libMesh::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_point (UnstructuredMesh &mesh, 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, const unsigned int n_smooth=2, const bool flat=true)
void build_extrusion (UnstructuredMesh &mesh, const MeshBase &cross_section, const unsigned int nz, RealVectorValue extrusion_vector)
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

Function Documentation

void libMesh::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 $ nx \times ny \times nz $ (elements) cube. Defaults to a unit cube (or line in 1D, square in 2D), but the dimensions can be specified through the optional arguments.

Boundary ids are set to be equal to the side indexing on a master hex

Definition at line 149 of file mesh_generation.C.

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

Referenced by build_line(), build_point(), and build_square().

00158 {
00159   START_LOG("build_cube()", "MeshTools::Generation");
00160 
00161   // Declare that we are using the indexing utility routine
00162   // in the "Private" part of our current namespace.  If this doesn't
00163   // work in GCC 2.95.3 we can either remove it or stop supporting
00164   // 2.95.3 altogether.
00165   // Changing this to import the whole namespace... just importing idx
00166   // causes an internal compiler error for Intel Compiler 11.0 on Linux
00167   // in debug mode.
00168   using namespace MeshTools::Generation::Private;
00169 
00170   // Clear the mesh and start from scratch
00171   mesh.clear();
00172 
00173   if (nz != 0)
00174     mesh.set_mesh_dimension(3);
00175   else if (ny != 0)
00176     mesh.set_mesh_dimension(2);
00177   else if (nx != 0)
00178     mesh.set_mesh_dimension(1);
00179   else
00180     mesh.set_mesh_dimension(0);
00181 
00182   switch (mesh.mesh_dimension())
00183     {
00184       //---------------------------------------------------------------------
00185       // Build a 0D point
00186     case 0:
00187       {
00188         libmesh_assert_equal_to (nx, 0);
00189         libmesh_assert_equal_to (ny, 0);
00190         libmesh_assert_equal_to (nz, 0);
00191 
00192         libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
00193 
00194         // Build one nodal element for the mesh
00195         mesh.add_point (Point(0, 0, 0), 0);
00196         Elem* elem = mesh.add_elem (new NodeElem);
00197         elem->set_node(0) = mesh.node_ptr(0);
00198 
00199         break;
00200       }
00201 
00202 
00203 
00204       //---------------------------------------------------------------------
00205       // Build a 1D line
00206     case 1:
00207       {
00208         libmesh_assert_not_equal_to (nx, 0);
00209         libmesh_assert_equal_to (ny, 0);
00210         libmesh_assert_equal_to (nz, 0);
00211         libmesh_assert_less (xmin, xmax);
00212 
00213         // Reserve elements
00214         switch (type)
00215           {
00216           case INVALID_ELEM:
00217           case EDGE2:
00218           case EDGE3:
00219           case EDGE4:
00220             {
00221               mesh.reserve_elem (nx);
00222               break;
00223             }
00224 
00225           default:
00226             {
00227               libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
00228               libmesh_error();
00229             }
00230           }
00231 
00232         // Reserve nodes
00233         switch (type)
00234           {
00235           case INVALID_ELEM:
00236           case EDGE2:
00237             {
00238               mesh.reserve_nodes(nx+1);
00239               break;
00240             }
00241 
00242           case EDGE3:
00243             {
00244               mesh.reserve_nodes(2*nx+1);
00245               break;
00246             }
00247 
00248           case EDGE4:
00249             {
00250               mesh.reserve_nodes(3*nx+1);
00251               break;
00252             }
00253 
00254           default:
00255             {
00256               libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
00257               libmesh_error();
00258             }
00259           }
00260 
00261 
00262         // Build the nodes, depends on whether we're using linears,
00263         // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
00264         unsigned int node_id = 0;
00265         switch(type)
00266         {
00267           case INVALID_ELEM:
00268           case EDGE2:
00269             {
00270               for (unsigned int i=0; i<=nx; i++)
00271               {
00272                 if (gauss_lobatto_grid)
00273                   mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0),
00274                         0,
00275                         0), node_id++);
00276                 else
00277                   mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
00278                         0,
00279                         0), node_id++);
00280               }
00281               break;
00282             }
00283 
00284           case EDGE3:
00285             {
00286               for (unsigned int i=0; i<=2*nx; i++)
00287               {
00288                 if (gauss_lobatto_grid)
00289                 {
00290                   // The x location of the point.
00291                   Real x=0.;
00292 
00293                   // Shortcut quantities (do not depend on i)
00294                   const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
00295 
00296                   // If i is even, compute a normal Gauss-Lobatto point
00297                   if (i%2 == 0)
00298                     x = 0.5*(1.0 - c);
00299 
00300                   // Otherwise, it is the average of the previous and next points
00301                   else
00302                   {
00303                     Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(2*nx) );
00304                     Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(2*nx) );
00305 
00306                     Real gl_xmin = 0.5*(1.0 - cmin);
00307                     Real gl_xmax = 0.5*(1.0 - cmax);
00308                     x = 0.5*(gl_xmin + gl_xmax);
00309                   }
00310 
00311                   mesh.add_point (Point(x,0.,0.), node_id++);
00312                 }
00313                 else
00314                   mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
00315                         0,
00316                         0), node_id++);
00317               }
00318               break;
00319             }
00320 
00321           case EDGE4:
00322             {
00323               for (unsigned int i=0; i<=3*nx; i++)
00324               {
00325                 if (gauss_lobatto_grid)
00326                 {
00327                   // The x location of the point
00328                   Real x=0.;
00329 
00330                   // Shortcut quantities
00331                   const Real c = std::cos( libMesh::pi*i / static_cast<Real>(3*nx) );
00332 
00333                   // If i is multiple of 3, compute a normal Gauss-Lobatto point
00334                   if (i%3 == 0)
00335                     x = 0.5*(1.0 - c);
00336 
00337                   // Otherwise, distribute points evenly within the element
00338                   else
00339                   {
00340                     if(i%3 == 1)
00341                     {
00342                       Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(3*nx) );
00343                       Real cmax = std::cos( libMesh::pi*(i+2) / static_cast<Real>(3*nx) );
00344 
00345                       Real gl_xmin = 0.5*(1.0 - cmin);
00346                       Real gl_xmax = 0.5*(1.0 - cmax);
00347 
00348                       x = (2.*gl_xmin + gl_xmax)/3.;
00349                     }
00350                     else
00351                     if(i%3 == 2)
00352                     {
00353                       Real cmin = std::cos( libMesh::pi*(i-2) / static_cast<Real>(3*nx) );
00354                       Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(3*nx) );
00355 
00356                       Real gl_xmin = 0.5*(1.0 - cmin);
00357                       Real gl_xmax = 0.5*(1.0 - cmax);
00358 
00359                       x = (gl_xmin + 2.*gl_xmax)/3.;
00360                     }
00361 
00362                   }
00363 
00364                   mesh.add_point (Point(x,0.,0.), node_id++);
00365                 }
00366                 else
00367                 mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(3*nx),
00368                         0,
00369                         0), node_id++);
00370               }
00371 
00372 
00373 
00374               break;
00375             }
00376 
00377           default:
00378             {
00379               libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
00380               libmesh_error();
00381             }
00382 
00383         }
00384 
00385         // Build the elements of the mesh
00386         switch(type)
00387           {
00388             case INVALID_ELEM:
00389             case EDGE2:
00390               {
00391                 for (unsigned int i=0; i<nx; i++)
00392                 {
00393                   Elem* elem = mesh.add_elem (new Edge2);
00394                   elem->set_node(0) = mesh.node_ptr(i);
00395                   elem->set_node(1) = mesh.node_ptr(i+1);
00396 
00397                   if (i == 0)
00398                     mesh.boundary_info->add_side(elem, 0, 0);
00399 
00400                   if (i == (nx-1))
00401                     mesh.boundary_info->add_side(elem, 1, 1);
00402                 }
00403               break;
00404               }
00405 
00406             case EDGE3:
00407               {
00408                 for (unsigned int i=0; i<nx; i++)
00409                 {
00410                   Elem* elem = mesh.add_elem (new Edge3);
00411                   elem->set_node(0) = mesh.node_ptr(2*i);
00412                   elem->set_node(2) = mesh.node_ptr(2*i+1);
00413                   elem->set_node(1) = mesh.node_ptr(2*i+2);
00414 
00415                   if (i == 0)
00416                     mesh.boundary_info->add_side(elem, 0, 0);
00417 
00418                   if (i == (nx-1))
00419                     mesh.boundary_info->add_side(elem, 1, 1);
00420                 }
00421               break;
00422               }
00423 
00424             case EDGE4:
00425               {
00426                 for (unsigned int i=0; i<nx; i++)
00427                 {
00428                   Elem* elem = mesh.add_elem (new Edge4);
00429                   elem->set_node(0) = mesh.node_ptr(3*i);
00430                   elem->set_node(2) = mesh.node_ptr(3*i+1);
00431                   elem->set_node(3) = mesh.node_ptr(3*i+2);
00432                   elem->set_node(1) = mesh.node_ptr(3*i+3);
00433 
00434                   if (i == 0)
00435                     mesh.boundary_info->add_side(elem, 0, 0);
00436 
00437                   if (i == (nx-1))
00438                     mesh.boundary_info->add_side(elem, 1, 1);
00439                 }
00440               break;
00441               }
00442 
00443             default:
00444               {
00445                 libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
00446                 libmesh_error();
00447               }
00448           }
00449 
00450         // Scale the nodal positions
00451         for (unsigned int p=0; p<mesh.n_nodes(); p++)
00452           mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
00453 
00454         // Add sideset names to boundary info
00455         mesh.boundary_info->sideset_name(0) = "left";
00456         mesh.boundary_info->sideset_name(1) = "right";
00457 
00458         // Add nodeset names to boundary info
00459         mesh.boundary_info->nodeset_name(0) = "left";
00460         mesh.boundary_info->nodeset_name(1) = "right";
00461 
00462         break;
00463       }
00464 
00465 
00466 
00467 
00468 
00469 
00470 
00471 
00472 
00473 
00474       //---------------------------------------------------------------------
00475       // Build a 2D quadrilateral
00476     case 2:
00477       {
00478         libmesh_assert_not_equal_to (nx, 0);
00479         libmesh_assert_not_equal_to (ny, 0);
00480         libmesh_assert_equal_to (nz, 0);
00481         libmesh_assert_less (xmin, xmax);
00482         libmesh_assert_less (ymin, ymax);
00483 
00484         // Reserve elements.  The TRI3 and TRI6 meshes
00485         // have twice as many elements...
00486         switch (type)
00487           {
00488           case INVALID_ELEM:
00489           case QUAD4:
00490           case QUAD8:
00491           case QUAD9:
00492             {
00493               mesh.reserve_elem (nx*ny);
00494               break;
00495             }
00496 
00497           case TRI3:
00498           case TRI6:
00499             {
00500               mesh.reserve_elem (2*nx*ny);
00501               break;
00502             }
00503 
00504           default:
00505             {
00506               libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
00507               libmesh_error();
00508             }
00509           }
00510 
00511 
00512 
00513         // Reserve nodes.  The quadratic element types
00514         // need to reserve more nodes than the linear types.
00515         switch (type)
00516           {
00517           case INVALID_ELEM:
00518           case QUAD4:
00519           case TRI3:
00520             {
00521               mesh.reserve_nodes( (nx+1)*(ny+1) );
00522               break;
00523             }
00524 
00525           case QUAD8:
00526           case QUAD9:
00527           case TRI6:
00528             {
00529               mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
00530               break;
00531             }
00532 
00533 
00534           default:
00535             {
00536               libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
00537               libmesh_error();
00538             }
00539           }
00540 
00541 
00542 
00543         // Build the nodes. Depends on whether you are using a linear
00544         // or quadratic element, and whether you are using a uniform
00545         // grid or the Gauss-Lobatto grid points.
00546         unsigned int node_id = 0;
00547         switch (type)
00548           {
00549           case INVALID_ELEM:
00550           case QUAD4:
00551           case TRI3:
00552             {
00553               for (unsigned int j=0; j<=ny; j++)
00554                 for (unsigned int i=0; i<=nx; i++)
00555                   {
00556                     if (gauss_lobatto_grid)
00557                       {
00558                         mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
00559                                               0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
00560                                               0.), node_id++);
00561                       }
00562 
00563                     else
00564                       mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
00565                                             static_cast<Real>(j)/static_cast<Real>(ny),
00566                                             0.), node_id++);
00567                   }
00568 
00569               break;
00570             }
00571 
00572           case QUAD8:
00573           case QUAD9:
00574           case TRI6:
00575             {
00576               for (unsigned int j=0; j<=(2*ny); j++)
00577                 for (unsigned int i=0; i<=(2*nx); i++)
00578                   {
00579                     if (gauss_lobatto_grid)
00580                       {
00581                         // The x,y locations of the point.
00582                         Real x=0., y=0.;
00583 
00584                         // Shortcut quantities (do not depend on i,j)
00585                         const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
00586                         const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
00587 
00588                         // Shortcut quantities (depend on i,j)
00589                         const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
00590                         const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
00591 
00592                         // If i is even, compute a normal Gauss-Lobatto point
00593                         if (i%2 == 0)
00594                           x = 0.5*(1.0 - c);
00595 
00596                         // Otherwise, it is the average of the previous and next points
00597                         else
00598                           x = 0.5*(1.0 - a*c);
00599 
00600                         // If j is even, compute a normal Gauss-Lobatto point
00601                         if (j%2 == 0)
00602                           y = 0.5*(1.0 - d);
00603 
00604                         // Otherwise, it is the average of the previous and next points
00605                         else
00606                           y = 0.5*(1.0 - b*d);
00607 
00608 
00609                         mesh.add_point (Point(x,y,0.), node_id++);
00610                       }
00611 
00612 
00613                     else
00614                       mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
00615                                             static_cast<Real>(j)/static_cast<Real>(2*ny),
00616                                             0), node_id++);
00617                 }
00618 
00619               break;
00620             }
00621 
00622 
00623           default:
00624             {
00625               libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
00626               libmesh_error();
00627             }
00628           }
00629 
00630 
00631 
00632 
00633 
00634 
00635         // Build the elements.  Each one is a bit different.
00636         switch (type)
00637           {
00638 
00639           case INVALID_ELEM:
00640           case QUAD4:
00641             {
00642               for (unsigned int j=0; j<ny; j++)
00643                 for (unsigned int i=0; i<nx; i++)
00644                   {
00645                     Elem* elem = mesh.add_elem(new Quad4);
00646 
00647                     elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
00648                     elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j)  );
00649                     elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
00650                     elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+1)  );
00651 
00652                     if (j == 0)
00653                       mesh.boundary_info->add_side(elem, 0, 0);
00654 
00655                     if (j == (ny-1))
00656                       mesh.boundary_info->add_side(elem, 2, 2);
00657 
00658                     if (i == 0)
00659                       mesh.boundary_info->add_side(elem, 3, 3);
00660 
00661                     if (i == (nx-1))
00662                       mesh.boundary_info->add_side(elem, 1, 1);
00663                   }
00664               break;
00665             }
00666 
00667 
00668           case TRI3:
00669             {
00670               for (unsigned int j=0; j<ny; j++)
00671                 for (unsigned int i=0; i<nx; i++)
00672                   {
00673                     Elem* elem = NULL;
00674 
00675                     // Add first Tri3
00676                     elem = mesh.add_elem(new Tri3);
00677 
00678                     elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
00679                     elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j)  );
00680                     elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
00681 
00682                     if (j == 0)
00683                       mesh.boundary_info->add_side(elem, 0, 0);
00684 
00685                     if (i == (nx-1))
00686                       mesh.boundary_info->add_side(elem, 1, 1);
00687 
00688                     // Add second Tri3
00689                     elem = mesh.add_elem(new Tri3);
00690 
00691                     elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
00692                     elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j+1));
00693                     elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+1)  );
00694 
00695                     if (j == (ny-1))
00696                       mesh.boundary_info->add_side(elem, 1, 2);
00697 
00698                     if (i == 0)
00699                       mesh.boundary_info->add_side(elem, 2, 3);
00700                   }
00701               break;
00702             }
00703 
00704 
00705 
00706           case QUAD8:
00707           case QUAD9:
00708             {
00709               for (unsigned int j=0; j<(2*ny); j += 2)
00710                 for (unsigned int i=0; i<(2*nx); i += 2)
00711                   {
00712                     Elem* elem = (type == QUAD8) ?
00713                       mesh.add_elem(new Quad8) :
00714                       mesh.add_elem(new Quad9);
00715 
00716 
00717                     elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
00718                     elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j)  );
00719                     elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
00720                     elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+2)  );
00721                     elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j)  );
00722                     elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+2,j+1));
00723                     elem->set_node(6) = mesh.node_ptr(idx(type,nx,i+1,j+2));
00724                     elem->set_node(7) = mesh.node_ptr(idx(type,nx,i,j+1)  );
00725                     if (type == QUAD9)
00726                       elem->set_node(8) = mesh.node_ptr(idx(type,nx,i+1,j+1));
00727 
00728 
00729                     if (j == 0)
00730                       mesh.boundary_info->add_side(elem, 0, 0);
00731 
00732                     if (j == 2*(ny-1))
00733                       mesh.boundary_info->add_side(elem, 2, 2);
00734 
00735                     if (i == 0)
00736                       mesh.boundary_info->add_side(elem, 3, 3);
00737 
00738                     if (i == 2*(nx-1))
00739                       mesh.boundary_info->add_side(elem, 1, 1);
00740                   }
00741               break;
00742             }
00743 
00744 
00745           case TRI6:
00746             {
00747               for (unsigned int j=0; j<(2*ny); j += 2)
00748                 for (unsigned int i=0; i<(2*nx); i += 2)
00749                   {
00750                     Elem* elem = NULL;
00751 
00752                     // Add first Tri6
00753                     elem = mesh.add_elem(new Tri6);
00754 
00755                     elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
00756                     elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j)  );
00757                     elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
00758                     elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j)  );
00759                     elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+2,j+1));
00760                     elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+1,j+1));
00761 
00762                     if (j == 0)
00763                       mesh.boundary_info->add_side(elem, 0, 0);
00764 
00765                     if (i == 2*(nx-1))
00766                       mesh.boundary_info->add_side(elem, 1, 1);
00767 
00768                     // Add second Tri6
00769                     elem = mesh.add_elem(new Tri6);
00770 
00771                     elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j)    );
00772                     elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j+2));
00773                     elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+2)  );
00774                     elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j+1));
00775                     elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j+2));
00776                     elem->set_node(5) = mesh.node_ptr(idx(type,nx,i,j+1)  );
00777 
00778                     if (j == 2*(ny-1))
00779                       mesh.boundary_info->add_side(elem, 1, 2);
00780 
00781                     if (i == 0)
00782                       mesh.boundary_info->add_side(elem, 2, 3);
00783 
00784                   }
00785               break;
00786             };
00787 
00788 
00789           default:
00790             {
00791               libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
00792               libmesh_error();
00793             }
00794           }
00795 
00796 
00797 
00798 
00799         // Scale the nodal positions
00800         for (unsigned int p=0; p<mesh.n_nodes(); p++)
00801           {
00802             mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
00803             mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
00804           }
00805 
00806         // Add sideset names to boundary info
00807         mesh.boundary_info->sideset_name(0) = "bottom";
00808         mesh.boundary_info->sideset_name(1) = "right";
00809         mesh.boundary_info->sideset_name(2) = "top";
00810         mesh.boundary_info->sideset_name(3) = "left";
00811 
00812         // Add nodeset names to boundary info
00813         mesh.boundary_info->nodeset_name(0) = "bottom";
00814         mesh.boundary_info->nodeset_name(1) = "right";
00815         mesh.boundary_info->nodeset_name(2) = "top";
00816         mesh.boundary_info->nodeset_name(3) = "left";
00817 
00818         break;
00819       }
00820 
00821 
00822 
00823 
00824 
00825 
00826 
00827 
00828 
00829 
00830 
00831       //---------------------------------------------------------------------
00832       // Build a 3D mesh using hexahedral or prismatic elements.
00833     case 3:
00834       {
00835         libmesh_assert_not_equal_to (nx, 0);
00836         libmesh_assert_not_equal_to (ny, 0);
00837         libmesh_assert_not_equal_to (nz, 0);
00838         libmesh_assert_less (xmin, xmax);
00839         libmesh_assert_less (ymin, ymax);
00840         libmesh_assert_less (zmin, zmax);
00841 
00842 
00843         // Reserve elements.  Meshes with prismatic elements require
00844         // twice as many elements.
00845         switch (type)
00846           {
00847           case INVALID_ELEM:
00848           case HEX8:
00849           case HEX20:
00850           case HEX27:
00851           case TET4:  // TET4's are created from an initial HEX27 discretization
00852           case TET10: // TET10's are created from an initial HEX27 discretization
00853           case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
00854             {
00855               mesh.reserve_elem(nx*ny*nz);
00856               break;
00857             }
00858 
00859           case PRISM6:
00860           case PRISM15:
00861           case PRISM18:
00862             {
00863               mesh.reserve_elem(2*nx*ny*nz);
00864               break;
00865             }
00866 
00867           default:
00868             {
00869               libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
00870               libmesh_error();
00871             }
00872           }
00873 
00874 
00875 
00876 
00877 
00878         // Reserve nodes.  Quadratic elements need twice as many nodes as linear elements.
00879         switch (type)
00880           {
00881           case INVALID_ELEM:
00882           case HEX8:
00883           case PRISM6:
00884             {
00885               mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
00886               break;
00887             }
00888 
00889           case HEX20:
00890           case HEX27:
00891           case TET4: // TET4's are created from an initial HEX27 discretization
00892           case TET10: // TET10's are created from an initial HEX27 discretization
00893           case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
00894           case PRISM15:
00895           case PRISM18:
00896             {
00897               // FYI: The resulting TET4 mesh will have exactly
00898               // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
00899               // nodes once the additional mid-edge nodes for the HEX27 discretization
00900               // have been deleted.
00901               mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
00902               break;
00903             }
00904 
00905           default:
00906             {
00907               libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
00908               libmesh_error();
00909             }
00910           }
00911 
00912 
00913 
00914 
00915         // Build the nodes.
00916         unsigned int node_id = 0;
00917         switch (type)
00918           {
00919           case INVALID_ELEM:
00920           case HEX8:
00921           case PRISM6:
00922             {
00923               for (unsigned int k=0; k<=nz; k++)
00924                 for (unsigned int j=0; j<=ny; j++)
00925                   for (unsigned int i=0; i<=nx; i++)
00926                     {
00927                       if (gauss_lobatto_grid)
00928                         {
00929                           mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
00930                                                 0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
00931                                                 0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(k)/static_cast<Real>(nz)))), node_id++);
00932                         }
00933 
00934                       else
00935                         mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(nx),
00936                                              static_cast<Real>(j)/static_cast<Real>(ny),
00937                                              static_cast<Real>(k)/static_cast<Real>(nz)), node_id++);
00938                     }
00939               break;
00940             }
00941 
00942           case HEX20:
00943           case HEX27:
00944           case TET4: // TET4's are created from an initial HEX27 discretization
00945           case TET10: // TET10's are created from an initial HEX27 discretization
00946           case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
00947           case PRISM15:
00948           case PRISM18:
00949             {
00950               for (unsigned int k=0; k<=(2*nz); k++)
00951                 for (unsigned int j=0; j<=(2*ny); j++)
00952                   for (unsigned int i=0; i<=(2*nx); i++)
00953                     {
00954                       if (gauss_lobatto_grid)
00955                         {
00956                           // The x,y locations of the point.
00957                           Real x=0., y=0., z=0.;
00958 
00959                           // Shortcut quantities (do not depend on i,j)
00960                           const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
00961                           const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
00962 
00963                           // Shortcut quantities (depend on i,j)
00964                           const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
00965                           const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
00966 
00967                           // Additional shortcut quantities (for 3D)
00968                           const Real e = std::cos( libMesh::pi / static_cast<Real>(2*nz) );
00969                           const Real f = std::cos( libMesh::pi*k / static_cast<Real>(2*nz) );
00970 
00971                           // If i is even, compute a normal Gauss-Lobatto point
00972                           if (i%2 == 0)
00973                             x = 0.5*(1.0 - c);
00974 
00975                           // Otherwise, it is the average of the previous and next points
00976                           else
00977                             x = 0.5*(1.0 - a*c);
00978 
00979                           // If j is even, compute a normal Gauss-Lobatto point
00980                           if (j%2 == 0)
00981                             y = 0.5*(1.0 - d);
00982 
00983                           // Otherwise, it is the average of the previous and next points
00984                           else
00985                             y = 0.5*(1.0 - b*d);
00986 
00987                           // If k is even, compute a normal Gauss-Lobatto point
00988                           if (k%2 == 0)
00989                             z = 0.5*(1.0 - f);
00990 
00991                           // Otherwise, it is the average of the previous and next points
00992                           else
00993                             z = 0.5*(1.0 - e*f);
00994 
00995 
00996                           mesh.add_point (Point(x,y,z), node_id++);
00997                         }
00998 
00999                       else
01000                         mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
01001                                              static_cast<Real>(j)/static_cast<Real>(2*ny),
01002                                              static_cast<Real>(k)/static_cast<Real>(2*nz)), node_id++);
01003                     }
01004               break;
01005             }
01006 
01007 
01008           default:
01009             {
01010               libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
01011               libmesh_error();
01012             }
01013           }
01014 
01015 
01016 
01017 
01018         // Build the elements.
01019         switch (type)
01020           {
01021           case INVALID_ELEM:
01022           case HEX8:
01023             {
01024               for (unsigned int k=0; k<nz; k++)
01025                 for (unsigned int j=0; j<ny; j++)
01026                   for (unsigned int i=0; i<nx; i++)
01027                     {
01028                       Elem* elem = mesh.add_elem(new Hex8);
01029 
01030                       elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k)      );
01031                       elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    );
01032                       elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
01033                       elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    );
01034                       elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1)    );
01035                       elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  );
01036                       elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
01037                       elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  );
01038 
01039                       if (k == 0)
01040                         mesh.boundary_info->add_side(elem, 0, 0);
01041 
01042                       if (k == (nz-1))
01043                         mesh.boundary_info->add_side(elem, 5, 5);
01044 
01045                       if (j == 0)
01046                         mesh.boundary_info->add_side(elem, 1, 1);
01047 
01048                       if (j == (ny-1))
01049                         mesh.boundary_info->add_side(elem, 3, 3);
01050 
01051                       if (i == 0)
01052                         mesh.boundary_info->add_side(elem, 4, 4);
01053 
01054                       if (i == (nx-1))
01055                         mesh.boundary_info->add_side(elem, 2, 2);
01056                     }
01057               break;
01058             }
01059 
01060 
01061 
01062 
01063           case PRISM6:
01064             {
01065               for (unsigned int k=0; k<nz; k++)
01066                 for (unsigned int j=0; j<ny; j++)
01067                   for (unsigned int i=0; i<nx; i++)
01068                     {
01069                       // First Prism
01070                       Elem* elem = NULL;
01071                       elem = mesh.add_elem(new Prism6);
01072 
01073                       elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k)      );
01074                       elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    );
01075                       elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    );
01076                       elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1)    );
01077                       elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  );
01078                       elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  );
01079 
01080                       // Add sides for first prism to boundary info object
01081                       if (i==0)
01082                         mesh.boundary_info->add_side(elem, 3, 4);
01083 
01084                       if (j==0)
01085                         mesh.boundary_info->add_side(elem, 1, 1);
01086 
01087                       if (k==0)
01088                         mesh.boundary_info->add_side(elem, 0, 0);
01089 
01090                       if (k == (nz-1))
01091                         mesh.boundary_info->add_side(elem, 4, 5);
01092 
01093                       // Second Prism
01094                       elem = mesh.add_elem(new Prism6);
01095 
01096                       elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k)    );
01097                       elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
01098                       elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k)    );
01099                       elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1)  );
01100                       elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
01101                       elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1)  );
01102 
01103                       // Add sides for second prism to boundary info object
01104                       if (i == (nx-1))
01105                         mesh.boundary_info->add_side(elem, 1, 2);
01106 
01107                       if (j == (ny-1))
01108                         mesh.boundary_info->add_side(elem, 2, 3);
01109 
01110                       if (k==0)
01111                         mesh.boundary_info->add_side(elem, 0, 0);
01112 
01113                       if (k == (nz-1))
01114                         mesh.boundary_info->add_side(elem, 4, 5);
01115                     }
01116               break;
01117             }
01118 
01119 
01120 
01121 
01122 
01123 
01124           case HEX20:
01125           case HEX27:
01126           case TET4: // TET4's are created from an initial HEX27 discretization
01127           case TET10: // TET10's are created from an initial HEX27 discretization
01128           case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
01129             {
01130               for (unsigned int k=0; k<(2*nz); k += 2)
01131                 for (unsigned int j=0; j<(2*ny); j += 2)
01132                   for (unsigned int i=0; i<(2*nx); i += 2)
01133                     {
01134                       Elem* elem = (type == HEX20) ?
01135                         mesh.add_elem(new Hex20) :
01136                         mesh.add_elem(new Hex27);
01137 
01138                       elem->set_node(0)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k)  );
01139                       elem->set_node(1)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k)  );
01140                       elem->set_node(2)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k)  );
01141                       elem->set_node(3)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k)  );
01142                       elem->set_node(4)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+2));
01143                       elem->set_node(5)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+2));
01144                       elem->set_node(6)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2));
01145                       elem->set_node(7)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+2));
01146                       elem->set_node(8)  = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k)  );
01147                       elem->set_node(9)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k)  );
01148                       elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k)  );
01149                       elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k)  );
01150                       elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+1));
01151                       elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+1));
01152                       elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
01153                       elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+1));
01154                       elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+2));
01155                       elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
01156                       elem->set_node(18) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
01157                       elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+2));
01158                       if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == PYRAMID5))
01159                         {
01160                           elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
01161                           elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+1));
01162                           elem->set_node(22) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
01163                           elem->set_node(23) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
01164                           elem->set_node(24) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+1));
01165                           elem->set_node(25) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
01166                           elem->set_node(26) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
01167                         }
01168 
01169 
01170                       if (k == 0)
01171                         mesh.boundary_info->add_side(elem, 0, 0);
01172 
01173                       if (k == 2*(nz-1))
01174                         mesh.boundary_info->add_side(elem, 5, 5);
01175 
01176                       if (j == 0)
01177                         mesh.boundary_info->add_side(elem, 1, 1);
01178 
01179                       if (j == 2*(ny-1))
01180                         mesh.boundary_info->add_side(elem, 3, 3);
01181 
01182                       if (i == 0)
01183                         mesh.boundary_info->add_side(elem, 4, 4);
01184 
01185                       if (i == 2*(nx-1))
01186                         mesh.boundary_info->add_side(elem, 2, 2);
01187                     }
01188               break;
01189             }
01190 
01191 
01192 
01193 
01194           case PRISM15:
01195           case PRISM18:
01196             {
01197               for (unsigned int k=0; k<(2*nz); k += 2)
01198                 for (unsigned int j=0; j<(2*ny); j += 2)
01199                   for (unsigned int i=0; i<(2*nx); i += 2)
01200                     {
01201                       // First Prism
01202                       Elem* elem = NULL;
01203                       elem = ((type == PRISM15) ?
01204                               mesh.add_elem(new Prism15) :
01205                               mesh.add_elem(new Prism18));
01206 
01207                       elem->set_node(0)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k)  );
01208                       elem->set_node(1)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k)  );
01209                       elem->set_node(2)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k)  );
01210                       elem->set_node(3)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+2));
01211                       elem->set_node(4)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+2));
01212                       elem->set_node(5)  = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+2));
01213                       elem->set_node(6)  = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k)  );
01214                       elem->set_node(7)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
01215                       elem->set_node(8)  = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k)  );
01216                       elem->set_node(9)  = mesh.node_ptr(idx(type,nx,ny,i,  j,  k+1));
01217                       elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j,  k+1));
01218                       elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,  j+2,k+1));
01219                       elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+2));
01220                       elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
01221                       elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+2));
01222                       if (type == PRISM18)
01223                         {
01224                           elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+1));
01225                           elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
01226                           elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+1));
01227                         }
01228 
01229                       // Add sides for first prism to boundary info object
01230                       if (i==0)
01231                         mesh.boundary_info->add_side(elem, 3, 4);
01232 
01233                       if (j==0)
01234                         mesh.boundary_info->add_side(elem, 1, 1);
01235 
01236                       if (k==0)
01237                         mesh.boundary_info->add_side(elem, 0, 0);
01238 
01239                       if (k == 2*(nz-1))
01240                         mesh.boundary_info->add_side(elem, 4, 5);
01241 
01242 
01243                       // Second Prism
01244                       elem = ((type == PRISM15) ?
01245                               mesh.add_elem(new Prism15) :
01246                               mesh.add_elem(new Prism18));
01247 
01248                       elem->set_node(0)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,k)     );
01249                       elem->set_node(1)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k)   );
01250                       elem->set_node(2)  = mesh.node_ptr(idx(type,nx,ny,i,j+2,k)     );
01251                       elem->set_node(3)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2)   );
01252                       elem->set_node(4)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) );
01253                       elem->set_node(5)  = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2)   );
01254                       elem->set_node(6)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k)  );
01255                       elem->set_node(7)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k)  );
01256                       elem->set_node(8)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
01257                       elem->set_node(9)  = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1)  );
01258                       elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
01259                       elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1)  );
01260                       elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
01261                       elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
01262                       elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
01263                       if (type == PRISM18)
01264                         {
01265                           elem->set_node(15)  = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
01266                           elem->set_node(16)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
01267                           elem->set_node(17)  = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
01268                         }
01269 
01270                       // Add sides for second prism to boundary info object
01271                       if (i == 2*(nx-1))
01272                         mesh.boundary_info->add_side(elem, 1, 2);
01273 
01274                       if (j == 2*(ny-1))
01275                         mesh.boundary_info->add_side(elem, 2, 3);
01276 
01277                       if (k==0)
01278                         mesh.boundary_info->add_side(elem, 0, 0);
01279 
01280                       if (k == 2*(nz-1))
01281                         mesh.boundary_info->add_side(elem, 4, 5);
01282 
01283                     }
01284               break;
01285             }
01286 
01287 
01288 
01289 
01290 
01291           default:
01292             {
01293               libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
01294               libmesh_error();
01295             }
01296           }
01297 
01298 
01299 
01300 
01301         //.......................................
01302         // Scale the nodal positions
01303         for (unsigned int p=0; p<mesh.n_nodes(); p++)
01304           {
01305             mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
01306             mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
01307             mesh.node(p)(2) = (mesh.node(p)(2))*(zmax-zmin) + zmin;
01308           }
01309 
01310 
01311 
01312 
01313         // Additional work for tets and pyramids: we take the existing
01314         // HEX27 discretization and split each element into 24
01315         // sub-tets or 6 sub-pyramids.
01316         //
01317         // 24 isn't the minimum-possible number of tets, but it
01318         // obviates any concerns about the edge orientations between
01319         // the various elements.
01320         if ((type == TET4) ||
01321             (type == TET10) ||
01322             (type == PYRAMID5))
01323           {
01324             // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
01325             std::vector<Elem*> new_elements;
01326 
01327             if ((type == TET4) || (type == TET10))
01328               new_elements.reserve(24*mesh.n_elem());
01329             else
01330               new_elements.reserve(6*mesh.n_elem());
01331 
01332             // Create tetrahedra or pyramids
01333             {
01334               MeshBase::element_iterator       el     = mesh.elements_begin();
01335               const MeshBase::element_iterator end_el = mesh.elements_end();
01336 
01337               for ( ; el != end_el;  ++el)
01338                 {
01339                   // Get a pointer to the HEX27 element.
01340                   Elem* base_hex = *el;
01341 
01342                   // Get a pointer to the node located at the HEX27 centroid
01343                   Node* apex_node = base_hex->get_node(26);
01344 
01345                   for (unsigned int s=0; s<base_hex->n_sides(); ++s)
01346                     {
01347                       // Get the boundary ID for this side
01348                       boundary_id_type b_id = mesh.boundary_info->boundary_id(*el, s);
01349 
01350                       // Need to build the full-ordered side!
01351                       AutoPtr<Elem> side = base_hex->build_side(s);
01352 
01353                       if ((type == TET4) || (type == TET10))
01354                         {
01355                           // Build 4 sub-tets per side
01356                           for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
01357                             {
01358                               new_elements.push_back( new Tet4 );
01359                               Elem* sub_elem = new_elements.back();
01360                               sub_elem->set_node(0) = side->get_node(sub_tet);
01361                               sub_elem->set_node(1) = side->get_node(8);                           // centroid of the face
01362                               sub_elem->set_node(2) = side->get_node(sub_tet==3 ? 0 : sub_tet+1 ); // wrap-around
01363                               sub_elem->set_node(3) = apex_node;                                   // apex node always used!
01364 
01365                               // If the original hex was a boundary hex, add the new sub_tet's side
01366                               // 0 with the same b_id.  Note: the tets are all aligned so that their
01367                               // side 0 is on the boundary.
01368                               if (b_id != BoundaryInfo::invalid_id)
01369                                 mesh.boundary_info->add_side(sub_elem, 0, b_id);
01370                             }
01371                         } // end if ((type == TET4) || (type == TET10))
01372 
01373                       else // type==PYRAMID5
01374                         {
01375                           // Build 1 sub-pyramid per side.
01376                           new_elements.push_back(new Pyramid5);
01377                           Elem* sub_elem = new_elements.back();
01378 
01379                           // Set the base.  Note that since the apex is *inside* the base_hex,
01380                           // and the pyramid uses a counter-clockwise base numbering, we need to
01381                           // reverse the [1] and [3] node indices.
01382                           sub_elem->set_node(0) = side->get_node(0);
01383                           sub_elem->set_node(1) = side->get_node(3);
01384                           sub_elem->set_node(2) = side->get_node(2);
01385                           sub_elem->set_node(3) = side->get_node(1);
01386 
01387                           // Set the apex
01388                           sub_elem->set_node(4) = apex_node;
01389 
01390                           // If the original hex was a boundary hex, add the new sub_pyr's side
01391                           // 4 (the square base) with the same b_id.
01392                           if (b_id != BoundaryInfo::invalid_id)
01393                             mesh.boundary_info->add_side(sub_elem, 4, b_id);
01394                         } // end else type==PYRAMID5
01395                     }
01396                 }
01397             }
01398 
01399 
01400             // Delete the original HEX27 elements from the mesh, and the boundary info structure.
01401             {
01402               MeshBase::element_iterator       el     = mesh.elements_begin();
01403               const MeshBase::element_iterator end_el = mesh.elements_end();
01404 
01405               for ( ; el != end_el;  ++el)
01406                 {
01407                   mesh.boundary_info->remove(*el); // Safe even if *el has no boundary info.
01408                   mesh.delete_elem(*el);
01409                 }
01410             }
01411 
01412             // Add the new elements
01413             for (unsigned int i=0; i<new_elements.size(); ++i)
01414               mesh.add_elem(new_elements[i]);
01415 
01416           } // end if (type == TET4,TET10,PYRAMID5
01417 
01418 
01419         // Use all_second_order to convert the TET4's to TET10's
01420         if (type == TET10)
01421           {
01422             mesh.all_second_order();
01423           }
01424 
01425         // Add sideset names to boundary info (Z axis out of the screen)
01426         mesh.boundary_info->sideset_name(0) = "back";
01427         mesh.boundary_info->sideset_name(1) = "bottom";
01428         mesh.boundary_info->sideset_name(2) = "right";
01429         mesh.boundary_info->sideset_name(3) = "top";
01430         mesh.boundary_info->sideset_name(4) = "left";
01431         mesh.boundary_info->sideset_name(5) = "front";
01432 
01433         // Add nodeset names to boundary info
01434         mesh.boundary_info->nodeset_name(0) = "back";
01435         mesh.boundary_info->nodeset_name(1) = "bottom";
01436         mesh.boundary_info->nodeset_name(2) = "right";
01437         mesh.boundary_info->nodeset_name(3) = "top";
01438         mesh.boundary_info->nodeset_name(4) = "left";
01439         mesh.boundary_info->nodeset_name(5) = "front";
01440 
01441         break;
01442       } // end case dim==3
01443 
01444     default:
01445       {
01446         libmesh_error();
01447       }
01448     }
01449 
01450   STOP_LOG("build_cube()", "MeshTools::Generation");
01451 
01452 
01453 
01454   // Done building the mesh.  Now prepare it for use.
01455   mesh.prepare_for_use (/*skip_renumber =*/ false);
01456 }

