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 <typename T>
class BaseExpr {};

template <std::size_t dim, std::size_t rank, typename... TRest, 
          template <std::size_t, std::size_t, typename...> typename T>
class BaseExpr<T<dim, rank, TRest...>>
{
// ...

public:
  template <typename... URest, template <std::size_t, std::size_t, typename...> typename U>
  SystemExpr<dim, rank, T<dim, rank, TRest...>, U<dim, rank, URest...>>
  operator ==(BaseExpr<U<dim, rank, URest...>> const &rhs) const
    { return {*this, rhs}; }

  SystemExpr<dim, rank, T<dim, rank, TRest...>, Tensor<dim, rank>>
  operator ==(Tensor<dim, rank> const &rhs) const
    { return {*this, rhs}; }
};

where SystemExpr<dim, rank, T<dim, rank, TRest...>, U<dim, rank, URest...>> 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<d, r1, T<d, r1, TRest ...>, VF::TTensor<d, r> > VF::TExprBase<T<d, r1, TRest ...> >::operator==(const VF::TTensor<d, r>&) const [with long unsigned int d = 3; long unsigned int r1 = 0; TRest = {VF::TExprBinaria<3, 0, VF::TExprd<3, 0>, double, std::divides<void> >, VF::TExprBinaria<3, 0, double, VF::TExprLap<3, 0, void>, std::multiplies<void> >, std::minus<void>}; 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<void> >, VF::TExprBinaria<3, 0, double, VF::TExprLap<3, 0, void>, std::multiplies<void> >, std::minus<void> >’

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

---

Here's a "minimal" example:

#include <cstddef>
#include <utility>
#include <functional>

template<typename>
class TExprBase;

template<std::size_t, std::size_t, typename, typename>
class TExprSistema;

template<std::size_t, std::size_t, typename, typename>
class TExprUnaria;

template<std::size_t, std::size_t, typename, typename, typename>
class TExprBinaria;

template<std::size_t, std::size_t>
class TCampo;

template<typename T>
class TExprBase {};

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

  template<typename... URest, template<std::size_t, std::size_t, typename...> typename U>
  TExprBinaria<d, r1, T<d, r1, TRest...>, U<d, r1, URest...>, std::minus<>>
  operator -(TExprBase<U<d, r1, URest...>> const &rhs) const
    { return {*this, rhs}; }

  TExprUnaria<d, r1, T<d, r1, TRest...>, std::negate<>>
  operator -() const
    { return {*this}; }

  template<std::size_t r2, typename... URest,
           template<std::size_t, std::size_t, typename...> typename U>
  TExprBinaria<d, r1 + r2, T<d, r1, TRest...>, U<d, r2, URest...>, std::multiplies<>>
  operator *(TExprBase<U<d, r2, URest...>> const &rhs) const
    { return {*this, rhs}; }

  TExprBinaria<d, r1, T<d, r1, TRest...>, double, std::multiplies<>>
  operator *(double const rhs) const
    { return {*this, rhs}; }

  template<typename... URest, template<std::size_t, std::size_t, typename...> typename U>
  TExprBinaria<d, r1, T<d, r1, TRest...>, U<d, 0u, URest...>, std::divides<>>
  operator /(TExprBase<U<d, 0u, URest...>> const &rhs) const
    { return {*this, rhs}; }

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

  template<typename... URest, template<std::size_t, std::size_t, typename...> typename U>
  TExprSistema<d, r1, T<d, r1, TRest...>, U<d, r1, URest...>>
  operator ==(TExprBase<U<d, r1, URest...>> const &rhs) const
    { return {*this, rhs}; }

  operator T<d, r1, TRest...> const &() const
    { return *static_cast<T<d, r1, TRest...> const *>(this); }
};

template<std::size_t d, std::size_t r, typename T, typename U>
class TExprSistema : public TExprBase<TExprSistema<d, r, T, U>>
{
private:
  TExprBase<T> const &lhs;
  TExprBase<U> const &rhs;

public:
  TExprSistema() = delete;

  TExprSistema(TExprBase<T> const &lhs_, TExprBase<U> const &rhs_) :
    lhs(lhs_), rhs(rhs_) {}

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

template<std::size_t d, std::size_t r, typename T, typename TOp>
class TExprUnaria : public TExprBase<TExprUnaria<d, r, T, TOp>>
{
private:
  T const &rhs;
  [[no_unique_address]] TOp const Op = {};

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

template<std::size_t d, std::size_t r, typename T, typename U, typename TOp>
class TExprBinaria : public TExprBase<TExprBinaria<d, r, T, U, TOp>>
{
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<std::size_t d, std::size_t r>
class TCampo : public TExprBase<TCampo<d, r>> {};

template<std::size_t d, std::size_t r, typename T>
class TExprDiv : public TExprBase<TExprDiv<d, r, T>> {};

template<std::size_t d, std::size_t r>
class TExprGrad : public TExprBase<TExprGrad<d, r>> {};

template<std::size_t d, std::size_t r, typename T>
class TExprLap : public TExprBase<TExprLap<d, r, T>> {};

template<std::size_t d, std::size_t r, typename... TRest,
         template<std::size_t, std::size_t, typename...> typename T>
TExprDiv<d, r - 1u, T<d, 1u, TRest...>>
inline div(TExprBinaria<d, r, T<d, 1u, TRest...>, TCampo<d, r - 1u>, std::multiplies<>> const &)
  { return {}; }

template<std::size_t d, std::size_t r>
TExprGrad<d, r + 1u>
inline grad(TCampo<d, r> const &)
  { return {}; }

template<std::size_t d, std::size_t r>
TExprLap<d, r, void>
inline lap(TCampo<d, r> const &)
  { return {}; }

template<std::size_t d, std::size_t r>
class TSistema
{
public:
  template<typename T, typename U>
  TSistema(TExprSistema<d, r, T, U> const &);

  void
  Solve() const;

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

template<std::size_t d, std::size_t r, typename T, typename U>
void
inline solve(TExprSistema<d, r, T, U> 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;
}

Answer 1

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

template <typename T>
struct Other { };

template <typename T>
struct Base {
    template <typename U>
    void operator==(Base<U> const&) const;

    template <typename U>
    void operator==(Other<U> const&) const;
};

struct A : Base<A> { };
template <typename T> struct B : Base<B<T>> { };

Basically, what you have is some CRTP base class template which is comparable to any other such thing, and also with any Other<U>. 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<T>, and Other<U> work?

---

The answer is: it probably does what you want.

Comparing A{} to B<int>{} or B<double>{} to B<char>{} 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<K> to an A in either direction. This is a symmetric comparison.

---

The other types are more interesting.

Comparing A{} to Other<T>{} 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<T>{} 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<T>{} 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<int>' and 'A')
    b == a;
    ~ ^  ~