mesh_generation.c

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

    }  

  STOP_LOG("build_cube()", "MeshTools::Generation");


  
  // Done building the mesh.  Now prepare it for use.
  mesh.prepare_for_use ();  
}



void MeshTools::Generation::build_line (UnstructuredMesh& mesh,
                                        const unsigned int nx,
                                        const Real xmin, const Real xmax,
                                        const ElemType type,
                                        const bool gauss_lobatto_grid)
{
    // This method only makes sense in 1D!
    libmesh_assert(mesh.mesh_dimension() == 1);

    build_cube(mesh,
               nx, 0, 0,
               xmin, xmax,
               0., 0.,
               0., 0.,
               type,
               gauss_lobatto_grid);
}
                                        


void MeshTools::Generation::build_square (UnstructuredMesh& mesh,
					  const unsigned int nx,
					  const unsigned int ny,
					  const Real xmin, const Real xmax,
					  const Real ymin, const Real ymax,
					  const ElemType type,
					  const bool gauss_lobatto_grid)
{
  // This method only makes sense in 2D!
  libmesh_assert (mesh.mesh_dimension() == 2);

  // Call the build_cube() member to actually do the work for us.
  build_cube (mesh,
	      nx, ny, 0,
	      xmin, xmax,
	      ymin, ymax,
	      0., 0.,
	      type,
	      gauss_lobatto_grid);
}









#ifndef ENABLE_AMR
void MeshTools::Generation::build_sphere (UnstructuredMesh&,
					  const Real,
					  const unsigned int,
					  const ElemType)
{
 	std::cout << "Building a circle/sphere only works with AMR." << std::endl;
 	libmesh_error();
}

#else
	
