mesh_generation.c

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

		      }
		  
		    else
		      mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
					    static_cast<Real>(j)/static_cast<Real>(ny),
					    0.), node_id++);
		  }

	      break;
	    }
	    
	  case QUAD8:
	  case QUAD9:
	  case TRI6:
	    {
	      for (unsigned int j=0; j<=(2*ny); j++)
		for (unsigned int i=0; i<=(2*nx); i++)
		  {
		    if (gauss_lobatto_grid)
		      {
			// The x,y locations of the point.
			Real x=0., y=0.;

			// Shortcut variable
			const Real pi = libMesh::pi;
			
			// Shortcut quantities (do not depend on i,j)
			const Real a = std::cos( pi / static_cast<Real>(2*nx) );
			const Real b = std::cos( pi / static_cast<Real>(2*ny) );

			// Shortcut quantities (depend on i,j)
			const Real c = std::cos( pi*i / static_cast<Real>(2*nx) );
			const Real d = std::cos( pi*j / static_cast<Real>(2*ny) );
			
			// If i is even, compute a normal Gauss-Lobatto point
			if (i%2 == 0)
			  x = 0.5*(1.0 - c);

			// Otherwise, it is the average of the previous and next points
			else
			  x = 0.5*(1.0 - a*c);
		      
			// If j is even, compute a normal Gauss-Lobatto point
			if (j%2 == 0)
			  y = 0.5*(1.0 - d);
		      
			// Otherwise, it is the average of the previous and next points
			else
			  y = 0.5*(1.0 - b*d);
		      

			mesh.add_point (Point(x,y,0.), 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),
					    0), node_id++);
		}

	      break;
	    }

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






	// Build the elements.  Each one is a bit different.
	switch (type)
	  {

	  case INVALID_ELEM:
	  case QUAD4:
	    {
	      for (unsigned int j=0; j<ny; j++)
		for (unsigned int i=0; i<nx; i++)
		  {
		    Elem* elem = mesh.add_elem(new Quad4);

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

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


	  case TRI3:
	    {
	      for (unsigned int j=0; j<ny; j++)
		for (unsigned int i=0; i<nx; i++)
		  {
		    Elem* elem = NULL;

		    // Add first Tri3
		    elem = mesh.add_elem(new Tri3);

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

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

		    if (i == (nx-1))
		      mesh.boundary_info->add_side(elem, 1, 1);

		    // Add second Tri3
		    elem = mesh.add_elem(new Tri3);

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

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

		    if (i == 0)
		      mesh.boundary_info->add_side(elem, 2, 3);
		  }
	      break;
	    }



	  case QUAD8:
	  case QUAD9:
	    {
	      for (unsigned int j=0; j<(2*ny); j += 2)
		for (unsigned int i=0; i<(2*nx); i += 2)
		  {
		    Elem* elem = (type == QUAD8) ?
		      mesh.add_elem(new Quad8) :
		      mesh.add_elem(new Quad9);
		  

		    elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
		    elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j)  );
		    elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
		    elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+2)  );
		    elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j)  );
		    elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+2,j+1));
		    elem->set_node(6) = mesh.node_ptr(idx(type,nx,i+1,j+2));
		    elem->set_node(7) = mesh.node_ptr(idx(type,nx,i,j+1)  );
		    if (type == QUAD9)
		      elem->set_node(8) = mesh.node_ptr(idx(type,nx,i+1,j+1));
		  

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


	  case TRI6:
	    {
	      for (unsigned int j=0; j<(2*ny); j += 2)
		for (unsigned int i=0; i<(2*nx); i += 2)
		  {
		    Elem* elem = NULL;

		    // Add first Tri6
		    elem = mesh.add_elem(new Tri6);

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

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

		    if (i == 2*(nx-1))
		      mesh.boundary_info->add_side(elem, 1, 1);

		    // Add second Tri6
		    elem = mesh.add_elem(new Tri6);

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

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

		    if (i == 0)
		      mesh.boundary_info->add_side(elem, 2, 3);

		  }
	      break;
	    };

	    
	  default:
	    {
	      std::cerr << "ERROR: Unrecognized 2D 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;
	  }

	
	break;
      }








      


      //---------------------------------------------------------------------
      // Build a 3D mesh using hexahedral or prismatic elements.
    case 3:
      {
	libmesh_assert (nx != 0);
	libmesh_assert (ny != 0);
	libmesh_assert (nz != 0);


	// Reserve elements.  Meshes with prismatic elements require
	// twice as many elements.
	switch (type)
	  {
	  case INVALID_ELEM:
	  case HEX8:
	  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
	    {
	      mesh.reserve_elem(nx*ny*nz);
	      break;
	    }

	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	    {
	      mesh.reserve_elem(2*nx*ny*nz);
	      break;
	    }

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




	
	// Reserve nodes.  Quadratic elements need twice as many nodes as linear elements.
	switch (type)
	  {
	  case INVALID_ELEM:
	  case HEX8:
	  case PRISM6:
	    {
	      mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
	      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
	  case PRISM15:
	  case PRISM18:
	    {
	      mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
	      break;
	    }

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




	// Build the nodes.
        unsigned int node_id = 0;
	switch (type)
	  {
	  case INVALID_ELEM:
	  case HEX8:
	  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++)
		    {
		      if (gauss_lobatto_grid)
			{
			  // Shortcut variable
			  const Real pi = libMesh::pi;

			  mesh.add_point (Point(0.5*(1.0 - std::cos(pi*static_cast<Real>(i)/static_cast<Real>(nx))),
						0.5*(1.0 - std::cos(pi*static_cast<Real>(j)/static_cast<Real>(ny))),
						0.5*(1.0 - std::cos(pi*static_cast<Real>(k)/static_cast<Real>(nz)))), node_id++);
			}
		      
		      else
			mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(nx),
					     static_cast<Real>(j)/static_cast<Real>(ny),
					     static_cast<Real>(k)/static_cast<Real>(nz)), node_id++);
		    }
	      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
	  case PRISM15:
	  case PRISM18:
	    {
	      for (unsigned int k=0; k<=(2*nz); k++)
		for (unsigned int j=0; j<=(2*ny); j++)
		  for (unsigned int i=0; i<=(2*nx); i++)
		    {
		      if (gauss_lobatto_grid)
			{
			  // Shortcut variable
			  const Real pi = libMesh::pi;
			  
			  // The x,y locations of the point.
			  Real x=0., y=0., z=0.;

			  // Shortcut quantities (do not depend on i,j)
			  const Real a = std::cos( pi / static_cast<Real>(2*nx) );
			  const Real b = std::cos( pi / static_cast<Real>(2*ny) );

			  // Shortcut quantities (depend on i,j)
			  const Real c = std::cos( pi*i / static_cast<Real>(2*nx) );
			  const Real d = std::cos( pi*j / static_cast<Real>(2*ny) );
			  
			  // Additional shortcut quantities (for 3D)
			  const Real e = std::cos( pi / static_cast<Real>(2*nz) );
			  const Real f = std::cos( pi*k / static_cast<Real>(2*nz) );
			  
			  // If i is even, compute a normal Gauss-Lobatto point
			  if (i%2 == 0)
			    x = 0.5*(1.0 - c);

			  // Otherwise, it is the average of the previous and next points
			  else
			    x = 0.5*(1.0 - a*c);
		      
			  // If j is even, compute a normal Gauss-Lobatto point
			  if (j%2 == 0)
			    y = 0.5*(1.0 - d);
		      
			  // Otherwise, it is the average of the previous and next points
			  else