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-2025 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 /// Vector representation of the imacro-th macro element
163 /// boundary idirect (L/R/D/U/B/F) at time level t
164 /// (t=0: present; t>0: previous): f(s)
165 //=================================================================
166 void TubeDomain::macro_element_boundary(const unsigned& t,
167 const unsigned& imacro,
168 const unsigned& idirect,
169 const Vector<double>& s,
170 Vector<double>& f)
171 {
172 using namespace OcTreeNames;
173
174#ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
175 // Warn about time argument being moved to the front
176 OomphLibWarning(
177 "Order of function arguments has changed between versions 0.8 and 0.85",
178 "TubeDomain::macro_element_boundary(...)",
179 OOMPH_EXCEPTION_LOCATION);
180#endif
181
182 // Get the number of the layer
183 unsigned ilayer = unsigned(imacro / 5);
184
185 // Get all required coordinates of the corners of the box
186 // at each end of the layer
187 Vector<Vector<Vector<double>>> Box(2);
188
189 // Get the centreline coordinates at the ends of the layer
190 Vector<double> zeta_centre(2);
191 // Storage for position in the volume
192 Vector<double> zeta(3);
193
194 // Loop over the layers
195 for (unsigned i = 0; i < 2; i++)
196 {
197 // Resize the storage
198 Box[i].resize(4);
199
200 // Get the centreline coordinate
201 zeta_centre[i] =
202 Centreline_limits[0] + (ilayer + i) *
203 (Centreline_limits[1] - Centreline_limits[0]) /
204 (double)(Nlayer);
205
206 // Now get the coordinates of each corner
207 zeta[0] = zeta_centre[i];
208
209 // Loop over the angles
210 for (unsigned j = 0; j < 4; j++)
211 {
212 // Set up the storage
213 Box[i][j].resize(3);
214
215 // Set the other values of zeta
216 zeta[1] = Theta_positions[j];
217 zeta[2] = Radius_box[j];
218
219 // Now get the values
220 Volume_pt->position(t, zeta, Box[i][j]);
221 }
222 }
223
224 // Local storage for positions on the boundaries
225 Vector<double> pos_1(3), pos_2(3);
226
227 const double pi = MathematicalConstants::Pi;
228
229 // Which macro element?
230 // --------------------
231 switch (imacro % 5)
232 {
233 // Macro element 0: Central box
234 case 0:
235
236 // Choose a direction
237 switch (idirect)
238 {
239 case L:
240 // Need to get the position from the domain
241 // Get the centreline position
242 zeta[0] = zeta_centre[0] +
243 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
244 // Get the lower corner
245 zeta[1] = Theta_positions[0];
246 zeta[2] = Radius_box[0];
247 Volume_pt->position(t, zeta, pos_1);
248
249 // Get the upper corner
250 zeta[1] = Theta_positions[3];
251 zeta[2] = Radius_box[3];
252 Volume_pt->position(t, zeta, pos_2);
253
254 // Now interpolate between the two corner positions
255 lin_interpolate(pos_1, pos_2, s[0], f);
256 break;
257
258 case R:
259 // Need to get the position from the domain
260 // Get the centreline position
261 zeta[0] = zeta_centre[0] +
262 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
263 // Get the lower corner
264 zeta[1] = Theta_positions[1];
265 zeta[2] = Radius_box[1];
266 Volume_pt->position(t, zeta, pos_1);
267
268 // Get the upper corner
269 zeta[1] = Theta_positions[2];
270 zeta[2] = Radius_box[2];
271 Volume_pt->position(t, zeta, pos_2);
272
273 // Now interpolate between the positions
274 lin_interpolate(pos_1, pos_2, s[0], f);
275 break;
276
277 case D:
278 // Need to get the position from the domain
279 // Get the centreline position
280 zeta[0] = zeta_centre[0] +
281 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
282 // Get the left corner
283 zeta[1] = Theta_positions[0];
284 zeta[2] = Radius_box[0];
285 Volume_pt->position(t, zeta, pos_1);
286
287 // Get the right corner
288 zeta[1] = Theta_positions[1];
289 zeta[2] = Radius_box[1];
290 Volume_pt->position(t, zeta, pos_2);
291 // Now interpolate between the positions
292 lin_interpolate(pos_1, pos_2, s[0], f);
293 break;
294
295 case U:
296 // Need to get the position from the domain
297 // Get the centreline position
298 zeta[0] = zeta_centre[0] +
299 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
300 // Get the left corner
301 zeta[1] = Theta_positions[3];
302 zeta[2] = Radius_box[3];
303 Volume_pt->position(t, zeta, pos_1);
304
305 // Get the right corner
306 zeta[1] = Theta_positions[2];
307 zeta[2] = Radius_box[2];
308 Volume_pt->position(t, zeta, pos_2);
309
310 // Now interpolate between the positions
311 lin_interpolate(pos_1, pos_2, s[0], f);
312 break;
313
314 case B:
315 // Linearly interpolate between lower left and lower right corners
316 lin_interpolate(Box[0][0], Box[0][1], s[0], pos_1);
317 // Linearly interpolate between upper left and upper right corners
318 lin_interpolate(Box[0][3], Box[0][2], s[0], pos_2);
319 // Finally, linearly interpolate between the upper and lower
320 // positions
321 lin_interpolate(pos_1, pos_2, s[1], f);
322 break;
323
324 case F:
325 // Linearly interpolate between lower left and lower right corners
326 lin_interpolate(Box[1][0], Box[1][1], s[0], pos_1);
327 // Linearly interpolate between upper left and upper right corners
328 lin_interpolate(Box[1][3], Box[1][2], s[0], pos_2);
329 // Finally, linearly interpolate between the upper and lower
330 // positions
331 lin_interpolate(pos_1, pos_2, s[1], f);
332 break;
333
334 default:
335 std::ostringstream error_stream;
336 error_stream << "idirect is " << idirect
337 << " not one of L, R, D, U, B, F" << std::endl;
338
339 throw OomphLibError(error_stream.str(),
340 OOMPH_CURRENT_FUNCTION,
341 OOMPH_EXCEPTION_LOCATION);
342 break;
343 }
344
345 break;
346
347 // Macro element 1: Bottom
348 case 1:
349
350 // Choose a direction
351 switch (idirect)
352 {
353 case L:
354 // Get the position on the wall
355 zeta[0] = zeta_centre[0] +
356 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
357 zeta[1] = Theta_positions[0];
358 zeta[2] = 1.0;
359 Volume_pt->position(t, zeta, pos_1);
360
361 // Get the position on the box
362 zeta[2] = Radius_box[0];
363 Volume_pt->position(t, zeta, pos_2);
364
365 // Now linearly interpolate between the two
366 lin_interpolate(pos_1, pos_2, s[0], f);
367 break;
368
369 case R:
370 // Get the position on the wall
371 zeta[0] = zeta_centre[0] +
372 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
373 zeta[1] = Theta_positions[1];
374 zeta[2] = 1.0;
375 Volume_pt->position(t, zeta, pos_1);
376
377 // Get the position from the box
378 zeta[2] = Radius_box[1];
379 Volume_pt->position(t, zeta, pos_2);
380
381 // Now linearly interpolate between the two
382 lin_interpolate(pos_1, pos_2, s[0], f);
383 break;
384
385 case D:
386 // This is entrirly on the wall
387 zeta[0] = zeta_centre[0] +
388 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
389 zeta[1] =
390 Theta_positions[0] +
391 (Theta_positions[1] - Theta_positions[0]) * 0.5 * (s[0] + 1.0);
392 zeta[2] = 1.0;
393 Volume_pt->position(t, zeta, f);
394 break;
395
396 case U:
397 // Need to get the position from the domain
398 // Get the centreline position
399 zeta[0] = zeta_centre[0] +
400 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
401 // Get the left corner
402 zeta[1] = Theta_positions[0];
403 zeta[2] = Radius_box[0];
404 Volume_pt->position(t, zeta, pos_1);
405
406 // Get the right corner
407 zeta[1] = Theta_positions[1];
408 zeta[2] = Radius_box[1];
409 Volume_pt->position(t, zeta, pos_2);
410 // Now interpolate between the positions
411 lin_interpolate(pos_1, pos_2, s[0], f);
412 break;
413
414 case B:
415 // Get the position on the wall
416 zeta[0] = zeta_centre[0];
417 zeta[1] =
418 Theta_positions[0] +
419 (Theta_positions[1] - Theta_positions[0]) * 0.5 * (s[0] + 1.0);
420 zeta[2] = 1.0;
421 Volume_pt->position(t, zeta, pos_1);
422
423 // Now linearly interpolate the position on the box
424 lin_interpolate(Box[0][0], Box[0][1], s[0], pos_2);
425
426 // Now linearly interpolate between the two
427 lin_interpolate(pos_1, pos_2, s[1], f);
428 break;
429
430 case F:
431 // Get the position on the wall
432 zeta[0] = zeta_centre[1];
433 zeta[1] =
434 Theta_positions[0] +
435 (Theta_positions[1] - Theta_positions[0]) * 0.5 * (s[0] + 1.0);
436 zeta[2] = 1.0;
437 Volume_pt->position(t, zeta, pos_1);
438
439 // Now linearly interpolate the position on the box
440 lin_interpolate(Box[1][0], Box[1][1], s[0], pos_2);
441
442 // Now linearly interpolate between the two
443 lin_interpolate(pos_1, pos_2, s[1], f);
444 break;
445
446 default:
447
448 std::ostringstream error_stream;
449 error_stream << "idirect is " << idirect
450 << " not one of L, R, D, U, B, F" << std::endl;
451
452 throw OomphLibError(error_stream.str(),
453 OOMPH_CURRENT_FUNCTION,
454 OOMPH_EXCEPTION_LOCATION);
455 break;
456 }
457
458 break;
459
460 // Macro element 2:Right
461 case 2:
462
463 // Which direction?
464 switch (idirect)
465 {
466 case L:
467 // Need to get the position from the domain
468 // Get the centreline position
469 zeta[0] = zeta_centre[0] +
470 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
471 // Get the lower corner
472 zeta[1] = Theta_positions[1];
473 zeta[2] = Radius_box[1];
474 Volume_pt->position(t, zeta, pos_1);
475
476 // Get the upper corner
477 zeta[1] = Theta_positions[2];
478 zeta[2] = Radius_box[2];
479 Volume_pt->position(t, zeta, pos_2);
480 // Now interpolate between the positions
481 lin_interpolate(pos_1, pos_2, s[0], f);
482 break;
483
484 case R:
485 // Entirely on the wall
486 zeta[0] = zeta_centre[0] +
487 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
488 zeta[1] =
489 Theta_positions[1] +
490 (Theta_positions[2] - Theta_positions[1]) * 0.5 * (s[0] + 1.0);
491 zeta[2] = 1.0;
492 Volume_pt->position(t, zeta, f);
493 break;
494
495 case U:
496 // Get the position on the wall
497 zeta[0] = zeta_centre[0] +
498 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
499 zeta[1] = Theta_positions[2];
500 zeta[2] = 1.0;
501 Volume_pt->position(t, zeta, pos_1);
502
503 // Get the position of the box
504 zeta[2] = Radius_box[2];
505 Volume_pt->position(t, zeta, pos_2);
506
507 // Now linearly interpolate between the two
508 lin_interpolate(pos_2, pos_1, s[0], f);
509 break;
510
511 case D:
512 // Get the position on the wall
513 zeta[0] = zeta_centre[0] +
514 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
515 zeta[1] = Theta_positions[1];
516 zeta[2] = 1.0;
517 Volume_pt->position(t, zeta, pos_1);
518
519 // Get the position of the box
520 zeta[2] = Radius_box[1];
521 Volume_pt->position(t, zeta, pos_2);
522
523 // Now linearly interpolate between the two
524 lin_interpolate(pos_2, pos_1, s[0], f);
525 break;
526
527 case B:
528 // Get the position on the wall
529 zeta[0] = zeta_centre[0];
530 zeta[1] =
531 Theta_positions[1] +
532 (Theta_positions[2] - Theta_positions[1]) * 0.5 * (s[1] + 1.0);
533 zeta[2] = 1.0;
534 Volume_pt->position(t, zeta, pos_1);
535
536 // Now linearly interpolate the position on the box
537 lin_interpolate(Box[0][1], Box[0][2], s[1], pos_2);
538
539 // Now linearly interpolate between the two
540 lin_interpolate(pos_2, pos_1, s[0], f);
541 break;
542
543 case F:
544 // Get the position on the wall
545 zeta[0] = zeta_centre[1];
546 zeta[1] =
547 Theta_positions[1] +
548 (Theta_positions[2] - Theta_positions[1]) * 0.5 * (s[1] + 1.0);
549 zeta[2] = 1.0;
550 Volume_pt->position(t, zeta, pos_1);
551
552 // Now linearly interpolate the position on the box
553 lin_interpolate(Box[1][1], Box[1][2], s[1], pos_2);
554
555 // Now linearly interpolate between the two
556 lin_interpolate(pos_2, pos_1, s[0], f);
557 break;
558
559 default:
560 std::ostringstream error_stream;
561 error_stream << "idirect is " << idirect
562 << " not one of L, R, D, U, B, F" << std::endl;
563
564 throw OomphLibError(error_stream.str(),
565 OOMPH_CURRENT_FUNCTION,
566 OOMPH_EXCEPTION_LOCATION);
567 }
568
569 break;
570
571 // Macro element 3: Top
572 case 3:
573
574 // Which direction?
575 switch (idirect)
576 {
577 case L:
578 // Get the position on the wall
579 zeta[0] = zeta_centre[0] +
580 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
581 zeta[1] = Theta_positions[3];
582 zeta[2] = 1.0;
583 Volume_pt->position(t, zeta, pos_1);
584
585 // Get the position on the box
586 zeta[2] = Radius_box[3];
587 Volume_pt->position(t, zeta, pos_2);
588
589 // Now linearly interpolate between the two
590 lin_interpolate(pos_2, pos_1, s[0], f);
591 break;
592
593 case R:
594 // Get the position on the wall
595 zeta[0] = zeta_centre[0] +
596 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
597 zeta[1] = Theta_positions[2];
598 zeta[2] = 1.0;
599 Volume_pt->position(t, zeta, pos_1);
600
601 // Get the position on the box
602 zeta[2] = Radius_box[2];
603 Volume_pt->position(t, zeta, pos_2);
604
605 // Now linearly interpolate between the two
606 lin_interpolate(pos_2, pos_1, s[0], f);
607 break;
608
609 case D:
610 // This is entirely on the box
611 // Need to get the position from the domain
612 // Get the centreline position
613 zeta[0] = zeta_centre[0] +
614 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
615 zeta[1] = Theta_positions[3];
616 zeta[2] = Radius_box[3];
617 // Get the lower corner
618 Volume_pt->position(t, zeta, pos_2);
619
620 // Get the upper corner
621 zeta[1] = Theta_positions[2];
622 zeta[2] = Radius_box[2];
623 // Get the upper corner
624 Volume_pt->position(t, zeta, pos_1);
625 // Now interpolate between the positions
626 lin_interpolate(pos_2, pos_1, s[0], f);
627 break;
628
629 case U:
630 // This is entirely on the wall
631 zeta[0] = zeta_centre[0] +
632 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
633 zeta[1] =
634 Theta_positions[3] +
635 (Theta_positions[2] - Theta_positions[3]) * 0.5 * (s[0] + 1.0);
636 zeta[2] = 1.0;
637 Volume_pt->position(t, zeta, f);
638 break;
639
640 case B:
641 // Get the position on the wall
642 zeta[0] = zeta_centre[0];
643 zeta[1] =
644 Theta_positions[3] +
645 (Theta_positions[2] - Theta_positions[3]) * 0.5 * (s[0] + 1.0);
646 zeta[2] = 1.0;
647 Volume_pt->position(t, zeta, pos_1);
648
649 // Now linearly interpolate the position on the box
650 lin_interpolate(Box[0][3], Box[0][2], s[0], pos_2);
651
652 // Now linearly interpolate between the two
653 lin_interpolate(pos_2, pos_1, s[1], f);
654 break;
655
656 case F:
657 // Get the position on the wall
658 zeta[0] = zeta_centre[1];
659 zeta[1] =
660 Theta_positions[3] +
661 (Theta_positions[2] - Theta_positions[3]) * 0.5 * (s[0] + 1.0);
662 zeta[2] = 1.0;
663 Volume_pt->position(t, zeta, pos_1);
664
665 // Now linearly interpolate the position on the box
666 lin_interpolate(Box[1][3], Box[1][2], s[0], pos_2);
667
668 // Now linearly interpolate between the two
669 lin_interpolate(pos_2, pos_1, s[1], f);
670 break;
671
672 default:
673 std::ostringstream error_stream;
674 error_stream << "idirect is " << idirect
675 << " not one of L, R, D, U, B, F" << std::endl;
676
677 throw OomphLibError(error_stream.str(),
678 OOMPH_CURRENT_FUNCTION,
679 OOMPH_EXCEPTION_LOCATION);
680 }
681
682 break;
683
684 // Macro element 4: Left
685 case 4:
686
687 // Which direction?
688 switch (idirect)
689 {
690 case L:
691 // Entirely on the wall, Need to be careful
692 // because this is the point at which theta passes through the
693 // branch cut
694 zeta[0] = zeta_centre[0] +
695 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
696 zeta[1] = Theta_positions[0] + 2.0 * pi +
697 (Theta_positions[3] - (Theta_positions[0] + 2.0 * pi)) *
698 0.5 * (s[0] + 1.0);
699 zeta[2] = 1.0;
700 Volume_pt->position(t, zeta, f);
701 break;
702
703 case R:
704 // Entirely on the box
705 // Need to get the position from the domain
706 // Get the centreline position
707 zeta[0] = zeta_centre[0] +
708 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
709 zeta[1] = Theta_positions[0];
710 zeta[2] = Radius_box[0];
711 // Get the lower corner
712 Volume_pt->position(t, zeta, pos_2);
713
714 // Get the upper corner
715 zeta[1] = Theta_positions[3];
716 zeta[2] = Radius_box[3];
717 // Get the upper corner
718 Volume_pt->position(t, zeta, pos_1);
719
720 lin_interpolate(pos_2, pos_1, s[0], f);
721 break;
722
723 case D:
724 // Get the position on the wall
725 zeta[0] = zeta_centre[0] +
726 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
727 zeta[1] = Theta_positions[0];
728 zeta[2] = 1.0;
729 Volume_pt->position(t, zeta, pos_1);
730
731 // Get the position on the box
732 zeta[2] = Radius_box[0];
733 Volume_pt->position(t, zeta, pos_2);
734
735 // Now linearly interpolate between the two
736 lin_interpolate(pos_1, pos_2, s[0], f);
737 break;
738
739 case U:
740 // Get the position on the wall
741 zeta[0] = zeta_centre[0] +
742 (zeta_centre[1] - zeta_centre[0]) * 0.5 * (s[1] + 1.0);
743 zeta[1] = Theta_positions[3];
744 zeta[2] = 1.0;
745 Volume_pt->position(t, zeta, pos_1);
746
747 // Get the position on the box
748 zeta[2] = Radius_box[3];
749 Volume_pt->position(t, zeta, pos_2);
750
751 // Now linearly interpolate between the two
752 lin_interpolate(pos_1, pos_2, s[0], f);
753 break;
754
755 case B:
756 // Get the position on the wall
757 // Again be careful of the branch cut
758 zeta[0] = zeta_centre[0];
759 zeta[1] = Theta_positions[0] + 2.0 * pi +
760 (Theta_positions[3] - (Theta_positions[0] + 2.0 * pi)) *
761 0.5 * (s[1] + 1.0);
762 zeta[2] = 1.0;
763 Volume_pt->position(t, zeta, pos_1);
764
765 // Now linearly interpolate the position on the box
766 lin_interpolate(Box[0][0], Box[0][3], s[1], pos_2);
767
768 // Now linearly interpolate between the two
769 lin_interpolate(pos_1, pos_2, s[0], f);
770 break;
771
772 case F:
773 // Get the position on the wall
774 // Again be careful of the branch cut
775 zeta[0] = zeta_centre[1];
776 zeta[1] = Theta_positions[0] + 2.0 * pi +
777 (Theta_positions[3] - (Theta_positions[0] + 2.0 * pi)) *
778 0.5 * (s[1] + 1.0);
779 zeta[2] = 1.0;
780 Volume_pt->position(t, zeta, pos_1);
781
782 // Now linearly interpolate the position on the box
783 lin_interpolate(Box[1][0], Box[1][3], s[1], pos_2);
784
785 // Now linearly interpolate between the two
786 lin_interpolate(pos_1, pos_2, s[0], f);
787 break;
788
789 default:
790 std::ostringstream error_stream;
791 error_stream << "idirect is " << idirect
792 << " not one of L, R, D, U, B, F" << std::endl;
793
794 throw OomphLibError(error_stream.str(),
795 OOMPH_CURRENT_FUNCTION,
796 OOMPH_EXCEPTION_LOCATION);
797 }
798 break;
799
800 default:
801 // Error
802 std::ostringstream error_stream;
803 error_stream << "Wrong imacro " << imacro << std::endl;
804 throw OomphLibError(
805 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
806 }
807 }
808
809} // namespace oomph
810
811#endif
s
static char t char * s
Definition cfortran.h:568
i
cstr elem_len * i
Definition cfortran.h:603
t
char t
Definition cfortran.h:568
oomph::Domain
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition domain.h:67
oomph::Domain::Macro_element_pt
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition domain.h:301
oomph::GeomObject
A geometric object is an object that provides a parametrised description of its shape via the functio...
Definition geom_objects.h:101
oomph::GeomObject::position
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
oomph::OomphLibError
An OomphLibError object which should be thrown when an run-time error is encountered....
Definition oomph_definitions.h:229
oomph::OomphLibWarning
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Definition oomph_definitions.h:274
oomph::TAdvectionDiffusionReactionElement
TAdvectionDiffusionReactionElement<NREAGENT,DIM,NNODE_1D> elements are isoparametric triangular DIM-d...
Definition Tadvection_diffusion_reaction_elements.h:66
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:166
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
geom_objects.h
oomph::MathematicalConstants::Pi
const double Pi
50 digits from maple
Definition oomph_utilities.h:167
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition advection_diffusion_elements.cc:30