cylinder_with_flag_domain.h
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-2025 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_DOMAIN_HEADER
27#define OOMPH_CYLINDER_WITH_FLAG_DOMAIN_HEADER
28
29// Generic includes
30#include "generic/geom_objects.h"
31#include "generic/macro_element.h"
32#include "generic/domain.h"
33
34namespace oomph
35{
36 //===========================================================
37 /// Domain for cylinder with flag as in Turek benchmark.
38 //===========================================================
39 class CylinderWithFlagDomain : public Domain
40 {
41 public:
42 /// Constructor. Pass the pointers to the GeomObjects that parametrise
43 /// the cylinder, the three edges of the flag, the length and height of the
44 /// domain, the length and height of the flag, the coordinates of the
45 /// centre of the cylinder and its radius.
46 CylinderWithFlagDomain(Circle* cylinder_pt,
47 GeomObject* top_flag_pt,
48 GeomObject* bottom_flag_pt,
49 GeomObject* tip_flag_pt,
50 const double& length,
51 const double& height,
52 const double& flag_length,
53 const double& flag_height,
54 const double& centre_x,
55 const double& centre_y,
56 const double& a);
57
58
59 /// Destructor: Emtpy because clean up happens in base class
60 /// as a service to the user!
61 ~CylinderWithFlagDomain() {}
62
63
64 /// Parametrisation of macro element boundaries: f(s) is the position
65 /// vector to macro-element m's boundary in the specified direction
66 /// [N/S/E/W] at the specfied discrete time level (time=0: present; time>0:
67 /// previous)
68 void macro_element_boundary(const unsigned& time,
69 const unsigned& m,
70 const unsigned& direction,
71 const Vector<double>& s,
72 Vector<double>& f);
73
74 /// Access fct to GeomObject (of type Circle)
75 /// that represents the cylinder
76 Circle* cylinder_pt()
77 {
78 return Cylinder_pt;
79 }
80
81 /// Access fct to GeomObjects for top, bottom and tip
82 GeomObject*& bottom_flag_pt()
83 {
84 return Bottom_flag_pt;
85 }
86 GeomObject*& top_flag_pt()
87 {
88 return Top_flag_pt;
89 }
90 GeomObject*& tip_flag_pt()
91 {
92 return Tip_flag_pt;
93 }
94
95 private:
96 /// Helper function to interpolate linearly between the
97 /// "right" and "left" points; \f$ s \in [-1,1] \f$
98 void linear_interpolate(const Vector<double>& left,
99 const Vector<double>& right,
100 const double& s,
101 Vector<double>& f)
102 {
103 for (unsigned i = 0; i < 2; i++)
104 {
105 f[i] = left[i] + (right[i] - left[i]) * 0.5 * (s + 1.0);
106 }
107 }
108
109 // Helper points
110 Vector<double> p1;
111 Vector<double> p2;
112 Vector<double> p3;
113 Vector<double> p4;
114 Vector<double> p5;
115 Vector<double> p6;
116 Vector<double> p7;
117 Vector<double> p8;
118 Vector<double> p9;
119 Vector<double> p10;
120 Vector<double> p11;
121 Vector<double> p12;
122 Vector<double> p13;
123 Vector<double> p14;
124 Vector<double> p15;
125 Vector<double> p16;
126 Vector<double> p17;
127 Vector<double> p18;
128 Vector<double> p19;
129 Vector<double> p20;
130 Vector<double> p21;
131 Vector<double> p22;
132 Vector<double> p23;
133 Vector<double> p24;
134 Vector<double> p25;
135 Vector<double> p26;
136 Vector<double> p27;
137 Vector<double> p28;
138 Vector<double> p29;
139 Vector<double> p30;
140 Vector<double> p31;
141 Vector<double> p32;
142 Vector<double> p33;
143 Vector<double> p34;
144 Vector<double> p35;
145 Vector<double> p36;
146 Vector<double> p37;
147 Vector<double> p38;
148 Vector<double> p39;
149 Vector<double> p40;
150 Vector<double> p41;
151 Vector<double> p42;
152 Vector<double> p43;
153 Vector<double> p44;
154 Vector<double> p45;
155 Vector<double> p46;
156 Vector<double> p47;
157 Vector<double> p48;
158 Vector<double> p49;
159 Vector<double> p50;
160
161 /// Pointer to geometric object that represents the central cylinder
162 Circle* Cylinder_pt;
163
164 /// Pointer to geometric object that represents the top of the flag
165 GeomObject* Top_flag_pt;
166
167 /// Pointer to geometric object that represents the bottom of the flag
168 GeomObject* Bottom_flag_pt;
169
170 /// Pointer to geometric object that represents the tip of the flag
171 GeomObject* Tip_flag_pt;
172
173 // Length of the flag
174 double Lx;
175
176 // Thickness of the flag
177 double Ly;
178
179 // Centre of the cylinder : x coordinate
180 double Centre_x;
181
182 // Centre of the cylinder : y coordinate
183 double Centre_y;
184
185 // Radius of the cylinder
186 double A;
187
188 }; // end of domain
189
190 //=======================================================================
191 /// Constructor, Pass the pointers to the GeomObjects that parametrise
192 /// the cylinder, the three edges of the flag, the length and height of the
193 /// domain, the length and height of the flag, the coordinates of the
194 /// centre of the cylinder and its radius.
195 //=======================================================================
196 inline CylinderWithFlagDomain::CylinderWithFlagDomain(
197 Circle* cylinder_pt,
198 GeomObject* top_flag_pt,
199 GeomObject* bottom_flag_pt,
200 GeomObject* tip_flag_pt,
201 const double& length,
202 const double& height,
203 const double& flag_length,
204 const double& flag_height,
205 const double& centre_x,
206 const double& centre_y,
207 const double& a)
208 : Cylinder_pt(cylinder_pt),
209 Top_flag_pt(top_flag_pt),
210 Bottom_flag_pt(bottom_flag_pt),
211 Tip_flag_pt(tip_flag_pt),
212 Lx(flag_length),
213 Ly(flag_height),
214 Centre_x(centre_x),
215 Centre_y(centre_y),
216 A(a)
217 {
218 // Vertices of rectangle
219 // Those are points of references of the domain
220 // to help create the macro_element_boundary sub_functions
221 p1.resize(2);
222 p1[0] = 0.0;
223 p1[1] = height;
224
225 p2.resize(2);
226 p2[0] = Centre_x;
227 p2[1] = height;
228
229 p3.resize(2);
230 p3[0] = 0.155596 * length;
231 p3[1] = height;
232
233 p4.resize(2);
234 p4[0] = 0.183596 * length;
235 p4[1] = height;
236
237 p5.resize(2);
238 p5[0] = 0.239596 * length;
239 p5[1] = height;
240
241 p6.resize(2);
242 p6[0] = 0.285123967 * length;
243 p6[1] = height;
244
245 p7.resize(2);
246 p7[0] = 0.433884298 * length;
247 p7[1] = height;
248
249 p8.resize(2);
250 p8[0] = 0.578512397 * length;
251 p8[1] = height;
252
253 p9.resize(2);
254 p9[0] = 0.789256198 * length;
255 p9[1] = height;
256
257 p10.resize(2);
258 p10[0] = length;
259 p10[1] = height;
260
261 p11.resize(2);
262 p11[0] = 0.127596 * length;
263 p11[1] = 0.778024390 * height;
264
265 p12.resize(2);
266 p12[0] = 0.155596 * length;
267 p12[1] = 0.778024390 * height;
268
269 p13.resize(2);
270 p13[0] = 0.183596 * length;
271 p13[1] = 0.778024390 * height;
272
273 p14.resize(2);
274 p14[0] = 0.211596 * length;
275 p14[1] = 0.778024390 * height;
276
277 p15.resize(2);
278 p15[0] = 0.285123967 * length;
279 p15[1] = 0.625 * height;
280
281 p16.resize(2);
282 p16[0] = 0.351239669 * length;
283 p16[1] = 0.625 * height;
284
285 p18.resize(2);
286 p18[0] = Centre_x;
287 p18[1] = Centre_y + A;
288
289 p33.resize(2);
290 p33[0] = Centre_x;
291 p33[1] = Centre_y - A;
292
293 p35.resize(2);
294 p35[0] = 0.285123967 * length;
295 p35[1] = 0.350609756 * height;
296
297 p36.resize(2);
298 p36[0] = 0.351239669 * length;
299 p36[1] = 0.350609756 * height;
300
301 p37.resize(2);
302 p37[0] = 0.127596 * length;
303 p37[1] = 0.197585366 * height;
304
305 p38.resize(2);
306 p38[0] = 0.155596 * length;
307 p38[1] = 0.197585366 * height;
308
309 p39.resize(2);
310 p39[0] = 0.183596 * length;
311 p39[1] = 0.197585366 * height;
312
313 p40.resize(2);
314 p40[0] = 0.211596 * length;
315 p40[1] = 0.197585366 * height;
316
317 p41.resize(2);
318 p41[0] = 0.0;
319 p41[1] = 0.;
320
321 p42.resize(2);
322 p42[0] = Centre_x;
323 p42[1] = 0.;
324
325 p43.resize(2);
326 p43[0] = 0.155596 * length;
327 p43[1] = 0.;
328
329 p44.resize(2);
330 p44[0] = 0.183596 * length;
331 p44[1] = 0.;
332
333 p45.resize(2);
334 p45[0] = 0.239596 * length;
335 p45[1] = 0.;
336
337 p46.resize(2);
338 p46[0] = 0.285123967 * length;
339 p46[1] = 0.;
340
341 p47.resize(2);
342 p47[0] = 0.433884298 * length;
343 p47[1] = 0.;
344
345 p48.resize(2);
346 p48[0] = 0.578512397 * length;
347 p48[1] = 0.;
348
349 p49.resize(2);
350 p49[0] = 0.789256198 * length;
351 p49[1] = 0.;
352
353 p50.resize(2);
354 p50[0] = length;
355 p50[1] = 0.;
356
357 // Allocate storage for variable points
358 p21.resize(2);
359 p22.resize(2);
360 p23.resize(2);
361 p24.resize(2);
362 p25.resize(2);
363 p27.resize(2);
364 p28.resize(2);
365 p29.resize(2);
366 p30.resize(2);
367 p31.resize(2);
368
369 // There are 31 macro elements
370 Macro_element_pt.resize(31);
371
372 // Build the 2D macro elements
373 for (unsigned i = 0; i < 31; i++)
374 {
375 Macro_element_pt[i] = new QMacroElement<2>(this, i);
376 }
377
378 } // end of constructor
379
380 //============================================================================
381 /// Parametrisation of macro element boundaries: f(s) is the position
382 /// vector to macro-element m's boundary in the specified direction [N/S/E/W]
383 /// at the specfied discrete time level (time=0: present; time>0: previous)
384 //============================================================================
385 inline void CylinderWithFlagDomain::macro_element_boundary(
386 const unsigned& time,
387 const unsigned& m,
388 const unsigned& direction,
389 const Vector<double>& s,
390 Vector<double>& f)
391 {
392 // Use Quadtree names for directions
393 using namespace QuadTreeNames;
394
395#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
396 // Warn about time argument being moved to the front
397 OomphLibWarning(
398 "Order of function arguments has changed between versions 0.8 and 0.85",
399 "CylinderWithFlagDomain::macro_element_boundary(...)",
400 OOMPH_EXCEPTION_LOCATION);
401#endif
402
403 // Lagrangian coordinate along surface of cylinder
404 Vector<double> xi(1);
405
406 // Point on circle
407 Vector<double> point_on_circle(2);
408
409 // Lagrangian coordinates on the flag
410 Vector<double> zeta(1);
411
412 // Definition of the points that depend on the shape of the flags
413 zeta[0] = 1. / 5. * Lx;
414 Top_flag_pt->position(time, zeta, p21);
415
416 zeta[0] = 2. / 5. * Lx;
417 Top_flag_pt->position(time, zeta, p22);
418
419 zeta[0] = 3. / 5. * Lx;
420 Top_flag_pt->position(time, zeta, p23);
421
422 zeta[0] = 4. / 5. * Lx;
423 Top_flag_pt->position(time, zeta, p24);
424
425 zeta[0] = Ly / 2.;
426 Tip_flag_pt->position(time, zeta, p25);
427
428 zeta[0] = 1. / 5. * Lx;
429 Bottom_flag_pt->position(time, zeta, p27);
430
431 zeta[0] = 2. / 5. * Lx;
432 Bottom_flag_pt->position(time, zeta, p28);
433
434 zeta[0] = 3. / 5. * Lx;
435 Bottom_flag_pt->position(time, zeta, p29);
436
437 zeta[0] = 4. / 5. * Lx;
438 Bottom_flag_pt->position(time, zeta, p30);
439
440 zeta[0] = -Ly / 2.;
441 Tip_flag_pt->position(time, zeta, p31);
442
443 std::ostringstream error_message;
444
445 // Switch on the macro element
446 switch (m)
447 {
448 // Macro element 0, is is immediately left of the cylinder
449 case 0:
450
451 switch (direction)
452 {
453 case N:
454 xi[0] = 3.0 * atan(1.0);
455 Cylinder_pt->position(time, xi, point_on_circle);
456 linear_interpolate(p1, point_on_circle, s[0], f);
457 break;
458
459 case S:
460 xi[0] = -3.0 * atan(1.0);
461 Cylinder_pt->position(time, xi, point_on_circle);
462 linear_interpolate(p41, point_on_circle, s[0], f);
463 break;
464
465 case W:
466 linear_interpolate(p41, p1, s[0], f);
467 break;
468
469 case E:
470 xi[0] = 5.0 * atan(1.0) - 2.0 * atan(1.0) * 0.5 * (1.0 + s[0]);
471 Cylinder_pt->position(time, xi, f);
472 break;
473
474 default:
475 error_message << "Direction is incorrect: " << direction
476 << std::endl;
477 throw OomphLibError(error_message.str(),
478 OOMPH_CURRENT_FUNCTION,
479 OOMPH_EXCEPTION_LOCATION);
480 }
481
482 break;
483
484 // Macro element 1, is immediately above the cylinder
485 case 1:
486
487 switch (direction)
488 {
489 case N:
490 linear_interpolate(p1, p2, s[0], f);
491 break;
492
493 case S:
494 xi[0] = 3.0 * atan(1.0) - atan(1.0) * 0.5 * (1.0 + s[0]);
495 Cylinder_pt->position(time, xi, f);
496 break;
497
498 case W:
499 xi[0] = 3.0 * atan(1.0);
500 Cylinder_pt->position(time, xi, point_on_circle);
501 linear_interpolate(point_on_circle, p1, s[0], f);
502 break;
503
504 case E:
505 linear_interpolate(p18, p2, s[0], f);
506 break;
507
508 default:
509 error_message << "Direction is incorrect: " << direction
510 << std::endl;
511 throw OomphLibError(error_message.str(),
512 OOMPH_CURRENT_FUNCTION,
513 OOMPH_EXCEPTION_LOCATION);
514 }
515
516 break;
517
518 // Macro element 2, is immediately right of the cylinder
519 case 2:
520
521 switch (direction)
522 {
523 case N:
524 linear_interpolate(p2, p11, s[0], f);
525 break;
526
527 case S:
528 xi[0] = 2.0 * atan(1.0) - atan(1.0) * 0.5 * (1.0 + s[0]);
529 Cylinder_pt->position(time, xi, f);
530 break;
531
532 case W:
533 linear_interpolate(p18, p2, s[0], f);
534 break;
535
536 case E:
537 xi[0] = atan(1.0);
538 Cylinder_pt->position(time, xi, point_on_circle);
539 linear_interpolate(point_on_circle, p11, s[0], f);
540 break;
541
542 default:
543 error_message << "Direction is incorrect: " << direction
544 << std::endl;
545 throw OomphLibError(error_message.str(),
546 OOMPH_CURRENT_FUNCTION,
547 OOMPH_EXCEPTION_LOCATION);
548 }
549
550 break;
551
552 // Macro element 3, is immediately below cylinder
553 case 3:
554
555 switch (direction)
556 {
557 case N:
558 xi[0] = atan(1.0);
559 Cylinder_pt->position(time, xi, point_on_circle);
560 linear_interpolate(point_on_circle, p11, s[0], f);
561 break;
562
563 case S:
564 xi[0] = (1. + s[0]) / 2. * 1. / 5. * Lx;
565 Top_flag_pt->position(time, xi, f);
566 break;
567
568 case W:
569 xi[0] = asin(Ly / A / 2.) +
570 (atan(1.0) - asin(Ly / A / 2.)) * 0.5 * (1.0 + s[0]);
571 Cylinder_pt->position(time, xi, f);
572 break;
573
574 case E:
575 linear_interpolate(p21, p11, s[0], f);
576 break;
577
578 default:
579 error_message << "Direction is incorrect: " << direction
580 << std::endl;
581 throw OomphLibError(error_message.str(),
582 OOMPH_CURRENT_FUNCTION,
583 OOMPH_EXCEPTION_LOCATION);
584 }
585
586 break;
587
588 // Macro element 4, is right hand block 1
589 case 4:
590
591 switch (direction)
592 {
593 case N:
594 xi[0] = (1. + s[0]) / 2. * 1. / 5. * Lx;
595 Bottom_flag_pt->position(time, xi, f);
596 break;
597
598 case S:
599 xi[0] = -atan(1.0);
600 Cylinder_pt->position(time, xi, point_on_circle);
601 linear_interpolate(point_on_circle, p37, s[0], f);
602 break;
603
604 case W:
605 xi[0] =
606 -atan(1.0) + (atan(1.0) - asin(Ly / A / 2.)) * 0.5 * (1.0 + s[0]);
607 Cylinder_pt->position(time, xi, f);
608 break;
609
610 case E:
611 linear_interpolate(p37, p27, s[0], f);
612 break;
613
614 default:
615 error_message << "Direction is incorrect: " << direction
616 << std::endl;
617 throw OomphLibError(error_message.str(),
618 OOMPH_CURRENT_FUNCTION,
619 OOMPH_EXCEPTION_LOCATION);
620 }
621
622 break;
623
624 // Macro element 5, is right hand block 2
625 case 5:
626
627 switch (direction)
628 {
629 case N:
630 xi[0] = 6 * atan(1.0) + atan(1.0) * 0.5 * (1.0 + s[0]);
631 Cylinder_pt->position(time, xi, f);
632 break;
633
634 case S:
635 linear_interpolate(p42, p37, s[0], f);
636 break;
637
638 case W:
639 linear_interpolate(p42, p33, s[0], f);
640 break;
641
642 case E:
643 xi[0] = -atan(1.0);
644 Cylinder_pt->position(time, xi, point_on_circle);
645 linear_interpolate(p37, point_on_circle, s[0], f);
646 break;
647
648 default:
649 error_message << "Direction is incorrect: " << direction
650 << std::endl;
651 throw OomphLibError(error_message.str(),
652 OOMPH_CURRENT_FUNCTION,
653 OOMPH_EXCEPTION_LOCATION);
654 }
655
656 break;
657
658 // Macro element 6, is right hand block 3
659 case 6:
660
661 switch (direction)
662 {
663 case N:
664 xi[0] = 5.0 * atan(1.0) + atan(1.0) * 0.5 * (1.0 + s[0]);
665 Cylinder_pt->position(time, xi, f);
666 break;
667
668 case S:
669 linear_interpolate(p41, p42, s[0], f);
670 break;
671
672 case W:
673 xi[0] = 5.0 * atan(1.0);
674 Cylinder_pt->position(time, xi, point_on_circle);
675 linear_interpolate(p41, point_on_circle, s[0], f);
676 break;
677
678 case E:
679 linear_interpolate(p42, p33, s[0], f);
680 break;
681
682 default:
683 error_message << "Direction is incorrect: " << direction
684 << std::endl;
685 throw OomphLibError(error_message.str(),
686 OOMPH_CURRENT_FUNCTION,
687 OOMPH_EXCEPTION_LOCATION);
688 }
689
690 break;
691
692 // Macro element 7, is right hand block 4
693 case 7:
694
695 switch (direction)
696 {
697 case N:
698 linear_interpolate(p2, p3, s[0], f);
699 break;
700
701 case S:
702 linear_interpolate(p11, p12, s[0], f);
703 break;
704
705 case W:
706 linear_interpolate(p11, p2, s[0], f);
707 break;
708
709 case E:
710 linear_interpolate(p12, p3, s[0], f);
711 break;
712
713 default:
714 error_message << "Direction is incorrect: " << direction
715 << std::endl;
716 throw OomphLibError(error_message.str(),
717 OOMPH_CURRENT_FUNCTION,
718 OOMPH_EXCEPTION_LOCATION);
719 }
720
721 break;
722
723 case 8:
724
725 switch (direction)
726 {
727 case N:
728 linear_interpolate(p11, p12, s[0], f);
729 break;
730
731 case S:
732 xi[0] = 1. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
733 Top_flag_pt->position(time, xi, f);
734 break;
735
736 case W:
737 linear_interpolate(p21, p11, s[0], f);
738 break;
739
740 case E:
741 linear_interpolate(p22, p12, s[0], f);
742 break;
743
744 default:
745 error_message << "Direction is incorrect: " << direction
746 << std::endl;
747 throw OomphLibError(error_message.str(),
748 OOMPH_CURRENT_FUNCTION,
749 OOMPH_EXCEPTION_LOCATION);
750 }
751
752 break;
753 case 9:
754
755 switch (direction)
756 {
757 case N:
758 xi[0] = 1. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
759 Bottom_flag_pt->position(time, xi, f);
760 break;
761
762 case S:
763 linear_interpolate(p37, p38, s[0], f);
764 break;
765
766 case W:
767 linear_interpolate(p37, p27, s[0], f);
768 break;
769
770 case E:
771 linear_interpolate(p38, p28, s[0], f);
772 break;
773
774 default:
775 error_message << "Direction is incorrect: " << direction
776 << std::endl;
777 throw OomphLibError(error_message.str(),
778 OOMPH_CURRENT_FUNCTION,
779 OOMPH_EXCEPTION_LOCATION);
780 }
781
782 break;
783
784 case 10:
785
786 switch (direction)
787 {
788 case N:
789 linear_interpolate(p37, p38, s[0], f);
790 break;
791
792 case S:
793 linear_interpolate(p42, p43, s[0], f);
794 break;
795
796 case W:
797 linear_interpolate(p42, p37, s[0], f);
798 break;
799
800 case E:
801 linear_interpolate(p43, p38, s[0], f);
802 break;
803
804 default:
805 error_message << "Direction is incorrect: " << direction
806 << std::endl;
807 throw OomphLibError(error_message.str(),
808 OOMPH_CURRENT_FUNCTION,
809 OOMPH_EXCEPTION_LOCATION);
810 }
811
812 break;
813 case 11:
814
815 switch (direction)
816 {
817 case N:
818 linear_interpolate(p3, p4, s[0], f);
819 break;
820
821 case S:
822 linear_interpolate(p12, p13, s[0], f);
823 break;
824
825 case W:
826 linear_interpolate(p12, p3, s[0], f);
827 break;
828
829 case E:
830 linear_interpolate(p13, p4, s[0], f);
831 break;
832
833 default:
834 error_message << "Direction is incorrect: " << direction
835 << std::endl;
836 throw OomphLibError(error_message.str(),
837 OOMPH_CURRENT_FUNCTION,
838 OOMPH_EXCEPTION_LOCATION);
839 }
840
841 break;
842
843 case 12:
844
845 switch (direction)
846 {
847 case N:
848 linear_interpolate(p12, p13, s[0], f);
849 break;
850
851 case S:
852 // linear_interpolate(p22,p23,s[0],f);
853 xi[0] = 2. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
854 Top_flag_pt->position(time, xi, f);
855 break;
856
857 case W:
858 linear_interpolate(p22, p12, s[0], f);
859 break;
860
861 case E:
862 linear_interpolate(p23, p13, s[0], f);
863 break;
864
865 default:
866 error_message << "Direction is incorrect: " << direction
867 << std::endl;
868 throw OomphLibError(error_message.str(),
869 OOMPH_CURRENT_FUNCTION,
870 OOMPH_EXCEPTION_LOCATION);
871 }
872
873 break;
874
875 case 13:
876
877 switch (direction)
878 {
879 case N:
880 xi[0] = 2. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
881 Bottom_flag_pt->position(time, xi, f);
882 break;
883
884 case S:
885 linear_interpolate(p38, p39, s[0], f);
886 break;
887
888 case W:
889 linear_interpolate(p38, p28, s[0], f);
890 break;
891
892 case E:
893 linear_interpolate(p39, p29, s[0], f);
894 break;
895
896 default:
897 error_message << "Direction is incorrect: " << direction
898 << std::endl;
899 throw OomphLibError(error_message.str(),
900 OOMPH_CURRENT_FUNCTION,
901 OOMPH_EXCEPTION_LOCATION);
902 }
903
904 break;
905
906 case 14:
907
908 switch (direction)
909 {
910 case N:
911 linear_interpolate(p38, p39, s[0], f);
912 break;
913
914 case S:
915 linear_interpolate(p43, p44, s[0], f);
916 break;
917
918 case W:
919 linear_interpolate(p43, p38, s[0], f);
920 break;
921
922 case E:
923 linear_interpolate(p44, p39, s[0], f);
924 break;
925
926 default:
927 error_message << "Direction is incorrect: " << direction
928 << std::endl;
929 throw OomphLibError(error_message.str(),
930 OOMPH_CURRENT_FUNCTION,
931 OOMPH_EXCEPTION_LOCATION);
932 }
933
934 break;
935
936 case 15:
937
938 switch (direction)
939 {
940 case N:
941 linear_interpolate(p4, p5, s[0], f);
942 break;
943
944 case S:
945 linear_interpolate(p13, p14, s[0], f);
946 break;
947
948 case W:
949 linear_interpolate(p13, p4, s[0], f);
950 break;
951
952 case E:
953 linear_interpolate(p14, p5, s[0], f);
954 break;
955
956 default:
957 error_message << "Direction is incorrect: " << direction
958 << std::endl;
959 throw OomphLibError(error_message.str(),
960 OOMPH_CURRENT_FUNCTION,
961 OOMPH_EXCEPTION_LOCATION);
962 }
963
964 break;
965
966 case 16:
967
968 switch (direction)
969 {
970 case N:
971 linear_interpolate(p13, p14, s[0], f);
972 break;
973
974 case S:
975 xi[0] = 3. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
976 Top_flag_pt->position(time, xi, f);
977 break;
978
979 case W:
980 linear_interpolate(p23, p13, s[0], f);
981 break;
982
983 case E:
984 linear_interpolate(p24, p14, s[0], f);
985 break;
986
987 default:
988 error_message << "Direction is incorrect: " << direction
989 << std::endl;
990 throw OomphLibError(error_message.str(),
991 OOMPH_CURRENT_FUNCTION,
992 OOMPH_EXCEPTION_LOCATION);
993 }
994
995 break;
996
997 case 17:
998
999 switch (direction)
1000 {
1001 case N:
1002 xi[0] = 3. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
1003 Bottom_flag_pt->position(time, xi, f);
1004 break;
1005
1006 case S:
1007 linear_interpolate(p39, p40, s[0], f);
1008 break;
1009
1010 case W:
1011 linear_interpolate(p39, p29, s[0], f);
1012 break;
1013
1014 case E:
1015 linear_interpolate(p40, p30, s[0], f);
1016 break;
1017
1018 default:
1019 error_message << "Direction is incorrect: " << direction
1020 << std::endl;
1021 throw OomphLibError(error_message.str(),
1022 OOMPH_CURRENT_FUNCTION,
1023 OOMPH_EXCEPTION_LOCATION);
1024 }
1025
1026 break;
1027
1028 case 18:
1029
1030 switch (direction)
1031 {
1032 case N:
1033 linear_interpolate(p39, p40, s[0], f);
1034 break;
1035
1036 case S:
1037 linear_interpolate(p44, p45, s[0], f);
1038 break;
1039
1040 case W:
1041 linear_interpolate(p44, p39, s[0], f);
1042 break;
1043
1044 case E:
1045 linear_interpolate(p45, p40, s[0], f);
1046 break;
1047
1048 default:
1049 error_message << "Direction is incorrect: " << direction
1050 << std::endl;
1051 throw OomphLibError(error_message.str(),
1052 OOMPH_CURRENT_FUNCTION,
1053 OOMPH_EXCEPTION_LOCATION);
1054 }
1055
1056 break;
1057
1058 case 19:
1059
1060 switch (direction)
1061 {
1062 case N:
1063 linear_interpolate(p14, p5, s[0], f);
1064 break;
1065
1066 case S:
1067 xi[0] = 4. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
1068 Top_flag_pt->position(time, xi, f);
1069 break;
1070
1071 case W:
1072 linear_interpolate(p24, p14, s[0], f);
1073 break;
1074
1075 case E:
1076 linear_interpolate(p25, p5, s[0], f);
1077 break;
1078
1079 default:
1080 error_message << "Direction is incorrect: " << direction
1081 << std::endl;
1082 throw OomphLibError(error_message.str(),
1083 OOMPH_CURRENT_FUNCTION,
1084 OOMPH_EXCEPTION_LOCATION);
1085 }
1086
1087 break;
1088
1089 case 20:
1090
1091 switch (direction)
1092 {
1093 case N:
1094 xi[0] = 4. / 5. * Lx + (1. + s[0]) / 2. * 1. / 5. * Lx;
1095 Bottom_flag_pt->position(time, xi, f);
1096 break;
1097
1098 case S:
1099 linear_interpolate(p40, p45, s[0], f);
1100 break;
1101
1102 case W:
1103 linear_interpolate(p40, p30, s[0], f);
1104 break;
1105
1106 case E:
1107 linear_interpolate(p45, p31, s[0], f);
1108 break;
1109
1110 default:
1111 error_message << "Direction is incorrect: " << direction
1112 << std::endl;
1113 throw OomphLibError(error_message.str(),
1114 OOMPH_CURRENT_FUNCTION,
1115 OOMPH_EXCEPTION_LOCATION);
1116 }
1117
1118 break;
1119
1120 case 21:
1121
1122 switch (direction)
1123 {
1124 case N:
1125 linear_interpolate(p5, p6, s[0], f);
1126 break;
1127
1128 case S:
1129 linear_interpolate(p25, p15, s[0], f);
1130 break;
1131
1132 case W:
1133 linear_interpolate(p25, p5, s[0], f);
1134 break;
1135
1136 case E:
1137 linear_interpolate(p15, p6, s[0], f);
1138 break;
1139
1140 default:
1141 error_message << "Direction is incorrect: " << direction
1142 << std::endl;
1143 throw OomphLibError(error_message.str(),
1144 OOMPH_CURRENT_FUNCTION,
1145 OOMPH_EXCEPTION_LOCATION);
1146 }
1147
1148 break;
1149
1150 case 22:
1151
1152 switch (direction)
1153 {
1154 case N:
1155 linear_interpolate(p25, p15, s[0], f);
1156 break;
1157
1158 case S:
1159 linear_interpolate(p31, p35, s[0], f);
1160 break;
1161
1162 case W:
1163 xi[0] = s[0] * Ly / 2.;
1164 Tip_flag_pt->position(time, xi, f);
1165 break;
1166
1167 case E:
1168 linear_interpolate(p35, p15, s[0], f);
1169 break;
1170
1171 default:
1172 error_message << "Direction is incorrect: " << direction
1173 << std::endl;
1174 throw OomphLibError(error_message.str(),
1175 OOMPH_CURRENT_FUNCTION,
1176 OOMPH_EXCEPTION_LOCATION);
1177 }
1178
1179 break;
1180
1181 case 23:
1182
1183 switch (direction)
1184 {
1185 case N:
1186 linear_interpolate(p31, p35, s[0], f);
1187 break;
1188
1189 case S:
1190 linear_interpolate(p45, p46, s[0], f);
1191 break;
1192
1193 case W:
1194 linear_interpolate(p45, p31, s[0], f);
1195 break;
1196
1197 case E:
1198 linear_interpolate(p46, p35, s[0], f);
1199 break;
1200
1201 default:
1202 error_message << "Direction is incorrect: " << direction
1203 << std::endl;
1204 throw OomphLibError(error_message.str(),
1205 OOMPH_CURRENT_FUNCTION,
1206 OOMPH_EXCEPTION_LOCATION);
1207 }
1208
1209 break;
1210
1211 case 24:
1212
1213 switch (direction)
1214 {
1215 case N:
1216 linear_interpolate(p6, p7, s[0], f);
1217 break;
1218
1219 case S:
1220 linear_interpolate(p15, p16, s[0], f);
1221 break;
1222
1223 case W:
1224 linear_interpolate(p15, p6, s[0], f);
1225 break;
1226
1227 case E:
1228 linear_interpolate(p16, p7, s[0], f);
1229 break;
1230
1231 default:
1232 error_message << "Direction is incorrect: " << direction
1233 << std::endl;
1234 throw OomphLibError(error_message.str(),
1235 OOMPH_CURRENT_FUNCTION,
1236 OOMPH_EXCEPTION_LOCATION);
1237 }
1238
1239 break;
1240
1241 case 25:
1242
1243 switch (direction)
1244 {
1245 case N:
1246 linear_interpolate(p15, p16, s[0], f);
1247 break;
1248
1249 case S:
1250 linear_interpolate(p35, p36, s[0], f);
1251 break;
1252
1253 case W:
1254 linear_interpolate(p35, p15, s[0], f);
1255 break;
1256
1257 case E:
1258 linear_interpolate(p36, p16, s[0], f);
1259 break;
1260
1261 default:
1262 error_message << "Direction is incorrect: " << direction
1263 << std::endl;
1264 throw OomphLibError(error_message.str(),
1265 OOMPH_CURRENT_FUNCTION,
1266 OOMPH_EXCEPTION_LOCATION);
1267 }
1268
1269 break;
1270
1271 case 26:
1272
1273 switch (direction)
1274 {
1275 case N:
1276 linear_interpolate(p35, p36, s[0], f);
1277 break;
1278
1279 case S:
1280 linear_interpolate(p46, p47, s[0], f);
1281 break;
1282
1283 case W:
1284 linear_interpolate(p46, p35, s[0], f);
1285 break;
1286
1287 case E:
1288 linear_interpolate(p47, p36, s[0], f);
1289 break;
1290
1291 default:
1292 error_message << "Direction is incorrect: " << direction
1293 << std::endl;
1294 throw OomphLibError(error_message.str(),
1295 OOMPH_CURRENT_FUNCTION,
1296 OOMPH_EXCEPTION_LOCATION);
1297 }
1298
1299 break;
1300
1301 case 27:
1302
1303 switch (direction)
1304 {
1305 case N:
1306 linear_interpolate(p16, p7, s[0], f);
1307 break;
1308
1309 case S:
1310 linear_interpolate(p36, p47, s[0], f);
1311 break;
1312
1313 case W:
1314 linear_interpolate(p36, p16, s[0], f);
1315 break;
1316
1317 case E:
1318 linear_interpolate(p47, p7, s[0], f);
1319 break;
1320
1321 default:
1322 error_message << "Direction is incorrect: " << direction
1323 << std::endl;
1324 throw OomphLibError(error_message.str(),
1325 OOMPH_CURRENT_FUNCTION,
1326 OOMPH_EXCEPTION_LOCATION);
1327 }
1328
1329 break;
1330
1331 case 28:
1332
1333 switch (direction)
1334 {
1335 case N:
1336 linear_interpolate(p7, p8, s[0], f);
1337 break;
1338
1339 case S:
1340 linear_interpolate(p47, p48, s[0], f);
1341 break;
1342
1343 case W:
1344 linear_interpolate(p47, p7, s[0], f);
1345 break;
1346
1347 case E:
1348 linear_interpolate(p48, p8, s[0], f);
1349 break;
1350
1351 default:
1352 error_message << "Direction is incorrect: " << direction
1353 << std::endl;
1354 throw OomphLibError(error_message.str(),
1355 OOMPH_CURRENT_FUNCTION,
1356 OOMPH_EXCEPTION_LOCATION);
1357 }
1358
1359 break;
1360
1361 case 29:
1362
1363 switch (direction)
1364 {
1365 case N:
1366 linear_interpolate(p8, p9, s[0], f);
1367 break;
1368
1369 case S:
1370 linear_interpolate(p48, p49, s[0], f);
1371 break;
1372
1373 case W:
1374 linear_interpolate(p48, p8, s[0], f);
1375 break;
1376
1377 case E:
1378 linear_interpolate(p49, p9, s[0], f);
1379 break;
1380
1381 default:
1382 error_message << "Direction is incorrect: " << direction
1383 << std::endl;
1384 throw OomphLibError(error_message.str(),
1385 OOMPH_CURRENT_FUNCTION,
1386 OOMPH_EXCEPTION_LOCATION);
1387 }
1388
1389 break;
1390
1391 case 30:
1392
1393 switch (direction)
1394 {
1395 case N:
1396 linear_interpolate(p9, p10, s[0], f);
1397 break;
1398
1399 case S:
1400 linear_interpolate(p49, p50, s[0], f);
1401 break;
1402
1403 case W:
1404 linear_interpolate(p49, p9, s[0], f);
1405 break;
1406
1407 case E:
1408 linear_interpolate(p50, p10, s[0], f);
1409 break;
1410
1411 default:
1412 error_message << "Direction is incorrect: " << direction
1413 << std::endl;
1414 throw OomphLibError(error_message.str(),
1415 OOMPH_CURRENT_FUNCTION,
1416 OOMPH_EXCEPTION_LOCATION);
1417 }
1418
1419 break;
1420
1421 default:
1422
1423 error_message << "Wrong macro element number" << m << std::endl;
1424 throw OomphLibError(error_message.str(),
1425 OOMPH_CURRENT_FUNCTION,
1426 OOMPH_EXCEPTION_LOCATION);
1427 }
1428
1429 } // end of macro element boundary
1430} // namespace oomph
1431
1432#endif
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
oomph::Circle
Circle in 2D space.
Definition geom_objects.h:873
oomph::Circle::position
void position(const Vector< double > &zeta, Vector< double > &r) const
Position Vector at Lagrangian coordinate zeta.
Definition geom_objects.h:1024
oomph::CylinderWithFlagDomain
Domain for cylinder with flag as in Turek benchmark.
Definition cylinder_with_flag_domain.h:40
oomph::CylinderWithFlagDomain::Tip_flag_pt
GeomObject * Tip_flag_pt
Pointer to geometric object that represents the tip of the flag.
Definition cylinder_with_flag_domain.h:171
oomph::CylinderWithFlagDomain::macro_element_boundary
void macro_element_boundary(const unsigned &time, const unsigned &m, const unsigned &direction, const Vector< double > &s, Vector< double > &f)
Parametrisation of macro element boundaries: f(s) is the position vector to macro-element m's boundar...
Definition cylinder_with_flag_domain.h:385
oomph::CylinderWithFlagDomain::cylinder_pt
Circle * cylinder_pt()
Access fct to GeomObject (of type Circle) that represents the cylinder.
Definition cylinder_with_flag_domain.h:76
oomph::CylinderWithFlagDomain::~CylinderWithFlagDomain
~CylinderWithFlagDomain()
Destructor: Emtpy because clean up happens in base class as a service to the user!
Definition cylinder_with_flag_domain.h:61
oomph::CylinderWithFlagDomain::bottom_flag_pt
GeomObject *& bottom_flag_pt()
Access fct to GeomObjects for top, bottom and tip.
Definition cylinder_with_flag_domain.h:82
oomph::CylinderWithFlagDomain::Top_flag_pt
GeomObject * Top_flag_pt
Pointer to geometric object that represents the top of the flag.
Definition cylinder_with_flag_domain.h:165
oomph::CylinderWithFlagDomain::linear_interpolate
void linear_interpolate(const Vector< double > &left, const Vector< double > &right, const double &s, Vector< double > &f)
Helper function to interpolate linearly between the "right" and "left" points; .
Definition cylinder_with_flag_domain.h:98
oomph::CylinderWithFlagDomain::CylinderWithFlagDomain
CylinderWithFlagDomain(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)
Constructor. Pass the pointers to the GeomObjects that parametrise the cylinder, the three edges of t...
Definition cylinder_with_flag_domain.h:196
oomph::CylinderWithFlagDomain::Bottom_flag_pt
GeomObject * Bottom_flag_pt
Pointer to geometric object that represents the bottom of the flag.
Definition cylinder_with_flag_domain.h:168
oomph::CylinderWithFlagDomain::Cylinder_pt
Circle * Cylinder_pt
Pointer to geometric object that represents the central cylinder.
Definition cylinder_with_flag_domain.h:162
oomph::Domain
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition domain.h:67
oomph::Domain::Macro_element_pt
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition domain.h:301
oomph::GeomObject
A geometric object is an object that provides a parametrised description of its shape via the functio...
Definition geom_objects.h:101
oomph::GeomObject::position
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::OomphLibWarning
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Definition oomph_definitions.h:274
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
geom_objects.h
macro_element.h
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30