libMesh
Functions
libMesh::MeshTools::Modification Namespace Reference

Tools for Mesh modification. More...

Functions

void distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
 Randomly perturb the nodal locations. More...
 
void permute_elements (MeshBase &mesh)
 Randomly permute the nodal ordering of each element (without twisting the element mapping). More...
 
void orient_elements (MeshBase &mesh)
 Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation. More...
 
void redistribute (MeshBase &mesh, const FunctionBase< Real > &mapfunc)
 Deterministically perturb the nodal locations. More...
 
void translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
 Translates the mesh. More...
 
RealTensorValue rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
 Rotates the mesh in the xy plane. More...
 
void scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
 Scales the mesh. More...
 
void all_tri (MeshBase &mesh)
 Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements. More...
 
void smooth (MeshBase &, unsigned int, Real)
 Smooth the mesh with a simple Laplace smoothing algorithm. More...
 
void flatten (MeshBase &mesh)
 Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements. More...
 
void change_boundary_id (MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id)
 Finds any boundary ids that are currently old_id, changes them to new_id. More...
 
void change_subdomain_id (MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)
 Finds any subdomain ids that are currently old_id, changes them to new_id. More...
 

Detailed Description

Tools for Mesh modification.

Author
Benjamin S. Kirk
Date
2004

Function Documentation

◆ all_tri()

void libMesh::MeshTools::Modification::all_tri ( MeshBase &  mesh)

Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.

Note
Only supports coarse / unrefined meshes. A uniformly refined mesh can be used only after a flatten() removes its coarser layers.

