libMesh
Namespaces | Classes | Functions
libMesh::MeshTools::Generation Namespace Reference

Tools for Mesh generation. More...

Namespaces

 Private
 

Classes

class  QueryElemSubdomainIDBase
 Class for receiving the callback during extrusion generation and providing user-defined subdomains based on the old (existing) element id and the current layer. More...
 

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)
 Builds a \( nx \times ny \times nz \) (elements) cube. More...
 
void build_point (UnstructuredMesh &mesh, const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 A specialized build_cube() for 0D meshes. More...
 
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)
 A specialized build_cube() for 1D meshes. More...
 
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)
 A specialized build_cube() for 2D meshes. More...
 
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)
 Meshes a spherical or mapped-spherical domain. More...
 
void build_extrusion (UnstructuredMesh &mesh, const MeshBase &cross_section, const unsigned int nz, RealVectorValue extrusion_vector, QueryElemSubdomainIDBase *elem_subdomain=nullptr)
 Meshes the tensor product of a 1D and a 1D-or-2D domain. More...
 
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=nullptr)
 Meshes a rectangular (2D) region (with or without holes) with a Delaunay triangulation. More...
 
void surface_octahedron (UnstructuredMesh &mesh, Real xmin, Real xmax, Real ymin, Real ymax, Real zmin, Real zmax, bool flip_tris=false)
 Meshes the surface of an octahedron with 8 Tri3 elements, with counter-clockwise (libMesh default) node ordering as viewed from the octahedron exterior if flip_tris is false or from the interior otherwise. More...
 

Detailed Description

Tools for Mesh generation.

Author
Benjamin S. Kirk
Date
2004

Function Documentation

◆ build_cube()

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 313 of file mesh_generation.C.