void MeshTools::Generation::build_sphere (UnstructuredMesh& mesh,
					  const Real rad,
					  const unsigned int nr,
					  const ElemType type)
{
  libmesh_assert (mesh.mesh_dimension() != 1);
  libmesh_assert (rad > 0.);
  //libmesh_assert (nr > 0); // must refine at least once otherwise will end up with a square/cube
  
  START_LOG("build_sphere()", "MeshTools::Generation");

  // Clear the mesh and start from scratch
  mesh.clear();
  
  // Sphere is centered at origin by default
  const Point cent;

  const Sphere sphere (cent, rad);
  
  switch (mesh.mesh_dimension())
    {
      //-----------------------------------------------------------------
      // Build a circle in two dimensions
    case 2:
      {
	const Real sqrt_2     = std::sqrt(2.);
        const Real rad_2      = .25*rad;
        const Real rad_sqrt_2 = rad/sqrt_2;

	// (Temporary) convenient storage for node pointers
	std::vector<Node*> nodes(8);
	    
	// Point 0
	nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.));
	    
	// Point 1
	nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.));
	    
	// Point 2
	nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.));

	// Point 3
	nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.));
	    
	// Point 4
	nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.));
	    
	// Point 5
	nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.));
	    
	// Point 6
	nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.));
	    
	// Point 7
	nodes[7] = mesh.add_point (Point(-rad_sqrt_2, rad_sqrt_2, 0.));

	// Build the elements & set node pointers
	    
	// Element 0
	Elem* elem0 = mesh.add_elem (new Quad4);
	elem0->set_node(0) = nodes[0];
	elem0->set_node(1) = nodes[1];
	elem0->set_node(2) = nodes[2];
	elem0->set_node(3) = nodes[3];
	    
	// Element 1
	Elem* elem1 = mesh.add_elem (new Quad4);
	elem1->set_node(0) = nodes[4];
	elem1->set_node(1) = nodes[0];
	elem1->set_node(2) = nodes[3];
	elem1->set_node(3) = nodes[7];
	    
	// Element 2
	Elem* elem2 = mesh.add_elem (new Quad4);
	elem2->set_node(0) = nodes[4];
	elem2->set_node(1) = nodes[5];
	elem2->set_node(2) = nodes[1];
	elem2->set_node(3) = nodes[0];
	    
	// Element 3
	Elem* elem3 = mesh.add_elem (new Quad4);
	elem3->set_node(0) = nodes[1];
	elem3->set_node(1) = nodes[5];
	elem3->set_node(2) = nodes[6];
	elem3->set_node(3) = nodes[2];
	    
	// Element 4
	Elem* elem4 = mesh.add_elem (new Quad4);
	elem4->set_node(0) = nodes[3];
	elem4->set_node(1) = nodes[2];
	elem4->set_node(2) = nodes[6];
	elem4->set_node(3) = nodes[7];

	break;
      } // end case 2


      

      
      //-----------------------------------------------------------------
      // Build a sphere in three dimensions
    case 3:
      {
	// (Currently) supported types
	if (!((type == HEX8) || (type == HEX27)))
	  {
	    // FIXME: We'd need an all_tet() routine (which could also be used by
	    // build_square()) to do Tets, or Prisms for that matter.
	    std::cerr << "Error: Only HEX8/27 currently supported."
		      << std::endl;
	    libmesh_error();
	  }

	
	// 3D analog of 2D initial grid:
	const Real
	  r_small = 0.25*rad,                      //  0.25 *radius
	  r_med   = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius 
	
	// (Temporary) convenient storage for node pointers
	std::vector<Node*> nodes(16);
	
	// Points 0-7 are the initial HEX8
	nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small));
	nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small));
	nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small));
	nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small));
	nodes[4] = mesh.add_point (Point(-r_small,-r_small,  r_small));
	nodes[5] = mesh.add_point (Point( r_small,-r_small,  r_small));
	nodes[6] = mesh.add_point (Point( r_small, r_small,  r_small));
	nodes[7] = mesh.add_point (Point(-r_small, r_small,  r_small));

	//  Points 8-15 are for the outer hexes, we number them in the same way
	nodes[8]  = mesh.add_point (Point(-r_med,-r_med, -r_med));
	nodes[9]  = mesh.add_point (Point( r_med,-r_med, -r_med));
	nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med));
	nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med));
	nodes[12] = mesh.add_point (Point(-r_med,-r_med,  r_med));
	nodes[13] = mesh.add_point (Point( r_med,-r_med,  r_med));
	nodes[14] = mesh.add_point (Point( r_med, r_med,  r_med));
	nodes[15] = mesh.add_point (Point(-r_med, r_med,  r_med));

	// Now create the elements and add them to the mesh
	// Element 0 - center element
	{
	  Elem* elem0 = mesh.add_elem (new Hex8);
	  elem0->set_node(0) = nodes[0];
	  elem0->set_node(1) = nodes[1];
	  elem0->set_node(2) = nodes[2];
	  elem0->set_node(3) = nodes[3];
	  elem0->set_node(4) = nodes[4];
	  elem0->set_node(5) = nodes[5];
	  elem0->set_node(6) = nodes[6];
	  elem0->set_node(7) = nodes[7];
	}
	
	// Element 1 - "bottom"
	{
	  Elem* elem1 = mesh.add_elem (new Hex8);
	  elem1->set_node(0) = nodes[8];
	  elem1->set_node(1) = nodes[9];
	  elem1->set_node(2) = nodes[10];
	  elem1->set_node(3) = nodes[11];
	  elem1->set_node(4) = nodes[0];
	  elem1->set_node(5) = nodes[1];
	  elem1->set_node(6) = nodes[2];
	  elem1->set_node(7) = nodes[3];
	}
	
	// Element 2 - "front"
	{
	  Elem* elem2 = mesh.add_elem (new Hex8);
	  elem2->set_node(0) = nodes[8];
	  elem2->set_node(1) = nodes[9];
	  elem2->set_node(2) = nodes[1];
	  elem2->set_node(3) = nodes[0];
	  elem2->set_node(4) = nodes[12];
	  elem2->set_node(5) = nodes[13];
	  elem2->set_node(6) = nodes[5];
	  elem2->set_node(7) = nodes[4];
	}

	// Element 3 - "right"
	{
	  Elem* elem3 = mesh.add_elem (new Hex8);
	  elem3->set_node(0) = nodes[1];
	  elem3->set_node(1) = nodes[9];
	  elem3->set_node(2) = nodes[10];
	  elem3->set_node(3) = nodes[2];
	  elem3->set_node(4) = nodes[5];
	  elem3->set_node(5) = nodes[13];
	  elem3->set_node(6) = nodes[14];
	  elem3->set_node(7) = nodes[6];
	}

	// Element 4 - "back"
	{
	  Elem* elem4 = mesh.add_elem (new Hex8);
	  elem4->set_node(0) = nodes[3];
	  elem4->set_node(1) = nodes[2];
	  elem4->set_node(2) = nodes[10];
	  elem4->set_node(3) = nodes[11];
	  elem4->set_node(4) = nodes[7];
	  elem4->set_node(5) = nodes[6];
	  elem4->set_node(6) = nodes[14];
	  elem4->set_node(7) = nodes[15];
	}

	// Element 5 - "left"
	{
	  Elem* elem5 = mesh.add_elem (new Hex8);
	  elem5->set_node(0) = nodes[8];
	  elem5->set_node(1) = nodes[0];
	  elem5->set_node(2) = nodes[3];
	  elem5->set_node(3) = nodes[11];
	  elem5->set_node(4) = nodes[12];
	  elem5->set_node(5) = nodes[4];
	  elem5->set_node(6) = nodes[7];
	  elem5->set_node(7) = nodes[15];
	}

	// Element 6 - "top"
	{
	  Elem* elem6 = mesh.add_elem (new Hex8);
	  elem6->set_node(0) = nodes[4];
	  elem6->set_node(1) = nodes[5];
	  elem6->set_node(2) = nodes[6];
	  elem6->set_node(3) = nodes[7];
	  elem6->set_node(4) = nodes[12];
	  elem6->set_node(5) = nodes[13];
	  elem6->set_node(6) = nodes[14];
	  elem6->set_node(7) = nodes[15];
	}
	
	break;
      } // end case 3
      
    default:
      libmesh_error();
      
    } // end switch (dim)

  
	
  // Now we have the beginnings of a sphere.
  // Add some more elements by doing uniform refinements and
  // popping nodes to the boundary.
  MeshRefinement mesh_refinement (mesh);

  // Loop over the elements, refine, pop nodes to boundary.
  for (unsigned int r=0; r<nr; r++)
    {
      mesh_refinement.uniformly_refine(1);

      MeshBase::element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::element_iterator end = mesh.active_elements_end(); 

      for (; it != end; ++it)
	{
	  Elem* elem = *it;
		
	  for (unsigned int s=0; s<elem->n_sides(); s++)
	    if (elem->neighbor(s) == NULL)
	      {
		AutoPtr<Elem> side(elem->build_side(s));

		// Pop each point to the sphere boundary
		for (unsigned int n=0; n<side->n_nodes(); n++)
		  side->point(n) =
		    sphere.closest_point(side->point(n));
	      }
	}
      }

  // The mesh now contains a refinement hierarchy due to the refinements
  // used to generate the grid.  In order to call other support functions
  // like all_tri() and all_second_order, you need a "flat" mesh file (with no
  // refinement trees) so
  MeshTools::Modification::flatten(mesh);

  // In 2D, convert all the quads to triangles if requested
  if (mesh.mesh_dimension()==2)
    {
      if ((type == TRI6) || (type == TRI3))
	{
	  MeshTools::Modification::all_tri(mesh);
	}
    }

  
  // Convert to second-order elements if the user requested it.
  if (Elem::second_order_equivalent_type(type) == INVALID_ELEM)
    {
      // type is already a second-order, determine if it is the
      // "full-ordered" second-order element, or the "serendipity"
      // second order element.  Note also that all_second_order
      // can't be called once the mesh has been refined.
      bool full_ordered = !((type==QUAD8) || (type==HEX20));
      mesh.all_second_order(full_ordered);

      // And pop to the boundary again...
      MeshBase::element_iterator       it  = mesh.active_elements_begin();
      const MeshBase::element_iterator end = mesh.active_elements_end(); 

      for (; it != end; ++it)
	{
	  Elem* elem = *it;
		
	  for (unsigned int s=0; s<elem->n_sides(); s++)
	    if (elem->neighbor(s) == NULL)
	      {
		AutoPtr<Elem> side(elem->build_side(s));

		// Pop each point to the sphere boundary
		for (unsigned int n=0; n<side->n_nodes(); n++)
		  side->point(n) =
		    sphere.closest_point(side->point(n));
	      }
	}
    }
  

  // The meshes could probably use some smoothing...
  LaplaceMeshSmoother smoother(mesh);
  smoother.smooth(2);
  
  STOP_LOG("build_sphere()", "MeshTools::Generation");

  
  // Done building the mesh.  Now prepare it for use.
  mesh.prepare_for_use();
}

#endif // #ifndef ENABLE_AMR