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=libmesh_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=libmesh_nullptr)
 Meshes a rectangular (2D) region (with or without holes) with a Delaunay triangulation. More...
 

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

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::UnstructuredMesh::all_second_order(), libMesh::BoundaryInfo::boundary_ids(), libMesh::MeshBase::clear(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::MeshBase::element_ptr_range(), libMesh::MeshTools::Generation::Private::GaussLobattoRedistributionFunction::GaussLobattoRedistributionFunction(), libMesh::MeshBase::get_boundary_info(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::MeshTools::Generation::Private::idx(), libMesh::INVALID_ELEM, libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_assert(), libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::MeshBase::node_ref(), libMesh::NODEELEM, libMesh::BoundaryInfo::nodeset_name(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::MeshTools::Modification::redistribute(), libMesh::BoundaryInfo::remove(), libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMesh::DofObject::set_id(), libMesh::MeshBase::set_mesh_dimension(), libMesh::Elem::set_node(), libMesh::MeshBase::set_spatial_dimension(), side, libMesh::BoundaryInfo::sideset_name(), libMesh::TET10, libMesh::TET4, libMesh::TRI3, and libMesh::TRI6.

Referenced by add_cube_convex_hull_to_mesh(), Biharmonic::Biharmonic(), build_line(), build_point(), build_square(), main(), ElemTest< elem_type >::setUp(), ParsedFEMFunctionTest::setUp(), setup(), FETest< order, family, elem_type >::setUp(), AllTriTest::test_helper_3D(), DefaultCouplingTest::testCoupling(), PointNeighborCouplingTest::testCoupling(), DofMapTest::testDofOwner(), BoundaryInfoTest::testEdgeBoundaryConditions(), PointLocatorTest::testLocator(), SystemsTest::testProjectCube(), and SystemsTest::testProjectCubeWithMeshFunction().

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

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

2311 {
2312  // Check for reasonable size
2313  libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2314  libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2315  libmesh_assert_less (xmin, xmax);
2316  libmesh_assert_less (ymin, ymax);
2317 
2318  // Clear out any data which may have been in the Mesh
2319  mesh.clear();
2320 
2321  BoundaryInfo & boundary_info = mesh.get_boundary_info();
2322 
2323  // Make sure the new Mesh will be 2D
2324  mesh.set_mesh_dimension(2);
2325 
2326  // The x and y spacing between boundary points
2327  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2328  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2329 
2330  // Bottom
2331  for (unsigned int p=0; p<=nx; ++p)
2332  mesh.add_point(Point(xmin + p*delta_x, ymin));
2333 
2334  // Right side
2335  for (unsigned int p=1; p<ny; ++p)
2336  mesh.add_point(Point(xmax, ymin + p*delta_y));
2337 
2338  // Top
2339  for (unsigned int p=0; p<=nx; ++p)
2340  mesh.add_point(Point(xmax - p*delta_x, ymax));
2341 
2342  // Left side
2343  for (unsigned int p=1; p<ny; ++p)
2344  mesh.add_point(Point(xmin, ymax - p*delta_y));
2345 
2346  // Be sure we added as many points as we thought we did
2347  libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2348 
2349  // Construct the Triangle Interface object
2350  TriangleInterface t(mesh);
2351 
2352  // Set custom variables for the triangulation
2353  t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2354  t.triangulation_type() = TriangleInterface::PSLG;
2355  t.elem_type() = type;
2356 
2357  if (holes != libmesh_nullptr)
2358  t.attach_hole_list(holes);
2359 
2360  // Triangulate!
2361  t.triangulate();
2362 
2363  // The mesh is now generated, but we still need to mark the boundaries
2364  // to be consistent with the other build_square routines. Note that all
2365  // hole boundary elements get the same ID, 4.
2366  for (auto & elem : mesh.element_ptr_range())
2367  for (auto s : elem->side_index_range())
2368  if (elem->neighbor_ptr(s) == libmesh_nullptr)
2369  {
2370  UniquePtr<const Elem> side (elem->build_side_ptr(s));
2371 
2372  // Check the location of the side's midpoint. Since
2373  // the square has straight sides, the midpoint is not
2374  // on the corner and thus it is uniquely on one of the
2375  // sides.
2376  Point side_midpoint= 0.5f*( side->point(0) + side->point(1) );
2377 
2378  // The boundary ids are set following the same convention as Quad4 sides
2379  // bottom = 0
2380  // right = 1
2381  // top = 2
2382  // left = 3
2383  // hole = 4
2385 
2386  // bottom
2387  if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2388  bc_id=0;
2389 
2390  // right
2391  else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2392  bc_id=1;
2393 
2394  // top
2395  else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2396  bc_id=2;
2397 
2398  // left
2399  else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2400  bc_id=3;
2401 
2402  // If the point is not on any of the external boundaries, it
2403  // is on one of the holes....
2404 
2405  // Finally, add this element's information to the boundary info object.
2406  boundary_info.add_side(elem->id(), s, bc_id);
2407  }
2408 
2409 } // end build_delaunay_square
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
static const Real TOLERANCE
boundary_id_type bc_id
Definition: xdr_io.C:50
int8_t boundary_id_type
Definition: id_types.h:51
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void libMesh::MeshTools::Generation::build_extrusion ( UnstructuredMesh mesh,
const MeshBase cross_section,
const unsigned int  nz,
RealVectorValue  extrusion_vector,
QueryElemSubdomainIDBase elem_subdomain = libmesh_nullptr 
)

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

Definition at line 1977 of file mesh_generation.C.

References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_node(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_remote_elements(), libMesh::Elem::dim(), libMesh::EDGE2, libMesh::EDGE3, libMesh::MeshBase::element_ptr_range(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::MeshBase::get_boundary_info(), libMesh::BoundaryInfo::get_side_boundary_ids(), libMesh::MeshTools::Generation::QueryElemSubdomainIDBase::get_subdomain_for_layer(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libMesh::Parallel::Communicator::max(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::MeshBase::node_ptr_range(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD9, libMesh::remote_elem, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMesh::SECOND, libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::DofObject::set_unique_id(), libMesh::Elem::subdomain_id(), libMesh::TRI3, libMesh::TRI6, and libMesh::DofObject::unique_id().

Referenced by MeshExtruderTest::testExtruder().

1982 {
1983  if (!cross_section.n_elem())
1984  return;
1985 
1986  START_LOG("build_extrusion()", "MeshTools::Generation");
1987 
1988  dof_id_type orig_elem = cross_section.n_elem();
1989  dof_id_type orig_nodes = cross_section.n_nodes();
1990 
1991 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1992  unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
1993 #endif
1994 
1995  unsigned int order = 1;
1996 
1997  BoundaryInfo & boundary_info = mesh.get_boundary_info();
1998  const BoundaryInfo & cross_section_boundary_info = cross_section.get_boundary_info();
1999 
2000  // If cross_section is distributed, so is its extrusion
2001  if (!cross_section.is_serial())
2002  mesh.delete_remote_elements();
2003 
2004  // We know a priori how many elements we'll need
2005  mesh.reserve_elem(nz*orig_elem);
2006 
2007  // For straightforward meshes we need one or two additional layers per
2008  // element.
2009  if (cross_section.elements_begin() != cross_section.elements_end() &&
2010  (*cross_section.elements_begin())->default_order() == SECOND)
2011  order = 2;
2012  mesh.comm().max(order);
2013 
2014  mesh.reserve_nodes((order*nz+1)*orig_nodes);
2015 
2016  // Container to catch the boundary IDs handed back by the BoundaryInfo object
2017  std::vector<boundary_id_type> ids_to_copy;
2018 
2019  for (const auto & node : cross_section.node_ptr_range())
2020  {
2021  for (unsigned int k=0; k != order*nz+1; ++k)
2022  {
2023  Node * new_node =
2024  mesh.add_point(*node +
2025  (extrusion_vector * k / nz / order),
2026  node->id() + (k * orig_nodes),
2027  node->processor_id());
2028 
2029 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2030  // Let's give the base of the extruded mesh the same
2031  // unique_ids as the source mesh, in case anyone finds that
2032  // a useful map to preserve.
2033  const unique_id_type uid = (k == 0) ?
2034  node->unique_id() :
2035  orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2036 
2037  new_node->set_unique_id() = uid;
2038 #endif
2039 
2040  cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2041  boundary_info.add_node(new_node, ids_to_copy);
2042  }
2043  }
2044 
2045  const std::set<boundary_id_type> & side_ids =
2046  cross_section_boundary_info.get_side_boundary_ids();
2047 
2048  boundary_id_type next_side_id = side_ids.empty() ?
2049  0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
2050 
2051  // side_ids may not include ids from remote elements, in which case
2052  // some processors may have underestimated the next_side_id; let's
2053  // fix that.
2054  cross_section.comm().max(next_side_id);
2055 
2056  for (const auto & elem : cross_section.element_ptr_range())
2057  {
2058  const ElemType etype = elem->type();
2059 
2060  // build_extrusion currently only works on coarse meshes
2061  libmesh_assert (!elem->parent());
2062 
2063  for (unsigned int k=0; k != nz; ++k)
2064  {
2065  Elem * new_elem;
2066  switch (etype)
2067  {
2068  case EDGE2:
2069  {
2070  new_elem = new Quad4;
2071  new_elem->set_node(0) = mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes));
2072  new_elem->set_node(1) = mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes));
2073  new_elem->set_node(2) = mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes));
2074  new_elem->set_node(3) = mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes));
2075 
2076  if (elem->neighbor_ptr(0) == remote_elem)
2077  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2078  if (elem->neighbor_ptr(1) == remote_elem)
2079  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2080 
2081  break;
2082  }
2083  case EDGE3:
2084  {
2085  new_elem = new Quad9;
2086  new_elem->set_node(0) = mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes));
2087  new_elem->set_node(1) = mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes));
2088  new_elem->set_node(2) = mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes));
2089  new_elem->set_node(3) = mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes));
2090  new_elem->set_node(4) = mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes));
2091  new_elem->set_node(5) = mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes));
2092  new_elem->set_node(6) = mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes));
2093  new_elem->set_node(7) = mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes));
2094  new_elem->set_node(8) = mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes));
2095 
2096  if (elem->neighbor_ptr(0) == remote_elem)
2097  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2098  if (elem->neighbor_ptr(1) == remote_elem)
2099  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2100 
2101  break;
2102  }
2103  case TRI3:
2104  {
2105  new_elem = new Prism6;
2106  new_elem->set_node(0) = mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes));
2107  new_elem->set_node(1) = mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes));
2108  new_elem->set_node(2) = mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes));
2109  new_elem->set_node(3) = mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes));
2110  new_elem->set_node(4) = mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes));
2111  new_elem->set_node(5) = mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes));
2112 
2113  if (elem->neighbor_ptr(0) == remote_elem)
2114  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2115  if (elem->neighbor_ptr(1) == remote_elem)
2116  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2117  if (elem->neighbor_ptr(2) == remote_elem)
2118  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2119 
2120  break;
2121  }
2122  case TRI6:
2123  {
2124  new_elem = new Prism18;
2125  new_elem->set_node(0) = mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes));
2126  new_elem->set_node(1) = mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes));
2127  new_elem->set_node(2) = mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes));
2128  new_elem->set_node(3) = mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes));
2129  new_elem->set_node(4) = mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes));
2130  new_elem->set_node(5) = mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes));
2131  new_elem->set_node(6) = mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes));
2132  new_elem->set_node(7) = mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes));
2133  new_elem->set_node(8) = mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes));
2134  new_elem->set_node(9) = mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes));
2135  new_elem->set_node(10) = mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes));
2136  new_elem->set_node(11) = mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes));
2137  new_elem->set_node(12) = mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes));
2138  new_elem->set_node(13) = mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes));
2139  new_elem->set_node(14) = mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes));
2140  new_elem->set_node(15) = mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes));
2141  new_elem->set_node(16) = mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes));
2142  new_elem->set_node(17) = mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes));
2143 
2144  if (elem->neighbor_ptr(0) == remote_elem)
2145  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2146  if (elem->neighbor_ptr(1) == remote_elem)
2147  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2148  if (elem->neighbor_ptr(2) == remote_elem)
2149  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2150 
2151  break;
2152  }
2153  case QUAD4:
2154  {
2155  new_elem = new Hex8;
2156  new_elem->set_node(0) = mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes));
2157  new_elem->set_node(1) = mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes));
2158  new_elem->set_node(2) = mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes));
2159  new_elem->set_node(3) = mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes));
2160  new_elem->set_node(4) = mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes));
2161  new_elem->set_node(5) = mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes));
2162  new_elem->set_node(6) = mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes));
2163  new_elem->set_node(7) = mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes));
2164 
2165  if (elem->neighbor_ptr(0) == remote_elem)
2166  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2167  if (elem->neighbor_ptr(1) == remote_elem)
2168  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2169  if (elem->neighbor_ptr(2) == remote_elem)
2170  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2171  if (elem->neighbor_ptr(3) == remote_elem)
2172  new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2173 
2174  break;
2175  }
2176  case QUAD9:
2177  {
2178  new_elem = new Hex27;
2179  new_elem->set_node(0) = mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes));
2180  new_elem->set_node(1) = mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes));
2181  new_elem->set_node(2) = mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes));
2182  new_elem->set_node(3) = mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes));
2183  new_elem->set_node(4) = mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes));
2184  new_elem->set_node(5) = mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes));
2185  new_elem->set_node(6) = mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes));
2186  new_elem->set_node(7) = mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes));
2187  new_elem->set_node(8) = mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes));
2188  new_elem->set_node(9) = mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes));
2189  new_elem->set_node(10) = mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes));
2190  new_elem->set_node(11) = mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes));
2191  new_elem->set_node(12) = mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes));
2192  new_elem->set_node(13) = mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes));
2193  new_elem->set_node(14) = mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes));
2194  new_elem->set_node(15) = mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes));
2195  new_elem->set_node(16) = mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes));
2196  new_elem->set_node(17) = mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes));
2197  new_elem->set_node(18) = mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes));
2198  new_elem->set_node(19) = mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes));
2199  new_elem->set_node(20) = mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes));
2200  new_elem->set_node(21) = mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes));
2201  new_elem->set_node(22) = mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes));
2202  new_elem->set_node(23) = mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes));
2203  new_elem->set_node(24) = mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes));
2204  new_elem->set_node(25) = mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes));
2205  new_elem->set_node(26) = mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes));
2206 
2207  if (elem->neighbor_ptr(0) == remote_elem)
2208  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2209  if (elem->neighbor_ptr(1) == remote_elem)
2210  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2211  if (elem->neighbor_ptr(2) == remote_elem)
2212  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2213  if (elem->neighbor_ptr(3) == remote_elem)
2214  new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2215 
2216  break;
2217  }
2218  default:
2219  {
2220  libmesh_not_implemented();
2221  break;
2222  }
2223  }
2224 
2225  new_elem->set_id(elem->id() + (k * orig_elem));
2226  new_elem->processor_id() = elem->processor_id();
2227 
2228 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2229  // Let's give the base of the extruded mesh the same
2230  // unique_ids as the source mesh, in case anyone finds that
2231  // a useful map to preserve.
2232  const unique_id_type uid = (k == 0) ?
2233  elem->unique_id() :
2234  orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + orig_nodes + elem->id();
2235 
2236  new_elem->set_unique_id() = uid;
2237 #endif
2238 
2239  if (!elem_subdomain)
2240  // maintain the subdomain_id
2241  new_elem->subdomain_id() = elem->subdomain_id();
2242  else
2243  // Allow the user to choose new subdomain_ids
2244  new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2245 
2246  new_elem = mesh.add_elem(new_elem);
2247 
2248  // Copy any old boundary ids on all sides
2249  for (auto s : elem->side_index_range())
2250  {
2251  cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2252 
2253  if (new_elem->dim() == 3)
2254  {
2255  // For 2D->3D extrusion, we give the boundary IDs
2256  // for side s on the old element to side s+1 on the
2257  // new element. This is just a happy coincidence as
2258  // far as I can tell...
2259  boundary_info.add_side(new_elem,
2260  cast_int<unsigned short>(s+1),
2261  ids_to_copy);
2262  }
2263  else
2264  {
2265  // For 1D->2D extrusion, the boundary IDs map as:
2266  // Old elem -> New elem
2267  // 0 -> 3
2268  // 1 -> 1
2269  libmesh_assert_less(s, 2);
2270  const unsigned short sidemap[2] = {3, 1};
2271  boundary_info.add_side(new_elem, sidemap[s], ids_to_copy);
2272  }
2273  }
2274 
2275  // Give new boundary ids to bottom and top
2276  if (k == 0)
2277  boundary_info.add_side(new_elem, 0, next_side_id);
2278  if (k == nz-1)
2279  {
2280  // For 2D->3D extrusion, the "top" ID is 1+the original
2281  // element's number of sides. For 1D->2D extrusion, the
2282  // "top" ID is side 2.
2283  const unsigned short top_id = new_elem->dim() == 3 ?
2284  cast_int<unsigned short>(elem->n_sides()+1) : 2;
2285  boundary_info.add_side
2286  (new_elem, top_id,
2287  cast_int<boundary_id_type>(next_side_id+1));
2288  }
2289  }
2290  }
2291 
2292  STOP_LOG("build_extrusion()", "MeshTools::Generation");
2293 
2294  // Done building the mesh. Now prepare it for use.
2295  mesh.prepare_for_use(/*skip_renumber =*/ false);
2296 }
ElemType
Defines an enum for geometric element types.
MeshBase & mesh
libmesh_assert(j)
int8_t boundary_id_type
Definition: id_types.h:51
uint8_t unique_id_type
Definition: id_types.h:79
uint8_t dof_id_type
Definition: id_types.h:64
const RemoteElem * remote_elem
Definition: remote_elem.C:57
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 1464 of file mesh_generation.C.