References libMesh::BoundaryInfo::add_node(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::Elem::build_with_id(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::Utility::enum_to_string(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::DofObject::id(), libMesh::MeshTools::Generation::Private::idx(), libMesh::index_range(), libMesh::INVALID_ELEM, libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_assert(), mesh, libMesh::NODEELEM, libMesh::BoundaryInfo::nodeset_name(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::QUADSHELL8, libMesh::QUADSHELL9, libMesh::Real, libMesh::MeshTools::Modification::redistribute(), libMesh::BoundaryInfo::remove(), libMesh::Elem::set_node(), libMesh::BoundaryInfo::sideset_name(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and libMesh::TRISHELL3.

Referenced by add_cube_convex_hull_to_mesh(), AllSecondOrderTest::allCompleteOrderMixedFixing3D(), AllSecondOrderTest::allSecondOrderMixedFixing3D(), Biharmonic::Biharmonic(), build_line(), ExtraIntegersTest::build_mesh(), build_point(), build_square(), CopyNodesAndElementsTest::collectMeshes(), main(), ParsedFEMFunctionTest::setUp(), PerElemTest< elem_type >::setUp(), RationalMapTest< elem_type >::setUp(), setup(), ParallelGhostSyncTest::setUp(), FETestBase< order, family, elem_type, 1 >::setUp(), SystemsTest::test3DProjectVectorFE(), DistortTest::test_helper_3D(), AllTriTest::test_helper_3D(), MeshFunctionTest::test_p_level(), VolumeTest::test_true_centroid_and_volume(), MeshStitchTest::testAmbiguousRemappingStitch(), SystemsTest::testAssemblyWithDgFemContext(), SystemsTest::testBlockRestrictedVarNDofs(), MeshStitchTest::testBoundaryInfo(), SystemsTest::testBoundaryProjectCube(), MeshGenerationTest::testBuildCube(), DefaultCouplingTest::testCoupling(), PointNeighborCouplingTest::testCoupling(), SystemsTest::testDofCouplingWithVarGroups(), DofMapTest::testDofOwner(), BoundaryInfoTest::testEdgeBoundaryConditions(), MeshInputTest::testExodusWriteElementDataFromDiscontinuousNodalData(), VolumeTest::testHex20PLevelTrueCentroid(), MeshTetTest::testHole(), InfFERadialTest::testInfQuants(), PointLocatorTest::testLocator(), MeshStitchTest::testMeshStitch(), MeshStitchTest::testMeshStitchElemsets(), PartitionerTest< PartitionerSubclass, MeshClass >::testPartition(), SystemsTest::testProjectCube(), SystemsTest::testProjectCubeWithMeshFunction(), SystemsTest::testProjectMatrix3D(), MeshStitchTest::testRemappingStitch(), SystemsTest::testSetSystemParameterOverEquationSystem(), InfFERadialTest::testSides(), InfFERadialTest::testSingleOrder(), VolumeTest::testTwistedVolume(), and WriteEdgesetData::testWriteImpl().

322 {
323  LOG_SCOPE("build_cube()", "MeshTools::Generation");
324 
325  // Declare that we are using the indexing utility routine
326  // in the "Private" part of our current namespace. If this doesn't
327  // work in GCC 2.95.3 we can either remove it or stop supporting
328  // 2.95.3 altogether.
329  // Changing this to import the whole namespace... just importing idx
330  // causes an internal compiler error for Intel Compiler 11.0 on Linux
331  // in debug mode.
332  using namespace MeshTools::Generation::Private;
333 
334  // Clear the mesh and start from scratch
335  mesh.clear();
336 
337  BoundaryInfo & boundary_info = mesh.get_boundary_info();
338 
339  if (nz != 0)
340  {
341  mesh.set_mesh_dimension(3);
342  mesh.set_spatial_dimension(3);
343  }
344  else if (ny != 0)
345  {
346  mesh.set_mesh_dimension(2);
347  mesh.set_spatial_dimension(2);
348  }
349  else if (nx != 0)
350  {
351  mesh.set_mesh_dimension(1);
352  mesh.set_spatial_dimension(1);
353  }
354  else
355  {
356  // Will we get here?
357  mesh.set_mesh_dimension(0);
358  mesh.set_spatial_dimension(0);
359  }
360 
361  switch (mesh.mesh_dimension())
362  {
363  //---------------------------------------------------------------------
364  // Build a 0D point
365  case 0:
366  {
367  libmesh_assert_equal_to (nx, 0);
368  libmesh_assert_equal_to (ny, 0);
369  libmesh_assert_equal_to (nz, 0);
370 
371  libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
372 
373  // Build one nodal element for the mesh
374  mesh.add_point (Point(0, 0, 0), 0);
375  Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
376  elem->set_node(0, mesh.node_ptr(0));
377 
378  break;
379  }
380 
381 
382 
383  //---------------------------------------------------------------------
384  // Build a 1D line
385  case 1:
386  {
387  libmesh_assert_not_equal_to (nx, 0);
388  libmesh_assert_equal_to (ny, 0);
389  libmesh_assert_equal_to (nz, 0);
390  libmesh_assert_less (xmin, xmax);
391 
392  // Reserve elements
393  switch (type)
394  {
395  case INVALID_ELEM:
396  case EDGE2:
397  case EDGE3:
398  case EDGE4:
399  {
400  mesh.reserve_elem (nx);
401  break;
402  }
403 
404  default:
405  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
406  }
407 
408  // Reserve nodes
409  switch (type)
410  {
411  case INVALID_ELEM:
412  case EDGE2:
413  {
414  mesh.reserve_nodes(nx+1);
415  break;
416  }
417 
418  case EDGE3:
419  {
420  mesh.reserve_nodes(2*nx+1);
421  break;
422  }
423 
424  case EDGE4:
425  {
426  mesh.reserve_nodes(3*nx+1);
427  break;
428  }
429 
430  default:
431  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
432  }
433 
434 
435  // Build the nodes, depends on whether we're using linears,
436  // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
437  unsigned int node_id = 0;
438  switch(type)
439  {
440  case INVALID_ELEM:
441  case EDGE2:
442  {
443  for (unsigned int i=0; i<=nx; i++)
444  {
445  const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/nx, 0, 0), node_id++);
446  if (i == 0)
447  boundary_info.add_node(node, 0);
448  if (i == nx)
449  boundary_info.add_node(node, 1);
450  }
451 
452  break;
453  }
454 
455  case EDGE3:
456  {
457  for (unsigned int i=0; i<=2*nx; i++)
458  {
459  const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(2*nx), 0, 0), node_id++);
460  if (i == 0)
461  boundary_info.add_node(node, 0);
462  if (i == 2*nx)
463  boundary_info.add_node(node, 1);
464  }
465  break;
466  }
467 
468  case EDGE4:
469  {
470  for (unsigned int i=0; i<=3*nx; i++)
471  {
472  const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(3*nx), 0, 0), node_id++);
473  if (i == 0)
474  boundary_info.add_node(node, 0);
475  if (i == 3*nx)
476  boundary_info.add_node(node, 1);
477  }
478 
479  break;
480  }
481 
482  default:
483  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
484 
485  }
486 
487  // Build the elements of the mesh
488  switch(type)
489  {
490  case INVALID_ELEM:
491  case EDGE2:
492  {
493  for (unsigned int i=0; i<nx; i++)
494  {
495  Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE2, i));
496  elem->set_node(0, mesh.node_ptr(i));
497  elem->set_node(1, mesh.node_ptr(i+1));
498 
499  if (i == 0)
500  boundary_info.add_side(elem, 0, 0);
501 
502  if (i == (nx-1))
503  boundary_info.add_side(elem, 1, 1);
504 
505  }
506  break;
507  }
508 
509  case EDGE3:
510  {
511  for (unsigned int i=0; i<nx; i++)
512  {
513  Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE3, i));
514  elem->set_node(0, mesh.node_ptr(2*i));
515  elem->set_node(2, mesh.node_ptr(2*i+1));
516  elem->set_node(1, mesh.node_ptr(2*i+2));
517 
518  if (i == 0)
519  boundary_info.add_side(elem, 0, 0);
520 
521  if (i == (nx-1))
522  boundary_info.add_side(elem, 1, 1);
523  }
524  break;
525  }
526 
527  case EDGE4:
528  {
529  for (unsigned int i=0; i<nx; i++)
530  {
531  Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE4, i));
532  elem->set_node(0, mesh.node_ptr(3*i));
533  elem->set_node(2, mesh.node_ptr(3*i+1));
534  elem->set_node(3, mesh.node_ptr(3*i+2));
535  elem->set_node(1, mesh.node_ptr(3*i+3));
536 
537  if (i == 0)
538  boundary_info.add_side(elem, 0, 0);
539 
540  if (i == (nx-1))
541  boundary_info.add_side(elem, 1, 1);
542  }
543  break;
544  }
545 
546  default:
547  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
548  }
549 
550  // Move the nodes to their final locations.
551  if (gauss_lobatto_grid)
552  {
553  GaussLobattoRedistributionFunction func(nx, xmin, xmax);
554  MeshTools::Modification::redistribute(mesh, func);
555  }
556  else // !gauss_lobatto_grid
557  {
558  for (Node * node : mesh.node_ptr_range())
559  (*node)(0) = (*node)(0)*(xmax-xmin) + xmin;
560  }
561 
562  // Add sideset names to boundary info
563  boundary_info.sideset_name(0) = "left";
564  boundary_info.sideset_name(1) = "right";
565 
566  // Add nodeset names to boundary info
567  boundary_info.nodeset_name(0) = "left";
568  boundary_info.nodeset_name(1) = "right";
569 
570  break;
571  }
572 
573 
574 
575 
576 
577 
578 
579 
580 
581 
582  //---------------------------------------------------------------------
583  // Build a 2D quadrilateral
584  case 2:
585  {
586  libmesh_assert_not_equal_to (nx, 0);
587  libmesh_assert_not_equal_to (ny, 0);
588  libmesh_assert_equal_to (nz, 0);
589  libmesh_assert_less (xmin, xmax);
590  libmesh_assert_less (ymin, ymax);
591 
592  // Reserve elements. The TRI3 and TRI6 meshes
593  // have twice as many elements...
594  switch (type)
595  {
596  case INVALID_ELEM:
597  case QUAD4:
598  case QUADSHELL4:
599  case QUAD8:
600  case QUADSHELL8:
601  case QUAD9:
602  case QUADSHELL9:
603  {
604  mesh.reserve_elem (nx*ny);
605  break;
606  }
607 
608  case TRI3:
609  case TRISHELL3:
610  case TRI6:
611  case TRI7:
612  {
613  mesh.reserve_elem (2*nx*ny);
614  break;
615  }
616 
617  default:
618  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
619  }
620 
621 
622 
623  // Reserve nodes. The quadratic element types
624  // need to reserve more nodes than the linear types.
625  switch (type)
626  {
627  case INVALID_ELEM:
628  case QUAD4:
629  case QUADSHELL4:
630  case TRI3:
631  case TRISHELL3:
632  {
633  mesh.reserve_nodes( (nx+1)*(ny+1) );
634  break;
635  }
636 
637  case QUAD8:
638  case QUADSHELL8:
639  case QUAD9:
640  case QUADSHELL9:
641  case TRI6:
642  {
643  mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
644  break;
645  }
646 
647  case TRI7:
648  {
649  mesh.reserve_nodes( (2*nx+1)*(2*ny+1) + 2*nx*ny );
650  break;
651  }
652 
653  default:
654  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
655  }
656 
657 
658 
659  // Build the nodes. Depends on whether you are using a linear
660  // or quadratic element, and whether you are using a uniform
661  // grid or the Gauss-Lobatto grid points.
662  unsigned int node_id = 0;
663  switch (type)
664  {
665  case INVALID_ELEM:
666  case QUAD4:
667  case QUADSHELL4:
668  case TRI3:
669  case TRISHELL3:
670  {
671  for (unsigned int j=0; j<=ny; j++)
672  for (unsigned int i=0; i<=nx; i++)
673  {
674  const Node * const node =
675  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
676  static_cast<Real>(j) / static_cast<Real>(ny),
677  0.),
678  node_id++);
679  if (j == 0)
680  boundary_info.add_node(node, 0);
681  if (j == ny)
682  boundary_info.add_node(node, 2);
683  if (i == 0)
684  boundary_info.add_node(node, 3);
685  if (i == nx)
686  boundary_info.add_node(node, 1);
687  }
688 
689  break;
690  }
691 
692  case QUAD8:
693  case QUADSHELL8:
694  case QUAD9:
695  case QUADSHELL9:
696  case TRI6:
697  case TRI7:
698  {
699  for (unsigned int j=0; j<=(2*ny); j++)
700  for (unsigned int i=0; i<=(2*nx); i++)
701  {
702  const Node * const node =
703  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
704  static_cast<Real>(j) / static_cast<Real>(2 * ny),
705  0),
706  node_id++);
707  if (j == 0)
708  boundary_info.add_node(node, 0);
709  if (j == 2*ny)
710  boundary_info.add_node(node, 2);
711  if (i == 0)
712  boundary_info.add_node(node, 3);
713  if (i == 2*nx)
714  boundary_info.add_node(node, 1);
715  }
716 
717  // We'll add any interior Tri7 nodes last, to keep from
718  // messing with our idx function
719  if (type == TRI7)
720  for (unsigned int j=0; j<(3*ny); j += 3)
721  for (unsigned int i=0; i<(3*nx); i += 3)
722  {
723  // The bottom-right triangle's center node
724  mesh.add_point(Point(static_cast<Real>(i+2) / static_cast<Real>(3 * nx),
725  static_cast<Real>(j+1) / static_cast<Real>(3 * ny),
726  0),
727  node_id++);
728  // The top-left triangle's center node
729  mesh.add_point(Point(static_cast<Real>(i+1) / static_cast<Real>(3 * nx),
730  static_cast<Real>(j+2) / static_cast<Real>(3 * ny),
731  0),
732  node_id++);
733  }
734 
735  break;
736  }
737 
738 
739  default:
740  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
741  }
742 
743 
744 
745 
746 
747 
748  // Build the elements. Each one is a bit different.
749  unsigned int elem_id = 0;
750  switch (type)
751  {
752 
753  case INVALID_ELEM:
754  case QUAD4:
755  case QUADSHELL4:
756  {
757  for (unsigned int j=0; j<ny; j++)
758  for (unsigned int i=0; i<nx; i++)
759  {
760  Elem * elem = mesh.add_elem(Elem::build_with_id(type == INVALID_ELEM ? QUAD4 : type, elem_id++));
761  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
762  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
763  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
764  elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+1) ));
765 
766  if (j == 0)
767  boundary_info.add_side(elem, 0, 0);
768 
769  if (j == (ny-1))
770  boundary_info.add_side(elem, 2, 2);
771 
772  if (i == 0)
773  boundary_info.add_side(elem, 3, 3);
774 
775  if (i == (nx-1))
776  boundary_info.add_side(elem, 1, 1);
777  }
778  break;
779  }
780 
781 
782  case TRI3:
783  case TRISHELL3:
784  {
785  for (unsigned int j=0; j<ny; j++)
786  for (unsigned int i=0; i<nx; i++)
787  {
788  // Add first Tri3
789  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
790  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
791  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
792  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
793 
794  if (j == 0)
795  boundary_info.add_side(elem, 0, 0);
796 
797  if (i == (nx-1))
798  boundary_info.add_side(elem, 1, 1);
799 
800  // Add second Tri3
801  elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
802  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
803  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j+1)));
804  elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+1) ));
805 
806  if (j == (ny-1))
807  boundary_info.add_side(elem, 1, 2);
808 
809  if (i == 0)
810  boundary_info.add_side(elem, 2, 3);
811  }
812  break;
813  }
814 
815 
816 
817  case QUAD8:
818  case QUADSHELL8:
819  case QUAD9:
820  case QUADSHELL9:
821  {
822  for (unsigned int j=0; j<(2*ny); j += 2)
823  for (unsigned int i=0; i<(2*nx); i += 2)
824  {
825  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
826  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
827  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
828  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
829  elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+2) ));
830  elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j) ));
831  elem->set_node(5, mesh.node_ptr(idx(type,nx,i+2,j+1)));
832  elem->set_node(6, mesh.node_ptr(idx(type,nx,i+1,j+2)));
833  elem->set_node(7, mesh.node_ptr(idx(type,nx,i,j+1) ));
834 
835  if (type == QUAD9 || type == QUADSHELL9)
836  elem->set_node(8, mesh.node_ptr(idx(type,nx,i+1,j+1)));
837 
838  if (j == 0)
839  boundary_info.add_side(elem, 0, 0);
840 
841  if (j == 2*(ny-1))
842  boundary_info.add_side(elem, 2, 2);
843 
844  if (i == 0)
845  boundary_info.add_side(elem, 3, 3);
846 
847  if (i == 2*(nx-1))
848  boundary_info.add_side(elem, 1, 1);
849  }
850  break;
851  }
852 
853 
854  case TRI6:
855  case TRI7:
856  {
857  for (unsigned int j=0; j<(2*ny); j += 2)
858  for (unsigned int i=0; i<(2*nx); i += 2)
859  {
860  // Add first Tri in the bottom-right of its quad
861  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
862  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
863  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
864  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
865  elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j) ));
866  elem->set_node(4, mesh.node_ptr(idx(type,nx,i+2,j+1)));
867  elem->set_node(5, mesh.node_ptr(idx(type,nx,i+1,j+1)));
868 
869  if (type == TRI7)
870  elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
871 
872  if (j == 0)
873  boundary_info.add_side(elem, 0, 0);
874 
875  if (i == 2*(nx-1))
876  boundary_info.add_side(elem, 1, 1);
877 
878  // Add second Tri in the top left of its quad
879  elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
880  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
881  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j+2)));
882  elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+2) ));
883  elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j+1)));
884  elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j+2)));
885  elem->set_node(5, mesh.node_ptr(idx(type,nx,i,j+1) ));
886 
887  if (type == TRI7)
888  elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
889 
890  if (j == 2*(ny-1))
891  boundary_info.add_side(elem, 1, 2);
892 
893  if (i == 0)
894  boundary_info.add_side(elem, 2, 3);
895  }
896  break;
897  };
898 
899 
900  default:
901  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
902  }
903 
904 
905 
906 
907  // Scale the nodal positions
908  if (gauss_lobatto_grid)
909  {
910  GaussLobattoRedistributionFunction func(nx, xmin, xmax,
911  ny, ymin, ymax);
912  MeshTools::Modification::redistribute(mesh, func);
913  }
914  else // !gauss_lobatto_grid
915  {
916  for (Node * node : mesh.node_ptr_range())
917  {
918  (*node)(0) = ((*node)(0))*(xmax-xmin) + xmin;
919  (*node)(1) = ((*node)(1))*(ymax-ymin) + ymin;
920  }
921  }
922 
923  // Add sideset names to boundary info
924  boundary_info.sideset_name(0) = "bottom";
925  boundary_info.sideset_name(1) = "right";
926  boundary_info.sideset_name(2) = "top";
927  boundary_info.sideset_name(3) = "left";
928 
929  // Add nodeset names to boundary info
930  boundary_info.nodeset_name(0) = "bottom";
931  boundary_info.nodeset_name(1) = "right";
932  boundary_info.nodeset_name(2) = "top";
933  boundary_info.nodeset_name(3) = "left";
934 
935  break;
936  }
937 
938 
939 
940 
941 
942 
943 
944 
945 
946 
947 
948  //---------------------------------------------------------------------
949  // Build a 3D mesh using hexes, tets, prisms, or pyramids.
950  case 3:
951  {
952  libmesh_assert_not_equal_to (nx, 0);
953  libmesh_assert_not_equal_to (ny, 0);
954  libmesh_assert_not_equal_to (nz, 0);
955  libmesh_assert_less (xmin, xmax);
956  libmesh_assert_less (ymin, ymax);
957  libmesh_assert_less (zmin, zmax);
958 
959 
960  // Reserve elements. Meshes with prismatic elements require
961  // twice as many elements.
962  switch (type)
963  {
964  case INVALID_ELEM:
965  case HEX8:
966  case HEX20:
967  case HEX27:
968  case TET4: // TET4's are created from an initial HEX27 discretization
969  case TET10: // TET10's are created from an initial HEX27 discretization
970  case TET14: // TET14's are created from an initial HEX27 discretization
971  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
972  case PYRAMID13:
973  case PYRAMID14:
974  case PYRAMID18:
975  {
976  mesh.reserve_elem(nx*ny*nz);
977  break;
978  }
979 
980  case PRISM6:
981  case PRISM15:
982  case PRISM18:
983  case PRISM20:
984  case PRISM21:
985  {
986  mesh.reserve_elem(2*nx*ny*nz);
987  break;
988  }
989 
990  default:
991  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
992  }
993 
994 
995 
996 
997 
998  // Reserve nodes. Quadratic elements need twice as many nodes as linear elements.
999  switch (type)
1000  {
1001  case INVALID_ELEM:
1002  case HEX8:
1003  case PRISM6:
1004  {
1005  mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
1006  break;
1007  }
1008 
1009  case HEX20:
1010  case HEX27:
1011  case TET4: // TET4's are created from an initial HEX27 discretization
1012  case TET10: // TET10's are created from an initial HEX27 discretization
1013  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1014  case PYRAMID13:
1015  case PYRAMID14:
1016  case PYRAMID18:
1017  case PRISM15:
1018  case PRISM18:
1019  {
1020  // FYI: The resulting TET4 mesh will have exactly
1021  // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
1022  // nodes once the additional mid-edge nodes for the HEX27 discretization
1023  // have been deleted.
1024  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
1025  break;
1026  }
1027 
1028  case TET14:
1029  {
1030  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1031  24*nx*ny*nz +
1032  4*(nx*ny + ny*nz + nx*nz) );
1033  break;
1034  }
1035 
1036  case PRISM20:
1037  {
1038  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1039  2*nx*ny*(nz+1) );
1040  break;
1041  }
1042 
1043  case PRISM21:
1044  {
1045  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1046  2*nx*ny*(2*nz+1) );
1047  break;
1048  }
1049 
1050  default:
1051  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1052  }
1053 
1054 
1055 
1056 
1057  // Build the nodes.
1058  unsigned int node_id = 0;
1059  switch (type)
1060  {
1061  case INVALID_ELEM:
1062  case HEX8:
1063  case PRISM6:
1064  {
1065  for (unsigned int k=0; k<=nz; k++)
1066  for (unsigned int j=0; j<=ny; j++)
1067  for (unsigned int i=0; i<=nx; i++)
1068  {
1069  const Node * const node =
1070  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
1071  static_cast<Real>(j) / static_cast<Real>(ny),
1072  static_cast<Real>(k) / static_cast<Real>(nz)),
1073  node_id++);
1074  if (k == 0)
1075  boundary_info.add_node(node, 0);
1076  if (k == nz)
1077  boundary_info.add_node(node, 5);
1078  if (j == 0)
1079  boundary_info.add_node(node, 1);
1080  if (j == ny)
1081  boundary_info.add_node(node, 3);
1082  if (i == 0)
1083  boundary_info.add_node(node, 4);
1084  if (i == nx)
1085  boundary_info.add_node(node, 2);
1086  }
1087 
1088  break;
1089  }
1090 
1091  case HEX20:
1092  case HEX27:
1093  case TET4: // TET4's are created from an initial HEX27 discretization
1094  case TET10: // TET10's are created from an initial HEX27 discretization
1095  case TET14: // TET14's are created from an initial HEX27 discretization
1096  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1097  case PYRAMID13:
1098  case PYRAMID14:
1099  case PYRAMID18:
1100  case PRISM15:
1101  case PRISM18:
1102  case PRISM20:
1103  case PRISM21:
1104  {
1105  for (unsigned int k=0; k<=(2*nz); k++)
1106  for (unsigned int j=0; j<=(2*ny); j++)
1107  for (unsigned int i=0; i<=(2*nx); i++)
1108  {
1109  const Node * const node =
1110  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
1111  static_cast<Real>(j) / static_cast<Real>(2 * ny),
1112  static_cast<Real>(k) / static_cast<Real>(2 * nz)),
1113  node_id++);
1114  if (k == 0)
1115  boundary_info.add_node(node, 0);
1116  if (k == 2*nz)
1117  boundary_info.add_node(node, 5);
1118  if (j == 0)
1119  boundary_info.add_node(node, 1);
1120  if (j == 2*ny)
1121  boundary_info.add_node(node, 3);
1122  if (i == 0)
1123  boundary_info.add_node(node, 4);
1124  if (i == 2*nx)
1125  boundary_info.add_node(node, 2);
1126  }
1127 
1128  if (type == PRISM20 ||
1129  type == PRISM21)
1130  {
1131  const unsigned int kmax = (type == PRISM20) ? nz : 2*nz;
1132  for (unsigned int k=0; k<=kmax; k++)
1133  for (unsigned int j=0; j<ny; j++)
1134  for (unsigned int i=0; i<nx; i++)
1135  {
1136  const Node * const node1 =
1137  mesh.add_point(Point((static_cast<Real>(i)+1/Real(3)) / static_cast<Real>(nx),
1138  (static_cast<Real>(j)+1/Real(3)) / static_cast<Real>(ny),
1139  static_cast<Real>(k) / static_cast<Real>(kmax)),
1140  node_id++);
1141  if (k == 0)
1142  boundary_info.add_node(node1, 0);
1143  if (k == kmax)
1144  boundary_info.add_node(node1, 5);
1145 
1146  const Node * const node2 =
1147  mesh.add_point(Point((static_cast<Real>(i)+2/Real(3)) / static_cast<Real>(nx),
1148  (static_cast<Real>(j)+2/Real(3)) / static_cast<Real>(ny),
1149  static_cast<Real>(k) / static_cast<Real>(kmax)),
1150  node_id++);
1151  if (k == 0)
1152  boundary_info.add_node(node2, 0);
1153  if (k == kmax)
1154  boundary_info.add_node(node2, 5);
1155  }
1156  }
1157 
1158  break;
1159  }
1160 
1161 
1162  default:
1163  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1164  }
1165 
1166 
1167 
1168 
1169  // Build the elements.
1170  unsigned int elem_id = 0;
1171  switch (type)
1172  {
1173  case INVALID_ELEM:
1174  case HEX8:
1175  {
1176  for (unsigned int k=0; k<nz; k++)
1177  for (unsigned int j=0; j<ny; j++)
1178  for (unsigned int i=0; i<nx; i++)
1179  {
1180  Elem * elem = mesh.add_elem(Elem::build_with_id(HEX8, elem_id++));
1181  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1182  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1183  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1184  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1185  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1186  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1187  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1188  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1189 
1190  if (k == 0)
1191  boundary_info.add_side(elem, 0, 0);
1192 
1193  if (k == (nz-1))
1194  boundary_info.add_side(elem, 5, 5);
1195 
1196  if (j == 0)
1197  boundary_info.add_side(elem, 1, 1);
1198 
1199  if (j == (ny-1))
1200  boundary_info.add_side(elem, 3, 3);
1201 
1202  if (i == 0)
1203  boundary_info.add_side(elem, 4, 4);
1204 
1205  if (i == (nx-1))
1206  boundary_info.add_side(elem, 2, 2);
1207  }
1208  break;
1209  }
1210 
1211 
1212 
1213 
1214  case PRISM6:
1215  {
1216  for (unsigned int k=0; k<nz; k++)
1217  for (unsigned int j=0; j<ny; j++)
1218  for (unsigned int i=0; i<nx; i++)
1219  {
1220  // First Prism
1221  Elem * elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1222  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1223  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1224  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1225  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1226  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1227  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1228 
1229  // Add sides for first prism to boundary info object
1230  if (i==0)
1231  boundary_info.add_side(elem, 3, 4);
1232 
1233  if (j==0)
1234  boundary_info.add_side(elem, 1, 1);
1235 
1236  if (k==0)
1237  boundary_info.add_side(elem, 0, 0);
1238 
1239  if (k == (nz-1))
1240  boundary_info.add_side(elem, 4, 5);
1241 
1242  // Second Prism
1243  elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1244  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1245  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1246  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1247  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1248  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1249  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1250 
1251  // Add sides for second prism to boundary info object
1252  if (i == (nx-1))
1253  boundary_info.add_side(elem, 1, 2);
1254 
1255  if (j == (ny-1))
1256  boundary_info.add_side(elem, 2, 3);
1257 
1258  if (k==0)
1259  boundary_info.add_side(elem, 0, 0);
1260 
1261  if (k == (nz-1))
1262  boundary_info.add_side(elem, 4, 5);
1263  }
1264  break;
1265  }
1266 
1267 
1268 
1269 
1270 
1271 
1272  case HEX20:
1273  case HEX27:
1274  case TET4: // TET4's are created from an initial HEX27 discretization
1275  case TET10: // TET10's are created from an initial HEX27 discretization
1276  case TET14: // TET14's are created from an initial HEX27 discretization
1277  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1278  case PYRAMID13:
1279  case PYRAMID14:
1280  case PYRAMID18:
1281  {
1282  for (unsigned int k=0; k<(2*nz); k += 2)
1283  for (unsigned int j=0; j<(2*ny); j += 2)
1284  for (unsigned int i=0; i<(2*nx); i += 2)
1285  {
1286  ElemType build_type = (type == HEX20) ? HEX20 : HEX27;
1287  Elem * elem = mesh.add_elem(Elem::build_with_id(build_type, elem_id++));
1288 
1289  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1290  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1291  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1292  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1293  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1294  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1295  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2)));
1296  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1297  elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1298  elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1299  elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1300  elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1301  elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1302  elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1303  elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1304  elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1305  elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1306  elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1307  elem->set_node(18, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1308  elem->set_node(19, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1309 
1310  if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == TET14) ||
1311  (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14) ||
1312  (type == PYRAMID18))
1313  {
1314  elem->set_node(20, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1315  elem->set_node(21, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1316  elem->set_node(22, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1317  elem->set_node(23, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1318  elem->set_node(24, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1319  elem->set_node(25, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1320  elem->set_node(26, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1321  }
1322 
1323  if (k == 0)
1324  boundary_info.add_side(elem, 0, 0);
1325 
1326  if (k == 2*(nz-1))
1327  boundary_info.add_side(elem, 5, 5);
1328 
1329  if (j == 0)
1330  boundary_info.add_side(elem, 1, 1);
1331 
1332  if (j == 2*(ny-1))
1333  boundary_info.add_side(elem, 3, 3);
1334 
1335  if (i == 0)
1336  boundary_info.add_side(elem, 4, 4);
1337 
1338  if (i == 2*(nx-1))
1339  boundary_info.add_side(elem, 2, 2);
1340  }
1341  break;
1342  }
1343 
1344 
1345 
1346 
1347  case PRISM15:
1348  case PRISM18:
1349  case PRISM20:
1350  case PRISM21:
1351  {
1352  for (unsigned int k=0; k<(2*nz); k += 2)
1353  for (unsigned int j=0; j<(2*ny); j += 2)
1354  for (unsigned int i=0; i<(2*nx); i += 2)
1355  {
1356  // First Prism
1357  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1358  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1359  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1360  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1361  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1362  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1363  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1364  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1365  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1366  elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1367  elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1368  elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1369  elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1370  elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1371  elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1372  elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1373 
1374  if (type == PRISM18 ||
1375  type == PRISM20 ||
1376  type == PRISM21)
1377  {
1378  elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1379  elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1380  elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1381  }
1382 
1383  if (type == PRISM20)
1384  {
1385  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1386  elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2));
1387  elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2));
1388  }
1389 
1390  if (type == PRISM21)
1391  {
1392  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1393  elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2));
1394  elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2));
1395  elem->set_node(20, mesh.node_ptr(base_idx+((k+1)*(nx*ny)+j/2*nx+i/2)*2));
1396  }
1397 
1398  // Add sides for first prism to boundary info object
1399  if (i==0)
1400  boundary_info.add_side(elem, 3, 4);
1401 
1402  if (j==0)
1403  boundary_info.add_side(elem, 1, 1);
1404 
1405  if (k==0)
1406  boundary_info.add_side(elem, 0, 0);
1407 
1408  if (k == 2*(nz-1))
1409  boundary_info.add_side(elem, 4, 5);
1410 
1411 
1412  // Second Prism
1413  elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1414  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+2,j,k) ));
1415  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1416  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+2,k) ));
1417  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) ));
1418  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) ));
1419  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) ));
1420  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1421  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1422  elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1423  elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) ));
1424  elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1425  elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) ));
1426  elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1427  elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1428  elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1429 
1430  if (type == PRISM18 ||
1431  type == PRISM20 ||
1432  type == PRISM21)
1433  {
1434  elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1435  elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1436  elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1437  }
1438 
1439  if (type == PRISM20)
1440  {
1441  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1442  elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2+1));
1443  elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1444  }
1445 
1446  if (type == PRISM21)
1447  {
1448  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1449  elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2+1));
1450  elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2+1));
1451  elem->set_node(20, mesh.node_ptr(base_idx+((k+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1452  }
1453 
1454  // Add sides for second prism to boundary info object
1455  if (i == 2*(nx-1))
1456  boundary_info.add_side(elem, 1, 2);
1457 
1458  if (j == 2*(ny-1))
1459  boundary_info.add_side(elem, 2, 3);
1460 
1461  if (k==0)
1462  boundary_info.add_side(elem, 0, 0);
1463 
1464  if (k == 2*(nz-1))
1465  boundary_info.add_side(elem, 4, 5);
1466 
1467  }
1468  break;
1469  }
1470 
1471 
1472 
1473 
1474 
1475  default:
1476  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1477  }
1478 
1479 
1480 
1481 
1482  //.......................................
1483  // Scale the nodal positions
1484  if (gauss_lobatto_grid)
1485  {
1486  GaussLobattoRedistributionFunction func(nx, xmin, xmax,
1487  ny, ymin, ymax,
1488  nz, zmin, zmax);
1489  MeshTools::Modification::redistribute(mesh, func);
1490  }
1491  else // !gauss_lobatto_grid
1492  {
1493  for (unsigned int p=0; p<mesh.n_nodes(); p++)
1494  {
1495  mesh.node_ref(p)(0) = (mesh.node_ref(p)(0))*(xmax-xmin) + xmin;
1496  mesh.node_ref(p)(1) = (mesh.node_ref(p)(1))*(ymax-ymin) + ymin;
1497  mesh.node_ref(p)(2) = (mesh.node_ref(p)(2))*(zmax-zmin) + zmin;
1498  }
1499  }
1500 
1501 
1502 
1503  // Additional work for tets and pyramids: we take the existing
1504  // HEX27 discretization and split each element into 24
1505  // sub-tets or 6 sub-pyramids.
1506  //
1507  // 24 isn't the minimum-possible number of tets, but it
1508  // obviates any concerns about the edge orientations between
1509  // the various elements.
1510  if ((type == TET4) ||
1511  (type == TET10) ||
1512  (type == TET14) ||
1513  (type == PYRAMID5) ||
1514  (type == PYRAMID13) ||
1515  (type == PYRAMID14) ||
1516  (type == PYRAMID18))
1517  {
1518  // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
1519  std::vector<std::unique_ptr<Elem>> new_elements;
1520 
1521  // For avoiding extraneous construction of element sides
1522  std::unique_ptr<Elem> side;
1523 
1524  if ((type == TET4) || (type == TET10) || (type == TET14))
1525  new_elements.reserve(24*mesh.n_elem());
1526  else
1527  new_elements.reserve(6*mesh.n_elem());
1528 
1529  // Create tetrahedra or pyramids
1530  for (auto & base_hex : mesh.element_ptr_range())
1531  {
1532  // Get a pointer to the node located at the HEX27 center
1533  Node * apex_node = base_hex->node_ptr(26);
1534 
1535  // Container to catch ids handed back from BoundaryInfo
1536  std::vector<boundary_id_type> ids;
1537 
1538  for (auto s : base_hex->side_index_range())
1539  {
1540  // Get the boundary ID(s) for this side
1541  boundary_info.boundary_ids(base_hex, s, ids);
1542 
1543  // We're creating this Mesh, so there should be 0 or 1 boundary IDs.
1544  libmesh_assert(ids.size() <= 1);
1545 
1546  // A convenient name for the side's ID.
1547  boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];
1548 
1549  // Need to build the full-ordered side!
1550  base_hex->build_side_ptr(side, s);
1551 
1552  if ((type == TET4) || (type == TET10) || (type == TET14))
1553  {
1554  // Build 4 sub-tets per side
1555  for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1556  {
1557  new_elements.push_back( Elem::build(TET4) );
1558  auto & sub_elem = new_elements.back();
1559  sub_elem->set_node(0, side->node_ptr(sub_tet));
1560  sub_elem->set_node(1, side->node_ptr(8)); // center of the face
1561  sub_elem->set_node(2, side->node_ptr(sub_tet==3 ? 0 : sub_tet+1 )); // wrap-around
1562  sub_elem->set_node(3, apex_node); // apex node always used!
1563 
1564  // If the original hex was a boundary hex, add the new sub_tet's side
1565  // 0 with the same b_id. Note: the tets are all aligned so that their
1566  // side 0 is on the boundary.
1567  if (b_id != BoundaryInfo::invalid_id)
1568  boundary_info.add_side(sub_elem.get(), 0, b_id);
1569  }
1570  } // end if ((type == TET4) || (type == TET10) || (type == TET14))
1571 
1572  else // type==PYRAMID*
1573  {
1574  // Build 1 sub-pyramid per side.
1575  new_elements.push_back( Elem::build(PYRAMID5) );
1576  auto & sub_elem = new_elements.back();
1577 
1578  // Set the base. Note that since the apex is *inside* the base_hex,
1579  // and the pyramid uses a counter-clockwise base numbering, we need to
1580  // reverse the [1] and [3] node indices.
1581  sub_elem->set_node(0, side->node_ptr(0));
1582  sub_elem->set_node(1, side->node_ptr(3));
1583  sub_elem->set_node(2, side->node_ptr(2));
1584  sub_elem->set_node(3, side->node_ptr(1));
1585 
1586  // Set the apex
1587  sub_elem->set_node(4, apex_node);
1588 
1589  // If the original hex was a boundary hex, add the new sub_pyr's side
1590  // 4 (the square base) with the same b_id.
1591  if (b_id != BoundaryInfo::invalid_id)
1592  boundary_info.add_side(sub_elem.get(), 4, b_id);
1593  } // end else type==PYRAMID*
1594  }
1595  }
1596 
1597 
1598  // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1599  for (auto & elem : mesh.element_ptr_range())
1600  {
1601  boundary_info.remove(elem); // Safe even if elem has no boundary info.
1602  mesh.delete_elem(elem);
1603  }
1604 
1605  // Add the new elements
1606  for (auto i : index_range(new_elements))
1607  {
1608  new_elements[i]->set_id(i);
1609  mesh.add_elem( std::move(new_elements[i]) );
1610  }
1611 
1612  } // end if (type == TET*,PYRAMID*)
1613 
1614 
1615  // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
1616  if ((type == TET10) || (type == PYRAMID14))
1617  mesh.all_second_order();
1618 
1619  else if (type == PYRAMID13)
1620  mesh.all_second_order(/*full_ordered=*/false);
1621 
1622  else if ((type == TET14) || (type == PYRAMID18))
1623  mesh.all_complete_order();
1624 
1625 
1626  // Add sideset names to boundary info (Z axis out of the screen)
1627  boundary_info.sideset_name(0) = "back";
1628  boundary_info.sideset_name(1) = "bottom";
1629  boundary_info.sideset_name(2) = "right";
1630  boundary_info.sideset_name(3) = "top";
1631  boundary_info.sideset_name(4) = "left";
1632  boundary_info.sideset_name(5) = "front";
1633 
1634  // Add nodeset names to boundary info
1635  boundary_info.nodeset_name(0) = "back";
1636  boundary_info.nodeset_name(1) = "bottom";
1637  boundary_info.nodeset_name(2) = "right";
1638  boundary_info.nodeset_name(3) = "top";
1639  boundary_info.nodeset_name(4) = "left";
1640  boundary_info.nodeset_name(5) = "front";
1641 
1642  break;
1643  } // end case dim==3
1644 
1645  default:
1646  libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension());
1647  }
1648 
1649  // Done building the mesh. Now prepare it for use.
1650  mesh.prepare_for_use ();
1651 }
libMesh::MeshTools::Generation::Private::idx
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int ny, const unsigned int i, const unsigned int j, const unsigned int k)
Definition: mesh_generation.C:115
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::boundary_id_type
int8_t boundary_id_type
Definition: id_types.h:51
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::MeshTools::Modification::redistribute
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
Definition: mesh_modification.C:280
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67

