tetgen_mesh.template.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 #ifndef OOMPH_TETGEN_MESH_TEMPLATE_CC
27 #define OOMPH_TETGEN_MESH_TEMPLATE_CC
28 
29 
30 #include <algorithm>
31 
32 #include "tetgen_mesh.template.h"
33 #include "../generic/Telements.h"
34 #include "../generic/map_matrix.h"
35 
36 
37 namespace oomph
38 {
39  /// ////////////////////////////////////////////////////////////////////
40  /// ////////////////////////////////////////////////////////////////////
41  /// ////////////////////////////////////////////////////////////////////
42 
43 
44  //========================================================================
45  /// Build unstructured tet mesh based on output from scaffold
46  //========================================================================
47  template<class ELEMENT>
48  void TetgenMesh<ELEMENT>::build_from_scaffold(TimeStepper* time_stepper_pt,
49  const bool& use_attributes)
50  {
51  // Mesh can only be built with 3D Telements.
52  MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
53 
54  // Create space for elements
55  unsigned nelem = Tmp_mesh_pt->nelement();
56  Element_pt.resize(nelem);
57 
58  // Create space for nodes
59  unsigned nnode_scaffold = Tmp_mesh_pt->nnode();
60  Node_pt.resize(nnode_scaffold);
61 
62  // Set number of boundaries
63  unsigned nbound = Tmp_mesh_pt->nboundary();
64  set_nboundary(nbound);
65  Boundary_element_pt.resize(nbound);
66  Face_index_at_boundary.resize(nbound);
67 
68  // If we have different regions, then resize the region
69  // information
70  if (use_attributes)
71  {
72  Boundary_region_element_pt.resize(nbound);
73  Face_index_region_at_boundary.resize(nbound);
74  }
75 
76  // Build elements
77  for (unsigned e = 0; e < nelem; e++)
78  {
79  Element_pt[e] = new ELEMENT;
80  }
81 
82  // Number of nodes per element
83  unsigned nnod_el = Tmp_mesh_pt->finite_element_pt(0)->nnode();
84 
85  // Setup map to check the (pseudo-)global node number
86  // Nodes whose number is zero haven't been copied across
87  // into the mesh yet.
88  std::map<Node*, unsigned> global_number;
89  unsigned global_count = 0;
90 
91  // Map of element attribute pairs
92  std::map<double, Vector<FiniteElement*>> element_attribute_map;
93 
94  // Loop over elements in scaffold mesh, visit their nodes
95  for (unsigned e = 0; e < nelem; e++)
96  {
97  // Loop over all nodes in element
98  for (unsigned j = 0; j < nnod_el; j++)
99  {
100  // Pointer to node in the scaffold mesh
101  Node* scaffold_node_pt = Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
102 
103  // Get the (pseudo-)global node number in scaffold mesh
104  // (It's zero [=default] if not visited this one yet)
105  unsigned j_global = global_number[scaffold_node_pt];
106 
107  // Haven't done this one yet
108  if (j_global == 0)
109  {
110  // Get pointer to set of mesh boundaries that this
111  // scaffold node occupies; NULL if the node is not on any boundary
112  std::set<unsigned>* boundaries_pt;
113  scaffold_node_pt->get_boundaries_pt(boundaries_pt);
114 
115  // Is it on boundaries?
116  if (boundaries_pt != 0)
117  {
118  // Create new boundary node
119  Node* new_node_pt =
120  finite_element_pt(e)->construct_boundary_node(j, time_stepper_pt);
121 
122  // Give it a number (not necessarily the global node
123  // number in the scaffold mesh -- we just need something
124  // to keep track...)
125  global_count++;
126  global_number[scaffold_node_pt] = global_count;
127 
128  // Add to boundaries
129  for (std::set<unsigned>::iterator it = boundaries_pt->begin();
130  it != boundaries_pt->end();
131  ++it)
132  {
133  add_boundary_node(*it, new_node_pt);
134  }
135  }
136  // Build normal node
137  else
138  {
139  // Create new normal node
140  finite_element_pt(e)->construct_node(j, time_stepper_pt);
141 
142  // Give it a number (not necessarily the global node
143  // number in the scaffold mesh -- we just need something
144  // to keep track...)
145  global_count++;
146  global_number[scaffold_node_pt] = global_count;
147  }
148 
149  // Copy new node, created using the NEW element's construct_node
150  // function into global storage, using the same global
151  // node number that we've just associated with the
152  // corresponding node in the scaffold mesh
153  Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
154 
155  // Assign coordinates
156  Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
157  Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
158  Node_pt[global_count - 1]->x(2) = scaffold_node_pt->x(2);
159  }
160  // This one has already been done: Copy across
161  else
162  {
163  finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
164  }
165  }
166 
167  // Store the attributes in the map
168  if (use_attributes)
169  {
170  element_attribute_map[Tmp_mesh_pt->element_attribute(e)].push_back(
171  finite_element_pt(e));
172  }
173  }
174 
175  // Now let's construct lists
176  // Find the number of attributes
177  if (use_attributes)
178  {
179  unsigned n_attribute = element_attribute_map.size();
180  // There are n_attribute different regions
181  Region_element_pt.resize(n_attribute);
182  Region_attribute.resize(n_attribute);
183  // Copy the vectors in the map over to our internal storage
184  unsigned count = 0;
185  for (std::map<double, Vector<FiniteElement*>>::iterator it =
186  element_attribute_map.begin();
187  it != element_attribute_map.end();
188  ++it)
189  {
190  Region_attribute[count] = it->first;
191  Region_element_pt[count] = it->second;
192  ++count;
193  }
194  }
195 
196  // At this point we've created all the elements and
197  // created their vertex nodes. Now we need to create
198  // the additional (midside and internal) nodes!
199  unsigned boundary_id;
200 
201  // Get number of nodes along element edge and dimension of element (3)
202  // from first element
203  unsigned n_node_1d = finite_element_pt(0)->nnode_1d();
204 
205  // At the moment we're only able to deal with nnode_1d=2 or 3.
206  if ((n_node_1d != 2) && (n_node_1d != 3))
207  {
208  std::ostringstream error_message;
209  error_message << "Mesh generation from tetgen currently only works\n";
210  error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
211  error_message << "for nnode_1d=" << n_node_1d << std::endl;
212 
213  throw OomphLibError(
214  error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
215  }
216 
217  // Spatial dimension of element = number of local coordinates
218  unsigned dim = finite_element_pt(0)->dim();
219 
220  // Storage for the local coordinate of the new node
221  Vector<double> s(dim);
222 
223  // Get number of nodes in the element from first element
224  unsigned n_node = finite_element_pt(0)->nnode();
225 
226  // Storage for each global edge of the mesh
227  unsigned n_global_edge = Tmp_mesh_pt->nglobal_edge();
228  Vector<Vector<Node*>> nodes_on_global_edge(n_global_edge);
229 
230  // Storage for each global face of the mesh
231  unsigned n_global_face = Tmp_mesh_pt->nglobal_face();
232  Vector<Vector<Node*>> nodes_on_global_face(n_global_face);
233 
234  // Map storing the mid-side of an edge; edge identified by
235  // pointers to vertex nodes in scaffold mesh
236  // MapMatrix<Node*,Node*> central_edge_node_pt;
237  // Node* edge_node1_pt=0;
238  // Node* edge_node2_pt=0;
239 
240  // Map storing the mid point of a face; face identified by
241  // set of pointers to vertex nodes in scaffold mesh
242  // std::map<std::set<Node*>,Node*> central_face_node_pt;
243  // std::set<Node*> face_nodes_pt;
244 
245  // Mapping of Tetgen faces to face nodes in the enriched element
246  unsigned face_map[4] = {1, 0, 2, 3};
247 
248  // Storage for the faces shared by the edges
249  const unsigned faces_on_edge[6][2] = {
250  {0, 1}, {0, 2}, {1, 2}, {0, 3}, {2, 3}, {1, 3}};
251 
252  // Loop over all elements
253  for (unsigned e = 0; e < nelem; e++)
254  {
255  // Cache pointers to the elements
256  FiniteElement* const elem_pt = this->finite_element_pt(e);
257  FiniteElement* const tmp_elem_pt = Tmp_mesh_pt->finite_element_pt(e);
258 
259  // The number of edge nodes is 4 + 6*(n_node1d-2)
260  unsigned n_edge_node = 4 + 6 * (n_node_1d - 2);
261 
262  // Now loop over the edge nodes
263  // assuming that the numbering scheme is the same as the triangles
264  // which puts edge nodes in ascending order.
265  // We don't have any higher than quadratic at the moment, so it's
266  // a bit academic really.
267 
268  // Only bother if we actually have some edge nodes
269  if (n_edge_node > 4)
270  {
271  // Start from node number 4
272  unsigned n = 4;
273 
274  // Loop over the edges
275  for (unsigned j = 0; j < 6; ++j)
276  {
277  // Find the global edge index
278  unsigned edge_index = Tmp_mesh_pt->edge_index(e, j);
279 
280  // Use the intersection of the appropriate faces to determine
281  // whether the boundaries on which an edge lies
282  std::set<unsigned> edge_boundaries;
283  for (unsigned i = 0; i < 2; ++i)
284  {
285  unsigned face_boundary_id =
286  Tmp_mesh_pt->face_boundary(e, faces_on_edge[j][i]);
287  if (face_boundary_id > 0)
288  {
289  edge_boundaries.insert(face_boundary_id);
290  }
291  }
292 
293  // If the nodes on the edge have not been allocated, construct them
294  if (nodes_on_global_edge[edge_index].size() == 0)
295  {
296  // Now loop over the nodes on the edge
297  for (unsigned j2 = 0; j2 < n_node_1d - 2; ++j2)
298  {
299  // Storage for the new node
300  Node* new_node_pt = 0;
301 
302  // If the edge is on a boundary, determined from the
303  // scaffold mesh, construct a boundary node
304  // The use of the scaffold mesh's edge_boundary data structure
305  // ensures that a boundary node is created, even if the faces of
306  // the current element do not lie on boundaries
307  if (Tmp_mesh_pt->edge_boundary(edge_index) == true)
308  {
309  new_node_pt =
310  elem_pt->construct_boundary_node(n, time_stepper_pt);
311  // Add it to the boundaries in the set,
312  // remembering to subtract one to get to the oomph-lib numbering
313  // scheme
314  for (std::set<unsigned>::iterator it = edge_boundaries.begin();
315  it != edge_boundaries.end();
316  ++it)
317  {
318  this->add_boundary_node((*it) - 1, new_node_pt);
319  }
320  }
321  // Otherwise construct a normal node
322  else
323  {
324  new_node_pt = elem_pt->construct_node(n, time_stepper_pt);
325  }
326 
327  // Find the local coordinates of the node
328  elem_pt->local_coordinate_of_node(n, s);
329 
330  // Find the coordinates of the new node from the exiting and
331  // fully-functional element in the scaffold mesh
332  for (unsigned i = 0; i < dim; ++i)
333  {
334  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
335  }
336 
337  // Add the newly created node to the global node list
338  Node_pt.push_back(new_node_pt);
339  // Add to the edge index
340  nodes_on_global_edge[edge_index].push_back(new_node_pt);
341  // Increment the local node number
342  ++n;
343  } // end of loop over edge nodes
344  }
345  // Otherwise just set the pointers (orientation assumed the same)
346  else
347  {
348  for (unsigned j2 = 0; j2 < n_node_1d - 2; ++j2)
349  {
350  elem_pt->node_pt(n) = nodes_on_global_edge[edge_index][j2];
351  // It is possible that the edge may be on additional boundaries
352  // through another element
353  // So add again (note that this function will not add to
354  // boundaries twice)
355  for (std::set<unsigned>::iterator it = edge_boundaries.begin();
356  it != edge_boundaries.end();
357  ++it)
358  {
359  this->add_boundary_node((*it) - 1, elem_pt->node_pt(n));
360  }
361  ++n;
362  }
363  }
364  } // End of loop over edges
365 
366  // Deal with enriched elements
367  if (n_node == 15)
368  {
369  // Now loop over the faces
370  for (unsigned j = 0; j < 4; ++j)
371  {
372  // Find the boundary id of the face (need to map between node
373  // numbers and the face)
374  boundary_id = Tmp_mesh_pt->face_boundary(e, face_map[j]);
375 
376  // Find the global face index (need to map between node numbers and
377  // the face)
378  unsigned face_index = Tmp_mesh_pt->face_index(e, face_map[j]);
379 
380  // If the nodes on the face have not been allocated
381  if (nodes_on_global_face[face_index].size() == 0)
382  {
383  // Storage for the new node
384  Node* new_node_pt = 0;
385 
386  // If the face is on a boundary, construct a boundary node
387  if (boundary_id > 0)
388  {
389  new_node_pt =
390  elem_pt->construct_boundary_node(n, time_stepper_pt);
391  // Add it to the boundary
392  this->add_boundary_node(boundary_id - 1, new_node_pt);
393  }
394  // Otherwise construct a normal node
395  else
396  {
397  new_node_pt = elem_pt->construct_node(n, time_stepper_pt);
398  }
399 
400  // Find the local coordinates of the node
401  elem_pt->local_coordinate_of_node(n, s);
402 
403  // Find the coordinates of the new node from the exiting and
404  // fully-functional element in the scaffold mesh
405  for (unsigned i = 0; i < dim; ++i)
406  {
407  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
408  }
409 
410  // Add the newly created node to the global node list
411  Node_pt.push_back(new_node_pt);
412  // Add to the face index
413  nodes_on_global_face[face_index].push_back(new_node_pt);
414  // Increment the local node number
415  ++n;
416  }
417  // Otherwise just set the single node from the face element
418  else
419  {
420  elem_pt->node_pt(n) = nodes_on_global_face[face_index][0];
421  ++n;
422  }
423  } // end of loop over faces
424 
425  // Construct the element's central node, which is not on a boundary
426  {
427  Node* new_node_pt =
428  finite_element_pt(e)->construct_node(n, time_stepper_pt);
429  Node_pt.push_back(new_node_pt);
430 
431  // Find the local coordinates of the node
432  elem_pt->local_coordinate_of_node(n, s);
433 
434  // Find the coordinates of the new node from the existing
435  // and fully-functional element in the scaffold mesh
436  for (unsigned i = 0; i < dim; i++)
437  {
438  new_node_pt->x(i) = tmp_elem_pt->interpolated_x(s, i);
439  }
440  }
441  } // End of enriched case
442 
443  } // End of case for edge nodes
444 
445 
446  // Now loop over the faces to setup the information about elements
447  // adjacent to the boundary
448  for (unsigned j = 0; j < 4; ++j)
449  {
450  // Find the boundary id of the face
451  boundary_id = Tmp_mesh_pt->face_boundary(e, j);
452 
453  if (boundary_id > 0)
454  {
455  Boundary_element_pt[boundary_id - 1].push_back(elem_pt);
456  // Need to put a shift in here because of an inconsistent naming
457  // convention between tetgen and our faces
458  // Tetgen Face 0 is our Face 3
459  // Tetgen Face 1 is our Face 2
460  // Tetgen Face 2 is our Face 1
461  // Tetgen Face 3 is our Face 0
462  Face_index_at_boundary[boundary_id - 1].push_back(3 - j);
463 
464  // If using regions set up the boundary information
465  if (use_attributes)
466  {
467  // Element adjacent to boundary
468  Boundary_region_element_pt[boundary_id - 1]
469  [static_cast<unsigned>(
470  Tmp_mesh_pt->element_attribute(e))]
471  .push_back(elem_pt);
472  // Need to put a shift in here because of an inconsistent naming
473  // convention between triangle and face elements
474  Face_index_region_at_boundary[boundary_id - 1]
475  [static_cast<unsigned>(
476  Tmp_mesh_pt->element_attribute(e))]
477  .push_back(3 - j);
478  }
479  }
480  } // End of loop over faces
481 
482 
483  // Lookup scheme has now been setup
484  Lookup_for_elements_next_boundary_is_setup = true;
485 
486 
487  // /*
488 
489  // //For standard quadratic elements all nodes are edge nodes
490  // unsigned n_vertex_and_edge_node = n_node;
491  // //If we have enriched elements, there are only 10 vertex and edge
492  // nodes if(n_node==15)
493  // {
494  // //There are only 10 vertex and edge nodes
495  // n_vertex_and_edge_node = 10;
496  // }
497 
498  // // Loop over the new (edge) nodes in the element and create them.
499  // for(unsigned j=4;j<n_vertex_and_edge_node;j++)
500  // {
501 
502  // // Figure out edge nodes
503  // switch(j)
504  // {
505 
506  // // Node 4 is between nodes 0 and 1
507  // case 4:
508 
509  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
510  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
511  // break;
512 
513 
514  // // Node 5 is between nodes 0 and 2
515  // case 5:
516 
517  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
518  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
519  // break;
520 
521  // // Node 6 is between nodes 0 and 3
522  // case 6:
523 
524  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(0);
525  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
526  // break;
527 
528  // // Node 7 is between nodes 1 and 2
529  // case 7:
530 
531  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
532  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
533  // break;
534 
535  // // Node 8 is between nodes 2 and 3
536  // case 8:
537 
538  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(2);
539  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
540  // break;
541 
542  // // Node 9 is between nodes 1 and 3
543  // case 9:
544 
545  // edge_node1_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(1);
546  // edge_node2_pt=Tmp_mesh_pt->finite_element_pt(e)->node_pt(3);
547  // break;
548 
549  // break;
550 
551  // //Error
552  // default:
553 
554  // throw OomphLibError("More than ten edge nodes in Tet Element",
555  // OOMPH_CURRENT_FUNCTION,
556  // OOMPH_EXCEPTION_LOCATION);
557  // }
558 
559 
560  // // Do we need a boundary node?
561  // bool need_boundary_node=false;
562 
563  // // hierher At some point fine tune via intersection of boundary
564  // sets if
565  // (edge_node1_pt->is_on_boundary()&&edge_node2_pt->is_on_boundary())
566  // {
567  // need_boundary_node=true;
568  // }
569 
570  // // Do we need a new node?
571  // if (central_edge_node_pt(edge_node1_pt,edge_node2_pt)==0)
572  // {
573  // Node* new_node_pt=0;
574 
575  // // Create new boundary node
576  // if (need_boundary_node)
577  // {
578  // new_node_pt=finite_element_pt(e)->
579  // construct_boundary_node(j,time_stepper_pt);
580  // }
581  // // Create new normal node
582  // else
583  // {
584  // new_node_pt=finite_element_pt(e)->
585  // construct_node(j,time_stepper_pt);
586  // }
587  // Node_pt.push_back(new_node_pt);
588 
589  // // Now indicate existence of newly created mideside node in map
590  // central_edge_node_pt(edge_node1_pt,edge_node2_pt)=new_node_pt;
591  // central_edge_node_pt(edge_node2_pt,edge_node1_pt)=new_node_pt;
592 
593  // // What are the node's local coordinates?
594  // finite_element_pt(e)->local_coordinate_of_node(j,s);
595 
596  // // Find the coordinates of the new node from the existing
597  // // and fully-functional element in the scaffold mesh
598  // for(unsigned i=0;i<dim;i++)
599  // {
600  // new_node_pt->x(i)=
601  // Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
602  // }
603  // }
604  // else
605  // {
606  // // Set pointer to the existing node
607  // finite_element_pt(e)->node_pt(j)=
608  // central_edge_node_pt(edge_node1_pt,edge_node2_pt);
609  // }
610 
611  // } // end of loop over new nodes
612 
613  // //Need to sort out the face nodes
614  // if(n_node==15)
615  // {
616  // // Loop over the new (face) nodes in the element and create them.
617  // for(unsigned j=10;j<14;j++)
618  // {
619  // //Empty the set of face nodes
620  // face_nodes_pt.clear();
621  // // Figure out face nodes
622  // switch(j)
623  // {
624 
625  // // Node 10 is between nodes 0 and 1 and 3
626  // case 10:
627 
628  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
629  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
630  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
631  // break;
632 
633  // // Node 11 is between nodes 0 and 1 and 2
634  // case 11:
635 
636  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
637  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
638  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
639  // break;
640 
641  // // Node 12 is between nodes 0 and 2 and 3
642  // case 12:
643 
644  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(0));
645  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
646  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
647  // break;
648 
649  // // Node 13 is between nodes 1 and 2 and 3
650  // case 13:
651 
652  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(1));
653  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(2));
654  // face_nodes_pt.insert(Tmp_mesh_pt->finite_element_pt(e)->node_pt(3));
655  // break;
656 
657 
658  // //Error
659  // default:
660 
661  // throw OomphLibError("More than four face nodes in Tet Element",
662  // OOMPH_CURRENT_FUNCTION,
663  // OOMPH_EXCEPTION_LOCATION);
664  // }
665 
666  // // Do we need a boundary node?
667  // bool need_boundary_node=false;
668 
669  // //Work it out from the face boundary
670  // boundary_id = Tmp_mesh_pt->face_boundary(e,face_map[j-10]);
671  // //If it's non-zero then we do need to create a boundary node
672  // if(boundary_id!=0) {need_boundary_node=true;}
673 
674  // // Do we need a new node?
675  // if (central_face_node_pt[face_nodes_pt]==0)
676  // {
677  // Node* new_node_pt=0;
678 
679  // // Create a new boundary node
680  // if (need_boundary_node)
681  // {
682  // new_node_pt=finite_element_pt(e)->
683  // construct_boundary_node(j,time_stepper_pt);
684  // //Add it to the boundary
685  // add_boundary_node(boundary_id-1,new_node_pt);
686  // }
687  // // Create new normal node
688  // else
689  // {
690  // new_node_pt=finite_element_pt(e)->
691  // construct_node(j,time_stepper_pt);
692  // }
693  // Node_pt.push_back(new_node_pt);
694 
695  // // Now indicate existence of newly created mideside node in map
696  // central_face_node_pt[face_nodes_pt]=new_node_pt;
697 
698  // // What are the node's local coordinates?
699  // finite_element_pt(e)->local_coordinate_of_node(j,s);
700 
701  // // Find the coordinates of the new node from the existing
702  // // and fully-functional element in the scaffold mesh
703  // for(unsigned i=0;i<dim;i++)
704  // {
705  // new_node_pt->x(i)=
706  // Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
707  // }
708  // }
709  // else
710  // {
711  // // Set pointer to the existing node
712  // finite_element_pt(e)->node_pt(j)=
713  // central_face_node_pt[face_nodes_pt];
714  // }
715  // } //End of loop over face nodes
716 
717  // //Construct the element's central node, which is not on a boundary
718  // {
719  // Node* new_node_pt=
720  // finite_element_pt(e)->construct_node(14,time_stepper_pt);
721  // Node_pt.push_back(new_node_pt);
722 
723  // // What are the node's local coordinates?
724  // finite_element_pt(e)->local_coordinate_of_node(14,s);
725  // // Find the coordinates of the new node from the existing
726  // // and fully-functional element in the scaffold mesh
727  // for(unsigned i=0;i<dim;i++)
728  // {
729  // new_node_pt->x(i)=
730  // Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s,i);
731  // }
732  // }
733  // } //End of enriched case
734 
735  // } //end of loop over elements
736 
737 
738  // //Boundary conditions
739 
740  // // Matrix map to check if a node has already been added to
741  // // the boundary number b
742  // MapMatrixMixed<Node*,unsigned,bool> bound_node_done;
743 
744  // // Loop over elements
745  // for (unsigned e=0;e<nelem;e++)
746  // {
747  // // Loop over new local nodes
748  // for(unsigned j=4;j<n_node;j++)
749  // {
750  // // Pointer to the element's local node
751  // Node* loc_node_pt=finite_element_pt(e)->node_pt(j);
752 
753  // // This will have to be changed for higher-order elements
754  // //=======================================================
755 
756  // // These are the face nodes on the element's face 0:
757  // if ( (j==4) || (j==5) || (j==7) )
758  // {
759  // boundary_id=Tmp_mesh_pt->face_boundary(e,0);
760  // if (boundary_id!=0)
761  // {
762  // if (!bound_node_done(loc_node_pt,boundary_id-1))
763  // {
764  // add_boundary_node(boundary_id-1,loc_node_pt);
765  // bound_node_done(loc_node_pt,boundary_id-1)=true;
766  // }
767  // }
768  // }
769 
770 
771  // // These are the face nodes on the element's face 1:
772  // if ( (j==4) || (j==6) || (j==9) )
773  // {
774  // boundary_id=Tmp_mesh_pt->face_boundary(e,1);
775  // if (boundary_id!=0)
776  // {
777  // if (!bound_node_done(loc_node_pt,boundary_id-1))
778  // {
779  // add_boundary_node(boundary_id-1,loc_node_pt);
780  // bound_node_done(loc_node_pt,boundary_id-1)=true;
781  // }
782  // }
783  // }
784 
785  // // These are the face nodes on the element's face 2:
786  // if ( (j==5) || (j==6) || (j==8) )
787  // {
788  // boundary_id=Tmp_mesh_pt->face_boundary(e,2);
789  // if (boundary_id!=0)
790  // {
791  // if (!bound_node_done(loc_node_pt,boundary_id-1))
792  // {
793  // add_boundary_node(boundary_id-1,loc_node_pt);
794  // bound_node_done(loc_node_pt,boundary_id-1)=true;
795  // }
796  // }
797  // }
798 
799 
800  // // These are the face nodes on the element's face 3:
801  // if ( (j==7) || (j==8) || (j==9) )
802  // {
803  // boundary_id=Tmp_mesh_pt->face_boundary(e,3);
804  // if (boundary_id!=0)
805  // {
806  // if (!bound_node_done(loc_node_pt,boundary_id-1))
807  // {
808  // add_boundary_node(boundary_id-1,loc_node_pt);
809  // bound_node_done(loc_node_pt,boundary_id-1)=true;
810  // }
811  // }
812  // }
813 
814  // } //end for j
815 
816  // */
817 
818  } // end for e
819 
820 
821  } // end function
822 
823  //=====================================================================
824  /// Helper function to set up the reverse look up scheme for facets.
825  /// This is used to set up boundary coordinates.
826  //=====================================================================
827  template<class ELEMENT>
829  TetMeshFacetedSurface* const& faceted_surface_pt)
830  {
831  // Set up reverse lookup scheme for a given faceted surface
832  // Assumption is that there there is one boundary ID per facet.
833  unsigned n_facet = faceted_surface_pt->nfacet();
834  for (unsigned f = 0; f < n_facet; f++)
835  {
836  unsigned b = faceted_surface_pt->one_based_facet_boundary_id(f);
837  if (b != 0)
838  {
839  this->Tet_mesh_faceted_surface_pt[b - 1] = faceted_surface_pt;
840  this->Tet_mesh_facet_pt[b - 1] = faceted_surface_pt->facet_pt(f);
841  }
842  else
843  {
844  std::ostringstream error_message;
845  error_message << "Boundary IDs have to be one-based. Yours is " << b
846  << "\n";
847  throw OomphLibError(error_message.str(),
848  OOMPH_CURRENT_FUNCTION,
849  OOMPH_EXCEPTION_LOCATION);
850  }
851  }
852  }
853 
854 
855  /// ////////////////////////////////////////////////////////////////////
856  /// ////////////////////////////////////////////////////////////////////
857  /// ////////////////////////////////////////////////////////////////////
858 
859 
860  //=========================================================================
861  /// Transfer tetgenio data from the input to the output
862  /// The output is assumed to have been constructed and "empty"
863  //=========================================================================
864  template<class ELEMENT>
865  void TetgenMesh<ELEMENT>::deep_copy_of_tetgenio(tetgenio* const& input_pt,
866  tetgenio*& output_pt)
867  {
868  // Copy data values
869  output_pt->firstnumber = input_pt->firstnumber;
870  output_pt->mesh_dim = input_pt->mesh_dim;
871  output_pt->useindex = input_pt->useindex;
872 
873  // Copy the number of points
874  output_pt->numberofpoints = input_pt->numberofpoints;
875  output_pt->numberofpointattributes = input_pt->numberofpointattributes;
876  output_pt->numberofpointmtrs = input_pt->numberofpointmtrs;
877 
878  const int n_point = output_pt->numberofpoints;
879  if (n_point > 0)
880  {
881  output_pt->pointlist = new double[n_point * 3];
882  // Copy the values
883  for (int n = 0; n < (n_point * 3); ++n)
884  {
885  output_pt->pointlist[n] = input_pt->pointlist[n];
886  }
887 
888  // If there are point markers copy as well
889  if (input_pt->pointmarkerlist != (int*)NULL)
890  {
891  output_pt->pointmarkerlist = new int[n_point];
892  for (int n = 0; n < n_point; ++n)
893  {
894  output_pt->pointmarkerlist[n] = input_pt->pointmarkerlist[n];
895  }
896  }
897 
898  const int n_attr = output_pt->numberofpointattributes;
899  if (n_attr > 0)
900  {
901  output_pt->pointattributelist = new double[n_point * n_attr];
902  for (int n = 0; n < (n_point * n_attr); ++n)
903  {
904  output_pt->pointattributelist[n] = input_pt->pointattributelist[n];
905  }
906  }
907  } // End of point case
908 
909  // Now add in metric tensors (if there are any)
910  const int n_point_mtr = output_pt->numberofpointmtrs;
911  if (n_point_mtr > 0)
912  {
913  output_pt->pointmtrlist = new double[n_point * n_point_mtr];
914  for (int n = 0; n < (n_point * n_point_mtr); ++n)
915  {
916  output_pt->pointmtrlist[n] = input_pt->pointmtrlist[n];
917  }
918  }
919 
920  // Next tetrahedrons
921  output_pt->numberoftetrahedra = input_pt->numberoftetrahedra;
922  output_pt->numberofcorners = input_pt->numberofcorners;
923  output_pt->numberoftetrahedronattributes =
924  input_pt->numberoftetrahedronattributes;
925 
926  const int n_tetra = output_pt->numberoftetrahedra;
927  const int n_corner = output_pt->numberofcorners;
928  if (n_tetra > 0)
929  {
930  output_pt->tetrahedronlist = new int[n_tetra * n_corner];
931  for (int n = 0; n < (n_tetra * n_corner); ++n)
932  {
933  output_pt->tetrahedronlist[n] = input_pt->tetrahedronlist[n];
934  }
935 
936  // Add in the volume constriants
937  if (input_pt->tetrahedronvolumelist != (double*)NULL)
938  {
939  output_pt->tetrahedronvolumelist = new double[n_tetra];
940  for (int n = 0; n < n_tetra; ++n)
941  {
942  output_pt->tetrahedronvolumelist[n] =
943  input_pt->tetrahedronvolumelist[n];
944  }
945  }
946 
947  // Add in the attributes
948  const int n_tetra_attr = output_pt->numberoftetrahedronattributes;
949  if (n_tetra_attr > 0)
950  {
951  output_pt->tetrahedronattributelist =
952  new double[n_tetra * n_tetra_attr];
953  for (int n = 0; n < (n_tetra * n_tetra_attr); ++n)
954  {
955  output_pt->tetrahedronattributelist[n] =
956  input_pt->tetrahedronattributelist[n];
957  }
958  }
959 
960  // Add in the neighbour list
961  if (input_pt->neighborlist != (int*)NULL)
962  {
963  output_pt->neighborlist = new int[n_tetra * 4];
964  for (int n = 0; n < (n_tetra * 4); ++n)
965  {
966  output_pt->neighborlist = input_pt->neighborlist;
967  }
968  }
969  } // End of tetrahedron section
970 
971  // Now arrange the facet list
972  output_pt->numberoffacets = input_pt->numberoffacets;
973  const int n_facet = output_pt->numberoffacets;
974  if (n_facet > 0)
975  {
976  output_pt->facetlist = new tetgenio::facet[n_facet];
977  for (int n = 0; n < n_facet; ++n)
978  {
979  tetgenio::facet* input_f_pt = &input_pt->facetlist[n];
980  tetgenio::facet* output_f_pt = &output_pt->facetlist[n];
981 
982  // Copy polygons and holes from the facets
983  output_f_pt->numberofpolygons = input_f_pt->numberofpolygons;
984 
985  // Loop over polygons
986  const int n_poly = output_f_pt->numberofpolygons;
987  if (n_poly > 0)
988  {
989  output_f_pt->polygonlist = new tetgenio::polygon[n_poly];
990  for (int p = 0; p < n_poly; ++p)
991  {
992  tetgenio::polygon* output_p_pt = &output_f_pt->polygonlist[p];
993  tetgenio::polygon* input_p_pt = &input_f_pt->polygonlist[p];
994  // Find out how many vertices each polygon has
995  output_p_pt->numberofvertices = input_p_pt->numberofvertices;
996  // Now copy of the vertices
997  const int n_vertex = output_p_pt->numberofvertices;
998  if (n_vertex > 0)
999  {
1000  output_p_pt->vertexlist = new int[n_vertex];
1001  for (int v = 0; v < n_vertex; ++v)
1002  {
1003  output_p_pt->vertexlist[v] = input_p_pt->vertexlist[v];
1004  }
1005  }
1006  }
1007  }
1008 
1009  // Hole information
1010  output_f_pt->numberofholes = input_f_pt->numberofholes;
1011  const int n_hole = output_f_pt->numberofholes;
1012  if (n_hole > 0)
1013  {
1014  throw OomphLibError("Don't know how to handle holes yet\n",
1015  OOMPH_CURRENT_FUNCTION,
1016  OOMPH_EXCEPTION_LOCATION);
1017  }
1018  } // end of loop over facets
1019 
1020  // Add the facetmarkers if there are any
1021  if (input_pt->facetmarkerlist != (int*)NULL)
1022  {
1023  output_pt->facetmarkerlist = new int[n_facet];
1024  for (int n = 0; n < n_facet; ++n)
1025  {
1026  output_pt->facetmarkerlist[n] = input_pt->facetmarkerlist[n];
1027  }
1028  }
1029  }
1030 
1031  // Now the holes
1032  output_pt->numberofholes = input_pt->numberofholes;
1033  const int n_hole = output_pt->numberofholes;
1034  if (n_hole > 0)
1035  {
1036  output_pt->holelist = new double[n_hole * 3];
1037  for (int n = 0; n < (n_hole * 3); ++n)
1038  {
1039  output_pt->holelist[n] = input_pt->holelist[n];
1040  }
1041  }
1042 
1043  // Now the regions
1044  output_pt->numberofregions = input_pt->numberofregions;
1045  const int n_region = output_pt->numberofregions;
1046  if (n_region > 0)
1047  {
1048  output_pt->regionlist = new double[n_region * 5];
1049  for (int n = 0; n < (n_region * 5); ++n)
1050  {
1051  output_pt->regionlist[n] = input_pt->regionlist[n];
1052  }
1053  }
1054 
1055  // Now the facet constraints
1056  output_pt->numberoffacetconstraints = input_pt->numberoffacetconstraints;
1057  const int n_facet_const = output_pt->numberoffacetconstraints;
1058  if (n_facet_const > 0)
1059  {
1060  output_pt->facetconstraintlist = new double[n_facet_const * 2];
1061  for (int n = 0; n < (n_facet_const * 2); ++n)
1062  {
1063  output_pt->facetconstraintlist[n] = input_pt->facetconstraintlist[n];
1064  }
1065  }
1066 
1067  // Now the segment constraints
1068  output_pt->numberofsegmentconstraints =
1069  input_pt->numberofsegmentconstraints;
1070  const int n_seg_const = output_pt->numberofsegmentconstraints;
1071  if (n_seg_const > 0)
1072  {
1073  output_pt->segmentconstraintlist = new double[n_seg_const * 2];
1074  for (int n = 0; n < (n_seg_const * 2); ++n)
1075  {
1076  output_pt->segmentconstraintlist[n] =
1077  input_pt->segmentconstraintlist[n];
1078  }
1079  }
1080 
1081  // Now the face list
1082  output_pt->numberoftrifaces = input_pt->numberoftrifaces;
1083  const int n_tri_face = output_pt->numberoftrifaces;
1084  if (n_tri_face > 0)
1085  {
1086  output_pt->trifacelist = new int[n_tri_face * 3];
1087  for (int n = 0; n < (n_tri_face * 3); ++n)
1088  {
1089  output_pt->trifacelist[n] = input_pt->trifacelist[n];
1090  }
1091 
1092  output_pt->trifacemarkerlist = new int[n_tri_face];
1093  for (int n = 0; n < n_tri_face; ++n)
1094  {
1095  output_pt->trifacemarkerlist[n] = input_pt->trifacemarkerlist[n];
1096  }
1097  }
1098 
1099  // Now the edge list
1100  output_pt->numberofedges = input_pt->numberofedges;
1101  const int n_edge = output_pt->numberofedges;
1102  if (n_edge > 0)
1103  {
1104  output_pt->edgelist = new int[n_edge * 2];
1105  for (int n = 0; n < (n_edge * 2); ++n)
1106  {
1107  output_pt->edgelist[n] = input_pt->edgelist[n];
1108  }
1109 
1110  output_pt->edgemarkerlist = new int[n_edge];
1111  for (int n = 0; n < n_edge; ++n)
1112  {
1113  output_pt->edgemarkerlist[n] = input_pt->edgemarkerlist[n];
1114  }
1115  }
1116  }
1117 
1118 
1119 } // namespace oomph
1120 #endif
void build_from_scaffold(TimeStepper *time_stepper_pt, const bool &use_attributes)
Build mesh from scaffold.
void setup_reverse_lookup_schemes_for_faceted_surface(TetMeshFacetedSurface *const &faceted_surface_pt)
Function to setup the reverse look-up schemes.
void deep_copy_of_tetgenio(tetgenio *const &input_pt, tetgenio *&output_pt)
Transfer tetgenio data from the input to the output The output is assumed to have been constructed an...
////////////////////////////////////////////////////////////////////// //////////////////////////////...