projector.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_PROJECTOR_HH
#define DUMUX_DISCRETIZATION_PROJECTOR_HH
 
#include <algorithm>
#include <string>
#include <utility>
#include <type_traits>
 
#include <dune/common/timer.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/promotiontraits.hh>
#include <dune/common/parametertree.hh>
 
#include <dune/geometry/quadraturerules.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
 
#include <dumux/io/format.hh>
#include <dumux/common/parameters.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/linearalgebratraits.hh>
#include <dumux/linear/istlsolvers.hh>
#include <dumux/assembly/jacobianpattern.hh>
#include <dumux/discretization/functionspacebasis.hh>
 
namespace Dumux {

template<class ScalarType>
class Projector
{
    using MatrixBlockType = Dune::FieldMatrix<ScalarType, 1, 1>;
 
public:
    using Scalar = ScalarType;
    using Matrix = Dune::BCRSMatrix< MatrixBlockType >;

    struct Params
    {
        std::size_t maxIterations{100};
        Scalar residualReduction{1e-13};
        int verbosity{0};
    };

    Projector() = delete;

    Projector(Matrix&& massMatrix, Matrix&& projectionMatrix)
    : massMat_(std::make_shared<Matrix>(std::move(massMatrix)))
    , projMat_(std::make_shared<Matrix>(std::move(projectionMatrix)))
    , numDofsTarget_(massMat_->N())
    {
        if (massMat_->N() != projMat_->N())
            DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_->N() << " vs " << projMat_->N());
 
        massMatrixSolver_.setMatrix(massMat_);
    }

    Projector(Matrix&& massMatrix,
              Matrix&& projectionMatrix,
              std::vector<std::size_t>&& indexMap,
              std::size_t numDofsTarget)
    : massMat_(std::make_shared<Matrix>(std::move(massMatrix)))
    , projMat_(std::make_shared<Matrix>(std::move(projectionMatrix)))
    , indexMapTarget_(std::move(indexMap))
    , numDofsTarget_(numDofsTarget)
    {
        if (indexMapTarget_.size() != massMat_->N())
            DUNE_THROW(Dune::InvalidStateException, "Target index map size mismatch: " << indexMapTarget_.size() << " vs " << massMat_->N());
 
        if (massMat_->N() != projMat_->N())
            DUNE_THROW(Dune::InvalidStateException, "Matrix row size mismatch: " << massMat_->N() << " vs " << projMat_->N());
 
        if (!indexMapTarget_.empty())
            if (*std::max_element(indexMapTarget_.begin(), indexMapTarget_.end()) > numDofsTarget_)
                DUNE_THROW(Dune::InvalidStateException, "Index map exceeds provided number of dofs in target domain!");
 
        massMatrixSolver_.setMatrix(massMat_);
    }

    template< class BlockType, std::enable_if_t<std::is_convertible<BlockType, ScalarType>::value, int> = 0 >
    Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const
    {
        // be picky about size of u
        if ( u.size() != projMat_->M())
            DUNE_THROW(Dune::InvalidStateException, "Vector size mismatch");
 
        Dune::BlockVector<BlockType> up(massMat_->N());
 
        auto rhs = up;
        projMat_->mv(u, rhs);
 
        Dune::ParameterTree solverParams;
        solverParams["maxit"] = std::to_string(params.maxIterations);
        solverParams["reduction"] = Fmt::format("{}", params.residualReduction);
        solverParams["verbose"] = std::to_string(params.verbosity);
        auto solver = massMatrixSolver_; // copy the solver to modify the parameters
        solver.setParams(solverParams);
        solver.solve(up, rhs);
 
        // if target space occupies a larger region, fill missing entries with zero
        if (!indexMapTarget_.empty())
        {
            Dune::BlockVector<BlockType> result(numDofsTarget_);
 
            result = 0.0;
            for (std::size_t i = 0; i < indexMapTarget_.size(); ++i)
                result[indexMapTarget_[i]] = up[i];
 
            return result;
        }
 
        return up;
    }

