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