Gradient of the normal field

dune/curvedgeometry/normalgradient.hh

+#ifndef DUNE_CURVEDGEOMETRY_NORMALGRADIENT_HH
+#define DUNE_CURVEDGEOMETRY_NORMALGRADIENT_HH
+
+#include <vector>
+
+#include <dune/common/fmatrix.hh>
+
+namespace Dune {
+
+template <class Geometry, class LFE>
+class NormalGradient
+{
+public:
+  /// type of the local coordinates in the geometry
+  using LocalCoordinate = typename Geometry::LocalCoordinates;
+
+  /// type of the global coordinates in the geometry
+  using GlobalCoordinate = typename Geometry::GlobalCoordinate;
+
+  /// type of the `J^{-T}` of the geometry
+  using JacobianInverseTransposed = typename Geometry::JacobianInverseTransposed;
+
+  /// type of the local finite-element
+  using FiniteElement = LFE;
+
+  /// traits types for the local basis
+  using BasisTraits = typename LFE::Traits::LocalBasisType::Traits;
+
+  /// type of the extended Weingarten map
+  using GlobalJacobian = FieldMatrix<typename Geometry::ctype, Geometry::coorddimension, Geometry::coorddimension>;
+
+public:
+  NormalGradient (const Geometry& geometry, const FiniteElement& localFE)
+    : geometry_(geometry)
+    , localFE_(localFE)
+  {
+    // create local discrete function of normal vectors by interpolation of the geometry normal
+    localFE_.localInterpolation().interpolate(
+      [&](const LocalCoordinate& local) { return geometry_.normalDirection(local); },
+      nCoefficients_);
+  }
+
+  GlobalJacobian operator() (const LocalCoordinate& local,
+                             const JacobianInverseTransposed& jiT) const
+  {
+    // Interpolated normal vector evaluated at local coordinate
+    localFE_.localBasis().evaluateFunction(local, nShapeValues_);
+    GlobalCoordinate nh(0);
+    for (std::size_t j = 0; j < nShapeValues_.size(); ++j)
+      nh.axpy(nShapeValues_[j], nCoefficients_[j]);
+    auto nh_nrm = nh.two_norm();
+    nh /= nh_nrm;
+
+    // P = I - n x n
+    GlobalJacobian Ph;
+    for (int r = 0; r < Geometry::coorddimension; ++r)
+      for (int s = 0; s < Geometry::coorddimension; ++s)
+        Ph[r][s] = (r == s ? 1 : 0) - nh[r]*nh[s];
+
+    // Compute the shape function gradients on the real element
+    localFE_.localBasis().evaluateJacobian(local, nShapeGradients_);
+    nGradients_.resize(nShapeGradients_.size());
+    for (std::size_t j = 0; j < nGradients_.size(); ++j)
+      jiT.mv(nShapeGradients_[j][0], nGradients_[j]);
+
+    // Normal gradient evaluated at local coordinate
+    GlobalJacobian H(0);
+    for (std::size_t j = 0; j < nGradients_.size(); ++j)
+      for (int r = 0; r < Geometry::coorddimension; ++r)
+        for (int s = 0; s < Geometry::coorddimension; ++s)
+          H[r][s] += nGradients_[j][s] * nCoefficients_[j][r];
+    H /= nh_nrm;
+    H.leftmultiply(Ph);
+    H.rightmultiply(Ph);
+
+    return H;
+  }
+
+  GlobalJacobian operator() (const LocalCoordinate& local) const
+  {
+    return (*this)(local, geometry_.jacobianInverseTransposed(local));
+  }
+
+private:
+  const Geometry& geometry_;
+  FiniteElement localFE_;
+
+  std::vector<GlobalCoordinate> nCoefficients_;
+  mutable std::vector<GlobalCoordinate> nGradients_;
+  mutable std::vector<typename BasisTraits::RangeType> nShapeValues_;
+  mutable std::vector<typename BasisTraits::JacobianType> nShapeGradients_;
+};
+
+
+template <class Geometry, class LFE>
+NormalGradient<Geometry,LFE> normalGradient (const Geometry& geometry, const LFE& localFE)
+{
+  return NormalGradient<Geometry,LFE>{geometry,localFE};
+}
+
+} // end namespace Dune
+
+#endif // DUNE_CURVEDGEOMETRY_NORMALGRADIENT_HH