    template< class BlockType, std::enable_if_t<!std::is_convertible<BlockType, ScalarType>::value, int> = 0 >
    Dune::BlockVector<BlockType> project(const Dune::BlockVector<BlockType>& u, const Params& params = Params{}) const
    {
        Dune::BlockVector<BlockType> result(numDofsTarget_);
 
        for (int pvIdx = 0; pvIdx < BlockType::size(); ++pvIdx)
        {
            Dune::BlockVector<Dune::FieldVector<Scalar, 1>> tmp(u.size());
            std::transform(u.begin(), u.end(), tmp.begin(), [pvIdx] (const auto& v) { return v[pvIdx]; });
 
            const auto p = project(tmp, params);
            for (std::size_t i = 0; i < p.size(); ++i)
                result[i][pvIdx] = p[i];
        }
 
        return result;
    }

    static Params defaultParams()
    { return {}; }
 
private:
    std::shared_ptr<Matrix> massMat_;
    std::shared_ptr<Matrix> projMat_;
 
    SSORCGIstlSolver<
        SeqLinearSolverTraits, LinearAlgebraTraits<Matrix, Dune::BlockVector<Dune::FieldVector<Scalar, 1>>>
    > massMatrixSolver_;
 
    std::vector<std::size_t> indexMapTarget_;
    std::size_t numDofsTarget_;
};

template<class FEBasisDomain, class FEBasisTarget>
class ProjectorTraits
{
    using FiniteElementDomain = typename FEBasisDomain::LocalView::Tree::FiniteElement;
    using FiniteElementTarget = typename FEBasisTarget::LocalView::Tree::FiniteElement;
    using ScalarDomain = typename FiniteElementDomain::Traits::LocalBasisType::Traits::RangeFieldType;
    using ScalarTarget = typename FiniteElementTarget::Traits::LocalBasisType::Traits::RangeFieldType;
    using Scalar = typename Dune::PromotionTraits<ScalarDomain, ScalarTarget>::PromotedType;
public:
    using Projector = Dumux::Projector< Scalar >;
};
 
 
// Projector construction details
namespace Detail {

template<class Matrix>
void setupReducedMatrices(const Matrix& massMatrix, const Matrix& projMatrix, const std::vector<bool>& dofIsVoid,
                          Matrix& reducedM, Matrix& reducedP, std::vector<std::size_t>& expansionMap)
{
    const std::size_t numNonVoidDofs = std::count_if(dofIsVoid.begin(), dofIsVoid.end(), [] (bool v) { return !v; });
 
    // reduce matrices to only dofs that take part and create index map
    std::vector<std::size_t> reductionMap(massMatrix.N());
    expansionMap.resize(numNonVoidDofs);
 
    std::size_t idxInReducedSpace = 0;
    for (std::size_t dofIdx = 0; dofIdx < dofIsVoid.size(); ++dofIdx)
        if (!dofIsVoid[dofIdx])
        {
            reductionMap[dofIdx] = idxInReducedSpace;
            expansionMap[idxInReducedSpace] = dofIdx;
            idxInReducedSpace++;
        }
 
    // create reduced M/P matrix
    Dune::MatrixIndexSet patternMReduced, patternPReduced;
    patternMReduced.resize(numNonVoidDofs, numNonVoidDofs);
    patternPReduced.resize(numNonVoidDofs, projMatrix.M());
    for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt)
        if (!dofIsVoid[rowIt.index()])
        {
            const auto reducedRowIdx = reductionMap[rowIt.index()];
            for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
                if (!dofIsVoid[colIt.index()])
                    patternMReduced.add(reducedRowIdx, reductionMap[colIt.index()]);
        }
 
