Add geometry file for moebiusstrip

dune/gfe/geometries/moebiusstrip.hh

+#ifndef DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH
+#define DUNE_GFE_MOEBIUSSTRIP_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 Möbius strip in 3D, where the base circle is the circle in the x-y-plane with radius_ around 0,0,0
+template <class T = double>
+class MoebiusStripProjection
+{
+  static const int dim = 3;
+  using Domain = FieldVector<T,dim>;
+  using Jacobian = FieldMatrix<T,dim,dim>;
+
+  T radius_;
+
+public:
+  MoebiusStripProjection (T radius)
+    : radius_(radius)
+  {}
+
+  /// \brief Project the coordinate to the MS using closest point projection
+  Domain operator() (const Domain& x) const
+  {
+    double nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    //center for the point x - it lies on the circle around (0,0,0) with radius radius_ in the x-y-plane
+    Domain center{x[0] * radius_/nrm, x[1] * radius_/nrm, 0}; 
+    double cosu = x[0]/nrm;
+    double sinu = x[1]/nrm;
+    double cosuhalf = std::sqrt((1+cosu)/2);
+    cosuhalf *= (sinu < 0) ? (-1) : 1 ; // u goes from 0 to 2*pi, so cosuhalf is negative if sinu < 0, we multiply that here
+    double sinuhalf = std::sqrt((1-cosu)/2); //u goes from 0 to 2*pi, so sinuhalf is always >= 0 
+
+    // now calculate line from center to the new point
+    // the direction is (cosuhalf*cosu,cosuhalf*sinu,sinuhalf), as can be seen from the formula to construct the MS, we normalize this vector
+    Domain centerToPointOnMS{cosuhalf*cosu,cosuhalf*sinu,sinuhalf};
+    centerToPointOnMS /= centerToPointOnMS.two_norm();
+
+    Domain centerToX = center - x;
+    // We need the length, let theta be the angle between the vector centerToX and center to centerToPointOnMS
+    // then cos(theta) = centerToX*centerToPointOnMS/(len(centerToX)*len(centerToPointOnMS))
+    // We want to project the point to MS, s.t. the angle between the MS and the projection is 90deg
+    // Then, length = cos(theta) * len(centerToX) = centerToX*centerToPointOnMS/len(centerToPointOnMS) = centerToX*centerToPointOnMS
+    double length = - centerToX * centerToPointOnMS;
+    centerToPointOnMS *= length;
+    return center + centerToPointOnMS;
+  }
+
+
+  /// \brief derivative of the projection to the MS
+  friend auto derivative (const MoebiusStripProjection& moebiusStrip)
+  {
+    DUNE_THROW(NotImplemented,"The derivative of the projection to the Möbius strip is not implemented yet!");
+    return [radius = moebiusStrip.radius_](const Domain& x)
+    {
+      Jacobian out;
+      return out;
+    };
+  }
+
+  /// \brief Normal Vector of the MS
+  Domain normal (const Domain& x) const
+  {
+    using std::sqrt;
+    Domain nVec = {x[0],x[1],0};
+    double nrm = std::sqrt(x[0]*x[0] + x[1]*x[1]);
+    double cosu = x[0]/nrm;
+    double sinu = x[1]/nrm;
+    double cosuhalf = std::sqrt((1+cosu)/2);
+    cosuhalf *= (sinu < 0) ? (-1) : 1 ; // u goes from 0 to 2*pi, so cosuhalf is negative if sinu < 0, we multiply that here
+    double sinuhalf = std::sqrt((1-cosu)/2);
+    nVec[2] = (-1)*cosuhalf*(cosu*x[0]+sinu*x[1])/sinuhalf;
+    nVec /= nVec.two_norm();
+    return nVec;
+  }
+
+  /// \brief The mean curvature of the MS
+  T mean_curvature (const Domain& /*x*/) const
+  {
+    DUNE_THROW(NotImplemented,"The mean curvature of the Möbius strip is not implemented yet!");
+    return 1;
+  }
+
+  /// \brief The area of the MS
+  T area () const
+  {
+    DUNE_THROW(NotImplemented,"The area of the Möbius strip is not implemented yet!");
+    return 1;
+  }
+};
+
+/// \brief construct a grid function representing the parametrization of a MS
+template <class Grid, class T>
+auto moebiusStripGridFunction (T radius)
+{
+  return analyticGridFunction<Grid>(MoebiusStripProjection<T>{radius});
+}
+
+} // end namespace Dune
+
+#endif // DUNE_GFE_MOEBIUSSTRIP_GRIDFUNCTION_HH