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-2022 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  //================================================================
39  //================================================================
40  class QuarterPipeDomain : public Domain
41  {
42  public:
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),
57  Axial_spacing_fct_pt(&default_axial_spacing_fct)
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 
78  QuarterPipeDomain(const QuarterPipeDomain&) = delete;
79 
81  void operator=(const QuarterPipeDomain&) = delete;
82 
84  ~QuarterPipeDomain()
85  {
86  // Note: macro elements are cleaned up in base class...
87  delete Outer_boundary_cross_section_pt;
88  delete Inner_boundary_cross_section_pt;
89  }
90 
93  typedef double (*AxialSpacingFctPt)(const double& xi);
94 
97  AxialSpacingFctPt& axial_spacing_fct_pt()
98  {
99  return Axial_spacing_fct_pt;
100  }
101 
104  double axial_spacing_fct(const double& xi)
105  {
106  return Axial_spacing_fct_pt(xi);
107  }
108 
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:
120  unsigned Ntheta;
121 
123  unsigned Nr;
124 
126  unsigned Nz;
127 
129  double Rmin;
130 
132  double Rmax;
133 
135  double Length;
136 
139  GeomObject* Outer_boundary_cross_section_pt;
140 
143  GeomObject* Inner_boundary_cross_section_pt;
144 
147  AxialSpacingFctPt Axial_spacing_fct_pt;
148 
151  static double default_axial_spacing_fct(const double& xi)
152  {
153  return xi;
154  }
155 
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 
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 
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 
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 
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 
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  //=================================================================
229  //=================================================================
230  void QuarterPipeDomain::macro_element_boundary(const unsigned& t,
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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  //=================================================================
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
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
D
dominoes D
Definition: Domino.cpp:55
L
MatrixXd L
Definition: LLT_example.cpp:6
R
@ R
Definition: StatisticsVector.h:21
oomph::Domain
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::Ellipse
Definition: geom_objects.h:644
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::Matrix
Definition: matrices.h:74
oomph::OomphLibError
Definition: oomph_definitions.h:222
oomph::QMacroElement< 3 >
Definition: macro_element.h:373
oomph::QuarterPipeDomain
Domain representing a quarter pipe.
Definition: quarter_pipe_domain.h:41
oomph::QuarterPipeDomain::axial_spacing_fct
double axial_spacing_fct(const double &xi)
Definition: quarter_pipe_domain.h:104
oomph::QuarterPipeDomain::Length
double Length
Length.
Definition: quarter_pipe_domain.h:135
oomph::QuarterPipeDomain::Nr
unsigned Nr
Number of elements radial direction.
Definition: quarter_pipe_domain.h:123
oomph::QuarterPipeDomain::r_L
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 .
Definition: quarter_pipe_domain.h:304
oomph::QuarterPipeDomain::Inner_boundary_cross_section_pt
GeomObject * Inner_boundary_cross_section_pt
Definition: quarter_pipe_domain.h:143
oomph::QuarterPipeDomain::default_axial_spacing_fct
static double default_axial_spacing_fct(const double &xi)
Definition: quarter_pipe_domain.h:151
oomph::QuarterPipeDomain::AxialSpacingFctPt
double(* AxialSpacingFctPt)(const double &xi)
Definition: quarter_pipe_domain.h:93
oomph::QuarterPipeDomain::r_U
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 .
Definition: quarter_pipe_domain.h:421
oomph::QuarterPipeDomain::Rmax
double Rmax
Outer radius.
Definition: quarter_pipe_domain.h:132
oomph::QuarterPipeDomain::QuarterPipeDomain
QuarterPipeDomain(const QuarterPipeDomain &)=delete
Broken copy constructor.
oomph::QuarterPipeDomain::r_F
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 .
Definition: quarter_pipe_domain.h:455
oomph::QuarterPipeDomain::~QuarterPipeDomain
~QuarterPipeDomain()
Destructor:
Definition: quarter_pipe_domain.h:84
oomph::QuarterPipeDomain::axial_spacing_fct_pt
AxialSpacingFctPt & axial_spacing_fct_pt()
Definition: quarter_pipe_domain.h:97
oomph::QuarterPipeDomain::Nz
unsigned Nz
Number of elements axial direction.
Definition: quarter_pipe_domain.h:126
oomph::QuarterPipeDomain::r_R
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 .
Definition: quarter_pipe_domain.h:345
oomph::QuarterPipeDomain::operator=
void operator=(const QuarterPipeDomain &)=delete
Broken assignment operator.
oomph::QuarterPipeDomain::r_B
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 .
Definition: quarter_pipe_domain.h:490
oomph::QuarterPipeDomain::Rmin
double Rmin
Inner radius.
Definition: quarter_pipe_domain.h:129
oomph::QuarterPipeDomain::Axial_spacing_fct_pt
AxialSpacingFctPt Axial_spacing_fct_pt
Definition: quarter_pipe_domain.h:147
oomph::QuarterPipeDomain::Outer_boundary_cross_section_pt
GeomObject * Outer_boundary_cross_section_pt
Definition: quarter_pipe_domain.h:139
oomph::QuarterPipeDomain::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)
Definition: quarter_pipe_domain.h:230
oomph::QuarterPipeDomain::r_D
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 .
Definition: quarter_pipe_domain.h:387
oomph::QuarterPipeDomain::Ntheta
unsigned Ntheta
Number of elements azimuthal direction.
Definition: quarter_pipe_domain.h:120
oomph::QuarterPipeDomain::QuarterPipeDomain
QuarterPipeDomain(const unsigned &ntheta, const unsigned &nr, const unsigned &nz, const double &rmin, const double &rmax, const double &length)
Definition: quarter_pipe_domain.h:45
f
static int f(const TensorMap< Tensor< int, 3 > > &tensor)
Definition: cxx11_tensor_map.cpp:237
s
RealScalar s
Definition: level1_cplx_impl.h:130
Eigen::zeta
EIGEN_STRONG_INLINE const Eigen::CwiseBinaryOp< Eigen::internal::scalar_zeta_op< typename DerivedX::Scalar >, const DerivedX, const DerivedQ > zeta(const Eigen::ArrayBase< DerivedX > &x, const Eigen::ArrayBase< DerivedQ > &q)
Definition: SpecialFunctionsArrayAPI.h:152
Global_Parameters::rmax
double rmax
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:70
Global_Parameters::zmax
double zmax
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:71
Global_Parameters::zmin
double zmin
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:69
Global_Parameters::rmin
double rmin
Definition: time_harmonic_fourier_decomposed_linear_elasticity/cylinder/cylinder.cc:68
Mesh_Parameters::nz
const unsigned nz
Definition: ConstraintElementsUnitTest.cpp:32
RachelsAdvectionDiffusion::U
double U
Swimming speed.
Definition: two_d_variable_diff_adapt.cc:53
constants::pi
const Mdouble pi
Definition: ExtendedMath.h:23
oomph::MathematicalConstants::Pi
const double Pi
50 digits from maple
Definition: oomph_utilities.h:157
oomph::OcTreeNames::F
@ F
Definition: octree.h:74
oomph
DRAIG: Change all instances of (SPATIAL_DIM) to (DIM-1).
Definition: AnisotropicHookean.h:10
plotDoE.x
list x
Definition: plotDoE.py:28
plotPSD.t
t
Definition: plotPSD.py:36
OOMPH_EXCEPTION_LOCATION
#define OOMPH_EXCEPTION_LOCATION
Definition: oomph_definitions.h:61
OOMPH_CURRENT_FUNCTION
#define OOMPH_CURRENT_FUNCTION
Definition: oomph_definitions.h:86