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-2024 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 
35 namespace 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
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
/////////////////////////////////////////////////////////////////////
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:222
oomph::OomphLibWarning
An OomphLibWarning object which should be created as a temporary object to issue a warning....
Definition: oomph_definitions.h:267
oomph::QMacroElement< 3 >
QMacroElement specialised to 3 spatial dimensions.
Definition: macro_element.h:373
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::MathematicalConstants::Pi
const double Pi
50 digits from maple
Definition: oomph_utilities.h:157
oomph
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
Definition: advection_diffusion_elements.cc:30