dune/elasticity/materials/neumannenergy.hh

@@ -168,8 +168,8 @@ public:
 
 private:
   /** \brief The Neumann boundary */
-  const std::shared_ptr<BoundaryPatch<GridView>>
-    DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::NeumannEnergy") neumannBoundary_;
+  [[deprecated("Use dune-functions powerBases with LocalView concept. See Elasticity::NeumannEnergy")]]
+  const std::shared_ptr<BoundaryPatch<GridView>> neumannBoundary_;
 
   /** \brief The function implementing the Neumann data */
   const std::shared_ptr<std::function<Dune::FieldVector<double,dim>(Dune::FieldVector<ctype,dim>)>> neumannFunction_;

dune/elasticity/materials/stvenantkirchhoffenergy.hh

@@ -7,7 +7,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class DUNE_DEPRECATED_MSG("Use Elasticity::LocalIntegralEnergy with Elasticity::StVenantKirchhoffDensity")
+class [[deprecated("Use Elasticity::LocalIntegralEnergy with Elasticity::StVenantKirchhoffDensity")]]
 StVenantKirchhoffEnergy
   : public Elasticity::LocalIntegralEnergy<GridView,LocalFiniteElement,field_type>
 {

dune/elasticity/materials/sumenergy.hh

@@ -90,9 +90,8 @@ public:
   }
 
 private:
-
-  std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > >
-    DUNE_DEPRECATED_MSG("Use dune-functions powerBases with LocalView concept. See Elasticity::SumEnergy") a_;
+  [[deprecated("Use dune-functions powerBases with LocalView concept. See Elasticity::SumEnergy")]]
+  std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > a_;
 
   std::shared_ptr<LocalEnergy<GridView,LocalFiniteElement,std::vector<FieldVector<field_type,dim> > > > b_;
 

problems/film-on-substrate-deformation.parset

