foamgridintersections.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:
#ifndef DUNE_FOAMGRID_INTERSECTIONS_HH
#define DUNE_FOAMGRID_INTERSECTIONS_HH
 
#include <memory>
#include <algorithm>
 
#include <dune/grid/common/intersection.hh>
 
#include <dune/foamgrid/foamgrid/foamgridentity.hh>
#include <dune/foamgrid/foamgrid/foamgridvertex.hh>
#include <dune/foamgrid/foamgrid/foamgridgeometry.hh>
#include <dune/foamgrid/foamgrid/foamgridnulliteratorfactory.hh>
 
namespace Dune {
 
template <class GridImp>
class FoamGridLevelIntersectionIterator;
 
template <class GridImp>
class FoamGridLeafIntersectionIterator;
 
template <class GridImp>
class FoamGridLevelIntersection;
 
template <class GridImp>
class FoamGridLeafIntersection;
 
 
template<class GridImp>
class FoamGridIntersection
{
    enum {dimgrid  = GridImp::dimension};
    enum {dimworld = GridImp::dimensionworld};
 
    // The type used to store coordinates
    typedef typename GridImp::ctype ctype;
    typedef typename GridImp::Traits::template Codim<1>::GeometryImpl GeometryImpl;
 
    // The iterators need access to all members
    friend class FoamGridLevelIntersectionIterator<GridImp>;
    friend class FoamGridLeafIntersectionIterator<GridImp>;
 
    // As weĺl as the implementations
    friend class FoamGridLevelIntersection<GridImp>;
    friend class FoamGridLeafIntersection<GridImp>;
 
    template<typename, typename>
    friend class Dune::Intersection;
 
    FoamGridIntersection()
      : center_(nullptr)
      , facetIndex_(-1)
      , neighbor_(FoamGridNullIteratorFactory<dimgrid, dimworld, ctype>::null())
    {}
 
    FoamGridIntersection(const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* center,
                         int facet)
        : center_(center), facetIndex_(facet), neighbor_(FoamGridNullIteratorFactory<dimgrid, dimworld, ctype>::null())
    {}
 
public:
 
    typedef typename GridImp::template Codim<0>::Entity Entity;
 
        Entity inside() const {
            return Entity(FoamGridEntity<0, dimgrid, GridImp> (center_));
        }
 
 
        Entity outside(/*std::size_t neighborIndex = 0*/) const {
            // Return the 'other' element on the current facet
            return Entity(FoamGridEntity<0, dimgrid, GridImp> ((*neighbor_)));
        }
 
        bool equals(const FoamGridIntersection<GridImp>& i) const {
            return center_==i.center_ && neighbor_ == i.neighbor_ && facetIndex_ == i.facetIndex_;
        }
 
        bool boundary () const {
           return center_->facet_[facetIndex_]->elements_.size()==1;
        }
 
        int boundarySegmentIndex () const
        {
            return center_->facet_[facetIndex_]->boundarySegmentIndex();
        }
 
        GeometryType type () const
        {
            return Dune::GeometryTypes::simplex(dimgrid-1);
        }
 
        int indexInInside () const
        {
            assert(facetIndex_ < center_->corners());
            return facetIndex_;
        }
 
        virtual int indexInOutside(std::size_t neighborIndex = 0) const=0;
 
