vertexmapper.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
 
#ifndef DUMUX_ENRICHED_VERTEX_DOF_MAPPER_HH
#define DUMUX_ENRICHED_VERTEX_DOF_MAPPER_HH
 
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/timer.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
#include "enrichmenthelper.hh"
 
namespace Dumux {
 
class EnrichmentIndicator
{
 
public:
    template< class GridView,
              class VertexMapper,
              class CodimOneGridView,
              class CodimOneGridAdapter >
    static void markVerticesForEnrichment(std::vector<bool>& vertexMarkers,
                                          const GridView& gridView,
                                          const VertexMapper& vertexMapper,
                                          const CodimOneGridView& codimOneGridView,
                                          const CodimOneGridAdapter& codimOneGridAdapter)
    {
        static constexpr int dim = GridView::dimension;
        static_assert(CodimOneGridView::dimension == dim-1, "Grid dimension mismatch");
 
        // reset the markers
        vertexMarkers.assign(gridView.size(dim), false);
 
        // first find all bulk grid vertices on the boundary
        std::vector<bool> isOnBoundary(gridView.size(dim), false);
        for (const auto& e : elements(gridView))
        {
            const auto refElem = referenceElement(e);
            for (const auto& is : intersections(gridView, e))
                if (is.boundary())
                    for (int i = 0; i < is.geometry().corners(); ++i)
                        isOnBoundary[ vertexMapper.subIndex( e,
                                                             refElem.subEntity(is.indexInInside(), 1, i, dim),
                                                             dim ) ] = true;
        }
 
        // mark all vertices on the lower-dimensional grid to be enriched (except immersed boundaries)
        std::vector<typename GridView::IndexSet::IndexType> vertexIndicesStorage;
        for (const auto& codimOneElement : elements(codimOneGridView))
        {
            // if a codimension one element has 2 or more embedments, we need to enrich
            if (codimOneGridAdapter.numEmbedments(codimOneElement) >= 2)
            {
                for (int i = 0; i < codimOneElement.subEntities(dim-1); ++i)
                    vertexMarkers[ codimOneGridAdapter.bulkGridVertexIndex(codimOneElement.template subEntity<dim-1>(i)) ] = true;
 
                // however, we have to exclude immersed boundaries
                const auto refElem = referenceElement(codimOneElement);
                for (const auto& intersection : intersections(codimOneGridView, codimOneElement))
                {
                    // skip if intersection is not on boundary
                    if (!intersection.boundary())
                        continue;
 
                    // obtain all bulk grid indices of the lower-dimensional intersection corners
                    const auto numCorners = intersection.geometry().corners();
                    vertexIndicesStorage.resize(numCorners);
                    for (int i = 0; i < numCorners; ++i)
                    {
                        const auto vIdxLocal = refElem.subEntity(intersection.indexInInside(), 1, i, dim-1);
                        vertexIndicesStorage[i] = codimOneGridAdapter.bulkGridVertexIndex( codimOneElement.template subEntity<dim-1>(vIdxLocal) );
                    }
 
                    // if any of the vertices is on an immersed boundary, we must not enrich any of them
                    if (std::any_of(vertexIndicesStorage.begin(), vertexIndicesStorage.end(), [&isOnBoundary] (auto idx) { return !isOnBoundary[idx]; }))
                        std::for_each(vertexIndicesStorage.begin(), vertexIndicesStorage.end(), [&vertexMarkers] (auto idx) { vertexMarkers[idx] = false; });
                }
            }
        }
    }
};
 
template<class GV>
class EnrichedVertexDofMapper
{
    static constexpr int dim = GV::dimension;
    static_assert(dim > 1, "Vertex dof enrichment mapper currently only works for dim > 1!");
 
    using GIType = typename GV::IndexSet::IndexType;
    using Vertex = typename GV::template Codim<dim>::Entity;
    using Element = typename GV::template Codim<0>::Entity;
    using MCMGMapper = Dune::MultipleCodimMultipleGeomTypeMapper<GV>;
 
public:
    using GridView = GV;
    using GridIndexType = GIType;
 
    EnrichedVertexDofMapper(const GV& gridView)
    : gridView_(gridView)
    , elementMapper_(gridView, Dune::mcmgElementLayout())
    , vertexMapper_(gridView, Dune::mcmgVertexLayout())
    {
        initialize_();
    }
 
