From d712478aacf93397e93d5691577280ab32fd0ef1 Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 9 Sep 2021 10:25:46 +0200
Subject: [PATCH] Construct compatible deformation gradient from incompatible
one
By an L^2-projection onto the space of deformation gradients.
---
dune/elasticity/materials/pdistanceenergy.hh | 118 ++++++++
src/quasiconvexity-test-micromorphic.cc | 303 +++++++++++++++++++
2 files changed, 421 insertions(+)
create mode 100644 dune/elasticity/materials/pdistanceenergy.hh
diff --git a/dune/elasticity/materials/pdistanceenergy.hh b/dune/elasticity/materials/pdistanceenergy.hh
new file mode 100644
index 0000000..3a0ac65
--- /dev/null
+++ b/dune/elasticity/materials/pdistanceenergy.hh
@@ -0,0 +1,118 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_ELASTICITY_MATERIALS_PDISTANCEENERGY_HH
+#define DUNE_ELASTICITY_MATERIALS_PDISTANCEENERGY_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/transpose.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/elasticity/assemblers/localenergy.hh>
+
+namespace Dune::Elasticity {
+
+template<class LocalView, class DataFunction, class field_type=double>
+class PDistanceEnergy
+ : 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 data term
+ */
+ PDistanceEnergy(const DataFunction& dataFunction,
+ FieldMatrix<double,2,2> F0,
+ double k)
+ : dataFunction_(dataFunction),
+ F0_(F0),
+ k_(k)
+ {}
+
+ /** \brief Virtual destructor */
+ virtual ~PDistanceEnergy()
+ {}
+
+ /** \brief Assemble the energy for a single element */
+ field_type energy(const LocalView& localView,
+ const std::vector<field_type>& localConfiguration) const override;
+
+ /** \brief Set scaling factor */
+ void setK(double k)
+ {
+ k_ = k;
+ }
+
+protected:
+ const DataFunction dataFunction_;
+
+ FieldMatrix<double,gridDim,gridDim> F0_;
+
+ /** \brief Scaling factor */
+ double k_;
+};
+
+template <class LocalView, class DataFunction, class field_type>
+field_type
+PDistanceEnergy<LocalView, DataFunction, 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();
+
+ auto localDataFunction = localFunction(dataFunction_);
+ localDataFunction.bind(element);
+
+ field_type energy = 0;
+
+ // store gradients of shape functions and base functions
+ using FiniteElement = typename LocalView::Tree::ChildType::FiniteElement;
+ using Jacobian = typename FiniteElement::Traits::LocalBasisType::Traits::JacobianType;
+
+ std::vector<Jacobian> jacobians(localFiniteElement.size());
+
+ int quadOrder = (element.type().isSimplex()) ? localFiniteElement.localBasis().order()
+ : localFiniteElement.localBasis().order() * gridDim;
+
+ const auto& quad = 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());
+
+ // Get gradients of shape functions
+ localFiniteElement.localBasis().evaluateJacobian(qp.position(), jacobians);
+
+ // compute gradients of base functions
+ for (size_t i=0; i<jacobians.size(); ++i)
+ jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
+
+ // Deformation gradient
+ FieldMatrix<field_type,gridDim,gridDim> deformationGradient(0);
+ for (size_t i=0; i<jacobians.size(); i++)
+ for (int j=0; j<gridDim; j++)
+ deformationGradient[j].axpy(localConfiguration[ localView.tree().child(j).localIndex(i) ], jacobians[i][0]);
+
+ deformationGradient += F0_;
+
+ FieldMatrix<double,2,2> P = localDataFunction(qp.position());
+
+ auto diff = deformationGradient - P;
+
+ // Integrate the micromorphic regularization
+ energy += qp.weight() * integrationElement * k_ * diff.frobenius_norm2();
+ }
+
+ return energy;
+}
+
+} // namespace Dune::Elasticity
+
+#endif //#ifndef DUNE_ELASTICITY_MATERIALS_PDISTANCEENERGY_HH
diff --git a/src/quasiconvexity-test-micromorphic.cc b/src/quasiconvexity-test-micromorphic.cc
index 98dbaf1..f4df370 100644
--- a/src/quasiconvexity-test-micromorphic.cc
+++ b/src/quasiconvexity-test-micromorphic.cc
@@ -38,6 +38,7 @@
#include <dune/elasticity/assemblers/localadolcstiffness.hh>
#include <dune/elasticity/assemblers/feassembler.hh>
#include <dune/elasticity/materials/micromorphicallyrelaxedenergy.hh>
+#include <dune/elasticity/materials/pdistanceenergy.hh>
#include <dune/elasticity/materials/quasiconvexitytestdensities.hh>
// grid dimension
@@ -108,6 +109,302 @@ void writeDeformationGradientInfo(const Basis& basis, const Vector& deformation,
vtkWriter.write(filename);
}
+template <class Basis, class Vector>
+void writeDisplacement(const Basis& basis, const Vector& x, const FieldMatrix<double,2,2>& F0, std::string filename)
+{
+ using GridView = typename Basis::GridView;
+
+ // Add the underlying affine displacement.
+ Vector affineDisplacement;
+ Dune::Functions::interpolate(basis, affineDisplacement, [&F0](FieldVector<double,dim> pos)
+ {
+ return FieldVector<double,2>{ (F0[0][0]-1)*pos[0] + F0[0][1]*pos[1], F0[1][0]*pos[0] + (F0[1][1]-1)*pos[1] };
+ });
+
+ Vector displacement = x;
+ displacement += affineDisplacement;
+
+ // Compute the determinant per element, for better understanding of the solution
+ Functions::LagrangeBasis<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();
+
+ Vector localConfiguration(localView.size());
+
+ for (size_t i=0; i<localConfiguration.size(); i++)
+ localConfiguration[i] = x[localView.index(i)[0]];
+
+ // store gradients of shape functions and base functions
+ std::vector<Dune::FieldMatrix<double,1,dim> > jacobians(localFiniteElement.size());
+
+ auto p = referenceElement(element).position(0,0); // Element center in local coordinates
+
+ const auto geometryJacobianIT = element.geometry().jacobianInverseTransposed(p);
+
+ // get gradients of shape functions
+ localFiniteElement.localBasis().evaluateJacobian(p, jacobians);
+
+ // compute Jacobians of basis functions
+ for (size_t i=0; i<jacobians.size(); ++i)
+ jacobians[i] = jacobians[i] * transpose(geometryJacobianIT);
+
+ FieldMatrix<double,dim,dim> derivative(0);
+ for (size_t i=0; i<jacobians.size(); i++)
+ for (int j=0; j<dim; j++)
+ derivative[j].axpy(localConfiguration[i][j], jacobians[i][0]);
+
+ derivative += F0;
+
+ deformationDeterminant[localP0View.index(0)] = derivative.determinant();
+ K[localP0View.index(0)] = 0.5 * derivative.frobenius_norm2() / derivative.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 displacement function
+ auto displacementFunction = Dune::Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis, displacement);
+
+ VTKWriter<GridView> vtkWriter(basis.gridView());
+ vtkWriter.addVertexData(displacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
+ 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);
+}
+
+
+template <typename Grid, typename IncompatibleDeformationGradientFunction>
+void constructDeformation(const std::shared_ptr<Grid>& grid,
+ const ParameterTree& parameterSet,
+ const IncompatibleDeformationGradientFunction& incompatibleDeformationGradient)
+{
+ using Vector = BlockVector<FieldVector<double,dim> >;
+
+ // 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
+ /////////////////////////////////////////
+
+#if HAVE_DUNE_PARMG
+ using GridView = typename Grid::LevelGridView;
+ GridView gridView = grid->levelGridView(0);
+#else
+ using GridView = typename Grid::LeafGridView;
+ GridView gridView = grid->leafGridView();
+#endif
+
+ // FE basis spanning the FE space that we are working in
+ using namespace Dune::Functions::BasisFactory;
+ auto feBasis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+ using FEBasis = decltype(feBasis);
+
+ std::cout << "Basis has " << feBasis.size() << " degrees of freedom" << std::endl;
+
+ ///////////////////////////////////////////
+ // Set Dirichlet values
+ ///////////////////////////////////////////
+
+ BoundaryPatch<GridView> dirichletBoundary(gridView, true);
+
+ 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));
+
+ // Output initial iterate
+ writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-deformation-initial");
+
+ //////////////////////////////////////////////////////////////
+ // Create an assembler for the energy functional
+ //////////////////////////////////////////////////////////////
+
+ const ParameterTree& materialParameters = parameterSet.sub("materialParameters");
+
+ std::cout << "Material parameters:" << std::endl;
+ materialParameters.report();
+
+ // Assembler using ADOL-C
+ auto energy = std::make_shared<Elasticity::PDistanceEnergy<typename FEBasis::LocalView,
+ IncompatibleDeformationGradientFunction,
+ adouble> >(incompatibleDeformationGradient,
+ F0,
+ parameterSet.get<double>("k"));
+
+ auto localADOLCStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<typename FEBasis::LocalView> >(energy);
+
+ auto assembler = std::make_shared<Elasticity::FEAssembler<FEBasis,Vector> >(feBasis, localADOLCStiffness);
+
+ ///////////////////////////////////////////////////
+ // Create a trust-region solver
+ ///////////////////////////////////////////////////
+
+ using Matrix = BCRSMatrix<FieldMatrix<double, dim, dim> >;
+
+#ifdef HAVE_IPOPT
+ // First create an IPOpt base solver
+ auto baseSolver = std::make_shared<QuadraticIPOptSolver<Matrix,Vector> >(
+ baseTolerance,
+ baseIterations,
+ NumProc::QUIET,
+ "mumps");
+#else
+ // First create a Gauss-seidel base solver
+ TrustRegionGSStep<Matrix, Vector>* baseSolverStep = new TrustRegionGSStep<Matrix, Vector>;
+
+ // Hack: the two-norm may not scale all that well, but it is fast!
+ TwoNorm<CorrectionType>* baseNorm = new TwoNorm<Vector>;
+
+ auto baseSolver = std::make_shared<::LoopSolver<Vector>>(baseSolverStep,
+ baseIterations,
+ baseTolerance,
+ baseNorm,
+ Solver::QUIET);
+#endif
+ // Make pre and postsmoothers
+ auto presmoother = std::make_shared< TrustRegionGSStep<Matrix, Vector> >();
+ auto postsmoother = std::make_shared< TrustRegionGSStep<Matrix, Vector> >();
+
+ auto mmgStep = std::make_shared<MonotoneMGStep<Matrix, Vector> >();
+
+ mmgStep->setVerbosity(NumProc::FULL);
+ mmgStep->setMGType(mu, nu1, nu2);
+ mmgStep->ignoreNodes_ = &dirichletDofs;
+ mmgStep->setBaseSolver(baseSolver);
+ mmgStep->setSmoother(presmoother, postsmoother);
+ mmgStep->setObstacleRestrictor(std::make_shared<MandelObstacleRestrictor<Vector> >());
+
+ //////////////////////////////////////
+ // Create the transfer operators
+ //////////////////////////////////////
+
+ using TransferOperatorType = typename TruncatedCompressedMGTransfer<Vector>::TransferOperatorType;
+ std::vector<std::shared_ptr<TruncatedCompressedMGTransfer<Vector>>> transferOperators(numLevels-1);
+
+ // For the restriction operators: FE bases on all levels
+ auto dummyLevelBasis = makeBasis(
+ grid->levelGridView(0),
+ power<dim>(
+ lagrange<order>(),
+ blockedInterleaved()
+ ));
+
+ using LevelBasis = decltype(dummyLevelBasis);
+
+ // Construct the actual level bases
+ std::vector<std::shared_ptr<LevelBasis> > levelBases(numLevels);
+ for (int j=0; j<numLevels; j++)
+ {
+ levelBases[j] = std::make_shared<LevelBasis>(makeBasis(
+ grid->levelGridView(j),
+ power<dim>(
+ lagrange<order>(),
+ blockedInterleaved()
+ )));
+ }
+
+ for (int j=0; j<numLevels-1; j++)
+ {
+ auto transferMatrix = std::make_shared<TransferOperatorType>();
+ assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases[j], *levelBases[j+1]);
+
+ // Construct the local multigrid transfer matrix
+ transferOperators[j] = std::make_shared<TruncatedCompressedMGTransfer<Vector> >();
+ transferOperators[j]->setMatrix(transferMatrix);
+ }
+
+ // What does this block do?
+ {
+ auto transferMatrix = std::make_shared<TransferOperatorType>();
+ assembleGlobalBasisTransferMatrix(*transferMatrix, *levelBases.back(), *levelBases.back());
+
+ // Construct the local multigrid transfer matrix
+ transferOperators.push_back(std::make_shared<TruncatedCompressedMGTransfer<Vector> >());
+ transferOperators.back()->setMatrix(transferMatrix);
+ }
+
+ mmgStep->setTransferOperators(transferOperators);
+
+ TrustRegionSolver<FEBasis,Vector> solver;
+ solver.setup(*grid,
+ assembler,
+ x,
+ dirichletDofs,
+ tolerance,
+ maxTrustRegionSteps,
+ initialTrustRegionRadius,
+ mmgStep,
+ mgTolerance,
+ multigridIterations);
+
+
+ ///////////////////////////////////////////////////////
+ // Solve!
+ // /////////////////////////////////////////////////////
+
+ std::cout << "2-Distance to incompatible strain: " << std::sqrt(assembler->computeEnergy(x)) << std::endl;
+
+ solver.setInitialIterate(x);
+ solver.solve();
+
+ x = solver.getSol();
+
+ /////////////////////////////////
+ // Output result
+ /////////////////////////////////
+
+ std::cout << "2-Distance to incompatible strain: " << std::sqrt(assembler->computeEnergy(x)) << std::endl;
+
+ writeDisplacement(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-deformation");
+}
+
int main (int argc, char *argv[]) try
{
// initialize MPI, finalize is done automatically on exit
@@ -404,6 +701,12 @@ int main (int argc, char *argv[]) try
writeDeformationGradientInfo(feBasis, x, F0, resultPath + "quasiconvexity-test-micromorphic-result_" + std::to_string(L_c) + "_" + std::to_string(x0));
+ //////////////////////////////////////////////////////////////////
+ // Construct a deformation field that has a deformation gradient
+ // close to the one we just constructed.
+ //////////////////////////////////////////////////////////////////
+ auto incompatibleDeformationGradient = Functions::makeDiscreteGlobalBasisFunction<FieldMatrix<double,2,2>>(feBasis, x);
+ constructDeformation(grid, parameterSet, incompatibleDeformationGradient);
} catch (Exception& e) {
std::cout << e.what() << std::endl;
}
--
GitLab