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