unstructured_mesh.c

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

	  if (name.rfind(".xda") < name.size())
	    {
	      xdr_io.binary() = false;
	      xdr_io.read (name);
	    }
	  else // .xdr* ==> true binary XDR file
	    {
	      xdr_io.binary() = true;
	      xdr_io.read (name);
	    }

	  // The xdr_io object gets constructed with legacy() == false.
	  // if legacy() == true then it means that a legacy file was detected and
	  // thus processor 0 performed the read. We therefore need to broadcast the
	  // mesh.  Further, for this flavor of mesh solution data ordering is tied
	  // to the node ordering, so we better not reorder the nodes!
	  if (xdr_io.legacy())
	    {
	      skip_renumber_nodes_and_elements = true;
	      MeshCommunication().broadcast(*this);
	    }
	}      
    }

  // Serial mesh formats
  else
    {
      START_LOG("read()", "Mesh");  
  
      // Read the file based on extension.  Only processor 0
      // needs to read the mesh.  It will then broadcast it and
      // the other processors will pick it up
      if (libMesh::processor_id() == 0)
	{
	  // Nasty hack for reading/writing zipped files
	  std::string new_name = name;
	  if (name.size() - name.rfind(".bz2") == 4)
	    {
	      new_name.erase(new_name.end() - 4, new_name.end());
	      std::string system_string = "bunzip2 -f -k ";
	      system_string += name;
	      START_LOG("system(bunzip2)", "Mesh");
	      std::system(system_string.c_str());
	      STOP_LOG("system(bunzip2)", "Mesh");
	    }

	  if (new_name.rfind(".mat") < new_name.size())
	    MatlabIO(*this).read(new_name);
	  
	  else if (new_name.rfind(".ucd") < new_name.size())
	    UCDIO(*this).read (new_name);
	  
	  else if (new_name.rfind(".exd") < new_name.size() ||
		   new_name.rfind(".e") < new_name.size())
	    ExodusII_IO(*this).read (new_name);
	  
	  else if ((new_name.rfind(".off")  < new_name.size()) ||
		   (new_name.rfind(".ogl")  < new_name.size()) ||
		   (new_name.rfind(".oogl") < new_name.size()))
	    OFFIO(*this).read (new_name);
     
	  else if (new_name.rfind(".mgf") < new_name.size())
	    LegacyXdrIO(*this,true).read_mgf (new_name);
      
	  else if (new_name.rfind(".unv") < new_name.size())
	    {
	      if (mesh_data == NULL)
		{
		  std::cerr << "Error! You must pass a "
			    << "valid MeshData pointer to "
			    << "read UNV files!" << std::endl;
		  libmesh_error();
		}
	      UNVIO(*this, *mesh_data).read (new_name);
	    }
      
	  else if ((new_name.rfind(".node")  < new_name.size()) ||
		   (new_name.rfind(".ele")   < new_name.size()))
	    TetGenIO(*this,mesh_data).read (new_name);

	  else if (new_name.rfind(".msh") < new_name.size())
	    GmshIO(*this).read (new_name);

	  else if (new_name.rfind(".gmv") < new_name.size())
	    GMVIO(*this).read (new_name);

	  else if (new_name.rfind(".vtu") < new_name.size())
	    VTKIO(*this).read(new_name);
      
	  else
	    {
	      std::cerr << " ERROR: Unrecognized file extension: " << name
			<< "\n   I understand the following:\n\n"
			<< "     *.mat  -- Matlab triangular ASCII file\n"
			<< "     *.ucd  -- AVS's ASCII UCD format\n"
			<< "     *.gmv  -- LANL's General Mesh Viewer format\n"
			<< "     *.off  -- OOGL OFF surface format\n"
			<< "     *.exd  -- Sandia's ExodusII format\n"
			<< "     *.e    -- Sandia's ExodusII format\n"
			<< "     *.xda  -- Internal ASCII format\n"
			<< "     *.xdr  -- Internal binary format,\n"
			<< "               compatible with XdrMGF\n"
			<< "     *.unv  -- I-deas Universal format\n"
			<< std::endl;
	      libmesh_error();	  
	    }    
	  
	  // If we temporarily decompressed a .bz2 file, remove the
	  // uncompressed version
	  if (name.size() - name.rfind(".bz2") == 4)
	    std::remove(new_name.c_str());
	}
      
  
      STOP_LOG("read()", "Mesh");

      // Send the mesh & bcs (which are now only on processor 0) to the other
      // processors
      MeshCommunication().broadcast (*this);
    }

  
  // Done reading the mesh.  Now prepare it for use.
  this->prepare_for_use(skip_renumber_nodes_and_elements);

}



