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