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-2023 Matthias Heil and Andrew Hazel
7 // LIC//
8 // LIC// This library is free software; you can redistribute it and/or
9 // LIC// modify it under the terms of the GNU Lesser General Public
10 // LIC// License as published by the Free Software Foundation; either
11 // LIC// version 2.1 of the License, or (at your option) any later version.
12 // LIC//
13 // LIC// This library is distributed in the hope that it will be useful,
14 // LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 // LIC// Lesser General Public License for more details.
17 // LIC//
18 // LIC// You should have received a copy of the GNU Lesser General Public
19 // LIC// License along with this library; if not, write to the Free Software
20 // LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 // LIC// 02110-1301 USA.
22 // LIC//
23 // LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24 // LIC//
25 // LIC//====================================================================
26 #include "cylinder_with_flag_domain.h"
27 
28 
29 namespace 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
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
cylinder_with_flag_domain.h
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition: annular_domain.h:35