geometryoperators.hh

Go to the documentation of this file.

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
 
// 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 DUNE_FUFEM_FORMS_GEOMETRYOPERATORS_HH
#define DUNE_FUFEM_FORMS_GEOMETRYOPERATORS_HH
 
#include <type_traits>
#include <utility>
 
#include <dune/common/typetraits.hh>
 
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
 
#include <dune/fufem/forms/baseclass.hh>
#include <dune/fufem/forms/tensors.hh>
#include <dune/fufem/forms/shapefunctioncache.hh>
 
 
 
namespace Dune::Fufem::Forms {
 
 
 
  template<class GV>
  class FaceNormalOperator : public MultilinearOperator<0>
  {
  public:
 
    using GridView = GV;
    using Element = typename GridView::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;
    using Range = typename Intersection::GlobalCoordinate;
 
    FaceNormalOperator()
      : cacheId_(*this)
    {}
 
    FaceNormalOperator(const GridView& gridView)
      : cacheId_(*this)
    {}
 
    class LocalOperator
    {
      static constexpr int dimension = FaceNormalOperator::Element::dimension;
      using ValueCache = SimpleCache<double, dimension, FaceNormalOperator::Range>;
    public:
 
      using Element = typename GridView::template Codim<0>::Entity;
      using Intersection = typename GridView::Intersection;
      using CacheManager = Dune::Fufem::Forms::CacheManager<typename Element::Geometry::ctype, Element::Geometry::mydimension>;
      using Range = typename Intersection::GlobalCoordinate;
 
      LocalOperator(UniqueCacheId cacheId)
        : cacheId_(cacheId)
      {}
 
      auto quadratureRuleKey() const
      {
        return quadratureRuleKey_;
      }
 
      template<class... LV>
      void registerLocalViews(const LV&... lvs)
      {}
 
      template<class... LV>
      void registerOutsideLocalViews(const LV&... lvs)
      {}
 
      void registerCaches(CacheManager& cacheManager)
      {
        cacheIndex_ = cacheManager.registerCache(cacheId_, ValueCache(true));
      }
 
      void registerOutsideCaches(CacheManager& cacheManager)
      {}
 
      template<class Dummy=void>
      void bind(const Element&)
      {
        static_assert(Dune::AlwaysFalse<Dummy>::value, "FaceNormalOperator can only be used for forms integrated over intersections");
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        intersection_ = &intersection;
        intersectionIsOriented_ = (intersection_->inside() == element);
        quadratureRuleKey_ = QuadratureRuleKey(element.type(), 0);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {
        valueCache_ = &cacheManager.template getCache<ValueCache>(cacheIndex_);
        if (valueCache_->isEmpty())
        {
          auto& values = valueCache_->getValues();
          const auto& rule = valueCache_->rule();
          if (intersectionIsOriented_)
            for (auto k : Dune::range(rule.size()))
            {
              auto xInIntersection = intersection_->geometryInInside().local(rule[k].position());
              values[k] = intersection_->unitOuterNormal(xInIntersection);
            }
          else
            for (auto k : Dune::range(rule.size()))
            {
              auto xInIntersection = intersection_->geometryInOutside().local(rule[k].position());
              values[k] = -intersection_->unitOuterNormal(xInIntersection);
            }
          valueCache_->setEmpty(false);
        }
      }
 
      auto operator()(std::size_t index) const
      {
        return Impl::Tensor::RankZeroTensor(valueCache_->getValues()[index]);
      }
 
    private:
      QuadratureRuleKey quadratureRuleKey_;
      const Intersection* intersection_ = nullptr;
      bool intersectionIsOriented_ = false;
      UniqueCacheId cacheId_;
      std::size_t cacheIndex_ = 0;
      ValueCache* valueCache_;
    };
 
    friend LocalOperator localOperator(const FaceNormalOperator& faceNormalOperator)
    {
      return LocalOperator(faceNormalOperator.cacheId_);
    }
 
    auto basis() const
    {
      return std::tuple<>();
    }
 
    auto treePath() const
    {
      return std::tuple<>();
    }
 
  private:
    UniqueCacheId cacheId_;
  };
 
 
 
  template<class GV>
  class MeshSizeOperator : public MultilinearOperator<0>
  {
  public:
 
    using GridView = GV;
    using Element = typename GridView::template Codim<0>::Entity;
    using Intersection = typename GridView::Intersection;
    using Range = double;
 
    MeshSizeOperator() = default;
 
    MeshSizeOperator(const GridView& gridView)
    {}
 
    class LocalOperator
    {
      static constexpr int dimension = MeshSizeOperator::Element::dimension;
    public:
 
      using Element = typename GridView::template Codim<0>::Entity;
      using Intersection = typename GridView::Intersection;
      using CacheManager = Dune::Fufem::Forms::CacheManager<typename Element::Geometry::ctype, Element::Geometry::mydimension>;
      using Range = double;
 
      auto quadratureRuleKey() const
      {
        return quadratureRuleKey_;
      }
 
      template<class... LV>
      void registerLocalViews(const LV&... lvs)
      {}
 
      template<class... LV>
      void registerOutsideLocalViews(const LV&... lvs)
      {}
 
      void registerCaches(CacheManager& cacheManager)
      {}
 
      void registerOutsideCaches(CacheManager& cacheManager)
      {}
 
      void bind(const Element& element)
      {
        const auto& geometry = element.geometry();
        auto&& re = Dune::referenceElement(geometry);
        h_ = std::pow(geometry.integrationElement(re.position(0, 0)), 1./dimension);
        quadratureRuleKey_ = QuadratureRuleKey(element.type(), 0);
      }
 
      void bind(const Intersection& intersection, const Element& element, const Element& otherElement)
      {
        const auto& geometry = intersection.geometry();
        auto&& re = Dune::referenceElement(geometry);
        h_ = std::pow(geometry.integrationElement(re.position(0, 0)), 1./(dimension-1));
        quadratureRuleKey_ = QuadratureRuleKey(element.type(), 0);
      }
 
      template<class... OutsideCacheManager>
      void bindToCaches(CacheManager& cacheManager, OutsideCacheManager&... outsideCacheManager)
      {}
 
      auto operator()(std::size_t index) const
      {
        return Impl::Tensor::RankZeroTensor(h_);
      }
 
    private:
      QuadratureRuleKey quadratureRuleKey_;
      double h_ = 0;
    };
 
    friend LocalOperator localOperator(const MeshSizeOperator& meshSizeOperator)
    {
      return LocalOperator();
    }
 
    auto basis() const
    {
      return std::tuple<>();
    }
 
    auto treePath() const
    {
      return std::tuple<>();
    }
  };
 
 
 
} // namespace Dune::Fufem::Forms
 
 
 
#endif // DUNE_FUFEM_FORMS_GEOMETRYOPERATORS_HH