legacy_xdr_io.c

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

            {
              elem = Elem::build(etypes[idx]).release();
              elem->set_id(_next_elem_id++);
              mesh.add_elem(elem);
            }
            
            // Add elements with the same id as in libMesh.  
            // Provided the data files that MeshData reads    
            // were only written with MeshData, then this      
            // should work properly.  This is an inline
            // function, so that for disabled MeshData, this
            // should not induce too much cost
            if (mesh_data != NULL)
              mesh_data->add_foreign_elem_id (elem, e);

            // Set the node pointers of the newly-created element
            for (unsigned int innd=0; innd < elem->n_nodes(); innd++)
            {
              elem->set_node(innd) = mesh.node_ptr(conn[innd+lastConnIndex]);
            }

            lastConnIndex += (m.get_orig_flag() == LegacyXdrIO::LIBM) ? (elem->n_nodes()+2) : elem->n_nodes();
          }
          lastFaceIndex += neeb[idx];
        }
        
      }

      if (m.get_orig_flag() == LegacyXdrIO::LIBM)
        {
          // Iterate in ascending elem ID order
          unsigned int _next_elem_id = 0;
          for (std::map<unsigned int, Elem *>::iterator i =
               parents.begin();
               i != parents.end(); ++i)
            {
              Elem *elem = i->second;
              if (elem)
                {
                  elem->set_id(_next_elem_id++);
                  mesh.add_elem(elem);
                }
              else
                // We can probably handle this, but we don't expect it
                libmesh_error();
            }
        }
    }
 
  // MGF-style (1) Hex27 mesh
  else if (m.get_orig_flag() == LegacyXdrIO::MGF) 
    {
      
#ifdef DEBUG
      if (mesh_data != NULL)
	if (mesh_data->active())
	  {
	    std::cerr << "ERROR: MeshData not implemented for MGF-style mesh."
		      << std::endl;
	    libmesh_error();
	  }
#endif
      
      for (int ielm=0; ielm < numElem; ++ielm)
	{
	  Elem* elem = new Hex27;
          elem->set_id(ielm);
	  mesh.add_elem(elem);
	  
	  for (int innd=0; innd < 27; ++innd)
	    elem->set_node(innd) = mesh.node_ptr(conn[innd+2+(27+2)*ielm]);	
	}
    }

  
  // tell the MeshData object that we are finished 
  // reading data
  if (mesh_data != NULL)
    mesh_data->close_foreign_id_maps ();
  
  // Free memory used in
  // the connectivity
  // vector.
  conn.clear();


  // If we are reading,
  // read in the BCs
  // from the mesh file,
  // otherwise write the
  // boundary conditions
  // if the BoundaryInfo
  // object exists.
  if (numBCs > 0)
    {
      std::vector<int> bcs(numBCs*3);

      // Read the BCs from the XDR file
      m.BC(&bcs[0], numBCs);
  
      // Add to the boundary_info
      for (int ibc=0; ibc < numBCs; ibc++)
	mesh.boundary_info->add_side(bcs[0+ibc*3], bcs[1+ibc*3], bcs[2+ibc*3]);
    }
}



