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 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...
 
void 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)
 Converts the 2D quadrilateral elements of a Mesh into triangular 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

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

Converts the 2D quadrilateral elements of a Mesh into triangular elements.

Note
Only works for 2D elements! 3D elements are ignored.
Probably won't do the right thing for meshes which have been refined previously.

Definition at line 740 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_side_ptr(), libMesh::ParallelObject::comm(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::MeshBase::element_ptr_range(), end, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::INFEDGE2, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::DofObject::invalid_id, libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libmesh_nullptr, libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::Parallel::Communicator::max(), libMesh::MeshBase::max_elem_id(), libMesh::MeshBase::mesh_dimension(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshBase::n_elem(), libMesh::Elem::neighbor_ptr(), libMesh::Elem::node_id(), libMesh::Elem::node_ptr(), libMesh::MeshBase::parallel_max_unique_id(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::point(), libMesh::MeshBase::prepare_for_use(), libMesh::PRISM18, libMesh::PRISM6, libMesh::ParallelObject::processor_id(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::remote_elem, libMesh::BoundaryInfo::remove(), libMesh::Elem::set_node(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::TET10, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by libMesh::MeshTools::Generation::build_sphere(), main(), AllTriTest::test_helper_2D(), and AllTriTest::test_helper_3D().

741 {
742  // The number of elements in the original mesh before any additions
743  // or deletions.
744  const dof_id_type n_orig_elem = mesh.n_elem();
745  const dof_id_type max_orig_id = mesh.max_elem_id();
746 
747  // We store pointers to the newly created elements in a vector
748  // until they are ready to be added to the mesh. This is because
749  // adding new elements on the fly can cause reallocation and invalidation
750  // of existing mesh element_iterators.
751  std::vector<Elem *> new_elements;
752 
753  unsigned int max_subelems = 1; // in 1D nothing needs to change
754  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
755  max_subelems = 2;
756  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
757  max_subelems = 6;
758 
759  new_elements.reserve (max_subelems*n_orig_elem);
760 
761  // If the original mesh has *side* boundary data, we carry that over
762  // to the new mesh with triangular elements. We currently only
763  // support bringing over side-based BCs to the all-tri mesh, but
764  // that could probably be extended to node and edge-based BCs as
765  // well.
766  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
767 
768  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
769  std::vector<Elem *> new_bndry_elements;
770  std::vector<unsigned short int> new_bndry_sides;
771  std::vector<boundary_id_type> new_bndry_ids;
772 
773  // We may need to add new points if we run into a 1.5th order
774  // element; if we do that on a DistributedMesh in a ghost element then
775  // we will need to fix their ids / unique_ids
776  bool added_new_ghost_point = false;
777 
778  // Iterate over the elements, splitting:
779  // QUADs into pairs of conforming triangles
780  // PYRAMIDs into pairs of conforming tets,
781  // PRISMs into triplets of conforming tets, and
782  // HEXs into quintets or sextets of conforming tets.
783  // We split on the shortest diagonal to give us better
784  // triangle quality in 2D, and we split based on node ids
785  // to guarantee consistency in 3D.
786 
787  // FIXME: This algorithm does not work on refined grids!
788  {
789 #ifdef LIBMESH_ENABLE_UNIQUE_ID
790  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
791 #endif
792 
793  for (auto & elem : mesh.element_ptr_range())
794  {
795  const ElemType etype = elem->type();
796 
797  // all_tri currently only works on coarse meshes
798  libmesh_assert (!elem->parent());
799 
800  // The new elements we will split the quad into.
801  // In 3D we may need as many as 6 tets per hex
802  Elem * subelem[6];
803 
804  for (unsigned int i = 0; i != max_subelems; ++i)
805  subelem[i] = libmesh_nullptr;
806 
807  switch (etype)
808  {
809  case QUAD4:
810  {
811  subelem[0] = new Tri3;
812  subelem[1] = new Tri3;
813 
814  // Check for possible edge swap
815  if ((elem->point(0) - elem->point(2)).norm() <
816  (elem->point(1) - elem->point(3)).norm())
817  {
818  subelem[0]->set_node(0) = elem->node_ptr(0);
819  subelem[0]->set_node(1) = elem->node_ptr(1);
820  subelem[0]->set_node(2) = elem->node_ptr(2);
821 
822  subelem[1]->set_node(0) = elem->node_ptr(0);
823  subelem[1]->set_node(1) = elem->node_ptr(2);
824  subelem[1]->set_node(2) = elem->node_ptr(3);
825  }
826 
827  else
828  {
829  subelem[0]->set_node(0) = elem->node_ptr(0);
830  subelem[0]->set_node(1) = elem->node_ptr(1);
831  subelem[0]->set_node(2) = elem->node_ptr(3);
832 
833  subelem[1]->set_node(0) = elem->node_ptr(1);
834  subelem[1]->set_node(1) = elem->node_ptr(2);
835  subelem[1]->set_node(2) = elem->node_ptr(3);
836  }
837 
838 
839  break;
840  }
841 
842  case QUAD8:
843  {
844  if (elem->processor_id() != mesh.processor_id())
845  added_new_ghost_point = true;
846 
847  subelem[0] = new Tri6;
848  subelem[1] = new Tri6;
849 
850  // Add a new node at the center (vertex average) of the element.
851  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
852  mesh.point(elem->node_id(1)) +
853  mesh.point(elem->node_id(2)) +
854  mesh.point(elem->node_id(3)))/4,
855  DofObject::invalid_id,
856  elem->processor_id());
857 
858  // Check for possible edge swap
859  if ((elem->point(0) - elem->point(2)).norm() <
860  (elem->point(1) - elem->point(3)).norm())
861  {
862  subelem[0]->set_node(0) = elem->node_ptr(0);
863  subelem[0]->set_node(1) = elem->node_ptr(1);
864  subelem[0]->set_node(2) = elem->node_ptr(2);
865  subelem[0]->set_node(3) = elem->node_ptr(4);
866  subelem[0]->set_node(4) = elem->node_ptr(5);
867  subelem[0]->set_node(5) = new_node;
868 
869  subelem[1]->set_node(0) = elem->node_ptr(0);
870  subelem[1]->set_node(1) = elem->node_ptr(2);
871  subelem[1]->set_node(2) = elem->node_ptr(3);
872  subelem[1]->set_node(3) = new_node;
873  subelem[1]->set_node(4) = elem->node_ptr(6);
874  subelem[1]->set_node(5) = elem->node_ptr(7);
875 
876  }
877 
878  else
879  {
880  subelem[0]->set_node(0) = elem->node_ptr(3);
881  subelem[0]->set_node(1) = elem->node_ptr(0);
882  subelem[0]->set_node(2) = elem->node_ptr(1);
883  subelem[0]->set_node(3) = elem->node_ptr(7);
884  subelem[0]->set_node(4) = elem->node_ptr(4);
885  subelem[0]->set_node(5) = new_node;
886 
887  subelem[1]->set_node(0) = elem->node_ptr(1);
888  subelem[1]->set_node(1) = elem->node_ptr(2);
889  subelem[1]->set_node(2) = elem->node_ptr(3);
890  subelem[1]->set_node(3) = elem->node_ptr(5);
891  subelem[1]->set_node(4) = elem->node_ptr(6);
892  subelem[1]->set_node(5) = new_node;
893  }
894 
895  break;
896  }
897 
898  case QUAD9:
899  {
900  subelem[0] = new Tri6;
901  subelem[1] = new Tri6;
902 
903  // Check for possible edge swap
904  if ((elem->point(0) - elem->point(2)).norm() <
905  (elem->point(1) - elem->point(3)).norm())
906  {
907  subelem[0]->set_node(0) = elem->node_ptr(0);
908  subelem[0]->set_node(1) = elem->node_ptr(1);
909  subelem[0]->set_node(2) = elem->node_ptr(2);
910  subelem[0]->set_node(3) = elem->node_ptr(4);
911  subelem[0]->set_node(4) = elem->node_ptr(5);
912  subelem[0]->set_node(5) = elem->node_ptr(8);
913 
914  subelem[1]->set_node(0) = elem->node_ptr(0);
915  subelem[1]->set_node(1) = elem->node_ptr(2);
916  subelem[1]->set_node(2) = elem->node_ptr(3);
917  subelem[1]->set_node(3) = elem->node_ptr(8);
918  subelem[1]->set_node(4) = elem->node_ptr(6);
919  subelem[1]->set_node(5) = elem->node_ptr(7);
920  }
921 
922  else
923  {
924  subelem[0]->set_node(0) = elem->node_ptr(0);
925  subelem[0]->set_node(1) = elem->node_ptr(1);
926  subelem[0]->set_node(2) = elem->node_ptr(3);
927  subelem[0]->set_node(3) = elem->node_ptr(4);
928  subelem[0]->set_node(4) = elem->node_ptr(8);
929  subelem[0]->set_node(5) = elem->node_ptr(7);
930 
931  subelem[1]->set_node(0) = elem->node_ptr(1);
932  subelem[1]->set_node(1) = elem->node_ptr(2);
933  subelem[1]->set_node(2) = elem->node_ptr(3);
934  subelem[1]->set_node(3) = elem->node_ptr(5);
935  subelem[1]->set_node(4) = elem->node_ptr(6);
936  subelem[1]->set_node(5) = elem->node_ptr(8);
937  }
938 
939  break;
940  }
941 
942  case PRISM6:
943  {
944  // Prisms all split into three tetrahedra
945  subelem[0] = new Tet4;
946  subelem[1] = new Tet4;
947  subelem[2] = new Tet4;
948 
949  // Triangular faces are not split.
950 
951  // On quad faces, we choose the node with the highest
952  // global id, and we split on the diagonal which
953  // includes that node. This ensures that (even in
954  // parallel, even on distributed meshes) the same
955  // diagonal split will be chosen for elements on either
956  // side of the same quad face. It also ensures that we
957  // always have a mix of "clockwise" and
958  // "counterclockwise" split faces (two of one and one
959  // of the other on each prism; this is useful since the
960  // alternative all-clockwise or all-counterclockwise
961  // face splittings can't be turned into tets without
962  // adding more nodes
963 
964  // Split on 0-4 diagonal
965  if (split_first_diagonal(elem, 0,4, 1,3))
966  {
967  // Split on 0-5 diagonal
968  if (split_first_diagonal(elem, 0,5, 2,3))
969  {
970  // Split on 1-5 diagonal
971  if (split_first_diagonal(elem, 1,5, 2,4))
972  {
973  subelem[0]->set_node(0) = elem->node_ptr(0);
974  subelem[0]->set_node(1) = elem->node_ptr(4);
975  subelem[0]->set_node(2) = elem->node_ptr(5);
976  subelem[0]->set_node(3) = elem->node_ptr(3);
977 
978  subelem[1]->set_node(0) = elem->node_ptr(0);
979  subelem[1]->set_node(1) = elem->node_ptr(4);
980  subelem[1]->set_node(2) = elem->node_ptr(1);
981  subelem[1]->set_node(3) = elem->node_ptr(5);
982 
983  subelem[2]->set_node(0) = elem->node_ptr(0);
984  subelem[2]->set_node(1) = elem->node_ptr(1);
985  subelem[2]->set_node(2) = elem->node_ptr(2);
986  subelem[2]->set_node(3) = elem->node_ptr(5);
987  }
988  else // Split on 2-4 diagonal
989  {
990  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
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[1]->set_node(0) = elem->node_ptr(0);
998  subelem[1]->set_node(1) = elem->node_ptr(4);
999  subelem[1]->set_node(2) = elem->node_ptr(2);
1000  subelem[1]->set_node(3) = elem->node_ptr(5);
1001 
1002  subelem[2]->set_node(0) = elem->node_ptr(0);
1003  subelem[2]->set_node(1) = elem->node_ptr(1);
1004  subelem[2]->set_node(2) = elem->node_ptr(2);
1005  subelem[2]->set_node(3) = elem->node_ptr(4);
1006  }
1007  }
1008  else // Split on 2-3 diagonal
1009  {
1010  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1011 
1012  // 0-4 and 2-3 split implies 2-4 split
1013  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1014 
1015  subelem[0]->set_node(0) = elem->node_ptr(0);
1016  subelem[0]->set_node(1) = elem->node_ptr(4);
1017  subelem[0]->set_node(2) = elem->node_ptr(2);
1018  subelem[0]->set_node(3) = elem->node_ptr(3);
1019 
1020  subelem[1]->set_node(0) = elem->node_ptr(3);
1021  subelem[1]->set_node(1) = elem->node_ptr(4);
1022  subelem[1]->set_node(2) = elem->node_ptr(2);
1023  subelem[1]->set_node(3) = elem->node_ptr(5);
1024 
1025  subelem[2]->set_node(0) = elem->node_ptr(0);
1026  subelem[2]->set_node(1) = elem->node_ptr(1);
1027  subelem[2]->set_node(2) = elem->node_ptr(2);
1028  subelem[2]->set_node(3) = elem->node_ptr(4);
1029  }
1030  }
1031  else // Split on 1-3 diagonal
1032  {
1033  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1034 
1035  // Split on 0-5 diagonal
1036  if (split_first_diagonal(elem, 0,5, 2,3))
1037  {
1038  // 1-3 and 0-5 split implies 1-5 split
1039  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1040 
1041  subelem[0]->set_node(0) = elem->node_ptr(1);
1042  subelem[0]->set_node(1) = elem->node_ptr(3);
1043  subelem[0]->set_node(2) = elem->node_ptr(4);
1044  subelem[0]->set_node(3) = elem->node_ptr(5);
1045 
1046  subelem[1]->set_node(0) = elem->node_ptr(1);
1047  subelem[1]->set_node(1) = elem->node_ptr(0);
1048  subelem[1]->set_node(2) = elem->node_ptr(3);
1049  subelem[1]->set_node(3) = elem->node_ptr(5);
1050 
1051  subelem[2]->set_node(0) = elem->node_ptr(0);
1052  subelem[2]->set_node(1) = elem->node_ptr(1);
1053  subelem[2]->set_node(2) = elem->node_ptr(2);
1054  subelem[2]->set_node(3) = elem->node_ptr(5);
1055  }
1056  else // Split on 2-3 diagonal
1057  {
1058  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1059 
1060  // Split on 1-5 diagonal
1061  if (split_first_diagonal(elem, 1,5, 2,4))
1062  {
1063  subelem[0]->set_node(0) = elem->node_ptr(0);
1064  subelem[0]->set_node(1) = elem->node_ptr(1);
1065  subelem[0]->set_node(2) = elem->node_ptr(2);
1066  subelem[0]->set_node(3) = elem->node_ptr(3);
1067 
1068  subelem[1]->set_node(0) = elem->node_ptr(3);
1069  subelem[1]->set_node(1) = elem->node_ptr(1);
1070  subelem[1]->set_node(2) = elem->node_ptr(2);
1071  subelem[1]->set_node(3) = elem->node_ptr(5);
1072 
1073  subelem[2]->set_node(0) = elem->node_ptr(1);
1074  subelem[2]->set_node(1) = elem->node_ptr(3);
1075  subelem[2]->set_node(2) = elem->node_ptr(4);
1076  subelem[2]->set_node(3) = elem->node_ptr(5);
1077  }
1078  else // Split on 2-4 diagonal
1079  {
1080  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1081 
1082  subelem[0]->set_node(0) = elem->node_ptr(0);
1083  subelem[0]->set_node(1) = elem->node_ptr(1);
1084  subelem[0]->set_node(2) = elem->node_ptr(2);
1085  subelem[0]->set_node(3) = elem->node_ptr(3);
1086 
1087  subelem[1]->set_node(0) = elem->node_ptr(2);
1088  subelem[1]->set_node(1) = elem->node_ptr(3);
1089  subelem[1]->set_node(2) = elem->node_ptr(4);
1090  subelem[1]->set_node(3) = elem->node_ptr(5);
1091 
1092  subelem[2]->set_node(0) = elem->node_ptr(3);
1093  subelem[2]->set_node(1) = elem->node_ptr(1);
1094  subelem[2]->set_node(2) = elem->node_ptr(2);
1095  subelem[2]->set_node(3) = elem->node_ptr(4);
1096  }
1097  }
1098  }
1099 
1100  break;
1101  }
1102 
1103  case PRISM18:
1104  {
1105  subelem[0] = new Tet10;
1106  subelem[1] = new Tet10;
1107  subelem[2] = new Tet10;
1108 
1109  // Split on 0-4 diagonal
1110  if (split_first_diagonal(elem, 0,4, 1,3))
1111  {
1112  // Split on 0-5 diagonal
1113  if (split_first_diagonal(elem, 0,5, 2,3))
1114  {
1115  // Split on 1-5 diagonal
1116  if (split_first_diagonal(elem, 1,5, 2,4))
1117  {
1118  subelem[0]->set_node(0) = elem->node_ptr(0);
1119  subelem[0]->set_node(1) = elem->node_ptr(4);
1120  subelem[0]->set_node(2) = elem->node_ptr(5);
1121  subelem[0]->set_node(3) = elem->node_ptr(3);
1122 
1123  subelem[0]->set_node(4) = elem->node_ptr(15);
1124  subelem[0]->set_node(5) = elem->node_ptr(13);
1125  subelem[0]->set_node(6) = elem->node_ptr(17);
1126  subelem[0]->set_node(7) = elem->node_ptr(9);
1127  subelem[0]->set_node(8) = elem->node_ptr(12);
1128  subelem[0]->set_node(9) = elem->node_ptr(14);
1129 
1130  subelem[1]->set_node(0) = elem->node_ptr(0);
1131  subelem[1]->set_node(1) = elem->node_ptr(4);
1132  subelem[1]->set_node(2) = elem->node_ptr(1);
1133  subelem[1]->set_node(3) = elem->node_ptr(5);
1134 
1135  subelem[1]->set_node(4) = elem->node_ptr(15);
1136  subelem[1]->set_node(5) = elem->node_ptr(10);
1137  subelem[1]->set_node(6) = elem->node_ptr(6);
1138  subelem[1]->set_node(7) = elem->node_ptr(17);
1139  subelem[1]->set_node(8) = elem->node_ptr(13);
1140  subelem[1]->set_node(9) = elem->node_ptr(16);
1141 
1142  subelem[2]->set_node(0) = elem->node_ptr(0);
1143  subelem[2]->set_node(1) = elem->node_ptr(1);
1144  subelem[2]->set_node(2) = elem->node_ptr(2);
1145  subelem[2]->set_node(3) = elem->node_ptr(5);
1146 
1147  subelem[2]->set_node(4) = elem->node_ptr(6);
1148  subelem[2]->set_node(5) = elem->node_ptr(7);
1149  subelem[2]->set_node(6) = elem->node_ptr(8);
1150  subelem[2]->set_node(7) = elem->node_ptr(17);
1151  subelem[2]->set_node(8) = elem->node_ptr(16);
1152  subelem[2]->set_node(9) = elem->node_ptr(11);
1153  }
1154  else // Split on 2-4 diagonal
1155  {
1156  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1157 
1158  subelem[0]->set_node(0) = elem->node_ptr(0);
1159  subelem[0]->set_node(1) = elem->node_ptr(4);
1160  subelem[0]->set_node(2) = elem->node_ptr(5);
1161  subelem[0]->set_node(3) = elem->node_ptr(3);
1162 
1163  subelem[0]->set_node(4) = elem->node_ptr(15);
1164  subelem[0]->set_node(5) = elem->node_ptr(13);
1165  subelem[0]->set_node(6) = elem->node_ptr(17);
1166  subelem[0]->set_node(7) = elem->node_ptr(9);
1167  subelem[0]->set_node(8) = elem->node_ptr(12);
1168  subelem[0]->set_node(9) = elem->node_ptr(14);
1169 
1170  subelem[1]->set_node(0) = elem->node_ptr(0);
1171  subelem[1]->set_node(1) = elem->node_ptr(4);
1172  subelem[1]->set_node(2) = elem->node_ptr(2);
1173  subelem[1]->set_node(3) = elem->node_ptr(5);
1174 
1175  subelem[1]->set_node(4) = elem->node_ptr(15);
1176  subelem[1]->set_node(5) = elem->node_ptr(16);
1177  subelem[1]->set_node(6) = elem->node_ptr(8);
1178  subelem[1]->set_node(7) = elem->node_ptr(17);
1179  subelem[1]->set_node(8) = elem->node_ptr(13);
1180  subelem[1]->set_node(9) = elem->node_ptr(11);
1181 
1182  subelem[2]->set_node(0) = elem->node_ptr(0);
1183  subelem[2]->set_node(1) = elem->node_ptr(1);
1184  subelem[2]->set_node(2) = elem->node_ptr(2);
1185  subelem[2]->set_node(3) = elem->node_ptr(4);
1186 
1187  subelem[2]->set_node(4) = elem->node_ptr(6);
1188  subelem[2]->set_node(5) = elem->node_ptr(7);
1189  subelem[2]->set_node(6) = elem->node_ptr(8);
1190  subelem[2]->set_node(7) = elem->node_ptr(15);
1191  subelem[2]->set_node(8) = elem->node_ptr(10);
1192  subelem[2]->set_node(9) = elem->node_ptr(16);
1193  }
1194  }
1195  else // Split on 2-3 diagonal
1196  {
1197  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1198 
1199  // 0-4 and 2-3 split implies 2-4 split
1200  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1201 
1202  subelem[0]->set_node(0) = elem->node_ptr(0);
1203  subelem[0]->set_node(1) = elem->node_ptr(4);
1204  subelem[0]->set_node(2) = elem->node_ptr(2);
1205  subelem[0]->set_node(3) = elem->node_ptr(3);
1206 
1207  subelem[0]->set_node(4) = elem->node_ptr(15);
1208  subelem[0]->set_node(5) = elem->node_ptr(16);
1209  subelem[0]->set_node(6) = elem->node_ptr(8);
1210  subelem[0]->set_node(7) = elem->node_ptr(9);
1211  subelem[0]->set_node(8) = elem->node_ptr(12);
1212  subelem[0]->set_node(9) = elem->node_ptr(17);
1213 
1214  subelem[1]->set_node(0) = elem->node_ptr(3);
1215  subelem[1]->set_node(1) = elem->node_ptr(4);
1216  subelem[1]->set_node(2) = elem->node_ptr(2);
1217  subelem[1]->set_node(3) = elem->node_ptr(5);
1218 
1219  subelem[1]->set_node(4) = elem->node_ptr(12);
1220  subelem[1]->set_node(5) = elem->node_ptr(16);
1221  subelem[1]->set_node(6) = elem->node_ptr(17);
1222  subelem[1]->set_node(7) = elem->node_ptr(14);
1223  subelem[1]->set_node(8) = elem->node_ptr(13);
1224  subelem[1]->set_node(9) = elem->node_ptr(11);
1225 
1226  subelem[2]->set_node(0) = elem->node_ptr(0);
1227  subelem[2]->set_node(1) = elem->node_ptr(1);
1228  subelem[2]->set_node(2) = elem->node_ptr(2);
1229  subelem[2]->set_node(3) = elem->node_ptr(4);
1230 
1231  subelem[2]->set_node(4) = elem->node_ptr(6);
1232  subelem[2]->set_node(5) = elem->node_ptr(7);
1233  subelem[2]->set_node(6) = elem->node_ptr(8);
1234  subelem[2]->set_node(7) = elem->node_ptr(15);
1235  subelem[2]->set_node(8) = elem->node_ptr(10);
1236  subelem[2]->set_node(9) = elem->node_ptr(16);
1237  }
1238  }
1239  else // Split on 1-3 diagonal
1240  {
1241  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1242 
1243  // Split on 0-5 diagonal
1244  if (split_first_diagonal(elem, 0,5, 2,3))
1245  {
1246  // 1-3 and 0-5 split implies 1-5 split
1247  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1248 
1249  subelem[0]->set_node(0) = elem->node_ptr(1);
1250  subelem[0]->set_node(1) = elem->node_ptr(3);
1251  subelem[0]->set_node(2) = elem->node_ptr(4);
1252  subelem[0]->set_node(3) = elem->node_ptr(5);
1253 
1254  subelem[0]->set_node(4) = elem->node_ptr(15);
1255  subelem[0]->set_node(5) = elem->node_ptr(12);
1256  subelem[0]->set_node(6) = elem->node_ptr(10);
1257  subelem[0]->set_node(7) = elem->node_ptr(16);
1258  subelem[0]->set_node(8) = elem->node_ptr(14);
1259  subelem[0]->set_node(9) = elem->node_ptr(13);
1260 
1261  subelem[1]->set_node(0) = elem->node_ptr(1);
1262  subelem[1]->set_node(1) = elem->node_ptr(0);
1263  subelem[1]->set_node(2) = elem->node_ptr(3);
1264  subelem[1]->set_node(3) = elem->node_ptr(5);
1265 
1266  subelem[1]->set_node(4) = elem->node_ptr(6);
1267  subelem[1]->set_node(5) = elem->node_ptr(9);
1268  subelem[1]->set_node(6) = elem->node_ptr(15);
1269  subelem[1]->set_node(7) = elem->node_ptr(16);
1270  subelem[1]->set_node(8) = elem->node_ptr(17);
1271  subelem[1]->set_node(9) = elem->node_ptr(14);
1272 
1273  subelem[2]->set_node(0) = elem->node_ptr(0);
1274  subelem[2]->set_node(1) = elem->node_ptr(1);
1275  subelem[2]->set_node(2) = elem->node_ptr(2);
1276  subelem[2]->set_node(3) = elem->node_ptr(5);
1277 
1278  subelem[2]->set_node(4) = elem->node_ptr(6);
1279  subelem[2]->set_node(5) = elem->node_ptr(7);
1280  subelem[2]->set_node(6) = elem->node_ptr(8);
1281  subelem[2]->set_node(7) = elem->node_ptr(17);
1282  subelem[2]->set_node(8) = elem->node_ptr(16);
1283  subelem[2]->set_node(9) = elem->node_ptr(11);
1284  }
1285  else // Split on 2-3 diagonal
1286  {
1287  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1288 
1289  // Split on 1-5 diagonal
1290  if (split_first_diagonal(elem, 1,5, 2,4))
1291  {
1292  subelem[0]->set_node(0) = elem->node_ptr(0);
1293  subelem[0]->set_node(1) = elem->node_ptr(1);
1294  subelem[0]->set_node(2) = elem->node_ptr(2);
1295  subelem[0]->set_node(3) = elem->node_ptr(3);
1296 
1297  subelem[0]->set_node(4) = elem->node_ptr(6);
1298  subelem[0]->set_node(5) = elem->node_ptr(7);
1299  subelem[0]->set_node(6) = elem->node_ptr(8);
1300  subelem[0]->set_node(7) = elem->node_ptr(9);
1301  subelem[0]->set_node(8) = elem->node_ptr(15);
1302  subelem[0]->set_node(9) = elem->node_ptr(17);
1303 
1304  subelem[1]->set_node(0) = elem->node_ptr(3);
1305  subelem[1]->set_node(1) = elem->node_ptr(1);
1306  subelem[1]->set_node(2) = elem->node_ptr(2);
1307  subelem[1]->set_node(3) = elem->node_ptr(5);
1308 
1309  subelem[1]->set_node(4) = elem->node_ptr(15);
1310  subelem[1]->set_node(5) = elem->node_ptr(7);
1311  subelem[1]->set_node(6) = elem->node_ptr(17);
1312  subelem[1]->set_node(7) = elem->node_ptr(14);
1313  subelem[1]->set_node(8) = elem->node_ptr(16);
1314  subelem[1]->set_node(9) = elem->node_ptr(11);
1315 
1316  subelem[2]->set_node(0) = elem->node_ptr(1);
1317  subelem[2]->set_node(1) = elem->node_ptr(3);
1318  subelem[2]->set_node(2) = elem->node_ptr(4);
1319  subelem[2]->set_node(3) = elem->node_ptr(5);
1320 
1321  subelem[2]->set_node(4) = elem->node_ptr(15);
1322  subelem[2]->set_node(5) = elem->node_ptr(12);
1323  subelem[2]->set_node(6) = elem->node_ptr(10);
1324  subelem[2]->set_node(7) = elem->node_ptr(16);
1325  subelem[2]->set_node(8) = elem->node_ptr(14);
1326  subelem[2]->set_node(9) = elem->node_ptr(13);
1327  }
1328  else // Split on 2-4 diagonal
1329  {
1330  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1331 
1332  subelem[0]->set_node(0) = elem->node_ptr(0);
1333  subelem[0]->set_node(1) = elem->node_ptr(1);
1334  subelem[0]->set_node(2) = elem->node_ptr(2);
1335  subelem[0]->set_node(3) = elem->node_ptr(3);
1336 
1337  subelem[0]->set_node(4) = elem->node_ptr(6);
1338  subelem[0]->set_node(5) = elem->node_ptr(7);
1339  subelem[0]->set_node(6) = elem->node_ptr(8);
1340  subelem[0]->set_node(7) = elem->node_ptr(9);
1341  subelem[0]->set_node(8) = elem->node_ptr(15);
1342  subelem[0]->set_node(9) = elem->node_ptr(17);
1343 
1344  subelem[1]->set_node(0) = elem->node_ptr(2);
1345  subelem[1]->set_node(1) = elem->node_ptr(3);
1346  subelem[1]->set_node(2) = elem->node_ptr(4);
1347  subelem[1]->set_node(3) = elem->node_ptr(5);
1348 
1349  subelem[1]->set_node(4) = elem->node_ptr(17);
1350  subelem[1]->set_node(5) = elem->node_ptr(12);
1351  subelem[1]->set_node(6) = elem->node_ptr(16);
1352  subelem[1]->set_node(7) = elem->node_ptr(11);
1353  subelem[1]->set_node(8) = elem->node_ptr(14);
1354  subelem[1]->set_node(9) = elem->node_ptr(13);
1355 
1356  subelem[2]->set_node(0) = elem->node_ptr(3);
1357  subelem[2]->set_node(1) = elem->node_ptr(1);
1358  subelem[2]->set_node(2) = elem->node_ptr(2);
1359  subelem[2]->set_node(3) = elem->node_ptr(4);
1360 
1361  subelem[2]->set_node(4) = elem->node_ptr(15);
1362  subelem[2]->set_node(5) = elem->node_ptr(7);
1363  subelem[2]->set_node(6) = elem->node_ptr(17);
1364  subelem[2]->set_node(7) = elem->node_ptr(12);
1365  subelem[2]->set_node(8) = elem->node_ptr(10);
1366  subelem[2]->set_node(9) = elem->node_ptr(16);
1367  }
1368  }
1369  }
1370 
1371  break;
1372  }
1373 
1374  // No need to split elements that are already simplicial:
1375  case EDGE2:
1376  case EDGE3:
1377  case EDGE4:
1378  case TRI3:
1379  case TRI6:
1380  case TET4:
1381  case TET10:
1382  case INFEDGE2:
1383  // No way to split infinite quad/prism elements, so
1384  // hopefully no need to
1385  case INFQUAD4:
1386  case INFQUAD6:
1387  case INFPRISM6:
1388  case INFPRISM12:
1389  continue;
1390  // If we're left with an unimplemented hex we're probably
1391  // out of luck. TODO: implement hexes
1392  default:
1393  {
1394  libMesh::err << "Error, encountered unimplemented element "
1395  << Utility::enum_to_string<ElemType>(etype)
1396  << " in MeshTools::Modification::all_tri()..."
1397  << std::endl;
1398  libmesh_not_implemented();
1399  }
1400  } // end switch (etype)
1401 
1402 
1403 
1404  // Be sure the correct IDs are also set for all subelems.
1405  for (unsigned int i=0; i != max_subelems; ++i)
1406  if (subelem[i]) {
1407  subelem[i]->processor_id() = elem->processor_id();
1408  subelem[i]->subdomain_id() = elem->subdomain_id();
1409  }
1410 
1411  // On a mesh with boundary data, we need to move that data to
1412  // the new elements.
1413 
1414  // On a mesh which is distributed, we need to move
1415  // remote_elem links to the new elements.
1416  bool mesh_is_serial = mesh.is_serial();
1417 
1418  if (mesh_has_boundary_data || mesh_is_serial)
1419  {
1420  // Container to key boundary IDs handed back by the BoundaryInfo object.
1421  std::vector<boundary_id_type> bc_ids;
1422 
1423  for (auto sn : elem->side_index_range())
1424  {
1425  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1426  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1427  {
1428  const boundary_id_type b_id = *id_it;
1429 
1430  if (mesh_is_serial && b_id == BoundaryInfo::invalid_id)
1431  continue;
1432 
1433  // Make a sorted list of node ids for elem->side(sn)
1434  UniquePtr<Elem> elem_side = elem->build_side_ptr(sn);
1435  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1436  for (std::size_t esn=0; esn<elem_side_nodes.size(); ++esn)
1437  elem_side_nodes[esn] = elem_side->node_id(esn);
1438  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1439 
1440  for (unsigned int i=0; i != max_subelems; ++i)
1441  if (subelem[i])
1442  {
1443  for (auto subside : subelem[i]->side_index_range())
1444  {
1445  UniquePtr<Elem> subside_elem = subelem[i]->build_side_ptr(subside);
1446 
1447  // Make a list of *vertex* node ids for this subside, see if they are all present
1448  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1449  // subelem[i]->key(subside) in the Prism cases, since the new side is
1450  // a different type. Note 2: we only use vertex nodes since, in the future,
1451  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1452  // original face will not contain the mid-edge node.
1453  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1454  for (std::size_t ssn=0; ssn<subside_nodes.size(); ++ssn)
1455  subside_nodes[ssn] = subside_elem->node_id(ssn);
1456  std::sort(subside_nodes.begin(), subside_nodes.end());
1457 
1458  // std::includes returns true if every element of the second sorted range is
1459  // contained in the first sorted range.
1460  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1461  subside_nodes.begin(), subside_nodes.end()))
1462  {
1463  if (b_id != BoundaryInfo::invalid_id)
1464  {
1465  new_bndry_ids.push_back(b_id);
1466  new_bndry_elements.push_back(subelem[i]);
1467  new_bndry_sides.push_back(subside);
1468  }
1469 
1470  // If the original element had a RemoteElem neighbor on side 'sn',
1471  // then the subelem has one on side 'subside'.
1472  if (elem->neighbor_ptr(sn) == remote_elem)
1473  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1474  }
1475  }
1476  }
1477  } // end for loop over boundary IDs
1478  } // end for loop over sides
1479 
1480  // Remove the original element from the BoundaryInfo structure.
1481  mesh.get_boundary_info().remove(elem);
1482 
1483  } // end if (mesh_has_boundary_data)
1484 
1485  // Determine new IDs for the split elements which will be
1486  // the same on all processors, therefore keeping the Mesh
1487  // in sync. Note: we offset the new IDs by max_orig_id to
1488  // avoid overwriting any of the original IDs.
1489  for (unsigned int i=0; i != max_subelems; ++i)
1490  if (subelem[i])
1491  {
1492  // Determine new IDs for the split elements which will be
1493  // the same on all processors, therefore keeping the Mesh
1494  // in sync. Note: we offset the new IDs by the max of the
1495  // pre-existing ids to avoid conflicting with originals.
1496  subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1497 
1498 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1499  subelem[i]->set_unique_id() = max_unique_id + 2*elem->unique_id() + i;
1500 #endif
1501 
1502  // Prepare to add the newly-created simplices
1503  new_elements.push_back(subelem[i]);
1504  }
1505 
1506  // Delete the original element
1507  mesh.delete_elem(elem);
1508  } // End for loop over elements
1509  } // end scope
1510 
1511 
1512  // Now, iterate over the new elements vector, and add them each to
1513  // the Mesh.
1514  {
1515  std::vector<Elem *>::iterator el = new_elements.begin();
1516  const std::vector<Elem *>::iterator end = new_elements.end();
1517  for (; el != end; ++el)
1518  mesh.add_elem(*el);
1519  }
1520 
1521  if (mesh_has_boundary_data)
1522  {
1523  // If the old mesh had boundary data, the new mesh better have
1524  // some. However, we can't assert that the size of
1525  // new_bndry_elements vector is > 0, since we may not have split
1526  // any elements actually on the boundary. We also can't assert
1527  // that the original number of boundary sides is equal to the
1528  // sum of the boundary sides currently in the mesh and the
1529  // newly-added boundary sides, since in 3D, we may have split a
1530  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1531  // too picky about the actual number of BCs, and just assert that
1532  // there are some, somewhere.
1533 #ifndef NDEBUG
1534  bool nbe_nonempty = new_bndry_elements.size();
1535  mesh.comm().max(nbe_nonempty);
1536  libmesh_assert(nbe_nonempty ||
1537  mesh.get_boundary_info().n_boundary_conds()>0);
1538 #endif
1539 
1540  // We should also be sure that the lengths of the new boundary data vectors
1541  // are all the same.
1542  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1543  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1544 
1545  // Add the new boundary info to the mesh
1546  for (std::size_t s=0; s<new_bndry_elements.size(); ++s)
1547  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1548  new_bndry_sides[s],
1549  new_bndry_ids[s]);
1550  }
1551 
1552  // In a DistributedMesh any newly added ghost node ids may be
1553  // inconsistent, and unique_ids of newly added ghost nodes remain
1554  // unset.
1555  // make_nodes_parallel_consistent() will fix all this.
1556  if (!mesh.is_serial())
1557  {
1558  mesh.comm().max(added_new_ghost_point);
1559 
1560  if (added_new_ghost_point)
1561  MeshCommunication().make_nodes_parallel_consistent (mesh);
1562  }
1563 
1564 
1565 
1566  // Prepare the newly created mesh for use.
1567  mesh.prepare_for_use(/*skip_renumber =*/ false);
1568 
1569  // Let the new_elements and new_bndry_elements vectors go out of scope.
1570 }
OStreamProxy err
ElemType
Defines an enum for geometric element types.
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
libmesh_assert(j)
int8_t boundary_id_type
Definition: id_types.h:51
uint8_t unique_id_type
Definition: id_types.h:79
uint8_t dof_id_type
Definition: id_types.h:64
const RemoteElem * remote_elem
Definition: remote_elem.C:57
void libMesh::MeshTools::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 1874 of file mesh_modification.C.