    for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt)
        if (!dofIsVoid[rowIt.index()])
        {
            const auto reducedRowIdx = reductionMap[rowIt.index()];
            for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
                patternPReduced.add(reducedRowIdx, colIt.index());
        }
 
    patternMReduced.exportIdx(reducedM);
    patternPReduced.exportIdx(reducedP);
 
    // fill matrix entries
    for (auto rowIt = massMatrix.begin(); rowIt != massMatrix.end(); ++rowIt)
        if (!dofIsVoid[rowIt.index()])
        {
            const auto reducedRowIdx = reductionMap[rowIt.index()];
            for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
                if (!dofIsVoid[colIt.index()])
                    reducedM[reducedRowIdx][reductionMap[colIt.index()]] = *colIt;
        }
 
    for (auto rowIt = projMatrix.begin(); rowIt != projMatrix.end(); ++rowIt)
        if (!dofIsVoid[rowIt.index()])
        {
            const auto reducedRowIdx = reductionMap[rowIt.index()];
            for (auto colIt = (*rowIt).begin(); colIt != (*rowIt).end(); ++colIt)
                reducedP[reducedRowIdx][colIt.index()] = *colIt;
        }
}

template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType>
auto createProjectionMatrices(const FEBasisDomain& feBasisDomain,
                              const FEBasisTarget& feBasisTarget,
                              const GlueType& glue,
                              bool treatDiagonalZeroes = true)
{
    // we assume that target dim <= domain dimension
    static constexpr int domainDim = FEBasisDomain::GridView::dimension;
    static constexpr int targetDim = FEBasisTarget::GridView::dimension;
    static_assert(targetDim <= domainDim, "This expects target dim < domain dim, please swap arguments");
 
    using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector;
    using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector;
 
    using ForwardProjectionMatrix = typename ForwardProjector::Matrix;
    using BackwardProjectionMatrix = typename BackwardProjector::Matrix;
 
    auto domainLocalView = feBasisDomain.localView();
    auto targetLocalView = feBasisTarget.localView();
 
    // determine mass matrix patterns (standard FE scheme pattern)
    Dune::MatrixIndexSet backwardPatternM, forwardPatternM;
    forwardPatternM = getFEJacobianPattern(feBasisTarget);
    if (doBidirectional) backwardPatternM = getFEJacobianPattern(feBasisDomain);
 
    // determine projection matrix patterns
    Dune::MatrixIndexSet backwardPatternP, forwardPatternP;
    forwardPatternP.resize(feBasisTarget.size(), feBasisDomain.size());
    if (doBidirectional) backwardPatternP.resize(feBasisDomain.size(), feBasisTarget.size());
 
    using std::max;
    unsigned int maxBasisOrder = 0;
    for (const auto& is : intersections(glue))
    {
        // since target dim <= domain dim there is maximum one!
        targetLocalView.bind( is.targetEntity(0) );
        const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis();
 
        for (unsigned int nIdx = 0; nIdx < is.numDomainNeighbors(); ++nIdx)
        {
            domainLocalView.bind( is.domainEntity(nIdx) );
            const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis();
 
            // keep track of maximum basis order (used in integration)
            maxBasisOrder = max(maxBasisOrder, max(domainLocalBasis.order(), targetLocalBasis.order()));
 
            for (unsigned int i = 0; i < domainLocalBasis.size(); ++i)
                for (unsigned int j = 0; j < targetLocalBasis.size(); ++j)
                {
                    forwardPatternP.add(targetLocalView.index(j), domainLocalView.index(i));
                    if (doBidirectional) backwardPatternP.add(domainLocalView.index(i), targetLocalView.index(j));
                }
        }
    }
 
    // assemble matrices
    ForwardProjectionMatrix forwardM, forwardP;
    forwardPatternM.exportIdx(forwardM); forwardM = 0.0;
    forwardPatternP.exportIdx(forwardP); forwardP = 0.0;
 
    BackwardProjectionMatrix backwardM, backwardP;
    if (doBidirectional)
    {
        backwardPatternM.exportIdx(backwardM); backwardM = 0.0;
        backwardPatternP.exportIdx(backwardP); backwardP = 0.0;
    }
 
    for (const auto& is : intersections(glue))
    {
        const auto& targetElement = is.targetEntity(0);
        const auto& targetElementGeometry = targetElement.geometry();
 
        targetLocalView.bind( targetElement );
        const auto& targetLocalBasis = targetLocalView.tree().finiteElement().localBasis();
 
        // perform integration over intersection geometry
        using IsGeometry = typename std::decay_t<decltype(is.geometry())>;
        using ctype = typename IsGeometry::ctype;
 
        const auto& isGeometry = is.geometry();
        const int intOrder = maxBasisOrder + 1;
        const auto& quad = Dune::QuadratureRules<ctype, IsGeometry::mydimension>::rule(isGeometry.type(), intOrder);
        for (auto&& qp : quad)
        {
            const auto weight = qp.weight();
            const auto ie = isGeometry.integrationElement(qp.position());
            const auto globalPos = isGeometry.global(qp.position());
 
            std::vector< Dune::FieldVector<ctype, 1> > targetShapeVals;
            targetLocalBasis.evaluateFunction(targetElementGeometry.local(globalPos), targetShapeVals);
 
            // mass matrix entries target domain
            for (unsigned int i = 0; i < targetLocalBasis.size(); ++i)
            {
                const auto dofIdxI = targetLocalView.index(i);
                forwardM[dofIdxI][dofIdxI][0][0] += ie*weight*targetShapeVals[i]*targetShapeVals[i];
 
                for (unsigned int j = i+1; j < targetLocalBasis.size(); ++j)
                {
                    const auto dofIdxJ = targetLocalView.index(j);
                    const auto value = ie*weight*targetShapeVals[i]*targetShapeVals[j];
                    forwardM[dofIdxI][dofIdxJ][0][0] += value;
                    forwardM[dofIdxJ][dofIdxI][0][0] += value;
                }
            }
 
            // If targetDim < domainDim, there can be several "neighbors" if
            // targetElement is aligned with a facet of domainElement. In
            // this case make sure the basis functions are not added
            // multiple times! (division by numNeighbors)
            const auto numNeighbors = is.numDomainNeighbors();
            for (unsigned int nIdx = 0; nIdx < numNeighbors; ++nIdx)
            {
                const auto& domainElement = is.domainEntity(nIdx);
                domainLocalView.bind( domainElement );
                const auto& domainLocalBasis = domainLocalView.tree().finiteElement().localBasis();
 
                std::vector< Dune::FieldVector<ctype, 1> > domainShapeVals;
                domainLocalBasis.evaluateFunction(domainElement.geometry().local(globalPos), domainShapeVals);
 
                // add entries in matrices
                for (unsigned int i = 0; i < domainLocalBasis.size(); ++i)
                {
                    const auto dofIdxDomain = domainLocalView.index(i);
                    const auto domainShapeVal = domainShapeVals[i];
                    if (doBidirectional)
                    {
                        backwardM[dofIdxDomain][dofIdxDomain][0][0] += ie*weight*domainShapeVal*domainShapeVal;
 
                        for (unsigned int j = i+1; j < domainLocalBasis.size(); ++j)
                        {
                            const auto dofIdxDomainJ = domainLocalView.index(j);
                            const auto value = ie*weight*domainShapeVal*domainShapeVals[j];
                            backwardM[dofIdxDomain][dofIdxDomainJ][0][0] += value;
                            backwardM[dofIdxDomainJ][dofIdxDomain][0][0] += value;
                        }
                    }
 
                    for (unsigned int j = 0; j < targetLocalBasis.size(); ++j)
                    {
                        const auto dofIdxTarget = targetLocalView.index(j);
                        const auto entry = ie*weight*domainShapeVal*targetShapeVals[j];
 
                        forwardP[dofIdxTarget][dofIdxDomain][0][0] += entry/numNeighbors;
                        if (doBidirectional)
                            backwardP[dofIdxDomain][dofIdxTarget][0][0] += entry;
                    }
                }
            }
        }
    }
 
    // maybe treat zeroes on the diagonal
    if (treatDiagonalZeroes)
    {
        for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget)
            if (forwardM[dofIdxTarget][dofIdxTarget][0][0] == 0.0)
                forwardM[dofIdxTarget][dofIdxTarget][0][0] = 1.0;
 
        if (doBidirectional)
        {
            for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain)
                if (backwardM[dofIdxDomain][dofIdxDomain][0][0] == 0.0)
                    backwardM[dofIdxDomain][dofIdxDomain][0][0] = 1.0;
        }
    }
 
    return std::make_pair( std::make_pair(std::move(forwardM), std::move(forwardP)),
                           std::make_pair(std::move(backwardM), std::move(backwardP)) );
}

