50 : airfoil_type(
"NACA")
52 , joukowski_center(-0.1, 0.14)
58 , incline_factor(0.35)
61 , n_subdivision_x_0(3)
62 , n_subdivision_x_1(2)
63 , n_subdivision_x_2(5)
65 , airfoil_sampling_factor(2)
68 airfoil_length <= height,
70 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
73 ExcMessage(
"incline_factor has to be in [0,1)!"));
86 "Mesh height measured from airfoil nose to horizontal boundaries");
90 "Length measured from airfoil leading edge to vertical outlet boundary");
94 "Define obliqueness of the vertical mesh around the airfoil");
103 "Type of airfoil geometry, either NACA or Joukowski airfoil",
118 "Joukowski circle center coordinates");
121 "Joukowski airfoil length leading to trailing edge");
129 "Number of global refinements");
131 "NumberSubdivisionX0",
133 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
135 "NumberSubdivisionX1",
137 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
139 "NumberSubdivisionX2",
141 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
144 "Number of subdivisions normal to airfoil");
148 "Factor to obtain a finer mesh at the airfoil surface");
174 : refinements(
data.refinements)
175 , n_subdivision_x_0(
data.n_subdivision_x_0)
176 , n_subdivision_x_1(
data.n_subdivision_x_1)
177 , n_subdivision_x_2(
data.n_subdivision_x_2)
178 , n_subdivision_y(
data.n_subdivision_y)
179 , height(
data.height)
180 , length_b2(
data.length_b2)
181 , incline_factor(
data.incline_factor)
182 , bias_factor(
data.bias_factor)
192 data.airfoil_type ==
"Joukowski" ?
195 data.airfoil_sampling_factor) :
196 (
data.airfoil_type ==
"NACA" ?
199 data.airfoil_sampling_factor) :
202 std::vector<Point<2>>{
206 data.airfoil_length))
228 , J(
tail_x + length_b2, 0)
233 data.airfoil_type ==
"NACA",
278 const unsigned int refinements;
281 const unsigned int n_subdivision_x_0;
284 const unsigned int n_subdivision_x_1;
287 const unsigned int n_subdivision_x_2;
291 const unsigned int n_subdivision_y;
298 const double length_b2;
302 const double incline_factor;
305 const double bias_factor;
328 const std::array<std::vector<Point<2>>, 2>
airfoil_1D;
380 const Point<2> A,
B,
C,
D,
E,
F, G,
H, I, J,
K,
L;
419 static std::array<std::vector<Point<2>>, 2>
422 const unsigned int factor)
424 std::array<std::vector<Point<2>>, 2>
airfoil_1D;
444 for (
unsigned int i = 0; i <
jouk_points.size(); ++i)
461 for (
unsigned int i = 0; i <=
nose_index; ++i)
477 for (
auto &i : vector)
511 static std::vector<Point<2>>
522 const double radius =
std::sqrt(center[1] * center[1] +
523 (1 - center[0]) * (1 - center[0]));
525 center[1] * center[1] + (1 + center[0]) * (1 + center[0]));
531 "Error creating lower circle: Circle for Joukowski-transform does"
532 " not enclose point zeta = -1! Choose different center "
540 center[1] - radius *
sin(i * theta));
553 static std::vector<Point<2>>
564 const std::complex<double> z =
zeta + 1. /
zeta;
587 static std::array<std::vector<Point<2>>, 2>
590 const unsigned int factor)
616 static std::vector<Point<2>>
622 ExcMessage(
"This NACA-serial number is not implemented!"));
643 static std::vector<Point<2>>
657 const double t =
static_cast<double>(
digit_23) / 100.0;
668 0.3516 * Utilities::fixed_power<2>(
x) +
669 0.2843 * Utilities::fixed_power<3>(
x) -
670 0.1036 * Utilities::fixed_power<4>(
681 const double m = 1.0 *
digit_0 / 100;
682 const double p = 1.0 *
digit_1 / 10;
687 m / Utilities::fixed_power<2>(p) *
688 (2 * p *
x - Utilities::fixed_power<2>(
x)) :
693 (
x <= p) ? 2 * m / Utilities::fixed_power<2>(p) * (p -
x) :
699 0.3516 * Utilities::fixed_power<2>(
x) +
700 0.2843 * Utilities::fixed_power<3>(
x) -
701 0.1036 * Utilities::fixed_power<4>(
727 static std::array<std::vector<Point<2>>, 2>
731 std::array<std::vector<Point<2>>, 2> output;
745 static std::vector<Point<2>>
749 std::vector<Point<2>> output = input;
752 desired_len / input.front().distance(input.back());
754 for (
auto &
x : output)
770 static std::vector<Point<2>>
775 const unsigned int n_points =
787 const auto airfoil_length =
801 for (
unsigned int j = 0, i = 1;
j < n_points - 1; ++
j)
833 std::vector<Triangulation<2>>
trias(10);
848 for (
auto it = tria.begin_vertex();
it < tria.end_vertex(); ++
it)
862 for (
auto cell : tria.active_cell_iterators())
870 {n_subdivision_y, n_subdivision_x_0},
874 {n_subdivision_y, n_subdivision_x_0},
878 {n_subdivision_x_1, n_subdivision_y},
882 {n_subdivision_x_1, n_subdivision_y},
886 {n_subdivision_x_2, n_subdivision_y},
890 {n_subdivision_x_2, n_subdivision_y},
917 for (
auto cell : tria.active_cell_iterators())
920 if (cell->face(f)->at_boundary() ==
false)
923 const auto mid = cell->material_id();
927 cell->face(f)->set_boundary_id(0);
930 cell->face(f)->set_boundary_id(1);
933 cell->face(f)->set_boundary_id(2);
936 cell->face(f)->set_boundary_id(3);
939 cell->face(f)->set_boundary_id(4);
942 cell->face(f)->set_boundary_id(5);
979 for (
const auto &cell : tria.cell_iterators())
989 auto &
node = cell->vertex(v);
1022 const double L = height /
std::sin(gamma);
1034 static_cast<int>(std::round((
node_[0] -
deltax) / dx));
1048 static_cast<int>(std::round(
node_[0] / dx));
1050 static_cast<int>(std::round(
node_[1] / dy));
1051 const double alpha =
1053 const double theta =
node_[0];
1060 (cell->material_id() ==
id_block_1) ? (0) : (1))]
1066 else if (cell->material_id() ==
id_block_2 ||
1075 "Points D,C,G and E,F,I are not defined symmetric to "
1076 "x-axis, which is required to interpolate block 2"
1077 " and 5 with same geometric computations."));
1078 const double l_y =
D[1] -
C[1];
1079 const double l_h =
D[1] -
l_y;
1082 const int iy =
static_cast<int>(
1085 const int ix =
static_cast<int>(
1093 incline_factor * length_b2 *
ix /
1101 (cell->material_id() ==
id_block_2) ? (0) : (1))]
1106 else if (cell->material_id() ==
id_block_3 ||
1112 const int ix =
static_cast<int>(
1123 const Point<2> p1(J[0] - (1 - incline_factor) * length_b2 *
1165 const AdditionalData &additional_data)
1196 const AdditionalData &)
1206 const AdditionalData &additional_data)
1219 const AdditionalData &additional_data)
1232 const AdditionalData &additional_data)
1236 (
void)additional_data;
1248 template <
int dim,
int spacedim>
1260 cell->face(f)->set_boundary_id(f);
1266 template <
int spacedim>
1277 if (cell->center()[0] > 0)
1278 cell->set_material_id(1);
1285 template <
int dim,
int spacedim>
1290 const double epsilon)
1303 for (; face !=
endface; ++face)
1304 if (face->at_boundary())
1305 if (face->boundary_id() == 0)
1310 face->set_boundary_id(0);
1311 else if (
std::abs(center[0] -
p2[0]) < epsilon)
1312 face->set_boundary_id(1);
1313 else if (dim > 1 &&
std::abs(center[1] -
p1[1]) < epsilon)
1314 face->set_boundary_id(2);
1315 else if (dim > 1 &&
std::abs(center[1] -
p2[1]) < epsilon)
1316 face->set_boundary_id(3);
1317 else if (dim > 2 &&
std::abs(center[2] -
p1[2]) < epsilon)
1318 face->set_boundary_id(4);
1319 else if (dim > 2 &&
std::abs(center[2] -
p2[2]) < epsilon)
1320 face->set_boundary_id(5);
1329 for (
const auto &cell : tria.cell_iterators())
1332 for (
unsigned int d = 0;
d < dim; ++
d)
1333 if (cell->center()[d] > 0)
1335 cell->set_material_id(
id);
1344 template <
int spacdim>
1357 for (
auto cell = tria.begin(); cell != tria.end(); ++cell)
1360 cell->face(2)->set_all_boundary_ids(1);
1379 if (tria.n_cells() == 6)
1384 cell->face(4)->set_all_boundary_ids(1);
1388 cell->face(2)->set_all_boundary_ids(1);
1392 cell->face(2)->set_all_boundary_ids(1);
1396 cell->face(0)->set_all_boundary_ids(1);
1400 cell->face(2)->set_all_boundary_ids(1);
1404 cell->face(0)->set_all_boundary_ids(1);
1406 else if (tria.n_cells() == 12)
1415 cell->face(5)->set_all_boundary_ids(1);
1418 else if (tria.n_cells() == 96)
1425 unsigned int count = 0;
1426 for (
const auto &cell : tria.cell_iterators())
1427 if (cell->face(5)->at_boundary())
1429 cell->face(5)->set_all_boundary_ids(1);
1449 const double inner_radius,
1452 if (tria.n_cells() != 3)
1459 for (; cell != tria.end(); ++cell)
1460 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
1462 if (!cell->face(f)->at_boundary())
1465 double radius = cell->face(f)->center().norm() - center.
norm();
1466 if (std::fabs(cell->face(f)->center()[0]) <
1469 cell->face(f)->set_boundary_id(2);
1472 if (cell->face(f)->line(
j)->at_boundary())
1473 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1474 cell->face(f)->line(
j)->vertex(1).norm()) >
1476 cell->face(f)->line(
j)->set_boundary_id(2);
1478 else if (std::fabs(cell->face(f)->center()[1]) <
1481 cell->face(f)->set_boundary_id(3);
1484 if (cell->face(f)->line(
j)->at_boundary())
1485 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1486 cell->face(f)->line(
j)->vertex(1).norm()) >
1488 cell->face(f)->line(
j)->set_boundary_id(3);
1490 else if (std::fabs(cell->face(f)->center()[2]) <
1493 cell->face(f)->set_boundary_id(4);
1496 if (cell->face(f)->line(
j)->at_boundary())
1497 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1498 cell->face(f)->line(
j)->vertex(1).norm()) >
1500 cell->face(f)->line(
j)->set_boundary_id(4);
1502 else if (radius <
middle)
1504 cell->face(f)->set_boundary_id(0);
1507 if (cell->face(f)->line(
j)->at_boundary())
1508 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1509 cell->face(f)->line(
j)->vertex(1).norm()) <
1511 cell->face(f)->line(
j)->set_boundary_id(0);
1513 else if (radius >
middle)
1515 cell->face(f)->set_boundary_id(1);
1518 if (cell->face(f)->line(
j)->at_boundary())
1519 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1520 cell->face(f)->line(
j)->vertex(1).norm()) <
1522 cell->face(f)->line(
j)->set_boundary_id(1);
1532 template <
int dim,
int spacedim>
1543 for (
unsigned int i = 0; i < dim; ++i)
1557 vertices[0] = vertices[1] =
p1;
1558 vertices[2] = vertices[3] =
p2;
1560 vertices[1][0] =
p2[0];
1561 vertices[2][0] =
p1[0];
1564 vertices[0] = vertices[1] = vertices[2] = vertices[3] =
p1;
1565 vertices[4] = vertices[5] = vertices[6] = vertices[7] =
p2;
1567 vertices[1][0] =
p2[0];
1568 vertices[2][1] =
p2[1];
1569 vertices[3][0] =
p2[0];
1570 vertices[3][1] =
p2[1];
1572 vertices[4][0] =
p1[0];
1573 vertices[4][1] =
p1[1];
1574 vertices[5][1] =
p1[1];
1575 vertices[6][0] =
p1[0];
1583 std::vector<CellData<dim>> cells(1);
1585 cells[0].vertices[i] = i;
1586 cells[0].material_id = 0;
1588 tria.create_triangulation(vertices, cells,
SubCellData());
1597 template <
int dim,
int spacedim>
1605 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1608 for (
unsigned int i = 0; i < dim; ++i)
1628 for (
unsigned int d = 0;
d < dim; ++
d)
1629 for (
unsigned int c = 1; c <= dim; ++c)
1632 ExcMessage(
"Vertices of simplex must form a right handed system"));
1636 std::vector<Point<dim>> points = vertices;
1640 for (
unsigned int i = 0; i <= dim; ++i)
1642 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1643 center += points[i];
1648 for (
unsigned int i = 1; i < dim; ++i)
1649 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1651 for (
unsigned int i = 0; i <= dim; ++i)
1652 points.push_back(1. / 3. *
1653 (points[i] + points[(i + 1) % (dim + 1)] +
1654 points[(i + 2) % (dim + 1)]));
1656 points.push_back((1. / (dim + 1)) * center);
1658 std::vector<CellData<dim>> cells(dim + 1);
1663 cells[0].vertices[0] = 0;
1664 cells[0].vertices[1] = 3;
1665 cells[0].vertices[2] = 5;
1666 cells[0].vertices[3] = 6;
1667 cells[0].material_id = 0;
1669 cells[1].vertices[0] = 3;
1670 cells[1].vertices[1] = 1;
1671 cells[1].vertices[2] = 6;
1672 cells[1].vertices[3] = 4;
1673 cells[1].material_id = 0;
1675 cells[2].vertices[0] = 5;
1676 cells[2].vertices[1] = 6;
1677 cells[2].vertices[2] = 2;
1678 cells[2].vertices[3] = 4;
1679 cells[2].material_id = 0;
1683 cells[0].vertices[0] = 0;
1684 cells[0].vertices[1] = 4;
1685 cells[0].vertices[2] = 8;
1686 cells[0].vertices[3] = 10;
1687 cells[0].vertices[4] = 7;
1688 cells[0].vertices[5] = 13;
1689 cells[0].vertices[6] = 12;
1690 cells[0].vertices[7] = 14;
1691 cells[0].material_id = 0;
1693 cells[1].vertices[0] = 4;
1694 cells[1].vertices[1] = 1;
1695 cells[1].vertices[2] = 10;
1696 cells[1].vertices[3] = 5;
1697 cells[1].vertices[4] = 13;
1698 cells[1].vertices[5] = 9;
1699 cells[1].vertices[6] = 14;
1700 cells[1].vertices[7] = 11;
1701 cells[1].material_id = 0;
1703 cells[2].vertices[0] = 8;
1704 cells[2].vertices[1] = 10;
1705 cells[2].vertices[2] = 2;
1706 cells[2].vertices[3] = 5;
1707 cells[2].vertices[4] = 12;
1708 cells[2].vertices[5] = 14;
1709 cells[2].vertices[6] = 6;
1710 cells[2].vertices[7] = 11;
1711 cells[2].material_id = 0;
1713 cells[3].vertices[0] = 7;
1714 cells[3].vertices[1] = 13;
1715 cells[3].vertices[2] = 12;
1716 cells[3].vertices[3] = 14;
1717 cells[3].vertices[4] = 3;
1718 cells[3].vertices[5] = 9;
1719 cells[3].vertices[6] = 6;
1720 cells[3].vertices[7] = 11;
1721 cells[3].material_id = 0;
1726 tria.create_triangulation(points, cells,
SubCellData());
1731 template <
int dim,
int spacedim>
1738 if (reference_cell == ReferenceCells::get_hypercube<dim>())
1747 std::vector<Point<spacedim>> vertices(
reference_cell.n_vertices());
1751 for (
unsigned int d = 0;
d < dim; ++
d)
1760 std::vector<CellData<dim>> cells(1);
1763 cells[0].vertices[v] = v;
1766 tria.create_triangulation(vertices, cells, {});
1772 const unsigned int n_cells,
1777 const unsigned int dim = 3;
1780 "More than 4 cells are needed to create a moebius grid."));
1782 ExcMessage(
"Outer and inner radius must be positive."));
1784 ExcMessage(
"Outer radius must be greater than inner radius."));
1787 std::vector<Point<dim>> vertices(4 * n_cells);
1791 for (
unsigned int i = 0; i <
n_cells; ++i)
1792 for (
unsigned int j = 0;
j < 4; ++
j)
1794 vertices[4 * i +
j][0] =
1798 vertices[4 * i +
j][1] =
1802 vertices[4 * i +
j][2] =
1806 unsigned int offset = 0;
1812 {0, 1, 5, 4, 2, 3, 7, 6}};
1813 std::vector<CellData<dim>> cells(n_cells);
1814 for (
unsigned int i = 0; i <
n_cells; ++i)
1816 for (
unsigned int j = 0;
j < 2; ++
j)
1828 cells[i].material_id = 0;
1842 tria.create_triangulation(vertices, cells,
SubCellData());
1850 const double centerline_radius,
1851 const double inner_radius,
1855 Assert(centerline_radius > inner_radius,
1856 ExcMessage(
"The centerline radius must be greater than the "
1858 Assert(inner_radius > 0.0,
1859 ExcMessage(
"The inner radius must be positive."));
1861 const unsigned int dim = 2;
1862 const unsigned int spacedim = 3;
1863 std::vector<Point<spacedim>> vertices(16);
1865 vertices[0] =
Point<spacedim>(centerline_radius - inner_radius, 0, 0);
1867 vertices[2] =
Point<spacedim>(centerline_radius + inner_radius, 0, 0);
1869 vertices[4] =
Point<spacedim>(0, 0, centerline_radius - inner_radius);
1871 vertices[6] =
Point<spacedim>(0, 0, centerline_radius + inner_radius);
1873 vertices[8] =
Point<spacedim>(-(centerline_radius - inner_radius), 0, 0);
1875 vertices[10] =
Point<spacedim>(-(centerline_radius + inner_radius), 0, 0);
1877 vertices[12] =
Point<spacedim>(0, 0, -(centerline_radius - inner_radius));
1879 vertices[14] =
Point<spacedim>(0, 0, -(centerline_radius + inner_radius));
1882 std::vector<CellData<dim>> cells(16);
1884 cells[0].vertices[0] = 0;
1885 cells[0].vertices[1] = 4;
1886 cells[0].vertices[2] = 3;
1887 cells[0].vertices[3] = 7;
1888 cells[0].material_id = 0;
1890 cells[1].vertices[0] = 1;
1891 cells[1].vertices[1] = 5;
1892 cells[1].vertices[2] = 0;
1893 cells[1].vertices[3] = 4;
1894 cells[1].material_id = 0;
1896 cells[2].vertices[0] = 2;
1897 cells[2].vertices[1] = 6;
1898 cells[2].vertices[2] = 1;
1899 cells[2].vertices[3] = 5;
1900 cells[2].material_id = 0;
1902 cells[3].vertices[0] = 3;
1903 cells[3].vertices[1] = 7;
1904 cells[3].vertices[2] = 2;
1905 cells[3].vertices[3] = 6;
1906 cells[3].material_id = 0;
1908 cells[4].vertices[0] = 4;
1909 cells[4].vertices[1] = 8;
1910 cells[4].vertices[2] = 7;
1911 cells[4].vertices[3] = 11;
1912 cells[4].material_id = 0;
1914 cells[5].vertices[0] = 5;
1915 cells[5].vertices[1] = 9;
1916 cells[5].vertices[2] = 4;
1917 cells[5].vertices[3] = 8;
1918 cells[5].material_id = 0;
1920 cells[6].vertices[0] = 6;
1921 cells[6].vertices[1] = 10;
1922 cells[6].vertices[2] = 5;
1923 cells[6].vertices[3] = 9;
1924 cells[6].material_id = 0;
1926 cells[7].vertices[0] = 7;
1927 cells[7].vertices[1] = 11;
1928 cells[7].vertices[2] = 6;
1929 cells[7].vertices[3] = 10;
1930 cells[7].material_id = 0;
1932 cells[8].vertices[0] = 8;
1933 cells[8].vertices[1] = 12;
1934 cells[8].vertices[2] = 11;
1935 cells[8].vertices[3] = 15;
1936 cells[8].material_id = 0;
1938 cells[9].vertices[0] = 9;
1939 cells[9].vertices[1] = 13;
1940 cells[9].vertices[2] = 8;
1941 cells[9].vertices[3] = 12;
1942 cells[9].material_id = 0;
1944 cells[10].vertices[0] = 10;
1945 cells[10].vertices[1] = 14;
1946 cells[10].vertices[2] = 9;
1947 cells[10].vertices[3] = 13;
1948 cells[10].material_id = 0;
1950 cells[11].vertices[0] = 11;
1951 cells[11].vertices[1] = 15;
1952 cells[11].vertices[2] = 10;
1953 cells[11].vertices[3] = 14;
1954 cells[11].material_id = 0;
1956 cells[12].vertices[0] = 12;
1957 cells[12].vertices[1] = 0;
1958 cells[12].vertices[2] = 15;
1959 cells[12].vertices[3] = 3;
1960 cells[12].material_id = 0;
1962 cells[13].vertices[0] = 13;
1963 cells[13].vertices[1] = 1;
1964 cells[13].vertices[2] = 12;
1965 cells[13].vertices[3] = 0;
1966 cells[13].material_id = 0;
1968 cells[14].vertices[0] = 14;
1969 cells[14].vertices[1] = 2;
1970 cells[14].vertices[2] = 13;
1971 cells[14].vertices[3] = 1;
1972 cells[14].material_id = 0;
1974 cells[15].vertices[0] = 15;
1975 cells[15].vertices[1] = 3;
1976 cells[15].vertices[2] = 14;
1977 cells[15].vertices[3] = 2;
1978 cells[15].material_id = 0;
1981 tria.create_triangulation(vertices, cells,
SubCellData());
1983 tria.set_all_manifold_ids(0);
2003 const double centerline_radius,
2004 const double inner_radius,
2008 Assert(centerline_radius > inner_radius,
2009 ExcMessage(
"The centerline radius must be greater than the "
2011 Assert(inner_radius > 0.0,
2012 ExcMessage(
"The inner radius must be positive."));
2014 ExcMessage(
"Number of cells in toroidal direction has "
2015 "to be at least 3 for a torus of polar extent 2*pi."));
2021 const double a = 1. / (1 +
std::sqrt(2.0));
2045 for (
unsigned int v = 0; v < 8; ++v)
2049 vertices[8 * c + v][1] = vertices[v][1];
2058 for (
unsigned int j = 0;
j < 2; ++
j)
2060 const unsigned int offset =
2064 for (
unsigned int c2 = 0;
c2 < 5; ++
c2)
2065 for (
unsigned int i = 0; i < 4; ++i)
2066 cells[5 * c +
c2].vertices[i +
j * 4] =
2070 cells[5 * c].material_id = 0;
2072 cells[5 * c + 1].material_id = 0;
2073 cells[5 * c + 2].material_id = 1;
2074 cells[5 * c + 3].material_id = 0;
2075 cells[5 * c + 4].material_id = 0;
2078 tria.create_triangulation(vertices, cells,
SubCellData());
2080 tria.reset_all_manifolds();
2081 tria.set_all_manifold_ids(0);
2083 for (
const auto &cell : tria.cell_iterators())
2086 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
2090 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2092 cell->face(f)->set_all_manifold_ids(1);
2097 if (cell->material_id() == 1)
2099 cell->set_all_manifold_ids(2);
2101 cell->set_material_id(0);
2106 tria.set_manifold(2,
2118 template <
int dim,
int spacedim>
2131 tria.begin_active();
2133 cell->vertex(i) = vertices[i];
2139 "The volume of the cell is not greater than zero. "
2140 "This could be due to the wrong ordering of the vertices."));
2171 std::array<Tensor<1, 2>, 2>
edges;
2187 unsigned int n_subdivisions[dim];
2188 for (
unsigned int i = 0; i < dim; ++i)
2189 n_subdivisions[i] = 1;
2198 const unsigned int n_subdivisions,
2205 for (
unsigned int i = 0; i < dim; ++i)
2216 const unsigned int (&n_subdivisions)[dim],
2218 const unsigned int *n_subdivisions,
2225 std::array<Tensor<1, dim>, dim>
edges;
2226 for (
unsigned int i = 0; i < dim; ++i)
2239 template <
int dim,
int spacedim>
2254 ExcMessage(
"One subdivision must be provided for each dimension."));
2256 for (
unsigned int i = 0; i < dim; ++i)
2263 "Edges in subdivided_parallelepiped() must not be degenerate."));
2297 "Edges in subdivided_parallelepiped() must point in"
2298 " different directions."));
2323 ExcInvalidInputOrientation(
2324 "The triangulation you are trying to create will consist of cells"
2325 " with negative measures. This is usually the result of input data"
2326 " that does not define a right-handed coordinate system. The usual"
2327 " fix for this is to ensure that in 1d the given point is to the"
2328 " right of the origin (or the given edge tensor is positive), in 2d"
2329 " that the two edges (and their cross product) obey the right-hand"
2330 " rule (which may usually be done by switching the order of the"
2331 " points or edge tensors), or in 3d that the edges form a"
2332 " right-handed coordinate system (which may also be accomplished by"
2333 " switching the order of the first two points or edge tensors)."));
2336 for (
unsigned int i = 0; i < dim; ++i)
2337 for (
unsigned int j = i + 1;
j < dim; ++
j)
2340 "Degenerate edges of subdivided_parallelepiped encountered."));
2343 std::vector<Point<spacedim>> points;
2363 points.push_back(
origin +
2375 for (
unsigned int i = 0; i < dim; ++i)
2377 std::vector<CellData<dim>> cells(n_cells);
2385 cells[
x].vertices[0] =
x;
2386 cells[
x].vertices[1] =
x + 1;
2389 cells[
x].material_id = 0;
2399 for (
unsigned int y = 0;
y <
n_dy; ++
y)
2400 for (
unsigned int x = 0;
x <
n_dx; ++
x)
2402 const unsigned int c =
y *
n_dx +
x;
2403 cells[c].vertices[0] =
y * (
n_dx + 1) +
x;
2404 cells[c].vertices[1] =
y * (
n_dx + 1) +
x + 1;
2405 cells[c].vertices[2] = (
y + 1) * (
n_dx + 1) +
x;
2406 cells[c].vertices[3] = (
y + 1) * (
n_dx + 1) +
x + 1;
2409 cells[c].material_id = 0;
2421 for (
unsigned int z = 0; z <
n_dz; ++z)
2422 for (
unsigned int y = 0;
y <
n_dy; ++
y)
2423 for (
unsigned int x = 0;
x <
n_dx; ++
x)
2427 cells[c].vertices[0] =
2429 cells[c].vertices[1] =
2431 cells[c].vertices[2] =
2433 cells[c].vertices[3] = z * (
n_dy + 1) * (
n_dx + 1) +
2434 (
y + 1) * (
n_dx + 1) +
x + 1;
2435 cells[c].vertices[4] =
2437 cells[c].vertices[5] = (z + 1) * (
n_dy + 1) * (
n_dx + 1) +
2439 cells[c].vertices[6] = (z + 1) * (
n_dy + 1) * (
n_dx + 1) +
2440 (
y + 1) * (
n_dx + 1) +
x;
2441 cells[c].vertices[7] = (z + 1) * (
n_dy + 1) * (
n_dx + 1) +
2442 (
y + 1) * (
n_dx + 1) +
x + 1;
2445 cells[c].material_id = 0;
2458 tria.create_triangulation(points, cells,
SubCellData());
2464 tria.begin_active(),
2466 for (; cell !=
endc; ++cell)
2468 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
2470 if (cell->face(face)->at_boundary())
2471 cell->face(face)->set_boundary_id(face);
2478 template <
int dim,
int spacedim>
2488 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2491 for (
unsigned int i = 0; i < dim; ++i)
2503 template <
int dim,
int spacedim>
2517 for (
unsigned int i = 0; i < dim; ++i)
2524 std::array<Point<spacedim>, dim> delta;
2525 for (
unsigned int i = 0; i < dim; ++i)
2533 "The first dim entries of coordinates of p1 and p2 need to be different."));
2537 std::vector<Point<spacedim>> points;
2542 points.push_back(
p1 +
x * delta[0]);
2548 points.push_back(
p1 +
x * delta[0] +
y * delta[1]);
2552 for (
unsigned int z = 0; z <=
repetitions[2]; ++z)
2555 points.push_back(
p1 +
x * delta[0] +
y * delta[1] +
2564 std::vector<CellData<dim>> cells;
2572 cells[
x].vertices[0] =
x;
2573 cells[
x].vertices[1] =
x + 1;
2574 cells[
x].material_id = 0;
2589 cells[c].vertices[3] = (
y + 1) * (
repetitions[0] + 1) +
x + 1;
2590 cells[c].material_id = 0;
2598 const unsigned int n_xy =
2608 cells[c].vertices[0] = z *
n_xy +
y *
n_x +
x;
2609 cells[c].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
2610 cells[c].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
2611 cells[c].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
2612 cells[c].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
2613 cells[c].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
2614 cells[c].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
2615 cells[c].vertices[7] =
2617 cells[c].material_id = 0;
2626 tria.create_triangulation(points, cells,
SubCellData());
2638 double epsilon = std::numeric_limits<double>::max();
2639 for (
unsigned int i = 0; i < dim; ++i)
2640 epsilon =
std::min(epsilon, 0.01 * delta[i][i]);
2643 "The distance between corner points must be positive."));
2656 const std::vector<std::vector<double>> &
step_sz,
2661 Assert(
step_sz.size() == dim, ExcInvalidRepetitionsDimension(dim));
2673 for (
unsigned int i = 0; i < dim; ++i)
2677 std::swap(
p1[i],
p2[i]);
2685 Assert(std::fabs(
x - (
p2[i] -
p1[i])) <= 1e-12 * std::fabs(
x),
2687 "The sequence of step sizes in coordinate direction " +
2689 " must be equal to the distance of the two given "
2690 "points in this coordinate direction."));
2697 std::vector<Point<dim>> points;
2703 for (
unsigned int i = 0;; ++i)
2723 for (
unsigned int j = 0;; ++
j)
2726 for (
unsigned int i = 0;; ++i)
2745 for (
unsigned int k = 0;; ++
k)
2748 for (
unsigned int j = 0;; ++
j)
2751 for (
unsigned int i = 0;; ++i)
2781 std::vector<CellData<dim>> cells;
2789 cells[
x].vertices[0] =
x;
2790 cells[
x].vertices[1] =
x + 1;
2791 cells[
x].material_id = 0;
2803 cells[c].vertices[0] =
y * (
step_sizes[0].size() + 1) +
x;
2804 cells[c].vertices[1] =
y * (
step_sizes[0].size() + 1) +
x + 1;
2805 cells[c].vertices[2] =
2807 cells[c].vertices[3] =
2809 cells[c].material_id = 0;
2817 const unsigned int n_xy =
2822 for (
unsigned int z = 0; z <
step_sizes[2].size(); ++z)
2826 const unsigned int c =
2829 cells[c].vertices[0] = z *
n_xy +
y *
n_x +
x;
2830 cells[c].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
2831 cells[c].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
2832 cells[c].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
2833 cells[c].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
2834 cells[c].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
2835 cells[c].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
2836 cells[c].vertices[7] =
2838 cells[c].material_id = 0;
2847 tria.create_triangulation(points, cells,
SubCellData());
2861 for (
unsigned int i = 1; i < dim; ++i)
2878 const std::vector<std::vector<double>> &
spacing,
2883 Assert(
spacing.size() == 1, ExcInvalidRepetitionsDimension(1));
2889 double delta = std::numeric_limits<double>::max();
2890 for (
unsigned int i = 0; i <
n_cells; ++i)
2897 std::vector<Point<1>> points;
2901 points.emplace_back(
ax);
2907 for (
unsigned int i = 0; i <
n_cells; ++i)
2912 unsigned int id = 0;
2916 cells[id].vertices[0] =
x;
2917 cells[id].vertices[1] =
x + 1;
2925 tria.create_triangulation(points, cells, t);
2936 const std::vector<std::vector<double>> &
spacing,
2941 Assert(
spacing.size() == 2, ExcInvalidRepetitionsDimension(2));
2944 double delta = std::numeric_limits<double>::max();
2945 for (
unsigned int i = 0; i < 2; ++i)
2954 ExcInvalidRepetitionsDimension(i));
2958 std::vector<Point<2>> points;
2965 points.emplace_back(
ax,
ay);
2975 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
2981 unsigned int id = 0;
2989 cells[id].vertices[3] = (
y + 1) * (
repetitions[0] + 1) +
x + 1;
2998 tria.create_triangulation(points, cells, t);
3003 double eps = 0.01 * delta;
3005 for (; cell !=
endc; ++cell)
3008 for (
const unsigned int f :
GeometryInfo<2>::face_indices())
3012 for (
unsigned int i = 0; i < 2; ++i)
3015 cell->face(f)->set_boundary_id(i * 2);
3017 cell->face(f)->set_boundary_id(i * 2 + 1);
3028 const std::vector<std::vector<double>> &
spacing,
3033 const unsigned int dim = 3;
3035 Assert(
spacing.size() == dim, ExcInvalidRepetitionsDimension(dim));
3038 double delta = std::numeric_limits<double>::max();
3039 for (
unsigned int i = 0; i < dim; ++i)
3048 ExcInvalidRepetitionsDimension(i));
3052 std::vector<Point<dim>> points;
3054 for (
unsigned int z = 0; z <=
repetitions[2]; ++z)
3062 points.emplace_back(
ax,
ay,
az);
3075 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
3082 unsigned int id = 0;
3090 cells[id].vertices[0] = z *
n_xy +
y *
n_x +
x;
3091 cells[id].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
3092 cells[id].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
3093 cells[id].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
3094 cells[id].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
3095 cells[id].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
3096 cells[id].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
3097 cells[id].vertices[7] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x + 1;
3106 tria.create_triangulation(points, cells, t);
3111 double eps = 0.01 * delta;
3114 for (; cell !=
endc; ++cell)
3121 for (
unsigned int i = 0; i < dim; ++i)
3124 cell->face(f)->set_boundary_id(i * 2);
3126 cell->face(f)->set_boundary_id(i * 2 + 1);
3133 template <
int dim,
int spacedim>
3136 const std::vector<unsigned int> &
holes)
3145 for (
unsigned int d = 0;
d < dim; ++
d)
3152 std::array<Point<spacedim>, dim> delta;
3154 for (
unsigned int i = 0; i < dim; ++i)
3157 ExcMessage(
"At least one hole needed in each direction"));
3159 delta[i][i] = (
p2[i] -
p1[i]);
3164 std::vector<Point<spacedim>> points;
3169 points.push_back(
p1 +
x * delta[0]);
3175 points.push_back(
p1 +
x * delta[0] +
y * delta[1]);
3179 for (
unsigned int z = 0; z <=
repetitions[2]; ++z)
3182 points.push_back(
p1 +
x * delta[0] +
y * delta[1] +
3192 std::vector<CellData<dim>> cells;
3202 if ((
x % 2 == 1) && (
y % 2 == 1))
3208 cells[c].vertices[3] = (
y + 1) * (
repetitions[0] + 1) +
x + 1;
3209 cells[c].material_id = 0;
3218 const unsigned int n_xy =
3229 cells[c].vertices[0] = z *
n_xy +
y *
n_x +
x;
3230 cells[c].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
3231 cells[c].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
3232 cells[c].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
3233 cells[c].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
3234 cells[c].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
3235 cells[c].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
3236 cells[c].vertices[7] =
3238 cells[c].material_id = 0;
3248 tria.create_triangulation(points, cells,
SubCellData());
3266 const unsigned int ,
3278 const unsigned int ,
3292 const double radius)
3294 return (
std::abs(p[0] - c[0]) < radius) &&
3303 template <
int dim,
int spacedim>
3307 double length = std::numeric_limits<double>::max();
3308 for (
const auto &cell :
triangulation.active_cell_iterators())
3310 length =
std::min(length, cell->line(n)->diameter());
3320 const double inner_radius,
3330 const unsigned int ,
3344 double length = std::numeric_limits<double>::max();
3345 for (
const auto &cell : tria.active_cell_iterators())
3346 for (unsigned
int n = 0; n < cell->n_lines(); ++n)
3347 length =
std::min(length, cell->line(n)->diameter());
3362 for (
const auto &cell :
cylinder_tria.active_cell_iterators())
3374 const double h) ->
void {
3378 static_cast<unsigned int>(std::round(
padding / h));
3382 for (
unsigned int i = 0; i <
num; ++i)
3386 std::vector<std::vector<double>>
step_sizes(2);
3411 for (
const auto &cell :
bulk_tria.active_cell_iterators())
3419 const double tolerance =
3429 for (
const auto &cell : tria.active_cell_iterators())
3435 for (
const unsigned int face_n : cell->face_indices())
3437 const auto &face = cell->face(
face_n);
3438 if (face->at_boundary() &&
3439 internal::point_in_2d_box(face->center(),
3455 static constexpr double tol =
3456 std::numeric_limits<double>::epsilon() * 10000;
3458 for (
const auto &cell : tria.active_cell_iterators())
3461 const auto face = cell->face(
face_n);
3462 if (face->at_boundary())
3464 const Point<2> center = face->center();
3467 face->set_boundary_id(0);
3470 face->set_boundary_id(1);
3473 face->set_boundary_id(2);
3476 face->set_boundary_id(3);
3482 face->set_boundary_id(4);
3503 const double inner_radius,
3562 ExcMessage(
"The width of the shell region must be less than 0.05 "
3563 "(and preferably close to 0.03)"));
3599 for (
const auto &cell :
bulk_tria.active_cell_iterators())
3601 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3613 2.0 * (cell->vertex(3) -
Point<2>());
3632 for (
const auto &cell :
cylinder_tria.active_cell_iterators())
3642 for (
const auto &cell :
cylinder_tria.active_cell_iterators())
3646 if (!cell->face(
face_n)->at_boundary())
3652 ExcMessage(
"If the shell region has positive width then "
3653 "there must be at least one shell."));
3698 const double shift =
3700 for (
const auto &cell : tria.active_cell_iterators())
3702 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3703 cell->vertex(v) =
Point<2>(0.2 -
shift, 0.205);
3704 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1e-10)
3705 cell->vertex(v) =
Point<2>(0.2 + shift, 0.205);
3706 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1e-10)
3707 cell->vertex(v) =
Point<2>(0.2, 0.2 - shift);
3708 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1e-10)
3709 cell->vertex(v) =
Point<2>(0.2, 0.2 + shift);
3714 for (
const auto &cell : tria.active_cell_iterators())
3720 for (
const auto &cell : tria.active_cell_iterators())
3728 for (
const auto &face : tria.active_face_iterators())
3747 *ptr +=
Point<2>(0.2, 0.2) - center;
3759 for (
const auto &face : tria.active_face_iterators())
3760 if (face->at_boundary())
3762 const Point<2> center = face->center();
3764 if (
std::abs(center[0] - 0.0) < 1e-10)
3765 face->set_boundary_id(0);
3767 else if (
std::abs(center[0] - 2.2) < 1e-10)
3768 face->set_boundary_id(1);
3771 face->set_boundary_id(2);
3776 std::abs(center[1] - 0.41) < 1.0e-10,
3778 face->set_boundary_id(3);
3820 for (
const auto &face : tria.active_face_iterators())
3822 face->set_boundary_id(3);
3827 template <
int dim,
int spacedim>
3830 const std::vector<unsigned int> &sizes,
3840 for (
unsigned int d = 0;
d < dim; ++
d)
3843 std::vector<Point<spacedim>> points;
3845 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
3848 std::vector<CellData<dim>> cells(n_cells);
3853 for (
unsigned int d = 0;
d < dim; ++
d)
3854 p[d] = 0.5 * dimensions[d] *
3857 points.push_back(p);
3858 cells[0].vertices[i] = i;
3860 cells[0].material_id = 0;
3865 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
3871 for (
unsigned int j = 0;
j < sizes[face]; ++
j, ++
cell_index)
3878 const unsigned int cellv =
3880 const unsigned int ocellv =
3892 points.push_back(p);
3897 tria.create_triangulation(points, cells,
SubCellData());
4045 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4047 std::vector<Point<2>> vertices(16);
4055 for (
const double y : 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;
4074 tria.create_triangulation(vertices,
4089 const double rl2 = (right + left) / 2;
4100 const int cell_vertices[4][4] = {{0, 1, 3, 2},
4105 for (
unsigned int i = 0; i < 4; ++i)
4107 for (
unsigned int j = 0;
j < 4; ++
j)
4108 cells[i].vertices[
j] = cell_vertices[i][
j];
4109 cells[i].material_id = 0;
4111 tria.create_triangulation(std::vector<
Point<2>>(std::begin(vertices),
4112 std::end(vertices)),
4119 cell->face(1)->set_boundary_id(1);
4121 cell->face(0)->set_boundary_id(2);
4141 const std::vector<Point<2>> vertices(std::begin(
vertices_tmp),
4146 cell_vertices[0][i] = i;
4151 cells[0].vertices[i] = cell_vertices[0][i];
4153 cells[0].material_id = 0;
4158 cell->face(0)->set_boundary_id(1);
4159 cell->face(1)->set_boundary_id(2);
4161 for (
unsigned int i = 2; i < 4; ++i)
4162 cell->face(i)->set_boundary_id(0);
4179 Point<2>((a + b) / 2, (a + b) / 2),
4183 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4187 for (
unsigned int i = 0; i < 3; ++i)
4189 for (
unsigned int j = 0;
j < 4; ++
j)
4190 cells[i].vertices[
j] = cell_vertices[i][
j];
4191 cells[i].material_id = 0;
4194 tria.create_triangulation(std::vector<
Point<2>>(std::begin(vertices),
4195 std::end(vertices)),
4203 cell->face(0)->set_boundary_id(0);
4204 cell->face(2)->set_boundary_id(1);
4207 cell->face(1)->set_boundary_id(2);
4208 cell->face(2)->set_boundary_id(1);
4209 cell->face(3)->set_boundary_id(3);
4212 cell->face(0)->set_boundary_id(0);
4213 cell->face(1)->set_boundary_id(4);
4214 cell->face(3)->set_boundary_id(5);
4220 template <
int dim,
int spacedim>
4232 for (
unsigned int d = 0;
d < dim; ++
d)
4235 ExcMessage(
"Attempting to cut away too many cells."));
4245 std::array<double, dim> h;
4247 for (
unsigned int d = 0;
d < dim; ++
d)
4267 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4269 for (
const auto &cell :
rectangle.active_cell_iterators())
4287 template <
int dim,
int spacedim>
4291 const double radius,
4294 if constexpr (dim == 2)
4298 if constexpr (spacedim == 2)
4300 else if constexpr (spacedim == 3)
4309 const double a = 1. / (1 +
std::sqrt(2.0));
4322 for (
unsigned int i = 0; i < 5; ++i)
4324 for (
unsigned int j = 0;
j < 4; ++
j)
4326 cells[i].material_id = 0;
4331 std::begin(vertices), std::end(vertices)),
4334 tria.set_all_manifold_ids_on_boundary(0);
4341 else if constexpr (dim == 3)
4347 const unsigned int n_vertices = 16;
4348 const Point<3> vertices[n_vertices] = {
4373 const unsigned int n_cells = 7;
4374 const int cell_vertices[
n_cells][8] = {
4375 {0, 1, 4, 5, 3, 2, 7, 6},
4376 {8, 9, 12, 13, 0, 1, 4, 5},
4377 {9, 13, 1, 5, 10, 14, 2, 6},
4378 {11, 10, 3, 2, 15, 14, 7, 6},
4379 {8, 0, 12, 4, 11, 3, 15, 7},
4380 {8, 9, 0, 1, 11, 10, 3, 2},
4381 {12, 4, 13, 5, 15, 7, 14, 6}};
4383 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
4385 for (
unsigned int i = 0; i <
n_cells; ++i)
4388 cells[i].vertices[
j] = cell_vertices[i][
j];
4389 cells[i].material_id = 0;
4393 tria.create_triangulation(std::vector<
Point<3>>(std::begin(vertices),
4394 std::end(vertices)),
4397 tria.set_all_manifold_ids_on_boundary(0);
4410 template <
int spacedim>
4414 const double inner_radius,
4416 const unsigned int n_cells,
4433 const unsigned int N =
4434 (
n_cells == 0 ?
static_cast<unsigned int>(
4445 std::vector<Point<spacedim>> vertices(2 * N);
4446 for (
unsigned int i = 0; i <
N; ++i)
4453 vertices[i +
N] = vertices[i] * (inner_radius /
outer_radius);
4455 vertices[i] += center;
4456 vertices[i +
N] += center;
4461 for (
unsigned int i = 0; i <
N; ++i)
4463 cells[i].vertices[0] = i;
4464 cells[i].vertices[1] = (i + 1) % N;
4465 cells[i].vertices[2] =
N + i;
4466 cells[i].vertices[3] =
N + ((i + 1) % N);
4468 cells[i].material_id = 0;
4471 tria.create_triangulation(vertices, cells,
SubCellData());
4476 tria.set_all_manifold_ids(0);
4486 const double inner_radius,
4488 const unsigned int n_cells,
4500 const double inner_radius,
4502 const unsigned int n_cells,
4515 const double inner_radius,
4517 const unsigned int n_cells)
4526 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4532 for (
const auto &face : tria.active_face_iterators())
4546 tria.set_all_manifold_ids(2);
4569 const double radius,
4581 switch (f->boundary_id())
4584 f->set_boundary_id(1);
4587 f->set_boundary_id(2);
4590 f->set_boundary_id(0);
4628 const double radius)
4630 const unsigned int dim = 2;
4644 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4648 for (
unsigned int i = 0; i < 3; ++i)
4650 for (
unsigned int j = 0;
j < 4; ++
j)
4651 cells[i].vertices[
j] = cell_vertices[i][
j];
4652 cells[i].material_id = 0;
4655 tria.create_triangulation(std::vector<
Point<dim>>(std::begin(vertices),
4656 std::end(vertices)),
4663 tria.set_all_manifold_ids_on_boundary(0);
4667 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
4669 if (cell->face(i)->boundary_id() ==
4675 if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius ||
4676 cell->face(i)->center()[1] < p[1] + 1.e-5 * radius)
4678 cell->face(i)->set_boundary_id(1);
4692 const double radius)
4697 const double a = 1. / (1 +
std::sqrt(2.0));
4708 const int cell_vertices[5][4] = {{0, 1, 2, 3},
4715 for (
unsigned int i = 0; i < 4; ++i)
4717 for (
unsigned int j = 0;
j < 4; ++
j)
4718 cells[i].vertices[
j] = cell_vertices[i][
j];
4719 cells[i].material_id = 0;
4722 tria.create_triangulation(std::vector<
Point<2>>(std::begin(vertices),
4723 std::end(vertices)),
4730 tria.set_all_manifold_ids_on_boundary(0);
4734 for (
const unsigned int i :
GeometryInfo<2>::face_indices())
4736 if (cell->face(i)->boundary_id() ==
4741 if (cell->face(i)->center()[0] < p[0] + 1.e-5 * radius)
4743 cell->face(i)->set_boundary_id(1);
4759 const double inner_radius,
4761 const unsigned int n_cells,
4777 const unsigned int N =
4778 (
n_cells == 0 ?
static_cast<unsigned int>(
4789 std::vector<Point<2>> vertices(2 * (N + 1));
4790 for (
unsigned int i = 0; i <=
N; ++i)
4802 vertices[i +
N + 1] = vertices[i] * (inner_radius /
outer_radius);
4804 vertices[i] += center;
4805 vertices[i +
N + 1] += center;
4811 for (
unsigned int i = 0; i <
N; ++i)
4813 cells[i].vertices[0] = i;
4814 cells[i].vertices[1] = (i + 1) % (N + 1);
4815 cells[i].vertices[2] =
N + 1 + i;
4816 cells[i].vertices[3] =
N + 1 + ((i + 1) % (N + 1));
4818 cells[i].material_id = 0;
4821 tria.create_triangulation(vertices, cells,
SubCellData());
4826 for (; cell != tria.end(); ++cell)
4828 cell->face(2)->set_boundary_id(1);
4830 tria.begin()->face(0)->set_boundary_id(3);
4832 tria.last()->face(1)->set_boundary_id(2);
4834 tria.set_all_manifold_ids(0);
4843 const double inner_radius,
4845 const unsigned int n_cells,
4861 const unsigned int N =
4862 (
n_cells == 0 ?
static_cast<unsigned int>(
4873 std::vector<Point<2>> vertices(2 * (N + 1));
4874 for (
unsigned int i = 0; i <=
N; ++i)
4884 vertices[i +
N + 1] = vertices[i] * (inner_radius /
outer_radius);
4886 vertices[i] += center;
4887 vertices[i +
N + 1] += center;
4893 for (
unsigned int i = 0; i <
N; ++i)
4895 cells[i].vertices[0] = i;
4896 cells[i].vertices[1] = (i + 1) % (N + 1);
4897 cells[i].vertices[2] =
N + 1 + i;
4898 cells[i].vertices[3] =
N + 1 + ((i + 1) % (N + 1));
4900 cells[i].material_id = 0;
4903 tria.create_triangulation(vertices, cells,
SubCellData());
4908 for (; cell != tria.end(); ++cell)
4910 cell->face(2)->set_boundary_id(1);
4912 tria.begin()->face(0)->set_boundary_id(3);
4914 tria.last()->face(1)->set_boundary_id(2);
4917 tria.set_all_manifold_ids(0);
4931 const double rl2 = (right + left) / 2;
4932 const double len = (right - left) / 2.;
4945 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4946 {9, 4, 2, 5, 19, 14, 12, 15},
4947 {3, 2, 7, 6, 13, 12, 17, 16},
4948 {2, 5, 6, 8, 12, 15, 16, 18}};
4950 for (
unsigned int i = 0; i < 4; ++i)
4952 for (
unsigned int j = 0;
j < 8; ++
j)
4953 cells[i].vertices[
j] = cell_vertices[i][
j];
4954 cells[i].material_id = 0;
4956 tria.create_triangulation(std::vector<
Point<3>>(std::begin(vertices),
4957 std::end(vertices)),
4964 cell->face(1)->set_boundary_id(1);
4966 cell->face(0)->set_boundary_id(2);
4982 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4984 std::vector<Point<3>> vertices(64);
4992 for (
const double z : coords)
4998 24, 26, 5, 4, 6, 1, 0,
4999 2, 9, 8, 10, 37, 36, 38,
5000 33, 32, 34, 41, 40, 42};
5002 std::vector<CellData<3>> cells(27);
5004 for (
unsigned int z = 0; z < 3; ++z)
5005 for (
unsigned int y = 0;
y < 3; ++
y)
5006 for (
unsigned int x = 0;
x < 3; ++
x)
5008 cells[
k].vertices[0] =
x + 4 *
y + 16 * z;
5009 cells[
k].vertices[1] =
x + 4 *
y + 16 * z + 1;
5010 cells[
k].vertices[2] =
x + 4 *
y + 16 * z + 4;
5011 cells[
k].vertices[3] =
x + 4 *
y + 16 * z + 5;
5012 cells[
k].vertices[4] =
x + 4 *
y + 16 * z + 16;
5013 cells[
k].vertices[5] =
x + 4 *
y + 16 * z + 17;
5014 cells[
k].vertices[6] =
x + 4 *
y + 16 * z + 20;
5015 cells[
k].vertices[7] =
x + 4 *
y + 16 * z + 21;
5020 tria.create_triangulation(vertices,
5035 ExcMessage(
"The output triangulation object needs to be empty."));
5040 const auto n_slices = 1 +
static_cast<unsigned int>(std::ceil(
5057 return Point<3>(p[0], factor * p[1], factor * p[2]);
5063 for (
const auto &face :
triangulation.active_face_iterators())
5064 if (face->at_boundary())
5067 face->set_boundary_id(1);
5070 face->set_boundary_id(2);
5072 face->set_all_manifold_ids(0);
5095 Point<3>((a + b) / 2, a, (a + b) / 2),
5102 Point<3>((a + b) / 2, (a + b) / 2, a),
5104 Point<3>(a, (a + b) / 2, (a + b) / 2),
5105 Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
5106 Point<3>(b, (a + b) / 2, (a + b) / 2),
5108 Point<3>((a + b) / 2, (a + b) / 2, b),
5116 Point<3>((a + b) / 2, b, (a + b) / 2),
5120 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5121 {1, 2, 10, 11, 4, 5, 13, 14},
5122 {3, 4, 12, 13, 6, 7, 15, 16},
5123 {4, 5, 13, 14, 7, 8, 16, 17},
5124 {9, 10, 18, 19, 12, 13, 21, 22},
5125 {10, 11, 19, 20, 13, 14, 22, 23},
5126 {12, 13, 21, 22, 15, 16, 24, 25}};
5130 for (
unsigned int i = 0; i < 7; ++i)
5132 for (
unsigned int j = 0;
j < 8; ++
j)
5133 cells[i].vertices[
j] = cell_vertices[i][
j];
5134 cells[i].material_id = 0;
5137 tria.create_triangulation(std::vector<
Point<3>>(std::begin(vertices),
5138 std::end(vertices)),
5155 ExcMessage(
"The number of rotation by pi/2 of the right square "
5156 "must be in the half-open range [0,4)."));
5158 constexpr unsigned int dim = 2;
5160 const unsigned int n_cells = 5;
5161 std::vector<CellData<dim>> cells(n_cells);
5164 const std::vector<Point<dim>> vertices = {
Point<dim>(0, 0),
5179 unsigned int cell_vertices[
n_cells][4] = {{0, 1, 2, 3},
5189 cell_vertices[1][0] = 4;
5190 cell_vertices[1][1] = 5;
5191 cell_vertices[1][2] = 1;
5192 cell_vertices[1][3] = 3;
5198 cell_vertices[1][0] = 5;
5199 cell_vertices[1][1] = 3;
5200 cell_vertices[1][2] = 4;
5201 cell_vertices[1][3] = 1;
5207 cell_vertices[1][0] = 3;
5208 cell_vertices[1][1] = 1;
5209 cell_vertices[1][2] = 5;
5210 cell_vertices[1][3] = 4;
5222 for (
const unsigned int vertex_index :
5231 tria.create_triangulation(vertices, cells,
SubCellData());
5237 const bool face_orientation,
5238 const bool face_flip,
5239 const bool face_rotation,
5242 constexpr unsigned int dim = 3;
5244 const unsigned int n_cells = 2;
5245 std::vector<CellData<dim>> cells(n_cells);
5248 const std::vector<Point<dim>> vertices = {
Point<dim>(0, 0, 0),
5261 unsigned int cell_vertices[
n_cells][8] = {
5262 {0, 1, 2, 3, 4, 5, 6, 7},
5263 {1, 8, 3, 9, 5, 10, 7, 11}};
5266 const unsigned int this_case = 4 *
static_cast<int>(face_orientation) +
5267 2 *
static_cast<int>(face_flip) +
5268 static_cast<int>(face_rotation);
5276 cell_vertices[0][0] = 1;
5277 cell_vertices[0][1] = 0;
5278 cell_vertices[0][2] = 5;
5279 cell_vertices[0][3] = 4;
5280 cell_vertices[0][4] = 3;
5281 cell_vertices[0][5] = 2;
5282 cell_vertices[0][6] = 7;
5283 cell_vertices[0][7] = 6;
5289 cell_vertices[0][0] = 5;
5290 cell_vertices[0][1] = 4;
5291 cell_vertices[0][2] = 7;
5292 cell_vertices[0][3] = 6;
5293 cell_vertices[0][4] = 1;
5294 cell_vertices[0][5] = 0;
5295 cell_vertices[0][6] = 3;
5296 cell_vertices[0][7] = 2;
5302 cell_vertices[0][0] = 7;
5303 cell_vertices[0][1] = 6;
5304 cell_vertices[0][2] = 3;
5305 cell_vertices[0][3] = 2;
5306 cell_vertices[0][4] = 5;
5307 cell_vertices[0][5] = 4;
5308 cell_vertices[0][6] = 1;
5309 cell_vertices[0][7] = 0;
5314 cell_vertices[0][0] = 3;
5315 cell_vertices[0][1] = 2;
5316 cell_vertices[0][2] = 1;
5317 cell_vertices[0][3] = 0;
5318 cell_vertices[0][4] = 7;
5319 cell_vertices[0][5] = 6;
5320 cell_vertices[0][6] = 5;
5321 cell_vertices[0][7] = 4;
5327 cell_vertices[0][0] = 0;
5328 cell_vertices[0][1] = 1;
5329 cell_vertices[0][2] = 2;
5330 cell_vertices[0][3] = 3;
5331 cell_vertices[0][4] = 4;
5332 cell_vertices[0][5] = 5;
5333 cell_vertices[0][6] = 6;
5334 cell_vertices[0][7] = 7;
5340 cell_vertices[0][0] = 2;
5341 cell_vertices[0][1] = 3;
5342 cell_vertices[0][2] = 6;
5343 cell_vertices[0][3] = 7;
5344 cell_vertices[0][4] = 0;
5345 cell_vertices[0][5] = 1;
5346 cell_vertices[0][6] = 4;
5347 cell_vertices[0][7] = 5;
5353 cell_vertices[0][0] = 6;
5354 cell_vertices[0][1] = 7;
5355 cell_vertices[0][2] = 4;
5356 cell_vertices[0][3] = 5;
5357 cell_vertices[0][4] = 2;
5358 cell_vertices[0][5] = 3;
5359 cell_vertices[0][6] = 0;
5360 cell_vertices[0][7] = 1;
5366 cell_vertices[0][0] = 4;
5367 cell_vertices[0][1] = 5;
5368 cell_vertices[0][2] = 0;
5369 cell_vertices[0][3] = 1;
5370 cell_vertices[0][4] = 6;
5371 cell_vertices[0][5] = 7;
5372 cell_vertices[0][6] = 2;
5373 cell_vertices[0][7] = 3;
5384 cell_vertices[1][0] = 8;
5385 cell_vertices[1][1] = 1;
5386 cell_vertices[1][2] = 10;
5387 cell_vertices[1][3] = 5;
5388 cell_vertices[1][4] = 9;
5389 cell_vertices[1][5] = 3;
5390 cell_vertices[1][6] = 11;
5391 cell_vertices[1][7] = 7;
5397 cell_vertices[1][0] = 10;
5398 cell_vertices[1][1] = 5;
5399 cell_vertices[1][2] = 11;
5400 cell_vertices[1][3] = 7;
5401 cell_vertices[1][4] = 8;
5402 cell_vertices[1][5] = 1;
5403 cell_vertices[1][6] = 9;
5404 cell_vertices[1][7] = 3;
5410 cell_vertices[1][0] = 11;
5411 cell_vertices[1][1] = 7;
5412 cell_vertices[1][2] = 9;
5413 cell_vertices[1][3] = 3;
5414 cell_vertices[1][4] = 10;
5415 cell_vertices[1][5] = 5;
5416 cell_vertices[1][6] = 8;
5417 cell_vertices[1][7] = 1;
5423 cell_vertices[1][0] = 9;
5424 cell_vertices[1][1] = 3;
5425 cell_vertices[1][2] = 8;
5426 cell_vertices[1][3] = 1;
5427 cell_vertices[1][4] = 11;
5428 cell_vertices[1][5] = 7;
5429 cell_vertices[1][6] = 10;
5430 cell_vertices[1][7] = 5;
5436 cell_vertices[1][0] = 1;
5437 cell_vertices[1][1] = 8;
5438 cell_vertices[1][2] = 3;
5439 cell_vertices[1][3] = 9;
5440 cell_vertices[1][4] = 5;
5441 cell_vertices[1][5] = 10;
5442 cell_vertices[1][6] = 7;
5443 cell_vertices[1][7] = 11;
5449 cell_vertices[1][0] = 5;
5450 cell_vertices[1][1] = 10;
5451 cell_vertices[1][2] = 1;
5452 cell_vertices[1][3] = 8;
5453 cell_vertices[1][4] = 7;
5454 cell_vertices[1][5] = 11;
5455 cell_vertices[1][6] = 3;
5456 cell_vertices[1][7] = 9;
5462 cell_vertices[1][0] = 7;
5463 cell_vertices[1][1] = 11;
5464 cell_vertices[1][2] = 5;
5465 cell_vertices[1][3] = 10;
5466 cell_vertices[1][4] = 3;
5467 cell_vertices[1][5] = 9;
5468 cell_vertices[1][6] = 1;
5469 cell_vertices[1][7] = 8;
5475 cell_vertices[1][0] = 3;
5476 cell_vertices[1][1] = 9;
5477 cell_vertices[1][2] = 7;
5478 cell_vertices[1][3] = 11;
5479 cell_vertices[1][4] = 1;
5480 cell_vertices[1][5] = 8;
5481 cell_vertices[1][6] = 5;
5482 cell_vertices[1][7] = 10;
5492 for (
const unsigned int vertex_index :
5501 tria.create_triangulation(vertices, cells,
SubCellData());
5506 template <
int spacedim>
5510 const double radius)
5514 const std::set<types::boundary_id> boundary_ids = {0};
5516 tria.set_all_manifold_ids(0);
5527 const double radius,
5535 std::vector<Point<3>> vertices;
5543 vertices.emplace_back(-d, height, -d);
5544 vertices.emplace_back(d, height, -d);
5545 vertices.emplace_back(-a, height, -a);
5546 vertices.emplace_back(a, height, -a);
5547 vertices.emplace_back(-a, height, a);
5548 vertices.emplace_back(a, height, a);
5549 vertices.emplace_back(-d, height, d);
5550 vertices.emplace_back(d, height, d);
5554 for (
auto &vertex : vertices)
5556 const double h = vertex[1];
5557 vertex[1] = -vertex[0];
5561 std::vector<std::vector<int>> cell_vertices;
5562 cell_vertices.push_back({0, 1, 8, 9, 2, 3, 10, 11});
5563 cell_vertices.push_back({0, 2, 8, 10, 6, 4, 14, 12});
5564 cell_vertices.push_back({2, 3, 10, 11, 4, 5, 12, 13});
5565 cell_vertices.push_back({1, 7, 9, 15, 3, 5, 11, 13});
5566 cell_vertices.push_back({6, 4, 14, 12, 7, 5, 15, 13});
5570 for (
unsigned int i = 0; i < 5; ++i)
5573 for (
unsigned int j = 0;
j < 8; ++
j)
5581 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5583 for (
unsigned int i = 0; i <
n_cells; ++i)
5585 for (
unsigned int j = 0;
j < 8; ++
j)
5586 cells[i].vertices[
j] = cell_vertices[i][
j];
5587 cells[i].material_id = 0;
5590 tria.create_triangulation(std::vector<
Point<3>>(std::begin(vertices),
5591 std::end(vertices)),
5605 tria.set_all_manifold_ids_on_boundary(0);
5611 for (
const auto &cell : tria.cell_iterators())
5613 if (cell->at_boundary(i))
5615 if (cell->face(i)->center()[0] >
half_length - tolerance)
5617 cell->face(i)->set_boundary_id(2);
5622 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5623 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5624 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5625 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5627 cell->face(i)->line(e)->set_boundary_id(2);
5628 cell->face(i)->line(e)->set_manifold_id(
5632 else if (cell->face(i)->center()[0] < -
half_length + tolerance)
5634 cell->face(i)->set_boundary_id(1);
5639 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5640 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5641 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5642 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5644 cell->face(i)->line(e)->set_boundary_id(1);
5645 cell->face(i)->line(e)->set_manifold_id(
5657 const double radius,
5667 const double radius)
5669 const unsigned int dim = 3;
5677 const double a = 0.528;
5678 const double b = 0.4533;
5679 const double c = 0.3752;
5683 center +
Point<dim>(+1, 0, 0) * (radius * a),
5684 center +
Point<dim>(0, +1, 0) * (radius * a),
5692 center +
Point<dim>(+1, +1, 1) * (radius * c),
5696 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5697 {1, 6, 2, 4, 8, 13, 9, 11},
5698 {5, 3, 6, 4, 12, 10, 13, 11},
5699 {7, 9, 10, 11, 14, 8, 12, 13}};
5703 for (
unsigned int i = 0; i < 4; ++i)
5705 for (
unsigned int j = 0;
j < 8; ++
j)
5706 cells[i].vertices[
j] = cell_vertices[i][
j];
5707 cells[i].material_id = 0;
5710 tria.create_triangulation(std::vector<
Point<dim>>(std::begin(vertices),
5711 std::end(vertices)),
5718 tria.set_all_manifold_ids_on_boundary(0);
5721 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
5723 if (cell->face(i)->boundary_id() ==
5728 if (cell->face(i)->center()[0] < center[0] + 1.e-5 * radius ||
5729 cell->face(i)->center()[1] < center[1] + 1.e-5 * radius ||
5730 cell->face(i)->center()[2] < center[2] + 1.e-5 * radius)
5732 cell->face(i)->set_boundary_id(1);
5739 const Point<3> line_vertices[2] = {
5740 cell->face(i)->line(
j)->vertex(0),
5741 cell->face(i)->line(
j)->vertex(1)};
5742 if ((std::fabs(line_vertices[0].distance(center) - radius) >
5744 (std::fabs(line_vertices[1].distance(center) - radius) >
5747 cell->face(i)->line(
j)->set_boundary_id(1);
5748 cell->face(i)->line(
j)->set_manifold_id(
5766 const double radius)
5772 const double b = a / 2.0;
5773 const double c =
d / 2.0;
5798 int cell_vertices[6][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
5799 {0, 2, 8, 10, 6, 4, 14, 12},
5800 {2, 3, 10, 11, 4, 5, 12, 13},
5801 {1, 7, 9, 15, 3, 5, 11, 13},
5802 {6, 4, 14, 12, 7, 5, 15, 13},
5803 {8, 10, 9, 11, 14, 12, 15, 13}};
5807 for (
unsigned int i = 0; i < 6; ++i)
5809 for (
unsigned int j = 0;
j < 8; ++
j)
5810 cells[i].vertices[
j] = cell_vertices[i][
j];
5811 cells[i].material_id = 0;
5814 tria.create_triangulation(std::vector<
Point<3>>(std::begin(vertices),
5815 std::end(vertices)),
5822 tria.set_all_manifold_ids_on_boundary(0);
5831 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
5833 if (!cell->at_boundary(i))
5839 if (cell->face(i)->center()[0] < center[0] + 1.e-5 * radius)
5841 cell->face(i)->set_boundary_id(1);
5846 const Point<3> line_vertices[2] = {
5847 cell->face(i)->line(
j)->vertex(0),
5848 cell->face(i)->line(
j)->vertex(1)};
5849 if ((std::fabs(line_vertices[0].distance(center) - radius) >
5851 (std::fabs(line_vertices[1].distance(center) - radius) >
5854 cell->face(i)->line(
j)->set_boundary_id(1);
5855 cell->face(i)->line(
j)->set_manifold_id(
5872 const double radius)
5892 for (
unsigned int v = 0; v <
tria_copy.n_vertices(); ++v)
5900 else if (
round == 1)
5902 for (
unsigned int v = 0; v <
tria_copy.n_vertices(); ++v)
5911 else if (
round == 2)
5912 for (
unsigned int v = 0; v <
tria_copy.n_vertices(); ++v)
5925 std::vector<CellData<dim>> cells;
5927 for (
const auto &cell :
tria_copy.cell_iterators())
5931 data.vertices[v] = cell->vertex_index(v);
5932 data.material_id = cell->material_id();
5933 data.manifold_id = cell->manifold_id();
5934 cells.push_back(
data);
5942 if (
round == dim - 1)
5951 for (
const auto &cell : tria.cell_iterators())
5952 if (cell->center().norm_square() > 0.4 * radius)
5953 cell->set_manifold_id(1);
5960 tria.set_all_manifold_ids_on_boundary(0);
5974 const double inner_radius,
5977 std::vector<Point<3>> vertices;
5978 std::vector<CellData<3>> cells;
5984 static const std::array<Point<3>, 8>
hexahedron = {{{-1, -1, -1},
5994 for (
unsigned int i = 0; i < 8; ++i)
5996 for (
unsigned int i = 0; i < 8; ++i)
5999 const unsigned int n_cells = 6;
6000 const int cell_vertices[
n_cells][8] = {
6001 {8, 9, 10, 11, 0, 1, 2, 3},
6002 {9, 11, 1, 3, 13, 15, 5, 7},
6003 {12, 13, 4, 5, 14, 15, 6, 7},
6004 {8, 0, 10, 2, 12, 4, 14, 6},
6005 {8, 9, 0, 1, 12, 13, 4, 5},
6006 {10, 2, 11, 3, 14, 6, 15, 7}};
6010 for (
unsigned int i = 0; i <
n_cells; ++i)
6013 cells[i].vertices[
j] = cell_vertices[i][
j];
6014 cells[i].material_id = 0;
6017 tria.create_triangulation(vertices, cells,
SubCellData());
6018 tria.set_all_manifold_ids(0);
6025 const double inner_radius,
6028 std::vector<Point<3>> vertices;
6029 std::vector<CellData<3>> cells;
6038 static const std::array<Point<3>, 6>
octahedron = {{{-1, 0, 0},
6046 static const std::array<Point<3>, 8>
hexahedron = {{{-1, -1, -1},
6055 for (
unsigned int i = 0; i < 8; ++i)
6057 for (
unsigned int i = 0; i < 6; ++i)
6058 vertices.push_back(p +
octahedron[i] * inner_radius);
6059 for (
unsigned int i = 0; i < 8; ++i)
6061 for (
unsigned int i = 0; i < 6; ++i)
6064 const unsigned int n_cells = 12;
6080 for (
unsigned int i = 0; i <
n_cells; ++i)
6082 for (
unsigned int j = 0;
j < 4; ++
j)
6084 cells[i].vertices[
j] =
rhombi[i][
j];
6085 cells[i].vertices[
j + 4] =
rhombi[i][
j] + 14;
6087 cells[i].material_id = 0;
6090 tria.create_triangulation(vertices, cells,
SubCellData());
6091 tria.set_all_manifold_ids(0);
6097 const unsigned int n,
6098 const unsigned int n_refinement_steps,
6100 const double inner_radius,
6129 for (
unsigned int r = 0;
r < n_refinement_steps; ++
r)
6131 tmp.refine_global(1);
6136 for (
const auto &cell : tmp.active_cell_iterators())
6140 if ((cell->vertex(v) - p).norm_square() <
6141 inner_radius * inner_radius * (1 + 1
e-12))
6148 if (
r == n_refinement_steps - 1)
6158 tmp = std::move(copy);
6159 tmp.set_all_manifold_ids(0);
6163 tria.set_all_manifold_ids(0);
6176 const double inner_radius,
6178 const unsigned int n_cells,
6184 unsigned int n_refinement_steps = 0;
6186 if (n_cells != 96 && n_cells > 12)
6189 ++n_refinement_steps;
6192 Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 ||
6193 (n_refinement_steps > 0 &&
6195 ExcMessage(
"Invalid number of coarse mesh cells"));
6197 const unsigned int n = n_refinement_steps > 0 ?
6204 internal::hyper_shell_6(tria, p, inner_radius,
outer_radius);
6207 internal::hyper_shell_12(tria, p, inner_radius,
outer_radius);
6211 internal::hyper_shell_24_48(
6212 tria, n, n_refinement_steps, p, inner_radius,
outer_radius);
6221 internal::hyper_shell_12(tmp, p, inner_radius,
outer_radius);
6222 tmp.refine_global(1);
6224 tria.set_all_manifold_ids(0);
6248 const double inner_radius,
6250 const unsigned int ,
6258 const double a = inner_radius /
std::sqrt(2.0);
6260 const double b = a / 2.0;
6261 const double c =
d / 2.0;
6263 const double hb = inner_radius *
std::sqrt(3.0) / 2.0;
6286 int cell_vertices[5][8] = {{0, 1, 8, 9, 2, 3, 10, 11},
6287 {0, 2, 8, 10, 6, 4, 14, 12},
6288 {1, 7, 9, 15, 3, 5, 11, 13},
6289 {6, 4, 14, 12, 7, 5, 15, 13},
6290 {8, 10, 9, 11, 14, 12, 15, 13}};
6294 for (
unsigned int i = 0; i < 5; ++i)
6296 for (
unsigned int j = 0;
j < 8; ++
j)
6297 cells[i].vertices[
j] = cell_vertices[i][
j];
6298 cells[i].material_id = 0;
6301 tria.create_triangulation(std::vector<
Point<3>>(std::begin(vertices),
6302 std::end(vertices)),
6312 for (; cell != tria.end(); ++cell)
6313 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6314 if (cell->at_boundary(i))
6315 cell->face(i)->set_all_boundary_ids(2);
6321 for (cell = tria.begin(); cell != tria.end(); ++cell)
6322 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6323 if (cell->at_boundary(i))
6333 face->set_all_boundary_ids(0);
6335 face->set_all_boundary_ids(1);
6339 tria.set_all_manifold_ids(0);
6349 const double inner_radius,
6351 const unsigned int n,
6356 if (n == 0 || n == 3)
6358 const double a = inner_radius *
std::sqrt(2.0) / 2
e0;
6363 const double h = inner_radius / 2
e0;
6365 std::vector<Point<3>> vertices;
6367 vertices.push_back(center +
Point<3>(0, inner_radius, 0));
6368 vertices.push_back(center +
Point<3>(a, a, 0));
6369 vertices.push_back(center +
Point<3>(b, b, 0));
6371 vertices.push_back(center +
Point<3>(0, a, a));
6372 vertices.push_back(center +
Point<3>(c, c, h));
6373 vertices.push_back(center +
Point<3>(d, d, e));
6374 vertices.push_back(center +
Point<3>(0, b, b));
6375 vertices.push_back(center +
Point<3>(inner_radius, 0, 0));
6377 vertices.push_back(center +
Point<3>(a, 0, a));
6378 vertices.push_back(center +
Point<3>(b, 0, b));
6379 vertices.push_back(center +
Point<3>(0, 0, inner_radius));
6382 const int cell_vertices[3][8] = {
6383 {0, 1, 3, 2, 4, 5, 7, 6},
6384 {1, 8, 2, 9, 5, 10, 6, 11},
6385 {4, 5, 7, 6, 12, 10, 13, 11},
6387 std::vector<CellData<3>> cells(3);
6389 for (
unsigned int i = 0; i < 3; ++i)
6391 for (
unsigned int j = 0;
j < 8; ++
j)
6392 cells[i].vertices[
j] = cell_vertices[i][
j];
6393 cells[i].material_id = 0;
6396 tria.create_triangulation(vertices,
6408 tria.set_all_manifold_ids(0);
6417 const double length,
6418 const double inner_radius,
6438 const unsigned int N_r =
6443 const unsigned int N_z =
6445 static_cast<unsigned int>(std::ceil(
6456 for (
unsigned int i = 0; i <
N_r; ++i)
6466 for (
unsigned int j = 0;
j <=
N_z; ++
j)
6467 for (
unsigned int i = 0; i < 2 *
N_r; ++i)
6477 for (
unsigned int j = 0;
j <
N_z; ++
j)
6478 for (
unsigned int i = 0; i <
N_r; ++i)
6480 cells[i +
j *
N_r].vertices[0] = i + (
j + 1) * 2 *
N_r;
6481 cells[i +
j *
N_r].vertices[1] = (i + 1) %
N_r + (
j + 1) * 2 *
N_r;
6482 cells[i +
j *
N_r].vertices[2] = i +
j * 2 *
N_r;
6483 cells[i +
j *
N_r].vertices[3] = (i + 1) %
N_r +
j * 2 *
N_r;
6485 cells[i +
j *
N_r].vertices[4] =
N_r + i + (
j + 1) * 2 *
N_r;
6486 cells[i +
j *
N_r].vertices[5] =
6488 cells[i +
j *
N_r].vertices[6] =
N_r + i +
j * 2 *
N_r;
6489 cells[i +
j *
N_r].vertices[7] =
N_r + ((i + 1) %
N_r) +
j * 2 *
N_r;
6491 cells[i +
j *
N_r].material_id = 0;
6495 tria.set_all_manifold_ids(0);
6510 double eps_z = 1
e-6 * length;
6521 for (
const auto &cell : tria.active_cell_iterators())
6524 if (!cell->face(f)->at_boundary())
6530 if ((std::fabs(z) >
eps_z) &&
6531 (std::fabs(z - length) >
eps_z))
6542 for (
const auto &cell : tria.active_cell_iterators())
6545 if (cell->face(f)->at_boundary())
6554 if (std::fabs(z) <
eps_z)
6556 cell->face(f)->set_boundary_id(2);
6558 else if (std::fabs(z - length) <
6561 cell->face(f)->set_boundary_id(3);
6566 cell->face(f)->set_boundary_id(1);
6571 cell->face(f)->set_boundary_id(0);
6580 template <
int dim,
int spacedim>
6589 std::vector<Point<spacedim>> vertices;
6590 std::vector<CellData<dim>> cells;
6597 ExcMessage(
"The input triangulations must be non-empty "
6598 "and must not be refined."));
6610 vertices.insert(vertices.end(),
6615 for (
unsigned int &
vertex_n : cell_data.vertices)
6617 cells.push_back(cell_data);
6659 if (std::all_of(cells.begin(), cells.end(), [](
const auto &cell) {
6660 return cell.vertices.size() ==
6661 ReferenceCells::get_hypercube<dim>().n_vertices();
6675 for (
const auto &cell : tria->cell_iterators())
6677 for (
const auto &f : cell->face_indices())
6693 template <
int dim,
int spacedim>
6746 template <
int structdim>
6750 static_assert(structdim == 1 || structdim == 2,
6751 "This function is only implemented for lines and "
6759 std::sort(std::begin(cell_data.vertices),
6760 std::end(cell_data.vertices));
6761 else if (structdim == 2)
6765 std::copy(std::begin(cell_data.vertices),
6766 std::end(cell_data.vertices),
6796 std::begin(cell_data.vertices));
6803 return std::lexicographical_compare(std::begin(a.vertices),
6804 std::end(a.vertices),
6805 std::begin(
b.vertices),
6806 std::end(
b.vertices));
6817 std::upper_bound(left,
subcell_data.end(), *left, compare);
6820 if (left + 1 != right)
6821 for (
auto it = left;
it != right; ++
it)
6824 Assert(
it->manifold_id == left->manifold_id,
6826 "In the process of grid generation a single "
6827 "line or quadrilateral has been assigned two "
6828 "different manifold ids. This can happen when "
6829 "a Triangulation is copied, e.g., via "
6830 "GridGenerator::replicate_triangulation() and "
6831 "not all external boundary faces have the same "
6832 "manifold id. Double check that all faces "
6833 "which you expect to be merged together have "
6834 "the same manifold id."));
6846 template <
int dim,
int spacedim>
6849 const std::vector<unsigned int> &
extents,
6856 ExcMessage(
"The Triangulation must be copied at least one time in "
6857 "each coordinate dimension."));
6860 const auto &
min =
bbox.get_boundary_points().first;
6861 const auto &
max =
bbox.get_boundary_points().second;
6863 std::array<Tensor<1, spacedim>, dim>
offsets;
6864 for (
unsigned int d = 0;
d < dim; ++
d)
6865 offsets[d][d] = max[d] - min[d];
6869 for (
unsigned int d = 0;
d < dim; ++
d)
6896 for (
unsigned int &vertex :
6904 for (
unsigned int &vertex :
6917 1e-6 * input.begin_active()->diameter());
6947 template <
int dim,
int spacedim>
6955 ExcMessage(
"The two input triangulations are not derived from "
6956 "the same coarse mesh as required."));
6963 ExcMessage(
"The source triangulations for this function must both "
6964 "be available entirely locally, and not be distributed "
6965 "triangulations."));
6999 result.execute_coarsening_and_refinement();
7005 template <
int dim,
int spacedim>
7019 std::vector<CellData<dim>> cells;
7023 Assert(
static_cast<unsigned int>(cell->level()) ==
7026 "Your input triangulation appears to have "
7027 "adaptively refined cells. This is not allowed. You can "
7028 "only call this function on a triangulation in which "
7029 "all cells are on the same refinement level."));
7033 this_cell.vertices[v] = cell->vertex_index(v);
7034 this_cell.material_id = cell->material_id();
7058 const double height,
7063 Assert(input.n_levels() == 1,
7065 "The input triangulation must be a coarse mesh, i.e., it must "
7066 "not have been refined."));
7068 ExcMessage(
"The output triangulation object needs to be empty."));
7070 ExcMessage(
"The given height for extrusion must be positive."));
7073 "The number of slices for extrusion must be at least 2."));
7077 for (
unsigned int i = 0; i <
n_slices; ++i)
7089 const double height,
7103 "GridTools::extrude_triangulation() is only available "
7104 "for Triangulation<3, 3> as output triangulation."));
7117 Assert(input.n_levels() == 1,
7119 "The input triangulation must be a coarse mesh, i.e., it must "
7120 "not have been refined."));
7122 ExcMessage(
"The output triangulation object needs to be empty."));
7125 "The number of slices for extrusion must be at least 2."));
7127 ExcMessage(
"Slice z-coordinates should be in ascending order"));
7128 Assert(input.all_reference_cells_are_hyper_cube(),
7130 "This function is only implemented for quadrilateral meshes."));
7132 const auto priorities = [&]() -> std::vector<types::manifold_id> {
7148 "The given vector of manifold ids may not contain any "
7149 "duplicated entries."));
7151 input.get_manifold_ids();
7156 message <<
"The given triangulation has manifold ids {";
7161 <<
" the given vector of manifold ids is {";
7167 <<
" These vectors should contain the same elements.\n";
7168 const std::string m =
message.str();
7177 input.get_manifold_ids();
7182 return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
7183 &input.get_manifold(id)) == nullptr;
7192 std::vector<Point<3>> points(
n_slices * input.n_vertices());
7193 std::vector<CellData<3>> cells;
7194 cells.reserve((
n_slices - 1) * input.n_active_cells());
7211 for (
const auto &cell : input.active_cell_iterators())
7224 (
slice_n + 1) * input.n_vertices();
7227 this_cell.material_id = cell->material_id();
7229 this_cell.manifold_id = cell->manifold_id();
7237 std::vector<CellData<2>> &quads =
subcell_data.boundary_quads;
7239 quads.reserve(input.n_active_lines() * (
n_slices - 1) +
7240 input.n_active_cells() * 2);
7241 for (
const auto &face : input.active_face_iterators())
7244 quad.boundary_id = face->boundary_id();
7245 if (face->at_boundary())
7248 quad.manifold_id = face->manifold_id();
7252 face->vertex_index(0) +
slice_n * input.n_vertices();
7254 face->vertex_index(1) +
slice_n * input.n_vertices();
7256 face->vertex_index(0) + (
slice_n + 1) * input.n_vertices();
7258 face->vertex_index(1) + (
slice_n + 1) * input.n_vertices();
7259 quads.push_back(quad);
7266 for (
const auto &cell : input.active_cell_iterators())
7270 quad.manifold_id = cell->manifold_id();
7274 cell->vertex_index(0) +
slice_n * input.n_vertices();
7276 cell->vertex_index(1) +
slice_n * input.n_vertices();
7278 cell->vertex_index(2) +
slice_n * input.n_vertices();
7280 cell->vertex_index(3) +
slice_n * input.n_vertices();
7281 quads.push_back(quad);
7292 "The input triangulation to this function is using boundary "
7293 "indicators in a range that do not allow using "
7294 "max_boundary_id+1 and max_boundary_id+2 as boundary "
7295 "indicators for the bottom and top faces of the "
7296 "extruded triangulation."));
7299 for (
const auto &cell : input.active_cell_iterators())
7303 quad.vertices[0] = cell->vertex_index(0);
7304 quad.vertices[1] = cell->vertex_index(1);
7305 quad.vertices[2] = cell->vertex_index(2);
7306 quad.vertices[3] = cell->vertex_index(3);
7308 quad.manifold_id = cell->manifold_id();
7309 quads.push_back(quad);
7312 for (
unsigned int &vertex : quad.vertices)
7313 vertex += (
n_slices - 1) * input.n_vertices();
7315 quad.manifold_id = cell->manifold_id();
7316 quads.push_back(quad);
7332 for (
const auto &face :
result.active_face_iterators())
7358 "GridTools::extrude_triangulation() is only available "
7359 "for Triangulation<3, 3> as output triangulation."));
7378 template <
int spacedim>
7381 const double inner_radius,
7390 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7400 for (
const auto &cell :
triangulation.active_cell_iterators())
7403 if (cell->face(f)->at_boundary())
7412 std::
fabs(cell->face(f)->vertex(v)[1])) <
7414 cell->face(f)->vertex(v) *=
std::
sqrt(2.);
7417 for (
const auto &cell :
triangulation.active_cell_iterators())
7419 for (
const unsigned int f : cell->face_indices())
7420 if (cell->face(f)->at_boundary())
7422 const double dx = cell->face(f)->center()[0] - center[0];
7423 const double dy = cell->face(f)->center()[1] - center[1];
7427 cell->face(f)->set_boundary_id(0);
7429 cell->face(f)->set_boundary_id(1);
7431 cell->face(f)->set_boundary_id(2);
7433 cell->face(f)->set_boundary_id(3);
7436 cell->face(f)->set_boundary_id(4);
7437 cell->face(f)->set_manifold_id(0);
7442 const double d = (cell->face(f)->center() - center).
norm();
7443 if (d - inner_radius < 0)
7445 cell->face(f)->set_boundary_id(1);
7446 cell->face(f)->set_manifold_id(0);
7449 cell->face(f)->set_boundary_id(0);
7461 const double inner_radius,
7480 const double inner_radius,
7500 const double inner_radius,
7504 const unsigned int n_cells,
7511 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7515 std::vector<double> radii;
7516 radii.push_back(inner_radius);
7520 radii.push_back(inner_radius +
7539 n_cells == 0 ? (dim == 2 ? 8 : 12) :
7566 100.0 * std::numeric_limits<double>::epsilon();
7570 const double radius) {
7585 for (
const auto &cell :
triangulation.active_cell_iterators())
7588 auto face = cell->face(
face_n);
7589 if (face->at_boundary())
7591 if (((face->vertex(0) - center).norm() - inner_radius) <
7596 face->set_all_boundary_ids(0);
7602 face->set_all_boundary_ids(1);
7613 const double inner_radius,
7616 const unsigned int Nz,
7622 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7637 for (
const auto &cell :
triangulation.active_cell_iterators())
7639 for (
const auto f : cell->face_indices())
7640 if (cell->face(f)->at_boundary())
7644 const unsigned int vv = cell->face(f)->vertex_index(v);
7657 cell->face(f)->vertex(v);
7665 cell->face(f)->vertex(v) =
7674 for (
const auto &cell :
triangulation.active_cell_iterators())
7676 for (
const unsigned int f : cell->face_indices())
7677 if (cell->face(f)->at_boundary())
7679 const double dx = cell->face(f)->center()[0];
7680 const double dy = cell->face(f)->center()[1];
7681 const double dz = cell->face(f)->center()[2];
7686 cell->face(f)->set_boundary_id(0);
7689 cell->face(f)->set_boundary_id(1);
7692 cell->face(f)->set_boundary_id(2);
7695 cell->face(f)->set_boundary_id(3);
7698 cell->face(f)->set_boundary_id(4);
7701 cell->face(f)->set_boundary_id(5);
7705 cell->face(f)->set_all_boundary_ids(6);
7706 cell->face(f)->set_all_manifold_ids(0);
7713 const double d = c.
norm();
7714 if (d - inner_radius < 0)
7716 cell->face(f)->set_all_boundary_ids(1);
7717 cell->face(f)->set_all_manifold_ids(0);
7720 cell->face(f)->set_boundary_id(0);
7729 template <
int dim,
int spacedim1,
int spacedim2>
7738 "This function cannot be used on "
7739 "parallel::distributed::Triangulation objects as inputs."));
7741 ExcMessage(
"This function does not work for meshes that have "
7750 std::vector<Point<spacedim2>> v(
in_vertices.size());
7751 for (
unsigned int i = 0; i <
in_vertices.size(); ++i)
7752 for (
unsigned int d = 0;
d < spacedim; ++
d)
7755 std::vector<CellData<dim>> cells(
in_tria.n_active_cells());
7756 for (
const auto &cell :
in_tria.active_cell_iterators())
7758 const unsigned int id = cell->active_cell_index();
7760 cells[id].vertices.resize(cell->n_vertices());
7762 cells[id].vertices[i] = cell->vertex_index(i);
7763 cells[id].material_id = cell->material_id();
7764 cells[id].manifold_id = cell->manifold_id();
7781 .clear_user_flags_line();
7785 for (
const auto &face :
in_tria.active_face_iterators())
7787 if (face->at_boundary())
7811 if ((face->user_flag_set() ==
false) &&
7826 face->set_user_flag();
7842 .clear_user_flags_line();
7847 .clear_user_flags_quad();
7851 for (
const auto &face :
in_tria.active_face_iterators())
7853 if (face->at_boundary())
7871 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7872 if (face->line(e)->user_flag_set() ==
false)
7876 edge = face->line(e);
7888 edge->set_user_flag();
7906 if (face->user_flag_set() ==
false)
7922 face->set_user_flag();
7931 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7932 if (face->line(e)->at_boundary() ==
false)
7933 if (face->line(e)->user_flag_set() ==
false)
7942 edge->vertex_index(i);
7950 edge->set_user_flag();
7969 for (
const auto i :
out_tria.get_manifold_ids())
7976 template <
int dim,
int spacedim>
7983 in_tria.all_reference_cells_are_hyper_cube(),
7985 "GridGenerator::convert_hypercube_to_simplex_mesh() expects a mesh that consists only of quads/hexes."));
7988 if (
in_tria.n_global_levels() > 1)
8017 {{{0, 1, 12, 10}}, {{2, 3, 11, 12}}, {{7, 6, 11, 13}},
8018 {{5, 4, 13, 10}}, {{0, 2, 8, 12}}, {{4, 6, 13, 8}},
8019 {{5, 13, 7, 9}}, {{1, 9, 3, 12}}, {{0, 8, 4, 10}},
8020 {{1, 5, 9, 10}}, {{3, 7, 11, 9}}, {{2, 6, 8, 11}},
8021 {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
8022 {{12, 13, 8, 10}}, {{13, 8, 10, 4}}, {{13, 10, 9, 5}},
8023 {{13, 9, 11, 7}}, {{13, 11, 8, 6}}, {{10, 12, 9, 1}},
8024 {{9, 12, 11, 3}}, {{11, 12, 8, 2}}, {{8, 12, 10, 0}}}};
8032 {{{{1, 5}}, {{5, 3}}}},
8033 {{{{0, 6}}, {{6, 1}}}},
8034 {{{{2, 7}}, {{7, 3}}}}}};
8042 {{{{{0, 4, 8}}, {{4, 8, 6}}, {{8, 6, 2}}, {{0, 2, 8}}}},
8043 {{{{1, 3, 9}}, {{3, 9, 7}}, {{9, 7, 5}}, {{1, 9, 5}}}},
8044 {{{{0, 1, 10}}, {{1, 10, 5}}, {{10, 5, 4}}, {{0, 10, 4}}}},
8045 {{{{2, 3, 11}}, {{3, 11, 7}}, {{11, 7, 6}}, {{2, 11, 6}}}},
8046 {{{{0, 1, 12}}, {{1, 12, 3}}, {{12, 3, 2}}, {{0, 12, 2}}}},
8047 {{{{4, 5, 13}}, {{5, 13, 7}}, {{13, 7, 6}}, {{4, 13, 6}}}}}};
8066 {{{0, 12, 10}}, {{12, 1, 10}}, {{12, 1, 9}}, {{12, 3, 9}},
8067 {{12, 2, 11}}, {{12, 3, 11}}, {{12, 0, 8}}, {{12, 2, 8}},
8068 {{9, 13, 5}}, {{13, 7, 9}}, {{11, 7, 13}}, {{11, 6, 13}},
8069 {{4, 8, 13}}, {{6, 8, 13}}, {{4, 13, 10}}, {{13, 5, 10}},
8070 {{10, 9, 5}}, {{10, 9, 1}}, {{11, 9, 7}}, {{11, 9, 3}},
8071 {{8, 11, 2}}, {{8, 11, 6}}, {{8, 10, 0}}, {{8, 10, 4}},
8072 {{12, 3, 9}}, {{12, 9, 11}}, {{12, 3, 11}}, {{3, 9, 11}},
8073 {{2, 12, 8}}, {{2, 12, 11}}, {{2, 11, 8}}, {{8, 12, 11}},
8074 {{0, 12, 10}}, {{0, 12, 8}}, {{0, 8, 10}}, {{8, 10, 12}},
8075 {{12, 1, 10}}, {{12, 1, 9}}, {{1, 10, 9}}, {{10, 9, 12}},
8076 {{10, 8, 4}}, {{10, 8, 13}}, {{4, 13, 8}}, {{4, 13, 10}},
8077 {{10, 9, 13}}, {{10, 9, 5}}, {{13, 5, 10}}, {{13, 5, 9}},
8078 {{13, 7, 9}}, {{13, 7, 11}}, {{9, 11, 13}}, {{9, 11, 7}},
8079 {{8, 11, 13}}, {{8, 11, 6}}, {{6, 13, 8}}, {{6, 13, 11}},
8080 {{12, 13, 10}}, {{12, 13, 8}}, {{8, 10, 13}}, {{8, 10, 12}},
8081 {{12, 13, 10}}, {{12, 13, 9}}, {{10, 9, 13}}, {{10, 9, 12}},
8082 {{12, 13, 9}}, {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
8083 {{12, 13, 11}}, {{12, 13, 8}}, {{8, 11, 13}}, {{8, 11, 12}}}};
8089 {{{12, 10}}, {{12, 9}}, {{12, 11}}, {{12, 8}}, {{9, 13}}, {{11, 13}},
8090 {{8, 13}}, {{10, 13}}, {{10, 9}}, {{9, 11}}, {{11, 8}}, {{8, 10}},
8091 {{12, 9}}, {{12, 11}}, {{11, 9}}, {{12, 8}}, {{12, 11}}, {{11, 8}},
8092 {{12, 8}}, {{12, 10}}, {{10, 8}}, {{12, 10}}, {{12, 9}}, {{9, 10}},
8093 {{13, 10}}, {{13, 8}}, {{8, 10}}, {{13, 10}}, {{13, 9}}, {{9, 10}},
8094 {{13, 11}}, {{13, 9}}, {{11, 9}}, {{13, 11}}, {{13, 8}}, {{11, 8}},
8095 {{12, 13}}, {{8, 10}}, {{8, 13}}, {{10, 13}}, {{8, 12}}, {{10, 12}},
8096 {{12, 13}}, {{10, 9}}, {{10, 13}}, {{9, 13}}, {{10, 12}}, {{9, 12}},
8097 {{12, 13}}, {{9, 11}}, {{9, 13}}, {{11, 13}}, {{9, 12}}, {{11, 12}},
8098 {{12, 13}}, {{11, 8}}, {{11, 13}}, {{8, 13}}, {{11, 12}}, {{8, 12}}}};
8111 {{{{{4, 8}}, {{6, 8}}, {{0, 8}}, {{2, 8}}}},
8112 {{{{5, 9}}, {{7, 9}}, {{1, 9}}, {{3, 9}}}},
8113 {{{{4, 10}}, {{5, 10}}, {{0, 10}}, {{1, 10}}}},
8114 {{{{6, 11}}, {{7, 11}}, {{2, 11}}, {{3, 11}}}},
8115 {{{{2, 12}}, {{3, 12}}, {{0, 12}}, {{1, 12}}}},
8116 {{{{6, 13}}, {{7, 13}}, {{4, 13}}, {{5, 13}}}}}};
8118 std::vector<Point<spacedim>> vertices;
8119 std::vector<CellData<dim>> cells;
8136 for (
const auto &cell :
ref_tria.cell_iterators())
8145 const auto v_global = cell->vertex_index(v);
8151 vertices.push_back(cell->vertex(v));
8159 for (
const auto f : cell->face_indices())
8161 const auto f_global = cell->face_index(f);
8168 cell->face(f)->center(
true));
8184 vertices.push_back(cell->center(
true));
8204 cell_data.material_id =
8211 cell_data.vertices[i] =
8214 cells.push_back(cell_data);
8316 for (
const auto f : cell->face_indices())
8318 const auto bid = cell->face(f)->boundary_id();
8319 const auto mid = cell->face(f)->manifold_id();
8349 for (
const auto e : cell->line_indices())
8351 auto edge = cell->line(e);
8371 for (
const auto i :
out_tria.get_manifold_ids())
8378 template <
int spacedim>
8389 template <
int dim,
int spacedim>
8395 if (
in_tria.n_global_levels() > 1)
8405 {{{0, 1, 3}}, {{1, 2, 3}}, {{2, 0, 3}}}};
8411 {{{{{0, 1}}}}, {{{{1, 2}}}}, {{{{2, 0}}}}}};
8415 {{{0, 1, 2, 4}}, {{1, 0, 3, 4}}, {{0, 2, 3, 4}}, {{2, 1, 3, 4}}}};
8421 {{{{{0, 1, 2}}}}, {{{{1, 0, 3}}}}, {{{{0, 2, 3}}}}, {{{{2, 1, 3}}}}}};
8434 std::vector<Point<spacedim>> vertices;
8435 std::vector<CellData<dim>> cells;
8448 for (
const auto &cell :
ref_tria.cell_iterators())
8451 cell->reference_cell().is_simplex(),
8453 "Cell with invalid ReferenceCell encountered. GridGenerator::alfeld_split_of_simplex_mesh() "
8454 "only supports simplex meshes as input."));
8464 const auto v_global = cell->vertex_index(v);
8470 vertices.push_back(cell->vertex(v));
8481 vertices.push_back(barycenter /
static_cast<double>(dim + 1));
8504 cell_data.vertices[i] =
8507 cell_data.material_id =
8510 cells.push_back(cell_data);
8581 for (
const auto f : cell->face_indices())
8583 const auto bid = cell->face(f)->boundary_id();
8584 const auto mid = cell->face(f)->manifold_id();
8606 for (
const auto l : cell->line_indices())
8608 const auto bid = cell->line(l)->boundary_id();
8609 const auto mid = cell->line(l)->manifold_id();
8611 for (
const auto &line_vertices :
8621 for (
const auto i :
out_tria.get_manifold_ids())
8628 template <
template <
int,
int>
class MeshType,
int dim,
int spacedim>
8632 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8639 const std::set<types::boundary_id> &boundary_ids)
8659 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>
8669 std::vector<CellData<boundary_dim>> cells;
8671 std::vector<Point<spacedim>> vertices;
8713 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
8718 if (face->at_boundary() &&
8719 (boundary_ids.empty() ||
8720 (boundary_ids.find(face->boundary_id()) != boundary_ids.end())))
8724 for (
const unsigned int j :
8727 const unsigned int v_index = face->vertex_index(
j);
8731 vertices.push_back(face->vertex(
j));
8740 c_data.manifold_id = face->manifold_id();
8750 for (
unsigned int e = 0;
e < 4; ++
e)
8783 edge.boundary_id = 0;
8784 edge.manifold_id = face->line(e)->manifold_id();
8803 for (
const auto &cell :
surface_mesh.active_cell_iterators())
8804 for (unsigned
int vertex = 0; vertex < 2; ++vertex)
8805 if (cell->face(vertex)->at_boundary())
8806 cell->face(vertex)->set_boundary_id(0);
8814 std::vector<std::pair<
8815 const typename MeshType<dim - 1, spacedim>::cell_iterator,
8816 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>>
8818 for (
const auto &cell :
surface_mesh.active_cell_iterators())
8849 .
second.first->has_children())
8854 .
second.first->refinement_case() ==
8858 .first->set_refine_flag();
8871 .execute_coarsening_and_refinement();
8877 const typename MeshType<dim - 1, spacedim>::cell_iterator
8886 for (
unsigned int child_n = 0;
8906 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8916 for (
const auto i :
volume_mesh.get_triangulation().get_manifold_ids())
8929 template <
int dim,
int spacedim>
8940 std::vector<Point<spacedim>> vertices;
8941 std::vector<CellData<dim>> cells;
8951 for (
unsigned int i = 0; i <=
repetitions[0]; ++i)
8960 std::array<unsigned int, 4> quad{{
8970 tri.vertices = {quad[0], quad[1], quad[2]};
8971 cells.push_back(
tri);
8977 tri.vertices = {quad[3], quad[2], quad[1]};
8978 cells.push_back(
tri);
8992 for (
unsigned int i = 0; i <=
repetitions[0]; ++i)
8995 p1[2] + dx[2] *
k));
9003 std::array<unsigned int, 8> quad{
9024 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
9025 cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}};
9027 cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}};
9029 cells.push_back(cell);
9035 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
9036 cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}};
9038 cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}};
9039 cells.push_back(cell);
9045 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
9046 cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}};
9048 cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}};
9049 cells.push_back(cell);
9055 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
9056 cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}};
9058 cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}};
9059 cells.push_back(cell);
9065 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
9066 cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}};
9068 cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}};
9069 cells.push_back(cell);
9079 tria.create_triangulation(vertices, cells,
SubCellData());
9091 double epsilon = std::numeric_limits<double>::max();
9092 for (
unsigned int i = 0; i < dim; ++i)
9097 "The distance between corner points must be positive."));
9107 template <
int dim,
int spacedim>
9137# include "grid/grid_generator.inst"
void add_parameter(const std::string &entry, ParameterType ¶meter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
numbers::NumberTraits< Number >::real_type norm() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcLowerRange(int arg1, int arg2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
void copy_boundary_to_manifold_id(Triangulation< dim, spacedim > &tria, const bool reset_boundary_ids=false)
void consistently_order_cells(std::vector< CellData< dim > > &cells)
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
std::vector< index_type > data
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
Expression fabs(const Expression &x)
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void subdivided_hyper_cube_with_simplices(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double p1=0.0, const double p2=1.0, const bool colorize=false)
void parallelepiped(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void hyper_cross(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &sizes, const bool colorize_cells=false)
A center cell with stacks of cell protruding from each surface.
void hyper_cube_with_cylindrical_hole(Triangulation< dim, spacedim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
void hyper_ball_balanced(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > ¢er=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
void enclosed_hyper_cube(Triangulation< dim > &tria, const double left=0., const double right=1., const double thickness=1., const bool colorize=false)
void replicate_triangulation(const Triangulation< dim, spacedim > &input, const std::vector< unsigned int > &extents, Triangulation< dim, spacedim > &result)
Replicate a given triangulation in multiple coordinate axes.
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void general_cell(Triangulation< dim, spacedim > &tria, const std::vector< Point< spacedim > > &vertices, const bool colorize=false)
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void hyper_shell(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)
void eccentric_hyper_shell(Triangulation< dim > &triangulation, const Point< dim > &inner_center, const Point< dim > &outer_center, const double inner_radius, const double outer_radius, const unsigned int n_cells)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void moebius(Triangulation< 3, 3 > &tria, const unsigned int n_cells, const unsigned int n_rotations, const double R, const double r)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void half_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void quarter_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
Rectangular domain with rectangular pattern of holes.
void create_union_triangulation(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result)
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > ¢er={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void subdivided_hyper_rectangle_with_simplices(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void non_standard_orientation_mesh(Triangulation< 2 > &tria, const unsigned int n_rotate_middle_square)
return_type extract_boundary_mesh(const MeshType< dim, spacedim > &volume_mesh, MeshType< dim - 1, spacedim > &surface_mesh, const std::set< types::boundary_id > &boundary_ids=std::set< types::boundary_id >())
void alfeld_split_of_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_parallelepiped(Triangulation< dim > &tria, const unsigned int n_subdivisions, const Point< dim >(&corners)[dim], const bool colorize=false)
void subdivided_cylinder(Triangulation< dim > &tria, const unsigned int x_subdivisions, const double radius=1., const double half_length=1.)
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
void subdivided_hyper_L(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &bottom_left, const Point< dim > &top_right, const std::vector< int > &n_cells_to_remove)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=Point< spacedim >(), const double radius=1.)
void concentric_hyper_shells(Triangulation< dim > &triangulation, const Point< dim > ¢er, const double inner_radius=0.125, const double outer_radius=0.25, const unsigned int n_shells=1, const double skewness=0.1, const unsigned int n_cells_per_shell=0, const bool colorize=false)
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
void truncated_cone(Triangulation< dim > &tria, const double radius_0=1.0, const double radius_1=0.5, const double half_length=1.0)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void half_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cylinder_shell(Triangulation< dim > &tria, const double length, const double inner_radius, const double outer_radius, const unsigned int n_radial_cells=0, const unsigned int n_axial_cells=0, const bool colorize=false)
void flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
constexpr T fixed_power(const T t)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
long double gamma(const unsigned int n)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy(const T *begin, const T *end, U *dest)
constexpr unsigned int invalid_unsigned_int
constexpr types::boundary_id internal_face_boundary_id
constexpr types::boundary_id invalid_boundary_id
constexpr types::manifold_id flat_manifold_id
constexpr types::material_id invalid_material_id
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
inline ::VectorizedArray< Number, width > tanh(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
inline ::VectorizedArray< Number, width > atan(const ::VectorizedArray< Number, width > &x)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
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)
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)