gmsh_io.c

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

                  // add node pointers to the elements
                  int nod = 0;
                  // if there is a node translation table, use it
                  if (eletype.nodes.size() > 0)
                    for (unsigned int i=0; i<nnodes; i++)
                      {
                        in >> nod;
                        elem->set_node(eletype.nodes[i]) = mesh.node_ptr(nodetrans[nod]);
                      }
                  else
                    {
                      for (unsigned int i=0; i<nnodes; i++)
                        {
                          in >> nod;
                          elem->set_node(i) = mesh.node_ptr(nodetrans[nod]);
                        }
                    }
                } // if element.dim == dim
              // if this is a boundary
              else if (eletype.dim == dim-1)
                {
                  /**
                   * add the boundary element nodes to the set of nodes
                   */

                  boundaryElementInfo binfo;
                  std::set<unsigned int>::iterator iter = binfo.nodes.begin();
                  int nod = 0;
                  for (unsigned int i=0; i<nnodes; i++)
                    {
                      in >> nod;
                      mesh.boundary_info->add_node(nodetrans[nod], physical);
                      binfo.nodes.insert(iter, nodetrans[nod]);
                    }
                  binfo.id = physical;
                  boundary_elem.push_back(binfo);
                }
              /**
               * If the element yet another dimension, just read in the nodes
               * and throw them away
               */
              else
                {
                  int nod = 0;
                  for (unsigned int i=0; i<nnodes; i++)
                    in >> nod;
                }
            }//element loop
          // read the $ENDELM delimiter
          in >> buf;

          /**
           * If any lower dimensional elements have been found in the file,
           * try to add them to the mesh.boundary_info as sides and nodes with
           * the respecitve id's (called "physical" in Gmsh).
           */
          if (boundary_elem.size() > 0)
            {
              // create a index of the boundary nodes to easily locate which
              // element might have that boundary
              std::map<unsigned int, std::vector<unsigned int> > node_index;
              for (unsigned int i=0; i<boundary_elem.size(); i++)
                {
                  boundaryElementInfo binfo = boundary_elem[i];
                  std::set<unsigned int>::iterator iter = binfo.nodes.begin();
                  for (;iter!= binfo.nodes.end(); iter++)
                    node_index[*iter].push_back(i);
                }
              
              MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
              const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
              
              // iterate over all elements and see which boundary element has
              // the same set of nodes as on of the boundary elements previously read
              for ( ; it != end; ++it)
                {
                  const Elem* elem = *it;
                  for (unsigned int s=0; s<elem->n_sides(); s++)
                    if (elem->neighbor(s) == NULL)
                      {
                        AutoPtr<Elem> side (elem->build_side(s));
                        std::set<unsigned int> side_nodes;
                        std::set<unsigned int>::iterator iter = side_nodes.begin();

                        // make a set with all nodes from this side
                        // this allows for easy comparison                        
                        for (unsigned int ns=0; ns<side->n_nodes(); ns++)
                          side_nodes.insert(iter, side->node(ns)); 

                        // See whether one of the side node occurs in the list
                        // of tagged nodes. If we would loop over all side
                        // nodes, we would just get multiple hits, so taking
                        // node 0 is enough to do the job
                        unsigned int sn = side->node(0);
                        if (node_index.count(sn) > 0)
                          {
                            // Loop over all tagged ("physical") "sides" which
                            // contain the node sn (typically just 1 to
                            // three). For each of these the set of nodes is
                            // compared to the current element's side nodes
                            for (unsigned int n=0; n<node_index[sn].size(); n++)
                              {
                                unsigned int bidx = node_index[sn][n];
                                if (boundary_elem[bidx].nodes == side_nodes)
                                  mesh.boundary_info->add_side(elem, s, boundary_elem[bidx].id);
                              }
                          }
                      } // if elem->neighbor(s) == NULL
                } // element loop
            } // if boundary_elem.size() > 0
        } // if $ELM

    } // while !in.eof()

  }

}



void GmshIO::write (const std::string& name)
{
  if (libMesh::processor_id() == 0)
    {
      std::ofstream out (name.c_str());
      this->write_mesh (out);
    }
}


void GmshIO::write_nodal_data (const std::string& fname,
                               const std::vector<Number>& soln,
                               const std::vector<std::string>& names)
{
  //this->_binary = true;
  if (libMesh::processor_id() == 0)
    this->write_post  (fname, &soln, &names);
}


