diff --git a/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp b/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp index 4bd10ddc02c1830ff05cc30341e48a10574f734f..9ed431f8aef454ed751fc52f6d4d84df97172c87 100644 --- a/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp +++ b/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&); diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index a86a752c2dbb3d7be4304e638a3f7c269311bcb1..14f798eb6bac43b77fd98367a009fa19d5dcfcbd 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/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"); } diff --git a/src/scheme/ODEGD1D.hpp b/src/scheme/ODEGD1D.hpp index c2d89c22ba3576b9b825ca15ec79b77ecb40d9e1..19198305f6d80f183c764503cdc7717e54ae9b0e 100644 --- a/src/scheme/ODEGD1D.hpp +++ b/src/scheme/ODEGD1D.hpp @@ -62,10 +62,14 @@ class ODEDG1D _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; + 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; }