quadraturerulecache.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef QUADRATURERULECACHE_HH
#define QUADRATURERULECACHE_HH
 
#include <mutex>
#include <map>
#include <algorithm>
 
#include <dune/geometry/type.hh>
 
#include <dune/geometry/quadraturerules.hh>
 
#include <dune/fufem/quadraturerules/refinedquadraturerule.hh>
#include <dune/fufem/quadraturerules/compositequadraturerule.hh>
#include <dune/fufem/quadraturerules/lumpingquadraturerule.hh>
 
#include <dune/fufem/refinedfehelper.hh>
 
 
 
class QuadratureRuleKey
{
    public:
 
        QuadratureRuleKey() :
            gt_(Dune::GeometryTypes::none(0)),
            order_(0),
            refinement_(0),
            lumping_(false)
        {}
 
        template<class LocalFiniteElement>
        QuadratureRuleKey(const LocalFiniteElement& fe)
        {
            gt_ = fe.type();
            order_ = fe.localBasis().order();
            refinement_ = IsRefinedLocalFiniteElement<LocalFiniteElement>::value(fe);
            lumping_ = false;
        }
 
        QuadratureRuleKey(const Dune::GeometryType& gt, const int order, const int refinement=0, bool lumping=false) :
            gt_(gt),
            order_(order),
            refinement_(refinement),
            lumping_(lumping)
        {}
 
        QuadratureRuleKey(const int dim, const int order, const int refinement=0, bool lumping=false) :
            gt_(Dune::GeometryTypes::none(dim)),
            order_(order),
            refinement_(refinement),
            lumping_(lumping)
        {}
 
        bool operator<(const QuadratureRuleKey& other) const
        {
            if (gt_ < other.gt_)
                return true;
            if (gt_ == other.gt_)
            {
                if (order_ < other.order_)
                    return true;
                if (order_ == other.order_)
                {
                    if (refinement_ < other.refinement_)
                        return true;
                    if (refinement_ == other.refinement_)
                        return lumping_ < other.lumping_;
                }
            }
            return false;
        }
 
        const Dune::GeometryType& geometryType() const
        {
            return gt_;
        }
 
        void setGeometryType(const Dune::GeometryType& gt)
        {
            gt_ = gt;
        }
 
        int order() const
        {
            return order_;
        }
 
        void setOrder(int order)
        {
            order_ = order;
        }
 
        int refinement() const
        {
            return refinement_;
        }
 
        void setRefinement(int refinement)
        {
            refinement_ = refinement;
        }
 
        bool lumping() const
        {
            return lumping_;
        }
 
        void setLumping(bool lumping)
        {
            lumping_ = lumping;
        }
 
 
 
        QuadratureRuleKey derivative() const
        {
            int derivativeOrder = gt_.isSimplex() ? std::max(0,order_-1) : order_;
            return QuadratureRuleKey(gt_, derivativeOrder, refinement_, lumping_);
        }
 
        QuadratureRuleKey sum(const QuadratureRuleKey& other) const
        {
            Dune::GeometryType resultGT = gt_.isNone() ? other.gt_ : gt_;
            return QuadratureRuleKey(resultGT, std::max(order_, other.order_), std::max(refinement_, other.refinement_), (lumping_ or other.lumping_));
        }
 
        QuadratureRuleKey product(const QuadratureRuleKey& other) const
        {
            Dune::GeometryType resultGT = gt_.isNone() ? other.gt_ : gt_;
            return QuadratureRuleKey(resultGT, order_ + other.order_, std::max(refinement_, other.refinement_), (lumping_ or other.lumping_));
        }
 
        QuadratureRuleKey square() const
        {
            return QuadratureRuleKey(gt_, 2*order_, refinement_, lumping_);
        }
 
    protected:
 
        Dune::GeometryType gt_;
        int order_;
        int refinement_;
        bool lumping_;
};
 
 
 
 
template <class coord_type, int dim>
class QuadratureRuleCache
{
        typedef std::map<QuadratureRuleKey, Dune::QuadratureRule<coord_type,dim> > ContainerType;
        static ContainerType  quadRules;
        static std::mutex mut;
 
    public:
        static const Dune::QuadratureRule<coord_type,dim>& rule(const Dune::GeometryType& gt, const int order, int refinement)
        {
            return rule(QuadratureRuleKey(gt, order, refinement));
        }
 
        static const Dune::QuadratureRule<coord_type,dim>& rule(const QuadratureRuleKey& index)
        {
            const Dune::GeometryType& gt = index.geometryType();
            const int order = index.order();
            const int refinement = index.refinement();
            const bool lumping = index.lumping();
 
            typename ContainerType::iterator qrule_it = quadRules.find(index);
 
            if (qrule_it != quadRules.end())
                return qrule_it->second;
 
            // if quadrature rule is not already cached, create it and insert it in the cache
            if (lumping)
            {
                auto lock = std::unique_lock(mut);
                qrule_it = (quadRules.insert(std::make_pair(QuadratureRuleKey(gt,order,refinement,true), LumpingQuadratureRule<coord_type,dim>(gt)))).first;
                return qrule_it->second;
            }
            if (refinement==0)
            {
                return Dune::template QuadratureRules<coord_type, dim>::rule(gt, order);
            }
            if ((refinement==1) and (gt.isSimplex()))
            {
                auto lock = std::unique_lock(mut);
                RefinedSimplexQuadratureRule<coord_type,dim> qrule( Dune::QuadratureRules<coord_type, dim>::rule(gt, order) );
                qrule_it = (quadRules.insert(std::make_pair(index, qrule))).first;
                return qrule_it->second;
            }
            if ((refinement>1) and (gt.isSimplex()))
            {
                auto lock = std::unique_lock(mut);
                CompositeQuadratureRule<coord_type,dim> qrule( Dune::QuadratureRules<coord_type, dim>::rule(gt, order), refinement);
                qrule_it = (quadRules.insert(std::make_pair(index, qrule))).first;
                return qrule_it->second;
            }
            DUNE_THROW(Dune::NotImplemented, "Quadrature rule for given QuadratureRuleKey(" << gt << "," << order << "," << refinement << " is not implemented!");
        }
};
// define static template members
template <class coord_type, int dim>  std::map<QuadratureRuleKey, Dune::QuadratureRule<coord_type,dim> > QuadratureRuleCache<coord_type,dim>::quadRules;
template <class coord_type, int dim>  std::mutex QuadratureRuleCache<coord_type,dim>::mut;
 
 
#endif