add the files to use Discontinuous Galerkin for 1D ODE

src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp

+#include <language/algorithms/ODEDiscontinuousGalerkin1D.hpp>
+#include <language/utils/InterpolateArray.hpp>
+#include <scheme/ODEGD1D.hpp>
+
+#include <fstream>
+
+template <size_t Dimension, size_t Degree>
+ODEDiscontinuousGalerkin1D<Dimension, Degree>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                                                                          std::shared_ptr<const IMesh> i_mesh,
+                                                                          const FunctionSymbolId& uex_id)
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+
+  std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+  static_assert(Dimension == 1, "Dimension have to be 1 for DiscontinuousGalerkin1D");
+  using R1 = TinyVector<Dimension>;
+  ODEDG1D<Degree> odedg1d(basis_type, mesh);
+  const NodeValue<const R1> xr    = mesh->xr();
+  const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+  Array<double> my_array;
+  for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+    const auto& cell_nodes          = cell_to_node_matrix[j];
+    const NodeId r1                 = cell_nodes[0];
+    const NodeId r2                 = cell_nodes[1];
+    const TinyVector<Dimension> xr1 = xr[r1];
+    const TinyVector<Dimension> xr2 = xr[r2];
+    double x13                      = (2 * xr1[0] + xr2[0]) / 3;
+    double x23                      = (xr1[0] + 2 * xr2[0]) / 3;
+    my_array[j * 4]                 = xr1[0];
+    my_array[j * 4 + 1]             = x13;
+    my_array[j * 4 + 2]             = x23;
+    my_array[j * 4 + 3]             = xr2[0];
+  }
+
+  Array<double> ej = InterpolateArray<double(double)>::interpolate(uex_id, my_array);
+
+  // InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(uex_id, mesh_data.xj());
+
+  CellValue<Polynomial<Degree>> U(mesh->connectivity());
+  odedg1d.globalSolve(U);
+  PolynomialBasis<Degree> B;
+  TinyVector<Degree + 1> zeros;
+  for (size_t i = 0; i <= Degree; i++) {
+    zeros[i] = i;
+  }
+  B.build(basis_type, 0.5, zeros);
+  std::string filename = "result-basis-";
+  filename += B.displayType();
+  filename += "-degree-";
+  filename += std::to_string(Degree);
+  filename += "-dof-";
+  filename += std::to_string(mesh->numberOfCells());
+  std::string filename2 = filename + "-polynomials";
+  std::ofstream sout(filename.c_str());
+  std::ofstream sout2(filename2.c_str());
+  for (CellId j = 0; j < mesh->numberOfCells(); ++j) {
+    const auto& cell_nodes          = cell_to_node_matrix[j];
+    const NodeId r1                 = cell_nodes[0];
+    const NodeId r2                 = cell_nodes[1];
+    const TinyVector<Dimension> xr1 = xr[r1];
+    const TinyVector<Dimension> xr2 = xr[r2];
+    const double deltax             = (xr2[0] - xr1[0]) / 100;
+    for (size_t i = 0; i <= 100; ++i) {
+      const double x = xr1[0] + deltax * i;
+      sout2 << x << " " << U[j](x) << '\n';
+    }
+    double x13 = (2 * xr1[0] + xr2[0]) / 3;
+    double x23 = (xr1[0] + 2 * xr2[0]) / 3;
+
+    double uex = std::exp(xr1[0]) + 3 * (std::exp(x13) + std::exp(x23)) + std::exp(xr2[0]);
+    uex /= 8;
+    double ucalc = integrate(U[j], xr1[0], xr2[0]) / (xr2[0] - xr1[0]);
+    sout << (0.5 * (xr1[0] + xr2[0])) << " " << uex << " " << ucalc << " " << std::abs(ucalc - uex) << '\n';
+  }
+}
+
+template ODEDiscontinuousGalerkin1D<1, 0>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                                                                      std::shared_ptr<const IMesh>,
+                                                                      const FunctionSymbolId&);
+template ODEDiscontinuousGalerkin1D<1, 1>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                                                                      std::shared_ptr<const IMesh>,
+                                                                      const FunctionSymbolId&);
+template ODEDiscontinuousGalerkin1D<1, 2>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                                                                      std::shared_ptr<const IMesh>,
+                                                                      const FunctionSymbolId&);

src/language/algorithms/ODEDiscontinuousGalerkin1D.hpp

+#ifndef ORDINARY_DIFFERENTIAL_EQUATION_DISCONTINUOUS_GALERKIN_1D_HPP
+#define ORDINARY_DIFFERENTIAL_EQUATION_DISCONTINUOUS_GALERKIN_1D_HPP
+
+#include <analysis/PolynomialBasis.hpp>
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+
+template <size_t Dimension, size_t Degree>
+struct ODEDiscontinuousGalerkin1D
+{
+  ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                             std::shared_ptr<const IMesh> i_mesh,
+                             const FunctionSymbolId& uex_id);
+};
+
+#endif   // ORDINARY_DIFFERENTIAL_EQUATION_DISCONTINUOUS_GALERKIN_1D_HPP

