Dear Eigen Users

I would like to implement a constitutive matrix occurring in plane stress conditions. I want this matrix to behave basically as a Eigen::Vector3d, except for the * operator, which I have overloaded in the following fashion:

Code: Select all

using TVector3=    Eigen::Matrix<double,3,1>;

class TMatTensPlana: public TVector3
{
public:
                TMatTensPlana()= default;
                TMatTensPlana(double C11,double C12,double C33): TVector3(C11,C12,C33) {}
  template <typename Derived>
                TMatTensPlana(const MatrixBase<Derived> &Other): TVector3(Other) {}
  REAL_         C11() const { return x(); }
  REAL_         C12() const { return y(); }
  REAL_         C21() const { return y(); }
  REAL_         C22() const { return x(); }
  REAL_         C33() const { return z(); }
  TVector3      operator *(const TVector3 &v) const
                  { return {C11()*v.x()+C12()*v.y(),C21()*v.x()+C22()*v.y(),C33()*v.z()}; }
};

I realize this implementation breaks Eigen’s lazy evaluation and vectorization, yet it works as I wanted. It is worth to implement a operator:

Code: Select all

typename ProductReturnType<Derived, OtherDerived>::Type TMatTensPlana::operator*(const MatrixBase<OtherDerived> &other) const

or whathever? Could somebody please provide advice on how to implement it?

Thanks so much

For small objects, it is perfectly fine to return by values. Since the values will stay in register, the compiler will easily optimize all temporaries away. An index-based lazy evaluator will even be slower because of the complex management of the shuffled indices.

BTW, be careful with such an implementation, because the following won't work as expected:

TVector3 a, b;

(a+a)*(b+b)

will call MatrixBase::operator*, not TVector3::operator*. In this case you'll get a compiler error because this matrix product is invalid, and since there is not matrix product taking 2 vectors and returning one vector, you should always get a compilation error in case of miss-use, but still, be aware of this possible loophole.

ggael wrote:will call MatrixBase::operator*, not TVector3::operator*.

Yes. This confused me.

I thought:

Code: Select all

TMatTensPlana C;
TVector3 epsT, epsV, sigma;

sigma= C*(epsT+epsV);


was going to fail because the result of epsT+epsV would have not been implicitly promoted to TVect3D, and thus calling the "incorrect" operator * and, consequently, get an assertion failed of incompatible matrix sizes.

However, it seem to work fine, at least with VC2015 Update 2. I'll try later with gcc and clang and post the result if relevant.

BTW, thanks for your prompt response.