void UnstructuredMesh::write (const std::string& name,
			      MeshData* mesh_data)
{
  // parallel formats are special -- they may choose to write
  // separate files, let's not try to handle the zipping here.
  if (is_parallel_file_format(name))
    {	
      // no need to handling bz2 files here -- the Xdr class does that.
      if (name.rfind(".xda") < name.size())
	XdrIO(*this).write(name);
	
      else if (name.rfind(".xdr") < name.size())
	XdrIO(*this,true).write(name);
    }

  // serial file formats
  else
    {
      START_LOG("write()", "Mesh");

      // Nasty hack for reading/writing zipped files
      std::string new_name = name;
      if (name.size() - name.rfind(".bz2") == 4)
	new_name.erase(new_name.end() - 4, new_name.end());
  
      // New scope so that io will close before we try to zip the file
      {
	// Write the file based on extension
	if (new_name.rfind(".dat") < new_name.size())
	  TecplotIO(*this).write (new_name);
	
	else if (new_name.rfind(".plt") < new_name.size())
	  TecplotIO(*this,true).write (new_name);
	
	else if (new_name.rfind(".ucd") < new_name.size())
	  UCDIO (*this).write (new_name);
	
	else if (new_name.rfind(".gmv") < new_name.size())
	  if (this->n_partitions() > 1)
	    GMVIO(*this).write (new_name);
	  else
	    {
	      GMVIO io(*this);
	      io.partitioning() = false;
	      io.write (new_name);
	    }
	
	else if (new_name.rfind(".ugrid") < new_name.size())
	  DivaIO(*this).write(new_name);
    
	else if (new_name.rfind(".mgf")  < new_name.size())
	  LegacyXdrIO(*this,true).write_mgf(new_name);
	
	else if (new_name.rfind(".unv") < new_name.size())
	  {
	    if (mesh_data == NULL)
	      {
		std::cerr << "Error! You must pass a "
			  << "valid MeshData pointer to "
			  << "write UNV files!" << std::endl;
		libmesh_error();
	      }
	    UNVIO(*this, *mesh_data).write (new_name);
	  }
	
	else if (new_name.rfind(".mesh") < new_name.size())
	  MEDITIO(*this).write (new_name);
	
	else if (new_name.rfind(".poly") < new_name.size())
	  TetGenIO(*this).write (new_name);
	
	else if (new_name.rfind(".msh") < new_name.size())
	  GmshIO(*this).write (new_name);
	
	else if (new_name.rfind(".fro") < new_name.size())
	  FroIO(*this).write (new_name);
	
	else if (new_name.rfind(".vtu") < new_name.size())
	  VTKIO(*this).write (new_name);
	
	else
	  {
	    std::cerr << " ERROR: Unrecognized file extension: " << name
		      << "\n   I understand the following:\n\n"
		      << "     *.dat   -- Tecplot ASCII file\n"
		      << "     *.plt   -- Tecplot binary file\n"
		      << "     *.ucd   -- AVS's ASCII UCD format\n"
		      << "     *.ugrid -- Kelly's DIVA ASCII format\n"
		      << "     *.gmv   -- LANL's GMV (General Mesh Viewer) format\n"
		      << "     *.xda   -- Internal ASCII format\n"
		      << "     *.xdr   -- Internal binary format,\n"
		      << "                compatible with XdrMGF\n"
		      << "     *.unv   -- I-deas Universal format\n"
		      << "     *.mesh  -- MEdit mesh format\n"
		      << "     *.poly  -- TetGen ASCII file\n"
		      << "     *.msh   -- GMSH ASCII file\n"
		      << "     *.fro   -- ACDL's surface triangulation file\n"
		      << std::endl
		      << "\n Exiting without writing output\n";
	  }    
      }
  
      // Nasty hack for reading/writing zipped files
      if (name.size() - name.rfind(".bz2") == 4)
	{
	  START_LOG("system(bzip2)", "Mesh");
	  if (libMesh::processor_id() == 0)
	    {
	      std::string system_string = "bzip2 -f ";
	      system_string += new_name;
	      std::system(system_string.c_str());
	    }
	  Parallel::barrier();
	  STOP_LOG("system(bzip2)", "Mesh");
	}
      
      STOP_LOG("write()", "Mesh");
    }  
}