◆ build_delaunay_square()

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 = nullptr 
)

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 2642 of file mesh_generation.C.

References libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::TriangulatorInterface::attach_hole_list(), libMesh::MeshBase::clear(), libMesh::TriangulatorInterface::desired_area(), libMesh::TriangulatorInterface::elem_type(), libMesh::MeshBase::get_boundary_info(), mesh, libMesh::MeshBase::n_nodes(), libMesh::TriangulatorInterface::PSLG, libMesh::Real, libMesh::MeshBase::set_mesh_dimension(), libMesh::TOLERANCE, libMesh::TriangleInterface::triangulate(), and libMesh::TriangulatorInterface::triangulation_type().

2649 {
2650  // Check for reasonable size
2651  libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2652  libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2653  libmesh_assert_less (xmin, xmax);
2654  libmesh_assert_less (ymin, ymax);
2655 
2656  // Clear out any data which may have been in the Mesh
2657  mesh.clear();
2658 
2659  BoundaryInfo & boundary_info = mesh.get_boundary_info();
2660 
2661  // Make sure the new Mesh will be 2D
2662  mesh.set_mesh_dimension(2);
2663 
2664  // The x and y spacing between boundary points
2665  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2666  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2667 
2668  // Bottom
2669  for (unsigned int p=0; p<=nx; ++p)
2670  mesh.add_point(Point(xmin + p*delta_x, ymin));
2671 
2672  // Right side
2673  for (unsigned int p=1; p<ny; ++p)
2674  mesh.add_point(Point(xmax, ymin + p*delta_y));
2675 
2676  // Top
2677  for (unsigned int p=0; p<=nx; ++p)
2678  mesh.add_point(Point(xmax - p*delta_x, ymax));
2679 
2680  // Left side
2681  for (unsigned int p=1; p<ny; ++p)
2682  mesh.add_point(Point(xmin, ymax - p*delta_y));
2683 
2684  // Be sure we added as many points as we thought we did
2685  libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2686 
2687  // Construct the Triangle Interface object
2688  TriangleInterface t(mesh);
2689 
2690  // Set custom variables for the triangulation
2691  t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2692  t.triangulation_type() = TriangleInterface::PSLG;
2693  t.elem_type() = type;
2694 
2695  if (holes != nullptr)
2696  t.attach_hole_list(holes);
2697 
2698  // Triangulate!
2699  t.triangulate();
2700 
2701  // For avoiding extraneous side element construction
2702  std::unique_ptr<const Elem> side;
2703 
2704  // The mesh is now generated, but we still need to mark the boundaries
2705  // to be consistent with the other build_square routines. Note that all
2706  // hole boundary elements get the same ID, 4.
2707  for (auto & elem : mesh.element_ptr_range())
2708  for (auto s : elem->side_index_range())
2709  if (elem->neighbor_ptr(s) == nullptr)
2710  {
2711  elem->build_side_ptr(side, s);
2712 
2713  // Check the location of the side's midpoint. Since
2714  // the square has straight sides, the midpoint is not
2715  // on the corner and thus it is uniquely on one of the
2716  // sides.
2717  Point side_midpoint= 0.5f*( side->point(0) + side->point(1) );
2718 
2719  // The boundary ids are set following the same convention as Quad4 sides
2720  // bottom = 0
2721  // right = 1
2722  // top = 2
2723  // left = 3
2724  // hole = 4
2725  boundary_id_type bc_id=4;
2726 
2727  // bottom
2728  if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2729  bc_id=0;
2730 
2731  // right
2732  else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2733  bc_id=1;
2734 
2735  // top
2736  else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2737  bc_id=2;
2738 
2739  // left
2740  else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2741  bc_id=3;
2742 
2743  // If the point is not on any of the external boundaries, it
2744  // is on one of the holes....
2745 
2746  // Finally, add this element's information to the boundary info object.
2747  boundary_info.add_side(elem->id(), s, bc_id);
2748  }
2749 
2750 } // end build_delaunay_square
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::boundary_id_type
int8_t boundary_id_type
Definition: id_types.h:51
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144

