26 #ifndef OOMPH_TRIANGLE_MESH_HEADER
27 #define OOMPH_TRIANGLE_MESH_HEADER
30 #include <oomph-lib-config.h>
39 #ifdef OOMPH_HAS_FPUCONTROLH
40 #include <fpu_control.h>
51 #include "../generic/problem.h"
52 #include "../generic/triangle_scaffold_mesh.h"
53 #include "../generic/triangle_mesh.h"
54 #include "../generic/refineable_mesh.h"
55 #include "../rigid_body/immersed_rigid_body_elements.h"
59 #ifdef OOMPH_HAS_TRIANGLE_LIB
74 struct oomph::TriangulateIO* in,
75 struct oomph::TriangulateIO* out,
76 struct oomph::TriangulateIO* vorout);
214 Vector<double>& region_coordinates)
219 std::ostringstream error_message;
221 <<
"Please use another region id different from zero.\n"
222 <<
"It is internally used as the default region number.\n";
223 throw OomphLibError(error_message.str(),
224 OOMPH_CURRENT_FUNCTION,
225 OOMPH_EXCEPTION_LOCATION);
229 std::map<unsigned, Vector<double>>::iterator it;
235 std::ostringstream error_message;
236 error_message <<
"The region id (" << i <<
") that you are using for"
238 <<
"your region is already in use. Use another\n"
239 <<
"region id and verify that you are not re-using\n"
240 <<
" previously defined regions ids\n"
242 OomphLibWarning(error_message.str(),
243 OOMPH_CURRENT_FUNCTION,
244 OOMPH_EXCEPTION_LOCATION);
422 template<
class ELEMENT>
429 #ifdef OOMPH_HAS_TRIANGLE_LIB
434 this->Allow_automatic_creation_of_vertices_on_boundaries =
true;
444 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
449 const std::string& node_file_name,
450 const std::string& element_file_name,
451 const std::string& poly_file_name,
452 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
453 const bool& allow_automatic_creation_of_vertices_on_boundaries =
true)
456 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
459 this->Allow_automatic_creation_of_vertices_on_boundaries =
460 allow_automatic_creation_of_vertices_on_boundaries;
473 bool should_use_attributes =
false;
475 #ifdef OOMPH_HAS_TRIANGLE_LIB
477 TriangleHelper::create_triangulateio_from_polyfiles(
482 should_use_attributes);
493 node_file_name, element_file_name, poly_file_name);
503 unsigned nb = nboundary();
504 for (
unsigned b = 0; b < nb; b++)
506 this->
template setup_boundary_coordinates<ELEMENT>(b);
511 #ifdef OOMPH_HAS_TRIANGLE_LIB
517 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
523 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
526 this->Allow_automatic_creation_of_vertices_on_boundaries =
527 triangle_mesh_parameters
547 unsigned max_boundary_id = 0;
554 Vector<TriangleMeshClosedCurve*> outer_boundary_pt =
559 if (outer_boundary_pt.size() == 0)
561 std::stringstream error_message;
563 <<
"There are no outer boundaries defined.\n"
564 <<
"Verify that you have specified the outer boundaries in the\n"
565 <<
"Triangle_mesh_parameter object\n\n";
566 throw OomphLibError(error_message.str(),
567 OOMPH_CURRENT_FUNCTION,
568 OOMPH_EXCEPTION_LOCATION);
573 unsigned n_outer_boundaries = outer_boundary_pt.size();
577 Vector<TriangleMeshPolygon*> outer_boundary_polygon_pt(
581 for (
unsigned i = 0; i < n_outer_boundaries; ++i)
587 outer_boundary_polygon_pt[i] =
588 closed_curve_to_polygon_helper(outer_boundary_pt[i], max_boundary_id);
596 Vector<TriangleMeshClosedCurve*> internal_closed_curve_pt =
600 unsigned n_internal_closed_curves = internal_closed_curve_pt.size();
605 Vector<TriangleMeshPolygon*> internal_polygon_pt(
606 n_internal_closed_curves);
609 for (
unsigned i = 0; i < n_internal_closed_curves; ++i)
613 internal_polygon_pt[i] = closed_curve_to_polygon_helper(
614 internal_closed_curve_pt[i], max_boundary_id);
622 Vector<TriangleMeshOpenCurve*> internal_open_curve_pt =
626 unsigned n_internal_open_curves = internal_open_curve_pt.size();
629 Vector<TriangleMeshOpenCurve*> internal_open_curve_poly_pt(
630 n_internal_open_curves);
633 for (
unsigned i = 0; i < n_internal_open_curves; i++)
638 internal_open_curve_poly_pt[i] = create_open_curve_with_polyline_helper(
639 internal_open_curve_pt[i], max_boundary_id);
649 for (
unsigned i = 0; i < n_outer_boundaries; i++)
651 set_geom_objects_and_coordinate_limits_for_close_curve(
652 outer_boundary_pt[i]);
658 for (
unsigned i = 0; i < n_internal_closed_curves; i++)
660 set_geom_objects_and_coordinate_limits_for_close_curve(
661 internal_closed_curve_pt[i]);
667 for (
unsigned i = 0; i < n_internal_open_curves; i++)
669 set_geom_objects_and_coordinate_limits_for_open_curve(
670 internal_open_curve_pt[i]);
679 const double element_area = triangle_mesh_parameters.
element_area();
682 Vector<Vector<double>> extra_holes_coordinates =
686 std::map<unsigned, Vector<double>> regions =
692 const bool refine_boundary =
695 const bool refine_internal_boundary =
698 if (!refine_internal_boundary && refine_boundary)
700 std::ostringstream error_stream;
702 <<
"You have specified that Triangle may refine the outer boundary, "
704 <<
"not internal boundaries. Triangle does not support this "
706 <<
"If you do not want Triangle to refine internal boundaries, it "
708 <<
"refine outer boundaries either!\n"
709 <<
"Please either disable all boundary refinement\n"
710 <<
"(call TriangleMeshParameters::disable_boundary_refinement()\n"
711 <<
"or enable internal boundary refinement (the default)\n";
713 throw OomphLibError(error_stream.str().c_str(),
714 OOMPH_CURRENT_FUNCTION,
715 OOMPH_EXCEPTION_LOCATION);
719 outer_boundary_polygon_pt,
721 internal_open_curve_poly_pt,
723 extra_holes_coordinates,
729 refine_internal_boundary);
732 unsigned nb = nboundary();
744 for (
unsigned b = 0; b < nb; b++)
746 this->
template setup_boundary_coordinates<ELEMENT>(b);
750 this->snap_nodes_onto_geometric_objects();
756 const std::string& poly_file_name,
757 const double& element_area,
758 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
759 const bool& allow_automatic_creation_of_vertices_on_boundaries =
true)
762 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
765 this->Allow_automatic_creation_of_vertices_on_boundaries =
766 allow_automatic_creation_of_vertices_on_boundaries;
775 std::string message =
776 "This constructor hasn't been tested since last cleanup.\n";
778 message,
"TriangleMesh::TriangleMesh()", OOMPH_EXCEPTION_LOCATION);
784 TriangulateIO triangle_in;
787 std::stringstream input_string_stream;
791 input_string_stream <<
"-pA -a -a" << element_area <<
"q30";
794 if (!this->is_creation_of_vertices_on_boundaries_allowed())
796 input_string_stream <<
" -YY";
800 char triswitches[100];
801 sprintf(triswitches,
"%s", input_string_stream.str().c_str());
806 bool use_attributes =
false;
815 triangulate(triswitches, &triangle_in, &Triangulateio, 0);
817 #ifdef OOMPH_HAS_FPUCONTROLH
819 fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
824 this->
Tmp_mesh_pt =
new TriangleScaffoldMesh(Triangulateio);
834 bool clear_hole_data =
false;
835 TriangleHelper::clear_triangulateio(triangle_in, clear_hole_data);
838 unsigned nb = nboundary();
839 for (
unsigned b = 0; b < nb; b++)
841 this->
template setup_boundary_coordinates<ELEMENT>(b);
856 #ifdef OOMPH_HAS_TRIANGLE_LIB
859 TriangleHelper::clear_triangulateio(Triangulateio);
862 std::set<TriangleMeshCurveSection*>::iterator it_polyline;
863 for (it_polyline = Free_curve_section_pt.begin();
864 it_polyline != Free_curve_section_pt.end();
867 delete (*it_polyline);
870 std::set<TriangleMeshPolygon*>::iterator it_polygon;
871 for (it_polygon = Free_polygon_pt.begin();
872 it_polygon != Free_polygon_pt.end();
875 delete (*it_polygon);
878 std::set<TriangleMeshOpenCurve*>::iterator it_open_polyline;
879 for (it_open_polyline = Free_open_curve_pt.begin();
880 it_open_polyline != Free_open_curve_pt.end();
883 delete (*it_open_polyline);
892 const bool& preserve_existing_data)
910 Vector<std::list<FiniteElement*>>& old_segment_sorted_ele_pt,
911 std::map<FiniteElement*, bool>& old_is_inverted);
925 Vector<FiniteElement*>& input_face_ele_pt,
926 const bool& is_internal_boundary,
927 std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt);
945 Vector<FiniteElement*>& face_el_pt,
947 bool& is_internal_boundary,
948 std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt);
954 return Boundary_segment_node_pt[b];
962 return Boundary_segment_node_pt[b][s];
970 return Boundary_segment_node_pt[b][s][n];
975 #ifdef OOMPH_HAS_TRIANGLE_LIB
983 unsigned nhole = Triangulateio.numberofholes;
984 unsigned count_coord = 0;
985 for (
unsigned ihole = 0; ihole < nhole; ihole++)
987 Triangulateio.holelist[count_coord] += internal_point[ihole][0];
988 Triangulateio.holelist[count_coord + 1] += internal_point[ihole][1];
1002 unsigned nnode = Triangulateio.numberofpoints;
1007 for (
unsigned inod = 0; inod < nnode; inod++)
1014 Node* mesh_node_pt = this->node_pt(inod);
1015 new_x = mesh_node_pt->x(0);
1016 new_y = mesh_node_pt->x(1);
1017 Triangulateio.pointlist[count * 2] = new_x;
1018 Triangulateio.pointlist[(count * 2) + 1] = new_y;
1022 #ifdef OOMPH_HAS_MPI
1028 std::string input_string;
1031 getline(read_file, input_string,
'#');
1034 read_file.ignore(200,
'\n');
1037 return std::atoi(input_string.c_str());
1046 OomphCommunicator* comm_pt, std::istream& restart_file)
1048 std::ostringstream error_stream;
1049 error_stream <<
"Empty default reestablish disributed info method "
1051 error_stream <<
"This should be overloaded in a specific "
1052 <<
"RefineableTriangleMesh\n";
1053 throw OomphLibError(
1054 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1063 this->remove_boundary_nodes();
1066 unsigned n_node = this->nnode();
1067 for (
unsigned n = n_node; n > 0; --n)
1069 delete this->Node_pt[n - 1];
1070 this->Node_pt[n - 1] = 0;
1073 unsigned n_element = this->nelement();
1074 for (
unsigned e = n_element; e > 0; --e)
1076 delete this->Element_pt[e - 1];
1077 this->Element_pt[e - 1] = 0;
1080 this->flush_element_and_node_storage();
1084 this->Boundary_element_pt.clear();
1085 this->Face_index_at_boundary.clear();
1086 this->Region_element_pt.clear();
1087 this->Region_attribute.clear();
1088 this->Boundary_region_element_pt.clear();
1089 this->Face_index_region_at_boundary.clear();
1090 this->Boundary_curve_section_pt.clear();
1091 this->Polygonal_vertex_arclength_info.clear();
1093 #ifdef OOMPH_HAS_MPI
1095 if (this->is_mesh_distributed())
1097 this->Halo_node_pt.clear();
1098 this->Root_halo_element_pt.clear();
1100 this->Haloed_node_pt.clear();
1101 this->Root_haloed_element_pt.clear();
1103 this->External_halo_node_pt.clear();
1104 this->External_halo_element_pt.clear();
1106 this->External_haloed_node_pt.clear();
1107 this->External_haloed_element_pt.clear();
1111 unsigned nbound = nboundary();
1112 Boundary_coordinate_exists.resize(nbound,
false);
1115 this->
Tmp_mesh_pt =
new TriangleScaffoldMesh(this->Triangulateio);
1127 #ifdef OOMPH_HAS_MPI
1128 if (!this->is_mesh_distributed())
1130 nbound = this->nboundary();
1143 nbound = this->nboundary();
1147 for (
unsigned b = 0; b < nbound; b++)
1149 this->
template setup_boundary_coordinates<ELEMENT>(b);
1156 if (!this->is_mesh_distributed())
1159 this->snap_nodes_onto_geometric_objects();
1192 const bool& use_attributes);
1194 #ifdef OOMPH_HAS_TRIANGLE_LIB
1199 TriangulateIO& triangulate_io,
1200 bool& use_attributes);
1206 Vector<TriangleMeshPolygon*>& outer_boundary_pt,
1207 Vector<TriangleMeshPolygon*>& internal_polygon_pt,
1208 Vector<TriangleMeshOpenCurve*>& open_polylines_pt,
1209 const double& element_area,
1210 Vector<Vector<double>>& extra_holes_coordinates,
1211 std::map<
unsigned, Vector<double>>& regions_coordinates,
1212 std::map<unsigned, double>& regions_areas,
1213 TimeStepper* time_stepper_pt,
1214 const bool& use_attributes,
1215 const bool& refine_boundary,
1216 const bool& refine_internal_boundary)
1219 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
1222 if (element_area < 10e-14)
1224 std::ostringstream warning_message;
1226 <<
"The current elements area was stated to (" << element_area
1227 <<
").\nThe current precision to generate the input to triangle "
1228 <<
"is fixed to 14 digits\n\n";
1229 OomphLibWarning(warning_message.str(),
1230 OOMPH_CURRENT_FUNCTION,
1231 OOMPH_EXCEPTION_LOCATION);
1242 Outer_boundary_pt = outer_boundary_pt;
1245 Internal_polygon_pt = internal_polygon_pt;
1248 Internal_open_curve_pt = open_polylines_pt;
1251 Extra_holes_coordinates = extra_holes_coordinates;
1254 Regions_coordinates = regions_coordinates;
1257 TriangulateIO triangulate_io;
1260 TriangleHelper::initialise_triangulateio(triangulate_io);
1264 UnstructuredTwoDMeshGeometryBase::build_triangulateio(
1266 internal_polygon_pt,
1268 extra_holes_coordinates,
1269 regions_coordinates,
1274 TriangleHelper::initialise_triangulateio(Triangulateio);
1280 std::stringstream input_string_stream;
1281 input_string_stream.precision(14);
1282 input_string_stream.setf(std::ios_base::fixed, std::ios_base::floatfield);
1288 input_string_stream <<
"-pA -a -a" << element_area <<
" -q30"
1292 if (!this->is_automatic_creation_of_vertices_on_boundaries_allowed())
1294 input_string_stream <<
" -YY";
1298 if (refine_boundary ==
false)
1300 input_string_stream <<
"-Y";
1302 if (refine_internal_boundary ==
false)
1304 input_string_stream <<
"Y";
1309 char triswitches[100];
1310 sprintf(triswitches,
"%s", input_string_stream.str().c_str());
1313 triangulate(triswitches, &triangulate_io, &Triangulateio, 0);
1315 #ifdef OOMPH_HAS_FPUCONTROLH
1317 fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
1322 this->
Tmp_mesh_pt =
new TriangleScaffoldMesh(Triangulateio);
1325 if (!regions_coordinates.empty())
1344 bool clear_hole_data =
false;
1345 TriangleHelper::clear_triangulateio(triangulate_io, clear_hole_data);
1360 #ifdef OOMPH_HAS_MPI
1383 std::set<unsigned> shared_boundaries_in_this_processor_set;
1386 const unsigned n_proc = this->communicator_pt()->nproc();
1388 const unsigned my_rank = this->communicator_pt()->my_rank();
1391 for (
unsigned iproc = 0; iproc < n_proc; iproc++)
1394 if (iproc != my_rank)
1398 const unsigned nshared_boundaries_with_iproc =
1403 if (nshared_boundaries_with_iproc > 0)
1406 Vector<unsigned> bound_ids_shared_with_iproc;
1407 bound_ids_shared_with_iproc =
1411 for (
unsigned bs = 0; bs < nshared_boundaries_with_iproc; bs++)
1413 const unsigned bnd_id = bound_ids_shared_with_iproc[bs];
1417 std::set<unsigned>::iterator it =
1418 shared_boundaries_in_this_processor_set.find(bnd_id);
1419 if (it != shared_boundaries_in_this_processor_set.end())
1421 std::stringstream error;
1422 error <<
"The current shared boundary (" << bnd_id <<
") was\n"
1423 <<
"already added by other pair of processors\n."
1424 <<
"This means that there are repeated shared boundaries "
1426 throw OomphLibError(error.str(),
1427 OOMPH_CURRENT_FUNCTION,
1428 OOMPH_EXCEPTION_LOCATION);
1430 shared_boundaries_in_this_processor_set.insert(bnd_id);
1444 const unsigned& q)
const
1470 const unsigned& q)
const
1483 const unsigned& i)
const
1494 const unsigned& c)
const
1500 const unsigned& p,
const unsigned& c)
1507 const unsigned& i)
const
1520 std::map<unsigned, Vector<FiniteElement*>>::iterator it =
1528 std::ostringstream error_stream;
1529 error_stream <<
"The shared boundary (" << b
1530 <<
") does not exist!!!\n\n";
1531 throw OomphLibError(
1532 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1544 std::map<unsigned, Vector<FiniteElement*>>::iterator it =
1552 std::ostringstream error_stream;
1553 error_stream <<
"The shared boundary (" << b
1554 <<
") does not exist!!!\n\n";
1555 throw OomphLibError(
1556 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1569 std::map<unsigned, Vector<FiniteElement*>>::iterator it =
1577 std::ostringstream error_stream;
1578 error_stream <<
"The shared boundary (" << b
1579 <<
") does not exist!!!\n\n";
1580 throw OomphLibError(
1581 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1598 std::map<unsigned, Vector<int>>::iterator it =
1606 std::ostringstream error_stream;
1607 error_stream <<
"The shared boundary (" << b
1608 <<
") does not exist!!!\n\n";
1609 throw OomphLibError(
1610 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1617 std::map<unsigned, Vector<Node*>>::iterator it =
1625 std::ostringstream error_stream;
1626 error_stream <<
"The shared boundary (" << b
1627 <<
") does not exist!!!\n\n";
1628 throw OomphLibError(
1629 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1650 bool node_already_on_this_boundary =
false;
1652 for (
unsigned n = 0; n < nbound_node; n++)
1657 node_already_on_this_boundary =
true;
1662 if (!node_already_on_this_boundary)
1671 std::map<unsigned, Vector<Node*>>::iterator it =
1679 std::ostringstream error_stream;
1680 error_stream <<
"The shared boundary (" << b
1681 <<
") does not exist!!!\n\n";
1682 throw OomphLibError(
1683 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1691 std::map<unsigned, Vector<Node*>>::iterator it =
1696 Vector<Node*>::iterator it_shd_nodes =
1712 std::ostringstream error_stream;
1713 error_stream <<
"The shared boundary (" << b
1714 <<
") does not exist!!!\n\n";
1715 throw OomphLibError(
1716 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1728 std::map<unsigned, Vector<unsigned>>::iterator it =
1733 std::ostringstream error_message;
1735 <<
"The boundary (" << b
1736 <<
") seems not to be shared by any processors,\n"
1737 <<
"it is possible that the boundary was created by the user an not\n"
1738 <<
"automatically by the common interfaces between "
1739 "processors-domains\n";
1740 throw OomphLibError(error_message.str(),
1741 OOMPH_CURRENT_FUNCTION,
1742 OOMPH_EXCEPTION_LOCATION);
1745 return (*it).second;
1757 const unsigned& shd_bnd_id)
1759 std::map<unsigned, unsigned>::iterator it =
1771 const unsigned& shd_bnd_id)
1773 std::map<unsigned, unsigned>::iterator it =
1778 std::ostringstream error_message;
1779 error_message <<
"The shared boundary (" << shd_bnd_id
1780 <<
") does not lie on an internal "
1782 <<
"Make sure to call this method just for shared "
1783 "boundaries that lie "
1784 <<
"on an internal boundary.\n\n";
1785 throw OomphLibError(error_message.str(),
1786 OOMPH_CURRENT_FUNCTION,
1787 OOMPH_EXCEPTION_LOCATION);
1790 return (*it).second;
1796 const unsigned& internal_bnd_id, Vector<unsigned>& shd_bnd_ids)
1799 shd_bnd_ids.clear();
1802 std::map<unsigned, unsigned>::iterator it =
1808 if ((*it).second == internal_bnd_id)
1811 shd_bnd_ids.push_back((*it).first);
1816 if (shd_bnd_ids.size() == 0)
1818 std::ostringstream error_message;
1820 <<
" The internal boundary (" << internal_bnd_id <<
") has no shared "
1821 <<
"boundaries overlapping it\n"
1822 <<
"Make sure to call this method just for internal boundaries that "
1823 <<
"are marked to as being\noverlaped by shared boundaries\n";
1824 throw OomphLibError(error_message.str(),
1825 OOMPH_CURRENT_FUNCTION,
1826 OOMPH_EXCEPTION_LOCATION);
1842 std::map<unsigned, bool>::iterator it;
1850 return (*it).second;
1858 std::map<unsigned, Vector<TriangleMeshPolyLine*>>::iterator it;
1863 std::ostringstream error_message;
1865 <<
"The boundary (" << b
1866 <<
") was marked as been splitted but there\n"
1867 <<
"are not registered polylines to represent the boundary.\n"
1868 <<
"The new polylines were not set up when the boundary was found "
1870 <<
"be splitted or the polylines have been explicitly deleted "
1873 throw OomphLibError(error_message.str(),
1874 OOMPH_CURRENT_FUNCTION,
1875 OOMPH_EXCEPTION_LOCATION);
1878 return (*it).second.size();
1886 std::map<unsigned, Vector<TriangleMeshPolyLine*>>::iterator it;
1890 std::ostringstream error_message;
1892 <<
"The boundary (" << b
1893 <<
") was marked as been splitted but there\n"
1894 <<
"are not registered polylines to represent the boundary.\n"
1895 <<
"The new polylines were not set up when the boundary was found "
1897 <<
"be splitted or the polylines have been explicitly deleted "
1900 throw OomphLibError(error_message.str(),
1901 OOMPH_CURRENT_FUNCTION,
1902 OOMPH_EXCEPTION_LOCATION);
1904 return (*it).second;
1911 const unsigned& isub)
1913 std::map<unsigned, std::vector<bool>>::iterator it;
1922 return (*it).second[isub];
1990 Vector<TriangleMeshPolygon*>& polygons_pt,
1991 Vector<TriangleMeshOpenCurve*>& open_curves_pt);
1996 Vector<TriangleMeshPolyLine*>& unsorted_polylines_pt,
1997 Vector<Vector<TriangleMeshPolyLine*>>& sorted_polylines_pt);
2002 Vector<Vector<TriangleMeshPolyLine*>>& polylines_pt,
2003 Vector<TriangleMeshPolygon*>& polygons_pt);
2009 Vector<Vector<TriangleMeshPolyLine*>>& sorted_open_curves_pt,
2010 Vector<TriangleMeshPolyLine*>& unsorted_shared_to_internal_poly_pt,
2011 Vector<TriangleMeshOpenCurve*>& open_curves_pt);
2023 Vector<Vector<double>>& output_holes_coordinates);
2031 Vector<TriangleMeshPolygon*>& polygons_pt,
2032 Vector<Vector<double>>& output_holes_coordinates);
2040 std::set<FiniteElement*>& element_in_processor_pt,
2041 const int& root_edge_bnd_id,
2042 std::map<std::pair<Node*, Node*>,
bool>& overlapped_face,
2043 std::map<
unsigned, std::map<Node*, bool>>&
2044 node_on_bnd_not_overlapped_by_shd_bnd,
2045 std::list<Node*>& current_polyline_nodes,
2046 std::map<
unsigned, std::list<Node*>>&
2047 shared_bnd_id_to_sorted_list_node_pt,
2048 const unsigned& node_degree,
2050 const bool called_from_load_balance =
false);
2059 OomphCommunicator* comm_pt,
2060 const Vector<unsigned>& element_domain,
2061 const Vector<GeneralisedElement*>& backed_up_el_pt,
2062 const Vector<FiniteElement*>& backed_up_f_el_pt,
2063 std::map<Data*, std::set<unsigned>>& processors_associated_with_data,
2064 const bool& overrule_keep_as_halo_element_status);
2070 const unsigned& nproc,
2071 const Vector<unsigned>& element_domain,
2072 const Vector<GeneralisedElement*>& backed_up_el_pt,
2073 std::map<Data*, std::set<unsigned>>& processors_associated_with_data,
2074 const bool& overrule_keep_as_halo_element_status,
2075 std::map<GeneralisedElement*, unsigned>& element_to_global_index,
2076 Vector<Vector<Vector<GeneralisedElement*>>>& output_halo_elements_pt);
2082 std::map<std::pair<Node*, Node*>,
unsigned>& element_edges_on_boundary);
2090 const Vector<unsigned>& element_domain,
2091 std::map<GeneralisedElement*, unsigned>& element_to_global_index,
2092 std::set<FiniteElement*>& element_in_processor_pt,
2093 Vector<Vector<Vector<GeneralisedElement*>>>& input_halo_elements,
2094 std::map<std::pair<Node*, Node*>,
unsigned>& elements_edges_on_boundary,
2095 Vector<Vector<Vector<TriangleMeshPolyLine*>>>& output_polylines_pt);
2100 const unsigned& initial_shd_bnd_id,
2101 std::list<Node*>& input_nodes,
2102 Vector<FiniteElement*>& input_boundary_element_pt,
2103 Vector<int>& input_face_index_element,
2104 const int& input_connect_to_the_left,
2105 const int& input_connect_to_the_right,
2106 Vector<std::list<Node*>>& output_sorted_nodes_pt,
2107 Vector<Vector<FiniteElement*>>& output_boundary_element_pt,
2108 Vector<Vector<int>>& output_face_index_element,
2109 Vector<int>& output_connect_to_the_left,
2110 Vector<int>& output_connect_to_the_right);
2116 const unsigned& initial_shd_bnd_id,
2117 std::list<Node*>& input_nodes,
2118 Vector<FiniteElement*>& input_boundary_element_pt,
2119 Vector<FiniteElement*>& input_boundary_face_element_pt,
2120 Vector<int>& input_face_index_element,
2121 const int& input_connect_to_the_left,
2122 const int& input_connect_to_the_right,
2123 Vector<std::list<Node*>>& output_sorted_nodes_pt,
2124 Vector<Vector<FiniteElement*>>& output_boundary_element_pt,
2125 Vector<Vector<FiniteElement*>>& output_boundary_face_element_pt,
2126 Vector<Vector<int>>& output_face_index_element,
2127 Vector<int>& output_connect_to_the_left,
2128 Vector<int>& output_connect_to_the_right);
2134 const unsigned& my_rank,
2135 const unsigned& shd_bnd_id,
2136 const unsigned& iproc,
2137 const unsigned& jproc,
2138 std::list<Node*>& sorted_nodes,
2139 const int& root_edge_bnd_id,
2140 Vector<FiniteElement*>& bulk_bnd_ele_pt,
2141 Vector<int>& face_index_ele,
2142 Vector<Vector<TriangleMeshPolyLine*>>& unsorted_polylines_pt,
2143 const int& connect_to_the_left_flag,
2144 const int& connect_to_the_right_flag);
2149 const Vector<unsigned>& target_domain_for_local_non_halo_element)
2151 std::ostringstream error_stream;
2152 error_stream <<
"Empty default load balancing function called.\n";
2153 error_stream <<
"This should be overloaded in a specific "
2154 <<
"RefineableTriangleMesh\n";
2155 throw OomphLibError(
2156 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
2162 Vector<unsigned>& ntmp_boundary_elements,
2163 Vector<Vector<unsigned>>& ntmp_boundary_elements_in_region,
2164 Vector<FiniteElement*>& deleted_elements);
2191 return x < p.
x || (x == p.
x && y < p.
y);
2201 return (A.
x - O.
x) * (B.
y - O.
y) - (A.
y - O.
y) * (B.
x - O.
x);
2209 int n = P.size(), k = 0;
2210 std::vector<Point>
H(2 * n);
2213 std::sort(P.begin(), P.end());
2216 for (
int i = 0; i < n; ++i)
2218 while (k >= 2 &&
cross(
H[k - 2],
H[k - 1], P[i]) <= 0) k--;
2223 for (
int i = n - 2, t = k + 1; i >= 0; i--)
2225 while (k >= t &&
cross(
H[k - 2],
H[k - 1], P[i]) <= 0) k--;
2246 template<
class ELEMENT>
2248 public virtual RefineableMeshBase
2255 typedef void (*MeshUpdateFctPt)(Mesh* mesh_pt);
2263 typedef void (*InternalHolePointUpdateFctPt)(
const unsigned& ihole,
2264 TriangleMeshPolygon* poly_pt);
2266 #ifdef OOMPH_HAS_TRIANGLE_LIB
2272 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2273 :
TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
2276 initialise_adaptation_data();
2279 initialise_boundary_refinement_data();
2286 const std::string& node_file_name,
2287 const std::string& element_file_name,
2288 const std::string& poly_file_name,
2289 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
2290 const bool& allow_automatic_creation_of_vertices_on_boundaries =
true)
2296 allow_automatic_creation_of_vertices_on_boundaries)
2300 create_polylines_from_polyfiles(node_file_name, poly_file_name);
2303 initialise_adaptation_data();
2306 initialise_boundary_refinement_data();
2310 #ifdef OOMPH_HAS_TRIANGLE_LIB
2317 const Vector<double>& target_area,
2318 TriangulateIO& triangulate_io,
2319 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
2320 const bool& use_attributes =
false,
2321 const bool& allow_automatic_creation_of_vertices_on_boundaries =
true,
2322 OomphCommunicator* comm_pt = 0)
2325 initialise_adaptation_data();
2328 initialise_boundary_refinement_data();
2331 this->Time_stepper_pt = time_stepper_pt;
2334 TriangulateIO triangle_refine;
2337 TriangleHelper::initialise_triangulateio(this->Triangulateio);
2340 this->Triangulateio_exists =
true;
2343 this->refine_triangulateio(triangulate_io, target_area, triangle_refine);
2346 std::stringstream input_string_stream;
2347 input_string_stream <<
"-pq30-ra";
2350 if (!allow_automatic_creation_of_vertices_on_boundaries)
2352 input_string_stream <<
" -YY";
2356 this->Allow_automatic_creation_of_vertices_on_boundaries =
2357 allow_automatic_creation_of_vertices_on_boundaries;
2360 this->Use_attributes = use_attributes;
2363 char triswitches[100];
2364 sprintf(triswitches,
"%s", input_string_stream.str().c_str());
2367 triangulate(triswitches, &triangle_refine, &this->Triangulateio, 0);
2369 #ifdef OOMPH_HAS_FPUCONTROLH
2371 fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
2376 this->Tmp_mesh_pt =
new TriangleScaffoldMesh(this->Triangulateio);
2379 this->build_from_scaffold(time_stepper_pt, use_attributes);
2382 delete this->Tmp_mesh_pt;
2383 this->Tmp_mesh_pt = 0;
2386 bool clear_hole_data =
false;
2387 TriangleHelper::clear_triangulateio(triangle_refine, clear_hole_data);
2389 #ifdef OOMPH_HAS_MPI
2395 this->set_communicator_pt(comm_pt);
2400 unsigned nb = nboundary();
2401 for (
unsigned b = 0; b < nb; b++)
2403 this->
template setup_boundary_coordinates<ELEMENT>(b);
2417 Print_timings_transfering_target_areas =
true;
2424 Print_timings_transfering_target_areas =
false;
2430 Disable_projection =
false;
2436 Disable_projection =
true;
2442 Print_timings_projection =
true;
2448 Print_timings_projection =
false;
2456 return Nbin_x_for_area_transfer;
2464 return Nbin_y_for_area_transfer;
2473 return Max_sample_points_for_limited_locate_zeta_during_target_area_transfer;
2479 return Max_element_size;
2485 return Min_element_size;
2491 return Min_permitted_angle;
2498 return Use_iterative_solver_for_projection;
2505 Use_iterative_solver_for_projection =
true;
2512 Use_iterative_solver_for_projection =
false;
2518 set_print_level_timings_adaptation(print_level);
2524 Print_timings_level_adaptation = 0;
2530 const unsigned max_print_level = 3;
2532 if (print_level > max_print_level)
2534 Print_timings_level_adaptation = max_print_level;
2538 Print_timings_level_adaptation = print_level;
2545 set_print_level_timings_load_balance(print_level);
2551 Print_timings_level_load_balance = 0;
2557 const unsigned max_print_level = 3;
2559 if (print_level > max_print_level)
2561 Print_timings_level_load_balance = max_print_level;
2565 Print_timings_level_load_balance = print_level;
2572 outfile << std::endl;
2573 outfile <<
"Targets for mesh adaptation: " << std::endl;
2574 outfile <<
"---------------------------- " << std::endl;
2575 outfile <<
"Target for max. error: " << Max_permitted_error << std::endl;
2576 outfile <<
"Target for min. error: " << Min_permitted_error << std::endl;
2577 outfile <<
"Target min angle: " << Min_permitted_angle << std::endl;
2578 outfile <<
"Min. allowed element size: " << Min_element_size << std::endl;
2579 outfile <<
"Max. allowed element size: " << Max_element_size << std::endl;
2580 outfile <<
"Don't unrefine if less than " << Max_keep_unrefined
2581 <<
" elements need unrefinement." << std::endl;
2582 outfile << std::endl;
2589 unsigned nelem = nelement();
2590 Vector<double> elem_error(nelem, DBL_MAX);
2593 double backup = Min_element_size;
2596 double orig_max_area, orig_min_area;
2597 this->max_and_min_element_size(orig_max_area, orig_min_area);
2600 Min_element_size = orig_min_area / 3.0;
2606 Min_element_size = backup;
2614 throw OomphLibError(
"unrefine_uniformly() not implemented yet",
2615 OOMPH_CURRENT_FUNCTION,
2616 OOMPH_EXCEPTION_LOCATION);
2622 void adapt(
const Vector<double>& elem_error);
2630 return Mesh_update_fct_pt;
2638 return Internal_hole_point_update_fct_pt;
2642 #ifdef OOMPH_HAS_MPI
2645 std::map<unsigned, Vector<Node*>>::iterator it =
2646 Sorted_shared_boundary_node_pt.find(b);
2647 if (it == Sorted_shared_boundary_node_pt.end())
2649 std::ostringstream error_message;
2650 error_message <<
"The boundary (" << b <<
") is not marked as shared\n";
2651 throw OomphLibError(error_message.str(),
2652 OOMPH_CURRENT_FUNCTION,
2653 OOMPH_EXCEPTION_LOCATION);
2655 return (*it).second.size();
2660 Sorted_shared_boundary_node_pt.clear();
2665 std::map<unsigned, Vector<Node*>>::iterator it =
2666 Sorted_shared_boundary_node_pt.find(b);
2667 if (it == Sorted_shared_boundary_node_pt.end())
2669 std::ostringstream error_message;
2670 error_message <<
"The boundary (" << b <<
") is not marked as shared\n";
2671 throw OomphLibError(error_message.str(),
2672 OOMPH_CURRENT_FUNCTION,
2673 OOMPH_EXCEPTION_LOCATION);
2675 return (*it).second[i];
2681 std::map<unsigned, Vector<Node*>>::iterator it =
2682 Sorted_shared_boundary_node_pt.find(b);
2683 if (it == Sorted_shared_boundary_node_pt.end())
2685 std::ostringstream error_message;
2686 error_message <<
"The boundary (" << b <<
") is not marked as shared\n";
2687 throw OomphLibError(error_message.str(),
2688 OOMPH_CURRENT_FUNCTION,
2689 OOMPH_EXCEPTION_LOCATION);
2691 return (*it).second;
2698 void create_polylines_from_polyfiles(
const std::string& node_file_name,
2699 const std::string& poly_file_name);
2701 #ifdef OOMPH_HAS_MPI
2704 void fill_boundary_elements_and_nodes_for_internal_boundaries();
2710 void fill_boundary_elements_and_nodes_for_internal_boundaries(
2711 std::ofstream& outfile);
2716 std::istream& restart_file)
2719 if (this->is_mesh_distributed())
2723 this->fill_boundary_elements_and_nodes_for_internal_boundaries();
2727 this->reset_shared_boundary_elements_and_nodes(comm_pt);
2731 this->sort_nodes_on_shared_boundaries();
2734 this->reset_halo_haloed_scheme();
2737 const unsigned noriginal_boundaries =
2738 this->initial_shared_boundary_id();
2739 this->set_nboundary(noriginal_boundaries);
2743 for (
unsigned b = 0; b < noriginal_boundaries; b++)
2748 this->identify_boundary_segments_and_assign_initial_zeta_values(b,
2751 if (this->boundary_geom_object_pt(b) != 0)
2754 this->
template setup_boundary_coordinates<ELEMENT>(b);
2759 this->snap_nodes_onto_geometric_objects();
2769 #ifdef OOMPH_HAS_MPI
2772 unsigned my_rank = 0;
2773 if (this->is_mesh_distributed())
2775 my_rank = this->communicator_pt()->my_rank();
2782 const unsigned ninternal = this->Internal_polygon_pt.size();
2783 for (
unsigned i_internal = 0; i_internal < ninternal; i_internal++)
2785 this->update_polygon_after_restart(
2786 this->Internal_polygon_pt[i_internal]);
2790 const unsigned nouter = this->Outer_boundary_pt.size();
2791 for (
unsigned i_outer = 0; i_outer < nouter; i_outer++)
2793 this->update_polygon_after_restart(this->Outer_boundary_pt[i_outer]);
2796 #ifdef OOMPH_HAS_MPI
2799 if (this->is_mesh_distributed())
2801 const unsigned ncurves = this->nshared_boundary_curves(my_rank);
2802 for (
unsigned nc = 0; nc < ncurves; nc++)
2805 this->update_shared_curve_after_restart(
2806 this->Shared_boundary_polyline_pt[my_rank][nc]
2813 const unsigned n_open_polyline = this->Internal_open_curve_pt.size();
2814 for (
unsigned i = 0; i < n_open_polyline; i++)
2816 this->update_open_curve_after_restart(this->Internal_open_curve_pt[i]);
2820 #ifdef OOMPH_HAS_MPI
2825 const Vector<unsigned>& input_target_domain_for_local_non_halo_element);
2830 void get_shared_boundary_elements_and_face_indexes(
2831 const Vector<FiniteElement*>& first_element_pt,
2832 const Vector<FiniteElement*>& second_element_pt,
2833 Vector<FiniteElement*>& first_shared_boundary_element_pt,
2834 Vector<unsigned>& first_shared_boundary_element_face_index,
2835 Vector<FiniteElement*>& second_shared_boundary_element_pt,
2836 Vector<unsigned>& second_shared_boundary_element_face_index);
2841 void create_new_shared_boundaries(
2842 std::set<FiniteElement*>& element_in_processor_pt,
2843 Vector<Vector<FiniteElement*>>& new_shared_boundary_element_pt,
2844 Vector<Vector<unsigned>>& new_shared_boundary_element_face_index);
2849 void compute_shared_node_degree_helper(
2850 Vector<Vector<FiniteElement*>>& unsorted_face_ele_pt,
2851 std::map<Node*, unsigned>& global_node_degree);
2857 void create_adjacency_matrix_new_shared_edges_helper(
2858 Vector<Vector<FiniteElement*>>& unsorted_face_ele_pt,
2859 Vector<Vector<Node*>>& tmp_sorted_shared_node_pt,
2860 std::map<Node*, Vector<Vector<unsigned>>>& node_alias,
2861 Vector<Vector<Vector<unsigned>>>& adjacency_matrix);
2865 void get_shared_boundary_segment_nodes_helper(
2866 const unsigned& shd_bnd_id, Vector<Vector<Node*>>& tmp_segment_nodes);
2874 void get_boundary_segment_nodes_helper(
2875 const unsigned& b, Vector<Vector<Node*>>& tmp_segment_nodes);
2881 Do_boundary_unrefinement_constrained_by_target_areas =
true;
2886 Do_boundary_unrefinement_constrained_by_target_areas =
false;
2891 Do_boundary_refinement_constrained_by_target_areas =
true;
2896 Do_boundary_refinement_constrained_by_target_areas =
false;
2903 Do_shared_boundary_unrefinement_constrained_by_target_areas =
true;
2908 Do_shared_boundary_unrefinement_constrained_by_target_areas =
false;
2913 Do_shared_boundary_refinement_constrained_by_target_areas =
true;
2918 Do_shared_boundary_refinement_constrained_by_target_areas =
false;
2934 std::set<Vector<double>>& vertices)
2937 std::map<unsigned, std::set<Vector<double>>>::iterator it =
2938 Boundary_connections_pt.find(b);
2940 if (it != Boundary_connections_pt.end())
2943 vertices = (*it).second;
2955 const void synchronize_shared_boundary_connections();
2960 void add_vertices_for_non_deletion();
2965 void add_non_delete_vertices_from_boundary_helper(
2966 Vector<Vector<Node*>> src_bound_segment_node_pt,
2967 Vector<Vector<Node*>> dst_bound_segment_node_pt,
2968 const unsigned& dst_bnd_id,
2969 const unsigned& dst_bnd_chunk);
2974 void create_temporary_boundary_connections(
2975 Vector<TriangleMeshPolygon*>& tmp_outer_polygons_pt,
2976 Vector<TriangleMeshOpenCurve*>& tmp_open_curves_pt);
2983 void restore_boundary_connections(
2984 Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
2985 Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
2993 void restore_polyline_connections_helper(
2994 TriangleMeshPolyLine* polyline_pt,
2995 Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
2996 Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
3004 void resume_boundary_connections(
3005 Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
3006 Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
3010 bool get_connected_vertex_number_on_dst_boundary(
3011 Vector<double>& vertex_coordinates,
3012 const unsigned& dst_b_id,
3013 unsigned& vertex_number);
3021 bool unrefine_boundary(
const unsigned& b,
3023 Vector<Vector<double>>& vector_bnd_vertices,
3024 double& unrefinement_tolerance,
3025 const bool& check_only =
false);
3033 bool refine_boundary(Mesh* face_mesh_pt,
3034 Vector<Vector<double>>& vector_bnd_vertices,
3035 double& refinement_tolerance,
3036 const bool& check_only =
false);
3042 bool apply_max_length_constraint(
3044 Vector<Vector<double>>& vector_bnd_vertices,
3045 double& max_length_constraint);
3052 bool unrefine_boundary_constrained_by_target_area(
3055 Vector<Vector<double>>& vector_bnd_vertices,
3056 double& unrefinement_tolerance,
3057 Vector<double>& area_constraint);
3064 bool refine_boundary_constrained_by_target_area(
3065 MeshAsGeomObject* mesh_geom_obj_pt,
3066 Vector<Vector<double>>& vector_bnd_vertices,
3067 double& refinement_tolerance,
3068 Vector<double>& area_constraint);
3075 bool unrefine_shared_boundary_constrained_by_target_area(
3078 Vector<Vector<double>>& vector_bnd_vertices,
3079 Vector<double>& area_constraint);
3086 bool refine_shared_boundary_constrained_by_target_area(
3087 Vector<Vector<double>>& vector_bnd_vertices,
3088 Vector<double>& area_constraint);
3107 Do_boundary_unrefinement_constrained_by_target_areas =
true;
3108 Do_boundary_refinement_constrained_by_target_areas =
true;
3109 Do_shared_boundary_unrefinement_constrained_by_target_areas =
true;
3110 Do_shared_boundary_refinement_constrained_by_target_areas =
true;
3113 #ifdef OOMPH_HAS_MPI
3123 void sort_nodes_on_shared_boundaries();
3128 void reset_shared_boundary_elements_and_nodes(
3129 const bool flush_elements =
true,
3130 const bool update_elements =
true,
3131 const bool flush_nodes =
true,
3132 const bool update_nodes =
true);
3138 void reset_halo_haloed_scheme();
3147 void compute_global_node_names_and_shared_nodes(
3148 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3149 other_proc_shd_bnd_node_pt,
3150 Vector<Vector<Vector<unsigned>>>& global_node_names,
3151 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3152 Vector<Node*>& global_shared_node_pt);
3158 void send_boundary_node_info_of_shared_nodes(
3159 Vector<Vector<Vector<unsigned>>>& global_node_names,
3160 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3161 Vector<Node*>& global_shared_node_pt);
3166 void reset_halo_haloed_scheme_helper(
3167 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3168 other_proc_shd_bnd_node_pt,
3169 Vector<Vector<Node*>>& iproc_currently_created_nodes_pt,
3170 Vector<Vector<Vector<unsigned>>>& global_node_names,
3171 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3172 Vector<Node*>& global_shared_node_pt);
3187 Vector<FiniteElement*>& new_elements_on_domain, FiniteElement*& ele_pt)
3190 const unsigned nnew_elements_on_domain = new_elements_on_domain.size();
3193 bool already_on_new_domain =
false;
3194 unsigned new_domain_ele_index = 0;
3196 for (
unsigned e = 0; e < nnew_elements_on_domain; e++)
3198 if (ele_pt == new_elements_on_domain[e])
3201 already_on_new_domain =
true;
3203 new_domain_ele_index = e;
3209 if (!already_on_new_domain)
3212 new_elements_on_domain.push_back(ele_pt);
3214 return nnew_elements_on_domain;
3219 return new_domain_ele_index;
3227 void get_required_elemental_information_load_balance_helper(
3229 Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
3230 FiniteElement* ele_pt);
3238 const unsigned nnew_nodes_on_domain = new_nodes_on_domain.size();
3241 bool already_on_new_domain =
false;
3242 unsigned new_domain_node_index = 0;
3244 for (
unsigned n = 0; n < nnew_nodes_on_domain; n++)
3246 if (node_pt == new_nodes_on_domain[n])
3249 already_on_new_domain =
true;
3251 new_domain_node_index = n;
3257 if (!already_on_new_domain)
3260 new_nodes_on_domain.push_back(node_pt);
3262 return nnew_nodes_on_domain;
3267 return new_domain_node_index;
3272 void add_node_load_balance_helper(
3274 Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
3275 Vector<Node*>& new_nodes_on_domain,
3286 void get_required_nodal_information_load_balance_helper(
3287 Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
3294 void create_element_load_balance_helper(
3296 Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
3297 Vector<Vector<std::map<unsigned, FiniteElement*>>>&
3298 received_old_haloed_element_pt,
3299 Vector<FiniteElement*>& new_elements_on_domain,
3300 Vector<Node*>& new_nodes_on_domain,
3301 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3302 other_proc_shd_bnd_node_pt,
3303 Vector<Vector<Vector<unsigned>>>& global_node_names,
3304 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3305 Vector<Node*>& global_shared_node_pt);
3312 void add_element_load_balance_helper(
3313 const unsigned& iproc,
3314 Vector<Vector<std::map<unsigned, FiniteElement*>>>&
3315 received_old_haloed_element_pt,
3316 FiniteElement* ele_pt);
3319 void add_received_node_load_balance_helper(
3321 Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
3322 Vector<Vector<std::map<unsigned, FiniteElement*>>>&
3323 received_old_haloed_element_pt,
3324 Vector<Node*>& new_nodes_on_domain,
3325 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3326 other_proc_shd_bnd_node_pt,
3328 unsigned& node_index,
3329 FiniteElement*
const& new_el_pt,
3330 Vector<Vector<Vector<unsigned>>>& global_node_names,
3331 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3332 Vector<Node*>& global_shared_node_pt);
3337 void construct_new_node_load_balance_helper(
3339 Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
3340 Vector<Vector<std::map<unsigned, FiniteElement*>>>&
3341 received_old_haloed_element_pt,
3342 Vector<Node*>& new_nodes_on_domain,
3343 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3344 other_proc_shd_bnd_node_pt,
3346 unsigned& node_index,
3347 FiniteElement*
const& new_el_pt,
3348 Vector<Vector<Vector<unsigned>>>& global_node_names,
3349 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3350 Vector<Node*>& global_shared_node_pt);
3375 #ifdef ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
3392 GeneralisedElement*& el_pt)
3395 unsigned n_haloed = this->nroot_haloed_element(p);
3398 bool already_haloed_element =
false;
3399 unsigned haloed_el_index = 0;
3400 for (
unsigned eh = 0; eh < n_haloed; eh++)
3402 if (el_pt == this->root_haloed_element_pt(p, eh))
3405 already_haloed_element =
true;
3407 haloed_el_index = eh;
3413 if (!already_haloed_element)
3416 this->add_root_haloed_element_pt(p, el_pt);
3423 return haloed_el_index;
3432 unsigned n_haloed_nod = this->nhaloed_node(p);
3435 bool is_an_haloed_node =
false;
3436 unsigned haloed_node_index = 0;
3437 for (
unsigned k = 0; k < n_haloed_nod; k++)
3439 if (nod_pt == this->haloed_node_pt(p, k))
3441 is_an_haloed_node =
true;
3442 haloed_node_index = k;
3448 if (!is_an_haloed_node)
3451 this->add_haloed_node_pt(p, nod_pt);
3453 return n_haloed_nod;
3458 return haloed_node_index;
3465 void get_required_elemental_information_helper(
unsigned& iproc,
3466 FiniteElement* ele_pt);
3471 void get_required_nodal_information_helper(
unsigned& iproc, Node* nod_pt);
3474 void add_haloed_node_helper(
unsigned& iproc, Node* nod_pt);
3477 void send_and_receive_elements_nodes_info(
int& send_proc,
int& recv_proc);
3481 void create_halo_element(
3483 Vector<Node*>& new_nodes_on_domain,
3484 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3485 other_proc_shd_bnd_node_pt,
3486 Vector<Vector<Vector<unsigned>>>& global_node_names,
3487 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3488 Vector<Node*>& global_shared_node_pt);
3494 void add_halo_element_helper(
unsigned& iproc, FiniteElement* ele_pt);
3497 void add_halo_node_helper(
3499 Vector<Node*>& new_nodes_on_domain,
3500 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3501 other_proc_shd_bnd_node_pt,
3503 unsigned& node_index,
3504 FiniteElement*
const& new_el_pt,
3505 Vector<Vector<Vector<unsigned>>>& global_node_names,
3506 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3507 Vector<Node*>& global_shared_node_pt);
3511 void construct_new_halo_node_helper(
3513 Vector<Node*>& new_nodes_on_domain,
3514 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3515 other_proc_shd_bnd_node_pt,
3517 unsigned& node_index,
3518 FiniteElement*
const& new_el_pt,
3519 Vector<Vector<Vector<unsigned>>>& global_node_names,
3520 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3521 Vector<Node*>& global_shared_node_pt);
3527 void update_other_proc_shd_bnd_node_helper(
3529 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3530 other_proc_shd_bnd_node_pt,
3531 Vector<unsigned>& other_processor_1,
3532 Vector<unsigned>& other_processor_2,
3533 Vector<unsigned>& other_shared_boundaries,
3534 Vector<unsigned>& other_indexes,
3535 Vector<Vector<Vector<unsigned>>>& global_node_names,
3536 std::map<Vector<unsigned>,
unsigned>& node_name_to_global_index,
3537 Vector<Node*>& global_shared_node_pt);
3552 bool update_polygon_using_face_mesh(TriangleMeshPolygon* polygon_pt,
3553 const bool& check_only =
false);
3562 bool update_open_curve_using_face_mesh(
3563 TriangleMeshOpenCurve* open_polyline_pt,
const bool& check_only =
false);
3570 virtual bool surface_remesh_for_inner_hole_boundaries(
3571 Vector<Vector<double>>& internal_point_coord,
3572 const bool& check_only =
false);
3582 void create_unsorted_face_mesh_representation(
const unsigned& boundary_id,
3583 Mesh* face_mesh_pt);
3589 void create_sorted_face_mesh_representation(
3590 const unsigned& boundary_id,
3592 std::map<FiniteElement*, bool>& is_inverted,
3593 bool& inverted_face_mesh);
3598 void get_face_mesh_representation(TriangleMeshPolygon* polygon_pt,
3599 Vector<Mesh*>& face_mesh_pt);
3603 void get_face_mesh_representation(TriangleMeshOpenCurve* open_polyline_pt,
3604 Vector<Mesh*>& face_mesh_pt);
3607 void update_polygon_after_restart(TriangleMeshPolygon*& polygon_pt);
3610 void update_open_curve_after_restart(TriangleMeshOpenCurve*& open_curve_pt);
3614 bool update_polygon_using_elements_area(TriangleMeshPolygon*& polygon_pt,
3615 const Vector<double>& target_area);
3619 bool update_open_curve_using_elements_area(
3620 TriangleMeshOpenCurve*& open_curve_pt,
const Vector<double>& target_area);
3622 #ifdef OOMPH_HAS_MPI
3625 bool update_shared_curve_using_elements_area(
3626 Vector<TriangleMeshPolyLine*>& vector_polyline_pt,
3627 const Vector<double>& target_areas);
3630 void update_shared_curve_after_restart(
3631 Vector<TriangleMeshPolyLine*>& vector_polyline_pt);
3641 this->Nbin_x_for_area_transfer = 100;
3646 this->Nbin_y_for_area_transfer = 100;
3650 Max_sample_points_for_limited_locate_zeta_during_target_area_transfer = 5;
3653 this->Max_element_size = 1.0;
3654 this->Min_element_size = 0.001;
3655 this->Min_permitted_angle = 15.0;
3658 this->Disable_projection =
false;
3661 this->Use_iterative_solver_for_projection =
true;
3664 this->Print_timings_level_adaptation = 0;
3667 this->Print_timings_level_load_balance = 0;
3671 this->Print_timings_transfering_target_areas =
false;
3674 this->Print_timings_projection =
false;
3680 Mesh_update_fct_pt = 0;
3684 Internal_hole_point_update_fct_pt = 0;
3687 #ifdef OOMPH_HAS_TRIANGLE_LIB
3691 void refine_triangulateio(TriangulateIO& triangulate_io,
3692 const Vector<double>& target_area,
3693 TriangulateIO& triangle_refine);
3700 Vector<double>& target_area)
3702 double min_angle = DBL_MAX;
3703 unsigned count_unrefined = 0;
3704 unsigned count_refined = 0;
3705 this->Nrefinement_overruled = 0;
3708 std::map<FiniteElement*, double> max_area_from_region;
3709 for (std::map<unsigned, double>::iterator it =
3710 this->Regions_areas.begin();
3711 it != this->Regions_areas.end();
3714 unsigned r = (*it).first;
3715 unsigned nel = this->nregion_element(r);
3716 for (
unsigned e = 0; e < nel; e++)
3718 max_area_from_region[this->region_element_pt(r, e)] = (*it).second;
3722 unsigned nel = this->nelement();
3723 for (
unsigned e = 0; e < nel; e++)
3726 FiniteElement* el_pt = this->finite_element_pt(e);
3729 double area = el_pt->size();
3733 double ax = el_pt->node_pt(0)->x(0);
3734 double ay = el_pt->node_pt(0)->x(1);
3736 double bx = el_pt->node_pt(1)->x(0);
3737 double by = el_pt->node_pt(1)->x(1);
3739 double cx = el_pt->node_pt(2)->x(0);
3740 double cy = el_pt->node_pt(2)->x(1);
3744 acos(((ax - cx) * (bx - cx) + (ay - cy) * (by - cy)) /
3745 (sqrt((ax - cx) * (ax - cx) + (ay - cy) * (ay - cy)) *
3746 sqrt((bx - cx) * (bx - cx) + (by - cy) * (by - cy)))) *
3747 180.0 / MathematicalConstants::Pi;
3748 min_angle = std::min(min_angle, angle0);
3751 acos(((ax - bx) * (cx - bx) + (ay - by) * (cy - by)) /
3752 (sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by)) *
3753 sqrt((cx - bx) * (cx - bx) + (cy - by) * (cy - by)))) *
3754 180.0 / MathematicalConstants::Pi;
3755 min_angle = std::min(min_angle, angle1);
3757 double angle2 = 180.0 - angle0 - angle1;
3758 min_angle = std::min(min_angle, angle2);
3763 double size_ratio = 3.0;
3764 if (elem_error[e] > max_permitted_error())
3767 target_area[e] = std::max(area / size_ratio, Min_element_size);
3770 target_area[e] = std::min(target_area[e], Max_element_size);
3772 if (target_area[e] != Min_element_size)
3778 this->Nrefinement_overruled++;
3781 else if (elem_error[e] < min_permitted_error())
3784 target_area[e] = std::min(size_ratio * area, Max_element_size);
3787 target_area[e] = std::max(target_area[e], Min_element_size);
3789 if (target_area[e] != Max_element_size)
3797 double area_leave_alone = std::max(area, Min_element_size);
3798 target_area[e] = std::min(area_leave_alone, Max_element_size);
3802 std::map<FiniteElement*, double>::iterator it =
3803 max_area_from_region.find(el_pt);
3804 if (it != max_area_from_region.end())
3806 target_area[e] = std::min(target_area[e], (*it).second);
3812 this->Nrefined = count_refined;
3813 this->Nunrefined = count_unrefined;
3815 if (this->Nrefinement_overruled != 0)
3818 <<
"\nNOTE: Refinement of " << this->Nrefinement_overruled
3820 <<
"overruled \nbecause the target area would have "
3821 <<
"been below \nthe minimum permitted area of " << Min_element_size
3822 <<
".\nYou can change the minimum permitted area with the\n"
3823 <<
"function RefineableTriangleMesh::min_element_size().\n\n";
3893 template<
class ELEMENT>
3895 public virtual SolidMesh
3898 #ifdef OOMPH_HAS_TRIANGLE_LIB
3904 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
3905 :
TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
3908 set_lagrangian_nodal_coordinates();
3914 const std::string& node_file_name,
3915 const std::string& element_file_name,
3916 const std::string& poly_file_name,
3917 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
3918 const bool& allow_automatic_creation_of_vertices_on_boundaries =
true)
3924 allow_automatic_creation_of_vertices_on_boundaries)
3927 set_lagrangian_nodal_coordinates();
3939 #ifdef OOMPH_HAS_TRIANGLE_LIB
3944 template<
class ELEMENT>
3947 public virtual SolidMesh
3954 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
3955 :
TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt),
3960 set_lagrangian_nodal_coordinates();
3966 const Vector<double>& target_area,
3967 TriangulateIO& triangulate_io,
3968 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper,
3969 const bool& use_attributes =
false,
3970 const bool& allow_automatic_creation_of_vertices_on_boundaries =
true,
3971 OomphCommunicator* comm_pt = 0)
3977 allow_automatic_creation_of_vertices_on_boundaries,
3981 set_lagrangian_nodal_coordinates();
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
RefineableSolidTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameter.
virtual ~RefineableSolidTriangleMesh()
Empty Destructor.
RefineableSolidTriangleMesh(const Vector< double > &target_area, TriangulateIO &triangulate_io, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false, const bool &allow_automatic_creation_of_vertices_on_boundaries=true, OomphCommunicator *comm_pt=0)
Build mesh from specified triangulation and associated target areas for elements in it.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
unsigned & nbin_y_for_area_transfer()
Read/write access to number of bins in the y-direction when transferring target areas by bin method....
void disable_iterative_solver_for_projection()
Enables the use of an iterative solver for the projection problem.
std::map< unsigned, std::set< Vector< double > > > Boundary_connections_pt
A map that stores the vertices that receive connections, they are identified by the boundary number t...
MeshUpdateFctPt & mesh_update_fct_pt()
Access to function pointer to function that updates the mesh following the snapping of boundary nodes...
unsigned Print_timings_level_adaptation
The printing level for adaptation.
unsigned try_to_add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Check if necessary to add the node as haloed or if it has been previously added to the haloed scheme.
double compute_area_target(const Vector< double > &elem_error, Vector< double > &target_area)
Compute target area based on the element's error and the error target; return minimum angle (in degre...
void enable_timings_tranfering_target_areas()
Enables info. and timings for tranferring of target areas.
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors.
virtual ~RefineableTriangleMesh()
Empty Destructor.
void set_print_level_timings_load_balance(const unsigned &print_level)
Sets the printing level of timings for load balance.
void disable_boundary_unrefinement_constrained_by_target_areas()
void initialise_adaptation_data()
Helper function to initialise data associated with adaptation.
void enable_print_timings_load_balance(const unsigned &print_level=1)
Enables printing of timings for load balance.
void update_polyline_representation_from_restart()
Method used to update the polylines representation after restart.
void initialise_boundary_refinement_data()
Set all the flags to true (the default values)
void enable_boundary_unrefinement_constrained_by_target_areas()
Enable/disable unrefinement/refinement methods for original boundaries.
void enable_timings_projection()
Enables info. and timings for projection.
bool Do_shared_boundary_unrefinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
double Max_element_size
Max permitted element size.
void set_print_level_timings_adaptation(const unsigned &print_level)
Sets the printing level of timings for adaptation.
void disable_print_timings_load_balance()
Disables printing of timings for load balance.
unsigned try_to_add_node_pt_load_balance(Vector< Node * > &new_nodes_on_domain, Node *&node_pt)
Check if necessary to add the node to the new domain or if it has been already added.
unsigned Nbin_x_for_area_transfer
Number of bins in the x-direction when transferring target areas by bin method. Only used if we don't...
void enable_shared_boundary_unrefinement_constrained_by_target_areas()
Enable/disable unrefinement/refinement methods for shared boundaries.
void reestablish_distribution_info_for_restart(OomphCommunicator *comm_pt, std::istream &restart_file)
Used to re-establish any additional info. related with the distribution after a re-starting for trian...
unsigned try_to_add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Check if necessary to add the element as haloed or if it has been previously added to the haloed sche...
void disable_timings_projection()
Disables info. and timings for projection.
void refine_uniformly(DocInfo &doc_info)
Refine mesh uniformly and doc process.
RefineableTriangleMesh(const Vector< double > &target_area, TriangulateIO &triangulate_io, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false, const bool &allow_automatic_creation_of_vertices_on_boundaries=true, OomphCommunicator *comm_pt=0)
Build mesh from specified triangulation and associated target areas for elements in it NOTE: This is ...
MeshUpdateFctPt Mesh_update_fct_pt
Function pointer to function that updates the mesh following the snapping of boundary nodes to the bo...
void enable_print_timings_adaptation(const unsigned &print_level=1)
Enables printing of timings for adaptation.
void disable_boundary_refinement_constrained_by_target_areas()
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds.
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles.
double & max_element_size()
Max element size allowed during adaptation.
InternalHolePointUpdateFctPt Internal_hole_point_update_fct_pt
Function pointer to function that can be set to update the position of the central point in internal ...
void disable_projection()
Disables the solution projection step during adaptation.
unsigned try_to_add_element_pt_load_balance(Vector< FiniteElement * > &new_elements_on_domain, FiniteElement *&ele_pt)
Check if necessary to add the element to the new domain or if it has been previously added.
void enable_boundary_refinement_constrained_by_target_areas()
bool Use_iterative_solver_for_projection
Flag to indicate whether to use or not an iterative solver (CG with diagonal preconditioned) for the ...
unsigned max_sample_points_for_limited_locate_zeta_during_target_area_transfer()
Read/write access to number of sample points from which we try to locate zeta by Newton method when t...
unsigned & nbin_x_for_area_transfer()
Read/write access to number of bins in the x-direction when transferring target areas by bin method....
void disable_shared_boundary_unrefinement_constrained_by_target_areas()
const bool boundary_connections(const unsigned &b, const unsigned &c, std::set< Vector< double >> &vertices)
Verifies if the given boundary receives a connection, and if that is the case then returns the list o...
RefineableTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
Build mesh, based on the polyfiles.
bool Disable_projection
Enable/disable solution projection during adaptation.
double Min_permitted_angle
Min angle before remesh gets triggered.
void disable_shared_boundary_refinement_constrained_by_target_areas()
RefineableTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
Node * sorted_shared_boundary_node_pt(unsigned &b, unsigned &i)
bool Do_boundary_unrefinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
double & min_permitted_angle()
Min angle before remesh gets triggered.
void doc_adaptivity_targets(std::ostream &outfile)
Doc the targets for mesh adaptation.
unsigned Max_sample_points_for_limited_locate_zeta_during_target_area_transfer
Default value for max. number of sample points used for locate_zeta when transferring target areas us...
Vector< Node * > sorted_shared_boundary_node_pt(unsigned &b)
void enable_projection()
Enables the solution projection step during adaptation.
void enable_iterative_solver_for_projection()
Enables the use of an iterative solver for the projection problem.
void disable_print_timings_adaptation()
Disables printing of timings for adaptation.
bool Do_boundary_refinement_constrained_by_target_areas
Flag that enables or disables boundary refinement (true by default)
std::map< unsigned, Vector< Node * > > Sorted_shared_boundary_node_pt
Stores the nodes in the boundaries in the same order in all the processors Sorted_shared_boundary_nod...
double Min_element_size
Min permitted element size.
bool Do_shared_boundary_refinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
bool Print_timings_projection
Enable/disable printing timings for projection.
void disable_timings_tranfering_target_areas()
Disables info. and timings for tranferring of target areas.
bool Print_timings_transfering_target_areas
Enable/disable printing timings for transfering target areas.
void flush_sorted_shared_boundary_node()
double & min_element_size()
Min element size allowed during adaptation.
unsigned Print_timings_level_load_balance
The printing level for load balance.
unsigned nsorted_shared_boundary_node(unsigned &b)
void enable_shared_boundary_refinement_constrained_by_target_areas()
bool use_iterative_solver_for_projection()
unsigned unrefine_uniformly()
Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarses...
unsigned Nbin_y_for_area_transfer
Number of bins in the y-direction when transferring target areas by bin method. Only used if we don't...
InternalHolePointUpdateFctPt & internal_hole_point_update_fct_pt()
Access to function pointer to can be used to generate the internal point for the ihole-th hole.
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
Vector< std::string > Flat_packed_unsigneds_string
Temporary vector of strings to enable full annotation of RefineableTriangleMesh comms.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
SolidTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on closed curve that specifies the outer boundary of the domain and any number of i...
SolidTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
virtual ~SolidTriangleMesh()
Empty Destructor.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Vector< TriangleMeshClosedCurve * > internal_closed_curve_pt() const
Helper function for getting the internal closed boundaries.
bool Boundary_refinement
Do not allow refinement of nodes on the boundary.
TriangleMeshClosedCurve * outer_boundary_pt(const unsigned &i) const
Helper function for getting the i-th outer boundary.
bool is_automatic_creation_of_vertices_on_boundaries_allowed()
Returns the status of the variable Allow_automatic_creation_of_vertices_on_boundaries.
TriangleMeshParameters(Vector< TriangleMeshClosedCurve * > &outer_boundary_pt)
Constructor: Only takes the outer boundary, all the other parameters are stated with the specific par...
Vector< Vector< double > > extra_holes_coordinates() const
Helper function for getting the extra holes.
virtual ~TriangleMeshParameters()
Empty destructor.
std::map< unsigned, double > & target_area_for_region()
Helper function for getting access to the region's target areas.
void set_target_area_for_region(const unsigned &i, const double &area)
Helper function to specify target area for region.
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
bool is_use_attributes() const
Helper function for getting the status of use_attributes variable.
double & element_area()
Helper function for getting access to the element area.
Vector< TriangleMeshClosedCurve * > outer_boundary_pt() const
Helper function for getting the outer boundary.
void enable_use_attributes()
Helper function for enabling the use of attributes.
void disable_use_attributes()
Helper function for disabling the use of attributes.
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed)
double element_area() const
Helper function for getting the element area.
void disable_automatic_creation_of_vertices_on_boundaries()
Disables the creation of points (by Triangle) on the outer and internal boundaries.
bool Internal_boundary_refinement
Do not allow refinement of nodes on the internal boundary.
void enable_automatic_creation_of_vertices_on_boundaries()
Enables the creation of points (by Triangle) on the outer and internal boundaries.
void disable_internal_boundary_refinement()
Helper function for disabling the use of boundary refinement.
TriangleMeshParameters(TriangleMeshClosedCurve *outer_boundary_pt)
Constructor: Only takes the outer boundary, all the other parameters are stated with the specific par...
std::map< unsigned, Vector< double > > Regions_coordinates
Store the coordinates for defining extra regions The key on the map is the region id.
Vector< Vector< double > > & extra_holes_coordinates()
Helper function for getting access to the extra holes.
Vector< TriangleMeshClosedCurve * > & internal_closed_curve_pt()
Helper function for getting access to the internal closed boundaries.
void add_region_coordinates(const unsigned &i, Vector< double > ®ion_coordinates)
Helper function for getting the extra regions.
TriangleMeshClosedCurve *& outer_boundary_pt(const unsigned &i)
Helper function for getting access to the i-th outer boundary.
TriangleMeshParameters()
Constructor: Takes nothing and initializes the other parameters to the default ones.
Vector< TriangleMeshClosedCurve * > Internal_closed_curve_pt
Internal closed boundaries.
void set_communicator_pt(OomphCommunicator *comm_pt)
Function to set communicator (mesh is then assumed to be distributed)
Vector< TriangleMeshClosedCurve * > & outer_boundary_pt()
Helper function for getting access to the outer boundary.
Vector< TriangleMeshOpenCurve * > & internal_open_curves_pt()
Helper function for getting access to the internal open boundaries.
Vector< Vector< double > > Extra_holes_coordinates
Store the coordinates for defining extra holes.
Vector< TriangleMeshClosedCurve * > Outer_boundary_pt
The outer boundary.
bool Allow_automatic_creation_of_vertices_on_boundaries
Allows automatic creation of vertices along boundaries by Triangle.
Vector< TriangleMeshOpenCurve * > internal_open_curves_pt() const
Helper function for getting the internal open boundaries.
std::map< unsigned, double > Regions_areas
Target areas for regions; defaults to 0.0 which (luckily) implies "no specific target area" for trian...
bool Use_attributes
Define the use of attributes (regions)
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed Required to pass it to new distribut...
void enable_internal_boundary_refinement()
Helper function for enabling the use of boundary refinement.
std::map< unsigned, Vector< double > > & regions_coordinates()
Helper function for getting access to the regions coordinates.
void disable_boundary_refinement()
Helper function for disabling the use of boundary refinement.
bool is_boundary_refinement_allowed() const
Helper function for getting the status of boundary refinement.
void enable_boundary_refinement()
Helper function for enabling the use of boundary refinement.
bool is_internal_boundary_refinement_allowed() const
Helper function for getting the status of boundary refinement.
Vector< TriangleMeshOpenCurve * > Internal_open_curves_pt
Internal boundaries.
double Element_area
The element are when calling triangulate external routine.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
void flush_shared_boundary_element()
TimeStepper * Time_stepper_pt
Timestepper used to build elements.
const unsigned nshared_boundary_element(const unsigned &b)
void re_assign_initial_zeta_values_for_internal_boundary(const unsigned &b, Vector< std::list< FiniteElement * >> &old_segment_sorted_ele_pt, std::map< FiniteElement *, bool > &old_is_inverted)
Re-assign the boundary segments initial zeta (arclength) value for those internal boundaries that wer...
Vector< Vector< Vector< unsigned > > > Shared_boundaries_ids
Stores the boundaries ids created by the interaction of two processors Shared_boundaries_ids[iproc][j...
TriangleScaffoldMesh * Tmp_mesh_pt
Temporary scaffold mesh.
std::map< unsigned, Vector< unsigned > > & shared_boundary_from_processors()
Return the association of the shared boundaries with the processors.
Vector< TriangleMeshPolyLine * > & boundary_subpolylines(const unsigned &b)
Gets the vector of auxiliar polylines that will represent the given boundary (useful only when the bo...
bool is_node_on_shared_boundary(const unsigned &b, Node *const &node_pt)
Is the node on the shared boundary.
bool Use_attributes
Boolean flag to indicate whether to use attributes or not (required for multidomain meshes)
void compute_boundary_segments_connectivity_and_initial_zeta_values(const unsigned &b)
Compute the boundary segments connectivity for those boundaries that were splited during the distribu...
std::map< unsigned, unsigned > Shared_boundary_overlaps_internal_boundary
Stores information about those shared boundaries that lie over or over a segment of an internal bound...
void build_triangulateio(const std::string &poly_file_name, TriangulateIO &triangulate_io, bool &use_attributes)
Helper function to create TriangulateIO object (return in triangulate_io) from the ....
virtual ~TriangleMesh()
Destructor.
std::map< unsigned, unsigned > & shared_boundary_overlaps_internal_boundary()
Gets the storage that indicates if a shared boundary is part of an internal boundary.
const bool boundary_was_splitted(const unsigned &b)
Helper function to verify if a given boundary was splitted in the distribution process.
void flush_shared_boundary_node(const unsigned &b)
Flush the boundary nodes associated to the shared boundary b.
void create_shared_polylines_connections()
Establish the connections of the polylines previously marked as having connections....
std::map< unsigned, double > Regions_areas
Target areas for regions; defaults to 0.0 which (luckily) implies "no specific target area" for trian...
void update_triangulateio()
Update the triangulateio object to the current nodal positions.
const unsigned initial_shared_boundary_id()
The initial boundary id for shared boundaries.
unsigned Initial_shared_boundary_id
The initial boundary id for shared boundaries.
const unsigned shared_boundaries_ids(const unsigned &p, const unsigned &q, const unsigned &i) const
Node *& boundary_segment_node_pt(const unsigned &b, const unsigned &s, const unsigned &n)
Return pointer to node n on boundary b.
void break_loops_on_shared_polyline_load_balance_helper(const unsigned &initial_shd_bnd_id, std::list< Node * > &input_nodes, Vector< FiniteElement * > &input_boundary_element_pt, Vector< FiniteElement * > &input_boundary_face_element_pt, Vector< int > &input_face_index_element, const int &input_connect_to_the_left, const int &input_connect_to_the_right, Vector< std::list< Node * >> &output_sorted_nodes_pt, Vector< Vector< FiniteElement * >> &output_boundary_element_pt, Vector< Vector< FiniteElement * >> &output_boundary_face_element_pt, Vector< Vector< int >> &output_face_index_element, Vector< int > &output_connect_to_the_left, Vector< int > &output_connect_to_the_right)
Break any possible loop created by the sorted list of nodes that is used to create a new shared polyl...
Vector< Vector< unsigned > > shared_boundaries_ids(const unsigned &p) const
void shared_boundaries_in_this_processor(Vector< unsigned > &shared_boundaries_in_this_processor)
Get the shared boundaries ids living in the current processor.
void synchronize_boundary_coordinates(const unsigned &b)
In charge of sinchronize the boundary coordinates for internal boundaries that were split as part of ...
void update_triangulateio(Vector< Vector< double >> &internal_point)
Update the TriangulateIO object to the current nodal position and the centre hole coordinates.
void compute_holes_left_by_halo_elements_helper(Vector< Vector< double >> &output_holes_coordinates)
Compute the holes left by the halo elements, those adjacent to the shared boundaries.
void identify_boundary_segments_and_assign_initial_zeta_values(const unsigned &b, Vector< FiniteElement * > &input_face_ele_pt, const bool &is_internal_boundary, std::map< FiniteElement *, FiniteElement * > &face_to_bulk_element_pt)
Identify the segments from the old mesh (original mesh) in the new mesh (this) and assign initial and...
std::map< unsigned, Vector< int > > Face_index_at_shared_boundary
For the e-th finite element on shared boundary b, this is the index of the face that lies along that ...
std::map< unsigned, Vector< unsigned > > Shared_boundary_from_processors
Stores the processors involved in the generation of a shared boundary, in 2D two processors give rise...
const bool shared_boundary_overlaps_internal_boundary(const unsigned &shd_bnd_id)
Checks if the shared boundary overlaps an internal boundary.
const unsigned nshared_boundary_node(const unsigned &b)
std::map< unsigned, Vector< TriangleMeshPolyLine * > > Boundary_subpolylines
The polylines that will temporary represent the boundary that was splitted in the distribution proces...
const unsigned nshared_boundary_overlaps_internal_boundary()
Get the number of shared boundaries overlaping internal boundaries.
void flush_shared_boundary_element(const unsigned &b)
virtual void reset_boundary_element_info(Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
Virtual function to perform the reset boundary elements info routines. Generally used after load bala...
const unsigned nshared_boundary_polyline(const unsigned &p, const unsigned &c) const
void add_shared_boundary_node(const unsigned &b, Node *node_pt)
Add the node the shared boundary.
Vector< Vector< unsigned > > & shared_boundaries_ids(const unsigned &p)
const unsigned nshared_boundary_curves(const unsigned &p) const
const int check_connections_of_polyline_nodes(std::set< FiniteElement * > &element_in_processor_pt, const int &root_edge_bnd_id, std::map< std::pair< Node *, Node * >, bool > &overlapped_face, std::map< unsigned, std::map< Node *, bool >> &node_on_bnd_not_overlapped_by_shd_bnd, std::list< Node * > ¤t_polyline_nodes, std::map< unsigned, std::list< Node * >> &shared_bnd_id_to_sorted_list_node_pt, const unsigned &node_degree, Node *&new_node_pt, const bool called_from_load_balance=false)
Check for any possible connections that the array of sorted nodes have with any previous boundaries o...
void break_loops_on_shared_polyline_helper(const unsigned &initial_shd_bnd_id, std::list< Node * > &input_nodes, Vector< FiniteElement * > &input_boundary_element_pt, Vector< int > &input_face_index_element, const int &input_connect_to_the_left, const int &input_connect_to_the_right, Vector< std::list< Node * >> &output_sorted_nodes_pt, Vector< Vector< FiniteElement * >> &output_boundary_element_pt, Vector< Vector< int >> &output_face_index_element, Vector< int > &output_connect_to_the_left, Vector< int > &output_connect_to_the_right)
Break any possible loop created by the sorted list of nodes that is used to create a new shared polyl...
void create_tmp_polygons_helper(Vector< Vector< TriangleMeshPolyLine * >> &polylines_pt, Vector< TriangleMeshPolygon * > &polygons_pt)
Take the polylines from the shared boundaries and create temporary polygon representations of the dom...
std::map< unsigned, bool > Boundary_was_splitted
Flag to indicate if a polyline has been splitted during the distribution process, the boundary id of ...
const unsigned read_unsigned_line_helper(std::istream &read_file)
bool Triangulateio_exists
Boolean defining if Triangulateio object has been built or not.
void re_scale_re_assigned_initial_zeta_values_for_internal_boundary(const unsigned &b)
Re-scale the re-assigned zeta values for the boundary nodes, apply only for internal boundaries.
const unsigned nshared_boundaries(const unsigned &p, const unsigned &q) const
Access functions to boundaries shared with processors.
Vector< Vector< double > > Original_extra_holes_coordinates
Backup the original extra holes coordinates.
virtual void load_balance(const Vector< unsigned > &target_domain_for_local_non_halo_element)
Virtual function to perform the load balance routines.
std::map< unsigned, Vector< Node * > > Shared_boundary_node_pt
Stores the boundary nodes adjacent to the shared boundaries, these nodes are a subset of the halo and...
Vector< TriangleMeshPolyLine * > & shared_boundary_polyline_pt(const unsigned &p, const unsigned &c)
void create_tmp_open_curves_helper(Vector< Vector< TriangleMeshPolyLine * >> &sorted_open_curves_pt, Vector< TriangleMeshPolyLine * > &unsorted_shared_to_internal_poly_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt)
Take the polylines from the original open curves and created new temporaly representations of open cu...
void sort_polylines_helper(Vector< TriangleMeshPolyLine * > &unsorted_polylines_pt, Vector< Vector< TriangleMeshPolyLine * >> &sorted_polylines_pt)
Sorts the polylines so they be continuous and then we can create a closed or open curve from them.
void create_distributed_domain_representation(Vector< TriangleMeshPolygon * > &polygons_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt)
Creates the distributed domain representation. Joins the original boundaires, shared boundaries and c...
Vector< Vector< Node * > > & boundary_segment_node_pt(const unsigned &b)
Return direct access to nodes associated with a boundary but sorted in segments.
void add_shared_boundary_element(const unsigned &b, FiniteElement *ele_pt)
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
coord2_t cross(const Point &O, const Point &A, const Point &B)
2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product....
void operator=(const TriangleMesh &)=delete
Broken assignment operator.
Vector< Vector< Vector< TriangleMeshPolyLine * > > > Shared_boundary_polyline_pt
Stores the polyline representation of the shared boundaries Shared_boundary_polyline_pt[iproc][ncurve...
Vector< Vector< Vector< unsigned > > > shared_boundaries_ids() const
TriangleMeshPolyLine * shared_boundary_polyline_pt(const unsigned &p, const unsigned &c, const unsigned &i) const
Vector< unsigned > & shared_boundary_from_processors(const unsigned &b)
void create_shared_boundaries(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, const Vector< FiniteElement * > &backed_up_f_el_pt, std::map< Data *, std::set< unsigned >> &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status)
Creates the shared boundaries.
TriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
Node * shared_boundary_node_pt(const unsigned &b, const unsigned &n)
std::vector< Point > convex_hull(std::vector< Point > P)
Returns a list of points on the convex hull in counter-clockwise order. Note: the last point in the r...
const unsigned nshared_boundaries() const
Vector< unsigned > oomph_vertex_nodes_id()
Return the vector that contains the oomph-lib node number for all vertex nodes in the TriangulateIO r...
void remesh_from_internal_triangulateio()
Completely regenerate the mesh from the trianglateio structure.
Vector< unsigned > Oomph_vertex_nodes_id
Vector storing oomph-lib node number for all vertex nodes in the TriangulateIO representation of the ...
TriangleMesh(const std::string &poly_file_name, const double &element_area, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
Build mesh from poly file, with specified target area for all elements.
const unsigned final_shared_boundary_id()
The final boundary id for shared boundaries.
virtual void reestablish_distribution_info_for_restart(OomphCommunicator *comm_pt, std::istream &restart_file)
Virtual function used to re-establish any additional info. related with the distribution after a re-s...
const unsigned nboundary_subpolylines(const unsigned &b)
Gets the number of subpolylines that create the boundarya (useful only when the boundary is marked as...
const unsigned shared_boundary_overlapping_internal_boundary(const unsigned &shd_bnd_id)
Gets the boundary id of the internal boundary that the shared boundary lies on.
TriangleMesh(const TriangleMesh &dummy)=delete
Broken copy constructor.
bool First_time_compute_holes_left_by_halo_elements
Flag to know if it is the first time we are going to compute the holes left by the halo elements.
void create_shared_polyline(const unsigned &my_rank, const unsigned &shd_bnd_id, const unsigned &iproc, const unsigned &jproc, std::list< Node * > &sorted_nodes, const int &root_edge_bnd_id, Vector< FiniteElement * > &bulk_bnd_ele_pt, Vector< int > &face_index_ele, Vector< Vector< TriangleMeshPolyLine * >> &unsorted_polylines_pt, const int &connect_to_the_left_flag, const int &connect_to_the_right_flag)
Create the shared polyline and fill the data structured that keep all the information associated with...
void dump_distributed_info_for_restart(std::ostream &dump_file)
Used to dump info. related with distributed triangle meshes.
void get_halo_elements_on_all_procs(const unsigned &nproc, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, std::map< Data *, std::set< unsigned >> &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status, std::map< GeneralisedElement *, unsigned > &element_to_global_index, Vector< Vector< Vector< GeneralisedElement * >>> &output_halo_elements_pt)
Creates the halo elements on all processors Gets the halo elements on all processors,...
const bool boundary_marked_as_shared_boundary(const unsigned &b, const unsigned &isub)
Returns the value that indicates if a subpolyline of a given boundary continues been used as internal...
void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Overload set_mesh_level_time_stepper so that the stored time stepper now corresponds to the new times...
unsigned Final_shared_boundary_id
The final boundary id for shared boundaries.
Vector< Vector< Vector< unsigned > > > & shared_boundaries_ids()
void read_distributed_info_for_restart(std::istream &restart_file)
Used to read info. related with distributed triangle meshes.
void select_boundary_face_elements(Vector< FiniteElement * > &face_el_pt, const unsigned &b, bool &is_internal_boundary, std::map< FiniteElement *, FiniteElement * > &face_to_bulk_element_pt)
Select face element from boundary using the criteria to decide which of the two face elements should ...
void output_boundary_coordinates(const unsigned &b, std::ostream &outfile)
Output the nodes on the boundary and their respective boundary coordinates(into separate tecplot zone...
void generic_constructor(Vector< TriangleMeshPolygon * > &outer_boundary_pt, Vector< TriangleMeshPolygon * > &internal_polygon_pt, Vector< TriangleMeshOpenCurve * > &open_polylines_pt, const double &element_area, Vector< Vector< double >> &extra_holes_coordinates, std::map< unsigned, Vector< double >> ®ions_coordinates, std::map< unsigned, double > ®ions_areas, TimeStepper *time_stepper_pt, const bool &use_attributes, const bool &refine_boundary, const bool &refine_internal_boundary)
A general-purpose construction function that builds the mesh once the different specific constructors...
TriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
Constructor with the input files.
std::map< unsigned, std::vector< bool > > Boundary_marked_as_shared_boundary
Flag to indicate if an internal boundary will be used as shared boundary because there is overlapping...
void flush_face_index_at_shared_boundary()
void add_face_index_at_shared_boundary(const unsigned &b, const unsigned &i)
Vector< unsigned > & shared_boundaries_ids(const unsigned &p, const unsigned &q)
void flush_shared_boundary_polyline_pt()
int face_index_at_shared_boundary(const unsigned &b, const unsigned &e)
void flush_shared_boundary_node()
Flush ALL the shared boundary nodes.
void get_shared_boundaries_overlapping_internal_boundary(const unsigned &internal_bnd_id, Vector< unsigned > &shd_bnd_ids)
Gets the shared boundaries ids that overlap the given internal boundary.
bool triangulateio_exists()
Boolean defining if Triangulateio object has been built or not.
FiniteElement * shared_boundary_element_pt(const unsigned &b, const unsigned &e)
void get_element_edges_on_boundary(std::map< std::pair< Node *, Node * >, unsigned > &element_edges_on_boundary)
Get the element edges (pair of nodes, edges) that lie on a boundary (used to mark shared boundaries t...
TriangleMesh()
Empty constructor.
void update_holes_information_helper(Vector< TriangleMeshPolygon * > &polygons_pt, Vector< Vector< double >> &output_holes_coordinates)
Keeps those vertices that define a hole, those that are inside closed internal boundaries in the new ...
std::map< unsigned, Vector< FiniteElement * > > Shared_boundary_element_pt
Stores the boundary elements adjacent to the shared boundaries, these elements are a subset of the ha...
Vector< Node * > & boundary_segment_node_pt(const unsigned &b, const unsigned &s)
Return direct access to nodes associated with a segment of a given boundary.
void create_polylines_from_halo_elements_helper(const Vector< unsigned > &element_domain, std::map< GeneralisedElement *, unsigned > &element_to_global_index, std::set< FiniteElement * > &element_in_processor_pt, Vector< Vector< Vector< GeneralisedElement * >>> &input_halo_elements, std::map< std::pair< Node *, Node * >, unsigned > &elements_edges_on_boundary, Vector< Vector< Vector< TriangleMeshPolyLine * >>> &output_polylines_pt)
Creates polylines from the intersection of halo elements on all processors. The new polylines define ...
Vector< unsigned > shared_boundaries_ids(const unsigned &p, const unsigned &q) const
double H
Non-dimensional wall thickness. As in Jensen & Heil (2003) paper.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void triangulate(char *triswitches, struct oomph::TriangulateIO *in, struct oomph::TriangulateIO *out, struct oomph::TriangulateIO *vorout)
bool operator<(const Point &p) const