Test geometry.normalGradients

dune/curvedgrid/test/convergence.hh

@@ -78,50 +78,22 @@ template <int order=2, class Grid, class Projection>
 typename Grid::ctype curvature_error (const Grid& grid, const Projection& projection, int quad_order = 8)
 {
   std::cout << "  curvature_error..." << std::endl;
-  using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
   using QuadProvider = Dune::QuadratureRules<typename Grid::ctype, 2>;
-  using FiniteElementType = Dune::PkLocalFiniteElement<typename Grid::ctype, typename Grid::ctype, 2, order+4>;
-  FiniteElementType fe;
-
-  using LBTraits = typename FiniteElementType::Traits::LocalBasisType::Traits;
-  std::vector<typename LBTraits::JacobianType> shapeGradients;
-  std::vector<GlobalCoordinate> gradients;
-  std::vector<GlobalCoordinate> normals;
 
   typename Grid::ctype dist = 0;
   for (const auto& element : elements(grid.leafGridView()))
   {
-    auto geometry = element.geometry();
-    const auto& localBasis = fe.localBasis();
-
-    { // get coefficients of normal vectors
-      const auto& localInterpolation = fe.localInterpolation();
-
-      auto n = [&](const auto& local) -> GlobalCoordinate
-      {
-        auto x = geometry.impl().normal(local);
-        return x;
-      };
-      localInterpolation.interpolate(n, normals);
-    }
-
-    shapeGradients.resize(localBasis.size());
-    gradients.resize(localBasis.size());
-
+    auto geometry = element.geometry().impl();
     const auto& quadRule = QuadProvider::rule(element.type(), quad_order);
-    for (const auto& qp : quadRule) {
-      auto jInvT = geometry.jacobianInverseTransposed(qp.position());
-      localBasis.evaluateJacobian(qp.position(), shapeGradients);
-
-      for (std::size_t i = 0; i < gradients.size(); ++i)
-        jInvT.mv(shapeGradients[i][0], gradients[i]);
-
-      typename Grid::ctype H = 0;
-      for (std::size_t i = 0; i < gradients.size(); ++i)
-        H += normals[i].dot(gradients[i]);
-
-      auto H_exact = projection.mean_curvature(geometry.global(qp.position()));
-      dist += std::abs(H/2 - H_exact) * qp.weight() * geometry.integrationElement(qp.position());
+    for (const auto& qp : quadRule)
+    {
+      auto H = geometry.normalGradient(qp.position());
+      typename Grid::ctype trH = 0;
+      for (int i = 0; i < H.N(); ++i)
+        trH += H[i][i];
+
+      auto trH_exact = projection.mean_curvature(geometry.global(qp.position()));
+      dist += std::abs(trH/2 - trH_exact) * qp.weight() * geometry.integrationElement(qp.position());
     }
   }