gmsh_io.c

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

// $Id: gmsh_io.C 2789 2008-04-13 02:24:40Z roystgnr $

// 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

// This file was massively overhauled and extended by Martin L黷hi, [e-mail]

// C++ includes
#include <fstream>
#include <set>
#include <cstring> // std::memcpy, std::strncmp

// Local includes
#include "libmesh_config.h"
#include "gmsh_io.h"
#include "elem.h"
#include "mesh_base.h"
#include "boundary_info.h"


// anonymous namespace to hold local data
namespace 
{
  
  /**
   * Defines a structure to hold boundary element information.
   *
   * We use a set because it keeps the nodes unique and ordered, and can be
   * easily compared to another set of nodes (the ones on the element side)
   */
  struct boundaryElementInfo {
    std::set<unsigned int> nodes;
    unsigned int id;
  };

  /**
   * Defines mapping from libMesh element types to Gmsh element types.
   */
  struct elementDefinition {
    std::string label;
    std::vector<unsigned int> nodes;
    ElemType type;
    unsigned int exptype;
    unsigned int dim;
    unsigned int nnodes;
  };


  // maps from a libMesh element type to the proper
  // Gmsh elementDefinition.  Placing the data structure
  // here in this anonymous namespace gives us the
  // benefits of a global variable without the nasty
  // side-effects
  std::map<ElemType, elementDefinition> eletypes_exp;
  std::map<unsigned int, elementDefinition> eletypes_imp;



  // ------------------------------------------------------------
  // helper function to initialize the eletypes map
  void init_eletypes ()
  {
    if (eletypes_exp.empty() && eletypes_imp.empty())
      {
	// This should happen only once.  The first time this method
	// is called the eletypes data struture will be empty, and
	// we will fill it.  Any subsequent calls will find an initialized
	// eletypes map and will do nothing.

	//==============================
	// setup the element definitions
	elementDefinition eledef;

	// use "swap trick" from Scott Meyer's "Effective STL" to initialize
	// eledef.nodes vector
	
	// POINT (only Gmsh)
	{
	  eledef.exptype = 15;
          eledef.dim     = 0;
          eledef.nnodes  = 1;
	  eledef.nodes.clear();
	
          // import only
          eletypes_imp[15] = eledef;
	}

	// EDGE2
	{
	  eledef.type    = EDGE2;
          eledef.dim     = 1;
          eledef.nnodes  = 2;
	  eledef.exptype = 1;
	  eledef.nodes.clear();
	
	  eletypes_exp[EDGE2] = eledef;
          eletypes_imp[1]     = eledef;
	}
  
	// EDGE3
	{
	  eledef.type    = EDGE3;
          eledef.dim     = 1;
          eledef.nnodes  = 3;
	  eledef.exptype = 8;
	  eledef.nodes.clear();
	  
	  eletypes_exp[EDGE3] = eledef;
          eletypes_imp[8]     = eledef;
	}
      
	// TRI3
	{
          eledef.type    = TRI3;
          eledef.dim     = 2;
          eledef.nnodes  = 3;
	  eledef.exptype = 2;
	  eledef.nodes.clear();
	  
	  eletypes_exp[TRI3] = eledef;
          eletypes_imp[2] = eledef;
	}
      
	// TRI6
	{
          eledef.type    = TRI6;
          eledef.dim     = 2;
          eledef.nnodes  = 6;
	  eledef.exptype = 9;
	  eledef.nodes.clear();

	  eletypes_exp[TRI6] = eledef;
          eletypes_imp[9]    = eledef;
	}
      
	// QUAD4
	{
          eledef.type    = QUAD4;
          eledef.dim     = 2;
          eledef.nnodes  = 4;
	  eledef.exptype = 3;
	  eledef.nodes.clear();

	  eletypes_exp[QUAD4] = eledef;
          eletypes_imp[3]     = eledef;
	}
      
	// QUAD8
        // TODO: what should be done with this on writing?
	{
          eledef.type    = QUAD8;
          eledef.dim     = 2;
          eledef.nnodes  = 8;
	  eledef.exptype = 100;
	  const unsigned int nodes[] = {1,2,3,4,5,6,7,8};
	  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);

	  eletypes_exp[QUAD8] = eledef;
          eletypes_imp[10]    = eledef;
	}
      
	// QUAD9
	{
          eledef.type    = QUAD9;
          eledef.dim     = 2;
          eledef.nnodes  = 9;
	  eledef.exptype = 10;
	  eledef.nodes.clear();

	  eletypes_exp[QUAD9] = eledef;
          eletypes_imp[10]    = eledef;
	}
      
	// HEX8
	{
          eledef.type    = HEX8;
          eledef.dim     = 3;
          eledef.nnodes  = 8;
	  eledef.exptype = 5;
	  eledef.nodes.clear();

	  eletypes_exp[HEX8] = eledef;
          eletypes_imp[5]    = eledef;
	}
      
	// HEX20
        // TODO: what should be done with this on writing?
	{
          eledef.type    = HEX20;
          eledef.dim     = 3;
          eledef.nnodes  = 20;
	  eledef.exptype = 101;
	  const unsigned int nodes[] = {1,2,3,4,5,6,7,8,9,10,11,16,17,18,19,12,13,14,15,16};
	  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);

	  eletypes_exp[HEX20] = eledef;
          eletypes_imp[12]    = eledef;
	}
      
