foamgridfactory.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_FACTORY_HH
#define DUNE_FOAMGRID_FACTORY_HH
 
#include <vector>
#include <memory>
#include <unordered_map>
#include <functional>
 
#include <dune/common/fvector.hh>
#include <dune/common/hash.hh>
 
#include <dune/grid/common/gridfactory.hh>
#include <dune/foamgrid/foamgrid.hh>
 
namespace Dune {
 
template <int dimgrid, int dimworld, class ct>
    class GridFactoryBase
        : public GridFactoryInterface<FoamGrid<dimgrid, dimworld, ct> >
    {
    using ctype = ct;
 
    public:
 
        using GridPtrType = std::unique_ptr<FoamGrid<dimgrid, dimworld,ctype>>;
 
        GridFactoryBase()
        : factoryOwnsGrid_(true)
        {
            grid_ = new FoamGrid<dimgrid, dimworld, ctype>;
            grid_->entityImps_.resize(1);
        }
 
        GridFactoryBase(FoamGrid<dimgrid, dimworld, ctype>* grid)
        : factoryOwnsGrid_(false)
        {
            grid_ = grid;
            grid_->entityImps_.resize(1);
        }
 
        ~GridFactoryBase() override {
            if (grid_ && factoryOwnsGrid_)
                delete grid_;
        }
 
        void insertVertex(const FieldVector<ctype,dimworld>& pos) override {
            std::get<0>(grid_->entityImps_[0]).emplace_back(
              0,    // level
              pos,  // position
              grid_->getNextFreeId()
            );
            vertexArray_.push_back(&*std::get<0>(grid_->entityImps_[0]).rbegin());
        }
 
        unsigned int
        insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::Traits::template Codim<0>::Entity &entity) const override
        {
            return entity.impl().target_->leafIndex_;
        }
 
        unsigned int
        insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::Traits::template Codim<dimgrid>::Entity &vertex) const override
        {
            return vertex.impl().target_->leafIndex_;
        }
 
        unsigned int
        insertionIndex(const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection& intersection) const override
        {
            return intersection.boundarySegmentIndex();
        }
 
    protected:
 
        // Pointer to the grid being built
        FoamGrid<dimgrid, dimworld, ctype>* grid_;
 
        // True if the factory allocated the grid itself, false if the
        // grid was handed over from the outside
        bool factoryOwnsGrid_;
 
        std::vector<FoamGridEntityImp<0, dimgrid, dimworld, ctype>*> vertexArray_;
 
        unsigned int boundarySegmentCounter_ = 0;
    };
 
template <int dimgrid, int dimworld, class ct>
    class GridFactory<FoamGrid<dimgrid, dimworld, ct> >
        : public GridFactoryBase<dimgrid, dimworld, ct>
    {};
 
template <int dimworld, class ct>
    class GridFactory<FoamGrid<1, dimworld, ct> >
        : public GridFactoryBase<1, dimworld, ct>
    {
        enum {dimgrid = 1};
 
        using ctype = ct;
 
        template<int mydim>
        using EntityImp = FoamGridEntityImp<mydim, dimgrid, dimworld, ctype>;
 
        using GridPtrType = typename GridFactoryBase<1, dimworld, ct>::GridPtrType;
 
    public:
 
        GridFactory() {}
 
        GridFactory(FoamGrid<1, dimworld, ctype>* grid):
            GridFactoryBase<1,dimworld,ctype>(grid)
        {}
 
        void insertBoundarySegment(const std::vector<unsigned int>& vertices) override
        {
            boundarySegmentIndices_[ vertices[0] ] = this->boundarySegmentCounter_++;
        }
 
        void insertBoundarySegment(const std::vector<unsigned int>& vertices,
                                   const std::shared_ptr<BoundarySegment<dimgrid, dimworld> >& boundarySegment) override
        {
            DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented");
        }
 
        bool wasInserted( const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection &intersection ) const override
        {
          if ( !intersection.boundary() || intersection.inside().level() != 0 )
            return false;
 
          const auto& vertex = intersection.inside().template subEntity<1>(intersection.indexInInside());
          const auto& it = boundarySegmentIndices_.find( this->insertionIndex(vertex) );
          return (it != boundarySegmentIndices_.end());
        }
 
        void insertElement(const GeometryType& type,
                           const std::vector<unsigned int>& vertices) override
        {
            assert(type.isLine());
            // insert new element
            std::get<1>(this->grid_->entityImps_[0]).emplace_back(
                this->vertexArray_[vertices[0]],
                this->vertexArray_[vertices[1]],
                0, // level
                this->grid_->getNextFreeId()
            );
        }
 
        void insertElement(const GeometryType& type,
                           const std::vector<unsigned int>& vertices,
                           std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dimgrid>)> elementParametrization)
        override
        {
            assert(type.isLine());
            EntityImp<1> newElement(this->vertexArray_[vertices[0]],
                                    this->vertexArray_[vertices[1]],
                                    0, // level
                                    this->grid_->getNextFreeId());
 
            newElement.elementParametrization_ = elementParametrization;
            std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
        }
 
        GridPtrType createGrid() override
        {
            // Prevent a crash when this method is called twice in a row
            // You never know who may do this...
            if (this->grid_==nullptr)
                return nullptr;
 
            // make facets (vertices) know about the element
            for (auto& element : std::get<dimgrid>(this->grid_->entityImps_[0]))
            {
                element.vertex_[0]->elements_.push_back(&element);
                element.vertex_[1]->elements_.push_back(&element);
            }
 
            // Create the index sets
            this->grid_->setIndices();
 
            // ////////////////////////////////////////////////
            //   Set the boundary ids
            // ////////////////////////////////////////////////
 
            for (auto& facet : std::get<0>(this->grid_->entityImps_[0]))
            {
                if (facet.elements_.size() == 1) // if boundary facet
                {
                    const auto it = boundarySegmentIndices_.find( facet.vertex_[0]->leafIndex_ );
                    if (it != boundarySegmentIndices_.end())
                        facet.boundarySegmentIndex_ = it->second;
                    else
                        facet.boundarySegmentIndex_ = this->boundarySegmentCounter_++;
                }
            }
 
            // ////////////////////////////////////////////////
            //   Hand over the new grid
            // ////////////////////////////////////////////////
 
            Dune::FoamGrid<dimgrid, dimworld, ctype>* tmp = this->grid_;
            tmp->numBoundarySegments_ = this->boundarySegmentCounter_;
            this->grid_ = nullptr;
            return GridPtrType(tmp);
        }
 
      private:
        std::unordered_map<unsigned int, unsigned int> boundarySegmentIndices_;
    };
 
