libMesh::MeshTools::Modification Namespace Reference

Functions

void distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
 
void translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
 
void rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
 
void scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
 
void all_tri (MeshBase &mesh)
 
void smooth (MeshBase &, unsigned int, Real)
 
void flatten (MeshBase &mesh)
 
void change_boundary_id (MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id)
 
void change_subdomain_id (MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)
 

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. Note: Probably won't do the right thing for meshes which have been refined previously.

Definition at line 686 of file mesh_modification.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::MeshBase::boundary_info, libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::err, libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_assert(), libMesh::libmesh_assert_greater(), libMesh::MeshBase::n_elem(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::Elem::node(), libMesh::MeshBase::node(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMesh::remote_elem, libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), libMeshEnums::TRI3, libMeshEnums::TRI6, and libMesh::Elem::type().

687 {
688  // The number of elements in the original mesh before any additions
689  // or deletions.
690  const dof_id_type n_orig_elem = mesh.n_elem();
691 
692  // We store pointers to the newly created elements in a vector
693  // until they are ready to be added to the mesh. This is because
694  // adding new elements on the fly can cause reallocation and invalidation
695  // of existing iterators.
696  std::vector<Elem*> new_elements;
697  new_elements.reserve (2*n_orig_elem);
698 
699  // If the original mesh has boundary data, we carry that over
700  // to the new mesh with triangular elements.
701  const bool mesh_has_boundary_data = (mesh.boundary_info->n_boundary_ids() > 0);
702 
703  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
704  std::vector<Elem*> new_bndry_elements;
705  std::vector<unsigned short int> new_bndry_sides;
706  std::vector<boundary_id_type> new_bndry_ids;
707 
708  // Iterate over the elements, splitting QUADS into
709  // pairs of conforming triangles.
710  // FIXME: This algorithm does not work on refined grids!
711  {
712  MeshBase::element_iterator el = mesh.elements_begin();
713  const MeshBase::element_iterator end = mesh.elements_end();
714 
715  for (; el!=end; ++el)
716  {
717  Elem* elem = *el;
718 
719  const ElemType etype = elem->type();
720 
721  // all_tri currently only works on coarse meshes
722  libmesh_assert (!elem->parent());
723 
724  // We split the quads using the shorter of the two diagonals
725  // to maintain the best angle properties.
726  bool edge_swap = false;
727 
728  // True if we actually split the current element.
729  bool split_elem = false;
730 
731  // The two new triangular elements we will split the quad into.
732  Elem* tri0 = NULL;
733  Elem* tri1 = NULL;
734 
735 
736  switch (etype)
737  {
738  case QUAD4:
739  {
740  split_elem = true;
741 
742  tri0 = new Tri3;
743  tri1 = new Tri3;
744 
745  // Check for possible edge swap
746  if ((elem->point(0) - elem->point(2)).size() <
747  (elem->point(1) - elem->point(3)).size())
748  {
749  tri0->set_node(0) = elem->get_node(0);
750  tri0->set_node(1) = elem->get_node(1);
751  tri0->set_node(2) = elem->get_node(2);
752 
753  tri1->set_node(0) = elem->get_node(0);
754  tri1->set_node(1) = elem->get_node(2);
755  tri1->set_node(2) = elem->get_node(3);
756  }
757 
758  else
759  {
760  edge_swap=true;
761 
762  tri0->set_node(0) = elem->get_node(0);
763  tri0->set_node(1) = elem->get_node(1);
764  tri0->set_node(2) = elem->get_node(3);
765 
766  tri1->set_node(0) = elem->get_node(1);
767  tri1->set_node(1) = elem->get_node(2);
768  tri1->set_node(2) = elem->get_node(3);
769  }
770 
771 
772  break;
773  }
774 
775  case QUAD8:
776  {
777  split_elem = true;
778 
779  tri0 = new Tri6;
780  tri1 = new Tri6;
781 
782  Node* new_node = mesh.add_point( (mesh.node(elem->node(0)) +
783  mesh.node(elem->node(1)) +
784  mesh.node(elem->node(2)) +
785  mesh.node(elem->node(3)) / 4)
786  );
787 
788  // Check for possible edge swap
789  if ((elem->point(0) - elem->point(2)).size() <
790  (elem->point(1) - elem->point(3)).size())
791  {
792  tri0->set_node(0) = elem->get_node(0);
793  tri0->set_node(1) = elem->get_node(1);
794  tri0->set_node(2) = elem->get_node(2);
795  tri0->set_node(3) = elem->get_node(4);
796  tri0->set_node(4) = elem->get_node(5);
797  tri0->set_node(5) = new_node;
798 
799  tri1->set_node(0) = elem->get_node(0);
800  tri1->set_node(1) = elem->get_node(2);
801  tri1->set_node(2) = elem->get_node(3);
802  tri1->set_node(3) = new_node;
803  tri1->set_node(4) = elem->get_node(6);
804  tri1->set_node(5) = elem->get_node(7);
805 
806  }
807 
808  else
809  {
810  edge_swap=true;
811 
812  tri0->set_node(0) = elem->get_node(3);
813  tri0->set_node(1) = elem->get_node(0);
814  tri0->set_node(2) = elem->get_node(1);
815  tri0->set_node(3) = elem->get_node(7);
816  tri0->set_node(4) = elem->get_node(4);
817  tri0->set_node(5) = new_node;
818 
819  tri1->set_node(0) = elem->get_node(1);
820  tri1->set_node(1) = elem->get_node(2);
821  tri1->set_node(2) = elem->get_node(3);
822  tri1->set_node(3) = elem->get_node(5);
823  tri1->set_node(4) = elem->get_node(6);
824  tri1->set_node(5) = new_node;
825  }
826 
827  break;
828  }
829 
830  case QUAD9:
831  {
832  split_elem = true;
833 
834  tri0 = new Tri6;
835  tri1 = new Tri6;
836 
837  // Check for possible edge swap
838  if ((elem->point(0) - elem->point(2)).size() <
839  (elem->point(1) - elem->point(3)).size())
840  {
841  tri0->set_node(0) = elem->get_node(0);
842  tri0->set_node(1) = elem->get_node(1);
843  tri0->set_node(2) = elem->get_node(2);
844  tri0->set_node(3) = elem->get_node(4);
845  tri0->set_node(4) = elem->get_node(5);
846  tri0->set_node(5) = elem->get_node(8);
847 
848  tri1->set_node(0) = elem->get_node(0);
849  tri1->set_node(1) = elem->get_node(2);
850  tri1->set_node(2) = elem->get_node(3);
851  tri1->set_node(3) = elem->get_node(8);
852  tri1->set_node(4) = elem->get_node(6);
853  tri1->set_node(5) = elem->get_node(7);
854  }
855 
856  else
857  {
858  edge_swap=true;
859 
860  tri0->set_node(0) = elem->get_node(0);
861  tri0->set_node(1) = elem->get_node(1);
862  tri0->set_node(2) = elem->get_node(3);
863  tri0->set_node(3) = elem->get_node(4);
864  tri0->set_node(4) = elem->get_node(8);
865  tri0->set_node(5) = elem->get_node(7);
866 
867  tri1->set_node(0) = elem->get_node(1);
868  tri1->set_node(1) = elem->get_node(2);
869  tri1->set_node(2) = elem->get_node(3);
870  tri1->set_node(3) = elem->get_node(5);
871  tri1->set_node(4) = elem->get_node(6);
872  tri1->set_node(5) = elem->get_node(8);
873  }
874 
875  break;
876  }
877  // No need to split elements that are already triangles
878  case TRI3:
879  case TRI6:
880  continue;
881  // Try to ignore non-2D elements for now
882  default:
883  {
884  libMesh::err << "Warning, encountered non-2D element "
885  << Utility::enum_to_string<ElemType>(etype)
886  << " in MeshTools::Modification::all_tri(), hope that's OK..."
887  << std::endl;
888  }
889  } // end switch (etype)
890 
891 
892 
893  if (split_elem)
894  {
895  // Be sure the correct ID's are also set for tri0 and
896  // tri1.
897  tri0->processor_id() = elem->processor_id();
898  tri0->subdomain_id() = elem->subdomain_id();
899  tri1->processor_id() = elem->processor_id();
900  tri1->subdomain_id() = elem->subdomain_id();
901 
902  if (mesh_has_boundary_data)
903  {
904  for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
905  {
906  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids(*el, sn);
907  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
908  {
909  const boundary_id_type b_id = *id_it;
910 
911  if (b_id != BoundaryInfo::invalid_id)
912  {
913  // Add the boundary ID to the list of new boundary ids
914  new_bndry_ids.push_back(b_id);
915 
916  // Convert the boundary side information of the old element to
917  // boundary side information for the new element.
918  if (!edge_swap)
919  {
920  switch (sn)
921  {
922  case 0:
923  {
924  // New boundary side is Tri 0, side 0
925  new_bndry_elements.push_back(tri0);
926  new_bndry_sides.push_back(0);
927  break;
928  }
929  case 1:
930  {
931  // New boundary side is Tri 0, side 1
932  new_bndry_elements.push_back(tri0);
933  new_bndry_sides.push_back(1);
934  break;
935  }
936  case 2:
937  {
938  // New boundary side is Tri 1, side 1
939  new_bndry_elements.push_back(tri1);
940  new_bndry_sides.push_back(1);
941  break;
942  }
943  case 3:
944  {
945  // New boundary side is Tri 1, side 2
946  new_bndry_elements.push_back(tri1);
947  new_bndry_sides.push_back(2);
948  break;
949  }
950 
951  default:
952  {
953  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
954  libmesh_error();
955  }
956  }
957  }
958 
959  else // edge_swap==true
960  {
961  switch (sn)
962  {
963  case 0:
964  {
965  // New boundary side is Tri 0, side 0
966  new_bndry_elements.push_back(tri0);
967  new_bndry_sides.push_back(0);
968  break;
969  }
970  case 1:
971  {
972  // New boundary side is Tri 1, side 0
973  new_bndry_elements.push_back(tri1);
974  new_bndry_sides.push_back(0);
975  break;
976  }
977  case 2:
978  {
979  // New boundary side is Tri 1, side 1
980  new_bndry_elements.push_back(tri1);
981  new_bndry_sides.push_back(1);
982  break;
983  }
984  case 3:
985  {
986  // New boundary side is Tri 0, side 2
987  new_bndry_elements.push_back(tri0);
988  new_bndry_sides.push_back(2);
989  break;
990  }
991 
992  default:
993  {
994  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
995  libmesh_error();
996  }
997  }
998  } // end edge_swap==true
999  } // end if (b_id != BoundaryInfo::invalid_id)
1000  } // end for loop over boundary IDs
1001  } // end for loop over sides
1002 
1003  // Remove the original element from the BoundaryInfo structure.
1004  mesh.boundary_info->remove(elem);
1005 
1006  } // end if (mesh_has_boundary_data)
1007 
1008 
1009  // On a distributed mesh, we need to preserve remote_elem
1010  // links, since prepare_for_use can't reconstruct them for
1011  // us.
1012  for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
1013  {
1014  if (elem->neighbor(sn) == remote_elem)
1015  {
1016  // Create a remote_elem link on one of the new
1017  // elements corresponding to the link from the old
1018  // element.
1019  if (!edge_swap)
1020  {
1021  switch (sn)
1022  {
1023  case 0:
1024  {
1025  // New remote side is Tri 0, side 0
1026  tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
1027  break;
1028  }
1029  case 1:
1030  {
1031  // New remote side is Tri 0, side 1
1032  tri0->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
1033  break;
1034  }
1035  case 2:
1036  {
1037  // New remote side is Tri 1, side 1
1038  tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
1039  break;
1040  }
1041  case 3:
1042  {
1043  // New remote side is Tri 1, side 2
1044  tri1->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
1045  break;
1046  }
1047 
1048  default:
1049  {
1050  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
1051  libmesh_error();
1052  }
1053  }
1054  }
1055 
1056  else // edge_swap==true
1057  {
1058  switch (sn)
1059  {
1060  case 0:
1061  {
1062  // New remote side is Tri 0, side 0
1063  tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
1064  break;
1065  }
1066  case 1:
1067  {
1068  // New remote side is Tri 1, side 0
1069  tri1->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
1070  break;
1071  }
1072  case 2:
1073  {
1074  // New remote side is Tri 1, side 1
1075  tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
1076  break;
1077  }
1078  case 3:
1079  {
1080  // New remote side is Tri 0, side 2
1081  tri0->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
1082  break;
1083  }
1084 
1085  default:
1086  {
1087  libMesh::err << "Quad4/8/9 cannot have more than 4 sides." << std::endl;
1088  libmesh_error();
1089  }
1090  }
1091  } // end edge_swap==true
1092  } // end if (elem->neighbor(sn) == remote_elem)
1093  } // end for loop over sides
1094 
1095  // Determine new IDs for the split elements which will be
1096  // the same on all processors, therefore keeping the Mesh
1097  // in sync. Note: we offset the new IDs by n_orig_elem to
1098  // avoid overwriting any of the original IDs, this assumes
1099  // they were contiguously-numbered to begin with...
1100  tri0->set_id( n_orig_elem + 2*elem->id() + 0 );
1101  tri1->set_id( n_orig_elem + 2*elem->id() + 1 );
1102 
1103  // Add the newly-created triangles to the temporary vector of new elements.
1104  new_elements.push_back(tri0);
1105  new_elements.push_back(tri1);
1106 
1107  // Delete the original element
1108  mesh.delete_elem(elem);
1109  } // end if (split_elem)
1110  } // End for loop over elements
1111  } // end scope
1112 
1113 
1114  // Now, iterate over the new elements vector, and add them each to
1115  // the Mesh.
1116  {
1117  std::vector<Elem*>::iterator el = new_elements.begin();
1118  const std::vector<Elem*>::iterator end = new_elements.end();
1119  for (; el != end; ++el)
1120  mesh.add_elem(*el);
1121  }
1122 
1123  if (mesh_has_boundary_data)
1124  {
1125  // By this time, we should have removed all of the original boundary sides
1126  // - except on a hybrid mesh, where we can't "start from a blank slate"! - RHS
1127  // libmesh_assert_equal_to (mesh.boundary_info->n_boundary_conds(), 0);
1128 
1129  // Clear the boundary info, to be sure and start from a blank slate.
1130  // mesh.boundary_info->clear();
1131 
1132  // If the old mesh had boundary data, the new mesh better have some.
1133  libmesh_assert_greater (new_bndry_elements.size(), 0);
1134 
1135  // We should also be sure that the lengths of the new boundary data vectors
1136  // are all the same.
1137  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1138  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1139 
1140  // Add the new boundary info to the mesh
1141  for (unsigned int s=0; s<new_bndry_elements.size(); ++s)
1142  mesh.boundary_info->add_side(new_bndry_elements[s],
1143  new_bndry_sides[s],
1144  new_bndry_ids[s]);
1145  }
1146 
1147 
1148  // Prepare the newly created mesh for use.
1149  mesh.prepare_for_use(/*skip_renumber =*/ false);
1150 
1151  // Let the new_elements and new_bndry_elements vectors go out of scope.
1152 }
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 1463 of file mesh_modification.C.

References libMesh::MeshBase::boundary_info, libMesh::Elem::get_node(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), libMesh::Elem::n_edges(), n_nodes, libMesh::Elem::n_nodes(), and libMesh::Elem::n_sides().

1466 {
1467  // Only level-0 elements store BCs. Loop over them.
1468  MeshBase::element_iterator el = mesh.level_elements_begin(0);
1469  const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
1470  for (; el != end_el; ++el)
1471  {
1472  Elem *elem = *el;
1473 
1474  unsigned int n_nodes = elem->n_nodes();
1475  for (unsigned int n=0; n != n_nodes; ++n)
1476  {
1477  const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->boundary_ids(elem->get_node(n));
1478  if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
1479  {
1480  std::vector<boundary_id_type> new_ids(old_ids);
1481  std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
1482  mesh.boundary_info->remove(elem->get_node(n));
1483  mesh.boundary_info->add_node(elem->get_node(n), new_ids);
1484  }
1485  }
1486 
1487  unsigned int n_edges = elem->n_edges();
1488  for (unsigned int edge=0; edge != n_edges; ++edge)
1489  {
1490  const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->edge_boundary_ids(elem, edge);
1491  if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
1492  {
1493  std::vector<boundary_id_type> new_ids(old_ids);
1494  std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
1495  mesh.boundary_info->remove_edge(elem, edge);
1496  mesh.boundary_info->add_edge(elem, edge, new_ids);
1497  }
1498  }
1499 
1500  unsigned int n_sides = elem->n_sides();
1501  for (unsigned int s=0; s != n_sides; ++s)
1502  {
1503  const std::vector<boundary_id_type>& old_ids = mesh.boundary_info->boundary_ids(elem, s);
1504  if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
1505  {
1506  std::vector<boundary_id_type> new_ids(old_ids);
1507  std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
1508  mesh.boundary_info->remove_side(elem, s);
1509  mesh.boundary_info->add_side(elem, s, new_ids);
1510  }
1511  }
1512  }
1513 }
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 1517 of file mesh_modification.C.

References libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), and libMesh::Elem::subdomain_id().

