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