References libMesh::BoundaryInfo::add_edge(), libMesh::BoundaryInfo::add_node(), libMesh::BoundaryInfo::add_shellface(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::BoundaryInfo::build_edge_list(), libMesh::BoundaryInfo::build_node_list(), libMesh::BoundaryInfo::build_shellface_list(), libMesh::BoundaryInfo::build_side_list(), libMesh::BoundaryInfo::edge_boundary_ids(), libMesh::MeshBase::elem_ptr(), libMesh::MeshBase::get_boundary_info(), libMesh::MeshTools::Generation::Private::idx(), libMesh::MeshBase::node(), libMesh::MeshBase::node_ptr(), libMesh::BoundaryInfo::remove(), libMesh::BoundaryInfo::remove_edge(), libMesh::BoundaryInfo::remove_id(), libMesh::BoundaryInfo::remove_shellface(), libMesh::BoundaryInfo::remove_side(), libMesh::BoundaryInfo::shellface_boundary_ids(), and side.

1877 {
1878  if (old_id == new_id)
1879  {
1880  // If the IDs are the same, this is a no-op.
1881  return;
1882  }
1883 
1884  // A reference to the Mesh's BoundaryInfo object, for convenience.
1885  BoundaryInfo & bi = mesh.get_boundary_info();
1886 
1887  {
1888  // Build a list of all nodes that have boundary IDs
1889  std::vector<dof_id_type> node_list;
1890  std::vector<boundary_id_type> bc_id_list;
1891  bi.build_node_list (node_list, bc_id_list);
1892 
1893  // Temporary vector to hold ids
1894  std::vector<boundary_id_type> bndry_ids;
1895 
1896  // For each node with the old_id...
1897  for (std::size_t idx=0; idx<node_list.size(); ++idx)
1898  if (bc_id_list[idx] == old_id)
1899  {
1900  // Get the node in question
1901  const Node * node = mesh.node_ptr(node_list[idx]);
1902 
1903  // Get all the current IDs for this node.
1904  bi.boundary_ids(node, bndry_ids);
1905 
1906  // Update the IDs accordingly
1907  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1908 
1909  // Remove all traces of that node from the BoundaryInfo object
1910  bi.remove(node);
1911 
1912  // Add it back with the updated IDs
1913  bi.add_node(node, bndry_ids);
1914  }
1915  }
1916 
1917  {
1918  // Build a list of all edges that have boundary IDs
1919  std::vector<dof_id_type> elem_list;
1920  std::vector<unsigned short int> edge_list;
1921  std::vector<boundary_id_type> bc_id_list;
1922  bi.build_edge_list (elem_list, edge_list, bc_id_list);
1923 
1924  // Temporary vector to hold ids
1925  std::vector<boundary_id_type> bndry_ids;
1926 
1927  // For each edge with the old_id...
1928  for (std::size_t idx=0; idx<elem_list.size(); ++idx)
1929  if (bc_id_list[idx] == old_id)
1930  {
1931  // Get the elem in question
1932  const Elem * elem = mesh.elem_ptr(elem_list[idx]);
1933 
1934  // The edge of the elem in question
1935  unsigned short int edge = edge_list[idx];
1936 
1937  // Get all the current IDs for the edge in question.
1938  bi.edge_boundary_ids(elem, edge, bndry_ids);
1939 
1940  // Update the IDs accordingly
1941  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1942 
1943  // Remove all traces of that edge from the BoundaryInfo object
1944  bi.remove_edge(elem, edge);
1945 
1946  // Add it back with the updated IDs
1947  bi.add_edge(elem, edge, bndry_ids);
1948  }
1949  }
1950 
1951  {
1952  // Build a list of all shell-faces that have boundary IDs
1953  std::vector<dof_id_type> elem_list;
1954  std::vector<unsigned short int> shellface_list;
1955  std::vector<boundary_id_type> bc_id_list;
1956  bi.build_shellface_list (elem_list, shellface_list, bc_id_list);
1957 
1958  // Temporary vector to hold ids
1959  std::vector<boundary_id_type> bndry_ids;
1960 
1961  // For each shellface with the old_id...
1962  for (std::size_t idx=0; idx<elem_list.size(); ++idx)
1963  if (bc_id_list[idx] == old_id)
1964  {
1965  // Get the elem in question
1966  const Elem * elem = mesh.elem_ptr(elem_list[idx]);
1967 
1968  // The shellface of the elem in question
1969  unsigned short int shellface = shellface_list[idx];
1970 
1971  // Get all the current IDs for the shellface in question.
1972  bi.shellface_boundary_ids(elem, shellface, bndry_ids);
1973 
1974  // Update the IDs accordingly
1975  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1976 
1977  // Remove all traces of that shellface from the BoundaryInfo object
1978  bi.remove_shellface(elem, shellface);
1979 
1980  // Add it back with the updated IDs
1981  bi.add_shellface(elem, shellface, bndry_ids);
1982  }
1983  }
1984 
1985  {
1986  // Build a list of all sides that have boundary IDs
1987  std::vector<dof_id_type> elem_list;
1988  std::vector<unsigned short int> side_list;
1989  std::vector<boundary_id_type> bc_id_list;
1990  bi.build_side_list (elem_list, side_list, bc_id_list);
1991 
1992  // Temporary vector to hold ids
1993  std::vector<boundary_id_type> bndry_ids;
1994 
1995  // For each side with the old_id...
1996  for (std::size_t idx=0; idx<elem_list.size(); ++idx)
1997  if (bc_id_list[idx] == old_id)
1998  {
1999  // Get the elem in question
2000  const Elem * elem = mesh.elem_ptr(elem_list[idx]);
2001 
2002  // The side of the elem in question
2003  unsigned short int side = side_list[idx];
2004 
2005  // Get all the current IDs for the side in question.
2006  bi.boundary_ids(elem, side, bndry_ids);
2007 
2008  // Update the IDs accordingly
2009  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
2010 
2011  // Remove all traces of that side from the BoundaryInfo object
2012  bi.remove_side(elem, side);
2013 
2014  // Add it back with the updated IDs
2015  bi.add_side(elem, side, bndry_ids);
2016  }
2017  }
2018 
2019  // Remove any remaining references to the old_id so it does not show
2020  // up in lists of boundary ids, etc.
2021  bi.remove_id(old_id);
2022 }
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
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 2026 of file mesh_modification.C.

References libMesh::MeshBase::element_ptr_range(), and libMesh::Elem::subdomain_id().

2029 {
2030  if (old_id == new_id)
2031  {
2032  // If the IDs are the same, this is a no-op.
2033  return;
2034  }
2035 
2036  for (auto & elem : mesh.element_ptr_range())
2037  {
2038  if (elem->subdomain_id() == old_id)
2039  elem->subdomain_id() = new_id;
2040  }
2041 }
MeshBase & mesh
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 66 of file mesh_modification.C.

References libMesh::MeshBase::active_element_ptr_range(), libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::hmin(), libMesh::libmesh_assert(), std::max(), libMesh::MeshBase::max_node_id(), libMesh::MeshBase::mesh_dimension(), std::min(), 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().

69 {
70  libmesh_assert (mesh.n_nodes());
71  libmesh_assert (mesh.n_elem());
72  libmesh_assert ((factor >= 0.) && (factor <= 1.));
73 
74  LOG_SCOPE("distort()", "MeshTools::Modification");
75 
76  // First find nodes on the boundary and flag them
77  // so that we don't move them
78  // on_boundary holds false (not on boundary) and true (on boundary)
79  std::vector<bool> on_boundary (mesh.max_node_id(), false);
80 
81  if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);
82 
83  // Now calculate the minimum distance to
84  // neighboring nodes for each node.
85  // hmin holds these distances.
86  std::vector<float> hmin (mesh.max_node_id(),
88 
89  for (const auto & elem : mesh.active_element_ptr_range())
90  for (auto & n : elem->node_ref_range())
91  hmin[n.id()] = std::min(hmin[n.id()],
92  static_cast<float>(elem->hmin()));
93 
94 
95  // Now actually move the nodes
96  {
97  const unsigned int seed = 123456;
98 
99  // seed the random number generator.
100  // We'll loop from 1 to n_nodes on every processor, even those
101  // that don't have a particular node, so that the pseudorandom
102  // numbers will be the same everywhere.
103  std::srand(seed);
104 
105  // If the node is on the boundary or
106  // the node is not used by any element (hmin[n]<1.e20)
107  // then we should not move it.
108  // [Note: Testing for (in)equality might be wrong
109  // (different types, namely float and double)]
110  for (unsigned int n=0; n<mesh.max_node_id(); n++)
111  if (!on_boundary[n] && (hmin[n] < 1.e20) )
112  {
113  // the direction, random but unit normalized
114 
115  Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
116  (mesh.mesh_dimension() > 1) ?
117  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
118  : 0.,
119  ((mesh.mesh_dimension() == 3) ?
120  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
121  : 0.)
122  );
123 
124  dir(0) = (dir(0)-.5)*2.;
125  if (mesh.mesh_dimension() > 1)
126  dir(1) = (dir(1)-.5)*2.;
127  if (mesh.mesh_dimension() == 3)
128  dir(2) = (dir(2)-.5)*2.;
129 
130  dir = dir.unit();
131 
132  Node * node = mesh.query_node_ptr(n);
133  if (!node)
134  continue;
135 
136  (*node)(0) += dir(0)*factor*hmin[n];
137  if (mesh.mesh_dimension() > 1)
138  (*node)(1) += dir(1)*factor*hmin[n];
139  if (mesh.mesh_dimension() == 3)
140  (*node)(2) += dir(2)*factor*hmin[n];
141  }
142  }
143 }
MeshBase & mesh
long double max(long double a, double b)
libmesh_assert(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
long double min(long double a, double b)
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:290
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 1749 of file mesh_modification.C.

References libMesh::MeshBase::active_element_ptr_range(), libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_side(), bc_id, libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::element_ptr_range(), libMesh::MeshBase::get_boundary_info(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), 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::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::DofObject::set_unique_id(), libMesh::Elem::side_index_range(), libMesh::Elem::subdomain_id(), libMesh::Elem::type(), and libMesh::DofObject::unique_id().

Referenced by libMesh::MeshTools::Generation::build_sphere(), and main().

1750 {
1751  // Algorithm:
1752  // .) For each active element in the mesh: construct a
1753  // copy which is the same in every way *except* it is
1754  // a level 0 element. Store the pointers to these in
1755  // a separate vector. Save any boundary information as well.
1756  // Delete the active element from the mesh.
1757  // .) Loop over all (remaining) elements in the mesh, delete them.
1758  // .) Add the level-0 copies back to the mesh
1759 
1760  // Temporary storage for new element pointers
1761  std::vector<Elem *> new_elements;
1762 
1763  // BoundaryInfo Storage for element ids, sides, and BC ids
1764  std::vector<Elem *> saved_boundary_elements;
1765  std::vector<boundary_id_type> saved_bc_ids;
1766  std::vector<unsigned short int> saved_bc_sides;
1767 
1768  // Container to catch boundary ids passed back by BoundaryInfo
1769  std::vector<boundary_id_type> bc_ids;
1770 
1771  // Reserve a reasonable amt. of space for each
1772  new_elements.reserve(mesh.n_active_elem());
1773  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1774  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1775  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1776 
1777  for (auto & elem : mesh.active_element_ptr_range())
1778  {
1779  // Make a new element of the same type
1780  Elem * copy = Elem::build(elem->type()).release();
1781 
1782  // Set node pointers (they still point to nodes in the original mesh)
1783  for (auto n : elem->node_index_range())
1784  copy->set_node(n) = elem->node_ptr(n);
1785 
1786  // Copy over ids
1787  copy->processor_id() = elem->processor_id();
1788  copy->subdomain_id() = elem->subdomain_id();
1789 
1790  // Retain the original element's ID(s) as well, otherwise
1791  // the Mesh may try to create them for you...
1792  copy->set_id( elem->id() );
1793 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1794  copy->set_unique_id() = elem->unique_id();
1795 #endif
1796 
1797  // This element could have boundary info or DistributedMesh
1798  // remote_elem links as well. We need to save the (elem,
1799  // side, bc_id) triples and those links
1800  for (auto s : elem->side_index_range())
1801  {
1802  if (elem->neighbor_ptr(s) == remote_elem)
1803  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1804 
1805  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1806  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1807  {
1808  const boundary_id_type bc_id = *id_it;
1809 
1810  if (bc_id != BoundaryInfo::invalid_id)
1811  {
1812  saved_boundary_elements.push_back(copy);
1813  saved_bc_ids.push_back(bc_id);
1814  saved_bc_sides.push_back(s);
1815  }
1816  }
1817  }
1818 
1819 
1820  // We're done with this element
1821  mesh.delete_elem(elem);
1822 
1823  // But save the copy
1824  new_elements.push_back(copy);
1825  }
1826 
1827  // Make sure we saved the same number of boundary conditions
1828  // in each vector.
1829  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1830  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1831 
1832  // Loop again, delete any remaining elements
1833  for (auto & elem : mesh.element_ptr_range())
1834  mesh.delete_elem(elem);
1835 
1836  // Add the copied (now level-0) elements back to the mesh
1837  {
1838  for (std::vector<Elem *>::iterator it = new_elements.begin();
1839  it != new_elements.end();
1840  ++it)
1841  {
1842 #ifndef NDEBUG
1843  dof_id_type orig_id = (*it)->id();
1844 
1845  // ugly mid-statement endif to avoid unused variable warnings
1846  Elem * added_elem =
1847 #endif
1848  mesh.add_elem(*it);
1849 
1850 #ifndef NDEBUG
1851  dof_id_type added_id = added_elem->id();
1852 #endif
1853 
1854  // If the Elem, as it was re-added to the mesh, now has a
1855  // different ID (this is unlikely, so it's just an assert)
1856  // the boundary information will no longer be correct.
1857  libmesh_assert_equal_to (orig_id, added_id);
1858  }
1859  }
1860 
1861  // Finally, also add back the saved boundary information
1862  for (std::size_t e=0; e<saved_boundary_elements.size(); ++e)
1863  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1864  saved_bc_sides[e],
1865  saved_bc_ids[e]);
1866 
1867  // Trim unused and renumber nodes and elements
1868  mesh.prepare_for_use(/*skip_renumber =*/ false);
1869 }
MeshBase & mesh
boundary_id_type bc_id
Definition: xdr_io.C:50
int8_t boundary_id_type
Definition: id_types.h:51
uint8_t dof_id_type
Definition: id_types.h:64
const RemoteElem * remote_elem
Definition: remote_elem.C:57
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 147 of file mesh_modification.C.

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

