WIP: Add test computation for bending isometries using reducedcubichermitetrianglebasis

src/bending-isometries.cc

@@ -164,6 +164,55 @@ auto nodewiseIsometryTest(DeformationFunction& deformationGridViewFunction)
 }
 
 
+/**
+ * @brief Compute L2/L1-Isometry Error over the grid. 
+ * 
+ */
+template<class LocalDiscreteKirchhoffFunction, class GridView>
+auto isometryError(LocalDiscreteKirchhoffFunction& deformationFunction,
+                   GridView& gridView)
+{
+    // QuadratureRule for the integral of the L^2 error
+    QuadratureRuleKey quadKey(dim,6);
+
+    // The error to be compute
+    double isometryL2Error = 0;
+    double isometryL1Error = 0;
+
+    FieldMatrix<double,2,2> IdentityMatrix = ScaledIdentityMatrix<double,2>(1);
+    for (auto&& element : elements(gridView))
+    {
+      deformationFunction.bind(element);
+
+      // Get quadrature formula
+      quadKey.setGeometryType(element.type());
+      const auto& quadRule = QuadratureRuleCache<double, dim>::rule(quadKey);
+
+
+      for (auto quadPoint : quadRule)
+      {
+        auto quadPos = quadPoint.position();
+        const auto integrationElement = element.geometry().integrationElement(quadPos);
+        const auto weight = quadPoint.weight();
+
+        auto numDir = deformationFunction.evaluateDerivative(quadPos);
+        auto isoValue = numDir.transposed() * numDir;
+        // printmatrix(std::cout, isoValue , "isoValue: ", "--");
+
+        // integrate error
+        isometryL2Error += (isoValue - IdentityMatrix).frobenius_norm2()  * weight * integrationElement;
+        isometryL1Error += (isoValue - IdentityMatrix).frobenius_norm()  * weight * integrationElement;
+      }
+    }
+    std::cout << std::setprecision(6);
+    std::cout << "elements: " << gridView.size(0)
+              << "      "
+              << "Isometry L2-error: " << std::sqrt(isometryL2Error)
+              << "      "
+              << "Isometry L1-error: " << isometryL1Error << std::endl;
+}
+
+
 
 /** \brief The identity function in R^d, and its derivative
  *
@@ -854,6 +903,10 @@ int main(int argc, char *argv[])
     // nodewiseIsometryTest(localDKFunctionD, gridView);
     // std::cout << "compute END isometry Error done." << std::endl;
 
+    // compute Isometry-Error
+    isometryError(localDKFunctionD, gridView);
+
+
 
     /**
      * @brief TEST: evaluate global discrete displacement function