quarter_tube_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// 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//====================================================================
30 //Include guards
31 #ifndef OOMPH_QUARTER_TUBE_DOMAIN_HEADER
32 #define OOMPH_QUARTER_TUBE_DOMAIN_HEADER
33 
34 // Generic oomph-lib includes
35 #include "../generic/quadtree.h"
36 #include "../generic/domain.h"
37 #include "../generic/geom_objects.h"
38 
39 namespace oomph
40 {
41 
42 
43 
44 
45 //=================================================================
46 /// \short Quarter tube as domain. Domain is bounded by
47 /// curved boundary which is represented by a GeomObject. Domain is
48 /// parametrised by three macro elements in each of the nlayer slices
49 //=================================================================
50 class QuarterTubeDomain : public Domain
51 {
52 
53 public:
54 
55  /// \short Constructor: Pass boundary object and start and end coordinates
56  /// and fraction along boundary object where outer ring is divided.
57  /// We form nlayer axial slices.
58  QuarterTubeDomain(GeomObject* boundary_geom_object_pt,
59  const Vector<double>& xi_lo,
60  const double& fract_mid,
61  const Vector<double>& xi_hi,
62  const unsigned& nlayer) :
63  Xi_lo(xi_lo), Fract_mid(fract_mid), Xi_hi(xi_hi), Nlayer(nlayer),
64  Wall_pt(boundary_geom_object_pt), BL_squash_fct_pt(&default_BL_squash_fct),
66  {
67  // There are three macro elements
68  unsigned nmacro=3*nlayer;
69 
70  // Resize
71  Macro_element_pt.resize(nmacro);
72 
73  // Create macro elements
74  for (unsigned i=0;i<nmacro;i++)
75  {
77  }
78  }
79 
80  /// Broken copy constructor
82  {
83  BrokenCopy::broken_copy("QuarterTubeDomain");
84  }
85 
86  /// Broken assignment operator
88  {
89  BrokenCopy::broken_assign("QuarterTubeDomain");
90  }
91 
92 
93  /// Destructor: Kill macro elements
95  {
96  for (unsigned i=0;i<3*Nlayer;i++)
97  {
98  delete Macro_element_pt[i];
99  }
100  }
101 
102  /// \short Typedef for function pointer for function that squashes
103  /// the outer two macro elements towards
104  /// the wall by mapping the input value of the "radial" macro element
105  /// coordinate to the return value
106  typedef double (*BLSquashFctPt)(const double& s);
107 
108 
109  /// \short Function pointer for function that squashes
110  /// the outer two macro elements towards
111  /// the wall by mapping the input value of the "radial" macro element
112  /// coordinate to the return value
114  {
115  return BL_squash_fct_pt;
116  }
117 
118 
119  /// \short Function that squashes the outer two macro elements towards
120  /// the wall by mapping the input value of the "radial" macro element
121  /// coordinate to the return value.
122  double s_squashed(const double& s)
123  {
124  return BL_squash_fct_pt(s);
125  }
126 
127 
128  /// \short Typedef for function pointer for function that implements
129  /// axial spacing of macro elements
130  typedef double (*AxialSpacingFctPt)(const double& xi);
131 
132 
133  /// \short Function pointer for function that implements
134  /// axial spacing of macro elements
136  {
137  return Axial_spacing_fct_pt;
138  }
139 
140  /// \short Function that implements
141  /// axial spacing of macro elements
142  double axial_spacing_fct(const double& xi)
143  {
144  return Axial_spacing_fct_pt(xi);
145  }
146 
147 
148  /// \short Vector representation of the i_macro-th macro element
149  /// boundary i_direct (L/R/D/U/B/F) at time level t
150  /// (t=0: present; t>0: previous):
151  /// f(s).
152  void macro_element_boundary(const unsigned& t,
153  const unsigned& i_macro,
154  const unsigned& i_direct,
155  const Vector<double>& s,
156  Vector<double>& f);
157 
158 private:
159 
160  /// Lower limit for the coordinates along the wall
162 
163  /// Fraction along wall where outer ring is to be divided
164  double Fract_mid;
165 
166  /// Upper limit for the coordinates along the wall
168 
169  /// Number of layers
170  unsigned Nlayer;
171 
172  /// Pointer to geometric object that represents the curved wall
174 
175 
176  /// \short Function pointer for function that squashes
177  /// the outer two macro elements towards
178  /// the wall by mapping the input value of the "radial" macro element
179  /// coordinate to the return value
181 
182 
183  /// \short Default for function that squashes
184  /// the outer two macro elements towards
185  /// the wall by mapping the input value of the "radial" macro element
186  /// coordinate to the return value: Identity.
187  static double default_BL_squash_fct(const double& s)
188  {
189  return s;
190  }
191 
192 
193  /// \short Function pointer for function that implements
194  /// axial spacing of macro elements
196 
197 
198  /// \short Default for function that implements
199  /// axial spacing of macro elements
200  static double default_axial_spacing_fct(const double& xi)
201  {
202  return xi;
203  }
204 
205 
206  /// \short Boundary of central box macro element in layer i_layer
207  /// zeta \f$ \in [-1,1]^2 \f$
208  void r_centr_L(const unsigned& t, const Vector<double>& zeta,
209  const unsigned& i_layer, Vector<double>& f);
210 
211  /// \short Boundary of central box macro element in layer i_layer
212  /// zeta \f$ \in [-1,1]^2 \f$
213  void r_centr_R(const unsigned& t, const Vector<double>& zeta,
214  const unsigned& i_layer, Vector<double>& f);
215 
216  /// \short Boundary of central box macro element in layer i_layer
217  /// zeta \f$ \in [-1,1]^2 \f$
218  void r_centr_D(const unsigned& t, const Vector<double>& zeta,
219  const unsigned& i_layer, Vector<double>& f);
220 
221  /// \short Boundary of central box macro element in layer i_layer
222  /// zeta \f$ \in [-1,1]^2 \f$
223  void r_centr_U(const unsigned& t, const Vector<double>& zeta,
224  const unsigned& i_layer, Vector<double>& f);
225 
226  /// \short Boundary of central box macro element in layer i_layer
227  /// zeta \f$ \in [-1,1]^2 \f$
228  void r_centr_B(const unsigned& t, const Vector<double>& zeta,
229  const unsigned& i_layer, Vector<double>& f);
230 
231  /// \short Boundary of central box macro element in layer i_layer
232  /// zeta \f$ \in [-1,1]^2 \f$
233  void r_centr_F(const unsigned& t, const Vector<double>& zeta,
234  const unsigned& i_layer, Vector<double>& f);
235 
236 
237 
238 
239 
240  /// \short Boundary of bottom right box macro element in layer i_layer
241  /// zeta \f$ \in [-1,1]^2 \f$
242  void r_bot_right_L(const unsigned& t, const Vector<double>& zeta,
243  const unsigned& i_layer, Vector<double>& f);
244 
245  /// \short Boundary of bottom right box macro element in layer i_layer
246  /// zeta \f$ \in [-1,1]^2 \f$
247  void r_bot_right_R(const unsigned& t, const Vector<double>& zeta,
248  const unsigned& i_layer, Vector<double>& f);
249 
250  /// \short Boundary of bottom right box macro element in layer i_layer
251  /// zeta \f$ \in [-1,1]^2 \f$
252  void r_bot_right_D(const unsigned& t, const Vector<double>& zeta,
253  const unsigned& i_layer, Vector<double>& f);
254 
255  /// \short Boundary of bottom right box macro element in layer i_layer
256  /// zeta \f$ \in [-1,1]^2 \f$
257  void r_bot_right_U(const unsigned& t, const Vector<double>& zeta,
258  const unsigned& i_layer, Vector<double>& f);
259 
260  /// \short Boundary of bottom right box macro element in layer i_layer
261  /// zeta \f$ \in [-1,1]^2 \f$
262  void r_bot_right_B(const unsigned& t, const Vector<double>& zeta,
263  const unsigned& i_layer, Vector<double>& f);
264 
265  /// \short Boundary of bottom right box macro element in layer i_layer
266  /// zeta \f$ \in [-1,1]^2 \f$
267  void r_bot_right_F(const unsigned& t, const Vector<double>& zeta,
268  const unsigned& i_layer, Vector<double>& f);
269 
270 
271 
272 
273 
274  /// \short Boundary of top left box macro element in layer i_layer
275  /// zeta \f$ \in [-1,1]^2 \f$
276  void r_top_left_L(const unsigned& t, const Vector<double>& zeta,
277  const unsigned& i_layer, Vector<double>& f);
278 
279  /// \short Boundary of top left box macro element in layer i_layer
280  /// zeta \f$ \in [-1,1]^2 \f$
281  void r_top_left_R(const unsigned& t, const Vector<double>& zeta,
282  const unsigned& i_layer, Vector<double>& f);
283 
284  /// \short Boundary of top left box macro element in layer i_layer
285  /// zeta \f$ \in [-1,1]^2 \f$
286  void r_top_left_D(const unsigned& t, const Vector<double>& zeta,
287  const unsigned& i_layer, Vector<double>& f);
288 
289  /// \short Boundary of top left box macro element in layer i_layer
290  /// zeta \f$ \in [-1,1]^2 \f$
291  void r_top_left_U(const unsigned& t, const Vector<double>& zeta,
292  const unsigned& i_layer, Vector<double>& f);
293 
294  /// \short Boundary of top left box macro element in layer i_layer
295  /// zeta \f$ \in [-1,1]^2 \f$
296  void r_top_left_B(const unsigned& t, const Vector<double>& zeta,
297  const unsigned& i_layer, Vector<double>& f);
298 
299  /// \short Boundary of top left box macro element in layer i_layer
300  /// zeta \f$ \in [-1,1]^2 \f$
301  void r_top_left_F(const unsigned& t, const Vector<double>& zeta,
302  const unsigned& i_layer, Vector<double>& f);
303 
304 };
305 
306 
307 /////////////////////////////////////////////////////////////////////////
308 /////////////////////////////////////////////////////////////////////////
309 /////////////////////////////////////////////////////////////////////////
310 
311 
312 
313 
314 //=================================================================
315 /// \short Vector representation of the imacro-th macro element
316 /// boundary idirect (L/R/D/U/B/F) at time level t
317 /// (t=0: present; t>0: previous): f(s)
318 //=================================================================
320 const unsigned& t,
321 const unsigned& imacro,
322 const unsigned& idirect,
323 const Vector<double>& s,
324 Vector<double>& f)
325 {
326 
327  using namespace OcTreeNames;
328 
329 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
330  // Warn about time argument being moved to the front
332  "Order of function arguments has changed between versions 0.8 and 0.85",
333  "QuarterTubeDomain::macro_element_boundary(...)",
334  OOMPH_EXCEPTION_LOCATION);
335 #endif
336 
337 
338  unsigned ilayer=unsigned(imacro/3);
339 
340  // Which macro element?
341  // --------------------
342  switch(imacro%3)
343  {
344 
345  // Macro element 0: Central box
346  case 0:
347 
348  // Which direction?
349  if (idirect==L)
350  {
351  r_centr_L(t,s,ilayer,f);
352  }
353  else if (idirect==R)
354  {
355  r_centr_R(t,s,ilayer,f);
356  }
357  else if (idirect==D)
358  {
359  r_centr_D(t,s,ilayer,f);
360  }
361  else if (idirect==U)
362  {
363  r_centr_U(t,s,ilayer,f);
364  }
365  else if (idirect==B)
366  {
367  r_centr_B(t,s,ilayer,f);
368  }
369  else if (idirect==F)
370  {
371  r_centr_F(t,s,ilayer,f);
372  }
373  else
374  {
375 
376  std::ostringstream error_stream;
377  error_stream << "idirect is " << idirect
378  << " not one of L, R, D, U, B, F" << std::endl;
379 
380  throw OomphLibError(
381  error_stream.str(),
382  OOMPH_CURRENT_FUNCTION,
383  OOMPH_EXCEPTION_LOCATION);
384  }
385 
386  break;
387 
388 
389  // Macro element 1: Bottom right
390  case 1:
391 
392  // Which direction?
393  if (idirect==L)
394  {
395  r_bot_right_L(t,s,ilayer,f);
396  }
397  else if (idirect==R)
398  {
399  r_bot_right_R(t,s,ilayer,f);
400  }
401  else if (idirect==D)
402  {
403  r_bot_right_D(t,s,ilayer,f);
404  }
405  else if (idirect==U)
406  {
407  r_bot_right_U(t,s,ilayer,f);
408  }
409  else if (idirect==B)
410  {
411  r_bot_right_B(t,s,ilayer,f);
412  }
413  else if (idirect==F)
414  {
415  r_bot_right_F(t,s,ilayer,f);
416  }
417  else
418  {
419 
420  std::ostringstream error_stream;
421  error_stream << "idirect is " << idirect
422  << " not one of L, R, D, U, B, F" << std::endl;
423 
424  throw OomphLibError(
425  error_stream.str(),
426  OOMPH_CURRENT_FUNCTION,
427  OOMPH_EXCEPTION_LOCATION);
428  }
429 
430 
431  break;
432 
433  // Macro element 2:Top left
434  case 2:
435 
436  // Which direction?
437  if (idirect==L)
438  {
439  r_top_left_L(t,s,ilayer,f);
440  }
441  else if (idirect==R)
442  {
443  r_top_left_R(t,s,ilayer,f);
444  }
445  else if (idirect==D)
446  {
447  r_top_left_D(t,s,ilayer,f);
448  }
449  else if (idirect==U)
450  {
451  r_top_left_U(t,s,ilayer,f);
452  }
453  else if (idirect==B)
454  {
455  r_top_left_B(t,s,ilayer,f);
456  }
457  else if (idirect==F)
458  {
459  r_top_left_F(t,s,ilayer,f);
460  }
461  else
462  {
463 
464  std::ostringstream error_stream;
465  error_stream << "idirect is " << idirect
466  << " not one of L, R, D, U, B, F" << std::endl;
467 
468  throw OomphLibError(
469  error_stream.str(),
470  OOMPH_CURRENT_FUNCTION,
471  OOMPH_EXCEPTION_LOCATION);
472  }
473 
474  break;
475 
476  default:
477 
478  // Error
479  std::ostringstream error_stream;
480  error_stream << "Wrong imacro " << imacro << std::endl;
481  throw OomphLibError(
482  error_stream.str(),
483  OOMPH_CURRENT_FUNCTION,
484  OOMPH_EXCEPTION_LOCATION);
485  }
486 
487 }
488 
489 
490 
491 //=======================================================================
492 /// \short Boundary of central box macro element in layer i_layer
493 /// zeta \f$ \in [-1,1]^2 \f$
494 //=======================================================================
495 void QuarterTubeDomain::r_centr_L(const unsigned& t,
496  const Vector<double>& zeta,
497  const unsigned& i_layer,
498  Vector<double>& f)
499 {
500 
501  // Wall coordinates along top edge of wall
502  Vector<double> x(2);
503  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
504  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
505  x[1]=Xi_hi[1];
506 
507  // Get position vector to upper edge of wall
508  Vector<double> r_top(3);
509  Wall_pt->position(t,x,r_top);
510 
511  // Scale it down to half the height
512  f[0]=r_top[0];
513  f[1]=r_top[1]*0.25*(1.0+zeta[0]);
514  // Warp it:
515  double rho=0.0; //0.25*(1.0+zeta[0]);
516  f[2]=x[0]+rho*(r_top[2]-x[0]);
517 
518  //f[2]=r_top[2];
519 
520 
521 }
522 
523 
524 
525  //=======================================================================
526  /// \short Boundary of central box macro element in layer i_layer
527  /// zeta \f$ \in [-1,1]^2 \f$
528  //=======================================================================
529  void QuarterTubeDomain::r_centr_R(const unsigned& t,
530  const Vector<double>& zeta,
531  const unsigned& i_layer,
532  Vector<double>& f)
533 {
534  // Note the repetition in the calculation, there is some scope
535  // for optimisation
536 
537  // Wall coordinates along top edge of wall
538  Vector<double> x(2);
539  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
540  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
541  x[1]=Xi_hi[1];
542 
543  // Get position vector to upper edge of wall
544  Vector<double> r_top(3);
545  Wall_pt->position(t,x,r_top);
546 
547  // Wall coordinates along bottom edge of wall
548  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
549  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
550  x[1]=Xi_lo[1];
551 
552  // Get position vector to bottom edge of wall
553  Vector<double> r_bottom(3);
554  Wall_pt->position(t,x,r_bottom);
555 
556  // Scale it down to half the height, halfway across width
557  f[0]=0.5*r_bottom[0];
558  f[1]=r_top[1]*0.25*(1.0+zeta[0]);
559 
560  // Warp it:
561  double rho=0.0; //0.25*(1.0+zeta[0]);
562  f[2]=x[0]+rho*(r_top[2]-x[0]);
563  //f[2]=r_top[2];
564 
565 
566 
567 
568 }
569 
570 
571 
572 
573  //=======================================================================
574  /// \short Boundary of central box macro element in layer i_layer
575  /// zeta \f$ \in [-1,1]^2 \f$
576  //=======================================================================
577  void QuarterTubeDomain::r_centr_D(const unsigned& t,
578  const Vector<double>& zeta,
579  const unsigned& i_layer,
580  Vector<double>& f)
581 {
582 
583  // Wall coordinates along bottom edge of wall
584  Vector<double> x(2);
585  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
586  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
587  x[1]=Xi_lo[1];
588 
589  // Get position vector to bottom edge of wall
590  Vector<double> r_bottom(3);
591  Wall_pt->position(t,x,r_bottom);
592 
593  // Scale it down to half the width
594  f[0]=r_bottom[0]*0.25*(1.0+zeta[0]);
595  f[1]=r_bottom[1];
596 
597  // Warp it:
598  double rho=0.0; //0.25*(1.0+zeta[0]);
599  f[2]=x[0]+rho*(r_bottom[2]-x[0]);
600  //f[2]=r_bottom[2];
601 
602 
603 }
604 
605 
606  //=======================================================================
607  /// \short Boundary of central box macro element in layer i_layer
608  /// zeta \f$ \in [-1,1]^2 \f$
609  //=======================================================================
610  void QuarterTubeDomain::r_centr_U(const unsigned& t,
611  const Vector<double>& zeta,
612  const unsigned& i_layer,
613  Vector<double>& f)
614 {
615  // Wall coordinates along top edge of wall
616  Vector<double> x(2);
617  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
618  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
619  x[1]=Xi_hi[1];
620 
621  // Get position vector to upper edge of wall
622  Vector<double> r_top(3);
623  Wall_pt->position(t,x,r_top);
624 
625  // Wall coordinates along bottom edge of wall
626  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
627  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
628  x[1]=Xi_lo[1];
629 
630  // Get position vector to bottom edge of wall
631  Vector<double> r_bottom(3);
632  Wall_pt->position(t,x,r_bottom);
633 
634  // Scale it down to half the width
635  f[0]=r_bottom[0]*0.25*(1.0+zeta[0]);
636  f[1]=0.5*r_top[1];
637 
638  // Warp it:
639  double rho=0.0; //0.25*(1.0+zeta[0]);
640  f[2]=x[0]+rho*(r_bottom[2]-x[0]);
641  //f[2]=r_bottom[2];
642 
643 
644 }
645 
646 
647 
648  //=======================================================================
649  /// \short Boundary of central box macro element in layer i_layer
650  /// zeta \f$ \in [-1,1]^2 \f$
651  //=======================================================================
652  void QuarterTubeDomain::r_centr_B(const unsigned& t,
653  const Vector<double>& zeta,
654  const unsigned& i_layer,
655  Vector<double>& f)
656 {
657  // Wall coordinates along bottom edge of wall
658  Vector<double> x(2);
659  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
660  axial_spacing_fct(double(i_layer)/double(Nlayer));
661  x[1]=Xi_lo[1];
662 
663  // Get position vector to bottom edge of wall
664  Vector<double> r_bottom(3);
665  Wall_pt->position(t,x,r_bottom);
666 
667 
668  // Wall coordinates along top edge of wall
669  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
670  axial_spacing_fct(double(i_layer)/double(Nlayer));
671  x[1]=Xi_hi[1];
672 
673  // Get position vector to top edge of wall
674  Vector<double> r_top(3);
675  Wall_pt->position(t,x,r_top);
676 
677  // Map it
678  f[0]=r_bottom[0]*0.25*(1.0+zeta[0]);
679  f[1]=r_top[1]*0.25*(1.0+zeta[1]);
680 
681  // Warp it:
682  double rho=0.0; // 0.25*(1.0+zeta[1]);
683  f[2]=x[0]+rho*(r_top[2]-x[0]);
684  //f[2]=r_top[2];
685 
686 
687 }
688 
689 
690 
691 
692  //=======================================================================
693  /// \short Boundary of central box macro element in layer i_layer
694  /// zeta \f$ \in [-1,1]^2 \f$
695  //=======================================================================
696  void QuarterTubeDomain::r_centr_F(const unsigned& t,
697  const Vector<double>& zeta,
698  const unsigned& i_layer,
699  Vector<double>& f)
700 {
701  // Wall coordinates along bottom edge of wall
702  Vector<double> x(2);
703  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
704  axial_spacing_fct(double(1+i_layer)/double(Nlayer));
705  x[1]=Xi_lo[1];
706 
707  // Get position vector to bottom edge of wall
708  Vector<double> r_bottom(3);
709  Wall_pt->position(t,x,r_bottom);
710 
711 
712  // Wall coordinates along top edge of wall
713  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
714  axial_spacing_fct(double(1+i_layer)/double(Nlayer));
715  x[1]=Xi_hi[1];
716 
717  // Get position vector to top edge of wall
718  Vector<double> r_top(3);
719  Wall_pt->position(t,x,r_top);
720 
721  // Map it
722  f[0]=r_bottom[0]*0.25*(1.0+zeta[0]);
723  f[1]=r_top[1]*0.25*(1.0+zeta[1]);
724 
725  // Warp it:
726  double rho=0.0; // 0.25*(1.0+zeta[1]);
727  f[2]=x[0]+rho*(r_top[2]-x[0]);
728  //f[2]=r_top[2];
729 
730 
731 }
732 
733 
734 
735 //#####################################################################
736 
737 
738 
739  //=======================================================================
740  /// \short Boundary of bottom right box macro element in layer i_layer
741  /// zeta \f$ \in [-1,1]^2 \f$
742  //=======================================================================
743  void QuarterTubeDomain::r_bot_right_L(const unsigned& t,
744  const Vector<double>& zeta,
745  const unsigned& i_layer,
746  Vector<double>& f)
747 {
748  r_centr_R(t,zeta,i_layer,f);
749 }
750 
751 
752 
753  //=======================================================================
754  /// \short Boundary of bottom right box macro element in layer i_layer
755  /// zeta \f$ \in [-1,1]^2 \f$
756  //=======================================================================
757  void QuarterTubeDomain::r_bot_right_R(const unsigned& t,
758  const Vector<double>& zeta,
759  const unsigned& i_layer,
760  Vector<double>& f)
761 {
762  // Wall coordinates
763  Vector<double> x(2);
764  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
765  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
766  x[1]=Xi_lo[1]+Fract_mid*(Xi_hi[1]-Xi_lo[1])*0.5*(1.0+zeta[0]);
767 
768  // Get position vector on wall
769  Wall_pt->position(t,x,f);
770 
771 
772 }
773 
774 
775 
776 
777  //=======================================================================
778  /// \short Boundary of bottom right box macro element in layer i_layer
779  /// zeta \f$ \in [-1,1]^2 \f$
780  //=======================================================================
781  void QuarterTubeDomain::r_bot_right_D(const unsigned& t,
782  const Vector<double>& zeta,
783  const unsigned& i_layer,
784  Vector<double>& f)
785 {
786 
787  // Wall coordinates along bottom edge of wall
788  Vector<double> x(2);
789  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
790  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
791  x[1]=Xi_lo[1];
792 
793  // Get position vector to bottom edge of wall
794  Vector<double> r_bottom(3);
795  Wall_pt->position(t,x,r_bottom);
796 
797  // Scale it down to half the width
798  f[0]=0.5*r_bottom[0]*(1.0+s_squashed(0.5*(1.0+zeta[0])));
799  f[1]=r_bottom[1];
800 
801  // Warp it:
802  double rho=s_squashed(0.5*(1.0+zeta[0]));
803  f[2]=x[0]+rho*(r_bottom[2]-x[0]);
804  //f[2]=r_bottom[2];
805 
806 
807 }
808 
809 
810  //=======================================================================
811  /// \short Boundary of bottom right box macro element in layer i_layer
812  /// zeta \f$ \in [-1,1]^2 \f$
813  //=======================================================================
814  void QuarterTubeDomain::r_bot_right_U(const unsigned& t,
815  const Vector<double>& zeta,
816  const unsigned& i_layer,
817  Vector<double>& f)
818 {
819 
820  // Wall coordinates of dividing line
821  Vector<double> x(2);
822  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
823  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
824  x[1]=Xi_lo[1]+Fract_mid*(Xi_hi[1]-Xi_lo[1]);
825 
826  // Get position vector on dividing line
827  Vector<double> r_div(3);
828  Wall_pt->position(t,x,r_div);
829 
830 
831  // Position vector to corner of central box
832  Vector<double> zeta_central(2);
833  Vector<double> r_central(3);
834  zeta_central[0]=1.0;
835  zeta_central[1]=zeta[1];
836  r_centr_R(t,zeta_central,i_layer,r_central);
837 
838 
839 
840  // Straight line across
841  f[0]=r_central[0]+(r_div[0]-r_central[0])*s_squashed(0.5*(1.0+zeta[0]));
842  f[1]=r_central[1]+(r_div[1]-r_central[1])*s_squashed(0.5*(1.0+zeta[0]));
843  f[2]=r_central[2]+(r_div[2]-r_central[2])*s_squashed(0.5*(1.0+zeta[0]));
844 
845 }
846 
847 
848 
849  //=======================================================================
850  /// \short Boundary of bottom right box macro element in layer i_layer
851  /// zeta \f$ \in [-1,1]^2 \f$
852  //=======================================================================
853  void QuarterTubeDomain::r_bot_right_B(const unsigned& t,
854  const Vector<double>& zeta,
855  const unsigned& i_layer,
856  Vector<double>& f)
857 {
858  // Wall coordinates
859  Vector<double> x(2);
860  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
861  axial_spacing_fct(double(i_layer)/double(Nlayer));
862  x[1]=Xi_lo[1]+Fract_mid*(Xi_hi[1]-Xi_lo[1])*0.5*(1.0+zeta[1]);
863 
864  // Get position vector to wall
865  Vector<double> r_wall(3);
866  Wall_pt->position(t,x,r_wall);
867 
868  // Get position vector on central box
869  Vector<double> zeta_central(2);
870  Vector<double> r_central(3);
871  zeta_central[0]=zeta[1];
872  zeta_central[1]=-1.0;
873  r_centr_R(t,zeta_central,i_layer,r_central);
874 
875 
876  // Straight line across
877  f[0]=r_central[0]+(r_wall[0]-r_central[0])*s_squashed(0.5*(1.0+zeta[0]));
878  f[1]=r_central[1]+(r_wall[1]-r_central[1])*s_squashed(0.5*(1.0+zeta[0]));
879  f[2]=r_central[2]+(r_wall[2]-r_central[2])*s_squashed(0.5*(1.0+zeta[0]));
880 
881 
882 
883 
884 
885 }
886 
887 
888 
889 
890  //=======================================================================
891  /// \short Boundary of bottom right box macro element in layer i_layer
892  /// zeta \f$ \in [-1,1]^2 \f$
893  //=======================================================================
894  void QuarterTubeDomain::r_bot_right_F(const unsigned& t,
895  const Vector<double>& zeta,
896  const unsigned& i_layer,
897  Vector<double>& f)
898 {
899 
900  // Wall coordinates
901  Vector<double> x(2);
902  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
903  axial_spacing_fct(double(1+i_layer)/double(Nlayer));
904  x[1]=Xi_lo[1]+Fract_mid*(Xi_hi[1]-Xi_lo[1])*0.5*(1.0+zeta[1]);
905 
906  // Get position vector to wall
907  Vector<double> r_wall(3);
908  Wall_pt->position(t,x,r_wall);
909 
910  // Get position vector on central box
911  Vector<double> zeta_central(2);
912  Vector<double> r_central(3);
913  zeta_central[0]=zeta[1];
914  zeta_central[1]=1.0;
915  r_centr_R(t,zeta_central,i_layer,r_central);
916 
917 
918  // Straight line across
919  f[0]=r_central[0]+(r_wall[0]-r_central[0])*s_squashed(0.5*(1.0+zeta[0]));
920  f[1]=r_central[1]+(r_wall[1]-r_central[1])*s_squashed(0.5*(1.0+zeta[0]));
921  f[2]=r_central[2]+(r_wall[2]-r_central[2])*s_squashed(0.5*(1.0+zeta[0]));
922 
923 
924 }
925 
926 
927 
928 //#####################################################################
929 
930 
931 
932  //=======================================================================
933  /// \short Boundary of top left box macro element in layer i_layer
934  /// zeta \f$ \in [-1,1]^2 \f$
935  //=======================================================================
936  void QuarterTubeDomain::r_top_left_L(const unsigned& t,
937  const Vector<double>& zeta,
938  const unsigned& i_layer,
939  Vector<double>& f)
940 {
941 
942  // Wall coordinates along top edge of wall
943  Vector<double> x(2);
944  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
945  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
946  x[1]=Xi_hi[1];
947 
948  // Get position vector to upper edge of wall
949  Vector<double> r_top(3);
950  Wall_pt->position(t,x,r_top);
951 
952  // Scale it down to half the height
953  f[0]=r_top[0];
954  f[1]=0.5*r_top[1]*(1.0+s_squashed(0.5*(1.0+zeta[0])));
955 
956  // Warp it:
957  double rho=s_squashed(0.5*(1.0+zeta[0]));
958  f[2]=x[0]+rho*(r_top[2]-x[0]);
959  //f[2]=r_top[2];
960 
961 }
962 
963 
964 
965  //=======================================================================
966  /// \short Boundary of top left box macro element in layer i_layer
967  /// zeta \f$ \in [-1,1]^2 \f$
968  //=======================================================================
969  void QuarterTubeDomain::r_top_left_R(const unsigned& t,
970  const Vector<double>& zeta,
971  const unsigned& i_layer,
972  Vector<double>& f)
973 {
974  // Swap coordinates
975  Vector<double> zeta_br(2);
976  zeta_br[0]=zeta[0];
977  zeta_br[1]=zeta[1];
978  r_bot_right_U(t,zeta_br,i_layer,f);
979 }
980 
981 
982 
983 
984  //=======================================================================
985  /// \short Boundary of top left box macro element in layer i_layer
986  /// zeta \f$ \in [-1,1]^2 \f$
987  //=======================================================================
988  void QuarterTubeDomain::r_top_left_D(const unsigned& t,
989  const Vector<double>& zeta,
990  const unsigned& i_layer,
991  Vector<double>& f)
992 {
993  r_centr_U(t,zeta,i_layer,f);
994 }
995 
996 
997  //=======================================================================
998  /// \short Boundary of top left box macro element in layer i_layer
999  /// zeta \f$ \in [-1,1]^2 \f$
1000  //=======================================================================
1001  void QuarterTubeDomain::r_top_left_U(const unsigned& t,
1002  const Vector<double>& zeta,
1003  const unsigned& i_layer,
1004  Vector<double>& f)
1005 {
1006 
1007 
1008  // Wall coordinates
1009  Vector<double> x(2);
1010  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
1011  axial_spacing_fct((0.5*(1.0+zeta[1])+double(i_layer))/double(Nlayer));
1012  x[1]=Xi_hi[1]+(Xi_lo[1]-Xi_hi[1])*(1-Fract_mid)*0.5*(1.0+zeta[0]);
1013 
1014  // Get position vector on wall
1015  Wall_pt->position(t,x,f);
1016 
1017 }
1018 
1019 
1020 
1021  //=======================================================================
1022  /// \short Boundary of top left box macro element in layer i_layer
1023  /// zeta \f$ \in [-1,1]^2 \f$
1024  //=======================================================================
1025  void QuarterTubeDomain::r_top_left_B(const unsigned& t,
1026  const Vector<double>& zeta,
1027  const unsigned& i_layer,
1028  Vector<double>& f)
1029 {
1030 
1031  // Wall coordinates
1032  Vector<double> x(2);
1033  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
1034  axial_spacing_fct(double(i_layer)/double(Nlayer));
1035  x[1]=Xi_hi[1]+(Xi_lo[1]-Xi_hi[1])*(1.0-Fract_mid)*0.5*(1.0+zeta[0]);
1036 
1037 
1038  // Get position vector to wall
1039  Vector<double> r_wall(3);
1040  Wall_pt->position(t,x,r_wall);
1041 
1042 
1043  // Get position vector on central box
1044  Vector<double> zeta_central(2);
1045  Vector<double> r_central(3);
1046  zeta_central[0]=zeta[0];
1047  zeta_central[1]=-1.0;
1048  r_centr_U(t,zeta_central,i_layer,r_central);
1049 
1050  // Straight line across
1051  f[0]=r_central[0]+(r_wall[0]-r_central[0])*s_squashed(0.5*(1.0+zeta[1]));
1052  f[1]=r_central[1]+(r_wall[1]-r_central[1])*s_squashed(0.5*(1.0+zeta[1]));
1053  f[2]=r_central[2]+(r_wall[2]-r_central[2])*s_squashed(0.5*(1.0+zeta[1]));
1054 
1055 
1056 }
1057 
1058 
1059 
1060 
1061  //=======================================================================
1062  /// \short Boundary of top left box macro element in layer i_layer
1063  /// zeta \f$ \in [-1,1]^2 \f$
1064  //=======================================================================
1065  void QuarterTubeDomain::r_top_left_F(const unsigned& t,
1066  const Vector<double>& zeta,
1067  const unsigned& i_layer,
1068  Vector<double>& f)
1069 {
1070 
1071  // Wall coordinates
1072  Vector<double> x(2);
1073  x[0]=Xi_lo[0]+(Xi_hi[0]-Xi_lo[0])*
1074  axial_spacing_fct(double(1+i_layer)/double(Nlayer));
1075  x[1]=Xi_hi[1]+(Xi_lo[1]-Xi_hi[1])*(1.0-Fract_mid)*0.5*(1.0+zeta[0]);
1076 
1077 
1078  // Get position vector to wall
1079  Vector<double> r_wall(3);
1080  Wall_pt->position(t,x,r_wall);
1081 
1082 
1083  // Get position vector on central box
1084  Vector<double> zeta_central(2);
1085  Vector<double> r_central(3);
1086  zeta_central[0]=zeta[0];
1087  zeta_central[1]=1.0;
1088  r_centr_U(t,zeta_central,i_layer,r_central);
1089 
1090  // Straight line across
1091  f[0]=r_central[0]+(r_wall[0]-r_central[0])*s_squashed(0.5*(1.0+zeta[1]));
1092  f[1]=r_central[1]+(r_wall[1]-r_central[1])*s_squashed(0.5*(1.0+zeta[1]));
1093  f[2]=r_central[2]+(r_wall[2]-r_central[2])*s_squashed(0.5*(1.0+zeta[1]));
1094 
1095 
1096 }
1097 
1098 
1099 
1100 
1101 }
1102 
1103 #endif
AxialSpacingFctPt & axial_spacing_fct_pt()
Function pointer for function that implements axial spacing of macro elements.
void r_bot_right_L(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of bottom right box macro element in layer i_layer zeta .
double axial_spacing_fct(const double &xi)
Function that implements axial spacing of macro elements.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
void r_top_left_R(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of top left box macro element in layer i_layer zeta .
cstr elem_len * i
Definition: cfortran.h:607
static double default_axial_spacing_fct(const double &xi)
Default for function that implements axial spacing of macro elements.
void r_centr_L(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of central box macro element in layer i_layer zeta .
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition: domain.h:237
void r_top_left_B(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of top left box macro element in layer i_layer zeta .
char t
Definition: cfortran.h:572
unsigned Nlayer
Number of layers.
void r_top_left_F(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of top left box macro element in layer i_layer zeta .
BLSquashFctPt BL_squash_fct_pt
Function pointer for function that squashes the outer two macro elements towards the wall by mapping ...
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
static double default_BL_squash_fct(const double &s)
Default for function that squashes the outer two macro elements towards the wall by mapping the input...
QuarterTubeDomain(GeomObject *boundary_geom_object_pt, const Vector< double > &xi_lo, const double &fract_mid, const Vector< double > &xi_hi, const unsigned &nlayer)
Constructor: Pass boundary object and start and end coordinates and fraction along boundary object wh...
double(* BLSquashFctPt)(const double &s)
Typedef for function pointer for function that squashes the outer two macro elements towards the wall...
void r_bot_right_B(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of bottom right box macro element in layer i_layer zeta .
void r_centr_R(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of central box macro element in layer i_layer zeta .
GeomObject * Wall_pt
Pointer to geometric object that represents the curved wall.
BLSquashFctPt & bl_squash_fct_pt()
Function pointer for function that squashes the outer two macro elements towards the wall by mapping ...
Vector< double > Xi_lo
Lower limit for the coordinates along the wall.
void r_centr_D(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of central box macro element in layer i_layer zeta .
static char t char * s
Definition: cfortran.h:572
double(* AxialSpacingFctPt)(const double &xi)
Typedef for function pointer for function that implements axial spacing of macro elements.
void r_centr_B(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of central box macro element in layer i_layer zeta .
AxialSpacingFctPt Axial_spacing_fct_pt
Function pointer for function that implements axial spacing of macro elements.
void r_top_left_D(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of top left box macro element in layer i_layer zeta .
void r_bot_right_R(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of bottom right box macro element in layer i_layer zeta .
Quarter tube as domain. Domain is bounded by curved boundary which is represented by a GeomObject...
double s_squashed(const double &s)
Function that squashes the outer two macro elements towards the wall by mapping the input value of th...
Vector< double > Xi_hi
Upper limit for the coordinates along the wall.
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void r_bot_right_D(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of bottom right box macro element in layer i_layer zeta .
void macro_element_boundary(const unsigned &t, const unsigned &i_macro, const unsigned &i_direct, const Vector< double > &s, Vector< double > &f)
Vector representation of the i_macro-th macro element boundary i_direct (L/R/D/U/B/F) at time level t...
void r_top_left_L(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of top left box macro element in layer i_layer zeta .
QuarterTubeDomain(const QuarterTubeDomain &)
Broken copy constructor.
void r_centr_F(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of central box macro element in layer i_layer zeta .
void operator=(const QuarterTubeDomain &)
Broken assignment operator.
double Fract_mid
Fraction along wall where outer ring is to be divided.
void r_bot_right_U(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of bottom right box macro element in layer i_layer zeta .
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition: domain.h:71
~QuarterTubeDomain()
Destructor: Kill macro elements.
void r_bot_right_F(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of bottom right box macro element in layer i_layer zeta .
void r_top_left_U(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of top left box macro element in layer i_layer zeta .
void r_centr_U(const unsigned &t, const Vector< double > &zeta, const unsigned &i_layer, Vector< double > &f)
Boundary of central box macro element in layer i_layer zeta .