boundaryLocator.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2017-2020 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef BOUNDARYLOCATOR_HH
#define BOUNDARYLOCATOR_HH
 
#include <memory>
#include <optional>
#include <tuple>
 
#include "dune/common/fvector.hh"
#include "dune/geometry/type.hh"
 
#include "fem/boundaryFace.hh"
#include "utilities/timing.hh"
 
namespace Kaskade
{
  // forward declaration for PIMPL idiome
  namespace BoundaryLocatorDetail
  {
    template <class GridView, class Function, int dimw>
    class SpatialIndex;    
    
    // Dummy implementation for FE function, just to be syntactically correct
    // TODO: replace this by void type and use constexpr if.
    struct ZeroDisplacement
    {
      template <class Cell>
      int order(Cell const&) const { return 0; }
 
      template <class Cell, class B>
      B value(Cell, B b) const { return B(0.0); }
 
      template <class Cell, class B>
      Dune::FieldMatrix<typename B::field_type,B::dimension,B::dimension> derivative(Cell, B) const { return 0.0; }
    };
 
  }
  
  // ----------------------------------------------------------------------------------------------
 
 
 
  template <class GridView, class Function=BoundaryLocatorDetail::ZeroDisplacement, int dimw=GridView::Grid::dimension>
  class BoundaryLocator
  {
    using Index = BoundaryLocatorDetail::SpatialIndex<GridView,Function,dimw>;
    using Grid = typename GridView::Grid;
 
  public:
    static int const dimension = GridView::Grid::dimension;
    static int const dimension_world = dimw;
    using ctype    = typename GridView::ctype;
    using Position = Dune::FieldVector<ctype,dimension_world>;
    using LocalPosition = Dune::FieldVector<ctype,dimension_world-1>;    
    using Face = std::conditional_t<dimension==dimension_world, 
                                    typename GridView::Intersection,
                                    typename GridView::template Codim<0>::Entity>;
    using DisplacedFace = BoundaryFace<Grid,Face,Function,dimw>;
 
    
    BoundaryLocator(GridView const& gridView);
 
    BoundaryLocator(GridView const& gridView, std::vector<bool> const& boundarySegments);
 
    template <class... Args>
    BoundaryLocator(GridView const& gridView, Args&&... args)
    : BoundaryLocator(gridView)
    {
      f = std::make_unique<Function>(std::forward<Args>(args)...);
    }
 
    // we need this because the default destructor doesn't know about the size of SpatialIndex when deallocating via unique_ptr
    ~BoundaryLocator();
    
    template <class Update>
    void updateDisplacement(Update const& update)
    {
      // Extract the faces from the spatial index. This is much faster than running through the whole
      // grid (O(n^(1-1/d)) vs O(n)). Moreover, faces on periodic boundaries and on not considered
      // boundary segments are already filtered out.
      auto& timer = Timings::instance();
      timer.start("boundary faces from rtree");
      auto fs = faces();
      timer.stop("boundary faces from rtree");
 
      if (f)
        *f = update;              // register new displacement
      else
        f = std::make_unique<Function>(update);
 
      enterFaces(fs);                   // recreate the spatial tree with newly displaced faces
    }
 
 
    std::optional<std::tuple<typename BoundaryLocator<GridView,Function,dimw>::DisplacedFace,
                             typename BoundaryLocator<GridView,Function,dimw>::LocalPosition,
                             typename GridView::ctype>>
    byLineSegment(Position const& from, Position const& to, ctype minAngle=-0.5) const;
    
    std::optional<std::tuple<typename BoundaryLocator<GridView,Function,dimw>::DisplacedFace,
                             typename BoundaryLocator<GridView,Function,dimw>::LocalPosition,
                             typename GridView::ctype>>
    byRay(Position const& from, Position const& direction, ctype minAngle=-0.5) const
    {
      return byLineSegment(from,from+(diam/direction.two_norm())*direction, minAngle);
    }
    
    ctype diameter() const 
    {
      return diam;
    }
    
    Position global(Face const& face, LocalPosition const& xi) const;
    
    Position unitOuterNormal(Face const& face, LocalPosition const& xi) const;
    
    Function const* displacement() const
    {
      return f.get();
    }
    
    std::vector<Face> faces() const;
 
    std::vector<DisplacedFace> displacedFaces() const;
 
  private:
    std::unique_ptr<BoundaryLocatorDetail::SpatialIndex<GridView,Function,dimension_world>> spatialIndex;
    std::unique_ptr<Function> f;
    ctype diam;
    
    using IntersectionSet = std::vector<std::tuple<BoundaryFace<Grid,Face,Function>,
                                                   LocalPosition,ctype,ctype>>;
 
    void init(GridView const& gridView, std::vector<bool> const& boundarySegments);
    
    IntersectionSet getIntersectingFaces(Position const& from, Position const& to, ctype minAngle) const;
    
    void enterFaces(std::vector<Face> const& fs);
  };
}
 
#endif