libMesh
Loading...
Searching...
No Matches
Namespaces | Classes | Functions
libMesh::MeshTools::Generation Namespace Reference

Tools for Mesh generation. More...

Namespaces

namespace  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.
 
void build_point (UnstructuredMesh &mesh, const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 A specialized build_cube() for 0D meshes.
 
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.
 
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.
 
void build_sphere (UnstructuredMesh &mesh, const Real radius=1, const unsigned int n_refinements=2, const ElemType type=INVALID_ELEM, const unsigned int n_smooth=2, const bool flat=true)
 Fills mesh with a mesh discretizing a ball (||x||<= radius) or sphere (||x|| = radius) domain.
 
void build_extrusion (UnstructuredMesh &mesh, const MeshBase &cross_section, const unsigned int nz, RealVectorValue extrusion_vector, QueryElemSubdomainIDBase *elem_subdomain=nullptr)
 Meshes the tensor product of a 1D and a 1D-or-2D domain.
 
void build_delaunay_square (UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin, const Real xmax, const Real ymin, const Real ymax, const ElemType type, const std::vector< TriangleInterface::Hole * > *holes=nullptr)
 Meshes a rectangular (2D) region (with or without holes) with a Delaunay triangulation.
 
void surface_octahedron (UnstructuredMesh &mesh, Real xmin, Real xmax, Real ymin, Real ymax, Real zmin, Real zmax, bool flip_tris=false)
 Meshes the surface of an octahedron with 8 Tri3 elements, with counter-clockwise (libMesh default) node ordering as viewed from the octahedron exterior if flip_tris is false or from the interior otherwise.
 

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

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

References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_node(), libMesh::MeshBase::add_node(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::MeshBase::all_complete_order(), libMesh::MeshBase::all_second_order(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::Elem::build_with_id(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::MeshBase::clear(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::Utility::enum_to_string(), libMesh::MeshBase::get_boundary_info(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::DofObject::id(), libMesh::index_range(), libMesh::INVALID_ELEM, libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_assert(), libMesh::make_range(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::node_ptr(), libMesh::NODEELEM, libMesh::BoundaryInfo::nodeset_name(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL4, libMesh::QUADSHELL8, libMesh::QUADSHELL9, libMesh::Real, libMesh::MeshTools::Modification::redistribute(), libMesh::BoundaryInfo::remove(), libMesh::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(), libMesh::BoundaryInfo::sideset_name(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, and libMesh::TRISHELL3.

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

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

2888{
2889 // Check for reasonable size
2890 libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2891 libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2892 libmesh_assert_less (xmin, xmax);
2893 libmesh_assert_less (ymin, ymax);
2894
2895 // Clear out any data which may have been in the Mesh
2896 mesh.clear();
2897
2898 BoundaryInfo & boundary_info = mesh.get_boundary_info();
2899
2900 // Make sure the new Mesh will be 2D
2901 mesh.set_mesh_dimension(2);
2902
2903 // The x and y spacing between boundary points
2904 const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2905 const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2906
2907 // Bottom
2908 for (unsigned int p=0; p<=nx; ++p)
2909 mesh.add_point(Point(xmin + p*delta_x, ymin));
2910
2911 // Right side
2912 for (unsigned int p=1; p<ny; ++p)
2913 mesh.add_point(Point(xmax, ymin + p*delta_y));
2914
2915 // Top
2916 for (unsigned int p=0; p<=nx; ++p)
2917 mesh.add_point(Point(xmax - p*delta_x, ymax));
2918
2919 // Left side
2920 for (unsigned int p=1; p<ny; ++p)
2921 mesh.add_point(Point(xmin, ymax - p*delta_y));
2922
2923 // Be sure we added as many points as we thought we did
2924 libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2925
2926 // Construct the Triangle Interface object
2927 TriangleInterface t(mesh);
2928
2929 // Set custom variables for the triangulation
2930 t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2931 t.triangulation_type() = TriangleInterface::PSLG;
2932 t.elem_type() = type;
2933
2934 if (holes != nullptr)
2935 t.attach_hole_list(holes);
2936
2937 // Triangulate!
2938 t.triangulate();
2939
2940 // For avoiding extraneous side element construction
2941 std::unique_ptr<const Elem> side;
2942
2943 // The mesh is now generated, but we still need to mark the boundaries
2944 // to be consistent with the other build_square routines. Note that all
2945 // hole boundary elements get the same ID, 4.
2946 for (auto & elem : mesh.element_ptr_range())
2947 for (auto s : elem->side_index_range())
2948 if (elem->neighbor_ptr(s) == nullptr)
2949 {
2950 elem->build_side_ptr(side, s);
2951
2952 // Check the location of the side's midpoint. Since
2953 // the square has straight sides, the midpoint is not
2954 // on the corner and thus it is uniquely on one of the
2955 // sides.
2956 Point side_midpoint= 0.5f*( side->point(0) + side->point(1) );
2957
2958 // The boundary ids are set following the same convention as Quad4 sides
2959 // bottom = 0
2960 // right = 1
2961 // top = 2
2962 // left = 3
2963 // hole = 4
2964 boundary_id_type bc_id=4;
2965
2966 // bottom
2967 if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2968 bc_id=0;
2969
2970 // right
2971 else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2972 bc_id=1;
2973
2974 // top
2975 else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2976 bc_id=2;
2977
2978 // left
2979 else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2980 bc_id=3;
2981
2982 // If the point is not on any of the external boundaries, it
2983 // is on one of the holes....
2984
2985 // Finally, add this element's information to the boundary info object.
2986 boundary_info.add_side(elem->id(), s, bc_id);
2987 }
2988
2989} // end build_delaunay_square
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::TriangleInterface
A C++ interface between LibMesh and the Triangle library written by J.R.
Definition mesh_triangle_interface.h:42

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

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

2513{
2514 LOG_SCOPE("build_extrusion()", "MeshTools::Generation");
2515
2516 if (!cross_section.n_elem())
2517 return;
2518
2519 dof_id_type orig_elem = cross_section.n_elem();
2520 dof_id_type orig_nodes = cross_section.n_nodes();
2521
2522#ifdef LIBMESH_ENABLE_UNIQUE_ID
2523 unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
2524#endif
2525
2526 unsigned int order = 1;
2527
2528 BoundaryInfo & boundary_info = mesh.get_boundary_info();
2529 const BoundaryInfo & cross_section_boundary_info = cross_section.get_boundary_info();
2530
2531 // Copy name maps from old to new boundary. We won't copy the whole
2532 // BoundaryInfo because that copies bc ids too, and we need to set
2533 // those more carefully.
2534 boundary_info.set_sideset_name_map() = cross_section_boundary_info.get_sideset_name_map();
2535 boundary_info.set_nodeset_name_map() = cross_section_boundary_info.get_nodeset_name_map();
2536 boundary_info.set_edgeset_name_map() = cross_section_boundary_info.get_edgeset_name_map();
2537
2538 // If cross_section is distributed, so is its extrusion
2539 if (!cross_section.is_serial())
2540 mesh.delete_remote_elements();
2541
2542 // We know a priori how many elements we'll need
2543 mesh.reserve_elem(nz*orig_elem);
2544
2545 // For straightforward meshes we need one or two additional layers per
2546 // element.
2547 if (cross_section.elements_begin() != cross_section.elements_end() &&
2548 (*cross_section.elements_begin())->default_order() == SECOND)
2549 order = 2;
2550 mesh.comm().max(order);
2551
2552 mesh.reserve_nodes((order*nz+1)*orig_nodes);
2553
2554 // Container to catch the boundary IDs handed back by the BoundaryInfo object
2555 std::vector<boundary_id_type> ids_to_copy;
2556
2557 for (const auto & node : cross_section.node_ptr_range())
2558 {
2559 for (unsigned int k=0; k != order*nz+1; ++k)
2560 {
2561 const dof_id_type new_node_id = node->id() + k * orig_nodes;
2562 Node * my_node = mesh.query_node_ptr(new_node_id);
2563 if (!my_node)
2564 {
2565 std::unique_ptr<Node> new_node = Node::build
2566 (*node + (extrusion_vector * k / nz / order),
2567 new_node_id);
2568 new_node->processor_id() = node->processor_id();
2569
2570#ifdef LIBMESH_ENABLE_UNIQUE_ID
2571 // Let's give the base of the extruded mesh the same
2572 // unique_ids as the source mesh, in case anyone finds that
2573 // a useful map to preserve.
2574 const unique_id_type uid = (k == 0) ?
2575 node->unique_id() :
2576 orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2577
2578 new_node->set_unique_id(uid);
2579#endif
2580
2581 cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2582 boundary_info.add_node(new_node.get(), ids_to_copy);
2583
2584 mesh.add_node(std::move(new_node));
2585 }
2586 }
2587 }
2588
2589 const std::set<boundary_id_type> & side_ids =
2590 cross_section_boundary_info.get_side_boundary_ids();
2591
2592 boundary_id_type next_side_id = side_ids.empty() ?
2593 0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
2594
2595 // side_ids may not include ids from remote elements, in which case
2596 // some processors may have underestimated the next_side_id; let's
2597 // fix that.
2598 cross_section.comm().max(next_side_id);
2599
2600 for (const auto & elem : cross_section.element_ptr_range())
2601 {
2602 const ElemType etype = elem->type();
2603
2604 // build_extrusion currently only works on coarse meshes
2605 libmesh_assert (!elem->parent());
2606
2607 for (unsigned int k=0; k != nz; ++k)
2608 {
2609 std::unique_ptr<Elem> new_elem;
2610 switch (etype)
2611 {
2612 case EDGE2:
2613 {
2614 new_elem = Elem::build(QUAD4);
2615 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2616 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2617 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2618 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2619
2620 if (elem->neighbor_ptr(0) == remote_elem)
2621 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2622 if (elem->neighbor_ptr(1) == remote_elem)
2623 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2624
2625 break;
2626 }
2627 case EDGE3:
2628 {
2629 new_elem = Elem::build(QUAD9);
2630 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2631 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2632 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2633 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2634 new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2635 new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2636 new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2637 new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2638 new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2639
2640 if (elem->neighbor_ptr(0) == remote_elem)
2641 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2642 if (elem->neighbor_ptr(1) == remote_elem)
2643 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2644
2645 break;
2646 }
2647 case TRI3:
2648 {
2649 new_elem = Elem::build(PRISM6);
2650 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2651 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2652 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2653 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2654 new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2655 new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2656
2657 if (elem->neighbor_ptr(0) == remote_elem)
2658 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2659 if (elem->neighbor_ptr(1) == remote_elem)
2660 new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2661 if (elem->neighbor_ptr(2) == remote_elem)
2662 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2663
2664 break;
2665 }
2666 case TRI6:
2667 {
2668 new_elem = Elem::build(PRISM18);
2669 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2670 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2671 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2672 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2673 new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2674 new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2675 new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2676 new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2677 new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2678 new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2679 new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2680 new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2681 new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2682 new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2683 new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2684 new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2685 new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2686 new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2687
2688 if (elem->neighbor_ptr(0) == remote_elem)
2689 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2690 if (elem->neighbor_ptr(1) == remote_elem)
2691 new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2692 if (elem->neighbor_ptr(2) == remote_elem)
2693 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2694
2695 break;
2696 }
2697 case TRI7:
2698 {
2699 new_elem = Elem::build(PRISM21);
2700 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2701 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2702 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2703 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2704 new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2705 new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2706 new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2707 new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2708 new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2709 new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2710 new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2711 new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2712 new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2713 new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2714 new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2715 new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2716 new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2717 new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2718
2719 new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2720 new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2721 new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2722
2723 if (elem->neighbor_ptr(0) == remote_elem)
2724 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2725 if (elem->neighbor_ptr(1) == remote_elem)
2726 new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2727 if (elem->neighbor_ptr(2) == remote_elem)
2728 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2729
2730 break;
2731 }
2732 case QUAD4:
2733 {
2734 new_elem = Elem::build(HEX8);
2735 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2736 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2737 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2738 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes)));
2739 new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2740 new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2741 new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2742 new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes)));
2743
2744 if (elem->neighbor_ptr(0) == remote_elem)
2745 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2746 if (elem->neighbor_ptr(1) == remote_elem)
2747 new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2748 if (elem->neighbor_ptr(2) == remote_elem)
2749 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2750 if (elem->neighbor_ptr(3) == remote_elem)
2751 new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2752
2753 break;
2754 }
2755 case QUAD9:
2756 {
2757 new_elem = Elem::build(HEX27);
2758 new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2759 new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2760 new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2761 new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2762 new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2763 new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2764 new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2765 new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2766 new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2767 new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2768 new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2769 new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes)));
2770 new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2771 new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2772 new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2773 new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2774 new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2775 new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2776 new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2777 new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes)));
2778 new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes)));
2779 new_elem->set_node(21, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2780 new_elem->set_node(22, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2781 new_elem->set_node(23, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2782 new_elem->set_node(24, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes)));
2783 new_elem->set_node(25, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes)));
2784 new_elem->set_node(26, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes)));
2785
2786 if (elem->neighbor_ptr(0) == remote_elem)
2787 new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2788 if (elem->neighbor_ptr(1) == remote_elem)
2789 new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2790 if (elem->neighbor_ptr(2) == remote_elem)
2791 new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2792 if (elem->neighbor_ptr(3) == remote_elem)
2793 new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2794
2795 break;
2796 }
2797 default:
2798 {
2799 libmesh_not_implemented();
2800 break;
2801 }
2802 }
2803
2804 new_elem->set_id(elem->id() + (k * orig_elem));
2805 new_elem->processor_id() = elem->processor_id();
2806
2807#ifdef LIBMESH_ENABLE_UNIQUE_ID
2808 // Let's give the base of the extruded mesh the same
2809 // unique_ids as the source mesh, in case anyone finds that
2810 // a useful map to preserve.
2811 const unique_id_type uid = (k == 0) ?
2812 elem->unique_id() :
2813 orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + orig_nodes + elem->id();
2814
2815 new_elem->set_unique_id(uid);
2816#endif
2817
2818 if (!elem_subdomain)
2819 // maintain the subdomain_id
2820 new_elem->subdomain_id() = elem->subdomain_id();
2821 else
2822 // Allow the user to choose new subdomain_ids
2823 new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2824
2825 Elem * added_elem = mesh.add_elem(std::move(new_elem));
2826
2827 // Copy any old boundary ids on all sides
2828 for (auto s : elem->side_index_range())
2829 {
2830 cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2831
2832 if (added_elem->dim() == 3)
2833 {
2834 // For 2D->3D extrusion, we give the boundary IDs
2835 // for side s on the old element to side s+1 on the
2836 // new element. This is just a happy coincidence as
2837 // far as I can tell...
2838 boundary_info.add_side(added_elem,
2839 cast_int<unsigned short>(s+1),
2840 ids_to_copy);
2841 }
2842 else
2843 {
2844 // For 1D->2D extrusion, the boundary IDs map as:
2845 // Old elem -> New elem
2846 // 0 -> 3
2847 // 1 -> 1
2848 libmesh_assert_less(s, 2);
2849 const unsigned short sidemap[2] = {3, 1};
2850 boundary_info.add_side(added_elem, sidemap[s], ids_to_copy);
2851 }
2852 }
2853
2854 // Give new boundary ids to bottom and top
2855 if (k == 0)
2856 boundary_info.add_side(added_elem, 0, next_side_id);
2857 if (k == nz-1)
2858 {
2859 // For 2D->3D extrusion, the "top" ID is 1+the original
2860 // element's number of sides. For 1D->2D extrusion, the
2861 // "top" ID is side 2.
2862 const unsigned short top_id = added_elem->dim() == 3 ?
2863 cast_int<unsigned short>(elem->n_sides()+1) : 2;
2864 boundary_info.add_side
2865 (added_elem, top_id,
2866 cast_int<boundary_id_type>(next_side_id+1));
2867 }
2868 }
2869 }
2870
2871 // Done building the mesh. Now prepare it for use.
2872 mesh.prepare_for_use();
2873}
libMesh::Parallel::Communicator::max
void max(const T &r, T &o, Request &req) const
libMesh::BoundaryInfo::set_sideset_name_map
std::map< boundary_id_type, std::string > & set_sideset_name_map()
Definition boundary_info.h:897
libMesh::BoundaryInfo::get_side_boundary_ids
const std::set< boundary_id_type > & get_side_boundary_ids() const
Definition boundary_info.h:817
libMesh::BoundaryInfo::set_edgeset_name_map
std::map< boundary_id_type, std::string > & set_edgeset_name_map()
Definition boundary_info.h:913
libMesh::BoundaryInfo::get_edgeset_name_map
const std::map< boundary_id_type, std::string > & get_edgeset_name_map() const
Definition boundary_info.h:915
libMesh::BoundaryInfo::get_nodeset_name_map
const std::map< boundary_id_type, std::string > & get_nodeset_name_map() const
Definition boundary_info.h:907
libMesh::BoundaryInfo::set_nodeset_name_map
std::map< boundary_id_type, std::string > & set_nodeset_name_map()
Definition boundary_info.h:905
libMesh::BoundaryInfo::get_sideset_name_map
const std::map< boundary_id_type, std::string > & get_sideset_name_map() const
Definition boundary_info.h:899
libMesh::MeshBase::is_serial
virtual bool is_serial() const
Definition mesh_base.h:357
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition mesh_base.h:170
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::MeshBase::parallel_max_unique_id
virtual unique_id_type parallel_max_unique_id() const =0
libMesh::MeshTools::Generation::QueryElemSubdomainIDBase::get_subdomain_for_layer
virtual subdomain_id_type get_subdomain_for_layer(const Elem *old_elem, unsigned int layer)=0
libMesh::Node
A Node is like a Point, but with more information.
Definition node.h:55
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition parallel_object.h:97
libMesh::unique_id_type
uint8_t unique_id_type
Definition id_types.h:86
top_id
const boundary_id_type top_id
Definition periodic_bc_test.C:156

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

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

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

