48 : airfoil_type(
"NACA")
50 , joukowski_center(-0.1, 0.14)
56 , incline_factor(0.35)
59 , n_subdivision_x_0(3)
60 , n_subdivision_x_1(2)
61 , n_subdivision_x_2(5)
63 , airfoil_sampling_factor(2)
66 airfoil_length <= height,
68 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
70 Assert(incline_factor < 1.0 && incline_factor >= 0.0,
71 ExcMessage(
"incline_factor has to be in [0,1)!"));
84 "Mesh height measured from airfoil nose to horizontal boundaries");
88 "Length measured from airfoil leading edge to vertical outlet boundary");
92 "Define obliqueness of the vertical mesh around the airfoil");
101 "Type of airfoil geometry, either NACA or Joukowski airfoil",
116 "Joukowski circle center coordinates");
119 "Joukowski airfoil length leading to trailing edge");
127 "Number of global refinements");
129 "NumberSubdivisionX0",
131 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
133 "NumberSubdivisionX1",
135 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
137 "NumberSubdivisionX2",
139 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
142 "Number of subdivisions normal to airfoil");
146 "Factor to obtain a finer mesh at the airfoil surface");
161 static const unsigned int id_block_1 = 1;
162 static const unsigned int id_block_2 = 2;
163 static const unsigned int id_block_3 = 3;
164 static const unsigned int id_block_4 = 4;
165 static const unsigned int id_block_5 = 5;
166 static const unsigned int id_block_6 = 6;
171 MeshGenerator(
const AdditionalData &data)
172 : refinements(data.refinements)
173 , n_subdivision_x_0(data.n_subdivision_x_0)
174 , n_subdivision_x_1(data.n_subdivision_x_1)
175 , n_subdivision_x_2(data.n_subdivision_x_2)
176 , n_subdivision_y(data.n_subdivision_y)
177 , height(data.height)
178 , length_b2(data.length_b2)
179 , incline_factor(data.incline_factor)
180 , bias_factor(data.bias_factor)
182 , n_cells_x_0(
Utilities::
pow(2, refinements) * n_subdivision_x_0)
183 , n_cells_x_1(
Utilities::
pow(2, refinements) * n_subdivision_x_1)
184 , n_cells_x_2(
Utilities::
pow(2, refinements) * n_subdivision_x_2)
185 , n_cells_y(
Utilities::
pow(2, refinements) * n_subdivision_y)
186 , n_points_on_each_side(n_cells_x_0 + n_cells_x_1 + 1)
188 , airfoil_1D(set_airfoil_length(
190 data.airfoil_type ==
"Joukowski" ?
191 joukowski(data.joukowski_center,
192 n_points_on_each_side,
193 data.airfoil_sampling_factor) :
194 (data.airfoil_type ==
"NACA" ?
196 n_points_on_each_side,
197 data.airfoil_sampling_factor) :
200 std::vector<Point<2>>{
204 data.airfoil_length))
205 , end_b0_x_u(airfoil_1D[0][n_cells_x_0](0))
206 , end_b0_x_l(airfoil_1D[1][n_cells_x_0](0))
207 , nose_x(airfoil_1D[0].front()(0))
208 , tail_x(airfoil_1D[0].back()(0))
209 , tail_y(airfoil_1D[0].back()(1))
210 , center_mesh(0.5 *
std::abs(end_b0_x_u + end_b0_x_l))
211 , length_b1_x(tail_x - center_mesh)
213 (edge_length +
std::abs(nose_x - center_mesh))))
217 ,
A(nose_x - edge_length, 0)
220 , D(center_mesh, height)
222 ,
F(center_mesh, -height)
226 , J(tail_x + length_b2, 0)
230 Assert(data.airfoil_type ==
"Joukowski" ||
231 data.airfoil_type ==
"NACA",
243 make_coarse_grid(tria_grid);
245 set_boundary_ids(tria_grid);
247 if (periodic_faces !=
nullptr)
250 tria_grid, 5, 4, 1, *periodic_faces);
267 (void)periodic_faces;
274 const unsigned int refinements;
277 const unsigned int n_subdivision_x_0;
280 const unsigned int n_subdivision_x_1;
283 const unsigned int n_subdivision_x_2;
287 const unsigned int n_subdivision_y;
294 const double length_b2;
298 const double incline_factor;
301 const double bias_factor;
304 const double edge_length;
307 const unsigned int n_cells_x_0;
310 const unsigned int n_cells_x_1;
313 const unsigned int n_cells_x_2;
317 const unsigned int n_cells_y;
320 const unsigned int n_points_on_each_side;
324 const std::array<std::vector<Point<2>>, 2> airfoil_1D;
328 const double end_b0_x_u;
332 const double end_b0_x_l;
346 const double center_mesh;
349 const double length_b1_x;
376 const Point<2> A, B,
C, D,
E,
F, G, H, I, J, K,
L;
415 static std::array<std::vector<Point<2>>, 2>
416 joukowski(
const Point<2> & centerpoint,
417 const unsigned int number_points,
418 const unsigned int factor)
420 std::array<std::vector<Point<2>>, 2> airfoil_1D;
421 const unsigned int total_points = 2 * number_points - 2;
422 const unsigned int n_airfoilpoints = factor * total_points;
424 const auto jouk_points =
425 joukowski_transform(joukowski_circle(centerpoint, n_airfoilpoints));
428 std::vector<Point<2>> upper_points;
429 std::vector<Point<2>> lower_points;
433 unsigned int nose_index = 0;
434 unsigned int tail_index = 0;
435 double nose_x_coordinate = 0;
436 double tail_x_coordinate = 0;
440 for (
unsigned int i = 0; i < jouk_points.size(); i++)
442 if (jouk_points[i](0) < nose_x_coordinate)
444 nose_x_coordinate = jouk_points[i](0);
447 if (jouk_points[i](0) > tail_x_coordinate)
449 tail_x_coordinate = jouk_points[i](0);
455 for (
unsigned int i = tail_index; i < jouk_points.size(); i++)
456 upper_points.emplace_back(jouk_points[i]);
457 for (
unsigned int i = 0; i <= nose_index; i++)
458 upper_points.emplace_back(jouk_points[i]);
459 std::reverse(upper_points.begin(), upper_points.end());
462 lower_points.insert(lower_points.end(),
463 jouk_points.begin() + nose_index,
464 jouk_points.begin() + tail_index + 1);
467 airfoil_1D[0] = make_points_equidistant(upper_points, number_points);
468 airfoil_1D[1] = make_points_equidistant(lower_points, number_points);
471 auto move_nose_to_origin = [](std::vector<Point<2>> &vector) {
472 const double nose_x_pos = vector.front()(0);
473 for (
auto &i : vector)
477 move_nose_to_origin(airfoil_1D[1]);
478 move_nose_to_origin(airfoil_1D[0]);
507 static std::vector<Point<2>>
509 const unsigned int number_points)
511 std::vector<Point<2>> circle_points;
525 radius_test < radius,
527 "Error creating lower circle: Circle for Joukowski-transform does"
528 " not enclose point zeta = -1! Choose different center "
532 const double theta = 2 *
numbers::PI / number_points;
534 for (
unsigned int i = 0; i < number_points; i++)
535 circle_points.emplace_back(
center[0] - radius *
cos(i * theta),
538 return circle_points;
549 static std::vector<Point<2>>
550 joukowski_transform(
const std::vector<
Point<2>> &circle_points)
552 std::vector<Point<2>> joukowski_points(circle_points.size());
555 for (
unsigned int i = 0; i < circle_points.size(); i++)
557 const double chi = circle_points[i](0);
558 const double eta = circle_points[i](1);
559 const std::complex<double> zeta(chi, eta);
560 const std::complex<double> z = zeta + 1. / zeta;
562 joukowski_points[i] = {real(z), imag(z)};
564 return joukowski_points;
583 static std::array<std::vector<Point<2>>, 2>
584 naca(
const std::string &serialnumber,
585 const unsigned int number_points,
586 const unsigned int factor)
590 const unsigned int n_airfoilpoints = factor * number_points;
593 return {{make_points_equidistant(
594 naca_create_points(serialnumber, n_airfoilpoints,
true),
596 make_points_equidistant(
597 naca_create_points(serialnumber, n_airfoilpoints,
false),
612 static std::vector<Point<2>>
613 naca_create_points(
const std::string &serialnumber,
614 const unsigned int number_points,
617 Assert(serialnumber.length() == 4,
618 ExcMessage(
"This NACA-serial number is not implemented!"));
620 return naca_create_points_4_digits(serialnumber,
639 static std::vector<Point<2>>
640 naca_create_points_4_digits(
const std::string &serialnumber,
641 const unsigned int number_points,
645 const unsigned int digit_0 = (serialnumber[0] -
'0');
646 const unsigned int digit_1 = (serialnumber[1] -
'0');
647 const unsigned int digit_2 = (serialnumber[2] -
'0');
648 const unsigned int digit_3 = (serialnumber[3] -
'0');
650 const unsigned int digit_23 = 10 * digit_2 + digit_3;
653 const double t =
static_cast<double>(digit_23) / 100.0;
655 std::vector<Point<2>> naca_points;
657 if (digit_0 == 0 && digit_1 == 0)
658 for (
unsigned int i = 0; i < number_points; i++)
660 const double x = i * 1 / (1.0 * number_points - 1);
663 (0.2969 *
std::pow(x, 0.5) - 0.126 * x -
668 naca_points.emplace_back(x, +y_t);
670 naca_points.emplace_back(x, -y_t);
673 for (
unsigned int i = 0; i < number_points; i++)
675 const double m = 1.0 * digit_0 / 100;
676 const double p = 1.0 * digit_1 / 10;
677 const double x = i * 1 / (1.0 * number_points - 1);
682 ((1 - 2 * p) + 2 * p * x -
std::pow(x, 2));
684 const double dy_c = (x <= p) ?
686 2 * m /
std::pow(1 - p, 2) * (p - x);
690 (0.2969 *
std::pow(x, 0.5) - 0.126 * x -
697 naca_points.emplace_back(x - y_t *
std::sin(theta),
700 naca_points.emplace_back(x + y_t *
std::sin(theta),
717 static std::array<std::vector<Point<2>>, 2>
718 set_airfoil_length(
const std::array<std::vector<
Point<2>>, 2> &input,
719 const double desired_len)
721 std::array<std::vector<Point<2>>, 2> output;
722 output[0] = set_airfoil_length(input[0], desired_len);
723 output[1] = set_airfoil_length(input[1], desired_len);
735 static std::vector<Point<2>>
736 set_airfoil_length(
const std::vector<
Point<2>> &input,
737 const double desired_len)
739 std::vector<Point<2>> output = input;
742 desired_len / input.front().distance(input.back());
744 for (
auto &x : output)
760 static std::vector<Point<2>>
761 make_points_equidistant(
762 const std::vector<
Point<2>> &non_equidistant_points,
763 const unsigned int number_points)
765 const unsigned int n_points =
766 non_equidistant_points
770 std::vector<double> arclength_L(non_equidistant_points.size(), 0);
771 for (
unsigned int i = 0; i < non_equidistant_points.size() - 1; i++)
774 non_equidistant_points[i + 1].distance(non_equidistant_points[i]);
777 const auto airfoil_length =
779 const auto deltaX = airfoil_length / (number_points - 1);
783 std::vector<Point<2>> equidist(
786 equidist[0] = non_equidistant_points[0];
787 equidist[number_points - 1] = non_equidistant_points[n_points - 1];
791 for (
unsigned int j = 0, i = 1; j < n_points - 1; j++)
794 const auto Lj = arclength_L[j];
795 const auto Ljp = arclength_L[j + 1];
797 while (Lj <= i * deltaX && i * deltaX <= Ljp &&
798 i < number_points - 1)
800 equidist[i] =
Point<2>((i * deltaX - Lj) / (Ljp - Lj) *
801 (non_equidistant_points[j + 1] -
802 non_equidistant_points[j]) +
803 non_equidistant_points[j]);
822 std::vector<Triangulation<2>> trias(10);
826 const std::vector<Point<2>> & corner_vertices,
827 const std::vector<unsigned int> &repetitions,
839 auto &
point = it->vertex();
840 const double xi =
point(0);
841 const double eta =
point(1);
844 point = 0.25 * ((1 - xi) * (1 - eta) * corner_vertices[0] +
845 (1 + xi) * (1 - eta) * corner_vertices[1] +
846 (1 - xi) * (1 + eta) * corner_vertices[2] +
847 (1 + xi) * (1 + eta) * corner_vertices[3]);
859 {n_subdivision_y, n_subdivision_x_0},
863 {n_subdivision_y, n_subdivision_x_0},
867 {n_subdivision_x_1, n_subdivision_y},
871 {n_subdivision_x_1, n_subdivision_y},
875 {n_subdivision_x_2, n_subdivision_y},
879 {n_subdivision_x_2, n_subdivision_y},
908 if (cell->face(f)->at_boundary() ==
false)
911 const auto mid = cell->material_id();
913 if ((mid == id_block_1 && f == 0) ||
914 (mid == id_block_4 && f == 0))
915 cell->face(f)->set_boundary_id(0);
916 else if ((mid == id_block_3 && f == 0) ||
917 (mid == id_block_6 && f == 2))
918 cell->face(f)->set_boundary_id(1);
919 else if ((mid == id_block_1 && f == 1) ||
920 (mid == id_block_2 && f == 1))
921 cell->face(f)->set_boundary_id(2);
922 else if ((mid == id_block_4 && f == 1) ||
923 (mid == id_block_5 && f == 3))
924 cell->face(f)->set_boundary_id(3);
925 else if ((mid == id_block_2 && f == 0) ||
926 (mid == id_block_3 && f == 2))
927 cell->face(f)->set_boundary_id(4);
928 else if ((mid == id_block_5 && f == 2) ||
929 (mid == id_block_6 && f == 0))
930 cell->face(f)->set_boundary_id(5);
949 std::vector<bool> vertex_processed(tria.
n_vertices(),
false);
959 rotation_matrix_2 =
transpose(rotation_matrix_1);
970 for (
auto &cell : tria)
974 if (vertex_processed[cell.vertex_index(v)])
978 vertex_processed[cell.vertex_index(v)] =
true;
980 auto &node = cell.vertex(v);
983 if (cell.material_id() == id_block_1 ||
984 cell.material_id() == id_block_4)
993 if (cell.material_id() == id_block_1)
996 (node - horizontal_offset) +
1002 else if (cell.material_id() == id_block_4)
1005 (node - horizontal_offset) -
1011 const double trapeze_height =
1017 const double x2 =
L - l_a - l_b;
1019 const double Dx = x1 + x2 + x3;
1020 const double deltax =
1022 const double dx = Dx / n_cells_x_0;
1023 const double dy = trapeze_height / n_cells_y;
1025 static_cast<int>(std::round((node_(0) - deltax) /
dx));
1027 static_cast<int>(std::round(
std::abs(node_(1)) / dy));
1029 node_(0) =
numbers::PI / 2 * (1.0 * ix) / n_cells_x_0;
1030 node_(1) = height * (1.0 * iy) / n_cells_y;
1037 const double dy = height / n_cells_y;
1039 static_cast<int>(std::round(node_(0) /
dx));
1041 static_cast<int>(std::round(node_(1) / dy));
1042 const double alpha =
1043 bias_alpha(1 - (1.0 * iy) / n_cells_y);
1044 const double theta = node_(0);
1046 ((cell.material_id() == id_block_1) ?
1052 (cell.material_id() == id_block_1) ? (0) : (1))][ix] *
1057 else if (cell.material_id() == id_block_2 ||
1058 cell.material_id() == id_block_5)
1066 "Points D,C,G and E,F,I are not defined symmetric to "
1067 "x-axis, which is required to interpolate block 2"
1068 " and 5 with same geometric computations."));
1069 const double l_y = D(1) -
C(1);
1070 const double l_h = D(1) - l_y;
1071 const double by = -l_h / length_b1_x * (node(0) - H(0));
1072 const double dy = (height - by) / n_cells_y;
1073 const int iy =
static_cast<int>(
1074 std::round((
std::abs(node(1)) - by) / dy));
1075 const double dx = length_b1_x / n_cells_x_1;
1076 const int ix =
static_cast<int>(
1077 std::round(
std::abs(node(0) - center_mesh) /
dx));
1079 const double alpha = bias_alpha(1 - (1.0 * iy) / n_cells_y);
1084 incline_factor * length_b2 * ix /
1086 ((cell.material_id() == id_block_2) ?
1092 (cell.material_id() == id_block_2) ? (0) : (1))]
1093 [n_cells_x_0 + ix] *
1097 else if (cell.material_id() == id_block_3 ||
1098 cell.material_id() == id_block_6)
1101 const double dx = length_b2 / n_cells_x_2;
1102 const double dy = height / n_cells_y;
1103 const int ix =
static_cast<int>(
1106 static_cast<int>(std::round(
std::abs(node(1)) / dy));
1108 const double alpha_y = bias_alpha(1 - 1.0 * iy / n_cells_y);
1109 const double alpha_x =
1110 bias_alpha(1 - (
static_cast<double>(ix)) / n_cells_x_2);
1114 const Point<2> p1(J(0) - (1 - incline_factor) * length_b2 *
1116 ((cell.material_id() == id_block_3) ?
1121 const Point<2> p2(J(0) - alpha_x * length_b2, tail_y);
1122 node = p1 * (1 - alpha_y) + p2 * alpha_y;
1142 bias_alpha(
double alpha)
const
1151 void internal_create_triangulation(
1155 const AdditionalData & additional_data)
1157 MeshGenerator mesh_generator(additional_data);
1160 if (
auto parallel_tria =
1163 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1164 else if (
auto parallel_tria =
dynamic_cast<
1167 mesh_generator.create_triangulation(*parallel_tria, periodic_faces);
1169 mesh_generator.create_triangulation(tria, periodic_faces);
1184 const AdditionalData &)
1193 const AdditionalData &additional_data)
1195 internal_create_triangulation(tria,
nullptr, additional_data);
1205 const AdditionalData & additional_data)
1207 internal_create_triangulation(tria, &periodic_faces, additional_data);
1217 const AdditionalData & additional_data)
1221 (void)additional_data;
1222 (void)periodic_faces;
1233 template <
int dim,
int spacedim>
1245 cell->face(f)->set_boundary_id(f);
1251 template <
int spacedim>
1261 if (cell->center()(0) > 0)
1262 cell->set_material_id(1);
1269 template <
int dim,
int spacedim>
1287 for (; face != endface; ++face)
1288 if (face->at_boundary())
1289 if (face->boundary_id() == 0)
1294 face->set_boundary_id(0);
1296 face->set_boundary_id(1);
1298 face->set_boundary_id(2);
1300 face->set_boundary_id(3);
1302 face->set_boundary_id(4);
1304 face->set_boundary_id(5);
1316 for (
unsigned int d = 0;
d < dim; ++
d)
1317 if (cell->center()(
d) > 0)
1319 cell->set_material_id(
id);
1344 cell->face(2)->set_all_boundary_ids(1);
1367 cell->face(4)->set_all_boundary_ids(1);
1371 cell->face(2)->set_all_boundary_ids(1);
1375 cell->face(2)->set_all_boundary_ids(1);
1379 cell->face(0)->set_all_boundary_ids(1);
1383 cell->face(2)->set_all_boundary_ids(1);
1387 cell->face(0)->set_all_boundary_ids(1);
1389 else if (tria.
n_cells() == 12)
1398 cell->face(5)->set_all_boundary_ids(1);
1401 else if (tria.
n_cells() == 96)
1414 unsigned int count = 0;
1418 if (cell->face(5)->at_boundary())
1420 cell->face(5)->set_all_boundary_ids(1);
1438 const double inner_radius,
1439 const double outer_radius)
1444 double middle = (outer_radius - inner_radius) / 2e0 + inner_radius;
1445 double eps = 1
e-3 * middle;
1448 for (; cell != tria.
end(); ++cell)
1451 if (!cell->face(f)->at_boundary())
1454 double radius = cell->face(f)->center().norm() -
center.
norm();
1455 if (
std::fabs(cell->face(f)->center()(0)) <
1458 cell->face(f)->set_boundary_id(2);
1459 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
1461 if (cell->face(f)->line(j)->at_boundary())
1462 if (
std::fabs(cell->face(f)->line(j)->vertex(0).norm() -
1463 cell->face(f)->line(j)->vertex(1).norm()) >
1465 cell->face(f)->line(j)->set_boundary_id(2);
1467 else if (
std::fabs(cell->face(f)->center()(1)) <
1470 cell->face(f)->set_boundary_id(3);
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(3);
1479 else if (
std::fabs(cell->face(f)->center()(2)) <
1482 cell->face(f)->set_boundary_id(4);
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(4);
1491 else if (radius < middle)
1493 cell->face(f)->set_boundary_id(0);
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(0);
1502 else if (radius > middle)
1504 cell->face(f)->set_boundary_id(1);
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(1);
1521 template <
int dim,
int spacedim>
1532 for (
unsigned int i = 0; i < dim; ++i)
1572 std::vector<CellData<dim>> cells(1);
1574 cells[0].vertices[i] = i;
1575 cells[0].material_id = 0;
1581 colorize_hyper_rectangle(tria);
1586 template <
int dim,
int spacedim>
1594 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1597 for (
unsigned int i = 0; i < dim; ++i)
1617 for (
unsigned int d = 0;
d < dim; ++
d)
1618 for (
unsigned int c = 1; c <= dim; ++c)
1621 ExcMessage(
"Vertices of simplex must form a right handed system"));
1625 std::vector<Point<dim>> points =
vertices;
1629 for (
unsigned int i = 0; i <= dim; ++i)
1631 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1637 for (
unsigned int i = 1; i < dim; ++i)
1638 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1640 for (
unsigned int i = 0; i <= dim; ++i)
1641 points.push_back(1. / 3. *
1642 (points[i] + points[(i + 1) % (dim + 1)] +
1643 points[(i + 2) % (dim + 1)]));
1645 points.push_back((1. / (dim + 1)) *
center);
1647 std::vector<CellData<dim>> cells(dim + 1);
1652 cells[0].vertices[0] = 0;
1653 cells[0].vertices[1] = 3;
1654 cells[0].vertices[2] = 5;
1655 cells[0].vertices[3] = 6;
1656 cells[0].material_id = 0;
1658 cells[1].vertices[0] = 3;
1659 cells[1].vertices[1] = 1;
1660 cells[1].vertices[2] = 6;
1661 cells[1].vertices[3] = 4;
1662 cells[1].material_id = 0;
1664 cells[2].vertices[0] = 5;
1665 cells[2].vertices[1] = 6;
1666 cells[2].vertices[2] = 2;
1667 cells[2].vertices[3] = 4;
1668 cells[2].material_id = 0;
1672 cells[0].vertices[0] = 0;
1673 cells[0].vertices[1] = 4;
1674 cells[0].vertices[2] = 8;
1675 cells[0].vertices[3] = 10;
1676 cells[0].vertices[4] = 7;
1677 cells[0].vertices[5] = 13;
1678 cells[0].vertices[6] = 12;
1679 cells[0].vertices[7] = 14;
1680 cells[0].material_id = 0;
1682 cells[1].vertices[0] = 4;
1683 cells[1].vertices[1] = 1;
1684 cells[1].vertices[2] = 10;
1685 cells[1].vertices[3] = 5;
1686 cells[1].vertices[4] = 13;
1687 cells[1].vertices[5] = 9;
1688 cells[1].vertices[6] = 14;
1689 cells[1].vertices[7] = 11;
1690 cells[1].material_id = 0;
1692 cells[2].vertices[0] = 8;
1693 cells[2].vertices[1] = 10;
1694 cells[2].vertices[2] = 2;
1695 cells[2].vertices[3] = 5;
1696 cells[2].vertices[4] = 12;
1697 cells[2].vertices[5] = 14;
1698 cells[2].vertices[6] = 6;
1699 cells[2].vertices[7] = 11;
1700 cells[2].material_id = 0;
1702 cells[3].vertices[0] = 7;
1703 cells[3].vertices[1] = 13;
1704 cells[3].vertices[2] = 12;
1705 cells[3].vertices[3] = 14;
1706 cells[3].vertices[4] = 3;
1707 cells[3].vertices[5] = 9;
1708 cells[3].vertices[6] = 6;
1709 cells[3].vertices[7] = 11;
1710 cells[3].material_id = 0;
1720 template <
int dim,
int spacedim>
1733 const std::vector<Point<spacedim>>
vertices = {
1739 std::vector<CellData<dim>> cells(1);
1740 cells[0].vertices = {0, 1, 2};
1748 static const std::vector<Point<spacedim>>
vertices = {
1749 {{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}}};
1751 std::vector<CellData<dim>> cells(1);
1752 cells[0].vertices = {0, 1, 2, 3};
1760 static const std::vector<Point<spacedim>>
vertices = {
1765 {+0.0, +0.0, 1.0}}};
1767 std::vector<CellData<dim>> cells(1);
1768 cells[0].vertices = {0, 1, 2, 3, 4};
1776 static const std::vector<Point<spacedim>>
vertices = {
1784 std::vector<CellData<dim>> cells(1);
1785 cells[0].vertices = {0, 1, 2, 3, 4, 5};
1797 const unsigned int n_rotations,
1801 const unsigned int dim = 3;
1804 "More than 4 cells are needed to create a moebius grid."));
1806 ExcMessage(
"Outer and inner radius must be positive."));
1808 ExcMessage(
"Outer radius must be greater than inner radius."));
1815 for (
unsigned int i = 0; i <
n_cells; ++i)
1816 for (
unsigned int j = 0; j < 4; ++j)
1830 unsigned int offset = 0;
1835 std::vector<CellData<dim>> cells(
n_cells);
1836 for (
unsigned int i = 0; i <
n_cells; ++i)
1838 for (
unsigned int j = 0; j < 2; ++j)
1850 cells[i].material_id = 0;
1855 (0 + n_rotations) % 4;
1857 (3 + n_rotations) % 4;
1859 (2 + n_rotations) % 4;
1861 (1 + n_rotations) % 4;
1877 ExcMessage(
"Outer radius R must be greater than the inner "
1881 const unsigned int dim = 2;
1882 const unsigned int spacedim = 3;
1883 std::vector<Point<spacedim>>
vertices(16);
1902 std::vector<CellData<dim>> cells(16);
1904 cells[0].vertices[0] = 0;
1905 cells[0].vertices[1] = 4;
1906 cells[0].vertices[2] = 3;
1907 cells[0].vertices[3] = 7;
1908 cells[0].material_id = 0;
1910 cells[1].vertices[0] = 1;
1911 cells[1].vertices[1] = 5;
1912 cells[1].vertices[2] = 0;
1913 cells[1].vertices[3] = 4;
1914 cells[1].material_id = 0;
1916 cells[2].vertices[0] = 2;
1917 cells[2].vertices[1] = 6;
1918 cells[2].vertices[2] = 1;
1919 cells[2].vertices[3] = 5;
1920 cells[2].material_id = 0;
1922 cells[3].vertices[0] = 3;
1923 cells[3].vertices[1] = 7;
1924 cells[3].vertices[2] = 2;
1925 cells[3].vertices[3] = 6;
1926 cells[3].material_id = 0;
1928 cells[4].vertices[0] = 4;
1929 cells[4].vertices[1] = 8;
1930 cells[4].vertices[2] = 7;
1931 cells[4].vertices[3] = 11;
1932 cells[4].material_id = 0;
1934 cells[5].vertices[0] = 5;
1935 cells[5].vertices[1] = 9;
1936 cells[5].vertices[2] = 4;
1937 cells[5].vertices[3] = 8;
1938 cells[5].material_id = 0;
1940 cells[6].vertices[0] = 6;
1941 cells[6].vertices[1] = 10;
1942 cells[6].vertices[2] = 5;
1943 cells[6].vertices[3] = 9;
1944 cells[6].material_id = 0;
1946 cells[7].vertices[0] = 7;
1947 cells[7].vertices[1] = 11;
1948 cells[7].vertices[2] = 6;
1949 cells[7].vertices[3] = 10;
1950 cells[7].material_id = 0;
1952 cells[8].vertices[0] = 8;
1953 cells[8].vertices[1] = 12;
1954 cells[8].vertices[2] = 11;
1955 cells[8].vertices[3] = 15;
1956 cells[8].material_id = 0;
1958 cells[9].vertices[0] = 9;
1959 cells[9].vertices[1] = 13;
1960 cells[9].vertices[2] = 8;
1961 cells[9].vertices[3] = 12;
1962 cells[9].material_id = 0;
1964 cells[10].vertices[0] = 10;
1965 cells[10].vertices[1] = 14;
1966 cells[10].vertices[2] = 9;
1967 cells[10].vertices[3] = 13;
1968 cells[10].material_id = 0;
1970 cells[11].vertices[0] = 11;
1971 cells[11].vertices[1] = 15;
1972 cells[11].vertices[2] = 10;
1973 cells[11].vertices[3] = 14;
1974 cells[11].material_id = 0;
1976 cells[12].vertices[0] = 12;
1977 cells[12].vertices[1] = 0;
1978 cells[12].vertices[2] = 15;
1979 cells[12].vertices[3] = 3;
1980 cells[12].material_id = 0;
1982 cells[13].vertices[0] = 13;
1983 cells[13].vertices[1] = 1;
1984 cells[13].vertices[2] = 12;
1985 cells[13].vertices[3] = 0;
1986 cells[13].material_id = 0;
1988 cells[14].vertices[0] = 14;
1989 cells[14].vertices[1] = 2;
1990 cells[14].vertices[2] = 13;
1991 cells[14].vertices[3] = 1;
1992 cells[14].material_id = 0;
1994 cells[15].vertices[0] = 15;
1995 cells[15].vertices[1] = 3;
1996 cells[15].vertices[2] = 14;
1997 cells[15].vertices[3] = 2;
1998 cells[15].material_id = 0;
2013 const unsigned int n_cells_toroidal,
2017 ExcMessage(
"Outer radius R must be greater than the inner "
2020 Assert(n_cells_toroidal > 2,
2021 ExcMessage(
"Number of cells in toroidal direction has "
2022 "to be at least 3."));
2028 double const a = 1. / (1 +
std::sqrt(2.0));
2031 const unsigned int additional_layer =
2035 const unsigned int n_point_layers_toroidal =
2036 n_cells_toroidal + additional_layer;
2037 std::vector<Point<3>>
vertices(8 * n_point_layers_toroidal);
2049 double const phi_cell = phi / n_cells_toroidal;
2050 for (
unsigned int c = 1; c < n_point_layers_toroidal; ++c)
2052 for (
unsigned int v = 0; v < 8; ++v)
2054 double const r_2d =
vertices[v][0];
2062 std::vector<CellData<3>> cells(5 * n_cells_toroidal);
2063 for (
unsigned int c = 0; c < n_cells_toroidal; ++c)
2065 for (
unsigned int j = 0; j < 2; ++j)
2067 const unsigned int offset =
2068 (8 * (c + j)) % (8 * n_point_layers_toroidal);
2071 cells[5 * c].vertices[0 + j * 4] = offset + 0;
2072 cells[5 * c].vertices[1 + j * 4] = offset + 1;
2073 cells[5 * c].vertices[2 + j * 4] = offset + 2;
2074 cells[5 * c].vertices[3 + j * 4] = offset + 3;
2076 cells[5 * c + 1].vertices[0 + j * 4] = offset + 2;
2077 cells[5 * c + 1].vertices[1 + j * 4] = offset + 3;
2078 cells[5 * c + 1].vertices[2 + j * 4] = offset + 4;
2079 cells[5 * c + 1].vertices[3 + j * 4] = offset + 5;
2081 cells[5 * c + 2].vertices[0 + j * 4] = offset + 4;
2082 cells[5 * c + 2].vertices[1 + j * 4] = offset + 5;
2083 cells[5 * c + 2].vertices[2 + j * 4] = offset + 6;
2084 cells[5 * c + 2].vertices[3 + j * 4] = offset + 7;
2086 cells[5 * c + 3].vertices[0 + j * 4] = offset + 0;
2087 cells[5 * c + 3].vertices[1 + j * 4] = offset + 2;
2088 cells[5 * c + 3].vertices[2 + j * 4] = offset + 6;
2089 cells[5 * c + 3].vertices[3 + j * 4] = offset + 4;
2091 cells[5 * c + 4].vertices[0 + j * 4] = offset + 3;
2092 cells[5 * c + 4].vertices[1 + j * 4] = offset + 1;
2093 cells[5 * c + 4].vertices[2 + j * 4] = offset + 5;
2094 cells[5 * c + 4].vertices[3 + j * 4] = offset + 7;
2097 cells[5 * c].material_id = 0;
2099 cells[5 * c + 1].material_id = 1;
2100 cells[5 * c + 2].material_id = 0;
2101 cells[5 * c + 3].material_id = 0;
2102 cells[5 * c + 4].material_id = 0;
2117 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2119 cell->face(f)->set_all_manifold_ids(1);
2124 if (cell->material_id() == 1)
2126 cell->set_all_manifold_ids(2);
2128 cell->set_material_id(0);
2143 template <
int dim,
int spacedim>
2164 "The volume of the cell is not greater than zero. "
2165 "This could be due to the wrong ordering of the vertices."));
2193 std::array<Tensor<1, 2>, 2> edges;
2194 edges[0] = corners[0];
2195 edges[1] = corners[1];
2196 std::vector<unsigned int> subdivisions;
2197 subdivided_parallelepiped<2, 2>(
2198 tria, origin, edges, subdivisions,
colorize);
2209 unsigned int n_subdivisions[dim];
2210 for (
unsigned int i = 0; i < dim; ++i)
2211 n_subdivisions[i] = 1;
2220 const unsigned int n_subdivisions,
2226 unsigned int n_subdivisions_[dim];
2227 for (
unsigned int i = 0; i < dim; ++i)
2228 n_subdivisions_[i] = n_subdivisions;
2238 const unsigned int (&n_subdivisions)[dim],
2240 const unsigned int *n_subdivisions,
2246 std::vector<unsigned int> subdivisions;
2247 std::array<Tensor<1, dim>, dim> edges;
2248 for (
unsigned int i = 0; i < dim; ++i)
2250 subdivisions.push_back(n_subdivisions[i]);
2251 edges[i] = corners[i];
2254 subdivided_parallelepiped<dim, dim>(
2255 tria, origin, edges, subdivisions,
colorize);
2261 template <
int dim,
int spacedim>
2266 const std::vector<unsigned int> &subdivisions,
2269 std::vector<unsigned int> compute_subdivisions = subdivisions;
2270 if (compute_subdivisions.size() == 0)
2272 compute_subdivisions.resize(dim, 1);
2275 Assert(compute_subdivisions.size() == dim,
2276 ExcMessage(
"One subdivision must be provided for each dimension."));
2278 for (
unsigned int i = 0; i < dim; ++i)
2280 Assert(compute_subdivisions[i] > 0,
2283 edges[i].
norm() > 0,
2285 "Edges in subdivided_parallelepiped() must not be degenerate."));
2293 bool twisted_data =
false;
2298 twisted_data = (edges[0][0] < 0);
2305 const double plane_normal =
2306 edges[0][0] * edges[1][1] - edges[0][1] * edges[1][0];
2307 twisted_data = (plane_normal < 0.0);
2316 (edges[0].
norm() * edges[1].
norm()) -
2319 "Edges in subdivided_parallelepiped() must point in"
2320 " different directions."));
2336 twisted_data = (plane_normal * edges[2] < 0.0);
2346 "The triangulation you are trying to create will consist of cells"
2347 " with negative measures. This is usually the result of input data"
2348 " that does not define a right-handed coordinate system. The usual"
2349 " fix for this is to ensure that in 1D the given point is to the"
2350 " right of the origin (or the given edge tensor is positive), in 2D"
2351 " that the two edges (and their cross product) obey the right-hand"
2352 " rule (which may usually be done by switching the order of the"
2353 " points or edge tensors), or in 3D that the edges form a"
2354 " right-handed coordinate system (which may also be accomplished by"
2355 " switching the order of the first two points or edge tensors)."));
2358 for (
unsigned int i = 0; i < dim; ++i)
2359 for (
unsigned int j = i + 1; j < dim; ++j)
2360 Assert((edges[i] != edges[j]),
2362 "Degenerate edges of subdivided_parallelepiped encountered."));
2365 std::vector<Point<spacedim>> points;
2370 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2371 points.push_back(origin + edges[0] / compute_subdivisions[0] * x);
2375 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2376 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2377 points.push_back(origin + edges[0] / compute_subdivisions[0] * x +
2378 edges[1] / compute_subdivisions[1] * y);
2382 for (
unsigned int z = 0; z <= compute_subdivisions[2]; ++z)
2383 for (
unsigned int y = 0; y <= compute_subdivisions[1]; ++y)
2384 for (
unsigned int x = 0; x <= compute_subdivisions[0]; ++x)
2385 points.push_back(origin +
2386 edges[0] / compute_subdivisions[0] * x +
2387 edges[1] / compute_subdivisions[1] * y +
2388 edges[2] / compute_subdivisions[2] * z);
2397 for (
unsigned int i = 0; i < dim; ++i)
2398 n_cells *= compute_subdivisions[i];
2399 std::vector<CellData<dim>> cells(
n_cells);
2405 for (
unsigned int x = 0; x < compute_subdivisions[0]; ++x)
2407 cells[x].vertices[0] = x;
2408 cells[x].vertices[1] = x + 1;
2411 cells[x].material_id = 0;
2418 const unsigned int n_dy = compute_subdivisions[1];
2419 const unsigned int n_dx = compute_subdivisions[0];
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 = y * n_dx + x;
2425 cells[c].vertices[0] = y * (n_dx + 1) + x;
2426 cells[c].vertices[1] = y * (n_dx + 1) + x + 1;
2427 cells[c].vertices[2] = (y + 1) * (n_dx + 1) + x;
2428 cells[c].vertices[3] = (y + 1) * (n_dx + 1) + x + 1;
2431 cells[c].material_id = 0;
2439 const unsigned int n_dz = compute_subdivisions[2];
2440 const unsigned int n_dy = compute_subdivisions[1];
2441 const unsigned int n_dx = compute_subdivisions[0];
2443 for (
unsigned int z = 0; z < n_dz; ++z)
2444 for (
unsigned int y = 0; y < n_dy; ++y)
2445 for (
unsigned int x = 0; x < n_dx; ++x)
2447 const unsigned int c = z * n_dy * n_dx + y * n_dx + x;
2449 cells[c].vertices[0] =
2450 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2451 cells[c].vertices[1] =
2452 z * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x + 1;
2453 cells[c].vertices[2] =
2454 z * (n_dy + 1) * (n_dx + 1) + (y + 1) * (n_dx + 1) + x;
2455 cells[c].vertices[3] = z * (n_dy + 1) * (n_dx + 1) +
2456 (y + 1) * (n_dx + 1) + x + 1;
2457 cells[c].vertices[4] =
2458 (z + 1) * (n_dy + 1) * (n_dx + 1) + y * (n_dx + 1) + x;
2459 cells[c].vertices[5] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2460 y * (n_dx + 1) + x + 1;
2461 cells[c].vertices[6] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2462 (y + 1) * (n_dx + 1) + x;
2463 cells[c].vertices[7] = (z + 1) * (n_dy + 1) * (n_dx + 1) +
2464 (y + 1) * (n_dx + 1) + x + 1;
2467 cells[c].material_id = 0;
2488 for (; cell != endc; ++cell)
2492 if (cell->face(face)->at_boundary())
2493 cell->face(face)->set_boundary_id(face);
2500 template <
int dim,
int spacedim>
2503 const unsigned int repetitions,
2510 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2513 for (
unsigned int i = 0; i < dim; ++i)
2519 std::vector<unsigned int> reps(dim, repetitions);
2525 template <
int dim,
int spacedim>
2528 const std::vector<unsigned int> &repetitions,
2539 for (
unsigned int i = 0; i < dim; ++i)
2546 std::vector<Point<spacedim>> delta(dim);
2547 for (
unsigned int i = 0; i < dim; ++i)
2551 delta[i][i] = (p2[i] - p1[i]) / repetitions[i];
2555 "The first dim entries of coordinates of p1 and p2 need to be different."));
2559 std::vector<Point<spacedim>> points;
2563 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2564 points.push_back(p1 + x * delta[0]);
2568 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2569 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2570 points.push_back(p1 + x * delta[0] + y * delta[1]);
2574 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
2575 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2576 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2577 points.push_back(p1 + x * delta[0] + y * delta[1] +
2586 std::vector<CellData<dim>> cells;
2591 cells.resize(repetitions[0]);
2592 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2594 cells[x].vertices[0] = x;
2595 cells[x].vertices[1] = x + 1;
2596 cells[x].material_id = 0;
2603 cells.resize(repetitions[1] * repetitions[0]);
2604 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2605 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2607 const unsigned int c = x + y * repetitions[0];
2608 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
2609 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
2610 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
2611 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
2612 cells[c].material_id = 0;
2619 const unsigned int n_x = (repetitions[0] + 1);
2620 const unsigned int n_xy =
2621 (repetitions[0] + 1) * (repetitions[1] + 1);
2623 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
2624 for (
unsigned int z = 0; z < repetitions[2]; ++z)
2625 for (
unsigned int y = 0; y < repetitions[1]; ++y)
2626 for (
unsigned int x = 0; x < repetitions[0]; ++x)
2628 const unsigned int c = x + y * repetitions[0] +
2629 z * repetitions[0] * repetitions[1];
2630 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2631 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2632 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2633 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2634 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2635 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2636 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2637 cells[c].vertices[7] =
2638 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2639 cells[c].material_id = 0;
2661 for (
unsigned int i = 0; i < dim; ++i)
2665 "The distance between corner points must be positive."))
2669 colorize_subdivided_hyper_rectangle(tria, p1, p2,
epsilon);
2678 const
std::vector<
std::vector<
double>> &step_sz,
2679 const
Point<dim> & p_1,
2680 const
Point<dim> & p_2,
2693 std::vector<std::vector<double>> step_sizes(step_sz);
2695 for (
unsigned int i = 0; i < dim; ++i)
2700 std::reverse(step_sizes[i].
begin(), step_sizes[i].
end());
2704 for (
unsigned int j = 0; j < step_sizes.at(i).size(); j++)
2705 x += step_sizes[i][j];
2708 "The sequence of step sizes in coordinate direction " +
2710 " must be equal to the distance of the two given "
2711 "points in this coordinate direction."));
2717 std::vector<Point<dim>> points;
2723 for (
unsigned int i = 0;; ++i)
2732 if (i == step_sizes[0].size())
2735 x += step_sizes[0][i];
2743 for (
unsigned int j = 0;; ++j)
2746 for (
unsigned int i = 0;; ++i)
2748 points.push_back(
Point<dim>(p1[0] + x, p1[1] + y));
2749 if (i == step_sizes[0].size())
2752 x += step_sizes[0][i];
2755 if (j == step_sizes[1].size())
2758 y += step_sizes[1][j];
2765 for (
unsigned int k = 0;; ++k)
2768 for (
unsigned int j = 0;; ++j)
2771 for (
unsigned int i = 0;; ++i)
2774 Point<dim>(p1[0] + x, p1[1] + y, p1[2] + z));
2775 if (i == step_sizes[0].size())
2778 x += step_sizes[0][i];
2781 if (j == step_sizes[1].size())
2784 y += step_sizes[1][j];
2787 if (k == step_sizes[2].size())
2790 z += step_sizes[2][k];
2801 std::vector<CellData<dim>> cells;
2806 cells.resize(step_sizes[0].size());
2807 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2809 cells[x].vertices[0] = x;
2810 cells[x].vertices[1] = x + 1;
2811 cells[x].material_id = 0;
2818 cells.resize(step_sizes[1].size() * step_sizes[0].size());
2819 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2820 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2822 const unsigned int c = x + y * step_sizes[0].size();
2823 cells[c].vertices[0] = y * (step_sizes[0].size() + 1) + x;
2824 cells[c].vertices[1] = y * (step_sizes[0].size() + 1) + x + 1;
2825 cells[c].vertices[2] =
2826 (y + 1) * (step_sizes[0].size() + 1) + x;
2827 cells[c].vertices[3] =
2828 (y + 1) * (step_sizes[0].size() + 1) + x + 1;
2829 cells[c].material_id = 0;
2836 const unsigned int n_x = (step_sizes[0].size() + 1);
2837 const unsigned int n_xy =
2838 (step_sizes[0].size() + 1) * (step_sizes[1].size() + 1);
2840 cells.resize(step_sizes[2].size() * step_sizes[1].size() *
2841 step_sizes[0].size());
2842 for (
unsigned int z = 0; z < step_sizes[2].size(); ++z)
2843 for (
unsigned int y = 0; y < step_sizes[1].size(); ++y)
2844 for (
unsigned int x = 0; x < step_sizes[0].size(); ++x)
2846 const unsigned int c =
2847 x + y * step_sizes[0].size() +
2848 z * step_sizes[0].size() * step_sizes[1].size();
2849 cells[c].vertices[0] = z * n_xy + y * n_x + x;
2850 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
2851 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
2852 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
2853 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
2854 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
2855 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
2856 cells[c].vertices[7] =
2857 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
2858 cells[c].material_id = 0;
2867 tria.create_triangulation(points, cells,
SubCellData());
2880 *std::min_element(step_sizes[0].
begin(), step_sizes[0].
end());
2881 for (
unsigned int i = 1; i < dim; ++i)
2883 *std::min_element(step_sizes[i].
begin(),
2884 step_sizes[i].
end()));
2885 const double epsilon = 0.01 * min_size;
2889 colorize_subdivided_hyper_rectangle(tria, p1, p2,
epsilon);
2898 const std::vector<std::vector<double>> &spacing,
2910 for (
unsigned int i = 0; i <
n_cells; i++)
2913 delta =
std::min(delta, spacing[0][i]);
2917 std::vector<Point<1>> points;
2919 for (
unsigned int x = 0; x <=
n_cells; ++x)
2921 points.emplace_back(ax);
2923 ax += spacing[0][x];
2926 unsigned int n_val_cells = 0;
2927 for (
unsigned int i = 0; i <
n_cells; i++)
2931 std::vector<CellData<1>> cells(n_val_cells);
2932 unsigned int id = 0;
2933 for (
unsigned int x = 0; x <
n_cells; ++x)
2936 cells[id].vertices[0] = x;
2937 cells[id].vertices[1] = x + 1;
2956 const std::vector<std::vector<double>> &spacing,
2963 std::vector<unsigned int> repetitions(2);
2966 for (
unsigned int i = 0; i < 2; i++)
2968 repetitions[i] = spacing[i].size();
2970 for (
unsigned int j = 0; j < repetitions[i]; j++)
2973 delta =
std::min(delta, spacing[i][j]);
2980 std::vector<Point<2>> points;
2982 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
2985 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
2987 points.emplace_back(ax, ay);
2988 if (x < repetitions[0])
2989 ax += spacing[0][x];
2991 if (y < repetitions[1])
2992 ay += spacing[1][y];
2996 unsigned int n_val_cells = 0;
2997 for (
unsigned int i = 0; i <
material_id.size(0); i++)
2998 for (
unsigned int j = 0; j <
material_id.size(1); j++)
3002 std::vector<CellData<2>> cells(n_val_cells);
3003 unsigned int id = 0;
3004 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3005 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3008 cells[id].vertices[0] = y * (repetitions[0] + 1) + x;
3009 cells[id].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3010 cells[id].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3011 cells[id].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3025 double eps = 0.01 * delta;
3027 for (; cell != endc; ++cell)
3029 Point<2> cell_center = cell->center();
3031 if (cell->face(f)->boundary_id() == 0)
3033 Point<2> face_center = cell->face(f)->center();
3034 for (
unsigned int i = 0; i < 2; ++i)
3036 if (face_center[i] < cell_center[i] -
eps)
3037 cell->face(f)->set_boundary_id(i * 2);
3038 if (face_center[i] > cell_center[i] +
eps)
3039 cell->face(f)->set_boundary_id(i * 2 + 1);
3050 const std::vector<std::vector<double>> &spacing,
3055 const unsigned int dim = 3;
3059 std::vector<unsigned int> repetitions(dim);
3062 for (
unsigned int i = 0; i < dim; i++)
3064 repetitions[i] = spacing[i].size();
3066 for (
unsigned int j = 0; j < repetitions[i]; j++)
3069 delta =
std::min(delta, spacing[i][j]);
3076 std::vector<Point<dim>> points;
3078 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3081 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3084 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3086 points.emplace_back(ax, ay, az);
3087 if (x < repetitions[0])
3088 ax += spacing[0][x];
3090 if (y < repetitions[1])
3091 ay += spacing[1][y];
3093 if (z < repetitions[2])
3094 az += spacing[2][z];
3098 unsigned int n_val_cells = 0;
3099 for (
unsigned int i = 0; i <
material_id.size(0); i++)
3100 for (
unsigned int j = 0; j <
material_id.size(1); j++)
3101 for (
unsigned int k = 0; k <
material_id.size(2); k++)
3105 std::vector<CellData<dim>> cells(n_val_cells);
3106 unsigned int id = 0;
3107 const unsigned int n_x = (repetitions[0] + 1);
3108 const unsigned int n_xy = (repetitions[0] + 1) * (repetitions[1] + 1);
3109 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3110 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3111 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3114 cells[id].vertices[0] = z * n_xy + y * n_x + x;
3115 cells[id].vertices[1] = z * n_xy + y * n_x + x + 1;
3116 cells[id].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3117 cells[id].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3118 cells[id].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3119 cells[id].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3120 cells[id].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3121 cells[id].vertices[7] = (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3135 double eps = 0.01 * delta;
3138 for (; cell != endc; ++cell)
3142 if (cell->face(f)->boundary_id() == 0)
3144 Point<dim> face_center = cell->face(f)->center();
3145 for (
unsigned int i = 0; i < dim; ++i)
3147 if (face_center[i] < cell_center[i] -
eps)
3148 cell->face(f)->set_boundary_id(i * 2);
3149 if (face_center[i] > cell_center[i] +
eps)
3150 cell->face(f)->set_boundary_id(i * 2 + 1);
3157 template <
int dim,
int spacedim>
3160 const std::vector<unsigned int> &holes)
3169 for (
unsigned int d = 0;
d < dim; ++
d)
3176 std::vector<Point<spacedim>> delta(dim);
3177 unsigned int repetitions[dim];
3178 for (
unsigned int i = 0; i < dim; ++i)
3181 ExcMessage(
"At least one hole needed in each direction"));
3182 repetitions[i] = 2 * holes[i] + 1;
3183 delta[i][i] = (p2[i] - p1[i]);
3188 std::vector<Point<spacedim>> points;
3192 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3193 points.push_back(p1 + x * delta[0]);
3197 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3198 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3199 points.push_back(p1 + x * delta[0] + y * delta[1]);
3203 for (
unsigned int z = 0; z <= repetitions[2]; ++z)
3204 for (
unsigned int y = 0; y <= repetitions[1]; ++y)
3205 for (
unsigned int x = 0; x <= repetitions[0]; ++x)
3206 points.push_back(p1 + x * delta[0] + y * delta[1] +
3216 std::vector<CellData<dim>> cells;
3221 cells.resize(repetitions[1] * repetitions[0] - holes[1] * holes[0]);
3223 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3224 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3226 if ((x % 2 == 1) && (y % 2 == 1))
3229 cells[c].vertices[0] = y * (repetitions[0] + 1) + x;
3230 cells[c].vertices[1] = y * (repetitions[0] + 1) + x + 1;
3231 cells[c].vertices[2] = (y + 1) * (repetitions[0] + 1) + x;
3232 cells[c].vertices[3] = (y + 1) * (repetitions[0] + 1) + x + 1;
3233 cells[c].material_id = 0;
3241 const unsigned int n_x = (repetitions[0] + 1);
3242 const unsigned int n_xy =
3243 (repetitions[0] + 1) * (repetitions[1] + 1);
3245 cells.resize(repetitions[2] * repetitions[1] * repetitions[0]);
3248 for (
unsigned int z = 0; z < repetitions[2]; ++z)
3249 for (
unsigned int y = 0; y < repetitions[1]; ++y)
3250 for (
unsigned int x = 0; x < repetitions[0]; ++x)
3253 cells[c].vertices[0] = z * n_xy + y * n_x + x;
3254 cells[c].vertices[1] = z * n_xy + y * n_x + x + 1;
3255 cells[c].vertices[2] = z * n_xy + (y + 1) * n_x + x;
3256 cells[c].vertices[3] = z * n_xy + (y + 1) * n_x + x + 1;
3257 cells[c].vertices[4] = (z + 1) * n_xy + y * n_x + x;
3258 cells[c].vertices[5] = (z + 1) * n_xy + y * n_x + x + 1;
3259 cells[c].vertices[6] = (z + 1) * n_xy + (y + 1) * n_x + x;
3260 cells[c].vertices[7] =
3261 (z + 1) * n_xy + (y + 1) * n_x + x + 1;
3262 cells[c].material_id = 0;
3289 const unsigned int ,
3300 const unsigned int ,
3312 bool inline point_in_2d_box(
const Point<2> &p,
3314 const double radius)
3316 return (
std::abs(p[0] - c[0]) < radius) &&
3325 template <
int dim,
int spacedim>
3330 for (
const auto &cell :
triangulation.active_cell_iterators())
3331 for (
unsigned int n = 0; n < GeometryInfo<dim>::lines_per_cell; ++n)
3332 length =
std::min(length, cell->line(n)->diameter());
3341 const double inner_radius,
3342 const double outer_radius,
3343 const double pad_bottom,
3344 const double pad_top,
3345 const double pad_left,
3346 const double pad_right,
3351 const unsigned int ,
3354 const bool with_padding =
3355 pad_bottom > 0 || pad_top > 0 || pad_left > 0 || pad_right > 0;
3367 for (
unsigned int n = 0; n < GeometryInfo<2>::lines_per_cell; ++n)
3368 length =
std::min(length, cell->line(n)->diameter());
3374 Triangulation<2> &cylinder_tria = with_padding ? cylinder_tria_maybe : tria;
3384 cell->set_manifold_id(tfi_manifold_id);
3386 const Point<2> bl(-outer_radius - pad_left, -outer_radius - pad_bottom);
3387 const Point<2> tr(outer_radius + pad_right, outer_radius + pad_top);
3393 auto add_sizes = [](std::vector<double> &step_sizes,
3394 const double padding,
3395 const double h) ->
void {
3398 const auto rounded =
3399 static_cast<unsigned int>(std::round(padding / h));
3402 const unsigned int num = (padding > 0. && rounded == 0) ? 1 : rounded;
3403 for (
unsigned int i = 0; i < num; ++i)
3404 step_sizes.push_back(padding / num);
3407 std::vector<std::vector<double>> step_sizes(2);
3410 add_sizes(step_sizes[0], pad_left, outer_radius);
3412 step_sizes[0].push_back(outer_radius);
3413 step_sizes[0].push_back(outer_radius);
3415 add_sizes(step_sizes[0], pad_right, outer_radius);
3418 add_sizes(step_sizes[1], pad_bottom, outer_radius);
3420 step_sizes[1].push_back(outer_radius);
3421 step_sizes[1].push_back(outer_radius);
3423 add_sizes(step_sizes[1], pad_top, outer_radius);
3428 bulk_tria, step_sizes, bl, tr,
colorize);
3431 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3433 if (internal::point_in_2d_box(cell->center(),
center, outer_radius))
3434 cells_to_remove.insert(cell);
3438 bulk_tria, cells_to_remove, tria_without_cylinder);
3440 const double tolerance =
3441 std::min(min_line_length(tria_without_cylinder),
3442 min_line_length(cylinder_tria)) /
3456 if (cell->manifold_id() == tfi_manifold_id)
3460 const auto &face = cell->face(face_n);
3461 if (face->at_boundary() &&
3462 internal::point_in_2d_box(face->center(),
3465 face->set_manifold_id(polar_manifold_id);
3467 face->set_manifold_id(tfi_manifold_id);
3478 static constexpr double tol =
3484 const auto face = cell->face(face_n);
3485 if (face->at_boundary())
3490 face->set_boundary_id(0);
3493 face->set_boundary_id(1);
3496 face->set_boundary_id(2);
3499 face->set_boundary_id(3);
3503 Assert(cell->manifold_id() == tfi_manifold_id,
3505 face->set_boundary_id(4);
3524 const double inner_radius,
3525 const double outer_radius,
3526 const double pad_bottom,
3527 const double pad_top,
3528 const double pad_left,
3529 const double pad_right,
3534 const unsigned int n_slices,
3545 Point<2>(new_center[0], new_center[1]),
3565 tria.
set_manifold(polar_manifold_id, cylindrical_manifold);
3573 const double shell_region_width,
3574 const unsigned int n_shells,
3575 const double skewness,
3578 Assert(0.0 <= shell_region_width && shell_region_width < 0.05,
3579 ExcMessage(
"The width of the shell region must be less than 0.05 "
3580 "(and preferably close to 0.03)"));
3614 std::set<Triangulation<2>::active_cell_iterator> cells_to_remove;
3618 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3619 cells_to_remove.insert(cell);
3623 for (
const unsigned int vertex_n :
3625 if (cell->vertex(vertex_n) ==
Point<2>())
3629 cylinder_triangulation_offset =
3630 2.0 * (cell->vertex(3) -
Point<2>());
3637 bulk_tria, cells_to_remove, tria_without_cylinder);
3644 0.05 + shell_region_width,
3652 if (
std::abs(cell->vertex(vertex_n)[0] - -0.41 / 4.0) < 1
e-10)
3653 cell->vertex(vertex_n)[0] = -0.1;
3654 else if (
std::abs(cell->vertex(vertex_n)[0] - 0.41 / 4.0) < 1
e-10)
3655 cell->vertex(vertex_n)[0] = 0.1;
3661 cell->set_manifold_id(tfi_manifold_id);
3663 if (!cell->face(face_n)->at_boundary())
3664 cell->face(face_n)->set_manifold_id(tfi_manifold_id);
3666 if (0.0 < shell_region_width)
3669 ExcMessage(
"If the shell region has positive width then "
3670 "there must be at least one shell."));
3675 0.05 + shell_region_width,
3682 const double vertex_tolerance =
3683 std::min(internal::minimal_vertex_distance(shell_tria),
3684 internal::minimal_vertex_distance(cylinder_tria)) *
3690 shell_tria, cylinder_tria, temp, vertex_tolerance,
true);
3691 cylinder_tria = std::move(temp);
3697 const double vertex_tolerance =
3698 std::min(internal::minimal_vertex_distance(tria_without_cylinder),
3699 internal::minimal_vertex_distance(cylinder_tria)) /
3702 tria_without_cylinder, cylinder_tria, tria, vertex_tolerance,
true);
3715 const double shift =
3716 std::min(0.125 + shell_region_width * 0.5, 0.1 * 4. / 3.);
3719 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3721 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1
e-10)
3723 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1
e-10)
3725 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1
e-10)
3732 if (cell->manifold_id() == polar_manifold_id)
3733 cell->set_all_manifold_ids(polar_manifold_id);
3738 if (cell->manifold_id() != polar_manifold_id &&
3739 cell->manifold_id() != tfi_manifold_id)
3744 std::vector<Point<2> *> cylinder_pointers;
3746 if (face->manifold_id() == polar_manifold_id)
3748 cylinder_pointers.push_back(&face->vertex(0));
3749 cylinder_pointers.push_back(&face->vertex(1));
3752 std::sort(cylinder_pointers.begin(), cylinder_pointers.end());
3753 cylinder_pointers.erase(std::unique(cylinder_pointers.begin(),
3754 cylinder_pointers.end()),
3755 cylinder_pointers.end());
3759 for (
const Point<2> *
const ptr : cylinder_pointers)
3760 center += *ptr / double(cylinder_pointers.size());
3763 for (
Point<2> *
const ptr : cylinder_pointers)
3775 if (face->at_boundary())
3780 face->set_boundary_id(0);
3783 face->set_boundary_id(1);
3785 else if (face->manifold_id() == polar_manifold_id)
3786 face->set_boundary_id(2);
3793 face->set_boundary_id(3);
3802 const double shell_region_width,
3803 const unsigned int n_shells,
3804 const double skewness,
3809 tria_2, shell_region_width, n_shells, skewness,
colorize);
3825 tria.
set_manifold(cylindrical_manifold_id, cylindrical_manifold);
3833 if (face->boundary_id() == 4 || face->boundary_id() == 5)
3834 face->set_boundary_id(3);
3839 template <
int dim,
int spacedim>
3842 const std::vector<unsigned int> &sizes,
3852 for (
unsigned int d = 0;
d < dim; ++
d)
3855 std::vector<Point<spacedim>> points;
3860 std::vector<CellData<dim>> cells(
n_cells);
3865 for (
unsigned int d = 0;
d < dim; ++
d)
3866 p(
d) = 0.5 * dimensions[
d] *
3869 points.push_back(p);
3870 cells[0].vertices[i] = i;
3872 cells[0].material_id = 0;
3883 for (
unsigned int j = 0; j < sizes[face]; ++j, ++
cell_index)
3885 const unsigned int last_cell = (j == 0) ? 0
U : (
cell_index - 1);
3887 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
3890 const unsigned int cellv =
3892 const unsigned int ocellv =
3896 cells[last_cell].vertices[cellv];
3899 cells[
cell_index].vertices[cellv] = points.size();
3904 points.push_back(p);
4041 const double thickness,
4045 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4047 std::vector<Point<2>>
vertices(16);
4049 coords[0] = left - thickness;
4052 coords[3] = right + thickness;
4055 for (
const double y : coords)
4056 for (
const double x : coords)
4061 std::vector<CellData<2>> cells(9);
4063 for (
unsigned int i0 = 0; i0 < 3; ++i0)
4064 for (
unsigned int i1 = 0; i1 < 3; ++i1)
4066 cells[k].vertices[0] = i1 + 4 * i0;
4067 cells[k].vertices[1] = i1 + 4 * i0 + 1;
4068 cells[k].vertices[2] = i1 + 4 * i0 + 4;
4069 cells[k].vertices[3] = i1 + 4 * i0 + 5;
4071 cells[k].material_id = materials[k];
4088 const double rl2 = (right + left) / 2;
4099 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4104 for (
unsigned int i = 0; i < 4; ++i)
4106 for (
unsigned int j = 0; j < 4; ++j)
4107 cells[i].
vertices[j] = cell_vertices[i][j];
4108 cells[i].material_id = 0;
4118 cell->face(1)->set_boundary_id(1);
4120 cell->face(0)->set_boundary_id(2);
4128 const double radius_0,
4129 const double radius_1,
4130 const double half_length)
4134 vertices_tmp[0] =
Point<2>(-half_length, -radius_0);
4135 vertices_tmp[1] =
Point<2>(half_length, -radius_1);
4136 vertices_tmp[2] =
Point<2>(-half_length, radius_0);
4137 vertices_tmp[3] =
Point<2>(half_length, radius_1);
4144 cell_vertices[0][i] = i;
4149 cells[0].vertices[i] = cell_vertices[0][i];
4151 cells[0].material_id = 0;
4156 cell->face(0)->set_boundary_id(1);
4157 cell->face(1)->set_boundary_id(2);
4159 for (
unsigned int i = 2; i < 4; ++i)
4160 cell->face(i)->set_boundary_id(0);
4180 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4184 for (
unsigned int i = 0; i < 3; ++i)
4186 for (
unsigned int j = 0; j < 4; ++j)
4187 cells[i].
vertices[j] = cell_vertices[i][j];
4188 cells[i].material_id = 0;
4200 cell->face(0)->set_boundary_id(0);
4201 cell->face(2)->set_boundary_id(1);
4204 cell->face(1)->set_boundary_id(2);
4205 cell->face(2)->set_boundary_id(1);
4206 cell->face(3)->set_boundary_id(3);
4209 cell->face(0)->set_boundary_id(0);
4210 cell->face(1)->set_boundary_id(4);
4211 cell->face(3)->set_boundary_id(5);
4217 template <
int dim,
int spacedim>
4220 const std::vector<unsigned int> &repetitions,
4223 const std::vector<int> & n_cells_to_remove)
4229 for (
unsigned int d = 0;
d < dim; ++
d)
4232 ExcMessage(
"Attempting to cut away too many cells."));
4242 std::array<double, dim> h;
4244 for (
unsigned int d = 0;
d < dim; ++
d)
4247 h[
d] = (top_right[
d] - bottom_left[
d]) / repetitions[
d];
4249 if (n_cells_to_remove[
d] >= 0)
4253 h[
d] *
std::fabs(n_cells_to_remove[
d]) + bottom_left[
d];
4258 cut_step[
d] = top_right[
d] - h[
d] *
std::fabs(n_cells_to_remove[
d]);
4264 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4269 std::inserter(cells_to_remove, cells_to_remove.end()),
4273 for (unsigned int d = 0; d < dim; ++d)
4274 if ((n_cells_to_remove[d] > 0 && cell->center()[d] >= cut_step[d]) ||
4275 (n_cells_to_remove[d] < 0 && cell->center()[d] <= cut_step[d]))
4292 const double radius,
4293 const bool internal_manifolds)
4298 const double a = 1. / (1 +
std::sqrt(2.0));
4309 const int cell_vertices[5][4] = {
4310 {0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}};
4314 for (
unsigned int i = 0; i < 5; ++i)
4316 for (
unsigned int j = 0; j < 4; ++j)
4317 cells[i].
vertices[j] = cell_vertices[i][j];
4318 cells[i].material_id = 0;
4328 if (internal_manifolds)
4337 const double inner_radius,
4338 const double outer_radius,
4342 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4356 const unsigned int N =
4357 (
n_cells == 0 ?
static_cast<unsigned int>(
4358 std::ceil((2 * pi * (outer_radius + inner_radius) / 2) /
4359 (outer_radius - inner_radius))) :
4369 for (
unsigned int i = 0; i <
N; ++i)
4382 for (
unsigned int i = 0; i <
N; ++i)
4384 cells[i].vertices[0] = i;
4385 cells[i].vertices[1] = (i + 1) %
N;
4386 cells[i].vertices[2] =
N + i;
4387 cells[i].vertices[3] =
N + ((i + 1) %
N);
4389 cells[i].material_id = 0;
4395 colorize_hyper_shell(tria,
center, inner_radius, outer_radius);
4408 const double inner_radius,
4409 const double outer_radius,
4413 tria, outer_center, inner_radius, outer_radius,
n_cells,
true);
4417 outer_radius - inner_radius > outer_center.
distance(inner_center),
4419 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4423 std::set<Point<dim> *> vertices_to_move;
4426 if (face->boundary_id() == 0)
4427 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
4428 vertices_to_move.insert(&face->vertex(v));
4430 const auto shift = inner_center - outer_center;
4431 for (
const auto &p : vertices_to_move)
4458 const double radius,
4459 const double half_length)
4461 Point<2> p1(-half_length, -radius);
4470 switch (f->boundary_id())
4473 f->set_boundary_id(1);
4476 f->set_boundary_id(2);
4479 f->set_boundary_id(0);
4513 const double radius)
4515 const unsigned int dim = 2;
4529 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4533 for (
unsigned int i = 0; i < 3; ++i)
4535 for (
unsigned int j = 0; j < 4; ++j)
4536 cells[i].
vertices[j] = cell_vertices[i][j];
4537 cells[i].material_id = 0;
4554 if (cell->face(i)->boundary_id() ==
4560 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius ||
4561 cell->face(i)->center()(1) < p(1) + 1.e-5 * radius)
4563 cell->face(i)->set_boundary_id(1);
4576 const double radius)
4581 const double a = 1. / (1 +
std::sqrt(2.0));
4592 const int cell_vertices[5][4] = {{0, 1, 2, 3},
4599 for (
unsigned int i = 0; i < 4; ++i)
4601 for (
unsigned int j = 0; j < 4; ++j)
4602 cells[i].
vertices[j] = cell_vertices[i][j];
4603 cells[i].material_id = 0;
4620 if (cell->face(i)->boundary_id() ==
4625 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius)
4627 cell->face(i)->set_boundary_id(1);
4642 const double inner_radius,
4643 const double outer_radius,
4647 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4660 const unsigned int N =
4661 (
n_cells == 0 ?
static_cast<unsigned int>(
4662 std::ceil((pi * (outer_radius + inner_radius) / 2) /
4663 (outer_radius - inner_radius))) :
4672 std::vector<Point<2>>
vertices(2 * (
N + 1));
4673 for (
unsigned int i = 0; i <=
N; ++i)
4694 for (
unsigned int i = 0; i <
N; ++i)
4696 cells[i].vertices[0] = i;
4697 cells[i].vertices[1] = (i + 1) % (
N + 1);
4698 cells[i].vertices[2] =
N + 1 + i;
4699 cells[i].vertices[3] =
N + 1 + ((i + 1) % (
N + 1));
4701 cells[i].material_id = 0;
4709 for (; cell != tria.
end(); ++cell)
4711 cell->face(2)->set_boundary_id(1);
4713 tria.
begin()->face(0)->set_boundary_id(3);
4715 tria.
last()->face(1)->set_boundary_id(2);
4725 const double inner_radius,
4726 const double outer_radius,
4730 Assert((inner_radius > 0) && (inner_radius < outer_radius),
4743 const unsigned int N =
4744 (
n_cells == 0 ?
static_cast<unsigned int>(
4745 std::ceil((pi * (outer_radius + inner_radius) / 4) /
4746 (outer_radius - inner_radius))) :
4755 std::vector<Point<2>>
vertices(2 * (
N + 1));
4756 for (
unsigned int i = 0; i <=
N; ++i)
4775 for (
unsigned int i = 0; i <
N; ++i)
4777 cells[i].vertices[0] = i;
4778 cells[i].vertices[1] = (i + 1) % (
N + 1);
4779 cells[i].vertices[2] =
N + 1 + i;
4780 cells[i].vertices[3] =
N + 1 + ((i + 1) % (
N + 1));
4782 cells[i].material_id = 0;
4790 for (; cell != tria.
end(); ++cell)
4792 cell->face(2)->set_boundary_id(1);
4794 tria.
begin()->face(0)->set_boundary_id(3);
4796 tria.
last()->face(1)->set_boundary_id(2);
4812 const double rl2 = (right + left) / 2;
4813 const double len = (right - left) / 2.;
4826 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4827 {9, 4, 2, 5, 19, 14, 12, 15},
4828 {3, 2, 7, 6, 13, 12, 17, 16},
4829 {2, 5, 6, 8, 12, 15, 16, 18}};
4831 for (
unsigned int i = 0; i < 4; ++i)
4833 for (
unsigned int j = 0; j < 8; ++j)
4834 cells[i].
vertices[j] = cell_vertices[i][j];
4835 cells[i].material_id = 0;
4845 cell->face(1)->set_boundary_id(1);
4847 cell->face(0)->set_boundary_id(2);
4858 const double thickness,
4862 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4864 std::vector<Point<3>>
vertices(64);
4866 coords[0] = left - thickness;
4869 coords[3] = right + thickness;
4872 for (
const double z : coords)
4873 for (
const double y : coords)
4874 for (
const double x : coords)
4878 24, 26, 5, 4, 6, 1, 0,
4879 2, 9, 8, 10, 37, 36, 38,
4880 33, 32, 34, 41, 40, 42};
4882 std::vector<CellData<3>> cells(27);
4884 for (
unsigned int z = 0; z < 3; ++z)
4885 for (
unsigned int y = 0; y < 3; ++y)
4886 for (
unsigned int x = 0; x < 3; ++x)
4888 cells[k].vertices[0] = x + 4 * y + 16 * z;
4889 cells[k].vertices[1] = x + 4 * y + 16 * z + 1;
4890 cells[k].vertices[2] = x + 4 * y + 16 * z + 4;
4891 cells[k].vertices[3] = x + 4 * y + 16 * z + 5;
4892 cells[k].vertices[4] = x + 4 * y + 16 * z + 16;
4893 cells[k].vertices[5] = x + 4 * y + 16 * z + 17;
4894 cells[k].vertices[6] = x + 4 * y + 16 * z + 20;
4895 cells[k].vertices[7] = x + 4 * y + 16 * z + 21;
4897 cells[k].material_id = materials[k];
4909 const double radius_0,
4910 const double radius_1,
4911 const double half_length)
4914 ExcMessage(
"The output triangulation object needs to be empty."));
4919 const auto n_slices = 1 +
static_cast<unsigned int>(
std::ceil(
4920 half_length /
std::max(radius_0, radius_1)));
4933 auto shift_radii = [=](
const Point<3> &p) {
4934 const double slope = (radius_1 / radius_0 - 1.0) / (2.0 * half_length);
4935 const double factor = slope * (p[0] - -half_length) + 1.0;
4936 return Point<3>(p[0], factor * p[1], factor * p[2]);
4942 for (
const auto &face :
triangulation.active_face_iterators())
4943 if (face->at_boundary())
4945 if (
std::abs(face->center()[0] - -half_length) < 1
e-8 * half_length)
4946 face->set_boundary_id(1);
4947 else if (
std::abs(face->center()[0] - half_length) <
4949 face->set_boundary_id(2);
4998 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
4999 {1, 2, 10, 11, 4, 5, 13, 14},
5000 {3, 4, 12, 13, 6, 7, 15, 16},
5001 {4, 5, 13, 14, 7, 8, 16, 17},
5002 {9, 10, 18, 19, 12, 13, 21, 22},
5003 {10, 11, 19, 20, 13, 14, 22, 23},
5004 {12, 13, 21, 22, 15, 16, 24, 25}};
5008 for (
unsigned int i = 0; i < 7; ++i)
5010 for (
unsigned int j = 0; j < 8; ++j)
5011 cells[i].
vertices[j] = cell_vertices[i][j];
5012 cells[i].material_id = 0;
5032 const double radius,
5033 const bool internal_manifold)
5039 const unsigned int n_vertices = 16;
5065 const unsigned int n_cells = 7;
5066 const int cell_vertices[
n_cells][8] = {
5067 {0, 1, 4, 5, 3, 2, 7, 6},
5068 {8, 9, 12, 13, 0, 1, 4, 5},
5069 {9, 13, 1, 5, 10, 14, 2, 6},
5070 {11, 10, 3, 2, 15, 14, 7, 6},
5071 {8, 0, 12, 4, 11, 3, 15, 7},
5072 {8, 9, 0, 1, 11, 10, 3, 2},
5073 {12, 4, 13, 5, 15, 7, 14, 6}};
5077 for (
unsigned int i = 0; i <
n_cells; ++i)
5080 cells[i].vertices[j] = cell_vertices[i][j];
5081 cells[i].material_id = 0;
5091 if (internal_manifold)
5097 const bool rotate_left_square,
5098 const bool rotate_right_square)
5100 constexpr unsigned int dim = 2;
5102 const unsigned int n_cells = 2;
5103 std::vector<CellData<dim>> cells(
n_cells);
5115 unsigned int cell_vertices[
n_cells][4] = {{0, 1, 2, 3},
5119 unsigned int this_case = 2 * rotate_left_square + rotate_right_square;
5125 cell_vertices[1][0] = 4;
5126 cell_vertices[1][1] = 5;
5127 cell_vertices[1][2] = 1;
5128 cell_vertices[1][3] = 3;
5134 cell_vertices[0][0] = 1;
5135 cell_vertices[0][1] = 3;
5136 cell_vertices[0][2] = 0;
5137 cell_vertices[0][3] = 2;
5143 cell_vertices[0][0] = 1;
5144 cell_vertices[0][1] = 3;
5145 cell_vertices[0][2] = 0;
5146 cell_vertices[0][3] = 2;
5148 cell_vertices[1][0] = 4;
5149 cell_vertices[1][1] = 5;
5150 cell_vertices[1][2] = 1;
5151 cell_vertices[1][3] = 3;
5163 for (
const unsigned int vertex_index :
5177 const bool face_orientation,
5178 const bool face_flip,
5179 const bool face_rotation,
5180 const bool manipulate_left_cube)
5182 constexpr unsigned int dim = 3;
5184 const unsigned int n_cells = 2;
5185 std::vector<CellData<dim>> cells(
n_cells);
5201 unsigned int cell_vertices[
n_cells][8] = {
5202 {0, 1, 2, 3, 4, 5, 6, 7},
5203 {1, 8, 3, 9, 5, 10, 7, 11}};
5206 const unsigned int this_case =
5207 4 * face_orientation + 2 * face_flip + face_rotation;
5209 if (manipulate_left_cube)
5215 cell_vertices[0][0] = 1;
5216 cell_vertices[0][1] = 0;
5217 cell_vertices[0][2] = 5;
5218 cell_vertices[0][3] = 4;
5219 cell_vertices[0][4] = 3;
5220 cell_vertices[0][5] = 2;
5221 cell_vertices[0][6] = 7;
5222 cell_vertices[0][7] = 6;
5228 cell_vertices[0][0] = 5;
5229 cell_vertices[0][1] = 4;
5230 cell_vertices[0][2] = 7;
5231 cell_vertices[0][3] = 6;
5232 cell_vertices[0][4] = 1;
5233 cell_vertices[0][5] = 0;
5234 cell_vertices[0][6] = 3;
5235 cell_vertices[0][7] = 2;
5241 cell_vertices[0][0] = 7;
5242 cell_vertices[0][1] = 6;
5243 cell_vertices[0][2] = 3;
5244 cell_vertices[0][3] = 2;
5245 cell_vertices[0][4] = 5;
5246 cell_vertices[0][5] = 4;
5247 cell_vertices[0][6] = 1;
5248 cell_vertices[0][7] = 0;
5253 cell_vertices[0][0] = 3;
5254 cell_vertices[0][1] = 2;
5255 cell_vertices[0][2] = 1;
5256 cell_vertices[0][3] = 0;
5257 cell_vertices[0][4] = 7;
5258 cell_vertices[0][5] = 6;
5259 cell_vertices[0][6] = 5;
5260 cell_vertices[0][7] = 4;
5266 cell_vertices[0][0] = 0;
5267 cell_vertices[0][1] = 1;
5268 cell_vertices[0][2] = 2;
5269 cell_vertices[0][3] = 3;
5270 cell_vertices[0][4] = 4;
5271 cell_vertices[0][5] = 5;
5272 cell_vertices[0][6] = 6;
5273 cell_vertices[0][7] = 7;
5279 cell_vertices[0][0] = 2;
5280 cell_vertices[0][1] = 3;
5281 cell_vertices[0][2] = 6;
5282 cell_vertices[0][3] = 7;
5283 cell_vertices[0][4] = 0;
5284 cell_vertices[0][5] = 1;
5285 cell_vertices[0][6] = 4;
5286 cell_vertices[0][7] = 5;
5292 cell_vertices[0][0] = 6;
5293 cell_vertices[0][1] = 7;
5294 cell_vertices[0][2] = 4;
5295 cell_vertices[0][3] = 5;
5296 cell_vertices[0][4] = 2;
5297 cell_vertices[0][5] = 3;
5298 cell_vertices[0][6] = 0;
5299 cell_vertices[0][7] = 1;
5305 cell_vertices[0][0] = 4;
5306 cell_vertices[0][1] = 5;
5307 cell_vertices[0][2] = 0;
5308 cell_vertices[0][3] = 1;
5309 cell_vertices[0][4] = 6;
5310 cell_vertices[0][5] = 7;
5311 cell_vertices[0][6] = 2;
5312 cell_vertices[0][7] = 3;
5323 cell_vertices[1][0] = 8;
5324 cell_vertices[1][1] = 1;
5325 cell_vertices[1][2] = 10;
5326 cell_vertices[1][3] = 5;
5327 cell_vertices[1][4] = 9;
5328 cell_vertices[1][5] = 3;
5329 cell_vertices[1][6] = 11;
5330 cell_vertices[1][7] = 7;
5336 cell_vertices[1][0] = 10;
5337 cell_vertices[1][1] = 5;
5338 cell_vertices[1][2] = 11;
5339 cell_vertices[1][3] = 7;
5340 cell_vertices[1][4] = 8;
5341 cell_vertices[1][5] = 1;
5342 cell_vertices[1][6] = 9;
5343 cell_vertices[1][7] = 3;
5349 cell_vertices[1][0] = 11;
5350 cell_vertices[1][1] = 7;
5351 cell_vertices[1][2] = 9;
5352 cell_vertices[1][3] = 3;
5353 cell_vertices[1][4] = 10;
5354 cell_vertices[1][5] = 5;
5355 cell_vertices[1][6] = 8;
5356 cell_vertices[1][7] = 1;
5362 cell_vertices[1][0] = 9;
5363 cell_vertices[1][1] = 3;
5364 cell_vertices[1][2] = 8;
5365 cell_vertices[1][3] = 1;
5366 cell_vertices[1][4] = 11;
5367 cell_vertices[1][5] = 7;
5368 cell_vertices[1][6] = 10;
5369 cell_vertices[1][7] = 5;
5375 cell_vertices[1][0] = 1;
5376 cell_vertices[1][1] = 8;
5377 cell_vertices[1][2] = 3;
5378 cell_vertices[1][3] = 9;
5379 cell_vertices[1][4] = 5;
5380 cell_vertices[1][5] = 10;
5381 cell_vertices[1][6] = 7;
5382 cell_vertices[1][7] = 11;
5388 cell_vertices[1][0] = 5;
5389 cell_vertices[1][1] = 10;
5390 cell_vertices[1][2] = 1;
5391 cell_vertices[1][3] = 8;
5392 cell_vertices[1][4] = 7;
5393 cell_vertices[1][5] = 11;
5394 cell_vertices[1][6] = 3;
5395 cell_vertices[1][7] = 9;
5401 cell_vertices[1][0] = 7;
5402 cell_vertices[1][1] = 11;
5403 cell_vertices[1][2] = 5;
5404 cell_vertices[1][3] = 10;
5405 cell_vertices[1][4] = 3;
5406 cell_vertices[1][5] = 9;
5407 cell_vertices[1][6] = 1;
5408 cell_vertices[1][7] = 8;
5414 cell_vertices[1][0] = 3;
5415 cell_vertices[1][1] = 9;
5416 cell_vertices[1][2] = 7;
5417 cell_vertices[1][3] = 11;
5418 cell_vertices[1][4] = 1;
5419 cell_vertices[1][5] = 8;
5420 cell_vertices[1][6] = 5;
5421 cell_vertices[1][7] = 10;
5431 for (
const unsigned int vertex_index :
5445 template <
int spacedim>
5448 const double radius)
5452 std::set<types::boundary_id> boundary_ids;
5453 boundary_ids.insert(0);
5464 const unsigned int x_subdivisions,
5465 const double radius,
5466 const double half_length)
5474 const double initial_height = -half_length;
5475 const double height_increment = 2. * half_length / x_subdivisions;
5477 for (
unsigned int rep = 0; rep < (x_subdivisions + 1); ++rep)
5479 const double height = initial_height + height_increment * rep;
5494 const double h = vertex(1);
5495 vertex(1) = -vertex(0);
5499 std::vector<std::vector<int>> cell_vertices;
5500 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5501 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5502 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5503 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5504 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5506 for (
unsigned int rep = 1; rep < x_subdivisions; ++rep)
5508 for (
unsigned int i = 0; i < 5; ++i)
5510 std::vector<int> new_cell_vertices(8);
5511 for (
unsigned int j = 0; j < 8; ++j)
5512 new_cell_vertices[j] = cell_vertices[i][j] + 8 * rep;
5513 cell_vertices.push_back(new_cell_vertices);
5517 unsigned int n_cells = x_subdivisions * 5;
5521 for (
unsigned int i = 0; i <
n_cells; ++i)
5523 for (
unsigned int j = 0; j < 8; ++j)
5524 cells[i].
vertices[j] = cell_vertices[i][j];
5525 cells[i].material_id = 0;
5547 if (cell->at_boundary(i))
5549 if (cell->face(i)->center()(0) > half_length - 1.e-5)
5551 cell->face(i)->set_boundary_id(2);
5554 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5556 if ((
std::fabs(cell->face(i)->line(
e)->vertex(0)[1]) == a) ||
5557 (
std::fabs(cell->face(i)->line(
e)->vertex(0)[2]) == a) ||
5558 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[1]) == a) ||
5559 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[2]) == a))
5561 cell->face(i)->line(
e)->set_boundary_id(2);
5562 cell->face(i)->line(
e)->set_manifold_id(
5566 else if (cell->face(i)->center()(0) < -half_length + 1.e-5)
5568 cell->face(i)->set_boundary_id(1);
5571 for (
unsigned int e = 0; e < GeometryInfo<3>::lines_per_face;
5573 if ((
std::fabs(cell->face(i)->line(
e)->vertex(0)[1]) == a) ||
5574 (
std::fabs(cell->face(i)->line(
e)->vertex(0)[2]) == a) ||
5575 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[1]) == a) ||
5576 (
std::fabs(cell->face(i)->line(
e)->vertex(1)[2]) == a))
5578 cell->face(i)->line(
e)->set_boundary_id(1);
5579 cell->face(i)->line(
e)->set_manifold_id(
5590 const double radius,
5591 const double half_length)
5599 const double radius)
5601 const unsigned int dim = 3;
5609 const double a = 0.528;
5610 const double b = 0.4533;
5611 const double c = 0.3752;
5628 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5629 {1, 6, 2, 4, 8, 13, 9, 11},
5630 {5, 3, 6, 4, 12, 10, 13, 11},
5631 {7, 9, 10, 11, 14, 8, 12, 13}};
5635 for (
unsigned int i = 0; i < 4; ++i)
5637 for (
unsigned int j = 0; j < 8; ++j)
5638 cells[i].
vertices[j] = cell_vertices[i][j];
5639 cells[i].material_id = 0;
5655 if (cell->face(i)->boundary_id() ==
5660 if (cell->face(i)->center()(0) <
center(0) + 1.e-5 * radius ||
5661 cell->face(i)->center()(1) <
center(1) + 1.e-5 * radius ||
5662 cell->face(i)->center()(2) <
center(2) + 1.e-5 * radius)
5664 cell->face(i)->set_boundary_id(1);
5668 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5671 const Point<3> line_vertices[2] = {
5672 cell->face(i)->line(j)->vertex(0),
5673 cell->face(i)->line(j)->vertex(1)};
5679 cell->face(i)->line(j)->set_boundary_id(1);
5680 cell->face(i)->line(j)->set_manifold_id(
5697 const double radius)
5703 const double b = a / 2.0;
5704 const double c =
d / 2.0;
5706 const double hb = radius *
std::sqrt(3.0) / 4.0;
5707 const double hc = radius *
std::sqrt(3.0) / 2.0;
5729 int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5730 {0, 2, 8, 10, 6, 4, 14, 12},
5731 {2, 3, 10, 11, 4, 5, 12, 13},
5732 {1, 7, 9, 15, 3, 5, 11, 13},
5733 {6, 4, 14, 12, 7, 5, 15, 13},
5734 {8, 10, 9, 11, 14, 12, 15, 13}};
5738 for (
unsigned int i = 0; i < 6; ++i)
5740 for (
unsigned int j = 0; j < 8; ++j)
5741 cells[i].
vertices[j] = cell_vertices[i][j];
5742 cells[i].material_id = 0;
5764 if (!cell->at_boundary(i))
5770 if (cell->face(i)->center()(0) <
center(0) + 1.e-5 * radius)
5772 cell->face(i)->set_boundary_id(1);
5774 for (
unsigned int j = 0; j < GeometryInfo<3>::lines_per_face;
5777 const Point<3> line_vertices[2] = {
5778 cell->face(i)->line(j)->vertex(0),
5779 cell->face(i)->line(j)->vertex(1)};
5785 cell->face(i)->line(j)->set_boundary_id(1);
5786 cell->face(i)->line(j)->set_manifold_id(
5803 const double radius)
5813 for (
unsigned int round = 0; round < dim; ++round)
5818 std::vector<Point<dim>> new_points(tria_copy.
n_vertices());
5820 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5828 else if (round == 1)
5830 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5839 else if (round == 2)
5840 for (
unsigned int v = 0; v < tria_copy.
n_vertices(); ++v)
5853 std::vector<CellData<dim>> cells;
5854 cells.reserve(tria_copy.
n_cells());
5859 data.
vertices[v] = cell->vertex_index(v);
5862 cells.push_back(data);
5870 if (round == dim - 1)
5880 if (cell->center().norm_square() > 0.4 * radius)
5881 cell->set_manifold_id(1);
5894 const double inner_radius,
5895 const double outer_radius,
5899 Assert((inner_radius > 0) && (inner_radius < outer_radius),
5902 unsigned int n_refinement_steps = 0;
5903 unsigned int n_cells_coarsened =
n_cells;
5905 while (n_cells_coarsened > 12 && n_cells_coarsened % 4 == 0)
5907 ++n_refinement_steps;
5908 n_cells_coarsened /= 4;
5911 (n_refinement_steps > 0 &&
5912 (n_cells_coarsened == 6 || n_cells_coarsened == 12)),
5913 ExcMessage(
"Invalid number of coarse mesh cells"));
5915 const unsigned int n = n_refinement_steps > 0 ?
5916 4 * n_cells_coarsened :
5919 const double irad = inner_radius /
std::sqrt(3.0);
5920 const double orad = outer_radius /
std::sqrt(3.0);
5922 std::vector<CellData<3>> cells;
5925 static const std::array<Point<3>, 8> hexahedron = {{{-1, -1, -1},
5939 for (
unsigned int i = 0; i < 8; ++i)
5940 vertices.push_back(p + hexahedron[i] * irad);
5941 for (
unsigned int i = 0; i < 8; ++i)
5942 vertices.push_back(p + hexahedron[i] * orad);
5944 const unsigned int n_cells = 6;
5945 const int cell_vertices[
n_cells][8] = {
5946 {8, 9, 10, 11, 0, 1, 2, 3},
5947 {9, 11, 1, 3, 13, 15, 5, 7},
5948 {12, 13, 4, 5, 14, 15, 6, 7},
5949 {8, 0, 10, 2, 12, 4, 14, 6},
5950 {8, 9, 0, 1, 12, 13, 4, 5},
5951 {10, 2, 11, 3, 14, 6, 15, 7}};
5955 for (
unsigned int i = 0; i <
n_cells; ++i)
5958 cells[i].vertices[j] = cell_vertices[i][j];
5959 cells[i].material_id = 0;
5971 static const std::array<Point<3>, 6> octahedron = {{{-1, 0, 0},
5978 for (
unsigned int i = 0; i < 8; ++i)
5979 vertices.push_back(p + hexahedron[i] * irad);
5980 for (
unsigned int i = 0; i < 6; ++i)
5981 vertices.push_back(p + octahedron[i] * inner_radius);
5982 for (
unsigned int i = 0; i < 8; ++i)
5983 vertices.push_back(p + hexahedron[i] * orad);
5984 for (
unsigned int i = 0; i < 6; ++i)
5985 vertices.push_back(p + octahedron[i] * outer_radius);
5987 const unsigned int n_cells = 12;
5988 const unsigned int rhombi[
n_cells][4] = {{10, 4, 0, 8},
6003 for (
unsigned int i = 0; i <
n_cells; ++i)
6005 for (
unsigned int j = 0; j < 4; ++j)
6007 cells[i].vertices[j] = rhombi[i][j];
6008 cells[i].vertices[j + 4] = rhombi[i][j] + 14;
6010 cells[i].material_id = 0;
6031 const unsigned int outer_radius_factor = 1 << n_refinement_steps;
6035 outer_radius_factor * outer_radius -
6036 (outer_radius_factor - 1) * inner_radius,
6038 for (
unsigned int r = 0; r < n_refinement_steps; ++r)
6041 std::set<Triangulation<3>::active_cell_iterator>
6048 unsigned int n_vertices_inside = 0;
6050 if ((cell->vertex(v) - p).norm_square() <
6051 inner_radius * inner_radius * (1 + 1
e-12))
6052 ++n_vertices_inside;
6053 if (n_vertices_inside < 4)
6054 cells_to_remove.insert(cell);
6059 if (r == n_refinement_steps - 1)
6069 tmp = std::move(
copy);
6083 hyper_shell(tmp, p, inner_radius, outer_radius, 12);
6098 colorize_hyper_shell(tria, p, inner_radius, outer_radius);
6109 const double inner_radius,
6110 const double outer_radius,
6111 const unsigned int ,
6114 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6118 const double d = outer_radius /
std::sqrt(2.0);
6119 const double a = inner_radius /
std::sqrt(2.0);
6121 const double b = a / 2.0;
6122 const double c =
d / 2.0;
6124 const double hb = inner_radius *
std::sqrt(3.0) / 2.0;
6125 const double hc = outer_radius *
std::sqrt(3.0) / 2.0;
6147 int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
6148 {0, 2, 8, 10, 6, 4, 14, 12},
6149 {1, 7, 9, 15, 3, 5, 11, 13},
6150 {6, 4, 14, 12, 7, 5, 15, 13},
6151 {8, 10, 9, 11, 14, 12, 15, 13}};
6155 for (
unsigned int i = 0; i < 5; ++i)
6157 for (
unsigned int j = 0; j < 8; ++j)
6158 cells[i].
vertices[j] = cell_vertices[i][j];
6159 cells[i].material_id = 0;
6173 for (; cell != tria.
end(); ++cell)
6175 if (cell->at_boundary(i))
6176 cell->face(i)->set_all_boundary_ids(2);
6182 for (cell = tria.
begin(); cell != tria.
end(); ++cell)
6184 if (cell->at_boundary(i))
6188 const Point<3> face_center(face->center());
6190 1.e-6 * face_center.
norm())
6194 face->set_all_boundary_ids(0);
6196 face->set_all_boundary_ids(1);
6209 const double inner_radius,
6210 const double outer_radius,
6211 const unsigned int n,
6214 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6216 if (n == 0 || n == 3)
6218 const double a = inner_radius *
std::sqrt(2.0) / 2e0;
6219 const double b = outer_radius *
std::sqrt(2.0) / 2e0;
6220 const double c = a *
std::sqrt(3.0) / 2e0;
6222 const double e = outer_radius / 2e0;
6223 const double h = inner_radius / 2e0;
6242 const int cell_vertices[3][8] = {
6243 {0, 1, 3, 2, 4, 5, 7, 6},
6244 {1, 8, 2, 9, 5, 10, 6, 11},
6245 {4, 5, 7, 6, 12, 10, 13, 11},
6247 std::vector<CellData<3>> cells(3);
6249 for (
unsigned int i = 0; i < 3; ++i)
6251 for (
unsigned int j = 0; j < 8; ++j)
6252 cells[i].
vertices[j] = cell_vertices[i][j];
6253 cells[i].material_id = 0;
6266 colorize_quarter_hyper_shell(tria,
center, inner_radius, outer_radius);
6276 const double length,
6277 const double inner_radius,
6278 const double outer_radius,
6279 const unsigned int n_radial_cells,
6280 const unsigned int n_axial_cells)
6282 Assert((inner_radius > 0) && (inner_radius < outer_radius),
6296 const unsigned int N_r =
6297 (n_radial_cells == 0 ?
static_cast<unsigned int>(
std::ceil(
6298 (2 * pi * (outer_radius + inner_radius) / 2) /
6299 (outer_radius - inner_radius))) :
6301 const unsigned int N_z =
6302 (n_axial_cells == 0 ?
6304 length / (2 * pi * (outer_radius + inner_radius) / 2 / N_r))) :
6313 std::vector<Point<2>> vertices_2d(2 * N_r);
6314 for (
unsigned int i = 0; i < N_r; ++i)
6319 vertices_2d[i + N_r] = vertices_2d[i] * (inner_radius / outer_radius);
6322 std::vector<Point<3>> vertices_3d;
6323 vertices_3d.reserve(2 * N_r * (N_z + 1));
6324 for (
unsigned int j = 0; j <= N_z; ++j)
6325 for (
unsigned int i = 0; i < 2 * N_r; ++i)
6327 const Point<3> v(vertices_2d[i][0],
6330 vertices_3d.push_back(v);
6333 std::vector<CellData<3>> cells(N_r * N_z,
CellData<3>());
6335 for (
unsigned int j = 0; j < N_z; ++j)
6336 for (
unsigned int i = 0; i < N_r; ++i)
6338 cells[i + j * N_r].vertices[0] = i + (j + 1) * 2 * N_r;
6339 cells[i + j * N_r].vertices[1] = (i + 1) % N_r + (j + 1) * 2 * N_r;
6340 cells[i + j * N_r].vertices[2] = i + j * 2 * N_r;
6341 cells[i + j * N_r].vertices[3] = (i + 1) % N_r + j * 2 * N_r;
6343 cells[i + j * N_r].vertices[4] = N_r + i + (j + 1) * 2 * N_r;
6344 cells[i + j * N_r].vertices[5] =
6345 N_r + ((i + 1) % N_r) + (j + 1) * 2 * N_r;
6346 cells[i + j * N_r].vertices[6] = N_r + i + j * 2 * N_r;
6347 cells[i + j * N_r].vertices[7] = N_r + ((i + 1) % N_r) + j * 2 * N_r;
6349 cells[i + j * N_r].material_id = 0;
6359 template <
int dim,
int spacedim>
6364 const double duplicated_vertex_tolerance,
6365 const bool copy_manifold_ids)
6367 std::vector<Point<spacedim>>
vertices;
6368 std::vector<CellData<dim>> cells;
6371 unsigned int n_accumulated_vertices = 0;
6375 ExcMessage(
"The input triangulations must be non-empty "
6376 "and must not be refined."));
6378 std::vector<Point<spacedim>> tria_vertices;
6379 std::vector<CellData<dim>> tria_cells;
6381 std::tie(tria_vertices, tria_cells, tria_subcell_data) =
6385 tria_vertices.begin(),
6386 tria_vertices.end());
6389 for (
unsigned int &vertex_n : cell_data.vertices)
6390 vertex_n += n_accumulated_vertices;
6391 cells.push_back(cell_data);
6395 if (copy_manifold_ids)
6401 for (
unsigned int &vertex_n : line_data.
vertices)
6402 vertex_n += n_accumulated_vertices;
6412 for (
unsigned int &vertex_n : quad_data.
vertices)
6413 vertex_n += n_accumulated_vertices;
6424 std::vector<unsigned int> considered_vertices;
6428 considered_vertices,
6429 duplicated_vertex_tolerance);
6440 template <
int dim,
int spacedim>
6445 const double duplicated_vertex_tolerance,
6446 const bool copy_manifold_ids)
6449 if (triangulation_1.
n_cells() == 0)
6454 if (triangulation_2.
n_cells() == 0)
6461 duplicated_vertex_tolerance,
6490 template <
int structdim>
6494 static_assert(structdim == 1 || structdim == 2,
6495 "This function is only implemented for lines and "
6505 else if (structdim == 2)
6508 std::array<unsigned int, 4> renumbering;
6511 renumbering.begin());
6522 std::swap(renumbering[2], renumbering[3]);
6524 std::min_element(renumbering.begin(),
6528 std::swap(renumbering[2], renumbering[3]);
6536 if (renumbering[1] > renumbering[2])
6537 std::swap(renumbering[1], renumbering[2]);
6552 std::sort(subcell_data.begin(), subcell_data.end(), compare);
6557 auto left = subcell_data.begin();
6558 while (left != subcell_data.end())
6561 std::upper_bound(left, subcell_data.end(), *left, compare);
6564 if (left + 1 != right)
6565 for (
auto it = left; it != right; ++it)
6568 Assert(it->manifold_id == left->manifold_id,
6570 "In the process of grid generation a single "
6571 "line or quadrilateral has been assigned two "
6572 "different manifold ids. This can happen when "
6573 "a Triangulation is copied, e.g., via "
6574 "GridGenerator::replicate_triangulation() and "
6575 "not all external boundary faces have the same "
6576 "manifold id. Double check that all faces "
6577 "which you expect to be merged together have "
6578 "the same manifold id."));
6583 subcell_data.erase(std::unique(subcell_data.begin(), subcell_data.end()),
6584 subcell_data.end());
6590 template <
int dim,
int spacedim>
6593 const std::vector<unsigned int> & extents,
6598 for (
const auto &extent : extents)
6600 ExcMessage(
"The Triangulation must be copied at least one time in "
6601 "each coordinate dimension."));
6604 const auto &
min = bbox.get_boundary_points().first;
6605 const auto &
max = bbox.get_boundary_points().second;
6607 std::array<Tensor<1, spacedim>, dim> offsets;
6608 for (
unsigned int d = 0;
d < dim; ++
d)
6613 for (
unsigned int d = 0;
d < dim; ++
d)
6615 std::vector<Point<spacedim>> input_vertices;
6616 std::vector<CellData<dim>> input_cell_data;
6618 std::tie(input_vertices, input_cell_data, input_subcell_data) =
6620 std::vector<Point<spacedim>> output_vertices = input_vertices;
6621 std::vector<CellData<dim>> output_cell_data = input_cell_data;
6622 SubCellData output_subcell_data = input_subcell_data;
6624 for (
unsigned int k = 1; k < extents[
d]; ++k)
6626 const std::size_t vertex_offset = k * input_vertices.size();
6629 output_vertices.push_back(
point +
double(k) * offsets[
d]);
6633 output_cell_data.push_back(cell_data);
6634 for (
unsigned int &vertex : output_cell_data.back().vertices)
6635 vertex += vertex_offset;
6642 for (
unsigned int &vertex :
6644 vertex += vertex_offset;
6650 for (
unsigned int &vertex :
6652 vertex += vertex_offset;
6657 std::vector<unsigned int> boundary_vertices;
6661 output_subcell_data,
6682 tria_to_replicate.
clear();
6685 output_subcell_data);
6693 template <
int dim,
int spacedim>
6701 ExcMessage(
"The two input triangulations are not derived from "
6702 "the same coarse mesh as required."));
6705 &triangulation_1) ==
nullptr) &&
6708 &triangulation_2) ==
nullptr),
6709 ExcMessage(
"The source triangulations for this function must both "
6710 "be available entirely locally, and not be distributed "
6711 "triangulations."));
6728 for (
unsigned int iteration = 0; iteration < triangulation_2.
n_levels();
6734 bool any_cell_flagged =
false;
6736 if (intergrid_map[result_cell]->has_children())
6738 any_cell_flagged =
true;
6739 result_cell->set_refine_flag();
6742 if (any_cell_flagged ==
false)
6751 template <
int dim,
int spacedim>
6765 std::vector<CellData<dim>> cells;
6767 if (cells_to_remove.find(cell) == cells_to_remove.end())
6769 Assert(
static_cast<unsigned int>(cell->level()) ==
6770 input_triangulation.
n_levels() - 1,
6772 "Your input triangulation appears to have "
6773 "adaptively refined cells. This is not allowed. You can "
6774 "only call this function on a triangulation in which "
6775 "all cells are on the same refinement level."));
6779 this_cell.
vertices[v] = cell->vertex_index(v);
6781 cells.push_back(this_cell);
6787 std::vector<unsigned int> considered_vertices;
6791 considered_vertices);
6803 const unsigned int n_slices,
6804 const double height,
6806 const bool copy_manifold_ids,
6807 const std::vector<types::manifold_id> &manifold_priorities)
6811 "The input triangulation must be a coarse mesh, i.e., it must "
6812 "not have been refined."));
6814 ExcMessage(
"The output triangulation object needs to be empty."));
6816 ExcMessage(
"The given height for extrusion must be positive."));
6819 "The number of slices for extrusion must be at least 2."));
6821 const double delta_h = height / (n_slices - 1);
6822 std::vector<double> slices_z_values;
6823 for (
unsigned int i = 0; i < n_slices; ++i)
6824 slices_z_values.push_back(i * delta_h);
6826 input, slices_z_values, result, copy_manifold_ids, manifold_priorities);
6834 const unsigned int n_slices,
6835 const double height,
6837 const bool copy_manifold_ids,
6838 const std::vector<types::manifold_id> &manifold_priorities)
6844 (void)copy_manifold_ids;
6845 (void)manifold_priorities;
6849 "GridTools::extrude_triangulation() is only available "
6850 "for Triangulation<3, 3> as output triangulation."));
6858 const std::vector<double> & slice_coordinates,
6860 const bool copy_manifold_ids,
6861 const std::vector<types::manifold_id> &manifold_priorities)
6865 "The input triangulation must be a coarse mesh, i.e., it must "
6866 "not have been refined."));
6868 ExcMessage(
"The output triangulation object needs to be empty."));
6869 Assert(slice_coordinates.size() >= 2,
6871 "The number of slices for extrusion must be at least 2."));
6872 Assert(std::is_sorted(slice_coordinates.begin(), slice_coordinates.end()),
6873 ExcMessage(
"Slice z-coordinates should be in ascending order"));
6875 const auto priorities = [&]() -> std::vector<types::manifold_id> {
6879 if (0 < manifold_priorities.size())
6883 std::vector<types::manifold_id> sorted_manifold_priorities =
6884 manifold_priorities;
6885 std::sort(sorted_manifold_priorities.begin(),
6886 sorted_manifold_priorities.end());
6887 Assert(std::unique(sorted_manifold_priorities.begin(),
6888 sorted_manifold_priorities.end()) ==
6889 sorted_manifold_priorities.end(),
6891 "The given vector of manifold ids may not contain any "
6892 "duplicated entries."));
6893 std::vector<types::manifold_id> sorted_manifold_ids =
6895 std::sort(sorted_manifold_ids.begin(), sorted_manifold_ids.end());
6896 if (sorted_manifold_priorities != sorted_manifold_ids)
6898 std::ostringstream message;
6899 message <<
"The given triangulation has manifold ids {";
6903 message << sorted_manifold_ids.back() <<
"}, but \n"
6904 <<
" the given vector of manifold ids is {";
6909 << manifold_priorities.back() <<
"}.\n"
6910 <<
" These vectors should contain the same elements.\n";
6911 const std::string m = message.str();
6915 return manifold_priorities;
6919 std::vector<types::manifold_id> default_priorities =
6922 default_priorities.begin(),
6923 default_priorities.end(),
6925 return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
6926 &input.get_manifold(id)) == nullptr;
6928 std::sort(default_priorities.begin(), first_tfi_it);
6929 std::sort(first_tfi_it, default_priorities.end());
6931 return default_priorities;
6934 const std::size_t n_slices = slice_coordinates.size();
6935 std::vector<Point<3>> points(n_slices * input.
n_vertices());
6936 std::vector<CellData<3>> cells;
6941 for (std::size_t slice_n = 0; slice_n < n_slices; ++slice_n)
6943 for (std::size_t vertex_n = 0; vertex_n < input.
n_vertices();
6947 points[slice_n * input.
n_vertices() + vertex_n] =
6948 Point<3>(vertex[0], vertex[1], slice_coordinates[slice_n]);
6956 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
6959 for (
const unsigned int vertex_n :
6963 cell->vertex_index(vertex_n) + slice_n * input.
n_vertices();
6966 cell->vertex_index(vertex_n) +
6971 if (copy_manifold_ids)
6973 cells.push_back(this_cell);
6988 if (face->at_boundary())
6990 if (copy_manifold_ids)
6992 for (std::size_t slice_n = 0; slice_n < n_slices - 1; ++slice_n)
6995 face->vertex_index(0) + slice_n * input.
n_vertices();
6997 face->vertex_index(1) + slice_n * input.
n_vertices();
6999 face->vertex_index(0) + (slice_n + 1) * input.
n_vertices();
7001 face->vertex_index(1) + (slice_n + 1) * input.
n_vertices();
7002 quads.push_back(quad);
7008 if (copy_manifold_ids)
7014 for (std::size_t slice_n = 1; slice_n < n_slices - 1; ++slice_n)
7017 cell->vertex_index(0) + slice_n * input.
n_vertices();
7019 cell->vertex_index(1) + slice_n * input.
n_vertices();
7021 cell->vertex_index(2) + slice_n * input.
n_vertices();
7023 cell->vertex_index(3) + slice_n * input.
n_vertices();
7024 quads.push_back(quad);
7035 "The input triangulation to this function is using boundary "
7036 "indicators in a range that do not allow using "
7037 "max_boundary_id+1 and max_boundary_id+2 as boundary "
7038 "indicators for the bottom and top faces of the "
7039 "extruded triangulation."));
7046 quad.
vertices[0] = cell->vertex_index(0);
7047 quad.
vertices[1] = cell->vertex_index(1);
7048 quad.
vertices[2] = cell->vertex_index(2);
7049 quad.
vertices[3] = cell->vertex_index(3);
7050 if (copy_manifold_ids)
7052 quads.push_back(quad);
7055 for (
unsigned int &vertex : quad.
vertices)
7056 vertex += (n_slices - 1) * input.
n_vertices();
7057 if (copy_manifold_ids)
7059 quads.push_back(quad);
7072 for (
auto manifold_id_it = priorities.rbegin();
7073 manifold_id_it != priorities.rend();
7076 if (face->manifold_id() == *manifold_id_it)
7077 for (
unsigned int line_n = 0;
7078 line_n < GeometryInfo<3>::lines_per_face;
7080 face->line(line_n)->set_manifold_id(*manifold_id_it);
7088 const std::vector<double> & slice_coordinates,
7090 const bool copy_manifold_ids,
7091 const std::vector<types::manifold_id> &manifold_priorities)
7094 (void)slice_coordinates;
7096 (void)copy_manifold_ids;
7097 (void)manifold_priorities;
7101 "GridTools::extrude_triangulation() is only available "
7102 "for Triangulation<3, 3> as output triangulation."));
7122 const double inner_radius,
7123 const double outer_radius,
7130 Assert(inner_radius < outer_radius,
7131 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7140 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
7141 for (; cell != endc; ++cell)
7144 if (cell->face(f)->at_boundary())
7146 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
7149 unsigned int vv = cell->face(f)->vertex_index(v);
7150 if (treated_vertices[vv] ==
false)
7152 treated_vertices[vv] =
true;
7156 cell->face(f)->vertex(v) =
7160 cell->face(f)->vertex(v) =
7164 cell->face(f)->vertex(v) =
7168 cell->face(f)->vertex(v) =
7178 double eps = 1
e-3 * outer_radius;
7180 for (; cell != endc; ++cell)
7183 if (cell->face(f)->at_boundary())
7185 double dx = cell->face(f)->center()(0) -
center(0);
7186 double dy = cell->face(f)->center()(1) -
center(1);
7190 cell->face(f)->set_boundary_id(0);
7192 cell->face(f)->set_boundary_id(1);
7194 cell->face(f)->set_boundary_id(2);
7196 cell->face(f)->set_boundary_id(3);
7199 cell->face(f)->set_boundary_id(4);
7200 cell->face(f)->set_manifold_id(0);
7205 double d = (cell->face(f)->center() -
center).
norm();
7206 if (
d - inner_radius < 0)
7208 cell->face(f)->set_boundary_id(1);
7209 cell->face(f)->set_manifold_id(0);
7212 cell->face(f)->set_boundary_id(0);
7225 const double inner_radius,
7226 const double outer_radius,
7227 const unsigned int n_shells,
7228 const double skewness,
7235 Assert(inner_radius < outer_radius,
7236 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7240 std::vector<double> radii;
7241 radii.push_back(inner_radius);
7242 for (
unsigned int shell_n = 1; shell_n < n_shells; ++shell_n)
7243 if (skewness == 0.0)
7245 radii.push_back(inner_radius +
7246 (outer_radius - inner_radius) *
7247 (1.0 - (1.0 -
double(shell_n) / n_shells)));
7251 (outer_radius - inner_radius) *
7252 (1.0 -
std::tanh(skewness * (1.0 -
double(shell_n) / n_shells)) /
7254 radii.push_back(outer_radius);
7256 double grid_vertex_tolerance = 0.0;
7257 for (
unsigned int shell_n = 0; shell_n < radii.size() - 1; ++shell_n)
7264 n_cells == 0 ? (dim == 2 ? 8 : 12) :
7269 if (grid_vertex_tolerance == 0.0)
7270 grid_vertex_tolerance =
7271 0.5 * internal::minimal_vertex_distance(current_shell);
7278 grid_vertex_tolerance);
7290 constexpr double radial_vertex_tolerance =
7292 auto assert_vertex_distance_within_tolerance =
7293 [
center, radial_vertex_tolerance](
7295 const double radius) {
7297 (void)radial_vertex_tolerance;
7300 for (
unsigned int vertex_n = 0;
7301 vertex_n < GeometryInfo<dim>::vertices_per_face;
7305 (
center.
norm() + radius) * radial_vertex_tolerance,
7310 for (
const auto &cell :
triangulation.active_cell_iterators())
7313 auto face = cell->face(face_n);
7314 if (face->at_boundary())
7316 if (((face->vertex(0) -
center).norm() - inner_radius) <
7317 (
center.
norm() + inner_radius) * radial_vertex_tolerance)
7320 assert_vertex_distance_within_tolerance(face, inner_radius);
7321 face->set_all_boundary_ids(0);
7326 assert_vertex_distance_within_tolerance(face, outer_radius);
7327 face->set_all_boundary_ids(1);
7337 const double inner_radius,
7338 const double outer_radius,
7340 const unsigned int Nz,
7345 Assert(inner_radius < outer_radius,
7346 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7356 std::vector<bool> treated_vertices(
triangulation.n_vertices(),
false);
7357 for (; cell != endc; ++cell)
7360 if (cell->face(f)->at_boundary())
7362 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face;
7365 unsigned int vv = cell->face(f)->vertex_index(v);
7366 if (treated_vertices[vv] ==
false)
7368 treated_vertices[vv] =
true;
7369 for (
unsigned int i = 0; i <= Nz; ++i)
7371 double d = i *
L / Nz;
7372 switch (vv - i * 16)
7375 cell->face(f)->vertex(v) =
7379 cell->face(f)->vertex(v) =
7383 cell->face(f)->vertex(v) =
7387 cell->face(f)->vertex(v) =
7398 double eps = 1
e-3 * outer_radius;
7400 for (; cell != endc; ++cell)
7403 if (cell->face(f)->at_boundary())
7405 double dx = cell->face(f)->center()(0);
7406 double dy = cell->face(f)->center()(1);
7407 double dz = cell->face(f)->center()(2);
7412 cell->face(f)->set_boundary_id(0);
7415 cell->face(f)->set_boundary_id(1);
7418 cell->face(f)->set_boundary_id(2);
7421 cell->face(f)->set_boundary_id(3);
7424 cell->face(f)->set_boundary_id(4);
7427 cell->face(f)->set_boundary_id(5);
7431 cell->face(f)->set_all_boundary_ids(6);
7432 cell->face(f)->set_all_manifold_ids(0);
7439 double d = c.
norm();
7440 if (
d - inner_radius < 0)
7442 cell->face(f)->set_all_boundary_ids(1);
7443 cell->face(f)->set_all_manifold_ids(0);
7446 cell->face(f)->set_boundary_id(0);
7453 template <
int dim,
int spacedim1,
int spacedim2>
7466 "Cannot use this function on parallel::distributed::Triangulation."));
7468 std::vector<Point<spacedim2>> v;
7469 std::vector<CellData<dim>> cells;
7472 const unsigned int spacedim =
std::min(spacedim1, spacedim2);
7473 const std::vector<Point<spacedim1>> &in_vertices = in_tria.
get_vertices();
7475 v.resize(in_vertices.size());
7476 for (
unsigned int i = 0; i < in_vertices.size(); ++i)
7477 for (
unsigned int d = 0;
d < spacedim; ++
d)
7478 v[i][
d] = in_vertices[i][
d];
7483 endc = in_tria.
end();
7485 for (
unsigned int id = 0; cell != endc; ++cell, ++id)
7487 cells[id].vertices.resize(cell->n_vertices());
7488 for (
const auto i : cell->vertex_indices())
7489 cells[id].vertices[i] = cell->vertex_index(i);
7490 cells[id].material_id = cell->material_id();
7491 cells[id].manifold_id = cell->manifold_id();
7507 for (; face != endf; ++face)
7508 if (face->at_boundary())
7511 face->n_vertices());
7512 for (
const auto i : face->vertex_indices())
7514 face->vertex_index(i);
7516 face->boundary_id();
7518 face->manifold_id();
7527 for (; face != endf; ++face)
7528 if (face->at_boundary())
7531 face->n_vertices());
7532 for (
const auto i : face->vertex_indices())
7534 face->vertex_index(i);
7536 face->boundary_id();
7538 face->manifold_id();
7553 template <
int dim,
int spacedim>
7559 ExcMessage(
"Number of global levels has to be 1."));
7587 {{{0, 1, 12, 10}}, {{2, 3, 11, 12}}, {{7, 6, 11, 13}},
7588 {{5, 4, 13, 10}}, {{0, 2, 8, 12}}, {{4, 6, 13, 8}},
7589 {{5, 13, 7, 9}}, {{1, 9, 3, 12}}, {{0, 8, 4, 10}},
7590 {{1, 5, 9, 10}}, {{3, 7, 11, 9}}, {{2, 6, 8, 11}},
7591 {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
7592 {{12, 13, 8, 10}}, {{13, 8, 10, 4}}, {{13, 10, 9, 5}},
7593 {{13, 9, 11, 7}}, {{13, 11, 8, 6}}, {{10, 12, 9, 1}},
7594 {{9, 12, 11, 3}}, {{11, 12, 8, 2}}, {{8, 12, 10, 0}}}};
7602 vertex_ids_for_boundary_faces_2d = {{{{{{0, 4}}, {{4, 2}}}},
7603 {{{{1, 5}}, {{5, 3}}}},
7604 {{{{0, 6}}, {{6, 1}}}},
7605 {{{{2, 7}}, {{7, 3}}}}}};
7613 vertex_ids_for_boundary_faces_3d = {
7614 {{{{{0, 4, 8}}, {{4, 8, 6}}, {{8, 6, 2}}, {{0, 2, 8}}}},
7615 {{{{1, 3, 9}}, {{3, 9, 7}}, {{9, 7, 5}}, {{1, 9, 5}}}},
7616 {{{{0, 1, 10}}, {{1, 10, 5}}, {{10, 5, 4}}, {{0, 10, 4}}}},
7617 {{{{2, 3, 11}}, {{3, 11, 7}}, {{11, 7, 6}}, {{2, 11, 6}}}},
7618 {{{{0, 1, 12}}, {{1, 12, 3}}, {{12, 3, 2}}, {{0, 12, 2}}}},
7619 {{{{4, 5, 13}}, {{5, 13, 7}}, {{13, 7, 6}}, {{4, 13, 6}}}}}};
7640 {{{0, 12, 10}}, {{12, 1, 10}}, {{12, 1, 9}}, {{12, 3, 9}},
7641 {{12, 2, 11}}, {{12, 3, 11}}, {{12, 0, 8}}, {{12, 2, 8}},
7642 {{9, 13, 5}}, {{13, 7, 9}}, {{11, 7, 13}}, {{11, 6, 13}},
7643 {{4, 8, 13}}, {{6, 8, 13}}, {{4, 13, 10}}, {{13, 5, 10}},
7644 {{10, 9, 5}}, {{10, 9, 1}}, {{11, 9, 7}}, {{11, 9, 3}},
7645 {{8, 11, 2}}, {{8, 11, 6}}, {{8, 10, 0}}, {{8, 10, 4}},
7646 {{12, 3, 9}}, {{12, 9, 11}}, {{12, 3, 11}}, {{3, 9, 11}},
7647 {{2, 12, 8}}, {{2, 12, 11}}, {{2, 11, 8}}, {{8, 12, 11}},
7648 {{0, 12, 10}}, {{0, 12, 8}}, {{0, 8, 10}}, {{8, 10, 12}},
7649 {{12, 1, 10}}, {{12, 1, 9}}, {{1, 10, 9}}, {{10, 9, 12}},
7650 {{10, 8, 4}}, {{10, 8, 13}}, {{4, 13, 8}}, {{4, 13, 10}},
7651 {{10, 9, 13}}, {{10, 9, 5}}, {{13, 5, 10}}, {{13, 5, 9}},
7652 {{13, 7, 9}}, {{13, 7, 11}}, {{9, 11, 13}}, {{9, 11, 7}},
7653 {{8, 11, 13}}, {{8, 11, 6}}, {{6, 13, 8}}, {{6, 13, 11}},
7654 {{12, 13, 10}}, {{12, 13, 8}}, {{8, 10, 13}}, {{8, 10, 12}},
7655 {{12, 13, 10}}, {{12, 13, 9}}, {{10, 9, 13}}, {{10, 9, 12}},
7656 {{12, 13, 9}}, {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
7657 {{12, 13, 11}}, {{12, 13, 8}}, {{8, 11, 13}}, {{8, 11, 12}}}};
7664 {{{12, 10}}, {{12, 9}}, {{12, 11}}, {{12, 8}}, {{9, 13}}, {{11, 13}},
7665 {{8, 13}}, {{10, 13}}, {{10, 9}}, {{9, 11}}, {{11, 8}}, {{8, 10}},
7666 {{12, 9}}, {{12, 11}}, {{11, 9}}, {{12, 8}}, {{12, 11}}, {{11, 8}},
7667 {{12, 8}}, {{12, 10}}, {{10, 8}}, {{12, 10}}, {{12, 9}}, {{9, 10}},
7668 {{13, 10}}, {{13, 8}}, {{8, 10}}, {{13, 10}}, {{13, 9}}, {{9, 10}},
7669 {{13, 11}}, {{13, 9}}, {{11, 9}}, {{13, 11}}, {{13, 8}}, {{11, 8}},
7670 {{12, 13}}, {{8, 10}}, {{8, 13}}, {{10, 13}}, {{8, 12}}, {{10, 12}},
7671 {{12, 13}}, {{10, 9}}, {{10, 13}}, {{9, 13}}, {{10, 12}}, {{9, 12}},
7672 {{12, 13}}, {{9, 11}}, {{9, 13}}, {{11, 13}}, {{9, 12}}, {{11, 12}},
7673 {{12, 13}}, {{11, 8}}, {{11, 13}}, {{8, 13}}, {{11, 12}}, {{8, 12}}}};
7681 vertex_ids_for_boundary_edges_3d = {{{{{{4, 6}},
7731 std::vector<Point<spacedim>>
vertices;
7732 std::vector<CellData<dim>> cells;
7737 std::vector<unsigned int> old_to_new_vertex_indices(
7739 std::vector<unsigned int> face_to_new_vertex_indices(
7749 for (
const auto &cell : in_tria)
7753 std::array<unsigned int, dim == 2 ? 9 : 14> local_vertex_indices;
7756 for (
const auto v : cell.vertex_indices())
7758 const auto v_global = cell.vertex_index(v);
7760 if (old_to_new_vertex_indices[v_global] ==
7763 old_to_new_vertex_indices[v_global] =
vertices.size();
7764 vertices.push_back(cell.vertex(v));
7767 local_vertex_indices[v] = old_to_new_vertex_indices[v_global];
7771 for (
const auto f : cell.face_indices())
7773 const auto f_global = cell.face_index(f);
7775 if (face_to_new_vertex_indices[f_global] ==
7778 face_to_new_vertex_indices[f_global] =
vertices.size();
7780 cell.face(f)->center(
true));
7783 local_vertex_indices[cell.n_vertices() + f] =
7784 face_to_new_vertex_indices[f_global];
7790 local_vertex_indices[cell.n_vertices() + cell.n_faces()] =
7792 vertices.push_back(cell.center(
true));
7796 const auto add_cell = [&](
const unsigned int struct_dim,
7797 const auto & index_vertices,
7798 const unsigned int material_or_boundary_id,
7802 if (struct_dim < dim &&
7807 if (struct_dim == dim)
7819 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
7821 cell_data.vertices[i] =
7822 local_vertex_indices[index_vertices[i]];
7823 cell_data.material_id =
7824 material_or_boundary_id;
7825 cell_data.manifold_id =
7828 cells.push_back(cell_data);
7830 else if (dim == 2 && struct_dim == 1)
7834 boundary_line.
boundary_id = material_or_boundary_id;
7836 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
7839 local_vertex_indices[index_vertices[i]];
7843 else if (dim == 3 && struct_dim == 2)
7847 boundary_quad.
material_id = material_or_boundary_id;
7849 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
7852 local_vertex_indices[index_vertices[i]];
7856 else if (dim == 3 && struct_dim == 1)
7861 for (
unsigned int i = 0; i < index_vertices.size(); ++i)
7864 local_vertex_indices[index_vertices[i]];
7874 const auto material_id_cell = cell.material_id();
7880 const auto manifold_id_cell = cell.manifold_id();
7882 for (
const auto &cell_vertices : table_2D_cell)
7883 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
7886 for (
const auto &face_vertices : vertex_ids_for_inner_faces_2d)
7898 const auto manifold_id_cell = cell.manifold_id();
7900 for (
const auto &cell_vertices : vertex_ids_for_cells_3d)
7901 add_cell(dim, cell_vertices, material_id_cell, manifold_id_cell);
7905 for (
const auto &face_vertices : vertex_ids_for_inner_faces_3d)
7913 for (
const auto &edge_vertices : vertex_ids_for_inner_edges_3d)
7923 for (
const auto f : cell.face_indices())
7925 const auto bid = cell.face(f)->boundary_id();
7926 const auto mid = cell.face(f)->manifold_id();
7930 for (
const auto &face_vertices :
7931 vertex_ids_for_boundary_faces_2d[f])
7932 add_cell(1, face_vertices, bid, mid);
7938 for (
const auto &face_vertices :
7939 vertex_ids_for_boundary_faces_3d[f])
7940 add_cell(2, face_vertices, bid, mid);
7943 for (
const auto &edge_vertices :
7944 vertex_ids_for_boundary_edges_3d[f])
7945 add_cell(1, edge_vertices, bid, mid);
7958 template <
int spacedim>
7969 template <
template <
int,
int>
class MeshType,
int dim,
int spacedim>
7971 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
7972 typename MeshType<dim, spacedim>::face_iterator>
7974 typename ExtractBoundaryMesh<MeshType, dim, spacedim>::return_type
7977 MeshType<dim - 1, spacedim> & surface_mesh,
7978 const std::set<types::boundary_id> &boundary_ids)
7982 &volume_mesh.get_triangulation()) ==
nullptr),
7991 const unsigned int boundary_dim = dim - 1;
7998 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>
7999 temporary_mapping_level0;
8004 std::vector<bool> touched(volume_mesh.get_triangulation().n_vertices(),
8008 std::vector<CellData<boundary_dim>> cells;
8010 std::vector<Point<spacedim>>
vertices;
8013 std::map<unsigned int, unsigned int> map_vert_index;
8026 for (
unsigned int i1 = 0; i1 < GeometryInfo<spacedim>::faces_per_cell; i1++)
8028 for (
unsigned int i2 = 0; i2 <
GeometryInfo<dim - 1>::vertices_per_cell;
8030 swap_matrix[i1][i2] = i2;
8036 std::swap(swap_matrix[0][1], swap_matrix[0][2]);
8037 std::swap(swap_matrix[2][1], swap_matrix[2][2]);
8038 std::swap(swap_matrix[4][1], swap_matrix[4][2]);
8042 std::swap(swap_matrix[1][0], swap_matrix[1][1]);
8043 std::swap(swap_matrix[2][0], swap_matrix[2][1]);
8048 for (
typename MeshType<dim, spacedim>::cell_iterator cell =
8049 volume_mesh.begin(0);
8050 cell != volume_mesh.end(0);
8054 const typename MeshType<dim, spacedim>::face_iterator face =
8057 if (face->at_boundary() &&
8058 (boundary_ids.empty() ||
8059 (boundary_ids.find(face->boundary_id()) != boundary_ids.end())))
8063 for (
const unsigned int j :
8066 const unsigned int v_index = face->vertex_index(j);
8068 if (!touched[v_index])
8070 vertices.push_back(face->vertex(j));
8071 map_vert_index[v_index] =
vertices.size() - 1;
8072 touched[v_index] =
true;
8075 c_data.
vertices[swap_matrix[i][j]] = map_vert_index[v_index];
8089 for (
unsigned int e = 0;
e < 4; ++
e)
8095 bool edge_found =
false;
8098 map_vert_index[face->line(
e)->vertex_index(0)]) &&
8100 map_vert_index[face->line(
e)->vertex_index(
8103 map_vert_index[face->line(
e)->vertex_index(1)]) &&
8105 map_vert_index[face->line(
e)->vertex_index(0)])))
8112 if (edge_found ==
true)
8119 map_vert_index[face->line(
e)->vertex_index(0)];
8121 map_vert_index[face->line(
e)->vertex_index(1)];
8128 cells.push_back(c_data);
8129 temporary_mapping_level0.push_back(std::make_pair(face, i));
8136 surface_mesh.get_triangulation())
8142 for (
const auto &cell : surface_mesh.active_cell_iterators())
8143 for (
unsigned int vertex = 0; vertex < 2; vertex++)
8144 if (cell->face(vertex)->at_boundary())
8145 cell->face(vertex)->set_boundary_id(0);
8153 std::vector<std::pair<
8154 const typename MeshType<dim - 1, spacedim>::cell_iterator,
8155 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>>
8156 temporary_map_boundary_cell_face;
8157 for (
const auto &cell : surface_mesh.active_cell_iterators())
8158 temporary_map_boundary_cell_face.push_back(
8159 std::make_pair(cell, temporary_mapping_level0.at(cell->index())));
8171 unsigned int index_cells_deepest_level = 0;
8174 bool changed =
false;
8178 std::vector<unsigned int> cells_refined;
8181 for (
unsigned int cell_n = index_cells_deepest_level;
8182 cell_n < temporary_map_boundary_cell_face.size();
8187 if (temporary_map_boundary_cell_face[cell_n]
8188 .
second.first->has_children())
8192 Assert(temporary_map_boundary_cell_face[cell_n]
8193 .
second.first->refinement_case() ==
8196 temporary_map_boundary_cell_face[cell_n]
8197 .first->set_refine_flag();
8198 cells_refined.push_back(cell_n);
8209 surface_mesh.get_triangulation())
8210 .execute_coarsening_and_refinement();
8213 index_cells_deepest_level = temporary_map_boundary_cell_face.size();
8214 for (
const auto &refined_cell_n : cells_refined)
8216 const typename MeshType<dim - 1, spacedim>::cell_iterator
8218 temporary_map_boundary_cell_face[refined_cell_n].first;
8219 const typename MeshType<dim,
8220 spacedim>::face_iterator refined_face =
8221 temporary_map_boundary_cell_face[refined_cell_n].second.first;
8222 const unsigned int refined_face_number =
8223 temporary_map_boundary_cell_face[refined_cell_n]
8225 for (
unsigned int child_n = 0;
8226 child_n < refined_cell->n_children();
8231 temporary_map_boundary_cell_face.push_back(
8232 std::make_pair(refined_cell->child(
8233 swap_matrix[refined_face_number][child_n]),
8234 std::make_pair(refined_face->child(child_n),
8235 refined_face_number)));
8245 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8246 typename MeshType<dim, spacedim>::face_iterator>
8247 surface_to_volume_mapping;
8248 for (
unsigned int i = 0; i < temporary_map_boundary_cell_face.size(); i++)
8249 surface_to_volume_mapping[temporary_map_boundary_cell_face[i].
first] =
8250 temporary_map_boundary_cell_face[i].
second.first;
8252 return surface_to_volume_mapping;
8257 template <
int dim,
int spacedim>
8261 const std::vector<unsigned int> &repetitions,
8270 std::vector<Point<spacedim>>
vertices;
8271 std::vector<CellData<dim>> cells;
8277 (p2[1] - p1[1]) / repetitions[1]);
8280 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
8281 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
8286 for (
unsigned int j = 0; j < repetitions[1]; ++j)
8287 for (
unsigned int i = 0; i < repetitions[0]; ++i)
8290 std::array<unsigned int, 4> quad{{
8291 (j + 0) * (repetitions[0] + 1) + i + 0,
8292 (j + 0) * (repetitions[0] + 1) + i + 1,
8293 (j + 1) * (repetitions[0] + 1) + i + 0,
8294 (j + 1) * (repetitions[0] + 1) + i + 1
8300 tri.
vertices = {quad[0], quad[1], quad[2]};
8301 cells.push_back(tri);
8307 tri.
vertices = {quad[3], quad[2], quad[1]};
8308 cells.push_back(tri);
8316 (p2[1] - p1[1]) / repetitions[1],
8317 (p2[2] - p1[2]) / repetitions[2]);
8320 for (
unsigned int k = 0; k <= repetitions[2]; ++k)
8321 for (
unsigned int j = 0; j <= repetitions[1]; ++j)
8322 for (
unsigned int i = 0; i <= repetitions[0]; ++i)
8325 p1[2] +
dx[2] * k));
8328 for (
unsigned int k = 0; k < repetitions[2]; ++k)
8329 for (
unsigned int j = 0; j < repetitions[1]; ++j)
8330 for (
unsigned int i = 0; i < repetitions[0]; ++i)
8333 std::array<unsigned int, 8> quad{
8334 {(k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8335 (j + 0) * (repetitions[0] + 1) + i + 0,
8336 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8337 (j + 0) * (repetitions[0] + 1) + i + 1,
8338 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8339 (j + 1) * (repetitions[0] + 1) + i + 0,
8340 (k + 0) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8341 (j + 1) * (repetitions[0] + 1) + i + 1,
8342 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8343 (j + 0) * (repetitions[0] + 1) + i + 0,
8344 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8345 (j + 0) * (repetitions[0] + 1) + i + 1,
8346 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8347 (j + 1) * (repetitions[0] + 1) + i + 0,
8348 (k + 1) * (repetitions[0] + 1) * (repetitions[1] + 1) +
8349 (j + 1) * (repetitions[0] + 1) + i + 1}};
8354 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8355 cell.
vertices = {{quad[0], quad[1], quad[2], quad[4]}};
8357 cell.
vertices = {{quad[0], quad[1], quad[3], quad[5]}};
8359 cells.push_back(cell);
8365 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8366 cell.
vertices = {{quad[2], quad[1], quad[3], quad[7]}};
8368 cell.
vertices = {{quad[0], quad[3], quad[2], quad[6]}};
8369 cells.push_back(cell);
8375 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8376 cell.
vertices = {{quad[1], quad[4], quad[5], quad[7]}};
8378 cell.
vertices = {{quad[0], quad[4], quad[5], quad[6]}};
8379 cells.push_back(cell);
8385 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8386 cell.
vertices = {{quad[2], quad[4], quad[7], quad[6]}};
8388 cell.
vertices = {{quad[3], quad[5], quad[7], quad[6]}};
8389 cells.push_back(cell);
8395 if (((i % 2) + (j % 2) + (k % 2)) % 2 == 0)
8396 cell.
vertices = {{quad[1], quad[2], quad[4], quad[7]}};
8398 cell.
vertices = {{quad[0], quad[3], quad[6], quad[5]}};
8399 cells.push_back(cell);
8414 template <
int dim,
int spacedim>
8417 const unsigned int repetitions,
8425 tria, {{repetitions, repetitions}}, {p1, p1}, {p2, p2},
colorize);
8431 {{repetitions, repetitions, repetitions}},
8445# 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)
numbers::NumberTraits< Number >::real_type distance(const Point< dim, Number > &p) const
static Point< dim, Number > unit_vector(const unsigned int i)
const Point< spacedim > center
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
numbers::NumberTraits< Number >::real_type norm() const
constexpr Tensor< 1, dim, typename ProductType< Number1, Number2 >::type > cross_product_3d(const Tensor< 1, dim, Number1 > &src1, const Tensor< 1, dim, Number2 > &src2)
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
face_iterator end_face() const
cell_iterator begin(const unsigned int level=0) const
unsigned int n_active_faces() 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
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
unsigned int n_vertices() const
active_face_iterator begin_active_face() const
active_cell_iterator begin_active(const unsigned int level=0) const
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
IteratorRange< active_face_iterator > active_face_iterators() const
IteratorRange< active_cell_iterator > active_cell_iterators() const
IteratorRange< cell_iterator > cell_iterators() const
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInvalidInputOrientation()
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
static ::ExceptionBase & ExcInvalidRepetitions(int arg1)
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcLowerRange(int arg1, int arg2)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcInvalidRadii()
static ::ExceptionBase & ExcInvalidRepetitionsDimension(int arg1)
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 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 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)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
Expression fabs(const Expression &x)
Expression ceil(const Expression &x)
Expression atan(const Expression &x)
Expression tanh(const Expression &x)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void parallelepiped(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void hyper_cross(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &sizes, const bool colorize_cells=false)
A center cell with stacks of cell protruding from each surface.
std::map< typename MeshType< dim - 1, spacedim >::cell_iterator, typename MeshType< dim, spacedim >::face_iterator > extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
void 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 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)
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
Rectangular domain with rectangular pattern of holes.
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)
void non_standard_orientation_mesh(Triangulation< 2 > &tria, const bool rotate_left_square, const bool rotate_right_square)
void create_union_triangulation(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result)
void subdivided_parallelepiped(Triangulation< dim > &tria, const unsigned int n_subdivisions, const Point< dim >(&corners)[dim], const bool colorize=false)
void subdivided_cylinder(Triangulation< dim > &tria, const unsigned int x_subdivisions, const double radius=1., const double half_length=1.)
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
void subdivided_hyper_L(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &bottom_left, const Point< dim > &top_right, const std::vector< int > &n_cells_to_remove)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=Point< spacedim >(), const double radius=1.)
void concentric_hyper_shells(Triangulation< dim > &triangulation, const Point< dim > ¢er, const double inner_radius=0.125, const double outer_radius=0.25, const unsigned int n_shells=1, const double skewness=0.1, const unsigned int n_cells_per_shell=0, const bool colorize=false)
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void 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 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 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)
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
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 > 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 const ReferenceCell Tetrahedron
constexpr const ReferenceCell Wedge
constexpr const ReferenceCell Pyramid
constexpr const ReferenceCell Triangle
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
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
const types::boundary_id invalid_boundary_id
static constexpr double E
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
::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