solidmechanics/hyperelastic_voliso/localresidual.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_SOLIDMECHANICS_HYPERELASTIC_VOLISO_LOCAL_RESIDUAL_HH
#define DUMUX_SOLIDMECHANICS_HYPERELASTIC_VOLISO_LOCAL_RESIDUAL_HH
 
#include <dune/common/fmatrix.hh>
#include <dune/grid/common/intersectioniterator.hh>
 
#include <dumux/common/math.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/common/typetraits/localdofs_.hh>
#include <dumux/discretization/defaultlocaloperator.hh>
#include <dumux/discretization/evalsolution.hh>
#include <dumux/discretization/cvfe/quadraturerules.hh>
 
namespace Dumux {

template<class TypeTag>
class HyperelasticVolIsoMomentumLocalResidual : public DiscretizationDefaultLocalOperator<TypeTag>
{
    using ParentType = DiscretizationDefaultLocalOperator<TypeTag>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using ElementDiscretization = typename GridGeometry::LocalView;
    using GridView = typename GridGeometry::GridView;
    using GlobalPosition = typename GridView::template Codim<0>::Entity::Geometry::GlobalCoordinate;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using Tensor = Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension>;
    static constexpr int dim = GridView::dimension;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using ElementVariables = typename GridVariables::GridVariablesCache::LocalView;
 
public:
    using ParentType::ParentType;
    using ElementResidualVector = typename ParentType::ElementResidualVector;
 
    NumEqVector computeStorage(const auto&, const auto&, const auto&) const
    { return NumEqVector(0.0); }
 
    NumEqVector storageIntegral(const ElementDiscretization&, const ElementVariables&,
                                const auto&, bool) const
    { return NumEqVector(0.0); }

    void addToElementFluxAndSourceResidual(ElementResidualVector& residual,
                                           const Problem& problem,
                                           const auto& element,
                                           const ElementDiscretization& elemDisc,
                                           const ElementVariables& elemVars) const
    {
        if constexpr (Dumux::Detail::LocalDofs::hasNonCVLocalDofsInterface<ElementDiscretization>())
        {
            if (nonCVLocalDofs(elemDisc).empty())
                return;
 
            for (const auto& qpData : CVFE::quadratureRule(elemDisc, element))
            {
                const auto& ipCache = cache(elemVars, qpData.ipData());
                const auto& ipGlobal = qpData.ipData().global();
 
                const auto F = deformationGradient_(elemDisc, elemVars, ipCache);
                const auto P = firstPiolaKirchhoffStressTensor_(problem, F, elemDisc, ipGlobal);
                const auto f = problem.source(elemDisc, elemVars, qpData.ipData());
                const auto& shapeValues = ipCache.shapeValues();
 
                for (const auto& nonCVdof : nonCVLocalDofs(elemDisc))
                {
                    const auto idx = nonCVdof.index();
                    for (int i = 0; i < dim; ++i)
                    {
                        for (int j = 0; j < dim; ++j)
                            residual[idx][Indices::momentum(i)] += P[i][j] * ipCache.gradN(idx)[j] * qpData.weight();
                        residual[idx][Indices::momentum(i)] -= f[Indices::momentum(i)] * shapeValues[idx][0] * qpData.weight();
                    }
                }
            }
        }
    }

    template<class SubControlVolumeFace>
    NumEqVector fluxIntegral(const ElementDiscretization& elemDisc,
                             const ElementVariables& elemVars,
                             const SubControlVolumeFace& scvf) const
    {
        const auto& problem = this->asImp().problem();
        const auto& n = scvf.unitOuterNormal();
        NumEqVector f(0.0);
 
        for (const auto& qpData : CVFE::quadratureRule(elemDisc, scvf))
        {
            const auto& ipCache = cache(elemVars, qpData.ipData());
            const auto F = deformationGradient_(elemDisc, elemVars, ipCache);
            const auto P = firstPiolaKirchhoffStressTensor_(problem, F, elemDisc, qpData.ipData().global());
 
            for (int i = 0; i < dim; ++i)
                for (int j = 0; j < dim; ++j)
                    f[Indices::momentum(i)] -= P[i][j] * n[j] * qpData.weight();
        }
        return f;
    }

    Scalar volumetricStressPressure(const Problem& problem,
                                    const ElementDiscretization& elemDisc,
                                    const GlobalPosition& ipGlobal) const
    { return problem.couplingManager().pressureAtPoint(elemDisc, ipGlobal); }
 
private:
    template<class IpCache>
    Tensor deformationGradient_(const ElementDiscretization& elemDisc,
                                const ElementVariables& elemVars,
                                const IpCache& ipCache) const
    {
        Tensor F(0.0);
        for (const auto& localDof : localDofs(elemDisc))
            for (int d = 0; d < dim; ++d)
                F[d].axpy(elemVars[localDof].displacement(d), ipCache.gradN(localDof.index()));
        for (int d = 0; d < dim; ++d)
            F[d][d] += 1.0;
        return F;
    }
 