◆ build_extrusion()

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

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

Definition at line 2269 of file mesh_generation.C.

References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_node(), libMesh::MeshBase::add_node(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Node::build(), libMesh::Elem::build(), libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_remote_elements(), libMesh::Elem::dim(), libMesh::EDGE2, libMesh::EDGE3, libMesh::MeshBase::get_boundary_info(), libMesh::BoundaryInfo::get_edgeset_name_map(), libMesh::BoundaryInfo::get_nodeset_name_map(), libMesh::BoundaryInfo::get_side_boundary_ids(), libMesh::BoundaryInfo::get_sideset_name_map(), libMesh::MeshTools::Generation::QueryElemSubdomainIDBase::get_subdomain_for_layer(), libMesh::HEX27, libMesh::HEX8, libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), TIMPI::Communicator::max(), mesh, libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM18, libMesh::PRISM21, libMesh::PRISM6, libMesh::QUAD4, libMesh::QUAD9, libMesh::MeshBase::query_node_ptr(), libMesh::remote_elem, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMesh::SECOND, libMesh::BoundaryInfo::set_edgeset_name_map(), libMesh::BoundaryInfo::set_nodeset_name_map(), libMesh::BoundaryInfo::set_sideset_name_map(), top_id, libMesh::TRI3, libMesh::TRI6, and libMesh::TRI7.

