diff --git a/src/quasiconvexity-test-micromorphic.cc b/src/quasiconvexity-test-micromorphic.cc index 8afc8bd327422acade11faad3ad11fe67fefac82..10a9b1f38f77528d7f2f352b17813395083a91e2 100644 --- a/src/quasiconvexity-test-micromorphic.cc +++ b/src/quasiconvexity-test-micromorphic.cc @@ -110,156 +110,6 @@ void writeDisplacement(const Basis& basis, const Vector& deformation, const Fiel } - -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