52 : airfoil_type(
"NACA")
54 , joukowski_center(-0.1, 0.14)
60 , incline_factor(0.35)
63 , n_subdivision_x_0(3)
64 , n_subdivision_x_1(2)
65 , n_subdivision_x_2(5)
67 , airfoil_sampling_factor(2)
70 airfoil_length <= height,
72 "Mesh is to small to enclose airfoil! Choose larger field or smaller"
75 ExcMessage(
"incline_factor has to be in [0,1)!"));
88 "Mesh height measured from airfoil nose to horizontal boundaries");
92 "Length measured from airfoil leading edge to vertical outlet boundary");
96 "Define obliqueness of the vertical mesh around the airfoil");
105 "Type of airfoil geometry, either NACA or Joukowski airfoil",
120 "Joukowski circle center coordinates");
123 "Joukowski airfoil length leading to trailing edge");
131 "Number of global refinements");
133 "NumberSubdivisionX0",
135 "Number of subdivisions along the airfoil in blocks with material ID 1 and 4");
137 "NumberSubdivisionX1",
139 "Number of subdivisions along the airfoil in blocks with material ID 2 and 5");
141 "NumberSubdivisionX2",
143 "Number of subdivisions in horizontal direction on the right of the trailing edge, i.e., blocks with material ID 3 and 6");
146 "Number of subdivisions normal to airfoil");
150 "Factor to obtain a finer mesh at the airfoil surface");
176 : refinements(
data.refinements)
177 , n_subdivision_x_0(
data.n_subdivision_x_0)
178 , n_subdivision_x_1(
data.n_subdivision_x_1)
179 , n_subdivision_x_2(
data.n_subdivision_x_2)
180 , n_subdivision_y(
data.n_subdivision_y)
181 , height(
data.height)
182 , length_b2(
data.length_b2)
183 , incline_factor(
data.incline_factor)
184 , bias_factor(
data.bias_factor)
194 data.airfoil_type ==
"Joukowski" ?
197 data.airfoil_sampling_factor) :
198 (
data.airfoil_type ==
"NACA" ?
201 data.airfoil_sampling_factor) :
204 std::vector<Point<2>>{
208 data.airfoil_length))
216 ,
gamma(std::atan(height /
230 , J(
tail_x + length_b2, 0)
235 data.airfoil_type ==
"NACA",
280 const unsigned int refinements;
283 const unsigned int n_subdivision_x_0;
286 const unsigned int n_subdivision_x_1;
289 const unsigned int n_subdivision_x_2;
293 const unsigned int n_subdivision_y;
300 const double length_b2;
304 const double incline_factor;
307 const double bias_factor;
330 const std::array<std::vector<Point<2>>, 2>
airfoil_1D;
382 const Point<2> A,
B,
C,
D,
E,
F, G, H, I, J,
K, L;
421 static std::array<std::vector<Point<2>>, 2>
424 const unsigned int factor)
426 std::array<std::vector<Point<2>>, 2>
airfoil_1D;
463 for (
unsigned int i = 0; i <=
nose_index; ++i)
479 for (
auto &i : vector)
513 static std::vector<Point<2>>
533 "Error creating lower circle: Circle for Joukowski-transform does"
534 " not enclose point zeta = -1! Choose different center "
555 static std::vector<Point<2>>
566 const std::complex<double> z =
zeta + 1. /
zeta;
589 static std::array<std::vector<Point<2>>, 2>
592 const unsigned int factor)
618 static std::vector<Point<2>>
624 ExcMessage(
"This NACA-serial number is not implemented!"));
645 static std::vector<Point<2>>
659 const double t =
static_cast<double>(
digit_23) / 100.0;
670 0.3516 * Utilities::fixed_power<2>(
x) +
671 0.2843 * Utilities::fixed_power<3>(
x) -
672 0.1036 * Utilities::fixed_power<4>(
683 const double m = 1.0 *
digit_0 / 100;
684 const double p = 1.0 *
digit_1 / 10;
689 m / Utilities::fixed_power<2>(p) *
690 (2 * p *
x - Utilities::fixed_power<2>(
x)) :
695 (
x <= p) ? 2 * m / Utilities::fixed_power<2>(p) * (p -
x) :
701 0.3516 * Utilities::fixed_power<2>(
x) +
702 0.2843 * Utilities::fixed_power<3>(
x) -
703 0.1036 * Utilities::fixed_power<4>(
706 const double theta = std::atan(
dy_c);
729 static std::array<std::vector<Point<2>>, 2>
733 std::array<std::vector<Point<2>>, 2> output;
747 static std::vector<Point<2>>
751 std::vector<Point<2>> output = input;
754 desired_len / input.front().distance(input.back());
756 for (
auto &
x : output)
772 static std::vector<Point<2>>
777 const unsigned int n_points =
789 const auto airfoil_length =
803 for (
unsigned int j = 0, i = 1;
j < n_points - 1; ++
j)
835 std::vector<Triangulation<2>>
trias(10);
864 for (
auto cell :
tria.active_cell_iterators())
872 {n_subdivision_y, n_subdivision_x_0},
876 {n_subdivision_y, n_subdivision_x_0},
880 {n_subdivision_x_1, n_subdivision_y},
884 {n_subdivision_x_1, n_subdivision_y},
888 {n_subdivision_x_2, n_subdivision_y},
892 {n_subdivision_x_2, n_subdivision_y},
919 for (
auto cell :
tria.active_cell_iterators())
922 if (cell->face(f)->at_boundary() ==
false)
925 const auto mid = cell->material_id();
929 cell->face(f)->set_boundary_id(0);
932 cell->face(f)->set_boundary_id(1);
935 cell->face(f)->set_boundary_id(2);
938 cell->face(f)->set_boundary_id(3);
941 cell->face(f)->set_boundary_id(4);
944 cell->face(f)->set_boundary_id(5);
981 for (
auto &cell :
tria)
991 auto &
node = cell.vertex(v);
1024 const double L = height /
std::sin(gamma);
1036 static_cast<int>(std::round((
node_(0) -
deltax) / dx));
1050 static_cast<int>(std::round(
node_(0) / dx));
1052 static_cast<int>(std::round(
node_(1) / dy));
1053 const double alpha =
1055 const double theta =
node_(0);
1077 "Points D,C,G and E,F,I are not defined symmetric to "
1078 "x-axis, which is required to interpolate block 2"
1079 " and 5 with same geometric computations."));
1080 const double l_y =
D(1) -
C(1);
1081 const double l_h =
D(1) -
l_y;
1084 const int iy =
static_cast<int>(
1087 const int ix =
static_cast<int>(
1095 incline_factor * length_b2 *
ix /
1103 (cell.material_id() ==
id_block_2) ? (0) : (1))]
1114 const int ix =
static_cast<int>(
1125 const Point<2> p1(J(0) - (1 - incline_factor) * length_b2 *
1155 return std::tanh(bias_factor *
alpha) / std::tanh(bias_factor);
1167 const AdditionalData & additional_data)
1198 const AdditionalData &)
1208 const AdditionalData &additional_data)
1221 const AdditionalData & additional_data)
1234 const AdditionalData & additional_data)
1238 (void)additional_data;
1250 template <
int dim,
int spacedim>
1262 cell->face(f)->set_boundary_id(f);
1268 template <
int spacedim>
1279 if (cell->center()(0) > 0)
1280 cell->set_material_id(1);
1287 template <
int dim,
int spacedim>
1292 const double epsilon)
1305 for (; face !=
endface; ++face)
1306 if (face->at_boundary())
1307 if (face->boundary_id() == 0)
1312 face->set_boundary_id(0);
1314 face->set_boundary_id(1);
1316 face->set_boundary_id(2);
1318 face->set_boundary_id(3);
1320 face->set_boundary_id(4);
1322 face->set_boundary_id(5);
1331 for (
const auto &cell :
tria.cell_iterators())
1334 for (
unsigned int d = 0;
d < dim; ++
d)
1335 if (cell->center()(d) > 0)
1337 cell->set_material_id(
id);
1363 cell->face(2)->set_all_boundary_ids(1);
1382 if (
tria.n_cells() == 6)
1387 cell->face(4)->set_all_boundary_ids(1);
1391 cell->face(2)->set_all_boundary_ids(1);
1395 cell->face(2)->set_all_boundary_ids(1);
1399 cell->face(0)->set_all_boundary_ids(1);
1403 cell->face(2)->set_all_boundary_ids(1);
1407 cell->face(0)->set_all_boundary_ids(1);
1409 else if (
tria.n_cells() == 12)
1418 cell->face(5)->set_all_boundary_ids(1);
1421 else if (
tria.n_cells() == 96)
1429 unsigned int count = 0;
1431 for (
const auto &cell :
tria.cell_iterators())
1432 if (cell->face(5)->at_boundary())
1434 cell->face(5)->set_all_boundary_ids(1);
1458 if (
tria.n_cells() != 3)
1465 for (; cell !=
tria.
end(); ++cell)
1466 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
1468 if (!cell->face(f)->at_boundary())
1471 double radius = cell->face(f)->center().norm() -
center.
norm();
1472 if (std::fabs(cell->face(f)->center()(0)) <
1475 cell->face(f)->set_boundary_id(2);
1478 if (cell->face(f)->line(
j)->at_boundary())
1479 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1480 cell->face(f)->line(
j)->vertex(1).norm()) >
1482 cell->face(f)->line(
j)->set_boundary_id(2);
1484 else if (std::fabs(cell->face(f)->center()(1)) <
1487 cell->face(f)->set_boundary_id(3);
1490 if (cell->face(f)->line(
j)->at_boundary())
1491 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1492 cell->face(f)->line(
j)->vertex(1).norm()) >
1494 cell->face(f)->line(
j)->set_boundary_id(3);
1496 else if (std::fabs(cell->face(f)->center()(2)) <
1499 cell->face(f)->set_boundary_id(4);
1502 if (cell->face(f)->line(
j)->at_boundary())
1503 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1504 cell->face(f)->line(
j)->vertex(1).norm()) >
1506 cell->face(f)->line(
j)->set_boundary_id(4);
1508 else if (radius <
middle)
1510 cell->face(f)->set_boundary_id(0);
1513 if (cell->face(f)->line(
j)->at_boundary())
1514 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1515 cell->face(f)->line(
j)->vertex(1).norm()) <
1517 cell->face(f)->line(
j)->set_boundary_id(0);
1519 else if (radius >
middle)
1521 cell->face(f)->set_boundary_id(1);
1524 if (cell->face(f)->line(
j)->at_boundary())
1525 if (std::fabs(cell->face(f)->line(
j)->vertex(0).norm() -
1526 cell->face(f)->line(
j)->vertex(1).norm()) <
1528 cell->face(f)->line(
j)->set_boundary_id(1);
1538 template <
int dim,
int spacedim>
1549 for (
unsigned int i = 0; i < dim; ++i)
1589 std::vector<CellData<dim>> cells(1);
1592 cells[0].material_id = 0;
1603 template <
int dim,
int spacedim>
1611 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
1614 for (
unsigned int i = 0; i < dim; ++i)
1634 for (
unsigned int d = 0;
d < dim; ++
d)
1635 for (
unsigned int c = 1; c <= dim; ++c)
1638 ExcMessage(
"Vertices of simplex must form a right handed system"));
1642 std::vector<Point<dim>> points =
vertices;
1646 for (
unsigned int i = 0; i <= dim; ++i)
1648 points.push_back(0.5 * (points[i] + points[(i + 1) % (dim + 1)]));
1654 for (
unsigned int i = 1; i < dim; ++i)
1655 points.push_back(0.5 * (points[i - 1] + points[i + 1]));
1657 for (
unsigned int i = 0; i <= dim; ++i)
1658 points.push_back(1. / 3. *
1659 (points[i] + points[(i + 1) % (dim + 1)] +
1660 points[(i + 2) % (dim + 1)]));
1662 points.push_back((1. / (dim + 1)) *
center);
1664 std::vector<CellData<dim>> cells(dim + 1);
1669 cells[0].vertices[0] = 0;
1670 cells[0].vertices[1] = 3;
1671 cells[0].vertices[2] = 5;
1672 cells[0].vertices[3] = 6;
1673 cells[0].material_id = 0;
1675 cells[1].vertices[0] = 3;
1676 cells[1].vertices[1] = 1;
1677 cells[1].vertices[2] = 6;
1678 cells[1].vertices[3] = 4;
1679 cells[1].material_id = 0;
1681 cells[2].vertices[0] = 5;
1682 cells[2].vertices[1] = 6;
1683 cells[2].vertices[2] = 2;
1684 cells[2].vertices[3] = 4;
1685 cells[2].material_id = 0;
1689 cells[0].vertices[0] = 0;
1690 cells[0].vertices[1] = 4;
1691 cells[0].vertices[2] = 8;
1692 cells[0].vertices[3] = 10;
1693 cells[0].vertices[4] = 7;
1694 cells[0].vertices[5] = 13;
1695 cells[0].vertices[6] = 12;
1696 cells[0].vertices[7] = 14;
1697 cells[0].material_id = 0;
1699 cells[1].vertices[0] = 4;
1700 cells[1].vertices[1] = 1;
1701 cells[1].vertices[2] = 10;
1702 cells[1].vertices[3] = 5;
1703 cells[1].vertices[4] = 13;
1704 cells[1].vertices[5] = 9;
1705 cells[1].vertices[6] = 14;
1706 cells[1].vertices[7] = 11;
1707 cells[1].material_id = 0;
1709 cells[2].vertices[0] = 8;
1710 cells[2].vertices[1] = 10;
1711 cells[2].vertices[2] = 2;
1712 cells[2].vertices[3] = 5;
1713 cells[2].vertices[4] = 12;
1714 cells[2].vertices[5] = 14;
1715 cells[2].vertices[6] = 6;
1716 cells[2].vertices[7] = 11;
1717 cells[2].material_id = 0;
1719 cells[3].vertices[0] = 7;
1720 cells[3].vertices[1] = 13;
1721 cells[3].vertices[2] = 12;
1722 cells[3].vertices[3] = 14;
1723 cells[3].vertices[4] = 3;
1724 cells[3].vertices[5] = 9;
1725 cells[3].vertices[6] = 6;
1726 cells[3].vertices[7] = 11;
1727 cells[3].material_id = 0;
1737 template <
int dim,
int spacedim>
1744 if (reference_cell == ReferenceCells::get_hypercube<dim>())
1757 for (
unsigned int d = 0;
d < dim; ++
d)
1766 std::vector<CellData<dim>> cells(1);
1778 const unsigned int n_cells,
1783 const unsigned int dim = 3;
1786 "More than 4 cells are needed to create a moebius grid."));
1788 ExcMessage(
"Outer and inner radius must be positive."));
1790 ExcMessage(
"Outer radius must be greater than inner radius."));
1793 std::vector<Point<dim>>
vertices(4 * n_cells);
1797 for (
unsigned int i = 0; i <
n_cells; ++i)
1798 for (
unsigned int j = 0;
j < 4; ++
j)
1812 unsigned int offset = 0;
1818 {0, 1, 5, 4, 2, 3, 7, 6}};
1819 std::vector<CellData<dim>> cells(n_cells);
1820 for (
unsigned int i = 0; i <
n_cells; ++i)
1822 for (
unsigned int j = 0;
j < 2; ++
j)
1834 cells[i].material_id = 0;
1862 ExcMessage(
"Outer radius R must be greater than the inner "
1866 const unsigned int dim = 2;
1867 const unsigned int spacedim = 3;
1868 std::vector<Point<spacedim>>
vertices(16);
1887 std::vector<CellData<dim>> cells(16);
1889 cells[0].vertices[0] = 0;
1890 cells[0].vertices[1] = 4;
1891 cells[0].vertices[2] = 3;
1892 cells[0].vertices[3] = 7;
1893 cells[0].material_id = 0;
1895 cells[1].vertices[0] = 1;
1896 cells[1].vertices[1] = 5;
1897 cells[1].vertices[2] = 0;
1898 cells[1].vertices[3] = 4;
1899 cells[1].material_id = 0;
1901 cells[2].vertices[0] = 2;
1902 cells[2].vertices[1] = 6;
1903 cells[2].vertices[2] = 1;
1904 cells[2].vertices[3] = 5;
1905 cells[2].material_id = 0;
1907 cells[3].vertices[0] = 3;
1908 cells[3].vertices[1] = 7;
1909 cells[3].vertices[2] = 2;
1910 cells[3].vertices[3] = 6;
1911 cells[3].material_id = 0;
1913 cells[4].vertices[0] = 4;
1914 cells[4].vertices[1] = 8;
1915 cells[4].vertices[2] = 7;
1916 cells[4].vertices[3] = 11;
1917 cells[4].material_id = 0;
1919 cells[5].vertices[0] = 5;
1920 cells[5].vertices[1] = 9;
1921 cells[5].vertices[2] = 4;
1922 cells[5].vertices[3] = 8;
1923 cells[5].material_id = 0;
1925 cells[6].vertices[0] = 6;
1926 cells[6].vertices[1] = 10;
1927 cells[6].vertices[2] = 5;
1928 cells[6].vertices[3] = 9;
1929 cells[6].material_id = 0;
1931 cells[7].vertices[0] = 7;
1932 cells[7].vertices[1] = 11;
1933 cells[7].vertices[2] = 6;
1934 cells[7].vertices[3] = 10;
1935 cells[7].material_id = 0;
1937 cells[8].vertices[0] = 8;
1938 cells[8].vertices[1] = 12;
1939 cells[8].vertices[2] = 11;
1940 cells[8].vertices[3] = 15;
1941 cells[8].material_id = 0;
1943 cells[9].vertices[0] = 9;
1944 cells[9].vertices[1] = 13;
1945 cells[9].vertices[2] = 8;
1946 cells[9].vertices[3] = 12;
1947 cells[9].material_id = 0;
1949 cells[10].vertices[0] = 10;
1950 cells[10].vertices[1] = 14;
1951 cells[10].vertices[2] = 9;
1952 cells[10].vertices[3] = 13;
1953 cells[10].material_id = 0;
1955 cells[11].vertices[0] = 11;
1956 cells[11].vertices[1] = 15;
1957 cells[11].vertices[2] = 10;
1958 cells[11].vertices[3] = 14;
1959 cells[11].material_id = 0;
1961 cells[12].vertices[0] = 12;
1962 cells[12].vertices[1] = 0;
1963 cells[12].vertices[2] = 15;
1964 cells[12].vertices[3] = 3;
1965 cells[12].material_id = 0;
1967 cells[13].vertices[0] = 13;
1968 cells[13].vertices[1] = 1;
1969 cells[13].vertices[2] = 12;
1970 cells[13].vertices[3] = 0;
1971 cells[13].material_id = 0;
1973 cells[14].vertices[0] = 14;
1974 cells[14].vertices[1] = 2;
1975 cells[14].vertices[2] = 13;
1976 cells[14].vertices[3] = 1;
1977 cells[14].material_id = 0;
1979 cells[15].vertices[0] = 15;
1980 cells[15].vertices[1] = 3;
1981 cells[15].vertices[2] = 14;
1982 cells[15].vertices[3] = 2;
1983 cells[15].material_id = 0;
1988 tria.set_all_manifold_ids(0);
2003 ExcMessage(
"Outer radius R must be greater than the inner "
2007 ExcMessage(
"Number of cells in toroidal direction has "
2008 "to be at least 3."));
2014 double const a = 1. / (1 +
std::sqrt(2.0));
2038 for (
unsigned int v = 0; v < 8; ++v)
2051 for (
unsigned int j = 0;
j < 2; ++
j)
2053 const unsigned int offset =
2057 cells[5 * c].vertices[0 +
j * 4] = offset + 0;
2058 cells[5 * c].vertices[1 +
j * 4] = offset + 1;
2059 cells[5 * c].vertices[2 +
j * 4] = offset + 2;
2060 cells[5 * c].vertices[3 +
j * 4] = offset + 3;
2062 cells[5 * c + 1].vertices[0 +
j * 4] = offset + 2;
2063 cells[5 * c + 1].vertices[1 +
j * 4] = offset + 3;
2064 cells[5 * c + 1].vertices[2 +
j * 4] = offset + 4;
2065 cells[5 * c + 1].vertices[3 +
j * 4] = offset + 5;
2067 cells[5 * c + 2].vertices[0 +
j * 4] = offset + 4;
2068 cells[5 * c + 2].vertices[1 +
j * 4] = offset + 5;
2069 cells[5 * c + 2].vertices[2 +
j * 4] = offset + 6;
2070 cells[5 * c + 2].vertices[3 +
j * 4] = offset + 7;
2072 cells[5 * c + 3].vertices[0 +
j * 4] = offset + 0;
2073 cells[5 * c + 3].vertices[1 +
j * 4] = offset + 2;
2074 cells[5 * c + 3].vertices[2 +
j * 4] = offset + 6;
2075 cells[5 * c + 3].vertices[3 +
j * 4] = offset + 4;
2077 cells[5 * c + 4].vertices[0 +
j * 4] = offset + 3;
2078 cells[5 * c + 4].vertices[1 +
j * 4] = offset + 1;
2079 cells[5 * c + 4].vertices[2 +
j * 4] = offset + 5;
2080 cells[5 * c + 4].vertices[3 +
j * 4] = offset + 7;
2083 cells[5 * c].material_id = 0;
2085 cells[5 * c + 1].material_id = 1;
2086 cells[5 * c + 2].material_id = 0;
2087 cells[5 * c + 3].material_id = 0;
2088 cells[5 * c + 4].material_id = 0;
2093 tria.reset_all_manifolds();
2094 tria.set_all_manifold_ids(0);
2096 for (
auto &cell :
tria.cell_iterators())
2099 for (
const unsigned int f :
GeometryInfo<3>::face_indices())
2103 if (cell->face(f)->at_boundary() && f != 4 && f != 5)
2105 cell->face(f)->set_all_manifold_ids(1);
2110 if (cell->material_id() == 1)
2112 cell->set_all_manifold_ids(2);
2114 cell->set_material_id(0);
2119 tria.set_manifold(2,
2129 template <
int dim,
int spacedim>
2142 tria.begin_active();
2150 "The volume of the cell is not greater than zero. "
2151 "This could be due to the wrong ordering of the vertices."));
2182 std::array<Tensor<1, 2>, 2>
edges;
2198 unsigned int n_subdivisions[dim];
2199 for (
unsigned int i = 0; i < dim; ++i)
2200 n_subdivisions[i] = 1;
2209 const unsigned int n_subdivisions,
2216 for (
unsigned int i = 0; i < dim; ++i)
2227 const unsigned int (&n_subdivisions)[dim],
2229 const unsigned int *n_subdivisions,
2236 std::array<Tensor<1, dim>, dim>
edges;
2237 for (
unsigned int i = 0; i < dim; ++i)
2250 template <
int dim,
int spacedim>
2265 ExcMessage(
"One subdivision must be provided for each dimension."));
2267 for (
unsigned int i = 0; i < dim; ++i)
2274 "Edges in subdivided_parallelepiped() must not be degenerate."));
2308 "Edges in subdivided_parallelepiped() must point in"
2309 " different directions."));
2334 ExcInvalidInputOrientation(
2335 "The triangulation you are trying to create will consist of cells"
2336 " with negative measures. This is usually the result of input data"
2337 " that does not define a right-handed coordinate system. The usual"
2338 " fix for this is to ensure that in 1d the given point is to the"
2339 " right of the origin (or the given edge tensor is positive), in 2d"
2340 " that the two edges (and their cross product) obey the right-hand"
2341 " rule (which may usually be done by switching the order of the"
2342 " points or edge tensors), or in 3d that the edges form a"
2343 " right-handed coordinate system (which may also be accomplished by"
2344 " switching the order of the first two points or edge tensors)."));
2347 for (
unsigned int i = 0; i < dim; ++i)
2348 for (
unsigned int j = i + 1;
j < dim; ++
j)
2351 "Degenerate edges of subdivided_parallelepiped encountered."));
2354 std::vector<Point<spacedim>> points;
2374 points.push_back(
origin +
2386 for (
unsigned int i = 0; i < dim; ++i)
2388 std::vector<CellData<dim>> cells(n_cells);
2396 cells[
x].vertices[0] =
x;
2397 cells[
x].vertices[1] =
x + 1;
2400 cells[
x].material_id = 0;
2410 for (
unsigned int y = 0;
y <
n_dy; ++
y)
2411 for (
unsigned int x = 0;
x <
n_dx; ++
x)
2413 const unsigned int c =
y *
n_dx +
x;
2414 cells[c].vertices[0] =
y * (
n_dx + 1) +
x;
2415 cells[c].vertices[1] =
y * (
n_dx + 1) +
x + 1;
2416 cells[c].vertices[2] = (
y + 1) * (
n_dx + 1) +
x;
2417 cells[c].vertices[3] = (
y + 1) * (
n_dx + 1) +
x + 1;
2420 cells[c].material_id = 0;
2432 for (
unsigned int z = 0; z <
n_dz; ++z)
2433 for (
unsigned int y = 0;
y <
n_dy; ++
y)
2434 for (
unsigned int x = 0;
x <
n_dx; ++
x)
2438 cells[c].vertices[0] =
2440 cells[c].vertices[1] =
2442 cells[c].vertices[2] =
2444 cells[c].vertices[3] = z * (
n_dy + 1) * (
n_dx + 1) +
2445 (
y + 1) * (
n_dx + 1) +
x + 1;
2446 cells[c].vertices[4] =
2448 cells[c].vertices[5] = (z + 1) * (
n_dy + 1) * (
n_dx + 1) +
2450 cells[c].vertices[6] = (z + 1) * (
n_dy + 1) * (
n_dx + 1) +
2451 (
y + 1) * (
n_dx + 1) +
x;
2452 cells[c].vertices[7] = (z + 1) * (
n_dy + 1) * (
n_dx + 1) +
2453 (
y + 1) * (
n_dx + 1) +
x + 1;
2456 cells[c].material_id = 0;
2475 tria.begin_active(),
2477 for (; cell !=
endc; ++cell)
2479 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
2481 if (cell->face(face)->at_boundary())
2482 cell->face(face)->set_boundary_id(face);
2489 template <
int dim,
int spacedim>
2499 ExcMessage(
"Invalid left-to-right bounds of hypercube"));
2502 for (
unsigned int i = 0; i < dim; ++i)
2514 template <
int dim,
int spacedim>
2528 for (
unsigned int i = 0; i < dim; ++i)
2535 std::array<Point<spacedim>, dim> delta;
2536 for (
unsigned int i = 0; i < dim; ++i)
2544 "The first dim entries of coordinates of p1 and p2 need to be different."));
2548 std::vector<Point<spacedim>> points;
2553 points.push_back(
p1 +
x * delta[0]);
2559 points.push_back(
p1 +
x * delta[0] +
y * delta[1]);
2563 for (
unsigned int z = 0; z <=
repetitions[2]; ++z)
2566 points.push_back(
p1 +
x * delta[0] +
y * delta[1] +
2575 std::vector<CellData<dim>> cells;
2583 cells[
x].vertices[0] =
x;
2584 cells[
x].vertices[1] =
x + 1;
2585 cells[
x].material_id = 0;
2600 cells[c].vertices[3] = (
y + 1) * (
repetitions[0] + 1) +
x + 1;
2601 cells[c].material_id = 0;
2609 const unsigned int n_xy =
2619 cells[c].vertices[0] = z *
n_xy +
y *
n_x +
x;
2620 cells[c].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
2621 cells[c].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
2622 cells[c].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
2623 cells[c].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
2624 cells[c].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
2625 cells[c].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
2626 cells[c].vertices[7] =
2628 cells[c].material_id = 0;
2650 for (
unsigned int i = 0; i < dim; ++i)
2651 epsilon =
std::min(epsilon, 0.01 * delta[i][i]);
2654 "The distance between corner points must be positive."))
2684 for (
unsigned int i = 0; i < dim; ++i)
2688 std::swap(
p1(i),
p2(i));
2696 Assert(std::fabs(
x - (
p2(i) -
p1(i))) <= 1e-12 * std::fabs(
x),
2698 "The sequence of step sizes in coordinate direction " +
2700 " must be equal to the distance of the two given "
2701 "points in this coordinate direction."));
2708 std::vector<Point<dim>> points;
2714 for (
unsigned int i = 0;; ++i)
2734 for (
unsigned int j = 0;; ++
j)
2737 for (
unsigned int i = 0;; ++i)
2756 for (
unsigned int k = 0;; ++
k)
2759 for (
unsigned int j = 0;; ++
j)
2762 for (
unsigned int i = 0;; ++i)
2792 std::vector<CellData<dim>> cells;
2800 cells[
x].vertices[0] =
x;
2801 cells[
x].vertices[1] =
x + 1;
2802 cells[
x].material_id = 0;
2816 cells[c].vertices[2] =
2818 cells[c].vertices[3] =
2820 cells[c].material_id = 0;
2828 const unsigned int n_xy =
2837 const unsigned int c =
2840 cells[c].vertices[0] = z *
n_xy +
y *
n_x +
x;
2841 cells[c].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
2842 cells[c].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
2843 cells[c].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
2844 cells[c].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
2845 cells[c].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
2846 cells[c].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
2847 cells[c].vertices[7] =
2849 cells[c].material_id = 0;
2872 for (
unsigned int i = 1; i < dim; ++i)
2889 const std::vector<std::vector<double>> &
spacing,
2898 Assert(
spacing[0].size() == n_cells, ExcInvalidRepetitionsDimension(1));
2900 double delta = std::numeric_limits<double>::max();
2901 for (
unsigned int i = 0; i <
n_cells; ++i)
2908 std::vector<Point<1>> points;
2912 points.emplace_back(
ax);
2918 for (
unsigned int i = 0; i <
n_cells; ++i)
2923 unsigned int id = 0;
2927 cells[id].vertices[0] =
x;
2928 cells[id].vertices[1] =
x + 1;
2936 tria.create_triangulation(points, cells, t);
2947 const std::vector<std::vector<double>> &
spacing,
2955 double delta = std::numeric_limits<double>::max();
2956 for (
unsigned int i = 0; i < 2; ++i)
2965 ExcInvalidRepetitionsDimension(i));
2969 std::vector<Point<2>> points;
2976 points.emplace_back(
ax,
ay);
2986 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
2992 unsigned int id = 0;
3000 cells[id].vertices[3] = (
y + 1) * (
repetitions[0] + 1) +
x + 1;
3009 tria.create_triangulation(points, cells, t);
3014 double eps = 0.01 * delta;
3016 for (; cell !=
endc; ++cell)
3019 for (
const unsigned int f :
GeometryInfo<2>::face_indices())
3023 for (
unsigned int i = 0; i < 2; ++i)
3026 cell->face(f)->set_boundary_id(i * 2);
3028 cell->face(f)->set_boundary_id(i * 2 + 1);
3039 const std::vector<std::vector<double>> &
spacing,
3044 const unsigned int dim = 3;
3049 double delta = std::numeric_limits<double>::max();
3050 for (
unsigned int i = 0; i < dim; ++i)
3059 ExcInvalidRepetitionsDimension(i));
3063 std::vector<Point<dim>> points;
3065 for (
unsigned int z = 0; z <=
repetitions[2]; ++z)
3073 points.emplace_back(
ax,
ay,
az);
3086 for (
unsigned int i = 0; i <
material_id.size(0); ++i)
3093 unsigned int id = 0;
3101 cells[id].vertices[0] = z *
n_xy +
y *
n_x +
x;
3102 cells[id].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
3103 cells[id].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
3104 cells[id].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
3105 cells[id].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
3106 cells[id].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
3107 cells[id].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
3108 cells[id].vertices[7] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x + 1;
3117 tria.create_triangulation(points, cells, t);
3122 double eps = 0.01 * delta;
3125 for (; cell !=
endc; ++cell)
3132 for (
unsigned int i = 0; i < dim; ++i)
3135 cell->face(f)->set_boundary_id(i * 2);
3137 cell->face(f)->set_boundary_id(i * 2 + 1);
3144 template <
int dim,
int spacedim>
3147 const std::vector<unsigned int> &
holes)
3156 for (
unsigned int d = 0;
d < dim; ++
d)
3163 std::array<Point<spacedim>, dim> delta;
3165 for (
unsigned int i = 0; i < dim; ++i)
3168 ExcMessage(
"At least one hole needed in each direction"));
3170 delta[i][i] = (
p2[i] -
p1[i]);
3175 std::vector<Point<spacedim>> points;
3180 points.push_back(
p1 +
x * delta[0]);
3186 points.push_back(
p1 +
x * delta[0] +
y * delta[1]);
3190 for (
unsigned int z = 0; z <=
repetitions[2]; ++z)
3193 points.push_back(
p1 +
x * delta[0] +
y * delta[1] +
3203 std::vector<CellData<dim>> cells;
3213 if ((
x % 2 == 1) && (
y % 2 == 1))
3219 cells[c].vertices[3] = (
y + 1) * (
repetitions[0] + 1) +
x + 1;
3220 cells[c].material_id = 0;
3229 const unsigned int n_xy =
3240 cells[c].vertices[0] = z *
n_xy +
y *
n_x +
x;
3241 cells[c].vertices[1] = z *
n_xy +
y *
n_x +
x + 1;
3242 cells[c].vertices[2] = z *
n_xy + (
y + 1) *
n_x +
x;
3243 cells[c].vertices[3] = z *
n_xy + (
y + 1) *
n_x +
x + 1;
3244 cells[c].vertices[4] = (z + 1) *
n_xy +
y *
n_x +
x;
3245 cells[c].vertices[5] = (z + 1) *
n_xy +
y *
n_x +
x + 1;
3246 cells[c].vertices[6] = (z + 1) *
n_xy + (
y + 1) *
n_x +
x;
3247 cells[c].vertices[7] =
3249 cells[c].material_id = 0;
3277 const unsigned int ,
3289 const unsigned int ,
3303 const double radius)
3305 return (
std::abs(p[0] - c[0]) < radius) &&
3314 template <
int dim,
int spacedim>
3318 double length = std::numeric_limits<double>::max();
3319 for (
const auto &cell :
triangulation.active_cell_iterators())
3321 length =
std::min(length, cell->line(n)->diameter());
3341 const unsigned int ,
3355 double length = std::numeric_limits<double>::max();
3356 for (
const auto &cell :
tria.active_cell_iterators())
3358 length =
std::min(length, cell->line(n)->diameter());
3373 for (
const auto &cell :
cylinder_tria.active_cell_iterators())
3385 const double h) ->
void {
3389 static_cast<unsigned int>(std::round(
padding / h));
3393 for (
unsigned int i = 0; i <
num; ++i)
3397 std::vector<std::vector<double>>
step_sizes(2);
3422 for (
const auto &cell :
bulk_tria.active_cell_iterators())
3430 const double tolerance =
3442 for (
const auto &cell :
tria.active_cell_iterators())
3450 const auto &face = cell->face(
face_n);
3451 if (face->at_boundary() &&
3452 internal::point_in_2d_box(face->center(),
3468 static constexpr double tol =
3469 std::numeric_limits<double>::epsilon() * 10000;
3471 for (
const auto &cell :
tria.active_cell_iterators())
3474 const auto face = cell->face(
face_n);
3475 if (face->at_boundary())
3480 face->set_boundary_id(0);
3483 face->set_boundary_id(1);
3486 face->set_boundary_id(2);
3489 face->set_boundary_id(3);
3495 face->set_boundary_id(4);
3571 ExcMessage(
"The width of the shell region must be less than 0.05 "
3572 "(and preferably close to 0.03)"));
3608 for (
const auto &cell :
bulk_tria.active_cell_iterators())
3610 if ((cell->center() -
Point<2>(0.2, 0.2)).norm() < 0.15)
3622 2.0 * (cell->vertex(3) -
Point<2>());
3641 for (
const auto &cell :
cylinder_tria.active_cell_iterators())
3651 for (
const auto &cell :
cylinder_tria.active_cell_iterators())
3655 if (!cell->face(
face_n)->at_boundary())
3661 ExcMessage(
"If the shell region has positive width then "
3662 "there must be at least one shell."));
3707 const double shift =
3709 for (
const auto &cell :
tria.active_cell_iterators())
3711 if (cell->vertex(v).distance(
Point<2>(0.1, 0.205)) < 1
e-10)
3712 cell->vertex(v) =
Point<2>(0.2 -
shift, 0.205);
3713 else if (cell->vertex(v).distance(
Point<2>(0.3, 0.205)) < 1e-10)
3714 cell->vertex(v) =
Point<2>(0.2 + shift, 0.205);
3715 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.1025)) < 1e-10)
3716 cell->vertex(v) =
Point<2>(0.2, 0.2 - shift);
3717 else if (cell->vertex(v).distance(
Point<2>(0.2, 0.3075)) < 1e-10)
3718 cell->vertex(v) =
Point<2>(0.2, 0.2 + shift);
3723 for (
const auto &cell :
tria.active_cell_iterators())
3729 for (
const auto &cell :
tria.active_cell_iterators())
3737 for (
const auto &face :
tria.active_face_iterators())
3766 for (
const auto &face :
tria.active_face_iterators())
3767 if (face->at_boundary())
3772 face->set_boundary_id(0);
3775 face->set_boundary_id(1);
3778 face->set_boundary_id(2);
3785 face->set_boundary_id(3);
3825 for (
const auto &face :
tria.active_face_iterators())
3827 face->set_boundary_id(3);
3832 template <
int dim,
int spacedim>
3835 const std::vector<unsigned int> &sizes,
3845 for (
unsigned int d = 0;
d < dim; ++
d)
3848 std::vector<Point<spacedim>> points;
3850 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
3853 std::vector<CellData<dim>> cells(n_cells);
3858 for (
unsigned int d = 0;
d < dim; ++
d)
3859 p(d) = 0.5 * dimensions[
d] *
3862 points.push_back(p);
3863 cells[0].vertices[i] = i;
3865 cells[0].material_id = 0;
3870 for (
const unsigned int face :
GeometryInfo<dim>::face_indices())
3876 for (
unsigned int j = 0;
j < sizes[face]; ++
j, ++
cell_index)
3883 const unsigned int cellv =
3885 const unsigned int ocellv =
3897 points.push_back(p);
4058 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4060 std::vector<Point<2>>
vertices(16);
4068 for (
const double y : coords)
4074 std::vector<CellData<2>> cells(9);
4076 for (
unsigned int i0 = 0;
i0 < 3; ++
i0)
4077 for (
unsigned int i1 = 0;
i1 < 3; ++
i1)
4079 cells[
k].vertices[0] =
i1 + 4 *
i0;
4080 cells[
k].vertices[1] =
i1 + 4 *
i0 + 1;
4081 cells[
k].vertices[2] =
i1 + 4 *
i0 + 4;
4082 cells[
k].vertices[3] =
i1 + 4 *
i0 + 5;
4102 const double rl2 = (right + left) / 2;
4118 for (
unsigned int i = 0; i < 4; ++i)
4120 for (
unsigned int j = 0;
j < 4; ++
j)
4122 cells[i].material_id = 0;
4132 cell->face(1)->set_boundary_id(1);
4134 cell->face(0)->set_boundary_id(2);
4166 cells[0].material_id = 0;
4171 cell->face(0)->set_boundary_id(1);
4172 cell->face(1)->set_boundary_id(2);
4174 for (
unsigned int i = 2; i < 4; ++i)
4175 cell->face(i)->set_boundary_id(0);
4192 Point<2>((a + b) / 2, (a + b) / 2),
4196 const int cell_vertices[3][4] = {{0, 1, 3, 4}, {1, 2, 4, 5}, {3, 4, 6, 7}};
4200 for (
unsigned int i = 0; i < 3; ++i)
4202 for (
unsigned int j = 0;
j < 4; ++
j)
4204 cells[i].material_id = 0;
4216 cell->face(0)->set_boundary_id(0);
4217 cell->face(2)->set_boundary_id(1);
4220 cell->face(1)->set_boundary_id(2);
4221 cell->face(2)->set_boundary_id(1);
4222 cell->face(3)->set_boundary_id(3);
4225 cell->face(0)->set_boundary_id(0);
4226 cell->face(1)->set_boundary_id(4);
4227 cell->face(3)->set_boundary_id(5);
4233 template <
int dim,
int spacedim>
4245 for (
unsigned int d = 0;
d < dim; ++
d)
4248 ExcMessage(
"Attempting to cut away too many cells."));
4258 std::array<double, dim> h;
4260 for (
unsigned int d = 0;
d < dim; ++
d)
4280 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>
4282 for (
const auto &cell :
rectangle.active_cell_iterators())
4305 const double radius,
4311 const double a = 1. / (1 +
std::sqrt(2.0));
4323 {0, 1, 2, 3}, {0, 2, 6, 4}, {2, 3, 4, 5}, {1, 7, 3, 5}, {6, 4, 7, 5}};
4327 for (
unsigned int i = 0; i < 5; ++i)
4329 for (
unsigned int j = 0;
j < 4; ++
j)
4331 cells[i].material_id = 0;
4339 tria.set_all_manifold_ids_on_boundary(0);
4353 const unsigned int n_cells,
4370 const unsigned int N =
4371 (
n_cells == 0 ?
static_cast<unsigned int>(
4382 std::vector<Point<2>>
vertices(2 * N);
4383 for (
unsigned int i = 0; i < N; ++i)
4396 for (
unsigned int i = 0; i < N; ++i)
4398 cells[i].vertices[0] = i;
4399 cells[i].vertices[1] = (i + 1) % N;
4400 cells[i].vertices[2] = N + i;
4401 cells[i].vertices[3] = N + ((i + 1) % N);
4403 cells[i].material_id = 0;
4411 tria.set_all_manifold_ids(0);
4424 const unsigned int n_cells)
4433 "The inner radius is greater than or equal to the outer radius plus eccentricity."));
4439 for (
const auto &face :
tria.active_face_iterators())
4453 tria.set_all_manifold_ids(2);
4473 const double radius,
4485 switch (f->boundary_id())
4488 f->set_boundary_id(1);
4491 f->set_boundary_id(2);
4494 f->set_boundary_id(0);
4531 const double radius)
4533 const unsigned int dim = 2;
4547 const int cell_vertices[3][4] = {{0, 2, 3, 4}, {1, 6, 2, 4}, {5, 3, 6, 4}};
4551 for (
unsigned int i = 0; i < 3; ++i)
4553 for (
unsigned int j = 0;
j < 4; ++
j)
4555 cells[i].material_id = 0;
4566 tria.set_all_manifold_ids_on_boundary(0);
4570 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
4572 if (cell->face(i)->boundary_id() ==
4578 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius ||
4579 cell->face(i)->center()(1) < p(1) + 1.e-5 * radius)
4581 cell->face(i)->set_boundary_id(1);
4595 const double radius)
4600 const double a = 1. / (1 +
std::sqrt(2.0));
4618 for (
unsigned int i = 0; i < 4; ++i)
4620 for (
unsigned int j = 0;
j < 4; ++
j)
4622 cells[i].material_id = 0;
4633 tria.set_all_manifold_ids_on_boundary(0);
4637 for (
const unsigned int i :
GeometryInfo<2>::face_indices())
4639 if (cell->face(i)->boundary_id() ==
4644 if (cell->face(i)->center()(0) < p(0) + 1.e-5 * radius)
4646 cell->face(i)->set_boundary_id(1);
4664 const unsigned int n_cells,
4680 const unsigned int N =
4681 (
n_cells == 0 ?
static_cast<unsigned int>(
4692 std::vector<Point<2>>
vertices(2 * (N + 1));
4693 for (
unsigned int i = 0; i <= N; ++i)
4714 for (
unsigned int i = 0; i < N; ++i)
4716 cells[i].vertices[0] = i;
4717 cells[i].vertices[1] = (i + 1) % (N + 1);
4718 cells[i].vertices[2] = N + 1 + i;
4719 cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4721 cells[i].material_id = 0;
4729 for (; cell !=
tria.
end(); ++cell)
4731 cell->face(2)->set_boundary_id(1);
4733 tria.
begin()->face(0)->set_boundary_id(3);
4735 tria.last()->face(1)->set_boundary_id(2);
4737 tria.set_all_manifold_ids(0);
4748 const unsigned int n_cells,
4764 const unsigned int N =
4765 (
n_cells == 0 ?
static_cast<unsigned int>(
4776 std::vector<Point<2>>
vertices(2 * (N + 1));
4777 for (
unsigned int i = 0; i <= N; ++i)
4796 for (
unsigned int i = 0; i < N; ++i)
4798 cells[i].vertices[0] = i;
4799 cells[i].vertices[1] = (i + 1) % (N + 1);
4800 cells[i].vertices[2] = N + 1 + i;
4801 cells[i].vertices[3] = N + 1 + ((i + 1) % (N + 1));
4803 cells[i].material_id = 0;
4811 for (; cell !=
tria.
end(); ++cell)
4813 cell->face(2)->set_boundary_id(1);
4815 tria.
begin()->face(0)->set_boundary_id(3);
4817 tria.last()->face(1)->set_boundary_id(2);
4820 tria.set_all_manifold_ids(0);
4834 const double rl2 = (right + left) / 2;
4835 const double len = (right - left) / 2.;
4848 const int cell_vertices[4][8] = {{0, 1, 3, 2, 10, 11, 13, 12},
4849 {9, 4, 2, 5, 19, 14, 12, 15},
4850 {3, 2, 7, 6, 13, 12, 17, 16},
4851 {2, 5, 6, 8, 12, 15, 16, 18}};
4853 for (
unsigned int i = 0; i < 4; ++i)
4855 for (
unsigned int j = 0;
j < 8; ++
j)
4857 cells[i].material_id = 0;
4867 cell->face(1)->set_boundary_id(1);
4869 cell->face(0)->set_boundary_id(2);
4885 ExcMessage(
"Invalid left-to-right bounds of enclosed hypercube"));
4887 std::vector<Point<3>>
vertices(64);
4895 for (
const double z : coords)
4901 24, 26, 5, 4, 6, 1, 0,
4902 2, 9, 8, 10, 37, 36, 38,
4903 33, 32, 34, 41, 40, 42};
4905 std::vector<CellData<3>> cells(27);
4907 for (
unsigned int z = 0; z < 3; ++z)
4908 for (
unsigned int y = 0;
y < 3; ++
y)
4909 for (
unsigned int x = 0;
x < 3; ++
x)
4911 cells[
k].vertices[0] =
x + 4 *
y + 16 * z;
4912 cells[
k].vertices[1] =
x + 4 *
y + 16 * z + 1;
4913 cells[
k].vertices[2] =
x + 4 *
y + 16 * z + 4;
4914 cells[
k].vertices[3] =
x + 4 *
y + 16 * z + 5;
4915 cells[
k].vertices[4] =
x + 4 *
y + 16 * z + 16;
4916 cells[
k].vertices[5] =
x + 4 *
y + 16 * z + 17;
4917 cells[
k].vertices[6] =
x + 4 *
y + 16 * z + 20;
4918 cells[
k].vertices[7] =
x + 4 *
y + 16 * z + 21;
4938 ExcMessage(
"The output triangulation object needs to be empty."));
4943 const auto n_slices = 1 +
static_cast<unsigned int>(std::ceil(
4960 return Point<3>(p[0], factor * p[1], factor * p[2]);
4966 for (
const auto &face :
triangulation.active_face_iterators())
4967 if (face->at_boundary())
4970 face->set_boundary_id(1);
4973 face->set_boundary_id(2);
4975 face->set_all_manifold_ids(0);
4998 Point<3>((a + b) / 2, a, (a + b) / 2),
5005 Point<3>((a + b) / 2, (a + b) / 2, a),
5007 Point<3>(a, (a + b) / 2, (a + b) / 2),
5008 Point<3>((a + b) / 2, (a + b) / 2, (a + b) / 2),
5009 Point<3>(b, (a + b) / 2, (a + b) / 2),
5011 Point<3>((a + b) / 2, (a + b) / 2, b),
5019 Point<3>((a + b) / 2, b, (a + b) / 2),
5023 const int cell_vertices[7][8] = {{0, 1, 9, 10, 3, 4, 12, 13},
5024 {1, 2, 10, 11, 4, 5, 13, 14},
5025 {3, 4, 12, 13, 6, 7, 15, 16},
5026 {4, 5, 13, 14, 7, 8, 16, 17},
5027 {9, 10, 18, 19, 12, 13, 21, 22},
5028 {10, 11, 19, 20, 13, 14, 22, 23},
5029 {12, 13, 21, 22, 15, 16, 24, 25}};
5033 for (
unsigned int i = 0; i < 7; ++i)
5035 for (
unsigned int j = 0;
j < 8; ++
j)
5037 cells[i].material_id = 0;
5058 const double radius,
5065 const unsigned int n_vertices = 16;
5091 const unsigned int n_cells = 7;
5093 {0, 1, 4, 5, 3, 2, 7, 6},
5094 {8, 9, 12, 13, 0, 1, 4, 5},
5095 {9, 13, 1, 5, 10, 14, 2, 6},
5096 {11, 10, 3, 2, 15, 14, 7, 6},
5097 {8, 0, 12, 4, 11, 3, 15, 7},
5098 {8, 9, 0, 1, 11, 10, 3, 2},
5099 {12, 4, 13, 5, 15, 7, 14, 6}};
5101 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5103 for (
unsigned int i = 0; i <
n_cells; ++i)
5107 cells[i].material_id = 0;
5115 tria.set_all_manifold_ids_on_boundary(0);
5128 ExcMessage(
"The number of rotation by pi/2 of the right square "
5129 "must be in the half-open range [0,4)."))
5133 const unsigned int n_cells = 5;
5195 for (
const unsigned int vertex_index :
5210 const bool face_orientation,
5211 const bool face_flip,
5212 const bool face_rotation,
5215 constexpr unsigned int dim = 3;
5217 const unsigned int n_cells = 2;
5218 std::vector<CellData<dim>> cells(n_cells);
5235 {0, 1, 2, 3, 4, 5, 6, 7},
5236 {1, 8, 3, 9, 5, 10, 7, 11}};
5239 const unsigned int this_case = 4 *
static_cast<int>(face_orientation) +
5240 2 *
static_cast<int>(face_flip) +
5241 static_cast<int>(face_rotation);
5465 for (
const unsigned int vertex_index :
5479 template <
int spacedim>
5483 const double radius)
5487 const std::set<types::boundary_id> boundary_ids = {0};
5489 tria.set_all_manifold_ids(0);
5500 const double radius,
5516 vertices.emplace_back(-d, height, -d);
5517 vertices.emplace_back(d, height, -d);
5518 vertices.emplace_back(-a, height, -a);
5519 vertices.emplace_back(a, height, -a);
5520 vertices.emplace_back(-a, height, a);
5521 vertices.emplace_back(a, height, a);
5522 vertices.emplace_back(-d, height, d);
5523 vertices.emplace_back(d, height, d);
5529 const double h = vertex(1);
5530 vertex(1) = -vertex(0);
5543 for (
unsigned int i = 0; i < 5; ++i)
5546 for (
unsigned int j = 0;
j < 8; ++
j)
5554 std::vector<CellData<3>> cells(n_cells,
CellData<3>());
5556 for (
unsigned int i = 0; i <
n_cells; ++i)
5558 for (
unsigned int j = 0;
j < 8; ++
j)
5560 cells[i].material_id = 0;
5578 tria.set_all_manifold_ids_on_boundary(0);
5584 for (
const auto &cell :
tria.cell_iterators())
5586 if (cell->at_boundary(i))
5588 if (cell->face(i)->center()(0) >
half_length - tolerance)
5590 cell->face(i)->set_boundary_id(2);
5595 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5596 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5597 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5598 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5600 cell->face(i)->line(e)->set_boundary_id(2);
5601 cell->face(i)->line(e)->set_manifold_id(
5605 else if (cell->face(i)->center()(0) < -
half_length + tolerance)
5607 cell->face(i)->set_boundary_id(1);
5612 if ((std::fabs(cell->face(i)->line(e)->vertex(0)[1]) == a) ||
5613 (std::fabs(cell->face(i)->line(e)->vertex(0)[2]) == a) ||
5614 (std::fabs(cell->face(i)->line(e)->vertex(1)[1]) == a) ||
5615 (std::fabs(cell->face(i)->line(e)->vertex(1)[2]) == a))
5617 cell->face(i)->line(e)->set_boundary_id(1);
5618 cell->face(i)->line(e)->set_manifold_id(
5630 const double radius,
5640 const double radius)
5642 const unsigned int dim = 3;
5650 const double a = 0.528;
5651 const double b = 0.4533;
5652 const double c = 0.3752;
5669 const int cell_vertices[4][8] = {{0, 2, 3, 4, 7, 9, 10, 11},
5670 {1, 6, 2, 4, 8, 13, 9, 11},
5671 {5, 3, 6, 4, 12, 10, 13, 11},
5672 {7, 9, 10, 11, 14, 8, 12, 13}};
5676 for (
unsigned int i = 0; i < 4; ++i)
5678 for (
unsigned int j = 0;
j < 8; ++
j)
5680 cells[i].material_id = 0;
5691 tria.set_all_manifold_ids_on_boundary(0);
5694 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
5696 if (cell->face(i)->boundary_id() ==
5701 if (cell->face(i)->center()(0) <
center(0) + 1.e-5 * radius ||
5702 cell->face(i)->center()(1) <
center(1) + 1.e-5 * radius ||
5703 cell->face(i)->center()(2) <
center(2) + 1.e-5 * radius)
5705 cell->face(i)->set_boundary_id(1);
5712 const Point<3> line_vertices[2] = {
5713 cell->face(i)->line(
j)->vertex(0),
5714 cell->face(i)->line(
j)->vertex(1)};
5715 if ((std::fabs(line_vertices[0].distance(
center) - radius) >
5717 (std::fabs(line_vertices[1].distance(
center) - radius) >
5720 cell->face(i)->line(
j)->set_boundary_id(1);
5721 cell->face(i)->line(
j)->set_manifold_id(
5739 const double radius)
5745 const double b = a / 2.0;
5746 const double c =
d / 2.0;
5772 {0, 2, 8, 10, 6, 4, 14, 12},
5773 {2, 3, 10, 11, 4, 5, 12, 13},
5774 {1, 7, 9, 15, 3, 5, 11, 13},
5775 {6, 4, 14, 12, 7, 5, 15, 13},
5776 {8, 10, 9, 11, 14, 12, 15, 13}};
5780 for (
unsigned int i = 0; i < 6; ++i)
5782 for (
unsigned int j = 0;
j < 8; ++
j)
5784 cells[i].material_id = 0;
5795 tria.set_all_manifold_ids_on_boundary(0);
5804 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
5806 if (!cell->at_boundary(i))
5812 if (cell->face(i)->center()(0) <
center(0) + 1.e-5 * radius)
5814 cell->face(i)->set_boundary_id(1);
5819 const Point<3> line_vertices[2] = {
5820 cell->face(i)->line(
j)->vertex(0),
5821 cell->face(i)->line(
j)->vertex(1)};
5822 if ((std::fabs(line_vertices[0].distance(
center) - radius) >
5824 (std::fabs(line_vertices[1].distance(
center) - radius) >
5827 cell->face(i)->line(
j)->set_boundary_id(1);
5828 cell->face(i)->line(
j)->set_manifold_id(
5845 const double radius)
5865 for (
unsigned int v = 0; v <
tria_copy.n_vertices(); ++v)
5873 else if (
round == 1)
5875 for (
unsigned int v = 0; v <
tria_copy.n_vertices(); ++v)
5884 else if (
round == 2)
5885 for (
unsigned int v = 0; v <
tria_copy.n_vertices(); ++v)
5898 std::vector<CellData<dim>> cells;
5900 for (
const auto &cell :
tria_copy.cell_iterators())
5904 data.
vertices[v] = cell->vertex_index(v);
5905 data.material_id = cell->material_id();
5906 data.manifold_id = cell->manifold_id();
5907 cells.push_back(data);
5915 if (
round == dim - 1)
5924 for (
const auto &cell :
tria.cell_iterators())
5925 if (cell->
center().norm_square() > 0.4 * radius)
5926 cell->set_manifold_id(1);
5931 tria.set_all_manifold_ids_on_boundary(0);
5949 std::vector<CellData<3>> cells;
5955 static const std::array<Point<3>, 8>
hexahedron = {{{-1, -1, -1},
5965 for (
unsigned int i = 0; i < 8; ++i)
5967 for (
unsigned int i = 0; i < 8; ++i)
5970 const unsigned int n_cells = 6;
5972 {8, 9, 10, 11, 0, 1, 2, 3},
5973 {9, 11, 1, 3, 13, 15, 5, 7},
5974 {12, 13, 4, 5, 14, 15, 6, 7},
5975 {8, 0, 10, 2, 12, 4, 14, 6},
5976 {8, 9, 0, 1, 12, 13, 4, 5},
5977 {10, 2, 11, 3, 14, 6, 15, 7}};
5981 for (
unsigned int i = 0; i <
n_cells; ++i)
5985 cells[i].material_id = 0;
5989 tria.set_all_manifold_ids(0);
6000 std::vector<CellData<3>> cells;
6009 static const std::array<Point<3>, 6>
octahedron = {{{-1, 0, 0},
6017 static const std::array<Point<3>, 8>
hexahedron = {{{-1, -1, -1},
6026 for (
unsigned int i = 0; i < 8; ++i)
6028 for (
unsigned int i = 0; i < 6; ++i)
6030 for (
unsigned int i = 0; i < 8; ++i)
6032 for (
unsigned int i = 0; i < 6; ++i)
6035 const unsigned int n_cells = 12;
6051 for (
unsigned int i = 0; i <
n_cells; ++i)
6053 for (
unsigned int j = 0;
j < 4; ++
j)
6055 cells[i].vertices[
j] =
rhombi[i][
j];
6056 cells[i].vertices[
j + 4] =
rhombi[i][
j] + 14;
6058 cells[i].material_id = 0;
6062 tria.set_all_manifold_ids(0);
6068 const unsigned int n,
6069 const unsigned int n_refinement_steps,
6100 for (
unsigned int r = 0; r < n_refinement_steps; ++r)
6102 tmp.refine_global(1);
6107 for (
const auto &cell : tmp.active_cell_iterators())
6111 if ((cell->vertex(v) - p).norm_square() <
6119 if (r == n_refinement_steps - 1)
6129 tmp = std::move(copy);
6130 tmp.set_all_manifold_ids(0);
6134 tria.set_all_manifold_ids(0);
6149 const unsigned int n_cells,
6155 unsigned int n_refinement_steps = 0;
6157 if (n_cells != 96 && n_cells > 12)
6160 ++n_refinement_steps;
6163 Assert(n_cells == 0 || n_cells == 6 || n_cells == 12 || n_cells == 96 ||
6164 (n_refinement_steps > 0 &&
6166 ExcMessage(
"Invalid number of coarse mesh cells"));
6168 const unsigned int n = n_refinement_steps > 0 ?
6182 internal::hyper_shell_24_48(
6193 tmp.refine_global(1);
6195 tria.set_all_manifold_ids(0);
6221 const unsigned int ,
6231 const double b = a / 2.0;
6232 const double c =
d / 2.0;
6258 {0, 2, 8, 10, 6, 4, 14, 12},
6259 {1, 7, 9, 15, 3, 5, 11, 13},
6260 {6, 4, 14, 12, 7, 5, 15, 13},
6261 {8, 10, 9, 11, 14, 12, 15, 13}};
6265 for (
unsigned int i = 0; i < 5; ++i)
6267 for (
unsigned int j = 0;
j < 8; ++
j)
6269 cells[i].material_id = 0;
6283 for (; cell !=
tria.
end(); ++cell)
6284 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6285 if (cell->at_boundary(i))
6286 cell->face(i)->set_all_boundary_ids(2);
6293 for (
const unsigned int i :
GeometryInfo<3>::face_indices())
6294 if (cell->at_boundary(i))
6304 face->set_all_boundary_ids(0);
6306 face->set_all_boundary_ids(1);
6310 tria.set_all_manifold_ids(0);
6322 const unsigned int n,
6327 if (n == 0 || n == 3)
6354 {0, 1, 3, 2, 4, 5, 7, 6},
6355 {1, 8, 2, 9, 5, 10, 6, 11},
6356 {4, 5, 7, 6, 12, 10, 13, 11},
6358 std::vector<CellData<3>> cells(3);
6360 for (
unsigned int i = 0; i < 3; ++i)
6362 for (
unsigned int j = 0;
j < 8; ++
j)
6364 cells[i].material_id = 0;
6379 tria.set_all_manifold_ids(0);
6388 const double length,
6408 const unsigned int N_r =
6413 const unsigned int N_z =
6415 static_cast<unsigned int>(std::ceil(
6426 for (
unsigned int i = 0; i <
N_r; ++i)
6436 for (
unsigned int j = 0;
j <=
N_z; ++
j)
6437 for (
unsigned int i = 0; i < 2 *
N_r; ++i)
6447 for (
unsigned int j = 0;
j <
N_z; ++
j)
6448 for (
unsigned int i = 0; i <
N_r; ++i)
6450 cells[i +
j *
N_r].vertices[0] = i + (
j + 1) * 2 *
N_r;
6451 cells[i +
j *
N_r].vertices[1] = (i + 1) %
N_r + (
j + 1) * 2 *
N_r;
6452 cells[i +
j *
N_r].vertices[2] = i +
j * 2 *
N_r;
6453 cells[i +
j *
N_r].vertices[3] = (i + 1) %
N_r +
j * 2 *
N_r;
6455 cells[i +
j *
N_r].vertices[4] =
N_r + i + (
j + 1) * 2 *
N_r;
6456 cells[i +
j *
N_r].vertices[5] =
6458 cells[i +
j *
N_r].vertices[6] =
N_r + i +
j * 2 *
N_r;
6459 cells[i +
j *
N_r].vertices[7] =
N_r + ((i + 1) %
N_r) +
j * 2 *
N_r;
6461 cells[i +
j *
N_r].material_id = 0;
6465 tria.set_all_manifold_ids(0);
6471 template <
int dim,
int spacedim>
6480 std::vector<Point<spacedim>>
vertices;
6481 std::vector<CellData<dim>> cells;
6488 ExcMessage(
"The input triangulations must be non-empty "
6489 "and must not be refined."));
6504 cells.push_back(cell_data);
6546 if (std::all_of(cells.begin(), cells.end(), [](
const auto &cell) {
6547 return cell.vertices.size() ==
6548 ReferenceCells::get_hypercube<dim>().n_vertices();
6559 for (
const auto &cell :
tria->cell_iterators())
6561 for (
auto const &f : cell->face_indices())
6577 template <
int dim,
int spacedim>
6629 template <
int structdim>
6633 static_assert(structdim == 1 || structdim == 2,
6634 "This function is only implemented for lines and "
6642 std::sort(std::begin(cell_data.vertices),
6643 std::end(cell_data.vertices));
6644 else if (structdim == 2)
6648 std::copy(std::begin(cell_data.vertices),
6649 std::end(cell_data.vertices),
6679 std::begin(cell_data.vertices));
6686 return std::lexicographical_compare(std::begin(a.vertices),
6687 std::end(a.vertices),
6688 std::begin(
b.vertices),
6689 std::end(
b.vertices));
6703 if (left + 1 != right)
6704 for (
auto it = left;
it != right; ++
it)
6707 Assert(
it->manifold_id == left->manifold_id,
6709 "In the process of grid generation a single "
6710 "line or quadrilateral has been assigned two "
6711 "different manifold ids. This can happen when "
6712 "a Triangulation is copied, e.g., via "
6713 "GridGenerator::replicate_triangulation() and "
6714 "not all external boundary faces have the same "
6715 "manifold id. Double check that all faces "
6716 "which you expect to be merged together have "
6717 "the same manifold id."));
6729 template <
int dim,
int spacedim>
6732 const std::vector<unsigned int> &
extents,
6739 ExcMessage(
"The Triangulation must be copied at least one time in "
6740 "each coordinate dimension."));
6743 const auto &
min =
bbox.get_boundary_points().first;
6744 const auto &
max =
bbox.get_boundary_points().second;
6746 std::array<Tensor<1, spacedim>, dim>
offsets;
6747 for (
unsigned int d = 0;
d < dim; ++
d)
6748 offsets[d][d] = max[d] - min[d];
6752 for (
unsigned int d = 0;
d < dim; ++
d)
6781 for (
unsigned int &vertex :
6789 for (
unsigned int &vertex :
6802 1e-6 * input.begin_active()->diameter());
6832 template <
int dim,
int spacedim>
6840 ExcMessage(
"The two input triangulations are not derived from "
6841 "the same coarse mesh as required."));
6848 ExcMessage(
"The source triangulations for this function must both "
6849 "be available entirely locally, and not be distributed "
6850 "triangulations."));
6884 result.execute_coarsening_and_refinement();
6890 template <
int dim,
int spacedim>
6904 std::vector<CellData<dim>> cells;
6908 Assert(
static_cast<unsigned int>(cell->level()) ==
6911 "Your input triangulation appears to have "
6912 "adaptively refined cells. This is not allowed. You can "
6913 "only call this function on a triangulation in which "
6914 "all cells are on the same refinement level."));
6919 this_cell.material_id = cell->material_id();
6943 const double height,
6948 Assert(input.n_levels() == 1,
6950 "The input triangulation must be a coarse mesh, i.e., it must "
6951 "not have been refined."));
6953 ExcMessage(
"The output triangulation object needs to be empty."));
6955 ExcMessage(
"The given height for extrusion must be positive."));
6958 "The number of slices for extrusion must be at least 2."));
6962 for (
unsigned int i = 0; i <
n_slices; ++i)
6974 const double height,
6988 "GridTools::extrude_triangulation() is only available "
6989 "for Triangulation<3, 3> as output triangulation."));
7002 Assert(input.n_levels() == 1,
7004 "The input triangulation must be a coarse mesh, i.e., it must "
7005 "not have been refined."));
7007 ExcMessage(
"The output triangulation object needs to be empty."));
7010 "The number of slices for extrusion must be at least 2."));
7012 ExcMessage(
"Slice z-coordinates should be in ascending order"));
7014 const auto priorities = [&]() -> std::vector<types::manifold_id> {
7030 "The given vector of manifold ids may not contain any "
7031 "duplicated entries."));
7033 input.get_manifold_ids();
7038 message <<
"The given triangulation has manifold ids {";
7043 <<
" the given vector of manifold ids is {";
7049 <<
" These vectors should contain the same elements.\n";
7050 const std::string m =
message.str();
7059 input.get_manifold_ids();
7064 return dynamic_cast<const TransfiniteInterpolationManifold<2, 2> *>(
7065 &input.get_manifold(id)) == nullptr;
7074 std::vector<Point<3>> points(
n_slices * input.n_vertices());
7075 std::vector<CellData<3>> cells;
7076 cells.reserve((
n_slices - 1) * input.n_active_cells());
7093 for (
const auto &cell : input.active_cell_iterators())
7106 (
slice_n + 1) * input.n_vertices();
7109 this_cell.material_id = cell->material_id();
7111 this_cell.manifold_id = cell->manifold_id();
7119 std::vector<CellData<2>> &quads =
subcell_data.boundary_quads;
7121 quads.reserve(input.n_active_lines() * (
n_slices - 1) +
7122 input.n_active_cells() * 2);
7123 for (
const auto &face : input.active_face_iterators())
7126 quad.boundary_id = face->boundary_id();
7127 if (face->at_boundary())
7130 quad.manifold_id = face->manifold_id();
7134 face->vertex_index(0) +
slice_n * input.n_vertices();
7136 face->vertex_index(1) +
slice_n * input.n_vertices();
7138 face->vertex_index(0) + (
slice_n + 1) * input.n_vertices();
7140 face->vertex_index(1) + (
slice_n + 1) * input.n_vertices();
7141 quads.push_back(quad);
7148 for (
const auto &cell : input.active_cell_iterators())
7152 quad.manifold_id = cell->manifold_id();
7156 cell->vertex_index(0) +
slice_n * input.n_vertices();
7158 cell->vertex_index(1) +
slice_n * input.n_vertices();
7160 cell->vertex_index(2) +
slice_n * input.n_vertices();
7162 cell->vertex_index(3) +
slice_n * input.n_vertices();
7163 quads.push_back(quad);
7174 "The input triangulation to this function is using boundary "
7175 "indicators in a range that do not allow using "
7176 "max_boundary_id+1 and max_boundary_id+2 as boundary "
7177 "indicators for the bottom and top faces of the "
7178 "extruded triangulation."));
7181 for (
const auto &cell : input.active_cell_iterators())
7185 quad.vertices[0] = cell->vertex_index(0);
7186 quad.vertices[1] = cell->vertex_index(1);
7187 quad.vertices[2] = cell->vertex_index(2);
7188 quad.vertices[3] = cell->vertex_index(3);
7190 quad.manifold_id = cell->manifold_id();
7191 quads.push_back(quad);
7194 for (
unsigned int &vertex : quad.
vertices)
7195 vertex += (
n_slices - 1) * input.n_vertices();
7197 quad.manifold_id = cell->manifold_id();
7198 quads.push_back(quad);
7214 for (
const auto &face :
result.active_face_iterators())
7240 "GridTools::extrude_triangulation() is only available "
7241 "for Triangulation<3, 3> as output triangulation."));
7272 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7282 for (; cell !=
endc; ++cell)
7285 if (cell->face(f)->at_boundary())
7290 unsigned int vv = cell->face(f)->vertex_index(v);
7297 cell->face(f)->vertex(v) =
7301 cell->face(f)->vertex(v) =
7305 cell->face(f)->vertex(v) =
7309 cell->face(f)->vertex(v) =
7321 for (; cell !=
endc; ++cell)
7324 if (cell->face(f)->at_boundary())
7326 double dx = cell->face(f)->center()(0) -
center(0);
7327 double dy = cell->face(f)->center()(1) -
center(1);
7331 cell->face(f)->set_boundary_id(0);
7333 cell->face(f)->set_boundary_id(1);
7335 cell->face(f)->set_boundary_id(2);
7337 cell->face(f)->set_boundary_id(3);
7340 cell->face(f)->set_boundary_id(4);
7341 cell->face(f)->set_manifold_id(0);
7346 double d = (cell->face(f)->center() -
center).
norm();
7349 cell->face(f)->set_boundary_id(1);
7350 cell->face(f)->set_manifold_id(0);
7353 cell->face(f)->set_boundary_id(0);
7370 const unsigned int n_cells,
7377 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7381 std::vector<double> radii;
7405 n_cells == 0 ? (dim == 2 ? 8 : 12) :
7432 100.0 * std::numeric_limits<double>::epsilon();
7436 const double radius) {
7451 for (
const auto &cell :
triangulation.active_cell_iterators())
7454 auto face = cell->face(
face_n);
7455 if (face->at_boundary())
7462 face->set_all_boundary_ids(0);
7468 face->set_all_boundary_ids(1);
7482 const unsigned int Nz,
7488 ExcMessage(
"outer_radius has to be bigger than inner_radius."));
7499 for (; cell !=
endc; ++cell)
7502 if (cell->face(f)->at_boundary())
7507 unsigned int vv = cell->face(f)->vertex_index(v);
7511 for (
unsigned int i = 0; i <=
Nz; ++i)
7513 double d = i * L /
Nz;
7514 switch (
vv - i * 16)
7517 cell->face(f)->vertex(v) =
7521 cell->face(f)->vertex(v) =
7525 cell->face(f)->vertex(v) =
7529 cell->face(f)->vertex(v) =
7542 for (; cell !=
endc; ++cell)
7545 if (cell->face(f)->at_boundary())
7547 double dx = cell->face(f)->center()(0);
7548 double dy = cell->face(f)->center()(1);
7549 double dz = cell->face(f)->center()(2);
7554 cell->face(f)->set_boundary_id(0);
7557 cell->face(f)->set_boundary_id(1);
7560 cell->face(f)->set_boundary_id(2);
7563 cell->face(f)->set_boundary_id(3);
7566 cell->face(f)->set_boundary_id(4);
7569 cell->face(f)->set_boundary_id(5);
7573 cell->face(f)->set_all_boundary_ids(6);
7574 cell->face(f)->set_all_manifold_ids(0);
7581 double d = c.norm();
7584 cell->face(f)->set_all_boundary_ids(1);
7585 cell->face(f)->set_all_manifold_ids(0);
7588 cell->face(f)->set_boundary_id(0);
7597 template <
int dim,
int spacedim1,
int spacedim2>
7606 "This function cannot be used on "
7607 "parallel::distributed::Triangulation objects as inputs."));
7609 ExcMessage(
"This function does not work for meshes that have "
7620 for (
unsigned int d = 0;
d < spacedim; ++
d)
7623 std::vector<CellData<dim>> cells(
in_tria.n_active_cells());
7624 for (
const auto &cell :
in_tria.active_cell_iterators())
7626 const unsigned int id = cell->active_cell_index();
7628 cells[id].vertices.resize(cell->n_vertices());
7630 cells[id].
vertices[i] = cell->vertex_index(i);
7631 cells[id].material_id = cell->material_id();
7632 cells[id].manifold_id = cell->manifold_id();
7649 .clear_user_flags_line();
7653 for (
const auto &face :
in_tria.active_face_iterators())
7655 if (face->at_boundary())
7679 if ((face->user_flag_set() ==
false) &&
7694 face->set_user_flag();
7710 .clear_user_flags_line();
7715 .clear_user_flags_quad();
7719 for (
const auto &face :
in_tria.active_face_iterators())
7721 if (face->at_boundary())
7739 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7740 if (face->line(e)->user_flag_set() ==
false)
7744 edge = face->line(e);
7756 edge->set_user_flag();
7774 if (face->user_flag_set() ==
false)
7790 face->set_user_flag();
7798 for (
unsigned int e = 0;
e < face->n_lines(); ++
e)
7799 if (face->line(e)->at_boundary() ==
false)
7800 if (face->line(e)->user_flag_set() ==
false)
7816 edge->set_user_flag();
7838 template <
int dim,
int spacedim>
7846 if (
in_tria.n_global_levels() > 1)
7878 {{{0, 1, 12, 10}}, {{2, 3, 11, 12}}, {{7, 6, 11, 13}},
7879 {{5, 4, 13, 10}}, {{0, 2, 8, 12}}, {{4, 6, 13, 8}},
7880 {{5, 13, 7, 9}}, {{1, 9, 3, 12}}, {{0, 8, 4, 10}},
7881 {{1, 5, 9, 10}}, {{3, 7, 11, 9}}, {{2, 6, 8, 11}},
7882 {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
7883 {{12, 13, 8, 10}}, {{13, 8, 10, 4}}, {{13, 10, 9, 5}},
7884 {{13, 9, 11, 7}}, {{13, 11, 8, 6}}, {{10, 12, 9, 1}},
7885 {{9, 12, 11, 3}}, {{11, 12, 8, 2}}, {{8, 12, 10, 0}}}};
7894 {{{{1, 5}}, {{5, 3}}}},
7895 {{{{0, 6}}, {{6, 1}}}},
7896 {{{{2, 7}}, {{7, 3}}}}}};
7905 {{{{{0, 4, 8}}, {{4, 8, 6}}, {{8, 6, 2}}, {{0, 2, 8}}}},
7906 {{{{1, 3, 9}}, {{3, 9, 7}}, {{9, 7, 5}}, {{1, 9, 5}}}},
7907 {{{{0, 1, 10}}, {{1, 10, 5}}, {{10, 5, 4}}, {{0, 10, 4}}}},
7908 {{{{2, 3, 11}}, {{3, 11, 7}}, {{11, 7, 6}}, {{2, 11, 6}}}},
7909 {{{{0, 1, 12}}, {{1, 12, 3}}, {{12, 3, 2}}, {{0, 12, 2}}}},
7910 {{{{4, 5, 13}}, {{5, 13, 7}}, {{13, 7, 6}}, {{4, 13, 6}}}}}};
7931 {{{0, 12, 10}}, {{12, 1, 10}}, {{12, 1, 9}}, {{12, 3, 9}},
7932 {{12, 2, 11}}, {{12, 3, 11}}, {{12, 0, 8}}, {{12, 2, 8}},
7933 {{9, 13, 5}}, {{13, 7, 9}}, {{11, 7, 13}}, {{11, 6, 13}},
7934 {{4, 8, 13}}, {{6, 8, 13}}, {{4, 13, 10}}, {{13, 5, 10}},
7935 {{10, 9, 5}}, {{10, 9, 1}}, {{11, 9, 7}}, {{11, 9, 3}},
7936 {{8, 11, 2}}, {{8, 11, 6}}, {{8, 10, 0}}, {{8, 10, 4}},
7937 {{12, 3, 9}}, {{12, 9, 11}}, {{12, 3, 11}}, {{3, 9, 11}},
7938 {{2, 12, 8}}, {{2, 12, 11}}, {{2, 11, 8}}, {{8, 12, 11}},
7939 {{0, 12, 10}}, {{0, 12, 8}}, {{0, 8, 10}}, {{8, 10, 12}},
7940 {{12, 1, 10}}, {{12, 1, 9}}, {{1, 10, 9}}, {{10, 9, 12}},
7941 {{10, 8, 4}}, {{10, 8, 13}}, {{4, 13, 8}}, {{4, 13, 10}},
7942 {{10, 9, 13}}, {{10, 9, 5}}, {{13, 5, 10}}, {{13, 5, 9}},
7943 {{13, 7, 9}}, {{13, 7, 11}}, {{9, 11, 13}}, {{9, 11, 7}},
7944 {{8, 11, 13}}, {{8, 11, 6}}, {{6, 13, 8}}, {{6, 13, 11}},
7945 {{12, 13, 10}}, {{12, 13, 8}}, {{8, 10, 13}}, {{8, 10, 12}},
7946 {{12, 13, 10}}, {{12, 13, 9}}, {{10, 9, 13}}, {{10, 9, 12}},
7947 {{12, 13, 9}}, {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
7948 {{12, 13, 11}}, {{12, 13, 8}}, {{8, 11, 13}}, {{8, 11, 12}}}};
7955 {{{12, 10}}, {{12, 9}}, {{12, 11}}, {{12, 8}}, {{9, 13}}, {{11, 13}},
7956 {{8, 13}}, {{10, 13}}, {{10, 9}}, {{9, 11}}, {{11, 8}}, {{8, 10}},
7957 {{12, 9}}, {{12, 11}}, {{11, 9}}, {{12, 8}}, {{12, 11}}, {{11, 8}},
7958 {{12, 8}}, {{12, 10}}, {{10, 8}}, {{12, 10}}, {{12, 9}}, {{9, 10}},
7959 {{13, 10}}, {{13, 8}}, {{8, 10}}, {{13, 10}}, {{13, 9}}, {{9, 10}},
7960 {{13, 11}}, {{13, 9}}, {{11, 9}}, {{13, 11}}, {{13, 8}}, {{11, 8}},
7961 {{12, 13}}, {{8, 10}}, {{8, 13}}, {{10, 13}}, {{8, 12}}, {{10, 12}},
7962 {{12, 13}}, {{10, 9}}, {{10, 13}}, {{9, 13}}, {{10, 12}}, {{9, 12}},
7963 {{12, 13}}, {{9, 11}}, {{9, 13}}, {{11, 13}}, {{9, 12}}, {{11, 12}},
7964 {{12, 13}}, {{11, 8}}, {{11, 13}}, {{8, 13}}, {{11, 12}}, {{8, 12}}}};
7978 {{{{{4, 8}}, {{6, 8}}, {{0, 8}}, {{2, 8}}}},
7979 {{{{5, 9}}, {{7, 9}}, {{1, 9}}, {{3, 9}}}},
7980 {{{{4, 10}}, {{5, 10}}, {{0, 10}}, {{1, 10}}}},
7981 {{{{6, 11}}, {{7, 11}}, {{2, 11}}, {{3, 11}}}},
7982 {{{{2, 12}}, {{3, 12}}, {{0, 12}}, {{1, 12}}}},
7983 {{{{6, 13}}, {{7, 13}}, {{4, 13}}, {{5, 13}}}}}};
7985 std::vector<Point<spacedim>>
vertices;
7986 std::vector<CellData<dim>> cells;
8012 const auto v_global = cell.vertex_index(v);
8018 vertices.push_back(cell.vertex(v));
8026 for (
const auto f : cell.face_indices())
8028 const auto f_global = cell.face_index(f);
8035 cell.face(f)->center(
true));
8051 vertices.push_back(cell.center(
true));
8082 cell_data.vertices[i] =
8084 cell_data.material_id =
8086 cell_data.manifold_id =
8089 cells.push_back(cell_data);
8191 for (
const auto f : cell.face_indices())
8193 const auto bid = cell.face(f)->boundary_id();
8194 const auto mid = cell.face(f)->manifold_id();
8224 for (
const auto e : cell.line_indices())
8226 auto edge = cell.line(e);
8248 template <
int spacedim>
8259 template <
template <
int,
int>
class MeshType,
int dim,
int spacedim>
8263 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8270 const std::set<types::boundary_id> &boundary_ids)
8290 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>
8300 std::vector<CellData<boundary_dim>> cells;
8302 std::vector<Point<spacedim>>
vertices;
8344 for (
const unsigned int i :
GeometryInfo<dim>::face_indices())
8349 if (face->at_boundary() &&
8350 (boundary_ids.empty() ||
8351 (boundary_ids.find(face->boundary_id()) != boundary_ids.end())))
8355 for (
const unsigned int j :
8358 const unsigned int v_index = face->vertex_index(
j);
8371 c_data.manifold_id = face->manifold_id();
8381 for (
unsigned int e = 0;
e < 4; ++
e)
8414 edge.boundary_id = 0;
8415 edge.manifold_id = face->line(e)->manifold_id();
8434 for (
const auto &cell :
surface_mesh.active_cell_iterators())
8435 for (unsigned
int vertex = 0; vertex < 2; ++vertex)
8436 if (cell->face(vertex)->at_boundary())
8437 cell->face(vertex)->set_boundary_id(0);
8445 std::vector<std::pair<
8446 const typename MeshType<dim - 1, spacedim>::cell_iterator,
8447 std::pair<typename MeshType<dim, spacedim>::face_iterator,
unsigned int>>>
8449 for (
const auto &cell :
surface_mesh.active_cell_iterators())
8480 .
second.first->has_children())
8485 .
second.first->refinement_case() ==
8489 .first->set_refine_flag();
8502 .execute_coarsening_and_refinement();
8508 const typename MeshType<dim - 1, spacedim>::cell_iterator
8517 for (
unsigned int child_n = 0;
8537 std::map<
typename MeshType<dim - 1, spacedim>::cell_iterator,
8549 template <
int dim,
int spacedim>
8562 std::vector<Point<spacedim>>
vertices;
8563 std::vector<CellData<dim>> cells;
8573 for (
unsigned int i = 0; i <=
repetitions[0]; ++i)
8582 std::array<unsigned int, 4> quad{{
8592 tri.vertices = {quad[0], quad[1], quad[2]};
8593 cells.push_back(
tri);
8599 tri.vertices = {quad[3], quad[2], quad[1]};
8600 cells.push_back(
tri);
8614 for (
unsigned int i = 0; i <=
repetitions[0]; ++i)
8617 p1[2] + dx[2] *
k));
8625 std::array<unsigned int, 8> quad{
8646 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
8647 cell.vertices = {{quad[0], quad[1], quad[2], quad[4]}};
8649 cell.vertices = {{quad[0], quad[1], quad[3], quad[5]}};
8651 cells.push_back(cell);
8657 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
8658 cell.vertices = {{quad[2], quad[1], quad[3], quad[7]}};
8660 cell.vertices = {{quad[0], quad[3], quad[2], quad[6]}};
8661 cells.push_back(cell);
8667 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
8668 cell.vertices = {{quad[1], quad[4], quad[5], quad[7]}};
8670 cell.vertices = {{quad[0], quad[4], quad[5], quad[6]}};
8671 cells.push_back(cell);
8677 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
8678 cell.vertices = {{quad[2], quad[4], quad[7], quad[6]}};
8680 cell.vertices = {{quad[3], quad[5], quad[7], quad[6]}};
8681 cells.push_back(cell);
8687 if (((i % 2) + (
j % 2) + (
k % 2)) % 2 == 0)
8688 cell.vertices = {{quad[1], quad[2], quad[4], quad[7]}};
8690 cell.vertices = {{quad[0], quad[3], quad[6], quad[5]}};
8691 cells.push_back(cell);
8706 template <
int dim,
int spacedim>
8736# include "grid_generator.inst"
value_type * data() const noexcept
void add_parameter(const std::string &entry, ParameterType ¶meter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void enter_subsection(const std::string &subsection)
numbers::NumberTraits< Number >::real_type 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)
const Mapping< dim, spacedim > & get_default_linear_mapping(const Triangulation< dim, spacedim > &triangulation)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void create_triangulation(Triangulation< dim, dim > &tria, const AdditionalData &additional_data=AdditionalData())
void subdivided_hyper_cube_with_simplices(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double p1=0.0, const double p2=1.0, const bool colorize=false)
void parallelepiped(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void hyper_cross(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &sizes, const bool colorize_cells=false)
A center cell with stacks of cell protruding from each surface.
void hyper_ball_balanced(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void plate_with_a_hole(Triangulation< dim > &tria, const double inner_radius=0.4, const double outer_radius=1., const double pad_bottom=2., const double pad_top=2., const double pad_left=1., const double pad_right=1., const Point< dim > ¢er=Point< dim >(), const types::manifold_id polar_manifold_id=0, const types::manifold_id tfi_manifold_id=1, const double L=1., const unsigned int n_slices=2, const bool colorize=false)
Rectangular plate with an (offset) cylindrical hole.
void enclosed_hyper_cube(Triangulation< dim > &tria, const double left=0., const double right=1., const double thickness=1., const bool colorize=false)
void replicate_triangulation(const Triangulation< dim, spacedim > &input, const std::vector< unsigned int > &extents, Triangulation< dim, spacedim > &result)
Replicate a given triangulation in multiple coordinate axes.
void parallelogram(Triangulation< dim > &tria, const Point< dim >(&corners)[dim], const bool colorize=false)
void general_cell(Triangulation< dim, spacedim > &tria, const std::vector< Point< spacedim > > &vertices, const bool colorize=false)
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void hyper_L(Triangulation< dim > &tria, const double left=-1., const double right=1., const bool colorize=false)
void hyper_cube_slit(Triangulation< dim > &tria, const double left=0., const double right=1., const bool colorize=false)
void hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void eccentric_hyper_shell(Triangulation< dim > &triangulation, const Point< dim > &inner_center, const Point< dim > &outer_center, const double inner_radius, const double outer_radius, const unsigned int n_cells)
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
void moebius(Triangulation< 3, 3 > &tria, const unsigned int n_cells, const unsigned int n_rotations, const double R, const double r)
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
void half_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void quarter_hyper_ball(Triangulation< dim > &tria, const Point< dim > ¢er=Point< dim >(), const double radius=1.)
void cylinder_shell(Triangulation< dim > &tria, const double length, const double inner_radius, const double outer_radius, const unsigned int n_radial_cells=0, const unsigned int n_axial_cells=0)
void cheese(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &holes)
Rectangular domain with rectangular pattern of holes.
void create_union_triangulation(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result)
void subdivided_hyper_rectangle_with_simplices(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void non_standard_orientation_mesh(Triangulation< 2 > &tria, const unsigned int n_rotate_middle_square)
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 subdivided_parallelepiped(Triangulation< dim > &tria, const unsigned int n_subdivisions, const Point< dim >(&corners)[dim], const bool colorize=false)
void subdivided_cylinder(Triangulation< dim > &tria, const unsigned int x_subdivisions, const double radius=1., const double half_length=1.)
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
void subdivided_hyper_L(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &bottom_left, const Point< dim > &top_right, const std::vector< int > &n_cells_to_remove)
void hyper_sphere(Triangulation< spacedim - 1, spacedim > &tria, const Point< spacedim > ¢er=Point< spacedim >(), const double radius=1.)
void concentric_hyper_shells(Triangulation< dim > &triangulation, const Point< dim > ¢er, const double inner_radius=0.125, const double outer_radius=0.25, const unsigned int n_shells=1, const double skewness=0.1, const unsigned int n_cells_per_shell=0, const bool colorize=false)
void convert_hypercube_to_simplex_mesh(const Triangulation< dim, spacedim > &in_tria, Triangulation< dim, spacedim > &out_tria)
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void quarter_hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, const bool colorize=false)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void hyper_shell(Triangulation< dim > &tria, const Point< dim > ¢er, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim > > &vertices)
void hyper_cube_with_cylindrical_hole(Triangulation< dim > &triangulation, const double inner_radius=.25, const double outer_radius=.5, const double L=.5, const unsigned int repetitions=1, const bool colorize=false)
void truncated_cone(Triangulation< dim > &tria, const double radius_0=1.0, const double radius_1=0.5, const double half_length=1.0)
void 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 flatten_triangulation(const Triangulation< dim, spacedim1 > &in_tria, Triangulation< dim, spacedim2 > &out_tria)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
long double gamma(const unsigned int n)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
void copy(const T *begin, const T *end, U *dest)
const types::material_id invalid_material_id
static constexpr double PI_2
const types::boundary_id invalid_boundary_id
static constexpr double PI
const types::boundary_id internal_face_boundary_id
static const unsigned int invalid_unsigned_int
const types::manifold_id flat_manifold_id
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
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 > &)
const ::Triangulation< dim, spacedim > & tria