Referenced by main(), and MeshExtruderTest::testExtruder().

2274 {
2275  LOG_SCOPE("build_extrusion()", "MeshTools::Generation");
2276 
2277  if (!cross_section.n_elem())
2278  return;
2279 
2280  dof_id_type orig_elem = cross_section.n_elem();
2281  dof_id_type orig_nodes = cross_section.n_nodes();
2282 
2283 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2284  unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
2285 #endif
2286 
2287  unsigned int order = 1;
2288 
2289  BoundaryInfo & boundary_info = mesh.get_boundary_info();
2290  const BoundaryInfo & cross_section_boundary_info = cross_section.get_boundary_info();
2291 
2292  // Copy name maps from old to new boundary. We won't copy the whole
2293  // BoundaryInfo because that copies bc ids too, and we need to set
2294  // those more carefully.
2295  boundary_info.set_sideset_name_map() = cross_section_boundary_info.get_sideset_name_map();
2296  boundary_info.set_nodeset_name_map() = cross_section_boundary_info.get_nodeset_name_map();
2297  boundary_info.set_edgeset_name_map() = cross_section_boundary_info.get_edgeset_name_map();
2298 
2299  // If cross_section is distributed, so is its extrusion
2300  if (!cross_section.is_serial())
2301  mesh.delete_remote_elements();
2302 
2303  // We know a priori how many elements we'll need
2304  mesh.reserve_elem(nz*orig_elem);
2305 
2306  // For straightforward meshes we need one or two additional layers per
2307  // element.
2308  if (cross_section.elements_begin() != cross_section.elements_end() &&
2309  (*cross_section.elements_begin())->default_order() == SECOND)
2310  order = 2;
2311  mesh.comm().max(order);
2312 
2313  mesh.reserve_nodes((order*nz+1)*orig_nodes);
2314 
2315  // Container to catch the boundary IDs handed back by the BoundaryInfo object
2316  std::vector<boundary_id_type> ids_to_copy;
2317 
2318  for (const auto & node : cross_section.node_ptr_range())
2319  {
2320  for (unsigned int k=0; k != order*nz+1; ++k)
2321  {
2322  const dof_id_type new_node_id = node->id() + k * orig_nodes;
2323  Node * my_node = mesh.query_node_ptr(new_node_id);
2324  if (!my_node)
2325  {
2326  std::unique_ptr<Node> new_node = Node::build
2327  (*node + (extrusion_vector * k / nz / order),
2328  new_node_id);
2329  new_node->processor_id() = node->processor_id();
2330 
2331 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2332  // Let's give the base of the extruded mesh the same
2333  // unique_ids as the source mesh, in case anyone finds that
2334  // a useful map to preserve.
2335  const unique_id_type uid = (k == 0) ?
2336  node->unique_id() :
2337  orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2338 
2339  new_node->set_unique_id(uid);
2340 #endif
2341 
2342  cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2343  boundary_info.add_node(new_node.get(), ids_to_copy);
2344 
2345  mesh.add_node(std::move(new_node));
2346  }
2347  }
2348  }
2349 
2350  const std::set<boundary_id_type> & side_ids =
2351  cross_section_boundary_info.get_side_boundary_ids();
2352 
2353  boundary_id_type next_side_id = side_ids.empty() ?
2354  0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
2355 
2356  // side_ids may not include ids from remote elements, in which case
2357  // some processors may have underestimated the next_side_id; let's
2358  // fix that.
2359  cross_section.comm().max(next_side_id);
2360 
2361  for (const auto & elem : cross_section.element_ptr_range())
2362  {
2363  const ElemType etype = elem->type();
2364 
2365  // build_extrusion currently only works on coarse meshes
2366  libmesh_assert (!elem->parent());
2367 
2368  for (unsigned int k=0; k != nz; ++k)
2369  {
2370  std::unique_ptr<Elem> new_elem;
2371  switch (etype)
2372  {
2373  case EDGE2:
2374  {
2375  new_elem = Elem::build(QUAD4);
2376  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2377  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2378  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2379  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2380 
2381  if (elem->neighbor_ptr(0) == remote_elem)
2382  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2383  if (elem->neighbor_ptr(1) == remote_elem)
2384  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2385 
2386  break;
2387  }
2388  case EDGE3:
2389  {
2390  new_elem = Elem::build(QUAD9);
2391  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2392  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2393  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2394  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2395  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2396  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2397  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2398  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2399  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2400 
2401  if (elem->neighbor_ptr(0) == remote_elem)
2402  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2403  if (elem->neighbor_ptr(1) == remote_elem)
2404  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2405 
2406  break;
2407  }
2408  case TRI3:
2409  {
2410  new_elem = Elem::build(PRISM6);
2411  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2412  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2413  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2414  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2415  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2416  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2417 
2418  if (elem->neighbor_ptr(0) == remote_elem)
2419  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2420  if (elem->neighbor_ptr(1) == remote_elem)
2421  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2422  if (elem->neighbor_ptr(2) == remote_elem)
2423  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2424 
2425  break;
2426  }
2427  case TRI6:
2428  {
2429  new_elem = Elem::build(PRISM18);
2430  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2431  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2432  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2433  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2434  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2435  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2436  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2437  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2438  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2439  new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2440  new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2441  new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2442  new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2443  new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2444  new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2445  new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2446  new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2447  new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2448 
2449  if (elem->neighbor_ptr(0) == remote_elem)
2450  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2451  if (elem->neighbor_ptr(1) == remote_elem)
2452  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2453  if (elem->neighbor_ptr(2) == remote_elem)
2454  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2455 
2456  break;
2457  }
2458  case TRI7:
2459  {
2460  new_elem = Elem::build(PRISM21);
2461  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2462  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2463  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2464  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2465  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2466  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2467  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2468  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2469  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2470  new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2471  new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2472  new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2473  new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2474  new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2475  new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2476  new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2477  new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2478  new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2479 
2480  new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2481  new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2482  new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2483 
2484  if (elem->neighbor_ptr(0) == remote_elem)
2485  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2486  if (elem->neighbor_ptr(1) == remote_elem)
2487  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2488  if (elem->neighbor_ptr(2) == remote_elem)
2489  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2490 
2491  break;
2492  }
2493  case QUAD4:
2494  {
2495  new_elem = Elem::build(HEX8);
2496  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2497  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2498  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2499  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes)));
2500  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2501  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2502  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2503  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes)));
2504 
2505  if (elem->neighbor_ptr(0) == remote_elem)
2506  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2507  if (elem->neighbor_ptr(1) == remote_elem)
2508  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2509  if (elem->neighbor_ptr(2) == remote_elem)
2510  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2511  if (elem->neighbor_ptr(3) == remote_elem)
2512  new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2513 
2514  break;
2515  }
2516  case QUAD9:
2517  {
2518  new_elem = Elem::build(HEX27);
2519  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2520  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2521  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2522  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2523  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2524  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2525  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2526  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2527  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2528  new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2529  new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2530  new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes)));
2531  new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2532  new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2533  new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2534  new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2535  new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2536  new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2537  new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2538  new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes)));
2539  new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes)));
2540  new_elem->set_node(21, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2541  new_elem->set_node(22, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2542  new_elem->set_node(23, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2543  new_elem->set_node(24, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes)));
2544  new_elem->set_node(25, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes)));
2545  new_elem->set_node(26, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes)));
2546 
2547  if (elem->neighbor_ptr(0) == remote_elem)
2548  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2549  if (elem->neighbor_ptr(1) == remote_elem)
2550  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2551  if (elem->neighbor_ptr(2) == remote_elem)
2552  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2553  if (elem->neighbor_ptr(3) == remote_elem)
2554  new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2555 
2556  break;
2557  }
2558  default:
2559  {
2560  libmesh_not_implemented();
2561  break;
2562  }
2563  }
2564 
2565  new_elem->set_id(elem->id() + (k * orig_elem));
2566  new_elem->processor_id() = elem->processor_id();
2567 
2568 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2569  // Let's give the base of the extruded mesh the same
2570  // unique_ids as the source mesh, in case anyone finds that
2571  // a useful map to preserve.
2572  const unique_id_type uid = (k == 0) ?
2573  elem->unique_id() :
2574  orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + orig_nodes + elem->id();
2575 
2576  new_elem->set_unique_id(uid);
2577 #endif
2578 
2579  if (!elem_subdomain)
2580  // maintain the subdomain_id
2581  new_elem->subdomain_id() = elem->subdomain_id();
2582  else
2583  // Allow the user to choose new subdomain_ids
2584  new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2585 
2586  Elem * added_elem = mesh.add_elem(std::move(new_elem));
2587 
2588  // Copy any old boundary ids on all sides
2589  for (auto s : elem->side_index_range())
2590  {
2591  cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2592 
2593  if (added_elem->dim() == 3)
2594  {
2595  // For 2D->3D extrusion, we give the boundary IDs
2596  // for side s on the old element to side s+1 on the
2597  // new element. This is just a happy coincidence as
2598  // far as I can tell...
2599  boundary_info.add_side(added_elem,
2600  cast_int<unsigned short>(s+1),
2601  ids_to_copy);
2602  }
2603  else
2604  {
2605  // For 1D->2D extrusion, the boundary IDs map as:
2606  // Old elem -> New elem
2607  // 0 -> 3
2608  // 1 -> 1
2609  libmesh_assert_less(s, 2);
2610  const unsigned short sidemap[2] = {3, 1};
2611  boundary_info.add_side(added_elem, sidemap[s], ids_to_copy);
2612  }
2613  }
2614 
2615  // Give new boundary ids to bottom and top
2616  if (k == 0)
2617  boundary_info.add_side(added_elem, 0, next_side_id);
2618  if (k == nz-1)
2619  {
2620  // For 2D->3D extrusion, the "top" ID is 1+the original
2621  // element's number of sides. For 1D->2D extrusion, the
2622  // "top" ID is side 2.
2623  const unsigned short top_id = added_elem->dim() == 3 ?
2624  cast_int<unsigned short>(elem->n_sides()+1) : 2;
2625  boundary_info.add_side
2626  (added_elem, top_id,
2627  cast_int<boundary_id_type>(next_side_id+1));
2628  }
2629  }
2630  }
2631 
2632  // Done building the mesh. Now prepare it for use.
2633  mesh.prepare_for_use();
2634 }
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::boundary_id_type
int8_t boundary_id_type
Definition: id_types.h:51
libMesh::libmesh_assert
libmesh_assert(ctx)
top_id
const boundary_id_type top_id
Definition: periodic_bc_test.C:156
libMesh::unique_id_type
uint8_t unique_id_type
Definition: id_types.h:86
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57

