ARTICLE AD BOX
I need to extend the Function class coming from the deal.II library.
The current approach is to have a RealFunction extend Function and then to have things like ConstantRealFunction, SineFunction, etc. extend RealFunction.
template <int dim> class RealFunction : public dealii::Function<dim> { public: RealFunction() : dealii::Function<dim>() {} virtual double value(const dealii::Point<dim> &p, const unsigned int component = 0) const override = 0; template <typename Number> Number value(const dealii::Point<dim, Number> &p, const unsigned int component = 0) const { throw std::logic_error("Class not extended correctly."); } }; template <int dim> class ConstantRealFunction : public RealFunction<dim> { private: double val; public: ConstantRealFunction(double v) : RealFunction<dim>(), val(v) {} virtual double value(const dealii::Point<dim> &p, const unsigned int component = 0) const override { return value<double>(p, component); } template <typename Number> Number value(const dealii::Point<dim, Number> & /*p*/, const unsigned int /*component*/ = 0) const { return Number(val); } };The main reason behind this choice of inheritance structure is to enable the possibility to evaluate the same functions both on a single coordinate (Point<dim, double>) or on multiple coordinates (Point<dim, VectorizedArray<double>>) at the same time with SIMD instructions without having to repeat two times the same code.
This idea is a reelaboration of the implementation of the Coefficient class from the deal.II step-37 tutorial, from which I copy here the relevant part.
// Code from deal.II step-37 tutorial template <int dim> class Coefficient : public Function<dim> { public: virtual double value(const Point<dim> &p, const unsigned int component = 0) const override; template <typename number> number value(const Point<dim, number> &p, const unsigned int component = 0) const; }; template <int dim> template <typename number> number Coefficient<dim>::value(const Point<dim, number> &p, const unsigned int /*component*/) const { return 1. / (0.05 + 2. * p.square()); } template <int dim> double Coefficient<dim>::value(const Point<dim> &p, const unsigned int component) const { return value<double>(p, component); }From my understanding of C++ inheritance, considering an instance of the Coefficient implementation from the tutorial, it is true that
a call to value passing an instance of Point<dim> would result in a call to the value method that does not have the template <typename number> which, in turn, calls the other one; a call to value passing an instance of Point<dim, VectorizedArray<double>> would result in a call directly to the value with template <typename number>.Now consider the case with RealFunction, ConstantRealFunction and these definitions
template <int dim> using DirichletBoundary = RealFunction<dim>; template <typename T> using Boundaries = std::unordered_map<types::boundary_id, T>; // boundary_id is just unsigned int template <int dim> using DirichletBoundaries = Boundaries<const DirichletBoundary<dim> *>; template <int dim> struct ProblemData { DirichletBoundaries<dim> dirichlet_boundaries; // ... } template <int dim> ProblemData<dim> example_problem_data() { ProblemData<dim> data; data.dirichlet_boundaries = DirichletBoundaries<dim>(); ConstantRealFunction<dim> crf(1.0); data.dirichlet_boundaries[0] = &crf; // ... return data; }Then, if we iterate over the Dirichlet boundaries like
for (const auto &[boundary_id, function] : this->problem.dirichlet_boundaries) { VectorTools::interpolate_boundary_values(mapping, dof_handler, boundary_id, *function, constraints); }even if this seems to be coherent with some working examples on the internet, whenever the VectorTools::interpolate_boundary_values tries to call the value of the function, the program throws a Class not extended correctly as if the value that gets called is not getting correctly overridden.
What am I doing wrong?