src/language/utils/InterpolateArray.hpp

+#ifndef INTERPOLATE_ARRAY_HPP
+#define INTERPOLATE_ARRAY_HPP
+
+#include <language/utils/PugsFunctionAdapter.hpp>
+
+template <typename T>
+class InterpolateArray;
+template <typename OutputType, typename InputType>
+class InterpolateArray<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+  static constexpr size_t Dimension = OutputType::Dimension;
+  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+  static inline Array<OutputType>
+  interpolate(const FunctionSymbolId& function_symbol_id, const Array<const InputType>& position)
+  {
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+    Array<OutputType> value(position.size());
+
+    parallel_for(position.size(), [=, &expression, &tokens](size_t i) {
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), position[i]);
+      auto result = expression.execute(execution_policy);
+      value[i]    = convert_result(std::move(result));
+
+      tokens.release(t);
+    });
+
+    return value;
+  }
+};
+
+#endif /* INTERPOLATE_ARRAY_HPP */

src/scheme/DiscontinuousGalerkin1DTools.hpp

+#ifndef DISCONTINUOUS_GALERKIN_1D_TOOLS
+#define DISCONTINUOUS_GALERKIN_1D_TOOLS
+
+#include <rang.hpp>
+
+#include <utils/ArrayUtils.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <analysis/Polynomial.hpp>
+#include <analysis/PolynomialBasis.hpp>
+
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+
+#include <utils/Exceptions.hpp>
+#include <utils/Messenger.hpp>
+
+#include <iostream>
+
+template <size_t Degree>
+class DiscontinuousGalerkin1DTools
+{
+  constexpr static size_t Dimension = 1;
+
+  using Rd  = TinyVector<Degree>;
+  using Rdd = TinyMatrix<Degree>;
+
+  using R1 = TinyVector<Dimension>;
+
+ public:
+  int
+  inverse() const
+  {
+    return 0;
+  }
+
+  void
+  elementarySolve(Polynomial<Degree>& U,
+                  const TinyVector<Degree + 1> moments,
+                  const PolynomialBasis<Degree> Basis,
+                  const double& xinf,
+                  const double& xsup) const
+  {
+    TinyMatrix<Degree + 1> M(zero);
+    for (size_t i = 0; i <= Degree; ++i) {
+      for (size_t j = 0; j <= Degree; ++j) {
+        Polynomial<2 * Degree> P = Basis.elements()[i] * Basis.elements()[j];
+        M(i, j)                  = integrate(P, xinf, xsup);
+      }
+    }
+    TinyMatrix inv_M = ::inverse(M);
+    U.coefficients() = inv_M * moments;
+  }
+
+  // void
+  // integrate() const
+  // {}
+
+ public:
+  // TODO add a flux manager to constructor
+  // TODO add a basis to constructor
+  DiscontinuousGalerkin1DTools()
+  {
+    ;
+  }
+};
+
+#endif   // DISCONTINUOUS_GALERKIN_1D_TOOLS

src/scheme/ODEGD1D.hpp

