Implement the L^2(Gamma_F) norm for the state

src/one-body-sample.cc

@@ -47,11 +47,13 @@
 #include <dune/istl/scaledidmatrix.hh>
 
 #include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/dunepython.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
@@ -134,13 +136,6 @@ template <class FunctionMap> void initPython(FunctionMap &functions) {
       .toC<typename FunctionMap::Base>(functions);
 }
 
-template <class VectorType>
-double diff_two_norm(VectorType const &v1, VectorType const &v2) {
-  VectorType tmp = v1;
-  tmp -= v2;
-  return tmp.two_norm();
-}
-
 int main(int argc, char *argv[]) {
   try {
     Dune::Timer timer;
@@ -152,7 +147,9 @@ int main(int argc, char *argv[]) {
 
     typedef Dune::FieldVector<double, dims> SmallVector;
     typedef Dune::FieldMatrix<double, dims, dims> SmallMatrix;
+    typedef Dune::FieldMatrix<double, 1, 1> SmallSingletonMatrix;
     typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
+    typedef Dune::BCRSMatrix<SmallSingletonMatrix> SingletonMatrixType;
     typedef Dune::BlockVector<SmallVector> VectorType;
     typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
 
@@ -289,7 +286,23 @@ int main(int argc, char *argv[]) {
       localStiffness(E, nu);
       OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis)
           .assemble(localStiffness, stiffnessMatrix);
-    };
+    }
+
+    SingletonMatrixType frictionalBoundaryMassMatrix;
+    EnergyNorm<SingletonMatrixType, SingletonVectorType> const stateEnergyNorm(
+        frictionalBoundaryMassMatrix);
+    {
+      BoundaryPatch<GridView> const frictionalBoundary(leafView,
+                                                       frictionalNodes);
+
+      BoundaryMassAssembler<
+          GridType, BoundaryPatch<GridView>, P1Basis::LocalFiniteElement,
+          P1Basis::LocalFiniteElement, SmallSingletonMatrix> const
+      frictionalBoundaryMassAssembler(frictionalBoundary);
+
+      OperatorAssembler<P1Basis, P1Basis>(p1Basis, p1Basis).assemble(
+          frictionalBoundaryMassAssembler, frictionalBoundaryMassMatrix);
+    }
 
     auto const nodalIntegrals =
         assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
@@ -412,7 +425,7 @@ int main(int argc, char *argv[]) {
         {
           SingletonVectorType computed_state;
           stateUpdater->extractState(computed_state);
-          double const correction = diff_two_norm(computed_state, alpha);
+          double const correction = stateEnergyNorm.diff(computed_state, alpha);
           if (state_fpi <= 2 // Let the first two steps pass through unchanged
               || correction < minimalCorrectionReduction * lastCorrection) {
             alpha = computed_state;