◆ build_line()

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 1673 of file mesh_generation.C.

References build_cube(), and mesh.

Referenced by Biharmonic::Biharmonic(), BoundaryMesh0DTest::build_mesh(), NodalNeighborsTest::do_test(), main(), DualShapeTest::setUp(), SystemsTest::simpleSetup(), SystemsTest::test100KVariables(), MeshSpatialDimensionTest::test1D(), ConstraintOperatorTest::test1DCoarseningNewNodes(), ConstraintOperatorTest::test1DCoarseningOperator(), MeshGenerationTest::testBuildLine(), ConnectedComponentsTest::testEdge(), VolumeTest::testEdge3Volume(), MeshSubdomainIDTest::testMultiple(), BoundaryInfoTest::testNameCopying(), EquationSystemsTest::testPostInitAddElem(), SystemsTest::testProjectLine(), SystemsTest::testProjectMatrix1D(), EquationSystemsTest::testReinitWithNodeElem(), DistributedMeshTest::testRemoteElemError(), and EquationSystemsTest::testSelectivePRefine().

1678 {
1679  // This method only makes sense in 1D!
1680  // But we now just turn a non-1D mesh into a 1D mesh
1681  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1682 
1683  build_cube(mesh,
1684  nx, 0, 0,
1685  xmin, xmax,
1686  0., 0.,
1687  0., 0.,
1688  type,
1689  gauss_lobatto_grid);
1690 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::MeshTools::Generation::build_cube
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)
Builds a (elements) cube.
Definition: mesh_generation.C:313

