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 409 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().

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

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

Referenced by MeshStitchTest::renameAndShift().

1722 {
1723  // This is just a shim around the member implementation, now
1724  mesh.get_boundary_info().renumber_id(old_id, new_id);
1725 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
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:1874

◆ 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 1729 of file mesh_modification.C.

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

1732 {
1733  if (old_id == new_id)
1734  {
1735  // If the IDs are the same, this is a no-op.
1736  return;
1737  }
1738 
1739  for (auto & elem : mesh.element_ptr_range())
1740  {
1741  if (elem->subdomain_id() == old_id)
1742  elem->subdomain_id() = new_id;
1743  }
1744 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476

◆ 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::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 }
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
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:517
libMesh::TypeVector::unit
TypeVector< T > unit() const
Definition: type_vector.h:1104
libMesh::MeshBase::query_node_ptr
virtual const Node * query_node_ptr(const dof_id_type i) const =0
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:140
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
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 1593 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().

1594 {
1595  libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
1596 
1597  // Algorithm:
1598  // .) For each active element in the mesh: construct a
1599  // copy which is the same in every way *except* it is
1600  // a level 0 element. Store the pointers to these in
1601  // a separate vector. Save any boundary information as well.
1602  // Delete the active element from the mesh.
1603  // .) Loop over all (remaining) elements in the mesh, delete them.
1604  // .) Add the level-0 copies back to the mesh
1605 
1606  // Temporary storage for new element pointers
1607  std::vector<std::unique_ptr<Elem>> new_elements;
1608 
1609  // BoundaryInfo Storage for element ids, sides, and BC ids
1610  std::vector<Elem *> saved_boundary_elements;
1611  std::vector<boundary_id_type> saved_bc_ids;
1612  std::vector<unsigned short int> saved_bc_sides;
1613 
1614  // Container to catch boundary ids passed back by BoundaryInfo
1615  std::vector<boundary_id_type> bc_ids;
1616 
1617  // Reserve a reasonable amt. of space for each
1618  new_elements.reserve(mesh.n_active_elem());
1619  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1620  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1621  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1622 
1623  for (auto & elem : mesh.active_element_ptr_range())
1624  {
1625  // Make a new element of the same type
1626  auto copy = Elem::build(elem->type());
1627 
1628  // Set node pointers (they still point to nodes in the original mesh)
1629  for (auto n : elem->node_index_range())
1630  copy->set_node(n, elem->node_ptr(n));
1631 
1632  // Copy over ids
1633  copy->processor_id() = elem->processor_id();
1634  copy->subdomain_id() = elem->subdomain_id();
1635 
1636  // Retain the original element's ID(s) as well, otherwise
1637  // the Mesh may try to create them for you...
1638  copy->set_id( elem->id() );
1639 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1640  copy->set_unique_id(elem->unique_id());
1641 #endif
1642 
1643  // This element could have boundary info or DistributedMesh
1644  // remote_elem links as well. We need to save the (elem,
1645  // side, bc_id) triples and those links
1646  for (auto s : elem->side_index_range())
1647  {
1648  if (elem->neighbor_ptr(s) == remote_elem)
1649  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1650 
1651  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1652  for (const auto & bc_id : bc_ids)
1653  if (bc_id != BoundaryInfo::invalid_id)
1654  {
1655  saved_boundary_elements.push_back(copy.get());
1656  saved_bc_ids.push_back(bc_id);
1657  saved_bc_sides.push_back(s);
1658  }
1659  }
1660 
1661  // Copy any extra element data. Since the copy hasn't been
1662  // added to the mesh yet any allocation has to be done manually.
1663  const unsigned int nei = elem->n_extra_integers();
1664  copy->add_extra_integers(nei);
1665  for (unsigned int i=0; i != nei; ++i)
1666  copy->set_extra_integer(i, elem->get_extra_integer(i));
1667 
1668  // Copy any mapping data.
1669  copy->set_mapping_type(elem->mapping_type());
1670  copy->set_mapping_data(elem->mapping_data());
1671 
1672  // We're done with this element
1673  mesh.delete_elem(elem);
1674 
1675  // But save the copy
1676  new_elements.push_back(std::move(copy));
1677  }
1678 
1679  // Make sure we saved the same number of boundary conditions
1680  // in each vector.
1681  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1682  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1683 
1684  // Loop again, delete any remaining elements
1685  for (auto & elem : mesh.element_ptr_range())
1686  mesh.delete_elem(elem);
1687 
1688  // Add the copied (now level-0) elements back to the mesh
1689  for (auto & new_elem : new_elements)
1690  {
1691  // Save the original ID, because the act of adding the Elem can
1692  // change new_elem's id!
1693  dof_id_type orig_id = new_elem->id();
1694 
1695  Elem * added_elem = mesh.add_elem(std::move(new_elem));
1696 
1697  // If the Elem, as it was re-added to the mesh, now has a
1698  // different ID (this is unlikely, so it's just an assert)
1699  // the boundary information will no longer be correct.
1700  libmesh_assert_equal_to (orig_id, added_elem->id());
1701 
1702  // Avoid compiler warnings in opt mode.
1703  libmesh_ignore(added_elem, orig_id);
1704  }
1705 
1706  // Finally, also add back the saved boundary information
1707  for (auto e : index_range(saved_boundary_elements))
1708  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1709  saved_bc_sides[e],
1710  saved_bc_ids[e]);
1711 
1712  // Trim unused and renumber nodes and elements
1713  mesh.prepare_for_use();
1714 }
libMesh::BoundaryInfo::n_boundary_conds
std::size_t n_boundary_conds() const
Definition: boundary_info.C:2240
libMesh::MeshBase::is_prepared
bool is_prepared() const
Definition: mesh_base.h:198
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 ecreated (or read) mesh for use.
Definition: mesh_base.C:759
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:1476
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:1333
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:515
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:828
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:233
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:1195
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
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 262 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().

