Rewritten comparison operators and expression templates

Question

I have a finite volume library, strongly influenced by openfoam, which enables the solution to continuum mechanics problems to be written in C++ similarly as one would in paper. For example, to solve the Navier-Stokes equation for incompressible, laminar flow I simply write:

solve(div(U * U) - nu * lap(U) == -grad(p));

This expression, of course, involves expression templates so that the system coefficients are calculated single-pass, and the system is solved in a matrix-free fashion.

Those expression templates are based on the CRTP, so that all of the necessary operators are defined in a base class. In particular, the equality operators are defined like this:

template 
class BaseExpr {};

template  typename T>
class BaseExpr>
{
// ...

public:
  template  typename U>
  SystemExpr, U>
  operator ==(BaseExpr> const &rhs) const
    { return {*this, rhs}; }

  SystemExpr, Tensor>
  operator ==(Tensor const &rhs) const
    { return {*this, rhs}; }
};

where SystemExpr, U> has functions to lazily evaluating the coefficients and calculating the solution.

In contrast with OpenFOAM, where you have to qualify, for example, fvm::lap, fvc::grad, etc. to signal that that term is to be evaluated explicitly or implicitly, I adopted the convention that I have always used in paper: everything in the left hand side is evaluated implicitly, while the right hand side is evaluated explicitly. Therefore, BaseExpr::operator == is not commutative. That becomes increasingly useful as the equations get longer. For example, the ε transport equation for compressible flow is:

solve(div(φ * ε) - div(µt * grad(ε)) / σε + ρ * Sp((C2 * ω + 2.0 / 3.0 * C1 * div(U)) * ε) == C1 * ω * G);

I'm concerned that under C++20 this design may broken because of the new "rewritten candidates" of operator ==. G++-10 compiles the library without any warning with -std=c++20 -Wall -Wextra -pedantic but I would like to double-check: is the above code well-formed under C++20?