+#############################################
+#  Grid parameters
+#############################################
+
+structuredGrid = true
+lower = 0 0 0
+
+# whole experiment: 45 mm x 10 mm x 2 mm, scaling with 10^7 such that the thickness, which is around 100 nm, so 100x10^-9 = 10^-7 is equal to 1.
+# upper = 45e4 10e4 2e4
+# using only a section of the whole experiment as deformed grid to start with for dune-gfe: use much smaller dimensions!
+upper = 600 200 200
+
+elements = 15 5 5
+
+# Number of grid levels
+numLevels = 3
+
+adaptiveRefinement = true
+
+#############################################
+#  Solver parameters
+#############################################
+
+# Number of homotopy steps for the Dirichlet boundary conditions
+numHomotopySteps = 1
+
+# Tolerance of the trust region solver
+tolerance = 1e-6
+
+# Max number of steps of the trust region solver
+maxTrustRegionSteps = 500
+
+# Initial trust-region radius
+initialTrustRegionRadius = 0.1
+
+# Number of multigrid iterations per trust-region step
+numIt = 1000
+
+# Number of presmoothing steps
+nu1 = 3
+
+# Number of postsmoothing steps
+nu2 = 3
+
+# Damping for the smoothers of the multigrid solver
+damping = 1.0
+
+# Number of coarse grid corrections
+mu = 1
+
+# Number of base solver iterations
+baseIt = 100
+
+# Tolerance of the multigrid solver
+mgTolerance = 1e-7
+
+# Tolerance of the base grid solver
+baseTolerance = 1e-8
+
+############################
+#   Material parameters
+############################
+
+energy = mooneyrivlin # stvenantkirchhoff, neohooke, hencky, exphencky or mooneyrivlin
+
+[materialParameters]
+
+## Lame parameters for stvenantkirchhoff, E = mu(3*lambda + 2*mu)/(lambda + mu)
+#mu = 2.7191e+4
+#lambda = 4.4364e+4
+
+#mooneyrivlin_a = 2.7191e+6
+#mooneyrivlin_b = 2.7191e+6
+#mooneyrivlin_c = 2.7191e+6
+
+#mooneyrivlin_10 = -7.28e+5 #182 20:1
+#mooneyrivlin_01 = 9.17e+5
+#mooneyrivlin_20 = 1.23e+5
+#mooneyrivlin_02 = 8.15e+5
+#mooneyrivlin_11 = -5.14e+5
+
+#mooneyrivlin_10 = -3.01e+6 #182 2:1
+#mooneyrivlin_01 = 3.36e+6
+#mooneyrivlin_20 = 5.03e+6
+#mooneyrivlin_02 = 13.1e+6
+#mooneyrivlin_11 = -15.2e+6
+
+mooneyrivlin_10 = -1.67e+6 #184 2:1
+mooneyrivlin_01 = 1.94e+6
+mooneyrivlin_20 = 2.42e+6
+mooneyrivlin_02 = 6.52e+6
+mooneyrivlin_11 = -7.34e+6
+
+mooneyrivlin_30 = 0
+mooneyrivlin_21 = 0
+mooneyrivlin_12 = 0
+mooneyrivlin_03 = 0
+
+# volume-preserving parameter 
+mooneyrivlin_k = 90e+6
+
+# How to choose the volume-preserving parameter?
+# We need a stretch of 30%  (45e4 10e4 2e4 in x-direction, so a stretch of 45e4*0.3 = 13.5e4)
+# 184 2:1, mooneyrivlin_k = 90e+6 approximately (depends also on the number of grid levels!) and mooneyrivlin_energy = square or log, neumannValues = 27e4 0 0
+
+mooneyrivlin_energy = square # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
+
+# ciarlet: Fomula from "Ciarlet: Three-Dimensional Elasticity", here no penalty term is 
+# log: Generalized Rivlin model or polynomial hyperelastic model, using  0.5*mooneyrivlin_k*log(det(∇φ)) as the volume-preserving penalty term
+# square: Generalized Rivlin model or polynomial hyperelastic model, using 0.5*mooneyrivlin_k*(det(∇φ)-1)² as the volume-preserving penalty term
+
+[]
+
+#############################################
+#  Boundary values
+#############################################
+
+problem = identity-deformation
+
+###  Python predicate specifying all Dirichlet grid vertices
+# x is the vertex coordinate
+dirichletVerticesPredicate = "x[0] < 0.01"
+
+###  Python predicate specifying all neumannValues grid vertices
+# x is the vertex coordinate
+neumannVerticesPredicate = "x[0] > 599.99"
+
+###  Neumann values
+neumannValues = 27e4 0 0
+
+# Initial deformation
+initialDeformation = "[x[0], x[1], x[2]]"

problems/finite-strain-elasticity-bending.parset

@@ -56,7 +56,6 @@ baseTolerance = 1e-8
 
 energy = stvenantkirchhoff
 
-## For the Wriggers L-shape example
 [materialParameters]
 
 # Lame parameters
@@ -79,7 +78,7 @@ mooneyrivlin_12 = 0
 mooneyrivlin_03 = 0
 
 mooneyrivlin_k = 1e+6	# volume-preserving parameter
-mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute the Mooney-Rivlin-Energy
+mooneyrivlin_energy = log # log, square or Ciarlet; different ways to compute the Mooney-Rivlin-Energy
 
 []
 
@@ -87,7 +86,7 @@ mooneyrivlin_energy = log # log, square or ciarlet; different ways to compute th
 #  Boundary values
 #############################################
 
-problem = wriggers-l-shape
+problem = identity-deformation
 
 ###  Python predicate specifying all Dirichlet grid vertices
 # x is the vertex coordinate

problems/identity-deformation-dirichlet-values.py

+class DirichletValues:
+    def __init__(self, homotopyParameter):
+        self.homotopyParameter = homotopyParameter
+
+    # The identity deformation of 3d classical materials
+    def dirichletValues(self, x):
+        # Clamp the shape in its reference configuration
+        return [x[0], x[1], x[2]]

src/CMakeLists.txt

-if(ADOLC_FOUND AND IPOPT_FOUND AND PYTHONLIBS_FOUND AND dune-uggrid_FOUND)
-  set(programs finite-strain-elasticity
-               linear-elasticity)
+# Note: PYTHONLIBS_FOUND is only for backwards compatibility with dune-fufem 2.7 and can be removed
+#       in the next release
+if(ADOLC_FOUND AND IPOPT_FOUND AND ( Python3_FOUND OR PYTHONLIBS_FOUND )  AND dune-uggrid_FOUND)
+  set(programs finite-strain-elasticity)
+              #linear-elasticity depends on Amiramesh which is dropped since DUNE 2.8
 
   foreach(_program ${programs})
     add_executable(${_program} ${_program}.cc)