        FieldVector<ctype, dimworld> outerNormal (const FieldVector<ctype, dimgrid-1>& local) const
        {
            // The intersection normal is a vector that is orthogonal to the element normal
            // and to the intersection itself.
 
            // only works for triangles and lines
            assert(center_->type().isTriangle() || center_->type().isLine());
 
            // Compute vertices
            const auto refElement = ReferenceElements<ctype, dimgrid>::general(center_->type());
 
            // facet vertices (vertex in 1d, edge in 2d)
            int v0 = refElement.subEntity(facetIndex_, 1, 0, dimgrid);
 
            if(dimgrid == 2) {
                //second vertex only exists for an edge
                int v1 = refElement.subEntity(facetIndex_, 1, 1, dimgrid);
                // opposite vertex
                int v2 = (v1+1)%3;
                if (v2==v0)
                    v2 = (v0+1)%3;
                assert(v2!=v0 and v2!=v1);
 
                FieldVector<ctype, dimworld> facet = center_->vertex_[v0]->pos_ - center_->vertex_[v1]->pos_;
                FieldVector<ctype, dimworld> otherEdge = center_->vertex_[v2]->pos_ - center_->vertex_[v1]->pos_;
 
                //Cross product of facet and otherEdge is a scaled element normal
                FieldVector<ctype, dimworld> scaledElementNormal;
 
                if(dimworld == 3) //dimworld==3
                {
                    scaledElementNormal[0] = facet[1]*otherEdge[2] - facet[2]*otherEdge[1];
                    scaledElementNormal[1] = facet[2]*otherEdge[0] - facet[0]*otherEdge[2];
                    scaledElementNormal[2] = facet[0]*otherEdge[1] - facet[1]*otherEdge[0];
                    outerNormal_[0] = facet[1]*scaledElementNormal[2] - facet[2]*scaledElementNormal[1];
                    outerNormal_[1] = facet[2]*scaledElementNormal[0] - facet[0]*scaledElementNormal[2];
                    outerNormal_[2] = facet[0]*scaledElementNormal[1] - facet[1]*scaledElementNormal[0];
                }
                else //dimworld==2
                {
                    outerNormal_[0] = facet[1];
                    outerNormal_[1] = -facet[0];
                }
 
                //Check if scaled EdgeNormal is inner normal, if yes flip
                otherEdge = center_->vertex_[v0]->pos_ - center_->vertex_[v2]->pos_;
                if(otherEdge*outerNormal_ < 0)
                    outerNormal_ *= -1.0;
 
                return outerNormal_;
            }
            else //dimgrid ==1
            {
                int v1 = 1-v0;
                assert(v0!=v1);
                //automatically has right orientation
                outerNormal_ = center_->vertex_[v0]->pos_ - center_->vertex_[v1]->pos_;
 
                return outerNormal_;
            }
            DUNE_THROW(GridError, "Non-existing grid dimension requested! Has to be 1 or 2!");
        }
 
        FieldVector<ctype, dimworld> integrationOuterNormal (const FieldVector<ctype, dimgrid-1>& local) const
        {
 
            if(dimgrid == 2)
            {
                const auto refElement = ReferenceElements<ctype, dimgrid>::general(center_->type());
 
                // facet vertices
                int v0 = refElement.subEntity(facetIndex_, 1, 0, dimgrid);
                int v1 = refElement.subEntity(facetIndex_, 1, 1, dimgrid);
 
                //facet length
                ctype facetLength = (center_->vertex_[v0]->pos_- center_->vertex_[v1]->pos_).two_norm();
 
                FieldVector<ctype, dimworld> integrationOuterNormal_ = unitOuterNormal(local);
                integrationOuterNormal_ *= facetLength;
                return integrationOuterNormal_;
            }
            else //dimgrid == 1
            {
                integrationOuterNormal_ = this->unitOuterNormal(local);
                return integrationOuterNormal_;
            }
            DUNE_THROW(GridError, "Non-existing grid dimension requested! Has to be 1 or 2!");
        }
 
        FieldVector<ctype, dimworld> unitOuterNormal (const FieldVector<ctype, dimgrid-1>& local) const
        {
            unitOuterNormal_ = this->outerNormal(local);
            unitOuterNormal_ /= unitOuterNormal_.two_norm();
            return unitOuterNormal_;
        }
 
        FieldVector<ctype, dimworld> centerUnitOuterNormal () const
        {
            return unitOuterNormal(FieldVector<ctype,dimgrid-1>(0.5));
        }
    private:
 
        mutable FieldVector<ctype, dimworld> outerNormal_;
        mutable FieldVector<ctype, dimworld> unitOuterNormal_;
        mutable FieldVector<ctype, dimworld> integrationOuterNormal_;
 
    protected:
 
        const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* center_;
 
        int facetIndex_;
 
