From dde70eb727a4cf5af07ca6c65d8c9ce65ebc397d Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Mon, 21 Oct 2019 17:18:56 +0200
Subject: [PATCH] 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 | 1 +
test/cosseratenergytest.cc | 127 ----------------------------------
test/rigidbodymotiontest.cc | 131 ++++++++++++++++++++++++++++++++++++
3 files changed, 132 insertions(+), 127 deletions(-)
create mode 100644 test/rigidbodymotiontest.cc
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 68a1dafe..5c737a11 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -9,6 +9,7 @@ set(TESTS
localgfetestfunctiontest
localprojectedfefunctiontest
orthogonalmatrixtest
+ rigidbodymotiontest
rodassemblertest
rotationtest
svdtest
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 05dc2011..256fba6f 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -1,7 +1,6 @@
#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
//////////////////////////////////////////////////////////////////////////////////////
diff --git a/test/rigidbodymotiontest.cc b/test/rigidbodymotiontest.cc
new file mode 100644
index 00000000..58d9f5c1
--- /dev/null
+++ b/test/rigidbodymotiontest.cc
@@ -0,0 +1,131 @@
+#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);
+ }
+}
--
GitLab