Fix gradientOf bug in DiscreteLocalFunction

amdis/gridfunctions/DiscreteLocalFunction.inc.hpp

@@ -162,7 +162,7 @@ private:
   // get the geometry of the bound element, that is constructed in the bind() method
   Geometry const& geometry() const
   {
-    assert(bound_);
+    assert(!!geometry_);
     return *geometry_;
   }
 
@@ -199,7 +199,7 @@ private:
   std::shared_ptr<Coefficients const> coefficients_;
   Type type_;
 
-  OptionalNoCopy<Geometry> geometry_;
+  std::optional<Geometry> geometry_;
   std::vector<Coefficient> localCoefficients_;
   bool bound_ = false;
 };

amdis/utility/LocalToGlobalAdapter.hpp

@@ -133,7 +133,7 @@ namespace AMDiS
                           std::vector<typename Traits::GradientRange>& out) const
     {
       auto const& localJacobian = cache_.localBasisJacobiansAt(x);
-      auto const& geoJacobian = geometry_.jacobianInverseTransposed(x);
+      auto&& geoJacobian = geometry_.jacobianInverseTransposed(x);
 
       out.resize(size_);
       for (std::size_t i = 0; i < size_; ++i)
@@ -159,7 +159,7 @@ namespace AMDiS
     {
       auto const& localJacobian = cache_.localBasisJacobiansAt(x);
       // NOTE: geoJacobian might be a Dune::DiagonalMatrix with limited interface!
-      auto const& geoJacobian = geometry_.jacobianInverseTransposed(x);
+      auto&& geoJacobian = geometry_.jacobianInverseTransposed(x);
 
       out.resize(size_);
       typename Traits::GradientRange grad;

test/GradientTest.cpp

@@ -91,6 +91,15 @@ void test(HostGrid& hostGrid)
   v_1 << evalAtQP([](auto const& x) { return std::sin(x[0])*std::sin(x[1]); });
   AMDIS_TEST((integrate(v_0 - partialDerivativeOf(u_0,0), gridView)) < 0.05);
   AMDIS_TEST((integrate(v_1 - partialDerivativeOf(u_1,1), gridView)) < 0.05);
+
+  {
+    using namespace Dune::Functions::BasisFactory;
+    DOFVector U0{gridView, lagrange<p>()};
+    DOFVector gradU0{gridView, power<HostGrid::dimensionworld>(lagrange<p>())};
+    valueOf(U0) << valueOf(uVector,0);
+    valueOf(gradU0) << gradientOf(U0);
+    [[maybe_unused]] auto integral1 = integrate(valueOf(gradU0,0), gridView);
+  }
 }
 
 int main(int argc, char** argv)