Definition at line 429 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::Elem::build_side_ptr(), libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::get_extra_integer(), libMesh::HEX8, libMesh::DofObject::id(), libMesh::index_range(), libMesh::INFEDGE2, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::Elem::mapping_data(), libMesh::Elem::mapping_type(), TIMPI::Communicator::max(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshBase::n_elem(), libMesh::DofObject::n_extra_integers(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_id(), libMesh::Elem::node_ptr(), libMesh::Elem::nodes_on_side(), libMesh::TensorTools::norm(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::point(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PRISM6, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::remote_elem, libMesh::BoundaryInfo::remove(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::TRI7, libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by OverlappingFunctorTest::checkCouplingFunctorTri(), OverlappingFunctorTest::checkCouplingFunctorTriUnifRef(), main(), AllTriTest::test_helper_2D(), AllTriTest::test_helper_3D(), and libMesh::MeshTetInterface::volume_to_surface_mesh().

430 {
431  if (!mesh.is_replicated() && !mesh.is_prepared())
432  mesh.prepare_for_use();
433 
434  // The number of elements in the original mesh before any additions
435  // or deletions.
436  const dof_id_type n_orig_elem = mesh.n_elem();
437  const dof_id_type max_orig_id = mesh.max_elem_id();
438 
439  // We store pointers to the newly created elements in a vector
440  // until they are ready to be added to the mesh. This is because
441  // adding new elements on the fly can cause reallocation and invalidation
442  // of existing mesh element_iterators.
443  std::vector<std::unique_ptr<Elem>> new_elements;
444 
445  unsigned int max_subelems = 1; // in 1D nothing needs to change
446  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
447  max_subelems = 2;
448  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
449  max_subelems = 6;
450 
451  new_elements.reserve (max_subelems*n_orig_elem);
452 
453  // If the original mesh has *side* boundary data, we carry that over
454  // to the new mesh with triangular elements. We currently only
455  // support bringing over side-based BCs to the all-tri mesh, but
456  // that could probably be extended to node and edge-based BCs as
457  // well.
458  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
459 
460  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
461  std::vector<Elem *> new_bndry_elements;
462  std::vector<unsigned short int> new_bndry_sides;
463  std::vector<boundary_id_type> new_bndry_ids;
464 
465  // We may need to add new points if we run into a 1.5th order
466  // element; if we do that on a DistributedMesh in a ghost element then
467  // we will need to fix their ids / unique_ids
468  bool added_new_ghost_point = false;
469 
470  // Iterate over the elements, splitting:
471  // QUADs into pairs of conforming triangles
472  // PYRAMIDs into pairs of conforming tets,
473  // PRISMs into triplets of conforming tets, and
474  // HEXs into quintets or sextets of conforming tets.
475  // We split on the shortest diagonal to give us better
476  // triangle quality in 2D, and we split based on node ids
477  // to guarantee consistency in 3D.
478 
479  // FIXME: This algorithm does not work on refined grids!
480  {
481 #ifdef LIBMESH_ENABLE_UNIQUE_ID
482  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
483 #endif
484 
485  // For avoiding extraneous allocation when building side elements
486  std::unique_ptr<const Elem> elem_side, subside_elem;
487 
488  for (auto & elem : mesh.element_ptr_range())
489  {
490  const ElemType etype = elem->type();
491 
492  // all_tri currently only works on coarse meshes
493  if (elem->parent())
494  libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
495 
496  // The new elements we will split the quad into.
497  // In 3D we may need as many as 6 tets per hex
498  std::array<std::unique_ptr<Elem>, 6> subelem {};
499 
500  switch (etype)
501  {
502  case QUAD4:
503  {
504  subelem[0] = Elem::build(TRI3);
505  subelem[1] = Elem::build(TRI3);
506 
507  // Check for possible edge swap
508  if ((elem->point(0) - elem->point(2)).norm() <
509  (elem->point(1) - elem->point(3)).norm())
510  {
511  subelem[0]->set_node(0, elem->node_ptr(0));
512  subelem[0]->set_node(1, elem->node_ptr(1));
513  subelem[0]->set_node(2, elem->node_ptr(2));
514 
515  subelem[1]->set_node(0, elem->node_ptr(0));
516  subelem[1]->set_node(1, elem->node_ptr(2));
517  subelem[1]->set_node(2, elem->node_ptr(3));
518  }
519 
520  else
521  {
522  subelem[0]->set_node(0, elem->node_ptr(0));
523  subelem[0]->set_node(1, elem->node_ptr(1));
524  subelem[0]->set_node(2, elem->node_ptr(3));
525 
526  subelem[1]->set_node(0, elem->node_ptr(1));
527  subelem[1]->set_node(1, elem->node_ptr(2));
528  subelem[1]->set_node(2, elem->node_ptr(3));
529  }
530 
531 
532  break;
533  }
534 
535  case QUAD8:
536  {
537  if (elem->processor_id() != mesh.processor_id())
538  added_new_ghost_point = true;
539 
540  subelem[0] = Elem::build(TRI6);
541  subelem[1] = Elem::build(TRI6);
542 
543  // Add a new node at the center (vertex average) of the element.
544  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
545  mesh.point(elem->node_id(1)) +
546  mesh.point(elem->node_id(2)) +
547  mesh.point(elem->node_id(3)))/4,
548  DofObject::invalid_id,
549  elem->processor_id());
550 
551  // Check for possible edge swap
552  if ((elem->point(0) - elem->point(2)).norm() <
553  (elem->point(1) - elem->point(3)).norm())
554  {
555  subelem[0]->set_node(0, elem->node_ptr(0));
556  subelem[0]->set_node(1, elem->node_ptr(1));
557  subelem[0]->set_node(2, elem->node_ptr(2));
558  subelem[0]->set_node(3, elem->node_ptr(4));
559  subelem[0]->set_node(4, elem->node_ptr(5));
560  subelem[0]->set_node(5, new_node);
561 
562  subelem[1]->set_node(0, elem->node_ptr(0));
563  subelem[1]->set_node(1, elem->node_ptr(2));
564  subelem[1]->set_node(2, elem->node_ptr(3));
565  subelem[1]->set_node(3, new_node);
566  subelem[1]->set_node(4, elem->node_ptr(6));
567  subelem[1]->set_node(5, elem->node_ptr(7));
568 
569  }
570 
571  else
572  {
573  subelem[0]->set_node(0, elem->node_ptr(3));
574  subelem[0]->set_node(1, elem->node_ptr(0));
575  subelem[0]->set_node(2, elem->node_ptr(1));
576  subelem[0]->set_node(3, elem->node_ptr(7));
577  subelem[0]->set_node(4, elem->node_ptr(4));
578  subelem[0]->set_node(5, new_node);
579 
580  subelem[1]->set_node(0, elem->node_ptr(1));
581  subelem[1]->set_node(1, elem->node_ptr(2));
582  subelem[1]->set_node(2, elem->node_ptr(3));
583  subelem[1]->set_node(3, elem->node_ptr(5));
584  subelem[1]->set_node(4, elem->node_ptr(6));
585  subelem[1]->set_node(5, new_node);
586  }
587 
588  break;
589  }
590 
591  case QUAD9:
592  {
593  subelem[0] = Elem::build(TRI6);
594  subelem[1] = Elem::build(TRI6);
595 
596  // Check for possible edge swap
597  if ((elem->point(0) - elem->point(2)).norm() <
598  (elem->point(1) - elem->point(3)).norm())
599  {
600  subelem[0]->set_node(0, elem->node_ptr(0));
601  subelem[0]->set_node(1, elem->node_ptr(1));
602  subelem[0]->set_node(2, elem->node_ptr(2));
603  subelem[0]->set_node(3, elem->node_ptr(4));
604  subelem[0]->set_node(4, elem->node_ptr(5));
605  subelem[0]->set_node(5, elem->node_ptr(8));
606 
607  subelem[1]->set_node(0, elem->node_ptr(0));
608  subelem[1]->set_node(1, elem->node_ptr(2));
609  subelem[1]->set_node(2, elem->node_ptr(3));
610  subelem[1]->set_node(3, elem->node_ptr(8));
611  subelem[1]->set_node(4, elem->node_ptr(6));
612  subelem[1]->set_node(5, elem->node_ptr(7));
613  }
614 
615  else
616  {
617  subelem[0]->set_node(0, elem->node_ptr(0));
618  subelem[0]->set_node(1, elem->node_ptr(1));
619  subelem[0]->set_node(2, elem->node_ptr(3));
620  subelem[0]->set_node(3, elem->node_ptr(4));
621  subelem[0]->set_node(4, elem->node_ptr(8));
622  subelem[0]->set_node(5, elem->node_ptr(7));
623 
624  subelem[1]->set_node(0, elem->node_ptr(1));
625  subelem[1]->set_node(1, elem->node_ptr(2));
626  subelem[1]->set_node(2, elem->node_ptr(3));
627  subelem[1]->set_node(3, elem->node_ptr(5));
628  subelem[1]->set_node(4, elem->node_ptr(6));
629  subelem[1]->set_node(5, elem->node_ptr(8));
630  }
631 
632  break;
633  }
634 
635  case HEX8:
636  {
637  BoundaryInfo & boundary_info = mesh.get_boundary_info();
638 
639  // Hexes all split into six tetrahedra
640  subelem[0] = Elem::build(TET4);
641  subelem[1] = Elem::build(TET4);
642  subelem[2] = Elem::build(TET4);
643  subelem[3] = Elem::build(TET4);
644  subelem[4] = Elem::build(TET4);
645  subelem[5] = Elem::build(TET4);
646 
647  // On faces, we choose the node with the highest
648  // global id, and we split on the diagonal which
649  // includes that node. This ensures that (even in
650  // parallel, even on distributed meshes) the same
651  // diagonal split will be chosen for elements on either
652  // side of the same quad face.
653  const unsigned int highest_n = highest_vertex_on(elem);
654 
655  // opposing_node[n] is the local node number of the node
656  // on the farthest corner of a hex8 from local node n
657  static const std::array<unsigned int, 8> opposing_node =
658  {6, 7, 4, 5, 2, 3, 0, 1};
659 
660  static const std::vector<std::vector<unsigned int>> sides_opposing_highest =
661  {{2,3,5},{3,4,5},{1,4,5},{1,2,5},{0,2,3},{0,3,4},{0,1,4},{0,1,2}};
662  static const std::vector<std::vector<unsigned int>> nodes_neighboring_highest =
663  {{1,3,4},{0,2,5},{1,3,6},{0,2,7},{0,5,7},{1,4,6},{2,5,7},{3,4,6}};
664 
665  // Start by looking in three directions away from the
666  // highest-id node. In each direction there will be two
667  // different possibilities for the split depending on
668  // how the opposing face nodes are numbered.
669  //
670  // This is tricky enough that I'm not going to worry
671  // about manually keeping tets oriented; we'll just call
672  // orient() on each as we go.
673 
674  unsigned int next_subelem = 0;
675  for (auto side : sides_opposing_highest[highest_n])
676  {
677  const std::vector<unsigned int> nodes_on_side =
678  elem->nodes_on_side(side);
679 
680  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
681 
682  unsigned int split_on_neighbor = false;
683  for (auto n : nodes_neighboring_highest[highest_n])
684  if (dn == n || dn2 == n)
685  {
686  split_on_neighbor = true;
687  break;
688  }
689 
690  // Add one or two elements for each opposing side,
691  // depending on whether the diagonal split there
692  // connects to the neighboring diagonal split or
693  // not.
694  if (split_on_neighbor)
695  {
696  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
697  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
698  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
699  for (auto n : nodes_on_side)
700  if (n != dn && n != dn2)
701  {
702  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
703  break;
704  }
705  subelem[next_subelem]->orient(&boundary_info);
706  ++next_subelem;
707 
708  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
709  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
710  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
711  for (auto n : reverse(nodes_on_side))
712  if (n != dn && n != dn2)
713  {
714  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
715  break;
716  }
717  subelem[next_subelem]->orient(&boundary_info);
718  ++next_subelem;
719  }
720  else
721  {
722  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
723  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
724  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
725  for (auto n : nodes_on_side)
726  for (auto n2 : nodes_neighboring_highest[highest_n])
727  if (n == n2)
728  {
729  subelem[next_subelem]->set_node(3, elem->node_ptr(n));
730  goto break_both_loops;
731  }
732 
733  break_both_loops:
734  subelem[next_subelem]->orient(&boundary_info);
735  ++next_subelem;
736  }
737  }
738 
739  // At this point we've created between 3 and 6 tets.
740  // What's left to do depends on how many.
741 
742  // If we just chopped off three vertices into three
743  // tets, then the best way to split this hex would be
744  // the symmetric five-split. Chop off the opposing
745  // vertex too, and then the remaining interior is our
746  // final tet.
747  if (next_subelem == 3)
748  {
749  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
750  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
751  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
752  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
753  subelem[next_subelem]->orient(&boundary_info);
754  ++next_subelem;
755 
756  subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
757  subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
758  subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
759  subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
760  subelem[next_subelem]->orient(&boundary_info);
761  ++next_subelem;
762 
763  // We don't need the 6th tet after all
764  subelem[next_subelem].reset();
765  ++next_subelem;
766  }
767 
768  // If we just chopped off one (or two) vertices into
769  // tets, then the remaining gap is best (or only) filled
770  // by pairing another tet with each.
771  if (next_subelem == 4 ||
772  next_subelem == 5)
773  {
774  for (auto side : sides_opposing_highest[highest_n])
775  {
776  const std::vector<unsigned int> nodes_on_side =
777  elem->nodes_on_side(side);
778 
779  auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
780 
781  unsigned int split_on_neighbor = false;
782  for (auto n : nodes_neighboring_highest[highest_n])
783  if (dn == n || dn2 == n)
784  {
785  split_on_neighbor = true;
786  break;
787  }
788 
789  // The two !split_on_neighbor sides are where we
790  // need the two remaining tets
791  if (!split_on_neighbor)
792  {
793  subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
794  subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
795  subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
796  subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
797  subelem[next_subelem]->orient(&boundary_info);
798  ++next_subelem;
799  }
800  }
801  }
802 
803  // Whether we got there by creating six tets from the
804  // first for loop or by patching up the split afterward,
805  // we should have considered six tets (possibly
806  // including one deleted one...) at this point.
807  libmesh_assert(next_subelem == 6);
808 
809  break;
810  }
811 
812  case PRISM6:
813  {
814  // Prisms all split into three tetrahedra
815  subelem[0] = Elem::build(TET4);
816  subelem[1] = Elem::build(TET4);
817  subelem[2] = Elem::build(TET4);
818 
819  // Triangular faces are not split.
820 
821  // On quad faces, we choose the node with the highest
822  // global id, and we split on the diagonal which
823  // includes that node. This ensures that (even in
824  // parallel, even on distributed meshes) the same
825  // diagonal split will be chosen for elements on either
826  // side of the same quad face. It also ensures that we
827  // always have a mix of "clockwise" and
828  // "counterclockwise" split faces (two of one and one
829  // of the other on each prism; this is useful since the
830  // alternative all-clockwise or all-counterclockwise
831  // face splittings can't be turned into tets without
832  // adding more nodes
833 
834  // Split on 0-4 diagonal
835  if (split_first_diagonal(elem, 0,4, 1,3))
836  {
837  // Split on 0-5 diagonal
838  if (split_first_diagonal(elem, 0,5, 2,3))
839  {
840  // Split on 1-5 diagonal
841  if (split_first_diagonal(elem, 1,5, 2,4))
842  {
843  subelem[0]->set_node(0, elem->node_ptr(0));
844  subelem[0]->set_node(1, elem->node_ptr(4));
845  subelem[0]->set_node(2, elem->node_ptr(5));
846  subelem[0]->set_node(3, elem->node_ptr(3));
847 
848  subelem[1]->set_node(0, elem->node_ptr(0));
849  subelem[1]->set_node(1, elem->node_ptr(4));
850  subelem[1]->set_node(2, elem->node_ptr(1));
851  subelem[1]->set_node(3, elem->node_ptr(5));
852 
853  subelem[2]->set_node(0, elem->node_ptr(0));
854  subelem[2]->set_node(1, elem->node_ptr(1));
855  subelem[2]->set_node(2, elem->node_ptr(2));
856  subelem[2]->set_node(3, elem->node_ptr(5));
857  }
858  else // Split on 2-4 diagonal
859  {
860  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
861 
862  subelem[0]->set_node(0, elem->node_ptr(0));
863  subelem[0]->set_node(1, elem->node_ptr(4));
864  subelem[0]->set_node(2, elem->node_ptr(5));
865  subelem[0]->set_node(3, elem->node_ptr(3));
866 
867  subelem[1]->set_node(0, elem->node_ptr(0));
868  subelem[1]->set_node(1, elem->node_ptr(4));
869  subelem[1]->set_node(2, elem->node_ptr(2));
870  subelem[1]->set_node(3, elem->node_ptr(5));
871 
872  subelem[2]->set_node(0, elem->node_ptr(0));
873  subelem[2]->set_node(1, elem->node_ptr(1));
874  subelem[2]->set_node(2, elem->node_ptr(2));
875  subelem[2]->set_node(3, elem->node_ptr(4));
876  }
877  }
878  else // Split on 2-3 diagonal
879  {
880  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
881 
882  // 0-4 and 2-3 split implies 2-4 split
883  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
884 
885  subelem[0]->set_node(0, elem->node_ptr(0));
886  subelem[0]->set_node(1, elem->node_ptr(4));
887  subelem[0]->set_node(2, elem->node_ptr(2));
888  subelem[0]->set_node(3, elem->node_ptr(3));
889 
890  subelem[1]->set_node(0, elem->node_ptr(3));
891  subelem[1]->set_node(1, elem->node_ptr(4));
892  subelem[1]->set_node(2, elem->node_ptr(2));
893  subelem[1]->set_node(3, elem->node_ptr(5));
894 
895  subelem[2]->set_node(0, elem->node_ptr(0));
896  subelem[2]->set_node(1, elem->node_ptr(1));
897  subelem[2]->set_node(2, elem->node_ptr(2));
898  subelem[2]->set_node(3, elem->node_ptr(4));
899  }
900  }
901  else // Split on 1-3 diagonal
902  {
903  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
904 
905  // Split on 0-5 diagonal
906  if (split_first_diagonal(elem, 0,5, 2,3))
907  {
908  // 1-3 and 0-5 split implies 1-5 split
909  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
910 
911  subelem[0]->set_node(0, elem->node_ptr(1));
912  subelem[0]->set_node(1, elem->node_ptr(3));
913  subelem[0]->set_node(2, elem->node_ptr(4));
914  subelem[0]->set_node(3, elem->node_ptr(5));
915 
916  subelem[1]->set_node(0, elem->node_ptr(1));
917  subelem[1]->set_node(1, elem->node_ptr(0));
918  subelem[1]->set_node(2, elem->node_ptr(3));
919  subelem[1]->set_node(3, elem->node_ptr(5));
920 
921  subelem[2]->set_node(0, elem->node_ptr(0));
922  subelem[2]->set_node(1, elem->node_ptr(1));
923  subelem[2]->set_node(2, elem->node_ptr(2));
924  subelem[2]->set_node(3, elem->node_ptr(5));
925  }
926  else // Split on 2-3 diagonal
927  {
928  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
929 
930  // Split on 1-5 diagonal
931  if (split_first_diagonal(elem, 1,5, 2,4))
932  {
933  subelem[0]->set_node(0, elem->node_ptr(0));
934  subelem[0]->set_node(1, elem->node_ptr(1));
935  subelem[0]->set_node(2, elem->node_ptr(2));
936  subelem[0]->set_node(3, elem->node_ptr(3));
937 
938  subelem[1]->set_node(0, elem->node_ptr(3));
939  subelem[1]->set_node(1, elem->node_ptr(1));
940  subelem[1]->set_node(2, elem->node_ptr(2));
941  subelem[1]->set_node(3, elem->node_ptr(5));
942 
943  subelem[2]->set_node(0, elem->node_ptr(1));
944  subelem[2]->set_node(1, elem->node_ptr(3));
945  subelem[2]->set_node(2, elem->node_ptr(4));
946  subelem[2]->set_node(3, elem->node_ptr(5));
947  }
948  else // Split on 2-4 diagonal
949  {
950  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
951 
952  subelem[0]->set_node(0, elem->node_ptr(0));
953  subelem[0]->set_node(1, elem->node_ptr(1));
954  subelem[0]->set_node(2, elem->node_ptr(2));
955  subelem[0]->set_node(3, elem->node_ptr(3));
956 
957  subelem[1]->set_node(0, elem->node_ptr(2));
958  subelem[1]->set_node(1, elem->node_ptr(3));
959  subelem[1]->set_node(2, elem->node_ptr(4));
960  subelem[1]->set_node(3, elem->node_ptr(5));
961 
962  subelem[2]->set_node(0, elem->node_ptr(3));
963  subelem[2]->set_node(1, elem->node_ptr(1));
964  subelem[2]->set_node(2, elem->node_ptr(2));
965  subelem[2]->set_node(3, elem->node_ptr(4));
966  }
967  }
968  }
969 
970  break;
971  }
972 
973  case PRISM20:
974  case PRISM21:
975  libmesh_experimental(); // We should upgrade this to TET14...
976  libmesh_fallthrough();
977  case PRISM18:
978  {
979  subelem[0] = Elem::build(TET10);
980  subelem[1] = Elem::build(TET10);
981  subelem[2] = Elem::build(TET10);
982 
983  // Split on 0-4 diagonal
984  if (split_first_diagonal(elem, 0,4, 1,3))
985  {
986  // Split on 0-5 diagonal
987  if (split_first_diagonal(elem, 0,5, 2,3))
988  {
989  // Split on 1-5 diagonal
990  if (split_first_diagonal(elem, 1,5, 2,4))
991  {
992  subelem[0]->set_node(0, elem->node_ptr(0));
993  subelem[0]->set_node(1, elem->node_ptr(4));
994  subelem[0]->set_node(2, elem->node_ptr(5));
995  subelem[0]->set_node(3, elem->node_ptr(3));
996 
997  subelem[0]->set_node(4, elem->node_ptr(15));
998  subelem[0]->set_node(5, elem->node_ptr(13));
999  subelem[0]->set_node(6, elem->node_ptr(17));
1000  subelem[0]->set_node(7, elem->node_ptr(9));
1001  subelem[0]->set_node(8, elem->node_ptr(12));
1002  subelem[0]->set_node(9, elem->node_ptr(14));
1003 
1004  subelem[1]->set_node(0, elem->node_ptr(0));
1005  subelem[1]->set_node(1, elem->node_ptr(4));
1006  subelem[1]->set_node(2, elem->node_ptr(1));
1007  subelem[1]->set_node(3, elem->node_ptr(5));
1008 
1009  subelem[1]->set_node(4, elem->node_ptr(15));
1010  subelem[1]->set_node(5, elem->node_ptr(10));
1011  subelem[1]->set_node(6, elem->node_ptr(6));
1012  subelem[1]->set_node(7, elem->node_ptr(17));
1013  subelem[1]->set_node(8, elem->node_ptr(13));
1014  subelem[1]->set_node(9, elem->node_ptr(16));
1015 
1016  subelem[2]->set_node(0, elem->node_ptr(0));
1017  subelem[2]->set_node(1, elem->node_ptr(1));
1018  subelem[2]->set_node(2, elem->node_ptr(2));
1019  subelem[2]->set_node(3, elem->node_ptr(5));
1020 
1021  subelem[2]->set_node(4, elem->node_ptr(6));
1022  subelem[2]->set_node(5, elem->node_ptr(7));
1023  subelem[2]->set_node(6, elem->node_ptr(8));
1024  subelem[2]->set_node(7, elem->node_ptr(17));
1025  subelem[2]->set_node(8, elem->node_ptr(16));
1026  subelem[2]->set_node(9, elem->node_ptr(11));
1027  }
1028  else // Split on 2-4 diagonal
1029  {
1030  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1031 
1032  subelem[0]->set_node(0, elem->node_ptr(0));
1033  subelem[0]->set_node(1, elem->node_ptr(4));
1034  subelem[0]->set_node(2, elem->node_ptr(5));
1035  subelem[0]->set_node(3, elem->node_ptr(3));
1036 
1037  subelem[0]->set_node(4, elem->node_ptr(15));
1038  subelem[0]->set_node(5, elem->node_ptr(13));
1039  subelem[0]->set_node(6, elem->node_ptr(17));
1040  subelem[0]->set_node(7, elem->node_ptr(9));
1041  subelem[0]->set_node(8, elem->node_ptr(12));
1042  subelem[0]->set_node(9, elem->node_ptr(14));
1043 
1044  subelem[1]->set_node(0, elem->node_ptr(0));
1045  subelem[1]->set_node(1, elem->node_ptr(4));
1046  subelem[1]->set_node(2, elem->node_ptr(2));
1047  subelem[1]->set_node(3, elem->node_ptr(5));
1048 
1049  subelem[1]->set_node(4, elem->node_ptr(15));
1050  subelem[1]->set_node(5, elem->node_ptr(16));
1051  subelem[1]->set_node(6, elem->node_ptr(8));
1052  subelem[1]->set_node(7, elem->node_ptr(17));
1053  subelem[1]->set_node(8, elem->node_ptr(13));
1054  subelem[1]->set_node(9, elem->node_ptr(11));
1055 
1056  subelem[2]->set_node(0, elem->node_ptr(0));
1057  subelem[2]->set_node(1, elem->node_ptr(1));
1058  subelem[2]->set_node(2, elem->node_ptr(2));
1059  subelem[2]->set_node(3, elem->node_ptr(4));
1060 
1061  subelem[2]->set_node(4, elem->node_ptr(6));
1062  subelem[2]->set_node(5, elem->node_ptr(7));
1063  subelem[2]->set_node(6, elem->node_ptr(8));
1064  subelem[2]->set_node(7, elem->node_ptr(15));
1065  subelem[2]->set_node(8, elem->node_ptr(10));
1066  subelem[2]->set_node(9, elem->node_ptr(16));
1067  }
1068  }
1069  else // Split on 2-3 diagonal
1070  {
1071  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1072 
1073  // 0-4 and 2-3 split implies 2-4 split
1074  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1075 
1076  subelem[0]->set_node(0, elem->node_ptr(0));
1077  subelem[0]->set_node(1, elem->node_ptr(4));
1078  subelem[0]->set_node(2, elem->node_ptr(2));
1079  subelem[0]->set_node(3, elem->node_ptr(3));
1080 
1081  subelem[0]->set_node(4, elem->node_ptr(15));
1082  subelem[0]->set_node(5, elem->node_ptr(16));
1083  subelem[0]->set_node(6, elem->node_ptr(8));
1084  subelem[0]->set_node(7, elem->node_ptr(9));
1085  subelem[0]->set_node(8, elem->node_ptr(12));
1086  subelem[0]->set_node(9, elem->node_ptr(17));
1087 
1088  subelem[1]->set_node(0, elem->node_ptr(3));
1089  subelem[1]->set_node(1, elem->node_ptr(4));
1090  subelem[1]->set_node(2, elem->node_ptr(2));
1091  subelem[1]->set_node(3, elem->node_ptr(5));
1092 
1093  subelem[1]->set_node(4, elem->node_ptr(12));
1094  subelem[1]->set_node(5, elem->node_ptr(16));
1095  subelem[1]->set_node(6, elem->node_ptr(17));
1096  subelem[1]->set_node(7, elem->node_ptr(14));
1097  subelem[1]->set_node(8, elem->node_ptr(13));
1098  subelem[1]->set_node(9, elem->node_ptr(11));
1099 
1100  subelem[2]->set_node(0, elem->node_ptr(0));
1101  subelem[2]->set_node(1, elem->node_ptr(1));
1102  subelem[2]->set_node(2, elem->node_ptr(2));
1103  subelem[2]->set_node(3, elem->node_ptr(4));
1104 
1105  subelem[2]->set_node(4, elem->node_ptr(6));
1106  subelem[2]->set_node(5, elem->node_ptr(7));
1107  subelem[2]->set_node(6, elem->node_ptr(8));
1108  subelem[2]->set_node(7, elem->node_ptr(15));
1109  subelem[2]->set_node(8, elem->node_ptr(10));
1110  subelem[2]->set_node(9, elem->node_ptr(16));
1111  }
1112  }
1113  else // Split on 1-3 diagonal
1114  {
1115  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1116 
1117  // Split on 0-5 diagonal
1118  if (split_first_diagonal(elem, 0,5, 2,3))
1119  {
1120  // 1-3 and 0-5 split implies 1-5 split
1121  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1122 
1123  subelem[0]->set_node(0, elem->node_ptr(1));
1124  subelem[0]->set_node(1, elem->node_ptr(3));
1125  subelem[0]->set_node(2, elem->node_ptr(4));
1126  subelem[0]->set_node(3, elem->node_ptr(5));
1127 
1128  subelem[0]->set_node(4, elem->node_ptr(15));
1129  subelem[0]->set_node(5, elem->node_ptr(12));
1130  subelem[0]->set_node(6, elem->node_ptr(10));
1131  subelem[0]->set_node(7, elem->node_ptr(16));
1132  subelem[0]->set_node(8, elem->node_ptr(14));
1133  subelem[0]->set_node(9, elem->node_ptr(13));
1134 
1135  subelem[1]->set_node(0, elem->node_ptr(1));
1136  subelem[1]->set_node(1, elem->node_ptr(0));
1137  subelem[1]->set_node(2, elem->node_ptr(3));
1138  subelem[1]->set_node(3, elem->node_ptr(5));
1139 
1140  subelem[1]->set_node(4, elem->node_ptr(6));
1141  subelem[1]->set_node(5, elem->node_ptr(9));
1142  subelem[1]->set_node(6, elem->node_ptr(15));
1143  subelem[1]->set_node(7, elem->node_ptr(16));
1144  subelem[1]->set_node(8, elem->node_ptr(17));
1145  subelem[1]->set_node(9, elem->node_ptr(14));
1146 
1147  subelem[2]->set_node(0, elem->node_ptr(0));
1148  subelem[2]->set_node(1, elem->node_ptr(1));
1149  subelem[2]->set_node(2, elem->node_ptr(2));
1150  subelem[2]->set_node(3, elem->node_ptr(5));
1151 
1152  subelem[2]->set_node(4, elem->node_ptr(6));
1153  subelem[2]->set_node(5, elem->node_ptr(7));
1154  subelem[2]->set_node(6, elem->node_ptr(8));
1155  subelem[2]->set_node(7, elem->node_ptr(17));
1156  subelem[2]->set_node(8, elem->node_ptr(16));
1157  subelem[2]->set_node(9, elem->node_ptr(11));
1158  }
1159  else // Split on 2-3 diagonal
1160  {
1161  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1162 
1163  // Split on 1-5 diagonal
1164  if (split_first_diagonal(elem, 1,5, 2,4))
1165  {
1166  subelem[0]->set_node(0, elem->node_ptr(0));
1167  subelem[0]->set_node(1, elem->node_ptr(1));
1168  subelem[0]->set_node(2, elem->node_ptr(2));
1169  subelem[0]->set_node(3, elem->node_ptr(3));
1170 
1171  subelem[0]->set_node(4, elem->node_ptr(6));
1172  subelem[0]->set_node(5, elem->node_ptr(7));
1173  subelem[0]->set_node(6, elem->node_ptr(8));
1174  subelem[0]->set_node(7, elem->node_ptr(9));
1175  subelem[0]->set_node(8, elem->node_ptr(15));
1176  subelem[0]->set_node(9, elem->node_ptr(17));
1177 
1178  subelem[1]->set_node(0, elem->node_ptr(3));
1179  subelem[1]->set_node(1, elem->node_ptr(1));
1180  subelem[1]->set_node(2, elem->node_ptr(2));
1181  subelem[1]->set_node(3, elem->node_ptr(5));
1182 
1183  subelem[1]->set_node(4, elem->node_ptr(15));
1184  subelem[1]->set_node(5, elem->node_ptr(7));
1185  subelem[1]->set_node(6, elem->node_ptr(17));
1186  subelem[1]->set_node(7, elem->node_ptr(14));
1187  subelem[1]->set_node(8, elem->node_ptr(16));
1188  subelem[1]->set_node(9, elem->node_ptr(11));
1189 
1190  subelem[2]->set_node(0, elem->node_ptr(1));
1191  subelem[2]->set_node(1, elem->node_ptr(3));
1192  subelem[2]->set_node(2, elem->node_ptr(4));
1193  subelem[2]->set_node(3, elem->node_ptr(5));
1194 
1195  subelem[2]->set_node(4, elem->node_ptr(15));
1196  subelem[2]->set_node(5, elem->node_ptr(12));
1197  subelem[2]->set_node(6, elem->node_ptr(10));
1198  subelem[2]->set_node(7, elem->node_ptr(16));
1199  subelem[2]->set_node(8, elem->node_ptr(14));
1200  subelem[2]->set_node(9, elem->node_ptr(13));
1201  }
1202  else // Split on 2-4 diagonal
1203  {
1204  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1205 
1206  subelem[0]->set_node(0, elem->node_ptr(0));
1207  subelem[0]->set_node(1, elem->node_ptr(1));
1208  subelem[0]->set_node(2, elem->node_ptr(2));
1209  subelem[0]->set_node(3, elem->node_ptr(3));
1210 
1211  subelem[0]->set_node(4, elem->node_ptr(6));
1212  subelem[0]->set_node(5, elem->node_ptr(7));
1213  subelem[0]->set_node(6, elem->node_ptr(8));
1214  subelem[0]->set_node(7, elem->node_ptr(9));
1215  subelem[0]->set_node(8, elem->node_ptr(15));
1216  subelem[0]->set_node(9, elem->node_ptr(17));
1217 
1218  subelem[1]->set_node(0, elem->node_ptr(2));
1219  subelem[1]->set_node(1, elem->node_ptr(3));
1220  subelem[1]->set_node(2, elem->node_ptr(4));
1221  subelem[1]->set_node(3, elem->node_ptr(5));
1222 
1223  subelem[1]->set_node(4, elem->node_ptr(17));
1224  subelem[1]->set_node(5, elem->node_ptr(12));
1225  subelem[1]->set_node(6, elem->node_ptr(16));
1226  subelem[1]->set_node(7, elem->node_ptr(11));
1227  subelem[1]->set_node(8, elem->node_ptr(14));
1228  subelem[1]->set_node(9, elem->node_ptr(13));
1229 
1230  subelem[2]->set_node(0, elem->node_ptr(3));
1231  subelem[2]->set_node(1, elem->node_ptr(1));
1232  subelem[2]->set_node(2, elem->node_ptr(2));
1233  subelem[2]->set_node(3, elem->node_ptr(4));
1234 
1235  subelem[2]->set_node(4, elem->node_ptr(15));
1236  subelem[2]->set_node(5, elem->node_ptr(7));
1237  subelem[2]->set_node(6, elem->node_ptr(17));
1238  subelem[2]->set_node(7, elem->node_ptr(12));
1239  subelem[2]->set_node(8, elem->node_ptr(10));
1240  subelem[2]->set_node(9, elem->node_ptr(16));
1241  }
1242  }
1243  }
1244 
1245  break;
1246  }
1247 
1248  // No need to split elements that are already simplicial:
1249  case EDGE2:
1250  case EDGE3:
1251  case EDGE4:
1252  case TRI3:
1253  case TRI6:
1254  case TRI7:
1255  case TET4:
1256  case TET10:
1257  case TET14:
1258  case INFEDGE2:
1259  // No way to split infinite quad/prism elements, so
1260  // hopefully no need to
1261  case INFQUAD4:
1262  case INFQUAD6:
1263  case INFPRISM6:
1264  case INFPRISM12:
1265  continue;
1266  // If we're left with an unimplemented hex we're probably
1267  // out of luck. TODO: implement hexes
1268  default:
1269  {
1270  libMesh::err << "Error, encountered unimplemented element "
1271  << Utility::enum_to_string<ElemType>(etype)
1272  << " in MeshTools::Modification::all_tri()..."
1273  << std::endl;
1274  libmesh_not_implemented();
1275  }
1276  } // end switch (etype)
1277 
1278  // Be sure the correct data is set for all subelems.
1279  const unsigned int nei = elem->n_extra_integers();
1280  for (unsigned int i=0; i != max_subelems; ++i)
1281  if (subelem[i]) {
1282  subelem[i]->processor_id() = elem->processor_id();
1283  subelem[i]->subdomain_id() = elem->subdomain_id();
1284 
1285  // Copy any extra element data. Since the subelements
1286  // haven't been added to the mesh yet any allocation has
1287  // to be done manually.
1288  subelem[i]->add_extra_integers(nei);
1289  for (unsigned int ei=0; ei != nei; ++ei)
1290  subelem[ei]->set_extra_integer(ei, elem->get_extra_integer(ei));
1291 
1292 
1293  // Copy any mapping data.
1294  subelem[i]->set_mapping_type(elem->mapping_type());
1295  subelem[i]->set_mapping_data(elem->mapping_data());
1296  }
1297 
1298  // On a mesh with boundary data, we need to move that data to
1299  // the new elements.
1300 
1301  // On a mesh which is distributed, we need to move
1302  // remote_elem links to the new elements.
1303  bool mesh_is_serial = mesh.is_serial();
1304 
1305  if (mesh_has_boundary_data || !mesh_is_serial)
1306  {
1307  // Container to key boundary IDs handed back by the BoundaryInfo object.
1308  std::vector<boundary_id_type> bc_ids;
1309 
1310  for (auto sn : elem->side_index_range())
1311  {
1312  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1313 
1314  if (bc_ids.empty() && elem->neighbor_ptr(sn) != remote_elem)
1315  continue;
1316 
1317  // Make a sorted list of node ids for elem->side(sn)
1318  elem->build_side_ptr(elem_side, sn);
1319  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1320  for (unsigned int esn=0,
1321  n_esn = cast_int<unsigned int>(elem_side_nodes.size());
1322  esn != n_esn; ++esn)
1323  elem_side_nodes[esn] = elem_side->node_id(esn);
1324  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1325 
1326  for (unsigned int i=0; i != max_subelems; ++i)
1327  if (subelem[i])
1328  {
1329  for (auto subside : subelem[i]->side_index_range())
1330  {
1331  subelem[i]->build_side_ptr(subside_elem, subside);
1332 
1333  // Make a list of *vertex* node ids for this subside, see if they are all present
1334  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1335  // subelem[i]->key(subside) in the Prism cases, since the new side is
1336  // a different type. Note 2: we only use vertex nodes since, in the future,
1337  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1338  // original face will not contain the mid-edge node.
1339  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1340  for (unsigned int ssn=0,
1341  n_ssn = cast_int<unsigned int>(subside_nodes.size());
1342  ssn != n_ssn; ++ssn)
1343  subside_nodes[ssn] = subside_elem->node_id(ssn);
1344  std::sort(subside_nodes.begin(), subside_nodes.end());
1345 
1346  // std::includes returns true if every element of the second sorted range is
1347  // contained in the first sorted range.
1348  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1349  subside_nodes.begin(), subside_nodes.end()))
1350  {
1351  for (const auto & b_id : bc_ids)
1352  if (b_id != BoundaryInfo::invalid_id)
1353  {
1354  new_bndry_ids.push_back(b_id);
1355  new_bndry_elements.push_back(subelem[i].get());
1356  new_bndry_sides.push_back(subside);
1357  }
1358 
1359  // If the original element had a RemoteElem neighbor on side 'sn',
1360  // then the subelem has one on side 'subside'.
1361  if (elem->neighbor_ptr(sn) == remote_elem)
1362  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1363  }
1364  }
1365  } // end for loop over subelem
1366  } // end for loop over sides
1367 
1368  // Remove the original element from the BoundaryInfo structure.
1369  mesh.get_boundary_info().remove(elem);
1370 
1371  } // end if (mesh_has_boundary_data)
1372 
1373  // Determine new IDs for the split elements which will be
1374  // the same on all processors, therefore keeping the Mesh
1375  // in sync. Note: we offset the new IDs by max_orig_id to
1376  // avoid overwriting any of the original IDs.
1377  for (unsigned int i=0; i != max_subelems; ++i)
1378  if (subelem[i])
1379  {
1380  // Determine new IDs for the split elements which will be
1381  // the same on all processors, therefore keeping the Mesh
1382  // in sync. Note: we offset the new IDs by the max of the
1383  // pre-existing ids to avoid conflicting with originals.
1384  subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1385 
1386 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1387  subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);
1388 #endif
1389 
1390  // Prepare to add the newly-created simplices
1391  new_elements.push_back(std::move(subelem[i]));
1392  }
1393 
1394  // Delete the original element
1395  mesh.delete_elem(elem);
1396  } // End for loop over elements
1397  } // end scope
1398 
1399 
1400  // Now, iterate over the new elements vector, and add them each to
1401  // the Mesh.
1402  for (auto & elem : new_elements)
1403  mesh.add_elem(std::move(elem));
1404 
1405  if (mesh_has_boundary_data)
1406  {
1407  // If the old mesh had boundary data, the new mesh better have
1408  // some. However, we can't assert that the size of
1409  // new_bndry_elements vector is > 0, since we may not have split
1410  // any elements actually on the boundary. We also can't assert
1411  // that the original number of boundary sides is equal to the
1412  // sum of the boundary sides currently in the mesh and the
1413  // newly-added boundary sides, since in 3D, we may have split a
1414  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1415  // too picky about the actual number of BCs, and just assert that
1416  // there are some, somewhere.
1417 #ifndef NDEBUG
1418  bool nbe_nonempty = new_bndry_elements.size();
1419  mesh.comm().max(nbe_nonempty);
1420  libmesh_assert(nbe_nonempty ||
1421  mesh.get_boundary_info().n_boundary_conds()>0);
1422 #endif
1423 
1424  // We should also be sure that the lengths of the new boundary data vectors
1425  // are all the same.
1426  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1427  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1428 
1429  // Add the new boundary info to the mesh
1430  for (auto s : index_range(new_bndry_elements))
1431  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1432  new_bndry_sides[s],
1433  new_bndry_ids[s]);
1434  }
1435 
1436  // In a DistributedMesh any newly added ghost node ids may be
1437  // inconsistent, and unique_ids of newly added ghost nodes remain
1438  // unset.
1439  // make_nodes_parallel_consistent() will fix all this.
1440  if (!mesh.is_serial())
1441  {
1442  mesh.comm().max(added_new_ghost_point);
1443 
1444  if (added_new_ghost_point)
1445  MeshCommunication().make_nodes_parallel_consistent (mesh);
1446  }
1447 
1448  // Prepare the newly created mesh for use.
1449  mesh.prepare_for_use();
1450 
1451  // Let the new_elements and new_bndry_elements vectors go out of scope.
1452 }
libMesh::err
OStreamProxy err
libMesh::BoundaryInfo::n_boundary_conds
std::size_t n_boundary_conds() const
Definition: boundary_info.C:2594
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::MeshBase::is_prepared
bool is_prepared() const
Definition: mesh_base.C:1042
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::MeshBase::parallel_max_unique_id
virtual unique_id_type parallel_max_unique_id() const =0
libMesh::MeshBase::prepare_for_use
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:811
libMesh::BoundaryInfo::remove
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
Definition: boundary_info.C:1786
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:97
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:1348
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:175
libMesh::MeshBase::add_point
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
libMesh::MeshBase::is_serial
virtual bool is_serial() const
Definition: mesh_base.h:353
libMesh::MeshBase::delete_elem
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
libMesh::MeshCommunication
This is the MeshCommunication class.
Definition: mesh_communication.h:53
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libMesh::MeshBase::max_elem_id
virtual dof_id_type max_elem_id() const =0
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::libmesh_assert
libmesh_assert(ctx)
TIMPI::Communicator::max
void max(const T &r, T &o, Request &req) const
libMesh::TensorTools::norm
auto norm(const T &a)
Definition: tensor_tools.h:74
libMesh::MeshCommunication::make_nodes_parallel_consistent
void make_nodes_parallel_consistent(MeshBase &)
Copy processor_ids and ids on ghost nodes from their local processors.
Definition: mesh_communication.C:2042
libMesh::MeshBase::is_replicated
virtual bool is_replicated() const
Definition: mesh_base.h:375
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:1210
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:422
libMesh::MeshBase::point
virtual const Point & point(const dof_id_type i) const =0
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::ParallelObject::processor_id
processor_id_type processor_id() const
Definition: parallel_object.h:114
libMesh::unique_id_type
uint8_t unique_id_type
Definition: id_types.h:86
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:150
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57