+#ifndef ODE_DG_1D
+#define ODE_DG_1D
+
+#include <rang.hpp>
+
+#include <utils/ArrayUtils.hpp>
+
+#include <utils/PugsAssert.hpp>
+
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+
+#include <analysis/Polynomial.hpp>
+#include <analysis/PolynomialBasis.hpp>
+
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/SubItemValuePerItem.hpp>
+
+#include <scheme/DiscontinuousGalerkin1DTools.hpp>
+#include <utils/Exceptions.hpp>
+#include <utils/Messenger.hpp>
+
+#include <iostream>
+enum class FluxType
+{
+  centered,
+  stable
+};
+
+template <size_t Degree>
+class ODEDG1D
+{
+  using MeshType                    = Mesh<Connectivity1D>;
+  constexpr static size_t Dimension = MeshType::Dimension;
+  static_assert(Dimension == 1, "Galerkin 1D only works in dimension 1");
+  using MeshDataType = MeshData<Dimension>;
+  //  using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
+
+  const BasisType& m_basis_type;
+  std::shared_ptr<const MeshType> m_mesh;
+  const typename MeshType::Connectivity& m_connectivity;
+
+  // class DirichletBoundaryCondition;
+
+  // using BoundaryCondition = std::variant<DirichletBoundaryCondition>;
+
+  // using BoundaryConditionList = std::vector<BoundaryCondition>;
+  // DiscontinuousGalerkin1DTools<Degree>& m_dgtools;
+  using Rd  = TinyVector<Degree>;
+  using Rdd = TinyMatrix<Degree>;
+
+  PolynomialBasis<Degree> m_basis;
+  using R1 = TinyVector<Dimension>;
+  // const BoundaryConditionList m_boundary_condition_list;
+  // NodeValue<R1> m_flux;
+  TinyVector<Degree + 1>
+  _buildZeros(double xmin, double xmax)
+  {
+    TinyVector<Degree + 1> zeros(zero);
+    zeros[0]  = xmin;
+    double dx = (xmax - xmin) / Degree;
+    for (size_t i = 1; i <= Degree; ++i) {
+      zeros[i] = xmin + i * dx;
+    }
+    return zeros;
+  }
+
+ public:
+  // void
+  // _applyBC()
+  // {
+  //   for (const auto& boundary_condition : m_boundary_condition_list) {
+  //     std::visit(
+  //       [&](auto&& bc) {
+  //         using T = std::decay_t<decltype(bc)>;
+  //         if constexpr (std::is_same_v<DirichletBoundaryCondition, T>) {
+  //           const auto& node_list  = bc.nodeList();
+  //           const auto& value_list = bc.valueList();
+
+  //           parallel_for(
+  //             node_list.size(), PUGS_LAMBDA(size_t i_node) {
+  //               NodeId node_id    = node_list[i_node];
+  //               const auto& value = value_list[i_node];
+
+  //               m_flux[node_id] = value;
+  //             });
+  //         }
+  //       },
+  //       boundary_condition);
+  //   }
+  // }
+
+  // NodeValue<R1>
+  // computeFluxes(FluxType flux_type, const PolynomialBasis<Degree> Basis, Polynomial<Degree>& U)
+  // {
+  //   switch (flux_type) {
+  //   // case FluxType::centered: {
+  //   //   return computeCenteredFlux(Basis, U);
+  //   //   break;
+  //   // }
+  //   case FluxType::stable: {
+  //     return computeStableFlux(Basis, U);
+  //     break;
+  //   }
+  //   default:
+  //     throw UnexpectedError("unknown flux type");
+  //   }
+  // }
+
+  double
+  computeStableFluxes(const CellValue<const Polynomial<Degree>>& U, NodeId r)
+  {
+    double flux;
+    const NodeValue<const R1> xr    = m_mesh->xr();
+    const auto& node_to_cell_matrix = m_mesh->connectivity().nodeToCellMatrix();
+
+    const auto& node_cells = node_to_cell_matrix[r];
+    if (node_cells.size() == 1) {
+      flux = 1.;
+    } else {
+      const CellId j1  = node_cells[0];
+      const double xr0 = xr[r][0];
+      flux             = U[j1](xr0);
+      std::cout << "U[" << j1 << "]=" << U[j1] << "\n";
+    }
+    return flux;
+  }
+
+  void
+  globalSolve(CellValue<Polynomial<Degree>>& U)
+  {
+    const NodeValue<const R1> xr    = m_mesh->xr();
+    const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
+    CellValue<TinyVector<Degree + 1>> moments(m_connectivity);
+    for (CellId j = 0; j < m_mesh->numberOfCells(); j++) {
+      const auto& cell_nodes             = cell_to_node_matrix[j];
+      const NodeId r1                    = cell_nodes[0];
+      const NodeId r2                    = cell_nodes[1];
+      const TinyVector<Dimension> xr1    = xr[r1];
+      const TinyVector<Dimension> xr2    = xr[r2];
+      const TinyVector<Degree + 1> zeros = _buildZeros(xr1[0], xr2[0]);
+      std::cout << "x1 " << xr1[0] << " x2 " << xr2[0] << " zeros " << zeros << "\n";
+      m_basis.build(m_basis_type, 0.5 * (xr1[0] + xr2[0]), zeros);
+      std::cout << " B[" << j << "]=" << m_basis.elements()[0] << "\n";
+      const double fluxg = computeStableFluxes(U, r1);
+      std::cout << j << " flux " << fluxg << "\n";
+      moments[j] = buildMoments(m_basis, fluxg, xr1[0]);
+      std::cout << " moments[" << j << "]=" << moments[j] << "\n";
+
+      U[j] = elementarySolve(moments[j], m_basis, xr1[0], xr2[0]);
+      std::cout << " U[" << j << "]=" << U[j] << "\n";
+    }
+  }
+
+  // void
+  // integrate() const
+  // {}
+  PUGS_INLINE constexpr TinyVector<Degree + 1>
+  buildMoments(PolynomialBasis<Degree> basis, double fluxg, double xr1)
+  {
+    TinyVector<Degree + 1> moments(zero);
+    for (size_t i = 0; i <= Degree; i++) {
+      // moments[i] = (basis.elements()[i](xr2) - basis.elements()[i](xr1)) * fluxg;
+      moments[i] = -basis.elements()[i](xr1) * fluxg;
+    }
+    return moments;
+  }
+
+  PUGS_INLINE constexpr Polynomial<Degree>
+  elementarySolve(const TinyVector<Degree + 1> moments,
+                  const PolynomialBasis<Degree> Basis,
+                  const double& xinf,
+                  const double& xsup) const
+  {
+    Polynomial<Degree> U;
+    U.coefficients() = zero;
+    TinyMatrix<Degree + 1> M(zero);
+    for (size_t i = 0; i <= Degree; ++i) {
+      for (size_t j = 0; j <= Degree; ++j) {
+        Polynomial<2 * Degree> P = Basis.elements()[j] * Basis.elements()[i];
+        P += derivative(Basis.elements()[i]) * Basis.elements()[j];
+        double fluxd = Basis.elements()[i](xsup) * Basis.elements()[j](xsup);
+        M(i, j)      = integrate(P, xinf, xsup) - fluxd;
+      }
+    }
+    std::cout << " M " << M << "\n";
+    TinyMatrix inv_M                    = ::inverse(M);
+    TinyVector<Degree + 1> coefficients = inv_M * moments;
+    for (size_t i = 0; i <= Degree; ++i) {
+      U += coefficients[i] * Basis.elements()[i];
+    }
+    std::cout << "U(" << xsup << ")=" << U(xsup) << "\n";
+    return U;
+  }
+
+ public:
+  // TODO add a flux manager to constructor
+  // TODO add a basis to constructor
+  ODEDG1D(const BasisType& basis_type, std::shared_ptr<const IMesh> i_mesh)   //, const BoundaryConditionList& bc_list)
+  : m_basis_type(basis_type),
+    m_mesh(std::dynamic_pointer_cast<const MeshType>(i_mesh)),
+    m_connectivity(m_mesh->connectivity())   //,
+  // m_boundary_condition_list(bc_list)
+  {
+    ;
+  }
+};
+
+#endif   // ODE_DG_1D