void libMesh::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 2206 of file mesh_generation.C.

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

02213 {
02214   // Check for reasonable size
02215   libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
02216   libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
02217   libmesh_assert_less (xmin, xmax);
02218   libmesh_assert_less (ymin, ymax);
02219 
02220   // Clear out any data which may have been in the Mesh
02221   mesh.clear();
02222 
02223   // Make sure the new Mesh will be 2D
02224   mesh.set_mesh_dimension(2);
02225 
02226   // The x and y spacing between boundary points
02227   const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
02228   const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
02229 
02230   // Bottom
02231   for (unsigned int p=0; p<=nx; ++p)
02232     mesh.add_point(Point(xmin + p*delta_x, ymin));
02233 
02234   // Right side
02235   for (unsigned int p=1; p<ny; ++p)
02236     mesh.add_point(Point(xmax, ymin + p*delta_y));
02237 
02238   // Top
02239   for (unsigned int p=0; p<=nx; ++p)
02240     mesh.add_point(Point(xmax - p*delta_x, ymax));
02241 
02242   // Left side
02243   for (unsigned int p=1; p<ny; ++p)
02244     mesh.add_point(Point(xmin,  ymax - p*delta_y));
02245 
02246   // Be sure we added as many points as we thought we did
02247   libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
02248 
02249   // Construct the Triangle Interface object
02250   TriangleInterface t(mesh);
02251 
02252   // Set custom variables for the triangulation
02253   t.desired_area()       = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
02254   t.triangulation_type() = TriangleInterface::PSLG;
02255   t.elem_type()          = type;
02256 
02257   if (holes != NULL)
02258     t.attach_hole_list(holes);
02259 
02260   // Triangulate!
02261   t.triangulate();
02262 
02263   // The mesh is now generated, but we still need to mark the boundaries
02264   // to be consistent with the other build_square routines.  Note that all
02265   // hole boundary elements get the same ID, 4.
02266   MeshBase::element_iterator       el     = mesh.elements_begin();
02267   const MeshBase::element_iterator end_el = mesh.elements_end();
02268   for ( ; el != end_el; ++el)
02269     {
02270       const Elem* elem = *el;
02271 
02272       for (unsigned int s=0; s<elem->n_sides(); s++)
02273         if (elem->neighbor(s) == NULL)
02274           {
02275             AutoPtr<Elem> side (elem->build_side(s));
02276 
02277             // Check the location of the side's midpoint.  Since
02278             // the square has straight sides, the midpoint is not
02279             // on the corner and thus it is uniquely on one of the
02280             // sides.
02281             Point side_midpoint= 0.5f*( (*side->get_node(0)) + (*side->get_node(1)) );
02282 
02283             // The boundary ids are set following the same convention as Quad4 sides
02284             // bottom = 0
02285             // right  = 1
02286             // top = 2
02287             // left = 3
02288             // hole = 4
02289             boundary_id_type bc_id=4;
02290 
02291             // bottom
02292             if      (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
02293               bc_id=0;
02294 
02295             // right
02296             else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
02297               bc_id=1;
02298 
02299             // top
02300             else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
02301               bc_id=2;
02302 
02303             // left
02304             else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
02305               bc_id=3;
02306 
02307             // If the point is not on any of the external boundaries, it
02308             // is on one of the holes....
02309 
02310             // Finally, add this element's information to the boundary info object.
02311             mesh.boundary_info->add_side(elem->id(), s, bc_id);
02312           }
02313     }
02314 
02315 } // end build_delaunay_square

void libMesh::MeshTools::Generation::build_extrusion ( UnstructuredMesh &  mesh,
const MeshBase &  cross_section,
const unsigned int  nz,
RealVectorValue  extrusion_vector 
)

Meshes the tensor product of a 1D and a 1D-or-2D domain.

Definition at line 1986 of file mesh_generation.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::MeshBase::boundary_info, libMesh::MeshBase::delete_remote_elements(), libMeshEnums::EDGE2, libMeshEnums::EDGE3, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::MeshBase::is_serial(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_sides(), libMesh::MeshBase::node_ptr(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::Elem::parent(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD9, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMeshEnums::SECOND, libMesh::DofObject::set_id(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), libMeshEnums::TRI3, libMeshEnums::TRI6, and libMesh::Elem::type().

01990 {
01991   if (!cross_section.n_elem())
01992     return;
01993 
01994   START_LOG("build_extrusion()", "MeshTools::Generation");
01995 
01996   dof_id_type orig_elem = cross_section.n_elem();
01997   dof_id_type orig_nodes = cross_section.n_nodes();
01998 
01999   unsigned int order = 1;
02000 
02001   // If cross_section is distributed, so is its extrusion
02002   if (!cross_section.is_serial())
02003     mesh.delete_remote_elements();
02004 
02005   // We know a priori how many elements we'll need
02006   mesh.reserve_elem(nz*orig_elem);
02007 
02008   // For straightforward meshes we need one or two additional layers per
02009   // element.
02010   if ((*cross_section.elements_begin())->default_order() == SECOND)
02011     order = 2;
02012 
02013   mesh.reserve_nodes((order*nz+1)*orig_nodes);
02014 
02015   MeshBase::const_node_iterator       nd  = cross_section.nodes_begin();
02016   const MeshBase::const_node_iterator nend = cross_section.nodes_end();
02017   for (; nd!=nend; ++nd)
02018     {
02019       const Node* node = *nd;
02020 
02021       for (unsigned int k=0; k != order*nz+1; ++k)
02022         {
02023           Node *new_node =
02024             mesh.add_point(*node +
02025                            (extrusion_vector * k / nz / order),
02026                            node->id() + (k * orig_nodes),
02027                            node->processor_id());
02028 
02029           const std::vector<boundary_id_type> ids_to_copy =
02030             cross_section.boundary_info->boundary_ids(node);
02031 
02032           mesh.boundary_info->add_node(new_node, ids_to_copy);
02033         }
02034     }
02035 
02036   const std::set<boundary_id_type> &side_ids =
02037     cross_section.boundary_info->get_side_boundary_ids();
02038   const boundary_id_type next_side_id = side_ids.empty() ?
02039     0 : *side_ids.rbegin() + 1;
02040 
02041   MeshBase::const_element_iterator       el  = cross_section.elements_begin();
02042   const MeshBase::const_element_iterator end = cross_section.elements_end();
02043   for (; el!=end; ++el)
02044     {
02045       const Elem* elem = *el;
02046       const ElemType etype = elem->type();
02047 
02048       // build_extrusion currently only works on coarse meshes
02049       libmesh_assert (!elem->parent());
02050 
02051       // We need a map from low-D to high-D sides for boundary id
02052       // setting
02053       std::vector<unsigned char> sidemap(4);
02054 
02055       for (unsigned int k=0; k != nz; ++k)
02056         {
02057           Elem *new_elem;
02058           switch (etype)
02059             {
02060             case EDGE2:
02061               {
02062                 new_elem = new Quad4;
02063                 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
02064                 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
02065                 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
02066                 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
02067                 break;
02068               }
02069             case EDGE3:
02070               {
02071                 new_elem = new Quad9;
02072                 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
02073                 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
02074                 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
02075                 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
02076                 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
02077                 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
02078                 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
02079                 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
02080                 new_elem->set_node(8) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
02081                 break;
02082               }
02083             case TRI3:
02084               {
02085                 new_elem = new Prism6;
02086                 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
02087                 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
02088                 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
02089                 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
02090                 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
02091                 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
02092                 break;
02093               }
02094             case TRI6:
02095               {
02096                 new_elem = new Prism18;
02097                 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
02098                 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
02099                 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
02100                 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
02101                 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
02102                 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
02103                 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
02104                 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
02105                 new_elem->set_node(8) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
02106                 new_elem->set_node(9) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
02107                 new_elem->set_node(10) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
02108                 new_elem->set_node(11) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
02109                 new_elem->set_node(12) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
02110                 new_elem->set_node(13) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
02111                 new_elem->set_node(14) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
02112                 new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
02113                 new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
02114                 new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
02115                 break;
02116               }
02117             case QUAD4:
02118               {
02119                 new_elem = new Hex8;
02120                 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
02121                 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
02122                 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
02123                 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (k * orig_nodes));
02124                 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
02125                 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
02126                 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
02127                 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((k+1) * orig_nodes));
02128                 break;
02129               }
02130             case QUAD9:
02131               {
02132                 new_elem = new Hex27;
02133                 new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
02134                 new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
02135                 new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
02136                 new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
02137                 new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
02138                 new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
02139                 new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
02140                 new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
02141                 new_elem->set_node(8) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
02142                 new_elem->set_node(9) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
02143                 new_elem->set_node(10) = mesh.node_ptr(elem->get_node(6)->id() + (2*k * orig_nodes));
02144                 new_elem->set_node(11) = mesh.node_ptr(elem->get_node(7)->id() + (2*k * orig_nodes));
02145                 new_elem->set_node(12) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
02146                 new_elem->set_node(13) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
02147                 new_elem->set_node(14) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
02148                 new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
02149                 new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
02150                 new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
02151                 new_elem->set_node(18) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+2) * orig_nodes));
02152                 new_elem->set_node(19) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+2) * orig_nodes));
02153                 new_elem->set_node(20) = mesh.node_ptr(elem->get_node(8)->id() + (2*k * orig_nodes));
02154                 new_elem->set_node(21) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
02155                 new_elem->set_node(22) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
02156                 new_elem->set_node(23) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+1) * orig_nodes));
02157                 new_elem->set_node(24) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+1) * orig_nodes));
02158                 new_elem->set_node(25) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+2) * orig_nodes));
02159                 new_elem->set_node(26) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+1) * orig_nodes));
02160                 break;
02161               }
02162             default:
02163               {
02164                 libmesh_not_implemented();
02165                 break;
02166               }
02167             }
02168 
02169           new_elem->set_id(elem->id() + (k * orig_elem));
02170           new_elem->processor_id() = elem->processor_id();
02171 
02172           // maintain the subdomain_id
02173           new_elem->subdomain_id() = elem->subdomain_id();
02174 
02175           new_elem = mesh.add_elem(new_elem);
02176 
02177           // Copy any old boundary ids on all sides
02178           for (unsigned int s = 0; s != elem->n_sides(); ++s)
02179             {
02180               const std::vector<boundary_id_type> ids_to_copy =
02181                 cross_section.boundary_info->boundary_ids(elem, s);
02182 
02183               mesh.boundary_info->add_side(new_elem, s+1, ids_to_copy);
02184             }
02185 
02186           // Give new boundary ids to bottom and top
02187           if (k == 0)
02188             mesh.boundary_info->add_side(new_elem, 0, next_side_id);
02189           if (k == nz-1)
02190             mesh.boundary_info->add_side(new_elem, elem->n_sides()+1, next_side_id+1);
02191         }
02192     }
02193 
02194   STOP_LOG("build_extrusion()", "MeshTools::Generation");
02195 
02196   // Done building the mesh.  Now prepare it for use.
02197   mesh.prepare_for_use(/*skip_renumber =*/ false);
02198 }

