New test for the RigidBodyMotion class

In particular, we move the test for the computeDR method from
cosseratenergytest.cc to here.  It is a method that does not
have anything to do with Cosserat mechanics directly.

The test still seems overly complicated.  Do I really need
a LocalGeodesicFEFunction?  Maybe I can simplify it later.

test/CMakeLists.txt

@@ -9,6 +9,7 @@ set(TESTS
   localgfetestfunctiontest
   localprojectedfefunctiontest
   orthogonalmatrixtest
+  rigidbodymotiontest
   rodassemblertest
   rotationtest
   svdtest

test/cosseratenergytest.cc

 #include "config.h"
 
 #include <dune/grid/uggrid.hh>
-#include <dune/grid/onedgrid.hh>
 
 #include <dune/geometry/type.hh>
 #include <dune/geometry/quadraturerules.hh>
@@ -23,19 +22,6 @@ typedef RigidBodyMotion<double,3> TargetSpace;
 using namespace Dune;
 
 
-/** \brief Computes the diameter of a set */
-template <class TargetSpace>
-double diameter(const std::vector<TargetSpace>& v)
-{
-    double d = 0;
-    for (size_t i=0; i<v.size(); i++)
-        for (size_t j=0; j<v.size(); j++)
-            d = std::max(d, Rotation<double,3>::distance(v[i].q,v[j].q));
-    return d;
-}
-
-
-
 // ////////////////////////////////////////////////////////
 //   Make a test grid consisting of a single simplex
 // ////////////////////////////////////////////////////////
@@ -64,83 +50,6 @@ std::unique_ptr<GridType> makeSingleSimplexGrid()
 }
 
 
-template <int domainDim,class LocalFiniteElement>
-Tensor3<double,3,3,domainDim> evaluateDerivativeFD(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace>& f,
-                                           const Dune::FieldVector<double,domainDim>& local)
-{
-    Tensor3<double,3,3,domainDim> result(0);
-
-    for (int i=0; i<domainDim; i++) {
-        
-        Dune::FieldVector<double, domainDim> forward  = local;
-        Dune::FieldVector<double, domainDim> backward = local;
-        
-        forward[i]  += eps;
-        backward[i] -= eps;
-
-        TargetSpace forwardValue = f.evaluate(forward);
-        TargetSpace backwardValue = f.evaluate(backward);
-        
-        FieldMatrix<double,3,3> forwardMatrix, backwardMatrix;
-        forwardValue.q.matrix(forwardMatrix);
-        backwardValue.q.matrix(backwardMatrix);
-        
-        FieldMatrix<double,3,3> fdDer = (forwardMatrix - backwardMatrix) / (2*eps);
-        
-        for (int j=0; j<3; j++)
-            for (int k=0; k<3; k++)
-                result[j][k][i] = fdDer[j][k];
-        
-    }
-
-    return result;
-
-}
-
-template <int domainDim>
-void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
-{
-    // Make local fe function to be tested
-    PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
-    typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
-
-    LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement, TargetSpace> f(feCache.get(GeometryTypes::simplex(domainDim)), corners);
-
-    // A quadrature rule as a set of test points
-    int quadOrder = 3;
-    
-    const auto& quad = Dune::QuadratureRules<double, domainDim>::rule(GeometryTypes::simplex(domainDim), quadOrder);
-    
-    for (size_t pt=0; pt<quad.size(); pt++) {
-        
-        const Dune::FieldVector<double,domainDim>& quadPos = quad[pt].position();
-
-        // evaluate actual derivative
-        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
-
-        Tensor3<double,3,3,domainDim> DR;
-        typedef Dune::Functions::LagrangeBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
-        CosseratEnergyLocalStiffness<FEBasis,3>::computeDR(f.evaluate(quadPos),derivative, DR);
-
-        // evaluate fd approximation of derivative
-        Tensor3<double,3,3,domainDim> DR_fd = evaluateDerivativeFD(f,quadPos);
-
-        double maxDiff = 0;
-        for (int i=0; i<3; i++)
-            for (int j=0; j<3; j++)
-                for (int k=0; k<domainDim; k++)
-                    maxDiff = std::max(maxDiff, std::abs(DR[i][j][k] - DR_fd[i][j][k]));
-        
-        if ( maxDiff > 100*eps ) {
-            std::cout << className(corners[0]) << ": Analytical gradient does not match fd approximation." << std::endl;
-            std::cout << "Analytical:\n " << DR << std::endl;
-            std::cout << "FD        :\n " << DR_fd << std::endl;
-            assert(false);
-        }
-
-    }
-}
-
 //////////////////////////////////////////////////////////////////////////////////////
 //   Test invariance of the energy functional under rotations
 //////////////////////////////////////////////////////////////////////////////////////
@@ -256,42 +165,6 @@ int main(int argc, char** argv)
 
     const int domainDim = 2;
 