◆ build_point()

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 1655 of file mesh_generation.C.

References build_cube(), and mesh.

Referenced by TimeSolverTestImplementation< NewmarkSolver >::run_test_with_exact_soln(), EquationSystemsTest::testInit(), EquationSystemsTest::testPostInitAddRealSystem(), and EquationSystemsTest::testPostInitAddSystem().

1658 {
1659  // This method only makes sense in 0D!
1660  // But we now just turn a non-0D mesh into a 0D mesh
1661  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1662 
1663  build_cube(mesh,
1664  0, 0, 0,
1665  0., 0.,
1666  0., 0.,
1667  0., 0.,
1668  type,
1669  gauss_lobatto_grid);
1670 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::MeshTools::Generation::build_cube
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)
Builds a (elements) cube.
Definition: mesh_generation.C:313

◆ build_sphere()

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 1725 of file mesh_generation.C.

Referenced by main(), MeshGenerationTest::testBuildSphere(), InfFERadialTest::testInfQuants_numericDeriv(), MeshTetTest::testSphereShell(), and MeshTriangulationTest::testTriangulatorRoundHole().

1731 {
1732  libmesh_error_msg("Building a circle/sphere only works with AMR.");
1733 }

◆ build_square()

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 1694 of file mesh_generation.C.

References build_cube(), and mesh.

