diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc index ecae2a3d5dcab0a88ca478aad8b9cf687dc9651f..a7a91f9422631016ffe4a1d07dab5a5381858202 100644 --- a/dune/elasticity/common/trustregionsolver.cc +++ b/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 diff --git a/dune/elasticity/common/trustregionsolver.hh b/dune/elasticity/common/trustregionsolver.hh index 3947b55b86ae3b32f2660eca8ce42b2ac2506365..7a3a1d6e751eb7af7d4c33805901a2fcf59f0ee0 100644 --- a/dune/elasticity/common/trustregionsolver.hh +++ b/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_; diff --git a/problems/quasiconvexity-circle-expacosh.parset b/problems/quasiconvexity-circle-expacosh.parset index 247656d5ce1e15f3d4d4c7156ed6d11ac4f03bee..2f0032bcdf6eda257f44b2da0804c6398c415824 100644 --- a/problems/quasiconvexity-circle-expacosh.parset +++ b/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 ############################################# diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 0d66716c1918b15fec869eb02d345b8afe3e48d7..84b207b05fb14a1e29383e71bfd1fe533679fcde 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,7 +1,8 @@ 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) diff --git a/src/quasiconvexity-test-micromorphic.cc b/src/quasiconvexity-test-micromorphic.cc index e20b21fad394a6abcbb4d7dfa05512053520636e..8afc8bd327422acade11faad3ad11fe67fefac82 100644 --- a/src/quasiconvexity-test-micromorphic.cc +++ b/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 diff --git a/src/quasiconvexity-test.cc b/src/quasiconvexity-test.cc index 103950f43c8a71b9f52ef2259b85c694c49aa096..94905a9e788e55f873f26a62af04aa268643b0d8 100644 --- a/src/quasiconvexity-test.cc +++ b/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