void LegacyXdrIO::write_mesh (const std::string& name,
			      const LegacyXdrIO::FileFormat originator)
{
  // get a read-only reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();

  // n_levels is a parallel-only method
  parallel_only();
  const unsigned int          n_levels = MeshTools::n_levels(mesh);

  // The Legacy Xdr IO code only works if we have a serialized mesh
  libmesh_assert (mesh.is_serial());

  // In which case only processor 0 needs to do any writing
  if (libMesh::processor_id() != 0)
    return;
  
  // Create an XdrMESH object.
  XdrMESH m;

  // Create a pointer
  // to an XdrMESH file
  // header.
  XdrMHEAD mh;

  // Open the XDR file for writing.
  // Note 1: Provide an additional argument
  // to specify the dimension.
  //
  // Note 2: Has to do the right thing for
  // both binary and ASCII files.
  m.set_orig_flag(originator);

  // From here on, things depend
  // on whether we are reading or
  // writing!  First, we define
  // header variables that may
  // be read OR written.
  std::vector<unsigned int> neeb;
  std::vector<ElemType> etypes;
  

  int n_non_subactive = 0;
  int non_subactive_weight = 0;

  // This map will associate 
  // the distance from the beginning of the set
  // to each node ID with the node ID itself.
  std::map<unsigned int, unsigned int> node_map;

  {
    // For each non-subactive element:
    // 1.) Increment the number of non subactive elements
    // 2.) Accumulate the total weight
    // 3.) Add the node ids to a set of non subactive node ids 
    std::set<unsigned int> not_subactive_node_ids;
    MeshBase::const_element_iterator el = mesh.elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.elements_end();
    for( ; el != end_el; ++el)
    {
      Elem* elem = (*el);
      if(!elem->subactive())
      {
        n_non_subactive++;
        non_subactive_weight += elem->n_nodes();

        for (unsigned int n=0; n<elem->n_nodes(); ++n)
          not_subactive_node_ids.insert(elem->node(n));
      }
    }

    // Now that the set is built, most of the hard work is done.  We build
    // the map next and let the set go out of scope.
    std::set<unsigned int>::iterator it = not_subactive_node_ids.begin();
    const std::set<unsigned int>::iterator end = not_subactive_node_ids.end();
    unsigned int cnt=0;
    for (; it!=end; ++it)
      node_map[*it] = cnt++;
  }


  const int                   numElem  = n_non_subactive;       
  const int                   numBCs   = mesh.boundary_info->n_boundary_conds();
  
  // Fill the etypes vector with all of the element types found in the mesh
  MeshTools::elem_types(mesh, etypes);
  
  // store number of elements in each block at each refinement level
  neeb.resize((n_levels+1)*etypes.size()); 

  // Store a variable for the number of element types                   
  const unsigned int n_el_types = etypes.size();
	
  m.set_num_levels(n_levels);

  // The last argument is zero because mesh files are always number 0 ...
  m.init((this->binary() ? XdrMGF::ENCODE : XdrMGF::W_ASCII), name.c_str(), 0); 

  // Loop over all levels and all element types to set the entries of neeb
  for(unsigned int level=0; level<=n_levels; level++)
    for (unsigned int el_type=0; el_type<n_el_types; el_type++)
      neeb[level*n_el_types + el_type] = 
        MeshTools::n_non_subactive_elem_of_type_at_level(mesh, etypes[el_type], level);
        // gotta change this function name!!!


  // Now we check to see if we're doing
  // MGF-style headers or libMesh-style
  // "augmented" headers.  An
  // augmented header contains 
  // information about mesh blocks,
  // allowing us to optimize storage
  // and minimize IO requirements
  // for these meshes.
  if ((m.get_orig_flag() == LegacyXdrIO::DEAL) || (m.get_orig_flag() == LegacyXdrIO::LIBM))
    {
      mh.set_n_blocks(etypes.size());
      mh.set_block_elt_types(etypes);
      mh.set_num_elem_each_block(neeb);
    }
  else
    libmesh_assert(etypes.size() == 1);
  
  mh.setNumEl(numElem);
  mh.setNumNodes(node_map.size());
  mh.setStrSize(65536);
 
  // set a local variable for the total weight of the mesh
  int totalWeight =0;

  if (m.get_orig_flag() == LegacyXdrIO::DEAL)  // old-style LibMesh
    totalWeight=MeshTools::total_weight(mesh);

  else if (m.get_orig_flag() == LegacyXdrIO::MGF) // MGF-style
    totalWeight = MeshTools::total_weight(mesh)+2*numElem;

  else if (m.get_orig_flag() == LegacyXdrIO::LIBM) // new-style LibMesh
    totalWeight = non_subactive_weight+2*numElem;

  else
    libmesh_error();
    
  // Set the total weight in the header
  mh.setSumWghts(totalWeight);
        
  mh.setNumBCs(numBCs);
  mh.setId("Id String");       // You can put whatever you want, it will be ignored 
  mh.setTitle("Title String"); // You can put whatever you want, it will be ignored 
  
  // Put the information
  // in the XDR file.
  m.header(&mh); 
  
  
  // Write the connectivity  
  {
    std::vector<int> conn;
    LegacyXdrIO::FileFormat orig_type = m.get_orig_flag();
   
    // Resize the connectivity vector to hold all the connectivity for the mesh
    conn.resize(totalWeight);
    
    unsigned int lastConnIndex = 0;
    unsigned int nn = 0;
   
    // Loop over levels and types again, write connectivity information to conn.
    for (unsigned int level=0; level<=n_levels; level++)
      for (unsigned int idx=0; idx<etypes.size(); idx++)
      {
        nn = lastConnIndex = 0;

        for (unsigned int e=0; e<mesh.n_elem(); e++)
          if ((mesh.elem(e)->type()  == etypes[idx]) && 
              (mesh.elem(e)->level() == level)       &&
              !mesh.elem(e)->subactive())
          {
            int nstart=0;
            
            if (orig_type == LegacyXdrIO::DEAL)
              nn = mesh.elem(e)->n_nodes();

            else if (orig_type == LegacyXdrIO::MGF)
            {
              nstart=2; // ignore the 27 and 0 entries
              nn = mesh.elem(e)->n_nodes()+2;
              conn[lastConnIndex + 0] = 27;
              conn[lastConnIndex + 1] = 0;
            }

            else if (orig_type == LegacyXdrIO::LIBM) // LIBMESH format
              nn = mesh.elem(e)->n_nodes() + 2;

            else
              libmesh_error();

            // Loop over the connectivity entries for this element and write to conn.
            START_LOG("set connectivity", "LegacyXdrIO::write_mesh");
            const unsigned int loopmax = (orig_type==LegacyXdrIO::LIBM) ? nn-2 : nn;
            for (unsigned int n=nstart; n<loopmax; n++)
            {
              unsigned int connectivity_value=0;

              // old-style Libmesh and MGF meshes
              if (orig_type != LegacyXdrIO::LIBM)
                connectivity_value = mesh.elem(e)->node(n-nstart);

              // new-style libMesh meshes: compress the connectivity entries to account for
              // subactive nodes that will not be in the mesh we write out.
              else
              {
                std::map<unsigned int, unsigned int>::iterator pos = 
                  node_map.find(mesh.elem(e)->node(n-nstart));

                libmesh_assert (pos != node_map.end());

                connectivity_value = (*pos).second;
              }
              conn[lastConnIndex + n] = connectivity_value;
            }
            STOP_LOG("set connectivity", "LegacyXdrIO::write_mesh");

            // In the case of an adaptive mesh, set last 2 entries to this ID and parent ID
            if (orig_type == LegacyXdrIO::LIBM)
            {
              int self_ID = mesh.elem(e)->id();
              int parent_ID = -1;
              if(level != 0)
                parent_ID = mesh.elem(e)->parent()->id();

              // Self ID is the second-to-last entry, Parent ID is the last
              // entry on each connectivity line
              conn[lastConnIndex+nn-2] = self_ID;
              conn[lastConnIndex+nn-1] = parent_ID;
            }

            lastConnIndex += nn;
          }

        // Send conn to the XDR file.  If there are no elements of this level and type,
        // then nn will be zero, and we there is no connectivity to write. 
        if (nn != 0)
          m.Icon(&conn[0], nn, lastConnIndex/nn);
      }
  }
    
  // create the vector of coords and send
  // it to the XDR file.
  {
    std::vector<Real> coords;
    
    coords.resize(mesh.spatial_dimension()*node_map.size());
    int lastIndex=0;

    std::map<unsigned int,unsigned int>::iterator it = node_map.begin();
    const std::map<unsigned int,unsigned int>::iterator end = node_map.end();
    for (; it != end; ++it)
      {
        const Point& p = mesh.node((*it).first);

        coords[lastIndex+0] = p(0);
        coords[lastIndex+1] = p(1);
        coords[lastIndex+2] = p(2);
        lastIndex += 3;
      }
   
    // Put the nodes in the XDR file
    m.coord(&coords[0], mesh.spatial_dimension(), node_map.size()); 
  }

  
  // write the
  // boundary conditions
  // if the BoundaryInfo
  // object exists.
  if (numBCs > 0)
    {
      std::vector<int> bcs(numBCs*3);
    
      //std::cout << "numBCs=" << numBCs << std::endl;
    
      //std::cout << "Preparing to write boundary conditions." << std::endl;
      std::vector<unsigned int> elem_list;
      std::vector<unsigned short int> side_list;
      std::vector<short int> elem_id_list;
      
      mesh.boundary_info->build_side_list (elem_list, side_list, elem_id_list);
    
      for (int ibc=0;  ibc<numBCs; ibc++)
	{
	  bcs[0+ibc*3] = elem_list[ibc];
	  bcs[1+ibc*3] = side_list[ibc];
	  bcs[2+ibc*3] = elem_id_list[ibc];
	}
    
      // Put the BCs in the XDR file
      m.BC(&bcs[0], numBCs);
    }
}



