unv_io.c

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

	case 116: // Solid Quadratic Brick
	  {
	    elem = new Hex20; // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=12;
	    assign_elem_nodes[3]=4;
	    assign_elem_nodes[4]=16;
	    assign_elem_nodes[5]=5;
	    assign_elem_nodes[6]=13;
	    assign_elem_nodes[7]=1;
	    assign_elem_nodes[8]=8;

	    assign_elem_nodes[9]=11;
	    assign_elem_nodes[10]=19;
	    assign_elem_nodes[11]=17;
	    assign_elem_nodes[12]=9;

	    assign_elem_nodes[13]=3;
	    assign_elem_nodes[14]=15;
	    assign_elem_nodes[15]=7;
	    assign_elem_nodes[16]=18;
	    assign_elem_nodes[17]=6;
	    assign_elem_nodes[18]=14;
	    assign_elem_nodes[19]=2;
	    assign_elem_nodes[20]=10;
	    break;
	  }	
	
	case 117: // Solid Cubic Brick
	  {
	    std::cerr << "Error: UNV-element type 117: Solid Cubic Brick"
		      << " not supported." 
		      << std::endl;
	    libmesh_error();	    
	    break;
	  }
	
	case 118: // Solid Quadratic Tetrahedron
	  {
	    elem = new Tet10; // create new element
	    
	    assign_elem_nodes[1]=0;
	    assign_elem_nodes[2]=4;
	    assign_elem_nodes[3]=1;
	    assign_elem_nodes[4]=5;
	    assign_elem_nodes[5]=2;
	    assign_elem_nodes[6]=6;
	    assign_elem_nodes[7]=7;
	    assign_elem_nodes[8]=8;
	    assign_elem_nodes[9]=9;
	    assign_elem_nodes[10]=3;
	    break;
	  }	
	
	default: // Unrecognized element type
	  {
	    std::cerr << "ERROR: UNV-element type " 
		      << fe_descriptor_id
		      << " not supported." 
		      << std::endl;
	    libmesh_error();	    
	    break;
	  }
	}

      // nodes are being stored in element
      for (unsigned int j=1; j<=n_nodes; j++)
	{
	  // Find the position of node_labels[j] in the _assign_nodes vector.
	  const std::pair<std::vector<unsigned int>::const_iterator,
	                  std::vector<unsigned int>::const_iterator>	    
	    it = std::equal_range (it_begin,
				   it_end,
				   node_labels[j]);

	  // it better be there, so libmesh_assert that it was found.
	  libmesh_assert (it.first  != it.second);
	  libmesh_assert (*(it.first) == node_labels[j]);

	  // Now, the distance between this UNV id and the beginning of
	  // the _assign_nodes vector will give us a unique id in the
	  // range [0,n_nodes) that we can use for defining a contiguous
	  // connectivity.
	  const unsigned int assigned_node = std::distance (it_begin,
							    it.first);

	  // Make sure we didn't get an out-of-bounds id
	  libmesh_assert (assigned_node < this->_n_nodes);
	  
	  elem->set_node(assign_elem_nodes[j]) =
	    mesh.node_ptr(assigned_node);
	}
      
      // add elem to the Mesh & 
      // tell the MeshData object the foreign elem id
      // (note that mesh.add_elem() returns a pointer to the new element)
      elem->set_id(i);
      this->_mesh_data.add_foreign_elem_id (mesh.add_elem(elem), element_lab);
    }

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






