From f193de2c909e763bd68e77f6cf4cb628528cc8a7 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 21 Dec 2017 12:36:24 +0100
Subject: [PATCH] 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 63e42ece72da4950160f61cdcc83eb38bea1d3c2 )
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...
---
test/cosseratenergytest.cc | 487 +++++--------------------------------
1 file changed, 62 insertions(+), 425 deletions(-)
diff --git a/test/cosseratenergytest.cc b/test/cosseratenergytest.cc
index 33e5c1fe..6881cb5d 100644
--- a/test/cosseratenergytest.cc
+++ b/test/cosseratenergytest.cc
@@ -6,7 +6,7 @@
#include <dune/geometry/type.hh>
#include <dune/geometry/quadraturerules.hh>
-#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#include <dune/functions/functionspacebases/pqknodalbasis.hh>
#include <dune/fufem/functions/constantfunction.hh>
@@ -16,10 +16,6 @@
#include "multiindex.hh"
#include "valuefactory.hh"
-#ifdef COSSERAT_ENERGY_FD_GRADIENT
-#warning Comparing two finite-difference approximations
-#endif
-
const double eps = 1e-4;
typedef RigidBodyMotion<double,3> TargetSpace;
@@ -45,7 +41,7 @@ double diameter(const std::vector<TargetSpace>& v)
// ////////////////////////////////////////////////////////
template <class GridType>
-std::auto_ptr<GridType> makeSingleSimplexGrid()
+std::unique_ptr<GridType> makeSingleSimplexGrid()
{
static const int domainDim = GridType::dimension;
GridFactory<GridType> factory;
@@ -62,17 +58,17 @@ std::auto_ptr<GridType> makeSingleSimplexGrid()
std::vector<unsigned int> v(domainDim+1);
for (int i=0; i<domainDim+1; 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>
-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)
{
- Tensor3<double,3,3,3> result(0);
+ Tensor3<double,3,3,domainDim> result(0);
for (int i=0; i<domainDim; i++) {
@@ -108,16 +104,12 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
PQkLocalFiniteElementCache<double,double,domainDim,1> feCache;
typedef typename PQkLocalFiniteElementCache<double,double,domainDim,1>::FiniteElementType LocalFiniteElement;
- GeometryType simplex;
- simplex.makeSimplex(domainDim);
-
- LocalGeodesicFEFunction<domainDim,double,LocalFiniteElement, TargetSpace> f(feCache.get(simplex), corners);
+ 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 Dune::QuadratureRule<double, domainDim>& quad
- = Dune::QuadratureRules<double, domainDim>::rule(GeometryType(GeometryType::simplex,domainDim), quadOrder);
+ const auto& quad = Dune::QuadratureRules<double, domainDim>::rule(GeometryTypes::simplex(domainDim), quadOrder);
for (size_t pt=0; pt<quad.size(); pt++) {
@@ -126,18 +118,17 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
// evaluate actual derivative
Dune::FieldMatrix<double, TargetSpace::EmbeddedTangentVector::dimension, domainDim> derivative = f.evaluateDerivative(quadPos);
- Tensor3<double,3,3,3> DR;
- CosseratEnergyLocalStiffness<typename UGGrid<domainDim>::LeafGridView,LocalFiniteElement,3>::computeDR(f.evaluate(quadPos),derivative, DR);
-
- //std::cout << "DR:\n" << DR << std::endl;
+ Tensor3<double,3,3,domainDim> DR;
+ typedef Dune::Functions::PQkNodalBasis<typename UGGrid<domainDim>::LeafGridView,1> FEBasis;
+ CosseratEnergyLocalStiffness<FEBasis,3>::computeDR(f.evaluate(quadPos),derivative, DR);
// 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;
for (int i=0; i<3; i++)
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]));
if ( maxDiff > 100*eps ) {
@@ -155,13 +146,8 @@ void testDerivativeOfRotationMatrix(const std::vector<TargetSpace>& corners)
//////////////////////////////////////////////////////////////////////////////////////
template <class GridType>
-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);
-
+void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficients)
+{
ParameterTree materialParameters;
materialParameters["thickness"] = "0.1";
materialParameters["mu"] = "3.8462e+05";
@@ -171,15 +157,19 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
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);
+ typedef Dune::Functions::PQkNodalBasis<typename GridType::LeafGridView,1> FEBasis;
+ FEBasis feBasis(grid->leafGridView());
+
+ CosseratEnergyLocalStiffness<FEBasis,3> assembler(materialParameters,
+ nullptr,
+ nullptr,
+ nullptr);
// compute reference energy
- double referenceEnergy = assembler.energy(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
+ auto localView = feBasis.localView();
+ localView.bind(*grid->leafGridView().template begin<0>());
+
+ double referenceEnergy = assembler.energy(localView,
coefficients);
// rotate the entire configuration
@@ -205,8 +195,7 @@ void testEnergy(const GridType* grid, const std::vector<TargetSpace>& coefficien
rotatedCoefficients[j].q = testRotations[i].mult(coefficients[j].q);
}
- double energy = assembler.energy(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
+ double energy = assembler.energy(localView,
rotatedCoefficients);
assert(std::fabs(energy-referenceEnergy) < 1e-4);
@@ -228,7 +217,7 @@ void testFrameInvariance()
// ////////////////////////////////////////////////////////
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
@@ -254,367 +243,12 @@ void testFrameInvariance()
}
-template <int domainDim, class LocalFiniteElement>
-void
-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()
+template <class Basis>
+void testEnergyGradient(Basis basis)
{
+ const int domainDim = Basis::GridView::dimension;
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;
materialParameters["thickness"] = "0.1";
materialParameters["mu"] = "3.8462e+05";
@@ -624,16 +258,19 @@ void testEnergyGradient()
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);
+ CosseratEnergyLocalStiffness<Basis,3> assembler(materialParameters,
+ nullptr,
+ nullptr,
+ nullptr);
// //////////////////////////////////////////////////////////////////////////////////////////
// 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;
ValueFactory<TargetSpace>::get(testPoints);
@@ -659,16 +296,13 @@ void testEnergyGradient()
}
// Compute the analytical gradient
- assembler.assembleGradient(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
+ assembler.assembleGradient(localView,
coefficients,
gradient);
// Compute the finite difference gradient
- assembler.LocalGeodesicFEStiffness<typename UGGrid<domainDim>::LeafGridView,
- LocalFiniteElement,
- RigidBodyMotion<double,3> >::assembleGradient(*grid->template leafbegin<0>(),
- feCache.get(grid->template leafbegin<0>()->type()),
+ assembler.LocalGeodesicFEStiffness<Basis,
+ RigidBodyMotion<double,3> >::assembleGradient(localView,
coefficients,
fdGradient);
@@ -700,9 +334,24 @@ void testEnergyGradient()
}
-int main(int argc, char** argv) try
+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::PQkNodalBasis<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;
@@ -717,13 +366,12 @@ int main(int argc, char** argv) try
MultiIndex index(domainDim+1, nTestPoints);
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++)
corners[j].q = testPoints[index[j]];
testDerivativeOfRotationMatrix<domainDim>(corners);
-
}
//////////////////////////////////////////////////////////////////////////////////////
@@ -732,17 +380,6 @@ int main(int argc, char** argv) try
testFrameInvariance<domainDim>();
- //////////////////////////////////////////////////////////////////////////////////////
- // Test the gradient of the energy functional
- //////////////////////////////////////////////////////////////////////////////////////
-
- testDerivativeOfGradientOfRotationMatrixWRTCoefficient<1>();
- testDerivativeOfGradientOfRotationMatrixWRTCoefficient<2>();
-
- testEnergyGradient<2>();
+ testEnergyGradient(basis);
-} catch (Exception e) {
-
- std::cout << e << std::endl;
-
- }
+}
--
GitLab