1520 {
1521  MeshBase::element_iterator el = mesh.elements_begin();
1522  const MeshBase::element_iterator end_el = mesh.elements_end();
1523 
1524  for (; el != end_el; ++el)
1525  {
1526  Elem *elem = *el;
1527 
1528  if (elem->subdomain_id() == old_id)
1529  elem->subdomain_id() = new_id;
1530  }
1531 }
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 47 of file mesh_modification.C.

References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::libmesh_assert(), std::max(), libMesh::MeshBase::mesh_dimension(), std::min(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::START_LOG(), libMesh::STOP_LOG(), and libMesh::TypeVector< T >::unit().

50 {
51  libmesh_assert (mesh.n_nodes());
52  libmesh_assert (mesh.n_elem());
53  libmesh_assert ((factor >= 0.) && (factor <= 1.));
54 
55  START_LOG("distort()", "MeshTools::Modification");
56 
57 
58 
59  // First find nodes on the boundary and flag them
60  // so that we don't move them
61  // on_boundary holds false (not on boundary) and true (on boundary)
62  std::vector<bool> on_boundary (mesh.n_nodes(), false);
63 
64  if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);
65 
66  // Now calculate the minimum distance to
67  // neighboring nodes for each node.
68  // hmin holds these distances.
69  std::vector<float> hmin (mesh.n_nodes(),
71 
72  MeshBase::element_iterator el = mesh.active_elements_begin();
73  const MeshBase::element_iterator end = mesh.active_elements_end();
74 
75  for (; el!=end; ++el)
76  for (unsigned int n=0; n<(*el)->n_nodes(); n++)
77  hmin[(*el)->node(n)] = std::min(hmin[(*el)->node(n)],
78  static_cast<float>((*el)->hmin()));
79 
80 
81  // Now actually move the nodes
82  {
83  const unsigned int seed = 123456;
84 
85  // seed the random number generator.
86  // We'll loop from 1 to n_nodes on every processor, even those
87  // that don't have a particular node, so that the pseudorandom
88  // numbers will be the same everywhere.
89  std::srand(seed);
90 
91  // If the node is on the boundary or
92  // the node is not used by any element (hmin[n]<1.e20)
93  // then we should not move it.
94  // [Note: Testing for (in)equality might be wrong
95  // (different types, namely float and double)]
96  for (unsigned int n=0; n<mesh.n_nodes(); n++)
97  if (!on_boundary[n] && (hmin[n] < 1.e20) )
98  {
99  // the direction, random but unit normalized
100 
101  Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
102  (mesh.mesh_dimension() > 1) ?
103  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
104  : 0.,
105  ((mesh.mesh_dimension() == 3) ?
106  static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
107  : 0.)
108  );
109 
110  dir(0) = (dir(0)-.5)*2.;
111  if (mesh.mesh_dimension() > 1)
112  dir(1) = (dir(1)-.5)*2.;
113  if (mesh.mesh_dimension() == 3)
114  dir(2) = (dir(2)-.5)*2.;
115 
116  dir = dir.unit();
117 
118  Node *node = mesh.node_ptr(n);
119  if (!node)
120  continue;
121 
122  (*node)(0) += dir(0)*factor*hmin[n];
123  if (mesh.mesh_dimension() > 1)
124  (*node)(1) += dir(1)*factor*hmin[n];
125  if (mesh.mesh_dimension() == 3)
126  (*node)(2) += dir(2)*factor*hmin[n];
127  }
128  }
129 
130 
131  // All done
132  STOP_LOG("distort()", "MeshTools::Modification");
133 }
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 that 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 1337 of file mesh_modification.C.