1965{
1966 // This method only makes sense in 1D!
1967 // But we now just turn a non-1D mesh into a 1D mesh
1968 //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1969
1970 build_cube(mesh,
1971 nx, 0, 0,
1972 xmin, xmax,
1973 0., 0.,
1974 0., 0.,
1975 type,
1976 gauss_lobatto_grid);
1977}
libMesh::MeshTools::Generation::build_cube
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
Definition mesh_generation.C:317

References build_cube(), and mesh.

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

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

1945{
1946 // This method only makes sense in 0D!
1947 // But we now just turn a non-0D mesh into a 0D mesh
1948 //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1949
1950 build_cube(mesh,
1951 0, 0, 0,
1952 0., 0.,
1953 0., 0.,
1954 0., 0.,
1955 type,
1956 gauss_lobatto_grid);
1957}

References build_cube(), and mesh.

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

◆ build_sphere()

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

Fills mesh with a mesh discretizing a ball (||x||<= radius) or sphere (||x|| = radius) domain.

A coarse grid of product elements (edge, quad, or hex) is refined, n_refinements times, snapping boundary nodes to the sphere surface each time; interior nodes are then smoothed n_smooth times.

If the selected element type is 1D, our 1D ball is a line.

If type is 2D, then if flat is true (the default), a 2D ball (a filled circle) is meshed.

