From 75fb1e6ec4ee6d3604d25cccd06558867360885a Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 14 May 2020 10:32:25 +0200
Subject: [PATCH] Use perturbation from an affine displacement as the unknown

Previously, the unknown was the displacement.  In this setting
it is difficult to implement periodic boundary conditions,
because the displacement itself is never periodic.
We therefore change the code to use the deviation from a
given affine transformation as the new variable.
That is surprisingly simple.  It is also what the ETH group does.
---
 dune/elasticity/common/trustregionsolver.cc |   1 +
 src/quasiconvexity-test.cc                  | 191 +++++++++++---------
 2 files changed, 107 insertions(+), 85 deletions(-)

diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index bdaaf12..1496cff 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -388,6 +388,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
             std::cout << "----------------------------------------------------" << std::endl;
             std::cout << "      Trust-Region Step Number: " << i
                       << ",     radius: " << trustRegion.radius()
+                      << std::setprecision(15)
                       << ",     energy: " << oldEnergy << std::endl;
             std::cout << "----------------------------------------------------" << std::endl;
         }
diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc
index 5491d41..7d2da5f 100644
--- a/src/quasiconvexity-test.cc
+++ b/src/quasiconvexity-test.cc
@@ -55,18 +55,19 @@ const int order = 1;
 
 using namespace Dune;
 
-template<int dim, class field_type = double>
+template<int dim, class field_type, class ctype=double>
 struct KPowerDensity final
