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