template<bool doBidirectional, class FEBasisDomain, class FEBasisTarget, class GlueType>
auto makeProjectorPair(const FEBasisDomain& feBasisDomain,
                       const FEBasisTarget& feBasisTarget,
                       const GlueType& glue)
{
    using ForwardProjector = typename ProjectorTraits<FEBasisDomain, FEBasisTarget>::Projector;
    using BackwardProjector = typename ProjectorTraits<FEBasisTarget, FEBasisDomain>::Projector;
 
    using ForwardProjectionMatrix = typename ForwardProjector::Matrix;
    using BackwardProjectionMatrix = typename BackwardProjector::Matrix;
 
    auto projectionMatrices = createProjectionMatrices<doBidirectional>(feBasisDomain, feBasisTarget, glue, false);
    auto& forwardMatrices = projectionMatrices.first;
    auto& backwardMatrices = projectionMatrices.second;
 
    auto& forwardM = forwardMatrices.first;
    auto& forwardP = forwardMatrices.second;
 
    auto& backwardM = backwardMatrices.first;
    auto& backwardP = backwardMatrices.second;
 
    // determine the dofs that do not take part in intersections
    std::vector<bool> isVoidTarget(forwardM.N(), false);
    for (std::size_t dofIdxTarget = 0; dofIdxTarget < forwardM.N(); ++dofIdxTarget)
        if (forwardM[dofIdxTarget][dofIdxTarget][0][0] == 0.0)
            isVoidTarget[dofIdxTarget] = true;
 
    std::vector<bool> isVoidDomain;
    if (doBidirectional)
    {
        isVoidDomain.resize(backwardM.N(), false);
        for (std::size_t dofIdxDomain = 0; dofIdxDomain < backwardM.N(); ++dofIdxDomain)
            if (backwardM[dofIdxDomain][dofIdxDomain][0][0] == 0.0)
                isVoidDomain[dofIdxDomain] = true;
    }
 
    const bool hasVoidTarget = std::any_of(isVoidTarget.begin(), isVoidTarget.end(), [] (bool v) { return v; });
    const bool hasVoidDomain = std::any_of(isVoidDomain.begin(), isVoidDomain.end(), [] (bool v) { return v; });
    if (!hasVoidDomain && !hasVoidTarget)
    {
        return std::make_pair(ForwardProjector(std::move(forwardM), std::move(forwardP)),
                              BackwardProjector(std::move(backwardM), std::move(backwardP)));
    }
    else if (!hasVoidDomain && hasVoidTarget)
    {
        std::vector<std::size_t> expansionMapTarget;
        ForwardProjectionMatrix forwardMReduced, forwardPReduced;
        setupReducedMatrices(forwardM, forwardP, isVoidTarget,
                             forwardMReduced, forwardPReduced, expansionMapTarget);
 
        return std::make_pair( ForwardProjector(std::move(forwardMReduced),
                                                std::move(forwardPReduced),
                                                std::move(expansionMapTarget),
                                                forwardM.N()),
                               BackwardProjector(std::move(backwardM), std::move(backwardP)) );
    }
    else if (hasVoidDomain && !hasVoidTarget)
    {
        if (doBidirectional)
        {
            std::vector<std::size_t> expansionMapDomain;
            BackwardProjectionMatrix backwardMReduced, backwardPReduced;
            setupReducedMatrices(backwardM, backwardP, isVoidDomain,
                                 backwardMReduced, backwardPReduced, expansionMapDomain);
 
            return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)),
                                   BackwardProjector(std::move(backwardMReduced),
                                                     std::move(backwardPReduced),
                                                     std::move(expansionMapDomain),
                                                     backwardM.N()) );
        }
        else
            return std::make_pair( ForwardProjector(std::move(forwardM), std::move(forwardP)),
                                   BackwardProjector(std::move(backwardM), std::move(backwardP)) );
    }
    else
    {
        std::vector<std::size_t> expansionMapTarget;
        ForwardProjectionMatrix forwardMReduced, forwardPReduced;
        setupReducedMatrices(forwardM, forwardP, isVoidTarget,
                             forwardMReduced, forwardPReduced, expansionMapTarget);
 
        if (doBidirectional)
        {
            std::vector<std::size_t> expansionMapDomain;
            BackwardProjectionMatrix backwardMReduced, backwardPReduced;
            setupReducedMatrices(backwardM, backwardP, isVoidDomain,
                                 backwardMReduced, backwardPReduced, expansionMapDomain);
 
            return std::make_pair( ForwardProjector(std::move(forwardMReduced),
                                                    std::move(forwardPReduced),
                                                    std::move(expansionMapTarget),
                                                    forwardM.N()),
                                   BackwardProjector(std::move(backwardMReduced),
                                                     std::move(backwardPReduced),
                                                     std::move(expansionMapDomain),
                                                     backwardM.N()) );
        }
        else
            return std::make_pair( ForwardProjector(std::move(forwardMReduced),
                                                    std::move(forwardPReduced),
                                                    std::move(expansionMapTarget),
                                                    forwardM.N()),
                                   BackwardProjector(std::move(backwardM), std::move(backwardP)) );
    }
}
 
} // end namespace Detail
 

