boundarypatch.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 BOUNDARY_PATCH_HH
#define BOUNDARY_PATCH_HH
 
#include <vector>
#include <set>
 
#include <dune/common/bitsetvector.hh>
 
#include <dune/geometry/referenceelements.hh>
#include <dune/grid/common/mcmgmapper.hh>
#include <dune/grid/common/gridview.hh>
#include <dune/grid/common/gridfactory.hh>
 
#include <dune/fufem/domains/boundary.hh>
#include <dune/fufem/domains/boundarypatch.hh>
 
 
 
template <class GridView>
class BoundaryPatchInsertionIndexProperty
{
    typedef typename Dune::GridFactory<typename GridView::Grid> GridFactory;
 
public:
    typedef typename GridView::Intersection Intersection;
 
    template<class T>
    BoundaryPatchInsertionIndexProperty(const GridFactory& factory, const T& t) : factory_(factory)
    {
        indices_.insert(t.begin(), t.end());
    }
 
    BoundaryPatchInsertionIndexProperty(const GridFactory& factory, int i) : factory_(factory)
    {
        indices_.insert(i);
    }
 
    bool operator() (const Intersection& i) const
    {
        return (factory_.wasInserted(i) and (indices_.find(factory_.insertionIndex(i))!=indices_.end()));
    }
private:
    const GridFactory& factory_;
    std::set<int> indices_;
};
 
 
 
template <class GridView>
class BoundaryPatchSegmentIndexProperty
{
public:
    typedef typename GridView::Intersection Intersection;
 
    template<class T>
    BoundaryPatchSegmentIndexProperty(const T& t)
    {
        indices_.insert(t.begin(), t.end());
    }
 
    BoundaryPatchSegmentIndexProperty(int i)
    {
        indices_.insert(i);
    }
 
    bool operator() (const Intersection& i) const
    {
        return (indices_.find(i.boundarySegmentIndex())!=indices_.end());
    }
private:
    std::set<int> indices_;
};
 
 
 
template <class GridView, int ncomp = 1>
class BoundaryPatchEnclosingVerticesProperty
{
    typedef typename GridView::IndexSet IndexSet;
    static const int dim = GridView::dimension;
public:
    typedef typename GridView::Intersection Intersection;
 
    BoundaryPatchEnclosingVerticesProperty(const GridView& gridView, const Dune::BitSetVector<ncomp>& vertices) :
        indexSet_(gridView.indexSet()),
        vertices_(vertices)
    {}
 
    bool operator() (const Intersection& i) const
    {
        const auto inside = i.inside();
        int localFaceIndex = i.indexInInside();
 
        auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
 
        // Get global node ids
        int n = refElement.size(localFaceIndex, 1, dim);
 
        // Using ReferenceElement::subEntity is OK here, because we loop
        // over all subEntities (i.e. sub-vertices) and just return false
        // if _any_ of them is not marked.
        for (int i=0; i<n; i++) {
            int localVertexIndex = refElement.subEntity(localFaceIndex, 1, i, dim);
            if (not(vertices_[indexSet_.subIndex(inside, localVertexIndex, dim)].any()))
                return false;
        }
        return true;
    }
 
private:
    const IndexSet& indexSet_;
    const Dune::BitSetVector<ncomp>& vertices_;
};
 
template <class GridView>
class BoundaryPatchTouchingVerticesProperty
{
    typedef typename GridView::IndexSet IndexSet;
    static const int dim = GridView::dimension;
public:
    typedef typename GridView::Intersection Intersection;
 
    BoundaryPatchTouchingVerticesProperty(const GridView& gridView, const Dune::BitSetVector<1>& vertices) :
        indexSet_(gridView.indexSet()),
        vertices_(vertices)
    {}
 
