From 481785637e47466613ab74a366c4ad6970479ebc Mon Sep 17 00:00:00 2001
From: LABOURASSE Emmanuel <labourassee@gmail.com>
Date: Thu, 17 Sep 2020 09:13:17 +0200
Subject: [PATCH] implementation of template class for polynomials of degree>2
for ODE with GD
---
.../algorithms/ODEDiscontinuousGalerkin1D.cpp | 32 +++++++++++++------
src/language/modules/SchemeModule.cpp | 8 +++++
src/scheme/ODEGD1D.hpp | 12 ++++---
3 files changed, 38 insertions(+), 14 deletions(-)
diff --git a/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp b/src/language/algorithms/ODEDiscontinuousGalerkin1D.cpp
index 4bd10ddc0..9ed431f8a 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 a86a752c2..14f798eb6 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 c2d89c22b..19198305f 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;
}
--
GitLab