FlatManifold<dim> chart_manifold;
};
+#ifndef DOXYGEN
+/**
+ * Specialization for the only properly implemented spacedim parameter.
+ */
+template <int dim>
+class CylindricalManifold<dim,3> : public ChartManifold<dim,3,3>
+{
+public:
+ /**
+ * Constructor. Using default values for the constructor arguments yields a
+ * cylinder along the x-axis (<tt>axis=0</tt>). Choose <tt>axis=1</tt> or
+ * <tt>axis=2</tt> for a tube along the y- or z-axis, respectively. The
+ * tolerance value is used to determine if a point is on the axis.
+ */
+ CylindricalManifold (const unsigned int axis = 0,
+ const double tolerance = 1e-10);
+
+ /**
+ * Constructor. If constructed with this constructor, the manifold described
+ * is a cylinder with an axis that points in direction #direction and goes
+ * through the given #point_on_axis. The direction may be arbitrarily
+ * scaled, and the given point may be any point on the axis. The tolerance
+ * value is used to determine if a point is on the axis.
+ */
+ CylindricalManifold (const Point<3> &direction,
+ const Point<3> &point_on_axis,
+ const double tolerance = 1e-10);
+
+ /**
+ * Compute the Cartesian coordinates for a point given in cylindrical
+ * coordinates.
+ */
+ virtual Point<3>
+ pull_back(const Point<3> &space_point) const override;
+
+ /**
+ * Compute the cylindrical coordinates $(r, \phi, \lambda)$ for the given
+ * point where $r$ denotes the distance from the axis,
+ * $\phi$ the angle between the given point and the computed normal
+ * direction and $\lambda$ the axial position.
+ */
+ virtual Point<3>
+ push_forward(const Point<3> &chart_point) const override;
+
+ /**
+ * Compute the derivatives of the mapping from cylindrical coordinates
+ * $(r, \phi, \lambda)$ to cartesian coordinates where $r$ denotes the
+ * distance from the axis, $\phi$ the angle between the given point and the
+ * computed normal direction and $\lambda$ the axial position.
+ */
+ virtual DerivativeForm<1, 3, 3>
+ push_forward_gradient(const Point<3> &chart_point) const override;
+
+ /**
+ * Compute new points on the CylindricalManifold. See the documentation of
+ * the base class for a detailed description of what this function does.
+ */
+ virtual Point<3>
+ get_new_point (const ArrayView<const Point<3>> &surrounding_points,
+ const ArrayView<const double> &weights) const override;
+
+protected:
+ /**
+ * A vector orthogonal to the normal direction.
+ */
+ const Tensor<1,3> normal_direction;
+
+ /**
+ * The direction vector of the axis.
+ */
+
+ const Tensor<1,3> direction;
+ /**
+ * An arbitrary point on the axis.
+ */
+ const Point<3> point_on_axis;
+
+private:
+ /**
+ * Relative tolerance to measure zero distances.
+ */
+ double tolerance;
+
+};
+#endif //DOXYGEN
+
DEAL_II_NAMESPACE_CLOSE
#endif
template <int dim, int spacedim>
CylindricalManifold<dim, spacedim>::CylindricalManifold(const unsigned int axis,
const double tolerance) :
- CylindricalManifold<dim, spacedim>(Point<3>::unit_vector(axis),
- Point<3>(),
+ CylindricalManifold<dim, spacedim>(Point<spacedim>::unit_vector(axis),
+ Point<spacedim>(),
tolerance)
{}
+template <int dim>
+CylindricalManifold<dim, 3>::CylindricalManifold(const unsigned int axis,
+ const double tolerance) :
+ CylindricalManifold<dim, 3>(Point<3>::unit_vector(axis),
+ Point<3>(),
+ tolerance)
+{}
+
+
+
template <int dim, int spacedim>
CylindricalManifold<dim, spacedim>::CylindricalManifold(const Point<spacedim> &direction_,
const Point<spacedim> &point_on_axis_,
const double tolerance) :
- ChartManifold<dim,3,3>(Tensor<1,3>({0,2.*numbers::PI,0})),
- normal_direction(internal::compute_normal(direction_, true)),
+ ChartManifold<dim,spacedim,3>(Tensor<1,3>({0,2.*numbers::PI,0})),
+ normal_direction(Tensor<1, spacedim>()),
direction (direction_/direction_.norm()),
point_on_axis (point_on_axis_),
tolerance(tolerance)
+template <int dim>
+CylindricalManifold<dim, 3>::CylindricalManifold(const Point<3> &direction_,
+ const Point<3> &point_on_axis_,
+ const double tolerance) :
+ ChartManifold<dim,3,3>(Tensor<1,3>({0,2.*numbers::PI,0})),
+ normal_direction(internal::compute_normal(direction_, true)),
+ direction (direction_/direction_.norm()),
+ point_on_axis (point_on_axis_),
+ tolerance(tolerance)
+{}
+
+
+
template <int dim, int spacedim>
Point<spacedim>
CylindricalManifold<dim,spacedim>::
get_new_point (const ArrayView<const Point<spacedim>> &surrounding_points,
const ArrayView<const double> &weights) const
+{
+ Assert (spacedim==3,
+ ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
+ return Point<spacedim>();
+}
+
+
+
+template <int dim>
+Point<3>
+CylindricalManifold<dim,3>::
+get_new_point (const ArrayView<const Point<3>> &surrounding_points,
+ const ArrayView<const double> &weights) const
{
// First check if the average in space lies on the axis.
- Point<spacedim> middle;
+ Point<3> middle;
double average_length = 0.;
for (unsigned int i=0; i<surrounding_points.size(); ++i)
{
const double lambda = middle*direction;
if ((middle-direction*lambda).square() < tolerance*average_length)
- return Point<spacedim>()+direction*lambda;
+ return Point<3>()+direction*lambda;
else // If not, using the ChartManifold should yield valid results.
- return ChartManifold<dim, spacedim, 3>::get_new_point(surrounding_points,
- weights);
+ return ChartManifold<dim, 3, 3>::get_new_point(surrounding_points,
+ weights);
}
template <int dim, int spacedim>
Point<3>
CylindricalManifold<dim, spacedim>::pull_back(const Point<spacedim> &space_point) const
+{
+ Assert (spacedim==3,
+ ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
+ return Point<3>();
+}
+
+
+
+template <int dim>
+Point<3>
+CylindricalManifold<dim, 3>::pull_back(const Point<3> &space_point) const
{
// First find the projection of the given point to the axis.
const Tensor<1,3> normalized_point = space_point-point_on_axis;
template <int dim, int spacedim>
Point<spacedim>
CylindricalManifold<dim, spacedim>::push_forward(const Point<3> &chart_point) const
+{
+ Assert (spacedim==3,
+ ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
+ return Point<spacedim>();
+}
+
+
+
+template<int dim>
+Point<3>
+CylindricalManifold<dim, 3>::push_forward(const Point<3> &chart_point) const
{
// Rotate the orthogonal direction by the given angle.
// Formula from Section 5.2 in
template<int dim, int spacedim>
DerivativeForm<1, 3, spacedim>
CylindricalManifold<dim, spacedim>::push_forward_gradient(const Point<3> &chart_point) const
+{
+ Assert (spacedim==3,
+ ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
+
+ return DerivativeForm<1,3,spacedim>();
+}
+
+
+
+template<int dim>
+DerivativeForm<1, 3, 3>
+CylindricalManifold<dim, 3>::push_forward_gradient(const Point<3> &chart_point) const
{
Tensor<2, 3> derivatives;
projects / dealii.git / commitdiff