Print error message if no parameter file is given

dune/elasticity/common/trustregionsolver.cc

@@ -120,6 +120,7 @@ setup(const typename BasisType::GridView::Grid& grid,
 
     auto mmgStep = std::make_shared< MonotoneMGStep<MatrixType, CorrectionType> >();
 
+    mmgStep->setVerbosity(NumProc::REDUCED);
     mmgStep->setMGType(mu, nu1, nu2);
     mmgStep->ignoreNodes_ = &dirichletNodes;
     mmgStep->setBaseSolver(baseSolver);

dune/elasticity/common/trustregionsolver.hh

@@ -113,6 +113,8 @@ public:
 
     SolutionType getSol() const {return x_;}
 
+    std::string iterateNamePrefix_;
+
 protected:
 
     /** \brief The grid */

dune/elasticity/materials/coshenergy.hh

@@ -11,7 +11,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class CosHEnergy
+class CosHEnergy final
   : public Elasticity::LocalEnergy<GridView,
                                    LocalFiniteElement,
                                    std::vector<Dune::FieldVector<field_type,GridView::dimension> > >

dune/elasticity/materials/expacoshenergy.hh

@@ -11,7 +11,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class ExpACosHEnergy
+class ExpACosHEnergy final
   : public Elasticity::LocalEnergy<GridView,
                                    LocalFiniteElement,
                                    std::vector<Dune::FieldVector<field_type,GridView::dimension> > >

dune/elasticity/materials/neffenergy.hh

@@ -11,7 +11,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffEnergy
+class NeffEnergy final
   : public Elasticity::LocalEnergy<GridView,
                                    LocalFiniteElement,
                                    std::vector<Dune::FieldVector<field_type,GridView::dimension> > >

dune/elasticity/materials/nefflogenergy.hh

@@ -11,7 +11,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffLogEnergy
+class NeffLogEnergy final
   : public Elasticity::LocalEnergy<GridView,
                                    LocalFiniteElement,
                                    std::vector<Dune::FieldVector<field_type,GridView::dimension> > >

dune/elasticity/materials/nefflogenergy2.hh

@@ -11,7 +11,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffLogEnergy2
+class NeffLogEnergy2 final
   : public Elasticity::LocalEnergy<GridView,
                                    LocalFiniteElement,
                                    std::vector<Dune::FieldVector<field_type,GridView::dimension> > >

dune/elasticity/materials/neffmagicenergy.hh

-#ifndef DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
-#define DUNE_ELASTICITY_MATERIALS_NEFFMAGICENERGY_HH
-
-#include <dune/common/fmatrix.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/geometry/quadraturerules.hh>
-
-#include <dune/elasticity/assemblers/localfestiffness.hh>
-
-namespace Dune::Elasticity {
-
-template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffMagicEnergy
-  : public Elasticity::LocalEnergy<GridView,
-                                   LocalFiniteElement,
-                                   std::vector<Dune::FieldVector<field_type,GridView::dimension> > >
-{
-  // grid types
-  typedef typename GridView::Grid::ctype DT;
-  typedef typename GridView::template Codim<0>::Entity Entity;
-
-  // some other sizes
-  enum {gridDim=GridView::dimension};
-  enum {dim=GridView::dimension};
-
-public:
-
-  /** \brief Constructor with a set of material parameters
-   * \param parameters The material parameters
-   */
-  NeffMagicEnergy(const Dune::ParameterTree& parameters)
-  {}
-
-  /** \brief Assemble the energy for a single element */
-  field_type energy (const Entity& e,
-                     const LocalFiniteElement& localFiniteElement,
-                     const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const;
-};
-
-template <class GridView, class LocalFiniteElement, class field_type>
-field_type
-NeffMagicEnergy<GridView, LocalFiniteElement, field_type>::
-energy(const Entity& element,
-       const LocalFiniteElement& localFiniteElement,
-       const std::vector<Dune::FieldVector<field_type, gridDim> >& localConfiguration) const
-{
-  assert(element.type() == localFiniteElement.type());
-
-  field_type energy = 0;
-
-  // store gradients of shape functions and base functions
-  std::vector<Dune::FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
-  std::vector<Dune::FieldVector<DT,gridDim> > gradients(localFiniteElement.size());
-
-  int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
-                                               : localFiniteElement.localBasis().order() * gridDim;
-
-  const auto& quad = Dune::QuadratureRules<DT, gridDim>::rule(element.type(), quadOrder);
-
-  for (const auto& qp : quad)
-  {
-    const DT integrationElement = element.geometry().integrationElement(qp.position());
-
-    const auto jacobianInverseTransposed = element.geometry().jacobianInverseTransposed(qp.position());
-
-    // get gradients of shape functions
-    localFiniteElement.localBasis().evaluateJacobian(qp.position(), referenceGradients);
-
-    // compute gradients of base functions
-    for (size_t i=0; i<gradients.size(); ++i)
-      jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-    Dune::FieldMatrix<field_type,gridDim,gridDim> derivative(0);
-    for (size_t i=0; i<gradients.size(); i++)
-      for (int j=0; j<gridDim; j++)
-        derivative[j].axpy(localConfiguration[i][j], gradients[i]);
-
-    auto det = derivative.determinant();
-    auto normSquared = derivative.frobenius_norm2();
-
-    auto K = normSquared / (2*det);
-
-    using std::acosh;
-
-    auto energyDensity = K + std::sqrt(K*K-1) - acosh(K) + log(det);
-
-    energy += qp.weight() * integrationElement * energyDensity;
-  }
-
-  return energy;
-}
-
-}  // namespace Dune::Elasticity
-
-#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_EXPACOSHENERGY_HH

dune/elasticity/materials/neffreducedenergy.hh

@@ -11,7 +11,7 @@
 namespace Dune {
 
 template<class GridView, class LocalFiniteElement, class field_type=double>
-class NeffReducedEnergy
+class NeffReducedEnergy final
   : public Elasticity::LocalEnergy<GridView,
                                    LocalFiniteElement,
                                    std::vector<Dune::FieldVector<field_type,GridView::dimension> > >

grids/circle2ndorder.msh

@@ -6,23 +6,23 @@ $Nodes
 1 -1 0 0
 2 1 0 0
 3 0 0 0
-4 -0.5000000000004422 -0.8660254037841834 0
-5 0.5000000000016158 -0.8660254037835057 0
-6 -0.8660254037845663 -0.4999999999997789 0
-7 6.776189039135044e-13 -1 0
-8 0.8660254037849052 -0.4999999999991919 0
-9 0.5000000000004422 0.8660254037841834 0
-10 -0.5000000000016158 0.8660254037835057 0
-11 0.8660254037845663 0.4999999999997789 0
-12 -6.776189039135044e-13 1 0
-13 -0.8660254037849052 0.4999999999991919 0
-14 -0.02222222222229404 0.03849001794593358 0
-15 -0.261111111111955 0.4522577108647196 0
-16 0.2388888888890741 0.4522577108650584 0
-17 0.488888888888853 0.01924500897296679 0
-18 -0.2611111111113681 -0.4137676929191249 0
-19 -0.511111111111147 0.01924500897296679 0
-20 0.2388888888896609 -0.4137676929187861 0
+4 -0.5 -0.8660254037844386 0
+5  0.5 -0.8660254037844386 0
+6 -0.8660254037844386 -0.5 0
+7  0 -1 0
+8  0.8660254037844386 -0.5 0
+9  0.5 0.8660254037844386 0
+10 -0.5 0.8660254037844386 0
+11 0.8660254037844386 0.5 0
+12 0 1 0
+13 -0.8660254037844386 0.5 0
+14 0 0 0
+15 -0.25 0.4330127018922193 0
+16 0.25 0.4330127018922193 0
+17 0.5 0 0
+18 -0.25 -0.4330127018922193 0
+19 -0.5 0 0
+20 0.25 -0.4330127018922193 0
 $EndNodes
 $Elements
 15

problems/quasiconvexity-circle-expacosh.parset

@@ -6,8 +6,13 @@ structuredGrid = false
 path = /home/sander/dune/dune-elasticity/grids/
 gridFile = circle2ndorder.msh
 
+#structuredGrid = true
+#lower = 0 0
+#upper = 1 1
+#elements = 2 2
+
 # Number of grid levels
-numLevels = 7
+numLevels = 6
 
 #############################################
 #  Solver parameters
@@ -72,7 +77,9 @@ dirichletVerticesPredicate = "True"
 dirichletValues = affine
 
 # Initial deformation is (x[0] * sqrt(x0), x[1] / sqrt(x0))
-x0 = 4.38
+#x0 = 4.38
+x0 = 12.0186
+
 # Setting initialIterateFile overrides the x0 parameter
-initialIterateFile = experiment-6.data
-oldX0 = 10
+#initialIterateFile = experiment-6.data
+#oldX0 = 10

src/CMakeLists.txt

@@ -10,7 +10,9 @@ if(ADOLC_FOUND AND IPOPT_FOUND AND PYTHONLIBS_FOUND AND dune-uggrid_FOUND)
     add_dune_pythonlibs_flags(${_program})
     add_dune_ug_flags(${_program})
     add_dune_mpi_flags(${_program})
-#    target_compile_options(${_program} PRIVATE "-fpermissive")
+    add_dune_superlu_flags(${_program})
+    target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
+#    target_compile_options(${_program} PRIVATE "-g")
   endforeach()
 
   if (AMIRAMESH_FOUND)

src/finite-strain-elasticity.cc

@@ -73,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);
 

src/quasiconvexity-test.cc

 #include <config.h>
 
+#undef HAVE_DUNE_PARMG
+
 // Includes for the ADOL-C automatic differentiation library
 // Need to come before (almost) all others.
 #include <adolc/adouble.h>
@@ -35,12 +37,16 @@
 #include <dune/elasticity/assemblers/localadolcstiffness.hh>
 #include <dune/elasticity/assemblers/feassembler.hh>
 #include <dune/elasticity/materials/neffenergy.hh>
+#include <dune/elasticity/materials/burkholderenergy.hh>
 #include <dune/elasticity/materials/coshenergy.hh>
+#include <dune/elasticity/materials/dacorognaenergy.hh>
 #include <dune/elasticity/materials/expacoshenergy.hh>
+#include <dune/elasticity/materials/fungenergy.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
 #include <dune/elasticity/materials/nefflogenergy.hh>
 #include <dune/elasticity/materials/nefflogenergy2.hh>
-#include <dune/elasticity/materials/neffmagicenergy.hh>
 #include <dune/elasticity/materials/neffreducedenergy.hh>
+#include <dune/elasticity/materials/neffw2energy.hh>
 
 // grid dimension
 const int dim = 2;
@@ -49,6 +55,103 @@ const int order = 1;
 
 using namespace Dune;
 
+template<int dim, class field_type = double>
+struct KPowerDensity final
+: public Elasticity::LocalDensity<dim,field_type>
+{
+  KPowerDensity(const Dune::ParameterTree& parameters)
+  {
+    p_ = parameters.template get<double>("p");
+  }
+
+#warning The type of 'x' should by ctype, not field_type
+  field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+  {
+    auto detF = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*detF);
+
+    return std::pow(K, p_);
+  }
+
+  double p_;
+};
+
+
+template<int dim, class field_type = double>
+struct MagicDensity final
+: public Elasticity::LocalDensity<dim,field_type>
+{
+  MagicDensity(const Dune::ParameterTree& parameters)
+  {
+    c_ = 1.0;
+    if (parameters.hasKey("c"))
+      c_ = parameters.template get<double>("c");
+  }
+
+  /** \brief Implement space-dependent coefficients */
+#warning The parameter should actually use ctype, not field_type
+  double c(const Dune::FieldVector<field_type,dim>& x) const
+  {
+    return c_;
+#if 0  // Three circles
+    bool inCircle0 = (x-FieldVector<double,dim>{-0.5,0}).two_norm() <0.2;
+    bool inCircle1 = (x-FieldVector<double,dim>{0.3535,0.3535}).two_norm() <0.2;
+    bool inCircle2 = (x-FieldVector<double,dim>{0.3535,-0.3535}).two_norm() <0.2;
+
+    return (inCircle0 or inCircle1 or inCircle2) ? c_ : 1;
+#endif
+//    return (x.two_norm() < 0.9) ? 1 : c_;
+  }
+
+  field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+  {
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto K = normSquared / (2*det);
+
+    using std::acosh;
+
+    return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
+  }
+
+  double c_;
+};
+
+
+template<int dim, class field_type = double>
+struct VossDensity final
+: public Elasticity::LocalDensity<dim,field_type>
+{
+  VossDensity(const Dune::ParameterTree& parameters)
+  {
+    c_ = parameters.template get<double>("c");
+  }
+
+  field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+  {
+    auto det = F.determinant();
+    auto normSquared = F.frobenius_norm2();
+
+    auto KBold = normSquared / (2*det);
+
+    using std::acosh;
+
+    auto K = KBold + std::sqrt(KBold*KBold - 1);
+
+    auto h = (1.0/(2*K)) * (1 + K*K + (1+K)*std::sqrt(1 + 6*K + K*K)) - std::log(K)
+           + std::log(1 + 4*K + K*K + (1+K)*std::sqrt(1 + 6*K + K*K));
+
+    return h + c_ * std::log( det );
+  }
+
+  double c_;
+};
+
+
+
 int main (int argc, char *argv[]) try
 {
   // initialize MPI, finalize is done automatically on exit
@@ -70,6 +173,8 @@ int main (int argc, char *argv[]) try
 
   // parse data file
   ParameterTree parameterSet;
+  if (argc < 2)
+    DUNE_THROW(Exception, "Usage: ./quasiconvexity-text <parameter file>");
 
   ParameterTreeParser::readINITree(argv[1], parameterSet);
 
@@ -142,6 +247,7 @@ 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));
 