References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::MeshBase::add_elem(), bc_id, libMesh::MeshBase::boundary_info, libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::n_active_elem(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), 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::Elem::subdomain_id(), and libMesh::Elem::type().

1338 {
1339  // Algorithm:
1340  // .) For each active element in the mesh: construct a
1341  // copy which is the same in every way *except* it is
1342  // a level 0 element. Store the pointers to these in
1343  // a separate vector. Save any boundary information as well.
1344  // Delete the active element from the mesh.
1345  // .) Loop over all (remaining) elements in the mesh, delete them.
1346  // .) Add the level-0 copies back to the mesh
1347 
1348  // Temporary storage for new element pointers
1349  std::vector<Elem*> new_elements;
1350 
1351  // BoundaryInfo Storage for element ids, sides, and BC ids
1352  std::vector<Elem*> saved_boundary_elements;
1353  std::vector<boundary_id_type> saved_bc_ids;
1354  std::vector<unsigned short int> saved_bc_sides;
1355 
1356  // Reserve a reasonable amt. of space for each
1357  new_elements.reserve(mesh.n_active_elem());
1358  saved_boundary_elements.reserve(mesh.boundary_info->n_boundary_conds());
1359  saved_bc_ids.reserve(mesh.boundary_info->n_boundary_conds());
1360  saved_bc_sides.reserve(mesh.boundary_info->n_boundary_conds());
1361  {
1362  MeshBase::element_iterator it = mesh.active_elements_begin();
1363  const MeshBase::element_iterator end = mesh.active_elements_end();
1364 
1365  for (; it != end; ++it)
1366  {
1367  Elem* elem = *it;
1368 
1369  // Make a new element of the same type
1370  Elem* copy = Elem::build(elem->type()).release();
1371 
1372  // Set node pointers (they still point to nodes in the original mesh)
1373  for(unsigned int n=0; n<elem->n_nodes(); n++)
1374  copy->set_node(n) = elem->get_node(n);
1375 
1376  // Copy over ids
1377  copy->processor_id() = elem->processor_id();
1378  copy->subdomain_id() = elem->subdomain_id();
1379 
1380  // Retain the original element's ID as well, otherwise ParallelMesh will
1381  // try to create one for you...
1382  copy->set_id( elem->id() );
1383 
1384  // This element could have boundary info or ParallelMesh
1385  // remote_elem links as well. We need to save the (elem,
1386  // side, bc_id) triples and those links
1387  for (unsigned int s=0; s<elem->n_sides(); s++)
1388  {
1389  if (elem->neighbor(s) == remote_elem)
1390  copy->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
1391 
1392  const std::vector<boundary_id_type>& bc_ids = mesh.boundary_info->boundary_ids(elem,s);
1393  for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
1394  {
1395  const boundary_id_type bc_id = *id_it;
1396 
1397  if (bc_id != BoundaryInfo::invalid_id)
1398  {
1399  saved_boundary_elements.push_back(copy);
1400  saved_bc_ids.push_back(bc_id);
1401  saved_bc_sides.push_back(s);
1402  }
1403  }
1404  }
1405 
1406 
1407  // We're done with this element
1408  mesh.delete_elem(elem);
1409 
1410  // But save the copy
1411  new_elements.push_back(copy);
1412  }
1413 
1414  // Make sure we saved the same number of boundary conditions
1415  // in each vector.
1416  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1417  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1418  }
1419 
1420 
1421  // Loop again, delete any remaining elements
1422  {
1423  MeshBase::element_iterator it = mesh.elements_begin();
1424  const MeshBase::element_iterator end = mesh.elements_end();
1425 
1426  for (; it != end; ++it)
1427  mesh.delete_elem( *it );
1428  }
1429 
1430 
1431  // Add the copied (now level-0) elements back to the mesh
1432  {
1433  for (std::vector<Elem*>::iterator it = new_elements.begin();
1434  it != new_elements.end();
1435  ++it)
1436  {
1437  dof_id_type orig_id = (*it)->id();
1438 
1439  Elem* added_elem = mesh.add_elem(*it);
1440 
1441  dof_id_type added_id = added_elem->id();
1442 
1443  // If the Elem, as it was re-added to the mesh, now has a
1444  // different ID (this is unlikely, so it's just an assert)
1445  // the boundary information will no longer be correct.
1446  libmesh_assert_equal_to (orig_id, added_id);
1447  }
1448  }
1449 
1450  // Finally, also add back the saved boundary information
1451  for (unsigned int e=0; e<saved_boundary_elements.size(); ++e)
1452  mesh.boundary_info->add_side(saved_boundary_elements[e],
1453  saved_bc_sides[e],
1454  saved_bc_ids[e]);
1455 
1456  // Trim unused and renumber nodes and elements
1457  mesh.prepare_for_use(/*skip_renumber =*/ false);
1458 }
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 174 of file mesh_modification.C.

References libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::pi, and libMesh::Real.

178 {
179  libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
180 
181  const Real p = -phi/180.*libMesh::pi;
182  const Real t = -theta/180.*libMesh::pi;
183  const Real s = -psi/180.*libMesh::pi;
184  const Real sp = std::sin(p), cp = std::cos(p);
185  const Real st = std::sin(t), ct = std::cos(t);
186  const Real ss = std::sin(s), cs = std::cos(s);
187 
188  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
189  // (equations 6-14 give the entries of the composite transformation matrix).
190  // The rotations are performed sequentially about the z, x, and z axes, in that order.
191  // A positive angle yields a counter-clockwise rotation about the axis in question.
192  const MeshBase::node_iterator nd_end = mesh.nodes_end();
193 
194  for (MeshBase::node_iterator nd = mesh.nodes_begin();
195  nd != nd_end; ++nd)
196  {
197  const Point pt = **nd;
198  const Real x = pt(0);
199  const Real y = pt(1);
200  const Real z = pt(2);
201  **nd = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
202  (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
203  ( sp*st)*x + (-cp*st)*y + (ct)*z );
204  }
205 }
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 208 of file mesh_modification.C.

References libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::Real, and libMesh::MeshBase::spatial_dimension().

Referenced by libMesh::DenseVector< T >::operator*=(), and libMesh::DenseMatrix< T >::operator*=().

