geometryinancestor.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUNCTIONS_COMMON_GEOMETRYINANCESTOR_HH
#define DUNE_FUNCTIONS_COMMON_GEOMETRYINANCESTOR_HH
 
#include <vector>
 
#include <dune/common/iteratorrange.hh>
 
namespace Dune::Functions {
 
 
 
template<class Element>
class GeometryInAncestor
{
  using GeometryInFather = typename Element::LocalGeometry;
 
  template<class V>
  static auto reverse_view(V&& v)
  {
    return Dune::IteratorRange(v.rbegin(), v.rend());
  }
 
public:
 
  using LocalCoordinate = typename Element::Geometry::LocalCoordinate;
  //
  using GlobalCoordinate = typename Element::Geometry::LocalCoordinate;
 
  GeometryInAncestor() = default;
  GeometryInAncestor(GeometryInAncestor&&) = delete;
 
  GeometryInAncestor(const GeometryInAncestor& other)
    : storedFineElement_(other.storedFineElement_)
    , storedCoarseElement_(other.storedCoarseElement_)
    , fineElementPtr_(&storedFineElement_)
    , coarseElementPtr_(&storedFineElement_)
    , geometryInFathersVector_(other.geometryInFathersVector_)
  {
    // If other points to external fine element,
    // we point to the same one.
    if (other.fineElementPtr_ != & other.storedFineElement_)
      fineElementPtr_ = other.fineElementPtr_;
    // If the other coarse element is not the stored one,
    // then it's the fine one.
    if (other.coarseElementPtr_ != & other.storedCoarseElement_)
      coarseElementPtr_ = fineElementPtr_;
  }
 
  GeometryInAncestor(const GeometryInAncestor& other, const Element& fineElement)
    : storedFineElement_(other.storedFineElement_)
    , storedCoarseElement_(other.storedCoarseElement_)
    , fineElementPtr_(&fineElement)
    , coarseElementPtr_(&storedFineElement_)
    , geometryInFathersVector_(other.geometryInFathersVector_)
  {
    // If the other coarse element is not the stored one,
    // then it's the fine one.
    if (other.coarseElementPtr_ != & other.storedCoarseElement_)
      coarseElementPtr_ = fineElementPtr_;
  }
 
  template<class F>
  requires (std::is_invocable_r_v<bool,F,Element>)
  const Element& bind(const Element& fineElement, F&& includeFather)
  {
    fineElementPtr_ = &fineElement;
    geometryInFathersVector_.clear();
    coarseElementPtr_ = fineElementPtr_;
    while (includeFather(*coarseElementPtr_))
    {
      geometryInFathersVector_.push_back(coarseElementPtr_->geometryInFather());
      storedCoarseElement_ = coarseElementPtr_->father();
      coarseElementPtr_ = &storedCoarseElement_;
    }
    return coarseElement();
  }
 
  template<class F>
  requires (std::is_invocable_r_v<bool,F,Element>)
  const Element& bind(Element&& fineElement, F&& includeFather)
  {
    storedFineElement_ = fineElement;
    bind(storedFineElement_, includeFather);
  }
 
  const Element& fineElement() const
  {
    return *fineElementPtr_;
  }
 
  const Element& coarseElement() const
  {
    return *coarseElementPtr_;
  }
 
  GlobalCoordinate global(LocalCoordinate x) const
  {
    for (const auto& g : geometryInFathersVector_)
      x = g.global(x);
    return x;
  }
 
  LocalCoordinate local(GlobalCoordinate x) const
  {
    // Here we can be instead use std::ranges::reverse_view once
    // all supported compilers provide this (gcc>=10, clang>=16).
    for (const auto& g : reverse_view(geometryInFathersVector_))
      x = g.local(x);
    return x;
  }
 
private:
  Element storedFineElement_;
  Element storedCoarseElement_;
  const Element* fineElementPtr_ = nullptr;
  const Element* coarseElementPtr_ = nullptr;
  std::vector<GeometryInFather> geometryInFathersVector_;
};
 
 
 
} // namespace Dune::Functions
 
 
 
 
#endif // DUNE_FUNCTIONS_COMMON_GEOMETRYINANCESTOR_HH