I know that some people may think that the above is a terrible design, but I love it to the point that I would prefer to stay in -std=c++17 mode instead of using a different operator (say, operator >> or whatever) to represent equality (yes, that's a form of equality).

With GCC-10.2 I'm getting the behaviour that I want. Consider the unsteady heat transfer equation:

solve(d(T) / dt - α * lap(T) == 0); // OK: implicit scheme
  
solve(d(T) / dt == α * lap(T)); // OK: explicit scheme

solve(d(T) / dt - 0.5 * α * lap(T) == 0.5 * α * lap(T)); // OK: Crank-Nicholson scheme
  
solve(0 == d(T) / dt - α * lap(T)); // OK: doesn't compile

It is perfectly fine that the last example does not compile, because it doen't make sense. The error I get is:

prova.cpp:51:11: error: return type of ‘VF::TExprSistema, VF::TTensor > VF::TExprBase >::operator==(const VF::TTensor&) const [with long unsigned int d = 3; long unsigned int r1 = 0; TRest = {VF::TExprBinaria<3, 0, VF::TExprd<3, 0>, double, std::divides >, VF::TExprBinaria<3, 0, double, VF::TExprLap<3, 0, void>, std::multiplies >, std::minus}; T = VF::TExprBinaria]’ is not ‘bool’
   51 |   solve(0 == d(T) / dt - α * lap(T));
      |         ~~^~~~~~~~~~~~~~~~~~~~~~~~~
prova.cpp:51:11: note: used as rewritten candidate for comparison of ‘int’ and ‘VF::TExprBinaria<3, 0, VF::TExprBinaria<3, 0, VF::TExprd<3, 0>, double, std::divides >, VF::TExprBinaria<3, 0, double, VF::TExprLap<3, 0, void>, std::multiplies >, std::minus >’

Therefore, GCC appears to behave exactly as I want even in C++20 mode.

Here's a "minimal" example:

#include 
#include 
#include 

template
class TExprBase;

template
class TExprSistema;

template
class TExprUnaria;

template
class TExprBinaria;

template
class TCampo;

template
class TExprBase {};

template typename T>
class TExprBase>
{
public:
  template typename U>
  TExprBinaria, U, std::plus<>>
  operator +(TExprBase> const &rhs) const
    { return {*this, rhs}; }

  template typename U>
  TExprBinaria, U, std::minus<>>
  operator -(TExprBase> const &rhs) const
    { return {*this, rhs}; }

  TExprUnaria, std::negate<>>
  operator -() const
    { return {*this}; }

  template typename U>
  TExprBinaria, U, std::multiplies<>>
  operator *(TExprBase> const &rhs) const
    { return {*this, rhs}; }

  TExprBinaria, double, std::multiplies<>>
  operator *(double const rhs) const
    { return {*this, rhs}; }

  template typename U>
  TExprBinaria, U, std::divides<>>
  operator /(TExprBase> const &rhs) const
    { return {*this, rhs}; }

  TExprBinaria, double, std::divides<>>
  operator /(double const rhs) const
    { return {*this, rhs}; }

  template typename U>
  TExprSistema, U>
  operator ==(TExprBase> const &rhs) const
    { return {*this, rhs}; }

  operator T const &() const
    { return *static_cast const *>(this); }
};

template
class TExprSistema : public TExprBase>
{
private:
  TExprBase const &lhs;
  TExprBase const &rhs;

public:
  TExprSistema() = delete;

  TExprSistema(TExprBase const &lhs_, TExprBase const &rhs_) :
    lhs(lhs_), rhs(rhs_) {}

  TExprSistema(TExprBase const &lhs_, U const &rhs_) :
    lhs(lhs_), rhs(rhs_) {}
};

template
class TExprUnaria : public TExprBase>
{
private:
  T const &rhs;
  [[no_unique_address]] TOp const Op = {};

public:
  TExprUnaria(T const &rhs_) :
    rhs(rhs_) {}
};

template
class TExprBinaria : public TExprBase>
{
private:
  T const &lhs;
  U const &rhs;
  [[no_unique_address]] TOp const Op = {};

public:
  TExprBinaria(T const &lhs_, U const &rhs_) :
    lhs(lhs_), rhs(rhs_) {}
};

template
class TCampo : public TExprBase> {};

template
class TExprDiv : public TExprBase> {};

template
class TExprGrad : public TExprBase> {};

template
class TExprLap : public TExprBase> {};

template typename T>
TExprDiv>
inline div(TExprBinaria, TCampo, std::multiplies<>> const &)
  { return {}; }

template
TExprGrad
inline grad(TCampo const &)
  { return {}; }

template
TExprLap
inline lap(TCampo const &)
  { return {}; }

template
class TSistema
{
public:
  template
  TSistema(TExprSistema const &);

  void
  Solve() const;

  void
  friend solve(TSistema const &Sistema)
    { Sistema.Solve(); }
};

template
void
inline solve(TExprSistema const &Expr)
  { solve(TSistema(Expr)); }

int main()
{
  TCampo<3u, 1u> U;
  TCampo<3u, 0u> nu, p;

  solve(div(U * U) - nu * lap(U) == -grad(p));

  return 0;
}

user_2069064 · Answer

I'm going to severely reduce your example as follows:

template 
struct Other { };

template 
struct Base {
    template 
    void operator==(Base const&) const;

    template 
    void operator==(Other const&) const;
};

struct A : Base { };
template  struct B : Base> { };

Basically, what you have is some CRTP base class template which is comparable to any other such thing, and also with any Other. I picked the simplest non-bool return type for these, which is to say... void.

So the question becomes: does the above work? Will every combination of comparing A, B, and Other work?

The answer is: it probably does what you want.

Comparing A{} to B{} or B{} to B{} is fine. In those kinds of scenarios, we have both the real candidate (from the left-hand side) and the rewritten candidate (from the right-hand side) and both candidates involve a derived-to-base conversion on both arguments, so the non-rewritten candidate is preferred.

This is true even if we directly compare a Base to an A in either direction. This is a symmetric comparison.

The other types are more interesting.

Comparing A{} to Other{} in that direction is fine. We invoke the member operator== as we would in C++17, that is the only candidate, and it does what you'd expect.

Comparing Other{} to A{} (i.e., in the other direction) was ill-formed in C++17 due to not having any candidates. But it's ill-formed in C++20 for a different reason. Now we do have a candidate: the reversed A{} == Other{} candidate. But we end up rejecting the candidate later due to [over.match.oper]/9:

If a rewritten operator== candidate is selected by overload resolution for an operator @, its return type shall be cv bool, and [...]

Note that the return type of bool requirement comes into play at the end of overload resolution (if a non-bool-returning rewritten candidate wins, it's ill-formed) rather than than at the beginning (as you alluded to in your comment, suggesting that non-bool-returning functions weren't even considered as rewrite candidates).

So we do consider this candidate (which is our only candidate), but we reject due to it being an invalid rewritten candidate. As a result, it's still ill-formed in C++20 as it was in C++17. It's just now ill-formed for a different reason, so you'll get a different error message. For instance, Clang gives me:

error: return type 'void' of selected 'operator==' function for rewritten '==' comparison is not 'bool'
    b == a;
    ~ ^  ~

whereas in C++17 it would've given:

error: invalid operands to binary expression ('Other' and 'A')
    b == a;
    ~ ^  ~

Rewritten comparison operators and expression templates

1 Answers