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