void UnstructuredMesh::write (const std::string& name,
		  const std::vector<Number>& v,
		  const std::vector<std::string>& vn)
{
  START_LOG("write()", "Mesh");

  // Write the file based on extension
  if (name.rfind(".dat") < name.size())
    TecplotIO(*this).write_nodal_data (name, v, vn);
  
  else if (name.rfind(".plt") < name.size())
    TecplotIO(*this,true).write_nodal_data (name, v, vn);
  
  else if (name.rfind(".gmv") < name.size())
    {
      if (n_subdomains() > 1)
	GMVIO(*this).write_nodal_data (name, v, vn);
      else
	{
	  GMVIO io(*this);
	  io.partitioning() = false;
	  io.write_nodal_data (name, v, vn);
	}
    }    
  else if (name.rfind(".pvtu") < name.size())
    {
      VTKIO(*this).write_nodal_data (name, v, vn);
    }
  else
    {
      std::cerr << " ERROR: Unrecognized file extension: " << name
		<< "\n   I understand the following:\n\n"
		<< "     *.dat  -- Tecplot ASCII file\n"
		<< "     *.plt  -- Tecplot binary file\n"
		<< "     *.pvtu -- Paraview VTK file\n"
		<< "     *.gmv  -- LANL's GMV (General Mesh Viewer) format\n"
		<< "\n Exiting without writing output\n";
    }

  STOP_LOG("write()", "Mesh");
}





void UnstructuredMesh::create_pid_mesh(UnstructuredMesh& pid_mesh,
			   const unsigned int pid) const
{

  // Issue a warning if the number the number of processors
  // currently available is less that that requested for
  // partitioning.  This is not necessarily an error since
  // you may run on one processor and still partition the
  // mesh into several partitions.
#ifdef DEBUG
  if (this->n_processors() < pid)
    {
      std::cout << "WARNING:  You are creating a "
		<< "mesh for a processor id (="
		<< pid
		<< ") greater than "
		<< "the number of processors available for "
		<< "the calculation. (="
		<< libMesh::n_processors()
		<< ")."
		<< std::endl;
    }
#endif
  
  // Create iterators to loop over the list of elements
//   const_active_pid_elem_iterator       it(this->elements_begin(),   pid);
//   const const_active_pid_elem_iterator it_end(this->elements_end(), pid);

  const_element_iterator       it     = this->active_pid_elements_begin(pid);
  const const_element_iterator it_end = this->active_pid_elements_end(pid);
    
  this->create_submesh (pid_mesh, it, it_end);
}