@@ -163,6 +269,20 @@ int main (int argc, char *argv[]) try
       for (int j=0; j<dim; j++)
         dirichletDofs[i][j] = true;
 
+#if 0
+  // HACK!
+#warning Teichmueller hack!
+  for (auto&& v : vertices(gridView))
+  {
+    if (v.geometry().corner(0).two_norm() < 1e-6)
+    {
+      auto idx = indexSet.index(v);
+      for (int j=0; j<dim; j++)
+        dirichletDofs[idx][j] = true;
+    }
+  }
+#endif
+
   // //////////////////////////
   //   Initial iterate
   // //////////////////////////
@@ -177,6 +297,8 @@ int main (int argc, char *argv[]) try
       blockedInterleaved()
   ));
 
+  if (parameterSet.hasKey("x0"))
+  {
     auto x0 = parameterSet.get<double>("x0");
 
     auto initialDeformation = [&x0](FieldVector<double,dim> x) -> FieldVector<double,dim>
@@ -185,10 +307,32 @@ int main (int argc, char *argv[]) try
     };
 
     Dune::Functions::interpolate(powerBasis, x, initialDeformation);
+  }
+
+  if (parameterSet.hasKey("initialIteratePython"))
+  {
+    // The python class that implements the Dirichlet values
+    Python::Module initialIterateModule = Python::import(parameterSet.get<std::string>("initialIteratePython"));
+
+    auto f = Python::Callable(initialIterateModule.get("initialIterate"));
+    PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > initialIteratePython(f);
+
+    Dune::Functions::interpolate(powerBasis, x, initialIteratePython);
+  }
 
   if (parameterSet.hasKey("initialIterateFile"))
   {
-    std::ifstream inFile(parameterSet.get<std::string>("initialIterateFile"), std::ios_base::binary);
+    std::string filename = parameterSet.get<std::string>("initialIterateFile");
+
+    // Guess the file format by looking at the file name suffix
+    auto dotPos = filename.rfind('.');
+    if (dotPos == std::string::npos)
+      DUNE_THROW(IOError, "Could not determine grid input file format");
+    std::string suffix = filename.substr(dotPos, filename.length()-dotPos);
+
+    if (suffix == ".data")
+    {
+      std::ifstream inFile(filename, std::ios_base::binary);
       if (not inFile)
         DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("initialIterateFile") << " could not be opened.");
       MatrixVector::Generic::readBinary(inFile, x);
@@ -198,6 +342,38 @@ int main (int argc, char *argv[]) try
       inFile.close();
     }
 
