quarter_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_QUARTER_TUBE_MESH_TEMPLATE_HEADER
27#define OOMPH_QUARTER_TUBE_MESH_TEMPLATE_HEADER
28
29#ifndef OOMPH_QUARTER_TUBE_MESH_HEADER
30#error __FILE__ should only be included from quarter_tube_mesh.h.
31#endif // OOMPH_QUARTER_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 boundary 3. Pass pointer to geometric object that
39 /// specifies the wall, start and end coordinates on the
40 /// geometric object, and the fraction along
41 /// which the dividing line is to be placed, and the timestepper.
42 /// Timestepper defaults to Static dummy timestepper.
43 //====================================================================
44 template<class ELEMENT>
45 QuarterTubeMesh<ELEMENT>::QuarterTubeMesh(GeomObject* wall_pt,
46 const Vector<double>& xi_lo,
47 const double& fract_mid,
48 const Vector<double>& xi_hi,
49 const unsigned& nlayer,
50 TimeStepper* time_stepper_pt)
51 : Wall_pt(wall_pt), Xi_lo(xi_lo), Fract_mid(fract_mid), Xi_hi(xi_hi)
52 {
53 // Mesh can only be built with 3D Qelements.
54 MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(3);
55
56 // Build macro element-based domain
57 Domain_pt = new QuarterTubeDomain(wall_pt, xi_lo, fract_mid, xi_hi, nlayer);
58
59 // Set the number of boundaries
60 set_nboundary(5);
61
62 // We have only bothered to parametrise boundary 3
63 set_boundary_coordinate_exists(3);
64
65 // Allocate the store for the elements
66 unsigned nelem = 3 * nlayer;
67 Element_pt.resize(nelem);
68
69 // Create dummy element so we can determine the number of nodes
70 ELEMENT* dummy_el_pt = new ELEMENT;
71
72 // Read out the number of linear points in the element
73 unsigned n_p = dummy_el_pt->nnode_1d();
74
75 // Kill the element
76 delete dummy_el_pt;
77
78 // Can now allocate the store for the nodes
79 unsigned nnodes_total =
80 (n_p * n_p + (n_p - 1) * n_p + (n_p - 1) * (n_p - 1)) *
81 (1 + nlayer * (n_p - 1));
82 Node_pt.resize(nnodes_total);
83
84 Vector<double> s(3);
85 Vector<double> r(3);
86
87 // Storage for the intrinsic boundary coordinate
88 Vector<double> zeta(2);
89
90 // Loop over elements and create all nodes
91 for (unsigned ielem = 0; ielem < nelem; ielem++)
92 {
93 // Create element
94 Element_pt[ielem] = new ELEMENT;
95
96 // Loop over rows in z/s_2-direction
97 for (unsigned i2 = 0; i2 < n_p; i2++)
98 {
99 // Loop over rows in y/s_1-direction
100 for (unsigned i1 = 0; i1 < n_p; i1++)
101 {
102 // Loop over rows in x/s_0-direction
103 for (unsigned i0 = 0; i0 < n_p; i0++)
104 {
105 // Local node number
106 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
107
108 // Create the node
109 Node* node_pt = finite_element_pt(ielem)->construct_node(
110 jnod_local, time_stepper_pt);
111
112 // Set the position of the node from macro element mapping
113 s[0] = -1.0 + 2.0 * double(i0) / double(n_p - 1);
114 s[1] = -1.0 + 2.0 * double(i1) / double(n_p - 1);
115 s[2] = -1.0 + 2.0 * double(i2) / double(n_p - 1);
116 Domain_pt->macro_element_pt(ielem)->macro_map(s, r);
117
118 node_pt->x(0) = r[0];
119 node_pt->x(1) = r[1];
120 node_pt->x(2) = r[2];
121 }
122 }
123 }
124 }
125
126 // Initialise number of global nodes
127 unsigned node_count = 0;
128
129 // Tolerance for node killing:
130 double node_kill_tol = 1.0e-12;
131
132 // Check for error in node killing
133 bool stopit = false;
134
135 // Loop over elements
136 for (unsigned ielem = 0; ielem < nelem; ielem++)
137 {
138 // Which layer?
139 unsigned ilayer = unsigned(ielem / 3);
140
141 // Which macro element?
142 switch (ielem % 3)
143 {
144 // Macro element 0: Central box
145 //-----------------------------
146 case 0:
147
148 // Loop over rows in z/s_2-direction
149 for (unsigned i2 = 0; i2 < n_p; i2++)
150 {
151 // Loop over rows in y/s_1-direction
152 for (unsigned i1 = 0; i1 < n_p; i1++)
153 {
154 // Loop over rows in x/s_0-direction
155 for (unsigned i0 = 0; i0 < n_p; i0++)
156 {
157 // Local node number
158 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
159
160 // Has the node been killed?
161 bool killed = false;
162
163 // First layer of all nodes in s_2 direction gets killed
164 // and re-directed to nodes in previous element layer
165 if ((i2 == 0) && (ilayer > 0))
166 {
167 // Neighbour element
168 unsigned ielem_neigh = ielem - 3;
169
170 // Node in neighbour element
171 unsigned i0_neigh = i0;
172 unsigned i1_neigh = i1;
173 unsigned i2_neigh = n_p - 1;
174 unsigned jnod_local_neigh =
175 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
176
177 // Check:
178 for (unsigned i = 0; i < 3; i++)
179 {
180 double error = std::fabs(
181 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
182 finite_element_pt(ielem_neigh)
183 ->node_pt(jnod_local_neigh)
184 ->x(i));
185 if (error > node_kill_tol)
186 {
187 oomph_info << "Error in node killing for i " << i << " "
188 << error << std::endl;
189 stopit = true;
190 }
191 }
192
193 // Kill node
194 delete finite_element_pt(ielem)->node_pt(jnod_local);
195 killed = true;
196
197 // Set pointer to neighbour:
198 finite_element_pt(ielem)->node_pt(jnod_local) =
199 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
200 }
201
202 // No duplicate node: Copy across to mesh
203 if (!killed)
204 {
205 Node_pt[node_count] =
206 finite_element_pt(ielem)->node_pt(jnod_local);
207
208 // Set boundaries:
209
210 // Back: Boundary 0
211 if ((i2 == 0) && (ilayer == 0))
212 {
213 this->convert_to_boundary_node(Node_pt[node_count]);
214 add_boundary_node(0, Node_pt[node_count]);
215 }
216
217 // Front: Boundary 4
218 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
219 {
220 this->convert_to_boundary_node(Node_pt[node_count]);
221 add_boundary_node(4, Node_pt[node_count]);
222 }
223
224 // Left symmetry plane: Boundary 1
225 if (i0 == 0)
226 {
227 this->convert_to_boundary_node(Node_pt[node_count]);
228 add_boundary_node(1, Node_pt[node_count]);
229 }
230
231 // Bottom symmetry plane: Boundary 2
232 if (i1 == 0)
233 {
234 this->convert_to_boundary_node(Node_pt[node_count]);
235 add_boundary_node(2, Node_pt[node_count]);
236 }
237
238 // Increment node counter
239 node_count++;
240 }
241 }
242 }
243 }
244
245 break;
246
247 // Macro element 1: Lower right box
248 //---------------------------------
249 case 1:
250
251 // Loop over rows in z/s_2-direction
252 for (unsigned i2 = 0; i2 < n_p; i2++)
253 {
254 // Loop over rows in y/s_1-direction
255 for (unsigned i1 = 0; i1 < n_p; i1++)
256 {
257 // Loop over rows in x/s_0-direction
258 for (unsigned i0 = 0; i0 < n_p; i0++)
259 {
260 // Local node number
261 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
262
263 // Has the node been killed?
264 bool killed = false;
265
266 // First layer of all nodes in s_2 direction gets killed
267 // and re-directed to nodes in previous element layer
268 if ((i2 == 0) && (ilayer > 0))
269 {
270 // Neighbour element
271 unsigned ielem_neigh = ielem - 3;
272
273 // Node in neighbour element
274 unsigned i0_neigh = i0;
275 unsigned i1_neigh = i1;
276 unsigned i2_neigh = n_p - 1;
277 unsigned jnod_local_neigh =
278 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
279
280 // Check:
281 for (unsigned i = 0; i < 3; i++)
282 {
283 double error = std::fabs(
284 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
285 finite_element_pt(ielem_neigh)
286 ->node_pt(jnod_local_neigh)
287 ->x(i));
288 if (error > node_kill_tol)
289 {
290 oomph_info << "Error in node killing for i " << i << " "
291 << error << std::endl;
292 stopit = true;
293 }
294 }
295
296 // Kill node
297 delete finite_element_pt(ielem)->node_pt(jnod_local);
298 killed = true;
299
300 // Set pointer to neighbour:
301 finite_element_pt(ielem)->node_pt(jnod_local) =
302 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
303 }
304 // Not in first layer:
305 else
306 {
307 // Duplicate node: kill and set pointer to central element
308 if (i0 == 0)
309 {
310 // Neighbour element
311 unsigned ielem_neigh = ielem - 1;
312
313 // Node in neighbour element
314 unsigned i0_neigh = n_p - 1;
315 unsigned i1_neigh = i1;
316 unsigned i2_neigh = i2;
317 unsigned jnod_local_neigh =
318 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
319
320 // Check:
321 for (unsigned i = 0; i < 3; i++)
322 {
323 double error = std::fabs(
324 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
325 finite_element_pt(ielem_neigh)
326 ->node_pt(jnod_local_neigh)
327 ->x(i));
328 if (error > node_kill_tol)
329 {
330 oomph_info << "Error in node killing for i " << i << " "
331 << error << std::endl;
332 stopit = true;
333 }
334 }
335
336 // Kill node
337 delete finite_element_pt(ielem)->node_pt(jnod_local);
338 killed = true;
339
340 // Set pointer to neighbour:
341 finite_element_pt(ielem)->node_pt(jnod_local) =
342 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
343 }
344 }
345
346 // No duplicate node: Copy across to mesh
347 if (!killed)
348 {
349 Node_pt[node_count] =
350 finite_element_pt(ielem)->node_pt(jnod_local);
351
352 // Set boundaries:
353
354 // Back: Boundary 0
355 if ((i2 == 0) && (ilayer == 0))
356 {
357 this->convert_to_boundary_node(Node_pt[node_count]);
358 add_boundary_node(0, Node_pt[node_count]);
359 }
360
361 // Front: Boundary 4
362 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
363 {
364 this->convert_to_boundary_node(Node_pt[node_count]);
365 add_boundary_node(4, Node_pt[node_count]);
366 }
367
368 // Bottom symmetry plane: Boundary 2
369 if (i1 == 0)
370 {
371 this->convert_to_boundary_node(Node_pt[node_count]);
372 add_boundary_node(2, Node_pt[node_count]);
373 }
374
375 // Tube wall: Boundary 3
376 if (i0 == n_p - 1)
377 {
378 this->convert_to_boundary_node(Node_pt[node_count]);
379 add_boundary_node(3, Node_pt[node_count]);
380
381 // Get axial boundary coordinate
382 zeta[0] = Xi_lo[0] +
383 (double(ilayer) + double(i2) / double(n_p - 1)) *
384 (Xi_hi[0] - Xi_lo[0]) / double(nlayer);
385
386 // Get azimuthal boundary coordinate
387 zeta[1] = Xi_lo[1] + double(i1) / double(n_p - 1) * 0.5 *
388 (Xi_hi[1] - Xi_lo[1]);
389
390 Node_pt[node_count]->set_coordinates_on_boundary(3, zeta);
391 }
392
393 // Increment node counter
394 node_count++;
395 }
396 }
397 }
398 }
399
400 break;
401
402 // Macro element 2: Top left box
403 //--------------------------------
404 case 2:
405
406 // Loop over rows in z/s_2-direction
407 for (unsigned i2 = 0; i2 < n_p; i2++)
408 {
409 // Loop over rows in y/s_1-direction
410 for (unsigned i1 = 0; i1 < n_p; i1++)
411 {
412 // Loop over rows in x/s_0-direction
413 for (unsigned i0 = 0; i0 < n_p; i0++)
414 {
415 // Local node number
416 unsigned jnod_local = i0 + i1 * n_p + i2 * n_p * n_p;
417
418 // Has the node been killed?
419 bool killed = false;
420
421 // First layer of all nodes in s_2 direction gets killed
422 // and re-directed to nodes in previous element layer
423 if ((i2 == 0) && (ilayer > 0))
424 {
425 // Neighbour element
426 unsigned ielem_neigh = ielem - 3;
427
428 // Node in neighbour element
429 unsigned i0_neigh = i0;
430 unsigned i1_neigh = i1;
431 unsigned i2_neigh = n_p - 1;
432 unsigned jnod_local_neigh =
433 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
434
435 // Check:
436 for (unsigned i = 0; i < 3; i++)
437 {
438 double error = std::fabs(
439 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
440 finite_element_pt(ielem_neigh)
441 ->node_pt(jnod_local_neigh)
442 ->x(i));
443 if (error > node_kill_tol)
444 {
445 oomph_info << "Error in node killing for i " << i << " "
446 << error << std::endl;
447 stopit = true;
448 }
449 }
450
451 // Kill node
452 delete finite_element_pt(ielem)->node_pt(jnod_local);
453 killed = true;
454
455 // Set pointer to neighbour:
456 finite_element_pt(ielem)->node_pt(jnod_local) =
457 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
458 }
459 // Not in first layer:
460 else
461 {
462 // Duplicate node: kill and set pointer to node in bottom
463 // right element
464 if (i0 == n_p - 1)
465 {
466 // Neighbour element
467 unsigned ielem_neigh = ielem - 1;
468
469 // Node in neighbour element
470 unsigned i0_neigh = i1;
471 unsigned i1_neigh = n_p - 1;
472 unsigned i2_neigh = i2;
473 unsigned jnod_local_neigh =
474 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
475
476 // Check:
477 for (unsigned i = 0; i < 3; i++)
478 {
479 double error = std::fabs(
480 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
481 finite_element_pt(ielem_neigh)
482 ->node_pt(jnod_local_neigh)
483 ->x(i));
484 if (error > node_kill_tol)
485 {
486 oomph_info << "Error in node killing for i " << i << " "
487 << error << std::endl;
488 stopit = true;
489 }
490 }
491
492 // Kill node
493 delete finite_element_pt(ielem)->node_pt(jnod_local);
494 killed = true;
495
496 // Set pointer to neighbour:
497 finite_element_pt(ielem)->node_pt(jnod_local) =
498 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
499 }
500
501 // Duplicate node: kill and set pointer to central element
502 if ((i1 == 0) && (i0 != n_p - 1))
503 {
504 // Neighbour element
505 unsigned ielem_neigh = ielem - 2;
506
507 // Node in neighbour element
508 unsigned i0_neigh = i0;
509 unsigned i1_neigh = n_p - 1;
510 unsigned i2_neigh = i2;
511 unsigned jnod_local_neigh =
512 i0_neigh + i1_neigh * n_p + i2_neigh * n_p * n_p;
513
514 // Check:
515 for (unsigned i = 0; i < 3; i++)
516 {
517 double error = std::fabs(
518 finite_element_pt(ielem)->node_pt(jnod_local)->x(i) -
519 finite_element_pt(ielem_neigh)
520 ->node_pt(jnod_local_neigh)
521 ->x(i));
522 if (error > node_kill_tol)
523 {
524 oomph_info << "Error in node killing for i " << i << " "
525 << error << std::endl;
526 stopit = true;
527 }
528 }
529
530 // Kill node
531 delete finite_element_pt(ielem)->node_pt(jnod_local);
532 killed = true;
533
534 // Set pointer to neighbour:
535 finite_element_pt(ielem)->node_pt(jnod_local) =
536 finite_element_pt(ielem_neigh)->node_pt(jnod_local_neigh);
537 }
538
539 // No duplicate node: Copy across to mesh
540 if (!killed)
541 {
542 Node_pt[node_count] =
543 finite_element_pt(ielem)->node_pt(jnod_local);
544
545 // Set boundaries:
546
547 // Back: Boundary 0
548 if ((i2 == 0) && (ilayer == 0))
549 {
550 this->convert_to_boundary_node(Node_pt[node_count]);
551 add_boundary_node(0, Node_pt[node_count]);
552 }
553
554 // Front: Boundary 4
555 if ((i2 == n_p - 1) && (ilayer == nlayer - 1))
556 {
557 this->convert_to_boundary_node(Node_pt[node_count]);
558 add_boundary_node(4, Node_pt[node_count]);
559 }
560
561 // Left symmetry plane: Boundary 1
562 if (i0 == 0)
563 {
564 this->convert_to_boundary_node(Node_pt[node_count]);
565 add_boundary_node(1, Node_pt[node_count]);
566 }
567
568 // Tube wall: Boundary 3
569 if (i1 == n_p - 1)
570 {
571 this->convert_to_boundary_node(Node_pt[node_count]);
572 add_boundary_node(3, Node_pt[node_count]);
573
574 // Get axial boundary coordinate
575 zeta[0] =
576 Xi_lo[0] +
577 (double(ilayer) + double(i2) / double(n_p - 1)) *
578 (Xi_hi[0] - Xi_lo[0]) / double(nlayer);
579
580 // Get azimuthal boundary coordinate
581 zeta[1] = Xi_hi[1] - double(i0) / double(n_p - 1) * 0.5 *
582 (Xi_hi[1] - Xi_lo[1]);
583
584 Node_pt[node_count]->set_coordinates_on_boundary(3, zeta);
585 }
586
587 // Increment node counter
588 node_count++;
589 }
590 }
591 }
592 }
593 }
594
595 break;
596 }
597 }
598
599 // Terminate if there's been an error
600 if (stopit)
601 {
602 std::ostringstream error_stream;
603 error_stream << "Error in killing nodes\n"
604 << "The most probable cause is that the domain is not\n"
605 << "compatible with the mesh.\n"
606 << "For the QuarterTubeMesh, the domain must be\n"
607 << "topologically consistent with a quarter tube with a\n"
608 << "non-curved centreline.\n";
609 throw OomphLibError(
610 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
611 }
612
613 // Setup boundary element lookup schemes
614 setup_boundary_element_info();
615 }
616
617 ///////////////////////////////////////////////////////////////////////
618 ///////////////////////////////////////////////////////////////////////
619 // Algebraic-mesh-version of RefineableQuarterTubeMesh
620 ///////////////////////////////////////////////////////////////////////
621 ///////////////////////////////////////////////////////////////////////
622
623 //======================================================================
624 /// Setup algebraic node update data, based on 3 regions, each
625 /// undergoing a different node update strategy. These regions are
626 /// defined by the three MacroElements in each of the nlayer slices
627 /// of the QuarterTubeDomain used to build the mesh.
628 /// The Mesh is suspended from the `wall' GeomObject pointed to
629 /// by wall_pt. The lower right edge of the mesh is located at the
630 /// wall's coordinate xi[1]==xi_lo[1], the upper left edge at
631 /// xi[1]=xi_hi[1], i.e. a view looking down the tube length.
632 /// The dividing line between the two outer regions is located
633 /// at the fraction fract_mid between these two coordinates.
634 /// Node updating strategy:
635 /// - the starting cross sectional shape along the tube length is
636 /// assumed to be uniform
637 /// - the cross sectional shape of the central region remains
638 /// rectangular; the position of its top right corner is located
639 /// at a fixed fraction of the starting width and height of the
640 /// domain measured at xi[1]==xi_lo[1] and xi[1]==xi_hi[1]
641 /// respectively. Nodes in this region are located at fixed
642 /// horizontal and vertical fractions of the region.
643 /// - Nodes in the two outer regions (bottom right and top left)
644 /// are located on straight lines running from the edges of the
645 /// central region to the outer wall.
646 //======================================================================
647 template<class ELEMENT>
648 void AlgebraicRefineableQuarterTubeMesh<
649 ELEMENT>::setup_algebraic_node_update()
650 {
651#ifdef PARANOID
652 /// Pointer to first algebraic element in central region
653 AlgebraicElementBase* algebraic_element_pt =
654 dynamic_cast<AlgebraicElementBase*>(Mesh::element_pt(0));
655
656 if (algebraic_element_pt == 0)
657 {
658 std::ostringstream error_message;
659 error_message
660 << "Element in AlgebraicRefineableQuarterTubeMesh must be\n ";
661 error_message << "derived from AlgebraicElementBase\n";
662 error_message << "but it is of type: "
663 << typeid(Mesh::element_pt(0)).name() << std::endl;
664 std::string function_name = "AlgebraicRefineableQuarterTubeMesh::";
665 function_name += "setup_algebraic_node_update()";
666 throw OomphLibError(
667 error_message.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
668 }
669#endif
670
671 // Find number of nodes in an element from the zeroth element
672 unsigned nnodes_3d = Mesh::finite_element_pt(0)->nnode();
673
674 // also find number of nodes in 1d line and 2d slice
675 unsigned nnodes_1d = Mesh::finite_element_pt(0)->nnode_1d();
676 unsigned nnodes_2d = nnodes_1d * nnodes_1d;
677
678 // find node number of a top-left and a bottom-right node in an element
679 // (orientation: looking down tube)
680 unsigned tl_node = nnodes_2d - nnodes_1d;
681 unsigned br_node = nnodes_1d - 1;
682
683 // find x & y distances to top-right node in element 0 - this is the same
684 // node as the top-left node of element 1
685 double x_c_element = Mesh::finite_element_pt(1)->node_pt(tl_node)->x(0);
686 double y_c_element = Mesh::finite_element_pt(1)->node_pt(tl_node)->x(1);
687
688 // Get x-distance to bottom-right edge of wall, i.e. coord of node
689 // at bottom-right of bottom-right of element 1
690 double x_wall = Mesh::finite_element_pt(1)->node_pt(br_node)->x(0);
691
692 // Get y-distance to top-left edge of wall, i.e. coord of node
693 // at top-left of element 2
694 double y_wall = Mesh::finite_element_pt(2)->node_pt(tl_node)->x(1);
695
696 // Establish fractional widths in central region
697 Lambda_x = Centre_box_size * x_c_element / x_wall;
698 Lambda_y = Centre_box_size * y_c_element / y_wall;
699
700 // how many elements are there?
701 unsigned nelements = Mesh::nelement();
702
703 // loop through the elements
704 for (unsigned e = 0; e < nelements; e++)
705 {
706 // get pointer to element
707 FiniteElement* el_pt = Mesh::finite_element_pt(e);
708
709 // set region id
710 unsigned region_id = e % 3;
711
712 // find the first node for which to set up node update info - must
713 // bypass the first nnodes_2d nodes after the first 3 elements
714 unsigned nstart = nnodes_2d;
715 if (e < 3)
716 {
717 nstart = 0;
718 }
719
720 // loop through the nodes,
721 for (unsigned n = nstart; n < nnodes_3d; n++)
722 {
723 // find z coordinate of node
724 // NOTE: to implement axial spacing replace z by z_spaced where
725 // z_spaced=axial_spacing_fct(z) when finding the GeomObjects
726 // and local coords below
727 // BUT store z as the third reference value since this is the value
728 // required by update_node_update()
729 double z = el_pt->node_pt(n)->x(2);
730
731 // Find wall GeomObject and the GeomObject local coordinates
732 // at bottom-right edge of wall
733 Vector<double> xi(2);
734 xi[0] = z;
735 xi[1] = QuarterTubeMesh<ELEMENT>::Xi_lo[1];
736 GeomObject* obj_br_pt;
737 Vector<double> s_br(2);
738 this->Wall_pt->locate_zeta(xi, obj_br_pt, s_br);
739
740 // Find wall GeomObject/local coordinate
741 // at top-left edge of wall
742 xi[1] = QuarterTubeMesh<ELEMENT>::Xi_hi[1];
743 GeomObject* obj_tl_pt;
744 Vector<double> s_tl(2);
745 this->Wall_pt->locate_zeta(xi, obj_tl_pt, s_tl);
746
747 // Element 0: central region
748 //--------------------------
749 if (region_id == 0)
750 {
751 // Nodal coordinates in x and y direction
752 double x = el_pt->node_pt(n)->x(0);
753 double y = el_pt->node_pt(n)->x(1);
754
755 // The update function requires two geometric objects
756 Vector<GeomObject*> geom_object_pt(2);
757
758 // Wall GeomObject at bottom right end of wall mesh:
759 geom_object_pt[0] = obj_br_pt;
760
761 // Wall GeomObject at top left end of wall mesh:
762 geom_object_pt[1] = obj_tl_pt;
763
764 // The update function requires seven parameters:
765 Vector<double> ref_value(7);
766
767 // First reference value: fractional x-position inside region
768 ref_value[0] = x / x_c_element;
769
770 // Second cfractional y-position inside region
771 ref_value[1] = y / y_c_element;
772
773 // Third reference value: initial z coordinate of node
774 ref_value[2] = z;
775
776 // Fourth and fifth reference values:
777 // local coordinate in GeomObject at bottom-right of wall.
778 // Note: must recompute this reference for new nodes
779 // since values interpolated from parent nodes will
780 // be wrong (this is true for all local coordinates
781 // within GeomObjects)
782 ref_value[3] = s_br[0];
783 ref_value[4] = s_br[1];
784
785 // Sixth and seventh reference values:
786 // local coordinate in GeomObject at top-left of wall.
787 ref_value[5] = s_tl[0];
788 ref_value[6] = s_tl[1];
789
790 // Setup algebraic update for node: Pass update information
791 static_cast<AlgebraicNode*>(el_pt->node_pt(n))
792 ->add_node_update_info(Central_region, // enumerated ID
793 this, // mesh
794 geom_object_pt, // vector of geom objects
795 ref_value); // vector of ref. values
796 }
797
798 // Element 1: Bottom right region
799 //-------------------------------
800 else if (region_id == 1)
801 {
802 // Fractional distance between nodes
803 double fract = 1.0 / double(nnodes_1d - 1);
804
805 // Fraction in the s_0-direction
806 double rho_0 = fract * double(n % nnodes_1d);
807
808 // Fraction in the s_1-direction
809 double rho_1 = fract * double((n / nnodes_1d) % nnodes_1d);
810
811 // Find coord of reference point on wall
812 xi[1] = QuarterTubeMesh<ELEMENT>::Xi_lo[1] +
813 rho_1 * QuarterTubeMesh<ELEMENT>::Fract_mid *
814 (QuarterTubeMesh<ELEMENT>::Xi_hi[1] -
815 QuarterTubeMesh<ELEMENT>::Xi_lo[1]);
816
817 // Identify wall GeomObject and local coordinate of
818 // reference point in GeomObject
819 GeomObject* obj_wall_pt;
820 Vector<double> s_wall(2);
821 this->Wall_pt->locate_zeta(xi, obj_wall_pt, s_wall);
822
823 // The update function requires three geometric objects
824 Vector<GeomObject*> geom_object_pt(3);
825
826 // Wall element at bottom-right end of wall mesh:
827 geom_object_pt[0] = obj_br_pt;
828
829 // Wall element at top left end of wall mesh:
830 geom_object_pt[1] = obj_tl_pt;
831
832 // Wall element at that contians reference point:
833 geom_object_pt[2] = obj_wall_pt;
834
835 // The update function requires nine parameters:
836 Vector<double> ref_value(9);
837
838 // First reference value: fractional s0-position inside region
839 ref_value[0] = rho_0;
840
841 // Second reference value: fractional s1-position inside region
842 ref_value[1] = rho_1;
843
844 // Thrid reference value: initial z coord of node
845 ref_value[2] = z;
846
847 // Fourth and fifth reference values:
848 // local coordinates at bottom-right of wall.
849 ref_value[3] = s_br[0];
850 ref_value[4] = s_br[1];
851
852 // Sixth and seventh reference values:
853 // local coordinates at top-left of wall.
854 ref_value[5] = s_tl[0];
855 ref_value[6] = s_tl[1];
856
857 // Eighth and ninth reference values:
858 // local coordinate of reference point.
859 ref_value[7] = s_wall[0];
860 ref_value[8] = s_wall[1];
861
862 // Setup algebraic update for node: Pass update information
863 static_cast<AlgebraicNode*>(el_pt->node_pt(n))
864 ->add_node_update_info(Lower_right_region, // enumerated ID
865 this, // mesh
866 geom_object_pt, // vector of geom objects
867 ref_value); // vector of ref. vals
868 }
869
870 // Element 2: Top left region
871 //---------------------------
872 else if (region_id == 2)
873 {
874 // Fractional distance between nodes
875 double fract = 1.0 / double(nnodes_1d - 1);
876
877 // Fraction in the s_0-direction
878 double rho_0 = fract * double(n % nnodes_1d);
879
880 // Fraction in the s_1-direction
881 double rho_1 = fract * double((n / nnodes_1d) % nnodes_1d);
882
883 // Find coord of reference point on wall
884 xi[1] = QuarterTubeMesh<ELEMENT>::Xi_hi[1] +
885 rho_0 * (1.0 - QuarterTubeMesh<ELEMENT>::Fract_mid) *
886 (QuarterTubeMesh<ELEMENT>::Xi_lo[1] -
887 QuarterTubeMesh<ELEMENT>::Xi_hi[1]);
888
889 // Identify GeomObject and local coordinate at
890 // reference point on wall
891 GeomObject* obj_wall_pt;
892 Vector<double> s_wall(2);
893 this->Wall_pt->locate_zeta(xi, obj_wall_pt, s_wall);
894
895 // The update function requires three geometric objects
896 Vector<GeomObject*> geom_object_pt(3);
897
898 // Wall element at bottom-right of wall mesh:
899 geom_object_pt[0] = obj_br_pt;
900
901 // Wall element at top-left of wall mesh:
902 geom_object_pt[1] = obj_tl_pt;
903
904 // Wall element contianing reference point:
905 geom_object_pt[2] = obj_wall_pt;
906
907 // The update function requires nine parameters:
908 Vector<double> ref_value(9);
909
910 // First reference value: fractional s0-position inside region
911 ref_value[0] = rho_0;
912
913 // Second reference value: fractional s1-position inside region
914 ref_value[1] = rho_1;
915
916 // Third reference value: initial z coord
917 ref_value[2] = z;
918
919 // Fourth and fifth reference values:
920 // local coordinates in bottom-right GeomObject
921 ref_value[3] = s_br[0];
922 ref_value[4] = s_br[1];
923
924 // Sixth and seventh reference values:
925 // local coordinates in top-left GeomObject
926 ref_value[5] = s_tl[0];
927 ref_value[6] = s_tl[1];
928
929 // Eighth and ninth reference values:
930 // local coordinates at reference point
931 ref_value[7] = s_wall[0];
932 ref_value[8] = s_wall[1];
933
934 // Setup algebraic update for node: Pass update information
935 static_cast<AlgebraicNode*>(el_pt->node_pt(n))
936 ->add_node_update_info(Upper_left_region, // Enumerated ID
937 this, // mesh
938 geom_object_pt, // vector of geom objects
939 ref_value); // vector of ref. vals
940 }
941 }
942 }
943 }
944
945 //======================================================================
946 /// Algebraic update function: Update in central region according
947 /// to wall shape at time level t (t=0: present; t>0: previous)
948 //======================================================================
949 template<class ELEMENT>
950 void AlgebraicRefineableQuarterTubeMesh<ELEMENT>::node_update_central_region(
951 const unsigned& t, AlgebraicNode*& node_pt)
952 {
953#ifdef PARANOID
954 // We're updating the nodal positions (!) at time level t
955 // and determine them by evaluating the wall GeomObject's
956 // position at that time level. I believe this only makes sense
957 // if the t-th history value in the positional timestepper
958 // actually represents previous values (rather than some
959 // generalised quantity). Hence if this function is called with
960 // t>nprev_values(), we issue a warning and terminate the execution.
961 // It *might* be possible that the code still works correctly
962 // even if this condition is violated (e.g. if the GeomObject's
963 // position() function returns the appropriate "generalised"
964 // position value that is required by the timestepping scheme but it's
965 // probably worth flagging this up and forcing the user to manually switch
966 // off this warning if he/she is 100% sure that this is kosher.
967 if (t > node_pt->position_time_stepper_pt()->nprev_values())
968 {
969 std::string error_message =
970 "Trying to update the nodal position at a time level\n";
971 error_message += "beyond the number of previous values in the nodes'\n";
972 error_message += "position timestepper. This seems highly suspect!\n";
973 error_message += "If you're sure the code behaves correctly\n";
974 error_message += "in your application, remove this warning \n";
975 error_message += "or recompile with PARNOID switched off.\n";
976
977 std::string function_name = "AlgebraicRefineableQuarterTubeMesh::";
978 function_name += "node_update_central_region()",
979 throw OomphLibError(
980 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
981 }
982#endif
983
984 // Extract references for update in central region by copy construction
985 Vector<double> ref_value(node_pt->vector_ref_value(Central_region));
986
987 // Extract geometric objects for update in central region by copy
988 // construction
989 Vector<GeomObject*> geom_object_pt(
990 node_pt->vector_geom_object_pt(Central_region));
991
992 // First reference value: fractional x-position of node inside region
993 double rho_x = ref_value[0];
994
995 // Second reference value: fractional y-position of node inside region
996 double rho_y = ref_value[1];
997
998 // Wall position in bottom right corner:
999
1000 // Pointer to GeomObject
1001 GeomObject* obj_br_pt = geom_object_pt[0];
1002
1003 // Local coordinate at bottom-right reference point:
1004 Vector<double> s_br(2);
1005 s_br[0] = ref_value[3];
1006 s_br[1] = ref_value[4];
1007
1008 // Get wall position
1009 unsigned n_dim = obj_br_pt->ndim();
1010 Vector<double> r_br(n_dim);
1011 obj_br_pt->position(t, s_br, r_br);
1012
1013 // Wall position in top left corner:
1014
1015 // Pointer to GeomObject
1016 GeomObject* obj_tl_pt = geom_object_pt[1];
1017
1018 // Local coordinate:
1019 Vector<double> s_tl(2);
1020 s_tl[0] = ref_value[5];
1021 s_tl[1] = ref_value[6];
1022
1023 // Get wall position
1024 Vector<double> r_tl(n_dim);
1025 obj_tl_pt->position(t, s_tl, r_tl);
1026
1027 // Assign new nodal coordinate
1028 node_pt->x(t, 0) = r_br[0] * Lambda_x * rho_x;
1029 node_pt->x(t, 1) = r_tl[1] * Lambda_y * rho_y;
1030 node_pt->x(t, 2) = (r_br[2] + r_tl[2]) * 0.5;
1031 }
1032
1033 //====================================================================
1034 /// Algebraic update function: Update in lower-right region according
1035 /// to wall shape at time level t (t=0: present; t>0: previous)
1036 //====================================================================
1037 template<class ELEMENT>
1038 void AlgebraicRefineableQuarterTubeMesh<
1039 ELEMENT>::node_update_lower_right_region(const unsigned& t,
1040 AlgebraicNode*& node_pt)
1041 {
1042#ifdef PARANOID
1043 // We're updating the nodal positions (!) at time level t
1044 // and determine them by evaluating the wall GeomObject's
1045 // position at that gime level. I believe this only makes sense
1046 // if the t-th history value in the positional timestepper
1047 // actually represents previous values (rather than some
1048 // generalised quantity). Hence if this function is called with
1049 // t>nprev_values(), we issue a warning and terminate the execution.
1050 // It *might* be possible that the code still works correctly
1051 // even if this condition is violated (e.g. if the GeomObject's
1052 // position() function returns the appropriate "generalised"
1053 // position value that is required by the timestepping scheme but it's
1054 // probably worth flagging this up and forcing the user to manually switch
1055 // off this warning if he/she is 100% sure that this is kosher.
1056 if (t > node_pt->position_time_stepper_pt()->nprev_values())
1057 {
1058 std::string error_message =
1059 "Trying to update the nodal position at a time level";
1060 error_message += "beyond the number of previous values in the nodes'";
1061 error_message += "position timestepper. This seems highly suspect!";
1062 error_message += "If you're sure the code behaves correctly";
1063 error_message += "in your application, remove this warning ";
1064 error_message += "or recompile with PARNOID switched off.";
1065
1066 std::string function_name = "AlgebraicRefineableQuarterTubeSectorMesh::";
1067 function_name += "node_update_lower_right_region()",
1068 throw OomphLibError(
1069 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1070 }
1071#endif
1072
1073 // Extract references for update in lower-right region by copy construction
1074 Vector<double> ref_value(node_pt->vector_ref_value(Lower_right_region));
1075
1076 // Extract geometric objects for update in lower-right region
1077 // by copy construction
1078 Vector<GeomObject*> geom_object_pt(
1079 node_pt->vector_geom_object_pt(Lower_right_region));
1080
1081 // First reference value: fractional s0-position of node inside region
1082 double rho_0 = ref_value[0];
1083
1084 // Use s_squashed to modify fractional s0 position
1085 rho_0 = this->Domain_pt->s_squashed(rho_0);
1086
1087 // Second reference value: fractional s1-position of node inside region
1088 double rho_1 = ref_value[1];
1089
1090 // Wall position in bottom right corner:
1091
1092 // Pointer to GeomObject
1093 GeomObject* obj_br_pt = geom_object_pt[0];
1094
1095 // Local coordinate
1096 Vector<double> s_br(2);
1097 s_br[0] = ref_value[3];
1098 s_br[1] = ref_value[4];
1099
1100 // Eulerian dimension
1101 unsigned n_dim = obj_br_pt->ndim();
1102
1103 // Get wall position
1104 Vector<double> r_br(n_dim);
1105 obj_br_pt->position(t, s_br, r_br);
1106
1107 // Wall position in top left corner:
1108
1109 // Pointer to GeomObject
1110 GeomObject* obj_tl_pt = geom_object_pt[1];
1111
1112 // Local coordinate
1113 Vector<double> s_tl(2);
1114 s_tl[0] = ref_value[5];
1115 s_tl[1] = ref_value[6];
1116
1117 Vector<double> r_tl(n_dim);
1118 obj_tl_pt->position(t, s_tl, r_tl);
1119
1120 // Wall position at reference point:
1121
1122 // Pointer to GeomObject
1123 GeomObject* obj_wall_pt = geom_object_pt[2];
1124
1125 // Local coordinate
1126 Vector<double> s_wall(2);
1127 s_wall[0] = ref_value[7];
1128 s_wall[1] = ref_value[8];
1129
1130 Vector<double> r_wall(n_dim);
1131 obj_wall_pt->position(t, s_wall, r_wall);
1132
1133 // Position Vector to corner of the central region
1134 Vector<double> r_corner(n_dim);
1135 r_corner[0] = Lambda_x * r_br[0];
1136 r_corner[1] = Lambda_y * r_tl[1];
1137 r_corner[2] = (r_br[2] + r_tl[2]) * 0.5;
1138
1139 // Position Vector to left edge of region
1140 Vector<double> r_left(n_dim);
1141 r_left[0] = r_corner[0];
1142 r_left[1] = rho_1 * r_corner[1];
1143 r_left[2] = r_corner[2];
1144
1145 // Assign new nodal coordinate
1146 node_pt->x(t, 0) = r_left[0] + rho_0 * (r_wall[0] - r_left[0]);
1147 node_pt->x(t, 1) = r_left[1] + rho_0 * (r_wall[1] - r_left[1]);
1148 node_pt->x(t, 2) = r_left[2] + rho_0 * (r_wall[2] - r_left[2]);
1149 }
1150 //====================================================================
1151 /// Algebraic update function: Update in upper left region according
1152 /// to wall shape at time level t (t=0: present; t>0: previous)
1153 //====================================================================
1154 template<class ELEMENT>
1155 void AlgebraicRefineableQuarterTubeMesh<
1156 ELEMENT>::node_update_upper_left_region(const unsigned& t,
1157 AlgebraicNode*& node_pt)
1158 {
1159#ifdef PARANOID
1160 // We're updating the nodal positions (!) at time level t
1161 // and determine them by evaluating the wall GeomObject's
1162 // position at that gime level. I believe this only makes sense
1163 // if the t-th history value in the positional timestepper
1164 // actually represents previous values (rather than some
1165 // generalised quantity). Hence if this function is called with
1166 // t>nprev_values(), we issue a warning and terminate the execution.
1167 // It *might* be possible that the code still works correctly
1168 // even if this condition is violated (e.g. if the GeomObject's
1169 // position() function returns the appropriate "generalised"
1170 // position value that is required by the timestepping scheme but it's
1171 // probably worth flagging this up and forcing the user to manually switch
1172 // off this warning if he/she is 100% sure that this is kosher.
1173 if (t > node_pt->position_time_stepper_pt()->nprev_values())
1174 {
1175 std::string error_message =
1176 "Trying to update the nodal position at a time level";
1177 error_message += "beyond the number of previous values in the nodes'";
1178 error_message += "position timestepper. This seems highly suspect!";
1179 error_message += "If you're sure the code behaves correctly";
1180 error_message += "in your application, remove this warning ";
1181 error_message += "or recompile with PARNOID switched off.";
1182
1183 std::string function_name = "AlgebraicRefineableQuarterTubeMesh::";
1184 function_name += "node_update_upper_left_region()";
1185
1186 throw OomphLibError(
1187 error_message, OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
1188 }
1189#endif
1190
1191 // Extract references for update in upper-left region by copy construction
1192 Vector<double> ref_value(node_pt->vector_ref_value(Upper_left_region));
1193
1194 // Extract geometric objects for update in upper-left region
1195 // by copy construction
1196 Vector<GeomObject*> geom_object_pt(
1197 node_pt->vector_geom_object_pt(Upper_left_region));
1198
1199 // First reference value: fractional s0-position of node inside region
1200 double rho_0 = ref_value[0];
1201
1202 // Second reference value: fractional s1-position of node inside region
1203 double rho_1 = ref_value[1];
1204
1205 // Use s_squashed to modify fractional s1 position
1206 rho_1 = this->Domain_pt->s_squashed(rho_1);
1207
1208 // Wall position in bottom right corner:
1209
1210 // Pointer to GeomObject
1211 GeomObject* obj_br_pt = geom_object_pt[0];
1212
1213 // Eulerian dimension
1214 unsigned n_dim = obj_br_pt->ndim();
1215
1216 // Local coordinate
1217 Vector<double> s_br(2);
1218 s_br[0] = ref_value[3];
1219 s_br[1] = ref_value[4];
1220
1221 // Get wall position
1222 Vector<double> r_br(n_dim);
1223 obj_br_pt->position(t, s_br, r_br);
1224
1225 // Wall position in top left corner:
1226
1227 // Pointer to GeomObject
1228 GeomObject* obj_tl_pt = node_pt->geom_object_pt(1);
1229
1230 // Local coordinate
1231 Vector<double> s_tl(2);
1232 s_tl[0] = ref_value[5];
1233 s_tl[1] = ref_value[6];
1234
1235 Vector<double> r_tl(n_dim);
1236 obj_tl_pt->position(t, s_tl, r_tl);
1237
1238 // Wall position at reference point:
1239
1240 // Pointer to GeomObject
1241 GeomObject* obj_wall_pt = node_pt->geom_object_pt(2);
1242
1243 // Local coordinate
1244 Vector<double> s_wall(2);
1245 s_wall[0] = ref_value[7];
1246 s_wall[1] = ref_value[8];
1247
1248 Vector<double> r_wall(n_dim);
1249 obj_wall_pt->position(t, s_wall, r_wall);
1250
1251 // Position vector to corner of the central region
1252 Vector<double> r_corner(n_dim);
1253 r_corner[0] = Lambda_x * r_br[0];
1254 r_corner[1] = Lambda_y * r_tl[1];
1255 r_corner[2] = (r_br[2] + r_tl[2]) * 0.5;
1256
1257 // Position Vector to top face of central region
1258 Vector<double> r_top(n_dim);
1259 r_top[0] = rho_0 * r_corner[0];
1260 r_top[1] = r_corner[1];
1261 r_top[2] = r_corner[2];
1262
1263 // Assign new nodal coordinate
1264 node_pt->x(t, 0) = r_top[0] + rho_1 * (r_wall[0] - r_top[0]);
1265 node_pt->x(t, 1) = r_top[1] + rho_1 * (r_wall[1] - r_top[1]);
1266 node_pt->x(t, 2) = r_top[2] + rho_1 * (r_wall[2] - r_top[2]);
1267 }
1268
1269 //======================================================================
1270 /// Update algebraic update function for nodes in region defined by
1271 /// region_id.
1272 //======================================================================
1273 template<class ELEMENT>
1274 void AlgebraicRefineableQuarterTubeMesh<
1275 ELEMENT>::update_node_update_in_region(AlgebraicNode*& node_pt,
1276 int& region_id)
1277 {
1278 // Extract references by copy construction
1279 Vector<double> ref_value(node_pt->vector_ref_value(region_id));
1280
1281 // Extract geometric objects for update by copy construction
1282 Vector<GeomObject*> geom_object_pt(
1283 node_pt->vector_geom_object_pt(region_id));
1284
1285 // Now remove the update info to allow overwriting below
1286 node_pt->kill_node_update_info(region_id);
1287
1288 // Fixed reference value: Start coordinate on wall
1289 double xi_lo = QuarterTubeMesh<ELEMENT>::Xi_lo[1];
1290
1291 // Fixed reference value: Fractional position of dividing line
1292 double fract_mid = QuarterTubeMesh<ELEMENT>::Fract_mid;
1293
1294 // Fixed reference value: End coordinate on wall
1295 double xi_hi = QuarterTubeMesh<ELEMENT>::Xi_hi[1];
1296
1297 // get initial z-coordinate of node
1298 // NOTE: use modified z found using axial_spacing_fct()
1299 // to implement axial spacing
1300 double z = ref_value[2];
1301
1302 // Update reference to wall point in bottom right corner
1303 //------------------------------------------------------
1304
1305 // Wall position in bottom right corner:
1306 Vector<double> xi(2);
1307 xi[0] = z;
1308 xi[1] = xi_lo;
1309
1310 // Find corresponding wall element/local coordinate
1311 GeomObject* obj_br_pt;
1312 Vector<double> s_br(2);
1313 this->Wall_pt->locate_zeta(xi, obj_br_pt, s_br);
1314
1315 // Wall element at bottom right end of wall mesh:
1316 geom_object_pt[0] = obj_br_pt;
1317
1318 // Local coordinate in this wall element.
1319 ref_value[3] = s_br[0];
1320 ref_value[4] = s_br[1];
1321
1322 // Update reference to wall point in upper left corner
1323 //----------------------------------------------------
1324
1325 // Wall position in top left corner
1326 xi[1] = xi_hi;
1327
1328 // Find corresponding wall element/local coordinate
1329 GeomObject* obj_tl_pt;
1330 Vector<double> s_tl(2);
1331 this->Wall_pt->locate_zeta(xi, obj_tl_pt, s_tl);
1332
1333 // Wall element at top left end of wall mesh:
1334 geom_object_pt[1] = obj_tl_pt;
1335
1336 // Local coordinate in this wall element.
1337 ref_value[5] = s_tl[0];
1338 ref_value[6] = s_tl[1];
1339
1340 if (region_id != Central_region)
1341 {
1342 // Update reference to reference point on wall
1343 //--------------------------------------------
1344
1345 // Reference point on wall
1346 if (region_id == Lower_right_region)
1347 {
1348 // Second reference value: fractional s1-position of node inside region
1349 double rho_1 = ref_value[1];
1350
1351 // position of reference point
1352 xi[1] = xi_lo + rho_1 * fract_mid * (xi_hi - xi_lo);
1353 }
1354 else // case of Upper_left region
1355 {
1356 // First reference value: fractional s0-position of node inside region
1357 double rho_0 = ref_value[0];
1358
1359 // position of reference point
1360 xi[1] = xi_hi - rho_0 * (1.0 - fract_mid) * (xi_hi - xi_lo);
1361 }
1362
1363 // Identify wall element number and local coordinate of
1364 // reference point on wall
1365 GeomObject* obj_wall_pt;
1366 Vector<double> s_wall(2);
1367 this->Wall_pt->locate_zeta(xi, obj_wall_pt, s_wall);
1368
1369 // Wall element at that contians reference point:
1370 geom_object_pt[2] = obj_wall_pt;
1371
1372 // Local coordinate in this wall element.
1373 ref_value[7] = s_wall[0];
1374 ref_value[8] = s_wall[1];
1375 }
1376
1377 // Setup algebraic update for node: Pass update information
1378 node_pt->add_node_update_info(region_id, // Enumerated ID
1379 this, // mesh
1380 geom_object_pt, // vector of geom objects
1381 ref_value); // vector of ref. vals
1382 }
1383
1384} // namespace oomph
1385#endif
oomph::AlgebraicRefineableQuarterTubeMesh
AlgebraicMesh version of RefineableQuarterTubeMesh.
Definition quarter_tube_mesh.h:433
oomph::AlgebraicRefineableQuarterTubeMesh::node_update_central_region
void node_update_central_region(const unsigned &t, AlgebraicNode *&node_pt)
Algebraic update function for a node that is located in the central region.
Definition quarter_tube_mesh.template.cc:950
oomph::QuarterTubeDomain
Quarter tube as domain. Domain is bounded by curved boundary which is represented by a GeomObject....
Definition quarter_tube_domain.h:43
oomph::QuarterTubeMesh
3D quarter tube mesh class. The domain is specified by the GeomObject that identifies boundary 3....
Definition quarter_tube_mesh.h:64
oomph::QuarterTubeMesh::wall_pt
GeomObject *& wall_pt()
Access function to GeomObject representing wall.
Definition quarter_tube_mesh.h:85
oomph::QuarterTubeMesh::Domain_pt
QuarterTubeDomain * Domain_pt
Pointer to domain.
Definition quarter_tube_mesh.h:120
oomph::QuarterTubeMesh::Xi_lo
Vector< double > Xi_lo
Lower limits for the coordinates along the wall.
Definition quarter_tube_mesh.h:126
oomph::QuarterTubeMesh::QuarterTubeMesh
QuarterTubeMesh(GeomObject *wall_pt, const Vector< double > &xi_lo, const double &fract_mid, const Vector< double > &xi_hi, const unsigned &nlayer, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor: Pass pointer to geometric object that specifies the wall, start and end coordinates on t...
Definition quarter_tube_mesh.template.cc:45
oomph::QuarterTubeMesh::Xi_hi
Vector< double > Xi_hi
Upper limits for the coordinates along the wall.
Definition quarter_tube_mesh.h:132
oomph::RefineableTriangleMesh
Unstructured refineable Triangle Mesh.
Definition triangle_mesh.h:2225