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