    EnrichedVertexDofMapper(const GV& gridView, Dune::MCMGLayout layout)
    : EnrichedVertexDofMapper(gridView)
    {
        if ( !( static_cast<bool>(layout(Dune::GeometryTypes::vertex, dim)) ) )
            DUNE_THROW(Dune::InvalidStateException, "Vertex mapper only makes sense for vertex layout!");
    }
 
    GridIndexType subIndex(const Element& e, unsigned int i, unsigned int codim) const
    {
        assert(codim == dim && "Only element corners can be mapped by this mapper");
        return indexMap_[elementMapper_.index(e)][i];
    }
 
    GridIndexType vertexIndex(const Element& e, unsigned int i, unsigned int codim) const
    {
        assert(codim == dim && "Only element corners can be mapped by this mapper");
        return vertexMapper_.subIndex(e, i, codim);
    }
 
    GridIndexType vertexIndex(const Vertex& v) const
    {
        assert(Vertex::Geometry::mydimension == 0 && "Only vertices can be mapped by this mapper");
        return vertexMapper_.index(v);
    }
 
    template< class EntityType >
    GridIndexType index(const EntityType& e) const
    {
        if (hasEnrichedVertices_)
            DUNE_THROW(Dune::InvalidStateException, "Index map contains enriched vertex dofs. Direct mapping from vertex to index not possible.");
 
        assert(EntityType::Geometry::mydimension == 0 && "Only vertices can be mapped by this mapper");
        return vertexMapper_.index(e);
    }
 
    std::size_t size() const
    { return size_; }
 
    bool isEnriched(const Vertex& v)
    { return isEnriched_[ vertexMapper_.index(v) ]; }
 
    void update(const GV& gridView)
    {
        gridView_ = gridView;
        initialize_();
    }
 
    void update(GV&& gridView)
    {
        gridView_ = std::move(gridView);
        initialize_();
    }
 
    template<class CodimOneGridView, class CodimOneGridAdapter>
    void enrich(const CodimOneGridView& codimOneGridView,
                const CodimOneGridAdapter& codimOneGridAdapter,
                bool verbose = false)
    {
        static const int codimOneDim = CodimOneGridView::dimension;
        static_assert(codimOneDim == dim-1, "Grid dimension mismatch!");
        static_assert(codimOneDim == 2 || codimOneDim == 1, "Inadmissible codimension one grid dimension");
        static_assert(int(CodimOneGridView::dimensionworld) == int(GV::dimensionworld), "Grid world dimension mismatch");
 
        // keep track of time
        Dune::Timer watch;
 
        // mark vertices for enrichment using the indicator
        EnrichmentIndicator::markVerticesForEnrichment(isEnriched_, gridView_, vertexMapper_, codimOneGridView, codimOneGridAdapter);
 
        // let the helper class do the enrichment of the index map
        size_ = VertexEnrichmentHelper< GridView, CodimOneGridView >::enrich(indexMap_,
                                                                             isEnriched_,
                                                                             gridView_,
                                                                             vertexMapper_,
                                                                             elementMapper_,
                                                                             codimOneGridView,
                                                                             codimOneGridAdapter);
 
        // check if new index map contains enriched dofs
        hasEnrichedVertices_ = std::any_of(isEnriched_.begin(), isEnriched_.end(), [] (bool isEnriched) { return isEnriched; });
 
        if (verbose)
            std::cout << "Vertex dof enrichment took " << watch.elapsed() << " seconds." << std::endl;
    }
 
private:
 
    void initialize_()
    {
        size_ = gridView_.size(dim);
        hasEnrichedVertices_ = false;
        indexMap_.resize(gridView_.size(0));
        isEnriched_.resize(gridView_.size(dim), false);
        for (const auto& e : elements(gridView_))
        {
            const auto numCorners = e.geometry().corners();
            const auto eIdxGlobal = elementMapper_.index(e);
            indexMap_[eIdxGlobal].resize(numCorners);
            for (unsigned int i = 0; i < numCorners; ++i)
                indexMap_[eIdxGlobal][i] = vertexMapper_.subIndex(e, i, dim);
        }
    }
 
    // data members
    std::size_t size_;                        
    GV gridView_;                             
    MCMGMapper elementMapper_;                
    MCMGMapper vertexMapper_;                 
    bool hasEnrichedVertices_;                
    std::vector<bool> isEnriched_;            
    std::vector< std::vector<GIType> > indexMap_; 
};
 
} // end namespace Dumux
 
#endif