tube_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-2024 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_TUBE_MESH_TEMPLATE_CC
27#define OOMPH_TUBE_MESH_TEMPLATE_CC
28
29#include "tube_mesh.template.h"
30
31
32namespace oomph
33{
34 //====================================================================
35 /// Constructor for deformable quarter tube mesh class.
36 /// The domain is specified by the GeomObject that
37 /// identifies the entire volume.
38 //====================================================================
39 template<class ELEMENT>
40 TubeMesh<ELEMENT>::TubeMesh(GeomObject* volume_pt,
41 const Vector<double>& centreline_limits,
42 const Vector<double>& theta_positions,
43 const Vector<double>& radius_box,
44 const unsigned& nlayer,
45 TimeStepper* time_stepper_pt)
46 : Volume_pt(volume_pt)
47 {
48 // Mesh can only be built with 3D Qelements.
49 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
50
51 // Build macro element-based domain
52 Domain_pt = new TubeDomain(
53 volume_pt, centreline_limits, theta_positions, radius_box, nlayer);
54
55 // Set the number of boundaries
56 set_nboundary(3);
57
58 // We have only bothered to parametrise boundary 1
59 Boundary_coordinate_exists[1] = true;
60
61 // Allocate the store for the elements
62 unsigned nelem = 5 * nlayer;
63 Element_pt.resize(nelem);
64
65 // Create dummy element so we can determine the number of nodes
66 ELEMENT* dummy_el_pt = new ELEMENT;
67
68 // Read out the number of linear points in the element
69 unsigned n_p = dummy_el_pt->nnode_1d();
70
71 // Kill the element
72 delete dummy_el_pt;
73
74 // Can now allocate the store for the nodes
75 unsigned nnodes_total =
76 (n_p * n_p + 4 * (n_p - 1) * (n_p - 1)) * (1 + nlayer * (n_p - 1));
77 Node_pt.resize(nnodes_total);
78
79 Vector<double> s(3);
80 Vector<double> r(3);
81
82 // Storage for the intrinsic boundary coordinate
83 Vector<double> zeta(2);
84
85 // Loop over elements and create all nodes
86 for (unsigned ielem = 0; ielem < nelem; ielem++)
87 {
88 // Create element
89 Element_pt[ielem] = new ELEMENT;
90
91 // Loop over rows in z/s_2-direction
92 for (unsigned i2 = 0; i2 < n_p; i2++)
93 {
94 // Loop over rows in y/s_1-direction
95 for (unsigned i1 = 0; i1 < n_p; i1++)
96 {
97 // Loop over rows in x/s_0-direction
98 for (unsigned i0 = 0; i0 < n_p; i0++)
99 {
100 // Local node number
101 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
102
103 // Create the node
104 Node* node_pt = finite_element_pt(ielem)->construct_node(
105 jnod_local, time_stepper_pt);
106
107 // Set the position of the node from macro element mapping
108 s[0] = -1.0 + 2.0 * double(i0) / double(n_p - 1);
109 s[1] = -1.0 + 2.0 * double(i1) / double(n_p - 1);
110 s[2] = -1.0 + 2.0 * double(i2) / double(n_p - 1);
111 Domain_pt->macro_element_pt(ielem)->macro_map(s, r);
112
113 node_pt->x(0) = r[0];
114 node_pt->x(1) = r[1];
115 node_pt->x(2) = r[2];
116 }
117 }
118 }
119 }
120
121 // Initialise number of global nodes
122 unsigned node_count = 0;
123
124 // Tolerance for node killing:
125 double node_kill_tol = 1.0e-12;
126
127 // Check for error in node killing
128 bool stopit = false;
129
130 // Define pine
131 const double pi = MathematicalConstants::Pi;
132
133 // Loop over elements
134 for (unsigned ielem = 0; ielem < nelem; ielem++)
135 {
136 // Which layer?
137 unsigned ilayer = unsigned(ielem / 5);
138
139 // Which macro element?
140 switch (ielem % 5)
141 {
142 // Macro element 0: Central box
143 //-----------------------------
144 case 0:
145
146 // Loop over rows in z/s_2-direction
147 for (unsigned i2 = 0; i2 < n_p; i2++)
148 {
149 // Loop over rows in y/s_1-direction
150 for (unsigned i1 = 0; i1 < n_p; i1++)
151 {
152 // Loop over rows in x/s_0-direction
153 for (unsigned i0 = 0; i0 < n_p; i0++)
154 {
155 // Local node number
156 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
157
158 // Has the node been killed?
159 bool killed = false;
160
161 // First layer of all nodes in s_2 direction gets killed
162 // and re-directed to nodes in previous element layer
163 if ((i2 == 0) && (ilayer > 0))
164 {
165 // Neighbour element
166 unsigned ielem_neigh = ielem - 5;
167
168 // Node in neighbour element
169 unsigned i0_neigh = i0;
170 unsigned i1_neigh = i1;
171 unsigned i2_neigh = n_p - 1;
172 unsigned jnod_local_neigh =
173 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
174
175 // Check:
176 for (unsigned i = 0; i < 3; i++)
177 {
178 double error = std::fabs(
179 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
180 finite_element_pt(ielem_neigh)
181 ->node_pt(jnod_local_neigh)
182 ->x(i));
183 if (error > node_kill_tol)
184 {
185 oomph_info << "Error in node killing for i " << i << " "
186 << error << std::endl;
187 stopit = true;
188 }
189 }
190
191 // Kill node
192 delete finite_element_pt(ielem)->node_pt(jnod_local);
193 killed = true;
194
195 // Set pointer to neighbour:
196 finite_element_pt(ielem)->node_pt(jnod_local) =
197 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
198 }
199
200 // No duplicate node: Copy across to mesh
201 if (!killed)
202 {
203 Node_pt[node_count] =
204 finite_element_pt(ielem)->node_pt(jnod_local);
205
206 // Set boundaries:
207
208 // Back: Boundary 0
209 if ((i2 == 0) && (ilayer == 0))
210 {
211 this->convert_to_boundary_node(Node_pt[node_count]);
212 add_boundary_node(0, Node_pt[node_count]);
213 }
214
215 // Front: Boundary 2
216 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
217 {
218 this->convert_to_boundary_node(Node_pt[node_count]);
219 add_boundary_node(2, Node_pt[node_count]);
220 }
221
222 // Increment node counter
223 node_count++;
224 }
225 }
226 }
227 }
228
229
230 break;
231
232 // Macro element 1: Bottom
233 //---------------------------------
234 case 1:
235
236
237 // Loop over rows in z/s_2-direction
238 for (unsigned i2 = 0; i2 < n_p; i2++)
239 {
240 // Loop over rows in y/s_1-direction
241 for (unsigned i1 = 0; i1 < n_p; i1++)
242 {
243 // Loop over rows in x/s_0-direction
244 for (unsigned i0 = 0; i0 < n_p; i0++)
245 {
246 // Local node number
247 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
248
249 // Has the node been killed?
250 bool killed = false;
251
252 // First layer of all nodes in s_2 direction gets killed
253 // and re-directed to nodes in previous element layer
254 if ((i2 == 0) && (ilayer > 0))
255 {
256 // Neighbour element
257 unsigned ielem_neigh = ielem - 5;
258
259 // Node in neighbour element
260 unsigned i0_neigh = i0;
261 unsigned i1_neigh = i1;
262 unsigned i2_neigh = n_p - 1;
263 unsigned jnod_local_neigh =
264 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
265
266
267 // Check:
268 for (unsigned i = 0; i < 3; i++)
269 {
270 double error = std::fabs(
271 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
272 finite_element_pt(ielem_neigh)
273 ->node_pt(jnod_local_neigh)
274 ->x(i));
275 if (error > node_kill_tol)
276 {
277 oomph_info << "Error in node killing for i " << i << " "
278 << error << std::endl;
279 stopit = true;
280 }
281 }
282
283 // Kill node
284 delete finite_element_pt(ielem)->node_pt(jnod_local);
285 killed = true;
286
287 // Set pointer to neighbour:
288 finite_element_pt(ielem)->node_pt(jnod_local) =
289 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
290 }
291 // Not in first layer:
292 else
293 {
294 // Duplicate node: kill and set pointer to central element
295 if (i1 == (n_p - 1))
296 {
297 // Neighbour element
298 unsigned ielem_neigh = ielem - 1;
299
300 // Node in neighbour element
301 unsigned i0_neigh = i0;
302 unsigned i1_neigh = 0;
303 unsigned i2_neigh = i2;
304 unsigned jnod_local_neigh =
305 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
306
307 // Check:
308 for (unsigned i = 0; i < 3; i++)
309 {
310 double error = std::fabs(
311 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
312 finite_element_pt(ielem_neigh)
313 ->node_pt(jnod_local_neigh)
314 ->x(i));
315 if (error > node_kill_tol)
316 {
317 oomph_info << "Error in node killing for i " << i << " "
318 << error << std::endl;
319 stopit = true;
320 }
321 }
322
323 // Kill node
324 delete finite_element_pt(ielem)->node_pt(jnod_local);
325 killed = true;
326
327 // Set pointer to neighbour:
328 finite_element_pt(ielem)->node_pt(jnod_local) =
329 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
330 }
331 }
332
333 // No duplicate node: Copy across to mesh
334 if (!killed)
335 {
336 Node_pt[node_count] =
337 finite_element_pt(ielem)->node_pt(jnod_local);
338
339 // Set boundaries:
340
341 // Back: Boundary 0
342 if ((i2 == 0) && (ilayer == 0))
343 {
344 this->convert_to_boundary_node(Node_pt[node_count]);
345 add_boundary_node(0, Node_pt[node_count]);
346 }
347
348 // Front: Boundary 2
349 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
350 {
351 this->convert_to_boundary_node(Node_pt[node_count]);
352 add_boundary_node(2, Node_pt[node_count]);
353 }
354
355 // Bottom: outer boundary
356 if (i1 == 0)
357 {
358 this->convert_to_boundary_node(Node_pt[node_count]);
359 add_boundary_node(1, Node_pt[node_count]);
360
361 // Get axial boundary coordinate
362 zeta[0] = centreline_limits[0] +
363 (double(ilayer) + double(i2) / double(n_p - 1)) *
364 (centreline_limits[1] - centreline_limits[0]) /
365 double(nlayer);
366
367 // Get azimuthal boundary coordinate
368 zeta[1] = theta_positions[0] +
369 double(i1) / double(n_p - 1) * 0.5 *
370 (theta_positions[1] - theta_positions[0]);
371
372 Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
373 }
374 // Increment node counter
375 node_count++;
376 }
377 }
378 }
379 } // End of loop over nodes
380
381 break;
382
383
384 // Macro element 2: Right element
385 //--------------------------------
386 case 2:
387
388 // Loop over rows in z/s_2-direction
389 for (unsigned i2 = 0; i2 < n_p; i2++)
390 {
391 // Loop over rows in y/s_1-direction
392 for (unsigned i1 = 0; i1 < n_p; i1++)
393 {
394 // Loop over rows in x/s_0-direction
395 for (unsigned i0 = 0; i0 < n_p; i0++)
396 {
397 // Local node number
398 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
399
400 // Has the node been killed?
401 bool killed = false;
402
403 // First layer of all nodes in s_2 direction gets killed
404 // and re-directed to nodes in previous element layer
405 if ((i2 == 0) && (ilayer > 0))
406 {
407 // Neighbour element
408 unsigned ielem_neigh = ielem - 5;
409
410 // Node in neighbour element
411 unsigned i0_neigh = i0;
412 unsigned i1_neigh = i1;
413 unsigned i2_neigh = n_p - 1;
414 unsigned jnod_local_neigh =
415 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
416
417 // Check:
418 for (unsigned i = 0; i < 3; i++)
419 {
420 double error = std::fabs(
421 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
422 finite_element_pt(ielem_neigh)
423 ->node_pt(jnod_local_neigh)
424 ->x(i));
425 if (error > node_kill_tol)
426 {
427 oomph_info << "Error in node killing for i " << i << " "
428 << error << std::endl;
429 stopit = true;
430 }
431 }
432
433 // Kill node
434 delete finite_element_pt(ielem)->node_pt(jnod_local);
435 killed = true;
436
437 // Set pointer to neighbour:
438 finite_element_pt(ielem)->node_pt(jnod_local) =
439 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
440 }
441 // Not in first layer:
442 else
443 {
444 // Duplicate node: kill and set pointer to previous element
445 if (i1 == 0)
446 {
447 // Neighbour element
448 unsigned ielem_neigh = ielem - 1;
449
450 // Node in neighbour element
451 unsigned i0_neigh = n_p - 1;
452 unsigned i1_neigh = n_p - 1 - i0;
453 unsigned i2_neigh = i2;
454 unsigned jnod_local_neigh =
455 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
456
457 // Check:
458 for (unsigned i = 0; i < 3; i++)
459 {
460 double error = std::fabs(
461 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
462 finite_element_pt(ielem_neigh)
463 ->node_pt(jnod_local_neigh)
464 ->x(i));
465 if (error > node_kill_tol)
466 {
467 oomph_info << "Error in node killing for i " << i << " "
468 << error << std::endl;
469 stopit = true;
470 }
471 }
472
473 // Kill node
474 delete finite_element_pt(ielem)->node_pt(jnod_local);
475 killed = true;
476
477 // Set pointer to neighbour:
478 finite_element_pt(ielem)->node_pt(jnod_local) =
479 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
480 }
481
482
483 // Duplicate node: kill and set pointer to central element
484 if ((i0 == 0) && (i1 != 0))
485 {
486 // Neighbour element
487 unsigned ielem_neigh = ielem - 2;
488
489 // Node in neighbour element
490 unsigned i0_neigh = n_p - 1;
491 unsigned i1_neigh = i1;
492 unsigned i2_neigh = i2;
493 unsigned jnod_local_neigh =
494 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
495
496 // Check:
497 for (unsigned i = 0; i < 3; i++)
498 {
499 double error = std::fabs(
500 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
501 finite_element_pt(ielem_neigh)
502 ->node_pt(jnod_local_neigh)
503 ->x(i));
504 if (error > node_kill_tol)
505 {
506 oomph_info << "Error in node killing for i " << i << " "
507 << error << std::endl;
508 stopit = true;
509 }
510 }
511
512 // Kill node
513 delete finite_element_pt(ielem)->node_pt(jnod_local);
514 killed = true;
515
516 // Set pointer to neighbour:
517 finite_element_pt(ielem)->node_pt(jnod_local) =
518 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
519 }
520
521 // No duplicate node: Copy across to mesh
522 if (!killed)
523 {
524 Node_pt[node_count] =
525 finite_element_pt(ielem)->node_pt(jnod_local);
526
527 // Set boundaries:
528
529 // Back: Boundary 0
530 if ((i2 == 0) && (ilayer == 0))
531 {
532 this->convert_to_boundary_node(Node_pt[node_count]);
533 add_boundary_node(0, Node_pt[node_count]);
534 }
535
536 // Front: Boundary 2
537 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
538 {
539 this->convert_to_boundary_node(Node_pt[node_count]);
540 add_boundary_node(2, Node_pt[node_count]);
541 }
542
543 // Tube boundary
544 if (i0 == n_p - 1)
545 {
546 this->convert_to_boundary_node(Node_pt[node_count]);
547 add_boundary_node(1, Node_pt[node_count]);
548
549 // Get axial boundary coordinate
550 zeta[0] =
551 centreline_limits[0] +
552 (double(ilayer) + double(i2) / double(n_p - 1)) *
553 (centreline_limits[1] - centreline_limits[0]) /
554 double(nlayer);
555
556 // Get azimuthal boundary coordinate
557 zeta[1] = theta_positions[1] +
558 double(i1) / double(n_p - 1) * 0.5 *
559 (theta_positions[2] - theta_positions[1]);
560
561 Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
562 }
563
564 // Increment node counter
565 node_count++;
566 }
567 }
568 }
569 }
570 }
571
572 break;
573
574 // Macro element 3: Top element
575 //--------------------------------
576 case 3:
577
578 // Loop over rows in z/s_2-direction
579 for (unsigned i2 = 0; i2 < n_p; i2++)
580 {
581 // Loop over rows in y/s_1-direction
582 for (unsigned i1 = 0; i1 < n_p; i1++)
583 {
584 // Loop over rows in x/s_0-direction
585 for (unsigned i0 = 0; i0 < n_p; i0++)
586 {
587 // Local node number
588 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
589
590 // Has the node been killed?
591 bool killed = false;
592
593 // First layer of all nodes in s_2 direction gets killed
594 // and re-directed to nodes in previous element layer
595 if ((i2 == 0) && (ilayer > 0))
596 {
597 // Neighbour element
598 unsigned ielem_neigh = ielem - 5;
599
600 // Node in neighbour element
601 unsigned i0_neigh = i0;
602 unsigned i1_neigh = i1;
603 unsigned i2_neigh = n_p - 1;
604 unsigned jnod_local_neigh =
605 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
606
607 // Check:
608 for (unsigned i = 0; i < 3; i++)
609 {
610 double error = std::fabs(
611 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
612 finite_element_pt(ielem_neigh)
613 ->node_pt(jnod_local_neigh)
614 ->x(i));
615 if (error > node_kill_tol)
616 {
617 oomph_info << "Error in node killing for i " << i << " "
618 << error << std::endl;
619 stopit = true;
620 }
621 }
622
623 // Kill node
624 delete finite_element_pt(ielem)->node_pt(jnod_local);
625 killed = true;
626
627 // Set pointer to neighbour:
628 finite_element_pt(ielem)->node_pt(jnod_local) =
629 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
630 }
631 // Not in first layer:
632 else
633 {
634 // Duplicate node: kill and set pointer to previous element
635 if (i0 == n_p - 1)
636 {
637 // Neighbour element
638 unsigned ielem_neigh = ielem - 1;
639
640 // Node in neighbour element
641 unsigned i0_neigh = i1;
642 unsigned i1_neigh = n_p - 1;
643 unsigned i2_neigh = i2;
644 unsigned jnod_local_neigh =
645 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
646
647 // Check:
648 for (unsigned i = 0; i < 3; i++)
649 {
650 double error = std::fabs(
651 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
652 finite_element_pt(ielem_neigh)
653 ->node_pt(jnod_local_neigh)
654 ->x(i));
655 if (error > node_kill_tol)
656 {
657 oomph_info << "Error in node killing for i " << i << " "
658 << error << std::endl;
659 stopit = true;
660 }
661 }
662
663 // Kill node
664 delete finite_element_pt(ielem)->node_pt(jnod_local);
665 killed = true;
666
667 // Set pointer to neighbour:
668 finite_element_pt(ielem)->node_pt(jnod_local) =
669 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
670 }
671
672 // Duplicate node: kill and set pointer to central element
673 if ((i1 == 0) && (i0 != n_p - 1))
674 {
675 // Neighbour element
676 unsigned ielem_neigh = ielem - 3;
677
678 // Node in neighbour element
679 unsigned i0_neigh = i0;
680 unsigned i1_neigh = n_p - 1;
681 unsigned i2_neigh = i2;
682 unsigned jnod_local_neigh =
683 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
684
685 // Check:
686 for (unsigned i = 0; i < 3; i++)
687 {
688 double error = std::fabs(
689 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
690 finite_element_pt(ielem_neigh)
691 ->node_pt(jnod_local_neigh)
692 ->x(i));
693 if (error > node_kill_tol)
694 {
695 oomph_info << "Error in node killing for i " << i << " "
696 << error << std::endl;
697 stopit = true;
698 }
699 }
700
701 // Kill node
702 delete finite_element_pt(ielem)->node_pt(jnod_local);
703 killed = true;
704
705 // Set pointer to neighbour:
706 finite_element_pt(ielem)->node_pt(jnod_local) =
707 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
708 }
709
710 // No duplicate node: Copy across to mesh
711 if (!killed)
712 {
713 Node_pt[node_count] =
714 finite_element_pt(ielem)->node_pt(jnod_local);
715
716 // Set boundaries:
717
718 // Back: Boundary 0
719 if ((i2 == 0) && (ilayer == 0))
720 {
721 this->convert_to_boundary_node(Node_pt[node_count]);
722 add_boundary_node(0, Node_pt[node_count]);
723 }
724
725 // Front: Boundary 2
726 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
727 {
728 this->convert_to_boundary_node(Node_pt[node_count]);
729 add_boundary_node(2, Node_pt[node_count]);
730 }
731
732 // Tube boundary
733 if (i1 == n_p - 1)
734 {
735 this->convert_to_boundary_node(Node_pt[node_count]);
736 add_boundary_node(1, Node_pt[node_count]);
737
738 // Get axial boundary coordinate
739 zeta[0] =
740 centreline_limits[0] +
741 (double(ilayer) + double(i2) / double(n_p - 1)) *
742 (centreline_limits[1] - centreline_limits[0]) /
743 double(nlayer);
744
745 // Get azimuthal boundary coordinate
746 zeta[1] = theta_positions[3] +
747 double(i0) / double(n_p - 1) * 0.5 *
748 (theta_positions[2] - theta_positions[3]);
749
750 Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
751 }
752
753 // Increment node counter
754 node_count++;
755 }
756 }
757 }
758 }
759 }
760 break;
761
762
763 // Macro element 4: Left element
764 //--------------------------------
765 case 4:
766
767 // Loop over rows in z/s_2-direction
768 for (unsigned i2 = 0; i2 < n_p; i2++)
769 {
770 // Loop over rows in y/s_1-direction
771 for (unsigned i1 = 0; i1 < n_p; i1++)
772 {
773 // Loop over rows in x/s_0-direction
774 for (unsigned i0 = 0; i0 < n_p; i0++)
775 {
776 // Local node number
777 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
778
779 // Has the node been killed?
780 bool killed = false;
781
782 // First layer of all nodes in s_2 direction gets killed
783 // and re-directed to nodes in previous element layer
784 if ((i2 == 0) && (ilayer > 0))
785 {
786 // Neighbour element
787 unsigned ielem_neigh = ielem - 5;
788
789 // Node in neighbour element
790 unsigned i0_neigh = i0;
791 unsigned i1_neigh = i1;
792 unsigned i2_neigh = n_p - 1;
793 unsigned jnod_local_neigh =
794 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
795
796 // Check:
797 for (unsigned i = 0; i < 3; i++)
798 {
799 double error = std::fabs(
800 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
801 finite_element_pt(ielem_neigh)
802 ->node_pt(jnod_local_neigh)
803 ->x(i));
804 if (error > node_kill_tol)
805 {
806 oomph_info << "Error in node killing for i " << i << " "
807 << error << std::endl;
808 stopit = true;
809 }
810 }
811
812 // Kill node
813 delete finite_element_pt(ielem)->node_pt(jnod_local);
814 killed = true;
815
816 // Set pointer to neighbour:
817 finite_element_pt(ielem)->node_pt(jnod_local) =
818 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
819 }
820 // Not in first layer:
821 else
822 {
823 // Duplicate node: kill and set pointer to previous element
824 if (i1 == n_p - 1)
825 {
826 // Neighbour element
827 unsigned ielem_neigh = ielem - 1;
828
829 // Node in neighbour element
830 unsigned i0_neigh = 0;
831 unsigned i1_neigh = n_p - 1 - i0;
832 unsigned i2_neigh = i2;
833 unsigned jnod_local_neigh =
834 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
835
836 // Check:
837 for (unsigned i = 0; i < 3; i++)
838 {
839 double error = std::fabs(
840 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
841 finite_element_pt(ielem_neigh)
842 ->node_pt(jnod_local_neigh)
843 ->x(i));
844 if (error > node_kill_tol)
845 {
846 oomph_info << "Error in node killing for i " << i << " "
847 << error << std::endl;
848 stopit = true;
849 }
850 }
851
852 // Kill node
853 delete finite_element_pt(ielem)->node_pt(jnod_local);
854 killed = true;
855
856 // Set pointer to neighbour:
857 finite_element_pt(ielem)->node_pt(jnod_local) =
858 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
859 }
860
861 // Duplicate node: kill and set pointer to central element
862 if ((i0 == n_p - 1) && (i1 != n_p - 1))
863 {
864 // Neighbour element
865 unsigned ielem_neigh = ielem - 4;
866
867 // Node in neighbour element
868 unsigned i0_neigh = 0;
869 unsigned i1_neigh = i1;
870 unsigned i2_neigh = i2;
871 unsigned jnod_local_neigh =
872 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
873
874 // Check:
875 for (unsigned i = 0; i < 3; i++)
876 {
877 double error = std::fabs(
878 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
879 finite_element_pt(ielem_neigh)
880 ->node_pt(jnod_local_neigh)
881 ->x(i));
882 if (error > node_kill_tol)
883 {
884 oomph_info << "Error in node killing for i " << i << " "
885 << error << std::endl;
886 stopit = true;
887 }
888 }
889
890 // Kill node
891 delete finite_element_pt(ielem)->node_pt(jnod_local);
892 killed = true;
893
894 // Set pointer to neighbour:
895 finite_element_pt(ielem)->node_pt(jnod_local) =
896 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
897 }
898
899 // Duplicate node: kill and set pointer to other ring element
900 if ((i1 == 0) && (i0 != n_p - 1))
901 {
902 // Neighbour element
903 unsigned ielem_neigh = ielem - 3;
904
905 // Node in neighbour element
906 unsigned i0_neigh = 0;
907 unsigned i1_neigh = i0;
908 unsigned i2_neigh = i2;
909 unsigned jnod_local_neigh =
910 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
911
912 // Check:
913 for (unsigned i = 0; i < 3; i++)
914 {
915 double error = std::fabs(
916 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
917 finite_element_pt(ielem_neigh)
918 ->node_pt(jnod_local_neigh)
919 ->x(i));
920 if (error > node_kill_tol)
921 {
922 oomph_info << "Error in node killing for i " << i << " "
923 << error << std::endl;
924 stopit = true;
925 }
926 }
927
928 // Kill node
929 delete finite_element_pt(ielem)->node_pt(jnod_local);
930 killed = true;
931
932 // Set pointer to neighbour:
933 finite_element_pt(ielem)->node_pt(jnod_local) =
934 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
935 }
936
937
938 // No duplicate node: Copy across to mesh
939 if (!killed)
940 {
941 Node_pt[node_count] =
942 finite_element_pt(ielem)->node_pt(jnod_local);
943
944 // Set boundaries:
945
946 // Back: Boundary 0
947 if ((i2 == 0) && (ilayer == 0))
948 {
949 this->convert_to_boundary_node(Node_pt[node_count]);
950 add_boundary_node(0, Node_pt[node_count]);
951 }
952
953 // Front: Boundary 2
954 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
955 {
956 this->convert_to_boundary_node(Node_pt[node_count]);
957 add_boundary_node(2, Node_pt[node_count]);
958 }
959
960 // Tube boundary
961 if (i0 == 0)
962 {
963 this->convert_to_boundary_node(Node_pt[node_count]);
964 add_boundary_node(1, Node_pt[node_count]);
965
966 // Get axial boundary coordinate
967 zeta[0] =
968 centreline_limits[0] +
969 (double(ilayer) + double(i2) / double(n_p - 1)) *
970 (centreline_limits[1] - centreline_limits[0]) /
971 double(nlayer);
972
973 // Get azimuthal boundary coordinate
974 zeta[1] =
975 theta_positions[0] + 2.0 * pi +
976 double(i1) / double(n_p - 1) * 0.5 *
977 (theta_positions[3] - theta_positions[0] + 2.0 * pi);
978
979 Node_pt[node_count]->set_coordinates_on_boundary(1, zeta);
980 }
981
982 // Increment node counter
983 node_count++;
984 }
985 }
986 }
987 }
988 }
989 break;
990 }
991 }
992
993 // Terminate if there's been an error
994 if (stopit)
995 {
996 std::ostringstream error_stream;
997 error_stream << "Error in killing nodes\n"
998 << "The most probable cause is that the domain is not\n"
999 << "compatible with the mesh.\n"
1000 << "For the TubeMesh, the domain must be\n"
1001 << "topologically consistent with a quarter tube with a\n"
1002 << "non-curved centreline.\n";
1003 throw OomphLibError(
1004 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1005 }
1006
1007 // Setup boundary element lookup schemes
1008 setup_boundary_element_info();
1009 }
1010
1011} // namespace oomph
1012#endif
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::BrickMeshBase::setup_boundary_element_info
void setup_boundary_element_info()
Setup lookup schemes which establish whic elements are located next to mesh's boundaries (wrapper to ...
Definition brick_mesh.h:195
oomph::Domain::macro_element_pt
MacroElement * macro_element_pt(const unsigned &i)
Access to i-th macro element.
Definition domain.h:116
oomph::FiniteElement::construct_node
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:2513
oomph::FiniteElement::node_pt
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition elements.h:2179
oomph::FiniteElement::nnode_1d
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition elements.h:2222
oomph::GeomObject
/////////////////////////////////////////////////////////////////////
Definition geom_objects.h:101
oomph::MacroElement::macro_map
void macro_map(const Vector< double > &s, Vector< double > &r)
The mapping from local to global coordinates at the current time : r(s)
Definition macro_element.h:126
oomph::Mesh::add_boundary_node
void add_boundary_node(const unsigned &b, Node *const &node_pt)
Add a (pointer to) a node to the b-th boundary.
Definition mesh.cc:243
oomph::Mesh::Node_pt
Vector< Node * > Node_pt
Vector of pointers to nodes.
Definition mesh.h:183
oomph::Mesh::Boundary_coordinate_exists
std::vector< bool > Boundary_coordinate_exists
Vector of boolean data that indicates whether the boundary coordinates have been set for the boundary...
Definition mesh.h:190
oomph::Mesh::finite_element_pt
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition mesh.h:473
oomph::Mesh::set_nboundary
void set_nboundary(const unsigned &nbound)
Set the number of boundaries in the mesh.
Definition mesh.h:505
oomph::Mesh::node_pt
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition mesh.h:436
oomph::Mesh::convert_to_boundary_node
void convert_to_boundary_node(Node *&node_pt, const Vector< FiniteElement * > &finite_element_pt)
A function that upgrades an ordinary node to a boundary node We shouldn't ever really use this,...
Definition mesh.cc:2590
oomph::Mesh::Element_pt
Vector< GeneralisedElement * > Element_pt
Vector of pointers to generalised elements.
Definition mesh.h:186
oomph::Node
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition nodes.h:906
oomph::Node::x
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition nodes.h:1060
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:222
oomph::TAdvectionDiffusionReactionElement
//////////////////////////////////////////////////////////////////////
Definition Tadvection_diffusion_reaction_elements.h:66
oomph::TimeStepper
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition timesteppers.h:231
oomph::TubeDomain
Tube as a domain. The entire domain must be defined by a GeomObject with the following convention: ze...
Definition tube_domain.h:70
oomph::TubeMesh::TubeMesh
TubeMesh(GeomObject *wall_pt, const Vector< double > &centreline_limits, const Vector< double > &theta_positions, const Vector< double > &radius_box, const unsigned &nlayer, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that specifies the volume, start and end coordinates fo...
Definition tube_mesh.template.cc:40
oomph::TubeMesh::Domain_pt
TubeDomain * Domain_pt
Pointer to domain.
Definition tube_mesh.template.h:99
oomph::TubeMesh::volume_pt
GeomObject *& volume_pt()
Access function to GeomObject representing wall.
Definition tube_mesh.template.h:80
oomph::MathematicalConstants::Pi
const double Pi
50 digits from maple
Definition oomph_utilities.h:157
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition advection_diffusion_elements.cc:30
oomph::oomph_info
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...
Definition oomph_definitions.cc:319
tube_mesh.template.h