gmv_io.c

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

  // a "fromfile" directive (?) But we currently don't make
  // any use of this feature in LibMesh.  Nonzero return val
  // from any function usually means an error has occurred.
  int ierr = GMV::gmvread_open_fromfileskip(const_cast<char*>(name.c_str()));
  if (ierr != 0)
    {
      std::cerr << "GMV::gmvread_open_fromfileskip failed!" << std::endl;
      libmesh_error();
    }

  
  // Loop through file until GMVEND.  
  int iend = 0;
  while (iend == 0)
    {
      GMV::gmvread_data();

      /*  Check for GMVEND.  */
      if (GMV::gmv_data.keyword == GMVEND)
        {
	  iend = 1;
	  GMV::gmvread_close();
	  break;
        }

      /*  Check for GMVERROR.  */
      if (GMV::gmv_data.keyword == GMVERROR)
        {
	  std::cerr << "Encountered GMVERROR while reading!" << std::endl;
	  libmesh_error();
        }

      /*  Process the data.  */
      switch (GMV::gmv_data.keyword)
        {
	case NODES:
	  {
	    //std::cout << "Reading nodes." << std::endl;

	    if (GMV::gmv_data.num2 == NODES)
	      this->_read_nodes();
	    
	    else if (GMV::gmv_data.num2 == NODE_V)
	      {
		std::cerr << "Unsupported GMV data type NODE_V!" << std::endl;
		libmesh_error();
	      }
	    break;
	  }
	  
	case CELLS:
	  {
	    // Read 1 cell at a time
	    // std::cout << "\nReading one cell." << std::endl;
	    this->_read_one_cell();
	    break;
	  }

	case MATERIAL:
	  {
	    // keyword == 6
	    // These are the materials, which we use to specify the mesh
	    // partitioning.
	    this->_read_materials();
	    break;
	  }
	  
	case VARIABLE:
	  {
	    // keyword == 8
	    // This is a field variable.  

	    // Check to see if we're done reading variables and break out.
	    if (GMV::gmv_data.datatype == ENDKEYWORD)
	      {
		// std::cout << "Done reading GMV variables." << std::endl;
		break;
	      }

	    if (GMV::gmv_data.datatype == NODE)
	      {
		// std::cout << "Reading node field data for variable "
		// 	  << GMV::gmv_data.name1 << std::endl;
		this->_read_var();
		break;
	      }
	    
	    else
	      {
		std::cerr << "Warning: Skipping variable: "
			  << GMV::gmv_data.name1
			  << " which is of unsupported GMV datatype "
			  << GMV::gmv_data.datatype
			  << ".  Nodal field data is currently the only type currently supported."
			  << std::endl;
		break;
	      }

	  }
	  
	default:
	  {
	    std::cerr << "Encountered unknown GMV keyword "
		      << GMV::gmv_data.keyword
		      << std::endl;
	    libmesh_error();
	  }
        } // end switch
    } // end while

  // Done reading in the mesh, now call find_neighbors, etc.
  // MeshInput<MeshBase>::mesh().find_neighbors();
  
  // Pass true flag to skip renumbering nodes and elements
  MeshInput<MeshBase>::mesh().prepare_for_use(true);
#endif
}




void GMVIO::_read_var()
{
#ifdef HAVE_GMV
  
  // Copy all the variable's values into a local storage vector.
  _nodal_data.insert ( std::make_pair(std::string(GMV::gmv_data.name1),
				      std::vector<Number>(GMV::gmv_data.doubledata1, GMV::gmv_data.doubledata1+GMV::gmv_data.num) ) );
#endif							 
}



void GMVIO::_read_materials()
{
#ifdef HAVE_GMV
  
  // LibMesh assigns materials on a per-cell basis
  libmesh_assert (GMV::gmv_data.datatype == CELL);

  //   // Material names: LibMesh has no use for these currently...
  //   std::cout << "Number of material names="
  // 	    << GMV::gmv_data.num
  // 	    << std::endl;
  
  //   for (int i = 0; i < GMV::gmv_data.num; i++)
  //     {
  //       // Build a 32-char string from the appropriate entries
  //       std::string mat_string(&GMV::gmv_data.chardata1[i*33], 32);

  //       std::cout << "Material name " << i << ": " << mat_string << std::endl;
  //     }

  //   // Material labels: These correspond to (1-based) CPU IDs, and
  //   // there should be 1 of these for each element.
  //   std::cout << "Number of material labels = "
  // 	    << GMV::gmv_data.nlongdata1
  // 	    << std::endl;
  
  for (int i = 0; i < GMV::gmv_data.nlongdata1; i++)
    {
      // Debugging Info
      // std::cout << "Material ID " << i << ": "
      // << GMV::gmv_data.longdata1[i]
      // << std::endl;
      
      MeshInput<MeshBase>::mesh().elem(i)->processor_id() =
	GMV::gmv_data.longdata1[i]-1;
    }
  
#endif
}