void LegacyXdrIO::read_soln (const std::string& name,
			     std::vector<Real>& soln,
			     std::vector<std::string>& var_names) const
{
  // Create an XdrSOLN object.
  XdrSOLN s;
  
  // Create an XdrSHEAD object.
  XdrSHEAD sh;
  
  // Open the XDR file for
  // reading or writing.
  // Note 1: Provide an additional argument
  // to specify the dimension.
  //
  // Note 2: Has to do the right thing for
  // both binary and ASCII files.
  s.init((this->binary() ? XdrMGF::DECODE : XdrMGF::R_ASCII), name.c_str(), 0); // mesh files are always number 0 ...
  
  // From here on, things depend
  // on whether we are reading or
  // writing!  First, we define
  // header variables that may
  // be read OR written.
  int         numVar      = 0;       
  int         numNodes    = 0;
  const char* varNames;
	
  // Get the information from
  // the header, and place it
  // in the header pointer.
  s.header(&sh);
	
  // Read information from the
  // file header.  This depends on
  // whether its a libMesh or MGF mesh.
  numVar   = sh.getWrtVar();
  numNodes = sh.getNumNodes();
  varNames = sh.getVarTitle();
	
  // Get the variable names
  {	  
    var_names.resize(numVar);
    
    const char* p = varNames;
    
    for (int i=0; i<numVar; i++)
      {
	var_names[i] = p;
	p += std::strlen(p) + 1;
      }
  }
  
  // Read the soln vector
  soln.resize(numVar*numNodes);
    
  s.values(&soln[0], numNodes);	
}  



