const double tolerance = 1e-10);
/**
- * Compute the cartesian coordinates for a point given in cylindrical
+ * Compute the Cartesian coordinates for a point given in cylindrical
* coordinates.
*/
virtual Point<3>
pull_back(const Point<spacedim> &space_point) const override;
/**
- * Compute the cylindrical coordinates \f$(r, \phi, \lambda)\f$ for the given
- * point where \f$r\f$ denotes the distance from the axis,
- * \f$\phi\f$ the angle between the given point and the computed normal
- * direction and \f$\lambda\f$ the axial position.
+ * 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<spacedim>
push_forward(const Point<3> &chart_point) const override;
get_new_point(const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const override;
-private:
- /**
- * Compute a vector that is orthogonal to the given one.
- * Used for initializing the member variable normal_direction.
- */
- Point<spacedim> compute_normal(const Tensor<1,spacedim> &vector) const;
-
+protected:
/**
* A vector orthogonal to direcetion.
*/
*/
const Point<spacedim> point_on_axis;
+private:
/**
* Relative tolerance to measure zero distances.
*/
double tolerance;
+
+ // explicitly check for sensible template arguments, but not on windows
+ // because MSVC creates bogus warnings during normal compilation
+#ifndef DEAL_II_MSVC
+ static_assert (spacedim==3,
+ "CylindricalManifold can only be used for spacedim==3!");
+#endif
+
};
// CylindricalManifold
// ============================================================
+namespace
+{
+ // helper function to compute a vector orthogonal to a given one.
+ template <int spacedim>
+ Point<spacedim>
+ compute_normal(const Tensor<1,spacedim> &vector)
+ {
+ AssertThrow(vector.norm() != 0.,
+ ExcMessage("The direction parameter must not be zero!"));
+ Point<3> normal;
+ if (std::abs(vector[0]) >= std::abs(vector[1])
+ && std::abs(vector[0]) >= std::abs(vector[2]))
+ {
+ normal[1]=-1.;
+ normal[2]=-1.;
+ normal[0]=(vector[1]+vector[2])/vector[0];
+ }
+ else if (std::abs(vector[1]) >= std::abs(vector[0])
+ && std::abs(vector[1]) >= std::abs(vector[2]))
+ {
+ normal[0]=-1.;
+ normal[2]=-1.;
+ normal[1]=(vector[0]+vector[2])/vector[1];
+ }
+ else
+ {
+ normal[0]=-1.;
+ normal[1]=-1.;
+ normal[2]=(vector[0]+vector[1])/vector[2];
+ }
+ return normal;
+ }
+}
+
+
+
template <int dim, int spacedim>
CylindricalManifold<dim, spacedim>::CylindricalManifold(const unsigned int axis,
const double tolerance) :
- ChartManifold<dim,spacedim,3>(Tensor<1,3>({0,2.*numbers::PI,0})),
- normal_direction(compute_normal(Point<3>::unit_vector(axis))),
- direction (Point<3>::unit_vector(axis)),
- point_on_axis (Point<3>()),
- tolerance(tolerance)
-{
- static_assert(spacedim==3,
- "CylindricalManifold can only be used for spacedim==3!");
-}
+ CylindricalManifold<dim, spacedim>(Point<3>::unit_vector(axis),
+ Point<3>(),
+ tolerance)
+{}
direction (direction_/direction_.norm()),
point_on_axis (point_on_axis_),
tolerance(tolerance)
-{
- static_assert(spacedim==3,
- "CylindricalManifold can only be used for spacedim==3!");
-}
-
-
-
-template <int dim, int spacedim>
-Point<spacedim>
-CylindricalManifold<dim, spacedim>::compute_normal(const Tensor<1,spacedim> &vector) const
-{
- AssertThrow(vector.norm() != 0.,
- ExcMessage("The direction parameter must not be zero!"));
- Point<3> normal;
- if (std::abs(vector[0]) >= std::abs(vector[1])
- && std::abs(vector[0]) >= std::abs(vector[2]))
- {
- normal[1]=-1.;
- normal[2]=-1.;
- normal[0]=(vector[1]+vector[2])/vector[0];
- }
- else if (std::abs(vector[1]) >= std::abs(vector[0])
- && std::abs(vector[1]) >= std::abs(vector[2]))
- {
- normal[0]=-1.;
- normal[2]=-1.;
- normal[1]=(vector[0]+vector[2])/vector[1];
- }
- else
- {
- normal[0]=-1.;
- normal[1]=-1.;
- normal[2]=(vector[0]+vector[1])/vector[2];
- }
- return normal;
-}
+{}
get_new_point (const std::vector<Point<spacedim> > &surrounding_points,
const std::vector<double> &weights) const
{
- // First check if the the average in space lies on the axis.
+ // First check if the average in space lies on the axis.
Point<spacedim> middle;
double average_length = 0.;
for (unsigned int i=0; i<surrounding_points.size(); ++i)
#include <string>
#include <sstream>
+template <int dim>
+struct ManifoldWrapper
+{
+ Manifold<dim> *operator()(const Point<dim> &direction,
+ const Point<dim> ¢er ) const;
+};
+template <>
+Manifold<2> *
+ManifoldWrapper<2>::operator()(const Point<2> &/*direction*/,
+ const Point<2> ¢er ) const
+{
+ return new SphericalManifold<2>(center);
+}
+
+template <>
+Manifold<3> *
+ManifoldWrapper<3>::operator()(const Point<3> &direction,
+ const Point<3> ¢er ) const
+{
+ return new CylindricalManifold<3>(direction, center);
+}
template <int dim>
void test()
direction[dim-1] = 1.;
std::shared_ptr<Manifold<dim> > cylinder_manifold
- (dim == 2 ? static_cast<Manifold<dim>*>(new SphericalManifold<dim>(Point<dim>())) :
- static_cast<Manifold<dim>*>(new CylindricalManifold<dim>(direction, Point<dim>())));
+ (ManifoldWrapper<dim>()(direction, Point<dim>()));
// create cube and shift it to position 0.7
Triangulation<dim> tria;
}
}
+template <int dim>
+struct ManifoldWrapper
+{
+ Manifold<dim> *operator()(const Point<dim> &direction,
+ const Point<dim> ¢er ) const;
+};
+
+template <>
+Manifold<2> *
+ManifoldWrapper<2>::operator()(const Point<2> &/*direction*/,
+ const Point<2> ¢er ) const
+{
+ return new SphericalManifold<2>(center);
+}
+template <>
+Manifold<3> *
+ManifoldWrapper<3>::operator()(const Point<3> &direction,
+ const Point<3> ¢er ) const
+{
+ return new CylindricalManifold<3>(direction, center);
+}
template <int dim>
void test()
direction[dim-1] = 1.;
std::shared_ptr<Manifold<dim> > cylinder_manifold
- (dim == 2 ? static_cast<Manifold<dim>*>(new SphericalManifold<dim>(center)) :
- static_cast<Manifold<dim>*>(new CylindricalManifold<dim>(direction, center)));
+ (ManifoldWrapper<dim>()(direction, center));
Triangulation<dim> tria;
create_triangulation(tria);
tria.set_manifold(MANIFOLD_ID, *cylinder_manifold);