MeshTools::Generation Namespace Reference

Namespaces
namespace	Private
Functions
void	build_cube (UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
void	build_line (UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
void	build_square (UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
void	build_sphere (UnstructuredMesh &mesh, const Real rad=1, const unsigned int nr=2, const ElemType type=INVALID_ELEM)
void	build_delaunay_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 std::vector< TriangleInterface::Hole * > *holes=NULL)

---

Detailed Description

Tools for Mesh generation.

Author:

Benjamin S. Kirk

Date:

2004

Version:

Revision

3391

---

Function Documentation

void MeshTools::Generation::build_cube	(	UnstructuredMesh &	mesh,
		const unsigned int	nx = 0,
		const unsigned int	ny = 0,
		const unsigned int	nz = 0,
		const Real	xmin = 0.,
		const Real	xmax = 1.,
		const Real	ymin = 0.,
		const Real	ymax = 1.,
		const Real	zmin = 0.,
		const Real	zmax = 1.,
		const ElemType	type = INVALID_ELEM,
		const bool	gauss_lobatto_grid = false
	)

Builds a (elements) cube. Defaults to a unit cube (or line in 1D, square in 2D), but the dimensions can be specified through the optional arguments.

Definition at line 55 of file mesh_generation.C.

References MeshBase::add_elem(), MeshBase::add_point(), UnstructuredMesh::all_second_order(), MeshBase::boundary_info, Elem::build_side(), MeshBase::clear(), MeshBase::delete_elem(), libMeshEnums::EDGE2, libMeshEnums::EDGE3, libMeshEnums::EDGE4, MeshBase::elements_begin(), MeshBase::elements_end(), Elem::get_node(), libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, MeshTools::Generation::Private::idx(), libMeshEnums::INVALID_ELEM, BoundaryInfo::invalid_id, MeshBase::mesh_dimension(), MeshBase::n_elem(), MeshBase::n_nodes(), Elem::n_sides(), MeshBase::node(), MeshBase::node_ptr(), libMesh::pi, MeshBase::prepare_for_use(), libMeshEnums::PRISM15, libMeshEnums::PRISM18, libMeshEnums::PRISM6, libMeshEnums::PYRAMID5, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, MeshBase::reserve_elem(), MeshBase::reserve_nodes(), MeshBase::set_mesh_dimension(), Elem::set_node(), side, libMeshEnums::TET10, libMeshEnums::TET4, libMeshEnums::TRI3, and libMeshEnums::TRI6.

Referenced by build_line(), and build_square().

{
  START_LOG("build_cube()", "MeshTools::Generation");

  // Declare that we are using the indexing utility routine
  // in the "Private" part of our current namespace.  If this doesn't
  // work in GCC 2.95.3 we can either remove it or stop supporting
  // 2.95.3 altogether.
  // Changing this to import the whole namespace... just importing idx
  // causes an internal compiler error for Intel Compiler 11.0 on Linux
  // in debug mode.
  using namespace MeshTools::Generation::Private;
  
  // Clear the mesh and start from scratch
  mesh.clear();

  if (nz != 0)
    mesh.set_mesh_dimension(3);
  else if (ny != 0)
    mesh.set_mesh_dimension(2);
  else
    mesh.set_mesh_dimension(1);
  
  switch (mesh.mesh_dimension())
    {

      
      //---------------------------------------------------------------------
      // Build a 1D line
    case 1:      
      {
        libmesh_assert (nx != 0);
        libmesh_assert (ny == 0);
        libmesh_assert (nz == 0);
        libmesh_assert (xmin < xmax);

        // Reserve elements
        switch (type)
          {
          case INVALID_ELEM:
          case EDGE2:
          case EDGE3:
          case EDGE4:
            {
              mesh.reserve_elem (nx);
              break;
            }

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

        // Reserve nodes
        switch (type)
          {
          case INVALID_ELEM:
          case EDGE2:
            {
              mesh.reserve_nodes(nx+1);
              break;
            }

          case EDGE3:
            {
              mesh.reserve_nodes(2*nx+1);
              break;
            }

          case EDGE4:
            {
              mesh.reserve_nodes(3*nx+1);
              break;
            }

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


        // Build the nodes, depends on whether we're using linears, 
        // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
        unsigned int node_id = 0;
        switch(type)
        {
          case INVALID_ELEM:
          case EDGE2:
            {
              for (unsigned int i=0; i<=nx; i++)
              {
                if (gauss_lobatto_grid)
                  mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0),
                        0, 
                        0), node_id++);
                else
                  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx), 
                        0, 
                        0), node_id++);
              }
              break;
            }

          case EDGE3:
            {
              for (unsigned int i=0; i<=2*nx; i++)
              {
                if (gauss_lobatto_grid)
                {
                  // The x location of the point.
                  Real x=0.;

                  // Shortcut variable
                  const Real pi = libMesh::pi;

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

                  // 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
                  {
                    Real cmin = std::cos( pi*(i-1) / static_cast<Real>(2*nx) );
                    Real cmax = std::cos( pi*(i+1) / static_cast<Real>(2*nx) );
                    
                    Real xmin = 0.5*(1.0 - cmin);
                    Real xmax = 0.5*(1.0 - cmax);
                    x = 0.5*(xmin + xmax);
                  }

                  mesh.add_point (Point(x,0.,0.), node_id++);
                }
                else
                  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
                        0,
                        0), node_id++);
              }
              break;
            }

          case EDGE4:
            {
              for (unsigned int i=0; i<=3*nx; i++)
              {
                if (gauss_lobatto_grid)
                {
                  // The x location of the point
                  Real x=0.;
                  
                  const Real pi = libMesh::pi;

                  // Shortcut quantities
                  const Real c = std::cos( pi*i / static_cast<Real>(3*nx) );

                  // If i is multiple of 3, compute a normal Gauss-Lobatto point
                  if (i%3 == 0)
                    x = 0.5*(1.0 - c);

                  // Otherwise, distribute points evenly within the element
                  else
                  {
                    if(i%3 == 1)
                    {
                      Real cmin = std::cos( pi*(i-1) / static_cast<Real>(3*nx) );
                      Real cmax = std::cos( pi*(i+2) / static_cast<Real>(3*nx) );

                      Real xmin = 0.5*(1.0 - cmin);
                      Real xmax = 0.5*(1.0 - cmax);

                      x = (2.*xmin + xmax)/3.;
                    }
                    else
                    if(i%3 == 2)
                    {
                      Real cmin = std::cos( pi*(i-2) / static_cast<Real>(3*nx) );
                      Real cmax = std::cos( pi*(i+1) / static_cast<Real>(3*nx) );

                      Real xmin = 0.5*(1.0 - cmin);
                      Real xmax = 0.5*(1.0 - cmax);

                      x = (xmin + 2.*xmax)/3.;
                    }

                  }

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


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

        // Build the elements of the mesh
        switch(type)
          {
            case INVALID_ELEM:
            case EDGE2:
              {
                for (unsigned int i=0; i<nx; i++)
                {
                  Elem* elem = mesh.add_elem (new Edge2);
                  elem->set_node(0) = mesh.node_ptr(i);
                  elem->set_node(1) = mesh.node_ptr(i+1);

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

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

            case EDGE3:
              {
                for (unsigned int i=0; i<nx; i++)
                {
                  Elem* elem = mesh.add_elem (new Edge3);
                  elem->set_node(0) = mesh.node_ptr(2*i);
                  elem->set_node(2) = mesh.node_ptr(2*i+1);
                  elem->set_node(1) = mesh.node_ptr(2*i+2);

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

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

            case EDGE4:
              {
                for (unsigned int i=0; i<nx; i++)
                {
                  Elem* elem = mesh.add_elem (new Edge4);
                  elem->set_node(0) = mesh.node_ptr(3*i);
                  elem->set_node(2) = mesh.node_ptr(3*i+1);
                  elem->set_node(3) = mesh.node_ptr(3*i+2);
                  elem->set_node(1) = mesh.node_ptr(3*i+3);

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

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

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







      
      
      
      //---------------------------------------------------------------------
      // Build a 2D quadrilateral
    case 2:
      {
        libmesh_assert (nx != 0);
        libmesh_assert (ny != 0);
        libmesh_assert (nz == 0);
        libmesh_assert (xmin < xmax);
        libmesh_assert (ymin < ymax);

        // Reserve elements.  The TRI3 and TRI6 meshes
        // have twice as many elements...
        switch (type)
          {
          case INVALID_ELEM:
          case QUAD4:
          case QUAD8:
          case QUAD9:
            {
              mesh.reserve_elem (nx*ny);
              break;
            }

          case TRI3:
          case TRI6:
            {
              mesh.reserve_elem (2*nx*ny);
              break;
            }

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


        
        // Reserve nodes.  The quadratic element types
        // need to reserve more nodes than the linear types.
        switch (type)
          {
          case INVALID_ELEM:
          case QUAD4:
          case TRI3:
            {
              mesh.reserve_nodes( (nx+1)*(ny+1) );
              break;
            }
            
          case QUAD8:
          case QUAD9:
          case TRI6:
            {
              mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
              break;
            }

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



        // Build the nodes. Depends on whether you are using a linear
        // or quadratic element, and whether you are using a uniform
        // grid or the Gauss-Lobatto grid points.
        unsigned int node_id = 0;
        switch (type)
          {
          case INVALID_ELEM:
          case QUAD4:
          case TRI3:
            {
              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.), node_id++);
                      }
                  
                    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);
        libmesh_assert (xmin < xmax);
        libmesh_assert (ymin < ymax);
        libmesh_assert (zmin < zmax);


        // 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
          case PYRAMID5: // PYRAMID5'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 PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
          case PRISM15:
          case PRISM18:
            {
              // FYI: The resulting TET4 mesh will have exactly
              // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
              // nodes once the additional mid-edge nodes for the HEX27 discretization
              // have been deleted.
              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 PYRAMID5: // PYRAMID5'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
                            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
          case PYRAMID5: // PYRAMID5'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) || (type == PYRAMID5))
                        {
                          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 tets and pyramids: we take the existing
        // HEX27 discretization and split each element into 24
        // sub-tets or 6 sub-pyramids.
        //
        // 24 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) ||
            (type == PYRAMID5)) 
          {
            // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
            std::vector<Elem*> new_elements;
            
            if ((type == TET4) || (type == TET10))
              new_elements.reserve(24*mesh.n_elem());
            else
              new_elements.reserve(6*mesh.n_elem());
              
            // Create tetrahedra or pyramids
            {
              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);

                      if ((type == TET4) || (type == TET10))
                        {
                          // Build 4 sub-tets per side
                          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);
                            }
                        } // end if ((type == TET4) || (type == TET10))

                      else // type==PYRAMID5
                        {
                          // Build 1 sub-pyramid per side.
                          new_elements.push_back(new Pyramid5);
                          Elem* sub_elem = new_elements.back();
                              
                          // Set the base.  Note that since the apex is *inside* the base_hex,
                          // and the pyramid uses a counter-clockwise base numbering, we need to
                          // reverse the [1] and [3] node indices.
                          sub_elem->set_node(0) = side->get_node(0);
                          sub_elem->set_node(1) = side->get_node(3);
                          sub_elem->set_node(2) = side->get_node(2);
                          sub_elem->set_node(3) = side->get_node(1);
                              
                          // Set the apex
                          sub_elem->set_node(4) = apex_node;
                              
                          // If the original hex was a boundary hex, add the new sub_pyr's side
                          // 4 (the square base) with the same b_id.  
                          if (b_id != BoundaryInfo::invalid_id)
                            mesh.boundary_info->add_side(sub_elem, 4, b_id);
                        } // end else type==PYRAMID5
                    }
                }
            }
            

            // 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,TET10,PYRAMID5


        // 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();
      }
    }  

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


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

void MeshTools::Generation::build_delaunay_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 std::vector< TriangleInterface::Hole * > *	holes = NULL
	)

Meshes a rectangular (2D) region (with or without holes) with a Delaunay triangulation. This function internally calls the triangle library written by J.R. Shewchuk.

Definition at line 34 of file mesh_triangle_support.C.

References MeshBase::add_point(), TriangleInterface::attach_hole_list(), bc_id, MeshBase::boundary_info, Elem::build_side(), MeshBase::clear(), TriangleInterface::desired_area(), TriangleInterface::elem_type(), MeshBase::elements_begin(), MeshBase::elements_end(), DofObject::id(), MeshBase::n_nodes(), Elem::n_sides(), Elem::neighbor(), TriangleInterface::PSLG, MeshBase::set_mesh_dimension(), side, TriangleInterface::triangulate(), and TriangleInterface::triangulation_type().

{
  // Check for reasonable size
  libmesh_assert (nx >= 1); // need at least 1 element in x-direction
  libmesh_assert (ny >= 1); // need at least 1 element in y-direction
  libmesh_assert (xmin < xmax);
  libmesh_assert (ymin < ymax);

  // Clear out any data which may have been in the Mesh
  mesh.clear();

  // Make sure the new Mesh will be 2D
  mesh.set_mesh_dimension(2);
  
  // The x and y spacing between boundary points
  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);  

  // Bottom
  for (unsigned int p=0; p<=nx; ++p)
    mesh.add_point(Point(xmin + p*delta_x, ymin));
  
  // Right side
  for (unsigned int p=1; p<ny; ++p)
    mesh.add_point(Point(xmax, ymin + p*delta_y)); 

  // Top
  for (unsigned int p=0; p<=nx; ++p)
    mesh.add_point(Point(xmax - p*delta_x, ymax));

  // Left side
  for (unsigned int p=1; p<ny; ++p)
    mesh.add_point(Point(xmin,  ymax - p*delta_y)); 
                   
  // Be sure we added as many points as we thought we did
  libmesh_assert (mesh.n_nodes() == 2*(nx+ny));

  // Construct the Triangle Interface object
  TriangleInterface t(mesh);

  // Set custom variables for the triangulation
  t.desired_area()       = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
  t.triangulation_type() = TriangleInterface::PSLG;
  t.elem_type()          = type;

  if (holes != NULL)
    t.attach_hole_list(holes);

  // Triangulate!
  t.triangulate();

  // The mesh is now generated, but we still need to mark the boundaries
  // to be consistent with the other build_square routines.  Note that only
  // the external boundaries of the square are given boundary ids, the holes
  // are not numbered.
  MeshBase::element_iterator       el     = mesh.elements_begin();
  const MeshBase::element_iterator end_el = mesh.elements_end();
  for ( ; el != end_el; ++el)
    {
      const Elem* elem = *el;
      
      for (unsigned int s=0; s<elem->n_sides(); s++)
        if (elem->neighbor(s) == NULL)
          {
            AutoPtr<Elem> side (elem->build_side(s));

            // Check the location of the side's midpoint.  Since
            // the square has straight sides, the midpoint is not
            // on the corner and thus it is uniquely on one of the
            // sides.
            Point side_midpoint= 0.5*( (*side->get_node(0)) + (*side->get_node(1)) );

            // The boundary ids are set following the same convention as Quad4 sides
            // bottom = 0
            // right  = 1
            // top = 2
            // left = 3
            // hole = 4
            short int bc_id=4;

            // bottom
            if      (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
              bc_id=0;

            // right
            else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
              bc_id=1;

            // top
            else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
              bc_id=2;

            // left
            else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
              bc_id=3;

            // If the point is not on any of the external boundaries, it
            // is on one of the holes....
            
            // Finally, add this element's information to the boundary info object.
            mesh.boundary_info->add_side(elem->id(), s, bc_id);
          }
    }
  
}

void MeshTools::Generation::build_line	(	UnstructuredMesh &	mesh,
		const unsigned int	nx,
		const Real	xmin = 0.,
		const Real	xmax = 1.,
		const ElemType	type = INVALID_ELEM,
		const bool	gauss_lobatto_grid = false
	)

A specialized build_cube() for 1D meshes

Definition at line 1329 of file mesh_generation.C.

References build_cube().

{
    // This method only makes sense in 1D!
    // But we now just turn a non-1D mesh into a 1D mesh
    //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_sphere	(	UnstructuredMesh &	mesh,
		const Real	rad = 1,
		const unsigned int	nr = 2,
		const ElemType	type = INVALID_ELEM
	)

Meshes a spherical or mapped-spherical domain.

Definition at line 1381 of file mesh_generation.C.

{
        std::cout << "Building a circle/sphere only works with AMR." << std::endl;
        libmesh_error();
}

void MeshTools::Generation::build_square	(	UnstructuredMesh &	mesh,
		const unsigned int	nx,
		const unsigned int	ny,
		const Real	xmin = 0.,
		const Real	xmax = 1.,
		const Real	ymin = 0.,
		const Real	ymax = 1.,
		const ElemType	type = INVALID_ELEM,
		const bool	gauss_lobatto_grid = false
	)

A specialized build_cube() for 2D meshes.

Definition at line 1350 of file mesh_generation.C.

References build_cube().

{
  // This method only makes sense in 2D!
  // But we now just turn a non-2D mesh into a 2D mesh
  //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);
}