prolongboundarypatch.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_PROLONG_BOUNDARY_PATCH_HH
#define DUNE_FUFEM_PROLONG_BOUNDARY_PATCH_HH
 
#include <vector>
 
#include <dune/grid/common/grid.hh>
#include <dune/grid/common/rangegenerators.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/facehierarchy.hh>
 
#if HAVE_DUNE_UGGRID
#include <dune/grid/uggrid.hh>
#endif
 
template <class GridType>
class PatchProlongator {
public:
 
    static void prolong(std::vector<BoundaryPatch<typename GridType::LevelGridView> >& patches)
    {
        using namespace Dune;
 
        const GridType& grid = patches[0].gridView().grid();
        if (!&grid)
            DUNE_THROW(Exception, "Coarsest boundary patch has not been set up correctly!");
 
        assert(patches.size()>=grid.maxLevel()+1);
 
        // Set array sizes correctly
        for (int i=1; i<=grid.maxLevel(); i++)
            patches[i].setup(grid.levelGridView(i));
 
        // Loop over the boundary on level 0
        typename GridType::LevelGridView levelView = grid.levelGridView(0);
        for (const auto& e : elements(levelView)) {
 
            // Loop over all neighbors
            for (const auto& i : intersections(levelView,e)) {
 
                // if the element is a boundary element
                if (patches[0].contains(i)) {
 
                    // ///////////////////////////////////////////////////////////////
                    //   This implementation doesn't use the UGGrid-specific method
                    //   getChildrenOfSubface.  Instead, it compares normals.  Therefore
                    //   it cannot be used with parametrized boundaries.
                    // //////////////////////////////////////////////////////////////
                    const int dim      = GridType::dimension;
                    const int dimworld = GridType::dimensionworld;
                    typedef typename GridType::ctype ctype;
                    FieldVector<ctype,dim-1> dummy;
                    FieldVector<ctype,dimworld> level0Normal = i.unitOuterNormal(dummy);
 
                    // Find all boundary segments of the descendants of this element
                    for (const auto& h : descendantElements(e,grid.maxLevel())) {
 
                        if (h.level() == e.level())
                            continue;
 
                        typename GridType::LevelGridView hlevelView = grid.levelGridView(h.level());
                        for (const auto& hi : intersections(hlevelView,h)) {
 
                            if (hi.boundary()) {
 
                                // Test whether this boundary segment is contained in the
                                // level 0 boundary segment by comparing the normals.
                                FieldVector<ctype,dim-1> dummy;
                                if (0.95 < level0Normal*hi.unitOuterNormal(dummy))
                                    patches[h.level()].addFace(hi);
 
 
                            }
 
                        }
 
                    }
 
                }
 
            }
 
        }
 
    }
 
    static void prolong(const BoundaryPatch<typename GridType::LevelGridView>& coarsePatch,
                        BoundaryPatch<typename GridType::LeafGridView>& finePatch)
    {
        const GridType& grid = coarsePatch.gridView().grid();
 
        if (!&grid)
            DUNE_THROW(Dune::Exception, "Coarsest boundary patch has not been set up correctly!");
 
        // Set array sizes correctly
        finePatch.setup(grid.leafGridView());
 
        for (const auto& pIt : coarsePatch) {
 
            const auto& inside = pIt.inside();
 
            if (inside.isLeaf())
                finePatch.addFace(inside, pIt.indexInInside());
            else
            {
                Face<GridType> face(grid, inside, pIt.indexInInside());
 
                typename Face<GridType>::HierarchicIterator it = face.hbegin(grid.maxLevel());
                typename Face<GridType>::HierarchicIterator end = face.hend(grid.maxLevel());
 
                for(; it!=end; ++it)
                {
                    if (it->element().isLeaf())
                        finePatch.addFace(it->element(), it->index());
                }
            }
        }
    }
 
};
 
 
#if HAVE_DUNE_UGGRID
template <int dim>
class PatchProlongator<Dune::UGGrid<dim> > {
public:
 
    typedef typename Dune::UGGrid<dim> GridType;
 
