pq1bubblelocalfiniteelement.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_PQ1BUBBLE_LOCAL_FINITE_ELEMENT_HH
#define DUMUX_DISCRETIZATION_PQ1BUBBLE_LOCAL_FINITE_ELEMENT_HH
 
#include <array>
#include <vector>
#include <numeric>
#include <algorithm>
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
 
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
 
#include <dune/localfunctions/common/localbasis.hh>
#include <dune/localfunctions/common/localfiniteelementtraits.hh>
#include <dune/localfunctions/common/localkey.hh>
 
#include <dune/localfunctions/lagrange/lagrangesimplex.hh>
#include <dune/localfunctions/lagrange/lagrangecube.hh>
 
namespace Dumux::Detail {

template<class D, class R, unsigned int dim, Dune::GeometryType::Id typeId, std::size_t numCubeBubbleDofs>
class PQ1BubbleLocalBasis
{
    using PQ1FiniteElement = std::conditional_t<
        Dune::GeometryType{ typeId } == Dune::GeometryTypes::cube(dim),
        Dune::LagrangeCubeLocalFiniteElement<D, R, dim, 1>,
        Dune::LagrangeSimplexLocalFiniteElement<D, R, dim, 1>
    >;
    static constexpr bool isCube = Dune::GeometryType{ typeId } == Dune::GeometryTypes::cube(dim);
    static constexpr std::size_t numVertexDofs = isCube ? (1<<dim) : (dim+1);
    static constexpr std::size_t numBubbleDofs = isCube ? numCubeBubbleDofs : 1;
    static constexpr std::size_t numDofs = numVertexDofs + numBubbleDofs;
public:
    using Traits = Dune::LocalBasisTraits<
        D, dim, Dune::FieldVector<D, dim>,
        R, 1, Dune::FieldVector<R, 1>,
        Dune::FieldMatrix<R, 1, dim>
    >;
 
    PQ1BubbleLocalBasis()
    : pq1FiniteElement_{}
    {
        // precompute center values for normalization
        const auto& p1Basis = pq1FiniteElement_.localBasis();
        const auto refElement = Dune::referenceElement<typename Traits::DomainFieldType, dim>(type());
        const auto& center = refElement.position(0, 0);
        p1Basis.evaluateFunction(center, pq1AtCenter_);
    }

    static constexpr unsigned int size()
    { return numDofs; }

    void evaluateFunction(const typename Traits::DomainType& x,
                          std::vector<typename Traits::RangeType>& out) const
    {
        out.reserve(size());
        const auto& p1Basis = pq1FiniteElement_.localBasis();
        p1Basis.evaluateFunction(x, out);
        const auto bubble = evaluateBubble_(x);
 
        out.resize(size());
        out[numVertexDofs] = bubble;
 
        if constexpr (numBubbleDofs == 2)
            out[numVertexDofs+1] = scaling_(x)*bubble;
 
        for (std::size_t i = 0; i < numVertexDofs; ++i)
            out[i] -= pq1AtCenter_[i]*out[numVertexDofs];
    }

    void evaluateJacobian(const typename Traits::DomainType& x,
                          std::vector<typename Traits::JacobianType>& out) const
    {
        out.reserve(size());
        const auto& p1Basis = pq1FiniteElement_.localBasis();
        p1Basis.evaluateJacobian(x, out);
 
        std::vector<typename Traits::RangeType> shapeValues;
        p1Basis.evaluateFunction(x, shapeValues);
 
        const auto bubble = evaluateBubble_(x);
        const auto bubbleJacobian = evaluateBubbleJacobian_(x);
 
        for (std::size_t i = 0; i < numVertexDofs; ++i)
            for (int k = 0; k < dim; ++k)
                out[i][0][k] -= pq1AtCenter_[i][0]*bubbleJacobian[0][k];
 
        out.resize(size());
        out[numVertexDofs] = bubbleJacobian;
        if constexpr (numBubbleDofs == 2)
            out[numVertexDofs+1] = scaling_(x)*bubbleJacobian + bubble[0]*scalingJacobian_(x);
    }

    void partial(const std::array<unsigned int, dim>& order,
                 const typename Traits::DomainType& in,
                 std::vector<typename Traits::RangeType>& out) const
    {
        DUNE_THROW(Dune::NotImplemented, "Partial derivatives");
    }

