fractionalmarking.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 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_FRACTIONAL_MARKING_HH
#define DUNE_FUFEM_FRACTIONAL_MARKING_HH
 
#include <vector>
#include <limits>
 
#include <dune/common/bitsetvector.hh>
#include <dune/common/shared_ptr.hh>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
 
#include <dune/fufem/estimators/refinementindicator.hh>
 
template <class GridType, class field_type = double>
class FractionalMarkingStrategy
{
    enum {dim = GridType::dimension};
 
    template<int localDim>
    struct P0Layout
    {
        bool contains (Dune::GeometryType gt) {
            return gt.dim() == localDim;
        }
    };
 
public:
    static void mark(RefinementIndicator<GridType>& refinementIndicator,
                     GridType& grid,
                     field_type fraction);
 
    // TODO: IMHO the implementation does not provide a mechanism that prevents
    // refining elements where all associated errors are zero.
    static void mark(const std::vector<std::shared_ptr<RefinementIndicator<GridType> > >& refinementIndicator,
                     const std::vector<GridType*>& grids,
                     field_type fraction);
 
};
 
template <class GridType, class field_type>
void FractionalMarkingStrategy<GridType,field_type>::
mark(RefinementIndicator<GridType>& refinementIndicator,
     GridType& grid,
     field_type fraction)
{
    std::vector<std::shared_ptr<RefinementIndicator<GridType> > > refinementIndicatorInVector(1);
    refinementIndicatorInVector[0] = Dune::stackobject_to_shared_ptr<RefinementIndicator<GridType> >(refinementIndicator);
 
    std::vector<GridType*> gridInVector(1);
    gridInVector[0] = &grid;
 
    mark(refinementIndicatorInVector, gridInVector, fraction);
}
 
template <class GridType, class field_type>
void FractionalMarkingStrategy<GridType,field_type>::
mark(const std::vector<std::shared_ptr<RefinementIndicator<GridType> > >& refinementIndicators,
     const std::vector<GridType*>& grids,
     field_type fraction)
{
    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.
    // /////////////////////////////////////////////////////////////////
    std::multimap<field_type, std::pair<int,int> > errorMap;
    using Iterator = typename std::multimap<field_type, std::pair<int,int> >::reverse_iterator;
 
    std::vector<field_type> minGridError(grids.size()), maxGridError(grids.size());
 
    for (size_t i=0; i<grids.size(); i++) {
 
        minGridError[i] =  std::numeric_limits<field_type>::max();
        maxGridError[i] = -std::numeric_limits<field_type>::max();
 
        Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> p0Mapper(grids[i]->leafGridView(), Dune::mcmgElementLayout());
 
        const auto& leafView = grids[i]->leafGridView();
 
        for (const auto& e : elements(leafView)) {
 
            // Get set of shape functions
            const auto& refElement = Dune::ReferenceElements<field_type, dim>::general(e.type());
 
            // Compute the maximum indicator of any subentity of this element
            field_type maxError = 0;
            for (int codim=0; codim<=dim; codim++)
                for (int entity=0; entity<refElement.size(codim); entity++)
                     maxError = std::max(maxError, (*refinementIndicators[i]).value(e, codim, entity));
 
            // Store the max error in a map ordered by the error
            errorMap.insert(std::pair<field_type, std::pair<int,int> >(maxError, std::make_pair(i, p0Mapper.index(e))));
 
            minGridError[i] = std::min(minGridError[i], maxError);
            maxGridError[i] = std::max(maxGridError[i], maxError);
 
        }
 
    }
 
    std::cout << "Errors per grid:" << std::endl;
    for (size_t i=0; i<grids.size(); i++)
        std::cout << i << ":  " << minGridError[i] << " -- " << maxGridError[i] << std::endl;
 
    // //////////////////////////////////////////////////////////////////////
    //   Mark the desired fraction of elements with the highest error in
    //   bitfields.  This is more memory efficient than storing EntityPointers
    //   in the map.
    // //////////////////////////////////////////////////////////////////////
 
    int numElementsToBeRefined = int(std::ceil(errorMap.size()*fraction));
 
    if (numElementsToBeRefined > int(errorMap.size()))
        DUNE_THROW(Dune::Exception, "fraction can't be more than one");
 
    std::vector<Dune::BitSetVector<1> > refinedElements(grids.size());
    for (size_t i=0; i<grids.size(); i++) {
        refinedElements[i].resize(grids[i]->leafIndexSet().size(0));
        refinedElements[i].unsetAll();
    }
 
    Iterator it    = errorMap.rbegin();
 
    for (int i=0; i<numElementsToBeRefined; i++, ++it) {
        int grid    = it->second.first;
        int element = it->second.second;
        refinedElements[grid][element] = true;
    }
 
    // /////////////////////////////////////////////////////////////////////
    //   Actually mark the elements for refinement
    // /////////////////////////////////////////////////////////////////////
 
    for (size_t i=0; i<grids.size(); i++) {
 
        Dune::MultipleCodimMultipleGeomTypeMapper<typename GridType::LeafGridView> p0Mapper(grids[i]->leafGridView(), Dune::mcmgElementLayout());
 
        for (const auto& e : elements(grids[i]->leafGridView()))
            if (refinedElements[i][p0Mapper.index(e)][0])
                grids[i]->mark(1,e);
 
    }
 
 
}
#endif