I'm writing a CFD application using Eigen for a lot of the calculations. Now I want to define a field class to hold the values for each variable. My current approach is to use a template class which is either instantiated for scalars and vectors or inherited by a ScalarField and a VectorField class and then mildly specialised. Scalars are stored in a 1 x n matrix, vectors in a 2 x n matrix (in 2D).
Now the problems arise when I want to access individual scalars or (column) vectors of those fields, as the return type of Eigen's access functions are expression templates, not the data type with which I instantiated the template. But I need to access them from inside the base class, otherwise I'd have a lot of duplicate code.
So in the end, I found two approaches that worked, which both involve overloading the subscript operator []:
pure virtual functions in the base class template, with the return type being decided by a std::conditional, based on the data type with which the template is being instantiated, and then having rather readable overrides in the derived classes and
using SFINAE to define all functions inside the base class template, using std::enable_if, thus avoiding inheritance altogether. But those template functions aren't very readable.
As the return type of operator[] for vectors I used Eigen::Block<...> by the way.
Here's a minimal working example, with the functions get() and get2() standing in for the virtual function approach and the SFINAE approach respectively, as well as int and double for scalars and vectors, so that it can be easily compiled.
#include <iostream>
#include <type_traits>
template<typename T>
class Base
{
protected:
T m_data;
public:
Base(T t) :
m_data {t}
{}
virtual typename std::conditional<
std::is_same<T, double>::value, int, T>::type
get() = 0;
template<typename U = T>
std::enable_if_t<std::is_integral<U>::value,T> get2(){
return m_data;
}
template<typename U = T>
std::enable_if_t<std::is_floating_point<U>::value,int> get2(){
return m_data;
}
void print() {
std::cout << this->get() << "\n";
std::cout << this->get2() << "\n";
}
};
class Int : public Base<int>
{
public:
Int(int i) :
Base(i)
{}
virtual int get() override
{ return m_data; }
};
class Double : public Base<double>
{
public:
Double(double d) :
Base<double>(d)
{}
virtual int get() override
{ return m_data; }
};
int main()
{
Int i { 1 };
i.print();
Double d { 1.1 };
d.print();
//Base<int> b1 { 1 };
//b1.print();
//Base<double> b2 { 1.1 };
//b2.print();
}
Now my question is: is there a more elegant and maintainable way to achieve my goal of template specialisation (and in the current case, accessing individual elements)? And if not, is one of the above ways preferable or can it be improved? Any feedback welcome, I'm rather new to C++.
EDIT: As requested, here's something more similar to my actual code.
I'm using the subscript operator in setBoundaryPatch() and setInternalField and their overloads to access m_data. Because for scalars, m_data[] works fine, but for vectors or blocks, I need to use m_data.col() or m_data.block().
I attempted just using m_data.col() or m_data.block() for both scalars and vectors, but there's no conversion from a scalar to the respective expression. I guess I could use m_data.block(...) = T { value }, but wouldn't that construct a temporary object for each call and be rather slow?
Also, overloading the subscript operator just makes it rather convenient to work with the class.
#include <Eigen/Dense>
#include <type_traits>
#include "mesh.h"
namespace fvm
{
constexpr auto dims { 2 }; // dimensions
using vector = Eigen::Matrix<double, dims, 1>;
using field = Eigen::Matrix<double, n, Eigen::Dynamic>;
template<typename T>
class Field
{
protected:
static constexpr int n {
(std::is_same<T, double>::value ? 1 : fvm::dims ) };
/* because of some more member variables than just the data, e.g. a
* reference to the mesh, I don't want to use an Eigen type directly. */
const fvm::Mesh& m_mesh;
fvm::field<n> m_data;
public:
Field( const fvm::Mesh& mesh ) :
m_mesh { mesh }
{
m_data.resize( Eigen::NoChange, m_mesh.cellCount() );
}
const fvm::field<n>& data() const
{ return m_data; }
fvm::field<n>& data()
{ return m_data; }
/* sets the field values belonging to the ghost cells of a boundary patch to
* a uniform value */
void setBoundaryPatch( const std::string& patchName, T )
{
for ( auto faceInd : m_mesh.boundaryPatch(patchName) )
(*this)[ m_mesh.face(faceInd).ghostCell().ID() ] = t;
}
/* and a couple overloads for non-uniform values */
/* sets the field values belonging to domain cells to a uniform value */
void setInternalField( T );
/* and a couple overloads for non-uniform values */
protected:
using col = Eigen::Block< fvm::field<n> >;
using constCol = const Eigen::Block< const fvm::field<n> >;
public:
/* using SFINAE to define subscript operator[] */
//template<typename U = T>
//std::enable_if_t<!std::is_same<U, double>::value, constCol>
//operator[] (int i) const {
// return m_data.block(0, i, n, 1);
//}
//template<typename U = T>
//std::enable_if_t<!std::is_same<U, double>::value, col>
//operator[] (int i) {
// return m_data.block(0, i, n, 1);
//}
//template<typename U = T>
//std::enable_if_t<std::is_same<U, double>::value, T>
//operator[] (int i) const {
// return m_data[i];
//}
//template<typename U = T>
//std::enable_if_t<std::is_same<U, double>::value, T&>
//operator[] (int i) {
// return m_data[i];
//}
/* using pure virtual functions to overload the subscript operator[] */
virtual typename std::conditional<
std::is_same<T, fvm::vector>::value, constCol, T>::type
operator[] (int) const = 0;
virtual typename std::conditional<
std::is_same<T, fvm::vector>::value, col, T&>::type
operator[] (int) = 0;
virtual void readFromFile( const std::string& path ) = 0;
virtual void writeToFile( const std::string& path ) const = 0;
};
/* if I defined everything in the template, I could just declare aliases
* -> SFINAE option*/
//using ScalarField = Field<fvm::scalar>;
//using VectorField = Field<fvm::vector>;
/* or I define them as classes and let them inherit from Field<> */
class ScalarField : public Field<fvm::scalar>
{
public:
ScalarField( const fvm::Mesh& mesh );
virtual fvm::scalar operator[] (int) const override;
virtual fvm::scalar& operator[] (int) override;
virtual void writeToFile( const std::string& path ) const override;
virtual void readFromFile( const std::string& path ) override;
};
class VectorField : public Field<fvm::vector>
{
public:
VectorField( const fvm::Mesh& );
virtual VectorField::constCol operator[] (int) const override;
virtual VectorField::col operator[] (int) override;
virtual void writeToFile( const std::string& path ) const override;
virtual void readFromFile( const std::string& path ) override;
};
} // end namespace fvm
I'm still very open to answers and feedback, however, I at least found a less verbose and more elegant solution to what I had before, by is using if constexpr (C++17 is pretty awesome!):
using col = Eigen::Block< fvm::field<n> >;
using constCol = const Eigen::Block< const fvm::field<n> >;
using elem = typename std::conditional<n==1, fvm::scalar&, col >::type;
using constElem = typename std::conditional<n==1, fvm::scalar , constCol>::type;
elem operator[] (int) {
if constexpr ( n==1 )
return m_data[i];
else
return m_data.block(0, i, n, 1);
}
constElem operator[] (int) const {
if constexpr ( n==1 )
return m_data[i];
else
return m_data.block(0, i, n, 1);
}
EDIT: What makes this solution even better, is that now, automatic type deduction actually works, reducing the code to
auto operator[] (int) {
if constexpr ( n==1 )
return m_data[i];
else
return m_data.block(0, i, n, 1);
}
auto operator[] (int) const {
if constexpr ( n==1 )
return m_data[i];
else
return m_data.block(0, i, n, 1);
}
EDIT2: I was wrong about auto being usable. The .obj-files of my unit test executable simply still had the correct type stored, but they couldn't be recreated using auto.