mesh_generation.c

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

			    y = 0.5*(1.0 - b*d);

			  // If k is even, compute a normal Gauss-Lobatto point
			  if (k%2 == 0)
			    z = 0.5*(1.0 - f);
		      
			  // Otherwise, it is the average of the previous and next points
			  else
			    z = 0.5*(1.0 - e*f);
		      
			  
			  mesh.add_point (Point(x,y,z), node_id++);
			}

		      else
			mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
					     static_cast<Real>(j)/static_cast<Real>(2*ny),
					     static_cast<Real>(k)/static_cast<Real>(2*nz)), node_id++);
		    }
	      break;
	    }


	  default:
	    {
	      std::cerr << "ERROR: Unrecognized 3D element type." << std::endl;
	      libmesh_error();
	    }
	  }

	    

	
	// Build the elements.
	switch (type)
	  {
	  case INVALID_ELEM:
	  case HEX8:
	    {
	      for (unsigned int k=0; k<nz; k++)
		for (unsigned int j=0; j<ny; j++)
		  for (unsigned int i=0; i<nx; i++)
		    {
		      Elem* elem = mesh.add_elem(new Hex8);
			
		      elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k)      );
		      elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    );
		      elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
		      elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    );
		      elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1)    );
		      elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  );
		      elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
		      elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  );
			
		      if (k == 0)
			mesh.boundary_info->add_side(elem, 0, 0);
			
		      if (k == (nz-1))
			mesh.boundary_info->add_side(elem, 5, 5);
			
		      if (j == 0)
			mesh.boundary_info->add_side(elem, 1, 1);

		      if (j == (ny-1))
			mesh.boundary_info->add_side(elem, 3, 3);
 			
		      if (i == 0)
			mesh.boundary_info->add_side(elem, 4, 4);
 			
		      if (i == (nx-1))
			mesh.boundary_info->add_side(elem, 2, 2);
		    }
	      break;
	    }




	  case PRISM6:
	    {
	      for (unsigned int k=0; k<nz; k++)
		for (unsigned int j=0; j<ny; j++)
		  for (unsigned int i=0; i<nx; i++)
		    {
		      // First Prism
		      Elem* elem = NULL;
		      elem = mesh.add_elem(new Prism6);

		      elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k)      );
		      elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    );
		      elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    );
		      elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1)    );
		      elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  );
		      elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  );

		      // Add sides for first prism to boundary info object
		      if (i==0)
			mesh.boundary_info->add_side(elem, 3, 4);

		      if (j==0)
			mesh.boundary_info->add_side(elem, 1, 1);

		      if (k==0)
			mesh.boundary_info->add_side(elem, 0, 0);

		      if (k == (nz-1))
			mesh.boundary_info->add_side(elem, 4, 5);

		      // Second Prism
		      elem = mesh.add_elem(new Prism6);

		      elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    );
		      elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
		      elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    );
		      elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  );
		      elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
		      elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  );

		      // Add sides for second prism to boundary info object
		      if (i == (nx-1))
			mesh.boundary_info->add_side(elem, 1, 2);

		      if (j == (ny-1))
			mesh.boundary_info->add_side(elem, 2, 3);

		      if (k==0)
			mesh.boundary_info->add_side(elem, 0, 0);

		      if (k == (nz-1))
			mesh.boundary_info->add_side(elem, 4, 5);
		    }
	      break;
	    }




	    
	    
	  case HEX20:
	  case HEX27:
	  case TET4: // TET4's are created from an initial HEX27 discretization
	  case TET10: // TET10's are created from an initial HEX27 discretization
	    {
	      for (unsigned int k=0; k<(2*nz); k += 2)
		for (unsigned int j=0; j<(2*ny); j += 2)
		  for (unsigned int i=0; i<(2*nx); i += 2)
		    {
		      Elem* elem = (type == HEX20) ?
			mesh.add_elem(new Hex20) :
			mesh.add_elem(new Hex27);
		    
		      elem->set_node(0)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k)  );
		      elem->set_node(1)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k)  );
		      elem->set_node(2)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k)  );
		      elem->set_node(3)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k)  );
		      elem->set_node(4)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+2));
		      elem->set_node(5)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+2));
		      elem->set_node(6)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2));
		      elem->set_node(7)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+2));
		      elem->set_node(8)  = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k)  );
		      elem->set_node(9)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k)  );
		      elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k)  );
		      elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k)  );
		      elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+1));
		      elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+1));
		      elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
		      elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+1));
		      elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+2));
		      elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
		      elem->set_node(18) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
		      elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+2));
		      if ((type == HEX27) || (type == TET4) || (type == TET10))
			{
			  elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
			  elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+1));
			  elem->set_node(22) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
			  elem->set_node(23) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
			  elem->set_node(24) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+1));
			  elem->set_node(25) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
			  elem->set_node(26) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
			}
			
			
		      if (k == 0)
			mesh.boundary_info->add_side(elem, 0, 0);
			
		      if (k == 2*(nz-1))
			mesh.boundary_info->add_side(elem, 5, 5);
			
		      if (j == 0)
			mesh.boundary_info->add_side(elem, 1, 1);

		      if (j == 2*(ny-1))
			mesh.boundary_info->add_side(elem, 3, 3);
 			
		      if (i == 0)
			mesh.boundary_info->add_side(elem, 4, 4);
 			
		      if (i == 2*(nx-1))
			mesh.boundary_info->add_side(elem, 2, 2);
		    }
	      break;
	    }



	    
	  case PRISM15:
	  case PRISM18:
	    {
	      for (unsigned int k=0; k<(2*nz); k += 2)
		for (unsigned int j=0; j<(2*ny); j += 2)
		  for (unsigned int i=0; i<(2*nx); i += 2)
		    {
		      // First Prism
		      Elem* elem = NULL;
		      elem = ((type == PRISM15) ?
			      mesh.add_elem(new Prism15) :
			      mesh.add_elem(new Prism18));

		      elem->set_node(0)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k)  );
		      elem->set_node(1)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k)  );
		      elem->set_node(2)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k)  );
		      elem->set_node(3)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+2));
		      elem->set_node(4)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+2));
		      elem->set_node(5)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+2));
		      elem->set_node(6)  = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k)  );
		      elem->set_node(7)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
		      elem->set_node(8)  = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k)  );
		      elem->set_node(9)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+1));
		      elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+1));
		      elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+1));
		      elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+2));
		      elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
		      elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+2));
		      if (type == PRISM18)
			{
			  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+1));
			  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
			  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+1));
			}

		      // Add sides for first prism to boundary info object
		      if (i==0)
			mesh.boundary_info->add_side(elem, 3, 4);

		      if (j==0)
			mesh.boundary_info->add_side(elem, 1, 1);

		      if (k==0)
			mesh.boundary_info->add_side(elem, 0, 0);

		      if (k == 2*(nz-1))
			mesh.boundary_info->add_side(elem, 4, 5);

		      
		      // Second Prism
		      elem = ((type == PRISM15) ?
			      mesh.add_elem(new Prism15) :
			      mesh.add_elem(new Prism18));
			
		      elem->set_node(0)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,k)     );
		      elem->set_node(1)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k)   );
		      elem->set_node(2)  = mesh.node_ptr(idx(type,nx,ny,i,j+2,k)     );
		      elem->set_node(3)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2)   );
		      elem->set_node(4)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) );
		      elem->set_node(5)  = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2)   );
		      elem->set_node(6)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k)  );
		      elem->set_node(7)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k)  );
		      elem->set_node(8)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
		      elem->set_node(9)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1)  );
		      elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
		      elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1)  );
		      elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
		      elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
		      elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
		      if (type == PRISM18)
			{
			  elem->set_node(15)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
			  elem->set_node(16)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
			  elem->set_node(17)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
			}
		      
		      // Add sides for second prism to boundary info object
		      if (i == 2*(nx-1))
			mesh.boundary_info->add_side(elem, 1, 2);

		      if (j == 2*(ny-1))
			mesh.boundary_info->add_side(elem, 2, 3);

		      if (k==0)
			mesh.boundary_info->add_side(elem, 0, 0);

		      if (k == 2*(nz-1))
			mesh.boundary_info->add_side(elem, 4, 5);

		    }
	      break;
	    }




	    
	  default:
	    {
	      std::cerr << "ERROR: Unrecognized 3D element type." << std::endl;
	      libmesh_error();
	    }
	  }

	    


	//.......................................
	// Scale the nodal positions
	for (unsigned int p=0; p<mesh.n_nodes(); p++)
	  {
	    mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
	    mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
	    mesh.node(p)(2) = (mesh.node(p)(2))*(zmax-zmin) + zmin;
	  }




	// Additional work for TET4 and TET10: we take the existing HEX27 discretization
	// and split each element into 24 sub tets.  This isn't the minimum-possible
	// number of tets, but it obviates any concerns about the edge orientations
	// between the various elements.
	if ((type == TET4) || (type == TET10)) 
	  {
	    // Temporary storage for new elements. (24 per hex)
	    std::vector<Elem*> new_elements;
	    new_elements.reserve(24*mesh.n_elem());

     
	    // Create tetrahedra
	    {
	      MeshBase::element_iterator       el     = mesh.elements_begin();
	      const MeshBase::element_iterator end_el = mesh.elements_end();

	      for ( ; el != end_el;  ++el)
		{
		  // Get a pointer to the HEX27 element.
		  Elem* base_hex = *el;

		  // Get a pointer to the node located at the HEX27 centroid
		  Node* apex_node = base_hex->get_node(26);
	 
		  for (unsigned int s=0; s<base_hex->n_sides(); ++s)
		    {
		      // Get the boundary ID for this side
		      short int b_id = mesh.boundary_info->boundary_id(*el, s);
		      
		      // Need to build the full-ordered side!
		      AutoPtr<Elem> side = base_hex->build_side(s);
	     
		      for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
			{
			  new_elements.push_back( new Tet4 );
			  Elem* sub_elem = new_elements.back();
			  sub_elem->set_node(0) = side->get_node(sub_tet);
			  sub_elem->set_node(1) = side->get_node(8);                           // centroid of the face
			  sub_elem->set_node(2) = side->get_node(sub_tet==3 ? 0 : sub_tet+1 ); // wrap-around
			  sub_elem->set_node(3) = apex_node;                                   // apex node always used!

			  // If the original hex was a boundary hex, add the new sub_tet's side
			  // 0 with the same b_id.  Note: the tets are all aligned so that their
			  // side 0 is on the boundary.
			  if (b_id != BoundaryInfo::invalid_id)
			    mesh.boundary_info->add_side(sub_elem, 0, b_id);
			}
		    }
		}
	    }
	    

	    // Delete the original HEX27 elements from the mesh, and the boundary info structure.
	    {
	      MeshBase::element_iterator       el     = mesh.elements_begin();
	      const MeshBase::element_iterator end_el = mesh.elements_end();
       
	      for ( ; el != end_el;  ++el)
		{
		  mesh.boundary_info->remove(*el); // Safe even if *el has no boundary info.
		  mesh.delete_elem(*el);
		}
	    }

	    // Add the new elements
	    for (unsigned int i=0; i<new_elements.size(); ++i)
	      mesh.add_elem(new_elements[i]);
	    
	  } // end if (type == TET4) || (type == TET10)


	// Use all_second_order to convert the TET4's to TET10's
	if (type == TET10)
	  {
	    mesh.all_second_order();
	  }

	
	break;
      } // end case dim==3

    default:
      {
	libmesh_error();
      }