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
 
 
 
#endif