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