test_DiscreteFunctionP0Vector.cpp

#include <catch2/catch_test_macros.hpp>
#include <catch2/matchers/catch_matchers_all.hpp>

#include <MeshDataBaseForTests.hpp>
#include <scheme/DiscreteFunctionP0Vector.hpp>

// clazy:excludeall=non-pod-global-static

TEST_CASE("DiscreteFunctionP0Vector", "[scheme]")
{
  auto same_values = [](const auto& f, const auto& g) {
    const size_t number_of_cells = f.cellArrays().numberOfItems();
    const size_t size_of_arrays  = f.cellArrays().sizeOfArrays();
    for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
      for (size_t i = 0; i < size_of_arrays; ++i) {
        if (f[cell_id][i] != g[cell_id][i]) {
          return false;
        }
      }
    }
    return true;
  };

  SECTION("constructors")
  {
    SECTION("1D")
    {
      const size_t size = 3;

      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh1D();
      constexpr size_t Dimension = 1;

      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0Vector);
      REQUIRE(f.size() == size);

      REQUIRE(f.mesh().get() == mesh.get());

      DiscreteFunctionP0Vector g{f};
      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0Vector);
      REQUIRE(g.size() == size);

      CellArray<double> h_arrays{mesh->connectivity(), size};
      h_arrays.fill(0);

      DiscreteFunctionP0Vector zero{mesh, [&] {
                                      CellArray<double> cell_array{mesh->connectivity(), size};
                                      cell_array.fill(0);
                                      return cell_array;
                                    }()};

      DiscreteFunctionP0Vector h{mesh, h_arrays};
      REQUIRE(same_values(h, zero));
      REQUIRE(same_values(h, h_arrays));

      h_arrays.fill(1);

      REQUIRE(same_values(h, h_arrays));
      REQUIRE(not same_values(h, zero));

      DiscreteFunctionP0Vector moved_h{std::move(h)};
      REQUIRE(same_values(moved_h, h_arrays));
    }

    SECTION("2D")
    {
      const size_t size = 3;

      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh2D();
      constexpr size_t Dimension = 2;

      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0Vector);
      REQUIRE(f.size() == size);

      REQUIRE(f.mesh().get() == mesh.get());

      DiscreteFunctionP0Vector g{f};
      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0Vector);
      REQUIRE(g.size() == size);

      CellArray<double> h_arrays{mesh->connectivity(), size};
      h_arrays.fill(0);

      DiscreteFunctionP0Vector zero{mesh, [&] {
                                      CellArray<double> cell_array{mesh->connectivity(), size};
                                      cell_array.fill(0);
                                      return cell_array;
                                    }()};

      DiscreteFunctionP0Vector h{mesh, h_arrays};
      REQUIRE(same_values(h, zero));
      REQUIRE(same_values(h, h_arrays));

      h_arrays.fill(1);

      REQUIRE(same_values(h, h_arrays));
      REQUIRE(not same_values(h, zero));

      DiscreteFunctionP0Vector moved_h{std::move(h)};
      REQUIRE(same_values(moved_h, h_arrays));
    }

    SECTION("3D")
    {
      const size_t size = 2;

      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh3D();
      constexpr size_t Dimension = 3;

      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
      REQUIRE(f.dataType() == ASTNodeDataType::double_t);
      REQUIRE(f.descriptor().type() == DiscreteFunctionType::P0Vector);
      REQUIRE(f.size() == size);

      REQUIRE(f.mesh().get() == mesh.get());

      DiscreteFunctionP0Vector g{f};
      REQUIRE(g.dataType() == ASTNodeDataType::double_t);
      REQUIRE(g.descriptor().type() == DiscreteFunctionType::P0Vector);
      REQUIRE(g.size() == size);

      CellArray<double> h_arrays{mesh->connectivity(), size};
      h_arrays.fill(0);

      DiscreteFunctionP0Vector zero{mesh, [&] {
                                      CellArray<double> cell_array{mesh->connectivity(), size};
                                      cell_array.fill(0);
                                      return cell_array;
                                    }()};

      DiscreteFunctionP0Vector h{mesh, h_arrays};
      REQUIRE(same_values(h, zero));
      REQUIRE(same_values(h, h_arrays));

      h_arrays.fill(1);

      REQUIRE(same_values(h, h_arrays));
      REQUIRE(not same_values(h, zero));

      DiscreteFunctionP0Vector moved_h{std::move(h)};
      REQUIRE(same_values(moved_h, h_arrays));
    }
  }

  SECTION("fill")
  {
    auto all_values_equal = [](const auto& f, const auto& g) {
      const size_t number_of_cells = f.cellArrays().numberOfItems();
      size_t size_of_arrays        = f.cellArrays().sizeOfArrays();
      for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
        for (size_t i = 0; i < size_of_arrays; ++i) {
          if (f[cell_id][i] != g) {
            return false;
          }
        }
      }
      return true;
    };

    SECTION("1D")
    {
      const size_t size = 3;

      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh1D();
      constexpr size_t Dimension = 1;

      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
      f.fill(3);

      REQUIRE(all_values_equal(f, 3));
    }

    SECTION("2D")
    {
      const size_t size = 3;

      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh2D();
      constexpr size_t Dimension = 2;

      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
      f.fill(2.3);

      REQUIRE(all_values_equal(f, 2.3));
    }

    SECTION("3D")
    {
      const size_t size = 2;

      std::shared_ptr mesh       = MeshDataBaseForTests::get().cartesianMesh3D();
      constexpr size_t Dimension = 3;

      DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
      f.fill(3.2);

      REQUIRE(all_values_equal(f, 3.2));
    }
  }

  SECTION("copies")
  {
    auto all_values_equal = [](const auto& f, const auto& g) {
      const size_t number_of_cells = f.cellArrays().numberOfItems();
      const size_t size_of_arrays  = f.cellArrays().sizeOfArrays();

      for (CellId cell_id = 0; cell_id < number_of_cells; ++cell_id) {
        for (size_t i = 0; i < size_of_arrays; ++i) {
          if (f[cell_id][i] != g) {
            return false;
          }
        }
      }
      return true;
    };

    SECTION("1D")
    {
      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();

      constexpr size_t Dimension = 1;

      const size_t size  = 3;
      const size_t value = parallel::rank() + 1;
      const size_t zero  = 0;

      DiscreteFunctionP0Vector<Dimension, size_t> f{mesh, size};
      f.fill(value);

      REQUIRE(all_values_equal(f, value));

      DiscreteFunctionP0Vector g = copy(f);
      f.fill(zero);

      REQUIRE(all_values_equal(f, zero));
      REQUIRE(all_values_equal(g, value));

      copy_to(g, f);
      g.fill(zero);

      DiscreteFunctionP0Vector<Dimension, const size_t> h = copy(f);

      DiscreteFunctionP0Vector<Dimension, size_t> shallow_g{mesh, size};
      shallow_g = g;

      REQUIRE(all_values_equal(f, value));
      REQUIRE(all_values_equal(g, zero));
      REQUIRE(all_values_equal(shallow_g, zero));
      REQUIRE(all_values_equal(h, value));

      copy_to(h, g);

      REQUIRE(all_values_equal(g, value));
      REQUIRE(all_values_equal(shallow_g, value));
    }

    SECTION("2D")
    {
      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();

      constexpr size_t Dimension = 2;

      const size_t size  = 3;
      const size_t value = parallel::rank() + 1;
      const size_t zero  = 0;

      DiscreteFunctionP0Vector<Dimension, size_t> f{mesh, size};
      f.fill(value);

      REQUIRE(all_values_equal(f, value));

      DiscreteFunctionP0Vector g = copy(f);
      f.fill(zero);

      REQUIRE(all_values_equal(f, zero));
      REQUIRE(all_values_equal(g, value));

      copy_to(g, f);
      g.fill(zero);

      DiscreteFunctionP0Vector<Dimension, const size_t> h = copy(f);

      DiscreteFunctionP0Vector<Dimension, size_t> shallow_g{mesh, size};
      shallow_g = g;

      REQUIRE(all_values_equal(f, value));
      REQUIRE(all_values_equal(g, zero));
      REQUIRE(all_values_equal(shallow_g, zero));
      REQUIRE(all_values_equal(h, value));

      copy_to(h, g);

      REQUIRE(all_values_equal(g, value));
      REQUIRE(all_values_equal(shallow_g, value));
    }

    SECTION("3D")
    {
      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();

      constexpr size_t Dimension = 3;

      const size_t size  = 3;
      const size_t value = parallel::rank() + 1;
      const size_t zero  = 0;

      DiscreteFunctionP0Vector<Dimension, size_t> f{mesh, size};
      f.fill(value);

      REQUIRE(all_values_equal(f, value));

      DiscreteFunctionP0Vector g = copy(f);
      f.fill(zero);

      REQUIRE(all_values_equal(f, zero));
      REQUIRE(all_values_equal(g, value));

      copy_to(g, f);
      g.fill(zero);

      DiscreteFunctionP0Vector<Dimension, const size_t> h = copy(f);

      DiscreteFunctionP0Vector<Dimension, size_t> shallow_g{mesh, size};
      shallow_g = g;

      REQUIRE(all_values_equal(f, value));
      REQUIRE(all_values_equal(g, zero));
      REQUIRE(all_values_equal(h, value));

      copy_to(h, g);

      REQUIRE(all_values_equal(g, value));
      REQUIRE(all_values_equal(shallow_g, value));
    }
  }

  SECTION("unary operators")
  {
    SECTION("1D")
    {
      const size_t size    = 3;
      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();

      constexpr size_t Dimension = 1;

      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();

      SECTION("unary minus")
      {
        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            const double x = xj[cell_id][0];
            for (size_t i = 0; i < size; ++i) {
              f[cell_id][i] = 2 * x + i;
            }
          });

        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;

        Table<double> minus_values{mesh->numberOfCells(), size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            for (size_t i = 0; i < size; ++i) {
              minus_values[cell_id][i] = -f[cell_id][i];
            }
          });

        REQUIRE(same_values(-f, minus_values));
        REQUIRE(same_values(-const_f, minus_values));
      }
    }

    SECTION("2D")
    {
      const size_t size    = 3;
      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();

      constexpr size_t Dimension = 2;

      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();

      SECTION("unary minus")
      {
        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            const double x = xj[cell_id][0];
            const double y = xj[cell_id][1];
            for (size_t i = 0; i < size; ++i) {
              f[cell_id][i] = 2 * x + i * y;
            }
          });

        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;

        Table<double> minus_values{mesh->numberOfCells(), size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            for (size_t i = 0; i < size; ++i) {
              minus_values[cell_id][i] = -f[cell_id][i];
            }
          });

        REQUIRE(same_values(-f, minus_values));
        REQUIRE(same_values(-const_f, minus_values));
      }
    }

    SECTION("3D")
    {
      const size_t size    = 2;
      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();

      constexpr size_t Dimension = 3;

      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();

      SECTION("unary minus")
      {
        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            const double x = xj[cell_id][0];
            const double y = xj[cell_id][1];
            const double z = xj[cell_id][2];
            for (size_t i = 0; i < size; ++i) {
              f[cell_id][i] = 2 * x + i * y - z;
            }
          });

        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;

        Table<double> minus_values{mesh->numberOfCells(), size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            for (size_t i = 0; i < size; ++i) {
              minus_values[cell_id][i] = -f[cell_id][i];
            }
          });

        REQUIRE(same_values(-f, minus_values));
        REQUIRE(same_values(-const_f, minus_values));
      }
    }
  }

  SECTION("binary operators")
  {
    SECTION("1D")
    {
      const size_t size = 3;

      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh1D();

      constexpr size_t Dimension = 1;

      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();

      SECTION("inner operators")
      {
        SECTION("scalar functions")
        {
          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
          parallel_for(
            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
              const double x = xj[cell_id][0];
              f[cell_id][0]  = 2 * x + 1;
              f[cell_id][1]  = x * x - 1;
              f[cell_id][2]  = 2 + x;
            });

          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
          parallel_for(
            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
              const double x = xj[cell_id][0];
              g[cell_id][0]  = (x + 1) * (x - 2) + 1;
              g[cell_id][1]  = 3 * (x + 2) - 1;
              g[cell_id][2]  = (x + 3) * 5;
            });

          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};

          SECTION("sum")
          {
            Table<double> sum_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  sum_values[cell_id][i] = f[cell_id][i] + g[cell_id][i];
                }
              });

            REQUIRE(same_values(f + g, sum_values));
            REQUIRE(same_values(const_f + g, sum_values));
            REQUIRE(same_values(f + const_g, sum_values));
            REQUIRE(same_values(const_f + const_g, sum_values));
          }

          SECTION("difference")
          {
            Table<double> difference_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  difference_values[cell_id][i] = f[cell_id][i] - g[cell_id][i];
                }
              });

            REQUIRE(same_values(f - g, difference_values));
            REQUIRE(same_values(const_f - g, difference_values));
            REQUIRE(same_values(f - const_g, difference_values));
            REQUIRE(same_values(const_f - const_g, difference_values));
          }
        }
      }

      SECTION("external operators")
      {
        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            const double x = xj[cell_id][0];
            for (size_t i = 0; i < size; ++i) {
              f[cell_id][i] = std::abs(2 * x) + i;
            }
          });

        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;

        SECTION("product")
        {
          SECTION("scalar lhs")
          {
            const double a = 3.2;
            Table<double> product_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  product_values[cell_id][i] = a * f[cell_id][i];
                }
              });

            REQUIRE(same_values(a * f, product_values));
            REQUIRE(same_values(a * const_f, product_values));
          }

          SECTION("DiscreteFunctionP0 lhs")
          {
            DiscreteFunctionP0<Dimension, double> a{mesh};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                const double x = xj[cell_id][0];
                a[cell_id]     = 2 * x + 1;
              });

            Table<double> product_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  product_values[cell_id][i] = a[cell_id] * f[cell_id][i];
                }
              });

            REQUIRE(same_values(a * f, product_values));
            REQUIRE(same_values(a * const_f, product_values));

            DiscreteFunctionP0<Dimension, const double> const_a = a;
            REQUIRE(same_values(const_a * f, product_values));
            REQUIRE(same_values(const_a * const_f, product_values));
          }
        }
      }
    }

    SECTION("2D")
    {
      const size_t size = 3;

      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh2D();

      constexpr size_t Dimension = 2;

      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();

      SECTION("inner operators")
      {
        SECTION("scalar functions")
        {
          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
          parallel_for(
            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
              const double x = xj[cell_id][0];
              const double y = xj[cell_id][1];
              f[cell_id][0]  = 2 * x + 1;
              f[cell_id][1]  = x * x - y;
              f[cell_id][2]  = 2 + x * y;
            });

          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
          parallel_for(
            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
              const double x = xj[cell_id][0];
              const double y = xj[cell_id][1];
              g[cell_id][0]  = (x + 1) * (y - 2) + 1;
              g[cell_id][1]  = 3 * (x + 2) - y;
              g[cell_id][2]  = (x + 3) + 5 * y;
            });

          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};

          SECTION("sum")
          {
            Table<double> sum_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  sum_values[cell_id][i] = f[cell_id][i] + g[cell_id][i];
                }
              });

            REQUIRE(same_values(f + g, sum_values));
            REQUIRE(same_values(const_f + g, sum_values));
            REQUIRE(same_values(f + const_g, sum_values));
            REQUIRE(same_values(const_f + const_g, sum_values));
          }

          SECTION("difference")
          {
            Table<double> difference_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  difference_values[cell_id][i] = f[cell_id][i] - g[cell_id][i];
                }
              });

            REQUIRE(same_values(f - g, difference_values));
            REQUIRE(same_values(const_f - g, difference_values));
            REQUIRE(same_values(f - const_g, difference_values));
            REQUIRE(same_values(const_f - const_g, difference_values));
          }
        }
      }

      SECTION("external operators")
      {
        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            const double x = xj[cell_id][0];
            const double y = xj[cell_id][1];
            for (size_t i = 0; i < size; ++i) {
              f[cell_id][i] = std::abs(2 * x) + i * y;
            }
          });

        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;

        SECTION("product")
        {
          SECTION("scalar lhs")
          {
            const double a = 3.2;
            Table<double> product_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  product_values[cell_id][i] = a * f[cell_id][i];
                }
              });

            REQUIRE(same_values(a * f, product_values));
            REQUIRE(same_values(a * const_f, product_values));
          }

          SECTION("DiscreteFunctionP0 lhs")
          {
            DiscreteFunctionP0<Dimension, double> a{mesh};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                const double x = xj[cell_id][0];
                const double y = xj[cell_id][1];
                a[cell_id]     = 2 * x + 1 - y;
              });

            Table<double> product_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  product_values[cell_id][i] = a[cell_id] * f[cell_id][i];
                }
              });

            REQUIRE(same_values(a * f, product_values));
            REQUIRE(same_values(a * const_f, product_values));

            DiscreteFunctionP0<Dimension, const double> const_a = a;
            REQUIRE(same_values(const_a * f, product_values));
            REQUIRE(same_values(const_a * const_f, product_values));
          }
        }
      }
    }

    SECTION("3D")
    {
      const size_t size = 2;

      std::shared_ptr mesh = MeshDataBaseForTests::get().cartesianMesh3D();

      constexpr size_t Dimension = 3;

      auto xj = MeshDataManager::instance().getMeshData(*mesh).xj();

      SECTION("inner operators")
      {
        SECTION("scalar functions")
        {
          DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
          parallel_for(
            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
              const double x = xj[cell_id][0];
              const double y = xj[cell_id][1];
              const double z = xj[cell_id][2];
              f[cell_id][0]  = 2 * x * z + 1;
              f[cell_id][1]  = x * z - y;
            });

          DiscreteFunctionP0Vector<Dimension, double> g{mesh, size};
          parallel_for(
            mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
              const double x = xj[cell_id][0];
              const double y = xj[cell_id][1];
              const double z = xj[cell_id][2];
              g[cell_id][0]  = (x + 1) * (y - 2) + 1 - z;
              g[cell_id][1]  = 3 * (x + 2) - y * z;
            });

          DiscreteFunctionP0Vector<Dimension, const double> const_f = f;
          DiscreteFunctionP0Vector<Dimension, const double> const_g{g};

          SECTION("sum")
          {
            Table<double> sum_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  sum_values[cell_id][i] = f[cell_id][i] + g[cell_id][i];
                }
              });

            REQUIRE(same_values(f + g, sum_values));
            REQUIRE(same_values(const_f + g, sum_values));
            REQUIRE(same_values(f + const_g, sum_values));
            REQUIRE(same_values(const_f + const_g, sum_values));
          }

          SECTION("difference")
          {
            Table<double> difference_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  difference_values[cell_id][i] = f[cell_id][i] - g[cell_id][i];
                }
              });

            REQUIRE(same_values(f - g, difference_values));
            REQUIRE(same_values(const_f - g, difference_values));
            REQUIRE(same_values(f - const_g, difference_values));
            REQUIRE(same_values(const_f - const_g, difference_values));
          }
        }
      }

      SECTION("external operators")
      {
        DiscreteFunctionP0Vector<Dimension, double> f{mesh, size};
        parallel_for(
          mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
            const double x = xj[cell_id][0];
            const double y = xj[cell_id][1];
            const double z = xj[cell_id][2];
            for (size_t i = 0; i < size; ++i) {
              f[cell_id][i] = std::abs(2 * x) + i * y + z;
            }
          });

        DiscreteFunctionP0Vector<Dimension, const double> const_f = f;

        SECTION("product")
        {
          SECTION("scalar lhs")
          {
            const double a = 3.2;
            Table<double> product_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  product_values[cell_id][i] = a * f[cell_id][i];
                }
              });

            REQUIRE(same_values(a * f, product_values));
            REQUIRE(same_values(a * const_f, product_values));
          }

          SECTION("DiscreteFunctionP0 lhs")
          {
            DiscreteFunctionP0<Dimension, double> a{mesh};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                const double x = xj[cell_id][0];
                const double y = xj[cell_id][1];
                const double z = xj[cell_id][2];
                a[cell_id]     = 2 * x + 1 - y * z;
              });

            Table<double> product_values{mesh->numberOfCells(), size};
            parallel_for(
              mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) {
                for (size_t i = 0; i < size; ++i) {
                  product_values[cell_id][i] = a[cell_id] * f[cell_id][i];
                }
              });

            REQUIRE(same_values(a * f, product_values));
            REQUIRE(same_values(a * const_f, product_values));

            DiscreteFunctionP0<Dimension, const double> const_a = a;
            REQUIRE(same_values(const_a * f, product_values));
            REQUIRE(same_values(const_a * const_f, product_values));
          }
        }
      }
    }
  }
}