+    // Read CSV files as provided by Sid Kumar
+    if (suffix == ".csv")
+    {
+      std::ifstream inFile(filename, std::ios_base::binary);
+      if (not inFile)
+        DUNE_THROW(IOError, "File " << parameterSet.get<std::string>("initialIterateFile") << " could not be opened.");
+
+      int vertex = 0;
+      for (std::string line; std::getline(inFile, line);)
+      {
+        // The first line starts with 'x' -- skip it
+        if (line[0]=='x')
+          continue;
+
+        std::vector<std::string> result;
+        std::stringstream s_stream(line); //create string stream from the string
+        while(s_stream.good()) {
+          std::string substr;
+          getline(s_stream, substr, ','); //get first string delimited by comma
+          result.push_back(substr);
+        }
+
+        x[vertex][0] = atof(result[2].c_str());
+        x[vertex][1] = atof(result[3].c_str());
+
+        vertex++;
+      }
+
+      inFile.close();
+    }
+  }
+
   ////////////////////////////////////////////////////////
   //   Main homotopy loop
   ////////////////////////////////////////////////////////
@@ -238,10 +414,18 @@ int main (int argc, char *argv[]) try
 
     // Assembler using ADOL-C
     std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
-    std::shared_ptr<LocalFEStiffness<GridView,
+    std::shared_ptr<Elasticity::LocalEnergy<GridView,
                                      FEBasis::LocalView::Tree::FiniteElement,
                                      std::vector<Dune::FieldVector<adouble, dim> > > > elasticEnergy;
 