    Tensor firstPiolaKirchhoffStressTensor_(const Problem& problem,
                                            const Tensor& F,
                                            const ElementDiscretization& elemDisc,
                                            const GlobalPosition& ipGlobal) const
    {
        auto P = problem.isochoricFirstPiolaKirchhoffStressTensor(F, elemDisc, ipGlobal);
 
        Tensor Finv(F);
        Finv.invert();
        const Scalar pJ = this->asImp().volumetricStressPressure(problem, elemDisc, ipGlobal)*F.determinant();
        for (int i = 0; i < dim; ++i)
            for (int j = 0; j < dim; ++j)
                P[i][j] += pJ*Finv[j][i];
 
        return P;
    }
};

template<class TypeTag>
class HyperelasticVolIsoPressureLocalResidual : public DiscretizationDefaultLocalOperator<TypeTag>
{
    using ParentType = DiscretizationDefaultLocalOperator<TypeTag>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using NumEqVector = Dumux::NumEqVector<GetPropType<TypeTag, Properties::PrimaryVariables>>;
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using ElementDiscretization = typename GridGeometry::LocalView;
    using GridView = typename GridGeometry::GridView;
    using Tensor = Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension>;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using ElementVariables = typename GridVariables::GridVariablesCache::LocalView;
 
public:
    using ParentType::ParentType;
    using ElementResidualVector = typename ParentType::ElementResidualVector;
 
    NumEqVector computeStorage(const Problem&, const auto&, const auto&) const
    { return NumEqVector(0.0); }
 
    NumEqVector storageIntegral(const ElementDiscretization&, const ElementVariables&,
                                const auto&, bool) const
    { return NumEqVector(0.0); }
 
    NumEqVector fluxIntegral(const ElementDiscretization&, const ElementVariables&,
                             const auto&) const
    { return NumEqVector(0.0); }

    template<class SubControlVolume_>
    NumEqVector sourceIntegral(const ElementDiscretization& elemDisc,
                               const ElementVariables& elemVars,
                               const SubControlVolume_& scv) const
    {
        NumEqVector source(0.0);
 
        const auto& problem = this->asImp().problem();
        const auto& localBasis = elemDisc.gridDiscretization().feCache().get(elemDisc.element().type()).localBasis();
        std::vector< Dune::FieldVector<Scalar, 1> > shapeValues;
 
        for (const auto& qpData : CVFE::quadratureRule(elemDisc, scv, QuadratureRules::DuneQuadrature<3>{}))
        {
            const Tensor F = problem.couplingManager().deformationGradientAtPoint(
                elemDisc, qpData.ipData().global());
            const Scalar J = F.determinant();
            const Scalar psEq = problem.volumetricPressure(J);
 
            localBasis.evaluateFunction(qpData.ipData().local(), shapeValues);
            Scalar ps = 0.0;
            for (const auto& localDof : localDofs(elemDisc))
                ps += elemVars[localDof].pressure() * shapeValues[localDof.index()][0];
 
            source[0] -= (psEq - ps) * qpData.weight();
        }
 
        return source;
    }

    void addToElementFluxAndSourceResidual(ElementResidualVector& residual,
                                           const Problem& problem,
                                           const auto& element,
                                           const ElementDiscretization& elemDisc,
                                           const ElementVariables& elemVars) const
    {
        if constexpr (Dumux::Detail::LocalDofs::hasNonCVLocalDofsInterface<ElementDiscretization>())
        {
            if (nonCVLocalDofs(elemDisc).empty())
                return;
 
            const auto& problem = this->asImp().problem();
            const auto& localBasis = elemDisc.gridDiscretization().feCache().get(elemDisc.element().type()).localBasis();
            std::vector< Dune::FieldVector<Scalar, 1> > shapeValues;
 
            for (const auto& qpData : CVFE::quadratureRule(elemDisc, element, QuadratureRules::DuneQuadrature<3>{}))
            {
                const Tensor F = problem.couplingManager().deformationGradientAtPoint(
                    elemDisc, qpData.ipData().global());
                const Scalar J = F.determinant();
                const Scalar psEq = problem.volumetricPressure(J);
 
                localBasis.evaluateFunction(qpData.ipData().local(), shapeValues);
                Scalar ps = 0.0;
                for (const auto& localDof : localDofs(elemDisc))
                    ps += elemVars[localDof].pressure() * shapeValues[localDof.index()][0];
 
                for (const auto& nonCVdof : nonCVLocalDofs(elemDisc))
                {
                    const auto idx = nonCVdof.index();
                    residual[idx][0] -= (psEq - ps) * shapeValues[idx][0] * qpData.weight();
                }
            }
 
        }
    }
 
};
 
} // end namespace Dumux
#endif