cylinder_with_flag_mesh.template.cc
Go to the documentation of this file.
1 // LIC// ====================================================================
2 // LIC// This file forms part of oomph-lib, the object-oriented,
3 // LIC// multi-physics finite-element library, available
4 // LIC// at http://www.oomph-lib.org.
5 // LIC//
6 // LIC// Copyright (C) 2006-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #ifndef OOMPH_CYLINDER_WITH_FLAG_MESH_TEMPLATE_CC
27 #define OOMPH_CYLINDER_WITH_FLAG_MESH_TEMPLATE_CC
28 
29 // Include the headers file
31 
32 
33 namespace oomph
34 {
35  //=============================================================
36  /// Constructor. Pass the pointers to the GeomObjects that parametrise
37  /// the cylinder, the three edges of the flag, the length and height of the
38  /// domain, the length and height of the flag, the coordinates of the
39  /// centre of the cylinder and its radius. Timestepper defaults to Steady
40  /// default timestepper.
41  //=============================================================
42  template<class ELEMENT>
44  Circle* cylinder_pt,
45  GeomObject* top_flag_pt,
46  GeomObject* bottom_flag_pt,
47  GeomObject* tip_flag_pt,
48  const double& length,
49  const double& height,
50  const double& flag_length,
51  const double& flag_height,
52  const double& centre_x,
53  const double& centre_y,
54  const double& a,
55  TimeStepper* time_stepper_pt)
56  {
57  // Mesh can only be built with 2D Qelements.
58  MeshChecker::assert_geometric_element<QElementGeometricBase, ELEMENT>(2);
59 
60  // Create the domain
61  Domain_pt = new CylinderWithFlagDomain(cylinder_pt,
62  top_flag_pt,
63  bottom_flag_pt,
64  tip_flag_pt,
65  length,
66  height,
67  flag_length,
68  flag_height,
69  centre_x,
70  centre_y,
71  a);
72 
73  // Initialise the node counter
74  unsigned long node_count = 0;
75 
76  // Vectors used to get data from domains
77  Vector<double> s(2);
78  Vector<double> r(2);
79 
80  // Setup temporary storage for the Node
81  Vector<Node*> tmp_node_pt;
82 
83  // Now blindly loop over the macro elements and associate and finite
84  // element with each
85  unsigned nmacro_element = Domain_pt->nmacro_element();
86  for (unsigned e = 0; e < nmacro_element; e++)
87  {
88  // Create the FiniteElement and add to the Element_pt Vector
89  Element_pt.push_back(new ELEMENT);
90 
91  // Read out the number of linear points in the element
92  unsigned np =
93  dynamic_cast<ELEMENT*>(this->finite_element_pt(e))->nnode_1d();
94 
95  // Loop over nodes in the column
96  for (unsigned l1 = 0; l1 < np; l1++)
97  {
98  // Loop over the nodes in the row
99  for (unsigned l2 = 0; l2 < np; l2++)
100  {
101  // Allocate the memory for the node
102  tmp_node_pt.push_back(
103  this->finite_element_pt(e)->construct_boundary_node(
104  l1 * np + l2, time_stepper_pt));
105 
106  // Read out the position of the node from the macro element
107  s[0] = -1.0 + 2.0 * (double)l2 / (double)(np - 1);
108  s[1] = -1.0 + 2.0 * (double)l1 / (double)(np - 1);
109  Domain_pt->macro_element_pt(e)->macro_map(s, r);
110 
111  // Set the position of the node
112  tmp_node_pt[node_count]->x(0) = r[0];
113  tmp_node_pt[node_count]->x(1) = r[1];
114 
115  // Increment the node number
116  node_count++;
117  }
118  }
119  } // End of loop over macro elements
120 
121  // Now the elements have been created, but there will be nodes in
122  // common, need to loop over the common edges and sort it, by reassigning
123  // pointers and the deleting excess nodes
124 
125  // Read out the number of linear points in the element
126  unsigned np =
127  dynamic_cast<ELEMENT*>(this->finite_element_pt(0))->nnode_1d();
128 
129  // Edge between Elements 0 and 1
130  for (unsigned n = 0; n < np; n++)
131  {
132  // Set the nodes in element 1 to be the same as in element 0
133  this->finite_element_pt(1)->node_pt(n * np) =
134  this->finite_element_pt(0)->node_pt((np - 1) * np + np - 1 - n);
135 
136  // Remove the nodes in element 1 from the temporaray node list
137  delete tmp_node_pt[np * np + n * np];
138  tmp_node_pt[np * np + n * np] = 0;
139  }
140 
141  // Edge between Elements 1 and 2
142  for (unsigned n = 0; n < np; n++)
143  {
144  // Set the nodes in element 2 to be the same as in element 1
145  this->finite_element_pt(2)->node_pt(n * np) =
146  this->finite_element_pt(1)->node_pt(np * n + np - 1);
147 
148  // Remove the nodes in element 2 from the temporaray node list
149  delete tmp_node_pt[2 * np * np + n * np];
150  tmp_node_pt[2 * np * np + n * np] = 0;
151  }
152 
153  // Edge between Element 2 and 3
154  for (unsigned n = 0; n < np; n++)
155  {
156  // Set the nodes in element 3 to be the same as in element 2
157  this->finite_element_pt(3)->node_pt(np * (np - 1) + n) =
158  this->finite_element_pt(2)->node_pt(np * n + np - 1);
159 
160  // Remove the nodes in element 3 from the temporaray node list
161  delete tmp_node_pt[3 * np * np + np * (np - 1) + n];
162  tmp_node_pt[3 * np * np + np * (np - 1) + n] = 0;
163  }
164 
165  // Edge between Element 5 and 4
166  for (unsigned n = 0; n < np; n++)
167  {
168  // Set the nodes in element 4 to be the same as in element 5
169  this->finite_element_pt(4)->node_pt(n) =
170  this->finite_element_pt(5)->node_pt(np * (np - n - 1) + np - 1);
171 
172  // Remove the nodes in element 4 from the temporaray node list
173  delete tmp_node_pt[4 * np * np + n];
174  tmp_node_pt[4 * np * np + n] = 0;
175  }
176 
177  // Edge between Elements 6 and 5
178  for (unsigned n = 0; n < np; n++)
179  {
180  // Set the nodes in element 5 to be the same as in element 6
181  this->finite_element_pt(5)->node_pt(n * np) =
182  this->finite_element_pt(6)->node_pt(np * n + np - 1);
183 
184  // Remove the nodes in element 5 from the temporaray node list
185  delete tmp_node_pt[5 * np * np + n * np];
186  tmp_node_pt[5 * np * np + n * np] = 0;
187  }
188 
189  // Edge between Elements 0 and 6
190  for (unsigned n = 0; n < np; n++)
191  {
192  // Set the nodes in element 6 to be the same as in element 0
193  this->finite_element_pt(6)->node_pt(n * np) =
194  this->finite_element_pt(0)->node_pt(n);
195 
196  // Remove the nodes in element 6 from the temporaray node list
197  delete tmp_node_pt[6 * np * np + n * np];
198  tmp_node_pt[6 * np * np + n * np] = 0;
199  }
200 
201  // Edge between Elements 2 and 7
202  for (unsigned n = 0; n < np; n++)
203  {
204  // Set the nodes in element 7 to be the same as in element 2
205  this->finite_element_pt(7)->node_pt(n * np) =
206  this->finite_element_pt(2)->node_pt((np - 1) * np + np - 1 - n);
207 
208  // Remove the nodes in element 7 from the temporaray node list
209  delete tmp_node_pt[7 * np * np + n * np];
210  tmp_node_pt[7 * np * np + n * np] = 0;
211  }
212 
213  // Edge between Elements 3 and 8
214  for (unsigned n = 0; n < np; n++)
215  {
216  // Set the nodes in element 8 to be the same as in element 3
217  this->finite_element_pt(8)->node_pt(n * np) =
218  this->finite_element_pt(3)->node_pt(np * n + np - 1);
219 
220  // Remove the nodes in element 8 from the temporaray node list
221  delete tmp_node_pt[8 * np * np + n * np];
222  tmp_node_pt[8 * np * np + n * np] = 0;
223  }
224 
225  // Edge between Elements 4 and 9
226  for (unsigned n = 0; n < np; n++)
227  {
228  // Set the nodes in element 9 to be the same as in element 4
229  this->finite_element_pt(9)->node_pt(n * np) =
230  this->finite_element_pt(4)->node_pt(np * n + np - 1);
231 
232  // Remove the nodes in element 9 from the temporaray node list
233  delete tmp_node_pt[9 * np * np + n * np];
234  tmp_node_pt[9 * np * np + n * np] = 0;
235  }
236 
237 
238  // Edge between Elements 5 and 10
239  for (unsigned n = 0; n < np; n++)
240  {
241  // Set the nodes in element 10 to be the same as in element 5
242  this->finite_element_pt(10)->node_pt(n * np) =
243  this->finite_element_pt(5)->node_pt(n);
244 
245  // Remove the nodes in element 10 from the temporaray node list
246  delete tmp_node_pt[10 * np * np + n * np];
247  tmp_node_pt[10 * np * np + n * np] = 0;
248  }
249 
250 
251  // Edge between Elements 7 and 11
252  for (unsigned n = 0; n < np; n++)
253  {
254  // Set the nodes in element 11 to be the same as in element 7
255  this->finite_element_pt(11)->node_pt(n * np) =
256  this->finite_element_pt(7)->node_pt(np * n + np - 1);
257 
258  // Remove the nodes in element 11 from the temporaray node list
259  delete tmp_node_pt[11 * np * np + n * np];
260  tmp_node_pt[11 * np * np + n * np] = 0;
261  }
262 
263  // Edge between Elements 8 and 12
264  for (unsigned n = 0; n < np; n++)
265  {
266  // Set the nodes in element 12 to be the same as in element 8
267  this->finite_element_pt(12)->node_pt(n * np) =
268  this->finite_element_pt(8)->node_pt(np * n + np - 1);
269 
270  // Remove the nodes in element 12 from the temporaray node list
271  delete tmp_node_pt[12 * np * np + n * np];
272  tmp_node_pt[12 * np * np + n * np] = 0;
273  }
274 
275  // Edge between Elements 9 and 13
276  for (unsigned n = 0; n < np; n++)
277  {
278  // Set the nodes in element 13 to be the same as in element 9
279  this->finite_element_pt(13)->node_pt(n * np) =
280  this->finite_element_pt(9)->node_pt(np * n + np - 1);
281 
282  // Remove the nodes in element 13 from the temporaray node list
283  delete tmp_node_pt[13 * np * np + n * np];
284  tmp_node_pt[13 * np * np + n * np] = 0;
285  }
286 
287  // Edge between Elements 10 and 14
288  for (unsigned n = 0; n < np; n++)
289  {
290  // Set the nodes in element 14 to be the same as in element 10
291  this->finite_element_pt(14)->node_pt(n * np) =
292  this->finite_element_pt(10)->node_pt(np * n + np - 1);
293 
294  // Remove the nodes in element 14 from the temporaray node list
295  delete tmp_node_pt[14 * np * np + n * np];
296  tmp_node_pt[14 * np * np + n * np] = 0;
297  }
298 
299 
300  // Edge between Elements 7 and 8
301  for (unsigned n = 0; n < np; n++)
302  {
303  // Set the nodes in element 8 to be the same as in element 7
304  this->finite_element_pt(8)->node_pt(np * (np - 1) + n) =
305  this->finite_element_pt(7)->node_pt(n);
306 
307  // Remove the nodes in element 8 from the temporaray node list
308  delete tmp_node_pt[8 * np * np + np * (np - 1) + n];
309  tmp_node_pt[8 * np * np + np * (np - 1) + n] = 0;
310  }
311 
312 
313  // Edge between Elements 9 and 10
314  for (unsigned n = 0; n < np; n++)
315  {
316  // Set the nodes in element 10 to be the same as in element 9
317  this->finite_element_pt(10)->node_pt(np * (np - 1) + n) =
318  this->finite_element_pt(9)->node_pt(n);
319 
320  // Remove the nodes in element 10 from the temporaray node list
321  delete tmp_node_pt[10 * np * np + np * (np - 1) + n];
322  tmp_node_pt[10 * np * np + np * (np - 1) + n] = 0;
323  }
324 
325 
326  // Edge between Elements 11 and 15
327  for (unsigned n = 0; n < np; n++)
328  {
329  // Set the nodes in element 15 to be the same as in element 11
330  this->finite_element_pt(15)->node_pt(n * np) =
331  this->finite_element_pt(11)->node_pt(np * n + np - 1);
332 
333  // Remove the nodes in element 15 from the temporaray node list
334  delete tmp_node_pt[15 * np * np + n * np];
335  tmp_node_pt[15 * np * np + n * np] = 0;
336  }
337 
338  // Edge between Elements 12 and 16
339  for (unsigned n = 0; n < np; n++)
340  {
341  // Set the nodes in element 16 to be the same as in element 12
342  this->finite_element_pt(16)->node_pt(n * np) =
343  this->finite_element_pt(12)->node_pt(np * n + np - 1);
344 
345  // Remove the nodes in element 16 from the temporaray node list
346  delete tmp_node_pt[16 * np * np + n * np];
347  tmp_node_pt[16 * np * np + n * np] = 0;
348  }
349 
350  // Edge between Elements 13 and 17
351  for (unsigned n = 0; n < np; n++)
352  {
353  // Set the nodes in element 17 to be the same as in element 13
354  this->finite_element_pt(17)->node_pt(n * np) =
355  this->finite_element_pt(13)->node_pt(np * n + np - 1);
356 
357  // Remove the nodes in element 17 from the temporaray node list
358  delete tmp_node_pt[17 * np * np + n * np];
359  tmp_node_pt[17 * np * np + n * np] = 0;
360  }
361 
362  // Edge between Elements 14 and 18
363  for (unsigned n = 0; n < np; n++)
364  {
365  // Set the nodes in element 18 to be the same as in element 14
366  this->finite_element_pt(18)->node_pt(n * np) =
367  this->finite_element_pt(14)->node_pt(np * n + np - 1);
368 
369  // Remove the nodes in element 18 from the temporaray node list
370  delete tmp_node_pt[18 * np * np + n * np];
371  tmp_node_pt[18 * np * np + n * np] = 0;
372  }
373 
374 
375  // Edge between Elements 11 and 12
376  for (unsigned n = 0; n < np; n++)
377  {
378  // Set the nodes in element 12 to be the same as in element 11
379  this->finite_element_pt(12)->node_pt(np * (np - 1) + n) =
380  this->finite_element_pt(11)->node_pt(n);
381 
382  // Remove the nodes in element 12 from the temporaray node list
383  delete tmp_node_pt[12 * np * np + np * (np - 1) + n];
384  tmp_node_pt[12 * np * np + np * (np - 1) + n] = 0;
385  }
386 
387 
388  // Edge between Elements 13 and 14
389  for (unsigned n = 0; n < np; n++)
390  {
391  // Set the nodes in element 14 to be the same as in element 13
392  this->finite_element_pt(14)->node_pt(np * (np - 1) + n) =
393  this->finite_element_pt(13)->node_pt(n);
394 
395  // Remove the nodes in element 14 from the temporaray node list
396  delete tmp_node_pt[14 * np * np + np * (np - 1) + n];
397  tmp_node_pt[14 * np * np + np * (np - 1) + n] = 0;
398  }
399 
400 
401  // Edge between Element 15 and 19
402  for (unsigned n = 0; n < np; n++)
403  {
404  // Set the nodes in element 19 to be the same as in element 15
405  this->finite_element_pt(19)->node_pt(np * (np - 1) + n) =
406  this->finite_element_pt(15)->node_pt(np * n + np - 1);
407 
408  // Remove the nodes in element 19 from the temporaray node list
409  delete tmp_node_pt[19 * np * np + np * (np - 1) + n];
410  tmp_node_pt[19 * np * np + np * (np - 1) + n] = 0;
411  }
412 
413 
414  // Edge between Elements 19 and 16
415  for (unsigned n = 0; n < np; n++)
416  {
417  // Set the nodes in element 19 to be the same as in element 16
418  this->finite_element_pt(16)->node_pt(np * n + np - 1) =
419  this->finite_element_pt(19)->node_pt(n * np);
420 
421  // Remove the nodes in element 16 from the temporaray node list
422  delete tmp_node_pt[16 * np * np + np * (np - 1) + n];
423  tmp_node_pt[16 * np * np + np * (np - 1) + n] = 0;
424  }
425 
426 
427  // Edge between Elements 15 and 16
428  for (unsigned n = 0; n < np; n++)
429  {
430  // Set the nodes in element 16 to be the same as in element 15
431  this->finite_element_pt(16)->node_pt(np * (np - 1) + n) =
432  this->finite_element_pt(15)->node_pt(n);
433 
434  // Remove the nodes in element 16 from the temporaray node list
435  delete tmp_node_pt[16 * np * np + np * (np - 1) + n];
436  tmp_node_pt[16 * np * np + np * (np - 1) + n] = 0;
437  }
438 
439 
440  // Edge between Element 18 and 20
441  for (unsigned n = 0; n < np; n++)
442  {
443  // Set the nodes in element 20 to be the same as in element 18
444  this->finite_element_pt(20)->node_pt(n) =
445  this->finite_element_pt(18)->node_pt(np * (np - n - 1) + np - 1);
446 
447  // Remove the nodes in element 20 from the temporaray node list
448  delete tmp_node_pt[20 * np * np + n];
449  tmp_node_pt[20 * np * np + n] = 0;
450  }
451 
452 
453  // Edge between Elements 17 and 20
454  for (unsigned n = 0; n < np; n++)
455  {
456  // Set the nodes in element 20 to be the same as in element 17
457  this->finite_element_pt(20)->node_pt(n * np) =
458  this->finite_element_pt(17)->node_pt(np * n + np - 1);
459 
460  // Remove the nodes in element 20 from the temporaray node list
461  delete tmp_node_pt[20 * np * np + n * np];
462  tmp_node_pt[20 * np * np + n * np] = 0;
463  }
464 
465 
466  // Edge between Elements 17 and 18
467  for (unsigned n = 0; n < np; n++)
468  {
469  // Set the nodes in element 18 to be the same as in element 17
470  this->finite_element_pt(18)->node_pt(np * (np - 1) + n) =
471  this->finite_element_pt(17)->node_pt(n);
472 
473  // Remove the nodes in element 18 from the temporaray node list
474  delete tmp_node_pt[18 * np * np + np * (np - 1) + n];
475  tmp_node_pt[18 * np * np + np * (np - 1) + n] = 0;
476  }
477 
478 
479  // Edge between Elements 19 and 21
480  for (unsigned n = 0; n < np; n++)
481  {
482  // Set the nodes in element 21 to be the same as in element 19
483  this->finite_element_pt(21)->node_pt(n * np) =
484  this->finite_element_pt(19)->node_pt(np * n + np - 1);
485 
486  // Remove the nodes in element 21 from the temporaray node list
487  delete tmp_node_pt[21 * np * np + n * np];
488  tmp_node_pt[21 * np * np + n * np] = 0;
489  }
490 
491 
492  // Edge between Elements 21 and 22
493  for (unsigned n = 0; n < np; n++)
494  {
495  // Set the nodes in element 22 to be the same as in element 21
496  this->finite_element_pt(22)->node_pt(np * (np - 1) + n) =
497  this->finite_element_pt(21)->node_pt(n);
498 
499  // Remove the nodes in element 22 from the temporaray node list
500  delete tmp_node_pt[22 * np * np + np * (np - 1) + n];
501  tmp_node_pt[22 * np * np + np * (np - 1) + n] = 0;
502  }
503 
504 
505  // Edge between Elements 20 and 23
506  for (unsigned n = 0; n < np; n++)
507  {
508  // Set the nodes in element 23 to be the same as in element 20
509  this->finite_element_pt(23)->node_pt(n * np) =
510  this->finite_element_pt(20)->node_pt(np * n + np - 1);
511 
512  // Remove the nodes in element 23 from the temporaray node list
513  delete tmp_node_pt[23 * np * np + n * np];
514  tmp_node_pt[23 * np * np + n * np] = 0;
515  }
516 
517 
518  // Edge between Elements 23 and 22
519  for (unsigned n = 0; n < np; n++)
520  {
521  // Set the nodes in element 22 to be the same as in element 23
522  this->finite_element_pt(22)->node_pt(n) =
523  this->finite_element_pt(23)->node_pt(np * (np - 1) + n);
524 
525  // Remove the nodes in element 22 from the temporaray node list
526  delete tmp_node_pt[22 * np * np + np * (np - 1) + n];
527  tmp_node_pt[22 * np * np + np * (np - 1) + n] = 0;
528  }
529 
530 
531  // Edge between Elements 21 and 24
532  for (unsigned n = 0; n < np; n++)
533  {
534  // Set the nodes in element 24 to be the same as in element 21
535  this->finite_element_pt(24)->node_pt(n * np) =
536  this->finite_element_pt(21)->node_pt(np * n + np - 1);
537 
538  // Remove the nodes in element 24 from the temporaray node list
539  delete tmp_node_pt[24 * np * np + n * np];
540  tmp_node_pt[24 * np * np + n * np] = 0;
541  }
542 
543 
544  // Edge between Elements 22 and 25
545  for (unsigned n = 0; n < np; n++)
546  {
547  // Set the nodes in element 25 to be the same as in element 22
548  this->finite_element_pt(25)->node_pt(n * np) =
549  this->finite_element_pt(22)->node_pt(np * n + np - 1);
550 
551  // Remove the nodes in element 25 from the temporaray node list
552  delete tmp_node_pt[25 * np * np + n * np];
553  tmp_node_pt[25 * np * np + n * np] = 0;
554  }
555 
556 
557  // Edge between Elements 23 and 26
558  for (unsigned n = 0; n < np; n++)
559  {
560  // Set the nodes in element 26 to be the same as in element 23
561  this->finite_element_pt(26)->node_pt(n * np) =
562  this->finite_element_pt(23)->node_pt(np * n + np - 1);
563 
564  // Remove the nodes in element 26 from the temporaray node list
565  delete tmp_node_pt[26 * np * np + n * np];
566  tmp_node_pt[26 * np * np + n * np] = 0;
567  }
568 
569 
570  // Edge between Elements 24 and 25
571  for (unsigned n = 0; n < np; n++)
572  {
573  // Set the nodes in element 25 to be the same as in element 24
574  this->finite_element_pt(25)->node_pt(np * (np - 1) + n) =
575  this->finite_element_pt(24)->node_pt(n);
576 
577  // Remove the nodes in element 25 from the temporaray node list
578  delete tmp_node_pt[25 * np * np + np * (np - 1) + n];
579  tmp_node_pt[25 * np * np + np * (np - 1) + n] = 0;
580  }
581 
582 
583  // Edge between Elements 25 and 26
584  for (unsigned n = 0; n < np; n++)
585  {
586  // Set the nodes in element 26 to be the same as in element 25
587  this->finite_element_pt(26)->node_pt(np * (np - 1) + n) =
588  this->finite_element_pt(25)->node_pt(n);
589 
590  // Remove the nodes in element 26 from the temporaray node list
591  delete tmp_node_pt[26 * np * np + np * (np - 1) + n];
592  tmp_node_pt[26 * np * np + np * (np - 1) + n] = 0;
593  }
594 
595 
596  // Edge between Element 24 and 27
597  for (unsigned n = 0; n < np; n++)
598  {
599  // Set the nodes in element 27 to be the same as in element 24
600  this->finite_element_pt(27)->node_pt(np * (np - 1) + n) =
601  this->finite_element_pt(24)->node_pt(np * n + np - 1);
602 
603  // Remove the nodes in element 27 from the temporaray node list
604  delete tmp_node_pt[27 * np * np + np * (np - 1) + n];
605  tmp_node_pt[27 * np * np + np * (np - 1) + n] = 0;
606  }
607 
608 
609  // Edge between Elements 25 and 27
610  for (unsigned n = 0; n < np; n++)
611  {
612  // Set the nodes in element 27 to be the same as in element 25
613  this->finite_element_pt(27)->node_pt(n * np) =
614  this->finite_element_pt(25)->node_pt(np * n + np - 1);
615 
616  // Remove the nodes in element 27 from the temporaray node list
617  delete tmp_node_pt[27 * np * np + n * np];
618  tmp_node_pt[27 * np * np + n * np] = 0;
619  }
620 
621 
622  // Edge between Element 26 and 27
623  for (unsigned n = 0; n < np; n++)
624  {
625  // Set the nodes in element 27 to be the same as in element 26
626  this->finite_element_pt(27)->node_pt(n) =
627  this->finite_element_pt(26)->node_pt(np * (np - n - 1) + np - 1);
628 
629  // Remove the nodes in element 27 from the temporaray node list
630  delete tmp_node_pt[27 * np * np + n];
631  tmp_node_pt[27 * np * np + n] = 0;
632  }
633 
634 
635  // Edge between Elements 27 and 28
636  for (unsigned n = 0; n < np; n++)
637  {
638  // Set the nodes in element 28 to be the same as in element 27
639  this->finite_element_pt(28)->node_pt(n * np) =
640  this->finite_element_pt(27)->node_pt(np * n + np - 1);
641 
642  // Remove the nodes in element 28 from the temporaray node list
643  delete tmp_node_pt[28 * np * np + n * np];
644  tmp_node_pt[28 * np * np + n * np] = 0;
645  }
646 
647 
648  // Edge between Elements 28 and 29
649  for (unsigned n = 0; n < np; n++)
650  {
651  // Set the nodes in element 29 to be the same as in element 28
652  this->finite_element_pt(29)->node_pt(n * np) =
653  this->finite_element_pt(28)->node_pt(np * n + np - 1);
654 
655  // Remove the nodes in element 29 from the temporaray node list
656  delete tmp_node_pt[29 * np * np + n * np];
657  tmp_node_pt[29 * np * np + n * np] = 0;
658  }
659 
660 
661  // Edge between Elements 29 and 30
662  for (unsigned n = 0; n < np; n++)
663  {
664  // Set the nodes in element 30 to be the same as in element 29
665  this->finite_element_pt(30)->node_pt(n * np) =
666  this->finite_element_pt(29)->node_pt(np * n + np - 1);
667 
668  // Remove the nodes in element 30 from the temporaray node list
669  delete tmp_node_pt[30 * np * np + n * np];
670  tmp_node_pt[30 * np * np + n * np] = 0;
671  }
672 
673  // Now set the actual true nodes
674  for (unsigned long n = 0; n < node_count; n++)
675  {
676  if (tmp_node_pt[n] != 0)
677  {
678  Node_pt.push_back(tmp_node_pt[n]);
679  }
680  }
681 
682  // Finally set the nodes on the boundaries
683  this->set_nboundary(8);
684 
685  for (unsigned n = 0; n < np; n++)
686  {
687  // Left hand side
688  this->add_boundary_node(3, this->finite_element_pt(0)->node_pt(n * np));
689  // Right hand side
690  this->add_boundary_node(
691  1, this->finite_element_pt(30)->node_pt(n * np + np - 1));
692 
693  // First part of lower boundary
694  this->add_boundary_node(0, this->finite_element_pt(6)->node_pt(n));
695 
696  // First part of upper boundary
697  this->add_boundary_node(
698  2, this->finite_element_pt(1)->node_pt(np * (np - 1) + n));
699 
700  // First part of hole boundary
701  this->add_boundary_node(4, this->finite_element_pt(3)->node_pt(np * n));
702 
703  // First part of lower flag
704  this->add_boundary_node(
705  5, this->finite_element_pt(4)->node_pt(np * (np - 1) + n));
706 
707  // First part of upper flag
708  this->add_boundary_node(6, this->finite_element_pt(3)->node_pt(n));
709 
710  // Right part of flag
711  this->add_boundary_node(7, this->finite_element_pt(22)->node_pt(n * np));
712  }
713 
714  for (unsigned n = 1; n < np; n++)
715  {
716  // Second part of lower boundary
717  this->add_boundary_node(0, this->finite_element_pt(10)->node_pt(n));
718 
719  // Second part of upper boundary
720  this->add_boundary_node(
721  2, this->finite_element_pt(7)->node_pt(np * (np - 1) + n));
722 
723  // Next part of lower flag
724  this->add_boundary_node(
725  5, this->finite_element_pt(9)->node_pt(np * (np - 1) + n));
726 
727  // Next part of upper flag
728  this->add_boundary_node(6, this->finite_element_pt(8)->node_pt(n));
729  }
730 
731  for (unsigned n = np - 2; n > 0; n--)
732  {
733  // Next part of hole
734  this->add_boundary_node(4, this->finite_element_pt(2)->node_pt(n));
735  }
736 
737  for (unsigned n = 1; n < np; n++)
738  {
739  // Next part of lower boundary
740  this->add_boundary_node(0, this->finite_element_pt(14)->node_pt(n));
741 
742  // Next part of upper boundary
743  this->add_boundary_node(
744  2, this->finite_element_pt(11)->node_pt(np * (np - 1) + n));
745 
746  // Next part of lower flag
747  this->add_boundary_node(
748  5, this->finite_element_pt(13)->node_pt(np * (np - 1) + n));
749 
750  // Next part of upper flag
751  this->add_boundary_node(6, this->finite_element_pt(12)->node_pt(n));
752  }
753 
754  for (unsigned n = np - 1; n > 0; n--)
755  {
756  // Next part of hole
757  this->add_boundary_node(4, this->finite_element_pt(1)->node_pt(n));
758  }
759 
760  for (unsigned n = 1; n < np; n++)
761  {
762  // Next part of lower boundary
763  this->add_boundary_node(0, this->finite_element_pt(18)->node_pt(n));
764  // Next part of upper boundary
765  this->add_boundary_node(
766  2, this->finite_element_pt(15)->node_pt(np * (np - 1) + n));
767 
768  // Next part of lower flag
769  this->add_boundary_node(
770  5, this->finite_element_pt(17)->node_pt(np * (np - 1) + n));
771 
772  // Next part of upper flag
773  this->add_boundary_node(6, this->finite_element_pt(16)->node_pt(n));
774  }
775 
776  for (unsigned n = np - 1; n > 0; n--)
777  {
778  // Next part of hole
779  this->add_boundary_node(
780  4, this->finite_element_pt(0)->node_pt(n * np + np - 1));
781  }
782 
783 
784  for (unsigned n = 1; n < np; n++)
785  {
786  // Next part of lower boundary
787  this->add_boundary_node(0, this->finite_element_pt(23)->node_pt(n));
788  // Next part of upper boundary
789  this->add_boundary_node(
790  2, this->finite_element_pt(21)->node_pt(np * (np - 1) + n));
791 
792  // Next part of hole
793  this->add_boundary_node(
794  4, this->finite_element_pt(6)->node_pt(np * (np - 1) + n));
795 
796  // Next part of lower flag
797  this->add_boundary_node(
798  5, this->finite_element_pt(20)->node_pt(np * (np - 1) + n));
799 
800  // Next part of upper flag
801  this->add_boundary_node(6, this->finite_element_pt(19)->node_pt(n));
802  }
803 
804  for (unsigned n = 0; n < np; n++)
805  {
806  // Next part of hole
807  this->add_boundary_node(
808  4, this->finite_element_pt(6)->node_pt(np * (np - 1) + n));
809  }
810 
811 
812  for (unsigned n = 1; n < np; n++)
813  {
814  // Next part of lower boundary
815  this->add_boundary_node(0, this->finite_element_pt(26)->node_pt(n));
816  // Next part of upper boundary
817  this->add_boundary_node(
818  2, this->finite_element_pt(24)->node_pt(np * (np - 1) + n));
819 
820  // Next part of hole
821  this->add_boundary_node(
822  4, this->finite_element_pt(5)->node_pt(np * (np - 1) + n));
823  }
824 
825 
826  for (unsigned n = 1; n < np; n++)
827  {
828  // Next part of lower boundary
829  this->add_boundary_node(0, this->finite_element_pt(28)->node_pt(n));
830  // Next part of upper boundary
831  this->add_boundary_node(
832  2, this->finite_element_pt(28)->node_pt(np * (np - 1) + n));
833 
834  // Next part of hole
835  this->add_boundary_node(4, this->finite_element_pt(4)->node_pt(np * n));
836  }
837 
838  for (unsigned n = 1; n < np; n++)
839  {
840  // Next part of lower boundary
841  this->add_boundary_node(0, this->finite_element_pt(29)->node_pt(n));
842  // Next part of upper boundary
843  this->add_boundary_node(
844  2, this->finite_element_pt(29)->node_pt(np * (np - 1) + n));
845  }
846 
847  for (unsigned n = 1; n < np; n++)
848  {
849  // Next part of lower boundary
850  this->add_boundary_node(0, this->finite_element_pt(30)->node_pt(n));
851  // Next part of upper boundary
852  this->add_boundary_node(
853  2, this->finite_element_pt(30)->node_pt(np * (np - 1) + n));
854  }
855 
856 
857  this->node_update();
858  setup_boundary_element_info();
859 
860  // Set boundary coordinates on the flag
861 
862  // Vector of Lagrangian coordinates used as boundary coordinate
863  Vector<double> zeta(1);
864 
865  // loop on nodes of boundary 5
866  unsigned nnode = this->nboundary_node(5);
867  for (unsigned k = 0; k < nnode; k++)
868  {
869  Node* nod_pt = this->boundary_node_pt(5, k);
870  zeta[0] = double(k) * flag_length / double(nnode - 1);
871  nod_pt->set_coordinates_on_boundary(5, zeta);
872  }
873 
874  // loop on nodes of boundary 6
875  nnode = this->nboundary_node(6);
876  for (unsigned k = 0; k < nnode; k++)
877  {
878  Node* nod_pt = this->boundary_node_pt(6, k);
879  zeta[0] = double(k) * flag_length / double(nnode - 1);
880  nod_pt->set_coordinates_on_boundary(6, zeta);
881  }
882 
883  // loop on nodes of boundary 7
884  nnode = this->nboundary_node(7);
885  for (unsigned k = 0; k < nnode; k++)
886  {
887  Node* nod_pt = this->boundary_node_pt(7, k);
888  zeta[0] = -flag_height / 2. + double(k) / double(nnode - 1) * flag_height;
889  nod_pt->set_coordinates_on_boundary(7, zeta);
890  }
891 
892  // We have parametrised boundary 5,6 and 7
893  this->Boundary_coordinate_exists[5] = true;
894  this->Boundary_coordinate_exists[6] = true;
895  this->Boundary_coordinate_exists[7] = true;
896 
897  // Loop over all elements and set macro element pointer
898  for (unsigned e = 0; e < 31; e++)
899  {
900  dynamic_cast<ELEMENT*>(this->element_pt(e))
901  ->set_macro_elem_pt(this->Domain_pt->macro_element_pt(e));
902  }
903  }
904 
905 
906  /// ///////////////////////////////////////////////////////////////////////
907  /// ///////////////////////////////////////////////////////////////////////
908  /// ///////////////////////////////////////////////////////////////////////
909 
910 
911  //=================================================================
912  /// Setup algebraic node update
913  //=================================================================
914  template<class ELEMENT>
916  {
917  // The update function requires six parameters in some cases:
918  Vector<double> ref_value(6);
919  for (unsigned i = 0; i < 5; i++)
920  {
921  ref_value[i] = 0.0;
922  }
923 
924  // Part I : macro elements 8,12,16
925  for (unsigned k = 0; k < 3; k++)
926  {
927  FiniteElement* el_pt = this->finite_element_pt(8 + k * 4);
928  unsigned nnode = el_pt->nnode();
929  for (unsigned i = 0; i < nnode; i++)
930  {
931  // Get local coordinates
932  Vector<double> coord_loc(2);
933  el_pt->local_coordinate_of_node(i, coord_loc);
934 
935  // First reference value : horizontal fraction
936  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
937 
938  // Second reference value : vertical fraction
939  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
940 
941  // Third reference value : zeta coordinate on flag
942  ref_value[2] = double(k + 1) / 5. * Flag_length +
943  ref_value[0] * 1. / 5. * Flag_length;
944 
945  // Sub-geomobject corresponding to the zeta coordinate on the flag
946  GeomObject* geom_obj_pt;
947  Vector<double> s(1);
948  Vector<double> zeta(1);
949  zeta[0] = ref_value[2];
950  Top_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
951 
952  // Create vector of geomobject for add_node_update_info()
953  Vector<GeomObject*> geom_object_pt(1);
954  geom_object_pt[0] = geom_obj_pt;
955 
956  // Fourth reference value : local coordinate in the sub geomobject
957  ref_value[3] = s[0];
958 
959  // Fifth reference value : x coordinate
960  ref_value[4] = el_pt->node_pt(i)->x(0);
961 
962  // Setup algebraic update for node: Pass update information
963  // to AlgebraicNode:
964  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
965  ->add_node_update_info(1, // ID
966  this, // mesh
967  geom_object_pt, // vector of geom objects
968  ref_value); // vector of ref. values
969  }
970  }
971 
972  // Part II : macro elements 9,13,17
973  for (unsigned k = 0; k < 3; k++)
974  {
975  FiniteElement* el_pt = this->finite_element_pt(9 + k * 4);
976  unsigned nnode = el_pt->nnode();
977  for (unsigned i = 0; i < nnode; i++)
978  {
979  // Get local coordinates
980  Vector<double> coord_loc(2);
981  el_pt->local_coordinate_of_node(i, coord_loc);
982 
983  // First reference value : horizontal fraction
984  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
985 
986  // Second reference value : vertical fraction
987  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
988 
989  // Third reference value : zeta coordinate on flag
990  ref_value[2] = double(k + 1) / 5. * Flag_length +
991  ref_value[0] * 1. / 5. * Flag_length;
992 
993  // Sub-geomobject corresponding to the zeta coordinate on the flag
994  GeomObject* geom_obj_pt;
995  Vector<double> s(1);
996  Vector<double> zeta(1);
997  zeta[0] = ref_value[2];
998  Bottom_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
999 
1000  // Create vector of geomobject for add_node_update_info()
1001  Vector<GeomObject*> geom_object_pt(1);
1002  geom_object_pt[0] = geom_obj_pt;
1003 
1004  // Fourth reference value : local coordinate in the sub geomobject
1005  ref_value[3] = s[0];
1006 
1007  // Fifth reference value : x coordinate
1008  ref_value[4] = el_pt->node_pt(i)->x(0);
1009 
1010  // Setup algebraic update for node: Pass update information
1011  // to AlgebraicNode:
1012  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1013  ->add_node_update_info(2, // ID
1014  this, // mesh
1015  geom_object_pt, // vector of geom objects
1016  ref_value); // vector of ref. values
1017  }
1018  }
1019 
1020  // Part III : macro element 22
1021  FiniteElement* el_pt = this->finite_element_pt(22);
1022  unsigned nnode = el_pt->nnode();
1023  for (unsigned i = 0; i < nnode; i++)
1024  {
1025  // Get local coordinates
1026  Vector<double> coord_loc(2);
1027  el_pt->local_coordinate_of_node(i, coord_loc);
1028 
1029  // First reference value : horizontal fraction
1030  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1031 
1032  // Second reference value : vertical fraction
1033  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1034 
1035  // Third reference value : zeta coordinate on flag
1036  ref_value[2] = coord_loc[1] * Flag_height / 2.;
1037 
1038  // Sub-geomobject corresponding to the zeta coordinate on the flag
1039  GeomObject* geom_obj_pt;
1040  Vector<double> s(1);
1041  Vector<double> zeta(1);
1042  zeta[0] = ref_value[2];
1043  Tip_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1044 
1045  // Create vector of geomobject for add_node_update_info()
1046  Vector<GeomObject*> geom_object_pt(1);
1047  geom_object_pt[0] = geom_obj_pt;
1048 
1049  // Fourth reference value : local coordinate in the sub geomobject
1050  ref_value[3] = s[0];
1051 
1052  // Setup algebraic update for node: Pass update information
1053  // to AlgebraicNode:
1054  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1055  ->add_node_update_info(3, // ID
1056  this, // mesh
1057  geom_object_pt, // vector of geom objects
1058  ref_value); // vector of ref. values
1059  }
1060 
1061  // Part IV : macro element 21
1062  el_pt = this->finite_element_pt(21);
1063  nnode = el_pt->nnode();
1064  for (unsigned i = 0; i < nnode; i++)
1065  {
1066  // Get local coordinates
1067  Vector<double> coord_loc(2);
1068  el_pt->local_coordinate_of_node(i, coord_loc);
1069 
1070  // First reference value : horizontal fraction
1071  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1072 
1073  // Second reference value : vertical fraction
1074  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1075 
1076  // Sub-geomobject corresponding to the tip of the Tip_flag
1077  GeomObject* geom_obj_pt;
1078  Vector<double> s(1);
1079  Vector<double> zeta(1);
1080  zeta[0] = Flag_height / 2.;
1081  Tip_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1082 
1083  // Create vector of geomobject for add_node_update_info()
1084  Vector<GeomObject*> geom_object_pt(1);
1085  geom_object_pt[0] = geom_obj_pt;
1086 
1087  // Third reference value : local coordinate in the sub geomobject
1088  ref_value[2] = s[0];
1089 
1090  // Setup algebraic update for node: Pass update information
1091  // to AlgebraicNode:
1092  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1093  ->add_node_update_info(4, // ID
1094  this, // mesh
1095  geom_object_pt, // vector of geom objects
1096  ref_value); // vector of ref. values
1097  }
1098 
1099  // Part V : macro element 23
1100  el_pt = this->finite_element_pt(23);
1101  nnode = el_pt->nnode();
1102  for (unsigned i = 0; i < nnode; i++)
1103  {
1104  // Get local coordinates
1105  Vector<double> coord_loc(2);
1106  el_pt->local_coordinate_of_node(i, coord_loc);
1107 
1108  // First reference value : horizontal fraction
1109  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1110 
1111  // Second reference value : vertical fraction
1112  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1113 
1114  // Sub-geomobject corresponding to the tip of the Tip_flag
1115  GeomObject* geom_obj_pt;
1116  Vector<double> s(1);
1117  Vector<double> zeta(1);
1118  zeta[0] = -Flag_height / 2.;
1119  Tip_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1120 
1121  // Create vector of geomobject for add_node_update_info()
1122  Vector<GeomObject*> geom_object_pt(1);
1123  geom_object_pt[0] = geom_obj_pt;
1124 
1125  // Third reference value : local coordinate in the sub geomobject
1126  ref_value[2] = s[0];
1127 
1128  // Setup algebraic update for node: Pass update information
1129  // to AlgebraicNode:
1130  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1131  ->add_node_update_info(5, // ID
1132  this, // mesh
1133  geom_object_pt, // vector of geom objects
1134  ref_value); // vector of ref. values
1135  }
1136 
1137  // Part VI = macro element 19
1138  el_pt = this->finite_element_pt(19);
1139  nnode = el_pt->nnode();
1140  for (unsigned i = 0; i < nnode; i++)
1141  {
1142  // Get local coordinates
1143  Vector<double> coord_loc(2);
1144  el_pt->local_coordinate_of_node(i, coord_loc);
1145 
1146  // First reference value : horizontal fraction
1147  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1148 
1149  // Second reference value : vertical fraction
1150  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1151 
1152  // Third reference value : zeta coordinate on the flag
1153  ref_value[2] =
1154  4. / 5. * Flag_length + ref_value[0] * 1. / 5. * Flag_length;
1155 
1156  // Sub-geomobject
1157  GeomObject* geom_obj_pt;
1158  Vector<double> s(1);
1159  Vector<double> zeta(1);
1160  zeta[0] = ref_value[2];
1161  Top_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1162 
1163  // Create vector of geomobject for add_node_update_info()
1164  Vector<GeomObject*> geom_object_pt(1);
1165  geom_object_pt[0] = geom_obj_pt;
1166 
1167  // Third reference value : local coordinate in the sub geomobject
1168  ref_value[3] = s[0];
1169 
1170  // Setup algebraic update for node: Pass update information
1171  // to AlgebraicNode:
1172  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1173  ->add_node_update_info(6, // ID
1174  this, // mesh
1175  geom_object_pt, // vector of geom objects
1176  ref_value); // vector of ref. values
1177  }
1178 
1179 
1180  // Part VII = macro element 20
1181  el_pt = this->finite_element_pt(20);
1182  nnode = el_pt->nnode();
1183  for (unsigned i = 0; i < nnode; i++)
1184  {
1185  // Get local coordinates
1186  Vector<double> coord_loc(2);
1187  el_pt->local_coordinate_of_node(i, coord_loc);
1188 
1189  // First reference value : horizontal fraction
1190  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1191 
1192  // Second reference value : vertical fraction
1193  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1194 
1195  // Third reference value : zeta coordinate on the flag
1196  ref_value[2] =
1197  4. / 5. * Flag_length + ref_value[0] * 1. / 5. * Flag_length;
1198 
1199  // Sub-geomobject
1200  GeomObject* geom_obj_pt;
1201  Vector<double> s(1);
1202  Vector<double> zeta(1);
1203  zeta[0] = ref_value[2];
1204  Bottom_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1205 
1206  // Create vector of geomobject for add_node_update_info()
1207  Vector<GeomObject*> geom_object_pt(1);
1208  geom_object_pt[0] = geom_obj_pt;
1209 
1210  // Third reference value : local coordinate in the sub geomobject
1211  ref_value[3] = s[0];
1212 
1213  // Setup algebraic update for node: Pass update information
1214  // to AlgebraicNode:
1215  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1216  ->add_node_update_info(7, // ID
1217  this, // mesh
1218  geom_object_pt, // vector of geom objects
1219  ref_value); // vector of ref. values
1220  }
1221 
1222  // Part VIII : macro element 3
1223  el_pt = this->finite_element_pt(3);
1224  nnode = el_pt->nnode();
1225  for (unsigned i = 0; i < nnode; i++)
1226  {
1227  // Get local coordinates
1228  Vector<double> coord_loc(2);
1229  el_pt->local_coordinate_of_node(i, coord_loc);
1230 
1231  // First reference value : horizontal fraction
1232  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1233 
1234  // Second reference value : vertical fraction
1235  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1236 
1237  // Third reference value : zeta coordinate on flag at reference point
1238  ref_value[2] = ref_value[0] * 1. / 5. * Flag_length;
1239 
1240  // Sub-geomobject corresponding to the zeta coordinate on the flag
1241  // at the reference point
1242  GeomObject* geom_obj_pt;
1243  Vector<double> s(1);
1244  Vector<double> zeta(1);
1245  zeta[0] = ref_value[2];
1246  Top_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1247 
1248  // Fourth reference value : local coordinate in the sub geomobject
1249  ref_value[3] = s[0];
1250 
1251  // Create vector of geomobject for add_node_update_info()
1252  Vector<GeomObject*> geom_object_pt(2);
1253  geom_object_pt[0] = geom_obj_pt;
1254 
1255  // Fifth reference value : zeta coordinate on flag at end of macro element
1256  ref_value[4] = 1. / 5. * Flag_length;
1257 
1258  // Sub-geomobject corresponding to the zeta coordinate on the flag
1259  // at the end of the macro element
1260  zeta[0] = ref_value[4];
1261  Top_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1262 
1263  // Add geom object
1264  geom_object_pt[1] = geom_obj_pt;
1265 
1266  // Sixth reference value : local coordinate in the sub geomobject
1267  ref_value[5] = s[0];
1268 
1269 
1270  // Setup algebraic update for node: Pass update information
1271  // to AlgebraicNode:
1272  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1273  ->add_node_update_info(8, // ID
1274  this, // mesh
1275  geom_object_pt, // vector of geom objects
1276  ref_value); // vector of ref. values
1277  }
1278 
1279 
1280  // Part IX : macro element 4
1281  el_pt = this->finite_element_pt(4);
1282  nnode = el_pt->nnode();
1283  for (unsigned i = 0; i < nnode; i++)
1284  {
1285  // Get local coordinates
1286  Vector<double> coord_loc(2); /**set the size ??*/
1287  el_pt->local_coordinate_of_node(i, coord_loc);
1288 
1289  // First reference value : horizontal fraction
1290  ref_value[0] = 0.5 * (coord_loc[0] + 1.0);
1291 
1292  // Second reference value : vertical fraction
1293  ref_value[1] = 0.5 * (coord_loc[1] + 1.0);
1294 
1295  // Third reference value : zeta coordinate on flag
1296  ref_value[2] = ref_value[0] * 1. / 5. * Flag_length;
1297 
1298  // Sub-geomobject corresponding to the zeta coordinate on the flag
1299  // at the reference point
1300  GeomObject* geom_obj_pt;
1301  Vector<double> s(1);
1302  Vector<double> zeta(1);
1303  zeta[0] = ref_value[2];
1304  Bottom_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1305 
1306  // Fourth reference value : local coordinate in the sub geomobject
1307  ref_value[3] = s[0];
1308 
1309  // Create vector of geomobject for add_node_update_info()
1310  Vector<GeomObject*> geom_object_pt(2);
1311  geom_object_pt[0] = geom_obj_pt;
1312 
1313  // Fifth reference value : zeta coordinate on flag at end of macro element
1314  ref_value[4] = 1. / 5. * Flag_length;
1315 
1316  // Sub-geomobject corresponding to the zeta coordinate on the flag
1317  // at the end of the macro element
1318  zeta[0] = ref_value[4];
1319  Bottom_flag_pt->locate_zeta(zeta, geom_obj_pt, s);
1320 
1321  // Add geom object
1322  geom_object_pt[1] = geom_obj_pt;
1323 
1324  // Sixth reference value : local coordinate in the sub geomobject
1325  ref_value[5] = s[0];
1326 
1327  // Setup algebraic update for node: Pass update information
1328  // to AlgebraicNode:
1329  dynamic_cast<AlgebraicNode*>(el_pt->node_pt(i))
1330  ->add_node_update_info(9, // ID
1331  this, // mesh
1332  geom_object_pt, // vector of geom objects
1333  ref_value); // vector of ref. values
1334  }
1335 
1336 
1337  } // end of setup_algebraic_node_update
1338 
1339 
1340  //=================================================================
1341  /// The algebraic node update function
1342  //=================================================================
1343  template<class ELEMENT>
1345  const unsigned& t, AlgebraicNode*& node_pt)
1346  {
1347  unsigned id = node_pt->node_update_fct_id();
1348 
1349  switch (id)
1350  {
1351  case 1:
1352  node_update_I(t, node_pt);
1353  break;
1354 
1355  case 2:
1356  node_update_II(t, node_pt);
1357  break;
1358 
1359  case 3:
1360  node_update_III(t, node_pt);
1361  break;
1362 
1363  case 4:
1364  node_update_IV(t, node_pt);
1365  break;
1366 
1367  case 5:
1368  node_update_V(t, node_pt);
1369  break;
1370 
1371  case 6:
1372  node_update_VI(t, node_pt);
1373  break;
1374 
1375  case 7:
1376  node_update_VII(t, node_pt);
1377  break;
1378 
1379  case 8:
1380  node_update_VIII(t, node_pt);
1381  break;
1382 
1383  case 9:
1384  node_update_IX(t, node_pt);
1385  break;
1386 
1387  default:
1388  std::ostringstream error_message;
1389  error_message << "Wrong id " << id << std::endl;
1390  throw OomphLibError(error_message.str(),
1391  OOMPH_CURRENT_FUNCTION,
1392  OOMPH_EXCEPTION_LOCATION);
1393  }
1394 
1395 
1396  } // end of algebraic_node_update()
1397 
1398  //=================================================================
1399  /// Node update for region I
1400  //=================================================================
1401  template<class ELEMENT>
1403  const unsigned& t, AlgebraicNode*& node_pt)
1404  {
1405  // Extract reference values for update by copy construction
1406  Vector<double> ref_value(node_pt->vector_ref_value());
1407 
1408  // Extract geometric objects for update by copy construction
1409  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1410 
1411  // Pointer to geom object
1412  GeomObject* flag_pt = geom_object_pt[0];
1413 
1414  // Point on the line y=p11[1] corresponding to the initial x.
1415  Vector<double> ref_point(2);
1416  ref_point[0] = ref_value[4];
1417  ref_point[1] = 0.778024390 * Height;
1418 
1419  // Point on the flag
1420  Vector<double> flag_point(2);
1421  Vector<double> zeta(1);
1422  zeta[0] = ref_value[3];
1423  flag_pt->position(t, zeta, flag_point);
1424 
1425  // Third reference value : fraction of the vertical line
1426  // between the straight line y = p11[1] and the flag
1427  double r = ref_value[1];
1428 
1429  // Assign new nodal coordinates
1430  node_pt->x(t, 0) =
1431  ref_point[0] + (1.0 - r) * (flag_point[0] - ref_point[0]);
1432  node_pt->x(t, 1) =
1433  ref_point[1] + (1.0 - r) * (flag_point[1] - ref_point[1]);
1434  }
1435 
1436 
1437  //=================================================================
1438  /// Node update for region II
1439  //=================================================================
1440  template<class ELEMENT>
1442  const unsigned& t, AlgebraicNode*& node_pt)
1443  {
1444  // Extract reference values for update by copy construction
1445  Vector<double> ref_value(node_pt->vector_ref_value());
1446 
1447  // Extract geometric objects for update by copy construction
1448  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1449 
1450  // Pointer to geom object
1451  GeomObject* flag_pt = geom_object_pt[0];
1452 
1453  // Point on the line y=p37[1] corresponding to the initial x.
1454  Vector<double> ref_point(2);
1455  ref_point[0] = ref_value[4];
1456  ref_point[1] = 0.197585366 * Height;
1457 
1458  // Point on the flag
1459  Vector<double> flag_point(2);
1460  Vector<double> zeta(1);
1461  zeta[0] = ref_value[3];
1462  flag_pt->position(t, zeta, flag_point);
1463 
1464  // Third reference value : fraction of the vertical line
1465  // between the straight line y = p11[1] and the flag
1466  double r = ref_value[1];
1467 
1468  // Assign new nodal coordinates
1469  node_pt->x(t, 0) = ref_point[0] + r * (flag_point[0] - ref_point[0]);
1470  node_pt->x(t, 1) = ref_point[1] + r * (flag_point[1] - ref_point[1]);
1471  }
1472 
1473  //=================================================================
1474  /// Node update for region III
1475  //=================================================================
1476  template<class ELEMENT>
1478  const unsigned& t, AlgebraicNode*& node_pt)
1479  {
1480  // useful points
1481  Vector<double> p15(2);
1482  Vector<double> p35(2);
1483 
1484  p15[0] = 0.285123967 * Length;
1485  p15[1] = 0.625 * Height;
1486 
1487  p35[0] = 0.285123967 * Length;
1488  p35[1] = 0.350609756 * Height;
1489 
1490  // Extract reference values for update by copy construction
1491  Vector<double> ref_value(node_pt->vector_ref_value());
1492 
1493  // Extract geometric objects for update by copy construction
1494  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1495 
1496  // Pointer to geom object
1497  GeomObject* flag_pt = geom_object_pt[0];
1498 
1499  // Point on the line x=p15[0]
1500  Vector<double> ref_point(2);
1501  ref_point[0] = p15[0];
1502  ref_point[1] = p35[1] + ref_value[1] * (p15[1] - p35[1]);
1503 
1504  // Point on the flag
1505  Vector<double> flag_point(2);
1506  Vector<double> zeta(1);
1507  zeta[0] = ref_value[3];
1508  flag_pt->position(t, zeta, flag_point);
1509 
1510  // Third reference value : fraction of the horizontal line
1511  // between the flag and the horizontal straight line in x=p15[0]
1512  double r = ref_value[0];
1513 
1514  // Assign new nodal coordinates
1515  node_pt->x(t, 0) = flag_point[0] + r * (ref_point[0] - flag_point[0]);
1516  node_pt->x(t, 1) = flag_point[1] + r * (ref_point[1] - flag_point[1]);
1517  }
1518 
1519  //=================================================================
1520  /// Node update for region IV
1521  //=================================================================
1522  template<class ELEMENT>
1524  const unsigned& t, AlgebraicNode*& node_pt)
1525  {
1526  // Useful points
1527  Vector<double> p15(2);
1528  Vector<double> p25(2);
1529  Vector<double> top_flag(2);
1530 
1531  p15[0] = 0.285123967 * Length;
1532  p15[1] = 0.625 * Height;
1533 
1534  p25[0] = Centre_x +
1535  A * sqrt(1.0 - Flag_height * Flag_height / (4.0 * A * A)) +
1536  Flag_length;
1537  p25[1] = Centre_y + Flag_height / 2.0;
1538 
1539  // Extract reference values for update by copy construction
1540  Vector<double> ref_value(node_pt->vector_ref_value());
1541 
1542  // Extract geometric objects for update by copy construction
1543  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1544 
1545  // Pointer to geom object
1546  GeomObject* flag_pt = geom_object_pt[0];
1547 
1548  Vector<double> zeta(1);
1549  zeta[0] = ref_value[2];
1550  flag_pt->position(t, zeta, top_flag);
1551 
1552  // point on the line linking p15 et top_flag
1553  Vector<double> p1(2);
1554  p1[0] = top_flag[0] + ref_value[0] * (p15[0] - top_flag[0]);
1555  p1[1] = top_flag[1] + ref_value[0] * (p15[1] - top_flag[1]);
1556 
1557  // Point on the line y = Height;
1558  Vector<double> p2(2);
1559  p2[0] = p25[0] + ref_value[0] * (p15[0] - p25[0]);
1560  p2[1] = Height;
1561 
1562  // Connect those points with the vertical fraction ref_value[1]
1563  node_pt->x(t, 0) = p1[0] + ref_value[1] * (p2[0] - p1[0]);
1564  node_pt->x(t, 1) = p1[1] + ref_value[1] * (p2[1] - p1[1]);
1565  }
1566 
1567  //=================================================================
1568  /// Node update for region V
1569  //=================================================================
1570  template<class ELEMENT>
1572  const unsigned& t, AlgebraicNode*& node_pt)
1573  {
1574  // Useful points
1575  Vector<double> p31(2);
1576  Vector<double> p35(2);
1577  Vector<double> top_flag(2);
1578 
1579  p31[0] = Centre_x +
1580  A * sqrt(1.0 - Flag_height * Flag_height / (4.0 * A * A)) +
1581  Flag_length;
1582  p31[1] = Centre_y - Flag_height / 2.;
1583 
1584  p35[0] = 0.285123967 * Length;
1585  p35[1] = 0.350609756 * Height;
1586 
1587  // Extract reference values for update by copy construction
1588  Vector<double> ref_value(node_pt->vector_ref_value());
1589 
1590  // Extract geometric objects for update by copy construction
1591  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1592 
1593  // Pointer to geom object
1594  GeomObject* flag_pt = geom_object_pt[0];
1595 
1596  Vector<double> zeta(1);
1597  zeta[0] = ref_value[2];
1598 
1599  flag_pt->position(t, zeta, top_flag);
1600 
1601  // point on the line linking p35 et top_flag
1602  Vector<double> p1(2);
1603  p1[0] = top_flag[0] + ref_value[0] * (p35[0] - top_flag[0]);
1604  p1[1] = top_flag[1] + ref_value[0] * (p35[1] - top_flag[1]);
1605 
1606  // Point on the line y = 0.0;
1607  Vector<double> p2(2);
1608  p2[0] = p31[0] + ref_value[0] * (p35[0] - p31[0]);
1609  p2[1] = 0.;
1610 
1611  // Connect those points with the vertical fraction ref_value[1]
1612  node_pt->x(t, 0) = p2[0] + ref_value[1] * (p1[0] - p2[0]);
1613  node_pt->x(t, 1) = p2[1] + ref_value[1] * (p1[1] - p2[1]);
1614  }
1615 
1616  //=================================================================
1617  /// Node update for region VI
1618  //=================================================================
1619  template<class ELEMENT>
1621  const unsigned& t, AlgebraicNode*& node_pt)
1622  {
1623  // Extract reference values for update by copy construction
1624  Vector<double> ref_value(node_pt->vector_ref_value());
1625 
1626  // Useful points
1627  Vector<double> p14(2);
1628  Vector<double> p5(2);
1629  Vector<double> point_flag(2);
1630 
1631  p14[0] = 0.211596 * Length;
1632  p14[1] = 0.778024390 * Height;
1633 
1634  p5[0] = 0.239596 * Length;
1635  p5[1] = Height;
1636 
1637  // Extract geometric objects for update by copy construction
1638  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1639 
1640  // Pointer to geom object
1641  GeomObject* flag_pt = geom_object_pt[0];
1642 
1643  Vector<double> zeta(1);
1644  zeta[0] = ref_value[3];
1645  flag_pt->position(t, zeta, point_flag);
1646 
1647  // point on the line linking p14 et p5
1648  Vector<double> p1(2);
1649  p1[0] = p14[0] + ref_value[0] * (p5[0] - p14[0]);
1650  p1[1] = p14[1] + ref_value[0] * (p5[1] - p14[1]);
1651 
1652 
1653  // Connect those points with the vertical fraction ref_value[1]
1654  node_pt->x(t, 0) = point_flag[0] + ref_value[1] * (p1[0] - point_flag[0]);
1655  node_pt->x(t, 1) = point_flag[1] + ref_value[1] * (p1[1] - point_flag[1]);
1656  }
1657 
1658  //=================================================================
1659  /// Node update for region VII
1660  //=================================================================
1661  template<class ELEMENT>
1663  const unsigned& t, AlgebraicNode*& node_pt)
1664  {
1665  // Extract reference values for update by copy construction
1666  Vector<double> ref_value(node_pt->vector_ref_value());
1667 
1668  // Useful points
1669  Vector<double> p40(2);
1670  Vector<double> p45(2);
1671  Vector<double> point_flag(2);
1672 
1673  p40[0] = 0.211596 * Length;
1674  p40[1] = 0.197585366 * Height;
1675 
1676  p45[0] = 0.239596 * Length;
1677  p45[1] = 0.0;
1678 
1679  // Extract geometric objects for update by copy construction
1680  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1681 
1682  // Pointer to geom object
1683  GeomObject* flag_pt = geom_object_pt[0];
1684 
1685  Vector<double> zeta(1);
1686  zeta[0] = ref_value[3];
1687  flag_pt->position(t, zeta, point_flag);
1688 
1689  // point on the line linking p40 et p45
1690  Vector<double> p1(2);
1691  p1[0] = p40[0] + ref_value[0] * (p45[0] - p40[0]);
1692  p1[1] = p40[1] + ref_value[0] * (p45[1] - p40[1]);
1693 
1694 
1695  // Connect those points with the vertical fraction ref_value[1]
1696  node_pt->x(t, 0) =
1697  point_flag[0] + (1 - ref_value[1]) * (p1[0] - point_flag[0]);
1698  node_pt->x(t, 1) =
1699  point_flag[1] + (1 - ref_value[1]) * (p1[1] - point_flag[1]);
1700  }
1701 
1702 
1703  //=================================================================
1704  /// Node update for region VIII
1705  //=================================================================
1706  template<class ELEMENT>
1708  const unsigned& t, AlgebraicNode*& node_pt)
1709  {
1710  // Useful point
1711  Vector<double> p11(2);
1712  p11[0] = 0.127596 * Length;
1713  p11[1] = 0.778024390 * Height;
1714 
1715  /// Extreme angles on circle
1716  double zeta_circle_top = atan(1.0);
1717  double zeta_circle_bot = asin(Flag_height / 2. / A);
1718 
1719  // Extract reference values for update by copy construction
1720  Vector<double> ref_value(node_pt->vector_ref_value());
1721 
1722  // Extract geometric objects for update by copy construction
1723  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1724 
1725  // Pointer to geom object containing the reference point
1726  GeomObject* flag_ref_pt = geom_object_pt[0];
1727 
1728  // Pointer to geom object containing the end of the macro element
1729  GeomObject* flag_end_pt = geom_object_pt[1];
1730 
1731  double omega_horizontal = ref_value[0];
1732  double omega_vertical = ref_value[1];
1733 
1734  // end of the macro element on the flag
1735  Vector<double> flag_end(2);
1736  Vector<double> zeta(1);
1737  zeta[0] = ref_value[5];
1738  flag_end_pt->position(t, zeta, flag_end);
1739 
1740  Vector<double> outer_point(p11);
1741 
1742  // Get reference point on circle
1743  Vector<double> ref_point_on_circle(2);
1744  zeta[0] =
1745  zeta_circle_bot + (zeta_circle_top - zeta_circle_bot) * omega_vertical;
1746  Cylinder_pt->position(zeta, ref_point_on_circle);
1747 
1748  // Get reference point on right line
1749  Vector<double> ref_point_on_right_line(2);
1750  ref_point_on_right_line[0] = // flag_end[0];
1751  outer_point[0] + (flag_end[0] - outer_point[0]) * (1.0 - omega_vertical);
1752  ref_point_on_right_line[1] =
1753  outer_point[1] + (flag_end[1] - outer_point[1]) * (1.0 - omega_vertical);
1754 
1755  // Get reference point on flag
1756  Vector<double> ref_point_on_flag(2);
1757  zeta[0] = ref_value[3];
1758  flag_ref_pt->position(t, zeta, ref_point_on_flag);
1759 
1760  // Get bottom-most point on circle
1761  Vector<double> circle_bot(2);
1762  zeta[0] = zeta_circle_bot;
1763  Cylinder_pt->position(zeta, circle_bot);
1764 
1765  // Get reference point on horizontal fraction of straight line
1766  // connecting the two bottom most reference points
1767  Vector<double> r_bot(2);
1768  r_bot[0] = circle_bot[0] + (flag_end[0] - circle_bot[0]) * omega_horizontal;
1769  r_bot[1] = circle_bot[1] + (flag_end[1] - circle_bot[1]) * omega_horizontal;
1770 
1771  // Place point on horizontal fraction of straight line
1772  // connecting reference points -- this won't match the
1773  // curved top boundary adjacent to the flag
1774  node_pt->x(t, 0) =
1775  ref_point_on_circle[0] +
1776  (ref_point_on_right_line[0] - ref_point_on_circle[0]) * omega_horizontal;
1777  node_pt->x(t, 1) =
1778  ref_point_on_circle[1] +
1779  (ref_point_on_right_line[1] - ref_point_on_circle[1]) * omega_horizontal;
1780 
1781  // Correct by scaled difference between bottom straight line
1782  // and bent flag
1783  node_pt->x(t, 0) +=
1784  (ref_point_on_flag[0] - r_bot[0]) * (1.0 - omega_vertical);
1785  node_pt->x(t, 1) +=
1786  (ref_point_on_flag[1] - r_bot[1]) * (1.0 - omega_vertical);
1787  }
1788 
1789  //=================================================================
1790  /// Node update for region IX
1791  //=================================================================
1792  template<class ELEMENT>
1794  const unsigned& t, AlgebraicNode*& node_pt)
1795  {
1796  // Useful point
1797  Vector<double> p37(2);
1798  p37[0] = 0.127596 * Length;
1799  p37[1] = 0.197585366 * Height;
1800 
1801  /// Extreme angles on circle
1802  double zeta_circle_top = -asin(Flag_height / 2. / A);
1803  double zeta_circle_bot = -atan(1.0);
1804 
1805  // Extract reference values for update by copy construction
1806  Vector<double> ref_value(node_pt->vector_ref_value());
1807 
1808  // Extract geometric objects for update by copy construction
1809  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1810 
1811  // Pointer to geom object containing the reference point
1812  GeomObject* flag_ref_pt = geom_object_pt[0];
1813 
1814  // Pointer to geom object containing the end of macro element
1815  GeomObject* flag_end_pt = geom_object_pt[1];
1816 
1817  double omega_horizontal = ref_value[0];
1818  double omega_vertical = ref_value[1];
1819 
1820  // end of the macro element on the flag
1821  Vector<double> flag_end(2);
1822  Vector<double> zeta(1);
1823  zeta[0] = ref_value[5];
1824  flag_end_pt->position(t, zeta, flag_end);
1825 
1826  Vector<double> outer_point(p37);
1827 
1828  // Get reference point on circle
1829  Vector<double> ref_point_on_circle(2);
1830  zeta[0] =
1831  zeta_circle_bot + (zeta_circle_top - zeta_circle_bot) * omega_vertical;
1832  Cylinder_pt->position(zeta, ref_point_on_circle);
1833 
1834  // Get reference point on right line
1835  Vector<double> ref_point_on_right_line(2);
1836  ref_point_on_right_line[0] = // flag_end[0];
1837  outer_point[0] + (flag_end[0] - outer_point[0]) * omega_vertical;
1838  ref_point_on_right_line[1] =
1839  outer_point[1] + (flag_end[1] - outer_point[1]) * omega_vertical;
1840 
1841  // Get reference point on flag
1842  Vector<double> ref_point_on_flag(2);
1843  zeta[0] = ref_value[3];
1844  flag_ref_pt->position(t, zeta, ref_point_on_flag);
1845 
1846  // Get top-most point on circle
1847  Vector<double> circle_top(2);
1848  zeta[0] = zeta_circle_top;
1849  Cylinder_pt->position(zeta, circle_top);
1850 
1851  // Get reference point on horizontal fraction of straight line
1852  // connecting the two top most reference points
1853  Vector<double> r_top(2);
1854  r_top[0] = circle_top[0] + (flag_end[0] - circle_top[0]) * omega_horizontal;
1855  r_top[1] = circle_top[1] + (flag_end[1] - circle_top[1]) * omega_horizontal;
1856 
1857  // Place point on horizontal fraction of straight line
1858  // connecting reference points -- this won't match the
1859  // curved top boundary adjacent to the flag
1860  node_pt->x(t, 0) =
1861  ref_point_on_circle[0] +
1862  (ref_point_on_right_line[0] - ref_point_on_circle[0]) * omega_horizontal;
1863  node_pt->x(t, 1) =
1864  ref_point_on_circle[1] +
1865  (ref_point_on_right_line[1] - ref_point_on_circle[1]) * omega_horizontal;
1866 
1867  // Correct by scaled difference between top straight line
1868  // and bent flag
1869  node_pt->x(t, 0) += (ref_point_on_flag[0] - r_top[0]) * omega_vertical;
1870  node_pt->x(t, 1) += (ref_point_on_flag[1] - r_top[1]) * omega_vertical;
1871  }
1872 
1873 
1874  /// ////////////////////////////////////////////////////////////////////////
1875  /// ////////////////////////////////////////////////////////////////////////
1876  /// ////////////////////////////////////////////////////////////////////////
1877 
1878 
1879  //===================================================================
1880  /// Update the node update functions
1881  //===================================================================
1882  template<class ELEMENT>
1884  AlgebraicNode*& node_pt)
1885  {
1886  // Extract ID
1887  unsigned id = node_pt->node_update_fct_id();
1888 
1889 
1890  if (id == 8)
1891  {
1892  // Extract geometric objects for update by copy construction
1893  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1894 
1895  // Extract reference values for node update by copy construction
1896  Vector<double> ref_value(node_pt->vector_ref_value());
1897 
1898  // Get zeta coordinate of reference point
1899  Vector<double> zeta_ref_flag(1);
1900  zeta_ref_flag[0] = ref_value[2];
1901 
1902  // Get the sub-geomobject and the local coordinate containing the
1903  // reference point
1904  Vector<double> s(1);
1905  GeomObject* geom_obj_pt;
1906  this->Top_flag_pt->locate_zeta(zeta_ref_flag, geom_obj_pt, s);
1907 
1908 
1909  // Update the pointer to the (sub-)GeomObject within which the
1910  // reference point is located.
1911  geom_object_pt[0] = geom_obj_pt;
1912 
1913  // Update second reference value: Reference local coordinate
1914  // in flag sub-element
1915  ref_value[3] = s[0];
1916 
1917 
1918  // Get zeta coordinate of point at end of macro element
1919  Vector<double> zeta_end_flag(1);
1920  zeta_end_flag[0] = ref_value[4];
1921 
1922  // Get the sub-geomobject and the local coordinate containing the
1923  // point at the end of the macro element
1924  this->Top_flag_pt->locate_zeta(zeta_end_flag, geom_obj_pt, s);
1925 
1926 
1927  // Update the pointer to the (sub-)GeomObject within which the
1928  // point at the end of the macro element is located
1929  geom_object_pt[1] = geom_obj_pt;
1930 
1931  // Update second reference value: Reference local coordinate
1932  // in flag sub-element
1933  ref_value[5] = s[0];
1934 
1935  // Kill the existing node update info
1936  node_pt->kill_node_update_info(8);
1937 
1938  // Setup algebraic update for node: Pass update information
1939  node_pt->add_node_update_info(8, // id
1940  this, // mesh
1941  geom_object_pt, // vector of geom objects
1942  ref_value); // vector of ref. values
1943  }
1944  else if (id == 9)
1945  {
1946  // Extract geometric objects for update by copy construction
1947  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
1948 
1949  // Extract reference values for node update by copy construction
1950  Vector<double> ref_value(node_pt->vector_ref_value());
1951 
1952  // Get zeta coordinate of reference point
1953  Vector<double> zeta_ref_flag(1);
1954  zeta_ref_flag[0] = ref_value[2];
1955 
1956  // Get the sub-geomobject and the local coordinate containing the
1957  // reference point
1958  Vector<double> s(1);
1959  GeomObject* geom_obj_pt;
1960  this->Bottom_flag_pt->locate_zeta(zeta_ref_flag, geom_obj_pt, s);
1961 
1962  // Update the pointer to the (sub-)GeomObject within which the
1963  // reference point is located.
1964  geom_object_pt[0] = geom_obj_pt;
1965 
1966  // Update second reference value: Reference local coordinate
1967  // in flag sub-element
1968  ref_value[3] = s[0];
1969 
1970 
1971  // Get zeta coordinate of point at end of macro element
1972  Vector<double> zeta_end_flag(1);
1973  zeta_end_flag[0] = ref_value[4];
1974 
1975  // Get the sub-geomobject and the local coordinate containing the
1976  // point at the end of the macro element
1977  this->Bottom_flag_pt->locate_zeta(zeta_end_flag, geom_obj_pt, s);
1978 
1979  // Update the pointer to the (sub-)GeomObject within which the
1980  // point at the end of the macro element is located
1981  geom_object_pt[1] = geom_obj_pt;
1982 
1983  // Update second reference value: Reference local coordinate
1984  // in flag sub-element
1985  ref_value[5] = s[0];
1986 
1987  // Kill the existing node update info
1988  node_pt->kill_node_update_info(9);
1989 
1990  // Setup algebraic update for node: Pass update information
1991  node_pt->add_node_update_info(9, // id
1992  this, // mesh
1993  geom_object_pt, // vector of geom objects
1994  ref_value); // vector of ref. values
1995  }
1996 
1997 
1998  if ((id == 1) || (id == 6))
1999  {
2000  // Extract reference values for node update by copy construction
2001  Vector<double> ref_value(node_pt->vector_ref_value());
2002 
2003  // Get zeta coordinate on flag
2004  Vector<double> zeta_flag(1);
2005  zeta_flag[0] = ref_value[2];
2006 
2007  // Get the sub-geomobject and the local coordinate
2008  Vector<double> s(1);
2009  GeomObject* geom_obj_pt;
2010  this->Top_flag_pt->locate_zeta(zeta_flag, geom_obj_pt, s);
2011 
2012  // Extract geometric objects for update by copy construction
2013  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
2014 
2015  // Update the pointer to the (sub-)GeomObject within which the
2016  // reference point is located. (If the flag is simple GeomObject
2017  // this is the same as Leaflet_pt; if it's a compound GeomObject
2018  // this points to the sub-object)
2019  geom_object_pt[0] = geom_obj_pt;
2020 
2021  // Update second reference value: Reference local coordinate
2022  // in flag sub-element
2023  ref_value[3] = s[0];
2024 
2025  if (id == 1)
2026  {
2027  // Kill the existing node update info
2028  node_pt->kill_node_update_info(1);
2029 
2030  // Setup algebraic update for node: Pass update information
2031  node_pt->add_node_update_info(1, // id
2032  this, // mesh
2033  geom_object_pt, // vector of geom objects
2034  ref_value); // vector of ref. values
2035  }
2036  else if (id == 6)
2037  {
2038  // Kill the existing node update info
2039  node_pt->kill_node_update_info(6);
2040 
2041  // Setup algebraic update for node: Pass update information
2042  node_pt->add_node_update_info(6, // id
2043  this, // mesh
2044  geom_object_pt, // vector of geom objects
2045  ref_value); // vector of ref. values
2046  }
2047  }
2048 
2049 
2050  if ((id == 2) || (id == 7))
2051  {
2052  // Extract reference values for node update by copy construction
2053  Vector<double> ref_value(node_pt->vector_ref_value());
2054 
2055  // Get zeta coordinate on flag
2056  Vector<double> zeta_flag(1);
2057  zeta_flag[0] = ref_value[2];
2058 
2059  // Get the sub-geomobject and the local coordinate
2060  Vector<double> s(1);
2061  GeomObject* geom_obj_pt;
2062  this->Bottom_flag_pt->locate_zeta(zeta_flag, geom_obj_pt, s);
2063 
2064  // Extract geometric objects for update by copy construction
2065  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
2066 
2067  // Update the pointer to the (sub-)GeomObject within which the
2068  // reference point is located. (If the flag is simple GeomObject
2069  // this is the same as Leaflet_pt; if it's a compound GeomObject
2070  // this points to the sub-object)
2071  geom_object_pt[0] = geom_obj_pt;
2072 
2073  // Update second reference value: Reference local coordinate
2074  // in flag sub-element
2075  ref_value[3] = s[0];
2076 
2077  if (id == 2)
2078  {
2079  // Kill the existing node update info
2080  node_pt->kill_node_update_info(2);
2081 
2082  // Setup algebraic update for node: Pass update information
2083  node_pt->add_node_update_info(2, // id
2084  this, // mesh
2085  geom_object_pt, // vector of geom objects
2086  ref_value); // vector of ref. values
2087  }
2088  else if (id == 7)
2089  {
2090  // Kill the existing node update info
2091  node_pt->kill_node_update_info(7);
2092 
2093  // Setup algebraic update for node: Pass update information
2094  node_pt->add_node_update_info(7, // id
2095  this, // mesh
2096  geom_object_pt, // vector of geom objects
2097  ref_value); // vector of ref. values
2098  }
2099  }
2100 
2101  if ((id == 3) || (id == 4) || (id == 5))
2102  {
2103  // Extract reference values for node update by copy construction
2104  Vector<double> ref_value(node_pt->vector_ref_value());
2105 
2106  // Get zeta coordinate on flag
2107  Vector<double> zeta_flag(1);
2108  if (id == 3)
2109  {
2110  zeta_flag[0] = ref_value[2];
2111  }
2112  else if (id == 4)
2113  {
2114  zeta_flag[0] = this->Flag_height / 2.;
2115  }
2116  else if (id == 5)
2117  {
2118  zeta_flag[0] = -this->Flag_height / 2.;
2119  }
2120 
2121  // Get the sub-geomobject and the local coordinate
2122  Vector<double> s(1);
2123  GeomObject* geom_obj_pt;
2124  this->Tip_flag_pt->locate_zeta(zeta_flag, geom_obj_pt, s);
2125 
2126  // Extract geometric objects for update by copy construction
2127  Vector<GeomObject*> geom_object_pt(node_pt->vector_geom_object_pt());
2128 
2129  // Update the pointer to the (sub-)GeomObject within which the
2130  // reference point is located. (If the flag is simple GeomObject
2131  // this is the same as Leaflet_pt; if it's a compound GeomObject
2132  // this points to the sub-object)
2133  geom_object_pt[0] = geom_obj_pt;
2134 
2135 
2136  if (id == 3)
2137  {
2138  // Update second reference value: Reference local coordinate
2139  // in flag sub-element
2140  ref_value[3] = s[0];
2141 
2142  // Kill the existing node update info
2143  node_pt->kill_node_update_info(3);
2144 
2145  // Setup algebraic update for node: Pass update information
2146  node_pt->add_node_update_info(3, // id
2147  this, // mesh
2148  geom_object_pt, // vector of geom objects
2149  ref_value); // vector of ref. values
2150  }
2151  else if (id == 4)
2152  {
2153  // Update second reference value: Reference local coordinate
2154  // in flag sub-element
2155  ref_value[2] = s[0];
2156 
2157  // Kill the existing node update info
2158  node_pt->kill_node_update_info(4);
2159 
2160  // Setup algebraic update for node: Pass update information
2161  node_pt->add_node_update_info(4, // id
2162  this, // mesh
2163  geom_object_pt, // vector of geom objects
2164  ref_value); // vector of ref. values
2165  }
2166  else if (id == 5)
2167  {
2168  // Update second reference value: Reference local coordinate
2169  // in flag sub-element
2170  ref_value[2] = s[0];
2171 
2172  // Kill the existing node update info
2173  node_pt->kill_node_update_info(5);
2174 
2175  // Setup algebraic update for node: Pass update information
2176  node_pt->add_node_update_info(5, // id
2177  this, // mesh
2178  geom_object_pt, // vector of geom objects
2179  ref_value); // vector of ref. values
2180  }
2181  }
2182  }
2183 
2184 
2185 } // namespace oomph
2186 
2187 #endif
void node_update_V(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void node_update_II(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void node_update_III(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void algebraic_node_update(const unsigned &t, AlgebraicNode *&node_pt)
Update nodal position at time level t (t=0: present; t>0: previous)
void node_update_IV(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void node_update_I(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void node_update_VIII(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void node_update_VI(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void node_update_VII(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
void setup_algebraic_node_update()
Function to setup the algebraic node update.
void node_update_IX(const unsigned &t, AlgebraicNode *&node_pt)
Helper function.
Domain for cylinder with flag as in Turek benchmark.
CylinderWithFlagMesh(Circle *cylinder_pt, GeomObject *top_flag_pt, GeomObject *bottom_flag_pt, GeomObject *tip_flag_pt, const double &length, const double &height, const double &flag_length, const double &flag_height, const double &centre_x, const double &centre_y, const double &a, TimeStepper *time_stepper_pt=&Mesh::Default_TimeStepper)
Constructor. Pass the pointers to the GeomObjects that parametrise the cylinder, the three edges of t...
void update_node_update(AlgebraicNode *&node_pt)
Update the node update data for specified node following any mesh adapation.
////////////////////////////////////////////////////////////////////// //////////////////////////////...