+    if (parameterSet.get<std::string>("energy") == "kpower")
+    {
+      auto density = std::make_shared<KPowerDensity<dim,adouble> >(materialParameters);
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+                                                                       FEBasis::LocalView::Tree::FiniteElement,
+                                                                       adouble> >(density);
+    }
+
     if (parameterSet.get<std::string>("energy") == "neff")
       elasticEnergy = std::make_shared<NeffEnergy<GridView,
                                                       FEBasis::LocalView::Tree::FiniteElement,
@@ -272,11 +456,42 @@ int main (int argc, char *argv[]) try
                                                          FEBasis::LocalView::Tree::FiniteElement,
                                                          adouble> >(materialParameters);
 
+    if (parameterSet.get<std::string>("energy") == "fung")
+      elasticEnergy = std::make_shared<Elasticity::FungEnergy<GridView,
+                                                  FEBasis::LocalView::Tree::FiniteElement,
+                                                  adouble> >(materialParameters);
+
     if (parameterSet.get<std::string>("energy") == "magic")
-      elasticEnergy = std::make_shared<Elasticity::NeffMagicEnergy<GridView,
+    {
+      auto density = std::make_shared<MagicDensity<dim,adouble> >(materialParameters);
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+                                                                       FEBasis::LocalView::Tree::FiniteElement,
+                                                                       adouble> >(density);
+    }
+
+    if (parameterSet.get<std::string>("energy") == "burkholder")
+      elasticEnergy = std::make_shared<Elasticity::BurkholderEnergy<GridView,
+                                                                    FEBasis::LocalView::Tree::FiniteElement,
+                                                                    adouble> >(materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "w2")
+      elasticEnergy = std::make_shared<Elasticity::NeffW2Energy<GridView,
+                                                                FEBasis::LocalView::Tree::FiniteElement,
+                                                                adouble> >(materialParameters);
+
+    if (parameterSet.get<std::string>("energy") == "dacorogna")
+      elasticEnergy = std::make_shared<Elasticity::DacorognaEnergy<GridView,
                                                                 FEBasis::LocalView::Tree::FiniteElement,
                                                                 adouble> >(materialParameters);
 
