From b71f5bf40caf60da026e3372be27238e15129011 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 3 Sep 2013 16:30:01 +0000
Subject: [PATCH] Use ValueFactory for the harmonic energy test
Hence many input configurations are tested now
[[Imported from SVN: r9418]]
---
test/adolctest.cc | 83 +++++++++++++++++++++++++++++++++--------------
1 file changed, 59 insertions(+), 24 deletions(-)
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 264cd620..27331135 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -33,8 +33,24 @@
#include <dune/gfe/cosseratenergystiffness.hh>
#include <dune/gfe/localgeodesicfeadolcstiffness.hh>
+#include "valuefactory.hh"
+#include "multiindex.hh"
+
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, TargetSpace::distance(v[i],v[j]));
+ return d;
+}
+
+
+
template <int blocksize>
void compareMatrices(const Matrix<FieldMatrix<double,blocksize,blocksize> >& fdMatrix,
const Matrix<FieldMatrix<double,blocksize,blocksize> >& adolcMatrix)
@@ -58,54 +74,73 @@ void compareMatrices(const Matrix<FieldMatrix<double,blocksize,blocksize> >& fdM
int testHarmonicEnergy() {
+ typedef UnitVector<double,3> TargetSpace;
const int dim = 1;
+
+ std::cout << " --- Testing " << className<TargetSpace>() << ", domain dimension: " << dim << " ---" << std::endl;
+
+ // set up a test grid
typedef YaspGrid<dim> GridType;
FieldVector<double,dim> l(1);
std::array<int,dim> elements;
std::fill(elements.begin(), elements.end(), 1);
GridType grid(l,elements);
- typedef Q1LocalFiniteElement<double,double,dim> LocalFE;
- LocalFE localFiniteElement;
-
- size_t nDofs = localFiniteElement.localBasis().size();
-
- typedef UnitVector<double,3> TargetSpace;
- std::vector<TargetSpace> localSolution(nDofs);
-
- for (size_t i=0; i<nDofs; i++)
- localSolution[i] = FieldVector<double,3>(1.0);
-
- for (size_t i=0; i<localSolution.size(); i++)
- std::cout << localSolution[i] << std::endl;
-
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
typedef Q1NodalBasis<GridType::LeafGridView,double> Q1Basis;
Q1Basis q1Basis(grid.leafView());
+ typedef Q1LocalFiniteElement<double,double,dim> LocalFE;
+ LocalFE localFiniteElement;
+
+ // Assembler using finite differences
HarmonicEnergyLocalStiffness<GridType::LeafGridView,
Q1Basis::LocalFiniteElement,
TargetSpace> harmonicEnergyLocalStiffness;
- // Assembler using ADOL-C
+ // Assembler using ADOL-C
typedef TargetSpace::rebind<adouble>::other ATargetSpace;
HarmonicEnergyLocalStiffness<GridType::LeafGridView,
- Q1Basis::LocalFiniteElement,
- ATargetSpace> harmonicEnergyADOLCLocalStiffness;
+ Q1Basis::LocalFiniteElement,
+ ATargetSpace> harmonicEnergyADOLCLocalStiffness;
LocalGeodesicFEADOLCStiffness<GridType::LeafGridView,
Q1Basis::LocalFiniteElement,
TargetSpace> localGFEADOLCStiffness(&harmonicEnergyADOLCLocalStiffness);
+ size_t nDofs = localFiniteElement.localBasis().size();
- harmonicEnergyLocalStiffness.assembleHessian(*grid.leafbegin<0>(),localFiniteElement, localSolution);
+ std::vector<TargetSpace> localSolution(nDofs);
- localGFEADOLCStiffness.assembleHessian(*grid.leafbegin<0>(),localFiniteElement, localSolution);
+ std::vector<TargetSpace> testPoints;
+ ValueFactory<TargetSpace>::get(testPoints);
- compareMatrices(harmonicEnergyLocalStiffness.A_, localGFEADOLCStiffness.A_);
+ int nTestPoints = testPoints.size();
- return 0;
+ MultiIndex index(nDofs, nTestPoints);
+ int numIndices = index.cycle();
+
+ for (int i=0; i<numIndices; i++, ++index) {
+
+ for (size_t j=0; j<nDofs; j++)
+ localSolution[j] = testPoints[index[j]];
+
+ if (diameter(localSolution) > 0.5*M_PI)
+ continue;
+
+ for (size_t j=0; j<localSolution.size(); j++)
+ std::cout << localSolution[j] << std::endl;
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ harmonicEnergyLocalStiffness.assembleHessian(*grid.leafbegin<0>(),localFiniteElement, localSolution);
+
+ localGFEADOLCStiffness.assembleHessian(*grid.leafbegin<0>(),localFiniteElement, localSolution);
+
+ compareMatrices(harmonicEnergyLocalStiffness.A_, localGFEADOLCStiffness.A_);
+
+ }
+
+ return 0;
}
int testCosseratEnergy() {
--
GitLab