simple_cubic_tet_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_SIMPLE_CUBIC_TET_MESH_TEMPLATE_CC
27#define OOMPH_SIMPLE_CUBIC_TET_MESH_TEMPLATE_CC
28
29#include <algorithm>
30
31// Simple 3D tetrahedral mesh class
33#include "../generic/map_matrix.h"
34#include <algorithm>
35
36
37namespace oomph
38{
39 //====================================================================
40 /// Simple tetrahedral mesh - with 24 tet elements constructed within a
41 /// "brick" form for each element block.
42 //====================================================================
43 template<class ELEMENT>
45 TimeStepper* time_stepper_pt)
46 {
47 // Mesh can only be built with 3D Telements.
48 MeshChecker::assert_geometric_element<TElementGeometricBase, ELEMENT>(3);
49
50 // Create space for elements
51 unsigned nelem = Tmp_mesh_pt->nelement();
52 Element_pt.resize(nelem);
53
54 // Create space for nodes
55 unsigned nnode_scaffold = Tmp_mesh_pt->nnode();
56 Node_pt.resize(nnode_scaffold);
57
58 // Set number of boundaries
59 unsigned nbound = Tmp_mesh_pt->nboundary();
60 set_nboundary(nbound);
61
62 // Loop over elements in scaffold mesh, visit their nodes
63 for (unsigned e = 0; e < nelem; e++)
64 {
65 Element_pt[e] = new ELEMENT;
66 }
67
68 // In the first instance build all nodes from within all the elements
69 unsigned nnod_el = Tmp_mesh_pt->finite_element_pt(0)->nnode();
70 // Loop over elements in scaffold mesh, visit their nodes
71 for (unsigned e = 0; e < nelem; e++)
72 {
73 // Loop over all nodes in element
74 for (unsigned j = 0; j < nnod_el; j++)
75 {
76 // Create new node, using the NEW element's construct_node
77 // member function
78 finite_element_pt(e)->construct_node(j, time_stepper_pt);
79 }
80 }
81
82
83 // Setup map to check the (pseudo-)global node number
84 // Nodes whose number is zero haven't been copied across
85 // into the mesh yet.
86 std::map<Node*, unsigned> global_number;
87 unsigned global_count = 0;
88 // Loop over elements in scaffold mesh, visit their nodes
89 for (unsigned e = 0; e < nelem; e++)
90 {
91 // Loop over all nodes in element
92 for (unsigned j = 0; j < nnod_el; j++)
93 {
94 // Pointer to node in the scaffold mesh
95 Node* scaffold_node_pt = Tmp_mesh_pt->finite_element_pt(e)->node_pt(j);
96
97 // Get the (pseudo-)global node number in scaffold mesh
98 // (It's zero [=default] if not visited this one yet)
99 unsigned j_global = global_number[scaffold_node_pt];
100
101 // Haven't done this one yet
102 if (j_global == 0)
103 {
104 // Give it a number (not necessarily the global node
105 // number in the scaffold mesh -- we just need something
106 // to keep track...)
107 global_count++;
108 global_number[scaffold_node_pt] = global_count;
109
110 // Copy new node, created using the NEW element's construct_node
111 // function into global storage, using the same global
112 // node number that we've just associated with the
113 // corresponding node in the scaffold mesh
114 Node_pt[global_count - 1] = finite_element_pt(e)->node_pt(j);
115
116 // Assign coordinates
117 Node_pt[global_count - 1]->x(0) = scaffold_node_pt->x(0);
118 Node_pt[global_count - 1]->x(1) = scaffold_node_pt->x(1);
119 Node_pt[global_count - 1]->x(2) = scaffold_node_pt->x(2);
120
121 // Get pointer to set of mesh boundaries that this
122 // scaffold node occupies; NULL if the node is not on any boundary
123 std::set<unsigned>* boundaries_pt;
124 scaffold_node_pt->get_boundaries_pt(boundaries_pt);
125
126 // Loop over the mesh boundaries that the node in the scaffold mesh
127 // occupies and assign new node to the same ones.
128 if (boundaries_pt != 0)
129 {
130 this->convert_to_boundary_node(Node_pt[global_count - 1]);
131 for (std::set<unsigned>::iterator it = boundaries_pt->begin();
132 it != boundaries_pt->end();
133 ++it)
134 {
135 add_boundary_node(*it, Node_pt[global_count - 1]);
136 }
137 }
138 }
139 // This one has already been done: Kill it
140 else
141 {
142 // Kill it
143 delete finite_element_pt(e)->node_pt(j);
144
145 // Overwrite the element's pointer to local node by
146 // pointer to the already existing node -- identified by
147 // the number kept in the map
148 finite_element_pt(e)->node_pt(j) = Node_pt[j_global - 1];
149 }
150 }
151 }
152
153
154 // At this point we've created all the elements and
155 // created their vertex nodes. Now we need to create
156 // the additional (midside and internal) nodes!
157
158
159 // We'll first create all local nodes for all elements
160 // and then delete the superfluous ones that have
161 // a matching node in an adjacent element.
162
163 // Get number of nodes along element edge and dimension of element (3)
164 // from first element
165 unsigned nnode_1d = finite_element_pt(0)->nnode_1d();
166
167 // At the moment we're only able to deal with nnode_1d=2 or 3.
168 if ((nnode_1d != 2) && (nnode_1d != 3))
169 {
170 std::ostringstream error_message;
171 error_message << "Mesh generation currently only works\n";
172 error_message << "for nnode_1d = 2 or 3. You're trying to use it\n";
173 error_message << "for nnode_1d=" << nnode_1d << std::endl;
174
175 throw OomphLibError(
176 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
177 }
178
179 // Spatial dimension of element = number of local coordinates
180 unsigned dim = finite_element_pt(0)->dim();
181
182 // Storage for the local coordinate of the new node
183 Vector<double> s(dim);
184
185 // Get number of nodes in the element from first element
186 unsigned nnode = finite_element_pt(0)->nnode();
187
188 // Loop over all elements
189 for (unsigned e = 0; e < nelem; e++)
190 {
191 // Loop over the new nodes in the element and create them.
192 for (unsigned j = 4; j < nnode; j++)
193 {
194 // Create new node
195 Node* new_node_pt =
196 finite_element_pt(e)->construct_node(j, time_stepper_pt);
197
198 // What are the node's local coordinates?
199 finite_element_pt(e)->local_coordinate_of_node(j, s);
200
201 // Find the coordinates of the new node from the existing
202 // and fully-functional element in the scaffold mesh
203 for (unsigned i = 0; i < dim; i++)
204 {
205 new_node_pt->x(i) =
206 Tmp_mesh_pt->finite_element_pt(e)->interpolated_x(s, i);
207 }
208 } // end of loop over new nodes
209 } // end of loop over elements
210
211
212 // Bracket this away so the edge map goes out of scope
213 // when we're done
214 {
215 // Storage for pointer to mid-edge node
216 MapMatrix<Node*, Node*> central_edge_node_pt;
217 Node* edge_node1_pt = 0;
218 Node* edge_node2_pt = 0;
219
220 // Loop over elements
221 for (unsigned e = 0; e < nelem; e++)
222 {
223 // Loop over new local nodes
224 for (unsigned j = 4; j < nnode; j++)
225 {
226 // Pointer to the element's local node
227 Node* node_pt = finite_element_pt(e)->node_pt(j);
228
229 // By default, we assume the node is not new
230 bool is_new = false;
231
232 // This will have to be changed for higher-order elements
233 //=======================================================
234
235 // Switch on local node number (all located on edges)
236 switch (j)
237 {
238 // Node 4 is located between nodes 0 and 1
239 case 4:
240
241 edge_node1_pt = finite_element_pt(e)->node_pt(0);
242 edge_node2_pt = finite_element_pt(e)->node_pt(1);
243 if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
244 {
245 is_new = true;
246 central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
247 central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
248 }
249 break;
250
251
252 // Node 5 is located between nodes 0 and 2
253 case 5:
254
255 edge_node1_pt = finite_element_pt(e)->node_pt(0);
256 edge_node2_pt = finite_element_pt(e)->node_pt(2);
257 if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
258 {
259 is_new = true;
260 central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
261 central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
262 }
263 break;
264
265
266 // Node 6 is located between nodes 0 and 3
267 case 6:
268
269 edge_node1_pt = finite_element_pt(e)->node_pt(0);
270 edge_node2_pt = finite_element_pt(e)->node_pt(3);
271 if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
272 {
273 is_new = true;
274 central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
275 central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
276 }
277 break;
278
279
280 // Node 7 is located between nodes 1 and 2
281 case 7:
282
283 edge_node1_pt = finite_element_pt(e)->node_pt(1);
284 edge_node2_pt = finite_element_pt(e)->node_pt(2);
285 if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
286 {
287 is_new = true;
288 central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
289 central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
290 }
291 break;
292
293
294 // Node 8 is located between nodes 2 and 3
295 case 8:
296
297 edge_node1_pt = finite_element_pt(e)->node_pt(2);
298 edge_node2_pt = finite_element_pt(e)->node_pt(3);
299 if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
300 {
301 is_new = true;
302 central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
303 central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
304 }
305 break;
306
307
308 // Node 9 is located between nodes 1 and 3
309 case 9:
310
311 edge_node1_pt = finite_element_pt(e)->node_pt(1);
312 edge_node2_pt = finite_element_pt(e)->node_pt(3);
313 if (central_edge_node_pt(edge_node1_pt, edge_node2_pt) == 0)
314 {
315 is_new = true;
316 central_edge_node_pt(edge_node1_pt, edge_node2_pt) = node_pt;
317 central_edge_node_pt(edge_node2_pt, edge_node1_pt) = node_pt;
318 }
319 break;
320
321 default:
322 throw OomphLibError("More than ten nodes in Tet Element",
323 OOMPH_CURRENT_FUNCTION,
324 OOMPH_EXCEPTION_LOCATION);
325 }
326
327 if (is_new)
328 {
329 // New node: Add it to mesh
330 Node_pt.push_back(node_pt);
331 }
332 else
333 {
334 // Delete local node in element...
335 delete finite_element_pt(e)->node_pt(j);
336
337 // ... and reset pointer to the existing node
338 finite_element_pt(e)->node_pt(j) =
339 central_edge_node_pt(edge_node1_pt, edge_node2_pt);
340 }
341 }
342 }
343 }
344
345
346 // Boundary conditions
347
348 // Matrix map to check if a node has already been added to
349 // the boundary number b (NOTE: Enumerated by pointer to ORIGINAL
350 // node before transfer to boundary node)
352
353 // Create lookup scheme for pointers to original local nodes
354 // in elements
355 Vector<Vector<Node*>> orig_node_pt(nelem);
356
357 // Loop over elements
358 for (unsigned e = 0; e < nelem; e++)
359 {
360 orig_node_pt[e].resize(nnode, 0);
361
362 // Loop over new local nodes
363 for (unsigned j = 4; j < nnode; j++)
364 {
365 // Pointer to the element's local node
366 orig_node_pt[e][j] = finite_element_pt(e)->node_pt(j);
367 }
368 }
369
370
371 // Loop over elements
372 for (unsigned e = 0; e < nelem; e++)
373 {
374 // Loop over new local nodes
375 for (unsigned j = 4; j < nnode; j++)
376 {
377 // Loop over the boundaries
378 for (unsigned bo = 0; bo < nbound; bo++)
379 {
380 // Pointer to the element's local node
381 Node* loc_node_pt = finite_element_pt(e)->node_pt(j);
382
383 // Pointer to original node
384 Node* orig_loc_node_pt = orig_node_pt[e][j];
385
386 // value of the map for the node and boundary specified
387 bool bound_test = bound_node_done(orig_loc_node_pt, bo);
388
389 if (bound_test == false)
390 {
391 bound_node_done(orig_loc_node_pt, bo) = true;
392
393 // This will have to be changed for higher-order elements
394 //=======================================================
395
396 // Switch on local node number (all located on edges)
397 switch (j)
398 {
399 // Node 4 is located between nodes 0 and 1
400 case 4:
401
402 if (finite_element_pt(e)->node_pt(0)->is_on_boundary(bo) &&
403 finite_element_pt(e)->node_pt(1)->is_on_boundary(bo))
404 {
405 this->convert_to_boundary_node(loc_node_pt);
406 add_boundary_node(bo, loc_node_pt);
407 }
408 break;
409
410
411 // Node 5 is located between nodes 0 and 2
412 case 5:
413
414 if (finite_element_pt(e)->node_pt(0)->is_on_boundary(bo) &&
415 finite_element_pt(e)->node_pt(2)->is_on_boundary(bo))
416 {
417 this->convert_to_boundary_node(loc_node_pt);
418 add_boundary_node(bo, loc_node_pt);
419 }
420 break;
421
422
423 // Node 6 is located between nodes 0 and 3
424 case 6:
425
426 if (finite_element_pt(e)->node_pt(0)->is_on_boundary(bo) &&
427 finite_element_pt(e)->node_pt(3)->is_on_boundary(bo))
428 {
429 this->convert_to_boundary_node(loc_node_pt);
430 add_boundary_node(bo, loc_node_pt);
431 }
432 break;
433
434
435 // Node 7 is located between nodes 1 and 2
436 case 7:
437
438 if (finite_element_pt(e)->node_pt(1)->is_on_boundary(bo) &&
439 finite_element_pt(e)->node_pt(2)->is_on_boundary(bo))
440 {
441 this->convert_to_boundary_node(loc_node_pt);
442 add_boundary_node(bo, loc_node_pt);
443 }
444 break;
445
446
447 // Node 8 is located between nodes 2 and 3
448 case 8:
449
450 if (finite_element_pt(e)->node_pt(2)->is_on_boundary(bo) &&
451 finite_element_pt(e)->node_pt(3)->is_on_boundary(bo))
452 {
453 this->convert_to_boundary_node(loc_node_pt);
454 add_boundary_node(bo, loc_node_pt);
455 }
456 break;
457
458
459 // Node 9 is located between nodes 1 and 3
460 case 9:
461
462 if (finite_element_pt(e)->node_pt(1)->is_on_boundary(bo) &&
463 finite_element_pt(e)->node_pt(3)->is_on_boundary(bo))
464 {
465 this->convert_to_boundary_node(loc_node_pt);
466 add_boundary_node(bo, loc_node_pt);
467 }
468 break;
469
470
471 default:
472 throw OomphLibError("More than ten nodes in Tet Element",
473 OOMPH_CURRENT_FUNCTION,
474 OOMPH_EXCEPTION_LOCATION);
475 }
476 }
477
478 } // end for bo
479 } // end for j
480 } // end for e
481 }
482
483} // namespace oomph
484
485#endif
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
MapMatrixMixed is a generalised, STL-map-based, sparse(ish) matrix class with mixed indices.
Definition: map_matrix.h:109
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: map_matrix.h:508
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....
void build_from_scaffold(TimeStepper *time_stepper_pt)
Build mesh from scaffold mesh.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
//////////////////////////////////////////////////////////////////// ////////////////////////////////...