diff --git a/test/adolcmaterialtest.cc b/test/adolcmaterialtest.cc
index c32f53190f4af36236ca090d5573e0bdf2f11817..a27f911e241ac21dc242c9fdb039898a85e1a2fc 100644
--- a/test/adolcmaterialtest.cc
+++ b/test/adolcmaterialtest.cc
@@ -14,25 +14,19 @@
#include <dune/istl/matrixindexset.hh>
#include <dune/grid/uggrid.hh>
-#include <dune/grid/geometrygrid/grid.hh>
-
#include <dune/fufem/makesphere.hh>
#include <dune/fufem/makering.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functions/deformationfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/adolclinearizationassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/adolchessianassembler.hh>
#define LAURSEN
#include <dune/elasticity/materials/adolcmaterial.hh>
#include <dune/elasticity/materials/neohookeanmaterial.hh>
-#include <dune/elasticity/materials/adolcneohookeanmaterial.hh>
-//#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
+#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
@@ -42,116 +36,85 @@ const int dim = 3;
using namespace Dune;
int main (int argc, char *argv[]) try
{
+ using FVector = FieldVector<double, dim>;
- double E = 1e5;
- if (argc==2)
- E = atof(argv[1]);
-
- typedef FieldVector<double,dim> FVector;
- typedef FieldMatrix<double,dim,dim> FMatrix;
- typedef BlockVector<FVector> VectorType;
- typedef BCRSMatrix<FMatrix> MatrixType;
-
- // ///////////////////////////////////////
// Create the grid
- // ///////////////////////////////////////
-
- typedef UGGrid<dim> GridType;
- auto grid = makeSphere<GridType>(FVector(0),0.2);
+ using Grid = UGGrid<dim>;
+ auto grid = makeSphere<Grid>(FVector(0), 0.2);
grid->globalRefine(3);
- typedef GridType::LeafGridView GridView;
- GridView gridView = grid->leafGridView();
- typedef P1NodalBasis<GridView> FEBasis;
- FEBasis feBasis(gridView);
+ using FEBasis = P1NodalBasis<typename Grid::LeafGridView>;
+ FEBasis feBasis(grid->leafGridView());
- // //////////////////////////
// Initial iterate
- // //////////////////////////
-
- VectorType x(feBasis.size());
-
double lower_bound = 0;
double upper_bound = 0.05;
- std::uniform_real_distribution<double> unif(lower_bound,upper_bound);
+ std::uniform_real_distribution<double> unif(lower_bound, upper_bound);
std::default_random_engine re;
- for (size_t i=0; i<x.size();i++)
- for (int j=0; j<dim;j++)
- x[i][j] = unif(re);
-
- auto xGridFunc = std::make_shared<BasisGridFunction<FEBasis,VectorType> >(feBasis,x);
-
- // ////////////////////////////////////////////////////////////
- // Create an assembler for the energy functional
- // ////////////////////////////////////////////////////////////
-
- double nu = 0.4;
-
- typedef FEBasis::LocalFiniteElement Lfe;
+ using BVector = BlockVector<FVector>;
+ BVector x(feBasis.size());
+ for (size_t i = 0; i < x.size(); ++i)
+ for (int j = 0; j < dim; ++j)
+ x[i][j] = unif(re);
+ auto xGridFunc = std::make_shared<BasisGridFunction<FEBasis,BVector> >(feBasis,x);
// make a material
- typedef NeoHookeanMaterial<FEBasis> MaterialType;
- MaterialType material(feBasis,E,nu);
+ using Material = NeoHookeanMaterial<FEBasis>;
+ Material material(feBasis, 1e5, 0.4);
// check against the Adol-C material
- LocalNeoHookeanEnergy<GridType,Lfe> localNeoHookeEnergy(E,nu);
- typedef AdolcMaterial<FEBasis> AdolcMaterialType;
- AdolcMaterialType adolcMaterial(feBasis,localNeoHookeEnergy,false);
+ AdolcMaterial<FEBasis> adolcMaterial(feBasis, material, false);
auto& linPaperAssembler = material.firstDerivative(xGridFunc);
auto& linAdolcAssembler = adolcMaterial.firstDerivative(xGridFunc);
FunctionalAssembler<FEBasis> functionalAssembler(feBasis);
- VectorType adolcGradient(feBasis.size());
- VectorType paperGrad(feBasis.size());
+ BVector adolcGradient(feBasis.size());
+ BVector paperGrad(feBasis.size());
adolcGradient = 0;
paperGrad = 0;
- Dune::Timer time;
- functionalAssembler.assemble(linAdolcAssembler,adolcGradient);
- std::cout<<"ADOL-C functional assembler needed "<<time.stop()<<std::endl;
+ Timer time;
+ functionalAssembler.assemble(linAdolcAssembler, adolcGradient);
+ std::cout << "ADOL-C functional assembler needed " << time.stop() << std::endl;
time.reset();
time.start();
- functionalAssembler.assemble(linPaperAssembler,paperGrad);
- std::cout<<"Paper&Pen functional assembler needed "<<time.stop()<<std::endl;
+ functionalAssembler.assemble(linPaperAssembler, paperGrad);
+ std::cout << "Paper&Pen functional assembler needed " << time.stop() << std::endl;
- for (size_t i=0; i<adolcGradient.size();i++) {
+ for (size_t i=0; i < adolcGradient.size(); ++i) {
- FVector diff = adolcGradient[i];
- diff -= paperGrad[i];
+ auto diff = adolcGradient[i] - paperGrad[i];
- if (diff.two_norm()>1e-8)
+ if (diff.two_norm()>1e-6)
DUNE_THROW(Dune::Exception,"Wrong local derivative, error is "<<diff.two_norm());
-
}
- // ////////////////////////
- // Test hessian assembler
- // ////////////////////////
-
+ // Test hessian assembler
auto& hessPaperAssembler = material.secondDerivative(xGridFunc);
auto& hessAdolcAssembler = adolcMaterial.secondDerivative(xGridFunc);
- OperatorAssembler<FEBasis,FEBasis> operatorAssembler(feBasis,feBasis);
+ OperatorAssembler<FEBasis, FEBasis> operatorAssembler(feBasis, feBasis);
- MatrixType adolcHessian, paperHessian;
+ using Matrix = BCRSMatrix<FieldMatrix<double, dim, dim> >;
+ Matrix adolcHessian, paperHessian;
time.reset();
time.start();
operatorAssembler.assemble(hessAdolcAssembler, adolcHessian);
- std::cout<<"ADOL-C operator assembler needed "<<time.stop()<<std::endl;
+ std::cout << "ADOL-C operator assembler needed " << time.stop() << std::endl;
time.reset();
time.start();
operatorAssembler.assemble(hessPaperAssembler, paperHessian);
- std::cout<<"Paper&Pen operator assembler needed "<<time.stop()<<std::endl;
+ std::cout << "Paper&Pen operator assembler needed " << time.stop() << std::endl;
- for (size_t i=0; i<adolcHessian.N();i++) {
+ for (size_t i=0; i < adolcHessian.N(); ++i) {
const auto& adolcRow = adolcHessian[i];
const auto& paperRow = paperHessian[i];
@@ -168,19 +131,18 @@ int main (int argc, char *argv[]) try
DUNE_THROW(Dune::Exception,"Not the same sparsity pattern!"<<adolcIt.index()
<<"!="<<paperIt.index());
- FMatrix diff = *adolcIt;
+ auto diff = *adolcIt;
diff -= *paperIt;
- if (diff.frobenius_norm()>1e-8)
+ if (diff.frobenius_norm() > 1e-6)
DUNE_THROW(Dune::Exception,"Wrong local hessian, error is "<<diff.frobenius_norm());
-
}
assert(paperIt==paperEndIt);
}
- // //////////////////////////////
-} catch (Exception e) {
+ return 0;
+} catch (Exception e) {
std::cout << e << std::endl;
-
+ return 1;
}
diff --git a/test/materialtest.cc b/test/materialtest.cc
index b2e7dd4e418f700320cd24f496eb62610392c561..4abaa888ddd98dae723091823a85cd5fbd08f7dc 100644
--- a/test/materialtest.cc
+++ b/test/materialtest.cc
@@ -1,246 +1,234 @@
-#include<config.h>
+#include <config.h>
#include <iomanip>
#include <random>
-#include <dune/common/fvector.hh>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
-#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
#include <dune/istl/matrixindexset.hh>
-#include <dune/grid/uggrid.hh>
#include <dune/grid/geometrygrid/grid.hh>
+#include <dune/grid/uggrid.hh>
-
-#include <dune/fufem/makesphere.hh>
-#include <dune/fufem/makering.hh>
+#include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/fufem/functions/deformationfunction.hh>
#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/assemblers/operatorassembler.hh>
-#include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/makering.hh>
+#include <dune/fufem/makesphere.hh>
#define LAURSEN
-#include <dune/elasticity/materials/neohookeanmaterial.hh>
-#include <dune/elasticity/materials/ogdenmaterial.hh>
#include <dune/elasticity/materials/geomexactstvenantkirchhoffmaterial.hh>
#include <dune/elasticity/materials/mooneyrivlinmaterial.hh>
+#include <dune/elasticity/materials/neohookeanmaterial.hh>
+#include <dune/elasticity/materials/adolcmaterial.hh>
+#include <dune/elasticity/materials/ogdenmaterial.hh>
-template<class Basis, class Material>
-bool checkHessian(const Basis& basis, Material& material,double eps)
-{
-
- const int dim = Basis::GridView::dimension;
-
- typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorType;
- typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,dim,dim> > MatrixType;
-
- VectorType conf(basis.getGridView().size(dim));
- conf = 0;
- double lower_bound = 0;
- double upper_bound = 0.005;
- std::uniform_real_distribution<double> unif(lower_bound,upper_bound);
- std::default_random_engine re;
+template <class Basis, class Material>
+bool checkHessian(const Basis& basis, Material& material, double eps) {
- for (size_t i=0; i<conf.size();i++)
- for (int j=0; j<dim;j++)
- conf[i][j] = unif(re);
+ constexpr int dim = Basis::GridView::dimension;
+ using Vector = Dune::BlockVector<Dune::FieldVector<double, dim>>;
+ double lower_bound{0}, upper_bound{0.005};
+ std::uniform_real_distribution<double> unif(lower_bound, upper_bound);
+ std::default_random_engine re;
- typedef BasisGridFunction<Basis,VectorType> GridFunc;
- auto confG = std::make_shared<GridFunc>(basis, conf);
+ Vector conf(basis.getGridView().size(dim));
+ for (size_t i = 0; i < conf.size(); ++i)
+ for (int j = 0; j < dim; ++j)
+ conf[i][j] = unif(re);
- // assemble Hessian using the material assemblers
- auto& locHessian = material.secondDerivative(confG);
+ using GridFunc = BasisGridFunction<Basis, Vector>;
+ auto confG = std::make_shared<GridFunc>(basis, conf);
- MatrixType exactHessian;
- OperatorAssembler<Basis,Basis> opAss(basis,basis);
- opAss.assemble(locHessian, exactHessian);
+ // assemble Hessian using the material assemblers
+ auto& locHessian = material.secondDerivative(confG);
- // compute finite difference approximation
+ OperatorAssembler<Basis, Basis> opAss(basis, basis);
- // dense matrix for the FD approximation
- Dune::MatrixIndexSet fdIndex(conf.size(),conf.size());
+ using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, dim, dim>>;
+ Matrix exactHessian;
+ opAss.assemble(locHessian, exactHessian);
- for (size_t i=0; i<conf.size();i++)
- for (size_t j=0; j<conf.size();j++)
- fdIndex.add(i,j);
+ // compute finite difference approximation
- // now the values
- MatrixType fdMat;
- fdIndex.exportIdx(fdMat);
- fdMat = 0;
+ // dense matrix for the FD approximation
+ Dune::MatrixIndexSet fdIndex(conf.size(), conf.size());
- FunctionalAssembler<Basis> funcAss(basis);
+ for (size_t i = 0; i < conf.size(); i++)
+ for (size_t j = 0; j < conf.size(); ++j)
+ fdIndex.add(i, j);
- for (size_t i=0; i<conf.size(); i++)
- for (int j=0; j<dim; j++) {
+ Matrix fdMat;
+ fdIndex.exportIdx(fdMat);
+ fdMat = 0;
+ FunctionalAssembler<Basis> funcAss(basis);
- VectorType forward = conf;
- VectorType backward = conf;
+ for (size_t i = 0; i < conf.size(); ++i)
+ for (int j = 0; j < dim; ++j) {
- forward[i][j] += eps;
- backward[i][j] -= eps;
+ auto forward = conf;
+ auto backward = conf;
- std::shared_ptr<GridFunc> fdConf(new GridFunc(basis,forward));
+ forward[i][j] += eps;
+ backward[i][j] -= eps;
- // assemble Hessian using the material assemblers
- typename Material::LocalLinearization& locLin = material.firstDerivative(fdConf);
+ auto fdConf = std::make_shared<GridFunc>(basis, forward);
- VectorType forw;
- funcAss.assemble(locLin,forw);
+ Vector forw;
+ funcAss.assemble(material.firstDerivative(fdConf), forw);
- std::shared_ptr<GridFunc> fdConfB(new GridFunc(basis,backward));
- auto& locLin2 = material.firstDerivative(fdConfB);
+ fdConf = std::make_shared<GridFunc>(basis, backward);
- VectorType backw;
- funcAss.assemble(locLin2,backw);
+ Vector backw;
+ funcAss.assemble(material.firstDerivative(fdConf), backw);
- for(size_t k=0; k<forw.size(); k++)
- for (int l=0; l<dim; l++)
- fdMat[k][i][l][j] = (forw[k][l]-backw[k][l])/(2*eps);
- }
+ for (size_t k = 0; k < forw.size(); ++k)
+ for (int l = 0; l < dim; ++l)
+ fdMat[k][i][l][j] = (forw[k][l] - backw[k][l]) / (2 * eps);
+ }
- MatrixType diff = fdMat;
- diff -= exactHessian;
- std::cout<<"Maximal difference "<<std::setprecision(12)<<diff.infinity_norm()<<"\n";
+ auto diff = fdMat;
+ diff -= exactHessian;
+ std::cout << "Maximal difference " << std::setprecision(12) << diff.infinity_norm() << "\n";
- return (diff.infinity_norm()<1e-3);
+ return (diff.infinity_norm() < 1e-3);
}
-template<class Basis, class Material>
-bool checkLinearization(const Basis& basis, Material& material,double eps)
-{
- const int dim = Basis::GridView::dimension;
-
- typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorType;
+template <class Basis, class Material>
+bool checkLinearization(const Basis& basis, Material& material, double eps) {
- VectorType conf(basis.getGridView().size(dim));
- conf = 0;
+ constexpr int dim = Basis::GridView::dimension;
+ using Vector = Dune::BlockVector<Dune::FieldVector<double, dim> >;
- double lower_bound = 0;
- double upper_bound = 0.05;
- std::uniform_real_distribution<double> unif(lower_bound,upper_bound);
- std::default_random_engine re;
+ double lower_bound{0}, upper_bound{0.05};
+ std::uniform_real_distribution<double> unif(lower_bound, upper_bound);
+ std::default_random_engine re;
- for (size_t i=0; i<conf.size();i++)
- for (int j=0; j<dim;j++)
- conf[i][j] = unif(re);
+ Vector conf(basis.getGridView().size(dim));
+ for (size_t i = 0; i < conf.size(); ++i)
+ for (int j = 0; j < dim; ++j)
+ conf[i][j] = unif(re);
- typedef BasisGridFunction<Basis,VectorType> GridFunc;
- auto confG = std::make_shared<GridFunc>(basis, conf);
+ using GridFunc = BasisGridFunction<Basis, Vector>;
+ auto confG = std::make_shared<GridFunc>(basis, conf);
- // assemble Hessian using the material assemblers
- auto& locLinear = material.firstDerivative(confG);
+ // assemble Hessian using the material assemblers
+ auto& locLinear = material.firstDerivative(confG);
- VectorType exactLinear;
- FunctionalAssembler<Basis> funcAss(basis);
- funcAss.assemble(locLinear, exactLinear);
+ FunctionalAssembler<Basis> funcAss(basis);
+ Vector exactLinear;
+ funcAss.assemble(locLinear, exactLinear);
- // compute finite difference approximation
+ // compute finite difference approximation
- VectorType fdLin = conf;
+ auto fdLin = conf;
- VectorType forward = conf;
- VectorType backward = conf;
+ auto forward = conf;
+ auto backward = conf;
- for (size_t i=0; i<conf.size(); i++)
- for (int j=0; j<dim; j++) {
+ for (size_t i = 0; i < conf.size(); ++i)
+ for (int j = 0; j < dim; ++j) {
- forward[i][j] += eps;
- backward[i][j] -= eps;
+ forward[i][j] += eps;
+ backward[i][j] -= eps;
- double forE = material.energy(std::make_shared<GridFunc>(basis,forward));
- double backE = material.energy(std::make_shared<GridFunc>(basis,backward));
+ auto forE = material.energy(std::make_shared<GridFunc>(basis, forward));
+ auto backE = material.energy(std::make_shared<GridFunc>(basis, backward));
- fdLin[i][j] = (forE-backE)/(2*eps);
+ fdLin[i][j] = (forE - backE) / (2 * eps);
- forward[i][j] -= eps;
- backward[i][j] += eps;
- }
+ forward[i][j] -= eps;
+ backward[i][j] += eps;
+ }
- VectorType diff = fdLin;
- diff -= exactLinear;
- std::cout<<"Two norm "<<std::setprecision(12)<<diff.two_norm()<<"\n";
+ auto diff = fdLin;
+ diff -= exactLinear;
+ std::cout << "Two norm " << std::setprecision(12) << diff.two_norm() << "\n";
- return (diff.two_norm()<1e-4);
+ return (diff.two_norm() < 1e-4);
}
-
using namespace Dune;
-int main (int argc, char *argv[]) try
-{
-
- // parse eps
- double eps = 1e-6;
- if (argc==2)
- eps = atof(argv[1]);
- std::cout<<"eps "<<eps<<std::endl;
-
- typedef UGGrid<3> GridType;
+int main(int argc, char* argv[]) try {
+ // parse eps
+ double eps = 1e-6;
+ if (argc == 2) eps = atof(argv[1]);
+ std::cout << "eps " << eps << std::endl;
- // make sphere grid
- FieldVector<double,3> center(0);
- double radius = 1;
+ using Grid = UGGrid<3>;
- auto grid = makeSphere<GridType>(center,radius);
+ // make sphere grid
+ FieldVector<double, 3> center(0);
+ double radius = 1;
- grid->setRefinementType(GridType::COPY);
- for (int i=0; i<1; i++)
- grid->globalRefine(1);
+ auto grid = makeSphere<Grid>(center, radius);
- // Create the materials
- typedef P1NodalBasis<GridType::LeafGridView> P1Basis;
- P1Basis p1Basis(grid->leafGridView());
+ grid->setRefinementType(Grid::COPY);
+ for (int i = 0; i < 1; i++)
+ grid->globalRefine(1);
- // make a material
- //typedef MooneyRivlinMaterial<P1Basis> MaterialType;
- typedef NeoHookeanMaterial<P1Basis> MaterialType;
- //typedef OgdenMaterial<P1Basis> MaterialType;
- //typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
+ // Create the materials
+ using P1Basis = P1NodalBasis<Grid::LeafGridView>;
+ P1Basis p1Basis(grid->leafGridView());
- MaterialType material;
- //material.setup(p1Basis, 1e5,0.3, 100);
- material.setup(p1Basis, 1e5,0.3);
+ // make a material
+ /*
+ using MoonRivlinMaterial = MooneyRivlinMaterial<P1Basis>;
+ using MaterialType = AdolcMaterial<P1Basis>;
+ MoonRivlinMaterial materialM;
+ materialM.setup(p1Basis, 1e5, 0.3);
+ MaterialType material(p1Basis, materialM);
+ */
- bool passed = checkLinearization(p1Basis,material,eps);
- passed = checkHessian(p1Basis,material,eps) and passed;
+ using MaterialType = NeoHookeanMaterial<P1Basis>;
+ // typedef OgdenMaterial<P1Basis> MaterialType;
+ // typedef GeomExactStVenantMaterial<P1Basis> MaterialType;
- std::cout<<"Checking 2D case!\n";
+ MaterialType material(p1Basis, 1e5, 0.3);
+ //materialM.setup(p1Basis, 1e5, 0.3);
- typedef UGGrid<2> GridType2;
- FieldVector<double,2> cent2(0);
+ bool passed = checkLinearization(p1Basis, material, eps);
+ passed = checkHessian(p1Basis, material, eps) and passed;
- auto grid2 = makeRingSegment2D<GridType2>(cent2, 0.2, 1, 0, 1.5, 10, false);
+ std::cout << "Checking 2D case!\n";
- grid2->setRefinementType(GridType2::COPY);
- for (int i=0; i<1; i++)
- grid2->globalRefine(1);
+ using Grid2 = UGGrid<2>;
+ FieldVector<double, 2> cent2(0);
- // Create the materials
- typedef P1NodalBasis<GridType2::LeafGridView> P1Basis2;
- P1Basis2 p1Basis2(grid2->leafGridView());
+ auto grid2 = makeRingSegment2D<Grid2>(cent2, 0.2, 1, 0, 1.5, 10, false);
- // make a material
- //typedef MooneyRivlinMaterial<P1Basis2> MaterialType2;
- typedef NeoHookeanMaterial<P1Basis2> MaterialType2;
- //typedef OgdenMaterial<P1Basis2> MaterialType2;
- //typedef GeomExactStVenantMaterial<P1Basis2> MaterialType2;
+ grid2->setRefinementType(Grid2::COPY);
+ for (int i = 0; i < 1; i++)
+ grid2->globalRefine(1);
- MaterialType2 material2;
- material2.setup(p1Basis2, 1e5, 0.3);
- //material2.setup(p1Basis2, 1e5, 0.3, 100);
+ // Create the materials
+ using P1Basis2 = P1NodalBasis<Grid2::LeafGridView>;
+ P1Basis2 p1Basis2(grid2->leafGridView());
- passed = checkLinearization(p1Basis2,material2,eps) and passed;
- passed = checkHessian(p1Basis2,material2,eps) and passed;
+ // make a material
+ //using MaterialType2 = MooneyRivlinMaterial<P1Basis2>;
+ using MaterialType2 = NeoHookeanMaterial<P1Basis2>;
+ // typedef OgdenMaterial<P1Basis2> MaterialType2;
+ // typedef GeomExactStVenantMaterial<P1Basis2> MaterialType2;
- return !passed;
+ //MaterialType2 material2(p1Basis2, 1e5, 0.3);
+ MaterialType2 material2(p1Basis2, 1e5, 0.3);
-} catch(Exception e) {
+ passed = checkLinearization(p1Basis2, material2, eps) and passed;
+ passed = checkHessian(p1Basis2, material2, eps) and passed;
-std::cout<<e<<std::endl;
+ return !passed;
+} catch (Exception e) {
+ std::cout << e << std::endl;
+ return 1;
}