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