src/finite-strain-elasticity.cc

@@ -55,6 +55,8 @@ int main (int argc, char *argv[]) try
   // initialize MPI, finalize is done automatically on exit
   Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
 
+  if (mpiHelper.rank()==0)
+    std::cout << "ORDER = " << order << std::endl;
   // Start Python interpreter
   Python::start();
   Python::Reference main = Python::import("__main__");
@@ -71,6 +73,8 @@ int main (int argc, char *argv[]) try
 
   // parse data file
   ParameterTree parameterSet;
+  if (argc < 2)
+    DUNE_THROW(Exception, "Usage: ./finite-strain-elasticity <parameter file>");
 
   ParameterTreeParser::readINITree(argv[1], parameterSet);
 
@@ -158,21 +162,16 @@ int main (int argc, char *argv[]) try
   // Make Python function that computes which vertices are on the Dirichlet boundary,
   // based on the vertex positions.
   std::string lambda = std::string("lambda x: (") + parameterSet.get<std::string>("dirichletVerticesPredicate") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonDirichletVertices(Python::evaluate(lambda));
+  auto pythonDirichletVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   // Same for the Neumann boundary
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("neumannVerticesPredicate", "0") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, bool> pythonNeumannVertices(Python::evaluate(lambda));
+  auto pythonNeumannVertices = Python::make_function<bool>(Python::evaluate(lambda));
 
   for (auto&& v : vertices(gridView))
   {
-    bool isDirichlet;
-    pythonDirichletVertices.evaluate(v.geometry().corner(0), isDirichlet);
-    dirichletVertices[indexSet.index(v)] = isDirichlet;
-
-    bool isNeumann;
-    pythonNeumannVertices.evaluate(v.geometry().corner(0), isNeumann);
-    neumannVertices[indexSet.index(v)] = isNeumann;
+    dirichletVertices[indexSet.index(v)] = pythonDirichletVertices(v.geometry().corner(0));
+    neumannVertices[indexSet.index(v)] = pythonNeumannVertices(v.geometry().corner(0));
   }
 
   BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
@@ -200,7 +199,7 @@ int main (int argc, char *argv[]) try
   SolutionType x(basis.size());
 
   lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
-  PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > pythonInitialDeformation(Python::evaluate(lambda));
+  auto pythonInitialDeformation = Python::make_function<FieldVector<double,dim> >(Python::evaluate(lambda));
 
   Dune::Functions::interpolate(basis, x, pythonInitialDeformation);
 
@@ -326,7 +325,7 @@ int main (int argc, char *argv[]) try
     Python::Reference dirichletValuesPythonObject = C(homotopyParameter);
 
     // Extract object member functions as Dune functions
-    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> >   dirichletValues(dirichletValuesPythonObject.get("dirichletValues"));
+    auto dirichletValues = Python::make_function<FieldVector<double, dim> >(dirichletValuesPythonObject.get("dirichletValues"));
 
     Dune::Functions::interpolate(basis, x, dirichletValues, dirichletDofs);
 

test/CMakeLists.txt