void GMVIO::_read_nodes()
{
#ifdef HAVE_GMV
  //   // Debug Info 
  //   std::cout << "gmv_data.datatype="
  // 	    <<  GMV::gmv_data.datatype
  // 	    << std::endl;

  // LibMesh writes UNSTRUCT=100 node data 
  libmesh_assert (GMV::gmv_data.datatype == UNSTRUCT);

  // The nodal data is stored in gmv_data.doubledata{1,2,3}
  // and is nnodes long
  for (int i = 0; i < GMV::gmv_data.num; i++)
    {
      //       std::cout << "(x,y,z)="
      // 		<< "("
      // 		<< GMV::gmv_data.doubledata1[i]
      // 		<< ","
      // 		<< GMV::gmv_data.doubledata2[i]
      // 		<< ","
      // 		<< GMV::gmv_data.doubledata3[i]
      // 		<< ")"
      // 		<< std::endl;
      
      // Add the point to the Mesh
      MeshInput<MeshBase>::mesh().add_point
                     ( Point(GMV::gmv_data.doubledata1[i],
			     GMV::gmv_data.doubledata2[i],
			     GMV::gmv_data.doubledata3[i]), i);
    }
#endif  
}


void GMVIO::_read_one_cell()
{
#ifdef HAVE_GMV
  //   // Debug Info 
  //   std::cout << "gmv_data.datatype="
  // 	    <<  GMV::gmv_data.datatype
  // 	    << std::endl;

  // This is either a REGULAR=111 cell or
  // the ENDKEYWORD=207 of the cells
#ifndef NDEBUG
    bool recognized =
      (GMV::gmv_data.datatype==REGULAR) ||
      (GMV::gmv_data.datatype==ENDKEYWORD);
#endif
    libmesh_assert (recognized);


  if (GMV::gmv_data.datatype == REGULAR)
    {
      //       std::cout << "Name of the cell is: "
      // 		<< GMV::gmv_data.name1
      // 		<< std::endl;

      //       std::cout << "Cell has "
      // 		<< GMV::gmv_data.num2
      // 		<< " vertices."
      // 		<< std::endl;

      // We need a mapping from GMV element types to LibMesh
      // ElemTypes.  Basically the reverse of the eletypes
      // std::map above.
      //
      // FIXME: Since Quad9's apparently don't exist for GMV, and since
      // In general we write linear sub-elements to GMV files, we need
      // to be careful to read back in exactly what we wrote out...
      ElemType type = this->_gmv_elem_to_libmesh_elem(GMV::gmv_data.name1);
      
      Elem* elem = Elem::build(type).release();
      elem->set_id(_next_elem_id++);

      // Print out the connectivity information for
      // this cell.
      for (int i = 0; i < GMV::gmv_data.num2; i++)
	{
	  // 	  // Debugging info
	  // 	  std::cout << "Vertex " << i << " is node "
	  // 		    << GMV::gmv_data.longdata1[i]
	  // 		    << std::endl;
	  
	  // Note: Node numbers are 1-based
	  elem->set_node(i) = MeshInput<MeshBase>::mesh().node_ptr(GMV::gmv_data.longdata1[i]-1);
	}

      // Add the newly-created element to the mesh
      MeshInput<MeshBase>::mesh().add_elem(elem);
    }


  if (GMV::gmv_data.datatype == ENDKEYWORD)
    {
      // There isn't a cell to read, so we just return
      return;
    }

#endif
}


ElemType GMVIO::_gmv_elem_to_libmesh_elem(const char* elemname)
{
  //
  // Linear Elements
  //
  if (!std::strncmp(elemname,"line",4))
    return EDGE2;
  
  if (!std::strncmp(elemname,"tri",3))
    return TRI3;

  if (!std::strncmp(elemname,"quad",4))
    return QUAD4;
  
  // FIXME: tet or ptet4?
  if ((!std::strncmp(elemname,"tet",3)) ||
      (!std::strncmp(elemname,"ptet4",5)))
    return TET4;

  // FIXME: hex or phex8?
  if ((!std::strncmp(elemname,"hex",3)) ||
      (!std::strncmp(elemname,"phex8",5)))
    return HEX8;

  // FIXME: prism or pprism6?
  if ((!std::strncmp(elemname,"prism",5)) ||
      (!std::strncmp(elemname,"pprism6",7)))
    return PRISM6;

  //
  // Quadratic Elements
  //
  if (!std::strncmp(elemname,"phex20",6))
    return HEX20;

  if (!std::strncmp(elemname,"phex27",6))
    return HEX27;

  if (!std::strncmp(elemname,"pprism15",8))
    return PRISM15;

  if (!std::strncmp(elemname,"ptet10",6))
    return TET10;

  if (!std::strncmp(elemname,"6tri",4))
    return TRI6;

  if (!std::strncmp(elemname,"8quad",5))
    return QUAD8;

  if (!std::strncmp(elemname,"3line",5))
    return EDGE3;
  
  // Unsupported/Unused types
  // if (!std::strncmp(elemname,"vface2d",7))
  // if (!std::strncmp(elemname,"vface3d",7))
  // if (!std::strncmp(elemname,"pyramid",7))
  // if (!std::strncmp(elemname,"ppyrmd5",7))
  // if (!std::strncmp(elemname,"ppyrmd13",8))

  // If we didn't return yet, then we didn't find the right cell!
  std::cerr << "Uknown/unsupported element: "
	    << elemname
	    << " was read."
	    << std::endl;
  libmesh_error();
}