In 2D, if flat is false, a 2D sphere (the surface of a 3D ball) is meshed, using TRI3 elements in an icosahedron for the initial mesh.

Definition at line 2012 of file mesh_generation.C.

2018{
2019 libmesh_error_msg("Building a circle/sphere only works with AMR.");
2020}

Referenced by main(), AllRBBTest::test_circle(), AllRBBTest::test_cylinder(), AllRBBTest::test_disk(), AllRBBTest::test_sphere(), MeshGenerationTest::testBuildSphere(), InfFERadialTest::testInfQuants_numericDeriv(), MeshTetTest::testSphereShell(), and MeshTriangulationTest::testTriangulatorRoundHole().

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

1988{
1989 // This method only makes sense in 2D!
1990 // But we now just turn a non-2D mesh into a 2D mesh
1991 //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1992
1993 // Call the build_cube() member to actually do the work for us.
1994 build_cube (mesh,
1995 nx, ny, 0,
1996 xmin, xmax,
1997 ymin, ymax,
1998 0., 0.,
1999 type,
2000 gauss_lobatto_grid);
2001}

References build_cube(), and mesh.

Referenced by AllSecondOrderTest::allCompleteOrder(), AllSecondOrderTest::allCompleteOrderDoNothing(), AllSecondOrderTest::allCompleteOrderMixed(), AllSecondOrderTest::allCompleteOrderMixedFixing(), AllSecondOrderTest::allCompleteOrderRange(), AllSecondOrderTest::allSecondOrder(), AllSecondOrderTest::allSecondOrderDoNothing(), AllSecondOrderTest::allSecondOrderMixed(), AllSecondOrderTest::allSecondOrderMixedFixing(), AllSecondOrderTest::allSecondOrderRange(), Biharmonic::Biharmonic(), MeshBaseTest::BrokenNeighborMesh(), BoundaryMeshTest::build_mesh(), NodalNeighborsTest::do_test(), main(), WriteVecAndScalar::setupTests(), MultiEvaluablePredTest::test(), MeshSpatialDimensionTest::test2D(), SystemsTest::test2DProjectVectorFE(), MeshFunctionTest::test_bad_gradient_var_with_out_of_mesh_value(), MeshFunctionTest::test_bad_hessian_var_with_out_of_mesh_value(), DistortTest::test_helper_2D(), ProjectSolutionTest::test_partial_project_solution(), MeshFunctionTest::test_subdomain_id_sets(), AllTriTest::testAllTriC0Polygon(), DofMapTest::testArrayDofIndicesWithType(), DofMapTest::testBadElemFECombo(), EquationSystemsTest::testBadVarNames(), BoundaryInfoTest::testBoundaryIDs(), BoundaryInfoTest::testBoundaryOnChildrenBoundaryIDs(), BoundaryInfoTest::testBoundaryOnChildrenBoundarySides(), BoundaryInfoTest::testBoundaryOnChildrenElementsRefineCoarsen(), BoundaryInfoTest::testBoundaryOnChildrenErrors(), BoundaryInfoTest::testBuildNodeListFromSideList(), BoundaryInfoTest::testBuildSideListFromNodeList(), MeshGenerationTest::testBuildSquare(), DofMapTest::testConstraintLoopDetection(), MeshAssignTest::testCopyConstruct(), MeshInputTest::testCopyElementSolutionImpl(), MeshInputTest::testCopyElementVectorImpl(), MeshInputTest::testCopyNodalSolutionImpl(), BoundaryMeshSubdomainTest::testDefaultSubdomain(), MeshDeletionsTest::testDeleteElem(), MeshInputTest::testExodusReadHeader(), MeshExtruderTest::testExtruder(), MixedOrderTest::testFindNeighbors(), BoundaryInfoTest::testInternalBoundary(), MeshSmootherTest::testLaplaceSmoother(), MappedSubdomainPartitionerTest::testMappedSubdomainPartitioner(), BoundaryInfoTest::testMesh(), MeshBaseTest::testMeshBaseVerifyHasCachedElemData(), MeshBaseTest::testMeshBaseVerifyRemovalPreparation(), MeshAssignTest::testMeshMoveAssign(), MeshInputTest::testNemesisReadImpl(), BoundaryMeshSubdomainTest::testPerBoundarySubdomain(), PointLocatorTest::testPlanar(), SystemsTest::testProjectMatrix2D(), SystemsTest::testProjectScalarCoarsening(), SystemsTest::testProjectSquare(), VolumeTest::testQuad4TrueCentroid(), EquationSystemsTest::testRefineThenReinitPreserveFlags(), BoundaryInfoTest::testRenumber(), EquationSystemsTest::testRepartitionThenReinit(), MeshInputTest::testSingleElementImpl(), BoundaryMeshSubdomainTest::testSingleSubdomain(), CheckpointIOTest::testSplitter(), MixedOrderTest::testStitch(), MeshTriangulationTest::testTriangulatorMeshedHoles(), SimplexRefinementTest::testTriRefinement(), MeshSmootherTest::testVariationalSmoother(), MeshInputTest::testVTKPreserveElemIds(), MeshInputTest::testVTKPreserveSubdomainIds(), WriteElemsetData::testWriteImpl(), WriteNodesetData::testWriteImpl(), and WriteSidesetData::testWriteImpl().