tests/test_DiscontinuousGalerkin1D.cpp

+#include <catch2/catch.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <utils/PugsAssert.hpp>
+#include <utils/Types.hpp>
+
+#include <algebra/TinyMatrix.hpp>
+#include <analysis/Polynomial.hpp>
+#include <analysis/PolynomialBasis.hpp>
+#include <mesh/CartesianMeshBuilder.hpp>
+#include <mesh/DiamondDualConnectivityManager.hpp>
+#include <mesh/DiamondDualMeshManager.hpp>
+#include <scheme/DiscontinuousGalerkin1DTools.hpp>
+#include <scheme/ODEGD1D.hpp>
+
+// Instantiate to ensure full coverage is performed
+template class Polynomial<0>;
+template class Polynomial<1>;
+template class Polynomial<2>;
+template class Polynomial<3>;
+template class Polynomial<4>;
+template class Polynomial<5>;
+
+// clazy:excludeall=non-pod-global-static
+
+TEST_CASE("Discontinuous-Galerkin-1D", "[discontinuous-galerkin-1d]")
+{
+  SECTION("Elementary-tests")
+  {
+    TinyVector<1> a{0};
+    TinyVector<1> b{1};
+    TinyVector<1, uint64_t> nbcells{10};
+    // CartesianMeshBuilder mesh_builder = CartesianMeshBuilder(a, b, nbcells);
+    // auto mesh                         = mesh_builder.mesh();
+    // DiscontinuousGalerkin1D<2>{mesh};
+    REQUIRE_NOTHROW(DiscontinuousGalerkin1DTools<2>());
+    DiscontinuousGalerkin1DTools dg = DiscontinuousGalerkin1DTools<1>();
+    Polynomial<1> U;
+    TinyVector<2> F({3, 2});
+    PolynomialBasis<1> Basis;
+    Basis.build(BasisType::canonical);
+    dg.elementarySolve(U, F, Basis, -1, 1);
+    REQUIRE(U == Polynomial<1>{{3. / 2, 3}});
+    Polynomial<2> U2;
+    TinyVector<3> F2({1, 1, 1});
+    PolynomialBasis<2> Basis2;
+    Basis2.build(BasisType::canonical);
+    DiscontinuousGalerkin1DTools dg2 = DiscontinuousGalerkin1DTools<2>();
+    dg2.elementarySolve(U2, F2, Basis2, -1, 1);
+    REQUIRE(integrate(U2, -1, 1) == Approx(1).epsilon(1E-14));
+    REQUIRE(integrate(U2 * Polynomial<2>{{0, 0, 1}}, -1, 1) == 1.);
+    dg2.elementarySolve(U2, F2, Basis2, 0, 1);
+    REQUIRE(integrate(U2 * Polynomial<2>{{0, 0, 1}}, 0, 1) == Approx(1).epsilon(1E-14));
+  }
+}