263 {
264  LOG_SCOPE("orient_elements()", "MeshTools::Modification");
265 
266  // We don't yet support doing orient() on a parent element, which
267  // would require us to consistently orient all its children and
268  // give them different local child numbers.
269  unsigned int n_levels = MeshTools::n_levels(mesh);
270  if (n_levels > 1)
271  libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
272 
273  BoundaryInfo & boundary_info = mesh.get_boundary_info();
274  for (auto elem : mesh.element_ptr_range())
275  elem->orient(&boundary_info);
276 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:802
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 222 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, 1 >::setUp().

223 {
224  LOG_SCOPE("permute_elements()", "MeshTools::Modification");
225 
226  // We don't yet support doing permute() on a parent element, which
227  // would require us to consistently permute all its children and
228  // give them different local child numbers.
229  unsigned int n_levels = MeshTools::n_levels(mesh);
230  if (n_levels > 1)
231  libmesh_error();
232 
233  const unsigned int seed = 123456;
234 
235  // seed the random number generator.
236  // We'll loop from 1 to max_elem_id on every processor, even those
237  // that don't have a particular element, so that the pseudorandom
238  // numbers will be the same everywhere.
239  std::srand(seed);
240 
241 
242  for (auto e_id : make_range(mesh.max_elem_id()))
243  {
244  int my_rand = std::rand();
245 
246  Elem * elem = mesh.query_elem_ptr(e_id);
247 
248  if (!elem)
249  continue;
250 
251  const unsigned int max_permutation = elem->n_permutations();
252  if (!max_permutation)
253  continue;
254 
255  const unsigned int perm = my_rand % max_permutation;
256 
257  elem->permute(perm);
258  }
259 }
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:1476
libMesh::MeshTools::n_levels
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:802
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:140
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 280 of file mesh_modification.C.

References 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, 1 >::setUp(), and MeshSmootherTest::testSmoother().

282 {
283  libmesh_assert (mesh.n_nodes());
284  libmesh_assert (mesh.n_elem());
285 
286  LOG_SCOPE("redistribute()", "MeshTools::Modification");
287 
288  DenseVector<Real> output_vec(LIBMESH_DIM);
289 
290  // FIXME - we should thread this later.
291  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
292 
293  for (auto & node : mesh.node_ptr_range())
294  {
295  (*myfunc)(*node, output_vec);
296 
297  (*node)(0) = output_vec(0);
298 #if LIBMESH_DIM > 1
299  (*node)(1) = output_vec(1);
300 #endif
301 #if LIBMESH_DIM > 2
302  (*node)(2) = output_vec(2);
303 #endif
304  }
305 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
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 344 of file mesh_modification.C.

References 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, 1 >::setUp(), BBoxTest::test_no_degenerate(), and BBoxTest::test_two_degenerate().

348 {
349 #if LIBMESH_DIM == 3
350  const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
351 
352  if (theta)
353  mesh.set_spatial_dimension(3);
354 
355  for (auto & node : mesh.node_ptr_range())
356  {
357  Point & pt = *node;
358  pt = R * pt;
359  }
360 
361  return R;
362 
363 #else
364  libmesh_ignore(mesh, phi, theta, psi);
365  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
366  // We'll never get here
367  return RealTensorValue();
368 #endif
369 }
libMesh::MeshBase::set_spatial_dimension
void set_spatial_dimension(unsigned char d)
Sets the "spatial dimension" of the Mesh.
Definition: mesh_base.C:550
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
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:515
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 372 of file mesh_modification.C.

References 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().

376 {
377  const Real x_scale = xs;
378  Real y_scale = ys;
379  Real z_scale = zs;
380 
381  if (ys == 0.)
382  {
383  libmesh_assert_equal_to (zs, 0.);
384 
385  y_scale = z_scale = x_scale;
386  }
387 
388  // Scale the x coordinate in all dimensions
389  for (auto & node : mesh.node_ptr_range())
390  (*node)(0) *= x_scale;
391 
392  // Only scale the y coordinate in 2 and 3D
393  if (LIBMESH_DIM < 2)
394  return;
395 
396  for (auto & node : mesh.node_ptr_range())
397  (*node)(1) *= y_scale;
398 
399  // Only scale the z coordinate in 3D
400  if (LIBMESH_DIM < 3)
401  return;
402 
403  for (auto & node : mesh.node_ptr_range())
404  (*node)(2) *= z_scale;
405 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144

◆ 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 1437 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::as_range(), 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().

1440 {
1444  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1445 
1446  /*
1447  * Create a quickly-searchable list of boundary nodes.
1448  */
1449  std::unordered_set<dof_id_type> boundary_node_ids =
1450  MeshTools::find_boundary_nodes (mesh);
1451 
1452  // For avoiding extraneous element side allocation
1453  ElemSideBuilder side_builder;
1454 
1455  for (unsigned int iter=0; iter<n_iterations; iter++)
1456  {
1457  /*
1458  * loop over the mesh refinement level
1459  */
1460  unsigned int n_levels = MeshTools::n_levels(mesh);
1461  for (unsigned int refinement_level=0; refinement_level != n_levels;
1462  refinement_level++)
1463  {
1464  // initialize the storage (have to do it on every level to get empty vectors
1465  std::vector<Point> new_positions;
1466  std::vector<Real> weight;
1467  new_positions.resize(mesh.n_nodes());
1468  weight.resize(mesh.n_nodes());
1469 
1470  {
1471  // Loop over the elements to calculate new node positions
1472  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1473  mesh.level_elements_end(refinement_level)))
1474  {
1475  /*
1476  * We relax all nodes on level 0 first
1477  * If the element is refined (level > 0), we interpolate the
1478  * parents nodes with help of the embedding matrix
1479  */
1480  if (refinement_level == 0)
1481  {
1482  for (auto s : elem->side_index_range())
1483  {
1484  /*
1485  * Only operate on sides which are on the
1486  * boundary or for which the current element's
1487  * id is greater than its neighbor's.
1488  * Sides get only built once.
1489  */
1490  if ((elem->neighbor_ptr(s) != nullptr) &&
1491  (elem->id() > elem->neighbor_ptr(s)->id()))
1492  {
1493  const Elem & side = side_builder(*elem, s);
1494  const Node & node0 = side.node_ref(0);
1495  const Node & node1 = side.node_ref(1);
1496 
1497  Real node_weight = 1.;
1498  // calculate the weight of the nodes
1499  if (power > 0)
1500  {
1501  Point diff = node0-node1;
1502  node_weight = std::pow(diff.norm(), power);
1503  }
1504 
1505  const dof_id_type id0 = node0.id(), id1 = node1.id();
1506  new_positions[id0].add_scaled( node1, node_weight );
1507  new_positions[id1].add_scaled( node0, node_weight );
1508  weight[id0] += node_weight;
1509  weight[id1] += node_weight;
1510  }
1511  } // element neighbor loop
1512  }
1513 #ifdef LIBMESH_ENABLE_AMR
1514  else // refinement_level > 0
1515  {
1516  /*
1517  * Find the positions of the hanging nodes of refined elements.
1518  * We do this by calculating their position based on the parent
1519  * (one level less refined) element, and the embedding matrix
1520  */
1521 
1522  const Elem * parent = elem->parent();
1523 
1524  /*
1525  * find out which child I am
1526  */
1527  unsigned int c = parent->which_child_am_i(elem);
1528  /*
1529  *loop over the childs (that is, the current elements) nodes
1530  */
1531  for (auto nc : elem->node_index_range())
1532  {
1533  /*
1534  * the new position of the node
1535  */
1536  Point point;
1537  for (auto n : parent->node_index_range())
1538  {
1539  /*
1540  * The value from the embedding matrix
1541  */
1542  const Real em_val = parent->embedding_matrix(c,nc,n);
1543 
1544  if (em_val != 0.)
1545  point.add_scaled (parent->point(n), em_val);
1546  }
1547 
1548  const dof_id_type id = elem->node_ptr(nc)->id();
1549  new_positions[id] = point;
1550  weight[id] = 1.;
1551  }
1552  } // if element refinement_level
1553 #endif // #ifdef LIBMESH_ENABLE_AMR
1554 
1555  } // element loop
1556 
1557  /*
1558  * finally reposition the vertex nodes
1559  */
1560  for (auto nid : make_range(mesh.n_nodes()))
1561  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1562  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1563  }
1564 
1565  // Now handle the additional second_order nodes by calculating
1566  // their position based on the vertex positions
1567  // we do a second loop over the level elements
1568  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1569  mesh.level_elements_end(refinement_level)))
1570  {
1571  const unsigned int son_begin = elem->n_vertices();
1572  const unsigned int son_end = elem->n_nodes();
1573  for (unsigned int n=son_begin; n<son_end; n++)
1574  {
1575  const unsigned int n_adjacent_vertices =
1576  elem->n_second_order_adjacent_vertices(n);
1577 
1578  Point point;
1579  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1580  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1581 
1582  const dof_id_type id = elem->node_ptr(n)->id();
1583  mesh.node_ref(id) = point/n_adjacent_vertices;
1584  }
1585  }
1586  } // refinement_level loop
1587  } // end iteration
1588 }
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:3030
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 -> decltype(std::norm(T()))
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:1476
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:517
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:328
libMesh::Elem::node_ref
const Node & node_ref(const unsigned int i) const
Definition: elem.h:2529
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:828
libMesh::MeshTools::weight
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:437
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:802
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:3192
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:140
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
libMesh::Elem::node_index_range
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2683
libMesh::MeshBase::node_ref
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
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:2453
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 309 of file mesh_modification.C.

References mesh.

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

313 {
314  const Point p(xt, yt, zt);
315 
316  for (auto & node : mesh.node_ptr_range())
317  *node += p;
318 }
mesh
MeshBase & mesh
Definition: mesh_communication.C:1476
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
Generated on Wed Jul 30 2025 07:12:06 for libMesh by   doxygen 1.8.14