rectangle_with_hole_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: 1194 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-05-27 08:44:43 +0100 (Fri, 27 May 2016) $
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 #ifndef OOMPH_RECTANGLE_WITH_HOLE_DOMAIN_HEADER
31 #define OOMPH_RECTANGLE_WITH_HOLE_DOMAIN_HEADER
32 
33 
34 // Generic includes
35 #include "../generic/quadtree.h"
36 #include "../generic/geom_objects.h"
37 #include "../generic/macro_element.h"
38 #include "../generic/domain.h"
39 
40 
41 namespace oomph
42 {
43 
44 
45 
46 //===========================================================
47 /// Rectangular domain with circular whole
48 //===========================================================
50 {
51 
52 public:
53 
54 
55  /// \short Constructor. Pass pointer to geometric object that
56  /// represents the cylinder, the length of the (square) domain.
57  /// The GeomObject must be parametrised such that
58  /// \f$\zeta \in [0,2\pi]\f$ sweeps around the circumference
59  /// in anticlockwise direction.
61  const double &length) :
62  Cylinder_pt(cylinder_pt)
63  {
64 
65  //Vertices of rectangle
66  Lower_left.resize(2);
67  Lower_left[0] = -0.5*length;
68  Lower_left[1] = -0.5*length;
69 
70  Upper_left.resize(2);
71  Upper_left[0] = -0.5*length;
72  Upper_left[1] = 0.5*length;
73 
74  Lower_right.resize(2);
75  Lower_right[0] = 0.5*length;
76  Lower_right[1] = -0.5*length;
77 
78  Upper_right.resize(2);
79  Upper_right[0] = 0.5*length;
80  Upper_right[1] = 0.5*length;
81 
82 
83  // Coordinates of points where the "radial" lines from central
84  // cylinder meet the upper and lower boundaries
85  Lower_mid_left.resize(2);
86  Lower_mid_left[0] = -0.5*length;
87  Lower_mid_left[1] = -0.5*length;
88 
89  Upper_mid_left.resize(2);
90  Upper_mid_left[0] = -0.5*length;
91  Upper_mid_left[1] = 0.5*length;
92 
93  Lower_mid_right.resize(2);
94  Lower_mid_right[0] = 0.5*length;
95  Lower_mid_right[1] = -0.5*length;
96 
97  Upper_mid_right.resize(2);
98  Upper_mid_right[0] = 0.5*length;
99  Upper_mid_right[1] = 0.5*length;
100 
101 
102  //There are four macro elements
103  Macro_element_pt.resize(4);
104 
105  // Build the 2D macro elements
106  for (unsigned i=0;i<4;i++)
107  {Macro_element_pt[i]= new QMacroElement<2>(this,i);}
108  }
109 
110 
111 
112  /// Destructor: Kill macro elements, why isn't this generic?
114  {
115  for (unsigned i=0;i<4;i++){delete Macro_element_pt[i];}
116  }
117 
118 
119  /// \short Helper function to interpolate linearly between the
120  /// "right" and "left" points; \f$ s \in [-1,1] \f$
122  const double &s, Vector<double> &f)
123  {
124  for(unsigned i=0;i<2;i++)
125  {
126  f[i] = left[i] + (right[i] - left[i])*0.5*(s+1.0);
127  }
128  }
129 
130 
131 
132  /// \short Parametrisation of macro element boundaries: f(s) is the position
133  /// vector to macro-element m's boundary in the specified direction [N/S/E/W]
134  /// at the specfied discrete time level (time=0: present; time>0: previous)
135  void macro_element_boundary(const unsigned &time,
136  const unsigned &m,
137  const unsigned &direction,
138  const Vector<double> &s,
139  Vector<double>& f)
140  {
141 
142 #ifdef WARN_ABOUT_SUBTLY_CHANGED_OOMPH_INTERFACES
143  // Warn about time argument being moved to the front
145  "Order of function arguments has changed between versions 0.8 and 0.85",
146  "RectangleWithHoleDomain::macro_element_boundary(...)",
147  OOMPH_EXCEPTION_LOCATION);
148 #endif
149 
150  // Lagrangian coordinate along surface of cylinder
151  Vector<double> xi(1);
152 
153  // Point on circle
154  Vector<double> point_on_circle(2);
155 
156  using namespace QuadTreeNames;
157 
158  //Switch on the macro element
159  switch(m)
160  {
161 
162  //Macro element 0, is is immediately left of the cylinder
163  case 0:
164 
165  switch(direction)
166  {
167  case N:
168  xi[0] = 3.0*atan(1.0);
169  Cylinder_pt->position(time,xi,point_on_circle);
170  linear_interpolate(Upper_mid_left,point_on_circle,s[0],f);
171  break;
172 
173  case S:
174  xi[0] = -3.0*atan(1.0);
175  Cylinder_pt->position(time,xi,point_on_circle);
176  linear_interpolate(Lower_mid_left,point_on_circle,s[0],f);
177  break;
178 
179  case W:
181  break;
182 
183  case E:
184  xi[0] = 5.0*atan(1.0) - 2.0*atan(1.0)*0.5*(1.0+s[0]);
185  Cylinder_pt->position(time,xi,f);
186  break;
187 
188  default:
189 
190  std::ostringstream error_stream;
191  error_stream << "Direction is incorrect: " << direction << std::endl;
192  throw OomphLibError(error_stream.str(),
193  OOMPH_CURRENT_FUNCTION,
194  OOMPH_EXCEPTION_LOCATION);
195  }
196 
197  break;
198 
199  //Macro element 1, is immediately above the cylinder
200  case 1:
201 
202  switch(direction)
203  {
204  case N:
206  break;
207 
208  case S:
209  xi[0] = 3.0*atan(1.0) - 2.0*atan(1.0)*0.5*(1.0+s[0]);
210  Cylinder_pt->position(time,xi,f);
211  break;
212 
213  case W:
214  xi[0] = 3.0*atan(1.0);
215  Cylinder_pt->position(time,xi,point_on_circle);
216  linear_interpolate(point_on_circle,Upper_mid_left,s[0],f);
217  break;
218 
219  case E:
220  xi[0] = 1.0*atan(1.0);
221  Cylinder_pt->position(time,xi,point_on_circle);
222  linear_interpolate(point_on_circle,Upper_mid_right,s[0],f);
223  break;
224 
225  default:
226 
227  std::ostringstream error_stream;
228  error_stream << "Direction is incorrect: " << direction << std::endl;
229  throw OomphLibError(error_stream.str(),
230  OOMPH_CURRENT_FUNCTION,
231  OOMPH_EXCEPTION_LOCATION);
232  }
233 
234  break;
235 
236  //Macro element 2, is immediately right of the cylinder
237  case 2:
238 
239  switch(direction)
240  {
241  case N:
242  xi[0] = 1.0*atan(1.0);
243  Cylinder_pt->position(time,xi,point_on_circle);
244  linear_interpolate(point_on_circle,Upper_mid_right,s[0],f);
245  break;
246 
247  case S:
248  xi[0] = -1.0*atan(1.0);
249  Cylinder_pt->position(time,xi,point_on_circle);
250  linear_interpolate(point_on_circle,Lower_mid_right,s[0],f);
251  break;
252 
253  case W:
254  xi[0] = -atan(1.0) + 2.0*atan(1.0)*0.5*(1.0+s[0]);
255  Cylinder_pt->position(time,xi,f);
256  break;
257 
258  case E:
260  break;
261 
262  default:
263 
264  std::ostringstream error_stream;
265  error_stream << "Direction is incorrect: " << direction << std::endl;
266  throw OomphLibError(error_stream.str(),
267  OOMPH_CURRENT_FUNCTION,
268  OOMPH_EXCEPTION_LOCATION);
269  }
270 
271  break;
272 
273  //Macro element 3, is immediately below cylinder
274  case 3:
275 
276  switch(direction)
277  {
278  case N:
279  xi[0] = -3.0*atan(1.0) + 2.0*atan(1.0)*0.5*(1.0+s[0]);
280  Cylinder_pt->position(time,xi,f);
281  break;
282 
283  case S:
285  break;
286 
287  case W:
288  xi[0] = -3.0*atan(1.0);
289  Cylinder_pt->position(time,xi,point_on_circle);
290  linear_interpolate(Lower_mid_left,point_on_circle,s[0],f);
291  break;
292 
293  case E:
294  xi[0] = -1.0*atan(1.0);
295  Cylinder_pt->position(time,xi,point_on_circle);
296  linear_interpolate(Lower_mid_right,point_on_circle,s[0],f);
297  break;
298 
299  default:
300 
301  std::ostringstream error_stream;
302  error_stream << "Direction is incorrect: " << direction << std::endl;
303  throw OomphLibError(error_stream.str(),
304  OOMPH_CURRENT_FUNCTION,
305  OOMPH_EXCEPTION_LOCATION);
306  }
307 
308  break;
309 
310  default:
311 
312  std::ostringstream error_stream;
313  error_stream << "Wrong macro element number" << m << std::endl;
314  throw OomphLibError(error_stream.str(),
315  OOMPH_CURRENT_FUNCTION,
316  OOMPH_EXCEPTION_LOCATION);
317  }
318  }
319 
320 
321 private:
322 
323  /// Lower left corner of rectangle
325 
326  /// Lower right corner of rectangle
328 
329  /// Where the "radial" line from circle meets lower boundary on left
331 
332  /// Where the "radial" line from circle meets lower boundary on right
334 
335  /// Upper left corner of rectangle
337 
338  /// Upper right corner of rectangle
340 
341  /// Where the "radial" line from circle meets upper boundary on left
343 
344  /// Where the "radial" line from circle meets upper boundary on right
346 
347  /// Pointer to geometric object that represents the central cylinder
349 
350 };
351 
352 
353 }
354 #endif
Rectangular domain with circular whole.
Vector< double > Lower_right
Lower right corner of rectangle.
cstr elem_len * i
Definition: cfortran.h:607
RectangleWithHoleDomain(GeomObject *cylinder_pt, const double &length)
Constructor. Pass pointer to geometric object that represents the cylinder, the length of the (square...
Vector< MacroElement * > Macro_element_pt
Vector of pointers to macro elements.
Definition: domain.h:237
virtual void position(const Vector< double > &zeta, Vector< double > &r) const =0
Parametrised position on object at current time: r(zeta).
Vector< double > Lower_mid_right
Where the "radial" line from circle meets lower boundary on right.
Vector< double > Lower_left
Lower left corner of rectangle.
~RectangleWithHoleDomain()
Destructor: Kill macro elements, why isn't this generic?
Vector< double > Upper_mid_left
Where the "radial" line from circle meets upper boundary on left.
Vector< double > Lower_mid_left
Where the "radial" line from circle meets lower boundary on left.
static char t char * s
Definition: cfortran.h:572
Vector< double > Upper_right
Upper right corner of rectangle.
Vector< double > Upper_left
Upper left corner of rectangle.
Vector< double > Upper_mid_right
Where the "radial" line from circle meets upper boundary on right.
GeomObject * Cylinder_pt
Pointer to geometric object that represents the central cylinder.
void linear_interpolate(Vector< double > left, Vector< double > right, const double &s, Vector< double > &f)
Helper function to interpolate linearly between the "right" and "left" points; .
Base class for Domains with curvilinear and/or time-dependent boundaries. Domain boundaries are typic...
Definition: domain.h:71
void macro_element_boundary(const unsigned &time, const unsigned &m, const unsigned &direction, const Vector< double > &s, Vector< double > &f)
Parametrisation of macro element boundaries: f(s) is the position vector to macro-element m's boundar...