Add Functor representing a cylinder in 3d

dune/curvedgrid/geometries/cylinder.hh

+#ifndef DUNE_CURVEDGRID_CYLINDER_GRIDFUNCTION_HH
+#define DUNE_CURVEDGRID_CYLINDER_GRIDFUNCTION_HH
+
+#include <cmath>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
+
+namespace Dune {
+
+/// \brief Functor representing a cylinder in 3D
+template <class T = double>
+class CylinderProjection
+{
+  static const int dim = 3;
+  using Domain = FieldVector<T,dim>;
+  using Jacobian = FieldMatrix<T,dim,dim>;
+
+  T radius_;
+
+public:
+  CylinderProjection (T radius)
+    : radius_(radius)
+  {}
+
+  /// \brief Project the coordinate to the cylinder at origin with radius
+  //         leaving the x[2] coordinate as it is
+  Domain operator() (const Domain& x) const
+  {
+    Domain point;
+    double radiusOverNorm = radius_/std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    point[0] = radiusOverNorm * x[0];
+    point[1] = radiusOverNorm * x[1];
+    point[2] = x[2];
+    return point;
+  }
+
+  /// \brief derivative of the projection to the cylinder
+  friend auto derivative (const CylinderProjection& cylinder)
+  {
+    return [r=cylinder.radius_](const Domain& x)
+    {
+      Jacobian out;
+      auto nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+      for (int i = 0; i < dim-1; ++i)
+        for (int j = 0; j < dim-1; ++j)
+          out[i][j] = r * ((i == j ? 1 : 0) - (x[i]/nrm) * (x[j]/nrm)) / nrm;
+      for (int i = 0; i < dim-1; ++i) {
+        out[2][i] = 0;
+        out[i][2] = 0;
+      }
+      out[2][2] = 1;
+      return out;
+    };
+  }
+
+  /// \brief normal vector, only points into x[0]-x[1]-direction, does not depend on x[2]
+  Domain normal (const Domain& x) const
+  {
+    Domain nVec = {0,0,0};
+    nVec[0] = x[0]/std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    nVec[1] = x[1]/std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    return nVec;
+  }
+
+  /// \brief mean curvature of the cylinder is 1/(2*radius)
+  T mean_curvature (const Domain& /*x*/) const
+  {
+    return T(1)/(2*radius_);
+  }
+
+  /// \brief area of the cylinder
+  T area (T height) const
+  {
+    return 2*M_PI*radius_*height;
+  }
+};
+
+/// \brief construct a grid function representing a cylinder parametrization
+template <class Grid, class T>
+auto cylinderGridFunction (T radius)
+{
+  return analyticGridFunction<Grid>(CylinderProjection<T>{radius});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_CURVEDGRID_CYLINDER_GRIDFUNCTION_HH