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