Skip to content
Snippets Groups Projects
Commit f193de2c authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Radical modernization of cosseratenergytest.cc 

In particular, this patch removes tests for the gradient
of the Cosserat energy and related things.  The
CosseratEnergyLocalStiffness class stopped computing
these things back in 2014 (in 63e42ece )
when it was decided that algorithmic differentiation was
the way to go.

With this patch, the test builds again, but it doesn't run yet.
I am afraid that fixing that may require a larger investment
of time...
parent 5e48a550
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
test/cosseratenergytest.cc
+ 62
− 425
View file @ f193de2c
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
#include <dune/geometry/type.hh> #include <dune/geometry/type.hh>
#include <dune/geometry/quadraturerules.hh> #include <dune/geometry/quadraturerules.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh> #include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/fufem/functions/constantfunction.hh> #include <dune/fufem/functions/constantfunction.hh>
...@@ -16,10 +16,6 @@ ...@@ -16,10 +16,6 @@
#include "multiindex.hh" #include "multiindex.hh"
#include "valuefactory.hh" #include "valuefactory.hh"
#ifdef COSSERAT_ENERGY_FD_GRADIENT
#warning Comparing two finite-difference approximations
#endif
const double eps = 1e-4; const double eps = 1e-4;
typedef RigidBodyMotion<double,3> TargetSpace; typedef RigidBodyMotion<double,3> TargetSpace;
...@@ -45,7 +41,7 @@ double diameter(const std::vector<TargetSpace>& v) ...@@ -45,7 +41,7 @@ double diameter(const std::vector<TargetSpace>& v)
// //////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////
template <class GridType> template <class GridType>
std::auto_ptr<GridType> makeSingleSimplexGrid() std::unique_ptr<GridType> makeSingleSimplexGrid()
{ {
static const int domainDim = GridType::dimension; static const int domainDim = GridType::dimension;
GridFactory<GridType> factory; GridFactory<GridType> factory;
...@@ -62,17 +58,17 @@ std::auto_ptr<GridType> makeSingleSimplexGrid() ...@@ -62,17 +58,17 @@ std::auto_ptr<GridType> makeSingleSimplexGrid()
std::vector<unsigned int> v(domainDim+1); std::vector<unsigned int> v(domainDim+1);
for (int i=0; i<domainDim+1; i++) for (int i=0; i<domainDim+1; i++)
v[i] = i; v[i] = i;
factory.insertElement(GeometryType(GeometryType::simplex,domainDim), v); factory.insertElement(GeometryTypes::simplex(domainDim), v);
return std::auto_ptr<GridType>(factory.createGrid()); return std::unique_ptr<GridType>(factory.createGrid());
} }
template <int domainDim,class LocalFiniteElement> template <int domainDim,class LocalFiniteElement>
Tensor3<double,3,3,3> evaluateDerivativeFD(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace>& f, Tensor3<double,3,3,domainDim> evaluateDerivativeFD(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace>& f,
const Dune::FieldVector<double,domainDim>& local) const Dune::FieldVector<double,domainDim>& local)
{ {
Tensor3<double,3,3,3> result(0); Tensor3<double,3,3,domainDim> result(0);
for (int i=0; i<domainDim; i++) { for (int i=0; i<domainDim; i++) {
...@@ -108,16 +104,12 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners) ...@@ -108,16 +104,12 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
PQkLocalFiniteElementCache<double,double,domainDim,1> feCache; PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement; typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
GeometryType simplex; LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement, TargetSpace> f(feCache.get(GeometryTypes::simplex(domainDim)), corners);
simplex.makeSimplex(domainDim);
LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement, TargetSpace> f(feCache.get(simplex), corners);
// A quadrature rule as a set of test points // A quadrature rule as a set of test points
int quadOrder = 3; int quadOrder = 3;
const Dune::QuadratureRule<double, domainDim>& quad const auto& quad = Dune::QuadratureRules<double, domainDim>::rule(GeometryTypes::simplex(domainDim), quadOrder);
= Dune::QuadratureRules<double, domainDim>::rule(GeometryType(GeometryType::simplex,domainDim), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) { for (size_t pt=0; pt<quad.size(); pt++) {
...@@ -126,18 +118,17 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners) ...@@ -126,18 +118,17 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
// evaluate actual derivative // evaluate actual derivative
Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos); Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
Tensor3<double,3,3,3> DR; Tensor3<double,3,3,domainDim> DR;
CosseratEnergyLocalStiffness<typename UGGrid<domainDim>::LeafGridView,LocalFiniteElement,3>::computeDR(f.evaluate(quadPos),derivative, DR); typedef Dune::Functions::PQkNodalBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
CosseratEnergyLocalStiffness<FEBasis,3>::computeDR(f.evaluate(quadPos),derivative, DR);
//std::cout << "DR:\n" << DR << std::endl;
// evaluate fd approximation of derivative // evaluate fd approximation of derivative
Tensor3<double,3,3,3> DR_fd = evaluateDerivativeFD(f,quadPos); Tensor3<double,3,3,domainDim> DR_fd = evaluateDerivativeFD(f,quadPos);
double maxDiff = 0; double maxDiff = 0;
for (int i=0; i<3; i++) for (int i=0; i<3; i++)
for (int j=0; j<3; j++) for (int j=0; j<3; j++)
for (int k=0; k<3; k++) for (int k=0; k<domainDim; k++)
maxDiff = std::max(maxDiff, std::abs(DR[i][j][k] - DR_fd[i][j][k])); maxDiff = std::max(maxDiff, std::abs(DR[i][j][k] - DR_fd[i][j][k]));
if ( maxDiff > 100*eps ) { if ( maxDiff > 100*eps ) {
...@@ -155,13 +146,8 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners) ...@@ -155,13 +146,8 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
template <class GridType> template <class GridType>
void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficients) { void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficients)
{
PQkLocalFiniteElementCache<double,double,GridType::dimension,1> feCache;
typedef typename PQkLocalFiniteElementCache<double,double,GridType::dimension,1>::FiniteElementType LocalFiniteElement;
//LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> f(feCache.get(element),corners);
ParameterTree materialParameters; ParameterTree materialParameters;
materialParameters["thickness"] = "0.1"; materialParameters["thickness"] = "0.1";
materialParameters["mu"] = "3.8462e+05"; materialParameters["mu"] = "3.8462e+05";
...@@ -171,15 +157,19 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien ...@@ -171,15 +157,19 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
materialParameters["q"] = "2.5"; materialParameters["q"] = "2.5";
materialParameters["kappa"] = "0.1"; materialParameters["kappa"] = "0.1";
ConstantFunction<Dune::FieldVector<double,GridType::dimension>, Dune::FieldVector<double,3> > zeroFunction(Dune::FieldVector<double,3>(0)); typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> FEBasis;
FEBasis feBasis(grid->leafGridView());
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3> assembler(materialParameters,
NULL, CosseratEnergyLocalStiffness<FEBasis,3> assembler(materialParameters,
&zeroFunction); nullptr,
nullptr,
nullptr);
// compute reference energy // compute reference energy
double referenceEnergy = assembler.energy(*grid->template leafbegin<0>(), auto localView = feBasis.localView();
feCache.get(grid->template leafbegin<0>()->type()), localView.bind(*grid->leafGridView().template begin<0>());
double referenceEnergy = assembler.energy(localView,
coefficients); coefficients);
// rotate the entire configuration // rotate the entire configuration
...@@ -205,8 +195,7 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien ...@@ -205,8 +195,7 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
rotatedCoefficients[j].q = testRotations[i].mult(coefficients[j].q); rotatedCoefficients[j].q = testRotations[i].mult(coefficients[j].q);
} }
double energy = assembler.energy(*grid->template leafbegin<0>(), double energy = assembler.energy(localView,
feCache.get(grid->template leafbegin<0>()->type()),
rotatedCoefficients); rotatedCoefficients);
assert(std::fabs(energy-referenceEnergy) < 1e-4); assert(std::fabs(energy-referenceEnergy) < 1e-4);
...@@ -228,7 +217,7 @@ void testFrameInvariance() ...@@ -228,7 +217,7 @@ void testFrameInvariance()
// //////////////////////////////////////////////////////// // ////////////////////////////////////////////////////////
typedef UGGrid<domainDim> GridType; typedef UGGrid<domainDim> GridType;
const std::auto_ptr<GridType> grid = makeSingleSimplexGrid<GridType>(); const std::unique_ptr<GridType> grid = makeSingleSimplexGrid<GridType>();
// ////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////
// Test whether the energy is invariant under isometries // Test whether the energy is invariant under isometries
...@@ -254,367 +243,12 @@ void testFrameInvariance() ...@@ -254,367 +243,12 @@ void testFrameInvariance()
} }
template <int domainDim, class LocalFiniteElement> template <class Basis>
void void testEnergyGradient(Basis basis)
evaluateFD_dDR_WRTCoefficient(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,RigidBodyMotion<double,3> >& f,
const Dune::FieldVector<double, domainDim>& local,
int coefficient,
Dune::array<Tensor3<double,3,3,4>, 3>& result)
{
typedef RigidBodyMotion<double,3> TargetSpace;
typedef double ctype;
double eps = 1e-3;
for (int j=0; j<4; j++) {
assert(f.type().isSimplex());
int ncoeff = domainDim+1;
std::vector<TargetSpace> cornersPlus(ncoeff);
std::vector<TargetSpace> cornersMinus(ncoeff);
for (int k=0; k<ncoeff; k++)
cornersMinus[k] = cornersPlus[k] = f.coefficient(k);
typename TargetSpace::CoordinateType aPlus = f.coefficient(coefficient).globalCoordinates();
typename TargetSpace::CoordinateType aMinus = f.coefficient(coefficient).globalCoordinates();
aPlus[j+3] += eps;
aMinus[j+3] -= eps;
cornersPlus[coefficient] = TargetSpace(aPlus);
cornersMinus[coefficient] = TargetSpace(aMinus);
LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> fPlus(f.localFiniteElement_,cornersPlus);
LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,TargetSpace> fMinus(f.localFiniteElement_,cornersMinus);
TargetSpace qPlus = fPlus.evaluate(local);
FieldMatrix<double,7,domainDim> qDerPlus = fPlus.evaluateDerivative(local);
TargetSpace qMinus = fMinus.evaluate(local);
FieldMatrix<double,7,domainDim> qDerMinus = fMinus.evaluateDerivative(local);
Tensor3<double,3,3,3> DRPlus,DRMinus;
typedef typename SelectType<domainDim==1,OneDGrid,UGGrid<domainDim> >::Type GridType;
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3>::computeDR(qPlus,qDerPlus,DRPlus);
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3>::computeDR(qMinus,qDerMinus,DRMinus);
for (int i=0; i<3; i++)
for (int k=0; k<3; k++)
for (int l=0; l<3; l++)
result[i][k][l][j] = (DRPlus[i][k][l] - DRMinus[i][k][l]) / (2*eps);
}
// Project onto the tangent space at M(q)
TargetSpace q = f.evaluate(local);
Dune::FieldMatrix<double,3,3> Mtmp;
q.q.matrix(Mtmp);
OrthogonalMatrix<double,3> M(Mtmp);
for (int k=0; k<domainDim; k++)
for (int v_i=0; v_i<4; v_i++) {
Dune::FieldMatrix<double,3,3> unprojected;
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
unprojected[i][j] = result[i][j][k][v_i];
Dune::FieldMatrix<double,3,3> projected = M.projectOntoTangentSpace(unprojected);
for (int i=0; i<3; i++)
for (int j=0; j<3; j++)
result[i][j][k][v_i] = projected[i][j];
}
}
template <int domainDim, class GridType, class LocalFiniteElement>
void testDerivativeOfGradientOfRotationMatrixWRTCoefficient(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,RigidBodyMotion<double,3> >& f,
const FieldVector<double,domainDim>& local,
unsigned int coeff)
{
RigidBodyMotion<double,3> value = f.evaluate(local);
FieldMatrix<double,7,domainDim> derOfValueWRTx = f.evaluateDerivative(local);
FieldMatrix<double,7,7> derOfValueWRTCoefficient;
f.evaluateDerivativeOfValueWRTCoefficient(local, coeff, derOfValueWRTCoefficient);
Tensor3<double,7,7,domainDim> derOfGradientWRTCoefficient;
f.evaluateDerivativeOfGradientWRTCoefficient(local, coeff, derOfGradientWRTCoefficient);
Tensor3<double,7,7,domainDim> FDderOfGradientWRTCoefficient;
f.evaluateFDDerivativeOfGradientWRTCoefficient(local, coeff, FDderOfGradientWRTCoefficient);
// Evaluate the analytical derivative
Dune::array<Tensor3<double,3,3,4>, 3> dDR_dv;
Tensor3<double,3,domainDim,4> dDR3_dv;
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3>::compute_dDR_dv(value,
derOfValueWRTx,
derOfValueWRTCoefficient,
derOfGradientWRTCoefficient,
dDR_dv,
dDR3_dv);
Dune::array<Tensor3<double,3,3,4>, 3> dDR_dv_FD;
evaluateFD_dDR_WRTCoefficient<domainDim,LocalFiniteElement>(f,local, coeff, dDR_dv_FD);
// Check whether the two are the same
double maxError = 0;
for (size_t j=0; j<3; j++)
for (size_t k=0; k<3; k++)
for (size_t l=0; l<3; l++)
for (size_t m=0; m<4; m++) {
double diff = std::fabs(dDR_dv[j][k][l][m] - dDR_dv_FD[j][k][l][m]);
maxError = std::max(maxError, diff);
}
if (maxError > 1e-3) {
std::cout << "Analytical and FD derivatives differ!"
<< " local: " << local
<< " coefficient: " << coeff << std::endl;
std::cout << "dDR_dv" << std::endl;
for (size_t i=0; i<dDR_dv.size(); i++)
std::cout << dDR_dv[i] << std::endl << std::endl;
std::cout << "finite differences: dDR_dv" << std::endl;
std::cout << dDR_dv_FD << std::endl;
abort();
}
}
template <int domainDim, class GridType, class LocalFiniteElement>
void testDerivativeOfBendingEnergy(const LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,RigidBodyMotion<double,3> >& f,
const FieldVector<double,domainDim>& local,
unsigned int coeff)
{
ParameterTree materialParameters;
materialParameters["thickness"] = "0.1";
materialParameters["mu"] = "3";
materialParameters["lambda"] = "2";
materialParameters["mu_c"] = "4";
materialParameters["L_c"] = "0.1";
materialParameters["q"] = "2.5";
materialParameters["kappa"] = "0.1";
ConstantFunction<Dune::FieldVector<double,GridType::dimension>, Dune::FieldVector<double,3> > zeroFunction(Dune::FieldVector<double,3>(0));
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3> assembler(materialParameters,
NULL,
&zeroFunction);
RigidBodyMotion<double,3> value = f.evaluate(local);
FieldMatrix<double,7,domainDim> derOfValueWRTx = f.evaluateDerivative(local);
FieldMatrix<double,7,7> derOfValueWRTCoefficient;
f.evaluateDerivativeOfValueWRTCoefficient(local, coeff, derOfValueWRTCoefficient);
Tensor3<double,7,7,domainDim> derOfGradientWRTCoefficient;
f.evaluateDerivativeOfGradientWRTCoefficient(local, coeff, derOfGradientWRTCoefficient);
Tensor3<double,7,7,domainDim> FDderOfGradientWRTCoefficient;
f.evaluateFDDerivativeOfGradientWRTCoefficient(local, coeff, FDderOfGradientWRTCoefficient);
// Evaluate the analytical derivative
Dune::array<Tensor3<double,3,3,4>, 3> dDR_dv;
Tensor3<double,3,domainDim,4> dDR3_dv;
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3>::compute_dDR_dv(value,
derOfValueWRTx,
derOfValueWRTCoefficient,
derOfGradientWRTCoefficient,
dDR_dv,
dDR3_dv);
//
Dune::FieldMatrix<double,3,3> R;
value.q.matrix(R);
Tensor3<double,3,3,3> DR;
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3>::computeDR(value, derOfValueWRTx, DR);
Tensor3<double,3,3,4> dR_dv;
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3>::computeDR_dv(value, derOfValueWRTCoefficient, dR_dv);
FieldVector<double,7> embeddedGradient;
assembler.bendingEnergyGradient(embeddedGradient,
R,
dR_dv,
DR,
dDR3_dv);
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
double eps = 1e-4;
FieldVector<double,4> embeddedFDGradient;
size_t nDofs = f.localFiniteElement_.localBasis().size();
std::vector<TargetSpace> forwardSolution(nDofs);
std::vector<TargetSpace> backwardSolution(nDofs);
for (size_t i=0; i<nDofs; i++)
forwardSolution[i] = backwardSolution[i] = f.coefficient(i);
for (int j=0; j<4; j++) {
FieldVector<double,7> forwardCorrection(0);
forwardCorrection[j+3] += eps;
FieldVector<double,7> backwardCorrection(0);
backwardCorrection[j+3] -= eps;
forwardSolution[coeff] = TargetSpace(f.coefficient(coeff).globalCoordinates() + forwardCorrection);
backwardSolution[coeff] = TargetSpace(f.coefficient(coeff).globalCoordinates() + backwardCorrection);
LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,RigidBodyMotion<double,3> > forwardF(f.localFiniteElement_,forwardSolution);
LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement,RigidBodyMotion<double,3> > backwardF(f.localFiniteElement_,backwardSolution);
RigidBodyMotion<double,3> forwardValue = forwardF.evaluate(local);
RigidBodyMotion<double,3> backwardValue = backwardF.evaluate(local);
FieldMatrix<double,3,3> forwardR, backwardR;
forwardValue.q.matrix(forwardR);
backwardValue.q.matrix(backwardR);
FieldMatrix<double,7,domainDim> forwardDerivative = forwardF.evaluateDerivative(local);
FieldMatrix<double,7,domainDim> backwardDerivative = backwardF.evaluateDerivative(local);
Tensor3<double,3,3,3> forwardDR;
Tensor3<double,3,3,3> backwardDR;
assembler.computeDR(forwardValue, forwardDerivative, forwardDR);
assembler.computeDR(backwardValue, backwardDerivative, backwardDR);
double forwardEnergy = assembler.bendingEnergy(forwardR, forwardDR);
double backwardEnergy = assembler.bendingEnergy(backwardR, backwardDR);
embeddedFDGradient[j] = (forwardEnergy - backwardEnergy) / (2*eps);
}
embeddedFDGradient = f.coefficient(coeff).q.projectOntoTangentSpace(embeddedFDGradient);
////////////////////////////////////////////////
// Check whether the two are the same
////////////////////////////////////////////////
double maxError = 0;
for (size_t i=0; i<4; i++) {
double diff = std::fabs(embeddedGradient[i+3] - embeddedFDGradient[i]);
maxError = std::max(maxError, diff);
}
if (maxError > 1e-3) {
std::cout << "Analytical and FD gradients differ!"
<< " local: " << local
<< " coefficient: " << coeff << std::endl;
std::cout << "embeddedGradient:" << std::endl;
std::cout << embeddedGradient[3] << " " << embeddedGradient[4] << " " << embeddedGradient[5] << " " << embeddedGradient[6] << std::endl << std::endl;
std::cout << "embeddedFDGradient:" << std::endl;
std::cout << embeddedFDGradient << std::endl;
abort();
}
}
template <int domainDim>
void testDerivativeOfGradientOfRotationMatrixWRTCoefficient()
{
std::cout << " --- test derivative of gradient of rotation matrix wrt coefficient ---" << std::endl;
// we just need the type
typedef typename SelectType<domainDim==1,OneDGrid,UGGrid<domainDim> >::Type GridType;
std::vector<TargetSpace> testPoints;
ValueFactory<TargetSpace>::get(testPoints);
// Set up elements of SE(3)
std::vector<TargetSpace> coefficients(domainDim+1);
MultiIndex index(domainDim+1, testPoints.size());
int numIndices = index.cycle();
for (int i=0; i<numIndices; i++, ++index) {
std::cout << "testing der of grad index: " << i << std::endl;
for (int j=0; j<domainDim+1; j++)
coefficients[j] = testPoints[index[j]];
if (diameter(coefficients) > M_PI-0.1) {
std::cout << "skipping, diameter = " << diameter(coefficients) << std::endl;
continue;
}
PQkLocalFiniteElementCache<double,double,GridType::dimension,1> feCache;
typedef typename PQkLocalFiniteElementCache<double,double,GridType::dimension,1>::FiniteElementType LocalFiniteElement;
GeometryType simplex;
simplex.makeSimplex(domainDim);
LocalGeodesicFEFunction<GridType::dimension,double,LocalFiniteElement,RigidBodyMotion<double,3> > f(feCache.get(simplex),coefficients);
// Loop over the coefficients -- we test derivation with respect to each of them
for (size_t j=0; j<coefficients.size(); j++) {
// A quadrature rule as a set of test points
int quadOrder = 3;
const Dune::QuadratureRule<double, domainDim>& quad
= Dune::QuadratureRules<double, domainDim>::rule(GeometryType(GeometryType::simplex,domainDim), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
FieldVector<double,domainDim> quadPos = quad[pt].position();
testDerivativeOfGradientOfRotationMatrixWRTCoefficient<domainDim,GridType,LocalFiniteElement>(f, quadPos, j);
testDerivativeOfBendingEnergy<domainDim,GridType,LocalFiniteElement>(f, quadPos, j);
}
}
}
}
template <int domainDim>
void testEnergyGradient()
{ {
const int domainDim = Basis::GridView::dimension;
std::cout << " --- Testing the gradient of the Cosserat energy functional, domain dimension: " << domainDim << " ---" << std::endl; std::cout << " --- Testing the gradient of the Cosserat energy functional, domain dimension: " << domainDim << " ---" << std::endl;
// ////////////////////////////////////////////////////////
// Make a test grid consisting of a single simplex
// ////////////////////////////////////////////////////////
typedef UGGrid<domainDim> GridType;
const std::auto_ptr<GridType> grid = makeSingleSimplexGrid<GridType>();
////////////////////////////////////////////////////////////////////////////
// Create a local assembler object
////////////////////////////////////////////////////////////////////////////
PQkLocalFiniteElementCache<double,double,GridType::dimension,1> feCache;
typedef typename PQkLocalFiniteElementCache<double,double,GridType::dimension,1>::FiniteElementType LocalFiniteElement;
ParameterTree materialParameters; ParameterTree materialParameters;
materialParameters["thickness"] = "0.1"; materialParameters["thickness"] = "0.1";
materialParameters["mu"] = "3.8462e+05"; materialParameters["mu"] = "3.8462e+05";
...@@ -624,16 +258,19 @@ void testEnergyGradient() ...@@ -624,16 +258,19 @@ void testEnergyGradient()
materialParameters["q"] = "2.5"; materialParameters["q"] = "2.5";
materialParameters["kappa"] = "0.1"; materialParameters["kappa"] = "0.1";
ConstantFunction<Dune::FieldVector<double,GridType::dimension>, Dune::FieldVector<double,3> > zeroFunction(Dune::FieldVector<double,3>(0)); CosseratEnergyLocalStiffness<Basis,3> assembler(materialParameters,
nullptr,
CosseratEnergyLocalStiffness<typename GridType::LeafGridView,LocalFiniteElement,3> assembler(materialParameters, nullptr,
NULL, nullptr);
&zeroFunction);
// ////////////////////////////////////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////////////////////////////////////
// Compare the gradient of the energy function with a finite difference approximation // Compare the gradient of the energy function with a finite difference approximation
// ////////////////////////////////////////////////////////////////////////////////////////// // //////////////////////////////////////////////////////////////////////////////////////////
auto element = *basis.gridView().template begin<0>();
auto localView = basis.localView();
localView.bind(element);
std::vector<TargetSpace> testPoints; std::vector<TargetSpace> testPoints;
ValueFactory<TargetSpace>::get(testPoints); ValueFactory<TargetSpace>::get(testPoints);
...@@ -659,16 +296,13 @@ void testEnergyGradient() ...@@ -659,16 +296,13 @@ void testEnergyGradient()
} }
// Compute the analytical gradient // Compute the analytical gradient
assembler.assembleGradient(*grid->template leafbegin<0>(), assembler.assembleGradient(localView,
feCache.get(grid->template leafbegin<0>()->type()),
coefficients, coefficients,
gradient); gradient);
// Compute the finite difference gradient // Compute the finite difference gradient
assembler.LocalGeodesicFEStiffness<typename UGGrid<domainDim>::LeafGridView, assembler.LocalGeodesicFEStiffness<Basis,
LocalFiniteElement, RigidBodyMotion<double,3> >::assembleGradient(localView,
RigidBodyMotion<double,3> >::assembleGradient(*grid->template leafbegin<0>(),
feCache.get(grid->template leafbegin<0>()->type()),
coefficients, coefficients,
fdGradient); fdGradient);
...@@ -700,9 +334,24 @@ void testEnergyGradient() ...@@ -700,9 +334,24 @@ void testEnergyGradient()
} }
int main(int argc, char** argv) try int main(int argc, char** argv)
{ {
const int domainDim = 2; 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::PQkNodalBasis<typename GridType::LeafGridView,1> Basis;
Basis basis(grid->leafGridView());
std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl; std::cout << " --- Testing derivative of rotation matrix, domain dimension: " << domainDim << " ---" << std::endl;
std::vector<Rotation<double,3> > testPoints; std::vector<Rotation<double,3> > testPoints;
...@@ -717,13 +366,12 @@ int main(int argc, char** argv) try ...@@ -717,13 +366,12 @@ int main(int argc, char** argv) try
MultiIndex index(domainDim+1, nTestPoints); MultiIndex index(domainDim+1, nTestPoints);
int numIndices = index.cycle(); int numIndices = index.cycle();
for (int i=0; i<numIndices; i++, ++index) { for (int i=0; i<numIndices; i++, ++index)
{
for (int j=0; j<domainDim+1; j++) for (int j=0; j<domainDim+1; j++)
corners[j].q = testPoints[index[j]]; corners[j].q = testPoints[index[j]];
testDerivativeOfRotationMatrix<domainDim>(corners); testDerivativeOfRotationMatrix<domainDim>(corners);
} }
////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////
...@@ -732,17 +380,6 @@ int main(int argc, char** argv) try ...@@ -732,17 +380,6 @@ int main(int argc, char** argv) try
testFrameInvariance<domainDim>(); testFrameInvariance<domainDim>();
////////////////////////////////////////////////////////////////////////////////////// testEnergyGradient(basis);
// Test the gradient of the energy functional
//////////////////////////////////////////////////////////////////////////////////////
testDerivativeOfGradientOfRotationMatrixWRTCoefficient<1>();
testDerivativeOfGradientOfRotationMatrixWRTCoefficient<2>();
testEnergyGradient<2>();
} catch (Exception e) { }
std::cout << e << std::endl;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment