deformationfunction.hh

Go to the documentation of this file.

// 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 DEFORMATION_FUNCTION_HH
#define DEFORMATION_FUNCTION_HH
 
#include <cstddef>
#include <iostream>
#include <limits>
#include <memory>
#include <vector>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/indices.hh>
 
#include <dune/istl/bvector.hh>
 
#include <dune/grid/geometrygrid/coordfunction.hh>
 
#include <dune/geometry/referenceelements.hh>
 
#include <dune/localfunctions/common/localkey.hh>
#include <dune/localfunctions/lagrange/lagrangelfecache.hh>
 
#include <dune/functions/functionspacebases/defaultglobalbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
 
 
 
namespace Dune::Fufem::Experimental {
 
 
 
template<class HG, std::size_t deformedDim, class ct = double>
class DeformationFunction
  : public Dune::DiscreteCoordFunction<ct, deformedDim, Dune::Fufem::Experimental::DeformationFunction<HG, deformedDim, ct>>
{
  static constexpr int dim = HG::dimension;
 
  using HostLeafGridView = typename HG::LeafGridView;
  using PreBasis = Dune::Functions::PowerPreBasis<
    Dune::Functions::BasisFactory::BlockedInterleaved,
    Dune::Functions::LagrangePreBasis<HostLeafGridView, 1>,
    deformedDim>;
 
public:
 
  using HostGrid = HG;
  using Domain = typename HostGrid::template Codim<dim>::Geometry::GlobalCoordinate;
  using ctype = ct;
  using Range = Dune::FieldVector<ct, deformedDim>;
  using Basis = Dune::Functions::DefaultGlobalBasis<PreBasis>;
  using CoefficientVector = Dune::BlockVector<Range>;
 
private:
 
  // Generic evaluate implementation. Since this is templated
  // is virtual, we have to redirect the virtual interface methods
  // here manually.
  template<unsigned int codim>
  void evaluateImpl(const typename HostGrid::template Codim<codim>::Entity& entity, unsigned int corner, Range& y) const
  {
    if (hostGrid_->leafIndexSet().contains(entity))
    {
      const auto& leafIndexSet = hostGrid_->leafIndexSet();
      y  = coefficients_[leafIndexSet.subIndex(entity, corner, dim)];
    }
    else
    {
      auto level = entity.level();
      const auto& levelIndexSet = hostGrid_->levelIndexSet(level);
      y  = coefficients_[leafIndices_[level][levelIndexSet.subIndex(entity, corner, dim)]];
    }
  }
 
public:
 
  DeformationFunction(const HostGrid& hostGrid, Dune::index_constant<deformedDim> deformedDimConstant = {})
    : hostGrid_(&hostGrid)
    , coefficients_(hostGrid_->leafGridView().size(dim), 0)
    , basis_(hostGrid_->leafGridView())
  {
    setupLeafIndices();
  }
 
  template<class F>
  DeformationFunction(const HostGrid& hostGrid, const F& f, Dune::index_constant<deformedDim> deformedDimConstant = {})
    : DeformationFunction(hostGrid)
  {
    interpolate(f);
  }
 
  virtual ~DeformationFunction()
  {}
 
  virtual void evaluate(const typename HostGrid::template Codim<0>::Entity& entity, unsigned int corner, Range& y) const
  {
    evaluateImpl<0>(entity, corner, y);
  }
 
  virtual void evaluate(const typename HostGrid::template Codim<dim>::Entity& entity, unsigned int corner, Range& y) const
  {
    evaluateImpl<dim>(entity, corner, y);
  }
 
  void adapt()
  {
    DUNE_THROW(Dune::NotImplemented, "DeformationFunction::adapt() is not implemented");
  }
 
  const Basis& basis() const
  {
    return basis_;
  }
 
  CoefficientVector& coefficients()
  {
    return coefficients_;
  }
 
  const CoefficientVector& coefficients() const
  {
    return coefficients_;
  }
 
  template<class F>
  void interpolate(const F& f)
  {
    Dune::Functions::interpolate(basis(), coefficients(), f);
  }
 
protected:
 
  void setupLeafIndices()
  {
    constexpr std::size_t unsetIndex = std::numeric_limits<std::size_t>::max();
 
    const auto& leafIndexSet = hostGrid_->leafIndexSet();
    const auto& idSet = hostGrid_->localIdSet();
 
    // Progress fine to coarse to compute all level-index to leaf-index map
    leafIndices_.resize(hostGrid_->maxLevel());
    for (int level=hostGrid_->maxLevel()-1; level>=0; --level)
    {
      const auto& levelGridView = hostGrid_->levelGridView(level);
      const auto& levelIndexSet = hostGrid_->levelIndexSet(level);
      const auto& fineLevelIndexSet = hostGrid_->levelIndexSet(level+1);
 
      leafIndices_[level].clear();
      leafIndices_[level].resize(levelIndexSet.size(dim), unsetIndex);
 
      for (const auto& element: elements(levelGridView))
      {
        const auto& re = Dune::referenceElement(element);
        // If the element is leaf, we can directly query the leaf-indexset
        if (element.isLeaf())
        {
          for (auto i : Dune::range(re.size(dim)))
            leafIndices_[level][levelIndexSet.subIndex(element, i, dim)] = leafIndexSet.subIndex(element, i, dim);
          continue;
        }
        else
        {
          // Loop over all direct children
          for(const auto& childElement : descendantElements(element, level+1))
          {
            if (childElement.level() == level)
              DUNE_THROW(Dune::RangeError, "HierarchicIterator visited element on the same level");
            // Loop over all vertices of the element
            const auto& childRe = Dune::referenceElement(childElement);
            for (auto i : Dune::range(re.size(dim)))
            {
              auto levelIndex = levelIndexSet.subIndex(element, i, dim);
              if (leafIndices_[level][levelIndex] == unsetIndex)
              {
                // Find child vertex with the same id
                auto id = idSet.subId(element, i, dim);
                for (auto j : Dune::range(childRe.size(dim)))
                {
                  if (idSet.subId(element, j, dim) == id)
                  {
                    if (childElement.isLeaf())
                      leafIndices_[level][levelIndex] = leafIndexSet.subIndex(childElement, j, dim);
                    else
                    {
                      auto fineLevelIndex = levelIndexSet.subIndex(childElement, j, dim);
                      leafIndices_[level][levelIndex] = leafIndices_[level+1][fineLevelIndex];
                    }
                    continue;
                  }
                }
              }
            }
          }
        }
      }
    }
  }
 
  std::vector<std::vector<std::size_t> > leafIndices_;
  const HostGrid* hostGrid_;
  Basis basis_;
  CoefficientVector coefficients_;
};
 
 
 
template<class HostGrid, class F>
DeformationFunction(const HostGrid& hostGrid, const F& f)
  -> DeformationFunction<HostGrid, std::decay_t<std::invoke_result_t<F, typename HostGrid::template Codim<HostGrid::dimension>::Geometry::GlobalCoordinate>>::size()>;
 
 
 
} // end namespace Dune::Fufem::Experimental
 
 
 
template <class GridView, class CoefficientVectorType>
class
DeformationFunction  :
 public Dune::DiscreteCoordFunction<typename GridView::ctype, GridView::dimension,
                                        DeformationFunction<GridView,CoefficientVectorType> >
{
    enum {dim = GridView::dimension};
 
    typedef typename GridView::ctype ctype;
    typedef DeformationFunction<GridView,CoefficientVectorType> This;
    typedef Dune::DiscreteCoordFunction<ctype, dim, This > Base;
public:
    typedef typename Base::RangeVector RangeVector;
 
    [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufen::Experimental::DeformationFunction instead.")]]
    DeformationFunction(const GridView& gridView)
    {
        gridView_ = std::unique_ptr<GridView>(new GridView(gridView));
        deformation_.resize(gridView_->size(dim));
        deformation_ = 0;
    }
 
 
    [[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufen::Experimental::DeformationFunction instead.")]]
    DeformationFunction(const GridView& gridView,
                        const CoefficientVectorType& deformation)
        : deformation_(deformation)
    {
        gridView_ = std::unique_ptr<GridView>(new GridView(gridView));
 
        if(gridView_->size(dim) != deformation_.size())
            DUNE_THROW(Dune::Exception,"The deformation coefficient vector doesn't match the gridview!");
    }
 
    virtual ~DeformationFunction() {}
 
    virtual void setGridView(const GridView& gridView) {
        gridView_ = std::unique_ptr<GridView>(new GridView(gridView));
    }
 
    void setDeformation(const CoefficientVectorType& deformation)
    {
        if (deformation_.size() != deformation.size())
            std::cout<<"Warning chanign deformation size "<<deformation_.size()<<" to "<<deformation.size()<<std::endl;
        deformation_ = deformation;
    }
 
    virtual void evaluate ( const typename GridView::template Codim<dim>::Entity& hostVertex, unsigned int corner,
                    Dune::FieldVector<ctype,dim> &y ) const
    {
        const typename GridView::IndexSet& indexSet = gridView_->indexSet();
        int idx = indexSet.index(hostVertex);
 
        y  = hostVertex.geometry().corner(0) + deformation_[idx];
    }
 
    virtual void evaluate (const typename GridView::template Codim<0>::Entity& hostEntity, unsigned int corner,
                    typename Base::RangeVector& y ) const
    {
        const typename GridView::IndexSet& indexSet = gridView_->indexSet();
        int idx = indexSet.subIndex(hostEntity, corner,dim);
        y  = hostEntity.geometry().corner(corner) + deformation_[idx];
    }
 
    void adapt() {
        if ((int) deformation_.size() != gridView_->size(dim)) {
            std::cout<<"Warning removing all data from DeformationFunction! \n";
            deformation_.resize(gridView_->size(dim));
            deformation_ = 0;
        }
    }
 
    const CoefficientVectorType& getDeformation() {return deformation_;}
protected:
    std::unique_ptr<GridView> gridView_;
 
    CoefficientVectorType deformation_;
};
 
 
template <class GridView, class CoefficientVectorType>
class
[[deprecated("This class is deprecated and will be removed after 2.11. Use Dune::Fufen::Experimental::DeformationFunction instead.")]]
DeformationHierarchyFunction  :
 public DeformationFunction<GridView,CoefficientVectorType>
{
    enum {dim = GridView::dimension};
 
    typedef typename GridView::ctype ctype;
    typedef typename CoefficientVectorType::field_type field_type;
    typedef DeformationFunction<GridView,CoefficientVectorType> Base;
public:
 
    DeformationHierarchyFunction(const GridView& gridView) :
        Base(gridView)
    {
        setupIndexMapping();
    }
 
 
    DeformationHierarchyFunction(const GridView& gridView,
                        const CoefficientVectorType& deformation) :
        Base(gridView, deformation)
    {
        setupIndexMapping();
    }
 
    virtual void setGridView(const GridView& gridView) {
        Base::setGridView(gridView);
        setupIndexMapping();
    }
    // TODO Implement evaluate method for corners...
 
    void evaluate (const typename GridView::template Codim<0>::Entity& hostEntity, unsigned int corner,
                    typename Base::RangeVector& y ) const
    {
        if (hostEntity.isLeaf()) {
            Base::evaluate(hostEntity,corner,y);
            return;
        }
        const auto& indexSet = this->gridView_->grid().levelIndexSet(hostEntity.level());
 
        int idx = indexSet.subIndex(hostEntity, corner,dim);
        y  = hostEntity.geometry().corner(corner) + this->deformation_[fineIndex_[hostEntity.level()][idx]];
    }
 
private:
    std::vector<std::vector<int> > fineIndex_;
 
    void setupIndexMapping() {
 
        int maxLevel = this->gridView_->grid().maxLevel();
 
        // all elements are leaf
        if (maxLevel==0)
            return;
 
        fineIndex_.resize(maxLevel+1);
 
        const typename GridView::IndexSet& leafIndex = this->gridView_->indexSet();
 
        // A factory for the P1 shape functions
        typedef typename Dune::LagrangeLocalFiniteElementCache<ctype, field_type, dim, 1> P1FECache;
        typedef typename P1FECache::FiniteElementType FEType;
        P1FECache p1FECache;
 
        const typename GridView::Grid& grid = this->gridView_->grid();
        typedef typename GridView::Grid::LevelGridView::template Codim<0>::Iterator LevelElementIterator;
        typedef typename LevelElementIterator::Entity EntityType;
        typedef typename EntityType::HierarchicIterator HierarchicIterator;
 
        for (int i=maxLevel; i>=0; i--) {
            fineIndex_[i].resize(grid.size(i,dim),-1);
 
            const auto& levelIndex = grid.levelIndexSet(i);
            int fineLevel = (i==maxLevel) ? i : (i+1);
            const auto& fineIndex = grid.levelIndexSet(fineLevel);
 
            const auto& lv = grid.levelGridView(i);
            for (const auto& e: elements(lv)) {
 
                // all vertices of leaf elements are leaf
                if (e.isLeaf()) {
 
                    auto ref
                        = Dune::ReferenceElements<ctype, dim>::general(e.type());
                    for (int j=0; j<ref.size(dim); j++)
                        fineIndex_[i][levelIndex.subIndex(e,j,dim)] = leafIndex.subIndex(e,j,dim);
                    continue;
                }
 
                // Get local finite element
                const FEType& coarseBaseSet = p1FECache.get(e.type());
 
                const size_t numCoarseBaseFct = coarseBaseSet.localBasis().size();
 
                // preallocate vector for function evaluations
                std::vector<Dune::FieldVector<field_type,1> > values(numCoarseBaseFct);
 
                HierarchicIterator fIt    = e.hbegin(i+1);
                HierarchicIterator fEndIt = e.hend(i+1);
 
                for (; fIt != fEndIt; ++fIt) {
 
                    if (fIt->level()==i)
                        continue;
 
                    auto fineRefElement
                        = Dune::ReferenceElements<ctype, dim>::general(fIt->type());
 
                    const auto& fGeometryInFather = fIt->geometryInFather();
 
                    // Get local finite element
                    const FEType& fineBaseSet = p1FECache.get(fIt->type());
 
                    const size_t numFineBaseFct = fineBaseSet.localBasis().size();
 
                    for (size_t j=0; j<numFineBaseFct; j++)
                    {
                        //compute position of child element node in father element
                        const Dune::LocalKey& jLocalKey = fineBaseSet.localCoefficients().localKey(j);
 
                        int globalLeaf = fineIndex_[fineLevel][fineIndex.subIndex(*fIt, jLocalKey.subEntity(), dim)];
 
                        Dune::FieldVector<ctype, dim> fineBasePosition = fineRefElement.position(jLocalKey.subEntity(), dim);
                        Dune::FieldVector<ctype, dim> local = fGeometryInFather.global(fineBasePosition);
 
                        // Evaluate coarse grid base functions
                        coarseBaseSet.localBasis().evaluateFunction(local, values);
 
                        for (size_t k=0; k<numCoarseBaseFct; k++)
                        {
                            if (values[k] > 0.9999)
                            {
                                const Dune::LocalKey& coarseKey = coarseBaseSet.localCoefficients().localKey(k);
                                int globalCoarse = levelIndex.subIndex(e, coarseKey.subEntity(), dim);
                                fineIndex_[i][globalCoarse] = globalLeaf;
                            }
                        }
                    }
 
                }
            }
        }
    }
};
 
 
 
 
#endif