    static void prolong(std::vector<BoundaryPatch<typename GridType::LevelGridView> >& patches)
    {
        using namespace Dune;
 
        typedef typename GridType::template Codim<0>::Entity EntityType;
 
        const GridType& grid = patches[0].gridView().grid();
        if (!&grid)
            DUNE_THROW(Exception, "Coarsest boundary patch has not been set up correctly!");
 
        assert(patches.size()>=grid.maxLevel()+1);
 
        // Set array sizes correctly
        for (int i=1; i<=grid.maxLevel(); i++)
            patches[i].setup(grid.levelGridView(i));
 
        // Loop over the boundary on level 0
        typename GridType::LevelGridView levelView = grid.levelGridView(0);
        for (const auto& e : elements(levelView)) {
 
            // Loop over all neighbors
            for (const auto& i : intersections(levelView,e)) {
 
                // if the element is a boundary element
                if (patches[0].contains(i)) {
 
                    // /////////////////////////////////////////////////////////////////
                    //   This version actually only works for UGGrids, but it works
                    //   purely topologically.  For the other grids we need to compare
                    //   surface normals.
                    // /////////////////////////////////////////////////////////////////
                    std::vector<EntityType> childElements(0,*levelView.template begin<0>());
                    std::vector<unsigned char> childElementSides;
                    grid.getChildrenOfSubface(e, i.indexInInside(), grid.maxLevel(),
                                              childElements, childElementSides);
 
                    for (size_t i=0; i<childElements.size(); i++) {
                        const int& level = childElements[i].level();
                        patches[level].addFace(childElements[i], childElementSides[i]);
                    }
 
                }
 
            }
 
        }
 
    }
 
 
 
#if ENABLE_GENERIC_PATCHPROLONGATOR_FOR_UG
    static void prolong(const BoundaryPatchBase<typename GridType::LevelGridView>& coarsePatch,
                        BoundaryPatchBase<typename GridType::LeafGridView>& finePatch)
    {
 
        const GridType& grid = coarsePatch.gridView().grid();
 
        if (!&grid)
            DUNE_THROW(Dune::Exception, "Coarsest boundary patch has not been set up correctly!");
 
        // Set array sizes correctly
        finePatch.setup(grid.leafGridView());
 
        for (const auto& pIt : coarsePatch) {
 
            const auto& inside = pIt.inside();
            if (inside.isLeaf())
                finePatch.addFace(inside, pIt.indexInInside());
            else
            {
                Face<GridType> face(grid, inside, pIt.indexInInside());
 
                typename Face<GridType>::HierarchicIterator it = face.hbegin(grid.maxLevel());
                typename Face<GridType>::HierarchicIterator end = face.hend(grid.maxLevel());
 
                for(; it!=end; ++it)
                {
                    if (it->element().isLeaf())
                        finePatch.addFace(it->element(), it->index());
                }
            }
        }
    }
 
 
 
    static void prolong_old(const BoundaryPatchBase<typename GridType::LevelGridView>& coarsePatch,
                        BoundaryPatchBase<typename GridType::LeafGridView>& finePatch)
#else
    static void prolong(const BoundaryPatch<typename GridType::LevelGridView>& coarsePatch,
                        BoundaryPatch<typename GridType::LeafGridView>& finePatch)
#endif
    {
        typedef typename GridType::template Codim<0>::Entity EntityType;
 
        const GridType& grid = coarsePatch.gridView().grid();
 
        if (!&grid)
            DUNE_THROW(Dune::Exception, "Coarsest boundary patch has not been set up correctly!");
 
        // Set array sizes correctly
        finePatch.setup(grid.leafGridView());
 
        // Loop over the boundary on level 0
        typename GridType::LevelGridView levelView = grid.levelGridView(0);
        for (const auto& e : elements(levelView)) {
 
            // Loop over all neighbors
            for (const auto& i : intersections(levelView,e)) {
 
                // if the element is a boundary element
                if (coarsePatch.contains(i)) {
 
                    // /////////////////////////////////////////////////////////////////
                    //   This version actually only works for UGGrids, but it works
                    //   purely topologically.  For the other grids we need to compare
                    //   surface normals.
                    // /////////////////////////////////////////////////////////////////
                    if (e.isLeaf())
                        finePatch.addFace(e, i.indexInInside());
 
                    std::vector<EntityType> childElements(0, *levelView.template begin<0>());
                    std::vector<unsigned char> childElementSides;
                    grid.getChildrenOfSubface(e, i.indexInInside(), grid.maxLevel(),
                                              childElements, childElementSides);
 
                    for (size_t i=0; i<childElements.size(); i++) {
 
                        if (childElements[i].isLeaf())
                            finePatch.addFace(childElements[i], childElementSides[i]);
 
                    }
 
                }
 
            }
 
        }
 
    }
 
};
#endif
 
#endif