transferoperatorassembler.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// 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 DUNE_FUFEM_ASSEMBLERS_TRANSFER_OPERATOR_ASSEMBLER_HH
#define DUNE_FUFEM_ASSEMBLERS_TRANSFER_OPERATOR_ASSEMBLER_HH
 
#include <vector>
#include <map>
#include <memory>
 
#include <dune/istl/matrixindexset.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/timer.hh>
#include <dune/geometry/referenceelements.hh>
 
#include <dune/functions/functionspacebases/lagrangebasis.hh>
 
#include <dune/fufem/utilities/asnumber.hh>
 
namespace Dune::Fufem::Impl {
 
template <class GridType>
class MultiLevelPQ1Indices
{
    private:
        static const int dim = GridType::dimension;
 
    public:
 
        MultiLevelPQ1Indices(const GridType& grid) :
            grid(grid)
        {
            const auto& localIdSet = grid.localIdSet();
            const auto& leafIndexSet = grid.leafIndexSet();
 
            int maxLevel = grid.maxLevel();
 
            idToIndex.resize(maxLevel+1);
            size_.resize(maxLevel+1);
 
            size_[maxLevel] = grid.size(dim);
 
            // build extended level indices as map (localId -> index)
            // 1. enumerate vertices on each level
            // 2. enumerate leaf vertices on each finer level to emulate copies of them
            // 3. on maxlevel take the leaf indices
            //
            // example: level=*, copy=(*), leaf=[*]
            // 
            //  [1] -  -  -  -  -  -  - [2]---[3]---[4] -  -  - [5]
            //   
            //  (4) -  -  -  -  -  -  -  1-----------2-----------3
            // 
            //   1-----------------------2-----------------------3
 
            // first use the level indices
            for(int level=0; level<maxLevel; ++level)
            {
                const auto& indexSet = grid.levelIndexSet(level);
                size_[level] = indexSet.size(dim);
                for (const auto& it : vertices(grid.levelGridView(level)))
                    idToIndex[level][localIdSet.id(it)] = indexSet.index(it);
            }
 
            for (const auto& it : vertices(grid.leafGridView()))
            {
                idToIndex[maxLevel][localIdSet.id(it)] = leafIndexSet.index(it);
 
                for(int level=it.level()+1; level<maxLevel; ++level)
                {
                    idToIndex[level][localIdSet.id(it)] = size_[level];
                    ++size_[level];
                }
            }
        }
 
        int size(int level)
        {
            return size_[level];
        }
 
        template<class LocalView>
        int index(const LocalView& localView, const int i, const int level) const
        {
            const Dune::LocalKey& localKey = localView.tree().finiteElement().localCoefficients().localKey(i);
            const IdType id = grid.localIdSet().subId(localView.element(), localKey.subEntity(), localKey.codim());
            return idToIndex[level].find(id)->second;
        }
 
 
    private:
        const GridType& grid;
 
        typedef typename GridType::Traits::LocalIdSet::IdType IdType;
        std::vector< std::map<IdType,int> > idToIndex;
        std::vector<int> size_;
};
 
} // namespace Dune::Fufem::Impl
 
 
template <class GridType>
class TransferOperatorAssembler {
 
    public:
        TransferOperatorAssembler(const GridType& grid) :
            grid(grid)
        {}
 
 
        template <class TransferOperator>
        void assembleOperatorHierarchy(std::vector<TransferOperator>& T) const
        {
            std::vector <typename TransferOperator::TransferOperatorType*> M;
 
            int maxLevel = grid.maxLevel();
 
            M.resize(maxLevel);
 
            for (int i=0; i<maxLevel; ++i)
                M[i] = &(const_cast<typename TransferOperator::TransferOperatorType&>(T[i].getMatrix()));
            assembleMatrixHierarchy(M);
        }
 
        template <class TransferOperator, class RealTransferOperator>
        void assembleDerivedOperatorPointerHierarchy(std::vector<TransferOperator*>& T) const
        {
            typedef typename RealTransferOperator::TransferOperatorType TransferOperatorType;
            std::vector <TransferOperatorType*> M;
 
            int maxLevel = grid.maxLevel();
 
            M.resize(maxLevel);
 
            for (int i=0; i<maxLevel; ++i)
            {
                RealTransferOperator* t = dynamic_cast<RealTransferOperator*>(T[i]);
                M[i] = &(const_cast<TransferOperatorType&>(t->getMatrix()));
            }
            assembleMatrixHierarchy(M);
        }
 
