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