unstructured_two_d_mesh_geometry_base.cc
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
28
29
30/// ////////////////////////////////////////////////////////////////////
31/// ////////////////////////////////////////////////////////////////////
32/// ////////////////////////////////////////////////////////////////////
33
34namespace oomph
35{
36#ifdef OOMPH_HAS_TRIANGLE_LIB
37
38 //==================================================================
39 /// Helper namespace for triangle meshes
40 //==================================================================
41 namespace TriangleHelper
42 {
43 /// Clear TriangulateIO structure
44 void clear_triangulateio(TriangulateIO& triangulate_io,
45 const bool& clear_hole_data)
46 {
47 // Clear the point,attribute and marker list
48 free(triangulate_io.pointlist);
49 free(triangulate_io.pointattributelist);
50 free(triangulate_io.pointmarkerlist);
51 triangulate_io.numberofpoints = 0;
52 triangulate_io.numberofpointattributes = 0;
53
54 // Clear the triangle, attribute,neighbor and area list
55 free(triangulate_io.trianglelist);
56 free(triangulate_io.triangleattributelist);
57 free(triangulate_io.trianglearealist);
58 free(triangulate_io.neighborlist);
59 triangulate_io.numberoftriangles = 0;
60 triangulate_io.numberofcorners = 0;
61 triangulate_io.numberoftriangleattributes = 0;
62
63 // Clear the segment and marker list
64 free(triangulate_io.segmentlist);
65 free(triangulate_io.segmentmarkerlist);
66 triangulate_io.numberofsegments = 0;
67
68 // Clear hole list
69 if (clear_hole_data) free(triangulate_io.holelist);
70 triangulate_io.numberofholes = 0;
71
72 // Clear region list
73 if (clear_hole_data)
74 {
75 free(triangulate_io.regionlist);
76 }
77 triangulate_io.numberofregions = 0;
78
79 // Clear edge, marker and norm list
80 free(triangulate_io.edgelist);
81 free(triangulate_io.edgemarkerlist);
82 free(triangulate_io.normlist);
83 triangulate_io.numberofedges = 0;
84
85 // Now null it all out again
86 initialise_triangulateio(triangulate_io);
87 }
88
89
90 /// Initialise TriangulateIO structure
92 {
93 // Initialize the point list
94 triangle_io.pointlist = (double*)NULL;
95 triangle_io.pointattributelist = (double*)NULL;
96 triangle_io.pointmarkerlist = (int*)NULL;
97 triangle_io.numberofpoints = 0;
98 triangle_io.numberofpointattributes = 0;
99
100 // Initialize the triangle list
101 triangle_io.trianglelist = (int*)NULL;
102 triangle_io.triangleattributelist = (double*)NULL;
103 triangle_io.trianglearealist = (double*)NULL;
104 triangle_io.neighborlist = (int*)NULL;
105 triangle_io.numberoftriangles = 0;
106 triangle_io.numberofcorners = 0;
107 triangle_io.numberoftriangleattributes = 0;
108
109 // Initialize the segment list
110 triangle_io.segmentlist = (int*)NULL;
111 triangle_io.segmentmarkerlist = (int*)NULL;
112 triangle_io.numberofsegments = 0;
113
114 // Initialise hole list
115 triangle_io.holelist = (double*)NULL;
116 triangle_io.numberofholes = 0;
117
118 // Initialize region list
119 triangle_io.regionlist = (double*)NULL;
120 triangle_io.numberofregions = 0;
121
122 // Initialize edge list
123 triangle_io.edgelist = (int*)NULL;
124 triangle_io.edgemarkerlist = (int*)NULL;
125 triangle_io.normlist = (double*)NULL;
126 triangle_io.numberofedges = 0;
127 }
128
129
130 /// Make (partial) deep copy of TriangulateIO object. We only copy
131 /// those items we need within oomph-lib's adaptation procedures.
132 /// Warnings are issued if triangulate_io contains data that is not
133 /// not copied, unless quiet=true;
135 TriangulateIO& triangle_io, const bool& quiet)
136 {
137 // Create the struct
138 TriangulateIO triangle_out;
139
140 // Initialise
141 initialise_triangulateio(triangle_out);
142
143 // Point data
144 triangle_out.numberofpoints = triangle_io.numberofpoints;
145 triangle_out.pointlist =
146 (double*)malloc(triangle_out.numberofpoints * 2 * sizeof(double));
147 for (int j = 0; j < triangle_out.numberofpoints * 2; j++)
148 {
149 triangle_out.pointlist[j] = triangle_io.pointlist[j];
150 }
151
152 triangle_out.pointmarkerlist =
153 (int*)malloc(triangle_out.numberofpoints * sizeof(int));
154 for (int j = 0; j < triangle_out.numberofpoints; j++)
155 {
156 triangle_out.pointmarkerlist[j] = triangle_io.pointmarkerlist[j];
157 }
158
159 // Warn about laziness...
160 if (!quiet)
161 {
162 if ((triangle_io.pointattributelist != 0) ||
163 (triangle_io.numberofpointattributes != 0))
164 {
166 "Point attributes are not currently copied across",
167 "TriangleHelper::deep_copy_of_triangulateio_representation",
168 OOMPH_EXCEPTION_LOCATION);
169 }
170 }
171
172
173 // Triangle data
174 triangle_out.numberoftriangles = triangle_io.numberoftriangles;
175 triangle_out.trianglelist =
176 (int*)malloc(triangle_out.numberoftriangles * 3 * sizeof(int));
177 for (int j = 0; j < triangle_out.numberoftriangles * 3; j++)
178 {
179 triangle_out.trianglelist[j] = triangle_io.trianglelist[j];
180 }
181
182
183 // Copy over the triangle attribute data
184 triangle_out.numberoftriangleattributes =
185 triangle_io.numberoftriangleattributes;
186 triangle_out.triangleattributelist = (double*)malloc(
187 triangle_out.numberoftriangles *
188 triangle_out.numberoftriangleattributes * sizeof(double));
189 for (int j = 0; j < (triangle_out.numberoftriangles *
190 triangle_out.numberoftriangleattributes);
191 ++j)
192 {
193 triangle_out.triangleattributelist[j] =
194 triangle_io.triangleattributelist[j];
195 }
196
197
198 // Warn about laziness...
199 if (!quiet)
200 {
201 /* if ((triangle_io.triangleattributelist!=0)||
202 (triangle_io.numberoftriangleattributes!=0))
203 {
204 OomphLibWarning(
205 "Triangle attributes are not currently copied across",
206 "TriangleHelper::deep_copy_of_triangulateio_representation",
207 OOMPH_EXCEPTION_LOCATION);
208 }*/
209
210 if ((triangle_io.trianglearealist != 0))
211 {
213 "Triangle areas are not currently copied across",
214 "TriangleHelper::deep_copy_of_triangulateio_representation",
215 OOMPH_EXCEPTION_LOCATION);
216 }
217
218 if ((triangle_io.neighborlist != 0))
219 {
221 "Triangle neighbours are not currently copied across",
222 "TriangleHelper::deep_copy_of_triangulateio_representation",
223 OOMPH_EXCEPTION_LOCATION);
224 }
225 }
226
227
228 triangle_out.numberofcorners = triangle_io.numberofcorners;
229
230 // Segment data
231 triangle_out.numberofsegments = triangle_io.numberofsegments;
232 triangle_out.segmentlist =
233 (int*)malloc(triangle_out.numberofsegments * 2 * sizeof(int));
234 for (int j = 0; j < triangle_out.numberofsegments * 2; j++)
235 {
236 triangle_out.segmentlist[j] = triangle_io.segmentlist[j];
237 }
238 triangle_out.segmentmarkerlist =
239 (int*)malloc(triangle_out.numberofsegments * sizeof(int));
240 for (int j = 0; j < triangle_out.numberofsegments; j++)
241 {
242 triangle_out.segmentmarkerlist[j] = triangle_io.segmentmarkerlist[j];
243 }
244
245
246 // Region data
247 triangle_out.numberofregions = triangle_io.numberofregions;
248 triangle_out.regionlist =
249 (double*)malloc(triangle_out.numberofregions * 4 * sizeof(double));
250 for (int j = 0; j < triangle_out.numberofregions * 4; ++j)
251 {
252 triangle_out.regionlist[j] = triangle_io.regionlist[j];
253 }
254
255 // Hole data
256 triangle_out.numberofholes = triangle_io.numberofholes;
257 triangle_out.holelist =
258 (double*)malloc(triangle_out.numberofholes * 2 * sizeof(double));
259 for (int j = 0; j < triangle_out.numberofholes * 2; j++)
260 {
261 triangle_out.holelist[j] = triangle_io.holelist[j];
262 }
263
264
265 // Warn about laziness...
266 if (!quiet)
267 {
268 /* if ((triangle_io.regionlist!=0)||
269 (triangle_io.numberofregions!=0))
270 {
271 OomphLibWarning(
272 "Regions are not currently copied across",
273 "TriangleHelper::deep_copy_of_triangulateio_representation",
274 OOMPH_EXCEPTION_LOCATION);
275 }*/
276
277 if ((triangle_io.edgelist != 0) || (triangle_io.numberofedges != 0))
278 {
280 "Edges are not currently copied across",
281 "TriangleHelper::deep_copy_of_triangulateio_representation",
282 OOMPH_EXCEPTION_LOCATION);
283 }
284
285 if ((triangle_io.edgemarkerlist != 0))
286 {
288 "Edge markers are not currently copied across",
289 "TriangleHelper::deep_copy_of_triangulateio_representation",
290 OOMPH_EXCEPTION_LOCATION);
291 }
292
293 if ((triangle_io.normlist != 0))
294 {
296 "Normals are not currently copied across",
297 "TriangleHelper::deep_copy_of_triangulateio_representation",
298 OOMPH_EXCEPTION_LOCATION);
299 }
300 }
301
302 // Give it back!
303 return triangle_out;
304 }
305
306 /// Write the triangulateio data to disk as a poly file,
307 /// mainly used for debugging
309 std::ostream& poly_file)
310 {
311 // Up the precision dramatiacally
312 poly_file.precision(20);
313
314 // Output the number of points and their attributes
315 // Store the number of attributes
316 const int n_attr = triangle_io.numberofpointattributes;
317 poly_file << triangle_io.numberofpoints << " " << 2 << " " << n_attr
318 << " ";
319 // Determine whether there are point markers
320 bool point_markers = true;
321 if (triangle_io.pointmarkerlist == NULL)
322 {
323 point_markers = false;
324 }
325 // Indicate this in the file
326 poly_file << point_markers << "\n";
327
328 // Now output the point data
329 poly_file << "#Points\n";
330 for (int n = 0; n < triangle_io.numberofpoints; ++n)
331 {
332 // Output the point number and x and y coordinates
333 poly_file << n + 1 << " " << triangle_io.pointlist[2 * n] << " "
334 << triangle_io.pointlist[2 * n + 1] << " ";
335 // Output any attributes
336 for (int i = 0; i < n_attr; ++i)
337 {
338 poly_file << triangle_io.pointattributelist[n_attr * n + i] << " ";
339 }
340 // Output the boundary marker
341 if (point_markers)
342 {
343 poly_file << triangle_io.pointmarkerlist[n] << " ";
344 }
345 poly_file << "\n";
346 }
347
348 // Now move onto the segments
349 poly_file << "#Lines\n";
350 poly_file << triangle_io.numberofsegments << " ";
351 // Determine whether there are segment markers
352 bool seg_markers = true;
353 if (triangle_io.segmentmarkerlist == NULL)
354 {
355 seg_markers = false;
356 }
357 // Output this info in the file
358 poly_file << seg_markers << "\n";
359
360 // Now output the segment data
361 for (int n = 0; n < triangle_io.numberofsegments; ++n)
362 {
363 poly_file << n + 1 << " " << triangle_io.segmentlist[2 * n] << " "
364 << triangle_io.segmentlist[2 * n + 1] << " ";
365 // If there is a boundary marker output
366 if (seg_markers)
367 {
368 poly_file << triangle_io.segmentmarkerlist[n] << " ";
369 }
370 poly_file << "\n";
371 }
372
373 // Now output the number of holes
374 poly_file << "#No holes\n";
375 poly_file << triangle_io.numberofholes << "\n";
376 // Output the hole data
377 for (int h = 0; h < triangle_io.numberofholes; ++h)
378 {
379 poly_file << h + 1 << " " << triangle_io.holelist[2 * h] << " "
380 << triangle_io.holelist[2 * h + 1] << "\n";
381 }
382
383 // Now output the number of regions
384 poly_file << "#Assignment of attributes to regions\n";
385 poly_file << triangle_io.numberofregions << "\n";
386
387 // Loop over the regions
388 for (int r = 0; r < triangle_io.numberofregions; ++r)
389 {
390 poly_file << r + 1 << " ";
391 for (unsigned i = 0; i < 4; i++)
392 {
393 poly_file << triangle_io.regionlist[4 * r + i] << " ";
394 }
395 poly_file << "\n";
396 }
397 }
398
399
400 /// Create a triangulateio data file from ele node and poly
401 /// files. This is used if the mesh is generated by using Triangle
402 /// externally. The triangulateio structure is required to dump the mesh
403 /// topology for restarts.
405 const std::string& node_file_name,
406 const std::string& element_file_name,
407 const std::string& poly_file_name,
408 TriangulateIO& triangle_io,
409 bool& use_attributes)
410 {
411 // Initialise the TriangulateIO data structure
412 initialise_triangulateio(triangle_io);
413
414 // Process element file
415 std::ifstream element_file(element_file_name.c_str(), std::ios_base::in);
416
417 // Check that the file actually opened correctly
418 if (!element_file.is_open())
419 {
420 std::string error_msg("Failed to open element file: ");
421 error_msg += "\"" + element_file_name + "\".";
422 throw OomphLibError(
423 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
424 }
425
426 // Read in the number of elements
427 element_file >> triangle_io.numberoftriangles;
428 const unsigned n_element =
429 static_cast<unsigned>(triangle_io.numberoftriangles);
430
431 // Read in the number of nodes per element
432 element_file >> triangle_io.numberofcorners;
433 const unsigned n_local_node =
434 static_cast<unsigned>(triangle_io.numberofcorners);
435
436 // Read in the element attributes
437 element_file >> triangle_io.numberoftriangleattributes;
438 const unsigned n_attributes =
439 static_cast<unsigned>(triangle_io.numberoftriangleattributes);
440
441 // Allocate storage in the data structure
442 triangle_io.trianglelist =
443 (int*)malloc(triangle_io.numberoftriangles *
444 triangle_io.numberofcorners * sizeof(int));
445
446 if (n_attributes > 0)
447 {
448 triangle_io.triangleattributelist = (double*)malloc(
449 triangle_io.numberoftriangles *
450 triangle_io.numberoftriangleattributes * sizeof(double));
451 }
452
453 // Dummy storage
454 int dummy_element_number;
455
456 // Initialise counter
457 unsigned counter = 0;
458 unsigned counter2 = 0;
459
460 // Read global node numbers for all elements
461 for (unsigned e = 0; e < n_element; e++)
462 {
463 element_file >> dummy_element_number;
464 for (unsigned j = 0; j < n_local_node; j++)
465 {
466 element_file >> triangle_io.trianglelist[counter];
467 ++counter;
468 }
469 for (unsigned j = 0; j < n_attributes; j++)
470 {
471 element_file >> triangle_io.triangleattributelist[counter2];
472 ++counter2;
473 }
474 }
475 // Close the element file
476 element_file.close();
477
478 // Process node file
479 // -----------------
480 std::ifstream node_file(node_file_name.c_str(), std::ios_base::in);
481
482 // Check that the file actually opened correctly
483 if (!node_file.is_open())
484 {
485 std::string error_msg("Failed to open node file: ");
486 error_msg += "\"" + node_file_name + "\".";
487 throw OomphLibError(
488 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
489 }
490
491 // Read number of nodes
492 node_file >> triangle_io.numberofpoints;
493 unsigned n_node = static_cast<unsigned>(triangle_io.numberofpoints);
494
495 // Spatial dimension of nodes
496 unsigned dimension;
497 node_file >> dimension;
498
499#ifdef PARANOID
500 if (dimension != 2)
501 {
502 throw OomphLibError("The dimension must be 2\n",
503 OOMPH_CURRENT_FUNCTION,
504 OOMPH_EXCEPTION_LOCATION);
505 }
506#endif
507
508 // Flag for attributes
509 node_file >> triangle_io.numberofpointattributes;
510 unsigned n_point_attributes =
511 static_cast<unsigned>(triangle_io.numberofpointattributes);
512
513 // Flag for boundary markers
514 unsigned boundary_markers_flag;
515 node_file >> boundary_markers_flag;
516
517 // Allocate storage
518 triangle_io.pointlist =
519 (double*)malloc(triangle_io.numberofpoints * 2 * sizeof(double));
520 triangle_io.pointattributelist =
521 (double*)malloc(triangle_io.numberofpoints *
522 triangle_io.numberofpointattributes * sizeof(double));
523 if (boundary_markers_flag)
524 {
525 triangle_io.pointmarkerlist =
526 (int*)malloc(triangle_io.numberofpoints * sizeof(int));
527 }
528
529 // Dummy for node number
530 unsigned dummy_node_number;
531
532 // Reset counter
533 counter = 0;
534 // Load in nodal posititions, point attributes
535 // and boundary markers
536 for (unsigned i = 0; i < n_node; i++)
537 {
538 node_file >> dummy_node_number;
539 node_file >> triangle_io.pointlist[2 * i];
540 node_file >> triangle_io.pointlist[2 * i + 1];
541 for (unsigned j = 0; j < n_point_attributes; ++j)
542 {
543 node_file >> triangle_io.pointattributelist[counter];
544 ++counter;
545 }
546 if (boundary_markers_flag)
547 {
548 node_file >> triangle_io.pointmarkerlist[i];
549 }
550 }
551 node_file.close();
552
553
554 // Process poly file to extract edges
555 //-----------------------------------
556
557 // Open poly file
558 std::ifstream poly_file(poly_file_name.c_str(), std::ios_base::in);
559
560 // Check that the file actually opened correctly
561 if (!poly_file.is_open())
562 {
563 std::string error_msg("Failed to open poly file: ");
564 error_msg += "\"" + poly_file_name + "\".";
565 throw OomphLibError(
566 error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
567 }
568
569 // Number of nodes in poly file --- these will be ignore
570 unsigned n_node_poly;
571 poly_file >> n_node_poly;
572
573 // Dimension
574 poly_file >> dimension;
575
576 // Attribute flag
577 unsigned attribute_flag;
578 poly_file >> attribute_flag;
579
580 // Boundary markers flag
581 poly_file >> boundary_markers_flag;
582
583
584 // Ignore node information: Note: No, we can't extract the
585 // actual nodes themselves from here!
586 unsigned dummy;
587 for (unsigned i = 0; i < n_node_poly; i++)
588 {
589 // Read in (and discard) node number and x and y coordinates
590 poly_file >> dummy;
591 poly_file >> dummy;
592 poly_file >> dummy;
593 // read in the attributes
594 for (unsigned j = 0; j < attribute_flag; ++j)
595 {
596 poly_file >> dummy;
597 }
598 // read in the boundary marker
599 if (boundary_markers_flag == 1)
600 {
601 poly_file >> dummy;
602 }
603 }
604
605 // Now extract the segment information
606 //------------------------------------
607
608 // Number of segments
609 poly_file >> triangle_io.numberofsegments;
610 unsigned n_segment = static_cast<unsigned>(triangle_io.numberofsegments);
611
612 // Boundary marker flag
613 poly_file >> boundary_markers_flag;
614
615 // Allocate storage
616 triangle_io.segmentlist =
617 (int*)malloc(triangle_io.numberofsegments * 2 * sizeof(int));
618 if (boundary_markers_flag)
619 {
620 triangle_io.segmentmarkerlist =
621 (int*)malloc(triangle_io.numberofsegments * sizeof(int));
622 }
623
624 // Dummy for global segment number
625 unsigned dummy_segment_number;
626
627 // Extract information for each segment
628 for (unsigned i = 0; i < n_segment; i++)
629 {
630 poly_file >> dummy_segment_number;
631 poly_file >> triangle_io.segmentlist[2 * i];
632 poly_file >> triangle_io.segmentlist[2 * i + 1];
633 if (boundary_markers_flag)
634 {
635 poly_file >> triangle_io.segmentmarkerlist[i];
636 }
637 }
638
639 // Extract hole center information
640 poly_file >> triangle_io.numberofholes;
641 unsigned n_hole = static_cast<unsigned>(triangle_io.numberofholes);
642
643 // Allocate memory
644 triangle_io.holelist =
645 (double*)malloc(triangle_io.numberofholes * 2 * sizeof(double));
646
647
648 // Dummy for hole number
649 unsigned dummy_hole;
650 // Loop over the holes to get centre coords
651 for (unsigned ihole = 0; ihole < n_hole; ihole++)
652 {
653 // Read the centre value
654 poly_file >> dummy_hole;
655 poly_file >> triangle_io.holelist[2 * ihole];
656 poly_file >> triangle_io.holelist[2 * ihole + 1];
657 }
658
659 // Extract regions information
660 poly_file >> triangle_io.numberofregions;
661 unsigned n_regions = static_cast<unsigned>(triangle_io.numberofregions);
662
663 // Allocate memory
664 triangle_io.regionlist =
665 (double*)malloc(triangle_io.numberofregions * 4 * sizeof(double));
666
667 // Check for using regions
668 if (n_regions > 0)
669 {
670 use_attributes = true;
671 }
672
673 // Dummy for regions number
674 unsigned dummy_region;
675
676 // Loop over the regions to get their coords
677 for (unsigned iregion = 0; iregion < n_regions; iregion++)
678 {
679 // Read the regions coordinates
680 poly_file >> dummy_region;
681 poly_file >> triangle_io.regionlist[4 * iregion];
682 poly_file >> triangle_io.regionlist[4 * iregion + 1];
683 poly_file >> triangle_io.regionlist[4 * iregion + 2];
684 triangle_io.regionlist[4 * iregion + 3] = 0.0;
685
686 // Skip the rest of the line, there is no need to read the size of
687 // the elements in the region since that value is no longer used
688 poly_file.ignore(80, '\n');
689
690 // Verify if not using the default region number (zero)
691 if (triangle_io.regionlist[4 * iregion + 2] == 0)
692 {
693 std::ostringstream error_message;
694 error_message
695 << "Please use another region id different from zero.\n"
696 << "It is internally used as the default region number.\n";
697 throw OomphLibError(error_message.str(),
698 OOMPH_CURRENT_FUNCTION,
699 OOMPH_EXCEPTION_LOCATION);
700 }
701 }
702
703 poly_file.close();
704 }
705
706
707 /// Write all the triangulateio data to disk in a dump file
708 /// that can then be used to restart simulations
709 void dump_triangulateio(TriangulateIO& triangle_io, std::ostream& dump_file)
710 {
711 // Dump the triangles first
712 dump_file << triangle_io.numberoftriangles
713 << " # number of elements in TriangulateIO" << std::endl;
714
715 dump_file << triangle_io.numberofcorners
716 << " # number of nodes in each triangle" << std::endl;
717
718 dump_file << triangle_io.numberoftriangleattributes
719 << " # number of triangle attributes" << std::endl;
720
721 // Loop over and dump the triangle information
722 const int n_element = triangle_io.numberoftriangles;
723 const int n_local_node = triangle_io.numberofcorners;
724 const int n_attribute = triangle_io.numberoftriangleattributes;
725 unsigned counter = 0, counter2 = 0;
726 for (int e = 0; e < n_element; ++e)
727 {
728 // Dump the corners
729 dump_file << e << " # element number " << std::endl;
730 for (int n = 0; n < n_local_node; ++n)
731 {
732 dump_file << triangle_io.trianglelist[counter] << std::endl;
733 ++counter;
734 }
735 // Dump the attributes
736 dump_file << n_attribute << " # number of attributes" << std::endl;
737 for (int n = 0; n < n_attribute; ++n)
738 {
739 dump_file << triangle_io.triangleattributelist[counter2] << std::endl;
740 ++counter2;
741 }
742 }
743
744
745 // Dump the points (nodes) next
746 dump_file << triangle_io.numberofpoints
747 << " # number of points in TriangulateIO" << std::endl;
748 dump_file << triangle_io.numberofpointattributes
749 << " # number of point attributes" << std::endl;
750 // Test whether there are point markers
751 bool point_marker_flag = true;
752 if (triangle_io.pointmarkerlist == NULL)
753 {
754 point_marker_flag = false;
755 }
756 dump_file << point_marker_flag << " # point marker flag" << std::endl;
757
758
759 // Now output the point data
760 const int n_nodes = triangle_io.numberofpoints;
761 const int n_point_attributes = triangle_io.numberofpointattributes;
762 counter = 0;
763 counter2 = 0;
764 for (int n = 0; n < n_nodes; ++n)
765 {
766 dump_file << n << " # point number " << std::endl;
767 for (int i = 0; i < 2; ++i)
768 {
769 dump_file << triangle_io.pointlist[counter] << std::endl;
770 ++counter;
771 }
772 dump_file << n_point_attributes << " # number of point attributes "
773 << std::endl;
774 // Output any attributes
775 for (int i = 0; i < n_point_attributes; ++i)
776 {
777 dump_file << triangle_io.pointattributelist[counter2] << std::endl;
778 ++counter2;
779 }
780 dump_file << point_marker_flag << " # point marker flag " << std::endl;
781 // Output the boundary marker
782 if (point_marker_flag)
783 {
784 dump_file << triangle_io.pointmarkerlist[n] << std::endl;
785 }
786 }
787
788 // Now move onto the segments
789 dump_file << triangle_io.numberofsegments
790 << " # Number of segments in TriangulateIO " << std::endl;
791
792 // Determine whether there are segment markers
793 bool seg_marker_flag = true;
794 if (triangle_io.segmentmarkerlist == NULL)
795 {
796 seg_marker_flag = false;
797 }
798 // Output this info in the file
799 dump_file << seg_marker_flag << " # segment marker flag " << std::endl;
800
801 const int n_segments = triangle_io.numberofsegments;
802 counter = 0;
803 // Now output the segment data
804 for (int n = 0; n < n_segments; ++n)
805 {
806 dump_file << n << " # segment number " << std::endl;
807 for (int i = 0; i < 2; ++i)
808 {
809 dump_file << triangle_io.segmentlist[counter] << std::endl;
810 ++counter;
811 }
812
813 // If there is a boundary marker output
814 dump_file << seg_marker_flag << " # segment marker flag " << std::endl;
815 if (seg_marker_flag)
816 {
817 dump_file << triangle_io.segmentmarkerlist[n] << std::endl;
818 }
819 }
820
821 // Now output the number of holes
822 dump_file << triangle_io.numberofholes << " # number of holes "
823 << std::endl;
824 const int n_hole = triangle_io.numberofholes;
825 // Output the hole data
826 for (int h = 0; h < n_hole; ++h)
827 {
828 dump_file << h << " # hole number " << std::endl;
829 dump_file << triangle_io.holelist[2 * h] << std::endl;
830 dump_file << triangle_io.holelist[2 * h + 1] << std::endl;
831 }
832
833 // Now output the number of regions
834 dump_file << triangle_io.numberofregions << " # number of regions "
835 << std::endl;
836
837 const int n_region = triangle_io.numberofregions;
838 // Loop over the regions
839 counter = 0;
840 for (int r = 0; r < n_region; ++r)
841 {
842 dump_file << r << " # region number " << std::endl;
843 for (unsigned i = 0; i < 4; i++)
844 {
845 dump_file << triangle_io.regionlist[counter] << std::endl;
846 ++counter;
847 }
848 }
849 }
850
851 /// Read the triangulateio data from a dump file on
852 /// disk, which can then be used to restart simulations
853 void read_triangulateio(std::istream& restart_file,
854 TriangulateIO& triangle_io)
855 {
856 // String for reading
857 std::string input_string;
858
859 // Initialise the triangulate data structure
860 initialise_triangulateio(triangle_io);
861
862 // Read the first line up to termination sign
863 getline(restart_file, input_string, '#');
864 // Ignore the rest of the line
865 restart_file.ignore(80, '\n');
866 // Convert the number
867 triangle_io.numberoftriangles = atoi(input_string.c_str());
868
869 // Read the next line up to termination sign
870 getline(restart_file, input_string, '#');
871 // Ignore the rest of the line
872 restart_file.ignore(80, '\n');
873 // Convert the number
874 triangle_io.numberofcorners = atoi(input_string.c_str());
875
876 // Read the next line up to termination sign
877 getline(restart_file, input_string, '#');
878 // Ignore the rest of the line
879 restart_file.ignore(80, '\n');
880 // Convert the number
881 triangle_io.numberoftriangleattributes = atoi(input_string.c_str());
882
883 // Convert numbers into register variables
884 const int n_element = triangle_io.numberoftriangles;
885 const int n_local_node = triangle_io.numberofcorners;
886 const int n_attribute = triangle_io.numberoftriangleattributes;
887
888 // Allocate storage in the data structure
889 triangle_io.trianglelist =
890 (int*)malloc(triangle_io.numberoftriangles *
891 triangle_io.numberofcorners * sizeof(int));
892
893 if (n_attribute > 0)
894 {
895 triangle_io.triangleattributelist = (double*)malloc(
896 triangle_io.numberoftriangles *
897 triangle_io.numberoftriangleattributes * sizeof(double));
898 }
899
900 // Loop over elements and load in data
901 unsigned counter = 0, counter2 = 0;
902 for (int e = 0; e < n_element; ++e)
903 {
904 // Read the next line and ignore it
905 getline(restart_file, input_string);
906 for (int n = 0; n < n_local_node; ++n)
907 {
908 getline(restart_file, input_string);
909 triangle_io.trianglelist[counter] = atoi(input_string.c_str());
910 ++counter;
911 }
912 // Read the attributes
913 getline(restart_file, input_string);
914 for (int n = 0; n < n_attribute; ++n)
915 {
916 getline(restart_file, input_string);
917 triangle_io.triangleattributelist[counter2] =
918 atof(input_string.c_str());
919 ++counter2;
920 }
921 }
922
923
924 // Read the points (nodes) next up to termination string
925 getline(restart_file, input_string, '#');
926 // ignore the rest
927 restart_file.ignore(80, '\n');
928 triangle_io.numberofpoints = atoi(input_string.c_str());
929
930 // Read the point attributes next up to termination string
931 getline(restart_file, input_string, '#');
932 // ignore the rest
933 restart_file.ignore(80, '\n');
934 triangle_io.numberofpointattributes = atoi(input_string.c_str());
935
936 // Read whether there are point markers
937 getline(restart_file, input_string, '#');
938 // ignore the rest
939 restart_file.ignore(80, '\n');
940 int point_marker_flag = atoi(input_string.c_str());
941
942 // Allocate storage
943 triangle_io.pointlist =
944 (double*)malloc(triangle_io.numberofpoints * 2 * sizeof(double));
945 triangle_io.pointattributelist =
946 (double*)malloc(triangle_io.numberofpoints *
947 triangle_io.numberofpointattributes * sizeof(double));
948 if (point_marker_flag)
949 {
950 triangle_io.pointmarkerlist =
951 (int*)malloc(triangle_io.numberofpoints * sizeof(int));
952 }
953
954
955 // Now read the point data
956 const int n_nodes = triangle_io.numberofpoints;
957 const int n_point_attributes = triangle_io.numberofpointattributes;
958 counter = 0;
959 counter2 = 0;
960 for (int n = 0; n < n_nodes; ++n)
961 {
962 // Ignore the first line
963 getline(restart_file, input_string);
964 // Get the positions
965 for (int i = 0; i < 2; ++i)
966 {
967 getline(restart_file, input_string);
968 triangle_io.pointlist[counter] = atof(input_string.c_str());
969 ++counter;
970 }
971
972 // Ignore the next line about point attributes
973 getline(restart_file, input_string);
974
975 // Read any attributes
976 for (int i = 0; i < n_point_attributes; ++i)
977 {
978 getline(restart_file, input_string);
979 triangle_io.pointattributelist[counter2] = atof(input_string.c_str());
980 ++counter2;
981 }
982
983 // Ignore the next line
984 getline(restart_file, input_string);
985 // Output the boundary marker
986 if (point_marker_flag)
987 {
988 getline(restart_file, input_string);
989 triangle_io.pointmarkerlist[n] = atoi(input_string.c_str());
990 }
991 }
992
993 // Next read the segments
994 getline(restart_file, input_string, '#');
995 restart_file.ignore(80, '\n');
996 triangle_io.numberofsegments = atoi(input_string.c_str());
997
998 // Determine whether there are segment markers
999 getline(restart_file, input_string, '#');
1000 // ignore the rest
1001 restart_file.ignore(80, '\n');
1002 int seg_marker_flag = atoi(input_string.c_str());
1003
1004 // Allocate storage
1005 triangle_io.segmentlist =
1006 (int*)malloc(triangle_io.numberofsegments * 2 * sizeof(int));
1007 if (seg_marker_flag)
1008 {
1009 triangle_io.segmentmarkerlist =
1010 (int*)malloc(triangle_io.numberofsegments * sizeof(int));
1011 }
1012
1013 const int n_segments = triangle_io.numberofsegments;
1014 counter = 0;
1015 // Now output the segment data
1016 for (int n = 0; n < n_segments; ++n)
1017 {
1018 getline(restart_file, input_string);
1019 // get input
1020 for (int i = 0; i < 2; ++i)
1021 {
1022 getline(restart_file, input_string);
1023 triangle_io.segmentlist[counter] = atoi(input_string.c_str());
1024 ++counter;
1025 }
1026
1027 // Ignore the next line
1028 getline(restart_file, input_string);
1029 if (seg_marker_flag)
1030 {
1031 getline(restart_file, input_string);
1032 triangle_io.segmentmarkerlist[n] = atoi(input_string.c_str());
1033 }
1034 }
1035
1036 // Now read the holes
1037 getline(restart_file, input_string, '#');
1038 restart_file.ignore(80, '\n');
1039 triangle_io.numberofholes = atoi(input_string.c_str());
1040
1041 // Allocate memory
1042 triangle_io.holelist =
1043 (double*)malloc(triangle_io.numberofholes * 2 * sizeof(double));
1044
1045 const int n_hole = triangle_io.numberofholes;
1046 // Output the hole data
1047 for (int h = 0; h < n_hole; ++h)
1048 {
1049 // Ignore the first line
1050 getline(restart_file, input_string);
1051 // get the centre
1052 getline(restart_file, input_string);
1053 triangle_io.holelist[2 * h] = atof(input_string.c_str());
1054 getline(restart_file, input_string);
1055 triangle_io.holelist[2 * h + 1] = atof(input_string.c_str());
1056 }
1057
1058 // Now read the number of regions
1059 getline(restart_file, input_string, '#');
1060 restart_file.ignore(80, '\n');
1061 triangle_io.numberofregions = atoi(input_string.c_str());
1062
1063 const int n_region = triangle_io.numberofregions;
1064
1065 // Allocate memory
1066 triangle_io.regionlist =
1067 (double*)malloc(triangle_io.numberofregions * 4 * sizeof(double));
1068
1069 // Loop over the regions
1070 counter = 0;
1071 for (int r = 0; r < n_region; ++r)
1072 {
1073 getline(restart_file, input_string);
1074 for (unsigned i = 0; i < 4; i++)
1075 {
1076 getline(restart_file, input_string);
1077 triangle_io.regionlist[counter] = atof(input_string.c_str());
1078 ++counter;
1079 }
1080 }
1081 }
1082 } // namespace TriangleHelper
1083
1084#endif
1085
1086
1087 /// /////////////////////////////////////////////////////////////////
1088 /// /////////////////////////////////////////////////////////////////
1089 /// /////////////////////////////////////////////////////////////////
1090
1091
1092 /// Namespace that allows the specification of a tolerance
1093 /// between vertices at the ends of polylines that are supposed
1094 /// to be at the same position.
1095 namespace ToleranceForVertexMismatchInPolygons
1096 {
1097 /// Acceptable discrepancy for mismatch in vertex coordinates.
1098 /// In paranoid mode, the code will die if the beginning/end of
1099 /// two adjacent polylines differ by more than that. If the
1100 /// discrepancy is smaller (but nonzero) one of the vertex coordinates
1101 /// get adjusted to match perfectly; without paranoia the vertex
1102 /// coordinates are taken as they come...
1103 double Tolerable_error = 1.0e-14;
1104
1105 } // namespace ToleranceForVertexMismatchInPolygons
1106
1107
1108 /// /////////////////////////////////////////////////////////////////
1109 /// /////////////////////////////////////////////////////////////////
1110 /// /////////////////////////////////////////////////////////////////
1111
1112
1113 // =======================================================================
1114 // Connects the initial vertex of the curve section to a desired
1115 /// target polyline by specifying the vertex number. There is a checking
1116 /// which verifies that the initial vertex is close enough to the
1117 /// destination vertex on the target polyline by no more than the specified
1118 /// tolerance
1119 // =======================================================================
1121 TriangleMeshPolyLine* polyline_pt,
1122 const unsigned& vertex_number,
1123 const double& tolerance_for_connection)
1124 {
1125#ifdef PARANOID
1126 unsigned n_vertices = polyline_pt->nvertex();
1127
1128 if (n_vertices <= vertex_number)
1129 {
1130 std::ostringstream error_stream;
1131 error_stream << "The vertex number you provided (" << vertex_number
1132 << ") is greater\n than the number of vertices ("
1133 << n_vertices << "in the specified TriangleMeshPolyLine.\n"
1134 << "Remember that the vertex index starts at 0" << std::endl
1135 << "Source boundary (" << boundary_id()
1136 << ") wants to connect "
1137 << "to destination boundary (" << polyline_pt->boundary_id()
1138 << ")" << std::endl;
1139 throw OomphLibError(
1140 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1141 }
1142
1143 // Verify if there is really a match point in the specified
1144 // connection values
1145 Vector<double> v_src(2);
1146 Vector<double> v_dst(2);
1147
1148 this->initial_vertex_coordinate(v_src);
1149 v_dst = polyline_pt->vertex_coordinate(vertex_number);
1150
1151 double error = sqrt((v_src[0] - v_dst[0]) * (v_src[0] - v_dst[0]) +
1152 (v_src[1] - v_dst[1]) * (v_src[1] - v_dst[1]));
1153
1154 if (error > tolerance_for_connection)
1155 {
1156 std::ostringstream error_stream;
1157 error_stream << "The associated vertices for the connection"
1158 << "\nare not close enough. Their respective values are:\n"
1159 << "Source boundary id:(" << this->boundary_id() << ")\n"
1160 << "Source vertex x:(" << v_src[0] << ") y:(" << v_src[1]
1161 << ")\n"
1162 << "Destination boundary id:(" << polyline_pt->boundary_id()
1163 << ")"
1164 << "\nAssociated vertex x:(" << v_dst[0] << ") y:("
1165 << v_dst[1] << ")"
1166 << "\nThe corresponding distance is: (" << error
1167 << ") but the "
1168 << "allowed\ntolerance is: (" << tolerance_for_connection
1169 << ")" << std::endl;
1170 throw OomphLibError(
1171 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1172 }
1173
1174#endif
1175
1178 Initial_vertex_connected_n_vertex = vertex_number;
1180 }
1181
1182 // =======================================================================
1183 // Connects the final vertex of the curve section to a desired
1184 /// target polyline by specifying the vertex number. There is a checking
1185 /// which verifies that the final vertex is close enough to the
1186 /// destination vertex on the target polyline by no more than the specified
1187 /// tolerance
1188 // =======================================================================
1190 TriangleMeshPolyLine* polyline_pt,
1191 const unsigned& vertex_number,
1192 const double& tolerance_for_connection)
1193 {
1194#ifdef PARANOID
1195 unsigned n_vertices = polyline_pt->nvertex();
1196
1197 if (n_vertices <= vertex_number)
1198 {
1199 std::ostringstream error_stream;
1200 error_stream << "The vertex number you provided (" << vertex_number
1201 << ") is greater\n than the number of vertices ("
1202 << n_vertices << "in the specified TriangleMeshPolyLine.\n"
1203 << "Remember that the vertex index starts at 0" << std::endl
1204 << "Source boundary (" << boundary_id()
1205 << ") wants to connect "
1206 << "to destination boundary (" << polyline_pt->boundary_id()
1207 << ")" << std::endl;
1208 throw OomphLibError(
1209 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1210 }
1211
1212 // Verify if there is really a match point in the specified
1213 // connection values
1214 Vector<double> v_src(2);
1215 Vector<double> v_dst(2);
1216
1217 this->final_vertex_coordinate(v_src);
1218 v_dst = polyline_pt->vertex_coordinate(vertex_number);
1219
1220 double error = sqrt((v_src[0] - v_dst[0]) * (v_src[0] - v_dst[0]) +
1221 (v_src[1] - v_dst[1]) * (v_src[1] - v_dst[1]));
1222
1223 if (error > tolerance_for_connection)
1224 {
1225 std::ostringstream error_stream;
1226 error_stream << "The associated vertices for the connection"
1227 << "\nare not close enough. Their respective values are:\n"
1228 << "Source boundary id:(" << this->boundary_id() << ")\n"
1229 << "Source vertex x:(" << v_src[0] << ") y:(" << v_src[1]
1230 << ")\n"
1231 << "Destination boundary id:(" << polyline_pt->boundary_id()
1232 << ")"
1233 << "\nAssociated vertex x:(" << v_dst[0] << ") y:("
1234 << v_dst[1] << ")"
1235 << "\nThe corresponding distance is: (" << error
1236 << ") but the "
1237 << "allowed\ntolerance is: (" << tolerance_for_connection
1238 << ")" << std::endl;
1239 throw OomphLibError(
1240 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1241 }
1242
1243#endif
1244
1247 Final_vertex_connected_n_vertex = vertex_number;
1249 }
1250
1251 // =======================================================================
1252 // Connects the initial vertex of the curve section to a desired
1253 /// target curviline by specifying the s value (intrinsic value on the
1254 /// geometric object of the curviline) where to connect on the target
1255 /// curviline. There is a checking which verifies that the initial vertex
1256 /// and the coordinates on the given s value are close enough by no more
1257 /// than the given tolerance
1258 // =======================================================================
1260 TriangleMeshCurviLine* curviline_pt,
1261 const double& s_value,
1262 const double& tolerance_for_connection)
1263 {
1264#ifdef PARANOID
1265 double z_initial = curviline_pt->zeta_start();
1266 double z_final = curviline_pt->zeta_end();
1267 double z_max = std::max(z_initial, z_final);
1268 double z_min = std::min(z_initial, z_final);
1269 if (s_value < z_min || z_max < s_value)
1270 {
1271 std::ostringstream error_stream;
1272 error_stream << "The s value you provided for connection (" << s_value
1273 << ") is out\nof the limits of the specified "
1274 << "TriangleMeshCurviLine.\nThe limits are [" << z_initial
1275 << ", " << z_final << "]" << std::endl
1276 << "Source boundary (" << boundary_id()
1277 << ") wants to connect "
1278 << "to destination boundary (" << curviline_pt->boundary_id()
1279 << ")" << std::endl;
1280 throw OomphLibError(
1281 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1282 }
1283
1284 // Verify if there is really a match point in the specified
1285 // connection values
1286 Vector<double> v_src(2);
1287 Vector<double> v_dst(2);
1288 Vector<double> z(1);
1289
1290 this->initial_vertex_coordinate(v_src);
1291 z[0] = s_value;
1292 curviline_pt->geom_object_pt()->position(z, v_dst);
1293 double error = sqrt((v_src[0] - v_dst[0]) * (v_src[0] - v_dst[0]) +
1294 (v_src[1] - v_dst[1]) * (v_src[1] - v_dst[1]));
1295 if (error >= tolerance_for_connection)
1296 {
1297 std::ostringstream error_stream;
1298 error_stream
1299 << "The associated vertex for the provided connection s value\n"
1300 << "are not close enough. Their respective values are:\n"
1301 << "Source boundary id:(" << this->boundary_id() << ")\n"
1302 << "Source vertex x:(" << v_src[0] << ") y:(" << v_src[1] << ")\n"
1303 << "Destination boundary id:(" << curviline_pt->boundary_id() << ")"
1304 << "\nDestination s value (" << s_value << ")\n"
1305 << "Associated vertex x:(" << v_dst[0] << ") y:(" << v_dst[1] << ")"
1306 << "\nThe corresponding distance is: (" << error << ") but the "
1307 << "allowed\ntolerance is: (" << tolerance_for_connection << ")"
1308 << std::endl;
1309 throw OomphLibError(
1310 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1311 }
1312
1313#endif
1314
1319
1320 // We are still not able to compute the vertex number but we can
1321 // at least store the corresponding s value
1322 // The corresponding vertex will be computed when the curviline be
1323 // changed into a polyline
1325 Tolerance_for_s_connection = tolerance_for_connection;
1326
1327 curviline_pt->add_connection_point(s_value, tolerance_for_connection);
1328 }
1329
1330 // =======================================================================
1331 // Connects the final vertex of the curve section to a desired
1332 /// target curviline by specifying the s value (intrinsic value on the
1333 /// geometric object of the curviline) where to connect on the target
1334 /// curviline. There is a checking which verifies that the final vertex
1335 /// and the coordinates on the given s value are close enough by no more
1336 /// than the given tolerance
1337 // =======================================================================
1339 TriangleMeshCurviLine* curviline_pt,
1340 const double& s_value,
1341 const double& tolerance_for_connection)
1342 {
1343#ifdef PARANOID
1344 double z_initial = curviline_pt->zeta_start();
1345 double z_final = curviline_pt->zeta_end();
1346 double z_max = std::max(z_initial, z_final);
1347 double z_min = std::min(z_initial, z_final);
1348 if (s_value < z_min || z_max < s_value)
1349 {
1350 std::ostringstream error_stream;
1351 error_stream << "The s value you provided for connection (" << s_value
1352 << ") is out\nof the limits of the specified "
1353 << "TriangleMeshCurviLine.\nThe limits are [" << z_initial
1354 << ", " << z_final << "]" << std::endl
1355 << "Source boundary (" << boundary_id()
1356 << ") wants to connect "
1357 << "to destination boundary (" << curviline_pt->boundary_id()
1358 << ")" << std::endl;
1359 throw OomphLibError(
1360 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1361 }
1362
1363 // Verify if there is really a match point in the specified
1364 // connection values
1365 Vector<double> v_src(2);
1366 Vector<double> v_dst(2);
1367 Vector<double> z(1);
1368
1369 this->final_vertex_coordinate(v_src);
1370 z[0] = s_value;
1371 curviline_pt->geom_object_pt()->position(z, v_dst);
1372
1373 double error = sqrt((v_src[0] - v_dst[0]) * (v_src[0] - v_dst[0]) +
1374 (v_src[1] - v_dst[1]) * (v_src[1] - v_dst[1]));
1375
1376 if (error >= tolerance_for_connection)
1377 {
1378 std::ostringstream error_stream;
1379 error_stream
1380 << "The associated vertex for the provided connection s value\n"
1381 << "are not close enough. Their respective values are:\n"
1382 << "Source boundary id:(" << this->boundary_id() << ")\n"
1383 << "Source vertex x:(" << v_src[0] << ") y:(" << v_src[1] << ")\n"
1384 << "Destination boundary id:(" << curviline_pt->boundary_id() << ")"
1385 << "\nDestination s value (" << s_value << ")\n"
1386 << "Associated vertex x:(" << v_dst[0] << ") y:(" << v_dst[1] << ")"
1387 << "\nThe corresponding distance is: (" << error << ") but the "
1388 << "allowed\ntolerance is: (" << tolerance_for_connection << ")"
1389 << std::endl;
1390 throw OomphLibError(
1391 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1392 }
1393
1394#endif
1395
1400
1401 // We are still not able to compute the vertex number but we can
1402 // at least store the corresponding s value.
1403 // The corresponding vertex will be computed when the curviline be
1404 // transformed into a polyline
1405 Final_s_connection_value = s_value;
1406 Tolerance_for_s_connection = tolerance_for_connection;
1407
1408 curviline_pt->add_connection_point(s_value, tolerance_for_connection);
1409 }
1410
1411
1412 /// ////////////////////////////////////////////////////////////////////
1413 /// ////////////////////////////////////////////////////////////////////
1414 /// ///////////////////////////////////////////////////////////////////
1415
1416
1417 //=====================================================================
1418 /// Class defining a closed curve for the Triangle mesh generation
1419 //=====================================================================
1421 const Vector<TriangleMeshCurveSection*>& curve_section_pt,
1422 const Vector<double>& internal_point_pt,
1423 const bool& is_internal_point_fixed)
1424 : TriangleMeshCurve(curve_section_pt),
1425 Internal_point_pt(internal_point_pt),
1426 Is_internal_point_fixed(is_internal_point_fixed)
1427 {
1428 // Matching of curve sections i.e. the last vertex of the i curve
1429 // section should match with the first vertex of the i+1 curve
1430 // section
1431
1432 // Total number of boundaries
1433 const unsigned n_boundaries = Curve_section_pt.size();
1434
1435 // Need at least two
1436 if (n_boundaries < 2)
1437 {
1438 std::ostringstream error_stream;
1439 error_stream << "Sorry -- I'm afraid we insist that a closed curve is\n"
1440 << "specified by at least two separate CurveSections.\n"
1441 << "You've only specified (" << n_boundaries << ")"
1442 << std::endl;
1443 throw OomphLibError(
1444 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1445 }
1446
1447 // Check last point of each boundary bit coincides with first point
1448 // on next one
1449 for (unsigned i = 0; i < n_boundaries - 1; i++)
1450 {
1451 // Auxiliary vertex for storing the vertex values of contiguous curves
1452 Vector<double> v1(2);
1453
1454 // This is for getting the final coordinates of the i curve section
1456
1457 // Auxiliary vertex for storing the vertex values of contiguous curves
1458 Vector<double> v2(2);
1459
1460 // This is for the start coordinates of the i+1 curve section
1462
1463 // Work out error
1464 double error = sqrt(pow(v1[0] - v2[0], 2) + pow(v1[1] - v2[1], 2));
1465
1467 {
1468 std::ostringstream error_stream;
1469 error_stream
1470 << "The start and end points of curve section boundary parts\n"
1471 << i << " and " << i + 1
1472 << " don't match when judged with the tolerance of "
1474 << " which\nis specified in the namespace variable\n"
1475 << "ToleranceForVertexMismatchInPolygons::Tolerable_error.\n\n"
1476 << "These are the vertices coordinates:\n"
1477 << "Curve section (" << i << ") final vertex: (" << v1[0] << ", "
1478 << v1[1] << ")\n"
1479 << "Curve section (" << i + 1 << ") initial vertex: (" << v2[0]
1480 << ", " << v2[1] << ")\n"
1481 << "The distance between the vertices is (" << error << ")\n"
1482 << "Feel free to adjust this or to recompile the code without\n"
1483 << "paranoia if you think this is OK...\n"
1484 << std::endl;
1485 throw OomphLibError(
1486 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1487 }
1488 else
1489 {
1490 // Aligns (only implemented for polylines)
1491 TriangleMeshPolyLine* current_polyline =
1492 dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i]);
1493 TriangleMeshPolyLine* next_polyline =
1494 dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i + 1]);
1495
1496 // Was it able to do the cast?
1497 if (current_polyline && next_polyline)
1498 {
1499 unsigned last_vertex = current_polyline->nvertex() - 1;
1500 next_polyline->vertex_coordinate(0) =
1501 current_polyline->vertex_coordinate(last_vertex);
1502 }
1503 }
1504
1505 } // For n_boundaries - 1
1506
1507 // Check wrap around
1508 // Auxiliary vertex for storing the vertex values of contiguous curves
1509 Vector<double> v1(2);
1510
1511 // This is for getting the first coordinates of the first curve section
1512 Curve_section_pt[0]->initial_vertex_coordinate(v1);
1513
1514 // Auxiliary vertex for storing the vertex values of contiguous curves
1515 Vector<double> v2(2);
1516
1517 // This is for getting the last coordinates of the last curve section
1518 Curve_section_pt[n_boundaries - 1]->final_vertex_coordinate(v2);
1519
1520 double error = sqrt(pow(v2[0] - v1[0], 2) + pow(v2[1] - v1[1], 2));
1521
1523 {
1524 std::ostringstream error_stream;
1525 error_stream
1526 << "The start and end points of the first and last curve segment\n"
1527 << "boundary parts don't match when judged \nwith the tolerance of "
1529 << " which is specified in the namespace \nvariable "
1530 << "ToleranceForVertexMismatchInPolygons::Tolerable_error.\n\n"
1531 << "Feel free to adjust this or to recompile the code without\n"
1532 << "paranoia if you think this is OK...\n"
1533 << std::endl;
1534 throw OomphLibError(
1535 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1536 }
1537 else
1538 {
1539 // Aligns (only implemented for polylines)
1540 TriangleMeshPolyLine* first_polyline =
1541 dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[0]);
1542 TriangleMeshPolyLine* last_polyline =
1543 dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[n_boundaries - 1]);
1544
1545 // Was it able to do the cast?
1546 if (first_polyline && last_polyline)
1547 {
1548 unsigned last_vertex = last_polyline->nvertex() - 1;
1549 first_polyline->vertex_coordinate(0) =
1550 last_polyline->vertex_coordinate(last_vertex);
1551 }
1552 }
1553 }
1554
1555
1556 /// //////////////////////////////////////////////////////////////////
1557 /// //////////////////////////////////////////////////////////////////
1558 /// //////////////////////////////////////////////////////////////////
1559
1560
1561 //=========================================================================
1562 /// Constructor: Specify vector of pointers to TriangleMeshCurveSection
1563 /// that define the boundary of the segments of the polygon.
1564 /// Each TriangleMeshCurveSection has its own boundary ID and can contain
1565 /// multiple (straight-line) segments. For consistency across the
1566 /// various uses of this class, we insist that the closed boundary
1567 /// is represented by at least two separate TriangleMeshCurveSection
1568 /// whose joint vertices must be specified in both.
1569 /// (This is to allow the setup of unique boundary coordinate(s)
1570 /// around the polygon.) This may seem slightly annoying
1571 /// in cases where a polygon really only represents a single
1572 /// boundary, but...
1573 /// Note: The specified vector of pointers must consist of only
1574 /// TriangleMeshPolyLine elements. There is a checking on the PARANOID
1575 /// mode for this constraint
1576 //=========================================================================
1578 const Vector<TriangleMeshCurveSection*>& boundary_polyline_pt,
1579 const Vector<double>& internal_point_pt,
1580 const bool& is_internal_point_fixed)
1581 : TriangleMeshCurve(boundary_polyline_pt),
1583 boundary_polyline_pt, internal_point_pt, is_internal_point_fixed),
1584 Enable_redistribution_of_segments_between_polylines(false),
1585 Can_update_configuration(false),
1586 Polygon_fixed(false)
1587 {
1588 // Get the number of polylines
1589 const unsigned n_bound = boundary_polyline_pt.size();
1590
1591 // Check that all the constituent TriangleMeshCurveSection are
1592 // instance of TriangleMeshPolyLine
1593 for (unsigned p = 0; p < n_bound; p++)
1594 {
1595 TriangleMeshPolyLine* tmp_polyline_pt =
1596 dynamic_cast<TriangleMeshPolyLine*>(boundary_polyline_pt[p]);
1597 if (tmp_polyline_pt == 0)
1598 {
1599 std::ostringstream error_stream;
1600 error_stream << "The (" << p << ") TriangleMeshCurveSection is not a "
1601 << "TriangleMeshPolyLine.\nThe TriangleMeshPolygon object"
1602 << "is constituent of only TriangleMeshPolyLine objects.\n"
1603 << "Verify that all the constituent elements of the "
1604 << "TriangleMeshPolygon\nare instantiated as "
1605 << "TriangleMeshPolyLines." << std::endl;
1606 throw OomphLibError(
1607 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1608 }
1609 }
1610
1611 // Check that the polylines are contiguous
1612 bool contiguous = true;
1613 unsigned i_offensive = 0;
1614
1615 // Need at least two
1616 if (n_bound < 2)
1617 {
1618 std::ostringstream error_stream;
1619 error_stream
1620 << "Sorry -- I'm afraid we insist that a closed curve is\n"
1621 << "specified by at least two separate TriangleMeshPolyLines.\n"
1622 << "You've only specified (" << n_bound << ")" << std::endl;
1623 throw OomphLibError(error_stream.str(),
1624 "TriangleMeshPolygon::TriangleMeshPolygon()",
1625 OOMPH_EXCEPTION_LOCATION);
1626 } // if (n_bound<2)
1627
1628 // Does the last node of the polyline connect to the first one of the
1629 // next one (only up the last but one!)
1630 for (unsigned i = 0; i < n_bound - 1; i++)
1631 {
1632 // Get vector last vertex in current polyline
1633 unsigned last_vertex = (polyline_pt(i)->nvertex()) - 1;
1634 Vector<double> v1 = polyline_pt(i)->vertex_coordinate(last_vertex);
1635
1636 // Get vector to first vertex in next polyline
1638
1639 // Work out error
1640 double error = sqrt(pow(v1[0] - v2[0], 2) + pow(v1[1] - v2[1], 2));
1641
1642 // Is error accetable?
1644 {
1645 contiguous = false;
1646 i_offensive = i;
1647 break;
1648 }
1649 // Align
1650 else
1651 {
1653 polyline_pt(i)->vertex_coordinate(last_vertex);
1654 }
1655 }
1656
1657 // Does the last one connect to the first one?
1658
1659 // Get vector last vertex last polyline
1660 unsigned last_vertex = (polyline_pt(n_bound - 1)->nvertex()) - 1;
1661 Vector<double> v1 =
1662 polyline_pt(n_bound - 1)->vertex_coordinate(last_vertex);
1663
1664 // Get vector first vertex first polyline
1666 double error = sqrt(pow(v1[0] - v2[0], 2) + pow(v1[1] - v2[1], 2));
1667
1669 {
1670 contiguous = false;
1671 i_offensive = n_bound - 1;
1672 }
1673 else
1674 {
1676 polyline_pt(n_bound - 1)->vertex_coordinate(last_vertex);
1677 }
1678
1679 if (!contiguous)
1680 {
1681 std::ostringstream error_stream;
1682 error_stream
1683 << "The polylines specified \n"
1684 << "should define a closed geometry, i.e. the first/last vertex of\n"
1685 << "adjacent polylines should match.\n\n"
1686 << "Your polyline number " << i_offensive
1687 << " has no contiguous neighbour, when judged \nwith the tolerance of "
1689 << " which is specified in the namespace \nvariable "
1690 << "ToleranceForVertexMismatchInPolygons::Tolerable_error.\n\n"
1691 << "Feel free to adjust this or to recompile the code without\n"
1692 << "paranoia if you think this is OK...\n"
1693 << std::endl;
1694 throw OomphLibError(error_stream.str(),
1695 "TriangleMeshPolygon::TriangleMeshPolygon()",
1696 OOMPH_EXCEPTION_LOCATION);
1697 }
1698
1699 // Check if internal point is actually located in bounding polygon
1700 // Reference: http://paulbourke.net/geometry/insidepoly/
1701
1702 // Only checked if there is an internal hole
1703 if (!Internal_point_pt.empty())
1704 {
1705 // Vertex coordinates
1706 Vector<Vector<double>> polygon_vertex;
1707
1708 // Total number of vertices
1709 unsigned nvertex = 0;
1710
1711 // Storage for first/last point on polyline for contiguousness check
1712 Vector<double> last_vertex(2);
1713 Vector<double> first_vertex(2);
1714
1715 // Get vertices
1716 unsigned npolyline = boundary_polyline_pt.size();
1717 for (unsigned i = 0; i < npolyline; i++)
1718 {
1719 // Number of vertices
1720 unsigned nvert = boundary_polyline_pt[i]->nvertex();
1721 for (unsigned j = 0; j < nvert; j++)
1722 {
1723 // Check contiguousness
1724 if ((i > 1) && (j == 0))
1725 {
1726 first_vertex = polyline_pt(i)->vertex_coordinate(j);
1727 double dist = sqrt(pow(first_vertex[0] - last_vertex[0], 2) +
1728 pow(first_vertex[1] - last_vertex[1], 2));
1730 {
1731 std::ostringstream error_stream;
1732 error_stream
1733 << "The start and end points of polylines " << i << " and "
1734 << i + 1 << " don't match when judged\n"
1735 << "with the tolerance ("
1737 << ") which is specified in the namespace \nvariable "
1738 << "ToleranceForVertexMismatchInPolygons::"
1739 << "Tolerable_error.\n\n"
1740 << "Feel free to adjust this or to recompile the "
1741 << "code without\n"
1742 << "paranoia if you think this is OK...\n"
1743 << std::endl;
1744 throw OomphLibError(error_stream.str(),
1745 "TriangleMeshPolygon::TriangleMeshPolygon()",
1746 OOMPH_EXCEPTION_LOCATION);
1747 }
1748 }
1749 // Get vertex (ignore end point)
1750 if (j < nvert - 1)
1751 {
1752 polygon_vertex.push_back(polyline_pt(i)->vertex_coordinate(j));
1753 }
1754 // Prepare for check of contiguousness
1755 else
1756 {
1757 last_vertex = polyline_pt(i)->vertex_coordinate(j);
1758 }
1759 }
1760 }
1761
1762 // Total number of vertices
1763 nvertex = polygon_vertex.size();
1764
1765 // Counter for number of intersections
1766 unsigned intersect_counter = 0;
1767
1768 // Get first vertex
1769 Vector<double> p1 = polygon_vertex[0];
1770 for (unsigned i = 1; i <= nvertex; i++)
1771 {
1772 // Get second vertex by wrap-around
1773 Vector<double> p2 = polygon_vertex[i % nvertex];
1774
1775 if (Internal_point_pt[1] > std::min(p1[1], p2[1]))
1776 {
1777 if (Internal_point_pt[1] <= std::max(p1[1], p2[1]))
1778 {
1779 if (Internal_point_pt[0] <= std::max(p1[0], p2[0]))
1780 {
1781 if (p1[1] != p2[1])
1782 {
1783 double xintersect = (Internal_point_pt[1] - p1[1]) *
1784 (p2[0] - p1[0]) / (p2[1] - p1[1]) +
1785 p1[0];
1786 if ((p1[0] == p2[0]) || (Internal_point_pt[0] <= xintersect))
1787 {
1788 intersect_counter++;
1789 }
1790 }
1791 }
1792 }
1793 }
1794 p1 = p2;
1795 }
1796
1797 // Even number of intersections: outside
1798 if (intersect_counter % 2 == 0)
1799 {
1800 std::ostringstream error_stream;
1801 error_stream
1802 << "The internal point at " << Internal_point_pt[0] << " "
1803 << Internal_point_pt[1]
1804 << " isn't in the polygon that describes the internal closed "
1805 << "curve!\nPolygon vertices are at: \n";
1806 for (unsigned i = 0; i < nvertex; i++)
1807 {
1808 error_stream << polygon_vertex[i][0] << " " << polygon_vertex[i][1]
1809 << "\n";
1810 }
1811 error_stream
1812 << "This may be because the internal point is defined by a\n"
1813 << "GeomObject that has deformed so much that it's \n"
1814 << "swept over the (initial) internal point.\n"
1815 << "If so, you should update the position of the internal point. \n"
1816 << "This could be done automatically by generating \n"
1817 << "an internal mesh inside the polygon and using one\n"
1818 << "of its internal nodes as the internal point. Actually not \n"
1819 << "why triangle doesn't do that automatically....\n";
1820 throw OomphLibError(error_stream.str(),
1821 "TriangleMeshPolygon::TriangleMeshPolygon()",
1822 OOMPH_EXCEPTION_LOCATION);
1823 }
1824 }
1825 }
1826
1827
1828 /// //////////////////////////////////////////////////////////////////
1829 /// //////////////////////////////////////////////////////////////////
1830 /// //////////////////////////////////////////////////////////////////
1831
1832
1833 //=====================================================================
1834 /// Class defining an open curve for the Triangle mesh generation
1835 //=====================================================================
1837 const Vector<TriangleMeshCurveSection*>& curve_section_pt)
1838 : TriangleMeshCurve(curve_section_pt)
1839 {
1840 // Matching of curve sections i.e. the last vertex of
1841 // the i curve section should match with the first
1842 // vertex of the i+1 curve section
1843
1844 // Total number of boundaries
1845 unsigned n_boundaries = Curve_section_pt.size();
1846
1847 // Check last point of each boundary bit coincides with first point
1848 // on next one
1849 for (unsigned i = 0; i < n_boundaries - 1; i++)
1850 {
1851 // Auxiliary vertex for storing the vertex values of contiguous curves
1852 Vector<double> v1(2);
1853 Vector<double> v2(2);
1854
1855 // This is for getting the final coordinates of the i curve section
1856 Curve_section_pt[i]->final_vertex_coordinate(v1);
1857
1858 // This is for the start coordinates of the i+1 curve section
1859 Curve_section_pt[i + 1]->initial_vertex_coordinate(v2);
1860
1861 // Work out error
1862 double error = sqrt(pow(v1[0] - v2[0], 2) + pow(v1[1] - v2[1], 2));
1864 {
1865 std::ostringstream error_stream;
1866 error_stream
1867 << "The start and end points of curve section boundary parts " << i
1868 << " and " << i + 1
1869 << " don't match when judged \nwith the tolerance of "
1871 << " which is specified in the namespace \nvariable "
1872 << "ToleranceForVertexMismatchInPolygons::Tolerable_error.\n\n"
1873 << "Feel free to adjust this or to recompile the code without\n"
1874 << "paranoia if you think this is OK...\n"
1875 << std::endl;
1876 throw OomphLibError(
1877 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1878 }
1879 else
1880 {
1881 // Aligns (only implemented for polylines)
1882 TriangleMeshPolyLine* current_polyline =
1883 dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i]);
1884 TriangleMeshPolyLine* next_polyline =
1885 dynamic_cast<TriangleMeshPolyLine*>(Curve_section_pt[i + 1]);
1886
1887 if (current_polyline && next_polyline)
1888 {
1889 unsigned last_vertex = current_polyline->nvertex() - 1;
1890 next_polyline->vertex_coordinate(0) =
1891 current_polyline->vertex_coordinate(last_vertex);
1892 }
1893 }
1894
1895 } // For n_boundaries - 1
1896 }
1897
1898#ifdef OOMPH_HAS_TRIANGLE_LIB
1899
1900 //======================================================================
1901 /// Create TriangulateIO object from outer boundaries, internal
1902 /// boundaries, and open curves. Add the holes and regions
1903 /// information as well
1904 //======================================================================
1906 Vector<TriangleMeshPolygon*>& outer_polygons_pt,
1907 Vector<TriangleMeshPolygon*>& internal_polygons_pt,
1908 Vector<TriangleMeshOpenCurve*>& open_curves_pt,
1909 Vector<Vector<double>>& extra_holes_coordinates,
1910 std::map<unsigned, Vector<double>>& regions_coordinates,
1911 std::map<unsigned, double>& regions_areas,
1912 TriangulateIO& triangulate_io)
1913 {
1914 // These are the general stages of the algorithm
1915 // --------------------------------------------------------------------
1916 // 1) Create the boundary connections structure, every entry states
1917 // if the initial or final vertex is connected to other vertex
1918 // --------------------------------------------------------------------
1919 // 2) Create the structure for base vertices, the base vertices are
1920 // those that are not connected
1921 // --------------------------------------------------------------------
1922 // 3) Assign a unique id to each vertex in the polylines (creates a
1923 // look up scheme used to create the segments)
1924 // --------------------------------------------------------------------
1925 // 4) Create the segments using the unique id of the vertices and
1926 // assign a segment id to each segment (the one from the
1927 // polyline-boundary)
1928 // ------------------------------------------------------------------
1929 // 5) Fill the triangulateio containers with the collected
1930 // information
1931 // --------------------------------------------------------------
1932
1933 // ------------------------------------------------------------------
1934 // 1st- stage
1935
1936 // 1) Create the boundary connections structure for the outer
1937 // boundaries, internal boundaries and open curves
1938
1939 // Store the number of vertices on each boundary (which is
1940 // represented by a polyline -- for quick access--). We may have
1941 // more than one polyline with the same boundary id, that is why we
1942 // need a vector to represent the set of polylines associated to a
1943 // boundary (the chunks). The mentioned case is quite common when
1944 // working in parallel, when a boundary may be split because of the
1945 // distribution strategy
1946 std::map<unsigned, std::map<unsigned, unsigned>> boundary_chunk_n_vertices;
1947
1948 // Note: For each polyline, we only consider (v-1) vertices since
1949 // the first vertex of the p-th polyline (when p > 1) is the same as
1950 // the last vertex of the (p-1)-th polyline (when p > 1). KEEP THIS
1951 // ---ALWAYS--- IN MIND WHEN REVIEWING THE CODE
1952
1953 // The connections matrix
1954
1955 // Stores the vertex_connection_info:
1956 // - is_connected
1957 // - boundary_id_to_connect
1958 // - boundary_chunk_to_connect
1959 // - vertex_number_to_connect
1960 // of the initial and final vertex of each polyline
1961 // -----------------------
1962 // (boundary, chunk#, vertex (initial[0] or final[1])) ->
1963 // vertex_connection_info
1964 // -----------------------
1965 // map[x][][] = boundary_id
1966 // map[][x][] = chunk_id
1967 // Vector[][][x] = vertex#, only initial or final (that is why only
1968 // two indexes)
1969 std::map<unsigned, std::map<unsigned, Vector<vertex_connection_info>>>
1970 connection_matrix;
1971
1972 // Initialize the base vertex structure (every vertex is a not base
1973 // vertex by default)
1974
1975 // The data structure that stores the base vertex information
1976 // Stores the base_vertex_info:
1977 // - done
1978 // - base_vertex
1979 // - boundary_id
1980 // - boundary_chunk
1981 // - vertex_number
1982 // of the initial and final vertex of each polyline
1983 // -----------------------
1984 // (boundary, chunk#, vertex (initial[0] or final[1])) -> base_vertex_info
1985 // -----------------------
1986 // map[x][][] = boundary_id
1987 // map[][x][] = chunk_id
1988 // Vector[][][x] = vertex#, only initial or final (that is why only
1989 // two indexes)
1990 std::map<unsigned, std::map<unsigned, Vector<base_vertex_info>>>
1991 base_vertices;
1992
1993 // Get the number of outer polygons
1994 const unsigned n_outer_polygons = outer_polygons_pt.size();
1995
1996 for (unsigned i = 0; i < n_outer_polygons; i++)
1997 {
1998 // The number of polylines in the current polygon
1999 const unsigned n_polylines = outer_polygons_pt[i]->npolyline();
2000
2001 for (unsigned p = 0; p < n_polylines; p++)
2002 {
2003 // Get a pointer to the current polyline
2004 TriangleMeshPolyLine* tmp_polyline_pt =
2005 outer_polygons_pt[i]->polyline_pt(p);
2006
2007 // Get the boundary id of the current polyline
2008 const unsigned bound_id = tmp_polyline_pt->boundary_id();
2009
2010 // The number of vertices in the current polyline
2011 const unsigned n_vertices = tmp_polyline_pt->nvertex();
2012
2013 // Get the chunk number associated with this boundary
2014 const unsigned bound_chunk = tmp_polyline_pt->boundary_chunk();
2015
2016 // Store the number of vertices of the polyline in the global
2017 // storage
2018 boundary_chunk_n_vertices[bound_id][bound_chunk] = n_vertices;
2019
2020 // Get the next polyline (or the initial polyline)
2021 TriangleMeshPolyLine* next_polyline_pt = 0;
2022
2023 // Is there next polyline
2024 if (p < n_polylines - 1)
2025 {
2026 // Set the next polyline
2027 next_polyline_pt = outer_polygons_pt[i]->polyline_pt(p + 1);
2028 }
2029 else
2030 {
2031 // The next polyline is the initial polyline
2032 next_polyline_pt = outer_polygons_pt[i]->polyline_pt(0);
2033 }
2034
2035 // Add the information to the connections matrix
2037 tmp_polyline_pt, connection_matrix, next_polyline_pt);
2038
2039 // Initialise the base vertices structure for the current
2040 // polyline
2041 initialise_base_vertex(tmp_polyline_pt, base_vertices);
2042
2043 } // for (p < n_polylines)
2044
2045 } // for (i < n_outer_polygons)
2046
2047 // Get the number of internal polygons
2048 const unsigned n_internal_polygons = internal_polygons_pt.size();
2049
2050 for (unsigned i = 0; i < n_internal_polygons; i++)
2051 {
2052 // The number of polylines in the current polygon
2053 const unsigned n_polylines = internal_polygons_pt[i]->npolyline();
2054
2055 for (unsigned p = 0; p < n_polylines; p++)
2056 {
2057 // Get a pointer to the current polyline
2058 TriangleMeshPolyLine* tmp_polyline_pt =
2059 internal_polygons_pt[i]->polyline_pt(p);
2060
2061 // Get the boundary id of the current polyline
2062 const unsigned bound_id = tmp_polyline_pt->boundary_id();
2063
2064 // The number of vertices in the current polyline
2065 const unsigned n_vertices = tmp_polyline_pt->nvertex();
2066
2067 // Get the chunk number associated with this boundary
2068 const unsigned bound_chunk = tmp_polyline_pt->boundary_chunk();
2069
2070 // Store the number of vertices of the polyline in the global
2071 // storage
2072 boundary_chunk_n_vertices[bound_id][bound_chunk] = n_vertices;
2073
2074 // Get the next polyline (or the initial polyline)
2075 TriangleMeshPolyLine* next_polyline_pt = 0;
2076
2077 // Is there next polyline
2078 if (p < n_polylines - 1)
2079 {
2080 // Set the next polyline
2081 next_polyline_pt = internal_polygons_pt[i]->polyline_pt(p + 1);
2082 }
2083 else
2084 {
2085 // The next polyline is the initial polyline
2086 next_polyline_pt = internal_polygons_pt[i]->polyline_pt(0);
2087 }
2088
2089 // Add the information to the connections matrix
2091 tmp_polyline_pt, connection_matrix, next_polyline_pt);
2092
2093 // Initialise the base vertices structure for the current
2094 // polyline
2095 initialise_base_vertex(tmp_polyline_pt, base_vertices);
2096
2097 } // for (p < n_polylines)
2098
2099 } // for (i < n_internal_polygons)
2100
2101 // Get the number of open curves
2102 const unsigned n_open_curves = open_curves_pt.size();
2103
2104 for (unsigned i = 0; i < n_open_curves; i++)
2105 {
2106 // The number of polylines in the current polygon
2107 const unsigned n_polylines = open_curves_pt[i]->ncurve_section();
2108
2109 for (unsigned p = 0; p < n_polylines; p++)
2110 {
2111 // Get a pointer to the current polyline
2112 TriangleMeshPolyLine* tmp_polyline_pt =
2113 open_curves_pt[i]->polyline_pt(p);
2114
2115 // Get the boundary id of the current polyline
2116 const unsigned bound_id = tmp_polyline_pt->boundary_id();
2117
2118 // The number of vertices in the current polyline
2119 const unsigned n_vertices = tmp_polyline_pt->nvertex();
2120
2121 // Get the chunk number associated with this boundary
2122 const unsigned bound_chunk = tmp_polyline_pt->boundary_chunk();
2123
2124 // Store the number of vertices of the polyline in the global
2125 // storage
2126 boundary_chunk_n_vertices[bound_id][bound_chunk] = n_vertices;
2127
2128 // Get the next polyline (or the initial polyline)
2129 TriangleMeshPolyLine* next_polyline_pt = 0;
2130
2131 // Is there next polyline
2132 if (p < n_polylines - 1)
2133 {
2134 // Set the next polyline
2135 next_polyline_pt = open_curves_pt[i]->polyline_pt(p + 1);
2136 }
2137 else
2138 {
2139 // If we are in the last polyline of the open curve there is
2140 // no actual next polyline
2141 }
2142
2143 // Add the information to the connections matrix
2145 tmp_polyline_pt, connection_matrix, next_polyline_pt);
2146
2147 // Initialise the base vertices structure for the current
2148 // polyline
2149 initialise_base_vertex(tmp_polyline_pt, base_vertices);
2150
2151 } // for (p < n_polylines)
2152
2153 } // for (i < n_open_curves)
2154
2155 // ------------------------------------------------------------------
2156
2157 // ------------------------------------------------------------------
2158 // 2) Create the structure for base vertices, the base vertices are
2159 // those that are not connected
2160 // ------------------------------------------------------------------
2161
2162 // Loop over the polylines in the outer polygons and indentify the
2163 // base vertices
2164 for (unsigned i = 0; i < n_outer_polygons; i++)
2165 {
2166 // The number of polylines in the current polygon
2167 const unsigned n_polylines = outer_polygons_pt[i]->npolyline();
2168
2169 for (unsigned p = 0; p < n_polylines; p++)
2170 {
2171 // Get a pointer to the current polyline
2172 TriangleMeshPolyLine* tmp_polyline_pt =
2173 outer_polygons_pt[i]->polyline_pt(p);
2174
2175 // Identify the base vertices in the current polyline
2176 add_base_vertex_info_helper(tmp_polyline_pt,
2177 base_vertices,
2178 connection_matrix,
2179 boundary_chunk_n_vertices);
2180
2181 } // for (p < n_polylines)
2182
2183 } // for (i < n_outer_polygons)
2184
2185 // Loop over the polylines in the internal polygons and indentify the
2186 // base vertices
2187 for (unsigned i = 0; i < n_internal_polygons; i++)
2188 {
2189 // The number of polylines in the current polygon
2190 const unsigned n_polylines = internal_polygons_pt[i]->npolyline();
2191
2192 for (unsigned p = 0; p < n_polylines; p++)
2193 {
2194 // Get a pointer to the current polyline
2195 TriangleMeshPolyLine* tmp_polyline_pt =
2196 internal_polygons_pt[i]->polyline_pt(p);
2197
2198 // Identify the base vertices in the current polyline
2199 add_base_vertex_info_helper(tmp_polyline_pt,
2200 base_vertices,
2201 connection_matrix,
2202 boundary_chunk_n_vertices);
2203
2204 } // for (p < n_polylines)
2205
2206 } // for (i < n_internal_polygons)
2207
2208 // Loop over the polylines in the open curves and indentify the base
2209 // vertices
2210 for (unsigned i = 0; i < n_open_curves; i++)
2211 {
2212 // The number of polylines in the current polygon
2213 const unsigned n_polylines = open_curves_pt[i]->ncurve_section();
2214
2215 for (unsigned p = 0; p < n_polylines; p++)
2216 {
2217 // Get a pointer to the current polyline
2218 TriangleMeshPolyLine* tmp_polyline_pt =
2219 open_curves_pt[i]->polyline_pt(p);
2220
2221 // Identify the base vertices in the current polyline
2222 add_base_vertex_info_helper(tmp_polyline_pt,
2223 base_vertices,
2224 connection_matrix,
2225 boundary_chunk_n_vertices);
2226
2227 } // for (p < n_polylines)
2228
2229 } // for (i < n_open_curves)
2230
2231 // ------------------------------------------------------------------
2232
2233 // ------------------------------------------------------------------
2234 // 3) Assign a unique id to each vertex in the polylines (creates a
2235 // look up scheme used to create the segments)
2236 // ------------------------------------------------------------------
2237
2238 // Create the storage for the look-up scheme
2239 // (boundary_local, chunk#, vertex#) -> global_vertex_id
2240 // map[x][][] = boundary
2241 // map[][x][] = chunk_id
2242 // Vector[][][x] = vertex#
2243 std::map<unsigned, std::map<unsigned, Vector<int>>>
2244 boundary_chunk_vertex_to_global_vertex_id;
2245
2246 // Create an entry in the map for each boundary, then do the same
2247 // for the chunks and finally resize the container (Vector) to store
2248 // the vertices
2249 for (std::map<unsigned, std::map<unsigned, unsigned>>::iterator it =
2250 boundary_chunk_n_vertices.begin();
2251 it != boundary_chunk_n_vertices.end();
2252 it++)
2253 {
2254 // Get the boundary id
2255 const unsigned b = (*it).first;
2256
2257 // Now loop over the chunks
2258 for (std::map<unsigned, unsigned>::iterator itt = (*it).second.begin();
2259 itt != (*it).second.end();
2260 itt++)
2261 {
2262 // Get the chunk id
2263 const unsigned c = (*itt).first;
2264
2265 // Get the number of vertices associated with the boundary-chunk
2266 // and resize the container
2267 const unsigned n_vertices = boundary_chunk_n_vertices[b][c];
2268
2269 // Now create storage in the container and resize the vector that
2270 // stores the vertices ids. Initialize the entries to -1
2271 boundary_chunk_vertex_to_global_vertex_id[b][c].resize(n_vertices, -1);
2272
2273 } // Loop over the chunks
2274
2275 } // Loop over the boundaries
2276
2277 // Counter for the numbering of the global vertices
2278 unsigned global_vertex_number = 0;
2279
2280 // Container for the vertices
2281 Vector<Vector<double>> global_vertices;
2282
2283 // Go through all the vertices in the polylines and assign a unique
2284 // id only to the base vertices, any other vertex copy the unique id
2285 // from its base vertex
2286
2287 // The total number of added vertices in the outer polygons
2288 unsigned n_vertices_outer_polygons = 0;
2289
2290 // Start with the outer polygons
2291 for (unsigned i = 0; i < n_outer_polygons; i++)
2292 {
2293 // The number of polylines in the current polygon
2294 const unsigned n_polylines = outer_polygons_pt[i]->npolyline();
2295
2296 for (unsigned p = 0; p < n_polylines; p++)
2297 {
2298 // Get a pointer to the current polyline
2299 TriangleMeshPolyLine* tmp_polyline_pt =
2300 outer_polygons_pt[i]->polyline_pt(p);
2301
2302 // Get the boundary id of the current polyline
2303 const unsigned bound_id = tmp_polyline_pt->boundary_id();
2304
2305 // The number of vertices in the current polyline
2306 const unsigned n_vertices = tmp_polyline_pt->nvertex();
2307
2308 // Get the current chunk number of the polyline
2309 const unsigned bnd_chunk = tmp_polyline_pt->boundary_chunk();
2310
2311 // Assign a global vertex id to the initial vertex
2312 // -----------------------------------------------
2313
2314 // Get the base vertex information of the initial vertex
2315 base_vertex_info initial_base_vertex_info =
2316 base_vertices[bound_id][bnd_chunk][0];
2317
2318#ifdef PARANOID
2319 if (!initial_base_vertex_info.has_base_vertex_assigned)
2320 {
2321 std::ostringstream error_message;
2322 std::string output_string =
2323 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2324 error_message
2325 << "The initial vertex of the current polyline has no base\n"
2326 << "vertex assigned\n"
2327 << "Outer polygon number: (" << i << ")\n\n"
2328 << "Polyline number: (" << p << ")\n"
2329 << "Boundary id: (" << bound_id << ")\n"
2330 << "Boundary chunk: (" << bnd_chunk << ")\n";
2331 throw OomphLibError(
2332 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2333 } // if (!initial_base_vertex_info.has_base_vertex_assigned)
2334#endif
2335
2336 // The base vertex boundary id
2337 unsigned bvbi = initial_base_vertex_info.boundary_id;
2338 // The base vertex boundary chunkx
2339 unsigned bvbc = initial_base_vertex_info.boundary_chunk;
2340 // The vertex number of the base vertex
2341 unsigned bvvn = initial_base_vertex_info.vertex_number;
2342
2343 // Get the global vertex id of the base vertex
2344 int global_vertex_id =
2345 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn];
2346
2347 // Check if the global vertex id has been already assigned
2348 if (global_vertex_id != -1)
2349 {
2350 // If that is the case then copy the global vertex id in the
2351 // current slot
2352 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][0] =
2353 global_vertex_id;
2354 } // if (global_vertex_id != -1)
2355 else
2356 {
2357 // Assign a global vertex id to the base vertex
2358 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn] =
2359 global_vertex_number;
2360
2361 // ... and copy the value to the initial vertex
2362 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][0] =
2363 global_vertex_number;
2364
2365 // ... increase the counter for the "global_vertex_number"
2366 global_vertex_number++;
2367
2368 // Get the vertex
2369 Vector<double> vertex = tmp_polyline_pt->vertex_coordinate(0);
2370
2371 // ... and add it to the global vertex container
2372 global_vertices.push_back(vertex);
2373 }
2374
2375 // Is the initial vertex a base vertex
2376 if (initial_base_vertex_info.is_base_vertex)
2377 {
2378 // Increase the independent vertices counter
2379 n_vertices_outer_polygons++;
2380 }
2381
2382 // ------------------------------------------------------------
2383 // Now loop over the intermediate vertices and assign a unique
2384 // vertex id (all intermediate vertices are base vertices)
2385 for (unsigned v = 1; v < n_vertices - 1; v++)
2386 {
2387 // Get the global vertex id
2388 global_vertex_id =
2389 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v];
2390
2391 // Check if the global vertex id has been already assigned.
2392 // We do nothing if it has been already assigned
2393 // (global_vertex_id!=-1).
2394
2395 // If it has not been already assigned (global_vertex_id==-1)
2396 // then set a new global vertex number, and add the vertex to
2397 // the global vertices container
2398 if (global_vertex_id == -1)
2399 {
2400 // Set a value for the global vertex
2401 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v] =
2402 global_vertex_number;
2403
2404 // ... increase the counter for the "global_vertex_number"
2405 global_vertex_number++;
2406
2407 // Add the vertex to the global vertex container
2408 Vector<double> vertex = tmp_polyline_pt->vertex_coordinate(v);
2409 // Add the vertex to the global vertex container
2410 global_vertices.push_back(vertex);
2411
2412 } // if (global_vertex_id == -1)
2413
2414 // Increase the independent vertices counter
2415 n_vertices_outer_polygons++;
2416
2417 } // for (n_vertices-1)
2418
2419 // Assign a global vertex id to the final vertex
2420 // -----------------------------------------------
2421
2422 // Get the base vertex information of the final vertex
2423 base_vertex_info final_base_vertex_info =
2424 base_vertices[bound_id][bnd_chunk][1];
2425
2426#ifdef PARANOID
2427 if (!final_base_vertex_info.has_base_vertex_assigned)
2428 {
2429 std::ostringstream error_message;
2430 std::string output_string =
2431 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2432 error_message
2433 << "The final vertex of the current polyline has no base\n"
2434 << "vertex assigned\n"
2435 << "Outer polygon number: (" << i << ")\n\n"
2436 << "Polyline number: (" << p << ")\n"
2437 << "Boundary id: (" << bound_id << ")\n"
2438 << "Boundary chunk: (" << bnd_chunk << ")\n";
2439 throw OomphLibError(
2440 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2441 } // if (!final_base_vertex_info.has_base_vertex_assigned)
2442#endif
2443
2444 // The base vertex boundary id
2445 bvbi = final_base_vertex_info.boundary_id;
2446 // The base vertex boundary chunkx
2447 bvbc = final_base_vertex_info.boundary_chunk;
2448 // The vertex number of the base vertex
2449 bvvn = final_base_vertex_info.vertex_number;
2450
2451 // Get the global vertex id of the base vertex
2452 global_vertex_id =
2453 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn];
2454
2455 // Check if the global vertex id has been already assigned
2456 if (global_vertex_id != -1)
2457 {
2458 // If that is the case then copy the global vertex id in the
2459 // current slot
2460 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
2461 [n_vertices - 1] =
2462 global_vertex_id;
2463 } // if (global_vertex_id != -1)
2464 else
2465 {
2466 // Assign a global vertex id to the base vertex
2467 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn] =
2468 global_vertex_number;
2469
2470 // ... and copy the value to the final vertex
2471 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
2472 [n_vertices - 1] =
2473 global_vertex_number;
2474
2475 // ... increase the counter for the "global_vertex_number"
2476 global_vertex_number++;
2477
2478 // Get the vertex
2479 Vector<double> vertex =
2480 tmp_polyline_pt->vertex_coordinate(n_vertices - 1);
2481
2482 // ... and add it to the global vertex container
2483 global_vertices.push_back(vertex);
2484 }
2485
2486 // Is the final vertex a base vertex
2487 if (final_base_vertex_info.is_base_vertex)
2488 {
2489 // Increase the independent vertices counter
2490 n_vertices_outer_polygons++;
2491 }
2492
2493 } // for (p < n_polylines)
2494
2495 } // for (i < n_outer_polygons)
2496
2497 // The total number of added vertices in the internal polygons
2498 unsigned n_vertices_internal_polygons = 0;
2499
2500 // Do the internal polygons
2501 for (unsigned i = 0; i < n_internal_polygons; i++)
2502 {
2503 // The number of polylines in the current polygon
2504 const unsigned n_polylines = internal_polygons_pt[i]->npolyline();
2505
2506 for (unsigned p = 0; p < n_polylines; p++)
2507 {
2508 // Get a pointer to the current polyline
2509 TriangleMeshPolyLine* tmp_polyline_pt =
2510 internal_polygons_pt[i]->polyline_pt(p);
2511
2512 // Get the boundary id of the current polyline
2513 const unsigned bound_id = tmp_polyline_pt->boundary_id();
2514
2515 // The number of vertices in the current polyline
2516 const unsigned n_vertices = tmp_polyline_pt->nvertex();
2517
2518 // Get the current chunk number of the polyline
2519 const unsigned bnd_chunk = tmp_polyline_pt->boundary_chunk();
2520
2521 // Assign a global vertex id to the initial vertex
2522 // -----------------------------------------------
2523
2524 // Get the base vertex information of the initial vertex
2525 base_vertex_info initial_base_vertex_info =
2526 base_vertices[bound_id][bnd_chunk][0];
2527
2528#ifdef PARANOID
2529 if (!initial_base_vertex_info.has_base_vertex_assigned)
2530 {
2531 std::ostringstream error_message;
2532 std::string output_string =
2533 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2534 error_message
2535 << "The initial vertex of the current polyline has no base\n"
2536 << "vertex assigned\n"
2537 << "Internal polygon number: (" << i << ")\n\n"
2538 << "Polyline number: (" << p << ")\n"
2539 << "Boundary id: (" << bound_id << ")\n"
2540 << "Boundary chunk: (" << bnd_chunk << ")\n";
2541 throw OomphLibError(
2542 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2543 } // if (!initial_base_vertex_info.has_base_vertex_assigned)
2544#endif
2545
2546 // The base vertex boundary id
2547 unsigned bvbi = initial_base_vertex_info.boundary_id;
2548 // The base vertex boundary chunkx
2549 unsigned bvbc = initial_base_vertex_info.boundary_chunk;
2550 // The vertex number of the base vertex
2551 unsigned bvvn = initial_base_vertex_info.vertex_number;
2552
2553 // Get the global vertex id of the base vertex
2554 int global_vertex_id =
2555 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn];
2556
2557 // Check if the global vertex id has been already assigned
2558 if (global_vertex_id != -1)
2559 {
2560 // If that is the case then copy the global vertex id in the
2561 // current slot
2562 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][0] =
2563 global_vertex_id;
2564 } // if (global_vertex_id != -1)
2565 else
2566 {
2567 // Assign a global vertex id to the base vertex
2568 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn] =
2569 global_vertex_number;
2570
2571 // ... and copy the value to the initial vertex
2572 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][0] =
2573 global_vertex_number;
2574
2575 // ... increase the counter for the "global_vertex_number"
2576 global_vertex_number++;
2577
2578 // Get the vertex
2579 Vector<double> vertex = tmp_polyline_pt->vertex_coordinate(0);
2580
2581 // ... and add it to the global vertex container
2582 global_vertices.push_back(vertex);
2583 }
2584
2585 // Is the initial vertex a base vertex
2586 if (initial_base_vertex_info.is_base_vertex)
2587 {
2588 // Increase the independent vertices counter
2589 n_vertices_internal_polygons++;
2590 }
2591
2592 // ------------------------------------------------------------
2593 // Now loop over the intermediate vertices and assign a unique
2594 // vertex id (all intermediate vertices are base vertices)
2595 for (unsigned v = 1; v < n_vertices - 1; v++)
2596 {
2597 // Get the global vertex id
2598 global_vertex_id =
2599 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v];
2600
2601 // Check if the global vertex id has been already assigned.
2602 // We do nothing if it has been already assigned
2603 // (global_vertex_id!=-1).
2604
2605 // If it has not been already assigned (global_vertex_id==-1)
2606 // then set a new global vertex number, and add the vertex to
2607 // the global vertices container
2608 if (global_vertex_id == -1)
2609 {
2610 // Set a value for the global vertex
2611 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v] =
2612 global_vertex_number;
2613
2614 // ... increase the counter for the "global_vertex_number"
2615 global_vertex_number++;
2616
2617 // Add the vertex to the global vertex container
2618 Vector<double> vertex = tmp_polyline_pt->vertex_coordinate(v);
2619 // Add the vertex to the global vertex container
2620 global_vertices.push_back(vertex);
2621
2622 } // if (global_vertex_id == -1)
2623
2624 // Increase the independent vertices counter
2625 n_vertices_internal_polygons++;
2626
2627 } // for (n_vertices-1)
2628
2629 // Assign a global vertex id to the final vertex
2630 // -----------------------------------------------
2631
2632 // Get the base vertex information of the final vertex
2633 base_vertex_info final_base_vertex_info =
2634 base_vertices[bound_id][bnd_chunk][1];
2635
2636#ifdef PARANOID
2637 if (!final_base_vertex_info.has_base_vertex_assigned)
2638 {
2639 std::ostringstream error_message;
2640 std::string output_string =
2641 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2642 error_message
2643 << "The final vertex of the current polyline has no base\n"
2644 << "vertex assigned\n"
2645 << "Internal polygon number: (" << i << ")\n\n"
2646 << "Polyline number: (" << p << ")\n"
2647 << "Boundary id: (" << bound_id << ")\n"
2648 << "Boundary chunk: (" << bnd_chunk << ")\n";
2649 throw OomphLibError(
2650 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2651 } // if (!final_base_vertex_info.has_base_vertex_assigned)
2652#endif
2653
2654 // The base vertex boundary id
2655 bvbi = final_base_vertex_info.boundary_id;
2656 // The base vertex boundary chunkx
2657 bvbc = final_base_vertex_info.boundary_chunk;
2658 // The vertex number of the base vertex
2659 bvvn = final_base_vertex_info.vertex_number;
2660
2661 // Get the global vertex id of the base vertex
2662 global_vertex_id =
2663 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn];
2664
2665 // Check if the global vertex id has been already assigned
2666 if (global_vertex_id != -1)
2667 {
2668 // If that is the case then copy the global vertex id in the
2669 // current slot
2670 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
2671 [n_vertices - 1] =
2672 global_vertex_id;
2673 } // if (global_vertex_id != -1)
2674 else
2675 {
2676 // Assign a global vertex id to the base vertex
2677 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn] =
2678 global_vertex_number;
2679
2680 // ... and copy the value to the final vertex
2681 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
2682 [n_vertices - 1] =
2683 global_vertex_number;
2684
2685 // ... increase the counter for the "global_vertex_number"
2686 global_vertex_number++;
2687
2688 // Get the vertex
2689 Vector<double> vertex =
2690 tmp_polyline_pt->vertex_coordinate(n_vertices - 1);
2691
2692 // ... and add it to the global vertex container
2693 global_vertices.push_back(vertex);
2694 }
2695
2696 // Is the final vertex a base vertex
2697 if (final_base_vertex_info.is_base_vertex)
2698 {
2699 // Increase the independent vertices counter
2700 n_vertices_internal_polygons++;
2701 }
2702
2703 } // for (p < n_polylines)
2704
2705 } // for (i < n_internal_polygons)
2706
2707 // The total number of added vertices in the open curves
2708 unsigned n_vertices_open_curves = 0;
2709
2710 // Do the open curves
2711 for (unsigned i = 0; i < n_open_curves; i++)
2712 {
2713 // The number of polylines in the current polygon
2714 const unsigned n_polylines = open_curves_pt[i]->ncurve_section();
2715
2716 for (unsigned p = 0; p < n_polylines; p++)
2717 {
2718 // Get a pointer to the current polyline
2719 TriangleMeshPolyLine* tmp_polyline_pt =
2720 open_curves_pt[i]->polyline_pt(p);
2721
2722 // Get the boundary id of the current polyline
2723 const unsigned bound_id = tmp_polyline_pt->boundary_id();
2724
2725 // The number of vertices in the current polyline
2726 const unsigned n_vertices = tmp_polyline_pt->nvertex();
2727
2728 // Get the current chunk number of the polyline
2729 const unsigned bnd_chunk = tmp_polyline_pt->boundary_chunk();
2730
2731 // Assign a global vertex id to the initial vertex
2732 // -----------------------------------------------
2733
2734 // Get the base vertex information of the initial vertex
2735 base_vertex_info initial_base_vertex_info =
2736 base_vertices[bound_id][bnd_chunk][0];
2737
2738#ifdef PARANOID
2739 if (!initial_base_vertex_info.has_base_vertex_assigned)
2740 {
2741 std::ostringstream error_message;
2742 std::string output_string =
2743 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2744 error_message
2745 << "The initial vertex of the current polyline has no base\n"
2746 << "vertex assigned\n"
2747 << "Open curve number: (" << i << ")\n\n"
2748 << "Polyline number: (" << p << ")\n"
2749 << "Boundary id: (" << bound_id << ")\n"
2750 << "Boundary chunk: (" << bnd_chunk << ")\n";
2751 throw OomphLibError(
2752 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2753 } // if (!initial_base_vertex_info.has_base_vertex_assigned)
2754#endif
2755
2756 // The base vertex boundary id
2757 unsigned bvbi = initial_base_vertex_info.boundary_id;
2758 // The base vertex boundary chunkx
2759 unsigned bvbc = initial_base_vertex_info.boundary_chunk;
2760 // The vertex number of the base vertex
2761 unsigned bvvn = initial_base_vertex_info.vertex_number;
2762
2763 // Get the global vertex id of the base vertex
2764 int global_vertex_id =
2765 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn];
2766
2767 // Check if the global vertex id has been already assigned
2768 if (global_vertex_id != -1)
2769 {
2770 // If that is the case then copy the global vertex id in the
2771 // current slot
2772 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][0] =
2773 global_vertex_id;
2774 } // if (global_vertex_id != -1)
2775 else
2776 {
2777 // Assign a global vertex id to the base vertex
2778 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn] =
2779 global_vertex_number;
2780
2781 // ... and copy the value to the initial vertex
2782 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][0] =
2783 global_vertex_number;
2784
2785 // ... increase the counter for the "global_vertex_number"
2786 global_vertex_number++;
2787
2788 // Get the vertex
2789 Vector<double> vertex = tmp_polyline_pt->vertex_coordinate(0);
2790
2791 // ... and add it to the global vertex container
2792 global_vertices.push_back(vertex);
2793 }
2794
2795 // Is the initial vertex a base vertex
2796 if (initial_base_vertex_info.is_base_vertex)
2797 {
2798 // Increase the independent vertices counter
2799 n_vertices_open_curves++;
2800 }
2801
2802 // ------------------------------------------------------------
2803 // Now loop over the intermediate vertices and assign a unique
2804 // vertex id (all intermediate vertices are base vertices)
2805 for (unsigned v = 1; v < n_vertices - 1; v++)
2806 {
2807 // Get the global vertex id
2808 global_vertex_id =
2809 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v];
2810
2811 // Check if the global vertex id has been already assigned.
2812 // We do nothing if it has been already assigned
2813 // (global_vertex_id!=-1).
2814
2815 // If it has not been already assigned (global_vertex_id==-1)
2816 // then set a new global vertex number, and add the vertex to
2817 // the global vertices container
2818 if (global_vertex_id == -1)
2819 {
2820 // Set a value for the global vertex
2821 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v] =
2822 global_vertex_number;
2823
2824 // ... increase the counter for the "global_vertex_number"
2825 global_vertex_number++;
2826
2827 // Add the vertex to the global vertex container
2828 Vector<double> vertex = tmp_polyline_pt->vertex_coordinate(v);
2829 // Add the vertex to the global vertex container
2830 global_vertices.push_back(vertex);
2831
2832 } // if (global_vertex_id == -1)
2833
2834 // Increase the independent vertices counter
2835 n_vertices_open_curves++;
2836
2837 } // for (n_vertices-1)
2838
2839 // Assign a global vertex id to the final vertex
2840 // -----------------------------------------------
2841
2842 // Get the base vertex information of the final vertex
2843 base_vertex_info final_base_vertex_info =
2844 base_vertices[bound_id][bnd_chunk][1];
2845
2846#ifdef PARANOID
2847 if (!final_base_vertex_info.has_base_vertex_assigned)
2848 {
2849 std::ostringstream error_message;
2850 std::string output_string =
2851 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2852 error_message
2853 << "The final vertex of the current polyline has no base\n"
2854 << "vertex assigned\n"
2855 << "Open curve number: (" << i << ")\n\n"
2856 << "Polyline number: (" << p << ")\n"
2857 << "Boundary id: (" << bound_id << ")\n"
2858 << "Boundary chunk: (" << bnd_chunk << ")\n";
2859 throw OomphLibError(
2860 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2861 } // if (!final_base_vertex_info.has_base_vertex_assigned)
2862#endif
2863
2864 // The base vertex boundary id
2865 bvbi = final_base_vertex_info.boundary_id;
2866 // The base vertex boundary chunkx
2867 bvbc = final_base_vertex_info.boundary_chunk;
2868 // The vertex number of the base vertex
2869 bvvn = final_base_vertex_info.vertex_number;
2870
2871 // Get the global vertex id of the base vertex
2872 global_vertex_id =
2873 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn];
2874
2875 // Check if the global vertex id has been already assigned
2876 if (global_vertex_id != -1)
2877 {
2878 // If that is the case then copy the global vertex id in the
2879 // current slot
2880 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
2881 [n_vertices - 1] =
2882 global_vertex_id;
2883 } // if (global_vertex_id != -1)
2884 else
2885 {
2886 // Assign a global vertex id to the base vertex
2887 boundary_chunk_vertex_to_global_vertex_id[bvbi][bvbc][bvvn] =
2888 global_vertex_number;
2889
2890 // ... and copy the value to the final vertex
2891 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
2892 [n_vertices - 1] =
2893 global_vertex_number;
2894
2895 // ... increase the counter for the "global_vertex_number"
2896 global_vertex_number++;
2897
2898 // Get the vertex
2899 Vector<double> vertex =
2900 tmp_polyline_pt->vertex_coordinate(n_vertices - 1);
2901
2902 // ... and add it to the global vertex container
2903 global_vertices.push_back(vertex);
2904 }
2905
2906 // Is the final vertex a base vertex
2907 if (final_base_vertex_info.is_base_vertex)
2908 {
2909 // Increase the independent vertices counter
2910 n_vertices_open_curves++;
2911 }
2912
2913 } // for (p < n_polylines)
2914
2915 } // for (i < n_open_curves)
2916
2917 // We have already assigned a unique id for each vertex in the
2918 // boundaries, and we have collected all the vertices in a container
2919
2920 // Store the global number of vertices. Add the number of vertices
2921 // of all the polygons (outer and internal), and the open curves to
2922 // get the total number of vertices. NB This is the number of NON
2923 // REPEATED vertices, the sum of all the vertices of each individual
2924 // polyline can be computed from the boundary_chunk_nvertices
2925 // container
2926 const unsigned n_global_vertices = n_vertices_outer_polygons +
2927 n_vertices_internal_polygons +
2928 n_vertices_open_curves;
2929
2930#ifdef PARANOID
2931 // Check that the total number of vertices be equal to the
2932 // pre-computed number of vertices
2933 const unsigned added_global_vertices = global_vertices.size();
2934 if (added_global_vertices != n_global_vertices)
2935 {
2936 std::ostringstream error_stream;
2937 std::string output_string =
2938 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2939 error_stream
2940 << "The number of added vertices to the global vertices container\n"
2941 << "is different from the pre-computed number of vertices\n"
2942 << "Added number of vertices: (" << added_global_vertices << ")\n"
2943 << "Pre-computed number of global vertices: (" << n_global_vertices
2944 << ")\n";
2945 throw OomphLibError(
2946 error_stream.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2947 } // if (added_global_vertices != n_global_vertices)
2948
2949 // Check that the counter for the global number of vertices is the same as
2950 // the pre-computed number of vertices
2951 if (global_vertex_number != n_global_vertices)
2952 {
2953 std::ostringstream error_stream;
2954 std::string output_string =
2955 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
2956 error_stream
2957 << "The counter for the global number of vertices is different from\n"
2958 << "the pre-computed number of vertices\n"
2959 << "Global vertices counter: (" << global_vertex_number << ")\n"
2960 << "Pre-computed number of global vertices: (" << n_global_vertices
2961 << ")\n";
2962 throw OomphLibError(
2963 error_stream.str(), output_string, OOMPH_EXCEPTION_LOCATION);
2964 } // if (global_vertex_number != n_global_vertices)
2965#endif
2966
2967
2968 // ------------------------------------------------------------------
2969
2970 // ------------------------------------------------------------------
2971 // 4th- stage
2972 // Create the segments using the unique id of the vertices and
2973 // assign a segment id to each segment (the one from the boundary)
2974
2975 // Create the segment storage, each segment is composed of two
2976 // vertices (the vertices id is obtained from the global vertex ids
2977 // container)
2978 Vector<Vector<int>> global_segments;
2979
2980 // Assign a segment marked to each segment (the boundary id to which
2981 // the segment belongs)
2982 Vector<int> global_segment_markers;
2983
2984 // Loop over the boundaries again, but know get the global vertex id
2985 // from the container and create the segments
2986
2987 // Start with the outer polygons
2988 for (unsigned i = 0; i < n_outer_polygons; i++)
2989 {
2990 // The number of polylines in the current polygon
2991 const unsigned n_polylines = outer_polygons_pt[i]->npolyline();
2992
2993 // Loop over the polylines
2994 for (unsigned p = 0; p < n_polylines; p++)
2995 {
2996 // Get the boundary id of the current polyline
2997 const unsigned bound_id =
2998 outer_polygons_pt[i]->polyline_pt(p)->boundary_id();
2999
3000 // The number of vertices in the current polyline
3001 const unsigned n_vertices =
3002 outer_polygons_pt[i]->polyline_pt(p)->nvertex();
3003
3004 // Get the current chunk number of the polyline
3005 const unsigned bnd_chunk =
3006 outer_polygons_pt[i]->polyline_pt(p)->boundary_chunk();
3007
3008 // Loop over the vertices in the polyline
3009 for (unsigned v = 1; v < n_vertices; v++)
3010 {
3011 // Get the global vertex id
3012 const int global_vertex_id_v_1 =
3013 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
3014 [v - 1];
3015
3016#ifdef PARANOID
3017 if (global_vertex_id_v_1 == -1)
3018 {
3019 // Get the current vertex
3020 Vector<double> current_vertex =
3021 outer_polygons_pt[i]->polyline_pt(p)->vertex_coordinate(v - 1);
3022 std::ostringstream error_message;
3023 std::string output_string =
3024 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
3025 error_message
3026 << "The global vertex number for the current vertex has not\n"
3027 << "been assigned\n."
3028 << "Outer polygon number: (" << i << ")\n\n"
3029 << "Polyline number: (" << p << ")\n"
3030 << "Boundary id: (" << bound_id << ")\n"
3031 << "Boundary chunk: (" << bnd_chunk << ")\n"
3032 << "Vertex[" << v - 1 << "]: (" << current_vertex[0] << ", "
3033 << current_vertex[1] << ")\n";
3034 throw OomphLibError(
3035 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3036 } // if (global_vertex_id_v_1 == -1)
3037#endif
3038
3039 // Get the global vertex id
3040 const int global_vertex_id_v =
3041 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v];
3042
3043#ifdef PARANOID
3044 if (global_vertex_id_v == -1)
3045 {
3046 // Get the current vertex
3047 Vector<double> current_vertex =
3048 outer_polygons_pt[i]->polyline_pt(p)->vertex_coordinate(v);
3049 std::ostringstream error_message;
3050 std::string output_string =
3051 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
3052 error_message
3053 << "The global vertex number for the current vertex has not\n"
3054 << "been assigned\n."
3055 << "Outer polygon number: (" << i << ")\n\n"
3056 << "Polyline number: (" << p << ")\n"
3057 << "Boundary id: (" << bound_id << ")\n"
3058 << "Boundary chunk: (" << bnd_chunk << ")\n"
3059 << "Vertex[" << v << "]: (" << current_vertex[0] << ", "
3060 << current_vertex[1] << ")\n";
3061 throw OomphLibError(
3062 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3063 } // if (global_vertex_id_v == -1)
3064#endif
3065
3066 // Add the vertices to the segments container
3067 Vector<int> current_segment(2);
3068 current_segment[0] = global_vertex_id_v_1;
3069 current_segment[1] = global_vertex_id_v;
3070 global_segments.push_back(current_segment);
3071
3072 // ... and associate the segments marker
3073 global_segment_markers.push_back(bound_id);
3074
3075 } // for (v < n_vertices)
3076
3077 } // for (p < n_polylines)
3078
3079 } // for (i < n_outer_polygons)
3080
3081 // Now work the internal polygons
3082 for (unsigned i = 0; i < n_internal_polygons; i++)
3083 {
3084 // The number of polylines in the current polygon
3085 const unsigned n_polylines = internal_polygons_pt[i]->npolyline();
3086
3087 // Loop over the polylines
3088 for (unsigned p = 0; p < n_polylines; p++)
3089 {
3090 // Get the boundary id of the current polyline
3091 const unsigned bound_id =
3092 internal_polygons_pt[i]->polyline_pt(p)->boundary_id();
3093
3094 // The number of vertices in the current polyline
3095 const unsigned n_vertices =
3096 internal_polygons_pt[i]->polyline_pt(p)->nvertex();
3097
3098 // Get the current chunk number of the polyline
3099 const unsigned bnd_chunk =
3100 internal_polygons_pt[i]->polyline_pt(p)->boundary_chunk();
3101
3102 // Loop over the vertices in the polyline
3103 for (unsigned v = 1; v < n_vertices; v++)
3104 {
3105 // Get the global vertex id
3106 const int global_vertex_id_v_1 =
3107 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
3108 [v - 1];
3109
3110#ifdef PARANOID
3111 if (global_vertex_id_v_1 == -1)
3112 {
3113 // Get the current vertex
3114 Vector<double> current_vertex =
3115 internal_polygons_pt[i]->polyline_pt(p)->vertex_coordinate(v - 1);
3116 std::ostringstream error_message;
3117 std::string output_string =
3118 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
3119 error_message
3120 << "The global vertex number for the current vertex has not\n"
3121 << "been assigned\n."
3122 << "Internal polygon number: (" << i << ")\n\n"
3123 << "Polyline number: (" << p << ")\n"
3124 << "Boundary id: (" << bound_id << ")\n"
3125 << "Boundary chunk: (" << bnd_chunk << ")\n"
3126 << "Vertex[" << v - 1 << "]: (" << current_vertex[0] << ", "
3127 << current_vertex[1] << ")\n";
3128 throw OomphLibError(
3129 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3130 } // if (global_vertex_id_v_1 == -1)
3131#endif
3132
3133 // Get the global vertex id
3134 const int global_vertex_id_v =
3135 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v];
3136
3137#ifdef PARANOID
3138 if (global_vertex_id_v == -1)
3139 {
3140 // Get the current vertex
3141 Vector<double> current_vertex =
3142 internal_polygons_pt[i]->polyline_pt(p)->vertex_coordinate(v);
3143 std::ostringstream error_message;
3144 std::string output_string =
3145 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
3146 error_message
3147 << "The global vertex number for the current vertex has not\n"
3148 << "been assigned\n."
3149 << "Internal polygon number: (" << i << ")\n\n"
3150 << "Polyline number: (" << p << ")\n"
3151 << "Boundary id: (" << bound_id << ")\n"
3152 << "Boundary chunk: (" << bnd_chunk << ")\n"
3153 << "Vertex[" << v << "]: (" << current_vertex[0] << ", "
3154 << current_vertex[1] << ")\n";
3155 throw OomphLibError(
3156 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3157 } // if (global_vertex_id_v == -1)
3158#endif
3159
3160 // Add the vertices to the segments container
3161 Vector<int> current_segment(2);
3162 current_segment[0] = global_vertex_id_v_1;
3163 current_segment[1] = global_vertex_id_v;
3164 global_segments.push_back(current_segment);
3165
3166 // ... and associate the segments marker
3167 global_segment_markers.push_back(bound_id);
3168
3169 } // for (v < n_vertices)
3170
3171 } // for (p < n_polylines)
3172
3173 } // for (i < n_internal_polygons)
3174
3175 // Now do the open curves
3176 for (unsigned i = 0; i < n_open_curves; i++)
3177 {
3178 // The number of polylines in the current polygon
3179 const unsigned n_polylines = open_curves_pt[i]->ncurve_section();
3180
3181 // Loop over the polylines
3182 for (unsigned p = 0; p < n_polylines; p++)
3183 {
3184 // Get the boundary id of the current polyline
3185 const unsigned bound_id =
3186 open_curves_pt[i]->polyline_pt(p)->boundary_id();
3187
3188 // The number of vertices in the current polyline
3189 const unsigned n_vertices =
3190 open_curves_pt[i]->polyline_pt(p)->nvertex();
3191
3192 // Get the current chunk number of the polyline
3193 const unsigned bnd_chunk =
3194 open_curves_pt[i]->polyline_pt(p)->boundary_chunk();
3195
3196 // Loop over the vertices in the polyline
3197 for (unsigned v = 1; v < n_vertices; v++)
3198 {
3199 // Get the global vertex id
3200 const int global_vertex_id_v_1 =
3201 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk]
3202 [v - 1];
3203
3204#ifdef PARANOID
3205 if (global_vertex_id_v_1 == -1)
3206 {
3207 // Get the current vertex
3208 Vector<double> current_vertex =
3209 open_curves_pt[i]->polyline_pt(p)->vertex_coordinate(v - 1);
3210 std::ostringstream error_message;
3211 std::string output_string =
3212 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
3213 error_message
3214 << "The global vertex number for the current vertex has not\n"
3215 << "been assigned\n."
3216 << "Open curve number: (" << i << ")\n\n"
3217 << "Polyline number: (" << p << ")\n"
3218 << "Boundary id: (" << bound_id << ")\n"
3219 << "Boundary chunk: (" << bnd_chunk << ")\n"
3220 << "Vertex[" << v - 1 << "]: (" << current_vertex[0] << ", "
3221 << current_vertex[1] << ")\n";
3222 throw OomphLibError(
3223 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3224 } // if (global_vertex_id_v_1 == -1)
3225#endif
3226
3227 // Get the global vertex id
3228 const int global_vertex_id_v =
3229 boundary_chunk_vertex_to_global_vertex_id[bound_id][bnd_chunk][v];
3230
3231#ifdef PARANOID
3232 if (global_vertex_id_v == -1)
3233 {
3234 // Get the current vertex
3235 Vector<double> current_vertex =
3236 open_curves_pt[i]->polyline_pt(p)->vertex_coordinate(v);
3237 std::ostringstream error_message;
3238 std::string output_string =
3239 "UnstructuredTwoDMeshGeometryBase::build_triangulateio()";
3240 error_message
3241 << "The global vertex number for the current vertex has not\n"
3242 << "been assigned\n."
3243 << "Open curve number: (" << i << ")\n\n"
3244 << "Polyline number: (" << p << ")\n"
3245 << "Boundary id: (" << bound_id << ")\n"
3246 << "Boundary chunk: (" << bnd_chunk << ")\n"
3247 << "Vertex[" << v << "]: (" << current_vertex[0] << ", "
3248 << current_vertex[1] << ")\n";
3249 throw OomphLibError(
3250 error_message.str(), output_string, OOMPH_EXCEPTION_LOCATION);
3251 } // if (global_vertex_id_v == -1)
3252#endif
3253
3254 // Add the vertices to the segments container
3255 Vector<int> current_segment(2);
3256 current_segment[0] = global_vertex_id_v_1;
3257 current_segment[1] = global_vertex_id_v;
3258 global_segments.push_back(current_segment);
3259
3260 // ... and associate the segments marker
3261 global_segment_markers.push_back(bound_id);
3262
3263 } // for (v < n_vertices)
3264
3265 } // for (p < n_polylines)
3266
3267 } // for (i < n_open_curves)
3268
3269 // ------------------------------------------------------------------
3270 // 5th- stage
3271 // Fill the triangulateio containers with the collected information
3272
3273 // Initialize the triangulateIO structure
3275
3276 // Set the number of vertices
3277 triangulate_io.numberofpoints = n_global_vertices;
3278
3279 // Get the number of segments
3280 const unsigned n_global_segments = global_segments.size();
3281 // Set the number of segments
3282 triangulate_io.numberofsegments = n_global_segments;
3283
3284 // Allocate memory in the triangulateIO structure to store the values
3285 triangulate_io.pointlist =
3286 (double*)malloc(triangulate_io.numberofpoints * 2 * sizeof(double));
3287 triangulate_io.segmentlist =
3288 (int*)malloc(triangulate_io.numberofsegments * 2 * sizeof(int));
3289 triangulate_io.segmentmarkerlist =
3290 (int*)malloc(triangulate_io.numberofsegments * sizeof(int));
3291
3292 // Fill triangulateIO data structure
3293 // ---------------------------------------
3294
3295 // Fill the vertices first
3296 // An external counter
3297 unsigned ii = 0;
3298 for (unsigned i = 0; i < n_global_vertices; i++, ii += 2)
3299 {
3300 triangulate_io.pointlist[ii] = global_vertices[i][0];
3301 triangulate_io.pointlist[ii + 1] = global_vertices[i][1];
3302 } // for (i < n_global_vertices)
3303
3304 // Then fill the segments, and segments markers
3305 // Reset the external counter
3306 ii = 0;
3307 for (unsigned i = 0; i < n_global_segments; i++, ii += 2)
3308 {
3309 // The segment marker should start in 1 (our enumeration started
3310 // in 0, therefore we add one to every entry)
3311 triangulate_io.segmentlist[ii] = global_segments[i][0] + 1;
3312 triangulate_io.segmentlist[ii + 1] = global_segments[i][1] + 1;
3313 // Set the segment marker as the boundary id + 1
3314 triangulate_io.segmentmarkerlist[i] = global_segment_markers[i] + 1;
3315 }
3316
3317 // Store any regions information
3318 // ---------------------------------------
3319
3320 // Get the number of regions
3321 unsigned n_regions = regions_coordinates.size();
3322
3323 // Check if there are regions to include
3324 if (n_regions > 0)
3325 {
3326 triangulate_io.numberofregions = n_regions;
3327 triangulate_io.regionlist =
3328 (double*)malloc(triangulate_io.numberofregions * 4 * sizeof(double));
3329
3330 // Loop over the regions map
3331 unsigned p = 1;
3332 for (std::map<unsigned, Vector<double>>::iterator it_regions =
3333 regions_coordinates.begin();
3334 it_regions != regions_coordinates.end();
3335 it_regions++)
3336 {
3337 // Get the region id
3338 unsigned region_id = (*it_regions).first;
3339 // Set the x-coordinate
3340 triangulate_io.regionlist[4 * p - 4] = ((*it_regions).second)[0];
3341 // Set the y-coordinate
3342 triangulate_io.regionlist[4 * p - 3] = ((*it_regions).second)[1];
3343 // Set the region id
3344 triangulate_io.regionlist[4 * p - 2] = static_cast<double>(region_id);
3345 // This is used to define a target area for the region, zero
3346 // means no target area defined
3347 triangulate_io.regionlist[4 * p - 1] = regions_areas[region_id];
3348 // Increase the auxiliary counter
3349 p++;
3350 } // Loop over the regions map
3351
3352 } // if (n_regions > 0)
3353
3354 // Holes information
3355 // ---------------------------------------
3356
3357 // Get the number of any extra defined holes
3358 const unsigned n_extra_holes = extra_holes_coordinates.size();
3359 // The internal polygon marked as a hole
3360 Vector<unsigned> internal_polygon_marked_as_hole;
3361 // Count how many internal polygons are holes
3362 for (unsigned h = 0; h < n_internal_polygons; h++)
3363 {
3364 if (!internal_polygons_pt[h]->internal_point().empty())
3365 {
3366 internal_polygon_marked_as_hole.push_back(h);
3367 } // Is internal polygon a hole?
3368 } // for (h < n_internal_polygons)
3369
3370 // Get the number of internal polygons that should be treated as
3371 // holes
3372 const unsigned n_internal_polygons_are_holes =
3373 internal_polygon_marked_as_hole.size();
3374
3375 // Set the number of holes in the triangulateIO structure
3376 triangulate_io.numberofholes =
3377 n_extra_holes + n_internal_polygons_are_holes;
3378
3379 // Allocate memory for the holes coordinates
3380 triangulate_io.holelist =
3381 (double*)malloc(triangulate_io.numberofholes * 2 * sizeof(double));
3382
3383 // Store the holes coordinates
3384 unsigned count_hole = 0;
3385 // Add first the internal polygons marked as holes
3386 for (; count_hole < n_internal_polygons_are_holes * 2; count_hole += 2)
3387 {
3388 const unsigned index_intenal_polygon_is_hole =
3389 internal_polygon_marked_as_hole[count_hole / 2];
3390 triangulate_io.holelist[count_hole] =
3391 internal_polygons_pt[index_intenal_polygon_is_hole]
3392 ->internal_point()[0];
3393 triangulate_io.holelist[count_hole + 1] =
3394 internal_polygons_pt[index_intenal_polygon_is_hole]
3395 ->internal_point()[1];
3396 } // for (count_hole < n_internal_polygons_are_holes*2)
3397
3398 // Add the extra holes coordinates
3399 for (unsigned i_eh = 0;
3400 count_hole < 2 * (n_extra_holes + n_internal_polygons_are_holes);
3401 count_hole += 2, i_eh++)
3402 {
3403 triangulate_io.holelist[count_hole] = extra_holes_coordinates[i_eh][0];
3404 triangulate_io.holelist[count_hole + 1] =
3405 extra_holes_coordinates[i_eh][1];
3406 } // for (count_hole < 2*(n_extra_holes + n_internal_polygons_are_holes))
3407
3408 // ------------------------------------------------------------------
3409 }
3410
3411 //========================================================================
3412 /// Helps to add information to the connection matrix of the
3413 /// given polyline
3414 //========================================================================
3416 TriangleMeshPolyLine* polyline_pt,
3417 std::map<unsigned, std::map<unsigned, Vector<vertex_connection_info>>>&
3418 connection_matrix,
3419 TriangleMeshPolyLine* next_polyline_pt)
3420 {
3421 // Get the boundary id of the current polyline
3422 const unsigned bound_id = polyline_pt->boundary_id();
3423
3424 // Get the chunk number associated with this boundary
3425 const unsigned bound_chunk = polyline_pt->boundary_chunk();
3426
3427 // Check if the ends of the polyline are connected
3428 const bool connected_init_end = polyline_pt->is_initial_vertex_connected();
3429
3430 // ... check for both connections
3431 const bool connected_final_end = polyline_pt->is_final_vertex_connected();
3432
3433 // Prepare the connection matrix (resize the vector to store the
3434 // initial and final vertices info.)
3435 connection_matrix[bound_id][bound_chunk].resize(2);
3436
3437 // The default information for the matrix
3438 vertex_connection_info default_vertex_connection_info;
3439 default_vertex_connection_info.is_connected = false;
3440 default_vertex_connection_info.boundary_id_to_connect = 0;
3441 default_vertex_connection_info.boundary_chunk_to_connect = 0;
3442 default_vertex_connection_info.vertex_number_to_connect = 0;
3443
3444 // Set the default information for the initial vertex
3445 connection_matrix[bound_id][bound_chunk][0] =
3446 default_vertex_connection_info;
3447 // Set the default information for the final vertex
3448 connection_matrix[bound_id][bound_chunk][1] =
3449 default_vertex_connection_info;
3450
3451 // Set the info. for the connection matrix
3452 if (connected_init_end)
3453 {
3454 // The boundary id to connect
3455 const unsigned bc = polyline_pt->initial_vertex_connected_bnd_id();
3456
3457 // The vertex number to which is connected
3458 const unsigned vc = polyline_pt->initial_vertex_connected_n_vertex();
3459
3460 // The boundary chunk to which is connected
3461 const unsigned cc = polyline_pt->initial_vertex_connected_n_chunk();
3462
3463 // Set the connection info
3464 vertex_connection_info initial_vertex_info;
3465 // Set as connected
3466 initial_vertex_info.is_connected = true;
3467 // The boundary id to connect
3468 initial_vertex_info.boundary_id_to_connect = bc;
3469 // The chunk number to connect
3470 initial_vertex_info.boundary_chunk_to_connect = cc;
3471 // The vertex number to connect
3472 initial_vertex_info.vertex_number_to_connect = vc;
3473
3474 // Add the connection information to the connection matrix
3475 connection_matrix[bound_id][bound_chunk][0] = initial_vertex_info;
3476
3477 } // if (connected_init_end)
3478
3479 if (connected_final_end)
3480 {
3481 // The boundary id to connect
3482 const unsigned bc = polyline_pt->final_vertex_connected_bnd_id();
3483
3484 // The vertex number to which is connected
3485 const unsigned vc = polyline_pt->final_vertex_connected_n_vertex();
3486
3487 // The boundary chunk to which is connected
3488 const unsigned cc = polyline_pt->final_vertex_connected_n_chunk();
3489
3490 // Set the connection info
3491 vertex_connection_info final_vertex_info;
3492 // Set as connected
3493 final_vertex_info.is_connected = true;
3494 // The boundary id to connect
3495 final_vertex_info.boundary_id_to_connect = bc;
3496 // The chunk number to connect
3497 final_vertex_info.boundary_chunk_to_connect = cc;
3498 // The vertex number to connect
3499 final_vertex_info.vertex_number_to_connect = vc;
3500
3501 // Add the connection information to the connection matrix
3502 connection_matrix[bound_id][bound_chunk][1] = final_vertex_info;
3503
3504 } // if (connected_final_end)
3505 else
3506 {
3507 // If the vertex is not connected but there is a next polyline in
3508 // the polygon (this only applies for polygons, for open curves
3509 // the next polyline is set to null)
3510
3511 if (next_polyline_pt != 0)
3512 {
3513 // Get the info of the next polyline
3514 const unsigned next_bound_id = next_polyline_pt->boundary_id();
3515
3516 // Get the chunk number associated with this boundary
3517 const unsigned next_bound_chunk = next_polyline_pt->boundary_chunk();
3518
3519 // Set the connection info
3520 vertex_connection_info final_vertex_info;
3521 // Set as connected
3522 final_vertex_info.is_connected = true;
3523 // The boundary id to connect
3524 final_vertex_info.boundary_id_to_connect = next_bound_id;
3525 // The chunk number to connect
3526 final_vertex_info.boundary_chunk_to_connect = next_bound_chunk;
3527 // The vertex number to connect
3528 final_vertex_info.vertex_number_to_connect = 0;
3529
3530 // Add the connection information to the connection matrix
3531 connection_matrix[bound_id][bound_chunk][1] = final_vertex_info;
3532
3533 } // if (next_polyline_pt != 0)
3534
3535 } // else if (connected_final_end)
3536 }
3537
3538 //========================================================================
3539 // Initialize the base vertex structure, set every vertex to
3540 // no visited and not being a base vertex
3541 //========================================================================
3543 TriangleMeshPolyLine* polyline_pt,
3544 std::map<unsigned, std::map<unsigned, Vector<base_vertex_info>>>&
3545 base_vertices)
3546 {
3547 // Get the boundary id of the current polyline
3548 const unsigned bound_id = polyline_pt->boundary_id();
3549
3550 // Get the chunk number associated with this boundary
3551 const unsigned bound_chunk = polyline_pt->boundary_chunk();
3552
3553 // Create the default data structure
3554 base_vertex_info default_base_vertex_info;
3555 // Set as not done
3556 default_base_vertex_info.has_base_vertex_assigned = false;
3557 // Set as not base vertex
3558 default_base_vertex_info.is_base_vertex = false;
3559 // Set the default base boundary id
3560 default_base_vertex_info.boundary_id = 0;
3561 // Set the default base boundary chunk
3562 default_base_vertex_info.boundary_chunk = 0;
3563 // Set the default base vertex number
3564 default_base_vertex_info.vertex_number = 0;
3565
3566 // Initialize the base vertex info. for the current polyline
3567 // Allocate memory for the initial and final vertex
3568 base_vertices[bound_id][bound_chunk].resize(2);
3569 // Set the initial vertex info.
3570 base_vertices[bound_id][bound_chunk][0] = default_base_vertex_info;
3571 // Set the final vertex info.
3572 base_vertices[bound_id][bound_chunk][1] = default_base_vertex_info;
3573 }
3574
3575 //========================================================================
3576 // Helps to identify the base vertex of the given polyline
3577 //========================================================================
3579 TriangleMeshPolyLine* polyline_pt,
3580 std::map<unsigned, std::map<unsigned, Vector<base_vertex_info>>>&
3581 base_vertices,
3582 std::map<unsigned, std::map<unsigned, Vector<vertex_connection_info>>>&
3583 connection_matrix,
3584 std::map<unsigned, std::map<unsigned, unsigned>>& boundary_chunk_nvertices)
3585 {
3586 // Get the general data of the polyline
3587
3588 // Get the boundary id of the current polyline
3589 const unsigned bound_id = polyline_pt->boundary_id();
3590
3591 // The number of vertices in the current polyline
3592 const unsigned n_vertices = polyline_pt->nvertex();
3593
3594 // Get the chunk number associated with this boundary
3595 const unsigned bound_chunk = polyline_pt->boundary_chunk();
3596
3597 // Keep track of the done vertices
3598 std::map<Vector<unsigned>, bool> done;
3599
3600 // Loop over the vertices in the polyline
3601 for (unsigned v = 0; v < n_vertices; v++)
3602 {
3603 // For each vertex find its base vertex
3604
3605 // Flag to state if repeat the search with the new vertex
3606 // reference
3607 bool repeat = false;
3608
3609 // Store the info. of the vertex currently serching for base
3610 // vertex
3611 unsigned ibnd_id = bound_id;
3612 unsigned ibnd_chunk = bound_chunk;
3613 unsigned ivertex_number = v;
3614
3615 // Store the info. of the base vertex of the current vertex
3616 unsigned base_bnd_id = 0;
3617 unsigned base_bnd_chunk = 0;
3618 unsigned base_vertex_number = 0;
3619
3620 // Store all the found references to the vertex
3621 Vector<unsigned> iter_bnd_id;
3622 Vector<unsigned> iter_bnd_chunk;
3623 Vector<unsigned> iter_vertex_number;
3624
3625 do
3626 {
3627 // Reset the flag
3628 repeat = false;
3629
3630 // Get the number of vertices of the polyline where the current
3631 // reference index lives on
3632 const unsigned i_n_vertices =
3633 boundary_chunk_nvertices[ibnd_id][ibnd_chunk];
3634 // Is the vertex an initial or final vertex on the (ibnd_id,
3635 // ibnd_chunk)
3636 bool is_initial = false;
3637 if (ivertex_number == 0)
3638 {
3639 is_initial = true;
3640 }
3641 bool is_final = false;
3642 if (ivertex_number == i_n_vertices - 1)
3643 {
3644 is_final = true;
3645 }
3646
3647 // Is initial or final?
3648 if (is_initial || is_final)
3649 {
3650 // Stores the base vertex info. of the current vertex
3651 base_vertex_info current_vertex_base_info;
3652
3653 if (is_initial)
3654 {
3655 // Get the base vertex info. of the current vertex
3656 current_vertex_base_info = base_vertices[ibnd_id][ibnd_chunk][0];
3657 }
3658 else if (is_final)
3659 {
3660 // Get the base vertex info. of the current vertex
3661 current_vertex_base_info = base_vertices[ibnd_id][ibnd_chunk][1];
3662 }
3663
3664 // Has the vertex assigned a base vertex?
3665 if (!current_vertex_base_info.has_base_vertex_assigned)
3666 {
3667 // Is the current vertex a base vertex?
3668 if (!current_vertex_base_info.is_base_vertex)
3669 {
3670 // Get the connection info. of the vertex
3671 vertex_connection_info vertex_connect_info =
3672 connection_matrix[ibnd_id][ibnd_chunk][0];
3673
3674 if (is_initial)
3675 {
3676 vertex_connect_info = connection_matrix[ibnd_id][ibnd_chunk][0];
3677 }
3678 else if (is_final)
3679 {
3680 vertex_connect_info = connection_matrix[ibnd_id][ibnd_chunk][1];
3681 }
3682
3683 // Get the complete connection information of the vertex
3684
3685 // The new connection info.
3686 const bool current_is_connected =
3687 vertex_connect_info.is_connected;
3688
3689 // Is the current vertex connected to other vertex
3690 if (current_is_connected)
3691 {
3692 // The new boundary id to connect
3693 const unsigned iibnd_id =
3694 vertex_connect_info.boundary_id_to_connect;
3695
3696 // The new boundary chunk to connect
3697 const unsigned iibnd_chunk =
3698 vertex_connect_info.boundary_chunk_to_connect;
3699
3700 // The new vertex number to connect
3701 const unsigned iivertex_number =
3702 vertex_connect_info.vertex_number_to_connect;
3703
3704 // Get the number of vertices in the new boundary to connect
3705 const unsigned ii_n_vertices =
3706 boundary_chunk_nvertices[iibnd_id][iibnd_chunk];
3707
3708 // Is the new vertex an initial or final vertex on the
3709 // (iibnd_id, iibnd_chunk)
3710 bool new_is_initial = false;
3711 if (iivertex_number == 0)
3712 {
3713 new_is_initial = true;
3714 }
3715 bool new_is_final = false;
3716 if (iivertex_number == ii_n_vertices - 1)
3717 {
3718 new_is_final = true;
3719 }
3720
3721 // Is new initial or final?
3722 if (new_is_initial || new_is_final)
3723 {
3724 // Stores the base vertex info. of the new vertex
3725 base_vertex_info new_vertex_base_info;
3726
3727 if (new_is_initial)
3728 {
3729 // Get the base vertex info. of the current vertex
3730 new_vertex_base_info =
3731 base_vertices[iibnd_id][iibnd_chunk][0];
3732 }
3733 else if (new_is_final)
3734 {
3735 // Get the base vertex info. of the current vertex
3736 new_vertex_base_info =
3737 base_vertices[iibnd_id][iibnd_chunk][1];
3738 }
3739
3740 // Verify if the new vertex has been done
3741 Vector<unsigned> check_done(3);
3742 check_done[0] = iibnd_id;
3743 check_done[1] = iibnd_chunk;
3744 check_done[2] = iivertex_number;
3745 // Has the new vertex been already visited?
3746 if (!done[check_done])
3747 {
3748 // Add the i-reference to the iteration vectors
3749 // since it needs to be assigned the base node when
3750 // it be found
3751 iter_bnd_id.push_back(ibnd_id);
3752 iter_bnd_chunk.push_back(ibnd_chunk);
3753 iter_vertex_number.push_back(ivertex_number);
3754
3755 Vector<unsigned> current_done(3);
3756 current_done[0] = ibnd_id;
3757 current_done[1] = ibnd_chunk;
3758 current_done[2] = ivertex_number;
3759
3760 // Mark the i-th reference of the vertex as done
3761 done[current_done] = true;
3762
3763 // Change the vertex reference and repeat
3764 ibnd_id = iibnd_id;
3765 ibnd_chunk = iibnd_chunk;
3766 ivertex_number = iivertex_number;
3767
3768 // Set the flag to repeat the search with the new
3769 // reference of the vertex
3770 repeat = true;
3771
3772 } // if (!done[check_done])
3773 else
3774 {
3775 // We have found a base vertex, we only need to
3776 // decide if it is the current vertex or the new
3777 // vertex
3778
3779 // Add the i-reference to the iteration vectors to
3780 // be assigned the found base vertex
3781 iter_bnd_id.push_back(ibnd_id);
3782 iter_bnd_chunk.push_back(ibnd_chunk);
3783 iter_vertex_number.push_back(ivertex_number);
3784
3785 // Has the new vertex a base vertes assigned
3786 if (!new_vertex_base_info.has_base_vertex_assigned)
3787 {
3788 // The current vertex is a base vertex (we are
3789 // looping over the current and new vertex)
3790
3791 Vector<unsigned> current_done(3);
3792 current_done[0] = ibnd_id;
3793 current_done[1] = ibnd_chunk;
3794 current_done[2] = ivertex_number;
3795
3796 // Mark the i-th reference of the vertex as done
3797 done[current_done] = true;
3798
3799 // Mark the i-th reference of the vertex as base
3800 // vertex
3801 if (is_initial)
3802 {
3803 // Mark the vertex as base vertex
3804 base_vertices[ibnd_id][ibnd_chunk][0].is_base_vertex =
3805 true;
3806 }
3807 else if (is_final)
3808 {
3809 // Mark the vertex as base vertex
3810 base_vertices[ibnd_id][ibnd_chunk][1].is_base_vertex =
3811 true;
3812 }
3813
3814 // Set as base vertex the current vertex info.
3815 // The base boundary id
3816 base_bnd_id = ibnd_id;
3817 // The base boundary chunk number
3818 base_bnd_chunk = ibnd_chunk;
3819 // The base vertex number
3820 base_vertex_number = ivertex_number;
3821
3822 } // if (!new_vertex_base_info.is_base_vertex)
3823 else
3824 {
3825 // The new vertex is a base vertex (we are looping
3826 // over the current and new vertex)
3827
3828 Vector<unsigned> current_done(3);
3829 current_done[0] = ibnd_id;
3830 current_done[1] = ibnd_chunk;
3831 current_done[2] = ivertex_number;
3832
3833 // Mark the i-th reference of the vertex as done
3834 done[current_done] = true;
3835
3836 // Set as base vertex the new vertex info.
3837 // The base boundary id
3838 base_bnd_id = iibnd_id;
3839 // The base boundary chunk number
3840 base_bnd_chunk = iibnd_chunk;
3841 // The base vertex number
3842 base_vertex_number = iivertex_number;
3843
3844 } // else if (!new_vertex_base_info.is_base_vertex)
3845
3846 } // else if (!new_vertex_base_info.done)
3847
3848 } // if (new_is_initial || new_is_final)
3849 else
3850 {
3851 // Add the i-reference to the iteration vectors since
3852 // it needs to be assigned the just found base vertex
3853 iter_bnd_id.push_back(ibnd_id);
3854 iter_bnd_chunk.push_back(ibnd_chunk);
3855 iter_vertex_number.push_back(ivertex_number);
3856
3857 Vector<unsigned> current_done(3);
3858 current_done[0] = ibnd_id;
3859 current_done[1] = ibnd_chunk;
3860 current_done[2] = ivertex_number;
3861
3862 // Mark the i-th reference of the vertex as done
3863 done[current_done] = true;
3864
3865 // Set as base node the new node info.
3866 // The base boundary id
3867 base_bnd_id = iibnd_id;
3868 // The base boundary chunk number
3869 base_bnd_chunk = iibnd_chunk;
3870 // The base vertex number
3871 base_vertex_number = iivertex_number;
3872 } // else if (new_is_initial || new_is_final)
3873
3874 } // if (current_is_connected)
3875 else
3876 {
3877 // The current vertex is a base vertex
3878
3879 // Add the i-reference to the iteration vectors to be
3880 // assigned the found base vertex
3881 iter_bnd_id.push_back(ibnd_id);
3882 iter_bnd_chunk.push_back(ibnd_chunk);
3883 iter_vertex_number.push_back(ivertex_number);
3884
3885 Vector<unsigned> current_done(3);
3886 current_done[0] = ibnd_id;
3887 current_done[1] = ibnd_chunk;
3888 current_done[2] = ivertex_number;
3889
3890 // Mark the i-th reference of the vertex as done
3891 done[current_done] = true;
3892
3893 // Mark the i-th reference of the vertex as base vertex
3894 if (is_initial)
3895 {
3896 // Mark the vertex as base vertex
3897 base_vertices[ibnd_id][ibnd_chunk][0].is_base_vertex = true;
3898 }
3899 else if (is_final)
3900 {
3901 // Mark the vertex as base vertex
3902 base_vertices[ibnd_id][ibnd_chunk][1].is_base_vertex = true;
3903 }
3904
3905 // Set as base vertex the current vertex info.
3906 // The base boundary id
3907 base_bnd_id = ibnd_id;
3908 // The base boundary chunk number
3909 base_bnd_chunk = ibnd_chunk;
3910 // The base vertex number
3911 base_vertex_number = ivertex_number;
3912
3913 } // else if (current_is_connected)
3914
3915 } // if (!current_vertex_base_info.is_base_vertex)
3916 else
3917 {
3918 // Copy the base vertex info. and assign it to all the
3919 // found vertex references
3920
3921 // The base boundary id
3922 base_bnd_id = ibnd_id;
3923 // The base boundary chunk number
3924 base_bnd_chunk = ibnd_chunk;
3925 // The base vertex number
3926 base_vertex_number = ivertex_number;
3927
3928 } // else if (!current_vertex_base_info.is_base_vertex)
3929
3930 } // if (!current_vertex_base_info.has_base_vertex_assigned)
3931 else
3932 {
3933 // Copy the base vertex info. and assign it to all the found
3934 // vertex references
3935
3936 // The base boundary id
3937 base_bnd_id = current_vertex_base_info.boundary_id;
3938 // The base boundary chunk number
3939 base_bnd_chunk = current_vertex_base_info.boundary_chunk;
3940 // The base vertex number
3941 base_vertex_number = current_vertex_base_info.vertex_number;
3942 } // else if (!current_vertex_base_info.has_base_vertex_assigned)
3943
3944 } // if (is_initial || is_final)
3945
3946 } while (repeat);
3947
3948 // Assign the base vertex to all the found references of the
3949 // vertex
3950
3951 // Get the number of found references
3952 const unsigned n_vertex_references = iter_bnd_id.size();
3953 // Loop over the found references and assign the base vertex
3954 for (unsigned r = 0; r < n_vertex_references; r++)
3955 {
3956 // Get the r-th reference of the vertex
3957 const unsigned rbnd_id = iter_bnd_id[r];
3958 const unsigned rbnd_chunk = iter_bnd_chunk[r];
3959 const unsigned rvertex_number = iter_vertex_number[r];
3960
3961 // Get the number of vertices in the r-th boundary r-th boundary
3962 // chunk
3963 const unsigned r_n_vertices =
3964 boundary_chunk_nvertices[rbnd_id][rbnd_chunk];
3965
3966 // Is the r-th vertex an initial or final vertex
3967 bool is_initial = false;
3968 if (rvertex_number == 0)
3969 {
3970 is_initial = true;
3971 }
3972 bool is_final = false;
3973 if (rvertex_number == r_n_vertices - 1)
3974 {
3975 is_final = true;
3976 }
3977
3978 if (is_initial)
3979 {
3980 // Set the base vertex info. of the r-th vertex reference
3981
3982 // Mark the vertex as having assigned a base vertex
3983 base_vertices[rbnd_id][rbnd_chunk][0].has_base_vertex_assigned = true;
3984 // The base boundary id
3985 base_vertices[rbnd_id][rbnd_chunk][0].boundary_id = base_bnd_id;
3986 // The base boundary chunk number
3987 base_vertices[rbnd_id][rbnd_chunk][0].boundary_chunk = base_bnd_chunk;
3988 // The base vertex number
3989 base_vertices[rbnd_id][rbnd_chunk][0].vertex_number =
3990 base_vertex_number;
3991 }
3992 else if (is_final)
3993 {
3994 // Set the base vertex info. of the r-th vertex reference
3995
3996 // Mark the vertex as having assigned a base vertex
3997 base_vertices[rbnd_id][rbnd_chunk][1].has_base_vertex_assigned = true;
3998 // The base boundary id
3999 base_vertices[rbnd_id][rbnd_chunk][1].boundary_id = base_bnd_id;
4000 // The base boundary chunk number
4001 base_vertices[rbnd_id][rbnd_chunk][1].boundary_chunk = base_bnd_chunk;
4002 // The base vertex number
4003 base_vertices[rbnd_id][rbnd_chunk][1].vertex_number =
4004 base_vertex_number;
4005 }
4006
4007 } // for (i < n_vertex_references)
4008
4009 } // for (v < n_vertices)
4010 }
4011
4012#endif
4013
4014
4015 //======================================================================
4016 /// Move the nodes on boundaries with associated Geometric Objects so
4017 /// that they exactly coincide with the geometric object. This requires
4018 /// that the boundary coordinates are set up consistently
4019 //======================================================================
4021 {
4022 // Loop over all boundaries
4023 unsigned n_bound = this->nboundary();
4024 for (unsigned b = 0; b < n_bound; b++)
4025 {
4026 // Find the geometric object
4027 GeomObject* const geom_object_pt = this->boundary_geom_object_pt(b);
4028
4029 // If there is one
4030 if (geom_object_pt != 0)
4031 {
4032 Vector<double> b_coord(1);
4033 Vector<double> new_x(2);
4034 const unsigned n_boundary_node = this->nboundary_node(b);
4035 for (unsigned n = 0; n < n_boundary_node; ++n)
4036 {
4037 // Get the boundary node and coordinates
4038 Node* const nod_pt = this->boundary_node_pt(b, n);
4039 nod_pt->get_coordinates_on_boundary(b, b_coord);
4040
4041 // Get the position and time history according to the underlying
4042 // geometric object, assuming that it has the same timestepper
4043 // as the nodes....
4044 unsigned n_tvalues =
4045 1 + nod_pt->position_time_stepper_pt()->nprev_values();
4046 for (unsigned t = 0; t < n_tvalues; ++t)
4047 {
4048 // Get the position according to the underlying geometric object
4049 geom_object_pt->position(t, b_coord, new_x);
4050
4051 // Move the node
4052 for (unsigned i = 0; i < 2; i++)
4053 {
4054 nod_pt->x(t, i) = new_x[i];
4055 }
4056 }
4057 }
4058 }
4059 }
4060 }
4061
4062 //======================================================================
4063 /// Helper function that checks if a given point is inside a polygon
4064 //======================================================================
4066 Vector<Vector<double>> polygon_vertices, Vector<double> point)
4067 {
4068 // Total number of vertices (the first and last vertex should be the same)
4069 const unsigned n_vertex = polygon_vertices.size();
4070
4071 // Check if internal point is actually located in bounding polygon
4072 // Reference: http://paulbourke.net/geometry/insidepoly/
4073
4074 // Counter for number of intersections
4075 unsigned intersect_counter = 0;
4076
4077 // Get first vertex
4078 Vector<double> p1 = polygon_vertices[0];
4079 for (unsigned i = 1; i <= n_vertex; i++)
4080 {
4081 // Get second vertex by wrap-around
4082 Vector<double> p2 = polygon_vertices[i % n_vertex];
4083
4084 if (point[1] > std::min(p1[1], p2[1]))
4085 {
4086 if (point[1] <= std::max(p1[1], p2[1]))
4087 {
4088 if (point[0] <= std::max(p1[0], p2[0]))
4089 {
4090 if (p1[1] != p2[1])
4091 {
4092 double xintersect =
4093 (point[1] - p1[1]) * (p2[0] - p1[0]) / (p2[1] - p1[1]) + p1[0];
4094 if ((p1[0] == p2[0]) || (point[0] <= xintersect))
4095 {
4096 intersect_counter++;
4097 }
4098 }
4099 }
4100 }
4101 }
4102 p1 = p2;
4103 }
4104
4105 // Even number of intersections: outside
4106 if (intersect_counter % 2 == 0)
4107 {
4108 return false;
4109 }
4110 else
4111 {
4112 return true;
4113 }
4114 }
4115
4116 //======================================================================
4117 /// Gets the vertex number on the destination polyline (used /
4118 /// to create the connections among shared boundaries)
4119 //======================================================================
4122 TriangleMeshPolyLine* dst_polyline_pt,
4123 Vector<double>& vertex_coordinates,
4124 unsigned& vertex_number)
4125 {
4126 // The returning flag indicating that the vertex was found in the
4127 // destination boundary
4128 bool found_vertex_number = false;
4129
4130 // Get the number of vertices in the destination polyline
4131 const unsigned n_vertices = dst_polyline_pt->nvertex();
4132
4133 // Loop over the vertices and return the closest vertex
4134 // to the given vertex coordinates
4135 for (unsigned i = 0; i < n_vertices; i++)
4136 {
4137 // Get the i-vertex on the polyline
4138 Vector<double> current_vertex = dst_polyline_pt->vertex_coordinate(i);
4139
4140 // Compute the distance to the input vertex
4141 double dist = (vertex_coordinates[0] - current_vertex[0]) *
4142 (vertex_coordinates[0] - current_vertex[0]) +
4143 (vertex_coordinates[1] - current_vertex[1]) *
4144 (vertex_coordinates[1] - current_vertex[1]);
4145 dist = sqrt(dist);
4146
4147 // Have we found the vertex?
4149 {
4150 // Set the vertex index
4151 vertex_number = i;
4152
4153 // Set the flag to found
4154 found_vertex_number = true;
4155
4156 // Break the serching
4157 break;
4158 }
4159 } // for (i < n_vertices)
4160
4161 // Return if the connection was found
4162 return found_vertex_number;
4163 }
4164
4165
4166#ifdef OOMPH_HAS_TRIANGLE_LIB
4167
4168 //======================================================================
4169 /// Helper function to copy the connection information from
4170 /// the input curve(polyline or curviline) to the output polyline
4171 //======================================================================
4173 TriangleMeshCurveSection* input_curve_pt,
4174 TriangleMeshCurveSection* output_curve_pt)
4175 {
4176 // Is there a connection to the initial vertex
4177 const bool initial_connection =
4178 input_curve_pt->is_initial_vertex_connected();
4179
4180 // Is there a connection to the initial vertex
4181 const bool final_connection = input_curve_pt->is_final_vertex_connected();
4182
4183 // If there are any connection at the ends then copy the information
4184 if (initial_connection || final_connection)
4185 {
4186 // Get the initial and final vertex from the input polyline
4187 Vector<double> input_initial_vertex(2);
4188 input_curve_pt->initial_vertex_coordinate(input_initial_vertex);
4189 Vector<double> input_final_vertex(2);
4190 input_curve_pt->final_vertex_coordinate(input_final_vertex);
4191
4192 // Get the initial and final vertex from the output polyline
4193 Vector<double> output_initial_vertex(2);
4194 output_curve_pt->initial_vertex_coordinate(output_initial_vertex);
4195 Vector<double> output_final_vertex(2);
4196 output_curve_pt->final_vertex_coordinate(output_final_vertex);
4197
4198 // Check if the polyline is in the same order or if it is
4199 // reversed. We need to know to copy the connection information
4200 // correctly
4201
4202 // The flag to state if the copy is in the same order or in
4203 // reversed order
4204 bool copy_reversed = false;
4205
4206 // Start with the initial vertices
4207 double dist_initial =
4208 ((output_initial_vertex[0] - input_initial_vertex[0]) *
4209 (output_initial_vertex[0] - input_initial_vertex[0])) +
4210 ((output_initial_vertex[1] - input_initial_vertex[1]) *
4211 (output_initial_vertex[1] - input_initial_vertex[1]));
4212 dist_initial = sqrt(dist_initial);
4213
4214 // Compute the distance of the final vertices
4215 double dist_final = ((output_final_vertex[0] - input_final_vertex[0]) *
4216 (output_final_vertex[0] - input_final_vertex[0])) +
4217 ((output_final_vertex[1] - input_final_vertex[1]) *
4218 (output_final_vertex[1] - input_final_vertex[1]));
4219 dist_final = sqrt(dist_final);
4220
4221 // Get the distace from the input vertices to the output vertices
4222 // in the same order
4223 const double dist = dist_initial + dist_final;
4224
4225 // Compute the distance between the input initial vertex and the
4226 // output final vertex
4227 double dist_initial_rev =
4228 ((input_initial_vertex[0] - output_final_vertex[0]) *
4229 (input_initial_vertex[0] - output_final_vertex[0])) +
4230 ((input_initial_vertex[1] - output_final_vertex[1]) *
4231 (input_initial_vertex[1] - output_final_vertex[1]));
4232 dist_initial_rev = sqrt(dist_initial_rev);
4233
4234 // Compute the distance between the input final vertex and the
4235 // output initial vertex
4236 double dist_final_rev =
4237 ((input_final_vertex[0] - output_initial_vertex[0]) *
4238 (input_final_vertex[0] - output_initial_vertex[0])) +
4239 ((input_final_vertex[1] - output_initial_vertex[1]) *
4240 (input_final_vertex[1] - output_initial_vertex[1]));
4241 dist_final_rev = sqrt(dist_final_rev);
4242
4243 // Get the distace from the input to the output vertices in
4244 // reversed order
4245 const double dist_rev = dist_initial_rev + dist_final_rev;
4246
4247 // If the distance is smaller than the reversed distance then the
4248 // order of the vertices is the same
4249 if (dist <= dist_rev)
4250 {
4251 // Do nothing
4252 copy_reversed = false;
4253 }
4254 else if (dist_rev < dist)
4255 {
4256 // The connection information is copied in reversed order
4257 copy_reversed = true;
4258 } // else if (dist_rev < dist)
4259
4260 // Copy the connection information
4261 // ----------------------------------------------------------------
4262 // Copy in reversed order? (copy normal)
4263 if (!copy_reversed)
4264 {
4265 // Initial vertex
4266 if (input_curve_pt->is_initial_vertex_connected())
4267 {
4268 output_curve_pt->set_initial_vertex_connected();
4269
4270 output_curve_pt->initial_vertex_connected_bnd_id() =
4271 input_curve_pt->initial_vertex_connected_bnd_id();
4272
4273 output_curve_pt->initial_vertex_connected_n_chunk() =
4274 input_curve_pt->initial_vertex_connected_n_chunk();
4275
4276 // We need to know if we have to compute the vertex number or
4277 // if we need just to copy it
4278 if (input_curve_pt->is_initial_vertex_connected_to_curviline())
4279 {
4280 double initial_s_connection =
4281 input_curve_pt->initial_s_connection_value();
4282
4283 unsigned bnd_id =
4284 output_curve_pt->initial_vertex_connected_bnd_id();
4285
4286 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4287
4288 output_curve_pt->initial_vertex_connected_n_vertex() =
4290 initial_s_connection, bnd_id, s_tolerance);
4291
4292 // The output polyline is not longer connected to a curviline
4293 // because we will be using the polyline representation of the
4294 // curviline
4296 }
4297 else
4298 {
4299 output_curve_pt->initial_vertex_connected_n_vertex() =
4300 input_curve_pt->initial_vertex_connected_n_vertex();
4301 }
4302
4303 } // if (input_curve_pt->is_initial_vertex_connected())
4304
4305 // Final vertex
4306 if (input_curve_pt->is_final_vertex_connected())
4307 {
4308 output_curve_pt->set_final_vertex_connected();
4309
4310 output_curve_pt->final_vertex_connected_bnd_id() =
4311 input_curve_pt->final_vertex_connected_bnd_id();
4312
4313 output_curve_pt->final_vertex_connected_n_chunk() =
4314 input_curve_pt->final_vertex_connected_n_chunk();
4315
4316 // We need to know if we have to compute the vertex number or
4317 // if we need just to copy it
4318 if (input_curve_pt->is_final_vertex_connected_to_curviline())
4319 {
4320 double final_s_connection =
4321 input_curve_pt->final_s_connection_value();
4322
4323 unsigned bnd_id = input_curve_pt->final_vertex_connected_bnd_id();
4324
4325 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4326
4327 output_curve_pt->final_vertex_connected_n_vertex() =
4329 final_s_connection, bnd_id, s_tolerance);
4330
4331 // The output polyline is not longer connected to a curviline
4332 // because we will be using the polyline representation of the
4333 // curviline
4335 }
4336 else
4337 {
4338 output_curve_pt->final_vertex_connected_n_vertex() =
4339 input_curve_pt->final_vertex_connected_n_vertex();
4340 }
4341
4342 } // if (input_curve_pt->is_final_vertex_connected())
4343
4344 } // if (!copy_reversed)
4345 else
4346 {
4347 // Copy the connections in reversed order
4348
4349 // Initial vertex
4350 if (input_curve_pt->is_initial_vertex_connected())
4351 {
4352 output_curve_pt->set_final_vertex_connected();
4353
4354 output_curve_pt->final_vertex_connected_bnd_id() =
4355 input_curve_pt->initial_vertex_connected_bnd_id();
4356
4357 output_curve_pt->final_vertex_connected_n_chunk() =
4358 input_curve_pt->initial_vertex_connected_n_chunk();
4359
4360 // We need to know if we have to compute the vertex number or
4361 // if we need just to copy it
4362 if (input_curve_pt->is_initial_vertex_connected_to_curviline())
4363 {
4364 double initial_s_connection =
4365 input_curve_pt->initial_s_connection_value();
4366
4367 unsigned bnd_id = input_curve_pt->initial_vertex_connected_bnd_id();
4368
4369 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4370
4371 output_curve_pt->final_vertex_connected_n_vertex() =
4373 initial_s_connection, bnd_id, s_tolerance);
4374
4375 // The output polyline is not longer connected to a curviline
4376 // because we will be using the polyline representation of the
4377 // curviline
4379 }
4380 else
4381 {
4382 output_curve_pt->final_vertex_connected_n_vertex() =
4383 input_curve_pt->initial_vertex_connected_n_vertex();
4384 }
4385
4386 } // if (input_curve_pt->is_initial_vertex_connected())
4387
4388 // Final vertex
4389 if (input_curve_pt->is_final_vertex_connected())
4390 {
4391 output_curve_pt->set_initial_vertex_connected();
4392
4393 output_curve_pt->initial_vertex_connected_bnd_id() =
4394 input_curve_pt->final_vertex_connected_bnd_id();
4395
4396 output_curve_pt->initial_vertex_connected_n_chunk() =
4397 input_curve_pt->final_vertex_connected_n_chunk();
4398
4399 // We need to know if we have to compute the vertex number or
4400 // if we need just to copy it
4401 if (input_curve_pt->is_final_vertex_connected_to_curviline())
4402 {
4403 double final_s_connection =
4404 input_curve_pt->final_s_connection_value();
4405
4406 unsigned bnd_id = input_curve_pt->final_vertex_connected_bnd_id();
4407
4408 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4409
4410 output_curve_pt->initial_vertex_connected_n_vertex() =
4412 final_s_connection, bnd_id, s_tolerance);
4413
4414 // The output polyline is not longer connected to a curviline
4415 // because we will be using the polyline representation of the
4416 // curviline
4418 }
4419 else
4420 {
4421 output_curve_pt->initial_vertex_connected_n_vertex() =
4422 input_curve_pt->final_vertex_connected_n_vertex();
4423 }
4424 } // if (input_curve_pt->is_final_vertex_connected())
4425 } // else if (!copy_reversed)
4426 } // if (initial_connection || final_connection)
4427 }
4428
4429 //======================================================================
4430 /// Helper function to copy the connection information from
4431 /// the input curve(polyline or curviline) to the output sub-polyline
4432 //======================================================================
4435 TriangleMeshCurveSection* input_curve_pt,
4436 TriangleMeshCurveSection* output_curve_pt)
4437 {
4438 // Is there a connection to the initial vertex
4439 const bool initial_connection =
4440 input_curve_pt->is_initial_vertex_connected();
4441
4442 // Is there a connection to the initial vertex
4443 const bool final_connection = input_curve_pt->is_final_vertex_connected();
4444
4445 // If there are any connection at the ends then copy the information
4446 if (initial_connection || final_connection)
4447 {
4448 // Get the initial and final vertex from the input polyline
4449 Vector<double> input_initial_vertex(2);
4450 input_curve_pt->initial_vertex_coordinate(input_initial_vertex);
4451 Vector<double> input_final_vertex(2);
4452 input_curve_pt->final_vertex_coordinate(input_final_vertex);
4453
4454 // Get the initial and final vertex from the output polyline
4455 Vector<double> output_initial_vertex(2);
4456 output_curve_pt->initial_vertex_coordinate(output_initial_vertex);
4457 Vector<double> output_final_vertex(2);
4458 output_curve_pt->final_vertex_coordinate(output_final_vertex);
4459
4460 // Check if the polyline is in the same order or if it is
4461 // reversed. We need to know to copy the connection information
4462 // correctly
4463
4464 // The flag to state if the copy is in the same order or in
4465 // reversed order
4466 bool copy_reversed = false;
4467
4468 // Start with the initial vertices
4469 double dist_initial =
4470 ((output_initial_vertex[0] - input_initial_vertex[0]) *
4471 (output_initial_vertex[0] - input_initial_vertex[0])) +
4472 ((output_initial_vertex[1] - input_initial_vertex[1]) *
4473 (output_initial_vertex[1] - input_initial_vertex[1]));
4474 dist_initial = sqrt(dist_initial);
4475
4476 // Compute the distance of the final vertices
4477 double dist_final = ((output_final_vertex[0] - input_final_vertex[0]) *
4478 (output_final_vertex[0] - input_final_vertex[0])) +
4479 ((output_final_vertex[1] - input_final_vertex[1]) *
4480 (output_final_vertex[1] - input_final_vertex[1]));
4481 dist_final = sqrt(dist_final);
4482
4483 // Get the distace from the input vertices to the output vertices
4484 // in the same order
4485 const double dist = dist_initial + dist_final;
4486
4487 // Compute the distance between the input initial vertex and the
4488 // output final vertex
4489 double dist_initial_rev =
4490 ((input_initial_vertex[0] - output_final_vertex[0]) *
4491 (input_initial_vertex[0] - output_final_vertex[0])) +
4492 ((input_initial_vertex[1] - output_final_vertex[1]) *
4493 (input_initial_vertex[1] - output_final_vertex[1]));
4494 dist_initial_rev = sqrt(dist_initial_rev);
4495
4496 // Compute the distance between the input final vertex and the
4497 // output initial vertex
4498 double dist_final_rev =
4499 ((input_final_vertex[0] - output_initial_vertex[0]) *
4500 (input_final_vertex[0] - output_initial_vertex[0])) +
4501 ((input_final_vertex[1] - output_initial_vertex[1]) *
4502 (input_final_vertex[1] - output_initial_vertex[1]));
4503 dist_final_rev = sqrt(dist_final_rev);
4504
4505 // Get the distace from the input to the output vertices in
4506 // reversed order
4507 const double dist_rev = dist_initial_rev + dist_final_rev;
4508
4509 // If the distance is smaller than the reversed distance then the
4510 // order of the vertices is the same
4511 if (dist <= dist_rev)
4512 {
4513 // Do nothing
4514 copy_reversed = false;
4515 }
4516 else if (dist_rev < dist)
4517 {
4518 // The connection information is copied in reversed order
4519 copy_reversed = true;
4520 } // else if (dist_rev < dist)
4521
4522 // Copy the connection information
4523 // ----------------------------------------------------------------
4524 // Copy in reversed order? (copy normal)
4525 if (!copy_reversed)
4526 {
4527 // Initial vertex. We need to double check that the initial
4528 // vertex of the input curve section correspond with the
4529 // initial vertex of the output sub-polyline. It may be the
4530 // case that this sub-polyline has not the initial vertex
4531 if (input_curve_pt->is_initial_vertex_connected() &&
4532 dist_initial <
4534 {
4535 output_curve_pt->set_initial_vertex_connected();
4536
4537 output_curve_pt->initial_vertex_connected_bnd_id() =
4538 input_curve_pt->initial_vertex_connected_bnd_id();
4539
4540 output_curve_pt->initial_vertex_connected_n_chunk() =
4541 input_curve_pt->initial_vertex_connected_n_chunk();
4542
4543 // We need to know if we have to compute the vertex
4544 // number or if we need just to copy it
4545 if (input_curve_pt->is_initial_vertex_connected_to_curviline())
4546 {
4547 double initial_s_connection =
4548 input_curve_pt->initial_s_connection_value();
4549
4550 unsigned bnd_id =
4551 output_curve_pt->initial_vertex_connected_bnd_id();
4552
4553 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4554
4555 output_curve_pt->initial_vertex_connected_n_vertex() =
4557 initial_s_connection, bnd_id, s_tolerance);
4558
4559 // The output polyline is not longer connected to a
4560 // curviline because we will be using the polyline
4561 // representation of the curviline
4563 }
4564 else
4565 {
4566 output_curve_pt->initial_vertex_connected_n_vertex() =
4567 input_curve_pt->initial_vertex_connected_n_vertex();
4568 }
4569 } // input_curve_pt->is_initial_vertex_connected() &&
4570 // dist_initial<Tolerance
4571
4572 // Final vertex. We need to double check that the final vertex
4573 // of the input curve section correspond with the final vertex
4574 // of the output sub-polyline. It may be the case that this
4575 // sub-polyline has not the final vertex
4576 if (input_curve_pt->is_final_vertex_connected() &&
4578 {
4579 output_curve_pt->set_final_vertex_connected();
4580
4581 output_curve_pt->final_vertex_connected_bnd_id() =
4582 input_curve_pt->final_vertex_connected_bnd_id();
4583
4584 output_curve_pt->final_vertex_connected_n_chunk() =
4585 input_curve_pt->final_vertex_connected_n_chunk();
4586
4587 // We need to know if we have to compute the vertex number or
4588 // if we need just to copy it
4589 if (input_curve_pt->is_final_vertex_connected_to_curviline())
4590 {
4591 double final_s_connection =
4592 input_curve_pt->final_s_connection_value();
4593
4594 unsigned bnd_id = input_curve_pt->final_vertex_connected_bnd_id();
4595
4596 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4597
4598 output_curve_pt->final_vertex_connected_n_vertex() =
4600 final_s_connection, bnd_id, s_tolerance);
4601
4602 // The output polyline is not longer connected to a
4603 // curviline because we will be using the polyline
4604 // representation of the curviline
4606 }
4607 else
4608 {
4609 output_curve_pt->final_vertex_connected_n_vertex() =
4610 input_curve_pt->final_vertex_connected_n_vertex();
4611 }
4612 } // input_curve_pt->is_final_vertex_connected() && dist_final<Tolerance
4613 } // if (!copy_reversed)
4614 else
4615 {
4616 // Copy the connections in reversed order
4617
4618 // Initial vertex. We need to double check that the initial
4619 // vertex of the input curve section correspond with the
4620 // initial vertex of the output sub-polyline. It may be the
4621 // case that this sub-polyline has not the initial vertex
4622 if (input_curve_pt->is_initial_vertex_connected() &&
4623 dist_initial_rev <
4625 {
4626 output_curve_pt->set_final_vertex_connected();
4627
4628 output_curve_pt->final_vertex_connected_bnd_id() =
4629 input_curve_pt->initial_vertex_connected_bnd_id();
4630
4631 output_curve_pt->final_vertex_connected_n_chunk() =
4632 input_curve_pt->initial_vertex_connected_n_chunk();
4633
4634 // We need to know if we have to compute the vertex number or
4635 // if we need just to copy it
4636 if (input_curve_pt->is_initial_vertex_connected_to_curviline())
4637 {
4638 double initial_s_connection =
4639 input_curve_pt->initial_s_connection_value();
4640
4641 unsigned bnd_id = input_curve_pt->initial_vertex_connected_bnd_id();
4642
4643 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4644
4645 output_curve_pt->final_vertex_connected_n_vertex() =
4647 initial_s_connection, bnd_id, s_tolerance);
4648
4649 // The output polyline is not longer connected to a
4650 // curviline because we will be using the polyline
4651 // representation of the curviline
4653 }
4654 else
4655 {
4656 output_curve_pt->final_vertex_connected_n_vertex() =
4657 input_curve_pt->initial_vertex_connected_n_vertex();
4658 }
4659 } // input_curve_pt->is_final_vertex_connected() && dist_final<Tolerance
4660
4661 // Final vertex. We need to double check that the final vertex
4662 // of the input curve section correspond with the final vertex
4663 // of the output sub-polyline. It may be the case that this
4664 // sub-polyline has not the final vertex
4665 if (input_curve_pt->is_final_vertex_connected() &&
4666 dist_final_rev <
4668 {
4669 output_curve_pt->set_initial_vertex_connected();
4670
4671 output_curve_pt->initial_vertex_connected_bnd_id() =
4672 input_curve_pt->final_vertex_connected_bnd_id();
4673
4674 output_curve_pt->initial_vertex_connected_n_chunk() =
4675 input_curve_pt->final_vertex_connected_n_chunk();
4676
4677 // We need to know if we have to compute the vertex number or
4678 // if we need just to copy it
4679 if (input_curve_pt->is_final_vertex_connected_to_curviline())
4680 {
4681 double final_s_connection =
4682 input_curve_pt->final_s_connection_value();
4683
4684 unsigned bnd_id = input_curve_pt->final_vertex_connected_bnd_id();
4685
4686 double s_tolerance = input_curve_pt->tolerance_for_s_connection();
4687
4688 output_curve_pt->initial_vertex_connected_n_vertex() =
4690 final_s_connection, bnd_id, s_tolerance);
4691
4692 // The output polyline is not longer connected to a
4693 // curviline because we will be using the polyline
4694 // representation of the curviline
4696 }
4697 else
4698 {
4699 output_curve_pt->initial_vertex_connected_n_vertex() =
4700 input_curve_pt->final_vertex_connected_n_vertex();
4701 }
4702 } // input_curve_pt->is_final_vertex_connected() && dist_final<Tolerance
4703 } // else if (!copy_reversed)
4704 } // if (initial_connection || final_connection)
4705 }
4706
4707
4708 /// Public static flag to suppress warning; defaults to true because
4709 /// it really clutters up the output. It's unlikely to ever be a
4710 /// genuine error...
4711 bool UnstructuredTwoDMeshGeometryBase::
4712 Suppress_warning_about_regions_and_boundaries = true;
4713
4714
4715#endif // OOMPH_HAS_TRIANGLE_LIB
4716
4717
4718} // namespace oomph
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
/////////////////////////////////////////////////////////////////////
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).
unsigned long nboundary_node(const unsigned &ibound) const
Return number of nodes on a particular boundary.
Definition: mesh.h:833
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:827
Node *& boundary_node_pt(const unsigned &b, const unsigned &n)
Return pointer to node n on boundary b.
Definition: mesh.h:493
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
TimeStepper *& position_time_stepper_pt()
Return a pointer to the position timestepper.
Definition: nodes.h:1022
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.cc:2379
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
virtual unsigned nprev_values() const =0
Number of previous values available: 0 for static, 1 for BDF<1>,...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
TriangleMeshClosedCurve(const Vector< TriangleMeshCurveSection * > &curve_section_pt, const Vector< double > &internal_point_pt=Vector< double >(0), const bool &is_internal_point_fixed=false)
Constructor prototype.
Vector< double > Internal_point_pt
Vector of vertex coordinates.
Base class for defining a triangle mesh boundary, this class has the methods that allow to connect th...
bool Initial_vertex_connected_to_curviline
States if the initial vertex is connected to a curviline.
virtual void initial_vertex_coordinate(Vector< double > &vertex)=0
Get first vertex coordinates.
bool is_initial_vertex_connected_to_curviline() const
Test whether the initial vertex is connected to a curviline.
bool Final_vertex_connected
Used for stating if the final end is connected to another boundary.
void unset_final_vertex_connected_to_curviline()
Sets the final vertex as non connected to a curviline.
unsigned Final_vertex_connected_n_vertex
Stores the vertex number used for connection with the final end.
double final_s_connection_value() const
Gets the s value to which the final end is connected.
bool Final_vertex_connected_to_curviline
States if the final vertex is connected to a curviline.
double Initial_s_connection_value
Stores the s value used for connecting the initial end with a curviline.
virtual unsigned boundary_id() const =0
Boundary id.
void connect_initial_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Connects the initial vertex of the curve section to a desired target polyline by specifying the verte...
double Final_s_connection_value
Stores the s value used for connecting the final end with a curviline.
double tolerance_for_s_connection() const
Gets the tolerance value for connections among curvilines.
unsigned Final_vertex_connected_bnd_id
Stores the id to which the initial end is connected.
unsigned Initial_vertex_connected_n_chunk
Stores the chunk number of the boundary to which is connected th initial end.
unsigned initial_vertex_connected_n_vertex() const
Gets the vertex number to which the initial end is connected.
void connect_final_vertex_to_polyline(TriangleMeshPolyLine *polyline_pt, const unsigned &vertex_number, const double &tolerance_for_connection=1.0e-14)
Connects the final vertex of the curve section to a desired target polyline by specifying the vertex ...
unsigned final_vertex_connected_n_chunk() const
Gets the boundary chunk to which the final end is connected.
bool Initial_vertex_connected
Used for stating if the initial end is connected to another boundary.
unsigned final_vertex_connected_n_vertex() const
Sets the vertex number to which the final end is connected.
unsigned Initial_vertex_connected_n_vertex
Stores the vertex number used for connection with the initial end.
bool is_final_vertex_connected_to_curviline() const
Test whether the final vertex is connected to a curviline.
bool is_final_vertex_connected() const
Test whether final vertex is connected or not.
double Tolerance_for_s_connection
Tolerance used for connecting the ends to a curviline.
unsigned initial_vertex_connected_bnd_id() const
Gets the id to which the initial end is connected.
virtual void final_vertex_coordinate(Vector< double > &vertex)=0
Get last vertex coordinates.
bool is_initial_vertex_connected() const
Test whether initial vertex is connected or not.
void set_initial_vertex_connected()
Sets the initial vertex as connected.
unsigned Final_vertex_connected_n_chunk
Stores the chunk number of the boundary to which is connected th initial end.
void set_final_vertex_connected()
Sets the final vertex as connected.
void connect_initial_vertex_to_curviline(TriangleMeshCurviLine *curviline_pt, const double &s_value, const double &tolerance_for_connection=1.0e-14)
Connects the initial vertex of the curve section to a desired target curviline by specifying the s va...
double initial_s_connection_value() const
Gets the s value to which the initial end is connected.
unsigned initial_vertex_connected_n_chunk() const
Gets the boundary chunk to which the initial end is connected.
void connect_final_vertex_to_curviline(TriangleMeshCurviLine *curviline_pt, const double &s_value, const double &tolerance_for_connection=1.0e-14)
Connects the final vertex of the curve section to a desired target curviline by specifying the s valu...
unsigned final_vertex_connected_bnd_id() const
Gets the id to which the final end is connected.
unsigned Initial_vertex_connected_bnd_id
Stores the id to which the initial end is connected.
void unset_initial_vertex_connected_to_curviline()
Sets the initial vertex as non connected to a curviline.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
virtual TriangleMeshCurveSection * curve_section_pt(const unsigned &i) const
Pointer to i-th constituent curve section.
Vector< TriangleMeshCurveSection * > Curve_section_pt
Vector of curve sections.
Class definining a curvilinear triangle mesh boundary in terms of a GeomObject. Curvlinear equivalent...
void add_connection_point(const double &z_value, const double &tol=1.0e-12)
Add the connection point (z_value) to the connection points that receive the curviline.
unsigned boundary_chunk() const
Boundary chunk (Used when a boundary is represented by more than one polyline.
double zeta_start()
Start coordinate in terms of the GeomObject's intrinsic coordinate.
double zeta_end()
End coordinate in terms of the GeomObject's intrinsic coordinate.
GeomObject * geom_object_pt()
Pointer to GeomObject that represents this part of the boundary.
TriangleMeshOpenCurve(const Vector< TriangleMeshCurveSection * > &curve_section_pt)
Constructor.
Class defining a polyline for use in Triangle Mesh generation.
Vector< double > vertex_coordinate(const unsigned &i) const
Coordinate vector of i-th vertex (const version)
unsigned boundary_chunk() const
Boundary chunk (Used when a boundary is represented by more than one polyline.
unsigned nvertex() const
Number of vertices.
unsigned npolyline() const
Number of constituent polylines.
TriangleMeshPolyLine * polyline_pt(const unsigned &i) const
Pointer to i-th constituent polyline.
TriangleMeshPolygon(const Vector< TriangleMeshCurveSection * > &boundary_polyline_pt, const Vector< double > &internal_point_pt=Vector< double >(0), const bool &is_internal_point_fixed=false)
Constructor: Specify vector of pointers to TriangleMeshCurveSection that define the boundary of the s...
const bool get_connected_vertex_number_on_destination_polyline(TriangleMeshPolyLine *dst_polyline_pt, Vector< double > &vertex_coordinates, unsigned &vertex_number)
Gets the vertex number on the destination polyline (used to create the connections among shared bound...
void snap_nodes_onto_geometric_objects()
Snap the boundary nodes onto any curvilinear boundaries defined by geometric objects.
unsigned get_associated_vertex_to_svalue(double &target_s_value, unsigned &bnd_id)
Get the associated vertex to the given s value by looking to the list of s values created when changi...
void copy_connection_information_to_sub_polylines(TriangleMeshCurveSection *input_curve_pt, TriangleMeshCurveSection *output_curve_pt)
Helper function to copy the connection information from the input curve(polyline or curviline) to the...
void copy_connection_information(TriangleMeshCurveSection *input_curve_pt, TriangleMeshCurveSection *output_curve_pt)
Helper function to copy the connection information from the input curve(polyline or curviline) to the...
void build_triangulateio(Vector< TriangleMeshPolygon * > &outer_polygons_pt, Vector< TriangleMeshPolygon * > &internal_polygons_pt, Vector< TriangleMeshOpenCurve * > &open_curves_pt, Vector< Vector< double > > &extra_holes_coordinates, std::map< unsigned, Vector< double > > &regions_coordinates, std::map< unsigned, double > &regions_areas, TriangulateIO &triangulate_io)
Create TriangulateIO object from outer boundaries, internal boundaries, and open curves....
bool is_point_inside_polygon_helper(Vector< Vector< double > > polygon_vertices, Vector< double > point)
Helper function that checks if a given point is inside a polygon (a set of sorted vertices that conne...
void add_connection_matrix_info_helper(TriangleMeshPolyLine *polyline_pt, std::map< unsigned, std::map< unsigned, Vector< vertex_connection_info > > > &connection_matrix, TriangleMeshPolyLine *next_polyline_pt=0)
Helps to add information to the connection matrix of the given polyline.
void initialise_base_vertex(TriangleMeshPolyLine *polyline_pt, std::map< unsigned, std::map< unsigned, Vector< base_vertex_info > > > &base_vertices)
Initialise the base vertex structure, set every vertex to no visited and not being a base vertex.
std::map< unsigned, GeomObject * > & boundary_geom_object_pt()
Return direct access to the geometric object storage.
void add_base_vertex_info_helper(TriangleMeshPolyLine *polyline_pt, std::map< unsigned, std::map< unsigned, Vector< base_vertex_info > > > &base_vertices, std::map< unsigned, std::map< unsigned, Vector< vertex_connection_info > > > &connection_matrix, std::map< unsigned, std::map< unsigned, unsigned > > &boundary_chunk_nvertices)
Helps to identify the base vertex of the given polyline.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
double Tolerable_error
Acceptable discrepancy for mismatch in vertex coordinates. In paranoid mode, the code will die if the...
void create_triangulateio_from_polyfiles(const std::string &node_file_name, const std::string &element_file_name, const std::string &poly_file_name, TriangulateIO &triangle_io, bool &use_attributes)
Create a triangulateio data file from ele node and poly files. This is used if the mesh is generated ...
void dump_triangulateio(TriangulateIO &triangle_io, std::ostream &dump_file)
Write all the triangulateio data to disk in a dump file that can then be used to restart simulations.
void initialise_triangulateio(TriangulateIO &triangle_io)
Initialise TriangulateIO structure.
void read_triangulateio(std::istream &restart_file, TriangulateIO &triangle_io)
Read the triangulateio data from a dump file on disk, which can then be used to restart simulations.
void write_triangulateio_to_polyfile(TriangulateIO &triangle_io, std::ostream &poly_file)
Write the triangulateio data to disk as a poly file, mainly used for debugging.
void clear_triangulateio(TriangulateIO &triangulate_io, const bool &clear_hole_data)
Clear TriangulateIO structure.
TriangulateIO deep_copy_of_triangulateio_representation(TriangulateIO &triangle_io, const bool &quiet)
Make (partial) deep copy of TriangulateIO object. We only copy those items we need within oomph-lib's...
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
The Triangle data structure, modified from the triangle.h header supplied with triangle 1....
double * pointlist
Pointer to list of points x coordinate followed by y coordinate.
int * pointmarkerlist
Pointer to list of point markers.
double * pointattributelist
Pointer to list of point attributes.
Data structure to store the base vertex info, initial or final vertex in the polylines have an associ...
Data structure filled when the connection matrix is created, for each polyline, there are two vertex_...