full_circle_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 
31 #ifndef OOMPH_FULL_CIRCLE_DOMAIN_HEADER
32 #define OOMPH_FULL_CIRCLE_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 Topologically circular domain, e.g. a tube cross section.
44 /// The entire domain must be defined by a GeomObject with the
45 /// following convention: zeta[0] is the radial coordinate and
46 /// zeta[1] is the theta coordinate around the cross-sectin.
47 /// The outer boundary must lie at zeta[0] = 1.
48 ///
49 /// The domain is
50 /// parametrised by five macro elements (a central box surrounded by
51 /// four curved elements). The labelling of the macro elements is shown
52 /// below.
53 ///
54 /// ----------------------------
55 /// |\ /|
56 /// | \ Macro / |
57 /// | 3 Element 3 2 |
58 /// | \ / |
59 /// | \----------------/ |
60 /// | | | |
61 /// | 4 | Macro | |
62 /// | | Element 0 | 2 |
63 /// | | | |
64 /// | ----------------- |
65 /// | / \ |
66 /// | 0 Macro 1 |
67 /// | / Element 1 \ |
68 /// | / \|
69 /// |/-------------------------|
70 ///
71 ///
72 //=================================================================
73 class FullCircleDomain : public Domain
74 {
75 
76 public:
77 
78  /// \short Constructor: Pass geometric object; the theta locations
79  /// marking the division between
80  /// the elements of the outer ring, labelled from the lower left to the
81  /// upper left in order, theta should be in the range \f$-\pi\f$ to
82  /// \f$\pi\f$; and the corresponding fractions of the
83  /// radius at which the central box is to be placed.
84  FullCircleDomain(GeomObject* area_geom_object_pt,
85  const Vector<double> &theta_positions,
86  const Vector<double> &radius_box) :
87  Theta_positions(theta_positions), Radius_box(radius_box),
88  Area_pt(area_geom_object_pt)
89  {
90  // There are five macro elements
91  const unsigned n_macro=5;
92  Macro_element_pt.resize(n_macro);
93 
94  // Create the macro elements
95  for (unsigned i=0;i<n_macro;i++)
96  {
98  }
99  }
100 
101  /// Broken copy constructor
103  {
104  BrokenCopy::broken_copy("FullCircleDomain");
105  }
106 
107  /// Broken assignment operator
109  {
110  BrokenCopy::broken_assign("FullCircleDomain");
111  }
112 
113 
114  /// Destructor: Kill all macro elements
116  {
117  for(unsigned i=0;i<5;i++)
118  {
119  delete Macro_element_pt[i];
120  }
121  }
122 
123 
124  /// \short Vector representation of the i_macro-th macro element
125  /// boundary i_direct (N/S/W/E) at time level t
126  /// (t=0: present; t>0: previous):
127  /// f(s).
128  void macro_element_boundary(const unsigned& t,
129  const unsigned& i_macro,
130  const unsigned& i_direct,
131  const Vector<double>& s,
132  Vector<double>& f);
133 
134 private:
135 
136  /// \short Storage for the dividing lines on the boundary
137  /// starting from the lower left and proceeding anticlockwise to
138  /// the upper left
140 
141  /// Storage for the fraction of the radius at which the central box
142  /// should be located corresponding to the chosen values of theta.
144 
145  /// Pointer to geometric object that represents the domain
147 
148  /// \short A very little linear interpolation helper.
149  /// Interpolate from the low
150  /// point to the high point using the coordinate s, which is
151  /// assumed to run from -1 to 1.
153  const Vector<double> &high,
154  const double &s,
155  Vector<double> &f)
156  {
157  //Loop over all coordinates (we are working in 2D)
158  for(unsigned i=0;i<2;i++)
159  {
160  f[i] = low[i] + (high[i] - low[i])*0.5*(s+1.0);
161  }
162  }
163 
164 };
165 
166 
167 /////////////////////////////////////////////////////////////////////////
168 /////////////////////////////////////////////////////////////////////////
169 /////////////////////////////////////////////////////////////////////////
170 
171 
172 
173 //=================================================================
174 /// \short Vector representation of the imacro-th macro element
175 /// boundary idirect (N/S/W/E) at time level t
176 /// (t=0: present; t>0: previous): f(s)
177 //=================================================================
179  const unsigned& t,
180  const unsigned& imacro,
181  const unsigned& idirect,
182  const Vector<double>& s,
183  Vector<double>& f)
184  {
185 
186  using namespace QuadTreeNames;
187 
188  //Get the coordinates of the corners of the box
189  Vector<Vector<double> > Box(4);
190  //Get the corresponding coordinates on the boundary
191  Vector<Vector<double> > Wall(4);
192 
193  //Storage for position in the area
194  Vector<double> zeta(2);
195 
196  //Loop over the angles
197  for(unsigned j=0;j<4;j++)
198  {
199  //Set up the storage
200  Box[j].resize(2);
201  Wall[j].resize(2);
202 
203  //Set the other values of zeta
204  zeta[0] = Radius_box[j];
205  zeta[1] = Theta_positions[j];
206  //Now get the values
207  Area_pt->position(t,zeta,Box[j]);
208 
209  //Now get the position on the boundary
210  zeta[0] = 1.0;
211  Area_pt->position(t,zeta,Wall[j]);
212  }
213 
214  //Define pi
215  const double pi = MathematicalConstants::Pi;
216 
217  // Which macro element?
218  // --------------------
219  switch(imacro)
220  {
221 
222  // Macro element 0: Central box
223  case 0:
224 
225  //Choose a direction
226  switch(idirect)
227  {
228 
229  case N:
230  //Linearly interpolate between the two corners of the box
231  lin_interpolate(Box[3],Box[2],s[0],f);
232  break;
233 
234  case S:
235  //Linearly interpolate between the two corners of the box
236  lin_interpolate(Box[0],Box[1],s[0],f);
237 
238  case W:
239  //Linearly interpolate between the two corners of the box
240  lin_interpolate(Box[0],Box[3],s[0],f);
241  break;
242 
243  case E:
244  //Linearly interpolate between the two corners of the box
245  lin_interpolate(Box[1],Box[2],s[0],f);
246  break;
247 
248  default:
249  std::ostringstream error_stream;
250  error_stream << "idirect is " << idirect
251  << " not one of N, S, E, W" << std::endl;
252 
253  throw OomphLibError(
254  error_stream.str(),
255  OOMPH_CURRENT_FUNCTION,
256  OOMPH_EXCEPTION_LOCATION);
257  break;
258  }
259 
260  break;
261 
262  // Macro element 1: Bottom
263  case 1:
264 
265  //Choose a direction
266  switch(idirect)
267  {
268 
269  case N:
270  //Linearly interpolate between the two corners of the box
271  lin_interpolate(Box[0],Box[1],s[0],f);
272  break;
273 
274  case S:
275  //Get the position on the wall by linearly interpolating in theta
276  zeta[0] = 1.0;
277  zeta[1] = Theta_positions[0] + (Theta_positions[1] - Theta_positions[0])
278  *0.5*(s[0]+1.0);
279 
280  Area_pt->position(t,zeta,f);
281  break;
282 
283  case W:
284  //Now linearly interpolate between the wall and the box
285  lin_interpolate(Wall[0],Box[0],s[0],f);
286  break;
287 
288  case E:
289  //Linearly interpolate between the wall and the box
290  lin_interpolate(Wall[1],Box[1],s[0],f);
291  break;
292 
293  default:
294 
295  std::ostringstream error_stream;
296  error_stream << "idirect is " << idirect
297  << " not one of N, S, E, W" << std::endl;
298 
299  throw OomphLibError(
300  error_stream.str(),
301  OOMPH_CURRENT_FUNCTION,
302  OOMPH_EXCEPTION_LOCATION);
303  break;
304  }
305 
306 
307  break;
308 
309  // Macro element 2:Right
310  case 2:
311 
312  // Which direction?
313  switch(idirect)
314  {
315  case N:
316  //Linearly interpolate between the box and the wall
317  lin_interpolate(Box[2],Wall[2],s[0],f);
318  break;
319 
320  case S:
321  //Linearly interpolate between the box and the wall
322  lin_interpolate(Box[1],Wall[1],s[0],f);
323  break;
324 
325  case W:
326  //Linearly interpolate between the two corners of the box
327  lin_interpolate(Box[1],Box[2],s[0],f);
328  break;
329 
330  case E:
331  //Get the position on the wall by linearly interpolating in theta
332  zeta[0] = 1.0;
333  zeta[1] = Theta_positions[1] + (Theta_positions[2] - Theta_positions[1])
334  *0.5*(s[0]+1.0);
335 
336  Area_pt->position(t,zeta,f);
337  break;
338 
339  default:
340  std::ostringstream error_stream;
341  error_stream << "idirect is " << idirect
342  << " not one of N, S, W, E" << std::endl;
343 
344  throw OomphLibError(
345  error_stream.str(),
346  OOMPH_CURRENT_FUNCTION,
347  OOMPH_EXCEPTION_LOCATION);
348  }
349 
350  break;
351 
352  // Macro element 3: Top
353  case 3:
354 
355  // Which direction?
356  switch(idirect)
357  {
358  case N:
359  //Get the position on the wall by linearly interpolating in theta
360  zeta[0] = 1.0;
361  zeta[1] = Theta_positions[3] + (Theta_positions[2] - Theta_positions[3])
362  *0.5*(s[0]+1.0);
363 
364  Area_pt->position(t,zeta,f);
365  break;
366 
367  case S:
368  //Linearly interpolate between the two corners of the box
369  lin_interpolate(Box[3],Box[1],s[0],f);
370  break;
371 
372  case W:
373  //Linearly interpolate between the box and the wall
374  lin_interpolate(Box[3],Wall[3],s[0],f);
375  break;
376 
377  case E:
378  //Linearly interpolate between the box and the wall
379  lin_interpolate(Box[2],Wall[2],s[0],f);
380  break;
381 
382  default:
383  std::ostringstream error_stream;
384  error_stream << "idirect is " << idirect
385  << " not one of N, S, E, W" << std::endl;
386 
387  throw OomphLibError(
388  error_stream.str(),
389  OOMPH_CURRENT_FUNCTION,
390  OOMPH_EXCEPTION_LOCATION);
391  }
392 
393  break;
394 
395 
396  // Macro element 4: Left
397  case 4:
398 
399  // Which direction?
400  switch(idirect)
401  {
402  case N:
403  //Linearly interpolate between the wall and the box
404  lin_interpolate(Wall[3],Box[3],s[0],f);
405  break;
406 
407  case S:
408  //Linearly interpolate between the wall and the box
409  lin_interpolate(Wall[0],Box[0],s[0],f);
410  break;
411 
412  case W:
413  //Entirely on the wall, Need to be careful
414  //because this is the point at which theta passes through the
415  //"branch cut"
416  zeta[0] = 1.0;
417  zeta[1] = Theta_positions[0] + 2.0*pi +
418  (Theta_positions[3] - (Theta_positions[0] + 2.0*pi))
419  *0.5*(s[0]+1.0);
420 
421  Area_pt->position(t,zeta,f);
422  break;
423 
424  case E:
425  //Linearly interpolate between the two corners of the box
426  lin_interpolate(Box[0],Box[3],s[0],f);
427  break;
428 
429  default:
430  std::ostringstream error_stream;
431  error_stream << "idirect is " << idirect
432  << " not one of N, S, W, E" << std::endl;
433 
434  throw OomphLibError(
435  error_stream.str(),
436  OOMPH_CURRENT_FUNCTION,
437  OOMPH_EXCEPTION_LOCATION);
438  }
439  break;
440 
441 
442  default:
443  // Error
444  std::ostringstream error_stream;
445  error_stream << "Wrong imacro " << imacro << std::endl;
446  throw OomphLibError(
447  error_stream.str(),
448  OOMPH_CURRENT_FUNCTION,
449  OOMPH_EXCEPTION_LOCATION);
450  }
451 
452 }
453 
454 }
455 
456 #endif
457 
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
FullCircleDomain(const FullCircleDomain &)
Broken copy constructor.
Topologically circular domain, e.g. a tube cross section. The entire domain must be defined by a Geom...
void operator=(const FullCircleDomain &)
Broken assignment operator.
cstr elem_len * i
Definition: cfortran.h:607
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
GeomObject * Area_pt
Pointer to geometric object that represents the domain.
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
~FullCircleDomain()
Destructor: Kill all macro elements.
Vector< double > Theta_positions
Storage for the dividing lines on the boundary starting from the lower left and proceeding anticlockw...
static char t char * s
Definition: cfortran.h:572
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...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
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 (N/S/W/E) at time level t (t=...
FullCircleDomain(GeomObject *area_geom_object_pt, const Vector< double > &theta_positions, const Vector< double > &radius_box)
Constructor: Pass geometric object; the theta locations marking the division between the elements of ...
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition: domain.h:71