Skip to content
Snippets Groups Projects
Select Git revision
  • d9f7577649f69b8d87daf2517a9a1e0f4c4bd0ed
  • develop default protected
  • feature/gmsh-reader
  • origin/stage/bouguettaia
  • feature/kinetic-schemes
  • feature/reconstruction
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • feature/serraille
  • feature/variational-hydro
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • master protected
  • feature/escobar-smoother
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

BuiltinFunctionProcessor.hpp

Blame
    • Stéphane Del Pino's avatar
      d9f75776
      Improve exception handling · d9f75776
      Stéphane Del Pino authored
      Backtrace is now only produced when '-p' option is provided
      
      This gives better user information and allows to debug more simply
      since exceptions are no more rethrown
      
      Going further would require to be able to indicate the current line of
      execution during a non exception crash
      d9f75776
      History
      Improve exception handling
      Stéphane Del Pino authored
      Backtrace is now only produced when '-p' option is provided
      
      This gives better user information and allows to debug more simply
      since exceptions are no more rethrown
      
      Going further would require to be able to indicate the current line of
      execution during a non exception crash
    test_DiscreteFunctionP0Vector.cpp 26.59 KiB
    #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));
              }
            }
          }
        }
      }
    }