medit_io.C

Go to the documentation of this file.

// The libMesh Finite Element Library.
// Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.

// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.

// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// C++ includes
#include <iomanip>
#include <fstream>

// Local includes
#include "libmesh/libmesh_config.h"
#include "libmesh/libmesh_logging.h"
#include "libmesh/mesh_base.h"
#include "libmesh/medit_io.h"
#include "libmesh/elem.h"

namespace libMesh
{


// ------------------------------------------------------------
// MEDITIO  members
void MEDITIO::write (const std::string & fname)
{
  if (this->mesh().processor_id() == 0)
    if (!this->binary())
      this->write_ascii  (fname);
}



void MEDITIO::write_nodal_data (const std::string & fname,
                                const std::vector<Number> & soln,
                                const std::vector<std::string> & names)
{
  LOG_SCOPE("write_nodal_data()", "MEDITIO");

  if (this->mesh().processor_id() == 0)
    if (!this->binary())
      this->write_ascii  (fname, &soln, &names);
}



void MEDITIO::write_ascii (const std::string & fname,
                           const std::vector<Number> * vec,
                           const std::vector<std::string> * solution_names)
{
  // Current lacks in implementation:
  //  (i)   only 3D meshes.
  //  (ii)  only QUAD4, TRI3, TET4 elements, others are omitted !
  //  (iii) no distinction between materials.
  //  (iv)  no vector output, just first scalar as output

  // libmesh_assert three dimensions (should be extended later)
  libmesh_assert_equal_to (MeshOutput<MeshBase>::mesh().mesh_dimension(), 3);

  // Open the output file stream
  std::ofstream out_stream (fname.c_str());

  // Make sure it opened correctly
  if (!out_stream.good())
    libmesh_file_error(fname.c_str());

  // Get a reference to the mesh
  const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();

  // Begin interfacing with the MEdit data file
  {
    // header:
    out_stream << "MeshVersionFormatted  1\n";
    out_stream << "Dimension  3\n";
    out_stream << "# Mesh generated by libmesh\n\n";

    // write the nodes:
    out_stream << "# Set of mesh vertices\n";
    out_stream << "Vertices\n";
    out_stream << the_mesh.n_nodes() << "\n";

    for (unsigned int v=0; v<the_mesh.n_nodes(); v++)
      out_stream << the_mesh.point(v)(0) << " " << the_mesh.point(v)(1) << " " << the_mesh.point(v)(2) << " 0\n";
  }

  {
    // write the connectivity:
    out_stream << "\n# Set of Polys\n\n";

    // count occurrences of output elements:
    int n_tri3  = 0;
    int n_quad4 = 0;
    int n_tet4  = 0;

    for (const auto & elem : the_mesh.active_element_ptr_range())
      {
        if (elem->type() == TRI3)  n_tri3++;
        if (elem->type() == QUAD4) n_quad4++;
        if (elem->type() == QUAD9) n_quad4+=4; // (QUAD9 is written as 4 QUAD4.)
        if (elem->type() == TET4)  n_tet4++;
      }

    // First: write out TRI3 elements:
    out_stream << "Triangles\n";
    out_stream << n_tri3 << "\n";

    for (const auto & elem : the_mesh.active_element_ptr_range())
      if (elem->type() == TRI3)
        out_stream << elem->node_id(0)+1  << " "
                   << elem->node_id(1)+1  << " "
                   << elem->node_id(2)+1  << " 0\n";

    // Second: write out QUAD4 elements:
    out_stream << "Quadrilaterals\n";
    out_stream << n_quad4 << "\n";

    for (const auto & elem : the_mesh.active_element_ptr_range())
      {
        if (elem->type() == QUAD4)
          {
            out_stream << elem->node_id(0)+1  << " "
                       << elem->node_id(1)+1  << " "
                       << elem->node_id(2)+1  << " "
                       << elem->node_id(3)+1  <<" 0\n";
          } // if
        else if (elem->type() == QUAD9)
          {
            out_stream << elem->node_id(0)+1  << " "
                       << elem->node_id(4)+1  << " "
                       << elem->node_id(8)+1  << " "
                       << elem->node_id(7)+1  <<" 0\n";
            out_stream << elem->node_id(7)+1  << " "
                       << elem->node_id(8)+1  << " "
                       << elem->node_id(6)+1  << " "
                       << elem->node_id(3)+1  <<" 0\n";
            out_stream << elem->node_id(4)+1  << " "
                       << elem->node_id(1)+1  << " "
                       << elem->node_id(5)+1  << " "
                       << elem->node_id(8)+1  <<" 0\n";
            out_stream << elem->node_id(8)+1  << " "
                       << elem->node_id(5)+1  << " "
                       << elem->node_id(2)+1  << " "
                       << elem->node_id(6)+1  <<" 0\n";
          }
      }

    // Third: write out TET4 elements:
    out_stream << "Tetrahedra\n";
    out_stream << n_tet4 << "\n";

    for (const auto & elem : the_mesh.active_element_ptr_range())
      if (elem->type() == TET4)
        {
          out_stream << elem->node_id(0)+1  << " "
                     << elem->node_id(1)+1  << " "
                     << elem->node_id(2)+1  << " "
                     << elem->node_id(3)+1  <<" 0\n";
        }
  }
  // end of the out file
  out_stream << '\n' << "# end of file\n";

  // optionally write the data
  if ((solution_names != nullptr) &&
      (vec != nullptr))
    {
      // Open the ".bb" file stream
      std::size_t idx = fname.find_last_of(".");
      std::string bbname = fname.substr(0,idx) + ".bb";

      std::ofstream bbout (bbname.c_str());

      // Make sure it opened correctly
      if (!bbout.good())
        libmesh_file_error(bbname.c_str());

      // Header: 3: 3D mesh, 1: scalar output, 2: node-indexed
      const std::size_t n_vars = solution_names->size();
      bbout << "3 1 " << the_mesh.n_nodes() << " 2\n";
      for (dof_id_type n=0; n<the_mesh.n_nodes(); n++)
        bbout << std::setprecision(10) << (*vec)[n*n_vars + scalar_idx] << " ";
      bbout << "\n";
    } // endif
}

} // namespace libMesh