        typename std::vector<const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>*>::const_iterator neighbor_;
 
};
 
template<class GridImp>
class FoamGridLevelIntersection
    : public FoamGridIntersection<GridImp>
{
    friend class FoamGridLevelIntersectionIterator<GridImp>;
 
    template<typename, typename>
    friend class Dune::Intersection;
 
    FoamGridLevelIntersection() : FoamGridIntersection<GridImp>() {}
 
    enum {dimgrid = GridImp::dimension};
    enum{ dimworld = GridImp::dimensionworld };
 
    // Geometry is a CachedMultiLinearGeometry
    typedef typename GridImp::template Codim<1>::Geometry Geometry;
    typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
 
    typedef typename GridImp::Traits::template Codim<1>::GeometryImpl GeometryImpl;
    typedef typename GridImp::Traits::template Codim<1>::LocalGeometryImpl LocalGeometryImpl;
 
    typedef typename GridImp::ctype ctype;
 
    FoamGridLevelIntersection(const FoamGridEntityImp<dimgrid, dimgrid ,dimworld, ctype>* center, std::size_t facet)
                              : FoamGridIntersection<GridImp>(center, facet)
    {}
 
public:
 
    bool conforming () const {
        // FoamGrid level intersections are always conforming
        return true;
    }
 
    int indexInOutside (std::size_t neighborIndex = 0) const {
        return std::find((*this->neighbor_)->facet_.begin(), (*this->neighbor_)->facet_.end(),
                         this->center_->facet_[this->facetIndex_])
            - (*this->neighbor_)->facet_.begin();
    }
 
    bool neighbor () const {
      return this->neighbor_!=neighborEnd_;
    }
 
    LocalGeometry geometryInInside () const
    {
        std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
 
        // Get reference element
        const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
 
        for (int idx = 0; idx < dimgrid; ++idx)
            coordinates[idx] = refElement.position(refElement.subEntity(this->facetIndex_, 1, idx, dimgrid), dimgrid);
 
        geometryInInside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
 
      return LocalGeometry(*geometryInInside_);
    }
 
    LocalGeometry geometryInOutside (std::size_t neighborIndex = 0) const {
 
        // Get two vertices of the intersection
        const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
 
        std::array<FoamGridEntityImp<0, dimgrid, dimworld, ctype>*, dimgrid> vtx;
 
        for (int idx = 0; idx < dimgrid; ++idx)
            vtx[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)];
 
        std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
 
        // Find the intersection vertices in local numbering of the outside element
        // That way we get the local orientation correctly.
        const auto refElementOther = ReferenceElements<ctype, dimgrid>::general((*this->neighbor_)->type());
 
        for (std::size_t j=0; j<dimgrid; j++)
          for (int i=0; i<refElementOther.size(dimgrid); i++)
             if (vtx[j] == (*this->neighbor_)->vertex_[refElementOther.subEntity(0, 0, i, dimgrid)])
              coordinates[j] = refElement.position(refElement.subEntity(0, 0, i, dimgrid), dimgrid);
 
        geometryInOutside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
 
        return LocalGeometry(*geometryInOutside_);
    }
 
    Geometry geometry () const {
 
        std::vector<FieldVector<ctype, dimworld> > coordinates(dimgrid);
 
        // Get two vertices of the intersection
        const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
 
        for (std::size_t idx = 0; idx < dimgrid; ++idx)
            coordinates[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)]->pos_;
 
        geometry_ = std::make_shared<GeometryImpl>(this->type(), coordinates);
 
