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