resultfunction.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers mueller@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <type_traits>
 
#include <dune/vtk/function.hh>
#include <dune/vtk/vtkwriter.hh>
 
#include <ikarus/finiteelements/ferequirements.hh>
#include <ikarus/utils/concepts.hh>
 
namespace Ikarus {
namespace Impl {
  struct DefaultUserFunction
  {
  };
} // namespace Impl
 
template <typename AS, template <typename, int, int> class RT, typename UserFunction = Impl::DefaultUserFunction>
requires(Concepts::FlatAssembler<AS> and Concepts::ResultType<RT>)
class ResultFunction : public Dune::VTKFunction<typename AS::GridView>
{
public:
  using Assembler         = AS;
  using GridView          = typename Assembler::GridView;
  using FERequirementType = typename Assembler::FERequirement;
  using FEContainer       = typename Assembler::FEContainer;
  using FiniteElement     = typename std::remove_cvref_t<FEContainer>::value_type;
 
  using ctype                  = typename GridView::ctype;
  constexpr static int griddim = GridView::dimension;
  using Entity                 = typename GridView::template Codim<0>::Entity;
 
  double evaluate(int comp, const Entity& e, const Dune::FieldVector<ctype, griddim>& local) const override {
    const auto index = gridView().indexSet().index(e);
    return evaluateComponent(index, local, comp);
  }
 
  [[nodiscard]] int ncomps() const override {
    if constexpr (std::is_same_v<UserFunction, Impl::DefaultUserFunction>) {
      Dune::FieldVector<ctype, griddim> val(0.0);
      auto sigma = finiteElements().at(0).template calculateAt<RT>(requirement(), val).asVec();
      return static_cast<int>(sigma.size());
    } else
      return userFunction_.ncomps();
  }
 
  [[nodiscard]] constexpr std::string name() const override {
    if constexpr (requires { userFunction_.name(); })
      return userFunction_.name();
    else
      return toString<RT>();
  }
 
  Dune::VTK::Precision precision() const override { return prec_; }
 
  ResultFunction(std::shared_ptr<Assembler> assembler, Dune::VTK::Precision prec = Dune::VTK::Precision::float64)
      : assembler_(assembler),
        prec_{prec},
        userFunction_{UserFunction{}} {}
 
private:
  double evaluateComponent(int eleID, const Dune::FieldVector<ctype, griddim>& local, int comp) const {
    auto result = finiteElements().at(eleID).template calculateAt<RT>(requirement(), local).asVec();
 
    if constexpr (!std::is_same_v<UserFunction, Impl::DefaultUserFunction>)
      return userFunction_(result, comp);
    else
      return result(comp);
  }
 
  const FEContainer& finiteElements() const { return assembler_->finiteElements(); }
  const FERequirementType& requirement() const { return assembler_->requirement(); }
  const GridView& gridView() const { return assembler_->gridView(); }
 
  std::shared_ptr<Assembler> assembler_;
 
  Dune::VTK::Precision prec_;
  [[no_unique_address]] std::string name_{};
  UserFunction userFunction_;
};
 
template <template <typename, int, int> class RT, typename UserFunction = Impl::DefaultUserFunction, typename AS>
auto makeResultFunction(std::shared_ptr<AS> assembler, Dune::VTK::Precision prec = Dune::VTK::Precision::float64) {
  return std::make_shared<ResultFunction<AS, RT, UserFunction>>(assembler, prec);
}
 
template <template <typename, int, int> class RT, typename UserFunction = Impl::DefaultUserFunction, typename AS>
auto makeResultVtkFunction(std::shared_ptr<AS> assembler) {
  return Dune::Vtk::Function<typename AS::GridView>(std::make_shared<ResultFunction<AS, RT, UserFunction>>(assembler));
}
 
} // namespace Ikarus