        return Geometry(*geometry_);
    }
 
 
    private:
    typename std::vector<const FoamGridEntityImp<dimgrid, dimgrid ,dimworld, ctype>*>::const_iterator neighborEnd_;
    mutable std::shared_ptr<GeometryImpl> geometry_;
    mutable std::shared_ptr<LocalGeometryImpl> geometryInInside_;
    mutable std::shared_ptr<LocalGeometryImpl> geometryInOutside_;
 
};
 
 
 
 
template<class GridImp>
class FoamGridLeafIntersection
    : public FoamGridIntersection<GridImp>
{
 
    friend class FoamGridLeafIntersectionIterator<GridImp>;
 
    template<typename, typename>
    friend class Dune::Intersection;
 
    FoamGridLeafIntersection() : FoamGridIntersection<GridImp>() {}
 
    enum {dimworld = GridImp::dimensionworld};
    enum {dimgrid  = GridImp::dimension};
 
    typedef typename GridImp::ctype ctype;
 
    typedef typename GridImp::template Codim<1>::Geometry Geometry;
    typedef typename GridImp::template Codim<1>::LocalGeometry LocalGeometry;
 
    typedef typename GridImp::Traits::template Codim<1>::GeometryImpl GeometryImpl;
    typedef typename GridImp::Traits::template Codim<1>::LocalGeometryImpl LocalGeometryImpl;
 
    FoamGridLeafIntersection(const FoamGridEntityImp<dimgrid, dimgrid, dimworld, ctype>* center,
                             std::size_t facet)
        : FoamGridIntersection<GridImp>(center, facet)
    {}
 
public:
 
    bool conforming () const {
        // FoamGrid leaf? intersections are always? conforming
        return true;
    }
 
    int indexInOutside (std::size_t neighborIndex = 0) const
    {
        // get our facet object
        auto* facet = this->center_->facet_[this->facetIndex_];
 
        // if the neighbor level is greater than the facet level
        // walk until our facet's father that has the same level as the neighbor
        // and look for the father facet in the neighbor to get the index
        while(facet->level() > (*this->neighbor_)->level())
        {
          assert(facet->father_ != nullptr);
          facet = facet->father_;
        }
 
        if (facet->level() == (*this->neighbor_)->level())
        {
          return std::distance((*this->neighbor_)->facet_.begin(),
                                std::find((*this->neighbor_)->facet_.begin(),
                                          (*this->neighbor_)->facet_.end(),
                                          facet)
                              );
 
        }
 
        // the neighbor level is smaller than the facet level
        // find the neighbor's facet that has the center element in it's element list
        else
        {
          return std::distance((*this->neighbor_)->facet_.begin(),
                                std::find_if((*this->neighbor_)->facet_.begin(),
                                             (*this->neighbor_)->facet_.end(),
                                             [this](const auto& facet) -> bool {
                                               for (const auto& e : facet->elements_)
                                                 if (this->center_ == e)
                                                   return true;
                                               return false;
                                             })
                              );
        }
    }
 
    LocalGeometry geometryInInside () const
    {
        std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
 
        // Get reference element
        const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
 
        for (std::size_t idx = 0; idx < dimgrid; ++idx)
            coordinates[idx] = refElement.position(refElement.subEntity(this->facetIndex_, 1, idx, dimgrid), dimgrid);
 
        geometryInInside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
 
      return LocalGeometry(*geometryInInside_);
    }
 
    LocalGeometry geometryInOutside (std::size_t neighborIndex = 0) const
    {
        // Get reference element
        const auto refElement = Dune::ReferenceElements<ctype, dimgrid>::general(this->center_->type());
 
        // Map the global position of the intersection vertices to the local position in the neighbor
        std::vector<FieldVector<ctype, dimgrid> > coordinates(dimgrid);
        for (std::size_t vIdx = 0; vIdx < dimgrid; ++vIdx)
        {
            const auto vIdxInElement = refElement.subEntity(this->facetIndex_, 1, vIdx, dimgrid);
            coordinates[vIdx] = (*this->neighbor_)->globalToLocal(this->center_->vertex_[vIdxInElement]->pos_);
        }
 
        geometryInOutside_ = std::make_shared<LocalGeometryImpl>(this->type(), coordinates);
 
        return LocalGeometry(*geometryInOutside_);
    }
 
    Geometry geometry () const {
 
        std::vector<FieldVector<ctype, dimworld> > coordinates(dimgrid);
 
        // Get two vertices of the intersection
        const auto refElement = ReferenceElements<ctype, dimgrid>::general(this->center_->type());
 
        for (std::size_t idx = 0; idx < dimgrid; ++idx)
            coordinates[idx] = this->center_->vertex_[refElement.subEntity(this->facetIndex_, 1, idx, dimgrid)]->pos_;
 
        geometry_ = std::make_shared<GeometryImpl>(this->type(), coordinates);
 
        return Geometry(*geometry_);
    }
 
    bool neighbor () const {
        return this->neighbor_!=FoamGridNullIteratorFactory<dimgrid, dimworld, ctype>::null();
    }
 
private:
 
    mutable std::shared_ptr<GeometryImpl> geometry_;
    mutable std::shared_ptr<LocalGeometryImpl> geometryInInside_;
    mutable std::shared_ptr<LocalGeometryImpl> geometryInOutside_;
};
 
 
}  // namespace Dune
 
#endif