212 {
213  const Real x_scale = xs;
214  Real y_scale = ys;
215  Real z_scale = zs;
216 
217  if (ys == 0.)
218  {
219  libmesh_assert_equal_to (zs, 0.);
220 
221  y_scale = z_scale = x_scale;
222  }
223 
224  // Scale the x coordinate in all dimensions
225  const MeshBase::node_iterator nd_end = mesh.nodes_end();
226 
227  for (MeshBase::node_iterator nd = mesh.nodes_begin();
228  nd != nd_end; ++nd)
229  (**nd)(0) *= x_scale;
230 
231 
232  // Only scale the y coordinate in 2 and 3D
233  if (mesh.spatial_dimension() < 2)
234  return;
235 
236  for (MeshBase::node_iterator nd = mesh.nodes_begin();
237  nd != nd_end; ++nd)
238  (**nd)(1) *= y_scale;
239 
240  // Only scale the z coordinate in 3D
241  if (mesh.spatial_dimension() < 3)
242  return;
243 
244  for (MeshBase::node_iterator nd = mesh.nodes_begin();
245  nd != nd_end; ++nd)
246  (**nd)(2) *= z_scale;
247 }
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 postition 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 Lüthi (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 1155 of file mesh_modification.C.

References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::build_side(), libMesh::Elem::child(), libMesh::Elem::embedding_matrix(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_children(), libMesh::MeshTools::n_levels(), libMesh::Elem::n_neighbors(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::neighbor(), libMesh::MeshBase::node(), libMesh::Elem::parent(), libMesh::Elem::point(), std::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), side, libMesh::TypeVector< T >::size(), and libMesh::MeshTools::weight().

1158 {
1162  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1163 
1164  /*
1165  * find the boundary nodes
1166  */
1167  std::vector<bool> on_boundary;
1168  MeshTools::find_boundary_nodes(mesh, on_boundary);
1169 
1170  for (unsigned int iter=0; iter<n_iterations; iter++)
1171 
1172  {
1173  /*
1174  * loop over the mesh refinement level
1175  */
1176  unsigned int n_levels = MeshTools::n_levels(mesh);
1177  for (unsigned int refinement_level=0; refinement_level != n_levels;
1178  refinement_level++)
1179  {
1180  // initialize the storage (have to do it on every level to get empty vectors
1181  std::vector<Point> new_positions;
1182  std::vector<Real> weight;
1183  new_positions.resize(mesh.n_nodes());
1184  weight.resize(mesh.n_nodes());
1185 
1186  {
1187  /*
1188  * Loop over the elements to calculate new node positions
1189  */
1190  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1191  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1192 
1193  for (; el != end; ++el)
1194  {
1195  /*
1196  * Constant handle for the element
1197  */
1198  const Elem* elem = *el;
1199 
1200  /*
1201  * We relax all nodes on level 0 first
1202  * If the element is refined (level > 0), we interpolate the
1203  * parents nodes with help of the embedding matrix
1204  */
1205  if (refinement_level == 0)
1206  {
1207  for (unsigned int s=0; s<elem->n_neighbors(); s++)
1208  {
1209  /*
1210  * Only operate on sides which are on the
1211  * boundary or for which the current element's
1212  * id is greater than its neighbor's.
1213  * Sides get only built once.
1214  */
1215  if ((elem->neighbor(s) != NULL) &&
1216  (elem->id() > elem->neighbor(s)->id()) )
1217  {
1218  AutoPtr<Elem> side(elem->build_side(s));
1219 
1220  Node* node0 = side->get_node(0);
1221  Node* node1 = side->get_node(1);
1222 
1223  Real node_weight = 1.;
1224  // calculate the weight of the nodes
1225  if (power > 0)
1226  {
1227  Point diff = (*node0)-(*node1);
1228  node_weight = std::pow( diff.size(), power );
1229  }
1230 
1231  const dof_id_type id0 = node0->id(), id1 = node1->id();
1232  new_positions[id0].add_scaled( *node1, node_weight );
1233  new_positions[id1].add_scaled( *node0, node_weight );
1234  weight[id0] += node_weight;
1235  weight[id1] += node_weight;
1236  }
1237  } // element neighbor loop
1238  }
1239 #ifdef LIBMESH_ENABLE_AMR
1240  else // refinement_level > 0
1241  {
1242  /*
1243  * Find the positions of the hanging nodes of refined elements.
1244  * We do this by calculating their position based on the parent
1245  * (one level less refined) element, and the embedding matrix
1246  */
1247 
1248  const Elem* parent = elem->parent();
1249 
1250  /*
1251  * find out which child I am
1252  */
1253  for (unsigned int c=0; c < parent->n_children(); c++)
1254  {
1255  if (parent->child(c) == elem)
1256  {
1257  /*
1258  *loop over the childs (that is, the current elements) nodes
1259  */
1260  for (unsigned int nc=0; nc < elem->n_nodes(); nc++)
1261  {
1262  /*
1263  * the new position of the node
1264  */
1265  Point point;
1266  for (unsigned int n=0; n<parent->n_nodes(); n++)
1267  {
1268  /*
1269  * The value from the embedding matrix
1270  */
1271  const float em_val = parent->embedding_matrix(c,nc,n);
1272 
1273  if (em_val != 0.)
1274  point.add_scaled (parent->point(n), em_val);
1275  }
1276 
1277  const dof_id_type id = elem->get_node(nc)->id();
1278  new_positions[id] = point;
1279  weight[id] = 1.;
1280  }
1281 
1282  } // if parent->child == elem
1283  } // for parent->n_children
1284  } // if element refinement_level
1285 #endif // #ifdef LIBMESH_ENABLE_AMR
1286 
1287  } // element loop
1288 
1289  /*
1290  * finally reposition the vertex nodes
1291  */
1292  for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
1293  if (!on_boundary[nid] && weight[nid] > 0.)
1294  mesh.node(nid) = new_positions[nid]/weight[nid];
1295  }
1296 
1297  {
1298  /*
1299  * Now handle the additional second_order nodes by calculating
1300  * their position based on the vertex postitions
1301  * we do a second loop over the level elements
1302  */
1303  MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
1304  const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
1305 
1306  for (; el != end; ++el)
1307  {
1308  /*
1309  * Constant handle for the element
1310  */
1311  const Elem* elem = *el;
1312  const unsigned int son_begin = elem->n_vertices();
1313  const unsigned int son_end = elem->n_nodes();
1314  for (unsigned int n=son_begin; n<son_end; n++)
1315  {
1316  const unsigned int n_adjacent_vertices =
1317  elem->n_second_order_adjacent_vertices(n);
1318 
1319  Point point;
1320  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1321  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1322 
1323  const dof_id_type id = elem->get_node(n)->id();
1324  mesh.node(id) = point/n_adjacent_vertices;
1325  }
1326  }
1327  }
1328 
1329  } // refinement_level loop
1330 
1331  } // end iteration
1332 }
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 137 of file mesh_modification.C.

References libMesh::MeshBase::nodes_begin(), and libMesh::MeshBase::nodes_end().

141 {
142  const Point p(xt, yt, zt);
143 
144  const MeshBase::node_iterator nd_end = mesh.nodes_end();
145 
146  for (MeshBase::node_iterator nd = mesh.nodes_begin();
147  nd != nd_end; ++nd)
148  **nd += p;
149 }

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:33 UTC

Hosted By:
SourceForge.net Logo