49 : airfoil_type(
"NACA")
51 , joukowski_center(-0.1, 0.14)
57 , incline_factor(0.35)
60 , n_subdivision_x_0(3)
61 , n_subdivision_x_1(2)
62 , n_subdivision_x_2(5)
64 , airfoil_sampling_factor(2)
69 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
71 Assert(incline_factor < 1.0 && incline_factor >= 0.0,
72 ExcMessage(
"incline_factor has to be in [0,1)!"));
85 "Mesh height measured from airfoil nose to horizontal boundaries");
89 "Length measured from airfoil leading edge to vertical outlet boundary");
93 "Define obliqueness of the vertical mesh around the airfoil");
102 "Type of airfoil geometry, either NACA or Joukowski airfoil",
117 "Joukowski circle center coordinates");
120 "Joukowski airfoil length leading to trailing edge");
128 "Number of global refinements");
130 "NumberSubdivisionX0",
132 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
134 "NumberSubdivisionX1",
136 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
138 "NumberSubdivisionX2",
140 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
143 "Number of subdivisions normal to airfoil");
147 "Factor to obtain a finer mesh at the airfoil surface");
162 static const unsigned int id_block_1 = 1;
163 static const unsigned int id_block_2 = 2;
164 static const unsigned int id_block_3 = 3;
165 static const unsigned int id_block_4 = 4;
166 static const unsigned int id_block_5 = 5;
167 static const unsigned int id_block_6 = 6;
172 MeshGenerator(
const AdditionalData &
data)
173 : refinements(
data.refinements)
174 , n_subdivision_x_0(
data.n_subdivision_x_0)
175 , n_subdivision_x_1(
data.n_subdivision_x_1)
176 , n_subdivision_x_2(
data.n_subdivision_x_2)
177 , n_subdivision_y(
data.n_subdivision_y)
179 , length_b2(
data.length_b2)
180 , incline_factor(
data.incline_factor)
181 , bias_factor(
data.bias_factor)
183 , n_cells_x_0(
Utilities::
pow(2, refinements) * n_subdivision_x_0)
184 , n_cells_x_1(
Utilities::
pow(2, refinements) * n_subdivision_x_1)
185 , n_cells_x_2(
Utilities::
pow(2, refinements) * n_subdivision_x_2)
186 , n_cells_y(
Utilities::
pow(2, refinements) * n_subdivision_y)
187 , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
189 , airfoil_1D(set_airfoil_length(
191 data.airfoil_type ==
"Joukowski" ?
192 joukowski(
data.joukowski_center,
193 n_points_on_each_side,
194 data.airfoil_sampling_factor) :
195 (
data.airfoil_type ==
"NACA" ?
197 n_points_on_each_side,
198 data.airfoil_sampling_factor) :
201 std::vector<Point<2>>{
205 data.airfoil_length))
206 , end_b0_x_u(airfoil_1D[0][n_cells_x_0][0])
207 , end_b0_x_l(airfoil_1D[1][n_cells_x_0][0])
208 , nose_x(airfoil_1D[0].front()[0])
209 , tail_x(airfoil_1D[0].back()[0])
210 , tail_y(airfoil_1D[0].back()[1])
211 , center_mesh(0.5 *
std::abs(end_b0_x_u + end_b0_x_l))
212 , length_b1_x(tail_x - center_mesh)
214 (edge_length +
std::abs(nose_x - center_mesh))))
218 ,
A(nose_x - edge_length, 0)
227 ,
J(tail_x + length_b2, 0)
232 data.airfoil_type ==
"NACA",
245 make_coarse_grid(tria_grid);
247 set_boundary_ids(tria_grid);
249 if (periodic_faces !=
nullptr)
252 tria_grid, 5, 4, 1, *periodic_faces);
270 (void)periodic_faces;
277 const unsigned int refinements;
280 const unsigned int n_subdivision_x_0;
283 const unsigned int n_subdivision_x_1;
286 const unsigned int n_subdivision_x_2;
290 const unsigned int n_subdivision_y;
297 const double length_b2;
301 const double incline_factor;
304 const double bias_factor;
307 const double edge_length;
310 const unsigned int n_cells_x_0;
313 const unsigned int n_cells_x_1;
316 const unsigned int n_cells_x_2;
320 const unsigned int n_cells_y;
323 const unsigned int n_points_on_each_side;
327 const std::array<std::vector<Point<2>>, 2> airfoil_1D;
331 const double end_b0_x_u;
335 const double end_b0_x_l;
349 const double center_mesh;
352 const double length_b1_x;
379 const Point<2> A, B,
C, D,
E,
F, G, H, I,
J,
K,
L;
418 static std::array<std::vector<Point<2>>, 2>
419 joukowski(
const Point<2> ¢erpoint,
420 const unsigned int number_points,
421 const unsigned int factor)
423 std::array<std::vector<Point<2>>, 2> airfoil_1D;
424 const unsigned int total_points = 2 * number_points - 2;
425 const unsigned int n_airfoilpoints = factor * total_points;
427 const auto jouk_points =
428 joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
431 std::vector<Point<2>> upper_points;
432 std::vector<Point<2>> lower_points;
436 unsigned int nose_index = 0;
437 unsigned int tail_index = 0;
438 double nose_x_coordinate = 0;
439 double tail_x_coordinate = 0;
443 for (
unsigned int i = 0; i < jouk_points.size(); ++i)
445 if (jouk_points[i][0] < nose_x_coordinate)
447 nose_x_coordinate = jouk_points[i][0];
450 if (jouk_points[i][0] > tail_x_coordinate)
452 tail_x_coordinate = jouk_points[i][0];
458 for (
unsigned int i = tail_index; i < jouk_points.size(); ++i)
459 upper_points.emplace_back(jouk_points[i]);
460 for (
unsigned int i = 0; i <= nose_index; ++i)
461 upper_points.emplace_back(jouk_points[i]);
462 std::reverse(upper_points.begin(), upper_points.end());
465 lower_points.insert(lower_points.end(),
466 jouk_points.begin() + nose_index,
467 jouk_points.begin() + tail_index + 1);
470 airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
471 airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
474 auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
475 const double nose_x_pos = vector.front()[0];
476 for (
auto &i : vector)
480 move_nose_to_origin(airfoil_1D[1]);
481 move_nose_to_origin(airfoil_1D[0]);
510 static std::vector<Point<2>>
511 joukowski_circle(
const Point<2> ¢er,
512 const unsigned int number_points)
514 std::vector<Point<2>> circle_points;
522 (1 - center[0]) * (1 - center[0]));
524 center[1] * center[1] + (1 + center[0]) * (1 + center[0]));
530 "Error creating lower circle: Circle for Joukowski-transform does"
531 " not enclose point zeta = -1! Choose different center "
537 circle_points.reserve(number_points);
538 for (
unsigned int i = 0; i < number_points; ++i)
539 circle_points.emplace_back(center[0] -
radius *
std::cos(i * theta),
543 return circle_points;
554 static std::vector<Point<2>>
555 joukowski_transform(
const std::vector<
Point<2>> &circle_points)
557 std::vector<Point<2>> joukowski_points(circle_points.size());
560 for (
unsigned int i = 0; i < circle_points.size(); ++i)
562 const double chi = circle_points[i][0];
563 const double eta = circle_points[i][1];
564 const std::complex<double> zeta(chi, eta);
565 const std::complex<double> z = zeta + 1. / zeta;
567 joukowski_points[i] = {std::real(z), std::imag(z)};
569 return joukowski_points;
588 static std::array<std::vector<Point<2>>, 2>
589 naca(
const std::string &serialnumber,
590 const unsigned int number_points,
591 const unsigned int factor)
595 const unsigned int n_airfoilpoints = factor * number_points;
598 return {{make_points_equidistant(
599 naca_create_points(serialnumber, n_airfoilpoints,
true),
601 make_points_equidistant(
602 naca_create_points(serialnumber, n_airfoilpoints,
false),
617 static std::vector<Point<2>>
618 naca_create_points(
const std::string &serialnumber,
619 const unsigned int number_points,
622 Assert(serialnumber.size() == 4,
623 ExcMessage(
"This NACA-serial number is not implemented!"));
625 return naca_create_points_4_digits(serialnumber,
644 static std::vector<Point<2>>
645 naca_create_points_4_digits(
const std::string &serialnumber,
646 const unsigned int number_points,
650 const unsigned int digit_0 = (serialnumber[0] -
'0');
651 const unsigned int digit_1 = (serialnumber[1] -
'0');
652 const unsigned int digit_2 = (serialnumber[2] -
'0');
653 const unsigned int digit_3 = (serialnumber[3] -
'0');
655 const unsigned int digit_23 = 10 * digit_2 + digit_3;
658 const double t =
static_cast<double>(digit_23) / 100.0;
660 std::vector<Point<2>> naca_points;
662 if (digit_0 == 0 && digit_1 == 0)
663 for (
unsigned int i = 0; i < number_points; ++i)
665 const double x = i * 1 / (1.0 * number_points - 1);
669 0.3516 * Utilities::fixed_power<2>(x) +
670 0.2843 * Utilities::fixed_power<3>(x) -
671 0.1036 * Utilities::fixed_power<4>(
675 naca_points.emplace_back(x, +y_t);
677 naca_points.emplace_back(x, -y_t);
680 for (
unsigned int i = 0; i < number_points; ++i)
682 const double m = 1.0 * digit_0 / 100;
683 const double p = 1.0 * digit_1 / 10;
684 const double x = i * 1 / (1.0 * number_points - 1);
688 m / Utilities::fixed_power<2>(p) *
689 (2 * p * x - Utilities::fixed_power<2>(x)) :
694 (x <= p) ? 2 * m / Utilities::fixed_power<2>(p) * (p - x) :
700 0.3516 * Utilities::fixed_power<2>(x) +
701 0.2843 * Utilities::fixed_power<3>(x) -
702 0.1036 * Utilities::fixed_power<4>(
708 naca_points.emplace_back(x - y_t *
std::sin(theta),
711 naca_points.emplace_back(x + y_t *
std::sin(theta),
728 static std::array<std::vector<Point<2>>, 2>
729 set_airfoil_length(
const std::array<std::vector<
Point<2>>, 2> &input,
730 const double desired_len)
732 std::array<std::vector<Point<2>>, 2> output;
733 output[0] = set_airfoil_length(input[0], desired_len);
734 output[1] = set_airfoil_length(input[1], desired_len);
746 static std::vector<Point<2>>
747 set_airfoil_length(
const std::vector<
Point<2>> &input,
748 const double desired_len)
750 std::vector<Point<2>> output = input;
753 desired_len / input.front().distance(input.back());
755 for (
auto &x : output)
771 static std::vector<Point<2>>
772 make_points_equidistant(
773 const std::vector<
Point<2>> &non_equidistant_points,
774 const unsigned int number_points)
776 const unsigned int n_points =
777 non_equidistant_points
781 std::vector<double> arclength_L(non_equidistant_points.size(), 0);
782 for (
unsigned int i = 0; i < non_equidistant_points.size() - 1; ++i)
785 non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
788 const auto airfoil_length =
790 const auto deltaX = airfoil_length / (number_points - 1);
794 std::vector<Point<2>> equidist(
797 equidist[0] = non_equidistant_points[0];
798 equidist[number_points - 1] = non_equidistant_points[n_points - 1];
802 for (
unsigned int j = 0, i = 1; j < n_points - 1; ++j)
805 const auto Lj = arclength_L[j];
806 const auto Ljp = arclength_L[j + 1];
808 while (Lj <= i * deltaX && i * deltaX <= Ljp &&
809 i < number_points - 1)
811 equidist[i] =
Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
812 (non_equidistant_points[j + 1] -
813 non_equidistant_points[j]) +
814 non_equidistant_points[j]);
834 std::vector<Triangulation<2>> trias(10);
838 const std::vector<Point<2>> &corner_vertices,
839 const std::vector<unsigned int> &repetitions,
851 auto &
point = it->vertex();
852 const double xi =
point[0];
853 const double eta =
point[1];
856 point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
857 (1 + xi) * (1 - eta) * corner_vertices[1] +
858 (1 - xi) * (1 + eta) * corner_vertices[2] +
859 (1 + xi) * (1 + eta) * corner_vertices[3]);
863 for (
auto cell : tria.active_cell_iterators())
871 {n_subdivision_y, n_subdivision_x_0},
875 {n_subdivision_y, n_subdivision_x_0},
879 {n_subdivision_x_1, n_subdivision_y},
883 {n_subdivision_x_1, n_subdivision_y},
887 {n_subdivision_x_2, n_subdivision_y},
891 {n_subdivision_x_2, n_subdivision_y},
918 for (
auto cell : tria.active_cell_iterators())
921 if (cell->face(f)->at_boundary() ==
false)
924 const auto mid = cell->material_id();
926 if ((mid == id_block_1 && f == 0) ||
927 (mid == id_block_4 && f == 0))
928 cell->face(f)->set_boundary_id(0);
929 else if ((mid == id_block_3 && f == 0) ||
930 (mid == id_block_6 && f == 2))
931 cell->face(f)->set_boundary_id(1);
932 else if ((mid == id_block_1 && f == 1) ||
933 (mid == id_block_2 && f == 1))
934 cell->face(f)->set_boundary_id(2);
935 else if ((mid == id_block_4 && f == 1) ||
936 (mid == id_block_5 && f == 3))
937 cell->face(f)->set_boundary_id(3);
938 else if ((mid == id_block_2 && f == 0) ||
939 (mid == id_block_3 && f == 2))
940 cell->face(f)->set_boundary_id(4);
941 else if ((mid == id_block_5 && f == 2) ||
942 (mid == id_block_6 && f == 0))
943 cell->face(f)->set_boundary_id(5);
963 std::vector<bool> vertex_processed(tria.
n_vertices(),
false);
980 for (
const auto &cell : tria.cell_iterators())
984 if (vertex_processed[cell->vertex_index(v)])
988 vertex_processed[cell->vertex_index(v)] =
true;
990 auto &node = cell->vertex(v);
993 if (cell->material_id() == id_block_1 ||
994 cell->material_id() == id_block_4)
1003 if (cell->material_id() == id_block_1)
1006 (node - horizontal_offset) +
1012 else if (cell->material_id() == id_block_4)
1015 (node - horizontal_offset) -
1021 const double trapeze_height =
1024 const double l_a =
std::cos(gamma) * edge_length;
1025 const double l_b = trapeze_height *
std::tan(gamma);
1027 const double x2 =
L - l_a - l_b;
1029 const double Dx = x1 + x2 + x3;
1030 const double deltax =
1032 const double dx = Dx / n_cells_x_0;
1033 const double dy = trapeze_height / n_cells_y;
1035 static_cast<int>(std::round((node_[0] - deltax) / dx));
1037 static_cast<int>(std::round(
std::abs(node_[1]) / dy));
1039 node_[0] =
numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1040 node_[1] =
height * (1.0 * iy) / n_cells_y;
1047 const double dy =
height / n_cells_y;
1049 static_cast<int>(std::round(node_[0] / dx));
1051 static_cast<int>(std::round(node_[1] / dy));
1052 const double alpha =
1053 bias_alpha(1 - (1.0 * iy) / n_cells_y);
1054 const double theta = node_[0];
1056 ((cell->material_id() == id_block_1) ?
1061 (cell->material_id() == id_block_1) ? (0) : (1))]
1067 else if (cell->material_id() == id_block_2 ||
1068 cell->material_id() == id_block_5)
1076 "Points D,C,G and E,F,I are not defined symmetric to "
1077 "x-axis, which is required to interpolate block 2"
1078 " and 5 with same geometric computations."));
1079 const double l_y = D[1] -
C[1];
1080 const double l_h = D[1] - l_y;
1081 const double by = -l_h / length_b1_x * (node[0] - H[0]);
1082 const double dy = (
height - by) / n_cells_y;
1083 const int iy =
static_cast<int>(
1084 std::round((
std::abs(node[1]) - by) / dy));
1085 const double dx = length_b1_x / n_cells_x_1;
1086 const int ix =
static_cast<int>(
1087 std::round(
std::abs(node[0] - center_mesh) / dx));
1089 const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1093 const Point<2> p(ix * dx + center_mesh +
1094 incline_factor * length_b2 * ix /
1096 ((cell->material_id() == id_block_2) ?
1102 (cell->material_id() == id_block_2) ? (0) : (1))]
1103 [n_cells_x_0 + ix] *
1107 else if (cell->material_id() == id_block_3 ||
1108 cell->material_id() == id_block_6)
1111 const double dx = length_b2 / n_cells_x_2;
1112 const double dy =
height / n_cells_y;
1113 const int ix =
static_cast<int>(
1114 std::round(
std::abs(node[0] - H[0]) / dx));
1116 static_cast<int>(std::round(
std::abs(node[1]) / dy));
1118 const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1119 const double alpha_x =
1120 bias_alpha(1 - (
static_cast<double>(ix)) / n_cells_x_2);
1124 const Point<2> p1(J[0] - (1 - incline_factor) * length_b2 *
1126 ((cell->material_id() == id_block_3) ?
1131 const Point<2> p2(J[0] - alpha_x * length_b2, tail_y);
1132 node = p1 * (1 - alpha_y) + p2 * alpha_y;
1152 bias_alpha(
double alpha)
const
1162 internal_create_triangulation(
1166 const AdditionalData &additional_data)
1168 MeshGenerator mesh_generator(additional_data);
1171 if (
auto *parallel_tria =
1174 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1175 else if (
auto *parallel_tria =
dynamic_cast<
1178 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1180 mesh_generator.create_triangulation(tria, periodic_faces);
1197 const AdditionalData &)
1207 const AdditionalData &additional_data)
1209 internal_create_triangulation(tria,
nullptr, additional_data);
1220 const AdditionalData &additional_data)
1222 internal_create_triangulation(tria, &periodic_faces, additional_data);
1233 const AdditionalData &additional_data)
1237 (void)additional_data;
1238 (void)periodic_faces;
1249 template <
int dim,
int spacedim>
1254 if constexpr (dim > 1)
1257 for (
const auto &f : tria.
begin()->face_indices())
1258 tria.
begin()->face(f)->set_boundary_id(f);
1264 template <
int spacedim>
1271 for (
const auto &cell : tria.cell_iterators())
1274 if (cell->center()[0] > 0)
1275 cell->set_material_id(1);
1280 template <
int dim,
int spacedim>
1285 const double epsilon)
1298 for (; face != endface; ++face)
1299 if (face->at_boundary())
1300 if (face->boundary_id() == 0)
1304 if (
std::abs(center[0] - p1[0]) < epsilon)
1305 face->set_boundary_id(0);
1306 else if (
std::abs(center[0] - p2[0]) < epsilon)
1307 face->set_boundary_id(1);
1308 else if (dim > 1 &&
std::abs(center[1] - p1[1]) < epsilon)
1309 face->set_boundary_id(2);
1310 else if (dim > 1 &&
std::abs(center[1] - p2[1]) < epsilon)
1311 face->set_boundary_id(3);
1312 else if (dim > 2 &&
std::abs(center[2] - p1[2]) < epsilon)
1313 face->set_boundary_id(4);
1314 else if (dim > 2 &&
std::abs(center[2] - p2[2]) < epsilon)
1315 face->set_boundary_id(5);
1324 for (
const auto &cell : tria.cell_iterators())
1327 for (
unsigned int d = 0;
d < dim; ++
d)
1328 if (cell->center()[d] > 0)
1330 cell->set_material_id(
id);
1339 template <
int spacdim>
1352 for (
auto cell = tria.
begin(); cell != tria.
end(); ++cell)
1355 cell->face(2)->set_all_boundary_ids(1);
1379 cell->face(4)->set_all_boundary_ids(1);
1383 cell->face(2)->set_all_boundary_ids(1);
1387 cell->face(2)->set_all_boundary_ids(1);
1391 cell->face(0)->set_all_boundary_ids(1);
1395 cell->face(2)->set_all_boundary_ids(1);
1399 cell->face(0)->set_all_boundary_ids(1);
1401 else if (tria.
n_cells() == 12)
1410 cell->face(5)->set_all_boundary_ids(1);
1413 else if (tria.
n_cells() == 96)
1420 unsigned int count = 0;
1421 for (
const auto &cell : tria.cell_iterators())
1422 if (cell->face(5)->at_boundary())
1424 cell->face(5)->set_all_boundary_ids(1);
1443 const double inner_radius,
1449 const double middle_radius =
1451 const double eps = 1
e-3 * middle_radius;
1453 for (
const auto &cell : tria.cell_iterators())
1456 const auto face = cell->face(f);
1457 if (!face->at_boundary())
1460 const double face_center_radius =
1461 (face->center(
true) - center).
norm();
1463 if (std::fabs(face->center()[0]) <
eps)
1465 face->set_boundary_id(2);
1466 for (
const unsigned int line_no : face->line_indices())
1467 if (face->line(line_no)->at_boundary())
1469 face->line(line_no)->vertex(1).
norm()) >
eps)
1470 face->line(line_no)->set_boundary_id(2);
1472 else if (std::fabs(face->center()[1]) < eps)
1474 face->set_boundary_id(3);
1475 for (
const unsigned int line_no : face->line_indices())
1476 if (face->line(line_no)->at_boundary())
1478 face->line(line_no)->vertex(1).
norm()) >
eps)
1479 face->line(line_no)->set_boundary_id(3);
1481 else if (std::fabs(face->center()[2]) < eps)
1483 face->set_boundary_id(4);
1484 for (
const unsigned int line_no : face->line_indices())
1485 if (face->line(line_no)->at_boundary())
1487 face->line(line_no)->vertex(1).
norm()) >
eps)
1488 face->line(line_no)->set_boundary_id(4);
1490 else if (face_center_radius <
1493 face->set_boundary_id(0);
1494 for (
const unsigned int line_no : face->line_indices())
1495 if (face->line(line_no)->at_boundary())
1497 face->line(line_no)->vertex(1).
norm()) <
eps)
1498 face->line(line_no)->set_boundary_id(0);
1500 else if (face_center_radius >
1503 face->set_boundary_id(1);
1504 for (
const unsigned int line_no : face->line_indices())
1505 if (face->line(line_no)->at_boundary())
1507 face->line(line_no)->vertex(1).
norm()) <
eps)
1508 face->line(line_no)->set_boundary_id(1);
1518 template <
int dim,
int spacedim>
1529 for (
unsigned int i = 0; i < dim; ++i)
1535 std::vector<Point<spacedim>> vertices(
1536 ReferenceCells::get_hypercube<dim>().n_vertices());
1544 vertices[0] = vertices[1] = p1;
1545 vertices[2] = vertices[3] = p2;
1547 vertices[1][0] = p2[0];
1548 vertices[2][0] = p1[0];
1551 vertices[0] = vertices[1] = vertices[2] = vertices[3] = p1;
1552 vertices[4] = vertices[5] = vertices[6] = vertices[7] = p2;
1554 vertices[1][0] = p2[0];
1555 vertices[2][1] = p2[1];
1556 vertices[3][0] = p2[0];
1557 vertices[3][1] = p2[1];
1559 vertices[4][0] = p1[0];
1560 vertices[4][1] = p1[1];
1561 vertices[5][1] = p1[1];
1562 vertices[6][0] = p1[0];
1570 std::vector<CellData<dim>> cells(1);
1571 for (
const unsigned int i :
1573 cells[0].vertices[i] = i;
1574 cells[0].material_id = 0;
1580 colorize_hyper_rectangle(tria);
1585 template <
int dim,
int spacedim>
1593 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1596 for (
unsigned int i = 0; i < dim; ++i)
1617 for (
unsigned int d = 0;
d < dim; ++
d)
1618 for (
unsigned int c = 1; c <= dim; ++c)
1619 vector_matrix[c - 1][d] = vertices[c][d] - vertices[0][d];
1622 "Vertices of simplex must form a right handed system"));
1626 std::vector<Point<dim>> points = vertices;
1630 for (
unsigned int i = 0; i <= dim; ++i)
1632 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1633 center += points[i];
1638 for (
unsigned int i = 1; i < dim; ++i)
1639 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1641 for (
unsigned int i = 0; i <= dim; ++i)
1642 points.push_back(1. / 3. *
1643 (points[i] + points[(i + 1) % (dim + 1)] +
1644 points[(i + 2) % (dim + 1)]));
1646 points.push_back((1. / (dim + 1)) * center);
1648 std::vector<CellData<dim>> cells(dim + 1);
1653 cells[0].vertices[0] = 0;
1654 cells[0].vertices[1] = 3;
1655 cells[0].vertices[2] = 5;
1656 cells[0].vertices[3] = 6;
1657 cells[0].material_id = 0;
1659 cells[1].vertices[0] = 3;
1660 cells[1].vertices[1] = 1;
1661 cells[1].vertices[2] = 6;
1662 cells[1].vertices[3] = 4;
1663 cells[1].material_id = 0;
1665 cells[2].vertices[0] = 5;
1666 cells[2].vertices[1] = 6;
1667 cells[2].vertices[2] = 2;
1668 cells[2].vertices[3] = 4;
1669 cells[2].material_id = 0;
1673 cells[0].vertices[0] = 0;
1674 cells[0].vertices[1] = 4;
1675 cells[0].vertices[2] = 8;
1676 cells[0].vertices[3] = 10;
1677 cells[0].vertices[4] = 7;
1678 cells[0].vertices[5] = 13;
1679 cells[0].vertices[6] = 12;
1680 cells[0].vertices[7] = 14;
1681 cells[0].material_id = 0;
1683 cells[1].vertices[0] = 4;
1684 cells[1].vertices[1] = 1;
1685 cells[1].vertices[2] = 10;
1686 cells[1].vertices[3] = 5;
1687 cells[1].vertices[4] = 13;
1688 cells[1].vertices[5] = 9;
1689 cells[1].vertices[6] = 14;
1690 cells[1].vertices[7] = 11;
1691 cells[1].material_id = 0;
1693 cells[2].vertices[0] = 8;
1694 cells[2].vertices[1] = 10;
1695 cells[2].vertices[2] = 2;
1696 cells[2].vertices[3] = 5;
1697 cells[2].vertices[4] = 12;
1698 cells[2].vertices[5] = 14;
1699 cells[2].vertices[6] = 6;
1700 cells[2].vertices[7] = 11;
1701 cells[2].material_id = 0;
1703 cells[3].vertices[0] = 7;
1704 cells[3].vertices[1] = 13;
1705 cells[3].vertices[2] = 12;
1706 cells[3].vertices[3] = 14;
1707 cells[3].vertices[4] = 3;
1708 cells[3].vertices[5] = 9;
1709 cells[3].vertices[6] = 6;
1710 cells[3].vertices[7] = 11;
1711 cells[3].material_id = 0;
1721 template <
int dim,
int spacedim>
1728 if (reference_cell == ReferenceCells::get_hypercube<dim>())
1737 std::vector<Point<spacedim>> vertices(
reference_cell.n_vertices());
1741 for (
unsigned int d = 0;
d < dim; ++
d)
1742 vertices[v][d] = this_vertex[d];
1750 std::vector<CellData<dim>> cells(1);
1753 cells[0].vertices[v] = v;
1762 const unsigned int n_cells,
1763 const unsigned int n_rotations,
1767 const unsigned int dim = 3;
1770 "More than 4 cells are needed to create a moebius grid."));
1772 ExcMessage(
"Outer and inner radius must be positive."));
1774 ExcMessage(
"Outer radius must be greater than inner radius."));
1777 std::vector<Point<dim>> vertices(4 * n_cells);
1781 for (
unsigned int i = 0; i <
n_cells; ++i)
1782 for (
unsigned int j = 0; j < 4; ++j)
1784 vertices[4 * i + j][0] =
1788 vertices[4 * i + j][1] =
1792 vertices[4 * i + j][2] =
1796 unsigned int offset = 0;
1801 static constexpr std::array<unsigned int, 8> local_vertex_numbering{
1802 {0, 1, 5, 4, 2, 3, 7, 6}};
1803 std::vector<CellData<dim>> cells(n_cells);
1804 for (
unsigned int i = 0; i <
n_cells; ++i)
1806 for (
unsigned int j = 0; j < 2; ++j)
1808 cells[i].vertices[local_vertex_numbering[0 + 4 * j]] =
1810 cells[i].vertices[local_vertex_numbering[1 + 4 * j]] =
1812 cells[i].vertices[local_vertex_numbering[2 + 4 * j]] =
1814 cells[i].vertices[local_vertex_numbering[3 + 4 * j]] =
1818 cells[i].material_id = 0;
1822 cells[
n_cells - 1].vertices[local_vertex_numbering[4]] =
1823 (0 + n_rotations) % 4;
1824 cells[
n_cells - 1].vertices[local_vertex_numbering[5]] =
1825 (3 + n_rotations) % 4;
1826 cells[
n_cells - 1].vertices[local_vertex_numbering[6]] =
1827 (2 + n_rotations) % 4;
1828 cells[
n_cells - 1].vertices[local_vertex_numbering[7]] =
1829 (1 + n_rotations) % 4;
1840 const double centerline_radius,
1841 const double inner_radius,
1845 Assert(centerline_radius > inner_radius,
1846 ExcMessage(
"The centerline radius must be greater than the "
1848 Assert(inner_radius > 0.0,
1849 ExcMessage(
"The inner radius must be positive."));
1851 const unsigned int dim = 2;
1852 const unsigned int spacedim = 3;
1853 std::vector<Point<spacedim>> vertices(16);
1855 vertices[0] =
Point<spacedim>(centerline_radius - inner_radius, 0, 0);
1857 vertices[2] =
Point<spacedim>(centerline_radius + inner_radius, 0, 0);
1859 vertices[4] =
Point<spacedim>(0, 0, centerline_radius - inner_radius);
1861 vertices[6] =
Point<spacedim>(0, 0, centerline_radius + inner_radius);
1863 vertices[8] =
Point<spacedim>(-(centerline_radius - inner_radius), 0, 0);
1865 vertices[10] =
Point<spacedim>(-(centerline_radius + inner_radius), 0, 0);
1867 vertices[12] =
Point<spacedim>(0, 0, -(centerline_radius - inner_radius));
1869 vertices[14] =
Point<spacedim>(0, 0, -(centerline_radius + inner_radius));
1872 std::vector<CellData<dim>> cells(16);
1874 cells[0].vertices[0] = 0;
1875 cells[0].vertices[1] = 4;
1876 cells[0].vertices[2] = 3;
1877 cells[0].vertices[3] = 7;
1878 cells[0].material_id = 0;
1880 cells[1].vertices[0] = 1;
1881 cells[1].vertices[1] = 5;
1882 cells[1].vertices[2] = 0;
1883 cells[1].vertices[3] = 4;
1884 cells[1].material_id = 0;
1886 cells[2].vertices[0] = 2;
1887 cells[2].vertices[1] = 6;
1888 cells[2].vertices[2] = 1;
1889 cells[2].vertices[3] = 5;
1890 cells[2].material_id = 0;
1892 cells[3].vertices[0] = 3;
1893 cells[3].vertices[1] = 7;
1894 cells[3].vertices[2] = 2;
1895 cells[3].vertices[3] = 6;
1896 cells[3].material_id = 0;
1898 cells[4].vertices[0] = 4;
1899 cells[4].vertices[1] = 8;
1900 cells[4].vertices[2] = 7;
1901 cells[4].vertices[3] = 11;
1902 cells[4].material_id = 0;
1904 cells[5].vertices[0] = 5;
1905 cells[5].vertices[1] = 9;
1906 cells[5].vertices[2] = 4;
1907 cells[5].vertices[3] = 8;
1908 cells[5].material_id = 0;
1910 cells[6].vertices[0] = 6;
1911 cells[6].vertices[1] = 10;
1912 cells[6].vertices[2] = 5;
1913 cells[6].vertices[3] = 9;
1914 cells[6].material_id = 0;
1916 cells[7].vertices[0] = 7;
1917 cells[7].vertices[1] = 11;
1918 cells[7].vertices[2] = 6;
1919 cells[7].vertices[3] = 10;
1920 cells[7].material_id = 0;
1922 cells[8].vertices[0] = 8;
1923 cells[8].vertices[1] = 12;
1924 cells[8].vertices[2] = 11;
1925 cells[8].vertices[3] = 15;
1926 cells[8].material_id = 0;
1928 cells[9].vertices[0] = 9;
1929 cells[9].vertices[1] = 13;
1930 cells[9].vertices[2] = 8;
1931 cells[9].vertices[3] = 12;
1932 cells[9].material_id = 0;
1934 cells[10].vertices[0] = 10;
1935 cells[10].vertices[1] = 14;
1936 cells[10].vertices[2] = 9;
1937 cells[10].vertices[3] = 13;
1938 cells[10].material_id = 0;
1940 cells[11].vertices[0] = 11;
1941 cells[11].vertices[1] = 15;
1942 cells[11].vertices[2] = 10;
1943 cells[11].vertices[3] = 14;
1944 cells[11].material_id = 0;
1946 cells[12].vertices[0] = 12;
1947 cells[12].vertices[1] = 0;
1948 cells[12].vertices[2] = 15;
1949 cells[12].vertices[3] = 3;
1950 cells[12].material_id = 0;
1952 cells[13].vertices[0] = 13;
1953 cells[13].vertices[1] = 1;
1954 cells[13].vertices[2] = 12;
1955 cells[13].vertices[3] = 0;
1956 cells[13].material_id = 0;
1958 cells[14].vertices[0] = 14;
1959 cells[14].vertices[1] = 2;
1960 cells[14].vertices[2] = 13;
1961 cells[14].vertices[3] = 1;
1962 cells[14].material_id = 0;
1964 cells[15].vertices[0] = 15;
1965 cells[15].vertices[1] = 3;
1966 cells[15].vertices[2] = 14;
1967 cells[15].vertices[3] = 2;
1968 cells[15].material_id = 0;
1981 static constexpr int circle_cell_vertices[5][4] = {{0, 1, 2, 3},
1993 const double centerline_radius,
1994 const double inner_radius,
1995 const unsigned int n_cells_toroidal,
1998 Assert(centerline_radius > inner_radius,
1999 ExcMessage(
"The centerline radius must be greater than the "
2001 Assert(inner_radius > 0.0,
2002 ExcMessage(
"The inner radius must be positive."));
2004 ExcMessage(
"Number of cells in toroidal direction has "
2005 "to be at least 3 for a torus of polar extent 2*pi."));
2011 const double a = 1. / (1 +
std::sqrt(2.0));
2014 const unsigned int additional_layer =
2018 const unsigned int n_point_layers_toroidal =
2019 n_cells_toroidal + additional_layer;
2020 std::vector<Point<3>> vertices(8 * n_point_layers_toroidal);
2032 const double phi_cell = phi / n_cells_toroidal;
2033 for (
unsigned int c = 1; c < n_point_layers_toroidal; ++c)
2035 for (
unsigned int v = 0; v < 8; ++v)
2037 const double inner_radius_2d = vertices[v][0];
2038 vertices[8 * c + v][0] = inner_radius_2d *
std::cos(phi_cell * c);
2039 vertices[8 * c + v][1] = vertices[v][1];
2040 vertices[8 * c + v][2] = inner_radius_2d *
std::sin(phi_cell * c);
2045 std::vector<CellData<3>> cells(5 * n_cells_toroidal);
2046 for (
unsigned int c = 0; c < n_cells_toroidal; ++c)
2048 for (
unsigned int j = 0; j < 2; ++j)
2050 const unsigned int offset =
2051 (8 * (c + j)) % (8 * n_point_layers_toroidal);
2054 for (
unsigned int c2 = 0; c2 < 5; ++c2)
2055 for (
unsigned int i = 0; i < 4; ++i)
2056 cells[5 * c + c2].vertices[i + j * 4] =
2057 offset + circle_cell_vertices[c2][i];
2060 cells[5 * c].material_id = 0;
2062 cells[5 * c + 1].material_id = 0;
2063 cells[5 * c + 2].material_id = 1;
2064 cells[5 * c + 3].material_id = 0;
2065 cells[5 * c + 4].material_id = 0;
2073 for (
const auto &cell : tria.cell_iterators())
2076 for (
const unsigned int f : cell->face_indices())
2080 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2082 cell->face(f)->set_all_manifold_ids(1);
2087 if (cell->material_id() == 1)
2089 cell->set_all_manifold_ids(2);
2091 cell->set_material_id(0);
2108 template <
int dim,
int spacedim>
2114 Assert(vertices.size() == ReferenceCells::get_hypercube<dim>().n_vertices(),
2124 cell->vertex(i) = vertices[i];
2130 "The volume of the cell is not greater than zero. "
2131 "This could be due to the wrong ordering of the vertices."));
2162 std::array<Tensor<1, 2>, 2> edges;
2163 edges[0] = corners[0];
2164 edges[1] = corners[1];
2165 std::vector<unsigned int> subdivisions;
2166 subdivided_parallelepiped<2, 2>(
2167 tria, origin, edges, subdivisions,
colorize);
2178 unsigned int n_subdivisions[dim];
2179 for (
unsigned int i = 0; i < dim; ++i)
2180 n_subdivisions[i] = 1;
2189 const unsigned int n_subdivisions,
2195 unsigned int n_subdivisions_[dim];
2196 for (
unsigned int i = 0; i < dim; ++i)
2197 n_subdivisions_[i] = n_subdivisions;
2207 const unsigned int (&n_subdivisions)[dim],
2209 const unsigned int *n_subdivisions,
2215 std::vector<unsigned int> subdivisions;
2216 std::array<Tensor<1, dim>, dim> edges;
2217 for (
unsigned int i = 0; i < dim; ++i)
2219 subdivisions.push_back(n_subdivisions[i]);
2220 edges[i] = corners[i];
2223 subdivided_parallelepiped<dim, dim>(
2224 tria, origin, edges, subdivisions,
colorize);
2230 template <
int dim,
int spacedim>
2235 const std::vector<unsigned int> &subdivisions,
2238 std::vector<unsigned int> compute_subdivisions = subdivisions;
2239 if (compute_subdivisions.empty())
2241 compute_subdivisions.resize(dim, 1);
2244 Assert(compute_subdivisions.size() == dim,
2245 ExcMessage(
"One subdivision must be provided for each dimension."));
2247 for (
unsigned int i = 0; i < dim; ++i)
2249 Assert(compute_subdivisions[i] > 0,
2250 ExcInvalidRepetitions(subdivisions[i]));
2252 edges[i].
norm() > 0,
2254 "Edges in subdivided_parallelepiped() must not be degenerate."));
2262 bool twisted_data =
false;
2267 twisted_data = (edges[0][0] < 0);
2274 const double plane_normal =
2275 edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2276 twisted_data = (plane_normal < 0.0);
2285 (edges[0].
norm() * edges[1].
norm()) -
2288 "Edges in subdivided_parallelepiped() must point in"
2289 " different directions."));
2291 cross_product_3d(edges[0], edges[1]);
2305 twisted_data = (plane_normal * edges[2] < 0.0);
2313 ExcInvalidInputOrientation(
2314 "The triangulation you are trying to create will consist of cells"
2315 " with negative measures. This is usually the result of input data"
2316 " that does not define a right-handed coordinate system. The usual"
2317 " fix for this is to ensure that in 1d the given point is to the"
2318 " right of the origin (or the given edge tensor is positive), in 2d"
2319 " that the two edges (and their cross product) obey the right-hand"
2320 " rule (which may usually be done by switching the order of the"
2321 " points or edge tensors), or in 3d that the edges form a"
2322 " right-handed coordinate system (which may also be accomplished by"
2323 " switching the order of the first two points or edge tensors)."));
2326 for (
unsigned int i = 0; i < dim; ++i)
2327 for (
unsigned int j = i + 1; j < dim; ++j)
2328 Assert((edges[i] != edges[j]),
2330 "Degenerate edges of subdivided_parallelepiped encountered."));
2333 std::vector<Point<spacedim>> points;
2338 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2339 points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2343 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2344 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2345 points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2346 edges[1] / compute_subdivisions[1] * y);
2350 for (
unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2351 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2352 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2353 points.push_back(origin +
2354 edges[0] / compute_subdivisions[0] * x +
2355 edges[1] / compute_subdivisions[1] * y +
2356 edges[2] / compute_subdivisions[2] * z);
2365 for (
unsigned int i = 0; i < dim; ++i)
2366 n_cells *= compute_subdivisions[i];
2367 std::vector<CellData<dim>> cells(n_cells);
2373 for (
unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2375 cells[x].vertices[0] = x;
2376 cells[x].vertices[1] = x + 1;
2379 cells[x].material_id = 0;
2386 const unsigned int n_dy = compute_subdivisions[1];
2387 const unsigned int n_dx = compute_subdivisions[0];
2389 for (
unsigned int y = 0; y < n_dy; ++y)
2390 for (
unsigned int x = 0; x < n_dx; ++x)
2392 const unsigned int c = y * n_dx + x;
2393 cells[c].vertices[0] = y * (n_dx + 1) + x;
2394 cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2395 cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2396 cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2399 cells[c].material_id = 0;
2407 const unsigned int n_dz = compute_subdivisions[2];
2408 const unsigned int n_dy = compute_subdivisions[1];
2409 const unsigned int n_dx = compute_subdivisions[0];
2411 for (
unsigned int z = 0; z < n_dz; ++z)
2412 for (
unsigned int y = 0; y < n_dy; ++y)
2413 for (
unsigned int x = 0; x < n_dx; ++x)
2415 const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2417 cells[c].vertices[0] =
2418 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2419 cells[c].vertices[1] =
2420 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2421 cells[c].vertices[2] =
2422 z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2423 cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2424 (y + 1) * (n_dx + 1) + x + 1;
2425 cells[c].vertices[4] =
2426 (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2427 cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2428 y * (n_dx + 1) + x + 1;
2429 cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2430 (y + 1) * (n_dx + 1) + x;
2431 cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2432 (y + 1) * (n_dx + 1) + x + 1;
2435 cells[c].material_id = 0;
2453 for (
const auto &cell : tria.active_cell_iterators())
2454 for (
const unsigned
int face : cell->face_indices())
2456 if (cell->face(face)->at_boundary())
2457 cell->face(face)->set_boundary_id(face);
2463 template <
int dim,
int spacedim>
2466 const unsigned int repetitions,
2471 Assert(repetitions >= 1, ExcInvalidRepetitions(repetitions));
2473 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2476 for (
unsigned int i = 0; i < dim; ++i)
2482 std::vector<unsigned int> reps(dim, repetitions);
2488 template <
int dim,
int spacedim>
2491 const std::vector<unsigned int> &repetitions,
2496 Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim));
2502 for (
unsigned int i = 0; i < dim; ++i)
2509 std::array<Point<spacedim>, dim> delta;
2510 for (
unsigned int i = 0; i < dim; ++i)
2512 Assert(repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
2514 delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2518 "The first dim entries of coordinates of p1 and p2 need to be different."));
2522 std::vector<Point<spacedim>> points;
2526 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2527 points.push_back(p1 + x * delta[0]);
2531 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2532 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2533 points.push_back(p1 + x * delta[0] + y * delta[1]);
2537 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
2538 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2539 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2540 points.push_back(p1 + x * delta[0] + y * delta[1] +
2549 std::vector<CellData<dim>> cells;
2554 cells.resize(repetitions[0]);
2555 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2557 cells[x].vertices[0] = x;
2558 cells[x].vertices[1] = x + 1;
2559 cells[x].material_id = 0;
2566 cells.resize(repetitions[1] * repetitions[0]);
2567 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2568 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2570 const unsigned int c = x + y * repetitions[0];
2571 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2572 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2573 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2574 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2575 cells[c].material_id = 0;
2582 const unsigned int n_x = (repetitions[0] + 1);
2583 const unsigned int n_xy =
2584 (repetitions[0] + 1) * (repetitions[1] + 1);
2586 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2587 for (
unsigned int z = 0; z < repetitions[2]; ++z)
2588 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2589 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2591 const unsigned int c = x + y * repetitions[0] +
2592 z * repetitions[0] * repetitions[1];
2593 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2594 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2595 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2596 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2597 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2598 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2599 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2600 cells[c].vertices[7] =
2601 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2602 cells[c].material_id = 0;
2623 double epsilon = std::numeric_limits<double>::max();
2624 for (
unsigned int i = 0; i < dim; ++i)
2625 epsilon =
std::min(epsilon, 0.01 * delta[i][i]);
2628 "The distance between corner points must be positive."));
2632 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2641 const std::vector<std::vector<double>> &step_sz,
2646 Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
2656 std::vector<std::vector<double>> step_sizes(step_sz);
2658 for (
unsigned int i = 0; i < dim; ++i)
2662 std::swap(p1[i], p2[i]);
2663 std::reverse(step_sizes[i].
begin(), step_sizes[i].
end());
2669 for (
unsigned int j = 0; j < step_sizes.at(i).size(); ++j)
2670 x += step_sizes[i][j];
2671 Assert(std::fabs(x - (p2[i] - p1[i])) <= 1e-12 * std::fabs(x),
2673 "The sequence of step sizes in coordinate direction " +
2675 " must be equal to the distance of the two given "
2676 "points in this coordinate direction."));
2683 std::vector<Point<dim>> points;
2689 for (
unsigned int i = 0;; ++i)
2698 if (i == step_sizes[0].
size())
2701 x += step_sizes[0][i];
2709 for (
unsigned int j = 0;; ++j)
2712 for (
unsigned int i = 0;; ++i)
2714 points.push_back(
Point<dim>(p1[0] + x, p1[1] + y));
2715 if (i == step_sizes[0].
size())
2718 x += step_sizes[0][i];
2721 if (j == step_sizes[1].
size())
2724 y += step_sizes[1][j];
2731 for (
unsigned int k = 0;; ++k)
2734 for (
unsigned int j = 0;; ++j)
2737 for (
unsigned int i = 0;; ++i)
2740 Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2741 if (i == step_sizes[0].
size())
2744 x += step_sizes[0][i];
2747 if (j == step_sizes[1].
size())
2750 y += step_sizes[1][j];
2753 if (k == step_sizes[2].
size())
2756 z += step_sizes[2][k];
2767 std::vector<CellData<dim>> cells;
2772 cells.resize(step_sizes[0].
size());
2773 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2775 cells[x].vertices[0] = x;
2776 cells[x].vertices[1] = x + 1;
2777 cells[x].material_id = 0;
2784 cells.resize(step_sizes[1].
size() * step_sizes[0].
size());
2785 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2786 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2788 const unsigned int c = x + y * step_sizes[0].size();
2789 cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2790 cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2791 cells[c].vertices[2] =
2792 (y + 1) * (step_sizes[0].
size() + 1) + x;
2793 cells[c].vertices[3] =
2794 (y + 1) * (step_sizes[0].
size() + 1) + x + 1;
2795 cells[c].material_id = 0;
2802 const unsigned int n_x = (step_sizes[0].size() + 1);
2803 const unsigned int n_xy =
2804 (step_sizes[0].size() + 1) * (step_sizes[1].
size() + 1);
2806 cells.resize(step_sizes[2].
size() * step_sizes[1].size() *
2807 step_sizes[0].size());
2808 for (
unsigned int z = 0; z < step_sizes[2].size(); ++z)
2809 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2810 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2812 const unsigned int c =
2813 x + y * step_sizes[0].size() +
2814 z * step_sizes[0].size() * step_sizes[1].size();
2815 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2816 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2817 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2818 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2819 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2820 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2821 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2822 cells[c].vertices[7] =
2823 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2824 cells[c].material_id = 0;
2846 *std::min_element(step_sizes[0].
begin(), step_sizes[0].
end());
2847 for (
unsigned int i = 1; i < dim; ++i)
2849 *std::min_element(step_sizes[i].
begin(),
2850 step_sizes[i].
end()));
2851 const double epsilon = 0.01 * min_size;
2855 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2864 const std::vector<std::vector<double>> &spacing,
2869 Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
2873 Assert(spacing[0].
size() == n_cells, ExcInvalidRepetitionsDimension(1));
2875 double delta = std::numeric_limits<double>::max();
2876 for (
unsigned int i = 0; i <
n_cells; ++i)
2878 Assert(spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
2879 delta =
std::min(delta, spacing[0][i]);
2883 std::vector<Point<1>> points;
2885 for (
unsigned int x = 0; x <=
n_cells; ++x)
2887 points.emplace_back(ax);
2889 ax += spacing[0][x];
2892 unsigned int n_val_cells = 0;
2893 for (
unsigned int i = 0; i <
n_cells; ++i)
2897 std::vector<CellData<1>> cells(n_val_cells);
2898 unsigned int id = 0;
2899 for (
unsigned int x = 0; x <
n_cells; ++x)
2902 cells[id].vertices[0] = x;
2903 cells[id].vertices[1] = x + 1;
2922 const std::vector<std::vector<double>> &spacing,
2927 Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
2929 std::vector<unsigned int> repetitions(2);
2930 double delta = std::numeric_limits<double>::max();
2931 for (
unsigned int i = 0; i < 2; ++i)
2933 repetitions[i] = spacing[i].size();
2934 for (
unsigned int j = 0; j < repetitions[i]; ++j)
2936 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
2937 delta =
std::min(delta, spacing[i][j]);
2940 ExcInvalidRepetitionsDimension(i));
2944 std::vector<Point<2>> points;
2946 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2949 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2951 points.emplace_back(ax, ay);
2952 if (x < repetitions[0])
2953 ax += spacing[0][x];
2955 if (y < repetitions[1])
2956 ay += spacing[1][y];
2960 unsigned int n_val_cells = 0;
2961 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
2962 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
2966 std::vector<CellData<2>> cells(n_val_cells);
2967 unsigned int id = 0;
2968 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2969 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2972 cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
2973 cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2974 cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2975 cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2989 double eps = 0.01 * delta;
2991 for (; cell != endc; ++cell)
2993 Point<2> cell_center = cell->center();
2994 for (
const unsigned int f : cell->face_indices())
2997 Point<2> face_center = cell->face(f)->center();
2998 for (
unsigned int i = 0; i < 2; ++i)
3000 if (face_center[i] < cell_center[i] - eps)
3001 cell->face(f)->set_boundary_id(i * 2);
3002 if (face_center[i] > cell_center[i] + eps)
3003 cell->face(f)->set_boundary_id(i * 2 + 1);
3014 const std::vector<std::vector<double>> &spacing,
3019 const unsigned int dim = 3;
3021 Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
3023 std::array<unsigned int, dim> repetitions;
3024 double delta = std::numeric_limits<double>::max();
3025 for (
unsigned int i = 0; i < dim; ++i)
3027 repetitions[i] = spacing[i].size();
3028 for (
unsigned int j = 0; j < repetitions[i]; ++j)
3030 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
3031 delta =
std::min(delta, spacing[i][j]);
3034 ExcInvalidRepetitionsDimension(i));
3038 std::vector<Point<dim>> points;
3040 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3043 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3046 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3048 points.emplace_back(ax, ay, az);
3049 if (x < repetitions[0])
3050 ax += spacing[0][x];
3052 if (y < repetitions[1])
3053 ay += spacing[1][y];
3055 if (z < repetitions[2])
3056 az += spacing[2][z];
3060 unsigned int n_val_cells = 0;
3061 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
3062 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
3063 for (
unsigned int k = 0; k <
material_id.size(2); ++k)
3067 std::vector<CellData<dim>> cells(n_val_cells);
3068 unsigned int id = 0;
3069 const unsigned int n_x = (repetitions[0] + 1);
3070 const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3071 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3072 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3073 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3076 cells[id].vertices[0] = z * n_xy + y * n_x + x;
3077 cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3078 cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3079 cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3080 cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3081 cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3082 cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3083 cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3097 double eps = 0.01 * delta;
3100 for (; cell != endc; ++cell)
3103 for (
auto f : cell->face_indices())
3106 Point<dim> face_center = cell->face(f)->center();
3107 for (
unsigned int i = 0; i < dim; ++i)
3109 if (face_center[i] < cell_center[i] - eps)
3110 cell->face(f)->set_boundary_id(i * 2);
3111 if (face_center[i] > cell_center[i] + eps)
3112 cell->face(f)->set_boundary_id(i * 2 + 1);
3121 template <
int dim,
int spacedim>
3124 const std::vector<unsigned int> &holes)
3135 std::array<Point<spacedim>, dim> delta;
3136 std::array<unsigned int, dim> repetitions;
3137 for (
unsigned int i = 0; i < dim; ++i)
3140 ExcMessage(
"At least one hole needed in each direction"));
3141 repetitions[i] = 2 * holes[i] + 1;
3147 std::vector<Point<spacedim>> points;
3151 points.reserve(repetitions[0] + 1);
3152 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3153 points.push_back(origin + x * delta[0]);
3157 points.reserve((repetitions[0] + 1) * (repetitions[1] + 1));
3158 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3159 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3160 points.push_back(origin + x * delta[0] + y * delta[1]);
3164 points.reserve((repetitions[0] + 1) * (repetitions[1] + 1) *
3165 (repetitions[2] + 1));
3166 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3167 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3168 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3169 points.push_back(origin + x * delta[0] + y * delta[1] +
3179 std::vector<CellData<dim>> cells;
3184 const unsigned int n_x = (repetitions[0] + 1);
3186 cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3188 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3189 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3193 if ((x % 2 == 1) && (y % 2 == 1))
3197 cells[c].vertices[0] = y * n_x + x;
3198 cells[c].vertices[1] = y * n_x + x + 1;
3199 cells[c].vertices[2] = (y + 1) * n_x + x;
3200 cells[c].vertices[3] = (y + 1) * n_x + x + 1;
3201 cells[c].material_id = 0;
3210 const unsigned int n_x = (repetitions[0] + 1);
3211 const unsigned int n_xy =
3212 (repetitions[0] + 1) * (repetitions[1] + 1);
3216 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3219 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3220 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3221 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3226 if ((x % 2 == 1) && (y % 2 == 1) && z % 2 == 1)
3230 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3231 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3232 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3233 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3234 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3235 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3236 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3237 cells[c].vertices[7] =
3238 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3239 cells[c].material_id = 0;
3259 template <
int dim,
int spacedim>
3263 const double voxel_edge_length)
3268 std::array<unsigned int, dim> voxels_per_direction;
3269 for (
unsigned int d = 0;
d < dim; ++
d)
3271 Assert(present_voxels.size()[d] > 0,
3273 "The 2d or 3d array that describes which pixels/voxels "
3274 "are present in the domain must have a non-zero extent in"
3275 "in coordinate direction."));
3276 voxels_per_direction[
d] = present_voxels.size()[
d];
3279 std::vector<Point<spacedim>> points;
3283 points.reserve(voxels_per_direction[0] + 1);
3284 for (
unsigned int x = 0; x < voxels_per_direction[0] + 1; ++x)
3285 points.emplace_back(x * voxel_edge_length);
3289 points.reserve((voxels_per_direction[0] + 1) *
3290 (voxels_per_direction[1] + 1));
3291 for (
unsigned int y = 0; y < voxels_per_direction[1] + 1; ++y)
3292 for (
unsigned int x = 0; x < voxels_per_direction[0] + 1; ++x)
3293 points.emplace_back(x * voxel_edge_length, y * voxel_edge_length);
3297 points.reserve((voxels_per_direction[0] + 1) *
3298 (voxels_per_direction[1] + 1) *
3299 (voxels_per_direction[2] + 1));
3300 for (
unsigned int z = 0; z < voxels_per_direction[2] + 1; ++z)
3301 for (
unsigned int y = 0; y < voxels_per_direction[1] + 1; ++y)
3302 for (
unsigned int x = 0; x < voxels_per_direction[0] + 1; ++x)
3303 points.emplace_back(x * voxel_edge_length,
3304 y * voxel_edge_length,
3305 z * voxel_edge_length);
3313 unsigned int n_voxels = 0;
3314 if constexpr (dim == 2)
3316 for (
unsigned int x = 0; x < voxels_per_direction[0]; ++x)
3317 for (
unsigned int y = 0; y < voxels_per_direction[1]; ++y)
3318 if (present_voxels[x][y] ==
true)
3321 else if constexpr (dim == 3)
3323 for (
unsigned int x = 0; x < voxels_per_direction[0]; ++x)
3324 for (
unsigned int y = 0; y < voxels_per_direction[1]; ++y)
3325 for (
unsigned int z = 0; z < voxels_per_direction[2]; ++z)
3326 if (present_voxels[x][y][z] ==
true)
3332 std::vector<CellData<dim>> cells(n_voxels);
3333 if constexpr (dim == 2)
3335 const unsigned int n_x = (voxels_per_direction[0] + 1);
3338 for (
unsigned int y = 0; y < voxels_per_direction[1]; ++y)
3339 for (
unsigned int x = 0; x < voxels_per_direction[0]; ++x)
3340 if (present_voxels[x][y] ==
true)
3342 cells[c].vertices[0] = y * n_x + x;
3343 cells[c].vertices[1] = y * n_x + x + 1;
3344 cells[c].vertices[2] = (y + 1) * n_x + x;
3345 cells[c].vertices[3] = (y + 1) * n_x + x + 1;
3346 cells[c].material_id = 0;
3351 else if constexpr (dim == 3)
3353 const unsigned int n_x = (voxels_per_direction[0] + 1);
3354 const unsigned int n_xy =
3355 (voxels_per_direction[0] + 1) * (voxels_per_direction[1] + 1);
3358 for (
unsigned int z = 0; z < voxels_per_direction[2]; ++z)
3359 for (
unsigned int y = 0; y < voxels_per_direction[1]; ++y)
3360 for (
unsigned int x = 0; x < voxels_per_direction[0]; ++x)
3361 if (present_voxels[x][y][z] ==
true)
3364 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3365 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3366 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3367 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3368 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3369 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3370 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3371 cells[c].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3372 cells[c].material_id = 0;
3403 const unsigned int ,
3415 const unsigned int ,
3439 template <
int dim,
int spacedim>
3443 double length = std::numeric_limits<double>::max();
3444 for (
const auto &cell : triangulation.active_cell_iterators())
3445 for (unsigned
int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3446 length =
std::min(length, cell->line(n)->diameter());
3456 const double inner_radius,
3458 const double pad_bottom,
3459 const double pad_top,
3460 const double pad_left,
3461 const double pad_right,
3466 const unsigned int ,
3469 const bool with_padding =
3470 pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3480 double length = std::numeric_limits<double>::max();
3481 for (
const auto &cell : tria.active_cell_iterators())
3482 for (unsigned
int n = 0; n < cell->n_lines(); ++n)
3483 length =
std::min(length, cell->line(n)->diameter());
3489 Triangulation<2> &cylinder_tria = with_padding ? cylinder_tria_maybe : tria;
3498 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3499 cell->set_manifold_id(tfi_manifold_id);
3508 auto add_sizes = [](std::vector<double> &step_sizes,
3509 const double padding,
3510 const double h) ->
void {
3513 const auto rounded =
3514 static_cast<unsigned int>(std::round(padding / h));
3517 const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3518 for (
unsigned int i = 0; i < num; ++i)
3519 step_sizes.push_back(padding / num);
3522 std::vector<std::vector<double>> step_sizes(2);
3543 bulk_tria, step_sizes, bl, tr,
colorize);
3546 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3547 for (
const auto &cell : bulk_tria.active_cell_iterators())
3549 cells_to_remove.
insert(cell);
3553 bulk_tria, cells_to_remove, tria_without_cylinder);
3555 const double tolerance =
3556 std::min(min_line_length(tria_without_cylinder),
3557 min_line_length(cylinder_tria)) /
3561 tria_without_cylinder, cylinder_tria, tria, tolerance,
true);
3565 for (
const auto &cell : tria.active_cell_iterators())
3569 if (cell->manifold_id() == tfi_manifold_id)
3571 for (
const unsigned int face_n : cell->face_indices())
3573 const auto &face = cell->face(face_n);
3574 if (face->at_boundary() &&
3575 internal::point_in_2d_box(face->center(),
3578 face->set_manifold_id(polar_manifold_id);
3580 face->set_manifold_id(tfi_manifold_id);
3591 static constexpr double tol =
3592 std::numeric_limits<double>::epsilon() * 10000;
3594 for (
const auto &cell : tria.active_cell_iterators())
3595 for (
const unsigned
int face_n : cell->face_indices())
3597 const auto face = cell->face(face_n);
3598 if (face->at_boundary())
3600 const Point<2> center = face->center();
3603 face->set_boundary_id(0);
3606 face->set_boundary_id(1);
3609 face->set_boundary_id(2);
3612 face->set_boundary_id(3);
3616 Assert(cell->manifold_id() == tfi_manifold_id,
3618 face->set_boundary_id(4);
3639 const double inner_radius,
3641 const double pad_bottom,
3642 const double pad_top,
3643 const double pad_left,
3644 const double pad_right,
3649 const unsigned int n_slices,
3660 Point<2>(new_center[0], new_center[1]),
3683 tria.
set_manifold(polar_manifold_id, cylindrical_manifold);
3692 const double shell_region_width,
3693 const unsigned int n_shells,
3694 const double skewness,
3697 Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3698 ExcMessage(
"The width of the shell region must be less than 0.05 "
3699 "(and preferably close to 0.03)"));
3733 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3735 for (
const auto &cell : bulk_tria.active_cell_iterators())
3737 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3738 cells_to_remove.insert(cell);
3742 for (
const unsigned int vertex_n :
3744 if (cell->vertex(vertex_n) ==
Point<2>())
3748 cylinder_triangulation_offset =
3749 2.0 * (cell->vertex(3) -
Point<2>());
3756 bulk_tria, cells_to_remove, tria_without_cylinder);
3763 0.05 + shell_region_width,
3768 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3771 if (
std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1e-10)
3772 cell->vertex(vertex_n)[0] = -0.1;
3773 else if (
std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1e-10)
3774 cell->vertex(vertex_n)[0] = 0.1;
3778 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3780 cell->set_manifold_id(tfi_manifold_id);
3781 for (
const unsigned int face_n :
GeometryInfo<2>::face_indices())
3782 if (!cell->face(face_n)->at_boundary())
3783 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3785 if (0.0 < shell_region_width)
3788 ExcMessage(
"If the shell region has positive width then "
3789 "there must be at least one shell."));
3794 0.05 + shell_region_width,
3801 const double vertex_tolerance =
3802 std::min(internal::minimal_vertex_distance(shell_tria),
3803 internal::minimal_vertex_distance(cylinder_tria)) *
3809 shell_tria, cylinder_tria, temp, vertex_tolerance,
true);
3810 cylinder_tria = std::move(temp);
3816 const double vertex_tolerance =
3817 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3818 internal::minimal_vertex_distance(cylinder_tria)) /
3821 tria_without_cylinder, cylinder_tria, tria, vertex_tolerance,
true);
3834 const double shift =
3835 std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3836 for (
const auto &cell : tria.active_cell_iterators())
3838 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3839 cell->vertex(v) =
Point<2>(0.2 -
shift, 0.205);
3840 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1e-10)
3841 cell->vertex(v) =
Point<2>(0.2 + shift, 0.205);
3842 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1e-10)
3843 cell->vertex(v) =
Point<2>(0.2, 0.2 - shift);
3844 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1e-10)
3845 cell->vertex(v) =
Point<2>(0.2, 0.2 + shift);
3850 for (
const auto &cell : tria.active_cell_iterators())
3852 cell->set_all_manifold_ids(polar_manifold_id);
3856 for (
const auto &cell : tria.active_cell_iterators())
3863 std::vector<Point<2> *> cylinder_pointers;
3864 for (
const auto &face : tria.active_face_iterators())
3867 cylinder_pointers.push_back(&face->vertex(0));
3868 cylinder_pointers.push_back(&face->vertex(1));
3871 std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3872 cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3873 cylinder_pointers.end()),
3874 cylinder_pointers.end());
3878 for (
const Point<2> *
const ptr : cylinder_pointers)
3879 center += *ptr / double(cylinder_pointers.
size());
3882 for (
Point<2> *
const ptr : cylinder_pointers)
3883 *ptr +=
Point<2>(0.2, 0.2) - center;
3895 for (
const auto &face : tria.active_face_iterators())
3896 if (face->at_boundary())
3898 const Point<2> center = face->center();
3900 if (
std::abs(center[0] - 0.0) < 1e-10)
3901 face->set_boundary_id(0);
3903 else if (
std::abs(center[0] - 2.2) < 1e-10)
3904 face->set_boundary_id(1);
3906 else if (face->manifold_id() == polar_manifold_id)
3907 face->set_boundary_id(2);
3912 std::abs(center[1] - 0.41) < 1.0e-10,
3914 face->set_boundary_id(3);
3924 const double shell_region_width,
3925 const unsigned int n_shells,
3926 const double skewness,
3931 tria_2, shell_region_width, n_shells, skewness,
colorize);
3951 tria.
set_manifold(cylindrical_manifold_id, cylindrical_manifold);
3958 for (
const auto &face : tria.active_face_iterators())
3960 face->set_boundary_id(3);
3966 const std::vector<unsigned int> &,
3982 const std::vector<unsigned int> &lengths_and_heights,
3985 const double shell_region_radius,
3986 const unsigned int n_shells,
3987 const double skewness,
3988 const bool use_transfinite_region,
3995 const double radius = 0.5;
3996 const double box_radius = 1;
4003 const unsigned int length_pre = lengths_and_heights[0];
4004 const unsigned int length_post = lengths_and_heights[1];
4005 const unsigned int height_below = lengths_and_heights[2];
4006 const unsigned int height_above = lengths_and_heights[3];
4008 const unsigned int length_repetitions = length_pre + length_post;
4009 const unsigned int height_repetitions = height_above + height_below;
4016 {(length_repetitions), height_repetitions},
4017 Point<2>(-
double(length_pre), -
double(height_below)),
4018 Point<2>(
double(length_post),
double(height_above)));
4036 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
4037 for (
const auto &cell : bulk_tria.active_cell_iterators())
4039 if ((cell->center() -
Point<2>(0., 0.)).norm() < 1.1 * box_radius)
4040 cells_to_remove.insert(cell);
4045 bulk_tria, cells_to_remove, tria_without_cylinder);
4052 shell_region_radius,
4056 for (
const auto &cell : cylinder_tria.active_cell_iterators())
4058 cell->set_manifold_id(tfi_manifold_id);
4059 for (
const unsigned int face_n :
GeometryInfo<2>::face_indices())
4060 if (!cell->face(face_n)->at_boundary())
4061 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
4066 if (
radius < shell_region_radius)
4069 ExcMessage(
"If the shell region has positive width then "
4070 "there must be at least one shell."));
4075 shell_region_radius,
4082 const double vertex_tolerance =
4083 std::min(internal::minimal_vertex_distance(shell_tria),
4084 internal::minimal_vertex_distance(cylinder_tria)) *
4090 shell_tria, cylinder_tria, temp, vertex_tolerance,
true);
4091 cylinder_tria = std::move(temp);
4096 const double vertex_tolerance =
4097 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
4098 internal::minimal_vertex_distance(cylinder_tria)) /
4102 tria_without_cylinder, cylinder_tria, tria, vertex_tolerance,
true);
4106 for (
const auto &cell : tria.active_cell_iterators())
4108 cell->set_all_manifold_ids(polar_manifold_id);
4112 for (
const auto &cell : tria.active_cell_iterators())
4121 if (use_transfinite_region)
4133 for (
const auto &face : tria.active_face_iterators())
4134 if (face->at_boundary())
4136 const Point<2> center = face->center();
4138 if (
std::abs(center[0] - (-
static_cast<double>(length_pre))) <
4140 face->set_boundary_id(0);
4142 else if (
std::abs(center[0] -
static_cast<double>(length_post)) <
4144 face->set_boundary_id(1);
4146 else if (face->manifold_id() == polar_manifold_id)
4147 face->set_boundary_id(2);
4150 (-
static_cast<double>(height_below))) < 1e-10)
4151 face->set_boundary_id(3);
4154 face->set_boundary_id(4);
4162 const std::vector<unsigned int> &lengths_and_heights,
4164 unsigned int depth_division,
4165 const double shell_region_radius,
4166 const unsigned int n_shells,
4167 const double skewness,
4168 const bool use_transfinite_region,
4173 lengths_and_heights,
4176 shell_region_radius,
4179 use_transfinite_region,
4198 tria.
set_manifold(cylindrical_manifold_id, cylindrical_manifold);
4200 if (use_transfinite_region)
4215 template <
int dim,
int spacedim>
4218 const std::vector<unsigned int> &sizes,
4228 for (
unsigned int d = 0;
d < dim; ++
d)
4231 std::vector<Point<spacedim>> points;
4233 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
4236 std::vector<CellData<dim>> cells(n_cells);
4241 for (
unsigned int d = 0;
d < dim; ++
d)
4242 p[d] = 0.5 * dimensions[d] *
4245 points.push_back(p);
4246 cells[0].vertices[i] = i;
4248 cells[0].material_id = 0;
4253 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
4259 for (
unsigned int j = 0; j < sizes[face]; ++j, ++
cell_index)
4261 const unsigned int last_cell = (j == 0) ? 0U : (
cell_index - 1);
4263 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
4266 const unsigned int cellv =
4268 const unsigned int ocellv =
4272 cells[last_cell].vertices[cellv];
4275 cells[
cell_index].vertices[cellv] = points.size();
4280 points.push_back(p);
4429 const double thickness,
4433 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4435 std::vector<Point<2>> vertices(16);
4437 coords[0] = left - thickness;
4440 coords[3] = right + thickness;
4443 for (
const double y : coords)
4445 vertices[k++] =
Point<2>(x, y);
4449 std::vector<CellData<2>> cells(9);
4451 for (
unsigned int i0 = 0; i0 < 3; ++i0)
4452 for (
unsigned int i1 = 0; i1 < 3; ++i1)
4454 cells[k].vertices[0] = i1 + 4 * i0;
4455 cells[k].vertices[1] = i1 + 4 * i0 + 1;
4456 cells[k].vertices[2] = i1 + 4 * i0 + 4;
4457 cells[k].vertices[3] = i1 + 4 * i0 + 5;
4459 cells[k].material_id = materials[k];
4477 const double rl2 = (right + left) / 2;
4488 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4493 for (
unsigned int i = 0; i < 4; ++i)
4495 for (
unsigned int j = 0; j < 4; ++j)
4496 cells[i].vertices[j] = cell_vertices[i][j];
4497 cells[i].material_id = 0;
4500 std::end(vertices)),
4507 cell->face(1)->set_boundary_id(1);
4509 cell->face(0)->set_boundary_id(2);
4518 const double radius_0,
4519 const double radius_1,
4520 const double half_length)
4524 vertices_tmp[0] =
Point<2>(-half_length, -radius_0);
4525 vertices_tmp[1] =
Point<2>(half_length, -radius_1);
4526 vertices_tmp[2] =
Point<2>(-half_length, radius_0);
4527 vertices_tmp[3] =
Point<2>(half_length, radius_1);
4529 const std::vector<Point<2>> vertices(std::begin(vertices_tmp),
4530 std::end(vertices_tmp));
4534 cell_vertices[0][i] = i;
4539 cells[0].vertices[i] = cell_vertices[0][i];
4541 cells[0].material_id = 0;
4546 cell->face(0)->set_boundary_id(1);
4547 cell->face(1)->set_boundary_id(2);
4549 for (
unsigned int i = 2; i < 4; ++i)
4550 cell->face(i)->set_boundary_id(0);
4567 Point<2>((a + b) / 2, (a + b) / 2),
4571 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4575 for (
unsigned int i = 0; i < 3; ++i)
4577 for (
unsigned int j = 0; j < 4; ++j)
4578 cells[i].vertices[j] = cell_vertices[i][j];
4579 cells[i].material_id = 0;
4583 std::end(vertices)),
4591 cell->face(0)->set_boundary_id(0);
4592 cell->face(2)->set_boundary_id(1);
4595 cell->face(1)->set_boundary_id(2);
4596 cell->face(2)->set_boundary_id(1);
4597 cell->face(3)->set_boundary_id(3);
4600 cell->face(0)->set_boundary_id(0);
4601 cell->face(1)->set_boundary_id(4);
4602 cell->face(3)->set_boundary_id(5);
4608 template <
int dim,
int spacedim>
4611 const std::vector<unsigned int> &repetitions,
4614 const std::vector<int> &n_cells_to_remove)
4620 for (
unsigned int d = 0;
d < dim; ++
d)
4622 Assert(std::fabs(n_cells_to_remove[d]) <= repetitions[d],
4623 ExcMessage(
"Attempting to cut away too many cells."));
4633 std::array<double, dim> h;
4635 for (
unsigned int d = 0;
d < dim; ++
d)
4638 h[
d] = (top_right[
d] - bottom_left[
d]) / repetitions[d];
4640 if (n_cells_to_remove[d] >= 0)
4644 h[
d] * std::fabs(n_cells_to_remove[d]) + bottom_left[
d];
4649 cut_step[
d] = top_right[
d] - h[
d] * std::fabs(n_cells_to_remove[d]);
4655 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4657 for (
const auto &cell : rectangle.active_cell_iterators())
4659 bool remove_cell =
true;
4660 for (
unsigned int d = 0;
d < dim && remove_cell; ++
d)
4661 if ((n_cells_to_remove[d] > 0 && cell->center()[d] >= cut_step[d]) ||
4662 (n_cells_to_remove[d] < 0 && cell->center()[d] <= cut_step[d]))
4663 remove_cell =
false;
4665 cells_to_remove.insert(cell);
4675 template <
int dim,
int spacedim>
4680 const bool internal_manifolds)
4682 if constexpr (dim == 2)
4684 const auto embed_point = [](
const double x,
4686 if constexpr (spacedim == 2)
4688 else if constexpr (spacedim == 3)
4697 const double a = 1. / (1 +
std::sqrt(2.0));
4710 for (
unsigned int i = 0; i < 5; ++i)
4712 for (
unsigned int j = 0; j < 4; ++j)
4713 cells[i].vertices[j] = circle_cell_vertices[i][j];
4714 cells[i].material_id = 0;
4719 std::begin(vertices), std::end(vertices)),
4724 if (internal_manifolds)
4729 else if constexpr (dim == 3)
4735 const unsigned int n_vertices = 16;
4736 const Point<3> vertices[n_vertices] = {
4761 const unsigned int n_cells = 7;
4762 const int cell_vertices[
n_cells][8] = {
4763 {0, 1, 4, 5, 3, 2, 7, 6},
4764 {8, 9, 12, 13, 0, 1, 4, 5},
4765 {9, 13, 1, 5, 10, 14, 2, 6},
4766 {11, 10, 3, 2, 15, 14, 7, 6},
4767 {8, 0, 12, 4, 11, 3, 15, 7},
4768 {8, 9, 0, 1, 11, 10, 3, 2},
4769 {12, 4, 13, 5, 15, 7, 14, 6}};
4771 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
4773 for (
unsigned int i = 0; i <
n_cells; ++i)
4776 cells[i].vertices[j] = cell_vertices[i][j];
4777 cells[i].material_id = 0;
4782 std::end(vertices)),
4787 if (internal_manifolds)
4798 template <
int spacedim>
4802 const double inner_radius,
4804 const unsigned int n_cells,
4821 const unsigned int N =
4822 (
n_cells == 0 ?
static_cast<unsigned int>(
4823 std::ceil((2 * pi * (
outer_radius + inner_radius) / 2) /
4833 std::vector<Point<spacedim>> vertices(2 * N);
4834 for (
unsigned int i = 0; i <
N; ++i)
4841 vertices[i +
N] = vertices[i] * (inner_radius /
outer_radius);
4843 vertices[i] += center;
4844 vertices[i +
N] += center;
4849 for (
unsigned int i = 0; i <
N; ++i)
4851 cells[i].vertices[0] = i;
4852 cells[i].vertices[1] = (i + 1) % N;
4853 cells[i].vertices[2] =
N + i;
4854 cells[i].vertices[3] =
N + ((i + 1) % N);
4856 cells[i].material_id = 0;
4862 colorize_hyper_shell(tria, center, inner_radius,
outer_radius);
4874 const double inner_radius,
4876 const unsigned int n_cells,
4888 const double inner_radius,
4890 const unsigned int n_cells,
4903 const double inner_radius,
4905 const unsigned int n_cells)
4908 tria, outer_center, inner_radius,
outer_radius, n_cells,
true);
4914 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4918 std::set<Point<dim> *> vertices_to_move;
4920 for (
const auto &face : tria.active_face_iterators())
4922 for (unsigned
int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4923 vertices_to_move.insert(&face->vertex(v));
4925 const auto shift = inner_center - outer_center;
4926 for (
const auto &p : vertices_to_move)
4958 const double half_length)
4969 switch (f->boundary_id())
4972 f->set_boundary_id(1);
4975 f->set_boundary_id(2);
4978 f->set_boundary_id(0);
5018 const unsigned int dim = 2;
5032 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
5036 for (
unsigned int i = 0; i < 3; ++i)
5038 for (
unsigned int j = 0; j < 4; ++j)
5039 cells[i].vertices[j] = cell_vertices[i][j];
5040 cells[i].material_id = 0;
5044 std::end(vertices)),
5055 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
5057 if (cell->face(i)->boundary_id() ==
5063 if (cell->face(i)->center()[0] < p[0] + 1.e-5 *
radius ||
5064 cell->face(i)->center()[1] < p[1] + 1.e-5 *
radius)
5066 cell->face(i)->set_boundary_id(1);
5085 const double a = 1. / (1 +
std::sqrt(2.0));
5096 const int cell_vertices[5][4] = {{0, 1, 2, 3},
5103 for (
unsigned int i = 0; i < 4; ++i)
5105 for (
unsigned int j = 0; j < 4; ++j)
5106 cells[i].vertices[j] = cell_vertices[i][j];
5107 cells[i].material_id = 0;
5111 std::end(vertices)),
5122 for (
const unsigned int i :
GeometryInfo<2>::face_indices())
5124 if (cell->face(i)->boundary_id() ==
5129 if (cell->face(i)->center()[0] < p[0] + 1.e-5 *
radius)
5131 cell->face(i)->set_boundary_id(1);
5147 const double inner_radius,
5149 const unsigned int n_cells,
5165 const unsigned int N =
5166 (
n_cells == 0 ?
static_cast<unsigned int>(
5177 std::vector<Point<2>> vertices(2 * (N + 1));
5178 for (
unsigned int i = 0; i <=
N; ++i)
5188 std::
sin(pi * i /
N - pi / 2)) *
5190 vertices[i +
N + 1] = vertices[i] * (inner_radius /
outer_radius);
5192 vertices[i] += center;
5193 vertices[i +
N + 1] += center;
5199 for (
unsigned int i = 0; i <
N; ++i)
5201 cells[i].vertices[0] = i;
5202 cells[i].vertices[1] = (i + 1) % (N + 1);
5203 cells[i].vertices[2] =
N + 1 + i;
5204 cells[i].vertices[3] =
N + 1 + ((i + 1) % (N + 1));
5206 cells[i].material_id = 0;
5214 for (; cell != tria.
end(); ++cell)
5216 cell->face(2)->set_boundary_id(1);
5218 tria.
begin()->face(0)->set_boundary_id(3);
5220 tria.
last()->face(1)->set_boundary_id(2);
5231 const double inner_radius,
5233 const unsigned int n_cells,
5249 const unsigned int N =
5250 (
n_cells == 0 ?
static_cast<unsigned int>(
5261 std::vector<Point<2>> vertices(2 * (N + 1));
5262 for (
unsigned int i = 0; i <=
N; ++i)
5272 vertices[i +
N + 1] = vertices[i] * (inner_radius /
outer_radius);
5274 vertices[i] += center;
5275 vertices[i +
N + 1] += center;
5281 for (
unsigned int i = 0; i <
N; ++i)
5283 cells[i].vertices[0] = i;
5284 cells[i].vertices[1] = (i + 1) % (N + 1);
5285 cells[i].vertices[2] =
N + 1 + i;
5286 cells[i].vertices[3] =
N + 1 + ((i + 1) % (N + 1));
5288 cells[i].material_id = 0;
5296 for (; cell != tria.
end(); ++cell)
5298 cell->face(2)->set_boundary_id(1);
5300 tria.
begin()->face(0)->set_boundary_id(3);
5302 tria.
last()->face(1)->set_boundary_id(2);
5319 const double rl2 = (right + left) / 2;
5320 const double len = (right - left) / 2.;
5333 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
5334 {9, 4, 2, 5, 19, 14, 12, 15},
5335 {3, 2, 7, 6, 13, 12, 17, 16},
5336 {2, 5, 6, 8, 12, 15, 16, 18}};
5338 for (
unsigned int i = 0; i < 4; ++i)
5340 for (
unsigned int j = 0; j < 8; ++j)
5341 cells[i].vertices[j] = cell_vertices[i][j];
5342 cells[i].material_id = 0;
5345 std::end(vertices)),
5352 cell->face(1)->set_boundary_id(1);
5354 cell->face(0)->set_boundary_id(2);
5366 const double thickness,
5370 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
5372 std::vector<Point<3>> vertices(64);
5374 coords[0] = left - thickness;
5377 coords[3] = right + thickness;
5380 for (
const double z : coords)
5383 vertices[k++] =
Point<3>(x, y, z);
5386 24, 26, 5, 4, 6, 1, 0,
5387 2, 9, 8, 10, 37, 36, 38,
5388 33, 32, 34, 41, 40, 42};
5390 std::vector<CellData<3>> cells(27);
5392 for (
unsigned int z = 0; z < 3; ++z)
5393 for (
unsigned int y = 0; y < 3; ++y)
5394 for (
unsigned int x = 0; x < 3; ++x)
5396 cells[k].vertices[0] = x + 4 * y + 16 * z;
5397 cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
5398 cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
5399 cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
5400 cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
5401 cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
5402 cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
5403 cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
5405 cells[k].material_id = materials[k];
5418 const double radius_0,
5419 const double radius_1,
5420 const double half_length)
5423 ExcMessage(
"The output triangulation object needs to be empty."));
5428 const auto n_slices = 1 +
static_cast<unsigned int>(std::ceil(
5429 half_length /
std::max(radius_0, radius_1)));
5442 auto shift_radii = [=](
const Point<3> &p) {
5443 const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
5444 const double factor = slope * (p[0] - -half_length) + 1.0;
5445 return Point<3>(p[0], factor * p[1], factor * p[2]);
5451 for (
const auto &face : triangulation.active_face_iterators())
5452 if (face->at_boundary())
5454 if (
std::abs(face->center()[0] - -half_length) < 1e-8 * half_length)
5455 face->set_boundary_id(1);
5456 else if (
std::abs(face->center()[0] - half_length) <
5458 face->set_boundary_id(2);
5460 face->set_all_manifold_ids(0);
5483 Point<3>((a + b) / 2, a, (a + b) / 2),
5490 Point<3>((a + b) / 2, (a + b) / 2, a),
5492 Point<3>(a, (a + b) / 2, (a + b) / 2),
5493 Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
5494 Point<3>(b, (a + b) / 2, (a + b) / 2),
5496 Point<3>((a + b) / 2, (a + b) / 2, b),
5504 Point<3>((a + b) / 2, b, (a + b) / 2),
5508 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5509 {1, 2, 10, 11, 4, 5, 13, 14},
5510 {3, 4, 12, 13, 6, 7, 15, 16},
5511 {4, 5, 13, 14, 7, 8, 16, 17},
5512 {9, 10, 18, 19, 12, 13, 21, 22},
5513 {10, 11, 19, 20, 13, 14, 22, 23},
5514 {12, 13, 21, 22, 15, 16, 24, 25}};
5518 for (
unsigned int i = 0; i < 7; ++i)
5520 for (
unsigned int j = 0; j < 8; ++j)
5521 cells[i].vertices[j] = cell_vertices[i][j];
5522 cells[i].material_id = 0;
5526 std::end(vertices)),
5540 const unsigned int n_rotate_middle_square)
5543 ExcMessage(
"The number of rotation by pi/2 of the right square "
5544 "must be in the half-open range [0,4)."));
5546 constexpr unsigned int dim = 2;
5548 const unsigned int n_cells = 5;
5549 std::vector<CellData<dim>> cells(n_cells);
5552 const std::vector<Point<dim>> vertices = {
Point<dim>(0, 0),
5567 unsigned int cell_vertices[
n_cells][4] = {{0, 1, 2, 3},
5573 switch (n_rotate_middle_square)
5577 cell_vertices[1][0] = 4;
5578 cell_vertices[1][1] = 5;
5579 cell_vertices[1][2] = 1;
5580 cell_vertices[1][3] = 3;
5586 cell_vertices[1][0] = 5;
5587 cell_vertices[1][1] = 3;
5588 cell_vertices[1][2] = 4;
5589 cell_vertices[1][3] = 1;
5595 cell_vertices[1][0] = 3;
5596 cell_vertices[1][1] = 1;
5597 cell_vertices[1][2] = 5;
5598 cell_vertices[1][3] = 4;
5610 for (
const unsigned int vertex_index :
5625 const bool face_orientation,
5626 const bool face_flip,
5627 const bool face_rotation,
5628 const bool manipulate_left_cube)
5630 constexpr unsigned int dim = 3;
5632 const unsigned int n_cells = 2;
5633 std::vector<CellData<dim>> cells(n_cells);
5636 const std::vector<Point<dim>> vertices = {
Point<dim>(0, 0, 0),
5649 unsigned int cell_vertices[
n_cells][8] = {
5650 {0, 1, 2, 3, 4, 5, 6, 7},
5651 {1, 8, 3, 9, 5, 10, 7, 11}};
5654 const unsigned int this_case = 4 *
static_cast<int>(face_orientation) +
5655 2 *
static_cast<int>(face_flip) +
5656 static_cast<int>(face_rotation);
5658 if (manipulate_left_cube)
5664 cell_vertices[0][0] = 1;
5665 cell_vertices[0][1] = 0;
5666 cell_vertices[0][2] = 5;
5667 cell_vertices[0][3] = 4;
5668 cell_vertices[0][4] = 3;
5669 cell_vertices[0][5] = 2;
5670 cell_vertices[0][6] = 7;
5671 cell_vertices[0][7] = 6;
5677 cell_vertices[0][0] = 5;
5678 cell_vertices[0][1] = 4;
5679 cell_vertices[0][2] = 7;
5680 cell_vertices[0][3] = 6;
5681 cell_vertices[0][4] = 1;
5682 cell_vertices[0][5] = 0;
5683 cell_vertices[0][6] = 3;
5684 cell_vertices[0][7] = 2;
5690 cell_vertices[0][0] = 7;
5691 cell_vertices[0][1] = 6;
5692 cell_vertices[0][2] = 3;
5693 cell_vertices[0][3] = 2;
5694 cell_vertices[0][4] = 5;
5695 cell_vertices[0][5] = 4;
5696 cell_vertices[0][6] = 1;
5697 cell_vertices[0][7] = 0;
5702 cell_vertices[0][0] = 3;
5703 cell_vertices[0][1] = 2;
5704 cell_vertices[0][2] = 1;
5705 cell_vertices[0][3] = 0;
5706 cell_vertices[0][4] = 7;
5707 cell_vertices[0][5] = 6;
5708 cell_vertices[0][6] = 5;
5709 cell_vertices[0][7] = 4;
5715 cell_vertices[0][0] = 0;
5716 cell_vertices[0][1] = 1;
5717 cell_vertices[0][2] = 2;
5718 cell_vertices[0][3] = 3;
5719 cell_vertices[0][4] = 4;
5720 cell_vertices[0][5] = 5;
5721 cell_vertices[0][6] = 6;
5722 cell_vertices[0][7] = 7;
5728 cell_vertices[0][0] = 2;
5729 cell_vertices[0][1] = 3;
5730 cell_vertices[0][2] = 6;
5731 cell_vertices[0][3] = 7;
5732 cell_vertices[0][4] = 0;
5733 cell_vertices[0][5] = 1;
5734 cell_vertices[0][6] = 4;
5735 cell_vertices[0][7] = 5;
5741 cell_vertices[0][0] = 6;
5742 cell_vertices[0][1] = 7;
5743 cell_vertices[0][2] = 4;
5744 cell_vertices[0][3] = 5;
5745 cell_vertices[0][4] = 2;
5746 cell_vertices[0][5] = 3;
5747 cell_vertices[0][6] = 0;
5748 cell_vertices[0][7] = 1;
5754 cell_vertices[0][0] = 4;
5755 cell_vertices[0][1] = 5;
5756 cell_vertices[0][2] = 0;
5757 cell_vertices[0][3] = 1;
5758 cell_vertices[0][4] = 6;
5759 cell_vertices[0][5] = 7;
5760 cell_vertices[0][6] = 2;
5761 cell_vertices[0][7] = 3;
5772 cell_vertices[1][0] = 8;
5773 cell_vertices[1][1] = 1;
5774 cell_vertices[1][2] = 10;
5775 cell_vertices[1][3] = 5;
5776 cell_vertices[1][4] = 9;
5777 cell_vertices[1][5] = 3;
5778 cell_vertices[1][6] = 11;
5779 cell_vertices[1][7] = 7;
5785 cell_vertices[1][0] = 10;
5786 cell_vertices[1][1] = 5;
5787 cell_vertices[1][2] = 11;
5788 cell_vertices[1][3] = 7;
5789 cell_vertices[1][4] = 8;
5790 cell_vertices[1][5] = 1;
5791 cell_vertices[1][6] = 9;
5792 cell_vertices[1][7] = 3;
5798 cell_vertices[1][0] = 11;
5799 cell_vertices[1][1] = 7;
5800 cell_vertices[1][2] = 9;
5801 cell_vertices[1][3] = 3;
5802 cell_vertices[1][4] = 10;
5803 cell_vertices[1][5] = 5;
5804 cell_vertices[1][6] = 8;
5805 cell_vertices[1][7] = 1;
5811 cell_vertices[1][0] = 9;
5812 cell_vertices[1][1] = 3;
5813 cell_vertices[1][2] = 8;
5814 cell_vertices[1][3] = 1;
5815 cell_vertices[1][4] = 11;
5816 cell_vertices[1][5] = 7;
5817 cell_vertices[1][6] = 10;
5818 cell_vertices[1][7] = 5;
5824 cell_vertices[1][0] = 1;
5825 cell_vertices[1][1] = 8;
5826 cell_vertices[1][2] = 3;
5827 cell_vertices[1][3] = 9;
5828 cell_vertices[1][4] = 5;
5829 cell_vertices[1][5] = 10;
5830 cell_vertices[1][6] = 7;
5831 cell_vertices[1][7] = 11;
5837 cell_vertices[1][0] = 5;
5838 cell_vertices[1][1] = 10;
5839 cell_vertices[1][2] = 1;
5840 cell_vertices[1][3] = 8;
5841 cell_vertices[1][4] = 7;
5842 cell_vertices[1][5] = 11;
5843 cell_vertices[1][6] = 3;
5844 cell_vertices[1][7] = 9;
5850 cell_vertices[1][0] = 7;
5851 cell_vertices[1][1] = 11;
5852 cell_vertices[1][2] = 5;
5853 cell_vertices[1][3] = 10;
5854 cell_vertices[1][4] = 3;
5855 cell_vertices[1][5] = 9;
5856 cell_vertices[1][6] = 1;
5857 cell_vertices[1][7] = 8;
5863 cell_vertices[1][0] = 3;
5864 cell_vertices[1][1] = 9;
5865 cell_vertices[1][2] = 7;
5866 cell_vertices[1][3] = 11;
5867 cell_vertices[1][4] = 1;
5868 cell_vertices[1][5] = 8;
5869 cell_vertices[1][6] = 5;
5870 cell_vertices[1][7] = 10;
5880 for (
const unsigned int vertex_index :
5894 template <
int spacedim>
5902 const std::set<types::boundary_id> boundary_ids = {0};
5914 const unsigned int x_subdivisions,
5916 const double half_length)
5923 std::vector<Point<3>> vertices;
5924 const double initial_height = -half_length;
5925 const double height_increment = 2. * half_length / x_subdivisions;
5927 for (
unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep)
5929 const double height = initial_height + height_increment * rep;
5931 vertices.emplace_back(-d,
height, -d);
5932 vertices.emplace_back(d,
height, -d);
5933 vertices.emplace_back(-a,
height, -a);
5934 vertices.emplace_back(a,
height, -a);
5935 vertices.emplace_back(-a,
height, a);
5936 vertices.emplace_back(a,
height, a);
5937 vertices.emplace_back(-d,
height, d);
5938 vertices.emplace_back(d,
height, d);
5942 for (
auto &vertex : vertices)
5944 const double h = vertex[1];
5945 vertex[1] = -vertex[0];
5949 std::vector<std::vector<int>> cell_vertices;
5950 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5951 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5952 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5953 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5954 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5956 for (
unsigned int rep = 1; rep < x_subdivisions; ++rep)
5958 for (
unsigned int i = 0; i < 5; ++i)
5960 std::vector<int> new_cell_vertices(8);
5961 for (
unsigned int j = 0; j < 8; ++j)
5962 new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep;
5963 cell_vertices.push_back(new_cell_vertices);
5967 unsigned int n_cells = x_subdivisions * 5;
5969 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5971 for (
unsigned int i = 0; i <
n_cells; ++i)
5973 for (
unsigned int j = 0; j < 8; ++j)
5974 cells[i].vertices[j] = cell_vertices[i][j];
5975 cells[i].material_id = 0;
5979 std::end(vertices)),
5999 for (
const auto &cell : tria.cell_iterators())
6001 if (cell->at_boundary(i))
6003 if (cell->face(i)->center()[0] > half_length - tolerance)
6005 cell->face(i)->set_boundary_id(2);
6008 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
6010 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
6011 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
6012 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
6013 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
6015 cell->face(i)->line(e)->set_boundary_id(2);
6016 cell->face(i)->line(e)->set_manifold_id(
6020 else if (cell->face(i)->center()[0] < -half_length + tolerance)
6022 cell->face(i)->set_boundary_id(1);
6025 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
6027 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
6028 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
6029 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
6030 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
6032 cell->face(i)->line(e)->set_boundary_id(1);
6033 cell->face(i)->line(e)->set_manifold_id(
6046 const double half_length)
6057 const unsigned int dim = 3;
6065 const double a = 0.528;
6066 const double b = 0.4533;
6067 const double c = 0.3752;
6084 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
6085 {1, 6, 2, 4, 8, 13, 9, 11},
6086 {5, 3, 6, 4, 12, 10, 13, 11},
6087 {7, 9, 10, 11, 14, 8, 12, 13}};
6091 for (
unsigned int i = 0; i < 4; ++i)
6093 for (
unsigned int j = 0; j < 8; ++j)
6094 cells[i].vertices[j] = cell_vertices[i][j];
6095 cells[i].material_id = 0;
6099 std::end(vertices)),
6109 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
6111 if (cell->face(i)->boundary_id() ==
6116 if (cell->face(i)->center()[0] < center[0] + 1.e-5 *
radius ||
6117 cell->face(i)->center()[1] < center[1] + 1.e-5 *
radius ||
6118 cell->face(i)->center()[2] < center[2] + 1.e-5 *
radius)
6120 cell->face(i)->set_boundary_id(1);
6124 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
6127 const Point<3> line_vertices[2] = {
6128 cell->face(i)->line(j)->vertex(0),
6129 cell->face(i)->line(j)->vertex(1)};
6130 if ((std::fabs(line_vertices[0].distance(center) -
radius) >
6132 (std::fabs(line_vertices[1].distance(center) -
radius) >
6135 cell->face(i)->line(j)->set_boundary_id(1);
6136 cell->face(i)->line(j)->set_manifold_id(
6160 const double b = a / 2.0;
6161 const double c =
d / 2.0;
6186 int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
6187 {0, 2, 8, 10, 6, 4, 14, 12},
6188 {2, 3, 10, 11, 4, 5, 12, 13},
6189 {1, 7, 9, 15, 3, 5, 11, 13},
6190 {6, 4, 14, 12, 7, 5, 15, 13},
6191 {8, 10, 9, 11, 14, 12, 15, 13}};
6195 for (
unsigned int i = 0; i < 6; ++i)
6197 for (
unsigned int j = 0; j < 8; ++j)
6198 cells[i].vertices[j] = cell_vertices[i][j];
6199 cells[i].material_id = 0;
6203 std::end(vertices)),
6219 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6221 if (!cell->at_boundary(i))
6227 if (cell->face(i)->center()[0] < center[0] + 1.e-5 *
radius)
6229 cell->face(i)->set_boundary_id(1);
6231 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
6234 const Point<3> line_vertices[2] = {
6235 cell->face(i)->line(j)->vertex(0),
6236 cell->face(i)->line(j)->vertex(1)};
6237 if ((std::fabs(line_vertices[0].distance(center) -
radius) >
6239 (std::fabs(line_vertices[1].distance(center) -
radius) >
6242 cell->face(i)->line(j)->set_boundary_id(1);
6243 cell->face(i)->line(j)->set_manifold_id(
6273 for (
unsigned int round = 0; round < dim; ++round)
6278 std::vector<Point<dim>> new_points(tria_copy.
n_vertices());
6280 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
6288 else if (round == 1)
6290 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
6299 else if (round == 2)
6300 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
6313 std::vector<CellData<dim>> cells;
6314 cells.reserve(tria_copy.
n_cells());
6315 for (
const auto &cell : tria_copy.cell_iterators())
6319 data.vertices[v] = cell->vertex_index(v);
6320 data.material_id = cell->material_id();
6321 data.manifold_id = cell->manifold_id();
6322 cells.push_back(
data);
6330 if (round == dim - 1)
6339 for (
const auto &cell : tria.cell_iterators())
6340 if (cell->center().norm_square() > 0.4 *
radius)
6341 cell->set_manifold_id(1);
6362 const double inner_radius,
6365 std::vector<Point<3>> vertices;
6366 std::vector<CellData<3>> cells;
6368 const double irad = inner_radius /
std::sqrt(3.0);
6372 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
6382 vertices.reserve(8);
6383 for (
unsigned int i = 0; i < 8; ++i)
6384 vertices.push_back(p + hexahedron[i] * irad);
6385 for (
unsigned int i = 0; i < 8; ++i)
6386 vertices.push_back(p + hexahedron[i] * orad);
6388 const unsigned int n_cells = 6;
6389 const int cell_vertices[
n_cells][8] = {
6390 {8, 9, 10, 11, 0, 1, 2, 3},
6391 {9, 11, 1, 3, 13, 15, 5, 7},
6392 {12, 13, 4, 5, 14, 15, 6, 7},
6393 {8, 0, 10, 2, 12, 4, 14, 6},
6394 {8, 9, 0, 1, 12, 13, 4, 5},
6395 {10, 2, 11, 3, 14, 6, 15, 7}};
6399 for (
unsigned int i = 0; i <
n_cells; ++i)
6402 cells[i].vertices[j] = cell_vertices[i][j];
6403 cells[i].material_id = 0;
6414 const double inner_radius,
6417 std::vector<Point<3>> vertices;
6418 std::vector<CellData<3>> cells;
6420 const double irad = inner_radius /
std::sqrt(3.0);
6427 static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0},
6435 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
6444 vertices.reserve(8);
6445 for (
unsigned int i = 0; i < 8; ++i)
6446 vertices.push_back(p + hexahedron[i] * irad);
6447 for (
unsigned int i = 0; i < 6; ++i)
6448 vertices.push_back(p + octahedron[i] * inner_radius);
6449 for (
unsigned int i = 0; i < 8; ++i)
6450 vertices.push_back(p + hexahedron[i] * orad);
6451 for (
unsigned int i = 0; i < 6; ++i)
6454 const unsigned int n_cells = 12;
6455 const unsigned int rhombi[
n_cells][4] = {{10, 4, 0, 8},
6470 for (
unsigned int i = 0; i <
n_cells; ++i)
6472 for (
unsigned int j = 0; j < 4; ++j)
6474 cells[i].vertices[j] = rhombi[i][j];
6475 cells[i].vertices[j + 4] = rhombi[i][j] + 14;
6477 cells[i].material_id = 0;
6487 const unsigned int n,
6488 const unsigned int n_refinement_steps,
6490 const double inner_radius,
6504 const unsigned int outer_radius_factor = 1 << n_refinement_steps;
6510 (outer_radius_factor - 1) * inner_radius);
6516 (outer_radius_factor - 1) * inner_radius);
6519 for (
unsigned int r = 0; r < n_refinement_steps; ++r)
6522 std::set<Triangulation<3>::active_cell_iterator> cells_to_remove;
6526 for (
const auto &cell : tmp.active_cell_iterators())
6528 unsigned int n_vertices_inside = 0;
6530 if ((cell->vertex(v) - p).norm_square() <
6531 inner_radius * inner_radius * (1 + 1
e-12))
6532 ++n_vertices_inside;
6533 if (n_vertices_inside < 4)
6534 cells_to_remove.insert(cell);
6538 if (r == n_refinement_steps - 1)
6548 tmp = std::move(copy);
6566 const double inner_radius,
6568 const unsigned int n_cells,
6574 unsigned int n_refinement_steps = 0;
6575 unsigned int n_cells_coarsened =
n_cells;
6576 if (n_cells != 96 && n_cells > 12)
6577 while (n_cells_coarsened > 12 && n_cells_coarsened % 4 == 0)
6579 ++n_refinement_steps;
6580 n_cells_coarsened /= 4;
6582 Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 ||
6583 (n_refinement_steps > 0 &&
6584 (n_cells_coarsened == 6 || n_cells_coarsened == 12)),
6585 ExcMessage(
"Invalid number of coarse mesh cells"));
6587 const unsigned int n = n_refinement_steps > 0 ?
6588 4 * n_cells_coarsened :
6594 internal::hyper_shell_6(tria, p, inner_radius,
outer_radius);
6597 internal::hyper_shell_12(tria, p, inner_radius,
outer_radius);
6601 internal::hyper_shell_24_48(
6602 tria, n, n_refinement_steps, p, inner_radius,
outer_radius);
6611 internal::hyper_shell_12(tmp, p, inner_radius,
outer_radius);
6628 colorize_hyper_shell(tria, p, inner_radius,
outer_radius);
6638 const double inner_radius,
6640 const unsigned int ,
6648 const double a = inner_radius /
std::sqrt(2.0);
6650 const double b = a / 2.0;
6651 const double c =
d / 2.0;
6653 const double hb = inner_radius *
std::sqrt(3.0) / 2.0;
6676 int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
6677 {0, 2, 8, 10, 6, 4, 14, 12},
6678 {1, 7, 9, 15, 3, 5, 11, 13},
6679 {6, 4, 14, 12, 7, 5, 15, 13},
6680 {8, 10, 9, 11, 14, 12, 15, 13}};
6684 for (
unsigned int i = 0; i < 5; ++i)
6686 for (
unsigned int j = 0; j < 8; ++j)
6687 cells[i].vertices[j] = cell_vertices[i][j];
6688 cells[i].material_id = 0;
6692 std::end(vertices)),
6704 for (
const auto &cell : tria.cell_iterators())
6706 if (cell->at_boundary(f))
6707 cell->face(f)->set_all_boundary_ids(2);
6713 for (
const auto &cell : tria.cell_iterators())
6715 if (cell->at_boundary(f))
6718 const Point<3> face_center(face->center(
true));
6720 if (
std::abs(face_center[0] - center[0]) >
6721 1.e-6 * face_center.
norm())
6723 if (
std::abs((face_center - center).
norm() - inner_radius) <
6725 face->set_all_boundary_ids(0);
6727 face->set_all_boundary_ids(1);
6739 const double inner_radius,
6741 const unsigned int n,
6746 if (n == 0 || n == 3)
6748 const double a = inner_radius *
std::sqrt(2.0) / 2e0;
6750 const double c = a *
std::sqrt(3.0) / 2e0;
6753 const double h = inner_radius / 2e0;
6755 std::vector<Point<3>> vertices;
6757 vertices.push_back(center +
Point<3>(0, inner_radius, 0));
6758 vertices.push_back(center +
Point<3>(a, a, 0));
6759 vertices.push_back(center +
Point<3>(b, b, 0));
6761 vertices.push_back(center +
Point<3>(0, a, a));
6762 vertices.push_back(center +
Point<3>(c, c, h));
6763 vertices.push_back(center +
Point<3>(d, d, e));
6764 vertices.push_back(center +
Point<3>(0, b, b));
6765 vertices.push_back(center +
Point<3>(inner_radius, 0, 0));
6767 vertices.push_back(center +
Point<3>(a, 0, a));
6768 vertices.push_back(center +
Point<3>(b, 0, b));
6769 vertices.push_back(center +
Point<3>(0, 0, inner_radius));
6772 const int cell_vertices[3][8] = {
6773 {0, 1, 3, 2, 4, 5, 7, 6},
6774 {1, 8, 2, 9, 5, 10, 6, 11},
6775 {4, 5, 7, 6, 12, 10, 13, 11},
6777 std::vector<CellData<3>> cells(3);
6779 for (
unsigned int i = 0; i < 3; ++i)
6781 for (
unsigned int j = 0; j < 8; ++j)
6782 cells[i].vertices[j] = cell_vertices[i][j];
6783 cells[i].material_id = 0;
6799 colorize_quarter_hyper_shell(tria, center, inner_radius,
outer_radius);
6807 const double length,
6808 const double inner_radius,
6810 const unsigned int n_radial_cells,
6811 const unsigned int n_axial_cells,
6828 const unsigned int N_r =
6829 (n_radial_cells == 0 ?
static_cast<unsigned int>(std::ceil(
6833 const unsigned int N_z =
6834 (n_axial_cells == 0 ?
6835 static_cast<unsigned int>(std::ceil(
6836 length / (2 * pi * (
outer_radius + inner_radius) / 2 / N_r))) :
6845 std::vector<Point<2>> vertices_2d(2 * N_r);
6846 for (
unsigned int i = 0; i < N_r; ++i)
6851 vertices_2d[i + N_r] = vertices_2d[i] * (inner_radius /
outer_radius);
6854 std::vector<Point<3>> vertices_3d;
6855 vertices_3d.reserve(2 * N_r * (N_z + 1));
6856 for (
unsigned int j = 0; j <= N_z; ++j)
6857 for (
unsigned int i = 0; i < 2 * N_r; ++i)
6859 const Point<3> v(vertices_2d[i][0],
6862 vertices_3d.push_back(v);
6865 std::vector<CellData<3>> cells(N_r * N_z,
CellData<3>());
6867 for (
unsigned int j = 0; j < N_z; ++j)
6868 for (
unsigned int i = 0; i < N_r; ++i)
6870 cells[i + j * N_r].vertices[0] = i + (j + 1) * 2 * N_r;
6871 cells[i + j * N_r].vertices[1] = (i + 1) % N_r + (j + 1) * 2 * N_r;
6872 cells[i + j * N_r].vertices[2] = i + j * 2 * N_r;
6873 cells[i + j * N_r].vertices[3] = (i + 1) % N_r + j * 2 * N_r;
6875 cells[i + j * N_r].vertices[4] = N_r + i + (j + 1) * 2 * N_r;
6876 cells[i + j * N_r].vertices[5] =
6877 N_r + ((i + 1) % N_r) + (j + 1) * 2 * N_r;
6878 cells[i + j * N_r].vertices[6] = N_r + i + j * 2 * N_r;
6879 cells[i + j * N_r].vertices[7] = N_r + ((i + 1) % N_r) + j * 2 * N_r;
6881 cells[i + j * N_r].material_id = 0;
6900 double eps_z = 1
e-6 * length;
6905 double face_inner_radius = std::numeric_limits<double>::max();
6906 double face_outer_radius = 0.;
6911 for (
const auto &cell : tria.active_cell_iterators())
6914 if (!cell->face(f)->at_boundary())
6917 const auto face_center = cell->face(f)->center();
6918 const double z = face_center[2];
6920 if ((std::fabs(z) > eps_z) &&
6921 (std::fabs(z - length) > eps_z))
6924 face_center[1] * face_center[1]);
6930 double mid_radial_distance = 0.5 * (face_outer_radius - face_inner_radius);
6932 for (
const auto &cell : tria.active_cell_iterators())
6935 if (cell->face(f)->at_boundary())
6937 const auto face_center = cell->face(f)->center();
6940 face_center[1] * face_center[1]);
6942 const double z = face_center[2];
6944 if (std::fabs(z) < eps_z)
6946 cell->face(f)->set_boundary_id(2);
6948 else if (std::fabs(z - length) <
6951 cell->face(f)->set_boundary_id(3);
6953 else if (std::fabs(
radius - face_inner_radius) >
6954 mid_radial_distance)
6956 cell->face(f)->set_boundary_id(1);
6958 else if (std::fabs(
radius - face_inner_radius) <
6959 mid_radial_distance)
6961 cell->face(f)->set_boundary_id(0);
6970 template <
int dim,
int spacedim>
6975 const double duplicated_vertex_tolerance,
6976 const bool copy_manifold_ids,
6977 const bool copy_boundary_ids)
6979 std::vector<Point<spacedim>> vertices;
6980 std::vector<CellData<dim>> cells;
6983 unsigned int n_accumulated_vertices = 0;
6984 for (
const auto triangulation : triangulations)
6987 ExcMessage(
"The input triangulations must be non-empty "
6988 "and must not be refined."));
6990 auto [tria_vertices, tria_cells, tria_subcell_data] =
7000 vertices.insert(vertices.end(),
7001 tria_vertices.begin(),
7002 tria_vertices.end());
7005 for (
unsigned int &vertex_n : cell_data.vertices)
7006 vertex_n += n_accumulated_vertices;
7007 cells.push_back(cell_data);
7011 if (copy_manifold_ids)
7013 for (
CellData<1> &line_data : tria_subcell_data.boundary_lines)
7017 for (
unsigned int &vertex_n : line_data.vertices)
7018 vertex_n += n_accumulated_vertices;
7019 line_data.boundary_id =
7024 for (
CellData<2> &quad_data : tria_subcell_data.boundary_quads)
7028 for (
unsigned int &vertex_n : quad_data.vertices)
7029 vertex_n += n_accumulated_vertices;
7030 quad_data.boundary_id =
7036 n_accumulated_vertices += triangulation->
n_vertices();
7040 std::vector<unsigned int> considered_vertices;
7044 considered_vertices,
7045 duplicated_vertex_tolerance);
7049 if (std::all_of(cells.begin(), cells.end(), [](
const auto &cell) {
7050 return cell.vertices.size() ==
7051 ReferenceCells::get_hypercube<dim>().n_vertices();
7057 if (!copy_manifold_ids)
7060 if (copy_boundary_ids)
7062 auto result_cell = result.
begin();
7063 for (
const auto &tria : triangulations)
7065 for (
const auto &cell : tria->cell_iterators())
7067 for (
const auto &f : cell->face_indices())
7068 if (result_cell->face(f)->at_boundary())
7069 result_cell->face(f)->set_boundary_id(
7076 Assert(duplicated_vertex_tolerance > 0.0 ||
7077 n_accumulated_vertices == result.
n_vertices(),
7083 template <
int dim,
int spacedim>
7088 const double duplicated_vertex_tolerance,
7089 const bool copy_manifold_ids,
7090 const bool copy_boundary_ids)
7093 if (triangulation_1.
n_cells() == 0)
7095 if (&result != &triangulation_2)
7098 else if (triangulation_2.
n_cells() == 0)
7100 if (&result != &triangulation_1)
7106 duplicated_vertex_tolerance,
7136 template <
int structdim>
7140 static_assert(structdim == 1 || structdim == 2,
7141 "This function is only implemented for lines and "
7149 std::sort(std::begin(cell_data.vertices),
7150 std::end(cell_data.vertices));
7151 else if (structdim == 2)
7154 std::array<unsigned int, 4> renumbering{};
7155 std::copy(std::begin(cell_data.vertices),
7156 std::end(cell_data.vertices),
7157 renumbering.begin());
7168 std::swap(renumbering[2], renumbering[3]);
7169 std::rotate(renumbering.begin(),
7170 std::min_element(renumbering.begin(),
7174 std::swap(renumbering[2], renumbering[3]);
7182 if (renumbering[1] > renumbering[2])
7183 std::swap(renumbering[1], renumbering[2]);
7184 std::copy(renumbering.begin(),
7186 std::begin(cell_data.vertices));
7193 return std::lexicographical_compare(std::begin(a.
vertices),
7195 std::begin(
b.vertices),
7196 std::end(
b.vertices));
7198 std::sort(subcell_data.begin(), subcell_data.end(), compare);
7203 auto left = subcell_data.begin();
7204 while (left != subcell_data.end())
7207 std::upper_bound(left, subcell_data.end(), *left, compare);
7210 if (left + 1 != right)
7211 for (
auto it = left; it != right; ++it)
7214 Assert(it->manifold_id == left->manifold_id,
7216 "In the process of grid generation a single "
7217 "line or quadrilateral has been assigned two "
7218 "different manifold ids. This can happen when "
7219 "a Triangulation is copied, e.g., via "
7220 "GridGenerator::replicate_triangulation() and "
7221 "not all external boundary faces have the same "
7222 "manifold id. Double check that all faces "
7223 "which you expect to be merged together have "
7224 "the same manifold id."));
7229 subcell_data.erase(std::unique(subcell_data.begin(), subcell_data.end()),
7230 subcell_data.end());
7236 template <
int dim,
int spacedim>
7239 const std::vector<unsigned int> &extents,
7245 for (
const auto &extent : extents)
7248 "The Triangulation must be copied at least one time in "
7249 "each coordinate dimension."));
7252 const auto &
min = bbox.get_boundary_points().first;
7253 const auto &
max = bbox.get_boundary_points().second;
7255 std::array<Tensor<1, spacedim>, dim> offsets;
7256 for (
unsigned int d = 0;
d < dim; ++
d)
7257 offsets[d][d] = max[d] - min[d];
7261 for (
unsigned int d = 0;
d < dim; ++
d)
7263 auto [input_vertices, input_cell_data, input_subcell_data] =
7266 std::vector<Point<spacedim>> output_vertices = input_vertices;
7267 std::vector<CellData<dim>> output_cell_data = input_cell_data;
7268 SubCellData output_subcell_data = input_subcell_data;
7270 for (
unsigned int k = 1; k < extents[
d]; ++k)
7272 const std::size_t vertex_offset = k * input_vertices.size();
7279 output_cell_data.push_back(cell_data);
7280 for (
unsigned int &vertex : output_cell_data.back().vertices)
7281 vertex += vertex_offset;
7285 input_subcell_data.boundary_lines)
7288 for (
unsigned int &vertex :
7289 output_subcell_data.boundary_lines.back().vertices)
7290 vertex += vertex_offset;
7293 input_subcell_data.boundary_quads)
7296 for (
unsigned int &vertex :
7297 output_subcell_data.boundary_quads.back().vertices)
7298 vertex += vertex_offset;
7303 std::vector<unsigned int> boundary_vertices;
7307 output_subcell_data,
7322 output_subcell_data.boundary_lines)
7328 tria_to_replicate.
clear();
7331 output_subcell_data);
7339 template <
int dim,
int spacedim>
7347 ExcMessage(
"The two input triangulations are not derived from "
7348 "the same coarse mesh as required."));
7351 &triangulation_1) ==
nullptr) &&
7354 &triangulation_2) ==
nullptr),
7355 ExcMessage(
"The source triangulations for this function must both "
7356 "be available entirely locally, and not be distributed "
7357 "triangulations."));
7374 for (
unsigned int iteration = 0; iteration < triangulation_2.
n_levels();
7380 bool any_cell_flagged =
false;
7381 for (
const auto &result_cell : result.active_cell_iterators())
7382 if (intergrid_map[result_cell]->has_children())
7384 any_cell_flagged =
true;
7385 result_cell->set_refine_flag();
7388 if (any_cell_flagged ==
false)
7397 template <
int dim,
int spacedim>
7407 std::vector<Point<spacedim>> vertices = input_triangulation.
get_vertices();
7411 std::vector<CellData<dim>> cells;
7412 for (
const auto &cell : input_triangulation.active_cell_iterators())
7413 if (cells_to_remove.find(cell) == cells_to_remove.
end())
7415 Assert(
static_cast<unsigned int>(cell->level()) ==
7416 input_triangulation.
n_levels() - 1,
7418 "Your input triangulation appears to have "
7419 "adaptively refined cells. This is not allowed. You can "
7420 "only call this function on a triangulation in which "
7421 "all cells are on the same refinement level."));
7425 this_cell.vertices[v] = cell->vertex_index(v);
7427 cells.push_back(this_cell);
7433 std::vector<unsigned int> considered_vertices;
7437 considered_vertices);
7449 const unsigned int n_slices,
7452 const bool copy_manifold_ids,
7453 const std::vector<types::manifold_id> &manifold_priorities)
7457 "The input triangulation must be a coarse mesh, i.e., it must "
7458 "not have been refined."));
7460 ExcMessage(
"The output triangulation object needs to be empty."));
7462 ExcMessage(
"The given height for extrusion must be positive."));
7465 "The number of slices for extrusion must be at least 2."));
7467 const double delta_h =
height / (n_slices - 1);
7468 std::vector<double> slices_z_values;
7469 slices_z_values.reserve(n_slices);
7470 for (
unsigned int i = 0; i < n_slices; ++i)
7471 slices_z_values.push_back(i * delta_h);
7473 input, slices_z_values, result, copy_manifold_ids, manifold_priorities);
7481 const unsigned int n_slices,
7484 const bool copy_manifold_ids,
7485 const std::vector<types::manifold_id> &manifold_priorities)
7491 (void)copy_manifold_ids;
7492 (void)manifold_priorities;
7496 "GridTools::extrude_triangulation() is only available "
7497 "for Triangulation<3, 3> as output triangulation."));
7505 const std::vector<double> &slice_coordinates,
7507 const bool copy_manifold_ids,
7508 const std::vector<types::manifold_id> &manifold_priorities)
7512 "The input triangulation must be a coarse mesh, i.e., it must "
7513 "not have been refined."));
7515 ExcMessage(
"The output triangulation object needs to be empty."));
7516 Assert(slice_coordinates.size() >= 2,
7518 "The number of slices for extrusion must be at least 2."));
7519 Assert(std::is_sorted(slice_coordinates.begin(), slice_coordinates.end()),
7520 ExcMessage(
"Slice z-coordinates should be in ascending order"));
7523 "This function is only implemented for quadrilateral meshes."));
7525 const auto priorities = [&]() -> std::vector<types::manifold_id> {
7529 if (0 < manifold_priorities.size())
7534 std::vector<types::manifold_id> sorted_manifold_priorities =
7535 manifold_priorities;
7536 std::sort(sorted_manifold_priorities.begin(),
7537 sorted_manifold_priorities.end());
7538 Assert(std::unique(sorted_manifold_priorities.begin(),
7539 sorted_manifold_priorities.end()) ==
7540 sorted_manifold_priorities.end(),
7542 "The given vector of manifold ids may not contain any "
7543 "duplicated entries."));
7544 std::vector<types::manifold_id> sorted_manifold_ids =
7546 std::sort(sorted_manifold_ids.begin(), sorted_manifold_ids.end());
7547 if (sorted_manifold_priorities != sorted_manifold_ids)
7549 std::ostringstream message;
7550 message <<
"The given triangulation has manifold ids {";
7552 sorted_manifold_ids)
7555 message << sorted_manifold_ids.back() <<
"}, but \n"
7556 <<
" the given vector of manifold ids is {";
7558 manifold_priorities)
7562 << manifold_priorities.back() <<
"}.\n"
7563 <<
" These vectors should contain the same elements.\n";
7564 const std::string m = message.str();
7568 return manifold_priorities;
7572 std::vector<types::manifold_id> default_priorities =
7574 const auto first_tfi_it = std::partition(
7575 default_priorities.begin(),
7576 default_priorities.end(),
7578 return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
7579 &input.get_manifold(id)) == nullptr;
7581 std::sort(default_priorities.begin(), first_tfi_it);
7582 std::sort(first_tfi_it, default_priorities.end());
7584 return default_priorities;
7587 const std::size_t n_slices = slice_coordinates.size();
7588 std::vector<Point<3>> points(n_slices * input.
n_vertices());
7589 std::vector<CellData<3>> cells;
7594 for (std::size_t slice_n = 0; slice_n < n_slices; ++slice_n)
7596 for (std::size_t vertex_n = 0; vertex_n < input.
n_vertices();
7600 points[slice_n * input.
n_vertices() + vertex_n] =
7601 Point<3>(vertex[0], vertex[1], slice_coordinates[slice_n]);
7607 for (
const auto &cell : input.active_cell_iterators())
7609 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
7612 for (
const unsigned int vertex_n :
7616 cell->vertex_index(vertex_n) + slice_n * input.
n_vertices();
7619 cell->vertex_index(vertex_n) +
7624 if (copy_manifold_ids)
7626 cells.push_back(this_cell);
7637 for (
const auto &face : input.active_face_iterators())
7641 if (face->at_boundary())
7643 if (copy_manifold_ids)
7645 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
7648 face->vertex_index(0) + slice_n * input.
n_vertices();
7650 face->vertex_index(1) + slice_n * input.
n_vertices();
7652 face->vertex_index(0) + (slice_n + 1) * input.
n_vertices();
7654 face->vertex_index(1) + (slice_n + 1) * input.
n_vertices();
7655 quads.push_back(quad);
7661 if (copy_manifold_ids)
7662 for (
const auto &cell : input.active_cell_iterators())
7667 for (std::size_t slice_n = 1; slice_n < n_slices - 1; ++slice_n)
7670 cell->vertex_index(0) + slice_n * input.
n_vertices();
7672 cell->vertex_index(1) + slice_n * input.
n_vertices();
7674 cell->vertex_index(2) + slice_n * input.
n_vertices();
7676 cell->vertex_index(3) + slice_n * input.
n_vertices();
7677 quads.push_back(quad);
7688 "The input triangulation to this function is using boundary "
7689 "indicators in a range that do not allow using "
7690 "max_boundary_id+1 and max_boundary_id+2 as boundary "
7691 "indicators for the bottom and top faces of the "
7692 "extruded triangulation."));
7695 for (
const auto &cell : input.active_cell_iterators())
7699 quad.
vertices[0] = cell->vertex_index(0);
7700 quad.
vertices[1] = cell->vertex_index(1);
7701 quad.
vertices[2] = cell->vertex_index(2);
7702 quad.
vertices[3] = cell->vertex_index(3);
7703 if (copy_manifold_ids)
7705 quads.push_back(quad);
7708 for (
unsigned int &vertex : quad.vertices)
7709 vertex += (n_slices - 1) * input.n_vertices();
7710 if (copy_manifold_ids)
7712 quads.push_back(quad);
7725 for (
auto manifold_id_it = priorities.rbegin();
7726 manifold_id_it != priorities.rend();
7728 for (
const auto &face : result.active_face_iterators())
7730 for (unsigned
int line_n = 0;
7731 line_n < GeometryInfo<3>::lines_per_face;
7733 face->line(line_n)->set_manifold_id(*manifold_id_it);
7741 const std::vector<double> &slice_coordinates,
7743 const bool copy_manifold_ids,
7744 const std::vector<types::manifold_id> &manifold_priorities)
7747 (void)slice_coordinates;
7749 (void)copy_manifold_ids;
7750 (void)manifold_priorities;
7754 "GridTools::extrude_triangulation() is only available "
7755 "for Triangulation<3, 3> as output triangulation."));
7774 template <
int spacedim>
7777 const double inner_radius,
7786 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7795 std::vector<bool> treated_vertices(triangulation.
n_vertices(),
false);
7796 for (
const auto &cell : triangulation.active_cell_iterators())
7799 if (cell->face(f)->at_boundary())
7808 std::
fabs(cell->face(f)->vertex(v)[1])) <
7810 cell->face(f)->vertex(v) *=
std::
sqrt(2.);
7813 for (
const auto &cell : triangulation.active_cell_iterators())
7815 for (
const unsigned int f : cell->face_indices())
7816 if (cell->face(f)->at_boundary())
7818 const double dx = cell->face(f)->center()[0] - center[0];
7819 const double dy = cell->face(f)->center()[1] - center[1];
7823 cell->face(f)->set_boundary_id(0);
7825 cell->face(f)->set_boundary_id(1);
7827 cell->face(f)->set_boundary_id(2);
7829 cell->face(f)->set_boundary_id(3);
7832 cell->face(f)->set_boundary_id(4);
7833 cell->face(f)->set_manifold_id(0);
7838 const double d = (cell->face(f)->center() - center).
norm();
7839 if (d - inner_radius < 0)
7841 cell->face(f)->set_boundary_id(1);
7842 cell->face(f)->set_manifold_id(0);
7845 cell->face(f)->set_boundary_id(0);
7857 const double inner_radius,
7860 const unsigned int width_repetition,
7863 hyper_cube_with_cylindrical_hole_2D(triangulation,
7876 const double inner_radius,
7879 const unsigned int width_repetition,
7882 hyper_cube_with_cylindrical_hole_2D(triangulation,
7896 const double inner_radius,
7898 const unsigned int n_shells,
7899 const double skewness,
7900 const unsigned int n_cells,
7905 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7909 std::vector<double> radii;
7910 radii.push_back(inner_radius);
7911 for (
unsigned int shell_n = 1; shell_n < n_shells; ++shell_n)
7912 if (skewness == 0.0)
7914 radii.push_back(inner_radius +
7916 (1.0 - (1.0 -
double(shell_n) / n_shells)));
7921 (1.0 -
std::tanh(skewness * (1.0 -
double(shell_n) / n_shells)) /
7925 double grid_vertex_tolerance = 0.0;
7926 for (
unsigned int shell_n = 0; shell_n < radii.size() - 1; ++shell_n)
7933 n_cells == 0 ? (dim == 2 ? 8 : 12) :
7938 if (grid_vertex_tolerance == 0.0)
7939 grid_vertex_tolerance =
7940 0.5 * internal::minimal_vertex_distance(current_shell);
7943 triangulation.
clear();
7947 grid_vertex_tolerance);
7959 constexpr double radial_vertex_tolerance =
7960 100.0 * std::numeric_limits<double>::epsilon();
7961 auto assert_vertex_distance_within_tolerance =
7962 [center, radial_vertex_tolerance](
7966 (void)radial_vertex_tolerance;
7967 for (
unsigned int vertex_n = 0;
7968 vertex_n < GeometryInfo<dim>::vertices_per_face;
7972 (center.
norm() +
radius) * radial_vertex_tolerance,
7977 for (
const auto &cell : triangulation.active_cell_iterators())
7980 auto face = cell->face(face_n);
7981 if (face->at_boundary())
7983 if (((face->vertex(0) - center).norm() - inner_radius) <
7984 (center.
norm() + inner_radius) * radial_vertex_tolerance)
7987 assert_vertex_distance_within_tolerance(face, inner_radius);
7988 face->set_all_boundary_ids(0);
7993 assert_vertex_distance_within_tolerance(face,
outer_radius);
7994 face->set_all_boundary_ids(1);
8005 const double inner_radius,
8008 const unsigned int Nz,
8014 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
8028 std::vector<bool> treated_vertices(triangulation.
n_vertices(),
false);
8029 for (
const auto &cell : triangulation.active_cell_iterators())
8031 for (
const auto f : cell->face_indices())
8032 if (cell->face(f)->at_boundary())
8036 const unsigned int vv = cell->face(f)->vertex_index(v);
8037 if (treated_vertices[vv] ==
false)
8039 treated_vertices[vv] =
true;
8049 cell->face(f)->vertex(v);
8050 if ((std::fabs(std::fabs(vertex_location[0]) -
8051 std::fabs(vertex_location[1])) <
8053 (std::fabs(vertex_location[0] * vertex_location[0] +
8054 vertex_location[1] * vertex_location[1] -
8057 cell->face(f)->vertex(v) =
8060 vertex_location[2]);
8066 for (
const auto &cell : triangulation.active_cell_iterators())
8068 for (
const unsigned int f : cell->face_indices())
8069 if (cell->face(f)->at_boundary())
8071 const double dx = cell->face(f)->center()[0];
8072 const double dy = cell->face(f)->center()[1];
8073 const double dz = cell->face(f)->center()[2];
8078 cell->face(f)->set_boundary_id(0);
8081 cell->face(f)->set_boundary_id(1);
8084 cell->face(f)->set_boundary_id(2);
8087 cell->face(f)->set_boundary_id(3);
8090 cell->face(f)->set_boundary_id(4);
8093 cell->face(f)->set_boundary_id(5);
8097 cell->face(f)->set_all_boundary_ids(6);
8098 cell->face(f)->set_all_manifold_ids(0);
8105 const double d = c.
norm();
8106 if (d - inner_radius < 0)
8108 cell->face(f)->set_all_boundary_ids(1);
8109 cell->face(f)->set_all_manifold_ids(0);
8112 cell->face(f)->set_boundary_id(0);
8121 template <
int dim,
int spacedim1,
int spacedim2>
8128 &in_tria) ==
nullptr),
8130 "This function cannot be used on "
8131 "parallel::distributed::Triangulation objects as inputs."));
8133 ExcMessage(
"This function does not work for meshes that have "
8137 const unsigned int spacedim =
std::min(spacedim1, spacedim2);
8138 const std::vector<Point<spacedim1>> &in_vertices = in_tria.
get_vertices();
8142 std::vector<Point<spacedim2>> v(in_vertices.size());
8143 for (
unsigned int i = 0; i < in_vertices.size(); ++i)
8144 for (
unsigned int d = 0;
d < spacedim; ++
d)
8145 v[i][d] = in_vertices[i][d];
8148 for (
const auto &cell : in_tria.active_cell_iterators())
8150 const unsigned int id = cell->active_cell_index();
8152 cells[id].vertices.resize(cell->n_vertices());
8154 cells[id].vertices[i] = cell->vertex_index(i);
8155 cells[id].material_id = cell->material_id();
8156 cells[id].manifold_id = cell->manifold_id();
8170 std::vector<bool> user_flags_line;
8173 .clear_user_flags_line();
8177 for (
const auto &face : in_tria.active_face_iterators())
8179 if (face->at_boundary())
8185 boundary_line.vertices[i] = face->vertex_index(i);
8190 std::move(boundary_line));
8203 if ((face->user_flag_set() ==
false) &&
8210 boundary_line.vertices[i] = face->vertex_index(i);
8216 std::move(boundary_line));
8218 face->set_user_flag();
8224 .load_user_flags_line(user_flags_line);
8231 std::vector<bool> user_flags_line;
8234 .clear_user_flags_line();
8236 std::vector<bool> user_flags_quad;
8239 .clear_user_flags_quad();
8243 for (
const auto &face : in_tria.active_face_iterators())
8245 if (face->at_boundary())
8251 boundary_face.vertices[i] = face->vertex_index(i);
8256 std::move(boundary_face));
8263 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
8264 if (face->line(e)->user_flag_set() ==
false)
8267 spacedim1>::line_iterator
8268 edge = face->line(e);
8273 boundary_edge.vertices[i] = edge->vertex_index(i);
8278 std::move(boundary_edge));
8280 edge->set_user_flag();
8298 if (face->user_flag_set() ==
false)
8306 boundary_face.vertices[i] = face->vertex_index(i);
8312 std::move(boundary_face));
8314 face->set_user_flag();
8323 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
8324 if (face->line(e)->at_boundary() ==
false)
8325 if (face->line(e)->user_flag_set() ==
false)
8328 line_iterator edge = face->line(e);
8333 boundary_edge.vertices[i] =
8334 edge->vertex_index(i);
8340 std::move(boundary_edge));
8342 edge->set_user_flag();
8349 .load_user_flags_line(user_flags_line);
8351 .load_user_flags_quad(user_flags_quad);
8361 for (
const auto i : out_tria.get_manifold_ids())
8363 out_tria.set_manifold(i,
FlatManifold<dim, spacedim2>());
8368 template <
int dim,
int spacedim>
8372 const unsigned int n_divisions)
8374 if constexpr (dim == 1)
8379 else if constexpr (dim == 2)
8381 n_divisions == 2 || n_divisions == 8,
8383 "Quadrilaterals must be split into either 2 or 8 triangles."));
8384 else if constexpr (dim == 3)
8385 AssertThrow(n_divisions == 6 || n_divisions == 24,
8387 "Hexahedra must be split into either 6 or 24 tetrahedra."));
8391 "GridGenerator::convert_hypercube_to_simplex_mesh() expects "
8392 "a Triangulation that consists only of quads/hexes."));
8410 {{{0, 1, 2}}, {{3, 2, 1}}}};
8420 [[maybe_unused]]
const auto vertex_ids_for_cells_2d =
8433 {{{0, 1, 12, 10}}, {{2, 3, 11, 12}}, {{7, 6, 11, 13}},
8434 {{5, 4, 13, 10}}, {{0, 2, 8, 12}}, {{4, 6, 13, 8}},
8435 {{5, 13, 7, 9}}, {{1, 9, 3, 12}}, {{0, 8, 4, 10}},
8436 {{1, 5, 9, 10}}, {{3, 7, 11, 9}}, {{2, 6, 8, 11}},
8437 {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
8438 {{12, 13, 8, 10}}, {{13, 8, 10, 4}}, {{13, 10, 9, 5}},
8439 {{13, 9, 11, 7}}, {{13, 11, 8, 6}}, {{10, 12, 9, 1}},
8440 {{9, 12, 11, 3}}, {{11, 12, 8, 2}}, {{8, 12, 10, 0}}}};
8442 [[maybe_unused]]
const auto vertex_ids_for_cells_3d =
8449 array<std::pair<unsigned int, std::array<unsigned int, 2>>, 4>
8450 vertex_ids_for_boundary_faces_2d_2 = {
8451 {{0, {{0, 2}}}, {1, {{1, 3}}}, {2, {{0, 1}}}, {3, {{2, 3}}}}};
8456 array<std::pair<unsigned int, std::array<unsigned int, 2>>, 8>
8457 vertex_ids_for_boundary_faces_2d_8 = {{{0, {{0, 4}}},
8465 [[maybe_unused]]
const auto vertex_ids_for_boundary_faces_2d =
8466 n_divisions == 2 ?
make_array_view(vertex_ids_for_boundary_faces_2d_2) :
8472 array<std::pair<unsigned int, std::array<unsigned int, 3>>, 12>
8473 vertex_ids_for_boundary_faces_3d_6 = {{{0, {{2, 6, 0}}},
8489 array<std::pair<unsigned int, std::array<unsigned int, 3>>, 24>
8490 vertex_ids_for_boundary_faces_3d_24 = {
8491 {{0, {{0, 4, 8}}}, {0, {{4, 8, 6}}}, {0, {{8, 6, 2}}},
8492 {0, {{0, 2, 8}}}, {1, {{1, 3, 9}}}, {1, {{3, 9, 7}}},
8493 {1, {{9, 7, 5}}}, {1, {{1, 9, 5}}}, {2, {{0, 1, 10}}},
8494 {2, {{1, 10, 5}}}, {2, {{10, 5, 4}}}, {2, {{0, 10, 4}}},
8495 {3, {{2, 3, 11}}}, {3, {{3, 11, 7}}}, {3, {{11, 7, 6}}},
8496 {3, {{2, 11, 6}}}, {4, {{0, 1, 12}}}, {4, {{1, 12, 3}}},
8497 {4, {{12, 3, 2}}}, {4, {{0, 12, 2}}}, {5, {{4, 5, 13}}},
8498 {5, {{5, 13, 7}}}, {5, {{13, 7, 6}}}, {5, {{4, 13, 6}}}}};
8499 [[maybe_unused]]
const auto vertex_ids_for_boundary_faces_3d =
8500 n_divisions == 6 ?
make_array_view(vertex_ids_for_boundary_faces_3d_6) :
8518 [[maybe_unused]]
const auto vertex_ids_for_inner_faces_2d =
8537 {{0, 12, 10}}, {{12, 1, 10}}, {{12, 1, 9}}, {{12, 3, 9}},
8538 {{12, 2, 11}}, {{12, 3, 11}}, {{12, 0, 8}}, {{12, 2, 8}},
8539 {{9, 13, 5}}, {{13, 7, 9}}, {{11, 7, 13}}, {{11, 6, 13}},
8540 {{4, 8, 13}}, {{6, 8, 13}}, {{4, 13, 10}}, {{13, 5, 10}},
8541 {{10, 9, 5}}, {{10, 9, 1}}, {{11, 9, 7}}, {{11, 9, 3}},
8542 {{8, 11, 2}}, {{8, 11, 6}}, {{8, 10, 0}}, {{8, 10, 4}},
8543 {{12, 3, 9}}, {{12, 9, 11}}, {{12, 3, 11}}, {{3, 9, 11}},
8544 {{2, 12, 8}}, {{2, 12, 11}}, {{2, 11, 8}}, {{8, 12, 11}},
8545 {{0, 12, 10}}, {{0, 12, 8}}, {{0, 8, 10}}, {{8, 10, 12}},
8546 {{12, 1, 10}}, {{12, 1, 9}}, {{1, 10, 9}}, {{10, 9, 12}},
8547 {{10, 8, 4}}, {{10, 8, 13}}, {{4, 13, 8}}, {{4, 13, 10}},
8548 {{10, 9, 13}}, {{10, 9, 5}}, {{13, 5, 10}}, {{13, 5, 9}},
8549 {{13, 7, 9}}, {{13, 7, 11}}, {{9, 11, 13}}, {{9, 11, 7}},
8550 {{8, 11, 13}}, {{8, 11, 6}}, {{6, 13, 8}}, {{6, 13, 11}},
8551 {{12, 13, 10}}, {{12, 13, 8}}, {{8, 10, 13}}, {{8, 10, 12}},
8552 {{12, 13, 10}}, {{12, 13, 9}}, {{10, 9, 13}}, {{10, 9, 12}},
8553 {{12, 13, 9}}, {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
8554 {{12, 13, 11}}, {{12, 13, 8}}, {{8, 11, 13}}, {{8, 11, 12}},
8556 [[maybe_unused]]
const auto vertex_ids_for_inner_faces_3d =
8567 {{{{12, 10}}, {{12, 9}}, {{12, 11}}, {{12, 8}}, {{9, 13}}, {{11, 13}},
8568 {{8, 13}}, {{10, 13}}, {{10, 9}}, {{9, 11}}, {{11, 8}}, {{8, 10}},
8569 {{12, 9}}, {{12, 11}}, {{11, 9}}, {{12, 8}}, {{12, 11}}, {{11, 8}},
8570 {{12, 8}}, {{12, 10}}, {{10, 8}}, {{12, 10}}, {{12, 9}}, {{9, 10}},
8571 {{13, 10}}, {{13, 8}}, {{8, 10}}, {{13, 10}}, {{13, 9}}, {{9, 10}},
8572 {{13, 11}}, {{13, 9}}, {{11, 9}}, {{13, 11}}, {{13, 8}}, {{11, 8}},
8573 {{12, 13}}, {{8, 10}}, {{8, 13}}, {{10, 13}}, {{8, 12}}, {{10, 12}},
8574 {{12, 13}}, {{10, 9}}, {{10, 13}}, {{9, 13}}, {{10, 12}}, {{9, 12}},
8575 {{12, 13}}, {{9, 11}}, {{9, 13}}, {{11, 13}}, {{9, 12}}, {{11, 12}},
8576 {{12, 13}}, {{11, 8}}, {{11, 13}}, {{8, 13}}, {{11, 12}}, {{8, 12}}}};
8577 [[maybe_unused]]
const auto vertex_ids_for_inner_edges_3d =
8589 array<std::pair<unsigned int, std::array<unsigned int, 2>>, 6>
8590 vertex_ids_for_new_boundary_edges_3d_6 = {{{0, {{0, 6}}},
8605 array<std::pair<unsigned int, std::array<unsigned int, 2>>, 24>
8606 vertex_ids_for_new_boundary_edges_3d_24 = {
8607 {{0, {{4, 8}}}, {0, {{6, 8}}}, {0, {{0, 8}}}, {0, {{2, 8}}},
8608 {1, {{5, 9}}}, {1, {{7, 9}}}, {1, {{1, 9}}}, {1, {{3, 9}}},
8609 {2, {{4, 10}}}, {2, {{5, 10}}}, {2, {{0, 10}}}, {2, {{1, 10}}},
8610 {3, {{6, 11}}}, {3, {{7, 11}}}, {3, {{2, 11}}}, {3, {{3, 11}}},
8611 {4, {{2, 12}}}, {4, {{3, 12}}}, {4, {{0, 12}}}, {4, {{1, 12}}},
8612 {5, {{6, 13}}}, {5, {{7, 13}}}, {5, {{4, 13}}}, {5, {{5, 13}}}}};
8613 [[maybe_unused]]
const auto vertex_ids_for_new_boundary_edges_3d =
8618 std::vector<Point<spacedim>> vertices;
8619 std::vector<CellData<dim>> cells;
8624 std::vector<unsigned int> old_to_new_vertex_indices(
8626 std::vector<unsigned int> face_to_new_vertex_indices(
8636 for (
const auto &cell : ref_tria.cell_iterators())
8640 std::array<unsigned int, dim == 2 ? 9 : 14> local_vertex_indices;
8646 const auto v_global = cell->vertex_index(v);
8648 if (old_to_new_vertex_indices[v_global] ==
8651 old_to_new_vertex_indices[v_global] = vertices.size();
8652 vertices.push_back(cell->vertex(v));
8656 local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
8660 if constexpr (dim > 1)
8661 for (
const auto f : cell->face_indices())
8663 const auto f_global = cell->face_index(f);
8665 if (face_to_new_vertex_indices[f_global] ==
8668 face_to_new_vertex_indices[f_global] = vertices.size();
8670 cell->face(f)->center(
true));
8674 local_vertex_indices.size());
8675 local_vertex_indices[cell->n_vertices() + f] =
8676 face_to_new_vertex_indices[f_global];
8683 local_vertex_indices.size());
8684 local_vertex_indices[cell->n_vertices() + cell->n_faces()] =
8686 vertices.push_back(cell->center(
true));
8690 [[maybe_unused]]
const auto add_cell =
8691 [&](
const unsigned int struct_dim,
8692 const auto &index_vertices,
8693 const unsigned int material_or_boundary_id,
8697 if (struct_dim < dim && (material_or_boundary_id ==
8702 if (struct_dim == dim)
8707 cell_data.material_id =
8708 material_or_boundary_id;
8710 for (
unsigned int i = 0; i < cell_data.vertices.size(); ++i)
8713 local_vertex_indices.size());
8714 cell_data.vertices[i] =
8715 local_vertex_indices[index_vertices[i]];
8717 cells.emplace_back(std::move(cell_data));
8719 else if (dim == 2 && struct_dim == 1)
8724 boundary_line.
boundary_id = material_or_boundary_id;
8726 for (
unsigned int i = 0; i < boundary_line.
vertices.
size(); ++i)
8729 local_vertex_indices.size());
8731 local_vertex_indices[index_vertices[i]];
8734 std::move(boundary_line));
8736 else if (dim == 3 && struct_dim == 2)
8741 boundary_quad.material_id = material_or_boundary_id;
8743 for (
unsigned int i = 0; i < boundary_quad.vertices.size(); ++i)
8746 local_vertex_indices.size());
8747 boundary_quad.vertices[i] =
8748 local_vertex_indices[index_vertices[i]];
8751 std::move(boundary_quad));
8753 else if (dim == 3 && struct_dim == 1)
8758 boundary_line.
boundary_id = material_or_boundary_id;
8760 for (
unsigned int i = 0; i < boundary_line.
vertices.
size(); ++i)
8763 local_vertex_indices.size());
8765 local_vertex_indices[index_vertices[i]];
8768 std::move(boundary_line));
8777 if constexpr (dim == 2)
8780 const auto manifold_id_cell = cell->manifold_id();
8781 const auto material_id_cell = cell->material_id();
8784 for (
const auto &cell_vertices : vertex_ids_for_cells_2d)
8785 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8788 for (
const auto &face_vertices : vertex_ids_for_inner_faces_2d)
8797 else if constexpr (dim == 3)
8800 const auto manifold_id_cell = cell->manifold_id();
8801 const auto material_id_cell = cell->material_id();
8803 for (
const auto &cell_vertices : vertex_ids_for_cells_3d)
8804 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8808 for (
const auto &face_vertices : vertex_ids_for_inner_faces_3d)
8816 for (
const auto &edge_vertices : vertex_ids_for_inner_edges_3d)
8826 if constexpr (dim == 2)
8828 for (
const auto &[quad_face_no, tri_face_vertices] :
8829 vertex_ids_for_boundary_faces_2d)
8831 const auto face = cell->face(quad_face_no);
8834 face->boundary_id(),
8835 face->manifold_id());
8838 else if constexpr (dim == 3)
8840 for (
const auto &[hex_face_no, tet_face_vertices] :
8841 vertex_ids_for_boundary_faces_3d)
8843 const auto face = cell->face(hex_face_no);
8846 face->boundary_id(),
8847 face->manifold_id());
8850 for (
const auto &[hex_face_no, tet_edge_vertices] :
8851 vertex_ids_for_new_boundary_edges_3d)
8853 const auto face = cell->face(hex_face_no);
8856 face->boundary_id(),
8857 face->manifold_id());
8865 if constexpr (dim == 3)
8867 for (
const auto e : cell->line_indices())
8869 auto edge = cell->line(e);
8875 old_to_new_vertex_indices[edge->vertex_index(0)];
8877 old_to_new_vertex_indices[edge->vertex_index(1)];
8889 for (
const auto i : out_tria.get_manifold_ids())
8891 out_tria.set_manifold(i,
FlatManifold<dim, spacedim>());
8896 template <
int dim,
int spacedim>
8901 if constexpr (dim == 1)
8911 "GridGenerator::convert_simplex_to_hypercube_mesh() expects "
8912 "a Triangulation that consists only of triangles/tetrahedra."));
8925 std::vector<Point<spacedim>> vertices = ref_tria.
get_vertices();
8927 (dim > 2 ? ref_tria.
n_quads() : 0) +
8929 std::map<typename Triangulation<dim, spacedim>::line_iterator,
8931 line_to_midpoint_vertex_map;
8932 std::map<typename Triangulation<dim, spacedim>::face_iterator,
8934 face_to_midpoint_vertex_map;
8935 std::map<typename Triangulation<dim, spacedim>::cell_iterator,
8937 cell_to_midpoint_vertex_map;
8938 for (
const auto &cell : ref_tria.active_cell_iterators())
8942 for (
const auto &face : cell->face_iterators())
8943 if (face_to_midpoint_vertex_map.find(face) ==
8944 face_to_midpoint_vertex_map.
end())
8949 vertices.push_back(new_vertex);
8953 face_to_midpoint_vertex_map[face] = vertices.size() - 1;
8957 if constexpr (dim == 3)
8959 for (
unsigned int l = 0;
l < face->n_lines(); ++
l)
8962 spacedim>::line_iterator
8963 line = face->line(l);
8964 if (line_to_midpoint_vertex_map.find(line) ==
8965 line_to_midpoint_vertex_map.end())
8971 vertices.push_back(new_vertex);
8976 line_to_midpoint_vertex_map[line] =
8977 vertices.size() - 1;
8987 vertices.push_back(new_vertex);
8991 cell_to_midpoint_vertex_map[cell] = vertices.size() - 1;
8994 Assert(vertices.size() ==
9002 std::vector<CellData<dim>> cell_data;
9003 cell_data.reserve(ref_tria.
n_cells() * (dim + 1));
9005 for (
const auto &cell : ref_tria.active_cell_iterators())
9006 for (unsigned
int v = 0; v < cell->n_vertices(); ++v)
9010 if constexpr (dim == 2)
9016 cell->vertex_index(v),
9017 face_to_midpoint_vertex_map[cell->line(v)],
9018 face_to_midpoint_vertex_map[cell->line((v + 2) % 3)],
9019 cell_to_midpoint_vertex_map[cell]};
9026 cell->vertex_index(0),
9027 line_to_midpoint_vertex_map[cell->line(0)],
9028 line_to_midpoint_vertex_map[cell->line(2)],
9029 face_to_midpoint_vertex_map[cell->face(0)],
9030 line_to_midpoint_vertex_map[cell->line(3)],
9031 face_to_midpoint_vertex_map[cell->face(1)],
9032 face_to_midpoint_vertex_map[cell->face(2)],
9033 cell_to_midpoint_vertex_map[cell]};
9038 line_to_midpoint_vertex_map[cell->line(0)],
9039 cell->vertex_index(1),
9040 face_to_midpoint_vertex_map[cell->face(0)],
9041 line_to_midpoint_vertex_map[cell->line(1)],
9042 face_to_midpoint_vertex_map[cell->face(1)],
9043 line_to_midpoint_vertex_map[cell->line(4)],
9044 cell_to_midpoint_vertex_map[cell],
9045 face_to_midpoint_vertex_map[cell->face(3)]};
9050 line_to_midpoint_vertex_map[cell->line(2)],
9051 face_to_midpoint_vertex_map[cell->face(0)],
9052 cell->vertex_index(2),
9053 line_to_midpoint_vertex_map[cell->line(1)],
9054 face_to_midpoint_vertex_map[cell->face(2)],
9055 cell_to_midpoint_vertex_map[cell],
9056 line_to_midpoint_vertex_map[cell->line(5)],
9057 face_to_midpoint_vertex_map[cell->face(3)]};
9062 line_to_midpoint_vertex_map[cell->line(3)],
9063 face_to_midpoint_vertex_map[cell->face(1)],
9064 face_to_midpoint_vertex_map[cell->face(2)],
9065 cell_to_midpoint_vertex_map[cell],
9066 cell->vertex_index(3),
9067 line_to_midpoint_vertex_map[cell->line(4)],
9068 line_to_midpoint_vertex_map[cell->line(5)],
9069 face_to_midpoint_vertex_map[cell->face(3)]};
9078 cell_data.emplace_back(std::move(new_cell));
9083 for (
const auto &cell : ref_tria.active_cell_iterators())
9084 for (
const auto &face : cell->face_iterators())
9089 if constexpr (dim == 2)
9091 CellData<dim - 1> child_face_1, child_face_2;
9092 child_face_1.
vertices = {face->vertex_index(0),
9093 face_to_midpoint_vertex_map[face]};
9094 child_face_1.boundary_id = face->boundary_id();
9095 child_face_2.vertices = {face_to_midpoint_vertex_map[face],
9096 face->vertex_index(1)};
9097 child_face_2.boundary_id = face->boundary_id();
9099 std::move(child_face_1));
9101 std::move(child_face_2));
9105 CellData<dim - 1> child_face_1, child_face_2, child_face_3;
9108 face->vertex_index(0),
9109 line_to_midpoint_vertex_map[face->line(0)],
9110 line_to_midpoint_vertex_map[face->line(2)],
9111 face_to_midpoint_vertex_map[face]};
9112 child_face_1.boundary_id = face->boundary_id();
9114 child_face_2.vertices = {
9115 line_to_midpoint_vertex_map[face->line(0)],
9116 face->vertex_index(1),
9117 face_to_midpoint_vertex_map[face],
9118 line_to_midpoint_vertex_map[face->line(1)]};
9120 child_face_2.boundary_id = face->boundary_id();
9122 child_face_3.vertices = {
9123 line_to_midpoint_vertex_map[face->line(2)],
9124 face_to_midpoint_vertex_map[face],
9125 face->vertex_index(2),
9126 line_to_midpoint_vertex_map[face->line(1)]};
9127 child_face_3.boundary_id = face->boundary_id();
9130 std::move(child_face_1));
9132 std::move(child_face_2));
9134 std::move(child_face_3));
9144 if constexpr (dim == 3)
9152 for (
const auto i : out_tria.get_manifold_ids())
9154 out_tria.set_manifold(i,
FlatManifold<dim, spacedim>());
9159 template <
int dim,
int spacedim>
9175 {{{0, 1, 3}}, {{1, 2, 3}}, {{2, 0, 3}}}};
9180 vertex_ids_for_boundary_faces_2d = {
9181 {{{{{0, 1}}}}, {{{{1, 2}}}}, {{{{2, 0}}}}}};
9185 {{{0, 1, 2, 4}}, {{1, 0, 3, 4}}, {{0, 2, 3, 4}}, {{2, 1, 3, 4}}}};
9190 vertex_ids_for_boundary_faces_3d = {
9191 {{{{{0, 1, 2}}}}, {{{{1, 0, 3}}}}, {{{{0, 2, 3}}}}, {{{{2, 1, 3}}}}}};
9197 vertex_ids_for_boundary_lines_3d = {{{{{{0, 1}}}},
9204 std::vector<Point<spacedim>> vertices;
9205 std::vector<CellData<dim>> cells;
9210 std::vector<unsigned int> old_to_new_vertex_indices(
9217 for (
const auto &cell : ref_tria.cell_iterators())
9220 cell->reference_cell().is_simplex(),
9222 "Cell with invalid ReferenceCell encountered. GridGenerator::alfeld_split_of_simplex_mesh() "
9223 "only supports simplex meshes as input."));
9227 std::array<
unsigned int, (dim == 3) ? 5 : 4> local_vertex_indices;
9233 const auto v_global = cell->vertex_index(v);
9235 if (old_to_new_vertex_indices[v_global] ==
9238 old_to_new_vertex_indices[v_global] = vertices.size();
9239 vertices.push_back(cell->vertex(v));
9243 local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
9249 local_vertex_indices[local_vertex_indices.size() - 1] = vertices.size();
9250 vertices.push_back(barycenter /
static_cast<double>(dim + 1));
9253 const auto add_cell = [&](
const unsigned int struct_dim,
9254 const auto &index_vertices,
9255 const unsigned int material_or_boundary_id,
9259 if (struct_dim < dim &&
9264 if (struct_dim == dim)
9269 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
9272 local_vertex_indices.size());
9273 cell_data.vertices[i] =
9274 local_vertex_indices[index_vertices[i]];
9276 cell_data.material_id =
9277 material_or_boundary_id;
9279 cells.emplace_back(std::move(cell_data));
9281 else if (dim == 2 && struct_dim == 1)
9286 boundary_line.
boundary_id = material_or_boundary_id;
9288 for (
unsigned int i = 0; i < boundary_line.
vertices.
size(); ++i)
9291 local_vertex_indices.size());
9293 local_vertex_indices[index_vertices[i]];
9296 std::move(boundary_line));
9298 else if (dim == 3 && struct_dim == 2)
9303 boundary_quad.material_id = material_or_boundary_id;
9305 for (
unsigned int i = 0; i < boundary_quad.vertices.size(); ++i)
9308 local_vertex_indices.size());
9309 boundary_quad.vertices[i] =
9310 local_vertex_indices[index_vertices[i]];
9313 std::move(boundary_quad));
9315 else if (dim == 3 && struct_dim == 1)
9320 boundary_line.
boundary_id = material_or_boundary_id;
9322 for (
unsigned int i = 0; i < boundary_line.
vertices.
size(); ++i)
9325 local_vertex_indices.size());
9327 local_vertex_indices[index_vertices[i]];
9330 std::move(boundary_line));
9338 const auto material_id_cell = cell->material_id();
9339 const auto manifold_id_cell = cell->manifold_id();
9344 for (
const auto &cell_vertices : vertex_ids_for_cells_2d)
9345 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
9349 for (
const auto &cell_vertices : vertex_ids_for_cells_3d)
9350 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
9356 for (
const auto f : cell->face_indices())
9358 const auto bid = cell->face(f)->boundary_id();
9359 const auto mid = cell->face(f)->manifold_id();
9364 for (
const auto &face_vertices :
9365 vertex_ids_for_boundary_faces_2d[f])
9366 add_cell(1, face_vertices, bid, mid);
9370 for (
const auto &face_vertices :
9371 vertex_ids_for_boundary_faces_3d[f])
9372 add_cell(2, face_vertices, bid, mid);
9381 for (
const auto l : cell->line_indices())
9383 const auto bid = cell->line(l)->boundary_id();
9384 const auto mid = cell->line(l)->manifold_id();
9386 for (
const auto &line_vertices :
9387 vertex_ids_for_boundary_lines_3d[
l])
9388 add_cell(1, line_vertices, bid, mid);
9396 for (
const auto i : out_tria.get_manifold_ids())
9398 out_tria.set_manifold(i,
FlatManifold<dim, spacedim>());
9403 template <
template <
int,
int>
class MeshType,
int dim,
int spacedim>
9407 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
9408 typename MeshType<dim, spacedim>::face_iterator>
9410 typename ExtractBoundaryMesh<MeshType, dim, spacedim>::return_type
9413 MeshType<dim - 1, spacedim> &surface_mesh,
9414 const std::set<types::boundary_id> &boundary_ids)
9418 &volume_mesh.get_triangulation()) ==
nullptr),
9427 constexpr unsigned int boundary_dim = dim - 1;
9434 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>
9435 temporary_mapping_level0;
9440 std::vector<bool> touched(volume_mesh.get_triangulation().n_vertices(),
9444 std::vector<CellData<boundary_dim>> cells;
9446 std::vector<Point<spacedim>> vertices;
9449 std::map<unsigned int, unsigned int> map_vert_index;
9462 (volume_mesh.get_triangulation().is_mixed_mesh() ==
false),
9464 "Boundary extraction for mixed meshes is currently not supported."));
9466 volume_mesh.begin(0)->reference_cell();
9471 const unsigned int max_n_face_vertices =
9475 max_n_face_vertices);
9477 for (
unsigned int j = 0; j < max_n_face_vertices; ++j)
9478 swap_matrix[f][j] = j;
9485 std::swap(swap_matrix[0][1], swap_matrix[0][2]);
9486 std::swap(swap_matrix[2][1], swap_matrix[2][2]);
9487 std::swap(swap_matrix[4][1], swap_matrix[4][2]);
9491 std::swap(swap_matrix[1][0], swap_matrix[1][1]);
9492 std::swap(swap_matrix[2][0], swap_matrix[2][1]);
9495 else if (reference_cell == ReferenceCells::get_simplex<dim>())
9502 std::swap(swap_matrix[f][0], swap_matrix[f][1]);
9512 for (
typename MeshType<dim, spacedim>::cell_iterator cell =
9513 volume_mesh.begin(0);
9514 cell != volume_mesh.end(0);
9516 for (
const unsigned int i : cell->
reference_cell().face_indices())
9518 const typename MeshType<dim, spacedim>::face_iterator face =
9521 if (face->at_boundary() &&
9522 (boundary_ids.empty() ||
9523 (boundary_ids.find(face->boundary_id()) != boundary_ids.end())))
9529 const unsigned int v_index = face->vertex_index(j);
9531 if (!touched[v_index])
9533 vertices.push_back(face->vertex(j));
9534 map_vert_index[v_index] = vertices.size() - 1;
9535 touched[v_index] =
true;
9538 c_data.vertices[swap_matrix[i][j]] = map_vert_index[v_index];
9540 c_data.material_id =
9542 c_data.manifold_id = face->manifold_id();
9552 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
9558 bool edge_found =
false;
9559 for (
auto &boundary_line : subcell_data.boundary_lines)
9560 if (((boundary_line.vertices[0] ==
9561 map_vert_index[face->line(
e)->vertex_index(0)]) &&
9562 (boundary_line.vertices[1] ==
9563 map_vert_index[face->line(
e)->vertex_index(
9565 ((boundary_line.vertices[0] ==
9566 map_vert_index[face->line(
e)->vertex_index(1)]) &&
9567 (boundary_line.vertices[1] ==
9568 map_vert_index[face->line(
e)->vertex_index(0)])))
9575 if (edge_found ==
true)
9582 map_vert_index[face->line(e)->vertex_index(0)];
9584 map_vert_index[face->line(e)->vertex_index(1)];
9591 cells.emplace_back(std::move(c_data));
9592 temporary_mapping_level0.push_back(std::make_pair(face, i));
9599 surface_mesh.get_triangulation())
9605 for (
const auto &cell : surface_mesh.active_cell_iterators())
9606 for (unsigned
int vertex = 0; vertex < 2; ++vertex)
9607 if (cell->face(vertex)->at_boundary())
9608 cell->face(vertex)->set_boundary_id(0);
9616 std::vector<std::pair<
9617 const typename MeshType<dim - 1, spacedim>::cell_iterator,
9618 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>>
9619 temporary_map_boundary_cell_face;
9620 for (
const auto &cell : surface_mesh.active_cell_iterators())
9621 temporary_map_boundary_cell_face.
push_back(
9622 std::make_pair(cell, temporary_mapping_level0.at(cell->
index())));
9634 unsigned int index_cells_deepest_level = 0;
9637 bool changed =
false;
9641 std::vector<unsigned int> cells_refined;
9645 for (
unsigned int cell_n = index_cells_deepest_level;
9646 cell_n < temporary_map_boundary_cell_face.size();
9651 if (temporary_map_boundary_cell_face[cell_n]
9652 .
second.first->has_children())
9656 Assert(temporary_map_boundary_cell_face[cell_n]
9657 .
second.first->refinement_case() ==
9660 temporary_map_boundary_cell_face[cell_n]
9661 .first->set_refine_flag();
9662 cells_refined.push_back(cell_n);
9673 surface_mesh.get_triangulation())
9674 .execute_coarsening_and_refinement();
9677 index_cells_deepest_level = temporary_map_boundary_cell_face.size();
9678 for (
const auto &refined_cell_n : cells_refined)
9680 const typename MeshType<dim - 1, spacedim>::cell_iterator
9682 temporary_map_boundary_cell_face[refined_cell_n].first;
9683 const typename MeshType<dim,
9684 spacedim>::face_iterator refined_face =
9685 temporary_map_boundary_cell_face[refined_cell_n].second.first;
9686 const unsigned int refined_face_number =
9687 temporary_map_boundary_cell_face[refined_cell_n]
9689 for (
unsigned int child_n = 0;
9690 child_n < refined_cell->n_children();
9695 temporary_map_boundary_cell_face.push_back(
9696 std::make_pair(refined_cell->child(
9697 swap_matrix[refined_face_number][child_n]),
9698 std::make_pair(refined_face->child(child_n),
9699 refined_face_number)));
9709 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
9710 typename MeshType<dim, spacedim>::face_iterator>
9711 surface_to_volume_mapping;
9712 for (
unsigned int i = 0; i < temporary_map_boundary_cell_face.size(); ++i)
9713 surface_to_volume_mapping[temporary_map_boundary_cell_face[i].
first] =
9714 temporary_map_boundary_cell_face[i].
second.first;
9718 const auto attached_mids =
9719 surface_mesh.get_triangulation().get_manifold_ids();
9720 for (
const auto i : volume_mesh.get_triangulation().get_manifold_ids())
9722 std::find(attached_mids.
begin(), attached_mids.
end(), i) ==
9723 attached_mids.
end())
9725 surface_mesh.get_triangulation())
9728 return surface_to_volume_mapping;
9733 template <
int dim,
int spacedim>
9737 const std::vector<unsigned int> &repetitions,
9744 std::vector<Point<spacedim>> vertices;
9745 std::vector<CellData<dim>> cells;
9751 (p2[1] - p1[1]) / repetitions[1]);
9754 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
9755 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
9760 for (
unsigned int j = 0; j < repetitions[1]; ++j)
9761 for (
unsigned int i = 0; i < repetitions[0]; ++i)
9764 std::array<unsigned int, 4> quad{{
9765 (j + 0) * (repetitions[0] + 1) + i + 0,
9766 (j + 0) * (repetitions[0] + 1) + i + 1,
9767 (j + 1) * (repetitions[0] + 1) + i + 0,
9768 (j + 1) * (repetitions[0] + 1) + i + 1
9774 tri.
vertices = {quad[0], quad[1], quad[2]};
9775 cells.emplace_back(std::move(tri));
9781 tri.
vertices = {quad[3], quad[2], quad[1]};
9782 cells.emplace_back(std::move(tri));
9790 (p2[1] - p1[1]) / repetitions[1],
9791 (p2[2] - p1[2]) / repetitions[2]);
9794 for (
unsigned int k = 0; k <= repetitions[2]; ++k)
9795 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
9796 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
9799 p1[2] + dx[2] * k));
9802 for (
unsigned int k = 0; k < repetitions[2]; ++k)
9803 for (
unsigned int j = 0; j < repetitions[1]; ++j)
9804 for (
unsigned int i = 0; i < repetitions[0]; ++i)
9807 std::array<unsigned int, 8> quad{
9808 {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9809 (j + 0) * (repetitions[0] + 1) + i + 0,
9810 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9811 (j + 0) * (repetitions[0] + 1) + i + 1,
9812 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9813 (j + 1) * (repetitions[0] + 1) + i + 0,
9814 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9815 (j + 1) * (repetitions[0] + 1) + i + 1,
9816 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9817 (j + 0) * (repetitions[0] + 1) + i + 0,
9818 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9819 (j + 0) * (repetitions[0] + 1) + i + 1,
9820 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9821 (j + 1) * (repetitions[0] + 1) + i + 0,
9822 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
9823 (j + 1) * (repetitions[0] + 1) + i + 1}};
9828 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9829 cell.
vertices = {{quad[0], quad[1], quad[2], quad[4]}};
9831 cell.
vertices = {{quad[0], quad[1], quad[3], quad[5]}};
9833 cells.emplace_back(std::move(cell));
9839 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9840 cell.
vertices = {{quad[2], quad[1], quad[3], quad[7]}};
9842 cell.
vertices = {{quad[0], quad[3], quad[2], quad[6]}};
9843 cells.emplace_back(std::move(cell));
9849 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9850 cell.
vertices = {{quad[1], quad[4], quad[5], quad[7]}};
9852 cell.
vertices = {{quad[0], quad[4], quad[5], quad[6]}};
9853 cells.emplace_back(std::move(cell));
9859 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9860 cell.
vertices = {{quad[2], quad[4], quad[7], quad[6]}};
9862 cell.
vertices = {{quad[3], quad[5], quad[7], quad[6]}};
9863 cells.emplace_back(std::move(cell));
9869 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
9870 cell.
vertices = {{quad[1], quad[2], quad[4], quad[7]}};
9872 cell.
vertices = {{quad[0], quad[3], quad[6], quad[5]}};
9873 cells.emplace_back(std::move(cell));
9895 double epsilon = std::numeric_limits<double>::max();
9896 for (
unsigned int i = 0; i < dim; ++i)
9898 0.01 * (
std::abs(p2[i] - p1[i]) / repetitions[i]));
9901 "The distance between corner points must be positive."));
9905 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
9911 template <
int dim,
int spacedim>
9914 const unsigned int repetitions,
9922 {{repetitions, repetitions}},
9931 {{repetitions, repetitions, repetitions}},
9933 {right, right, right},
9947 const double inner_radius,
9953 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
9958 for (
const auto &cell : tria.active_cell_iterators())
9961 if (cell->face(f)->at_boundary())
9966 const double ax = std::fabs(cell->face(f)->vertex(v)[0]);
9967 const double ay = std::fabs(cell->face(f)->vertex(v)[1]);
9968 const double az = std::fabs(cell->face(f)->vertex(v)[2]);
9969 const double max_coord =
std::max({ax, ay, az});
9978 for (
const auto &cell : tria.active_cell_iterators())
9980 if (cell->face(f)->at_boundary())
9982 bool vertex_on_sphere =
true;
9984 if (
std::
abs(cell->face(f)->vertex(v).
norm() - inner_radius) >
9987 vertex_on_sphere =
false;
9991 if (vertex_on_sphere)
9993 cell->face(f)->set_all_manifold_ids(0);
9994 cell->face(f)->set_boundary_id(
colorize ? 6 : 1);
9998 const Point<3> fc = cell->face(f)->center();
10000 cell->face(f)->set_boundary_id(0);
10002 cell->face(f)->set_boundary_id(1);
10004 cell->face(f)->set_boundary_id(2);
10006 cell->face(f)->set_boundary_id(3);
10008 cell->face(f)->set_boundary_id(4);
10010 cell->face(f)->set_boundary_id(5);
10044 const std::vector<unsigned int> &lengths_heights_widths,
10045 const double inner_radius,
10047 const bool use_transfinite_region,
10052 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
10055 const double length_pre = lengths_heights_widths[0] * 2 *
outer_radius;
10056 const double length_post = lengths_heights_widths[1] * 2 *
outer_radius;
10057 const double height_below = lengths_heights_widths[2] * 2 *
outer_radius;
10058 const double height_above = lengths_heights_widths[3] * 2 *
outer_radius;
10059 const double width_front = lengths_heights_widths[4] * 2 *
outer_radius;
10060 const double width_back = lengths_heights_widths[5] * 2 *
outer_radius;
10063 const unsigned int length_pre_repetitions = lengths_heights_widths[0];
10064 const unsigned int length_post_repetitions = lengths_heights_widths[1];
10065 const unsigned int height_below_repetitions = lengths_heights_widths[2];
10066 const unsigned int height_above_repetitions = lengths_heights_widths[3];
10067 const unsigned int width_front_repetitions = lengths_heights_widths[4];
10068 const unsigned int width_back_repetitions = lengths_heights_widths[5];
10073 {length_pre_repetitions + length_post_repetitions + 1,
10074 height_below_repetitions + height_above_repetitions + 1,
10075 width_front_repetitions + width_back_repetitions + 1},
10084 std::set<Triangulation<3>::active_cell_iterator> cells_to_remove;
10085 for (
const auto &cell : bulk_tria.active_cell_iterators())
10087 const Point<3> center = cell->center();
10091 cells_to_remove.insert(cell);
10097 tria_without_sphere);
10100 const double shell_radius = 0.5 * (inner_radius +
outer_radius);
10108 shell_tria,
Point<3>(), inner_radius, shell_radius, 1, 0.1, 6,
false);
10110 const double vertex_tol =
10111 std::min(internal::minimal_vertex_distance(shell_tria),
10112 internal::minimal_vertex_distance(sphere_tria)) *
10117 shell_tria, sphere_tria, temp, vertex_tol,
true);
10118 sphere_tria = std::move(temp);
10120 const double vertex_tolerance =
10121 std::min(internal::minimal_vertex_distance(tria_without_sphere),
10122 internal::minimal_vertex_distance(sphere_tria)) /
10126 tria_without_sphere, sphere_tria, tria, vertex_tolerance,
true);
10131 for (
const auto &cell : tria.cell_iterators())
10133 for (
const auto &face : cell->face_iterators())
10135 bool face_spherical =
true;
10138 if (face->vertex(v).norm_square() >
std::sqrt(inner_radius))
10139 face_spherical =
false;
10141 if (face_spherical)
10142 face->set_all_manifold_ids(1);
10148 if (use_transfinite_region)
10160 const double tol = 1
e-6;
10161 for (
const auto &face : tria.active_face_iterators())
10162 if (face->at_boundary())
10164 const Point<3> center = face->center();
10166 if (
std::abs(center[0] - (-
static_cast<double>(
10168 face->set_boundary_id(0);
10170 else if (
std::abs(center[0] -
static_cast<double>(length_post +
10173 face->set_boundary_id(1);
10175 else if (face->manifold_id() == 1)
10176 face->set_boundary_id(2);
10178 else if (
std::abs(center[1] - (-
static_cast<double>(
10181 face->set_boundary_id(3);
10183 else if (
std::abs(center[1] -
static_cast<double>(height_above +
10186 face->set_boundary_id(4);
10188 else if (
std::abs(center[2] - (-
static_cast<double>(
10191 face->set_boundary_id(5);
10193 else if (
std::abs(center[2] -
static_cast<double>(
10195 face->set_boundary_id(6);
10205 const std::vector<unsigned int> &,
10219 const std::vector<unsigned int> &,
10230# include "grid/grid_generator.inst"
* * for(const auto &cell :triangulation.active_cell_iterators())
* * const Number outer_radius
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
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)
static constexpr Point< dim, Number > unit_vector(const unsigned int i)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
const Point< spacedim > & get_center() const
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
unsigned int n_quads() 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
unsigned int n_lines() const
std::vector< Point< spacedim > > vertices
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
bool all_reference_cells_are_simplex() const
active_cell_iterator begin_active(const unsigned int level=0) const
constexpr size_type size() const noexcept
#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 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 uniform_channel_with_sphere(Triangulation< dim > &tria, const std::vector< unsigned int > &lengths_heights_widths, const double inner_radius, const double outer_radius, const bool use_transfinite_region=false, 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 subdivided_hyper_cube_with_simplices(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0.0, const double right=1.0, const bool colorize=false)
void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)
void hyper_cube_with_spherical_hole(Triangulation< dim > &tria, const double inner_radius=0.25, const double outer_radius=0.5, 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 uniform_channel_with_cylinder(Triangulation< dim > &tria, const std::vector< unsigned int > &lengths_and_heights, const double depth=1, const unsigned int depth_division=1, const double shell_region_radius=0.75, const unsigned int n_shells=2, const double skewness=2.0, const bool use_transfinite_region=false, 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 convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria, const unsigned int n_divisions=(dim==2 ? 8u :24u))
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 convert_simplex_to_hypercube_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
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 reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell< dim > &reference_cell)
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 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)
* * if(update_pressure &update_flags) * compute_pressure(constitutive_request
* * * * std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters const
constexpr const ReferenceCell< dim > & get_hypercube()
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
types::manifold_id manifold_id
std_cxx26::inplace_vector< unsigned int, ReferenceCells::max_n_vertices< structdim >()> vertices
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
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)