gsAdaptiveRefUtils.h

Go to the documentation of this file.

#pragma once


#include <iostream>

namespace gismo
{

enum MarkingStrategy
{
    GARU=1,
    PUCA=2,
    BULK=3,
    PBULK=4
};

template <class T>
void gsMarkThreshold( const std::vector<T> & elError, T refParameter,
                      std::vector<bool> & elMarked)
{
    // First, conduct a brutal search for the maximum local error
    const T maxErr = *std::max_element(elError.begin(), elError.end() );

    // Compute the threshold:
    const T Thr = refParameter * maxErr;

    elMarked.resize( elError.size() );
    // Now just check for each element, whether the local error
    // is above the computed threshold or not, and mark accordingly.

    typename std::vector<T>::const_iterator err = elError.begin();
    for(std::vector<bool>::iterator i = elMarked.begin(); i!=  elMarked.end(); ++i, ++err)
        *i = ( *err >= Thr );
}

template <class T>
void gsMarkPercentage( const std::vector<T> & elError, T refParameter,
                       std::vector<bool> & elMarked)
{
    T Thr = (T)(0);

    // Total number of elements:
    size_t NE = elError.size();
    // The vector of local errors will need to be sorted,
    // which will be done on a copy:
    std::vector<T> elErrCopy = elError;

    // Compute the index from which the refinement should start,
    // once the vector is sorted.
    size_t idxRefineStart = cast<T,size_t>( math::floor( refParameter * (T)(NE) ) );
    // ...and just to be sure we are in range:
    if( idxRefineStart == elErrCopy.size() )
    {
        GISMO_ASSERT(idxRefineStart >= 1, "idxRefineStart can't get negative");
        idxRefineStart -= 1;
    }

    // Sort the list using bubblesort.
    // After each loop, the largest elements are at the end
    // of the list. Since we are only interested in the largest elements,
    // it is enough to run the sorting until enough "largest" elements
    // have been found, i.e., until we have reached indexRefineStart
    size_t lastSwapDone = elErrCopy.size() - 1;
    size_t lastCheckIdx = lastSwapDone;

    bool didSwap;
    T tmp;
    do{
        didSwap = false;
        lastCheckIdx = lastSwapDone;
        for( size_t i=0; i < lastCheckIdx; i++)
            if( elErrCopy[i] > elErrCopy[i+1] )
            {
                tmp = elErrCopy[i];
                elErrCopy[i] = elErrCopy[i+1];
                elErrCopy[i+1] = tmp;

                didSwap = true;
                lastSwapDone = i;
            }
    }while( didSwap && (lastSwapDone+1 >= idxRefineStart ) );

    // Compute the threshold:
    Thr = elErrCopy[ idxRefineStart ];
    elMarked.resize( elError.size() );
    // Now just check for each element, whether the local error
    // is above the computed threshold or not, and mark accordingly.
    for( size_t i=0; i < elError.size(); i++)
        ( elError[i] >= Thr ? elMarked[i] = true : elMarked[i] = false );
}


template <class T>
void gsMarkFraction( const std::vector<T> & elError, T refParameter,
                     std::vector<bool> & elMarked)
{
    T Thr = (T)(0);

    // The vector of local errors will need to be sorted,
    // which will be done on a copy:
    std::vector<T> elErrCopy = elError;

    // Compute the sum, i.e., the global/total error
    T totalError = (T)(0);
    for( size_t i = 0; i < elErrCopy.size(); ++i)
        totalError += elErrCopy[i];

    // We want to mark just enough cells such that their
    // cummulated errors add up to a certain fraction
    // of the total error.
    T errorMarkSum = (1-refParameter) * totalError;
    T cummulErrMarked = 0;

    T tmp;
    GISMO_ASSERT(elErrCopy.size() >= 1, "elErrCopy needs at least 1 element");
    size_t lastSwapDone = elErrCopy.size() - 1;
    do{
        for( size_t i=0; i < lastSwapDone; i++)
            if( elErrCopy[i] > elErrCopy[i+1] )
            {
                tmp = elErrCopy[i];
                elErrCopy[i] = elErrCopy[i+1];
                elErrCopy[i+1] = tmp;
            }

        cummulErrMarked += elErrCopy[ lastSwapDone  ];
        lastSwapDone -= 1;

    }while( cummulErrMarked < errorMarkSum && lastSwapDone > 0 );

    // Compute the threshold:
    Thr = elErrCopy[ lastSwapDone + 1 ];
    elMarked.resize( elError.size() );
    // Now just check for each element, whether the local error
    // is above the computed threshold or not, and mark accordingly.
    for( size_t i=0; i < elError.size(); i++)
        ( elError[i] >= Thr ? elMarked[i] = true : elMarked[i] = false );
}


template <class T>
void gsMarkElementsForRef( const std::vector<T> & elError, int refCriterion,
                           T refParameter, std::vector<bool> & elMarked)
{
    switch (refCriterion)
    {
    case GARU:
        gsMarkThreshold(elError,refParameter,elMarked);
        break;
    case PUCA:
        gsMarkPercentage(elError,refParameter,elMarked);
        break;
    case BULK:
        gsMarkFraction(elError,refParameter,elMarked);
        break;
    default:
        GISMO_ERROR("unknown marking strategy");
    }
}

template <class T>
std::vector<gsMatrix<T>> gsGetRefinementBoxes(  gsMultiBasis<T> & basis,
                                                const std::vector<bool> & elMarked)
{
    GISMO_ASSERT(basis.totalElements()==elMarked.size(),
                 "Incorrect input, num_el="<<basis.totalElements()
                 <<", mark_vec="<<elMarked.size() );

    std::vector<gsMatrix<T>> result(basis.nBases());
    const short_t dim = basis.dim();

    // numMarked: Number of marked cells on current patch, also currently marked cell
    // poffset  : offset index for the first element on a patch
    // globalCount: counter for the current global element index
    index_t numMarked, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes;

    for (size_t pn=0; pn < basis.nBases(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const size_t numEl = basis[pn].numElements();
        numMarked = std::count_if(elMarked.begin() + poffset,
                                  elMarked.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );
        poffset += numEl;
        refBoxes.resize(dim, 2*numMarked);
        numMarked = 0;// counting current patch element to be refined

        // for all elements in patch pn
        typename gsBasis<T>::domainIter domIt = basis.basis(pn).makeDomainIterator();
        for (; domIt->good(); domIt->next())
        {
            if( elMarked[ globalCount++ ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                refBoxes.col(2*numMarked  ) = domIt->lowerCorner();
                refBoxes.col(2*numMarked+1) = domIt->upperCorner();

                // Advance marked cells counter
                numMarked++;
            }
        }

        result[pn] = refBoxes;
    }
    return result;
}

template <class T>
void gsRefineMarkedElements(gsMultiBasis<T> & basis,
                            const std::vector<bool> & elMarked,
                            index_t refExtension = 0)
{
    GISMO_ASSERT(basis.totalElements()==elMarked.size(),
                 "Incorrect input, num_el="<<basis.totalElements()
                 <<", mark_vec="<<elMarked.size() );

    std::vector<gsMatrix<T>> boxes = gsGetRefinementBoxes(basis,elMarked);
    for (size_t pn=0; pn < basis.nBases(); ++pn )// for all patches
    {
        // Refine all of the found refBoxes in this patch
        basis.refine( pn, boxes[pn], refExtension );
    }
}

template <class T>
void gsRefineWithBoxes( gsMultiBasis<T> & basis,
                        const std::vector<gsMatrix<T>> & boxes,
                        index_t refExtension = 0)
{
    for (size_t pn=0; pn < basis.nBases(); ++pn )// for all patches
    {
        // Refine all of the found refBoxes in this patch
        basis.refine( pn, boxes[pn], refExtension );
    }
}

template <class T>
void gsRefineMarkedElements(gsMultiPatch<T> & mp,
                            const std::vector<bool> & elMarked,
                            int refExtension = 0)
{
    const int dim = mp.dim();

    // numMarked: Number of marked cells on current patch, also currently marked cell
    // poffset  : offset index for the first element on a patch
    // globalCount: counter for the current global element index
    int numMarked, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes;

    for (size_t pn=0; pn < mp.nPatches(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const int numEl = mp[pn].basis().numElements();
        numMarked = std::count_if(elMarked.begin() + poffset,
                                  elMarked.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );

        poffset += numEl;
        refBoxes.resize(dim, 2*numMarked);
        numMarked = 0;// counting current patch element to be refined

        // for all elements in patch pn
        typename gsBasis<T>::domainIter domIt = mp.patch(pn).basis().makeDomainIterator();
        for (; domIt->good(); domIt->next())
        {
            if( elMarked[ globalCount++ ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                refBoxes.col(2*numMarked  ) =
                        refBoxes.col(2*numMarked+1) = domIt->centerPoint();

                // Advance marked cells counter
                numMarked++;
            }
        }
        // Refine all of the found refBoxes in this patch
        std::vector<index_t> elements = mp.patch(pn).basis().asElements(refBoxes, refExtension);
        mp.patch(pn).refineElements( elements );
    }
}

template <class T>
void gsRefineMarkedFunctions(gsMultiPatch<T> & mp,
                            const std::vector<bool> & funMarked,
                            int refExtension = 0)
{
    GISMO_ASSERT(funMarked.size()==gsMultiBasis<T>(mp).totalSize(),"Vector of marked functions must have the same size as the total number of basis functions (over all patches)!");

    const int dim = mp.dim();

    int numMarked, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes;

    for (size_t pn=0; pn < mp.nPatches(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const int numFun = mp[pn].basis().size();
        numMarked = std::count_if(funMarked.begin() + poffset,
                                  funMarked.begin() + poffset + numFun,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );

        poffset += numFun;
        refBoxes.resize(dim, 2*numMarked);
        numMarked = 0;// counting current patch element to be refined
        for (index_t i = 0; i != mp.basis(pn).size(); ++i)
            if (funMarked[globalCount++])
            {
                refBoxes.block(0,2*numMarked,2,2) = mp.basis(pn).support(i);
                numMarked++;
            }

        // Refine all of the found refBoxes in this patch
        std::vector<index_t> elements = mp.patch(pn).basis().asElements(refBoxes, refExtension);

        mp.patch(pn).refineElements( elements );
    }
}

template <class T>
void gsUnrefineMarkedElements(gsMultiBasis<T> & basis,
                            const std::vector<bool> & elMarked,
                            index_t refExtension = 0)
{
    GISMO_ASSERT(basis.numElements()==elMarked.size(),
                 "Incorrect input, num_el="<<basis.numElements()
                 <<", mark_vec="<<elMarked.size() );

    const short_t dim = basis.dim();

    // numMarked: Number of marked cells on current patch, also currently marked cell
    // poffset  : offset index for the first element on a patch
    // globalCount: counter for the current global element index
    index_t numMarked, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes;

    for (size_t pn=0; pn < basis.nBases(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const size_t numEl = basis[pn].numElements();
        numMarked = std::count_if(elMarked.begin() + poffset,
                                  elMarked.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );
        poffset += numEl;
        refBoxes.resize(dim, 2*numMarked);
        numMarked = 0;// counting current patch element to be refined

        // for all elements in patch pn
        typename gsBasis<T>::domainIter domIt = basis.basis(pn).makeDomainIterator();
        for (; domIt->good(); domIt->next())
        {
            if( elMarked[ globalCount++ ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                refBoxes.col(2*numMarked  ) = domIt->lowerCorner();
                refBoxes.col(2*numMarked+1) = domIt->upperCorner();

                // Advance marked cells counter
                numMarked++;
            }
        }

        // Refine all of the found refBoxes in this patch
        basis.unrefine( pn, refBoxes, refExtension );
    }
}

template <class T>
void gsProcessMarkedElements(gsMultiBasis<T> & basis,
                            const std::vector<bool> & elRefined,
                            const std::vector<bool> & elCoarsened,
                            index_t refExtension = 0,
                            index_t crsExtension = 0)
{
    const short_t dim = basis.dim();

    // numMarked: Number of marked cells on current patch, also currently marked cell
    // poffset  : offset index for the first element on a patch
    // globalCount: counter for the current global element index
    index_t numRefined, numCoarsened, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes, crsBoxes;

    for (size_t pn=0; pn < basis.nBases(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const size_t numEl = basis[pn].numElements();
        numRefined = std::count_if(elRefined.begin() + poffset,
                                  elRefined.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );

        numCoarsened = std::count_if(elCoarsened.begin() + poffset,
                                  elCoarsened.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );


        poffset += numEl;

        refBoxes.resize(dim, 2*numRefined);
        crsBoxes.resize(dim, 2*numCoarsened);
        numRefined = 0;// counting current patch element to be refined
        numCoarsened = 0;// counting current patch element to be refined

        // for all elements in patch pn
        typename gsBasis<T>::domainIter domIt = basis.basis(pn).makeDomainIterator();
        for (; domIt->good(); domIt->next())
        {
            if( elRefined[ globalCount ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                refBoxes.col(2*numRefined  ) = domIt->lowerCorner();
                refBoxes.col(2*numRefined+1) = domIt->upperCorner();

                // Advance marked cells counter
                numRefined++;
            }
            if( elCoarsened[ globalCount ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                crsBoxes.col(2*numCoarsened  ) = domIt->lowerCorner();
                crsBoxes.col(2*numCoarsened+1) = domIt->upperCorner();

                // Advance marked cells counter
                numCoarsened++;
            }

            globalCount++;
        }
        // Refine all of the found refBoxes in this patch
        basis.unrefine( pn, crsBoxes, crsExtension);
        basis.refine( pn, refBoxes, refExtension );
    }
}

template <class T>
void gsProcessMarkedElements(gsMultiPatch<T> & mp,
                            const std::vector<bool> & elRefined,
                            const std::vector<bool> & elCoarsened,
                            index_t refExtension = 0,
                            index_t crsExtension = 0)
{
    const int dim = mp.dim();

    // numMarked: Number of marked cells on current patch, also currently marked cell
    // poffset  : offset index for the first element on a patch
    // globalCount: counter for the current global element index
    index_t numRefined, numCoarsened, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes, crsBoxes;

    for (size_t pn=0; pn < mp.nPatches(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const size_t numEl = mp[pn].basis().numElements();
        numRefined = std::count_if(elRefined.begin() + poffset,
                                  elRefined.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );

        numCoarsened = std::count_if(elCoarsened.begin() + poffset,
                                  elCoarsened.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );


        poffset += numEl;
        refBoxes.resize(dim, 2*numRefined);
        crsBoxes.resize(dim, 2*numCoarsened);
        numRefined = 0;// counting current patch element to be refined
        numCoarsened = 0;// counting current patch element to be refined

        // for all elements in patch pn
        typename gsBasis<T>::domainIter domIt = mp.patch(pn).basis().makeDomainIterator();
        for (; domIt->good(); domIt->next())
        {
            if( elRefined[ globalCount ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                refBoxes.col(2*numRefined  ) = domIt->lowerCorner();
                refBoxes.col(2*numRefined+1) = domIt->upperCorner();

                // Advance marked cells counter
                numRefined++;
            }
            if( elCoarsened[ globalCount ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                crsBoxes.col(2*numCoarsened  ) = domIt->lowerCorner();
                crsBoxes.col(2*numCoarsened+1) = domIt->upperCorner();

                // Advance marked cells counter
                numCoarsened++;
            }

            globalCount++;
        }

        std::vector<index_t> elements;
        // Unrefine all of the found refBoxes in this patch
        elements = mp.patch(pn).basis().asElementsUnrefine(crsBoxes, crsExtension);
        mp.patch(pn).unrefineElements( elements );

        // Refine all of the found refBoxes in this patch
        elements = mp.patch(pn).basis().asElements(refBoxes, refExtension);
        mp.patch(pn).refineElements( elements );
    }
}


template <class T>
void gsUnrefineMarkedElements(gsMultiPatch<T> & mp,
                            const std::vector<bool> & elMarked,
                            int refExtension = 0)
{
    const int dim = mp.dim();

    // numMarked: Number of marked cells on current patch, also currently marked cell
    // poffset  : offset index for the first element on a patch
    // globalCount: counter for the current global element index
    int numMarked, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes;

    for (size_t pn=0; pn < mp.nPatches(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const int numEl = mp[pn].basis().numElements();
        numMarked = std::count_if(elMarked.begin() + poffset,
                                  elMarked.begin() + poffset + numEl,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );

        poffset += numEl;
        refBoxes.resize(dim, 2*numMarked);
        numMarked = 0;// counting current patch element to be refined

        // for all elements in patch pn
        typename gsBasis<T>::domainIter domIt = mp.patch(pn).basis().makeDomainIterator();
        for (; domIt->good(); domIt->next())
        {
            if( elMarked[ globalCount++ ] ) // refine this element ?
            {
                // Construct degenerate box by setting both
                // corners equal to the center
                refBoxes.col(2*numMarked  ) =
                        refBoxes.col(2*numMarked+1) = domIt->centerPoint();

                // Advance marked cells counter
                numMarked++;
            }
        }
        // Refine all of the found refBoxes in this patch
        std::vector<index_t> elements = mp.patch(pn).basis().asElementsUnrefine(refBoxes, refExtension);
        mp.patch(pn).unrefineElements( elements );
    }
}

template <class T>
void gsUnrefineMarkedFunctions(gsMultiPatch<T> & mp,
                            const std::vector<bool> & funMarked,
                            int refExtension = 0)
{
    GISMO_ASSERT(funMarked.size()==gsMultiBasis<T>(mp).totalSize(),"Vector of marked functions must have the same size as the total number of basis functions (over all patches)!");

    const int dim = mp.dim();

    int numMarked, poffset = 0, globalCount = 0;

    // refBoxes: contains marked boxes on a given patch
    gsMatrix<T> refBoxes;

    for (size_t pn=0; pn < mp.nPatches(); ++pn )// for all patches
    {
        // Get number of elements to be refined on this patch
        const int numFun = mp[pn].basis().size();
        numMarked = std::count_if(funMarked.begin() + poffset,
                                  funMarked.begin() + poffset + numFun,
                                  GS_BIND2ND(std::equal_to<bool>(), true) );

        poffset += numFun;
        refBoxes.resize(dim, 2*numMarked);
        numMarked = 0;// counting current patch element to be refined
        for (index_t i = 0; i != mp.basis(pn).size(); ++i)
            if (funMarked[globalCount++])
            {
                refBoxes.block(0,2*numMarked,2,2) = mp.basis(pn).support(i);
                numMarked++;
            }

        // Refine all of the found refBoxes in this patch
        std::vector<index_t> elements = mp.patch(pn).basis().asElements(refBoxes, refExtension);

        mp.patch(pn).unrefineElements( elements );
    }
}


} // namespace gismo