I am implementing a finite element code.
In finite element methods, we need an integrator and an interpolator. An integrator is an object that performs numerical integration on a geometrical object, for example, a quadrilateral, a triangle, etc. The integrator places multiple integration points, or abscissae, inside the geometrical objects, and then uses an interpolator to approximate the value of a function at these integration points.
For example, a quadrilateral object can use an integrator that uses 4 integration points.
* ------------- *
| |
| @ @ |
| |
| |
| @ @ |
| |
* ------------- *
where @ represents the location of the integration points. The interpolator approximates the value of a function at these integration points by using the values at the corner nodes, represented by *. You can think of it as each value of @ being a sort of average of the values of all *'s.
The following chart shows here, for convenience, the connection between the different classes used in this question:
Interpolator
|
| is a base for:
v
Interpolator_Template
|
| is a base for:
|
------------------------------------------
| | |
| | More shapes
V V
Interpolator_Quad Interpolator_Tria
| | |
| ... More QUADs More TRIAs
v
Interpolator_Quad_04
Integrator
|
| has member variable
v
vector<Abscissa*> abscissa
|
| Abscissa is a base class for:
|
--------------------------------
| | |
| More shapes |
V V
Abscissa_Quad Abscissar_Tria
Each geometrical shape has a different coordinate system, so my integrator and abscissae look like this:
class Integrator {
public:
void
integrate();
private:
/**
* An interpolator class.
*/
Interpolator * _interpolator;
/**
* List of abscissae (for the quadrilateral shown above, _abscissae.size() == 4).
*/
std::vector< Abscissa * > _abscissae;
};
Base class for all natural coordinates abscissae.
class Abscissa {
};
A quadrilateral abscissa operates on ξ and η natural coordinates.
class Abscissa_Quad final : public Abscissa {
public:
const double xi;
const double eta;
};
A triangle abscissa operates on ζ1, ζ2, and ζ3 natural coordinates.
class Abscissa_Tria final : public Abscissa {
public:
const double zeta_1;
const double zeta_2;
const double zeta_3;
};
The integrator implementation would integrate somewhat like this:
void
Integrator::integrate()
{
for ( Abscissa * abscissa : _abscissae ) {
_intepolator->eval_det_J( abscissa );
}
}
So far, so good. Let me show you my interpolator class.
class Interpolator {
public:
/**
* Evaluate the determinant of the Jacobian (important for numerical integration).
*
* @note Common interface for all abscissa types.
*
* @param[in] abscissa Integration abscissa.
* @return Shape functions.
*/
virtual double
eval_det_J(
const Abscissa * abscissa ) const = 0;
};
From the interpolator, I derive classes for all geometrical shapes. You may notice that I use an Interpolator_Template class as the base. Ignore it for now, I will explain the details in a second. This class contains functions common to all quadrilaterals.
class Interpolator_Quad : public Interpolator_Template< Abscissa_Quad > {
public:
// ... More functions common to all quadrilaterals.
};
This derived class corresponds to the quadrilateral drawn at the beginning of this question. The reason it is derived is because there might be quadrilaterals with more interpolation nodes. This class implements a QUAD_04 element (a quadrilateral with 4 interpolation nodes), but in finite elements we also have QUAD_08, QUAD_09, etc.
class Interpolator_Quad_04 final : public Interpolator_Quad {
public:
double
eval_det_J(
const Abscissa_Quad * abscissa ) const;
};
double
Interpolator_Quad_04::eval_det_J(
const Abscissa_Quad * abscissa ) const
{
// Important! Get coordinates from an Abscissa_Quad object.
const double xi = abscissa.xi;
const double eta = abscissa.eta;
double det_J = ...
// ... Perform some computations and return the determinant of the Jacobian.
return det_J;
}
Let me circle back to the Interpolator_Template class that I missed to explain before. At some point in my code, I perform a downcast from an Abscissa * to an Abscissa_Quad * object. I have achieved this by using a template class in combination with a non-virtual interface pattern.
template< class Abscissa_Derived >
class Interpolator_Template : public Interpolator {
public:
/**
* Implements Interpolator::eval_det_J.
*/
double
eval_det_J(
const Abscissa * abscissa ) const;
protected:
/**
* Implemented by Interpolator_Quad_04 in this example.
*/
virtual double
eval_det_J(
const Abscissa_Derived * abscissa ) const = 0;
private:
Abscissa_Derived *
eval_abscissa(
const Abscissa * abscissa ) const;
};
template< class Abscissa_Derived >
double
Interpolator_Template< Abscissa_Derived >::eval_det_J(
const Abscissa * abscissa ) const
{
Abscissa_Derived * abscissa_type = this->eval_abscissa( abscissa );
double det_J = this->eval_det_J( abscissa_type );
return det_J;
}
template< class Abscissa_Derived >
Abscissa_Derived *
Interpolator_Template< Abscissa_Derived >::eval_abscissa(
const Abscissa * abscissa ) const
{
// Dynamic cast occurs here.
// I will include some check later to check for nullptr.
return dynamic_cast< Abscissa_Derived * >( abscissa )
}
I'm sure this code contains errors because I had to copy and paste only what I thought was necessary to get my point across, as well as perform modifications. I hope, however, that my idea gets through correctly.
I know downcasting is usually a code smell, so before I went along and started implementing integrators and interpolators for all geometrical shapes in finite elements I wanted to get your opinion.
This is the last design pattern I have implemented. I will explain other designs I have attempted below; however, you can skip reading this section.
A double dispatch design pattern (visitor pattern specifically) where the integrator is the one being derived rather than the interpolator. For example, I had an Integrator_Quad_04 rather than an Interpolator_Quad_04. The Integrator_Quad_04 had a Abscissa_Quad as a member variable, as the abscissae were not longer derived.
class Integrator_Quad_04 final : public Integrator {
private:
std::vector< Abscissa_Quad * > _abscissae;
public:
double
eval_det_J(
const std::size_t & index,
const Interpolator * interpolator ) const
{
// The interpolator acts as the visitor.
interpolator->eval_det_J( _abscissa[index] );
}
}
/// Abscissa_Quad is no longer derived from Abscissa.
class Abscissa_Quad {
public:
const double xi;
const double eta;
};
The interpolator then becomes a visitor to the integrator class, and accesses its _abscissae member variable. I decided to not go along this design because then the interpolator would have had to be implemented based on the operations, rather than the shape.
class Interpolator {
// ...
};
class Eval_Det_J : public Interpolator {
double
eval_det_J(
const Abscissa_Quad * abscissa ) const;
double
eval_det_J(
const Abscissa_Tria * abscissa ) const;
};
I tried doing something with multiple dispatch, but the number of functions necessary for all shapes grew pretty fast.
Multiple variations of double dispatch + templating.
I found the current design pattern I'm using here:
Object Oriented Design Problem, Liskov Substitution Principle
As you may have inferred from the code, I'm compiling using C++11.
You may wonder why I don't simply combine the integrator and the interpolator into a single class, and the answer is because the integrator may operate on a sub-domain of the quadrilateral. For example, I could introduce a fictitious triangle inside the quadrilateral and place integration points inside the triangle, but I would still use the quadrilateral interpolation to approximate the solution inside the triangle points. Of course, I would need to implement some mapping between the triangle and the quadrilateral coordinates, but that's a problem for another day.
I want to know if you think downcasting is not a bad solution to this problem, or if I'm missing something. Perhaps I don't have the knowledge of a design pattern that fixes this issue, or perhaps my polymorphism architecture is incorrect.
Any feedback is appreciated. Thanks!
I'm posting this as an answer because it's too much for a comment, and maybe it helps.
You could make Integrator a template class:
template<class Shape>
class Integrator {
typedef typename Shape::AbscissaType AbscissaType;
public:
void integrate() const {
for (const AbscissaType& abscissa : _abscissae) {
_intepolator->eval_det_J(abscissa);
}
}
template<class OtherShape>
Integrator<OtherShape> convert(/* maybe pass range of abscissae indices */) const {
// Abscissae classes must have converting constructors like AbscissaTria(const AbscissaQuad&);
std::vector<typename OtherShape::AbscissaType> newAbscissae(_abscissae.begin(), _abscissae.end());
// initialize resulting integrator with newAbscissae
}
Interpolator_Template<AbscissaType> _interpolator;
std::vector<AbscissaType> _abscissae;
};
Specific integrators are inherited from appropriate base template, if it's still needed:
class Integrator_Quad_04 : public Integrator<Quad> {
};
So you don't need base Interpolator class and eval_det_J is a non-virtual function accepting appropriate type of abscissae:
template<class AbscissaType>
class Interpolator_Template {
public:
double eval_det_J(const AbscissaType& abscissa) const;
}
I added Shape here to emphasize that abscissa type depends on the shape, and you may have other things that depend on the shape too:
struct Tria {
typedef AbscissaTria AbscissaType;
static const size_t NUM_EDGES = 3; // just example
};
struct Quad {
typedef AbscissaQuad AbscissaType;
static const size_t NUM_EDGES = 4; // just example
};
I think you already have such classes in your code.
Also note that you also don't need base Abscissa anymore, all classes for abscissae become independent.
EDIT: if you need to convert abscissae to different type you can implement the following function in the Integrator class:
template<class OtherShape>
Integrator<OtherShape> convert(/* maybe pass range of abscissae indices */) const {
}
EDIT 2: Added sample implementation of convert into Integrator class.