	// HEX27
	{
          eledef.type    = HEX27;
          eledef.dim     = 3;
          eledef.nnodes  = 27;
	  eledef.exptype = 12;
          const unsigned int nodes[] = {0,1,2,3,4,5,6,7,8,11,12,9,13,10,14,
                                        15,16,19,17,18,20,21,24,22,23,25,26};
	  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);

	  eletypes_exp[HEX27] = eledef;
          eletypes_imp[12]    = eledef;
	}
      
	// TET4
	{
          eledef.type    = TET4;
          eledef.dim     = 3;
          eledef.nnodes  = 4;
	  eledef.exptype = 4;
	  eledef.nodes.clear();

          eletypes_exp[TET4] = eledef;
          eletypes_imp[4]    = eledef;
	}
      
	// TET10
	{
          eledef.type    = TET10;
          eledef.dim     = 3;
          eledef.nnodes  = 10;
	  eledef.exptype = 11;
          const unsigned int nodes[] = {0,1,2,3,4,5,6,7,9,8};
	  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);
	  eletypes_exp[TET10] = eledef;
          eletypes_imp[11]    = eledef;
	}
      
	// PRISM6
	{
          eledef.type    = PRISM6;
          eledef.dim     = 3;
          eledef.nnodes  = 6;
	  eledef.exptype = 6;
	  eledef.nodes.clear();
	  
	  eletypes_exp[PRISM6] = eledef;
          eletypes_imp[6]      = eledef;
	}
      
	// PRISM15
        // TODO: what should be done with this on writing?
	{
          eledef.type    = PRISM15;
          eledef.dim     = 3;
          eledef.nnodes  = 15;
	  eledef.exptype = 103;
	  eledef.nodes.clear();
	  
	  eletypes_exp[PRISM15] = eledef;
          eletypes_imp[13] = eledef;
	}

	// PRISM18
	{
          eledef.type    = PRISM18;
          eledef.dim     = 3;
          eledef.nnodes  = 18;
	  eledef.exptype = 13;
          const unsigned int nodes[] = {0,1,2,3,4,5,6,8,9,7,10,11,
                                        12,14,13,15,17,16};
	  std::vector<unsigned int>(nodes, nodes+eledef.nnodes).swap(eledef.nodes);

	  eletypes_exp[PRISM18] = eledef;
          eletypes_imp[13]      = eledef;
	}

	// PYRAMID5
	{
          eledef.type    = PYRAMID5;
          eledef.dim     = 3;
          eledef.nnodes  = 5;
	  eledef.exptype = 7;
	  eledef.nodes.clear();
	  
	  eletypes_exp[PYRAMID5] = eledef;
          eletypes_imp[7]        = eledef;
	}

	//==============================      
      }
  }
  
} // end anonymous namespace


// ------------------------------------------------------------
// GmshIO  members
void GmshIO::read (const std::string& name)
{
  std::ifstream in (name.c_str());
  this->read_mesh (in);
}


