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