diff --git a/test/Makefile.am b/test/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..8e88968d8a341cb5c265d56e0042305040387be3
--- /dev/null
+++ b/test/Makefile.am
@@ -0,0 +1,17 @@
+# $Id: Makefile.am 1867 2007-12-31 15:45:00Z sander@PCPOOL.MI.FU-BERLIN.DE $
+
+# possible options
+LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS)
+AM_CPPFLAGS += $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) -Wall
+
+check_PROGRAMS = frameinvariancetest rotationtest fdcheck
+
+frameinvariancetest_SOURCES = frameinvariancetest.cc
+
+rotationtest_SOURCES = rotationtest.cc
+
+fdcheck_SOURCES = fdcheck.cc
+
+# don't follow the full GNU-standard
+# we need automake 1.5
+AUTOMAKE_OPTIONS = foreign 1.5
diff --git a/test/fdcheck.cc b/test/fdcheck.cc
new file mode 100644
index 0000000000000000000000000000000000000000..079caa45b5ef85188ca94d04951b4adcaf6b5cb9
--- /dev/null
+++ b/test/fdcheck.cc
@@ -0,0 +1,326 @@
+#include <config.h>
+
+#include <dune/common/bitfield.hh>
+#include <dune/common/configparser.hh>
+
+#include <dune/grid/onedgrid.hh>
+
+#include <dune/istl/io.hh>
+
+
+#include "../common/iterativesolver.hh"
+
+#include "src/configuration.hh"
+#include "src/quaternion.hh"
+#include "src/rodassembler.hh"
+
+#include "fdcheck.hh"
+
+// Number of degrees of freedom:
+// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
+const int blocksize = 6;
+
+using namespace Dune;
+
+void testDDExp()
+{
+ std::tr1::array<FieldVector<double,3>, 125> v;
+ int ct = 0;
+ double eps = 1e-4;
+
+ for (int i=-2; i<3; i++)
+ for (int j=-2; j<3; j++)
+ for (int k=-2; k<3; k++) {
+ v[ct][0] = i;
+ v[ct][1] = j;
+ v[ct][2] = k;
+ ct++;
+ }
+
+ for (size_t i=0; i<v.size(); i++) {
+
+ // Compute FD approximation of second derivative of exp
+ Dune::array<Dune::FieldMatrix<double,3,3>, 4> fdDDExp;
+
+ for (int j=0; j<3; j++) {
+
+ for (int k=0; k<3; k++) {
+
+ if (j==k) {
+
+ Quaternion<double> forwardQ = Quaternion<double>::exp(v[i][0] + (j==0)*eps,
+ v[i][1] + (j==1)*eps,
+ v[i][2] + (j==2)*eps);
+ Quaternion<double> centerQ = Quaternion<double>::exp(v[i][0],v[i][1],v[i][2]);
+ Quaternion<double> backwardQ = Quaternion<double>::exp(v[i][0] - (j==0)*eps,
+ v[i][1] - (j==1)*eps,
+ v[i][2] - (j==2)*eps);
+
+ for (int l=0; l<4; l++)
+ fdDDExp[l][j][j] = (forwardQ[l] - 2*centerQ[l] + backwardQ[l]) / (eps*eps);
+
+
+ } else {
+
+ FieldVector<double,3> ffV = v[i]; ffV[j] += eps; ffV[k] += eps;
+ FieldVector<double,3> fbV = v[i]; fbV[j] += eps; fbV[k] -= eps;
+ FieldVector<double,3> bfV = v[i]; bfV[j] -= eps; bfV[k] += eps;
+ FieldVector<double,3> bbV = v[i]; bbV[j] -= eps; bbV[k] -= eps;
+
+ Quaternion<double> forwardForwardQ = Quaternion<double>::exp(ffV);
+ Quaternion<double> forwardBackwardQ = Quaternion<double>::exp(fbV);
+ Quaternion<double> backwardForwardQ = Quaternion<double>::exp(bfV);
+ Quaternion<double> backwardBackwardQ = Quaternion<double>::exp(bbV);
+
+ for (int l=0; l<4; l++)
+ fdDDExp[l][j][k] = (forwardForwardQ[l] + backwardBackwardQ[l]
+ - forwardBackwardQ[l] - backwardForwardQ[l]) / (4*eps*eps);
+
+ }
+
+ }
+
+ }
+
+ // Compute analytical second derivative of exp
+ Dune::array<Dune::FieldMatrix<double,3,3>, 4> ddExp;
+ Quaternion<double>::DDexp(v[i], ddExp);
+
+ for (int m=0; m<4; m++)
+ for (int j=0; j<3; j++)
+ for (int k=0; k<3; k++)
+ if ( std::abs(fdDDExp[m][j][k] - ddExp[m][j][k]) > eps) {
+ std::cout << "Error at v = " << v[i]
+ << "[" << m << ", " << j << ", " << k << "] "
+ << " fd: " << fdDDExp[m][j][k]
+ << " analytical: " << ddExp[m][j][k] << std::endl;
+ }
+ }
+}
+
+void testDerivativeOfInterpolatedPosition()
+{
+ std::tr1::array<Quaternion<double>, 6> q;
+
+ FieldVector<double,3> xAxis(0); xAxis[0] = 1;
+ FieldVector<double,3> yAxis(0); yAxis[1] = 1;
+ FieldVector<double,3> zAxis(0); zAxis[2] = 1;
+
+ q[0] = Quaternion<double>(xAxis, 0);
+ q[1] = Quaternion<double>(xAxis, M_PI/2);
+ q[2] = Quaternion<double>(yAxis, 0);
+ q[3] = Quaternion<double>(yAxis, M_PI/2);
+ q[4] = Quaternion<double>(zAxis, 0);
+ q[5] = Quaternion<double>(zAxis, M_PI/2);
+
+ double eps = 1e-7;
+
+ for (int i=0; i<6; i++) {
+
+ for (int j=0; j<6; j++) {
+
+ for (int k=0; k<7; k++) {
+
+ double s = k/6.0;
+
+ std::tr1::array<Quaternion<double>,6> fdGrad;
+
+ // ///////////////////////////////////////////////////////////
+ // First: test the interpolated position
+ // ///////////////////////////////////////////////////////////
+ fdGrad[0] = Quaternion<double>::interpolate(q[i].mult(Quaternion<double>::exp(eps,0,0)), q[j], s);
+ fdGrad[0] -= Quaternion<double>::interpolate(q[i].mult(Quaternion<double>::exp(-eps,0,0)), q[j], s);
+ fdGrad[0] /= 2*eps;
+
+ fdGrad[1] = Quaternion<double>::interpolate(q[i].mult(Quaternion<double>::exp(0,eps,0)), q[j], s);
+ fdGrad[1] -= Quaternion<double>::interpolate(q[i].mult(Quaternion<double>::exp(0,-eps,0)), q[j], s);
+ fdGrad[1] /= 2*eps;
+
+ fdGrad[2] = Quaternion<double>::interpolate(q[i].mult(Quaternion<double>::exp(0,0,eps)), q[j], s);
+ fdGrad[2] -= Quaternion<double>::interpolate(q[i].mult(Quaternion<double>::exp(0,0,-eps)), q[j], s);
+ fdGrad[2] /= 2*eps;
+
+ fdGrad[3] = Quaternion<double>::interpolate(q[i], q[j].mult(Quaternion<double>::exp(eps,0,0)), s);
+ fdGrad[3] -= Quaternion<double>::interpolate(q[i], q[j].mult(Quaternion<double>::exp(-eps,0,0)), s);
+ fdGrad[3] /= 2*eps;
+
+ fdGrad[4] = Quaternion<double>::interpolate(q[i], q[j].mult(Quaternion<double>::exp(0,eps,0)), s);
+ fdGrad[4] -= Quaternion<double>::interpolate(q[i], q[j].mult(Quaternion<double>::exp(0,-eps,0)), s);
+ fdGrad[4] /= 2*eps;
+
+ fdGrad[5] = Quaternion<double>::interpolate(q[i], q[j].mult(Quaternion<double>::exp(0,0,eps)), s);
+ fdGrad[5] -= Quaternion<double>::interpolate(q[i], q[j].mult(Quaternion<double>::exp(0,0,-eps)), s);
+ fdGrad[5] /= 2*eps;
+
+ // Compute analytical gradient
+ std::tr1::array<Quaternion<double>,6> grad;
+ RodLocalStiffness<OneDGrid,double>::interpolationDerivative(q[i], q[j], s, grad);
+
+ for (int l=0; l<6; l++) {
+ Quaternion<double> diff = fdGrad[l];
+ diff -= grad[l];
+ if (diff.two_norm() > 1e-6) {
+ std::cout << "Error in position " << l << ": fd: " << fdGrad[l]
+ << " analytical: " << grad[l] << std::endl;
+ }
+
+ }
+
+ // ///////////////////////////////////////////////////////////
+ // Second: test the interpolated velocity vector
+ // ///////////////////////////////////////////////////////////
+
+ for (int l=1; l<7; l++) {
+
+ double intervalLength = l/(double(3));
+
+ fdGrad[0] = Quaternion<double>::interpolateDerivative(q[i].mult(Quaternion<double>::exp(eps,0,0)),
+ q[j], s, intervalLength);
+ fdGrad[0] -= Quaternion<double>::interpolateDerivative(q[i].mult(Quaternion<double>::exp(-eps,0,0)),
+ q[j], s, intervalLength);
+ fdGrad[0] /= 2*eps;
+
+ fdGrad[1] = Quaternion<double>::interpolateDerivative(q[i].mult(Quaternion<double>::exp(0,eps,0)),
+ q[j], s, intervalLength);
+ fdGrad[1] -= Quaternion<double>::interpolateDerivative(q[i].mult(Quaternion<double>::exp(0,-eps,0)),
+ q[j], s, intervalLength);
+ fdGrad[1] /= 2*eps;
+
+ fdGrad[2] = Quaternion<double>::interpolateDerivative(q[i].mult(Quaternion<double>::exp(0,0,eps)),
+ q[j], s, intervalLength);
+ fdGrad[2] -= Quaternion<double>::interpolateDerivative(q[i].mult(Quaternion<double>::exp(0,0,-eps)),
+ q[j], s, intervalLength);
+ fdGrad[2] /= 2*eps;
+
+ fdGrad[3] = Quaternion<double>::interpolateDerivative(q[i], q[j].mult(Quaternion<double>::exp(eps,0,0)), s, intervalLength);
+ fdGrad[3] -= Quaternion<double>::interpolateDerivative(q[i], q[j].mult(Quaternion<double>::exp(-eps,0,0)), s, intervalLength);
+ fdGrad[3] /= 2*eps;
+
+ fdGrad[4] = Quaternion<double>::interpolateDerivative(q[i], q[j].mult(Quaternion<double>::exp(0,eps,0)), s, intervalLength);
+ fdGrad[4] -= Quaternion<double>::interpolateDerivative(q[i], q[j].mult(Quaternion<double>::exp(0,-eps,0)), s, intervalLength);
+ fdGrad[4] /= 2*eps;
+
+ fdGrad[5] = Quaternion<double>::interpolateDerivative(q[i], q[j].mult(Quaternion<double>::exp(0,0,eps)), s, intervalLength);
+ fdGrad[5] -= Quaternion<double>::interpolateDerivative(q[i], q[j].mult(Quaternion<double>::exp(0,0,-eps)), s, intervalLength);
+ fdGrad[5] /= 2*eps;
+
+ // Compute analytical velocity vector gradient
+ RodLocalStiffness<OneDGrid,double>::interpolationVelocityDerivative(q[i], q[j], s, intervalLength, grad);
+
+ for (int m=0; m<6; m++) {
+ Quaternion<double> diff = fdGrad[m];
+ diff -= grad[m];
+ if (diff.two_norm() > 1e-6) {
+ std::cout << "Error in velocity " << m
+ << ": s = " << s << " of (" << intervalLength << ")"
+ << " fd: " << fdGrad[m] << " analytical: " << grad[m] << std::endl;
+ }
+
+ }
+
+ }
+
+ }
+
+ }
+
+ }
+
+}
+
+
+int main (int argc, char *argv[]) try
+{
+ // //////////////////////////////////////////////
+ // Test second derivative of exp
+ // //////////////////////////////////////////////
+ testDDExp();
+
+ // //////////////////////////////////////////////
+ // Test derivative of interpolated position
+ // //////////////////////////////////////////////
+ testDerivativeOfInterpolatedPosition();
+ exit(0);
+ typedef std::vector<Configuration> SolutionType;
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ typedef OneDGrid GridType;
+ GridType grid(1, 0, 1);
+
+ SolutionType x(grid.size(1));
+
+ // //////////////////////////
+ // Initial solution
+ // //////////////////////////
+ FieldVector<double,3> xAxis(0);
+ xAxis[0] = 1;
+ FieldVector<double,3> zAxis(0);
+ zAxis[2] = 1;
+
+ for (size_t i=0; i<x.size(); i++) {
+ x[i].r[0] = 0; // x
+ x[i].r[1] = 0; // y
+ x[i].r[2] = double(i)/(x.size()-1); // z
+ //x[i].r[2] = i+5;
+ x[i].q = Quaternion<double>::identity();
+ //x[i].q = Quaternion<double>(zAxis, M_PI/2 * double(i)/(x.size()-1));
+ }
+
+ //x.back().r[1] = 0.1;
+ //x.back().r[2] = 2;
+ //x.back().q = Quaternion<double>(zAxis, M_PI/4);
+
+ std::cout << "Left boundary orientation:" << std::endl;
+ std::cout << "director 0: " << x[0].q.director(0) << std::endl;
+ std::cout << "director 1: " << x[0].q.director(1) << std::endl;
+ std::cout << "director 2: " << x[0].q.director(2) << std::endl;
+ std::cout << std::endl;
+ std::cout << "Right boundary orientation:" << std::endl;
+ std::cout << "director 0: " << x[x.size()-1].q.director(0) << std::endl;
+ std::cout << "director 1: " << x[x.size()-1].q.director(1) << std::endl;
+ std::cout << "director 2: " << x[x.size()-1].q.director(2) << std::endl;
+// exit(0);
+
+ //x[0].r[2] = -1;
+
+ // ///////////////////////////////////////////
+ // Create a solver for the rod problem
+ // ///////////////////////////////////////////
+ RodAssembler<GridType> rodAssembler(grid);
+ //rodAssembler.setShapeAndMaterial(0.01, 0.0001, 0.0001, 2.5e5, 0.3);
+ //rodAssembler.setParameters(0,0,0,0,1,0);
+ rodAssembler.setParameters(0,0,100,0,0,0);
+
+ std::cout << "Energy: " << rodAssembler.computeEnergy(x) << std::endl;
+
+ double pos = (argc==2) ? atof(argv[1]) : 0.5;
+
+ FieldVector<double,1> shapeGrad[2];
+ shapeGrad[0] = -1;
+ shapeGrad[1] = 1;
+
+ FieldVector<double,1> shapeFunction[2];
+ shapeFunction[0] = 1-pos;
+ shapeFunction[1] = pos;
+
+ exit(0);
+ BlockVector<FieldVector<double,6> > rhs(x.size());
+ BCRSMatrix<FieldMatrix<double,6,6> > hessianMatrix;
+ MatrixIndexSet indices(grid.size(1), grid.size(1));
+ rodAssembler.getNeighborsPerVertex(indices);
+ indices.exportIdx(hessianMatrix);
+
+ rodAssembler.assembleGradient(x, rhs);
+ rodAssembler.assembleMatrix(x, hessianMatrix);
+
+ gradientFDCheck(x, rhs, rodAssembler);
+ hessianFDCheck(x, hessianMatrix, rodAssembler);
+
+ // //////////////////////////////
+ } catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+ }
diff --git a/test/fdcheck.hh b/test/fdcheck.hh
new file mode 100644
index 0000000000000000000000000000000000000000..09159862202a925217a4dbad7e623813f17b0e23
--- /dev/null
+++ b/test/fdcheck.hh
@@ -0,0 +1,479 @@
+#ifndef ASSEMBLER_FINITE_DIFFERENCE_CHECK
+#define ASSEMBLER_FINITE_DIFFERENCE_CHECK
+
+#define ABORT_ON_ERROR
+
+void infinitesimalVariation(Configuration& c, double eps, int i)
+{
+ if (i<3)
+ c.r[i] += eps;
+ else
+ c.q = c.q.mult(Quaternion<double>::exp((i==3)*eps,
+ (i==4)*eps,
+ (i==5)*eps));
+}
+
+template <class GridType>
+void strainFD(const std::vector<Configuration>& x,
+ double pos,
+ Dune::array<Dune::FieldMatrix<double,2,6>, 6>& fdStrainDerivatives,
+ const RodAssembler<GridType>& assembler)
+{
+ assert(x.size()==2);
+ double eps = 1e-8;
+
+ typename GridType::template Codim<0>::EntityPointer element = assembler.grid_->template leafbegin<0>();
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<Configuration> forwardSolution = x;
+ std::vector<Configuration> backwardSolution = x;
+
+ for (size_t i=0; i<x.size(); i++) {
+
+ Dune::FieldVector<double,6> fdGradient;
+
+ for (int j=0; j<6; j++) {
+
+ infinitesimalVariation(forwardSolution[i], eps, j);
+ infinitesimalVariation(backwardSolution[i], -eps, j);
+
+// fdGradient[j] = (assembler.computeEnergy(forwardSolution) - assembler.computeEnergy(backwardSolution))
+// / (2*eps);
+ Dune::FieldVector<double,6> strain;
+ strain = assembler.getStrain(forwardSolution, element, pos);
+ strain -= assembler.getStrain(backwardSolution, element, pos);
+ strain /= 2*eps;
+
+ for (int m=0; m<6; m++)
+ fdStrainDerivatives[m][i][j] = strain[m];
+
+ forwardSolution[i] = x[i];
+ backwardSolution[i] = x[i];
+ }
+
+ }
+
+}
+
+
+template <class GridType>
+void strain2ndOrderFD(const std::vector<Configuration>& x,
+ double pos,
+ Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer,
+ Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer,
+ const RodAssembler<GridType>& assembler)
+{
+ assert(x.size()==2);
+ double eps = 1e-3;
+
+ typename GridType::template Codim<0>::EntityPointer element = assembler.grid_->template leafbegin<0>();
+
+ for (int m=0; m<3; m++) {
+
+ translationDer[m].setSize(2,2);
+ translationDer[m] = 0;
+
+ rotationDer[m].setSize(2,2);
+ rotationDer[m] = 0;
+
+ }
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<Configuration> forwardSolution = x;
+ std::vector<Configuration> backwardSolution = x;
+
+ std::vector<Configuration> forwardForwardSolution = x;
+ std::vector<Configuration> forwardBackwardSolution = x;
+ std::vector<Configuration> backwardForwardSolution = x;
+ std::vector<Configuration> backwardBackwardSolution = x;
+
+ for (int i=0; i<2; i++) {
+
+ for (int j=0; j<3; j++) {
+
+ for (int k=0; k<2; k++) {
+
+ for (int l=0; l<3; l++) {
+
+ if (i==k && j==l) {
+
+ infinitesimalVariation(forwardSolution[i], eps, j+3);
+ infinitesimalVariation(backwardSolution[i], -eps, j+3);
+
+ // Second derivative
+// fdHessian[j][k] = (assembler.computeEnergy(forwardSolution)
+// - 2*assembler.computeEnergy(x)
+// + assembler.computeEnergy(backwardSolution)) / (eps*eps);
+
+ Dune::FieldVector<double,6> strain;
+ strain = assembler.getStrain(forwardSolution, element, pos);
+ strain += assembler.getStrain(backwardSolution, element, pos);
+ strain.axpy(-2, assembler.getStrain(x, element, pos));
+ strain /= eps*eps;
+
+ for (int m=0; m<3; m++)
+ rotationDer[m][i][k][j][l] = strain[3+m];
+
+ forwardSolution = x;
+ backwardSolution = x;
+
+ } else {
+
+ infinitesimalVariation(forwardForwardSolution[i], eps, j+3);
+ infinitesimalVariation(forwardForwardSolution[k], eps, l+3);
+ infinitesimalVariation(forwardBackwardSolution[i], eps, j+3);
+ infinitesimalVariation(forwardBackwardSolution[k], -eps, l+3);
+ infinitesimalVariation(backwardForwardSolution[i], -eps, j+3);
+ infinitesimalVariation(backwardForwardSolution[k], eps, l+3);
+ infinitesimalVariation(backwardBackwardSolution[i],-eps, j+3);
+ infinitesimalVariation(backwardBackwardSolution[k],-eps, l+3);
+
+
+// fdHessian[j][k] = (assembler.computeEnergy(forwardForwardSolution)
+// + assembler.computeEnergy(backwardBackwardSolution)
+// - assembler.computeEnergy(forwardBackwardSolution)
+// - assembler.computeEnergy(backwardForwardSolution)) / (4*eps*eps);
+
+ Dune::FieldVector<double,6> strain;
+ strain = assembler.getStrain(forwardForwardSolution, element, pos);
+ strain += assembler.getStrain(backwardBackwardSolution, element, pos);
+ strain -= assembler.getStrain(forwardBackwardSolution, element, pos);
+ strain -= assembler.getStrain(backwardForwardSolution, element, pos);
+ strain /= 4*eps*eps;
+
+
+ for (int m=0; m<3; m++)
+ rotationDer[m][i][k][j][l] = strain[3+m];
+
+ forwardForwardSolution = x;
+ forwardBackwardSolution = x;
+ backwardForwardSolution = x;
+ backwardBackwardSolution = x;
+
+
+ }
+
+ }
+
+ }
+ }
+
+ }
+
+}
+
+
+template <class GridType>
+void strain2ndOrderFD2(const std::vector<Configuration>& x,
+ double pos,
+ Dune::FieldVector<double,1> shapeGrad[2],
+ Dune::FieldVector<double,1> shapeFunction[2],
+ Dune::array<Dune::Matrix<Dune::FieldMatrix<double,6,6> >, 3>& translationDer,
+ Dune::array<Dune::Matrix<Dune::FieldMatrix<double,3,3> >, 3>& rotationDer,
+ const RodAssembler<GridType>& assembler)
+{
+ assert(x.size()==2);
+ double eps = 1e-3;
+
+ for (int m=0; m<3; m++) {
+
+ translationDer[m].setSize(2,2);
+ translationDer[m] = 0;
+
+ rotationDer[m].setSize(2,2);
+ rotationDer[m] = 0;
+
+ }
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<Configuration> forwardSolution = x;
+ std::vector<Configuration> backwardSolution = x;
+
+ for (int k=0; k<2; k++) {
+
+ for (int l=0; l<3; l++) {
+
+ infinitesimalVariation(forwardSolution[k], eps, l+3);
+ infinitesimalVariation(backwardSolution[k], -eps, l+3);
+
+ Dune::array<Dune::FieldMatrix<double,2,6>, 6> forwardStrainDer;
+ assembler.strainDerivative(forwardSolution, pos, shapeGrad, shapeFunction, forwardStrainDer);
+
+ Dune::array<Dune::FieldMatrix<double,2,6>, 6> backwardStrainDer;
+ assembler.strainDerivative(backwardSolution, pos, shapeGrad, shapeFunction, backwardStrainDer);
+
+ for (int i=0; i<2; i++) {
+ for (int j=0; j<3; j++) {
+ for (int m=0; m<3; m++) {
+ rotationDer[m][i][k][j][l] = (forwardStrainDer[m][i][j]-backwardStrainDer[m][i][j]) / (2*eps);
+ }
+ }
+ }
+
+ forwardSolution = x;
+ backwardSolution = x;
+
+ }
+
+ }
+
+}
+
+
+
+
+
+template <class GridType>
+void expHessianFD()
+{
+ using namespace Dune;
+ double eps = 1e-3;
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ FieldVector<double,3> forward;
+ FieldVector<double,3> backward;
+
+ FieldVector<double,3> forwardForward;
+ FieldVector<double,3> forwardBackward;
+ FieldVector<double,3> backwardForward;
+ FieldVector<double,3> backwardBackward;
+
+ for (int i=0; i<3; i++) {
+
+ for (int j=0; j<3; j++) {
+
+ Quaternion<double> hessian;
+
+ if (i==j) {
+
+ forward = backward = 0;
+ forward[i] += eps;
+ backward[i] -= eps;
+
+ // Second derivative
+ // fdHessian[j][k] = (assembler.computeEnergy(forward)
+ // - 2*assembler.computeEnergy(x)
+ // + assembler.computeEnergy(backward)) / (eps*eps);
+
+ hessian = Quaternion<double>::exp(forward);
+ hessian += Quaternion<double>::exp(backward);
+ hessian.axpy(-2, Quaternion<double>::exp(0,0,0));
+ hessian /= eps*eps;
+
+ } else {
+
+ forwardForward = forwardBackward = 0;
+ backwardForward = backwardBackward = 0;
+
+ forwardForward[i] += eps;
+ forwardForward[j] += eps;
+ forwardBackward[i] += eps;
+ forwardBackward[j] -= eps;
+ backwardForward[i] -= eps;
+ backwardForward[j] += eps;
+ backwardBackward[i] -= eps;
+ backwardBackward[j] -= eps;
+
+
+ hessian = Quaternion<double>::exp(forwardForward);
+ hessian += Quaternion<double>::exp(backwardBackward);
+ hessian -= Quaternion<double>::exp(forwardBackward);
+ hessian -= Quaternion<double>::exp(backwardForward);
+ hessian /= 4*eps*eps;
+
+ }
+
+ printf("i: %d, j: %d ", i, j);
+ std::cout << hessian << std::endl;
+
+ }
+
+ }
+}
+
+
+template <class GridType>
+void gradientFDCheck(const std::vector<Configuration>& x,
+ const Dune::BlockVector<Dune::FieldVector<double,6> >& gradient,
+ const RodAssembler<GridType>& assembler)
+{
+#ifndef ABORT_ON_ERROR
+ static int gradientError = 0;
+ double maxError = -1;
+#endif
+ double eps = 1e-6;
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+
+ std::vector<Configuration> forwardSolution = x;
+ std::vector<Configuration> backwardSolution = x;
+
+ for (size_t i=0; i<x.size(); i++) {
+
+ Dune::FieldVector<double,6> fdGradient;
+
+ for (int j=0; j<6; j++) {
+
+ infinitesimalVariation(forwardSolution[i], eps, j);
+ infinitesimalVariation(backwardSolution[i], -eps, j);
+
+ fdGradient[j] = (assembler.computeEnergy(forwardSolution) - assembler.computeEnergy(backwardSolution))
+ / (2*eps);
+
+ forwardSolution[i] = x[i];
+ backwardSolution[i] = x[i];
+ }
+
+ if ( (fdGradient-gradient[i]).two_norm() > 1e-6 ) {
+#ifdef ABORT_ON_ERROR
+ std::cout << "Differing gradients at " << i
+ << "! (" << (fdGradient-gradient[i]).two_norm() / gradient[i].two_norm()
+ << ") Analytical: " << gradient[i] << ", fd: " << fdGradient << std::endl;
+ //std::cout << "Current configuration " << std::endl;
+ for (size_t j=0; j<x.size(); j++)
+ std::cout << x[j] << std::endl;
+ //abort();
+#else
+ gradientError++;
+ maxError = std::max(maxError, (fdGradient-gradient[i]).two_norm());
+#endif
+ }
+
+ }
+
+#ifndef ABORT_ON_ERROR
+ std::cout << gradientError << " errors in the gradient corrected (max: " << maxError << ")!" << std::endl;
+#endif
+
+}
+
+
+template <class GridType>
+void hessianFDCheck(const std::vector<Configuration>& x,
+ const Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >& hessian,
+ const RodAssembler<GridType>& assembler)
+{
+#ifndef ABORT_ON_ERROR
+ static int hessianError = 0;
+#endif
+ double eps = 1e-2;
+
+ typedef typename Dune::BCRSMatrix<Dune::FieldMatrix<double,6,6> >::row_type::const_iterator ColumnIterator;
+
+ // ///////////////////////////////////////////////////////////
+ // Compute gradient by finite-difference approximation
+ // ///////////////////////////////////////////////////////////
+ std::vector<Configuration> forwardSolution = x;
+ std::vector<Configuration> backwardSolution = x;
+
+ std::vector<Configuration> forwardForwardSolution = x;
+ std::vector<Configuration> forwardBackwardSolution = x;
+ std::vector<Configuration> backwardForwardSolution = x;
+ std::vector<Configuration> backwardBackwardSolution = x;
+
+
+ // ///////////////////////////////////////////////////////////////
+ // Loop over all blocks of the outer matrix
+ // ///////////////////////////////////////////////////////////////
+ for (int i=0; i<hessian.N(); i++) {
+
+ ColumnIterator cIt = hessian[i].begin();
+ ColumnIterator cEndIt = hessian[i].end();
+
+ for (; cIt!=cEndIt; ++cIt) {
+
+ // ////////////////////////////////////////////////////////////////////////////
+ // Compute a finite-difference approximation of this hessian matrix block
+ // ////////////////////////////////////////////////////////////////////////////
+
+ Dune::FieldMatrix<double,6,6> fdHessian;
+
+ for (int j=0; j<6; j++) {
+
+ for (int k=0; k<6; k++) {
+
+ if (i==cIt.index() && j==k) {
+
+ infinitesimalVariation(forwardSolution[i], eps, j);
+ infinitesimalVariation(backwardSolution[i], -eps, j);
+
+ // Second derivative
+ fdHessian[j][k] = (assembler.computeEnergy(forwardSolution)
+ + assembler.computeEnergy(backwardSolution)
+ - 2*assembler.computeEnergy(x)
+ ) / (eps*eps);
+
+ forwardSolution[i] = x[i];
+ backwardSolution[i] = x[i];
+
+ } else {
+
+ infinitesimalVariation(forwardForwardSolution[i], eps, j);
+ infinitesimalVariation(forwardForwardSolution[cIt.index()], eps, k);
+ infinitesimalVariation(forwardBackwardSolution[i], eps, j);
+ infinitesimalVariation(forwardBackwardSolution[cIt.index()], -eps, k);
+ infinitesimalVariation(backwardForwardSolution[i], -eps, j);
+ infinitesimalVariation(backwardForwardSolution[cIt.index()], eps, k);
+ infinitesimalVariation(backwardBackwardSolution[i], -eps, j);
+ infinitesimalVariation(backwardBackwardSolution[cIt.index()],-eps, k);
+
+
+ fdHessian[j][k] = (assembler.computeEnergy(forwardForwardSolution)
+ + assembler.computeEnergy(backwardBackwardSolution)
+ - assembler.computeEnergy(forwardBackwardSolution)
+ - assembler.computeEnergy(backwardForwardSolution)) / (4*eps*eps);
+
+ forwardForwardSolution[i] = x[i];
+ forwardForwardSolution[cIt.index()] = x[cIt.index()];
+ forwardBackwardSolution[i] = x[i];
+ forwardBackwardSolution[cIt.index()] = x[cIt.index()];
+ backwardForwardSolution[i] = x[i];
+ backwardForwardSolution[cIt.index()] = x[cIt.index()];
+ backwardBackwardSolution[i] = x[i];
+ backwardBackwardSolution[cIt.index()] = x[cIt.index()];
+
+ }
+
+ }
+
+ }
+ //exit(0);
+ // /////////////////////////////////////////////////////////////////////////////
+ // Compare analytical and fd Hessians
+ // /////////////////////////////////////////////////////////////////////////////
+
+ Dune::FieldMatrix<double,6,6> diff = fdHessian;
+ diff -= *cIt;
+
+ if ( diff.frobenius_norm() > 1e-5 ) {
+#ifdef ABORT_ON_ERROR
+ std::cout << "Differing hessians at [(" << i << ", "<< cIt.index() << ")]!" << std::endl
+ << "Analytical: " << std::endl << *cIt << std::endl
+ << "fd: " << std::endl << fdHessian << std::endl;
+ abort();
+#else
+ hessianError++;
+#endif
+ }
+
+ }
+
+ }
+
+#ifndef ABORT_ON_ERROR
+ std::cout << hessianError << " errors in the Hessian corrected!" << std::endl;
+#endif
+
+}
+
+#endif
diff --git a/test/frameinvariancetest.cc b/test/frameinvariancetest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3f8b7f848b0d590809fc2635e2c2235951c7a06e
--- /dev/null
+++ b/test/frameinvariancetest.cc
@@ -0,0 +1,107 @@
+#include <config.h>
+
+//#define DUNE_EXPRESSIONTEMPLATES
+#include <dune/grid/onedgrid.hh>
+
+#include <dune/istl/io.hh>
+
+#include <dune/common/bitfield.hh>
+#include "src/quaternion.hh"
+
+#include "src/rodassembler.hh"
+
+#include "src/configuration.hh"
+#include "src/rodwriter.hh"
+
+// Number of degrees of freedom:
+// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
+const int blocksize = 6;
+
+using namespace Dune;
+using std::string;
+
+
+int main (int argc, char *argv[]) try
+{
+ // Some types that I need
+ typedef BCRSMatrix<FieldMatrix<double, blocksize, blocksize> > MatrixType;
+ typedef BlockVector<FieldVector<double, blocksize> > CorrectionType;
+ typedef std::vector<Configuration> SolutionType;
+
+ // Problem settings
+ const int numRodBaseElements = 1;
+
+ // ///////////////////////////////////////
+ // Create the grid
+ // ///////////////////////////////////////
+ typedef OneDGrid GridType;
+ GridType grid(numRodBaseElements, 0, 1);
+
+ SolutionType x(grid.size(1));
+
+ // //////////////////////////
+ // Initial solution
+ // //////////////////////////
+
+ for (size_t i=0; i<x.size(); i++) {
+ x[i].r[0] = 0; // x
+ x[i].r[1] = 0;
+ x[i].r[2] = double(i)/(x.size()-1); // z
+ x[i].q = Quaternion<double>::identity();
+ //x[i].q = Quaternion<double>(zAxis, (double(i)*M_PI)/(2*(x.size()-1)) );
+ }
+
+ FieldVector<double,3> zAxis(0); zAxis[2]=1;
+ x.back().q = Quaternion<double>(zAxis, M_PI/4);
+
+ // /////////////////////////////////////////////////////////////////////
+ // Create a second, rotated copy of the configuration
+ // /////////////////////////////////////////////////////////////////////
+
+ FieldVector<double,3> displacement;
+ displacement[0] = 0;
+ displacement[1] = 0;
+ displacement[2] = 0;
+
+ FieldVector<double,3> axis(0); axis[0]=1;
+ Quaternion<double> rotation(axis,M_PI/2);
+// std::cout << "Rotation:" << std::endl;
+// std::cout << "director 0: " << rotation.director(0) << std::endl;
+// std::cout << "director 1: " << rotation.director(1) << std::endl;
+// std::cout << "director 2: " << rotation.director(2) << std::endl;
+
+ SolutionType rotatedX = x;
+
+ for (size_t i=0; i<rotatedX.size(); i++) {
+
+ rotatedX[i].r = rotation.rotate(x[i].r);
+ rotatedX[i].r += displacement;
+
+ rotatedX[i].q = rotation.mult(x[i].q);
+
+// std::cout << "Vertex " << i << ": (" << rotatedX[i].r << ")" << std::endl;
+// // std::cout << "Right boundary orientation:" << std::endl;
+// std::cout << "director 0: " << rotatedX[i].q.director(0) << std::endl;
+// std::cout << "director 1: " << rotatedX[i].q.director(1) << std::endl;
+// std::cout << "director 2: " << rotatedX[i].q.director(2) << std::endl;
+ }
+
+ writeRod(x,"rod");
+ writeRod(rotatedX, "rotated");
+
+ RodLocalStiffness<GridType,double> assembler;
+
+ for (int i=1; i<2; i++) {
+
+ double p = double(i)/2;
+
+ assembler.getStrain(x,grid.lbegin<0>(0), p);
+ assembler.getStrain(rotatedX,grid.lbegin<0>(0), p);
+
+ }
+
+ } catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+ }
diff --git a/test/rotationtest.cc b/test/rotationtest.cc
new file mode 100644
index 0000000000000000000000000000000000000000..f781aa2eb84844adfcca91c2128801ea848cd075
--- /dev/null
+++ b/test/rotationtest.cc
@@ -0,0 +1,124 @@
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/fmatrix.hh>
+
+#include "src/quaternion.hh"
+#include "src/svd.hh"
+
+using namespace Dune;
+
+void testUnitQuaternion(Quaternion<double> q)
+{
+ // Make sure it really is a unit quaternion
+ q.normalize();
+
+ assert(std::abs(1-q.two_norm()) < 1e-12);
+
+ // Turn it into a matrix
+ FieldMatrix<double,3,3> matrix;
+ q.matrix(matrix);
+
+ // make sure it is an orthogonal matrix
+ if (std::abs(1-matrix.determinant()) > 1e-12 )
+ DUNE_THROW(Exception, "Expected determinant 1, but the computed value is " << matrix.determinant());
+
+ assert( std::abs( matrix[0]*matrix[1] ) < 1e-12 );
+ assert( std::abs( matrix[0]*matrix[2] ) < 1e-12 );
+ assert( std::abs( matrix[1]*matrix[2] ) < 1e-12 );
+
+ // Turn the matrix back into a quaternion, and check whether it is the same one
+ // Since the quaternions form a double covering of SO(3), we may either get q back
+ // or -q. We have to check both.
+ Quaternion<double> newQ;
+ newQ.set(matrix);
+
+ Quaternion<double> diff = newQ;
+ diff -= q;
+
+ Quaternion<double> sum = newQ;
+ sum += q;
+
+ if (diff.infinity_norm() > 1e-12 && sum.infinity_norm() > 1e-12)
+ DUNE_THROW(Exception, "Backtransformation failed for " << q << ". ");
+
+ // //////////////////////////////////////////////////////
+ // Check the director vectors
+ // //////////////////////////////////////////////////////
+
+ for (int i=0; i<3; i++)
+ for (int j=0; j<3; j++)
+ assert( std::abs(matrix[i][j] - q.director(j)[i]) < 1e-12 );
+
+ // //////////////////////////////////////////////////////
+ // Check multiplication with another unit quaternion
+ // //////////////////////////////////////////////////////
+
+ for (int i=-2; i<2; i++)
+ for (int j=-2; j<2; j++)
+ for (int k=-2; k<2; k++)
+ for (int l=-2; l<2; l++)
+ if (i!=0 || j!=0 || k!=0 || l!=0) {
+
+ Quaternion<double> q2(i,j,k,l);
+ q2.normalize();
+
+ // set up corresponding rotation matrix
+ FieldMatrix<double,3,3> q2Matrix;
+ q2.matrix(q2Matrix);
+
+ // q2 = q2 * q Quaternion multiplication
+ q2 = q2.mult(q);
+
+ // q2 = q2 * q Matrix multiplication
+ q2Matrix.rightmultiply(matrix);
+
+ FieldMatrix<double,3,3> productMatrix;
+ q2.matrix(productMatrix);
+
+ // Make sure we got identical results
+ productMatrix -= q2Matrix;
+ assert(productMatrix.infinity_norm() < 1e-10);
+
+ }
+
+ // ////////////////////////////////////////////////////////////////
+ // Check the operators 'B' that create an orthonormal basis of H
+ // ////////////////////////////////////////////////////////////////
+
+ Quaternion<double> Bq[4];
+ Bq[0] = q;
+ Bq[1] = q.B(0);
+ Bq[2] = q.B(1);
+ Bq[3] = q.B(2);
+
+ for (int i=0; i<4; i++) {
+
+ for (int j=0; j<4; j++) {
+
+ double prod = Bq[i]*Bq[j];
+ assert( std::abs( prod - (i==j) ) < 1e-6 );
+
+ }
+
+ }
+
+}
+
+int main (int argc, char *argv[]) try
+{
+ // Make a few random unit quaternions and pipe them into the test
+ for (int i=-5; i<5; i++)
+ for (int j=-5; j<5; j++)
+ for (int k=-5; k<5; k++)
+ for (int l=-5; l<5; l++)
+ if (i!=0 || j!=0 || k!=0 || l!=0)
+ testUnitQuaternion(Quaternion<double>(i,j,k,l));
+
+
+ } catch (Exception e) {
+
+ std::cout << e << std::endl;
+
+ }