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-2023 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26// Include guards
27#ifndef OOMPH_TUBE_DOMAIN_HEADER
28#define OOMPH_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 /// Tube as a domain. The entire domain must be defined by
39 /// a GeomObject with the following convention: zeta[0] is the coordinate
40 /// along the centreline, zeta[1] is the theta coordinate around the tube
41 /// wall and zeta[2] is the radial coordinate.
42 /// The outer boundary must lie at zeta[2] = 1.
43 ///
44 /// The domain is
45 /// parametrised by five macro elements (a central box surrounded by
46 /// four curved elements) in each of the nlayer slices. The labelling
47 /// of the macro elements is as follows with the zeta[0] coordinate
48 /// coming out of the page.
49 ///
50 /// ----------------------------
51 /// |\ /|
52 /// | \ Macro / |
53 /// | 3 Element 3 2 |
54 /// | \ / |
55 /// | \----------------/ |
56 /// | | | |
57 /// | 4 | Macro | |
58 /// | | Element 0 | 2 |
59 /// | | | |
60 /// | ----------------- |
61 /// | / \ |
62 /// | 0 Macro 1 |
63 /// | / Element 1 \ |
64 /// | / \|
65 /// |/-------------------------|
66 ///
67 ///
68 //=================================================================
69 class TubeDomain : public Domain
70 {
71 public:
72 /// Constructor: Pass geometric object; start and end limit of the
73 /// centreline coordinate; the theta locations marking the division between
74 /// the elements of the outer ring, labelled from the lower left to the
75 /// upper left in order, theta should be in the range \f$-\pi\f$ to
76 /// \f$\pi\f$; the corresponding fractions of the
77 /// radius at which the central box is to be placed; and the number of
78 /// layers in the domain
79 TubeDomain(GeomObject* volume_geom_object_pt,
80 const Vector<double>& centreline_limits,
81 const Vector<double>& theta_positions,
82 const Vector<double>& radius_box,
83 const unsigned& nlayer)
84 : Centreline_limits(centreline_limits),
85 Theta_positions(theta_positions),
86 Radius_box(radius_box),
87 Nlayer(nlayer),
88 Volume_pt(volume_geom_object_pt)
89 {
90 // There are five macro elements
91 const unsigned n_macro = 5 * nlayer;
92 Macro_element_pt.resize(n_macro);
93
94 // Create the macro elements
95 for (unsigned i = 0; i < n_macro; i++)
96 {
97 Macro_element_pt[i] = new QMacroElement<3>(this, i);
98 }
99 }
100
101 /// Broken copy constructor
102 TubeDomain(const TubeDomain&) = delete;
103
104 /// Broken assignment operator
105 void operator=(const TubeDomain&) = delete;
106
107 /// Destructor: Empty; cleanup done in base class
108 ~TubeDomain() {}
109
110 /// Vector representation of the i_macro-th macro element
111 /// boundary i_direct (L/R/D/U/B/F) at time level t
112 /// (t=0: present; t>0: previous):
113 /// f(s).
114 void macro_element_boundary(const unsigned& t,
115 const unsigned& i_macro,
116 const unsigned& i_direct,
117 const Vector<double>& s,
118 Vector<double>& f);
119
120 private:
121 /// Storage for the limits of the centreline coordinate
122 Vector<double> Centreline_limits;
123
124 /// Storage for the dividing lines on the boundary
125 /// starting from the lower left and proceeding anticlockwise to
126 /// the upper left
127 Vector<double> Theta_positions;
128
129 /// Storage for the fraction of the radius at which the central box
130 /// should be located corresponding to the chosen values of theta.
131 Vector<double> Radius_box;
132
133 /// Number of axial layers
134 unsigned Nlayer;
135
136 /// Pointer to geometric object that represents the domain
137 GeomObject* Volume_pt;
138
139 /// A very little linear interpolation helper.
140 /// Interpolate from the low
141 /// point to the high point using the coordinate s which is
142 /// assumed to run from -1 to 1.
143 void lin_interpolate(const Vector<double>& low,
144 const Vector<double>& high,
145 const double& s,
146 Vector<double>& f)
147 {
148 // Loop over all coordinates
149 for (unsigned i = 0; i < 3; i++)
150 {
151 f[i] = low[i] + (high[i] - low[i]) * 0.5 * (s + 1.0);
152 }
153 }
154 };
155
156
157 /// //////////////////////////////////////////////////////////////////////
158 /// //////////////////////////////////////////////////////////////////////
159 /// //////////////////////////////////////////////////////////////////////
160
161
162 //=================================================================
163 /// Vector representation of the imacro-th macro element
164 /// boundary idirect (L/R/D/U/B/F) at time level t
165 /// (t=0: present; t>0: previous): f(s)
166 //=================================================================
167 void TubeDomain::macro_element_boundary(const unsigned& t,
168 const unsigned& imacro,
169 const unsigned& idirect,
170 const Vector<double>& s,
171 Vector<double>& f)
172 {
173 using namespace OcTreeNames;
174
175#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
176 // Warn about time argument being moved to the front
177 OomphLibWarning(
178 "Order of function arguments has changed between versions 0.8 and 0.85",
179 "TubeDomain::macro_element_boundary(...)",
180 OOMPH_EXCEPTION_LOCATION);
181#endif
182
183 // Get the number of the layer
184 unsigned ilayer = unsigned(imacro / 5);
185
186 // Get all required coordinates of the corners of the box
187 // at each end of the layer
188 Vector<Vector<Vector<double>>> Box(2);
189
190 // Get the centreline coordinates at the ends of the layer
191 Vector<double> zeta_centre(2);
192 // Storage for position in the volume
193 Vector<double> zeta(3);
194
195 // Loop over the layers
196 for (unsigned i = 0; i < 2; i++)
197 {
198 // Resize the storage
199 Box[i].resize(4);
200
201 // Get the centreline coordinate
202 zeta_centre[i] =
203 Centreline_limits[0] + (ilayer + i) *
204 (Centreline_limits[1] - Centreline_limits[0]) /
205 (double)(Nlayer);
206
207 // Now get the coordinates of each corner
208 zeta[0] = zeta_centre[i];
209
210 // Loop over the angles
211 for (unsigned j = 0; j < 4; j++)
212 {
213 // Set up the storage
214 Box[i][j].resize(3);
215
216 // Set the other values of zeta
217 zeta[1] = Theta_positions[j];
218 zeta[2] = Radius_box[j];
219
220 // Now get the values
221 Volume_pt->position(t, zeta, Box[i][j]);
222 }
223 }
224
225 // Local storage for positions on the boundaries
226 Vector<double> pos_1(3), pos_2(3);
227
228 const double pi = MathematicalConstants::Pi;
229
230 // Which macro element?
231 // --------------------
232 switch (imacro % 5)
233 {
234 // Macro element 0: Central box
235 case 0:
236
237 // Choose a direction
238 switch (idirect)
239 {
240 case L:
241 // Need to get the position from the domain
242 // Get the centreline position
243 zeta[0] = zeta_centre[0] +
244 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
245 // Get the lower corner
246 zeta[1] = Theta_positions[0];
247 zeta[2] = Radius_box[0];
248 Volume_pt->position(t, zeta, pos_1);
249
250 // Get the upper corner
251 zeta[1] = Theta_positions[3];
252 zeta[2] = Radius_box[3];
253 Volume_pt->position(t, zeta, pos_2);
254
255 // Now interpolate between the two corner positions
256 lin_interpolate(pos_1, pos_2, s[0], f);
257 break;
258
259 case R:
260 // Need to get the position from the domain
261 // Get the centreline position
262 zeta[0] = zeta_centre[0] +
263 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
264 // Get the lower corner
265 zeta[1] = Theta_positions[1];
266 zeta[2] = Radius_box[1];
267 Volume_pt->position(t, zeta, pos_1);
268
269 // Get the upper corner
270 zeta[1] = Theta_positions[2];
271 zeta[2] = Radius_box[2];
272 Volume_pt->position(t, zeta, pos_2);
273
274 // Now interpolate between the positions
275 lin_interpolate(pos_1, pos_2, s[0], f);
276 break;
277
278 case D:
279 // Need to get the position from the domain
280 // Get the centreline position
281 zeta[0] = zeta_centre[0] +
282 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
283 // Get the left corner
284 zeta[1] = Theta_positions[0];
285 zeta[2] = Radius_box[0];
286 Volume_pt->position(t, zeta, pos_1);
287
288 // Get the right corner
289 zeta[1] = Theta_positions[1];
290 zeta[2] = Radius_box[1];
291 Volume_pt->position(t, zeta, pos_2);
292 // Now interpolate between the positions
293 lin_interpolate(pos_1, pos_2, s[0], f);
294 break;
295
296 case U:
297 // Need to get the position from the domain
298 // Get the centreline position
299 zeta[0] = zeta_centre[0] +
300 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
301 // Get the left corner
302 zeta[1] = Theta_positions[3];
303 zeta[2] = Radius_box[3];
304 Volume_pt->position(t, zeta, pos_1);
305
306 // Get the right corner
307 zeta[1] = Theta_positions[2];
308 zeta[2] = Radius_box[2];
309 Volume_pt->position(t, zeta, pos_2);
310
311 // Now interpolate between the positions
312 lin_interpolate(pos_1, pos_2, s[0], f);
313 break;
314
315 case B:
316 // Linearly interpolate between lower left and lower right corners
317 lin_interpolate(Box[0][0], Box[0][1], s[0], pos_1);
318 // Linearly interpolate between upper left and upper right corners
319 lin_interpolate(Box[0][3], Box[0][2], s[0], pos_2);
320 // Finally, linearly interpolate between the upper and lower
321 // positions
322 lin_interpolate(pos_1, pos_2, s[1], f);
323 break;
324
325 case F:
326 // Linearly interpolate between lower left and lower right corners
327 lin_interpolate(Box[1][0], Box[1][1], s[0], pos_1);
328 // Linearly interpolate between upper left and upper right corners
329 lin_interpolate(Box[1][3], Box[1][2], s[0], pos_2);
330 // Finally, linearly interpolate between the upper and lower
331 // positions
332 lin_interpolate(pos_1, pos_2, s[1], f);
333 break;
334
335 default:
336 std::ostringstream error_stream;
337 error_stream << "idirect is " << idirect
338 << " not one of L, R, D, U, B, F" << std::endl;
339
340 throw OomphLibError(error_stream.str(),
341 OOMPH_CURRENT_FUNCTION,
342 OOMPH_EXCEPTION_LOCATION);
343 break;
344 }
345
346 break;
347
348
349 // Macro element 1: Bottom
350 case 1:
351
352 // Choose a direction
353 switch (idirect)
354 {
355 case L:
356 // Get the position on the wall
357 zeta[0] = zeta_centre[0] +
358 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
359 zeta[1] = Theta_positions[0];
360 zeta[2] = 1.0;
361 Volume_pt->position(t, zeta, pos_1);
362
363 // Get the position on the box
364 zeta[2] = Radius_box[0];
365 Volume_pt->position(t, zeta, pos_2);
366
367 // Now linearly interpolate between the two
368 lin_interpolate(pos_1, pos_2, s[0], f);
369 break;
370
371 case R:
372 // Get the position on the wall
373 zeta[0] = zeta_centre[0] +
374 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
375 zeta[1] = Theta_positions[1];
376 zeta[2] = 1.0;
377 Volume_pt->position(t, zeta, pos_1);
378
379 // Get the position from the box
380 zeta[2] = Radius_box[1];
381 Volume_pt->position(t, zeta, pos_2);
382
383 // Now linearly interpolate between the two
384 lin_interpolate(pos_1, pos_2, s[0], f);
385 break;
386
387 case D:
388 // This is entrirly on the wall
389 zeta[0] = zeta_centre[0] +
390 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
391 zeta[1] =
392 Theta_positions[0] +
393 (Theta_positions[1] - Theta_positions[0]) * 0.5 * (s[0] + 1.0);
394 zeta[2] = 1.0;
395 Volume_pt->position(t, zeta, f);
396 break;
397
398 case U:
399 // Need to get the position from the domain
400 // Get the centreline position
401 zeta[0] = zeta_centre[0] +
402 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
403 // Get the left corner
404 zeta[1] = Theta_positions[0];
405 zeta[2] = Radius_box[0];
406 Volume_pt->position(t, zeta, pos_1);
407
408 // Get the right corner
409 zeta[1] = Theta_positions[1];
410 zeta[2] = Radius_box[1];
411 Volume_pt->position(t, zeta, pos_2);
412 // Now interpolate between the positions
413 lin_interpolate(pos_1, pos_2, s[0], f);
414 break;
415
416 case B:
417 // Get the position on the wall
418 zeta[0] = zeta_centre[0];
419 zeta[1] =
420 Theta_positions[0] +
421 (Theta_positions[1] - Theta_positions[0]) * 0.5 * (s[0] + 1.0);
422 zeta[2] = 1.0;
423 Volume_pt->position(t, zeta, pos_1);
424
425 // Now linearly interpolate the position on the box
426 lin_interpolate(Box[0][0], Box[0][1], s[0], pos_2);
427
428 // Now linearly interpolate between the two
429 lin_interpolate(pos_1, pos_2, s[1], f);
430 break;
431
432 case F:
433 // Get the position on the wall
434 zeta[0] = zeta_centre[1];
435 zeta[1] =
436 Theta_positions[0] +
437 (Theta_positions[1] - Theta_positions[0]) * 0.5 * (s[0] + 1.0);
438 zeta[2] = 1.0;
439 Volume_pt->position(t, zeta, pos_1);
440
441 // Now linearly interpolate the position on the box
442 lin_interpolate(Box[1][0], Box[1][1], s[0], pos_2);
443
444 // Now linearly interpolate between the two
445 lin_interpolate(pos_1, pos_2, s[1], f);
446 break;
447
448
449 default:
450
451 std::ostringstream error_stream;
452 error_stream << "idirect is " << idirect
453 << " not one of L, R, D, U, B, F" << std::endl;
454
455 throw OomphLibError(error_stream.str(),
456 OOMPH_CURRENT_FUNCTION,
457 OOMPH_EXCEPTION_LOCATION);
458 break;
459 }
460
461
462 break;
463
464 // Macro element 2:Right
465 case 2:
466
467 // Which direction?
468 switch (idirect)
469 {
470 case L:
471 // Need to get the position from the domain
472 // Get the centreline position
473 zeta[0] = zeta_centre[0] +
474 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
475 // Get the lower corner
476 zeta[1] = Theta_positions[1];
477 zeta[2] = Radius_box[1];
478 Volume_pt->position(t, zeta, pos_1);
479
480 // Get the upper corner
481 zeta[1] = Theta_positions[2];
482 zeta[2] = Radius_box[2];
483 Volume_pt->position(t, zeta, pos_2);
484 // Now interpolate between the positions
485 lin_interpolate(pos_1, pos_2, s[0], f);
486 break;
487
488 case R:
489 // Entirely on the wall
490 zeta[0] = zeta_centre[0] +
491 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
492 zeta[1] =
493 Theta_positions[1] +
494 (Theta_positions[2] - Theta_positions[1]) * 0.5 * (s[0] + 1.0);
495 zeta[2] = 1.0;
496 Volume_pt->position(t, zeta, f);
497 break;
498
499 case U:
500 // Get the position on the wall
501 zeta[0] = zeta_centre[0] +
502 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
503 zeta[1] = Theta_positions[2];
504 zeta[2] = 1.0;
505 Volume_pt->position(t, zeta, pos_1);
506
507 // Get the position of the box
508 zeta[2] = Radius_box[2];
509 Volume_pt->position(t, zeta, pos_2);
510
511 // Now linearly interpolate between the two
512 lin_interpolate(pos_2, pos_1, s[0], f);
513 break;
514
515 case D:
516 // Get the position on the wall
517 zeta[0] = zeta_centre[0] +
518 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
519 zeta[1] = Theta_positions[1];
520 zeta[2] = 1.0;
521 Volume_pt->position(t, zeta, pos_1);
522
523 // Get the position of the box
524 zeta[2] = Radius_box[1];
525 Volume_pt->position(t, zeta, pos_2);
526
527 // Now linearly interpolate between the two
528 lin_interpolate(pos_2, pos_1, s[0], f);
529 break;
530
531 case B:
532 // Get the position on the wall
533 zeta[0] = zeta_centre[0];
534 zeta[1] =
535 Theta_positions[1] +
536 (Theta_positions[2] - Theta_positions[1]) * 0.5 * (s[1] + 1.0);
537 zeta[2] = 1.0;
538 Volume_pt->position(t, zeta, pos_1);
539
540 // Now linearly interpolate the position on the box
541 lin_interpolate(Box[0][1], Box[0][2], s[1], pos_2);
542
543 // Now linearly interpolate between the two
544 lin_interpolate(pos_2, pos_1, s[0], f);
545 break;
546
547 case F:
548 // Get the position on the wall
549 zeta[0] = zeta_centre[1];
550 zeta[1] =
551 Theta_positions[1] +
552 (Theta_positions[2] - Theta_positions[1]) * 0.5 * (s[1] + 1.0);
553 zeta[2] = 1.0;
554 Volume_pt->position(t, zeta, pos_1);
555
556 // Now linearly interpolate the position on the box
557 lin_interpolate(Box[1][1], Box[1][2], s[1], pos_2);
558
559 // Now linearly interpolate between the two
560 lin_interpolate(pos_2, pos_1, s[0], f);
561 break;
562
563
564 default:
565 std::ostringstream error_stream;
566 error_stream << "idirect is " << idirect
567 << " not one of L, R, D, U, B, F" << std::endl;
568
569 throw OomphLibError(error_stream.str(),
570 OOMPH_CURRENT_FUNCTION,
571 OOMPH_EXCEPTION_LOCATION);
572 }
573
574 break;
575
576 // Macro element 3: Top
577 case 3:
578
579 // Which direction?
580 switch (idirect)
581 {
582 case L:
583 // Get the position on the wall
584 zeta[0] = zeta_centre[0] +
585 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
586 zeta[1] = Theta_positions[3];
587 zeta[2] = 1.0;
588 Volume_pt->position(t, zeta, pos_1);
589
590 // Get the position on the box
591 zeta[2] = Radius_box[3];
592 Volume_pt->position(t, zeta, pos_2);
593
594
595 // Now linearly interpolate between the two
596 lin_interpolate(pos_2, pos_1, s[0], f);
597 break;
598
599 case R:
600 // Get the position on the wall
601 zeta[0] = zeta_centre[0] +
602 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
603 zeta[1] = Theta_positions[2];
604 zeta[2] = 1.0;
605 Volume_pt->position(t, zeta, pos_1);
606
607 // Get the position on the box
608 zeta[2] = Radius_box[2];
609 Volume_pt->position(t, zeta, pos_2);
610
611
612 // Now linearly interpolate between the two
613 lin_interpolate(pos_2, pos_1, s[0], f);
614 break;
615
616 case D:
617 // This is entirely on the box
618 // Need to get the position from the domain
619 // Get the centreline position
620 zeta[0] = zeta_centre[0] +
621 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
622 zeta[1] = Theta_positions[3];
623 zeta[2] = Radius_box[3];
624 // Get the lower corner
625 Volume_pt->position(t, zeta, pos_2);
626
627 // Get the upper corner
628 zeta[1] = Theta_positions[2];
629 zeta[2] = Radius_box[2];
630 // Get the upper corner
631 Volume_pt->position(t, zeta, pos_1);
632 // Now interpolate between the positions
633 lin_interpolate(pos_2, pos_1, s[0], f);
634 break;
635
636 case U:
637 // This is entirely on the wall
638 zeta[0] = zeta_centre[0] +
639 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
640 zeta[1] =
641 Theta_positions[3] +
642 (Theta_positions[2] - Theta_positions[3]) * 0.5 * (s[0] + 1.0);
643 zeta[2] = 1.0;
644 Volume_pt->position(t, zeta, f);
645 break;
646
647 case B:
648 // Get the position on the wall
649 zeta[0] = zeta_centre[0];
650 zeta[1] =
651 Theta_positions[3] +
652 (Theta_positions[2] - Theta_positions[3]) * 0.5 * (s[0] + 1.0);
653 zeta[2] = 1.0;
654 Volume_pt->position(t, zeta, pos_1);
655
656 // Now linearly interpolate the position on the box
657 lin_interpolate(Box[0][3], Box[0][2], s[0], pos_2);
658
659 // Now linearly interpolate between the two
660 lin_interpolate(pos_2, pos_1, s[1], f);
661 break;
662
663 case F:
664 // Get the position on the wall
665 zeta[0] = zeta_centre[1];
666 zeta[1] =
667 Theta_positions[3] +
668 (Theta_positions[2] - Theta_positions[3]) * 0.5 * (s[0] + 1.0);
669 zeta[2] = 1.0;
670 Volume_pt->position(t, zeta, pos_1);
671
672 // Now linearly interpolate the position on the box
673 lin_interpolate(Box[1][3], Box[1][2], s[0], pos_2);
674
675 // Now linearly interpolate between the two
676 lin_interpolate(pos_2, pos_1, s[1], f);
677 break;
678
679
680 default:
681 std::ostringstream error_stream;
682 error_stream << "idirect is " << idirect
683 << " not one of L, R, D, U, B, F" << std::endl;
684
685 throw OomphLibError(error_stream.str(),
686 OOMPH_CURRENT_FUNCTION,
687 OOMPH_EXCEPTION_LOCATION);
688 }
689
690 break;
691
692
693 // Macro element 4: Left
694 case 4:
695
696 // Which direction?
697 switch (idirect)
698 {
699 case L:
700 // Entirely on the wall, Need to be careful
701 // because this is the point at which theta passes through the
702 // branch cut
703 zeta[0] = zeta_centre[0] +
704 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
705 zeta[1] = Theta_positions[0] + 2.0 * pi +
706 (Theta_positions[3] - (Theta_positions[0] + 2.0 * pi)) *
707 0.5 * (s[0] + 1.0);
708 zeta[2] = 1.0;
709 Volume_pt->position(t, zeta, f);
710 break;
711
712 case R:
713 // Entirely on the box
714 // Need to get the position from the domain
715 // Get the centreline position
716 zeta[0] = zeta_centre[0] +
717 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
718 zeta[1] = Theta_positions[0];
719 zeta[2] = Radius_box[0];
720 // Get the lower corner
721 Volume_pt->position(t, zeta, pos_2);
722
723 // Get the upper corner
724 zeta[1] = Theta_positions[3];
725 zeta[2] = Radius_box[3];
726 // Get the upper corner
727 Volume_pt->position(t, zeta, pos_1);
728
729 lin_interpolate(pos_2, pos_1, s[0], f);
730 break;
731
732 case D:
733 // Get the position on the wall
734 zeta[0] = zeta_centre[0] +
735 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
736 zeta[1] = Theta_positions[0];
737 zeta[2] = 1.0;
738 Volume_pt->position(t, zeta, pos_1);
739
740
741 // Get the position on the box
742 zeta[2] = Radius_box[0];
743 Volume_pt->position(t, zeta, pos_2);
744
745
746 // Now linearly interpolate between the two
747 lin_interpolate(pos_1, pos_2, s[0], f);
748 break;
749
750 case U:
751 // Get the position on the wall
752 zeta[0] = zeta_centre[0] +
753 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
754 zeta[1] = Theta_positions[3];
755 zeta[2] = 1.0;
756 Volume_pt->position(t, zeta, pos_1);
757
758 // Get the position on the box
759 zeta[2] = Radius_box[3];
760 Volume_pt->position(t, zeta, pos_2);
761
762 // Now linearly interpolate between the two
763 lin_interpolate(pos_1, pos_2, s[0], f);
764 break;
765
766 case B:
767 // Get the position on the wall
768 // Again be careful of the branch cut
769 zeta[0] = zeta_centre[0];
770 zeta[1] = Theta_positions[0] + 2.0 * pi +
771 (Theta_positions[3] - (Theta_positions[0] + 2.0 * pi)) *
772 0.5 * (s[1] + 1.0);
773 zeta[2] = 1.0;
774 Volume_pt->position(t, zeta, pos_1);
775
776 // Now linearly interpolate the position on the box
777 lin_interpolate(Box[0][0], Box[0][3], s[1], pos_2);
778
779 // Now linearly interpolate between the two
780 lin_interpolate(pos_1, pos_2, s[0], f);
781 break;
782
783
784 case F:
785 // Get the position on the wall
786 // Again be careful of the branch cut
787 zeta[0] = zeta_centre[1];
788 zeta[1] = Theta_positions[0] + 2.0 * pi +
789 (Theta_positions[3] - (Theta_positions[0] + 2.0 * pi)) *
790 0.5 * (s[1] + 1.0);
791 zeta[2] = 1.0;
792 Volume_pt->position(t, zeta, pos_1);
793
794 // Now linearly interpolate the position on the box
795 lin_interpolate(Box[1][0], Box[1][3], s[1], pos_2);
796
797 // Now linearly interpolate between the two
798 lin_interpolate(pos_1, pos_2, s[0], f);
799 break;
800
801
802 default:
803 std::ostringstream error_stream;
804 error_stream << "idirect is " << idirect
805 << " not one of L, R, D, U, B, F" << std::endl;
806
807 throw OomphLibError(error_stream.str(),
808 OOMPH_CURRENT_FUNCTION,
809 OOMPH_EXCEPTION_LOCATION);
810 }
811 break;
812
813
814 default:
815 // Error
816 std::ostringstream error_stream;
817 error_stream << "Wrong imacro " << imacro << std::endl;
818 throw OomphLibError(
819 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
820 }
821 }
822
823} // namespace oomph
824
825#endif
oomph::RefineableTriangleMesh
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition triangle_mesh.template.h:2249
oomph::TubeDomain
Tube as a domain. The entire domain must be defined by a GeomObject with the following convention: ze...
Definition tube_domain.h:70
oomph::TubeDomain::TubeDomain
TubeDomain(GeomObject *volume_geom_object_pt, const Vector< double > &centreline_limits, const Vector< double > &theta_positions, const Vector< double > &radius_box, const unsigned &nlayer)
Constructor: Pass geometric object; start and end limit of the centreline coordinate; the theta locat...
Definition tube_domain.h:79
oomph::TubeDomain::macro_element_boundary
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...
Definition tube_domain.h:167
oomph::TubeDomain::TubeDomain
TubeDomain(const TubeDomain &)=delete
Broken copy constructor.
oomph::TubeDomain::Centreline_limits
Vector< double > Centreline_limits
Storage for the limits of the centreline coordinate.
Definition tube_domain.h:122
oomph::TubeDomain::operator=
void operator=(const TubeDomain &)=delete
Broken assignment operator.
oomph::TubeDomain::lin_interpolate
void lin_interpolate(const Vector< double > &low, const Vector< double > &high, const double &s, Vector< double > &f)
A very little linear interpolation helper. Interpolate from the low point to the high point using the...
Definition tube_domain.h:143
oomph::TubeDomain::Volume_pt
GeomObject * Volume_pt
Pointer to geometric object that represents the domain.
Definition tube_domain.h:137
oomph::TubeDomain::Nlayer
unsigned Nlayer
Number of axial layers.
Definition tube_domain.h:134
oomph::TubeDomain::Theta_positions
Vector< double > Theta_positions
Storage for the dividing lines on the boundary starting from the lower left and proceeding anticlockw...
Definition tube_domain.h:127
oomph::TubeDomain::Radius_box
Vector< double > Radius_box
Storage for the fraction of the radius at which the central box should be located corresponding to th...
Definition tube_domain.h:131
oomph::TubeDomain::~TubeDomain
~TubeDomain()
Destructor: Empty; cleanup done in base class.
Definition tube_domain.h:108
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition annular_domain.h:35