WIP

dune/elasticity/common/trustregionsolver.cc

@@ -12,13 +12,11 @@
 #include <dune/fufem/assemblers/dunefunctionsoperatorassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/laplaceassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/massassembler.hh>
-#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
 #include <dune/fufem/assemblers/istlbackend.hh>
 
 // Using a monotone multigrid as the inner solver
 #include <dune/solvers/iterationsteps/trustregiongsstep.hh>
 #include <dune/solvers/iterationsteps/mmgstep.hh>
-#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
 #include <dune/solvers/transferoperators/mandelobsrestrictor.hh>
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/norms/twonorm.hh>
@@ -85,10 +83,11 @@ setup(const typename BasisType::GridView::Grid& grid,
     if (grid.comm().rank()==0)
         std::cout << "Parallel multigrid setup took " << setupTimer.stop() << " seconds." << std::endl;
 #else
-    std::conditional_t< // do we have a dune-functions basis?
-      Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
-        BasisType,
-        DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
+//     std::conditional_t< // do we have a dune-functions basis?
+//       Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+//         BasisType,
+//         DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
+    const auto& basis = assembler_->basis_;
     // ////////////////////////////////
     //   Create a multigrid solver
     // ////////////////////////////////
@@ -258,7 +257,7 @@ setup(const typename BasisType::GridView::Grid& grid,
     // ////////////////////////////////////
     //   Create the transfer operators
     // ////////////////////////////////////
-
+#if 0
     ////////////////////////////////////////////////////////////////////////
     //  The P1 space (actually P1/Q1, depending on the grid) is special:
     //  If we work in such a space, then the multigrid hierarchy of spaces
@@ -333,6 +332,9 @@ setup(const typename BasisType::GridView::Grid& grid,
     }
 
     mmgStep->setTransferOperators(transferOperators);
+#else
+    mmgStep->setTransferOperators(transferOperators_);
+#endif
 
     // //////////////////////////////////////////////////////////
     //   Create obstacles

dune/elasticity/common/trustregionsolver.hh

@@ -12,9 +12,11 @@
 #include <dune/istl/bvector.hh>
 
 #include <dune/solvers/common/boxconstraint.hh>
+#include <dune/solvers/iterationsteps/mmgstep.hh>
 #include <dune/solvers/norms/h1seminorm.hh>
 #include <dune/solvers/solvers/iterativesolver.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/transferoperators/truncatedcompressedmgtransfer.hh>
 
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 #include <dune/fufem/functionspacebases/p2nodalbasis.hh>
@@ -182,7 +184,10 @@ protected:
 #endif
     /** \brief The solver for the quadratic inner problems */
     std::shared_ptr<Solver> innerSolver_;
+public:
+    std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<VectorType> > > transferOperators_;
 
+protected:
 #if ! HAVE_DUNE_PARMG
     /** \brief Contains 'true' everywhere -- the trust-region is bounded */
     Dune::BitSetVector<blocksize> hasObstacle_;

problems/quasiconvexity-circle-expacosh.parset

@@ -12,7 +12,7 @@ gridFile = circle2ndorder.msh
 #elements = 2 2
 
 # Number of grid levels
-numLevels = 8
+numLevels = 6
 
 #############################################
 #  Solver parameters
@@ -65,6 +65,9 @@ c2 = 0
 
 []
 
+# Micromorphic regularization
+L_c = 10
+
 #############################################
 #  Boundary values
 #############################################

src/CMakeLists.txt

 if(ADOLC_FOUND AND IPOPT_FOUND AND PYTHONLIBS_FOUND AND dune-uggrid_FOUND)
   set(programs finite-strain-elasticity
                linear-elasticity
-               quasiconvexity-test)
+               quasiconvexity-test
+               quasiconvexity-test-micromorphic)
 
   foreach(_program ${programs})
     add_executable(${_program} ${_program}.cc)
@@ -11,8 +12,8 @@ if(ADOLC_FOUND AND IPOPT_FOUND AND PYTHONLIBS_FOUND AND dune-uggrid_FOUND)
     add_dune_ug_flags(${_program})
     add_dune_mpi_flags(${_program})
     add_dune_superlu_flags(${_program})
-#    target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
-    target_compile_options(${_program} PRIVATE "-g")
+    target_compile_options(${_program} PRIVATE "-fpermissive" -O3 -funroll-loops -DNDEBUG)
+#   target_compile_options(${_program} PRIVATE "-g")
   endforeach()
 
   if (AMIRAMESH_FOUND)

src/quasiconvexity-test-micromorphic.cc

@@ -109,6 +109,157 @@ void writeDisplacement(const Basis& basis, const Vector& deformation, const Fiel
   vtkWriter.write(filename);
 }
 
+
+
+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 setRegularization(double regularization)
+  {
+    regularization_ = regularization;
+  }
+
+protected:
+  const std::shared_ptr<LocalDensity<gridDim,field_type,DT>> localDensity_ = nullptr;
+
+  /** \brief Strength of the micromorphic regularization term */
+  double regularization_;
+};
+
+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;
+
+#if 0
+  // store gradients of shape functions and base functions
+  std::vector<FieldMatrix<DT,1,gridDim> > referenceGradients(localFiniteElement.size());
+  std::vector<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());
+
+    // Global position
+    auto x = element.geometry().global(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]);
+
+    // Deformation gradient
+    FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
+    for (size_t i=0; i<gradients.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], gradients[i]);
+
+    // Prescribed deformation gradient
+    // See the mail by Jendrik Voss from 13.1.2020
+    // TODO: This constructs the value of P from a basis and a vector of coefficients.
+    // It would be easier to pass a discrete functions.
+    FieldMatrix<field_type,gridDim,gridDim> P(0);
+    for (size_t i=0; i<values.size(); i++)
+      for (int j=0; j<gridDim; j++)
+        P[j].axpy(PCoefficients_[ localView.tree().child(j).index(i) ], PValues[i]);
+
+    // Integrate the energy density
+    energy += qp.weight() * integrationElement * (*localDensity_)(x, deformationGradient);
+
+    //
+    energy += q.weight() * integrationElement * 0.5 * k * (deformationGradient - P).frobenius_norm2();
+  }
+
+  return energy;
+#else
+  // 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
+
+    // Integrate the micromorphic regularization
+    energy += qp.weight() * integrationElement * 0.5 * power(L_c_, 2) * curl.two_norm2();
+  }
+
+  return energy;
+#endif
+}
+
+
 int main (int argc, char *argv[]) try
 {
   // initialize MPI, finalize is done automatically on exit

src/quasiconvexity-test.cc

@@ -27,6 +27,7 @@
 #include <dune/functions/functionspacebases/periodicbasis.hh>
 #include <dune/functions/functionspacebases/powerbasis.hh>
 
+#include <dune/fufem/assemblers/basisinterpolationmatrixassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/functiontools/boundarydofs.hh>
 #include <dune/fufem/dunepython.hh>
@@ -534,7 +535,7 @@ int main (int argc, char *argv[]) try
     auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,SolutionType> >(
         baseTolerance,
         baseIterations,
-        NumProc::REDUCED,
+        NumProc::QUIET,
         "mumps");
 #else
     // First create a Gauss-seidel base solver