+    if (parameterSet.get<std::string>("energy") == "voss")
+    {
+      auto density = std::make_shared<VossDensity<dim,adouble> >(materialParameters);
+      elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
+                                                                       FEBasis::LocalView::Tree::FiniteElement,
+                                                                       adouble> >(density);
+    }
+
     if(!elasticEnergy)
       DUNE_THROW(Exception, "Error: Selected energy not available!");
 
@@ -306,7 +521,8 @@ int main (int argc, char *argv[]) try
                  mgTolerance,
                  mu, nu1, nu2,
                  baseIterations,
-                 baseTolerance
+                 baseTolerance,
+                 1.0
                 );
 
     ////////////////////////////////////////////////////////
@@ -328,9 +544,9 @@ int main (int argc, char *argv[]) try
 #else
 //     lambda = std::string("lambda x: (") + parameterSet.get<std::string>("initialDeformation") + std::string(")");
 //     PythonFunction<FieldVector<double,dim>, FieldVector<double,dim> > dirichletValues(Python::evaluate(lambda));
-    Python::Module module = Python::import(parameterSet.get<std::string>("dirichletValues"));
-    auto dirichletValues = Python::makeFunction<FieldVector<double,dim>(const FieldVector<double,dim>&)>(module.get("f"));
-    Dune::Functions::interpolate(powerBasis, x, dirichletValues, dirichletDofs);
+//     Python::Module module = Python::import(parameterSet.get<std::string>("dirichletValues"));
+//     auto dirichletValues = Python::makeFunction<FieldVector<double,dim>(const FieldVector<double,dim>&)>(module.get("f"));
+//     Dune::Functions::interpolate(powerBasis, x, dirichletValues, dirichletDofs);
 #endif
 
     //::Functions::interpolate(fufemFEBasis, x, dirichletValues, dirichletDofs);