◆ surface_octahedron()

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

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

Definition at line 2994 of file mesh_generation.C.

3000{
3001 const Real xavg = (xmin + xmax)/2;
3002 const Real yavg = (ymin + ymax)/2;
3003 const Real zavg = (zmin + zmax)/2;
3004 mesh.add_point(Point(xavg,yavg,zmin), 0);
3005 mesh.add_point(Point(xmax,yavg,zavg), 1);
3006 mesh.add_point(Point(xavg,ymax,zavg), 2);
3007 mesh.add_point(Point(xmin,yavg,zavg), 3);
3008 mesh.add_point(Point(xavg,ymin,zavg), 4);
3009 mesh.add_point(Point(xavg,yavg,zmax), 5);
3010
3011 auto add_tri = [&mesh, flip_tris](std::array<dof_id_type,3> nodes)
3012 {
3013 auto elem = mesh.add_elem(Elem::build(TRI3));
3014 elem->set_node(0, mesh.node_ptr(nodes[0]));
3015 elem->set_node(1, mesh.node_ptr(nodes[1]));
3016 elem->set_node(2, mesh.node_ptr(nodes[2]));
3017 if (flip_tris)
3018 elem->flip(&mesh.get_boundary_info());
3019 };
3020
3021 add_tri({0,2,1});
3022 add_tri({0,3,2});
3023 add_tri({0,4,3});
3024 add_tri({0,1,4});
3025 add_tri({5,4,1});
3026 add_tri({5,3,4});
3027 add_tri({5,2,3});
3028 add_tri({5,1,2});
3029
3030 mesh.prepare_for_use();
3031}

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

Referenced by SimplexRefinementTest::test3DTriRefinement().

Generated on Mon Aug 24 2026 19:01:13 for libMesh by doxygen 1.9.8