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
}
// helper function to compute a vector orthogonal to a given one.
+ // does nothing unless spacedim == 3.
+ template <int spacedim>
+ Point<spacedim>
+ compute_normal(const Tensor<1, spacedim> &/*vector*/,
+ bool /*normalize*/=false)
+ {
+ return {};
+ }
+
Point<3>
compute_normal(const Tensor<1,3> &vector, bool normalize=false)
{
-template<int dim, int spacedim>
+template <int dim, int spacedim>
Tensor<1, spacedim>
SphericalManifold<dim, spacedim>::
normal_vector (const typename Triangulation<dim, spacedim >::face_iterator &face,
// ============================================================
// CylindricalManifold
// ============================================================
-
template <int dim, int spacedim>
CylindricalManifold<dim, spacedim>::CylindricalManifold(const unsigned int axis,
- const double tolerance) :
+ const double tolerance)
+ :
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)
-{}
+{
+ // do not use static_assert to make dimension-independent programming
+ // easier.
+ Assert (spacedim==3,
+ ExcMessage("CylindricalManifold can only be used for spacedim==3!"));
+}
const Point<spacedim> &point_on_axis_,
const double tolerance) :
ChartManifold<dim,spacedim,3>(Tensor<1,3>({0,2.*numbers::PI,0})),
- normal_direction(Tensor<1, spacedim>()),
+ normal_direction(internal::compute_normal(direction_, true)),
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
+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<3> middle;
+ Point<spacedim> 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<3>()+direction*lambda;
+ return Point<spacedim>()+direction*lambda;
else // If not, using the ChartManifold should yield valid results.
- return ChartManifold<dim, 3, 3>::get_new_point(surrounding_points,
- weights);
+ return ChartManifold<dim, spacedim, 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
+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;
+ const Tensor<1,spacedim> normalized_point = space_point-point_on_axis;
const double lambda = normalized_point*direction;
- const Point<3> projection = point_on_axis + direction*lambda;
- const Tensor<1,3> p_diff = space_point - projection;
+ const Point<spacedim> projection = point_on_axis + direction*lambda;
+ const Tensor<1,spacedim> p_diff = space_point - projection;
// Then compute the angle between the projection direction and
// another vector orthogonal to the direction vector.
template <int dim, int spacedim>
Point<spacedim>
-CylindricalManifold<dim, spacedim>::push_forward(const Point<3> &/*chart_point*/) const
+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
// http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
// and unit vectors.
const double sine_r = std::sin(chart_point(1))*chart_point(0);
const double cosine_r = std::cos(chart_point(1))*chart_point(0);
- const Tensor<1,3> dxn = cross_product_3d(direction, normal_direction);
- const Tensor<1,3> intermediate = normal_direction*cosine_r+dxn*sine_r;
+ const Tensor<1, spacedim> dxn = cross_product_3d(direction, normal_direction);
+ const Tensor<1, spacedim> intermediate = normal_direction*cosine_r+dxn*sine_r;
// Finally, put everything together.
return point_on_axis+direction*chart_point(2)+intermediate;
-template<int dim, int spacedim>
+template <int dim, int spacedim>
DerivativeForm<1, 3, spacedim>
-CylindricalManifold<dim, spacedim>::push_forward_gradient(const Point<3> &/*chart_point*/) const
+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;
// Rotate the orthogonal direction by the given angle.
// and unit vectors.
const double sine = std::sin(chart_point(1));
const double cosine = std::cos(chart_point(1));
- const Tensor<1,3> dxn = cross_product_3d(direction, normal_direction);
- const Tensor<1,3> intermediate = normal_direction*cosine+dxn*sine;
+ const Tensor<1,spacedim> dxn = cross_product_3d(direction, normal_direction);
+ const Tensor<1,spacedim> intermediate = normal_direction*cosine+dxn*sine;
// derivative w.r.t the radius
derivatives[0][0] = intermediate[0];
projects / dealii.git / commitdiff