-: public Elasticity::LocalDensity<dim,field_type>
+: public Elasticity::LocalDensity<dim,field_type,ctype>
 {
-  KPowerDensity(const Dune::ParameterTree& parameters)
+  KPowerDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
   {
     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
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
   {
+    auto F = derivative + F0_;
     auto detF = F.determinant();
     auto normSquared = F.frobenius_norm2();
 
@@ -75,15 +76,17 @@ struct KPowerDensity final
     return std::pow(K, p_);
   }
 
+  FieldMatrix<ctype,dim,dim> F0_;
   double p_;
 };
 
 
-template<int dim, class field_type = double>
+template<int dim, class field_type, class ctype=double>
 struct MagicDensity final
-: public Elasticity::LocalDensity<dim,field_type>
+: public Elasticity::LocalDensity<dim,field_type,ctype>
 {
-  MagicDensity(const Dune::ParameterTree& parameters)
+  MagicDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
   {
     c_ = 1.0;
     if (parameters.hasKey("c"))
@@ -91,8 +94,7 @@ struct MagicDensity final
   }
 
   /** \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
+  double c(const Dune::FieldVector<ctype,dim>& x) const
   {
     return c_;
 #if 0  // Three circles
@@ -105,8 +107,9 @@ struct MagicDensity final
 //    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
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
   {
+    auto F = derivative + F0_;
     auto det = F.determinant();
     auto normSquared = F.frobenius_norm2();
 
@@ -117,21 +120,24 @@ struct MagicDensity final
     return K + std::sqrt(K*K-1) - acosh(K) + c(x)*log(det);
   }
 
+  FieldMatrix<ctype,dim,dim> F0_;
   double c_;
 };
 
 
-template<int dim, class field_type = double>
+template<int dim, class field_type, class ctype=double>
 struct VossDensity final
-: public Elasticity::LocalDensity<dim,field_type>
+: public Elasticity::LocalDensity<dim,field_type,ctype>
 {
-  VossDensity(const Dune::ParameterTree& parameters)
+  VossDensity(const FieldMatrix<ctype,dim,dim>& F0, const Dune::ParameterTree& parameters)
+  : F0_(F0)
   {
     c_ = parameters.template get<double>("c");
   }
 
-  field_type operator()(const FieldVector<field_type,dim>& x, const FieldMatrix<field_type,dim,dim>& F) const
+  field_type operator()(const FieldVector<ctype,dim>& x, const FieldMatrix<field_type,dim,dim>& derivative) const
   {
+    auto F = derivative + F0_;
     auto det = F.determinant();
     auto normSquared = F.frobenius_norm2();
 
@@ -147,6 +153,7 @@ struct VossDensity final
     return h + c_ * std::log( det );
   }
 
+  FieldMatrix<ctype,dim,dim> F0_;
   double c_;
 };
 
@@ -210,7 +217,7 @@ int main (int argc, char *argv[]) try
     upper = parameterSet.get<FieldVector<double,dim> >("upper");
 
     std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
-    grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+    grid = StructuredGridFactory<GridType>::createSimplexGrid(lower, upper, elements);
 
   } else {
     std::string path                = parameterSet.get<std::string>("path");
@@ -297,17 +304,15 @@ int main (int argc, char *argv[]) try
       blockedInterleaved()
   ));
 
-  if (parameterSet.hasKey("x0"))
-  {
-    auto x0 = parameterSet.get<double>("x0");
+  // Set the underlying affine deformation
+  auto x0 = parameterSet.get<double>("x0");
 
-    auto initialDeformation = [&x0](FieldVector<double,dim> x) -> FieldVector<double,dim>
-    {
-      return {x[0]*sqrt(x0), x[1] / sqrt(x0)};
-    };
+  FieldMatrix<double,2,2> F0 = { {sqrt(x0), 0},
+                                 {       0, 1.0 / sqrt(x0)} };
 
-    Dune::Functions::interpolate(powerBasis, x, initialDeformation);
-  }
+  // The actual unknown is a perturbation of the affine transformation
+  // Start here with the zero function
+  std::fill(x.begin(), x.end(), FieldVector<double,2>(0));
 
   if (parameterSet.hasKey("initialIteratePython"))
   {
@@ -364,9 +369,23 @@ int main (int argc, char *argv[]) try
           result.push_back(substr);
         }
 
-        x[vertex][0] = atof(result[2].c_str());
-        x[vertex][1] = atof(result[3].c_str());
+        FieldVector<double,2> pos{atof(result[0].c_str()), atof(result[1].c_str())};
+        FieldVector<double,2> u{atof(result[2].c_str()), atof(result[3].c_str())};
+
+        //x[vertex] = u - (F0-Id)*pos;
+        FieldVector<double,2> phi = { u[0] - (F0[0][0]-1)*pos[0] -  F0[0][1]   *pos[1],
+                                      u[1] -  F0[1][0]   *pos[0] - (F0[1][1]-1)*pos[1] };
+
+        for (const auto& vert : vertices(gridView))
+          if ( (vert.geometry().corner(0)-pos).two_norm() < 1e-8)
+          {
+            x[gridView.indexSet().index(vert)] = phi;
+            std::cout << "pos: " << pos << ",  vertex: " << gridView.indexSet().index(vert) << std::endl;
+            break;
+          }
+        //x[vertex] = phi;
 
+        //std::cout << "pos: " << pos << ",  phi: " << x[vertex] << std::endl;
         vertex++;
       }
 
@@ -380,7 +399,10 @@ int main (int argc, char *argv[]) try
 
   // Output initial iterate (of homotopy loop)
   SolutionType identity;
-  Dune::Functions::interpolate(powerBasis, identity, [](FieldVector<double,dim> x){ return x; });
+  Dune::Functions::interpolate(powerBasis, identity, [&F0](FieldVector<double,dim> x)
+  {
+    return FieldVector<double,2>{ (F0[0][0]-1)*x[0] + F0[0][1]*x[1], F0[1][0]*x[0] + (F0[1][1]-1)*x[1] };
+  });
 
   SolutionType displacement = x;
   displacement -= identity;
@@ -388,8 +410,8 @@ int main (int argc, char *argv[]) try
   auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
   auto localDisplacementFunction = localFunction(displacementFunction);
 
-  //  We need to subsample, because VTK cannot natively display real second-order functions
-  SubsamplingVTKWriter<GridView> vtkWriter(gridView, Dune::refinementLevels(0));
+  //  Write the initial iterate
+  VTKWriter<GridView> vtkWriter(gridView);
   vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
   vtkWriter.write(resultPath + "quasiconvexity-test_homotopy_0");
 
@@ -420,73 +442,73 @@ int main (int argc, char *argv[]) try
 
     if (parameterSet.get<std::string>("energy") == "kpower")
     {
-      auto density = std::make_shared<KPowerDensity<dim,adouble> >(materialParameters);
+      auto density = std::make_shared<KPowerDensity<dim,adouble> >(F0, 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,
-                                                      adouble> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "nefflog")
-      elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
-                                                      FEBasis::LocalView::Tree::FiniteElement,
-                                                      adouble> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "nefflog2")
-      elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
-                                                      FEBasis::LocalView::Tree::FiniteElement,
-                                                      adouble> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "neffreduced")
-      elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
-                                                         FEBasis::LocalView::Tree::FiniteElement,
-                                                         adouble> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "cosh")
-      elasticEnergy = std::make_shared<CosHEnergy<GridView,
-                                                         FEBasis::LocalView::Tree::FiniteElement,
-                                                         adouble> >(materialParameters);
-
-    if (parameterSet.get<std::string>("energy") == "expacosh")
-      elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
-                                                         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") == "neff")
+//       elasticEnergy = std::make_shared<NeffEnergy<GridView,
+//                                                       FEBasis::LocalView::Tree::FiniteElement,
+//                                                       adouble> >(materialParameters);
+//
+//     if (parameterSet.get<std::string>("energy") == "nefflog")
+//       elasticEnergy = std::make_shared<NeffLogEnergy<GridView,
+//                                                       FEBasis::LocalView::Tree::FiniteElement,
+//                                                       adouble> >(materialParameters);
+//
+//     if (parameterSet.get<std::string>("energy") == "nefflog2")
+//       elasticEnergy = std::make_shared<NeffLogEnergy2<GridView,
+//                                                       FEBasis::LocalView::Tree::FiniteElement,
+//                                                       adouble> >(materialParameters);
+//
+//     if (parameterSet.get<std::string>("energy") == "neffreduced")
+//       elasticEnergy = std::make_shared<NeffReducedEnergy<GridView,
+//                                                          FEBasis::LocalView::Tree::FiniteElement,
+//                                                          adouble> >(materialParameters);
+//
+//     if (parameterSet.get<std::string>("energy") == "cosh")
+//       elasticEnergy = std::make_shared<CosHEnergy<GridView,
+//                                                          FEBasis::LocalView::Tree::FiniteElement,
+//                                                          adouble> >(materialParameters);
+//
+//     if (parameterSet.get<std::string>("energy") == "expacosh")
+//       elasticEnergy = std::make_shared<ExpACosHEnergy<GridView,
+//                                                          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")
     {
-      auto density = std::make_shared<MagicDensity<dim,adouble> >(materialParameters);
+      auto density = std::make_shared<MagicDensity<dim,adouble> >(F0,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") == "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);
+      auto density = std::make_shared<VossDensity<dim,adouble> >(F0,materialParameters);
       elasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<GridView,
                                                                        FEBasis::LocalView::Tree::FiniteElement,
                                                                        adouble> >(density);
@@ -654,7 +676,7 @@ int main (int argc, char *argv[]) try
         for (int j=0; j<dim; j++)
           derivative[j].axpy(localConfiguration[i][j], gradients[i]);
 
-    //auto det = derivative.determinant();
+      derivative += F0;
 
       //std::cout << derivative << std::endl;
       deformationDeterminant[localP0View.index(0)] = derivative.determinant();
@@ -741,8 +763,7 @@ int main (int argc, char *argv[]) try
   displacement -= identity;
   auto displacementFunction2 = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(feBasis, displacement);
 
-  //  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.clear();
   vtkWriter.addVertexData(displacementFunction2, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
   vtkWriter.write(resultPath + "final-result");
-- 
GitLab