    bool operator() (const Intersection& i) const
    {
        const auto inside = i.inside();
        int localFaceIndex = i.indexInInside();
 
        auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
 
        // Get global node ids
        int n = refElement.size(localFaceIndex, 1, dim);
 
        // Using ReferenceElement::subEntity is OK here, because we loop
        // over all subEntities (i.e. sub-vertices) and return true as soon
        // as _any_ of them is marked.
        for (int i=0; i<n; i++) {
            int localVertexIndex = refElement.subEntity(localFaceIndex, 1, i, dim);
            if (vertices_[indexSet_.subIndex(inside, localVertexIndex, dim)][0])
                return true;
        }
        return false;
    }
 
private:
    const IndexSet& indexSet_;
    const Dune::BitSetVector<1>& vertices_;
};
 
 
 
 
template <class GV>
class BoundaryPatch
    : public Dune::Fufem::BoundaryPatch<GV>
{
protected:
 
    using Base = typename Dune::Fufem::BoundaryPatch<GV>;
    using Element = typename Base::Element;
    using Intersection = typename Base::Intersection;
    using IntersectionIterator = typename GV::IntersectionIterator;
 
public:
 
    using iterator = typename Base::iterator;
 
    using GridView = GV;
 
    using Base::dim;
    using Base::gridView;
    using Base::begin;
    using Base::end;
    using Base::contains;
    using Base::isInitialized;
    using Base::insertFace;
    using Base::size;
 
 
    BoundaryPatch(const GridView& gridView, bool initialValue = false)
        : Base()
        , verticesAreUpToDate_(false)
    {
        setup(gridView, initialValue);
    }
 
    template <int blocksize>
    BoundaryPatch(const GridView& gridView, const Dune::BitSetVector<blocksize>& vertices)
        : Base()
        , verticesAreUpToDate_(false)
    {
        setup(gridView, vertices);
    }
 
    BoundaryPatch() = default;
 
    BoundaryPatch(const BoundaryPatch&) = default;
 
    BoundaryPatch(BoundaryPatch&&) = default;
 
    BoundaryPatch& operator=(const BoundaryPatch&) = default;
 
    BoundaryPatch& operator=(BoundaryPatch&&) = default;
 
    void setup(const GridView& gridView, bool insertAllBoundaryFaces = false)
    {
        Base::setGridView(gridView);
        if (insertAllBoundaryFaces)
            for (const auto& intersection : Dune::Fufem::Boundary(this->gridView()))
                insertFace(intersection);
        verticesAreUpToDate_ = false;
    }
 
    template <int blocksize>
    void setup(const GridView& gridView, const Dune::BitSetVector<blocksize>& vertices)
    {
        // setup empty BoundaryPatch
        setup(gridView);
 
        // create property for insertion by vertex vector
        BoundaryPatchEnclosingVerticesProperty<GridView, blocksize> prop(this->gridView(), vertices);
 
        // insert by property
        insertFacesByProperty(prop);
    }
 
    void clear() {
        Base::clear();
        vertices_.unsetAll();
        verticesAreUpToDate_ = false;
    }
 
    void getFaces(Dune::BitSetVector<1>& faces) const {
        faces = Base::faces_;
    }
 
    int numFaces() const {
        return Base::size();
    }
 
    void addFace(const IntersectionIterator& nIt) {
        insertFace(*nIt);
    }
 
    void addFace(const Intersection& nIt) {
        insertFace(nIt);
    }
 
    void addFace(const Element& en, int fIdx) {
        insertFace(en, fIdx);
    }
 
    template<class InsertionProperty>
    void insertFacesByProperty(const InsertionProperty& property)
    {
        for(const auto& element : elements(gridView()))
            if (element.hasBoundaryIntersections())
                for (const auto& is : intersections(gridView(), element))
                    if (is.boundary())
                        if (property(is))
                            insertFace(is);
        this->verticesAreUpToDate_ = false;
    }
 
    bool contains(const IntersectionIterator& nIt) const
    {
        return contains(*nIt);
    }
 
 
    bool containsVertex(size_t v) const {
        if (!verticesAreUpToDate_)
            makeVertices();
 
#ifdef DUNE_RANGE_CHECKING
        if (v<0 || v>=vertices_.size())
            DUNE_THROW(Dune::RangeError, "Trying to check for vertex " << v
                       << ", but the vertices are only in the range [0..." << int(vertices_.size())-1
                       << "]!");
#endif
 
        return vertices_[v][0];
    }
 
    typename GridView::Grid::ctype area() const {
        auto patchArea = typename GridView::Grid::ctype(0);
        for(const auto& intersection : *this)
            patchArea += intersection.geometry().volume();
        return patchArea;
    }
 
    void getVertices(Dune::BitSetVector<1>& vertices) const {
        if (!this->verticesAreUpToDate_)
            this->makeVertices();
 
        vertices = this->vertices_;
    }
 
    const Dune::BitSetVector<1>* getVertices() const {
        if (!this->verticesAreUpToDate_)
            this->makeVertices();
 
        return &this->vertices_;
    }
 
    void getPatchBoundaryVertices(Dune::BitSetVector<1>& patchBoundary) const {
 
        const auto& indexSet = gridView().indexSet();
 
        patchBoundary.resize(indexSet.size(dim));
        patchBoundary.unsetAll();
 
        // loop over all elements
        for (const auto& e : elements(gridView())) {
 
            auto refElement = Dune::ReferenceElements<double, dim>::general(e.type());
 
            // Loop over all neighbors
            for (const auto& is : intersections(gridView(),e)) {
 
                // if the element is a boundary element, but _not_ a patch element
                if (is.boundary() && !contains(is)) {
 
                    // Get global node ids
                    int n = refElement.size(is.indexInInside(),1,dim);
 
                    // ATTENTION: This implementation is broken!
                    //
                    // It will not work in general, because there's no guarantee that:
                    //
                    //   subIndex(e, refElement.subEntity(face,1,i,dim)) == subIndex(e.subEntity<1>(face), i, dim)
                    //
                    // The latter is what you really want because it uses the correct embedding
                    // of the face into the element (which is not the same as the one in the
                    // reference element)
                    for (int i=0; i<n; i++) {
                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
                        int globalIdx = indexSet.subIndex(e, faceIdxi,dim);
                        if (this->containsVertex(globalIdx))
                            patchBoundary[globalIdx] = true;
                    }
 
                }
 
            }
 
        }
 
    }
 
    void getPatchBoundaryEdges(Dune::BitSetVector<1>& patchBoundaryEdges) const {
 
        assert(dim==3);
        const auto& indexSet = gridView().indexSet();
 
        patchBoundaryEdges.resize(indexSet.size(dim-1));
        patchBoundaryEdges.unsetAll();
 
        // loop over all elements
        for (const auto& e : elements(gridView())) {
 
            auto refElement = Dune::ReferenceElements<double, dim>::general(e.type());
 
            // Loop over all neighbors
            for (const auto& is : intersections(gridView(),e)) {
 
                // if the element is a boundary element, but _not_ a patch element
                if (is.boundary() && !contains(is)) {
 
                    // Loop over the face edges
                    int n = refElement.size(is.indexInInside(),1,dim-1);
 
                    for (int i=0; i<n; i++) {
 
                        // ATTENTION: This implementation is broken!
                        //
                        // It will not work in general, because there's no guarantee that:
                        //
                        //   subIndex(e, refElement.subEntity(face,1,i,dim-1)) == subIndex(e.subEntity<1>(face), i, dim-1)
                        //
                        // The latter is what you really want because it uses the correct embedding
                        // of the face into the element (which is not the same as the one in the
                        // reference element).
                        int edge = refElement.subEntity(is.indexInInside(), 1, i, dim-1);
 
                        // If 'edge' would be correct, the following use of ReferenceElement::subEntity
                        // would be OK, because we just do a check for all vertices we pass, i.e. 0 and 1.
                        int v0 = refElement.subEntity(edge, dim-1, 0, dim);
                        int v1 = refElement.subEntity(edge, dim-1, 1, dim);
 
                        if (containsVertex(indexSet.subIndex(e, v0,dim))
                            && containsVertex(indexSet.subIndex(e, v1,dim)))
                            patchBoundaryEdges[indexSet.subIndex(e, edge,dim-1)] = true;
                    }
 
                }
 
            }
 
        }
 
    }
 
    int numVertices() const {
        if (not(isInitialized()))
            return 0;
        if (!this->verticesAreUpToDate_)
            this->makeVertices();
 
        return this->vertices_.count();
    }
 
    void getNormals(std::vector<Dune::FieldVector<double, dim> >& normals) const {
        normals = getNormals();
    }
 
    std::vector<Dune::FieldVector<double,dim> > getNormals() const {
 
        const auto& indexSet = gridView().indexSet();
 
        if (!this->verticesAreUpToDate_)
            this->makeVertices();
 
        typedef Dune::FieldVector<double, dim> Vector;
        std::vector<Vector> normals(this->vertices_.size(), Vector(0.0));
 
        // Loop over all boundary segments
        iterator it    = begin();
        iterator endIt = end();
 
        for (; it!=endIt; ++it) {
            const auto inside = it->inside();
 
            auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
 
            // Compute actual surface normals
            int n = refElement.size(it->indexInInside(),1,dim);
 
            for (int i=0; i<n; i++) {
                // ATTENTION: This implementation does only work if unitOuterNormal is constant on the intersection.
                // It will not work in general, because there's no guarantee that:
                //
                //   subIndex(e, refElement.subEntity(inside,1,i,dim)) == subIndex(e.subEntity<1>(inside), i, dim)
                //
                // The latter is what you really want because it uses the correct embedding
                // of the face into the element (which is not the same as the one in the
                // reference element).
                //
                // If unitOuterNormal is constant we add the same faceNormal to all vertices
                // of this face, hence the proper order does not matter. Otherwise this
                // implementation may compute the faceNormal at the wrong posInBoundary
                // because posInElement and vertexIdx don't correspond to each other.
                int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
 
                // local position of vertex in element
                const Dune::FieldVector<double, dim>& posInElement = refElement.position(faceIdxi, dim);
 
                // local position in boundary segment
                Dune::FieldVector<double, dim-1> posInBoundary = it->geometryInInside().local(posInElement);
 
                // Get normal there
                const Dune::FieldVector<double, dim>& faceNormal = it->unitOuterNormal(posInBoundary);
 
                int vertexIdx = indexSet.subIndex(inside, faceIdxi,dim);
                normals[vertexIdx] += faceNormal;
            }
 
        }
 
        // Normalize vertex normals
        for (size_t i=0; i<normals.size(); i++)
            if (this->vertices_[i][0])
                normals[i] /= normals[i].two_norm();
 
        return normals;
    }
 
    std::vector<Dune::FieldVector<double,dim> > getLocalNormals() const {
 
        const auto& indexSet = gridView().indexSet();
 
        if (!this->verticesAreUpToDate_)
            this->makeVertices();
 
        auto globalToLocal = makeGlobalToLocal();
 
        using Vector = Dune::FieldVector<double, dim>;
        std::vector<Vector> normals(this->vertices_.count(), Vector(0.0));
 
        // Loop over all boundary segments
        for (const auto& it : *this) {
 
            const auto inside = it.inside();
 
            auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
 
            // Compute actual surface normals
            int n = refElement.size(it.indexInInside(),1,dim);
 
            const auto& geomInInside = it.geometryInInside();
            for (int i=0; i<n; i++) {
 
                // See the Warning in the method above!
                int faceIdxi = refElement.subEntity(it.indexInInside(), 1, i, dim);
 
                // local position of vertex in element
                auto& posInElement = refElement.position(faceIdxi, dim);
 
                // local position in boundary segment
                auto posInBoundary = geomInInside.local(posInElement);
 
                // Get normal there
                const auto& faceNormal = it.unitOuterNormal(posInBoundary);
 
                int vertexIdx = indexSet.subIndex(inside, faceIdxi,dim);
                normals[globalToLocal[vertexIdx]] += faceNormal;
            }
 
        }
 
        // Normalize vertex normals
        for (size_t i=0; i<normals.size(); i++)
            normals[i] /= normals[i].two_norm();
 
        return normals;
    }
 
    void makeGlobalToLocal(std::vector<int>& globalToLocal) const {
        globalToLocal = makeGlobalToLocal();
    }
 
    std::vector<int> makeGlobalToLocal() const {
        if (!this->verticesAreUpToDate_)
            this->makeVertices();
 
        std::vector<int> globalToLocal(this->vertices_.size());
        int idx = 0;
        for (size_t i=0; i< globalToLocal.size(); i++)
            globalToLocal[i] = (this->vertices_[i][0]) ? idx++ : -1;
 
        return globalToLocal;
    }
 
 
    BoundaryPatch<GridView> addPatch(const BoundaryPatch<GridView>& patch) {
 
        /* test whether both patches live on the same grid
         * for LevelBoundaryPatches it is not tested whether both patches refer to the same level*/
        if (&this->gridView().grid() != &patch.gridView().grid())
            DUNE_THROW(Dune::NotImplemented, "addPatch only works for patches defined on the same grid");
 
        Dune::BitSetVector<1> vertices, additionalVertices;
        this->getVertices(vertices);
        patch.getVertices(additionalVertices);
 
        // create combined BitSetVector
        for (size_t i=0; i<vertices.size(); i++)
            if (additionalVertices[i][0])
                vertices[i] = true;
 
        GridView gridView = this->gridView();
        this->setup(gridView, vertices);
 
        return *this;
    }
 
 
    BoundaryPatch<GridView> intersect(const BoundaryPatch<GridView>& patch) {
 
        /* test whether both patches live on the same grid
         * for LevelBoundaryPatches it is not tested whether both patches refer to the same level*/
        if (&this->gridView().grid() != &patch.gridView().grid())
            DUNE_THROW(Dune::NotImplemented, "addPatch only works for patches defined on the same grid");
 
        Dune::BitSetVector<1> vertices, intersectionVertices;
        this->getVertices(vertices);
        patch.getVertices(intersectionVertices);
 
        // create intersection BitSetVector
        for (size_t i=0; i<vertices.size(); i++)
            if (vertices[i][0] && intersectionVertices[i][0])
                vertices[i][0] = true;
            else
                vertices[i][0] = false;
 
        GridView gridView = this->gridView();
        this->setup(gridView, vertices);
 
        return *this;
    }
 
protected:
 
    void makeVertices() const {
 
        const auto& indexSet = gridView().indexSet();
 
        this->vertices_.resize(gridView().size(dim));
        this->vertices_.unsetAll();
 
        // loop over all elements
        iterator it    = begin();
        iterator endIt = end();
 
        for (; it!=endIt; ++it) {
            const auto inside = it->inside();
 
            auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
 
            // Get global node ids
            int n = refElement.size(it->indexInInside(),1,dim);
 
            // This use of ReferenceElement::subEntity is OK, because
            // we mark the vertices corresponding to all returned indices.
            for (int i=0; i<n; i++) {
                int faceIdxi = refElement.subEntity(it->indexInInside(), 1, i, dim);
                vertices_[indexSet.subIndex(inside,faceIdxi,dim)] = true;
            }
 
        }
 
        verticesAreUpToDate_ = true;
 
    }
 
 
 
    mutable Dune::BitSetVector<1> vertices_;
 
    mutable bool verticesAreUpToDate_ = false;
};
 
 
 
 
#endif