Add test for the automatic differentiation assemblers comparing them to

local assemblers of a nonlinear elasticity functional, derived by hand.

dune/fufem/test/CMakeLists.txt

 # Put your test in here if it needs access to external grids
 set(GRID_BASED_TESTS
+  adolctest
   basisgridfunctiontest
   basisinterpolatortest
   boundarypatchtest
@@ -66,7 +67,9 @@ if(PYTHONLIBS_FOUND)
   target_compile_definitions(dunepythontest PRIVATE -DDUNE_FUFEM_TEST_MODULE_PATH=\"${CMAKE_CURRENT_BINARY_DIR}\")
 endif()
 
-
+if (ADOLC_FOUND)
+    add_dune_adolc_flags(adolctest)
+endif()
 
 
 add_directory_test_target(_test_target)

dune/fufem/test/Makefile.am

 
 # list of tests to run
 TESTS = arithmetictest\
+		adolctest \
         basisgridfunctiontest \
         basisinterpolatortest \
         boundarypatchtest \
@@ -48,6 +49,11 @@ arithmetictest_CPPFLAGS = $(COMMON_CPPFLAGS)
 arithmetictest_LDADD = $(COMMON_LDADD)
 arithmetictest_LDFLAGS = $(COMMON_LDFLAGS)
 
+adolctest_SOURCES = adolctest.cc
+adolctest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS) $(ADOLC_CPPFLAGS)
+adolctest_LDADD = $(COMMON_LDADD) $(GRID_LDADD) $(ADOLC_LDFLAGS)
+adolctest_LDFLAGS = $(COMMON_LDFLAGS) $(GRID_LDFLAGS) $(ADOLC_LDFLAGS)
+
 basisgridfunctiontest_SOURCES = basisgridfunctiontest.cc
 basisgridfunctiontest_CPPFLAGS = $(COMMON_CPPFLAGS) $(GRID_CPPFLAGS)
 basisgridfunctiontest_LDADD = $(COMMON_LDADD) $(GRID_LDADD)

dune/fufem/test/adolctest.cc

+#include <config.h>
+#include <fstream>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/grid/uggrid.hh>
+//#include <dune/grid/io/file/dgfparser/gridptr.hh>
+//#include <dune/grid/io/file/dgfparser/dgfug.hh>
+
+#include <dune/istl/io.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/fufem/assemblers/localassemblers/adolclocalenergy.hh>
+#include <dune/fufem/assemblers/localassemblers/adolclinearizationassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/adolchessianassembler.hh>
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/assemblers/functionalassembler.hh>
+#include <dune/fufem/assemblers/operatorassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/geomexactstvenantkirchhofffunctionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/geomexactstvenantkirchhoffoperatorassembler.hh>
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
+#include <dune/fufem/symmetrictensor.hh>
+#include <dune/fufem/makesphere.hh>
+
+
+//! Local energy for a geometric exact St. Venant--Kirchhoff material
+template <class GridType, class LocalFiniteElement>
+class LocalElasticity : public Adolc::LocalEnergy<GridType, LocalFiniteElement,GridType::dimension>
+{
+  public:
+
+    enum {dim = GridType::dimension};
+    typedef typename GridType::ctype ctype;
+    typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::RangeFieldType field_type;
+
+    typedef Adolc::LocalEnergy<GridType, LocalFiniteElement, dim> Base;
+    typedef typename Base::CoefficientVectorType CoefficientVectorType;
+    typedef typename Base::ReturnType ReturnType;
+
+    typedef typename GridType::template Codim<0>::Entity Element;
+    typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits::JacobianType JacobianType;
+
+    LocalElasticity(field_type E, field_type nu) : E_(E), nu_(nu) {}
+
+    ReturnType energy(const Element& element, const LocalFiniteElement& lfe, 
+            const CoefficientVectorType& localCoeff) const
+    {
+        ReturnType energy(0);
+
+        // get quadrature rule
+        QuadratureRuleKey quad1(lfe);
+        QuadratureRuleKey quadKey = quad1.derivative().square().square();
+
+        const Dune::QuadratureRule<ctype, dim>& quad = QuadratureRuleCache<ctype, dim>::rule(quadKey);
+
+        // store gradients of shape functions and base functions
+        std::vector<JacobianType> referenceGradients(lfe.localBasis().size());
+
+        // the element geometry mapping
+        const typename Element::Geometry geometry = element.geometry();
+
+        // loop over quadrature points
+        for (size_t pt=0; pt < quad.size(); ++pt) {
+
+            // get quadrature point
+            const auto& quadPos = quad[pt].position();
+
+            // evaluate displacement gradient at the quadrature point
+
+            // get transposed inverse of Jacobian of transformation
+            const auto& invJacobian = geometry.jacobianInverseTransposed(quadPos);
+
+            // get integration factor
+            const ctype integrationElement = geometry.integrationElement(quadPos);
+
+            // get gradients of shape functions
+            lfe.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+            // compute gradient of the configuration
+            Dune::FieldMatrix<ReturnType,dim,dim> localGradient(0);
+            for (size_t k=0; k<referenceGradients.size(); ++k) {
+                Dune::FieldVector<ReturnType,dim> gradient(0);
+                invJacobian.umv(referenceGradients[k][0], gradient);
+                for (int i=0; i<dim; ++i)
+                    for (int j=0; j<dim; ++j)
+                        localGradient[i][j] += localCoeff[k][i] * gradient[j];
+            }
+
+            // Compute the nonlinear strain tensor from the deformation gradient
+            SymmetricTensor<dim,ReturnType> strain(0),stress(0);
+
+            computeNonlinearStrain(localGradient,strain);
+
+            // and the stress
+            stress = hookeTimesStrain(strain);
+
+            ReturnType z = quad[pt].weight()*integrationElement;
+            energy += (stress*strain)*z;
+        }
+
+        return 0.5*energy; 
+    }
+
+  private:
+    /** \brief Compute nonlinear strain tensor from a displacement gradient.
+     *
+     *  \param grad The gradient of the direction in which the linearisation is computed.
+     *  \param strain The tensor to store the strain in. 
+     */
+    static void computeNonlinearStrain(const Dune::FieldMatrix<ReturnType, dim, dim>& grad, 
+                                SymmetricTensor<dim, ReturnType>& strain)  {
+        strain = ReturnType(0);
+        for (int i=0; i<dim ; ++i) {
+            strain(i,i) +=grad[i][i];
+
+            for (int k=0;k<dim;++k)
+                strain(i,i) += 0.5*grad[k][i]*grad[k][i];
+
+            for (int j=i+1; j<dim; ++j) {
+                strain(i,j) += 0.5*(grad[i][j] + grad[j][i]);
+                for (int k=0;k<dim;++k)
+                    strain(i,j) += 0.5*grad[k][i]*grad[k][j];
+            }
+        }
+    }
+
+    //! Compute linear elastic stress from the strain
+    SymmetricTensor<dim,ReturnType> hookeTimesStrain(const SymmetricTensor<dim,ReturnType>& strain) const {
+
+        SymmetricTensor<dim,ReturnType> stress = strain;
+        stress *= E_/(1+nu_);
+
+        ReturnType f = E_*nu_ / ((1+nu_)*(1-2*nu_)) * strain.trace();
+        stress.addToDiag(f);
+
+        return stress;
+    }
+
+    field_type E_;
+    field_type nu_;
+
+};
+
+
+// grid dimension
+const int dim = 3;
+
+using namespace Dune;
+int main (int argc, char *argv[]) try
+{
+
+    typedef FieldVector<double,dim> FVector;
+    typedef FieldMatrix<double,dim,dim> FMatrix;
+    typedef BlockVector<FVector> VectorType;
+    typedef BCRSMatrix<FMatrix> MatrixType;
+    typedef std::vector<FVector> CoeffType;
+
+    // ///////////////////////////////////////
+    //    Create the grid
+    // ///////////////////////////////////////
+
+    typedef UGGrid<dim> GridType;
+    GridType* grid = new GridType;
+    makeSphere(*grid,FVector(0),0.2);
+    grid->globalRefine(3);
+    typedef GridType::LeafGridView GridView;
+    GridView gridView = grid->leafGridView();
+
+    typedef P1NodalBasis<GridView> FEBasis;
+    FEBasis feBasis(gridView);
+
+    // //////////////////////////
+    //   Initial iterate
+    // //////////////////////////
+
+    std::ifstream infile("adolctest.sol");
+    std::string line;
+
+    VectorType x(feBasis.size());
+    int i=0;
+    while (std::getline(infile, line))
+    {
+        std::istringstream iss(line);
+        for (int j=0; j<dim;j++)
+            iss >> x[i][j];
+        i++;
+    }
+
+    auto xGridFunc = make_shared<BasisGridFunction<FEBasis,VectorType> >(feBasis,x); 
+
+    // ////////////////////////////////////////////////////////////
+    //   Create an assembler for the energy functional
+    // ////////////////////////////////////////////////////////////
+
+    double E = 1e5;
+    double nu = 0.4;
+
+    typedef FEBasis::LocalFiniteElement Lfe;
+    LocalElasticity<GridType,Lfe> localEnergy(E,nu);
+    AdolcLinearizationAssembler<GridType,Lfe,FVector> linAdolcAssembler(localEnergy,*xGridFunc);
+
+    GeomExactStVenantKirchhoffFunctionalAssembler<GridType, Lfe> linPaperAssembler(E,nu);
+    linPaperAssembler.setConfiguration(xGridFunc);
+
+    FunctionalAssembler<FEBasis> functionalAssembler(feBasis);
+ 
+    VectorType adolcGradient(feBasis.size());
+    VectorType 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;
+
+    time.reset();
+    time.start();
+    functionalAssembler.assemble(linPaperAssembler,paperGrad);
+    std::cout<<"Paper&Pen functional assembler needed "<<time.stop()<<std::endl;
+
+    for (size_t i=0; i<adolcGradient.size();i++) {
+
+        FVector diff = adolcGradient[i];
+        diff -= paperGrad[i];
+
+        if (diff.two_norm()>1e-9)
+            DUNE_THROW(Dune::Exception,"Wrong local derivative, error is "<<diff.two_norm()); 
+
+    }
+    // ////////////////////////
+    // Test hessian assembler
+    // ////////////////////////
+
+    AdolcHessianAssembler<GridType,Lfe,Lfe,FMatrix> hessAdolcAssembler(localEnergy,*xGridFunc);
+    GeomExactStVenantKirchhoffOperatorAssembler<GridType, Lfe,Lfe> hessPaperAssembler(E,nu);
+
+    hessPaperAssembler.setConfiguration(xGridFunc);
+
+    OperatorAssembler<FEBasis,FEBasis> operatorAssembler(feBasis,feBasis);
+
+    MatrixType adolcHessian, paperHessian;
+
+    time.reset();
+    time.start();
+    operatorAssembler.assemble(hessAdolcAssembler, adolcHessian);
+    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;
+
+    for (size_t i=0; i<adolcHessian.N();i++) {
+
+        const auto& adolcRow = adolcHessian[i];
+        const auto& paperRow = paperHessian[i];
+
+        auto adolcIt = adolcRow.begin();
+        auto adolcEndIt = adolcRow.end();
+
+        auto paperIt = paperRow.begin();
+        auto paperEndIt = paperRow.end();
+
+        for (; adolcIt != adolcEndIt; ++adolcIt, ++paperIt) {
+
+            if (adolcIt.index() != paperIt.index())
+                DUNE_THROW(Dune::Exception,"Not the same sparsity pattern!"<<adolcIt.index()
+                                <<"!="<<paperIt.index());
+
+            FMatrix diff = *adolcIt;
+            diff -= *paperIt;
+
+            if (diff.frobenius_norm()>1e-8)
+                DUNE_THROW(Dune::Exception,"Wrong local hessian, error is "<<diff.frobenius_norm()); 
+
+        }
+        assert(paperIt==paperEndIt);
+    }
+
+    // //////////////////////////////
+} catch (Exception e) {
+
+    std::cout << e << std::endl;
+
+}