void GmshIO::write_mesh (std::ostream& out)
{
  // Be sure that the stream is valid.
  libmesh_assert (out.good());
  
  // initialize the map with element types
  init_eletypes();

  // Get a const reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();
  
  // Note: we are using version 2.0 of the gmsh output format.
  
  {
    // Write the file header.
    out << "$MeshFormat\n";
    out << "2.0 0 " << sizeof(Real) << '\n';
    out << "$EndMeshFormat\n";
  }

  {
    // write the nodes in (n x y z) format
    out << "$Nodes\n";
    out << mesh.n_nodes() << '\n';
    
    for (unsigned int v=0; v<mesh.n_nodes(); v++)
      out << mesh.node(v).id()+1 << " "
	  << mesh.node(v)(0) << " "
	  << mesh.node(v)(1) << " "
	  << mesh.node(v)(2) << '\n';
    out << "$EndNodes\n";
  }
  
  {
    // write the connectivity
    out << "$Elements\n";
    out << mesh.n_active_elem() << '\n';

    MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
    const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
    
    // loop over the elements
    for ( ; it != end; ++it)
      {
        const Elem* elem = *it;

	// Make sure we have a valid entry for
	// the current element type.
	libmesh_assert (eletypes_exp.count(elem->type()));

        // consult the export element table 
        const elementDefinition& eletype = eletypes_exp[elem->type()];

	// The element mapper better not require any more nodes
	// than are present in the current element!
        libmesh_assert (eletype.nodes.size() <= elem->n_nodes());
	
        // elements ids are 1 based in Gmsh
        out << elem->id()+1 << " ";

        // element type
        out << eletype.exptype;

        // write the number of tags and
        // tag1 (physical entity), and tag2 (geometric entity)
        out << " 3 1 1 ";

        // write the partition the element belongs to
        out << elem->processor_id()+1 << " ";

        // if there is a node translation table, use it
        if (eletype.nodes.size() > 0)
          for (unsigned int i=0; i < elem->n_nodes(); i++)
            out << elem->node(eletype.nodes[i])+1 << " "; // gmsh is 1-based 
        // otherwise keep the same node order
        else
          for (unsigned int i=0; i < elem->n_nodes(); i++)
            out << elem->node(i)+1 << " ";                  // gmsh is 1-based 
        out << "\n";
      } // element loop
    out << "$EndElements\n";
  }
}


