errorfractionmarking.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 ERROR_FRACTION_MARKING_HH
#define ERROR_FRACTION_MARKING_HH
 
#include <vector>
 
#include <dune/common/bitsetvector.hh>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/rangegenerators.hh>
 
#include <dune/fufem/estimators/refinementindicator.hh>
 
template <class GridType>
class ErrorFractionMarkingStrategy
{
    private:
 
        enum {dim = GridType::dimension};
 
        struct AllCodimLayout {
            bool operator() ([[maybe_unused]] Dune::GeometryType gt, int) const {
                return true;
            }
        };
 
    public:
        static void mark(RefinementIndicator<GridType>& refinementIndicator,
                         GridType& grid,
                         double fraction);
 
        static void mark(const std::vector<RefinementIndicator<GridType>* >& refinementIndicator,
                         const std::vector<GridType*>& grids,
                         double fraction);
 
};
 
template <class GridType>
void ErrorFractionMarkingStrategy<GridType>::
mark(RefinementIndicator<GridType>& refinementIndicator,
     GridType& grid,
     double fraction)
{
    std::vector<RefinementIndicator<GridType>* > refinementIndicatorInVector(1);
    refinementIndicatorInVector[0] = &refinementIndicator;
 
    std::vector<GridType*> gridInVector(1);
    gridInVector[0] = &grid;
 
    mark(refinementIndicatorInVector, gridInVector, fraction);
}
 
template <class GridType>
void ErrorFractionMarkingStrategy<GridType>::
mark(const std::vector<RefinementIndicator<GridType>*>& refinementIndicators,
     const std::vector<GridType*>& grids,
     double fraction)
{
    typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> AllCodimMapper;
    typedef typename std::multimap<double, std::pair<int,int>, std::greater<double> > ErrorMap;
 
    if (grids.size() != refinementIndicators.size())
        DUNE_THROW(Dune::Exception, "The number of grids must match the number of error sets");
 
    // /////////////////////////////////////////////////////////////////
    //   Create a map which stores the error for all elements of
    //   all grids and which can be traversed in sorted order.
    // /////////////////////////////////////////////////////////////////
    ErrorMap errorMap;
 
    std::vector<Dune::BitSetVector<1> > refineIndex(grids.size());
 
    // we need to sum the error indicators up
    double totalError = 0.0;
 
    for (size_t i=0; i<grids.size(); ++i)
    {
        // mapper that maps all entities to indices
        AllCodimMapper mapper(*grids[i], AllCodimLayout());
 
        // initialize bitfield
        refineIndex[i].resize(mapper.size(), false);
 
        // abuse bitfield for processed flags
        Dune::BitSetVector<1>& entityProcessed = refineIndex[i];
 
        for (const auto& e : elements(grids[i]->leafGridView()))
        {
            // Get set of shape functions
            auto refElement
                = Dune::ReferenceElements<double, dim>::general(e.type());
 
            // Compute the indicator of any subentity of this element
            for (int codim=0; codim<=dim; ++codim)
            {
                for (int entity=0; entity<refElement.size(codim); ++entity)
                {
                    int entityIndex = mapper.subIndex(e, entity, codim);
                    if (not(entityProcessed[entityIndex][0]))
                    {
                        // map error indicator values to entity index and sum them up
                        double errorIndicator = (*refinementIndicators[i]).value(e, codim, entity);
                        totalError += errorIndicator;
                        errorMap.insert(std::pair<double, std::pair<int,int> >(errorIndicator, std::pair<int,int>(i, entityIndex)));
                        entityProcessed[entityIndex][0] = true;
                    }
                }
            }
        }
 
        // reset bitfields
        entityProcessed.unsetAll();
    }
 
    // mark indicators by value fraction
    {
        double refinedError = 0.0;
        for (const auto& e : errorMap)
        {
            refineIndex[e.second.first][e.second.second][0] = true;
            refinedError += e.first;
            if (refinedError >= totalError*fraction)
              break;
        }
    }
 
    for (size_t i=0; i<grids.size(); ++i)
    {
        AllCodimMapper mapper(*grids[i], AllCodimLayout());
 
        for (const auto& e : elements(grids[i]->leafGridView()))
        {
            // Get set of shape functions
            auto refElement
                = Dune::ReferenceElements<double, dim>::general(e.type());
 
            // Determine if the element is marked to refinement by any of its subentities
            bool refineElement = false;
            for (int codim=0; (codim<=dim) and not(refineElement); ++codim)
            {
                for (int entity=0; (entity<refElement.size(codim)) and not(refineElement); ++entity)
                {
                    int entityIndex = mapper.subIndex(e, entity, codim);
                    if (refineIndex[i][entityIndex][0])
                        refineElement = true;
                }
            }
            if (refineElement)
                grids[i]->mark(1,e);
        }
    }
}
#endif