ryujin 2.1.1 revision 284c1ee75f38e11ba18808dd0d991d905cf7c040
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(
907 typename Geometry<dim>::Triangulation &triangulation) final
908 {
909 /*
910 * Step 1: Create parametrization:
911 *
912 * Runtime parameters: airfoil_type_, airfoil_length_, psi_center_,
913 * psi_samples_
914 */
915
916 const auto [x_upper, y_upper, x_lower, y_lower] = [&]() {
917 if (airfoil_type_.rfind("NACA ", 0) == 0) {
918 return naca_4digit_points(airfoil_type_.substr(5), psi_samples_);
919 } else if (airfoil_type_.rfind("NASA SC(2) ", 0) == 0) {
920 return nasa_sc2(airfoil_type_.substr(11));
921 } else if (airfoil_type_.rfind("ONERA ", 0) == 0) {
922 return onera(airfoil_type_.substr(6));
923 } else if (airfoil_type_.rfind("BELL ", 0) == 0) {
924 return bell(airfoil_type_.substr(5));
925 }
926 AssertThrow(false, ExcMessage("Unknown airfoil type"));
927 }();
928
929 const auto [psi_front, psi_upper, psi_lower] =
930 create_psi(x_upper,
931 y_upper,
932 x_lower,
933 y_lower,
934 psi_center_[0],
935 psi_center_[1],
936 airfoil_length_);
937
938 /*
939 * Step 2: Create coarse mesh.
940 *
941 * Runtime parameters: airfoil_center_, height_,
942 */
943
944 /* The radius of the radial front part of the mesh: */
945 const auto outer_radius = 0.5 * height_;
946
947 /* by convention, psi_front(0.) returns the "back length" */
948 const auto back_length = psi_front(0.);
949
950 /* sharp trailing edge? */
951 const bool sharp_trailing_edge =
952 std::abs(psi_upper(back_length) - psi_lower(back_length)) < 1.0e-10;
953 AssertThrow(
954 sharp_trailing_edge ||
955 std::abs(psi_upper(back_length) - psi_lower(back_length)) >
956 0.001 * back_length,
957 dealii::ExcMessage("Blunt trailing edge has a width of less than "
958 "0.1% of the trailing airfoil length."));
959
960 /* Front part: */
961 dealii::Triangulation<2> tria_front;
962
963 {
964 std::vector<dealii::Point<2>> vertices{
965 {-outer_radius, 0.0}, // 0
966 {airfoil_center_[0] - psi_front(M_PI), airfoil_center_[1]}, // 1
967 {-0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 2
968 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 3
969 {0., airfoil_center_[1] + psi_lower(-airfoil_center_[0])}, // 4
970 {airfoil_center_[0] + back_length, //
971 airfoil_center_[1] + psi_lower(back_length)}, // 5
972 {0., airfoil_center_[1] + psi_upper(-airfoil_center_[0])}, // 6
973 {-0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 7
974 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 8
975 };
976
977 std::vector<dealii::CellData<2>> cells(4);
978 assign(cells[0].vertices, {2, 3, 4, 5});
979 assign(cells[1].vertices, {0, 2, 1, 4});
980 assign(cells[2].vertices, {7, 0, 6, 1});
981 if (sharp_trailing_edge) {
982 assign(cells[3].vertices, {8, 7, 5, 6});
983 } else {
984 vertices.push_back({airfoil_center_[0] + back_length,
985 airfoil_center_[1] + psi_upper(back_length)});
986 assign(cells[3].vertices, {8, 7, 9, 6});
987 }
988
989 tria_front.create_triangulation(
990 vertices, cells, dealii::SubCellData());
991 }
992
993 /* Back part: */
994 dealii::Triangulation<2> tria_back;
995
996 if (sharp_trailing_edge) {
997 /* Back part for sharp trailing edge: */
998
999 const std::vector<dealii::Point<2>> vertices{
1000 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 0
1001 {airfoil_center_[0] + back_length,
1002 airfoil_center_[1] + psi_lower(back_length)}, // 1
1003 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 2
1004 {outer_radius, -0.5 * outer_radius}, // 3
1005 {outer_radius, 0.0}, // 4
1006 {outer_radius, 0.5 * outer_radius}, // 5
1007 };
1008
1009 std::vector<dealii::CellData<2>> cells(2);
1010 assign(cells[0].vertices, {0, 3, 1, 4});
1011 assign(cells[1].vertices, {1, 4, 2, 5});
1012
1013 tria_back.create_triangulation(
1014 vertices, cells, dealii::SubCellData());
1015
1016 } else {
1017 /* Back part for blunt trailing edge: */
1018
1019 /* Good width for the anisotropically refined center trailing cell: */
1020 double trailing_height =
1021 0.5 / (0.5 + std::pow(2., n_anisotropic_refinements_airfoil_)) *
1022 0.5 * outer_radius;
1023
1024 const std::vector<dealii::Point<2>> vertices{
1025 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius}, // 0
1026 {airfoil_center_[0] + back_length,
1027 airfoil_center_[1] + psi_lower(back_length)}, // 1
1028 {airfoil_center_[0] + back_length,
1029 airfoil_center_[1] + psi_upper(back_length)}, // 2
1030 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius}, // 3
1031 {outer_radius, -0.5 * outer_radius}, // 4
1032 {outer_radius, -trailing_height}, // 5
1033 {outer_radius, trailing_height}, // 6
1034 {outer_radius, 0.5 * outer_radius}, // 7
1035 };
1036
1037 std::vector<dealii::CellData<2>> cells(3);
1038 assign(cells[0].vertices, {0, 4, 1, 5});
1039 assign(cells[1].vertices, {1, 5, 2, 6});
1040 assign(cells[2].vertices, {2, 6, 3, 7});
1041
1042 tria_back.create_triangulation(
1043 vertices, cells, dealii::SubCellData());
1044 }
1045
1046 dealii::Triangulation<2> coarse_triangulation;
1047 GridGenerator::merge_triangulations(
1048 {&tria_front, &tria_back}, coarse_triangulation, 1.e-12, true);
1049
1050 /*
1051 * Step 3: Set manifold IDs and attach manifolds to preliminary
1052 * coarse triangulation:
1053 *
1054 * Curvature for boundaries:
1055 * 1 -> upper airfoil (inner boundary)
1056 * 2 -> lower airfoil (inner boundary)
1057 * 3 -> spherical manifold (outer boundary)
1058 *
1059 * Transfinite interpolation with grading on coarse cells:
1060 *
1061 * 10 -> bottom center cell
1062 * 11 -> bottom front cell
1063 * 12 -> top front cell
1064 * 13 -> top center cell
1065 * 14 -> bottom trailing cell
1066 * 15 -> top trailing cell (sharp), center trailing cell (blunt)
1067 * 16 -> top trailing cell (blunt)
1068 */
1069
1070 /*
1071 * Colorize boundary faces and add manifolds for curvature
1072 * information on boundaries:
1073 */
1074
1075 for (auto cell : coarse_triangulation.active_cell_iterators()) {
1076 for (auto f : dealii::GeometryInfo<2>::face_indices()) {
1077 const auto face = cell->face(f);
1078 if (!face->at_boundary())
1079 continue;
1080
1081 bool airfoil = true;
1082 bool spherical_boundary = true;
1083 for (const auto v : dealii::GeometryInfo<1>::vertex_indices())
1084 if (std::abs((face->vertex(v)).norm() - outer_radius) < 1.0e-10)
1085 airfoil = false;
1086 else
1087 spherical_boundary = false;
1088
1089 if (spherical_boundary) {
1090 face->set_manifold_id(3);
1091 } else if (airfoil) {
1092 if (face->center()[0] <
1093 airfoil_center_[0] + back_length - 1.e-6) {
1094 if (face->center()[1] >= airfoil_center_[1]) {
1095 face->set_manifold_id(1);
1096 } else {
1097 face->set_manifold_id(2);
1098 }
1099 }
1100 }
1101 } /* f */
1102 } /* cell */
1103
1104 Manifolds::AirfoilManifold airfoil_manifold_upper{
1105 airfoil_center_, psi_front, psi_upper, psi_lower, true, psi_ratio_};
1106 coarse_triangulation.set_manifold(1, airfoil_manifold_upper);
1107
1108 Manifolds::AirfoilManifold airfoil_manifold_lower{airfoil_center_,
1109 psi_front,
1110 psi_upper,
1111 psi_lower,
1112 false,
1113 psi_ratio_};
1114 coarse_triangulation.set_manifold(2, airfoil_manifold_lower);
1115
1116 dealii::SphericalManifold<2> spherical_manifold;
1117 coarse_triangulation.set_manifold(3, spherical_manifold);
1118
1119 /*
1120 * Create transfinite interpolation manifolds for the interior of
1121 * the 2D coarse cells:
1122 */
1123
1124 Assert(!sharp_trailing_edge || (coarse_triangulation.n_cells() == 6),
1125 dealii::ExcInternalError());
1126 Assert(sharp_trailing_edge || (coarse_triangulation.n_cells() == 7),
1127 dealii::ExcInternalError());
1128
1129 std::vector<std::unique_ptr<dealii::Manifold<2, 2>>> manifolds;
1130 manifolds.resize(sharp_trailing_edge ? 6 : 7);
1131
1132 /* FIXME: Remove workaround - mark cells as off limit: */
1133 // WORKAROUND
1134 const auto first_cell = coarse_triangulation.begin_active();
1135 std::next(first_cell, 4)->set_material_id(42);
1136 std::next(first_cell, sharp_trailing_edge ? 5 : 6)->set_material_id(42);
1137 // end WORKAROUND
1138
1139 for (auto i : {0, 1, 2, 3, 5}) {
1140 const auto index = 10 + i;
1141
1142 dealii::Point<2> center;
1143 dealii::Tensor<1, 2> direction;
1144 if (i < 4) {
1145 /* cells: bottom center, bottom front, top front, top center */
1146 direction[1] = 1.;
1147 } else {
1148 Assert(i == 5, dealii::ExcInternalError());
1149 /* cell: center trailing (blunt) */
1150 center[0] = 1.;
1151 direction[0] = -1.;
1152 }
1153
1154 Manifolds::GradingManifold<2> grading{
1155 center,
1156 direction,
1157 grading_,
1158 i == 5 ? grading_epsilon_trailing_ : grading_epsilon_};
1159
1160 auto transfinite =
1161 std::make_unique<TransfiniteInterpolationManifold<2>>();
1162 transfinite->initialize(coarse_triangulation, grading);
1163
1164 coarse_triangulation.set_manifold(index, *transfinite);
1165 manifolds[i] = std::move(transfinite);
1166 }
1167
1168 /* Remove erroneous manifold: */
1169 if (sharp_trailing_edge)
1170 coarse_triangulation.reset_manifold(5);
1171
1172 /*
1173 * Remove unneeded manifolds now. Our custom
1174 * TransfiniteInterpolationManifolds did copy all necessary
1175 * geometry information from the coarse grid already. The boundary
1176 * manifolds are thus not needed any more.
1177 */
1178
1179 coarse_triangulation.reset_manifold(1);
1180 coarse_triangulation.reset_manifold(2);
1181 coarse_triangulation.reset_manifold(3);
1182
1183 /* We can set the final sequence of manifold ids: */
1184 for (unsigned int i = 0; i < (sharp_trailing_edge ? 6 : 7); ++i) {
1185 const auto &cell = std::next(coarse_triangulation.begin_active(), i);
1186 const auto index = 10 + i;
1187 if (i == 4 || i == (sharp_trailing_edge ? 5 : 6)) {
1188 cell->set_manifold_id(index);
1189 } else {
1190 cell->set_all_manifold_ids(index);
1191 }
1192 }
1193
1194 /*
1195 * Attach separate transfinite interpolation manifolds (without a
1196 * grading) to the top and bottom trailing cells:
1197 */
1198
1199 /* FIXME: Remove workaround - mark cells as off limit: */
1200 // WORKAROUND
1201 for (auto cell : coarse_triangulation.active_cell_iterators())
1202 cell->set_material_id(42);
1203 // const auto first_cell = coarse_triangulation.begin_active();
1204 std::next(first_cell, 4)->set_material_id(0);
1205 std::next(first_cell, sharp_trailing_edge ? 5 : 6)->set_material_id(0);
1206 // end WORKAROUND
1207
1208 for (auto i : {4, sharp_trailing_edge ? 5 : 6}) {
1209 const auto index = 10 + i;
1210 auto transfinite =
1211 std::make_unique<ryujin::TransfiniteInterpolationManifold<2>>();
1212 transfinite->initialize(coarse_triangulation);
1213 coarse_triangulation.set_manifold(index, *transfinite);
1214 manifolds[i] = std::move(transfinite);
1215 }
1216
1217 /*
1218 * For good measure, also set material ids. We will need those
1219 * in a minute to reconstruct material ids...
1220 */
1221
1222 for (unsigned int i = 0; i < (sharp_trailing_edge ? 6 : 7); ++i) {
1223 const auto &cell = std::next(coarse_triangulation.begin_active(), i);
1224 const auto index = 10 + i;
1225 cell->set_material_id(index);
1226 }
1227
1228 /*
1229 * Step 4: Anisotropic pre refinement.
1230 *
1231 * Runtime parameters: n_anisotropic_refinements_airfoil_,
1232 * n_anisotropic_refinements_trailing_
1233 */
1234
1235 /* Additional radials in upper and lower cell on airfoil (material
1236 * ids 10 and 13): */
1237 for (unsigned int i = 0; i < n_anisotropic_refinements_airfoil_; ++i) {
1238 for (auto cell : coarse_triangulation.active_cell_iterators()) {
1239 const auto id = cell->material_id();
1240 if (id == 10 || id == 13)
1241 cell->set_refine_flag(dealii::RefinementCase<2>::cut_axis(0));
1242 }
1243
1244 coarse_triangulation.execute_coarsening_and_refinement();
1245 }
1246
1247 /* Anisotropic refinement into trailing cell (material id 15): */
1248 if (!sharp_trailing_edge)
1249 for (unsigned i = 0; i < n_anisotropic_refinements_trailing_; ++i) {
1250 for (auto cell : coarse_triangulation.active_cell_iterators())
1251 if (cell->material_id() == 15)
1252 cell->set_refine_flag(dealii::RefinementCase<2>::cut_axis(0));
1253 else
1254 cell->set_refine_flag();
1255 coarse_triangulation.execute_coarsening_and_refinement();
1256 }
1257
1258 /*
1259 * Step 5: Flatten triangulation, create distributed coarse
1260 * triangulation, and reattach manifolds
1261 *
1262 * Runtime parameters: width_, subdivisions_z_ (for dim == 3)
1263 */
1264
1265 if constexpr (dim == 1) {
1266 AssertThrow(false, dealii::ExcNotImplemented());
1267 __builtin_trap();
1268
1269 } else if constexpr (dim == 2) {
1270 /* Flatten manifold: */
1271 dealii::Triangulation<2> tria3;
1272 tria3.set_mesh_smoothing(triangulation.get_mesh_smoothing());
1273 GridGenerator::flatten_triangulation(coarse_triangulation, tria3);
1274
1275 triangulation.copy_triangulation(tria3);
1276
1277 } else {
1278 static_assert(dim == 3);
1279
1280 /* Flatten manifold: */
1281 dealii::Triangulation<2> tria3;
1282 GridGenerator::flatten_triangulation(coarse_triangulation, tria3);
1283
1284 /* extrude mesh: */
1285 dealii::Triangulation<3, 3> tria4;
1286 tria4.set_mesh_smoothing(triangulation.get_mesh_smoothing());
1287 GridGenerator::extrude_triangulation(
1288 tria3, subdivisions_z_, width_, tria4);
1289
1290 triangulation.copy_triangulation(tria4);
1291 }
1292
1293 /*
1294 * Somewhere during flattening the triangulation, extruding and
1295 * copying, all manifold ids got lost. Reconstruct manifold IDs
1296 * from the material ids we set earlier:
1297 */
1298
1299 for (auto &cell : triangulation.active_cell_iterators()) {
1300 const auto id = cell->material_id();
1301 cell->set_all_manifold_ids(id);
1302 }
1303
1304 /*
1305 * Reattach manifolds:
1306 */
1307 if constexpr (dim == 1) {
1308 AssertThrow(false, dealii::ExcNotImplemented());
1309 __builtin_trap();
1310
1311 } else if constexpr (dim == 2) {
1312 unsigned int index = 10;
1313 for (const auto &manifold : manifolds)
1314 triangulation.set_manifold(index++, *manifold);
1315
1316 } else {
1317 static_assert(dim == 3);
1318
1319 unsigned int index = 10;
1320 for (const auto &manifold : manifolds)
1321 triangulation.set_manifold(
1322 index++, Manifolds::ExtrudedManifold<3>(*manifold));
1323 }
1324
1325 /* Set boundary ids: */
1326
1327 for (auto cell : triangulation.active_cell_iterators()) {
1328 for (auto f : dealii::GeometryInfo<dim>::face_indices()) {
1329 auto face = cell->face(f);
1330
1331 /* Handle boundary faces: */
1332 if (!face->at_boundary())
1333 continue;
1334
1335 bool airfoil = true;
1336 bool spherical_boundary = true;
1337
1338 const auto &indices =
1339 dealii::GeometryInfo<dim - 1>::vertex_indices();
1340 for (const auto v : indices) {
1341 const auto vert = face->vertex(v);
1342 const auto radius_sqr = vert[0] * vert[0] + vert[1] * vert[1];
1343 if (radius_sqr >= outer_radius * outer_radius - 1.0e-10 ||
1344 vert[0] > airfoil_center_[0] + 1.001 * back_length)
1345 airfoil = false;
1346 else
1347 spherical_boundary = false;
1348 }
1349
1350 bool periodic_face = (dim == 3);
1351
1352 if constexpr (dim == 3) {
1353 const auto &indices =
1354 dealii::GeometryInfo<dim - 1>::vertex_indices();
1355 bool not_left = false;
1356 bool not_right = false;
1357 for (const auto v : indices) {
1358 const auto vert = face->vertex(v);
1359 if (vert[2] > 1.0e-10)
1360 not_left = true;
1361 if (vert[2] < width_ - 1.0e-10)
1362 not_right = true;
1363 if (not_left && not_right) {
1364 periodic_face = false;
1365 break;
1366 }
1367 }
1368 }
1369
1370 if (periodic_face) {
1371 face->set_boundary_id(Boundary::periodic);
1372 } else if (spherical_boundary) {
1373 face->set_boundary_id(Boundary::dynamic);
1374 } else if (airfoil) {
1375 face->set_boundary_id(Boundary::no_slip);
1376 } else {
1377 Assert(false, dealii::ExcInternalError());
1378 __builtin_trap();
1379 }
1380 }
1381 }
1382
1383 /* Add periodicity: */
1384
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 }
1398
1399 private:
1400 dealii::Point<2> airfoil_center_;
1401 double airfoil_length_;
1402 std::string airfoil_type_;
1403 dealii::Point<2> psi_center_;
1404 double psi_ratio_;
1405 unsigned int psi_samples_;
1406 double height_;
1407 double width_;
1408 double grading_;
1409 double grading_epsilon_;
1410 double grading_epsilon_trailing_;
1411 unsigned int n_anisotropic_refinements_airfoil_;
1412 unsigned int n_anisotropic_refinements_trailing_;
1413 unsigned int subdivisions_z_;
1414 };
1415 } /* namespace Geometries */
1416} /* namespace ryujin */
ryujin::Geometries::Airfoil::Airfoil
Airfoil(const std::string subsection)
Definition: geometry_airfoil.h:826
ryujin::Geometries::Airfoil::create_triangulation
void create_triangulation(typename Geometry< dim >::Triangulation &triangulation) final
Definition: geometry_airfoil.h:906
ryujin::Geometry::Triangulation
typename Discretization< dim >::Triangulation Triangulation
Definition: geometry.h:38
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:85
ryujin::dynamic
@ dynamic
Definition: discretization.h:105
ryujin::periodic
@ periodic
Definition: discretization.h:64
ryujin::pow
T pow(const T x, const T b)
transfinite_interpolation.template.h
Generated by doxygen 1.9.4