template <>
void
torus<2, 3>(Triangulation<2, 3> &tria,
- const double R,
- const double r,
+ const double centerline_radius,
+ const double inner_radius,
const unsigned int,
const double)
{
- Assert(R > r,
- ExcMessage("Outer radius R must be greater than the inner "
- "radius r."));
- Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
+ Assert(centerline_radius > inner_radius,
+ ExcMessage("The centerline radius must be greater than the "
+ "inner radius."));
+ Assert(inner_radius > 0.0,
+ ExcMessage("The inner radius must be positive."));
const unsigned int dim = 2;
const unsigned int spacedim = 3;
std::vector<Point<spacedim>> vertices(16);
- vertices[0] = Point<spacedim>(R - r, 0, 0);
- vertices[1] = Point<spacedim>(R, -r, 0);
- vertices[2] = Point<spacedim>(R + r, 0, 0);
- vertices[3] = Point<spacedim>(R, r, 0);
- vertices[4] = Point<spacedim>(0, 0, R - r);
- vertices[5] = Point<spacedim>(0, -r, R);
- vertices[6] = Point<spacedim>(0, 0, R + r);
- vertices[7] = Point<spacedim>(0, r, R);
- vertices[8] = Point<spacedim>(-(R - r), 0, 0);
- vertices[9] = Point<spacedim>(-R, -r, 0);
- vertices[10] = Point<spacedim>(-(R + r), 0, 0);
- vertices[11] = Point<spacedim>(-R, r, 0);
- vertices[12] = Point<spacedim>(0, 0, -(R - r));
- vertices[13] = Point<spacedim>(0, -r, -R);
- vertices[14] = Point<spacedim>(0, 0, -(R + r));
- vertices[15] = Point<spacedim>(0, r, -R);
+ vertices[0] = Point<spacedim>(centerline_radius - inner_radius, 0, 0);
+ vertices[1] = Point<spacedim>(centerline_radius, -inner_radius, 0);
+ vertices[2] = Point<spacedim>(centerline_radius + inner_radius, 0, 0);
+ vertices[3] = Point<spacedim>(centerline_radius, inner_radius, 0);
+ vertices[4] = Point<spacedim>(0, 0, centerline_radius - inner_radius);
+ vertices[5] = Point<spacedim>(0, -inner_radius, centerline_radius);
+ vertices[6] = Point<spacedim>(0, 0, centerline_radius + inner_radius);
+ vertices[7] = Point<spacedim>(0, inner_radius, centerline_radius);
+ vertices[8] = Point<spacedim>(-(centerline_radius - inner_radius), 0, 0);
+ vertices[9] = Point<spacedim>(-centerline_radius, -inner_radius, 0);
+ vertices[10] = Point<spacedim>(-(centerline_radius + inner_radius), 0, 0);
+ vertices[11] = Point<spacedim>(-centerline_radius, inner_radius, 0);
+ vertices[12] = Point<spacedim>(0, 0, -(centerline_radius - inner_radius));
+ vertices[13] = Point<spacedim>(0, -inner_radius, -centerline_radius);
+ vertices[14] = Point<spacedim>(0, 0, -(centerline_radius + inner_radius));
+ vertices[15] = Point<spacedim>(0, inner_radius, -centerline_radius);
std::vector<CellData<dim>> cells(16);
// Right Hand Orientation
tria.create_triangulation(vertices, cells, SubCellData());
tria.set_all_manifold_ids(0);
- tria.set_manifold(0, TorusManifold<2>(R, r));
+ tria.set_manifold(0, TorusManifold<2>(centerline_radius, inner_radius));
}
template <>
void
torus<3, 3>(Triangulation<3, 3> &tria,
- const double R,
- const double r,
+ const double centerline_radius,
+ const double inner_radius,
const unsigned int n_cells_toroidal,
const double phi)
{
- Assert(R > r,
- ExcMessage("Outer radius R must be greater than the inner "
- "radius r."));
- Assert(r > 0.0, ExcMessage("The inner radius r must be positive."));
+ Assert(centerline_radius > inner_radius,
+ ExcMessage("The centerline radius must be greater than the "
+ "inner radius."));
+ Assert(inner_radius > 0.0,
+ ExcMessage("The inner radius must be positive."));
Assert(n_cells_toroidal > static_cast<unsigned int>(phi / numbers::PI),
ExcMessage("Number of cells in toroidal direction has "
"to be at least 3 for a torus of polar extent 2*pi."));
ExcMessage("Invalid angle phi specified."));
// the first 8 vertices are in the x-y-plane
- const Point<3> p = Point<3>(R, 0.0, 0.0);
+ const Point<3> p = Point<3>(centerline_radius, 0.0, 0.0);
const double a = 1. / (1 + std::sqrt(2.0));
// A value of 1 indicates "open" torus with angle < 2*pi, which
// means that we need an additional layer of vertices
const unsigned int n_point_layers_toroidal =
n_cells_toroidal + additional_layer;
std::vector<Point<3>> vertices(8 * n_point_layers_toroidal);
- vertices[0] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0)),
- vertices[1] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0)),
- vertices[2] = p + Point<3>(-1, -1, 0) * (r / std::sqrt(2.0) * a),
- vertices[3] = p + Point<3>(+1, -1, 0) * (r / std::sqrt(2.0) * a),
- vertices[4] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0) * a),
- vertices[5] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0) * a),
- vertices[6] = p + Point<3>(-1, +1, 0) * (r / std::sqrt(2.0)),
- vertices[7] = p + Point<3>(+1, +1, 0) * (r / std::sqrt(2.0));
+ vertices[0] = p + Point<3>(-1, -1, 0) * (inner_radius / std::sqrt(2.0)),
+ vertices[1] = p + Point<3>(+1, -1, 0) * (inner_radius / std::sqrt(2.0)),
+ vertices[2] = p + Point<3>(-1, -1, 0) * (inner_radius / std::sqrt(2.0) * a),
+ vertices[3] = p + Point<3>(+1, -1, 0) * (inner_radius / std::sqrt(2.0) * a),
+ vertices[4] = p + Point<3>(-1, +1, 0) * (inner_radius / std::sqrt(2.0) * a),
+ vertices[5] = p + Point<3>(+1, +1, 0) * (inner_radius / std::sqrt(2.0) * a),
+ vertices[6] = p + Point<3>(-1, +1, 0) * (inner_radius / std::sqrt(2.0)),
+ vertices[7] = p + Point<3>(+1, +1, 0) * (inner_radius / std::sqrt(2.0));
// create remaining vertices by rotating around negative y-axis (the
// direction is to ensure positive cell measures)
{
for (unsigned int v = 0; v < 8; ++v)
{
- const double r_2d = vertices[v][0];
- vertices[8 * c + v][0] = r_2d * std::cos(phi_cell * c);
+ const double inner_radius_2d = vertices[v][0];
+ vertices[8 * c + v][0] = inner_radius_2d * std::cos(phi_cell * c);
vertices[8 * c + v][1] = vertices[v][1];
- vertices[8 * c + v][2] = r_2d * std::sin(phi_cell * c);
+ vertices[8 * c + v][2] = inner_radius_2d * std::sin(phi_cell * c);
}
}
}
}
- tria.set_manifold(1, TorusManifold<3>(R, r));
+ tria.set_manifold(1, TorusManifold<3>(centerline_radius, inner_radius));
tria.set_manifold(2,
CylindricalManifold<3>(Tensor<1, 3>({0., 1., 0.}),
Point<3>()));
double z = p[1];
double y = p[2];
double phi = std::atan2(y, x);
- double theta = std::atan2(z, std::sqrt(x * x + y * y) - R);
+ double theta = std::atan2(z, std::sqrt(x * x + y * y) - centerline_radius);
double w =
- std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * R) +
- Utilities::fixed_power<2>(x - std::cos(phi) * R) + z * z) /
- r;
+ std::sqrt(Utilities::fixed_power<2>(y - std::sin(phi) * centerline_radius) +
+ Utilities::fixed_power<2>(x - std::cos(phi) * centerline_radius) +
+ z * z) /
+ inner_radius;
return {phi, theta, w};
}
double theta = chart_point[1];
double w = chart_point[2];
- return {std::cos(phi) * R + r * w * std::cos(theta) * std::cos(phi),
- r * w * std::sin(theta),
- std::sin(phi) * R + r * w * std::cos(theta) * std::sin(phi)};
+ return {std::cos(phi) * centerline_radius +
+ inner_radius * w * std::cos(theta) * std::cos(phi),
+ inner_radius * w * std::sin(theta),
+ std::sin(phi) * centerline_radius +
+ inner_radius * w * std::cos(theta) * std::sin(phi)};
}
template <int dim>
-TorusManifold<dim>::TorusManifold(const double R, const double r)
+TorusManifold<dim>::TorusManifold(const double centerline_radius,
+ const double inner_radius)
: ChartManifold<dim, 3, 3>(Point<3>(2 * numbers::PI, 2 * numbers::PI, 0.0))
- , r(r)
- , R(R)
+ , centerline_radius(centerline_radius)
+ , inner_radius(inner_radius)
{
- Assert(R > r,
- ExcMessage("Outer radius R must be greater than the inner "
- "radius r."));
- Assert(r > 0.0, ExcMessage("inner radius must be positive."));
+ Assert(centerline_radius > inner_radius,
+ ExcMessage("The centerline radius must be greater than the "
+ "inner radius."));
+ Assert(inner_radius > 0.0, ExcMessage("The inner radius must be positive."));
}
std::unique_ptr<Manifold<dim, 3>>
TorusManifold<dim>::clone() const
{
- return std::make_unique<TorusManifold<dim>>(R, r);
+ return std::make_unique<TorusManifold<dim>>(centerline_radius, inner_radius);
}
double theta = chart_point[1];
double w = chart_point[2];
- DX[0][0] = -std::sin(phi) * R - r * w * std::cos(theta) * std::sin(phi);
- DX[0][1] = -r * w * std::sin(theta) * std::cos(phi);
- DX[0][2] = r * std::cos(theta) * std::cos(phi);
+ DX[0][0] = -std::sin(phi) * centerline_radius -
+ inner_radius * w * std::cos(theta) * std::sin(phi);
+ DX[0][1] = -inner_radius * w * std::sin(theta) * std::cos(phi);
+ DX[0][2] = inner_radius * std::cos(theta) * std::cos(phi);
DX[1][0] = 0;
- DX[1][1] = r * w * std::cos(theta);
- DX[1][2] = r * std::sin(theta);
+ DX[1][1] = inner_radius * w * std::cos(theta);
+ DX[1][2] = inner_radius * std::sin(theta);
- DX[2][0] = std::cos(phi) * R + r * w * std::cos(theta) * std::cos(phi);
- DX[2][1] = -r * w * std::sin(theta) * std::sin(phi);
- DX[2][2] = r * std::cos(theta) * std::sin(phi);
+ DX[2][0] = std::cos(phi) * centerline_radius +
+ inner_radius * w * std::cos(theta) * std::cos(phi);
+ DX[2][1] = -inner_radius * w * std::sin(theta) * std::sin(phi);
+ DX[2][2] = inner_radius * std::cos(theta) * std::sin(phi);
return DX;
}