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