triangle_scaffold_mesh.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #include "triangle_scaffold_mesh.h"
27 
28 
29 namespace oomph
30 {
31  //=====================================================================
32  /// Check mesh integrity -- performs some internal consistency checks
33  /// and throws error if violated.
34  //=====================================================================
36  {
37  // Check that all nodes are attached to elements
38  std::map<Node*, bool> done;
39 
40  // Number of nodes per element
41  unsigned nnod_el = finite_element_pt(0)->nnode();
42 
43  // Loop over elements to visit their nodes
44  unsigned nelem = nelement();
45  for (unsigned e = 0; e < nelem; e++)
46  {
47  // Loop over all nodes in element
48  for (unsigned j = 0; j < nnod_el; j++)
49  {
50  // Pointer to node in the scaffold mesh
52 
53  done[node_pt] = true;
54  }
55  }
56 
57  // Now check if all nodes have been visited
58  std::ostringstream error_stream;
59  bool broken = false;
60  unsigned nnod = nnode();
61  for (unsigned j = 0; j < nnod; j++)
62  {
63  if (!done[Node_pt[j]])
64  {
65  error_stream << "Scaffold node " << j
66  << " does not seem to be attached to any element";
67  broken = true;
68  }
69  }
70  if (broken)
71  {
72  throw OomphLibError(
73  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
74  }
75  }
76 
77  //=====================================================================
78  /// Constructor: Pass the filenames of the triangle files
79  /// The assumptions are that the nodes have been assigned boundary
80  /// information which is used in the nodal construction to ensure that
81  /// BoundaryNodes are indeed constructed when necessary. Additional
82  /// boundary information is added from the segment boundaries.
83  //=====================================================================
85  const std::string& ele_file_name,
86  const std::string& poly_file_name)
87  {
88  // Process element file
89  //---------------------
90  std::ifstream element_file(ele_file_name.c_str(), std::ios_base::in);
91 
92  // Check that the file actually opened correctly
93  if (!element_file.is_open())
94  {
95  std::string error_msg("Failed to open element file: ");
96  error_msg += "\"" + ele_file_name + "\".";
97  throw OomphLibError(
98  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
99  }
100 
101  // Number of elements
102  unsigned n_element;
103  element_file >> n_element;
104 
105  // Number of nodes per element
106  unsigned n_local_node;
107  element_file >> n_local_node;
108  if (n_local_node != 3)
109  {
110  std::ostringstream error_stream;
111  error_stream
112  << "Triangle should only be used to generate 3-noded triangles!\n"
113  << "Your triangle input file, contains data for " << n_local_node
114  << "-noded triangles" << std::endl;
115 
116  throw OomphLibError(
117  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
118  }
119 
120  // Dummy nodal attributes
121  unsigned dummy_attribute;
122 
123  // Element attributes may be used if we have internal boundaries
124  Element_attribute.resize(n_element, 0.0);
125 
126  // Dummy storage for element numbers
127  unsigned dummy_element_number;
128 
129  // Resize stoorage for global node numbers listed element-by-element
130  Global_node.resize(n_element * n_local_node);
131 
132 #ifdef PARANOID
133  std::map<unsigned, bool> global_node_done;
134 #endif
135 
136  // Initialise counter
137  unsigned k = 0;
138 
139  // Are attributes specified?
140  unsigned attribute_flag;
141  element_file >> attribute_flag;
142 
143  // Read global node numbers for all elements
144  if (attribute_flag == 0)
145  {
146  for (unsigned i = 0; i < n_element; i++)
147  {
148  element_file >> dummy_element_number;
149  for (unsigned j = 0; j < n_local_node; j++)
150  {
151  element_file >> Global_node[k];
152 #ifdef PARANOID
153  global_node_done[Global_node[k] - 1] = true;
154 #endif
155  k++;
156  }
157  }
158  }
159  else
160  {
161  for (unsigned i = 0; i < n_element; i++)
162  {
163  element_file >> dummy_element_number;
164  for (unsigned j = 0; j < n_local_node; j++)
165  {
166  element_file >> Global_node[k];
167 #ifdef PARANOID
168  global_node_done[Global_node[k] - 1] = true;
169 #endif
170  k++;
171  }
172  element_file >> Element_attribute[i];
173  }
174  }
175  element_file.close();
176 
177 #ifdef PARANOID
178  // Determine if the node numbering starts at 0 or 1 (triangle can do
179  // either depending on input arguments). We can't (currently) handle
180  // 0-based indexing so throw an error if it is being used.
181  if (*std::min_element(Global_node.begin(), Global_node.end()) != 1)
182  {
183  std::string err("Triangle is using 0-based indexing, ");
184  err += "however the oomph-lib interface can only handle 1-based indexing "
185  "in Triangle.\n";
186  err +=
187  "This is likely to be due to the input file using 0-based indexing.";
188  err += "Alternatively you may have specified the -z flag when running "
189  "Triangle.";
190  throw OomphLibError(
191  err, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
192  }
193 #endif
194 
195  // Resize the Element vector
196  Element_pt.resize(n_element);
197 
198  // Process node file
199  // -----------------
200  std::ifstream node_file(node_file_name.c_str(), std::ios_base::in);
201 
202  // Check that the file actually opened correctly
203  if (!node_file.is_open())
204  {
205  std::string error_msg("Failed to open node file: ");
206  error_msg += "\"" + node_file_name + "\".";
207  throw OomphLibError(
208  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
209  }
210 
211  // Read number of nodes
212  unsigned n_node;
213  node_file >> n_node;
214 
215  // Create a vector of boolean so as not to create the same node twice
216  std::vector<bool> done(n_node, false);
217 
218  // Resize the Node vector
219  Node_pt.resize(n_node);
220 
221  // Spatial dimension of nodes
222  unsigned dimension;
223  node_file >> dimension;
224 
225 #ifdef PARANOID
226  if (dimension != 2)
227  {
228  throw OomphLibError("The dimension must be 2\n",
229  OOMPH_CURRENT_FUNCTION,
230  OOMPH_EXCEPTION_LOCATION);
231  }
232 #endif
233 
234  // Flag for attributes
235  node_file >> attribute_flag;
236 
237  // Flag for boundary markers
238  unsigned boundary_markers_flag;
239  node_file >> boundary_markers_flag;
240 
241  // Dummy for node number
242  unsigned dummy_node_number;
243 
244  // Create storage for nodal posititions and boundary markers
245  Vector<double> x_node(n_node);
246  Vector<double> y_node(n_node);
247  Vector<unsigned> bound(n_node);
248 
249  // Check if there are attributes
250  if (attribute_flag == 0)
251  {
252  if (boundary_markers_flag == 1)
253  {
254  for (unsigned i = 0; i < n_node; i++)
255  {
256  node_file >> dummy_node_number;
257  node_file >> x_node[i];
258  node_file >> y_node[i];
259  node_file >> bound[i];
260  }
261  node_file.close();
262  }
263  else
264  {
265  for (unsigned i = 0; i < n_node; i++)
266  {
267  node_file >> dummy_node_number;
268  node_file >> x_node[i];
269  node_file >> y_node[i];
270  bound[i] = 0;
271  }
272  node_file.close();
273  }
274  }
275  else
276  {
277  if (boundary_markers_flag == 1)
278  {
279  for (unsigned i = 0; i < n_node; i++)
280  {
281  node_file >> dummy_node_number;
282  node_file >> x_node[i];
283  node_file >> y_node[i];
284  node_file >> dummy_attribute;
285  node_file >> bound[i];
286  }
287  node_file.close();
288  }
289  else
290  {
291  for (unsigned i = 0; i < n_node; i++)
292  {
293  node_file >> dummy_node_number;
294  node_file >> x_node[i];
295  node_file >> y_node[i];
296  node_file >> dummy_attribute;
297  bound[i] = 0;
298  }
299  node_file.close();
300  }
301  } // end
302 
303 
304  // Determine highest boundary index
305  // --------------------------------
306  unsigned n_bound = 0;
307  if (boundary_markers_flag == 1)
308  {
309  n_bound = bound[0];
310  for (unsigned i = 1; i < n_node; i++)
311  {
312  if (bound[i] > n_bound)
313  {
314  n_bound = bound[i];
315  }
316  }
317  }
318 
319  // Process poly file to extract edges
320  //-----------------------------------
321 
322  // Open poly file
323  std::ifstream poly_file(poly_file_name.c_str(), std::ios_base::in);
324 
325  // Check that the file actually opened correctly
326  if (!poly_file.is_open())
327  {
328  std::string error_msg("Failed to open poly file: ");
329  error_msg += "\"" + poly_file_name + "\".";
330  throw OomphLibError(
331  error_msg, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
332  }
333 
334  // Number of nodes in poly file
335  unsigned n_node_poly;
336  poly_file >> n_node_poly;
337 
338  // Dimension
339  poly_file >> dimension;
340 
341  // Attribute flag
342  poly_file >> attribute_flag;
343 
344  // Boundary markers flag
345  poly_file >> boundary_markers_flag;
346 
347 
348  // Ignore node information: Note: No, we can't extract the
349  // actual nodes themselves from here!
350  unsigned dummy;
351  if (n_node_poly > 0)
352  {
353  if (attribute_flag == 0)
354  {
355  if (boundary_markers_flag == 1)
356  {
357  for (unsigned i = 0; i < n_node_poly; i++)
358  {
359  poly_file >> dummy;
360  poly_file >> dummy;
361  poly_file >> dummy;
362  poly_file >> dummy;
363  }
364  }
365  else
366  {
367  for (unsigned i = 0; i < n_node_poly; i++)
368  {
369  poly_file >> dummy;
370  poly_file >> dummy;
371  poly_file >> dummy;
372  }
373  }
374  }
375  else
376  {
377  if (boundary_markers_flag == 1)
378  {
379  for (unsigned i = 0; i < n_node_poly; i++)
380  {
381  poly_file >> dummy;
382  poly_file >> dummy;
383  poly_file >> dummy;
384  poly_file >> dummy;
385  poly_file >> dummy;
386  }
387  }
388  else
389  {
390  for (unsigned i = 0; i < n_node_poly; i++)
391  {
392  poly_file >> dummy;
393  poly_file >> dummy;
394  poly_file >> dummy;
395  poly_file >> dummy;
396  }
397  }
398  }
399  }
400 
401  // Now extract the segment information
402  // Segements are lines that lie on boundaries of the domain
403  //----------------------------------------------------------
404 
405  // Number of segments
406  unsigned n_segment;
407  poly_file >> n_segment;
408 
409  // Boundary marker flag
410  poly_file >> boundary_markers_flag;
411 
412  // Storage for the global node numbers (in the triangle 1-based
413  // numbering scheme!) of the first and second node in each segment
414  Vector<unsigned> first_node(n_segment);
415  Vector<unsigned> second_node(n_segment);
416 
417  // Storage for the boundary marker for each segment
418  Vector<unsigned> segment_boundary(n_segment);
419 
420  // Dummy for global segment number
421  unsigned dummy_segment_number;
422 
423  // Storage for the edges associated with each node. Nodes are indexed
424  // using the Triangle 1-based index which is why there is a +1 here.
425  Vector<std::set<unsigned>> node_on_edges(n_node + 1);
426 
427  // Extract information for each segment
428  for (unsigned i = 0; i < n_segment; i++)
429  {
430  poly_file >> dummy_segment_number;
431  poly_file >> first_node[i];
432  poly_file >> second_node[i];
433  poly_file >> segment_boundary[i];
434  // Check that we don't have a higher segment boundary number
435  if (segment_boundary[i] > n_bound)
436  {
437  n_bound = segment_boundary[i];
438  }
439  // Add the segment index to each node
440  node_on_edges[first_node[i]].insert(i);
441  node_on_edges[second_node[i]].insert(i);
442  }
443 
444  // Extract hole center information
445  unsigned nhole = 0;
446  poly_file >> nhole;
447  if (nhole != 0)
448  {
449  Hole_centre.resize(nhole);
450 
451  // Dummy for hole number
452  unsigned dummy_hole;
453  // Loop over the holes to get centre coords
454  for (unsigned ihole = 0; ihole < nhole; ihole++)
455  {
456  Hole_centre[ihole].resize(2);
457 
458  // Read the centre value
459  poly_file >> dummy_hole;
460  poly_file >> Hole_centre[ihole][0];
461  poly_file >> Hole_centre[ihole][1];
462  }
463  }
464  else
465  {
466  Hole_centre.resize(0);
467  }
468  poly_file.close();
469 
470  // Set the number of boundaries
471  if (n_bound > 0)
472  {
473  this->set_nboundary(n_bound);
474  }
475 
476 
477  // Create the elements
478  //--------------------
479 
480  // Counter for nodes in the vector that lists
481  // the global node numbers of the elements' local nodes
482  unsigned counter = 0;
483  for (unsigned e = 0; e < n_element; e++)
484  {
485  Element_pt[e] = new TElement<2, 2>;
486  for (unsigned j = 0; j < n_local_node; j++)
487  {
488  unsigned global_node_number = Global_node[counter];
489  if (done[global_node_number - 1] == false) //... -1 because node number
490  // begins at 1 in triangle
491  {
492  // If we are on the boundary
493  if ((boundary_markers_flag == 1) &&
494  (bound[global_node_number - 1] > 0))
495  {
496  // Construct a boundary node
499  // Add to the boundary node look-up scheme
500  add_boundary_node(bound[global_node_number - 1] - 1,
502  }
503  // Otherwise make an ordinary node
504  else
505  {
508  }
509  done[global_node_number - 1] = true;
510  Node_pt[global_node_number - 1]->x(0) =
511  x_node[global_node_number - 1];
512  Node_pt[global_node_number - 1]->x(1) =
513  y_node[global_node_number - 1];
514  }
515  else
516  {
518  }
519  counter++;
520  }
521  }
522 
523  // Resize the "matrix" that stores the boundary id for each
524  // edge in each element.
525  Edge_boundary.resize(n_element);
526  Edge_index.resize(n_element);
527 
528  // Storage for the global node numbers (in triangle's 1-based
529  // numbering scheme) for the zero-th, 1st, and 2nd node in each
530  // triangle.
531  unsigned glob_num[3] = {0, 0, 0};
532 
533  // 0-based index used to construct a global index-based lookup scheme
534  // for each edge that will be used to uniquely construct mid-side
535  // nodes.
536  // The segments (edges that lie on boundaries) have already
537  // been added to the scheme, so we start with the number of segments.
538  Nglobal_edge = n_segment;
539 
540  // Loop over the elements
541  for (unsigned e = 0; e < n_element; e++)
542  {
543  // Each element has three edges
544  Edge_boundary[e].resize(3);
545  Edge_index[e].resize(3);
546  // By default each edge is NOT on a boundary
547  for (unsigned i = 0; i < 3; i++)
548  {
549  Edge_boundary[e][i] = 0;
550  }
551 
552  // Read out the global node numbers from the triangle data structure
553  const unsigned element_offset = e * n_local_node;
554  for (unsigned i = 0; i < 3; i++)
555  {
556  glob_num[i] = Global_node[element_offset + i];
557  }
558 
559  // Now we know the global node numbers of the elements' three nodes
560  // in triangle's 1-based numbering.
561 
562  // Determine whether any of the elements edges have already been
563  // allocated an index. This may be because they are on boundaries
564  // (segments) or because they have already occured.
565  // The global index for the i-th edge will be stored in edge_index[i]
566  for (unsigned i = 0; i < 3; i++)
567  {
568  std::vector<unsigned> local_edge_index;
569 
570  // Find the common global edge index for the nodes on
571  // the i-th element edge (note the use of moular arithmetic here)
572  std::set_intersection(node_on_edges[glob_num[i]].begin(),
573  node_on_edges[glob_num[i]].end(),
574  node_on_edges[glob_num[(i + 1) % 3]].begin(),
575  node_on_edges[glob_num[(i + 1) % 3]].end(),
576  std::insert_iterator<std::vector<unsigned>>(
577  local_edge_index, local_edge_index.begin()));
578 
579  // If the nodes share more than one global edge index, then
580  // we have a problem
581  if (local_edge_index.size() > 1)
582  {
583  throw OomphLibError(
584  "Nodes in scaffold mesh share more than one global edge",
585  OOMPH_CURRENT_FUNCTION,
586  OOMPH_EXCEPTION_LOCATION);
587  }
588 
589  // If the element's edge is not already allocated, the intersection
590  // will be empty
591  if (local_edge_index.size() == 0)
592  {
593  // Allocate the next global index
595  // Associate the new edge index with the nodes
596  node_on_edges[glob_num[i]].insert(Nglobal_edge);
597  node_on_edges[glob_num[(i + 1) % 3]].insert(Nglobal_edge);
598  // Increment the global edge index
599  ++Nglobal_edge;
600  }
601  // Otherwise we already have an edge
602  else if (local_edge_index.size() == 1)
603  {
604  // Set the edge index
605  Edge_index[e][i] = local_edge_index[0];
606  // Allocate the boundary index, if it is a segment
607  if (local_edge_index[0] < n_segment)
608  {
609  Edge_boundary[e][i] = segment_boundary[local_edge_index[0]];
610  // Add the nodes to the boundary look-up scheme in
611  // oomph-lib (0-based) index
612  add_boundary_node(segment_boundary[local_edge_index[0]] - 1,
613  Node_pt[glob_num[i] - 1]);
614  add_boundary_node(segment_boundary[local_edge_index[0]] - 1,
615  Node_pt[glob_num[(i + 1) % 3] - 1]);
616  }
617  }
618  }
619  }
620 
621 #ifdef PARANOID
622 
623  std::ostringstream error_stream;
624  bool broken = false;
625  unsigned nnod = nnode();
626  error_stream << "Checking presence of " << nnod << " global nodes\n";
627  for (unsigned j = 0; j < nnod; j++)
628  {
629  if (!global_node_done[j])
630  {
631  error_stream << "Global node " << j
632  << " was not listed in *.ele file\n";
633  broken = true;
634  }
635  }
636  if (broken)
637  {
638  throw OomphLibError(
639  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
640  }
641 
642  // Check things and throw if mesh is broken...
644 #endif
645  }
646 
647 #ifdef OOMPH_HAS_TRIANGLE_LIB
648 
649  //=====================================================================
650  /// Constructor: Pass a data structure obtained from the triangulate
651  /// function
652  //=====================================================================
654  {
655  // Number of elements
656  unsigned n_element = static_cast<unsigned>(triangle_data.numberoftriangles);
657 
658  // Number of nodes per element
659  unsigned n_local_node =
660  static_cast<unsigned>(triangle_data.numberofcorners);
661  if (n_local_node != 3)
662  {
663  std::ostringstream error_stream;
664  error_stream
665  << "Triangle should only be used to generate 3-noded triangles!\n"
666  << "Your triangle input file, contains data for " << n_local_node
667  << "-noded triangles" << std::endl;
668 
669  throw OomphLibError(
670  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
671  }
672 
673  // Element attributes may be used if we have internal boundaries
674  Element_attribute.resize(n_element, 0.0);
675 
676  // Resizestoorage for global node numbers listed element-by-element
677  Global_node.resize(n_element * n_local_node);
678 
679 #ifdef PARANOID
680  std::map<unsigned, bool> global_node_done;
681 #endif
682 
683  // Initialise counter
684  unsigned k = 0;
685 
686  // Are attributes specified?
687  unsigned attribute_flag =
688  static_cast<unsigned>(triangle_data.numberoftriangleattributes);
689 
690  // Read global node numbers for all elements
691  if (attribute_flag == 0)
692  {
693  for (unsigned i = 0; i < n_element; i++)
694  {
695  for (unsigned j = 0; j < n_local_node; j++)
696  {
697  Global_node[k] = static_cast<unsigned>(triangle_data.trianglelist[k]);
698 #ifdef PARANOID
699  global_node_done[Global_node[k] - 1] = true;
700 #endif
701  k++;
702  }
703  }
704  }
705  else
706  {
707  for (unsigned i = 0; i < n_element; i++)
708  {
709  for (unsigned j = 0; j < n_local_node; j++)
710  {
711  Global_node[k] = static_cast<unsigned>(triangle_data.trianglelist[k]);
712 #ifdef PARANOID
713  global_node_done[Global_node[k] - 1] = true;
714 #endif
715  k++;
716  }
717  Element_attribute[i] = triangle_data.triangleattributelist[i];
718  }
719  }
720 
721  // Resize the Element vector
722  Element_pt.resize(n_element);
723 
724  // Read number of nodes
725  unsigned n_node = triangle_data.numberofpoints;
726 
727  // Create a vector of boolean so as not to create the same node twice
728  std::vector<bool> done(n_node, false);
729 
730  // Resize the Node vector
731  Node_pt.resize(n_node);
732 
733  // Flag for boundary markers
734  unsigned boundary_markers_flag = 0;
735  if (triangle_data.pointmarkerlist != 0)
736  {
737  boundary_markers_flag = 1;
738  }
739 
740  // Create storage for nodal posititions and boundary markers
741  Vector<double> x_node(n_node);
742  Vector<double> y_node(n_node);
743  Vector<unsigned> bound(n_node);
744 
745  // We shall ingnore all point attributes
746  if (boundary_markers_flag == 1)
747  {
748  for (unsigned i = 0; i < n_node; i++)
749  {
750  x_node[i] = triangle_data.pointlist[2 * i];
751  y_node[i] = triangle_data.pointlist[2 * i + 1];
752  bound[i] = static_cast<unsigned>(triangle_data.pointmarkerlist[i]);
753  }
754  }
755  else
756  {
757  for (unsigned i = 0; i < n_node; i++)
758  {
759  x_node[i] = triangle_data.pointlist[2 * i];
760  y_node[i] = triangle_data.pointlist[2 * i + 1];
761  bound[i] = 0;
762  }
763  }
764 
765  // Determine highest boundary index
766  // --------------------------------
767  unsigned n_bound = 0;
768  if (boundary_markers_flag == 1)
769  {
770  n_bound = bound[0];
771  for (unsigned i = 1; i < n_node; i++)
772  {
773  if (bound[i] > n_bound)
774  {
775  n_bound = bound[i];
776  }
777  }
778  }
779 
780 
781  // Now extract the segment information
782  //------------------------------------
783 
784  // Number of segments
785  unsigned n_segment = triangle_data.numberofsegments;
786 
787  // Storage for the global node numbers (in the triangle 1-based
788  // numbering scheme!) of the first and second node in each segment
789  Vector<unsigned> first_node(n_segment);
790  Vector<unsigned> second_node(n_segment);
791 
792  // Storage for the boundary marker for each segment
793  Vector<unsigned> segment_boundary(n_segment);
794 
795  // Storage for the edges associated with each node. Nodes are indexed
796  // using the Triangle 1-based index which is why there is a +1 here.
797  Vector<std::set<unsigned>> node_on_edges(n_node + 1);
798 
799  // Extract information for each segment
800  for (unsigned i = 0; i < n_segment; i++)
801  {
802  first_node[i] = static_cast<unsigned>(triangle_data.segmentlist[2 * i]);
803  second_node[i] =
804  static_cast<unsigned>(triangle_data.segmentlist[2 * i + 1]);
805  segment_boundary[i] =
806  static_cast<unsigned>(triangle_data.segmentmarkerlist[i]);
807  // Check that we don't have a higher segment boundary number
808  if (segment_boundary[i] > n_bound)
809  {
810  n_bound = segment_boundary[i];
811  }
812  // Add the segment index to each node
813  node_on_edges[first_node[i]].insert(i);
814  node_on_edges[second_node[i]].insert(i);
815  }
816 
817  // Extract hole center information
818  unsigned nhole = triangle_data.numberofholes;
819  if (nhole != 0)
820  {
821  Hole_centre.resize(nhole);
822 
823  // Coords counter
824  unsigned count_coords = 0;
825 
826  // Loop over the holes to get centre coords
827  for (unsigned ihole = 0; ihole < nhole; ihole++)
828  {
829  Hole_centre[ihole].resize(2);
830 
831  // Read the centre value
832  Hole_centre[ihole][0] = triangle_data.holelist[count_coords];
833  Hole_centre[ihole][1] = triangle_data.holelist[count_coords + 1];
834 
835  // Increment coords counter
836  count_coords += 2;
837  }
838  }
839  else
840  {
841  Hole_centre.resize(0);
842  }
843 
844  // Set number of boundaries
845  if (n_bound > 0)
846  {
847  set_nboundary(n_bound);
848  }
849 
850 
851  // Create the elements
852  //--------------------
853 
854  // Counter for nodes in the vector that lists
855  // the global node numbers of the elements' local nodes
856  unsigned counter = 0;
857  for (unsigned e = 0; e < n_element; e++)
858  {
859  Element_pt[e] = new TElement<2, 2>;
860  for (unsigned j = 0; j < n_local_node; j++)
861  {
862  unsigned global_node_number = Global_node[counter];
863  if (done[global_node_number - 1] == false) //... -1 because node number
864  // begins at 1 in triangle
865  {
866  // If we are on the boundary
867  if ((boundary_markers_flag == 1) &&
868  (bound[global_node_number - 1] > 0))
869  {
870  // Construct a boundary node
873 
874  // Add to the boundary node look-up scheme
875  add_boundary_node(bound[global_node_number - 1] - 1,
877  }
878  // Otherwise make an ordinary node
879  else
880  {
883  }
884  done[global_node_number - 1] = true;
885  Node_pt[global_node_number - 1]->x(0) =
886  x_node[global_node_number - 1];
887  Node_pt[global_node_number - 1]->x(1) =
888  y_node[global_node_number - 1];
889  }
890  else
891  {
893  }
894  counter++;
895  }
896  }
897 
898  // Resize the "matrix" that stores the boundary id for each
899  // edge in each element.
900  Edge_boundary.resize(n_element);
901  Edge_index.resize(n_element);
902 
903  // Storage for the global node numbers (in triangle's 1-based
904  // numbering scheme) for the zero-th, 1st, and 2nd node in each
905  // triangle.
906  unsigned glob_num[3] = {0, 0, 0};
907 
908  // 0-based index used to construct a global index-based lookup scheme
909  // for each edge that will be used to uniquely construct mid-side
910  // nodes.
911  // The segments (edges that lie on boundaries) have already
912  // been added to the scheme, so we start with the number of segments.
913  Nglobal_edge = n_segment;
914 
915  // Loop over the elements
916  for (unsigned e = 0; e < n_element; e++)
917  {
918  // Each element has three edges
919  Edge_boundary[e].resize(3);
920  Edge_index[e].resize(3);
921  // By default each edge is NOT on a boundary
922  for (unsigned i = 0; i < 3; i++)
923  {
924  Edge_boundary[e][i] = 0;
925  }
926 
927  // Read out the global node numbers from the triangle data structure
928  const unsigned element_offset = e * n_local_node;
929  for (unsigned i = 0; i < 3; i++)
930  {
931  glob_num[i] = Global_node[element_offset + i];
932  }
933 
934  // Now we know the global node numbers of the elements' three nodes
935  // in triangle's 1-based numbering.
936 
937  // Determine whether any of the elements edges have already been
938  // allocated an index. This may be because they are on boundaries
939  // (segments) or because they have already occured.
940  // The global index for the i-th edge will be stored in edge_index[i]
941  for (unsigned i = 0; i < 3; i++)
942  {
943  std::vector<unsigned> local_edge_index;
944 
945  // Find the common global edge index for the nodes on
946  // the i-th element edge (note the use of moular arithmetic here)
947  std::set_intersection(node_on_edges[glob_num[i]].begin(),
948  node_on_edges[glob_num[i]].end(),
949  node_on_edges[glob_num[(i + 1) % 3]].begin(),
950  node_on_edges[glob_num[(i + 1) % 3]].end(),
951  std::insert_iterator<std::vector<unsigned>>(
952  local_edge_index, local_edge_index.begin()));
953 
954  // If the nodes share more than one global edge index, then
955  // we have a problem
956  if (local_edge_index.size() > 1)
957  {
958  throw OomphLibError(
959  "Nodes in scaffold mesh share more than one global edge",
960  OOMPH_CURRENT_FUNCTION,
961  OOMPH_EXCEPTION_LOCATION);
962  }
963 
964 
965  // If the element's edge is not already allocated, the intersection
966  // will be empty
967  if (local_edge_index.size() == 0)
968  {
969  // Allocate the next global index
971  // Associate the new edge index with the nodes
972  node_on_edges[glob_num[i]].insert(Nglobal_edge);
973  node_on_edges[glob_num[(i + 1) % 3]].insert(Nglobal_edge);
974  // Increment the global edge index
975  ++Nglobal_edge;
976  }
977  // Otherwise we already have an edge
978  else if (local_edge_index.size() == 1)
979  {
980  // Set the edge index
981  Edge_index[e][i] = local_edge_index[0];
982  // Allocate the boundary index, if it is a segment
983  if (local_edge_index[0] < n_segment)
984  {
985  Edge_boundary[e][i] = segment_boundary[local_edge_index[0]];
986  // Add the nodes to the boundary look-up scheme in
987  // oomph-lib (0-based) index
988  add_boundary_node(segment_boundary[local_edge_index[0]] - 1,
989  Node_pt[glob_num[i] - 1]);
990  add_boundary_node(segment_boundary[local_edge_index[0]] - 1,
991  Node_pt[glob_num[(i + 1) % 3] - 1]);
992  }
993  }
994  }
995  }
996 
997 #ifdef PARANOID
998 
999  std::ostringstream error_stream;
1000  bool broken = false;
1001  unsigned nnod = nnode();
1002  error_stream << "Checking presence of " << nnod << " global nodes\n";
1003  for (unsigned j = 0; j < nnod; j++)
1004  {
1005  if (!global_node_done[j])
1006  {
1007  error_stream << "Global node " << j
1008  << " was not listed in *.ele file\n";
1009  broken = true;
1010  }
1011  }
1012  if (broken)
1013  {
1014  error_stream
1015  << "This error means that some of the nodes are not connected \n"
1016  << " to (bulk) elements. This can happen if there is an isolated\n"
1017  << " boundary line in the mesh. One possible cause for this is\n"
1018  << " specifying a hole coordinate in the wrong place so that there\n"
1019  << " a gap between the mesh and the outer boundary.\n";
1020  throw OomphLibError(
1021  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1022  }
1023 
1024  // Check things and throw if mesh is broken...
1026 #endif
1027  }
1028 
1029 #endif
1030 
1031 } // namespace oomph
e
Definition: cfortran.h:571
cstr elem_len * i
Definition: cfortran.h:603
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition: elements.h:2509
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2538
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition: mesh.cc:243
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:183
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:473
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition: mesh.h:505
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:596
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:436
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition: mesh.h:186
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:590
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
An OomphLibError object which should be thrown when an run-time error is encountered....
General TElement class.
Definition: Telements.h:1208
unsigned global_node_number(const unsigned &i)
Return the global node of each local node listed element-by-element e*n_local_node + n_local Note tha...
Vector< unsigned > Global_node
Storage for global node numbers listed element-by-element.
void check_mesh_integrity()
Check mesh integrity – performs some internal consistency checks and throws error if violated.
unsigned Nglobal_edge
Number of internal edges.
Vector< double > Element_attribute
Vector of double attributes for each element.
TriangleScaffoldMesh()
Empty constructor.
Vector< Vector< unsigned > > Edge_boundary
Vector of vectors containing the boundary ids of the elements' edges.
Vector< Vector< double > > Hole_centre
Vectors of hole centre coordinates.
Vector< Vector< unsigned > > Edge_index
Vector of vectors containing the global edge index of.
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
The Triangle data structure, modified from the triangle.h header supplied with triangle 1....
double * pointlist
Pointer to list of points x coordinate followed by y coordinate.
int * pointmarkerlist
Pointer to list of point markers.