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