void UnstructuredMesh::create_submesh (UnstructuredMesh& new_mesh,
			   const_element_iterator& it,
			   const const_element_iterator& it_end) const
{
  // Just in case the subdomain_mesh already has some information
  // in it, get rid of it.
  new_mesh.clear();

  // Fail if (*this == new_mesh), we cannot create a submesh inside ourself!
  // This may happen if the user accidently passes the original mesh into
  // this function!  We will check this by making sure we did not just
  // clear ourself.
  libmesh_assert (this->n_nodes() != 0);
  libmesh_assert (this->n_elem()  != 0); 

  // How the nodes on this mesh will be renumbered to nodes
  // on the new_mesh.  
  std::vector<unsigned int> new_node_numbers (this->n_nodes());

  std::fill (new_node_numbers.begin(),
	     new_node_numbers.end(),
	     libMesh::invalid_uint);

  
  
  // the number of nodes on the new mesh, will be incremented
  unsigned int n_new_nodes = 0;
  unsigned int n_new_elem  = 0;
    
  for (; it != it_end; ++it)
    {
      // increment the new element counter
      n_new_elem++;
	
      const Elem* old_elem = *it;

      // Add an equivalent element type to the new_mesh
      Elem* new_elem = 
	new_mesh.add_elem (Elem::build(old_elem->type()).release());

      libmesh_assert (new_elem != NULL);
	
      // Loop over the nodes on this element.  
      for (unsigned int n=0; n<old_elem->n_nodes(); n++)
	{
	  libmesh_assert (old_elem->node(n) < new_node_numbers.size());

	  if (new_node_numbers[old_elem->node(n)] == libMesh::invalid_uint)
	    {
	      new_node_numbers[old_elem->node(n)] = n_new_nodes;

	      // Add this node to the new mesh
	      new_mesh.add_point (old_elem->point(n));

	      // Increment the new node counter
	      n_new_nodes++;
	    }

	  // Define this element's connectivity on the new mesh
	  libmesh_assert (new_node_numbers[old_elem->node(n)] < new_mesh.n_nodes());
	    
	  new_elem->set_node(n) = new_mesh.node_ptr (new_node_numbers[old_elem->node(n)]);
	}

      // Maybe add boundary conditions for this element
      for (unsigned int s=0; s<old_elem->n_sides(); s++)
	if (old_elem->neighbor(s) == NULL)
	  if (this->boundary_info->boundary_id (old_elem, s) !=
	      this->boundary_info->invalid_id)
	    new_mesh.boundary_info->add_side (new_elem,
					     s,
					     this->boundary_info->boundary_id (old_elem, s));
    } // end loop over elements
  

  // Prepare the new_mesh for use
  new_mesh.prepare_for_use();
  
}



#ifdef ENABLE_AMR
bool UnstructuredMesh::contract ()
{
  START_LOG ("contract()", "Mesh");

  // Flag indicating if this call actually changes the mesh
  bool mesh_changed = false;

  element_iterator in        = elements_begin();
  element_iterator out       = elements_begin();
  const element_iterator end = elements_end();

#ifdef DEBUG
  for ( ; in != end; ++in)
    if (*in != NULL)
      {
	Elem* elem = *in;
	libmesh_assert(elem->active() || elem->subactive() || elem->ancestor());
      }
  in = elements_begin();
#endif
  
  // Loop over the elements.   
  for ( ; in != end; ++in)
    if (*in != NULL)
      {
	Elem* elem = *in;

	// Delete all the subactive ones
	if (elem->subactive())
	  {
	    // Huh?  no level-0 element should be subactive
	    libmesh_assert (elem->level() != 0);

	    // Delete the element
	    // This just sets a pointer to NULL, and doesn't
	    // invalidate any iterators
	    this->delete_elem(elem);
	    
	    // the mesh has certainly changed
	    mesh_changed = true;
	  }
	else
	  {
	    // Compress all the active ones
	    if (elem->active())
	      elem->contract();
	    else
	      libmesh_assert (elem->ancestor());
	  }
      }

  // Strip any newly-created NULL voids out of the element array
  this->renumber_nodes_and_elements();

  STOP_LOG ("contract()", "Mesh");
  
  return mesh_changed;
}
#endif // #ifdef ENABLE_AMR