◆ change_boundary_id()

void libMesh::MeshTools::Modification::change_boundary_id ( MeshBase &  mesh,
const boundary_id_type  old_id,
const boundary_id_type  new_id 
)

Finds any boundary ids that are currently old_id, changes them to new_id.

Definition at line 1741 of file mesh_modification.C.

References libMesh::MeshBase::get_boundary_info(), mesh, and libMesh::BoundaryInfo::renumber_id().

Referenced by MeshStitchTest::renameAndShift().

1744 {
1745  // This is just a shim around the member implementation, now
1746  mesh.get_boundary_info().renumber_id(old_id, new_id);
1747 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:175
libMesh::BoundaryInfo::renumber_id
void renumber_id(boundary_id_type old_id, boundary_id_type new_id)
Changes all entities (nodes, sides, edges, shellfaces) with boundary id old_id to instead be labeled ...
Definition: boundary_info.C:2070

◆ change_subdomain_id()

void libMesh::MeshTools::Modification::change_subdomain_id ( MeshBase &  mesh,
const subdomain_id_type  old_id,
const subdomain_id_type  new_id 
)

Finds any subdomain ids that are currently old_id, changes them to new_id.

Definition at line 1751 of file mesh_modification.C.

References mesh, libMesh::Elem::subdomain_id(), and libMesh::MeshBase::unset_has_cached_elem_data().

1754 {
1755  if (old_id == new_id)
1756  {
1757  // If the IDs are the same, this is a no-op.
1758  return;
1759  }
1760 
1761  for (auto & elem : mesh.element_ptr_range())
1762  {
1763  if (elem->subdomain_id() == old_id)
1764  elem->subdomain_id() = new_id;
1765  }
1766 
1767  // We just invalidated mesh.get_subdomain_ids(), but it might not be
1768  // efficient to fix that here.
1769  mesh.unset_has_cached_elem_data();
1770 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::unset_has_cached_elem_data
void unset_has_cached_elem_data()
Tells this we have done some operation (e.g.
Definition: mesh_base.h:287

◆ distort()

void libMesh::MeshTools::Modification::distort ( MeshBase &  mesh,
const Real  factor,
const bool  perturb_boundary = false 
)

Randomly perturb the nodal locations.

This function will move each node factor fraction of its minimum neighboring node separation distance. Nodes on the boundary are not moved by default, however they may be by setting the flag perturb_boundary true.

Definition at line 141 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::hmin(), libMesh::libmesh_assert(), libMesh::make_range(), libMesh::MeshBase::max_node_id(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::Elem::node_ref_range(), libMesh::MeshBase::query_node_ptr(), and libMesh::TypeVector< T >::unit().

Referenced by main(), DistortTest::perturb_and_check(), VolumeTest::test_true_centroid_and_volume(), and VolumeTest::testQuad4TrueCentroid().

144 {
145  libmesh_assert (mesh.n_nodes());
146  libmesh_assert (mesh.n_elem());
147  libmesh_assert ((factor >= 0.) && (factor <= 1.));
148 
149  LOG_SCOPE("distort()", "MeshTools::Modification");
150 
151  // If we are not perturbing boundary nodes, make a
152  // quickly-searchable list of node ids we can check against.
153  std::unordered_set<dof_id_type> boundary_node_ids;
154  if (!perturb_boundary)
155  boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
156 
157  // Now calculate the minimum distance to
158  // neighboring nodes for each node.
159  // hmin holds these distances.
160  std::vector<float> hmin (mesh.max_node_id(),
161  std::numeric_limits<float>::max());
162 
163  for (const auto & elem : mesh.active_element_ptr_range())
164  for (auto & n : elem->node_ref_range())
165  hmin[n.id()] = std::min(hmin[n.id()],
166  static_cast<float>(elem->hmin()));
167 
168  // Now actually move the nodes
169  {
170  const unsigned int seed = 123456;
171 
172  // seed the random number generator.
173  // We'll loop from 1 to n_nodes on every processor, even those
174  // that don't have a particular node, so that the pseudorandom
175  // numbers will be the same everywhere.
176  std::srand(seed);
177 
178  // If the node is on the boundary or
179  // the node is not used by any element (hmin[n]<1.e20)
180  // then we should not move it.
181  // [Note: Testing for (in)equality might be wrong
182  // (different types, namely float and double)]
183  for (auto n : make_range(mesh.max_node_id()))
184  if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
185  {
186  // the direction, random but unit normalized
187  Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
188  (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
189  ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
190 
191  dir(0) = (dir(0)-.5)*2.;
192 #if LIBMESH_DIM > 1
193  if (mesh.mesh_dimension() > 1)
194  dir(1) = (dir(1)-.5)*2.;
195 #endif
196 #if LIBMESH_DIM > 2
197  if (mesh.mesh_dimension() == 3)
198  dir(2) = (dir(2)-.5)*2.;
199 #endif
200 
201  dir = dir.unit();
202 
203  Node * node = mesh.query_node_ptr(n);
204  if (!node)
205  continue;
206 
207  (*node)(0) += dir(0)*factor*hmin[n];
208 #if LIBMESH_DIM > 1
209  if (mesh.mesh_dimension() > 1)
210  (*node)(1) += dir(1)*factor*hmin[n];
211 #endif
212 #if LIBMESH_DIM > 2
213  if (mesh.mesh_dimension() == 3)
214  (*node)(2) += dir(2)*factor*hmin[n];
215 #endif
216  }
217  }
218 
219  // We haven't changed any topology, but just changing geometry could
220  // have invalidated a point locator.
221  mesh.clear_point_locator();
222 }
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshTools::find_boundary_nodes
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:524
libMesh::TypeVector::unit
TypeVector< T > unit() const
Definition: type_vector.h:1109
libMesh::MeshBase::query_node_ptr
virtual const Node * query_node_ptr(const dof_id_type i) const =0
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1809
libMesh::libmesh_assert
libmesh_assert(ctx)
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:173
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:422
libMesh::MeshBase::max_node_id
virtual dof_id_type max_node_id() const =0
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0

◆ flatten()

void libMesh::MeshTools::Modification::flatten ( MeshBase &  mesh)

Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements.

This is useful when you want to write out a uniformly-refined grid to be treated later as an initial mesh.

Note
Many functions in LibMesh assume a conforming (with no hanging nodes) grid exists at some level, so you probably only want to do this on meshes which have been uniformly refined.

Definition at line 1615 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::get_extra_integer(), libMesh::DofObject::id(), libMesh::index_range(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_prepared(), libMesh::MeshBase::is_replicated(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::Elem::mapping_data(), libMesh::Elem::mapping_type(), mesh, libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::DofObject::n_extra_integers(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::remote_elem, libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by main().

1616 {
1617  libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
1618 
1619  // Algorithm:
1620  // .) For each active element in the mesh: construct a
1621  // copy which is the same in every way *except* it is
1622  // a level 0 element. Store the pointers to these in
1623  // a separate vector. Save any boundary information as well.
1624  // Delete the active element from the mesh.
1625  // .) Loop over all (remaining) elements in the mesh, delete them.
1626  // .) Add the level-0 copies back to the mesh
1627 
1628  // Temporary storage for new element pointers
1629  std::vector<std::unique_ptr<Elem>> new_elements;
1630 
1631  // BoundaryInfo Storage for element ids, sides, and BC ids
1632  std::vector<Elem *> saved_boundary_elements;
1633  std::vector<boundary_id_type> saved_bc_ids;
1634  std::vector<unsigned short int> saved_bc_sides;
1635 
1636  // Container to catch boundary ids passed back by BoundaryInfo
1637  std::vector<boundary_id_type> bc_ids;
1638 
1639  // Reserve a reasonable amt. of space for each
1640  new_elements.reserve(mesh.n_active_elem());
1641  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1642  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1643  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1644 
1645  for (auto & elem : mesh.active_element_ptr_range())
1646  {
1647  // Make a new element of the same type
1648  auto copy = Elem::build(elem->type());
1649 
1650  // Set node pointers (they still point to nodes in the original mesh)
1651  for (auto n : elem->node_index_range())
1652  copy->set_node(n, elem->node_ptr(n));
1653 
1654  // Copy over ids
1655  copy->processor_id() = elem->processor_id();
1656  copy->subdomain_id() = elem->subdomain_id();
1657 
1658  // Retain the original element's ID(s) as well, otherwise
1659  // the Mesh may try to create them for you...
1660  copy->set_id( elem->id() );
1661 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1662  copy->set_unique_id(elem->unique_id());
1663 #endif
1664 
1665  // This element could have boundary info or DistributedMesh
1666  // remote_elem links as well. We need to save the (elem,
1667  // side, bc_id) triples and those links
1668  for (auto s : elem->side_index_range())
1669  {
1670  if (elem->neighbor_ptr(s) == remote_elem)
1671  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1672 
1673  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1674  for (const auto & bc_id : bc_ids)
1675  if (bc_id != BoundaryInfo::invalid_id)
1676  {
1677  saved_boundary_elements.push_back(copy.get());
1678  saved_bc_ids.push_back(bc_id);
1679  saved_bc_sides.push_back(s);
1680  }
1681  }
1682 
1683  // Copy any extra element data. Since the copy hasn't been
1684  // added to the mesh yet any allocation has to be done manually.
1685  const unsigned int nei = elem->n_extra_integers();
1686  copy->add_extra_integers(nei);
1687  for (unsigned int i=0; i != nei; ++i)
1688  copy->set_extra_integer(i, elem->get_extra_integer(i));
1689 
1690  // Copy any mapping data.
1691  copy->set_mapping_type(elem->mapping_type());
1692  copy->set_mapping_data(elem->mapping_data());
1693 
1694  // We're done with this element
1695  mesh.delete_elem(elem);
1696 
1697  // But save the copy
1698  new_elements.push_back(std::move(copy));
1699  }
1700 
1701  // Make sure we saved the same number of boundary conditions
1702  // in each vector.
1703  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1704  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1705 
1706  // Loop again, delete any remaining elements
1707  for (auto & elem : mesh.element_ptr_range())
1708  mesh.delete_elem(elem);
1709 
1710  // Add the copied (now level-0) elements back to the mesh
1711  for (auto & new_elem : new_elements)
1712  {
1713  // Save the original ID, because the act of adding the Elem can
1714  // change new_elem's id!
1715  dof_id_type orig_id = new_elem->id();
1716 
1717  Elem * added_elem = mesh.add_elem(std::move(new_elem));
1718 
1719  // If the Elem, as it was re-added to the mesh, now has a
1720  // different ID (this is unlikely, so it's just an assert)
1721  // the boundary information will no longer be correct.
1722  libmesh_assert_equal_to (orig_id, added_elem->id());
1723 
1724  // Avoid compiler warnings in opt mode.
1725  libmesh_ignore(added_elem, orig_id);
1726  }
1727 
1728  // Finally, also add back the saved boundary information
1729  for (auto e : index_range(saved_boundary_elements))
1730  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1731  saved_bc_sides[e],
1732  saved_bc_ids[e]);
1733 
1734  // Trim unused and renumber nodes and elements
1735  mesh.prepare_for_use();
1736 }
libMesh::BoundaryInfo::n_boundary_conds
std::size_t n_boundary_conds() const
Definition: boundary_info.C:2594
libMesh::MeshBase::is_prepared
bool is_prepared() const
Definition: mesh_base.C:1042
libMesh::MeshBase::n_active_elem
virtual dof_id_type n_active_elem() const =0
libMesh::MeshBase::prepare_for_use
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:811
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
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:1348
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:175
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:535
libMesh::MeshBase::delete_elem
virtual void delete_elem(Elem *e)=0
Removes element e from the mesh.
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:819
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::MeshBase::is_replicated
virtual bool is_replicated() const
Definition: mesh_base.h:375
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:1210
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:150
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57

◆ orient_elements()

void libMesh::MeshTools::Modification::orient_elements ( MeshBase &  mesh)

Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientation.

This function does not currently handle meshes with any element refinement.

Definition at line 266 of file mesh_modification.C.

References libMesh::MeshBase::get_boundary_info(), mesh, libMesh::MeshTools::n_levels(), and libMesh::Elem::orient().

Referenced by ElemTest< elem_type >::test_orient_elements().

267 {
268  LOG_SCOPE("orient_elements()", "MeshTools::Modification");
269 
270  // We don't yet support doing orient() on a parent element, which
271  // would require us to consistently orient all its children and
272  // give them different local child numbers.
273  unsigned int n_levels = MeshTools::n_levels(mesh);
274  if (n_levels > 1)
275  libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
276 
277  BoundaryInfo & boundary_info = mesh.get_boundary_info();
278  for (auto elem : mesh.element_ptr_range())
279  elem->orient(&boundary_info);
280 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:175
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:809
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57

◆ permute_elements()

void libMesh::MeshTools::Modification::permute_elements ( MeshBase &  mesh)

Randomly permute the nodal ordering of each element (without twisting the element mapping).

This is useful for regression testing with a variety of element orientations.

This function does not currently handle meshes with any element refinement.

This function does not currently permute BoundaryInfo data associated with element sides, which will likely be scrambled by the permutation.

Definition at line 226 of file mesh_modification.C.

References libMesh::make_range(), libMesh::MeshBase::max_elem_id(), mesh, libMesh::MeshTools::n_levels(), libMesh::Elem::n_permutations(), libMesh::Elem::permute(), and libMesh::MeshBase::query_elem_ptr().

Referenced by main(), and FETestBase< order, family, elem_type, N_x >::setUp().

227 {
228  LOG_SCOPE("permute_elements()", "MeshTools::Modification");
229 
230  // We don't yet support doing permute() on a parent element, which
231  // would require us to consistently permute all its children and
232  // give them different local child numbers.
233  unsigned int n_levels = MeshTools::n_levels(mesh);
234  if (n_levels > 1)
235  libmesh_error();
236 
237  const unsigned int seed = 123456;
238 
239  // seed the random number generator.
240  // We'll loop from 1 to max_elem_id on every processor, even those
241  // that don't have a particular element, so that the pseudorandom
242  // numbers will be the same everywhere.
243  std::srand(seed);
244 
245 
246  for (auto e_id : make_range(mesh.max_elem_id()))
247  {
248  int my_rand = std::rand();
249 
250  Elem * elem = mesh.query_elem_ptr(e_id);
251 
252  if (!elem)
253  continue;
254 
255  const unsigned int max_permutation = elem->n_permutations();
256  if (!max_permutation)
257  continue;
258 
259  const unsigned int perm = my_rand % max_permutation;
260 
261  elem->permute(perm);
262  }
263 }
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:809
libMesh::MeshBase::max_elem_id
virtual dof_id_type max_elem_id() const =0
libMesh::MeshBase::query_elem_ptr
virtual const Elem * query_elem_ptr(const dof_id_type i) const =0
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:173
libMesh::Elem::permute
virtual void permute(unsigned int perm_num)=0
Permutes the element (by swapping node and neighbor pointers) according to the specified index...
libMesh::Elem::n_permutations
virtual unsigned int n_permutations() const =0
Returns the number of independent permutations of element nodes - e.g.

◆ redistribute()

void libMesh::MeshTools::Modification::redistribute ( MeshBase &  mesh,
const FunctionBase< Real > &  mapfunc 
)

Deterministically perturb the nodal locations.

This function will move each node from it's current x/y/z coordinates to a new x/y/z coordinate given by the first LIBMESH_DIM components of the specified function mapfunc

Nodes on the boundary are also moved.

Currently, non-vertex nodes are moved in the same way as vertex nodes, according to (newx,newy,newz) = mapfunc(x,y,z). This behavior is often suboptimal for higher order geometries and may be subject to change in future libMesh versions.

Definition at line 284 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), libMesh::FunctionBase< Output >::clone(), libMesh::libmesh_assert(), mesh, libMesh::MeshBase::n_elem(), and libMesh::MeshBase::n_nodes().

Referenced by libMesh::MeshTools::Generation::build_cube(), FETestBase< order, family, elem_type, N_x >::setUp(), MeshSmootherTest::testLaplaceSmoother(), MeshSmootherTest::testVariationalSmoother(), and MeshSmootherTest::testVariationalSmootherRegression().

286 {
287  libmesh_assert (mesh.n_nodes());
288  libmesh_assert (mesh.n_elem());
289 
290  LOG_SCOPE("redistribute()", "MeshTools::Modification");
291 
292  DenseVector<Real> output_vec(LIBMESH_DIM);
293 
294  // FIXME - we should thread this later.
295  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
296 
297  for (auto & node : mesh.node_ptr_range())
298  {
299  (*myfunc)(*node, output_vec);
300 
301  (*node)(0) = output_vec(0);
302 #if LIBMESH_DIM > 1
303  (*node)(1) = output_vec(1);
304 #endif
305 #if LIBMESH_DIM > 2
306  (*node)(2) = output_vec(2);
307 #endif
308  }
309 
310  // We haven't changed any topology, but just changing geometry could
311  // have invalidated a point locator.
312  mesh.clear_point_locator();
313 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1809
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::FunctionBase::clone
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0

◆ rotate()

RealTensorValue libMesh::MeshTools::Modification::rotate ( MeshBase &  mesh,
const Real  phi,
const Real  theta = 0.,
const Real  psi = 0. 
)

Rotates the mesh in the xy plane.

The rotation is counter-clock-wise (mathematical definition). The angle is in degrees (360 make a full circle) Rotates the mesh in 3D space. Here the standard Euler angles are adopted (http://mathworld.wolfram.com/EulerAngles.html) The angles are in degrees (360 make a full circle)

Returns
the 3x3 rotation matrix implied by (phi, theta, psi)

Definition at line 356 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), libMesh::TensorValue< T >::intrinsic_rotation_matrix(), libMesh::libmesh_ignore(), mesh, and libMesh::MeshBase::set_spatial_dimension().

Referenced by main(), libMesh::Prism6::permute(), libMesh::Hex8::permute(), libMesh::Hex20::permute(), libMesh::Tet10::permute(), libMesh::Tet14::permute(), libMesh::Prism15::permute(), libMesh::Tet4::permute(), libMesh::Hex27::permute(), libMesh::Prism18::permute(), libMesh::Prism20::permute(), libMesh::Prism21::permute(), FETestBase< order, family, elem_type, N_x >::setUp(), BBoxTest::test_no_degenerate(), and BBoxTest::test_two_degenerate().

360 {
361  // We won't change any topology, but just changing geometry could
362  // invalidate a point locator.
363  mesh.clear_point_locator();
364 
365 #if LIBMESH_DIM == 3
366  const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
367 
368  if (theta)
369  mesh.set_spatial_dimension(3);
370 
371  for (auto & node : mesh.node_ptr_range())
372  {
373  Point & pt = *node;
374  pt = R * pt;
375  }
376 
377  return R;
378 
379 #else
380  libmesh_ignore(mesh, phi, theta, psi);
381  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
382  // We'll never get here
383  return RealTensorValue();
384 #endif
385 }
libMesh::MeshBase::set_spatial_dimension
void set_spatial_dimension(unsigned char d)
Sets the "spatial dimension" of the Mesh.
Definition: mesh_base.C:600
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::RealTensorValue
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:48
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:535
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1809
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ scale()

void libMesh::MeshTools::Modification::scale ( MeshBase &  mesh,
const Real  xs,
const Real  ys = 0.,
const Real  zs = 0. 
)

Scales the mesh.

The grid points are scaled in the x direction by xs, in the y direction by ys, etc... If only xs is specified then the scaling is assumed uniform in all directions.

Definition at line 388 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), mesh, and libMesh::Real.

Referenced by libMesh::PetscVector< libMesh::Number >::add(), main(), libMesh::DenseVector< Output >::operator*=(), libMesh::DenseMatrix< Real >::operator*=(), libMesh::PetscMatrix< T >::scale(), and libMesh::SparseMatrix< ValOut >::scale().

392 {
393  const Real x_scale = xs;
394  Real y_scale = ys;
395  Real z_scale = zs;
396 
397  if (ys == 0.)
398  {
399  libmesh_assert_equal_to (zs, 0.);
400 
401  y_scale = z_scale = x_scale;
402  }
403 
404  // Scale the x coordinate in all dimensions
405  for (auto & node : mesh.node_ptr_range())
406  (*node)(0) *= x_scale;
407 
408  // Only scale the y coordinate in 2 and 3D
409  if (LIBMESH_DIM < 2)
410  return;
411 
412  for (auto & node : mesh.node_ptr_range())
413  (*node)(1) *= y_scale;
414 
415  // Only scale the z coordinate in 3D
416  if (LIBMESH_DIM < 3)
417  return;
418 
419  for (auto & node : mesh.node_ptr_range())
420  (*node)(2) *= z_scale;
421 
422  // We haven't changed any topology, but just changing geometry could
423  // have invalidated a point locator.
424  mesh.clear_point_locator();
425 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1809
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144

◆ smooth()

void libMesh::MeshTools::Modification::smooth ( MeshBase &  mesh,
unsigned int  n_iterations,
Real  power 
)

Smooth the mesh with a simple Laplace smoothing algorithm.

The mesh is smoothed n_iterations times. If the parameter power is 0, each node is moved to the average position of the neighboring connected nodes. If power > 0, the node positions are weighted by their distance. The positions of higher order nodes, and nodes living in refined elements, are calculated from the vertex positions of their parent nodes. Only works in 2D.

Author
Martin Luthi (luthi.nosp@m.@gi..nosp@m.alask.nosp@m.a.ed.nosp@m.u)
Date
2005

This implementation assumes every element "side" has only 2 nodes.

Definition at line 1455 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::as_range(), libMesh::MeshBase::clear_point_locator(), libMesh::Elem::embedding_matrix(), libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), libMesh::make_range(), mesh, libMesh::MeshBase::mesh_dimension(), libMesh::MeshTools::n_levels(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_index_range(), libMesh::Elem::node_ptr(), libMesh::Elem::node_ref(), libMesh::MeshBase::node_ref(), libMesh::TypeVector< T >::norm(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::Utility::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), libMesh::Elem::side_index_range(), libMesh::MeshTools::weight(), and libMesh::Elem::which_child_am_i().

1458 {
1462  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1463 
1464  /*
1465  * Create a quickly-searchable list of boundary nodes.
1466  */
1467  std::unordered_set<dof_id_type> boundary_node_ids =
1468  MeshTools::find_boundary_nodes (mesh);
1469 
1470  // For avoiding extraneous element side allocation
1471  ElemSideBuilder side_builder;
1472 
1473  for (unsigned int iter=0; iter<n_iterations; iter++)
1474  {
1475  /*
1476  * loop over the mesh refinement level
1477  */
1478  unsigned int n_levels = MeshTools::n_levels(mesh);
1479  for (unsigned int refinement_level=0; refinement_level != n_levels;
1480  refinement_level++)
1481  {
1482  // initialize the storage (have to do it on every level to get empty vectors
1483  std::vector<Point> new_positions;
1484  std::vector<Real> weight;
1485  new_positions.resize(mesh.n_nodes());
1486  weight.resize(mesh.n_nodes());
1487 
1488  {
1489  // Loop over the elements to calculate new node positions
1490  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1491  mesh.level_elements_end(refinement_level)))
1492  {
1493  /*
1494  * We relax all nodes on level 0 first
1495  * If the element is refined (level > 0), we interpolate the
1496  * parents nodes with help of the embedding matrix
1497  */
1498  if (refinement_level == 0)
1499  {
1500  for (auto s : elem->side_index_range())
1501  {
1502  /*
1503  * Only operate on sides which are on the
1504  * boundary or for which the current element's
1505  * id is greater than its neighbor's.
1506  * Sides get only built once.
1507  */
1508  if ((elem->neighbor_ptr(s) != nullptr) &&
1509  (elem->id() > elem->neighbor_ptr(s)->id()))
1510  {
1511  const Elem & side = side_builder(*elem, s);
1512  const Node & node0 = side.node_ref(0);
1513  const Node & node1 = side.node_ref(1);
1514 
1515  Real node_weight = 1.;
1516  // calculate the weight of the nodes
1517  if (power > 0)
1518  {
1519  Point diff = node0-node1;
1520  node_weight = std::pow(diff.norm(), power);
1521  }
1522 
1523  const dof_id_type id0 = node0.id(), id1 = node1.id();
1524  new_positions[id0].add_scaled( node1, node_weight );
1525  new_positions[id1].add_scaled( node0, node_weight );
1526  weight[id0] += node_weight;
1527  weight[id1] += node_weight;
1528  }
1529  } // element neighbor loop
1530  }
1531 #ifdef LIBMESH_ENABLE_AMR
1532  else // refinement_level > 0
1533  {
1534  /*
1535  * Find the positions of the hanging nodes of refined elements.
1536  * We do this by calculating their position based on the parent
1537  * (one level less refined) element, and the embedding matrix
1538  */
1539 
1540  const Elem * parent = elem->parent();
1541 
1542  /*
1543  * find out which child I am
1544  */
1545  unsigned int c = parent->which_child_am_i(elem);
1546  /*
1547  *loop over the childs (that is, the current elements) nodes
1548  */
1549  for (auto nc : elem->node_index_range())
1550  {
1551  /*
1552  * the new position of the node
1553  */
1554  Point point;
1555  for (auto n : parent->node_index_range())
1556  {
1557  /*
1558  * The value from the embedding matrix
1559  */
1560  const Real em_val = parent->embedding_matrix(c,nc,n);
1561 
1562  if (em_val != 0.)
1563  point.add_scaled (parent->point(n), em_val);
1564  }
1565 
1566  const dof_id_type id = elem->node_ptr(nc)->id();
1567  new_positions[id] = point;
1568  weight[id] = 1.;
1569  }
1570  } // if element refinement_level
1571 #endif // #ifdef LIBMESH_ENABLE_AMR
1572 
1573  } // element loop
1574 
1575  /*
1576  * finally reposition the vertex nodes
1577  */
1578  for (auto nid : make_range(mesh.n_nodes()))
1579  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1580  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1581  }
1582 
1583  // Now handle the additional second_order nodes by calculating
1584  // their position based on the vertex positions
1585  // we do a second loop over the level elements
1586  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1587  mesh.level_elements_end(refinement_level)))
1588  {
1589  const unsigned int son_begin = elem->n_vertices();
1590  const unsigned int son_end = elem->n_nodes();
1591  for (unsigned int n=son_begin; n<son_end; n++)
1592  {
1593  const unsigned int n_adjacent_vertices =
1594  elem->n_second_order_adjacent_vertices(n);
1595 
1596  Point point;
1597  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1598  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1599 
1600  const dof_id_type id = elem->node_ptr(n)->id();
1601  mesh.node_ref(id) = point/n_adjacent_vertices;
1602  }
1603  }
1604  } // refinement_level loop
1605  } // end iteration
1606 
1607  // We haven't changed any topology, but just changing geometry could
1608  // have invalidated a point locator.
1609  mesh.clear_point_locator();
1610 }
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:3044
libMesh::TypeVector::add_scaled
void add_scaled(const TypeVector< T2 > &, const T &)
Add a scaled value to this vector without creating a temporary.
Definition: type_vector.h:629
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::TypeVector::norm
auto norm() const
Definition: type_vector.h:907
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::TypeVector::add
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:605
libMesh::MeshTools::find_boundary_nodes
std::unordered_set< dof_id_type > find_boundary_nodes(const MeshBase &mesh)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:524
libMesh::Elem::embedding_matrix
virtual Real embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:296
libMesh::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2535
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:819
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:444
libMesh::as_range
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
Helper function that allows us to treat a homogenous pair as a range.
Definition: simple_range.h:57
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:809
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1809
libMesh::ElemSideBuilder
Helper for building element sides that minimizes the construction of new elements.
Definition: elem_side_builder.h:50
libMesh::Elem::which_child_am_i
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:3206
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:173
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:422
libMesh::Elem::node_index_range
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2697
libMesh::MeshBase::node_ref
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:741
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:2459
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67

◆ translate()

void libMesh::MeshTools::Modification::translate ( MeshBase &  mesh,
const Real  xt = 0.,
const Real  yt = 0.,
const Real  zt = 0. 
)

Translates the mesh.

The grid points are translated in the x direction by xt, in the y direction by yt, etc...

Definition at line 317 of file mesh_modification.C.

References libMesh::MeshBase::clear_point_locator(), and mesh.

Referenced by MeshStitchTest::testMeshStitchElemsets(), and MeshTriangulationTest::testTriangulatorRoundHole().

321 {
322  const Point p(xt, yt, zt);
323 
324  for (auto & node : mesh.node_ptr_range())
325  *node += p;
326 
327  // We haven't changed any topology, but just changing geometry could
328  // have invalidated a point locator.
329  mesh.clear_point_locator();
330 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::MeshBase::clear_point_locator
void clear_point_locator()
Releases the current PointLocator object.
Definition: mesh_base.C:1809
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
Generated on Wed Jun 3 2026 21:24:58 for libMesh by   doxygen 1.8.14