ryujin 2.1.1 revision 350e54cc11f3d21282bddcf3e6517944dbc508bf
geometry_airfoil.h
Go to the documentation of this file.
1//
2// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
3// Copyright (C) 2020 - 2023 by the ryujin authors
4//
5
6#pragma once
7
8#include "geometry_common_includes.h"
9
10#include "cubic_spline.h"
11#include "transfinite_interpolation.template.h"
12
13namespace ryujin
14{
15 namespace
16 {
23 template <typename T, typename T2>
24 inline void assign(T &array,
25 const std::initializer_list<T2> &initializer_list)
26 {
27 /* accomodate for a possible c-style array... */
28 Assert(std::size(array) == std::size(initializer_list),
29 dealii::ExcMessage(
30 "size of initializer list and array does not match"));
31 std::copy(
32 initializer_list.begin(), initializer_list.end(), std::begin(array));
33 }
34 } // namespace
35
36 namespace Manifolds
37 {
41 template <int dim>
42 class AirfoilManifold : public dealii::ChartManifold<dim>
43 {
44 static_assert(dim == 2, "not implemented for dim != 2");
45
46 public:
47 AirfoilManifold(const dealii::Point<dim> airfoil_center,
48 const std::function<double(const double)> &psi_front,
49 const std::function<double(const double)> &psi_upper,
50 const std::function<double(const double)> &psi_lower,
51 const bool upper_side,
52 const double psi_ratio = 1.0)
53 : airfoil_center(airfoil_center)
54 , psi_front(psi_front)
55 , psi_upper(psi_upper)
56 , psi_lower(psi_lower)
57 , upper_side(upper_side)
58 , ratio_(psi_ratio * psi_front(0.) / psi_front(M_PI))
59 , polar_manifold()
60 {
61 Assert(std::abs(psi_upper(0.) - psi_front(0.5 * M_PI)) < 1.0e-10,
62 dealii::ExcInternalError());
63 Assert(std::abs(psi_lower(0.) + psi_front(1.5 * M_PI)) < 1.0e-10,
64 dealii::ExcInternalError());
65 }
66
67 dealii::Point<dim>
68 pull_back(const dealii::Point<dim> &space_point) const final
69 {
70 auto coordinate = dealii::Point<dim>() + (space_point - airfoil_center);
71
72 /* transform: */
73
74 dealii::Point<dim> chart_point;
75 if (coordinate[0] > 0.) {
76 if (upper_side) {
77 /* upper back airfoil part */
78 chart_point[0] = 1. + coordinate[1] - psi_upper(coordinate[0]);
79 chart_point[1] = 0.5 * M_PI - ratio_ * coordinate[0];
80 } else {
81 /* lower back airfoil part */
82 chart_point[0] = 1. - coordinate[1] + psi_lower(coordinate[0]);
83 chart_point[1] = 1.5 * M_PI + ratio_ * coordinate[0];
84 }
85 } else {
86 /* front part */
87 chart_point = polar_manifold.pull_back(coordinate);
88 chart_point[0] = 1. + chart_point[0] - psi_front(chart_point[1]);
89 }
90
91 return chart_point;
92 }
93
94 dealii::Point<dim>
95 push_forward(const dealii::Point<dim> &point) const final
96 {
97 auto chart_point = point;
98
99 /* transform back */
100
101 dealii::Point<dim> coordinate;
102 if (chart_point[1] < 0.5 * M_PI) {
103 Assert(upper_side, dealii::ExcInternalError());
104 /* upper back airfoil part */
105 coordinate[0] = (0.5 * M_PI - chart_point[1]) / ratio_;
106 Assert(coordinate[0] >= -1.0e-10, dealii::ExcInternalError());
107 coordinate[1] = chart_point[0] - 1. + psi_upper(coordinate[0]);
108 } else if (chart_point[1] > 1.5 * M_PI) {
109 Assert(!upper_side, dealii::ExcInternalError());
110 /* lower back airfoil part */
111 coordinate[0] = (chart_point[1] - 1.5 * M_PI) / ratio_;
112 Assert(coordinate[0] >= -1.0e-10, dealii::ExcInternalError());
113 coordinate[1] = 1. - chart_point[0] + psi_lower(coordinate[0]);
114 } else {
115 /* front part */
116 chart_point[0] = chart_point[0] - 1. + psi_front(chart_point[1]);
117 coordinate = polar_manifold.push_forward(chart_point);
118 }
119
120 return dealii::Point<dim>() + (coordinate + airfoil_center);
121 }
122
123 std::unique_ptr<dealii::Manifold<dim, dim>> clone() const final
124 {
125 const double psi_ratio = ratio_ * psi_front(M_PI) / psi_front(0.);
126 return std::make_unique<AirfoilManifold<dim>>(airfoil_center,
127 psi_front,
128 psi_upper,
129 psi_lower,
130 upper_side,
131 psi_ratio);
132 }
133
134 private:
135 const dealii::Point<dim> airfoil_center;
136 const std::function<double(const double)> psi_front;
137 const std::function<double(const double)> psi_upper;
138 const std::function<double(const double)> psi_lower;
139 const bool upper_side;
140
141 const double ratio_;
142
143 dealii::PolarManifold<dim> polar_manifold;
144 };
145
146
150 template <int dim>
151 class GradingManifold : public dealii::ChartManifold<dim>
152 {
153 public:
154 GradingManifold(const dealii::Point<dim> center,
155 const dealii::Tensor<1, dim> direction,
156 const double grading,
157 const double epsilon)
158 : center(center)
159 , direction(direction)
160 , grading(grading)
161 , epsilon(epsilon)
162 {
163 }
164
165 /* FIXME: find out why weights are not normalized. */
166 Point<dim>
167 get_new_point(const ArrayView<const Point<dim>> &surrounding_points,
168 const ArrayView<const double> &weights) const override
169 {
170 if (weights[0] > 1.0)
171 return surrounding_points[0];
172
173 if (weights[1] > 1.0)
174 return surrounding_points[1];
175
176 return dealii::ChartManifold<dim>::get_new_point(surrounding_points,
177 weights);
178 }
179
180 dealii::Point<dim>
181 pull_back(const dealii::Point<dim> &space_point) const final
182 {
183 auto chart_point = space_point - center;
184
185 for (unsigned int d = 0; d < dim; ++d) {
186 if (std::abs(direction[d]) > 1.0e-10) {
187 const double x = chart_point[d] * std::copysign(1., direction[d]);
188 Assert(x + epsilon > 0, dealii::ExcInternalError());
189 const double x_hat = std::pow(x + epsilon, 1. / grading) -
190 std::pow(epsilon, 1. / grading) + 1.e-14;
191 chart_point[d] += (x_hat - x) * std::copysign(1., direction[d]);
192 }
193 }
194
195 return dealii::Point<dim>() + chart_point;
196 }
197
198 dealii::Point<dim>
199 push_forward(const dealii::Point<dim> &chart_point) const final
200 {
201 auto space_point = chart_point;
202
203 for (unsigned int d = 0; d < dim; ++d) {
204 if (std::abs(direction[d]) > 1.0e-10) {
205 const double x_hat =
206 space_point[d] * std::copysign(1., direction[d]);
207 Assert(x_hat + std::pow(epsilon, 1. / grading) > 0,
208 dealii::ExcInternalError());
209 const double x =
210 std::pow(x_hat + std::pow(epsilon, 1. / grading), grading) -
211 epsilon + 1.e-14;
212 space_point[d] += (x - x_hat) * std::copysign(1., direction[d]);
213 }
214 }
215
216 return center + (space_point - dealii::Point<dim>());
217 }
218
219 std::unique_ptr<dealii::Manifold<dim, dim>> clone() const final
220 {
221 return std::make_unique<GradingManifold<dim>>(
222 center, direction, grading, epsilon);
223 }
224
225 private:
226 const dealii::Point<dim> center;
227 const dealii::Tensor<1, dim> direction;
228 const double grading;
229 const double epsilon;
230 };
231
232
236 template <int dim>
237 class ExtrudedManifold : public dealii::Manifold<dim>
238 {
239 public:
240 ExtrudedManifold(const dealii::Manifold<dim - 1> &manifold)
241 : manifold(manifold.clone())
242 {
243 }
244
245 std::unique_ptr<Manifold<dim>> clone() const override
246 {
247 return std::make_unique<ExtrudedManifold<dim>>(*manifold);
248 }
249
250 Point<dim>
251 get_new_point(const ArrayView<const Point<dim>> &surrounding_points,
252 const ArrayView<const double> &weights) const override
253 {
254 Assert(surrounding_points.size() == weights.size(),
255 dealii::ExcInternalError());
256
257 boost::container::small_vector<dealii::Point<dim - 1>, 100>
258 surrounding_points_projected;
259 std::transform(surrounding_points.begin(),
260 surrounding_points.end(),
261 surrounding_points_projected.begin(),
262 [](const dealii::Point<dim> &source) {
263 dealii::Point<dim - 1> result;
264 for (unsigned int d = 0; d < dim - 1; ++d)
265 result[d] = source[d];
266 return result;
267 });
268
269 const auto projected = manifold->get_new_point(
270 ArrayView<const Point<dim - 1>>{surrounding_points_projected.data(),
271 weights.size()},
272 weights);
273
274 dealii::Point<dim> result;
275
276 for (unsigned int d = 0; d < dim - 1; ++d)
277 result[d] = projected[d];
278
279 for (unsigned int i = 0; i < weights.size(); ++i)
280 result[dim - 1] += weights[i] * surrounding_points[i][dim - 1];
281
282 return result;
283 }
284
285 private:
286 std::unique_ptr<const dealii::Manifold<dim - 1>> manifold;
287 };
288
289 } // namespace Manifolds
290
291
292 namespace
293 {
297 std::array<std::vector<double>, 4>
298 naca_4digit_points(const std::string &serial_number,
299 const unsigned int n_samples)
300 {
301 AssertThrow(serial_number.size() == 4,
302 dealii::ExcMessage("Invalid NACA 4 digit serial number"));
303 std::array<unsigned int, 4> digit;
304 std::transform(serial_number.begin(),
305 serial_number.end(),
306 digit.begin(),
307 [](auto it) { return it - '0'; });
308
309 /* thickness */
310 const double t = 0.1 * digit[2] + 0.01 * digit[3];
311 AssertThrow(t > 0.,
312 dealii::ExcMessage("Invalid NACA 4 digit serial number"));
313
314 /* maximal chamber */
315 const double m = 0.01 * digit[0];
316 /* x position of maximal chamber */
317 const double p = 0.1 * digit[1];
318
319 std::vector<double> x_upper;
320 std::vector<double> y_upper;
321 std::vector<double> x_lower;
322 std::vector<double> y_lower;
323
324 for (unsigned int i = 0; i < n_samples; i++) {
325 const double x = 1. * i / (n_samples - 1);
326 const double y =
327 5. * t *
328 (0.2969 * std::sqrt(x) +
329 x * (-0.126 + x * (-0.3516 + x * (0.2843 + x * (-0.1036)))));
330
331 const double y_c = (x < p) ? m / (p * p) * (2. * p * x - x * x)
332 : m / ((1. - p) * (1. - p)) *
333 (1. - 2. * p + 2. * p * x - x * x);
334
335 const double dy_c = (x < p) ? 2. * m / (p * p) * (p - x)
336 : 2. * m / ((1. - p) * (1. - p)) * (p - x);
337
338 const double theta = std::atan(dy_c);
339
340 x_upper.emplace_back(x - y * std::sin(theta));
341 y_upper.emplace_back(y_c + y * std::cos(theta));
342 x_lower.emplace_back(x + y * std::sin(theta));
343 y_lower.emplace_back(y_c - y * std::cos(theta));
344 }
345
346 /* Fix up roundoff errors: */
347 y_upper.front() = 0.;
348 y_upper.back() = 0.;
349 y_lower.front() = 0.;
350 y_lower.back() = 0.;
351
352 return {{x_upper, y_upper, x_lower, y_lower}};
353 }
354
355
365 std::array<std::vector<double>, 4>
366 nasa_sc2(const std::string &serial_number)
367 {
368 if (serial_number == "0714") {
369 std::vector<double> x_upper{
370 .0, .002, .005, .01, .02, .03, .04, .05, .07, .1, .12, .15,
371 .17, .2, .22, .25, .27, .3, .33, .35, .38, .4, .43, .45,
372 .48, .50, .53, .55, .57, .6, .62, .65, .68, .7, .72, .75,
373 .77, .8, .82, .85, .87, .9, .92, .95, .97, .98, .99, 1.};
374
375 std::vector<double> y_upper{
376 .0, .0095, .0158, .0219, .0293, .0343, .0381, .0411,
377 .0462, .0518, .0548, .0585, .0606, .0632, .0646, .0664,
378 .0673, .0685, .0692, .0696, .0698, .0697, .0695, .0692,
379 .0684, .0678, .0666, .0656, .0645, .0625, .0610, .0585,
380 .0555, .0533, .0509, .0469, .0439, .0389, .0353, .0294,
381 .0251, .0181, .0131, .0049, -.0009, -.0039, -.0071, -.0104};
382
383 std::vector<double> x_lower{
384 .0, .002, .005, .01, .02, .03, .04, .05, .07, .1, .12, .15, .17,
385 .20, .22, .25, .28, .3, .32, .35, .37, .4, .42, .45, .48, .5,
386 .53, .55, .58, .6, .63, .65, .68, .70, .73, .75, .77, .80, .83,
387 .85, .87, .89, .92, .94, .95, .96, .97, .98, .99, 1.};
388
389 std::vector<double> y_lower{
390 .0, -.0093, -.016, -.0221, -.0295, -.0344, -.0381, -.0412,
391 -.0462, -.0517, -.0547, -.0585, -.0606, -.0633, -.0647, -.0666,
392 -.068, -.0687, -.0692, -.0696, -.0696, -.0692, -.0688, -.0676,
393 -.0657, -.0644, -.0614, -.0588, -.0543, -.0509, -.0451, -.041,
394 -.0346, -.0302, -.0235, -.0192, -.0150, -.0093, -.0048, -.0024,
395 -.0013, -.0008, -.0016, -.0035, -.0049, -.0066, -.0085, -.0109,
396 -.0137, -.0163};
397
398 return {{x_upper, y_upper, x_lower, y_lower}};
399
400 } else {
401
402 AssertThrow(false,
403 dealii::ExcMessage("Invalid NASA SC(2) serial number"));
404 }
405 }
406
407
416 std::array<std::vector<double>, 4> onera(const std::string &serial_number)
417 {
418 if (serial_number == "OAT15a") {
419 std::vector<double> x_upper{
420 0., 2.95888e-05, 0.000117865, 0.000263239, 0.000464245,
421 0.000719821, 0.00103013, 0.00139605, 0.00181884, 0.00230024,
422 0.00284243, 0.00344764, 0.00411805, 0.00485595, 0.00566349,
423 0.00654241, 0.00749421, 0.0085197, 0.0096197, 0.0107945,
424 0.0120442, 0.013369, 0.0147688, 0.0162438, 0.0177939,
425 0.0194194, 0.0211207, 0.0228982, 0.0247526, 0.0266845,
426 0.028695, 0.0307852, 0.0329562, 0.0352094, 0.0375463,
427 0.0399687, 0.0424782, 0.0450769, 0.0477668, 0.05055,
428 0.0534291, 0.0564063, 0.0594842, 0.0626655, 0.0659531,
429 0.0693498, 0.0728588, 0.0764831, 0.0802261, 0.0840914,
430 0.0880824, 0.0922027, 0.0964564, 0.100847, 0.10538,
431 0.110058, 0.114885, 0.119868, 0.125009, 0.130314,
432 0.135789, 0.14139, 0.147074, 0.152839, 0.158682,
433 0.164603, 0.170599, 0.17667, 0.182814, 0.189028,
434 0.195312, 0.201665, 0.208083, 0.214567, 0.221115,
435 0.227724, 0.234394, 0.241123, 0.24791, 0.254753,
436 0.26165, 0.2686, 0.275601, 0.282653, 0.289753,
437 0.2969, 0.304093, 0.311329, 0.318609, 0.325929,
438 0.333289, 0.340686, 0.348121, 0.35559, 0.363093,
439 0.370629, 0.378194, 0.385789, 0.393412, 0.40106,
440 0.408734, 0.41643, 0.424148, 0.431886, 0.439643,
441 0.447417, 0.455207, 0.46301, 0.470827, 0.478654,
442 0.486491, 0.494337, 0.502189, 0.510046, 0.517907,
443 0.525769, 0.533632, 0.541494, 0.549354, 0.557209,
444 0.565059, 0.572902, 0.580736, 0.588559, 0.596371,
445 0.60417, 0.611953, 0.61972, 0.627469, 0.635198,
446 0.642906, 0.65059, 0.658251, 0.665886, 0.673493,
447 0.68107, 0.688617, 0.69613, 0.70361, 0.711054,
448 0.71846, 0.725827, 0.733154, 0.740438, 0.747679,
449 0.754875, 0.762025, 0.769127, 0.776181, 0.783185,
450 0.790139, 0.79704, 0.80389, 0.810685, 0.817426,
451 0.82411, 0.830738, 0.837307, 0.843817, 0.850265,
452 0.856652, 0.862974, 0.869233, 0.87538, 0.881373,
453 0.887216, 0.892913, 0.898467, 0.903883, 0.909163,
454 0.914311, 0.91933, 0.924224, 0.928996, 0.933648,
455 0.938183, 0.942606, 0.946917, 0.95112, 0.955217,
456 0.959212, 0.963107, 0.966904, 0.970605, 0.974213,
457 0.977731, 0.98116, 0.984503, 0.987762, 0.990939,
458 0.994036, 0.997056, 1.};
459
460 std::vector<double> y_upper{
461 0., 0.000899353, 0.0018231, 0.00276894, 0.00373508,
462 0.00472011, 0.0057226, 0.00674103, 0.0077738, 0.00881906,
463 0.00987467, 0.0109383, 0.0120074, 0.0130793, 0.0141511,
464 0.01522, 0.0162832, 0.0173387, 0.0183841, 0.0194179,
465 0.0204389, 0.021446, 0.0224386, 0.0234164, 0.0243794,
466 0.0253276, 0.0262612, 0.0271805, 0.028086, 0.0289783,
467 0.0298578, 0.0307252, 0.0315811, 0.0324262, 0.0332611,
468 0.0340861, 0.0349022, 0.0357098, 0.0365095, 0.0373016,
469 0.0380867, 0.0388652, 0.0396375, 0.0404039, 0.041165,
470 0.0419211, 0.0426722, 0.0434189, 0.0441614, 0.0448995,
471 0.0456336, 0.0463636, 0.0470894, 0.047811, 0.0485286,
472 0.0492423, 0.0499518, 0.0506574, 0.0513591, 0.0520569,
473 0.0527506, 0.0534343, 0.0541023, 0.054755, 0.0553921,
474 0.0560138, 0.05662, 0.0572108, 0.0577861, 0.0583462,
475 0.0588909, 0.0594202, 0.0599341, 0.0604325, 0.0609153,
476 0.0613826, 0.0618341, 0.06227, 0.06269, 0.0630941,
477 0.0634823, 0.0638544, 0.0642103, 0.06455, 0.0648734,
478 0.0651806, 0.0654713, 0.0657454, 0.0660031, 0.0662442,
479 0.0664685, 0.066676, 0.0668664, 0.0670398, 0.067196,
480 0.0673349, 0.0674562, 0.0675598, 0.0676456, 0.0677134,
481 0.0677629, 0.067794, 0.0678065, 0.0678, 0.0677743,
482 0.0677293, 0.0676646, 0.0675798, 0.0674748, 0.0673492,
483 0.0672027, 0.0670349, 0.0668456, 0.0666344, 0.066401,
484 0.066145, 0.0658661, 0.065564, 0.0652385, 0.064889,
485 0.0645151, 0.0641169, 0.0636938, 0.0632454, 0.0627715,
486 0.0622718, 0.061746, 0.0611937, 0.0606145, 0.0600083,
487 0.0593747, 0.0587136, 0.0580244, 0.0573069, 0.0565607,
488 0.0557853, 0.0549807, 0.0541461, 0.0532814, 0.0523863,
489 0.0514606, 0.0505046, 0.0495188, 0.0485042, 0.047462,
490 0.0463943, 0.0453031, 0.0441914, 0.0430618, 0.0419174,
491 0.0407612, 0.0395961, 0.038425, 0.0372503, 0.0360742,
492 0.034899, 0.0337262, 0.0325572, 0.0313935, 0.030236,
493 0.029086, 0.0279442, 0.0268114, 0.0256966, 0.0246079,
494 0.023545, 0.0225073, 0.0214947, 0.0205065, 0.0195422,
495 0.0186011, 0.0176826, 0.0167863, 0.0159112, 0.0150567,
496 0.0142221, 0.0134066, 0.0126095, 0.0118301, 0.0110678,
497 0.0103219, 0.00959177, 0.00887695, 0.00817697, 0.00749135,
498 0.00681977, 0.0061621, 0.00551806, 0.00488739, 0.00427007,
499 0.00366612, 0.00307588, 0.0024997};
500
501 std::vector<double> x_lower{
502 0., 3.22311e-05, 0.000136327, 0.000324365, 0.000606007,
503 0.000986654, 0.00146626, 0.00204126, 0.00270571, 0.00345312,
504 0.00427708, 0.00517234, 0.00613428, 0.00715943, 0.00824517,
505 0.00938973, 0.0105917, 0.0118504, 0.0131652, 0.0145362,
506 0.0159635, 0.0174476, 0.0189891, 0.0205889, 0.0222479,
507 0.0239675, 0.0257488, 0.0275934, 0.0295028, 0.031479,
508 0.0335237, 0.035639, 0.0378271, 0.04009, 0.0424303,
509 0.0448502, 0.0473523, 0.049939, 0.0526132, 0.0553774,
510 0.0582343, 0.0611868, 0.0642377, 0.0673899, 0.0706465,
511 0.0740105, 0.077485, 0.0810733, 0.0847788, 0.088605,
512 0.0925553, 0.0966336, 0.100844, 0.10519, 0.109675,
513 0.114305, 0.119084, 0.124016, 0.129106, 0.134358,
514 0.139779, 0.145328, 0.150961, 0.156676, 0.162473,
515 0.168349, 0.174303, 0.180333, 0.186437, 0.192615,
516 0.198863, 0.205181, 0.211568, 0.21802, 0.224537,
517 0.231117, 0.237759, 0.24446, 0.251219, 0.258035,
518 0.264905, 0.271829, 0.278804, 0.285828, 0.292901,
519 0.30002, 0.307184, 0.314391, 0.321639, 0.328927,
520 0.336253, 0.343615, 0.351011, 0.358441, 0.365902,
521 0.373392, 0.38091, 0.388455, 0.396024, 0.403616,
522 0.41123, 0.418864, 0.426517, 0.434187, 0.441871,
523 0.44957, 0.457282, 0.465005, 0.472737, 0.480477,
524 0.488225, 0.495977, 0.503733, 0.511493, 0.519252,
525 0.527012, 0.53477, 0.542525, 0.550276, 0.55802,
526 0.565759, 0.573488, 0.581208, 0.588918, 0.596615,
527 0.604298, 0.611967, 0.61962, 0.627257, 0.634874,
528 0.642473, 0.65005, 0.657605, 0.665137, 0.672645,
529 0.680128, 0.687584, 0.695012, 0.702411, 0.70978,
530 0.717118, 0.724424, 0.731697, 0.738935, 0.746137,
531 0.753303, 0.76043, 0.767518, 0.774565, 0.78157,
532 0.788531, 0.795447, 0.802316, 0.809136, 0.815905,
533 0.822623, 0.829286, 0.835893, 0.842441, 0.84893,
534 0.855357, 0.86172, 0.868018, 0.874204, 0.880239,
535 0.886123, 0.891863, 0.89746, 0.902919, 0.908242,
536 0.913433, 0.918495, 0.923431, 0.928244, 0.932938,
537 0.937515, 0.941978, 0.94633, 0.950574, 0.954712,
538 0.958747, 0.962681, 0.966518, 0.970259, 0.973907,
539 0.977464, 0.980932, 0.984314, 0.987612, 0.990827,
540 0.993963, 0.997019, 1.};
541
542 std::vector<double> y_lower{
543 0., -0.000899234, -0.00182108, -0.00275889, -0.00370397,
544 -0.00464681, -0.00557909, -0.00649452, -0.0073895, -0.00826265,
545 -0.00911444, -0.00994611, -0.0107598, -0.0115578, -0.0123422,
546 -0.0131152, -0.0138788, -0.0146348, -0.0153851, -0.0161313,
547 -0.0168749, -0.0176173, -0.0183597, -0.0191032, -0.0198489,
548 -0.0205974, -0.0213498, -0.0221064, -0.0228678, -0.0236341,
549 -0.0244053, -0.0251817, -0.0259627, -0.0267482, -0.0275378,
550 -0.0283309, -0.029127, -0.0299255, -0.0307258, -0.0315273,
551 -0.0323295, -0.0331322, -0.0339346, -0.0347366, -0.035538,
552 -0.0363383, -0.0371378, -0.037936, -0.0387331, -0.0395288,
553 -0.0403228, -0.0411152, -0.0419054, -0.0426933, -0.0434778,
554 -0.0442587, -0.0450347, -0.045805, -0.0465681, -0.0473225,
555 -0.0480666, -0.0487929, -0.0494941, -0.0501695, -0.0508179,
556 -0.0514387, -0.052031, -0.0525942, -0.0531277, -0.0536309,
557 -0.0541035, -0.054545, -0.0549547, -0.0553323, -0.0556773,
558 -0.0559891, -0.056267, -0.0565107, -0.0567191, -0.0568918,
559 -0.0570281, -0.0571272, -0.0571886, -0.0572116, -0.0571958,
560 -0.0571405, -0.0570453, -0.0569098, -0.0567339, -0.0565171,
561 -0.0562594, -0.0559608, -0.055621, -0.0552405, -0.0548193,
562 -0.0543577, -0.0538564, -0.0533157, -0.0527365, -0.0521197,
563 -0.0514661, -0.0507767, -0.0500525, -0.0492947, -0.0485041,
564 -0.0476819, -0.0468294, -0.0459474, -0.0450371, -0.0440995,
565 -0.0431356, -0.0421465, -0.0411333, -0.0400971, -0.039039,
566 -0.0379601, -0.0368617, -0.0357451, -0.0346114, -0.0334621,
567 -0.0322983, -0.0311217, -0.0299336, -0.0287354, -0.0275285,
568 -0.0263145, -0.025095, -0.0238717, -0.0226459, -0.0214194,
569 -0.020194, -0.0189714, -0.0177534, -0.016542, -0.015339,
570 -0.0141466, -0.0129671, -0.0118026, -0.0106558, -0.00952898,
571 -0.00842491, -0.00734634, -0.00629613, -0.00527714, -0.0042922,
572 -0.00334424, -0.00243619, -0.00157086, -0.000750868, 2.13187e-05,
573 0.000743489, 0.00141373, 0.00203031, 0.00259174, 0.00309711,
574 0.00354599, 0.00393819, 0.00427381, 0.00455352, 0.00477839,
575 0.00494994, 0.00506989, 0.00514012, 0.0051629, 0.00514197,
576 0.0050815, 0.0049854, 0.00485738, 0.0047009, 0.00451919,
577 0.0043152, 0.00409157, 0.00385071, 0.00359475, 0.00332561,
578 0.00304503, 0.00275452, 0.00245546, 0.00214904, 0.00183633,
579 0.00151827, 0.00119574, 0.00086945, 0.00054002, 0.000208013,
580 -0.000126036, -0.000461602, -0.000798383, -0.0011363, -0.00147529,
581 -0.00181444, -0.00215832, -0.0024967};
582
583 return {{x_upper, y_upper, x_lower, y_lower}};
584
585 } else {
586
587 AssertThrow(false, dealii::ExcMessage("Invalid ONERA serial number"));
588 }
589 }
590
600 std::array<std::vector<double>, 4> bell(const std::string &serial_number)
601 {
602 if (serial_number == "NLR-1T") {
603 std::vector<double> x_upper{
604 .0, .00259, .00974, .02185, .03796, .05675, .07753,
605 .09845, .12341, .15412, .18767, .22313, .26054, .29979,
606 .34064, .38269, .42528, .46849, .51162, .55383, .59596,
607 .63728, .67732, .71079, .73905, .76946, .80263, .84055,
608 .87846, .90845, .93589, .96199, 1.};
609
610 std::vector<double> y_upper{
611 .0, .00704, .01524, .02296, .02972, .03588, .04098,
612 .04469, .04741, .04986, .05188, .05345, .05459, .05531,
613 .05565, .0556, .05518, .05438, .05323, .05175, .04992,
614 .04774, .04524, .04291, .04017, .03644, .0314, .02533,
615 .01901, .01421, .0102, .00651, .00104};
616
617 std::vector<double> x_lower{
618 .0, .00259, .00974, .02185, .03796, .05675, .07753,
619 .09845, .12341, .15412, .18767, .22313, .26054, .29979,
620 .34064, .38269, .42528, .46849, .51162, .55383, .59596,
621 .63728, .67732, .71079, .73905, .76946, .80263, .84055,
622 .87846, .90845, .93589, .96199, 1.};
623 std::vector<double> y_lower{
624 .0, -.00512, -.00867, -.0118, -.01465, -.01713, -.01929,
625 -.02112, -.02299, -.02494, -.02671, -.02821, -.02944, -.0304,
626 -.03104, -.03142, -.0315, -.03132, -.0308, -.02992, -.02867,
627 -.02734, -.0258, -.02432, -.02305, -.02164, -.01996, -.01794,
628 -.01571, -.01364, -.01087, -.00711, -.00104};
629
630 return {{x_upper, y_upper, x_lower, y_lower}};
631
632 } else {
633
634 AssertThrow(false, dealii::ExcMessage("Invalid BELL serial number"));
635 }
636 }
637
638
642 std::array<std::function<double(const double)>, 3>
643 create_psi(const std::vector<double> &x_upper [[maybe_unused]],
644 const std::vector<double> &y_upper [[maybe_unused]],
645 const std::vector<double> &x_lower [[maybe_unused]],
646 const std::vector<double> &y_lower [[maybe_unused]],
647 const double x_center [[maybe_unused]],
648 const double y_center [[maybe_unused]],
649 const double scaling [[maybe_unused]] = 1.)
650 {
651 Assert(x_upper.size() >= 2, dealii::ExcInternalError());
652 Assert(x_upper.front() == 0. && x_upper.back() == 1.,
653 dealii::ExcInternalError());
654 Assert(std::is_sorted(x_upper.begin(), x_upper.end()),
655 dealii::ExcInternalError());
656
657 Assert(x_lower.size() >= 2, dealii::ExcInternalError());
658 Assert(x_lower.front() == 0. && x_lower.back() == 1.,
659 dealii::ExcInternalError());
660 Assert(std::is_sorted(x_lower.begin(), x_lower.end()),
661 dealii::ExcInternalError());
662
663 Assert(y_upper.size() == x_upper.size(), dealii::ExcInternalError());
664 Assert(y_upper.front() == 0., dealii::ExcInternalError());
665
666 Assert(y_lower.size() == x_lower.size(), dealii::ExcInternalError());
667 Assert(y_lower.front() == 0., dealii::ExcInternalError());
668
669 Assert(y_lower.back() < y_upper.back(), dealii::ExcInternalError());
670
671 Assert(0. < x_center && x_center < 1., dealii::ExcInternalError());
672
673#ifdef DEAL_II_WITH_GSL
674 CubicSpline upper_airfoil(x_upper, y_upper);
675 auto psi_upper =
676 [upper_airfoil, x_center, y_center, scaling](const double x_hat) {
677 /* Past the trailing edge return the the final upper y position: */
678 const double x = x_hat / scaling;
679 if (x > 1. - x_center)
680 return scaling * (upper_airfoil.eval(1.0) - y_center);
681 return scaling * (upper_airfoil.eval(x + x_center) - y_center);
682 };
683
684 CubicSpline lower_airfoil(x_lower, y_lower);
685
686 auto psi_lower =
687 [lower_airfoil, x_center, y_center, scaling](const double x_hat) {
688 /* Past the trailing edge return the the final lower y position: */
689 const double x = x_hat / scaling;
690 if (x > 1. - x_center)
691 return scaling * (lower_airfoil.eval(1.0) - y_center);
692 return scaling * (lower_airfoil.eval(x + x_center) - y_center);
693 };
694
695 /*
696 * Create a combined point set for psi_front:
697 */
698
699 std::vector<double> x_combined;
700 std::vector<double> y_combined;
701
702 for (std::size_t i = 0; i < x_upper.size(); ++i) {
703 if (x_upper[i] >= x_center)
704 break;
705 x_combined.push_back(x_upper[i]);
706 y_combined.push_back(y_upper[i]);
707 }
708
709 /*
710 * We are about to create a spline interpolation in polar coordinates
711 * for the front part. In order to blend this interpolation with the
712 * two splines for the upper and lower part that we have just created
713 * we have to add some additional sample points around the
714 * coordinates were we glue together
715 */
716 for (double x : {x_center, x_center + 0.01, x_center + 0.02}) {
717 x_combined.push_back(x);
718 y_combined.push_back(upper_airfoil.eval(x));
719 }
720
721 std::reverse(x_combined.begin(), x_combined.end());
722 std::reverse(y_combined.begin(), y_combined.end());
723 x_combined.pop_back();
724 y_combined.pop_back();
725
726 for (std::size_t i = 0; i < x_lower.size(); ++i) {
727 if (x_lower[i] >= x_center)
728 break;
729 x_combined.push_back(x_lower[i]);
730 y_combined.push_back(y_lower[i]);
731 }
732
733 for (double x : {x_center, x_center + 0.01, x_center + 0.02}) {
734 x_combined.push_back(x);
735 y_combined.push_back(lower_airfoil.eval(x));
736 }
737
738 /* Translate into polar coordinates: */
739
740 for (unsigned int i = 0; i < y_combined.size(); ++i) {
741 const auto x = x_combined[i] - x_center;
742 const auto y = y_combined[i] - y_center;
743
744 const auto rho = std::sqrt(x * x + y * y);
745 auto phi = std::atan2(y, x);
746 if (phi < 0)
747 phi += 2 * dealii::numbers::PI;
748
749 x_combined[i] = phi;
750 y_combined[i] = rho;
751 }
752
753 /* Ensure that x_combined is monotonically increasing: */
754 if (x_combined.back() == 0.)
755 x_combined.back() = 2. * dealii::numbers::PI;
756 Assert(std::is_sorted(x_combined.begin(), x_combined.end()),
757 dealii::ExcInternalError());
758
759 CubicSpline front_airfoil(x_combined, y_combined);
760 auto psi_front = [front_airfoil, x_center, scaling](const double phi) {
761 /* By convention we return the "back length" for phi == 0.: */
762 if (phi == 0.)
763 return scaling * (1. - x_center);
764
765 return scaling * front_airfoil.eval(phi);
766 };
767
768 return {{psi_front, psi_upper, psi_lower}};
769#else
770 AssertThrow(
771 false,
772 dealii::ExcNotImplemented("Airfoil grid needs deal.II with GSL"));
773 return {};
774#endif
775 }
776
777
778 } // namespace
779
780
781 namespace Geometries
782 {
822 template <int dim>
823 class Airfoil : public Geometry<dim>
824 {
825 public:
826 Airfoil(const std::string subsection)
827 : Geometry<dim>("airfoil", subsection)
828 {
829 /* Parameters affecting parameterization: */
830
831 airfoil_type_ = "NASA SC(2) 0714";
832 this->add_parameter(
833 "airfoil type", airfoil_type_, "airfoil type and serial number");
834
835 airfoil_length_ = 2.;
836 this->add_parameter("airfoil length",
837 airfoil_length_,
838 "length of airfoil (leading to trailing edge)");
839
840 psi_samples_ = 100;
841 this->add_parameter("psi samples",
842 psi_samples_,
843 "number of samples used for generating spline psi");
844
845 psi_center_[0] = 0.05;
846 this->add_parameter("psi center",
847 psi_center_,
848 "center position of airfoil for sampling psi");
849
850 psi_ratio_ = 0.30;
851 this->add_parameter(
852 "psi ratio",
853 psi_ratio_,
854 "Scaling parameter for averages in curved nose region, can be "
855 "adjusted by hand to equliabrate the size of faces at the nose "
856 "part of the airfoil");
857
858
859 airfoil_center_[0] = -.5;
860 this->add_parameter("airfoil center",
861 airfoil_center_,
862 "position of airfoil center in the mesh");
863
864 /* Parameters affecting mesh generation: */
865
866 grading_ = 5.5;
867 this->add_parameter(
868 "grading exponent", grading_, "graded mesh: exponent");
869
870 grading_epsilon_ = 0.02;
871 this->add_parameter("grading epsilon",
872 grading_epsilon_,
873 "graded mesh: regularization parameter");
874
875 grading_epsilon_trailing_ = 0.01;
876 this->add_parameter(
877 "grading epsilon trailing",
878 grading_epsilon_trailing_,
879 "graded mesh: regularization parameter for trailing cells");
880
881 height_ = 6.;
882 this->add_parameter(
883 "height", height_, "height of computational domain");
884
885 width_ = 1.;
886 this->add_parameter("width", width_, "width of computational domain");
887
888 n_anisotropic_refinements_airfoil_ = 1;
889 this->add_parameter(
890 "anisotropic pre refinement airfoil",
891 n_anisotropic_refinements_airfoil_,
892 "number of anisotropic pre refinement steps for the airfoil");
893
894 n_anisotropic_refinements_trailing_ = 3;
895 this->add_parameter("anisotropic pre refinement trailing",
896 n_anisotropic_refinements_trailing_,
897 "number of anisotropic pre refinement steps for "
898 "the blunt trailing edge cell");
899
900 subdivisions_z_ = 2;
901 this->add_parameter("subdivisions z",
902 subdivisions_z_,
903 "number of subdivisions in z direction");
904 }
905
906 void create_triangulation(dealii::Triangulation<dim> &triangulation) final
907 {
908 /*
909 * Step 1: Create parametrization:
910 *
911 * Runtime parameters: airfoil_type_, airfoil_length_, psi_center_,
912 * psi_samples_
913 */
914
915 const auto [x_upper, y_upper, x_lower, y_lower] = [&]() {
916 if (airfoil_type_.rfind("NACA ", 0) == 0) {
917 return naca_4digit_points(airfoil_type_.substr(5), psi_samples_);
918 } else if (airfoil_type_.rfind("NASA SC(2) ", 0) == 0) {
919 return nasa_sc2(airfoil_type_.substr(11));
920 } else if (airfoil_type_.rfind("ONERA ", 0) == 0) {
921 return onera(airfoil_type_.substr(6));
922 } else if (airfoil_type_.rfind("BELL ", 0) == 0) {
923 return bell(airfoil_type_.substr(5));
924 }
925 AssertThrow(false, ExcMessage("Unknown airfoil type"));
926 }();
927
928 const auto [psi_front, psi_upper, psi_lower] =
929 create_psi(x_upper,
930 y_upper,
931 x_lower,
932 y_lower,
933 psi_center_[0],
934 psi_center_[1],
935 airfoil_length_);
936
937 /*
938 * Step 2: Create coarse mesh.
939 *
940 * Runtime parameters: airfoil_center_, height_,
941 */
942
943 /* The radius of the radial front part of the mesh: */
944 const auto outer_radius = 0.5 * height_;
945
946 /* by convention, psi_front(0.) returns the "back length" */
947 const auto back_length = psi_front(0.);
948
949 /* sharp trailing edge? */
950 const bool sharp_trailing_edge =
951 std::abs(psi_upper(back_length) - psi_lower(back_length)) < 1.0e-10;
952 AssertThrow(
953 sharp_trailing_edge ||
954 std::abs(psi_upper(back_length) - psi_lower(back_length)) >
955 0.001 * back_length,
956 dealii::ExcMessage("Blunt trailing edge has a width of less than "
957 "0.1% of the trailing airfoil length."));
958
959 /* Front part: */
960 dealii::Triangulation<2> tria_front;
961
962 {
963 std::vector<dealii::Point<2>> vertices{
964 {-outer_radius, 0.0}, // 0
965 {airfoil_center_[0] - psi_front(M_PI), airfoil_center_[1]}, // 1
966 {-0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 2
967 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 3
968 {0., airfoil_center_[1] + psi_lower(-airfoil_center_[0])}, // 4
969 {airfoil_center_[0] + back_length, //
970 airfoil_center_[1] + psi_lower(back_length)}, // 5
971 {0., airfoil_center_[1] + psi_upper(-airfoil_center_[0])}, // 6
972 {-0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 7
973 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 8
974 };
975
976 std::vector<dealii::CellData<2>> cells(4);
977 assign(cells[0].vertices, {2, 3, 4, 5});
978 assign(cells[1].vertices, {0, 2, 1, 4});
979 assign(cells[2].vertices, {7, 0, 6, 1});
980 if (sharp_trailing_edge) {
981 assign(cells[3].vertices, {8, 7, 5, 6});
982 } else {
983 vertices.push_back({airfoil_center_[0] + back_length,
984 airfoil_center_[1] + psi_upper(back_length)});
985 assign(cells[3].vertices, {8, 7, 9, 6});
986 }
987
988 tria_front.create_triangulation(
989 vertices, cells, dealii::SubCellData());
990 }
991
992 /* Back part: */
993 dealii::Triangulation<2> tria_back;
994
995 if (sharp_trailing_edge) {
996 /* Back part for sharp trailing edge: */
997
998 const std::vector<dealii::Point<2>> vertices{
999 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 0
1000 {airfoil_center_[0] + back_length,
1001 airfoil_center_[1] + psi_lower(back_length)}, // 1
1002 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 2
1003 {outer_radius, -0.5 * outer_radius}, // 3
1004 {outer_radius, 0.0}, // 4
1005 {outer_radius, 0.5 * outer_radius}, // 5
1006 };
1007
1008 std::vector<dealii::CellData<2>> cells(2);
1009 assign(cells[0].vertices, {0, 3, 1, 4});
1010 assign(cells[1].vertices, {1, 4, 2, 5});
1011
1012 tria_back.create_triangulation(
1013 vertices, cells, dealii::SubCellData());
1014
1015 } else {
1016 /* Back part for blunt trailing edge: */
1017
1018 /* Good width for the anisotropically refined center trailing cell: */
1019 double trailing_height =
1020 0.5 / (0.5 + std::pow(2., n_anisotropic_refinements_airfoil_)) *
1021 0.5 * outer_radius;
1022
1023 const std::vector<dealii::Point<2>> vertices{
1024 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 0
1025 {airfoil_center_[0] + back_length,
1026 airfoil_center_[1] + psi_lower(back_length)}, // 1
1027 {airfoil_center_[0] + back_length,
1028 airfoil_center_[1] + psi_upper(back_length)}, // 2
1029 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 3
1030 {outer_radius, -0.5 * outer_radius}, // 4
1031 {outer_radius, -trailing_height}, // 5
1032 {outer_radius, trailing_height}, // 6
1033 {outer_radius, 0.5 * outer_radius}, // 7
1034 };
1035
1036 std::vector<dealii::CellData<2>> cells(3);
1037 assign(cells[0].vertices, {0, 4, 1, 5});
1038 assign(cells[1].vertices, {1, 5, 2, 6});
1039 assign(cells[2].vertices, {2, 6, 3, 7});
1040
1041 tria_back.create_triangulation(
1042 vertices, cells, dealii::SubCellData());
1043 }
1044
1045 dealii::Triangulation<2> coarse_triangulation;
1046 GridGenerator::merge_triangulations(
1047 {&tria_front, &tria_back}, coarse_triangulation, 1.e-12, true);
1048
1049 /*
1050 * Step 3: Set manifold IDs and attach manifolds to preliminary
1051 * coarse triangulation:
1052 *
1053 * Curvature for boundaries:
1054 * 1 -> upper airfoil (inner boundary)
1055 * 2 -> lower airfoil (inner boundary)
1056 * 3 -> spherical manifold (outer boundary)
1057 *
1058 * Transfinite interpolation with grading on coarse cells:
1059 *
1060 * 10 -> bottom center cell
1061 * 11 -> bottom front cell
1062 * 12 -> top front cell
1063 * 13 -> top center cell
1064 * 14 -> bottom trailing cell
1065 * 15 -> top trailing cell (sharp), center trailing cell (blunt)
1066 * 16 -> top trailing cell (blunt)
1067 */
1068
1069 /*
1070 * Colorize boundary faces and add manifolds for curvature
1071 * information on boundaries:
1072 */
1073
1074 for (auto cell : coarse_triangulation.active_cell_iterators()) {
1075 for (auto f : cell->face_indices()) {
1076 const auto face = cell->face(f);
1077 if (!face->at_boundary())
1078 continue;
1079
1080 bool airfoil = true;
1081 bool spherical_boundary = true;
1082 for (const auto v : dealii::GeometryInfo<1>::vertex_indices())
1083 if (std::abs((face->vertex(v)).norm() - outer_radius) < 1.0e-10)
1084 airfoil = false;
1085 else
1086 spherical_boundary = false;
1087
1088 if (spherical_boundary) {
1089 face->set_manifold_id(3);
1090 } else if (airfoil) {
1091 if (face->center()[0] <
1092 airfoil_center_[0] + back_length - 1.e-6) {
1093 if (face->center()[1] >= airfoil_center_[1]) {
1094 face->set_manifold_id(1);
1095 } else {
1096 face->set_manifold_id(2);
1097 }
1098 }
1099 }
1100 } /* f */
1101 } /* cell */
1102
1103 Manifolds::AirfoilManifold airfoil_manifold_upper{
1104 airfoil_center_, psi_front, psi_upper, psi_lower, true, psi_ratio_};
1105 coarse_triangulation.set_manifold(1, airfoil_manifold_upper);
1106
1107 Manifolds::AirfoilManifold airfoil_manifold_lower{airfoil_center_,
1108 psi_front,
1109 psi_upper,
1110 psi_lower,
1111 false,
1112 psi_ratio_};
1113 coarse_triangulation.set_manifold(2, airfoil_manifold_lower);
1114
1115 dealii::SphericalManifold<2> spherical_manifold;
1116 coarse_triangulation.set_manifold(3, spherical_manifold);
1117
1118 /*
1119 * Create transfinite interpolation manifolds for the interior of
1120 * the 2D coarse cells:
1121 */
1122
1123 Assert(!sharp_trailing_edge || (coarse_triangulation.n_cells() == 6),
1124 dealii::ExcInternalError());
1125 Assert(sharp_trailing_edge || (coarse_triangulation.n_cells() == 7),
1126 dealii::ExcInternalError());
1127
1128 std::vector<std::unique_ptr<dealii::Manifold<2, 2>>> manifolds;
1129 manifolds.resize(sharp_trailing_edge ? 6 : 7);
1130
1131 /* FIXME: Remove workaround - mark cells as off limit: */
1132 // WORKAROUND
1133 const auto first_cell = coarse_triangulation.begin_active();
1134 std::next(first_cell, 4)->set_material_id(42);
1135 std::next(first_cell, sharp_trailing_edge ? 5 : 6)->set_material_id(42);
1136 // end WORKAROUND
1137
1138 for (auto i : {0, 1, 2, 3, 5}) {
1139 const auto index = 10 + i;
1140
1141 dealii::Point<2> center;
1142 dealii::Tensor<1, 2> direction;
1143 if (i < 4) {
1144 /* cells: bottom center, bottom front, top front, top center */
1145 direction[1] = 1.;
1146 } else {
1147 Assert(i == 5, dealii::ExcInternalError());
1148 /* cell: center trailing (blunt) */
1149 center[0] = 1.;
1150 direction[0] = -1.;
1151 }
1152
1153 Manifolds::GradingManifold<2> grading{
1154 center,
1155 direction,
1156 grading_,
1157 i == 5 ? grading_epsilon_trailing_ : grading_epsilon_};
1158
1159 auto transfinite =
1160 std::make_unique<TransfiniteInterpolationManifold<2>>();
1161 transfinite->initialize(coarse_triangulation, grading);
1162
1163 coarse_triangulation.set_manifold(index, *transfinite);
1164 manifolds[i] = std::move(transfinite);
1165 }
1166
1167 /* Remove erroneous manifold: */
1168 if (sharp_trailing_edge)
1169 coarse_triangulation.reset_manifold(5);
1170
1171 /*
1172 * Remove unneeded manifolds now. Our custom
1173 * TransfiniteInterpolationManifolds did copy all necessary
1174 * geometry information from the coarse grid already. The boundary
1175 * manifolds are thus not needed any more.
1176 */
1177
1178 coarse_triangulation.reset_manifold(1);
1179 coarse_triangulation.reset_manifold(2);
1180 coarse_triangulation.reset_manifold(3);
1181
1182 /* We can set the final sequence of manifold ids: */
1183 for (unsigned int i = 0; i < (sharp_trailing_edge ? 6 : 7); ++i) {
1184 const auto &cell = std::next(coarse_triangulation.begin_active(), i);
1185 const auto index = 10 + i;
1186 if (i == 4 || i == (sharp_trailing_edge ? 5 : 6)) {
1187 cell->set_manifold_id(index);
1188 } else {
1189 cell->set_all_manifold_ids(index);
1190 }
1191 }
1192
1193 /*
1194 * Attach separate transfinite interpolation manifolds (without a
1195 * grading) to the top and bottom trailing cells:
1196 */
1197
1198 /* FIXME: Remove workaround - mark cells as off limit: */
1199 // WORKAROUND
1200 for (auto cell : coarse_triangulation.active_cell_iterators())
1201 cell->set_material_id(42);
1202 // const auto first_cell = coarse_triangulation.begin_active();
1203 std::next(first_cell, 4)->set_material_id(0);
1204 std::next(first_cell, sharp_trailing_edge ? 5 : 6)->set_material_id(0);
1205 // end WORKAROUND
1206
1207 for (auto i : {4, sharp_trailing_edge ? 5 : 6}) {
1208 const auto index = 10 + i;
1209 auto transfinite =
1210 std::make_unique<ryujin::TransfiniteInterpolationManifold<2>>();
1211 transfinite->initialize(coarse_triangulation);
1212 coarse_triangulation.set_manifold(index, *transfinite);
1213 manifolds[i] = std::move(transfinite);
1214 }
1215
1216 /*
1217 * For good measure, also set material ids. We will need those
1218 * in a minute to reconstruct material ids...
1219 */
1220
1221 for (unsigned int i = 0; i < (sharp_trailing_edge ? 6 : 7); ++i) {
1222 const auto &cell = std::next(coarse_triangulation.begin_active(), i);
1223 const auto index = 10 + i;
1224 cell->set_material_id(index);
1225 }
1226
1227 /*
1228 * Step 4: Anisotropic pre refinement.
1229 *
1230 * Runtime parameters: n_anisotropic_refinements_airfoil_,
1231 * n_anisotropic_refinements_trailing_
1232 */
1233
1234 /* Additional radials in upper and lower cell on airfoil (material
1235 * ids 10 and 13): */
1236 for (unsigned int i = 0; i < n_anisotropic_refinements_airfoil_; ++i) {
1237 for (auto cell : coarse_triangulation.active_cell_iterators()) {
1238 const auto id = cell->material_id();
1239 if (id == 10 || id == 13)
1240 cell->set_refine_flag(dealii::RefinementCase<2>::cut_axis(0));
1241 }
1242
1243 coarse_triangulation.execute_coarsening_and_refinement();
1244 }
1245
1246 /* Anisotropic refinement into trailing cell (material id 15): */
1247 if (!sharp_trailing_edge)
1248 for (unsigned i = 0; i < n_anisotropic_refinements_trailing_; ++i) {
1249 for (auto cell : coarse_triangulation.active_cell_iterators())
1250 if (cell->material_id() == 15)
1251 cell->set_refine_flag(dealii::RefinementCase<2>::cut_axis(0));
1252 else
1253 cell->set_refine_flag();
1254 coarse_triangulation.execute_coarsening_and_refinement();
1255 }
1256
1257 /*
1258 * Step 5: Flatten triangulation, create distributed coarse
1259 * triangulation, and reattach manifolds
1260 *
1261 * Runtime parameters: width_, subdivisions_z_ (for dim == 3)
1262 */
1263
1264 if constexpr (dim == 1) {
1265 AssertThrow(false, dealii::ExcNotImplemented());
1266 __builtin_trap();
1267
1268 } else if constexpr (dim == 2) {
1269 /* Flatten manifold: */
1270 dealii::Triangulation<2> tria3;
1271 tria3.set_mesh_smoothing(triangulation.get_mesh_smoothing());
1272 GridGenerator::flatten_triangulation(coarse_triangulation, tria3);
1273
1274 triangulation.copy_triangulation(tria3);
1275
1276 } else {
1277 static_assert(dim == 3);
1278
1279 /* Flatten manifold: */
1280 dealii::Triangulation<2> tria3;
1281 GridGenerator::flatten_triangulation(coarse_triangulation, tria3);
1282
1283 /* extrude mesh: */
1284 dealii::Triangulation<3, 3> tria4;
1285 tria4.set_mesh_smoothing(triangulation.get_mesh_smoothing());
1286 GridGenerator::extrude_triangulation(
1287 tria3, subdivisions_z_, width_, tria4);
1288
1289 triangulation.copy_triangulation(tria4);
1290 }
1291
1292 /*
1293 * Somewhere during flattening the triangulation, extruding and
1294 * copying, all manifold ids got lost. Reconstruct manifold IDs
1295 * from the material ids we set earlier:
1296 */
1297
1298 for (auto &cell : triangulation.active_cell_iterators()) {
1299 const auto id = cell->material_id();
1300 cell->set_all_manifold_ids(id);
1301 }
1302
1303 /*
1304 * Reattach manifolds:
1305 */
1306 if constexpr (dim == 1) {
1307 AssertThrow(false, dealii::ExcNotImplemented());
1308 __builtin_trap();
1309
1310 } else if constexpr (dim == 2) {
1311 unsigned int index = 10;
1312 for (const auto &manifold : manifolds)
1313 triangulation.set_manifold(index++, *manifold);
1314
1315 } else {
1316 static_assert(dim == 3);
1317
1318 unsigned int index = 10;
1319 for (const auto &manifold : manifolds)
1320 triangulation.set_manifold(
1321 index++, Manifolds::ExtrudedManifold<3>(*manifold));
1322 }
1323
1324 /* Set boundary ids: */
1325
1326 for (auto cell : triangulation.active_cell_iterators()) {
1327 for (auto f : cell->face_indices()) {
1328 auto face = cell->face(f);
1329
1330 /* Handle boundary faces: */
1331 if (!face->at_boundary())
1332 continue;
1333
1334 bool airfoil = true;
1335 bool spherical_boundary = true;
1336
1337 const auto &indices =
1338 dealii::GeometryInfo<dim - 1>::vertex_indices();
1339 for (const auto v : indices) {
1340 const auto vert = face->vertex(v);
1341 const auto radius_sqr = vert[0] * vert[0] + vert[1] * vert[1];
1342 if (radius_sqr >= outer_radius * outer_radius - 1.0e-10 ||
1343 vert[0] > airfoil_center_[0] + 1.001 * back_length)
1344 airfoil = false;
1345 else
1346 spherical_boundary = false;
1347 }
1348
1349 bool periodic_face = (dim == 3);
1350
1351 if constexpr (dim == 3) {
1352 const auto &indices =
1353 dealii::GeometryInfo<dim - 1>::vertex_indices();
1354 bool not_left = false;
1355 bool not_right = false;
1356 for (const auto v : indices) {
1357 const auto vert = face->vertex(v);
1358 if (vert[2] > 1.0e-10)
1359 not_left = true;
1360 if (vert[2] < width_ - 1.0e-10)
1361 not_right = true;
1362 if (not_left && not_right) {
1363 periodic_face = false;
1364 break;
1365 }
1366 }
1367 }
1368
1369 if (periodic_face) {
1370 face->set_boundary_id(Boundary::periodic);
1371 } else if (spherical_boundary) {
1372 face->set_boundary_id(Boundary::dynamic);
1373 } else if (airfoil) {
1374 face->set_boundary_id(Boundary::no_slip);
1375 } else {
1376 Assert(false, dealii::ExcInternalError());
1377 __builtin_trap();
1378 }
1379 }
1380 }
1381
1382 /* Add periodicity: */
1383
1384#ifndef BUG_COLLECT_PERIODIC_FACES_INSTANTIATION
1385 if constexpr (dim == 3) {
1386 std::vector<dealii::GridTools::PeriodicFacePair<
1387 typename dealii::Triangulation<dim>::cell_iterator>>
1388 periodic_faces;
1389
1390 GridTools::collect_periodic_faces(triangulation,
1391 /*b_id */ Boundary::periodic,
1392 /*direction*/ 2,
1393 periodic_faces);
1394
1395 triangulation.add_periodicity(periodic_faces);
1396 }
1397#endif
1398 }
1399
1400 private:
1401 dealii::Point<2> airfoil_center_;
1402 double airfoil_length_;
1403 std::string airfoil_type_;
1404 dealii::Point<2> psi_center_;
1405 double psi_ratio_;
1406 unsigned int psi_samples_;
1407 double height_;
1408 double width_;
1409 double grading_;
1410 double grading_epsilon_;
1411 double grading_epsilon_trailing_;
1412 unsigned int n_anisotropic_refinements_airfoil_;
1413 unsigned int n_anisotropic_refinements_trailing_;
1414 unsigned int subdivisions_z_;
1415 };
1416 } /* namespace Geometries */
1417} /* namespace ryujin */
ryujin::Geometries::Airfoil::Airfoil
Airfoil(const std::string subsection)
Definition: geometry_airfoil.h:826
ryujin::Geometries::Airfoil::create_triangulation
void create_triangulation(dealii::Triangulation< dim > &triangulation) final
Definition: geometry_airfoil.h:906
ryujin::Manifolds::AirfoilManifold::pull_back
dealii::Point< dim > pull_back(const dealii::Point< dim > &space_point) const final
Definition: geometry_airfoil.h:68
ryujin::Manifolds::AirfoilManifold::clone
std::unique_ptr< dealii::Manifold< dim, dim > > clone() const final
Definition: geometry_airfoil.h:123
ryujin::Manifolds::AirfoilManifold::push_forward
dealii::Point< dim > push_forward(const dealii::Point< dim > &point) const final
Definition: geometry_airfoil.h:95
ryujin::Manifolds::AirfoilManifold::AirfoilManifold
AirfoilManifold(const dealii::Point< dim > airfoil_center, const std::function< double(const double)> &psi_front, const std::function< double(const double)> &psi_upper, const std::function< double(const double)> &psi_lower, const bool upper_side, const double psi_ratio=1.0)
Definition: geometry_airfoil.h:47
ryujin::Manifolds::ExtrudedManifold::get_new_point
Point< dim > get_new_point(const ArrayView< const Point< dim > > &surrounding_points, const ArrayView< const double > &weights) const override
Definition: geometry_airfoil.h:251
ryujin::Manifolds::ExtrudedManifold::clone
std::unique_ptr< Manifold< dim > > clone() const override
Definition: geometry_airfoil.h:245
ryujin::Manifolds::ExtrudedManifold::ExtrudedManifold
ExtrudedManifold(const dealii::Manifold< dim - 1 > &manifold)
Definition: geometry_airfoil.h:240
ryujin::Manifolds::GradingManifold::push_forward
dealii::Point< dim > push_forward(const dealii::Point< dim > &chart_point) const final
Definition: geometry_airfoil.h:199
ryujin::Manifolds::GradingManifold::GradingManifold
GradingManifold(const dealii::Point< dim > center, const dealii::Tensor< 1, dim > direction, const double grading, const double epsilon)
Definition: geometry_airfoil.h:154
ryujin::Manifolds::GradingManifold::clone
std::unique_ptr< dealii::Manifold< dim, dim > > clone() const final
Definition: geometry_airfoil.h:219
ryujin::Manifolds::GradingManifold::pull_back
dealii::Point< dim > pull_back(const dealii::Point< dim > &space_point) const final
Definition: geometry_airfoil.h:181
ryujin::Manifolds::GradingManifold::get_new_point
Point< dim > get_new_point(const ArrayView< const Point< dim > > &surrounding_points, const ArrayView< const double > &weights) const override
Definition: geometry_airfoil.h:167
cubic_spline.h
geometry_common_includes.h
ryujin::no_slip
@ no_slip
Definition: discretization.h:90
ryujin::dynamic
@ dynamic
Definition: discretization.h:118
ryujin::periodic
@ periodic
Definition: discretization.h:69
ryujin::pow
T pow(const T x, const T b)
transfinite_interpolation.template.h
Generated by doxygen 1.9.4