-    // ////////////////////////////////////////////////////////
-    //   Make a test grid consisting of a single simplex
-    // ////////////////////////////////////////////////////////
-
-    typedef UGGrid<domainDim> GridType;
-    const std::unique_ptr<GridType> grid = makeSingleSimplexGrid<GridType>();
-
-    ////////////////////////////////////////////////////////////////////////////
-    //  Create a local assembler object
-    ////////////////////////////////////////////////////////////////////////////
-
-    typedef Dune::Functions::LagrangeBasis<typename GridType::LeafGridView,1> Basis;
-    Basis basis(grid->leafGridView());
-
-    std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl;
-    
-    std::vector<Rotation<double,3> > testPoints;
-    
-    ValueFactory<Rotation<double,3> >::get(testPoints);
-    
-    int nTestPoints = testPoints.size();
-
-    // Set up elements of SO(3)
-    std::vector<TargetSpace> corners(domainDim+1);
-
-    ::MultiIndex index(domainDim+1, nTestPoints);
-    int numIndices = index.cycle();
-
-    for (int i=0; i<numIndices; i++, ++index)
-    {
-        for (int j=0; j<domainDim+1; j++)
-            corners[j].q = testPoints[index[j]];
-
-        testDerivativeOfRotationMatrix<domainDim>(corners);
-    }
-    
     //////////////////////////////////////////////////////////////////////////////////////
     //   Test invariance of the energy functional under rotations
     //////////////////////////////////////////////////////////////////////////////////////

test/rigidbodymotiontest.cc

+#include "config.h"
+
+#include <dune/geometry/type.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/uggrid.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/gfe/rigidbodymotion.hh>
+#include <dune/gfe/cosseratenergystiffness.hh>
+
+#include "multiindex.hh"
+#include "valuefactory.hh"
+
+const double eps = 1e-4;
+
+typedef RigidBodyMotion<double,3> TargetSpace;
+
+using namespace Dune;
+
+
+template <int domainDim,class LocalFiniteElement>
+Tensor3<double,3,3,domainDim> evaluateDerivativeFD(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace>& f,
+                                           const Dune::FieldVector<double,domainDim>& local)
+{
+  const double stepSize = std::sqrt(std::numeric_limits<double>::epsilon());
+    Tensor3<double,3,3,domainDim> result(0);
+
+    for (int i=0; i<domainDim; i++) {
+        
+        Dune::FieldVector<double, domainDim> forward  = local;
+        Dune::FieldVector<double, domainDim> backward = local;
+        
+        forward[i]  += stepSize;
+        backward[i] -= stepSize;
+
+        TargetSpace forwardValue = f.evaluate(forward);
+        TargetSpace backwardValue = f.evaluate(backward);
+        
+        FieldMatrix<double,3,3> forwardMatrix, backwardMatrix;
+        forwardValue.q.matrix(forwardMatrix);
+        backwardValue.q.matrix(backwardMatrix);
+        
+        FieldMatrix<double,3,3> fdDer = (forwardMatrix - backwardMatrix) / (2*stepSize);
+        
+        for (int j=0; j<3; j++)
+            for (int k=0; k<3; k++)
+                result[j][k][i] = fdDer[j][k];
+        
+    }
+
+    return result;
+
+}
+
+template <int domainDim>
+void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
+{
+    // Make local fe function to be tested
+    PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
+    typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
+
+    LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement, TargetSpace> f(feCache.get(GeometryTypes::simplex(domainDim)), corners);
+
+    // A quadrature rule as a set of test points
+    int quadOrder = 3;
+    
+    const auto& quad = Dune::QuadratureRules<double, domainDim>::rule(GeometryTypes::simplex(domainDim), quadOrder);
+    
+    for (size_t pt=0; pt<quad.size(); pt++) {
+        
+        const Dune::FieldVector<double,domainDim>& quadPos = quad[pt].position();
+
+        // evaluate actual derivative
+        Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
+
+        Tensor3<double,3,3,domainDim> DR;
+        typedef Dune::Functions::LagrangeBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
+        CosseratEnergyLocalStiffness<FEBasis,3>::computeDR(f.evaluate(quadPos),derivative, DR);
+
+        // evaluate fd approximation of derivative
+        Tensor3<double,3,3,domainDim> DR_fd = evaluateDerivativeFD(f,quadPos);
+
+        double maxDiff = 0;
+        for (int i=0; i<3; i++)
+            for (int j=0; j<3; j++)
+                for (int k=0; k<domainDim; k++)
+                    maxDiff = std::max(maxDiff, std::abs(DR[i][j][k] - DR_fd[i][j][k]));
+        
+        if ( maxDiff > 100*eps ) {
+            std::cout << className(corners[0]) << ": Analytical gradient does not match fd approximation." << std::endl;
+            std::cout << "Analytical:\n " << DR << std::endl;
+            std::cout << "FD        :\n " << DR_fd << std::endl;
+            assert(false);
+        }
+
+    }
+}
+
+int main(int argc, char** argv)
+{
+  MPIHelper::instance(argc, argv);
+
+  const int domainDim = 2;
+
+  ////////////////////////////////////////////////////////////////////////////
+  //  Create a local assembler object
+  ////////////////////////////////////////////////////////////////////////////
+  std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl;
+
+  std::vector<Rotation<double,3> > testPoints;
+
+  ValueFactory<Rotation<double,3> >::get(testPoints);
+
+  int nTestPoints = testPoints.size();
+
+  // Set up elements of SO(3)
+  std::vector<TargetSpace> corners(domainDim+1);
+
+  ::MultiIndex index(domainDim+1, nTestPoints);
+  int numIndices = index.cycle();
+
+  for (int i=0; i<numIndices; i++, ++index)
+  {
+      for (int j=0; j<domainDim+1; j++)
+          corners[j].q = testPoints[index[j]];
+
+      testDerivativeOfRotationMatrix<domainDim>(corners);
+  }
+}