void UNVIO::node_out (std::ostream& out_file)
{
  
  libmesh_assert (this->_mesh_data.active() ||
	  this->_mesh_data.compatibility_mode());


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

  // Write beginning of dataset
  out_file << "    -1\n"
	   << "  " 
	   << _label_dataset_nodes
	   << '\n';


  unsigned int exp_coord_sys_dummy  = 0; // export coordinate sys. (not supported yet)
  unsigned int disp_coord_sys_dummy = 0; // displacement coordinate sys. (not supp. yet)
  unsigned int color_dummy          = 0; // color(not supported yet)

  // A reference to the parent class's mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();

  MeshBase::const_node_iterator       nd  = mesh.nodes_begin();
  const MeshBase::const_node_iterator end = mesh.nodes_end();

  for (; nd != end; ++nd)
    {
      const Node* current_node = *nd;
      
      char buf[78];
      std::sprintf(buf, "%10d%10d%10d%10d\n", 
		   this->_mesh_data.node_to_foreign_id(current_node),
		   exp_coord_sys_dummy,
		   disp_coord_sys_dummy,
		   color_dummy);
      out_file << buf;

      // the coordinates
      if (mesh.spatial_dimension() == 3)
	std::sprintf(buf, "%25.16E%25.16E%25.16E\n", 
		     static_cast<double>((*current_node)(0)),
		     static_cast<double>((*current_node)(1)),
		     static_cast<double>((*current_node)(2)));
      else if (mesh.spatial_dimension() == 2)
	std::sprintf(buf, "%25.16E%25.16E\n", 
		     static_cast<double>((*current_node)(0)),
		     static_cast<double>((*current_node)(1)));
      else
	std::sprintf(buf, "%25.16E\n", 
		     static_cast<double>((*current_node)(0)));
      
      out_file << buf;
    }
  
  
  // Write end of dataset
  out_file << "    -1\n";
}






