The quasiconvexity test program with micromorphic relaxation

dune/elasticity/materials/micromorphicallyrelaxedenergy.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_ELASTICITY_MATERIALS_MICROMORPHICALLYRELAXEDENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_MICROMORPHICALLYRELAXEDENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localenergy.hh>
+#include <dune/elasticity/materials/localdensity.hh>
+
+namespace Dune::Elasticity {
+
+template<class LocalView, class field_type=double>
+class MicromorphicallyRelaxedEnergy
+  : public Elasticity::LocalEnergy<LocalView,field_type>
+{
+  using GridView = typename LocalView::GridView;
+  using DT = typename GridView::Grid::ctype;
+
+  enum {gridDim=GridView::dimension};
+
+public:
+
+  /** \brief Constructor with a local energy density
+    */
+  MicromorphicallyRelaxedEnergy(const std::shared_ptr<LocalDensity<gridDim,field_type,DT>>& ld,
+                                double L_c)
+  : localDensity_(ld),
+    L_c_(L_c)
+  {}
+
+  /** \brief Virtual destructor */
+  virtual ~MicromorphicallyRelaxedEnergy()
+  {}
+
+  /** \brief Assemble the energy for a single element */
+  field_type energy(const LocalView& localView,
+                    const std::vector<field_type>& localConfiguration) const override;
+
+  void setLc(double L_c)
+  {
+    L_c_ = L_c;
+  }
+
+protected:
+  const std::shared_ptr<LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
+
+  /** \brief Strength of the micromorphic regularization term */
+  double L_c_;
+};
+
+template <class LocalView, class field_type>
+field_type
+MicromorphicallyRelaxedEnergy<LocalView, field_type>::
+energy(const LocalView& localView,
+       const std::vector<field_type>& localConfiguration) const
+{
+  // powerBasis: grab the finite element of the first child
+  const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+  const auto& element = localView.element();
+
+  field_type energy = 0;
+
+  // store values and gradients of basis functions
+  using RangeType    = typename std::decay_t<decltype(localFiniteElement)>::Traits::LocalBasisType::Traits::RangeType;
+  using JacobianType = typename std::decay_t<decltype(localFiniteElement)>::Traits::LocalBasisType::Traits::JacobianType;
+  std::vector<RangeType> values(localFiniteElement.size());
+  std::vector<JacobianType> jacobians(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 geometryJacobianIT = element.geometry().jacobianInverseTransposed(qp.position());
+
+    // Global position
+    auto x = element.geometry().global(qp.position());
+
+    // Get values of the basis functions
+    localFiniteElement.localBasis().evaluateFunction(qp.position(), values);
+    localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
+    for (std::size_t i=0; i<jacobians.size(); i++)
+      jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
+
+    // Deformation gradient
+    // This class is supposed to work with discretizations that approximate the deformation gradient
+    // directly.  Therefore we are evaluation shape function values, but the result is still
+    // a deformation *gradient*.
+    FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
+    for (size_t i=0; i<values.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], values[i]);
+
+    // Integrate the energy density
+    energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+
+    // Compute curl
+
+    static_assert(gridDim==2, "curl is only implemented for dim==2");
+
+    field_type P12x = 0;
+    field_type P11y = 0;
+    field_type P22x = 0;
+    field_type P21y = 0;
+
+    for (std::size_t i=0; i<jacobians.size(); i++)
+    {
+      P12x += localConfiguration[localView.tree().child(0).localIndex(i)] * jacobians[i][1][0];
+      P11y += localConfiguration[localView.tree().child(0).localIndex(i)] * jacobians[i][0][1];
+      P22x += localConfiguration[localView.tree().child(1).localIndex(i)] * jacobians[i][1][0];
+      P21y += localConfiguration[localView.tree().child(1).localIndex(i)] * jacobians[i][0][1];
+    }
+
+    FieldVector<field_type,2> curl;
+    curl[0] = P12x - P11y;
+    curl[1] = P22x - P21y;
+
+
+    // Integrate the micromorphic regularization
+    energy += qp.weight() * integrationElement * 0.5 * power(L_c_, 2) * curl.two_norm2();
+  }
+
+  return energy;
+}
+
+}  // namespace Dune::Elasticity
+
+#endif   //#ifndef DUNE_ELASTICITY_MATERIALS_MICROMORPHICALLYRELAXEDENERGY_HH

src/quasiconvexity-test-micromorphic.cc