Referenced by libMesh::MeshTools::Generation::build_cube().

149 {
150  libmesh_assert (mesh.n_nodes());
151  libmesh_assert (mesh.n_elem());
152 
153  LOG_SCOPE("redistribute()", "MeshTools::Modification");
154 
155  DenseVector<Real> output_vec(LIBMESH_DIM);
156 
157  // FIXME - we should thread this later.
158  UniquePtr<FunctionBase<Real>> myfunc = mapfunc.clone();
159 
160  for (auto & node : mesh.node_ptr_range())
161  {
162  (*myfunc)(*node, output_vec);
163 
164  (*node)(0) = output_vec(0);
165 #if LIBMESH_DIM > 1
166  (*node)(1) = output_vec(1);
167 #endif
168 #if LIBMESH_DIM > 2
169  (*node)(2) = output_vec(2);
170 #endif
171  }
172 }
MeshBase & mesh
libmesh_assert(j)
virtual UniquePtr< FunctionBase< Output > > clone() const =0
void 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)

Definition at line 210 of file mesh_modification.C.

References libMesh::MeshBase::node_ptr_range(), libMesh::pi, libMesh::Real, and libMesh::x.

214 {
215 #if LIBMESH_DIM == 3
216  const Real p = -phi/180.*libMesh::pi;
217  const Real t = -theta/180.*libMesh::pi;
218  const Real s = -psi/180.*libMesh::pi;
219  const Real sp = std::sin(p), cp = std::cos(p);
220  const Real st = std::sin(t), ct = std::cos(t);
221  const Real ss = std::sin(s), cs = std::cos(s);
222 
223  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
224  // (equations 6-14 give the entries of the composite transformation matrix).
225  // The rotations are performed sequentially about the z, x, and z axes, in that order.
226  // A positive angle yields a counter-clockwise rotation about the axis in question.
227  for (auto & node : mesh.node_ptr_range())
228  {
229  const Point pt = *node;
230  const Real x = pt(0);
231  const Real y = pt(1);
232  const Real z = pt(2);
233  *node = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
234  (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
235  ( sp*st)*x + (-cp*st)*y + (ct)*z );
236  }
237 #else
238  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
239 #endif
240 }
MeshBase & mesh
PetscErrorCode Vec x
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Real pi
.
Definition: libmesh.h:172
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 243 of file mesh_modification.C.