void GmshIO::read_mesh(std::istream& in)
{
  // This is a serial-only process for now;
  // the Mesh should be read on processor 0 and
  // broadcast later
  libmesh_assert(libMesh::processor_id() == 0);

  libmesh_assert(in.good());

  // initialize the map with element types
  init_eletypes();

  // clear any data in the mesh
  MeshBase& mesh = MeshInput<MeshBase>::mesh();
  mesh.clear();

  const unsigned int dim = mesh.mesh_dimension();

  // some variables
  const int  bufLen = 256;
  char       buf[bufLen+1];
  int        format=0, size=0;
  Real       version = 1.0;
  
  // map to hold the node numbers for translation
  // note the the nodes can be non-consecutive
  std::map<unsigned int, unsigned int> nodetrans;
  
  {
    while (!in.eof()) {
      in >> buf;

      if (!std::strncmp(buf,"$MeshFormat",11))
        {
          in >> version >> format >> size;
          if(version != 2.0){
            std::cerr << "Error: Wrong msh file version " << version << "\n";
            libmesh_error();
          }
          if(format){
            std::cerr << "Error: Unknown data format for mesh\n";
            libmesh_error();
          }
        }
      
      // read the node block
      else if (!std::strncmp(buf,"$NOD",4) ||
	       !std::strncmp(buf,"$NOE",4) ||
	       !std::strncmp(buf,"$Nodes",6) 
	       )
        {
          unsigned int numNodes = 0;
          in >> numNodes;
          mesh.reserve_nodes (numNodes);

          // read in the nodal coordinates and form points.
          Real x, y, z;
          unsigned int id;
        
          // add the nodal coordinates to the mesh
          for (unsigned int i=0; i<numNodes; ++i)      
            {
              in >> id >> x >> y >> z;
              mesh.add_point (Point(x, y, z), i);
              nodetrans[id] = i;
            }
          // read the $ENDNOD delimiter
          in >> buf;
        }

      /**
       * Read the element block
       *
       * If the element dimension is smaller than the mesh dimension, this is a
       * boundary element and will be added to mesh.boundary_info.
       *
       * Because the elements might not yet exist, the sides are put on hold
       * until the elements are created, and inserted once reading elements is
       * finished
       */
      else if (!std::strncmp(buf,"$ELM",4) ||
               !std::strncmp(buf,"$Elements",9)
               )
        {
          unsigned int numElem = 0;
          std::vector< boundaryElementInfo > boundary_elem;

          // read how many elements are there, and reserve space in the mesh
          in >> numElem;
          mesh.reserve_elem (numElem);

          // read the elements
	  unsigned int elem_id_counter = 0;
          for (unsigned int iel=0; iel<numElem; ++iel)      
            {
              unsigned int id, type, physical, elementary,
                /* partition = 1,*/ nnodes, ntags;
	      // note - partition was assigned but never used - BSK
              if(version <= 1.0)
                {
                  in >> id >> type >> physical >> elementary >> nnodes;
                }
              else
                {
                  in >> id >> type >> ntags;
                  elementary = physical = /* partition = */ 1;
                  for(unsigned int j = 0; j < ntags; j++)
                    {
                      int tag;
                      in >> tag;
                      if(j == 0)
                        physical = tag;
                      else if(j == 1)
                        elementary = tag;
                      // else if(j == 2)
                      //  partition = tag;
                      // ignore any other tags for now
                    }
                }

              // consult the import element table which element to build
              const elementDefinition& eletype = eletypes_imp[type];
              nnodes = eletype.nnodes;

              // only elements that match the mesh dimension are added
              // if the element dimension is one less than dim, the nodes and
              // sides are added to the mesh.boundary_info
              if (eletype.dim == dim)
                {
                  // add the elements to the mesh
                  Elem* elem = Elem::build(eletype.type).release();
                  elem->set_id(elem_id_counter);
                  mesh.add_elem(elem);

		  // different to iel, lower dimensional elems aren't added
		  elem_id_counter++;

                  // check number of nodes. We cannot do that for version 2.0
                  if (version <= 1.0) 
                    {
                      if (elem->n_nodes() != nnodes)
                        {
                          std::cerr << "Number of nodes for element " << id
                                    << " of type " << eletypes_imp[type].type
                                    << " (Gmsh type " << type  
                                    << ") does not match Libmesh definition. "
                                    << "I expected " << elem->n_nodes()
                                    << " nodes, but got " << nnodes << "\n";
                          libmesh_error();
                        }
                    }