Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : // C++ includes
      21             : #include <iomanip>
      22             : #include <fstream>
      23             : 
      24             : // Local includes
      25             : #include "libmesh/libmesh_config.h"
      26             : #include "libmesh/libmesh_logging.h"
      27             : #include "libmesh/mesh_base.h"
      28             : #include "libmesh/medit_io.h"
      29             : #include "libmesh/elem.h"
      30             : 
      31             : namespace libMesh
      32             : {
      33             : 
      34             : 
      35             : // ------------------------------------------------------------
      36             : // MEDITIO  members
      37           0 : void MEDITIO::write (const std::string & fname)
      38             : {
      39           0 :   if (this->mesh().processor_id() == 0)
      40           0 :     if (!this->binary())
      41           0 :       this->write_ascii  (fname);
      42           0 : }
      43             : 
      44             : 
      45             : 
      46           0 : void MEDITIO::write_nodal_data (const std::string & fname,
      47             :                                 const std::vector<Number> & soln,
      48             :                                 const std::vector<std::string> & names)
      49             : {
      50           0 :   LOG_SCOPE("write_nodal_data()", "MEDITIO");
      51             : 
      52           0 :   if (this->mesh().processor_id() == 0)
      53           0 :     if (!this->binary())
      54           0 :       this->write_ascii  (fname, &soln, &names);
      55           0 : }
      56             : 
      57             : 
      58             : 
      59           0 : void MEDITIO::write_ascii (const std::string & fname,
      60             :                            const std::vector<Number> * vec,
      61             :                            const std::vector<std::string> * solution_names)
      62             : {
      63             :   // Current lacks in implementation:
      64             :   //  (i)   only 3D meshes.
      65             :   //  (ii)  only QUAD4, TRI3, TET4 elements, others are omitted !
      66             :   //  (iii) no distinction between materials.
      67             :   //  (iv)  no vector output, just first scalar as output
      68             : 
      69             :   // libmesh_assert three dimensions (should be extended later)
      70           0 :   libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().mesh_dimension(), 3);
      71             : 
      72             :   // Open the output file stream
      73           0 :   std::ofstream out_stream (fname.c_str());
      74             : 
      75             :   // Make sure it opened correctly
      76           0 :   if (!out_stream.good())
      77           0 :     libmesh_file_error(fname.c_str());
      78             : 
      79             :   // Get a reference to the mesh
      80           0 :   const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
      81             : 
      82             :   // Begin interfacing with the MEdit data file
      83             :   {
      84             :     // header:
      85           0 :     out_stream << "MeshVersionFormatted  1\n";
      86           0 :     out_stream << "Dimension  3\n";
      87           0 :     out_stream << "# Mesh generated by libmesh\n\n";
      88             : 
      89             :     // write the nodes:
      90           0 :     out_stream << "# Set of mesh vertices\n";
      91           0 :     out_stream << "Vertices\n";
      92           0 :     out_stream << the_mesh.n_nodes() << "\n";
      93             : 
      94           0 :     for (auto v : make_range(the_mesh.n_nodes()))
      95           0 :       out_stream << the_mesh.point(v)(0) << " " << the_mesh.point(v)(1) << " " << the_mesh.point(v)(2) << " 0\n";
      96             :   }
      97             : 
      98             :   {
      99             :     // write the connectivity:
     100           0 :     out_stream << "\n# Set of Polys\n\n";
     101             : 
     102             :     // count occurrences of output elements:
     103           0 :     int n_tri3  = 0;
     104           0 :     int n_quad4 = 0;
     105           0 :     int n_tet4  = 0;
     106             : 
     107           0 :     for (const auto & elem : the_mesh.active_element_ptr_range())
     108             :       {
     109           0 :         if (elem->type() == TRI3)  n_tri3++;
     110           0 :         if (elem->type() == QUAD4) n_quad4++;
     111           0 :         if (elem->type() == QUAD9) n_quad4+=4; // (QUAD9 is written as 4 QUAD4.)
     112           0 :         if (elem->type() == TET4)  n_tet4++;
     113           0 :       }
     114             : 
     115             :     // First: write out TRI3 elements:
     116           0 :     out_stream << "Triangles\n";
     117           0 :     out_stream << n_tri3 << "\n";
     118             : 
     119           0 :     for (const auto & elem : the_mesh.active_element_ptr_range())
     120           0 :       if (elem->type() == TRI3)
     121           0 :         out_stream << elem->node_id(0)+1  << " "
     122           0 :                    << elem->node_id(1)+1  << " "
     123           0 :                    << elem->node_id(2)+1  << " 0\n";
     124             : 
     125             :     // Second: write out QUAD4 elements:
     126           0 :     out_stream << "Quadrilaterals\n";
     127           0 :     out_stream << n_quad4 << "\n";
     128             : 
     129           0 :     for (const auto & elem : the_mesh.active_element_ptr_range())
     130             :       {
     131           0 :         if (elem->type() == QUAD4)
     132             :           {
     133           0 :             out_stream << elem->node_id(0)+1  << " "
     134           0 :                        << elem->node_id(1)+1  << " "
     135           0 :                        << elem->node_id(2)+1  << " "
     136           0 :                        << elem->node_id(3)+1  <<" 0\n";
     137             :           } // if
     138           0 :         else if (elem->type() == QUAD9)
     139             :           {
     140           0 :             out_stream << elem->node_id(0)+1  << " "
     141           0 :                        << elem->node_id(4)+1  << " "
     142           0 :                        << elem->node_id(8)+1  << " "
     143           0 :                        << elem->node_id(7)+1  <<" 0\n";
     144           0 :             out_stream << elem->node_id(7)+1  << " "
     145           0 :                        << elem->node_id(8)+1  << " "
     146           0 :                        << elem->node_id(6)+1  << " "
     147           0 :                        << elem->node_id(3)+1  <<" 0\n";
     148           0 :             out_stream << elem->node_id(4)+1  << " "
     149           0 :                        << elem->node_id(1)+1  << " "
     150           0 :                        << elem->node_id(5)+1  << " "
     151           0 :                        << elem->node_id(8)+1  <<" 0\n";
     152           0 :             out_stream << elem->node_id(8)+1  << " "
     153           0 :                        << elem->node_id(5)+1  << " "
     154           0 :                        << elem->node_id(2)+1  << " "
     155           0 :                        << elem->node_id(6)+1  <<" 0\n";
     156             :           }
     157           0 :       }
     158             : 
     159             :     // Third: write out TET4 elements:
     160           0 :     out_stream << "Tetrahedra\n";
     161           0 :     out_stream << n_tet4 << "\n";
     162             : 
     163           0 :     for (const auto & elem : the_mesh.active_element_ptr_range())
     164           0 :       if (elem->type() == TET4)
     165             :         {
     166           0 :           out_stream << elem->node_id(0)+1  << " "
     167           0 :                      << elem->node_id(1)+1  << " "
     168           0 :                      << elem->node_id(2)+1  << " "
     169           0 :                      << elem->node_id(3)+1  <<" 0\n";
     170           0 :         }
     171             :   }
     172             :   // end of the out file
     173           0 :   out_stream << '\n' << "# end of file\n";
     174             : 
     175             :   // optionally write the data
     176           0 :   if ((solution_names != nullptr) &&
     177           0 :       (vec != nullptr))
     178             :     {
     179             :       // Open the ".bb" file stream
     180           0 :       std::size_t idx = fname.find_last_of(".");
     181           0 :       std::string bbname = fname.substr(0,idx) + ".bb";
     182             : 
     183           0 :       std::ofstream bbout (bbname.c_str());
     184             : 
     185             :       // Make sure it opened correctly
     186           0 :       if (!bbout.good())
     187           0 :         libmesh_file_error(bbname.c_str());
     188             : 
     189             :       // Header: 3: 3D mesh, 1: scalar output, 2: node-indexed
     190           0 :       const std::size_t n_vars = solution_names->size();
     191           0 :       bbout << "3 1 " << the_mesh.n_nodes() << " 2\n";
     192           0 :       for (auto n : make_range(the_mesh.n_nodes()))
     193           0 :         bbout << std::setprecision(this->ascii_precision()) << (*vec)[n*n_vars + scalar_idx] << " ";
     194           0 :       bbout << "\n";
     195           0 :     } // endif
     196           0 : }
     197             : 
     198             : } // namespace libMesh