void libMesh::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

Boundary ids are set to be equal to the side indexing on a master edge

Definition at line 1478 of file mesh_generation.C.

References build_cube().

01483 {
01484     // This method only makes sense in 1D!
01485     // But we now just turn a non-1D mesh into a 1D mesh
01486     //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
01487 
01488     build_cube(mesh,
01489                nx, 0, 0,
01490                xmin, xmax,
01491                0., 0.,
01492                0., 0.,
01493                type,
01494                gauss_lobatto_grid);
01495 }

void libMesh::MeshTools::Generation::build_point ( UnstructuredMesh &  mesh,
const ElemType  type = INVALID_ELEM,
const bool  gauss_lobatto_grid = false 
)

A specialized build_cube() for 0D meshes. The resulting mesh is a single NodeElem suitable for ODE tests

Definition at line 1460 of file mesh_generation.C.

References build_cube().

01463 {
01464     // This method only makes sense in 0D!
01465     // But we now just turn a non-0D mesh into a 0D mesh
01466     //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
01467 
01468     build_cube(mesh,
01469                0, 0, 0,
01470                0., 0.,
01471                0., 0.,
01472                0., 0.,
01473                type,
01474                gauss_lobatto_grid);
01475 }

void libMesh::MeshTools::Generation::build_sphere ( UnstructuredMesh &  mesh,
const Real  rad = 1,
const unsigned int  nr = 2,
const ElemType  type = INVALID_ELEM,
const unsigned int  n_smooth = 2,
const bool  flat = true 
)

Meshes a spherical or mapped-spherical domain.

Definition at line 1530 of file mesh_generation.C.

References libMesh::out.

01536 {
01537         libMesh::out << "Building a circle/sphere only works with AMR." << std::endl;
01538         libmesh_error();
01539 }

void libMesh::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.

Boundary ids are set to be equal to the side indexing on a master quad

Definition at line 1499 of file mesh_generation.C.

References build_cube().

01506 {
01507   // This method only makes sense in 2D!
01508   // But we now just turn a non-2D mesh into a 2D mesh
01509   //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
01510 
01511   // Call the build_cube() member to actually do the work for us.
01512   build_cube (mesh,
01513               nx, ny, 0,
01514               xmin, xmax,
01515               ymin, ymax,
01516               0., 0.,
01517               type,
01518               gauss_lobatto_grid);
01519 }

Site Created By: libMesh Developers
Last modified: February 05 2013 19:55:44 UTC

Hosted By:
SourceForge.net Logo