References libMesh::MeshBase::node_ptr_range(), and libMesh::Real.

Referenced by libMesh::libmesh_final< T >::add_vector_transpose(), and main().

247 {
248  const Real x_scale = xs;
249  Real y_scale = ys;
250  Real z_scale = zs;
251 
252  if (ys == 0.)
253  {
254  libmesh_assert_equal_to (zs, 0.);
255 
256  y_scale = z_scale = x_scale;
257  }
258 
259  // Scale the x coordinate in all dimensions
260  for (auto & node : mesh.node_ptr_range())
261  (*node)(0) *= x_scale;
262 
263  // Only scale the y coordinate in 2 and 3D
264  if (LIBMESH_DIM < 2)
265  return;
266 
267  for (auto & node : mesh.node_ptr_range())
268  (*node)(1) *= y_scale;
269 
270  // Only scale the z coordinate in 3D
271  if (LIBMESH_DIM < 3)
272  return;
273 
274  for (auto & node : mesh.node_ptr_range())
275  (*node)(2) *= z_scale;
276 }
MeshBase & mesh
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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 1573 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::build_side_ptr(), libMesh::Elem::embedding_matrix(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::DofObject::id(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), libmesh_nullptr, 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::MeshBase::node_ref(), libMesh::TypeVector< T >::norm(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::point(), std::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), side, libMesh::Elem::side_index_range(), libMesh::MeshTools::weight(), and libMesh::Elem::which_child_am_i().

1576 {
1580  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1581 
1582  /*
1583  * find the boundary nodes
1584  */
1585  std::vector<bool> on_boundary;
1586  MeshTools::find_boundary_nodes(mesh, on_boundary);
1587 
1588  for (unsigned int iter=0; iter<n_iterations; iter++)
1589 
1590  {
1591  /*
1592  * loop over the mesh refinement level
1593  */
1594  unsigned int n_levels = MeshTools::n_levels(mesh);
1595  for (unsigned int refinement_level=0; refinement_level != n_levels;
1596  refinement_level++)
1597  {
1598  // initialize the storage (have to do it on every level to get empty vectors
1599  std::vector<Point> new_positions;
1600  std::vector<Real> weight;
1601  new_positions.resize(mesh.n_nodes());
1602  weight.resize(mesh.n_nodes());
1603 
1604  {
1605  /*
1606  * Loop over the elements to calculate new node positions
1607  */
1608  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1609  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1610 
1611  for (; el != end; ++el)
1612  {
1613  /*
1614  * Constant handle for the element
1615  */
1616  const Elem * elem = *el;
1617 
1618  /*
1619  * We relax all nodes on level 0 first
1620  * If the element is refined (level > 0), we interpolate the
1621  * parents nodes with help of the embedding matrix
1622  */
1623  if (refinement_level == 0)
1624  {
1625  for (auto s : elem->side_index_range())
1626  {
1627  /*
1628  * Only operate on sides which are on the
1629  * boundary or for which the current element's
1630  * id is greater than its neighbor's.
1631  * Sides get only built once.
1632  */
1633  if ((elem->neighbor_ptr(s) != libmesh_nullptr) &&
1634  (elem->id() > elem->neighbor_ptr(s)->id()))
1635  {
1636  UniquePtr<const Elem> side(elem->build_side_ptr(s));
1637 
1638  const Node & node0 = side->node_ref(0);
1639  const Node & node1 = side->node_ref(1);
1640 
1641  Real node_weight = 1.;
1642  // calculate the weight of the nodes
1643  if (power > 0)
1644  {
1645  Point diff = node0-node1;
1646  node_weight = std::pow(diff.norm(), power);
1647  }
1648 
1649  const dof_id_type id0 = node0.id(), id1 = node1.id();
1650  new_positions[id0].add_scaled( node1, node_weight );
1651  new_positions[id1].add_scaled( node0, node_weight );
1652  weight[id0] += node_weight;
1653  weight[id1] += node_weight;
1654  }
1655  } // element neighbor loop
1656  }
1657 #ifdef LIBMESH_ENABLE_AMR
1658  else // refinement_level > 0
1659  {
1660  /*
1661  * Find the positions of the hanging nodes of refined elements.
1662  * We do this by calculating their position based on the parent
1663  * (one level less refined) element, and the embedding matrix
1664  */
1665 
1666  const Elem * parent = elem->parent();
1667 
1668  /*
1669  * find out which child I am
1670  */
1671  unsigned int c = parent->which_child_am_i(elem);
1672  /*
1673  *loop over the childs (that is, the current elements) nodes
1674  */
1675  for (auto nc : elem->node_index_range())
1676  {
1677  /*
1678  * the new position of the node
1679  */
1680  Point point;
1681  for (auto n : parent->node_index_range())
1682  {
1683  /*
1684  * The value from the embedding matrix
1685  */
1686  const float em_val = parent->embedding_matrix(c,nc,n);
1687 
1688  if (em_val != 0.)
1689  point.add_scaled (parent->point(n), em_val);
1690  }
1691 
1692  const dof_id_type id = elem->node_ptr(nc)->id();
1693  new_positions[id] = point;
1694  weight[id] = 1.;
1695  }
1696  } // if element refinement_level
1697 #endif // #ifdef LIBMESH_ENABLE_AMR
1698 
1699  } // element loop
1700 
1701  /*
1702  * finally reposition the vertex nodes
1703  */
1704  for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
1705  if (!on_boundary[nid] && weight[nid] > 0.)
1706  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1707  }
1708 
1709  {
1710  /*
1711  * Now handle the additional second_order nodes by calculating
1712  * their position based on the vertex positions
1713  * we do a second loop over the level elements
1714  */
1715  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1716  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1717 
1718  for (; el != end; ++el)
1719  {
1720  /*
1721  * Constant handle for the element
1722  */
1723  const Elem * elem = *el;
1724  const unsigned int son_begin = elem->n_vertices();
1725  const unsigned int son_end = elem->n_nodes();
1726  for (unsigned int n=son_begin; n<son_end; n++)
1727  {
1728  const unsigned int n_adjacent_vertices =
1729  elem->n_second_order_adjacent_vertices(n);
1730 
1731  Point point;
1732  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1733  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1734 
1735  const dof_id_type id = elem->node_ptr(n)->id();
1736  mesh.node_ref(id) = point/n_adjacent_vertices;
1737  }
1738  }
1739  }
1740 
1741  } // refinement_level loop
1742 
1743  } // end iteration
1744 }
unsigned short int side
Definition: xdr_io.C:49
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:245
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:603
double pow(double a, int b)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Calling this function on a 2D mesh will convert all the elements to triangles.
Definition: mesh_tools.C:290
uint8_t dof_id_type
Definition: id_types.h:64
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 176 of file mesh_modification.C.

References libMesh::MeshBase::node_ptr_range().

180 {
181  const Point p(xt, yt, zt);
182 
183  for (auto & node : mesh.node_ptr_range())
184  *node += p;
185 }
MeshBase & mesh