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