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-2023 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 
57 namespace 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
98  TriangleMeshParameters(Vector<TriangleMeshClosedCurve*>& outer_boundary_pt)
100  Element_area(0.2),
101  Use_attributes(false),
102  Boundary_refinement(true),
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
111  TriangleMeshParameters(TriangleMeshClosedCurve* outer_boundary_pt)
112  : Element_area(0.2),
113  Use_attributes(false),
114  Boundary_refinement(true),
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),
128  Boundary_refinement(true),
131  Comm_pt(0)
132  {
133  }
134 
135  /// Empty destructor
137 
138  /// Helper function for getting the outer boundary
139  Vector<TriangleMeshClosedCurve*> outer_boundary_pt() const
140  {
141  return Outer_boundary_pt;
142  }
143 
144  /// Helper function for getting access to the outer boundary
145  Vector<TriangleMeshClosedCurve*>& outer_boundary_pt()
146  {
147  return Outer_boundary_pt;
148  }
149 
150  /// Helper function for getting the i-th outer boundary
151  TriangleMeshClosedCurve* outer_boundary_pt(const unsigned& i) const
152  {
153  return Outer_boundary_pt[i];
154  }
155 
156  /// Helper function for getting access to the i-th outer boundary
157  TriangleMeshClosedCurve*& outer_boundary_pt(const unsigned& i)
158  {
159  return Outer_boundary_pt[i];
160  }
161 
162  /// Helper function for getting the internal closed boundaries
163  Vector<TriangleMeshClosedCurve*> internal_closed_curve_pt() const
164  {
166  }
167 
168  /// Helper function for getting access to the internal
169  /// closed boundaries
170  Vector<TriangleMeshClosedCurve*>& internal_closed_curve_pt()
171  {
173  }
174 
175  /// Helper function for getting the internal open boundaries
176  Vector<TriangleMeshOpenCurve*> internal_open_curves_pt() const
177  {
179  }
180 
181  /// Helper function for getting access to the internal
182  /// open boundaries
183  Vector<TriangleMeshOpenCurve*>& internal_open_curves_pt()
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
201  Vector<Vector<double>> extra_holes_coordinates() const
202  {
204  }
205 
206  /// Helper function for getting access to the extra holes
207  Vector<Vector<double>>& extra_holes_coordinates()
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
298  bool is_mesh_distributed() const
299  {
300  return (Comm_pt != 0);
301  }
302 
303  /// Function to set communicator (mesh is then assumed to be distributed)
304  void set_communicator_pt(OomphCommunicator* comm_pt)
305  {
306  Comm_pt = comm_pt;
307  }
308 
309  /// Read-only access fct to communicator (Null if mesh is not distributed)
310  OomphCommunicator* communicator_pt() const
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
368  Vector<TriangleMeshClosedCurve*> Outer_boundary_pt;
369 
370  /// Internal closed boundaries
371  Vector<TriangleMeshClosedCurve*> Internal_closed_curve_pt;
372 
373  /// Internal boundaries
374  Vector<TriangleMeshOpenCurve*> Internal_open_curves_pt;
375 
376  /// The element are when calling triangulate external routine
377  double Element_area;
378 
379  /// Store the coordinates for defining extra holes
380  Vector<Vector<double>> Extra_holes_coordinates;
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
406  OomphCommunicator* Comm_pt;
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
434  this->Allow_automatic_creation_of_vertices_on_boundaries = true;
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
459  this->Allow_automatic_creation_of_vertices_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
477  TriangleHelper::create_triangulateio_from_polyfiles(
478  node_file_name,
479  element_file_name,
480  poly_file_name,
481  Triangulateio,
482  should_use_attributes);
483 
484  // Record that the triangulateio object has been created
485  Triangulateio_exists = true;
486 #endif
487 
488  // Store the attributes
489  Use_attributes = should_use_attributes;
490 
491  // Build scaffold
492  this->Tmp_mesh_pt = new TriangleScaffoldMesh(
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
526  this->Allow_automatic_creation_of_vertices_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  {
651  set_geom_objects_and_coordinate_limits_for_close_curve(
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  {
660  set_geom_objects_and_coordinate_limits_for_close_curve(
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  {
669  set_geom_objects_and_coordinate_limits_for_open_curve(
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 
718  this->generic_constructor(
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!
750  this->snap_nodes_onto_geometric_objects();
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
765  this->Allow_automatic_creation_of_vertices_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";
777  OomphLibWarning(
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
824  this->Tmp_mesh_pt = new TriangleScaffoldMesh(Triangulateio);
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
854  virtual ~TriangleMesh()
855  {
856 #ifdef OOMPH_HAS_TRIANGLE_LIB
858  {
859  TriangleHelper::clear_triangulateio(Triangulateio);
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
952  Vector<Vector<Node*>>& boundary_segment_node_pt(const unsigned& b)
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
959  Vector<Node*>& boundary_segment_node_pt(const unsigned& b,
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.
979  void update_triangulateio(Vector<Vector<double>>& internal_point)
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
1002  unsigned nnode = Triangulateio.numberofpoints;
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
1080  this->flush_element_and_node_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
1115  this->Tmp_mesh_pt = new TriangleScaffoldMesh(this->Triangulateio);
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
1159  this->snap_nodes_onto_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
1173  Vector<unsigned> oomph_vertex_nodes_id()
1174  {
1175  return Oomph_vertex_nodes_id;
1176  }
1177 
1178  /// Timestepper used to build elements
1179  TimeStepper* Time_stepper_pt;
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
1260  TriangleHelper::initialise_triangulateio(triangulate_io);
1261 
1262  // Convert TriangleMeshPolyLine and TriangleMeshClosedCurvePolyLine
1263  // to a triangulateio object
1264  UnstructuredTwoDMeshGeometryBase::build_triangulateio(
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
1274  TriangleHelper::initialise_triangulateio(Triangulateio);
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
1292  if (!this->is_automatic_creation_of_vertices_on_boundaries_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
1322  this->Tmp_mesh_pt = new TriangleScaffoldMesh(Triangulateio);
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
1354  TriangleScaffoldMesh* Tmp_mesh_pt;
1355 
1356  /// Vector storing oomph-lib node number
1357  /// for all vertex nodes in the TriangulateIO representation of the mesh
1358  Vector<unsigned> Oomph_vertex_nodes_id;
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
1370  const unsigned final_shared_boundary_id()
1371  {
1372  return Final_shared_boundary_id;
1373  }
1374 
1375 
1376  protected:
1377  /// Get the shared boundaries ids living in the current processor
1379  Vector<unsigned>& shared_boundaries_in_this_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 
1449  Vector<Vector<Vector<unsigned>>> shared_boundaries_ids() const
1450  {
1451  return Shared_boundaries_ids;
1452  }
1453 
1454  Vector<Vector<Vector<unsigned>>>& shared_boundaries_ids()
1455  {
1456  return Shared_boundaries_ids;
1457  }
1458 
1459  Vector<Vector<unsigned>> shared_boundaries_ids(const unsigned& p) const
1460  {
1461  return Shared_boundaries_ids[p];
1462  }
1463 
1464  Vector<Vector<unsigned>>& shared_boundaries_ids(const unsigned& p)
1465  {
1466  return Shared_boundaries_ids[p];
1467  }
1468 
1469  Vector<unsigned> shared_boundaries_ids(const unsigned& p,
1470  const unsigned& q) const
1471  {
1472  return Shared_boundaries_ids[p][q];
1473  }
1474 
1475  Vector<unsigned>& shared_boundaries_ids(const unsigned& p,
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 
1499  Vector<TriangleMeshPolyLine*>& shared_boundary_polyline_pt(
1500  const unsigned& p, const unsigned& c)
1501  {
1502  return Shared_boundary_polyline_pt[p][c];
1503  }
1504 
1505  TriangleMeshPolyLine* shared_boundary_polyline_pt(const unsigned& p,
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 
1565  FiniteElement* shared_boundary_element_pt(const unsigned& b,
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  {
1602  return Face_index_at_shared_boundary[b][e];
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 =
1618  Shared_boundary_node_pt.find(b);
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  {
1636  Shared_boundary_node_pt.clear();
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 =
1672  Shared_boundary_node_pt.find(b);
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 =
1692  Shared_boundary_node_pt.find(b);
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 
1726  Vector<unsigned>& shared_boundary_from_processors(const unsigned& b)
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)
1884  Vector<TriangleMeshPolyLine*>& boundary_subpolylines(const unsigned& b)
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
1935  Vector<Vector<Vector<unsigned>>> Shared_boundaries_ids;
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)
1945  std::map<unsigned, unsigned> Shared_boundary_overlaps_internal_boundary;
1946 
1947  /// Stores the polyline representation of the shared boundaries
1948  /// Shared_boundary_polyline_pt[iproc][ncurve][npolyline] = polyline_pt
1949  Vector<Vector<Vector<TriangleMeshPolyLine*>>> Shared_boundary_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
1995  void sort_polylines_helper(
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
2002  Vector<Vector<TriangleMeshPolyLine*>>& polylines_pt,
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
2018  Vector<Vector<double>> 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
2586  void refine_uniformly(DocInfo& doc_info)
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 
2679  Vector<Node*> sorted_shared_boundary_node_pt(unsigned& b)
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
2715  void reestablish_distribution_info_for_restart(OomphCommunicator* comm_pt,
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
3361  Vector<double> Flat_packed_doubles;
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
3369  Vector<unsigned> Flat_packed_unsigneds;
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
3378  Vector<std::string> Flat_packed_unsigneds_string;
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)))) *
3747  180.0 / MathematicalConstants::Pi;
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)))) *
3754  180.0 / MathematicalConstants::Pi;
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  {
3817  oomph_info
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.
3903  SolidTriangleMesh(TriangleMeshParameters& triangle_mesh_parameters,
3904  TimeStepper* time_stepper_pt = &Mesh::Default_TimeStepper)
3905  : TriangleMesh<ELEMENT>(triangle_mesh_parameters, time_stepper_pt)
3906  {
3907  // Assign the Lagrangian coordinates
3908  set_lagrangian_nodal_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
3927  set_lagrangian_nodal_coordinates();
3928  }
3929 
3930  /// Empty Destructor
3931  virtual ~SolidTriangleMesh() {}
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
3960  set_lagrangian_nodal_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
3981  set_lagrangian_nodal_coordinates();
3982  }
3983 
3984  /// Empty Destructor
3986  };
3987 
3988 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3989 
3990 } // namespace oomph
3991 
3992 #endif
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
RefineableSolidTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameter.
virtual ~RefineableSolidTriangleMesh()
Empty Destructor.
RefineableSolidTriangleMesh(const Vector< double > &target_area, TriangulateIO &triangulate_io, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false, const bool &allow_automatic_creation_of_vertices_on_boundaries=true, OomphCommunicator *comm_pt=0)
Build mesh from specified triangulation and associated target areas for elements in it.
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
unsigned & nbin_y_for_area_transfer()
Read/write access to number of bins in the y-direction when transferring target areas by bin method....
void disable_iterative_solver_for_projection()
Enables the use of an iterative solver for the projection problem.
std::map< unsigned, std::set< Vector< double > > > Boundary_connections_pt
A map that stores the vertices that receive connections, they are identified by the boundary number t...
MeshUpdateFctPt & mesh_update_fct_pt()
Access to function pointer to function that updates the mesh following the snapping of boundary nodes...
unsigned Print_timings_level_adaptation
The printing level for adaptation.
unsigned try_to_add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
Check if necessary to add the node as haloed or if it has been previously added to the haloed scheme.
double compute_area_target(const Vector< double > &elem_error, Vector< double > &target_area)
Compute target area based on the element's error and the error target; return minimum angle (in degre...
void enable_timings_tranfering_target_areas()
Enables info. and timings for tranferring of target areas.
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors.
virtual ~RefineableTriangleMesh()
Empty Destructor.
void set_print_level_timings_load_balance(const unsigned &print_level)
Sets the printing level of timings for load balance.
void initialise_adaptation_data()
Helper function to initialise data associated with adaptation.
void enable_print_timings_load_balance(const unsigned &print_level=1)
Enables printing of timings for load balance.
void update_polyline_representation_from_restart()
Method used to update the polylines representation after restart.
void initialise_boundary_refinement_data()
Set all the flags to true (the default values)
void enable_boundary_unrefinement_constrained_by_target_areas()
Enable/disable unrefinement/refinement methods for original boundaries.
void enable_timings_projection()
Enables info. and timings for projection.
bool Do_shared_boundary_unrefinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
double Max_element_size
Max permitted element size.
void set_print_level_timings_adaptation(const unsigned &print_level)
Sets the printing level of timings for adaptation.
void disable_print_timings_load_balance()
Disables printing of timings for load balance.
unsigned try_to_add_node_pt_load_balance(Vector< Node * > &new_nodes_on_domain, Node *&node_pt)
Check if necessary to add the node to the new domain or if it has been already added.
unsigned Nbin_x_for_area_transfer
Number of bins in the x-direction when transferring target areas by bin method. Only used if we don't...
void enable_shared_boundary_unrefinement_constrained_by_target_areas()
Enable/disable unrefinement/refinement methods for shared boundaries.
void reestablish_distribution_info_for_restart(OomphCommunicator *comm_pt, std::istream &restart_file)
Used to re-establish any additional info. related with the distribution after a re-starting for trian...
unsigned try_to_add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
Check if necessary to add the element as haloed or if it has been previously added to the haloed sche...
void disable_timings_projection()
Disables info. and timings for projection.
void refine_uniformly(DocInfo &doc_info)
Refine mesh uniformly and doc process.
RefineableTriangleMesh(const Vector< double > &target_area, TriangulateIO &triangulate_io, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false, const bool &allow_automatic_creation_of_vertices_on_boundaries=true, OomphCommunicator *comm_pt=0)
Build mesh from specified triangulation and associated target areas for elements in it NOTE: This is ...
MeshUpdateFctPt Mesh_update_fct_pt
Function pointer to function that updates the mesh following the snapping of boundary nodes to the bo...
void enable_print_timings_adaptation(const unsigned &print_level=1)
Enables printing of timings for adaptation.
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds.
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles.
double & max_element_size()
Max element size allowed during adaptation.
InternalHolePointUpdateFctPt Internal_hole_point_update_fct_pt
Function pointer to function that can be set to update the position of the central point in internal ...
void disable_projection()
Disables the solution projection step during adaptation.
unsigned try_to_add_element_pt_load_balance(Vector< FiniteElement * > &new_elements_on_domain, FiniteElement *&ele_pt)
Check if necessary to add the element to the new domain or if it has been previously added.
bool Use_iterative_solver_for_projection
Flag to indicate whether to use or not an iterative solver (CG with diagonal preconditioned) for the ...
unsigned max_sample_points_for_limited_locate_zeta_during_target_area_transfer()
Read/write access to number of sample points from which we try to locate zeta by Newton method when t...
unsigned & nbin_x_for_area_transfer()
Read/write access to number of bins in the x-direction when transferring target areas by bin method....
void disable_shared_boundary_unrefinement_constrained_by_target_areas()
const bool boundary_connections(const unsigned &b, const unsigned &c, std::set< Vector< double >> &vertices)
Verifies if the given boundary receives a connection, and if that is the case then returns the list o...
RefineableTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
Build mesh, based on the polyfiles.
bool Disable_projection
Enable/disable solution projection during adaptation.
double Min_permitted_angle
Min angle before remesh gets triggered.
RefineableTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
Node * sorted_shared_boundary_node_pt(unsigned &b, unsigned &i)
bool Do_boundary_unrefinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
double & min_permitted_angle()
Min angle before remesh gets triggered.
void doc_adaptivity_targets(std::ostream &outfile)
Doc the targets for mesh adaptation.
unsigned Max_sample_points_for_limited_locate_zeta_during_target_area_transfer
Default value for max. number of sample points used for locate_zeta when transferring target areas us...
Vector< Node * > sorted_shared_boundary_node_pt(unsigned &b)
void enable_projection()
Enables the solution projection step during adaptation.
void enable_iterative_solver_for_projection()
Enables the use of an iterative solver for the projection problem.
void disable_print_timings_adaptation()
Disables printing of timings for adaptation.
bool Do_boundary_refinement_constrained_by_target_areas
Flag that enables or disables boundary refinement (true by default)
std::map< unsigned, Vector< Node * > > Sorted_shared_boundary_node_pt
Stores the nodes in the boundaries in the same order in all the processors Sorted_shared_boundary_nod...
double Min_element_size
Min permitted element size.
bool Do_shared_boundary_refinement_constrained_by_target_areas
Flag that enables or disables boundary unrefinement (true by default)
bool Print_timings_projection
Enable/disable printing timings for projection.
void disable_timings_tranfering_target_areas()
Disables info. and timings for tranferring of target areas.
bool Print_timings_transfering_target_areas
Enable/disable printing timings for transfering target areas.
double & min_element_size()
Min element size allowed during adaptation.
unsigned Print_timings_level_load_balance
The printing level for load balance.
unsigned nsorted_shared_boundary_node(unsigned &b)
unsigned unrefine_uniformly()
Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarses...
unsigned Nbin_y_for_area_transfer
Number of bins in the y-direction when transferring target areas by bin method. Only used if we don't...
InternalHolePointUpdateFctPt & internal_hole_point_update_fct_pt()
Access to function pointer to can be used to generate the internal point for the ihole-th hole.
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
Vector< std::string > Flat_packed_unsigneds_string
Temporary vector of strings to enable full annotation of RefineableTriangleMesh comms.
////////////////////////////////////////////////////////////////// //////////////////////////////////...
SolidTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on closed curve that specifies the outer boundary of the domain and any number of i...
SolidTriangleMesh(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
virtual ~SolidTriangleMesh()
Empty Destructor.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Vector< TriangleMeshClosedCurve * > internal_closed_curve_pt() const
Helper function for getting the internal closed boundaries.
bool Boundary_refinement
Do not allow refinement of nodes on the boundary.
TriangleMeshClosedCurve * outer_boundary_pt(const unsigned &i) const
Helper function for getting the i-th outer boundary.
bool is_automatic_creation_of_vertices_on_boundaries_allowed()
Returns the status of the variable Allow_automatic_creation_of_vertices_on_boundaries.
TriangleMeshParameters(Vector< TriangleMeshClosedCurve * > &outer_boundary_pt)
Constructor: Only takes the outer boundary, all the other parameters are stated with the specific par...
Vector< Vector< double > > extra_holes_coordinates() const
Helper function for getting the extra holes.
virtual ~TriangleMeshParameters()
Empty destructor.
std::map< unsigned, double > & target_area_for_region()
Helper function for getting access to the region's target areas.
void set_target_area_for_region(const unsigned &i, const double &area)
Helper function to specify target area for region.
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
bool is_use_attributes() const
Helper function for getting the status of use_attributes variable.
double & element_area()
Helper function for getting access to the element area.
Vector< TriangleMeshClosedCurve * > outer_boundary_pt() const
Helper function for getting the outer boundary.
void enable_use_attributes()
Helper function for enabling the use of attributes.
void disable_use_attributes()
Helper function for disabling the use of attributes.
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed)
double element_area() const
Helper function for getting the element area.
void disable_automatic_creation_of_vertices_on_boundaries()
Disables the creation of points (by Triangle) on the outer and internal boundaries.
bool Internal_boundary_refinement
Do not allow refinement of nodes on the internal boundary.
void enable_automatic_creation_of_vertices_on_boundaries()
Enables the creation of points (by Triangle) on the outer and internal boundaries.
void disable_internal_boundary_refinement()
Helper function for disabling the use of boundary refinement.
TriangleMeshParameters(TriangleMeshClosedCurve *outer_boundary_pt)
Constructor: Only takes the outer boundary, all the other parameters are stated with the specific par...
std::map< unsigned, Vector< double > > Regions_coordinates
Store the coordinates for defining extra regions The key on the map is the region id.
Vector< Vector< double > > & extra_holes_coordinates()
Helper function for getting access to the extra holes.
Vector< TriangleMeshClosedCurve * > & internal_closed_curve_pt()
Helper function for getting access to the internal closed boundaries.
void add_region_coordinates(const unsigned &i, Vector< double > &region_coordinates)
Helper function for getting the extra regions.
TriangleMeshClosedCurve *& outer_boundary_pt(const unsigned &i)
Helper function for getting access to the i-th outer boundary.
TriangleMeshParameters()
Constructor: Takes nothing and initializes the other parameters to the default ones.
Vector< TriangleMeshClosedCurve * > Internal_closed_curve_pt
Internal closed boundaries.
void set_communicator_pt(OomphCommunicator *comm_pt)
Function to set communicator (mesh is then assumed to be distributed)
Vector< TriangleMeshClosedCurve * > & outer_boundary_pt()
Helper function for getting access to the outer boundary.
Vector< TriangleMeshOpenCurve * > & internal_open_curves_pt()
Helper function for getting access to the internal open boundaries.
Vector< Vector< double > > Extra_holes_coordinates
Store the coordinates for defining extra holes.
Vector< TriangleMeshClosedCurve * > Outer_boundary_pt
The outer boundary.
bool Allow_automatic_creation_of_vertices_on_boundaries
Allows automatic creation of vertices along boundaries by Triangle.
Vector< TriangleMeshOpenCurve * > internal_open_curves_pt() const
Helper function for getting the internal open boundaries.
std::map< unsigned, double > Regions_areas
Target areas for regions; defaults to 0.0 which (luckily) implies "no specific target area" for trian...
bool Use_attributes
Define the use of attributes (regions)
OomphCommunicator * Comm_pt
Pointer to communicator – set to NULL if mesh is not distributed Required to pass it to new distribut...
void enable_internal_boundary_refinement()
Helper function for enabling the use of boundary refinement.
std::map< unsigned, Vector< double > > & regions_coordinates()
Helper function for getting access to the regions coordinates.
void disable_boundary_refinement()
Helper function for disabling the use of boundary refinement.
bool is_boundary_refinement_allowed() const
Helper function for getting the status of boundary refinement.
void enable_boundary_refinement()
Helper function for enabling the use of boundary refinement.
bool is_internal_boundary_refinement_allowed() const
Helper function for getting the status of boundary refinement.
Vector< TriangleMeshOpenCurve * > Internal_open_curves_pt
Internal boundaries.
double Element_area
The element are when calling triangulate external routine.
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
TimeStepper * Time_stepper_pt
Timestepper used to build elements.
const unsigned nshared_boundary_element(const unsigned &b)
void re_assign_initial_zeta_values_for_internal_boundary(const unsigned &b, Vector< std::list< FiniteElement * >> &old_segment_sorted_ele_pt, std::map< FiniteElement *, bool > &old_is_inverted)
Re-assign the boundary segments initial zeta (arclength) value for those internal boundaries that wer...
Vector< Vector< Vector< unsigned > > > Shared_boundaries_ids
Stores the boundaries ids created by the interaction of two processors Shared_boundaries_ids[iproc][j...
TriangleScaffoldMesh * Tmp_mesh_pt
Temporary scaffold mesh.
std::map< unsigned, Vector< unsigned > > & shared_boundary_from_processors()
Return the association of the shared boundaries with the processors.
Vector< TriangleMeshPolyLine * > & boundary_subpolylines(const unsigned &b)
Gets the vector of auxiliar polylines that will represent the given boundary (useful only when the bo...
bool is_node_on_shared_boundary(const unsigned &b, Node *const &node_pt)
Is the node on the shared boundary.
bool Use_attributes
Boolean flag to indicate whether to use attributes or not (required for multidomain meshes)
void compute_boundary_segments_connectivity_and_initial_zeta_values(const unsigned &b)
Compute the boundary segments connectivity for those boundaries that were splited during the distribu...
std::map< unsigned, unsigned > Shared_boundary_overlaps_internal_boundary
Stores information about those shared boundaries that lie over or over a segment of an internal bound...
void build_triangulateio(const std::string &poly_file_name, TriangulateIO &triangulate_io, bool &use_attributes)
Helper function to create TriangulateIO object (return in triangulate_io) from the ....
virtual ~TriangleMesh()
Destructor.
std::map< unsigned, unsigned > & shared_boundary_overlaps_internal_boundary()
Gets the storage that indicates if a shared boundary is part of an internal boundary.
const bool boundary_was_splitted(const unsigned &b)
Helper function to verify if a given boundary was splitted in the distribution process.
void flush_shared_boundary_node(const unsigned &b)
Flush the boundary nodes associated to the shared boundary b.
void create_shared_polylines_connections()
Establish the connections of the polylines previously marked as having connections....
std::map< unsigned, double > Regions_areas
Target areas for regions; defaults to 0.0 which (luckily) implies "no specific target area" for trian...
void update_triangulateio()
Update the triangulateio object to the current nodal positions.
const unsigned initial_shared_boundary_id()
The initial boundary id for shared boundaries.
unsigned Initial_shared_boundary_id
The initial boundary id for shared boundaries.
const unsigned shared_boundaries_ids(const unsigned &p, const unsigned &q, const unsigned &i) const
Node *& boundary_segment_node_pt(const unsigned &b, const unsigned &s, const unsigned &n)
Return pointer to node n on boundary b.
void break_loops_on_shared_polyline_load_balance_helper(const unsigned &initial_shd_bnd_id, std::list< Node * > &input_nodes, Vector< FiniteElement * > &input_boundary_element_pt, Vector< FiniteElement * > &input_boundary_face_element_pt, Vector< int > &input_face_index_element, const int &input_connect_to_the_left, const int &input_connect_to_the_right, Vector< std::list< Node * >> &output_sorted_nodes_pt, Vector< Vector< FiniteElement * >> &output_boundary_element_pt, Vector< Vector< FiniteElement * >> &output_boundary_face_element_pt, Vector< Vector< int >> &output_face_index_element, Vector< int > &output_connect_to_the_left, Vector< int > &output_connect_to_the_right)
Break any possible loop created by the sorted list of nodes that is used to create a new shared polyl...
Vector< Vector< unsigned > > shared_boundaries_ids(const unsigned &p) const
void shared_boundaries_in_this_processor(Vector< unsigned > &shared_boundaries_in_this_processor)
Get the shared boundaries ids living in the current processor.
void synchronize_boundary_coordinates(const unsigned &b)
In charge of sinchronize the boundary coordinates for internal boundaries that were split as part of ...
void update_triangulateio(Vector< Vector< double >> &internal_point)
Update the TriangulateIO object to the current nodal position and the centre hole coordinates.
void compute_holes_left_by_halo_elements_helper(Vector< Vector< double >> &output_holes_coordinates)
Compute the holes left by the halo elements, those adjacent to the shared boundaries.
void identify_boundary_segments_and_assign_initial_zeta_values(const unsigned &b, Vector< FiniteElement * > &input_face_ele_pt, const bool &is_internal_boundary, std::map< FiniteElement *, FiniteElement * > &face_to_bulk_element_pt)
Identify the segments from the old mesh (original mesh) in the new mesh (this) and assign initial and...
std::map< unsigned, Vector< int > > Face_index_at_shared_boundary
For the e-th finite element on shared boundary b, this is the index of the face that lies along that ...
std::map< unsigned, Vector< unsigned > > Shared_boundary_from_processors
Stores the processors involved in the generation of a shared boundary, in 2D two processors give rise...
const bool shared_boundary_overlaps_internal_boundary(const unsigned &shd_bnd_id)
Checks if the shared boundary overlaps an internal boundary.
const unsigned nshared_boundary_node(const unsigned &b)
std::map< unsigned, Vector< TriangleMeshPolyLine * > > Boundary_subpolylines
The polylines that will temporary represent the boundary that was splitted in the distribution proces...
const unsigned nshared_boundary_overlaps_internal_boundary()
Get the number of shared boundaries overlaping internal boundaries.
void flush_shared_boundary_element(const unsigned &b)
virtual void reset_boundary_element_info(Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned >> &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
Virtual function to perform the reset boundary elements info routines. Generally used after load bala...
const unsigned nshared_boundary_polyline(const unsigned &p, const unsigned &c) const
void add_shared_boundary_node(const unsigned &b, Node *node_pt)
Add the node the shared boundary.
Vector< Vector< unsigned > > & shared_boundaries_ids(const unsigned &p)
const unsigned nshared_boundary_curves(const unsigned &p) const
const int check_connections_of_polyline_nodes(std::set< FiniteElement * > &element_in_processor_pt, const int &root_edge_bnd_id, std::map< std::pair< Node *, Node * >, bool > &overlapped_face, std::map< unsigned, std::map< Node *, bool >> &node_on_bnd_not_overlapped_by_shd_bnd, std::list< Node * > &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...
void break_loops_on_shared_polyline_helper(const unsigned &initial_shd_bnd_id, std::list< Node * > &input_nodes, Vector< FiniteElement * > &input_boundary_element_pt, Vector< int > &input_face_index_element, const int &input_connect_to_the_left, const int &input_connect_to_the_right, Vector< std::list< Node * >> &output_sorted_nodes_pt, Vector< Vector< FiniteElement * >> &output_boundary_element_pt, Vector< Vector< int >> &output_face_index_element, Vector< int > &output_connect_to_the_left, Vector< int > &output_connect_to_the_right)
Break any possible loop created by the sorted list of nodes that is used to create a new shared polyl...
void create_tmp_polygons_helper(Vector< Vector< TriangleMeshPolyLine * >> &polylines_pt, Vector< TriangleMeshPolygon * > &polygons_pt)
Take the polylines from the shared boundaries and create temporary polygon representations of the dom...
std::map< unsigned, bool > Boundary_was_splitted
Flag to indicate if a polyline has been splitted during the distribution process, the boundary id of ...
const unsigned read_unsigned_line_helper(std::istream &read_file)
bool Triangulateio_exists
Boolean defining if Triangulateio object has been built or not.
void re_scale_re_assigned_initial_zeta_values_for_internal_boundary(const unsigned &b)
Re-scale the re-assigned zeta values for the boundary nodes, apply only for internal boundaries.
const unsigned nshared_boundaries(const unsigned &p, const unsigned &q) const
Access functions to boundaries shared with processors.
Vector< Vector< double > > Original_extra_holes_coordinates
Backup the original extra holes coordinates.
virtual void load_balance(const Vector< unsigned > &target_domain_for_local_non_halo_element)
Virtual function to perform the load balance routines.
std::map< unsigned, Vector< Node * > > Shared_boundary_node_pt
Stores the boundary nodes adjacent to the shared boundaries, these nodes are a subset of the halo and...
Vector< TriangleMeshPolyLine * > & shared_boundary_polyline_pt(const unsigned &p, const unsigned &c)
void create_tmp_open_curves_helper(Vector< Vector< TriangleMeshPolyLine * >> &sorted_open_curves_pt, Vector< TriangleMeshPolyLine * > &unsorted_shared_to_internal_poly_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt)
Take the polylines from the original open curves and created new temporaly representations of open cu...
void sort_polylines_helper(Vector< TriangleMeshPolyLine * > &unsorted_polylines_pt, Vector< Vector< TriangleMeshPolyLine * >> &sorted_polylines_pt)
Sorts the polylines so they be continuous and then we can create a closed or open curve from them.
void create_distributed_domain_representation(Vector< TriangleMeshPolygon * > &polygons_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt)
Creates the distributed domain representation. Joins the original boundaires, shared boundaries and c...
Vector< Vector< Node * > > & boundary_segment_node_pt(const unsigned &b)
Return direct access to nodes associated with a boundary but sorted in segments.
void add_shared_boundary_element(const unsigned &b, FiniteElement *ele_pt)
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
coord2_t cross(const Point &O, const Point &A, const Point &B)
2D cross product of OA and OB vectors, i.e. z-component of their 3D cross product....
void operator=(const TriangleMesh &)=delete
Broken assignment operator.
Vector< Vector< Vector< TriangleMeshPolyLine * > > > Shared_boundary_polyline_pt
Stores the polyline representation of the shared boundaries Shared_boundary_polyline_pt[iproc][ncurve...
Vector< Vector< Vector< unsigned > > > shared_boundaries_ids() const
TriangleMeshPolyLine * shared_boundary_polyline_pt(const unsigned &p, const unsigned &c, const unsigned &i) const
Vector< unsigned > & shared_boundary_from_processors(const unsigned &b)
void create_shared_boundaries(OomphCommunicator *comm_pt, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, const Vector< FiniteElement * > &backed_up_f_el_pt, std::map< Data *, std::set< unsigned >> &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status)
Creates the shared boundaries.
TriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
Node * shared_boundary_node_pt(const unsigned &b, const unsigned &n)
std::vector< Point > convex_hull(std::vector< Point > P)
Returns a list of points on the convex hull in counter-clockwise order. Note: the last point in the r...
const unsigned nshared_boundaries() const
Vector< unsigned > oomph_vertex_nodes_id()
Return the vector that contains the oomph-lib node number for all vertex nodes in the TriangulateIO r...
void remesh_from_internal_triangulateio()
Completely regenerate the mesh from the trianglateio structure.
Vector< unsigned > Oomph_vertex_nodes_id
Vector storing oomph-lib node number for all vertex nodes in the TriangulateIO representation of the ...
TriangleMesh(const std::string &poly_file_name, const double &element_area, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
Build mesh from poly file, with specified target area for all elements.
const unsigned final_shared_boundary_id()
The final boundary id for shared boundaries.
virtual void reestablish_distribution_info_for_restart(OomphCommunicator *comm_pt, std::istream &restart_file)
Virtual function used to re-establish any additional info. related with the distribution after a re-s...
const unsigned nboundary_subpolylines(const unsigned &b)
Gets the number of subpolylines that create the boundarya (useful only when the boundary is marked as...
const unsigned shared_boundary_overlapping_internal_boundary(const unsigned &shd_bnd_id)
Gets the boundary id of the internal boundary that the shared boundary lies on.
TriangleMesh(const TriangleMesh &dummy)=delete
Broken copy constructor.
bool First_time_compute_holes_left_by_halo_elements
Flag to know if it is the first time we are going to compute the holes left by the halo elements.
void create_shared_polyline(const unsigned &my_rank, const unsigned &shd_bnd_id, const unsigned &iproc, const unsigned &jproc, std::list< Node * > &sorted_nodes, const int &root_edge_bnd_id, Vector< FiniteElement * > &bulk_bnd_ele_pt, Vector< int > &face_index_ele, Vector< Vector< TriangleMeshPolyLine * >> &unsorted_polylines_pt, const int &connect_to_the_left_flag, const int &connect_to_the_right_flag)
Create the shared polyline and fill the data structured that keep all the information associated with...
void dump_distributed_info_for_restart(std::ostream &dump_file)
Used to dump info. related with distributed triangle meshes.
void get_halo_elements_on_all_procs(const unsigned &nproc, const Vector< unsigned > &element_domain, const Vector< GeneralisedElement * > &backed_up_el_pt, std::map< Data *, std::set< unsigned >> &processors_associated_with_data, const bool &overrule_keep_as_halo_element_status, std::map< GeneralisedElement *, unsigned > &element_to_global_index, Vector< Vector< Vector< GeneralisedElement * >>> &output_halo_elements_pt)
Creates the halo elements on all processors Gets the halo elements on all processors,...
const bool boundary_marked_as_shared_boundary(const unsigned &b, const unsigned &isub)
Returns the value that indicates if a subpolyline of a given boundary continues been used as internal...
void set_mesh_level_time_stepper(TimeStepper *const &time_stepper_pt, const bool &preserve_existing_data)
Overload set_mesh_level_time_stepper so that the stored time stepper now corresponds to the new times...
unsigned Final_shared_boundary_id
The final boundary id for shared boundaries.
Vector< Vector< Vector< unsigned > > > & shared_boundaries_ids()
void read_distributed_info_for_restart(std::istream &restart_file)
Used to read info. related with distributed triangle meshes.
void select_boundary_face_elements(Vector< FiniteElement * > &face_el_pt, const unsigned &b, bool &is_internal_boundary, std::map< FiniteElement *, FiniteElement * > &face_to_bulk_element_pt)
Select face element from boundary using the criteria to decide which of the two face elements should ...
void output_boundary_coordinates(const unsigned &b, std::ostream &outfile)
Output the nodes on the boundary and their respective boundary coordinates(into separate tecplot zone...
void generic_constructor(Vector< TriangleMeshPolygon * > &outer_boundary_pt, Vector< TriangleMeshPolygon * > &internal_polygon_pt, Vector< TriangleMeshOpenCurve * > &open_polylines_pt, const double &element_area, Vector< Vector< double >> &extra_holes_coordinates, std::map< unsigned, Vector< double >> &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 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 add_face_index_at_shared_boundary(const unsigned &b, const unsigned &i)
Vector< unsigned > & shared_boundaries_ids(const unsigned &p, const unsigned &q)
int face_index_at_shared_boundary(const unsigned &b, const unsigned &e)
void flush_shared_boundary_node()
Flush ALL the shared boundary nodes.
void get_shared_boundaries_overlapping_internal_boundary(const unsigned &internal_bnd_id, Vector< unsigned > &shd_bnd_ids)
Gets the shared boundaries ids that overlap the given internal boundary.
bool triangulateio_exists()
Boolean defining if Triangulateio object has been built or not.
FiniteElement * shared_boundary_element_pt(const unsigned &b, const unsigned &e)
void get_element_edges_on_boundary(std::map< std::pair< Node *, Node * >, unsigned > &element_edges_on_boundary)
Get the element edges (pair of nodes, edges) that lie on a boundary (used to mark shared boundaries t...
TriangleMesh()
Empty constructor.
void update_holes_information_helper(Vector< TriangleMeshPolygon * > &polygons_pt, Vector< Vector< double >> &output_holes_coordinates)
Keeps those vertices that define a hole, those that are inside closed internal boundaries in the new ...
std::map< unsigned, Vector< FiniteElement * > > Shared_boundary_element_pt
Stores the boundary elements adjacent to the shared boundaries, these elements are a subset of the ha...
Vector< Node * > & boundary_segment_node_pt(const unsigned &b, const unsigned &s)
Return direct access to nodes associated with a segment of a given boundary.
void create_polylines_from_halo_elements_helper(const Vector< unsigned > &element_domain, std::map< GeneralisedElement *, unsigned > &element_to_global_index, std::set< FiniteElement * > &element_in_processor_pt, Vector< Vector< Vector< GeneralisedElement * >>> &input_halo_elements, std::map< std::pair< Node *, Node * >, unsigned > &elements_edges_on_boundary, Vector< Vector< Vector< TriangleMeshPolyLine * >>> &output_polylines_pt)
Creates polylines from the intersection of halo elements on all processors. The new polylines define ...
Vector< unsigned > shared_boundaries_ids(const unsigned &p, const unsigned &q) const
////////////////////////////////////////////////////////////////////// //////////////////////////////...
void triangulate(char *triswitches, struct oomph::TriangulateIO *in, struct oomph::TriangulateIO *out, struct oomph::TriangulateIO *vorout)
bool operator<(const Point &p) const