Toggle navigation
Documentation
Big picture
The finite element method
The data structure
Not-so-quick guide
Optimisation
Order of action functions
Example codes and tutorials
List of example codes and tutorials
Meshing
Solvers
MPI parallel processing
Post-processing/visualisation
Other
Change log
Creating documentation
Coding conventions
Index
FAQ
Installation
Installation guide
Copyright
About
People
Contact/Get involved
Publications
Acknowledgements
Picture show
Go
src
meshes
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-2023 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
35
namespace
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
{
91
Macro_element_pt
[
i
] =
new
QMacroElement<2>
(
this
,
i
);
92
}
93
}
94
95
/// Broken copy constructor
96
FullCircleDomain
(
const
FullCircleDomain
&) =
delete
;
97
98
/// Broken assignment operator
99
void
operator=
(
const
FullCircleDomain
&) =
delete
;
100
101
102
/// Destructor: Empty; cleanup done in base class
103
~FullCircleDomain
() {}
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
120
Vector<double>
Theta_positions
;
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.
124
Vector<double>
Radius_box
;
125
126
/// Pointer to geometric object that represents the domain
127
GeomObject
*
Area_pt
;
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.
133
void
lin_interpolate
(
const
Vector<double>
&
low
,
134
const
Vector<double>
&
high
,
135
const
double
&
s
,
136
Vector<double>
&
f
)
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
//=================================================================
157
void
FullCircleDomain::macro_element_boundary
(
const
unsigned
&
t
,
158
const
unsigned
&
imacro
,
159
const
unsigned
&
idirect
,
160
const
Vector<double>
&
s
,
161
Vector<double>
&
f
)
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
oomph::FullCircleDomain
Topologically circular domain, e.g. a tube cross section. The entire domain must be defined by a Geom...
Definition
full_circle_domain.h:69
oomph::FullCircleDomain::FullCircleDomain
FullCircleDomain(const FullCircleDomain &)=delete
Broken copy constructor.
oomph::FullCircleDomain::FullCircleDomain
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 ...
Definition
full_circle_domain.h:77
oomph::FullCircleDomain::Area_pt
GeomObject * Area_pt
Pointer to geometric object that represents the domain.
Definition
full_circle_domain.h:127
oomph::FullCircleDomain::lin_interpolate
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...
Definition
full_circle_domain.h:133
oomph::FullCircleDomain::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 (N/S/W/E) at time level t (t=...
Definition
full_circle_domain.h:157
oomph::FullCircleDomain::Radius_box
Vector< double > Radius_box
Storage for the fraction of the radius at which the central box should be located corresponding to th...
Definition
full_circle_domain.h:124
oomph::FullCircleDomain::operator=
void operator=(const FullCircleDomain &)=delete
Broken assignment operator.
oomph::FullCircleDomain::Theta_positions
Vector< double > Theta_positions
Storage for the dividing lines on the boundary starting from the lower left and proceeding anticlockw...
Definition
full_circle_domain.h:120
oomph::FullCircleDomain::~FullCircleDomain
~FullCircleDomain()
Destructor: Empty; cleanup done in base class.
Definition
full_circle_domain.h:103
oomph::GeompackQuadMesh
Quadrilateral mesh generator; Uses input from Geompack++. See: http://members.shaw....
Definition
geompack_mesh.template.h:42
oomph
////////////////////////////////////////////////////////////////////// //////////////////////////////...
Definition
annular_domain.h:35