@@ -8,5 +8,8 @@ dune_add_test(SOURCES materialtest
 add_dune_ug_flags(materialtest)
 add_dune_adolc_flags(materialtest)
 
+dune_add_test(SOURCES mooneyrivlintest
+              CMAKE_GUARD ADOLC_FOUND)
+
 dune_add_test(SOURCES localhyperdualstiffnesstest
               CMAKE_GUARD ADOLC_FOUND)

test/localhyperdualstiffnesstest.cc

@@ -85,10 +85,10 @@ TestSuite hyperdualEqualsADOL_C()
 
   Timer timer;
   hyperdualAssembler.assembleGradientAndHessian(x,hyperdualGradient,hyperdualHessian);
-  std::cout <<  dim <<"D: hyperdual assembler took " << timer.elapsed() << " sec." << std::endl;
+  std::cout <<  dim <<"D: hyperdual hessian and gradient assembler took " << timer.elapsed() << " sec." << std::endl;
   timer.reset();
   adolcAssembler.assembleGradientAndHessian(x,adolcGradient,adolcHessian);
-  std::cout << dim << "D: ADOL-C assembler took " << timer.elapsed() << " sec." << std::endl;
+  std::cout << dim << "D: ADOL-C hessian and gradient assembler took " << timer.elapsed() << " sec." << std::endl;
 
 
   hyperdualHessian -= adolcHessian;
@@ -97,6 +97,18 @@ TestSuite hyperdualEqualsADOL_C()
   // check the relative error
   t.check(hyperdualHessian.frobenius_norm()/adolcHessian.frobenius_norm() < 1e-12 && hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
 
+
+  timer.reset();
+  hyperdualAssembler.assembleGradient(x,hyperdualGradient);
+  std::cout <<  dim <<"D: hyperdual gradient assembler took " << timer.elapsed() << " sec." << std::endl;
+  timer.reset();
+  adolcAssembler.assembleGradient(x,adolcGradient);
+  std::cout << dim << "D: ADOL-C gradient assembler took " << timer.elapsed() << " sec." << std::endl;
+  hyperdualGradient -= adolcGradient;
+
+  // check the relative error
+  t.check(hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
+
   return t;
 
 }

test/mooneyrivlintest.cc

+#include <config.h>
+#include <iomanip>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dune/elasticity/assemblers/localadolcstiffness.hh>
+#include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
+#include <dune/elasticity/materials/mooneyrivlindensity.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/grid/utility/structuredgridfactory.hh>
+#include <dune/grid/uggrid.hh>
+
+// grid dimension
+const int dim = 3;
+
+//order of integration
+const int order = 1;
+
+using namespace Dune;
+
+//differentiation method: ADOL-C
+typedef adouble ValueType;
+
+using namespace Dune;
+
+typedef BlockVector<FieldVector<double,dim> > VectorType;
+typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
+
+template <class Basis>
+int assembleAndCompare (const Basis basis, const ParameterTree parameters, const VectorType x,
+                        double expectedEnergy, double expectedGradientTwoNorm, double expectedGradientInfinityNorm, double expectedMatrixFrobeniusNorm) {
+  using LocalView = typename Basis::LocalView;
+
+  std::string mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
+
+  auto elasticDensity = std::make_shared<Elasticity::MooneyRivlinDensity<dim,ValueType>>(parameters);
+  auto elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(elasticDensity);
+  auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(elasticEnergy);
+  Elasticity::FEAssembler<Basis,VectorType> assembler(basis, localADOLCStiffness);
+
+  //////////////////////////////////////////////////////////////
+  //  Compute the energy, assemble and compare!
+  //////////////////////////////////////////////////////////////
+
+  VectorType gradient;
+  MatrixType hessianMatrix;
+  double energy = assembler.computeEnergy(x);
+
+  if ( std::abs(energy - expectedEnergy)/expectedEnergy > 1e-3)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated energy is " << energy << " but '" << expectedEnergy << "' was expected!" << std::endl;
+      return 1;
+  }
+
+  assembler.assembleGradientAndHessian(x, gradient, hessianMatrix, true);
+
+  double gradientTwoNorm = gradient.two_norm();
+  double gradientInfinityNorm = gradient.infinity_norm();
+  double matrixFrobeniusNorm = hessianMatrix.frobenius_norm();
+  if ( std::abs(gradientTwoNorm - expectedGradientTwoNorm)/expectedGradientTwoNorm > 1e-4 || std::abs(gradientInfinityNorm - expectedGradientInfinityNorm)/expectedGradientInfinityNorm > 1e-8)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated gradient max norm is " << gradientInfinityNorm << " but '" << expectedGradientInfinityNorm << "' was expected!" << std::endl;
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated gradient norm is " << gradientTwoNorm << " but '" << expectedGradientTwoNorm << "' was expected!" << std::endl;
+      return 1;
+  }
+
+  if ( std::abs(matrixFrobeniusNorm - expectedMatrixFrobeniusNorm)/expectedMatrixFrobeniusNorm > 1e-4)
+  {
+      std::cerr << std::setprecision(9);
+      std::cerr << "In " << mooneyrivlin_energy << ": Calculated matrix norm is " << matrixFrobeniusNorm << " but '" << expectedMatrixFrobeniusNorm << "' was expected!" << std::endl;
+      return 1;
+  }
+
+  return 0;
+}
+
+int main (int argc, char *argv[])
+{
+
+  [[maybe_unused]] Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+  // ///////////////////////////////////////
+  //    Create the grid
+  // ///////////////////////////////////////
+  using GridType = UGGrid<dim>;
+  auto grid = StructuredGridFactory<GridType>::createCubeGrid({0,0,0}, {1,1,1}, {3,3,3});
+  int refine = 3;
+  while(refine > 0) {
+    for (auto&& e : elements(grid->leafGridView())){
+      bool refine = false;
+      for (int i = 0; i < e.geometry().corners(); i++) {
+          refine = refine || (e.geometry().corner(i)[0] > 0.99);
+      }
+      grid->mark(refine ? 1 : 0, e);
+    }
+    grid->adapt();
+    refine--;
+  }
+  grid->loadBalance();
+
+  using GridView = GridType::LeafGridView;
+  GridView gridView = grid->leafGridView();
+
+  /////////////////////////////////////////////////////////////////////////
+  //  Create a power basis
+  /////////////////////////////////////////////////////////////////////////
+
+  using namespace Functions::BasisFactory;
+  auto basis = makeBasis(
+    gridView,
+    power<dim>(
+      lagrange<order>(),
+      blockedInterleaved()
+  ));
+
+  //////////////////////////////////////////////////////////////
+  //  Prepare the iterate x where we want to assemble
+  //////////////////////////////////////////////////////////////
+  VectorType x(basis.size());
+  Functions::interpolate(basis, x, [](FieldVector<double,dim> x){
+    FieldVector<double,dim> y = {0.5*x[0]*x[0], 0.5*x[1] + 0.2*x[0], 4*std::abs(std::sin(x[2]))};
+    return y;
+  });
+
+  ParameterTree parameters;
+  parameters["mooneyrivlin_10"] = "1.67e+6";
+  parameters["mooneyrivlin_01"] = "1.94e+6";
+  parameters["mooneyrivlin_20"] = "2.42e+6";
+  parameters["mooneyrivlin_02"] = "6.52e+6";
+  parameters["mooneyrivlin_11"] = "-7.34e+6";
+  parameters["mooneyrivlin_30"] = "0";
+  parameters["mooneyrivlin_03"] = "0";
+  parameters["mooneyrivlin_21"] = "0";
+  parameters["mooneyrivlin_12"] = "0";
+  parameters["mooneyrivlin_k"] = "75e+6";
+
+  parameters["mooneyrivlin_energy"] = "square";
+  double expectedEnergy = 167636683;
+  double expectedGradientTwoNorm = 2.10685704e+09;
+  double expectedGradientInfinityNorm = 589695526;
+  double expectedMatrixFrobeniusNorm = 1.61941167e+11;
+  int testSquare = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+  parameters["mooneyrivlin_energy"] = "log";
+  expectedEnergy = 181180273;
+  expectedGradientTwoNorm = 2.29399362e+09;
+  expectedGradientInfinityNorm = 648551554;
+  expectedMatrixFrobeniusNorm = 1.67152642e+11;
+  int testLog = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+  parameters["mooneyrivlin_energy"] = "ciarlet";
+  parameters["mooneyrivlin_a"] = "2.42e+6";
+  parameters["mooneyrivlin_b"] = "6.52e+6";
+  parameters["mooneyrivlin_c"] = "-7.34e+6";
+  expectedEnergy = 76302830.4;
+  expectedGradientTwoNorm = 40670527.3;
+  expectedGradientInfinityNorm = 11116511.8;
+  expectedMatrixFrobeniusNorm = 2.1978108e+09;
+  int testCiarlet = assembleAndCompare(basis, parameters, x, expectedEnergy, expectedGradientTwoNorm, expectedGradientInfinityNorm, expectedMatrixFrobeniusNorm);
+
+  return testCiarlet + testLog + testSquare;
+}