52 : airfoil_type(
"NACA")
54 , joukowski_center(-0.1, 0.14)
60 , incline_factor(0.35)
63 , n_subdivision_x_0(3)
64 , n_subdivision_x_1(2)
65 , n_subdivision_x_2(5)
67 , airfoil_sampling_factor(2)
70 airfoil_length <= height,
72 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
74 Assert(incline_factor < 1.0 && incline_factor >= 0.0,
75 ExcMessage(
"incline_factor has to be in [0,1)!"));
88 "Mesh height measured from airfoil nose to horizontal boundaries");
92 "Length measured from airfoil leading edge to vertical outlet boundary");
96 "Define obliqueness of the vertical mesh around the airfoil");
105 "Type of airfoil geometry, either NACA or Joukowski airfoil",
120 "Joukowski circle center coordinates");
123 "Joukowski airfoil length leading to trailing edge");
131 "Number of global refinements");
133 "NumberSubdivisionX0",
135 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
137 "NumberSubdivisionX1",
139 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
141 "NumberSubdivisionX2",
143 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
146 "Number of subdivisions normal to airfoil");
150 "Factor to obtain a finer mesh at the airfoil surface");
165 static const unsigned int id_block_1 = 1;
166 static const unsigned int id_block_2 = 2;
167 static const unsigned int id_block_3 = 3;
168 static const unsigned int id_block_4 = 4;
169 static const unsigned int id_block_5 = 5;
170 static const unsigned int id_block_6 = 6;
175 MeshGenerator(
const AdditionalData &data)
176 : refinements(data.refinements)
177 , n_subdivision_x_0(data.n_subdivision_x_0)
178 , n_subdivision_x_1(data.n_subdivision_x_1)
179 , n_subdivision_x_2(data.n_subdivision_x_2)
180 , n_subdivision_y(data.n_subdivision_y)
181 , height(data.height)
182 , length_b2(data.length_b2)
183 , incline_factor(data.incline_factor)
184 , bias_factor(data.bias_factor)
186 , n_cells_x_0(
Utilities::
pow(2, refinements) * n_subdivision_x_0)
187 , n_cells_x_1(
Utilities::
pow(2, refinements) * n_subdivision_x_1)
188 , n_cells_x_2(
Utilities::
pow(2, refinements) * n_subdivision_x_2)
189 , n_cells_y(
Utilities::
pow(2, refinements) * n_subdivision_y)
190 , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
192 , airfoil_1D(set_airfoil_length(
194 data.airfoil_type ==
"Joukowski" ?
195 joukowski(data.joukowski_center,
196 n_points_on_each_side,
197 data.airfoil_sampling_factor) :
198 (data.airfoil_type ==
"NACA" ?
200 n_points_on_each_side,
201 data.airfoil_sampling_factor) :
202 std::array<std::vector<
Point<2>>, 2>{
204 std::vector<Point<2>>{
208 data.airfoil_length))
209 , end_b0_x_u(airfoil_1D[0][n_cells_x_0](0))
210 , end_b0_x_l(airfoil_1D[1][n_cells_x_0](0))
211 , nose_x(airfoil_1D[0].front()(0))
212 , tail_x(airfoil_1D[0].back()(0))
213 , tail_y(airfoil_1D[0].back()(1))
214 , center_mesh(0.5 * std::abs(end_b0_x_u + end_b0_x_l))
215 , length_b1_x(tail_x - center_mesh)
217 (edge_length + std::abs(nose_x - center_mesh))))
221 ,
A(nose_x - edge_length, 0)
223 ,
C(center_mesh, +std::abs(nose_x - center_mesh) * std::tan(
gamma))
224 , D(center_mesh, height)
225 ,
E(center_mesh, -std::abs(nose_x - center_mesh) * std::tan(
gamma))
226 ,
F(center_mesh, -height)
230 , J(tail_x + length_b2, 0)
234 Assert(data.airfoil_type ==
"Joukowski" ||
235 data.airfoil_type ==
"NACA",
248 make_coarse_grid(tria_grid);
250 set_boundary_ids(tria_grid);
252 if (periodic_faces !=
nullptr)
255 tria_grid, 5, 4, 1, *periodic_faces);
273 (void)periodic_faces;
280 const unsigned int refinements;
283 const unsigned int n_subdivision_x_0;
286 const unsigned int n_subdivision_x_1;
289 const unsigned int n_subdivision_x_2;
293 const unsigned int n_subdivision_y;
300 const double length_b2;
304 const double incline_factor;
307 const double bias_factor;
310 const double edge_length;
313 const unsigned int n_cells_x_0;
316 const unsigned int n_cells_x_1;
319 const unsigned int n_cells_x_2;
323 const unsigned int n_cells_y;
326 const unsigned int n_points_on_each_side;
330 const std::array<std::vector<Point<2>>, 2> airfoil_1D;
334 const double end_b0_x_u;
338 const double end_b0_x_l;
352 const double center_mesh;
355 const double length_b1_x;
382 const Point<2> A, B,
C, D,
E,
F, G, H, I, J,
K,
L;
421 static std::array<std::vector<Point<2>>, 2>
422 joukowski(
const Point<2> & centerpoint,
423 const unsigned int number_points,
424 const unsigned int factor)
426 std::array<std::vector<Point<2>>, 2> airfoil_1D;
427 const unsigned int total_points = 2 * number_points - 2;
428 const unsigned int n_airfoilpoints = factor * total_points;
430 const auto jouk_points =
431 joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
434 std::vector<Point<2>> upper_points;
435 std::vector<Point<2>> lower_points;
439 unsigned int nose_index = 0;
440 unsigned int tail_index = 0;
441 double nose_x_coordinate = 0;
442 double tail_x_coordinate = 0;
446 for (
unsigned int i = 0; i < jouk_points.size(); ++i)
448 if (jouk_points[i](0) < nose_x_coordinate)
450 nose_x_coordinate = jouk_points[i](0);
453 if (jouk_points[i](0) > tail_x_coordinate)
455 tail_x_coordinate = jouk_points[i](0);
461 for (
unsigned int i = tail_index; i < jouk_points.size(); ++i)
462 upper_points.emplace_back(jouk_points[i]);
463 for (
unsigned int i = 0; i <= nose_index; ++i)
464 upper_points.emplace_back(jouk_points[i]);
465 std::reverse(upper_points.begin(), upper_points.end());
468 lower_points.insert(lower_points.end(),
469 jouk_points.begin() + nose_index,
470 jouk_points.begin() + tail_index + 1);
473 airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
474 airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
477 auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
478 const double nose_x_pos = vector.front()(0);
479 for (
auto &i : vector)
483 move_nose_to_origin(airfoil_1D[1]);
484 move_nose_to_origin(airfoil_1D[0]);
513 static std::vector<Point<2>>
515 const unsigned int number_points)
517 std::vector<Point<2>> circle_points;
531 radius_test < radius,
533 "Error creating lower circle: Circle for Joukowski-transform does"
534 " not enclose point zeta = -1! Choose different center "
538 const double theta = 2 *
numbers::PI / number_points;
540 for (
unsigned int i = 0; i < number_points; ++i)
541 circle_points.emplace_back(
center[0] - radius *
cos(i * theta),
544 return circle_points;
555 static std::vector<Point<2>>
556 joukowski_transform(
const std::vector<
Point<2>> &circle_points)
558 std::vector<Point<2>> joukowski_points(circle_points.size());
561 for (
unsigned int i = 0; i < circle_points.size(); ++i)
563 const double chi = circle_points[i](0);
564 const double eta = circle_points[i](1);
565 const std::complex<double> zeta(chi, eta);
566 const std::complex<double> z = zeta + 1. / zeta;
568 joukowski_points[i] = {real(z), imag(z)};
570 return joukowski_points;
589 static std::array<std::vector<Point<2>>, 2>
590 naca(
const std::string &serialnumber,
591 const unsigned int number_points,
592 const unsigned int factor)
596 const unsigned int n_airfoilpoints = factor * number_points;
599 return {{make_points_equidistant(
600 naca_create_points(serialnumber, n_airfoilpoints,
true),
602 make_points_equidistant(
603 naca_create_points(serialnumber, n_airfoilpoints,
false),
618 static std::vector<Point<2>>
619 naca_create_points(
const std::string &serialnumber,
620 const unsigned int number_points,
623 Assert(serialnumber.size() == 4,
624 ExcMessage(
"This NACA-serial number is not implemented!"));
626 return naca_create_points_4_digits(serialnumber,
645 static std::vector<Point<2>>
646 naca_create_points_4_digits(
const std::string &serialnumber,
647 const unsigned int number_points,
651 const unsigned int digit_0 = (serialnumber[0] -
'0');
652 const unsigned int digit_1 = (serialnumber[1] -
'0');
653 const unsigned int digit_2 = (serialnumber[2] -
'0');
654 const unsigned int digit_3 = (serialnumber[3] -
'0');
656 const unsigned int digit_23 = 10 * digit_2 + digit_3;
659 const double t =
static_cast<double>(digit_23) / 100.0;
661 std::vector<Point<2>> naca_points;
663 if (digit_0 == 0 && digit_1 == 0)
664 for (
unsigned int i = 0; i < number_points; ++i)
666 const double x = i * 1 / (1.0 * number_points - 1);
669 (0.2969 *
std::pow(x, 0.5) - 0.126 * x -
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);
686 (x <= p) ? m / std::pow(p, 2) * (2 * p * x -
std::pow(x, 2)) :
687 m / std::pow(1 - p, 2) *
688 ((1 - 2 * p) + 2 * p * x - std::pow(x, 2));
690 const double dy_c = (x <= p) ?
691 2 * m / std::pow(p, 2) * (p - x) :
692 2 * m / std::pow(1 - p, 2) * (p - x);
696 (0.2969 *
std::pow(x, 0.5) - 0.126 * x -
703 naca_points.emplace_back(x - y_t * std::sin(theta),
704 y_c + y_t * std::cos(theta));
706 naca_points.emplace_back(x + y_t * std::sin(theta),
707 y_c - y_t * std::cos(theta));
723 static std::array<std::vector<Point<2>>, 2>
724 set_airfoil_length(
const std::array<std::vector<
Point<2>>, 2> &input,
725 const double desired_len)
727 std::array<std::vector<Point<2>>, 2> output;
728 output[0] = set_airfoil_length(input[0], desired_len);
729 output[1] = set_airfoil_length(input[1], desired_len);
741 static std::vector<Point<2>>
742 set_airfoil_length(
const std::vector<
Point<2>> &input,
743 const double desired_len)
745 std::vector<Point<2>> output = input;
748 desired_len / input.front().distance(input.back());
750 for (
auto &x : output)
766 static std::vector<Point<2>>
767 make_points_equidistant(
768 const std::vector<
Point<2>> &non_equidistant_points,
769 const unsigned int number_points)
771 const unsigned int n_points =
772 non_equidistant_points
776 std::vector<double> arclength_L(non_equidistant_points.size(), 0);
777 for (
unsigned int i = 0; i < non_equidistant_points.size() - 1; ++i)
780 non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
783 const auto airfoil_length =
785 const auto deltaX = airfoil_length / (number_points - 1);
789 std::vector<Point<2>> equidist(
792 equidist[0] = non_equidistant_points[0];
793 equidist[number_points - 1] = non_equidistant_points[n_points - 1];
797 for (
unsigned int j = 0, i = 1; j < n_points - 1; ++j)
800 const auto Lj = arclength_L[j];
801 const auto Ljp = arclength_L[j + 1];
803 while (Lj <= i * deltaX && i * deltaX <= Ljp &&
804 i < number_points - 1)
806 equidist[i] =
Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
807 (non_equidistant_points[j + 1] -
808 non_equidistant_points[j]) +
809 non_equidistant_points[j]);
829 std::vector<Triangulation<2>> trias(10);
833 const std::vector<Point<2>> & corner_vertices,
834 const std::vector<unsigned int> &repetitions,
846 auto &
point = it->vertex();
847 const double xi =
point(0);
848 const double eta =
point(1);
851 point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
852 (1 + xi) * (1 - eta) * corner_vertices[1] +
853 (1 - xi) * (1 + eta) * corner_vertices[2] +
854 (1 + xi) * (1 + eta) * corner_vertices[3]);
866 {n_subdivision_y, n_subdivision_x_0},
870 {n_subdivision_y, n_subdivision_x_0},
874 {n_subdivision_x_1, n_subdivision_y},
878 {n_subdivision_x_1, n_subdivision_y},
882 {n_subdivision_x_2, n_subdivision_y},
886 {n_subdivision_x_2, n_subdivision_y},
916 if (cell->face(f)->at_boundary() ==
false)
919 const auto mid = cell->material_id();
921 if ((mid == id_block_1 && f == 0) ||
922 (mid == id_block_4 && f == 0))
923 cell->face(f)->set_boundary_id(0);
924 else if ((mid == id_block_3 && f == 0) ||
925 (mid == id_block_6 && f == 2))
926 cell->face(f)->set_boundary_id(1);
927 else if ((mid == id_block_1 && f == 1) ||
928 (mid == id_block_2 && f == 1))
929 cell->face(f)->set_boundary_id(2);
930 else if ((mid == id_block_4 && f == 1) ||
931 (mid == id_block_5 && f == 3))
932 cell->face(f)->set_boundary_id(3);
933 else if ((mid == id_block_2 && f == 0) ||
934 (mid == id_block_3 && f == 2))
935 cell->face(f)->set_boundary_id(4);
936 else if ((mid == id_block_5 && f == 2) ||
937 (mid == id_block_6 && f == 0))
938 cell->face(f)->set_boundary_id(5);
972 std::sin(
gamma) * edge_length);
975 for (
auto &cell :
tria)
979 if (vertex_processed[cell.vertex_index(v)])
983 vertex_processed[cell.vertex_index(v)] =
true;
985 auto &node = cell.vertex(v);
988 if (cell.material_id() == id_block_1 ||
989 cell.material_id() == id_block_4)
998 if (cell.material_id() == id_block_1)
1001 (node - horizontal_offset) +
1007 else if (cell.material_id() == id_block_4)
1010 (node - horizontal_offset) -
1016 const double trapeze_height =
1022 const double x2 =
L - l_a - l_b;
1024 const double Dx = x1 + x2 + x3;
1025 const double deltax =
1027 const double dx = Dx / n_cells_x_0;
1028 const double dy = trapeze_height / n_cells_y;
1030 static_cast<int>(std::round((node_(0) - deltax) /
dx));
1032 static_cast<int>(std::round(std::abs(node_(1)) / dy));
1034 node_(0) =
numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1035 node_(1) = height * (1.0 * iy) / n_cells_y;
1042 const double dy = height / n_cells_y;
1044 static_cast<int>(std::round(node_(0) /
dx));
1046 static_cast<int>(std::round(node_(1) / dy));
1047 const double alpha =
1048 bias_alpha(1 - (1.0 * iy) / n_cells_y);
1049 const double theta = node_(0);
1050 const Point<2> p(-height * std::cos(theta) + center_mesh,
1051 ((cell.material_id() == id_block_1) ?
1057 (cell.material_id() == id_block_1) ? (0) : (1))][ix] *
1062 else if (cell.material_id() == id_block_2 ||
1063 cell.material_id() == id_block_5)
1067 (std::abs(D(1) -
C(1)) == std::abs(
F(1) -
E(1))) &&
1068 (std::abs(
C(1)) == std::abs(
E(1))) &&
1069 (std::abs(G(1)) == std::abs(I(1))),
1071 "Points D,C,G and E,F,I are not defined symmetric to "
1072 "x-axis, which is required to interpolate block 2"
1073 " and 5 with same geometric computations."));
1074 const double l_y = D(1) -
C(1);
1075 const double l_h = D(1) - l_y;
1076 const double by = -l_h / length_b1_x * (node(0) - H(0));
1077 const double dy = (height - by) / n_cells_y;
1078 const int iy =
static_cast<int>(
1079 std::round((std::abs(node(1)) - by) / dy));
1080 const double dx = length_b1_x / n_cells_x_1;
1081 const int ix =
static_cast<int>(
1082 std::round(std::abs(node(0) - center_mesh) /
dx));
1084 const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1089 incline_factor * length_b2 * ix /
1091 ((cell.material_id() == id_block_2) ?
1097 (cell.material_id() == id_block_2) ? (0) : (1))]
1098 [n_cells_x_0 + ix] *
1102 else if (cell.material_id() == id_block_3 ||
1103 cell.material_id() == id_block_6)
1106 const double dx = length_b2 / n_cells_x_2;
1107 const double dy = height / n_cells_y;
1108 const int ix =
static_cast<int>(
1109 std::round(std::abs(node(0) - H(0)) /
dx));
1111 static_cast<int>(std::round(std::abs(node(1)) / dy));
1113 const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1114 const double alpha_x =
1115 bias_alpha(1 - (
static_cast<double>(ix)) / n_cells_x_2);
1119 const Point<2> p1(J(0) - (1 - incline_factor) * length_b2 *
1121 ((cell.material_id() == id_block_3) ?
1126 const Point<2> p2(J(0) - alpha_x * length_b2, tail_y);
1127 node = p1 * (1 - alpha_y) + p2 * alpha_y;
1147 bias_alpha(
double alpha)
const
1157 internal_create_triangulation(
1161 const AdditionalData & additional_data)
1163 MeshGenerator mesh_generator(additional_data);
1166 if (
auto parallel_tria =
1169 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1170 else if (
auto parallel_tria =
dynamic_cast<
1173 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1175 mesh_generator.create_triangulation(
tria, periodic_faces);
1192 const AdditionalData &)
1202 const AdditionalData &additional_data)
1204 internal_create_triangulation(
tria,
nullptr, additional_data);
1215 const AdditionalData & additional_data)
1217 internal_create_triangulation(
tria, &periodic_faces, additional_data);
1228 const AdditionalData & additional_data)
1232 (void)additional_data;
1233 (void)periodic_faces;
1244 template <
int dim,
int spacedim>
1256 cell->face(f)->set_boundary_id(f);
1262 template <
int spacedim>
1273 if (cell->center()(0) > 0)
1274 cell->set_material_id(1);
1281 template <
int dim,
int spacedim>
1299 for (; face != endface; ++face)
1300 if (face->at_boundary())
1301 if (face->boundary_id() == 0)
1306 face->set_boundary_id(0);
1308 face->set_boundary_id(1);
1310 face->set_boundary_id(2);
1312 face->set_boundary_id(3);
1314 face->set_boundary_id(4);
1316 face->set_boundary_id(5);
1328 for (
unsigned int d = 0;
d < dim; ++
d)
1329 if (cell->center()(
d) > 0)
1331 cell->set_material_id(
id);
1357 cell->face(2)->set_all_boundary_ids(1);
1381 cell->face(4)->set_all_boundary_ids(1);
1385 cell->face(2)->set_all_boundary_ids(1);
1389 cell->face(2)->set_all_boundary_ids(1);
1393 cell->face(0)->set_all_boundary_ids(1);
1397 cell->face(2)->set_all_boundary_ids(1);
1401 cell->face(0)->set_all_boundary_ids(1);
1412 cell->face(5)->set_all_boundary_ids(1);
1423 unsigned int count = 0;
1426 if (cell->face(5)->at_boundary())
1428 cell->face(5)->set_all_boundary_ids(1);
1449 const double inner_radius,
1450 const double outer_radius)
1455 double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
1456 double eps = 1
e-3 * middle;
1459 for (; cell !=
tria.
end(); ++cell)
1462 if (!cell->face(f)->at_boundary())
1465 double radius = cell->face(f)->center().norm() -
center.
norm();
1466 if (
std::fabs(cell->face(f)->center()(0)) <
1469 cell->face(f)->set_boundary_id(2);
1470 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1472 if (cell->face(f)->line(j)->at_boundary())
1473 if (
std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1474 cell->face(f)->line(j)->vertex(1).norm()) >
1476 cell->face(f)->line(j)->set_boundary_id(2);
1478 else if (
std::fabs(cell->face(f)->center()(1)) <
1481 cell->face(f)->set_boundary_id(3);
1482 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1484 if (cell->face(f)->line(j)->at_boundary())
1485 if (
std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1486 cell->face(f)->line(j)->vertex(1).norm()) >
1488 cell->face(f)->line(j)->set_boundary_id(3);
1490 else if (
std::fabs(cell->face(f)->center()(2)) <
1493 cell->face(f)->set_boundary_id(4);
1494 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1496 if (cell->face(f)->line(j)->at_boundary())
1497 if (
std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1498 cell->face(f)->line(j)->vertex(1).norm()) >
1500 cell->face(f)->line(j)->set_boundary_id(4);
1502 else if (radius < middle)
1504 cell->face(f)->set_boundary_id(0);
1505 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1507 if (cell->face(f)->line(j)->at_boundary())
1508 if (
std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1509 cell->face(f)->line(j)->vertex(1).norm()) <
1511 cell->face(f)->line(j)->set_boundary_id(0);
1513 else if (radius > middle)
1515 cell->face(f)->set_boundary_id(1);
1516 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1518 if (cell->face(f)->line(j)->at_boundary())
1519 if (
std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1520 cell->face(f)->line(j)->vertex(1).norm()) <
1522 cell->face(f)->line(j)->set_boundary_id(1);
1532 template <
int dim,
int spacedim>
1543 for (
unsigned int i = 0; i < dim; ++i)
1583 std::vector<CellData<dim>> cells(1);
1585 cells[0].vertices[i] = i;
1586 cells[0].material_id = 0;
1592 colorize_hyper_rectangle(
tria);
1597 template <
int dim,
int spacedim>
1605 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1608 for (
unsigned int i = 0; i < dim; ++i)
1628 for (
unsigned int d = 0;
d < dim; ++
d)
1629 for (
unsigned int c = 1; c <= dim; ++c)
1632 ExcMessage(
"Vertices of simplex must form a right handed system"));
1636 std::vector<Point<dim>> points =
vertices;
1640 for (
unsigned int i = 0; i <= dim; ++i)
1642 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1648 for (
unsigned int i = 1; i < dim; ++i)
1649 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1651 for (
unsigned int i = 0; i <= dim; ++i)
1652 points.push_back(1. / 3. *
1653 (points[i] + points[(i + 1) % (dim + 1)] +
1654 points[(i + 2) % (dim + 1)]));
1656 points.push_back((1. / (dim + 1)) *
center);
1658 std::vector<CellData<dim>> cells(dim + 1);
1663 cells[0].vertices[0] = 0;
1664 cells[0].vertices[1] = 3;
1665 cells[0].vertices[2] = 5;
1666 cells[0].vertices[3] = 6;
1667 cells[0].material_id = 0;
1669 cells[1].vertices[0] = 3;
1670 cells[1].vertices[1] = 1;
1671 cells[1].vertices[2] = 6;
1672 cells[1].vertices[3] = 4;
1673 cells[1].material_id = 0;
1675 cells[2].vertices[0] = 5;
1676 cells[2].vertices[1] = 6;
1677 cells[2].vertices[2] = 2;
1678 cells[2].vertices[3] = 4;
1679 cells[2].material_id = 0;
1683 cells[0].vertices[0] = 0;
1684 cells[0].vertices[1] = 4;
1685 cells[0].vertices[2] = 8;
1686 cells[0].vertices[3] = 10;
1687 cells[0].vertices[4] = 7;
1688 cells[0].vertices[5] = 13;
1689 cells[0].vertices[6] = 12;
1690 cells[0].vertices[7] = 14;
1691 cells[0].material_id = 0;
1693 cells[1].vertices[0] = 4;
1694 cells[1].vertices[1] = 1;
1695 cells[1].vertices[2] = 10;
1696 cells[1].vertices[3] = 5;
1697 cells[1].vertices[4] = 13;
1698 cells[1].vertices[5] = 9;
1699 cells[1].vertices[6] = 14;
1700 cells[1].vertices[7] = 11;
1701 cells[1].material_id = 0;
1703 cells[2].vertices[0] = 8;
1704 cells[2].vertices[1] = 10;
1705 cells[2].vertices[2] = 2;
1706 cells[2].vertices[3] = 5;
1707 cells[2].vertices[4] = 12;
1708 cells[2].vertices[5] = 14;
1709 cells[2].vertices[6] = 6;
1710 cells[2].vertices[7] = 11;
1711 cells[2].material_id = 0;
1713 cells[3].vertices[0] = 7;
1714 cells[3].vertices[1] = 13;
1715 cells[3].vertices[2] = 12;
1716 cells[3].vertices[3] = 14;
1717 cells[3].vertices[4] = 3;
1718 cells[3].vertices[5] = 9;
1719 cells[3].vertices[6] = 6;
1720 cells[3].vertices[7] = 11;
1721 cells[3].material_id = 0;
1731 template <
int dim,
int spacedim>
1751 for (
unsigned int d = 0;
d < dim; ++
d)
1760 std::vector<CellData<dim>> cells(1);
1763 cells[0].vertices[v] = v;
1773 const unsigned int n_rotations,
1777 const unsigned int dim = 3;
1780 "More than 4 cells are needed to create a moebius grid."));
1782 ExcMessage(
"Outer and inner radius must be positive."));
1784 ExcMessage(
"Outer radius must be greater than inner radius."));
1791 for (
unsigned int i = 0; i <
n_cells; ++i)
1792 for (
unsigned int j = 0; j < 4; ++j)
1806 unsigned int offset = 0;
1811 static constexpr std::array<unsigned int, 8> local_vertex_numbering{
1812 {0, 1, 5, 4, 2, 3, 7, 6}};
1813 std::vector<CellData<dim>> cells(
n_cells);
1814 for (
unsigned int i = 0; i <
n_cells; ++i)
1816 for (
unsigned int j = 0; j < 2; ++j)
1818 cells[i].vertices[local_vertex_numbering[0 + 4 * j]] =
1820 cells[i].vertices[local_vertex_numbering[1 + 4 * j]] =
1822 cells[i].vertices[local_vertex_numbering[2 + 4 * j]] =
1824 cells[i].vertices[local_vertex_numbering[3 + 4 * j]] =
1828 cells[i].material_id = 0;
1832 cells[
n_cells - 1].vertices[local_vertex_numbering[4]] =
1833 (0 + n_rotations) % 4;
1834 cells[
n_cells - 1].vertices[local_vertex_numbering[5]] =
1835 (3 + n_rotations) % 4;
1836 cells[
n_cells - 1].vertices[local_vertex_numbering[6]] =
1837 (2 + n_rotations) % 4;
1838 cells[
n_cells - 1].vertices[local_vertex_numbering[7]] =
1839 (1 + n_rotations) % 4;
1856 ExcMessage(
"Outer radius R must be greater than the inner "
1860 const unsigned int dim = 2;
1861 const unsigned int spacedim = 3;
1862 std::vector<Point<spacedim>>
vertices(16);
1881 std::vector<CellData<dim>> cells(16);
1883 cells[0].vertices[0] = 0;
1884 cells[0].vertices[1] = 4;
1885 cells[0].vertices[2] = 3;
1886 cells[0].vertices[3] = 7;
1887 cells[0].material_id = 0;
1889 cells[1].vertices[0] = 1;
1890 cells[1].vertices[1] = 5;
1891 cells[1].vertices[2] = 0;
1892 cells[1].vertices[3] = 4;
1893 cells[1].material_id = 0;
1895 cells[2].vertices[0] = 2;
1896 cells[2].vertices[1] = 6;
1897 cells[2].vertices[2] = 1;
1898 cells[2].vertices[3] = 5;
1899 cells[2].material_id = 0;
1901 cells[3].vertices[0] = 3;
1902 cells[3].vertices[1] = 7;
1903 cells[3].vertices[2] = 2;
1904 cells[3].vertices[3] = 6;
1905 cells[3].material_id = 0;
1907 cells[4].vertices[0] = 4;
1908 cells[4].vertices[1] = 8;
1909 cells[4].vertices[2] = 7;
1910 cells[4].vertices[3] = 11;
1911 cells[4].material_id = 0;
1913 cells[5].vertices[0] = 5;
1914 cells[5].vertices[1] = 9;
1915 cells[5].vertices[2] = 4;
1916 cells[5].vertices[3] = 8;
1917 cells[5].material_id = 0;
1919 cells[6].vertices[0] = 6;
1920 cells[6].vertices[1] = 10;
1921 cells[6].vertices[2] = 5;
1922 cells[6].vertices[3] = 9;
1923 cells[6].material_id = 0;
1925 cells[7].vertices[0] = 7;
1926 cells[7].vertices[1] = 11;
1927 cells[7].vertices[2] = 6;
1928 cells[7].vertices[3] = 10;
1929 cells[7].material_id = 0;
1931 cells[8].vertices[0] = 8;
1932 cells[8].vertices[1] = 12;
1933 cells[8].vertices[2] = 11;
1934 cells[8].vertices[3] = 15;
1935 cells[8].material_id = 0;
1937 cells[9].vertices[0] = 9;
1938 cells[9].vertices[1] = 13;
1939 cells[9].vertices[2] = 8;
1940 cells[9].vertices[3] = 12;
1941 cells[9].material_id = 0;
1943 cells[10].vertices[0] = 10;
1944 cells[10].vertices[1] = 14;
1945 cells[10].vertices[2] = 9;
1946 cells[10].vertices[3] = 13;
1947 cells[10].material_id = 0;
1949 cells[11].vertices[0] = 11;
1950 cells[11].vertices[1] = 15;
1951 cells[11].vertices[2] = 10;
1952 cells[11].vertices[3] = 14;
1953 cells[11].material_id = 0;
1955 cells[12].vertices[0] = 12;
1956 cells[12].vertices[1] = 0;
1957 cells[12].vertices[2] = 15;
1958 cells[12].vertices[3] = 3;
1959 cells[12].material_id = 0;
1961 cells[13].vertices[0] = 13;
1962 cells[13].vertices[1] = 1;
1963 cells[13].vertices[2] = 12;
1964 cells[13].vertices[3] = 0;
1965 cells[13].material_id = 0;
1967 cells[14].vertices[0] = 14;
1968 cells[14].vertices[1] = 2;
1969 cells[14].vertices[2] = 13;
1970 cells[14].vertices[3] = 1;
1971 cells[14].material_id = 0;
1973 cells[15].vertices[0] = 15;
1974 cells[15].vertices[1] = 3;
1975 cells[15].vertices[2] = 14;
1976 cells[15].vertices[3] = 2;
1977 cells[15].material_id = 0;
1993 const unsigned int n_cells_toroidal,
1997 ExcMessage(
"Outer radius R must be greater than the inner "
2000 Assert(n_cells_toroidal > 2,
2001 ExcMessage(
"Number of cells in toroidal direction has "
2002 "to be at least 3."));
2008 double const a = 1. / (1 +
std::sqrt(2.0));
2011 const unsigned int additional_layer =
2015 const unsigned int n_point_layers_toroidal =
2016 n_cells_toroidal + additional_layer;
2017 std::vector<Point<3>>
vertices(8 * n_point_layers_toroidal);
2029 double const phi_cell = phi / n_cells_toroidal;
2030 for (
unsigned int c = 1; c < n_point_layers_toroidal; ++c)
2032 for (
unsigned int v = 0; v < 8; ++v)
2034 double const r_2d =
vertices[v][0];
2042 std::vector<CellData<3>> cells(5 * n_cells_toroidal);
2043 for (
unsigned int c = 0; c < n_cells_toroidal; ++c)
2045 for (
unsigned int j = 0; j < 2; ++j)
2047 const unsigned int offset =
2048 (8 * (c + j)) % (8 * n_point_layers_toroidal);
2051 cells[5 * c].vertices[0 + j * 4] = offset + 0;
2052 cells[5 * c].vertices[1 + j * 4] = offset + 1;
2053 cells[5 * c].vertices[2 + j * 4] = offset + 2;
2054 cells[5 * c].vertices[3 + j * 4] = offset + 3;
2056 cells[5 * c + 1].vertices[0 + j * 4] = offset + 2;
2057 cells[5 * c + 1].vertices[1 + j * 4] = offset + 3;
2058 cells[5 * c + 1].vertices[2 + j * 4] = offset + 4;
2059 cells[5 * c + 1].vertices[3 + j * 4] = offset + 5;
2061 cells[5 * c + 2].vertices[0 + j * 4] = offset + 4;
2062 cells[5 * c + 2].vertices[1 + j * 4] = offset + 5;
2063 cells[5 * c + 2].vertices[2 + j * 4] = offset + 6;
2064 cells[5 * c + 2].vertices[3 + j * 4] = offset + 7;
2066 cells[5 * c + 3].vertices[0 + j * 4] = offset + 0;
2067 cells[5 * c + 3].vertices[1 + j * 4] = offset + 2;
2068 cells[5 * c + 3].vertices[2 + j * 4] = offset + 6;
2069 cells[5 * c + 3].vertices[3 + j * 4] = offset + 4;
2071 cells[5 * c + 4].vertices[0 + j * 4] = offset + 3;
2072 cells[5 * c + 4].vertices[1 + j * 4] = offset + 1;
2073 cells[5 * c + 4].vertices[2 + j * 4] = offset + 5;
2074 cells[5 * c + 4].vertices[3 + j * 4] = offset + 7;
2077 cells[5 * c].material_id = 0;
2079 cells[5 * c + 1].material_id = 1;
2080 cells[5 * c + 2].material_id = 0;
2081 cells[5 * c + 3].material_id = 0;
2082 cells[5 * c + 4].material_id = 0;
2097 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2099 cell->face(f)->set_all_manifold_ids(1);
2104 if (cell->material_id() == 1)
2106 cell->set_all_manifold_ids(2);
2108 cell->set_material_id(0);
2123 template <
int dim,
int spacedim>
2144 "The volume of the cell is not greater than zero. "
2145 "This could be due to the wrong ordering of the vertices."));
2176 std::array<Tensor<1, 2>, 2> edges;
2177 edges[0] = corners[0];
2178 edges[1] = corners[1];
2179 std::vector<unsigned int> subdivisions;
2180 subdivided_parallelepiped<2, 2>(
2192 unsigned int n_subdivisions[dim];
2193 for (
unsigned int i = 0; i < dim; ++i)
2194 n_subdivisions[i] = 1;
2203 const unsigned int n_subdivisions,
2209 unsigned int n_subdivisions_[dim];
2210 for (
unsigned int i = 0; i < dim; ++i)
2211 n_subdivisions_[i] = n_subdivisions;
2221 const unsigned int (&n_subdivisions)[dim],
2223 const unsigned int *n_subdivisions,
2229 std::vector<unsigned int> subdivisions;
2230 std::array<Tensor<1, dim>, dim> edges;
2231 for (
unsigned int i = 0; i < dim; ++i)
2233 subdivisions.push_back(n_subdivisions[i]);
2234 edges[i] = corners[i];
2237 subdivided_parallelepiped<dim, dim>(
2244 template <
int dim,
int spacedim>
2249 const std::vector<unsigned int> &subdivisions,
2252 std::vector<unsigned int> compute_subdivisions = subdivisions;
2253 if (compute_subdivisions.size() == 0)
2255 compute_subdivisions.resize(dim, 1);
2258 Assert(compute_subdivisions.size() == dim,
2259 ExcMessage(
"One subdivision must be provided for each dimension."));
2261 for (
unsigned int i = 0; i < dim; ++i)
2263 Assert(compute_subdivisions[i] > 0,
2266 edges[i].
norm() > 0,
2268 "Edges in subdivided_parallelepiped() must not be degenerate."));
2276 bool twisted_data =
false;
2281 twisted_data = (edges[0][0] < 0);
2288 const double plane_normal =
2289 edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2290 twisted_data = (plane_normal < 0.0);
2298 Assert(std::abs(edges[0] * edges[1] /
2299 (edges[0].
norm() * edges[1].
norm()) -
2302 "Edges in subdivided_parallelepiped() must point in"
2303 " different directions."));
2305 cross_product_3d(edges[0], edges[1]);
2319 twisted_data = (plane_normal * edges[2] < 0.0);
2329 "The triangulation you are trying to create will consist of cells"
2330 " with negative measures. This is usually the result of input data"
2331 " that does not define a right-handed coordinate system. The usual"
2332 " fix for this is to ensure that in 1D the given point is to the"
2333 " right of the origin (or the given edge tensor is positive), in 2D"
2334 " that the two edges (and their cross product) obey the right-hand"
2335 " rule (which may usually be done by switching the order of the"
2336 " points or edge tensors), or in 3D that the edges form a"
2337 " right-handed coordinate system (which may also be accomplished by"
2338 " switching the order of the first two points or edge tensors)."));
2341 for (
unsigned int i = 0; i < dim; ++i)
2342 for (
unsigned int j = i + 1; j < dim; ++j)
2343 Assert((edges[i] != edges[j]),
2345 "Degenerate edges of subdivided_parallelepiped encountered."));
2348 std::vector<Point<spacedim>> points;
2353 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2354 points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2358 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2359 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2360 points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2361 edges[1] / compute_subdivisions[1] * y);
2365 for (
unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2366 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2367 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2368 points.push_back(origin +
2369 edges[0] / compute_subdivisions[0] * x +
2370 edges[1] / compute_subdivisions[1] * y +
2371 edges[2] / compute_subdivisions[2] * z);
2380 for (
unsigned int i = 0; i < dim; ++i)
2381 n_cells *= compute_subdivisions[i];
2382 std::vector<CellData<dim>> cells(
n_cells);
2388 for (
unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2390 cells[x].vertices[0] = x;
2391 cells[x].vertices[1] = x + 1;
2394 cells[x].material_id = 0;
2401 const unsigned int n_dy = compute_subdivisions[1];
2402 const unsigned int n_dx = compute_subdivisions[0];
2404 for (
unsigned int y = 0; y < n_dy; ++y)
2405 for (
unsigned int x = 0; x < n_dx; ++x)
2407 const unsigned int c = y * n_dx + x;
2408 cells[c].vertices[0] = y * (n_dx + 1) + x;
2409 cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2410 cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2411 cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2414 cells[c].material_id = 0;
2422 const unsigned int n_dz = compute_subdivisions[2];
2423 const unsigned int n_dy = compute_subdivisions[1];
2424 const unsigned int n_dx = compute_subdivisions[0];
2426 for (
unsigned int z = 0; z < n_dz; ++z)
2427 for (
unsigned int y = 0; y < n_dy; ++y)
2428 for (
unsigned int x = 0; x < n_dx; ++x)
2430 const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2432 cells[c].vertices[0] =
2433 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2434 cells[c].vertices[1] =
2435 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2436 cells[c].vertices[2] =
2437 z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2438 cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2439 (y + 1) * (n_dx + 1) + x + 1;
2440 cells[c].vertices[4] =
2441 (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2442 cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2443 y * (n_dx + 1) + x + 1;
2444 cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2445 (y + 1) * (n_dx + 1) + x;
2446 cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2447 (y + 1) * (n_dx + 1) + x + 1;
2450 cells[c].material_id = 0;
2471 for (; cell != endc; ++cell)
2475 if (cell->face(face)->at_boundary())
2476 cell->face(face)->set_boundary_id(face);
2483 template <
int dim,
int spacedim>
2486 const unsigned int repetitions,
2493 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2496 for (
unsigned int i = 0; i < dim; ++i)
2502 std::vector<unsigned int> reps(dim, repetitions);
2508 template <
int dim,
int spacedim>
2511 const std::vector<unsigned int> &repetitions,
2522 for (
unsigned int i = 0; i < dim; ++i)
2529 std::array<Point<spacedim>, dim> delta;
2530 for (
unsigned int i = 0; i < dim; ++i)
2534 delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2538 "The first dim entries of coordinates of p1 and p2 need to be different."));
2542 std::vector<Point<spacedim>> points;
2546 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2547 points.push_back(p1 + x * delta[0]);
2551 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2552 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2553 points.push_back(p1 + x * delta[0] + y * delta[1]);
2557 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
2558 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2559 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2560 points.push_back(p1 + x * delta[0] + y * delta[1] +
2569 std::vector<CellData<dim>> cells;
2574 cells.resize(repetitions[0]);
2575 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2577 cells[x].vertices[0] = x;
2578 cells[x].vertices[1] = x + 1;
2579 cells[x].material_id = 0;
2586 cells.resize(repetitions[1] * repetitions[0]);
2587 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2588 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2590 const unsigned int c = x + y * repetitions[0];
2591 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2592 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2593 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2594 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2595 cells[c].material_id = 0;
2602 const unsigned int n_x = (repetitions[0] + 1);
2603 const unsigned int n_xy =
2604 (repetitions[0] + 1) * (repetitions[1] + 1);
2606 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2607 for (
unsigned int z = 0; z < repetitions[2]; ++z)
2608 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2609 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2611 const unsigned int c = x + y * repetitions[0] +
2612 z * repetitions[0] * repetitions[1];
2613 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2614 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2615 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2616 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2617 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2618 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2619 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2620 cells[c].vertices[7] =
2621 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2622 cells[c].material_id = 0;
2644 for (
unsigned int i = 0; i < dim; ++i)
2648 "The distance between corner points must be positive."))
2652 colorize_subdivided_hyper_rectangle(
tria, p1, p2,
epsilon);
2661 const std::vector<std::vector<
double>> &step_sz,
2662 const
Point<dim> & p_1,
2663 const
Point<dim> & p_2,
2676 std::vector<std::vector<double>> step_sizes(step_sz);
2678 for (
unsigned int i = 0; i < dim; ++i)
2683 std::reverse(step_sizes[i].
begin(), step_sizes[i].
end());
2688 for (
unsigned int j = 0; j < step_sizes.at(i).size(); ++j)
2689 x += step_sizes[i][j];
2692 "The sequence of step sizes in coordinate direction " +
2694 " must be equal to the distance of the two given "
2695 "points in this coordinate direction."));
2702 std::vector<Point<dim>> points;
2708 for (
unsigned int i = 0;; ++i)
2717 if (i == step_sizes[0].size())
2720 x += step_sizes[0][i];
2728 for (
unsigned int j = 0;; ++j)
2731 for (
unsigned int i = 0;; ++i)
2733 points.push_back(
Point<dim>(p1[0] + x, p1[1] + y));
2734 if (i == step_sizes[0].size())
2737 x += step_sizes[0][i];
2740 if (j == step_sizes[1].size())
2743 y += step_sizes[1][j];
2750 for (
unsigned int k = 0;; ++k)
2753 for (
unsigned int j = 0;; ++j)
2756 for (
unsigned int i = 0;; ++i)
2759 Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2760 if (i == step_sizes[0].size())
2763 x += step_sizes[0][i];
2766 if (j == step_sizes[1].size())
2769 y += step_sizes[1][j];
2772 if (k == step_sizes[2].size())
2775 z += step_sizes[2][k];
2786 std::vector<CellData<dim>> cells;
2791 cells.resize(step_sizes[0].size());
2792 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2794 cells[x].vertices[0] = x;
2795 cells[x].vertices[1] = x + 1;
2796 cells[x].material_id = 0;
2803 cells.resize(step_sizes[1].size() * step_sizes[0].size());
2804 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2805 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2807 const unsigned int c = x + y * step_sizes[0].size();
2808 cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2809 cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2810 cells[c].vertices[2] =
2811 (y + 1) * (step_sizes[0].size() + 1) + x;
2812 cells[c].vertices[3] =
2813 (y + 1) * (step_sizes[0].size() + 1) + x + 1;
2814 cells[c].material_id = 0;
2821 const unsigned int n_x = (step_sizes[0].size() + 1);
2822 const unsigned int n_xy =
2823 (step_sizes[0].size() + 1) * (step_sizes[1].size() + 1);
2825 cells.resize(step_sizes[2].size() * step_sizes[1].size() *
2826 step_sizes[0].size());
2827 for (
unsigned int z = 0; z < step_sizes[2].size(); ++z)
2828 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2829 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2831 const unsigned int c =
2832 x + y * step_sizes[0].size() +
2833 z * step_sizes[0].size() * step_sizes[1].size();
2834 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2835 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2836 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2837 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2838 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2839 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2840 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2841 cells[c].vertices[7] =
2842 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2843 cells[c].material_id = 0;
2865 *std::min_element(step_sizes[0].
begin(), step_sizes[0].
end());
2866 for (
unsigned int i = 1; i < dim; ++i)
2868 *std::min_element(step_sizes[i].
begin(),
2869 step_sizes[i].
end()));
2870 const double epsilon = 0.01 * min_size;
2874 colorize_subdivided_hyper_rectangle(
tria, p1, p2,
epsilon);
2883 const std::vector<std::vector<double>> &spacing,
2895 for (
unsigned int i = 0; i <
n_cells; ++i)
2898 delta =
std::min(delta, spacing[0][i]);
2902 std::vector<Point<1>> points;
2904 for (
unsigned int x = 0; x <=
n_cells; ++x)
2906 points.emplace_back(ax);
2908 ax += spacing[0][x];
2911 unsigned int n_val_cells = 0;
2912 for (
unsigned int i = 0; i <
n_cells; ++i)
2916 std::vector<CellData<1>> cells(n_val_cells);
2917 unsigned int id = 0;
2918 for (
unsigned int x = 0; x <
n_cells; ++x)
2921 cells[id].vertices[0] = x;
2922 cells[id].vertices[1] = x + 1;
2941 const std::vector<std::vector<double>> &spacing,
2948 std::vector<unsigned int> repetitions(2);
2950 for (
unsigned int i = 0; i < 2; ++i)
2952 repetitions[i] = spacing[i].size();
2953 for (
unsigned int j = 0; j < repetitions[i]; ++j)
2956 delta =
std::min(delta, spacing[i][j]);
2963 std::vector<Point<2>> points;
2965 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2968 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2970 points.emplace_back(ax, ay);
2971 if (x < repetitions[0])
2972 ax += spacing[0][x];
2974 if (y < repetitions[1])
2975 ay += spacing[1][y];
2979 unsigned int n_val_cells = 0;
2980 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
2981 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
2985 std::vector<CellData<2>> cells(n_val_cells);
2986 unsigned int id = 0;
2987 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2988 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2991 cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
2992 cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2993 cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2994 cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3008 double eps = 0.01 * delta;
3010 for (; cell != endc; ++cell)
3012 Point<2> cell_center = cell->center();
3014 if (cell->face(f)->boundary_id() == 0)
3016 Point<2> face_center = cell->face(f)->center();
3017 for (
unsigned int i = 0; i < 2; ++i)
3019 if (face_center[i] < cell_center[i] -
eps)
3020 cell->face(f)->set_boundary_id(i * 2);
3021 if (face_center[i] > cell_center[i] +
eps)
3022 cell->face(f)->set_boundary_id(i * 2 + 1);
3033 const std::vector<std::vector<double>> &spacing,
3038 const unsigned int dim = 3;
3042 std::array<unsigned int, dim> repetitions;
3044 for (
unsigned int i = 0; i < dim; ++i)
3046 repetitions[i] = spacing[i].size();
3047 for (
unsigned int j = 0; j < repetitions[i]; ++j)
3050 delta =
std::min(delta, spacing[i][j]);
3057 std::vector<Point<dim>> points;
3059 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3062 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3065 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3067 points.emplace_back(ax, ay, az);
3068 if (x < repetitions[0])
3069 ax += spacing[0][x];
3071 if (y < repetitions[1])
3072 ay += spacing[1][y];
3074 if (z < repetitions[2])
3075 az += spacing[2][z];
3079 unsigned int n_val_cells = 0;
3080 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
3081 for (
unsigned int j = 0; j <
material_id.size(1); ++j)
3082 for (
unsigned int k = 0; k <
material_id.size(2); ++k)
3086 std::vector<CellData<dim>> cells(n_val_cells);
3087 unsigned int id = 0;
3088 const unsigned int n_x = (repetitions[0] + 1);
3089 const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3090 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3091 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3092 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3095 cells[id].vertices[0] = z * n_xy + y * n_x + x;
3096 cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3097 cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3098 cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3099 cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3100 cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3101 cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3102 cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3116 double eps = 0.01 * delta;
3119 for (; cell != endc; ++cell)
3123 if (cell->face(f)->boundary_id() == 0)
3125 Point<dim> face_center = cell->face(f)->center();
3126 for (
unsigned int i = 0; i < dim; ++i)
3128 if (face_center[i] < cell_center[i] -
eps)
3129 cell->face(f)->set_boundary_id(i * 2);
3130 if (face_center[i] > cell_center[i] +
eps)
3131 cell->face(f)->set_boundary_id(i * 2 + 1);
3138 template <
int dim,
int spacedim>
3141 const std::vector<unsigned int> &holes)
3150 for (
unsigned int d = 0;
d < dim; ++
d)
3157 std::array<Point<spacedim>, dim> delta;
3158 std::array<unsigned int, dim> repetitions;
3159 for (
unsigned int i = 0; i < dim; ++i)
3162 ExcMessage(
"At least one hole needed in each direction"));
3163 repetitions[i] = 2 * holes[i] + 1;
3164 delta[i][i] = (p2[i] - p1[i]);
3169 std::vector<Point<spacedim>> points;
3173 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3174 points.push_back(p1 + x * delta[0]);
3178 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3179 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3180 points.push_back(p1 + x * delta[0] + y * delta[1]);
3184 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3185 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3186 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3187 points.push_back(p1 + x * delta[0] + y * delta[1] +
3197 std::vector<CellData<dim>> cells;
3202 cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3204 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3205 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3207 if ((x % 2 == 1) && (y % 2 == 1))
3210 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
3211 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3212 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3213 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3214 cells[c].material_id = 0;
3222 const unsigned int n_x = (repetitions[0] + 1);
3223 const unsigned int n_xy =
3224 (repetitions[0] + 1) * (repetitions[1] + 1);
3226 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3229 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3230 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3231 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3234 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3235 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3236 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3237 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3238 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3239 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3240 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3241 cells[c].vertices[7] =
3242 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3243 cells[c].material_id = 0;
3271 const unsigned int ,
3283 const unsigned int ,
3295 bool inline point_in_2d_box(
const Point<2> &p,
3297 const double radius)
3299 return (std::abs(p[0] - c[0]) < radius) &&
3308 template <
int dim,
int spacedim>
3313 for (
const auto &cell :
triangulation.active_cell_iterators())
3314 for (
unsigned int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3315 length =
std::min(length, cell->line(n)->diameter());
3325 const double inner_radius,
3326 const double outer_radius,
3327 const double pad_bottom,
3328 const double pad_top,
3329 const double pad_left,
3330 const double pad_right,
3335 const unsigned int ,
3338 const bool with_padding =
3339 pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3351 for (
unsigned int n = 0; n < GeometryInfo<2>::lines_per_cell; ++n)
3352 length =
std::min(length, cell->line(n)->diameter());
3368 cell->set_manifold_id(tfi_manifold_id);
3370 const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
3371 const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
3377 auto add_sizes = [](std::vector<double> &step_sizes,
3378 const double padding,
3379 const double h) ->
void {
3382 const auto rounded =
3383 static_cast<unsigned int>(std::round(padding / h));
3386 const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3387 for (
unsigned int i = 0; i < num; ++i)
3388 step_sizes.push_back(padding / num);
3391 std::vector<std::vector<double>> step_sizes(2);
3394 add_sizes(step_sizes[0], pad_left, outer_radius);
3396 step_sizes[0].push_back(outer_radius);
3397 step_sizes[0].push_back(outer_radius);
3399 add_sizes(step_sizes[0], pad_right, outer_radius);
3402 add_sizes(step_sizes[1], pad_bottom, outer_radius);
3404 step_sizes[1].push_back(outer_radius);
3405 step_sizes[1].push_back(outer_radius);
3407 add_sizes(step_sizes[1], pad_top, outer_radius);
3412 bulk_tria, step_sizes, bl, tr,
colorize);
3415 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3417 if (internal::point_in_2d_box(cell->center(),
center, outer_radius))
3418 cells_to_remove.insert(cell);
3422 bulk_tria, cells_to_remove, tria_without_cylinder);
3424 const double tolerance =
3425 std::min(min_line_length(tria_without_cylinder),
3426 min_line_length(cylinder_tria)) /
3440 if (cell->manifold_id() == tfi_manifold_id)
3444 const auto &face = cell->face(face_n);
3445 if (face->at_boundary() &&
3446 internal::point_in_2d_box(face->center(),
3449 face->set_manifold_id(polar_manifold_id);
3451 face->set_manifold_id(tfi_manifold_id);
3462 static constexpr
double tol =
3468 const auto face = cell->face(face_n);
3469 if (face->at_boundary())
3473 if (std::abs(
center[0] - bl[0]) < tol * std::abs(bl[0]))
3474 face->set_boundary_id(0);
3476 else if (std::abs(
center[0] - tr[0]) < tol * std::abs(tr[0]))
3477 face->set_boundary_id(1);
3479 else if (std::abs(
center[1] - bl[1]) < tol * std::abs(bl[1]))
3480 face->set_boundary_id(2);
3482 else if (std::abs(
center[1] - tr[1]) < tol * std::abs(tr[1]))
3483 face->set_boundary_id(3);
3487 Assert(cell->manifold_id() == tfi_manifold_id,
3489 face->set_boundary_id(4);
3509 const double inner_radius,
3510 const double outer_radius,
3511 const double pad_bottom,
3512 const double pad_top,
3513 const double pad_left,
3514 const double pad_right,
3519 const unsigned int n_slices,
3530 Point<2>(new_center[0], new_center[1]),
3559 const double shell_region_width,
3560 const unsigned int n_shells,
3561 const double skewness,
3564 Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3565 ExcMessage(
"The width of the shell region must be less than 0.05 "
3566 "(and preferably close to 0.03)"));
3600 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3604 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3605 cells_to_remove.insert(cell);
3609 for (
const unsigned int vertex_n :
3611 if (cell->vertex(vertex_n) ==
Point<2>())
3615 cylinder_triangulation_offset =
3616 2.0 * (cell->vertex(3) -
Point<2>());
3623 bulk_tria, cells_to_remove, tria_without_cylinder);
3630 0.05 + shell_region_width,
3638 if (std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1
e-10)
3639 cell->vertex(vertex_n)[0] = -0.1;
3640 else if (std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1
e-10)
3641 cell->vertex(vertex_n)[0] = 0.1;
3647 cell->set_manifold_id(tfi_manifold_id);
3649 if (!cell->face(face_n)->at_boundary())
3650 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3652 if (0.0 < shell_region_width)
3655 ExcMessage(
"If the shell region has positive width then "
3656 "there must be at least one shell."));
3661 0.05 + shell_region_width,
3668 const double vertex_tolerance =
3669 std::min(internal::minimal_vertex_distance(shell_tria),
3670 internal::minimal_vertex_distance(cylinder_tria)) *
3676 shell_tria, cylinder_tria, temp, vertex_tolerance,
true);
3677 cylinder_tria = std::move(temp);
3683 const double vertex_tolerance =
3684 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3685 internal::minimal_vertex_distance(cylinder_tria)) /
3688 tria_without_cylinder, cylinder_tria,
tria, vertex_tolerance,
true);
3701 const double shift =
3702 std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3705 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3707 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1
e-10)
3709 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1
e-10)
3711 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1
e-10)
3718 if (cell->manifold_id() == polar_manifold_id)
3719 cell->set_all_manifold_ids(polar_manifold_id);
3724 if (cell->manifold_id() != polar_manifold_id &&
3725 cell->manifold_id() != tfi_manifold_id)
3730 std::vector<Point<2> *> cylinder_pointers;
3732 if (face->manifold_id() == polar_manifold_id)
3734 cylinder_pointers.push_back(&face->vertex(0));
3735 cylinder_pointers.push_back(&face->vertex(1));
3738 std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3739 cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3740 cylinder_pointers.end()),
3741 cylinder_pointers.end());
3745 for (
const Point<2> *
const ptr : cylinder_pointers)
3749 for (
Point<2> *
const ptr : cylinder_pointers)
3761 if (face->at_boundary())
3765 if (std::abs(
center[0] - 0.0) < 1
e-10)
3766 face->set_boundary_id(0);
3768 else if (std::abs(
center[0] - 2.2) < 1
e-10)
3769 face->set_boundary_id(1);
3771 else if (face->manifold_id() == polar_manifold_id)
3772 face->set_boundary_id(2);
3777 std::abs(
center[1] - 0.41) < 1.0e-10,
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);
3820 if (face->boundary_id() == 4 || face->boundary_id() == 5)
3821 face->set_boundary_id(3);
3826 template <
int dim,
int spacedim>
3829 const std::vector<unsigned int> &sizes,
3839 for (
unsigned int d = 0;
d < dim; ++
d)
3842 std::vector<Point<spacedim>> points;
3847 std::vector<CellData<dim>> cells(
n_cells);
3852 for (
unsigned int d = 0;
d < dim; ++
d)
3853 p(
d) = 0.5 * dimensions[
d] *
3856 points.push_back(p);
3857 cells[0].vertices[i] = i;
3859 cells[0].material_id = 0;
3870 for (
unsigned int j = 0; j < sizes[face]; ++j, ++
cell_index)
3872 const unsigned int last_cell = (j == 0) ? 0
U : (
cell_index - 1);
3874 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
3877 const unsigned int cellv =
3879 const unsigned int ocellv =
3883 cells[last_cell].vertices[cellv];
3886 cells[
cell_index].vertices[cellv] = points.size();
3891 points.push_back(p);
4048 const double thickness,
4052 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4054 std::vector<Point<2>>
vertices(16);
4056 coords[0] = left - thickness;
4059 coords[3] = right + thickness;
4062 for (
const double y : coords)
4063 for (
const double x : coords)
4068 std::vector<CellData<2>> cells(9);
4070 for (
unsigned int i0 = 0; i0 < 3; ++i0)
4071 for (
unsigned int i1 = 0; i1 < 3; ++i1)
4073 cells[k].vertices[0] = i1 + 4 * i0;
4074 cells[k].vertices[1] = i1 + 4 * i0 + 1;
4075 cells[k].vertices[2] = i1 + 4 * i0 + 4;
4076 cells[k].vertices[3] = i1 + 4 * i0 + 5;
4078 cells[k].material_id = materials[k];
4096 const double rl2 = (right + left) / 2;
4107 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4112 for (
unsigned int i = 0; i < 4; ++i)
4114 for (
unsigned int j = 0; j < 4; ++j)
4115 cells[i].
vertices[j] = cell_vertices[i][j];
4116 cells[i].material_id = 0;
4126 cell->face(1)->set_boundary_id(1);
4128 cell->face(0)->set_boundary_id(2);
4137 const double radius_0,
4138 const double radius_1,
4139 const double half_length)
4143 vertices_tmp[0] =
Point<2>(-half_length, -radius_0);
4144 vertices_tmp[1] =
Point<2>(half_length, -radius_1);
4145 vertices_tmp[2] =
Point<2>(-half_length, radius_0);
4146 vertices_tmp[3] =
Point<2>(half_length, radius_1);
4153 cell_vertices[0][i] = i;
4158 cells[0].vertices[i] = cell_vertices[0][i];
4160 cells[0].material_id = 0;
4165 cell->face(0)->set_boundary_id(1);
4166 cell->face(1)->set_boundary_id(2);
4168 for (
unsigned int i = 2; i < 4; ++i)
4169 cell->face(i)->set_boundary_id(0);
4190 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4194 for (
unsigned int i = 0; i < 3; ++i)
4196 for (
unsigned int j = 0; j < 4; ++j)
4197 cells[i].
vertices[j] = cell_vertices[i][j];
4198 cells[i].material_id = 0;
4210 cell->face(0)->set_boundary_id(0);
4211 cell->face(2)->set_boundary_id(1);
4214 cell->face(1)->set_boundary_id(2);
4215 cell->face(2)->set_boundary_id(1);
4216 cell->face(3)->set_boundary_id(3);
4219 cell->face(0)->set_boundary_id(0);
4220 cell->face(1)->set_boundary_id(4);
4221 cell->face(3)->set_boundary_id(5);
4227 template <
int dim,
int spacedim>
4230 const std::vector<unsigned int> &repetitions,
4233 const std::vector<int> & n_cells_to_remove)
4239 for (
unsigned int d = 0;
d < dim; ++
d)
4242 ExcMessage(
"Attempting to cut away too many cells."));
4252 std::array<double, dim> h;
4254 for (
unsigned int d = 0;
d < dim; ++
d)
4257 h[
d] = (top_right[
d] - bottom_left[
d]) / repetitions[
d];
4259 if (n_cells_to_remove[
d] >= 0)
4263 h[
d] *
std::fabs(n_cells_to_remove[
d]) + bottom_left[
d];
4268 cut_step[
d] = top_right[
d] - h[
d] *
std::fabs(n_cells_to_remove[
d]);
4274 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4278 bool remove_cell =
true;
4279 for (
unsigned int d = 0;
d < dim && remove_cell; ++
d)
4280 if ((n_cells_to_remove[
d] > 0 && cell->center()[
d] >= cut_step[
d]) ||
4281 (n_cells_to_remove[
d] < 0 && cell->center()[
d] <= cut_step[
d]))
4282 remove_cell =
false;
4284 cells_to_remove.insert(cell);
4299 const double radius,
4300 const bool internal_manifolds)
4305 const double a = 1. / (1 +
std::sqrt(2.0));
4316 const int cell_vertices[5][4] = {
4317 {0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}};
4321 for (
unsigned int i = 0; i < 5; ++i)
4323 for (
unsigned int j = 0; j < 4; ++j)
4324 cells[i].
vertices[j] = cell_vertices[i][j];
4325 cells[i].material_id = 0;
4335 if (internal_manifolds)
4345 const double inner_radius,
4346 const double outer_radius,
4350 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4364 const unsigned int N =
4365 (
n_cells == 0 ?
static_cast<unsigned int>(
4366 std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
4367 (outer_radius - inner_radius))) :
4377 for (
unsigned int i = 0; i <
N; ++i)
4380 Point<2>(std::cos(2 * pi * i /
N), std::sin(2 * pi * i /
N)) *
4390 for (
unsigned int i = 0; i <
N; ++i)
4392 cells[i].vertices[0] = i;
4393 cells[i].vertices[1] = (i + 1) %
N;
4394 cells[i].vertices[2] =
N + i;
4395 cells[i].vertices[3] =
N + ((i + 1) %
N);
4397 cells[i].material_id = 0;
4403 colorize_hyper_shell(
tria,
center, inner_radius, outer_radius);
4416 const double inner_radius,
4417 const double outer_radius,
4421 tria, outer_center, inner_radius, outer_radius,
n_cells,
true);
4425 outer_radius - inner_radius > outer_center.
distance(inner_center),
4427 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4431 std::set<Point<dim> *> vertices_to_move;
4434 if (face->boundary_id() == 0)
4435 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4436 vertices_to_move.insert(&face->vertex(v));
4438 const auto shift = inner_center - outer_center;
4439 for (
const auto &p : vertices_to_move)
4467 const double radius,
4468 const double half_length)
4470 Point<2> p1(-half_length, -radius);
4479 switch (f->boundary_id())
4482 f->set_boundary_id(1);
4485 f->set_boundary_id(2);
4488 f->set_boundary_id(0);
4525 const double radius)
4527 const unsigned int dim = 2;
4541 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4545 for (
unsigned int i = 0; i < 3; ++i)
4547 for (
unsigned int j = 0; j < 4; ++j)
4548 cells[i].
vertices[j] = cell_vertices[i][j];
4549 cells[i].material_id = 0;
4566 if (cell->face(i)->boundary_id() ==
4572 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius ||
4573 cell->face(i)->center()(1) < p(1) + 1.e-5 * radius)
4575 cell->face(i)->set_boundary_id(1);
4589 const double radius)
4594 const double a = 1. / (1 +
std::sqrt(2.0));
4605 const int cell_vertices[5][4] = {{0, 1, 2, 3},
4612 for (
unsigned int i = 0; i < 4; ++i)
4614 for (
unsigned int j = 0; j < 4; ++j)
4615 cells[i].
vertices[j] = cell_vertices[i][j];
4616 cells[i].material_id = 0;
4633 if (cell->face(i)->boundary_id() ==
4638 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius)
4640 cell->face(i)->set_boundary_id(1);
4656 const double inner_radius,
4657 const double outer_radius,
4661 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4674 const unsigned int N =
4675 (
n_cells == 0 ?
static_cast<unsigned int>(
4676 std::ceil((pi * (outer_radius + inner_radius) / 2) /
4677 (outer_radius - inner_radius))) :
4686 std::vector<Point<2>>
vertices(2 * (
N + 1));
4687 for (
unsigned int i = 0; i <=
N; ++i)
4696 Point<2>(((i == 0) || (i ==
N) ? 0 : std::cos(pi * i /
N - pi / 2)),
4697 std::sin(pi * i /
N - pi / 2)) *
4708 for (
unsigned int i = 0; i <
N; ++i)
4710 cells[i].vertices[0] = i;
4711 cells[i].vertices[1] = (i + 1) % (
N + 1);
4712 cells[i].vertices[2] =
N + 1 + i;
4713 cells[i].vertices[3] =
N + 1 + ((i + 1) % (
N + 1));
4715 cells[i].material_id = 0;
4723 for (; cell !=
tria.
end(); ++cell)
4725 cell->face(2)->set_boundary_id(1);
4727 tria.
begin()->face(0)->set_boundary_id(3);
4729 tria.
last()->face(1)->set_boundary_id(2);
4740 const double inner_radius,
4741 const double outer_radius,
4745 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4758 const unsigned int N =
4759 (
n_cells == 0 ?
static_cast<unsigned int>(
4760 std::ceil((pi * (outer_radius + inner_radius) / 4) /
4761 (outer_radius - inner_radius))) :
4770 std::vector<Point<2>>
vertices(2 * (
N + 1));
4771 for (
unsigned int i = 0; i <=
N; ++i)
4779 std::sin(pi * i /
N / 2)) *
4790 for (
unsigned int i = 0; i <
N; ++i)
4792 cells[i].vertices[0] = i;
4793 cells[i].vertices[1] = (i + 1) % (
N + 1);
4794 cells[i].vertices[2] =
N + 1 + i;
4795 cells[i].vertices[3] =
N + 1 + ((i + 1) % (
N + 1));
4797 cells[i].material_id = 0;
4805 for (; cell !=
tria.
end(); ++cell)
4807 cell->face(2)->set_boundary_id(1);
4809 tria.
begin()->face(0)->set_boundary_id(3);
4811 tria.
last()->face(1)->set_boundary_id(2);
4828 const double rl2 = (right + left) / 2;
4829 const double len = (right - left) / 2.;
4842 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4843 {9, 4, 2, 5, 19, 14, 12, 15},
4844 {3, 2, 7, 6, 13, 12, 17, 16},
4845 {2, 5, 6, 8, 12, 15, 16, 18}};
4847 for (
unsigned int i = 0; i < 4; ++i)
4849 for (
unsigned int j = 0; j < 8; ++j)
4850 cells[i].
vertices[j] = cell_vertices[i][j];
4851 cells[i].material_id = 0;
4861 cell->face(1)->set_boundary_id(1);
4863 cell->face(0)->set_boundary_id(2);
4875 const double thickness,
4879 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4881 std::vector<Point<3>>
vertices(64);
4883 coords[0] = left - thickness;
4886 coords[3] = right + thickness;
4889 for (
const double z : coords)
4890 for (
const double y : coords)
4891 for (
const double x : coords)
4895 24, 26, 5, 4, 6, 1, 0,
4896 2, 9, 8, 10, 37, 36, 38,
4897 33, 32, 34, 41, 40, 42};
4899 std::vector<CellData<3>> cells(27);
4901 for (
unsigned int z = 0; z < 3; ++z)
4902 for (
unsigned int y = 0; y < 3; ++y)
4903 for (
unsigned int x = 0; x < 3; ++x)
4905 cells[k].vertices[0] = x + 4 * y + 16 * z;
4906 cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
4907 cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
4908 cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
4909 cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
4910 cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
4911 cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
4912 cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
4914 cells[k].material_id = materials[k];
4927 const double radius_0,
4928 const double radius_1,
4929 const double half_length)
4932 ExcMessage(
"The output triangulation object needs to be empty."));
4937 const auto n_slices = 1 +
static_cast<unsigned int>(
std::ceil(
4938 half_length /
std::max(radius_0, radius_1)));
4951 auto shift_radii = [=](
const Point<3> &p) {
4952 const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
4953 const double factor = slope * (p[0] - -half_length) + 1.0;
4954 return Point<3>(p[0], factor * p[1], factor * p[2]);
4960 for (
const auto &face :
triangulation.active_face_iterators())
4961 if (face->at_boundary())
4963 if (std::abs(face->center()[0] - -half_length) < 1
e-8 * half_length)
4964 face->set_boundary_id(1);
4965 else if (std::abs(face->center()[0] - half_length) <
4967 face->set_boundary_id(2);
5017 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5018 {1, 2, 10, 11, 4, 5, 13, 14},
5019 {3, 4, 12, 13, 6, 7, 15, 16},
5020 {4, 5, 13, 14, 7, 8, 16, 17},
5021 {9, 10, 18, 19, 12, 13, 21, 22},
5022 {10, 11, 19, 20, 13, 14, 22, 23},
5023 {12, 13, 21, 22, 15, 16, 24, 25}};
5027 for (
unsigned int i = 0; i < 7; ++i)
5029 for (
unsigned int j = 0; j < 8; ++j)
5030 cells[i].
vertices[j] = cell_vertices[i][j];
5031 cells[i].material_id = 0;
5052 const double radius,
5053 const bool internal_manifold)
5059 const unsigned int n_vertices = 16;
5085 const unsigned int n_cells = 7;
5086 const int cell_vertices[
n_cells][8] = {
5087 {0, 1, 4, 5, 3, 2, 7, 6},
5088 {8, 9, 12, 13, 0, 1, 4, 5},
5089 {9, 13, 1, 5, 10, 14, 2, 6},
5090 {11, 10, 3, 2, 15, 14, 7, 6},
5091 {8, 0, 12, 4, 11, 3, 15, 7},
5092 {8, 9, 0, 1, 11, 10, 3, 2},
5093 {12, 4, 13, 5, 15, 7, 14, 6}};
5097 for (
unsigned int i = 0; i <
n_cells; ++i)
5100 cells[i].vertices[j] = cell_vertices[i][j];
5101 cells[i].material_id = 0;
5111 if (internal_manifold)
5119 const unsigned int n_rotate_middle_square)
5122 ExcMessage(
"The number of rotation by pi/2 of the right square "
5123 "must be in the half-open range [0,4)."))
5125 constexpr
unsigned int dim = 2;
5127 const
unsigned int n_cells = 5;
5146 unsigned int cell_vertices[
n_cells][4] = {{0, 1, 2, 3},
5152 switch (n_rotate_middle_square)
5156 cell_vertices[1][0] = 4;
5157 cell_vertices[1][1] = 5;
5158 cell_vertices[1][2] = 1;
5159 cell_vertices[1][3] = 3;
5165 cell_vertices[1][0] = 5;
5166 cell_vertices[1][1] = 3;
5167 cell_vertices[1][2] = 4;
5168 cell_vertices[1][3] = 1;
5174 cell_vertices[1][0] = 3;
5175 cell_vertices[1][1] = 1;
5176 cell_vertices[1][2] = 5;
5177 cell_vertices[1][3] = 4;
5189 for (
const unsigned int vertex_index :
5204 const bool face_orientation,
5205 const bool face_flip,
5206 const bool face_rotation,
5207 const bool manipulate_left_cube)
5209 constexpr
unsigned int dim = 3;
5211 const unsigned int n_cells = 2;
5212 std::vector<CellData<dim>> cells(
n_cells);
5228 unsigned int cell_vertices[
n_cells][8] = {
5229 {0, 1, 2, 3, 4, 5, 6, 7},
5230 {1, 8, 3, 9, 5, 10, 7, 11}};
5233 const unsigned int this_case = 4 *
static_cast<int>(face_orientation) +
5234 2 *
static_cast<int>(face_flip) +
5235 static_cast<int>(face_rotation);
5237 if (manipulate_left_cube)
5243 cell_vertices[0][0] = 1;
5244 cell_vertices[0][1] = 0;
5245 cell_vertices[0][2] = 5;
5246 cell_vertices[0][3] = 4;
5247 cell_vertices[0][4] = 3;
5248 cell_vertices[0][5] = 2;
5249 cell_vertices[0][6] = 7;
5250 cell_vertices[0][7] = 6;
5256 cell_vertices[0][0] = 5;
5257 cell_vertices[0][1] = 4;
5258 cell_vertices[0][2] = 7;
5259 cell_vertices[0][3] = 6;
5260 cell_vertices[0][4] = 1;
5261 cell_vertices[0][5] = 0;
5262 cell_vertices[0][6] = 3;
5263 cell_vertices[0][7] = 2;
5269 cell_vertices[0][0] = 7;
5270 cell_vertices[0][1] = 6;
5271 cell_vertices[0][2] = 3;
5272 cell_vertices[0][3] = 2;
5273 cell_vertices[0][4] = 5;
5274 cell_vertices[0][5] = 4;
5275 cell_vertices[0][6] = 1;
5276 cell_vertices[0][7] = 0;
5281 cell_vertices[0][0] = 3;
5282 cell_vertices[0][1] = 2;
5283 cell_vertices[0][2] = 1;
5284 cell_vertices[0][3] = 0;
5285 cell_vertices[0][4] = 7;
5286 cell_vertices[0][5] = 6;
5287 cell_vertices[0][6] = 5;
5288 cell_vertices[0][7] = 4;
5294 cell_vertices[0][0] = 0;
5295 cell_vertices[0][1] = 1;
5296 cell_vertices[0][2] = 2;
5297 cell_vertices[0][3] = 3;
5298 cell_vertices[0][4] = 4;
5299 cell_vertices[0][5] = 5;
5300 cell_vertices[0][6] = 6;
5301 cell_vertices[0][7] = 7;
5307 cell_vertices[0][0] = 2;
5308 cell_vertices[0][1] = 3;
5309 cell_vertices[0][2] = 6;
5310 cell_vertices[0][3] = 7;
5311 cell_vertices[0][4] = 0;
5312 cell_vertices[0][5] = 1;
5313 cell_vertices[0][6] = 4;
5314 cell_vertices[0][7] = 5;
5320 cell_vertices[0][0] = 6;
5321 cell_vertices[0][1] = 7;
5322 cell_vertices[0][2] = 4;
5323 cell_vertices[0][3] = 5;
5324 cell_vertices[0][4] = 2;
5325 cell_vertices[0][5] = 3;
5326 cell_vertices[0][6] = 0;
5327 cell_vertices[0][7] = 1;
5333 cell_vertices[0][0] = 4;
5334 cell_vertices[0][1] = 5;
5335 cell_vertices[0][2] = 0;
5336 cell_vertices[0][3] = 1;
5337 cell_vertices[0][4] = 6;
5338 cell_vertices[0][5] = 7;
5339 cell_vertices[0][6] = 2;
5340 cell_vertices[0][7] = 3;
5351 cell_vertices[1][0] = 8;
5352 cell_vertices[1][1] = 1;
5353 cell_vertices[1][2] = 10;
5354 cell_vertices[1][3] = 5;
5355 cell_vertices[1][4] = 9;
5356 cell_vertices[1][5] = 3;
5357 cell_vertices[1][6] = 11;
5358 cell_vertices[1][7] = 7;
5364 cell_vertices[1][0] = 10;
5365 cell_vertices[1][1] = 5;
5366 cell_vertices[1][2] = 11;
5367 cell_vertices[1][3] = 7;
5368 cell_vertices[1][4] = 8;
5369 cell_vertices[1][5] = 1;
5370 cell_vertices[1][6] = 9;
5371 cell_vertices[1][7] = 3;
5377 cell_vertices[1][0] = 11;
5378 cell_vertices[1][1] = 7;
5379 cell_vertices[1][2] = 9;
5380 cell_vertices[1][3] = 3;
5381 cell_vertices[1][4] = 10;
5382 cell_vertices[1][5] = 5;
5383 cell_vertices[1][6] = 8;
5384 cell_vertices[1][7] = 1;
5390 cell_vertices[1][0] = 9;
5391 cell_vertices[1][1] = 3;
5392 cell_vertices[1][2] = 8;
5393 cell_vertices[1][3] = 1;
5394 cell_vertices[1][4] = 11;
5395 cell_vertices[1][5] = 7;
5396 cell_vertices[1][6] = 10;
5397 cell_vertices[1][7] = 5;
5403 cell_vertices[1][0] = 1;
5404 cell_vertices[1][1] = 8;
5405 cell_vertices[1][2] = 3;
5406 cell_vertices[1][3] = 9;
5407 cell_vertices[1][4] = 5;
5408 cell_vertices[1][5] = 10;
5409 cell_vertices[1][6] = 7;
5410 cell_vertices[1][7] = 11;
5416 cell_vertices[1][0] = 5;
5417 cell_vertices[1][1] = 10;
5418 cell_vertices[1][2] = 1;
5419 cell_vertices[1][3] = 8;
5420 cell_vertices[1][4] = 7;
5421 cell_vertices[1][5] = 11;
5422 cell_vertices[1][6] = 3;
5423 cell_vertices[1][7] = 9;
5429 cell_vertices[1][0] = 7;
5430 cell_vertices[1][1] = 11;
5431 cell_vertices[1][2] = 5;
5432 cell_vertices[1][3] = 10;
5433 cell_vertices[1][4] = 3;
5434 cell_vertices[1][5] = 9;
5435 cell_vertices[1][6] = 1;
5436 cell_vertices[1][7] = 8;
5442 cell_vertices[1][0] = 3;
5443 cell_vertices[1][1] = 9;
5444 cell_vertices[1][2] = 7;
5445 cell_vertices[1][3] = 11;
5446 cell_vertices[1][4] = 1;
5447 cell_vertices[1][5] = 8;
5448 cell_vertices[1][6] = 5;
5449 cell_vertices[1][7] = 10;
5459 for (
const unsigned int vertex_index :
5473 template <
int spacedim>
5477 const double radius)
5481 const std::set<types::boundary_id> boundary_ids = {0};
5493 const unsigned int x_subdivisions,
5494 const double radius,
5495 const double half_length)
5503 const double initial_height = -half_length;
5504 const double height_increment = 2. * half_length / x_subdivisions;
5506 for (
unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep)
5508 const double height = initial_height + height_increment * rep;
5512 vertices.emplace_back(-a, height, -a);
5513 vertices.emplace_back(a, height, -a);
5514 vertices.emplace_back(-a, height, a);
5515 vertices.emplace_back(a, height, a);
5523 const double h = vertex(1);
5524 vertex(1) = -vertex(0);
5528 std::vector<std::vector<int>> cell_vertices;
5529 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5530 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5531 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5532 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5533 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5535 for (
unsigned int rep = 1; rep < x_subdivisions; ++rep)
5537 for (
unsigned int i = 0; i < 5; ++i)
5539 std::vector<int> new_cell_vertices(8);
5540 for (
unsigned int j = 0; j < 8; ++j)
5541 new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep;
5542 cell_vertices.push_back(new_cell_vertices);
5546 unsigned int n_cells = x_subdivisions * 5;
5550 for (
unsigned int i = 0; i <
n_cells; ++i)
5552 for (
unsigned int j = 0; j < 8; ++j)
5553 cells[i].
vertices[j] = cell_vertices[i][j];
5554 cells[i].material_id = 0;
5576 const double tolerance = 1
e-5 *
std::min(radius, half_length);
5580 if (cell->at_boundary(i))
5582 if (cell->face(i)->center()(0) > half_length - tolerance)
5584 cell->face(i)->set_boundary_id(2);
5587 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5589 if ((
std::fabs(cell->face(i)->line(
e)->vertex(0)[1]) == a) ||
5590 (
std::fabs(cell->face(i)->line(
e)->vertex(0)[2]) == a) ||
5591 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[1]) == a) ||
5592 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[2]) == a))
5594 cell->face(i)->line(
e)->set_boundary_id(2);
5595 cell->face(i)->line(
e)->set_manifold_id(
5599 else if (cell->face(i)->center()(0) < -half_length + tolerance)
5601 cell->face(i)->set_boundary_id(1);
5604 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5606 if ((
std::fabs(cell->face(i)->line(
e)->vertex(0)[1]) == a) ||
5607 (
std::fabs(cell->face(i)->line(
e)->vertex(0)[2]) == a) ||
5608 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[1]) == a) ||
5609 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[2]) == a))
5611 cell->face(i)->line(
e)->set_boundary_id(1);
5612 cell->face(i)->line(
e)->set_manifold_id(
5624 const double radius,
5625 const double half_length)
5634 const double radius)
5636 const unsigned int dim = 3;
5644 const double a = 0.528;
5645 const double b = 0.4533;
5646 const double c = 0.3752;
5663 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5664 {1, 6, 2, 4, 8, 13, 9, 11},
5665 {5, 3, 6, 4, 12, 10, 13, 11},
5666 {7, 9, 10, 11, 14, 8, 12, 13}};
5670 for (
unsigned int i = 0; i < 4; ++i)
5672 for (
unsigned int j = 0; j < 8; ++j)
5673 cells[i].
vertices[j] = cell_vertices[i][j];
5674 cells[i].material_id = 0;
5690 if (cell->face(i)->boundary_id() ==
5695 if (cell->face(i)->center()(0) <
center(0) + 1.e-5 * radius ||
5696 cell->face(i)->center()(1) <
center(1) + 1.e-5 * radius ||
5697 cell->face(i)->center()(2) <
center(2) + 1.e-5 * radius)
5699 cell->face(i)->set_boundary_id(1);
5703 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5706 const Point<3> line_vertices[2] = {
5707 cell->face(i)->line(j)->vertex(0),
5708 cell->face(i)->line(j)->vertex(1)};
5714 cell->face(i)->line(j)->set_boundary_id(1);
5715 cell->face(i)->line(j)->set_manifold_id(
5733 const double radius)
5739 const double b = a / 2.0;
5740 const double c =
d / 2.0;
5742 const double hb = radius *
std::sqrt(3.0) / 4.0;
5743 const double hc = radius *
std::sqrt(3.0) / 2.0;
5765 int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5766 {0, 2, 8, 10, 6, 4, 14, 12},
5767 {2, 3, 10, 11, 4, 5, 12, 13},
5768 {1, 7, 9, 15, 3, 5, 11, 13},
5769 {6, 4, 14, 12, 7, 5, 15, 13},
5770 {8, 10, 9, 11, 14, 12, 15, 13}};
5774 for (
unsigned int i = 0; i < 6; ++i)
5776 for (
unsigned int j = 0; j < 8; ++j)
5777 cells[i].
vertices[j] = cell_vertices[i][j];
5778 cells[i].material_id = 0;
5800 if (!cell->at_boundary(i))
5806 if (cell->face(i)->center()(0) <
center(0) + 1.e-5 * radius)
5808 cell->face(i)->set_boundary_id(1);
5810 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5813 const Point<3> line_vertices[2] = {
5814 cell->face(i)->line(j)->vertex(0),
5815 cell->face(i)->line(j)->vertex(1)};
5821 cell->face(i)->line(j)->set_boundary_id(1);
5822 cell->face(i)->line(j)->set_manifold_id(
5839 const double radius)
5852 for (
unsigned int round = 0; round < dim; ++round)
5857 std::vector<Point<dim>> new_points(tria_copy.
n_vertices());
5859 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5867 else if (round == 1)
5869 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5878 else if (round == 2)
5879 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5892 std::vector<CellData<dim>> cells;
5893 cells.reserve(tria_copy.
n_cells());
5898 data.
vertices[v] = cell->vertex_index(v);
5901 cells.push_back(data);
5909 if (round == dim - 1)
5919 if (cell->center().norm_square() > 0.4 * radius)
5920 cell->set_manifold_id(1);
5939 const double inner_radius,
5940 const double outer_radius)
5943 std::vector<CellData<3>> cells;
5945 const double irad = inner_radius /
std::sqrt(3.0);
5946 const double orad = outer_radius /
std::sqrt(3.0);
5949 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
5959 for (
unsigned int i = 0; i < 8; ++i)
5960 vertices.push_back(p + hexahedron[i] * irad);
5961 for (
unsigned int i = 0; i < 8; ++i)
5962 vertices.push_back(p + hexahedron[i] * orad);
5964 const unsigned int n_cells = 6;
5965 const int cell_vertices[
n_cells][8] = {
5966 {8, 9, 10, 11, 0, 1, 2, 3},
5967 {9, 11, 1, 3, 13, 15, 5, 7},
5968 {12, 13, 4, 5, 14, 15, 6, 7},
5969 {8, 0, 10, 2, 12, 4, 14, 6},
5970 {8, 9, 0, 1, 12, 13, 4, 5},
5971 {10, 2, 11, 3, 14, 6, 15, 7}};
5975 for (
unsigned int i = 0; i <
n_cells; ++i)
5978 cells[i].vertices[j] = cell_vertices[i][j];
5979 cells[i].material_id = 0;
5990 const double inner_radius,
5991 const double outer_radius)
5994 std::vector<CellData<3>> cells;
5996 const double irad = inner_radius /
std::sqrt(3.0);
5997 const double orad = outer_radius /
std::sqrt(3.0);
6003 static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0},
6011 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
6020 for (
unsigned int i = 0; i < 8; ++i)
6021 vertices.push_back(p + hexahedron[i] * irad);
6022 for (
unsigned int i = 0; i < 6; ++i)
6023 vertices.push_back(p + octahedron[i] * inner_radius);
6024 for (
unsigned int i = 0; i < 8; ++i)
6025 vertices.push_back(p + hexahedron[i] * orad);
6026 for (
unsigned int i = 0; i < 6; ++i)
6027 vertices.push_back(p + octahedron[i] * outer_radius);
6029 const unsigned int n_cells = 12;
6030 const unsigned int rhombi[
n_cells][4] = {{10, 4, 0, 8},
6045 for (
unsigned int i = 0; i <
n_cells; ++i)
6047 for (
unsigned int j = 0; j < 4; ++j)
6049 cells[i].vertices[j] = rhombi[i][j];
6050 cells[i].vertices[j + 4] = rhombi[i][j] + 14;
6052 cells[i].material_id = 0;
6062 const unsigned int n,
6063 const unsigned int n_refinement_steps,
6065 const double inner_radius,
6066 const double outer_radius)