    static constexpr unsigned int order()
    {
        return dim;
    }

    static constexpr Dune::GeometryType type()
    {
        return { typeId };
    }
private:
    // evaluate bubble function at x
    typename Traits::RangeType evaluateBubble_(const typename Traits::DomainType& x) const
    {
        if constexpr (type() == Dune::GeometryTypes::simplex(dim))
        {
            if constexpr (dim == 2)
                return 27*x[0]*x[1]*(1-x[0]-x[1]);
            else if constexpr (dim == 3)
                return 256*x[0]*x[1]*x[2]*(1-x[0]-x[1]-x[2]);
        }
        else if constexpr (type() == Dune::GeometryTypes::cube(dim))
        {
            if constexpr (dim == 2)
                return 16*x[0]*x[1]*(1-x[0])*(1-x[1]);
            else if constexpr (dim == 3)
                return 64*x[0]*x[1]*x[2]*(1-x[0])*(1-x[1])*(1-x[2]);
        }
        else
            DUNE_THROW(Dune::NotImplemented, "Bubble function for " << type());
    }
 
    // evaluate bubble function at x
    typename Traits::JacobianType evaluateBubbleJacobian_(const typename Traits::DomainType& x) const
    {
        if constexpr (type() == Dune::GeometryTypes::simplex(dim))
        {
            if constexpr (dim == 2)
                return {{27*(x[1]*(1-x[0]-x[1]) - x[0]*x[1]),
                         27*(x[0]*(1-x[0]-x[1]) - x[0]*x[1])}};
            else if constexpr (dim == 3)
                return {{256*(x[1]*x[2]*(1-x[0]-x[1]-x[2]) - x[0]*x[1]*x[2]),
                         256*(x[0]*x[2]*(1-x[0]-x[1]-x[2]) - x[0]*x[1]*x[2]),
                         256*(x[0]*x[1]*(1-x[0]-x[1]-x[2]) - x[0]*x[1]*x[2])}};
        }
        else if constexpr (type() == Dune::GeometryTypes::cube(dim))
        {
            if constexpr (dim == 2)
                return {{16*(x[1]*(1-x[0])*(1-x[1]) - x[0]*x[1]*(1-x[1])),
                         16*(x[0]*(1-x[0])*(1-x[1]) - x[0]*x[1]*(1-x[0]))}};
            else if constexpr (dim == 3)
                // by the product rule,
                //   64*(x[1]*x[2]*(1-x[0])*(1-x[1])*(1-x[2]) - x[0]*x[1]*x[2]*(1-x[1])*(1-x[2]))
                //   64*(x[0]*x[2]*(1-x[0])*(1-x[1])*(1-x[2]) - x[0]*x[1]*x[2]*(1-x[0])*(1-x[2]))
                //   64*(x[0]*x[1]*(1-x[0])*(1-x[1])*(1-x[2]) - x[0]*x[1]*x[2]*(1-x[0])*(1-x[1]))
                // which simplifies to:
                return {{64*x[1]*x[2]*(1-x[1])*(1-x[2])*(1-2*x[0]),
                         64*x[0]*x[2]*(1-x[0])*(1-x[2])*(1-2*x[1]),
                         64*x[0]*x[1]*(1-x[0])*(1-x[1])*(1-2*x[2])}};
        }
        else
            DUNE_THROW(Dune::NotImplemented, "Bubble function for " << type() << " dim = " << dim);
    }
 
    auto scaling_(const typename Traits::DomainType& x) const
    {
        if constexpr (dim == 2)
            return (x[0]-0.5) + (x[1]-0.5);
        else if constexpr (dim == 3)
            return (x[0]-0.5) + (x[1]-0.5) + (x[2]-0.5);
        else
            DUNE_THROW(Dune::NotImplemented, "Bubble function function for dim = " << dim);
    }
 
    typename Traits::JacobianType scalingJacobian_(const typename Traits::DomainType& x) const
    {
        if constexpr (dim == 2)
            return {{1.0, 1.0}};
        else if constexpr (dim == 3)
            return {{1.0, 1.0, 1.0}};
        else
            DUNE_THROW(Dune::NotImplemented, "Bubble function function for dim = " << dim);
    }
 
