libMesh: diva_io.C Source File

00001 // The libMesh Finite Element Library.
00002 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
00003 
00004 // This library is free software; you can redistribute it and/or
00005 // modify it under the terms of the GNU Lesser General Public
00006 // License as published by the Free Software Foundation; either
00007 // version 2.1 of the License, or (at your option) any later version.
00008 
00009 // This library is distributed in the hope that it will be useful,
00010 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00011 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00012 // Lesser General Public License for more details.
00013 
00014 // You should have received a copy of the GNU Lesser General Public
00015 // License along with this library; if not, write to the Free Software
00016 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00017 
00018 
00019 // C++ includes
00020 #include <fstream>
00021 
00022 // Local includes
00023 #include "libmesh/diva_io.h"
00024 #include "libmesh/boundary_mesh.h"
00025 #include "libmesh/mesh_tools.h"
00026 #include "libmesh/elem.h"
00027 #include "libmesh/boundary_info.h"
00028 
00029 namespace libMesh
00030 {
00031 
00032 // ------------------------------------------------------------
00033 // DivaIO class members
00034 void DivaIO::write (const std::string& fname)
00035 {
00036   // We may need to gather a ParallelMesh to output it, making that
00037   // const qualifier in our constructor a dirty lie
00038   MeshSerializer serialize(const_cast<MeshBase&>(this->mesh()), !_is_parallel_format);
00039 
00040   // Open the output file stream
00041   std::ofstream out_file(fname.c_str());
00042 
00043   // Make sure it opened correctly
00044   if (!out_file.good())
00045     libmesh_file_error(fname.c_str());
00046 
00047   this->write_stream (out_file);
00048 }
00049 
00050 
00051 
00052 
00053 void DivaIO::write_stream (std::ostream& out_file)
00054 {
00055   /*
00056     From Kelly: (kelly@tacc.utexas.edu)
00057 
00058     Ok, the following is the format:
00059 
00060     #points #triangles #quads #tets #prisms #pyramids #hexs
00061     loop over all points (written out x y z x y z ...)
00062     loop over all triangles (written out i1 i2 i3) (These are indices into
00063     the points going from
00064     1 to #points)
00065     loop over all quads (written out i1 i2 i3 i4) (Same numbering scheme)
00066     loop over all triangles and quads (write out b1) (This is a boundary
00067     condition for each
00068     triangle and each
00069     hex. You can put
00070     anything you want
00071     here)
00072     loop over all tets (written out i1 i2 i3 i4) (Same)
00073     loop over all pyramids (written out i1 i2 i3 i4 i5) (Same)
00074     loop over all prisms (written out i1 i2 i3 i4 i5 i6) (Same)
00075     loop over all hexs (written out i1 i2 i3 i4 i5 i6 i7 i8) (Same)
00076 
00077   */
00078 
00079   // Be sure the stream has been created successfully.
00080   libmesh_assert (out_file.good());
00081 
00082   // Can't use a constant mesh reference since we have to
00083   // sync the boundary info.
00084   libmesh_here();
00085   libMesh::err << "WARNING...  Sure you want to do this?"
00086                << std::endl;
00087   MeshBase& the_mesh = const_cast<MeshBase&>
00088     (MeshOutput<MeshBase>::mesh());
00089 
00090   if (the_mesh.mesh_dimension() < 3)
00091     {
00092       libMesh::err << "WARNING: DIVA only supports 3D meshes.\n\n"
00093                    << "Exiting without producing output.\n";
00094       return;
00095     }
00096 
00097 
00098 
00099   BoundaryMesh boundary_mesh
00100     (the_mesh.comm(),
00101      cast_int<unsigned char>(the_mesh.mesh_dimension()-1));
00102   the_mesh.get_boundary_info().sync(boundary_mesh);
00103 
00104 
00108   out_file << the_mesh.n_nodes()
00109            << ' '
00110            << (MeshTools::n_active_elem_of_type(boundary_mesh,TRI3) +
00111                MeshTools::n_active_elem_of_type(boundary_mesh,TRI6)*4)
00112            << ' '
00113            << (MeshTools::n_active_elem_of_type(boundary_mesh, QUAD4) +
00114                MeshTools::n_active_elem_of_type(boundary_mesh, QUAD8) +
00115                MeshTools::n_active_elem_of_type(boundary_mesh, QUAD9)*4)
00116            << ' '
00117            << (MeshTools::n_active_elem_of_type(the_mesh, TET4) +
00118                MeshTools::n_active_elem_of_type(the_mesh, TET10)*8)
00119            << ' '
00120            << MeshTools::n_active_elem_of_type(the_mesh, PYRAMID5)
00121            << ' '
00122            << (MeshTools::n_active_elem_of_type(the_mesh, PRISM6) +
00123                MeshTools::n_active_elem_of_type(the_mesh, PRISM18)*8)
00124            << ' '
00125            << (MeshTools::n_active_elem_of_type(the_mesh, HEX8)  +
00126                MeshTools::n_active_elem_of_type(the_mesh, HEX20) +
00127                MeshTools::n_active_elem_of_type(the_mesh, HEX27)*8)
00128            << ' '
00129            << '\n';
00130 
00131   boundary_mesh.clear();
00132 
00133 
00137   for (unsigned int v=0; v<the_mesh.n_nodes(); v++)
00138     the_mesh.point(v).write_unformatted(out_file);
00139 
00140 
00144   {
00148     for(unsigned int e=0; e<the_mesh.n_elem(); e++)
00149       if (the_mesh.elem(e)->active())
00150         for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
00151           if (the_mesh.elem(e)->neighbor(s) == NULL)
00152             {
00153               const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
00154 
00155               if (side->type() == TRI3)
00156                 {
00157                   out_file << side->node(0)+1 << " "
00158                            << side->node(1)+1 << " "
00159                            << side->node(2)+1 << '\n';
00160                 }
00161               else if (side->type() == TRI6)
00162                 {
00163                   out_file << side->node(0)+1 << " "
00164                            << side->node(3)+1 << " "
00165                            << side->node(5)+1 << '\n'
00166 
00167                            << side->node(3)+1 << " "
00168                            << side->node(1)+1 << " "
00169                            << side->node(4)+1 << '\n'
00170 
00171                            << side->node(5)+1 << " "
00172                            << side->node(4)+1 << " "
00173                            << side->node(2)+1 << '\n'
00174 
00175                            << side->node(3)+1 << " "
00176                            << side->node(4)+1 << " "
00177                            << side->node(5)+1 << '\n';
00178                 }
00179             }
00180 
00181 
00185     for(unsigned int e=0; e<the_mesh.n_elem(); e++)
00186       if (the_mesh.elem(e)->active())
00187         for (unsigned int s=0; s<the_mesh.elem(e)->n_sides(); s++)
00188           if (the_mesh.elem(e)->neighbor(s) == NULL)
00189             {
00190               const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
00191 
00192               if ((side->type() == QUAD4) ||
00193                   (side->type() == QUAD8)  )
00194                 {
00195                   out_file << side->node(0)+1 << " "
00196                            << side->node(1)+1 << " "
00197                            << side->node(2)+1 << " "
00198                            << side->node(3)+1 << '\n';
00199                 }
00200               else if (side->type() == QUAD9)
00201                 {
00202                   out_file << side->node(0)+1 << " "
00203                            << side->node(4)+1 << " "
00204                            << side->node(8)+1 << " "
00205                            << side->node(7)+1 << '\n'
00206 
00207                            << side->node(4)+1 << " "
00208                            << side->node(1)+1 << " "
00209                            << side->node(5)+1 << " "
00210                            << side->node(8)+1 << '\n'
00211 
00212                            << side->node(7)+1 << " "
00213                            << side->node(8)+1 << " "
00214                            << side->node(6)+1 << " "
00215                            << side->node(3)+1 << '\n'
00216 
00217                            << side->node(8)+1 << " "
00218                            << side->node(5)+1 << " "
00219                            << side->node(2)+1 << " "
00220                            << side->node(6)+1 << '\n';
00221                 }
00222             }
00223   }
00224 
00225 
00226 
00230   {
00234     for(unsigned int e=0; e<the_mesh.n_elem(); e++)
00235       if (the_mesh.elem(e)->active())
00236         for (unsigned short s=0; s<the_mesh.elem(e)->n_sides(); s++)
00237           if (the_mesh.elem(e)->neighbor(s) == NULL)
00238             {
00239               const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
00240 
00241               if ((side->type() == TRI3) ||
00242                   (side->type() == TRI6)  )
00243 
00244                 out_file << the_mesh.get_boundary_info().boundary_id(the_mesh.elem(e), s)
00245                          << '\n';
00246             }
00247 
00248 
00252     for(unsigned int e=0; e<the_mesh.n_elem(); e++)
00253       if (the_mesh.elem(e)->active())
00254         for (unsigned short s=0; s<the_mesh.elem(e)->n_sides(); s++)
00255           if (the_mesh.elem(e)->neighbor(s) == NULL)
00256             {
00257               const UniquePtr<Elem> side(the_mesh.elem(e)->build_side(s));
00258 
00259               if ((side->type() == QUAD4) ||
00260                   (side->type() == QUAD8) ||
00261                   (side->type() == QUAD9))
00262 
00263                 out_file << the_mesh.get_boundary_info().boundary_id(the_mesh.elem(e), s);
00264             }
00265   }
00266 
00267 
00268 
00272   for (unsigned int e=0; e<the_mesh.n_elem(); e++)
00273     if (the_mesh.elem(e)->active())
00274       {
00275         if (the_mesh.elem(e)->type() == TET4)
00276           {
00277             out_file << the_mesh.elem(e)->node(0)+1 << " "
00278                      << the_mesh.elem(e)->node(1)+1 << " "
00279                      << the_mesh.elem(e)->node(2)+1 << " "
00280                      << the_mesh.elem(e)->node(3)+1 << '\n';
00281           }
00282         else if (the_mesh.elem(e)->type() == TET10)
00283           {
00284             out_file << the_mesh.elem(e)->node(0)+1 << " "
00285                      << the_mesh.elem(e)->node(4)+1 << " "
00286                      << the_mesh.elem(e)->node(6)+1 << " "
00287                      << the_mesh.elem(e)->node(7)+1 << '\n';
00288 
00289             out_file << the_mesh.elem(e)->node(4)+1 << " "
00290                      << the_mesh.elem(e)->node(1)+1 << " "
00291                      << the_mesh.elem(e)->node(5)+1 << " "
00292                      << the_mesh.elem(e)->node(8)+1 << '\n';
00293 
00294             out_file << the_mesh.elem(e)->node(6)+1 << " "
00295                      << the_mesh.elem(e)->node(5)+1 << " "
00296                      << the_mesh.elem(e)->node(2)+1 << " "
00297                      << the_mesh.elem(e)->node(9)+1 << '\n';
00298 
00299             out_file << the_mesh.elem(e)->node(7)+1 << " "
00300                      << the_mesh.elem(e)->node(8)+1 << " "
00301                      << the_mesh.elem(e)->node(9)+1 << " "
00302                      << the_mesh.elem(e)->node(3)+1 << '\n';
00303 
00304             out_file << the_mesh.elem(e)->node(4)+1 << " "
00305                      << the_mesh.elem(e)->node(8)+1 << " "
00306                      << the_mesh.elem(e)->node(6)+1 << " "
00307                      << the_mesh.elem(e)->node(7)+1 << '\n';
00308 
00309             out_file << the_mesh.elem(e)->node(4)+1 << " "
00310                      << the_mesh.elem(e)->node(5)+1 << " "
00311                      << the_mesh.elem(e)->node(6)+1 << " "
00312                      << the_mesh.elem(e)->node(8)+1 << '\n';
00313 
00314             out_file << the_mesh.elem(e)->node(6)+1 << " "
00315                      << the_mesh.elem(e)->node(5)+1 << " "
00316                      << the_mesh.elem(e)->node(9)+1 << " "
00317                      << the_mesh.elem(e)->node(8)+1 << '\n';
00318 
00319             out_file << the_mesh.elem(e)->node(6)+1 << " "
00320                      << the_mesh.elem(e)->node(8)+1 << " "
00321                      << the_mesh.elem(e)->node(9)+1 << " "
00322                      << the_mesh.elem(e)->node(7)+1 << '\n';
00323           }
00324       }
00325 
00326 
00330   for (unsigned int e=0; e<the_mesh.n_elem(); e++)
00331     if (the_mesh.elem(e)->active())
00332       if (the_mesh.elem(e)->type() == PYRAMID5)
00333         {
00334           out_file << the_mesh.elem(e)->node(0)+1 << " "
00335                    << the_mesh.elem(e)->node(1)+1 << " "
00336                    << the_mesh.elem(e)->node(2)+1 << " "
00337                    << the_mesh.elem(e)->node(3)+1 << " "
00338                    << the_mesh.elem(e)->node(4)+1 << '\n';
00339         }
00340 
00341 
00342 
00346   for (unsigned int e=0; e<the_mesh.n_elem(); e++)
00347     if (the_mesh.elem(e)->active())
00348       {
00349         if (the_mesh.elem(e)->type() == PRISM6)
00350           {
00351             out_file << the_mesh.elem(e)->node(0)+1 << " "
00352                      << the_mesh.elem(e)->node(1)+1 << " "
00353                      << the_mesh.elem(e)->node(2)+1 << " "
00354                      << the_mesh.elem(e)->node(3)+1 << " "
00355                      << the_mesh.elem(e)->node(4)+1 << " "
00356                      << the_mesh.elem(e)->node(5)+1 << '\n';
00357           }
00358         else if (the_mesh.elem(e)->type() == PRISM18)
00359           libmesh_error_msg("PRISM18 element type not supported.");
00360       }
00361 
00362 
00366   for (unsigned int e=0; e<the_mesh.n_elem(); e++)
00367     if (the_mesh.elem(e)->active())
00368       if ((the_mesh.elem(e)->type() == HEX8)   ||
00369           (the_mesh.elem(e)->type() == HEX20) ||
00370           (the_mesh.elem(e)->type() == HEX27)   )
00371         {
00372           std::vector<dof_id_type> conn;
00373           for (unsigned int se=0; se<the_mesh.elem(e)->n_sub_elem(); se++)
00374             {
00375               the_mesh.elem(e)->connectivity(se, TECPLOT, conn);
00376 
00377               out_file << conn[0] << ' '
00378                        << conn[1] << ' '
00379                        << conn[2] << ' '
00380                        << conn[3] << ' '
00381                        << conn[4] << ' '
00382                        << conn[5] << ' '
00383                        << conn[6] << ' '
00384                        << conn[7] << '\n';
00385             }
00386         }
00387 }
00388 
00389 } // namespace libMesh