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-2022 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
33namespace 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
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
e
Definition: cfortran.h:571
static char t char * s
Definition: cfortran.h:568
cstr elem_len * i
Definition: cfortran.h:603
char t
Definition: cfortran.h:568
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.
////////////////////////////////////////////////////////////////////
Vector< GeomObject * > & vector_geom_object_pt(const int &id)
Return vector of geometric objects involved in id-th update function.
int node_update_fct_id()
Default (usually first if there are multiple ones) node update fct id.
void kill_node_update_info(const int &id=0)
Erase algebraic node update information for id-th node update function. Id defaults to 0.
Vector< double > & vector_ref_value()
Return vector of reference values involved in default (usually first) update function.
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 ...
////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:873
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...
A general Finite Element class.
Definition: elements.h:1313
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:1842
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2210
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2175
/////////////////////////////////////////////////////////////////////
Definition: geom_objects.h:101
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
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:381
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:906
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:2394
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:1060
An OomphLibError object which should be thrown when an run-time error is encountered....
void update_node_update(AlgebraicNode *&node_pt)
Update the node update data for specified node following any mesh adapation.
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: timesteppers.h:231
//////////////////////////////////////////////////////////////////// ////////////////////////////////...