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

Namespaces

 Private
 

Classes

class  QueryElemSubdomainIDBase
 

Functions

void build_cube (UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 
void build_point (UnstructuredMesh &mesh, const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 
void build_line (UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 
void build_square (UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 
void build_sphere (UnstructuredMesh &mesh, const Real rad=1, const unsigned int nr=2, const ElemType type=INVALID_ELEM, const unsigned int n_smooth=2, const bool flat=true)
 
void build_extrusion (UnstructuredMesh &mesh, const MeshBase &cross_section, const unsigned int nz, RealVectorValue extrusion_vector, QueryElemSubdomainIDBase *elem_subdomain=nullptr)
 
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)
 

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

References libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::MeshTools::Generation::Private::idx(), libMesh::INVALID_ELEM, libMesh::BoundaryInfo::invalid_id, mesh, libMesh::NODEELEM, libMesh::BoundaryInfo::nodeset_name(), 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::DofObject::set_id(), libMesh::Elem::set_node(), side, libMesh::BoundaryInfo::sideset_name(), libMesh::TET10, libMesh::TET4, libMesh::TRI3, and libMesh::TRI6.

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

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

◆ 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 2319 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(), mesh, 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().

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

◆ 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 1992 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::Parallel::Communicator::max(), mesh, 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().

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

References build_cube(), and mesh.

1483 {
1484  // This method only makes sense in 1D!
1485  // But we now just turn a non-1D mesh into a 1D mesh
1486  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1487 
1488  build_cube(mesh,
1489  nx, 0, 0,
1490  xmin, xmax,
1491  0., 0.,
1492  0., 0.,
1493  type,
1494  gauss_lobatto_grid);
1495 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
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)
Definition: mesh_generation.C:296

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

References build_cube(), and mesh.

1463 {
1464  // This method only makes sense in 0D!
1465  // But we now just turn a non-0D mesh into a 0D mesh
1466  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1467 
1468  build_cube(mesh,
1469  0, 0, 0,
1470  0., 0.,
1471  0., 0.,
1472  0., 0.,
1473  type,
1474  gauss_lobatto_grid);
1475 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
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)
Definition: mesh_generation.C:296

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

1536 {
1537  libmesh_error_msg("Building a circle/sphere only works with AMR.");
1538 }

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

References build_cube(), and mesh.

1506 {
1507  // This method only makes sense in 2D!
1508  // But we now just turn a non-2D mesh into a 2D mesh
1509  //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1510 
1511  // Call the build_cube() member to actually do the work for us.
1512  build_cube (mesh,
1513  nx, ny, 0,
1514  xmin, xmax,
1515  ymin, ymax,
1516  0., 0.,
1517  type,
1518  gauss_lobatto_grid);
1519 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1245
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)
Definition: mesh_generation.C:296
generated by doxygen