template <int dimworld, class ct>
    class GridFactory<FoamGrid<2, dimworld, ct> >
        : public GridFactoryBase<2, dimworld, ct>
    {
        enum {dimgrid = 2};
 
        using ctype = ct;
 
        template<int mydim>
        using EntityImp = FoamGridEntityImp<mydim, dimgrid, dimworld, ctype>;
 
        using GridPtrType = typename GridFactoryBase<2, dimworld, ct>::GridPtrType;
 
    public:
 
        GridFactory() {}
 
        GridFactory(FoamGrid<2, dimworld, ctype>* grid):
            GridFactoryBase<2, dimworld, ctype>(grid)
        {}
 
        void insertBoundarySegment(const std::vector<unsigned int>& vertices) override
        {
            std::array<unsigned int, 2> vertexIndices {{ vertices[0], vertices[1] }};
 
            // sort the indices
            if ( vertexIndices[0] > vertexIndices[1] )
              std::swap( vertexIndices[0], vertexIndices[1] );
 
            boundarySegmentIndices_[ vertexIndices ] = this->boundarySegmentCounter_++;
        }
 
        void insertBoundarySegment(const std::vector<unsigned int>& vertices,
                                   const std::shared_ptr<BoundarySegment<dimgrid, dimworld> >& boundarySegment) override
        {
            DUNE_THROW(Dune::NotImplemented, "Parameterized boundary segments are not implemented");
        }
 
        bool wasInserted( const typename FoamGrid<dimgrid, dimworld, ctype>::LeafIntersection &intersection ) const override
        {
          if ( !intersection.boundary() || intersection.inside().level() != 0 )
            return false;
 
          // Get the vertices of the intersection
          const auto refElement = ReferenceElements<ctype, dimgrid>::general(intersection.inside().type());
 
          const int subIdx0 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/0, dimgrid);
          const auto vertex0 = intersection.inside().template subEntity<2>( subIdx0 );
          const int subIdx1 = refElement.subEntity(intersection.indexInInside(), 1, /*idx*/1, dimgrid);
          const auto vertex1 = intersection.inside().template subEntity<2>( subIdx1 );
 
          std::array<unsigned int, 2> vertexIndices {{
            this->insertionIndex( vertex0 ),
            this->insertionIndex( vertex1 )
          }};
 
          // sort the indices
          if ( vertexIndices[0] > vertexIndices[1] )
            std::swap( vertexIndices[0], vertexIndices[1] );
 
          const auto& it = boundarySegmentIndices_.find( vertexIndices );
          return (it != boundarySegmentIndices_.end());
        }
 
        void insertElement(const GeometryType& type,
                           const std::vector<unsigned int>& vertices) override {
 
            assert(type.isTriangle());
 
            EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
            newElement.vertex_[0] = this->vertexArray_[vertices[0]];
            newElement.vertex_[1] = this->vertexArray_[vertices[1]];
            newElement.vertex_[2] = this->vertexArray_[vertices[2]];
 
            std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
        }
 
        void insertElement(const GeometryType& type,
                           const std::vector<unsigned int>& vertices,
                           std::function<FieldVector<ctype,dimworld>(FieldVector<ctype,dimgrid>)> elementParametrization)
        override
        {
            assert(type.isTriangle());
            EntityImp<dimgrid> newElement(/*level=*/0, this->grid_->getNextFreeId());
            newElement.vertex_[0] = this->vertexArray_[vertices[0]];
            newElement.vertex_[1] = this->vertexArray_[vertices[1]];
            newElement.vertex_[2] = this->vertexArray_[vertices[2]];
            newElement.elementParametrization_ = elementParametrization;
            std::get<dimgrid>(this->grid_->entityImps_[0]).push_back(newElement);
        }
 
        GridPtrType createGrid() override
        {
            // Prevent a crash when this method is called twice in a row
            // You never know who may do this...
            if (this->grid_==nullptr)
                return nullptr;
 
            // ////////////////////////////////////////////////////
            //   Create the edges
            // ////////////////////////////////////////////////////
 
            // For fast retrieval: a map from pairs of vertices to the edge that connects them
            std::unordered_map<std::pair<const EntityImp<0>*, const EntityImp<0>*>, EntityImp<1>*,
                               HashPair<const EntityImp<0>*>> edgeMap;
            for (auto& element : std::get<dimgrid>(this->grid_->entityImps_[0]))
            {
                const auto refElement = ReferenceElements<ctype, dimgrid>::general(element.type());
 
                // Loop over all edges of this element
                for (std::size_t i=0; i<element.facet_.size(); ++i)
                {
                    // get pointer to the edge (insert edge if it's not in the entity list yet)
                    auto edge = [&](){
                      // Get two vertices of the potential edge
                      auto v0 = element.vertex_[refElement.subEntity(i, 1, 0, 2)];
                      auto v1 = element.vertex_[refElement.subEntity(i, 1, 1, 2)];
                      // sort pointers
                      if ( v0 > v1 )
                        std::swap( v0, v1 );
 
                      // see if the edge was already inserted
                      // the pointer pair hash is symmetric so we don't have to check for pair(v1,v0)
                      auto e = edgeMap.find(std::make_pair(v0, v1));
                      if (e != edgeMap.end())
                        return e->second;
 
                      // edge was not inserted yet: insert the edge now
                      std::get<1>(this->grid_->entityImps_[0]).emplace_back(v0, v1, /*level=*/0, this->grid_->getNextFreeId());
                      // insert it into the map of inserted edges
                      auto newEdge = &*std::get<1>(this->grid_->entityImps_[0]).rbegin();
                      edgeMap.insert(std::make_pair(std::make_pair(v0, v1), newEdge));
                      return newEdge;
                    }();
 
                    // make element know about the edge
                    element.facet_[i] = edge;
 
                    // make edge know about the element
                    edge->elements_.push_back(&element);
                }
            }
 
            // Create the index sets
            this->grid_->setIndices();
 
 
            // ////////////////////////////////////////////////
            //   Set the boundary ids
            // ////////////////////////////////////////////////
 
            for (auto& facet : std::get<1>(this->grid_->entityImps_[0]))
            {
                if (facet.elements_.size() == 1) // if boundary facet
                {
                    std::array<unsigned int, 2> vertexIndices {{ facet.vertex_[0]->leafIndex_, facet.vertex_[1]->leafIndex_ }};
                    // sort the indices
                    if ( vertexIndices[0] > vertexIndices[1] )
                      std::swap( vertexIndices[0], vertexIndices[1] );
 
                    const auto it = boundarySegmentIndices_.find( vertexIndices );
                    if (it != boundarySegmentIndices_.end())
                        facet.boundarySegmentIndex_ = it->second;
                    else
                        facet.boundarySegmentIndex_ = this->boundarySegmentCounter_++;
                }
            }
 
            // ////////////////////////////////////////////////
            //   Hand over the new grid
            // ////////////////////////////////////////////////
 
            Dune::FoamGrid<dimgrid, dimworld, ctype>* tmp = this->grid_;
            tmp->numBoundarySegments_ = this->boundarySegmentCounter_;
            this->grid_ = nullptr;
            return GridPtrType(tmp);
        }
 
      private:
        template<class T, class U = T>
        struct HashPair {
          std::size_t operator() (const std::pair<T, U>& a) const {
            std::size_t seed = 0;
            hash_combine(seed, a.first);
            hash_combine(seed, a.second);
            return seed;
          }
        };
 
        struct HashUIntArray {
          std::size_t operator() (const std::array<unsigned int, 2>& a) const {
            return hash_range(a.begin(), a.end());
          }
        };
 
        std::unordered_map<std::array<unsigned int, 2>, unsigned int, HashUIntArray> boundarySegmentIndices_;
    };
 
} // end namespace Dune
 
#endif