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-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_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
35namespace 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),
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
77 /// Broken copy constructor
78 QuarterPipeDomain(const QuarterPipeDomain&) = delete;
79
80 /// Broken assignment operator
81 void operator=(const QuarterPipeDomain&) = delete;
82
83 /// Destructor:
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
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
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
139 GeomObject* Outer_boundary_cross_section_pt;
140
141 /// Geom object representing the inner boundary of
142 /// the cross section
143 GeomObject* Inner_boundary_cross_section_pt;
144
145 /// Function pointer for function that implements
146 /// axial spacing of macro elements
147 AxialSpacingFctPt Axial_spacing_fct_pt;
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 //=================================================================
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 //=================================================================
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
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)
Function that implements axial spacing of macro elements.
Definition quarter_pipe_domain.h:104
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
Geom object representing the inner boundary of the cross section.
Definition quarter_pipe_domain.h:143
oomph::QuarterPipeDomain::default_axial_spacing_fct
static double default_axial_spacing_fct(const double &xi)
Default for function that implements axial spacing of macro elements.
Definition quarter_pipe_domain.h:151
oomph::QuarterPipeDomain::AxialSpacingFctPt
double(* AxialSpacingFctPt)(const double &xi)
Typedef for function pointer for function that implements axial spacing of macro elements.
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::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::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::axial_spacing_fct_pt
AxialSpacingFctPt & axial_spacing_fct_pt()
Function pointer for function that implements axial spacing of macro elements.
Definition quarter_pipe_domain.h:97
oomph::QuarterPipeDomain::Axial_spacing_fct_pt
AxialSpacingFctPt Axial_spacing_fct_pt
Function pointer for function that implements axial spacing of macro elements.
Definition quarter_pipe_domain.h:147
oomph::QuarterPipeDomain::Outer_boundary_cross_section_pt
GeomObject * Outer_boundary_cross_section_pt
Geom object representing the outer boundary of the cross section.
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)
Vector representation of the i_macro-th macro element boundary i_direct (U/D/L/R/F/B) at time level t...
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)
Constructor: Pass number of elements in various directions, the inner and outer radius and the length...
Definition quarter_pipe_domain.h:45
oomph::RefineableTriangleMesh
/////////////////////////////////////////////////////////////////// /////////////////////////////////...
Definition triangle_mesh.template.h:2249
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition annular_domain.h:35