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