medit_io.c

一个用来实现偏微分方程中网格的计算库

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// 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_config.h"
#include "mesh_base.h"
#include "medit_io.h"
#include "elem.h"


// ------------------------------------------------------------
// MEDITIO  members
void MEDITIO::write (const std::string& fname)
{
  if (libMesh::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)
{
  if (libMesh::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 (MeshOutput<MeshBase>::mesh().mesh_dimension() == 3);

  // Open the output file stream
  std::ofstream out (fname.c_str());
  
  libmesh_assert (out.good());

  // Get a reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();
  
  // Begin interfacing with the MEdit data file
  {
    // header:
    out << "MeshVersionFormatted  1\n";
    out << "Dimension  3\n";
    out << "# Mesh generated by libmesh\n\n";

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

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

  {
    // write the connectivity:
    out << "\n# Set of Polys\n\n";
    
    // count occurrences of output elements:
    int n_tri3  = 0;
    int n_quad4 = 0;
    int n_tet4  = 0;

    {
//       const_active_elem_iterator       it (mesh.elements_begin());
//       const const_active_elem_iterator end(mesh.elements_end());

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
      
      for ( ; it != end; ++it)
	{
	  if ((*it)->type() == TRI3)  n_tri3++;
	  if ((*it)->type() == QUAD4) n_quad4++;
	  if ((*it)->type() == QUAD9) n_quad4+=4; // (QUAD9 is written as 4 QUAD4.)
	  if ((*it)->type() == TET4)  n_tet4++;
	} // for
    }
    
    // First: write out TRI3 elements:
    out << "Triangles\n";
    out << n_tri3 << "\n";

    {
//       const_active_elem_iterator       it (mesh.elements_begin());
//       const const_active_elem_iterator end(mesh.elements_end());

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
    
      for ( ; it != end; ++it) 
	if ((*it)->type() == TRI3)
	  out << (*it)->node(0)+1  << " " << (*it)->node(1)+1  << " " << (*it)->node(2)+1  << " 0\n";
    }
    
    // Second: write out QUAD4 elements:
    out << "Quadrilaterals\n";
    out << n_quad4 << "\n";

    {
//       const_active_elem_iterator       it (mesh.elements_begin());
//       const const_active_elem_iterator end(mesh.elements_end());

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 

      for ( ; it != end; ++it) 
	if ((*it)->type() == QUAD4)
	  {
	    out << (*it)->node(0)+1  << " " 
		<< (*it)->node(1)+1  << " " 
		<< (*it)->node(2)+1  << " " 
		<< (*it)->node(3)+1  <<" 0\n";
	  } // if
	else if ((*it)->type() == QUAD9)
	  {
	    out << (*it)->node(0)+1  << " " 
		<< (*it)->node(4)+1  << " " 
		<< (*it)->node(8)+1  << " " 
		<< (*it)->node(7)+1  <<" 0\n";
	    out << (*it)->node(7)+1  << " " 
		<< (*it)->node(8)+1  << " " 
		<< (*it)->node(6)+1  << " " 
		<< (*it)->node(3)+1  <<" 0\n";
	    out << (*it)->node(4)+1  << " " 
		<< (*it)->node(1)+1  << " " 
		<< (*it)->node(5)+1  << " " 
		<< (*it)->node(8)+1  <<" 0\n";
	    out << (*it)->node(8)+1  << " " 
		<< (*it)->node(5)+1  << " " 
		<< (*it)->node(2)+1  << " " 
		<< (*it)->node(6)+1  <<" 0\n";
	  } // if
    }

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

    {
//       const_active_elem_iterator       it (mesh.elements_begin());
//       const const_active_elem_iterator end(mesh.elements_end());

      MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
      
      for ( ; it != end; ++it) 
	if ((*it)->type() == TET4)
	  {
	    out << (*it)->node(0)+1  << " " 
		<< (*it)->node(1)+1  << " " 
		<< (*it)->node(2)+1  << " " 
		<< (*it)->node(3)+1  <<" 0\n";
	  } // if
    }
    
  }  
  // end of the out file
  out << '\n' << "# end of file\n";

  // optionally write the data
  if ((solution_names != NULL) &&
      (vec != NULL))
    {
      // Open the ".bb" file stream
      int idx = fname.find_last_of(".");
      std::string bbname = fname.substr(0,idx) + ".bb";
      
      std::ofstream bbout (bbname.c_str());
      libmesh_assert (bbout.good());      

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