51 : airfoil_type(
"NACA")
53 , joukowski_center(-0.1, 0.14)
59 , incline_factor(0.35)
62 , n_subdivision_x_0(3)
63 , n_subdivision_x_1(2)
64 , n_subdivision_x_2(5)
66 , airfoil_sampling_factor(2)
69 airfoil_length <= height,
71 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
73 Assert(incline_factor < 1.0 && incline_factor >= 0.0,
74 ExcMessage(
"incline_factor has to be in [0,1)!"));
87 "Mesh height measured from airfoil nose to horizontal boundaries");
91 "Length measured from airfoil leading edge to vertical outlet boundary");
95 "Define obliqueness of the vertical mesh around the airfoil");
104 "Type of airfoil geometry, either NACA or Joukowski airfoil",
119 "Joukowski circle center coordinates");
122 "Joukowski airfoil length leading to trailing edge");
130 "Number of global refinements");
132 "NumberSubdivisionX0",
134 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
136 "NumberSubdivisionX1",
138 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
140 "NumberSubdivisionX2",
142 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
145 "Number of subdivisions normal to airfoil");
149 "Factor to obtain a finer mesh at the airfoil surface");
164 static const unsigned int id_block_1 = 1;
165 static const unsigned int id_block_2 = 2;
166 static const unsigned int id_block_3 = 3;
167 static const unsigned int id_block_4 = 4;
168 static const unsigned int id_block_5 = 5;
169 static const unsigned int id_block_6 = 6;
174 MeshGenerator(
const AdditionalData &
data)
175 : refinements(
data.refinements)
176 , n_subdivision_x_0(
data.n_subdivision_x_0)
177 , n_subdivision_x_1(
data.n_subdivision_x_1)
178 , n_subdivision_x_2(
data.n_subdivision_x_2)
179 , n_subdivision_y(
data.n_subdivision_y)
180 , height(
data.height)
181 , length_b2(
data.length_b2)
182 , incline_factor(
data.incline_factor)
183 , bias_factor(
data.bias_factor)
185 , n_cells_x_0(
Utilities::
pow(2, refinements) * n_subdivision_x_0)
186 , n_cells_x_1(
Utilities::
pow(2, refinements) * n_subdivision_x_1)
187 , n_cells_x_2(
Utilities::
pow(2, refinements) * n_subdivision_x_2)
188 , n_cells_y(
Utilities::
pow(2, refinements) * n_subdivision_y)
189 , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
191 , airfoil_1D(set_airfoil_length(
193 data.airfoil_type ==
"Joukowski" ?
194 joukowski(
data.joukowski_center,
195 n_points_on_each_side,
196 data.airfoil_sampling_factor) :
197 (
data.airfoil_type ==
"NACA" ?
199 n_points_on_each_side,
200 data.airfoil_sampling_factor) :
203 std::vector<Point<2>>{
207 data.airfoil_length))
208 , end_b0_x_u(airfoil_1D[0][n_cells_x_0][0])
209 , end_b0_x_l(airfoil_1D[1][n_cells_x_0][0])
210 , nose_x(airfoil_1D[0].front()[0])
211 , tail_x(airfoil_1D[0].back()[0])
212 , tail_y(airfoil_1D[0].back()[1])
213 , center_mesh(0.5 *
std::abs(end_b0_x_u + end_b0_x_l))
214 , length_b1_x(tail_x - center_mesh)
216 (edge_length +
std::abs(nose_x - center_mesh))))
220 ,
A(nose_x - edge_length, 0)
223 , D(center_mesh, height)
225 ,
F(center_mesh, -height)
229 , J(tail_x + length_b2, 0)
234 data.airfoil_type ==
"NACA",
247 make_coarse_grid(tria_grid);
249 set_boundary_ids(tria_grid);
251 if (periodic_faces !=
nullptr)
254 tria_grid, 5, 4, 1, *periodic_faces);
272 (void)periodic_faces;
279 const unsigned int refinements;
282 const unsigned int n_subdivision_x_0;
285 const unsigned int n_subdivision_x_1;
288 const unsigned int n_subdivision_x_2;
292 const unsigned int n_subdivision_y;
299 const double length_b2;
303 const double incline_factor;
306 const double bias_factor;
309 const double edge_length;
312 const unsigned int n_cells_x_0;
315 const unsigned int n_cells_x_1;
318 const unsigned int n_cells_x_2;
322 const unsigned int n_cells_y;
325 const unsigned int n_points_on_each_side;
329 const std::array<std::vector<Point<2>>, 2> airfoil_1D;
333 const double end_b0_x_u;
337 const double end_b0_x_l;
351 const double center_mesh;
354 const double length_b1_x;
381 const Point<2> A, B,
C, D,
E,
F, G, H, I, J,
K,
L;
420 static std::array<std::vector<Point<2>>, 2>
421 joukowski(
const Point<2> ¢erpoint,
422 const unsigned int number_points,
423 const unsigned int factor)
425 std::array<std::vector<Point<2>>, 2> airfoil_1D;
426 const unsigned int total_points = 2 * number_points - 2;
427 const unsigned int n_airfoilpoints = factor * total_points;
429 const auto jouk_points =
430 joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
433 std::vector<Point<2>> upper_points;
434 std::vector<Point<2>> lower_points;
438 unsigned int nose_index = 0;
439 unsigned int tail_index = 0;
440 double nose_x_coordinate = 0;
441 double tail_x_coordinate = 0;
445 for (
unsigned int i = 0; i < jouk_points.size(); ++i)
447 if (jouk_points[i][0] < nose_x_coordinate)
449 nose_x_coordinate = jouk_points[i][0];
452 if (jouk_points[i][0] > tail_x_coordinate)
454 tail_x_coordinate = jouk_points[i][0];
460 for (
unsigned int i = tail_index; i < jouk_points.size(); ++i)
461 upper_points.emplace_back(jouk_points[i]);
462 for (
unsigned int i = 0; i <= nose_index; ++i)
463 upper_points.emplace_back(jouk_points[i]);
464 std::reverse(upper_points.begin(), upper_points.end());
467 lower_points.insert(lower_points.end(),
468 jouk_points.begin() + nose_index,
469 jouk_points.begin() + tail_index + 1);
472 airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
473 airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
476 auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
477 const double nose_x_pos = vector.front()[0];
478 for (
auto &i : vector)
482 move_nose_to_origin(airfoil_1D[1]);
483 move_nose_to_origin(airfoil_1D[0]);
512 static std::vector<Point<2>>
513 joukowski_circle(
const Point<2> ¢er,
514 const unsigned int number_points)
516 std::vector<Point<2>> circle_points;
523 const double radius =
std::sqrt(center[1] * center[1] +
524 (1 - center[0]) * (1 - center[0]));
526 center[1] * center[1] + (1 + center[0]) * (1 + center[0]));
530 radius_test < radius,
532 "Error creating lower circle: Circle for Joukowski-transform does"
533 " not enclose point zeta = -1! Choose different center "
537 const double theta = 2 *
numbers::PI / number_points;
539 circle_points.reserve(number_points);
540 for (
unsigned int i = 0; i < number_points; ++i)
541 circle_points.emplace_back(center[0] - radius *
std::cos(i * theta),
545 return circle_points;
556 static std::vector<Point<2>>
557 joukowski_transform(
const std::vector<
Point<2>> &circle_points)
559 std::vector<Point<2>> joukowski_points(circle_points.size());
562 for (
unsigned int i = 0; i < circle_points.size(); ++i)
564 const double chi = circle_points[i][0];
565 const double eta = circle_points[i][1];
566 const std::complex<double> zeta(chi, eta);
567 const std::complex<double> z = zeta + 1. / zeta;
569 joukowski_points[i] = {std::real(z), std::imag(z)};
571 return joukowski_points;
590 static std::array<std::vector<Point<2>>, 2>
591 naca(
const std::string &serialnumber,
592 const unsigned int number_points,
593 const unsigned int factor)
597 const unsigned int n_airfoilpoints = factor * number_points;
600 return {{make_points_equidistant(
601 naca_create_points(serialnumber, n_airfoilpoints,
true),
603 make_points_equidistant(
604 naca_create_points(serialnumber, n_airfoilpoints,
false),
619 static std::vector<Point<2>>
620 naca_create_points(
const std::string &serialnumber,
621 const unsigned int number_points,
624 Assert(serialnumber.size() == 4,
625 ExcMessage(
"This NACA-serial number is not implemented!"));
627 return naca_create_points_4_digits(serialnumber,
646 static std::vector<Point<2>>
647 naca_create_points_4_digits(
const std::string &serialnumber,
648 const unsigned int number_points,
652 const unsigned int digit_0 = (serialnumber[0] -
'0');
653 const unsigned int digit_1 = (serialnumber[1] -
'0');
654 const unsigned int digit_2 = (serialnumber[2] -
'0');
655 const unsigned int digit_3 = (serialnumber[3] -
'0');
657 const unsigned int digit_23 = 10 * digit_2 + digit_3;
660 const double t =
static_cast<double>(digit_23) / 100.0;
662 std::vector<Point<2>> naca_points;
664 if (digit_0 == 0 && digit_1 == 0)
665 for (
unsigned int i = 0; i < number_points; ++i)
667 const double x = i * 1 / (1.0 * number_points - 1);
671 0.3516 * Utilities::fixed_power<2>(x) +
672 0.2843 * Utilities::fixed_power<3>(x) -
673 0.1036 * Utilities::fixed_power<4>(
677 naca_points.emplace_back(x, +y_t);
679 naca_points.emplace_back(x, -y_t);
682 for (
unsigned int i = 0; i < number_points; ++i)
684 const double m = 1.0 * digit_0 / 100;
685 const double p = 1.0 * digit_1 / 10;
686 const double x = i * 1 / (1.0 * number_points - 1);
690 m / Utilities::fixed_power<2>(p) *
691 (2 * p * x - Utilities::fixed_power<2>(x)) :
696 (x <= p) ? 2 * m / Utilities::fixed_power<2>(p) * (p - x) :
702 0.3516 * Utilities::fixed_power<2>(x) +
703 0.2843 * Utilities::fixed_power<3>(x) -
704 0.1036 * Utilities::fixed_power<4>(
710 naca_points.emplace_back(x - y_t *
std::sin(theta),
713 naca_points.emplace_back(x + y_t *
std::sin(theta),
730 static std::array<std::vector<Point<2>>, 2>
731 set_airfoil_length(
const std::array<std::vector<
Point<2>>, 2> &input,
732 const double desired_len)
734 std::array<std::vector<Point<2>>, 2> output;
735 output[0] = set_airfoil_length(input[0], desired_len);
736 output[1] = set_airfoil_length(input[1], desired_len);
748 static std::vector<Point<2>>
749 set_airfoil_length(
const std::vector<
Point<2>> &input,
750 const double desired_len)
752 std::vector<Point<2>> output = input;
755 desired_len / input.front().distance(input.back());
757 for (
auto &x : output)
773 static std::vector<Point<2>>
774 make_points_equidistant(
775 const std::vector<
Point<2>> &non_equidistant_points,
776 const unsigned int number_points)
778 const unsigned int n_points =
779 non_equidistant_points
783 std::vector<double> arclength_L(non_equidistant_points.size(), 0);
784 for (
unsigned int i = 0; i < non_equidistant_points.size() - 1; ++i)
787 non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
790 const auto airfoil_length =
792 const auto deltaX = airfoil_length / (number_points - 1);
796 std::vector<Point<2>> equidist(
799 equidist[0] = non_equidistant_points[0];
800 equidist[number_points - 1] = non_equidistant_points[n_points - 1];
804 for (
unsigned int j = 0, i = 1; j < n_points - 1; ++j)
807 const auto Lj = arclength_L[j];
808 const auto Ljp = arclength_L[j + 1];
810 while (Lj <= i * deltaX && i * deltaX <= Ljp &&
811 i < number_points - 1)
813 equidist[i] =
Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
814 (non_equidistant_points[j + 1] -
815 non_equidistant_points[j]) +
816 non_equidistant_points[j]);
836 std::vector<Triangulation<2>> trias(10);
840 const std::vector<Point<2>> &corner_vertices,
841 const std::vector<unsigned int> &repetitions,
853 auto &
point = it->vertex();
854 const double xi =
point[0];
855 const double eta =
point[1];
858 point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
859 (1 + xi) * (1 - eta) * corner_vertices[1] +
860 (1 - xi) * (1 + eta) * corner_vertices[2] +
861 (1 + xi) * (1 + eta) * corner_vertices[3]);
865 for (
auto cell : tria.active_cell_iterators())
873 {n_subdivision_y, n_subdivision_x_0},
877 {n_subdivision_y, n_subdivision_x_0},
881 {n_subdivision_x_1, n_subdivision_y},
885 {n_subdivision_x_1, n_subdivision_y},
889 {n_subdivision_x_2, n_subdivision_y},
893 {n_subdivision_x_2, n_subdivision_y},
920 for (
auto cell : tria.active_cell_iterators())
923 if (cell->face(f)->at_boundary() ==
false)
926 const auto mid = cell->material_id();
928 if ((mid == id_block_1 && f == 0) ||
929 (mid == id_block_4 && f == 0))
930 cell->face(f)->set_boundary_id(0);
931 else if ((mid == id_block_3 && f == 0) ||
932 (mid == id_block_6 && f == 2))
933 cell->face(f)->set_boundary_id(1);
934 else if ((mid == id_block_1 && f == 1) ||
935 (mid == id_block_2 && f == 1))
936 cell->face(f)->set_boundary_id(2);
937 else if ((mid == id_block_4 && f == 1) ||
938 (mid == id_block_5 && f == 3))
939 cell->face(f)->set_boundary_id(3);
940 else if ((mid == id_block_2 && f == 0) ||
941 (mid == id_block_3 && f == 2))
942 cell->face(f)->set_boundary_id(4);
943 else if ((mid == id_block_5 && f == 2) ||
944 (mid == id_block_6 && f == 0))
945 cell->face(f)->set_boundary_id(5);
965 std::vector<bool> vertex_processed(tria.
n_vertices(),
false);
982 for (
const auto &cell : tria.cell_iterators())
986 if (vertex_processed[cell->vertex_index(v)])
990 vertex_processed[cell->vertex_index(v)] =
true;
992 auto &node = cell->vertex(v);
995 if (cell->material_id() == id_block_1 ||
996 cell->material_id() == id_block_4)
1005 if (cell->material_id() == id_block_1)
1008 (node - horizontal_offset) +
1014 else if (cell->material_id() == id_block_4)
1017 (node - horizontal_offset) -
1023 const double trapeze_height =
1025 const double L = height /
std::sin(gamma);
1026 const double l_a =
std::cos(gamma) * edge_length;
1027 const double l_b = trapeze_height *
std::tan(gamma);
1029 const double x2 =
L - l_a - l_b;
1031 const double Dx = x1 + x2 + x3;
1032 const double deltax =
1034 const double dx = Dx / n_cells_x_0;
1035 const double dy = trapeze_height / n_cells_y;
1037 static_cast<int>(std::round((node_[0] - deltax) / dx));
1039 static_cast<int>(std::round(
std::abs(node_[1]) / dy));
1041 node_[0] =
numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1042 node_[1] = height * (1.0 * iy) / n_cells_y;
1049 const double dy = height / n_cells_y;
1051 static_cast<int>(std::round(node_[0] / dx));
1053 static_cast<int>(std::round(node_[1] / dy));
1054 const double alpha =
1055 bias_alpha(1 - (1.0 * iy) / n_cells_y);
1056 const double theta = node_[0];
1058 ((cell->material_id() == id_block_1) ?
1063 (cell->material_id() == id_block_1) ? (0) : (1))]
1069 else if (cell->material_id() == id_block_2 ||
1070 cell->material_id() == id_block_5)
1078 "Points D,C,G and E,F,I are not defined symmetric to "
1079 "x-axis, which is required to interpolate block 2"
1080 " and 5 with same geometric computations."));
1081 const double l_y = D[1] -
C[1];
1082 const double l_h = D[1] - l_y;
1083 const double by = -l_h / length_b1_x * (node[0] - H[0]);
1084 const double dy = (height - by) / n_cells_y;
1085 const int iy =
static_cast<int>(
1086 std::round((
std::abs(node[1]) - by) / dy));
1087 const double dx = length_b1_x / n_cells_x_1;
1088 const int ix =
static_cast<int>(
1089 std::round(
std::abs(node[0] - center_mesh) / dx));
1091 const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1095 const Point<2> p(ix * dx + center_mesh +
1096 incline_factor * length_b2 * ix /
1098 ((cell->material_id() == id_block_2) ?
1104 (cell->material_id() == id_block_2) ? (0) : (1))]
1105 [n_cells_x_0 + ix] *
1109 else if (cell->material_id() == id_block_3 ||
1110 cell->material_id() == id_block_6)
1113 const double dx = length_b2 / n_cells_x_2;
1114 const double dy = height / n_cells_y;
1115 const int ix =
static_cast<int>(
1116 std::round(
std::abs(node[0] - H[0]) / dx));
1118 static_cast<int>(std::round(
std::abs(node[1]) / dy));
1120 const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1121 const double alpha_x =
1122 bias_alpha(1 - (
static_cast<double>(ix)) / n_cells_x_2);
1126 const Point<2> p1(J[0] - (1 - incline_factor) * length_b2 *
1128 ((cell->material_id() == id_block_3) ?
1133 const Point<2> p2(J[0] - alpha_x * length_b2, tail_y);
1134 node = p1 * (1 - alpha_y) + p2 * alpha_y;
1154 bias_alpha(
double alpha)
const
1164 internal_create_triangulation(
1168 const AdditionalData &additional_data)
1170 MeshGenerator mesh_generator(additional_data);
1173 if (
auto *parallel_tria =
1176 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1177 else if (
auto *parallel_tria =
dynamic_cast<
1180 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1182 mesh_generator.create_triangulation(tria, periodic_faces);
1199 const AdditionalData &)
1209 const AdditionalData &additional_data)
1211 internal_create_triangulation(tria,
nullptr, additional_data);
1222 const AdditionalData &additional_data)
1224 internal_create_triangulation(tria, &periodic_faces, additional_data);
1235 const AdditionalData &additional_data)
1239 (void)additional_data;
1240 (void)periodic_faces;
1251 template <
int dim,
int spacedim>
1263 cell->face(f)->set_boundary_id(f);
1269 template <
int spacedim>
1280 if (cell->center()[0] > 0)
1281 cell->set_material_id(1);
1288 template <
int dim,
int spacedim>
1293 const double epsilon)
1306 for (; face != endface; ++face)
1307 if (face->at_boundary())
1308 if (face->boundary_id() == 0)
1312 if (
std::abs(center[0] - p1[0]) < epsilon)
1313 face->set_boundary_id(0);
1314 else if (
std::abs(center[0] - p2[0]) < epsilon)
1315 face->set_boundary_id(1);
1316 else if (dim > 1 &&
std::abs(center[1] - p1[1]) < epsilon)
1317 face->set_boundary_id(2);
1318 else if (dim > 1 &&
std::abs(center[1] - p2[1]) < epsilon)
1319 face->set_boundary_id(3);
1320 else if (dim > 2 &&
std::abs(center[2] - p1[2]) < epsilon)
1321 face->set_boundary_id(4);
1322 else if (dim > 2 &&
std::abs(center[2] - p2[2]) < epsilon)
1323 face->set_boundary_id(5);
1332 for (
const auto &cell : tria.cell_iterators())
1335 for (
unsigned int d = 0;
d < dim; ++
d)
1336 if (cell->center()[d] > 0)
1338 cell->set_material_id(
id);
1347 template <
int spacdim>
1360 for (
auto cell = tria.
begin(); cell != tria.
end(); ++cell)
1363 cell->face(2)->set_all_boundary_ids(1);
1387 cell->face(4)->set_all_boundary_ids(1);
1391 cell->face(2)->set_all_boundary_ids(1);
1395 cell->face(2)->set_all_boundary_ids(1);
1399 cell->face(0)->set_all_boundary_ids(1);
1403 cell->face(2)->set_all_boundary_ids(1);
1407 cell->face(0)->set_all_boundary_ids(1);
1409 else if (tria.
n_cells() == 12)
1418 cell->face(5)->set_all_boundary_ids(1);
1421 else if (tria.
n_cells() == 96)
1428 unsigned int count = 0;
1429 for (
const auto &cell : tria.cell_iterators())
1430 if (cell->face(5)->at_boundary())
1432 cell->face(5)->set_all_boundary_ids(1);
1452 const double inner_radius,
1453 const double outer_radius)
1458 double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
1459 double eps = 1
e-3 * middle;
1462 for (; cell != tria.
end(); ++cell)
1463 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
1465 if (!cell->face(f)->at_boundary())
1468 double radius = cell->face(f)->center().norm() - center.
norm();
1469 if (std::fabs(cell->face(f)->center()[0]) <
1472 cell->face(f)->set_boundary_id(2);
1473 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1475 if (cell->face(f)->line(j)->at_boundary())
1476 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1477 cell->face(f)->line(j)->vertex(1).norm()) >
1479 cell->face(f)->line(j)->set_boundary_id(2);
1481 else if (std::fabs(cell->face(f)->center()[1]) <
1484 cell->face(f)->set_boundary_id(3);
1485 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1487 if (cell->face(f)->line(j)->at_boundary())
1488 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1489 cell->face(f)->line(j)->vertex(1).norm()) >
1491 cell->face(f)->line(j)->set_boundary_id(3);
1493 else if (std::fabs(cell->face(f)->center()[2]) <
1496 cell->face(f)->set_boundary_id(4);
1497 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1499 if (cell->face(f)->line(j)->at_boundary())
1500 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1501 cell->face(f)->line(j)->vertex(1).norm()) >
1503 cell->face(f)->line(j)->set_boundary_id(4);
1505 else if (radius < middle)
1507 cell->face(f)->set_boundary_id(0);
1508 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1510 if (cell->face(f)->line(j)->at_boundary())
1511 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1512 cell->face(f)->line(j)->vertex(1).norm()) <
1514 cell->face(f)->line(j)->set_boundary_id(0);
1516 else if (radius > middle)
1518 cell->face(f)->set_boundary_id(1);
1519 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1521 if (cell->face(f)->line(j)->at_boundary())
1522 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1523 cell->face(f)->line(j)->vertex(1).norm()) <
1525 cell->face(f)->line(j)->set_boundary_id(1);
1535 template <
int dim,
int spacedim>
1546 for (
unsigned int i = 0; i < dim; ++i)
1560 vertices[0] = vertices[1] = p1;
1561 vertices[2] = vertices[3] = p2;
1563 vertices[1][0] = p2[0];
1564 vertices[2][0] = p1[0];
1567 vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
1568 vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
1570 vertices[1][0] = p2[0];
1571 vertices[2][1] = p2[1];
1572 vertices[3][0] = p2[0];
1573 vertices[3][1] = p2[1];
1575 vertices[4][0] = p1[0];
1576 vertices[4][1] = p1[1];
1577 vertices[5][1] = p1[1];
1578 vertices[6][0] = p1[0];
1586 std::vector<CellData<dim>> cells(1);
1588 cells[0].vertices[i] = i;
1589 cells[0].material_id = 0;
1595 colorize_hyper_rectangle(tria);
1600 template <
int dim,
int spacedim>
1608 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1611 for (
unsigned int i = 0; i < dim; ++i)
1632 for (
unsigned int d = 0;
d < dim; ++
d)
1633 for (
unsigned int c = 1; c <= dim; ++c)
1634 vector_matrix[c - 1][d] = vertices[c][d] - vertices[0][d];
1637 "Vertices of simplex must form a right handed system"));
1641 std::vector<Point<dim>> points = vertices;
1645 for (
unsigned int i = 0; i <= dim; ++i)
1647 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1648 center += points[i];
1653 for (
unsigned int i = 1; i < dim; ++i)
1654 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1656 for (
unsigned int i = 0; i <= dim; ++i)
1657 points.push_back(1. / 3. *
1658 (points[i] + points[(i + 1) % (dim + 1)] +
1659 points[(i + 2) % (dim + 1)]));
1661 points.push_back((1. / (dim + 1)) * center);
1663 std::vector<CellData<dim>> cells(dim + 1);
1668 cells[0].vertices[0] = 0;
1669 cells[0].vertices[1] = 3;
1670 cells[0].vertices[2] = 5;
1671 cells[0].vertices[3] = 6;
1672 cells[0].material_id = 0;
1674 cells[1].vertices[0] = 3;
1675 cells[1].vertices[1] = 1;
1676 cells[1].vertices[2] = 6;
1677 cells[1].vertices[3] = 4;
1678 cells[1].material_id = 0;
1680 cells[2].vertices[0] = 5;
1681 cells[2].vertices[1] = 6;
1682 cells[2].vertices[2] = 2;
1683 cells[2].vertices[3] = 4;
1684 cells[2].material_id = 0;
1688 cells[0].vertices[0] = 0;
1689 cells[0].vertices[1] = 4;
1690 cells[0].vertices[2] = 8;
1691 cells[0].vertices[3] = 10;
1692 cells[0].vertices[4] = 7;
1693 cells[0].vertices[5] = 13;
1694 cells[0].vertices[6] = 12;
1695 cells[0].vertices[7] = 14;
1696 cells[0].material_id = 0;
1698 cells[1].vertices[0] = 4;
1699 cells[1].vertices[1] = 1;
1700 cells[1].vertices[2] = 10;
1701 cells[1].vertices[3] = 5;
1702 cells[1].vertices[4] = 13;
1703 cells[1].vertices[5] = 9;
1704 cells[1].vertices[6] = 14;
1705 cells[1].vertices[7] = 11;
1706 cells[1].material_id = 0;
1708 cells[2].vertices[0] = 8;
1709 cells[2].vertices[1] = 10;
1710 cells[2].vertices[2] = 2;
1711 cells[2].vertices[3] = 5;
1712 cells[2].vertices[4] = 12;
1713 cells[2].vertices[5] = 14;
1714 cells[2].vertices[6] = 6;
1715 cells[2].vertices[7] = 11;
1716 cells[2].material_id = 0;
1718 cells[3].vertices[0] = 7;
1719 cells[3].vertices[1] = 13;
1720 cells[3].vertices[2] = 12;
1721 cells[3].vertices[3] = 14;
1722 cells[3].vertices[4] = 3;
1723 cells[3].vertices[5] = 9;
1724 cells[3].vertices[6] = 6;
1725 cells[3].vertices[7] = 11;
1726 cells[3].material_id = 0;
1736 template <
int dim,
int spacedim>
1743 if (reference_cell == ReferenceCells::get_hypercube<dim>())
1752 std::vector<Point<spacedim>> vertices(
reference_cell.n_vertices());
1756 for (
unsigned int d = 0;
d < dim; ++
d)
1757 vertices[v][d] = this_vertex[d];
1765 std::vector<CellData<dim>> cells(1);
1768 cells[0].vertices[v] = v;
1777 const unsigned int n_cells,
1778 const unsigned int n_rotations,
1782 const unsigned int dim = 3;
1785 "More than 4 cells are needed to create a moebius grid."));
1787 ExcMessage(
"Outer and inner radius must be positive."));
1789 ExcMessage(
"Outer radius must be greater than inner radius."));
1792 std::vector<Point<dim>> vertices(4 * n_cells);
1796 for (
unsigned int i = 0; i <
n_cells; ++i)
1797 for (
unsigned int j = 0; j < 4; ++j)
1799 vertices[4 * i + j][0] =
1803 vertices[4 * i + j][1] =
1807 vertices[4 * i + j][2] =
1811 unsigned int offset = 0;
1816 static constexpr std::array<unsigned int, 8> local_vertex_numbering{
1817 {0, 1, 5, 4, 2, 3, 7, 6}};
1818 std::vector<CellData<dim>> cells(n_cells);
1819 for (
unsigned int i = 0; i <
n_cells; ++i)
1821 for (
unsigned int j = 0; j < 2; ++j)
1823 cells[i].vertices[local_vertex_numbering[0 + 4 * j]] =
1825 cells[i].vertices[local_vertex_numbering[1 + 4 * j]] =
1827 cells[i].vertices[local_vertex_numbering[2 + 4 * j]] =
1829 cells[i].vertices[local_vertex_numbering[3 + 4 * j]] =
1833 cells[i].material_id = 0;
1837 cells[
n_cells - 1].vertices[local_vertex_numbering[4]] =
1838 (0 + n_rotations) % 4;
1839 cells[
n_cells - 1].vertices[local_vertex_numbering[5]] =
1840 (3 + n_rotations) % 4;
1841 cells[
n_cells - 1].vertices[local_vertex_numbering[6]] =
1842 (2 + n_rotations) % 4;
1843 cells[
n_cells - 1].vertices[local_vertex_numbering[7]] =
1844 (1 + n_rotations) % 4;
1855 const double centerline_radius,
1856 const double inner_radius,
1860 Assert(centerline_radius > inner_radius,
1861 ExcMessage(
"The centerline radius must be greater than the "
1863 Assert(inner_radius > 0.0,
1864 ExcMessage(
"The inner radius must be positive."));
1866 const unsigned int dim = 2;
1867 const unsigned int spacedim = 3;
1868 std::vector<Point<spacedim>> vertices(16);
1870 vertices[0] =
Point<spacedim>(centerline_radius - inner_radius, 0, 0);
1872 vertices[2] =
Point<spacedim>(centerline_radius + inner_radius, 0, 0);
1874 vertices[4] =
Point<spacedim>(0, 0, centerline_radius - inner_radius);
1876 vertices[6] =
Point<spacedim>(0, 0, centerline_radius + inner_radius);
1878 vertices[8] =
Point<spacedim>(-(centerline_radius - inner_radius), 0, 0);
1880 vertices[10] =
Point<spacedim>(-(centerline_radius + inner_radius), 0, 0);
1882 vertices[12] =
Point<spacedim>(0, 0, -(centerline_radius - inner_radius));
1884 vertices[14] =
Point<spacedim>(0, 0, -(centerline_radius + inner_radius));
1887 std::vector<CellData<dim>> cells(16);
1889 cells[0].vertices[0] = 0;
1890 cells[0].vertices[1] = 4;
1891 cells[0].vertices[2] = 3;
1892 cells[0].vertices[3] = 7;
1893 cells[0].material_id = 0;
1895 cells[1].vertices[0] = 1;
1896 cells[1].vertices[1] = 5;
1897 cells[1].vertices[2] = 0;
1898 cells[1].vertices[3] = 4;
1899 cells[1].material_id = 0;
1901 cells[2].vertices[0] = 2;
1902 cells[2].vertices[1] = 6;
1903 cells[2].vertices[2] = 1;
1904 cells[2].vertices[3] = 5;
1905 cells[2].material_id = 0;
1907 cells[3].vertices[0] = 3;
1908 cells[3].vertices[1] = 7;
1909 cells[3].vertices[2] = 2;
1910 cells[3].vertices[3] = 6;
1911 cells[3].material_id = 0;
1913 cells[4].vertices[0] = 4;
1914 cells[4].vertices[1] = 8;
1915 cells[4].vertices[2] = 7;
1916 cells[4].vertices[3] = 11;
1917 cells[4].material_id = 0;
1919 cells[5].vertices[0] = 5;
1920 cells[5].vertices[1] = 9;
1921 cells[5].vertices[2] = 4;
1922 cells[5].vertices[3] = 8;
1923 cells[5].material_id = 0;
1925 cells[6].vertices[0] = 6;
1926 cells[6].vertices[1] = 10;
1927 cells[6].vertices[2] = 5;
1928 cells[6].vertices[3] = 9;
1929 cells[6].material_id = 0;
1931 cells[7].vertices[0] = 7;
1932 cells[7].vertices[1] = 11;
1933 cells[7].vertices[2] = 6;
1934 cells[7].vertices[3] = 10;
1935 cells[7].material_id = 0;
1937 cells[8].vertices[0] = 8;
1938 cells[8].vertices[1] = 12;
1939 cells[8].vertices[2] = 11;
1940 cells[8].vertices[3] = 15;
1941 cells[8].material_id = 0;
1943 cells[9].vertices[0] = 9;
1944 cells[9].vertices[1] = 13;
1945 cells[9].vertices[2] = 8;
1946 cells[9].vertices[3] = 12;
1947 cells[9].material_id = 0;
1949 cells[10].vertices[0] = 10;
1950 cells[10].vertices[1] = 14;
1951 cells[10].vertices[2] = 9;
1952 cells[10].vertices[3] = 13;
1953 cells[10].material_id = 0;
1955 cells[11].vertices[0] = 11;
1956 cells[11].vertices[1] = 15;
1957 cells[11].vertices[2] = 10;
1958 cells[11].vertices[3] = 14;
1959 cells[11].material_id = 0;
1961 cells[12].vertices[0] = 12;
1962 cells[12].vertices[1] = 0;
1963 cells[12].vertices[2] = 15;
1964 cells[12].vertices[3] = 3;
1965 cells[12].material_id = 0;
1967 cells[13].vertices[0] = 13;
1968 cells[13].vertices[1] = 1;
1969 cells[13].vertices[2] = 12;
1970 cells[13].vertices[3] = 0;
1971 cells[13].material_id = 0;
1973 cells[14].vertices[0] = 14;
1974 cells[14].vertices[1] = 2;
1975 cells[14].vertices[2] = 13;
1976 cells[14].vertices[3] = 1;
1977 cells[14].material_id = 0;
1979 cells[15].vertices[0] = 15;
1980 cells[15].vertices[1] = 3;
1981 cells[15].vertices[2] = 14;
1982 cells[15].vertices[3] = 2;
1983 cells[15].material_id = 0;
1996 static constexpr int circle_cell_vertices[5][4] = {{0, 1, 2, 3},
2008 const double centerline_radius,
2009 const double inner_radius,
2010 const unsigned int n_cells_toroidal,
2013 Assert(centerline_radius > inner_radius,
2014 ExcMessage(
"The centerline radius must be greater than the "
2016 Assert(inner_radius > 0.0,
2017 ExcMessage(
"The inner radius must be positive."));
2019 ExcMessage(
"Number of cells in toroidal direction has "
2020 "to be at least 3 for a torus of polar extent 2*pi."));
2026 const double a = 1. / (1 +
std::sqrt(2.0));
2029 const unsigned int additional_layer =
2033 const unsigned int n_point_layers_toroidal =
2034 n_cells_toroidal + additional_layer;
2035 std::vector<Point<3>> vertices(8 * n_point_layers_toroidal);
2047 const double phi_cell = phi / n_cells_toroidal;
2048 for (
unsigned int c = 1; c < n_point_layers_toroidal; ++c)
2050 for (
unsigned int v = 0; v < 8; ++v)
2052 const double inner_radius_2d = vertices[v][0];
2053 vertices[8 * c + v][0] = inner_radius_2d *
std::cos(phi_cell * c);
2054 vertices[8 * c + v][1] = vertices[v][1];
2055 vertices[8 * c + v][2] = inner_radius_2d *
std::sin(phi_cell * c);
2060 std::vector<CellData<3>> cells(5 * n_cells_toroidal);
2061 for (
unsigned int c = 0; c < n_cells_toroidal; ++c)
2063 for (
unsigned int j = 0; j < 2; ++j)
2065 const unsigned int offset =
2066 (8 * (c + j)) % (8 * n_point_layers_toroidal);
2069 for (
unsigned int c2 = 0; c2 < 5; ++c2)
2070 for (
unsigned int i = 0; i < 4; ++i)
2071 cells[5 * c + c2].vertices[i + j * 4] =
2072 offset + circle_cell_vertices[c2][i];
2075 cells[5 * c].material_id = 0;
2077 cells[5 * c + 1].material_id = 0;
2078 cells[5 * c + 2].material_id = 1;
2079 cells[5 * c + 3].material_id = 0;
2080 cells[5 * c + 4].material_id = 0;
2088 for (
const auto &cell : tria.cell_iterators())
2091 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
2095 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2097 cell->face(f)->set_all_manifold_ids(1);
2102 if (cell->material_id() == 1)
2104 cell->set_all_manifold_ids(2);
2106 cell->set_material_id(0);
2123 template <
int dim,
int spacedim>
2138 cell->vertex(i) = vertices[i];
2144 "The volume of the cell is not greater than zero. "
2145 "This could be due to the wrong ordering of the vertices."));
2176 std::array<Tensor<1, 2>, 2> edges;
2177 edges[0] = corners[0];
2178 edges[1] = corners[1];
2179 std::vector<unsigned int> subdivisions;
2180 subdivided_parallelepiped<2, 2>(
2181 tria, origin, edges, subdivisions,
colorize);
2192 unsigned int n_subdivisions[dim];
2193 for (
unsigned int i = 0; i < dim; ++i)
2194 n_subdivisions[i] = 1;
2203 const unsigned int n_subdivisions,
2209 unsigned int n_subdivisions_[dim];
2210 for (
unsigned int i = 0; i < dim; ++i)
2211 n_subdivisions_[i] = n_subdivisions;
2221 const unsigned int (&n_subdivisions)[dim],
2223 const unsigned int *n_subdivisions,
2229 std::vector<unsigned int> subdivisions;
2230 std::array<Tensor<1, dim>, dim> edges;
2231 for (
unsigned int i = 0; i < dim; ++i)
2233 subdivisions.push_back(n_subdivisions[i]);
2234 edges[i] = corners[i];
2237 subdivided_parallelepiped<dim, dim>(
2238 tria, origin, edges, subdivisions,
colorize);
2244 template <
int dim,
int spacedim>
2249 const std::vector<unsigned int> &subdivisions,
2252 std::vector<unsigned int> compute_subdivisions = subdivisions;
2253 if (compute_subdivisions.empty())
2255 compute_subdivisions.resize(dim, 1);
2258 Assert(compute_subdivisions.size() == dim,
2259 ExcMessage(
"One subdivision must be provided for each dimension."));
2261 for (
unsigned int i = 0; i < dim; ++i)
2263 Assert(compute_subdivisions[i] > 0,
2264 ExcInvalidRepetitions(subdivisions[i]));
2266 edges[i].
norm() > 0,
2268 "Edges in subdivided_parallelepiped() must not be degenerate."));
2276 bool twisted_data =
false;
2281 twisted_data = (edges[0][0] < 0);
2288 const double plane_normal =
2289 edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2290 twisted_data = (plane_normal < 0.0);
2299 (edges[0].
norm() * edges[1].
norm()) -
2302 "Edges in subdivided_parallelepiped() must point in"
2303 " different directions."));
2305 cross_product_3d(edges[0], edges[1]);
2319 twisted_data = (plane_normal * edges[2] < 0.0);
2328 ExcInvalidInputOrientation(
2329 "The triangulation you are trying to create will consist of cells"
2330 " with negative measures. This is usually the result of input data"
2331 " that does not define a right-handed coordinate system. The usual"
2332 " fix for this is to ensure that in 1d the given point is to the"
2333 " right of the origin (or the given edge tensor is positive), in 2d"
2334 " that the two edges (and their cross product) obey the right-hand"
2335 " rule (which may usually be done by switching the order of the"
2336 " points or edge tensors), or in 3d that the edges form a"
2337 " right-handed coordinate system (which may also be accomplished by"
2338 " switching the order of the first two points or edge tensors)."));
2341 for (
unsigned int i = 0; i < dim; ++i)
2342 for (
unsigned int j = i + 1; j < dim; ++j)
2343 Assert((edges[i] != edges[j]),
2345 "Degenerate edges of subdivided_parallelepiped encountered."));
2348 std::vector<Point<spacedim>> points;
2353 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2354 points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2358 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2359 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2360 points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2361 edges[1] / compute_subdivisions[1] * y);
2365 for (
unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2366 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2367 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2368 points.push_back(origin +
2369 edges[0] / compute_subdivisions[0] * x +
2370 edges[1] / compute_subdivisions[1] * y +
2371 edges[2] / compute_subdivisions[2] * z);
2380 for (
unsigned int i = 0; i < dim; ++i)
2381 n_cells *= compute_subdivisions[i];
2382 std::vector<CellData<dim>> cells(n_cells);
2388 for (
unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2390 cells[x].vertices[0] = x;
2391 cells[x].vertices[1] = x + 1;
2394 cells[x].material_id = 0;
2401 const unsigned int n_dy = compute_subdivisions[1];
2402 const unsigned int n_dx = compute_subdivisions[0];
2404 for (
unsigned int y = 0; y < n_dy; ++y)
2405 for (
unsigned int x = 0; x < n_dx; ++x)
2407 const unsigned int c = y * n_dx + x;
2408 cells[c].vertices[0] = y * (n_dx + 1) + x;
2409 cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2410 cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2411 cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2414 cells[c].material_id = 0;
2422 const unsigned int n_dz = compute_subdivisions[2];
2423 const unsigned int n_dy = compute_subdivisions[1];
2424 const unsigned int n_dx = compute_subdivisions[0];
2426 for (
unsigned int z = 0; z < n_dz; ++z)
2427 for (
unsigned int y = 0; y < n_dy; ++y)
2428 for (
unsigned int x = 0; x < n_dx; ++x)
2430 const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2432 cells[c].vertices[0] =
2433 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2434 cells[c].vertices[1] =
2435 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2436 cells[c].vertices[2] =
2437 z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2438 cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2439 (y + 1) * (n_dx + 1) + x + 1;
2440 cells[c].vertices[4] =
2441 (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2442 cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2443 y * (n_dx + 1) + x + 1;
2444 cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2445 (y + 1) * (n_dx + 1) + x;
2446 cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2447 (y + 1) * (n_dx + 1) + x + 1;
2450 cells[c].material_id = 0;
2471 for (; cell != endc; ++cell)
2473 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
2475 if (cell->face(face)->at_boundary())
2476 cell->face(face)->set_boundary_id(face);
2483 template <
int dim,
int spacedim>
2486 const unsigned int repetitions,
2491 Assert(repetitions >= 1, ExcInvalidRepetitions(repetitions));
2493 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2496 for (
unsigned int i = 0; i < dim; ++i)
2502 std::vector<unsigned int> reps(dim, repetitions);
2508 template <
int dim,
int spacedim>
2511 const std::vector<unsigned int> &repetitions,
2516 Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim));
2522 for (
unsigned int i = 0; i < dim; ++i)
2529 std::array<Point<spacedim>, dim> delta;
2530 for (
unsigned int i = 0; i < dim; ++i)
2532 Assert(repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
2534 delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2538 "The first dim entries of coordinates of p1 and p2 need to be different."));
2542 std::vector<Point<spacedim>> points;
2546 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2547 points.push_back(p1 + x * delta[0]);
2551 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2552 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2553 points.push_back(p1 + x * delta[0] + y * delta[1]);
2557 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
2558 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2559 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2560 points.push_back(p1 + x * delta[0] + y * delta[1] +
2569 std::vector<CellData<dim>> cells;
2574 cells.resize(repetitions[0]);
2575 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2577 cells[x].vertices[0] = x;
2578 cells[x].vertices[1] = x + 1;
2579 cells[x].material_id = 0;
2586 cells.resize(repetitions[1] * repetitions[0]);
2587 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2588 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2590 const unsigned int c = x + y * repetitions[0];
2591 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2592 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2593 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2594 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2595 cells[c].material_id = 0;
2602 const unsigned int n_x = (repetitions[0] + 1);
2603 const unsigned int n_xy =
2604 (repetitions[0] + 1) * (repetitions[1] + 1);
2606 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2607 for (
unsigned int z = 0; z < repetitions[2]; ++z)
2608 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2609 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2611 const unsigned int c = x + y * repetitions[0] +
2612 z * repetitions[0] * repetitions[1];
2613 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2614 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2615 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2616 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2617 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2618 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2619 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2620 cells[c].vertices[7] =
2621 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2622 cells[c].material_id = 0;
2643 double epsilon = std::numeric_limits<double>::max();
2644 for (
unsigned int i = 0; i < dim; ++i)
2645 epsilon =
std::min(epsilon, 0.01 * delta[i][i]);
2648 "The distance between corner points must be positive."));
2652 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2661 const std::vector<std::vector<double>> &step_sz,
2666 Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
2676 std::vector<std::vector<double>> step_sizes(step_sz);
2678 for (
unsigned int i = 0; i < dim; ++i)
2682 std::swap(p1[i], p2[i]);
2683 std::reverse(step_sizes[i].
begin(), step_sizes[i].
end());
2689 for (
unsigned int j = 0; j < step_sizes.at(i).size(); ++j)
2690 x += step_sizes[i][j];
2691 Assert(std::fabs(x - (p2[i] - p1[i])) <= 1e-12 * std::fabs(x),
2693 "The sequence of step sizes in coordinate direction " +
2695 " must be equal to the distance of the two given "
2696 "points in this coordinate direction."));
2703 std::vector<Point<dim>> points;
2709 for (
unsigned int i = 0;; ++i)
2718 if (i == step_sizes[0].
size())
2721 x += step_sizes[0][i];
2729 for (
unsigned int j = 0;; ++j)
2732 for (
unsigned int i = 0;; ++i)
2734 points.push_back(
Point<dim>(p1[0] + x, p1[1] + y));
2735 if (i == step_sizes[0].
size())
2738 x += step_sizes[0][i];
2741 if (j == step_sizes[1].
size())
2744 y += step_sizes[1][j];
2751 for (
unsigned int k = 0;; ++k)
2754 for (
unsigned int j = 0;; ++j)
2757 for (
unsigned int i = 0;; ++i)
2760 Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2761 if (i == step_sizes[0].
size())
2764 x += step_sizes[0][i];
2767 if (j == step_sizes[1].
size())
2770 y += step_sizes[1][j];
2773 if (k == step_sizes[2].
size())
2776 z += step_sizes[2][k];
2787 std::vector<CellData<dim>> cells;
2792 cells.resize(step_sizes[0].
size());
2793 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2795 cells[x].vertices[0] = x;
2796 cells[x].vertices[1] = x + 1;
2797 cells[x].material_id = 0;
2804 cells.resize(step_sizes[1].
size() * step_sizes[0].
size());
2805 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2806 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2808 const unsigned int c = x + y * step_sizes[0].size();
2809 cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2810 cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2811 cells[c].vertices[2] =
2812 (y + 1) * (step_sizes[0].
size() + 1) + x;
2813 cells[c].vertices[3] =
2814 (y + 1) * (step_sizes[0].
size() + 1) + x + 1;
2815 cells[c].material_id = 0;
2822 const unsigned int n_x = (step_sizes[0].size() + 1);
2823 const unsigned int n_xy =
2824 (step_sizes[0].size() + 1) * (step_sizes[1].
size() + 1);
2826 cells.resize(step_sizes[2].
size() * step_sizes[1].size() *
2827 step_sizes[0].size());
2828 for (
unsigned int z = 0; z < step_sizes[2].size(); ++z)
2829 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2830 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2832 const unsigned int c =
2833 x + y * step_sizes[0].size() +
2834 z * step_sizes[0].size() * step_sizes[1].size();
2835 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2836 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2837 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2838 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2839 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2840 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2841 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2842 cells[c].vertices[7] =
2843 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2844 cells[c].material_id = 0;
2866 *std::min_element(step_sizes[0].
begin(), step_sizes[0].
end());
2867 for (
unsigned int i = 1; i < dim; ++i)
2869 *std::min_element(step_sizes[i].
begin(),
2870 step_sizes[i].
end()));
2871 const double epsilon = 0.01 * min_size;
2875 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2884 const std::vector<std::vector<double>> &spacing,
2889 Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
2893 Assert(spacing[0].
size() == n_cells, ExcInvalidRepetitionsDimension(1));
2895 double delta = std::numeric_limits<double>::max();
2896 for (
unsigned int i = 0; i <
n_cells; ++i)
2898 Assert(spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
2899 delta =
std::min(delta, spacing[0][i]);
2903 std::vector<Point<1>> points;
2905 for (
unsigned int x = 0; x <=
n_cells; ++x)
2907 points.emplace_back(ax);
2909 ax += spacing[0][x];
2912 unsigned int n_val_cells = 0;
2913 for (
unsigned int i = 0; i <
n_cells; ++i)
2917 std::vector<CellData<1>> cells(n_val_cells);
2918 unsigned int id = 0;
2919 for (
unsigned int x = 0; x <
n_cells; ++x)
2922 cells[id].vertices[0] = x;
2923 cells[id].vertices[1] = x + 1;
2942 const std::vector<std::vector<double>> &spacing,
2947 Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
2949 std::vector<unsigned int> repetitions(2);
2950 double delta = std::numeric_limits<double>::max();
2951 for (
unsigned int i = 0; i < 2; ++i)
2953 repetitions[i] = spacing[i].size();
2954 for (
unsigned int j = 0; j < repetitions[i]; ++j)
2956 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
2957 delta =
std::min(delta, spacing[i][j]);
2960 ExcInvalidRepetitionsDimension(i));
2964 std::vector<Point<2>> points;
2966 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2969 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2971 points.emplace_back(ax, ay);
2972 if (x < repetitions[0])
2973 ax += spacing[0][x];
2975 if (y < repetitions[1])
2976 ay += spacing[1][y];
2980 unsigned int n_val_cells = 0;
2981 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
2982 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
2986 std::vector<CellData<2>> cells(n_val_cells);
2987 unsigned int id = 0;
2988 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2989 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2992 cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
2993 cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2994 cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2995 cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3009 double eps = 0.01 * delta;
3011 for (; cell != endc; ++cell)
3013 Point<2> cell_center = cell->center();
3014 for (
const unsigned int f :
GeometryInfo<2>::face_indices())
3017 Point<2> face_center = cell->face(f)->center();
3018 for (
unsigned int i = 0; i < 2; ++i)
3020 if (face_center[i] < cell_center[i] - eps)
3021 cell->face(f)->set_boundary_id(i * 2);
3022 if (face_center[i] > cell_center[i] + eps)
3023 cell->face(f)->set_boundary_id(i * 2 + 1);
3034 const std::vector<std::vector<double>> &spacing,
3039 const unsigned int dim = 3;
3041 Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
3043 std::array<unsigned int, dim> repetitions;
3044 double delta = std::numeric_limits<double>::max();
3045 for (
unsigned int i = 0; i < dim; ++i)
3047 repetitions[i] = spacing[i].size();
3048 for (
unsigned int j = 0; j < repetitions[i]; ++j)
3050 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
3051 delta =
std::min(delta, spacing[i][j]);
3054 ExcInvalidRepetitionsDimension(i));
3058 std::vector<Point<dim>> points;
3060 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3063 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3066 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3068 points.emplace_back(ax, ay, az);
3069 if (x < repetitions[0])
3070 ax += spacing[0][x];
3072 if (y < repetitions[1])
3073 ay += spacing[1][y];
3075 if (z < repetitions[2])
3076 az += spacing[2][z];
3080 unsigned int n_val_cells = 0;
3081 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
3082 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
3083 for (
unsigned int k = 0; k <
material_id.size(2); ++k)
3087 std::vector<CellData<dim>> cells(n_val_cells);
3088 unsigned int id = 0;
3089 const unsigned int n_x = (repetitions[0] + 1);
3090 const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3091 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3092 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3093 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3096 cells[id].vertices[0] = z * n_xy + y * n_x + x;
3097 cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3098 cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3099 cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3100 cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3101 cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3102 cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3103 cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3117 double eps = 0.01 * delta;
3120 for (; cell != endc; ++cell)
3126 Point<dim> face_center = cell->face(f)->center();
3127 for (
unsigned int i = 0; i < dim; ++i)
3129 if (face_center[i] < cell_center[i] - eps)
3130 cell->face(f)->set_boundary_id(i * 2);
3131 if (face_center[i] > cell_center[i] + eps)
3132 cell->face(f)->set_boundary_id(i * 2 + 1);
3139 template <
int dim,
int spacedim>
3142 const std::vector<unsigned int> &holes)
3151 for (
unsigned int d = 0;
d < dim; ++
d)
3158 std::array<Point<spacedim>, dim> delta;
3159 std::array<unsigned int, dim> repetitions;
3160 for (
unsigned int i = 0; i < dim; ++i)
3163 ExcMessage(
"At least one hole needed in each direction"));
3164 repetitions[i] = 2 * holes[i] + 1;
3165 delta[i][i] = (p2[i] - p1[i]);
3170 std::vector<Point<spacedim>> points;
3174 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3175 points.push_back(p1 + x * delta[0]);
3179 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3180 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3181 points.push_back(p1 + x * delta[0] + y * delta[1]);
3185 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3186 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3187 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3188 points.push_back(p1 + x * delta[0] + y * delta[1] +
3198 std::vector<CellData<dim>> cells;
3203 cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3205 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3206 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3208 if ((x % 2 == 1) && (y % 2 == 1))
3211 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
3212 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3213 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3214 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3215 cells[c].material_id = 0;
3223 const unsigned int n_x = (repetitions[0] + 1);
3224 const unsigned int n_xy =
3225 (repetitions[0] + 1) * (repetitions[1] + 1);
3227 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3230 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3231 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3232 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3235 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3236 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3237 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3238 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3239 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3240 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3241 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3242 cells[c].vertices[7] =
3243 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3244 cells[c].material_id = 0;
3272 const unsigned int ,
3284 const unsigned int ,
3296 bool inline point_in_2d_box(
const Point<2> &p,
3298 const double radius)
3300 return (
std::abs(p[0] - c[0]) < radius) &&
3309 template <
int dim,
int spacedim>
3313 double length = std::numeric_limits<double>::max();
3314 for (
const auto &cell :
triangulation.active_cell_iterators())
3315 for (unsigned
int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3316 length =
std::min(length, cell->line(n)->diameter());
3326 const double inner_radius,
3327 const double outer_radius,
3328 const double pad_bottom,
3329 const double pad_top,
3330 const double pad_left,
3331 const double pad_right,
3336 const unsigned int ,
3339 const bool with_padding =
3340 pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3350 double length = std::numeric_limits<double>::max();
3351 for (
const auto &cell : tria.active_cell_iterators())
3352 for (unsigned
int n = 0; n < cell->n_lines(); ++n)
3353 length =
std::min(length, cell->line(n)->diameter());
3359 Triangulation<2> &cylinder_tria = with_padding ? cylinder_tria_maybe : tria;
3368 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3369 cell->set_manifold_id(tfi_manifold_id);
3371 const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
3372 const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
3378 auto add_sizes = [](std::vector<double> &step_sizes,
3379 const double padding,
3380 const double h) ->
void {
3383 const auto rounded =
3384 static_cast<unsigned int>(std::round(padding / h));
3387 const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3388 for (
unsigned int i = 0; i < num; ++i)
3389 step_sizes.push_back(padding / num);
3392 std::vector<std::vector<double>> step_sizes(2);
3395 add_sizes(step_sizes[0], pad_left, outer_radius);
3397 step_sizes[0].push_back(outer_radius);
3398 step_sizes[0].push_back(outer_radius);
3400 add_sizes(step_sizes[0], pad_right, outer_radius);
3403 add_sizes(step_sizes[1], pad_bottom, outer_radius);
3405 step_sizes[1].push_back(outer_radius);
3406 step_sizes[1].push_back(outer_radius);
3408 add_sizes(step_sizes[1], pad_top, outer_radius);
3413 bulk_tria, step_sizes, bl, tr,
colorize);
3416 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3417 for (
const auto &cell : bulk_tria.active_cell_iterators())
3418 if (
internal::point_in_2d_box(cell->center(), center, outer_radius))
3419 cells_to_remove.
insert(cell);
3423 bulk_tria, cells_to_remove, tria_without_cylinder);
3425 const double tolerance =
3426 std::min(min_line_length(tria_without_cylinder),
3427 min_line_length(cylinder_tria)) /
3431 tria_without_cylinder, cylinder_tria, tria, tolerance,
true);
3435 for (
const auto &cell : tria.active_cell_iterators())
3439 if (cell->manifold_id() == tfi_manifold_id)
3441 for (
const unsigned int face_n : cell->face_indices())
3443 const auto &face = cell->face(face_n);
3444 if (face->at_boundary() &&
3445 internal::point_in_2d_box(face->center(),
3447 outer_radius * (1. - 1e-12)))
3448 face->set_manifold_id(polar_manifold_id);
3450 face->set_manifold_id(tfi_manifold_id);
3461 static constexpr double tol =
3462 std::numeric_limits<double>::epsilon() * 10000;
3464 for (
const auto &cell : tria.active_cell_iterators())
3465 for (const unsigned
int face_n : cell->face_indices())
3467 const auto face = cell->face(face_n);
3468 if (face->at_boundary())
3470 const Point<2> center = face->center();
3473 face->set_boundary_id(0);
3476 face->set_boundary_id(1);
3479 face->set_boundary_id(2);
3482 face->set_boundary_id(3);
3486 Assert(cell->manifold_id() == tfi_manifold_id,
3488 face->set_boundary_id(4);
3509 const double inner_radius,
3510 const double outer_radius,
3511 const double pad_bottom,
3512 const double pad_top,
3513 const double pad_left,
3514 const double pad_right,
3519 const unsigned int n_slices,
3530 Point<2>(new_center[0], new_center[1]),
3553 tria.
set_manifold(polar_manifold_id, cylindrical_manifold);
3562 const double shell_region_width,
3563 const unsigned int n_shells,
3564 const double skewness,
3567 Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3568 ExcMessage(
"The width of the shell region must be less than 0.05 "
3569 "(and preferably close to 0.03)"));
3603 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3605 for (
const auto &cell : bulk_tria.active_cell_iterators())
3607 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3608 cells_to_remove.insert(cell);
3612 for (
const unsigned int vertex_n :
3614 if (cell->vertex(vertex_n) ==
Point<2>())
3618 cylinder_triangulation_offset =
3619 2.0 * (cell->vertex(3) -
Point<2>());
3626 bulk_tria, cells_to_remove, tria_without_cylinder);
3633 0.05 + shell_region_width,
3638 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3641 if (
std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1e-10)
3642 cell->vertex(vertex_n)[0] = -0.1;
3643 else if (
std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1e-10)
3644 cell->vertex(vertex_n)[0] = 0.1;
3648 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3650 cell->set_manifold_id(tfi_manifold_id);
3651 for (
const unsigned int face_n :
GeometryInfo<2>::face_indices())
3652 if (!cell->face(face_n)->at_boundary())
3653 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3655 if (0.0 < shell_region_width)
3658 ExcMessage(
"If the shell region has positive width then "
3659 "there must be at least one shell."));
3664 0.05 + shell_region_width,
3671 const double vertex_tolerance =
3672 std::min(internal::minimal_vertex_distance(shell_tria),
3673 internal::minimal_vertex_distance(cylinder_tria)) *
3679 shell_tria, cylinder_tria, temp, vertex_tolerance,
true);
3680 cylinder_tria = std::move(temp);
3686 const double vertex_tolerance =
3687 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3688 internal::minimal_vertex_distance(cylinder_tria)) /
3691 tria_without_cylinder, cylinder_tria, tria, vertex_tolerance,
true);
3704 const double shift =
3705 std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3706 for (
const auto &cell : tria.active_cell_iterators())
3708 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3709 cell->vertex(v) =
Point<2>(0.2 -
shift, 0.205);
3710 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1e-10)
3711 cell->vertex(v) =
Point<2>(0.2 + shift, 0.205);
3712 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1e-10)
3713 cell->vertex(v) =
Point<2>(0.2, 0.2 - shift);
3714 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1e-10)
3715 cell->vertex(v) =
Point<2>(0.2, 0.2 + shift);
3720 for (
const auto &cell : tria.active_cell_iterators())
3722 cell->set_all_manifold_ids(polar_manifold_id);
3726 for (
const auto &cell : tria.active_cell_iterators())
3733 std::vector<Point<2> *> cylinder_pointers;
3734 for (
const auto &face : tria.active_face_iterators())
3737 cylinder_pointers.push_back(&face->vertex(0));
3738 cylinder_pointers.push_back(&face->vertex(1));
3741 std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3742 cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3743 cylinder_pointers.end()),
3744 cylinder_pointers.end());
3748 for (
const Point<2> *
const ptr : cylinder_pointers)
3749 center += *ptr / double(cylinder_pointers.
size());
3752 for (
Point<2> *
const ptr : cylinder_pointers)
3753 *ptr +=
Point<2>(0.2, 0.2) - center;
3765 for (
const auto &face : tria.active_face_iterators())
3766 if (face->at_boundary())
3768 const Point<2> center = face->center();
3770 if (
std::abs(center[0] - 0.0) < 1e-10)
3771 face->set_boundary_id(0);
3773 else if (
std::abs(center[0] - 2.2) < 1e-10)
3774 face->set_boundary_id(1);
3776 else if (face->manifold_id() == polar_manifold_id)
3777 face->set_boundary_id(2);
3782 std::abs(center[1] - 0.41) < 1.0e-10,
3784 face->set_boundary_id(3);
3794 const double shell_region_width,
3795 const unsigned int n_shells,
3796 const double skewness,
3801 tria_2, shell_region_width, n_shells, skewness,
colorize);
3819 tria.
set_manifold(cylindrical_manifold_id, cylindrical_manifold);
3826 for (
const auto &face : tria.active_face_iterators())
3828 face->set_boundary_id(3);
3833 template <
int dim,
int spacedim>
3836 const std::vector<unsigned int> &sizes,
3846 for (
unsigned int d = 0;
d < dim; ++
d)
3849 std::vector<Point<spacedim>> points;
3851 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
3854 std::vector<CellData<dim>> cells(n_cells);
3859 for (
unsigned int d = 0;
d < dim; ++
d)
3860 p[d] = 0.5 * dimensions[d] *
3863 points.push_back(p);
3864 cells[0].vertices[i] = i;
3866 cells[0].material_id = 0;
3871 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
3877 for (
unsigned int j = 0; j < sizes[face]; ++j, ++
cell_index)
3879 const unsigned int last_cell = (j == 0) ? 0U : (
cell_index - 1);
3881 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
3884 const unsigned int cellv =
3886 const unsigned int ocellv =
3890 cells[last_cell].vertices[cellv];
3893 cells[
cell_index].vertices[cellv] = points.size();
3898 points.push_back(p);
4047 const double thickness,
4051 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4053 std::vector<Point<2>> vertices(16);
4055 coords[0] = left - thickness;
4058 coords[3] = right + thickness;
4061 for (
const double y : coords)
4062 for (const double x : coords)
4063 vertices[k++] =
Point<2>(x, y);
4067 std::vector<CellData<2>> cells(9);
4069 for (
unsigned int i0 = 0; i0 < 3; ++i0)
4070 for (
unsigned int i1 = 0; i1 < 3; ++i1)
4072 cells[k].vertices[0] = i1 + 4 * i0;
4073 cells[k].vertices[1] = i1 + 4 * i0 + 1;
4074 cells[k].vertices[2] = i1 + 4 * i0 + 4;
4075 cells[k].vertices[3] = i1 + 4 * i0 + 5;
4077 cells[k].material_id = materials[k];
4095 const double rl2 = (right + left) / 2;
4106 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4111 for (
unsigned int i = 0; i < 4; ++i)
4113 for (
unsigned int j = 0; j < 4; ++j)
4114 cells[i].vertices[j] = cell_vertices[i][j];
4115 cells[i].material_id = 0;
4118 std::end(vertices)),
4125 cell->face(1)->set_boundary_id(1);
4127 cell->face(0)->set_boundary_id(2);
4136 const double radius_0,
4137 const double radius_1,
4138 const double half_length)
4142 vertices_tmp[0] =
Point<2>(-half_length, -radius_0);
4143 vertices_tmp[1] =
Point<2>(half_length, -radius_1);
4144 vertices_tmp[2] =
Point<2>(-half_length, radius_0);
4145 vertices_tmp[3] =
Point<2>(half_length, radius_1);
4147 const std::vector<Point<2>> vertices(std::begin(vertices_tmp),
4148 std::end(vertices_tmp));
4152 cell_vertices[0][i] = i;
4157 cells[0].vertices[i] = cell_vertices[0][i];
4159 cells[0].material_id = 0;
4164 cell->face(0)->set_boundary_id(1);
4165 cell->face(1)->set_boundary_id(2);
4167 for (
unsigned int i = 2; i < 4; ++i)
4168 cell->face(i)->set_boundary_id(0);
4185 Point<2>((a + b) / 2, (a + b) / 2),
4189 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4193 for (
unsigned int i = 0; i < 3; ++i)
4195 for (
unsigned int j = 0; j < 4; ++j)
4196 cells[i].vertices[j] = cell_vertices[i][j];
4197 cells[i].material_id = 0;
4201 std::end(vertices)),
4209 cell->face(0)->set_boundary_id(0);
4210 cell->face(2)->set_boundary_id(1);
4213 cell->face(1)->set_boundary_id(2);
4214 cell->face(2)->set_boundary_id(1);
4215 cell->face(3)->set_boundary_id(3);
4218 cell->face(0)->set_boundary_id(0);
4219 cell->face(1)->set_boundary_id(4);
4220 cell->face(3)->set_boundary_id(5);
4226 template <
int dim,
int spacedim>
4229 const std::vector<unsigned int> &repetitions,
4232 const std::vector<int> &n_cells_to_remove)
4238 for (
unsigned int d = 0;
d < dim; ++
d)
4240 Assert(std::fabs(n_cells_to_remove[d]) <= repetitions[d],
4241 ExcMessage(
"Attempting to cut away too many cells."));
4251 std::array<double, dim> h;
4253 for (
unsigned int d = 0;
d < dim; ++
d)
4256 h[
d] = (top_right[
d] - bottom_left[
d]) / repetitions[d];
4258 if (n_cells_to_remove[d] >= 0)
4262 h[
d] * std::fabs(n_cells_to_remove[d]) + bottom_left[
d];
4267 cut_step[
d] = top_right[
d] - h[
d] * std::fabs(n_cells_to_remove[d]);
4273 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4275 for (
const auto &cell : rectangle.active_cell_iterators())
4277 bool remove_cell =
true;
4278 for (
unsigned int d = 0;
d < dim && remove_cell; ++
d)
4279 if ((n_cells_to_remove[d] > 0 && cell->center()[d] >= cut_step[d]) ||
4280 (n_cells_to_remove[d] < 0 && cell->center()[d] <= cut_step[d]))
4281 remove_cell =
false;
4283 cells_to_remove.insert(cell);
4293 template <
int dim,
int spacedim>
4297 const double radius,
4298 const bool internal_manifolds)
4300 if constexpr (dim == 2)
4302 const auto embed_point = [](
const double x,
4304 if constexpr (spacedim == 2)
4306 else if constexpr (spacedim == 3)
4315 const double a = 1. / (1 +
std::sqrt(2.0));
4317 p + embed_point(-1., -1.) * (radius /
std::sqrt(2.0)),
4318 p + embed_point(+1., -1.) * (radius /
std::sqrt(2.0)),
4319 p + embed_point(-1., -1.) * (radius /
std::sqrt(2.0) * a),
4320 p + embed_point(+1., -1.) * (radius /
std::sqrt(2.0) * a),
4321 p + embed_point(-1., +1.) * (radius /
std::sqrt(2.0) * a),
4322 p + embed_point(+1., +1.) * (radius /
std::sqrt(2.0) * a),
4323 p + embed_point(-1., +1.) * (radius /
std::sqrt(2.0)),
4324 p + embed_point(+1., +1.) * (radius /
std::sqrt(2.0))};
4328 for (
unsigned int i = 0; i < 5; ++i)
4330 for (
unsigned int j = 0; j < 4; ++j)
4331 cells[i].vertices[j] = circle_cell_vertices[i][j];
4332 cells[i].material_id = 0;
4337 std::begin(vertices), std::end(vertices)),
4342 if (internal_manifolds)
4347 else if constexpr (dim == 3)
4353 const unsigned int n_vertices = 16;
4354 const Point<3> vertices[n_vertices] = {
4379 const unsigned int n_cells = 7;
4380 const int cell_vertices[
n_cells][8] = {
4381 {0, 1, 4, 5, 3, 2, 7, 6},
4382 {8, 9, 12, 13, 0, 1, 4, 5},
4383 {9, 13, 1, 5, 10, 14, 2, 6},
4384 {11, 10, 3, 2, 15, 14, 7, 6},
4385 {8, 0, 12, 4, 11, 3, 15, 7},
4386 {8, 9, 0, 1, 11, 10, 3, 2},
4387 {12, 4, 13, 5, 15, 7, 14, 6}};
4389 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
4391 for (
unsigned int i = 0; i <
n_cells; ++i)
4394 cells[i].vertices[j] = cell_vertices[i][j];
4395 cells[i].material_id = 0;
4400 std::end(vertices)),
4405 if (internal_manifolds)
4416 template <
int spacedim>
4420 const double inner_radius,
4421 const double outer_radius,
4422 const unsigned int n_cells,
4425 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4439 const unsigned int N =
4440 (
n_cells == 0 ?
static_cast<unsigned int>(
4441 std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
4442 (outer_radius - inner_radius))) :
4451 std::vector<Point<spacedim>> vertices(2 * N);
4452 for (
unsigned int i = 0; i <
N; ++i)
4458 vertices[i] =
point * outer_radius;
4459 vertices[i +
N] = vertices[i] * (inner_radius / outer_radius);
4461 vertices[i] += center;
4462 vertices[i +
N] += center;
4467 for (
unsigned int i = 0; i <
N; ++i)
4469 cells[i].vertices[0] = i;
4470 cells[i].vertices[1] = (i + 1) % N;
4471 cells[i].vertices[2] =
N + i;
4472 cells[i].vertices[3] =
N + ((i + 1) % N);
4474 cells[i].material_id = 0;
4480 colorize_hyper_shell(tria, center, inner_radius, outer_radius);
4492 const double inner_radius,
4493 const double outer_radius,
4494 const unsigned int n_cells,
4497 hyper_shell_2D(tria, center, inner_radius, outer_radius, n_cells,
colorize);
4506 const double inner_radius,
4507 const double outer_radius,
4508 const unsigned int n_cells,
4511 hyper_shell_2D(tria, center, inner_radius, outer_radius, n_cells,
colorize);
4521 const double inner_radius,
4522 const double outer_radius,
4523 const unsigned int n_cells)
4526 tria, outer_center, inner_radius, outer_radius, n_cells,
true);
4530 outer_radius - inner_radius > outer_center.
distance(inner_center),
4532 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4536 std::set<Point<dim> *> vertices_to_move;
4538 for (
const auto &face : tria.active_face_iterators())
4540 for (unsigned
int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4541 vertices_to_move.insert(&face->vertex(v));
4543 const auto shift = inner_center - outer_center;
4544 for (
const auto &p : vertices_to_move)
4575 const double radius,
4576 const double half_length)
4578 Point<2> p1(-half_length, -radius);
4587 switch (f->boundary_id())
4590 f->set_boundary_id(1);
4593 f->set_boundary_id(2);
4596 f->set_boundary_id(0);
4634 const double radius)
4636 const unsigned int dim = 2;
4650 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4654 for (
unsigned int i = 0; i < 3; ++i)
4656 for (
unsigned int j = 0; j < 4; ++j)
4657 cells[i].vertices[j] = cell_vertices[i][j];
4658 cells[i].material_id = 0;
4662 std::end(vertices)),
4673 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
4675 if (cell->face(i)->boundary_id() ==
4681 if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius ||
4682 cell->face(i)->center()[1] < p[1] + 1.e-5 * radius)
4684 cell->face(i)->set_boundary_id(1);
4698 const double radius)
4703 const double a = 1. / (1 +
std::sqrt(2.0));
4714 const int cell_vertices[5][4] = {{0, 1, 2, 3},
4721 for (
unsigned int i = 0; i < 4; ++i)
4723 for (
unsigned int j = 0; j < 4; ++j)
4724 cells[i].vertices[j] = cell_vertices[i][j];
4725 cells[i].material_id = 0;
4729 std::end(vertices)),
4740 for (
const unsigned int i :
GeometryInfo<2>::face_indices())
4742 if (cell->face(i)->boundary_id() ==
4747 if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius)
4749 cell->face(i)->set_boundary_id(1);
4765 const double inner_radius,
4766 const double outer_radius,
4767 const unsigned int n_cells,
4770 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4783 const unsigned int N =
4784 (
n_cells == 0 ?
static_cast<unsigned int>(
4785 std::ceil((pi * (outer_radius + inner_radius) / 2) /
4786 (outer_radius - inner_radius))) :
4795 std::vector<Point<2>> vertices(2 * (N + 1));
4796 for (
unsigned int i = 0; i <=
N; ++i)
4806 std::
sin(pi * i /
N - pi / 2)) *
4808 vertices[i +
N + 1] = vertices[i] * (inner_radius / outer_radius);
4810 vertices[i] += center;
4811 vertices[i +
N + 1] += center;
4817 for (
unsigned int i = 0; i <
N; ++i)
4819 cells[i].vertices[0] = i;
4820 cells[i].vertices[1] = (i + 1) % (N + 1);
4821 cells[i].vertices[2] =
N + 1 + i;
4822 cells[i].vertices[3] =
N + 1 + ((i + 1) % (N + 1));
4824 cells[i].material_id = 0;
4832 for (; cell != tria.
end(); ++cell)
4834 cell->face(2)->set_boundary_id(1);
4836 tria.
begin()->face(0)->set_boundary_id(3);
4838 tria.
last()->face(1)->set_boundary_id(2);
4849 const double inner_radius,
4850 const double outer_radius,
4851 const unsigned int n_cells,
4854 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4867 const unsigned int N =
4868 (
n_cells == 0 ?
static_cast<unsigned int>(
4869 std::ceil((pi * (outer_radius + inner_radius) / 4) /
4870 (outer_radius - inner_radius))) :
4879 std::vector<Point<2>> vertices(2 * (N + 1));
4880 for (
unsigned int i = 0; i <=
N; ++i)
4890 vertices[i +
N + 1] = vertices[i] * (inner_radius / outer_radius);
4892 vertices[i] += center;
4893 vertices[i +
N + 1] += center;
4899 for (
unsigned int i = 0; i <
N; ++i)
4901 cells[i].vertices[0] = i;
4902 cells[i].vertices[1] = (i + 1) % (N + 1);
4903 cells[i].vertices[2] =
N + 1 + i;
4904 cells[i].vertices[3] =
N + 1 + ((i + 1) % (N + 1));
4906 cells[i].material_id = 0;
4914 for (; cell != tria.
end(); ++cell)
4916 cell->face(2)->set_boundary_id(1);
4918 tria.
begin()->face(0)->set_boundary_id(3);
4920 tria.
last()->face(1)->set_boundary_id(2);
4937 const double rl2 = (right + left) / 2;
4938 const double len = (right - left) / 2.;
4951 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4952 {9, 4, 2, 5, 19, 14, 12, 15},
4953 {3, 2, 7, 6, 13, 12, 17, 16},
4954 {2, 5, 6, 8, 12, 15, 16, 18}};
4956 for (
unsigned int i = 0; i < 4; ++i)
4958 for (
unsigned int j = 0; j < 8; ++j)
4959 cells[i].vertices[j] = cell_vertices[i][j];
4960 cells[i].material_id = 0;
4963 std::end(vertices)),
4970 cell->face(1)->set_boundary_id(1);
4972 cell->face(0)->set_boundary_id(2);
4984 const double thickness,
4988 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4990 std::vector<Point<3>> vertices(64);
4992 coords[0] = left - thickness;
4995 coords[3] = right + thickness;
4998 for (
const double z : coords)
4999 for (const double y : coords)
5000 for (const double x : coords)
5001 vertices[k++] =
Point<3>(x, y, z);
5004 24, 26, 5, 4, 6, 1, 0,
5005 2, 9, 8, 10, 37, 36, 38,
5006 33, 32, 34, 41, 40, 42};
5008 std::vector<CellData<3>> cells(27);
5010 for (
unsigned int z = 0; z < 3; ++z)
5011 for (
unsigned int y = 0; y < 3; ++y)
5012 for (
unsigned int x = 0; x < 3; ++x)
5014 cells[k].vertices[0] = x + 4 * y + 16 * z;
5015 cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
5016 cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
5017 cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
5018 cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
5019 cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
5020 cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
5021 cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
5023 cells[k].material_id = materials[k];
5036 const double radius_0,
5037 const double radius_1,
5038 const double half_length)
5041 ExcMessage(
"The output triangulation object needs to be empty."));
5046 const auto n_slices = 1 +
static_cast<unsigned int>(std::ceil(
5047 half_length /
std::max(radius_0, radius_1)));
5060 auto shift_radii = [=](
const Point<3> &p) {
5061 const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
5062 const double factor = slope * (p[0] - -half_length) + 1.0;
5063 return Point<3>(p[0], factor * p[1], factor * p[2]);
5069 for (
const auto &face :
triangulation.active_face_iterators())
5070 if (face->at_boundary())
5072 if (
std::abs(face->center()[0] - -half_length) < 1e-8 * half_length)
5073 face->set_boundary_id(1);
5074 else if (
std::abs(face->center()[0] - half_length) <
5076 face->set_boundary_id(2);
5078 face->set_all_manifold_ids(0);
5101 Point<3>((a + b) / 2, a, (a + b) / 2),
5108 Point<3>((a + b) / 2, (a + b) / 2, a),
5110 Point<3>(a, (a + b) / 2, (a + b) / 2),
5111 Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
5112 Point<3>(b, (a + b) / 2, (a + b) / 2),
5114 Point<3>((a + b) / 2, (a + b) / 2, b),
5122 Point<3>((a + b) / 2, b, (a + b) / 2),
5126 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5127 {1, 2, 10, 11, 4, 5, 13, 14},
5128 {3, 4, 12, 13, 6, 7, 15, 16},
5129 {4, 5, 13, 14, 7, 8, 16, 17},
5130 {9, 10, 18, 19, 12, 13, 21, 22},
5131 {10, 11, 19, 20, 13, 14, 22, 23},
5132 {12, 13, 21, 22, 15, 16, 24, 25}};
5136 for (
unsigned int i = 0; i < 7; ++i)
5138 for (
unsigned int j = 0; j < 8; ++j)
5139 cells[i].vertices[j] = cell_vertices[i][j];
5140 cells[i].material_id = 0;
5144 std::end(vertices)),
5158 const unsigned int n_rotate_middle_square)
5161 ExcMessage(
"The number of rotation by pi/2 of the right square "
5162 "must be in the half-open range [0,4)."));
5164 constexpr unsigned int dim = 2;
5166 const unsigned int n_cells = 5;
5167 std::vector<CellData<dim>> cells(n_cells);
5170 const std::vector<Point<dim>> vertices = {
Point<dim>(0, 0),
5185 unsigned int cell_vertices[
n_cells][4] = {{0, 1, 2, 3},
5191 switch (n_rotate_middle_square)
5195 cell_vertices[1][0] = 4;
5196 cell_vertices[1][1] = 5;
5197 cell_vertices[1][2] = 1;
5198 cell_vertices[1][3] = 3;
5204 cell_vertices[1][0] = 5;
5205 cell_vertices[1][1] = 3;
5206 cell_vertices[1][2] = 4;
5207 cell_vertices[1][3] = 1;
5213 cell_vertices[1][0] = 3;
5214 cell_vertices[1][1] = 1;
5215 cell_vertices[1][2] = 5;
5216 cell_vertices[1][3] = 4;
5228 for (
const unsigned int vertex_index :
5243 const bool face_orientation,
5244 const bool face_flip,
5245 const bool face_rotation,
5246 const bool manipulate_left_cube)
5248 constexpr unsigned int dim = 3;
5250 const unsigned int n_cells = 2;
5251 std::vector<CellData<dim>> cells(n_cells);
5254 const std::vector<Point<dim>> vertices = {
Point<dim>(0, 0, 0),
5267 unsigned int cell_vertices[
n_cells][8] = {
5268 {0, 1, 2, 3, 4, 5, 6, 7},
5269 {1, 8, 3, 9, 5, 10, 7, 11}};
5272 const unsigned int this_case = 4 *
static_cast<int>(face_orientation) +
5273 2 *
static_cast<int>(face_flip) +
5274 static_cast<int>(face_rotation);
5276 if (manipulate_left_cube)
5282 cell_vertices[0][0] = 1;
5283 cell_vertices[0][1] = 0;
5284 cell_vertices[0][2] = 5;
5285 cell_vertices[0][3] = 4;
5286 cell_vertices[0][4] = 3;
5287 cell_vertices[0][5] = 2;
5288 cell_vertices[0][6] = 7;
5289 cell_vertices[0][7] = 6;
5295 cell_vertices[0][0] = 5;
5296 cell_vertices[0][1] = 4;
5297 cell_vertices[0][2] = 7;
5298 cell_vertices[0][3] = 6;
5299 cell_vertices[0][4] = 1;
5300 cell_vertices[0][5] = 0;
5301 cell_vertices[0][6] = 3;
5302 cell_vertices[0][7] = 2;
5308 cell_vertices[0][0] = 7;
5309 cell_vertices[0][1] = 6;
5310 cell_vertices[0][2] = 3;
5311 cell_vertices[0][3] = 2;
5312 cell_vertices[0][4] = 5;
5313 cell_vertices[0][5] = 4;
5314 cell_vertices[0][6] = 1;
5315 cell_vertices[0][7] = 0;
5320 cell_vertices[0][0] = 3;
5321 cell_vertices[0][1] = 2;
5322 cell_vertices[0][2] = 1;
5323 cell_vertices[0][3] = 0;
5324 cell_vertices[0][4] = 7;
5325 cell_vertices[0][5] = 6;
5326 cell_vertices[0][6] = 5;
5327 cell_vertices[0][7] = 4;
5333 cell_vertices[0][0] = 0;
5334 cell_vertices[0][1] = 1;
5335 cell_vertices[0][2] = 2;
5336 cell_vertices[0][3] = 3;
5337 cell_vertices[0][4] = 4;
5338 cell_vertices[0][5] = 5;
5339 cell_vertices[0][6] = 6;
5340 cell_vertices[0][7] = 7;
5346 cell_vertices[0][0] = 2;
5347 cell_vertices[0][1] = 3;
5348 cell_vertices[0][2] = 6;
5349 cell_vertices[0][3] = 7;
5350 cell_vertices[0][4] = 0;
5351 cell_vertices[0][5] = 1;
5352 cell_vertices[0][6] = 4;
5353 cell_vertices[0][7] = 5;
5359 cell_vertices[0][0] = 6;
5360 cell_vertices[0][1] = 7;
5361 cell_vertices[0][2] = 4;
5362 cell_vertices[0][3] = 5;
5363 cell_vertices[0][4] = 2;
5364 cell_vertices[0][5] = 3;
5365 cell_vertices[0][6] = 0;
5366 cell_vertices[0][7] = 1;
5372 cell_vertices[0][0] = 4;
5373 cell_vertices[0][1] = 5;
5374 cell_vertices[0][2] = 0;
5375 cell_vertices[0][3] = 1;
5376 cell_vertices[0][4] = 6;
5377 cell_vertices[0][5] = 7;
5378 cell_vertices[0][6] = 2;
5379 cell_vertices[0][7] = 3;
5390 cell_vertices[1][0] = 8;
5391 cell_vertices[1][1] = 1;
5392 cell_vertices[1][2] = 10;
5393 cell_vertices[1][3] = 5;
5394 cell_vertices[1][4] = 9;
5395 cell_vertices[1][5] = 3;
5396 cell_vertices[1][6] = 11;
5397 cell_vertices[1][7] = 7;
5403 cell_vertices[1][0] = 10;
5404 cell_vertices[1][1] = 5;
5405 cell_vertices[1][2] = 11;
5406 cell_vertices[1][3] = 7;
5407 cell_vertices[1][4] = 8;
5408 cell_vertices[1][5] = 1;
5409 cell_vertices[1][6] = 9;
5410 cell_vertices[1][7] = 3;
5416 cell_vertices[1][0] = 11;
5417 cell_vertices[1][1] = 7;
5418 cell_vertices[1][2] = 9;
5419 cell_vertices[1][3] = 3;
5420 cell_vertices[1][4] = 10;
5421 cell_vertices[1][5] = 5;
5422 cell_vertices[1][6] = 8;
5423 cell_vertices[1][7] = 1;
5429 cell_vertices[1][0] = 9;
5430 cell_vertices[1][1] = 3;
5431 cell_vertices[1][2] = 8;
5432 cell_vertices[1][3] = 1;
5433 cell_vertices[1][4] = 11;
5434 cell_vertices[1][5] = 7;
5435 cell_vertices[1][6] = 10;
5436 cell_vertices[1][7] = 5;
5442 cell_vertices[1][0] = 1;
5443 cell_vertices[1][1] = 8;
5444 cell_vertices[1][2] = 3;
5445 cell_vertices[1][3] = 9;
5446 cell_vertices[1][4] = 5;
5447 cell_vertices[1][5] = 10;
5448 cell_vertices[1][6] = 7;
5449 cell_vertices[1][7] = 11;
5455 cell_vertices[1][0] = 5;
5456 cell_vertices[1][1] = 10;
5457 cell_vertices[1][2] = 1;
5458 cell_vertices[1][3] = 8;
5459 cell_vertices[1][4] = 7;
5460 cell_vertices[1][5] = 11;
5461 cell_vertices[1][6] = 3;
5462 cell_vertices[1][7] = 9;
5468 cell_vertices[1][0] = 7;
5469 cell_vertices[1][1] = 11;
5470 cell_vertices[1][2] = 5;
5471 cell_vertices[1][3] = 10;
5472 cell_vertices[1][4] = 3;
5473 cell_vertices[1][5] = 9;
5474 cell_vertices[1][6] = 1;
5475 cell_vertices[1][7] = 8;
5481 cell_vertices[1][0] = 3;
5482 cell_vertices[1][1] = 9;
5483 cell_vertices[1][2] = 7;
5484 cell_vertices[1][3] = 11;
5485 cell_vertices[1][4] = 1;
5486 cell_vertices[1][5] = 8;
5487 cell_vertices[1][6] = 5;
5488 cell_vertices[1][7] = 10;
5498 for (
const unsigned int vertex_index :
5512 template <
int spacedim>
5516 const double radius)
5520 const std::set<types::boundary_id> boundary_ids = {0};
5532 const unsigned int x_subdivisions,
5533 const double radius,
5534 const double half_length)
5541 std::vector<Point<3>> vertices;
5542 const double initial_height = -half_length;
5543 const double height_increment = 2. * half_length / x_subdivisions;
5545 for (
unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep)
5547 const double height = initial_height + height_increment * rep;
5549 vertices.emplace_back(-d, height, -d);
5550 vertices.emplace_back(d, height, -d);
5551 vertices.emplace_back(-a, height, -a);
5552 vertices.emplace_back(a, height, -a);
5553 vertices.emplace_back(-a, height, a);
5554 vertices.emplace_back(a, height, a);
5555 vertices.emplace_back(-d, height, d);
5556 vertices.emplace_back(d, height, d);
5560 for (
auto &vertex : vertices)
5562 const double h = vertex[1];
5563 vertex[1] = -vertex[0];
5567 std::vector<std::vector<int>> cell_vertices;
5568 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5569 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5570 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5571 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5572 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5574 for (
unsigned int rep = 1; rep < x_subdivisions; ++rep)
5576 for (
unsigned int i = 0; i < 5; ++i)
5578 std::vector<int> new_cell_vertices(8);
5579 for (
unsigned int j = 0; j < 8; ++j)
5580 new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep;
5581 cell_vertices.push_back(new_cell_vertices);
5585 unsigned int n_cells = x_subdivisions * 5;
5587 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5589 for (
unsigned int i = 0; i <
n_cells; ++i)
5591 for (
unsigned int j = 0; j < 8; ++j)
5592 cells[i].vertices[j] = cell_vertices[i][j];
5593 cells[i].material_id = 0;
5597 std::end(vertices)),
5615 const double tolerance = 1
e-5 *
std::min(radius, half_length);
5617 for (
const auto &cell : tria.cell_iterators())
5619 if (cell->at_boundary(i))
5621 if (cell->face(i)->center()[0] > half_length - tolerance)
5623 cell->face(i)->set_boundary_id(2);
5626 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5628 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5629 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5630 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5631 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5633 cell->face(i)->line(e)->set_boundary_id(2);
5634 cell->face(i)->line(e)->set_manifold_id(
5638 else if (cell->face(i)->center()[0] < -half_length + tolerance)
5640 cell->face(i)->set_boundary_id(1);
5643 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5645 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5646 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5647 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5648 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5650 cell->face(i)->line(e)->set_boundary_id(1);
5651 cell->face(i)->line(e)->set_manifold_id(
5663 const double radius,
5664 const double half_length)
5673 const double radius)
5675 const unsigned int dim = 3;
5683 const double a = 0.528;
5684 const double b = 0.4533;
5685 const double c = 0.3752;
5689 center +
Point<dim>(+1, 0, 0) * (radius * a),
5690 center +
Point<dim>(0, +1, 0) * (radius * a),
5698 center +
Point<dim>(+1, +1, 1) * (radius * c),
5702 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5703 {1, 6, 2, 4, 8, 13, 9, 11},
5704 {5, 3, 6, 4, 12, 10, 13, 11},
5705 {7, 9, 10, 11, 14, 8, 12, 13}};
5709 for (
unsigned int i = 0; i < 4; ++i)
5711 for (
unsigned int j = 0; j < 8; ++j)
5712 cells[i].vertices[j] = cell_vertices[i][j];
5713 cells[i].material_id = 0;
5717 std::end(vertices)),
5727 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
5729 if (cell->face(i)->boundary_id() ==
5734 if (cell->face(i)->center()[0] < center[0] + 1.e-5 * radius ||
5735 cell->face(i)->center()[1] < center[1] + 1.e-5 * radius ||
5736 cell->face(i)->center()[2] < center[2] + 1.e-5 * radius)
5738 cell->face(i)->set_boundary_id(1);
5742 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5745 const Point<3> line_vertices[2] = {
5746 cell->face(i)->line(j)->vertex(0),
5747 cell->face(i)->line(j)->vertex(1)};
5748 if ((std::fabs(line_vertices[0].distance(center) - radius) >
5750 (std::fabs(line_vertices[1].distance(center) - radius) >
5753 cell->face(i)->line(j)->set_boundary_id(1);
5754 cell->face(i)->line(j)->set_manifold_id(
5772 const double radius)
5778 const double b = a / 2.0;
5779 const double c =
d / 2.0;
5781 const double hb = radius *
std::sqrt(3.0) / 4.0;
5782 const double hc = radius *
std::sqrt(3.0) / 2.0;
5804 int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5805 {0, 2, 8, 10, 6, 4, 14, 12},
5806 {2, 3, 10, 11, 4, 5, 12, 13},
5807 {1, 7, 9, 15, 3, 5, 11, 13},
5808 {6, 4, 14, 12, 7, 5, 15, 13},
5809 {8, 10, 9, 11, 14, 12, 15, 13}};
5813 for (
unsigned int i = 0; i < 6; ++i)
5815 for (
unsigned int j = 0; j < 8; ++j)
5816 cells[i].vertices[j] = cell_vertices[i][j];
5817 cells[i].material_id = 0;
5821 std::end(vertices)),
5837 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
5839 if (!cell->at_boundary(i))
5845 if (cell->face(i)->center()[0] < center[0] + 1.e-5 * radius)
5847 cell->face(i)->set_boundary_id(1);
5849 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5852 const Point<3> line_vertices[2] = {
5853 cell->face(i)->line(j)->vertex(0),
5854 cell->face(i)->line(j)->vertex(1)};
5855 if ((std::fabs(line_vertices[0].distance(center) - radius) >
5857 (std::fabs(line_vertices[1].distance(center) - radius) >
5860 cell->face(i)->line(j)->set_boundary_id(1);
5861 cell->face(i)->line(j)->set_manifold_id(
5878 const double radius)
5891 for (
unsigned int round = 0; round < dim; ++round)
5896 std::vector<Point<dim>> new_points(tria_copy.
n_vertices());
5898 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5906 else if (round == 1)
5908 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5917 else if (round == 2)
5918 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5931 std::vector<CellData<dim>> cells;
5932 cells.reserve(tria_copy.
n_cells());
5933 for (
const auto &cell : tria_copy.cell_iterators())
5937 data.vertices[v] = cell->vertex_index(v);
5938 data.material_id = cell->material_id();
5939 data.manifold_id = cell->manifold_id();
5940 cells.push_back(
data);
5948 if (round == dim - 1)
5957 for (
const auto &cell : tria.cell_iterators())
5958 if (cell->center().norm_square() > 0.4 * radius)
5959 cell->set_manifold_id(1);
5980 const double inner_radius,
5981 const double outer_radius)
5983 std::vector<Point<3>> vertices;
5984 std::vector<CellData<3>> cells;
5986 const double irad = inner_radius /
std::sqrt(3.0);
5987 const double orad = outer_radius /
std::sqrt(3.0);
5990 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
6000 vertices.reserve(8);
6001 for (
unsigned int i = 0; i < 8; ++i)
6002 vertices.push_back(p + hexahedron[i] * irad);
6003 for (
unsigned int i = 0; i < 8; ++i)
6004 vertices.push_back(p + hexahedron[i] * orad);
6006 const unsigned int n_cells = 6;
6007 const int cell_vertices[
n_cells][8] = {
6008 {8, 9, 10, 11, 0, 1, 2, 3},
6009 {9, 11, 1, 3, 13, 15, 5, 7},
6010 {12, 13, 4, 5, 14, 15, 6, 7},
6011 {8, 0, 10, 2, 12, 4, 14, 6},
6012 {8, 9, 0, 1, 12, 13, 4, 5},
6013 {10, 2, 11, 3, 14, 6, 15, 7}};
6017 for (
unsigned int i = 0; i <
n_cells; ++i)
6020 cells[i].vertices[j] = cell_vertices[i][j];
6021 cells[i].material_id = 0;
6032 const double inner_radius,
6033 const double outer_radius)
6035 std::vector<Point<3>> vertices;
6036 std::vector<CellData<3>> cells;
6038 const double irad = inner_radius /
std::sqrt(3.0);
6039 const double orad = outer_radius /
std::sqrt(3.0);
6045 static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0},
6053 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
6062 vertices.reserve(8);
6063 for (
unsigned int i = 0; i < 8; ++i)
6064 vertices.push_back(p + hexahedron[i] * irad);
6065 for (
unsigned int i = 0; i < 6; ++i)
6066 vertices.push_back(p + octahedron[i] * inner_radius);
6067 for (
unsigned int i = 0; i < 8; ++i)
6068 vertices.push_back(p + hexahedron[i] * orad);
6069 for (
unsigned int i = 0; i < 6; ++i)
6070 vertices.push_back(p + octahedron[i] * outer_radius);
6072 const unsigned int n_cells = 12;
6073 const unsigned int rhombi[
n_cells][4] = {{10, 4, 0, 8},
6088 for (
unsigned int i = 0; i <
n_cells; ++i)
6090 for (
unsigned int j = 0; j < 4; ++j)
6092 cells[i].vertices[j] = rhombi[i][j];
6093 cells[i].vertices[j + 4] = rhombi[i][j] + 14;
6095 cells[i].material_id = 0;
6105 const unsigned int n,
6106 const unsigned int n_refinement_steps,
6108 const double inner_radius,
6109 const double outer_radius)
6122 const unsigned int outer_radius_factor = 1 << n_refinement_steps;
6127 outer_radius_factor * outer_radius -
6128 (outer_radius_factor - 1) * inner_radius);
6133 outer_radius_factor * outer_radius -
6134 (outer_radius_factor - 1) * inner_radius);
6137 for (
unsigned int r = 0; r < n_refinement_steps; ++r)
6140 std::set<Triangulation<3>::active_cell_iterator> cells_to_remove;
6144 for (
const auto &cell : tmp.active_cell_iterators())
6146 unsigned int n_vertices_inside = 0;
6148 if ((cell->vertex(v) - p).norm_square() <
6149 inner_radius * inner_radius * (1 + 1
e-12))
6150 ++n_vertices_inside;
6151 if (n_vertices_inside < 4)
6152 cells_to_remove.insert(cell);
6156 if (r == n_refinement_steps - 1)
6166 tmp = std::move(copy);
6184 const double inner_radius,
6185 const double outer_radius,
6186 const unsigned int n_cells,
6189 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6192 unsigned int n_refinement_steps = 0;
6193 unsigned int n_cells_coarsened =
n_cells;
6194 if (n_cells != 96 && n_cells > 12)
6195 while (n_cells_coarsened > 12 && n_cells_coarsened % 4 == 0)
6197 ++n_refinement_steps;
6198 n_cells_coarsened /= 4;
6200 Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 ||
6201 (n_refinement_steps > 0 &&
6202 (n_cells_coarsened == 6 || n_cells_coarsened == 12)),
6203 ExcMessage(
"Invalid number of coarse mesh cells"));
6205 const unsigned int n = n_refinement_steps > 0 ?
6206 4 * n_cells_coarsened :
6212 internal::hyper_shell_6(tria, p, inner_radius, outer_radius);
6215 internal::hyper_shell_12(tria, p, inner_radius, outer_radius);
6219 internal::hyper_shell_24_48(
6220 tria, n, n_refinement_steps, p, inner_radius, outer_radius);
6229 internal::hyper_shell_12(tmp, p, inner_radius, outer_radius);
6246 colorize_hyper_shell(tria, p, inner_radius, outer_radius);
6256 const double inner_radius,
6257 const double outer_radius,
6258 const unsigned int ,
6261 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6265 const double d = outer_radius /
std::sqrt(2.0);
6266 const double a = inner_radius /
std::sqrt(2.0);
6268 const double b = a / 2.0;
6269 const double c =
d / 2.0;
6271 const double hb = inner_radius *
std::sqrt(3.0) / 2.0;
6272 const double hc = outer_radius *
std::sqrt(3.0) / 2.0;
6294 int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
6295 {0, 2, 8, 10, 6, 4, 14, 12},
6296 {1, 7, 9, 15, 3, 5, 11, 13},
6297 {6, 4, 14, 12, 7, 5, 15, 13},
6298 {8, 10, 9, 11, 14, 12, 15, 13}};
6302 for (
unsigned int i = 0; i < 5; ++i)
6304 for (
unsigned int j = 0; j < 8; ++j)
6305 cells[i].vertices[j] = cell_vertices[i][j];
6306 cells[i].material_id = 0;
6310 std::end(vertices)),
6320 for (; cell != tria.
end(); ++cell)
6321 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6322 if (cell->at_boundary(i))
6323 cell->face(i)->set_all_boundary_ids(2);
6329 for (cell = tria.
begin(); cell != tria.
end(); ++cell)
6330 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6331 if (cell->at_boundary(i))
6335 const Point<3> face_center(face->center());
6336 if (
std::abs(face_center[0] - center[0]) >
6337 1.e-6 * face_center.
norm())
6339 if (
std::abs((face_center - center).
norm() - inner_radius) <
6340 std::abs((face_center - center).
norm() - outer_radius))
6341 face->set_all_boundary_ids(0);
6343 face->set_all_boundary_ids(1);
6357 const double inner_radius,
6358 const double outer_radius,
6359 const unsigned int n,
6362 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6364 if (n == 0 || n == 3)
6366 const double a = inner_radius *
std::sqrt(2.0) / 2e0;
6367 const double b = outer_radius *
std::sqrt(2.0) / 2e0;
6368 const double c = a *
std::sqrt(3.0) / 2e0;
6370 const double e = outer_radius / 2e0;
6371 const double h = inner_radius / 2e0;
6373 std::vector<Point<3>> vertices;
6375 vertices.push_back(center +
Point<3>(0, inner_radius, 0));
6376 vertices.push_back(center +
Point<3>(a, a, 0));
6377 vertices.push_back(center +
Point<3>(b, b, 0));
6378 vertices.push_back(center +
Point<3>(0, outer_radius, 0));
6379 vertices.push_back(center +
Point<3>(0, a, a));
6380 vertices.push_back(center +
Point<3>(c, c, h));
6381 vertices.push_back(center +
Point<3>(d, d, e));
6382 vertices.push_back(center +
Point<3>(0, b, b));
6383 vertices.push_back(center +
Point<3>(inner_radius, 0, 0));
6384 vertices.push_back(center +
Point<3>(outer_radius, 0, 0));
6385 vertices.push_back(center +
Point<3>(a, 0, a));
6386 vertices.push_back(center +
Point<3>(b, 0, b));
6387 vertices.push_back(center +
Point<3>(0, 0, inner_radius));
6388 vertices.push_back(center +
Point<3>(0, 0, outer_radius));
6390 const int cell_vertices[3][8] = {
6391 {0, 1, 3, 2, 4, 5, 7, 6},
6392 {1, 8, 2, 9, 5, 10, 6, 11},
6393 {4, 5, 7, 6, 12, 10, 13, 11},
6395 std::vector<CellData<3>> cells(3);
6397 for (
unsigned int i = 0; i < 3; ++i)
6399 for (
unsigned int j = 0; j < 8; ++j)
6400 cells[i].vertices[j] = cell_vertices[i][j];
6401 cells[i].material_id = 0;
6414 colorize_quarter_hyper_shell(tria, center, inner_radius, outer_radius);
6425 const double length,
6426 const double inner_radius,
6427 const double outer_radius,
6428 const unsigned int n_radial_cells,
6429 const unsigned int n_axial_cells,
6432 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6446 const unsigned int N_r =
6447 (n_radial_cells == 0 ?
static_cast<unsigned int>(std::ceil(
6448 (2 * pi * (outer_radius + inner_radius) / 2) /
6449 (outer_radius - inner_radius))) :
6451 const unsigned int N_z =
6452 (n_axial_cells == 0 ?
6453 static_cast<unsigned int>(std::ceil(
6454 length / (2 * pi * (outer_radius + inner_radius) / 2 / N_r))) :
6463 std::vector<Point<2>> vertices_2d(2 * N_r);
6464 for (
unsigned int i = 0; i < N_r; ++i)
6469 vertices_2d[i + N_r] = vertices_2d[i] * (inner_radius / outer_radius);
6472 std::vector<Point<3>> vertices_3d;
6473 vertices_3d.reserve(2 * N_r * (N_z + 1));
6474 for (
unsigned int j = 0; j <= N_z; ++j)
6475 for (
unsigned int i = 0; i < 2 * N_r; ++i)
6477 const Point<3> v(vertices_2d[i][0],
6480 vertices_3d.push_back(v);
6483 std::vector<CellData<3>> cells(N_r * N_z,
CellData<3>());
6485 for (
unsigned int j = 0; j < N_z; ++j)
6486 for (
unsigned int i = 0; i < N_r; ++i)
6488 cells[i + j * N_r].vertices[0] = i + (j + 1) * 2 * N_r;
6489 cells[i + j * N_r].vertices[1] = (i + 1) % N_r + (j + 1) * 2 * N_r;
6490 cells[i + j * N_r].vertices[2] = i + j * 2 * N_r;
6491 cells[i + j * N_r].vertices[3] = (i + 1) % N_r + j * 2 * N_r;
6493 cells[i + j * N_r].vertices[4] = N_r + i + (j + 1) * 2 * N_r;
6494 cells[i + j * N_r].vertices[5] =
6495 N_r + ((i + 1) % N_r) + (j + 1) * 2 * N_r;
6496 cells[i + j * N_r].vertices[6] = N_r + i + j * 2 * N_r;
6497 cells[i + j * N_r].vertices[7] = N_r + ((i + 1) % N_r) + j * 2 * N_r;
6499 cells[i + j * N_r].material_id = 0;
6518 double eps_z = 1
e-6 * length;
6523 double face_inner_radius = std::numeric_limits<double>::max();
6524 double face_outer_radius = 0.;
6529 for (
const auto &cell : tria.active_cell_iterators())
6532 if (!cell->face(f)->at_boundary())
6535 const auto face_center = cell->face(f)->center();
6536 const double z = face_center[2];
6538 if ((std::fabs(z) > eps_z) &&
6539 (std::fabs(z - length) > eps_z))
6541 const double radius =
std::sqrt(face_center[0] * face_center[0] +
6542 face_center[1] * face_center[1]);
6543 face_inner_radius =
std::min(face_inner_radius, radius);
6544 face_outer_radius =
std::max(face_outer_radius, radius);
6548 double mid_radial_distance = 0.5 * (face_outer_radius - face_inner_radius);
6550 for (
const auto &cell : tria.active_cell_iterators())
6553 if (cell->face(f)->at_boundary())
6555 const auto face_center = cell->face(f)->center();
6557 const double radius =
std::sqrt(face_center[0] * face_center[0] +
6558 face_center[1] * face_center[1]);
6560 const double z = face_center[2];
6562 if (std::fabs(z) < eps_z)
6564 cell->face(f)->set_boundary_id(2);
6566 else if (std::fabs(z - length) <
6569 cell->face(f)->set_boundary_id(3);
6571 else if (std::fabs(radius - face_inner_radius) >
6572 mid_radial_distance)
6574 cell->face(f)->set_boundary_id(1);
6576 else if (std::fabs(radius - face_inner_radius) <
6577 mid_radial_distance)
6579 cell->face(f)->set_boundary_id(0);
6588 template <
int dim,
int spacedim>
6593 const double duplicated_vertex_tolerance,
6594 const bool copy_manifold_ids,
6595 const bool copy_boundary_ids)
6597 std::vector<Point<spacedim>> vertices;
6598 std::vector<CellData<dim>> cells;
6601 unsigned int n_accumulated_vertices = 0;
6605 ExcMessage(
"The input triangulations must be non-empty "
6606 "and must not be refined."));
6608 auto [tria_vertices, tria_cells, tria_subcell_data] =
6618 vertices.insert(vertices.end(),
6619 tria_vertices.begin(),
6620 tria_vertices.end());
6623 for (
unsigned int &vertex_n : cell_data.vertices)
6624 vertex_n += n_accumulated_vertices;
6625 cells.push_back(cell_data);
6629 if (copy_manifold_ids)
6631 for (
CellData<1> &line_data : tria_subcell_data.boundary_lines)
6635 for (
unsigned int &vertex_n : line_data.vertices)
6636 vertex_n += n_accumulated_vertices;
6637 line_data.boundary_id =
6642 for (
CellData<2> &quad_data : tria_subcell_data.boundary_quads)
6646 for (
unsigned int &vertex_n : quad_data.vertices)
6647 vertex_n += n_accumulated_vertices;
6648 quad_data.boundary_id =
6658 std::vector<unsigned int> considered_vertices;
6662 considered_vertices,
6663 duplicated_vertex_tolerance);
6667 if (std::all_of(cells.begin(), cells.end(), [](
const auto &cell) {
6668 return cell.vertices.size() ==
6669 ReferenceCells::get_hypercube<dim>().n_vertices();
6675 if (!copy_manifold_ids)
6678 if (copy_boundary_ids)
6680 auto result_cell = result.
begin();
6681 for (
const auto &tria : triangulations)
6683 for (
const auto &cell : tria->cell_iterators())
6685 for (
const auto &f : cell->face_indices())
6686 if (result_cell->face(f)->at_boundary())
6687 result_cell->face(f)->set_boundary_id(
6694 Assert(duplicated_vertex_tolerance > 0.0 ||
6695 n_accumulated_vertices == result.
n_vertices(),
6701 template <
int dim,
int spacedim>
6706 const double duplicated_vertex_tolerance,
6707 const bool copy_manifold_ids,
6708 const bool copy_boundary_ids)
6711 if (triangulation_1.
n_cells() == 0)
6713 if (&result != &triangulation_2)
6716 else if (triangulation_2.
n_cells() == 0)
6718 if (&result != &triangulation_1)
6724 duplicated_vertex_tolerance,
6754 template <
int structdim>
6758 static_assert(structdim == 1 || structdim == 2,
6759 "This function is only implemented for lines and "
6767 std::sort(std::begin(cell_data.vertices),
6768 std::end(cell_data.vertices));
6769 else if (structdim == 2)
6772 std::array<unsigned int, 4> renumbering{};
6773 std::copy(std::begin(cell_data.vertices),
6774 std::end(cell_data.vertices),
6775 renumbering.begin());
6786 std::swap(renumbering[2], renumbering[3]);
6787 std::rotate(renumbering.begin(),
6788 std::min_element(renumbering.begin(),
6792 std::swap(renumbering[2], renumbering[3]);
6800 if (renumbering[1] > renumbering[2])
6801 std::swap(renumbering[1], renumbering[2]);
6802 std::copy(renumbering.begin(),
6804 std::begin(cell_data.vertices));
6811 return std::lexicographical_compare(std::begin(a.
vertices),
6813 std::begin(
b.vertices),
6814 std::end(
b.vertices));
6816 std::sort(subcell_data.begin(), subcell_data.end(), compare);
6821 auto left = subcell_data.begin();
6822 while (left != subcell_data.end())
6825 std::upper_bound(left, subcell_data.end(), *left, compare);
6828 if (left + 1 != right)
6829 for (
auto it = left; it != right; ++it)
6832 Assert(it->manifold_id == left->manifold_id,
6834 "In the process of grid generation a single "
6835 "line or quadrilateral has been assigned two "
6836 "different manifold ids. This can happen when "
6837 "a Triangulation is copied, e.g., via "
6838 "GridGenerator::replicate_triangulation() and "
6839 "not all external boundary faces have the same "
6840 "manifold id. Double check that all faces "
6841 "which you expect to be merged together have "
6842 "the same manifold id."));
6847 subcell_data.erase(std::unique(subcell_data.begin(), subcell_data.end()),
6848 subcell_data.end());
6854 template <
int dim,
int spacedim>
6857 const std::vector<unsigned int> &extents,
6863 for (
const auto &extent : extents)
6866 "The Triangulation must be copied at least one time in "
6867 "each coordinate dimension."));
6870 const auto &
min = bbox.get_boundary_points().first;
6871 const auto &
max = bbox.get_boundary_points().second;
6873 std::array<Tensor<1, spacedim>, dim> offsets;
6874 for (
unsigned int d = 0;
d < dim; ++
d)
6875 offsets[d][d] = max[d] - min[d];
6879 for (
unsigned int d = 0;
d < dim; ++
d)
6881 auto [input_vertices, input_cell_data, input_subcell_data] =
6884 std::vector<Point<spacedim>> output_vertices = input_vertices;
6885 std::vector<CellData<dim>> output_cell_data = input_cell_data;
6886 SubCellData output_subcell_data = input_subcell_data;
6888 for (
unsigned int k = 1; k < extents[
d]; ++k)
6890 const std::size_t vertex_offset = k * input_vertices.size();
6893 output_vertices.push_back(
point + double(k) * offsets[
d]);
6897 output_cell_data.push_back(cell_data);
6898 for (
unsigned int &vertex : output_cell_data.back().vertices)
6899 vertex += vertex_offset;
6903 input_subcell_data.boundary_lines)
6906 for (
unsigned int &vertex :
6907 output_subcell_data.boundary_lines.back().vertices)
6908 vertex += vertex_offset;
6911 input_subcell_data.boundary_quads)
6914 for (
unsigned int &vertex :
6915 output_subcell_data.boundary_quads.back().vertices)
6916 vertex += vertex_offset;
6921 std::vector<unsigned int> boundary_vertices;
6925 output_subcell_data,
6940 output_subcell_data.boundary_lines)
6946 tria_to_replicate.
clear();
6949 output_subcell_data);
6957 template <
int dim,
int spacedim>
6965 ExcMessage(
"The two input triangulations are not derived from "
6966 "the same coarse mesh as required."));
6969 &triangulation_1) ==
nullptr) &&
6972 &triangulation_2) ==
nullptr),
6973 ExcMessage(
"The source triangulations for this function must both "
6974 "be available entirely locally, and not be distributed "
6975 "triangulations."));
6992 for (
unsigned int iteration = 0; iteration < triangulation_2.
n_levels();
6998 bool any_cell_flagged =
false;
6999 for (
const auto &result_cell : result.active_cell_iterators())
7000 if (intergrid_map[result_cell]->has_children())
7002 any_cell_flagged =
true;
7003 result_cell->set_refine_flag();
7006 if (any_cell_flagged ==
false)
7015 template <
int dim,
int spacedim>
7025 std::vector<Point<spacedim>> vertices = input_triangulation.
get_vertices();
7029 std::vector<CellData<dim>> cells;
7030 for (
const auto &cell : input_triangulation.active_cell_iterators())
7031 if (cells_to_remove.find(cell) == cells_to_remove.
end())
7033 Assert(
static_cast<unsigned int>(cell->level()) ==
7034 input_triangulation.
n_levels() - 1,
7036 "Your input triangulation appears to have "
7037 "adaptively refined cells. This is not allowed. You can "
7038 "only call this function on a triangulation in which "
7039 "all cells are on the same refinement level."));
7043 this_cell.vertices[v] = cell->vertex_index(v);
7045 cells.push_back(this_cell);
7051 std::vector<unsigned int> considered_vertices;
7055 considered_vertices);
7067 const unsigned int n_slices,
7068 const double height,
7070 const bool copy_manifold_ids,
7071 const std::vector<types::manifold_id> &manifold_priorities)
7075 "The input triangulation must be a coarse mesh, i.e., it must "
7076 "not have been refined."));
7078 ExcMessage(
"The output triangulation object needs to be empty."));
7080 ExcMessage(
"The given height for extrusion must be positive."));
7083 "The number of slices for extrusion must be at least 2."));
7085 const double delta_h = height / (n_slices - 1);
7086 std::vector<double> slices_z_values;
7087 slices_z_values.reserve(n_slices);
7088 for (
unsigned int i = 0; i < n_slices; ++i)
7089 slices_z_values.push_back(i * delta_h);
7091 input, slices_z_values, result, copy_manifold_ids, manifold_priorities);
7099 const unsigned int n_slices,
7100 const double height,
7102 const bool copy_manifold_ids,
7103 const std::vector<types::manifold_id> &manifold_priorities)
7109 (void)copy_manifold_ids;
7110 (void)manifold_priorities;
7114 "GridTools::extrude_triangulation() is only available "
7115 "for Triangulation<3, 3> as output triangulation."));
7123 const std::vector<double> &slice_coordinates,
7125 const bool copy_manifold_ids,
7126 const std::vector<types::manifold_id> &manifold_priorities)
7130 "The input triangulation must be a coarse mesh, i.e., it must "
7131 "not have been refined."));
7133 ExcMessage(
"The output triangulation object needs to be empty."));
7134 Assert(slice_coordinates.size() >= 2,
7136 "The number of slices for extrusion must be at least 2."));
7137 Assert(std::is_sorted(slice_coordinates.begin(), slice_coordinates.end()),
7138 ExcMessage(
"Slice z-coordinates should be in ascending order"));
7141 "This function is only implemented for quadrilateral meshes."));
7143 const auto priorities = [&]() -> std::vector<types::manifold_id> {
7147 if (0 < manifold_priorities.size())
7152 std::vector<types::manifold_id> sorted_manifold_priorities =
7153 manifold_priorities;
7154 std::sort(sorted_manifold_priorities.begin(),
7155 sorted_manifold_priorities.end());
7156 Assert(std::unique(sorted_manifold_priorities.begin(),
7157 sorted_manifold_priorities.end()) ==
7158 sorted_manifold_priorities.end(),
7160 "The given vector of manifold ids may not contain any "
7161 "duplicated entries."));
7162 std::vector<types::manifold_id> sorted_manifold_ids =
7164 std::sort(sorted_manifold_ids.begin(), sorted_manifold_ids.end());
7165 if (sorted_manifold_priorities != sorted_manifold_ids)
7167 std::ostringstream message;
7168 message <<
"The given triangulation has manifold ids {";
7170 sorted_manifold_ids)
7173 message << sorted_manifold_ids.back() <<
"}, but \n"
7174 <<
" the given vector of manifold ids is {";
7176 manifold_priorities)
7180 << manifold_priorities.back() <<
"}.\n"
7181 <<
" These vectors should contain the same elements.\n";
7182 const std::string m = message.str();
7186 return manifold_priorities;
7190 std::vector<types::manifold_id> default_priorities =
7192 const auto first_tfi_it = std::partition(
7193 default_priorities.begin(),
7194 default_priorities.end(),
7196 return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
7197 &input.get_manifold(id)) == nullptr;
7199 std::sort(default_priorities.begin(), first_tfi_it);
7200 std::sort(first_tfi_it, default_priorities.end());
7202 return default_priorities;
7205 const std::size_t n_slices = slice_coordinates.size();
7206 std::vector<Point<3>> points(n_slices * input.
n_vertices());
7207 std::vector<CellData<3>> cells;
7212 for (std::size_t slice_n = 0; slice_n < n_slices; ++slice_n)
7214 for (std::size_t vertex_n = 0; vertex_n < input.
n_vertices();
7218 points[slice_n * input.
n_vertices() + vertex_n] =
7219 Point<3>(vertex[0], vertex[1], slice_coordinates[slice_n]);
7225 for (
const auto &cell : input.active_cell_iterators())
7227 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
7230 for (
const unsigned int vertex_n :
7234 cell->vertex_index(vertex_n) + slice_n * input.
n_vertices();
7237 cell->vertex_index(vertex_n) +
7242 if (copy_manifold_ids)
7244 cells.push_back(this_cell);
7255 for (
const auto &face : input.active_face_iterators())
7259 if (face->at_boundary())
7261 if (copy_manifold_ids)
7263 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
7266 face->vertex_index(0) + slice_n * input.
n_vertices();
7268 face->vertex_index(1) + slice_n * input.
n_vertices();
7270 face->vertex_index(0) + (slice_n + 1) * input.
n_vertices();
7272 face->vertex_index(1) + (slice_n + 1) * input.
n_vertices();
7273 quads.push_back(quad);
7279 if (copy_manifold_ids)
7280 for (
const auto &cell : input.active_cell_iterators())
7285 for (std::size_t slice_n = 1; slice_n < n_slices - 1; ++slice_n)
7288 cell->vertex_index(0) + slice_n * input.
n_vertices();
7290 cell->vertex_index(1) + slice_n * input.
n_vertices();
7292 cell->vertex_index(2) + slice_n * input.
n_vertices();
7294 cell->vertex_index(3) + slice_n * input.
n_vertices();
7295 quads.push_back(quad);
7306 "The input triangulation to this function is using boundary "
7307 "indicators in a range that do not allow using "
7308 "max_boundary_id+1 and max_boundary_id+2 as boundary "
7309 "indicators for the bottom and top faces of the "
7310 "extruded triangulation."));
7313 for (
const auto &cell : input.active_cell_iterators())
7317 quad.
vertices[0] = cell->vertex_index(0);
7318 quad.
vertices[1] = cell->vertex_index(1);
7319 quad.
vertices[2] = cell->vertex_index(2);
7320 quad.
vertices[3] = cell->vertex_index(3);
7321 if (copy_manifold_ids)
7323 quads.push_back(quad);
7326 for (
unsigned int &vertex : quad.vertices)
7327 vertex += (n_slices - 1) * input.n_vertices();
7328 if (copy_manifold_ids)
7330 quads.push_back(quad);
7343 for (
auto manifold_id_it = priorities.rbegin();
7344 manifold_id_it != priorities.rend();
7346 for (
const auto &face : result.active_face_iterators())
7348 for (unsigned
int line_n = 0;
7349 line_n < GeometryInfo<3>::lines_per_face;
7351 face->line(line_n)->set_manifold_id(*manifold_id_it);
7359 const std::vector<double> &slice_coordinates,
7361 const bool copy_manifold_ids,
7362 const std::vector<types::manifold_id> &manifold_priorities)
7365 (void)slice_coordinates;
7367 (void)copy_manifold_ids;
7368 (void)manifold_priorities;
7372 "GridTools::extrude_triangulation() is only available "
7373 "for Triangulation<3, 3> as output triangulation."));
7392 template <
int spacedim>
7395 const double inner_radius,
7396 const double outer_radius,
7403 Assert(inner_radius < outer_radius,
7404 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7413 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
7414 for (
const auto &cell :
triangulation.active_cell_iterators())
7417 if (cell->face(f)->at_boundary())
7420 (
std::
fabs(cell->face(f)->vertex(v).
norm() - outer_radius) <
7421 1
e-12 * outer_radius)
7426 std::
fabs(cell->face(f)->vertex(v)[1])) <
7427 1
e-12 * outer_radius))
7428 cell->face(f)->vertex(v) *=
std::
sqrt(2.);
7430 const double eps = 1
e-3 * outer_radius;
7431 for (
const auto &cell :
triangulation.active_cell_iterators())
7433 for (
const unsigned int f : cell->face_indices())
7434 if (cell->face(f)->at_boundary())
7436 const double dx = cell->face(f)->center()[0] - center[0];
7437 const double dy = cell->face(f)->center()[1] - center[1];
7440 if (
std::abs(dx + outer_radius) < eps)
7441 cell->face(f)->set_boundary_id(0);
7442 else if (
std::abs(dx - outer_radius) < eps)
7443 cell->face(f)->set_boundary_id(1);
7444 else if (
std::abs(dy + outer_radius) < eps)
7445 cell->face(f)->set_boundary_id(2);
7446 else if (
std::abs(dy - outer_radius) < eps)
7447 cell->face(f)->set_boundary_id(3);
7450 cell->face(f)->set_boundary_id(4);
7451 cell->face(f)->set_manifold_id(0);
7456 const double d = (cell->face(f)->center() - center).
norm();
7457 if (d - inner_radius < 0)
7459 cell->face(f)->set_boundary_id(1);
7460 cell->face(f)->set_manifold_id(0);
7463 cell->face(f)->set_boundary_id(0);
7475 const double inner_radius,
7476 const double outer_radius,
7478 const unsigned int width_repetition,
7494 const double inner_radius,
7495 const double outer_radius,
7497 const unsigned int width_repetition,
7514 const double inner_radius,
7515 const double outer_radius,
7516 const unsigned int n_shells,
7517 const double skewness,
7518 const unsigned int n_cells,
7524 Assert(inner_radius < outer_radius,
7525 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7529 std::vector<double> radii;
7530 radii.push_back(inner_radius);
7531 for (
unsigned int shell_n = 1; shell_n < n_shells; ++shell_n)
7532 if (skewness == 0.0)
7534 radii.push_back(inner_radius +
7535 (outer_radius - inner_radius) *
7536 (1.0 - (1.0 -
double(shell_n) / n_shells)));
7540 (outer_radius - inner_radius) *
7541 (1.0 -
std::tanh(skewness * (1.0 -
double(shell_n) / n_shells)) /
7543 radii.push_back(outer_radius);
7545 double grid_vertex_tolerance = 0.0;
7546 for (
unsigned int shell_n = 0; shell_n < radii.size() - 1; ++shell_n)
7553 n_cells == 0 ? (dim == 2 ? 8 : 12) :
7558 if (grid_vertex_tolerance == 0.0)
7559 grid_vertex_tolerance =
7560 0.5 * internal::minimal_vertex_distance(current_shell);
7567 grid_vertex_tolerance);
7579 constexpr double radial_vertex_tolerance =
7580 100.0 * std::numeric_limits<double>::epsilon();
7581 auto assert_vertex_distance_within_tolerance =
7582 [center, radial_vertex_tolerance](
7584 const double radius) {
7586 (void)radial_vertex_tolerance;
7589 for (
unsigned int vertex_n = 0;
7590 vertex_n < GeometryInfo<dim>::vertices_per_face;
7593 Assert(
std::abs((face->vertex(vertex_n) - center).norm() - radius) <
7594 (center.
norm() + radius) * radial_vertex_tolerance,
7599 for (
const auto &cell :
triangulation.active_cell_iterators())
7602 auto face = cell->face(face_n);
7603 if (face->at_boundary())
7605 if (((face->vertex(0) - center).norm() - inner_radius) <
7606 (center.
norm() + inner_radius) * radial_vertex_tolerance)
7609 assert_vertex_distance_within_tolerance(face, inner_radius);
7610 face->set_all_boundary_ids(0);
7615 assert_vertex_distance_within_tolerance(face, outer_radius);
7616 face->set_all_boundary_ids(1);
7627 const double inner_radius,
7628 const double outer_radius,
7630 const unsigned int Nz,
7635 Assert(inner_radius < outer_radius,
7636 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7650 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
7651 for (
const auto &cell :
triangulation.active_cell_iterators())
7653 for (
const auto f : cell->face_indices())
7654 if (cell->face(f)->at_boundary())
7658 const unsigned int vv = cell->face(f)->vertex_index(v);
7659 if (treated_vertices[vv] ==
false)
7661 treated_vertices[vv] =
true;
7671 cell->face(f)->vertex(v);
7672 if ((std::fabs(std::fabs(vertex_location[0]) -
7673 std::fabs(vertex_location[1])) <
7674 1e-12 * outer_radius) &&
7675 (std::fabs(vertex_location[0] * vertex_location[0] +
7676 vertex_location[1] * vertex_location[1] -
7677 outer_radius * outer_radius) <
7678 1e-12 * outer_radius))
7679 cell->face(f)->vertex(v) =
7682 vertex_location[2]);
7687 double eps = 1
e-3 * outer_radius;
7688 for (
const auto &cell :
triangulation.active_cell_iterators())
7690 for (
const unsigned int f : cell->face_indices())
7691 if (cell->face(f)->at_boundary())
7693 const double dx = cell->face(f)->center()[0];
7694 const double dy = cell->face(f)->center()[1];
7695 const double dz = cell->face(f)->center()[2];
7699 if (
std::abs(dx + outer_radius) < eps)
7700 cell->face(f)->set_boundary_id(0);
7702 else if (
std::abs(dx - outer_radius) < eps)
7703 cell->face(f)->set_boundary_id(1);
7705 else if (
std::abs(dy + outer_radius) < eps)
7706 cell->face(f)->set_boundary_id(2);
7708 else if (
std::abs(dy - outer_radius) < eps)
7709 cell->face(f)->set_boundary_id(3);
7712 cell->face(f)->set_boundary_id(4);
7715 cell->face(f)->set_boundary_id(5);
7719 cell->face(f)->set_all_boundary_ids(6);
7720 cell->face(f)->set_all_manifold_ids(0);
7727 const double d = c.
norm();
7728 if (d - inner_radius < 0)
7730 cell->face(f)->set_all_boundary_ids(1);
7731 cell->face(f)->set_all_manifold_ids(0);
7734 cell->face(f)->set_boundary_id(0);
7743 template <
int dim,
int spacedim1,
int spacedim2>
7750 &in_tria) ==
nullptr),
7752 "This function cannot be used on "
7753 "parallel::distributed::Triangulation objects as inputs."));
7755 ExcMessage(
"This function does not work for meshes that have "
7759 const unsigned int spacedim =
std::min(spacedim1, spacedim2);
7760 const std::vector<Point<spacedim1>> &in_vertices = in_tria.
get_vertices();
7764 std::vector<Point<spacedim2>> v(in_vertices.size());
7765 for (
unsigned int i = 0; i < in_vertices.size(); ++i)
7766 for (
unsigned int d = 0;
d < spacedim; ++
d)
7767 v[i][d] = in_vertices[i][d];
7770 for (
const auto &cell : in_tria.active_cell_iterators())
7772 const unsigned int id = cell->active_cell_index();
7774 cells[id].vertices.resize(cell->n_vertices());
7776 cells[id].vertices[i] = cell->vertex_index(i);
7777 cells[id].material_id = cell->material_id();
7778 cells[id].manifold_id = cell->manifold_id();
7792 std::vector<bool> user_flags_line;
7795 .clear_user_flags_line();
7799 for (
const auto &face : in_tria.active_face_iterators())
7801 if (face->at_boundary())
7805 boundary_line.
vertices.resize(face->n_vertices());
7807 boundary_line.vertices[i] = face->vertex_index(i);
7812 std::move(boundary_line));
7825 if ((face->user_flag_set() ==
false) &&
7830 boundary_line.
vertices.resize(face->n_vertices());
7832 boundary_line.vertices[i] = face->vertex_index(i);
7838 std::move(boundary_line));
7840 face->set_user_flag();
7846 .load_user_flags_line(user_flags_line);
7853 std::vector<bool> user_flags_line;
7856 .clear_user_flags_line();
7858 std::vector<bool> user_flags_quad;
7861 .clear_user_flags_quad();
7865 for (
const auto &face : in_tria.active_face_iterators())
7867 if (face->at_boundary())
7871 boundary_face.
vertices.resize(face->n_vertices());
7873 boundary_face.vertices[i] = face->vertex_index(i);
7878 std::move(boundary_face));
7885 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7886 if (face->line(e)->user_flag_set() ==
false)
7889 spacedim1>::line_iterator
7890 edge = face->line(e);
7893 boundary_edge.
vertices.resize(edge->n_vertices());
7895 boundary_edge.vertices[i] = edge->vertex_index(i);
7900 std::move(boundary_edge));
7902 edge->set_user_flag();
7920 if (face->user_flag_set() ==
false)
7926 boundary_face.
vertices.resize(face->n_vertices());
7928 boundary_face.vertices[i] = face->vertex_index(i);
7934 std::move(boundary_face));
7936 face->set_user_flag();
7945 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7946 if (face->line(e)->at_boundary() ==
false)
7947 if (face->line(e)->user_flag_set() ==
false)
7950 line_iterator edge = face->line(e);
7955 boundary_edge.vertices[i] =
7956 edge->vertex_index(i);
7962 std::move(boundary_edge));
7964 edge->set_user_flag();
7971 .load_user_flags_line(user_flags_line);
7973 .load_user_flags_quad(user_flags_quad);
7983 for (
const auto i : out_tria.get_manifold_ids())
7985 out_tria.set_manifold(i,
FlatManifold<dim, spacedim2>());
7990 template <
int dim,
int spacedim>
7999 "GridGenerator::convert_hypercube_to_simplex_mesh() expects a mesh that consists only of quads/hexes."));
8031 {{{0, 1, 12, 10}}, {{2, 3, 11, 12}}, {{7, 6, 11, 13}},
8032 {{5, 4, 13, 10}}, {{0, 2, 8, 12}}, {{4, 6, 13, 8}},
8033 {{5, 13, 7, 9}}, {{1, 9, 3, 12}}, {{0, 8, 4, 10}},
8034 {{1, 5, 9, 10}}, {{3, 7, 11, 9}}, {{2, 6, 8, 11}},
8035 {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
8036 {{12, 13, 8, 10}}, {{13, 8, 10, 4}}, {{13, 10, 9, 5}},
8037 {{13, 9, 11, 7}}, {{13, 11, 8, 6}}, {{10, 12, 9, 1}},
8038 {{9, 12, 11, 3}}, {{11, 12, 8, 2}}, {{8, 12, 10, 0}}}};
8045 vertex_ids_for_boundary_faces_2d = {{{{{{0, 4}}, {{4, 2}}}},
8046 {{{{1, 5}}, {{5, 3}}}},
8047 {{{{0, 6}}, {{6, 1}}}},
8048 {{{{2, 7}}, {{7, 3}}}}}};
8055 vertex_ids_for_boundary_faces_3d = {
8056 {{{{{0, 4, 8}}, {{4, 8, 6}}, {{8, 6, 2}}, {{0, 2, 8}}}},
8057 {{{{1, 3, 9}}, {{3, 9, 7}}, {{9, 7, 5}}, {{1, 9, 5}}}},
8058 {{{{0, 1, 10}}, {{1, 10, 5}}, {{10, 5, 4}}, {{0, 10, 4}}}},
8059 {{{{2, 3, 11}}, {{3, 11, 7}}, {{11, 7, 6}}, {{2, 11, 6}}}},
8060 {{{{0, 1, 12}}, {{1, 12, 3}}, {{12, 3, 2}}, {{0, 12, 2}}}},
8061 {{{{4, 5, 13}}, {{5, 13, 7}}, {{13, 7, 6}}, {{4, 13, 6}}}}}};
8080 {{{0, 12, 10}}, {{12, 1, 10}}, {{12, 1, 9}}, {{12, 3, 9}},
8081 {{12, 2, 11}}, {{12, 3, 11}}, {{12, 0, 8}}, {{12, 2, 8}},
8082 {{9, 13, 5}}, {{13, 7, 9}}, {{11, 7, 13}}, {{11, 6, 13}},
8083 {{4, 8, 13}}, {{6, 8, 13}}, {{4, 13, 10}}, {{13, 5, 10}},
8084 {{10, 9, 5}}, {{10, 9, 1}}, {{11, 9, 7}}, {{11, 9, 3}},
8085 {{8, 11, 2}}, {{8, 11, 6}}, {{8, 10, 0}}, {{8, 10, 4}},
8086 {{12, 3, 9}}, {{12, 9, 11}}, {{12, 3, 11}}, {{3, 9, 11}},
8087 {{2, 12, 8}}, {{2, 12, 11}}, {{2, 11, 8}}, {{8, 12, 11}},
8088 {{0, 12, 10}}, {{0, 12, 8}}, {{0, 8, 10}}, {{8, 10, 12}},
8089 {{12, 1, 10}}, {{12, 1, 9}}, {{1, 10, 9}}, {{10, 9, 12}},
8090 {{10, 8, 4}}, {{10, 8, 13}}, {{4, 13, 8}}, {{4, 13, 10}},
8091 {{10, 9, 13}}, {{10, 9, 5}}, {{13, 5, 10}}, {{13, 5, 9}},
8092 {{13, 7, 9}}, {{13, 7, 11}}, {{9, 11, 13}}, {{9, 11, 7}},
8093 {{8, 11, 13}}, {{8, 11, 6}}, {{6, 13, 8}}, {{6, 13, 11}},
8094 {{12, 13, 10}}, {{12, 13, 8}}, {{8, 10, 13}}, {{8, 10, 12}},
8095 {{12, 13, 10}}, {{12, 13, 9}}, {{10, 9, 13}}, {{10, 9, 12}},
8096 {{12, 13, 9}}, {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
8097 {{12, 13, 11}}, {{12, 13, 8}}, {{8, 11, 13}}, {{8, 11, 12}}}};
8103 {{{12, 10}}, {{12, 9}}, {{12, 11}}, {{12, 8}}, {{9, 13}}, {{11, 13}},
8104 {{8, 13}}, {{10, 13}}, {{10, 9}}, {{9, 11}}, {{11, 8}}, {{8, 10}},
8105 {{12, 9}}, {{12, 11}}, {{11, 9}}, {{12, 8}}, {{12, 11}}, {{11, 8}},
8106 {{12, 8}}, {{12, 10}}, {{10, 8}}, {{12, 10}}, {{12, 9}}, {{9, 10}},
8107 {{13, 10}}, {{13, 8}}, {{8, 10}}, {{13, 10}}, {{13, 9}}, {{9, 10}},
8108 {{13, 11}}, {{13, 9}}, {{11, 9}}, {{13, 11}}, {{13, 8}}, {{11, 8}},
8109 {{12, 13}}, {{8, 10}}, {{8, 13}}, {{10, 13}}, {{8, 12}}, {{10, 12}},
8110 {{12, 13}}, {{10, 9}}, {{10, 13}}, {{9, 13}}, {{10, 12}}, {{9, 12}},
8111 {{12, 13}}, {{9, 11}}, {{9, 13}}, {{11, 13}}, {{9, 12}}, {{11, 12}},
8112 {{12, 13}}, {{11, 8}}, {{11, 13}}, {{8, 13}}, {{11, 12}}, {{8, 12}}}};
8124 vertex_ids_for_new_boundary_edges_3d = {
8125 {{{{{4, 8}}, {{6, 8}}, {{0, 8}}, {{2, 8}}}},
8126 {{{{5, 9}}, {{7, 9}}, {{1, 9}}, {{3, 9}}}},
8127 {{{{4, 10}}, {{5, 10}}, {{0, 10}}, {{1, 10}}}},
8128 {{{{6, 11}}, {{7, 11}}, {{2, 11}}, {{3, 11}}}},
8129 {{{{2, 12}}, {{3, 12}}, {{0, 12}}, {{1, 12}}}},
8130 {{{{6, 13}}, {{7, 13}}, {{4, 13}}, {{5, 13}}}}}};
8132 std::vector<Point<spacedim>> vertices;
8133 std::vector<CellData<dim>> cells;
8138 std::vector<unsigned int> old_to_new_vertex_indices(
8140 std::vector<unsigned int> face_to_new_vertex_indices(
8150 for (
const auto &cell : ref_tria.cell_iterators())
8154 std::array<unsigned int, dim == 2 ? 9 : 14> local_vertex_indices;
8159 const auto v_global = cell->vertex_index(v);
8161 if (old_to_new_vertex_indices[v_global] ==
8164 old_to_new_vertex_indices[v_global] = vertices.size();
8165 vertices.push_back(cell->vertex(v));
8169 local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
8173 for (
const auto f : cell->face_indices())
8175 const auto f_global = cell->face_index(f);
8177 if (face_to_new_vertex_indices[f_global] ==
8180 face_to_new_vertex_indices[f_global] = vertices.size();
8182 cell->face(f)->center(
true));
8186 local_vertex_indices.size());
8187 local_vertex_indices[cell->n_vertices() + f] =
8188 face_to_new_vertex_indices[f_global];
8195 local_vertex_indices.size());
8196 local_vertex_indices[cell->n_vertices() + cell->n_faces()] =
8198 vertices.push_back(cell->center(
true));
8202 const auto add_cell = [&](
const unsigned int struct_dim,
8203 const auto &index_vertices,
8204 const unsigned int material_or_boundary_id,
8208 if (struct_dim < dim &&
8213 if (struct_dim == dim)
8218 cell_data.material_id =
8219 material_or_boundary_id;
8221 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8224 local_vertex_indices.size());
8225 cell_data.vertices[i] =
8226 local_vertex_indices[index_vertices[i]];
8228 cells.push_back(cell_data);
8230 else if (dim == 2 && struct_dim == 1)
8234 boundary_line.
boundary_id = material_or_boundary_id;
8236 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8239 local_vertex_indices.size());
8241 local_vertex_indices[index_vertices[i]];
8245 else if (dim == 3 && struct_dim == 2)
8249 boundary_quad.material_id = material_or_boundary_id;
8251 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8254 local_vertex_indices.size());
8255 boundary_quad.vertices[i] =
8256 local_vertex_indices[index_vertices[i]];
8260 else if (dim == 3 && struct_dim == 1)
8264 boundary_line.
boundary_id = material_or_boundary_id;
8266 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8269 local_vertex_indices.size());
8271 local_vertex_indices[index_vertices[i]];
8281 const auto material_id_cell = cell->material_id();
8287 const auto manifold_id_cell = cell->manifold_id();
8289 for (
const auto &cell_vertices : table_2D_cell)
8290 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8293 for (
const auto &face_vertices : vertex_ids_for_inner_faces_2d)
8305 const auto manifold_id_cell = cell->manifold_id();
8307 for (
const auto &cell_vertices : vertex_ids_for_cells_3d)
8308 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8312 for (
const auto &face_vertices : vertex_ids_for_inner_faces_3d)
8320 for (
const auto &edge_vertices : vertex_ids_for_inner_edges_3d)
8330 for (
const auto f : cell->face_indices())
8332 const auto bid = cell->face(f)->boundary_id();
8333 const auto mid = cell->face(f)->manifold_id();
8338 for (
const auto &face_vertices :
8339 vertex_ids_for_boundary_faces_2d[f])
8340 add_cell(1, face_vertices, bid, mid);
8346 for (
const auto &face_vertices :
8347 vertex_ids_for_boundary_faces_3d[f])
8348 add_cell(2, face_vertices, bid, mid);
8351 for (
const auto &edge_vertices :
8352 vertex_ids_for_new_boundary_edges_3d[f])
8353 add_cell(1, edge_vertices, bid, mid);
8363 for (
const auto e : cell->line_indices())
8365 auto edge = cell->line(e);
8371 old_to_new_vertex_indices[edge->vertex_index(0)];
8373 old_to_new_vertex_indices[edge->vertex_index(1)];
8385 for (
const auto i : out_tria.get_manifold_ids())
8387 out_tria.set_manifold(i,
FlatManifold<dim, spacedim>());
8392 template <
int spacedim>
8403 template <
int dim,
int spacedim>
8419 {{{0, 1, 3}}, {{1, 2, 3}}, {{2, 0, 3}}}};
8424 vertex_ids_for_boundary_faces_2d = {
8425 {{{{{0, 1}}}}, {{{{1, 2}}}}, {{{{2, 0}}}}}};
8429 {{{0, 1, 2, 4}}, {{1, 0, 3, 4}}, {{0, 2, 3, 4}}, {{2, 1, 3, 4}}}};
8434 vertex_ids_for_boundary_faces_3d = {
8435 {{{{{0, 1, 2}}}}, {{{{1, 0, 3}}}}, {{{{0, 2, 3}}}}, {{{{2, 1, 3}}}}}};
8441 vertex_ids_for_boundary_lines_3d = {{{{{{0, 1}}}},
8448 std::vector<Point<spacedim>> vertices;
8449 std::vector<CellData<dim>> cells;
8454 std::vector<unsigned int> old_to_new_vertex_indices(
8462 for (
const auto &cell : ref_tria.cell_iterators())
8465 cell->reference_cell().is_simplex(),
8467 "Cell with invalid ReferenceCell encountered. GridGenerator::alfeld_split_of_simplex_mesh() "
8468 "only supports simplex meshes as input."));
8472 std::array<
unsigned int, (dim == 3) ? 5 : 4> local_vertex_indices;
8478 const auto v_global = cell->vertex_index(v);
8480 if (old_to_new_vertex_indices[v_global] ==
8483 old_to_new_vertex_indices[v_global] = vertices.size();
8484 vertices.push_back(cell->vertex(v));
8488 local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
8494 local_vertex_indices[local_vertex_indices.size() - 1] = vertices.size();
8495 vertices.push_back(barycenter /
static_cast<double>(dim + 1));
8498 const auto add_cell = [&](
const unsigned int struct_dim,
8499 const auto &index_vertices,
8500 const unsigned int material_or_boundary_id,
8504 if (struct_dim < dim &&
8509 if (struct_dim == dim)
8514 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8517 local_vertex_indices.size());
8518 cell_data.vertices[i] =
8519 local_vertex_indices[index_vertices[i]];
8521 cell_data.material_id =
8522 material_or_boundary_id;
8524 cells.push_back(cell_data);
8526 else if (dim == 2 && struct_dim == 1)
8530 boundary_line.
boundary_id = material_or_boundary_id;
8532 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8535 local_vertex_indices.size());
8537 local_vertex_indices[index_vertices[i]];
8541 else if (dim == 3 && struct_dim == 2)
8545 boundary_quad.material_id = material_or_boundary_id;
8547 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8550 local_vertex_indices.size());
8551 boundary_quad.vertices[i] =
8552 local_vertex_indices[index_vertices[i]];
8556 else if (dim == 3 && struct_dim == 1)
8560 boundary_line.
boundary_id = material_or_boundary_id;
8562 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8565 local_vertex_indices.size());
8567 local_vertex_indices[index_vertices[i]];
8577 const auto material_id_cell = cell->material_id();
8578 const auto manifold_id_cell = cell->manifold_id();
8583 for (
const auto &cell_vertices : table_2D_cell)
8584 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8588 for (
const auto &cell_vertices : table_3D_cell)
8589 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8595 for (
const auto f : cell->face_indices())
8597 const auto bid = cell->face(f)->boundary_id();
8598 const auto mid = cell->face(f)->manifold_id();
8603 for (
const auto &face_vertices :
8604 vertex_ids_for_boundary_faces_2d[f])
8605 add_cell(1, face_vertices, bid, mid);
8609 for (
const auto &face_vertices :
8610 vertex_ids_for_boundary_faces_3d[f])
8611 add_cell(2, face_vertices, bid, mid);
8620 for (
const auto l : cell->line_indices())
8622 const auto bid = cell->line(l)->boundary_id();
8623 const auto mid = cell->line(l)->manifold_id();
8625 for (
const auto &line_vertices :
8626 vertex_ids_for_boundary_lines_3d[
l])
8627 add_cell(1, line_vertices, bid, mid);
8635 for (
const auto i : out_tria.get_manifold_ids())
8637 out_tria.set_manifold(i,
FlatManifold<dim, spacedim>());
8642 template <
template <
int,
int>
class MeshType,
int dim,
int spacedim>
8646 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8647 typename MeshType<dim, spacedim>::face_iterator>
8649 typename ExtractBoundaryMesh<MeshType, dim, spacedim>::return_type
8652 MeshType<dim - 1, spacedim> &surface_mesh,
8653 const std::set<types::boundary_id> &boundary_ids)
8657 &volume_mesh.get_triangulation()) ==
nullptr),
8666 const unsigned int boundary_dim = dim - 1;
8673 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>
8674 temporary_mapping_level0;
8679 std::vector<bool> touched(volume_mesh.get_triangulation().n_vertices(),
8683 std::vector<CellData<boundary_dim>> cells;
8685 std::vector<Point<spacedim>> vertices;
8688 std::map<unsigned int, unsigned int> map_vert_index;
8701 for (
unsigned int i1 = 0; i1 < GeometryInfo<spacedim>::faces_per_cell; ++i1)
8703 for (
unsigned int i2 = 0; i2 <
GeometryInfo<dim - 1>::vertices_per_cell;
8705 swap_matrix[i1][i2] = i2;
8711 std::swap(swap_matrix[0][1], swap_matrix[0][2]);
8712 std::swap(swap_matrix[2][1], swap_matrix[2][2]);
8713 std::swap(swap_matrix[4][1], swap_matrix[4][2]);
8717 std::swap(swap_matrix[1][0], swap_matrix[1][1]);
8718 std::swap(swap_matrix[2][0], swap_matrix[2][1]);
8723 for (
typename MeshType<dim, spacedim>::cell_iterator cell =
8724 volume_mesh.begin(0);
8725 cell != volume_mesh.end(0);
8727 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
8729 const typename MeshType<dim, spacedim>::face_iterator face =
8732 if (face->at_boundary() &&
8733 (boundary_ids.empty() ||
8734 (boundary_ids.find(face->boundary_id()) != boundary_ids.end())))
8738 for (
const unsigned int j :
8741 const unsigned int v_index = face->vertex_index(j);
8743 if (!touched[v_index])
8745 vertices.push_back(face->vertex(j));
8746 map_vert_index[v_index] = vertices.size() - 1;
8747 touched[v_index] =
true;
8750 c_data.
vertices[swap_matrix[i][j]] = map_vert_index[v_index];
8764 for (
unsigned int e = 0;
e < 4; ++
e)
8770 bool edge_found =
false;
8771 for (
auto &boundary_line : subcell_data.boundary_lines)
8772 if (((boundary_line.vertices[0] ==
8773 map_vert_index[face->line(
e)->vertex_index(0)]) &&
8774 (boundary_line.vertices[1] ==
8775 map_vert_index[face->line(
e)->vertex_index(
8777 ((boundary_line.vertices[0] ==
8778 map_vert_index[face->line(
e)->vertex_index(1)]) &&
8779 (boundary_line.vertices[1] ==
8780 map_vert_index[face->line(
e)->vertex_index(0)])))
8787 if (edge_found ==
true)
8794 map_vert_index[face->line(e)->vertex_index(0)];
8796 map_vert_index[face->line(e)->vertex_index(1)];
8803 cells.push_back(c_data);
8804 temporary_mapping_level0.push_back(std::make_pair(face, i));
8811 surface_mesh.get_triangulation())
8817 for (
const auto &cell : surface_mesh.active_cell_iterators())
8818 for (unsigned
int vertex = 0; vertex < 2; ++vertex)
8819 if (cell->face(vertex)->at_boundary())
8820 cell->face(vertex)->set_boundary_id(0);
8828 std::vector<std::pair<
8829 const typename MeshType<dim - 1, spacedim>::cell_iterator,
8830 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>>
8831 temporary_map_boundary_cell_face;
8832 for (
const auto &cell : surface_mesh.active_cell_iterators())
8833 temporary_map_boundary_cell_face.push_back(
8834 std::make_pair(cell, temporary_mapping_level0.at(cell->
index())));
8846 unsigned int index_cells_deepest_level = 0;
8849 bool changed =
false;
8853 std::vector<unsigned int> cells_refined;
8856 for (
unsigned int cell_n = index_cells_deepest_level;
8857 cell_n < temporary_map_boundary_cell_face.size();
8862 if (temporary_map_boundary_cell_face[cell_n]
8863 .
second.first->has_children())
8867 Assert(temporary_map_boundary_cell_face[cell_n]
8868 .
second.first->refinement_case() ==
8871 temporary_map_boundary_cell_face[cell_n]
8872 .first->set_refine_flag();
8873 cells_refined.push_back(cell_n);
8884 surface_mesh.get_triangulation())
8885 .execute_coarsening_and_refinement();
8888 index_cells_deepest_level = temporary_map_boundary_cell_face.size();
8889 for (
const auto &refined_cell_n : cells_refined)
8891 const typename MeshType<dim - 1, spacedim>::cell_iterator
8893 temporary_map_boundary_cell_face[refined_cell_n].first;
8894 const typename MeshType<dim,
8895 spacedim>::face_iterator refined_face =
8896 temporary_map_boundary_cell_face[refined_cell_n].second.first;
8897 const unsigned int refined_face_number =
8898 temporary_map_boundary_cell_face[refined_cell_n]
8900 for (
unsigned int child_n = 0;
8901 child_n < refined_cell->n_children();
8906 temporary_map_boundary_cell_face.push_back(
8907 std::make_pair(refined_cell->child(
8908 swap_matrix[refined_face_number][child_n]),
8909 std::make_pair(refined_face->child(child_n),
8910 refined_face_number)));
8920 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8921 typename MeshType<dim, spacedim>::face_iterator>
8922 surface_to_volume_mapping;
8923 for (
unsigned int i = 0; i < temporary_map_boundary_cell_face.size(); ++i)
8924 surface_to_volume_mapping[temporary_map_boundary_cell_face[i].
first] =
8925 temporary_map_boundary_cell_face[i].
second.first;
8928 const auto attached_mids =
8929 surface_mesh.get_triangulation().get_manifold_ids();
8930 for (
const auto i : volume_mesh.get_triangulation().get_manifold_ids())
8932 std::find(attached_mids.
begin(), attached_mids.
end(), i) ==
8933 attached_mids.
end())
8935 surface_mesh.get_triangulation())
8938 return surface_to_volume_mapping;
8943 template <
int dim,
int spacedim>
8947 const std::vector<unsigned int> &repetitions,
8954 std::vector<Point<spacedim>> vertices;
8955 std::vector<CellData<dim>> cells;
8961 (p2[1] - p1[1]) / repetitions[1]);
8964 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
8965 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
8970 for (
unsigned int j = 0; j < repetitions[1]; ++j)
8971 for (
unsigned int i = 0; i < repetitions[0]; ++i)
8974 std::array<unsigned int, 4> quad{{
8975 (j + 0) * (repetitions[0] + 1) + i + 0,
8976 (j + 0) * (repetitions[0] + 1) + i + 1,
8977 (j + 1) * (repetitions[0] + 1) + i + 0,
8978 (j + 1) * (repetitions[0] + 1) + i + 1
8984 tri.
vertices = {quad[0], quad[1], quad[2]};
8985 cells.push_back(tri);
8991 tri.
vertices = {quad[3], quad[2], quad[1]};
8992 cells.push_back(tri);
9000 (p2[1] - p1[1]) / repetitions[1],
9001 (p2[2] - p1[2]) / repetitions[2]);
9004 for (
unsigned int k = 0; k <= repetitions[2]; ++k)
9005 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
9006 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
9009 p1[2] + dx[2] * k));
9012 for (
unsigned int k = 0; k < repetitions[2]; ++k)
9013 for (
unsigned int j = 0; j < repetitions[1]; ++j)
9014 for (
unsigned int i = 0; i < repetitions[0]; ++i)
9017 std::array<unsigned int, 8> quad{
9018 {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9019 (j + 0) * (repetitions[0] + 1) + i + 0,
9020 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9021 (j + 0) * (repetitions[0] + 1) + i + 1,
9022 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9023 (j + 1) * (repetitions[0] + 1) + i + 0,
9024 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9025 (j + 1) * (repetitions[0] + 1) + i + 1,
9026 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9027 (j + 0) * (repetitions[0] + 1) + i + 0,
9028 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9029 (j + 0) * (repetitions[0] + 1) + i + 1,
9030 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9031 (j + 1) * (repetitions[0] + 1) + i + 0,
9032 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9033 (j + 1) * (repetitions[0] + 1) + i + 1}};
9038 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9039 cell.
vertices = {{quad[0], quad[1], quad[2], quad[4]}};
9041 cell.
vertices = {{quad[0], quad[1], quad[3], quad[5]}};
9043 cells.push_back(cell);
9049 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9050 cell.
vertices = {{quad[2], quad[1], quad[3], quad[7]}};
9052 cell.
vertices = {{quad[0], quad[3], quad[2], quad[6]}};
9053 cells.push_back(cell);
9059 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9060 cell.
vertices = {{quad[1], quad[4], quad[5], quad[7]}};
9062 cell.
vertices = {{quad[0], quad[4], quad[5], quad[6]}};
9063 cells.push_back(cell);
9069 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9070 cell.
vertices = {{quad[2], quad[4], quad[7], quad[6]}};
9072 cell.
vertices = {{quad[3], quad[5], quad[7], quad[6]}};
9073 cells.push_back(cell);
9079 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9080 cell.
vertices = {{quad[1], quad[2], quad[4], quad[7]}};
9082 cell.
vertices = {{quad[0], quad[3], quad[6], quad[5]}};
9083 cells.push_back(cell);
9105 double epsilon = std::numeric_limits<double>::max();
9106 for (
unsigned int i = 0; i < dim; ++i)
9108 0.01 * (
std::abs(p2[i] - p1[i]) / repetitions[i]));
9111 "The distance between corner points must be positive."));
9115 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
9121 template <
int dim,
int spacedim>
9124 const unsigned int repetitions,
9132 tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2},
colorize);
9138 {{repetitions, repetitions, repetitions}},
9151# include "grid/grid_generator.inst"
void make_mapping(const MeshType &source_grid, const MeshType &destination_grid)
void add_parameter(const std::string &entry, ParameterType ¶meter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
const Point< spacedim > center
numbers::NumberTraits< Number >::real_type norm() const
void initialize(const Triangulation< dim, spacedim > &triangulation)
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator > > &)
virtual types::global_cell_index n_global_active_cells() const
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
unsigned int n_faces() const
bool all_reference_cells_are_hyper_cube() const
void save_user_flags_line(std::ostream &out) const
face_iterator end_face() const
cell_iterator begin(const unsigned int level=0) const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_active_lines() const
unsigned int n_levels() const
cell_iterator end() const
virtual bool has_hanging_nodes() const
vertex_iterator begin_vertex() const
vertex_iterator end_vertex() const
virtual void execute_coarsening_and_refinement()
virtual unsigned int n_global_levels() const
cell_iterator last() const
face_iterator begin_face() const
unsigned int n_cells() const
void save_user_flags_quad(std::ostream &out) const
unsigned int n_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcLowerRange(int arg1, int arg2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void set_all_manifold_ids_on_boundary(const types::manifold_id number)
const Manifold< dim, spacedim > & get_manifold(const types::manifold_id number) const
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
virtual std::vector< types::manifold_id > get_manifold_ids() const
void set_manifold(const types::manifold_id number, const Manifold< dim, spacedim > &manifold_object)
void reset_all_manifolds()
void set_all_manifold_ids(const types::manifold_id number)
void consistently_order_cells(std::vector< CellData< dim > > &cells)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
std::vector< index_type > data
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
Expression fabs(const Expression &x)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void subdivided_hyper_cube_with_simplices(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double p1=0.0, const double p2=1.0, const bool colorize=false)
void parallelepiped(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void hyper_cross(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &sizes, const bool colorize_cells=false)
A center cell with stacks of cell protruding from each surface.
void hyper_cube_with_cylindrical_hole(Triangulation< dim, spacedim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
void hyper_ball_balanced(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > ¢er=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
void enclosed_hyper_cube(Triangulation< dim > &tria, const double left=0., const double right=1., const double thickness=1., const bool colorize=false)
void replicate_triangulation(const Triangulation< dim, spacedim > &input, const std::vector< unsigned int > &extents, Triangulation< dim, spacedim > &result)
Replicate a given triangulation in multiple coordinate axes.
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void general_cell(Triangulation< dim, spacedim > &tria, const std::vector< Point< spacedim > > &vertices, const bool colorize=false)
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void hyper_shell(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)
void eccentric_hyper_shell(Triangulation< dim > &triangulation, const Point< dim > &inner_center, const Point< dim > &outer_center, const double inner_radius, const double outer_radius, const unsigned int n_cells)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void moebius(Triangulation< 3, 3 > &tria, const unsigned int n_cells, const unsigned int n_rotations, const double R, const double r)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void half_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void quarter_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
Rectangular domain with rectangular pattern of holes.
void create_union_triangulation(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result)
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void subdivided_hyper_rectangle_with_simplices(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void non_standard_orientation_mesh(Triangulation< 2 > &tria, const unsigned int n_rotate_middle_square)
return_type extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
void alfeld_split_of_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_parallelepiped(Triangulation< dim > &tria, const unsigned int n_subdivisions, const Point< dim >(&corners)[dim], const bool colorize=false)
void subdivided_cylinder(Triangulation< dim > &tria, const unsigned int x_subdivisions, const double radius=1., const double half_length=1.)
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
void subdivided_hyper_L(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &bottom_left, const Point< dim > &top_right, const std::vector< int > &n_cells_to_remove)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=Point< spacedim >(), const double radius=1.)
void concentric_hyper_shells(Triangulation< dim > &triangulation, const Point< dim > ¢er, const double inner_radius=0.125, const double outer_radius=0.25, const unsigned int n_shells=1, const double skewness=0.1, const unsigned int n_cells_per_shell=0, const bool colorize=false)
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
void truncated_cone(Triangulation< dim > &tria, const double radius_0=1.0, const double radius_1=0.5, const double half_length=1.0)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cylinder_shell(Triangulation< dim > &tria, const double length, const double inner_radius, const double outer_radius, const unsigned int n_radial_cells=0, const unsigned int n_axial_cells=0, const bool colorize=false)
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
constexpr T fixed_power(const T t)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
long double gamma(const unsigned int n)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy(const T *begin, const T *end, U *dest)
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id internal_face_boundary_id
constexpr types::boundary_id invalid_boundary_id
constexpr types::manifold_id flat_manifold_id
constexpr types::material_id invalid_material_id
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > tanh(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< unsigned int > vertices
types::manifold_id manifold_id
types::material_id material_id
types::boundary_id boundary_id
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
std::vector< CellData< 2 > > boundary_quads
std::vector< CellData< 1 > > boundary_lines
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)