    PQ1FiniteElement pq1FiniteElement_;
    std::vector<typename Traits::RangeType> pq1AtCenter_;
};

template<int dim, Dune::GeometryType::Id typeId, std::size_t numCubeBubbleDofs>
class PQ1BubbleLocalCoefficients
{
    static constexpr bool isCube = Dune::GeometryType{ typeId } == Dune::GeometryTypes::cube(dim);
    static constexpr std::size_t numVertexDofs = isCube ? (1<<dim) : (dim+1);
    static constexpr std::size_t numBubbleDofs = isCube ? numCubeBubbleDofs : 1;
    static constexpr std::size_t numDofs = numVertexDofs + numBubbleDofs;
public:
 
    PQ1BubbleLocalCoefficients()
    {
        // Dune::LocalKey is constructed with
        // - the local index with respect to the subentity type
        // - the codim of the subentity the dof is attached to
        // - the local number of the function to evaluate
 
        // vertex dofs
        for (std::size_t i = 0; i < numVertexDofs; ++i)
            localKeys_[i] = Dune::LocalKey(i, dim, 0);
 
        // cell dof(s)
        for (std::size_t i = 0; i < numBubbleDofs; ++i)
            localKeys_[numVertexDofs + i] = Dune::LocalKey(0, 0, i);
    }

    static constexpr std::size_t size ()
    { return numDofs; }

    const Dune::LocalKey& localKey (std::size_t i) const
    { return localKeys_[i]; }
 
private:
    std::array<Dune::LocalKey, numDofs> localKeys_;
};

template<class LocalBasis>
class PQ1BubbleLocalInterpolation
{
public:
    template<typename F, typename C>
    void interpolate (const F& f, std::vector<C>& out) const
    {
        constexpr auto dim = LocalBasis::Traits::dimDomain;
        if constexpr (LocalBasis::type() == Dune::GeometryTypes::cube(dim))
            DUNE_THROW(Dune::NotImplemented, "Interpolation for Q1 bubble not implemented yet");
 
        out.resize(LocalBasis::size());
 
        const auto refElement = Dune::referenceElement<typename LocalBasis::Traits::DomainFieldType,dim>(LocalBasis::type());
 
        // f may return a scalar or a vector with a single entry
        const auto coefficient = [](const auto& value) -> C
        {
            if constexpr (Dune::IsNumber<std::decay_t<decltype(value)>>::value)
                return value;
            else
                return value[0];
        };
 
        // Evaluate at the vertices and at the center
        for (int i = 0; i < refElement.size(dim); ++i)
            out[i] = coefficient(f(refElement.position(i, dim)));
        out.back() = coefficient(f(refElement.position(0, 0)));
    }
};
 
} // end namespace Dumux::Detail
 
namespace Dumux {

template<class D, class R, int dim, Dune::GeometryType::Id typeId, std::size_t numCubeBubbleDofs = 1>
class PQ1BubbleLocalFiniteElement
{
    static_assert(numCubeBubbleDofs == 1 || numCubeBubbleDofs == 2,
                    "P1/Q1 bubble FE supports numCubeBubbleDofs = 1 or 2");
    using Basis = Detail::PQ1BubbleLocalBasis<D, R, dim, typeId, numCubeBubbleDofs>;
    using Coefficients = Detail::PQ1BubbleLocalCoefficients<dim, typeId, numCubeBubbleDofs>;
    using Interpolation = Detail::PQ1BubbleLocalInterpolation<Basis>;
 
    static constexpr Dune::GeometryType gt = Dune::GeometryType{ typeId };
    static_assert(
        gt == Dune::GeometryTypes::cube(dim) || gt == Dune::GeometryTypes::simplex(dim),
        "Only implemented for cubes and simplices"
    );
 
public:
    using Traits = Dune::LocalFiniteElementTraits<Basis, Coefficients, Interpolation>;

    const typename Traits::LocalBasisType& localBasis() const
    {
        return basis_;
    }

    const typename Traits::LocalCoefficientsType& localCoefficients() const
    {
        return coefficients_;
    }

    const typename Traits::LocalInterpolationType& localInterpolation() const
    {
        return interpolation_;
    }

    static constexpr std::size_t size()
    {
        return Basis::size();
    }

    static constexpr Dune::GeometryType type()
    {
        return Traits::LocalBasisType::type();
    }
 
private:
    Basis basis_;
    Coefficients coefficients_;
    Interpolation interpolation_;
};
 
} // end namespace Dumux
 
#endif