EllipticalManifold(const Point<spacedim> & center,
const Tensor<1, spacedim> &major_axis_direction,
const double c_parameter);
- virtual ~EllipticalManifold() = default;
virtual std::unique_ptr<Manifold<dim, spacedim>>
clone() const override;
+ /**
+ * @copydoc ChartManifold::pull_back()
+ */
virtual Point<spacedim>
pull_back(const Point<spacedim> &space_point) const override;
+ /**
+ * @copydoc ChartManifold::push_forward()
+ */
virtual Point<spacedim>
push_forward(const Point<spacedim> &chart_point) const override;
+ /**
+ * @copydoc ChartManifold::push_forward_gradient()
+ */
virtual DerivativeForm<1, spacedim, spacedim>
push_forward_gradient(const Point<spacedim> &chart_point) const override;
private:
/**
- * Return the periodicity of the coordinate variables of the manifold.
+ * @copydoc ChartManifold::get_periodicity()
*
- * For $dim=2$ and $spacedim=2$, first coordinate is non-periodic, while
- * second coordinate has a periodicity of $2\pi$.
+ * For $\text{dim}=2$ and $\text{spacedim}=2$, the first coordinate is
+ * non-periodic, while the second coordinate has a periodicity of $2\pi$.
*/
static Tensor<1, spacedim>
get_periodicity();
}
+
// ============================================================
// EllipticalManifold
// ============================================================
direction /= direction_norm;
}
+
+
template <int dim, int spacedim>
std::unique_ptr<Manifold<dim, spacedim>>
EllipticalManifold<dim, spacedim>::clone() const
c_parameter);
}
+
+
template <int dim, int spacedim>
Tensor<1, spacedim>
EllipticalManifold<dim, spacedim>::get_periodicity()
{
Tensor<1, spacedim> periodicity;
// The second elliptical coordinate is periodic, while the first is not.
- // Enforce periodicity on last variable.
+ // Enforce periodicity on the last variable.
periodicity[spacedim - 1] = 2.0 * numbers::PI;
return periodicity;
}
+
+
template <int dim, int spacedim>
Point<spacedim>
EllipticalManifold<dim, spacedim>::push_forward(const Point<spacedim> &) const
Assert(false, ExcNotImplemented());
return Point<spacedim>();
}
+
+
+
template <>
Point<2>
EllipticalManifold<2, 2>::push_forward(const Point<2> &chart_point) const
direction[1] * x + direction[0] * y);
return p + center;
}
+
+
+
template <int dim, int spacedim>
Point<spacedim>
EllipticalManifold<dim, spacedim>::pull_back(const Point<spacedim> &) const
Assert(false, ExcNotImplemented());
return Point<spacedim>();
}
+
+
+
template <>
Point<2>
EllipticalManifold<2, 2>::pull_back(const Point<2> &space_point) const
const double y0 = space_point[1] - center[1];
const double x = direction[0] * x0 + direction[1] * y0;
const double y = -direction[1] * x0 + direction[0] * y0;
- // From here we try to find solutions to equation
- // x^2/(1-q)-y^2/q = c^2 for q, and q = -sinh^2(pt[0]),
- // and to equation
- // x^2/(1-p)-y^2/p = c^2 for p, and p = sin^2(pt[1]).
- // Note that p and q are solutions to the same quadratic equation.
+ // From here we try to find solutions of two equations:
+ // x^2/(1-q)-y^2/q = c^2 for q, given q = -sinh^2(pt[0]),
+ // and
+ // x^2/(1-p)-y^2/p = c^2 for p, given p = sin^2(pt[1]).
+ // Note that, in the end, p and q are solutions to the same quadratic
+ // equation.
const double x2 = x * x;
const double y2 = y * y;
//
p = 0.0;
if (p > 1.0)
p = 1.0;
- const double q = (-b - srdelta) / (2.0 * c2);
- const bool x_pos = std::signbit(x) == 0;
- const bool y_pos = std::signbit(y) == 0;
- const double eta0 = std::asin(std::sqrt(p));
+ const double q = (-b - srdelta) / (2.0 * c2);
+ const bool x_is_positive = !std::signbit(x);
+ const bool y_is_positive = !std::signbit(y);
+ const double eta0 = std::asin(std::sqrt(p));
// Given q = -sinh^2(pt[0]),
// pt[0] is calculated by straight inversion since q <= 0 for any
// pt[0].
Point<2> pt(std::log(std::sqrt(-q) + std::sqrt(1.0 - q)), eta0);
// Unfolding pt[1] according to the quadrant.
- if (x_pos && !y_pos)
+ if (x_is_positive && !y_is_positive)
pt[1] = 2.0 * numbers::PI - eta0;
- else if (!x_pos && y_pos)
+ else if (!x_is_positive && y_is_positive)
pt[1] = numbers::PI - eta0;
- else if (!x_pos && !y_pos)
+ else if (!x_is_positive && !y_is_positive)
pt[1] = numbers::PI + eta0;
return pt;
}
+
+
template <int dim, int spacedim>
DerivativeForm<1, spacedim, spacedim>
EllipticalManifold<dim, spacedim>::push_forward_gradient(
- const Point<spacedim> &chart_point) const
+ const Point<spacedim> &) const
{
- DerivativeForm<1, spacedim, spacedim> dX;
- const double ch = std::cosh(chart_point[0]);
- const double sh = std::sinh(chart_point[0]);
- const double cs = std::cos(chart_point[1]);
- const double sn = std::sin(chart_point[1]);
- switch (spacedim)
- {
- case 2:
- {
- dX[0][0] = c_parameter * sh * cs;
- dX[0][1] = -c_parameter * ch * sn;
- dX[1][0] = c_parameter * ch * sn;
- dX[1][1] = c_parameter * sh * cs;
- }
- break;
- default:
- Assert(false, ExcNotImplemented());
- }
+ Assert(false, ExcNotImplemented());
+ return DerivativeForm<1, spacedim, spacedim>();
+}
+
+
+
+template <>
+DerivativeForm<1, 2, 2>
+EllipticalManifold<2, 2>::push_forward_gradient(
+ const Point<2> &chart_point) const
+{
+ const double ch = std::cosh(chart_point[0]);
+ const double sh = std::sinh(chart_point[0]);
+ const double cs = std::cos(chart_point[1]);
+ const double sn = std::sin(chart_point[1]);
+ DerivativeForm<1, 2, 2> dX;
+ dX[0][0] = c_parameter * sh * cs;
+ dX[0][1] = -c_parameter * ch * sn;
+ dX[1][0] = c_parameter * ch * sn;
+ dX[1][1] = c_parameter * sh * cs;
return dX;
}
+
+
// ============================================================
// FunctionManifold
// ============================================================
// hyper_shell made of 8 elements.
namespace
{
+ /* The following system of equations
+ // /
+ // |x = ellPar*cosh(u)*cos(v)
+ // |y = ellPar*sinh(u)*sin(v)
+ // \
+ // allows to transform a point (u,v) in the chart space to a
+ // point (x,y) in the cartesian space.
+ //
+ // Setting v=0, the system of equations becomes
+ // /
+ // |x = x_0 = ellPar*cosh(u)
+ // |y = 0
+ // \
+ // which represents the crossing point of an ellipsis, non-rotated and
+ // centered at the origin of the cartesian system, with its major axis.
+ //
+ // For the sake of completeness, by setting v=pi/2, one finds y_0 =
+ // ellPar*sinh(u) as the crossing point on the minor axis.
+ //
+ // This function takes in input a point pt={x_0,v}
+ // and transforms it into (x,y) in accordance with the first system of
+ // equations above.
+ //
+ // On a final note, the condition x_0 >= ellPar has to be satisfied. This is
+ // required because min(cosh(u)) = 1 is obtained for u = 0, and imposing v = 0
+ // the minimum value that x can assume is x = ellPar.
+ // For the sake of completeness, it's worth to notice that there would be no
+ // restrictions on the input if we passed y_0, instead of x_0, as pt[0].
+ */
Point<2>
chart_to_cartesian(const Point<2> &pt, const double ellPar)
{
c[1] = ellPar * sh * std::sin(pt[1]);
return c;
}
+
+
+
std::vector<Point<2>>
generate_shell_points(const Point<2> ¢er,
const double radius0,
}
return points;
}
+
+
+
+ // Generate an hyper_shell over an EllipticalManifold having an abitrary
+ // center, and the major axis oriented in the direction of the x-axis.
+ //
+ // inner_radius and outer_radius parameters correspond to the
+ // distances from the center of the manifold of the
+ // crossing points of two concentric ellipsis with their major axis.
+ //
+ // Input parameters must respect the following constrains:
+ // c_parameter < inner_radius < outer_radius.
+ //
+ // To understand the constrain on the c_parameter refer to the
+ // documentation of function chart_to_cartesian() whithin this namespace.
void build_simple_hyper_shell(Triangulation<2, 2> &grid,
const Point<2> & center,
const double inner_radius,
grid.set_manifold(0, EllipticalManifold<2, 2>(center, axis, c_param));
grid.set_all_manifold_ids(0);
}
+
+
+
} // namespace
+
+
+
// Helper function
// Generate a simple hyper_shell over an elliptical manifold centered at the
// origin. Major axis is the x-axis.
gridout.write_msh(tria, deallog.get_file_stream());
}
+
+
int
main()
{