implementation of template class for polynomials of degree>2 for ODE with GD

src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp

@@ -18,22 +18,24 @@ ODEDiscontinuousGalerkin1D<Dimension, Degree>::ODEDiscontinuousGalerkin1D(const
   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;
+  Array<double> my_array(100 * mesh->numberOfCells() + 1);
+
   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];
+    const double deltax             = (xr2[0] - xr1[0]) / 100;
+    bool start                      = (j == 0);
+    for (size_t i = (!start); i <= 100; ++i) {
+      const double x        = xr1[0] + deltax * i;
+      my_array[100 * j + i] = x;
+    }
   }
 
   Array<double> ej = InterpolateArray<double(double)>::interpolate(uex_id, my_array);
+  // Array<TinyVector<Dimension>> ej = InterpolateArray<TinyVector<Dimension>(double)>::interpolate(uex_id, my_array);
 
   // InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(uex_id, mesh_data.xj());
 
@@ -61,9 +63,12 @@ ODEDiscontinuousGalerkin1D<Dimension, Degree>::ODEDiscontinuousGalerkin1D(const
     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) {
+    bool start                      = (j == 0);
+    double norm1error               = 0.;
+    for (size_t i = (!start); i <= 100; ++i) {
       const double x = xr1[0] + deltax * i;
-      sout2 << x << " " << U[j](x) << '\n';
+      sout2 << x << " " << U[j](x) << " " << ej[100 * j + i] << " " << std::abs(U[j](x) - ej[100 * j + i]) << '\n';
+      norm1error += std::abs(U[j](x) - ej[100 * j + i]);
     }
     double x13 = (2 * xr1[0] + xr2[0]) / 3;
     double x23 = (xr1[0] + 2 * xr2[0]) / 3;
@@ -71,7 +76,8 @@ ODEDiscontinuousGalerkin1D<Dimension, Degree>::ODEDiscontinuousGalerkin1D(const
     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';
+    sout << (0.5 * (xr1[0] + xr2[0])) << " " << uex << " " << ucalc << " " << std::abs(ucalc - uex) << " " << norm1error
+         << '\n';
   }
 }
 
@@ -84,3 +90,9 @@ template ODEDiscontinuousGalerkin1D<1, 1>::ODEDiscontinuousGalerkin1D(const Basi
 template ODEDiscontinuousGalerkin1D<1, 2>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
                                                                       std::shared_ptr<const IMesh>,
                                                                       const FunctionSymbolId&);
+template ODEDiscontinuousGalerkin1D<1, 3>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                                                                      std::shared_ptr<const IMesh>,
+                                                                      const FunctionSymbolId&);
+template ODEDiscontinuousGalerkin1D<1, 4>::ODEDiscontinuousGalerkin1D(const BasisType& basis_type,
+                                                                      std::shared_ptr<const IMesh>,
+                                                                      const FunctionSymbolId&);

src/language/modules/SchemeModule.cpp

@@ -223,6 +223,14 @@ SchemeModule::SchemeModule()
                                   ODEDiscontinuousGalerkin1D<1, 2>(*basis_type, p_mesh, uex_id);
                                   break;
                                 }
+                                case 3: {
+                                  ODEDiscontinuousGalerkin1D<1, 3>(*basis_type, p_mesh, uex_id);
+                                  break;
+                                }
+                                case 4: {
+                                  ODEDiscontinuousGalerkin1D<1, 4>(*basis_type, p_mesh, uex_id);
+                                  break;
+                                }
                                 default: {
                                   throw UnexpectedError("invalid polynomial degree");
                                 }

src/scheme/ODEGD1D.hpp

@@ -62,11 +62,15 @@ class ODEDG1D
   _buildZeros(double xmin, double xmax)
   {
     TinyVector<Degree + 1> zeros(zero);
+    if constexpr (Degree == 0) {
+      zeros[0] = 0.5 * (xmax - xmin);
+    } else {
       zeros[0]  = xmin;
       double dx = (xmax - xmin) / Degree;
       for (size_t i = 1; i <= Degree; ++i) {
         zeros[i] = xmin + i * dx;
       }
+    }
     return zeros;
   }