quarter_pipe_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_QUARTER_PIPE_DOMAIN_HEADER
28 #define OOMPH_QUARTER_PIPE_DOMAIN_HEADER
29 
30 // Generic oomph-lib includes
31 #include "../generic/octree.h"
32 #include "../generic/domain.h"
33 #include "../generic/geom_objects.h"
34 
35 namespace oomph
36 {
37  //================================================================
38  /// Domain representing a quarter pipe
39  //================================================================
40  class QuarterPipeDomain : public Domain
41  {
42  public:
43  /// Constructor: Pass number of elements in various directions,
44  /// the inner and outer radius and the length of the tube
45  QuarterPipeDomain(const unsigned& ntheta,
46  const unsigned& nr,
47  const unsigned& nz,
48  const double& rmin,
49  const double& rmax,
50  const double& length)
51  : Ntheta(ntheta),
52  Nr(nr),
53  Nz(nz),
54  Rmin(rmin),
55  Rmax(rmax),
56  Length(length),
58  {
59  // Number of macroelements
60  unsigned nmacro = nr * ntheta * nz;
61 
62  // Resize
63  Macro_element_pt.resize(nmacro);
64 
65  // Create macro elements
66  for (unsigned i = 0; i < nmacro; i++)
67  {
68  Macro_element_pt[i] = new QMacroElement<3>(this, i);
69  }
70 
71  // Make geom object representing the outer and inner boundaries of
72  // the cross section
73  Inner_boundary_cross_section_pt = new Ellipse(rmin, rmin);
74  Outer_boundary_cross_section_pt = new Ellipse(rmax, rmax);
75  }
76 
77  /// Broken copy constructor
79 
80  /// Broken assignment operator
81  void operator=(const QuarterPipeDomain&) = delete;
82 
83  /// Destructor:
85  {
86  // Note: macro elements are cleaned up in base class...
89  }
90 
91  /// Typedef for function pointer for function that implements
92  /// axial spacing of macro elements
93  typedef double (*AxialSpacingFctPt)(const double& xi);
94 
95  /// Function pointer for function that implements
96  /// axial spacing of macro elements
98  {
99  return Axial_spacing_fct_pt;
100  }
101 
102  /// Function that implements
103  /// axial spacing of macro elements
104  double axial_spacing_fct(const double& xi)
105  {
106  return Axial_spacing_fct_pt(xi);
107  }
108 
109  /// Vector representation of the i_macro-th macro element
110  /// boundary i_direct (U/D/L/R/F/B) at time level t
111  /// (t=0: present; t>0: previous): f(s).
112  void macro_element_boundary(const unsigned& t,
113  const unsigned& i_macro,
114  const unsigned& i_direct,
115  const Vector<double>& s,
116  Vector<double>& f);
117 
118  private:
119  /// Number of elements azimuthal direction
120  unsigned Ntheta;
121 
122  /// Number of elements radial direction
123  unsigned Nr;
124 
125  /// Number of elements axial direction
126  unsigned Nz;
127 
128  /// Inner radius
129  double Rmin;
130 
131  /// Outer radius
132  double Rmax;
133 
134  /// Length
135  double Length;
136 
137  /// Geom object representing the outer boundary of
138  /// the cross section
140 
141  /// Geom object representing the inner boundary of
142  /// the cross section
144 
145  /// Function pointer for function that implements
146  /// axial spacing of macro elements
148 
149  /// Default for function that implements
150  /// axial spacing of macro elements
151  static double default_axial_spacing_fct(const double& xi)
152  {
153  return xi;
154  }
155 
156  /// Boundary of macro element zeta \f$ \in [-1,1]x[-1,1] \f$
157  void r_U(const unsigned& t,
158  const Vector<double>& zeta,
159  Vector<double>& f,
160  const double& rmin,
161  const double& rmax,
162  const double& thetamin,
163  const double& thetamax,
164  const double& zmin,
165  const double& zmax);
166 
167  /// Boundary of macro element zeta \f$ \in [-1,1]x[-1,1] \f$
168  void r_L(const unsigned& t,
169  const Vector<double>& zeta,
170  Vector<double>& f,
171  const double& rmin,
172  const double& rmax,
173  const double& thetamin,
174  const double& thetamax,
175  const double& zmin,
176  const double& zmax);
177 
178  /// Boundary of macro element zeta \f$ \in [-1,1]x[-1,1] \f$
179  void r_D(const unsigned& t,
180  const Vector<double>& zeta,
181  Vector<double>& f,
182  const double& rmin,
183  const double& rmax,
184  const double& thetamin,
185  const double& thetamax,
186  const double& zmin,
187  const double& zmax);
188 
189  /// Boundary of macro element zeta \f$ \in [-1,1]x[-1,1] \f$
190  void r_R(const unsigned& t,
191  const Vector<double>& zeta,
192  Vector<double>& f,
193  const double& rmin,
194  const double& rmax,
195  const double& thetamin,
196  const double& thetamax,
197  const double& zmin,
198  const double& zmax);
199 
200  /// Boundary of macro element zeta \f$ \in [-1,1]x[-1,1] \f$
201  void r_F(const unsigned& t,
202  const Vector<double>& zeta,
203  Vector<double>& f,
204  const double& rmin,
205  const double& rmax,
206  const double& thetamin,
207  const double& thetamax,
208  const double& zmin,
209  const double& zmax);
210 
211  /// Boundary of macro element zeta \f$ \in [-1,1]x[-1,1] \f$
212  void r_B(const unsigned& t,
213  const Vector<double>& zeta,
214  Vector<double>& f,
215  const double& rmin,
216  const double& rmax,
217  const double& thetamin,
218  const double& thetamax,
219  const double& zmin,
220  const double& zmax);
221 
222  }; // endofclass
223 
224 
225  //=================================================================
226  /// Vector representation of the imacro-th macro element
227  /// boundary idirect (U/D/L/R/F/B) at time level t:
228  /// f(s)
229  //=================================================================
231  const unsigned& imacro,
232  const unsigned& idirect,
233  const Vector<double>& s,
234  Vector<double>& f)
235  {
236  using namespace OcTreeNames;
237 
238  const double pi = MathematicalConstants::Pi;
239 
240  // Match the elements number with the position
241  unsigned num_z = imacro / (Nr * Ntheta);
242  unsigned num_y = (imacro % (Nr * Ntheta)) / Ntheta;
243  unsigned num_x = imacro % Ntheta;
244 
245  // Define the extreme coordinates
246 
247  // radial direction
248  double rmin = Rmin + (Rmax - Rmin) * double(num_y) / double(Nr);
249  double rmax = Rmin + (Rmax - Rmin) * double(num_y + 1) / double(Nr);
250 
251  // theta direction
252  double thetamin = (pi / 2.0) * (1.0 - double(num_x + 1) / double(Ntheta));
253  double thetamax = (pi / 2.0) * (1.0 - double(num_x) / double(Ntheta));
254 
255  // zdirection (tube)
256  double zmin = double(num_z) * Length / double(Nz);
257  double zmax = double(num_z + 1) * Length / double(Nz);
258 
259 
260  // Which direction?
261  if (idirect == U)
262  {
263  r_U(t, s, f, rmin, rmax, thetamin, thetamax, zmin, zmax);
264  }
265  else if (idirect == D)
266  {
267  r_D(t, s, f, rmin, rmax, thetamin, thetamax, zmin, zmax);
268  }
269  else if (idirect == L)
270  {
271  r_L(t, s, f, rmin, rmax, thetamin, thetamax, zmin, zmax);
272  }
273  else if (idirect == R)
274  {
275  r_R(t, s, f, rmin, rmax, thetamin, thetamax, zmin, zmax);
276  }
277  else if (idirect == F)
278  {
279  r_F(t, s, f, rmin, rmax, thetamin, thetamax, zmin, zmax);
280  }
281  else if (idirect == B)
282  {
283  r_B(t, s, f, rmin, rmax, thetamin, thetamax, zmin, zmax);
284  }
285  else
286  {
287  std::ostringstream error_stream;
288  error_stream << "idirect is " << idirect << " not one of U, D, L, R, F, B"
289  << std::endl;
290 
291  throw OomphLibError(
292  error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
293  }
294 
295  // Now redistribute points in the axial direction
296  double z_frac = f[2] / Length;
297  f[2] = Length * axial_spacing_fct(z_frac);
298  }
299 
300 
301  //=================================================================
302  /// Left face of a macro element \f$ s \in [-1,1]*[-1,1] \f$
303  //=================================================================
304  void QuarterPipeDomain::r_L(const unsigned& t,
305  const Vector<double>& s,
306  Vector<double>& f,
307  const double& rmin,
308  const double& rmax,
309  const double& thetamin,
310  const double& thetamax,
311  const double& zmin,
312  const double& zmax)
313  {
314  Vector<double> x(1);
315  x[0] = thetamax;
316 
317  // Point on outer wall
318  Vector<double> r_outer(2);
319  Outer_boundary_cross_section_pt->position(t, x, r_outer);
320 
321  // Point on inner wall
322  Vector<double> r_inner(2);
323  Inner_boundary_cross_section_pt->position(t, x, r_inner);
324 
325  // Get layer boundaries
326  Vector<double> r_top(2);
327  Vector<double> r_bot(2);
328  for (unsigned i = 0; i < 2; i++)
329  {
330  r_top[i] =
331  r_inner[i] + (r_outer[i] - r_inner[i]) * (rmax - Rmin) / (Rmax - Rmin);
332  r_bot[i] =
333  r_inner[i] + (r_outer[i] - r_inner[i]) * (rmin - Rmin) / (Rmax - Rmin);
334  }
335 
336  // Compute coordinates
337  f[0] = r_bot[0] + (0.5 * (s[0] + 1.0)) * (r_top[0] - r_bot[0]);
338  f[1] = r_bot[1] + (0.5 * (s[0] + 1.0)) * (r_top[1] - r_bot[1]);
339  f[2] = zmin + (zmax - zmin) * (0.5 * (s[1] + 1.0));
340  }
341 
342  //=================================================================
343  /// Right face of a macro element \f$ s \in [-1,1]*[-1,1] \f$
344  //=================================================================
345  void QuarterPipeDomain::r_R(const unsigned& t,
346  const Vector<double>& s,
347  Vector<double>& f,
348  const double& rmin,
349  const double& rmax,
350  const double& thetamin,
351  const double& thetamax,
352  const double& zmin,
353  const double& zmax)
354  {
355  Vector<double> x(1);
356  x[0] = thetamin;
357 
358  // Point on outer wall
359  Vector<double> r_outer(2);
360  Outer_boundary_cross_section_pt->position(t, x, r_outer);
361 
362  // Point on inner wall
363  Vector<double> r_inner(2);
364  Inner_boundary_cross_section_pt->position(t, x, r_inner);
365 
366  // Get layer boundaries
367  Vector<double> r_top(2);
368  Vector<double> r_bot(2);
369  for (unsigned i = 0; i < 2; i++)
370  {
371  r_top[i] =
372  r_inner[i] + (r_outer[i] - r_inner[i]) * (rmax - Rmin) / (Rmax - Rmin);
373  r_bot[i] =
374  r_inner[i] + (r_outer[i] - r_inner[i]) * (rmin - Rmin) / (Rmax - Rmin);
375  }
376 
377  // Compute coordinates
378  f[0] = r_bot[0] + (0.5 * (s[0] + 1.0)) * (r_top[0] - r_bot[0]);
379  f[1] = r_bot[1] + (0.5 * (s[0] + 1.0)) * (r_top[1] - r_bot[1]);
380  f[2] = zmin + (zmax - zmin) * (0.5 * (s[1] + 1.0));
381  }
382 
383 
384  //=================================================================
385  /// Left face of a macro element \f$s \in [-1,1]*[-1,1] \f$
386  //=================================================================
387  void QuarterPipeDomain::r_D(const unsigned& t,
388  const Vector<double>& s,
389  Vector<double>& f,
390  const double& rmin,
391  const double& rmax,
392  const double& thetamin,
393  const double& thetamax,
394  const double& zmin,
395  const double& zmax)
396  {
397  Vector<double> x(1);
398  x[0] = thetamax + (0.5 * (s[0] + 1.0)) * (thetamin - thetamax);
399 
400  // Point on outer wall
401  Vector<double> r_outer(2);
402  Outer_boundary_cross_section_pt->position(t, x, r_outer);
403 
404  // Point on inner wall
405  Vector<double> r_inner(2);
406  Inner_boundary_cross_section_pt->position(t, x, r_inner);
407 
408  // Get layer
409  for (unsigned i = 0; i < 2; i++)
410  {
411  f[i] =
412  r_inner[i] + (r_outer[i] - r_inner[i]) * (rmin - Rmin) / (Rmax - Rmin);
413  }
414  f[2] = zmin + (zmax - zmin) * (0.5 * (s[1] + 1.0));
415  }
416 
417 
418  //=================================================================
419  /// Right face of a macro element \f$ s \in [-1,1]*[-1,1] \f$
420  //=================================================================
421  void QuarterPipeDomain::r_U(const unsigned& t,
422  const Vector<double>& s,
423  Vector<double>& f,
424  const double& rmin,
425  const double& rmax,
426  const double& thetamin,
427  const double& thetamax,
428  const double& zmin,
429  const double& zmax)
430  {
431  Vector<double> x(1);
432  x[0] = thetamax + (0.5 * (s[0] + 1.0)) * (thetamin - thetamax);
433 
434  // Point on outer wall
435  Vector<double> r_outer(2);
436  Outer_boundary_cross_section_pt->position(t, x, r_outer);
437 
438  // Point on inner wall
439  Vector<double> r_inner(2);
440  Inner_boundary_cross_section_pt->position(t, x, r_inner);
441 
442  // Get layer
443  for (unsigned i = 0; i < 2; i++)
444  {
445  f[i] =
446  r_inner[i] + (r_outer[i] - r_inner[i]) * (rmax - Rmin) / (Rmax - Rmin);
447  }
448  f[2] = zmin + (zmax - zmin) * (0.5 * (s[1] + 1.0));
449  }
450 
451 
452  //=================================================================
453  /// Front face of a macro element \f$ s \in [-1,1]*[-1,1] \f$
454  //=================================================================
455  void QuarterPipeDomain::r_F(const unsigned& t,
456  const Vector<double>& s,
457  Vector<double>& f,
458  const double& rmin,
459  const double& rmax,
460  const double& thetamin,
461  const double& thetamax,
462  const double& zmin,
463  const double& zmax)
464  {
465  Vector<double> x(1);
466  x[0] = thetamax + (0.5 * (s[0] + 1.0)) * (thetamin - thetamax);
467 
468  // Point on outer wall
469  Vector<double> r_outer(2);
470  Outer_boundary_cross_section_pt->position(t, x, r_outer);
471 
472  // Point on inner wall
473  Vector<double> r_inner(2);
474  Inner_boundary_cross_section_pt->position(t, x, r_inner);
475 
476  // Get layer
477  double rad = rmin + (0.5 * (s[1] + 1.0)) * (rmax - rmin);
478  for (unsigned i = 0; i < 2; i++)
479  {
480  f[i] =
481  r_inner[i] + (r_outer[i] - r_inner[i]) * (rad - Rmin) / (Rmax - Rmin);
482  }
483  f[2] = zmax;
484  }
485 
486 
487  //=================================================================
488  /// Back face of a macro element \f$ s \in [-1,1]*[-1,1] \f$
489  //=================================================================
490  void QuarterPipeDomain::r_B(const unsigned& t,
491  const Vector<double>& s,
492  Vector<double>& f,
493  const double& rmin,
494  const double& rmax,
495  const double& thetamin,
496  const double& thetamax,
497  const double& zmin,
498  const double& zmax)
499  {
500  Vector<double> x(1);
501  x[0] = thetamax + (0.5 * (s[0] + 1.0)) * (thetamin - thetamax);
502 
503  // Point on outer wall
504  Vector<double> r_outer(2);
505  Outer_boundary_cross_section_pt->position(t, x, r_outer);
506 
507  // Point on inner wall
508  Vector<double> r_inner(2);
509  Inner_boundary_cross_section_pt->position(t, x, r_inner);
510 
511  // Get layer
512  double rad = rmin + (0.5 * (s[1] + 1.0)) * (rmax - rmin);
513  for (unsigned i = 0; i < 2; i++)
514  {
515  f[i] =
516  r_inner[i] + (r_outer[i] - r_inner[i]) * (rad - Rmin) / (Rmax - Rmin);
517  }
518  f[2] = zmin;
519  }
520 
521 
522 } // namespace oomph
523 
524 #endif
Domain representing a quarter pipe.
double axial_spacing_fct(const double &xi)
Function that implements axial spacing of macro elements.
unsigned Nr
Number of elements radial direction.
void r_L(const unsigned &t, const Vector< double > &zeta, Vector< double > &f, const double &rmin, const double &rmax, const double &thetamin, const double &thetamax, const double &zmin, const double &zmax)
Boundary of macro element zeta .
GeomObject * Inner_boundary_cross_section_pt
Geom object representing the inner boundary of the cross section.
static double default_axial_spacing_fct(const double &xi)
Default for function that implements axial spacing of macro elements.
double(* AxialSpacingFctPt)(const double &xi)
Typedef for function pointer for function that implements axial spacing of macro elements.
void r_U(const unsigned &t, const Vector< double > &zeta, Vector< double > &f, const double &rmin, const double &rmax, const double &thetamin, const double &thetamax, const double &zmin, const double &zmax)
Boundary of macro element zeta .
QuarterPipeDomain(const QuarterPipeDomain &)=delete
Broken copy constructor.
void r_F(const unsigned &t, const Vector< double > &zeta, Vector< double > &f, const double &rmin, const double &rmax, const double &thetamin, const double &thetamax, const double &zmin, const double &zmax)
Boundary of macro element zeta .
AxialSpacingFctPt & axial_spacing_fct_pt()
Function pointer for function that implements axial spacing of macro elements.
unsigned Nz
Number of elements axial direction.
void r_R(const unsigned &t, const Vector< double > &zeta, Vector< double > &f, const double &rmin, const double &rmax, const double &thetamin, const double &thetamax, const double &zmin, const double &zmax)
Boundary of macro element zeta .
void operator=(const QuarterPipeDomain &)=delete
Broken assignment operator.
void r_B(const unsigned &t, const Vector< double > &zeta, Vector< double > &f, const double &rmin, const double &rmax, const double &thetamin, const double &thetamax, const double &zmin, const double &zmax)
Boundary of macro element zeta .
AxialSpacingFctPt Axial_spacing_fct_pt
Function pointer for function that implements axial spacing of macro elements.
GeomObject * Outer_boundary_cross_section_pt
Geom object representing the outer boundary of the cross section.
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 (U/D/L/R/F/B) at time level t...
void r_D(const unsigned &t, const Vector< double > &zeta, Vector< double > &f, const double &rmin, const double &rmax, const double &thetamin, const double &thetamax, const double &zmin, const double &zmax)
Boundary of macro element zeta .
unsigned Ntheta
Number of elements azimuthal direction.
QuarterPipeDomain(const unsigned &ntheta, const unsigned &nr, const unsigned &nz, const double &rmin, const double &rmax, const double &length)
Constructor: Pass number of elements in various directions, the inner and outer radius and the length...
////////////////////////////////////////////////////////////////////// //////////////////////////////...