unstructured_two_d_mesh_geometry_base.h
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 // Contains the definition of a TriangulateIO object. This is used to
27 // define the complex geometry of a two-dimensional mesh which is why
28 // it resides here. The definition of things like TriangleMeshPolygons
29 // and other classes which define geometrical aspects of a 2D mesh can
30 // also be found here. The class UnstructuredTwoDMeshGeometryBase is
31 // defined here. It forms the base class for QuadFromTriangleMesh and
32 // TriangleMeshBase. This class makes use of previously written code
33 // to create TriangleScaffoldMeshes and avoid a large amount of code
34 // duplication.
35 #ifndef OOMPH_UNSTRUCTURED_TWO_D_MESH_GEOMETRY_BASE_HEADER
36 #define OOMPH_UNSTRUCTURED_TWO_D_MESH_GEOMETRY_BASE_HEADER
37 
38 // The scaffold mesh
39 #include "mesh.h"
40 
41 /// /////////////////////////////////////////////////////////////////
42 /// /////////////////////////////////////////////////////////////////
43 /// /////////////////////////////////////////////////////////////////
44 
45 
46 namespace oomph
47 {
48 #ifdef OOMPH_HAS_TRIANGLE_LIB
49 
50  //=====================================================================
51  /// The Triangle data structure, modified from the triangle.h header
52  /// supplied with triangle 1.6. by J. R. Schewchuk. We need to define
53  /// this here separately because we can't include a c header directly
54  /// into C++ code!
55  //=====================================================================
57  {
58  /// Pointer to list of points x coordinate followed by y coordinate
59  double* pointlist;
60 
61  /// Pointer to list of point attributes
63 
64  /// Pointer to list of point markers
68 
76 
80 
81  double* holelist;
83 
84  double* regionlist;
86 
87  int* edgelist;
88  int* edgemarkerlist; // <---- contains boundary ID (offset by one)
89  double* normlist;
91  };
92 
93 
94  /// /////////////////////////////////////////////////////////////////
95  /// /////////////////////////////////////////////////////////////////
96  /// /////////////////////////////////////////////////////////////////
97 
98 
99  //==================================================================
100  /// Helper namespace for triangle meshes
101  //==================================================================
102  namespace TriangleHelper
103  {
104  /// Clear TriangulateIO structure
105  extern void clear_triangulateio(TriangulateIO& triangulate_io,
106  const bool& clear_hole_data = true);
107 
108  /// Initialise TriangulateIO structure
109  extern void initialise_triangulateio(TriangulateIO& triangle_io);
110 
111  /// Make (partial) deep copy of TriangulateIO object. We only copy
112  /// those items we need within oomph-lib's adaptation procedures.
113  /// Warnings are issued if triangulate_io contains data that is not
114  /// not copied, unless quiet=true;
116  TriangulateIO& triangle_io, const bool& quiet);
117 
118  /// Write the triangulateio data to disk as a poly file,
119  /// mainly used for debugging
120  extern void write_triangulateio_to_polyfile(TriangulateIO& triangle_io,
121  std::ostream& poly_file);
122 
123  /// Create a triangulateio data file from ele node and poly
124  /// files. This is used if the mesh is generated by using Triangle
125  /// externally. The triangulateio structure is required to dump the mesh
126  /// topology for restarts.
128  const std::string& node_file_name,
129  const std::string& element_file_name,
130  const std::string& poly_file_name,
131  TriangulateIO& triangle_io,
132  bool& use_attributes);
133 
134  /// Write all the triangulateio data to disk in a dump file
135  /// that can then be used to restart simulations
136  extern void dump_triangulateio(TriangulateIO& triangle_io,
137  std::ostream& dump_file);
138 
139  /// Read the triangulateio data from a dump file on
140  /// disk, which can then be used to restart simulations
141  extern void read_triangulateio(std::istream& restart_file,
142  TriangulateIO& triangle_io);
143 
144  } // namespace TriangleHelper
145 
146 #endif
147 
148 
149  /// //////////////////////////////////////////////////////////////////////
150  /// //////////////////////////////////////////////////////////////////////
151  /// //////////////////////////////////////////////////////////////////////
152 
153 
154  class TriangleMeshPolyLine;
155  class TriangleMeshCurviLine;
156 
157  //=====================================================================
158  /// Base class for defining a triangle mesh boundary, this class has the
159  /// methods that allow to connect the initial and final ends to other
160  /// triangle mesh boundaries
161  //=====================================================================
163  {
164  public:
165  /// Empty constructor. Initialises the curve section as non connected
167  : Initial_vertex_connected(false),
168  Final_vertex_connected(false),
173  Refinement_tolerance(0.08),
175  Maximum_length(-1.0)
176  {
177  }
178 
179  /// Empty destructor
181 
182  /// Number of segments that this part of the
183  /// boundary is to be represented by. This corresponds
184  /// to the number of straight-line segments in triangle
185  /// representation.
186  virtual unsigned nsegment() const = 0;
187 
188  /// Boundary id
189  virtual unsigned boundary_id() const = 0;
190 
191  /// Boundary chunk (Used when a boundary is represented by more
192  /// than one polyline
193  virtual unsigned boundary_chunk() const = 0;
194 
195  /// Number of vertices
196  virtual unsigned nvertex() const = 0;
197 
198  /// Output the curve_section
199  virtual void output(std::ostream& outfile,
200  const unsigned& n_sample = 50) = 0;
201 
202  /// Enable refinement of curve section to create a better
203  /// representation of curvilinear boundaries (e.g. in free-surface
204  /// problems). See tutorial for
205  /// interpretation of the optional argument which specifies the
206  /// refinement tolerance. It defaults to 0.08 and the smaller the
207  /// number the finer the surface representation.
208  void enable_refinement_tolerance(const double& tolerance = 0.08)
209  {
210  Refinement_tolerance = tolerance;
211  }
212 
213  /// Set tolerance for refinement of curve sections to create a better
214  /// representation of curvilinear boundaries (e.g. in free-surface
215  /// problems). See tutorial for
216  /// interpretation of the refinement tolerance. (The smaller the
217  /// number the finer the surface representation). If set to
218  /// a negative value, we're switching off refinement --
219  /// equivalent to calling disable_polyline_refinement()
220  void set_refinement_tolerance(const double& tolerance)
221  {
222  Refinement_tolerance = tolerance;
223  }
224 
225  /// Get tolerance for refinement of curve sections to create a better
226  /// representation of curvilinear boundaries (e.g. in free-surface
227  /// problems). See tutorial for
228  /// interpretation. If it's negative refinement is disabled.
230  {
231  return Refinement_tolerance;
232  }
233 
234  /// Disable refinement of curve section
236  {
237  Refinement_tolerance = -1.0;
238  }
239 
240  /// Enable unrefinement of curve sections to avoid unnecessarily
241  /// large numbers of elements on of curvilinear boundaries (e.g. in
242  /// free-surface problems). See tutorial for interpretation of the optional
243  /// argument which specifies the unrefinement tolerance. It defaults to 0.04
244  /// and the larger the number the more agressive we are when removing
245  /// unnecessary vertices on gently curved polylines.
246  void enable_unrefinement_tolerance(const double& tolerance = 0.04)
247  {
248  Unrefinement_tolerance = tolerance;
249  }
250 
251  /// Set tolerance for unrefinement of curve sections
252  /// to avoid unnecessarily large
253  /// numbers of elements on of curvilinear boundaries (e.g. in free-surface
254  /// problems). See tutorial for
255  /// interpretation of the optional argument which specifies the
256  /// unrefinement tolerance. It defaults to 0.04 and the larger the number
257  /// the more agressive we are when removing unnecessary vertices on
258  /// gently curved polylines. If set to
259  /// a negative value, we're switching off unrefinement --
260  /// equivalent to calling disable_curve_section_unrefinement()
261  void set_unrefinement_tolerance(const double& tolerance)
262  {
263  Unrefinement_tolerance = tolerance;
264  }
265 
266  /// Get tolerance for unrefinement of curve section to create a
267  /// better representation of curvilinear boundaries (e.g. in free-surface
268  /// problems). See tutorial for
269  /// interpretation. If it's negative unrefinement is disabled.
271  {
272  return Unrefinement_tolerance;
273  }
274 
275  /// Disable unrefinement of curve sections
277  {
278  Unrefinement_tolerance = -1.0;
279  }
280 
281  /// Allows to specify the maximum distance between two vertices
282  /// that define the associated polyline of the curve section, it only
283  /// takes effect on the unrefinement and refinement steps
284  void set_maximum_length(const double& maximum_length)
285  {
287  }
288 
289  /// Disables the use of the maximum length criteria on the
290  /// unrefinement or refinement steps
292  {
293  Maximum_length = -1.0;
294  }
295 
296  /// Gets access to the maximum length variable
297  double maximum_length()
298  {
299  return Maximum_length;
300  }
301 
302  /// Get first vertex coordinates
303  virtual void initial_vertex_coordinate(Vector<double>& vertex) = 0;
304 
305  /// Get last vertex coordinates
306  virtual void final_vertex_coordinate(Vector<double>& vertex) = 0;
307 
308  /// Connects the initial vertex of the curve section to a desired
309  /// target polyline by specifying the vertex number. There is a checking
310  /// which verifies that the initial vertex is close enough to the
311  /// destination vertex on the target polyline by no more than the specified
312  /// tolerance
314  TriangleMeshPolyLine* polyline_pt,
315  const unsigned& vertex_number,
316  const double& tolerance_for_connection = 1.0e-14);
317 
318  /// Connects the final vertex of the curve section to a desired
319  /// target polyline by specifying the vertex number. There is a checking
320  /// which verifies that the final vertex is close enough to the
321  /// destination vertex on the target polyline by no more than the specified
322  /// tolerance
324  TriangleMeshPolyLine* polyline_pt,
325  const unsigned& vertex_number,
326  const double& tolerance_for_connection = 1.0e-14);
327 
328  /// Connects the initial vertex of the curve section to a desired
329  /// target curviline by specifying the s value (intrinsic value on the
330  /// geometric object of the curviline) where to connect on the target
331  /// curviline. There is a checking which verifies that the initial vertex
332  /// and the coordinates on the given s value are close enough by no more
333  /// than the given tolerance
335  TriangleMeshCurviLine* curviline_pt,
336  const double& s_value,
337  const double& tolerance_for_connection = 1.0e-14);
338 
339  /// Connects the final vertex of the curve section to a desired
340  /// target curviline by specifying the s value (intrinsic value on the
341  /// geometric object of the curviline) where to connect on the target
342  /// curviline. There is a checking which verifies that the final vertex
343  /// and the coordinates on the given s value are close enough by no more
344  /// than the given tolerance
346  TriangleMeshCurviLine* curviline_pt,
347  const double& s_value,
348  const double& tolerance_for_connection = 1.0e-14);
349 
350  /// Test whether initial vertex is connected or not
352  {
354  }
355 
356  /// Sets the initial vertex as connected
358  {
360  }
361 
362  /// Sets the initial vertex as non connected
364  {
365  Initial_vertex_connected = false;
366  }
367 
368  /// Set the initial vertex connection as suspended, it will be
369  /// resumed when the method to resume the connections is called
370  /// This method is only used in a distributed context, when the
371  /// boundary to connect is no longer part of the domain in the
372  /// processor
374  {
376  {
377  Initial_vertex_connected = false;
379  }
380  }
381 
382  /// Resumes the initial vertex connection, it may be that after load
383  /// balancing the boundary to which the connection was suspended be part
384  /// of the domain
386  {
388  {
391  }
392  }
393 
394  /// Test whether final vertex is connected or not
396  {
397  return Final_vertex_connected;
398  }
399 
400  /// Sets the final vertex as connected
402  {
403  Final_vertex_connected = true;
404  }
405 
406  /// Sets the final vertex as non connected
408  {
409  Final_vertex_connected = false;
410  }
411 
412  /// Set the final vertex connection as suspended, it will be
413  /// resumed when the method to resume the connections is called
414  /// This method is only used in a distributed context, when the
415  /// boundary to connect is no longer part of the domain in the
416  /// processor
418  {
420  {
421  Final_vertex_connected = false;
423  }
424  }
425 
426  /// Resumes the final vertex connection, it may be that after load
427  /// balancing the boundary to which the connection was suspended be part
428  /// of the domain
430  {
432  {
433  Final_vertex_connected = true;
435  }
436  }
437 
438  /// Gets the id to which the initial end is connected
440  {
442  }
443 
444  /// Sets the id to which the initial end is connected
446  {
448  }
449 
450  /// Gets the vertex number to which the initial end is connected
452  {
454  }
455 
456  /// Sets the vertex number to which the initial end is connected
458  {
460  }
461 
462  /// Gets the boundary chunk to which the initial end is connected
464  {
466  }
467 
468  /// Sets the boundary chunk to which the initial end is connected
470  {
472  }
473 
474  /// Gets the id to which the final end is connected
476  {
478  }
479 
480  /// Sets the id to which the final end is connected
482  {
484  }
485 
486  /// Sets the vertex number to which the final end is connected
488  {
490  }
491 
492  /// Gets the vertex number to which the final end is connected
494  {
496  }
497 
498  /// Gets the boundary chunk to which the final end is connected
500  {
502  }
503 
504  /// Sets the boundary chunk to which the final end is connected
506  {
508  }
509 
510  /// Test whether the initial vertex is connected to a curviline
512  {
514  }
515 
516  /// Sets the initial vertex as connected to a curviline
518  {
520  }
521 
522  /// Sets the initial vertex as non connected to a curviline
524  {
526  }
527 
528  /// Test whether the final vertex is connected to a curviline
530  {
532  }
533 
534  /// Sets the final vertex as connected to a curviline
536  {
538  }
539 
540  /// Sets the final vertex as non connected to a curviline
542  {
544  }
545 
546  /// Gets the s value to which the initial end is connected
548  {
550  }
551 
552  /// Sets the s value to which the initial end is connected
554  {
556  }
557 
558  /// Gets the s value to which the final end is connected
560  {
562  }
563 
564  /// Sets the s value to which the final end is connected
566  {
568  }
569 
570  /// Gets the tolerance value for connections among
571  /// curvilines
573  {
575  }
576 
577  /// Sets the tolerance value for connections among
578  /// curvilines
580  {
582  }
583 
584  protected:
585  /// Used for stating if the initial end is connected
586  /// to another boundary
588 
589  /// Used for stating if the final end is connected
590  /// to another boundary
592 
593  /// Indicates if the connection is suspended because the
594  /// boundary to connect is no longer part of the domain (only used in
595  /// a distributed context)
597 
598  /// Indicates if the connection is suspended because the
599  /// boundary to connect is no longer part of the domain (only used in
600  /// a distributed context)
602 
603  /// Stores the id to which the initial end is connected
605 
606  /// Stores the vertex number used for connection with
607  /// the initial end
609 
610  /// Stores the chunk number of the boundary to which is
611  /// connected th initial end
613 
614  /// Stores the id to which the initial end is connected
616 
617  /// Stores the vertex number used for connection with
618  /// the final end
620 
621  /// Stores the chunk number of the boundary to which is
622  /// connected th initial end
624 
625  /// States if the initial vertex is connected to a curviline
627 
628  /// States if the final vertex is connected to a curviline
630 
631  /// Stores the s value used for connecting the
632  /// initial end with a curviline
634 
635  /// Stores the s value used for connecting the
636  /// final end with a curviline
638 
639  /// Tolerance used for connecting the ends to a curviline
641 
642  private:
643  /// Tolerance for refinement of curve sections (neg if refinement is
644  /// disabled)
646 
647  /// Tolerance for unrefinement of curve sections (neg if refinement
648  /// is disabled)
650 
651  /// Maximum allowed distance between two vertices on the polyline
652  /// representation of the curve section
654  };
655 
656 
657  //=====================================================================
658  /// Class definining a curvilinear triangle mesh boundary in terms
659  /// of a GeomObject. Curvlinear equivalent of PolyLine.
660  //=====================================================================
662  {
663  public:
664  /// Constructor: Specify GeomObject, the start and end coordinates
665  /// of the relevant boundary in terms of the GeomObject's intrinsic
666  /// coordinate, the number of (initially straight-line) segments that
667  /// this GeomObject is to be split up into, and the boundary ID.
668  /// The final optional boolean argument specifies if vertices in
669  /// polygonhal represenation are spaced
670  /// (approximately) evenly in arclength along the GeomObject [true,
671  /// default] or in equal increments in zeta.
672  /// This is the curvlinear equivalent of PolyLine.
674  const double& zeta_start,
675  const double& zeta_end,
676  const unsigned& nsegment,
677  const unsigned& boundary_id,
678  const bool& space_vertices_evenly_in_arclength = true,
679  const unsigned& boundary_chunk = 0)
688  {
689  }
690 
691 
692  /// Empty Destuctor
694 
695  /// Pointer to GeomObject that represents this part of the boundary
697  {
698  return Geom_object_pt;
699  }
700 
701  /// Start coordinate in terms of the GeomObject's intrinsic coordinate
702  double zeta_start()
703  {
704  return Zeta_start;
705  }
706 
707  /// End coordinate in terms of the GeomObject's intrinsic coordinate
708  double zeta_end()
709  {
710  return Zeta_end;
711  }
712 
713  /// Number of (initially straight-line) segments that this part of
714  /// the boundary is to be represented by
715  unsigned nsegment() const
716  {
717  return Nsegment;
718  }
719 
720  /// Number of (initially straight-line) segments that this part of
721  /// the boundary is to be represented by. This version allows the change of
722  /// the number of segments
723  unsigned& nsegment()
724  {
725  return Nsegment;
726  }
727 
728  /// Boundary ID
729  unsigned boundary_id() const
730  {
731  return Boundary_id;
732  }
733 
734  /// Boundary chunk (Used when a boundary is represented by more
735  /// than one polyline
736  unsigned boundary_chunk() const
737  {
738  return Boundary_chunk;
739  }
740 
741  /// Output curvilinear boundary at n_sample (default: 50) points
742  void output(std::ostream& outfile, const unsigned& n_sample = 50)
743  {
744  outfile << "ZONE T=\"Boundary " << Boundary_id << "\"\n";
745  Vector<double> zeta(1);
746  Vector<double> r(2);
747  for (unsigned i = 0; i < n_sample; i++)
748  {
749  zeta[0] = Zeta_start +
750  (Zeta_end - Zeta_start) * double(i) / double(n_sample - 1);
751  Geom_object_pt->position(zeta, r);
752  outfile << r[0] << " " << r[1] << std::endl;
753  }
754  }
755 
756  /// Boolean to indicate if vertices in polygonal representation
757  /// of the Curvline are spaced (approximately) evenly in arclength
758  /// along the GeomObject [true] or in equal increments in zeta [false]
760  {
762  }
763 
764  /// Number of vertices
765  unsigned nvertex() const
766  {
767  return 2;
768  }
769 
770  /// Get first vertex coordinates
772  {
773  Vector<double> z(1);
774  z[0] = Zeta_start;
775  Geom_object_pt->position(z, vertex);
776  }
777 
778  /// Get last vertex coordinates
780  {
781  Vector<double> z(1);
782  z[0] = Zeta_end;
783  Geom_object_pt->position(z, vertex);
784  }
785 
786  /// Does the vector for storing connections has elements?
788  {
789  return !Connection_points_pt.empty();
790  }
791 
792  /// Returns the connection points vector
794  {
795  return &Connection_points_pt;
796  }
797 
798  /// Add the connection point (z_value) to the connection
799  /// points that receive the curviline
800  void add_connection_point(const double& z_value,
801  const double& tol = 1.0e-12)
802  {
803  // If we are trying to connect to the initial or final
804  // point then it is not necessary to add this point
805  // to the list since it will explicitly be created when
806  // transforming the curviline to polyline
807  if (std::fabs(z_value - Zeta_start) < tol ||
808  std::fabs(z_value - Zeta_end) < tol)
809  {
810  return;
811  }
812 
813  // We need to deal with repeated connection points,
814  // if the connection point is already in the list then
815  // we will not add it!!!
816  // Search for repeated elements
817  unsigned n_size = Connection_points_pt.size();
818  for (unsigned i = 0; i < n_size; i++)
819  {
820  if (fabs(z_value - Connection_points_pt[i]) < tol)
821  {
822  return;
823  }
824  }
825 
826  // Only add the connection point if it is not the
827  // initial or final zeta value and it is not already on the
828  // list
829  Connection_points_pt.push_back(z_value);
830  }
831 
832  private:
833  /// Pointer to GeomObject that represents this part of the boundary
835 
836  /// Start coordinate in terms of the GeomObject's intrinsic coordinate
837  double Zeta_start;
838 
839  /// End coordinate in terms of the GeomObject's intrinsic coordinate
840  double Zeta_end;
841 
842  /// Number of (initially straight-line) segments that this part of the
843  /// boundary is to be represented by
844  unsigned Nsegment;
845 
846  /// Boundary ID
847  unsigned Boundary_id;
848 
849  /// Boolean to indicate if vertices in polygonal representation
850  /// of the Curviline are spaced (approximately) evenly in arclength
851  /// along the GeomObject [true] or in equal increments in zeta [false]
853 
854  /// Boundary chunk (Used when a boundary is represented by more
855  /// than one polyline
856  unsigned Boundary_chunk;
857 
858  /// Stores the information for connections received on the
859  /// curviline. Used when converting to polyline
861  };
862 
863 
864  //=====================================================================
865  /// Class defining a polyline for use in Triangle Mesh generation
866  //=====================================================================
868  {
869  public:
870  /// Constructor: Takes vectors of vertex coordinates in order
871  /// Also allows the optional specification of a boundary ID -- useful
872  /// in a mesh generation context. If not specified it defaults to zero.
874  const unsigned& boundary_id,
875  const unsigned& boundary_chunk = 0)
880  {
881 #ifdef PARANOID
882  unsigned nvert = Vertex_coordinate.size();
883  for (unsigned i = 0; i < nvert; i++)
884  {
885  if (Vertex_coordinate[i].size() != 2)
886  {
887  std::ostringstream error_stream;
888  error_stream << "TriangleMeshPolyLine should only be used in 2D!\n"
889  << "Your Vector of coordinates, contains data for "
890  << Vertex_coordinate[i].size()
891  << "-dimensional coordinates." << std::endl;
892  throw OomphLibError(error_stream.str(),
893  OOMPH_CURRENT_FUNCTION,
894  OOMPH_EXCEPTION_LOCATION);
895  }
896  }
897 #endif
898  }
899 
900  /// Empty destructor
902 
903  /// Number of vertices
904  unsigned nvertex() const
905  {
906  return Vertex_coordinate.size();
907  }
908 
909  /// Number of segments
910  unsigned nsegment() const
911  {
912  return Vertex_coordinate.size() - 1;
913  }
914 
915  /// Boundary id
916  unsigned boundary_id() const
917  {
918  return Boundary_id;
919  }
920 
921  /// Boundary chunk (Used when a boundary is represented by more
922  /// than one polyline
923  unsigned boundary_chunk() const
924  {
925  return Boundary_chunk;
926  }
927 
928  /// Coordinate vector of i-th vertex (const version)
929  Vector<double> vertex_coordinate(const unsigned& i) const
930  {
931  return Vertex_coordinate[i];
932  }
933 
934  /// Coordinate vector of i-th vertex
936  {
937  return Vertex_coordinate[i];
938  }
939 
940  /// Get first vertex coordinates
942  {
943  vertex = Vertex_coordinate[0];
944  }
945 
946  /// Get last vertex coordinates
948  {
949  vertex = Vertex_coordinate[nvertex() - 1];
950  }
951 
952  /// Output the polyline -- n_sample is ignored
953  void output(std::ostream& outfile, const unsigned& n_sample = 50)
954  {
955  outfile << "ZONE T=\"TriangleMeshPolyLine with boundary ID" << Boundary_id
956  << "\"" << std::endl;
957  unsigned nvert = Vertex_coordinate.size();
958  for (unsigned i = 0; i < nvert; i++)
959  {
960  outfile << Vertex_coordinate[i][0] << " " << Vertex_coordinate[i][1]
961  << std::endl;
962  }
963  }
964 
965  /// Reverse the polyline, this includes the connection information
966  /// and the vertices order
967  void reverse()
968  {
969  // Do the reversing of the connection information
970 
971  // Is there a connection to the initial vertex
972  const bool initial_connection = is_initial_vertex_connected();
973 
974  // Is there a connection to the initial vertex
975  const bool final_connection = is_final_vertex_connected();
976 
977  // If there are any connection at the ends that info. needs to be
978  // reversed
979  if (initial_connection || final_connection)
980  {
981  // Backup the connection information
982 
983  // -------------------------------------------------------------------
984  // Backup the initial vertex connection information
985  // The boundary id
986  const unsigned backup_initial_vertex_connected_bnd_id =
988  // The chunk number
989  const unsigned backup_initial_vertex_connected_chunk =
991  // The vertex number
992  const unsigned backup_initial_vertex_connected_n_vertex =
994  // Is it connected to a curviline
995  const bool backup_initial_vertex_connected_to_curviline =
997  // The s value for the curviline connection
998  const double backup_initial_s_connection = initial_s_connection_value();
999  // The tolerance
1000  const double backup_initial_s_tolerance = tolerance_for_s_connection();
1001 
1002  // -------------------------------------------------------------------
1003  // Backup the final vertex connection information
1004  // The boundary id
1005  const unsigned backup_final_vertex_connected_bnd_id =
1007  // The chunk number
1008  const unsigned backup_final_vertex_connected_chunk =
1010  // The vertex number
1011  const unsigned backup_final_vertex_connected_n_vertex =
1013  // Is it connected to a curviline
1014  const bool backup_final_vertex_connected_to_curviline =
1016  // The s value for the curviline connection
1017  const double backup_final_s_connection = final_s_connection_value();
1018  // The tolerance
1019  const double backup_final_s_tolerance = tolerance_for_s_connection();
1020  // -------------------------------------------------------------------
1021 
1022  // Disconnect the polyline
1023 
1024  // Disconnect the initial vertex
1027 
1028  // Disconnect the final vertex
1031 
1032  // Now reconnected but in inverted order
1033  if (initial_connection)
1034  {
1035  // Set the final vertex as connected
1037  // Copy the boundary id
1039  backup_initial_vertex_connected_bnd_id;
1040  // Copy the chunk number
1042  backup_initial_vertex_connected_chunk;
1043  // Copy the vertex number
1045  backup_initial_vertex_connected_n_vertex;
1046  // Is it connected to a curviline
1047  if (backup_initial_vertex_connected_to_curviline)
1048  {
1049  // Set the final vertex as connected to curviline
1051  // Copy the s value to connected
1052  final_s_connection_value() = backup_initial_s_connection;
1053  // Copy the tolerance
1054  tolerance_for_s_connection() = backup_initial_s_tolerance;
1055  } // if (backup_initial_vertex_connected_to_curviline)
1056 
1057  } // if (initial_connection)
1058 
1059  if (final_connection)
1060  {
1061  // Set the initial vertex as connected
1063  // Copy the boundary id
1065  backup_final_vertex_connected_bnd_id;
1066  // Copy the chunk number
1068  backup_final_vertex_connected_chunk;
1069  // Copy the vertex number
1071  backup_final_vertex_connected_n_vertex;
1072  // Is it connected to a curviline
1073  if (backup_final_vertex_connected_to_curviline)
1074  {
1075  // Set the initial vertex as connected to curviline
1077  // Copy the s value to connected
1078  initial_s_connection_value() = backup_final_s_connection;
1079  // Copy the tolerance
1080  tolerance_for_s_connection() = backup_final_s_tolerance;
1081  } // if (backup_final_vertex_connected_to_curviline)
1082 
1083  } // if (final_connection)
1084 
1085  } // if (initial_connection || final_connection)
1086 
1087  // Do the reversing of the vertices
1088  std::reverse(Vertex_coordinate.begin(), Vertex_coordinate.end());
1089  }
1090 
1091  private:
1092  /// Vector of Vector of vertex coordinates
1094 
1095  /// Boundary ID
1096  unsigned Boundary_id;
1097 
1098  /// Boundary chunk (Used when a boundary is represented by more
1099  /// than one polyline
1100  unsigned Boundary_chunk;
1101  };
1102 
1103 
1104  /// ////////////////////////////////////////////////////////////////////
1105  /// ////////////////////////////////////////////////////////////////////
1106  /// ////////////////////////////////////////////////////////////////////
1107 
1108 
1109  //===================================================================
1110  /// Namespace that allows the specification of a tolerance
1111  /// between vertices at the ends of polylines that are supposed
1112  /// to be at the same position.
1113  //===================================================================
1114  namespace ToleranceForVertexMismatchInPolygons
1115  {
1116  /// Acceptable discrepancy for mismatch in vertex coordinates.
1117  /// In paranoid mode, the code will die if the beginning/end of
1118  /// two adjacent polylines differ by more than that. If the
1119  /// discrepancy is smaller (but nonzero) one of the vertex coordinates
1120  /// get adjusted to match perfectly; without paranoia the vertex
1121  /// coordinates are taken as they come...
1122  extern double Tolerable_error;
1123  } // namespace ToleranceForVertexMismatchInPolygons
1124 
1125  /// ////////////////////////////////////////////////////////////////////
1126  /// ////////////////////////////////////////////////////////////////////
1127  /// ////////////////////////////////////////////////////////////////////
1128 
1129 
1130  //=====================================================================
1131  // Class defining triangle mesh curves. Abstract class for
1132  /// closed curves and open curves. All TriangleMeshCurves are composed
1133  /// of a Vector of TriangleMeshCurveSections
1134  //=====================================================================
1136  {
1137  public:
1138  /// Empty constructor
1143  {
1144  }
1145 
1146  /// Empty destructor
1147  virtual ~TriangleMeshCurve() {}
1148 
1149  /// Number of vertices
1150  virtual unsigned nvertices() const = 0;
1151 
1152  /// Total number of segments
1153  virtual unsigned nsegments() const = 0;
1154 
1155  /// Return max boundary id of associated curves
1156  unsigned max_boundary_id()
1157  {
1158  unsigned max = 0;
1159  unsigned n_curve_section = ncurve_section();
1160  for (unsigned i = 0; i < n_curve_section; i++)
1161  {
1162  unsigned boundary_id = Curve_section_pt[i]->boundary_id();
1163  if (boundary_id > max)
1164  {
1165  max = boundary_id;
1166  }
1167  }
1168  return max;
1169  }
1170 
1171  /// Number of constituent curves
1172  virtual unsigned ncurve_section() const
1173  {
1174  return Curve_section_pt.size();
1175  }
1176 
1177  /// Enable refinement of polylines to create a better
1178  /// representation of curvilinear boundaries (e.g. in free-surface
1179  /// problems). See tutorial for
1180  /// interpretation of the optional argument which specifies the
1181  /// refinement tolerance. It defaults to 0.08 and the smaller the
1182  /// number the finer the surface representation.
1183  void enable_polyline_refinement(const double& tolerance = 0.08)
1184  {
1185  Polyline_refinement_tolerance = tolerance;
1186  // Establish the refinement tolerance for all the
1187  // curve sections on the TriangleMeshCurve
1188  unsigned n_curve_sections = Curve_section_pt.size();
1189  for (unsigned i = 0; i < n_curve_sections; i++)
1190  {
1191  Curve_section_pt[i]->set_refinement_tolerance(
1193  }
1194  }
1195 
1196  /// Set tolerance for refinement of polylines to create a better
1197  /// representation of curvilinear boundaries (e.g. in free-surface
1198  /// problems). See tutorial for
1199  /// interpretation of the refinement tolerance. (The smaller the
1200  /// number the finer the surface representation). If set to
1201  /// a negative value, we're switching off refinement --
1202  /// equivalent to calling disable_polyline_refinement()
1203  void set_polyline_refinement_tolerance(const double& tolerance)
1204  {
1205  Polyline_refinement_tolerance = tolerance;
1206  // Establish the refinement tolerance for all the
1207  // curve sections on the TriangleMeshCurve
1208  unsigned n_curve_sections = Curve_section_pt.size();
1209  for (unsigned i = 0; i < n_curve_sections; i++)
1210  {
1211  Curve_section_pt[i]->set_refinement_tolerance(
1213  }
1214  }
1215 
1216  /// Get tolerance for refinement of polylines to create a better
1217  /// representation of curvilinear boundaries (e.g. in free-surface
1218  /// problems). See tutorial for
1219  /// interpretation. If it's negative refinement is disabled.
1221  {
1223  }
1224 
1225  /// Disable refinement of polylines
1227  {
1229  // Disable the refinement tolerance for all the
1230  // curve sections on the TriangleMeshCurve
1231  unsigned n_curve_sections = Curve_section_pt.size();
1232  for (unsigned i = 0; i < n_curve_sections; i++)
1233  {
1234  Curve_section_pt[i]->disable_refinement_tolerance();
1235  }
1236  }
1237 
1238  /// Enable unrefinement of polylines to avoid unnecessarily large
1239  /// numbers of elements on of curvilinear boundaries (e.g. in free-surface
1240  /// problems). See tutorial for
1241  /// interpretation of the optional argument which specifies the
1242  /// unrefinement tolerance. It defaults to 0.04 and the larger the number
1243  /// the more agressive we are when removing unnecessary vertices on
1244  /// gently curved polylines.
1245  void enable_polyline_unrefinement(const double& tolerance = 0.04)
1246  {
1247  Polyline_unrefinement_tolerance = tolerance;
1248  // Establish the unrefinement tolerance for all the
1249  // curve sections on the TriangleMeshCurve
1250  unsigned n_curve_sections = Curve_section_pt.size();
1251  for (unsigned i = 0; i < n_curve_sections; i++)
1252  {
1253  Curve_section_pt[i]->set_unrefinement_tolerance(
1255  }
1256  }
1257 
1258  /// Set tolerance for unrefinement of polylines
1259  /// to avoid unnecessarily large
1260  /// numbers of elements on of curvilinear boundaries (e.g. in free-surface
1261  /// problems). See tutorial for
1262  /// interpretation of the optional argument which specifies the
1263  /// unrefinement tolerance. It defaults to 0.04 and the larger the number
1264  /// the more agressive we are when removing unnecessary vertices on
1265  /// gently curved polylines. If set to
1266  /// a negative value, we're switching off unrefinement --
1267  /// equivalent to calling disable_polyline_unrefinement()
1268  void set_polyline_unrefinement_tolerance(const double& tolerance)
1269  {
1270  Polyline_unrefinement_tolerance = tolerance;
1271  // Establish the unrefinement tolerance for all the
1272  // curve sections on the TriangleMeshCurve
1273  unsigned n_curve_sections = Curve_section_pt.size();
1274  for (unsigned i = 0; i < n_curve_sections; i++)
1275  {
1276  Curve_section_pt[i]->set_unrefinement_tolerance(
1278  }
1279  }
1280 
1281  /// Get tolerance for unrefinement of polylines to create a better
1282  /// representation of curvilinear boundaries (e.g. in free-surface
1283  /// problems). See tutorial for
1284  /// interpretation. If it's negative unrefinement is disabled.
1286  {
1288  }
1289 
1290  /// Disable unrefinement of polylines
1292  {
1294  // Disable the unrefinement tolerance for all the
1295  // curve sections on the TriangleMeshCurve
1296  unsigned n_curve_sections = Curve_section_pt.size();
1297  for (unsigned i = 0; i < n_curve_sections; i++)
1298  {
1299  Curve_section_pt[i]->disable_unrefinement_tolerance();
1300  }
1301  }
1302 
1303  /// Output each sub-boundary at n_sample (default: 50) points
1304  virtual void output(std::ostream& outfile,
1305  const unsigned& n_sample = 50) = 0;
1306 
1307  /// Pointer to i-th constituent curve section
1308  virtual TriangleMeshCurveSection* curve_section_pt(const unsigned& i) const
1309  {
1310  return Curve_section_pt[i];
1311  }
1312 
1313  /// Pointer to i-th constituent curve section
1314  virtual TriangleMeshCurveSection*& curve_section_pt(const unsigned& i)
1315  {
1316  return Curve_section_pt[i];
1317  }
1318 
1319  protected:
1320  /// Vector of curve sections
1322 
1323  private:
1324  /// Tolerance for refinement of polylines (neg if refinement is disabled)
1326 
1327  /// Tolerance for unrefinement of polylines (neg if refinement is disabled)
1329  };
1330 
1331  /// ////////////////////////////////////////////////////////////////////
1332  /// ////////////////////////////////////////////////////////////////////
1333  /// ////////////////////////////////////////////////////////////////////
1334 
1335  //=====================================================================
1336  /// Base class defining a closed curve for the Triangle mesh generation
1337  //=====================================================================
1339  {
1340  public:
1341  /// Constructor prototype
1344  const Vector<double>& internal_point_pt = Vector<double>(0),
1345  const bool& is_internal_point_fixed = false);
1346 
1347  /// Empty destructor
1349 
1350  /// Number of vertices
1351  unsigned nvertices() const
1352  {
1353  unsigned n_curve_section = ncurve_section();
1354  unsigned n_vertices = 0;
1355  for (unsigned j = 0; j < n_curve_section; j++)
1356  {
1357  // Storing the number of the vertices
1358  n_vertices += Curve_section_pt[j]->nvertex() - 1;
1359  }
1360  // If there's just one boundary. All the vertices should be counted
1361  if (n_curve_section == 1)
1362  {
1363  n_vertices += 1;
1364  }
1365  return n_vertices;
1366  }
1367 
1368  /// Total number of segments
1369  unsigned nsegments() const
1370  {
1371  unsigned n_curve_section = ncurve_section();
1372  unsigned nseg = 0;
1373  for (unsigned j = 0; j < n_curve_section; j++)
1374  {
1375  nseg += Curve_section_pt[j]->nsegment();
1376  }
1377  // If there's just one boundary poly line we have another segment
1378  // connecting the last vertex to the first one
1379  if (n_curve_section == 1)
1380  {
1381  nseg += 1;
1382  }
1383  return nseg;
1384  }
1385 
1386  /// Output each sub-boundary at n_sample (default: 50) points
1387  void output(std::ostream& outfile, const unsigned& n_sample = 50)
1388  {
1389  unsigned nb = Curve_section_pt.size();
1390  for (unsigned i = 0; i < nb; i++)
1391  {
1392  Curve_section_pt[i]->output(outfile, n_sample);
1393  }
1394 
1395  if (!Internal_point_pt.empty())
1396  {
1397  outfile << "ZONE T=\"Internal point\"\n";
1398  outfile << Internal_point_pt[0] << " " << Internal_point_pt[1] << "\n";
1399  }
1400  }
1401 
1402  /// Coordinates of the internal point
1404  {
1405  return Internal_point_pt;
1406  }
1407 
1408  /// Coordinates of the internal point
1410  {
1411  return Internal_point_pt;
1412  }
1413 
1414  /// Fix the internal point (i.e. do not allow our automatic machinery
1415  /// to update it)
1417  {
1418  Is_internal_point_fixed = true;
1419  }
1420 
1421  /// Unfix the internal point (i.e. allow our automatic machinery
1422  /// to update it)
1424  {
1425  Is_internal_point_fixed = false;
1426  }
1427 
1428  /// Test whether the internal point is fixed
1430  {
1431  return Is_internal_point_fixed;
1432  }
1433 
1434  protected:
1435  /// Vector of vertex coordinates
1437 
1438  /// Indicate whether the internal point should be updated automatically
1440  };
1441 
1442  /// ////////////////////////////////////////////////////////////////////
1443  /// ////////////////////////////////////////////////////////////////////
1444  /// ////////////////////////////////////////////////////////////////////
1445 
1446 
1447  //=====================================================================
1448  /// Class defining a closed polygon for the Triangle mesh generation
1449  //=====================================================================
1451  {
1452  public:
1453  /// Constructor: Specify vector of pointers to
1454  /// TriangleMeshCurveSection that define the boundary of the segments of the
1455  /// polygon. Each TriangleMeshCurveSection has its own boundary ID and can
1456  /// contain multiple (straight-line) segments. For consistency across the
1457  /// various uses of this class, we insist that the closed boundary
1458  /// is represented by at least two separate TriangleMeshCurveSection
1459  /// whose joint vertices must be specified in both.
1460  /// (This is to allow the setup of unique boundary coordinate(s)
1461  /// around the polygon.) This may seem slightly annoying
1462  /// in cases where a polygon really only represents a single
1463  /// boundary, but...
1464  /// Note: The specified vector of pointers must consist of only
1465  /// TriangleMeshPolyLine elements. There is a checking on the PARANOID
1466  /// mode for this constraint
1468  const Vector<TriangleMeshCurveSection*>& boundary_polyline_pt,
1469  const Vector<double>& internal_point_pt = Vector<double>(0),
1470  const bool& is_internal_point_fixed = false);
1471 
1472  /// Empty virtual destructor
1474 
1475  /// Number of constituent curves
1476  unsigned ncurve_section() const
1477  {
1478  return npolyline();
1479  }
1480 
1481  /// Number of constituent polylines
1482  unsigned npolyline() const
1483  {
1484  return Curve_section_pt.size();
1485  }
1486 
1487  /// Pointer to i-th constituent polyline
1488  TriangleMeshPolyLine* polyline_pt(const unsigned& i) const
1489  {
1490  TriangleMeshPolyLine* tmp_polyline =
1491  dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i]);
1492 #ifdef PARANOID
1493  if (tmp_polyline == NULL)
1494  {
1495  std::ostringstream error_stream;
1496  error_stream
1497  << "The (" << i << ") TriangleMeshCurveSection is not a "
1498  << "TriangleMeshPolyLine\nThe TriangleMeshPolygon object"
1499  << "is constituent of only TriangleMeshPolyLine objects.\n"
1500  << "The problem could be generated when changing the constituent "
1501  << "objects of the TriangleMeshPolygon.\nCheck where you got "
1502  << "access to this objects and review that you are not introducing "
1503  << "any other objects than TriangleMeshPolyLines" << std::endl;
1504  throw OomphLibError(
1505  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1506  }
1507 #endif
1508  return tmp_polyline;
1509  }
1510 
1511  /// Pointer to i-th constituent polyline
1513  {
1514  TriangleMeshPolyLine* tmp_polyline =
1515  dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i]);
1516 #ifdef PARANOID
1517  if (tmp_polyline == NULL)
1518  {
1519  std::ostringstream error_stream;
1520  error_stream
1521  << "The (" << i << ") TriangleMeshCurveSection is not a "
1522  << "TriangleMeshPolyLine\nThe TriangleMeshPolygon object"
1523  << "is constituent of only TriangleMeshPolyLine objects.\n"
1524  << "The problem could be generated when changing the constituent "
1525  << "objects of the TriangleMeshPolygon.\nCheck where you got "
1526  << "access to this objects and review that you are not introducing "
1527  << "any other objects than TriangleMeshPolyLines" << std::endl;
1528  throw OomphLibError(
1529  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1530  }
1531 #endif
1532  return tmp_polyline;
1533  }
1534 
1535  /// Return vector of boundary ids of associated polylines
1537  {
1538  // Get the number of polylines
1539  unsigned nline = npolyline();
1540  Vector<unsigned> boundary_id(nline);
1541 
1542  // Loop over the polyline to get the id
1543  for (unsigned iline = 0; iline < nline; iline++)
1544  {
1545  boundary_id[iline] = Curve_section_pt[iline]->boundary_id();
1546  }
1547  return boundary_id;
1548  }
1549 
1550  /// Is re-distribution of polyline segments in the curve
1551  /// between different boundaries during adaptation enabled?
1553  {
1555  }
1556 
1557  /// Enable re-distribution of polyline segments in the curve
1558  /// between different boundaries during adaptation
1560  {
1562  }
1563 
1564  /// Disable re-distribution of polyline segments in the curve
1565  /// between different boundaries during adaptation
1567  {
1569  }
1570 
1571  /// Test whether curve can update reference
1573  {
1574  return Can_update_configuration;
1575  }
1576 
1577  /// Virtual function that should be overloaded to update the polygons
1578  /// reference configuration
1580  {
1581  std::ostringstream error_stream;
1582  error_stream
1583  << "Broken Default Called\n"
1584  << "This function should be overloaded if Can_update_configuration = "
1585  "true,"
1586  << "\nwhich indicates that the curve in it polylines parts can update "
1587  "its "
1588  << "own position (i.e. it\n"
1589  << "may be a rigid body. Otherwise the update will be via a FaceMesh \n"
1590  << "representation of the boundary which is appropriate for \n"
1591  << "general deforming surfaces\n";
1592 
1593  throw OomphLibError(
1594  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1595  }
1596 
1597  /// Test whether the polygon is fixed or not
1598  bool is_fixed() const
1599  {
1600  return Polygon_fixed;
1601  }
1602 
1603  /// Set the polygon to be fixed
1604  void set_fixed()
1605  {
1606  Polygon_fixed = true;
1607  }
1608 
1609  /// Set the polygon to be allowed to move (default)
1611  {
1612  Polygon_fixed = false;
1613  }
1614 
1615  protected:
1616  /// Is re-distribution of polyline segments between different
1617  /// boundaries during adaptation enabled? (Default: false)
1619 
1620  /// Boolean flag to indicate whether the polygon can update its
1621  /// own reference configuration after it has moved i.e. if it is
1622  /// upgraded to a rigid body rather than being a free surface (default
1623  /// false)
1625 
1626  private:
1627  /// Boolean flag to indicate whether the polygon can move
1628  /// (default false because if it doesn't move this will just lead to
1629  /// wasted work)
1631  };
1632 
1633  /// //////////////////////////////////////////////////////////////////////
1634  /// //////////////////////////////////////////////////////////////////////
1635  /// //////////////////////////////////////////////////////////////////////
1636 
1637  //=====================================================================
1638  /// Base class defining an open curve for the Triangle mesh generation
1639  /// Basically used to define internal boundaries on the mesh
1640  //=====================================================================
1642  {
1643  public:
1644  /// Constructor
1647 
1648  /// Empty destructor
1650 
1651  /// Number of vertices
1652  unsigned nvertices() const
1653  {
1654  unsigned n_vertices = 0;
1655  unsigned n_curve_section = ncurve_section();
1656  for (unsigned i = 0; i < n_curve_section; i++)
1657  n_vertices += Curve_section_pt[i]->nvertex();
1658  // If there's just one boundary. All the vertices should be counted
1659  if (n_curve_section == 1)
1660  {
1661  n_vertices += 1;
1662  }
1663  return n_vertices;
1664  }
1665 
1666  /// Total number of segments
1667  unsigned nsegments() const
1668  {
1669  unsigned n_curve_section = ncurve_section();
1670  unsigned nseg = 0;
1671  for (unsigned j = 0; j < n_curve_section; j++)
1672  {
1673  nseg += Curve_section_pt[j]->nsegment();
1674  }
1675  return nseg;
1676  }
1677 
1678  /// Output each sub-boundary at n_sample (default: 50) points
1679  void output(std::ostream& outfile, const unsigned& n_sample = 50)
1680  {
1681  unsigned nb = Curve_section_pt.size();
1682  for (unsigned i = 0; i < nb; i++)
1683  {
1684  Curve_section_pt[i]->output(outfile, n_sample);
1685  }
1686  }
1687 
1688  /// Pointer to i-th constituent polyline
1689  TriangleMeshPolyLine* polyline_pt(const unsigned& i) const
1690  {
1691  TriangleMeshPolyLine* tmp_polyline =
1692  dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i]);
1693 #ifdef PARANOID
1694  if (tmp_polyline == NULL)
1695  {
1696  std::ostringstream error_stream;
1697  error_stream << "The (" << i
1698  << ")-th TriangleMeshCurveSection is not a "
1699  << "TriangleMeshPolyLine.\nPlease make sure that when you"
1700  << "first created this object the (" << i << ")-th\n"
1701  << "TriangleCurveSection is a TriangleMeshPolyLine."
1702  << std::endl;
1703  throw OomphLibError(
1704  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1705  }
1706 #endif
1707  return tmp_polyline;
1708  }
1709 
1710  /// Pointer to i-th constituent polyline
1712  {
1713  TriangleMeshPolyLine* tmp_polyline =
1714  dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i]);
1715 #ifdef PARANOID
1716  if (tmp_polyline == NULL)
1717  {
1718  std::ostringstream error_stream;
1719  error_stream << "The (" << i
1720  << ")-th TriangleMeshCurveSection is not a "
1721  << "TriangleMeshPolyLine.\nPlease make sure that when you"
1722  << "first created this object the (" << i << ")-th\n"
1723  << "TriangleCurveSection is a TriangleMeshPolyLine."
1724  << std::endl;
1725  throw OomphLibError(
1726  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1727  }
1728 #endif
1729  return tmp_polyline;
1730  }
1731  };
1732 
1733  //==============start_of_geometry_helper_functions_class================
1734  /// Contains functions which define the geometry of the mesh, i.e.
1735  /// regions, boundaries, etc.
1736  //======================================================================
1738  {
1739  public:
1740  /// Public static flag to suppress warning; defaults to false
1742 
1743  /// Empty constructor
1745 
1746  /// Broken copy constructor
1748  const UnstructuredTwoDMeshGeometryBase& dummy) = delete;
1749 
1750  /// Broken assignment operator
1752 
1753  /// Empty destructor
1755 
1756  /// Return the number of regions specified by attributes
1757  unsigned nregion()
1758  {
1759  return Region_element_pt.size();
1760  }
1761 
1762  /// Return the number of elements in the i-th region
1763  unsigned nregion_element(const unsigned& i)
1764  {
1765  // Create an iterator to iterate over Region_element_pt
1766  std::map<unsigned, Vector<FiniteElement*>>::iterator it;
1767 
1768  // Find the entry of Region_element_pt associated with the i-th region
1769  it = Region_element_pt.find(i);
1770 
1771  // If there is an entry associated with the i-th region
1772  if (it != Region_element_pt.end())
1773  {
1774  return (it->second).size();
1775  }
1776  else
1777  {
1778  return 0;
1779  }
1780  } // End of nregion_element
1781 
1782  /// Return the e-th element in the i-th region
1783  FiniteElement* region_element_pt(const unsigned& i, const unsigned& e)
1784  {
1785  // Create an iterator to iterate over Region_element_pt
1786  std::map<unsigned, Vector<FiniteElement*>>::iterator it;
1787 
1788  // Find the entry of Region_element_pt associated with the i-th region
1789  it = Region_element_pt.find(i);
1790 
1791  // If there is an entry associated with the i-th region
1792  if (it != Region_element_pt.end())
1793  {
1794  // Return a pointer to the e-th element in the i-th region
1795  return (it->second)[e];
1796  }
1797  else
1798  {
1799  // Create a stringstream object
1800  std::stringstream error_message;
1801 
1802  // Create the error message
1803  error_message << "There are no regions associated with "
1804  << "region ID " << i << ".";
1805 
1806  // Throw an error
1807  throw OomphLibError(error_message.str(),
1808  OOMPH_CURRENT_FUNCTION,
1809  OOMPH_EXCEPTION_LOCATION);
1810  }
1811  } // End of region_element_pt
1812 
1813  /// Return the number of attributes used in the mesh
1815  {
1816  return Region_attribute.size();
1817  }
1818 
1819  /// Return the attribute associated with region i
1820  double region_attribute(const unsigned& i)
1821  {
1822  return Region_attribute[i];
1823  }
1824 
1825  /// Return the geometric object associated with the b-th boundary or
1826  /// null if the boundary has associated geometric object.
1828  {
1829  std::map<unsigned, GeomObject*>::iterator it;
1830  it = Boundary_geom_object_pt.find(b);
1831  if (it == Boundary_geom_object_pt.end())
1832  {
1833  return 0;
1834  }
1835  else
1836  {
1837  return (*it).second;
1838  }
1839  }
1840 
1841  /// Return direct access to the geometric object storage
1842  std::map<unsigned, GeomObject*>& boundary_geom_object_pt()
1843  {
1844  return Boundary_geom_object_pt;
1845  }
1846 
1847  /// Return access to the vector of boundary coordinates associated
1848  /// with each geometric object
1849  std::map<unsigned, Vector<double>>& boundary_coordinate_limits()
1850  {
1852  }
1853 
1854  /// Return access to the coordinate limits associated with
1855  /// the geometric object associated with boundary b
1857  {
1858  std::map<unsigned, Vector<double>>::iterator it;
1859  it = Boundary_coordinate_limits.find(b);
1860  if (it == Boundary_coordinate_limits.end())
1861  {
1862  throw OomphLibError(
1863  "No coordinate limits associated with this boundary\n",
1864  OOMPH_CURRENT_FUNCTION,
1865  OOMPH_EXCEPTION_LOCATION);
1866  }
1867  else
1868  {
1869  return (*it).second;
1870  }
1871  }
1872 
1873  /// Return the number of elements adjacent to boundary b in region r
1874  inline unsigned nboundary_element_in_region(const unsigned& b,
1875  const unsigned& r) const
1876  {
1877  // Need to use a constant iterator here to keep the function "const".
1878  // Return an iterator to the appropriate entry, if we find it
1879  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
1880  Boundary_region_element_pt[b].find(r);
1881  if (it != Boundary_region_element_pt[b].end())
1882  {
1883  return (it->second).size();
1884  }
1885  // Otherwise there are no elements adjacent to boundary b in the region r
1886  else
1887  {
1888  return 0;
1889  }
1890  }
1891 
1892  /// Return pointer to the e-th element adjacent to boundary b in region r
1894  const unsigned& r,
1895  const unsigned& e) const
1896  {
1897  // Use a constant iterator here to keep function "const" overall
1898  std::map<unsigned, Vector<FiniteElement*>>::const_iterator it =
1899  Boundary_region_element_pt[b].find(r);
1900  if (it != Boundary_region_element_pt[b].end())
1901  {
1902  return (it->second)[e];
1903  }
1904  else
1905  {
1906  return 0;
1907  }
1908  }
1909 
1910  /// Return face index of the e-th element adjacent to boundary b in region r
1911  int face_index_at_boundary_in_region(const unsigned& b,
1912  const unsigned& r,
1913  const unsigned& e) const
1914  {
1915  // Use a constant iterator here to keep function "const" overall
1916  std::map<unsigned, Vector<int>>::const_iterator it =
1917  Face_index_region_at_boundary[b].find(r);
1918  if (it != Face_index_region_at_boundary[b].end())
1919  {
1920  return (it->second)[e];
1921  }
1922  else
1923  {
1924  std::ostringstream error_message;
1925  error_message << "Face indices not set up for boundary " << b
1926  << " in region " << r << "\n";
1927  error_message << "This probably means that the boundary is not "
1928  "adjacent to region\n";
1929  throw OomphLibError(error_message.str(),
1930  OOMPH_CURRENT_FUNCTION,
1931  OOMPH_EXCEPTION_LOCATION);
1932  }
1933  }
1934 
1935  /// Return pointer to the current polyline that describes
1936  /// the b-th mesh boundary
1938  {
1939  std::map<unsigned, TriangleMeshCurveSection*>::iterator it =
1940  Boundary_curve_section_pt.find(b);
1941  // Search for the polyline associated with the given boundary
1942  if (it != Boundary_curve_section_pt.end())
1943  {
1944  return dynamic_cast<TriangleMeshPolyLine*>(
1946  }
1947  // If the boundary was not established then return 0, or a null pointer
1948  return 0;
1949  }
1950 
1951  /// Gets a pointer to a set with all the nodes related with a
1952  /// boundary
1953  std::map<unsigned, std::set<Node*>>& nodes_on_boundary_pt()
1954  {
1955  return Nodes_on_boundary_pt;
1956  }
1957 
1958  /// Gets the vertex number on the destination polyline (used
1959  /// to create the connections among shared boundaries)
1961  TriangleMeshPolyLine* dst_polyline_pt,
1962  Vector<double>& vertex_coordinates,
1963  unsigned& vertex_number);
1964 
1965  /// Sort the polylines coming from joining them. Check whether
1966  /// it is necessary to reverse them or not. Used when joining two curve
1967  /// polylines in order to create a polygon
1969  Vector<TriangleMeshPolyLine*>& polylines_pt, unsigned& index);
1970 
1971  /// Sort the polylines coming from joining them. Check whether
1972  /// it is necessary to reverse them or not. Used when joining polylines
1973  /// and they still do not create a polygon
1975  Vector<TriangleMeshPolyLine*>& polylines_pt,
1976  unsigned& index_halo_start,
1977  unsigned& index_halo_end);
1978 
1979  /// Helper function that checks if a given point is inside a polygon
1980  /// (a set of sorted vertices that connected create a polygon)
1981  bool is_point_inside_polygon_helper(Vector<Vector<double>> polygon_vertices,
1982  Vector<double> point);
1983 
1984  /// Enables the creation of points (by Triangle) on the outer and
1985  /// internal boundaries
1987  {
1989  }
1990 
1991  /// Disables the creation of points (by Triangle) on the outer and
1992  /// internal boundaries
1994  {
1996  }
1997 
1998  /// Returns the status of the variable
1999  /// Allow_automatic_creation_of_vertices_on_boundaries
2001  {
2003  }
2004 
2005 #ifdef OOMPH_HAS_MPI
2006 
2007  /// Flush the boundary segment node storage
2008  void flush_boundary_segment_node(const unsigned& b)
2009  {
2010  Boundary_segment_node_pt[b].clear();
2011  }
2012 
2013  /// Set the number of segments associated with a boundary
2014  void set_nboundary_segment_node(const unsigned& b, const unsigned& s)
2015  {
2016  Boundary_segment_node_pt[b].resize(s);
2017  }
2018 
2019  /// Return the number of segments associated with a boundary
2020  unsigned nboundary_segment(const unsigned& b)
2021  {
2022  return Boundary_segment_node_pt[b].size();
2023  }
2024 
2025  /// Return the number of segments associated with a boundary
2026  unsigned long nboundary_segment_node(const unsigned& b)
2027  {
2028  unsigned ntotal_nodes = 0;
2029  unsigned nsegments = Boundary_segment_node_pt[b].size();
2030  for (unsigned is = 0; is < nsegments; is++)
2031  {
2032  ntotal_nodes += nboundary_segment_node(b, is);
2033  }
2034  return ntotal_nodes;
2035  }
2036 
2037  /// Return the number of nodes associated with a given segment of a
2038  /// boundary
2039  unsigned long nboundary_segment_node(const unsigned& b, const unsigned& s)
2040  {
2041  return Boundary_segment_node_pt[b][s].size();
2042  }
2043 
2044  /// Add the node node_pt to the b-th boundary and the s-th segment of
2045  /// the mesh
2046  void add_boundary_segment_node(const unsigned& b,
2047  const unsigned& s,
2048  Node* const& node_pt)
2049  {
2050  // Get the size of the Boundary_node_pt vector
2051  unsigned nbound_seg_node = nboundary_segment_node(b, s);
2052  bool node_already_on_this_boundary_segment = false;
2053 
2054  // Loop over the vector
2055  for (unsigned n = 0; n < nbound_seg_node; n++)
2056  {
2057  // Is the current node here already?
2058  if (node_pt == Boundary_segment_node_pt[b][s][n])
2059  {
2060  node_already_on_this_boundary_segment = true;
2061  }
2062  }
2063 
2064  // Add the base node pointer to the vector if it's not there already
2065  if (!node_already_on_this_boundary_segment)
2066  {
2067  Boundary_segment_node_pt[b][s].push_back(node_pt);
2068  }
2069  }
2070 
2071  /// Flag used at the setup_boundary_coordinate function to know
2072  /// if initial zeta values for segments have been assigned
2073  std::map<unsigned, bool> Assigned_segments_initial_zeta_values;
2074 
2075  /// Return direct access to the initial coordinates of a boundary
2076  std::map<unsigned, Vector<double>>& boundary_initial_coordinate()
2077  {
2079  }
2080 
2081  /// Return direct access to the final coordinates of a boundary
2082  std::map<unsigned, Vector<double>>& boundary_final_coordinate()
2083  {
2085  }
2086 
2087  /// Return direct access to the initial zeta coordinate of a
2088  /// boundary
2089  std::map<unsigned, Vector<double>>& boundary_initial_zeta_coordinate()
2090  {
2092  }
2093 
2094  /// Return direct access to the final zeta coordinates of a
2095  /// boundary
2096  std::map<unsigned, Vector<double>>& boundary_final_zeta_coordinate()
2097  {
2099  }
2100 
2101  /// Return the info. to know if it is necessary to reverse the
2102  /// segment based on a previous mesh
2103  std::map<unsigned, Vector<unsigned>>& boundary_segment_inverted()
2104  {
2106  }
2107 
2108  /// Return direct access to the initial coordinates for the
2109  /// segments that are part of a boundary
2110  std::map<unsigned, Vector<Vector<double>>>& boundary_segment_initial_coordinate()
2111  {
2113  }
2114 
2115  /// Return direct access to the final coordinates for the
2116  /// segments that are part of a boundary
2117  std::map<unsigned, Vector<Vector<double>>>& boundary_segment_final_coordinate()
2118  {
2120  }
2121 
2122  /// Return direct access to the initial arclength for the
2123  /// segments that are part of a boundary
2124  std::map<unsigned, Vector<double>>& boundary_segment_initial_arclength()
2125  {
2127  }
2128 
2129  /// Return direct access to the final arclength for the
2130  /// segments that are part of a boundary
2131  std::map<unsigned, Vector<double>>& boundary_segment_final_arclength()
2132  {
2134  }
2135 
2136  /// Return direct access to the initial zeta for the
2137  /// segments that are part of a boundary
2138  std::map<unsigned, Vector<double>>& boundary_segment_initial_zeta()
2139  {
2141  }
2142 
2143  /// Return direct access to the final zeta for the
2144  /// segments that are part of a boundary
2145  std::map<unsigned, Vector<double>>& boundary_segment_final_zeta()
2146  {
2148  }
2149 
2150  /// Return the initial zeta for the segments that are
2151  /// part of a boundary
2153  {
2154  std::map<unsigned, Vector<double>>::iterator it =
2156 
2157 #ifdef PARANOID
2158 
2159  if (it == Boundary_segment_initial_zeta.end())
2160  {
2161  std::stringstream error_message;
2162  error_message << "The boundary (" << b
2163  << ") has no segments associated with it!!\n\n";
2164  throw OomphLibError(error_message.str(),
2165  OOMPH_CURRENT_FUNCTION,
2166  OOMPH_EXCEPTION_LOCATION);
2167  }
2168 
2169 #endif // PARANOID
2170 
2171  return (*it).second;
2172  }
2173 
2174  /// Return the final zeta for the segments that are
2175  /// part of a boundary
2177  {
2178  std::map<unsigned, Vector<double>>::iterator it =
2180 
2181 #ifdef PARANOID
2182 
2183  if (it == Boundary_segment_final_zeta.end())
2184  {
2185  std::stringstream error_message;
2186  error_message << "The boundary (" << b
2187  << ") has no segments associated with it!!\n\n";
2188  throw OomphLibError(error_message.str(),
2189  OOMPH_CURRENT_FUNCTION,
2190  OOMPH_EXCEPTION_LOCATION);
2191  }
2192 
2193 #endif // PARANOID
2194 
2195  return (*it).second;
2196  }
2197 
2198  /// Return the initial coordinates for the boundary
2200  {
2201  std::map<unsigned, Vector<double>>::iterator it =
2203 
2204 #ifdef PARANOID
2205 
2206  if (it == Boundary_initial_coordinate.end())
2207  {
2208  std::stringstream error_message;
2209  error_message << "The boundary (" << b
2210  << ") has not established initial coordinates\n";
2211  throw OomphLibError(error_message.str(),
2212  OOMPH_CURRENT_FUNCTION,
2213  OOMPH_EXCEPTION_LOCATION);
2214  }
2215 
2216 #endif
2217  return (*it).second;
2218  }
2219 
2220  /// Return the final coordinates for the boundary
2222  {
2223  std::map<unsigned, Vector<double>>::iterator it =
2224  Boundary_final_coordinate.find(b);
2225 
2226 #ifdef PARANOID
2227 
2228  if (it == Boundary_final_coordinate.end())
2229  {
2230  std::stringstream error_message;
2231  error_message << "The boundary (" << b
2232  << ") has not established final coordinates\n";
2233  throw OomphLibError(error_message.str(),
2234  OOMPH_CURRENT_FUNCTION,
2235  OOMPH_EXCEPTION_LOCATION);
2236  }
2237 
2238 #endif
2239 
2240  return (*it).second;
2241  }
2242 
2243  /// Return the info. to know if it is necessary to reverse the
2244  /// segment based on a previous mesh
2245  const Vector<unsigned> boundary_segment_inverted(const unsigned& b) const
2246  {
2247  std::map<unsigned, Vector<unsigned>>::const_iterator it =
2248  Boundary_segment_inverted.find(b);
2249 
2250 #ifdef PARANOID
2251 
2252  if (it == Boundary_segment_inverted.end())
2253  {
2254  std::stringstream error_message;
2255  error_message << "The boundary (" << b
2256  << ") has not established inv. segments info\n";
2257  throw OomphLibError(error_message.str(),
2258  OOMPH_CURRENT_FUNCTION,
2259  OOMPH_EXCEPTION_LOCATION);
2260  }
2261 
2262 #endif
2263 
2264  return (*it).second;
2265  }
2266 
2267  /// Return the info. to know if it is necessary to reverse the
2268  /// segment based on a previous mesh
2270  {
2271  std::map<unsigned, Vector<unsigned>>::iterator it =
2272  Boundary_segment_inverted.find(b);
2273 
2274 #ifdef PARANOID
2275 
2276  if (it == Boundary_segment_inverted.end())
2277  {
2278  std::stringstream error_message;
2279  error_message << "The boundary (" << b
2280  << ") has not established inv. segments info\n";
2281  throw OomphLibError(error_message.str(),
2282  OOMPH_CURRENT_FUNCTION,
2283  OOMPH_EXCEPTION_LOCATION);
2284  }
2285 
2286 #endif
2287 
2288  return (*it).second;
2289  }
2290 
2291  /// Return the initial zeta coordinate for the boundary
2293  {
2294  std::map<unsigned, Vector<double>>::iterator it =
2296 
2297 #ifdef PARANOID
2298 
2299  if (it == Boundary_initial_zeta_coordinate.end())
2300  {
2301  std::stringstream error_message;
2302  error_message << "The boundary (" << b
2303  << ") has not established initial zeta "
2304  << "coordinate\n";
2305  throw OomphLibError(error_message.str(),
2306  OOMPH_CURRENT_FUNCTION,
2307  OOMPH_EXCEPTION_LOCATION);
2308  }
2309 
2310 #endif
2311 
2312  return (*it).second;
2313  }
2314 
2315  /// Return the final zeta coordinate for the boundary
2317  {
2318  std::map<unsigned, Vector<double>>::iterator it =
2320 
2321 #ifdef PARANOID
2322 
2323  if (it == Boundary_final_zeta_coordinate.end())
2324  {
2325  std::stringstream error_message;
2326  error_message << "The boundary (" << b
2327  << ") has not established final zeta coordinate\n";
2328  throw OomphLibError(error_message.str(),
2329  OOMPH_CURRENT_FUNCTION,
2330  OOMPH_EXCEPTION_LOCATION);
2331  }
2332 
2333 #endif
2334 
2335  return (*it).second;
2336  }
2337 
2338  /// Return the initial arclength for the segments that are
2339  /// part of a boundary
2341  {
2342  std::map<unsigned, Vector<double>>::iterator it =
2344 
2345 #ifdef PARANOID
2346 
2347  if (it == Boundary_segment_initial_arclength.end())
2348  {
2349  std::stringstream error_message;
2350  error_message << "The boundary (" << b
2351  << ") has no segments associated with it!!\n\n";
2352  throw OomphLibError(error_message.str(),
2353  OOMPH_CURRENT_FUNCTION,
2354  OOMPH_EXCEPTION_LOCATION);
2355  }
2356 
2357 #endif
2358 
2359  return (*it).second;
2360  }
2361 
2362  /// Return the final arclength for the segments that are
2363  /// part of a boundary
2365  {
2366  std::map<unsigned, Vector<double>>::iterator it =
2368 
2369 #ifdef PARANOID
2370 
2371  if (it == Boundary_segment_final_arclength.end())
2372  {
2373  std::stringstream error_message;
2374  error_message << "The boundary (" << b
2375  << ") has no segments associated with it!!\n\n";
2376  throw OomphLibError(error_message.str(),
2377  OOMPH_CURRENT_FUNCTION,
2378  OOMPH_EXCEPTION_LOCATION);
2379  }
2380 
2381 #endif
2382 
2383  return (*it).second;
2384  }
2385 
2386  /// Return the initial coordinates for the segments that are
2387  /// part of a boundary
2389  const unsigned& b)
2390  {
2391  std::map<unsigned, Vector<Vector<double>>>::iterator it =
2393 
2394 #ifdef PARANOID
2395 
2396  if (it == Boundary_segment_initial_coordinate.end())
2397  {
2398  std::stringstream error_message;
2399  error_message << "The boundary (" << b
2400  << ") has no segments associated with it!!\n\n";
2401  throw OomphLibError(error_message.str(),
2402  OOMPH_CURRENT_FUNCTION,
2403  OOMPH_EXCEPTION_LOCATION);
2404  }
2405 
2406 #endif
2407 
2408  return (*it).second;
2409  }
2410 
2411  /// Return the final coordinates for the segments that are
2412  /// part of a boundary
2414  {
2415  std::map<unsigned, Vector<Vector<double>>>::iterator it =
2417 
2418 #ifdef PARANOID
2419 
2420  if (it == Boundary_segment_final_coordinate.end())
2421  {
2422  std::stringstream error_message;
2423  error_message << "The boundary (" << b
2424  << ") has no segments associated with it!!\n\n";
2425  throw OomphLibError(error_message.str(),
2426  OOMPH_CURRENT_FUNCTION,
2427  OOMPH_EXCEPTION_LOCATION);
2428  }
2429 
2430 #endif
2431 
2432  return (*it).second;
2433  }
2434 
2435 #endif // OOMPH_HAS_MPI
2436 
2437  /// Setup boundary coordinate on boundary b.
2438  /// Boundary coordinate increases continously along
2439  /// polygonal boundary. It's zero at the lowest left
2440  /// node on the boundary.
2441  template<class ELEMENT>
2442  void setup_boundary_coordinates(const unsigned& b)
2443  {
2444  // Dummy file
2445  std::ofstream some_file;
2446  setup_boundary_coordinates<ELEMENT>(b, some_file);
2447  }
2448 
2449  /// Setup boundary coordinate on boundary b. Doc Faces
2450  /// in outfile.
2451  /// Boundary coordinate increases continously along
2452  /// polygonal boundary. It's zero at the lowest left
2453  /// node on the boundary.
2454  template<class ELEMENT>
2455  void setup_boundary_coordinates(const unsigned& b, std::ofstream& outfile);
2456 
2457 
2458  protected:
2459 #ifdef OOMPH_HAS_TRIANGLE_LIB
2460 
2461  /// Create TriangulateIO object from outer boundaries,
2462  /// internal boundaries, and open curves. Add the holes and regions
2463  /// information as well
2464  void build_triangulateio(
2465  Vector<TriangleMeshPolygon*>& outer_polygons_pt,
2466  Vector<TriangleMeshPolygon*>& internal_polygons_pt,
2467  Vector<TriangleMeshOpenCurve*>& open_curves_pt,
2468  Vector<Vector<double>>& extra_holes_coordinates,
2469  std::map<unsigned, Vector<double>>& regions_coordinates,
2470  std::map<unsigned, double>& regions_areas,
2471  TriangulateIO& triangulate_io);
2472 
2473  /// Data structure filled when the connection matrix is created, for
2474  /// each polyline, there are two vertex_connection_info structures,
2475  /// one for each end
2477  {
2482  }; // vertex_connection_info
2483 
2484  /// Data structure to store the base vertex info, initial or final
2485  /// vertex in the polylines have an associated base vertex
2487  {
2490  unsigned boundary_id;
2491  unsigned boundary_chunk;
2492  unsigned vertex_number;
2493  }; // base_vertex_info
2494 
2495  /// Helps to add information to the connection matrix of the
2496  /// given polyline
2498  TriangleMeshPolyLine* polyline_pt,
2499  std::map<unsigned, std::map<unsigned, Vector<vertex_connection_info>>>&
2500  connection_matrix,
2501  TriangleMeshPolyLine* next_polyline_pt = 0);
2502 
2503  /// Initialise the base vertex structure, set every vertex to
2504  /// no visited and not being a base vertex
2506  TriangleMeshPolyLine* polyline_pt,
2507  std::map<unsigned, std::map<unsigned, Vector<base_vertex_info>>>&
2508  base_vertices);
2509 
2510  /// Helps to identify the base vertex of the given polyline
2512  TriangleMeshPolyLine* polyline_pt,
2513  std::map<unsigned, std::map<unsigned, Vector<base_vertex_info>>>&
2514  base_vertices,
2515  std::map<unsigned, std::map<unsigned, Vector<vertex_connection_info>>>&
2516  connection_matrix,
2517  std::map<unsigned, std::map<unsigned, unsigned>>&
2518  boundary_chunk_nvertices);
2519 
2520 #endif
2521 
2522 #ifdef OOMPH_HAS_MPI
2523 
2524  /// Used to store the nodes associated to a boundary and to an
2525  /// specific segment (this only applies in distributed meshes where the
2526  /// boundary is splitted in segments)
2527  std::map<unsigned, Vector<Vector<Node*>>> Boundary_segment_node_pt;
2528 
2529  /// Stores the initial zeta coordinate for the segments that
2530  /// appear when a boundary is splited among processors
2531  std::map<unsigned, Vector<double>> Boundary_segment_initial_zeta;
2532 
2533  /// Stores the final zeta coordinate for the segments that
2534  /// appear when a boundary is splited among processors
2535  std::map<unsigned, Vector<double>> Boundary_segment_final_zeta;
2536 
2537  /// Stores the initial coordinates for the boundary
2538  std::map<unsigned, Vector<double>> Boundary_initial_coordinate;
2539 
2540  /// Stores the final coordinates for the boundary
2541  std::map<unsigned, Vector<double>> Boundary_final_coordinate;
2542 
2543  /// Stores the info. to know if it is necessary to reverse the
2544  /// segment based on a previous mesh
2545  std::map<unsigned, Vector<unsigned>> Boundary_segment_inverted;
2546 
2547  /// Stores the initial zeta coordinate for the boundary
2548  std::map<unsigned, Vector<double>> Boundary_initial_zeta_coordinate;
2549 
2550  /// Stores the final zeta coordinate for the boundary
2551  std::map<unsigned, Vector<double>> Boundary_final_zeta_coordinate;
2552 
2553  /// Stores the initial arclength for the segments that appear when
2554  /// a boundary is splited among processors
2555  std::map<unsigned, Vector<double>> Boundary_segment_initial_arclength;
2556 
2557  /// Stores the final arclength for the segments that appear when
2558  /// a boundary is splited among processors
2559  std::map<unsigned, Vector<double>> Boundary_segment_final_arclength;
2560 
2561  /// Stores the initial coordinates for the segments that appear
2562  /// when a boundary is splited among processors
2563  std::map<unsigned, Vector<Vector<double>>>
2565 
2566  /// Stores the final coordinates for the segments that appear
2567  /// when a boundary is splited among processors
2568  std::map<unsigned, Vector<Vector<double>>>
2570 
2571 #endif
2572 
2573  /// Flag to indicate whether the automatic creation of vertices
2574  /// along the boundaries by Triangle is allowed
2576 
2577  /// Snap the boundary nodes onto any curvilinear boundaries
2578  /// defined by geometric objects
2580 
2581  /// Vector of elements in each region differentiated by attribute
2582  /// (the key of the map is the attribute)
2583  std::map<unsigned, Vector<FiniteElement*>> Region_element_pt;
2584 
2585  /// Vector of attributes associated with the elements in each region
2587 
2588  /// Storage for the geometric objects associated with any boundaries
2589  std::map<unsigned, GeomObject*> Boundary_geom_object_pt;
2590 
2591  /// Storage for the limits of the boundary coordinates defined by the use
2592  /// of geometric objects. Only used for curvilinear boundaries.
2593  std::map<unsigned, Vector<double>> Boundary_coordinate_limits;
2594 
2595  /// Polygon that defines outer boundaries
2597 
2598  /// Vector of polygons that define internal polygons
2600 
2601  /// Vector of open polylines that define internal curves
2603 
2604  /// Storage for extra coordinates for holes
2606 
2607  /// Storage for extra coordinates for regions. The key on the map
2608  /// is the region id
2609  std::map<unsigned, Vector<double>> Regions_coordinates;
2610 
2611  /// A map that stores the associated curve section of the specified boundary
2612  /// id
2613  std::map<unsigned, TriangleMeshCurveSection*> Boundary_curve_section_pt;
2614 
2615  /// Storage for elements adjacent to a boundary in a particular region
2618 
2619  /// Storage for the face index adjacent to a boundary in a particular region
2621 
2622  /// Storage for pairs of doubles representing:
2623  /// .first: the arclength along the polygonal representation of
2624  /// the curviline
2625  /// .second: the corresponding intrinsic coordinate on the associated
2626  /// geometric object
2627  /// at which the vertices on the specified boundary are located.
2628  /// Only used for boundaries represented by geom objects.
2629  std::map<unsigned, Vector<std::pair<double, double>>>
2631 
2632  /// Stores a pointer to a set with all the nodes
2633  /// related with a boundary
2634  std::map<unsigned, std::set<Node*>> Nodes_on_boundary_pt;
2635 
2636  /// A set that contains the curve sections created by this object
2637  /// therefore it is necessary to free their associated memory
2638  std::set<TriangleMeshCurveSection*> Free_curve_section_pt;
2639 
2640  /// A set that contains the polygons created by this object
2641  /// therefore it is necessary to free their associated memory
2642  std::set<TriangleMeshPolygon*> Free_polygon_pt;
2643 
2644  /// A set that contains the open curves created by this
2645  /// object therefore it is necessary to free their associated memory
2646  std::set<TriangleMeshOpenCurve*> Free_open_curve_pt;
2647 
2648  /// Helper function to copy the connection information from
2649  /// the input curve(polyline or curviline) to the output polyline
2651  TriangleMeshCurveSection* output_curve_pt);
2652 
2653  /// Helper function to copy the connection information from
2654  /// the input curve(polyline or curviline) to the output sub-polyline
2656  TriangleMeshCurveSection* input_curve_pt,
2657  TriangleMeshCurveSection* output_curve_pt);
2658 
2659 
2660 #ifdef PARANOID
2661 
2662  // Used to verify if any of the polygons (closedcurves) that define
2663  // the mesh are of type ImmersedRigidBodyTriangleMeshPolygon, if
2664  // that is the case it may lead to problems in case of using load
2665  // balance
2667 
2668 #endif
2669 
2670 #ifdef OOMPH_HAS_TRIANGLE_LIB
2671 
2672  /// Helper function to create polyline vertex coordinates for
2673  /// curvilinear boundary specified by boundary_pt, using either
2674  /// equal increments in zeta or in (approximate) arclength
2675  /// along the curviline. vertex_coord[i_vertex][i_dim] stores
2676  /// i_dim-th coordinate of i_vertex-th vertex.
2677  /// polygonal_vertex_arclength_info[i_vertex] contains the pair of doubles
2678  /// made of the arclength of the i_vertex-th vertex along the polygonal
2679  /// representation (.first), and the corresponding coordinate on the
2680  /// GeomObject (.second)
2682  TriangleMeshCurviLine* boundary_pt,
2683  Vector<Vector<double>>& vertex_coord,
2684  Vector<std::pair<double, double>>& polygonal_vertex_arclength_info)
2685  {
2686  // Intrinsic coordinate along GeomObjects
2687  Vector<double> zeta(1);
2688 
2689  // Position vector to point on GeomObject
2690  Vector<double> posn(2);
2691 
2692  // Start coordinate
2693  double zeta_initial = boundary_pt->zeta_start();
2694 
2695  // How many segments do we want on this polyline?
2696  unsigned n_seg = boundary_pt->nsegment();
2697  vertex_coord.resize(n_seg + 1);
2698  polygonal_vertex_arclength_info.resize(n_seg + 1);
2699  polygonal_vertex_arclength_info[0].first = 0.0;
2700  polygonal_vertex_arclength_info[0].second = zeta_initial;
2701 
2702  // Vertices placed in equal zeta increments
2703  if (!(boundary_pt->space_vertices_evenly_in_arclength()))
2704  {
2705  // Read the values of the limiting coordinates, assuming equal
2706  // spacing of the nodes
2707  double zeta_increment =
2708  (boundary_pt->zeta_end() - boundary_pt->zeta_start()) /
2709  (double(n_seg));
2710 
2711  // Loop over the n_seg+1 points bounding the segments
2712  for (unsigned s = 0; s < n_seg + 1; s++)
2713  {
2714  // Get the coordinates
2715  zeta[0] = zeta_initial + zeta_increment * double(s);
2716  boundary_pt->geom_object_pt()->position(zeta, posn);
2717  vertex_coord[s] = posn;
2718 
2719  // Bump up the polygonal arclength
2720  if (s > 0)
2721  {
2722  polygonal_vertex_arclength_info[s].first =
2723  polygonal_vertex_arclength_info[s - 1].first +
2724  sqrt(pow(vertex_coord[s][0] - vertex_coord[s - 1][0], 2) +
2725  pow(vertex_coord[s][1] - vertex_coord[s - 1][1], 2));
2726  polygonal_vertex_arclength_info[s].second = zeta[0];
2727  }
2728  }
2729  }
2730  // Vertices placed in equal increments in (approximate) arclength
2731  else
2732  {
2733  // Number of sampling points to compute arclength and
2734  // arclength increments
2735  unsigned nsample_per_segment = 100;
2736  unsigned nsample = nsample_per_segment * n_seg;
2737 
2738  // Work out start and increment
2739  double zeta_increment =
2740  (boundary_pt->zeta_end() - boundary_pt->zeta_start()) /
2741  (double(nsample));
2742 
2743  // Get coordinate of first point
2744  Vector<double> start_point(2);
2745  zeta[0] = zeta_initial;
2746 
2747  boundary_pt->geom_object_pt()->position(zeta, start_point);
2748 
2749  // Storage for coordinates of end point
2750  Vector<double> end_point(2);
2751 
2752  // Compute total arclength
2753  double total_arclength = 0.0;
2754  for (unsigned i = 1; i < nsample; i++)
2755  {
2756  // Next point
2757  zeta[0] += zeta_increment;
2758 
2759  // Get coordinate of end point
2760  boundary_pt->geom_object_pt()->position(zeta, end_point);
2761 
2762  // Increment arclength
2763  total_arclength += sqrt(pow(end_point[0] - start_point[0], 2) +
2764  pow(end_point[1] - start_point[1], 2));
2765 
2766  // Shift back
2767  start_point = end_point;
2768  }
2769 
2770  // Desired arclength increment
2771  double target_s_increment = total_arclength / (double(n_seg));
2772 
2773  // Get coordinate of first point again
2774  zeta[0] = zeta_initial;
2775  boundary_pt->geom_object_pt()->position(zeta, start_point);
2776 
2777  // Assign as coordinate
2778  vertex_coord[0] = start_point;
2779 
2780  // Start sampling point
2781  unsigned i_lo = 1;
2782 
2783  // Loop over the n_seg-1 internal points bounding the segments
2784  for (unsigned s = 1; s < n_seg; s++)
2785  {
2786  // Visit potentially all sample points until we've found
2787  // the one at which we exceed the target arclength increment
2788  double arclength_increment = 0.0;
2789  for (unsigned i = i_lo; i < nsample; i++)
2790  {
2791  // Next point
2792  zeta[0] += zeta_increment;
2793 
2794  // Get coordinate of end point
2795  boundary_pt->geom_object_pt()->position(zeta, end_point);
2796 
2797  // Increment arclength increment
2798  arclength_increment += sqrt(pow(end_point[0] - start_point[0], 2) +
2799  pow(end_point[1] - start_point[1], 2));
2800 
2801  // Shift back
2802  start_point = end_point;
2803 
2804  // Are we there yet?
2805  if (arclength_increment > target_s_increment)
2806  {
2807  // Remember how far we've got
2808  i_lo = i;
2809 
2810  // And bail out
2811  break;
2812  }
2813  }
2814 
2815  // Store the coordinates
2816  vertex_coord[s] = end_point;
2817 
2818  // Bump up the polygonal arclength
2819  if (s > 0)
2820  {
2821  polygonal_vertex_arclength_info[s].first =
2822  polygonal_vertex_arclength_info[s - 1].first +
2823  sqrt(pow(vertex_coord[s][0] - vertex_coord[s - 1][0], 2) +
2824  pow(vertex_coord[s][1] - vertex_coord[s - 1][1], 2));
2825  polygonal_vertex_arclength_info[s].second = zeta[0];
2826  }
2827  }
2828 
2829  // Final point
2830  unsigned s = n_seg;
2831  zeta[0] = boundary_pt->zeta_end();
2832  boundary_pt->geom_object_pt()->position(zeta, end_point);
2833  vertex_coord[s] = end_point;
2834  polygonal_vertex_arclength_info[s].first =
2835  polygonal_vertex_arclength_info[s - 1].first +
2836  sqrt(pow(vertex_coord[s][0] - vertex_coord[s - 1][0], 2) +
2837  pow(vertex_coord[s][1] - vertex_coord[s - 1][1], 2));
2838  polygonal_vertex_arclength_info[s].second = zeta[0];
2839  }
2840  }
2841 
2842  /// Helper function to create polyline vertex coordinates for
2843  /// curvilinear boundary specified by boundary_pt, using either
2844  /// equal increments in zeta or in (approximate) arclength
2845  /// along the curviline. vertex_coord[i_vertex][i_dim] stores
2846  /// i_dim-th coordinate of i_vertex-th vertex.
2847  /// polygonal_vertex_arclength_info[i_vertex] contains the pair of doubles
2848  /// made of the arclength of the i_vertex-th vertex along the polygonal
2849  /// representation (.first), and the corresponding coordinate on the
2850  /// GeomObject (.second)
2852  TriangleMeshCurviLine* boundary_pt,
2853  Vector<Vector<double>>& vertex_coord,
2854  Vector<std::pair<double, double>>& polygonal_vertex_arclength_info)
2855  {
2856  // Start coordinate
2857  double zeta_initial = boundary_pt->zeta_start();
2858  // Final coordinate
2859  double zeta_final = boundary_pt->zeta_end();
2860 
2861  Vector<double>* connection_points_pt =
2862  boundary_pt->connection_points_pt();
2863 
2864  unsigned n_connections = connection_points_pt->size();
2865 
2866  // We need to sort the connection points
2867  if (n_connections > 1)
2868  {
2869  std::sort(connection_points_pt->begin(), connection_points_pt->end());
2870  }
2871 
2872 #ifdef PARANOID
2873  // Are the connection points out of range of the polyline
2874  bool out_of_range_connection_points = false;
2875  std::ostringstream error_message;
2876  // Check if the curviline should be created on a reversed way
2877  bool reversed = false;
2878  if (zeta_final < zeta_initial)
2879  {
2880  reversed = true;
2881  }
2882  if (!reversed)
2883  {
2884  if (zeta_initial > (*connection_points_pt)[0])
2885  {
2886  error_message
2887  << "One of the specified connection points is out of the\n"
2888  << "curviline limits. We found that the point ("
2889  << (*connection_points_pt)[0] << ") is\n"
2890  << "less than the"
2891  << "initial s value which is (" << zeta_initial << ").\n"
2892  << "Initial value: (" << zeta_initial << ")\n"
2893  << "Final value: (" << zeta_final << ")\n"
2894  << std::endl;
2895  out_of_range_connection_points = true;
2896  }
2897 
2898  if (zeta_final < (*connection_points_pt)[n_connections - 1])
2899  {
2900  error_message
2901  << "One of the specified connection points is out of the\n"
2902  << "curviline limits. We found that the point ("
2903  << (*connection_points_pt)[n_connections - 1] << ") is\n"
2904  << "greater than the final s value which is (" << zeta_final
2905  << ").\n"
2906  << "Initial value: (" << zeta_initial << ")\n"
2907  << "Final value: (" << zeta_final << ")\n"
2908  << std::endl;
2909  out_of_range_connection_points = true;
2910  }
2911  }
2912  else
2913  {
2914  if (zeta_initial < (*connection_points_pt)[0])
2915  {
2916  error_message
2917  << "One of the specified connection points is out of the\n"
2918  << "curviline limits. We found that the point ("
2919  << (*connection_points_pt)[0] << ") is\n"
2920  << "greater than the"
2921  << "initial s value which is (" << zeta_initial << ").\n"
2922  << "Initial value: (" << zeta_initial << ")\n"
2923  << "Final value: (" << zeta_final << ")\n"
2924  << std::endl;
2925  out_of_range_connection_points = true;
2926  }
2927 
2928  if (zeta_final > (*connection_points_pt)[n_connections - 1])
2929  {
2930  error_message
2931  << "One of the specified connection points is out of the\n"
2932  << "curviline limits. We found that the point ("
2933  << (*connection_points_pt)[n_connections - 1] << ") is\n"
2934  << "less than the final s value which is (" << zeta_final << ").\n"
2935  << "Initial value: (" << zeta_initial << ")\n"
2936  << "Final value: (" << zeta_final << ")\n"
2937  << std::endl;
2938  out_of_range_connection_points = true;
2939  }
2940  }
2941 
2942  if (out_of_range_connection_points)
2943  {
2944  throw OomphLibError(error_message.str(),
2945  OOMPH_CURRENT_FUNCTION,
2946  OOMPH_EXCEPTION_LOCATION);
2947  }
2948 
2949 #endif // PARANOID
2950 
2951  // Intrinsic coordinate along GeomObjects
2952  Vector<double> zeta(1);
2953 
2954  // Position vector to point on GeomObject
2955  Vector<double> posn(2);
2956 
2957  // How many segments do we want on this polyline?
2958  unsigned n_seg = boundary_pt->nsegment();
2959 
2960  // How many connection vertices have we already created
2961  unsigned i_connection = 0;
2962  Vector<double> zeta_connection(1);
2963 
2964  // If we have more connection points than the generated
2965  // by the number of segments then we have to change the
2966  // number of segments and create all the vertices
2967  // according to the connection points list
2968  if (n_connections >= n_seg - 1)
2969  {
2970  std::ostringstream warning_message;
2971  std::string output_string = "UnstructuredTwoDMeshGeometryBase::";
2972  output_string += "create_vertex_coordinates_for_polyline_connections()";
2973 
2974  warning_message
2975  << "The number of segments specified for the curviline with\n"
2976  << "boundary id (" << boundary_pt->boundary_id() << ") is less "
2977  << "(or equal) than the ones that will be\ngenerated by using "
2978  << "the specified number of connection points.\n"
2979  << "You specified (" << n_seg << ") segments but ("
2980  << n_connections + 1 << ") segments\nwill be generated." << std::endl;
2982  warning_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2983 
2984  // We have to explicitly change the number of segments
2985  boundary_pt->nsegment() = n_connections + 1;
2986  n_seg = boundary_pt->nsegment();
2987  vertex_coord.resize(n_seg + 1);
2988 
2989  // Initial coordinate and initial values
2990  zeta[0] = zeta_initial;
2991  boundary_pt->geom_object_pt()->position(zeta, posn);
2992  vertex_coord[0] = posn;
2993 
2994  polygonal_vertex_arclength_info.resize(n_seg + 1);
2995  polygonal_vertex_arclength_info[0].first = 0.0;
2996  polygonal_vertex_arclength_info[0].second = zeta_initial;
2997 
2998  // Loop over the n_connections points bounding the segments
2999  for (i_connection = 0; i_connection < n_connections; i_connection++)
3000  {
3001  // Get the coordinates
3002  zeta[0] = (*connection_points_pt)[i_connection];
3003  boundary_pt->geom_object_pt()->position(zeta, posn);
3004  vertex_coord[i_connection + 1] = posn;
3005 
3006  // Bump up the polygonal arclength
3007  polygonal_vertex_arclength_info[i_connection + 1].first =
3008  polygonal_vertex_arclength_info[i_connection].first +
3009  sqrt(pow(vertex_coord[i_connection + 1][0] -
3010  vertex_coord[i_connection][0],
3011  2) +
3012  pow(vertex_coord[i_connection + 1][1] -
3013  vertex_coord[i_connection][1],
3014  2));
3015  polygonal_vertex_arclength_info[i_connection + 1].second = zeta[0];
3016  }
3017 
3018  // Final coordinate and final values
3019  zeta[0] = zeta_final;
3020  boundary_pt->geom_object_pt()->position(zeta, posn);
3021  vertex_coord[n_seg] = posn;
3022 
3023  polygonal_vertex_arclength_info[n_seg].first =
3024  polygonal_vertex_arclength_info[n_seg - 1].first +
3025  sqrt(pow(vertex_coord[n_seg][0] - vertex_coord[n_seg - 1][0], 2) +
3026  pow(vertex_coord[n_seg][1] - vertex_coord[n_seg - 1][1], 2));
3027  polygonal_vertex_arclength_info[n_seg].second = zeta_final;
3028  }
3029  else
3030  {
3031  // Total number of vertices
3032  unsigned n_t_vertices = n_seg + 1;
3033 
3034  // Number of vertices left for creation
3035  unsigned l_vertices = n_t_vertices;
3036 
3037  // Total number of already created vertices
3038  unsigned n_assigned_vertices = 0;
3039 
3040  // Stores the distance between current vertices in the list
3041  // Edge vertices + Connection points - 1
3042  Vector<double> delta_z(2 + n_connections - 1);
3043 
3044  std::list<double> zeta_values_pt;
3045  zeta_values_pt.push_back(zeta_initial);
3046  for (unsigned s = 0; s < n_connections; s++)
3047  {
3048  zeta_values_pt.push_back((*connection_points_pt)[s]);
3049  }
3050  zeta_values_pt.push_back(zeta_final);
3051 
3052  l_vertices -= 2; // Edge vertices
3053  l_vertices -= n_connections; // Connection points
3054  n_assigned_vertices += 2; // Edge vertices
3055  n_assigned_vertices += n_connections; // Connection points
3056 
3057  // Vertices placed in equal zeta increments
3058  if (!(boundary_pt->space_vertices_evenly_in_arclength()))
3059  {
3060  double local_zeta_initial;
3061  double local_zeta_final;
3062  double local_zeta_increment;
3063  double local_zeta_insert;
3064 
3065  // How many vertices for each section
3066  unsigned local_n_vertices;
3067 
3068  std::list<double>::iterator l_it = zeta_values_pt.begin();
3069  std::list<double>::iterator r_it = zeta_values_pt.begin();
3070  r_it++;
3071 
3072  for (unsigned h = 0; r_it != zeta_values_pt.end();
3073  l_it++, r_it++, h++)
3074  {
3075  delta_z[h] = *r_it - *l_it;
3076  }
3077 
3078  l_it = r_it = zeta_values_pt.begin();
3079  r_it++;
3080 
3081  for (unsigned h = 0; r_it != zeta_values_pt.end(); h++)
3082  {
3083  local_n_vertices =
3084  static_cast<unsigned>(((double)n_t_vertices * delta_z[h]) /
3085  std::fabs(zeta_final - zeta_initial));
3086 
3087  local_zeta_initial = *l_it;
3088  local_zeta_final = *r_it;
3089  local_zeta_increment = (local_zeta_final - local_zeta_initial) /
3090  (double)(local_n_vertices + 1);
3091 
3092  for (unsigned s = 0; s < local_n_vertices; s++)
3093  {
3094  local_zeta_insert =
3095  local_zeta_initial + local_zeta_increment * double(s + 1);
3096  zeta_values_pt.insert(r_it, local_zeta_insert);
3097  n_assigned_vertices++;
3098  }
3099  // Moving to the next segment
3100  l_it = r_it;
3101  r_it++;
3102  }
3103 
3104  // Finishing it ...!!!
3105 #ifdef PARANOID
3106  // Counting the vertices number and the total of
3107  // assigned vertices values
3108  unsigned s = zeta_values_pt.size();
3109 
3110  if (s != n_assigned_vertices)
3111  {
3112  error_message
3113  << "The total number of assigned vertices is different from\n"
3114  << "the number of elements in the z_values list. The number"
3115  << "of\nelements in the z_values list is (" << s << ") but "
3116  << "the number\n"
3117  << "of assigned vertices is (" << n_assigned_vertices << ")."
3118  << std::endl
3119  << std::endl;
3120  throw OomphLibError(error_message.str(),
3121  OOMPH_CURRENT_FUNCTION,
3122  OOMPH_EXCEPTION_LOCATION);
3123  }
3124 #endif // PARANOID
3125 
3126  vertex_coord.resize(n_assigned_vertices);
3127  polygonal_vertex_arclength_info.resize(n_assigned_vertices);
3128  polygonal_vertex_arclength_info[0].first = 0.0;
3129  polygonal_vertex_arclength_info[0].second = zeta_initial;
3130 
3131  // Creating the vertices with the corresponding z_values
3132  l_it = zeta_values_pt.begin();
3133  for (unsigned s = 0; l_it != zeta_values_pt.end(); s++, l_it++)
3134  {
3135  // Get the coordinates
3136  zeta[0] = *l_it;
3137  boundary_pt->geom_object_pt()->position(zeta, posn);
3138  vertex_coord[s] = posn;
3139 
3140  // Bump up the polygonal arclength
3141  if (s > 0)
3142  {
3143  polygonal_vertex_arclength_info[s].first =
3144  polygonal_vertex_arclength_info[s - 1].first +
3145  sqrt(pow(vertex_coord[s][0] - vertex_coord[s - 1][0], 2) +
3146  pow(vertex_coord[s][1] - vertex_coord[s - 1][1], 2));
3147  }
3148  }
3149  }
3150  // Vertices placed in equal increments in (approximate) arclength
3151  else
3152  {
3153  // Compute the total arclength
3154  // Number of sampling points to compute arclength and
3155  // arclength increments
3156  unsigned nsample_per_segment = 100;
3157  unsigned nsample = nsample_per_segment * n_seg;
3158 
3159  // Work out start and increment
3160  double zeta_increment =
3161  (zeta_final - zeta_initial) / (double(nsample));
3162 
3163  // Get coordinate of first point
3164  Vector<double> start_point(2);
3165  zeta[0] = zeta_initial;
3166  boundary_pt->geom_object_pt()->position(zeta, start_point);
3167 
3168  // Storage for coordinates of end point
3169  Vector<double> end_point(2);
3170 
3171  // Compute total arclength
3172  double total_arclength = 0.0;
3173  for (unsigned i = 1; i < nsample; i++)
3174  {
3175  // Next point
3176  zeta[0] += zeta_increment;
3177 
3178  // Get coordinate of end point
3179  boundary_pt->geom_object_pt()->position(zeta, end_point);
3180 
3181  // Increment arclength
3182  total_arclength += sqrt(pow(end_point[0] - start_point[0], 2) +
3183  pow(end_point[1] - start_point[1], 2));
3184 
3185  // Shift back
3186  start_point = end_point;
3187  }
3188 
3189  double local_zeta_initial;
3190  double local_zeta_final;
3191  double local_zeta_increment;
3192 
3193  // How many vertices per section
3194  unsigned local_n_vertices;
3195 
3196  std::list<double>::iterator l_it = zeta_values_pt.begin();
3197  std::list<double>::iterator r_it = zeta_values_pt.begin();
3198  r_it++;
3199 
3200  for (unsigned h = 0; r_it != zeta_values_pt.end(); h++)
3201  {
3202  // There is no need to move the r_it iterator since it is
3203  // moved at the final of this loop
3204  local_zeta_initial = *l_it;
3205  local_zeta_final = *r_it;
3206  local_zeta_increment =
3207  (local_zeta_final - local_zeta_initial) / (double)(nsample);
3208 
3209  // Compute local arclength
3210  // Get coordinate of first point
3211  zeta[0] = local_zeta_initial;
3212  boundary_pt->geom_object_pt()->position(zeta, start_point);
3213 
3214  delta_z[h] = 0.0;
3215 
3216  for (unsigned i = 1; i < nsample; i++)
3217  {
3218  // Next point
3219  zeta[0] += local_zeta_increment;
3220 
3221  // Get coordinate of end point
3222  boundary_pt->geom_object_pt()->position(zeta, end_point);
3223 
3224  // Increment arclength
3225  delta_z[h] += sqrt(pow(end_point[0] - start_point[0], 2) +
3226  pow(end_point[1] - start_point[1], 2));
3227 
3228  // Shift back
3229  start_point = end_point;
3230  }
3231 
3232  local_n_vertices = static_cast<unsigned>(
3233  ((double)n_t_vertices * delta_z[h]) / (total_arclength));
3234 
3235  // Desired arclength increment
3236  double local_target_s_increment =
3237  delta_z[h] / double(local_n_vertices + 1);
3238 
3239  // Get coordinate of first point again
3240  zeta[0] = local_zeta_initial;
3241  boundary_pt->geom_object_pt()->position(zeta, start_point);
3242 
3243  // Start sampling point
3244  unsigned i_lo = 1;
3245 
3246  // Loop over the n_seg-1 internal points bounding the segments
3247  for (unsigned s = 0; s < local_n_vertices; s++)
3248  {
3249  // Visit potentially all sample points until we've found
3250  // the one at which we exceed the target arclength increment
3251  double local_arclength_increment = 0.0;
3252  for (unsigned i = i_lo; i < nsample; i++)
3253  // for (unsigned i=i_lo;i<nsample_per_segment;i++)
3254  {
3255  // Next point
3256  zeta[0] += local_zeta_increment;
3257 
3258  // Get coordinate of end point
3259  boundary_pt->geom_object_pt()->position(zeta, end_point);
3260 
3261  // Increment arclength increment
3262  local_arclength_increment +=
3263  sqrt(pow(end_point[0] - start_point[0], 2) +
3264  pow(end_point[1] - start_point[1], 2));
3265 
3266  // Shift back
3267  start_point = end_point;
3268 
3269  // Are we there yet?
3270  if (local_arclength_increment > local_target_s_increment)
3271  {
3272  // Remember how far we've got
3273  i_lo = i;
3274 
3275  // And bail out
3276  break;
3277  }
3278  }
3279 
3280  zeta_values_pt.insert(r_it, zeta[0]);
3281  n_assigned_vertices++;
3282  }
3283  // Moving to the next segments
3284  l_it = r_it;
3285  r_it++;
3286  }
3287 
3288  // Finishing it ... !!!
3289 #ifdef PARANOID
3290  // Counting the vertices number and the total of
3291  // assigned vertices values
3292  unsigned h = zeta_values_pt.size();
3293 
3294  if (h != n_assigned_vertices)
3295  {
3296  error_message
3297  << "The total number of assigned vertices is different from\n"
3298  << "the number of elements in the z_values list. The number of\n"
3299  << "elements in the z_values list is (" << h
3300  << ") but the number\n"
3301  << "of assigned vertices is (" << n_assigned_vertices << ")."
3302  << std::endl
3303  << std::endl;
3304  throw OomphLibError(error_message.str(),
3305  OOMPH_CURRENT_FUNCTION,
3306  OOMPH_EXCEPTION_LOCATION);
3307  }
3308 #endif // PARANOID
3309 
3310  vertex_coord.resize(n_assigned_vertices);
3311  polygonal_vertex_arclength_info.resize(n_assigned_vertices);
3312  polygonal_vertex_arclength_info[0].first = 0.0;
3313  polygonal_vertex_arclength_info[0].second = zeta_initial;
3314 
3315  // Creating the vertices with the corresponding z_values
3316  l_it = zeta_values_pt.begin();
3317  for (unsigned s = 0; l_it != zeta_values_pt.end(); s++, l_it++)
3318  {
3319  // Get the coordinates
3320  zeta[0] = *l_it;
3321  boundary_pt->geom_object_pt()->position(zeta, posn);
3322  vertex_coord[s] = posn;
3323 
3324  // Bump up the polygonal arclength
3325  if (s > 0)
3326  {
3327  polygonal_vertex_arclength_info[s].first =
3328  polygonal_vertex_arclength_info[s - 1].first +
3329  sqrt(pow(vertex_coord[s][0] - vertex_coord[s - 1][0], 2) +
3330  pow(vertex_coord[s][1] - vertex_coord[s - 1][1], 2));
3331  polygonal_vertex_arclength_info[s].second = zeta[0];
3332  }
3333  }
3334  } // Arclength uniformly spaced
3335  } // Less number of insertion points than vertices
3336  }
3337 
3338  /// Helper function that returns a polygon representation for
3339  /// the given closed curve, it also computes the maximum boundary id of
3340  /// the constituent curves.
3341  /// If the TriangleMeshClosedCurve is already a TriangleMeshPolygon
3342  /// we simply return a pointer to it. Otherwise a new TrilangleMeshPolygon
3343  /// is created -- this is deleted automatically when the TriangleMesh
3344  /// destructor is called, so no external book-keeping is required.
3346  TriangleMeshClosedCurve* closed_curve_pt, unsigned& max_bnd_id_local)
3347  {
3348  // How many separate boundaries do we have
3349  unsigned nb = closed_curve_pt->ncurve_section();
3350 
3351 #ifdef PARANOID
3352  if (nb < 2)
3353  {
3354  std::ostringstream error_message;
3355  error_message << "TriangleMeshClosedCurve that defines outer boundary\n"
3356  << "must be made up of at least two "
3357  << "TriangleMeshCurveSections\n"
3358  << "to allow the automatic set up boundary coordinates.\n"
3359  << "Yours only has (" << nb << ")" << std::endl;
3360  throw OomphLibError(error_message.str(),
3361  OOMPH_CURRENT_FUNCTION,
3362  OOMPH_EXCEPTION_LOCATION);
3363  }
3364 #endif
3365 
3366  // Provide storage for accompanying polylines
3367  Vector<TriangleMeshCurveSection*> my_boundary_polyline_pt(nb);
3368 
3369  // Store refinement tolerance
3370  Vector<double> refinement_tolerance(nb);
3371 
3372  // Store unrefinement tolerance
3373  Vector<double> unrefinement_tolerance(nb);
3374 
3375  // Store max. length
3376  Vector<double> max_length(nb);
3377 
3378  // Loop over boundaries that make up this boundary
3379  for (unsigned b = 0; b < nb; b++)
3380  {
3381  // Get pointer to the curve segment boundary that makes up
3382  // this part of the boundary
3383  TriangleMeshCurviLine* curviline_pt =
3384  dynamic_cast<TriangleMeshCurviLine*>(
3385  closed_curve_pt->curve_section_pt(b));
3386 
3387  TriangleMeshPolyLine* polyline_pt = dynamic_cast<TriangleMeshPolyLine*>(
3388  closed_curve_pt->curve_section_pt(b));
3389 
3390  if (curviline_pt != 0)
3391  {
3392  // Boundary id
3393  unsigned bnd_id = curviline_pt->boundary_id();
3394 
3395  // Build associated polyline
3396  my_boundary_polyline_pt[b] =
3397  curviline_to_polyline(curviline_pt, bnd_id);
3398 
3399  // Copy the unrefinement tolerance
3400  unrefinement_tolerance[b] = curviline_pt->unrefinement_tolerance();
3401 
3402  // Copy the refinement tolerance
3403  refinement_tolerance[b] = curviline_pt->refinement_tolerance();
3404 
3405  // Copy the maximum length
3406  max_length[b] = curviline_pt->maximum_length();
3407 
3408  // Updates bnd_id<--->curve section map
3409  Boundary_curve_section_pt[bnd_id] = my_boundary_polyline_pt[b];
3410 
3411  // Keep track of curve sections that need to be deleted!!!
3412  Free_curve_section_pt.insert(my_boundary_polyline_pt[b]);
3413 
3414  // Keep track...
3415  if (bnd_id > max_bnd_id_local)
3416  {
3417  max_bnd_id_local = bnd_id;
3418  }
3419  }
3420  else if (polyline_pt != 0)
3421  {
3422  // Boundary id
3423  unsigned bnd_id = polyline_pt->boundary_id();
3424 
3425  // Pass the pointer of the already existing polyline
3426  my_boundary_polyline_pt[b] = polyline_pt;
3427 
3428  // Copy the unrefinement tolerance
3429  unrefinement_tolerance[b] = polyline_pt->unrefinement_tolerance();
3430 
3431  // Copy the refinement tolerance
3432  refinement_tolerance[b] = polyline_pt->refinement_tolerance();
3433 
3434  // Copy the maximum length
3435  max_length[b] = polyline_pt->maximum_length();
3436 
3437  // Updates bnd_id<--->curve section map
3438  Boundary_curve_section_pt[bnd_id] = my_boundary_polyline_pt[b];
3439 
3440  // Keep track...
3441  if (bnd_id > max_bnd_id_local)
3442  {
3443  max_bnd_id_local = bnd_id;
3444  }
3445  }
3446  else
3447  {
3448  std::ostringstream error_stream;
3449  error_stream << "The 'curve_segment' is not a curviline neither a\n "
3450  << "polyline: What is it?\n"
3451  << std::endl;
3452  throw OomphLibError(error_stream.str(),
3453  OOMPH_CURRENT_FUNCTION,
3454  OOMPH_EXCEPTION_LOCATION);
3455  }
3456 
3457  } // end of loop over boundaries
3458 
3459  // Create a new polygon by using the new created polylines
3460  TriangleMeshPolygon* output_polygon_pt =
3461  new TriangleMeshPolygon(my_boundary_polyline_pt,
3462  closed_curve_pt->internal_point(),
3463  closed_curve_pt->is_internal_point_fixed());
3464 
3465  // Keep track of new created polygons that need to be deleted!!!
3466  Free_polygon_pt.insert(output_polygon_pt);
3467 
3468  // Pass on refinement information
3469  output_polygon_pt->set_polyline_refinement_tolerance(
3470  closed_curve_pt->polyline_refinement_tolerance());
3471  output_polygon_pt->set_polyline_unrefinement_tolerance(
3472  closed_curve_pt->polyline_unrefinement_tolerance());
3473 
3474  // Loop over boundaries that make up this boundary and copy
3475  // refinement, unrefinement and max length information
3476  for (unsigned b = 0; b < nb; b++)
3477  {
3478  // Set the unrefinement and refinement information
3479  my_boundary_polyline_pt[b]->set_unrefinement_tolerance(
3480  unrefinement_tolerance[b]);
3481 
3482  my_boundary_polyline_pt[b]->set_refinement_tolerance(
3483  refinement_tolerance[b]);
3484 
3485  // Copy the maximum length constraint
3486  my_boundary_polyline_pt[b]->set_maximum_length(max_length[b]);
3487  }
3488  return output_polygon_pt;
3489  }
3490 
3491  /// Helper function that creates and returns an open curve with
3492  /// the polyline representation of its constituent curve sections. The
3493  /// new created open curve is deleted when the TriangleMesh destructor
3494  /// is called
3496  TriangleMeshOpenCurve* open_curve_pt, unsigned& max_bnd_id_local)
3497  {
3498  unsigned nb = open_curve_pt->ncurve_section();
3499 
3500  // Provide storage for accompanying polylines
3501  Vector<TriangleMeshCurveSection*> my_boundary_polyline_pt(nb);
3502 
3503  // Store refinement tolerance
3504  Vector<double> refinement_tolerance(nb);
3505 
3506  // Store unrefinement tolerance
3507  Vector<double> unrefinement_tolerance(nb);
3508 
3509  // Store max. length
3510  Vector<double> max_length(nb);
3511 
3512  // Loop over the number of curve sections on the open curve
3513  for (unsigned i = 0; i < nb; i++)
3514  {
3515  // Get pointer to the curve segment boundary that makes up
3516  // this part of the boundary
3517  TriangleMeshCurviLine* curviline_pt =
3518  dynamic_cast<TriangleMeshCurviLine*>(
3519  open_curve_pt->curve_section_pt(i));
3520  TriangleMeshPolyLine* polyline_pt = dynamic_cast<TriangleMeshPolyLine*>(
3521  open_curve_pt->curve_section_pt(i));
3522 
3523  if (curviline_pt != 0)
3524  {
3525  // Boundary id
3526  unsigned bnd_id = curviline_pt->boundary_id();
3527 
3528  // Build associated polyline
3529  my_boundary_polyline_pt[i] =
3530  curviline_to_polyline(curviline_pt, bnd_id);
3531 
3532  // Copy the unrefinement tolerance
3533  unrefinement_tolerance[i] = curviline_pt->unrefinement_tolerance();
3534 
3535  // Copy the refinement tolerance
3536  refinement_tolerance[i] = curviline_pt->refinement_tolerance();
3537 
3538  // Copy the maximum length
3539  max_length[i] = curviline_pt->maximum_length();
3540 
3541  // Pass the connection information to the polyline representation
3542  copy_connection_information(curviline_pt, my_boundary_polyline_pt[i]);
3543 
3544  // Updates bnd_id<--->curve section map
3545  Boundary_curve_section_pt[bnd_id] = my_boundary_polyline_pt[i];
3546 
3547  // Keep track of curve sections that need to be deleted!!!
3548  Free_curve_section_pt.insert(my_boundary_polyline_pt[i]);
3549 
3550  // Keep track...
3551  if (bnd_id > max_bnd_id_local)
3552  {
3553  max_bnd_id_local = bnd_id;
3554  }
3555  }
3556  else if (polyline_pt != 0)
3557  {
3558  // Boundary id
3559  unsigned bnd_id = polyline_pt->boundary_id();
3560 
3561  // Storage pointer
3562  my_boundary_polyline_pt[i] = polyline_pt;
3563 
3564  // Copy the unrefinement tolerance
3565  unrefinement_tolerance[i] = polyline_pt->unrefinement_tolerance();
3566 
3567  // Copy the refinement tolerance
3568  refinement_tolerance[i] = polyline_pt->refinement_tolerance();
3569 
3570  // Copy the maximum length
3571  max_length[i] = polyline_pt->maximum_length();
3572 
3573  // Pass the connection information to the polyline representation
3574  copy_connection_information(polyline_pt, my_boundary_polyline_pt[i]);
3575 
3576  // Updates bnd_id<--->curve section map
3577  Boundary_curve_section_pt[bnd_id] = my_boundary_polyline_pt[i];
3578 
3579  // Keep track...
3580  if (bnd_id > max_bnd_id_local)
3581  {
3582  max_bnd_id_local = bnd_id;
3583  }
3584  }
3585  else
3586  {
3587  std::ostringstream error_stream;
3588  error_stream
3589  << "The 'curve_segment' (open) is not a curviline neither a\n "
3590  << "polyline: What is it?\n"
3591  << std::endl;
3592  throw OomphLibError(error_stream.str(),
3593  OOMPH_CURRENT_FUNCTION,
3594  OOMPH_EXCEPTION_LOCATION);
3595  }
3596  } // end of loop over boundaries
3597 
3598  // Create open curve with polylines boundaries
3599  TriangleMeshOpenCurve* output_open_polyline_pt =
3600  new TriangleMeshOpenCurve(my_boundary_polyline_pt);
3601 
3602  // Keep track of open polylines that need to be deleted!!!
3603  Free_open_curve_pt.insert(output_open_polyline_pt);
3604 
3605  // Pass on refinement information
3606  output_open_polyline_pt->set_polyline_refinement_tolerance(
3607  open_curve_pt->polyline_refinement_tolerance());
3608  output_open_polyline_pt->set_polyline_unrefinement_tolerance(
3609  open_curve_pt->polyline_unrefinement_tolerance());
3610 
3611  // Loop over boundaries that make up this boundary and copy
3612  // refinement, unrefinement and max length information
3613  for (unsigned b = 0; b < nb; b++)
3614  {
3615  // Set the unrefinement and refinement information
3616  my_boundary_polyline_pt[b]->set_unrefinement_tolerance(
3617  unrefinement_tolerance[b]);
3618 
3619  my_boundary_polyline_pt[b]->set_refinement_tolerance(
3620  refinement_tolerance[b]);
3621 
3622  // Copy the maximum length constraint
3623  my_boundary_polyline_pt[b]->set_maximum_length(max_length[b]);
3624  }
3625  return output_open_polyline_pt;
3626  }
3627 
3628  /// Stores the geometric objects associated to the
3629  /// curve sections that compound the closed curve. It also
3630  /// stores the limits defined by these geometric objects
3632  TriangleMeshClosedCurve* input_closed_curve_pt)
3633  {
3634  unsigned nb = input_closed_curve_pt->ncurve_section();
3635 
3636 #ifdef PARANOID
3637 
3638  if (nb < 2)
3639  {
3640  std::ostringstream error_message;
3641  error_message << "TriangleMeshCurve that defines closed boundary\n"
3642  << "must be made up of at least two "
3643  << "TriangleMeshCurveSection\n"
3644  << "to allow the automatic set up boundary coordinates.\n"
3645  << "Yours only has " << nb << std::endl;
3646  throw OomphLibError(error_message.str(),
3647  OOMPH_CURRENT_FUNCTION,
3648  OOMPH_EXCEPTION_LOCATION);
3649  }
3650 
3651 #endif
3652 
3653  // TODO: Used for the ImmersedRigidBodyTriangleMeshPolygon objects only
3654  // ImmersedRigidBodyTriangleMeshPolygon* bound_geom_obj_pt
3655  //= dynamic_cast<ImmersedRigidBodyTriangleMeshPolygon*>
3656  // (input_closed_curve_pt);
3657  GeomObject* bound_geom_obj_pt =
3658  dynamic_cast<GeomObject*>(input_closed_curve_pt);
3659 
3660  // If cast successful set up the coordinates
3661  if (bound_geom_obj_pt != 0)
3662  {
3663  unsigned n_poly = input_closed_curve_pt->ncurve_section();
3664  for (unsigned p = 0; p < n_poly; p++)
3665  {
3666  // Read out the index of the boundary from the polyline
3667  unsigned b_index =
3668  input_closed_curve_pt->curve_section_pt(p)->boundary_id();
3669 
3670  // Set the geometric object
3671  Boundary_geom_object_pt[b_index] = bound_geom_obj_pt;
3672 
3673  // The coordinates along each polyline boundary are scaled to
3674  // of unit length so the total coordinate limits are simply
3675  // (p,p+1)
3676  Boundary_coordinate_limits[b_index].resize(2);
3677  Boundary_coordinate_limits[b_index][0] = p;
3678  Boundary_coordinate_limits[b_index][1] = p + 1.0;
3679  }
3680 
3681 #ifdef PARANOID
3682  // If we are using parallel mesh adaptation and load balancing,
3683  // we are going to need to check for the use of this type of
3684  // Polygon at this stage, so switch on the flag
3686 #endif
3687  }
3688  else
3689  {
3690  // Loop over curve sections that make up this boundary
3691  for (unsigned b = 0; b < nb; b++)
3692  {
3693  TriangleMeshCurviLine* curviline_pt =
3694  dynamic_cast<TriangleMeshCurviLine*>(
3695  input_closed_curve_pt->curve_section_pt(b));
3696 
3697  if (curviline_pt != 0)
3698  {
3699  // Read the values of the limiting coordinates
3700  Vector<double> zeta_bound(2);
3701  zeta_bound[0] = curviline_pt->zeta_start();
3702  zeta_bound[1] = curviline_pt->zeta_end();
3703 
3704  // Boundary id
3705  unsigned bnd_id = curviline_pt->boundary_id();
3706 
3707  // Set the boundary geometric object and limits
3708  Boundary_geom_object_pt[bnd_id] = curviline_pt->geom_object_pt();
3709  Boundary_coordinate_limits[bnd_id] = zeta_bound;
3710  }
3711  } // for
3712  } // else
3713  } // function
3714 
3715  /// Stores the geometric objects associated to the
3716  /// curve sections that compound the open curve. It also
3717  /// stores the limits defined by these geometric objects
3719  TriangleMeshOpenCurve* input_open_curve_pt)
3720  {
3721  unsigned nb = input_open_curve_pt->ncurve_section();
3722 
3723  // Loop over curve sections that make up this boundary
3724  for (unsigned b = 0; b < nb; b++)
3725  {
3726  TriangleMeshCurviLine* curviline_pt =
3727  dynamic_cast<TriangleMeshCurviLine*>(
3728  input_open_curve_pt->curve_section_pt(b));
3729 
3730  if (curviline_pt != 0)
3731  {
3732  // ead the values of the limiting coordinates
3733  Vector<double> zeta_bound(2);
3734  zeta_bound[0] = curviline_pt->zeta_start();
3735  zeta_bound[1] = curviline_pt->zeta_end();
3736 
3737  // Boundary id
3738  unsigned bnd_id = curviline_pt->boundary_id();
3739 
3740  // Set the boundary geometric object and limits
3741  Boundary_geom_object_pt[bnd_id] = curviline_pt->geom_object_pt();
3742  Boundary_coordinate_limits[bnd_id] = zeta_bound;
3743  }
3744  } // for
3745  } // function
3746 
3747 #endif // OOMPH_HAS_TRIANGLE_LIB
3748 
3749  private:
3750 #ifdef OOMPH_HAS_TRIANGLE_LIB
3751 
3752  /// Helper function that creates the associated polyline
3753  /// representation for curvilines
3755  TriangleMeshCurviLine*& curviline_pt, unsigned& bnd_id)
3756  {
3757  // Create vertex coordinates for polygonal representation
3758  Vector<Vector<double>> bound;
3759  Vector<std::pair<double, double>> polygonal_vertex_arclength;
3760 
3761  if (curviline_pt->are_there_connection_points())
3762  {
3764  curviline_pt, bound, polygonal_vertex_arclength);
3765  }
3766  else
3767  {
3769  curviline_pt, bound, polygonal_vertex_arclength);
3770  }
3771 
3772  // Store the vertex-arclength information
3773  Polygonal_vertex_arclength_info[bnd_id] = polygonal_vertex_arclength;
3774 
3775  // Build associated polyline
3776  return new TriangleMeshPolyLine(bound, bnd_id);
3777  }
3778 
3779  /// Get the associated vertex to the given s value by looking to
3780  /// the list of s values created when changing from curviline to polyline
3781  unsigned get_associated_vertex_to_svalue(double& target_s_value,
3782  unsigned& bnd_id)
3783  {
3784  double s_tolerance = 1.0e-14;
3786  target_s_value, bnd_id, s_tolerance);
3787  }
3788 
3789  /// Get the associated vertex to the given s value by looking to
3790  /// the list of s values created when changing from curviline to polyline
3791  unsigned get_associated_vertex_to_svalue(double& target_s_value,
3792  unsigned& bnd_id,
3793  double& s_tolerance)
3794  {
3795  // Create a pointer to the list of s coordinates and arclength values
3796  // associated with a vertex
3797  Vector<std::pair<double, double>>* vertex_info =
3799 
3800  // Total vertex number
3801  unsigned vector_size = vertex_info->size();
3802 
3803  // Counter for current vertex number
3804  unsigned n_vertex = 0;
3805 
3806  // Find the associated value to the given s value
3807  do
3808  {
3809  // Store the current zeta value
3810  double s = (*vertex_info)[n_vertex].second;
3811 
3812  // When find it save the vertex number and return it
3813  if (std::fabs(s - target_s_value) < s_tolerance)
3814  {
3815  break;
3816  }
3817 
3818  // Increment n_vertex
3819  n_vertex++;
3820  } while (n_vertex < vector_size);
3821 
3822 #ifdef PARANOID
3823 
3824  if (n_vertex >= vector_size)
3825  {
3826  std::ostringstream error_message;
3827  error_message << "Could not find the associated vertex number in\n"
3828  << "boundary " << bnd_id << " with the given s\n"
3829  << "connection value (" << target_s_value << ") using\n"
3830  << "this tolerance: " << s_tolerance << std::endl;
3831  throw OomphLibError(error_message.str(),
3832  OOMPH_CURRENT_FUNCTION,
3833  OOMPH_EXCEPTION_LOCATION);
3834  }
3835 
3836 #endif
3837  return n_vertex;
3838  }
3839 
3840 #endif // OOMPH_HAS_TRIANGLE_LIB
3841  };
3842 
3843 
3844  //======================================================================
3845  /// Setup boundary coordinate on boundary b. Doc Faces
3846  /// in outfile. Boundary coordinate increases continously along
3847  /// polygonal boundary. It's zero at the lexicographically
3848  /// smallest node on the boundary.
3849  //======================================================================
3850  template<class ELEMENT>
3852  const unsigned& b, std::ofstream& outfile)
3853  {
3854  // Temporary storage for face elements
3855  Vector<FiniteElement*> face_el_pt;
3856 
3857  // Temporary storage for number of elements adjacent to the boundary
3858  unsigned nel = 0;
3859 
3860  // =================================================================
3861  // BEGIN: Get face elements from boundary elements
3862  // =================================================================
3863 
3864  // Temporary storage for elements adjacent to the boundary that have
3865  // an common edge (related with internal boundaries)
3866  unsigned n_repeated_ele = 0;
3867 
3868  unsigned n_regions = this->nregion();
3869 
3870 #ifdef OOMPH_HAS_MPI
3871  // map to associate the face element to the bulk element, this info.
3872  // is only necessary for the setup of boundary coordinates in a
3873  // distributed mesh when we need to extract the halo/haloed info.
3874  std::map<FiniteElement*, FiniteElement*> face_to_bulk_element_pt;
3875 #endif
3876 
3877  // Temporary storage for already done nodes
3878  Vector<std::pair<Node*, Node*>> done_nodes_pt;
3879 
3880  // If there is more than one region then only use boundary
3881  // coordinates from the bulk side (region 0)
3882  if (n_regions > 1)
3883  {
3884  for (unsigned rr = 0; rr < n_regions; rr++)
3885  {
3886  unsigned region_id = static_cast<unsigned>(this->Region_attribute[rr]);
3887 
3888 #ifdef PARANOID
3889  double diff =
3890  fabs(Region_attribute[rr] - static_cast<double>(static_cast<unsigned>(
3891  this->Region_attribute[rr])));
3892  if (diff > 0.0)
3893  {
3894  std::ostringstream error_message;
3895  error_message << "Region attributes should be unsigneds because we \n"
3896  << "only use them to set region ids\n";
3897  throw OomphLibError(error_message.str(),
3898  OOMPH_CURRENT_FUNCTION,
3899  OOMPH_EXCEPTION_LOCATION);
3900  }
3901 #endif
3902 
3903  // Loop over all elements on boundaries in region rr
3904  unsigned nel_in_region =
3905  this->nboundary_element_in_region(b, region_id);
3906  unsigned nel_repeated_in_region = 0;
3907 
3908 #ifdef PARANOID
3910  {
3911  if (nel_in_region == 0)
3912  {
3913  std::ostringstream warning_message;
3914  std::string output_string =
3915  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
3916  warning_message
3917  << "There are no elements associated with boundary (" << b
3918  << ")\n"
3919  << "in region (" << region_id << "). This could happen because:\n"
3920  << "1) You did not specify boundaries with this boundary id.\n"
3921  << "---- Review carefully the indexing of your boundaries.\n"
3922  << "2) The boundary (" << b << ") is not associated with region ("
3923  << region_id << ").\n"
3924  << "---- The boundary does not touch the region.\n"
3925  << "You can suppress this warning by setting the static public "
3926  "bool\n\n"
3927  << " "
3928  "UnstructuredTwoDMeshGeometryBase::Suppress_warning_about_"
3929  "regions_and_boundaries\n\n"
3930  << "to true.\n";
3932  warning_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3933  }
3934  }
3935 #endif
3936 
3937  // Only bother to do anything else, if there are elements
3938  // associated with the boundary and the current region
3939  if (nel_in_region > 0)
3940  {
3941  // Flag that activates when a repeated face element is found,
3942  // possibly because we are dealing with an internal boundary
3943  bool repeated = false;
3944 
3945  // Loop over the bulk elements adjacent to boundary b
3946  for (unsigned e = 0; e < nel_in_region; e++)
3947  {
3948  // Get pointer to the bulk element that is adjacent to boundary b
3949  FiniteElement* bulk_elem_pt =
3950  this->boundary_element_in_region_pt(b, region_id, e);
3951 
3952 #ifdef OOMPH_HAS_MPI
3953  // In a distributed mesh only work with nonhalo elements
3954  if (this->is_mesh_distributed() && bulk_elem_pt->is_halo())
3955  {
3956  // Increase the number of repeated elements
3957  n_repeated_ele++;
3958  // Skip this element and go for the next one
3959  continue;
3960  }
3961 #endif
3962 
3963  // Find the index of the face of element e along boundary b
3964  int face_index =
3965  this->face_index_at_boundary_in_region(b, region_id, e);
3966 
3967  // Before adding the new element we need to be sure that
3968  // the edge that this element represent has not been
3969  // already added
3970  FiniteElement* tmp_ele_pt =
3971  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
3972 
3973  const unsigned n_nodes = tmp_ele_pt->nnode();
3974 
3975  std::pair<Node*, Node*> tmp_pair = std::make_pair(
3976  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
3977 
3978  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
3979  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
3980 
3981  // Search for repeated nodes
3982  const unsigned n_done_nodes = done_nodes_pt.size();
3983  for (unsigned l = 0; l < n_done_nodes; l++)
3984  {
3985  if (tmp_pair == done_nodes_pt[l] ||
3986  tmp_pair_inverse == done_nodes_pt[l])
3987  {
3988  nel_repeated_in_region++;
3989  repeated = true;
3990  break;
3991  }
3992  }
3993 
3994  // Create new face element?
3995  if (!repeated)
3996  {
3997  // Add the pair of nodes (edge) to the node dones
3998  done_nodes_pt.push_back(tmp_pair);
3999 #ifdef OOMPH_HAS_MPI
4000  // If the mesh is distributed then create a map from the
4001  // temporary face element to the bulk element, further
4002  // info. will be extracted from the bulk element for setup
4003  // of boundary coordinates in a distributed mesh
4004  if (this->is_mesh_distributed())
4005  {
4006  face_to_bulk_element_pt[tmp_ele_pt] = bulk_elem_pt;
4007  }
4008 #endif
4009  // Add the face element to the storage
4010  face_el_pt.push_back(tmp_ele_pt);
4011  }
4012  else
4013  {
4014  // Clean up
4015  delete tmp_ele_pt;
4016  tmp_ele_pt = 0;
4017  }
4018 
4019  // Re-start
4020  repeated = false;
4021 
4022  // Output faces?
4023  if (outfile.is_open())
4024  {
4025  face_el_pt[face_el_pt.size() - 1]->output(outfile);
4026  }
4027  } // for nel
4028 
4029  nel += nel_in_region;
4030 
4031  n_repeated_ele += nel_repeated_in_region;
4032 
4033  } // if (nel_in_region > 0)
4034 
4035  } // for (rr < n_regions)
4036 
4037  } // if (n_regions > 1)
4038  // Otherwise it's just the normal boundary functions
4039  else
4040  {
4041  // Loop over all elements on boundaries
4042  nel = this->nboundary_element(b);
4043 
4044 #ifdef PARANOID
4046  {
4047  if (nel == 0)
4048  {
4049  std::ostringstream warning_message;
4050  std::string output_string =
4051  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4052  warning_message
4053  << "There are no elements associated with boundary (" << b << ").\n"
4054  << "This could happen because you did not specify boundaries with\n"
4055  << "this boundary id. Review carefully the indexing of your\n"
4056  << "boundaries.";
4058  warning_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4059  }
4060  }
4061 #endif
4062 
4063  // Only bother to do anything else, if there are elements
4064  if (nel > 0)
4065  {
4066  // Flag that activates when a repeated face element is found,
4067  // possibly because we are dealing with an internal boundary
4068  bool repeated = false;
4069 
4070  // Loop over the bulk elements adjacent to boundary b
4071  for (unsigned e = 0; e < nel; e++)
4072  {
4073  // Get pointer to the bulk element that is adjacent to boundary b
4074  FiniteElement* bulk_elem_pt = this->boundary_element_pt(b, e);
4075 
4076 #ifdef OOMPH_HAS_MPI
4077 
4078  // In a distributed mesh only work with nonhalo elements
4079  if (this->is_mesh_distributed() && bulk_elem_pt->is_halo())
4080  {
4081  // Increase the number of repeated elements
4082  n_repeated_ele++;
4083  // Skip this element and go for the next one
4084  continue;
4085  }
4086 
4087 #endif
4088 
4089  // Find the index of the face of element e along boundary b
4090  int face_index = this->face_index_at_boundary(b, e);
4091 
4092  // Before adding the new element we need to be sure that the
4093  // edge that this element represent has not been already added
4094  FiniteElement* tmp_ele_pt =
4095  new DummyFaceElement<ELEMENT>(bulk_elem_pt, face_index);
4096 
4097  const unsigned n_nodes = tmp_ele_pt->nnode();
4098 
4099  std::pair<Node*, Node*> tmp_pair = std::make_pair(
4100  tmp_ele_pt->node_pt(0), tmp_ele_pt->node_pt(n_nodes - 1));
4101 
4102  std::pair<Node*, Node*> tmp_pair_inverse = std::make_pair(
4103  tmp_ele_pt->node_pt(n_nodes - 1), tmp_ele_pt->node_pt(0));
4104 
4105  // Search for repeated nodes
4106  const unsigned n_done_nodes = done_nodes_pt.size();
4107  for (unsigned l = 0; l < n_done_nodes; l++)
4108  {
4109  if (tmp_pair == done_nodes_pt[l] ||
4110  tmp_pair_inverse == done_nodes_pt[l])
4111  {
4112  n_repeated_ele++;
4113  repeated = true;
4114  break;
4115  }
4116  }
4117 
4118  // Create new face element
4119  if (!repeated)
4120  {
4121  // Add the pair of nodes (edge) to the node dones
4122  done_nodes_pt.push_back(tmp_pair);
4123 #ifdef OOMPH_HAS_MPI
4124  // Create a map from the temporary face element to the bulk
4125  // element, further info. will be extracted from the bulk
4126  // element for setup of boundary coordinates in a
4127  // distributed mesh
4128  if (this->is_mesh_distributed())
4129  {
4130  face_to_bulk_element_pt[tmp_ele_pt] = bulk_elem_pt;
4131  }
4132 #endif
4133  face_el_pt.push_back(tmp_ele_pt);
4134  }
4135  else
4136  {
4137  // Free the repeated bulk element!!
4138  delete tmp_ele_pt;
4139  tmp_ele_pt = 0;
4140  }
4141 
4142  // Re-start
4143  repeated = false;
4144 
4145  // Output faces?
4146  if (outfile.is_open())
4147  {
4148  face_el_pt[face_el_pt.size() - 1]->output(outfile);
4149  }
4150 
4151  } // for (e < nel)
4152 
4153  } // if (nel > 0)
4154 
4155  } // else if (n_regions > 1)
4156 
4157  // Do not consider the repeated elements
4158  nel -= n_repeated_ele;
4159 
4160 #ifdef PARANOID
4161  if (nel != face_el_pt.size())
4162  {
4163  std::ostringstream error_message;
4164  error_message
4165  << "The independent counting of face elements (" << nel << ") for "
4166  << "boundary (" << b << ") is different\n"
4167  << "from the real number of face elements in the container ("
4168  << face_el_pt.size() << ")\n";
4169  throw OomphLibError(
4170  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
4171  }
4172 #endif
4173 
4174  // =================================================================
4175  // END: Get face elements from boundary elements
4176  // =================================================================
4177 
4178  // Only bother to do anything else, if there are elements
4179  if (nel > 0)
4180  {
4181  // A flag vector to mark those face elements that are considered
4182  // as halo in the current processor
4183  std::vector<bool> is_halo_face_element(nel, false);
4184 
4185  // Count the total number of non halo face elements
4186  unsigned nnon_halo_face_elements = 0;
4187 
4188 #ifdef OOMPH_HAS_MPI
4189  // Only mark the face elements as halo if the mesh is marked as
4190  // distributed
4191  if (this->is_mesh_distributed())
4192  {
4193  for (unsigned ie = 0; ie < nel; ie++)
4194  {
4195  FiniteElement* face_ele_pt = face_el_pt[ie];
4196  // Get the bulk element
4197  FiniteElement* tmp_bulk_ele_pt = face_to_bulk_element_pt[face_ele_pt];
4198  // Check if the bulk element is halo
4199  if (!tmp_bulk_ele_pt->is_halo())
4200  {
4201  // Mark the face element as nonhalo
4202  is_halo_face_element[ie] = false;
4203  // Increase the counter for nonhalo elements
4204  nnon_halo_face_elements++;
4205  }
4206  else
4207  {
4208  // Mark the face element as halo
4209  is_halo_face_element[ie] = true;
4210  }
4211  } // for (ie < nel)
4212  } // if (this->is_mesh_distributed())
4213  else
4214  {
4215 #endif // OOMPH_HAS_MPI
4216 
4217  // If the mesh is not distributed then the number of non halo
4218  // elements is the same as the number of elements
4219  nnon_halo_face_elements = nel;
4220 
4221 #ifdef OOMPH_HAS_MPI
4222  } // else if (this->is_mesh_distributed())
4223 #endif
4224 
4225 #ifdef PARANOID
4226  // Get the total number of halo face elements
4227  const unsigned nhalo_face_element = nel - nnon_halo_face_elements;
4228 
4229  if (nhalo_face_element > 0)
4230  {
4231  std::ostringstream error_message;
4232  error_message
4233  << "There should not be halo face elements since they were not\n"
4234  << "considered when computing the face elements.\n"
4235  << "The number of found halo face elements is: ("
4236  << nhalo_face_element << ").\n\n";
4237  throw OomphLibError(error_message.str(),
4238  OOMPH_CURRENT_FUNCTION,
4239  OOMPH_EXCEPTION_LOCATION);
4240  }
4241 #endif
4242 
4243  // =================================================================
4244  // BEGIN: Sort face elements
4245  // =================================================================
4246 
4247  // The vector of lists to store the "segments" that compound the
4248  // boundaries (segments may appear only in a distributed mesh
4249  // because the boundary may have been split across multiple
4250  // processors)
4251  Vector<std::list<FiniteElement*>> segment_sorted_ele_pt;
4252 
4253  // Number of already sorted face elements (only nonhalo face
4254  // elements for a distributed mesh)
4255  unsigned nsorted_face_elements = 0;
4256 
4257  // Keep track of who's done (in a distributed mesh this apply to
4258  // nonhalo only)
4259  std::map<FiniteElement*, bool> done_el;
4260 
4261  // Keep track of which element is inverted (in distributed mesh
4262  // the elements may be inverted with respect to the segment they
4263  // belong)
4264  std::map<FiniteElement*, bool> is_inverted;
4265 
4266  // Iterate until all possible segments have been created. In a
4267  // non distributed mesh there is only one segment which defines
4268  // the complete boundary
4269  while (nsorted_face_elements < nnon_halo_face_elements)
4270  {
4271  // The sorted list of face elements (in a distributed mesh a
4272  // collection of continuous face elements define a segment)
4273  std::list<FiniteElement*> sorted_el_pt;
4274 
4275  FiniteElement* ele_face_pt = 0;
4276 
4277 #ifdef PARANOID
4278  // Select an initial element for the segment
4279  bool found_initial_face_element = false;
4280 #endif
4281 
4282  // Store the index of the initial face element
4283  unsigned iface = 0;
4284 #ifdef OOMPH_HAS_MPI
4285  if (this->is_mesh_distributed())
4286  {
4287  for (iface = 0; iface < nel; iface++)
4288  {
4289  // Only work with nonhalo face elements
4290  if (!is_halo_face_element[iface])
4291  {
4292  ele_face_pt = face_el_pt[iface];
4293  // If not done then take it as initial face element
4294  if (!done_el[ele_face_pt])
4295  {
4296 #ifdef PARANOID
4297  // Set the flag to indicate the initial element was
4298  // found
4299  found_initial_face_element = true;
4300 #endif
4301  // Increase the number of sorted face elements
4302  nsorted_face_elements++;
4303  // Set the index to the next face element
4304  iface++;
4305  // Add the face element in the container
4306  sorted_el_pt.push_back(ele_face_pt);
4307  // Mark as done
4308  done_el[ele_face_pt] = true;
4309  break;
4310  } // if (!done_el[ele_face_pt])
4311  } // if (!is_halo_face_element[iface])
4312  } // for (iface < nel)
4313  } // if (this->is_mesh_distributed())
4314  else
4315  {
4316 #endif // #ifdef OOMPH_HAS_MPI
4317 
4318  // When the mesh is not distributed just take the first
4319  // element and put it in the sorted list
4320  ele_face_pt = face_el_pt[0];
4321 #ifdef PARANOID
4322  // Set the flag to indicate the initial element was found
4323  found_initial_face_element = true;
4324 #endif
4325  // Increase the number of sorted face elements
4326  nsorted_face_elements++;
4327  // Set the index to the next face element
4328  iface = 1;
4329  // Add the face element in the container
4330  sorted_el_pt.push_back(ele_face_pt);
4331  // Mark as done
4332  done_el[ele_face_pt] = true;
4333 
4334 #ifdef OOMPH_HAS_MPI
4335  } // else if (this->is_mesh_distributed())
4336 #endif
4337 
4338 #ifdef PARANOID
4339  if (!found_initial_face_element)
4340  {
4341  std::ostringstream error_message;
4342  error_message << "Could not find an initial face element for the "
4343  "current segment\n";
4344  throw OomphLibError(error_message.str(),
4345  OOMPH_CURRENT_FUNCTION,
4346  OOMPH_EXCEPTION_LOCATION);
4347  }
4348 #endif
4349 
4350  // Number of nodes of the initial face element
4351  const unsigned nnod = ele_face_pt->nnode();
4352 
4353  // Left and rightmost nodes (the left and right nodes of the
4354  // current face element)
4355  Node* left_node_pt = ele_face_pt->node_pt(0);
4356  Node* right_node_pt = ele_face_pt->node_pt(nnod - 1);
4357 
4358  // Continue iterating if a new face element has been added to
4359  // the list
4360  bool face_element_added = false;
4361 
4362  // While a new face element has been added to the set of sorted
4363  // face elements continue iterating
4364  do
4365  {
4366  // Start from the next face element since we have already
4367  // added the previous one as the initial face element (any
4368  // previous face element had to be added on previous
4369  // iterations)
4370  for (unsigned iiface = iface; iiface < nel; iiface++)
4371  {
4372  // Re-start flag
4373  face_element_added = false;
4374 
4375  // Get the candidate element
4376  ele_face_pt = face_el_pt[iiface];
4377 
4378  // Check that the candidate element has not been done and
4379  // is not a halo element
4380  if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
4381  {
4382  // Get the left and right nodes of the current element
4383  Node* local_left_node_pt = ele_face_pt->node_pt(0);
4384  Node* local_right_node_pt = ele_face_pt->node_pt(nnod - 1);
4385 
4386  // New element fits at the left of segment and is not inverted
4387  if (left_node_pt == local_right_node_pt)
4388  {
4389  left_node_pt = local_left_node_pt;
4390  sorted_el_pt.push_front(ele_face_pt);
4391  is_inverted[ele_face_pt] = false;
4392  face_element_added = true;
4393  }
4394  // New element fits at the left of segment and is inverted
4395  else if (left_node_pt == local_left_node_pt)
4396  {
4397  left_node_pt = local_right_node_pt;
4398  sorted_el_pt.push_front(ele_face_pt);
4399  is_inverted[ele_face_pt] = true;
4400  face_element_added = true;
4401  }
4402  // New element fits on the right of segment and is not inverted
4403  else if (right_node_pt == local_left_node_pt)
4404  {
4405  right_node_pt = local_right_node_pt;
4406  sorted_el_pt.push_back(ele_face_pt);
4407  is_inverted[ele_face_pt] = false;
4408  face_element_added = true;
4409  }
4410  // New element fits on the right of segment and is inverted
4411  else if (right_node_pt == local_right_node_pt)
4412  {
4413  right_node_pt = local_left_node_pt;
4414  sorted_el_pt.push_back(ele_face_pt);
4415  is_inverted[ele_face_pt] = true;
4416  face_element_added = true;
4417  }
4418 
4419  if (face_element_added)
4420  {
4421  done_el[ele_face_pt] = true;
4422  nsorted_face_elements++;
4423  break;
4424  }
4425 
4426  } // if (!(done_el[ele_face_pt] || is_halo_face_element[iiface]))
4427  } // for (iiface<nnon_halo_face_element)
4428  } while (face_element_added &&
4429  (nsorted_face_elements < nnon_halo_face_elements));
4430 
4431  // Store the created segment in the vector of segments
4432  segment_sorted_ele_pt.push_back(sorted_el_pt);
4433 
4434  } // while(nsorted_face_elements < nnon_halo_face_elements);
4435 
4436 #ifdef OOMPH_HAS_MPI
4437  if (!this->is_mesh_distributed())
4438  {
4439 #endif
4440  // Are we done?
4441  if (nsorted_face_elements != nel || segment_sorted_ele_pt.size() != 1)
4442  {
4443  std::ostringstream error_message;
4444  error_message << "Was only able to setup boundary coordinate on "
4445  << "boundary " << b << "\nfor " << nsorted_face_elements
4446  << " of " << nel
4447  << " face elements. This usually means\n"
4448  << "that the boundary is not simply connected.\n\n"
4449  << "Re-run the setup_boundary_coordintes() function\n"
4450  << "with an output file specified "
4451  << "as the second argument.\n"
4452  << "This file will contain FaceElements that\n"
4453  << "oomph-lib believes to be located on the boundary.\n"
4454  << std::endl;
4455  throw OomphLibError(error_message.str(),
4456  OOMPH_CURRENT_FUNCTION,
4457  OOMPH_EXCEPTION_LOCATION);
4458  }
4459 #ifdef OOMPH_HAS_MPI
4460  } // if (!this->is_mesh_distributed())
4461 #endif
4462 
4463  // =================================================================
4464  // END: Sort face elements
4465  // =================================================================
4466 
4467  // ----------------------------------------------------------------
4468 
4469  // =================================================================
4470  // BEGIN: Assign global/local (non distributed mesh/distributed
4471  // mesh) boundary coordinates to nodes
4472  // =================================================================
4473 
4474  // Compute the (local) boundary coordinates of the nodes in the
4475  // segments. In a distributed mesh this info. will be used by a
4476  // root processor to compute the (global) boundary coordinates.
4477 
4478  // Vector of sets that stores the nodes of each segment based on
4479  // a lexicographically order starting from the bottom left node
4480  // of each segment
4481  Vector<std::set<Node*>> segment_all_nodes_pt;
4482 
4483  // The number of segments in this processor
4484  const unsigned nsegments = segment_sorted_ele_pt.size();
4485 
4486 #ifdef PARANOID
4487  if (nnon_halo_face_elements > 0 && nsegments == 0)
4488  {
4489  std::ostringstream error_message;
4490  error_message
4491  << "The number of segments is zero, but the number of nonhalo\n"
4492  << "elements is: (" << nnon_halo_face_elements << ")\n";
4493  throw OomphLibError(error_message.str(),
4494  OOMPH_CURRENT_FUNCTION,
4495  OOMPH_EXCEPTION_LOCATION);
4496  } // if (nnon_halo_face_elements > 0 && nsegments == 0)
4497 
4498 #endif
4499 
4500  // Store the arclength of each segment in the current processor
4501  Vector<double> segment_arclength(nsegments);
4502 
4503 #ifdef OOMPH_HAS_MPI
4504 
4505  // Clear the boundary segment nodes storage
4506  this->flush_boundary_segment_node(b);
4507 
4508  // Set the number of segments for the current boundary
4509  this->set_nboundary_segment_node(b, nsegments);
4510 
4511 #endif // #ifdef OOMPH_HAS_MPI
4512 
4513  // Go through all the segments and compute their (local) boundary
4514  // coordinates
4515  for (unsigned is = 0; is < nsegments; is++)
4516  {
4517 #ifdef PARANOID
4518 
4519  if (segment_sorted_ele_pt[is].size() == 0)
4520  {
4521  std::ostringstream error_message;
4522  std::string output_string =
4523  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4524  error_message << "The (" << is << ")-th segment has no elements\n";
4525  throw OomphLibError(
4526  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4527  } // if (segment_sorted_ele_pt[is].size() == 0)
4528 
4529 #endif
4530 
4531  // Get access to the first element on the segment
4532  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
4533 
4534  // Number of nodes
4535  unsigned nnod = first_ele_pt->nnode();
4536 
4537  // Get the first node of the current segment
4538  Node* first_node_pt = first_ele_pt->node_pt(0);
4539  if (is_inverted[first_ele_pt])
4540  {
4541  first_node_pt = first_ele_pt->node_pt(nnod - 1);
4542  }
4543 
4544  // Coordinates of left node
4545  double x_left = first_node_pt->x(0);
4546  double y_left = first_node_pt->x(1);
4547 
4548  // Initialise boundary coordinate (local boundary coordinate
4549  // for boundaries with more than one segment)
4550  Vector<double> zeta(1, 0.0);
4551 
4552  // Set boundary coordinate
4553  first_node_pt->set_coordinates_on_boundary(b, zeta);
4554 
4555  // Lexicographically bottom left node
4556  std::set<Node*> local_nodes_pt;
4557  local_nodes_pt.insert(first_node_pt);
4558 
4559 #ifdef OOMPH_HAS_MPI
4560 
4561  // Insert the node in the look-up scheme for nodes in segments
4562  this->add_boundary_segment_node(b, is, first_node_pt);
4563 
4564 #endif // #ifdef OOMPH_HAS_MPI
4565 
4566  // Now loop over nodes in order
4567  for (std::list<FiniteElement*>::iterator it =
4568  segment_sorted_ele_pt[is].begin();
4569  it != segment_sorted_ele_pt[is].end();
4570  it++)
4571  {
4572  // Get element
4573  FiniteElement* el_pt = *it;
4574 
4575  // Start node and increment
4576  unsigned k_nod = 1;
4577  int nod_diff = 1;
4578  if (is_inverted[el_pt])
4579  {
4580  k_nod = nnod - 2;
4581  nod_diff = -1;
4582  }
4583 
4584  // Loop over nodes
4585  for (unsigned j = 1; j < nnod; j++)
4586  {
4587  Node* nod_pt = el_pt->node_pt(k_nod);
4588  k_nod += nod_diff;
4589 
4590  // Coordinates of right node
4591  double x_right = nod_pt->x(0);
4592  double y_right = nod_pt->x(1);
4593 
4594  // Increment boundary coordinate
4595  zeta[0] += sqrt((x_right - x_left) * (x_right - x_left) +
4596  (y_right - y_left) * (y_right - y_left));
4597 
4598  // Set boundary coordinate
4599  nod_pt->set_coordinates_on_boundary(b, zeta);
4600 
4601  // Increment reference coordinate
4602  x_left = x_right;
4603  y_left = y_right;
4604 
4605  // Get lexicographically bottom left node but only
4606  // use vertex nodes as candidates
4607  local_nodes_pt.insert(nod_pt);
4608 
4609 #ifdef OOMPH_HAS_MPI
4610 
4611  // Insert the node in the look-up scheme for nodes in segments
4612  this->add_boundary_segment_node(b, is, nod_pt);
4613 
4614 #endif // #ifdef OOMPH_HAS_MPI
4615 
4616  } // for (j < nnod)
4617  } // iterator over the elements in the segment
4618 
4619  // Assign the arclength of the current segment
4620  segment_arclength[is] = zeta[0];
4621 
4622  // Add the nodes for the corresponding segment in the container
4623  segment_all_nodes_pt.push_back(local_nodes_pt);
4624 
4625  } // for (is < nsegments)
4626 
4627  // =================================================================
4628  // END: Assign global/local (non distributed mesh/distributed
4629  // mesh) boundary coordinates to nodes
4630  // =================================================================
4631 
4632  // Store the initial arclength for each segment of boundary in
4633  // the current processor, initialise to zero in case we have a
4634  // non distributed mesh
4635  Vector<double> initial_segment_arclength(nsegments, 0.0);
4636 
4637  // Store the final arclength for each segment of boundary in the
4638  // current processor, initalise to zero in case we have a non
4639  // distributed mesh
4640  Vector<double> final_segment_arclength(nsegments, 0.0);
4641 
4642  // Store the initial zeta for each segment of boundary in the
4643  // current processor, initalise to zero in case we have a non
4644  // distributed mesh
4645  Vector<double> initial_segment_zeta(nsegments, 0.0);
4646 
4647  // Store the final zeta for each segment of boundary in the
4648  // current processor, initalise to zero in case we have a non
4649  // distributed mesh
4650  Vector<double> final_segment_zeta(nsegments, 0.0);
4651 
4652  // --------------------------------------------------------------
4653  // DISTRIBUTED MESH: BEGIN - Check orientation of boundaries and
4654  // assign boundary coordinates accordingly
4655  // --------------------------------------------------------------
4656 
4657 #ifdef OOMPH_HAS_MPI
4658 
4659  // Check that the mesh is distributed and that the initial zeta
4660  // values for the boundary segments have been already
4661  // assigned. When the method is called by the first time (when
4662  // the mesh is just created) the initial zeta values for the
4663  // boundary segments are not available
4664  if (this->is_mesh_distributed() &&
4666  {
4667  // Get the initial and final coordinates of each segment
4668 
4669  // For each segment in the processor check whether it was
4670  // inverted in the root processor, if that is the case then it
4671  // is necessary to re-sort the face elements and invert them
4672  for (unsigned is = 0; is < nsegments; is++)
4673  {
4674  // Check if we need/or not to invert the current zeta values
4675  // on the segment (only apply for boundaries with GeomObject
4676  // associated)
4677 
4678  // Does the boundary HAS a GeomObject associated?
4679  if (this->boundary_geom_object_pt(b) != 0)
4680  {
4681  // Get the first and last node of the current segment and
4682  // their zeta values
4683 
4684  // Get access to the first element on the segment
4685  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
4686 
4687  // Number of nodes
4688  const unsigned nnod = first_ele_pt->nnode();
4689 
4690  // Get the first node of the current segment
4691  Node* first_node_pt = first_ele_pt->node_pt(0);
4692  if (is_inverted[first_ele_pt])
4693  {
4694  first_node_pt = first_ele_pt->node_pt(nnod - 1);
4695  }
4696 
4697  // Get access to the last element on the segment
4698  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
4699 
4700  // Get the last node of the current segment
4701  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
4702  if (is_inverted[last_ele_pt])
4703  {
4704  last_node_pt = last_ele_pt->node_pt(0);
4705  }
4706 
4707  // Get the zeta coordinates for the first and last node
4708  Vector<double> current_segment_initial_zeta(1);
4709  Vector<double> current_segment_final_zeta(1);
4710 
4711  first_node_pt->get_coordinates_on_boundary(
4712  b, current_segment_initial_zeta);
4713  last_node_pt->get_coordinates_on_boundary(
4714  b, current_segment_final_zeta);
4715 
4716 #ifdef PARANOID
4717  // Check whether the zeta values in the segment are in
4718  // increasing or decreasing order
4719  if (current_segment_initial_zeta[0] < current_segment_final_zeta[0])
4720  {
4721  // do nothing
4722  }
4723  else if (current_segment_initial_zeta[0] >
4724  current_segment_final_zeta[0])
4725  {
4726  std::stringstream error_message;
4727  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
4728  "setup_boundary_coordinates()";
4729  error_message
4730  << "The zeta values are in decreasing order, this is weird\n"
4731  << "since they have just been set-up in increasing order few\n"
4732  << "lines above\n"
4733  << "Boundary: (" << b << ")\n"
4734  << "Segment: (" << is << ")\n"
4735  << "Initial zeta value: (" << current_segment_initial_zeta[0]
4736  << ")\n"
4737  << "Initial coordinate: (" << first_node_pt->x(0) << ", "
4738  << first_node_pt->x(1) << ")\n"
4739  << "Final zeta value: (" << current_segment_final_zeta[0]
4740  << ")\n"
4741  << "Initial coordinate: (" << last_node_pt->x(0) << ", "
4742  << last_node_pt->x(1) << ")\n";
4743  throw OomphLibError(
4744  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4745  }
4746  else
4747  {
4748  std::stringstream error_message;
4749  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
4750  "setup_boundary_coordinates()";
4751  error_message
4752  << "It was not possible to determine whether the zeta values on"
4753  << "the current segment\nof the boundary are in"
4754  << "increasing or decreasing order\n\n"
4755  << "Boundary: (" << b << ")\n"
4756  << "Segment: (" << is << ")\n"
4757  << "Arclength: (" << segment_arclength[is] << "\n"
4758  << "Initial zeta value: (" << current_segment_initial_zeta[0]
4759  << ")\n"
4760  << "Initial coordinate: (" << first_node_pt->x(0) << ", "
4761  << first_node_pt->x(1) << ")\n"
4762  << "Final zeta value: (" << current_segment_final_zeta[0]
4763  << ")\n"
4764  << "Initial coordinate: (" << last_node_pt->x(0) << ", "
4765  << last_node_pt->x(1) << ")\n";
4766  throw OomphLibError(
4767  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4768  }
4769 #endif // #ifdef PARANOID
4770 
4771  // Now get the original zeta values and check if they are in
4772  // increasing or decreasing order
4773  const double original_segment_initial_zeta =
4775  const double original_segment_final_zeta =
4777 
4778  // Now check if the zeta values go in increase or decrease
4779  // order
4780  if (original_segment_final_zeta > original_segment_initial_zeta)
4781  {
4782  // The original values go in increasing order, only need
4783  // to change the values if the original segment is marked
4784  // as inverted
4785  if (boundary_segment_inverted(b)[is])
4786  {
4787  // The original segment is inverted, then we need to
4788  // reverse the boundary coordinates. Go through all the
4789  // nodes and reverse its value
4790  std::set<Node*> all_nodes_pt = segment_all_nodes_pt[is];
4791  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
4792  it != all_nodes_pt.end();
4793  it++)
4794  {
4795  Vector<double> zeta(1);
4796  // Get the node
4797  Node* nod_pt = (*it);
4798  // Get the boundary coordinate associated to the node
4799  nod_pt->get_coordinates_on_boundary(b, zeta);
4800  // Compute its new value
4801  zeta[0] = segment_arclength[is] - zeta[0];
4802  // Set the new boundary coordinate value
4803  nod_pt->set_coordinates_on_boundary(b, zeta);
4804  } // Loop over the nodes in the segment
4805  } // if (boundary_segment_inverted(b)[is])
4806  else
4807  {
4808  // The boundary is not inverted, do nothing!!!
4809  }
4810 
4811  } // original zeta values in increasing order
4812  else if (original_segment_final_zeta <
4813  original_segment_initial_zeta)
4814  {
4815  // The original values go in decreasing order, only need
4816  // to change the values if the original segment is NOT
4817  // marked as inverted
4818  if (boundary_segment_inverted(b)[is])
4819  {
4820  // The boundary is inverted, do nothing!!!
4821  }
4822  else
4823  {
4824  // The original segment is NOT inverted, then we need
4825  // to reverse the boundary coordinates. Go through all
4826  // the nodes and reverse its value
4827  std::set<Node*> all_nodes_pt = segment_all_nodes_pt[is];
4828  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
4829  it != all_nodes_pt.end();
4830  it++)
4831  {
4832  Vector<double> zeta(1);
4833  // Get the node
4834  Node* nod_pt = (*it);
4835  // Get the boundary coordinate associated to the node
4836  nod_pt->get_coordinates_on_boundary(b, zeta);
4837  // Compute its new value
4838  zeta[0] = segment_arclength[is] - zeta[0];
4839  // Set the new boundary coordinate value
4840  nod_pt->set_coordinates_on_boundary(b, zeta);
4841  } // Loop over the nodes in the segment
4842  } // else if (boundary_segment_inverted(b)[is])
4843  } // original zeta values in decreasing order
4844 #ifdef PARANOID
4845  else
4846  {
4847  std::stringstream error_message;
4848  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
4849  "setup_boundary_coordinates()";
4850  error_message
4851  << "It was not possible to identify if the zeta values on the\n"
4852  << "current segment in the boundary should go in increasing\n"
4853  << "or decreasing order.\n\n"
4854  << "Boundary: (" << b << ")\n"
4855  << "Segment: (" << is << ")\n"
4856  << "Initial zeta value: (" << original_segment_initial_zeta
4857  << ")\n"
4858  << "Final zeta value: (" << original_segment_final_zeta
4859  << ")\n";
4860  throw OomphLibError(
4861  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4862  }
4863 #endif
4864 
4865  } // if (this->boundary_geom_object_pt(b)!=0)
4866  else
4867  {
4868  // No GeomObject associated, do nothing!!!
4869 
4870  } // else if (this->boundary_geom_object_pt(b)!=0)
4871  } // for (is < nsegments)
4872  } // if (this->is_mesh_distributed() &&
4873  // Assigned_segments_initial_zeta_values[b])
4874 
4875 #endif // #ifdef OOMPH_HAS_MPI
4876 
4877  // Get the number of sets for nodes
4878 #ifdef PARANOID
4879  if (segment_all_nodes_pt.size() != nsegments)
4880  {
4881  std::ostringstream error_message;
4882  std::string output_string =
4883  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4884  error_message << "The number of segments (" << nsegments
4885  << ") and the number of "
4886  << "sets of nodes (" << segment_all_nodes_pt.size()
4887  << ") representing\n"
4888  << "the\nsegments is different!!!\n\n";
4889  throw OomphLibError(
4890  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4891  }
4892 #endif
4893 
4894  // --------------------------------------------------------------
4895  // DISTRIBUTED MESH: END - Check orientation of boundaries and
4896  // assign boundary coordinates accordingly
4897  // --------------------------------------------------------------
4898 
4899  // =================================================================
4900  // BEGIN: Get the lenght of the boundaries or segments of the
4901  // boundary
4902  // =================================================================
4903 
4904  // The nodes have been assigned arc-length coordinates from one
4905  // end or the other of the segments. If the mesh is distributed
4906  // the values are set so that they agree with the increasing or
4907  // decreasing order of the zeta values for the segments
4908 
4909  // Storage for the coordinates of the first and last nodes
4910  Vector<double> first_coordinate(2);
4911  Vector<double> last_coordinate(2);
4912  // Storage for the zeta coordinate of the first and last nodes
4913  Vector<double> first_node_zeta_coordinate(1, 0.0);
4914  Vector<double> last_node_zeta_coordinate(1, 0.0);
4915 
4916  // Store the accumulated arclength, used at the end to check the
4917  // max arclength of the boundary
4918  double boundary_arclength = 0.0;
4919 
4920  // If the mesh is marked as distributed and the initial zeta
4921  // values for the segments have been computed then get the
4922  // info. regarding the initial and final nodes coordinates, same
4923  // as the zeta boundary values for those nodes
4924 
4925 #ifdef OOMPH_HAS_MPI
4926  if (this->is_mesh_distributed() &&
4928  {
4929  // --------------------------------------------------------------
4930  // DISTRIBUTED MESH: BEGIN
4931  // --------------------------------------------------------------
4932 
4933  // Get the initial and final coordinates for the complete boundary
4934  first_coordinate = boundary_initial_coordinate(b);
4935  last_coordinate = boundary_final_coordinate(b);
4936  // Get the initial and final zeta values for the boundary
4937  // (arclength)
4938  first_node_zeta_coordinate = boundary_initial_zeta_coordinate(b);
4939  last_node_zeta_coordinate = boundary_final_zeta_coordinate(b);
4940 
4941  // The total arclength is the maximum between the initial and
4942  // final zeta coordinate
4943  boundary_arclength =
4944  std::max(first_node_zeta_coordinate[0], last_node_zeta_coordinate[0]);
4945 
4946 #ifdef PARANOID
4947  if (boundary_arclength == 0)
4948  {
4949  std::ostringstream error_message;
4950  std::string output_string =
4951  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
4952  error_message << "The boundary arclength is zero for boundary (" << b
4953  << ")\n";
4954  throw OomphLibError(
4955  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
4956  }
4957 #endif
4958 
4959  // Check if there is a GeomObject associated to the boundary,
4960  // and assign the corresponding zeta (arclength) values
4961  if (this->boundary_geom_object_pt(b) != 0)
4962  {
4963  initial_segment_zeta = boundary_segment_initial_zeta(b);
4964  final_segment_zeta = boundary_segment_final_zeta(b);
4965  }
4966  else
4967  {
4968  initial_segment_arclength = boundary_segment_initial_arclength(b);
4969  final_segment_arclength = boundary_segment_final_arclength(b);
4970  }
4971 
4972  // --------------------------------------------------------------
4973  // DISTRIBUTED MESH: END
4974  // --------------------------------------------------------------
4975 
4976  } // if (this->is_mesh_distributed() &&
4977  // Assigned_segments_initial_zeta_values[b])
4978  else
4979  {
4980 #endif // #ifdef OOMPH_HAS_MPI
4981 
4982  // If the mesh is NOT distributed or the initial and final zeta
4983  // values of the segments have not been assigned then perform
4984  // as in the serial case
4985  if (this->boundary_geom_object_pt(b) != 0)
4986  {
4987  initial_segment_zeta[0] = this->boundary_coordinate_limits(b)[0];
4988  final_segment_zeta[0] = this->boundary_coordinate_limits(b)[1];
4989  }
4990  else
4991  {
4992  // When the mesh is or not distributed but the initial
4993  // segment's zeta values HAVE NOT been established then
4994  // initalize the initial segment to zero
4995  initial_segment_arclength[0] = 0.0;
4996  }
4997 #ifdef OOMPH_HAS_MPI
4998  } // The mesh is NOT distributed or the zeta values for the
4999  // segments have not been established
5000 
5001 #endif // #ifdef OOMPH_HAS_MPI
5002 
5003  // =================================================================
5004  // END: Get the lenght of the boundaries or segments of the
5005  // boundary
5006  // =================================================================
5007 
5008  // =================================================================
5009  // BEGIN: Scale the boundary coordinates based on whether the
5010  // boundary has an associated GeomObj or not
5011  // =================================================================
5012 
5013  // Go through all the segments and assign the scaled boundary
5014  // coordinates
5015  for (unsigned is = 0; is < nsegments; is++)
5016  {
5017 #ifdef PARANOID
5018  if (segment_sorted_ele_pt[is].size() == 0)
5019  {
5020  std::ostringstream error_message;
5021  std::string output_string =
5022  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
5023  error_message << "The (" << is << ")-th segment has no elements\n";
5024  throw OomphLibError(
5025  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5026  } // if (segment_sorted_ele_pt[is].size() == 0)x
5027 #endif
5028 
5029  // Get the first and last nodes coordinates for the current
5030  // segment
5031 
5032 #ifdef OOMPH_HAS_MPI
5033  // Check if the initial zeta values of the segments have been
5034  // established, if they have not then get the first and last
5035  // coordinates from the current segments, same as the zeta
5036  // values
5038  {
5039 #endif // #ifdef OOMPH_HAS_MPI
5040 
5041  // Get access to the first element on the segment
5042  FiniteElement* first_ele_pt = segment_sorted_ele_pt[is].front();
5043 
5044  // Number of nodes
5045  const unsigned nnod = first_ele_pt->nnode();
5046 
5047  // Get the first node of the current segment
5048  Node* first_node_pt = first_ele_pt->node_pt(0);
5049  if (is_inverted[first_ele_pt])
5050  {
5051  first_node_pt = first_ele_pt->node_pt(nnod - 1);
5052  }
5053 
5054  // Get access to the last element on the segment
5055  FiniteElement* last_ele_pt = segment_sorted_ele_pt[is].back();
5056 
5057  // Get the last node of the current segment
5058  Node* last_node_pt = last_ele_pt->node_pt(nnod - 1);
5059  if (is_inverted[last_ele_pt])
5060  {
5061  last_node_pt = last_ele_pt->node_pt(0);
5062  }
5063 
5064  // Get the coordinates for the first and last node
5065  for (unsigned i = 0; i < 2; i++)
5066  {
5067  first_coordinate[i] = first_node_pt->x(i);
5068  last_coordinate[i] = last_node_pt->x(i);
5069  }
5070 
5071  // Get the zeta coordinates for the first and last node
5072  first_node_pt->get_coordinates_on_boundary(
5073  b, first_node_zeta_coordinate);
5074  last_node_pt->get_coordinates_on_boundary(b,
5075  last_node_zeta_coordinate);
5076 
5077 #ifdef OOMPH_HAS_MPI
5078  } // if (!Assigned_segments_initial_zeta_values[b])
5079 #endif // #ifdef OOMPH_HAS_MPI
5080 
5081  // Get the nodes of the current segment
5082  std::set<Node*> all_nodes_pt = segment_all_nodes_pt[is];
5083 
5084  // If the boundary has a geometric object representation then
5085  // scale the coordinates to match those of the geometric object
5086  GeomObject* const geom_object_pt = this->boundary_geom_object_pt(b);
5087  if (geom_object_pt != 0)
5088  {
5089  Vector<double> bound_coord_limits =
5090  this->boundary_coordinate_limits(b);
5091  // Get the position of the ends of the geometric object
5092  Vector<double> zeta(1);
5093  Vector<double> first_geom_object_location(2);
5094  Vector<double> last_geom_object_location(2);
5095 
5096  // Get the zeta value for the initial coordinates of the
5097  // GeomObject
5098  zeta[0] = bound_coord_limits[0];
5099  // zeta[0] = initial_segment_zeta[is];
5100 
5101  // Get the coordinates from the initial zeta value from the
5102  // GeomObject
5103  geom_object_pt->position(zeta, first_geom_object_location);
5104 
5105  // Get the zeta value for the final coordinates of the
5106  // GeomObject
5107  zeta[0] = bound_coord_limits[1];
5108  // zeta[0] = final_segment_zeta[is];
5109 
5110  // Get the coordinates from the final zeta value from the
5111  // GeomObject
5112  geom_object_pt->position(zeta, last_geom_object_location);
5113 
5114  // Calculate the errors in position between the first and
5115  // last nodes and the endpoints of the geometric object
5116  double error = 0.0;
5117  double tmp_error = 0.0;
5118  for (unsigned i = 0; i < 2; i++)
5119  {
5120  const double dist =
5121  first_geom_object_location[i] - first_coordinate[i];
5122  tmp_error += dist * dist;
5123  }
5124  error += sqrt(tmp_error);
5125  tmp_error = 0.0;
5126  for (unsigned i = 0; i < 2; i++)
5127  {
5128  const double dist =
5129  last_geom_object_location[i] - last_coordinate[i];
5130  tmp_error += dist * dist;
5131  }
5132  error += sqrt(tmp_error);
5133 
5134  // Calculate the errors in position between the first and
5135  // last nodes and the endpoints of the geometric object if
5136  // reversed
5137  double rev_error = 0.0;
5138  tmp_error = 0.0;
5139  for (unsigned i = 0; i < 2; i++)
5140  {
5141  const double dist =
5142  first_geom_object_location[i] - last_coordinate[i];
5143  tmp_error += dist * dist;
5144  }
5145  rev_error += sqrt(tmp_error);
5146  tmp_error = 0.0;
5147  for (unsigned i = 0; i < 2; i++)
5148  {
5149  const double dist =
5150  last_geom_object_location[i] - first_coordinate[i];
5151  tmp_error += dist * dist;
5152  }
5153  rev_error += sqrt(tmp_error);
5154 
5155  // Number of polyline vertices along this boundary
5156  const unsigned n_vertex = Polygonal_vertex_arclength_info[b].size();
5157 
5158  // Get polygonal vertex data
5159  Vector<std::pair<double, double>> polygonal_vertex_arclength =
5161 
5162  // If the (normal) error is small than reversed then we have
5163  // the coordinate direction correct. If not then we must
5164  // reverse it bool reversed = false;
5165  if (error < rev_error)
5166  {
5167  // Coordinates are aligned (btw: don't delete this block --
5168  // there's a final else below to catch errors!)
5169 
5170  // reversed = false;
5171  }
5172  else if (error > rev_error)
5173  {
5174  // reversed = true;
5175 
5176  // Reverse the limits of the boundary coordinates along the
5177  // geometric object
5178  double temp = bound_coord_limits[0];
5179  bound_coord_limits[0] = bound_coord_limits[1];
5180  bound_coord_limits[1] = temp;
5181 
5182 #ifdef OOMPH_HAS_MPI
5183  // If we are working with a NON distributed mesh then
5184  // re-assign the initial and final zeta values
5185  if (!this->is_mesh_distributed())
5186  {
5187 #endif
5188  temp = initial_segment_zeta[is];
5189  initial_segment_zeta[is] = final_segment_zeta[is];
5190  final_segment_zeta[is] = temp;
5191 #ifdef OOMPH_HAS_MPI
5192  }
5193 #endif
5194  // Reverse the vertices information
5195  for (unsigned v = 0; v < n_vertex; v++)
5196  {
5197  polygonal_vertex_arclength[v].first =
5198  Polygonal_vertex_arclength_info[b][v].first;
5199 
5200  polygonal_vertex_arclength[v].second =
5201  Polygonal_vertex_arclength_info[b][n_vertex - v - 1].second;
5202  } // for (v < n_vertex)
5203  }
5204  else
5205  {
5206  std::ostringstream error_stream;
5207  std::string output_string =
5208  "UnstructuredTwoDMeshGeometryBase::setup_boundary_coordinates()";
5209  error_stream
5210  << "Something very strange has happened.\n"
5211  << "The error between the endpoints of the geometric object\n"
5212  << "and the first and last nodes on the boundary is the same\n"
5213  << "irrespective of the direction of the coordinate.\n"
5214  << "This probably means that things are way off.\n"
5215  << "The errors are " << error << " and " << rev_error << "\n";
5216  std::cout << error_stream.str();
5217  throw OomphLibError(
5218  error_stream.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5219  }
5220 
5221  // Get the total arclength of the edge
5222  // last_node_pt->get_coordinates_on_boundary(b, zeta);
5223  // double zeta_old_range=zeta[0];
5224 
5225  // Store the arclength of the segment (not yet mapped to the
5226  // boundary coordinates defined by the GeomObject)
5227  const double zeta_old_range = segment_arclength[is];
5228 
5229  // double zeta_new_range=bound_coord_limits[1]-bound_coord_limits[0];
5230 
5231  // The range to map the zeta values
5232  double zeta_new_range =
5233  final_segment_zeta[is] - initial_segment_zeta[is];
5234  // The initial zeta value for the current segment
5235  double initial_local_segment_zeta = initial_segment_zeta[is];
5236 
5237 #ifdef OOMPH_HAS_MPI
5238 
5239  // If we are working with a distributed mes and the initial
5240  // and final zeta values for the segments have been
5241  // established then reset the range to map
5242  if (this->is_mesh_distributed() &&
5244  {
5245  // Re-set the range to map the zeta values
5246  zeta_new_range =
5247  std::fabs(final_segment_zeta[is] - initial_segment_zeta[is]);
5248 
5249  // Re-set the initial zeta value of the segment
5250  initial_local_segment_zeta =
5251  std::min(initial_segment_zeta[is], final_segment_zeta[is]);
5252  }
5253 
5254 #endif
5255 
5256  // Re-assign boundary coordinate for the case where boundary
5257  // is represented by polygon
5258  unsigned use_old = false;
5259  if (n_vertex == 0) use_old = true;
5260 
5261  // Now scale the coordinates accordingly
5262  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
5263  it != all_nodes_pt.end();
5264  it++)
5265  {
5266  // Get the node
5267  Node* nod_pt = (*it);
5268 
5269  // Get coordinate based on arclength along polygonal repesentation
5270  nod_pt->get_coordinates_on_boundary(b, zeta);
5271 
5272  if (use_old)
5273  {
5274  // Boundary is actually a polygon -- simply rescale
5275  zeta[0] = initial_local_segment_zeta +
5276  (zeta_new_range / zeta_old_range) * (zeta[0]);
5277  }
5278  else
5279  {
5280  // Scale such that vertex nodes stay where they were
5281  bool found = false;
5282 
5283  // Loop over vertex nodes
5284  for (unsigned v = 1; v < n_vertex; v++)
5285  {
5286  if ((zeta[0] >= polygonal_vertex_arclength[v - 1].first) &&
5287  (zeta[0] <= polygonal_vertex_arclength[v].first))
5288  {
5289  // Increment in intrinsic coordinate along geom object
5290  double delta_zeta =
5291  (polygonal_vertex_arclength[v].second -
5292  polygonal_vertex_arclength[v - 1].second);
5293  // Increment in arclength along segment
5294  double delta_polyarc =
5295  (polygonal_vertex_arclength[v].first -
5296  polygonal_vertex_arclength[v - 1].first);
5297 
5298  // Mapped arclength coordinate
5299  double zeta_new =
5300  polygonal_vertex_arclength[v - 1].second +
5301  delta_zeta *
5302  (zeta[0] - polygonal_vertex_arclength[v - 1].first) /
5303  delta_polyarc;
5304  zeta[0] = zeta_new;
5305 
5306  // Success!
5307  found = true;
5308 
5309  // Bail out
5310  break;
5311  }
5312  } // for (v < n_vertex)
5313 
5314  // If we still haven't found it's probably the last point along
5315  if (!found)
5316  {
5317 #ifdef PARANOID
5318  double diff = std::fabs(
5319  zeta[0] - polygonal_vertex_arclength[n_vertex - 1].first);
5320  if (diff >
5322  {
5323  std::ostringstream error_stream;
5324  error_stream
5325  << "Wasn't able to locate the polygonal arclength exactly\n"
5326  << "during re-setup of boundary coordinates and have\n"
5327  << "assumed that we're dealing with the final point along\n"
5328  << "the curvilinear segment and encountered some roundoff\n"
5329  << "However,the difference in the polygonal zeta "
5330  "coordinates\n"
5331  << "between zeta[0] " << zeta[0]
5332  << " and the originallly stored value "
5333  << polygonal_vertex_arclength[n_vertex - 1].first << "\n"
5334  << "is " << diff
5335  << " which exceeds the threshold specified\n"
5336  << "in the publically modifiable variable\n"
5337  << "ToleranceForVertexMismatchInPolygons::Tolerable_error\n"
5338  << "whose current value is: "
5340  << "\nPlease check your mesh carefully and increase the\n"
5341  << "threshold if you're sure this is appropriate\n";
5342  throw OomphLibError(error_stream.str(),
5343  OOMPH_CURRENT_FUNCTION,
5344  OOMPH_EXCEPTION_LOCATION);
5345  }
5346 #endif
5347  zeta[0] = polygonal_vertex_arclength[n_vertex - 1].second;
5348  }
5349  }
5350 
5351  // Assign updated coordinate
5352  nod_pt->set_coordinates_on_boundary(b, zeta);
5353  }
5354  } // if (geom_object_pt != 0)
5355  else
5356  {
5357  // Establish the initial zeta value for this segment
5358  double z_initial = initial_segment_arclength[is];
5359 
5360  // Only use end points of the whole boundary and pick the
5361  // bottom left node
5362 
5363  // Set the bottom left coordinate as the first coordinates
5364  Vector<double> bottom_left_coordinate(2);
5365  bottom_left_coordinate = first_coordinate;
5366 
5367  // ... do the same with the zeta value
5368  Vector<double> bottom_left_zeta_coordinate(1);
5369  bottom_left_zeta_coordinate = first_node_zeta_coordinate;
5370 
5371  // Is the last y-coordinate smaller than y-coordinate of the
5372  // current bottom left coordinate
5373  if (last_coordinate[1] < bottom_left_coordinate[1])
5374  {
5375  // Re-set the bottom-left coordinate as the last coordinate
5376  bottom_left_coordinate = last_coordinate;
5377 
5378  // ... do the same with the zeta value
5379  bottom_left_zeta_coordinate = last_node_zeta_coordinate;
5380  }
5381  // The y-coordinates are the same
5382  else if (last_coordinate[1] == bottom_left_coordinate[1])
5383  {
5384  // Then check for the x-coordinate, which is the most-left
5385  if (last_coordinate[0] < bottom_left_coordinate[0])
5386  {
5387  // Re-set the bottom-left coordinate as the last
5388  // coordinate
5389  bottom_left_coordinate = last_coordinate;
5390 
5391  // ... do the same with the zeta value
5392  bottom_left_zeta_coordinate = last_node_zeta_coordinate;
5393  }
5394  } // else (The y-coordinates are the same)
5395 
5396  Vector<double> zeta(1, 0.0);
5397 
5398  // Now adjust boundary coordinate so that the bottom left
5399  // node has a boundary coordinate of 0.0 and that zeta
5400  // increases away from that point
5401  zeta = bottom_left_zeta_coordinate;
5402  const double zeta_ref = zeta[0];
5403 
5404  // Also get the maximum zeta value
5405  double zeta_max = 0.0;
5406  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
5407  it != all_nodes_pt.end();
5408  it++)
5409  {
5410  Node* nod_pt = (*it);
5411  nod_pt->get_coordinates_on_boundary(b, zeta);
5412 
5413 #ifdef OOMPH_HAS_MPI
5414 
5415  // If the mesh is distributed and the initial and final
5416  // zeta values for the segment have been assigned then
5417  // check if the segment is inverted, we need to consider
5418  // that to correctly assgin the zeta values for the segment
5419  if (this->is_mesh_distributed() &&
5421  {
5422  // Is the segment inverted, if that is the case then
5423  // invert the zeta values
5424  if (boundary_segment_inverted(b)[is])
5425  {
5426  zeta[0] = segment_arclength[is] - zeta[0];
5427  } // if (boundary_segment_inverted(b)[is])
5428  }
5429 
5430 #endif // #ifdef OOMPH_HAS_MPI
5431 
5432  // Set the zeta value
5433  zeta[0] += z_initial;
5434  // Adjust the value based on the bottom-left condition
5435  zeta[0] -= zeta_ref;
5436  // If direction is reversed, then take absolute value
5437  if (zeta[0] < 0.0)
5438  {
5439  zeta[0] = std::fabs(zeta[0]);
5440  }
5441  // Get the max zeta value (we will use it to scale the
5442  // values to [0,1])
5443  if (zeta[0] > zeta_max)
5444  {
5445  zeta_max = zeta[0];
5446  }
5447  nod_pt->set_coordinates_on_boundary(b, zeta);
5448  } // Loop through the nodes in the segment (boundary)
5449 
5450 #ifdef OOMPH_HAS_MPI
5451  // After assigning boundary coordinates, BUT before scaling,
5452  // copy the initial and final segment arclengths so that we
5453  // know if the values go in increasing or decreasing
5454  // order. This will be used to identify the correct direction
5455  // of the segments in the new meshes created in the
5456  // adaptation method.
5457  if (this->is_mesh_distributed() &&
5459  {
5460  // Get the first face element
5461  FiniteElement* first_seg_ele_pt = segment_sorted_ele_pt[is].front();
5462 
5463 #ifdef PARANOID
5464  // Check if the face element is nonhalo, it shouldn't, but
5465  // better check
5466  if (first_seg_ele_pt->is_halo())
5467  {
5468  std::ostringstream error_message;
5469  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
5470  "setup_boundary_coordinates()";
5471  error_message << "The first face element in the (" << is
5472  << ")-th segment"
5473  << " is halo\n";
5474  throw OomphLibError(
5475  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5476  } // if (first_seg_ele_pt->is_halo())
5477 #endif
5478 
5479  // Number of nodes
5480  const unsigned nnod = first_seg_ele_pt->nnode();
5481 
5482  // Get the first node of the current segment
5483  Node* first_seg_node_pt = first_seg_ele_pt->node_pt(0);
5484  if (is_inverted[first_seg_ele_pt])
5485  {
5486  first_seg_node_pt = first_seg_ele_pt->node_pt(nnod - 1);
5487  }
5488 
5489  // Get the last face element
5490  FiniteElement* last_seg_ele_pt = segment_sorted_ele_pt[is].back();
5491 
5492 #ifdef PARANOID
5493  // Check if the face element is nonhalo, it shouldn't, but
5494  // better check
5495  if (last_seg_ele_pt->is_halo())
5496  {
5497  std::ostringstream error_message;
5498  std::string output_string = "UnstructuredTwoDMeshGeometryBase::"
5499  "setup_boundary_coordinates()";
5500  error_message << "The last face element in the (" << is
5501  << ")-th segment"
5502  << " is halo\n";
5503  throw OomphLibError(
5504  error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
5505  } // if (last_seg_ele_pt->is_halo())
5506 #endif
5507 
5508  // Get the last node of the current segment
5509  Node* last_seg_node_pt = last_seg_ele_pt->node_pt(nnod - 1);
5510  if (is_inverted[last_seg_ele_pt])
5511  {
5512  last_seg_node_pt = last_seg_ele_pt->node_pt(0);
5513  }
5514 
5515  // Now get the first and last node boundary coordinate values
5516  Vector<double> first_seg_arclen(1);
5517  Vector<double> last_seg_arclen(1);
5518 
5519  first_seg_node_pt->get_coordinates_on_boundary(b, first_seg_arclen);
5520  last_seg_node_pt->get_coordinates_on_boundary(b, last_seg_arclen);
5521 
5522  // Update the initial and final segments arclength
5523  boundary_segment_initial_arclength(b)[is] = first_seg_arclen[0];
5524  boundary_segment_final_arclength(b)[is] = last_seg_arclen[0];
5525 
5526  // Update the initial and final coordinates
5527  Vector<double> updated_segment_initial_coord(2);
5528  Vector<double> updated_segment_final_coord(2);
5529  for (unsigned k = 0; k < 2; k++)
5530  {
5531  updated_segment_initial_coord[k] = first_seg_node_pt->x(k);
5532  updated_segment_final_coord[k] = last_seg_node_pt->x(k);
5533  }
5534 
5535  // Set the updated initial coordinate
5537  updated_segment_initial_coord;
5538 
5539  // Set the updated final coordinate
5541  updated_segment_final_coord;
5542 
5543  } // if (this->is_mesh_distributed() &&
5544  // Assigned_segments_initial_zeta_values[b])
5545 #endif // #ifdef OOMPH_HAS_MPI
5546 
5547  // The max. value will be incorrect if we are working with
5548  // distributed meshes where the current boundary has been
5549  // split by the distribution process. Here correct the
5550  // maximum value
5551  if (zeta_max < boundary_arclength)
5552  {
5553  zeta_max = boundary_arclength;
5554  }
5555 
5556  // Scale all surface coordinates so that all the points be on
5557  // the range [0, 1]
5558  for (std::set<Node*>::iterator it = all_nodes_pt.begin();
5559  it != all_nodes_pt.end();
5560  it++)
5561  {
5562  // Get the node
5563  Node* nod_pt = (*it);
5564 
5565  // Get the boundary coordinate
5566  nod_pt->get_coordinates_on_boundary(b, zeta);
5567 
5568  // Scale the value of the current node
5569  zeta[0] /= zeta_max;
5570 
5571  // Set the new scaled value
5572  nod_pt->set_coordinates_on_boundary(b, zeta);
5573  } // Loop over the nodes
5574 
5575  } // else if (geom_object_pt != 0)
5576 
5577  } // for (is < nsegments)
5578 
5579  // =================================================================
5580  // END: Scale the boundary coordinates based on whether the
5581  // boundary has an associated GeomObj or not
5582  // =================================================================
5583 
5584  // Cleanup
5585  for (unsigned e = 0; e < nel; e++)
5586  {
5587  delete face_el_pt[e];
5588  face_el_pt[e] = 0;
5589  }
5590 
5591  } // if (nel > 0), from the beginning of the method
5592 
5593  // Indicate that boundary coordinate has been set up
5594  Boundary_coordinate_exists[b] = true;
5595  }
5596 
5597 
5598  /// ///////////////////////////////////////////////////////////////////////
5599  /// ///////////////////////////////////////////////////////////////////////
5600  /// ///////////////////////////////////////////////////////////////////////
5601 
5602 
5603  //=================================================
5604  /// Helper namespace for BCInfo object used
5605  /// in the identification of boundary elements.
5606  //=================================================
5607  namespace TriangleBoundaryHelper
5608  {
5609  /// Structure for Boundary Informations
5610  struct BCInfo
5611  {
5612  /// Face ID
5613  unsigned Face_id;
5614 
5615  /// Boundary ID
5616  unsigned Boundary;
5617 
5618  /// Pointer to bulk finite element
5620  };
5621 
5622  } // namespace TriangleBoundaryHelper
5623 
5624 } // namespace oomph
5625 
5626 #endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: elements.h:5014
A general Finite Element class.
Definition: elements.h:1313
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
bool is_halo() const
Is this element a halo?
Definition: elements.h:1163
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
A general mesh class.
Definition: mesh.h:67
bool is_mesh_distributed() const
Boolean to indicate if Mesh has been distributed.
Definition: mesh.h:1588
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition: mesh.h:190
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:896
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:878
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:840
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
virtual void set_coordinates_on_boundary(const unsigned &b, const unsigned &k, const Vector< double > &boundary_zeta)
Set the vector of the k-th generalised boundary coordinates on mesh boundary b. Broken virtual interf...
Definition: nodes.cc:2394
virtual void get_coordinates_on_boundary(const unsigned &b, const unsigned &k, Vector< double > &boundary_zeta)
Return the vector of the k-th generalised boundary coordinates on mesh boundary b....
Definition: nodes.cc:2379
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
void output(std::ostream &outfile, const unsigned &n_sample=50)
Output each sub-boundary at n_sample (default: 50) points.
Vector< double > internal_point() const
Coordinates of the internal point.
TriangleMeshClosedCurve(const Vector< TriangleMeshCurveSection * > &curve_section_pt, const Vector< double > &internal_point_pt=Vector< double >(0), const bool &is_internal_point_fixed=false)
Constructor prototype.
unsigned nvertices() const
Number of vertices.
Vector< double > Internal_point_pt
Vector of vertex coordinates.
bool is_internal_point_fixed() const
Test whether the internal point is fixed.
bool Is_internal_point_fixed
Indicate whether the internal point should be updated automatically.
unsigned nsegments() const
Total number of segments.
void unfix_internal_point()
Unfix the internal point (i.e. allow our automatic machinery to update it)
Vector< double > & internal_point()
Coordinates of the internal point.
void fix_internal_point()
Fix the internal point (i.e. do not allow our automatic machinery to update it)
Base class for defining a triangle mesh boundary, this class has the methods that allow to connect th...
bool Final_vertex_connected_suspended
Indicates if the connection is suspended because the boundary to connect is no longer part of the dom...
bool Initial_vertex_connected_to_curviline
States if the initial vertex is connected to a curviline.
virtual void initial_vertex_coordinate(Vector< double > &vertex)=0
Get first vertex coordinates.
bool is_initial_vertex_connected_to_curviline() const
Test whether the initial vertex is connected to a curviline.
virtual unsigned boundary_chunk() const =0
Boundary chunk (Used when a boundary is represented by more than one polyline.
virtual void output(std::ostream &outfile, const unsigned &n_sample=50)=0
Output the curve_section.
bool Final_vertex_connected
Used for stating if the final end is connected to another boundary.
double unrefinement_tolerance()
Get tolerance for unrefinement of curve section to create a better representation of curvilinear boun...
void unset_final_vertex_connected_to_curviline()
Sets the final vertex as non connected to a curviline.
unsigned Final_vertex_connected_n_vertex
Stores the vertex number used for connection with the final end.
double final_s_connection_value() const
Gets the s value to which the final end is connected.
double Unrefinement_tolerance
Tolerance for unrefinement of curve sections (neg if refinement is disabled)
unsigned & final_vertex_connected_n_vertex()
Gets the vertex number to which the final end is connected.
bool Final_vertex_connected_to_curviline
States if the final vertex is connected to a curviline.
double Initial_s_connection_value
Stores the s value used for connecting the initial end with a curviline.
void set_unrefinement_tolerance(const double &tolerance)
Set tolerance for unrefinement of curve sections to avoid unnecessarily large numbers of elements on ...
virtual unsigned boundary_id() const =0
Boundary id.
unsigned & final_vertex_connected_bnd_id()
Sets the id to which the final end is connected.
void connect_initial_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Connects the initial vertex of the curve section to a desired target polyline by specifying the verte...
double & tolerance_for_s_connection()
Sets the tolerance value for connections among curvilines.
void resume_initial_vertex_connected()
Resumes the initial vertex connection, it may be that after load balancing the boundary to which the ...
double Final_s_connection_value
Stores the s value used for connecting the final end with a curviline.
double tolerance_for_s_connection() const
Gets the tolerance value for connections among curvilines.
unsigned Final_vertex_connected_bnd_id
Stores the id to which the initial end is connected.
void suspend_final_vertex_connected()
Set the final vertex connection as suspended, it will be resumed when the method to resume the connec...
unsigned Initial_vertex_connected_n_chunk
Stores the chunk number of the boundary to which is connected th initial end.
unsigned initial_vertex_connected_n_vertex() const
Gets the vertex number to which the initial end is connected.
void resume_final_vertex_connected()
Resumes the final vertex connection, it may be that after load balancing the boundary to which the co...
double Refinement_tolerance
Tolerance for refinement of curve sections (neg if refinement is disabled)
void set_refinement_tolerance(const double &tolerance)
Set tolerance for refinement of curve sections to create a better representation of curvilinear bound...
virtual unsigned nvertex() const =0
Number of vertices.
virtual unsigned nsegment() const =0
Number of segments that this part of the boundary is to be represented by. This corresponds to the nu...
double Maximum_length
Maximum allowed distance between two vertices on the polyline representation of the curve section.
void connect_final_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Connects the final vertex of the curve section to a desired target polyline by specifying the vertex ...
unsigned final_vertex_connected_n_chunk() const
Gets the boundary chunk to which the final end is connected.
void set_maximum_length(const double &maximum_length)
Allows to specify the maximum distance between two vertices that define the associated polyline of th...
void unset_initial_vertex_connected()
Sets the initial vertex as non connected.
bool Initial_vertex_connected
Used for stating if the initial end is connected to another boundary.
unsigned final_vertex_connected_n_vertex() const
Sets the vertex number to which the final end is connected.
double refinement_tolerance()
Get tolerance for refinement of curve sections to create a better representation of curvilinear bound...
unsigned & final_vertex_connected_n_chunk()
Sets the boundary chunk to which the final end is connected.
unsigned Initial_vertex_connected_n_vertex
Stores the vertex number used for connection with the initial end.
bool Initial_vertex_connected_suspended
Indicates if the connection is suspended because the boundary to connect is no longer part of the dom...
TriangleMeshCurveSection()
Empty constructor. Initialises the curve section as non connected.
double & initial_s_connection_value()
Sets the s value to which the initial end is connected.
void disable_refinement_tolerance()
Disable refinement of curve section.
bool is_final_vertex_connected_to_curviline() const
Test whether the final vertex is connected to a curviline.
bool is_final_vertex_connected() const
Test whether final vertex is connected or not.
double Tolerance_for_s_connection
Tolerance used for connecting the ends to a curviline.
void disable_unrefinement_tolerance()
Disable unrefinement of curve sections.
unsigned & initial_vertex_connected_n_chunk()
Sets the boundary chunk to which the initial end is connected.
unsigned & initial_vertex_connected_bnd_id()
Sets the id to which the initial end is connected.
unsigned initial_vertex_connected_bnd_id() const
Gets the id to which the initial end is connected.
virtual void final_vertex_coordinate(Vector< double > &vertex)=0
Get last vertex coordinates.
bool is_initial_vertex_connected() const
Test whether initial vertex is connected or not.
double & final_s_connection_value()
Sets the s value to which the final end is connected.
void enable_refinement_tolerance(const double &tolerance=0.08)
Enable refinement of curve section to create a better representation of curvilinear boundaries (e....
void set_initial_vertex_connected()
Sets the initial vertex as connected.
unsigned Final_vertex_connected_n_chunk
Stores the chunk number of the boundary to which is connected th initial end.
void set_final_vertex_connected()
Sets the final vertex as connected.
void enable_unrefinement_tolerance(const double &tolerance=0.04)
Enable unrefinement of curve sections to avoid unnecessarily large numbers of elements on of curvilin...
void unset_final_vertex_connected()
Sets the final vertex as non connected.
void connect_initial_vertex_to_curviline(TriangleMeshCurviLine *curviline_pt, const double &s_value, const double &tolerance_for_connection=1.0e-14)
Connects the initial vertex of the curve section to a desired target curviline by specifying the s va...
void set_initial_vertex_connected_to_curviline()
Sets the initial vertex as connected to a curviline.
unsigned & initial_vertex_connected_n_vertex()
Sets the vertex number to which the initial end is connected.
double initial_s_connection_value() const
Gets the s value to which the initial end is connected.
void suspend_initial_vertex_connected()
Set the initial vertex connection as suspended, it will be resumed when the method to resume the conn...
unsigned initial_vertex_connected_n_chunk() const
Gets the boundary chunk to which the initial end is connected.
double maximum_length()
Gets access to the maximum length variable.
void connect_final_vertex_to_curviline(TriangleMeshCurviLine *curviline_pt, const double &s_value, const double &tolerance_for_connection=1.0e-14)
Connects the final vertex of the curve section to a desired target curviline by specifying the s valu...
void set_final_vertex_connected_to_curviline()
Sets the final vertex as connected to a curviline.
unsigned final_vertex_connected_bnd_id() const
Gets the id to which the final end is connected.
unsigned Initial_vertex_connected_bnd_id
Stores the id to which the initial end is connected.
void disable_use_maximum_length()
Disables the use of the maximum length criteria on the unrefinement or refinement steps.
void unset_initial_vertex_connected_to_curviline()
Sets the initial vertex as non connected to a curviline.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
double Polyline_unrefinement_tolerance
Tolerance for unrefinement of polylines (neg if refinement is disabled)
void set_polyline_refinement_tolerance(const double &tolerance)
Set tolerance for refinement of polylines to create a better representation of curvilinear boundaries...
void set_polyline_unrefinement_tolerance(const double &tolerance)
Set tolerance for unrefinement of polylines to avoid unnecessarily large numbers of elements on of cu...
void disable_polyline_unrefinement()
Disable unrefinement of polylines.
Vector< TriangleMeshCurveSection * > Curve_section_pt
Vector of curve sections.
TriangleMeshCurve(const Vector< TriangleMeshCurveSection * > &curve_section_pt)
Empty constructor.
virtual void output(std::ostream &outfile, const unsigned &n_sample=50)=0
Output each sub-boundary at n_sample (default: 50) points.
void disable_polyline_refinement()
Disable refinement of polylines.
virtual TriangleMeshCurveSection *& curve_section_pt(const unsigned &i)
Pointer to i-th constituent curve section.
virtual TriangleMeshCurveSection * curve_section_pt(const unsigned &i) const
Pointer to i-th constituent curve section.
double polyline_unrefinement_tolerance()
Get tolerance for unrefinement of polylines to create a better representation of curvilinear boundari...
void enable_polyline_refinement(const double &tolerance=0.08)
Enable refinement of polylines to create a better representation of curvilinear boundaries (e....
virtual unsigned nsegments() const =0
Total number of segments.
double Polyline_refinement_tolerance
Tolerance for refinement of polylines (neg if refinement is disabled)
double polyline_refinement_tolerance()
Get tolerance for refinement of polylines to create a better representation of curvilinear boundaries...
virtual unsigned nvertices() const =0
Number of vertices.
void enable_polyline_unrefinement(const double &tolerance=0.04)
Enable unrefinement of polylines to avoid unnecessarily large numbers of elements on of curvilinear b...
unsigned max_boundary_id()
Return max boundary id of associated curves.
virtual unsigned ncurve_section() const
Number of constituent curves.
Class definining a curvilinear triangle mesh boundary in terms of a GeomObject. Curvlinear equivalent...
GeomObject * geom_object_pt()
Pointer to GeomObject that represents this part of the boundary.
bool are_there_connection_points()
Does the vector for storing connections has elements?
void initial_vertex_coordinate(Vector< double > &vertex)
Get first vertex coordinates.
unsigned & nsegment()
Number of (initially straight-line) segments that this part of the boundary is to be represented by....
Vector< double > Connection_points_pt
Stores the information for connections received on the curviline. Used when converting to polyline.
void add_connection_point(const double &z_value, const double &tol=1.0e-12)
Add the connection point (z_value) to the connection points that receive the curviline.
bool Space_vertices_evenly_in_arclength
Boolean to indicate if vertices in polygonal representation of the Curviline are spaced (approximatel...
unsigned boundary_chunk() const
Boundary chunk (Used when a boundary is represented by more than one polyline.
void output(std::ostream &outfile, const unsigned &n_sample=50)
Output curvilinear boundary at n_sample (default: 50) points.
unsigned nvertex() const
Number of vertices.
double Zeta_end
End coordinate in terms of the GeomObject's intrinsic coordinate.
void final_vertex_coordinate(Vector< double > &vertex)
Get last vertex coordinates.
double zeta_start()
Start coordinate in terms of the GeomObject's intrinsic coordinate.
GeomObject * Geom_object_pt
Pointer to GeomObject that represents this part of the boundary.
TriangleMeshCurviLine(GeomObject *geom_object_pt, const double &zeta_start, const double &zeta_end, const unsigned &nsegment, const unsigned &boundary_id, const bool &space_vertices_evenly_in_arclength=true, const unsigned &boundary_chunk=0)
Constructor: Specify GeomObject, the start and end coordinates of the relevant boundary in terms of t...
unsigned Boundary_chunk
Boundary chunk (Used when a boundary is represented by more than one polyline.
double Zeta_start
Start coordinate in terms of the GeomObject's intrinsic coordinate.
Vector< double > * connection_points_pt()
Returns the connection points vector.
double zeta_end()
End coordinate in terms of the GeomObject's intrinsic coordinate.
unsigned Nsegment
Number of (initially straight-line) segments that this part of the boundary is to be represented by.
bool space_vertices_evenly_in_arclength() const
Boolean to indicate if vertices in polygonal representation of the Curvline are spaced (approximately...
unsigned nsegment() const
Number of (initially straight-line) segments that this part of the boundary is to be represented by.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
TriangleMeshOpenCurve(const Vector< TriangleMeshCurveSection * > &curve_section_pt)
Constructor.
TriangleMeshPolyLine * polyline_pt(const unsigned &i)
Pointer to i-th constituent polyline.
unsigned nsegments() const
Total number of segments.
unsigned nvertices() const
Number of vertices.
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
void output(std::ostream &outfile, const unsigned &n_sample=50)
Output each sub-boundary at n_sample (default: 50) points.
Class defining a polyline for use in Triangle Mesh generation.
unsigned boundary_chunk() const
Boundary chunk (Used when a boundary is represented by more than one polyline.
unsigned Boundary_chunk
Boundary chunk (Used when a boundary is represented by more than one polyline.
unsigned nsegment() const
Number of segments.
TriangleMeshPolyLine(const Vector< Vector< double >> &vertex_coordinate, const unsigned &boundary_id, const unsigned &boundary_chunk=0)
Constructor: Takes vectors of vertex coordinates in order Also allows the optional specification of a...
void final_vertex_coordinate(Vector< double > &vertex)
Get last vertex coordinates.
unsigned nvertex() const
Number of vertices.
void output(std::ostream &outfile, const unsigned &n_sample=50)
Output the polyline – n_sample is ignored.
Vector< double > & vertex_coordinate(const unsigned &i)
Coordinate vector of i-th vertex.
Vector< Vector< double > > Vertex_coordinate
Vector of Vector of vertex coordinates.
void reverse()
Reverse the polyline, this includes the connection information and the vertices order.
void initial_vertex_coordinate(Vector< double > &vertex)
Get first vertex coordinates.
Vector< double > vertex_coordinate(const unsigned &i) const
Coordinate vector of i-th vertex (const version)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
unsigned npolyline() const
Number of constituent polylines.
TriangleMeshPolyLine * polyline_pt(const unsigned &i)
Pointer to i-th constituent polyline.
bool is_fixed() const
Test whether the polygon is fixed or not.
bool can_update_reference_configuration() const
Test whether curve can update reference.
void set_fixed()
Set the polygon to be fixed.
Vector< unsigned > polygon_boundary_id()
Return vector of boundary ids of associated polylines.
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
bool is_redistribution_of_segments_between_polylines_enabled()
Is re-distribution of polyline segments in the curve between different boundaries during adaptation e...
TriangleMeshPolygon(const Vector< TriangleMeshCurveSection * > &boundary_polyline_pt, const Vector< double > &internal_point_pt=Vector< double >(0), const bool &is_internal_point_fixed=false)
Constructor: Specify vector of pointers to TriangleMeshCurveSection that define the boundary of the s...
unsigned ncurve_section() const
Number of constituent curves.
bool Can_update_configuration
Boolean flag to indicate whether the polygon can update its own reference configuration after it has ...
bool Polygon_fixed
Boolean flag to indicate whether the polygon can move (default false because if it doesn't move this ...
void set_unfixed()
Set the polygon to be allowed to move (default)
virtual ~TriangleMeshPolygon()
Empty virtual destructor.
virtual void reset_reference_configuration()
Virtual function that should be overloaded to update the polygons reference configuration.
void disable_redistribution_of_segments_between_polylines()
Disable re-distribution of polyline segments in the curve between different boundaries during adaptat...
bool Enable_redistribution_of_segments_between_polylines
Is re-distribution of polyline segments between different boundaries during adaptation enabled?...
void enable_redistribution_of_segments_between_polylines()
Enable re-distribution of polyline segments in the curve between different boundaries during adaptati...
Contains functions which define the geometry of the mesh, i.e. regions, boundaries,...
std::map< unsigned, GeomObject * > Boundary_geom_object_pt
Storage for the geometric objects associated with any boundaries.
std::map< unsigned, Vector< double > > & boundary_segment_final_zeta()
Return direct access to the final zeta for the segments that are part of a boundary.
void check_contiguousness_on_polylines_helper(Vector< TriangleMeshPolyLine * > &polylines_pt, unsigned &index)
Sort the polylines coming from joining them. Check whether it is necessary to reverse them or not....
Vector< TriangleMeshOpenCurve * > Internal_open_curve_pt
Vector of open polylines that define internal curves.
const Vector< unsigned > boundary_segment_inverted(const unsigned &b) const
Return the info. to know if it is necessary to reverse the segment based on a previous mesh.
std::map< unsigned, Vector< double > > & boundary_initial_zeta_coordinate()
Return direct access to the initial zeta coordinate of a boundary.
double region_attribute(const unsigned &i)
Return the attribute associated with region i.
Vector< double > & boundary_segment_final_arclength(const unsigned &b)
Return the final arclength for the segments that are part of a boundary.
static bool Suppress_warning_about_regions_and_boundaries
Public static flag to suppress warning; defaults to false.
std::map< unsigned, Vector< Vector< double > > > Boundary_segment_initial_coordinate
Stores the initial coordinates for the segments that appear when a boundary is splited among processo...
void disable_automatic_creation_of_vertices_on_boundaries()
Disables the creation of points (by Triangle) on the outer and internal boundaries.
std::map< unsigned, Vector< double > > & boundary_final_coordinate()
Return direct access to the final coordinates of a boundary.
bool is_automatic_creation_of_vertices_on_boundaries_allowed()
Returns the status of the variable Allow_automatic_creation_of_vertices_on_boundaries.
const bool get_connected_vertex_number_on_destination_polyline(TriangleMeshPolyLine *dst_polyline_pt, Vector< double > &vertex_coordinates, unsigned &vertex_number)
Gets the vertex number on the destination polyline (used to create the connections among shared bound...
Vector< double > & boundary_final_zeta_coordinate(const unsigned &b)
Return the final zeta coordinate for the boundary.
Vector< TriangleMeshPolygon * > Outer_boundary_pt
Polygon that defines outer boundaries.
Vector< Vector< double > > & boundary_segment_final_coordinate(const unsigned &b)
Return the final coordinates for the segments that are part of a boundary.
void initialise_base_vertex(TriangleMeshPolyLine *polyline_pt, std::map< unsigned, std::map< unsigned, Vector< base_vertex_info >>> &base_vertices)
Initialise the base vertex structure, set every vertex to no visited and not being a base vertex.
std::map< unsigned, Vector< double > > Boundary_segment_initial_zeta
Stores the initial zeta coordinate for the segments that appear when a boundary is splited among proc...
void add_connection_matrix_info_helper(TriangleMeshPolyLine *polyline_pt, std::map< unsigned, std::map< unsigned, Vector< vertex_connection_info >>> &connection_matrix, TriangleMeshPolyLine *next_polyline_pt=0)
Helps to add information to the connection matrix of the given polyline.
void set_nboundary_segment_node(const unsigned &b, const unsigned &s)
Set the number of segments associated with a boundary.
void check_contiguousness_on_polylines_helper(Vector< TriangleMeshPolyLine * > &polylines_pt, unsigned &index_halo_start, unsigned &index_halo_end)
Sort the polylines coming from joining them. Check whether it is necessary to reverse them or not....
Vector< TriangleMeshPolygon * > Internal_polygon_pt
Vector of polygons that define internal polygons.
std::map< unsigned, Vector< unsigned > > & boundary_segment_inverted()
Return the info. to know if it is necessary to reverse the segment based on a previous mesh.
std::map< unsigned, Vector< Vector< double > > > Boundary_segment_final_coordinate
Stores the final coordinates for the segments that appear when a boundary is splited among processors...
std::map< unsigned, Vector< double > > & boundary_segment_initial_zeta()
Return direct access to the initial zeta for the segments that are part of a boundary.
Vector< double > & boundary_segment_initial_arclength(const unsigned &b)
Return the initial arclength for the segments that are part of a boundary.
std::map< unsigned, bool > Assigned_segments_initial_zeta_values
Flag used at the setup_boundary_coordinate function to know if initial zeta values for segments have ...
unsigned nboundary_segment(const unsigned &b)
Return the number of segments associated with a boundary.
void snap_nodes_onto_geometric_objects()
Snap the boundary nodes onto any curvilinear boundaries defined by geometric objects.
std::map< unsigned, Vector< double > > Boundary_initial_zeta_coordinate
Stores the initial zeta coordinate for the boundary.
std::map< unsigned, std::set< Node * > > Nodes_on_boundary_pt
Stores a pointer to a set with all the nodes related with a boundary.
unsigned long nboundary_segment_node(const unsigned &b)
Return the number of segments associated with a boundary.
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...
FiniteElement * boundary_element_in_region_pt(const unsigned &b, const unsigned &r, const unsigned &e) const
Return pointer to the e-th element adjacent to boundary b in region r.
std::map< unsigned, Vector< Vector< Node * > > > Boundary_segment_node_pt
Used to store the nodes associated to a boundary and to an specific segment (this only applies in dis...
std::map< unsigned, Vector< unsigned > > Boundary_segment_inverted
Stores the info. to know if it is necessary to reverse the segment based on a previous mesh.
std::map< unsigned, Vector< double > > & boundary_segment_final_arclength()
Return direct access to the final arclength for the segments that are part of a boundary.
std::map< unsigned, Vector< double > > Regions_coordinates
Storage for extra coordinates for regions. The key on the map is the region id.
void enable_automatic_creation_of_vertices_on_boundaries()
Enables the creation of points (by Triangle) on the outer and internal boundaries.
Vector< Vector< double > > Extra_holes_coordinates
Storage for extra coordinates for holes.
std::map< unsigned, GeomObject * > & boundary_geom_object_pt()
Return direct access to the geometric object storage.
bool Allow_automatic_creation_of_vertices_on_boundaries
Flag to indicate whether the automatic creation of vertices along the boundaries by Triangle is allow...
UnstructuredTwoDMeshGeometryBase(const UnstructuredTwoDMeshGeometryBase &dummy)=delete
Broken copy constructor.
int face_index_at_boundary_in_region(const unsigned &b, const unsigned &r, const unsigned &e) const
Return face index of the e-th element adjacent to boundary b in region r.
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< Vector< double > > > & boundary_segment_initial_coordinate()
Return direct access to the initial coordinates for the segments that are part of a boundary.
std::map< unsigned, Vector< double > > Boundary_coordinate_limits
Storage for the limits of the boundary coordinates defined by the use of geometric objects....
std::map< unsigned, Vector< double > > Boundary_segment_final_arclength
Stores the final arclength for the segments that appear when a boundary is splited among processors.
std::map< unsigned, Vector< double > > & boundary_segment_initial_arclength()
Return direct access to the initial arclength for the segments that are part of a boundary.
unsigned get_associated_vertex_to_svalue(double &target_s_value, unsigned &bnd_id)
Get the associated vertex to the given s value by looking to the list of s values created when changi...
Vector< double > Region_attribute
Vector of attributes associated with the elements in each region.
Vector< std::map< unsigned, Vector< int > > > Face_index_region_at_boundary
Storage for the face index adjacent to a boundary in a particular region.
void copy_connection_information_to_sub_polylines(TriangleMeshCurveSection *input_curve_pt, TriangleMeshCurveSection *output_curve_pt)
Helper function to copy the connection information from the input curve(polyline or curviline) to the...
unsigned nboundary_element_in_region(const unsigned &b, const unsigned &r) const
Return the number of elements adjacent to boundary b in region r.
void add_boundary_segment_node(const unsigned &b, const unsigned &s, Node *const &node_pt)
Add the node node_pt to the b-th boundary and the s-th segment of the mesh.
Vector< double > & boundary_initial_coordinate(const unsigned &b)
Return the initial coordinates for the boundary.
unsigned nregion()
Return the number of regions specified by attributes.
TriangleMeshPolyLine * boundary_polyline_pt(const unsigned &b)
Return pointer to the current polyline that describes the b-th mesh boundary.
Vector< double > & boundary_segment_final_zeta(const unsigned &b)
Return the final zeta for the segments that are part of a boundary.
void add_base_vertex_info_helper(TriangleMeshPolyLine *polyline_pt, std::map< unsigned, std::map< unsigned, Vector< base_vertex_info >>> &base_vertices, std::map< unsigned, std::map< unsigned, Vector< vertex_connection_info >>> &connection_matrix, std::map< unsigned, std::map< unsigned, unsigned >> &boundary_chunk_nvertices)
Helps to identify the base vertex of the given polyline.
std::map< unsigned, std::set< Node * > > & nodes_on_boundary_pt()
Gets a pointer to a set with all the nodes related with a boundary.
void setup_boundary_coordinates(const unsigned &b)
Setup boundary coordinate on boundary b. Boundary coordinate increases continously along polygonal bo...
void copy_connection_information(TriangleMeshCurveSection *input_curve_pt, TriangleMeshCurveSection *output_curve_pt)
Helper function to copy the connection information from the input curve(polyline or curviline) to the...
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...
Vector< unsigned > & boundary_segment_inverted(const unsigned &b)
Return the info. to know if it is necessary to reverse the segment based on a previous mesh.
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...
Vector< Vector< double > > & boundary_segment_initial_coordinate(const unsigned &b)
Return the initial coordinates for the segments that are part of a boundary.
Vector< double > & boundary_initial_zeta_coordinate(const unsigned &b)
Return the initial zeta coordinate for the boundary.
Vector< double > & boundary_final_coordinate(const unsigned &b)
Return the final coordinates for the boundary.
unsigned nregion_attribute()
Return the number of attributes used in the mesh.
unsigned get_associated_vertex_to_svalue(double &target_s_value, unsigned &bnd_id, double &s_tolerance)
Get the associated vertex to the given s value by looking to the list of s values created when changi...
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....
FiniteElement * region_element_pt(const unsigned &i, const unsigned &e)
Return the e-th element in the i-th region.
GeomObject * boundary_geom_object_pt(const unsigned &b)
Return the geometric object associated with the b-th boundary or null if the boundary has associated ...
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...
Vector< double > & boundary_coordinate_limits(const unsigned &b)
Return access to the coordinate limits associated with the geometric object associated with boundary ...
std::map< unsigned, Vector< std::pair< double, double > > > Polygonal_vertex_arclength_info
Storage for pairs of doubles representing: .first: the arclength along the polygonal representation o...
std::map< unsigned, Vector< double > > Boundary_final_coordinate
Stores the final coordinates for the boundary.
std::map< unsigned, Vector< double > > Boundary_initial_coordinate
Stores the initial coordinates for the boundary.
std::set< TriangleMeshPolygon * > Free_polygon_pt
A set that contains the polygons created by this object therefore it is necessary to free their assoc...
std::map< unsigned, Vector< FiniteElement * > > Region_element_pt
Vector of elements in each region differentiated by attribute (the key of the map is the attribute)
std::map< unsigned, Vector< double > > Boundary_final_zeta_coordinate
Stores the final zeta coordinate for the boundary.
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....
std::map< unsigned, Vector< double > > Boundary_segment_initial_arclength
Stores the initial arclength for the segments that appear when a boundary is splited among processors...
std::map< unsigned, Vector< double > > & boundary_final_zeta_coordinate()
Return direct access to the final zeta coordinates of a boundary.
TriangleMeshCurveSection * curviline_to_polyline(TriangleMeshCurviLine *&curviline_pt, unsigned &bnd_id)
Helper function that creates the associated polyline representation for curvilines.
void set_geom_objects_and_coordinate_limits_for_close_curve(TriangleMeshClosedCurve *input_closed_curve_pt)
Stores the geometric objects associated to the curve sections that compound the closed curve....
std::map< unsigned, Vector< double > > Boundary_segment_final_zeta
Stores the final zeta coordinate for the segments that appear when a boundary is splited among proces...
std::map< unsigned, TriangleMeshCurveSection * > Boundary_curve_section_pt
A map that stores the associated curve section of the specified boundary id.
bool is_point_inside_polygon_helper(Vector< Vector< double >> polygon_vertices, Vector< double > point)
Helper function that checks if a given point is inside a polygon (a set of sorted vertices that conne...
void create_vertex_coordinates_for_polyline_connections(TriangleMeshCurviLine *boundary_pt, Vector< Vector< double >> &vertex_coord, Vector< std::pair< double, double >> &polygonal_vertex_arclength_info)
Helper function to create polyline vertex coordinates for curvilinear boundary specified by boundary_...
unsigned nregion_element(const unsigned &i)
Return the number of elements in the i-th region.
void create_vertex_coordinates_for_polyline_no_connections(TriangleMeshCurviLine *boundary_pt, Vector< Vector< double >> &vertex_coord, Vector< std::pair< double, double >> &polygonal_vertex_arclength_info)
Helper function to create polyline vertex coordinates for curvilinear boundary specified by boundary_...
std::map< unsigned, Vector< double > > & boundary_coordinate_limits()
Return access to the vector of boundary coordinates associated with each geometric object.
std::map< unsigned, Vector< double > > & boundary_initial_coordinate()
Return direct access to the initial coordinates of a boundary.
void flush_boundary_segment_node(const unsigned &b)
Flush the boundary segment node storage.
std::map< unsigned, Vector< Vector< double > > > & boundary_segment_final_coordinate()
Return direct access to the final coordinates for the segments that are part of a boundary.
unsigned long nboundary_segment_node(const unsigned &b, const unsigned &s)
Return the number of nodes associated with a given segment of a boundary.
void operator=(const UnstructuredTwoDMeshGeometryBase &)=delete
Broken assignment operator.
Vector< double > & boundary_segment_initial_zeta(const unsigned &b)
Return the initial zeta for the segments that are part of a boundary.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double Tolerable_error
Acceptable discrepancy for mismatch in vertex coordinates. In paranoid mode, the code will die if the...
void create_triangulateio_from_polyfiles(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TriangulateIO &triangle_io, bool &use_attributes)
Create a triangulateio data file from ele node and poly files. This is used if the mesh is generated ...
void dump_triangulateio(TriangulateIO &triangle_io, std::ostream &dump_file)
Write all the triangulateio data to disk in a dump file that can then be used to restart simulations.
void initialise_triangulateio(TriangulateIO &triangle_io)
Initialise TriangulateIO structure.
void read_triangulateio(std::istream &restart_file, TriangulateIO &triangle_io)
Read the triangulateio data from a dump file on disk, which can then be used to restart simulations.
void write_triangulateio_to_polyfile(TriangulateIO &triangle_io, std::ostream &poly_file)
Write the triangulateio data to disk as a poly file, mainly used for debugging.
void clear_triangulateio(TriangulateIO &triangulate_io, const bool &clear_hole_data)
Clear TriangulateIO structure.
TriangulateIO deep_copy_of_triangulateio_representation(TriangulateIO &triangle_io, const bool &quiet)
Make (partial) deep copy of TriangulateIO object. We only copy those items we need within oomph-lib's...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
FiniteElement * FE_pt
Pointer to bulk finite element.
The Triangle data structure, modified from the triangle.h header supplied with triangle 1....
double * pointlist
Pointer to list of points x coordinate followed by y coordinate.
int * pointmarkerlist
Pointer to list of point markers.
double * pointattributelist
Pointer to list of point attributes.
Data structure to store the base vertex info, initial or final vertex in the polylines have an associ...
Data structure filled when the connection matrix is created, for each polyline, there are two vertex_...