50 : airfoil_type(
"NACA")
52 , joukowski_center(-0.1, 0.14)
58 , incline_factor(0.35)
61 , n_subdivision_x_0(3)
62 , n_subdivision_x_1(2)
63 , n_subdivision_x_2(5)
65 , airfoil_sampling_factor(2)
68 airfoil_length <= height,
70 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
72 Assert(incline_factor < 1.0 && incline_factor >= 0.0,
73 ExcMessage(
"incline_factor has to be in [0,1)!"));
86 "Mesh height measured from airfoil nose to horizontal boundaries");
90 "Length measured from airfoil leading edge to vertical outlet boundary");
94 "Define obliqueness of the vertical mesh around the airfoil");
103 "Type of airfoil geometry, either NACA or Joukowski airfoil",
118 "Joukowski circle center coordinates");
121 "Joukowski airfoil length leading to trailing edge");
129 "Number of global refinements");
131 "NumberSubdivisionX0",
133 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
135 "NumberSubdivisionX1",
137 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
139 "NumberSubdivisionX2",
141 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
144 "Number of subdivisions normal to airfoil");
148 "Factor to obtain a finer mesh at the airfoil surface");
163 static const unsigned int id_block_1 = 1;
164 static const unsigned int id_block_2 = 2;
165 static const unsigned int id_block_3 = 3;
166 static const unsigned int id_block_4 = 4;
167 static const unsigned int id_block_5 = 5;
168 static const unsigned int id_block_6 = 6;
173 MeshGenerator(
const AdditionalData &data)
174 : refinements(data.refinements)
175 , n_subdivision_x_0(data.n_subdivision_x_0)
176 , n_subdivision_x_1(data.n_subdivision_x_1)
177 , n_subdivision_x_2(data.n_subdivision_x_2)
178 , n_subdivision_y(data.n_subdivision_y)
179 , height(data.height)
180 , length_b2(data.length_b2)
181 , incline_factor(data.incline_factor)
182 , bias_factor(data.bias_factor)
184 , n_cells_x_0(
Utilities::
pow(2, refinements) * n_subdivision_x_0)
185 , n_cells_x_1(
Utilities::
pow(2, refinements) * n_subdivision_x_1)
186 , n_cells_x_2(
Utilities::
pow(2, refinements) * n_subdivision_x_2)
187 , n_cells_y(
Utilities::
pow(2, refinements) * n_subdivision_y)
188 , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
190 , airfoil_1D(set_airfoil_length(
192 data.airfoil_type ==
"Joukowski" ?
193 joukowski(data.joukowski_center,
194 n_points_on_each_side,
195 data.airfoil_sampling_factor) :
196 (data.airfoil_type ==
"NACA" ?
198 n_points_on_each_side,
199 data.airfoil_sampling_factor) :
202 std::vector<Point<2>>{
206 data.airfoil_length))
207 , end_b0_x_u(airfoil_1D[0][n_cells_x_0][0])
208 , end_b0_x_l(airfoil_1D[1][n_cells_x_0][0])
209 , nose_x(airfoil_1D[0].front()[0])
210 , tail_x(airfoil_1D[0].back()[0])
211 , tail_y(airfoil_1D[0].back()[1])
212 , center_mesh(0.5 *
std::abs(end_b0_x_u + end_b0_x_l))
213 , length_b1_x(tail_x - center_mesh)
214 ,
gamma(std::atan(height /
215 (edge_length +
std::abs(nose_x - center_mesh))))
219 , A(nose_x - edge_length, 0)
222 , D(center_mesh, height)
224 ,
F(center_mesh, -height)
228 , J(tail_x + length_b2, 0)
232 Assert(data.airfoil_type ==
"Joukowski" ||
233 data.airfoil_type ==
"NACA",
246 make_coarse_grid(tria_grid);
248 set_boundary_ids(tria_grid);
250 if (periodic_faces !=
nullptr)
253 tria_grid, 5, 4, 1, *periodic_faces);
271 (void)periodic_faces;
278 const unsigned int refinements;
281 const unsigned int n_subdivision_x_0;
284 const unsigned int n_subdivision_x_1;
287 const unsigned int n_subdivision_x_2;
291 const unsigned int n_subdivision_y;
298 const double length_b2;
302 const double incline_factor;
305 const double bias_factor;
308 const double edge_length;
311 const unsigned int n_cells_x_0;
314 const unsigned int n_cells_x_1;
317 const unsigned int n_cells_x_2;
321 const unsigned int n_cells_y;
324 const unsigned int n_points_on_each_side;
328 const std::array<std::vector<Point<2>>, 2> airfoil_1D;
332 const double end_b0_x_u;
336 const double end_b0_x_l;
350 const double center_mesh;
353 const double length_b1_x;
380 const Point<2> A, B,
C, D,
E,
F, G, H, I, J,
K, L;
419 static std::array<std::vector<Point<2>>, 2>
420 joukowski(
const Point<2> ¢erpoint,
421 const unsigned int number_points,
422 const unsigned int factor)
424 std::array<std::vector<Point<2>>, 2> airfoil_1D;
425 const unsigned int total_points = 2 * number_points - 2;
426 const unsigned int n_airfoilpoints = factor * total_points;
428 const auto jouk_points =
429 joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
432 std::vector<Point<2>> upper_points;
433 std::vector<Point<2>> lower_points;
437 unsigned int nose_index = 0;
438 unsigned int tail_index = 0;
439 double nose_x_coordinate = 0;
440 double tail_x_coordinate = 0;
444 for (
unsigned int i = 0; i < jouk_points.size(); ++i)
446 if (jouk_points[i][0] < nose_x_coordinate)
448 nose_x_coordinate = jouk_points[i][0];
451 if (jouk_points[i][0] > tail_x_coordinate)
453 tail_x_coordinate = jouk_points[i][0];
459 for (
unsigned int i = tail_index; i < jouk_points.size(); ++i)
460 upper_points.emplace_back(jouk_points[i]);
461 for (
unsigned int i = 0; i <= nose_index; ++i)
462 upper_points.emplace_back(jouk_points[i]);
463 std::reverse(upper_points.begin(), upper_points.end());
466 lower_points.insert(lower_points.end(),
467 jouk_points.begin() + nose_index,
468 jouk_points.begin() + tail_index + 1);
471 airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
472 airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
475 auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
476 const double nose_x_pos = vector.front()[0];
477 for (
auto &i : vector)
481 move_nose_to_origin(airfoil_1D[1]);
482 move_nose_to_origin(airfoil_1D[0]);
511 static std::vector<Point<2>>
513 const unsigned int number_points)
515 std::vector<Point<2>> circle_points;
529 radius_test < radius,
531 "Error creating lower circle: Circle for Joukowski-transform does"
532 " not enclose point zeta = -1! Choose different center "
536 const double theta = 2 *
numbers::PI / number_points;
538 for (
unsigned int i = 0; i < number_points; ++i)
539 circle_points.emplace_back(
center[0] - radius *
cos(i * theta),
542 return circle_points;
553 static std::vector<Point<2>>
554 joukowski_transform(
const std::vector<
Point<2>> &circle_points)
556 std::vector<Point<2>> joukowski_points(circle_points.size());
559 for (
unsigned int i = 0; i < circle_points.size(); ++i)
561 const double chi = circle_points[i][0];
562 const double eta = circle_points[i][1];
563 const std::complex<double> zeta(chi, eta);
564 const std::complex<double> z = zeta + 1. / zeta;
566 joukowski_points[i] = {real(z), imag(z)};
568 return joukowski_points;
587 static std::array<std::vector<Point<2>>, 2>
588 naca(
const std::string &serialnumber,
589 const unsigned int number_points,
590 const unsigned int factor)
594 const unsigned int n_airfoilpoints = factor * number_points;
597 return {{make_points_equidistant(
598 naca_create_points(serialnumber, n_airfoilpoints,
true),
600 make_points_equidistant(
601 naca_create_points(serialnumber, n_airfoilpoints,
false),
616 static std::vector<Point<2>>
617 naca_create_points(
const std::string &serialnumber,
618 const unsigned int number_points,
621 Assert(serialnumber.size() == 4,
622 ExcMessage(
"This NACA-serial number is not implemented!"));
624 return naca_create_points_4_digits(serialnumber,
643 static std::vector<Point<2>>
644 naca_create_points_4_digits(
const std::string &serialnumber,
645 const unsigned int number_points,
649 const unsigned int digit_0 = (serialnumber[0] -
'0');
650 const unsigned int digit_1 = (serialnumber[1] -
'0');
651 const unsigned int digit_2 = (serialnumber[2] -
'0');
652 const unsigned int digit_3 = (serialnumber[3] -
'0');
654 const unsigned int digit_23 = 10 * digit_2 + digit_3;
657 const double t =
static_cast<double>(digit_23) / 100.0;
659 std::vector<Point<2>> naca_points;
661 if (digit_0 == 0 && digit_1 == 0)
662 for (
unsigned int i = 0; i < number_points; ++i)
664 const double x = i * 1 / (1.0 * number_points - 1);
668 0.3516 * Utilities::fixed_power<2>(x) +
669 0.2843 * Utilities::fixed_power<3>(x) -
670 0.1036 * Utilities::fixed_power<4>(
674 naca_points.emplace_back(x, +y_t);
676 naca_points.emplace_back(x, -y_t);
679 for (
unsigned int i = 0; i < number_points; ++i)
681 const double m = 1.0 * digit_0 / 100;
682 const double p = 1.0 * digit_1 / 10;
683 const double x = i * 1 / (1.0 * number_points - 1);
687 m / Utilities::fixed_power<2>(p) *
688 (2 *
p * x - Utilities::fixed_power<2>(x)) :
693 (x <=
p) ? 2 * m / Utilities::fixed_power<2>(p) * (
p - x) :
699 0.3516 * Utilities::fixed_power<2>(x) +
700 0.2843 * Utilities::fixed_power<3>(x) -
701 0.1036 * Utilities::fixed_power<4>(
704 const double theta = std::atan(dy_c);
707 naca_points.emplace_back(x - y_t *
std::sin(theta),
710 naca_points.emplace_back(x + y_t *
std::sin(theta),
727 static std::array<std::vector<Point<2>>, 2>
728 set_airfoil_length(
const std::array<std::vector<
Point<2>>, 2> &input,
729 const double desired_len)
731 std::array<std::vector<Point<2>>, 2> output;
732 output[0] = set_airfoil_length(input[0], desired_len);
733 output[1] = set_airfoil_length(input[1], desired_len);
745 static std::vector<Point<2>>
746 set_airfoil_length(
const std::vector<
Point<2>> &input,
747 const double desired_len)
749 std::vector<Point<2>> output = input;
752 desired_len / input.front().distance(input.back());
754 for (
auto &x : output)
770 static std::vector<Point<2>>
771 make_points_equidistant(
772 const std::vector<
Point<2>> &non_equidistant_points,
773 const unsigned int number_points)
775 const unsigned int n_points =
776 non_equidistant_points
780 std::vector<double> arclength_L(non_equidistant_points.size(), 0);
781 for (
unsigned int i = 0; i < non_equidistant_points.size() - 1; ++i)
784 non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
787 const auto airfoil_length =
789 const auto deltaX = airfoil_length / (number_points - 1);
793 std::vector<Point<2>> equidist(
796 equidist[0] = non_equidistant_points[0];
797 equidist[number_points - 1] = non_equidistant_points[n_points - 1];
801 for (
unsigned int j = 0, i = 1; j < n_points - 1; ++j)
804 const auto Lj = arclength_L[j];
805 const auto Ljp = arclength_L[j + 1];
807 while (Lj <= i * deltaX && i * deltaX <= Ljp &&
808 i < number_points - 1)
810 equidist[i] =
Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
811 (non_equidistant_points[j + 1] -
812 non_equidistant_points[j]) +
813 non_equidistant_points[j]);
833 std::vector<Triangulation<2>> trias(10);
837 const std::vector<Point<2>> &corner_vertices,
838 const std::vector<unsigned int> &repetitions,
851 const double xi =
point[0];
852 const double eta =
point[1];
855 point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
856 (1 + xi) * (1 - eta) * corner_vertices[1] +
857 (1 - xi) * (1 + eta) * corner_vertices[2] +
858 (1 + xi) * (1 + eta) * corner_vertices[3]);
862 for (
auto cell :
tria.active_cell_iterators())
870 {n_subdivision_y, n_subdivision_x_0},
874 {n_subdivision_y, n_subdivision_x_0},
878 {n_subdivision_x_1, n_subdivision_y},
882 {n_subdivision_x_1, n_subdivision_y},
886 {n_subdivision_x_2, n_subdivision_y},
890 {n_subdivision_x_2, n_subdivision_y},
917 for (
auto cell :
tria.active_cell_iterators())
920 if (cell->face(f)->at_boundary() ==
false)
923 const auto mid = cell->material_id();
925 if ((mid == id_block_1 && f == 0) ||
926 (mid == id_block_4 && f == 0))
927 cell->face(f)->set_boundary_id(0);
928 else if ((mid == id_block_3 && f == 0) ||
929 (mid == id_block_6 && f == 2))
930 cell->face(f)->set_boundary_id(1);
931 else if ((mid == id_block_1 && f == 1) ||
932 (mid == id_block_2 && f == 1))
933 cell->face(f)->set_boundary_id(2);
934 else if ((mid == id_block_4 && f == 1) ||
935 (mid == id_block_5 && f == 3))
936 cell->face(f)->set_boundary_id(3);
937 else if ((mid == id_block_2 && f == 0) ||
938 (mid == id_block_3 && f == 2))
939 cell->face(f)->set_boundary_id(4);
940 else if ((mid == id_block_5 && f == 2) ||
941 (mid == id_block_6 && f == 0))
942 cell->face(f)->set_boundary_id(5);
979 for (
const auto &cell :
tria.cell_iterators())
983 if (vertex_processed[cell->vertex_index(v)])
987 vertex_processed[cell->vertex_index(v)] =
true;
989 auto &node = cell->vertex(v);
992 if (cell->material_id() == id_block_1 ||
993 cell->material_id() == id_block_4)
1002 if (cell->material_id() == id_block_1)
1005 (node - horizontal_offset) +
1011 else if (cell->material_id() == id_block_4)
1014 (node - horizontal_offset) -
1020 const double trapeze_height =
1022 const double L = height /
std::sin(gamma);
1023 const double l_a =
std::cos(gamma) * edge_length;
1024 const double l_b = trapeze_height *
std::tan(gamma);
1026 const double x2 = L - l_a - l_b;
1028 const double Dx = x1 + x2 + x3;
1029 const double deltax =
1031 const double dx = Dx / n_cells_x_0;
1032 const double dy = trapeze_height / n_cells_y;
1034 static_cast<int>(std::round((node_[0] - deltax) / dx));
1036 static_cast<int>(std::round(
std::abs(node_[1]) / dy));
1038 node_[0] =
numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1039 node_[1] = height * (1.0 * iy) / n_cells_y;
1046 const double dy = height / n_cells_y;
1048 static_cast<int>(std::round(node_[0] / dx));
1050 static_cast<int>(std::round(node_[1] / dy));
1051 const double alpha =
1052 bias_alpha(1 - (1.0 * iy) / n_cells_y);
1053 const double theta = node_[0];
1055 ((cell->material_id() == id_block_1) ?
1060 (cell->material_id() == id_block_1) ? (0) : (1))]
1066 else if (cell->material_id() == id_block_2 ||
1067 cell->material_id() == id_block_5)
1075 "Points D,C,G and E,F,I are not defined symmetric to "
1076 "x-axis, which is required to interpolate block 2"
1077 " and 5 with same geometric computations."));
1078 const double l_y = D[1] -
C[1];
1079 const double l_h = D[1] - l_y;
1080 const double by = -l_h / length_b1_x * (node[0] - H[0]);
1081 const double dy = (height - by) / n_cells_y;
1082 const int iy =
static_cast<int>(
1083 std::round((
std::abs(node[1]) - by) / dy));
1084 const double dx = length_b1_x / n_cells_x_1;
1085 const int ix =
static_cast<int>(
1086 std::round(
std::abs(node[0] - center_mesh) / dx));
1088 const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1092 const Point<2> p(ix * dx + center_mesh +
1093 incline_factor * length_b2 * ix /
1095 ((cell->material_id() == id_block_2) ?
1101 (cell->material_id() == id_block_2) ? (0) : (1))]
1102 [n_cells_x_0 + ix] *
1106 else if (cell->material_id() == id_block_3 ||
1107 cell->material_id() == id_block_6)
1110 const double dx = length_b2 / n_cells_x_2;
1111 const double dy = height / n_cells_y;
1112 const int ix =
static_cast<int>(
1113 std::round(
std::abs(node[0] - H[0]) / dx));
1115 static_cast<int>(std::round(
std::abs(node[1]) / dy));
1117 const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1118 const double alpha_x =
1119 bias_alpha(1 - (
static_cast<double>(ix)) / n_cells_x_2);
1123 const Point<2> p1(J[0] - (1 - incline_factor) * length_b2 *
1125 ((cell->material_id() == id_block_3) ?
1130 const Point<2> p2(J[0] - alpha_x * length_b2, tail_y);
1131 node = p1 * (1 - alpha_y) + p2 * alpha_y;
1151 bias_alpha(
double alpha)
const
1153 return std::tanh(bias_factor * alpha) / std::tanh(bias_factor);
1161 internal_create_triangulation(
1165 const AdditionalData &additional_data)
1167 MeshGenerator mesh_generator(additional_data);
1170 if (
auto *parallel_tria =
1173 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1174 else if (
auto *parallel_tria =
dynamic_cast<
1177 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1179 mesh_generator.create_triangulation(tria, periodic_faces);
1196 const AdditionalData &)
1206 const AdditionalData &additional_data)
1208 internal_create_triangulation(tria,
nullptr, additional_data);
1219 const AdditionalData &additional_data)
1221 internal_create_triangulation(tria, &periodic_faces, additional_data);
1232 const AdditionalData &additional_data)
1236 (void)additional_data;
1237 (void)periodic_faces;
1248 template <
int dim,
int spacedim>
1260 cell->face(f)->set_boundary_id(f);
1266 template <
int spacedim>
1277 if (cell->center()[0] > 0)
1278 cell->set_material_id(1);
1285 template <
int dim,
int spacedim>
1290 const double epsilon)
1303 for (; face != endface; ++face)
1304 if (face->at_boundary())
1305 if (face->boundary_id() == 0)
1310 face->set_boundary_id(0);
1312 face->set_boundary_id(1);
1314 face->set_boundary_id(2);
1316 face->set_boundary_id(3);
1318 face->set_boundary_id(4);
1320 face->set_boundary_id(5);
1329 for (
const auto &cell :
tria.cell_iterators())
1332 for (
unsigned int d = 0;
d < dim; ++
d)
1333 if (cell->center()[d] > 0)
1335 cell->set_material_id(
id);
1361 cell->face(2)->set_all_boundary_ids(1);
1385 cell->face(4)->set_all_boundary_ids(1);
1389 cell->face(2)->set_all_boundary_ids(1);
1393 cell->face(2)->set_all_boundary_ids(1);
1397 cell->face(0)->set_all_boundary_ids(1);
1401 cell->face(2)->set_all_boundary_ids(1);
1405 cell->face(0)->set_all_boundary_ids(1);
1416 cell->face(5)->set_all_boundary_ids(1);
1426 unsigned int count = 0;
1427 for (
const auto &cell :
tria.cell_iterators())
1428 if (cell->face(5)->at_boundary())
1430 cell->face(5)->set_all_boundary_ids(1);
1450 const double inner_radius,
1451 const double outer_radius)
1456 double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
1457 double eps = 1
e-3 * middle;
1460 for (; cell !=
tria.
end(); ++cell)
1461 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
1463 if (!cell->face(f)->at_boundary())
1466 double radius = cell->face(f)->center().norm() -
center.
norm();
1467 if (std::fabs(cell->face(f)->center()[0]) <
1470 cell->face(f)->set_boundary_id(2);
1471 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1473 if (cell->face(f)->line(j)->at_boundary())
1474 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1475 cell->face(f)->line(j)->vertex(1).norm()) >
1477 cell->face(f)->line(j)->set_boundary_id(2);
1479 else if (std::fabs(cell->face(f)->center()[1]) <
1482 cell->face(f)->set_boundary_id(3);
1483 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1485 if (cell->face(f)->line(j)->at_boundary())
1486 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1487 cell->face(f)->line(j)->vertex(1).norm()) >
1489 cell->face(f)->line(j)->set_boundary_id(3);
1491 else if (std::fabs(cell->face(f)->center()[2]) <
1494 cell->face(f)->set_boundary_id(4);
1495 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1497 if (cell->face(f)->line(j)->at_boundary())
1498 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1499 cell->face(f)->line(j)->vertex(1).norm()) >
1501 cell->face(f)->line(j)->set_boundary_id(4);
1503 else if (radius < middle)
1505 cell->face(f)->set_boundary_id(0);
1506 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1508 if (cell->face(f)->line(j)->at_boundary())
1509 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1510 cell->face(f)->line(j)->vertex(1).norm()) <
1512 cell->face(f)->line(j)->set_boundary_id(0);
1514 else if (radius > middle)
1516 cell->face(f)->set_boundary_id(1);
1517 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1519 if (cell->face(f)->line(j)->at_boundary())
1520 if (std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1521 cell->face(f)->line(j)->vertex(1).norm()) <
1523 cell->face(f)->line(j)->set_boundary_id(1);
1533 template <
int dim,
int spacedim>
1544 for (
unsigned int i = 0; i < dim; ++i)
1584 std::vector<CellData<dim>> cells(1);
1587 cells[0].material_id = 0;
1593 colorize_hyper_rectangle(tria);
1598 template <
int dim,
int spacedim>
1606 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1609 for (
unsigned int i = 0; i < dim; ++i)
1629 for (
unsigned int d = 0;
d < dim; ++
d)
1630 for (
unsigned int c = 1; c <= dim; ++c)
1633 ExcMessage(
"Vertices of simplex must form a right handed system"));
1637 std::vector<Point<dim>> points =
vertices;
1641 for (
unsigned int i = 0; i <= dim; ++i)
1643 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1649 for (
unsigned int i = 1; i < dim; ++i)
1650 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1652 for (
unsigned int i = 0; i <= dim; ++i)
1653 points.push_back(1. / 3. *
1654 (points[i] + points[(i + 1) % (dim + 1)] +
1655 points[(i + 2) % (dim + 1)]));
1657 points.push_back((1. / (dim + 1)) *
center);
1659 std::vector<CellData<dim>> cells(dim + 1);
1664 cells[0].vertices[0] = 0;
1665 cells[0].vertices[1] = 3;
1666 cells[0].vertices[2] = 5;
1667 cells[0].vertices[3] = 6;
1668 cells[0].material_id = 0;
1670 cells[1].vertices[0] = 3;
1671 cells[1].vertices[1] = 1;
1672 cells[1].vertices[2] = 6;
1673 cells[1].vertices[3] = 4;
1674 cells[1].material_id = 0;
1676 cells[2].vertices[0] = 5;
1677 cells[2].vertices[1] = 6;
1678 cells[2].vertices[2] = 2;
1679 cells[2].vertices[3] = 4;
1680 cells[2].material_id = 0;
1684 cells[0].vertices[0] = 0;
1685 cells[0].vertices[1] = 4;
1686 cells[0].vertices[2] = 8;
1687 cells[0].vertices[3] = 10;
1688 cells[0].vertices[4] = 7;
1689 cells[0].vertices[5] = 13;
1690 cells[0].vertices[6] = 12;
1691 cells[0].vertices[7] = 14;
1692 cells[0].material_id = 0;
1694 cells[1].vertices[0] = 4;
1695 cells[1].vertices[1] = 1;
1696 cells[1].vertices[2] = 10;
1697 cells[1].vertices[3] = 5;
1698 cells[1].vertices[4] = 13;
1699 cells[1].vertices[5] = 9;
1700 cells[1].vertices[6] = 14;
1701 cells[1].vertices[7] = 11;
1702 cells[1].material_id = 0;
1704 cells[2].vertices[0] = 8;
1705 cells[2].vertices[1] = 10;
1706 cells[2].vertices[2] = 2;
1707 cells[2].vertices[3] = 5;
1708 cells[2].vertices[4] = 12;
1709 cells[2].vertices[5] = 14;
1710 cells[2].vertices[6] = 6;
1711 cells[2].vertices[7] = 11;
1712 cells[2].material_id = 0;
1714 cells[3].vertices[0] = 7;
1715 cells[3].vertices[1] = 13;
1716 cells[3].vertices[2] = 12;
1717 cells[3].vertices[3] = 14;
1718 cells[3].vertices[4] = 3;
1719 cells[3].vertices[5] = 9;
1720 cells[3].vertices[6] = 6;
1721 cells[3].vertices[7] = 11;
1722 cells[3].material_id = 0;
1732 template <
int dim,
int spacedim>
1739 if (reference_cell == ReferenceCells::get_hypercube<dim>())
1752 for (
unsigned int d = 0;
d < dim; ++
d)
1761 std::vector<CellData<dim>> cells(1);
1773 const unsigned int n_cells,
1774 const unsigned int n_rotations,
1778 const unsigned int dim = 3;
1781 "More than 4 cells are needed to create a moebius grid."));
1783 ExcMessage(
"Outer and inner radius must be positive."));
1785 ExcMessage(
"Outer radius must be greater than inner radius."));
1788 std::vector<Point<dim>>
vertices(4 * n_cells);
1792 for (
unsigned int i = 0; i <
n_cells; ++i)
1793 for (
unsigned int j = 0; j < 4; ++j)
1807 unsigned int offset = 0;
1812 static constexpr std::array<unsigned int, 8> local_vertex_numbering{
1813 {0, 1, 5, 4, 2, 3, 7, 6}};
1814 std::vector<CellData<dim>> cells(n_cells);
1815 for (
unsigned int i = 0; i <
n_cells; ++i)
1817 for (
unsigned int j = 0; j < 2; ++j)
1819 cells[i].vertices[local_vertex_numbering[0 + 4 * j]] =
1821 cells[i].vertices[local_vertex_numbering[1 + 4 * j]] =
1823 cells[i].vertices[local_vertex_numbering[2 + 4 * j]] =
1825 cells[i].vertices[local_vertex_numbering[3 + 4 * j]] =
1829 cells[i].material_id = 0;
1833 cells[
n_cells - 1].vertices[local_vertex_numbering[4]] =
1834 (0 + n_rotations) % 4;
1835 cells[
n_cells - 1].vertices[local_vertex_numbering[5]] =
1836 (3 + n_rotations) % 4;
1837 cells[
n_cells - 1].vertices[local_vertex_numbering[6]] =
1838 (2 + n_rotations) % 4;
1839 cells[
n_cells - 1].vertices[local_vertex_numbering[7]] =
1840 (1 + n_rotations) % 4;
1857 ExcMessage(
"Outer radius R must be greater than the inner "
1861 const unsigned int dim = 2;
1862 const unsigned int spacedim = 3;
1863 std::vector<Point<spacedim>>
vertices(16);
1882 std::vector<CellData<dim>> cells(16);
1884 cells[0].vertices[0] = 0;
1885 cells[0].vertices[1] = 4;
1886 cells[0].vertices[2] = 3;
1887 cells[0].vertices[3] = 7;
1888 cells[0].material_id = 0;
1890 cells[1].vertices[0] = 1;
1891 cells[1].vertices[1] = 5;
1892 cells[1].vertices[2] = 0;
1893 cells[1].vertices[3] = 4;
1894 cells[1].material_id = 0;
1896 cells[2].vertices[0] = 2;
1897 cells[2].vertices[1] = 6;
1898 cells[2].vertices[2] = 1;
1899 cells[2].vertices[3] = 5;
1900 cells[2].material_id = 0;
1902 cells[3].vertices[0] = 3;
1903 cells[3].vertices[1] = 7;
1904 cells[3].vertices[2] = 2;
1905 cells[3].vertices[3] = 6;
1906 cells[3].material_id = 0;
1908 cells[4].vertices[0] = 4;
1909 cells[4].vertices[1] = 8;
1910 cells[4].vertices[2] = 7;
1911 cells[4].vertices[3] = 11;
1912 cells[4].material_id = 0;
1914 cells[5].vertices[0] = 5;
1915 cells[5].vertices[1] = 9;
1916 cells[5].vertices[2] = 4;
1917 cells[5].vertices[3] = 8;
1918 cells[5].material_id = 0;
1920 cells[6].vertices[0] = 6;
1921 cells[6].vertices[1] = 10;
1922 cells[6].vertices[2] = 5;
1923 cells[6].vertices[3] = 9;
1924 cells[6].material_id = 0;
1926 cells[7].vertices[0] = 7;
1927 cells[7].vertices[1] = 11;
1928 cells[7].vertices[2] = 6;
1929 cells[7].vertices[3] = 10;
1930 cells[7].material_id = 0;
1932 cells[8].vertices[0] = 8;
1933 cells[8].vertices[1] = 12;
1934 cells[8].vertices[2] = 11;
1935 cells[8].vertices[3] = 15;
1936 cells[8].material_id = 0;
1938 cells[9].vertices[0] = 9;
1939 cells[9].vertices[1] = 13;
1940 cells[9].vertices[2] = 8;
1941 cells[9].vertices[3] = 12;
1942 cells[9].material_id = 0;
1944 cells[10].vertices[0] = 10;
1945 cells[10].vertices[1] = 14;
1946 cells[10].vertices[2] = 9;
1947 cells[10].vertices[3] = 13;
1948 cells[10].material_id = 0;
1950 cells[11].vertices[0] = 11;
1951 cells[11].vertices[1] = 15;
1952 cells[11].vertices[2] = 10;
1953 cells[11].vertices[3] = 14;
1954 cells[11].material_id = 0;
1956 cells[12].vertices[0] = 12;
1957 cells[12].vertices[1] = 0;
1958 cells[12].vertices[2] = 15;
1959 cells[12].vertices[3] = 3;
1960 cells[12].material_id = 0;
1962 cells[13].vertices[0] = 13;
1963 cells[13].vertices[1] = 1;
1964 cells[13].vertices[2] = 12;
1965 cells[13].vertices[3] = 0;
1966 cells[13].material_id = 0;
1968 cells[14].vertices[0] = 14;
1969 cells[14].vertices[1] = 2;
1970 cells[14].vertices[2] = 13;
1971 cells[14].vertices[3] = 1;
1972 cells[14].material_id = 0;
1974 cells[15].vertices[0] = 15;
1975 cells[15].vertices[1] = 3;
1976 cells[15].vertices[2] = 14;
1977 cells[15].vertices[3] = 2;
1978 cells[15].material_id = 0;
1991 static constexpr int circle_cell_vertices[5][4] = {{0, 1, 2, 3},
2005 const unsigned int n_cells_toroidal,
2009 ExcMessage(
"Outer radius R must be greater than the inner "
2013 ExcMessage(
"Number of cells in toroidal direction has "
2014 "to be at least 3 for a torus of polar extent 2*pi."));
2020 const double a = 1. / (1 +
std::sqrt(2.0));
2023 const unsigned int additional_layer =
2027 const unsigned int n_point_layers_toroidal =
2028 n_cells_toroidal + additional_layer;
2029 std::vector<Point<3>>
vertices(8 * n_point_layers_toroidal);
2041 const double phi_cell = phi / n_cells_toroidal;
2042 for (
unsigned int c = 1; c < n_point_layers_toroidal; ++c)
2044 for (
unsigned int v = 0; v < 8; ++v)
2046 const double r_2d =
vertices[v][0];
2054 std::vector<CellData<3>> cells(5 * n_cells_toroidal);
2055 for (
unsigned int c = 0; c < n_cells_toroidal; ++c)
2057 for (
unsigned int j = 0; j < 2; ++j)
2059 const unsigned int offset =
2060 (8 * (c + j)) % (8 * n_point_layers_toroidal);
2063 for (
unsigned int c2 = 0; c2 < 5; ++c2)
2064 for (
unsigned int i = 0; i < 4; ++i)
2065 cells[5 * c + c2].
vertices[i + j * 4] =
2066 offset + circle_cell_vertices[c2][i];
2069 cells[5 * c].material_id = 0;
2071 cells[5 * c + 1].material_id = 0;
2072 cells[5 * c + 2].material_id = 1;
2073 cells[5 * c + 3].material_id = 0;
2074 cells[5 * c + 4].material_id = 0;
2082 for (
const auto &cell :
tria.cell_iterators())
2085 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
2089 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2091 cell->face(f)->set_all_manifold_ids(1);
2096 if (cell->material_id() == 1)
2098 cell->set_all_manifold_ids(2);
2100 cell->set_material_id(0);
2117 template <
int dim,
int spacedim>
2138 "The volume of the cell is not greater than zero. "
2139 "This could be due to the wrong ordering of the vertices."));
2170 std::array<Tensor<1, 2>, 2> edges;
2171 edges[0] = corners[0];
2172 edges[1] = corners[1];
2173 std::vector<unsigned int> subdivisions;
2174 subdivided_parallelepiped<2, 2>(
2175 tria, origin, edges, subdivisions,
colorize);
2186 unsigned int n_subdivisions[dim];
2187 for (
unsigned int i = 0; i < dim; ++i)
2188 n_subdivisions[i] = 1;
2197 const unsigned int n_subdivisions,
2203 unsigned int n_subdivisions_[dim];
2204 for (
unsigned int i = 0; i < dim; ++i)
2205 n_subdivisions_[i] = n_subdivisions;
2215 const unsigned int (&n_subdivisions)[dim],
2217 const unsigned int *n_subdivisions,
2223 std::vector<unsigned int> subdivisions;
2224 std::array<Tensor<1, dim>, dim> edges;
2225 for (
unsigned int i = 0; i < dim; ++i)
2227 subdivisions.push_back(n_subdivisions[i]);
2228 edges[i] = corners[i];
2231 subdivided_parallelepiped<dim, dim>(
2232 tria, origin, edges, subdivisions,
colorize);
2238 template <
int dim,
int spacedim>
2243 const std::vector<unsigned int> &subdivisions,
2246 std::vector<unsigned int> compute_subdivisions = subdivisions;
2247 if (compute_subdivisions.empty())
2249 compute_subdivisions.resize(dim, 1);
2252 Assert(compute_subdivisions.size() == dim,
2253 ExcMessage(
"One subdivision must be provided for each dimension."));
2255 for (
unsigned int i = 0; i < dim; ++i)
2257 Assert(compute_subdivisions[i] > 0,
2258 ExcInvalidRepetitions(subdivisions[i]));
2260 edges[i].
norm() > 0,
2262 "Edges in subdivided_parallelepiped() must not be degenerate."));
2270 bool twisted_data =
false;
2275 twisted_data = (edges[0][0] < 0);
2282 const double plane_normal =
2283 edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2284 twisted_data = (plane_normal < 0.0);
2293 (edges[0].
norm() * edges[1].
norm()) -
2296 "Edges in subdivided_parallelepiped() must point in"
2297 " different directions."));
2299 cross_product_3d(edges[0], edges[1]);
2313 twisted_data = (plane_normal * edges[2] < 0.0);
2322 ExcInvalidInputOrientation(
2323 "The triangulation you are trying to create will consist of cells"
2324 " with negative measures. This is usually the result of input data"
2325 " that does not define a right-handed coordinate system. The usual"
2326 " fix for this is to ensure that in 1d the given point is to the"
2327 " right of the origin (or the given edge tensor is positive), in 2d"
2328 " that the two edges (and their cross product) obey the right-hand"
2329 " rule (which may usually be done by switching the order of the"
2330 " points or edge tensors), or in 3d that the edges form a"
2331 " right-handed coordinate system (which may also be accomplished by"
2332 " switching the order of the first two points or edge tensors)."));
2335 for (
unsigned int i = 0; i < dim; ++i)
2336 for (
unsigned int j = i + 1; j < dim; ++j)
2337 Assert((edges[i] != edges[j]),
2339 "Degenerate edges of subdivided_parallelepiped encountered."));
2342 std::vector<Point<spacedim>> points;
2347 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2348 points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2352 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2353 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2354 points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2355 edges[1] / compute_subdivisions[1] * y);
2359 for (
unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2360 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2361 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2362 points.push_back(origin +
2363 edges[0] / compute_subdivisions[0] * x +
2364 edges[1] / compute_subdivisions[1] * y +
2365 edges[2] / compute_subdivisions[2] * z);
2374 for (
unsigned int i = 0; i < dim; ++i)
2375 n_cells *= compute_subdivisions[i];
2376 std::vector<CellData<dim>> cells(n_cells);
2382 for (
unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2384 cells[x].vertices[0] = x;
2385 cells[x].vertices[1] = x + 1;
2388 cells[x].material_id = 0;
2395 const unsigned int n_dy = compute_subdivisions[1];
2396 const unsigned int n_dx = compute_subdivisions[0];
2398 for (
unsigned int y = 0; y < n_dy; ++y)
2399 for (
unsigned int x = 0; x < n_dx; ++x)
2401 const unsigned int c = y * n_dx + x;
2402 cells[c].vertices[0] = y * (n_dx + 1) + x;
2403 cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2404 cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2405 cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2408 cells[c].material_id = 0;
2416 const unsigned int n_dz = compute_subdivisions[2];
2417 const unsigned int n_dy = compute_subdivisions[1];
2418 const unsigned int n_dx = compute_subdivisions[0];
2420 for (
unsigned int z = 0; z < n_dz; ++z)
2421 for (
unsigned int y = 0; y < n_dy; ++y)
2422 for (
unsigned int x = 0; x < n_dx; ++x)
2424 const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2426 cells[c].vertices[0] =
2427 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2428 cells[c].vertices[1] =
2429 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2430 cells[c].vertices[2] =
2431 z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2432 cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2433 (y + 1) * (n_dx + 1) + x + 1;
2434 cells[c].vertices[4] =
2435 (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2436 cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2437 y * (n_dx + 1) + x + 1;
2438 cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2439 (y + 1) * (n_dx + 1) + x;
2440 cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2441 (y + 1) * (n_dx + 1) + x + 1;
2444 cells[c].material_id = 0;
2465 for (; cell != endc; ++cell)
2467 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
2469 if (cell->face(face)->at_boundary())
2470 cell->face(face)->set_boundary_id(face);
2477 template <
int dim,
int spacedim>
2480 const unsigned int repetitions,
2485 Assert(repetitions >= 1, ExcInvalidRepetitions(repetitions));
2487 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2490 for (
unsigned int i = 0; i < dim; ++i)
2496 std::vector<unsigned int> reps(dim, repetitions);
2502 template <
int dim,
int spacedim>
2505 const std::vector<unsigned int> &repetitions,
2510 Assert(repetitions.size() == dim, ExcInvalidRepetitionsDimension(dim));
2516 for (
unsigned int i = 0; i < dim; ++i)
2523 std::array<Point<spacedim>, dim> delta;
2524 for (
unsigned int i = 0; i < dim; ++i)
2526 Assert(repetitions[i] >= 1, ExcInvalidRepetitions(repetitions[i]));
2528 delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2532 "The first dim entries of coordinates of p1 and p2 need to be different."));
2536 std::vector<Point<spacedim>> points;
2540 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2541 points.push_back(p1 + x * delta[0]);
2545 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2546 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2547 points.push_back(p1 + x * delta[0] + y * delta[1]);
2551 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
2552 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2553 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2554 points.push_back(p1 + x * delta[0] + y * delta[1] +
2563 std::vector<CellData<dim>> cells;
2568 cells.resize(repetitions[0]);
2569 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2571 cells[x].vertices[0] = x;
2572 cells[x].vertices[1] = x + 1;
2573 cells[x].material_id = 0;
2580 cells.resize(repetitions[1] * repetitions[0]);
2581 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2582 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2584 const unsigned int c = x + y * repetitions[0];
2585 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2586 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2587 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2588 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2589 cells[c].material_id = 0;
2596 const unsigned int n_x = (repetitions[0] + 1);
2597 const unsigned int n_xy =
2598 (repetitions[0] + 1) * (repetitions[1] + 1);
2600 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2601 for (
unsigned int z = 0; z < repetitions[2]; ++z)
2602 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2603 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2605 const unsigned int c = x + y * repetitions[0] +
2606 z * repetitions[0] * repetitions[1];
2607 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2608 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2609 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2610 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2611 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2612 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2613 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2614 cells[c].vertices[7] =
2615 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2616 cells[c].material_id = 0;
2637 double epsilon = std::numeric_limits<double>::max();
2638 for (
unsigned int i = 0; i < dim; ++i)
2639 epsilon =
std::min(epsilon, 0.01 * delta[i][i]);
2642 "The distance between corner points must be positive."));
2646 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2655 const std::vector<std::vector<double>> &step_sz,
2660 Assert(step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
2670 std::vector<std::vector<double>> step_sizes(step_sz);
2672 for (
unsigned int i = 0; i < dim; ++i)
2676 std::swap(p1[i], p2[i]);
2677 std::reverse(step_sizes[i].
begin(), step_sizes[i].
end());
2682 for (
unsigned int j = 0; j < step_sizes.at(i).size(); ++j)
2683 x += step_sizes[i][j];
2684 Assert(std::fabs(x - (p2[i] - p1[i])) <= 1e-12 * std::fabs(x),
2686 "The sequence of step sizes in coordinate direction " +
2688 " must be equal to the distance of the two given "
2689 "points in this coordinate direction."));
2696 std::vector<Point<dim>> points;
2702 for (
unsigned int i = 0;; ++i)
2711 if (i == step_sizes[0].size())
2714 x += step_sizes[0][i];
2722 for (
unsigned int j = 0;; ++j)
2725 for (
unsigned int i = 0;; ++i)
2727 points.push_back(
Point<dim>(p1[0] + x, p1[1] + y));
2728 if (i == step_sizes[0].size())
2731 x += step_sizes[0][i];
2734 if (j == step_sizes[1].size())
2737 y += step_sizes[1][j];
2744 for (
unsigned int k = 0;; ++k)
2747 for (
unsigned int j = 0;; ++j)
2750 for (
unsigned int i = 0;; ++i)
2753 Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2754 if (i == step_sizes[0].size())
2757 x += step_sizes[0][i];
2760 if (j == step_sizes[1].size())
2763 y += step_sizes[1][j];
2766 if (k == step_sizes[2].size())
2769 z += step_sizes[2][k];
2780 std::vector<CellData<dim>> cells;
2785 cells.resize(step_sizes[0].size());
2786 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2788 cells[x].vertices[0] = x;
2789 cells[x].vertices[1] = x + 1;
2790 cells[x].material_id = 0;
2797 cells.resize(step_sizes[1].size() * step_sizes[0].size());
2798 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2799 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2801 const unsigned int c = x + y * step_sizes[0].size();
2802 cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2803 cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2804 cells[c].vertices[2] =
2805 (y + 1) * (step_sizes[0].size() + 1) + x;
2806 cells[c].vertices[3] =
2807 (y + 1) * (step_sizes[0].size() + 1) + x + 1;
2808 cells[c].material_id = 0;
2815 const unsigned int n_x = (step_sizes[0].size() + 1);
2816 const unsigned int n_xy =
2817 (step_sizes[0].size() + 1) * (step_sizes[1].size() + 1);
2819 cells.resize(step_sizes[2].size() * step_sizes[1].size() *
2820 step_sizes[0].size());
2821 for (
unsigned int z = 0; z < step_sizes[2].size(); ++z)
2822 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2823 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2825 const unsigned int c =
2826 x + y * step_sizes[0].size() +
2827 z * step_sizes[0].size() * step_sizes[1].size();
2828 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2829 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2830 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2831 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2832 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2833 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2834 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2835 cells[c].vertices[7] =
2836 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2837 cells[c].material_id = 0;
2859 *std::min_element(step_sizes[0].
begin(), step_sizes[0].
end());
2860 for (
unsigned int i = 1; i < dim; ++i)
2862 *std::min_element(step_sizes[i].
begin(),
2863 step_sizes[i].
end()));
2864 const double epsilon = 0.01 * min_size;
2868 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
2877 const std::vector<std::vector<double>> &spacing,
2882 Assert(spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
2886 Assert(spacing[0].size() == n_cells, ExcInvalidRepetitionsDimension(1));
2888 double delta = std::numeric_limits<double>::max();
2889 for (
unsigned int i = 0; i <
n_cells; ++i)
2891 Assert(spacing[0][i] >= 0, ExcInvalidRepetitions(-1));
2892 delta =
std::min(delta, spacing[0][i]);
2896 std::vector<Point<1>> points;
2898 for (
unsigned int x = 0; x <=
n_cells; ++x)
2900 points.emplace_back(ax);
2902 ax += spacing[0][x];
2905 unsigned int n_val_cells = 0;
2906 for (
unsigned int i = 0; i <
n_cells; ++i)
2910 std::vector<CellData<1>> cells(n_val_cells);
2911 unsigned int id = 0;
2912 for (
unsigned int x = 0; x <
n_cells; ++x)
2915 cells[id].vertices[0] = x;
2916 cells[id].vertices[1] = x + 1;
2935 const std::vector<std::vector<double>> &spacing,
2940 Assert(spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
2942 std::vector<unsigned int> repetitions(2);
2943 double delta = std::numeric_limits<double>::max();
2944 for (
unsigned int i = 0; i < 2; ++i)
2946 repetitions[i] = spacing[i].size();
2947 for (
unsigned int j = 0; j < repetitions[i]; ++j)
2949 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
2950 delta =
std::min(delta, spacing[i][j]);
2953 ExcInvalidRepetitionsDimension(i));
2957 std::vector<Point<2>> points;
2959 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2962 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2964 points.emplace_back(ax, ay);
2965 if (x < repetitions[0])
2966 ax += spacing[0][x];
2968 if (y < repetitions[1])
2969 ay += spacing[1][y];
2973 unsigned int n_val_cells = 0;
2974 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
2975 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
2979 std::vector<CellData<2>> cells(n_val_cells);
2980 unsigned int id = 0;
2981 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2982 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2985 cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
2986 cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2987 cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2988 cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3002 double eps = 0.01 * delta;
3004 for (; cell != endc; ++cell)
3006 Point<2> cell_center = cell->center();
3007 for (
const unsigned int f :
GeometryInfo<2>::face_indices())
3010 Point<2> face_center = cell->face(f)->center();
3011 for (
unsigned int i = 0; i < 2; ++i)
3013 if (face_center[i] < cell_center[i] - eps)
3014 cell->face(f)->set_boundary_id(i * 2);
3015 if (face_center[i] > cell_center[i] + eps)
3016 cell->face(f)->set_boundary_id(i * 2 + 1);
3027 const std::vector<std::vector<double>> &spacing,
3032 const unsigned int dim = 3;
3034 Assert(spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
3036 std::array<unsigned int, dim> repetitions;
3037 double delta = std::numeric_limits<double>::max();
3038 for (
unsigned int i = 0; i < dim; ++i)
3040 repetitions[i] = spacing[i].size();
3041 for (
unsigned int j = 0; j < repetitions[i]; ++j)
3043 Assert(spacing[i][j] >= 0, ExcInvalidRepetitions(-1));
3044 delta =
std::min(delta, spacing[i][j]);
3047 ExcInvalidRepetitionsDimension(i));
3051 std::vector<Point<dim>> points;
3053 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3056 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3059 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3061 points.emplace_back(ax, ay, az);
3062 if (x < repetitions[0])
3063 ax += spacing[0][x];
3065 if (y < repetitions[1])
3066 ay += spacing[1][y];
3068 if (z < repetitions[2])
3069 az += spacing[2][z];
3073 unsigned int n_val_cells = 0;
3074 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
3075 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
3076 for (
unsigned int k = 0; k <
material_id.size(2); ++k)
3080 std::vector<CellData<dim>> cells(n_val_cells);
3081 unsigned int id = 0;
3082 const unsigned int n_x = (repetitions[0] + 1);
3083 const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3084 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3085 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3086 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3089 cells[id].vertices[0] = z * n_xy + y * n_x + x;
3090 cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3091 cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3092 cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3093 cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3094 cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3095 cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3096 cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3110 double eps = 0.01 * delta;
3113 for (; cell != endc; ++cell)
3119 Point<dim> face_center = cell->face(f)->center();
3120 for (
unsigned int i = 0; i < dim; ++i)
3122 if (face_center[i] < cell_center[i] - eps)
3123 cell->face(f)->set_boundary_id(i * 2);
3124 if (face_center[i] > cell_center[i] + eps)
3125 cell->face(f)->set_boundary_id(i * 2 + 1);
3132 template <
int dim,
int spacedim>
3135 const std::vector<unsigned int> &holes)
3144 for (
unsigned int d = 0;
d < dim; ++
d)
3151 std::array<Point<spacedim>, dim> delta;
3152 std::array<unsigned int, dim> repetitions;
3153 for (
unsigned int i = 0; i < dim; ++i)
3156 ExcMessage(
"At least one hole needed in each direction"));
3157 repetitions[i] = 2 * holes[i] + 1;
3158 delta[i][i] = (p2[i] - p1[i]);
3163 std::vector<Point<spacedim>> points;
3167 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3168 points.push_back(p1 + x * delta[0]);
3172 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3173 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3174 points.push_back(p1 + x * delta[0] + y * delta[1]);
3178 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3179 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3180 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3181 points.push_back(p1 + x * delta[0] + y * delta[1] +
3191 std::vector<CellData<dim>> cells;
3196 cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3198 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3199 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3201 if ((x % 2 == 1) && (y % 2 == 1))
3204 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
3205 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3206 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3207 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3208 cells[c].material_id = 0;
3216 const unsigned int n_x = (repetitions[0] + 1);
3217 const unsigned int n_xy =
3218 (repetitions[0] + 1) * (repetitions[1] + 1);
3220 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3223 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3224 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3225 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3228 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3229 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3230 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3231 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3232 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3233 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3234 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3235 cells[c].vertices[7] =
3236 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3237 cells[c].material_id = 0;
3265 const unsigned int ,
3277 const unsigned int ,
3289 bool inline point_in_2d_box(
const Point<2> &p,
3291 const double radius)
3293 return (
std::abs(p[0] - c[0]) < radius) &&
3302 template <
int dim,
int spacedim>
3306 double length = std::numeric_limits<double>::max();
3307 for (
const auto &cell :
triangulation.active_cell_iterators())
3308 for (unsigned
int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3309 length =
std::min(length, cell->line(n)->diameter());
3319 const double inner_radius,
3320 const double outer_radius,
3321 const double pad_bottom,
3322 const double pad_top,
3323 const double pad_left,
3324 const double pad_right,
3329 const unsigned int ,
3332 const bool with_padding =
3333 pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3343 double length = std::numeric_limits<double>::max();
3344 for (
const auto &cell :
tria.active_cell_iterators())
3345 for (unsigned
int n = 0; n < cell->n_lines(); ++n)
3346 length =
std::min(length, cell->line(n)->diameter());
3361 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3362 cell->set_manifold_id(tfi_manifold_id);
3364 const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
3365 const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
3371 auto add_sizes = [](std::vector<double> &step_sizes,
3372 const double padding,
3373 const double h) ->
void {
3376 const auto rounded =
3377 static_cast<unsigned int>(std::round(padding / h));
3380 const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3381 for (
unsigned int i = 0; i < num; ++i)
3382 step_sizes.push_back(padding / num);
3385 std::vector<std::vector<double>> step_sizes(2);
3388 add_sizes(step_sizes[0], pad_left, outer_radius);
3390 step_sizes[0].push_back(outer_radius);
3391 step_sizes[0].push_back(outer_radius);
3393 add_sizes(step_sizes[0], pad_right, outer_radius);
3396 add_sizes(step_sizes[1], pad_bottom, outer_radius);
3398 step_sizes[1].push_back(outer_radius);
3399 step_sizes[1].push_back(outer_radius);
3401 add_sizes(step_sizes[1], pad_top, outer_radius);
3406 bulk_tria, step_sizes, bl, tr,
colorize);
3409 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3410 for (
const auto &cell : bulk_tria.active_cell_iterators())
3412 cells_to_remove.
insert(cell);
3416 bulk_tria, cells_to_remove, tria_without_cylinder);
3418 const double tolerance =
3419 std::min(min_line_length(tria_without_cylinder),
3420 min_line_length(cylinder_tria)) /
3430 for (
const auto &cell :
tria.active_cell_iterators())
3434 if (cell->manifold_id() == tfi_manifold_id)
3436 for (
const unsigned int face_n : cell->face_indices())
3438 const auto &face = cell->face(face_n);
3439 if (face->at_boundary() &&
3440 internal::point_in_2d_box(face->center(),
3442 outer_radius * (1. - 1e-12)))
3443 face->set_manifold_id(polar_manifold_id);
3445 face->set_manifold_id(tfi_manifold_id);
3456 static constexpr double tol =
3457 std::numeric_limits<double>::epsilon() * 10000;
3459 for (
const auto &cell :
tria.active_cell_iterators())
3460 for (const unsigned
int face_n : cell->face_indices())
3462 const auto face = cell->face(face_n);
3463 if (face->at_boundary())
3468 face->set_boundary_id(0);
3471 face->set_boundary_id(1);
3474 face->set_boundary_id(2);
3477 face->set_boundary_id(3);
3481 Assert(cell->manifold_id() == tfi_manifold_id,
3483 face->set_boundary_id(4);
3504 const double inner_radius,
3505 const double outer_radius,
3506 const double pad_bottom,
3507 const double pad_top,
3508 const double pad_left,
3509 const double pad_right,
3514 const unsigned int n_slices,
3525 Point<2>(new_center[0], new_center[1]),
3557 const double shell_region_width,
3558 const unsigned int n_shells,
3559 const double skewness,
3562 Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3563 ExcMessage(
"The width of the shell region must be less than 0.05 "
3564 "(and preferably close to 0.03)"));
3598 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3600 for (
const auto &cell : bulk_tria.active_cell_iterators())
3602 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3603 cells_to_remove.insert(cell);
3607 for (
const unsigned int vertex_n :
3609 if (cell->vertex(vertex_n) ==
Point<2>())
3613 cylinder_triangulation_offset =
3614 2.0 * (cell->vertex(3) -
Point<2>());
3621 bulk_tria, cells_to_remove, tria_without_cylinder);
3628 0.05 + shell_region_width,
3633 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3636 if (
std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1e-10)
3637 cell->vertex(vertex_n)[0] = -0.1;
3638 else if (
std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1e-10)
3639 cell->vertex(vertex_n)[0] = 0.1;
3643 for (
const auto &cell : cylinder_tria.active_cell_iterators())
3645 cell->set_manifold_id(tfi_manifold_id);
3646 for (
const unsigned int face_n :
GeometryInfo<2>::face_indices())
3647 if (!cell->face(face_n)->at_boundary())
3648 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3650 if (0.0 < shell_region_width)
3653 ExcMessage(
"If the shell region has positive width then "
3654 "there must be at least one shell."));
3659 0.05 + shell_region_width,
3666 const double vertex_tolerance =
3667 std::min(internal::minimal_vertex_distance(shell_tria),
3668 internal::minimal_vertex_distance(cylinder_tria)) *
3674 shell_tria, cylinder_tria, temp, vertex_tolerance,
true);
3675 cylinder_tria = std::move(temp);
3681 const double vertex_tolerance =
3682 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3683 internal::minimal_vertex_distance(cylinder_tria)) /
3686 tria_without_cylinder, cylinder_tria, tria, vertex_tolerance,
true);
3699 const double shift =
3700 std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3701 for (
const auto &cell :
tria.active_cell_iterators())
3703 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3704 cell->vertex(v) =
Point<2>(0.2 -
shift, 0.205);
3705 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1e-10)
3706 cell->vertex(v) =
Point<2>(0.2 + shift, 0.205);
3707 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1e-10)
3708 cell->vertex(v) =
Point<2>(0.2, 0.2 - shift);
3709 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1e-10)
3710 cell->vertex(v) =
Point<2>(0.2, 0.2 + shift);
3715 for (
const auto &cell :
tria.active_cell_iterators())
3717 cell->set_all_manifold_ids(polar_manifold_id);
3721 for (
const auto &cell :
tria.active_cell_iterators())
3728 std::vector<Point<2> *> cylinder_pointers;
3729 for (
const auto &face :
tria.active_face_iterators())
3732 cylinder_pointers.push_back(&face->vertex(0));
3733 cylinder_pointers.push_back(&face->vertex(1));
3736 std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3737 cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3738 cylinder_pointers.end()),
3739 cylinder_pointers.end());
3743 for (
const Point<2> *
const ptr : cylinder_pointers)
3744 center += *ptr / double(cylinder_pointers.size());
3747 for (
Point<2> *
const ptr : cylinder_pointers)
3760 for (
const auto &face :
tria.active_face_iterators())
3761 if (face->at_boundary())
3766 face->set_boundary_id(0);
3769 face->set_boundary_id(1);
3771 else if (face->manifold_id() == polar_manifold_id)
3772 face->set_boundary_id(2);
3779 face->set_boundary_id(3);
3789 const double shell_region_width,
3790 const unsigned int n_shells,
3791 const double skewness,
3796 tria_2, shell_region_width, n_shells, skewness,
colorize);
3821 for (
const auto &face :
tria.active_face_iterators())
3823 face->set_boundary_id(3);
3828 template <
int dim,
int spacedim>
3831 const std::vector<unsigned int> &sizes,
3841 for (
unsigned int d = 0;
d < dim; ++
d)
3844 std::vector<Point<spacedim>> points;
3846 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
3849 std::vector<CellData<dim>> cells(n_cells);
3854 for (
unsigned int d = 0;
d < dim; ++
d)
3855 p[d] = 0.5 * dimensions[d] *
3858 points.push_back(p);
3859 cells[0].vertices[i] = i;
3861 cells[0].material_id = 0;
3866 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
3872 for (
unsigned int j = 0; j < sizes[face]; ++j, ++
cell_index)
3874 const unsigned int last_cell = (j == 0) ? 0U : (
cell_index - 1);
3876 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
3879 const unsigned int cellv =
3881 const unsigned int ocellv =
3885 cells[last_cell].vertices[cellv];
3888 cells[
cell_index].vertices[cellv] = points.size();
3893 points.push_back(p);
4051 const double thickness,
4055 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4057 std::vector<Point<2>>
vertices(16);
4059 coords[0] = left - thickness;
4062 coords[3] = right + thickness;
4065 for (
const double y : coords)
4066 for (const double x : coords)
4071 std::vector<CellData<2>> cells(9);
4073 for (
unsigned int i0 = 0; i0 < 3; ++i0)
4074 for (
unsigned int i1 = 0; i1 < 3; ++i1)
4076 cells[k].vertices[0] = i1 + 4 * i0;
4077 cells[k].vertices[1] = i1 + 4 * i0 + 1;
4078 cells[k].vertices[2] = i1 + 4 * i0 + 4;
4079 cells[k].vertices[3] = i1 + 4 * i0 + 5;
4081 cells[k].material_id = materials[k];
4099 const double rl2 = (right + left) / 2;
4110 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4115 for (
unsigned int i = 0; i < 4; ++i)
4117 for (
unsigned int j = 0; j < 4; ++j)
4118 cells[i].
vertices[j] = cell_vertices[i][j];
4119 cells[i].material_id = 0;
4129 cell->face(1)->set_boundary_id(1);
4131 cell->face(0)->set_boundary_id(2);
4140 const double radius_0,
4141 const double radius_1,
4142 const double half_length)
4146 vertices_tmp[0] =
Point<2>(-half_length, -radius_0);
4147 vertices_tmp[1] =
Point<2>(half_length, -radius_1);
4148 vertices_tmp[2] =
Point<2>(-half_length, radius_0);
4149 vertices_tmp[3] =
Point<2>(half_length, radius_1);
4151 const std::vector<Point<2>>
vertices(std::begin(vertices_tmp),
4152 std::end(vertices_tmp));
4156 cell_vertices[0][i] = i;
4161 cells[0].
vertices[i] = cell_vertices[0][i];
4163 cells[0].material_id = 0;
4168 cell->face(0)->set_boundary_id(1);
4169 cell->face(1)->set_boundary_id(2);
4171 for (
unsigned int i = 2; i < 4; ++i)
4172 cell->face(i)->set_boundary_id(0);
4189 Point<2>((a + b) / 2, (a + b) / 2),
4193 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4197 for (
unsigned int i = 0; i < 3; ++i)
4199 for (
unsigned int j = 0; j < 4; ++j)
4200 cells[i].
vertices[j] = cell_vertices[i][j];
4201 cells[i].material_id = 0;
4213 cell->face(0)->set_boundary_id(0);
4214 cell->face(2)->set_boundary_id(1);
4217 cell->face(1)->set_boundary_id(2);
4218 cell->face(2)->set_boundary_id(1);
4219 cell->face(3)->set_boundary_id(3);
4222 cell->face(0)->set_boundary_id(0);
4223 cell->face(1)->set_boundary_id(4);
4224 cell->face(3)->set_boundary_id(5);
4230 template <
int dim,
int spacedim>
4233 const std::vector<unsigned int> &repetitions,
4236 const std::vector<int> &n_cells_to_remove)
4242 for (
unsigned int d = 0;
d < dim; ++
d)
4244 Assert(std::fabs(n_cells_to_remove[d]) <= repetitions[d],
4245 ExcMessage(
"Attempting to cut away too many cells."));
4255 std::array<double, dim> h;
4257 for (
unsigned int d = 0;
d < dim; ++
d)
4260 h[
d] = (top_right[
d] - bottom_left[
d]) / repetitions[d];
4262 if (n_cells_to_remove[d] >= 0)
4266 h[
d] * std::fabs(n_cells_to_remove[d]) + bottom_left[
d];
4271 cut_step[
d] = top_right[
d] - h[
d] * std::fabs(n_cells_to_remove[d]);
4277 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4279 for (
const auto &cell : rectangle.active_cell_iterators())
4281 bool remove_cell =
true;
4282 for (
unsigned int d = 0;
d < dim && remove_cell; ++
d)
4283 if ((n_cells_to_remove[d] > 0 && cell->center()[d] >= cut_step[d]) ||
4284 (n_cells_to_remove[d] < 0 && cell->center()[d] <= cut_step[d]))
4285 remove_cell =
false;
4287 cells_to_remove.insert(cell);
4302 const double radius,
4303 const bool internal_manifolds)
4308 const double a = 1. / (1 +
std::sqrt(2.0));
4321 for (
unsigned int i = 0; i < 5; ++i)
4323 for (
unsigned int j = 0; j < 4; ++j)
4324 cells[i].
vertices[j] = circle_cell_vertices[i][j];
4325 cells[i].material_id = 0;
4335 if (internal_manifolds)
4347 const double inner_radius,
4348 const double outer_radius,
4349 const unsigned int n_cells,
4352 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4366 const unsigned int N =
4367 (
n_cells == 0 ?
static_cast<unsigned int>(
4368 std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
4369 (outer_radius - inner_radius))) :
4378 std::vector<Point<2>>
vertices(2 * N);
4379 for (
unsigned int i = 0; i < N; ++i)
4392 for (
unsigned int i = 0; i < N; ++i)
4394 cells[i].vertices[0] = i;
4395 cells[i].vertices[1] = (i + 1) % N;
4396 cells[i].vertices[2] = N + i;
4397 cells[i].vertices[3] = N + ((i + 1) % N);
4399 cells[i].material_id = 0;
4405 colorize_hyper_shell(tria,
center, inner_radius, outer_radius);
4418 const double inner_radius,
4419 const double outer_radius,
4420 const unsigned int n_cells)
4423 tria, outer_center, inner_radius, outer_radius, n_cells,
true);
4427 outer_radius - inner_radius > outer_center.
distance(inner_center),
4429 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4433 std::set<Point<dim> *> vertices_to_move;
4435 for (
const auto &face :
tria.active_face_iterators())
4437 for (unsigned
int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4438 vertices_to_move.insert(&face->vertex(v));
4440 const auto shift = inner_center - outer_center;
4441 for (
const auto &p : vertices_to_move)
4472 const double radius,
4473 const double half_length)
4475 Point<2> p1(-half_length, -radius);
4484 switch (f->boundary_id())
4487 f->set_boundary_id(1);
4490 f->set_boundary_id(2);
4493 f->set_boundary_id(0);
4531 const double radius)
4533 const unsigned int dim = 2;
4547 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4551 for (
unsigned int i = 0; i < 3; ++i)
4553 for (
unsigned int j = 0; j < 4; ++j)
4554 cells[i].
vertices[j] = cell_vertices[i][j];
4555 cells[i].material_id = 0;
4570 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
4572 if (cell->face(i)->boundary_id() ==
4578 if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius ||
4579 cell->face(i)->center()[1] < p[1] + 1.e-5 * radius)
4581 cell->face(i)->set_boundary_id(1);
4595 const double radius)
4600 const double a = 1. / (1 +
std::sqrt(2.0));
4611 const int cell_vertices[5][4] = {{0, 1, 2, 3},
4618 for (
unsigned int i = 0; i < 4; ++i)
4620 for (
unsigned int j = 0; j < 4; ++j)
4621 cells[i].
vertices[j] = cell_vertices[i][j];
4622 cells[i].material_id = 0;
4637 for (
const unsigned int i :
GeometryInfo<2>::face_indices())
4639 if (cell->face(i)->boundary_id() ==
4644 if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius)
4646 cell->face(i)->set_boundary_id(1);
4662 const double inner_radius,
4663 const double outer_radius,
4664 const unsigned int n_cells,
4667 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4680 const unsigned int N =
4681 (
n_cells == 0 ?
static_cast<unsigned int>(
4682 std::ceil((pi * (outer_radius + inner_radius) / 2) /
4683 (outer_radius - inner_radius))) :
4692 std::vector<Point<2>>
vertices(2 * (N + 1));
4693 for (
unsigned int i = 0; i <= N; ++i)
4702 Point<2>(((i == 0) || (i == N) ? 0 :
std::
cos(pi * i / N - pi / 2)),
4703 std::
sin(pi * i / N - pi / 2)) *
4714 for (
unsigned int i = 0; i < N; ++i)
4716 cells[i].vertices[0] = i;
4717 cells[i].vertices[1] = (i + 1) % (N + 1);
4718 cells[i].vertices[2] = N + 1 + i;
4719 cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4721 cells[i].material_id = 0;
4729 for (; cell !=
tria.
end(); ++cell)
4731 cell->face(2)->set_boundary_id(1);
4733 tria.
begin()->face(0)->set_boundary_id(3);
4735 tria.
last()->face(1)->set_boundary_id(2);
4746 const double inner_radius,
4747 const double outer_radius,
4748 const unsigned int n_cells,
4751 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4764 const unsigned int N =
4765 (
n_cells == 0 ?
static_cast<unsigned int>(
4766 std::ceil((pi * (outer_radius + inner_radius) / 4) /
4767 (outer_radius - inner_radius))) :
4776 std::vector<Point<2>>
vertices(2 * (N + 1));
4777 for (
unsigned int i = 0; i <= N; ++i)
4785 std::
sin(pi * i / N / 2)) *
4796 for (
unsigned int i = 0; i < N; ++i)
4798 cells[i].vertices[0] = i;
4799 cells[i].vertices[1] = (i + 1) % (N + 1);
4800 cells[i].vertices[2] = N + 1 + i;
4801 cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4803 cells[i].material_id = 0;
4811 for (; cell !=
tria.
end(); ++cell)
4813 cell->face(2)->set_boundary_id(1);
4815 tria.
begin()->face(0)->set_boundary_id(3);
4817 tria.
last()->face(1)->set_boundary_id(2);
4834 const double rl2 = (right + left) / 2;
4835 const double len = (right - left) / 2.;
4848 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4849 {9, 4, 2, 5, 19, 14, 12, 15},
4850 {3, 2, 7, 6, 13, 12, 17, 16},
4851 {2, 5, 6, 8, 12, 15, 16, 18}};
4853 for (
unsigned int i = 0; i < 4; ++i)
4855 for (
unsigned int j = 0; j < 8; ++j)
4856 cells[i].
vertices[j] = cell_vertices[i][j];
4857 cells[i].material_id = 0;
4867 cell->face(1)->set_boundary_id(1);
4869 cell->face(0)->set_boundary_id(2);
4881 const double thickness,
4885 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4887 std::vector<Point<3>>
vertices(64);
4889 coords[0] = left - thickness;
4892 coords[3] = right + thickness;
4895 for (
const double z : coords)
4896 for (const double y : coords)
4897 for (const double x : coords)
4901 24, 26, 5, 4, 6, 1, 0,
4902 2, 9, 8, 10, 37, 36, 38,
4903 33, 32, 34, 41, 40, 42};
4905 std::vector<CellData<3>> cells(27);
4907 for (
unsigned int z = 0; z < 3; ++z)
4908 for (
unsigned int y = 0; y < 3; ++y)
4909 for (
unsigned int x = 0; x < 3; ++x)
4911 cells[k].vertices[0] = x + 4 * y + 16 * z;
4912 cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
4913 cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
4914 cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
4915 cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
4916 cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
4917 cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
4918 cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
4920 cells[k].material_id = materials[k];
4933 const double radius_0,
4934 const double radius_1,
4935 const double half_length)
4938 ExcMessage(
"The output triangulation object needs to be empty."));
4943 const auto n_slices = 1 +
static_cast<unsigned int>(std::ceil(
4944 half_length /
std::max(radius_0, radius_1)));
4957 auto shift_radii = [=](
const Point<3> &
p) {
4958 const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
4959 const double factor = slope * (
p[0] - -half_length) + 1.0;
4960 return Point<3>(p[0], factor * p[1], factor * p[2]);
4966 for (
const auto &face :
triangulation.active_face_iterators())
4967 if (face->at_boundary())
4969 if (
std::abs(face->center()[0] - -half_length) < 1e-8 * half_length)
4970 face->set_boundary_id(1);
4971 else if (
std::abs(face->center()[0] - half_length) <
4973 face->set_boundary_id(2);
4975 face->set_all_manifold_ids(0);
4998 Point<3>((a + b) / 2, a, (a + b) / 2),
5005 Point<3>((a + b) / 2, (a + b) / 2, a),
5007 Point<3>(a, (a + b) / 2, (a + b) / 2),
5008 Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
5009 Point<3>(b, (a + b) / 2, (a + b) / 2),
5011 Point<3>((a + b) / 2, (a + b) / 2, b),
5019 Point<3>((a + b) / 2, b, (a + b) / 2),
5023 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5024 {1, 2, 10, 11, 4, 5, 13, 14},
5025 {3, 4, 12, 13, 6, 7, 15, 16},
5026 {4, 5, 13, 14, 7, 8, 16, 17},
5027 {9, 10, 18, 19, 12, 13, 21, 22},
5028 {10, 11, 19, 20, 13, 14, 22, 23},
5029 {12, 13, 21, 22, 15, 16, 24, 25}};
5033 for (
unsigned int i = 0; i < 7; ++i)
5035 for (
unsigned int j = 0; j < 8; ++j)
5036 cells[i].
vertices[j] = cell_vertices[i][j];
5037 cells[i].material_id = 0;
5058 const double radius,
5059 const bool internal_manifold)
5065 const unsigned int n_vertices = 16;
5091 const unsigned int n_cells = 7;
5092 const int cell_vertices[
n_cells][8] = {
5093 {0, 1, 4, 5, 3, 2, 7, 6},
5094 {8, 9, 12, 13, 0, 1, 4, 5},
5095 {9, 13, 1, 5, 10, 14, 2, 6},
5096 {11, 10, 3, 2, 15, 14, 7, 6},
5097 {8, 0, 12, 4, 11, 3, 15, 7},
5098 {8, 9, 0, 1, 11, 10, 3, 2},
5099 {12, 4, 13, 5, 15, 7, 14, 6}};
5101 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5103 for (
unsigned int i = 0; i <
n_cells; ++i)
5106 cells[i].
vertices[j] = cell_vertices[i][j];
5107 cells[i].material_id = 0;
5117 if (internal_manifold)
5127 const unsigned int n_rotate_middle_square)
5130 ExcMessage(
"The number of rotation by pi/2 of the right square "
5131 "must be in the half-open range [0,4)."));
5133 constexpr unsigned int dim = 2;
5135 const unsigned int n_cells = 5;
5136 std::vector<CellData<dim>> cells(n_cells);
5154 unsigned int cell_vertices[
n_cells][4] = {{0, 1, 2, 3},
5160 switch (n_rotate_middle_square)
5164 cell_vertices[1][0] = 4;
5165 cell_vertices[1][1] = 5;
5166 cell_vertices[1][2] = 1;
5167 cell_vertices[1][3] = 3;
5173 cell_vertices[1][0] = 5;
5174 cell_vertices[1][1] = 3;
5175 cell_vertices[1][2] = 4;
5176 cell_vertices[1][3] = 1;
5182 cell_vertices[1][0] = 3;
5183 cell_vertices[1][1] = 1;
5184 cell_vertices[1][2] = 5;
5185 cell_vertices[1][3] = 4;
5197 for (
const unsigned int vertex_index :
5212 const bool face_orientation,
5213 const bool face_flip,
5214 const bool face_rotation,
5215 const bool manipulate_left_cube)
5217 constexpr unsigned int dim = 3;
5219 const unsigned int n_cells = 2;
5220 std::vector<CellData<dim>> cells(n_cells);
5236 unsigned int cell_vertices[
n_cells][8] = {
5237 {0, 1, 2, 3, 4, 5, 6, 7},
5238 {1, 8, 3, 9, 5, 10, 7, 11}};
5241 const unsigned int this_case = 4 *
static_cast<int>(face_orientation) +
5242 2 *
static_cast<int>(face_flip) +
5243 static_cast<int>(face_rotation);
5245 if (manipulate_left_cube)
5251 cell_vertices[0][0] = 1;
5252 cell_vertices[0][1] = 0;
5253 cell_vertices[0][2] = 5;
5254 cell_vertices[0][3] = 4;
5255 cell_vertices[0][4] = 3;
5256 cell_vertices[0][5] = 2;
5257 cell_vertices[0][6] = 7;
5258 cell_vertices[0][7] = 6;
5264 cell_vertices[0][0] = 5;
5265 cell_vertices[0][1] = 4;
5266 cell_vertices[0][2] = 7;
5267 cell_vertices[0][3] = 6;
5268 cell_vertices[0][4] = 1;
5269 cell_vertices[0][5] = 0;
5270 cell_vertices[0][6] = 3;
5271 cell_vertices[0][7] = 2;
5277 cell_vertices[0][0] = 7;
5278 cell_vertices[0][1] = 6;
5279 cell_vertices[0][2] = 3;
5280 cell_vertices[0][3] = 2;
5281 cell_vertices[0][4] = 5;
5282 cell_vertices[0][5] = 4;
5283 cell_vertices[0][6] = 1;
5284 cell_vertices[0][7] = 0;
5289 cell_vertices[0][0] = 3;
5290 cell_vertices[0][1] = 2;
5291 cell_vertices[0][2] = 1;
5292 cell_vertices[0][3] = 0;
5293 cell_vertices[0][4] = 7;
5294 cell_vertices[0][5] = 6;
5295 cell_vertices[0][6] = 5;
5296 cell_vertices[0][7] = 4;
5302 cell_vertices[0][0] = 0;
5303 cell_vertices[0][1] = 1;
5304 cell_vertices[0][2] = 2;
5305 cell_vertices[0][3] = 3;
5306 cell_vertices[0][4] = 4;
5307 cell_vertices[0][5] = 5;
5308 cell_vertices[0][6] = 6;
5309 cell_vertices[0][7] = 7;
5315 cell_vertices[0][0] = 2;
5316 cell_vertices[0][1] = 3;
5317 cell_vertices[0][2] = 6;
5318 cell_vertices[0][3] = 7;
5319 cell_vertices[0][4] = 0;
5320 cell_vertices[0][5] = 1;
5321 cell_vertices[0][6] = 4;
5322 cell_vertices[0][7] = 5;
5328 cell_vertices[0][0] = 6;
5329 cell_vertices[0][1] = 7;
5330 cell_vertices[0][2] = 4;
5331 cell_vertices[0][3] = 5;
5332 cell_vertices[0][4] = 2;
5333 cell_vertices[0][5] = 3;
5334 cell_vertices[0][6] = 0;
5335 cell_vertices[0][7] = 1;
5341 cell_vertices[0][0] = 4;
5342 cell_vertices[0][1] = 5;
5343 cell_vertices[0][2] = 0;
5344 cell_vertices[0][3] = 1;
5345 cell_vertices[0][4] = 6;
5346 cell_vertices[0][5] = 7;
5347 cell_vertices[0][6] = 2;
5348 cell_vertices[0][7] = 3;
5359 cell_vertices[1][0] = 8;
5360 cell_vertices[1][1] = 1;
5361 cell_vertices[1][2] = 10;
5362 cell_vertices[1][3] = 5;
5363 cell_vertices[1][4] = 9;
5364 cell_vertices[1][5] = 3;
5365 cell_vertices[1][6] = 11;
5366 cell_vertices[1][7] = 7;
5372 cell_vertices[1][0] = 10;
5373 cell_vertices[1][1] = 5;
5374 cell_vertices[1][2] = 11;
5375 cell_vertices[1][3] = 7;
5376 cell_vertices[1][4] = 8;
5377 cell_vertices[1][5] = 1;
5378 cell_vertices[1][6] = 9;
5379 cell_vertices[1][7] = 3;
5385 cell_vertices[1][0] = 11;
5386 cell_vertices[1][1] = 7;
5387 cell_vertices[1][2] = 9;
5388 cell_vertices[1][3] = 3;
5389 cell_vertices[1][4] = 10;
5390 cell_vertices[1][5] = 5;
5391 cell_vertices[1][6] = 8;
5392 cell_vertices[1][7] = 1;
5398 cell_vertices[1][0] = 9;
5399 cell_vertices[1][1] = 3;
5400 cell_vertices[1][2] = 8;
5401 cell_vertices[1][3] = 1;
5402 cell_vertices[1][4] = 11;
5403 cell_vertices[1][5] = 7;
5404 cell_vertices[1][6] = 10;
5405 cell_vertices[1][7] = 5;
5411 cell_vertices[1][0] = 1;
5412 cell_vertices[1][1] = 8;
5413 cell_vertices[1][2] = 3;
5414 cell_vertices[1][3] = 9;
5415 cell_vertices[1][4] = 5;
5416 cell_vertices[1][5] = 10;
5417 cell_vertices[1][6] = 7;
5418 cell_vertices[1][7] = 11;
5424 cell_vertices[1][0] = 5;
5425 cell_vertices[1][1] = 10;
5426 cell_vertices[1][2] = 1;
5427 cell_vertices[1][3] = 8;
5428 cell_vertices[1][4] = 7;
5429 cell_vertices[1][5] = 11;
5430 cell_vertices[1][6] = 3;
5431 cell_vertices[1][7] = 9;
5437 cell_vertices[1][0] = 7;
5438 cell_vertices[1][1] = 11;
5439 cell_vertices[1][2] = 5;
5440 cell_vertices[1][3] = 10;
5441 cell_vertices[1][4] = 3;
5442 cell_vertices[1][5] = 9;
5443 cell_vertices[1][6] = 1;
5444 cell_vertices[1][7] = 8;
5450 cell_vertices[1][0] = 3;
5451 cell_vertices[1][1] = 9;
5452 cell_vertices[1][2] = 7;
5453 cell_vertices[1][3] = 11;
5454 cell_vertices[1][4] = 1;
5455 cell_vertices[1][5] = 8;
5456 cell_vertices[1][6] = 5;
5457 cell_vertices[1][7] = 10;
5467 for (
const unsigned int vertex_index :
5481 template <
int spacedim>
5485 const double radius)
5490 GridGenerator::extract_boundary_mesh(volume_mesh, tria, boundary_ids);
5501 const unsigned int x_subdivisions,
5502 const double radius,
5503 const double half_length)
5511 const double initial_height = -half_length;
5512 const double height_increment = 2. * half_length / x_subdivisions;
5514 for (
unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep)
5516 const double height = initial_height + height_increment * rep;
5518 vertices.emplace_back(-d, height, -d);
5519 vertices.emplace_back(d, height, -d);
5520 vertices.emplace_back(-a, height, -a);
5521 vertices.emplace_back(a, height, -a);
5522 vertices.emplace_back(-a, height, a);
5523 vertices.emplace_back(a, height, a);
5524 vertices.emplace_back(-d, height, d);
5525 vertices.emplace_back(d, height, d);
5531 const double h = vertex[1];
5532 vertex[1] = -vertex[0];
5536 std::vector<std::vector<int>> cell_vertices;
5537 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5538 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5539 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5540 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5541 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5543 for (
unsigned int rep = 1; rep < x_subdivisions; ++rep)
5545 for (
unsigned int i = 0; i < 5; ++i)
5547 std::vector<int> new_cell_vertices(8);
5548 for (
unsigned int j = 0; j < 8; ++j)
5549 new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep;
5550 cell_vertices.push_back(new_cell_vertices);
5554 unsigned int n_cells = x_subdivisions * 5;
5556 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5558 for (
unsigned int i = 0; i <
n_cells; ++i)
5560 for (
unsigned int j = 0; j < 8; ++j)
5561 cells[i].
vertices[j] = cell_vertices[i][j];
5562 cells[i].material_id = 0;
5584 const double tolerance = 1
e-5 *
std::min(radius, half_length);
5586 for (
const auto &cell :
tria.cell_iterators())
5588 if (cell->at_boundary(i))
5590 if (cell->face(i)->center()[0] > half_length - tolerance)
5592 cell->face(i)->set_boundary_id(2);
5595 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5597 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5598 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5599 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5600 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5602 cell->face(i)->line(e)->set_boundary_id(2);
5603 cell->face(i)->line(e)->set_manifold_id(
5607 else if (cell->face(i)->center()[0] < -half_length + tolerance)
5609 cell->face(i)->set_boundary_id(1);
5612 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5614 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5615 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5616 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5617 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5619 cell->face(i)->line(e)->set_boundary_id(1);
5620 cell->face(i)->line(e)->set_manifold_id(
5632 const double radius,
5633 const double half_length)
5642 const double radius)
5644 const unsigned int dim = 3;
5652 const double a = 0.528;
5653 const double b = 0.4533;
5654 const double c = 0.3752;
5671 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5672 {1, 6, 2, 4, 8, 13, 9, 11},
5673 {5, 3, 6, 4, 12, 10, 13, 11},
5674 {7, 9, 10, 11, 14, 8, 12, 13}};
5678 for (
unsigned int i = 0; i < 4; ++i)
5680 for (
unsigned int j = 0; j < 8; ++j)
5681 cells[i].
vertices[j] = cell_vertices[i][j];
5682 cells[i].material_id = 0;
5696 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
5698 if (cell->face(i)->boundary_id() ==
5703 if (cell->face(i)->center()[0] <
center[0] + 1.e-5 * radius ||
5704 cell->face(i)->center()[1] <
center[1] + 1.e-5 * radius ||
5705 cell->face(i)->center()[2] <
center[2] + 1.e-5 * radius)
5707 cell->face(i)->set_boundary_id(1);
5711 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5714 const Point<3> line_vertices[2] = {
5715 cell->face(i)->line(j)->vertex(0),
5716 cell->face(i)->line(j)->vertex(1)};
5717 if ((std::fabs(line_vertices[0].distance(
center) - radius) >
5719 (std::fabs(line_vertices[1].distance(
center) - radius) >
5722 cell->face(i)->line(j)->set_boundary_id(1);
5723 cell->face(i)->line(j)->set_manifold_id(
5741 const double radius)
5747 const double b = a / 2.0;
5748 const double c =
d / 2.0;
5750 const double hb = radius *
std::sqrt(3.0) / 4.0;
5751 const double hc = radius *
std::sqrt(3.0) / 2.0;
5773 int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5774 {0, 2, 8, 10, 6, 4, 14, 12},
5775 {2, 3, 10, 11, 4, 5, 12, 13},
5776 {1, 7, 9, 15, 3, 5, 11, 13},
5777 {6, 4, 14, 12, 7, 5, 15, 13},
5778 {8, 10, 9, 11, 14, 12, 15, 13}};
5782 for (
unsigned int i = 0; i < 6; ++i)
5784 for (
unsigned int j = 0; j < 8; ++j)
5785 cells[i].
vertices[j] = cell_vertices[i][j];
5786 cells[i].material_id = 0;
5806 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
5808 if (!cell->at_boundary(i))
5814 if (cell->face(i)->center()[0] <
center[0] + 1.e-5 * radius)
5816 cell->face(i)->set_boundary_id(1);
5818 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5821 const Point<3> line_vertices[2] = {
5822 cell->face(i)->line(j)->vertex(0),
5823 cell->face(i)->line(j)->vertex(1)};
5824 if ((std::fabs(line_vertices[0].distance(
center) - radius) >
5826 (std::fabs(line_vertices[1].distance(
center) - radius) >
5829 cell->face(i)->line(j)->set_boundary_id(1);
5830 cell->face(i)->line(j)->set_manifold_id(
5847 const double radius)
5860 for (
unsigned int round = 0; round < dim; ++round)
5865 std::vector<Point<dim>> new_points(tria_copy.
n_vertices());
5867 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5875 else if (round == 1)
5877 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5886 else if (round == 2)
5887 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5900 std::vector<CellData<dim>> cells;
5901 cells.reserve(tria_copy.
n_cells());
5902 for (
const auto &cell : tria_copy.cell_iterators())
5906 data.
vertices[v] = cell->vertex_index(v);
5909 cells.push_back(data);
5917 if (round == dim - 1)
5926 for (
const auto &cell :
tria.cell_iterators())
5927 if (cell->
center().norm_square() > 0.4 * radius)
5928 cell->set_manifold_id(1);
5949 const double inner_radius,
5950 const double outer_radius)
5953 std::vector<CellData<3>> cells;
5955 const double irad = inner_radius /
std::sqrt(3.0);
5956 const double orad = outer_radius /
std::sqrt(3.0);
5959 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
5969 for (
unsigned int i = 0; i < 8; ++i)
5970 vertices.push_back(p + hexahedron[i] * irad);
5971 for (
unsigned int i = 0; i < 8; ++i)
5972 vertices.push_back(p + hexahedron[i] * orad);
5974 const unsigned int n_cells = 6;
5975 const int cell_vertices[
n_cells][8] = {
5976 {8, 9, 10, 11, 0, 1, 2, 3},
5977 {9, 11, 1, 3, 13, 15, 5, 7},
5978 {12, 13, 4, 5, 14, 15, 6, 7},
5979 {8, 0, 10, 2, 12, 4, 14, 6},
5980 {8, 9, 0, 1, 12, 13, 4, 5},
5981 {10, 2, 11, 3, 14, 6, 15, 7}};
5985 for (
unsigned int i = 0; i <
n_cells; ++i)
5988 cells[i].
vertices[j] = cell_vertices[i][j];
5989 cells[i].material_id = 0;
6000 const double inner_radius,
6001 const double outer_radius)
6004 std::vector<CellData<3>> cells;
6006 const double irad = inner_radius /
std::sqrt(3.0);
6007 const double orad = outer_radius /
std::sqrt(3.0);
6013 static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0},
6021 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
6030 for (
unsigned int i = 0; i < 8; ++i)
6031 vertices.push_back(p + hexahedron[i] * irad);
6032 for (
unsigned int i = 0; i < 6; ++i)
6033 vertices.push_back(p + octahedron[i] * inner_radius);
6034 for (
unsigned int i = 0; i < 8; ++i)
6035 vertices.push_back(p + hexahedron[i] * orad);
6036 for (
unsigned int i = 0; i < 6; ++i)
6037 vertices.push_back(p + octahedron[i] * outer_radius);
6039 const unsigned int n_cells = 12;
6040 const unsigned int rhombi[
n_cells][4] = {{10, 4, 0, 8},
6055 for (
unsigned int i = 0; i <
n_cells; ++i)
6057 for (
unsigned int j = 0; j < 4; ++j)
6059 cells[i].vertices[j] = rhombi[i][j];
6060 cells[i].vertices[j + 4] = rhombi[i][j] + 14;
6062 cells[i].material_id = 0;
6072 const unsigned int n,
6073 const unsigned int n_refinement_steps,
6075 const double inner_radius,
6076 const double outer_radius)
6089 const unsigned int outer_radius_factor = 1 << n_refinement_steps;
6094 outer_radius_factor * outer_radius -
6095 (outer_radius_factor - 1) * inner_radius);
6100 outer_radius_factor * outer_radius -
6101 (outer_radius_factor - 1) * inner_radius);
6104 for (
unsigned int r = 0; r < n_refinement_steps; ++r)
6107 std::set<Triangulation<3>::active_cell_iterator> cells_to_remove;
6111 for (
const auto &cell : tmp.active_cell_iterators())
6113 unsigned int n_vertices_inside = 0;
6115 if ((cell->vertex(v) -
p).norm_square() <
6116 inner_radius * inner_radius * (1 + 1
e-12))
6117 ++n_vertices_inside;
6118 if (n_vertices_inside < 4)
6119 cells_to_remove.insert(cell);
6123 if (r == n_refinement_steps - 1)
6133 tmp = std::move(copy);
6151 const double inner_radius,
6152 const double outer_radius,
6153 const unsigned int n_cells,
6156 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6159 unsigned int n_refinement_steps = 0;
6160 unsigned int n_cells_coarsened =
n_cells;
6161 if (n_cells != 96 && n_cells > 12)
6162 while (n_cells_coarsened > 12 && n_cells_coarsened % 4 == 0)
6164 ++n_refinement_steps;
6165 n_cells_coarsened /= 4;
6167 Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 ||
6168 (n_refinement_steps > 0 &&
6169 (n_cells_coarsened == 6 || n_cells_coarsened == 12)),
6170 ExcMessage(
"Invalid number of coarse mesh cells"));
6172 const unsigned int n = n_refinement_steps > 0 ?
6173 4 * n_cells_coarsened :
6179 internal::hyper_shell_6(tria, p, inner_radius, outer_radius);
6182 internal::hyper_shell_12(tria, p, inner_radius, outer_radius);
6186 internal::hyper_shell_24_48(
6187 tria, n, n_refinement_steps, p, inner_radius, outer_radius);
6196 internal::hyper_shell_12(tmp, p, inner_radius, outer_radius);
6213 colorize_hyper_shell(tria, p, inner_radius, outer_radius);
6223 const double inner_radius,
6224 const double outer_radius,
6225 const unsigned int ,
6228 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6232 const double d = outer_radius /
std::sqrt(2.0);
6233 const double a = inner_radius /
std::sqrt(2.0);
6235 const double b = a / 2.0;
6236 const double c =
d / 2.0;
6238 const double hb = inner_radius *
std::sqrt(3.0) / 2.0;
6239 const double hc = outer_radius *
std::sqrt(3.0) / 2.0;
6261 int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
6262 {0, 2, 8, 10, 6, 4, 14, 12},
6263 {1, 7, 9, 15, 3, 5, 11, 13},
6264 {6, 4, 14, 12, 7, 5, 15, 13},
6265 {8, 10, 9, 11, 14, 12, 15, 13}};
6269 for (
unsigned int i = 0; i < 5; ++i)
6271 for (
unsigned int j = 0; j < 8; ++j)
6272 cells[i].
vertices[j] = cell_vertices[i][j];
6273 cells[i].material_id = 0;
6287 for (; cell !=
tria.
end(); ++cell)
6288 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6289 if (cell->at_boundary(i))
6290 cell->face(i)->set_all_boundary_ids(2);
6297 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6298 if (cell->at_boundary(i))
6302 const Point<3> face_center(face->center());
6304 1.e-6 * face_center.
norm())
6308 face->set_all_boundary_ids(0);
6310 face->set_all_boundary_ids(1);
6324 const double inner_radius,
6325 const double outer_radius,
6326 const unsigned int n,
6329 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6331 if (n == 0 || n == 3)
6333 const double a = inner_radius *
std::sqrt(2.0) / 2e0;
6334 const double b = outer_radius *
std::sqrt(2.0) / 2e0;
6335 const double c = a *
std::sqrt(3.0) / 2e0;
6337 const double e = outer_radius / 2e0;
6338 const double h = inner_radius / 2e0;
6357 const int cell_vertices[3][8] = {
6358 {0, 1, 3, 2, 4, 5, 7, 6},
6359 {1, 8, 2, 9, 5, 10, 6, 11},
6360 {4, 5, 7, 6, 12, 10, 13, 11},
6362 std::vector<CellData<3>> cells(3);
6364 for (
unsigned int i = 0; i < 3; ++i)
6366 for (
unsigned int j = 0; j < 8; ++j)
6367 cells[i].
vertices[j] = cell_vertices[i][j];
6368 cells[i].material_id = 0;
6381 colorize_quarter_hyper_shell(tria,
center, inner_radius, outer_radius);
6392 const double length,
6393 const double inner_radius,
6394 const double outer_radius,
6395 const unsigned int n_radial_cells,
6396 const unsigned int n_axial_cells,
6399 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6413 const unsigned int N_r =
6414 (n_radial_cells == 0 ?
static_cast<unsigned int>(std::ceil(
6415 (2 * pi * (outer_radius + inner_radius) / 2) /
6416 (outer_radius - inner_radius))) :
6418 const unsigned int N_z =
6419 (n_axial_cells == 0 ?
6420 static_cast<unsigned int>(std::ceil(
6421 length / (2 * pi * (outer_radius + inner_radius) / 2 / N_r))) :
6430 std::vector<Point<2>> vertices_2d(2 * N_r);
6431 for (
unsigned int i = 0; i < N_r; ++i)
6436 vertices_2d[i + N_r] = vertices_2d[i] * (inner_radius / outer_radius);
6439 std::vector<Point<3>> vertices_3d;
6440 vertices_3d.reserve(2 * N_r * (N_z + 1));
6441 for (
unsigned int j = 0; j <= N_z; ++j)
6442 for (
unsigned int i = 0; i < 2 * N_r; ++i)
6444 const Point<3> v(vertices_2d[i][0],
6447 vertices_3d.push_back(v);
6450 std::vector<CellData<3>> cells(N_r * N_z,
CellData<3>());
6452 for (
unsigned int j = 0; j < N_z; ++j)
6453 for (
unsigned int i = 0; i < N_r; ++i)
6455 cells[i + j * N_r].vertices[0] = i + (j + 1) * 2 * N_r;
6456 cells[i + j * N_r].vertices[1] = (i + 1) % N_r + (j + 1) * 2 * N_r;
6457 cells[i + j * N_r].vertices[2] = i + j * 2 * N_r;
6458 cells[i + j * N_r].vertices[3] = (i + 1) % N_r + j * 2 * N_r;
6460 cells[i + j * N_r].vertices[4] = N_r + i + (j + 1) * 2 * N_r;
6461 cells[i + j * N_r].vertices[5] =
6462 N_r + ((i + 1) % N_r) + (j + 1) * 2 * N_r;
6463 cells[i + j * N_r].vertices[6] = N_r + i + j * 2 * N_r;
6464 cells[i + j * N_r].vertices[7] = N_r + ((i + 1) % N_r) + j * 2 * N_r;
6466 cells[i + j * N_r].material_id = 0;
6485 double eps_z = 1
e-6 * length;
6490 double face_inner_radius = std::numeric_limits<double>::max();
6491 double face_outer_radius = 0.;
6496 for (
const auto &cell :
tria.active_cell_iterators())
6499 if (!cell->face(f)->at_boundary())
6502 const auto face_center = cell->face(f)->center();
6503 const double z = face_center[2];
6505 if ((std::fabs(z) > eps_z) &&
6506 (std::fabs(z - length) > eps_z))
6508 const double radius =
std::sqrt(face_center[0] * face_center[0] +
6509 face_center[1] * face_center[1]);
6510 face_inner_radius =
std::min(face_inner_radius, radius);
6511 face_outer_radius =
std::max(face_outer_radius, radius);
6515 double mid_radial_distance = 0.5 * (face_outer_radius - face_inner_radius);
6517 for (
const auto &cell :
tria.active_cell_iterators())
6520 if (cell->face(f)->at_boundary())
6522 const auto face_center = cell->face(f)->center();
6524 const double radius =
std::sqrt(face_center[0] * face_center[0] +
6525 face_center[1] * face_center[1]);
6527 const double z = face_center[2];
6529 if (std::fabs(z) < eps_z)
6531 cell->face(f)->set_boundary_id(2);
6533 else if (std::fabs(z - length) <
6536 cell->face(f)->set_boundary_id(3);
6538 else if (std::fabs(radius - face_inner_radius) >
6539 mid_radial_distance)
6541 cell->face(f)->set_boundary_id(1);
6543 else if (std::fabs(radius - face_inner_radius) <
6544 mid_radial_distance)
6546 cell->face(f)->set_boundary_id(0);
6555 template <
int dim,
int spacedim>
6560 const double duplicated_vertex_tolerance,
6561 const bool copy_manifold_ids,
6562 const bool copy_boundary_ids)
6564 std::vector<Point<spacedim>>
vertices;
6565 std::vector<CellData<dim>> cells;
6568 unsigned int n_accumulated_vertices = 0;
6572 ExcMessage(
"The input triangulations must be non-empty "
6573 "and must not be refined."));
6575 auto [tria_vertices, tria_cells, tria_subcell_data] =
6586 tria_vertices.begin(),
6587 tria_vertices.end());
6590 for (
unsigned int &vertex_n : cell_data.
vertices)
6591 vertex_n += n_accumulated_vertices;
6592 cells.push_back(cell_data);
6596 if (copy_manifold_ids)
6598 for (
CellData<1> &line_data : tria_subcell_data.boundary_lines)
6602 for (
unsigned int &vertex_n : line_data.
vertices)
6603 vertex_n += n_accumulated_vertices;
6604 line_data.boundary_id =
6609 for (
CellData<2> &quad_data : tria_subcell_data.boundary_quads)
6613 for (
unsigned int &vertex_n : quad_data.
vertices)
6614 vertex_n += n_accumulated_vertices;
6615 quad_data.boundary_id =
6625 std::vector<unsigned int> considered_vertices;
6629 considered_vertices,
6630 duplicated_vertex_tolerance);
6634 if (std::all_of(cells.begin(), cells.end(), [](
const auto &cell) {
6635 return cell.vertices.size() ==
6636 ReferenceCells::get_hypercube<dim>().n_vertices();
6642 if (copy_boundary_ids)
6644 auto result_cell = result.
begin();
6645 for (
const auto &tria : triangulations)
6647 for (
const auto &cell :
tria->cell_iterators())
6649 for (
const auto &f : cell->face_indices())
6650 if (result_cell->face(f)->at_boundary())
6651 result_cell->face(f)->set_boundary_id(
6658 Assert(duplicated_vertex_tolerance > 0.0 ||
6659 n_accumulated_vertices == result.
n_vertices(),
6665 template <
int dim,
int spacedim>
6670 const double duplicated_vertex_tolerance,
6671 const bool copy_manifold_ids,
6672 const bool copy_boundary_ids)
6675 if (triangulation_1.
n_cells() == 0)
6677 if (&result != &triangulation_2)
6680 else if (triangulation_2.
n_cells() == 0)
6682 if (&result != &triangulation_1)
6688 duplicated_vertex_tolerance,
6718 template <
int structdim>
6722 static_assert(structdim == 1 || structdim == 2,
6723 "This function is only implemented for lines and "
6731 std::sort(std::begin(cell_data.vertices),
6732 std::end(cell_data.vertices));
6733 else if (structdim == 2)
6736 std::array<unsigned int, 4> renumbering{};
6737 std::copy(std::begin(cell_data.vertices),
6738 std::end(cell_data.vertices),
6739 renumbering.begin());
6750 std::swap(renumbering[2], renumbering[3]);
6751 std::rotate(renumbering.begin(),
6752 std::min_element(renumbering.begin(),
6756 std::swap(renumbering[2], renumbering[3]);
6764 if (renumbering[1] > renumbering[2])
6765 std::swap(renumbering[1], renumbering[2]);
6766 std::copy(renumbering.begin(),
6768 std::begin(cell_data.vertices));
6775 return std::lexicographical_compare(std::begin(a.
vertices),
6777 std::begin(
b.vertices),
6778 std::end(
b.vertices));
6780 std::sort(subcell_data.begin(), subcell_data.end(), compare);
6785 auto left = subcell_data.begin();
6786 while (left != subcell_data.end())
6789 std::upper_bound(left, subcell_data.end(), *left, compare);
6792 if (left + 1 != right)
6793 for (
auto it = left;
it != right; ++
it)
6796 Assert(
it->manifold_id == left->manifold_id,
6798 "In the process of grid generation a single "
6799 "line or quadrilateral has been assigned two "
6800 "different manifold ids. This can happen when "
6801 "a Triangulation is copied, e.g., via "
6802 "GridGenerator::replicate_triangulation() and "
6803 "not all external boundary faces have the same "
6804 "manifold id. Double check that all faces "
6805 "which you expect to be merged together have "
6806 "the same manifold id."));
6811 subcell_data.erase(std::unique(subcell_data.begin(), subcell_data.end()),
6812 subcell_data.end());
6818 template <
int dim,
int spacedim>
6821 const std::vector<unsigned int> &extents,
6826 for (
const auto &extent : extents)
6828 ExcMessage(
"The Triangulation must be copied at least one time in "
6829 "each coordinate dimension."));
6832 const auto &
min = bbox.get_boundary_points().first;
6833 const auto &
max = bbox.get_boundary_points().second;
6835 std::array<Tensor<1, spacedim>, dim> offsets;
6836 for (
unsigned int d = 0;
d < dim; ++
d)
6837 offsets[d][d] = max[d] - min[d];
6841 for (
unsigned int d = 0;
d < dim; ++
d)
6843 auto [input_vertices, input_cell_data, input_subcell_data] =
6846 std::vector<Point<spacedim>> output_vertices = input_vertices;
6847 std::vector<CellData<dim>> output_cell_data = input_cell_data;
6848 SubCellData output_subcell_data = input_subcell_data;
6850 for (
unsigned int k = 1; k < extents[
d]; ++k)
6852 const std::size_t vertex_offset = k * input_vertices.size();
6855 output_vertices.push_back(
point + double(k) * offsets[
d]);
6859 output_cell_data.push_back(cell_data);
6860 for (
unsigned int &vertex : output_cell_data.back().
vertices)
6861 vertex += vertex_offset;
6865 input_subcell_data.boundary_lines)
6868 for (
unsigned int &vertex :
6869 output_subcell_data.boundary_lines.back().
vertices)
6870 vertex += vertex_offset;
6873 input_subcell_data.boundary_quads)
6876 for (
unsigned int &vertex :
6877 output_subcell_data.boundary_quads.back().
vertices)
6878 vertex += vertex_offset;
6883 std::vector<unsigned int> boundary_vertices;
6887 output_subcell_data,
6902 output_subcell_data.boundary_lines)
6908 tria_to_replicate.
clear();
6911 output_subcell_data);
6919 template <
int dim,
int spacedim>
6927 ExcMessage(
"The two input triangulations are not derived from "
6928 "the same coarse mesh as required."));
6931 &triangulation_1) ==
nullptr) &&
6934 &triangulation_2) ==
nullptr),
6935 ExcMessage(
"The source triangulations for this function must both "
6936 "be available entirely locally, and not be distributed "
6937 "triangulations."));
6954 for (
unsigned int iteration = 0; iteration < triangulation_2.
n_levels();
6960 bool any_cell_flagged =
false;
6961 for (
const auto &result_cell : result.active_cell_iterators())
6962 if (intergrid_map[result_cell]->has_children())
6964 any_cell_flagged =
true;
6965 result_cell->set_refine_flag();
6968 if (any_cell_flagged ==
false)
6977 template <
int dim,
int spacedim>
6991 std::vector<CellData<dim>> cells;
6992 for (
const auto &cell : input_triangulation.active_cell_iterators())
6993 if (cells_to_remove.find(cell) == cells_to_remove.
end())
6995 Assert(
static_cast<unsigned int>(cell->level()) ==
6996 input_triangulation.
n_levels() - 1,
6998 "Your input triangulation appears to have "
6999 "adaptively refined cells. This is not allowed. You can "
7000 "only call this function on a triangulation in which "
7001 "all cells are on the same refinement level."));
7005 this_cell.
vertices[v] = cell->vertex_index(v);
7007 cells.push_back(this_cell);
7013 std::vector<unsigned int> considered_vertices;
7017 considered_vertices);
7029 const unsigned int n_slices,
7030 const double height,
7032 const bool copy_manifold_ids,
7033 const std::vector<types::manifold_id> &manifold_priorities)
7037 "The input triangulation must be a coarse mesh, i.e., it must "
7038 "not have been refined."));
7040 ExcMessage(
"The output triangulation object needs to be empty."));
7042 ExcMessage(
"The given height for extrusion must be positive."));
7045 "The number of slices for extrusion must be at least 2."));
7047 const double delta_h = height / (n_slices - 1);
7048 std::vector<double> slices_z_values;
7049 for (
unsigned int i = 0; i < n_slices; ++i)
7050 slices_z_values.push_back(i * delta_h);
7052 input, slices_z_values, result, copy_manifold_ids, manifold_priorities);
7060 const unsigned int n_slices,
7061 const double height,
7063 const bool copy_manifold_ids,
7064 const std::vector<types::manifold_id> &manifold_priorities)
7070 (void)copy_manifold_ids;
7071 (void)manifold_priorities;
7075 "GridTools::extrude_triangulation() is only available "
7076 "for Triangulation<3, 3> as output triangulation."));
7084 const std::vector<double> &slice_coordinates,
7086 const bool copy_manifold_ids,
7087 const std::vector<types::manifold_id> &manifold_priorities)
7091 "The input triangulation must be a coarse mesh, i.e., it must "
7092 "not have been refined."));
7094 ExcMessage(
"The output triangulation object needs to be empty."));
7095 Assert(slice_coordinates.size() >= 2,
7097 "The number of slices for extrusion must be at least 2."));
7098 Assert(std::is_sorted(slice_coordinates.begin(), slice_coordinates.end()),
7099 ExcMessage(
"Slice z-coordinates should be in ascending order"));
7102 "This function is only implemented for quadrilateral meshes."));
7104 const auto priorities = [&]() -> std::vector<types::manifold_id> {
7108 if (0 < manifold_priorities.size())
7112 std::vector<types::manifold_id> sorted_manifold_priorities =
7113 manifold_priorities;
7114 std::sort(sorted_manifold_priorities.begin(),
7115 sorted_manifold_priorities.end());
7116 Assert(std::unique(sorted_manifold_priorities.begin(),
7117 sorted_manifold_priorities.end()) ==
7118 sorted_manifold_priorities.end(),
7120 "The given vector of manifold ids may not contain any "
7121 "duplicated entries."));
7122 std::vector<types::manifold_id> sorted_manifold_ids =
7124 std::sort(sorted_manifold_ids.begin(), sorted_manifold_ids.end());
7125 if (sorted_manifold_priorities != sorted_manifold_ids)
7127 std::ostringstream message;
7128 message <<
"The given triangulation has manifold ids {";
7132 message << sorted_manifold_ids.back() <<
"}, but \n"
7133 <<
" the given vector of manifold ids is {";
7138 << manifold_priorities.back() <<
"}.\n"
7139 <<
" These vectors should contain the same elements.\n";
7140 const std::string m = message.str();
7144 return manifold_priorities;
7148 std::vector<types::manifold_id> default_priorities =
7150 const auto first_tfi_it = std::partition(
7151 default_priorities.begin(),
7152 default_priorities.end(),
7154 return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
7155 &input.get_manifold(id)) == nullptr;
7157 std::sort(default_priorities.begin(), first_tfi_it);
7158 std::sort(first_tfi_it, default_priorities.end());
7160 return default_priorities;
7163 const std::size_t n_slices = slice_coordinates.size();
7164 std::vector<Point<3>> points(n_slices * input.
n_vertices());
7165 std::vector<CellData<3>> cells;
7170 for (std::size_t slice_n = 0; slice_n < n_slices; ++slice_n)
7172 for (std::size_t vertex_n = 0; vertex_n < input.
n_vertices();
7176 points[slice_n * input.
n_vertices() + vertex_n] =
7177 Point<3>(vertex[0], vertex[1], slice_coordinates[slice_n]);
7183 for (
const auto &cell : input.active_cell_iterators())
7185 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
7188 for (
const unsigned int vertex_n :
7192 cell->vertex_index(vertex_n) + slice_n * input.
n_vertices();
7195 cell->vertex_index(vertex_n) +
7200 if (copy_manifold_ids)
7202 cells.push_back(this_cell);
7213 for (
const auto &face : input.active_face_iterators())
7217 if (face->at_boundary())
7219 if (copy_manifold_ids)
7221 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
7224 face->vertex_index(0) + slice_n * input.
n_vertices();
7226 face->vertex_index(1) + slice_n * input.
n_vertices();
7228 face->vertex_index(0) + (slice_n + 1) * input.
n_vertices();
7230 face->vertex_index(1) + (slice_n + 1) * input.
n_vertices();
7231 quads.push_back(quad);
7237 if (copy_manifold_ids)
7238 for (
const auto &cell : input.active_cell_iterators())
7243 for (std::size_t slice_n = 1; slice_n < n_slices - 1; ++slice_n)
7246 cell->vertex_index(0) + slice_n * input.
n_vertices();
7248 cell->vertex_index(1) + slice_n * input.
n_vertices();
7250 cell->vertex_index(2) + slice_n * input.
n_vertices();
7252 cell->vertex_index(3) + slice_n * input.
n_vertices();
7253 quads.push_back(quad);
7264 "The input triangulation to this function is using boundary "
7265 "indicators in a range that do not allow using "
7266 "max_boundary_id+1 and max_boundary_id+2 as boundary "
7267 "indicators for the bottom and top faces of the "
7268 "extruded triangulation."));
7271 for (
const auto &cell : input.active_cell_iterators())
7275 quad.
vertices[0] = cell->vertex_index(0);
7276 quad.
vertices[1] = cell->vertex_index(1);
7277 quad.
vertices[2] = cell->vertex_index(2);
7278 quad.
vertices[3] = cell->vertex_index(3);
7279 if (copy_manifold_ids)
7281 quads.push_back(quad);
7284 for (
unsigned int &vertex : quad.
vertices)
7285 vertex += (n_slices - 1) * input.n_vertices();
7286 if (copy_manifold_ids)
7288 quads.push_back(quad);
7301 for (
auto manifold_id_it = priorities.rbegin();
7302 manifold_id_it != priorities.rend();
7304 for (
const auto &face : result.active_face_iterators())
7306 for (unsigned
int line_n = 0;
7307 line_n < GeometryInfo<3>::lines_per_face;
7309 face->line(line_n)->set_manifold_id(*manifold_id_it);
7317 const std::vector<double> &slice_coordinates,
7319 const bool copy_manifold_ids,
7320 const std::vector<types::manifold_id> &manifold_priorities)
7323 (void)slice_coordinates;
7325 (void)copy_manifold_ids;
7326 (void)manifold_priorities;
7330 "GridTools::extrude_triangulation() is only available "
7331 "for Triangulation<3, 3> as output triangulation."));
7353 const double inner_radius,
7354 const double outer_radius,
7361 Assert(inner_radius < outer_radius,
7362 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7371 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
7372 for (
const auto &cell :
triangulation.active_cell_iterators())
7375 if (cell->face(f)->at_boundary())
7378 (
std::
fabs(cell->face(f)->vertex(v).
norm() - outer_radius) <
7379 1
e-12 * outer_radius)
7384 std::
fabs(cell->face(f)->vertex(v)[1])) <
7385 1
e-12 * outer_radius))
7386 cell->face(f)->vertex(v) *=
std::
sqrt(2.);
7388 const double eps = 1
e-3 * outer_radius;
7389 for (
const auto &cell :
triangulation.active_cell_iterators())
7391 for (
const unsigned int f : cell->face_indices())
7392 if (cell->face(f)->at_boundary())
7394 const double dx = cell->face(f)->center()[0] -
center[0];
7395 const double dy = cell->face(f)->center()[1] -
center[1];
7398 if (
std::abs(dx + outer_radius) < eps)
7399 cell->face(f)->set_boundary_id(0);
7400 else if (
std::abs(dx - outer_radius) < eps)
7401 cell->face(f)->set_boundary_id(1);
7402 else if (
std::abs(dy + outer_radius) < eps)
7403 cell->face(f)->set_boundary_id(2);
7404 else if (
std::abs(dy - outer_radius) < eps)
7405 cell->face(f)->set_boundary_id(3);
7408 cell->face(f)->set_boundary_id(4);
7409 cell->face(f)->set_manifold_id(0);
7414 const double d = (cell->face(f)->center() -
center).
norm();
7415 if (d - inner_radius < 0)
7417 cell->face(f)->set_boundary_id(1);
7418 cell->face(f)->set_manifold_id(0);
7421 cell->face(f)->set_boundary_id(0);
7434 const double inner_radius,
7435 const double outer_radius,
7436 const unsigned int n_shells,
7437 const double skewness,
7438 const unsigned int n_cells,
7444 Assert(inner_radius < outer_radius,
7445 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7449 std::vector<double> radii;
7450 radii.push_back(inner_radius);
7451 for (
unsigned int shell_n = 1; shell_n < n_shells; ++shell_n)
7452 if (skewness == 0.0)
7454 radii.push_back(inner_radius +
7455 (outer_radius - inner_radius) *
7456 (1.0 - (1.0 -
double(shell_n) / n_shells)));
7460 (outer_radius - inner_radius) *
7461 (1.0 - std::tanh(skewness * (1.0 -
double(shell_n) / n_shells)) /
7462 std::tanh(skewness)));
7463 radii.push_back(outer_radius);
7465 double grid_vertex_tolerance = 0.0;
7466 for (
unsigned int shell_n = 0; shell_n < radii.size() - 1; ++shell_n)
7473 n_cells == 0 ? (dim == 2 ? 8 : 12) :
7478 if (grid_vertex_tolerance == 0.0)
7479 grid_vertex_tolerance =
7480 0.5 * internal::minimal_vertex_distance(current_shell);
7487 grid_vertex_tolerance);
7499 constexpr double radial_vertex_tolerance =
7500 100.0 * std::numeric_limits<double>::epsilon();
7501 auto assert_vertex_distance_within_tolerance =
7502 [
center, radial_vertex_tolerance](
7504 const double radius) {
7506 (void)radial_vertex_tolerance;
7509 for (
unsigned int vertex_n = 0;
7510 vertex_n < GeometryInfo<dim>::vertices_per_face;
7514 (
center.
norm() + radius) * radial_vertex_tolerance,
7519 for (
const auto &cell :
triangulation.active_cell_iterators())
7522 auto face = cell->face(face_n);
7523 if (face->at_boundary())
7525 if (((face->vertex(0) -
center).norm() - inner_radius) <
7526 (
center.
norm() + inner_radius) * radial_vertex_tolerance)
7529 assert_vertex_distance_within_tolerance(face, inner_radius);
7530 face->set_all_boundary_ids(0);
7535 assert_vertex_distance_within_tolerance(face, outer_radius);
7536 face->set_all_boundary_ids(1);
7547 const double inner_radius,
7548 const double outer_radius,
7550 const unsigned int Nz,
7555 Assert(inner_radius < outer_radius,
7556 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7570 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
7571 for (
const auto &cell :
triangulation.active_cell_iterators())
7573 for (
const auto f : cell->face_indices())
7574 if (cell->face(f)->at_boundary())
7578 const unsigned int vv = cell->face(f)->vertex_index(v);
7579 if (treated_vertices[vv] ==
false)
7581 treated_vertices[vv] =
true;
7591 cell->face(f)->vertex(v);
7592 if ((std::fabs(std::fabs(vertex_location[0]) -
7593 std::fabs(vertex_location[1])) <
7594 1e-12 * outer_radius) &&
7595 (std::fabs(vertex_location[0] * vertex_location[0] +
7596 vertex_location[1] * vertex_location[1] -
7597 outer_radius * outer_radius) <
7598 1e-12 * outer_radius))
7599 cell->face(f)->vertex(v) =
7602 vertex_location[2]);
7607 double eps = 1
e-3 * outer_radius;
7608 for (
const auto &cell :
triangulation.active_cell_iterators())
7610 for (
const unsigned int f : cell->face_indices())
7611 if (cell->face(f)->at_boundary())
7613 const double dx = cell->face(f)->center()[0];
7614 const double dy = cell->face(f)->center()[1];
7615 const double dz = cell->face(f)->center()[2];
7619 if (
std::abs(dx + outer_radius) < eps)
7620 cell->face(f)->set_boundary_id(0);
7622 else if (
std::abs(dx - outer_radius) < eps)
7623 cell->face(f)->set_boundary_id(1);
7625 else if (
std::abs(dy + outer_radius) < eps)
7626 cell->face(f)->set_boundary_id(2);
7628 else if (
std::abs(dy - outer_radius) < eps)
7629 cell->face(f)->set_boundary_id(3);
7632 cell->face(f)->set_boundary_id(4);
7635 cell->face(f)->set_boundary_id(5);
7639 cell->face(f)->set_all_boundary_ids(6);
7640 cell->face(f)->set_all_manifold_ids(0);
7647 const double d = c.
norm();
7648 if (d - inner_radius < 0)
7650 cell->face(f)->set_all_boundary_ids(1);
7651 cell->face(f)->set_all_manifold_ids(0);
7654 cell->face(f)->set_boundary_id(0);
7663 template <
int dim,
int spacedim1,
int spacedim2>
7670 &in_tria) ==
nullptr),
7672 "This function cannot be used on "
7673 "parallel::distributed::Triangulation objects as inputs."));
7675 ExcMessage(
"This function does not work for meshes that have "
7679 const unsigned int spacedim =
std::min(spacedim1, spacedim2);
7680 const std::vector<Point<spacedim1>> &in_vertices = in_tria.
get_vertices();
7684 std::vector<Point<spacedim2>> v(in_vertices.size());
7685 for (
unsigned int i = 0; i < in_vertices.size(); ++i)
7686 for (
unsigned int d = 0;
d < spacedim; ++
d)
7687 v[i][d] = in_vertices[i][d];
7690 for (
const auto &cell : in_tria.active_cell_iterators())
7692 const unsigned int id = cell->active_cell_index();
7694 cells[id].vertices.resize(cell->n_vertices());
7696 cells[id].
vertices[i] = cell->vertex_index(i);
7697 cells[id].material_id = cell->material_id();
7698 cells[id].manifold_id = cell->manifold_id();
7712 std::vector<bool> user_flags_line;
7715 .clear_user_flags_line();
7719 for (
const auto &face : in_tria.active_face_iterators())
7721 if (face->at_boundary())
7725 boundary_line.
vertices.resize(face->n_vertices());
7727 boundary_line.
vertices[i] = face->vertex_index(i);
7732 std::move(boundary_line));
7745 if ((face->user_flag_set() ==
false) &&
7750 boundary_line.
vertices.resize(face->n_vertices());
7752 boundary_line.
vertices[i] = face->vertex_index(i);
7758 std::move(boundary_line));
7760 face->set_user_flag();
7766 .load_user_flags_line(user_flags_line);
7773 std::vector<bool> user_flags_line;
7776 .clear_user_flags_line();
7778 std::vector<bool> user_flags_quad;
7781 .clear_user_flags_quad();
7785 for (
const auto &face : in_tria.active_face_iterators())
7787 if (face->at_boundary())
7791 boundary_face.
vertices.resize(face->n_vertices());
7793 boundary_face.
vertices[i] = face->vertex_index(i);
7798 std::move(boundary_face));
7805 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7806 if (face->line(e)->user_flag_set() ==
false)
7809 spacedim1>::line_iterator
7810 edge = face->line(e);
7813 boundary_edge.
vertices.resize(edge->n_vertices());
7815 boundary_edge.
vertices[i] = edge->vertex_index(i);
7820 std::move(boundary_edge));
7822 edge->set_user_flag();
7840 if (face->user_flag_set() ==
false)
7846 boundary_face.
vertices.resize(face->n_vertices());
7848 boundary_face.
vertices[i] = face->vertex_index(i);
7854 std::move(boundary_face));
7856 face->set_user_flag();
7865 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7866 if (face->line(e)->at_boundary() ==
false)
7867 if (face->line(e)->user_flag_set() ==
false)
7870 line_iterator edge = face->line(e);
7876 edge->vertex_index(i);
7882 std::move(boundary_edge));
7884 edge->set_user_flag();
7891 .load_user_flags_line(user_flags_line);
7893 .load_user_flags_quad(user_flags_quad);
7903 for (
const auto i : out_tria.get_manifold_ids())
7905 out_tria.set_manifold(i,
FlatManifold<dim, spacedim2>());
7910 template <
int dim,
int spacedim>
7950 {{{0, 1, 12, 10}}, {{2, 3, 11, 12}}, {{7, 6, 11, 13}},
7951 {{5, 4, 13, 10}}, {{0, 2, 8, 12}}, {{4, 6, 13, 8}},
7952 {{5, 13, 7, 9}}, {{1, 9, 3, 12}}, {{0, 8, 4, 10}},
7953 {{1, 5, 9, 10}}, {{3, 7, 11, 9}}, {{2, 6, 8, 11}},
7954 {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
7955 {{12, 13, 8, 10}}, {{13, 8, 10, 4}}, {{13, 10, 9, 5}},
7956 {{13, 9, 11, 7}}, {{13, 11, 8, 6}}, {{10, 12, 9, 1}},
7957 {{9, 12, 11, 3}}, {{11, 12, 8, 2}}, {{8, 12, 10, 0}}}};
7965 vertex_ids_for_boundary_faces_2d = {{{{{{0, 4}}, {{4, 2}}}},
7966 {{{{1, 5}}, {{5, 3}}}},
7967 {{{{0, 6}}, {{6, 1}}}},
7968 {{{{2, 7}}, {{7, 3}}}}}};
7976 vertex_ids_for_boundary_faces_3d = {
7977 {{{{{0, 4, 8}}, {{4, 8, 6}}, {{8, 6, 2}}, {{0, 2, 8}}}},
7978 {{{{1, 3, 9}}, {{3, 9, 7}}, {{9, 7, 5}}, {{1, 9, 5}}}},
7979 {{{{0, 1, 10}}, {{1, 10, 5}}, {{10, 5, 4}}, {{0, 10, 4}}}},
7980 {{{{2, 3, 11}}, {{3, 11, 7}}, {{11, 7, 6}}, {{2, 11, 6}}}},
7981 {{{{0, 1, 12}}, {{1, 12, 3}}, {{12, 3, 2}}, {{0, 12, 2}}}},
7982 {{{{4, 5, 13}}, {{5, 13, 7}}, {{13, 7, 6}}, {{4, 13, 6}}}}}};
8003 {{{0, 12, 10}}, {{12, 1, 10}}, {{12, 1, 9}}, {{12, 3, 9}},
8004 {{12, 2, 11}}, {{12, 3, 11}}, {{12, 0, 8}}, {{12, 2, 8}},
8005 {{9, 13, 5}}, {{13, 7, 9}}, {{11, 7, 13}}, {{11, 6, 13}},
8006 {{4, 8, 13}}, {{6, 8, 13}}, {{4, 13, 10}}, {{13, 5, 10}},
8007 {{10, 9, 5}}, {{10, 9, 1}}, {{11, 9, 7}}, {{11, 9, 3}},
8008 {{8, 11, 2}}, {{8, 11, 6}}, {{8, 10, 0}}, {{8, 10, 4}},
8009 {{12, 3, 9}}, {{12, 9, 11}}, {{12, 3, 11}}, {{3, 9, 11}},
8010 {{2, 12, 8}}, {{2, 12, 11}}, {{2, 11, 8}}, {{8, 12, 11}},
8011 {{0, 12, 10}}, {{0, 12, 8}}, {{0, 8, 10}}, {{8, 10, 12}},
8012 {{12, 1, 10}}, {{12, 1, 9}}, {{1, 10, 9}}, {{10, 9, 12}},
8013 {{10, 8, 4}}, {{10, 8, 13}}, {{4, 13, 8}}, {{4, 13, 10}},
8014 {{10, 9, 13}}, {{10, 9, 5}}, {{13, 5, 10}}, {{13, 5, 9}},
8015 {{13, 7, 9}}, {{13, 7, 11}}, {{9, 11, 13}}, {{9, 11, 7}},
8016 {{8, 11, 13}}, {{8, 11, 6}}, {{6, 13, 8}}, {{6, 13, 11}},
8017 {{12, 13, 10}}, {{12, 13, 8}}, {{8, 10, 13}}, {{8, 10, 12}},
8018 {{12, 13, 10}}, {{12, 13, 9}}, {{10, 9, 13}}, {{10, 9, 12}},
8019 {{12, 13, 9}}, {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
8020 {{12, 13, 11}}, {{12, 13, 8}}, {{8, 11, 13}}, {{8, 11, 12}}}};
8027 {{{12, 10}}, {{12, 9}}, {{12, 11}}, {{12, 8}}, {{9, 13}}, {{11, 13}},
8028 {{8, 13}}, {{10, 13}}, {{10, 9}}, {{9, 11}}, {{11, 8}}, {{8, 10}},
8029 {{12, 9}}, {{12, 11}}, {{11, 9}}, {{12, 8}}, {{12, 11}}, {{11, 8}},
8030 {{12, 8}}, {{12, 10}}, {{10, 8}}, {{12, 10}}, {{12, 9}}, {{9, 10}},
8031 {{13, 10}}, {{13, 8}}, {{8, 10}}, {{13, 10}}, {{13, 9}}, {{9, 10}},
8032 {{13, 11}}, {{13, 9}}, {{11, 9}}, {{13, 11}}, {{13, 8}}, {{11, 8}},
8033 {{12, 13}}, {{8, 10}}, {{8, 13}}, {{10, 13}}, {{8, 12}}, {{10, 12}},
8034 {{12, 13}}, {{10, 9}}, {{10, 13}}, {{9, 13}}, {{10, 12}}, {{9, 12}},
8035 {{12, 13}}, {{9, 11}}, {{9, 13}}, {{11, 13}}, {{9, 12}}, {{11, 12}},
8036 {{12, 13}}, {{11, 8}}, {{11, 13}}, {{8, 13}}, {{11, 12}}, {{8, 12}}}};
8049 vertex_ids_for_new_boundary_edges_3d = {
8050 {{{{{4, 8}}, {{6, 8}}, {{0, 8}}, {{2, 8}}}},
8051 {{{{5, 9}}, {{7, 9}}, {{1, 9}}, {{3, 9}}}},
8052 {{{{4, 10}}, {{5, 10}}, {{0, 10}}, {{1, 10}}}},
8053 {{{{6, 11}}, {{7, 11}}, {{2, 11}}, {{3, 11}}}},
8054 {{{{2, 12}}, {{3, 12}}, {{0, 12}}, {{1, 12}}}},
8055 {{{{6, 13}}, {{7, 13}}, {{4, 13}}, {{5, 13}}}}}};
8057 std::vector<Point<spacedim>>
vertices;
8058 std::vector<CellData<dim>> cells;
8063 std::vector<unsigned int> old_to_new_vertex_indices(
8065 std::vector<unsigned int> face_to_new_vertex_indices(
8075 for (
const auto &cell : ref_tria.cell_iterators())
8079 std::array<unsigned int, dim == 2 ? 9 : 14> local_vertex_indices;
8084 const auto v_global = cell->vertex_index(v);
8086 if (old_to_new_vertex_indices[v_global] ==
8089 old_to_new_vertex_indices[v_global] =
vertices.size();
8090 vertices.push_back(cell->vertex(v));
8094 local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
8098 for (
const auto f : cell->face_indices())
8100 const auto f_global = cell->face_index(f);
8102 if (face_to_new_vertex_indices[f_global] ==
8105 face_to_new_vertex_indices[f_global] =
vertices.size();
8107 cell->face(f)->center(
true));
8111 local_vertex_indices.size());
8112 local_vertex_indices[cell->n_vertices() + f] =
8113 face_to_new_vertex_indices[f_global];
8120 local_vertex_indices.size());
8121 local_vertex_indices[cell->n_vertices() + cell->n_faces()] =
8123 vertices.push_back(cell->center(
true));
8127 const auto add_cell = [&](
const unsigned int struct_dim,
8128 const auto &index_vertices,
8129 const unsigned int material_or_boundary_id,
8133 if (struct_dim < dim &&
8138 if (struct_dim == dim)
8150 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8153 local_vertex_indices.size());
8154 cell_data.vertices[i] =
8155 local_vertex_indices[index_vertices[i]];
8156 cell_data.material_id =
8157 material_or_boundary_id;
8158 cell_data.manifold_id =
8161 cells.push_back(cell_data);
8163 else if (dim == 2 && struct_dim == 1)
8167 boundary_line.
boundary_id = material_or_boundary_id;
8169 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8172 local_vertex_indices.size());
8174 local_vertex_indices[index_vertices[i]];
8178 else if (dim == 3 && struct_dim == 2)
8182 boundary_quad.material_id = material_or_boundary_id;
8184 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8187 local_vertex_indices.size());
8188 boundary_quad.vertices[i] =
8189 local_vertex_indices[index_vertices[i]];
8193 else if (dim == 3 && struct_dim == 1)
8197 boundary_line.
boundary_id = material_or_boundary_id;
8199 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
8202 local_vertex_indices.size());
8204 local_vertex_indices[index_vertices[i]];
8214 const auto material_id_cell = cell->material_id();
8220 const auto manifold_id_cell = cell->manifold_id();
8222 for (
const auto &cell_vertices : table_2D_cell)
8223 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8226 for (
const auto &face_vertices : vertex_ids_for_inner_faces_2d)
8238 const auto manifold_id_cell = cell->manifold_id();
8240 for (
const auto &cell_vertices : vertex_ids_for_cells_3d)
8241 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
8245 for (
const auto &face_vertices : vertex_ids_for_inner_faces_3d)
8253 for (
const auto &edge_vertices : vertex_ids_for_inner_edges_3d)
8263 for (
const auto f : cell->face_indices())
8265 const auto bid = cell->face(f)->boundary_id();
8266 const auto mid = cell->face(f)->manifold_id();
8271 for (
const auto &face_vertices :
8272 vertex_ids_for_boundary_faces_2d[f])
8273 add_cell(1, face_vertices, bid, mid);
8279 for (
const auto &face_vertices :
8280 vertex_ids_for_boundary_faces_3d[f])
8281 add_cell(2, face_vertices, bid, mid);
8284 for (
const auto &edge_vertices :
8285 vertex_ids_for_new_boundary_edges_3d[f])
8286 add_cell(1, edge_vertices, bid, mid);
8296 for (
const auto e : cell->line_indices())
8298 auto edge = cell->line(e);
8304 old_to_new_vertex_indices[edge->vertex_index(0)];
8306 old_to_new_vertex_indices[edge->vertex_index(1)];
8318 for (
const auto i : out_tria.get_manifold_ids())
8320 out_tria.set_manifold(i,
FlatManifold<dim, spacedim>());
8325 template <
int spacedim>
8336 template <
template <
int,
int>
class MeshType,
int dim,
int spacedim>
8340 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8341 typename MeshType<dim, spacedim>::face_iterator>
8343 typename ExtractBoundaryMesh<MeshType, dim, spacedim>::return_type
8345 extract_boundary_mesh(
const MeshType<dim, spacedim> &volume_mesh,
8346 MeshType<dim - 1, spacedim> &surface_mesh,
8347 const std::set<types::boundary_id> &boundary_ids)
8360 const unsigned int boundary_dim = dim - 1;
8367 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>
8368 temporary_mapping_level0;
8373 std::vector<bool> touched(
volume_mesh.get_triangulation().n_vertices(),
8377 std::vector<CellData<boundary_dim>> cells;
8379 std::vector<Point<spacedim>>
vertices;
8382 std::map<unsigned int, unsigned int> map_vert_index;
8395 for (
unsigned int i1 = 0; i1 < GeometryInfo<spacedim>::faces_per_cell; ++i1)
8397 for (
unsigned int i2 = 0; i2 <
GeometryInfo<dim - 1>::vertices_per_cell;
8399 swap_matrix[i1][i2] = i2;
8405 std::swap(swap_matrix[0][1], swap_matrix[0][2]);
8406 std::swap(swap_matrix[2][1], swap_matrix[2][2]);
8407 std::swap(swap_matrix[4][1], swap_matrix[4][2]);
8411 std::swap(swap_matrix[1][0], swap_matrix[1][1]);
8412 std::swap(swap_matrix[2][0], swap_matrix[2][1]);
8417 for (
typename MeshType<dim, spacedim>::cell_iterator cell =
8421 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
8423 const typename MeshType<dim, spacedim>::face_iterator face =
8426 if (face->at_boundary() &&
8432 for (
const unsigned int j :
8435 const unsigned int v_index = face->vertex_index(j);
8437 if (!touched[v_index])
8439 vertices.push_back(face->vertex(j));
8440 map_vert_index[v_index] =
vertices.size() - 1;
8441 touched[v_index] =
true;
8444 c_data.
vertices[swap_matrix[i][j]] = map_vert_index[v_index];
8458 for (
unsigned int e = 0;
e < 4; ++
e)
8464 bool edge_found =
false;
8465 for (
auto &boundary_line : subcell_data.boundary_lines)
8467 map_vert_index[face->line(
e)->vertex_index(0)]) &&
8469 map_vert_index[face->line(
e)->vertex_index(
8472 map_vert_index[face->line(
e)->vertex_index(1)]) &&
8474 map_vert_index[face->line(
e)->vertex_index(0)])))
8481 if (edge_found ==
true)
8488 map_vert_index[face->line(e)->vertex_index(0)];
8490 map_vert_index[face->line(e)->vertex_index(1)];
8497 cells.push_back(c_data);
8498 temporary_mapping_level0.push_back(std::make_pair(face, i));
8511 for (
const auto &cell :
surface_mesh.active_cell_iterators())
8512 for (unsigned
int vertex = 0; vertex < 2; ++vertex)
8513 if (cell->face(vertex)->at_boundary())
8514 cell->face(vertex)->set_boundary_id(0);
8522 std::vector<std::pair<
8523 const typename MeshType<dim - 1, spacedim>::cell_iterator,
8524 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>>
8525 temporary_map_boundary_cell_face;
8526 for (
const auto &cell :
surface_mesh.active_cell_iterators())
8527 temporary_map_boundary_cell_face.push_back(
8528 std::make_pair(cell, temporary_mapping_level0.at(cell->
index())));
8540 unsigned int index_cells_deepest_level = 0;
8543 bool changed =
false;
8547 std::vector<unsigned int> cells_refined;
8550 for (
unsigned int cell_n = index_cells_deepest_level;
8551 cell_n < temporary_map_boundary_cell_face.size();
8556 if (temporary_map_boundary_cell_face[cell_n]
8557 .
second.first->has_children())
8561 Assert(temporary_map_boundary_cell_face[cell_n]
8562 .
second.first->refinement_case() ==
8565 temporary_map_boundary_cell_face[cell_n]
8566 .first->set_refine_flag();
8567 cells_refined.push_back(cell_n);
8579 .execute_coarsening_and_refinement();
8582 index_cells_deepest_level = temporary_map_boundary_cell_face.size();
8583 for (
const auto &refined_cell_n : cells_refined)
8585 const typename MeshType<dim - 1, spacedim>::cell_iterator
8587 temporary_map_boundary_cell_face[refined_cell_n].first;
8588 const typename MeshType<dim,
8589 spacedim>::face_iterator refined_face =
8590 temporary_map_boundary_cell_face[refined_cell_n].second.first;
8591 const unsigned int refined_face_number =
8592 temporary_map_boundary_cell_face[refined_cell_n]
8594 for (
unsigned int child_n = 0;
8595 child_n < refined_cell->n_children();
8600 temporary_map_boundary_cell_face.push_back(
8601 std::make_pair(refined_cell->child(
8602 swap_matrix[refined_face_number][child_n]),
8603 std::make_pair(refined_face->child(child_n),
8604 refined_face_number)));
8614 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8615 typename MeshType<dim, spacedim>::face_iterator>
8616 surface_to_volume_mapping;
8617 for (
unsigned int i = 0; i < temporary_map_boundary_cell_face.size(); ++i)
8618 surface_to_volume_mapping[temporary_map_boundary_cell_face[i].
first] =
8619 temporary_map_boundary_cell_face[i].
second.first;
8622 const auto attached_mids =
8624 for (
const auto i :
volume_mesh.get_triangulation().get_manifold_ids())
8626 std::find(attached_mids.
begin(), attached_mids.
end(), i) ==
8627 attached_mids.
end())
8632 return surface_to_volume_mapping;
8637 template <
int dim,
int spacedim>
8641 const std::vector<unsigned int> &repetitions,
8648 std::vector<Point<spacedim>>
vertices;
8649 std::vector<CellData<dim>> cells;
8655 (p2[1] - p1[1]) / repetitions[1]);
8658 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
8659 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
8664 for (
unsigned int j = 0; j < repetitions[1]; ++j)
8665 for (
unsigned int i = 0; i < repetitions[0]; ++i)
8668 std::array<unsigned int, 4> quad{{
8669 (j + 0) * (repetitions[0] + 1) + i + 0,
8670 (j + 0) * (repetitions[0] + 1) + i + 1,
8671 (j + 1) * (repetitions[0] + 1) + i + 0,
8672 (j + 1) * (repetitions[0] + 1) + i + 1
8678 tri.
vertices = {quad[0], quad[1], quad[2]};
8679 cells.push_back(tri);
8685 tri.
vertices = {quad[3], quad[2], quad[1]};
8686 cells.push_back(tri);
8694 (p2[1] - p1[1]) / repetitions[1],
8695 (p2[2] - p1[2]) / repetitions[2]);
8698 for (
unsigned int k = 0; k <= repetitions[2]; ++k)
8699 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
8700 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
8703 p1[2] + dx[2] * k));
8706 for (
unsigned int k = 0; k < repetitions[2]; ++k)
8707 for (
unsigned int j = 0; j < repetitions[1]; ++j)
8708 for (
unsigned int i = 0; i < repetitions[0]; ++i)
8711 std::array<unsigned int, 8> quad{
8712 {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8713 (j + 0) * (repetitions[0] + 1) + i + 0,
8714 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8715 (j + 0) * (repetitions[0] + 1) + i + 1,
8716 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8717 (j + 1) * (repetitions[0] + 1) + i + 0,
8718 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8719 (j + 1) * (repetitions[0] + 1) + i + 1,
8720 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8721 (j + 0) * (repetitions[0] + 1) + i + 0,
8722 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8723 (j + 0) * (repetitions[0] + 1) + i + 1,
8724 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8725 (j + 1) * (repetitions[0] + 1) + i + 0,
8726 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8727 (j + 1) * (repetitions[0] + 1) + i + 1}};
8732 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8733 cell.
vertices = {{quad[0], quad[1], quad[2], quad[4]}};
8735 cell.
vertices = {{quad[0], quad[1], quad[3], quad[5]}};
8737 cells.push_back(cell);
8743 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8744 cell.
vertices = {{quad[2], quad[1], quad[3], quad[7]}};
8746 cell.
vertices = {{quad[0], quad[3], quad[2], quad[6]}};
8747 cells.push_back(cell);
8753 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8754 cell.
vertices = {{quad[1], quad[4], quad[5], quad[7]}};
8756 cell.
vertices = {{quad[0], quad[4], quad[5], quad[6]}};
8757 cells.push_back(cell);
8763 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8764 cell.
vertices = {{quad[2], quad[4], quad[7], quad[6]}};
8766 cell.
vertices = {{quad[3], quad[5], quad[7], quad[6]}};
8767 cells.push_back(cell);
8773 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8774 cell.
vertices = {{quad[1], quad[2], quad[4], quad[7]}};
8776 cell.
vertices = {{quad[0], quad[3], quad[6], quad[5]}};
8777 cells.push_back(cell);
8799 double epsilon = std::numeric_limits<double>::max();
8800 for (
unsigned int i = 0; i < dim; ++i)
8802 0.01 * (
std::abs(p2[i] - p1[i]) / repetitions[i]));
8805 "The distance between corner points must be positive."));
8809 colorize_subdivided_hyper_rectangle(tria, p1, p2, epsilon);
8815 template <
int dim,
int spacedim>
8818 const unsigned int repetitions,
8826 tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2},
colorize);
8832 {{repetitions, repetitions, repetitions}},
8845# include "grid_generator.inst"
void make_mapping(const MeshType &source_grid, const MeshType &destination_grid)
void add_parameter(const std::string &entry, ParameterType ¶meter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
const Point< spacedim > center
numbers::NumberTraits< Number >::real_type norm() const
void initialize(const Triangulation< dim, spacedim > &triangulation)
virtual void add_periodicity(const std::vector< GridTools::PeriodicFacePair< cell_iterator > > &)
virtual types::global_cell_index n_global_active_cells() const
virtual void copy_triangulation(const Triangulation< dim, spacedim > &other_tria)
unsigned int n_faces() const
bool all_reference_cells_are_hyper_cube() const
void save_user_flags_line(std::ostream &out) const
face_iterator end_face() const
cell_iterator begin(const unsigned int level=0) const
virtual void create_triangulation(const std::vector< Point< spacedim > > &vertices, const std::vector< CellData< dim > > &cells, const SubCellData &subcelldata)
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
const std::vector< Point< spacedim > > & get_vertices() const
unsigned int n_active_lines() const
unsigned int n_levels() const
cell_iterator end() const
virtual bool has_hanging_nodes() const
vertex_iterator begin_vertex() const
vertex_iterator end_vertex() const
virtual void execute_coarsening_and_refinement()
virtual unsigned int n_global_levels() const
cell_iterator last() const
face_iterator begin_face() const
unsigned int n_cells() const
void save_user_flags_quad(std::ostream &out) const
unsigned int n_vertices() const
active_cell_iterator begin_active(const unsigned int level=0) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
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)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
Expression fabs(const Expression &x)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void subdivided_hyper_cube_with_simplices(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double p1=0.0, const double p2=1.0, const bool colorize=false)
void parallelepiped(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void hyper_cross(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &sizes, const bool colorize_cells=false)
A center cell with stacks of cell protruding from each surface.
void hyper_ball_balanced(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > ¢er=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
void enclosed_hyper_cube(Triangulation< dim > &tria, const double left=0., const double right=1., const double thickness=1., const bool colorize=false)
void replicate_triangulation(const Triangulation< dim, spacedim > &input, const std::vector< unsigned int > &extents, Triangulation< dim, spacedim > &result)
Replicate a given triangulation in multiple coordinate axes.
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void general_cell(Triangulation< dim, spacedim > &tria, const std::vector< Point< spacedim > > &vertices, const bool colorize=false)
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)
void hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void eccentric_hyper_shell(Triangulation< dim > &triangulation, const Point< dim > &inner_center, const Point< dim > &outer_center, const double inner_radius, const double outer_radius, const unsigned int n_cells)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void moebius(Triangulation< 3, 3 > &tria, const unsigned int n_cells, const unsigned int n_rotations, const double R, const double r)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void half_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void quarter_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
Rectangular domain with rectangular pattern of holes.
void create_union_triangulation(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result)
void 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)
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.)
spacedim MeshType< dim - 1, spacedim > const std::set< types::boundary_id > & boundary_ids
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
void subdivided_hyper_L(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &bottom_left, const Point< dim > &top_right, const std::vector< int > &n_cells_to_remove)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=Point< spacedim >(), const double radius=1.)
void concentric_hyper_shells(Triangulation< dim > &triangulation, const Point< dim > ¢er, const double inner_radius=0.125, const double outer_radius=0.25, const unsigned int n_shells=1, const double skewness=0.1, const unsigned int n_cells_per_shell=0, const bool colorize=false)
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
spacedim MeshType< dim - 1, spacedim > & surface_mesh
void hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, 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 hyper_cube_with_cylindrical_hole(Triangulation< dim > &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 truncated_cone(Triangulation< dim > &tria, const double radius_0=1.0, const double radius_1=0.5, const double half_length=1.0)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cylinder_shell(Triangulation< dim > &tria, const double length, const double inner_radius, const double outer_radius, const unsigned int n_radial_cells=0, const unsigned int n_axial_cells=0, const bool colorize=false)
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
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)
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)
const types::material_id invalid_material_id
static constexpr double PI_2
const types::boundary_id invalid_boundary_id
static constexpr double PI
const types::boundary_id internal_face_boundary_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
const Iterator const std_cxx20::type_identity_t< Iterator > & end
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::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)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
typename internal::ndarray::HelperArray< T, Ns... >::type ndarray
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< unsigned int > vertices
types::manifold_id manifold_id
types::material_id material_id
types::boundary_id boundary_id
static unsigned int face_to_cell_vertices(const unsigned int face, const unsigned int vertex, const bool face_orientation=true, const bool face_flip=false, const bool face_rotation=false)
std::vector< CellData< 2 > > boundary_quads
std::vector< CellData< 1 > > boundary_lines
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)