triangle_mesh.template.h
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26#ifndef OOMPH_TRIANGLE_MESH_HEADER
27#define OOMPH_TRIANGLE_MESH_HEADER
28// Config header generated by autoconfig
29#ifdef HAVE_CONFIG_H
30#include <oomph-lib-config.h>
31#endif
32
33#ifdef OOMPH_HAS_MPI
34
35// MPI headers
36#include <mpi.h>
37#endif
38
39#ifdef OOMPH_HAS_FPUCONTROLH
40#include <fpu_control.h>
41#endif
42
43
44// Standards
45#include <float.h>
46#include <iostream>
47#include <fstream>
48#include <string.h>
49#include <iomanip>
50
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"
56
57namespace oomph
58{
59#ifdef OOMPH_HAS_TRIANGLE_LIB
60
61 // Interface to triangulate function
62 //
63 // NOTE: POSTFIX ANY CALLS TO THIS FUNCTION BY
64 //--------------------------------------------
65 // #ifdef OOMPH_HAS_FPUCONTROLH
66 // // Reset flags that are tweaked by triangle; can cause nasty crashes
67 // fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
68 // _FPU_SETCW(cw);
69 // #endif
70 //
71 extern "C"
72 {
73 void triangulate(char* triswitches,
74 struct oomph::TriangulateIO* in,
75 struct oomph::TriangulateIO* out,
76 struct oomph::TriangulateIO* vorout);
77 }
78
79#endif
80
81
82 /// /////////////////////////////////////////////////////////////////////
83 /// /////////////////////////////////////////////////////////////////////
84 /// /////////////////////////////////////////////////////////////////////
85 /// /////////////////////////////////////////////////////////////////////
86 /// /////////////////////////////////////////////////////////////////////
87
88
89 //=========================================================================
90 /// Helper object for dealing with the parameters used for the
91 /// TriangleMesh objects
92 //=========================================================================
94 {
95 public:
96 /// Constructor: Only takes the outer boundary, all the other parameters
97 /// are stated with the specific parameters
100 Element_area(0.2),
101 Use_attributes(false),
105 Comm_pt(0)
106 {
107 }
108
109 /// Constructor: Only takes the outer boundary, all the other parameters
110 /// are stated with the specific parameters
112 : Element_area(0.2),
113 Use_attributes(false),
117 Comm_pt(0)
118 {
119 Outer_boundary_pt.resize(1);
121 }
122
123 /// Constructor: Takes nothing and initializes the other parameters to
124 /// the default ones
126 : Element_area(0.2),
127 Use_attributes(false),
131 Comm_pt(0)
132 {
133 }
134
135 /// Empty destructor
137
138 /// Helper function for getting the outer boundary
140 {
141 return Outer_boundary_pt;
142 }
143
144 /// Helper function for getting access to the outer boundary
146 {
147 return Outer_boundary_pt;
148 }
149
150 /// Helper function for getting the i-th outer boundary
152 {
153 return Outer_boundary_pt[i];
154 }
155
156 /// Helper function for getting access to the i-th outer boundary
158 {
159 return Outer_boundary_pt[i];
160 }
161
162 /// Helper function for getting the internal closed boundaries
164 {
166 }
167
168 /// Helper function for getting access to the internal
169 /// closed boundaries
171 {
173 }
174
175 /// Helper function for getting the internal open boundaries
177 {
179 }
180
181 /// Helper function for getting access to the internal
182 /// open boundaries
184 {
186 }
187
188 /// Helper function for getting the element area
189 double element_area() const
190 {
191 return Element_area;
192 }
193
194 /// Helper function for getting access to the element area
195 double& element_area()
196 {
197 return Element_area;
198 }
199
200 /// Helper function for getting the extra holes
202 {
204 }
205
206 /// Helper function for getting access to the extra holes
208 {
210 }
211
212 /// Helper function for getting the extra regions
213 void add_region_coordinates(const unsigned& i,
214 Vector<double>& region_coordinates)
215 {
216 // Verify if not using the default region number (zero)
217 if (i == 0)
218 {
219 std::ostringstream error_message;
220 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);
226 }
227
228 // First check if the region with the specified id does not already exist
229 std::map<unsigned, Vector<double>>::iterator it;
230 it = Regions_coordinates.find(i);
231
232 // If it is already a region defined with that id throw an error
233 if (it != Regions_coordinates.end())
234 {
235 std::ostringstream error_message;
236 error_message << "The region id (" << i << ") that you are using for"
237 << "defining\n"
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"
241 << std::endl;
242 OomphLibWarning(error_message.str(),
243 OOMPH_CURRENT_FUNCTION,
244 OOMPH_EXCEPTION_LOCATION);
245 }
246
247 // If it does not exist then create the map
248 Regions_coordinates[i] = region_coordinates;
249
250 // Automatically set the using of attributes to enable
252 }
253
254 /// Helper function for getting access to the regions coordinates
255 std::map<unsigned, Vector<double>>& regions_coordinates()
256 {
257 return Regions_coordinates;
258 }
259
260 /// Helper function to specify target area for region
261 void set_target_area_for_region(const unsigned& i, const double& area)
262 {
263 Regions_areas[i] = area;
264 }
265
266 /// Helper function for getting access to the region's target areas
267 std::map<unsigned, double>& target_area_for_region()
268 {
269 return Regions_areas;
270 }
271
272 /// Helper function for enabling the use of attributes
274 {
275 Use_attributes = true;
276 }
277
278 /// Helper function for disabling the use of attributes
280 {
281 Use_attributes = false;
282 }
283
284 /// Helper function for getting the status of use_attributes
285 /// variable
286 bool is_use_attributes() const
287 {
288 return Use_attributes;
289 }
290
291 /// Helper function for enabling the use of boundary refinement
293 {
294 Boundary_refinement = true;
295 }
296
297 /// Boolean to indicate if Mesh has been distributed
299 {
300 return (Comm_pt != 0);
301 }
302
303 /// Function to set communicator (mesh is then assumed to be distributed)
305 {
306 Comm_pt = comm_pt;
307 }
308
309 /// Read-only access fct to communicator (Null if mesh is not distributed)
311 {
312 return Comm_pt;
313 }
314
315 /// Helper function for disabling the use of boundary refinement
317 {
318 Boundary_refinement = false;
319 }
320
321 /// Helper function for getting the status of boundary refinement
323 {
324 return Boundary_refinement;
325 }
326
327 /// Helper function for enabling the use of boundary refinement
329 {
331 }
332
333 /// Helper function for disabling the use of boundary refinement
335 {
337 }
338
339 /// Helper function for getting the status of boundary refinement
341 {
343 }
344
345 /// Enables the creation of points (by Triangle) on the outer and
346 /// internal boundaries
348 {
350 }
351
352 /// Disables the creation of points (by Triangle) on the outer and
353 /// internal boundaries
355 {
357 }
358
359 /// Returns the status of the variable
360 /// Allow_automatic_creation_of_vertices_on_boundaries
362 {
364 }
365
366 protected:
367 /// The outer boundary
369
370 /// Internal closed boundaries
372
373 /// Internal boundaries
375
376 /// The element are when calling triangulate external routine
378
379 /// Store the coordinates for defining extra holes
381
382 /// Store the coordinates for defining extra regions
383 /// The key on the map is the region id
384 std::map<unsigned, Vector<double>> Regions_coordinates;
385
386 /// Target areas for regions; defaults to 0.0 which (luckily)
387 /// implies "no specific target area" for triangle!
388 std::map<unsigned, double> Regions_areas;
389
390 /// Define the use of attributes (regions)
392
393 /// Do not allow refinement of nodes on the boundary
395
396 /// Do not allow refinement of nodes on the internal boundary
398
399 /// Allows automatic creation of vertices along boundaries by
400 /// Triangle
402
403 /// Pointer to communicator -- set to NULL if mesh is not distributed
404 /// Required to pass it to new distributed meshes created at the
405 /// adaptation stage
407 };
408
409
410 /// /////////////////////////////////////////////////////////////////////
411 /// /////////////////////////////////////////////////////////////////////
412 /// /////////////////////////////////////////////////////////////////////
413 /// /////////////////////////////////////////////////////////////////////
414 /// /////////////////////////////////////////////////////////////////////
415
416
417 //============start_of_triangle_class===================================
418 /// Triangle mesh build with the help of the scaffold mesh coming
419 /// from the triangle mesh generator Triangle.
420 /// http://www.cs.cmu.edu/~quake/triangle.html
421 //======================================================================
422 template<class ELEMENT>
423 class TriangleMesh : public virtual TriangleMeshBase
424 {
425 public:
426 /// Empty constructor
428 {
429#ifdef OOMPH_HAS_TRIANGLE_LIB
430 // Using this constructor no Triangulateio object is built
431 Triangulateio_exists = false;
432 // By default allow the automatic creation of vertices along the
433 // boundaries by Triangle
435#ifdef OOMPH_HAS_MPI
436 // Initialize the flag to indicate this is the first time to
437 // compute the holes left by the halo elements
439#endif // #ifdef OOMPH_HAS_MPI
440
441#endif
442
443 // Mesh can only be built with 2D Telements.
444 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
445 }
446
447 /// Constructor with the input files
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)
454 {
455 // Mesh can only be built with 2D Telements.
456 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
457
458 // Initialize the value for allowing creation of points on boundaries
460 allow_automatic_creation_of_vertices_on_boundaries;
461
462#ifdef OOMPH_HAS_MPI
463 // Initialize the flag to indicate this is the first time to
464 // compute the holes left by the halo elements
466#endif // #ifdef OOMPH_HAS_MPI
467
468 // Store Timestepper used to build elements
469 Time_stepper_pt = time_stepper_pt;
470
471 // Check if we should use attributes. This is set to true if the .poly
472 // file specifies regions
473 bool should_use_attributes = false;
474
475#ifdef OOMPH_HAS_TRIANGLE_LIB
476 // Using this constructor build the triangulatio
478 node_file_name,
479 element_file_name,
480 poly_file_name,
482 should_use_attributes);
483
484 // Record that the triangulateio object has been created
486#endif
487
488 // Store the attributes
489 Use_attributes = should_use_attributes;
490
491 // Build scaffold
493 node_file_name, element_file_name, poly_file_name);
494
495 // Convert mesh from scaffold to actual mesh
496 build_from_scaffold(time_stepper_pt, should_use_attributes);
497
498 // kill the scaffold
499 delete this->Tmp_mesh_pt;
500 this->Tmp_mesh_pt = 0;
501
502 // Setup boundary coordinates for boundaries
503 unsigned nb = nboundary();
504 for (unsigned b = 0; b < nb; b++)
505 {
506 this->template setup_boundary_coordinates<ELEMENT>(b);
507 }
508 }
509
510 protected:
511#ifdef OOMPH_HAS_TRIANGLE_LIB
512
513 public:
514 /// Build mesh, based on the specifications on
515 /// TriangleMeshParameters
516 TriangleMesh(TriangleMeshParameters& triangle_mesh_parameters,
517 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
518 {
519 // Store the region target areas
520 Regions_areas = triangle_mesh_parameters.target_area_for_region();
521
522 // Mesh can only be built with 2D Telements.
523 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
524
525 // Initialize the value for allowing creation of points on boundaries
527 triangle_mesh_parameters
529
530 // Store Timestepper used to build elements
531 Time_stepper_pt = time_stepper_pt;
532
533#ifdef OOMPH_HAS_MPI
534 // Initialize the flag to indicate this is the first time to
535 // compute the holes left by the halo elements
537#endif // #ifdef OOMPH_HAS_MPI
538
539 // ********************************************************************
540 // First part - Get polylines representations
541 // ********************************************************************
542
543 // Create the polyline representation of all the boundaries and
544 // then create the mesh by calling to "generic_constructor()"
545
546 // Initialise highest boundary id
547 unsigned max_boundary_id = 0;
548
549 // *****************************************************************
550 // Part 1.1 - Outer boundary
551 // *****************************************************************
552 // Get the representation of the outer boundaries from the
553 // TriangleMeshParameters object
554 Vector<TriangleMeshClosedCurve*> outer_boundary_pt =
555 triangle_mesh_parameters.outer_boundary_pt();
556
557#ifdef PARANOID
558 // Verify that the outer_boundary_object_pt has been set
559 if (outer_boundary_pt.size() == 0)
560 {
561 std::stringstream error_message;
562 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);
569 } // if (outer_boundary_pt!=0)
570#endif
571
572 // Find the number of outer closed curves
573 unsigned n_outer_boundaries = outer_boundary_pt.size();
574
575 // Create the storage for the polygons that define the outer
576 // boundaries
577 Vector<TriangleMeshPolygon*> outer_boundary_polygon_pt(
578 n_outer_boundaries);
579
580 // Loop over the number of outer boundaries
581 for (unsigned i = 0; i < n_outer_boundaries; ++i)
582 {
583 // Get the polygon representation and compute the max boundary_id on
584 // each outer polygon. Does nothing (i.e. just returns a pointer to
585 // the outer boundary that was input) if the outer boundary is
586 // already a polygon
587 outer_boundary_polygon_pt[i] =
588 closed_curve_to_polygon_helper(outer_boundary_pt[i], max_boundary_id);
589 }
590
591 // *****************************************************************
592 // Part 1.2 - Internal closed boundaries (possible holes)
593 // *****************************************************************
594 // Get the representation of the internal closed boundaries from the
595 // TriangleMeshParameters object
596 Vector<TriangleMeshClosedCurve*> internal_closed_curve_pt =
597 triangle_mesh_parameters.internal_closed_curve_pt();
598
599 // Find the number of internal closed curves
600 unsigned n_internal_closed_curves = internal_closed_curve_pt.size();
601
602 // Create the storage for the polygons that define the internal closed
603 // boundaries (again nothing happens (as above) if an internal closed
604 // curve is already a polygon)
605 Vector<TriangleMeshPolygon*> internal_polygon_pt(
606 n_internal_closed_curves);
607
608 // Loop over the number of internal closed curves
609 for (unsigned i = 0; i < n_internal_closed_curves; ++i)
610 {
611 // Get the polygon representation and compute the max boundary_id on
612 // each internal polygon
613 internal_polygon_pt[i] = closed_curve_to_polygon_helper(
614 internal_closed_curve_pt[i], max_boundary_id);
615 }
616
617 // *****************************************************************
618 // Part 1.3 - Internal open boundaries
619 // *****************************************************************
620 // Get the representation of open boundaries from the
621 // TriangleMeshParameteres object
622 Vector<TriangleMeshOpenCurve*> internal_open_curve_pt =
623 triangle_mesh_parameters.internal_open_curves_pt();
624
625 // Find the number of internal open curves
626 unsigned n_internal_open_curves = internal_open_curve_pt.size();
627
628 // Create the storage for the polylines that define the open boundaries
629 Vector<TriangleMeshOpenCurve*> internal_open_curve_poly_pt(
630 n_internal_open_curves);
631
632 // Loop over the number of internal open curves
633 for (unsigned i = 0; i < n_internal_open_curves; i++)
634 {
635 // Get the open polyline representation and compute the max boundary_id
636 // on each open polyline (again, nothing happens if there are curve
637 // sections on the current internal open curve)
638 internal_open_curve_poly_pt[i] = create_open_curve_with_polyline_helper(
639 internal_open_curve_pt[i], max_boundary_id);
640 }
641
642 // ********************************************************************
643 // Second part - Get associated geom objects and coordinate limits
644 // ********************************************************************
645
646 // ***************************************************************
647 // Part 2.1 Outer boundary
648 // ***************************************************************
649 for (unsigned i = 0; i < n_outer_boundaries; i++)
650 {
652 outer_boundary_pt[i]);
653 }
654
655 // ***************************************************************
656 // Part 2.2 - Internal closed boundaries (possible holes)
657 // ***************************************************************
658 for (unsigned i = 0; i < n_internal_closed_curves; i++)
659 {
661 internal_closed_curve_pt[i]);
662 }
663
664 // ********************************************************************
665 // Part 2.3 - Internal open boundaries
666 // ********************************************************************
667 for (unsigned i = 0; i < n_internal_open_curves; i++)
668 {
670 internal_open_curve_pt[i]);
671 }
672
673 // ********************************************************************
674 // Third part - Creates the TriangulateIO object by calling the
675 // "generic_constructor()" function
676 // ********************************************************************
677 // Get all the other parameters from the TriangleMeshParameters object
678 // The maximum element area
679 const double element_area = triangle_mesh_parameters.element_area();
680
681 // The holes coordinates
682 Vector<Vector<double>> extra_holes_coordinates =
683 triangle_mesh_parameters.extra_holes_coordinates();
684
685 // The regions coordinates
686 std::map<unsigned, Vector<double>> regions =
687 triangle_mesh_parameters.regions_coordinates();
688
689 // If we use regions then we use attributes
690 const bool use_attributes = triangle_mesh_parameters.is_use_attributes();
691
692 const bool refine_boundary =
693 triangle_mesh_parameters.is_boundary_refinement_allowed();
694
695 const bool refine_internal_boundary =
696 triangle_mesh_parameters.is_internal_boundary_refinement_allowed();
697
698 if (!refine_internal_boundary && refine_boundary)
699 {
700 std::ostringstream error_stream;
701 error_stream
702 << "You have specified that Triangle may refine the outer boundary, "
703 "but\n"
704 << "not internal boundaries. Triangle does not support this "
705 "combination.\n"
706 << "If you do not want Triangle to refine internal boundaries, it "
707 "can't\n"
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";
712
713 throw OomphLibError(error_stream.str().c_str(),
714 OOMPH_CURRENT_FUNCTION,
715 OOMPH_EXCEPTION_LOCATION);
716 }
717
719 outer_boundary_polygon_pt,
720 internal_polygon_pt,
721 internal_open_curve_poly_pt,
722 element_area,
723 extra_holes_coordinates,
724 regions,
725 triangle_mesh_parameters.target_area_for_region(),
726 time_stepper_pt,
727 use_attributes,
728 refine_boundary,
729 refine_internal_boundary);
730
731 // Setup boundary coordinates for boundaries
732 unsigned nb = nboundary();
733
734#ifdef OOMPH_HAS_MPI
735 // Before calling setup boundary coordinates check if the mesh is
736 // marked as distrbuted
737 if (triangle_mesh_parameters.is_mesh_distributed())
738 {
739 // Set the mesh as distributed by passing the communicator
740 this->set_communicator_pt(triangle_mesh_parameters.communicator_pt());
741 }
742#endif
743
744 for (unsigned b = 0; b < nb; b++)
745 {
746 this->template setup_boundary_coordinates<ELEMENT>(b);
747 }
748
749 // Snap it!
751 }
752
753 /// Build mesh from poly file, with specified target
754 /// area for all elements.
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)
760 {
761 // Mesh can only be built with 2D Telements.
762 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
763
764 // Initialize the value for allowing creation of points on boundaries
766 allow_automatic_creation_of_vertices_on_boundaries;
767
768#ifdef OOMPH_HAS_MPI
769 // Initialize the flag to indicate this is the first time to
770 // compute the holes left by the halo elements
772#endif // #ifdef OOMPH_HAS_MPI
773
774 // Disclaimer
775 std::string message =
776 "This constructor hasn't been tested since last cleanup.\n";
778 message, "TriangleMesh::TriangleMesh()", OOMPH_EXCEPTION_LOCATION);
779
780 // Store Timestepper used to build elements
781 Time_stepper_pt = time_stepper_pt;
782
783 // Create the data structures required to call the triangulate function
784 TriangulateIO triangle_in;
785
786 // Input string for triangle
787 std::stringstream input_string_stream;
788
789 // MH: Like everything else, this hasn't been tested!
790 // used to be input_string_stream<<"-pA -a" << element_area << "q30";
791 input_string_stream << "-pA -a -a" << element_area << "q30";
792
793 // Verify if creation of new points on boundaries is allowed
794 if (!this->is_creation_of_vertices_on_boundaries_allowed())
795 {
796 input_string_stream << " -YY";
797 }
798
799 // Convert to a *char required by the triangulate function
800 char triswitches[100];
801 sprintf(triswitches, "%s", input_string_stream.str().c_str());
802
803 // Create a boolean to decide whether or not to use attributes.
804 // The value of this will only be changed in build_triangulateio
805 // depending on whether or not the .poly file contains regions
806 bool use_attributes = false;
807
808 // Build the input triangulateio object from the .poly file
809 build_triangulateio(poly_file_name, triangle_in, use_attributes);
810
811 // Store the attributes flag
812 Use_attributes = use_attributes;
813
814 // Build the triangulateio out object
815 triangulate(triswitches, &triangle_in, &Triangulateio, 0);
816
817#ifdef OOMPH_HAS_FPUCONTROLH
818 // Reset flags that are tweaked by triangle; can cause nasty crashes
819 fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
820 _FPU_SETCW(cw);
821#endif
822
823 // Build scaffold
825
826 // Convert mesh from scaffold to actual mesh
827 build_from_scaffold(time_stepper_pt, use_attributes);
828
829 // Kill the scaffold
830 delete this->Tmp_mesh_pt;
831 this->Tmp_mesh_pt = 0;
832
833 // Cleanup but leave hole alone
834 bool clear_hole_data = false;
835 TriangleHelper::clear_triangulateio(triangle_in, clear_hole_data);
836
837 // Setup boundary coordinates for boundaries
838 unsigned nb = nboundary();
839 for (unsigned b = 0; b < nb; b++)
840 {
841 this->template setup_boundary_coordinates<ELEMENT>(b);
842 }
843 }
844
845#endif
846
847 /// Broken copy constructor
848 TriangleMesh(const TriangleMesh& dummy) = delete;
849
850 /// Broken assignment operator
851 void operator=(const TriangleMesh&) = delete;
852
853 /// Destructor
855 {
856#ifdef OOMPH_HAS_TRIANGLE_LIB
858 {
860 }
861
862 std::set<TriangleMeshCurveSection*>::iterator it_polyline;
863 for (it_polyline = Free_curve_section_pt.begin();
864 it_polyline != Free_curve_section_pt.end();
865 it_polyline++)
866 {
867 delete (*it_polyline);
868 }
869
870 std::set<TriangleMeshPolygon*>::iterator it_polygon;
871 for (it_polygon = Free_polygon_pt.begin();
872 it_polygon != Free_polygon_pt.end();
873 it_polygon++)
874 {
875 delete (*it_polygon);
876 }
877
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();
881 it_open_polyline++)
882 {
883 delete (*it_open_polyline);
884 }
885
886#endif
887 }
888
889 /// Overload set_mesh_level_time_stepper so that the stored
890 /// time stepper now corresponds to the new timestepper
891 void set_mesh_level_time_stepper(TimeStepper* const& time_stepper_pt,
892 const bool& preserve_existing_data)
893 {
894 this->Time_stepper_pt = time_stepper_pt;
895 }
896
897#ifdef OOMPH_HAS_MPI
898
899 /// Compute the boundary segments connectivity for those
900 /// boundaries that were splited during the distribution process
902 const unsigned& b);
903
904 /// Re-assign the boundary segments initial zeta (arclength)
905 /// value for those internal boundaries that were splited during the
906 /// distribution process. Those boundaries that have one face element
907 /// at each side of the boundary
909 const unsigned& b,
910 Vector<std::list<FiniteElement*>>& old_segment_sorted_ele_pt,
911 std::map<FiniteElement*, bool>& old_is_inverted);
912
913 /// Re-scale the re-assigned zeta values for the boundary
914 /// nodes, apply only for internal boundaries
916 const unsigned& b);
917
918 /// Identify the segments from the old mesh (original mesh)
919 /// in the new mesh (this) and assign initial and final boundary
920 /// coordinates for the segments that create the boundary. (This is
921 /// the version called from the original mesh to identify its own
922 /// segments)
924 const unsigned& b,
925 Vector<FiniteElement*>& input_face_ele_pt,
926 const bool& is_internal_boundary,
927 std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt);
928
929 /// Identify the segments from the old mesh (original mesh)
930 /// in the new mesh (this) and assign initial and final boundary
931 /// coordinates for the segments that create the boundary
933 const unsigned& b, TriangleMesh<ELEMENT>* original_mesh_pt);
934
935 /// In charge of sinchronize the boundary coordinates for
936 /// internal boundaries that were split as part of the distribution
937 /// process. Called after setup_boundary_coordinates() for the
938 /// original mesh only
939 void synchronize_boundary_coordinates(const unsigned& b);
940
941 /// Select face element from boundary using the criteria to
942 /// decide which of the two face elements should be used on internal
943 /// boundaries
945 Vector<FiniteElement*>& face_el_pt,
946 const unsigned& b,
947 bool& is_internal_boundary,
948 std::map<FiniteElement*, FiniteElement*>& face_to_bulk_element_pt);
949
950 /// Return direct access to nodes associated with a boundary but
951 /// sorted in segments
953 {
954 return Boundary_segment_node_pt[b];
955 }
956
957 /// Return direct access to nodes associated with a segment of
958 /// a given boundary
960 const unsigned& s)
961 {
962 return Boundary_segment_node_pt[b][s];
963 }
964
965 /// Return pointer to node n on boundary b
966 Node*& boundary_segment_node_pt(const unsigned& b,
967 const unsigned& s,
968 const unsigned& n)
969 {
970 return Boundary_segment_node_pt[b][s][n];
971 }
972
973#endif // OOMPH_HAS_MPI
974
975#ifdef OOMPH_HAS_TRIANGLE_LIB
976
977 /// Update the TriangulateIO object to the current nodal position
978 /// and the centre hole coordinates.
980 {
981 // Move the hole center
982 // Get number of holes
983 unsigned nhole = Triangulateio.numberofholes;
984 unsigned count_coord = 0;
985 for (unsigned ihole = 0; ihole < nhole; ihole++)
986 {
987 Triangulateio.holelist[count_coord] += internal_point[ihole][0];
988 Triangulateio.holelist[count_coord + 1] += internal_point[ihole][1];
989
990 // Increment counter
991 count_coord += 2;
992 }
993
994 // Do the update
996 }
997
998 /// Update the triangulateio object to the current nodal positions
1000 {
1001 // Get number of points
1003 double new_x = 0;
1004 double new_y = 0;
1005
1006 // Loop over the points
1007 for (unsigned inod = 0; inod < nnode; inod++)
1008 {
1009 // Get the node Id to be updated
1010 unsigned count = Oomph_vertex_nodes_id[inod];
1011
1012 // Update vertices using the vertex_node_id giving for the TriangulateIO
1013 // vertex enumeration the corresponding oomphlib mesh enumeration
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;
1019 }
1020 }
1021
1022#ifdef OOMPH_HAS_MPI
1023 /// Used to dump info. related with distributed triangle meshes
1024 void dump_distributed_info_for_restart(std::ostream& dump_file);
1025
1026 const unsigned read_unsigned_line_helper(std::istream& read_file)
1027 {
1028 std::string input_string;
1029
1030 // Read line up to termination sign
1031 getline(read_file, input_string, '#');
1032
1033 // Ignore rest of line
1034 read_file.ignore(200, '\n');
1035
1036 // Convert
1037 return std::atoi(input_string.c_str());
1038 }
1039
1040 /// Used to read info. related with distributed triangle meshes
1041 void read_distributed_info_for_restart(std::istream& restart_file);
1042
1043 /// Virtual function used to re-establish any additional info. related with
1044 /// the distribution after a re-starting for triangle meshes
1046 OomphCommunicator* comm_pt, std::istream& restart_file)
1047 {
1048 std::ostringstream error_stream;
1049 error_stream << "Empty default reestablish disributed info method "
1050 << "called.\n";
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);
1055 }
1056
1057#endif // #ifdef OOMPH_HAS_MPI
1058
1059 /// Completely regenerate the mesh from the trianglateio structure
1061 {
1062 // Remove all the boundary node information
1063 this->remove_boundary_nodes();
1064
1065 // Delete exisiting nodes
1066 unsigned n_node = this->nnode();
1067 for (unsigned n = n_node; n > 0; --n)
1068 {
1069 delete this->Node_pt[n - 1];
1070 this->Node_pt[n - 1] = 0;
1071 }
1072 // Delete exisiting elements
1073 unsigned n_element = this->nelement();
1074 for (unsigned e = n_element; e > 0; --e)
1075 {
1076 delete this->Element_pt[e - 1];
1077 this->Element_pt[e - 1] = 0;
1078 }
1079 // Flush the storage
1081
1082 // Delete all boundary element information
1083 // ALH: Kick up the object hierarchy?
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();
1092
1093#ifdef OOMPH_HAS_MPI
1094 // Delete Halo(ed) information in the old mesh
1095 if (this->is_mesh_distributed())
1096 {
1097 this->Halo_node_pt.clear();
1098 this->Root_halo_element_pt.clear();
1099
1100 this->Haloed_node_pt.clear();
1101 this->Root_haloed_element_pt.clear();
1102
1103 this->External_halo_node_pt.clear();
1104 this->External_halo_element_pt.clear();
1105
1106 this->External_haloed_node_pt.clear();
1107 this->External_haloed_element_pt.clear();
1108 }
1109#endif
1110
1111 unsigned nbound = nboundary();
1112 Boundary_coordinate_exists.resize(nbound, false);
1113
1114 // Now build the new scaffold
1116
1117 // Triangulation has been created -- remember to wipe it!
1118 Triangulateio_exists = true;
1119
1120 // Convert mesh from scaffold to actual mesh
1122
1123 // Kill the scaffold
1124 delete this->Tmp_mesh_pt;
1125 this->Tmp_mesh_pt = 0;
1126
1127#ifdef OOMPH_HAS_MPI
1128 if (!this->is_mesh_distributed())
1129 {
1130 nbound = this->nboundary(); // The original number of boundaries
1131 }
1132 else
1133 {
1134 nbound = this->initial_shared_boundary_id();
1135 // NOTE: The total number of boundaries is the number of
1136 // original bondaries plus the number of shared boundaries, but
1137 // here we only establish boundary coordinates for the original
1138 // boundaries. Once all the info. related with the distribution
1139 // has been established then the number of boundaries is reset
1140 // to the correct one (after reset the halo/haloed scheme)
1141 }
1142#else
1143 nbound = this->nboundary(); // The original number of boundaries
1144#endif
1145
1146 // Setup boundary coordinates for boundaries
1147 for (unsigned b = 0; b < nbound; b++)
1148 {
1149 this->template setup_boundary_coordinates<ELEMENT>(b);
1150 }
1151
1152 // Snap nodes only if the mesh is not distributed, if the mesh is
1153 // distributed it will be called after the re-establishment of the
1154 // halo/haloed scheme, and the proper identification of the segments
1155 // in the boundary
1156 if (!this->is_mesh_distributed())
1157 {
1158 // Deform the boundary onto any geometric objects
1160 }
1161 }
1162
1163 /// Boolean defining if Triangulateio object has been built or not
1165 {
1166 return Triangulateio_exists;
1167 }
1168
1169#endif
1170
1171 /// Return the vector that contains the oomph-lib node number
1172 /// for all vertex nodes in the TriangulateIO representation of the mesh
1174 {
1175 return Oomph_vertex_nodes_id;
1176 }
1177
1178 /// Timestepper used to build elements
1180
1181 /// Boolean flag to indicate whether to use attributes or not (required
1182 /// for multidomain meshes)
1184
1185 protected:
1186 /// Target areas for regions; defaults to 0.0 which (luckily)
1187 /// implies "no specific target area" for triangle!
1188 std::map<unsigned, double> Regions_areas;
1189
1190 /// Build mesh from scaffold
1191 void build_from_scaffold(TimeStepper* time_stepper_pt,
1192 const bool& use_attributes);
1193
1194#ifdef OOMPH_HAS_TRIANGLE_LIB
1195
1196 /// Helper function to create TriangulateIO object (return in
1197 /// triangulate_io) from the .poly file
1198 void build_triangulateio(const std::string& poly_file_name,
1199 TriangulateIO& triangulate_io,
1200 bool& use_attributes);
1201
1202 /// A general-purpose construction function that builds the
1203 /// mesh once the different specific constructors have assembled the
1204 /// appropriate information.
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)
1217 {
1218 // Mesh can only be built with 2D Telements.
1219 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(2);
1220
1221#ifdef PARANOID
1222 if (element_area < 10e-14)
1223 {
1224 std::ostringstream warning_message;
1225 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);
1232 }
1233#endif
1234
1235 // Store the attribute flag
1236 Use_attributes = use_attributes;
1237
1238 // Store Timestepper used to build elements
1239 Time_stepper_pt = time_stepper_pt;
1240
1241 // Store outer polygon
1242 Outer_boundary_pt = outer_boundary_pt;
1243
1244 // Store internal polygons by copy constructor
1245 Internal_polygon_pt = internal_polygon_pt;
1246
1247 // Store internal polylines by copy constructor
1248 Internal_open_curve_pt = open_polylines_pt;
1249
1250 // Store the extra holes coordinates
1251 Extra_holes_coordinates = extra_holes_coordinates;
1252
1253 // Store the extra regions coordinates
1254 Regions_coordinates = regions_coordinates;
1255
1256 // Create the data structures required to call the triangulate function
1257 TriangulateIO triangulate_io;
1258
1259 // Initialize TriangulateIO structure
1261
1262 // Convert TriangleMeshPolyLine and TriangleMeshClosedCurvePolyLine
1263 // to a triangulateio object
1265 outer_boundary_pt,
1266 internal_polygon_pt,
1267 open_polylines_pt,
1268 extra_holes_coordinates,
1269 regions_coordinates,
1270 regions_areas,
1271 triangulate_io);
1272
1273 // Initialize TriangulateIO structure
1275
1276 // Triangulation has been created -- remember to wipe it!
1277 Triangulateio_exists = true;
1278
1279 // Input string for triangle
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);
1283
1284 // MH: Used to be:
1285 // input_string_stream<<"-pA -a" << element_area << " -q30" << std::fixed;
1286 // The repeated -a allows the specification of areas for different
1287 // regions (if any)
1288 input_string_stream << "-pA -a -a" << element_area << " -q30"
1289 << std::fixed;
1290
1291 // Verify if creation of new points on boundaries is allowed
1293 {
1294 input_string_stream << " -YY";
1295 }
1296
1297 // Suppress insertion of additional points on outer boundary
1298 if (refine_boundary == false)
1299 {
1300 input_string_stream << "-Y";
1301 // Add the extra flag to suppress additional points on interior segments
1302 if (refine_internal_boundary == false)
1303 {
1304 input_string_stream << "Y";
1305 }
1306 }
1307
1308 // Convert the Input string in *char required by the triangulate function
1309 char triswitches[100];
1310 sprintf(triswitches, "%s", input_string_stream.str().c_str());
1311
1312 // Build the mesh using triangulate function
1313 triangulate(triswitches, &triangulate_io, &Triangulateio, 0);
1314
1315#ifdef OOMPH_HAS_FPUCONTROLH
1316 // Reset flags that are tweaked by triangle; can cause nasty crashes
1317 fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
1318 _FPU_SETCW(cw);
1319#endif
1320
1321 // Build scaffold
1323
1324 // If we have filled holes then we must use the attributes
1325 if (!regions_coordinates.empty())
1326 {
1327 // Convert mesh from scaffold to actual mesh
1328 build_from_scaffold(time_stepper_pt, true);
1329 // Record the attribute flag
1330 Use_attributes = true;
1331 }
1332 // Otherwise use what was asked
1333 else
1334 {
1335 // Convert mesh from scaffold to actual mesh
1336 build_from_scaffold(time_stepper_pt, use_attributes);
1337 }
1338
1339 // Kill the scaffold
1340 delete this->Tmp_mesh_pt;
1341 this->Tmp_mesh_pt = 0;
1342
1343 // Cleanup but leave hole and regions alone since it's still used
1344 bool clear_hole_data = false;
1345 TriangleHelper::clear_triangulateio(triangulate_io, clear_hole_data);
1346 }
1347
1348 /// Boolean defining if Triangulateio object has been built or not
1350
1351#endif // OOMPH_HAS_TRIANGLE_LIB
1352
1353 /// Temporary scaffold mesh
1355
1356 /// Vector storing oomph-lib node number
1357 /// for all vertex nodes in the TriangulateIO representation of the mesh
1359
1360#ifdef OOMPH_HAS_MPI
1361
1362 public:
1363 /// The initial boundary id for shared boundaries
1365 {
1367 }
1368
1369 /// The final boundary id for shared boundaries
1371 {
1373 }
1374
1375
1376 protected:
1377 /// Get the shared boundaries ids living in the current processor
1380 {
1381#ifdef PARANOID
1382 // Used to check if there are repeated shared boundaries
1383 std::set<unsigned> shared_boundaries_in_this_processor_set;
1384#endif
1385 // Get the number of processors
1386 const unsigned n_proc = this->communicator_pt()->nproc();
1387 // Get the current processor
1388 const unsigned my_rank = this->communicator_pt()->my_rank();
1389 // Loop over all the processor and get the shared boundaries ids
1390 // associated with each processor
1391 for (unsigned iproc = 0; iproc < n_proc; iproc++)
1392 {
1393 // Work with other processors only
1394 if (iproc != my_rank)
1395 {
1396 // Get the number of boundaries shared with the "iproc"-th
1397 // processor
1398 const unsigned nshared_boundaries_with_iproc =
1399 this->nshared_boundaries(my_rank, iproc);
1400
1401 // If there are shared boundaries associated with the current
1402 // processor then add them
1403 if (nshared_boundaries_with_iproc > 0)
1404 {
1405 // Get the boundaries ids shared with "iproc"-th processor
1406 Vector<unsigned> bound_ids_shared_with_iproc;
1407 bound_ids_shared_with_iproc =
1408 this->shared_boundaries_ids(my_rank, iproc);
1409
1410 // Loop over shared boundaries with "iproc"-th processor
1411 for (unsigned bs = 0; bs < nshared_boundaries_with_iproc; bs++)
1412 {
1413 const unsigned bnd_id = bound_ids_shared_with_iproc[bs];
1414#ifdef PARANOID
1415 // Check that the current shared boundary id has not been
1416 // previously added
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())
1420 {
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 "
1425 "ids\n";
1426 throw OomphLibError(error.str(),
1427 OOMPH_CURRENT_FUNCTION,
1428 OOMPH_EXCEPTION_LOCATION);
1429 } // if (it != shared_boundaries_in_this_processor_set.end())
1430 shared_boundaries_in_this_processor_set.insert(bnd_id);
1431#endif
1432 shared_boundaries_in_this_processor.push_back(bnd_id);
1433 } // for (bs < nshared_boundaries_with_iproc)
1434
1435 } // if (nshared_boundaries_with_iproc > 0)
1436
1437 } // if (iproc != my_rank)
1438
1439 } // for (iproc < nproc)
1440 }
1441
1442 /// Access functions to boundaries shared with processors
1443 const unsigned nshared_boundaries(const unsigned& p,
1444 const unsigned& q) const
1445 {
1446 return Shared_boundaries_ids[p][q].size();
1447 }
1448
1450 {
1451 return Shared_boundaries_ids;
1452 }
1453
1455 {
1456 return Shared_boundaries_ids;
1457 }
1458
1460 {
1461 return Shared_boundaries_ids[p];
1462 }
1463
1465 {
1466 return Shared_boundaries_ids[p];
1467 }
1468
1470 const unsigned& q) const
1471 {
1472 return Shared_boundaries_ids[p][q];
1473 }
1474
1476 const unsigned& q)
1477 {
1478 return Shared_boundaries_ids[p][q];
1479 }
1480
1481 const unsigned shared_boundaries_ids(const unsigned& p,
1482 const unsigned& q,
1483 const unsigned& i) const
1484 {
1485 return Shared_boundaries_ids[p][q][i];
1486 }
1487
1488 const unsigned nshared_boundary_curves(const unsigned& p) const
1489 {
1490 return Shared_boundary_polyline_pt[p].size();
1491 }
1492
1493 const unsigned nshared_boundary_polyline(const unsigned& p,
1494 const unsigned& c) const
1495 {
1496 return Shared_boundary_polyline_pt[p][c].size();
1497 }
1498
1500 const unsigned& p, const unsigned& c)
1501 {
1502 return Shared_boundary_polyline_pt[p][c];
1503 }
1504
1506 const unsigned& c,
1507 const unsigned& i) const
1508 {
1509 return Shared_boundary_polyline_pt[p][c][i];
1510 }
1511
1512 const unsigned nshared_boundaries() const
1513 {
1514 return Shared_boundary_element_pt.size();
1515 }
1516
1517 const unsigned nshared_boundary_element(const unsigned& b)
1518 {
1519 // First check if the boundary exist
1520 std::map<unsigned, Vector<FiniteElement*>>::iterator it =
1522 if (it != Shared_boundary_element_pt.end())
1523 {
1524 return Shared_boundary_element_pt[b].size();
1525 }
1526 else
1527 {
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);
1533 }
1534 }
1535
1537 {
1539 }
1540
1541 void flush_shared_boundary_element(const unsigned& b)
1542 {
1543 // First check if the boundary exist
1544 std::map<unsigned, Vector<FiniteElement*>>::iterator it =
1546 if (it != Shared_boundary_element_pt.end())
1547 {
1548 Shared_boundary_element_pt[b].clear();
1549 }
1550 else
1551 {
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);
1557 }
1558 }
1559
1560 void add_shared_boundary_element(const unsigned& b, FiniteElement* ele_pt)
1561 {
1562 Shared_boundary_element_pt[b].push_back(ele_pt);
1563 }
1564
1566 const unsigned& e)
1567 {
1568 // First check if the boundary exist
1569 std::map<unsigned, Vector<FiniteElement*>>::iterator it =
1571 if (it != Shared_boundary_element_pt.end())
1572 {
1573 return Shared_boundary_element_pt[b][e];
1574 }
1575 else
1576 {
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);
1582 }
1583 }
1584
1586 {
1588 }
1589
1590 void add_face_index_at_shared_boundary(const unsigned& b, const unsigned& i)
1591 {
1592 Face_index_at_shared_boundary[b].push_back(i);
1593 }
1594
1595 int face_index_at_shared_boundary(const unsigned& b, const unsigned& e)
1596 {
1597 // First check if the boundary exist
1598 std::map<unsigned, Vector<int>>::iterator it =
1600 if (it != Face_index_at_shared_boundary.end())
1601 {
1603 }
1604 else
1605 {
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);
1611 }
1612 }
1613
1614 const unsigned nshared_boundary_node(const unsigned& b)
1615 {
1616 // First check if the boundary exist
1617 std::map<unsigned, Vector<Node*>>::iterator it =
1619 if (it != Shared_boundary_node_pt.end())
1620 {
1621 return Shared_boundary_node_pt[b].size();
1622 }
1623 else
1624 {
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);
1630 }
1631 }
1632
1633 /// Flush ALL the shared boundary nodes
1635 {
1637 }
1638
1639 /// Flush the boundary nodes associated to the shared boundary b
1640 void flush_shared_boundary_node(const unsigned& b)
1641 {
1642 Shared_boundary_node_pt[b].clear();
1643 }
1644
1645 /// Add the node the shared boundary
1646 void add_shared_boundary_node(const unsigned& b, Node* node_pt)
1647 {
1648 // Get the size of the Shared_boundary_node_pt vector
1649 const unsigned nbound_node = Shared_boundary_node_pt[b].size();
1650 bool node_already_on_this_boundary = false;
1651 // Loop over the vector
1652 for (unsigned n = 0; n < nbound_node; n++)
1653 {
1654 // is the current node here already?
1655 if (node_pt == Shared_boundary_node_pt[b][n])
1656 {
1657 node_already_on_this_boundary = true;
1658 }
1659 }
1660
1661 // Add the base node pointer to the vector if it's not there already
1662 if (!node_already_on_this_boundary)
1663 {
1664 Shared_boundary_node_pt[b].push_back(node_pt);
1665 }
1666 }
1667
1668 Node* shared_boundary_node_pt(const unsigned& b, const unsigned& n)
1669 {
1670 // First check if the boundary exist
1671 std::map<unsigned, Vector<Node*>>::iterator it =
1673 if (it != Shared_boundary_node_pt.end())
1674 {
1675 return Shared_boundary_node_pt[b][n];
1676 }
1677 else
1678 {
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);
1684 }
1685 }
1686
1687 /// Is the node on the shared boundary
1688 bool is_node_on_shared_boundary(const unsigned& b, Node* const& node_pt)
1689 {
1690 // First check if the boundary exist
1691 std::map<unsigned, Vector<Node*>>::iterator it =
1693 if (it != Shared_boundary_node_pt.end())
1694 {
1695 // Now check if the node lives on the shared boundary
1696 Vector<Node*>::iterator it_shd_nodes =
1697 std::find(Shared_boundary_node_pt[b].begin(),
1698 Shared_boundary_node_pt[b].end(),
1699 node_pt);
1700 // If the node is on this boundary
1701 if (it_shd_nodes != Shared_boundary_node_pt[b].end())
1702 {
1703 return true;
1704 }
1705 else // The node is not on the boundary
1706 {
1707 return false;
1708 }
1709 }
1710 else
1711 {
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);
1717 }
1718 }
1719
1720 /// Return the association of the shared boundaries with the processors
1721 std::map<unsigned, Vector<unsigned>>& shared_boundary_from_processors()
1722 {
1724 }
1725
1727 {
1728 std::map<unsigned, Vector<unsigned>>::iterator it =
1730#ifdef PARANOID
1731 if (it == Shared_boundary_from_processors.end())
1732 {
1733 std::ostringstream error_message;
1734 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);
1743 }
1744#endif
1745 return (*it).second;
1746 }
1747
1748 /// Get the number of shared boundaries overlaping internal
1749 /// boundaries
1751 {
1753 }
1754
1755 /// Checks if the shared boundary overlaps an internal boundary
1757 const unsigned& shd_bnd_id)
1758 {
1759 std::map<unsigned, unsigned>::iterator it =
1762 {
1763 return true;
1764 }
1765 return false;
1766 }
1767
1768 /// Gets the boundary id of the internal boundary that the
1769 /// shared boundary lies on
1771 const unsigned& shd_bnd_id)
1772 {
1773 std::map<unsigned, unsigned>::iterator it =
1775#ifdef PARANOID
1777 {
1778 std::ostringstream error_message;
1779 error_message << "The shared boundary (" << shd_bnd_id
1780 << ") does not lie on an internal "
1781 << "boundary!!!.\n"
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);
1788 }
1789#endif
1790 return (*it).second;
1791 }
1792
1793 /// Gets the shared boundaries ids that overlap the given
1794 /// internal boundary
1796 const unsigned& internal_bnd_id, Vector<unsigned>& shd_bnd_ids)
1797 {
1798 // Clear any data in the output storage
1799 shd_bnd_ids.clear();
1800 // Loop over the map and store in the output vector the shared
1801 // boundaries ids that overlap the internal boundary
1802 std::map<unsigned, unsigned>::iterator it =
1804 for (; it != Shared_boundary_overlaps_internal_boundary.end(); it++)
1805 {
1806 // If the second entry is the internal boundary, then add the
1807 // first entry to the output vector
1808 if ((*it).second == internal_bnd_id)
1809 {
1810 // Add the first entry
1811 shd_bnd_ids.push_back((*it).first);
1812 }
1813 } // loop over the map entries
1814
1815#ifdef PARANOID
1816 if (shd_bnd_ids.size() == 0)
1817 {
1818 std::ostringstream error_message;
1819 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);
1827 }
1828#endif
1829 }
1830
1831 /// Gets the storage that indicates if a shared boundary is part
1832 /// of an internal boundary
1833 std::map<unsigned, unsigned>& shared_boundary_overlaps_internal_boundary()
1834 {
1836 }
1837
1838 /// Helper function to verify if a given boundary was splitted
1839 /// in the distribution process
1840 const bool boundary_was_splitted(const unsigned& b)
1841 {
1842 std::map<unsigned, bool>::iterator it;
1843 it = Boundary_was_splitted.find(b);
1844 if (it == Boundary_was_splitted.end())
1845 {
1846 return false;
1847 }
1848 else
1849 {
1850 return (*it).second;
1851 }
1852 }
1853
1854 /// Gets the number of subpolylines that create the boundarya
1855 /// (useful only when the boundary is marked as split)
1856 const unsigned nboundary_subpolylines(const unsigned& b)
1857 {
1858 std::map<unsigned, Vector<TriangleMeshPolyLine*>>::iterator it;
1859 it = Boundary_subpolylines.find(b);
1860#ifdef PARANOID
1861 if (it == Boundary_subpolylines.end())
1862 {
1863 std::ostringstream error_message;
1864 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 "
1869 "to\n"
1870 << "be splitted or the polylines have been explicitly deleted "
1871 "before\n"
1872 << "being used.";
1873 throw OomphLibError(error_message.str(),
1874 OOMPH_CURRENT_FUNCTION,
1875 OOMPH_EXCEPTION_LOCATION);
1876 }
1877#endif
1878 return (*it).second.size();
1879 }
1880
1881 /// Gets the vector of auxiliar polylines that will represent
1882 /// the given boundary (useful only when the boundaries were
1883 /// split)
1885 {
1886 std::map<unsigned, Vector<TriangleMeshPolyLine*>>::iterator it;
1887 it = Boundary_subpolylines.find(b);
1888 if (it == Boundary_subpolylines.end())
1889 {
1890 std::ostringstream error_message;
1891 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 "
1896 "to\n"
1897 << "be splitted or the polylines have been explicitly deleted "
1898 "before\n"
1899 << "being used.";
1900 throw OomphLibError(error_message.str(),
1901 OOMPH_CURRENT_FUNCTION,
1902 OOMPH_EXCEPTION_LOCATION);
1903 }
1904 return (*it).second;
1905 }
1906
1907 /// Returns the value that indicates if a subpolyline of a
1908 /// given boundary continues been used as internal boundary or should
1909 /// be changed as shared boundary
1910 const bool boundary_marked_as_shared_boundary(const unsigned& b,
1911 const unsigned& isub)
1912 {
1913 std::map<unsigned, std::vector<bool>>::iterator it;
1915 if (it == Boundary_marked_as_shared_boundary.end())
1916 {
1917 // If no info. was found for the shared boundary then it may be
1918 // a non internal boundary, so no shared boundaries are
1919 // overlaping it
1920 return false;
1921 }
1922 return (*it).second[isub];
1923 }
1924
1925 /// The initial boundary id for shared boundaries
1927
1928 /// The final boundary id for shared boundaries
1930
1931 /// Stores the boundaries ids created by the interaction of two
1932 /// processors Shared_boundaries_ids[iproc][jproc] = Vector of shared
1933 /// boundaries ids "iproc" processor shares boundaries with "jproc"
1934 /// processor
1936
1937 /// Stores the processors involved in the generation of a shared
1938 /// boundary, in 2D two processors give rise to the creation of a
1939 /// shared boundary
1940 std::map<unsigned, Vector<unsigned>> Shared_boundary_from_processors;
1941
1942 /// Stores information about those shared boundaries that lie over or
1943 /// over a segment of an internal boundary (only used when using
1944 /// internal boundaries in the domain)
1946
1947 /// Stores the polyline representation of the shared boundaries
1948 /// Shared_boundary_polyline_pt[iproc][ncurve][npolyline] = polyline_pt
1950
1952 {
1954 }
1955
1956 /// Stores the boundary elements adjacent to the shared boundaries,
1957 /// these
1958 /// elements are a subset of the halo and haloed elements
1959 std::map<unsigned, Vector<FiniteElement*>> Shared_boundary_element_pt;
1960
1961 /// For the e-th finite element on shared boundary b, this is
1962 /// the index of the face that lies along that boundary
1963 std::map<unsigned, Vector<int>> Face_index_at_shared_boundary;
1964
1965 /// Stores the boundary nodes adjacent to the shared boundaries,
1966 /// these nodes are a subset of the halo and haloed nodes
1967 std::map<unsigned, Vector<Node*>> Shared_boundary_node_pt;
1968
1969 /// Flag to indicate if a polyline has been splitted during the
1970 /// distribution process, the boundary id of the polyline is used to
1971 /// indicate if spplited
1972 std::map<unsigned, bool> Boundary_was_splitted;
1973
1974 /// The polylines that will temporary represent the boundary that was
1975 /// splitted in the distribution process. Used to ease the sending of
1976 /// info. to Triangle during the adaptation process.
1977 std::map<unsigned, Vector<TriangleMeshPolyLine*>> Boundary_subpolylines;
1978
1979 /// Flag to indicate if an internal boundary will be used as shared
1980 /// boundary
1981 /// because there is overlapping of the internal boundary with the shared
1982 /// boundary
1983 std::map<unsigned, std::vector<bool>> Boundary_marked_as_shared_boundary;
1984
1985 /// Creates the distributed domain representation. Joins the
1986 /// original boundaires, shared boundaries and creates connections among
1987 /// them to create the new polygons that represent the distributed
1988 /// domain
1990 Vector<TriangleMeshPolygon*>& polygons_pt,
1991 Vector<TriangleMeshOpenCurve*>& open_curves_pt);
1992
1993 /// Sorts the polylines so they be continuous and then we can
1994 /// create a closed or open curve from them
1996 Vector<TriangleMeshPolyLine*>& unsorted_polylines_pt,
1997 Vector<Vector<TriangleMeshPolyLine*>>& sorted_polylines_pt);
1998
1999 /// Take the polylines from the shared boundaries and create
2000 /// temporary polygon representations of the domain
2003 Vector<TriangleMeshPolygon*>& polygons_pt);
2004
2005 /// Take the polylines from the original open curves and created
2006 /// new temporaly representations of open curves with the bits of
2007 /// original curves not overlapped by shared boundaries
2009 Vector<Vector<TriangleMeshPolyLine*>>& sorted_open_curves_pt,
2010 Vector<TriangleMeshPolyLine*>& unsorted_shared_to_internal_poly_pt,
2011 Vector<TriangleMeshOpenCurve*>& open_curves_pt);
2012
2013 /// Flag to know if it is the first time we are going to compute the
2014 /// holes left by the halo elements
2016
2017 /// Backup the original extra holes coordinates
2019
2020 /// Compute the holes left by the halo elements, those
2021 /// adjacent to the shared boundaries
2023 Vector<Vector<double>>& output_holes_coordinates);
2024
2025 /// Keeps those vertices that define a hole, those that are
2026 /// inside closed internal boundaries in the new polygons that define the
2027 /// domain. Delete those outside/inside the outer polygons (this is
2028 /// required since Triangle can not deal with vertices that define
2029 /// holes outside the new outer polygons of the domain)
2031 Vector<TriangleMeshPolygon*>& polygons_pt,
2032 Vector<Vector<double>>& output_holes_coordinates);
2033
2034 /// Check for any possible connections that the array of
2035 /// sorted nodes have with any previous boundaries or with
2036 /// itself. Return -1 if no connection was found, return -2 if the
2037 /// connection is with the same polyline, return the boundary id of
2038 /// the boundary to which the connection is performed
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,
2049 Node*& new_node_pt,
2050 const bool called_from_load_balance = false);
2051
2052 /// Establish the connections of the polylines previously marked
2053 /// as having connections. This connections were marked in the function
2054 /// TriangleMesh::create_polylines_from_halo_elements_helper().
2056
2057 /// Creates the shared boundaries
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);
2065
2066 /// Creates the halo elements on all processors
2067 /// Gets the halo elements on all processors, these elements are then used
2068 /// on the function that computes the shared boundaries among the processors
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);
2077
2078 /// Get the element edges (pair of nodes, edges) that lie
2079 /// on a boundary (used to mark shared boundaries that lie on
2080 /// internal boundaries)
2082 std::map<std::pair<Node*, Node*>, unsigned>& element_edges_on_boundary);
2083
2084 /// Creates polylines from the intersection of halo elements
2085 /// on all processors. The new polylines define the shared boundaries
2086 /// in the domain This get the polylines on ALL processors, that is
2087 /// why the three dimensions
2088 /// output_polylines_pt[iproc][ncurve][npolyline]
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);
2096
2097 /// Break any possible loop created by the sorted list of
2098 /// nodes that is used to create a new shared polyline
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);
2111
2112 /// Break any possible loop created by the sorted list of
2113 /// nodes that is used to create a new shared polyline (modified
2114 /// version for load balance)
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);
2129
2130 /// Create the shared polyline and fill the data structured
2131 /// that keep all the information associated with the creationg of the
2132 /// shared boundary
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);
2145
2146 public:
2147 /// Virtual function to perform the load balance routines
2148 virtual void load_balance(
2149 const Vector<unsigned>& target_domain_for_local_non_halo_element)
2150 {
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);
2157 }
2158
2159 /// Virtual function to perform the reset boundary elements info
2160 /// routines. Generally used after load balance.
2161 virtual void reset_boundary_element_info(
2162 Vector<unsigned>& ntmp_boundary_elements,
2163 Vector<Vector<unsigned>>& ntmp_boundary_elements_in_region,
2164 Vector<FiniteElement*>& deleted_elements);
2165
2166#endif // #ifdef OOMPH_HAS_MPI
2167
2168
2169 public:
2170 /// Output the nodes on the boundary and their respective boundary
2171 /// coordinates(into separate tecplot zones)
2172 void output_boundary_coordinates(const unsigned& b, std::ostream& outfile);
2173
2174
2175 private:
2176 // Reference :
2177 // http://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#C.2B.2B
2178
2179 // Monotone chain 2D convex hull algorithm.
2180 // Asymptotic complexity: O(n log n).
2181 // Practical performance: 0.5-1.0 seconds for n=1000000 on a 1GHz machine.
2182 typedef double coord_t; // coordinate type
2183 typedef double coord2_t; // must be big enough to hold 2*max(|coordinate|)^2
2184
2185 struct Point
2186 {
2188
2189 bool operator<(const Point& p) const
2190 {
2191 return x < p.x || (x == p.x && y < p.y);
2192 }
2193 };
2194
2195 /// 2D cross product of OA and OB vectors, i.e. z-component of their
2196 /// 3D cross product. Returns a positive value, if OAB makes a
2197 /// counter-clockwise turn, negative for clockwise turn, and zero if the
2198 /// points are collinear.
2199 coord2_t cross(const Point& O, const Point& A, const Point& B)
2200 {
2201 return (A.x - O.x) * (B.y - O.y) - (A.y - O.y) * (B.x - O.x);
2202 }
2203
2204 /// Returns a list of points on the convex hull in counter-clockwise
2205 /// order. Note: the last point in the returned list is the same as the
2206 /// first one.
2207 std::vector<Point> convex_hull(std::vector<Point> P)
2208 {
2209 int n = P.size(), k = 0;
2210 std::vector<Point> H(2 * n);
2211
2212 // Sort points lexicographically
2213 std::sort(P.begin(), P.end());
2214
2215 // Build lower hull
2216 for (int i = 0; i < n; ++i)
2217 {
2218 while (k >= 2 && cross(H[k - 2], H[k - 1], P[i]) <= 0) k--;
2219 H[k++] = P[i];
2220 }
2221
2222 // Build upper hull
2223 for (int i = n - 2, t = k + 1; i >= 0; i--)
2224 {
2225 while (k >= t && cross(H[k - 2], H[k - 1], P[i]) <= 0) k--;
2226 H[k++] = P[i];
2227 }
2228
2229 H.resize(k);
2230 return H;
2231 }
2232 };
2233
2234
2235 /// ///////////////////////////////////////////////////////////////////
2236 /// ///////////////////////////////////////////////////////////////////
2237 /// ///////////////////////////////////////////////////////////////////
2238
2239
2240 //=========================================================================
2241 /// Unstructured refineable Triangle Mesh
2242 //=========================================================================
2243 // Temporary flag to enable full annotation of RefineableTriangleMesh
2244 // comms
2245 //#define ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
2246 template<class ELEMENT>
2248 public virtual RefineableMeshBase
2249 {
2250 public:
2251 /// Function pointer to function that updates the
2252 /// mesh following the snapping of boundary nodes to the
2253 /// boundaries (e.g. to move boundary nodes very slightly
2254 /// to satisfy volume constraints)
2255 typedef void (*MeshUpdateFctPt)(Mesh* mesh_pt);
2256
2257 /// Function pointer to a function that can generate
2258 /// a point within the ihole-th hole, so that this can be
2259 /// overloaded by the user if they have a better way of
2260 /// doing it than our clunky default. The function
2261 /// should update the components of the
2262 /// Vector poly_pt->internal_point()
2263 typedef void (*InternalHolePointUpdateFctPt)(const unsigned& ihole,
2264 TriangleMeshPolygon* poly_pt);
2265
2266#ifdef OOMPH_HAS_TRIANGLE_LIB
2267
2268 /// Build mesh, based on the specifications on
2269 /// TriangleMeshParameters
2271 TriangleMeshParameters& triangle_mesh_parameters,
2272 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
2273 : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
2274 {
2275 // Initialise the data associated with adaptation
2276 initialise_adaptation_data();
2277
2278 // Initialise the data associated with boundary refinements
2279 initialise_boundary_refinement_data();
2280 }
2281
2282#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
2283
2284 /// Build mesh, based on the polyfiles
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)
2291 : TriangleMesh<ELEMENT>(
2292 node_file_name,
2293 element_file_name,
2294 poly_file_name,
2295 time_stepper_pt,
2296 allow_automatic_creation_of_vertices_on_boundaries)
2297 {
2298 // Create and fill the data structures to give rise to polylines so that
2299 // the mesh can use the adapt methods
2300 create_polylines_from_polyfiles(node_file_name, poly_file_name);
2301
2302 // Initialise the data associated with adaptation
2303 initialise_adaptation_data();
2304
2305 // Initialise the data associated with boundary refinements
2306 initialise_boundary_refinement_data();
2307 }
2308
2309 protected:
2310#ifdef OOMPH_HAS_TRIANGLE_LIB
2311
2312 /// Build mesh from specified triangulation and
2313 /// associated target areas for elements in it
2314 /// NOTE: This is used ONLY during adaptation and should not be used
2315 /// as a method of constructing a TriangleMesh object in demo drivers!
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)
2323 {
2324 // Initialise the data associated with adaptation
2325 initialise_adaptation_data();
2326
2327 // Initialise the data associated with boundary refinements
2328 initialise_boundary_refinement_data();
2329
2330 // Store Timestepper used to build elements
2331 this->Time_stepper_pt = time_stepper_pt;
2332
2333 // Create triangulateio object to refine
2334 TriangulateIO triangle_refine;
2335
2336 // Initialize triangulateio structure
2337 TriangleHelper::initialise_triangulateio(this->Triangulateio);
2338
2339 // Triangulation has been created -- remember to wipe it!
2340 this->Triangulateio_exists = true;
2341
2342 // Create refined TriangulateIO structure based on target areas
2343 this->refine_triangulateio(triangulate_io, target_area, triangle_refine);
2344
2345 // Input string for triangle
2346 std::stringstream input_string_stream;
2347 input_string_stream << "-pq30-ra";
2348
2349 // Verify if creation of new points on boundaries is allowed
2350 if (!allow_automatic_creation_of_vertices_on_boundaries)
2351 {
2352 input_string_stream << " -YY";
2353 }
2354
2355 // Copy the allowing of creation of points on the boundaries status
2356 this->Allow_automatic_creation_of_vertices_on_boundaries =
2357 allow_automatic_creation_of_vertices_on_boundaries;
2358
2359 // Store the attribute flag
2360 this->Use_attributes = use_attributes;
2361
2362 // Convert to a *char required by the triangulate function
2363 char triswitches[100];
2364 sprintf(triswitches, "%s", input_string_stream.str().c_str());
2365
2366 // Build triangulateio refined object
2367 triangulate(triswitches, &triangle_refine, &this->Triangulateio, 0);
2368
2369#ifdef OOMPH_HAS_FPUCONTROLH
2370 // Reset flags that are tweaked by triangle; can cause nasty crashes
2371 fpu_control_t cw = (_FPU_DEFAULT & ~_FPU_EXTENDED) | _FPU_DOUBLE;
2372 _FPU_SETCW(cw);
2373#endif
2374
2375 // Build scaffold
2376 this->Tmp_mesh_pt = new TriangleScaffoldMesh(this->Triangulateio);
2377
2378 // Convert mesh from scaffold to actual mesh
2379 this->build_from_scaffold(time_stepper_pt, use_attributes);
2380
2381 // Kill the scaffold
2382 delete this->Tmp_mesh_pt;
2383 this->Tmp_mesh_pt = 0;
2384
2385 // Cleanup but leave hole alone as it's still used
2386 bool clear_hole_data = false;
2387 TriangleHelper::clear_triangulateio(triangle_refine, clear_hole_data);
2388
2389#ifdef OOMPH_HAS_MPI
2390 // Before calling setup boundary coordinates check if the mesh
2391 // should be treated as a distributed mesh
2392 if (comm_pt != 0)
2393 {
2394 // Set the communicator which will mark the mesh as distributed
2395 this->set_communicator_pt(comm_pt);
2396 }
2397#endif
2398
2399 // Setup boundary coordinates for boundaries
2400 unsigned nb = nboundary();
2401 for (unsigned b = 0; b < nb; b++)
2402 {
2403 this->template setup_boundary_coordinates<ELEMENT>(b);
2404 }
2405 }
2406
2407 public:
2408#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
2409
2410 /// Empty Destructor
2412
2413 /// Enables info. and timings for tranferring of target
2414 /// areas
2416 {
2417 Print_timings_transfering_target_areas = true;
2418 }
2419
2420 /// Disables info. and timings for tranferring of target
2421 /// areas
2423 {
2424 Print_timings_transfering_target_areas = false;
2425 }
2426
2427 /// Enables the solution projection step during adaptation
2429 {
2430 Disable_projection = false;
2431 }
2432
2433 /// Disables the solution projection step during adaptation
2435 {
2436 Disable_projection = true;
2437 }
2438
2439 /// Enables info. and timings for projection
2441 {
2442 Print_timings_projection = true;
2443 }
2444
2445 /// Disables info. and timings for projection
2447 {
2448 Print_timings_projection = false;
2449 }
2450
2451 /// Read/write access to number of bins in the x-direction
2452 /// when transferring target areas by bin method. Only used if we
2453 /// don't have CGAL!
2455 {
2456 return Nbin_x_for_area_transfer;
2457 }
2458
2459 /// Read/write access to number of bins in the y-direction
2460 /// when transferring target areas by bin method. Only used if we
2461 /// don't have CGAL!
2463 {
2464 return Nbin_y_for_area_transfer;
2465 }
2466
2467 /// Read/write access to number of sample points from which
2468 /// we try to locate zeta by Newton method when transferring target areas
2469 /// using cgal-based sample point container. If Newton method doesn't
2470 /// converge from any of these we use the nearest sample point.
2472 {
2473 return Max_sample_points_for_limited_locate_zeta_during_target_area_transfer;
2474 }
2475
2476 /// Max element size allowed during adaptation
2478 {
2479 return Max_element_size;
2480 }
2481
2482 /// Min element size allowed during adaptation
2484 {
2485 return Min_element_size;
2486 }
2487
2488 /// Min angle before remesh gets triggered
2490 {
2491 return Min_permitted_angle;
2492 }
2493
2494 // Returns the status of using an iterative solver for the
2495 // projection problem
2497 {
2498 return Use_iterative_solver_for_projection;
2499 }
2500
2501 /// Enables the use of an iterative solver for the projection
2502 /// problem
2504 {
2505 Use_iterative_solver_for_projection = true;
2506 }
2507
2508 /// Enables the use of an iterative solver for the projection
2509 /// problem
2511 {
2512 Use_iterative_solver_for_projection = false;
2513 }
2514
2515 /// Enables printing of timings for adaptation
2516 void enable_print_timings_adaptation(const unsigned& print_level = 1)
2517 {
2518 set_print_level_timings_adaptation(print_level);
2519 }
2520
2521 /// Disables printing of timings for adaptation
2523 {
2524 Print_timings_level_adaptation = 0;
2525 }
2526
2527 /// Sets the printing level of timings for adaptation
2528 void set_print_level_timings_adaptation(const unsigned& print_level)
2529 {
2530 const unsigned max_print_level = 3;
2531 // If printing level is greater than max. printing level
2532 if (print_level > max_print_level)
2533 {
2534 Print_timings_level_adaptation = max_print_level;
2535 }
2536 else
2537 {
2538 Print_timings_level_adaptation = print_level;
2539 }
2540 }
2541
2542 /// Enables printing of timings for load balance
2543 void enable_print_timings_load_balance(const unsigned& print_level = 1)
2544 {
2545 set_print_level_timings_load_balance(print_level);
2546 }
2547
2548 /// Disables printing of timings for load balance
2550 {
2551 Print_timings_level_load_balance = 0;
2552 }
2553
2554 /// Sets the printing level of timings for load balance
2555 void set_print_level_timings_load_balance(const unsigned& print_level)
2556 {
2557 const unsigned max_print_level = 3;
2558 // If printing level is greater than max. printing level
2559 if (print_level > max_print_level)
2560 {
2561 Print_timings_level_load_balance = max_print_level;
2562 }
2563 else
2564 {
2565 Print_timings_level_load_balance = print_level;
2566 }
2567 }
2568
2569 /// Doc the targets for mesh adaptation
2570 void doc_adaptivity_targets(std::ostream& outfile)
2571 {
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;
2583 }
2584
2585 /// Refine mesh uniformly and doc process
2587 {
2588 // Set the element error to something big
2589 unsigned nelem = nelement();
2590 Vector<double> elem_error(nelem, DBL_MAX);
2591
2592 // Limit the min element size to 1/3 of the current minimum
2593 double backup = Min_element_size;
2594
2595 // Get current max and min element size
2596 double orig_max_area, orig_min_area;
2597 this->max_and_min_element_size(orig_max_area, orig_min_area);
2598
2599 // Limit
2600 Min_element_size = orig_min_area / 3.0;
2601
2602 // Do it...
2603 adapt(elem_error);
2604
2605 // Reset
2606 Min_element_size = backup;
2607 }
2608
2609 /// Unrefine mesh uniformly: Return 0 for success,
2610 /// 1 for failure (if unrefinement has reached the coarsest permitted
2611 /// level)
2613 {
2614 throw OomphLibError("unrefine_uniformly() not implemented yet",
2615 OOMPH_CURRENT_FUNCTION,
2616 OOMPH_EXCEPTION_LOCATION);
2617 // dummy return
2618 return 0;
2619 }
2620
2621 /// Adapt mesh, based on elemental error provided
2622 void adapt(const Vector<double>& elem_error);
2623
2624 /// Access to function pointer to function that updates the
2625 /// mesh following the snapping of boundary nodes to the
2626 /// boundaries (e.g. to move boundary nodes very slightly
2627 /// to satisfy volume constraints)
2628 MeshUpdateFctPt& mesh_update_fct_pt()
2629 {
2630 return Mesh_update_fct_pt;
2631 }
2632
2633 /// Access to function pointer to can be
2634 /// used to generate the internal point for the ihole-th
2635 /// hole.
2636 InternalHolePointUpdateFctPt& internal_hole_point_update_fct_pt()
2637 {
2638 return Internal_hole_point_update_fct_pt;
2639 }
2640
2641
2642#ifdef OOMPH_HAS_MPI
2643 unsigned nsorted_shared_boundary_node(unsigned& b)
2644 {
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())
2648 {
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);
2654 }
2655 return (*it).second.size();
2656 }
2657
2659 {
2660 Sorted_shared_boundary_node_pt.clear();
2661 }
2662
2663 Node* sorted_shared_boundary_node_pt(unsigned& b, unsigned& i)
2664 {
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())
2668 {
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);
2674 }
2675 return (*it).second[i];
2676 }
2677
2678
2680 {
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())
2684 {
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);
2690 }
2691 return (*it).second;
2692 }
2693
2694#endif // #ifdef OOMPH_HAS_MPI
2695
2696 /// Helper function to create polylines and fill associate data
2697 // structures, used when creating from a mesh from polyfiles
2698 void create_polylines_from_polyfiles(const std::string& node_file_name,
2699 const std::string& poly_file_name);
2700
2701#ifdef OOMPH_HAS_MPI
2702 // Fill the boundary elements structures when dealing with
2703 // shared boundaries that overlap internal boundaries
2704 void fill_boundary_elements_and_nodes_for_internal_boundaries();
2705
2706 // Fill the boundary elements structures when dealing with
2707 // shared boundaries that overlap internal boundaries. Document the
2708 // number of elements on the shared boundaries that go to internal
2709 // boundaries
2710 void fill_boundary_elements_and_nodes_for_internal_boundaries(
2711 std::ofstream& outfile);
2712
2713 /// Used to re-establish any additional info. related with
2714 /// the distribution after a re-starting for triangle meshes
2716 std::istream& restart_file)
2717 {
2718 // Ensure that the mesh is distributed
2719 if (this->is_mesh_distributed())
2720 {
2721 // Fill the structures for the boundary elements and face indexes
2722 // of the boundary elements
2723 this->fill_boundary_elements_and_nodes_for_internal_boundaries();
2724
2725 // Re-establish the shared boundary elements and nodes scheme
2726 // before re-establish halo and haloed elements
2727 this->reset_shared_boundary_elements_and_nodes(comm_pt);
2728
2729 // Sort the nodes on the boundaries so that they have the same order
2730 // on all the boundaries
2731 this->sort_nodes_on_shared_boundaries();
2732
2733 // Re-establish the halo(ed) scheme
2734 this->reset_halo_haloed_scheme();
2735
2736 // Re-set the number of boundaries to the original one
2737 const unsigned noriginal_boundaries =
2738 this->initial_shared_boundary_id();
2739 this->set_nboundary(noriginal_boundaries);
2740
2741 // Go through the original boudaries and re-establish the
2742 // boundary coordinates
2743 for (unsigned b = 0; b < noriginal_boundaries; b++)
2744 {
2745 // Identify the segment boundaries and re-call the
2746 // setup_boundary_coordinates() method for the new boundaries
2747 // from restart
2748 this->identify_boundary_segments_and_assign_initial_zeta_values(b,
2749 this);
2750
2751 if (this->boundary_geom_object_pt(b) != 0)
2752 {
2753 // Re-set the boundary coordinates
2754 this->template setup_boundary_coordinates<ELEMENT>(b);
2755 }
2756 }
2757
2758 // Deform the boundary onto any geometric objects
2759 this->snap_nodes_onto_geometric_objects();
2760
2761 } // if (this->is_mesh_distributed())
2762 }
2763
2764#endif // #ifdef OOMPH_HAS_MPI
2765
2766 /// Method used to update the polylines representation after restart
2768 {
2769#ifdef OOMPH_HAS_MPI
2770 // If the mesh is distributed then also update the shared
2771 // boundaries
2772 unsigned my_rank = 0;
2773 if (this->is_mesh_distributed())
2774 {
2775 my_rank = this->communicator_pt()->my_rank();
2776 }
2777#endif
2778
2779 // Update the polyline representation after restart
2780
2781 // First update all internal boundaries
2782 const unsigned ninternal = this->Internal_polygon_pt.size();
2783 for (unsigned i_internal = 0; i_internal < ninternal; i_internal++)
2784 {
2785 this->update_polygon_after_restart(
2786 this->Internal_polygon_pt[i_internal]);
2787 }
2788
2789 // then update the external boundaries
2790 const unsigned nouter = this->Outer_boundary_pt.size();
2791 for (unsigned i_outer = 0; i_outer < nouter; i_outer++)
2792 {
2793 this->update_polygon_after_restart(this->Outer_boundary_pt[i_outer]);
2794 }
2795
2796#ifdef OOMPH_HAS_MPI
2797 // If the mesh is distributed then also update the shared
2798 // boundaries
2799 if (this->is_mesh_distributed())
2800 {
2801 const unsigned ncurves = this->nshared_boundary_curves(my_rank);
2802 for (unsigned nc = 0; nc < ncurves; nc++)
2803 {
2804 // Update the shared polyline
2805 this->update_shared_curve_after_restart(
2806 this->Shared_boundary_polyline_pt[my_rank][nc] // shared_curve
2807 );
2808 }
2809
2810 } // if (this->is_mesh_distributed())
2811#endif // #ifdef OOMPH_HAS_MPI
2812
2813 const unsigned n_open_polyline = this->Internal_open_curve_pt.size();
2814 for (unsigned i = 0; i < n_open_polyline; i++)
2815 {
2816 this->update_open_curve_after_restart(this->Internal_open_curve_pt[i]);
2817 }
2818 }
2819
2820#ifdef OOMPH_HAS_MPI
2821
2822 /// Performs the load balancing for unstructured meshes, the
2823 /// load balancing strategy is based on mesh migration
2824 void load_balance(
2825 const Vector<unsigned>& input_target_domain_for_local_non_halo_element);
2826
2827 /// Use the first and second group of elements to find the
2828 /// intersection between them to get the shared boundary
2829 /// elements from the first and second group
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);
2837
2838 /// Creates the new shared boundaries, this method is also in
2839 /// charge of computing the shared boundaries ids of each processor
2840 /// and send that info. to all the processors
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);
2845
2846 /// Computes the degree of the nodes on the shared boundaries, the
2847 /// degree of the node is computed from the global graph created by the
2848 /// shared boundaries of all processors
2849 void compute_shared_node_degree_helper(
2850 Vector<Vector<FiniteElement*>>& unsorted_face_ele_pt,
2851 std::map<Node*, unsigned>& global_node_degree);
2852
2853 /// Sort the nodes on the new shared boundaries (after load
2854 /// balancing), computes the alias of the nodes and creates the adjacency
2855 /// matrix that represent the graph created by the shared edges between each
2856 /// pair of processors
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);
2862
2863 /// Get the nodes on the shared boundary (b), these are stored
2864 /// in the segment they belong
2865 void get_shared_boundary_segment_nodes_helper(
2866 const unsigned& shd_bnd_id, Vector<Vector<Node*>>& tmp_segment_nodes);
2867
2868#endif // #ifdef OOMPH_HAS_MPI
2869
2870 /// Get the nodes on the boundary (b), these are stored in
2871 /// the segment they belong (also used by the load balance method
2872 /// to re-set the number of segments per boundary after load
2873 /// balance has taken place)
2874 void get_boundary_segment_nodes_helper(
2875 const unsigned& b, Vector<Vector<Node*>>& tmp_segment_nodes);
2876
2877 /// Enable/disable unrefinement/refinement methods for original
2878 /// boundaries
2880 {
2881 Do_boundary_unrefinement_constrained_by_target_areas = true;
2882 }
2883
2885 {
2886 Do_boundary_unrefinement_constrained_by_target_areas = false;
2887 }
2888
2890 {
2891 Do_boundary_refinement_constrained_by_target_areas = true;
2892 }
2893
2895 {
2896 Do_boundary_refinement_constrained_by_target_areas = false;
2897 }
2898
2899 /// Enable/disable unrefinement/refinement methods for shared
2900 /// boundaries
2902 {
2903 Do_shared_boundary_unrefinement_constrained_by_target_areas = true;
2904 }
2905
2907 {
2908 Do_shared_boundary_unrefinement_constrained_by_target_areas = false;
2909 }
2910
2912 {
2913 Do_shared_boundary_refinement_constrained_by_target_areas = true;
2914 }
2915
2917 {
2918 Do_shared_boundary_refinement_constrained_by_target_areas = false;
2919 }
2920
2921
2922 protected:
2923 /// A map that stores the vertices that receive connections, they
2924 /// are identified by the boundary number that receive the connection
2925 /// This is necessary for not erasing them on the adaptation process,
2926 /// specifically for the un-refinement process
2927 std::map<unsigned, std::set<Vector<double>>> Boundary_connections_pt;
2928
2929 /// Verifies if the given boundary receives a connection, and
2930 /// if that is the case then returns the list of vertices that
2931 /// receive the connections
2932 const bool boundary_connections(const unsigned& b,
2933 const unsigned& c,
2934 std::set<Vector<double>>& vertices)
2935 {
2936 // Search for the given boundary
2937 std::map<unsigned, std::set<Vector<double>>>::iterator it =
2938 Boundary_connections_pt.find(b);
2939 // Was the boundary found?
2940 if (it != Boundary_connections_pt.end())
2941 {
2942 // Return the set of vertices that receive the connection
2943 vertices = (*it).second;
2944 return true;
2945 }
2946 else
2947 {
2948 return false;
2949 }
2950 }
2951
2952 /// Synchronise the vertices that are marked for non deletion
2953 // on the shared boundaries. Unrefinement of shared boundaries is
2954 // performed only if the candidate node is not marked for non deletion
2955 const void synchronize_shared_boundary_connections();
2956
2957 /// Mark the vertices that are not allowed for deletion by
2958 /// the unrefienment/refinement polyline methods. In charge of
2959 /// filling the Boundary_chunk_connections_pt structure
2960 void add_vertices_for_non_deletion();
2961
2962 /// Adds the vertices from the sources boundary that are
2963 /// repeated in the destination boundary to the list of non
2964 /// delete-able vertices in the destination boundary
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);
2970
2971 /// After unrefinement and refinement has taken place compute
2972 /// the new vertices numbers of the temporary representation of the
2973 // boundaries to connect.
2974 void create_temporary_boundary_connections(
2975 Vector<TriangleMeshPolygon*>& tmp_outer_polygons_pt,
2976 Vector<TriangleMeshOpenCurve*>& tmp_open_curves_pt);
2977
2978 /// After unrefinement and refinement has taken place compute
2979 /// the new vertices numbers of the boundaries to connect (in a
2980 /// distributed scheme it may be possible that the destination boundary
2981 /// does no longer exist, therefore the connection is suspended and
2982 /// resumed after the adaptation processor
2983 void restore_boundary_connections(
2984 Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
2985 Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
2986
2987 /// Restore the connections of the specific polyline
2988 /// The vertices numbering on the destination boundaries may have
2989 /// change because of (un)refinement in the destination boundaries.
2990 /// Also deals with connection that do not longer exist because the
2991 /// destination boundary does no longer exist because of the distribution
2992 /// process
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);
2997
2998 /// Resume the boundary connections that may have been
2999 /// suspended because the destination boundary is no part of the
3000 /// domain. The connections are no permanently suspended because if
3001 /// load balance takes place the destination boundary may be part of
3002 /// the new domain representation therefore the connection would
3003 /// exist
3004 void resume_boundary_connections(
3005 Vector<TriangleMeshPolyLine*>& resume_initial_connection_polyline_pt,
3006 Vector<TriangleMeshPolyLine*>& resume_final_connection_polyline_pt);
3007
3008 /// Computes the associated vertex number on the destination
3009 /// boundary
3010 bool get_connected_vertex_number_on_dst_boundary(
3011 Vector<double>& vertex_coordinates,
3012 const unsigned& dst_b_id,
3013 unsigned& vertex_number);
3014
3015 /// Helper function that performs the unrefinement process
3016 // on the specified boundary by using the provided vertices
3017 /// representation. Optional boolean is used to run it as test only (if
3018 /// true is specified as input) in which case vertex coordinates aren't
3019 /// actually modified. Returned boolean indicates if polyline was (or
3020 /// would have been -- if called with check_only=false) changed.
3021 bool unrefine_boundary(const unsigned& b,
3022 const unsigned& c,
3023 Vector<Vector<double>>& vector_bnd_vertices,
3024 double& unrefinement_tolerance,
3025 const bool& check_only = false);
3026
3027 /// Helper function that performs the refinement process
3028 /// on the specified boundary by using the provided vertices
3029 /// representation. Optional boolean is used to run it as test only (if
3030 /// true is specified as input) in which case vertex coordinates aren't
3031 /// actually modified. Returned boolean indicates if polyline was (or
3032 /// would have been -- if called with check_only=false) changed.
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);
3037
3038 // Helper function that applies the maximum length constraint
3039 // when it was specified. This will increase the number of points in
3040 // the current curve section in case that any segment on it does not
3041 // fulfils the requirement
3042 bool apply_max_length_constraint(
3043 Mesh* face_mesh_pt,
3044 Vector<Vector<double>>& vector_bnd_vertices,
3045 double& max_length_constraint);
3046
3047 /// Helper function that performs the unrefinement process on
3048 /// the specified boundary by using the provided vertices
3049 /// representation and the associated target area. Used only when the
3050 /// 'allow_automatic_creation_of_vertices_on_boundaries' flag is set to
3051 /// true.
3052 bool unrefine_boundary_constrained_by_target_area(
3053 const unsigned& b,
3054 const unsigned& c,
3055 Vector<Vector<double>>& vector_bnd_vertices,
3056 double& unrefinement_tolerance,
3057 Vector<double>& area_constraint);
3058
3059 /// Helper function that performs the refinement process on
3060 /// the specified boundary by using the provided vertices
3061 /// representation and the associated elements target area. Used
3062 /// only when the 'allow_automatic_creation_of_vertices_on_boundaries'
3063 /// flag is set to true.
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);
3069
3070 /// Helper function that performs the unrefinement process
3071 /// on the specified boundary by using the provided vertices
3072 /// representation and the associated target area.
3073 /// NOTE: This is the version that applies unrefinement to shared
3074 /// boundaries
3075 bool unrefine_shared_boundary_constrained_by_target_area(
3076 const unsigned& b,
3077 const unsigned& c,
3078 Vector<Vector<double>>& vector_bnd_vertices,
3079 Vector<double>& area_constraint);
3080
3081 /// Helper function that performs the refinement process
3082 /// on the specified boundary by using the provided vertices
3083 /// representation and the associated elements target area.
3084 /// NOTE: This is the version that applies refinement to shared
3085 /// boundaries
3086 bool refine_shared_boundary_constrained_by_target_area(
3087 Vector<Vector<double>>& vector_bnd_vertices,
3088 Vector<double>& area_constraint);
3089
3090 /// Flag that enables or disables boundary unrefinement (true by default)
3092
3093 /// Flag that enables or disables boundary refinement (true by default)
3095
3096 /// Flag that enables or disables boundary unrefinement (true by default)
3098
3099 /// Flag that enables or disables boundary unrefinement (true by default)
3101
3102 /// Set all the flags to true (the default values)
3104 {
3105 // All boundaries refinement and unrefinement are allowed by
3106 // default
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;
3111 }
3112
3113#ifdef OOMPH_HAS_MPI
3114 /// Stores the nodes in the boundaries in the same order in all the
3115 /// processors
3116 /// Sorted_shared_boundary_node_pt[bnd_id][i-th node] = Node*
3117 /// It is a map since the boundary id may not start at zero
3118 std::map<unsigned, Vector<Node*>> Sorted_shared_boundary_node_pt;
3119
3120 /// Sort the nodes on shared boundaries so that the processors
3121 /// that share a boundary agree with the order of the nodes on the
3122 /// boundary
3123 void sort_nodes_on_shared_boundaries();
3124
3125 /// Re-establish the shared boundary elements after the
3126 /// adaptation process (the updating of shared nodes is optional and
3127 /// performed by default)
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);
3133
3134 /// In charge of. re-establish the halo(ed) scheme on all processors.
3135 /// Sends info. to create halo elements and nodes on the processors
3136 /// that need it. It uses and all to all communication strategy therefore
3137 /// must be called on all processors.
3138 void reset_halo_haloed_scheme();
3139
3140 /// Compute the names of the nodes on shared boundaries in
3141 /// this (my_rank) processor with other processors. Also compute the
3142 /// names of nodes on shared boundaries of other processors with
3143 /// other processors (useful when there is an element that requires
3144 /// to be sent to this (my_rank) processor because there is a shared
3145 /// node between this (my_rank) and other processors BUT there is
3146 /// not a shared boundary between this and the other processor
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);
3153
3154 /// Get the original boundaries to which is associated each
3155 /// shared node, and send the info. to the related processors. We
3156 /// need to do this so that at the reset of halo(ed) info. stage,
3157 /// the info. be already updated.
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);
3162
3163 /// In charge of creating additional halo(ed) elements on
3164 /// those processors that have no shared boundaries in common but have
3165 /// shared nodes
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);
3173
3174 // ====================================================================
3175 // Methods for load balancing
3176 // ====================================================================
3177
3178 //#define ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION_LOAD_BALANCE
3179
3180 // *********************************************************************
3181 // BEGIN: Methods to perform load balance
3182 // *********************************************************************
3183
3184 /// Check if necessary to add the element to the new domain or if it
3185 /// has been previously added
3187 Vector<FiniteElement*>& new_elements_on_domain, FiniteElement*& ele_pt)
3188 {
3189 // Get the number of elements currently added to the new domain
3190 const unsigned nnew_elements_on_domain = new_elements_on_domain.size();
3191
3192 // Flag to state if has been added or not
3193 bool already_on_new_domain = false;
3194 unsigned new_domain_ele_index = 0;
3195
3196 for (unsigned e = 0; e < nnew_elements_on_domain; e++)
3197 {
3198 if (ele_pt == new_elements_on_domain[e])
3199 {
3200 // It's already there, so...
3201 already_on_new_domain = true;
3202 // ...set the index of this element
3203 new_domain_ele_index = e;
3204 break;
3205 }
3206 }
3207
3208 // Has it been found?
3209 if (!already_on_new_domain)
3210 {
3211 // Not found, so add it:
3212 new_elements_on_domain.push_back(ele_pt);
3213 // Return the index where it's just been added
3214 return nnew_elements_on_domain;
3215 }
3216 else
3217 {
3218 // Return the index where it was found
3219 return new_domain_ele_index;
3220 }
3221 }
3222
3223 /// Helper function to get the required elemental information from
3224 /// the element to be sent. This info. involves the association of
3225 /// the element to a boundary or region, and if its part of the
3226 /// halo(ed) elements within a processor
3227 void get_required_elemental_information_load_balance_helper(
3228 unsigned& iproc,
3229 Vector<Vector<FiniteElement*>>& f_haloed_ele_pt,
3230 FiniteElement* ele_pt);
3231
3232 /// Check if necessary to add the node to the new domain or if it has
3233 /// been already added
3234 unsigned try_to_add_node_pt_load_balance(Vector<Node*>& new_nodes_on_domain,
3235 Node*& node_pt)
3236 {
3237 // Get the number of nodes currently added to the new domain
3238 const unsigned nnew_nodes_on_domain = new_nodes_on_domain.size();
3239
3240 // Flag to state if has been added or not
3241 bool already_on_new_domain = false;
3242 unsigned new_domain_node_index = 0;
3243
3244 for (unsigned n = 0; n < nnew_nodes_on_domain; n++)
3245 {
3246 if (node_pt == new_nodes_on_domain[n])
3247 {
3248 // It's already there, so...
3249 already_on_new_domain = true;
3250 // ...set the index of this element
3251 new_domain_node_index = n;
3252 break;
3253 }
3254 }
3255
3256 // Has it been found?
3257 if (!already_on_new_domain)
3258 {
3259 // Not found, so add it:
3260 new_nodes_on_domain.push_back(node_pt);
3261 // Return the index where it's just been added
3262 return nnew_nodes_on_domain;
3263 }
3264 else
3265 {
3266 // Return the index where it was found
3267 return new_domain_node_index;
3268 }
3269 }
3270
3271 /// Helper function to add haloed node
3272 void add_node_load_balance_helper(
3273 unsigned& iproc,
3274 Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
3275 Vector<Node*>& new_nodes_on_domain,
3276 Node* nod_pt);
3277
3278 /// Helper function to get the required nodal information
3279 /// from an haloed node so that a fully-functional node (and
3280 /// therefore element) can be created on the receiving process
3281 /// (this is the specific version for the load balance strategy,
3282 /// the difference with the original method is that it checks if
3283 /// the node is on a shared boundary no associated with the current
3284 /// processor --my_rank--, or in a haloed element from other
3285 /// processors
3286 void get_required_nodal_information_load_balance_helper(
3287 Vector<Vector<FiniteElement*>>& f_halo_ele_pt,
3288 unsigned& iproc,
3289 Node* nod_pt);
3290
3291 /// Helper function to create elements on the loop
3292 /// process based on the info received in
3293 /// send_and_received_elements_nodes_info
3294 void create_element_load_balance_helper(
3295 unsigned& iproc,
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);
3306
3307 /// Helper function to create elements on the loop
3308 /// process based on the info received in
3309 /// send_and_received_elements_nodes_info
3310 /// This function is in charge of verify if the element is associated
3311 /// to a boundary and associate to it if that is the case
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);
3317
3318 /// Helper function to add a new node from load balance
3319 void add_received_node_load_balance_helper(
3320 Node*& new_nod_pt,
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,
3327 unsigned& iproc,
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);
3333
3334 /// Helper function which constructs a new node (on an
3335 /// element) with the information sent from the load balance
3336 /// process
3337 void construct_new_node_load_balance_helper(
3338 Node*& new_nod_pt,
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,
3345 unsigned& iproc,
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);
3351
3352 // *********************************************************************
3353 // END: Methods to perform load balance
3354 // *********************************************************************
3355
3356 // *********************************************************************
3357 // Start communication variables
3358 // *********************************************************************
3359 /// Vector of flat-packed doubles to be communicated with
3360 /// other processors
3362
3363 /// Counter used when processing vector of flat-packed
3364 /// doubles
3366
3367 /// Vector of flat-packed unsigneds to be communicated with
3368 /// other processors
3370
3371 /// Counter used when processing vector of flat-packed
3372 /// unsigneds
3374
3375#ifdef ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
3376 /// Temporary vector of strings to enable full annotation of
3377 /// RefineableTriangleMesh comms
3379#endif
3380
3381 // *********************************************************************
3382 // End communication variables
3383 // *********************************************************************
3384
3385 // *********************************************************************
3386 // Start communication functions
3387 // *********************************************************************
3388
3389 /// Check if necessary to add the element as haloed or if it has been
3390 /// previously added to the haloed scheme
3391 unsigned try_to_add_root_haloed_element_pt(const unsigned& p,
3392 GeneralisedElement*& el_pt)
3393 {
3394 // Loop over current storage
3395 unsigned n_haloed = this->nroot_haloed_element(p);
3396
3397 // Is this already an haloed element?
3398 bool already_haloed_element = false;
3399 unsigned haloed_el_index = 0;
3400 for (unsigned eh = 0; eh < n_haloed; eh++)
3401 {
3402 if (el_pt == this->root_haloed_element_pt(p, eh))
3403 {
3404 // It's already there, so...
3405 already_haloed_element = true;
3406 // ...set the index of this element
3407 haloed_el_index = eh;
3408 break;
3409 }
3410 }
3411
3412 // Has it been found?
3413 if (!already_haloed_element)
3414 {
3415 // Not found, so add it:
3416 this->add_root_haloed_element_pt(p, el_pt);
3417 // Return the index where it's just been added
3418 return n_haloed;
3419 }
3420 else
3421 {
3422 // Return the index where it was found
3423 return haloed_el_index;
3424 }
3425 }
3426
3427 /// Check if necessary to add the node as haloed or if it has been
3428 /// previously added to the haloed scheme
3429 unsigned try_to_add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
3430 {
3431 // Loop over current storage
3432 unsigned n_haloed_nod = this->nhaloed_node(p);
3433
3434 // Is this already an haloed node?
3435 bool is_an_haloed_node = false;
3436 unsigned haloed_node_index = 0;
3437 for (unsigned k = 0; k < n_haloed_nod; k++)
3438 {
3439 if (nod_pt == this->haloed_node_pt(p, k))
3440 {
3441 is_an_haloed_node = true;
3442 haloed_node_index = k;
3443 break;
3444 }
3445 }
3446
3447 // Has it been found?
3448 if (!is_an_haloed_node)
3449 {
3450 // Not found, so add it
3451 this->add_haloed_node_pt(p, nod_pt);
3452 // Return the index where it's just been added
3453 return n_haloed_nod;
3454 }
3455 else
3456 {
3457 // Return the index where it was found
3458 return haloed_node_index;
3459 }
3460 }
3461
3462 /// Helper function to get the required elemental information from
3463 /// an haloed element. This info. involves the association of the element
3464 /// to a boundary or region.
3465 void get_required_elemental_information_helper(unsigned& iproc,
3466 FiniteElement* ele_pt);
3467
3468 /// Helper function to get the required nodal information
3469 /// from a haloed node so that a fully-functional halo node (and
3470 /// therefore element) can be created on the receiving process
3471 void get_required_nodal_information_helper(unsigned& iproc, Node* nod_pt);
3472
3473 /// Helper function to add haloed node
3474 void add_haloed_node_helper(unsigned& iproc, Node* nod_pt);
3475
3476 /// Helper function to send back halo and haloed information
3477 void send_and_receive_elements_nodes_info(int& send_proc, int& recv_proc);
3478
3479 /// Helper function to create (halo) elements on the loop
3480 /// process based on the info received in send_and_received_located_info
3481 void create_halo_element(
3482 unsigned& iproc,
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);
3489
3490 /// Helper function to create (halo) elements on the loop
3491 /// process based on the info received in send_and_received_located_info
3492 /// This function is in charge of verify if the element is associated to
3493 /// a boundary
3494 void add_halo_element_helper(unsigned& iproc, FiniteElement* ele_pt);
3495
3496 /// Helper function to add halo node
3497 void add_halo_node_helper(
3498 Node*& new_nod_pt,
3499 Vector<Node*>& new_nodes_on_domain,
3500 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3501 other_proc_shd_bnd_node_pt,
3502 unsigned& iproc,
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);
3508
3509 /// Helper function which constructs a new halo node
3510 /// (on an element) with the information sent from the haloed process
3511 void construct_new_halo_node_helper(
3512 Node*& new_nod_pt,
3513 Vector<Node*>& new_nodes_on_domain,
3514 Vector<Vector<Vector<std::map<unsigned, Node*>>>>&
3515 other_proc_shd_bnd_node_pt,
3516 unsigned& iproc,
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);
3522
3523 /// Helper function that assigns/updates the references to the node
3524 /// so that it can be found with any other reference. The return
3525 /// value indicates whether or not a node was found on the same
3526 /// reference
3527 void update_other_proc_shd_bnd_node_helper(
3528 Node*& new_nod_pt,
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);
3538
3539 // *********************************************************************
3540 // End Communication funtions
3541 // *********************************************************************
3542
3543#endif // #ifdef OOMPH_HAS_MPI
3544
3545 /// Helper function that updates the input polygon's PSLG
3546 /// by using the end-points of elements from FaceMesh(es) that are
3547 /// constructed for the boundaries associated with the segments of the
3548 /// polygon. Optional boolean is used to run it as test only (if
3549 /// true is specified as input) in which case polygon isn't actually
3550 /// modified. Returned boolean indicates if polygon was (or would have
3551 /// been -- if called with check_only=false) changed.
3552 bool update_polygon_using_face_mesh(TriangleMeshPolygon* polygon_pt,
3553 const bool& check_only = false);
3554
3555 /// Helper function that updates the input open curve by using
3556 /// end-points of elements from FaceMesh(es) that are constructed for the
3557 /// boundaries associated with the polylines. Optional boolean is used to
3558 /// run it as test only (if true is specified as input) in which case the
3559 /// polylines are not actually modified. Returned boolean indicates if
3560 /// polylines were (or would have been -- if called with check_only=false)
3561 /// changed.
3562 bool update_open_curve_using_face_mesh(
3563 TriangleMeshOpenCurve* open_polyline_pt, const bool& check_only = false);
3564
3565 /// Generate a new PSLG representation of the inner hole
3566 /// boundaries. Optional boolean is used to run it as test only (if
3567 /// true is specified as input) in which case PSLG isn't actually
3568 /// modified. Returned boolean indicates if PSLG was (or would have
3569 /// been -- if called with check_only=false) changed.
3570 virtual bool surface_remesh_for_inner_hole_boundaries(
3571 Vector<Vector<double>>& internal_point_coord,
3572 const bool& check_only = false);
3573
3574 /// Snap the boundary nodes onto any curvilinear boundaries
3575 void snap_nodes_onto_boundary(RefineableTriangleMesh<ELEMENT>*& new_mesh_pt,
3576 const unsigned& b);
3577
3578 /// Helper function
3579 /// Creates an unsorted face mesh representation from the specified
3580 /// boundary id. It means that the elements are not sorted along the
3581 /// boundary
3582 void create_unsorted_face_mesh_representation(const unsigned& boundary_id,
3583 Mesh* face_mesh_pt);
3584
3585 /// Helper function
3586 /// Creates a sorted face mesh representation of the specified PolyLine
3587 /// It means that the elements are sorted along the boundary
3588 /// It also returns a map that indicated the inverted elements
3589 void create_sorted_face_mesh_representation(
3590 const unsigned& boundary_id,
3591 Mesh* face_mesh_pt,
3592 std::map<FiniteElement*, bool>& is_inverted,
3593 bool& inverted_face_mesh);
3594
3595 /// Helper function to construct face mesh representation of all
3596 /// polylines, possibly with segments re-distributed between polylines
3597 /// to maintain an appxroximately even sub-division of the polygon
3598 void get_face_mesh_representation(TriangleMeshPolygon* polygon_pt,
3599 Vector<Mesh*>& face_mesh_pt);
3600
3601 /// Helper function to construct face mesh representation of
3602 /// open curves
3603 void get_face_mesh_representation(TriangleMeshOpenCurve* open_polyline_pt,
3604 Vector<Mesh*>& face_mesh_pt);
3605
3606 /// Updates the polylines representation after restart
3607 void update_polygon_after_restart(TriangleMeshPolygon*& polygon_pt);
3608
3609 /// Updates the open curve representation after restart
3610 void update_open_curve_after_restart(TriangleMeshOpenCurve*& open_curve_pt);
3611
3612 /// Updates the polylines using the elements area as constraint for
3613 /// the number of points along the boundaries
3614 bool update_polygon_using_elements_area(TriangleMeshPolygon*& polygon_pt,
3615 const Vector<double>& target_area);
3616
3617 /// Updates the open curve but using the elements area instead
3618 /// of the default refinement and unrefinement methods
3619 bool update_open_curve_using_elements_area(
3620 TriangleMeshOpenCurve*& open_curve_pt, const Vector<double>& target_area);
3621
3622#ifdef OOMPH_HAS_MPI
3623 /// Updates the polylines using the elements area as
3624 /// constraint for the number of points along the boundaries
3625 bool update_shared_curve_using_elements_area(
3626 Vector<TriangleMeshPolyLine*>& vector_polyline_pt,
3627 const Vector<double>& target_areas);
3628
3629 /// Updates the shared polylines representation after restart
3630 void update_shared_curve_after_restart(
3631 Vector<TriangleMeshPolyLine*>& vector_polyline_pt);
3632
3633#endif // #ifdef OOMPH_HAS_MPI
3634
3635 /// Helper function to initialise data associated with adaptation
3637 {
3638 // Number of bins in the x-direction
3639 // when transferring target areas by bin method. Only used if we
3640 // don't have CGAL!
3641 this->Nbin_x_for_area_transfer = 100;
3642
3643 // Number of bins in the y-direction
3644 // when transferring target areas by bin method. Only used if we
3645 // don't have CGAL!
3646 this->Nbin_y_for_area_transfer = 100;
3647
3648 // Initialise "what it says" -- used when transferring target areas
3649 // using cgal-based sample point container
3650 Max_sample_points_for_limited_locate_zeta_during_target_area_transfer = 5;
3651
3652 // Set max and min targets for adaptation
3653 this->Max_element_size = 1.0;
3654 this->Min_element_size = 0.001;
3655 this->Min_permitted_angle = 15.0;
3656
3657 // By default we want to do projection
3658 this->Disable_projection = false;
3659
3660 // Use by default an iterative solver for the projection problem
3661 this->Use_iterative_solver_for_projection = true;
3662
3663 // Set the defaul value for printing level adaptation (default 0)
3664 this->Print_timings_level_adaptation = 0;
3665
3666 // Set the defaul value for printing level load balance (default 0)
3667 this->Print_timings_level_load_balance = 0;
3668
3669 // By default we want no info. about timings for transferring of
3670 // target areas
3671 this->Print_timings_transfering_target_areas = false;
3672
3673 // By default we want no info. about timings for projection
3674 this->Print_timings_projection = false;
3675
3676 // Initialise function pointer to function that updates the
3677 // mesh following the snapping of boundary nodes to the
3678 // boundaries (e.g. to move boundary nodes very slightly
3679 // to satisfy volume constraints)
3680 Mesh_update_fct_pt = 0;
3681
3682 // Initialise function point for update of internal hole
3683 // point to null
3684 Internal_hole_point_update_fct_pt = 0;
3685 }
3686
3687#ifdef OOMPH_HAS_TRIANGLE_LIB
3688
3689 /// Build a new TriangulateIO object from previous TriangulateIO
3690 /// based on target area for each element
3691 void refine_triangulateio(TriangulateIO& triangulate_io,
3692 const Vector<double>& target_area,
3693 TriangulateIO& triangle_refine);
3694
3695#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3696
3697 /// Compute target area based on the element's error and the
3698 /// error target; return minimum angle (in degrees)
3699 double compute_area_target(const Vector<double>& elem_error,
3700 Vector<double>& target_area)
3701 {
3702 double min_angle = DBL_MAX;
3703 unsigned count_unrefined = 0;
3704 unsigned count_refined = 0;
3705 this->Nrefinement_overruled = 0;
3706
3707 // Record max. area constraint set by region
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();
3712 it++)
3713 {
3714 unsigned r = (*it).first;
3715 unsigned nel = this->nregion_element(r);
3716 for (unsigned e = 0; e < nel; e++)
3717 {
3718 max_area_from_region[this->region_element_pt(r, e)] = (*it).second;
3719 }
3720 }
3721
3722 unsigned nel = this->nelement();
3723 for (unsigned e = 0; e < nel; e++)
3724 {
3725 // Get element
3726 FiniteElement* el_pt = this->finite_element_pt(e);
3727
3728 // Area
3729 double area = el_pt->size();
3730
3731 // Min angle based on vertex coordinates
3732 // (vertices are enumerated first)
3733 double ax = el_pt->node_pt(0)->x(0);
3734 double ay = el_pt->node_pt(0)->x(1);
3735
3736 double bx = el_pt->node_pt(1)->x(0);
3737 double by = el_pt->node_pt(1)->x(1);
3738
3739 double cx = el_pt->node_pt(2)->x(0);
3740 double cy = el_pt->node_pt(2)->x(1);
3741
3742 // Min angle
3743 double angle0 =
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)))) *
3748 min_angle = std::min(min_angle, angle0);
3749
3750 double angle1 =
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)))) *
3755 min_angle = std::min(min_angle, angle1);
3756
3757 double angle2 = 180.0 - angle0 - angle1;
3758 min_angle = std::min(min_angle, angle2);
3759
3760 // Mimick refinement in tree-based procedure: Target areas
3761 // for elements that exceed permitted error is 1/3 of their
3762 // current area, corresponding to a uniform sub-division.
3763 double size_ratio = 3.0;
3764 if (elem_error[e] > max_permitted_error())
3765 {
3766 // Reduce area
3767 target_area[e] = std::max(area / size_ratio, Min_element_size);
3768
3769 //...but also make sure we're below the max element size
3770 target_area[e] = std::min(target_area[e], Max_element_size);
3771
3772 if (target_area[e] != Min_element_size)
3773 {
3774 count_refined++;
3775 }
3776 else
3777 {
3778 this->Nrefinement_overruled++;
3779 }
3780 }
3781 else if (elem_error[e] < min_permitted_error())
3782 {
3783 // Increase the area
3784 target_area[e] = std::min(size_ratio * area, Max_element_size);
3785
3786 //...but also make sure we're above the min element size
3787 target_area[e] = std::max(target_area[e], Min_element_size);
3788
3789 if (target_area[e] != Max_element_size)
3790 {
3791 count_unrefined++;
3792 }
3793 }
3794 else
3795 {
3796 // Leave it alone but enforce size limits
3797 double area_leave_alone = std::max(area, Min_element_size);
3798 target_area[e] = std::min(area_leave_alone, Max_element_size);
3799 }
3800
3801 // Enforce max areas from regions
3802 std::map<FiniteElement*, double>::iterator it =
3803 max_area_from_region.find(el_pt);
3804 if (it != max_area_from_region.end())
3805 {
3806 target_area[e] = std::min(target_area[e], (*it).second);
3807 }
3808 }
3809
3810
3811 // Tell everybody
3812 this->Nrefined = count_refined;
3813 this->Nunrefined = count_unrefined;
3814
3815 if (this->Nrefinement_overruled != 0)
3816 {
3818 << "\nNOTE: Refinement of " << this->Nrefinement_overruled
3819 << " elements was "
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";
3824 }
3825 return min_angle;
3826 }
3827
3828 /// Number of bins in the x-direction
3829 /// when transferring target areas by bin method. Only used if we
3830 /// don't have CGAL!
3832
3833 /// Number of bins in the y-direction
3834 /// when transferring target areas by bin method. Only used if we
3835 /// don't have CGAL!
3837
3838 /// Default value for max. number of sample points used for
3839 /// locate_zeta when transferring target areas using cgal-based sample point
3840 /// container
3841 unsigned
3843
3844 /// Max permitted element size
3846
3847 /// Min permitted element size
3849
3850 /// Min angle before remesh gets triggered
3852
3853 /// Enable/disable solution projection during adaptation
3855
3856 /// Flag to indicate whether to use or not an iterative solver (CG
3857 /// with diagonal preconditioned) for the projection problem
3859
3860 /// Enable/disable printing timings for transfering target areas
3862
3863 /// Enable/disable printing timings for projection
3865
3866 /// The printing level for adaptation
3868
3869 /// The printing level for load balance
3871
3872 /// Function pointer to function that updates the
3873 /// mesh following the snapping of boundary nodes to the
3874 /// boundaries (e.g. to move boundary nodes very slightly
3875 /// to satisfy volume constraints)
3876 MeshUpdateFctPt Mesh_update_fct_pt;
3877
3878 /// Function pointer to function that can be set
3879 /// to update the position of the central point in internal
3880 /// holes
3881 InternalHolePointUpdateFctPt Internal_hole_point_update_fct_pt;
3882 };
3883
3884
3885 /// //////////////////////////////////////////////////////////////////
3886 /// //////////////////////////////////////////////////////////////////
3887 /// //////////////////////////////////////////////////////////////////
3888
3889
3890 //=========================================================================
3891 /// Unstructured Triangle Mesh upgraded to solid mesh
3892 //=========================================================================
3893 template<class ELEMENT>
3894 class SolidTriangleMesh : public virtual TriangleMesh<ELEMENT>,
3895 public virtual SolidMesh
3896 {
3897 public:
3898#ifdef OOMPH_HAS_TRIANGLE_LIB
3899
3900 /// Build mesh, based on closed curve that specifies
3901 /// the outer boundary of the domain and any number of internal
3902 /// clsed curves. Specify target area for uniform element size.
3904 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
3905 : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
3906 {
3907 // Assign the Lagrangian coordinates
3909 }
3910
3911#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3912
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)
3919 : TriangleMesh<ELEMENT>(
3920 node_file_name,
3921 element_file_name,
3922 poly_file_name,
3923 time_stepper_pt,
3924 allow_automatic_creation_of_vertices_on_boundaries)
3925 {
3926 // Assign the Lagrangian coordinates
3928 }
3929
3930 /// Empty Destructor
3932 };
3933
3934
3935 /// /////////////////////////////////////////////////////////////////////
3936 /// /////////////////////////////////////////////////////////////////////
3937 /// /////////////////////////////////////////////////////////////////////
3938
3939#ifdef OOMPH_HAS_TRIANGLE_LIB
3940
3941 //=========================================================================
3942 /// Unstructured refineable Triangle Mesh upgraded to solid mesh
3943 //=========================================================================
3944 template<class ELEMENT>
3946 : public virtual RefineableTriangleMesh<ELEMENT>,
3947 public virtual SolidMesh
3948 {
3949 public:
3950 /// Build mesh, based on the specifications on
3951 /// TriangleMeshParameter
3953 TriangleMeshParameters& triangle_mesh_parameters,
3954 TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
3955 : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt),
3956 RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters,
3957 time_stepper_pt)
3958 {
3959 // Assign the Lagrangian coordinates
3961 }
3962
3963 /// Build mesh from specified triangulation and
3964 /// associated target areas for elements in it.
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)
3972 : RefineableTriangleMesh<ELEMENT>(
3973 target_area,
3974 triangulate_io,
3975 time_stepper_pt,
3976 use_attributes,
3977 allow_automatic_creation_of_vertices_on_boundaries,
3978 comm_pt)
3979 {
3980 // Assign the Lagrangian coordinates
3982 }
3983
3984 /// Empty Destructor
3986 };
3987
3988#endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3989
3990} // namespace oomph
3991
3992#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:86
Information for documentation of results: Directory and file number to enable output in the form RESL...
A general Finite Element class.
Definition: elements.h:1313
double size() const
Calculate the size of the element (length, area, volume,...) in Eulerian computational coordinates....
Definition: elements.cc:4290
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
A Generalised Element class.
Definition: elements.h:73
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
A general mesh class.
Definition: mesh.h:67
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors.
Definition: mesh.h:75
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
std::map< unsigned, Vector< GeneralisedElement * > > Root_haloed_element_pt
Map of vectors holding the pointers to the root haloed elements.
Definition: mesh.h:103
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,...
Definition: mesh.h:91
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them....
Definition: mesh.h:407
std::map< unsigned, Vector< GeneralisedElement * > > External_haloed_element_pt
Map of vectors holding the pointers to the external haloed elements.
Definition: mesh.h:129
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:190
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:95
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
std::map< unsigned, Vector< Node * > > External_haloed_node_pt
Map of vectors holding the pointers to the external haloed nodes.
Definition: mesh.h:138
std::map< unsigned, Vector< GeneralisedElement * > > Root_halo_element_pt
Map of vectors holding the pointers to the root halo elements.
Definition: mesh.h:100
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:204
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
std::map< unsigned, Vector< Node * > > Halo_node_pt
Map of vectors holding the pointers to the halo nodes.
Definition: mesh.h:106
void set_communicator_pt(OomphCommunicator *comm_pt)
Function to set communicator (mesh is assumed to be distributed if the communicator pointer is non-nu...
Definition: mesh.h:1615
std::map< unsigned, Vector< Node * > > Haloed_node_pt
Map of vectors holding the pointers to the haloed nodes.
Definition: mesh.h:109
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed, i.e. if we don't have mpi).
Definition: mesh.h:1600
std::map< unsigned, Vector< Node * > > External_halo_node_pt
Map of vectors holding the pointers to the external halo nodes.
Definition: mesh.h:135
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
std::map< unsigned, Vector< GeneralisedElement * > > External_halo_element_pt
External halo(ed) elements are created as and when they are needed to act as source elements for the ...
Definition: mesh.h:126
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:54
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Base class for refineable meshes. Provides standardised interfaces for the following standard mesh ad...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
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.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
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...
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 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.
double & max_element_size()
Max element size allowed during adaptation.
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...
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_timings_projection()
Disables info. and timings for projection.
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...
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.
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.
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.
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...
void disable_shared_boundary_unrefinement_constrained_by_target_areas()
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.
double & min_element_size()
Min element size allowed during adaptation.
bool Disable_projection
Enable/disable solution projection during adaptation.
double Min_permitted_angle
Min angle before remesh gets triggered.
RefineableTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
bool Do_boundary_unrefinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
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.
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...
Node * sorted_shared_boundary_node_pt(unsigned &b, unsigned &i)
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)
MeshUpdateFctPt & mesh_update_fct_pt()
Access to function pointer to function that updates the mesh following the snapping of boundary nodes...
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.
double & min_permitted_angle()
Min angle before remesh gets triggered.
Vector< Node * > sorted_shared_boundary_node_pt(unsigned &b)
bool Print_timings_transfering_target_areas
Enable/disable printing timings for transfering target areas.
unsigned Print_timings_level_load_balance
The printing level for load balance.
unsigned nsorted_shared_boundary_node(unsigned &b)
unsigned & nbin_x_for_area_transfer()
Read/write access to number of bins in the x-direction when transferring target areas by bin method....
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...
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.
General SolidMesh class.
Definition: mesh.h:2562
void set_lagrangian_nodal_coordinates()
Make the current configuration the undeformed one by setting the nodal Lagrangian coordinates to thei...
Definition: mesh.cc:9564
////////////////////////////////////////////////////////////////// //////////////////////////////////...
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.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
Base class for triangle meshes (meshes made of 2D triangle elements). Note: we choose to template Tri...
Definition: triangle_mesh.h:52
TriangulateIO Triangulateio
TriangulateIO representation of the mesh.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
TriangleMeshClosedCurve * outer_boundary_pt(const unsigned &i) const
Helper function for getting the i-th outer boundary.
bool Boundary_refinement
Do not allow refinement of nodes on the boundary.
bool is_automatic_creation_of_vertices_on_boundaries_allowed()
Returns the status of the variable Allow_automatic_creation_of_vertices_on_boundaries.
Vector< TriangleMeshClosedCurve * > outer_boundary_pt() const
Helper function for getting the outer boundary.
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()
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.
virtual ~TriangleMeshParameters()
Empty destructor.
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.
void enable_use_attributes()
Helper function for enabling the use of attributes.
void disable_use_attributes()
Helper function for disabling the use of attributes.
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, double > & target_area_for_region()
Helper function for getting access to the region's target areas.
std::map< unsigned, Vector< double > > Regions_coordinates
Store the coordinates for defining extra regions The key on the map is the region id.
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed)
void add_region_coordinates(const unsigned &i, Vector< double > &region_coordinates)
Helper function for getting the extra regions.
std::map< unsigned, Vector< double > > & regions_coordinates()
Helper function for getting access to the regions coordinates.
Vector< TriangleMeshClosedCurve * > internal_closed_curve_pt() const
Helper function for getting the internal closed boundaries.
TriangleMeshParameters()
Constructor: Takes nothing and initializes the other parameters to the default ones.
Vector< TriangleMeshOpenCurve * > & internal_open_curves_pt()
Helper function for getting access to the internal open boundaries.
Vector< Vector< double > > extra_holes_coordinates() const
Helper function for getting the extra holes.
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< Vector< double > > Extra_holes_coordinates
Store the coordinates for defining extra holes.
Vector< TriangleMeshClosedCurve * > & outer_boundary_pt()
Helper function for getting access to the outer boundary.
double & element_area()
Helper function for getting access to the element area.
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.
TriangleMeshClosedCurve *& outer_boundary_pt(const unsigned &i)
Helper function for getting access to the i-th outer boundary.
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.
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() const
Helper function for getting the internal open boundaries.
Vector< TriangleMeshOpenCurve * > Internal_open_curves_pt
Internal boundaries.
double Element_area
The element are when calling triangulate external routine.
Class defining a polyline for use in Triangle Mesh generation.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
TimeStepper * Time_stepper_pt
Timestepper used to build elements.
const unsigned nshared_boundary_element(const unsigned &b)
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 * > &current_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...
std::map< unsigned, Vector< unsigned > > & shared_boundary_from_processors()
Return the association of the shared boundaries with the processors.
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.
Vector< Vector< Vector< unsigned > > > & shared_boundaries_ids()
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 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 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.
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.
const bool boundary_was_splitted(const unsigned &b)
Helper function to verify if a given boundary was splitted in the distribution process.
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 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.
Vector< Vector< unsigned > > & shared_boundaries_ids(const unsigned &p)
unsigned Initial_shared_boundary_id
The initial boundary id for shared boundaries.
FiniteElement * shared_boundary_element_pt(const unsigned &b, const unsigned &e)
const unsigned shared_boundaries_ids(const unsigned &p, const unsigned &q, const unsigned &i) const
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...
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 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...
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...
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.
Vector< Vector< Node * > > & boundary_segment_node_pt(const unsigned &b)
Return direct access to nodes associated with a boundary but sorted in segments.
void flush_shared_boundary_element(const unsigned &b)
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.
const unsigned nshared_boundary_curves(const unsigned &p) const
Vector< unsigned > & shared_boundary_from_processors(const unsigned &b)
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.
Vector< Vector< Vector< unsigned > > > shared_boundaries_ids() const
std::map< unsigned, unsigned > & shared_boundary_overlaps_internal_boundary()
Gets the storage that indicates if a shared boundary is part of an internal boundary.
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.
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...
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.
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< Node * > > Shared_boundary_node_pt
Stores the boundary nodes adjacent to the shared boundaries, these nodes are a subset of the halo and...
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...
Node * shared_boundary_node_pt(const unsigned &b, const unsigned &n)
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.
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.
Vector< Vector< Vector< TriangleMeshPolyLine * > > > Shared_boundary_polyline_pt
Stores the polyline representation of the shared boundaries Shared_boundary_polyline_pt[iproc][ncurve...
Vector< Vector< unsigned > > shared_boundaries_ids(const unsigned &p) const
TriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
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...
const unsigned nshared_boundaries() const
void remesh_from_internal_triangulateio()
Completely regenerate the mesh from the trianglateio structure.
Node *& boundary_segment_node_pt(const unsigned &b, const unsigned &s, const unsigned &n)
Return pointer to node n on boundary b.
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...
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.
Vector< TriangleMeshPolyLine * > & shared_boundary_polyline_pt(const unsigned &p, const unsigned &c)
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...
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 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.
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 > > &regions_coordinates, std::map< unsigned, double > &regions_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 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 dump_distributed_info_for_restart(std::ostream &dump_file)
Used to dump info. related with distributed triangle meshes.
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...
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...
Vector< unsigned > shared_boundaries_ids(const unsigned &p, const unsigned &q) const
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.
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 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 ...
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...
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 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...
void add_face_index_at_shared_boundary(const unsigned &b, const unsigned &i)
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.
void update_triangulateio(Vector< Vector< double > > &internal_point)
Update the TriangulateIO object to the current nodal position and the centre hole coordinates.
Vector< TriangleMeshPolyLine * > & boundary_subpolylines(const unsigned &b)
Gets the vector of auxiliar polylines that will represent the given boundary (useful only when the bo...
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.
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.
Vector< unsigned > & shared_boundaries_ids(const unsigned &p, const unsigned &q)
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...
TriangleMeshPolyLine * shared_boundary_polyline_pt(const unsigned &p, const unsigned &c, const unsigned &i) const
Triangle Mesh that is based on input files generated by the triangle mesh generator Triangle.
Vector< TriangleMeshOpenCurve * > Internal_open_curve_pt
Vector of open polylines that define internal curves.
bool is_automatic_creation_of_vertices_on_boundaries_allowed()
Returns the status of the variable Allow_automatic_creation_of_vertices_on_boundaries.
Vector< TriangleMeshPolygon * > Outer_boundary_pt
Polygon that defines outer boundaries.
Vector< TriangleMeshPolygon * > Internal_polygon_pt
Vector of polygons that define internal polygons.
void snap_nodes_onto_geometric_objects()
Snap the boundary nodes onto any curvilinear boundaries defined by geometric objects.
TriangleMeshOpenCurve * create_open_curve_with_polyline_helper(TriangleMeshOpenCurve *open_curve_pt, unsigned &max_bnd_id_local)
Helper function that creates and returns an open curve with the polyline representation of its consti...
std::map< unsigned, Vector< Vector< Node * > > > Boundary_segment_node_pt
Used to store the nodes associated to a boundary and to an specific segment (this only applies in dis...
std::map< unsigned, Vector< double > > Regions_coordinates
Storage for extra coordinates for regions. The key on the map is the region id.
Vector< Vector< double > > Extra_holes_coordinates
Storage for extra coordinates for holes.
bool Allow_automatic_creation_of_vertices_on_boundaries
Flag to indicate whether the automatic creation of vertices along the boundaries by Triangle is allow...
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region.
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
TriangleMeshPolygon * closed_curve_to_polygon_helper(TriangleMeshClosedCurve *closed_curve_pt, unsigned &max_bnd_id_local)
Helper function that returns a polygon representation for the given closed curve, it also computes th...
std::set< TriangleMeshOpenCurve * > Free_open_curve_pt
A set that contains the open curves created by this object therefore it is necessary to free their as...
std::set< TriangleMeshCurveSection * > Free_curve_section_pt
A set that contains the curve sections created by this object therefore it is necessary to free their...
std::map< unsigned, Vector< std::pair< double, double > > > Polygonal_vertex_arclength_info
Storage for pairs of doubles representing: .first: the arclength along the polygonal representation o...
std::set< TriangleMeshPolygon * > Free_polygon_pt
A set that contains the polygons created by this object therefore it is necessary to free their assoc...
std::map< unsigned, Vector< FiniteElement * > > Region_element_pt
Vector of elements in each region differentiated by attribute (the key of the map is the attribute)
void build_triangulateio(Vector< TriangleMeshPolygon * > &outer_polygons_pt, Vector< TriangleMeshPolygon * > &internal_polygons_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt, Vector< Vector< double > > &extra_holes_coordinates, std::map< unsigned, Vector< double > > &regions_coordinates, std::map< unsigned, double > &regions_areas, TriangulateIO &triangulate_io)
Create TriangulateIO object from outer boundaries, internal boundaries, and open curves....
void set_geom_objects_and_coordinate_limits_for_open_curve(TriangleMeshOpenCurve *input_open_curve_pt)
Stores the geometric objects associated to the curve sections that compound the open curve....
void set_geom_objects_and_coordinate_limits_for_close_curve(TriangleMeshClosedCurve *input_closed_curve_pt)
Stores the geometric objects associated to the curve sections that compound the closed curve....
std::map< unsigned, TriangleMeshCurveSection * > Boundary_curve_section_pt
A map that stores the associated curve section of the specified boundary id.
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
const double Pi
50 digits from maple
void get_required_nodal_information_helper(int &iproc, Node *nod_pt, Mesh *const &mesh_pt, int &n_cont_inter_values, Vector< unsigned > &send_unsigneds, Vector< double > &send_doubles)
Helper function to get the required nodal information from an external haloed node so that a fully-fu...
void create_triangulateio_from_polyfiles(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TriangulateIO &triangle_io, bool &use_attributes)
Create a triangulateio data file from ele node and poly files. This is used if the mesh is generated ...
void initialise_triangulateio(TriangulateIO &triangle_io)
Initialise TriangulateIO structure.
void clear_triangulateio(TriangulateIO &triangulate_io, const bool &clear_hole_data)
Clear TriangulateIO structure.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void triangulate(char *triswitches, struct oomph::TriangulateIO *in, struct oomph::TriangulateIO *out, struct oomph::TriangulateIO *vorout)
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...
bool operator<(const Point &p) const
The Triangle data structure, modified from the triangle.h header supplied with triangle 1....
double * pointlist
Pointer to list of points x coordinate followed by y coordinate.