void UNVIO::element_out(std::ostream& out_file)
{
  libmesh_assert (this->_mesh_data.active() ||
	  this->_mesh_data.compatibility_mode());

  if (this->verbose())
    std::cout << "  Writing elements" << std::endl;
  
  // Write beginning of dataset
  out_file << "    -1\n"
	   << "  " 
	   << _label_dataset_elements
	   << "\n";

  unsigned long int fe_descriptor_id = 0;    // FE descriptor id 
  unsigned long int phys_prop_tab_dummy = 2; // physical property (not supported yet)
  unsigned long int mat_prop_tab_dummy = 1;  // material property (not supported yet)
  unsigned long int color_dummy = 7;         // color (not supported yet)


  // vector that assigns element nodes to their correct position
  // currently only elements with up to 20 nodes
  //  
  // Example:
  // QUAD4               | 44:plane stress
  //                     | linear quad
  // position in libMesh | UNV numbering
  // (note: 0-based)     | (note: 1-based)
  //     
  // assign_elem_node[0]  = 0
  // assign_elem_node[1]  = 3
  // assign_elem_node[2]  = 2
  // assign_elem_node[3]  = 1
  unsigned int assign_elem_nodes[20];

  unsigned int n_elem_written=0;

  // A reference to the parent class's mesh
  const MeshBase& mesh = MeshOutput<MeshBase>::mesh();

  MeshBase::const_element_iterator it  = mesh.elements_begin();
  const MeshBase::const_element_iterator end = mesh.elements_end();

  for (; it != end; ++it)
    {
      const Elem* elem = *it;

      elem->n_nodes();
      
      switch (elem->type())
	{
	
	case TRI3:
	  {
	    fe_descriptor_id = 41; // Plane Stress Linear Triangle
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 2;
	    assign_elem_nodes[2] = 1;
	    break;
	  }

	case TRI6:
	  {
	    fe_descriptor_id = 42; // Plane Stress Quadratic Triangle
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 5;
	    assign_elem_nodes[2] = 2;
	    assign_elem_nodes[3] = 4;
	    assign_elem_nodes[4] = 1;
	    assign_elem_nodes[5] = 3;
	    break;
	  }
	  
	case QUAD4:
	  {
	    fe_descriptor_id = 44; // Plane Stress Linear Quadrilateral
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 3;
	    assign_elem_nodes[2] = 2;
	    assign_elem_nodes[3] = 1;
	    break;
	  }
	
	case QUAD8:
	  {
	    fe_descriptor_id = 45; // Plane Stress Quadratic Quadrilateral
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 7;
	    assign_elem_nodes[2] = 3;
	    assign_elem_nodes[3] = 6;
	    assign_elem_nodes[4] = 2;
	    assign_elem_nodes[5] = 5;
	    assign_elem_nodes[6] = 1;
	    assign_elem_nodes[7] = 4;
	    break;
	  }
	
	case QUAD9:
	  {
	    fe_descriptor_id = 300; // Plane Stress Quadratic Quadrilateral
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 7;
	    assign_elem_nodes[2] = 3;
	    assign_elem_nodes[3] = 6;
	    assign_elem_nodes[4] = 2;
	    assign_elem_nodes[5] = 5;
	    assign_elem_nodes[6] = 1;
	    assign_elem_nodes[7] = 4;
	    assign_elem_nodes[8] = 8;
	    break;
	  }
	
	case TET4:
	  {
	    fe_descriptor_id = 111; // Solid Linear Tetrahedron
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 1;
	    assign_elem_nodes[2] = 2;
	    assign_elem_nodes[3] = 3;
	    break;
	  }

	case PRISM6:
	  {
	    fe_descriptor_id = 112; // Solid Linear Prism
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 1;
	    assign_elem_nodes[2] = 2;
	    assign_elem_nodes[3] = 3;
	    assign_elem_nodes[4] = 4;
	    assign_elem_nodes[5] = 5;
	    break;
	  }

	case HEX8:
	  {
	    fe_descriptor_id = 115; // Solid Linear Brick
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 4;
	    assign_elem_nodes[2] = 5;
	    assign_elem_nodes[3] = 1;
	    assign_elem_nodes[4] = 3;
	    assign_elem_nodes[5] = 7;
	    assign_elem_nodes[6] = 6;
	    assign_elem_nodes[7] = 2;
	    break;
	  }
	
	case HEX20:
	  {
	    fe_descriptor_id = 116; // Solid Quadratic Brick
	    assign_elem_nodes[ 0] = 0;
	    assign_elem_nodes[ 1] = 12;
	    assign_elem_nodes[ 2] = 4;
	    assign_elem_nodes[ 3] = 16;
	    assign_elem_nodes[ 4] = 5;
	    assign_elem_nodes[ 5] = 13;
	    assign_elem_nodes[ 6] = 1;
	    assign_elem_nodes[ 7] = 8;

	    assign_elem_nodes[ 8] = 11;
	    assign_elem_nodes[ 9] = 19;
	    assign_elem_nodes[10] = 17;
	    assign_elem_nodes[11] = 9;

	    assign_elem_nodes[12] = 3;
	    assign_elem_nodes[13] = 15;
	    assign_elem_nodes[14] = 7;
	    assign_elem_nodes[15] = 18;
	    assign_elem_nodes[16] = 6;
	    assign_elem_nodes[17] = 14;
	    assign_elem_nodes[18] = 2;
	    assign_elem_nodes[19] = 10;


	    break;
	  }
		
	case TET10:
	  {
	    fe_descriptor_id = 118; // Solid Quadratic Tetrahedron
	    assign_elem_nodes[0] = 0;
	    assign_elem_nodes[1] = 4;
	    assign_elem_nodes[2] = 1;
	    assign_elem_nodes[3] = 5;
	    assign_elem_nodes[4] = 2;
	    assign_elem_nodes[5] = 6;
	    assign_elem_nodes[6] = 7;
	    assign_elem_nodes[7] = 8;
	    assign_elem_nodes[8] = 9;
	    assign_elem_nodes[9] = 3;
	    break;
	  }
		
	default:
	  {
	    std::cerr << "ERROR: Element type = " 
		      << elem->type() 
		      << " not supported in "
		      << "UNVIO!"
		      << std::endl;
	    libmesh_error();	
	    break;
	  }
	}


      out_file << std::setw(10) << this->_mesh_data.elem_to_foreign_id(elem)  // element ID
	       << std::setw(10) << fe_descriptor_id                           // type of element
	       << std::setw(10) << phys_prop_tab_dummy                        // not supported 
	       << std::setw(10) << mat_prop_tab_dummy                         // not supported 
	       << std::setw(10) << color_dummy                                // not supported 
	       << std::setw(10) << elem->n_nodes()                            // No. of nodes per element
	       << '\n';

      for (unsigned int j=0; j<elem->n_nodes(); j++)
	{
	  // assign_elem_nodes[j]-th node: i.e., j loops over the
	  // libMesh numbering, and assign_elem_nodes[j] over the
	  // UNV numbering.
	  const Node* node_in_unv_order = elem->get_node(assign_elem_nodes[j]);

	  // new record after 8 id entries
	  if (j==8 || j==16)
	    out_file << '\n';

	  // write foreign label for this node
	  out_file << std::setw(10) << this->_mesh_data.node_to_foreign_id(node_in_unv_order);


	}

      out_file << '\n';

      n_elem_written++;
    }

  if (this->verbose())
    std::cout << "  Finished writing " << n_elem_written << " elements" << std::endl;

  // Write end of dataset
  out_file << "    -1\n";
}