pq3/fegriddiscretization.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_DISCRETIZATION_PQ3_FE_GRID_DISCRETIZATION_HH
#define DUMUX_DISCRETIZATION_PQ3_FE_GRID_DISCRETIZATION_HH
 
#include <ranges>
 
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
#include <dumux/discretization/method.hh>
#include <dumux/common/defaultmappertraits.hh>
#include <dumux/discretization/cvfe/quadraturerules.hh>
#include <dumux/discretization/fem/feelementdiscretization.hh>
#include <dumux/discretization/fem/fegriddiscretization.hh>
#include <dumux/discretization/periodicdofmap.hh>
#include <dumux/discretization/pq3/fvgridgeometry.hh>
 
#include "dofhelper.hh"
 
namespace Dumux::Experimental {

template<class GridView,
         class ElementRule = Dumux::QuadratureRules::DuneQuadrature<6>,
         class IntersectionRule = Dumux::QuadratureRules::DuneQuadrature<6>,
         class BoundaryFaceRule = Dumux::QuadratureRules::DuneQuadrature<6>>
struct PQ3FEQuadratureTraits
{
    using ElementQuadratureRule = ElementRule;
    using IntersectionQuadratureRule = IntersectionRule;
    using BoundaryFaceQuadratureRule = BoundaryFaceRule;
};

template<class GridView,
         class Scalar,
         class MapperTraits = Dumux::PQ3MapperTraits<GridView>,
         class QuadratureTraits = PQ3FEQuadratureTraits<GridView>>
struct PQ3FEDefaultGridDiscretizationTraits
: public MapperTraits, public QuadratureTraits
{
    using FeCache = Dune::LagrangeLocalFiniteElementCache<typename GridView::ctype, Scalar, GridView::dimension, 3>;
    using DofHelper = PQ3LagrangeDofHelper<GridView>;
 
    template<class GridDiscretization, bool enableCache>
    using LocalView = FEElementDiscretization<GridDiscretization, enableCache>;
 
    // Maximum number of element-local DOFs (for cubes, order 3)
    static constexpr std::size_t maxNumElementDofs = []()
        {
            if constexpr (GridView::dimension == 1)
                return 4;   // Q3 interval: 2 vertices + 2 edge interior
            else if constexpr (GridView::dimension == 2)
                return 16;  // Q3 quad: 4 vertices + 8 edge + 4 interior
            else if constexpr (GridView::dimension == 3)
                return 64;  // Q3 hex: 8 vertices + 24 edge + 24 face + 8 interior
        }();
};

template<class Scalar,
         class GV,
         bool enableCaching = true,
         class Traits = PQ3FEDefaultGridDiscretizationTraits<GV, Scalar>>
class PQ3FEGridDiscretization
: public BaseGridGeometry<GV, Traits>
{
    using ThisType = PQ3FEGridDiscretization<Scalar, GV, enableCaching, Traits>;
    using ParentType = BaseGridGeometry<GV, Traits>;
    using GridIndexType = typename IndexTraits<GV>::GridIndex;
    using LocalIndexType = typename IndexTraits<GV>::LocalIndex;
    static const int dim = GV::dimension;
 
public:
    using DiscretizationMethod = typename DiscretizationMethods::PQ3;
    static constexpr DiscretizationMethod discMethod{};
 
    static constexpr std::size_t maxNumElementDofs = Traits::maxNumElementDofs;

    using BoundaryFace = Experimental::BoundaryFace<GV>;

    using BasicGridGeometry = BasicGridGeometry_t<GV, Traits>;
    using LocalView = typename Traits::template LocalView<ThisType, true>;
    using Extrusion = Extrusion_t<Traits>;
    using DofMapper = typename Traits::DofMapper;
    using FeCache = typename Traits::FeCache;
    using GridView = GV;
    using SupportsPeriodicity = typename PeriodicGridTraits<typename GV::Grid>::SupportsPeriodicity;
    using ElementQuadratureRule = typename Traits::ElementQuadratureRule;
    using IntersectionQuadratureRule = typename Traits::IntersectionQuadratureRule;
    using BoundaryFaceQuadratureRule = typename Traits::BoundaryFaceQuadratureRule;

    PQ3FEGridDiscretization(std::shared_ptr<BasicGridGeometry> gd)
    : ParentType(std::move(gd))
    , dofMapper_(this->gridView(), Traits::layout())
    , cache_(*this)
    , periodicGridTraits_(this->gridView().grid())
    {
        update_();
    }

    PQ3FEGridDiscretization(const GridView& gridView)
    : PQ3FEGridDiscretization(std::make_shared<BasicGridGeometry>(gridView))
    {}

    const DofMapper& dofMapper() const
    { return dofMapper_; }

    std::size_t numDofs() const
    { return this->dofMapper().size(); }

    void update(const GridView& gridView)
    {
        ParentType::update(gridView);
        update_();
    }

    void update(GridView&& gridView)
    {
        ParentType::update(std::move(gridView));
        update_();
    }

    const FeCache& feCache() const
    { return feCache_; }

    template<class LocalKey>
    auto dofIndex(const typename GV::template Codim<0>::Entity& element,
                  const LocalKey& localKey) const
    { return Traits::DofHelper::dofIndex(dofMapper_, element, localKey, this->gridView().grid().globalIdSet()); }

    bool dofOnBoundary(GridIndexType dofIdx) const
    { return boundaryDofIndices_[dofIdx]; }

    bool dofOnPeriodicBoundary(GridIndexType dofIdx) const
    { return periodicDofMap_.count(dofIdx); }

    [[deprecated("Will be removed after release 3.11. Use periodicallyMappedDofs, returning a range of dofs")]]
    GridIndexType periodicallyMappedDof(GridIndexType dofIdx) const
    { return periodicDofMap_.at(dofIdx)[0]; }

    const std::ranges::range auto periodicallyMappedDofs(GridIndexType dofIdx) const
    { return std::views::all(periodicDofMap_.at(dofIdx)); }

    const std::unordered_map<GridIndexType, std::vector<GridIndexType>>& periodicDofMap() const
    { return periodicDofMap_; }

    friend inline LocalView localView(const PQ3FEGridDiscretization& gd)
    { return { gd.cache_ }; }
 
private:
    class PQ3FEGridDiscretizationCache
    {
        friend class PQ3FEGridDiscretization;
    public:
        using DofHelper = typename Traits::DofHelper;
 
        explicit PQ3FEGridDiscretizationCache(const PQ3FEGridDiscretization& gd)
        : gridDiscretization_(&gd)
        {}
 
        const PQ3FEGridDiscretization& gridDiscretization() const
        { return *gridDiscretization_; }

        auto boundaryFaces(GridIndexType eIdx) const -> std::span<const BoundaryFace>
        {
            if (auto it = boundaryFaces_.find(eIdx); it != boundaryFaces_.end())
                return {it->second};
            return {};
        }
 
    private:
        void clear_()
        { boundaryFaces_.clear(); }
 
        std::unordered_map<GridIndexType, Dune::ReservedVector<BoundaryFace, 2*dim>> boundaryFaces_;
 
        const PQ3FEGridDiscretization* gridDiscretization_;
    };
 
public:
    using Cache = PQ3FEGridDiscretizationCache;
 
private:
    void update_()
    {
        using DofHelper = typename Cache::DofHelper;
 
        dofMapper_.update(this->gridView());
        cache_.clear_();
        boundaryDofIndices_.assign(numDofs(), false);
        periodicDofMap_.clear();
        const auto& idSet = this->gridView().grid().globalIdSet();
 
        for (const auto& element : elements(this->gridView()))
        {
            auto elementGeometry = element.geometry();
            const auto& localCoefficients = this->feCache().get(element.type()).localCoefficients();
            const auto eIdx = this->elementMapper().index(element);
 
            LocalIndexType numBoundaryFaces = 0;
            for (const auto& intersection : intersections(this->gridView(), element))
            {
                if (intersection.boundary() && !periodicGridTraits_.isPeriodic(intersection))
                {
                    // add one boundary face per boundary intersection
                    const auto isGeometry = intersection.geometry();
                    cache_.boundaryFaces_[eIdx].push_back(BoundaryFace{
                        isGeometry.center(),
                        isGeometry.volume(),
                        intersection.centerUnitOuterNormal(),
                        numBoundaryFaces++,
                        static_cast<LocalIndexType>(intersection.indexInInside()),
                        typename BoundaryFace::Traits::BoundaryFlag{intersection}
                    });
 
                    // mark all dofs on this intersection as boundary dofs
                    for (LocalIndexType keyIdx = 0; keyIdx < localCoefficients.size(); ++keyIdx)
                    {
                        if (DofHelper::localDofOnIntersection(elementGeometry.type(),
                                                              intersection.indexInInside(),
                                                              localCoefficients.localKey(keyIdx)))
                        {
                            const auto dofIdxGlobal = DofHelper::dofIndex(
                                this->dofMapper(), element, localCoefficients.localKey(keyIdx), idSet);
                            boundaryDofIndices_[dofIdxGlobal] = true;
                        }
                    }
                }
 
                // inform the grid discretization if we have periodic boundaries
                else if (periodicGridTraits_.isPeriodic(intersection))
                {
                    this->setPeriodic();
 
                    const auto eps = 1e-7*(elementGeometry.corner(1) - elementGeometry.corner(0)).two_norm();
                    for (int localDofIdx = 0; localDofIdx < localCoefficients.size(); ++localDofIdx)
                    {
                        if (!DofHelper::localDofOnIntersection(elementGeometry.type(),
                                                               intersection.indexInInside(),
                                                               localCoefficients.localKey(localDofIdx)))
                            continue;
 
                        const auto dofIdxGlobal = DofHelper::dofIndex(
                            this->dofMapper(), element, localCoefficients.localKey(localDofIdx), idSet);
                        const auto dofPos = DofHelper::dofPosition(elementGeometry, localCoefficients.localKey(localDofIdx));
 
                        const auto& outside = intersection.outside();
                        const auto outsideGeometry = outside.geometry();
                        const auto& localCoefficientsOut = this->feCache().get(outsideGeometry.type()).localCoefficients();
                        for (const auto& isOutside : intersections(this->gridView(), outside))
                        {
                            if (periodicGridTraits_.isPeriodic(isOutside))
                            {
                                for (int localDofIdxOut = 0; localDofIdxOut < localCoefficientsOut.size(); ++localDofIdxOut)
                                {
                                    const auto& localKeyOut = localCoefficientsOut.localKey(localDofIdxOut);
                                    if (!DofHelper::localDofOnIntersection(outsideGeometry.type(),
                                                                           isOutside.indexInInside(),
                                                                           localKeyOut))
                                        continue;
 
                                    const auto dofIdxGlobalOut = DofHelper::dofIndex(this->dofMapper(), outside, localKeyOut, idSet);
                                    const auto dofPosOutside = DofHelper::dofPosition(outsideGeometry, localKeyOut);
                                    const auto shift = std::abs((this->bBoxMax()-this->bBoxMin())*intersection.centerUnitOuterNormal());
                                    if ((dofPosOutside - dofPos + shift*intersection.centerUnitOuterNormal()).two_norm() < eps)
                                        Dumux::Detail::addPeriodicallyMappedDof(periodicDofMap_, dofIdxGlobal, dofIdxGlobalOut);
                                }
                            }
                        }
                    }
                }
            }
        }
 
        Dumux::Detail::closePeriodicDofMap(periodicDofMap_);
 
        // error check: periodic boundaries currently don't work in parallel
        if (this->isPeriodic() && this->gridView().comm().size() > 1)
            DUNE_THROW(Dune::NotImplemented, "Periodic boundaries for FE PQ3 method for parallel simulations!");
    }
 
    DofMapper dofMapper_;
    const FeCache feCache_;
    // dofs on the boundary
    std::vector<bool> boundaryDofIndices_;
    // a map for periodic boundary dofs
    std::unordered_map<GridIndexType, std::vector<GridIndexType>> periodicDofMap_;
    Cache cache_;
    PeriodicGridTraits<typename GridView::Grid> periodicGridTraits_;
};
 
} // end namespace Dumux::Experimental
 
#endif