23 template <
typename T,
typename T2>
24 inline void assign(T &array,
25 const std::initializer_list<T2> &initializer_list)
28 Assert(std::size(array) == std::size(initializer_list),
30 "size of initializer list and array does not match"));
32 initializer_list.begin(), initializer_list.end(), std::begin(array));
44 static_assert(dim == 2,
"not implemented for dim != 2");
48 const std::function<
double(
const double)> &psi_front,
49 const std::function<
double(
const double)> &psi_upper,
50 const std::function<
double(
const double)> &psi_lower,
51 const bool upper_side,
52 const double psi_ratio = 1.0)
53 : airfoil_center(airfoil_center)
54 , psi_front(psi_front)
55 , psi_upper(psi_upper)
56 , psi_lower(psi_lower)
57 , upper_side(upper_side)
58 , ratio_(psi_ratio * psi_front(0.) / psi_front(M_PI))
61 Assert(std::abs(psi_upper(0.) - psi_front(0.5 * M_PI)) < 1.0e-10,
62 dealii::ExcInternalError());
63 Assert(std::abs(psi_lower(0.) + psi_front(1.5 * M_PI)) < 1.0e-10,
64 dealii::ExcInternalError());
68 pull_back(
const dealii::Point<dim> &space_point)
const final
70 auto coordinate = dealii::Point<dim>() + (space_point - airfoil_center);
74 dealii::Point<dim> chart_point;
75 if (coordinate[0] > 0.) {
78 chart_point[0] = 1. + coordinate[1] - psi_upper(coordinate[0]);
79 chart_point[1] = 0.5 * M_PI - ratio_ * coordinate[0];
82 chart_point[0] = 1. - coordinate[1] + psi_lower(coordinate[0]);
83 chart_point[1] = 1.5 * M_PI + ratio_ * coordinate[0];
87 chart_point = polar_manifold.pull_back(coordinate);
88 chart_point[0] = 1. + chart_point[0] - psi_front(chart_point[1]);
97 auto chart_point = point;
101 dealii::Point<dim> coordinate;
102 if (chart_point[1] < 0.5 * M_PI) {
103 Assert(upper_side, dealii::ExcInternalError());
105 coordinate[0] = (0.5 * M_PI - chart_point[1]) / ratio_;
106 Assert(coordinate[0] >= -1.0e-10, dealii::ExcInternalError());
107 coordinate[1] = chart_point[0] - 1. + psi_upper(coordinate[0]);
108 }
else if (chart_point[1] > 1.5 * M_PI) {
109 Assert(!upper_side, dealii::ExcInternalError());
111 coordinate[0] = (chart_point[1] - 1.5 * M_PI) / ratio_;
112 Assert(coordinate[0] >= -1.0e-10, dealii::ExcInternalError());
113 coordinate[1] = 1. - chart_point[0] + psi_lower(coordinate[0]);
116 chart_point[0] = chart_point[0] - 1. + psi_front(chart_point[1]);
117 coordinate = polar_manifold.push_forward(chart_point);
120 return dealii::Point<dim>() + (coordinate + airfoil_center);
123 std::unique_ptr<dealii::Manifold<dim, dim>>
clone() const final
125 const double psi_ratio = ratio_ * psi_front(M_PI) / psi_front(0.);
126 return std::make_unique<AirfoilManifold<dim>>(airfoil_center,
135 const dealii::Point<dim> airfoil_center;
136 const std::function<double(
const double)> psi_front;
137 const std::function<double(
const double)> psi_upper;
138 const std::function<double(
const double)> psi_lower;
139 const bool upper_side;
143 dealii::PolarManifold<dim> polar_manifold;
155 const dealii::Tensor<1, dim> direction,
156 const double grading,
157 const double epsilon)
159 , direction(direction)
168 const ArrayView<const double> &weights)
const override
170 if (weights[0] > 1.0)
171 return surrounding_points[0];
173 if (weights[1] > 1.0)
174 return surrounding_points[1];
176 return dealii::ChartManifold<dim>::get_new_point(surrounding_points,
181 pull_back(
const dealii::Point<dim> &space_point)
const final
183 auto chart_point = space_point - center;
185 for (
unsigned int d = 0; d < dim; ++d) {
186 if (std::abs(direction[d]) > 1.0e-10) {
187 const double x = chart_point[d] * std::copysign(1., direction[d]);
188 Assert(x + epsilon > 0, dealii::ExcInternalError());
189 const double x_hat =
std::pow(x + epsilon, 1. / grading) -
190 std::pow(epsilon, 1. / grading) + 1.e-14;
191 chart_point[d] += (x_hat - x) * std::copysign(1., direction[d]);
195 return dealii::Point<dim>() + chart_point;
201 auto space_point = chart_point;
203 for (
unsigned int d = 0; d < dim; ++d) {
204 if (std::abs(direction[d]) > 1.0e-10) {
206 space_point[d] * std::copysign(1., direction[d]);
207 Assert(x_hat +
std::pow(epsilon, 1. / grading) > 0,
208 dealii::ExcInternalError());
212 space_point[d] += (x - x_hat) * std::copysign(1., direction[d]);
216 return center + (space_point - dealii::Point<dim>());
219 std::unique_ptr<dealii::Manifold<dim, dim>>
clone() const final
221 return std::make_unique<GradingManifold<dim>>(
222 center, direction, grading, epsilon);
226 const dealii::Point<dim> center;
227 const dealii::Tensor<1, dim> direction;
228 const double grading;
229 const double epsilon;
241 : manifold(manifold.
clone())
245 std::unique_ptr<Manifold<dim>>
clone()
const override
247 return std::make_unique<ExtrudedManifold<dim>>(*manifold);
252 const ArrayView<const double> &weights)
const override
254 Assert(surrounding_points.size() == weights.size(),
255 dealii::ExcInternalError());
257 boost::container::small_vector<dealii::Point<dim - 1>, 100>
258 surrounding_points_projected;
259 std::transform(surrounding_points.begin(),
260 surrounding_points.end(),
261 surrounding_points_projected.begin(),
262 [](
const dealii::Point<dim> &source) {
263 dealii::Point<dim - 1> result;
264 for (unsigned int d = 0; d < dim - 1; ++d)
265 result[d] = source[d];
269 const auto projected = manifold->get_new_point(
270 ArrayView<
const Point<dim - 1>>{surrounding_points_projected.data(),
274 dealii::Point<dim> result;
276 for (
unsigned int d = 0; d < dim - 1; ++d)
277 result[d] = projected[d];
279 for (
unsigned int i = 0; i < weights.size(); ++i)
280 result[dim - 1] += weights[i] * surrounding_points[i][dim - 1];
286 std::unique_ptr<
const dealii::Manifold<dim - 1>> manifold;
297 std::array<std::vector<double>, 4>
298 naca_4digit_points(
const std::string &serial_number,
299 const unsigned int n_samples)
301 AssertThrow(serial_number.size() == 4,
302 dealii::ExcMessage(
"Invalid NACA 4 digit serial number"));
303 std::array<unsigned int, 4> digit;
304 std::transform(serial_number.begin(),
307 [](
auto it) { return it -
'0'; });
310 const double t = 0.1 * digit[2] + 0.01 * digit[3];
312 dealii::ExcMessage(
"Invalid NACA 4 digit serial number"));
315 const double m = 0.01 * digit[0];
317 const double p = 0.1 * digit[1];
319 std::vector<double> x_upper;
320 std::vector<double> y_upper;
321 std::vector<double> x_lower;
322 std::vector<double> y_lower;
324 for (
unsigned int i = 0; i < n_samples; i++) {
325 const double x = 1. * i / (n_samples - 1);
328 (0.2969 * std::sqrt(x) +
329 x * (-0.126 + x * (-0.3516 + x * (0.2843 + x * (-0.1036)))));
331 const double y_c = (x < p) ? m / (p * p) * (2. * p * x - x * x)
332 : m / ((1. - p) * (1. - p)) *
333 (1. - 2. * p + 2. * p * x - x * x);
335 const double dy_c = (x < p) ? 2. * m / (p * p) * (p - x)
336 : 2. * m / ((1. - p) * (1. - p)) * (p - x);
338 const double theta = std::atan(dy_c);
340 x_upper.emplace_back(x - y * std::sin(theta));
341 y_upper.emplace_back(y_c + y * std::cos(theta));
342 x_lower.emplace_back(x + y * std::sin(theta));
343 y_lower.emplace_back(y_c - y * std::cos(theta));
347 y_upper.front() = 0.;
349 y_lower.front() = 0.;
352 return {{x_upper, y_upper, x_lower, y_lower}};
365 std::array<std::vector<double>, 4>
366 nasa_sc2(
const std::string &serial_number)
368 if (serial_number ==
"0714") {
369 std::vector<double> x_upper{
370 .0, .002, .005, .01, .02, .03, .04, .05, .07, .1, .12, .15,
371 .17, .2, .22, .25, .27, .3, .33, .35, .38, .4, .43, .45,
372 .48, .50, .53, .55, .57, .6, .62, .65, .68, .7, .72, .75,
373 .77, .8, .82, .85, .87, .9, .92, .95, .97, .98, .99, 1.};
375 std::vector<double> y_upper{
376 .0, .0095, .0158, .0219, .0293, .0343, .0381, .0411,
377 .0462, .0518, .0548, .0585, .0606, .0632, .0646, .0664,
378 .0673, .0685, .0692, .0696, .0698, .0697, .0695, .0692,
379 .0684, .0678, .0666, .0656, .0645, .0625, .0610, .0585,
380 .0555, .0533, .0509, .0469, .0439, .0389, .0353, .0294,
381 .0251, .0181, .0131, .0049, -.0009, -.0039, -.0071, -.0104};
383 std::vector<double> x_lower{
384 .0, .002, .005, .01, .02, .03, .04, .05, .07, .1, .12, .15, .17,
385 .20, .22, .25, .28, .3, .32, .35, .37, .4, .42, .45, .48, .5,
386 .53, .55, .58, .6, .63, .65, .68, .70, .73, .75, .77, .80, .83,
387 .85, .87, .89, .92, .94, .95, .96, .97, .98, .99, 1.};
389 std::vector<double> y_lower{
390 .0, -.0093, -.016, -.0221, -.0295, -.0344, -.0381, -.0412,
391 -.0462, -.0517, -.0547, -.0585, -.0606, -.0633, -.0647, -.0666,
392 -.068, -.0687, -.0692, -.0696, -.0696, -.0692, -.0688, -.0676,
393 -.0657, -.0644, -.0614, -.0588, -.0543, -.0509, -.0451, -.041,
394 -.0346, -.0302, -.0235, -.0192, -.0150, -.0093, -.0048, -.0024,
395 -.0013, -.0008, -.0016, -.0035, -.0049, -.0066, -.0085, -.0109,
398 return {{x_upper, y_upper, x_lower, y_lower}};
403 dealii::ExcMessage(
"Invalid NASA SC(2) serial number"));
416 std::array<std::vector<double>, 4> onera(
const std::string &serial_number)
418 if (serial_number ==
"OAT15a") {
419 std::vector<double> x_upper{
420 0., 2.95888e-05, 0.000117865, 0.000263239, 0.000464245,
421 0.000719821, 0.00103013, 0.00139605, 0.00181884, 0.00230024,
422 0.00284243, 0.00344764, 0.00411805, 0.00485595, 0.00566349,
423 0.00654241, 0.00749421, 0.0085197, 0.0096197, 0.0107945,
424 0.0120442, 0.013369, 0.0147688, 0.0162438, 0.0177939,
425 0.0194194, 0.0211207, 0.0228982, 0.0247526, 0.0266845,
426 0.028695, 0.0307852, 0.0329562, 0.0352094, 0.0375463,
427 0.0399687, 0.0424782, 0.0450769, 0.0477668, 0.05055,
428 0.0534291, 0.0564063, 0.0594842, 0.0626655, 0.0659531,
429 0.0693498, 0.0728588, 0.0764831, 0.0802261, 0.0840914,
430 0.0880824, 0.0922027, 0.0964564, 0.100847, 0.10538,
431 0.110058, 0.114885, 0.119868, 0.125009, 0.130314,
432 0.135789, 0.14139, 0.147074, 0.152839, 0.158682,
433 0.164603, 0.170599, 0.17667, 0.182814, 0.189028,
434 0.195312, 0.201665, 0.208083, 0.214567, 0.221115,
435 0.227724, 0.234394, 0.241123, 0.24791, 0.254753,
436 0.26165, 0.2686, 0.275601, 0.282653, 0.289753,
437 0.2969, 0.304093, 0.311329, 0.318609, 0.325929,
438 0.333289, 0.340686, 0.348121, 0.35559, 0.363093,
439 0.370629, 0.378194, 0.385789, 0.393412, 0.40106,
440 0.408734, 0.41643, 0.424148, 0.431886, 0.439643,
441 0.447417, 0.455207, 0.46301, 0.470827, 0.478654,
442 0.486491, 0.494337, 0.502189, 0.510046, 0.517907,
443 0.525769, 0.533632, 0.541494, 0.549354, 0.557209,
444 0.565059, 0.572902, 0.580736, 0.588559, 0.596371,
445 0.60417, 0.611953, 0.61972, 0.627469, 0.635198,
446 0.642906, 0.65059, 0.658251, 0.665886, 0.673493,
447 0.68107, 0.688617, 0.69613, 0.70361, 0.711054,
448 0.71846, 0.725827, 0.733154, 0.740438, 0.747679,
449 0.754875, 0.762025, 0.769127, 0.776181, 0.783185,
450 0.790139, 0.79704, 0.80389, 0.810685, 0.817426,
451 0.82411, 0.830738, 0.837307, 0.843817, 0.850265,
452 0.856652, 0.862974, 0.869233, 0.87538, 0.881373,
453 0.887216, 0.892913, 0.898467, 0.903883, 0.909163,
454 0.914311, 0.91933, 0.924224, 0.928996, 0.933648,
455 0.938183, 0.942606, 0.946917, 0.95112, 0.955217,
456 0.959212, 0.963107, 0.966904, 0.970605, 0.974213,
457 0.977731, 0.98116, 0.984503, 0.987762, 0.990939,
458 0.994036, 0.997056, 1.};
460 std::vector<double> y_upper{
461 0., 0.000899353, 0.0018231, 0.00276894, 0.00373508,
462 0.00472011, 0.0057226, 0.00674103, 0.0077738, 0.00881906,
463 0.00987467, 0.0109383, 0.0120074, 0.0130793, 0.0141511,
464 0.01522, 0.0162832, 0.0173387, 0.0183841, 0.0194179,
465 0.0204389, 0.021446, 0.0224386, 0.0234164, 0.0243794,
466 0.0253276, 0.0262612, 0.0271805, 0.028086, 0.0289783,
467 0.0298578, 0.0307252, 0.0315811, 0.0324262, 0.0332611,
468 0.0340861, 0.0349022, 0.0357098, 0.0365095, 0.0373016,
469 0.0380867, 0.0388652, 0.0396375, 0.0404039, 0.041165,
470 0.0419211, 0.0426722, 0.0434189, 0.0441614, 0.0448995,
471 0.0456336, 0.0463636, 0.0470894, 0.047811, 0.0485286,
472 0.0492423, 0.0499518, 0.0506574, 0.0513591, 0.0520569,
473 0.0527506, 0.0534343, 0.0541023, 0.054755, 0.0553921,
474 0.0560138, 0.05662, 0.0572108, 0.0577861, 0.0583462,
475 0.0588909, 0.0594202, 0.0599341, 0.0604325, 0.0609153,
476 0.0613826, 0.0618341, 0.06227, 0.06269, 0.0630941,
477 0.0634823, 0.0638544, 0.0642103, 0.06455, 0.0648734,
478 0.0651806, 0.0654713, 0.0657454, 0.0660031, 0.0662442,
479 0.0664685, 0.066676, 0.0668664, 0.0670398, 0.067196,
480 0.0673349, 0.0674562, 0.0675598, 0.0676456, 0.0677134,
481 0.0677629, 0.067794, 0.0678065, 0.0678, 0.0677743,
482 0.0677293, 0.0676646, 0.0675798, 0.0674748, 0.0673492,
483 0.0672027, 0.0670349, 0.0668456, 0.0666344, 0.066401,
484 0.066145, 0.0658661, 0.065564, 0.0652385, 0.064889,
485 0.0645151, 0.0641169, 0.0636938, 0.0632454, 0.0627715,
486 0.0622718, 0.061746, 0.0611937, 0.0606145, 0.0600083,
487 0.0593747, 0.0587136, 0.0580244, 0.0573069, 0.0565607,
488 0.0557853, 0.0549807, 0.0541461, 0.0532814, 0.0523863,
489 0.0514606, 0.0505046, 0.0495188, 0.0485042, 0.047462,
490 0.0463943, 0.0453031, 0.0441914, 0.0430618, 0.0419174,
491 0.0407612, 0.0395961, 0.038425, 0.0372503, 0.0360742,
492 0.034899, 0.0337262, 0.0325572, 0.0313935, 0.030236,
493 0.029086, 0.0279442, 0.0268114, 0.0256966, 0.0246079,
494 0.023545, 0.0225073, 0.0214947, 0.0205065, 0.0195422,
495 0.0186011, 0.0176826, 0.0167863, 0.0159112, 0.0150567,
496 0.0142221, 0.0134066, 0.0126095, 0.0118301, 0.0110678,
497 0.0103219, 0.00959177, 0.00887695, 0.00817697, 0.00749135,
498 0.00681977, 0.0061621, 0.00551806, 0.00488739, 0.00427007,
499 0.00366612, 0.00307588, 0.0024997};
501 std::vector<double> x_lower{
502 0., 3.22311e-05, 0.000136327, 0.000324365, 0.000606007,
503 0.000986654, 0.00146626, 0.00204126, 0.00270571, 0.00345312,
504 0.00427708, 0.00517234, 0.00613428, 0.00715943, 0.00824517,
505 0.00938973, 0.0105917, 0.0118504, 0.0131652, 0.0145362,
506 0.0159635, 0.0174476, 0.0189891, 0.0205889, 0.0222479,
507 0.0239675, 0.0257488, 0.0275934, 0.0295028, 0.031479,
508 0.0335237, 0.035639, 0.0378271, 0.04009, 0.0424303,
509 0.0448502, 0.0473523, 0.049939, 0.0526132, 0.0553774,
510 0.0582343, 0.0611868, 0.0642377, 0.0673899, 0.0706465,
511 0.0740105, 0.077485, 0.0810733, 0.0847788, 0.088605,
512 0.0925553, 0.0966336, 0.100844, 0.10519, 0.109675,
513 0.114305, 0.119084, 0.124016, 0.129106, 0.134358,
514 0.139779, 0.145328, 0.150961, 0.156676, 0.162473,
515 0.168349, 0.174303, 0.180333, 0.186437, 0.192615,
516 0.198863, 0.205181, 0.211568, 0.21802, 0.224537,
517 0.231117, 0.237759, 0.24446, 0.251219, 0.258035,
518 0.264905, 0.271829, 0.278804, 0.285828, 0.292901,
519 0.30002, 0.307184, 0.314391, 0.321639, 0.328927,
520 0.336253, 0.343615, 0.351011, 0.358441, 0.365902,
521 0.373392, 0.38091, 0.388455, 0.396024, 0.403616,
522 0.41123, 0.418864, 0.426517, 0.434187, 0.441871,
523 0.44957, 0.457282, 0.465005, 0.472737, 0.480477,
524 0.488225, 0.495977, 0.503733, 0.511493, 0.519252,
525 0.527012, 0.53477, 0.542525, 0.550276, 0.55802,
526 0.565759, 0.573488, 0.581208, 0.588918, 0.596615,
527 0.604298, 0.611967, 0.61962, 0.627257, 0.634874,
528 0.642473, 0.65005, 0.657605, 0.665137, 0.672645,
529 0.680128, 0.687584, 0.695012, 0.702411, 0.70978,
530 0.717118, 0.724424, 0.731697, 0.738935, 0.746137,
531 0.753303, 0.76043, 0.767518, 0.774565, 0.78157,
532 0.788531, 0.795447, 0.802316, 0.809136, 0.815905,
533 0.822623, 0.829286, 0.835893, 0.842441, 0.84893,
534 0.855357, 0.86172, 0.868018, 0.874204, 0.880239,
535 0.886123, 0.891863, 0.89746, 0.902919, 0.908242,
536 0.913433, 0.918495, 0.923431, 0.928244, 0.932938,
537 0.937515, 0.941978, 0.94633, 0.950574, 0.954712,
538 0.958747, 0.962681, 0.966518, 0.970259, 0.973907,
539 0.977464, 0.980932, 0.984314, 0.987612, 0.990827,
540 0.993963, 0.997019, 1.};
542 std::vector<double> y_lower{
543 0., -0.000899234, -0.00182108, -0.00275889, -0.00370397,
544 -0.00464681, -0.00557909, -0.00649452, -0.0073895, -0.00826265,
545 -0.00911444, -0.00994611, -0.0107598, -0.0115578, -0.0123422,
546 -0.0131152, -0.0138788, -0.0146348, -0.0153851, -0.0161313,
547 -0.0168749, -0.0176173, -0.0183597, -0.0191032, -0.0198489,
548 -0.0205974, -0.0213498, -0.0221064, -0.0228678, -0.0236341,
549 -0.0244053, -0.0251817, -0.0259627, -0.0267482, -0.0275378,
550 -0.0283309, -0.029127, -0.0299255, -0.0307258, -0.0315273,
551 -0.0323295, -0.0331322, -0.0339346, -0.0347366, -0.035538,
552 -0.0363383, -0.0371378, -0.037936, -0.0387331, -0.0395288,
553 -0.0403228, -0.0411152, -0.0419054, -0.0426933, -0.0434778,
554 -0.0442587, -0.0450347, -0.045805, -0.0465681, -0.0473225,
555 -0.0480666, -0.0487929, -0.0494941, -0.0501695, -0.0508179,
556 -0.0514387, -0.052031, -0.0525942, -0.0531277, -0.0536309,
557 -0.0541035, -0.054545, -0.0549547, -0.0553323, -0.0556773,
558 -0.0559891, -0.056267, -0.0565107, -0.0567191, -0.0568918,
559 -0.0570281, -0.0571272, -0.0571886, -0.0572116, -0.0571958,
560 -0.0571405, -0.0570453, -0.0569098, -0.0567339, -0.0565171,
561 -0.0562594, -0.0559608, -0.055621, -0.0552405, -0.0548193,
562 -0.0543577, -0.0538564, -0.0533157, -0.0527365, -0.0521197,
563 -0.0514661, -0.0507767, -0.0500525, -0.0492947, -0.0485041,
564 -0.0476819, -0.0468294, -0.0459474, -0.0450371, -0.0440995,
565 -0.0431356, -0.0421465, -0.0411333, -0.0400971, -0.039039,
566 -0.0379601, -0.0368617, -0.0357451, -0.0346114, -0.0334621,
567 -0.0322983, -0.0311217, -0.0299336, -0.0287354, -0.0275285,
568 -0.0263145, -0.025095, -0.0238717, -0.0226459, -0.0214194,
569 -0.020194, -0.0189714, -0.0177534, -0.016542, -0.015339,
570 -0.0141466, -0.0129671, -0.0118026, -0.0106558, -0.00952898,
571 -0.00842491, -0.00734634, -0.00629613, -0.00527714, -0.0042922,
572 -0.00334424, -0.00243619, -0.00157086, -0.000750868, 2.13187e-05,
573 0.000743489, 0.00141373, 0.00203031, 0.00259174, 0.00309711,
574 0.00354599, 0.00393819, 0.00427381, 0.00455352, 0.00477839,
575 0.00494994, 0.00506989, 0.00514012, 0.0051629, 0.00514197,
576 0.0050815, 0.0049854, 0.00485738, 0.0047009, 0.00451919,
577 0.0043152, 0.00409157, 0.00385071, 0.00359475, 0.00332561,
578 0.00304503, 0.00275452, 0.00245546, 0.00214904, 0.00183633,
579 0.00151827, 0.00119574, 0.00086945, 0.00054002, 0.000208013,
580 -0.000126036, -0.000461602, -0.000798383, -0.0011363, -0.00147529,
581 -0.00181444, -0.00215832, -0.0024967};
583 return {{x_upper, y_upper, x_lower, y_lower}};
587 AssertThrow(
false, dealii::ExcMessage(
"Invalid ONERA serial number"));
600 std::array<std::vector<double>, 4> bell(
const std::string &serial_number)
602 if (serial_number ==
"NLR-1T") {
603 std::vector<double> x_upper{
604 .0, .00259, .00974, .02185, .03796, .05675, .07753,
605 .09845, .12341, .15412, .18767, .22313, .26054, .29979,
606 .34064, .38269, .42528, .46849, .51162, .55383, .59596,
607 .63728, .67732, .71079, .73905, .76946, .80263, .84055,
608 .87846, .90845, .93589, .96199, 1.};
610 std::vector<double> y_upper{
611 .0, .00704, .01524, .02296, .02972, .03588, .04098,
612 .04469, .04741, .04986, .05188, .05345, .05459, .05531,
613 .05565, .0556, .05518, .05438, .05323, .05175, .04992,
614 .04774, .04524, .04291, .04017, .03644, .0314, .02533,
615 .01901, .01421, .0102, .00651, .00104};
617 std::vector<double> x_lower{
618 .0, .00259, .00974, .02185, .03796, .05675, .07753,
619 .09845, .12341, .15412, .18767, .22313, .26054, .29979,
620 .34064, .38269, .42528, .46849, .51162, .55383, .59596,
621 .63728, .67732, .71079, .73905, .76946, .80263, .84055,
622 .87846, .90845, .93589, .96199, 1.};
623 std::vector<double> y_lower{
624 .0, -.00512, -.00867, -.0118, -.01465, -.01713, -.01929,
625 -.02112, -.02299, -.02494, -.02671, -.02821, -.02944, -.0304,
626 -.03104, -.03142, -.0315, -.03132, -.0308, -.02992, -.02867,
627 -.02734, -.0258, -.02432, -.02305, -.02164, -.01996, -.01794,
628 -.01571, -.01364, -.01087, -.00711, -.00104};
630 return {{x_upper, y_upper, x_lower, y_lower}};
634 AssertThrow(
false, dealii::ExcMessage(
"Invalid BELL serial number"));
642 std::array<std::function<double(
const double)>, 3>
643 create_psi(
const std::vector<double> &x_upper [[maybe_unused]],
644 const std::vector<double> &y_upper [[maybe_unused]],
645 const std::vector<double> &x_lower [[maybe_unused]],
646 const std::vector<double> &y_lower [[maybe_unused]],
647 const double x_center [[maybe_unused]],
648 const double y_center [[maybe_unused]],
649 const double scaling [[maybe_unused]] = 1.)
651 Assert(x_upper.size() >= 2, dealii::ExcInternalError());
652 Assert(x_upper.front() == 0. && x_upper.back() == 1.,
653 dealii::ExcInternalError());
654 Assert(std::is_sorted(x_upper.begin(), x_upper.end()),
655 dealii::ExcInternalError());
657 Assert(x_lower.size() >= 2, dealii::ExcInternalError());
658 Assert(x_lower.front() == 0. && x_lower.back() == 1.,
659 dealii::ExcInternalError());
660 Assert(std::is_sorted(x_lower.begin(), x_lower.end()),
661 dealii::ExcInternalError());
663 Assert(y_upper.size() == x_upper.size(), dealii::ExcInternalError());
664 Assert(y_upper.front() == 0., dealii::ExcInternalError());
666 Assert(y_lower.size() == x_lower.size(), dealii::ExcInternalError());
667 Assert(y_lower.front() == 0., dealii::ExcInternalError());
669 Assert(y_lower.back() < y_upper.back(), dealii::ExcInternalError());
671 Assert(0. < x_center && x_center < 1., dealii::ExcInternalError());
673#ifdef DEAL_II_WITH_GSL
674 CubicSpline upper_airfoil(x_upper, y_upper);
676 [upper_airfoil, x_center, y_center, scaling](
const double x_hat) {
678 const double x = x_hat / scaling;
679 if (x > 1. - x_center)
680 return scaling * (upper_airfoil.eval(1.0) - y_center);
681 return scaling * (upper_airfoil.eval(x + x_center) - y_center);
684 CubicSpline lower_airfoil(x_lower, y_lower);
687 [lower_airfoil, x_center, y_center, scaling](
const double x_hat) {
689 const double x = x_hat / scaling;
690 if (x > 1. - x_center)
691 return scaling * (lower_airfoil.eval(1.0) - y_center);
692 return scaling * (lower_airfoil.eval(x + x_center) - y_center);
699 std::vector<double> x_combined;
700 std::vector<double> y_combined;
702 for (std::size_t i = 0; i < x_upper.size(); ++i) {
703 if (x_upper[i] >= x_center)
705 x_combined.push_back(x_upper[i]);
706 y_combined.push_back(y_upper[i]);
716 for (
double x : {x_center, x_center + 0.01, x_center + 0.02}) {
717 x_combined.push_back(x);
718 y_combined.push_back(upper_airfoil.eval(x));
721 std::reverse(x_combined.begin(), x_combined.end());
722 std::reverse(y_combined.begin(), y_combined.end());
723 x_combined.pop_back();
724 y_combined.pop_back();
726 for (std::size_t i = 0; i < x_lower.size(); ++i) {
727 if (x_lower[i] >= x_center)
729 x_combined.push_back(x_lower[i]);
730 y_combined.push_back(y_lower[i]);
733 for (
double x : {x_center, x_center + 0.01, x_center + 0.02}) {
734 x_combined.push_back(x);
735 y_combined.push_back(lower_airfoil.eval(x));
740 for (
unsigned int i = 0; i < y_combined.size(); ++i) {
741 const auto x = x_combined[i] - x_center;
742 const auto y = y_combined[i] - y_center;
744 const auto rho = std::sqrt(x * x + y * y);
745 auto phi = std::atan2(y, x);
747 phi += 2 * dealii::numbers::PI;
754 if (x_combined.back() == 0.)
755 x_combined.back() = 2. * dealii::numbers::PI;
756 Assert(std::is_sorted(x_combined.begin(), x_combined.end()),
757 dealii::ExcInternalError());
759 CubicSpline front_airfoil(x_combined, y_combined);
760 auto psi_front = [front_airfoil, x_center, scaling](
const double phi) {
763 return scaling * (1. - x_center);
765 return scaling * front_airfoil.eval(phi);
768 return {{psi_front, psi_upper, psi_lower}};
772 dealii::ExcNotImplemented(
"Airfoil grid needs deal.II with GSL"));
827 :
Geometry<dim>(
"airfoil", subsection)
831 airfoil_type_ =
"NASA SC(2) 0714";
833 "airfoil type", airfoil_type_,
"airfoil type and serial number");
835 airfoil_length_ = 2.;
836 this->add_parameter(
"airfoil length",
838 "length of airfoil (leading to trailing edge)");
841 this->add_parameter(
"psi samples",
843 "number of samples used for generating spline psi");
845 psi_center_[0] = 0.05;
846 this->add_parameter(
"psi center",
848 "center position of airfoil for sampling psi");
854 "Scaling parameter for averages in curved nose region, can be "
855 "adjusted by hand to equliabrate the size of faces at the nose "
856 "part of the airfoil");
859 airfoil_center_[0] = -.5;
860 this->add_parameter(
"airfoil center",
862 "position of airfoil center in the mesh");
868 "grading exponent", grading_,
"graded mesh: exponent");
870 grading_epsilon_ = 0.02;
871 this->add_parameter(
"grading epsilon",
873 "graded mesh: regularization parameter");
875 grading_epsilon_trailing_ = 0.01;
877 "grading epsilon trailing",
878 grading_epsilon_trailing_,
879 "graded mesh: regularization parameter for trailing cells");
883 "height", height_,
"height of computational domain");
886 this->add_parameter(
"width", width_,
"width of computational domain");
888 n_anisotropic_refinements_airfoil_ = 1;
890 "anisotropic pre refinement airfoil",
891 n_anisotropic_refinements_airfoil_,
892 "number of anisotropic pre refinement steps for the airfoil");
894 n_anisotropic_refinements_trailing_ = 3;
895 this->add_parameter(
"anisotropic pre refinement trailing",
896 n_anisotropic_refinements_trailing_,
897 "number of anisotropic pre refinement steps for "
898 "the blunt trailing edge cell");
901 this->add_parameter(
"subdivisions z",
903 "number of subdivisions in z direction");
915 const auto [x_upper, y_upper, x_lower, y_lower] = [&]() {
916 if (airfoil_type_.rfind(
"NACA ", 0) == 0) {
917 return naca_4digit_points(airfoil_type_.substr(5), psi_samples_);
918 }
else if (airfoil_type_.rfind(
"NASA SC(2) ", 0) == 0) {
919 return nasa_sc2(airfoil_type_.substr(11));
920 }
else if (airfoil_type_.rfind(
"ONERA ", 0) == 0) {
921 return onera(airfoil_type_.substr(6));
922 }
else if (airfoil_type_.rfind(
"BELL ", 0) == 0) {
923 return bell(airfoil_type_.substr(5));
925 AssertThrow(
false, ExcMessage(
"Unknown airfoil type"));
928 const auto [psi_front, psi_upper, psi_lower] =
944 const auto outer_radius = 0.5 * height_;
947 const auto back_length = psi_front(0.);
950 const bool sharp_trailing_edge =
951 std::abs(psi_upper(back_length) - psi_lower(back_length)) < 1.0e-10;
953 sharp_trailing_edge ||
954 std::abs(psi_upper(back_length) - psi_lower(back_length)) >
956 dealii::ExcMessage(
"Blunt trailing edge has a width of less than "
957 "0.1% of the trailing airfoil length."));
960 dealii::Triangulation<2> tria_front;
963 std::vector<dealii::Point<2>> vertices{
964 {-outer_radius, 0.0},
965 {airfoil_center_[0] - psi_front(M_PI), airfoil_center_[1]},
966 {-0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius},
967 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius},
968 {0., airfoil_center_[1] + psi_lower(-airfoil_center_[0])},
969 {airfoil_center_[0] + back_length,
970 airfoil_center_[1] + psi_lower(back_length)},
971 {0., airfoil_center_[1] + psi_upper(-airfoil_center_[0])},
972 {-0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius},
973 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius},
976 std::vector<dealii::CellData<2>> cells(4);
977 assign(cells[0].vertices, {2, 3, 4, 5});
978 assign(cells[1].vertices, {0, 2, 1, 4});
979 assign(cells[2].vertices, {7, 0, 6, 1});
980 if (sharp_trailing_edge) {
981 assign(cells[3].vertices, {8, 7, 5, 6});
983 vertices.push_back({airfoil_center_[0] + back_length,
984 airfoil_center_[1] + psi_upper(back_length)});
985 assign(cells[3].vertices, {8, 7, 9, 6});
988 tria_front.create_triangulation(
989 vertices, cells, dealii::SubCellData());
993 dealii::Triangulation<2> tria_back;
995 if (sharp_trailing_edge) {
998 const std::vector<dealii::Point<2>> vertices{
999 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius},
1000 {airfoil_center_[0] + back_length,
1001 airfoil_center_[1] + psi_lower(back_length)},
1002 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius},
1003 {outer_radius, -0.5 * outer_radius},
1004 {outer_radius, 0.0},
1005 {outer_radius, 0.5 * outer_radius},
1008 std::vector<dealii::CellData<2>> cells(2);
1009 assign(cells[0].vertices, {0, 3, 1, 4});
1010 assign(cells[1].vertices, {1, 4, 2, 5});
1012 tria_back.create_triangulation(
1013 vertices, cells, dealii::SubCellData());
1019 double trailing_height =
1020 0.5 / (0.5 +
std::pow(2., n_anisotropic_refinements_airfoil_)) *
1023 const std::vector<dealii::Point<2>> vertices{
1024 {0.5 * outer_radius, -std::sqrt(3.) / 2. * outer_radius},
1025 {airfoil_center_[0] + back_length,
1026 airfoil_center_[1] + psi_lower(back_length)},
1027 {airfoil_center_[0] + back_length,
1028 airfoil_center_[1] + psi_upper(back_length)},
1029 {0.5 * outer_radius, std::sqrt(3.) / 2. * outer_radius},
1030 {outer_radius, -0.5 * outer_radius},
1031 {outer_radius, -trailing_height},
1032 {outer_radius, trailing_height},
1033 {outer_radius, 0.5 * outer_radius},
1036 std::vector<dealii::CellData<2>> cells(3);
1037 assign(cells[0].vertices, {0, 4, 1, 5});
1038 assign(cells[1].vertices, {1, 5, 2, 6});
1039 assign(cells[2].vertices, {2, 6, 3, 7});
1041 tria_back.create_triangulation(
1042 vertices, cells, dealii::SubCellData());
1045 dealii::Triangulation<2> coarse_triangulation;
1046 GridGenerator::merge_triangulations(
1047 {&tria_front, &tria_back}, coarse_triangulation, 1.e-12,
true);
1074 for (
auto cell : coarse_triangulation.active_cell_iterators()) {
1075 for (
auto f : cell->face_indices()) {
1076 const auto face = cell->face(f);
1077 if (!face->at_boundary())
1080 bool airfoil =
true;
1081 bool spherical_boundary =
true;
1082 for (
const auto v : dealii::GeometryInfo<1>::vertex_indices())
1083 if (std::abs((face->vertex(v)).norm() - outer_radius) < 1.0e-10)
1086 spherical_boundary =
false;
1088 if (spherical_boundary) {
1089 face->set_manifold_id(3);
1090 }
else if (airfoil) {
1091 if (face->center()[0] <
1092 airfoil_center_[0] + back_length - 1.e-6) {
1093 if (face->center()[1] >= airfoil_center_[1]) {
1094 face->set_manifold_id(1);
1096 face->set_manifold_id(2);
1104 airfoil_center_, psi_front, psi_upper, psi_lower,
true, psi_ratio_};
1105 coarse_triangulation.set_manifold(1, airfoil_manifold_upper);
1113 coarse_triangulation.set_manifold(2, airfoil_manifold_lower);
1115 dealii::SphericalManifold<2> spherical_manifold;
1116 coarse_triangulation.set_manifold(3, spherical_manifold);
1123 Assert(!sharp_trailing_edge || (coarse_triangulation.n_cells() == 6),
1124 dealii::ExcInternalError());
1125 Assert(sharp_trailing_edge || (coarse_triangulation.n_cells() == 7),
1126 dealii::ExcInternalError());
1128 std::vector<std::unique_ptr<dealii::Manifold<2, 2>>> manifolds;
1129 manifolds.resize(sharp_trailing_edge ? 6 : 7);
1133 const auto first_cell = coarse_triangulation.begin_active();
1134 std::next(first_cell, 4)->set_material_id(42);
1135 std::next(first_cell, sharp_trailing_edge ? 5 : 6)->set_material_id(42);
1138 for (
auto i : {0, 1, 2, 3, 5}) {
1139 const auto index = 10 + i;
1141 dealii::Point<2> center;
1142 dealii::Tensor<1, 2> direction;
1147 Assert(i == 5, dealii::ExcInternalError());
1157 i == 5 ? grading_epsilon_trailing_ : grading_epsilon_};
1160 std::make_unique<TransfiniteInterpolationManifold<2>>();
1161 transfinite->initialize(coarse_triangulation, grading);
1163 coarse_triangulation.set_manifold(index, *transfinite);
1164 manifolds[i] = std::move(transfinite);
1168 if (sharp_trailing_edge)
1169 coarse_triangulation.reset_manifold(5);
1178 coarse_triangulation.reset_manifold(1);
1179 coarse_triangulation.reset_manifold(2);
1180 coarse_triangulation.reset_manifold(3);
1183 for (
unsigned int i = 0; i < (sharp_trailing_edge ? 6 : 7); ++i) {
1184 const auto &cell = std::next(coarse_triangulation.begin_active(), i);
1185 const auto index = 10 + i;
1186 if (i == 4 || i == (sharp_trailing_edge ? 5 : 6)) {
1187 cell->set_manifold_id(index);
1189 cell->set_all_manifold_ids(index);
1200 for (
auto cell : coarse_triangulation.active_cell_iterators())
1201 cell->set_material_id(42);
1203 std::next(first_cell, 4)->set_material_id(0);
1204 std::next(first_cell, sharp_trailing_edge ? 5 : 6)->set_material_id(0);
1207 for (
auto i : {4, sharp_trailing_edge ? 5 : 6}) {
1208 const auto index = 10 + i;
1210 std::make_unique<ryujin::TransfiniteInterpolationManifold<2>>();
1211 transfinite->initialize(coarse_triangulation);
1212 coarse_triangulation.set_manifold(index, *transfinite);
1213 manifolds[i] = std::move(transfinite);
1221 for (
unsigned int i = 0; i < (sharp_trailing_edge ? 6 : 7); ++i) {
1222 const auto &cell = std::next(coarse_triangulation.begin_active(), i);
1223 const auto index = 10 + i;
1224 cell->set_material_id(index);
1236 for (
unsigned int i = 0; i < n_anisotropic_refinements_airfoil_; ++i) {
1237 for (
auto cell : coarse_triangulation.active_cell_iterators()) {
1238 const auto id = cell->material_id();
1239 if (
id == 10 ||
id == 13)
1240 cell->set_refine_flag(dealii::RefinementCase<2>::cut_axis(0));
1243 coarse_triangulation.execute_coarsening_and_refinement();
1247 if (!sharp_trailing_edge)
1248 for (
unsigned i = 0; i < n_anisotropic_refinements_trailing_; ++i) {
1249 for (
auto cell : coarse_triangulation.active_cell_iterators())
1250 if (cell->material_id() == 15)
1251 cell->set_refine_flag(dealii::RefinementCase<2>::cut_axis(0));
1253 cell->set_refine_flag();
1254 coarse_triangulation.execute_coarsening_and_refinement();
1264 if constexpr (dim == 1) {
1265 AssertThrow(
false, dealii::ExcNotImplemented());
1268 }
else if constexpr (dim == 2) {
1270 dealii::Triangulation<2> tria3;
1271 tria3.set_mesh_smoothing(triangulation.get_mesh_smoothing());
1272 GridGenerator::flatten_triangulation(coarse_triangulation, tria3);
1274 triangulation.copy_triangulation(tria3);
1277 static_assert(dim == 3);
1280 dealii::Triangulation<2> tria3;
1281 GridGenerator::flatten_triangulation(coarse_triangulation, tria3);
1284 dealii::Triangulation<3, 3> tria4;
1285 tria4.set_mesh_smoothing(triangulation.get_mesh_smoothing());
1286 GridGenerator::extrude_triangulation(
1287 tria3, subdivisions_z_, width_, tria4);
1289 triangulation.copy_triangulation(tria4);
1298 for (
auto &cell : triangulation.active_cell_iterators()) {
1299 const auto id = cell->material_id();
1300 cell->set_all_manifold_ids(
id);
1306 if constexpr (dim == 1) {
1307 AssertThrow(
false, dealii::ExcNotImplemented());
1310 }
else if constexpr (dim == 2) {
1311 unsigned int index = 10;
1312 for (
const auto &manifold : manifolds)
1313 triangulation.set_manifold(index++, *manifold);
1316 static_assert(dim == 3);
1318 unsigned int index = 10;
1319 for (
const auto &manifold : manifolds)
1320 triangulation.set_manifold(
1326 for (
auto cell : triangulation.active_cell_iterators()) {
1327 for (
auto f : cell->face_indices()) {
1328 auto face = cell->face(f);
1331 if (!face->at_boundary())
1334 bool airfoil =
true;
1335 bool spherical_boundary =
true;
1337 const auto &indices =
1338 dealii::GeometryInfo<dim - 1>::vertex_indices();
1339 for (
const auto v : indices) {
1340 const auto vert = face->vertex(v);
1341 const auto radius_sqr = vert[0] * vert[0] + vert[1] * vert[1];
1342 if (radius_sqr >= outer_radius * outer_radius - 1.0e-10 ||
1343 vert[0] > airfoil_center_[0] + 1.001 * back_length)
1346 spherical_boundary =
false;
1349 bool periodic_face = (dim == 3);
1351 if constexpr (dim == 3) {
1352 const auto &indices =
1353 dealii::GeometryInfo<dim - 1>::vertex_indices();
1354 bool not_left =
false;
1355 bool not_right =
false;
1356 for (
const auto v : indices) {
1357 const auto vert = face->vertex(v);
1358 if (vert[2] > 1.0e-10)
1360 if (vert[2] < width_ - 1.0e-10)
1362 if (not_left && not_right) {
1363 periodic_face =
false;
1369 if (periodic_face) {
1371 }
else if (spherical_boundary) {
1373 }
else if (airfoil) {
1376 Assert(
false, dealii::ExcInternalError());
1384#ifndef BUG_COLLECT_PERIODIC_FACES_INSTANTIATION
1385 if constexpr (dim == 3) {
1386 std::vector<dealii::GridTools::PeriodicFacePair<
1387 typename dealii::Triangulation<dim>::cell_iterator>>
1390 GridTools::collect_periodic_faces(triangulation,
1395 triangulation.add_periodicity(periodic_faces);
1401 dealii::Point<2> airfoil_center_;
1402 double airfoil_length_;
1403 std::string airfoil_type_;
1404 dealii::Point<2> psi_center_;
1406 unsigned int psi_samples_;
1410 double grading_epsilon_;
1411 double grading_epsilon_trailing_;
1412 unsigned int n_anisotropic_refinements_airfoil_;
1413 unsigned int n_anisotropic_refinements_trailing_;
1414 unsigned int subdivisions_z_;
Airfoil(const std::string subsection)
void create_triangulation(dealii::Triangulation< dim > &triangulation) final
dealii::Point< dim > pull_back(const dealii::Point< dim > &space_point) const final
std::unique_ptr< dealii::Manifold< dim, dim > > clone() const final
dealii::Point< dim > push_forward(const dealii::Point< dim > &point) const final
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)
Point< dim > get_new_point(const ArrayView< const Point< dim > > &surrounding_points, const ArrayView< const double > &weights) const override
std::unique_ptr< Manifold< dim > > clone() const override
ExtrudedManifold(const dealii::Manifold< dim - 1 > &manifold)
dealii::Point< dim > push_forward(const dealii::Point< dim > &chart_point) const final
GradingManifold(const dealii::Point< dim > center, const dealii::Tensor< 1, dim > direction, const double grading, const double epsilon)
std::unique_ptr< dealii::Manifold< dim, dim > > clone() const final
dealii::Point< dim > pull_back(const dealii::Point< dim > &space_point) const final
Point< dim > get_new_point(const ArrayView< const Point< dim > > &surrounding_points, const ArrayView< const double > &weights) const override
T pow(const T x, const T b)