References build_cube().

Referenced by Biharmonic::Biharmonic(), build_sphere(), NodalNeighborsTest::do_test(), main(), MeshSpatialDimensionTest::test1D(), EquationSystemsTest::testPostInitAddElem(), and SystemsTest::testProjectLine().

1469 {
1470  // This method only makes sense in 1D!
1471  // But we now just turn a non-1D mesh into a 1D mesh
1472  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1473 
1474  build_cube(mesh,
1475  nx, 0, 0,
1476  xmin, xmax,
1477  0., 0.,
1478  0., 0.,
1479  type,
1480  gauss_lobatto_grid);
1481 }
MeshBase & mesh
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.
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 1446 of file mesh_generation.C.

References build_cube().

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

1449 {
1450  // This method only makes sense in 0D!
1451  // But we now just turn a non-0D mesh into a 0D mesh
1452  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1453 
1454  build_cube(mesh,
1455  0, 0, 0,
1456  0., 0.,
1457  0., 0.,
1458  0., 0.,
1459  type,
1460  gauss_lobatto_grid);
1461 }
MeshBase & mesh
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.
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 
)
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 1485 of file mesh_generation.C.

References build_cube().

Referenced by Biharmonic::Biharmonic(), BoundaryMeshTest::build_mesh(), main(), MeshSpatialDimensionTest::test2D(), AllTriTest::test_helper_2D(), MeshExtruderTest::testExtruder(), MappedSubdomainPartitionerTest::testMappedSubdomainPartitioner(), BoundaryInfoTest::testMesh(), SystemsTest::testProjectSquare(), EquationSystemsTest::testRefineThenReinitPreserveFlags(), CheckpointIOTest::testSplitter(), and WriteVecAndScalar::testWrite().

1492 {
1493  // This method only makes sense in 2D!
1494  // But we now just turn a non-2D mesh into a 2D mesh
1495  //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1496 
1497  // Call the build_cube() member to actually do the work for us.
1498  build_cube (mesh,
1499  nx, ny, 0,
1500  xmin, xmax,
1501  ymin, ymax,
1502  0., 0.,
1503  type,
1504  gauss_lobatto_grid);
1505 }
MeshBase & mesh
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.