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// Version 1.0; svn revision $LastChangedRevision: 1299 $
7 //LIC//
8 //LIC// $LastChangedDate: 2017-10-06 07:40:18 +0100 (Fri, 06 Oct 2017) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 #ifndef OOMPH_TRIANGLE_MESH_HEADER
31 #define OOMPH_TRIANGLE_MESH_HEADER
32 // Config header generated by autoconfig
33 #ifdef HAVE_CONFIG_H
34 #include <oomph-lib-config.h>
35 #endif
36 
37 #ifdef OOMPH_HAS_MPI
38 
39 // MPI headers
40 #include <mpi.h>
41 #endif
42 
43 //Standards
44 #include<float.h>
45 #include <iostream>
46 #include <fstream>
47 #include <string.h>
48 #include <iomanip>
49 
50 #include "../generic/problem.h"
51 #include "../generic/triangle_scaffold_mesh.h"
52 #include "../generic/triangle_mesh.h"
53 #include "../generic/refineable_mesh.h"
54 #include "../rigid_body/immersed_rigid_body_elements.h"
55 
56 namespace oomph
57 {
58 
59 #ifdef OOMPH_HAS_TRIANGLE_LIB
60 
61 // Interface to triangulate function
62  extern "C" {
63  void triangulate(char *triswitches, struct oomph::TriangulateIO *in,
64  struct oomph::TriangulateIO *out,
65  struct oomph::TriangulateIO *vorout);
66  }
67 
68 #endif
69 
70 
71 
72 
73 ////////////////////////////////////////////////////////////////////////
74 ////////////////////////////////////////////////////////////////////////
75 ////////////////////////////////////////////////////////////////////////
76 ////////////////////////////////////////////////////////////////////////
77 ////////////////////////////////////////////////////////////////////////
78 
79 
80 
81 
82 //=========================================================================
83 /// \short Helper object for dealing with the parameters used for the
84 /// TriangleMesh objects
85 //=========================================================================
87  {
88 
89  public:
90 
91  /// Constructor: Only takes the outer boundary, all the other parameters
92  /// are stated with the specific parameters
94  : Outer_boundary_pt(outer_boundary_pt),
95  Element_area(0.2),
96  Use_attributes(false),
97  Boundary_refinement(true),
100  Comm_pt(0)
101  { }
102 
103  /// Constructor: Only takes the outer boundary, all the other parameters
104  /// are stated with the specific parameters
106  : Element_area(0.2),
107  Use_attributes(false),
108  Boundary_refinement(true),
111  Comm_pt(0)
112  {
113  Outer_boundary_pt.resize(1);
115  }
116 
117  /// Constructor: Takes nothing and initialize the other paremeters to
118  /// the default ones
120  : Element_area(0.2),
121  Use_attributes(false),
122  Boundary_refinement(true),
125  Comm_pt(0)
126  { }
127 
128  /// Empty destructor
130 
131  /// Helper function for getting the outer boundary
133  {return Outer_boundary_pt;}
134 
135  /// Helper function for getting access to the outer boundary
137  {return Outer_boundary_pt;}
138 
139  /// Helper function for getting the i-th outer boundary
141  {return Outer_boundary_pt[i];}
142 
143  /// Helper function for getting access to the i-th outer boundary
145  {return Outer_boundary_pt[i];}
146 
147  /// Helper function for getting the internal closed boundaries
149  {return Internal_closed_curve_pt;}
150 
151  /// \short Helper function for getting access to the internal
152  /// closed boundaries
154  {return Internal_closed_curve_pt;}
155 
156  /// Helper function for getting the internal open boundaries
158  {return Internal_open_curves_pt;}
159 
160  /// \short Helper function for getting access to the internal
161  /// open boundaries
163  {return Internal_open_curves_pt;}
164 
165  /// Helper function for getting the element area
166  double element_area() const {return Element_area;}
167 
168  /// Helper function for getting access to the element area
169  double &element_area(){return Element_area;}
170 
171  /// Helper function for getting the extra holes
173  {return Extra_holes_coordinates;}
174 
175  /// Helper function for getting access to the extra holes
177  {return Extra_holes_coordinates;}
178 
179  /// Helper function for getting the extra regions
180  void add_region_coordinates(const unsigned &i,
181  Vector<double> &region_coordinates)
182  {
183  // Verify if not using the default region number (zero)
184  if (i == 0)
185  {
186  std::ostringstream error_message;
187  error_message << "Please use another region id different from zero.\n"
188  << "It is internally used as the default region number.\n";
189  throw OomphLibError(error_message.str(),
190  OOMPH_CURRENT_FUNCTION,
191  OOMPH_EXCEPTION_LOCATION);
192  }
193 
194  // First check if the region with the specified id does not already exist
195  std::map<unsigned, Vector<double> >::iterator it;
196  it = Regions_coordinates.find(i);
197 
198  // If it is already a region defined with that id throw an error
199  if (it != Regions_coordinates.end())
200  {
201  std::ostringstream error_message;
202  error_message << "The region id ("<<i<<") that you are using for"
203  << "defining\n"
204  << "your region is already in use. Use another\n"
205  << "region id and verify that you are not re-using\n"
206  <<" previously defined regions ids\n"<<std::endl;
207  OomphLibWarning(error_message.str(),
208  OOMPH_CURRENT_FUNCTION,
209  OOMPH_EXCEPTION_LOCATION);
210  }
211 
212  // If it does not exist then create the map
213  Regions_coordinates[i] = region_coordinates;
214 
215  // Automatically set the using of attributes to enable
217  }
218 
219  /// Helper function for getting access to the regions coordinates
220  std::map<unsigned, Vector<double> >&regions_coordinates()
221  {return Regions_coordinates;}
222 
223  /// Helper function to specify target area for region
224  void set_target_area_for_region(const unsigned &i, const double& area)
225  {
226  Regions_areas[i] = area;
227  }
228 
229  /// Helper function for getting access to the region's target areas
230  std::map<unsigned, double >&target_area_for_region()
231  {return Regions_areas;}
232 
233  /// \short Helper function for enabling the use of attributes
235 
236  /// \short Helper function for disabling the use of attributes
238 
239  /// \short Helper function for getting the status of use_attributes
240  /// variable
241  bool is_use_attributes() const {return Use_attributes;}
242 
243  /// \short Helper function for enabling the use of boundary refinement
245 
246  /// Boolean to indicate if Mesh has been distributed
247  bool is_mesh_distributed() const {return (Comm_pt!=0);}
248 
249  /// Function to set communicator (mesh is then assumed to be distributed)
251  {Comm_pt=comm_pt;}
252 
253  /// Read-only access fct to communicator (Null if mesh is not distributed)
255  {return Comm_pt;}
256 
257  /// \short Helper function for disabling the use of boundary refinement
259 
260  /// \short Helper function for getting the status of boundary refinement
262 
263  /// \short Helper function for enabling the use of boundary refinement
266 
267  /// \short Helper function for disabling the use of boundary refinement
270 
271  /// \short Helper function for getting the status of boundary refinement
274 
275  /// Enables the creation of points (by Triangle) on the outer and
276  /// internal boundaries
278  {
280  }
281 
282  /// Disables the creation of points (by Triangle) on the outer and
283  /// internal boundaries
285  {
287  }
288 
289  /// Returns the status of the variable
290  /// Allow_automatic_creation_of_vertices_on_boundaries
292  {
294  }
295 
296  protected:
297 
298  /// The outer boundary
300 
301  /// Internal closed boundaries
303 
304  /// Internal boundaries
306 
307  /// The element are when calling triangulate external routine
308  double Element_area;
309 
310  /// Store the coordinates for defining extra holes
312 
313  /// \short Store the coordinates for defining extra regions
314  /// The key on the map is the region id
315  std::map<unsigned, Vector<double> > Regions_coordinates;
316 
317  /// \short Target areas for regions; defaults to 0.0 which (luckily)
318  /// implies "no specific target area" for triangle!
319  std::map<unsigned, double> Regions_areas;
320 
321  /// Define the use of attributes (regions)
323 
324  /// Do not allow refinement of nodes on the boundary
326 
327  /// Do not allow refinement of nodes on the internal boundary
329 
330  /// Allows automatic creation of vertices along boundaries by
331  /// Triangle
333 
334  /// Pointer to communicator -- set to NULL if mesh is not distributed
335  /// Required to pass it to new distributed meshes created at the
336  /// adaptation stage
338 
339  };
340 
341 
342 ////////////////////////////////////////////////////////////////////////
343 ////////////////////////////////////////////////////////////////////////
344 ////////////////////////////////////////////////////////////////////////
345 ////////////////////////////////////////////////////////////////////////
346 ////////////////////////////////////////////////////////////////////////
347 
348 
349 //============start_of_triangle_class===================================
350 /// Triangle mesh build with the help of the scaffold mesh coming
351 /// from the triangle mesh generator Triangle.
352 /// http://www.cs.cmu.edu/~quake/triangle.html
353 //======================================================================
354  template<class ELEMENT>
355  class TriangleMesh : public virtual TriangleMeshBase
356  {
357  public:
358 
359  /// \short Empty constructor
361  {
362 #ifdef OOMPH_HAS_TRIANGLE_LIB
363  // Using this constructor no Triangulateio object is built
364  Triangulateio_exists=false;
365  // By default allow the automatic creation of vertices along the
366  // boundaries by Triangle
368 #ifdef OOMPH_HAS_MPI
369  // Initialize the flag to indicate this is the first time to
370  // compute the holes left by the halo elements
372 #endif // #ifdef OOMPH_HAS_MPI
373 
374 #endif
375 
376  // Mesh can only be built with 2D Telements.
377  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(2);
378  }
379 
380  /// \short Constructor with the input files
381  TriangleMesh(const std::string& node_file_name,
382  const std::string& element_file_name,
383  const std::string& poly_file_name,
384  TimeStepper* time_stepper_pt=
386  const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
387  {
388  // Mesh can only be built with 2D Telements.
389  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(2);
390 
391  // Initialize the value for allowing creation of points on boundaries
393  allow_automatic_creation_of_vertices_on_boundaries;
394 
395 #ifdef OOMPH_HAS_MPI
396  // Initialize the flag to indicate this is the first time to
397  // compute the holes left by the halo elements
399 #endif // #ifdef OOMPH_HAS_MPI
400 
401  // Store Timestepper used to build elements
402  Time_stepper_pt=time_stepper_pt;
403 
404  // Check if we should use attributes. This is set to true if the .poly file
405  // specifies regions
406  bool should_use_attributes=false;
407 
408 #ifdef OOMPH_HAS_TRIANGLE_LIB
409  // Using this constructor build the triangulatio
411  element_file_name,
412  poly_file_name,
414  should_use_attributes);
415 
416  // Record that the triangulateio object has been created
418 #endif
419 
420  // Store the attributes
421  Use_attributes=should_use_attributes;
422 
423  // Build scaffold
424  this->Tmp_mesh_pt= new
425  TriangleScaffoldMesh(node_file_name,
426  element_file_name,
427  poly_file_name);
428 
429  // Convert mesh from scaffold to actual mesh
430  build_from_scaffold(time_stepper_pt,should_use_attributes);
431 
432  // kill the scaffold
433  delete this->Tmp_mesh_pt;
434  this->Tmp_mesh_pt=0;
435 
436  // Setup boundary coordinates for boundaries
437  unsigned nb=nboundary();
438  for (unsigned b=0;b<nb;b++)
439  {
440  this->template setup_boundary_coordinates<ELEMENT>(b);
441  }
442  }
443 
444  protected:
445 
446 #ifdef OOMPH_HAS_TRIANGLE_LIB
447 
448  public:
449 
450  /// \short Build mesh, based on the specifications on
451  /// TriangleMeshParameters
452  TriangleMesh(TriangleMeshParameters &triangle_mesh_parameters,
453  TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper)
454  {
455 
456  // Store the region target areas
457  Regions_areas=triangle_mesh_parameters.target_area_for_region();
458 
459  // Mesh can only be built with 2D Telements.
460  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(2);
461 
462  // Initialize the value for allowing creation of points on boundaries
465 
466  // Store Timestepper used to build elements
467  Time_stepper_pt=time_stepper_pt;
468 
469 #ifdef OOMPH_HAS_MPI
470  // Initialize the flag to indicate this is the first time to
471  // compute the holes left by the halo elements
473 #endif // #ifdef OOMPH_HAS_MPI
474 
475  // ********************************************************************
476  // First part - Get polylines representations
477  // ********************************************************************
478 
479  // Create the polyline representation of all the boundaries and
480  // then create the mesh by calling to "generic_constructor()"
481 
482  // Initialise highest boundary id
483  unsigned max_boundary_id = 0;
484 
485  // *****************************************************************
486  // Part 1.1 - Outer boundary
487  // *****************************************************************
488  // Get the representation of the outer boundaries from the
489  // TriangleMeshParameters object
490  Vector<TriangleMeshClosedCurve *> outer_boundary_pt =
491  triangle_mesh_parameters.outer_boundary_pt();
492 
493 #ifdef PARANOID
494  // Verify that the outer_boundary_object_pt has been set
495  if (outer_boundary_pt.size()==0)
496  {
497  std::stringstream error_message;
498  error_message
499  << "There are no outer boundaries defined.\n"
500  << "Verify that you have specified the outer boundaries in the\n"
501  << "Triangle_mesh_parameter object\n\n";
502  throw OomphLibError(error_message.str(),
503  OOMPH_CURRENT_FUNCTION,
504  OOMPH_EXCEPTION_LOCATION);
505  } // if (outer_boundary_pt!=0)
506 #endif
507 
508  // Find the number of outer closed curves
509  unsigned n_outer_boundaries = outer_boundary_pt.size();
510 
511  // Create the storage for the polygons that define the outer
512  // boundaries
513  Vector<TriangleMeshPolygon*> outer_boundary_polygon_pt(n_outer_boundaries);
514 
515  // Loop over the number of outer boundaries
516  for(unsigned i=0;i<n_outer_boundaries;++i)
517  {
518  // Get the polygon representation and compute the max boundary_id on
519  // each outer polygon. Does nothing (i.e. just returns a pointer to
520  // the outer boundary that was input) if the outer boundary is
521  // already a polygon
522  outer_boundary_polygon_pt[i] =
523  closed_curve_to_polygon_helper(outer_boundary_pt[i], max_boundary_id);
524  }
525 
526  // *****************************************************************
527  // Part 1.2 - Internal closed boundaries (possible holes)
528  // *****************************************************************
529  // Get the representation of the internal closed boundaries from the
530  // TriangleMeshParameters object
531  Vector<TriangleMeshClosedCurve *> internal_closed_curve_pt =
532  triangle_mesh_parameters.internal_closed_curve_pt();
533 
534  // Find the number of internal closed curves
535  unsigned n_internal_closed_curves = internal_closed_curve_pt.size();
536 
537  // Create the storage for the polygons that define the internal closed
538  // boundaries (again nothing happens (as above) if an internal closed
539  // curve is already a polygon)
540  Vector<TriangleMeshPolygon*> internal_polygon_pt(n_internal_closed_curves);
541 
542  // Loop over the number of internal closed curves
543  for(unsigned i=0;i<n_internal_closed_curves;++i)
544  {
545  // Get the polygon representation and compute the max boundary_id on
546  // each internal polygon
547  internal_polygon_pt[i] =
549  internal_closed_curve_pt[i], max_boundary_id);
550  }
551 
552  // *****************************************************************
553  // Part 1.3 - Internal open boundaries
554  // *****************************************************************
555  // Get the representation of open boundaries from the
556  // TriangleMeshParameteres object
557  Vector<TriangleMeshOpenCurve*> internal_open_curve_pt =
558  triangle_mesh_parameters.internal_open_curves_pt();
559 
560  //Find the number of internal open curves
561  unsigned n_internal_open_curves = internal_open_curve_pt.size();
562 
563  // Create the storage for the polylines that define the open boundaries
564  Vector<TriangleMeshOpenCurve*> internal_open_curve_poly_pt(
565  n_internal_open_curves);
566 
567  // Loop over the number of internal open curves
568  for (unsigned i = 0; i < n_internal_open_curves; i++)
569  {
570  // Get the open polyline representation and compute the max boundary_id
571  // on each open polyline (again, nothing happens if there are curve
572  // sections on the current internal open curve)
573  internal_open_curve_poly_pt[i] =
575  internal_open_curve_pt[i], max_boundary_id);
576  }
577 
578  // ********************************************************************
579  // Second part - Get associated geom objects and coordinate limits
580  // ********************************************************************
581 
582  // ***************************************************************
583  // Part 2.1 Outer boundary
584  // ***************************************************************
585  for (unsigned i = 0; i < n_outer_boundaries; i++)
586  {
588  outer_boundary_pt[i]);
589  }
590 
591  // ***************************************************************
592  // Part 2.2 - Internal closed boundaries (possible holes)
593  // ***************************************************************
594  for (unsigned i = 0; i < n_internal_closed_curves; i++)
595  {
597  internal_closed_curve_pt[i]);
598  }
599 
600  // ********************************************************************
601  // Part 2.3 - Internal open boundaries
602  // ********************************************************************
603  for (unsigned i = 0; i < n_internal_open_curves; i++)
604  {
606  internal_open_curve_pt[i]);
607  }
608 
609  // ********************************************************************
610  // Third part - Creates the TriangulateIO object by calling the
611  // "generic_constructor()" function
612  // ********************************************************************
613  // Get all the other parameters from the TriangleMeshParameters object
614  // The maximum element area
615  const double element_area =
616  triangle_mesh_parameters.element_area();
617 
618  // The holes coordinates
619  Vector<Vector<double> > extra_holes_coordinates =
620  triangle_mesh_parameters.extra_holes_coordinates();
621 
622  // The regions coordinates
623  std::map<unsigned, Vector<double> > regions =
624  triangle_mesh_parameters.regions_coordinates();
625 
626  // If we use regions then we use attributes
627  const bool use_attributes = triangle_mesh_parameters.is_use_attributes();
628 
629  const bool refine_boundary =
630  triangle_mesh_parameters.is_boundary_refinement_allowed();
631 
632  const bool refine_internal_boundary =
633  triangle_mesh_parameters.is_internal_boundary_refinement_allowed();
634 
635  if(!refine_internal_boundary && refine_boundary)
636  {
637  std::ostringstream error_stream;
638  error_stream
639  <<
640  "You have specified that Triangle may refine the outer boundary, but\n"
641  <<
642  "not internal boundaries. Triangle does not support this combination.\n"
643  <<
644  "If you do not want Triangle to refine internal boundaries, it can't\n"
645  <<
646  "refine outer boundaries either!\n"
647  << "Please either disable all boundary refinement\n"
648  << "(call TriangleMeshParameters::disable_boundary_refinement()\n"
649  << "or enable internal boundary refinement (the default)\n";
650 
651  throw OomphLibError(error_stream.str().c_str(),
652  OOMPH_CURRENT_FUNCTION,
653  OOMPH_EXCEPTION_LOCATION);
654  }
655 
656  this->generic_constructor(outer_boundary_polygon_pt,
657  internal_polygon_pt,
658  internal_open_curve_poly_pt,
659  element_area,
660  extra_holes_coordinates,
661  regions,
662  triangle_mesh_parameters.target_area_for_region(),
663  time_stepper_pt,
664  use_attributes,
665  refine_boundary,
666  refine_internal_boundary);
667 
668  // Setup boundary coordinates for boundaries
669  unsigned nb=nboundary();
670 
671 #ifdef OOMPH_HAS_MPI
672  // Before calling setup boundary coordinates check if the mesh is
673  // marked as distrbuted
674  if (triangle_mesh_parameters.is_mesh_distributed())
675  {
676  // Set the mesh as distributed by passing the communicator
677  this->set_communicator_pt(triangle_mesh_parameters.communicator_pt());
678  }
679 #endif
680 
681  for (unsigned b=0;b<nb;b++)
682  {
683  this->template setup_boundary_coordinates<ELEMENT>(b);
684  }
685 
686  // Snap it!
688  }
689 
690  /// \short Build mesh from poly file, with specified target
691  /// area for all elements.
692  TriangleMesh(const std::string& poly_file_name,
693  const double& element_area,
694  TimeStepper* time_stepper_pt=
696  const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
697  {
698  // Mesh can only be built with 2D Telements.
699  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(2);
700 
701  // Initialize the value for allowing creation of points on boundaries
703  allow_automatic_creation_of_vertices_on_boundaries;
704 
705 #ifdef OOMPH_HAS_MPI
706  // Initialize the flag to indicate this is the first time to
707  // compute the holes left by the halo elements
709 #endif // #ifdef OOMPH_HAS_MPI
710 
711  // Disclaimer
712  std::string message=
713  "This constructor hasn't been tested since last cleanup.\n";
714  OomphLibWarning(message,
715  "TriangleMesh::TriangleMesh()",
716  OOMPH_EXCEPTION_LOCATION);
717 
718  // Store Timestepper used to build elements
719  Time_stepper_pt=time_stepper_pt;
720 
721  // Create the data structures required to call the triangulate function
722  TriangulateIO triangle_in;
723 
724  // Input string for triangle
725  std::stringstream input_string_stream;
726 
727  // MH: Like everything else, this hasn't been tested!
728  // used to be input_string_stream<<"-pA -a" << element_area << "q30";
729  input_string_stream<< "-pA -a -a" << element_area << "q30";
730 
731  // Verify if creation of new points on boundaries is allowed
732  if (!this->is_creation_of_vertices_on_boundaries_allowed())
733  {
734  input_string_stream<<" -YY";
735  }
736 
737  // Convert to a *char required by the triangulate function
738  char triswitches[100];
739  sprintf(triswitches,"%s",input_string_stream.str().c_str());
740 
741  // Create a boolean to decide whether or not to use attributes.
742  // The value of this will only be changed in build_triangulateio
743  // depending on whether or not the .poly file contains regions
744  bool use_attributes=false;
745 
746  // Build the input triangulateio object from the .poly file
747  build_triangulateio(poly_file_name,triangle_in,use_attributes);
748 
749  //Store the attributes flag
750  Use_attributes=use_attributes;
751 
752  // Build the triangulateio out object
753  triangulate(triswitches,&triangle_in,&Triangulateio,0);
754 
755  // Build scaffold
757 
758  // Convert mesh from scaffold to actual mesh
759  build_from_scaffold(time_stepper_pt,use_attributes);
760 
761  // Kill the scaffold
762  delete this->Tmp_mesh_pt;
763  this->Tmp_mesh_pt=0;
764 
765  // Cleanup but leave hole alone
766  bool clear_hole_data=false;
767  TriangleHelper::clear_triangulateio(triangle_in,clear_hole_data);
768 
769  // Setup boundary coordinates for boundaries
770  unsigned nb=nboundary();
771  for (unsigned b=0;b<nb;b++)
772  {
773  this->template setup_boundary_coordinates<ELEMENT>(b);
774  }
775  }
776 
777 #endif
778 
779  /// Broken copy constructor
780  TriangleMesh(const TriangleMesh& dummy)
781  {
782  BrokenCopy::broken_copy("TriangleMesh");
783  }
784 
785  /// Broken assignment operator
786  void operator=(const TriangleMesh&)
787  {
788  BrokenCopy::broken_assign("TriangleMesh");
789  }
790 
791  /// Destructor
792  virtual ~TriangleMesh()
793  {
794 #ifdef OOMPH_HAS_TRIANGLE_LIB
796  {
798  }
799 
800  std::set<TriangleMeshCurveSection*>::iterator it_polyline;
801  for (it_polyline = Free_curve_section_pt.begin();
802  it_polyline != Free_curve_section_pt.end();
803  it_polyline++)
804  {
805  delete (*it_polyline);
806  }
807 
808  std::set<TriangleMeshPolygon*>::iterator it_polygon;
809  for (it_polygon = Free_polygon_pt.begin();
810  it_polygon != Free_polygon_pt.end();
811  it_polygon++)
812  {
813  delete (*it_polygon);
814  }
815 
816  std::set<TriangleMeshOpenCurve*>::iterator it_open_polyline;
817  for (it_open_polyline = Free_open_curve_pt.begin();
818  it_open_polyline != Free_open_curve_pt.end();
819  it_open_polyline++)
820  {
821  delete (*it_open_polyline);
822  }
823 
824 #endif
825  }
826 
827  /// \short Overload set_mesh_level_time_stepper so that the stored
828  /// time stepper now corresponds to the new timestepper
829  void set_mesh_level_time_stepper(TimeStepper* const &time_stepper_pt,
830  const bool &preserve_existing_data)
831  {
832  this->Time_stepper_pt = time_stepper_pt;
833  }
834 
835 #ifdef OOMPH_HAS_MPI
836 
837  /// \short Compute the boundary segments connectivity for those
838  /// boundaries that were splited during the distribution process
840  const unsigned& b);
841 
842  /// \short Re-assign the boundary segments initial zeta (arclength)
843  /// value for those internal boundaries that were splited during the
844  /// distribution process. Those boundaries that have one face element
845  /// at each side of the boundary
847  const unsigned& b,
848  Vector<std::list<FiniteElement*> > &old_segment_sorted_ele_pt,
849  std::map<FiniteElement*, bool> &old_is_inverted);
850 
851  /// \short Re-scale the re-assigned zeta values for the boundary
852  /// nodes, apply only for internal boundaries
854  const unsigned& b);
855 
856  /// \short Identify the segments from the old mesh (original mesh)
857  /// in the new mesh (this) and assign initial and final boundary
858  /// coordinates for the segments that create the boundary. (This is
859  /// the version called from the original mesh to identify its own
860  /// segments)
862  const unsigned& b, Vector<FiniteElement*> &input_face_ele_pt,
863  const bool &is_internal_boundary,
864  std::map<FiniteElement*,FiniteElement*> &face_to_bulk_element_pt);
865 
866  /// \short Identify the segments from the old mesh (original mesh)
867  /// in the new mesh (this) and assign initial and final boundary
868  /// coordinates for the segments that create the boundary
870  const unsigned& b, TriangleMesh<ELEMENT>* original_mesh_pt);
871 
872  /// \short In charge of sinchronize the boundary coordinates for
873  /// internal boundaries that were split as part of the distribution
874  /// process. Called after setup_boundary_coordinates() for the
875  /// original mesh only
876  void synchronize_boundary_coordinates(const unsigned& b);
877 
878  /// \short Select face element from boundary using the criteria to
879  /// decide which of the two face elements should be used on internal
880  /// boundaries
882  Vector<FiniteElement*> &face_el_pt, const unsigned &b,
883  bool &is_internal_boundary,
884  std::map<FiniteElement*,FiniteElement*> &face_to_bulk_element_pt);
885 
886  /// \short Return direct access to nodes associated with a boundary but
887  /// sorted in segments
889  {return Boundary_segment_node_pt[b];}
890 
891  /// \short Return direct access to nodes associated with a segment of
892  /// a given boundary
893  Vector<Node*> &boundary_segment_node_pt(const unsigned &b,const unsigned &s)
894  {return Boundary_segment_node_pt[b][s];}
895 
896  /// Return pointer to node n on boundary b
897  Node* &boundary_segment_node_pt(const unsigned &b, const unsigned &s,
898  const unsigned &n)
899  {return Boundary_segment_node_pt[b][s][n];}
900 
901 #endif // OOMPH_HAS_MPI
902 
903 #ifdef OOMPH_HAS_TRIANGLE_LIB
904 
905  /// \short Update the TriangulateIO object to the current nodal position
906  /// and the centre hole coordinates.
908  {
909  // Move the hole center
910  // Get number of holes
911  unsigned nhole=Triangulateio.numberofholes;
912  unsigned count_coord=0;
913  for(unsigned ihole=0;ihole<nhole;ihole++)
914  {
915  Triangulateio.holelist[count_coord]+=internal_point[ihole][0];
916  Triangulateio.holelist[count_coord+1]+=internal_point[ihole][1];
917 
918  // Increment counter
919  count_coord+=2;
920  }
921 
922  // Do the update
924  }
925 
926  /// \short Update the triangulateio object to the current nodal positions
928  {
929  // Get number of points
931  double new_x=0;
932  double new_y=0;
933 
934  // Loop over the points
935  for(unsigned inod=0;inod<nnode;inod++)
936  {
937  // Get the node Id to be updated
938  unsigned count=Oomph_vertex_nodes_id[inod];
939 
940  // Update vertices using the vertex_node_id giving for the TriangulateIO
941  // vertex enumeration the corresponding oomphlib mesh enumeration
942  Node* mesh_node_pt=this->node_pt(inod);
943  new_x=mesh_node_pt->x(0);
944  new_y=mesh_node_pt->x(1);
945  Triangulateio.pointlist[count*2] = new_x;
946  Triangulateio.pointlist[(count*2)+1] = new_y;
947  }
948  }
949 
950 #ifdef OOMPH_HAS_MPI
951  /// Used to dump info. related with distributed triangle meshes
952  void dump_distributed_info_for_restart(std::ostream &dump_file);
953 
954  const unsigned read_unsigned_line_helper(std::istream &read_file)
955  {
956  std::string input_string;
957 
958  // Read line up to termination sign
959  getline(read_file,input_string,'#');
960 
961  // Ignore rest of line
962  read_file.ignore(200,'\n');
963 
964  // Convert
965  return std::atoi(input_string.c_str());
966  }
967 
968  /// Used to read info. related with distributed triangle meshes
969  void read_distributed_info_for_restart(std::istream &restart_file);
970 
971  /// Virtual function used to re-establish any additional info. related with
972  /// the distribution after a re-starting for triangle meshes
974  OomphCommunicator* comm_pt, std::istream &restart_file)
975  {
976  std::ostringstream error_stream;
977  error_stream << "Empty default reestablish disributed info method "
978  << "called.\n";
979  error_stream << "This should be overloaded in a specific "
980  << "RefineableTriangleMesh\n";
981  throw OomphLibError(error_stream.str(),
982  OOMPH_CURRENT_FUNCTION,
983  OOMPH_EXCEPTION_LOCATION);
984  }
985 
986 #endif // #ifdef OOMPH_HAS_MPI
987 
988  /// \short Completely regenerate the mesh from the trianglateio structure
990  {
991  //Remove all the boundary node information
992  this->remove_boundary_nodes();
993 
994  //Delete exisiting nodes
995  unsigned n_node = this->nnode();
996  for(unsigned n=n_node;n>0;--n)
997  {
998  delete this->Node_pt[n-1];
999  this->Node_pt[n-1] = 0;
1000  }
1001  //Delete exisiting elements
1002  unsigned n_element = this->nelement();
1003  for(unsigned e=n_element;e>0;--e)
1004  {
1005  delete this->Element_pt[e-1];
1006  this->Element_pt[e-1] = 0;
1007  }
1008  //Flush the storage
1010 
1011  //Delete all boundary element information
1012  //ALH: Kick up the object hierarchy?
1013  this->Boundary_element_pt.clear();
1014  this->Face_index_at_boundary.clear();
1015  this->Region_element_pt.clear();
1016  this->Region_attribute.clear();
1017  this->Boundary_region_element_pt.clear();
1018  this->Face_index_region_at_boundary.clear();
1019  this->Boundary_curve_section_pt.clear();
1020  this->Polygonal_vertex_arclength_info.clear();
1021 
1022 #ifdef OOMPH_HAS_MPI
1023  // Delete Halo(ed) information in the old mesh
1024  if (this->is_mesh_distributed())
1025  {
1026  this->Halo_node_pt.clear();
1027  this->Root_halo_element_pt.clear();
1028 
1029  this->Haloed_node_pt.clear();
1030  this->Root_haloed_element_pt.clear();
1031 
1032  this->External_halo_node_pt.clear();
1033  this->External_halo_element_pt.clear();
1034 
1035  this->External_haloed_node_pt.clear();
1036  this->External_haloed_element_pt.clear();
1037  }
1038 #endif
1039 
1040  unsigned nbound=nboundary();
1041  Boundary_coordinate_exists.resize(nbound,false);
1042 
1043  //Now build the new scaffold
1045 
1046  // Triangulation has been created -- remember to wipe it!
1047  Triangulateio_exists=true;
1048 
1049  // Convert mesh from scaffold to actual mesh
1051 
1052  // Kill the scaffold
1053  delete this->Tmp_mesh_pt;
1054  this->Tmp_mesh_pt=0;
1055 
1056 #ifdef OOMPH_HAS_MPI
1057  if (!this->is_mesh_distributed())
1058  {
1059  nbound = this->nboundary(); // The original number of boundaries
1060  }
1061  else
1062  {
1063  nbound = this->initial_shared_boundary_id();
1064  // NOTE: The total number of boundaries is the number of
1065  // original bondaries plus the number of shared boundaries, but
1066  // here we only establish boundary coordinates for the original
1067  // boundaries. Once all the info. related with the distribution
1068  // has been established then the number of boundaries is reset
1069  // to the correct one (after reset the halo/haloed scheme)
1070  }
1071 #else
1072  nbound = this->nboundary(); // The original number of boundaries
1073 #endif
1074 
1075  // Setup boundary coordinates for boundaries
1076  for (unsigned b=0;b<nbound;b++)
1077  {
1078  this->template setup_boundary_coordinates<ELEMENT>(b);
1079  }
1080 
1081  // Snap nodes only if the mesh is not distributed, if the mesh is
1082  // distributed it will be called after the re-establishment of the
1083  // halo/haloed scheme, and the proper identification of the segments
1084  // in the boundary
1085  if (!this->is_mesh_distributed())
1086  {
1087  //Deform the boundary onto any geometric objects
1089  }
1090  }
1091 
1092  /// Boolean defining if Triangulateio object has been built or not
1094  {
1095  return Triangulateio_exists;
1096  }
1097 
1098 #endif
1099 
1100  /// \short Return the vector that contains the oomph-lib node number
1101  /// for all vertex nodes in the TriangulateIO representation of the mesh
1103  {
1104  return Oomph_vertex_nodes_id;
1105  }
1106 
1107  /// Timestepper used to build elements
1109 
1110  /// Boolean flag to indicate whether to use attributes or not (required
1111  /// for multidomain meshes)
1113 
1114  protected:
1115 
1116  /// Build mesh from scaffold
1117  void build_from_scaffold(TimeStepper* time_stepper_pt,
1118  const bool &use_attributes);
1119 
1120 #ifdef OOMPH_HAS_TRIANGLE_LIB
1121 
1122  /// \short Helper function to create TriangulateIO object (return in
1123  /// triangulate_io) from the .poly file
1124  void build_triangulateio(const std::string& poly_file_name,
1125  TriangulateIO& triangulate_io,
1126  bool& use_attributes);
1127 
1128  /// \short A general-purpose construction function that builds the
1129  /// mesh once the different specific constructors have assembled the
1130  /// appropriate information.
1132  Vector<TriangleMeshPolygon*> &internal_polygon_pt,
1134  &open_polylines_pt,
1135  const double &element_area,
1136  Vector<Vector<double> > &extra_holes_coordinates,
1137  std::map<unsigned, Vector<double> >
1138  &regions_coordinates,
1139  std::map<unsigned, double >
1140  &regions_areas,
1141  TimeStepper* time_stepper_pt,
1142  const bool &use_attributes,
1143  const bool &refine_boundary,
1144  const bool &refine_internal_boundary)
1145  {
1146  // Mesh can only be built with 2D Telements.
1147  MeshChecker::assert_geometric_element<TElementGeometricBase,ELEMENT>(2);
1148 
1149 #ifdef PARANOID
1150  if (element_area < 10e-14)
1151  {
1152  std::ostringstream warning_message;
1153  warning_message
1154  << "The current elements area was stated to (" << element_area
1155  << ").\nThe current precision to generate the input to triangle "
1156  << "is fixed to 14 digits\n\n";
1157  OomphLibWarning(warning_message.str(),
1158  OOMPH_CURRENT_FUNCTION,
1159  OOMPH_EXCEPTION_LOCATION);
1160  }
1161 #endif
1162 
1163  //Store the attribute flag
1164  Use_attributes = use_attributes;
1165 
1166  // Store Timestepper used to build elements
1167  Time_stepper_pt=time_stepper_pt;
1168 
1169  // Store outer polygon
1170  Outer_boundary_pt=outer_boundary_pt;
1171 
1172  // Store internal polygons by copy constructor
1173  Internal_polygon_pt=internal_polygon_pt;
1174 
1175  // Store internal polylines by copy constructor
1176  Internal_open_curve_pt = open_polylines_pt;
1177 
1178  // Store the extra holes coordinates
1179  Extra_holes_coordinates = extra_holes_coordinates;
1180 
1181  // Store the extra regions coordinates
1182  Regions_coordinates = regions_coordinates;
1183 
1184  // Create the data structures required to call the triangulate function
1185  TriangulateIO triangulate_io;
1186 
1187  // Initialize TriangulateIO structure
1189 
1190  // Convert TriangleMeshPolyLine and TriangleMeshClosedCurvePolyLine
1191  // to a triangulateio object
1193  internal_polygon_pt,
1194  open_polylines_pt,
1195  extra_holes_coordinates,
1196  regions_coordinates,
1197  regions_areas,
1198  triangulate_io);
1199 
1200  // Initialize TriangulateIO structure
1202 
1203  // Triangulation has been created -- remember to wipe it!
1204  Triangulateio_exists=true;
1205 
1206  // Input string for triangle
1207  std::stringstream input_string_stream;
1208  input_string_stream.precision(14);
1209  input_string_stream.setf(std::ios_base::fixed,
1210  std::ios_base::floatfield);
1211 
1212  // MH: Used to be:
1213  // input_string_stream<<"-pA -a" << element_area << " -q30" << std::fixed;
1214  // The repeated -a allows the specification of areas for different
1215  // regions (if any)
1216  input_string_stream <<"-pA -a -a" << element_area << " -q30" << std::fixed;
1217 
1218  // Verify if creation of new points on boundaries is allowed
1220  {input_string_stream<<" -YY";}
1221 
1222  //Suppress insertion of additional points on outer boundary
1223  if(refine_boundary==false)
1224  {
1225  input_string_stream << "-Y";
1226  //Add the extra flag to suppress additional points on interior segments
1227  if(refine_internal_boundary==false) {input_string_stream << "Y";}
1228  }
1229 
1230  // Convert the Input string in *char required by the triangulate function
1231  char triswitches[100];
1232  sprintf(triswitches,"%s",input_string_stream.str().c_str());
1233 
1234  // Build the mesh using triangulate function
1235  triangulate(triswitches, &triangulate_io, &Triangulateio, 0);
1236 
1237  // Build scaffold
1239 
1240  //If we have filled holes then we must use the attributes
1241  if(!regions_coordinates.empty())
1242  {
1243  // Convert mesh from scaffold to actual mesh
1244  build_from_scaffold(time_stepper_pt,true);
1245  //Record the attribute flag
1246  Use_attributes=true;
1247  }
1248  //Otherwise use what was asked
1249  else
1250  {
1251  // Convert mesh from scaffold to actual mesh
1252  build_from_scaffold(time_stepper_pt,use_attributes);
1253  }
1254 
1255  // Kill the scaffold
1256  delete this->Tmp_mesh_pt;
1257  this->Tmp_mesh_pt=0;
1258 
1259  // Cleanup but leave hole and regions alone since it's still used
1260  bool clear_hole_data=false;
1261  TriangleHelper::clear_triangulateio(triangulate_io,clear_hole_data);
1262  }
1263 
1264  /// Boolean defining if Triangulateio object has been built or not
1266 
1267 #endif // OOMPH_HAS_TRIANGLE_LIB
1268 
1269  /// \short Target areas for regions; defaults to 0.0 which (luckily)
1270  /// implies "no specific target area" for triangle!
1271  std::map<unsigned, double> Regions_areas;
1272 
1273  // Stores the ids of the internal boundaries
1274  std::set<unsigned> Internal_boundaries;
1275 
1276  /// Temporary scaffold mesh
1278 
1279  /// \short Vector storing oomph-lib node number
1280  /// for all vertex nodes in the TriangulateIO representation of the mesh
1282 
1283 #ifdef OOMPH_HAS_MPI
1284 
1285  public:
1286 
1287  // short The initial boundary id for shared boundaries
1289  {return Initial_shared_boundary_id;}
1290 
1291  // short The final boundary id for shared boundaries
1292  const unsigned final_shared_boundary_id()
1293  {return Final_shared_boundary_id;}
1294 
1295 
1296  protected:
1297 
1298 /* // \short Stores the initial arclength for the segments that appear when */
1299 /* // a boundary is splited among processors (data from previous mesh) */
1300 /* std::map<unsigned, Vector<double> > Boundary_segment_initial_arclength_previous_mesh; */
1301 
1302  // \short Stores the initial number of vertices for the segments that
1303  // appear when a boundary is splited among processors
1304  //std::map<unsigned, Vector<unsigned> > Boundary_segment_initial_vertices;
1305 
1306  // Get the shared boundaries ids living in the current processor
1309  {
1310 #ifdef PARANOID
1311  // Used to check if there are repeated shared boundaries
1312  std::set<unsigned> shared_boundaries_in_this_processor_set;
1313 #endif
1314  // Get the number of processors
1315  const unsigned n_proc = this->communicator_pt()->nproc();
1316  // Get the current processor
1317  const unsigned my_rank = this->communicator_pt()->my_rank();
1318  // Loop over all the processor and get the shared boundaries ids
1319  // associated with each processor
1320  for (unsigned iproc = 0; iproc < n_proc; iproc++)
1321  {
1322  // Work with other processors only
1323  if (iproc != my_rank)
1324  {
1325  // Get the number of boundaries shared with the "iproc"-th
1326  // processor
1327  const unsigned nshared_boundaries_with_iproc =
1328  this->nshared_boundaries(my_rank, iproc);
1329 
1330  // If there are shared boundaries associated with the current
1331  // processor then add them
1332  if (nshared_boundaries_with_iproc > 0)
1333  {
1334  // Get the boundaries ids shared with "iproc"-th processor
1335  Vector<unsigned> bound_ids_shared_with_iproc;
1336  bound_ids_shared_with_iproc =
1337  this->shared_boundaries_ids(my_rank, iproc);
1338 
1339  // Loop over shared boundaries with "iproc"-th processor
1340  for (unsigned bs = 0; bs < nshared_boundaries_with_iproc; bs++)
1341  {
1342  const unsigned bnd_id = bound_ids_shared_with_iproc[bs];
1343 #ifdef PARANOID
1344  // Check that the current shared boundary id has not been
1345  // previously added
1346  std::set<unsigned>::iterator it =
1347  shared_boundaries_in_this_processor_set.find(bnd_id);
1348  if (it != shared_boundaries_in_this_processor_set.end())
1349  {
1350  std::stringstream error;
1351  error
1352  << "The current shared boundary (" << bnd_id << ") was\n"
1353  << "already added by other pair of processors\n."
1354  << "This means that there are repeated shared boundaries ids\n";
1355  throw OomphLibError(error.str(),
1356  OOMPH_CURRENT_FUNCTION,
1357  OOMPH_EXCEPTION_LOCATION);
1358  } // if (it != shared_boundaries_in_this_processor_set.end())
1359  shared_boundaries_in_this_processor_set.insert(bnd_id);
1360 #endif
1361  shared_boundaries_in_this_processor.push_back(bnd_id);
1362  } // for (bs < nshared_boundaries_with_iproc)
1363 
1364  } // if (nshared_boundaries_with_iproc > 0)
1365 
1366  } // if (iproc != my_rank)
1367 
1368  } // for (iproc < nproc)
1369 
1370  }
1371 
1372  // Access functions to boundaries shared with processors
1373  const unsigned nshared_boundaries(const unsigned &p,
1374  const unsigned &q) const
1375  {return Shared_boundaries_ids[p][q].size();}
1376 
1378  {return Shared_boundaries_ids;}
1379 
1381  {return Shared_boundaries_ids;}
1382 
1384  {return Shared_boundaries_ids[p];}
1385 
1387  {return Shared_boundaries_ids[p];}
1388 
1390  const unsigned &q) const
1391  {return Shared_boundaries_ids[p][q];}
1392 
1393  Vector<unsigned> &shared_boundaries_ids(const unsigned &p, const unsigned &q)
1394  {return Shared_boundaries_ids[p][q];}
1395 
1396  const unsigned shared_boundaries_ids(const unsigned &p,
1397  const unsigned &q,
1398  const unsigned &i) const
1399  {return Shared_boundaries_ids[p][q][i];}
1400 
1401  const unsigned nshared_boundary_curves(const unsigned &p) const
1402  {return Shared_boundary_polyline_pt[p].size();}
1403 
1404  const unsigned nshared_boundary_polyline(const unsigned &p,
1405  const unsigned &c) const
1406  {return Shared_boundary_polyline_pt[p][c].size();}
1407 
1409  const unsigned &c)
1410  {return Shared_boundary_polyline_pt[p][c];}
1411 
1413  const unsigned &c,
1414  const unsigned &i) const
1415  {return Shared_boundary_polyline_pt[p][c][i];}
1416 
1417  const unsigned nshared_boundaries() const
1418  {return Shared_boundary_element_pt.size();}
1419 
1420  const unsigned nshared_boundary_element(const unsigned &b)
1421  {
1422  // First check if the boundary exist
1423  std::map<unsigned, Vector<FiniteElement*> >::iterator it =
1425  if (it != Shared_boundary_element_pt.end())
1426  {
1427  return Shared_boundary_element_pt[b].size();
1428  }
1429  else
1430  {
1431  std::ostringstream error_stream;
1432  error_stream
1433  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1434  throw OomphLibError(error_stream.str(),
1435  OOMPH_CURRENT_FUNCTION,
1436  OOMPH_EXCEPTION_LOCATION);
1437  }
1438  }
1439 
1441  {Shared_boundary_element_pt.clear();}
1442 
1443  void flush_shared_boundary_element(const unsigned &b)
1444  {
1445  // First check if the boundary exist
1446  std::map<unsigned, Vector<FiniteElement*> >::iterator it =
1448  if (it != Shared_boundary_element_pt.end())
1449  {
1450  Shared_boundary_element_pt[b].clear();
1451  }
1452  else
1453  {
1454  std::ostringstream error_stream;
1455  error_stream
1456  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1457  throw OomphLibError(error_stream.str(),
1458  OOMPH_CURRENT_FUNCTION,
1459  OOMPH_EXCEPTION_LOCATION);
1460  }
1461  }
1462 
1463  void add_shared_boundary_element(const unsigned &b,
1464  FiniteElement* ele_pt)
1465  {Shared_boundary_element_pt[b].push_back(ele_pt);}
1466 
1468  const unsigned &e)
1469  {
1470  // First check if the boundary exist
1471  std::map<unsigned, Vector<FiniteElement*> >::iterator it =
1473  if (it != Shared_boundary_element_pt.end())
1474  {
1475  return Shared_boundary_element_pt[b][e];
1476  }
1477  else
1478  {
1479  std::ostringstream error_stream;
1480  error_stream
1481  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1482  throw OomphLibError(error_stream.str(),
1483  OOMPH_CURRENT_FUNCTION,
1484  OOMPH_EXCEPTION_LOCATION);
1485  }
1486  }
1487 
1490 
1491  void add_face_index_at_shared_boundary(const unsigned &b,
1492  const unsigned &i)
1493  {Face_index_at_shared_boundary[b].push_back(i);}
1494 
1495  int face_index_at_shared_boundary(const unsigned &b,
1496  const unsigned &e)
1497  {
1498  // First check if the boundary exist
1499  std::map<unsigned, Vector<int> >::iterator it =
1501  if (it != Face_index_at_shared_boundary.end())
1502  {
1503  return Face_index_at_shared_boundary[b][e];
1504  }
1505  else
1506  {
1507  std::ostringstream error_stream;
1508  error_stream
1509  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1510  throw OomphLibError(error_stream.str(),
1511  OOMPH_CURRENT_FUNCTION,
1512  OOMPH_EXCEPTION_LOCATION);
1513  }
1514  }
1515 
1516  const unsigned nshared_boundary_node(const unsigned &b)
1517  {
1518  // First check if the boundary exist
1519  std::map<unsigned, Vector<Node*> >::iterator it =
1520  Shared_boundary_node_pt.find(b);
1521  if (it != Shared_boundary_node_pt.end())
1522  {
1523  return Shared_boundary_node_pt[b].size();
1524  }
1525  else
1526  {
1527  std::ostringstream error_stream;
1528  error_stream
1529  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1530  throw OomphLibError(error_stream.str(),
1531  OOMPH_CURRENT_FUNCTION,
1532  OOMPH_EXCEPTION_LOCATION);
1533  }
1534  }
1535 
1536  // Flush ALL the shared boundary nodes
1538  {Shared_boundary_node_pt.clear();}
1539 
1540  // Flush the boundary nodes associated to the shared boundary b
1541  void flush_shared_boundary_node(const unsigned &b)
1542  {Shared_boundary_node_pt[b].clear();}
1543 
1544  // Add the node the shared boundary
1545  void add_shared_boundary_node(const unsigned &b, Node* node_pt)
1546  {
1547  //Get the size of the Shared_boundary_node_pt vector
1548  const unsigned nbound_node=Shared_boundary_node_pt[b].size();
1549  bool node_already_on_this_boundary=false;
1550  //Loop over the vector
1551  for (unsigned n=0; n<nbound_node; n++)
1552  {
1553  // is the current node here already?
1554  if (node_pt==Shared_boundary_node_pt[b][n])
1555  {
1556  node_already_on_this_boundary=true;
1557  }
1558  }
1559 
1560  //Add the base node pointer to the vector if it's not there already
1561  if (!node_already_on_this_boundary)
1562  {
1563  Shared_boundary_node_pt[b].push_back(node_pt);
1564  }
1565 
1566  }
1567 
1568  Node* shared_boundary_node_pt(const unsigned &b, const unsigned &n)
1569  {
1570  // First check if the boundary exist
1571  std::map<unsigned, Vector<Node*> >::iterator it =
1572  Shared_boundary_node_pt.find(b);
1573  if (it != Shared_boundary_node_pt.end())
1574  {
1575  return Shared_boundary_node_pt[b][n];
1576  }
1577  else
1578  {
1579  std::ostringstream error_stream;
1580  error_stream
1581  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1582  throw OomphLibError(error_stream.str(),
1583  OOMPH_CURRENT_FUNCTION,
1584  OOMPH_EXCEPTION_LOCATION);
1585  }
1586  }
1587 
1588  // Is the node on the shared boundary
1589  bool is_node_on_shared_boundary(const unsigned &b, Node* const &node_pt)
1590  {
1591  // First check if the boundary exist
1592  std::map<unsigned, Vector<Node*> >::iterator it =
1593  Shared_boundary_node_pt.find(b);
1594  if (it != Shared_boundary_node_pt.end())
1595  {
1596  // Now check if the node lives on the shared boundary
1597  Vector<Node*>::iterator it_shd_nodes =
1598  std::find(Shared_boundary_node_pt[b].begin(),
1599  Shared_boundary_node_pt[b].end(),
1600  node_pt);
1601  //If the node is on this boundary
1602  if(it_shd_nodes!=Shared_boundary_node_pt[b].end())
1603  {return true;}
1604  else // The node is not on the boundary
1605  {return false;}
1606  }
1607  else
1608  {
1609  std::ostringstream error_stream;
1610  error_stream
1611  << "The shared boundary ("<< b << ") does not exist!!!\n\n";
1612  throw OomphLibError(error_stream.str(),
1613  OOMPH_CURRENT_FUNCTION,
1614  OOMPH_EXCEPTION_LOCATION);
1615  }
1616  }
1617 
1618  // Return the association of the shared boundaries with the processors
1619  std::map<unsigned, Vector<unsigned> > &shared_boundary_from_processors()
1621 
1623  {
1624  std::map<unsigned, Vector<unsigned> >::iterator it =
1626 #ifdef PARANOID
1627  if (it == Shared_boundary_from_processors.end())
1628  {
1629  std::ostringstream error_message;
1630  error_message
1631  <<"The boundary ("<<b<<") seems not to be shared by any processors,\n"
1632  <<"it is possible that the boundary was created by the user an not\n"
1633  <<"automatically by the common interfaces between processors-domains\n";
1634  throw OomphLibError(error_message.str(),
1635  OOMPH_CURRENT_FUNCTION,
1636  OOMPH_EXCEPTION_LOCATION);
1637  }
1638 #endif
1639  return (*it).second;
1640  }
1641 
1642  // \short Get the number of shared boundaries overlaping internal
1643  // boundaries
1645  {
1647  }
1648 
1649  // \short Checks if the shared boundary overlaps an internal boundary
1651  const unsigned &shd_bnd_id)
1652  {
1653  std::map<unsigned, unsigned>::iterator it =
1656  {
1657  return true;
1658  }
1659  return false;
1660  }
1661 
1662  // \short Gets the boundary id of the internal boundary that the
1663  // shared boundary lies on
1665  const unsigned &shd_bnd_id)
1666  {
1667  std::map<unsigned, unsigned>::iterator it =
1669 #ifdef PARANOID
1671  {
1672  std::ostringstream error_message;
1673  error_message
1674  <<"The shared boundary ("<<shd_bnd_id<<") does not lie on an internal "
1675  << "boundary!!!.\n"
1676  << "Make sure to call this method just for shared boundaries that lie "
1677  << "on an internal boundary.\n\n";
1678  throw OomphLibError(error_message.str(),
1679  OOMPH_CURRENT_FUNCTION,
1680  OOMPH_EXCEPTION_LOCATION);
1681  }
1682 #endif
1683  return (*it).second;
1684  }
1685 
1686  // \short Gets the shared boundaries ids that overlap the given
1687  // internal boundary
1689  const unsigned &internal_bnd_id,
1690  Vector<unsigned> &shd_bnd_ids)
1691  {
1692  // Clear any data in the output storage
1693  shd_bnd_ids.clear();
1694  // Loop over the map and store in the output vector the shared
1695  // boundaries ids that overlap the internal boundary
1696  std::map<unsigned, unsigned>::iterator it =
1698  for (; it !=Shared_boundary_overlaps_internal_boundary.end(); it++)
1699  {
1700  // If the second entry is the internal boundary, then add the
1701  // first entry to the output vector
1702  if ((*it).second == internal_bnd_id)
1703  {
1704  // Add the first entry
1705  shd_bnd_ids.push_back((*it).first);
1706  }
1707  } // loop over the map entries
1708 
1709 #ifdef PARANOID
1710  if (shd_bnd_ids.size() == 0)
1711  {
1712  std::ostringstream error_message;
1713  error_message
1714  <<" The internal boundary (" << internal_bnd_id << ") has no shared "
1715  << "boundaries overlapping it\n"
1716  << "Make sure to call this method just for internal boundaries that "
1717  << "are marked to as being\noverlaped by shared boundaries\n";
1718  throw OomphLibError(error_message.str(),
1719  OOMPH_CURRENT_FUNCTION,
1720  OOMPH_EXCEPTION_LOCATION);
1721  }
1722 #endif
1723 
1724  }
1725 
1726  // \short Gets the storage that indicates if a shared boundary is part
1727  // of an internal boundary
1728  std::map<unsigned, unsigned> &shared_boundary_overlaps_internal_boundary()
1730 
1731  // \short Helper function to verify if a given boundary was splitted
1732  // in the distribution process
1733  const bool boundary_was_splitted(const unsigned &b)
1734  {
1735  std::map<unsigned, bool>::iterator it;
1736  it = Boundary_was_splitted.find(b);
1737  if (it == Boundary_was_splitted.end())
1738  {
1739  return false;
1740  }
1741  else
1742  {
1743  return (*it).second;
1744  }
1745  }
1746 
1747  // \short Gets the number of subpolylines that create the boundarya
1748  // (useful only when the boundary is marked as split)
1749  const unsigned nboundary_subpolylines(const unsigned &b)
1750  {
1751  std::map<unsigned, Vector<TriangleMeshPolyLine*> >::iterator it;
1752  it = Boundary_subpolylines.find(b);
1753 #ifdef PARANOID
1754  if (it == Boundary_subpolylines.end())
1755  {
1756  std::ostringstream error_message;
1757  error_message
1758  <<"The boundary ("<<b<<") was marked as been splitted but there\n"
1759  <<"are not registered polylines to represent the boundary.\n"
1760  <<"The new polylines were not set up when the boundary was found to\n"
1761  <<"be splitted or the polylines have been explicitly deleted before\n"
1762  <<"being used.";
1763  throw OomphLibError(error_message.str(),
1764  OOMPH_CURRENT_FUNCTION,
1765  OOMPH_EXCEPTION_LOCATION);
1766  }
1767 #endif
1768  return (*it).second.size();
1769  }
1770 
1771  // \short Gets the vector of auxiliar polylines that will represent
1772  // the given boundary (useful only when the boundaries were
1773  // splitted)
1775  {
1776  std::map<unsigned, Vector<TriangleMeshPolyLine*> >::iterator it;
1777  it = Boundary_subpolylines.find(b);
1778  if (it == Boundary_subpolylines.end())
1779  {
1780  std::ostringstream error_message;
1781  error_message
1782  <<"The boundary ("<<b<<") was marked as been splitted but there\n"
1783  <<"are not registered polylines to represent the boundary.\n"
1784  <<"The new polylines were not set up when the boundary was found to\n"
1785  <<"be splitted or the polylines have been explicitly deleted before\n"
1786  <<"being used.";
1787  throw OomphLibError(error_message.str(),
1788  OOMPH_CURRENT_FUNCTION,
1789  OOMPH_EXCEPTION_LOCATION);
1790  }
1791  return (*it).second;
1792  }
1793 
1794  // \short Returns the value that indicates if a subpolyline of a
1795  // given boundary continues been used as internal boundary or should
1796  // be changed as shared boundary
1797  const bool boundary_marked_as_shared_boundary(const unsigned &b,
1798  const unsigned &isub)
1799  {
1800  std::map<unsigned, std::vector<bool> >::iterator it;
1802  if (it == Boundary_marked_as_shared_boundary.end())
1803  {
1804  // If no info. was found for the shared boundary then it may be
1805  // a non internal boundary, so no shared boundaries are
1806  // overlaping it
1807  return false;
1808  }
1809  return (*it).second[isub];
1810  }
1811 
1812  // short The initial boundary id for shared boundaries
1814 
1815  // short The final boundary id for shared boundaries
1817 
1818  // Stores the boundaries ids created by the interaction of two processors
1819  // Shared_boundaries_ids[iproc][jproc] = Vector of shared boundaries ids
1820  // "iproc" processor shares boundaries with "jproc" processor
1822 
1823  // Stores the processors involved in the generation of a shared
1824  // boundary, in 2D two processors give rise to the creation of a
1825  // shared boundary
1826  std::map<unsigned, Vector<unsigned> >Shared_boundary_from_processors;
1827 
1828  // Stores information about those shared boundaries that lie over or
1829  // over a segment of an internal boundary (only used when using
1830  // internal boundaries in the domain)
1831  std::map<unsigned, unsigned> Shared_boundary_overlaps_internal_boundary;
1832 
1833  // Stores the polyline representation of the shared boundaries
1834  // Shared_boundary_polyline_pt[iproc][ncurve][npolyline] = polyline_pt
1836 
1838  {Shared_boundary_polyline_pt.clear();}
1839 
1840  // Stores the boundary elements adjacent to the shared boundaries, these
1841  // elements are a subset of the halo and haloed elements
1842  std::map<unsigned, Vector<FiniteElement*> > Shared_boundary_element_pt;
1843 
1844  /// \short For the e-th finite element on shared boundary b, this is
1845  /// the index of the face that lies along that boundary
1846  std::map<unsigned, Vector<int> > Face_index_at_shared_boundary;
1847 
1848  // Stores the boundary nodes adjacent to the shared boundaries, these
1849  // nodes are a subset of the halo and haloed nodes
1850  std::map<unsigned, Vector<Node*> > Shared_boundary_node_pt;
1851 
1852  // Flag to indicate if a polyline has been splitted during the distribution
1853  // process, the boundary id of the polyline is used to indicate if spplited
1854  std::map<unsigned, bool> Boundary_was_splitted;
1855 
1856  // The polylines that will temporary represent the boundary that was
1857  // splitted in the distribution process. Used to ease the sending of
1858  // info. to Triangle during the adaptation process.
1859  std::map<unsigned, Vector<TriangleMeshPolyLine*> > Boundary_subpolylines;
1860 
1861  // Flag to indicate if an internal boundary will be used as shared boundary
1862  // because there is overlapping of the internal boundary with the shared
1863  // boundary
1864  std::map<unsigned, std::vector<bool> > Boundary_marked_as_shared_boundary;
1865 
1866  // \short Creates the distributed domain representation. Joins the
1867  // original boundaires, shared boundaries and creates connections among
1868  // them to create the new polygons that represent the distributed
1869  // domain
1871  &polygons_pt,
1873  &open_curves_pt);
1874 
1875  // \short Sorts the polylines so they be continuous and then we can
1876  // create a closed or open curve from them
1878  &unsorted_polylines_pt,
1880  &sorted_polylines_pt);
1881 
1882  // \short Take the polylines from the shared boundaries and create
1883  // temporary polygon representations of the domain
1885  &polylines_pt,
1886  Vector<TriangleMeshPolygon *> &polygons_pt);
1887 
1888  //\short Take the polylines from the original open curves and created
1889  //new temporaly representations of open curves with the bits of
1890  //original curves not overlapped by shared boundaries
1892  &sorted_open_curves_pt,
1894  &unsorted_shared_to_internal_poly_pt,
1896  &open_curves_pt);
1897 
1898  // Flag to know if it is the first time we are going to compute the
1899  // holes left by the halo elements
1901 
1902  // Backup the original extra holes coordinates
1904 
1905  // \short Compute the holes left by the halo elements, those
1906  // adjacent to the shared boundaries
1908  Vector<Vector<double> > &output_holes_coordinates);
1909 
1910  // \short Keeps those vertices that define a hole, those that are
1911  // inside closed internal boundaries in the new polygons that define the
1912  // domain. Delete those outside/inside the outer polygons (this is
1913  // required since Triangle can not deal with vertices that define
1914  // holes outside the new outer polygons of the domain)
1916  Vector<TriangleMeshPolygon *> &polygons_pt,
1917  Vector<Vector<double> > &output_holes_coordinates);
1918 
1919  // \short Check for any possible connections that the array of
1920  // sorted nodes have with any previous boundaries or with
1921  // itself. Return -1 if no connection was found, return -2 if the
1922  // connection is with the same polyline, return the boundary id of
1923  // the boundary to which the connection is performed
1925  std::set<FiniteElement*> &element_in_processor_pt,
1926  const int &root_edge_bnd_id,
1927  std::map<std::pair<Node*,Node*>, bool> &overlapped_face,
1928  std::map<unsigned, std::map<Node*, bool> >
1929  &node_on_bnd_not_overlapped_by_shd_bnd,
1930  std::list<Node*> &current_polyline_nodes,
1931  std::map<unsigned, std::list<Node*> >
1932  &shared_bnd_id_to_sorted_list_node_pt,
1933  const unsigned &node_degree,
1934  Node* &new_node_pt,
1935  const bool called_from_load_balance=false);
1936 
1937  // \short Establish the connections of the polylines previously marked
1938  // as having connections. This connections were marked in the function
1939  // TriangleMesh::create_polylines_from_halo_elements_helper().
1941 
1942  /// \short Creates the shared boundaries
1944  const Vector<unsigned> &element_domain,
1946  &backed_up_el_pt,
1948  &backed_up_f_el_pt,
1949  std::map<Data*,std::set<unsigned> >
1950  &processors_associated_with_data,
1951  const bool&
1952  overrule_keep_as_halo_element_status);
1953 
1954  // \short Creates the halo elements on all processors
1955  // Gets the halo elements on all processors, these elements are then used
1956  // on the function that computes the shared boundaries among the processors
1957  void get_halo_elements_on_all_procs(const unsigned &nproc,
1958  const Vector<unsigned>
1959  &element_domain,
1961  &backed_up_el_pt,
1962  std::map<Data*,std::set<unsigned> >
1963  &processors_associated_with_data,
1964  const bool&
1965  overrule_keep_as_halo_element_status,
1966  std::map<GeneralisedElement*, unsigned>
1967  &element_to_global_index,
1968  Vector<Vector<
1970  output_halo_elements_pt);
1971 
1972  // \short Get the element edges (pair of nodes, edges) that lie
1973  // on a boundary (used to mark shared boundaries that lie on
1974  // internal boundaries)
1975  void get_element_edges_on_boundary(std::map<std::pair<Node*,Node*>,
1976  unsigned> &element_edges_on_boundary);
1977 
1978  // \short Creates polylines from the intersection of halo elements
1979  // on all processors. The new polylines define the shared boundaries
1980  // in the domain This get the polylines on ALL processors, that is
1981  // why the three dimensions
1982  // output_polylines_pt[iproc][ncurve][npolyline]
1984  const Vector<unsigned> &element_domain,
1985  std::map<GeneralisedElement*, unsigned> &element_to_global_index,
1986  std::set<FiniteElement*> &element_in_processor_pt,
1987  Vector<Vector<Vector<GeneralisedElement*> > > &input_halo_elements,
1988  std::map<std::pair<Node*,Node*>, unsigned> &elements_edges_on_boundary,
1989  Vector<Vector<Vector<TriangleMeshPolyLine *> > > &output_polylines_pt);
1990 
1991  // \short Break any possible loop created by the sorted list of
1992  // nodes that is used to create a new shared polyline
1994  const unsigned &initial_shd_bnd_id,
1995  std::list<Node*> &input_nodes,
1996  Vector<FiniteElement*> &input_boundary_element_pt,
1997  Vector<int> &input_face_index_element,
1998  const int &input_connect_to_the_left,
1999  const int &input_connect_to_the_right,
2000  Vector<std::list<Node*> > &output_sorted_nodes_pt,
2001  Vector<Vector<FiniteElement*> > &output_boundary_element_pt,
2002  Vector<Vector<int> > &output_face_index_element,
2003  Vector<int> &output_connect_to_the_left,
2004  Vector<int> &output_connect_to_the_right);
2005 
2006  // \short Break any possible loop created by the sorted list of
2007  // nodes that is used to create a new shared polyline (modified
2008  // version for load balance)
2010  const unsigned &initial_shd_bnd_id,
2011  std::list<Node*> &input_nodes,
2012  Vector<FiniteElement*> &input_boundary_element_pt,
2013  Vector<FiniteElement*> &input_boundary_face_element_pt,
2014  Vector<int> &input_face_index_element,
2015  const int &input_connect_to_the_left,
2016  const int &input_connect_to_the_right,
2017  Vector<std::list<Node*> > &output_sorted_nodes_pt,
2018  Vector<Vector<FiniteElement*> > &output_boundary_element_pt,
2019  Vector<Vector<FiniteElement*> > &output_boundary_face_element_pt,
2020  Vector<Vector<int> > &output_face_index_element,
2021  Vector<int> &output_connect_to_the_left,
2022  Vector<int> &output_connect_to_the_right);
2023 
2024  // \short Create the shared polyline and fill the data structured
2025  // that keep all the information associated with the creationg of the
2026  // shared boundary
2027  void create_shared_polyline(const unsigned &my_rank,
2028  const unsigned &shd_bnd_id,
2029  const unsigned &iproc,
2030  const unsigned &jproc,
2031  std::list<Node*> &sorted_nodes,
2032  const int &root_edge_bnd_id,
2034  &bulk_bnd_ele_pt,
2035  Vector<int> &face_index_ele,
2037  &unsorted_polylines_pt,
2038  const int &connect_to_the_left_flag,
2039  const int &connect_to_the_right_flag);
2040 
2041  public:
2042 
2043  /// Virtual function to perform the load balance routines
2044  virtual void load_balance(const Vector<unsigned>&
2045  target_domain_for_local_non_halo_element)
2046  {
2047  std::ostringstream error_stream;
2048  error_stream << "Empty default load balancing function called.\n";
2049  error_stream << "This should be overloaded in a specific "
2050  << "RefineableTriangleMesh\n";
2051  throw OomphLibError(error_stream.str(),
2052  OOMPH_CURRENT_FUNCTION,
2053  OOMPH_EXCEPTION_LOCATION);
2054  }
2055 
2056  /// Virtual function to perform the reset boundary elements info
2057  /// routines. Generally used after load balance.
2058  virtual void reset_boundary_element_info(
2059  Vector<unsigned> &ntmp_boundary_elements,
2060  Vector<Vector<unsigned> > &ntmp_boundary_elements_in_region,
2061  Vector<FiniteElement*> &deleted_elements);
2062 
2063 #endif // #ifdef OOMPH_HAS_MPI
2064 
2065 
2066  public:
2067 
2068  /// Output the nodes on the boundary and their respective boundary
2069  /// coordinates(into separate tecplot zones)
2070  void output_boundary_coordinates(const unsigned &b,
2071  std::ostream &outfile);
2072 
2073 
2074  private:
2075 
2076  // Reference :
2077  // http://en.wikibooks.org/wiki/Algorithm_Implementation/Geometry/Convex_hull/Monotone_chain#C.2B.2B
2078 
2079  // Monotone chain 2D convex hull algorithm.
2080  // Asymptotic complexity: O(n log n).
2081  // Practical performance: 0.5-1.0 seconds for n=1000000 on a 1GHz machine.
2082  typedef double coord_t; // coordinate type
2083  typedef double coord2_t; // must be big enough to hold 2*max(|coordinate|)^2
2084 
2085  struct Point {
2087 
2088  bool operator <(const Point &p) const {
2089  return x < p.x || (x == p.x && y < p.y);
2090  }
2091  };
2092 
2093  // 2D cross product of OA and OB vectors, i.e. z-component of their 3D
2094  // cross product.
2095  // Returns a positive value, if OAB makes a counter-clockwise turn,
2096  // negative for clockwise turn, and zero if the points are collinear.
2097  coord2_t cross(const Point &O, const Point &A, const Point &B)
2098  {
2099  return (A.x - O.x) * (B.y - O.y) - (A.y - O.y) * (B.x - O.x);
2100  }
2101 
2102  // Returns a list of points on the convex hull in counter-clockwise order.
2103  // Note: the last point in the returned list is the same as the first one.
2104  std::vector<Point> convex_hull(std::vector<Point> P)
2105  {
2106  int n = P.size(), k = 0;
2107  std::vector<Point> H(2*n);
2108 
2109  // Sort points lexicographically
2110  std::sort(P.begin(), P.end());
2111 
2112  // Build lower hull
2113  for (int i = 0; i < n; ++i) {
2114  while (k >= 2 && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
2115  H[k++] = P[i];
2116  }
2117 
2118  // Build upper hull
2119  for (int i = n-2, t = k+1; i >= 0; i--) {
2120  while (k >= t && cross(H[k-2], H[k-1], P[i]) <= 0) k--;
2121  H[k++] = P[i];
2122  }
2123 
2124  H.resize(k);
2125  return H;
2126  }
2127 
2128  };
2129 
2130 
2131 //////////////////////////////////////////////////////////////////////
2132 //////////////////////////////////////////////////////////////////////
2133 //////////////////////////////////////////////////////////////////////
2134 
2135 
2136 //=========================================================================
2137 /// Unstructured refineable Triangle Mesh
2138 //=========================================================================
2139 // Temporary flag to enable full annotation of RefineableTriangleMesh
2140 // comms
2141 //#define ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
2142  template<class ELEMENT>
2144  public virtual RefineableMeshBase
2145  {
2146 
2147  public:
2148 
2149  /// \short Function pointer to function that updates the
2150  /// mesh following the snapping of boundary nodes to the
2151  /// boundaries (e.g. to move boundary nodes very slightly
2152  /// to satisfy volume constraints)
2153  typedef void (*MeshUpdateFctPt)(Mesh* mesh_pt);
2154 
2155  /// \short Function pointer to a function that can generate
2156  /// a point within the ihole-th hole, so that this can be
2157  /// overloaded by the user if they have a better way of
2158  /// doing it than our clunky default. The function
2159  /// should update the components of the
2160  /// Vector poly_pt->internal_point()
2161  typedef void (*InternalHolePointUpdateFctPt)(
2162  const unsigned &ihole, TriangleMeshPolygon* poly_pt);
2163 
2164 #ifdef OOMPH_HAS_TRIANGLE_LIB
2165 
2166  /// \short Build mesh, based on the specifications on
2167  /// TriangleMeshParameters
2169  TimeStepper* time_stepper_pt=
2171  : TriangleMesh<ELEMENT>(triangle_mesh_parameters,
2172  time_stepper_pt)
2173  {
2174  // Initialise the data associated with adaptation
2175  initialise_adaptation_data();
2176 
2177  // Initialise the data associated with boundary refinements
2178  initialise_boundary_refinement_data();
2179  }
2180 
2181 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
2182 
2183  /// \short Build mesh, based on the polyfiles
2184  RefineableTriangleMesh(const std::string& node_file_name,
2185  const std::string& element_file_name,
2186  const std::string& poly_file_name,
2187  TimeStepper* time_stepper_pt=
2189  const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
2190  : TriangleMesh<ELEMENT>(node_file_name,
2191  element_file_name,
2192  poly_file_name,
2193  time_stepper_pt,
2194  allow_automatic_creation_of_vertices_on_boundaries)
2195  {
2196  // Create and fill the data structures to give rise to polylines so that
2197  // the mesh can use the adapt methods
2198  create_polylines_from_polyfiles(node_file_name, poly_file_name);
2199 
2200  // Initialise the data associated with adaptation
2201  initialise_adaptation_data();
2202 
2203  // Initialise the data associated with boundary refinements
2204  initialise_boundary_refinement_data();
2205  }
2206 
2207  protected:
2208 
2209 #ifdef OOMPH_HAS_TRIANGLE_LIB
2210 
2211  /// \short Build mesh from specified triangulation and
2212  /// associated target areas for elements in it
2213  /// NOTE: This is used ONLY during adaptation and should not be used
2214  /// as a method of constructing a TriangleMesh object in demo drivers!
2216  TriangulateIO& triangulate_io,
2217  TimeStepper* time_stepper_pt=
2219  const bool &use_attributes=false,
2220  const bool
2221  &allow_automatic_creation_of_vertices_on_boundaries=true,
2222  OomphCommunicator* comm_pt = 0)
2223  {
2224  // Initialise the data associated with adaptation
2225  initialise_adaptation_data();
2226 
2227  // Initialise the data associated with boundary refinements
2228  initialise_boundary_refinement_data();
2229 
2230  // Store Timestepper used to build elements
2231  this->Time_stepper_pt=time_stepper_pt;
2232 
2233  // Create triangulateio object to refine
2234  TriangulateIO triangle_refine;
2235 
2236  // Initialize triangulateio structure
2237  TriangleHelper::initialise_triangulateio(this->Triangulateio);
2238 
2239  // Triangulation has been created -- remember to wipe it!
2240  this->Triangulateio_exists=true;
2241 
2242  // Create refined TriangulateIO structure based on target areas
2243  this->refine_triangulateio(triangulate_io,
2244  target_area,
2245  triangle_refine);
2246 
2247  // Input string for triangle
2248  std::stringstream input_string_stream;
2249  input_string_stream<<"-pq30-ra";
2250 
2251  // Verify if creation of new points on boundaries is allowed
2252  if (!allow_automatic_creation_of_vertices_on_boundaries)
2253  {input_string_stream<<" -YY";}
2254 
2255  // Copy the allowing of creation of points on the boundaries status
2256  this->Allow_automatic_creation_of_vertices_on_boundaries =
2257  allow_automatic_creation_of_vertices_on_boundaries;
2258 
2259  //Store the attribute flag
2260  this->Use_attributes=use_attributes;
2261 
2262  // Convert to a *char required by the triangulate function
2263  char triswitches[100];
2264  sprintf(triswitches,"%s",input_string_stream.str().c_str());
2265 
2266  // Build triangulateio refined object
2267  triangulate(triswitches,&triangle_refine,&this->Triangulateio,0);
2268 
2269  // Build scaffold
2270  this->Tmp_mesh_pt=new TriangleScaffoldMesh(this->Triangulateio);
2271 
2272  // Convert mesh from scaffold to actual mesh
2273  this->build_from_scaffold(time_stepper_pt,use_attributes);
2274 
2275  // Kill the scaffold
2276  delete this->Tmp_mesh_pt;
2277  this->Tmp_mesh_pt=0;
2278 
2279  // Cleanup but leave hole alone as it's still used
2280  bool clear_hole_data=false;
2281  TriangleHelper::clear_triangulateio(triangle_refine,clear_hole_data);
2282 
2283 #ifdef OOMPH_HAS_MPI
2284  // Before calling setup boundary coordinates check if the mesh
2285  // should be treated as a distributed mesh
2286  if (comm_pt!=0)
2287  {
2288  // Set the communicator which will mark the mesh as distributed
2289  this->set_communicator_pt(comm_pt);
2290  }
2291 #endif
2292 
2293  // Setup boundary coordinates for boundaries
2294  unsigned nb=nboundary();
2295  for (unsigned b=0;b<nb;b++)
2296  {
2297  this->template setup_boundary_coordinates<ELEMENT>(b);
2298  }
2299  }
2300 
2301  public:
2302 
2303 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
2304 
2305  /// Empty Destructor
2307 
2308  /// \short Enables info. and timings for tranferring of target
2309  /// areas
2311  {Print_timings_transfering_target_areas=true;}
2312 
2313  /// \short Disables info. and timings for tranferring of target
2314  /// areas
2316  {Print_timings_transfering_target_areas=false;}
2317 
2318  /// \short Enables the solution projection step during adaptation
2320  {Disable_projection=false;}
2321 
2322  /// \short Disables the solution projection step during adaptation
2324  {Disable_projection=true;}
2325 
2326  /// \short Enables info. and timings for projection
2328  {Print_timings_projection=true;}
2329 
2330  /// \short Disables info. and timings for projection
2332  {Print_timings_projection=false;}
2333 
2334  /// \short Read/write access to number of bins in the x-direction
2335  /// when transferring target areas by bin method
2336  unsigned& nbin_x_for_area_transfer(){return Nbin_x_for_area_transfer;}
2337 
2338  /// \short Read/write access to number of bins in the y-direction
2339  /// when transferring target areas by bin method
2340  unsigned& nbin_y_for_area_transfer(){return Nbin_y_for_area_transfer;}
2341 
2342  /// \short Read/write access to number of bins in the x-direction
2343  /// when projecting old solution onto new mesh
2344  unsigned& nbin_x_for_projection(){return Nbin_x_for_projection;}
2345 
2346  /// \short Read/write access to number of bins in the y-direction
2347  /// when projecting old solution onto new mesh
2348  unsigned& nbin_y_for_projection(){return Nbin_y_for_projection;}
2349 
2350  /// Max element size allowed during adaptation
2351  double& max_element_size(){return Max_element_size;}
2352 
2353  /// Min element size allowed during adaptation
2354  double& min_element_size(){return Min_element_size;}
2355 
2356  /// Min angle before remesh gets triggered
2357  double& min_permitted_angle(){return Min_permitted_angle;}
2358 
2359  // Returns the status of using an iterative solver for the
2360  // projection problem
2362  {return Use_iterative_solver_for_projection;}
2363 
2364  /// Enables the use of an iterative solver for the projection
2365  /// problem
2367  {Use_iterative_solver_for_projection=true;}
2368 
2369  /// Enables the use of an iterative solver for the projection
2370  /// problem
2372  {Use_iterative_solver_for_projection=false;}
2373 
2374  /// Enables printing of timings for adaptation
2375  void enable_print_timings_adaptation(const unsigned &print_level=1)
2376  {set_print_level_timings_adaptation(print_level);}
2377 
2378  /// Disables printing of timings for adaptation
2380  {Print_timings_level_adaptation=0;}
2381 
2382  /// Sets the printing level of timings for adaptation
2383  void set_print_level_timings_adaptation(const unsigned &print_level)
2384  {
2385  const unsigned max_print_level = 3;
2386  // If printing level is greater than max. printing level
2387  if (print_level > max_print_level)
2388  {
2389  Print_timings_level_adaptation=max_print_level;
2390  }
2391  else
2392  {
2393  Print_timings_level_adaptation=print_level;
2394  }
2395  }
2396 
2397  /// Enables printing of timings for load balance
2398  void enable_print_timings_load_balance(const unsigned &print_level=1)
2399  {set_print_level_timings_load_balance(print_level);}
2400 
2401  /// Disables printing of timings for load balance
2403  {Print_timings_level_load_balance=0;}
2404 
2405  /// Sets the printing level of timings for load balance
2406  void set_print_level_timings_load_balance(const unsigned &print_level)
2407  {
2408  const unsigned max_print_level = 3;
2409  // If printing level is greater than max. printing level
2410  if (print_level > max_print_level)
2411  {
2412  Print_timings_level_load_balance=max_print_level;
2413  }
2414  else
2415  {
2416  Print_timings_level_load_balance=print_level;
2417  }
2418  }
2419 
2420  /// Doc the targets for mesh adaptation
2421  void doc_adaptivity_targets(std::ostream &outfile)
2422  {
2423  outfile << std::endl;
2424  outfile << "Targets for mesh adaptation: " << std::endl;
2425  outfile << "---------------------------- " << std::endl;
2426  outfile << "Target for max. error: " << Max_permitted_error << std::endl;
2427  outfile << "Target for min. error: " << Min_permitted_error << std::endl;
2428  outfile << "Target min angle: " << Min_permitted_angle << std::endl;
2429  outfile << "Min. allowed element size: " << Min_element_size << std::endl;
2430  outfile << "Max. allowed element size: " << Max_element_size << std::endl;
2431  outfile << "Don't unrefine if less than " << Max_keep_unrefined
2432  << " elements need unrefinement." << std::endl;
2433  outfile << std::endl;
2434  }
2435 
2436  /// Refine mesh uniformly and doc process
2437  void refine_uniformly(DocInfo& doc_info)
2438  {
2439  // Set the element error to something big
2440  unsigned nelem=nelement();
2441  Vector<double> elem_error(nelem,DBL_MAX);
2442 
2443  // Limit the min element size to 1/3 of the current minimum
2444  double backup=Min_element_size;
2445 
2446  // Get current max and min element size
2447  double orig_max_area, orig_min_area;
2448  this->max_and_min_element_size(orig_max_area, orig_min_area);
2449 
2450  // Limit
2451  Min_element_size=orig_min_area/3.0;
2452 
2453  // Do it...
2454  adapt(elem_error);
2455 
2456  // Reset
2457  Min_element_size=backup;
2458 
2459 /* throw OomphLibError("refine_uniformly() not implemented yet", */
2460 /* OOMPH_CURRENT_FUNCTION, */
2461 /* OOMPH_EXCEPTION_LOCATION); */
2462  }
2463 
2464  /// \short Unrefine mesh uniformly: Return 0 for success,
2465  /// 1 for failure (if unrefinement has reached the coarsest permitted
2466  /// level)
2468  {
2469  throw OomphLibError("unrefine_uniformly() not implemented yet",
2470  OOMPH_CURRENT_FUNCTION,
2471  OOMPH_EXCEPTION_LOCATION);
2472  // dummy return
2473  return 0;
2474  }
2475 
2476  /// Adapt mesh, based on elemental error provided
2477  void adapt(const Vector<double>& elem_error);
2478 
2479  /// \short Access to function pointer to function that updates the
2480  /// mesh following the snapping of boundary nodes to the
2481  /// boundaries (e.g. to move boundary nodes very slightly
2482  /// to satisfy volume constraints)
2483  MeshUpdateFctPt& mesh_update_fct_pt()
2484  {
2485  return Mesh_update_fct_pt;
2486  }
2487 
2488  /// \short Access to function pointer to can be
2489  /// used to generate the internal point for the ihole-th
2490  ///hole.
2491  InternalHolePointUpdateFctPt& internal_hole_point_update_fct_pt()
2492  {
2493  return Internal_hole_point_update_fct_pt;
2494  }
2495 
2496 
2497 
2498 
2499 #ifdef OOMPH_HAS_MPI
2500  unsigned nsorted_shared_boundary_node(unsigned &b)
2501  {
2502  std::map<unsigned, Vector<Node*> >::iterator it =
2503  Sorted_shared_boundary_node_pt.find(b);
2504  if (it == Sorted_shared_boundary_node_pt.end())
2505  {
2506  std::ostringstream error_message;
2507  error_message
2508  <<"The boundary ("<<b<<") is not marked as shared\n";
2509  throw OomphLibError(error_message.str(),
2510  OOMPH_CURRENT_FUNCTION,
2511  OOMPH_EXCEPTION_LOCATION);
2512  }
2513  return (*it).second.size();
2514  }
2515 
2517  {Sorted_shared_boundary_node_pt.clear();}
2518 
2519  Node* sorted_shared_boundary_node_pt(unsigned &b, unsigned &i)
2520  {
2521  std::map<unsigned, Vector<Node*> >::iterator it =
2522  Sorted_shared_boundary_node_pt.find(b);
2523  if (it == Sorted_shared_boundary_node_pt.end())
2524  {
2525  std::ostringstream error_message;
2526  error_message
2527  <<"The boundary ("<<b<<") is not marked as shared\n";
2528  throw OomphLibError(error_message.str(),
2529  OOMPH_CURRENT_FUNCTION,
2530  OOMPH_EXCEPTION_LOCATION);
2531  }
2532  return (*it).second[i];
2533  }
2534 
2535 
2537  {
2538  std::map<unsigned, Vector<Node*> >::iterator it =
2539  Sorted_shared_boundary_node_pt.find(b);
2540  if (it == Sorted_shared_boundary_node_pt.end())
2541  {
2542  std::ostringstream error_message;
2543  error_message
2544  <<"The boundary ("<<b<<") is not marked as shared\n";
2545  throw OomphLibError(error_message.str(),
2546  OOMPH_CURRENT_FUNCTION,
2547  OOMPH_EXCEPTION_LOCATION);
2548  }
2549  return (*it).second;
2550  }
2551 
2552 #endif // #ifdef OOMPH_HAS_MPI
2553 
2554  /// Helper function to create polylines and fill associate data
2555  // structures, used when creating from a mesh from polyfiles
2556  void create_polylines_from_polyfiles(const std::string& node_file_name,
2557  const std::string& poly_file_name);
2558 
2559 #ifdef OOMPH_HAS_MPI
2560  // \short Fill the boundary elements structures when dealing with
2561  // shared boundaries that overlap internal boundaries
2562  void fill_boundary_elements_and_nodes_for_internal_boundaries();
2563 
2564  // \short Fill the boundary elements structures when dealing with
2565  // shared boundaries that overlap internal boundaries. Document the
2566  // number of elements on the shared boundaries that go to internal
2567  // boundaries
2568  void fill_boundary_elements_and_nodes_for_internal_boundaries(
2569  std::ofstream& outfile);
2570 
2571  /// Used to re-establish any additional info. related with
2572  /// the distribution after a re-starting for triangle meshes
2574  OomphCommunicator* comm_pt, std::istream &restart_file)
2575  {
2576  // Ensure that the mesh is distributed
2577  if (this->is_mesh_distributed())
2578  {
2579  // Fill the structures for the boundary elements and face indexes
2580  // of the boundary elements
2581  this->fill_boundary_elements_and_nodes_for_internal_boundaries();
2582 
2583  // Re-establish the shared boundary elements and nodes scheme
2584  // before re-establish halo and haloed elements
2585  this->reset_shared_boundary_elements_and_nodes(comm_pt);
2586 
2587  // Sort the nodes on the boundaries so that they have the same order
2588  // on all the boundaries
2589  this->sort_nodes_on_shared_boundaries();
2590 
2591  // Re-establish the halo(ed) scheme
2592  this->reset_halo_haloed_scheme();
2593 
2594  // Re-set the number of boundaries to the original one
2595  const unsigned noriginal_boundaries =
2596  this->initial_shared_boundary_id();
2597  this->set_nboundary(noriginal_boundaries);
2598 
2599  // Go through the original boudaries and re-establish the
2600  // boundary coordinates
2601  for (unsigned b = 0; b < noriginal_boundaries; b++)
2602  {
2603  // Identify the segment boundaries and re-call the
2604  // setup_boundary_coordinates() method for the new boundaries
2605  // from restart
2606  this->identify_boundary_segments_and_assign_initial_zeta_values(b, this);
2607 
2608  if (this->boundary_geom_object_pt(b) != 0)
2609  {
2610  // Re-set the boundary coordinates
2611  this->template setup_boundary_coordinates<ELEMENT>(b);
2612  }
2613  }
2614 
2615  //Deform the boundary onto any geometric objects
2616  this->snap_nodes_onto_geometric_objects();
2617 
2618  } // if (this->is_mesh_distributed())
2619 
2620  }
2621 
2622 #endif // #ifdef OOMPH_HAS_MPI
2623 
2624  /// Method used to update the polylines representation after restart
2626  {
2627 #ifdef OOMPH_HAS_MPI
2628  // If the mesh is distributed then also update the shared
2629  // boundaries
2630  unsigned my_rank = 0;
2631  if (this->is_mesh_distributed())
2632  {
2633  my_rank = this->communicator_pt()->my_rank();
2634  }
2635 #endif
2636 
2637  // Update the polyline representation after restart
2638 
2639  // First update all internal boundaries
2640  const unsigned ninternal = this->Internal_polygon_pt.size();
2641  for (unsigned i_internal = 0; i_internal < ninternal; i_internal++)
2642  {
2643  this->update_polygon_after_restart(
2644  this->Internal_polygon_pt[i_internal]);
2645  }
2646 
2647  // then update the external boundaries
2648  const unsigned nouter = this->Outer_boundary_pt.size();
2649  for (unsigned i_outer = 0; i_outer < nouter; i_outer++)
2650  {
2651  this->update_polygon_after_restart(this->Outer_boundary_pt[i_outer]);
2652  }
2653 
2654 #ifdef OOMPH_HAS_MPI
2655  // If the mesh is distributed then also update the shared
2656  // boundaries
2657  if (this->is_mesh_distributed())
2658  {
2659  const unsigned ncurves = this->nshared_boundary_curves(my_rank);
2660  for (unsigned nc = 0; nc < ncurves; nc ++)
2661  {
2662  // Update the shared polyline
2663  this->update_shared_curve_after_restart(
2664  this->Shared_boundary_polyline_pt[my_rank][nc]//shared_curve
2665  );
2666  }
2667 
2668  } // if (this->is_mesh_distributed())
2669 #endif // #ifdef OOMPH_HAS_MPI
2670 
2671  const unsigned n_open_polyline = this->Internal_open_curve_pt.size();
2672  for (unsigned i = 0; i < n_open_polyline; i++)
2673  {
2674  this->update_open_curve_after_restart(
2675  this->Internal_open_curve_pt[i]);
2676  }
2677 
2678  }
2679 
2680 #ifdef OOMPH_HAS_MPI
2681 
2682  /// \short Performs the load balancing for unstructured meshes, the
2683  /// load balancing strategy is based on mesh migration
2684  void load_balance(const Vector<unsigned>&
2685  input_target_domain_for_local_non_halo_element);
2686 
2687  /// Use the first and second group of elements to find the
2688  /// intersection between them to get the shared boundary
2689  /// elements from the first and second group
2690  void get_shared_boundary_elements_and_face_indexes(
2691  const Vector<FiniteElement*> &first_element_pt,
2692  const Vector<FiniteElement*> &second_element_pt,
2693  Vector<FiniteElement*> &first_shared_boundary_element_pt,
2694  Vector<unsigned> &first_shared_boundary_element_face_index,
2695  Vector<FiniteElement*> &second_shared_boundary_element_pt,
2696  Vector<unsigned> &second_shared_boundary_element_face_index);
2697 
2698  /// \short Creates the new shared boundaries, this method is also in
2699  /// charge of computing the shared boundaries ids of each processor
2700  /// and send that info. to all the processors
2701  void create_new_shared_boundaries(std::set<FiniteElement*>
2702  &element_in_processor_pt,
2704  &new_shared_boundary_element_pt,
2706  &new_shared_boundary_element_face_index);
2707 
2708  // Computes the degree of the nodes on the shared boundaries, the
2709  // degree of the node is computed from the global graph created by the
2710  // shared boundaries of all processors
2711  void compute_shared_node_degree_helper(Vector<Vector<FiniteElement*> >
2712  &unsorted_face_ele_pt,
2713  std::map<Node*, unsigned>
2714  &global_node_degree);
2715 
2716  // Sort the nodes on the new shared boundaries (after load balancing),
2717  // computes the alias of the nodes and creates the adjacency matrix
2718  // that represent the graph created by the shared edges between each
2719  // pair of processors
2720  void create_adjacency_matrix_new_shared_edges_helper(
2721  Vector<Vector<FiniteElement*> > &unsorted_face_ele_pt,
2722  Vector<Vector<Node*> > &tmp_sorted_shared_node_pt,
2723  std::map<Node*, Vector<Vector<unsigned> > > &node_alias,
2724  Vector<Vector<Vector<unsigned> > > &adjacency_matrix);
2725 
2726  /// \short Get the nodes on the shared boundary (b), these are stored
2727  /// in the segment they belong
2728  void get_shared_boundary_segment_nodes_helper(
2729  const unsigned &shd_bnd_id, Vector<Vector<Node*> > &tmp_segment_nodes);
2730 
2731 #endif // #ifdef OOMPH_HAS_MPI
2732 
2733  /// \short Get the nodes on the boundary (b), these are stored in
2734  /// the segment they belong (also used by the load balance method
2735  /// to re-set the number of segments per boundary after load
2736  /// balance has taken place)
2737  void get_boundary_segment_nodes_helper(
2738  const unsigned &b, Vector<Vector<Node*> > &tmp_segment_nodes);
2739 
2740  // Enable/disable unrefinement/refinement methods for original
2741  // boundaries
2743  {Do_boundary_unrefinement_constrained_by_target_areas = true;}
2744 
2746  {Do_boundary_unrefinement_constrained_by_target_areas = false;}
2747 
2749  {Do_boundary_refinement_constrained_by_target_areas = true;}
2750 
2752  {Do_boundary_refinement_constrained_by_target_areas = false;}
2753 
2754  // Enable/disable unrefinement/refinement methods for shared
2755  // boundaries
2757  {Do_shared_boundary_unrefinement_constrained_by_target_areas = true;}
2758 
2760  {Do_shared_boundary_unrefinement_constrained_by_target_areas = false;}
2761 
2763  {Do_shared_boundary_refinement_constrained_by_target_areas = true;}
2764 
2766  {Do_shared_boundary_refinement_constrained_by_target_areas = false;}
2767 
2768  // From here to the end of the class everything is protected
2769  protected:
2770 
2771  /// A map that stores the vertices that receive connections, they
2772  /// are identified by the boundary number that receive the connection
2773  /// This is necessary for not erasing them on the adaptation process,
2774  /// specifically for the un-refinement process
2775  std::map<unsigned, std::set<Vector<double> > > Boundary_connections_pt;
2776 
2777  /// \short Verifies if the given boundary receives a connection, and
2778  /// if that is the case then returns the list of vertices that
2779  /// receive the connections
2780  const bool boundary_connections(const unsigned &b,
2781  const unsigned &c,
2782  std::set<Vector<double> > &vertices)
2783  {
2784  // Search for the given boundary
2785  std::map<unsigned, std::set<Vector<double> > >::
2786  iterator it = Boundary_connections_pt.find(b);
2787  // Was the boundary found?
2788  if (it != Boundary_connections_pt.end())
2789  {
2790  // Return the set of vertices that receive the connection
2791  vertices = (*it).second;
2792  return true;
2793  }
2794  else
2795  {
2796  return false;
2797  }
2798 
2799  }
2800 
2801  /// \short Synchronise the vertices that are marked for non deletion
2802  // on the shared boundaries. Unrefinement of shared boundaries is
2803  // performed only if the candidate node is not marked for non deletion
2804  const void synchronize_shared_boundary_connections();
2805 
2806  /// \short Mark the vertices that are not allowed for deletion by
2807  /// the unrefienment/refinement polyline methods. In charge of
2808  /// filling the Boundary_chunk_connections_pt structure
2809  void add_vertices_for_non_deletion();
2810 
2811  /// \short Adds the vertices from the sources boundary that are
2812  /// repeated in the destination boundary to the list of non
2813  /// delete-able vertices in the destination boundary
2814  void add_non_delete_vertices_from_boundary_helper(
2815  Vector<Vector<Node*> > src_bound_segment_node_pt,
2816  Vector<Vector<Node*> > dst_bound_segment_node_pt,
2817  const unsigned &dst_bnd_id, const unsigned &dst_bnd_chunk);
2818 
2819  /// \short After unrefinement and refinement has taken place compute
2820  /// the new vertices numbers of the temporary representation of the
2821  // boundaries to connect.
2822  void create_temporary_boundary_connections(
2823  Vector<TriangleMeshPolygon *> &tmp_outer_polygons_pt,
2824  Vector<TriangleMeshOpenCurve *> &tmp_open_curves_pt);
2825 
2826  /// \short After unrefinement and refinement has taken place compute
2827  /// the new vertices numbers of the boundaries to connect (in a
2828  /// distributed scheme it may be possible that the destination boundary
2829  /// does no longer exist, therefore the connection is suspended and
2830  /// resumed after the adaptation processor
2831  void restore_boundary_connections(
2832  Vector<TriangleMeshPolyLine*> &resume_initial_connection_polyline_pt,
2833  Vector<TriangleMeshPolyLine*> &resume_final_connection_polyline_pt);
2834 
2835  /// \short Restore the connections of the specific polyline
2836  /// The vertices numbering on the destination boundaries may have
2837  /// change because of (un)refinement in the destination boundaries.
2838  /// Also deals with connection that do not longer exist because the
2839  /// destination boundary does no longer exist because of the distribution
2840  /// process
2841  void restore_polyline_connections_helper(
2842  TriangleMeshPolyLine* polyline_pt,
2843  Vector<TriangleMeshPolyLine*> &resume_initial_connection_polyline_pt,
2844  Vector<TriangleMeshPolyLine*> &resume_final_connection_polyline_pt);
2845 
2846  /// \short Resume the boundary connections that may have been
2847  /// suspended because the destination boundary is no part of the
2848  /// domain. The connections are no permanently suspended because if
2849  /// load balance takes place the destination boundary may be part of
2850  /// the new domain representation therefore the connection would
2851  /// exist
2852  void resume_boundary_connections(
2853  Vector<TriangleMeshPolyLine*> &resume_initial_connection_polyline_pt,
2854  Vector<TriangleMeshPolyLine*> &resume_final_connection_polyline_pt);
2855 
2856  /// \short Computes the associated vertex number on the destination
2857  /// boundary
2858  bool get_connected_vertex_number_on_dst_boundary(
2859  Vector<double> &vertex_coordinates,
2860  const unsigned &dst_b_id, unsigned &vertex_number);
2861 
2862  /// \short Helper function that performs the unrefinement process
2863  // on the specified boundary by using the provided vertices
2864  /// representation. Optional boolean is used to run it as test only (if
2865  /// true is specified as input) in which case vertex coordinates aren't
2866  /// actually modified. Returned boolean indicates if polyline was (or
2867  /// would have been -- if called with check_only=false) changed.
2868  bool unrefine_boundary(const unsigned &b,
2869  const unsigned &c,
2870  Vector<Vector<double> > &vector_bnd_vertices,
2871  double &unrefinement_tolerance,
2872  const bool &check_only = false);
2873 
2874  /// \short Helper function that performs the refinement process
2875  /// on the specified boundary by using the provided vertices
2876  /// representation. Optional boolean is used to run it as test only (if
2877  /// true is specified as input) in which case vertex coordinates aren't
2878  /// actually modified. Returned boolean indicates if polyline was (or
2879  /// would have been -- if called with check_only=false) changed.
2880  bool refine_boundary(Mesh* face_mesh_pt,
2881  Vector<Vector<double> > &vector_bnd_vertices,
2882  double &refinement_tolerance,
2883  const bool &check_only = false);
2884 
2885  // \short Helper function that applies the maximum length constraint
2886  // when it was specified. This will increase the number of points in
2887  // the current curve section in case that any segment on it does not
2888  // fulfils the requirement
2889  bool apply_max_length_constraint(Mesh* face_mesh_pt,
2891  &vector_bnd_vertices,
2892  double &max_length_constraint);
2893 
2894  /// \short Helper function that performs the unrefinement process on
2895  /// the specified boundary by using the provided vertices
2896  /// representation and the associated target area. Used only when the
2897  /// 'allow_automatic_creation_of_vertices_on_boundaries' flag is set to
2898  /// true.
2899  bool unrefine_boundary_constrained_by_target_area(
2900  const unsigned &b, const unsigned &c,
2901  Vector<Vector<double> > &vector_bnd_vertices,
2902  double &unrefinement_tolerance, Vector<double> &area_constraint);
2903 
2904  /// \short Helper function that performs the refinement process on
2905  /// the specified boundary by using the provided vertices
2906  /// representation and the associated elements target area. Used
2907  /// only when the 'allow_automatic_creation_of_vertices_on_boundaries'
2908  /// flag is set to true.
2909  bool refine_boundary_constrained_by_target_area(
2910  MeshAsGeomObject* mesh_geom_obj_pt,
2911  Vector<Vector<double> > &vector_bnd_vertices,
2912  double &refinement_tolerance, Vector<double> &area_constraint);
2913 
2914  /// \short Helper function that performs the unrefinement process
2915  /// on the specified boundary by using the provided vertices
2916  /// representation and the associated target area.
2917  /// NOTE: This is the version that applies unrefinement to shared
2918  /// boundaries
2919  bool unrefine_shared_boundary_constrained_by_target_area(
2920  const unsigned &b, const unsigned &c,
2921  Vector<Vector<double> > &vector_bnd_vertices,
2922  Vector<double> &area_constraint);
2923 
2924  /// \short Helper function that performs the refinement process
2925  /// on the specified boundary by using the provided vertices
2926  /// representation and the associated elements target area.
2927  /// NOTE: This is the version that applies refinement to shared
2928  /// boundaries
2929  bool refine_shared_boundary_constrained_by_target_area(
2930  Vector<Vector<double> > &vector_bnd_vertices,
2931  Vector<double> &area_constraint);
2932 
2933  // Flag that enables or disables boundary unrefinement (true by default)
2935  // Flag that enables or disables boundary refinement (true by default)
2937  // Flag that enables or disables boundary unrefinement (true by default)
2939  // Flag that enables or disables boundary unrefinement (true by default)
2941 
2942  // Set all the flags to true (the default values)
2944  {
2945  // All boundaries refinement and unrefinement are allowed by
2946  // default
2947  Do_boundary_unrefinement_constrained_by_target_areas = true;
2948  Do_boundary_refinement_constrained_by_target_areas = true;
2949  Do_shared_boundary_unrefinement_constrained_by_target_areas = true;
2950  Do_shared_boundary_refinement_constrained_by_target_areas = true;
2951  }
2952 
2953 #ifdef OOMPH_HAS_MPI
2954  // \short Stores the nodes in the boundaries in the same order in all the
2955  // processors
2956  // Sorted_shared_boundary_node_pt[bnd_id][i-th node] = Node*
2957  // It is a map since the boundary id may not start at zero
2958  std::map<unsigned, Vector<Node*> > Sorted_shared_boundary_node_pt;
2959 
2960  // \short Sort the nodes on shared boundaries so that the processors
2961  // that share a boundary agree with the order of the nodes on the
2962  // boundary
2963  void sort_nodes_on_shared_boundaries();
2964 
2965  // \short Re-establish the shared boundary elements after the
2966  // adaptation process (the updating of shared nodes is optional and
2967  // performed by default)
2968  void reset_shared_boundary_elements_and_nodes(const bool flush_elements
2969  = true,
2970  const bool update_elements
2971  = true,
2972  const bool flush_nodes
2973  = true,
2974  const bool update_nodes
2975  = true);
2976 
2977  // \short In charge of. re-establish the halo(ed) scheme on all processors.
2978  // Sends info. to create halo elements and nodes on the processors
2979  // that need it. It uses and all to all communication strategy therefore
2980  // must be called on all processors.
2981  void reset_halo_haloed_scheme();
2982 
2983  // \short Compute the names of the nodes on shared boundaries in
2984  // this (my_rank) processor with other processors. Also compute the
2985  // names of nodes on shared boundaries of other processors with
2986  // other processors (useful when there is an element that requires
2987  // to be sent to this (my_rank) processor because there is a shared
2988  // node between this (my_rank) and other processors BUT there is
2989  // not a shared boundary between this and the other processor
2990  void compute_global_node_names_and_shared_nodes(
2991  Vector<Vector<Vector<std::map<unsigned, Node*> > > >
2992  &other_proc_shd_bnd_node_pt,
2993  Vector<Vector<Vector<unsigned> > > &global_node_names,
2994  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
2995  Vector<Node*> &global_shared_node_pt);
2996 
2997  // \short Get the original boundaries to which is associated each
2998  // shared node, and send the info. to the related processors. We
2999  // need to do this so that at the reset of halo(ed) info. stage,
3000  // the info. be already updated.
3001  void send_boundary_node_info_of_shared_nodes(
3002  Vector<Vector<Vector<unsigned> > > &global_node_names,
3003  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
3004  Vector<Node*> &global_shared_node_pt);
3005 
3006  // \short In charge of creating additional halo(ed) elements on
3007  // those processors that have no shared boundaries in common but have
3008  // shared nodes
3009  void reset_halo_haloed_scheme_helper(
3010  Vector<Vector<Vector<std::map<unsigned, Node*> > > >
3011  &other_proc_shd_bnd_node_pt,
3012  Vector<Vector<Node *> > &iproc_currently_created_nodes_pt,
3013  Vector<Vector<Vector<unsigned> > > &global_node_names,
3014  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
3015  Vector<Node*> &global_shared_node_pt);
3016 
3017  // ====================================================================
3018  // Methods for load balancing
3019  // ====================================================================
3020 
3021 //#define ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION_LOAD_BALANCE
3022 
3023  // *********************************************************************
3024  // BEGIN: Methods to perform load balance
3025  // *********************************************************************
3026 
3027  // Check if necessary to add the element to the new domain or if it has
3028  // been previously added
3030  Vector<FiniteElement*> &new_elements_on_domain,
3031  FiniteElement* &ele_pt)
3032  {
3033  // Get the number of elements currently added to the new domain
3034  const unsigned nnew_elements_on_domain = new_elements_on_domain.size();
3035 
3036  // Flag to state if has been added or not
3037  bool already_on_new_domain = false;
3038  unsigned new_domain_ele_index = 0;
3039 
3040  for (unsigned e = 0; e < nnew_elements_on_domain; e++)
3041  {
3042  if (ele_pt == new_elements_on_domain[e])
3043  {
3044  // It's already there, so...
3045  already_on_new_domain = true;
3046  // ...set the index of this element
3047  new_domain_ele_index = e;
3048  break;
3049  }
3050  }
3051 
3052  // Has it been found?
3053  if (!already_on_new_domain)
3054  {
3055  // Not found, so add it:
3056  new_elements_on_domain.push_back(ele_pt);
3057  // Return the index where it's just been added
3058  return nnew_elements_on_domain;
3059  }
3060  else
3061  {
3062  // Return the index where it was found
3063  return new_domain_ele_index;
3064  }
3065 
3066  }
3067 
3068  /// \short Helper function to get the required elemental information from
3069  /// the element to be sent. This info. involves the association of
3070  /// the element to a boundary or region, and if its part of the
3071  /// halo(ed) elements within a processor
3072  void get_required_elemental_information_load_balance_helper(
3073  unsigned& iproc,
3074  Vector<Vector<FiniteElement*> > &f_haloed_ele_pt,
3075  FiniteElement* ele_pt);
3076 
3077  // Check if necessary to add the node to the new domain or if it has been
3078  // already added
3080  Vector<Node*> &new_nodes_on_domain,
3081  Node*& node_pt)
3082  {
3083  // Get the number of nodes currently added to the new domain
3084  const unsigned nnew_nodes_on_domain = new_nodes_on_domain.size();
3085 
3086  // Flag to state if has been added or not
3087  bool already_on_new_domain = false;
3088  unsigned new_domain_node_index = 0;
3089 
3090  for (unsigned n = 0; n < nnew_nodes_on_domain; n++)
3091  {
3092  if (node_pt == new_nodes_on_domain[n])
3093  {
3094  // It's already there, so...
3095  already_on_new_domain = true;
3096  // ...set the index of this element
3097  new_domain_node_index = n;
3098  break;
3099  }
3100  }
3101 
3102  // Has it been found?
3103  if (!already_on_new_domain)
3104  {
3105  // Not found, so add it:
3106  new_nodes_on_domain.push_back(node_pt);
3107  // Return the index where it's just been added
3108  return nnew_nodes_on_domain;
3109  }
3110  else
3111  {
3112  // Return the index where it was found
3113  return new_domain_node_index;
3114  }
3115 
3116  }
3117 
3118  /// \short Helper function to add haloed node
3119  void add_node_load_balance_helper(unsigned& iproc,
3121  &f_halo_ele_pt,
3122  Vector<Node*> &new_nodes_on_domain,
3123  Node* nod_pt);
3124 
3125  /// \short Helper function to get the required nodal information
3126  /// from an haloed node so that a fully-functional node (and
3127  /// therefore element) can be created on the receiving process
3128  /// (this is the specific version for the load balance strategy,
3129  /// the difference with the original method is that it checks if
3130  /// the node is on a shared boundary no associated with the current
3131  /// processor --my_rank--, or in a haloed element from other
3132  /// processors
3133  void get_required_nodal_information_load_balance_helper(
3134  Vector<Vector<FiniteElement*> > &f_halo_ele_pt,
3135  unsigned& iproc,
3136  Node* nod_pt);
3137 
3138  /// \short Helper function to create elements on the loop
3139  /// process based on the info received in
3140  /// send_and_received_elements_nodes_info
3141  void create_element_load_balance_helper(unsigned& iproc,
3142  Vector<Vector<FiniteElement*> > &f_haloed_ele_pt,
3143  Vector<Vector<std::map<unsigned,FiniteElement*> > >
3144  &received_old_haloed_element_pt,
3145  Vector<FiniteElement*> &new_elements_on_domain,
3146  Vector<Node*> &new_nodes_on_domain,
3147  Vector<Vector<Vector<std::map<unsigned, Node*> > > >
3148  &other_proc_shd_bnd_node_pt,
3149  Vector<Vector<Vector<unsigned> > > &global_node_names,
3150  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
3151  Vector<Node*> &global_shared_node_pt);
3152 
3153  /// \short Helper function to create elements on the loop
3154  /// process based on the info received in
3155  /// send_and_received_elements_nodes_info
3156  /// This function is in charge of verify if the element is associated
3157  /// to a boundary and associate to it if that is the case
3158  void add_element_load_balance_helper(const unsigned &iproc,
3159  Vector<Vector<std::map<
3160  unsigned,FiniteElement*> > >
3161  &received_old_haloed_element_pt,
3162  FiniteElement* ele_pt);
3163 
3164  /// \short Helper function to add a new node from load balance
3165  void add_received_node_load_balance_helper(Node* &new_nod_pt,
3166  Vector<Vector<FiniteElement*> > &f_haloed_ele_pt,
3167  Vector<Vector<std::map<unsigned,FiniteElement*> > >
3168  &received_old_haloed_element_pt,
3169  Vector<Node*> &new_nodes_on_domain,
3170  Vector<Vector<Vector<std::map<unsigned, Node*> > > >
3171  &other_proc_shd_bnd_node_pt,
3172  unsigned& iproc, unsigned& node_index,
3173  FiniteElement* const &new_el_pt,
3174  Vector<Vector<Vector<unsigned> > > &global_node_names,
3175  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
3176  Vector<Node*> &global_shared_node_pt);
3177 
3178  /// \short Helper function which constructs a new node (on an
3179  /// element) with the information sent from the load balance
3180  /// process
3181  void construct_new_node_load_balance_helper(Node* &new_nod_pt,
3182  Vector<Vector<FiniteElement*> > &f_haloed_ele_pt,
3183  Vector<Vector<std::map<unsigned,FiniteElement*> > >
3184  &received_old_haloed_element_pt,
3185  Vector<Node*> &new_nodes_on_domain,
3186  Vector<Vector<Vector<std::map<unsigned, Node*> > > >
3187  &other_proc_shd_bnd_node_pt,
3188  unsigned& iproc, unsigned& node_index,
3189  FiniteElement* const &new_el_pt,
3190  Vector<Vector<Vector<unsigned> > > &global_node_names,
3191  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
3192  Vector<Node*> &global_shared_node_pt);
3193 
3194  // *********************************************************************
3195  // END: Methods to perform load balance
3196  // *********************************************************************
3197 
3198  // *********************************************************************
3199  // Start communication variables
3200  // *********************************************************************
3201  /// \short Vector of flat-packed doubles to be communicated with
3202  /// other processors
3204 
3205  /// \short Counter used when processing vector of flat-packed
3206  /// doubles
3208 
3209  /// \short Vector of flat-packed unsigneds to be communicated with
3210  /// other processors
3212 
3213  /// \short Counter used when processing vector of flat-packed
3214  /// unsigneds
3216 
3217 #ifdef ANNOTATE_REFINEABLE_TRIANGLE_MESH_COMMUNICATION
3218  // Temporary vector of strings to enable full annotation of
3219  // RefineableTriangleMesh comms
3221 #endif
3222 
3223  // *********************************************************************
3224  // End communication variables
3225  // *********************************************************************
3226 
3227  // *********************************************************************
3228  // Start communication functions
3229  // *********************************************************************
3230 
3231  // Check if necessary to add the element as haloed or if it has been
3232  // previously added to the haloed scheme
3233  unsigned try_to_add_root_haloed_element_pt(const unsigned& p,
3234  GeneralisedElement*& el_pt)
3235  {
3236  // Loop over current storage
3237  unsigned n_haloed=this->nroot_haloed_element(p);
3238 
3239  // Is this already an haloed element?
3240  bool already_haloed_element=false;
3241  unsigned haloed_el_index=0;
3242  for (unsigned eh=0;eh<n_haloed;eh++)
3243  {
3244  if (el_pt==this->root_haloed_element_pt(p, eh))
3245  {
3246  // It's already there, so...
3247  already_haloed_element=true;
3248  // ...set the index of this element
3249  haloed_el_index=eh;
3250  break;
3251  }
3252  }
3253 
3254  // Has it been found?
3255  if (!already_haloed_element)
3256  {
3257  // Not found, so add it:
3258  this->add_root_haloed_element_pt(p, el_pt);
3259  // Return the index where it's just been added
3260  return n_haloed;
3261  }
3262  else
3263  {
3264  // Return the index where it was found
3265  return haloed_el_index;
3266  }
3267  }
3268 
3269  // Check if necessary to add the node as haloed or if it has been
3270  // previously added to the haloed scheme
3271  unsigned try_to_add_haloed_node_pt(const unsigned& p, Node*& nod_pt)
3272  {
3273  // Loop over current storage
3274  unsigned n_haloed_nod=this->nhaloed_node(p);
3275 
3276  // Is this already an haloed node?
3277  bool is_an_haloed_node=false;
3278  unsigned haloed_node_index=0;
3279  for (unsigned k=0;k<n_haloed_nod;k++)
3280  {
3281  if (nod_pt==this->haloed_node_pt(p,k))
3282  {
3283  is_an_haloed_node=true;
3284  haloed_node_index=k;
3285  break;
3286  }
3287  }
3288 
3289  // Has it been found?
3290  if (!is_an_haloed_node)
3291  {
3292  // Not found, so add it
3293  this->add_haloed_node_pt(p, nod_pt);
3294  // Return the index where it's just been added
3295  return n_haloed_nod;
3296  }
3297  else
3298  {
3299  // Return the index where it was found
3300  return haloed_node_index;
3301  }
3302  }
3303 
3304  /// \short Helper function to get the required elemental information from
3305  /// an haloed element. This info. involves the association of the element
3306  /// to a boundary or region.
3307  void get_required_elemental_information_helper(unsigned& iproc,
3308  FiniteElement* ele_pt);
3309 
3310  /// \short Helper function to get the required nodal information
3311  /// from a haloed node so that a fully-functional halo node (and
3312  /// therefore element) can be created on the receiving process
3313  void get_required_nodal_information_helper(unsigned& iproc, Node* nod_pt);
3314 
3315  /// \short Helper function to add haloed node
3316  void add_haloed_node_helper(unsigned& iproc, Node* nod_pt);
3317 
3318  /// \short Helper function to send back halo and haloed information
3319  void send_and_receive_elements_nodes_info(int& send_proc, int &recv_proc);
3320 
3321  /// \short Helper function to create (halo) elements on the loop
3322  /// process based on the info received in send_and_received_located_info
3323  void create_halo_element(unsigned &iproc,
3324  Vector<Node*> &new_nodes_on_domain,
3325  Vector<Vector<Vector<std::map<unsigned,
3326  Node*> > > > &other_proc_shd_bnd_node_pt,
3328  &global_node_names,
3329  std::map<Vector<unsigned>, unsigned>
3330  &node_name_to_global_index,
3331  Vector<Node*> &global_shared_node_pt);
3332 
3333  /// \short Helper function to create (halo) elements on the loop
3334  /// process based on the info received in send_and_received_located_info
3335  /// This function is in charge of verify if the element is associated to
3336  /// a boundary
3337  void add_halo_element_helper(unsigned& iproc, FiniteElement* ele_pt);
3338 
3339  /// \short Helper function to add halo node
3340  void add_halo_node_helper(Node* &new_nod_pt,
3341  Vector<Node*> &new_nodes_on_domain,
3342  Vector<Vector<Vector<std::map<unsigned,
3343  Node*> > > >
3344  &other_proc_shd_bnd_node_pt,
3345  unsigned& iproc,
3346  unsigned& node_index,
3347  FiniteElement* const &new_el_pt,
3349  &global_node_names,
3350  std::map<Vector<unsigned>, unsigned>
3351  &node_name_to_global_index,
3352  Vector<Node*> &global_shared_node_pt);
3353 
3354  /// \short Helper function which constructs a new halo node
3355  /// (on an element) with the information sent from the haloed process
3356  void construct_new_halo_node_helper(Node* &new_nod_pt,
3357  Vector<Node*> &new_nodes_on_domain,
3358  Vector<Vector<Vector<std::map<unsigned,
3359  Node*> > > >
3360  &other_proc_shd_bnd_node_pt,
3361  unsigned& iproc,
3362  unsigned& node_index,
3363  FiniteElement* const &new_el_pt,
3365  &global_node_names,
3366  std::map<Vector<unsigned>, unsigned>
3367  &node_name_to_global_index,
3368  Vector<Node*> &global_shared_node_pt);
3369 
3370  //Helper function that assigns/updates the references to the node
3371  //so that it can be found with any other reference. The return
3372  //value indicates whether or not a node was found on the same
3373  //reference
3374  void update_other_proc_shd_bnd_node_helper
3375  (Node* &new_nod_pt,
3376  Vector<Vector<Vector<std::map<unsigned, Node*> > > >
3377  &other_proc_shd_bnd_node_pt,
3378  Vector<unsigned> &other_processor_1,
3379  Vector<unsigned> &other_processor_2,
3380  Vector<unsigned> &other_shared_boundaries,
3381  Vector<unsigned> &other_indexes,
3382  Vector<Vector<Vector<unsigned> > > &global_node_names,
3383  std::map<Vector<unsigned>, unsigned> &node_name_to_global_index,
3384  Vector<Node*> &global_shared_node_pt);
3385 
3386  // *********************************************************************
3387  // End Communication funtions
3388  // *********************************************************************
3389 
3390 #endif // #ifdef OOMPH_HAS_MPI
3391 
3392  /// \short Helper function that updates the input polygon's PSLG
3393  /// by using the end-points of elements from FaceMesh(es) that are
3394  /// constructed for the boundaries associated with the segments of the
3395  /// polygon. Optional boolean is used to run it as test only (if
3396  /// true is specified as input) in which case polygon isn't actually
3397  /// modified. Returned boolean indicates if polygon was (or would have
3398  /// been -- if called with check_only=false) changed.
3399  bool update_polygon_using_face_mesh(TriangleMeshPolygon* polygon_pt,
3400  const bool& check_only=false);
3401 
3402  /// \short Helper function that updates the input open curve by using
3403  /// end-points of elements from FaceMesh(es) that are constructed for the
3404  /// boundaries associated with the polylines. Optional boolean is used to
3405  /// run it as test only (if true is specified as input) in which case the
3406  /// polylines are not actually modified. Returned boolean indicates if
3407  /// polylines were (or would have been -- if called with check_only=false)
3408  /// changed.
3409  bool update_open_curve_using_face_mesh(
3410  TriangleMeshOpenCurve* open_polyline_pt,
3411  const bool& check_only=false);
3412 
3413  /// \short Generate a new PSLG representation of the inner hole
3414  /// boundaries. Optional boolean is used to run it as test only (if
3415  /// true is specified as input) in which case PSLG isn't actually
3416  /// modified. Returned boolean indicates if PSLG was (or would have
3417  /// been -- if called with check_only=false) changed.
3418  virtual bool surface_remesh_for_inner_hole_boundaries(
3419  Vector<Vector<double> > &internal_point_coord,
3420  const bool& check_only=false);
3421 
3422  /// \short Snap the boundary nodes onto any curvilinear boundaries
3423  void snap_nodes_onto_boundary(RefineableTriangleMesh<ELEMENT>* &new_mesh_pt,
3424  const unsigned &b);
3425 
3426  /// \short Helper function
3427  /// Creates an unsorted face mesh representation from the specified
3428  /// boundary id. It means that the elements are not sorted along the
3429  /// boundary
3430  void create_unsorted_face_mesh_representation(
3431  const unsigned &boundary_id,
3432  Mesh* face_mesh_pt);
3433 
3434  /// \short Helper function
3435  /// Creates a sorted face mesh representation of the specified PolyLine
3436  /// It means that the elements are sorted along the boundary
3437  /// It also returns a map that indicated the inverted elements
3438  void create_sorted_face_mesh_representation(
3439  const unsigned &boundary_id,
3440  Mesh* face_mesh_pt,
3441  std::map<FiniteElement*, bool> &is_inverted,
3442  bool &inverted_face_mesh);
3443 
3444  /// \short Helper function to construct face mesh representation of all
3445  /// polylines, possibly with segments re-distributed between polylines
3446  /// to maintain an appxroximately even sub-division of the polygon
3447  void get_face_mesh_representation(TriangleMeshPolygon* polygon_pt,
3448  Vector<Mesh*>& face_mesh_pt);
3449 
3450  /// \short Helper function to construct face mesh representation of
3451  /// open curves
3452  void get_face_mesh_representation(
3453  TriangleMeshOpenCurve* open_polyline_pt,
3454  Vector<Mesh*>& face_mesh_pt);
3455 
3456  /// \short Updates the polylines representation after restart
3457  void update_polygon_after_restart(TriangleMeshPolygon* &polygon_pt);
3458 
3459  /// \short Updates the open curve representation after restart
3460  void update_open_curve_after_restart(TriangleMeshOpenCurve* &open_curve_pt);
3461 
3462  /// \short Updates the polylines using the elements area as constraint for
3463  /// the number of points along the boundaries
3464  bool update_polygon_using_elements_area(
3465  TriangleMeshPolygon* &polygon_pt, const Vector<double> &target_area);
3466 
3467  /// \short Updates the open curve but using the elements area instead
3468  /// of the default refinement and unrefinement methods
3469  bool update_open_curve_using_elements_area(
3470  TriangleMeshOpenCurve* &open_curve_pt, const Vector<double> &target_area);
3471 
3472 #ifdef OOMPH_HAS_MPI
3473  /// \short Updates the polylines using the elements area as
3474  /// constraint for the number of points along the boundaries
3475  bool update_shared_curve_using_elements_area(
3476  Vector<TriangleMeshPolyLine*> &vector_polyline_pt,
3477  const Vector<double> &target_areas);
3478 
3479  /// \short Updates the shared polylines representation after restart
3480  void update_shared_curve_after_restart(Vector<TriangleMeshPolyLine*>
3481  &vector_polyline_pt);
3482 
3483 #endif // #ifdef OOMPH_HAS_MPI
3484 
3485  /// Helper function to initialise data associated with adaptation
3487  {
3488  // Number of bins in the x-direction
3489  // when transferring target areas by bin method
3490  this->Nbin_x_for_area_transfer=100;
3491 
3492  // Number of bins in the y-direction
3493  // when transferring target areas by bin method
3494  this->Nbin_y_for_area_transfer=100;
3495 
3496  /// Number of bins in the x-direction when projecting the solution
3497  /// from the old mesh into the new mesh
3498  this->Nbin_x_for_projection=100;
3499 
3500  /// Number of bins in the y-direction when projecting the solution
3501  /// from the old mesh into the new mesh
3502  this->Nbin_y_for_projection=100;
3503 
3504  // Set max and min targets for adaptation
3505  this->Max_element_size=1.0;
3506  this->Min_element_size=0.001;
3507  this->Min_permitted_angle=15.0;
3508 
3509  // By default we want to do projection
3510  this->Disable_projection=false;
3511 
3512  // Use by default an iterative solver for the projection problem
3513  this->Use_iterative_solver_for_projection=true;
3514 
3515  // Set the defaul value for printing level adaptation (default 0)
3516  this->Print_timings_level_adaptation=0;
3517 
3518  // Set the defaul value for printing level load balance (default 0)
3519  this->Print_timings_level_load_balance=0;
3520 
3521  // By default we want no info. about timings for transferring of
3522  // target areas
3523  this->Print_timings_transfering_target_areas=false;
3524 
3525  // By default we want no info. about timings for projection
3526  this->Print_timings_projection=false;
3527 
3528  // Initialise function pointer to function that updates the
3529  // mesh following the snapping of boundary nodes to the
3530  // boundaries (e.g. to move boundary nodes very slightly
3531  // to satisfy volume constraints)
3532  Mesh_update_fct_pt=0;
3533 
3534  // Initialise function point for update of internal hole
3535  // point to null
3536  Internal_hole_point_update_fct_pt=0;
3537 
3538  }
3539 
3540 #ifdef OOMPH_HAS_TRIANGLE_LIB
3541 
3542  /// \short Build a new TriangulateIO object from previous TriangulateIO
3543  /// based on target area for each element
3544  void refine_triangulateio(TriangulateIO& triangulate_io,
3545  const Vector<double> &target_area,
3546  TriangulateIO &triangle_refine);
3547 
3548 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3549 
3550  /// \short Compute target area based on the element's error and the
3551  /// error target; return minimum angle (in degrees)
3552  double compute_area_target(const Vector<double> &elem_error,
3553  Vector<double> &target_area)
3554  {
3555  double min_angle=DBL_MAX;
3556  unsigned count_unrefined=0;
3557  unsigned count_refined=0;
3558  this->Nrefinement_overruled=0;
3559 
3560  // Record max. area constraint set by region
3561  std::map<FiniteElement*,double> max_area_from_region;
3562  for (std::map<unsigned, double>::iterator it=this->Regions_areas.begin();
3563  it!=this->Regions_areas.end();it++)
3564  {
3565  unsigned r=(*it).first;
3566  unsigned nel=this->nregion_element(r);
3567  for(unsigned e=0;e<nel;e++)
3568  {
3569  max_area_from_region[this->region_element_pt(r,e)]=(*it).second;
3570  }
3571  }
3572 
3573  unsigned nel=this->nelement();
3574  for (unsigned e=0;e<nel;e++)
3575  {
3576  // Get element
3577  FiniteElement* el_pt=this->finite_element_pt(e);
3578 
3579  // Area
3580  double area=el_pt->size();
3581 
3582  // Min angle based on vertex coordinates
3583  // (vertices are enumerated first)
3584  double ax=el_pt->node_pt(0)->x(0);
3585  double ay=el_pt->node_pt(0)->x(1);
3586 
3587  double bx=el_pt->node_pt(1)->x(0);
3588  double by=el_pt->node_pt(1)->x(1);
3589 
3590  double cx=el_pt->node_pt(2)->x(0);
3591  double cy=el_pt->node_pt(2)->x(1);
3592 
3593  // Min angle
3594  double angle0=
3595  acos(((ax-cx)*(bx-cx)+(ay-cy)*(by-cy))/
3596  (sqrt((ax-cx)*(ax-cx)+(ay-cy)*(ay-cy))*
3597  sqrt((bx-cx)*(bx-cx)+(by-cy)*(by-cy))))*
3599  min_angle=std::min(min_angle,angle0);
3600 
3601  double angle1=
3602  acos(((ax-bx)*(cx-bx)+(ay-by)*(cy-by))/
3603  (sqrt((ax-bx)*(ax-bx)+(ay-by)*(ay-by))*
3604  sqrt((cx-bx)*(cx-bx)+(cy-by)*(cy-by))))*
3606  min_angle=std::min(min_angle,angle1);
3607 
3608  double angle2=180.0-angle0-angle1;
3609  min_angle=std::min(min_angle,angle2);
3610 
3611  // Mimick refinement in tree-based procedure: Target areas
3612  // for elements that exceed permitted error is 1/3 of their
3613  // current area, corresponding to a uniform sub-division.
3614  double size_ratio=3.0;
3615  if (elem_error[e]>max_permitted_error())
3616  {
3617  // Reduce area
3618  target_area[e]=std::max(area/size_ratio,Min_element_size);
3619 
3620  //...but also make sure we're below the max element size
3621  target_area[e]=std::min(target_area[e],Max_element_size);
3622 
3623  if (target_area[e]!=Min_element_size)
3624  {
3625  count_refined++;
3626  }
3627  else
3628  {
3629  this->Nrefinement_overruled++;
3630  }
3631  }
3632  else if (elem_error[e]<min_permitted_error())
3633  {
3634  // Increase the area
3635  target_area[e]=std::min(size_ratio*area,Max_element_size);
3636 
3637  //...but also make sure we're above the min element size
3638  target_area[e]=std::max(target_area[e],Min_element_size);
3639 
3640  if (target_area[e]!=Max_element_size)
3641  {
3642  count_unrefined++;
3643  }
3644  }
3645  else
3646  {
3647  // Leave it alone but enforce size limits
3648  double area_leave_alone = std::max(area,Min_element_size);
3649  target_area[e] = std::min(area_leave_alone,Max_element_size);
3650  }
3651 
3652  // Enforce max areas from regions
3653  std::map<FiniteElement*,double>::iterator it=
3654  max_area_from_region.find(el_pt);
3655  if (it!=max_area_from_region.end())
3656  {
3657  target_area[e]=std::min(target_area[e],(*it).second);
3658  }
3659 
3660  }
3661 
3662 
3663  // Tell everybody
3664  this->Nrefined=count_refined;
3665  this->Nunrefined=count_unrefined;
3666 
3667  if (this->Nrefinement_overruled!=0)
3668  {
3669  oomph_info
3670  << "\nNOTE: Refinement of "
3671  << this->Nrefinement_overruled << " elements was "
3672  << "overruled \nbecause the target area would have "
3673  << "been below \nthe minimum permitted area of "
3674  << Min_element_size
3675  << ".\nYou can change the minimum permitted area with the\n"
3676  << "function RefineableTriangleMesh::min_element_size().\n\n";
3677  }
3678  return min_angle;
3679  }
3680 
3681  /// \short Number of bins in the x-direction
3682  /// when transferring target areas by bin method.
3684 
3685  /// \short Number of bins in the y-direction
3686  /// when transferring target areas by bin method.
3688 
3689  /// \short Number of bins in the x-direction when projecting the
3690  /// solution from the old mesh into the new mesh
3692 
3693  /// \short Number of bins in the y-direction when projecting the
3694  /// solution from the old mesh into the new mesh
3696 
3697  /// Max permitted element size
3699 
3700  /// Min permitted element size
3702 
3703  /// Min angle before remesh gets triggered
3705 
3706  /// Enable/disable solution projection during adaptation
3708 
3709  /// Flag to indicate whether to use or not an iterative solver (CG
3710  /// with diagonal preconditioned) for the projection problem
3712 
3713  /// Enable/disable printing timings for transfering target areas
3715 
3716  /// Enable/disable printing timings for projection
3718 
3719  /// The printing level for adaptation
3721 
3722  /// The printing level for load balance
3724 
3725  /// \short Function pointer to function that updates the
3726  /// mesh following the snapping of boundary nodes to the
3727  /// boundaries (e.g. to move boundary nodes very slightly
3728  /// to satisfy volume constraints)
3729  MeshUpdateFctPt Mesh_update_fct_pt;
3730 
3731  /// \short Function pointer to function that can be set
3732  /// to update the position of the central point in internal
3733  /// holes
3734  InternalHolePointUpdateFctPt Internal_hole_point_update_fct_pt;
3735 
3736  };
3737 
3738 
3739 /////////////////////////////////////////////////////////////////////
3740 /////////////////////////////////////////////////////////////////////
3741 /////////////////////////////////////////////////////////////////////
3742 
3743 
3744 //=========================================================================
3745 /// Unstructured Triangle Mesh upgraded to solid mesh
3746 //=========================================================================
3747  template<class ELEMENT>
3748  class SolidTriangleMesh : public virtual TriangleMesh<ELEMENT>,
3749  public virtual SolidMesh
3750  {
3751 
3752  public:
3753 
3754 #ifdef OOMPH_HAS_TRIANGLE_LIB
3755 
3756  /// \short Build mesh, based on closed curve that specifies
3757  /// the outer boundary of the domain and any number of internal
3758  /// clsed curves. Specify target area for uniform element size.
3759  SolidTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters,
3760  TimeStepper* time_stepper_pt=&Mesh::Default_TimeStepper)
3761  : TriangleMesh<ELEMENT>(triangle_mesh_parameters,
3762  time_stepper_pt)
3763  {
3764  //Assign the Lagrangian coordinates
3766  }
3767 
3768 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3769 
3770  SolidTriangleMesh(const std::string& node_file_name,
3771  const std::string& element_file_name,
3772  const std::string& poly_file_name,
3773  TimeStepper* time_stepper_pt=
3775  const bool &allow_automatic_creation_of_vertices_on_boundaries=true)
3776  : TriangleMesh<ELEMENT>(node_file_name,
3777  element_file_name,
3778  poly_file_name,
3779  time_stepper_pt,
3780  allow_automatic_creation_of_vertices_on_boundaries)
3781  {
3782  //Assign the Lagrangian coordinates
3784  }
3785 
3786  /// Empty Destructor
3787  virtual ~SolidTriangleMesh() { }
3788  };
3789 
3790 
3791 ////////////////////////////////////////////////////////////////////////
3792 ////////////////////////////////////////////////////////////////////////
3793 ////////////////////////////////////////////////////////////////////////
3794 
3795 #ifdef OOMPH_HAS_TRIANGLE_LIB
3796 
3797 //=========================================================================
3798 /// Unstructured refineable Triangle Mesh upgraded to solid mesh
3799 //=========================================================================
3800  template<class ELEMENT>
3802  public virtual RefineableTriangleMesh<ELEMENT>,
3803  public virtual SolidMesh
3804  {
3805 
3806  public:
3807 
3808  /// \short Build mesh, based on the specifications on
3809  /// TriangleMeshParameter
3811  TimeStepper* time_stepper_pt=
3813  : TriangleMesh<ELEMENT>(triangle_mesh_parameters,
3814  time_stepper_pt),
3815  RefineableTriangleMesh<ELEMENT>(triangle_mesh_parameters,
3816  time_stepper_pt)
3817  {
3818  //Assign the Lagrangian coordinates
3820  }
3821 
3822  /// \short Build mesh from specified triangulation and
3823  /// associated target areas for elements in it.
3825  TriangulateIO& triangulate_io,
3826  TimeStepper* time_stepper_pt=
3828  const bool &use_attributes=false,
3829  const bool
3830  &allow_automatic_creation_of_vertices_on_boundaries=true,
3831  OomphCommunicator* comm_pt = 0) :
3832  RefineableTriangleMesh<ELEMENT>(target_area,
3833  triangulate_io,
3834  time_stepper_pt,
3835  use_attributes,
3836  allow_automatic_creation_of_vertices_on_boundaries,
3837  comm_pt)
3838  {
3839  //Assign the Lagrangian coordinates
3841  }
3842 
3843  /// Empty Destructor
3845  };
3846 
3847 #endif // #ifdef OOMPH_HAS_TRIANGLE_LIB
3848 
3849 }
3850 
3851 #endif
const unsigned nshared_boundaries(const unsigned &p, const unsigned &q) const
bool Print_timings_transfering_target_areas
Enable/disable printing timings for transfering target areas.
A Generalised Element class.
Definition: elements.h:76
RefineableSolidTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameter.
void set_communicator_pt(OomphCommunicator *comm_pt)
Definition: mesh.h:1332
TriangulateIO Triangulateio
TriangulateIO representation of the mesh.
void set_geom_objects_and_coordinate_limits_for_close_curve(TriangleMeshClosedCurve *input_closed_curve_pt)
Stores the geometric objects associated to the curve sections that compound the closed curve...
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...
void disable_timings_projection()
Disables info. and timings for projection.
Vector< Vector< unsigned > > shared_boundaries_ids(const unsigned &p) const
Vector< unsigned > & shared_boundaries_ids(const unsigned &p, const unsigned &q)
Node *& boundary_segment_node_pt(const unsigned &b, const unsigned &s, const unsigned &n)
Return pointer to node n on boundary b.
std::map< unsigned, Vector< FiniteElement * > > Region_element_pt
Vector of elements in each region differentiated by attribute (the key of the map is the attribute) ...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
double Max_element_size
Max permitted element size.
void set_geom_objects_and_coordinate_limits_for_open_curve(TriangleMeshOpenCurve *input_open_curve_pt)
Stores the geometric objects associated to the curve sections that compound the open curve...
Triangle Mesh that is based on input files generated by the triangle mesh generator Triangle...
double * pointlist
Pointer to list of points x coordinate followed by y coordinate.
double & min_element_size()
Min element size allowed during adaptation.
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)
virtual void load_balance(const Vector< unsigned > &target_domain_for_local_non_halo_element)
Virtual function to perform the load balance routines.
Vector< TriangleMeshOpenCurve * > Internal_open_curves_pt
Internal boundaries.
const unsigned nshared_boundary_node(const unsigned &b)
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:194
void enable_timings_projection()
Enables info. and timings for projection.
Unstructured refineable Triangle Mesh upgraded to solid mesh.
std::map< unsigned, Vector< GeneralisedElement * > > Root_halo_element_pt
Map of vectors holding the pointers to the root halo elements.
Definition: mesh.h:111
std::map< unsigned, Vector< double > > & regions_coordinates()
Helper function for getting access to the regions coordinates.
Vector< TriangleMeshClosedCurve * > Internal_closed_curve_pt
Internal closed boundaries.
std::map< unsigned, Vector< Node * > > Haloed_node_pt
Map of vectors holding the pointers to the haloed nodes.
Definition: mesh.h:120
TriangleMeshPolygon * closed_curve_to_polygon_helper(TriangleMeshClosedCurve *closed_curve_pt, unsigned &max_bnd_id_local)
Helper function that returns a polygon representation for the given closed curve, it also computes th...
std::map< unsigned, double > Regions_areas
Target areas for regions; defaults to 0.0 which (luckily) implies "no specific target area" for trian...
Vector< TriangleMeshOpenCurve * > Internal_open_curve_pt
Vector of open polylines that define internal curves.
std::map< unsigned, Vector< GeneralisedElement * > > External_halo_element_pt
Map of vectors holding the pointers to the external halo elements.
Definition: mesh.h:137
void enable_projection()
Enables the solution projection step during adaptation.
virtual ~RefineableTriangleMesh()
Empty Destructor.
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...
void initialise_adaptation_data()
Helper function to initialise data associated with adaptation.
static Steady< 0 > Default_TimeStepper
Default Steady Timestepper, to be used in default arguments to Mesh constructors. ...
Definition: mesh.h:85
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
void disable_shared_boundary_unrefinement_constrained_by_target_areas()
bool Disable_projection
Enable/disable solution projection during adaptation.
MeshUpdateFctPt & mesh_update_fct_pt()
Access to function pointer to function that updates the mesh following the snapping of boundary nodes...
std::map< unsigned, unsigned > Shared_boundary_overlaps_internal_boundary
Vector< TriangleMeshClosedCurve * > & outer_boundary_pt()
Helper function for getting access to the outer boundary.
void re_scale_re_assigned_initial_zeta_values_for_internal_boundary(const unsigned &b)
Re-scale the re-assigned zeta values for the boundary nodes, apply only for internal boundaries...
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:201
void set_target_area_for_region(const unsigned &i, const double &area)
Helper function to specify target area for region.
Information for documentation of results: Directory and file number to enable output in the form RESL...
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 ...
std::map< unsigned, TriangleMeshCurveSection * > Boundary_curve_section_pt
A map that stores the associated curve section of the specified boundary id.
void update_polyline_representation_from_restart()
Method used to update the polylines representation after restart.
void reestablish_distribution_info_for_restart(OomphCommunicator *comm_pt, std::istream &restart_file)
cstr elem_len * i
Definition: cfortran.h:607
virtual void reestablish_distribution_info_for_restart(OomphCommunicator *comm_pt, std::istream &restart_file)
std::map< unsigned, Vector< TriangleMeshPolyLine * > > Boundary_subpolylines
void disable_internal_boundary_refinement()
Helper function for disabling the use of boundary refinement.
unsigned Counter_for_flat_packed_doubles
Counter used when processing vector of flat-packed doubles.
std::map< unsigned, Vector< Node * > > External_haloed_node_pt
Map of vectors holding the pointers to the external haloed nodes.
Definition: mesh.h:149
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)
const double Pi
50 digits from maple
Vector< TriangleMeshOpenCurve * > & internal_open_curves_pt()
Helper function for getting access to the internal open boundaries.
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 ...
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1305
void operator=(const TriangleMesh &)
Broken assignment operator.
unsigned try_to_add_root_haloed_element_pt(const unsigned &p, GeneralisedElement *&el_pt)
void add_shared_boundary_node(const unsigned &b, Node *node_pt)
void add_region_coordinates(const unsigned &i, Vector< double > &region_coordinates)
Helper function for getting the extra regions.
A general Finite Element class.
Definition: elements.h:1271
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.
Vector< TriangleMeshPolygon * > Internal_polygon_pt
Vector of polygons that define internal polygons.
char t
Definition: cfortran.h:572
void enable_use_attributes()
Helper function for enabling the use of attributes.
void get_element_edges_on_boundary(std::map< std::pair< Node *, Node * >, unsigned > &element_edges_on_boundary)
Helper object for dealing with the parameters used for the TriangleMesh objects.
std::map< unsigned, double > Regions_areas
Target areas for regions; defaults to 0.0 which (luckily) implies "no specific target area" for trian...
const unsigned shared_boundaries_ids(const unsigned &p, const unsigned &q, const unsigned &i) const
void disable_boundary_refinement()
Helper function for disabling the use of boundary refinement.
void update_triangulateio()
Update the triangulateio object to the current nodal positions.
bool Triangulateio_exists
Boolean defining if Triangulateio object has been built or not.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void get_required_nodal_information_helper(int &iproc, Node *nod_pt, Mesh *const &mesh_pt, int &n_cont_inter_values, Vector< unsigned > &send_unsigneds, Vector< double > &send_doubles)
Helper function to get the required nodal information from an external haloed node so that a fully-fu...
void disable_print_timings_adaptation()
Disables printing of timings for adaptation.
std::map< unsigned, Vector< unsigned > > & shared_boundary_from_processors()
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.
std::map< unsigned, Vector< unsigned > > Shared_boundary_from_processors
void update_holes_information_helper(Vector< TriangleMeshPolygon * > &polygons_pt, Vector< Vector< double > > &output_holes_coordinates)
Vector< unsigned > Flat_packed_unsigneds
Vector of flat-packed unsigneds to be communicated with other processors.
void refine_uniformly(DocInfo &doc_info)
Refine mesh uniformly and doc process.
coord2_t cross(const Point &O, const Point &A, const Point &B)
OomphInfo oomph_info
Vector< Vector< double > > & extra_holes_coordinates()
Helper function for getting access to the extra holes.
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 ...
unsigned & nbin_x_for_projection()
Read/write access to number of bins in the x-direction when projecting old solution onto new mesh...
std::set< unsigned > Internal_boundaries
e
Definition: cfortran.h:575
std::map< unsigned, Vector< Vector< Node * > > > Boundary_segment_node_pt
Used to store the nodes associated to a boundary and to an specific segment (this only applies in dis...
void output_boundary_coordinates(const unsigned &b, std::ostream &outfile)
std::map< unsigned, Vector< double > > Regions_coordinates
Store the coordinates for defining extra regions The key on the map is the region id...
TriangleMesh()
Empty constructor.
std::set< TriangleMeshCurveSection * > Free_curve_section_pt
A set that contains the curve sections created by this object therefore it is necessary to free their...
void disable_projection()
Disables the solution projection step during adaptation.
unsigned Nbin_y_for_area_transfer
Number of bins in the y-direction when transferring target areas by bin method.
std::map< unsigned, unsigned > & shared_boundary_overlaps_internal_boundary()
TriangleMeshParameters(TriangleMeshClosedCurve *outer_boundary_pt)
std::set< TriangleMeshPolygon * > Free_polygon_pt
A set that contains the polygons created by this object therefore it is necessary to free their assoc...
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region.
Vector< TriangleMeshClosedCurve * > & internal_closed_curve_pt()
Helper function for getting access to the internal closed boundaries.
void remesh_from_internal_triangulateio()
Completely regenerate the mesh from the trianglateio structure.
void flush_shared_boundary_node(const unsigned &b)
Vector< Vector< Node * > > & boundary_segment_node_pt(const unsigned &b)
Return direct access to nodes associated with a boundary but sorted in segments.
void set_print_level_timings_adaptation(const unsigned &print_level)
Sets the printing level of timings for adaptation.
Vector< TriangleMeshClosedCurve * > internal_closed_curve_pt() const
Helper function for getting the internal closed boundaries.
const unsigned initial_shared_boundary_id()
const unsigned nshared_boundary_polyline(const unsigned &p, const unsigned &c) const
void dump_distributed_info_for_restart(std::ostream &dump_file)
Used to dump info. related with distributed triangle meshes.
void set_communicator_pt(OomphCommunicator *comm_pt)
Function to set communicator (mesh is then assumed to be distributed)
Vector< TriangleMeshClosedCurve * > outer_boundary_pt() const
Helper function for getting the outer boundary.
double & max_element_size()
Max element size allowed during adaptation.
bool Internal_boundary_refinement
Do not allow refinement of nodes on the internal boundary.
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
virtual ~TriangleMeshParameters()
Empty destructor.
void create_triangulateio_from_polyfiles(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TriangulateIO &triangle_io, bool &use_attributes)
Vector< std::map< unsigned, Vector< FiniteElement * > > > Boundary_region_element_pt
Storage for elements adjacent to a boundary in a particular region.
std::map< unsigned, Vector< Node * > > Shared_boundary_node_pt
virtual ~TriangleMesh()
Destructor.
Vector< double > Flat_packed_doubles
Vector of flat-packed doubles to be communicated with other processors.
std::vector< Point > convex_hull(std::vector< Point > P)
TriangleMeshOpenCurve * create_open_curve_with_polyline_helper(TriangleMeshOpenCurve *open_curve_pt, unsigned &max_bnd_id_local)
Helper function that creates and returns an open curve with the polyline representation of its consti...
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...
TriangleMeshPolyLine * shared_boundary_polyline_pt(const unsigned &p, const unsigned &c, const unsigned &i) 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)
double element_area() const
Helper function for getting the element area.
unsigned nsorted_shared_boundary_node(unsigned &b)
const unsigned nshared_boundaries() const
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
OomphCommunicator * communicator_pt() const
Definition: mesh.h:1317
unsigned Nbin_x_for_projection
Number of bins in the x-direction when projecting the solution from the old mesh into the new mesh...
double size() const
Definition: elements.cc:4121
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)
Vector< Vector< Vector< unsigned > > > shared_boundaries_ids() const
void synchronize_boundary_coordinates(const unsigned &b)
In charge of sinchronize the boundary coordinates for internal boundaries that were split as part of ...
const bool shared_boundary_overlaps_internal_boundary(const unsigned &shd_bnd_id)
Vector< Vector< double > > extra_holes_coordinates() const
Helper function for getting the extra holes.
void enable_shared_boundary_unrefinement_constrained_by_target_areas()
void get_shared_boundaries_overlapping_internal_boundary(const unsigned &internal_bnd_id, Vector< unsigned > &shd_bnd_ids)
Class defining a polyline for use in Triangle Mesh generation.
std::map< unsigned, Vector< GeneralisedElement * > > Root_haloed_element_pt
Map of vectors holding the pointers to the root haloed elements.
Definition: mesh.h:114
bool Boundary_refinement
Do not allow refinement of nodes on the boundary.
std::map< unsigned, Vector< GeneralisedElement * > > External_haloed_element_pt
Map of vectors holding the pointers to the external haloed elements.
Definition: mesh.h:140
Vector< unsigned > Oomph_vertex_nodes_id
Vector storing oomph-lib node number for all vertex nodes in the TriangulateIO representation of the ...
unsigned & nbin_y_for_projection()
Read/write access to number of bins in the y-direction when projecting old solution onto new mesh...
Vector< Vector< int > > Face_index_at_boundary
For the e-th finite element on boundary b, this is the index of the face that lies along that boundar...
Definition: mesh.h:106
double Element_area
The element are when calling triangulate external routine.
std::map< unsigned, bool > Boundary_was_splitted
TriangleMeshClosedCurve * outer_boundary_pt(const unsigned &i) const
Helper function for getting the i-th outer boundary.
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.
Vector< Vector< double > > Extra_holes_coordinates
Storage for extra coordinates for holes.
void add_face_index_at_shared_boundary(const unsigned &b, const unsigned &i)
void clear_triangulateio(TriangulateIO &triangulate_io, const bool &clear_hole_data)
Clear TriangulateIO structure.
void add_shared_boundary_element(const unsigned &b, FiniteElement *ele_pt)
void create_distributed_domain_representation(Vector< TriangleMeshPolygon * > &polygons_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt)
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...
MeshUpdateFctPt Mesh_update_fct_pt
Function pointer to function that updates the mesh following the snapping of boundary nodes to the bo...
const unsigned read_unsigned_line_helper(std::istream &read_file)
bool Print_timings_projection
Enable/disable printing timings for projection.
double Min_element_size
Min permitted element size.
FiniteElement * shared_boundary_element_pt(const unsigned &b, const unsigned &e)
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
static char t char * s
Definition: cfortran.h:572
void enable_timings_tranfering_target_areas()
Enables info. and timings for tranferring of target areas.
unsigned unrefine_uniformly()
Unrefine mesh uniformly: Return 0 for success, 1 for failure (if unrefinement has reached the coarses...
TriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
std::map< unsigned, std::vector< bool > > Boundary_marked_as_shared_boundary
bool Allow_automatic_creation_of_vertices_on_boundaries
Flag to indicate whether the automatic creation of vertices along the boundaries by Triangle is allow...
Vector< unsigned > & shared_boundary_from_processors(const unsigned &b)
void enable_print_timings_adaptation(const unsigned &print_level=1)
Enables printing of timings for adaptation.
Vector< unsigned > shared_boundaries_ids(const unsigned &p, const unsigned &q) const
Node * shared_boundary_node_pt(const unsigned &b, const unsigned &n)
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 ...
double & element_area()
Helper function for getting access to the element area.
void disable_print_timings_load_balance()
Disables printing of timings for load balance.
std::map< unsigned, Vector< std::pair< double, double > > > Polygonal_vertex_arclength_info
Storage for pairs of doubles representing: .first: the arclength along the polygonal representation o...
void disable_use_attributes()
Helper function for disabling the use of attributes.
void compute_holes_left_by_halo_elements_helper(Vector< Vector< double > > &output_holes_coordinates)
const unsigned nboundary_subpolylines(const unsigned &b)
A class that represents a collection of data; each Data object may contain many different individual ...
Definition: nodes.h:89
Vector< TriangleMeshClosedCurve * > Outer_boundary_pt
The outer boundary.
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:197
virtual ~RefineableSolidTriangleMesh()
Empty Destructor.
std::map< unsigned, std::set< Vector< double > > > Boundary_connections_pt
TriangleScaffoldMesh * Tmp_mesh_pt
Temporary scaffold mesh.
void triangulate(char *triswitches, struct oomph::TriangulateIO *in, struct oomph::TriangulateIO *out, struct oomph::TriangulateIO *vorout)
TimeStepper * Time_stepper_pt
Timestepper used to build elements.
bool is_node_on_shared_boundary(const unsigned &b, Node *const &node_pt)
void snap_nodes_onto_geometric_objects()
Snap the boundary nodes onto any curvilinear boundaries defined by geometric objects.
const unsigned nshared_boundary_overlaps_internal_boundary()
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
void enable_boundary_refinement()
Helper function for enabling the use of boundary refinement.
Base class defining a closed curve for the Triangle mesh generation.
RefineableTriangleMesh(TriangleMeshParameters &triangle_mesh_parameters, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Build mesh, based on the specifications on TriangleMeshParameters.
void doc_adaptivity_targets(std::ostream &outfile)
Doc the targets for mesh adaptation.
void initialise_triangulateio(TriangulateIO &triangle_io)
Initialise TriangulateIO structure.
std::map< unsigned, Vector< Node * > > Halo_node_pt
Map of vectors holding the pointers to the halo nodes.
Definition: mesh.h:117
virtual void reset_boundary_element_info(Vector< unsigned > &ntmp_boundary_elements, Vector< Vector< unsigned > > &ntmp_boundary_elements_in_region, Vector< FiniteElement * > &deleted_elements)
Reset the boundary elements info. after load balance have taken place.
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
bool is_internal_boundary_refinement_allowed() const
Helper function for getting the status of boundary refinement.
TriangleMeshParameters(Vector< TriangleMeshClosedCurve * > &outer_boundary_pt)
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
Vector< Vector< FiniteElement * > > Boundary_element_pt
Vector of Vector of pointers to elements on the boundaries: Boundary_element_pt(b,e)
Definition: mesh.h:102
std::map< unsigned, Vector< Node * > > Sorted_shared_boundary_node_pt
OomphCommunicator * communicator_pt() const
Read-only access fct to communicator (Null if mesh is not distributed)
double Min_permitted_angle
Min angle before remesh gets triggered.
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.
unsigned Print_timings_level_adaptation
The printing level for adaptation.
std::map< unsigned, Vector< Node * > > External_halo_node_pt
Map of vectors holding the pointers to the external halo nodes.
Definition: mesh.h:146
General SolidMesh class.
Definition: mesh.h:2252
Vector< std::string > Flat_packed_unsigneds_string
virtual ~SolidTriangleMesh()
Empty Destructor.
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...
unsigned Nbin_y_for_projection
Number of bins in the y-direction when projecting the solution from the old mesh into the new mesh...
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, these elements are then used on the function that computes the shared boundaries among the processors.
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)
std::map< unsigned, double > & target_area_for_region()
Helper function for getting access to the region's target areas.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void update_triangulateio(Vector< Vector< double > > &internal_point)
Update the TriangulateIO object to the current nodal position and the centre hole coordinates...
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)
std::set< TriangleMeshOpenCurve * > Free_open_curve_pt
A set that contains the open curves created by this object therefore it is necessary to free their as...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void disable_timings_tranfering_target_areas()
Disables info. and timings for tranferring of target areas.
Unstructured Triangle Mesh upgraded to solid mesh.
bool is_boundary_refinement_allowed() const
Helper function for getting the status of boundary refinement.
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 flush_shared_boundary_element(const unsigned &b)
Vector< Vector< double > > Extra_holes_coordinates
Store the coordinates for defining extra holes.
const unsigned nshared_boundary_element(const unsigned &b)
Vector< TriangleMeshOpenCurve * > internal_open_curves_pt() const
Helper function for getting the internal open boundaries.
bool Use_attributes
Define the use of attributes (regions)
const bool boundary_was_splitted(const unsigned &b)
const unsigned shared_boundary_overlapping_internal_boundary(const unsigned &shd_bnd_id)
const unsigned nshared_boundary_curves(const unsigned &p) const
unsigned try_to_add_node_pt_load_balance(Vector< Node * > &new_nodes_on_domain, Node *&node_pt)
RefineableSolidTriangleMesh(const Vector< double > &target_area, TriangulateIO &triangulate_io, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper, const bool &use_attributes=false, const bool &allow_automatic_creation_of_vertices_on_boundaries=true, OomphCommunicator *comm_pt=0)
Build mesh from specified triangulation and associated target areas for elements in it...
void enable_print_timings_load_balance(const unsigned &print_level=1)
Enables printing of timings for load balance.
void enable_internal_boundary_refinement()
Helper function for enabling the use of boundary refinement.
void set_print_level_timings_load_balance(const unsigned &print_level)
Sets the printing level of timings for load balance.
Vector< Vector< double > > Original_extra_holes_coordinates
const bool boundary_marked_as_shared_boundary(const unsigned &b, const unsigned &isub)
std::map< unsigned, Vector< double > > Regions_coordinates
Node * sorted_shared_boundary_node_pt(unsigned &b, unsigned &i)
Vector< Vector< Vector< unsigned > > > Shared_boundaries_ids
unsigned Counter_for_flat_packed_unsigneds
Counter used when processing vector of flat-packed unsigneds.
void sort_polylines_helper(Vector< TriangleMeshPolyLine * > &unsorted_polylines_pt, Vector< Vector< TriangleMeshPolyLine * > > &sorted_polylines_pt)
Vector< Node * > sorted_shared_boundary_node_pt(unsigned &b)
Vector< Vector< Vector< unsigned > > > & shared_boundaries_ids()
unsigned Print_timings_level_load_balance
The printing level for load balance.
Vector< TriangleMeshPolyLine * > & shared_boundary_polyline_pt(const unsigned &p, const unsigned &c)
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
TriangleMesh(const TriangleMesh &dummy)
Broken copy constructor.
Vector< Vector< unsigned > > & shared_boundaries_ids(const unsigned &p)
const unsigned final_shared_boundary_id()
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...
int face_index_at_shared_boundary(const unsigned &b, const unsigned &e)
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.
void flush_element_and_node_storage()
Flush storage for elements and nodes by emptying the vectors that store the pointers to them...
Definition: mesh.h:427
unsigned Nbin_x_for_area_transfer
Number of bins in the x-direction when transferring target areas by bin method.
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...
void create_tmp_polygons_helper(Vector< Vector< TriangleMeshPolyLine * > > &polylines_pt, Vector< TriangleMeshPolygon * > &polygons_pt)
void read_distributed_info_for_restart(std::istream &restart_file)
Used to read info. related with distributed triangle meshes.
double & min_permitted_angle()
Min angle before remesh gets triggered.
Base class for time-stepping schemes. Timestepper provides an approximation of the temporal derivativ...
Definition: timesteppers.h:219
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 set_lagrangian_nodal_coordinates()
Make the current configuration the undeformed one by setting the nodal Lagrangian coordinates to thei...
Definition: mesh.cc:9196
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...
Unstructured refineable Triangle Mesh.
std::map< unsigned, Vector< FiniteElement * > > Shared_boundary_element_pt
bool triangulateio_exists()
Boolean defining if Triangulateio object has been built or not.
Vector< unsigned > oomph_vertex_nodes_id()
Return the vector that contains the oomph-lib node number for all vertex nodes in the TriangulateIO r...
TriangleMeshClosedCurve *& outer_boundary_pt(const unsigned &i)
Helper function for getting access to the i-th outer boundary.
unsigned & nbin_x_for_area_transfer()
Read/write access to number of bins in the x-direction when transferring target areas by bin method...
Vector< TriangleMeshPolyLine * > & boundary_subpolylines(const unsigned &b)
A general mesh class.
Definition: mesh.h:74
Vector< TriangleMeshPolygon * > Outer_boundary_pt
Polygon that defines outer boundaries.
unsigned try_to_add_haloed_node_pt(const unsigned &p, Node *&nod_pt)
void shared_boundaries_in_this_processor(Vector< unsigned > &shared_boundaries_in_this_processor)
bool is_use_attributes() const
Helper function for getting the status of use_attributes variable.
void remove_boundary_nodes()
Clear all pointers to boundary nodes.
Definition: mesh.cc:208
Vector< Vector< Vector< TriangleMeshPolyLine * > > > Shared_boundary_polyline_pt
void build_triangulateio(Vector< TriangleMeshPolygon * > &outer_polygons_pt, Vector< TriangleMeshPolygon * > &internal_polygons_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt, Vector< Vector< double > > &extra_holes_coordinates, std::map< unsigned, Vector< double > > &regions_coordinates, std::map< unsigned, double > &regions_areas, TriangulateIO &triangulate_io)
Create TriangulateIO object from outer boundaries, internal boundaries, and open curves. Add the holes and regions information as well.
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)
unsigned try_to_add_element_pt_load_balance(Vector< FiniteElement * > &new_elements_on_domain, FiniteElement *&ele_pt)
unsigned & nbin_y_for_area_transfer()
Read/write access to number of bins in the y-direction when transferring target areas by bin method...
An oomph-lib wrapper to the MPI_Comm communicator object. Just contains an MPI_Comm object (which is ...
Definition: communicator.h:57
Class defining a closed polygon for the Triangle mesh generation.