Merge branch 'fix/projection-based' into 'master'

Fix the use of projection-based finite elements in cosserat-continuum

See merge request !78

dune/gfe/linearalgebra.hh

@@ -3,6 +3,7 @@
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/version.hh>
+#include <dune/istl/scaledidmatrix.hh>
 
 
 ///////////////////////////////////////////////////////////////////////////////////////////
@@ -13,6 +14,20 @@ namespace Dune {
 
   namespace GFE {
 
+  /** \brief Multiplication of a ScalecIdentityMatrix with another FieldMatrix */
+  template <class T, int N, int otherCols>
+  Dune::FieldMatrix<T,N,otherCols> operator* ( const Dune::ScaledIdentityMatrix<T, N>& diagonalMatrix,
+                        const Dune::FieldMatrix<T, N, otherCols>& matrix)
+  {
+      Dune::FieldMatrix<T,N,otherCols> result(0);
+
+      for (size_t i = 0; i < N; ++i)
+        for (size_t j = 0; j < otherCols; ++j)
+          result[i][j] = diagonalMatrix[i][i]*matrix[i][j];
+
+      return result;
+  }
+
   /** \brief Return the trace of a matrix */
   template <class T, int n>
   static T trace(const FieldMatrix<T,n,n>& A)

dune/gfe/localprojectedfefunction.hh

@@ -152,17 +152,7 @@ namespace Dune {
 
       auto derivativeOfProjection = TargetSpace::derivativeOfProjection(embeddedInterpolation);
 
-      typename LocalProjectedFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType result;
-
-      for (size_t i=0; i<result.N(); i++)
-        for (size_t j=0; j<result.M(); j++)
-        {
-          result[i][j] = 0;
-          for (size_t k=0; k<derivativeOfProjection.M(); k++)
-            result[i][j] += derivativeOfProjection[i][k]*derivative[k][j];
-        }
-
-      return result;
+      return derivativeOfProjection*derivative;
     }
 
     /** \brief Interpolate in an embedding Euclidean space, and project back onto the Riemannian manifold -- specialization for SO(3)

src/cosserat-continuum.cc

 #define MIXED_SPACE 0
+//#define PROJECTED_INTERPOLATION
 
 #include <config.h>
 
@@ -39,6 +40,8 @@
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/energynorm.hh>
 
+#include <dune/gfe/localgeodesicfefunction.hh>
+#include <dune/gfe/localprojectedfefunction.hh>
 #include <dune/gfe/localgeodesicfeadolcstiffness.hh>
 #include <dune/gfe/mixedlocalgfeadolcstiffness.hh>
 #include <dune/gfe/cosseratenergystiffness.hh>
@@ -318,7 +321,7 @@ int main (int argc, char *argv[]) try
       InitialBasis initialBasis(initialIterateGrid->leafGridView());
 
 #ifdef PROJECTED_INTERPOLATION
-      using LocalInterpolationRule  = LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
+      using LocalInterpolationRule  = GFE::LocalProjectedFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #else
       using LocalInterpolationRule  = LocalGeodesicFEFunction<dim, double, DeformationFEBasis::LocalView::Tree::FiniteElement, TargetSpace>;
 #endif

test/localprojectedfefunctiontest.cc

@@ -64,9 +64,9 @@ evaluateDerivativeFD(const LocalFunction& f, const Dune::FieldVector<ctype, dim>
 }
 
 
-template <int domainDim>
-void testDerivativeTangentiality(const RealTuple<double,1>& x,
-                                 const FieldMatrix<double,1,domainDim>& derivative)
+template <int domainDim, int dim>
+void testDerivativeTangentiality(const RealTuple<double,dim>& x,
+                                 const FieldMatrix<double,dim,domainDim>& derivative)
 {
     // By construction, derivatives of RealTuples are always tangent
 }
@@ -191,6 +191,7 @@ void testDerivative(const GFE::LocalProjectedFEFunction<domainDim,double,typenam
             std::cout << className<TargetSpace>() << ": Analytical gradient does not match fd approximation." << std::endl;
             std::cout << "Analytical: " << derivative << std::endl;
             std::cout << "FD        : " << fdDerivative << std::endl;
+            assert(false);
         }
 
         testDerivativeTangentiality(f.evaluate(quadPos), derivative);
@@ -262,6 +263,7 @@ int main()
 
     test<RealTuple<double,1>,2>(GeometryTypes::triangle);
     test<UnitVector<double,2>,2>(GeometryTypes::triangle);
+    test<RealTuple<double,3>,2>(GeometryTypes::triangle);
     test<UnitVector<double,3>,2>(GeometryTypes::triangle);
     test<Rotation<double,3>,2>(GeometryTypes::triangle);
     test<RigidBodyMotion<double,3>,2>(GeometryTypes::triangle);