void GmshIO::write_post (const std::string& fname,
                         const std::vector<Number>* v,
                         const std::vector<std::string>* solution_names)
{

  // Should only do this on processor 0!
  libmesh_assert (libMesh::processor_id() == 0);
  
  // Create an output stream
  std::ofstream out(fname.c_str());

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

  if (!out.good())
    {
      std::cerr << "ERROR: opening output file " << fname
		<< std::endl;
      libmesh_error();
    }

  // create a character buffer
  char buf[80];

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

  //  write the data
  if ((solution_names != NULL) && (v != NULL))
    {      
      const unsigned int n_vars = solution_names->size();
    
      if (!(v->size() == mesh.n_nodes()*n_vars))
        std::cerr << "ERROR: v->size()=" << v->size()
                  << ", mesh.n_nodes()=" << mesh.n_nodes()
                  << ", n_vars=" << n_vars
                  << ", mesh.n_nodes()*n_vars=" << mesh.n_nodes()*n_vars
                  << "\n";
      
      libmesh_assert (v->size() == mesh.n_nodes()*n_vars);

      // write the header
      out << "$PostFormat\n";
      if (this->binary())
        out << "1.2 1 " << sizeof(double) << "\n";
      else
        out << "1.2 0 " << sizeof(double) << "\n";
      out << "$EndPostFormat\n";

      // Loop over the elements to see how much of each type there are
      unsigned int n_points=0, n_lines=0, n_triangles=0, n_quadrangles=0,
        n_tetrahedra=0, n_hexahedra=0, n_prisms=0, n_pyramids=0;
      unsigned int n_scalar=0, n_vector=0, n_tensor=0;
      unsigned int nb_text2d=0, nb_text2d_chars=0, nb_text3d=0, nb_text3d_chars=0;

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


        for ( ; it != end; ++it)
          {
            const ElemType elemtype = (*it)->type();
            
            switch (elemtype)
              {
              case EDGE2:
              case EDGE3:
              case EDGE4:
                {
                  n_lines += 1;
                  break;
                }
              case TRI3:
              case TRI6:
                {
                  n_triangles += 1;
                  break;
                }
              case QUAD4:
              case QUAD8:
              case QUAD9:
                {
                  n_quadrangles += 1;
                  break;
                }
              case TET4:
              case TET10:
                {
                  n_tetrahedra += 1;
                  break;
                }
              case HEX8:
              case HEX20:
              case HEX27:
                {
                  n_hexahedra += 1;
                  break;
                }
              case PRISM6:
              case PRISM15:
              case PRISM18:
                {
                  n_prisms += 1;
                  break;
                }
              case PYRAMID5:
                {
                  n_pyramids += 1;
                  break;
                }
              default:
                {
                  std::cerr << "ERROR: Not existant element type "
                            << (*it)->type() << std::endl;
                  libmesh_error();
                }
              }
          }
      }

      // create a view for each variable
      for (unsigned int ivar=0; ivar < n_vars; ivar++)
        {
          std::string varname = (*solution_names)[ivar];

          // at the moment, we just write out scalar quantities
          // later this should be made configurable through
          // options to the writer class
          n_scalar = 1;
      
          // write the variable as a view, and the number of time steps
          out << "$View\n" << varname << " " << 1 << "\n";
      
          // write how many of each geometry type are written
          out << n_points * n_scalar << " "
              << n_points * n_vector << " "
              << n_points * n_tensor << " "
              << n_lines * n_scalar << " "
              << n_lines * n_vector << " "
              << n_lines * n_tensor << " "
              << n_triangles * n_scalar << " "
              << n_triangles * n_vector << " "
              << n_triangles * n_tensor << " "
              << n_quadrangles * n_scalar << " "
              << n_quadrangles * n_vector << " "
              << n_quadrangles * n_tensor << " "
              << n_tetrahedra * n_scalar << " "
              << n_tetrahedra * n_vector << " "
              << n_tetrahedra * n_tensor << " "
              << n_hexahedra * n_scalar << " "
              << n_hexahedra * n_vector << " "
              << n_hexahedra * n_tensor << " "
              << n_prisms * n_scalar << " "
              << n_prisms * n_vector << " "
              << n_prisms * n_tensor << " "
              << n_pyramids * n_scalar << " "
              << n_pyramids * n_vector << " "
              << n_pyramids * n_tensor << " "
              << nb_text2d << " "
              << nb_text2d_chars << " "
              << nb_text3d << " "
              << nb_text3d_chars << "\n";
      
          // if binary, write a marker to identify the endianness of the file
          if (this->binary())
            {
              const int one = 1;
              std::memcpy(buf, &one, sizeof(int));
              out.write(buf, sizeof(int));
            }

          // the time steps (there is just 1 at the moment)
          if (this->binary())
            {
              double one = 1;
              std::memcpy(buf, &one, sizeof(double));
              out.write(buf, sizeof(double));
            }
          else
            out << "1\n";

          // Loop over the elements and write out the data
          MeshBase::const_element_iterator       it  = mesh.active_elements_begin();
          const MeshBase::const_element_iterator end = mesh.active_elements_end(); 
    
          for ( ; it != end; ++it)
            {
              const Elem* elem = *it;
	
              // this is quite crappy, but I did not invent that file format!
              for (unsigned int d=0; d<3; d++)  // loop over the dimensions 
                {
                  for (unsigned int n=0; n < elem->n_vertices(); n++)   // loop over vertices
                    {
                      const Point vertex = elem->point(n);
                      if (this->binary())
                        {
                          double tmp = vertex(d);
                          std::memcpy(buf, &tmp, sizeof(double));
                          out.write(reinterpret_cast<char *>(buf), sizeof(double));
                        }
                      else
                        out << vertex(d) << " ";
                    }
                  if (!this->binary())
                    out << "\n";
                }

              // now finally write out the data
              for (unsigned int i=0; i < elem->n_vertices(); i++)   // loop over vertices
                if (this->binary())
                  {
#ifdef USE_COMPLEX_NUMBERS
		    std::cout << "WARNING: Gmsh::write_post does not fully support "
			      << "complex numbers. Will only write the real part of "
			      << "variable " << varname << std::endl;
#endif
                    double tmp = libmesh_real((*v)[elem->node(i)*n_vars + ivar]);
                    std::memcpy(buf, &tmp, sizeof(double));
                    out.write(reinterpret_cast<char *>(buf), sizeof(double));
                  }
                else
		  {
#ifdef USE_COMPLEX_NUMBERS
		    std::cout << "WARNING: Gmsh::write_post does not fully support "
			      << "complex numbers. Will only write the real part of "
			      << "variable " << varname << std::endl;
#endif
		    out << libmesh_real((*v)[elem->node(i)*n_vars + ivar]) << "\n";
		  }
            }
          if (this->binary())
            out << "\n";
          out << "$EndView\n";

        } // end variable loop (writing the views)
    }

}