unv_io.c

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

  // this was our first node
  this->_n_nodes++;



  // proceed _counting_ the remaining
  // nodes.
  while (in_file.good())
    {
      // read the node label
      in_file >> data; 
	   	
      if (data == "-1")
	// end of dataset is reached
	break;
      
      // ignore the remaining data (coord_sys_label, color etc)
      in_file.ignore (256, '\n');
      // ignore the coordinates
      in_file.ignore (256, '\n');

      this->_n_nodes++;
    }


  if (in_file.eof())
    {
      std::cerr << "ERROR: File ended before end of node dataset!"
		<< std::endl;
      libmesh_error();
    }

  if (this->verbose())
    std::cout << "  Nodes   : " << this->_n_nodes << std::endl;

  STOP_LOG("count_nodes()","UNVIO");
}






void UNVIO::count_elements (std::istream& in_file)
{
  START_LOG("count_elements()","UNVIO");

  if (this->_n_elements != 0)
    {
      std::cerr << "Error: Trying to scan elements twice!" 
		<< std::endl;
      libmesh_error();
    }


  // Simply read the element
  // dataset for the @e only
  // purpose to count nodes!
  
  std::string data;
  unsigned int fe_id;

  while (!in_file.eof())
    {
      // read element label
      in_file >> data;
      
      // end of dataset?
      if (data == "-1") 
	break;
	
      // read fe_id
      in_file >> fe_id;

      // Skip related data,
      // and node number list
      in_file.ignore (256,'\n');
      in_file.ignore (256,'\n');

      // For some elements the node numbers
      // are given more than one record

      // TET10 or QUAD9
      if (fe_id == 118 || fe_id == 300)
	  in_file.ignore (256,'\n');

      // HEX20
      if (fe_id == 116)
	{
	  in_file.ignore (256,'\n');
	  in_file.ignore (256,'\n');
	}
     
      this->_n_elements++;
    }


  if (in_file.eof())
    {
      std::cerr << "ERROR: File ended before end of element dataset!"
		<< std::endl;
      libmesh_error();
    }

  if (this->verbose())
    std::cout << "  Elements: " << this->_n_elements << std::endl;

  STOP_LOG("count_elements()","UNVIO");
}



void UNVIO::node_in (std::istream& in_file)
{
  START_LOG("node_in()","UNVIO");

  if (this->verbose())
    std::cout << "  Reading nodes" << std::endl;

  // adjust the \p istream to our position
  const bool ok = this->beginning_of_dataset(in_file, _label_dataset_nodes);

  if (!ok)
    {
      std::cerr << "ERROR: Could not find node dataset!" << std::endl;
      libmesh_error();
    }

  MeshBase& mesh = MeshInput<MeshBase>::mesh();
  
  unsigned int node_lab;           // label of the node
  unsigned int exp_coord_sys_num,  // export coordinate system number       (not supported yet)
               disp_coord_sys_num, // displacement coordinate system number (not supported yet)
               color;              // color                                 (not supported yet)

  // allocate the correct amount 
  // of memory for the node vector
  this->_assign_nodes.reserve (this->_n_nodes);


  // always 3 coordinates in the UNV file, no matter
  // which dimensionality libMesh is in
  //std::vector<Real> xyz (3);
  Point xyz;
      
  // depending on whether we have to convert each
  // coordinate (float), we offer two versions.
  // Note that \p count_nodes() already verified
  // whether this file uses "D" of "e"
  if (this->_need_D_to_e)
    {
      // ok, convert...
      std::string num_buf;

      for(unsigned int i=0; i<this->_n_nodes; i++)
        {
	  libmesh_assert (!in_file.eof());

	  in_file >> node_lab                // read the node label
		  >> exp_coord_sys_num       // (not supported yet)
		  >> disp_coord_sys_num      // (not supported yet)
		  >> color;                  // (not supported yet)

	  // take care of the
	  // floating-point data
	  for (unsigned int d=0; d<3; d++)
	    {
	      in_file >> num_buf;
	      xyz(d) = this->D_to_e (num_buf);
	    }

	  // set up the id map
	  this->_assign_nodes.push_back (node_lab);
	  
	  // add node to the Mesh & 
	  // tell the MeshData object the foreign node id
	  // (note that mesh.add_point() returns a pointer to the new node)
	  this->_mesh_data.add_foreign_node_id (mesh.add_point(xyz,i), node_lab);
	}
    }

  else
    {
      // very well, no need to convert anything,
      // just plain import.
      for (unsigned int i=0;i<this->_n_nodes;i++)
        {
	  libmesh_assert (!in_file.eof());

	  in_file >> node_lab                // read the node label
		  >> exp_coord_sys_num       // (not supported yet)
		  >> disp_coord_sys_num      // (not supported yet)
		  >> color                   // (not supported yet)
		  >> xyz(0)                  // read x-coordinate
		  >> xyz(1)                  // read y-coordinate
		  >> xyz(2);                 // read z-coordinate

	  // set up the id map
	  this->_assign_nodes.push_back (node_lab);
	  
	  // add node to the Mesh & 
	  // tell the MeshData object the foreign node id
	  // (note that mesh.add_point() returns a pointer to the new node)
	  this->_mesh_data.add_foreign_node_id (mesh.add_point(xyz,i), node_lab);
	}
    }

  // now we need to sort the _assign_nodes vector so we can
  // search it efficiently like a map
  std::sort (this->_assign_nodes.begin(),
	     this->_assign_nodes.end());

  STOP_LOG("node_in()","UNVIO");
}