template< class FEBasisDomain, class FEBasisTarget, class GlueType >
auto makeProjectorPair(const FEBasisDomain& feBasisDomain,
                       const FEBasisTarget& feBasisTarget,
                       GlueType glue)
{
    // we assume that target dim <= domain dimension
    static constexpr int domainDim = FEBasisDomain::GridView::dimension;
    static constexpr int targetDim = FEBasisTarget::GridView::dimension;
    static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments");
 
    return Detail::makeProjectorPair<true>(feBasisDomain, feBasisTarget, glue);
}

template< class FEBasisDomain, class FEBasisTarget, class GlueType >
auto makeProjector(const FEBasisDomain& feBasisDomain,
                   const FEBasisTarget& feBasisTarget,
                   GlueType glue)
{
    // we assume that target dim <= domain dimension
    static constexpr int domainDim = FEBasisDomain::GridView::dimension;
    static constexpr int targetDim = FEBasisTarget::GridView::dimension;
    static_assert(targetDim <= domainDim, "makeProjectorPair() expects targetDim < domainDim, please swap arguments");
 
    return Detail::makeProjectorPair<false>(feBasisDomain, feBasisTarget, glue).first;
}

template< class FEBasisDomain, class FEBasisTarget, class GlueType >
auto makeProjectionMatricesPair(const FEBasisDomain& feBasisDomain,
                                const FEBasisTarget& feBasisTarget,
                                GlueType glue)
{
    // we assume that target dim <= domain dimension
    static constexpr int domainDim = FEBasisDomain::GridView::dimension;
    static constexpr int targetDim = FEBasisTarget::GridView::dimension;
    static_assert(targetDim <= domainDim, "makeProjectionMatrixPair() expects targetDim < domainDim, please swap arguments");
 
    return Detail::createProjectionMatrices<true>(feBasisDomain, feBasisTarget, glue);
}

template< class FEBasisDomain, class FEBasisTarget, class GlueType >
auto makeProjectionMatrices(const FEBasisDomain& feBasisDomain,
                            const FEBasisTarget& feBasisTarget,
                            GlueType glue)
{
    // we assume that target dim <= domain dimension
    static constexpr int domainDim = FEBasisDomain::GridView::dimension;
    static constexpr int targetDim = FEBasisTarget::GridView::dimension;
    static_assert(targetDim <= domainDim, "makeProjectionMatrixPair() expects targetDim < domainDim, please swap arguments");
 
    return Detail::createProjectionMatrices<false>(feBasisDomain, feBasisTarget, glue).first;
}
 
} // end namespace Dumux
 
#endif