tetgen_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: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
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 "mesh.h"
31 #include "Telements.h"
32 #include "tetgen_scaffold_mesh.h"
33 
34 namespace oomph
35 {
36 
37 //======================================================================
38 /// Constructor: Pass the filename of the tetrahedra file
39 /// The assumptions are that the nodes have been assigned boundary
40 /// information which is used in the nodal construction to make sure
41 /// that BoundaryNodes are constructed. Any additional boundaries are
42 /// determined from the face boundary information.
43 //======================================================================
45  const std::string& element_file_name,
46  const std::string& face_file_name)
47 {
48 
49  // Process the element file
50  // --------------------------
51  std::ifstream element_file(element_file_name.c_str(),std::ios_base::in);
52 
53  // Check that the file actually opened correctly
54 #ifdef PARANOID
55  if(!element_file.is_open())
56  {
57  std::string error_msg("Failed to open element file: ");
58  error_msg += "\"" + element_file_name + "\".";
59  throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
60  OOMPH_EXCEPTION_LOCATION);
61  }
62 #endif
63 
64 
65  //Read in number of elements
66  unsigned n_element;
67  element_file>>n_element;
68 
69  //Read in number of nodes per element
70  unsigned n_local_node;
71  element_file>>n_local_node;
72 
73  //Throw an error if we have anything but linear simplices
74  if (n_local_node!=4)
75  {
76  std::ostringstream error_stream;
77  error_stream
78  << "Tetgen should only be used to generate 4-noded tetrahedra!\n"
79  << "Your tetgen input file, contains data for "
80  << n_local_node << "-noded tetrahedra" << std::endl;
81 
82  throw OomphLibError(error_stream.str(),
83  OOMPH_CURRENT_FUNCTION,
84  OOMPH_EXCEPTION_LOCATION);
85  }
86 
87  // Dummy nodal attributes
88  unsigned dummy_attribute;
89 
90  // Element attributes may be used to distinguish internal regions
91  Element_attribute.resize(n_element,0.0);
92 
93  // Dummy storage for element numbers
94  unsigned dummy_element_number;
95 
96  // Resize storage for the global node numbers listed element-by-element
97  Global_node.resize(n_element*n_local_node);
98 
99  //Initialise (global) node counter
100  unsigned k=0;
101  //Are there attributes?
102  unsigned attribute_flag;
103  element_file >> attribute_flag;
104 
105  //If no attributes
106  if(attribute_flag==0)
107  {
108  //Read in global node numbers
109  for(unsigned i=0;i<n_element;i++)
110  {
111  element_file>>dummy_element_number;
112  for(unsigned j=0;j<n_local_node;j++)
113  {
114  element_file >> Global_node[k];
115  k++;
116  }
117  }
118  }
119  //Otherwise read in the attributes as well
120  else
121  {
122  for(unsigned i=0;i<n_element;i++)
123  {
124  element_file>>dummy_element_number;
125  for(unsigned j=0;j<n_local_node;j++)
126  {
127  element_file >> Global_node[k];
128  k++;
129  }
130  element_file >> Element_attribute[i];
131  }
132  }
133  element_file.close();
134 
135  // Resize the Element vector
136  Element_pt.resize(n_element);
137 
138  //Process node file
139  //--------------------
140  std::ifstream node_file(node_file_name.c_str(),std::ios_base::in);
141 
142  // Check that the file actually opened correctly
143 #ifdef PARANOID
144  if(!node_file.is_open())
145  {
146  std::string error_msg("Failed to open node file: ");
147  error_msg += "\"" + node_file_name + "\".";
148  throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
149  OOMPH_EXCEPTION_LOCATION);
150  }
151 #endif
152 
153 
154  //Read in the number of nodes
155  unsigned n_node;
156  node_file>>n_node;
157 
158  // Create a vector of boolean so as not to create the same node twice
159  std::vector<bool> done (n_node,false);
160 
161  // Resize the Node vector
162  Node_pt.resize(n_node);
163 
164  //Set the spatial dimension of the nodes
165  unsigned dimension;
166  node_file>>dimension;
167 
168 #ifdef PARANOID
169  if(dimension!=3)
170  {
171  throw OomphLibError("The dimesion of the nodes must be 3\n",
172  OOMPH_CURRENT_FUNCTION,
173  OOMPH_EXCEPTION_LOCATION);
174  }
175 #endif
176 
177  // Flag for attributes
178  node_file>>attribute_flag;
179 
180  //Flag for boundary markers
181  unsigned boundary_markers_flag;
182  node_file>>boundary_markers_flag;
183 
184  // Dummy storage for the node number
185  unsigned dummy_node_number;
186 
187  // Create storage for nodal positions and boundary markers
188  Vector<double> x_node(n_node);
189  Vector<double> y_node(n_node);
190  Vector<double> z_node(n_node);
191  Vector<unsigned> bound(n_node);
192 
193  // Check if there are attributes
194  //If not
195  if(attribute_flag==0)
196  {
197  //Are there boundary markers
198  if(boundary_markers_flag==1)
199  {
200  for(unsigned i=0;i<n_node;i++)
201  {
202  node_file>>dummy_node_number;
203  node_file>>x_node[i];
204  node_file>>y_node[i];
205  node_file>>z_node[i];
206  node_file>>bound[i];
207  }
208  node_file.close();
209  }
210  //Otherwise no boundary markers
211  else
212  {
213  for(unsigned i=0;i<n_node;i++)
214  {
215  node_file>>dummy_node_number;
216  node_file>>x_node[i];
217  node_file>>y_node[i];
218  node_file>>z_node[i];
219  bound[i] = 0;
220  }
221  node_file.close();
222  }
223  }
224  //Otherwise there are attributes
225  else
226  {
227  if(boundary_markers_flag==1)
228  {
229  for(unsigned i=0;i<n_node;i++)
230  {
231  node_file>>dummy_node_number;
232  node_file>>x_node[i];
233  node_file>>y_node[i];
234  node_file>>z_node[i];
235  node_file>>dummy_attribute;
236  node_file>>bound[i];
237  }
238  node_file.close();
239  }
240  else
241  {
242  for(unsigned i=0;i<n_node;i++)
243  {
244  node_file>>dummy_node_number;
245  node_file>>x_node[i];
246  node_file>>y_node[i];
247  node_file>>z_node[i];
248  node_file>>dummy_attribute;
249  bound[i] = 0;
250  }
251  node_file.close();
252  }
253  }
254 
255  // Determine highest boundary index
256  //------------------------------------
257  unsigned n_bound=0;
258  if(boundary_markers_flag==1)
259  {
260  n_bound=bound[0];
261  for(unsigned i=1;i<n_node;i++)
262  {
263  if(bound[i] > n_bound)
264  {
265  n_bound = bound[i];
266  }
267  }
268  }
269 
270  // Process face file to extract boundary faces
271  //--------------------------------------------
272 
273  // Open face file
274  std::ifstream face_file(face_file_name.c_str(),std::ios_base::in);
275 
276  // Check that the file actually opened correctly
277 #ifdef PARANOID
278  if(!face_file.is_open())
279  {
280  std::string error_msg("Failed to open face file: ");
281  error_msg += "\"" + face_file_name + "\".";
282  throw OomphLibError(error_msg, OOMPH_CURRENT_FUNCTION,
283  OOMPH_EXCEPTION_LOCATION);
284  }
285 #endif
286 
287  // Number of faces in face file
288  unsigned n_face;
289  face_file>>n_face;
290 
291  // Boundary markers flag
292  face_file>>boundary_markers_flag;
293 
294  // Storage for the global node numbers (in the tetgen 1-based
295  // numbering scheme!) of the first, second and third node in
296  // each segment
297  Vector<unsigned> first_node(n_face);
298  Vector<unsigned> second_node(n_face);
299  Vector<unsigned> third_node(n_face);
300 
301  // Storage for the boundary marker for each face
303 
304  // Dummy for global face number
305  unsigned dummy_face_number;
306 
307  // Storage for the (boundary) faces associated with each node.
308  // Nodes are indexed using Tetgen's 1-based scheme, which is why
309  // there is a +1 here
310  Vector<std::set<unsigned> > node_on_faces(n_node+1);
311 
312  // Extract information for each segment
313  for(unsigned i=0;i<n_face;i++)
314  {
315  face_file >> dummy_face_number;
316  face_file >> first_node[i];
317  face_file >> second_node[i];
318  face_file >> third_node[i];
319  face_file >> face_boundary[i];
320  if (face_boundary[i]>n_bound) {n_bound=face_boundary[i];}
321  //Add the face index to each node
322  node_on_faces[first_node[i]].insert(i);
323  node_on_faces[second_node[i]].insert(i);
324  node_on_faces[third_node[i]].insert(i);
325  }
326  face_file.close();
327 
328  // Set number of boundaries
329  if(n_bound > 0)
330  {
331  this->set_nboundary(n_bound);
332  }
333 
334  // Create the elements
335  unsigned counter=0;
336  for(unsigned e=0;e<n_element;e++)
337  {
339  unsigned global_node_number = Global_node[counter];
340  if(done[global_node_number-1]==false)
341  // ... -1 because node number begins at 1 in tetgen
342  {
343  //If the node is on a boundary, construct a boundary node
344  if((boundary_markers_flag==1) &&
345  (bound[global_node_number-1] > 0))
346  {
347  //Construct the boundary ndoe
348  Node_pt[global_node_number-1] =
350  //Add to the boundary lookup scheme
351  add_boundary_node(bound[global_node_number-1]-1,
352  Node_pt[global_node_number-1]);
353  }
354  //Otherwise just construct a normal node
355  else
356  {
357  Node_pt[global_node_number-1] =
359  }
360 
361  done[global_node_number-1]=true;
362  Node_pt[global_node_number-1]->x(0)=x_node[global_node_number-1];
363  Node_pt[global_node_number-1]->x(1)=y_node[global_node_number-1];
364  Node_pt[global_node_number-1]->x(2)=z_node[global_node_number-1];
365  }
366  //Otherwise just copy the node numbr accross
367  else
368  {
369  finite_element_pt(e)->node_pt(3) = Node_pt[global_node_number-1];
370  }
371  counter++;
372 
373  //Loop over the other nodes
374  for(unsigned j=0;j<(n_local_node-1);j++)
375  {
376  global_node_number=Global_node[counter];
377  if(done[global_node_number-1]==false)
378  // ... -1 because node number begins at 1 in tetgen
379  {
380  //If we're on a boundary
381  if((boundary_markers_flag==1) &&
382  (bound[global_node_number-1] > 0))
383  {
384  //Construct the boundary ndoe
385  Node_pt[global_node_number-1] =
387  //Add to the boundary lookup scheme
388  add_boundary_node(bound[global_node_number-1]-1,
389  Node_pt[global_node_number-1]);
390  }
391  else
392  {
393  Node_pt[global_node_number-1]=finite_element_pt(e)->construct_node(j);
394  }
395  done[global_node_number-1]=true;
396  Node_pt[global_node_number-1]->x(0)=x_node[global_node_number-1];
397  Node_pt[global_node_number-1]->x(1)=y_node[global_node_number-1];
398  Node_pt[global_node_number-1]->x(2)=z_node[global_node_number-1];
399  }
400  //Otherwise copy the pointer over
401  else
402  {
403  finite_element_pt(e)->node_pt(j) = Node_pt[global_node_number-1];
404  }
405  counter++;
406  }
407  }
408 
409 
410  // Resize the "matrix" that stores the boundary id for each
411  // face in each element.
412  Face_boundary.resize(n_element);
413  Face_index.resize(n_element);
414  Edge_index.resize(n_element);
415 
416 
417  //0-based index scheme used to construct a global lookup for
418  //each face that will be used to uniquely construct interior
419  //face nodes.
420  //The faces that lie on boundaries will have already been allocated so
421  //we can start from there
422  Nglobal_face = n_face;
423 
424  // Storage for the edges associated with each node.
425  // Nodes are indexed using Tetgen's 1-based scheme, which is why
426  // there is a +1 here
427  Vector<std::set<unsigned> > node_on_edges(n_node+1);
428 
429  //0-based index scheme used to construct a global lookup for each
430  //edge that will be used to uniquely construct interior edge nodes
431  Nglobal_edge = 0;
432 
433  // Conversion from the edge numbers to the nodes at the end
434  // of each each edge
435  const unsigned first_local_edge_node[6] = {0,0,0,1,2,1};
436  const unsigned second_local_edge_node[6] = {1,2,3,2,3,3};
437 
438  // Loop over the elements
439  for(unsigned e=0;e<n_element;e++)
440  {
441  // Each element has four faces
442  Face_boundary[e].resize(4);
443  Face_index[e].resize(4);
444  // By default each face is NOT on a boundary
445  for(unsigned i=0;i<4;i++)
446  {
447  Face_boundary[e][i]=0;
448  }
449 
450  Edge_index[e].resize(6);
451 
452  // Read out the global node numbers of the nodes from
453  // tetgen's 1-based numbering.
454  // The offset is to match the offset used above
455  const unsigned element_offset = e*n_local_node;
456  unsigned glob_num[4] = {0,0,0,0};
457  for(unsigned i=0;i<4;++i)
458  {
459  glob_num[i] = Global_node[element_offset + ((i+1)%4)];
460  }
461 
462  // Now we know the global node numbers of the elements' four nodes
463  // in tetgen's 1-based numbering.
464 
465  //Determine whether any of the element faces have already been allocated an
466  //index. This may be because they are located on boundaries, or
467  //have already been visited
468  //The global index for the i-th face will be stored in Face_index[i]
469 
470  //Loop over the local faces in the element
471  for(unsigned i=0;i<4;++i)
472  {
473  //On the i-th face, our numbering convention is such that
474  //it is the (3-i)th node of the element that is omitted
475  const unsigned omitted_node = 3-i;
476 
477  //Start with the set of global faces associated with the i-th node
478  //which is always on the i-th face
479  std::set<unsigned> input = node_on_faces[glob_num[i]];
480 
481  //If there is no data yet, not point doing intersections
482  //the face has not been set up
483  if(input.size() > 0)
484  {
485  //Loop over the other nodes on the face
486  unsigned local_node = i;
487  for(unsigned i2=0;i2<3;++i2)
488  {
489  //Create the next node index (mod 4)
490  local_node = (local_node+1)%4;
491  //Only take the intersection of nodes on the face
492  //i.e. not 3-i
493  if(local_node != omitted_node)
494  {
495  //Local storage for the intersection
496  std::set<unsigned> intersection_result;
497  //Calculate the intersection of the current input
498  //and next node
499  std::set_intersection(input.begin(),input.end(),
500  node_on_faces[glob_num[local_node]].begin(),
501  node_on_faces[glob_num[local_node]].end(),
502  std::insert_iterator<std::set<unsigned> >(
503  intersection_result,
504  intersection_result.begin()));
505  //Let the result be the next input
506  input = intersection_result;
507  }
508  }
509  }
510 
511  //If the nodes share more than one global face index, then
512  //we have a problem
513  if(input.size() > 1)
514  {
515  throw OomphLibError(
516  "Nodes in scaffold mesh share more than one global face",
517  OOMPH_CURRENT_FUNCTION,
518  OOMPH_EXCEPTION_LOCATION);
519  }
520 
521  //If the element's face is not already allocated, the intersection
522  //will be empty
523  if(input.size()==0)
524  {
525  //Allocate the next global index
527  //Associate the new face index with the nodes
528  for(unsigned i2=0;i2<4;++i2)
529  {
530  //Don't add the omitted node
531  if(i2 != omitted_node)
532  {
533  node_on_faces[glob_num[i2]].insert(Nglobal_face);
534  }
535  }
536  //Increment the global face index
537  ++Nglobal_face;
538  }
539  //Otherwise we already have a face
540  else if(input.size()==1)
541  {
542  const unsigned global_face_index = *input.begin();
543  //Set the face index
544  Face_index[e][i] = global_face_index;
545  //Allocate the boundary index, if it's a boundary
546  if(global_face_index < n_face)
547  {
548  Face_boundary[e][i] = face_boundary[global_face_index];
549  //Add the nodes to the boundary look-up scheme in
550  //oomph-lib (0-based) index
551  for(unsigned i2=0;i2<4;++i2)
552  {
553  //Don't add the omitted node
554  if(i2 != omitted_node)
555  {
556  add_boundary_node(face_boundary[global_face_index]-1,
557  Node_pt[glob_num[i2]-1]);
558  }
559  }
560  }
561  }
562  } //end of loop over the faces
563 
564 
565  //Loop over the element edges and assign global edge numbers
566  for(unsigned i=0;i<6;++i)
567  {
568  //Storage for the result of the intersection
569  std::vector<unsigned> local_edge_index;
570 
571  const unsigned first_global_num = glob_num[first_local_edge_node[i]];
572  const unsigned second_global_num = glob_num[second_local_edge_node[i]];
573 
574  //Find the common global edge index for the nodes on
575  //the i-th element edge
576  std::set_intersection(node_on_edges[first_global_num].begin(),
577  node_on_edges[first_global_num].end(),
578  node_on_edges[second_global_num].begin(),
579  node_on_edges[second_global_num].end(),
580  std::insert_iterator<std::vector<unsigned> >(
581  local_edge_index,local_edge_index.begin()));
582 
583  //If the nodes share more than one global edge index, then
584  //we have a problem
585  if(local_edge_index.size() > 1)
586  {
587  throw OomphLibError(
588  "Nodes in scaffold mesh share more than one global edge",
589  OOMPH_CURRENT_FUNCTION,
590  OOMPH_EXCEPTION_LOCATION);
591  }
592 
593  //If the element's edge is not already allocated, the intersection
594  //will be empty
595  if(local_edge_index.size()==0)
596  {
597  //Allocate the next global index
599  //Associate the new edge index with the nodes
600  node_on_edges[first_global_num].insert(Nglobal_edge);
601  node_on_edges[second_global_num].insert(Nglobal_edge);
602  //Increment the global edge index
603  ++Nglobal_edge;
604  }
605  //Otherwise we already have an edge
606  else if(local_edge_index.size()==1)
607  {
608  //Set the edge index from the result of the intersection
609  Edge_index[e][i] = local_edge_index[0];
610  }
611  }
612 
613  } // end for e
614 
615 
616 
617  //Now determine whether any edges lie on boundaries by using the
618  //face boundary scheme and
619 
620  //Resize the storage
621  Edge_boundary.resize(Nglobal_edge,false);
622 
623  //Now loop over all the boundary faces and mark that all edges
624  //must also lie on the bounadry
625  for(unsigned i=0;i<n_face;++i)
626  {
627  {
628  //Storage for the result of the intersection
629  std::vector<unsigned> local_edge_index;
630 
631  //Find the common global edge index for first edge of the face
632  std::set_intersection(node_on_edges[first_node[i]].begin(),
633  node_on_edges[first_node[i]].end(),
634  node_on_edges[second_node[i]].begin(),
635  node_on_edges[second_node[i]].end(),
636  std::insert_iterator<std::vector<unsigned> >(
637  local_edge_index,local_edge_index.begin()));
638 
639  //If the nodes do not share exactly one global edge index, then
640  //we have a problem
641  if(local_edge_index.size() != 1)
642  {
643  throw OomphLibError(
644  "Nodes in scaffold mesh face do not share exactly one global edge",
645  OOMPH_CURRENT_FUNCTION,
646  OOMPH_EXCEPTION_LOCATION);
647  }
648  else
649  {
650  Edge_boundary[local_edge_index[0]] = true;
651  }
652  }
653 
654  {
655  //Storage for the result of the intersection
656  std::vector<unsigned> local_edge_index;
657 
658  //Find the common global edge index for second edge of the face
659  std::set_intersection(node_on_edges[second_node[i]].begin(),
660  node_on_edges[second_node[i]].end(),
661  node_on_edges[third_node[i]].begin(),
662  node_on_edges[third_node[i]].end(),
663  std::insert_iterator<std::vector<unsigned> >(
664  local_edge_index,local_edge_index.begin()));
665 
666  //If the nodes do not share exactly one global edge index, then
667  //we have a problem
668  if(local_edge_index.size() != 1)
669  {
670  throw OomphLibError(
671  "Nodes in scaffold mesh face do not share exactly one global edge",
672  OOMPH_CURRENT_FUNCTION,
673  OOMPH_EXCEPTION_LOCATION);
674  }
675  else
676  {
677  Edge_boundary[local_edge_index[0]] = true;
678  }
679  }
680 
681  {
682  //Storage for the result of the intersection
683  std::vector<unsigned> local_edge_index;
684 
685  //Find the common global edge index for third edge of the face
686  std::set_intersection(node_on_edges[first_node[i]].begin(),
687  node_on_edges[first_node[i]].end(),
688  node_on_edges[third_node[i]].begin(),
689  node_on_edges[third_node[i]].end(),
690  std::insert_iterator<std::vector<unsigned> >(
691  local_edge_index,local_edge_index.begin()));
692 
693  //If the nodes do not share exactly one global edge index, then
694  //we have a problem
695  if(local_edge_index.size() != 1)
696  {
697  throw OomphLibError(
698  "Nodes in scaffold mesh face do not share exactly one global edge",
699  OOMPH_CURRENT_FUNCTION,
700  OOMPH_EXCEPTION_LOCATION);
701  }
702  else
703  {
704  Edge_boundary[local_edge_index[0]] = true;
705  }
706  }
707  }
708 
709 
710 } //end of constructor
711 
712 
713 //======================================================================
714 /// Constructor: Pass a tetgenio data structure that represents
715 /// the input data of the mesh.
716 /// The assumptions are that the nodes have been assigned boundary
717 /// information which is used in the nodal construction to make sure
718 /// that BoundaryNodes are constructed. Any additional boundaries are
719 /// determined from the face boundary information.
720 //======================================================================
722 {
723  //Find the number of elements
724  unsigned n_element = static_cast<unsigned>(tetgen_data.numberoftetrahedra);
725 
726  //Read in number of nodes per element
727  unsigned n_local_node = static_cast<unsigned>(tetgen_data.numberofcorners);
728 
729  //Throw an error if we have anything but linear simplices
730  if (n_local_node!=4)
731  {
732  std::ostringstream error_stream;
733  error_stream
734  << "Tetgen should only be used to generate 4-noded tetrahedra!\n"
735  << "Your tetgen input data, contains data for "
736  << n_local_node << "-noded tetrahedra" << std::endl;
737 
738  throw OomphLibError(error_stream.str(),
739  OOMPH_CURRENT_FUNCTION,
740  OOMPH_EXCEPTION_LOCATION);
741  }
742 
743  // Element attributes may be used to distinguish internal regions
744  Element_attribute.resize(n_element,0.0);
745 
746  // Resize storage for the global node numbers listed element-by-element
747  Global_node.resize(n_element*n_local_node);
748 
749  //Initialise (global) node counter
750  unsigned k=0;
751  //Are there attributes?
752  unsigned attribute_flag =
753  static_cast<unsigned>(tetgen_data.numberoftetrahedronattributes);
754 
755  //Read global node numbers for all elements
756  if(attribute_flag==0)
757  {
758  for(unsigned i=0;i<n_element;i++)
759  {
760  for(unsigned j=0;j<n_local_node;j++)
761  {
762  Global_node[k] = static_cast<unsigned>(tetgen_data.tetrahedronlist[k]);
763  k++;
764  }
765  }
766  }
767  //Otherwise read in the attributes as well
768  else
769  {
770  for(unsigned i=0;i<n_element;i++)
771  {
772  for(unsigned j=0;j<n_local_node;j++)
773  {
774  Global_node[k] = static_cast<unsigned>(tetgen_data.tetrahedronlist[k]);
775  k++;
776  }
777  Element_attribute[i] = tetgen_data.tetrahedronattributelist[i];
778  }
779  }
780 
781  // Resize the Element vector
782  Element_pt.resize(n_element);
783 
784  //Read in the number of nodes
785  unsigned n_node = tetgen_data.numberofpoints;
786 
787  // Create a vector of boolean so as not to create the same node twice
788  std::vector<bool> done (n_node,false);
789 
790  // Resize the Node vector
791  Node_pt.resize(n_node);
792 
793  //Flag for boundary markers
794  unsigned boundary_markers_flag = 0;
795  if(tetgen_data.pointmarkerlist!=0) {boundary_markers_flag=1;}
796 
797  // Create storage for nodal positions and boundary markers
798  Vector<double> x_node(n_node);
799  Vector<double> y_node(n_node);
800  Vector<double> z_node(n_node);
801  Vector<unsigned> bound(n_node);
802 
803  //We shall ignore all point attributes
804  if(boundary_markers_flag==1)
805  {
806  for(unsigned i=0;i<n_node;i++)
807  {
808  x_node[i] = tetgen_data.pointlist[3*i];
809  y_node[i] = tetgen_data.pointlist[3*i+1];
810  z_node[i] = tetgen_data.pointlist[3*i+2];
811  bound[i] = static_cast<unsigned>(tetgen_data.pointmarkerlist[i]);
812  }
813  }
814  //Otherwise no boundary markers
815  else
816  {
817  for(unsigned i=0;i<n_node;i++)
818  {
819  x_node[i] = tetgen_data.pointlist[3*i];
820  y_node[i] = tetgen_data.pointlist[3*i+1];
821  z_node[i] = tetgen_data.pointlist[3*i+2];
822  bound[i] = 0;
823  }
824  }
825 
826 
827  // Determine highest boundary index
828  //------------------------------------
829  unsigned n_bound=0;
830  if(boundary_markers_flag==1)
831  {
832  n_bound=bound[0];
833  for(unsigned i=1;i<n_node;i++)
834  {
835  if(bound[i] > n_bound)
836  {
837  n_bound = bound[i];
838  }
839  }
840  }
841 
842  //Now extract face information
843  //---------------------------------
844 
845  // Number of faces in face file
846  unsigned n_face = tetgen_data.numberoftrifaces;
847 
848  // Storage for the global node numbers (in the tetgen 1-based
849  // numbering scheme!) of the first, second and third node in
850  // each segment
851  Vector<unsigned> first_node(n_face);
852  Vector<unsigned> second_node(n_face);
853  Vector<unsigned> third_node(n_face);
854 
855  // Storage for the boundary marker for each face
857 
858  // Storage for the (boundary) faces associated with each node.
859  // Nodes are indexed using Tetgen's 1-based scheme, which is why
860  // there is a +1 here
861  Vector<std::set<unsigned> > node_on_faces(n_node+1);
862 
863  // Extract information for each segment
864  for(unsigned i=0;i<n_face;i++)
865  {
866  first_node[i] = static_cast<unsigned>(tetgen_data.trifacelist[3*i]);
867  second_node[i] = static_cast<unsigned>(tetgen_data.trifacelist[3*i+1]);
868  third_node[i] = static_cast<unsigned>(tetgen_data.trifacelist[3*i+2]);
869  face_boundary[i] =
870  static_cast<unsigned>(tetgen_data.trifacemarkerlist[i]);
871  if (face_boundary[i]>n_bound) {n_bound=face_boundary[i];}
872  //Add the face index to each node
873  node_on_faces[first_node[i]].insert(i);
874  node_on_faces[second_node[i]].insert(i);
875  node_on_faces[third_node[i]].insert(i);
876  }
877 
878  // Extract hole center information
879  //unsigned nhole=tetgen_data.numberofholes;
880  /*if(nhole!=0)
881  {
882  Hole_centre.resize(nhole);
883 
884  // Coords counter
885  unsigned count_coords=0;
886 
887  // Loop over the holes to get centre coords
888  for(unsigned ihole=0;ihole<nhole;ihole++)
889  {
890  Hole_centre[ihole].resize(3);
891 
892  // Read the centre value
893  Hole_centre[ihole][0]=triangle_data.holelist[count_coords];
894  Hole_centre[ihole][1]=triangle_data.holelist[count_coords+1];
895  Hole_centre[ihole][2]=triangle_data.holelist[count_coords+2];
896 
897  // Increment coords counter
898  count_coords+=3;
899  }
900  }
901  else
902  {
903  Hole_centre.resize(0);
904  }*/
905 
906 
907  // Set number of boundaries
908  if(n_bound > 0)
909  {
910  this->set_nboundary(n_bound);
911  }
912 
913  // Create the elements
914  unsigned counter=0;
915  for(unsigned e=0;e<n_element;e++)
916  {
918  unsigned global_node_number = Global_node[counter];
919  if(done[global_node_number-1]==false)
920  // ... -1 because node number begins at 1 in tetgen
921  {
922  //If the node is on a boundary, construct a boundary node
923  if((boundary_markers_flag==1) &&
924  (bound[global_node_number-1] > 0))
925  {
926  //Construct the boundary ndoe
927  Node_pt[global_node_number-1] =
929  //Add to the boundary lookup scheme
930  add_boundary_node(bound[global_node_number-1]-1,
931  Node_pt[global_node_number-1]);
932  }
933  //Otherwise just construct a normal node
934  else
935  {
936  Node_pt[global_node_number-1] =
938  }
939 
940  done[global_node_number-1]=true;
941  Node_pt[global_node_number-1]->x(0)=x_node[global_node_number-1];
942  Node_pt[global_node_number-1]->x(1)=y_node[global_node_number-1];
943  Node_pt[global_node_number-1]->x(2)=z_node[global_node_number-1];
944  }
945  //Otherwise just copy the node numbr accross
946  else
947  {
948  finite_element_pt(e)->node_pt(3) = Node_pt[global_node_number-1];
949  }
950  counter++;
951 
952  //Loop over the other nodes
953  for(unsigned j=0;j<(n_local_node-1);j++)
954  {
955  global_node_number=Global_node[counter];
956  if(done[global_node_number-1]==false)
957  // ... -1 because node number begins at 1 in tetgen
958  {
959  //If we're on a boundary
960  if((boundary_markers_flag==1) &&
961  (bound[global_node_number-1] > 0))
962  {
963  //Construct the boundary ndoe
964  Node_pt[global_node_number-1] =
966  //Add to the boundary lookup scheme
967  add_boundary_node(bound[global_node_number-1]-1,
968  Node_pt[global_node_number-1]);
969  }
970  else
971  {
972  Node_pt[global_node_number-1]=finite_element_pt(e)->construct_node(j);
973  }
974  done[global_node_number-1]=true;
975  Node_pt[global_node_number-1]->x(0)=x_node[global_node_number-1];
976  Node_pt[global_node_number-1]->x(1)=y_node[global_node_number-1];
977  Node_pt[global_node_number-1]->x(2)=z_node[global_node_number-1];
978  }
979  //Otherwise copy the pointer over
980  else
981  {
982  finite_element_pt(e)->node_pt(j) = Node_pt[global_node_number-1];
983  }
984  counter++;
985  }
986  }
987 
988 
989  // Resize the "matrix" that stores the boundary id for each
990  // face in each element.
991  Face_boundary.resize(n_element);
992  Face_index.resize(n_element);
993  Edge_index.resize(n_element);
994 
995 
996  //0-based index scheme used to construct a global lookup for
997  //each face that will be used to uniquely construct interior
998  //face nodes.
999  //The faces that lie on boundaries will have already been allocated so
1000  //we can start from there
1001  Nglobal_face = n_face;
1002 
1003  // Storage for the edges associated with each node.
1004  // Nodes are indexed using Tetgen's 1-based scheme, which is why
1005  // there is a +1 here
1006  Vector<std::set<unsigned> > node_on_edges(n_node+1);
1007 
1008  //0-based index scheme used to construct a global lookup for each
1009  //edge that will be used to uniquely construct interior edge nodes
1010  Nglobal_edge = 0;
1011 
1012  // Conversion from the edge numbers to the nodes at the end
1013  // of each each edge
1014  const unsigned first_local_edge_node[6] = {0,0,0,1,2,1};
1015  const unsigned second_local_edge_node[6] = {1,2,3,2,3,3};
1016 
1017  // Loop over the elements
1018  for(unsigned e=0;e<n_element;e++)
1019  {
1020  // Each element has four faces
1021  Face_boundary[e].resize(4);
1022  Face_index[e].resize(4);
1023  // By default each face is NOT on a boundary
1024  for(unsigned i=0;i<4;i++)
1025  {
1026  Face_boundary[e][i]=0;
1027  }
1028 
1029  Edge_index[e].resize(6);
1030 
1031  // Read out the global node numbers of the nodes from
1032  // tetgen's 1-based numbering.
1033  // The offset is to match the offset used above
1034  const unsigned element_offset = e*n_local_node;
1035  unsigned glob_num[4] = {0,0,0,0};
1036  for(unsigned i=0;i<4;++i)
1037  {
1038  glob_num[i] = Global_node[element_offset + ((i+1)%4)];
1039  }
1040 
1041  // Now we know the global node numbers of the elements' four nodes
1042  // in tetgen's 1-based numbering.
1043 
1044  //Determine whether any of the element faces have already been allocated an
1045  //index. This may be because they are located on boundaries, or
1046  //have already been visited
1047  //The global index for the i-th face will be stored in Face_index[i]
1048 
1049  //Loop over the local faces in the element
1050  for(unsigned i=0;i<4;++i)
1051  {
1052  //On the i-th face, our numbering convention is such that
1053  //it is the (3-i)th node of the element that is omitted
1054  const unsigned omitted_node = 3-i;
1055 
1056  //Start with the set of global faces associated with the i-th node
1057  //which is always on the i-th face
1058  std::set<unsigned> input = node_on_faces[glob_num[i]];
1059 
1060  //If there is no data yet, not point doing intersections
1061  //the face has not been set up
1062  if(input.size() > 0)
1063  {
1064  //Loop over the other nodes on the face
1065  unsigned local_node = i;
1066  for(unsigned i2=0;i2<3;++i2)
1067  {
1068  //Create the next node index (mod 4)
1069  local_node = (local_node+1)%4;
1070  //Only take the intersection of nodes on the face
1071  //i.e. not 3-i
1072  if(local_node != omitted_node)
1073  {
1074  //Local storage for the intersection
1075  std::set<unsigned> intersection_result;
1076  //Calculate the intersection of the current input
1077  //and next node
1078  std::set_intersection(input.begin(),input.end(),
1079  node_on_faces[glob_num[local_node]].begin(),
1080  node_on_faces[glob_num[local_node]].end(),
1081  std::insert_iterator<std::set<unsigned> >(
1082  intersection_result,
1083  intersection_result.begin()));
1084  //Let the result be the next input
1085  input = intersection_result;
1086  }
1087  }
1088  }
1089 
1090  //If the nodes share more than one global face index, then
1091  //we have a problem
1092  if(input.size() > 1)
1093  {
1094  throw OomphLibError(
1095  "Nodes in scaffold mesh share more than one global face",
1096  OOMPH_CURRENT_FUNCTION,
1097  OOMPH_EXCEPTION_LOCATION);
1098  }
1099 
1100  //If the element's face is not already allocated, the intersection
1101  //will be empty
1102  if(input.size()==0)
1103  {
1104  //Allocate the next global index
1105  Face_index[e][i] = Nglobal_face;
1106  //Associate the new face index with the nodes
1107  for(unsigned i2=0;i2<4;++i2)
1108  {
1109  //Don't add the omitted node
1110  if(i2 != omitted_node)
1111  {
1112  node_on_faces[glob_num[i2]].insert(Nglobal_face);
1113  }
1114  }
1115  //Increment the global face index
1116  ++Nglobal_face;
1117  }
1118  //Otherwise we already have a face
1119  else if(input.size()==1)
1120  {
1121  const unsigned global_face_index = *input.begin();
1122  //Set the face index
1123  Face_index[e][i] = global_face_index;
1124  //Allocate the boundary index, if it's a boundary
1125  if(global_face_index < n_face)
1126  {
1127  Face_boundary[e][i] = face_boundary[global_face_index];
1128  //Add the nodes to the boundary look-up scheme in
1129  //oomph-lib (0-based) index
1130  for(unsigned i2=0;i2<4;++i2)
1131  {
1132  //Don't add the omitted node
1133  if(i2 != omitted_node)
1134  {
1135  add_boundary_node(face_boundary[global_face_index]-1,
1136  Node_pt[glob_num[i2]-1]);
1137  }
1138  }
1139  }
1140  }
1141  } //end of loop over the faces
1142 
1143 
1144  //Loop over the element edges and assign global edge numbers
1145  for(unsigned i=0;i<6;++i)
1146  {
1147  //Storage for the result of the intersection
1148  std::vector<unsigned> local_edge_index;
1149 
1150  const unsigned first_global_num = glob_num[first_local_edge_node[i]];
1151  const unsigned second_global_num = glob_num[second_local_edge_node[i]];
1152 
1153  //Find the common global edge index for the nodes on
1154  //the i-th element edge
1155  std::set_intersection(node_on_edges[first_global_num].begin(),
1156  node_on_edges[first_global_num].end(),
1157  node_on_edges[second_global_num].begin(),
1158  node_on_edges[second_global_num].end(),
1159  std::insert_iterator<std::vector<unsigned> >(
1160  local_edge_index,local_edge_index.begin()));
1161 
1162  //If the nodes share more than one global edge index, then
1163  //we have a problem
1164  if(local_edge_index.size() > 1)
1165  {
1166  throw OomphLibError(
1167  "Nodes in scaffold mesh share more than one global edge",
1168  OOMPH_CURRENT_FUNCTION,
1169  OOMPH_EXCEPTION_LOCATION);
1170  }
1171 
1172  //If the element's edge is not already allocated, the intersection
1173  //will be empty
1174  if(local_edge_index.size()==0)
1175  {
1176  //Allocate the next global index
1177  Edge_index[e][i] = Nglobal_edge;
1178  //Associate the new edge index with the nodes
1179  node_on_edges[first_global_num].insert(Nglobal_edge);
1180  node_on_edges[second_global_num].insert(Nglobal_edge);
1181  //Increment the global edge index
1182  ++Nglobal_edge;
1183  }
1184  //Otherwise we already have an edge
1185  else if(local_edge_index.size()==1)
1186  {
1187  //Set the edge index from the result of the intersection
1188  Edge_index[e][i] = local_edge_index[0];
1189  }
1190  }
1191 
1192  } // end for e
1193 
1194 
1195  //Now determine whether any edges lie on boundaries by using the
1196  //face boundary scheme and
1197 
1198  //Resize the storage
1199  Edge_boundary.resize(Nglobal_edge,false);
1200 
1201  //Now loop over all the boundary faces and mark that all edges
1202  //must also lie on the bounadry
1203  for(unsigned i=0;i<n_face;++i)
1204  {
1205  {
1206  //Storage for the result of the intersection
1207  std::vector<unsigned> local_edge_index;
1208 
1209  //Find the common global edge index for first edge of the face
1210  std::set_intersection(node_on_edges[first_node[i]].begin(),
1211  node_on_edges[first_node[i]].end(),
1212  node_on_edges[second_node[i]].begin(),
1213  node_on_edges[second_node[i]].end(),
1214  std::insert_iterator<std::vector<unsigned> >(
1215  local_edge_index,local_edge_index.begin()));
1216 
1217  //If the nodes do not share exactly one global edge index, then
1218  //we have a problem
1219  if(local_edge_index.size() != 1)
1220  {
1221  throw OomphLibError(
1222  "Nodes in scaffold mesh face do not share exactly one global edge",
1223  OOMPH_CURRENT_FUNCTION,
1224  OOMPH_EXCEPTION_LOCATION);
1225  }
1226  else
1227  {
1228  Edge_boundary[local_edge_index[0]] = true;
1229  }
1230  }
1231 
1232  {
1233  //Storage for the result of the intersection
1234  std::vector<unsigned> local_edge_index;
1235 
1236  //Find the common global edge index for second edge of the face
1237  std::set_intersection(node_on_edges[second_node[i]].begin(),
1238  node_on_edges[second_node[i]].end(),
1239  node_on_edges[third_node[i]].begin(),
1240  node_on_edges[third_node[i]].end(),
1241  std::insert_iterator<std::vector<unsigned> >(
1242  local_edge_index,local_edge_index.begin()));
1243 
1244  //If the nodes do not share exactly one global edge index, then
1245  //we have a problem
1246  if(local_edge_index.size() != 1)
1247  {
1248  throw OomphLibError(
1249  "Nodes in scaffold mesh face do not share exactly one global edge",
1250  OOMPH_CURRENT_FUNCTION,
1251  OOMPH_EXCEPTION_LOCATION);
1252  }
1253  else
1254  {
1255  Edge_boundary[local_edge_index[0]] = true;
1256  }
1257  }
1258 
1259  {
1260  //Storage for the result of the intersection
1261  std::vector<unsigned> local_edge_index;
1262 
1263  //Find the common global edge index for third edge of the face
1264  std::set_intersection(node_on_edges[first_node[i]].begin(),
1265  node_on_edges[first_node[i]].end(),
1266  node_on_edges[third_node[i]].begin(),
1267  node_on_edges[third_node[i]].end(),
1268  std::insert_iterator<std::vector<unsigned> >(
1269  local_edge_index,local_edge_index.begin()));
1270 
1271  //If the nodes do not share exactly one global edge index, then
1272  //we have a problem
1273  if(local_edge_index.size() != 1)
1274  {
1275  throw OomphLibError(
1276  "Nodes in scaffold mesh face do not share exactly one global edge",
1277  OOMPH_CURRENT_FUNCTION,
1278  OOMPH_EXCEPTION_LOCATION);
1279  }
1280  else
1281  {
1282  Edge_boundary[local_edge_index[0]] = true;
1283  }
1284  }
1285  }
1286 
1287 
1288 } //end of constructor
1289 
1290 
1291 }
1292 
1293 
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< Node * > Node_pt
Vector of pointers to nodes.
Definition: mesh.h:194
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 Nglobal_face
Storage for the number of global faces.
Vector< Vector< unsigned > > Face_boundary
Vector of vectors containing the boundary ids of the elements' faces.
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< 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
e
Definition: cfortran.h:575
TetgenScaffoldMesh()
Empty constructor.
unsigned Nglobal_edge
Storage for the number of global edges.
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
Vector< double > Element_attribute
Vector of double attributes for each element.
Vector< unsigned > Global_node
Storage for global node numbers listed element-by-element.
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
unsigned face_boundary(const unsigned &e, const unsigned &i) const
Return the boundary id of the i-th face in the e-th element: This is zero-based as in tetgen...
Vector< Vector< unsigned > > Face_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.
std::vector< bool > Edge_boundary
Vector of booleans to indicate whether a global edge lies on a boundary.