+#include <config.h>
+
+// Includes for the ADOL-C automatic differentiation library
+// Need to come before (almost) all others.
+#include <adolc/adouble.h>
+#include <adolc/drivers/drivers.h>    // use of "Easy to Use" drivers
+#include <adolc/taping.h>
+
+#include <dune/fufem/utilities/adolcnamespaceinjections.hh>
+
+#undef HAVE_DUNE_PARMG
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/transpose.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/grid/io/file/gmshreader.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/functionspacebases/nedelecbasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+#include <dune/fufem/dunepython.hh>
+
+#include <dune/solvers/iterationsteps/mmgstep.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/norms/energynorm.hh>
+
+#include <dune/elasticity/common/trustregionsolver.hh>
+#include <dune/elasticity/assemblers/localadolcstiffness.hh>
+#include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/materials/micromorphicallyrelaxedenergy.hh>
+#include <dune/elasticity/materials/quasiconvexitytestdensities.hh>
+
+// grid dimension
+const int dim = 2;
+
+const int order = 1;
+
+using namespace Dune;
+
+template <class Basis, class Vector>
+void writeDisplacement(const Basis& basis, const Vector& deformation, const FieldMatrix<double,2,2>& F0, std::string filename)
+{
+  // Compute the determinant per element, for better understanding of the solution
+  Functions::LagrangeBasis<typename Basis::GridView,0> p0Basis(basis.gridView());
+  std::vector<double> deformationDeterminant(p0Basis.size());
+  // K = 0.5 F.frobenius_norm2() / F.determinant();
+  std::vector<double> K(p0Basis.size());
+
+  auto localView = basis.localView();
+  auto localP0View = p0Basis.localView();
+
+  for (const auto& element : elements(basis.gridView()))
+  {
+    localView.bind(element);
+    localP0View.bind(element);
+
+    const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+    std::vector<double> localConfiguration(localView.size());
+
+    auto deformationBackend = Dune::Functions::istlVectorBackend(deformation);
+    for (size_t i=0; i<localConfiguration.size(); i++)
+      localConfiguration[i] = deformationBackend[localView.index(i)];
+
+    // store values of shape functions
+    std::vector<Dune::FieldVector<double,dim> > shapeFunctionValues(localFiniteElement.size());
+
+    auto p = element.geometry().center();
+
+    // get gradients of shape functions
+    localFiniteElement.localBasis().evaluateFunction(p, shapeFunctionValues);
+
+    FieldMatrix<double,dim,dim> deformationGradient(0);
+    for (size_t i=0; i<shapeFunctionValues.size(); i++)
+      for (int j=0; j<dim; j++)
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ],
+                                    shapeFunctionValues[i]);
+
+    deformationGradient += F0;
+
+    //std::cout << derivative << std::endl;
+    deformationDeterminant[localP0View.index(0)] = deformationGradient.determinant();
+    K[localP0View.index(0)] = 0.5 * deformationGradient.frobenius_norm2() / deformationGradient.determinant();
+  }
+
+  std::cout << "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);
+
+  // Actually write the functions
+  VTKWriter<typename Basis::GridView> vtkWriter(basis.gridView());
+  vtkWriter.addCellData(deformationDeterminantFunction,
+                        VTK::FieldInfo("determinant", VTK::FieldInfo::Type::scalar, 1));
+  vtkWriter.addCellData(KFunction,
+                        VTK::FieldInfo("K", VTK::FieldInfo::Type::scalar, 1));
+  vtkWriter.write(filename);
+}
+
+int main (int argc, char *argv[]) try
+{
+  // initialize MPI, finalize is done automatically on exit
+  Dune::MPIHelper& mpiHelper = MPIHelper::instance(argc, argv);
+
+  // Start Python interpreter
+  Python::start();
+  Python::Reference main = Python::import("__main__");
+  Python::run("import math");
+
+  //feenableexcept(FE_INVALID);
+  Python::runStream()
+        << std::endl << "import sys"
+        << std::endl << "import os"
+        << std::endl << "sys.path.append(os.getcwd() + '/../../problems/')"
+        << std::endl;
+
+  using Vector = BlockVector<FieldVector<double,dim> >;
+
+  // parse data file
+  ParameterTree parameterSet;
+  if (argc < 2)
+    DUNE_THROW(Exception, "Usage: ./quasiconvexity-test-micromorphic <parameter file>");
+
+  ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+  ParameterTreeParser::readOptions(argc, argv, parameterSet);
+
+  // read solver settings
+  const int numLevels                   = parameterSet.get<int>("numLevels");
+  const double tolerance                = parameterSet.get<double>("tolerance");
+  const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
+  const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+  const int multigridIterations         = parameterSet.get<int>("numIt");
+  const int nu1                         = parameterSet.get<int>("nu1");
+  const int nu2                         = parameterSet.get<int>("nu2");
+  const int mu                          = parameterSet.get<int>("mu");
+  const int baseIterations              = parameterSet.get<int>("baseIt");
+  const double mgTolerance              = parameterSet.get<double>("mgTolerance");
+  const double baseTolerance            = parameterSet.get<double>("baseTolerance");
+  std::string resultPath                = parameterSet.get("resultPath", "");
+
+  // ///////////////////////////////////////
+  //    Create the grid
+  // ///////////////////////////////////////
+  using Grid = UGGrid<dim>;
+
+  std::shared_ptr<Grid> grid;
+
+  FieldVector<double,dim> lower(0), upper(1);
+
+  if (parameterSet.get<bool>("structuredGrid")) {
+
+    lower = parameterSet.get<FieldVector<double,dim> >("lower");
+    upper = parameterSet.get<FieldVector<double,dim> >("upper");
+
+    std::array<unsigned int,dim> elements = parameterSet.get<std::array<unsigned int,dim> >("elements");
+    grid = StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, elements);
+
+  } else {
+    std::string path                = parameterSet.get<std::string>("path");
+    std::string gridFile            = parameterSet.get<std::string>("gridFile");
+    grid = std::shared_ptr<Grid>(GmshReader<Grid>::read(path + "/" + gridFile));
+  }
+
+  grid->globalRefine(numLevels-1);
+
+#if HAVE_DUNE_PARMG
+  typedef Grid::LevelGridView GridView;
+  GridView gridView = grid->levelGridView(0);
+#else
+  typedef Grid::LeafGridView GridView;
+  GridView gridView = grid->leafGridView();
+#endif
+
+  // Extract all boundary vertices
+//   std::vector<std::pair<std::size_t, FieldVector<double,dim> > > boundaryVertices;
+//   BoundaryPatch<GridView> domainBoundary(gridView, true);
+
+  using namespace Dune::Functions::BasisFactory;
+
+  // FE basis spanning the FE space that we are working in
+  auto feBasis = makeBasis(
+      gridView,
+      power<dim>(
+        nedelec<1,1>(),
+        blockedInterleaved()
+    ));
+  using FEBasis = decltype(feBasis);
+
+  std::cout << "Basis has " << feBasis.size() << " degrees of freedom" << std::endl;
+
+  ///////////////////////////////////////////
+  //   Set Dirichlet values
+  ///////////////////////////////////////////
+
+  // The entire boundary is Dirichlet boundary
+  BitSetVector<1> dirichletVertices(gridView.size(dim), true);
+  BoundaryPatch<GridView> dirichletBoundary(gridView, dirichletVertices);
+
+  // TODO: Do we still need to distinguish between dirichletNodes and dirichletDofs?
+  BitSetVector<dim> dirichletDofs(feBasis.size(), false);
+  constructBoundaryDofs(dirichletBoundary,feBasis,dirichletDofs);
+
+  ////////////////////////////
+  //   Initial iterate
+  ////////////////////////////
+
+  Vector x(feBasis.size());
+
+  // Set the underlying affine deformation
+  auto x0 = parameterSet.get<double>("x0");
+
+  FieldMatrix<double,2,2> F0 = { {sqrt(x0), 0},
+                                 {       0, 1.0 / sqrt(x0)} };
+
+  // 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"))
+  {
+    DUNE_THROW(NotImplemented, "initialIteratePython is not implemented");
+  }
+
+  if (parameterSet.hasKey("initialIterateFile"))
+  {
+    DUNE_THROW(NotImplemented, "initialIterateFile is not implemented");
+  }
+
+  // Output initial iterate
+  writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-initial");
+
+  //////////////////////////////////////////////////////////////
+  //   Create an assembler for the energy functional
+  //////////////////////////////////////////////////////////////
+
+  const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+  const double L_c = parameterSet.get<double>("L_c");   // Micromorphic regularization strength
+
+  if (mpiHelper.rank() == 0)
+  {
+    std::cout << "Material parameters:" << std::endl;
+    materialParameters.report();
+  }
+
+  // Assembler using ADOL-C
+  std::cout << "Selected energy is: " << parameterSet.get<std::string>("energy") << std::endl;
+//   std::shared_ptr<Elasticity::LocalEnergy<FEBasis::LocalView,adouble> > elasticEnergy;
+  std::shared_ptr<Elasticity::MicromorphicallyRelaxedEnergy<FEBasis::LocalView,adouble> > elasticEnergy;
+
+  if (parameterSet.get<std::string>("energy") == "kpower")
+  {
+    auto density = std::make_shared<Elasticity::KPowerDensity<dim,adouble> >(F0, materialParameters);
+    elasticEnergy = std::make_shared<Elasticity::MicromorphicallyRelaxedEnergy<FEBasis::LocalView,
+                                                                     adouble> >(density,L_c);
+  }
+
+  if (parameterSet.get<std::string>("energy") == "expacosh")
+  {
+    auto density = std::make_shared<Elasticity::ExpACosHDensity<dim,adouble> >(F0,materialParameters);
+    elasticEnergy = std::make_shared<Elasticity::MicromorphicallyRelaxedEnergy<FEBasis::LocalView,
+                                                                     adouble> >(density,L_c);
+  }
+
+  if (parameterSet.get<std::string>("energy") == "magic")
+  {
+    auto density = std::make_shared<Elasticity::MagicDensity<dim,adouble> >(F0,materialParameters);
+    elasticEnergy = std::make_shared<Elasticity::MicromorphicallyRelaxedEnergy<FEBasis::LocalView,
+                                                                     adouble> >(density,L_c);
+  }
+
+  if (parameterSet.get<std::string>("energy") == "burkholder")
+  {
+    auto density = std::make_shared<Elasticity::BurkholderDensity<dim,adouble> >(F0,materialParameters);
+    elasticEnergy = std::make_shared<Elasticity::MicromorphicallyRelaxedEnergy<FEBasis::LocalView,
+                                                                     adouble> >(density,L_c);
+  }
+
+  if (parameterSet.get<std::string>("energy") == "voss")
+  {
+    auto density = std::make_shared<Elasticity::VossDensity<dim,adouble> >(F0,materialParameters);
+    elasticEnergy = std::make_shared<Elasticity::MicromorphicallyRelaxedEnergy<FEBasis::LocalView,
+                                                                     adouble> >(density,L_c);
+  }
+
+  if(!elasticEnergy)
+    DUNE_THROW(Exception, "Error: Selected energy not available!");
+
+  auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<FEBasis::LocalView> >(elasticEnergy);
+
+  Elasticity::FEAssembler<FEBasis,Vector> assembler(feBasis, localADOLCStiffness);
+
+  ///////////////////////////////////////////////////
+  //   Create a Riemannian trust-region solver
+  ///////////////////////////////////////////////////
+
+  // Multigrid transfer operators
+  using TransferOperatorType = typename TruncatedCompressedMGTransfer<Vector>::TransferOperatorType;
+  std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<Vector> > > transferOperators(numLevels-1);
+
+  //using Matrix = BCRSMatrix<FieldMatrix<double,1,1> >;
+  using Matrix = TransferOperatorType;
+
+  for (int i=0; i<grid->maxLevel(); i++)
+  {
+    // Coarse and fine level basis
+    auto coarseLevelBasis = makeBasis(
+      grid->levelGridView(i),
+      power<dim>(
+        nedelec<1,1>(),
+        blockedInterleaved()
+    ));
+
+    auto fineLevelBasis = makeBasis(
+      grid->levelGridView(i+1),
+      power<dim>(
+        nedelec<1,1>(),
+        blockedInterleaved()
+    ));
+
+    // Assemble the transfer operator matrix
+    auto transferMatrix = std::make_shared<Matrix>();
+    assembleGlobalBasisTransferMatrix(*transferMatrix,coarseLevelBasis,fineLevelBasis);
+
+    transferOperators[i] = std::make_shared<TruncatedCompressedMGTransfer<Vector> >();
+    transferOperators[i]->setMatrix(transferMatrix);
+  }
+
+  // The actual solver
+  TrustRegionSolver<FEBasis,Vector> solver;
+  solver.setup(*grid,
+               &assembler,
+               x,
+               dirichletDofs,
+               tolerance,
+               maxTrustRegionSteps,
+               initialTrustRegionRadius,
+               multigridIterations,
+               mgTolerance,
+               mu, nu1, nu2,
+               baseIterations,
+               baseTolerance,
+               1.0
+              );
+
+   solver.transferOperators_ = transferOperators;
+
+  ///////////////////////////////////////////////////////
+  //   Solve!
+  ///////////////////////////////////////////////////////
+
+  solver.iterateNamePrefix_ = "stage1-";
+  solver.setInitialIterate(x);
+  solver.solve();
+
+  x = solver.getSol();
+
+  /////////////////////////////////
+  //   Output result
+  /////////////////////////////////
+
+  elasticEnergy->setLc(0);
+  std::cout << "Energy without regularization: " << assembler.computeEnergy(x) << std::endl;
+  elasticEnergy->setLc(parameterSet.get<double>("L_c"));
+  std::cout << "Energy with regularization: " << assembler.computeEnergy(x) << std::endl;
+
+  writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-result");
+
+} catch (Exception& e) {
+    std::cout << e.what() << std::endl;
+}