Referenced by AllSecondOrderTest::allCompleteOrder(), AllSecondOrderTest::allCompleteOrderDoNothing(), AllSecondOrderTest::allCompleteOrderMixed(), AllSecondOrderTest::allCompleteOrderMixedFixing(), AllSecondOrderTest::allCompleteOrderRange(), AllSecondOrderTest::allSecondOrder(), AllSecondOrderTest::allSecondOrderDoNothing(), AllSecondOrderTest::allSecondOrderMixed(), AllSecondOrderTest::allSecondOrderMixedFixing(), AllSecondOrderTest::allSecondOrderRange(), Biharmonic::Biharmonic(), BoundaryMeshTest::build_mesh(), main(), WriteVecAndScalar::setupTests(), MultiEvaluablePredTest::test(), MeshSpatialDimensionTest::test2D(), SystemsTest::test2DProjectVectorFE(), DistortTest::test_helper_2D(), AllTriTest::test_helper_2D(), MeshFunctionTest::test_subdomain_id_sets(), DofMapTest::testBadElemFECombo(), EquationSystemsTest::testBadVarNames(), BoundaryInfoTest::testBoundaryOnChildrenBoundaryIDs(), BoundaryInfoTest::testBoundaryOnChildrenBoundarySides(), BoundaryInfoTest::testBoundaryOnChildrenElementsRefineCoarsen(), BoundaryInfoTest::testBoundaryOnChildrenErrors(), MeshGenerationTest::testBuildSquare(), DofMapTest::testConstraintLoopDetection(), MeshAssignTest::testCopyConstruct(), MeshInputTest::testCopyElementSolutionImpl(), MeshInputTest::testCopyElementVectorImpl(), MeshInputTest::testCopyNodalSolutionImpl(), MeshDeletionsTest::testDeleteElem(), MeshInputTest::testExodusReadHeader(), MeshExtruderTest::testExtruder(), MixedOrderTest::testFindNeighbors(), MappedSubdomainPartitionerTest::testMappedSubdomainPartitioner(), BoundaryInfoTest::testMesh(), MeshBaseTest::testMeshBaseVerifyIsPrepared(), MeshAssignTest::testMeshMoveAssign(), MeshInputTest::testNemesisReadImpl(), PointLocatorTest::testPlanar(), SystemsTest::testProjectMatrix2D(), SystemsTest::testProjectSquare(), VolumeTest::testQuad4TrueCentroid(), EquationSystemsTest::testRefineThenReinitPreserveFlags(), BoundaryInfoTest::testRenumber(), EquationSystemsTest::testRepartitionThenReinit(), MeshInputTest::testSingleElementImpl(), MeshSmootherTest::testSmoother(), CheckpointIOTest::testSplitter(), MixedOrderTest::testStitch(), MeshTriangulationTest::testTriangulatorMeshedHoles(), SimplexRefinementTest::testTriRefinement(), MeshInputTest::testVTKPreserveElemIds(), MeshInputTest::testVTKPreserveSubdomainIds(), WriteNodesetData::testWriteImpl(), WriteSidesetData::testWriteImpl(), and WriteElemsetData::testWriteImpl().

1701 {
1702  // This method only makes sense in 2D!
1703  // But we now just turn a non-2D mesh into a 2D mesh
1704  //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1705 
1706  // Call the build_cube() member to actually do the work for us.
1707  build_cube (mesh,
1708  nx, ny, 0,
1709  xmin, xmax,
1710  ymin, ymax,
1711  0., 0.,
1712  type,
1713  gauss_lobatto_grid);
1714 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::MeshTools::Generation::build_cube
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)
Builds a (elements) cube.
Definition: mesh_generation.C:313

◆ surface_octahedron()

void libMesh::MeshTools::Generation::surface_octahedron ( UnstructuredMesh &  mesh,
Real  xmin,
Real  xmax,
Real  ymin,
Real  ymax,
Real  zmin,
Real  zmax,
bool  flip_tris = false 
)

Meshes the surface of an octahedron with 8 Tri3 elements, with counter-clockwise (libMesh default) node ordering as viewed from the octahedron exterior if flip_tris is false or from the interior otherwise.

Definition at line 2756 of file mesh_generation.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::Elem::build(), libMesh::MeshBase::get_boundary_info(), mesh, libMesh::MeshBase::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::Real, libMesh::Elem::set_node(), and libMesh::TRI3.

Referenced by SimplexRefinementTest::test3DTriRefinement().

2761 {
2762  const Real xavg = (xmin + xmax)/2;
2763  const Real yavg = (ymin + ymax)/2;
2764  const Real zavg = (zmin + zmax)/2;
2765  mesh.add_point(Point(xavg,yavg,zmin), 0);
2766  mesh.add_point(Point(xmax,yavg,zavg), 1);
2767  mesh.add_point(Point(xavg,ymax,zavg), 2);
2768  mesh.add_point(Point(xmin,yavg,zavg), 3);
2769  mesh.add_point(Point(xavg,ymin,zavg), 4);
2770  mesh.add_point(Point(xavg,yavg,zmax), 5);
2771 
2772  auto add_tri = [&mesh, flip_tris](std::array<dof_id_type,3> nodes)
2773  {
2774  auto elem = mesh.add_elem(Elem::build(TRI3));
2775  elem->set_node(0, mesh.node_ptr(nodes[0]));
2776  elem->set_node(1, mesh.node_ptr(nodes[1]));
2777  elem->set_node(2, mesh.node_ptr(nodes[2]));
2778  if (flip_tris)
2779  elem->flip(&mesh.get_boundary_info());
2780  };
2781 
2782  add_tri({0,2,1});
2783  add_tri({0,3,2});
2784  add_tri({0,4,3});
2785  add_tri({0,1,4});
2786  add_tri({5,4,1});
2787  add_tri({5,3,4});
2788  add_tri({5,2,3});
2789  add_tri({5,1,2});
2790 
2791  mesh.prepare_for_use();
2792 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
Generated on Wed Jul 23 2025 20:03:49 for libMesh by   doxygen 1.8.14