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-2023 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 
32 namespace 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::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::QuarterTubeMesh::wall_pt
GeomObject *& wall_pt()
Access function to GeomObject representing wall.
Definition: quarter_tube_mesh.template.h:86
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: annular_domain.h:35
quarter_tube_mesh.template.h