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-2022 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
33#include "../generic/Telements.h"
34#include "../generic/map_matrix.h"
35
36
37namespace oomph
38{
39 /// ////////////////////////////////////////////////////////////////////
40 /// ////////////////////////////////////////////////////////////////////
41 /// ////////////////////////////////////////////////////////////////////
42
43
44 //========================================================================
45 /// Build unstructured tet mesh based on output from scaffold
46 //========================================================================
47 template<class ELEMENT>
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
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
A general Finite Element class.
Definition: elements.h:1313
virtual void local_coordinate_of_node(const unsigned &j, Vector< double > &s) const
Get local coordinates of node j in the element; vector sets its own size (broken virtual)
Definition: elements.h:1842
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object.
Definition: elements.h:2509
virtual double interpolated_x(const Vector< double > &s, const unsigned &i) const
Return FE interpolated coordinate x[i] at local coordinate s.
Definition: elements.cc:3962
virtual Node * construct_boundary_node(const unsigned &n)
Construct the local node n as a boundary node; that is a node that MAY be placed on a mesh boundary a...
Definition: elements.h:2538
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
virtual void get_boundaries_pt(std::set< unsigned > *&boundaries_pt)
Return a pointer to set of mesh boundaries that this node occupies; this will be overloaded by Bounda...
Definition: nodes.h:1365
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
///////////////////////////////////////////////////////////////////// ///////////////////////////////...
Definition: tet_mesh.h:306
TetMeshFacet * facet_pt(const unsigned &j) const
Pointer to j-th facet.
Definition: tet_mesh.h:373
unsigned nfacet() const
Number of facets.
Definition: tet_mesh.h:325
unsigned one_based_facet_boundary_id(const unsigned &j) const
One-based boundary id of j-th facet.
Definition: tet_mesh.h:331
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...
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
A slight extension to the standard template vector class so that we can include "graceful" array rang...
Definition: Vector.h:58
//////////////////////////////////////////////////////////////////// ////////////////////////////////...