@@ -339,6 +555,7 @@ int main (int argc, char *argv[]) try
     //   Solve!
     // /////////////////////////////////////////////////////
 
+    solver.iterateNamePrefix_ = "stage1-";
     solver.setInitialIterate(x);
     solver.solve();
 
@@ -348,6 +565,7 @@ int main (int argc, char *argv[]) try
     //   Output result
     /////////////////////////////////
 
+#if 0
     std::cout << "Total energy: " << assembler.computeEnergy(x) << std::endl;
 
     elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
@@ -363,7 +581,7 @@ int main (int argc, char *argv[]) try
     localADOLCStiffness.localEnergy_ = elasticEnergy.get();
 
     std::cout << "Energy part 2: " << assembler.computeEnergy(x) << std::endl;
-
+#endif
     // Compute the displacement
     auto displacement = x;
     displacement -= identity;
@@ -445,11 +663,14 @@ int main (int argc, char *argv[]) try
     }
 #endif
 
+    std::cout << "p: " << parameterSet.sub("materialParameters")["p"] << ",  "
+              << "K range: " << (*std::min_element(K.begin(), K.end()))
+              << ",  " << (*std::max_element(K.begin(), K.end())) << std::endl;
+
     auto deformationDeterminantFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, deformationDeterminant);
     auto KFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,1>>(p0Basis, K);
 
-    //  We need to subsample, because VTK cannot natively display real second-order functions
-    SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+    VTKWriter<GridView> vtkWriter(gridView);
     vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
     vtkWriter.addCellData(deformationDeterminantFunction,
                           VTK::FieldInfo("determinant", VTK::FieldInfo::Type::scalar, 1));
@@ -463,22 +684,24 @@ int main (int argc, char *argv[]) try
   MatrixVector::Generic::writeBinary(outFile, x);
   outFile.close();
 
+  exit(0);
+#if 0
   //////////////////////////////////////////////////////////////////////////////////
   //  Use the solution as initial iterate for the 'real' energy
   //////////////////////////////////////////////////////////////////////////////////
 
   ParameterTree neffMaterialParameters;
-  neffMaterialParameters["scaling"] = "1.15";
+  neffMaterialParameters["c"] = "1";
 
   // Careful: The following code computes the initial iterate given by the Python file,
   // and *assumes* that that file describes the homogeneous deformation!
   SolutionType homogeneous(feBasis.size());
     Dune::Functions::interpolate(powerBasis, homogeneous, initialDeformation);
 
-  std::shared_ptr<LocalFEStiffness<GridView,
+  std::shared_ptr<Elasticity::LocalEnergy<GridView,
                                           FEBasis::LocalView::Tree::FiniteElement,
                                           std::vector<Dune::FieldVector<adouble, dim> > > >
-    energy = std::make_shared<NeffEnergy<GridView,
+    energy = std::make_shared<Elasticity::NeffW2Energy<GridView,
                               FEBasis::LocalView::Tree::FiniteElement,
                               adouble> >(neffMaterialParameters);
 
@@ -510,6 +733,7 @@ int main (int argc, char *argv[]) try
                baseTolerance
               );
 
+  solver.iterateNamePrefix_ = "stage2-";
   solver.setInitialIterate(x);
   solver.solve();
 
@@ -522,7 +746,7 @@ int main (int argc, char *argv[]) try
   vtkWriter.clear();
   vtkWriter.addVertexData(displacementFunction2, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
   vtkWriter.write(resultPath + "final-result");
-
+#endif
 } catch (Exception& e) {
     std::cout << e.what() << std::endl;
 }