void GMVIO::copy_nodal_solution(EquationSystems& es)
{
  // Check for easy return if there isn't any nodal data
  if (_nodal_data.empty())
    {
      std::cerr << "Unable to copy nodal solution: No nodal "
		<< "solution has been read in from file." << std::endl;
      return;
    }

  // Be sure there is at least one system
  libmesh_assert (es.n_systems());

  // Keep track of variable names which have been found and
  // copied already.  This could be used to prevent us from
  // e.g. copying the same var into 2 different systems ...
  // but this seems unlikely.  Also, it is used to tell if
  // any variables which were read in were not successfully
  // copied to the EquationSystems.
  std::set<std::string> vars_copied;
  
  // For each entry in the nodal data map, try to find a system
  // that has the same variable key name.
  for (unsigned int sys=0; sys<es.n_systems(); ++sys)
    {
      // Get a generic refernence to the current System
      System& system = es.get_system(sys);

      // And a reference to that system's dof_map
      // const DofMap & dof_map = system.get_dof_map();
	
      // For each var entry in the _nodal_data map, try to find
      // that var in the system
      std::map<std::string, std::vector<Number> >::iterator it = _nodal_data.begin();
      const std::map<std::string, std::vector<Number> >::iterator end = _nodal_data.end();
      for (; it != end; ++it)
	{
	  std::string var_name = (*it).first;
	  // std::cout << "Searching for var " << var_name << " in system " << sys << std::endl;

	  if (system.has_variable(var_name))
	    {
	      // Check if there are as many nodes in the mesh as there are entries
	      // in the stored nodal data vector
	      libmesh_assert ( (*it).second.size() == MeshInput<MeshBase>::mesh().n_nodes() );
	      
	      const unsigned int var_num = system.variable_number(var_name);
	      
	      // std::cout << "Variable "
	      // 			<< var_name
	      // 			<< " is variable "
	      // 			<< var_num 
	      // 			<< " in system " << sys << std::endl;

	      // The only type of nodal data we can read in from GMV is for
	      // linear LAGRANGE type elements.
	      const FEType& fe_type = system.variable_type(var_num);
	      if ((fe_type.order != FIRST) || (fe_type.family != LAGRANGE))
		{
		  std::cerr << "Only FIRST-order LAGRANGE variables can be read from GMV files. "
			    << "Skipping variable " << var_name << std::endl;
		  //libmesh_error();
		  break;
		}

	      
	      // Loop over the stored vector's entries, inserting them into
	      // the System's solution if appropriate.
	      for (unsigned int i=0; i<(*it).second.size(); ++i)
		{
		  // Since this var came from a GMV file, the index i corresponds to
		  // the (single) DOF value of the current variable for node i.
		  const unsigned int dof_index =
		    MeshInput<MeshBase>::mesh().node_ptr(i)->dof_number(sys,      /*system #*/
									var_num,  /*var # */
									0);       /*component #, always zero for LAGRANGE */

		  // std::cout << "Value " << i << ": "
		  // 			    << (*it).second [i]
		  // 			    << ", dof index="
		  // 			    << dof_index << std::endl;

		  // If the dof_index is local to this processor, set the value
		  if ((dof_index >= system.solution->first_local_index()) &&
		      (dof_index <  system.solution->last_local_index()))
		    system.solution->set (dof_index, (*it).second [i]);
		} // end loop over my GMVIO's copy of the solution

	      // Add the most recently copied var to the set of copied vars
	      vars_copied.insert (var_name);
	    } // end if (system.has_variable)
	} // end for loop over _nodal_data

      // Communicate parallel values before going to the next system.
      system.update();

    } // end loop over all systems


  
  // Warn the user if any GMV variables were not successfully copied over to the EquationSystems object
  {
    std::map<std::string, std::vector<Number> >::iterator it = _nodal_data.begin();
    const std::map<std::string, std::vector<Number> >::iterator end = _nodal_data.end();

    for (; it != end; ++it)
      {
	if (vars_copied.find( (*it).first ) == vars_copied.end())
	  {
	    std::cerr << "Warning: Variable "
		      << (*it).first
		      << " was not copied to the EquationSystems object."
		      <