        template <class TransferOperator>
        void assembleOperatorPointerHierarchy(std::vector<TransferOperator*>& T) const
        {
            std::vector <typename TransferOperator::TransferOperatorType*> M;
 
            int maxLevel = grid.maxLevel();
 
            M.resize(maxLevel);
 
            for (int i=0; i<maxLevel; ++i)
                M[i] = &(const_cast<typename TransferOperator::TransferOperatorType&>(T[i]->getMatrix()));
            assembleMatrixHierarchy(M);
        }
 
        template <class Matrix>
        void assembleMatrixHierarchy(std::vector<std::shared_ptr<Matrix> >& M) const
        {
            while(M.size() < uint(grid.maxLevel()))
                M.push_back(std::make_shared<Matrix>());
 
            std::vector <Matrix*> Mraw(M.size());
            for (size_t i=0; i<M.size(); ++i)
                Mraw[i] = M[i].get();
 
            assembleMatrixHierarchy(Mraw);
        }
 
        template <class Matrix>
        void assembleMatrixHierarchy(std::vector<Matrix*>& T) const
        {
            typedef std::map<int, double> LinearCombination;
            typedef std::vector<LinearCombination> BaseTransformation;
            typedef std::vector<BaseTransformation> TransformationHierarchy;
            typedef Dune::Functions::LagrangeBasis<typename GridType::LevelGridView, 1> LevelBasis;
            typedef typename LevelBasis::LocalView::Tree::FiniteElement LFE;
            typedef typename LFE::Traits::LocalBasisType::Traits::RangeType FERange;
 
            Dune::Timer timer;
 
            int maxLevel = grid.maxLevel();
 
            std::vector<LevelBasis> levelBasis;
            for(int level=0; level<=maxLevel; ++level)
                levelBasis.emplace_back(grid.levelGridView(level));
 
            Dune::Fufem::Impl::MultiLevelPQ1Indices<GridType> multiLevelIndexSet(grid);
 
#ifdef FE_VERBOSE
            std::cout << "FE:" << "localId -> index maps build in " << timer.elapsed() << " seconds." << std::endl;
#endif
 
 
            TransformationHierarchy transformationHierarchy(maxLevel);
            for (int level=0; level<maxLevel; ++level)
                transformationHierarchy[level].resize(multiLevelIndexSet.size(level+1));
 
            // set all nodes as not processed
            std::vector< Dune::BitSetVector<1> > processed(maxLevel+1);
            for (int level=0; level<=maxLevel; ++level)
                processed[level].resize(multiLevelIndexSet.size(level), false);
 
            // loop over all levels
            timer.reset();
            for(int level=0; level<grid.maxLevel(); ++level)
            {
                auto coarseLocalView = levelBasis[level].localView();
                auto fineLocalView = levelBasis[level+1].localView();
 
                // loop over all elements of current level
                for (const auto& cIt : elements(grid.levelGridView(level)))
                {
                    coarseLocalView.bind(cIt);
                    const LFE& coarseFE = coarseLocalView.tree().finiteElement();
 
                    // if element is leaf the transfer to the next level is locally the identity
                    if (cIt.isLeaf())
                    {
                        for (int coarseLevel=level; coarseLevel<maxLevel; ++coarseLevel)
                        {
                            int fineLevel = coarseLevel+1;
 
                            for(size_t j=0; j<coarseFE.localBasis().size(); ++j)
                            {
                                int fineIndex = multiLevelIndexSet.index(coarseLocalView, j, fineLevel);
 
                                // visit each child node only once
                                if (processed[fineLevel][fineIndex][0])
                                    continue;
 
                                int coarseIndex = multiLevelIndexSet.index(coarseLocalView, j, coarseLevel);
                                transformationHierarchy[coarseLevel][fineIndex][coarseIndex] = 1.0;
                                processed[fineLevel][fineIndex][0] = true;
                            }
                        }
                    }
                    else
                    {
                        int coarseLevel = level;
                        int fineLevel = level+1;
 
                        // store coarse node indices since we need them often
                        std::vector<int> coarseIndex(coarseFE.localBasis().size());
                        for(size_t i=0; i<coarseFE.localBasis().size(); ++i)
                            coarseIndex[i] = multiLevelIndexSet.index(coarseLocalView, i, coarseLevel);
 
                        std::vector<FERange> valuesAtPosition(coarseFE.localBasis().size());
 
                        // loop over all children on next level
                        for (const auto& fIt : descendantElements(cIt,level+1))
                        {
                            fineLocalView.bind(fIt);
                            const LFE& fineFE = fineLocalView.tree().finiteElement();
 
                            // we need the reference element to get the local position of the subentities corresponding to fine basis functions
                            auto fineRefElement = Dune::ReferenceElements<double, dim>::general(fIt.type());
 
                            // loop over all child nodes
                            for(size_t j=0; j<fineFE.localBasis().size(); ++j)
                            {
                                int fineIndex = multiLevelIndexSet.index(fineLocalView, j, fineLevel);
 
                                // visit each child node only once
                                if (processed[fineLevel][fineIndex][0])
                                    continue;
 
                                // get local coordinates of subentity in fine element
                                const Dune::LocalKey& localKey = fineFE.localCoefficients().localKey(j);
                                auto localPositionFine = fineRefElement.position(localKey.subEntity(), localKey.codim());
 
                                // compute local coordinates of subentity in coarse element
                                auto localPositionCoarse = fIt.geometryInFather().global(localPositionFine);
 
                                // evaluate coarse basis functions at the position of the subentity corresponding to the fine basis function
                                coarseFE.localBasis().evaluateFunction(localPositionCoarse, valuesAtPosition);
 
                                for(size_t i=0; i<coarseFE.localBasis().size(); ++i)
                                {
                                    if (Dune::Fufem::asNumber(valuesAtPosition[i]) > 1e-5)
                                        transformationHierarchy[coarseLevel][fineIndex][coarseIndex[i]] = Dune::Fufem::asNumber(valuesAtPosition[i]);
                                }
 
                                processed[fineLevel][fineIndex][0] = true;
                            } // loop over all child nodes
                        } // end of loop over all children on next level
                    }
                } // end of loop over all elements of current level
            } // end of loop over all levels
#ifdef FE_VERBOSE
            std::cout << "FE:" << "Grid traversed for transfer operators in " << timer.elapsed() << " seconds." << std::endl; 
#endif
 
            // setup transfer operator matrices
            timer.reset();
            //  T.resize(maxLevel);
            for(int level=0; level<maxLevel; ++level)
            {
                Dune::MatrixIndexSet indices(multiLevelIndexSet.size(level+1), multiLevelIndexSet.size(level));
 
                for(size_t row=0; row<transformationHierarchy[level].size(); ++row)
                    for(const auto& colIt : transformationHierarchy[level][row])
                        indices.add(row,colIt.first);
 
                indices.exportIdx(*T[level]);
 
                for(size_t row=0; row<transformationHierarchy[level].size(); ++row)
                {
                    for(const auto& colIt : transformationHierarchy[level][row])
                    {
                        //                if (colIt->second != 0.0)
                        (*T[level])[row][colIt.first] = 0.0;
                        for(size_t i=0; i<Matrix::block_type::rows; ++i)
                            (*T[level])[row][colIt.first][i][i] = colIt.second;
                    }
                }
 
                transformationHierarchy[level].clear();
            }
 
#ifdef FE_VERBOSE
            for(int level=0; level<maxLevel; ++level)
                std::cout << "FE:" << "Transfer " << level << " -> " << level+1 << " is a "<< (*T[level]).N() << " x " << (*T[level]).M() << " matrix." << std::endl;
            std::cout << "FE:" << "Transfer operator matrices set up in " << timer.elapsed() << " seconds." << std::endl;
#endif
        }
 
 
    private:
        static const int dim = GridType::dimension;
        const GridType& grid;
};
 
#endif