From 2f81ce9ae0f422077e438f9364b3bcf08dfcc7b4 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 3 Sep 2013 16:30:16 +0000
Subject: [PATCH] Test the Cosserat energy with all points from the test
function set
This still fails -- there must still be a bug in the Hessian computation
[[Imported from SVN: r9427]]
---
test/adolctest.cc | 87 +++++++++++++++++++++++++----------------------
1 file changed, 46 insertions(+), 41 deletions(-)
diff --git a/test/adolctest.cc b/test/adolctest.cc
index 83cb327d..85deaf7c 100644
--- a/test/adolctest.cc
+++ b/test/adolctest.cc
@@ -144,59 +144,62 @@ int testHarmonicEnergy() {
int testCosseratEnergy() {
- size_t nDofs = 4;
+ typedef RigidBodyMotion<double,3> TargetSpace;
+ const int gridDim = 2;
+
+ std::cout << " --- Testing " << className<TargetSpace>() << ", domain dimension: " << gridDim << " ---" << std::endl;
+
+ ParameterTree parameterSet;
+ ParameterTreeParser::readINITree("../cosserat-continuum.parset", parameterSet);
+
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+ std::cout << "Material parameters:" << std::endl;
+ materialParameters.report();
- const int dim = 2;
- typedef YaspGrid<dim> GridType;
- FieldVector<double,dim> l(1);
- std::array<int,dim> elements = {{1, 1}};
+ // set up a test grid
+ typedef YaspGrid<gridDim> GridType;
+ FieldVector<double,gridDim> l(1);
+ std::array<int,gridDim> elements;
+ std::fill(elements.begin(), elements.end(), 1);
GridType grid(l,elements);
- typedef Q1LocalFiniteElement<double,double,dim> LocalFE;
- LocalFE localFiniteElement;
+ typedef Q1NodalBasis<typename GridType::LeafGridView,double> Q1Basis;
+ Q1Basis q1Basis(grid.leafView());
- //typedef RealTuple<double,1> TargetSpace;
- typedef RigidBodyMotion<double,3> TargetSpace;
- std::vector<TargetSpace> localSolution(nDofs);
- FieldVector<double,7> identity(0);
- identity[6] = 1;
-
- for (auto vIt = grid.leafbegin<dim>(); vIt != grid.leafend<dim>(); ++vIt) {
- localSolution[grid.leafView().indexSet().index(*vIt)].r = 0;
- for (int i=0; i<dim; i++)
- localSolution[grid.leafView().indexSet().index(*vIt)].r[i] = 2*vIt->geometry().corner(0)[i];
- localSolution[grid.leafView().indexSet().index(*vIt)].q = Rotation<double,3>::identity();
- }
+ typedef Q1LocalFiniteElement<double,double,gridDim> LocalFE;
+ LocalFE localFiniteElement;
- for (size_t i=0; i<localSolution.size(); i++)
- std::cout << localSolution[i] << std::endl;
+ // Assembler using finite differences
+ CosseratEnergyLocalStiffness<GridType::LeafGridView,
+ Q1Basis::LocalFiniteElement,
+ 3> cosseratEnergyLocalStiffness(materialParameters,NULL,NULL);
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Assembler using ADOL-C
+ CosseratEnergyLocalStiffness<GridType::LeafGridView,
+ Q1Basis::LocalFiniteElement,
+ 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters, NULL, NULL);
+ LocalGeodesicFEADOLCStiffness<GridType::LeafGridView,
+ Q1Basis::LocalFiniteElement,
+ TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
- typedef Q1NodalBasis<GridType::LeafGridView,double> Q1Basis;
- Q1Basis q1Basis(grid.leafView());
+ size_t nDofs = localFiniteElement.localBasis().size();
- ParameterTree parameterSet;
- ParameterTreeParser::readINITree("../cosserat-continuum.parset", parameterSet);
+ std::vector<TargetSpace> localSolution(nDofs);
+ std::vector<TargetSpace> testPoints;
+ ValueFactory<TargetSpace>::get(testPoints);
- const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
- std::cout << "Material parameters:" << std::endl;
- materialParameters.report();
+ int nTestPoints = testPoints.size();
+ MultiIndex index(nDofs, nTestPoints);
+ int numIndices = index.cycle();
- CosseratEnergyLocalStiffness<GridType::LeafGridView,
- Q1Basis::LocalFiniteElement,
- 3> cosseratEnergyLocalStiffness(materialParameters,NULL,NULL);
+ for (int i=0; i<numIndices; i++, ++index) {
- // Assembler using ADOL-C
- CosseratEnergyLocalStiffness<GridType::LeafGridView,
- Q1Basis::LocalFiniteElement,
- 3,adouble> cosseratEnergyADOLCLocalStiffness(materialParameters, NULL, NULL);
- LocalGeodesicFEADOLCStiffness<GridType::LeafGridView,
- Q1Basis::LocalFiniteElement,
- TargetSpace> localGFEADOLCStiffness(&cosseratEnergyADOLCLocalStiffness);
+ for (size_t j=0; j<nDofs; j++)
+ localSolution[j] = testPoints[index[j]];
+ if (diameter(localSolution) > TargetSpace::convexityRadius)
+ continue;
cosseratEnergyLocalStiffness.assembleHessian(*grid.leafbegin<0>(),localFiniteElement, localSolution);
@@ -204,7 +207,9 @@ int testCosseratEnergy() {
compareMatrices(cosseratEnergyLocalStiffness.A_, localGFEADOLCStiffness.A_, localSolution);
- return 0;
+ }
+
+ return 0;
}
--
GitLab