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