void LegacyXdrIO::write_soln (const std::string& name,
			      std::vector<Real>& soln,
			      std::vector<std::string>& var_names) const
{
  // get a read-only reference to the mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();
  
  // Create an XdrSOLN object.
  XdrSOLN s;
  
  // Create an XdrSHEAD object.
  XdrSHEAD sh;
  
  // Open the XDR file for
  // reading or writing.
  // Note 1: Provide an additional argument
  // to specify the dimension.
  //
  // Note 2: Has to do the right thing for
  // both binary and ASCII files.
  s.init((this->binary() ? XdrMGF::ENCODE : XdrMGF::W_ASCII), name.c_str(), 0); // mesh files are always number 0 ...

  // Build the header
  sh.setWrtVar(var_names.size());
  sh.setNumVar(var_names.size());
  sh.setNumNodes(mesh.n_nodes());
  sh.setNumBCs(mesh.boundary_info->n_boundary_conds());
  sh.setMeshCnt(0);
  sh.setKstep(0);
  sh.setTime(0.);
  sh.setStrSize(65536);
  sh.setId("Id String");	               // Ignored
  sh.setTitle("Title String");          // Ignored
  sh.setUserTitle("User Title String"); // Ignored
  
  // create the variable array
  {
    std::string var_title;
    
    for (unsigned int var=0; var<var_names.size(); var++)
      {
	for (unsigned int c=0; c<var_names[var].size(); c++)
	  var_title += var_names[var][c];
	
	var_title += '\0';
      }
    
    sh.setVarTitle(var_title.c_str(), var_title.size());
  }
  
  // Put the informationin the XDR file.
  s.header(&sh); // Needs to work for both types of file
  
  // Write the solution vector
  libmesh_assert (soln.size() == var_names.size()*mesh.n_nodes());
  
  s.values(&soln[0], mesh.n_nodes());
}