void UNVIO::element_in (std::istream& in_file)
{
  START_LOG("element_in()","UNVIO");

  if (this->verbose())
    std::cout << "  Reading elements" << std::endl;

  MeshBase& mesh = MeshInput<MeshBase>::mesh();

  // adjust the \p istream to our
  // position
  const bool ok = this->beginning_of_dataset(in_file, _label_dataset_elements);

  if (!ok)
    {
      std::cerr << "ERROR: Could not find element dataset!" << std::endl;
      libmesh_error();
    }


  unsigned int      element_lab,       // element label (not supported yet)
                    n_nodes;           // number of nodes on element
  unsigned long int fe_descriptor_id,  // FE descriptor id
                    phys_prop_tab_num, // physical property table number (not supported yet)
                    mat_prop_tab_num,  // material property table number (not supported yet)
                    color;             // color (not supported yet)


  // vector that temporarily holds the node labels defining element
  std::vector<unsigned int> node_labels (21);


  // vector that assigns element nodes to their correct position
  // for example:
  // 44:plane stress      | QUAD4
  // linear quadrilateral |
  // position in UNV-file | position in libmesh
  // assign_elem_node[1]   = 0
  // assign_elem_node[2]   = 3
  // assign_elem_node[3]   = 2
  // assign_elem_node[4]   = 1
  //
  // UNV is 1-based, we leave the 0th element of the vectors unused in order
  // to prevent confusion, this way we can store elements with up to 20 nodes
  unsigned int assign_elem_nodes[21];
  

  // Get the beginning and end of the _assign_nodes vector
  // to eliminate repeated function calls
  const std::vector<unsigned int>::const_iterator it_begin =
    this->_assign_nodes.begin();
  
  const std::vector<unsigned int>::const_iterator it_end   =
    this->_assign_nodes.end();


      
  // read from the virtual file
  for (unsigned int i=0; i<this->_n_elements; i++)
    {
      in_file >> element_lab             // read element label
	      >> fe_descriptor_id        // read FE descriptor id
	      >> phys_prop_tab_num       // (not supported yet)
	      >> mat_prop_tab_num        // (not supported yet)
	      >> color                   // (not supported yet)
	      >> n_nodes;                // read number of nodes on element

      for (unsigned int j=1; j<=n_nodes; j++)
	in_file >> node_labels[j];       // read node labels

      Elem* elem = NULL;                 // element pointer
      
      switch (fe_descriptor_id)
	{
	
	case 41: // Plane Stress Linear Triangle
	case 91: // Thin Shell   Linear Triangle
	  {
	    elem = new Tri3;  // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=2;
	    assign_elem_nodes[3]=1;
	    break;
	  }
	
	case 42: // Plane Stress Quadratic Triangle
	case 92: // Thin Shell   Quadratic Triangle
	  {
	    elem = new Tri6;  // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=5;
	    assign_elem_nodes[3]=2;
	    assign_elem_nodes[4]=4;
	    assign_elem_nodes[5]=1;
	    assign_elem_nodes[6]=3;
	    break;
	  }
	
	case 43: // Plane Stress Cubic Triangle
	  {
	    std::cerr << "ERROR: UNV-element type 43: Plane Stress Cubic Triangle"
		      << " not supported." 
		      << std::endl;
	    libmesh_error();	    
	    break;
	  }
	
	case 44: // Plane Stress Linear Quadrilateral
	case 94: // Thin Shell   Linear Quadrilateral
	  {
	    elem = new Quad4; // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=3;
	    assign_elem_nodes[3]=2;
	    assign_elem_nodes[4]=1;
	    break;
	  }
	
	case 45: // Plane Stress Quadratic Quadrilateral
	case 95: // Thin Shell   Quadratic Quadrilateral
	  {
	    elem = new Quad8; // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=7;
	    assign_elem_nodes[3]=3;
	    assign_elem_nodes[4]=6;
	    assign_elem_nodes[5]=2;
	    assign_elem_nodes[6]=5;
	    assign_elem_nodes[7]=1;
	    assign_elem_nodes[8]=4;
	    break;
	  }
	
	case 300: // Thin Shell   Quadratic Quadrilateral (nine nodes)
	  {
	    elem = new Quad9; // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=7;
	    assign_elem_nodes[3]=3;
	    assign_elem_nodes[4]=6;
	    assign_elem_nodes[5]=2;
	    assign_elem_nodes[6]=5;
	    assign_elem_nodes[7]=1;
	    assign_elem_nodes[8]=4;
	    assign_elem_nodes[9]=8;
	    break;
	  }
	
	case 46: // Plane Stress Cubic Quadrilateral
	  {
	    std::cerr << "ERROR: UNV-element type 46: Plane Stress Cubic Quadrilateral"
		      << " not supported." 
		      << std::endl;
	    libmesh_error();
	    break;
	  }
	
	case 111: // Solid Linear Tetrahedron
	  {
	    elem = new Tet4;  // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=1;
	    assign_elem_nodes[3]=2;
	    assign_elem_nodes[4]=3;
	    break;
	  }

	case 112: // Solid Linear Prism
	  {
	    elem = new Prism6;  // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=1;
	    assign_elem_nodes[3]=2;
	    assign_elem_nodes[4]=3;
	    assign_elem_nodes[5]=4;
	    assign_elem_nodes[6]=5;
	    break;
	  }
	
	case 115: // Solid Linear Brick
	  {
	    elem = new Hex8;  // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=4;
	    assign_elem_nodes[3]=5;
	    assign_elem_nodes[4]=1;
	    assign_elem_nodes[5]=3;
	    assign_elem_nodes[6]=7;
	    assign_elem_nodes[7]=6;
	    assign_elem_nodes[8]=2;
	    break;
	  }