From 6c591d2084ac23172480a3ef386e768cea901ee5 Mon Sep 17 00:00:00 2001 From: Jan Breuer <74359367+jbreue16@users.noreply.github.com> Date: Tue, 16 Jan 2024 14:27:20 +0100 Subject: [PATCH] Move DG toolbox into a separate file --- src/libcadet/CMakeLists.txt | 3 +- src/libcadet/model/GeneralRateModelDG.cpp | 12 +- src/libcadet/model/GeneralRateModelDG.hpp | 184 +------- .../parts/ConvectionDispersionOperatorDG.cpp | 8 +- .../parts/ConvectionDispersionOperatorDG.hpp | 347 +--------------- src/libcadet/model/parts/DGToolbox.cpp | 393 ++++++++++++++++++ src/libcadet/model/parts/DGToolbox.hpp | 84 ++++ 7 files changed, 495 insertions(+), 536 deletions(-) create mode 100644 src/libcadet/model/parts/DGToolbox.cpp create mode 100644 src/libcadet/model/parts/DGToolbox.hpp diff --git a/src/libcadet/CMakeLists.txt b/src/libcadet/CMakeLists.txt index 015318228..6a2b7ca79 100644 --- a/src/libcadet/CMakeLists.txt +++ b/src/libcadet/CMakeLists.txt @@ -139,7 +139,6 @@ set(LIBCADET_SOURCES ${CMAKE_SOURCE_DIR}/src/libcadet/model/GeneralRateModelBuilder.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/ParameterMultiplexing.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperator.cpp - ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/AxialConvectionDispersionKernel.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/RadialConvectionDispersionKernel.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/reaction/ReactionModelBase.cpp @@ -196,6 +195,8 @@ configure_file("${CMAKE_CURRENT_SOURCE_DIR}/linalg/SparseSolverInterface.hpp.in" if (ENABLE_DG) list (APPEND LIBCADET_SOURCES + ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/DGToolbox.cpp + ${CMAKE_SOURCE_DIR}/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithPoresDG.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/LumpedRateModelWithoutPoresDG.cpp ${CMAKE_SOURCE_DIR}/src/libcadet/model/GeneralRateModelDG.cpp diff --git a/src/libcadet/model/GeneralRateModelDG.cpp b/src/libcadet/model/GeneralRateModelDG.cpp index 5d5557c65..46688a30a 100644 --- a/src/libcadet/model/GeneralRateModelDG.cpp +++ b/src/libcadet/model/GeneralRateModelDG.cpp @@ -2045,26 +2045,26 @@ void GeneralRateModelDG::updateRadialDisc() // compute mass matrices for exact integration based on particle geometry, via transformation to normalized Jacobi polynomials with weight function w if (_parGeomSurfToVol[parType] == SurfVolRatioSphere) { // w = (1 + \xi)^2 - _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.invMMatrix(_disc.nParNode[parType], _disc.parNodes[parType], 0.0, 2.0).inverse() * pow((static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0), 2.0); + _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 2.0).inverse() * pow((static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0), 2.0); if(cell > 0 || _parCoreRadius[parType] != 0.0) // following contributions are zero for first cell when R_c = 0 (no particle core) - _disc.parInvMM[_disc.offsetMetric[parType] + cell] += _disc.invMMatrix(_disc.nParNode[parType], _disc.parNodes[parType], 0.0, 1.0).inverse() * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) * static_cast(r_L)) - + _disc.invMMatrix(_disc.nParNode[parType], _disc.parNodes[parType], 0.0, 0.0).inverse() * pow(static_cast(r_L), 2.0); + _disc.parInvMM[_disc.offsetMetric[parType] + cell] += parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 1.0).inverse() * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) * static_cast(r_L)) + + parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0).inverse() * pow(static_cast(r_L), 2.0); _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); _disc.minus_InvMM_ST[_disc.offsetMetric[parType] + cell] = - _disc.parInvMM[_disc.offsetMetric[parType] + cell] * _disc.parPolyDerM[parType].transpose() * _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); } else if (_parGeomSurfToVol[parType] == SurfVolRatioCylinder) { // w = (1 + \xi) - _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.invMMatrix(_disc.nParNode[parType], _disc.parNodes[parType], 0.0, 1.0).inverse() * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0); + _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 1.0).inverse() * (static_cast(_disc.deltaR[_disc.offsetMetric[parType] + cell]) / 2.0); if (cell > 0 || _parCoreRadius[parType] != 0.0) // following contribution is zero for first cell when R_c = 0 (no particle core) - _disc.parInvMM[_disc.offsetMetric[parType] + cell] += _disc.invMMatrix(_disc.nParNode[parType], _disc.parNodes[parType], 0.0, 0.0).inverse() * static_cast(r_L); + _disc.parInvMM[_disc.offsetMetric[parType] + cell] += parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0).inverse() * static_cast(r_L); _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); _disc.minus_InvMM_ST[_disc.offsetMetric[parType] + cell] = -_disc.parInvMM[_disc.offsetMetric[parType] + cell] * _disc.parPolyDerM[parType].transpose() * _disc.parInvMM[_disc.offsetMetric[parType] + cell].inverse(); } else if (_parGeomSurfToVol[parType] == SurfVolRatioSlab) { // w = 1 - _disc.parInvMM[_disc.offsetMetric[parType] + cell] = _disc.invMMatrix(_disc.nParNode[parType], _disc.parNodes[parType], 0.0, 0.0); + _disc.parInvMM[_disc.offsetMetric[parType] + cell] = parts::dgtoolbox::invMMatrix(_disc.parPolyDeg[parType], _disc.parNodes[parType], 0.0, 0.0); } } } diff --git a/src/libcadet/model/GeneralRateModelDG.hpp b/src/libcadet/model/GeneralRateModelDG.hpp index a98c735b5..b593afe67 100644 --- a/src/libcadet/model/GeneralRateModelDG.hpp +++ b/src/libcadet/model/GeneralRateModelDG.hpp @@ -362,9 +362,9 @@ class GeneralRateModelDG : public UnitOperationBase for (int parType = 0; parType < nParType; parType++) { - lglNodesWeights(parPolyDeg[parType], parNodes[parType], parInvWeights[parType]); - derivativeMatrix(parPolyDeg[parType], parPolyDerM[parType], parNodes[parType]); - parInvMM_Leg[parType] = invMMatrix(nParNode[parType], parNodes[parType], 0.0, 0.0); + parts::dgtoolbox::lglNodesWeights(parPolyDeg[parType], parNodes[parType], parInvWeights[parType], true); + parPolyDerM[parType] = parts::dgtoolbox::derivativeMatrix(parPolyDeg[parType], parNodes[parType]); + parInvMM_Leg[parType] = parts::dgtoolbox::invMMatrix(parPolyDeg[parType], parNodes[parType], 0.0, 0.0); } } @@ -382,175 +382,6 @@ class GeneralRateModelDG : public UnitOperationBase private: - /* =================================================================================== - * Polynomial Basis operators and auxiliary functions - * =================================================================================== */ - - /** - * @brief computes the Legendre polynomial L_N and q = L_N+1 - L_N-2 and q' at point x - * @param [in] polyDeg polynomial degree of spatial Discretization - * @param [in] x evaluation point - * @param [in] L <- L(x) - * @param [in] q <- q(x) = L_N+1 (x) - L_N-2(x) - * @param [in] qder <- q'(x) = [L_N+1 (x) - L_N-2(x)]' - */ - void qAndL(const int _polyDeg, const double x, double& L, double& q, double& qder) { - - double L_2 = 1.0; - double L_1 = x; - double Lder_2 = 0.0; - double Lder_1 = 1.0; - double Lder = 0.0; - - for (double k = 2; k <= _polyDeg; k++) { // note that this function is only called for polyDeg >= 2. - L = ((2 * k - 1) * x * L_1 - (k - 1) * L_2) / k; - Lder = Lder_2 + (2 * k - 1) * L_1; - L_2 = L_1; - L_1 = L; - Lder_2 = Lder_1; - Lder_1 = Lder; - } - q = ((2.0 * _polyDeg + 1) * x * L - _polyDeg * L_2) / (_polyDeg + 1.0) - L_2; - qder = Lder_1 + (2.0 * _polyDeg + 1) * L_1 - Lder_2; - } - - /** - * @brief computes the Legendre-Gauss-Lobatto nodes and (inverse) quadrature weights - * @detail inexact LGL-quadrature leads to a diagonal mass matrix (mass lumping), defined by the quadrature weights - */ - void lglNodesWeights(const int _polyDeg, Eigen::VectorXd& _nodes, Eigen::VectorXd& _invWeights) { - - const double pi = 3.1415926535897932384626434; - - // tolerance and max #iterations for Newton iteration - int nIterations = 10; - double tolerance = 1e-15; - - // Legendre polynomial and derivative - double L = 0; - double q = 0; - double qder = 0; - - switch (_polyDeg) { - case 0: - throw std::invalid_argument("Polynomial degree must be at least 1 !"); - break; - case 1: - _nodes[0] = -1; - _invWeights[0] = 1; - _nodes[1] = 1; - _invWeights[1] = 1; - break; - default: - _nodes[0] = -1; - _nodes[_polyDeg] = 1; - _invWeights[0] = 2.0 / (_polyDeg * (_polyDeg + 1.0)); - _invWeights[_polyDeg] = _invWeights[0]; - - // use symmetrie, only compute half of points and weights - for (unsigned int j = 1; j <= floor((_polyDeg + 1) / 2) - 1; j++) { - // first guess for Newton iteration - _nodes[j] = -cos(pi * (j + 0.25) / _polyDeg - 3 / (8.0 * _polyDeg * pi * (j + 0.25))); - // Newton iteration to find roots of Legendre Polynomial - for (unsigned int k = 0; k <= nIterations; k++) { - qAndL(_polyDeg, _nodes[j], L, q, qder); - _nodes[j] = _nodes[j] - q / qder; - if (abs(q / qder) <= tolerance * abs(_nodes[j])) { - break; - } - } - // calculate weights - qAndL(_polyDeg, _nodes[j], L, q, qder); - _invWeights[j] = 2.0 / (_polyDeg * (_polyDeg + 1.0) * pow(L, 2.0)); - _nodes[_polyDeg - j] = -_nodes[j]; // copy to second half of points and weights - _invWeights[_polyDeg - j] = _invWeights[j]; - } - } - - if (_polyDeg % 2 == 0) { // for even polyDeg we have an odd number of points which include 0.0 - qAndL(_polyDeg, 0.0, L, q, qder); - _nodes[_polyDeg / 2] = 0; - _invWeights[_polyDeg / 2] = 2.0 / (_polyDeg * (_polyDeg + 1.0) * pow(L, 2.0)); - } - - _invWeights = _invWeights.cwiseInverse(); // we need the inverse of the weights - } - - /** - * @brief computation of barycentric weights for fast polynomial evaluation - * @param [in] baryWeights vector to store barycentric weights. Must already be initialized with ones! - */ - void barycentricWeights(Eigen::VectorXd& baryWeights, const Eigen::VectorXd& _nodes, const int _polyDeg) { - - for (unsigned int j = 1; j <= _polyDeg; j++) { - for (unsigned int k = 0; k <= j - 1; k++) { - baryWeights[k] = baryWeights[k] * (_nodes[k] - _nodes[j]) * 1.0; - baryWeights[j] = baryWeights[j] * (_nodes[j] - _nodes[k]) * 1.0; - } - } - - for (unsigned int j = 0; j <= _polyDeg; j++) { - baryWeights[j] = 1 / baryWeights[j]; - } - } - - /** - * @brief computation of nodal (lagrange) polynomial derivative matrix - */ - void derivativeMatrix(const int _polyDeg, Eigen::MatrixXd& _polyDerM, const Eigen::VectorXd& _nodes) { - - Eigen::VectorXd baryWeights = Eigen::VectorXd::Ones(_polyDeg + 1u); - barycentricWeights(baryWeights, _nodes, _polyDeg); - - for (unsigned int i = 0; i <= _polyDeg; i++) { - for (unsigned int j = 0; j <= _polyDeg; j++) { - if (i != j) { - _polyDerM(i, j) = baryWeights[j] / (baryWeights[i] * (_nodes[i] - _nodes[j])); - _polyDerM(i, i) += -_polyDerM(i, j); - } - } - } - } - - /** - * @brief factor to normalize legendre polynomials - */ - double orthonFactor(const int _polyDeg, double a, double b) { - - double n = static_cast (_polyDeg); - return std::sqrt(((2.0 * n + a + b + 1.0) * std::tgamma(n + 1.0) * std::tgamma(n + a + b + 1.0)) - / (std::pow(2.0, a + b + 1.0) * std::tgamma(n + a + 1.0) * std::tgamma(n + b + 1.0))); - } - /** - * @brief calculates the Vandermonde matrix of the normalized legendre polynomials - */ - Eigen::MatrixXd getVandermonde_JACOBI(const int _nNodes, const Eigen::VectorXd _nodes, double a, double b) { - - Eigen::MatrixXd V(_nNodes, _nNodes); - - // degree 0 - V.block(0, 0, _nNodes, 1) = VectorXd::Ones(_nNodes) * orthonFactor(0, a, b); - // degree 1 - for (int node = 0; node < static_cast(_nNodes); node++) { - V(node, 1) = ((_nodes[node] - 1.0) / 2.0 * (a + b + 2.0) + (a + 1.0)) * orthonFactor(1, a, b); - } - - for (int deg = 2; deg <= static_cast(_nNodes - 1); deg++) { - - for (int node = 0; node < static_cast(_nNodes); node++) { - - double orthn_1 = orthonFactor(deg, a, b) / orthonFactor(deg - 1, a, b); - double orthn_2 = orthonFactor(deg, a, b) / orthonFactor(deg - 2, a, b); - - // recurrence relation - V(node, deg) = orthn_1 * ((2.0 * deg + a + b - 1.0) * ((2.0 * deg + a + b) * (2.0 * deg + a + b - 2.0) * _nodes[node] + a * a - b * b) * V(node, deg - 1)); - V(node, deg) -= orthn_2 * (2.0 * (deg + a - 1.0) * (deg + b - 1.0) * (2.0 * deg + a + b) * V(node, deg - 2)); - V(node, deg) /= 2.0 * deg * (deg + a + b) * (2.0 * deg + a + b - 2.0); - } - } - - return V; - } /** * @brief calculates the DG Jacobian auxiliary block * @param [in] exInt true if exact integration DG scheme @@ -663,15 +494,6 @@ class GeneralRateModelDG : public UnitOperationBase return -dispBlock; // *-1 for residual } - public: - /** - * @brief calculates the inverse mass matrix for exact integration - * @detail calculation via transformation of Lagrange to Jacobi polynomial coefficients - */ - MatrixXd invMMatrix(const int _nnodes, const Eigen::VectorXd _nodes, double alpha, double beta) { - return (getVandermonde_JACOBI(_nnodes, _nodes, alpha, beta) * (getVandermonde_JACOBI(_nnodes, _nodes, alpha, beta).transpose())); - } - }; enum class ParticleDiscretizationMode : int diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp index ce0a6f629..3a0666faa 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.cpp @@ -102,12 +102,12 @@ bool AxialConvectionDispersionOperatorBaseDG::configureModelDiscretization(IPara _newStaticJac = true; - lglNodesWeights(); - invMMatrix(); - derivativeMatrix(); + dgtoolbox::lglNodesWeights(_polyDeg, _nodes, _invWeights, true); + _invMM = dgtoolbox::invMMatrix(_polyDeg, _nodes); + _polyDerM = dgtoolbox::derivativeMatrix(_polyDeg, _nodes); if(polynomial_integration_mode == 2) // use Gauss quadrature for exact integration - _invMM = gaussQuadratureMMatrix(_nodes, _nNodes).inverse(); + _invMM = dgtoolbox::gaussQuadratureMMatrix(_nodes, _nNodes).inverse(); if (paramProvider.exists("COL_DISPERSION_DEP")) { diff --git a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp index f5a4a1b8a..b98e8b296 100644 --- a/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp +++ b/src/libcadet/model/parts/ConvectionDispersionOperatorDG.hpp @@ -25,6 +25,7 @@ #include "SimulationTypes.hpp" #include #include "linalg/BandedEigenSparseRowIterator.hpp" +#include "model/parts/DGToolbox.hpp" #include #include @@ -32,11 +33,6 @@ #include #include -#ifndef _USE_MATH_DEFINES -#define _USE_MATH_DEFINES -#include -#endif - using namespace Eigen; namespace cadet @@ -197,342 +193,9 @@ namespace cadet IParameterParameterDependence* _dispersionDep; /* =================================================================================== - * Polynomial Basis operators and auxiliary functions - * =================================================================================== */ - - /** - * @brief computes the Legendre polynomial L_N and q = L_N+1 - L_N-2 and q' at point x - * @param [in] x evaluation point - * @param [in] L <- L(x) - * @param [in] q <- q(x) = L_N+1 (x) - L_N-2(x) - * @param [in] qder <- q'(x) = [L_N+1 (x) - L_N-2(x)]' - */ - void qAndL(const double x, double& L, double& q, double& qder) { - // auxiliary variables (Legendre polynomials) - double L_2 = 1.0; - double L_1 = x; - double Lder_2 = 0.0; - double Lder_1 = 1.0; - double Lder = 0.0; - for (double k = 2; k <= _polyDeg; k++) { // note that this function is only called for polyDeg >= 2. - L = ((2 * k - 1) * x * L_1 - (k - 1) * L_2) / k; - Lder = Lder_2 + (2 * k - 1) * L_1; - L_2 = L_1; - L_1 = L; - Lder_2 = Lder_1; - Lder_1 = Lder; - } - q = ((2.0 * _polyDeg + 1) * x * L - _polyDeg * L_2) / (_polyDeg + 1.0) - L_2; - qder = Lder_1 + (2.0 * _polyDeg + 1) * L_1 - Lder_2; - } - - /** - * @brief computes the Legendre-Gauss-Lobatto nodes and (inverse) quadrature weights - */ - void lglNodesWeights() { - // tolerance and max #iterations for Newton iteration - int nIterations = 10; - double tolerance = 1e-15; - // Legendre polynomial and derivative - double L = 0; - double q = 0; - double qder = 0; - switch (_polyDeg) { - case 0: - throw std::invalid_argument("Polynomial degree must be at least 1 !"); - break; - case 1: - _nodes[0] = -1; - _invWeights[0] = 1; - _nodes[1] = 1; - _invWeights[1] = 1; - break; - default: - _nodes[0] = -1; - _nodes[_polyDeg] = 1; - _invWeights[0] = 2.0 / (_polyDeg * (_polyDeg + 1.0)); - _invWeights[_polyDeg] = _invWeights[0]; - // use symmetrie, only compute half of points and weights - for (unsigned int j = 1; j <= floor((_polyDeg + 1) / 2) - 1; j++) { - // first guess for Newton iteration - _nodes[j] = -cos(M_PI * (j + 0.25) / _polyDeg - 3 / (8.0 * _polyDeg * M_PI * (j + 0.25))); - // Newton iteration to find roots of Legendre Polynomial - for (unsigned int k = 0; k <= nIterations; k++) { - qAndL(_nodes[j], L, q, qder); - _nodes[j] = _nodes[j] - q / qder; - if (abs(q / qder) <= tolerance * abs(_nodes[j])) { - break; - } - } - // calculate weights - qAndL(_nodes[j], L, q, qder); - _invWeights[j] = 2.0 / (_polyDeg * (_polyDeg + 1.0) * pow(L, 2.0)); - _nodes[_polyDeg - j] = -_nodes[j]; // copy to second half of points and weights - _invWeights[_polyDeg - j] = _invWeights[j]; - } - } - if (_polyDeg % 2 == 0) { // for even _polyDeg we have an odd number of points which include 0.0 - qAndL(0.0, L, q, qder); - _nodes[_polyDeg / 2] = 0; - _invWeights[_polyDeg / 2] = 2.0 / (_polyDeg * (_polyDeg + 1.0) * pow(L, 2.0)); - } - // inverse the weights - _invWeights = _invWeights.cwiseInverse(); - } - - /** - * @brief computes the Legendre polynomial and its derivative - * @param [in, out] leg Legendre polynomial - * @param [in, out] leg Legendre polynomial derivative - * @param [in] N polynomial degree - * @param [in] x evaluation point - */ - void legendrePolynomialAndDerivative(double& leg, double& legDer, const int N, const double x) { - - switch (N) { - case 0: - leg = 1.0; - legDer = 0.0; - break; - case 1: - leg = x; - legDer = 1.0; - break; - default: - double leg_2 = 1.0; - double leg_1 = x; - double legDer_2 = 0.0; - double legDer_1 = 1.0; - - for (int k = 2; k <= N; k++) { - leg = (2.0 * k - 1.0) / k * x * leg_1 - (k - 1.0) / k * leg_2; - legDer = legDer_2 + (2.0 * k - 1.0) * leg_1; - leg_2 = leg_1; - leg_1 = leg; - legDer_2 = legDer_1; - legDer_1 = legDer; - } - } - } - - /** - * @brief computes the Legendre-Gauss nodes and quadrature weights - * @param [in, out] nodes Legendre Gauss nodes - * @param [in, out] weights Legendre Gauss quadrature weights - * @param [in] N polynomial degree - */ - void lgNodesWeights(Eigen::VectorXd& nodes, Eigen::VectorXd& weights, const int N) { - // tolerance and max #iterations for Newton iteration - int nIterations = 10; - double tolerance = 1e-15; - - switch (N) { - case 0: - nodes[0] = 0.0; - weights[0] = 2.0; - break; - case 1: - nodes[0] = -std::sqrt(1.0 / 3.0); - weights[0] = 1; - nodes[1] = -nodes[0]; - weights[1] = weights[0]; - break; - default: - - double leg = 0.0; - double legDer = 0.0; - double delta = 0.0; - - for (int j = 0; j <= std::floor((N + 1) / 2) - 1; j++) - { - nodes[j] = -std::cos((2.0 * j + 1.0) / (2.0 * N + 2.0) * M_PI); - for (int k = 0; k <= nIterations; k++) - { - legendrePolynomialAndDerivative(leg, legDer, N + 1, nodes[j]); - delta = -leg / legDer; - nodes[j] = nodes[j] + delta; - if (std::abs(delta) <= tolerance * std::abs(nodes[j])) - break; - } - legendrePolynomialAndDerivative(leg, legDer, N + 1, nodes[j]); - nodes[N - j] = -nodes[j]; - weights[j] = 2.0 / ((1.0 - std::pow(nodes[j], 2.0)) * std::pow(legDer, 2.0)); - weights[N - j] = weights[j]; - } - - if (N % 2 == 0) - { - legendrePolynomialAndDerivative(leg, legDer, N + 1, 0.0); - nodes[N / 2] = 0.0; - weights[N / 2] = 2.0 / std::pow(legDer, 2.0); - } - } - } - - /** - * @brief evaluates a Lagrange polynomial built on input nodes at a set of points - * @detail can be used to establish quadrature rules - * @param [in] j index of Lagrange basis function - * @param [in] intNodes interpolation nodes the Lagrange basis is constructed with - * @param [in] evalNodes nodes the Lagrange basis is evaluated at - */ - Eigen::VectorXd evalLagrangeBasis(const int j, const Eigen::VectorXd intNodes, const Eigen::VectorXd evalNodes) { - - const int nIntNodes = intNodes.size(); - const int nEvalNodes = evalNodes.size(); - VectorXd evalEll = VectorXd::Zero(nEvalNodes); - - double nominator = 1.0; - double denominator = 1.0; - - for (int i = 0; i < nIntNodes; i++) - if (i != j) - denominator *= (intNodes[j] - intNodes[i]); - - for (int k = 0; k < nEvalNodes; k++) - { - for (int i = 0; i < nIntNodes; i++) - { - if (i != j) - { - if (std::abs(evalNodes[k] - intNodes[i]) < std::numeric_limits::epsilon()) - { - nominator = denominator; - break; - } - else - nominator *= (evalNodes[k] - intNodes[i]); - } - } - evalEll[k] = nominator / denominator; - nominator = 1.0; - } - - return evalEll; - } - - /** - * @brief calculates the Gauss quadrature mass matrix for LGL interpolation polynomial on LG points - * @detail exact integration of polynomials up to order 2N - 1 - * @param [in] LGLnodes Lagrange Gauss Lobatto nodes - * @param [in] nGaussNodes number of Gauss quadrature nodes - */ - Eigen::MatrixXd gaussQuadratureMMatrix(const Eigen::VectorXd LGLnodes, const int nLGNodes) { - - const int Ldegree = nLGNodes - 1; // Legendre polynomial degree - const int nLGLnodes = LGLnodes.size(); - - MatrixXd evalEll = MatrixXd::Zero(nLGNodes, nLGNodes); - MatrixXd massMatrix = MatrixXd::Zero(nLGNodes, nLGNodes); - - VectorXd LGnodes = VectorXd::Zero(nLGNodes); - VectorXd LGweigths = VectorXd::Zero(nLGNodes); - lgNodesWeights(LGnodes, LGweigths, Ldegree); - - for (int i = 0; i < nLGLnodes; i++) - evalEll.row(i) = evalLagrangeBasis(i, LGLnodes, LGnodes); - - VectorXd aux = VectorXd::Zero(nLGNodes); - - for (int i = 0; i < nLGLnodes; i++) - { - for (int j = 0; j < nLGLnodes; j++) - { - aux = evalEll.row(i).array() * evalEll.row(j).array(); - massMatrix(i, j) = (aux.array() * LGweigths.array()).sum(); - } - } - - return massMatrix; - } - - /** - * @brief computation of barycentric weights for fast polynomial evaluation - * @param [in] baryWeights vector to store barycentric weights. Must already be initialized with ones! - */ - void barycentricWeights(Eigen::VectorXd& baryWeights) { - for (unsigned int j = 1; j <= _polyDeg; j++) { - for (unsigned int k = 0; k <= j - 1; k++) { - baryWeights[k] = baryWeights[k] * (_nodes[k] - _nodes[j]) * 1.0; - baryWeights[j] = baryWeights[j] * (_nodes[j] - _nodes[k]) * 1.0; - } - } - for (unsigned int j = 0; j <= _polyDeg; j++) { - baryWeights[j] = 1 / baryWeights[j]; - } - } - - /** - * @brief computation of nodal (lagrange) polynomial derivative matrix - */ - void derivativeMatrix() { - Eigen::VectorXd baryWeights = Eigen::VectorXd::Ones(_polyDeg + 1u); - barycentricWeights(baryWeights); - for (unsigned int i = 0; i <= _polyDeg; i++) { - for (unsigned int j = 0; j <= _polyDeg; j++) { - if (i != j) { - _polyDerM(i, j) = baryWeights[j] / (baryWeights[i] * (_nodes[i] - _nodes[j])); - _polyDerM(i, i) += -_polyDerM(i, j); - } - } - } - } - - /** - * @brief factor to normalize legendre polynomials - */ - double orthonFactor(int _polyDeg) { - - double n = static_cast (_polyDeg); - // alpha = beta = 0 to get legendre polynomials as special case from jacobi polynomials. - double a = 0.0; - double b = 0.0; - return std::sqrt(((2.0 * n + a + b + 1.0) * std::tgamma(n + 1.0) * std::tgamma(n + a + b + 1.0)) - / (std::pow(2.0, a + b + 1.0) * std::tgamma(n + a + 1.0) * std::tgamma(n + b + 1.0))); - } + * Functions to calculate Jacobian blocks + * =================================================================================== */ - /** - * @brief calculates the Vandermonde matrix of the normalized legendre polynomials - */ - Eigen::MatrixXd getVandermonde_LEGENDRE() { - - Eigen::MatrixXd V(_nodes.size(), _nodes.size()); - - double alpha = 0.0; - double beta = 0.0; - - // degree 0 - V.block(0, 0, _nNodes, 1) = Eigen::VectorXd::Ones(_nNodes) * orthonFactor(0); - // degree 1 - for (int node = 0; node < static_cast(_nNodes); node++) { - V(node, 1) = _nodes[node] * orthonFactor(1); - } - - for (int deg = 2; deg <= static_cast(_polyDeg); deg++) { - - for (int node = 0; node < static_cast(_nNodes); node++) { - - double orthn_1 = orthonFactor(deg) / orthonFactor(deg - 1); - double orthn_2 = orthonFactor(deg) / orthonFactor(deg - 2); - - double fac_1 = ((2.0 * deg - 1.0) * 2.0 * deg * (2.0 * deg - 2.0) * _nodes[node]) / (2.0 * deg * deg * (2.0 * deg - 2.0)); - double fac_2 = (2.0 * (deg - 1.0) * (deg - 1.0) * 2.0 * deg) / (2.0 * deg * deg * (2.0 * deg - 2.0)); - - V(node, deg) = orthn_1 * fac_1 * V(node, deg - 1) - orthn_2 * fac_2 * V(node, deg - 2); - - } - - } - - return V; - } - /** - * @brief calculates mass matrix via transformation to orthonormal Legendre (modal) basis - * @detail exact integration for integrals of the form \int_E \ell_i \ell_j d\xi - */ - void invMMatrix() { - _invMM = (getVandermonde_LEGENDRE() * (getVandermonde_LEGENDRE().transpose())); - } /** * @brief calculates the convection part of the DG jacobian */ @@ -570,10 +233,6 @@ namespace cadet return -convBlock; // *-1 for residual } - /* =================================================================================== - * Functions to calculate Jacobian blocks - * =================================================================================== */ - /** * @brief calculates the DG Jacobian auxiliary block * @param [in] exInt true if exact integration DG scheme diff --git a/src/libcadet/model/parts/DGToolbox.cpp b/src/libcadet/model/parts/DGToolbox.cpp new file mode 100644 index 000000000..fe9f03140 --- /dev/null +++ b/src/libcadet/model/parts/DGToolbox.cpp @@ -0,0 +1,393 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2022: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +/** + * @file + * Defines the convection dispersion transport operator according to the discontinuous Galerkin discretization. + */ + +#include "model/parts/DGToolbox.hpp" + +using namespace Eigen; + +namespace cadet +{ + +namespace model +{ + +namespace parts +{ + +namespace dgtoolbox +{ + +/** + * @brief computes the Legendre polynomial L_N and q = L_N+1 - L_N-2 and q' at point x + * @param [in] polyDeg polynomial degree + * @param [in] x evaluation point + * @param [in] L <- L(x) + * @param [in] q <- q(x) = L_N+1 (x) - L_N-2(x) + * @param [in] qder <- q'(x) = [L_N+1 (x) - L_N-2(x)]' + */ +void qAndL(const unsigned int polyDeg, const double x, double& L, double& q, double& qder) { + // auxiliary variables (Legendre polynomials) + double L_2 = 1.0; + double L_1 = x; + double Lder_2 = 0.0; + double Lder_1 = 1.0; + double Lder = 0.0; + for (double k = 2; k <= polyDeg; k++) { // note that this function is only called for polyDeg >= 2. + L = ((2 * k - 1) * x * L_1 - (k - 1) * L_2) / k; + Lder = Lder_2 + (2 * k - 1) * L_1; + L_2 = L_1; + L_1 = L; + Lder_2 = Lder_1; + Lder_1 = Lder; + } + q = ((2.0 * polyDeg + 1) * x * L - polyDeg * L_2) / (polyDeg + 1.0) - L_2; + qder = Lder_1 + (2.0 * polyDeg + 1) * L_1 - Lder_2; +} +/** + * @brief computes the Legendre-Gauss-Lobatto nodes and (inverse) quadrature weights + * @param [in] polyDeg polynomial degree + * @param [in, out] nodes Legendre Gauss Lobatto nodes + * @param [in, out] invWeights Legendre Gauss quadrature weights + * @param [in] invertWeights specifies if weights should be inverted + */ +void lglNodesWeights(const unsigned int polyDeg, VectorXd& nodes, VectorXd& invWeights, bool invertWeights) { + + const double pi = 3.1415926535897932384626434; + + // tolerance and max #iterations for Newton iteration + int nIterations = 10; + double tolerance = 1e-15; + // Legendre polynomial and derivative + double L = 0; + double q = 0; + double qder = 0; + switch (polyDeg) { + case 0: + throw std::invalid_argument("Polynomial degree must be at least 1 !"); + break; + case 1: + nodes[0] = -1; + invWeights[0] = 1; + nodes[1] = 1; + invWeights[1] = 1; + break; + default: + nodes[0] = -1; + nodes[polyDeg] = 1; + invWeights[0] = 2.0 / (polyDeg * (polyDeg + 1.0)); + invWeights[polyDeg] = invWeights[0]; + // use symmetrie, only compute half of points and weights + for (unsigned int j = 1; j <= floor((polyDeg + 1) / 2) - 1; j++) { + // first guess for Newton iteration + nodes[j] = -cos(pi * (j + 0.25) / polyDeg - 3 / (8.0 * polyDeg * pi * (j + 0.25))); + // Newton iteration to find roots of Legendre Polynomial + for (unsigned int k = 0; k <= nIterations; k++) { + qAndL(polyDeg, nodes[j], L, q, qder); + nodes[j] = nodes[j] - q / qder; + if (abs(q / qder) <= tolerance * abs(nodes[j])) { + break; + } + } + // calculate weights + qAndL(polyDeg, nodes[j], L, q, qder); + invWeights[j] = 2.0 / (polyDeg * (polyDeg + 1.0) * pow(L, 2.0)); + nodes[polyDeg - j] = -nodes[j]; // copy to second half of points and weights + invWeights[polyDeg - j] = invWeights[j]; + } + } + if (polyDeg % 2 == 0) { // for even polyDeg we have an odd number of points which include 0.0 + qAndL(polyDeg, 0.0, L, q, qder); + nodes[polyDeg / 2] = 0; + invWeights[polyDeg / 2] = 2.0 / (polyDeg * (polyDeg + 1.0) * pow(L, 2.0)); + } + // inverse the weights + invWeights = invWeights.cwiseInverse(); +} +/** + * @brief computes the Legendre polynomial and its derivative + * @param [in] polyDeg polynomial degree + * @param [in, out] leg Legendre polynomial + * @param [in, out] legDer Legendre polynomial derivative + * @param [in] x evaluation point + */ +void legendrePolynomialAndDerivative(const int polyDeg, double& leg, double& legDer, const double x) { + + switch (polyDeg) { + case 0: + leg = 1.0; + legDer = 0.0; + break; + case 1: + leg = x; + legDer = 1.0; + break; + default: + double leg_2 = 1.0; + double leg_1 = x; + double legDer_2 = 0.0; + double legDer_1 = 1.0; + + for (int k = 2; k <= polyDeg; k++) { + leg = (2.0 * k - 1.0) / k * x * leg_1 - (k - 1.0) / k * leg_2; + legDer = legDer_2 + (2.0 * k - 1.0) * leg_1; + leg_2 = leg_1; + leg_1 = leg; + legDer_2 = legDer_1; + legDer_1 = legDer; + } + } +} +/** + * @brief computes the Legendre-Gauss nodes and quadrature weights + * @param [in] polyDeg polynomial degree + * @param [in, out] nodes Legendre Gauss nodes + * @param [in, out] weights Legendre Gauss quadrature weights + * @param [in] invertWeights specifies if weights should be inverted + */ +void lgNodesWeights(const unsigned int polyDeg, VectorXd& nodes, VectorXd& weights, bool invertWeights = true) { + + const double pi = 3.1415926535897932384626434; + + // tolerance and max #iterations for Newton iteration + int nIterations = 10; + double tolerance = 1e-15; + + switch (polyDeg) { + case 0: + nodes[0] = 0.0; + weights[0] = 2.0; + break; + case 1: + nodes[0] = -std::sqrt(1.0 / 3.0); + weights[0] = 1; + nodes[1] = -nodes[0]; + weights[1] = weights[0]; + break; + default: + + double leg = 0.0; + double legDer = 0.0; + double delta = 0.0; + + for (int j = 0; j <= std::floor((polyDeg + 1) / 2) - 1; j++) + { + nodes[j] = -std::cos((2.0 * j + 1.0) / (2.0 * polyDeg + 2.0) * pi); + for (int k = 0; k <= nIterations; k++) + { + legendrePolynomialAndDerivative(polyDeg + 1, leg, legDer, nodes[j]); + delta = -leg / legDer; + nodes[j] = nodes[j] + delta; + if (std::abs(delta) <= tolerance * std::abs(nodes[j])) + break; + } + legendrePolynomialAndDerivative(polyDeg + 1, leg, legDer, nodes[j]); + nodes[polyDeg - j] = -nodes[j]; + weights[j] = 2.0 / ((1.0 - std::pow(nodes[j], 2.0)) * std::pow(legDer, 2.0)); + weights[polyDeg - j] = weights[j]; + } + + if (polyDeg % 2 == 0) + { + legendrePolynomialAndDerivative(polyDeg + 1, leg, legDer, 0.0); + nodes[polyDeg / 2] = 0.0; + weights[polyDeg / 2] = 2.0 / std::pow(legDer, 2.0); + } + } +} +/** + * @brief evaluates a Lagrange polynomial built on input nodes at a set of points + * @detail can be used to establish quadrature rules + * @param [in] j index of Lagrange basis function + * @param [in] intNodes interpolation nodes the Lagrange basis is constructed with + * @param [in] evalNodes nodes the Lagrange basis is evaluated at + */ +VectorXd evalLagrangeBasis(const int j, const VectorXd intNodes, const VectorXd evalNodes) { + + const int nIntNodes = intNodes.size(); + const int nEvalNodes = evalNodes.size(); + VectorXd evalEll = VectorXd::Zero(nEvalNodes); + + double nominator = 1.0; + double denominator = 1.0; + + for (int i = 0; i < nIntNodes; i++) + if (i != j) + denominator *= (intNodes[j] - intNodes[i]); + + for (int k = 0; k < nEvalNodes; k++) + { + for (int i = 0; i < nIntNodes; i++) + { + if (i != j) + { + if (std::abs(evalNodes[k] - intNodes[i]) < std::numeric_limits::epsilon()) + { + nominator = denominator; + break; + } + else + nominator *= (evalNodes[k] - intNodes[i]); + } + } + evalEll[k] = nominator / denominator; + nominator = 1.0; + } + + return evalEll; +} +/** + * @brief calculates the Gauss quadrature mass matrix for LGL interpolation polynomial on LG points + * @detail exact integration of polynomials up to order 2N - 1 + * @param [in] LGLnodes Legendre Gauss Lobatto nodes + * @param [in] nLGNodes number of Gauss quadrature nodes + */ +MatrixXd gaussQuadratureMMatrix(const VectorXd LGLnodes, const int nLGNodes) { + + const int Ldegree = nLGNodes - 1; // Legendre polynomial degree + const int nLGLnodes = LGLnodes.size(); + + MatrixXd evalEll = MatrixXd::Zero(nLGNodes, nLGNodes); + MatrixXd massMatrix = MatrixXd::Zero(nLGNodes, nLGNodes); + + VectorXd LGnodes = VectorXd::Zero(nLGNodes); + VectorXd LGweigths = VectorXd::Zero(nLGNodes); + lgNodesWeights(Ldegree, LGnodes, LGweigths, false); + + for (int i = 0; i < nLGLnodes; i++) + evalEll.row(i) = evalLagrangeBasis(i, LGLnodes, LGnodes); + + VectorXd aux = VectorXd::Zero(nLGNodes); + + for (int i = 0; i < nLGLnodes; i++) + { + for (int j = 0; j < nLGLnodes; j++) + { + aux = evalEll.row(i).array() * evalEll.row(j).array(); + massMatrix(i, j) = (aux.array() * LGweigths.array()).sum(); + } + } + + return massMatrix; +} +/** + * @brief computation of barycentric weights for fast polynomial evaluation + * @param [in] polyDeg polynomial degree + * @param [in, out] baryWeights vector to store barycentric weights. Must already be initialized with ones! + */ +VectorXd barycentricWeights(const unsigned int polyDeg, const VectorXd nodes) { + + VectorXd baryWeights = VectorXd::Ones(polyDeg + 1u); + + for (unsigned int j = 1; j <= polyDeg; j++) { + for (unsigned int k = 0; k <= j - 1; k++) { + baryWeights[k] = baryWeights[k] * (nodes[k] - nodes[j]) * 1.0; + baryWeights[j] = baryWeights[j] * (nodes[j] - nodes[k]) * 1.0; + } + } + for (unsigned int j = 0; j <= polyDeg; j++) { + baryWeights[j] = 1 / baryWeights[j]; + } + + return baryWeights; +} +/** + * @brief computation of nodal (lagrange) polynomial derivative matrix + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +MatrixXd derivativeMatrix(const unsigned int polyDeg, const VectorXd nodes) { + + MatrixXd polyDerM = MatrixXd::Zero(polyDeg + 1u, polyDeg + 1u); + VectorXd baryWeights = barycentricWeights(polyDeg, nodes); + + for (unsigned int i = 0; i <= polyDeg; i++) { + for (unsigned int j = 0; j <= polyDeg; j++) { + if (i != j) { + polyDerM(i, j) = baryWeights[j] / (baryWeights[i] * (nodes[i] - nodes[j])); + polyDerM(i, i) += -polyDerM(i, j); + } + } + } + + return polyDerM; +} +/** + * @brief factor to normalize Jacobi polynomials + */ +double orthonFactor(const int polyDeg, double a, double b) { + + double n = static_cast (polyDeg); + return std::sqrt(((2.0 * n + a + b + 1.0) * std::tgamma(n + 1.0) * std::tgamma(n + a + b + 1.0)) + / (std::pow(2.0, a + b + 1.0) * std::tgamma(n + a + 1.0) * std::tgamma(n + b + 1.0))); +} +/** + * @brief calculates the Vandermonde matrix of the normalized Jacobi polynomials + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + * @param [in] a Jacobi polynomial parameter + * @param [in] b Jacobi polynomial parameter + */ +MatrixXd getVandermonde_JACOBI(const unsigned int polyDeg, const VectorXd nodes, const double a, const double b) { + + const unsigned int nNodes = polyDeg + 1u; + MatrixXd V(nNodes, nNodes); + + // degree 0 + V.block(0, 0, nNodes, 1) = VectorXd::Ones(nNodes) * orthonFactor(0, a, b); + // degree 1 + for (int node = 0; node < static_cast(nNodes); node++) { + V(node, 1) = ((nodes[node] - 1.0) / 2.0 * (a + b + 2.0) + (a + 1.0)) * orthonFactor(1, a, b); + } + + for (int deg = 2; deg <= static_cast(nNodes - 1); deg++) { + + for (int node = 0; node < static_cast(nNodes); node++) { + + double orthn_1 = orthonFactor(deg, a, b) / orthonFactor(deg - 1, a, b); + double orthn_2 = orthonFactor(deg, a, b) / orthonFactor(deg - 2, a, b); + + // recurrence relation + V(node, deg) = orthn_1 * ((2.0 * deg + a + b - 1.0) * ((2.0 * deg + a + b) * (2.0 * deg + a + b - 2.0) * nodes[node] + a * a - b * b) * V(node, deg - 1)); + V(node, deg) -= orthn_2 * (2.0 * (deg + a - 1.0) * (deg + b - 1.0) * (2.0 * deg + a + b) * V(node, deg - 2)); + V(node, deg) /= 2.0 * deg * (deg + a + b) * (2.0 * deg + a + b - 2.0); + } + } + + return V; +} +/** + * @brief calculates the Vandermonde matrix of the normalized Legendre polynomials + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +MatrixXd getVandermonde_LEGENDRE(const unsigned int polyDeg, const VectorXd nodes) { + return getVandermonde_JACOBI(polyDeg, nodes, 0.0, 0.0); +} +/** + * @brief calculates mass matrix via transformation to orthonormal Jacobi (modal) basis + * @detail exact integration for integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +Eigen::MatrixXd invMMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes, const double alpha, const double beta) { + return (getVandermonde_JACOBI(polyDeg, nodes, alpha, beta) * (getVandermonde_JACOBI(polyDeg, nodes, alpha, beta).transpose())); +} + +} // namespace dgtoolbox +} // namespace parts +} // namespace model +} // namespace cadet \ No newline at end of file diff --git a/src/libcadet/model/parts/DGToolbox.hpp b/src/libcadet/model/parts/DGToolbox.hpp new file mode 100644 index 000000000..0dd9c2d56 --- /dev/null +++ b/src/libcadet/model/parts/DGToolbox.hpp @@ -0,0 +1,84 @@ +// ============================================================================= +// CADET +// +// Copyright © 2008-2022: The CADET Authors +// Please see the AUTHORS and CONTRIBUTORS file. +// +// All rights reserved. This program and the accompanying materials +// are made available under the terms of the GNU Public License v3.0 (or, at +// your option, any later version) which accompanies this distribution, and +// is available at http://www.gnu.org/licenses/gpl.html +// ============================================================================= + +/** + * @file + * Defines a discrete calculus toolbox for nodal spectral elements. + * @details See @cite Kopriva2009 + */ + +#ifndef LIBCADET_DGTOOLBOX_HPP_ +#define LIBCADET_DGTOOLBOX_HPP_ + +#include + +namespace cadet +{ + +namespace model +{ + +namespace parts +{ + +namespace dgtoolbox +{ + +/** + * @brief computes the Legendre-Gauss-Lobatto nodes and (inverse) quadrature weights + * @param [in] polyDeg polynomial degree + * @param [in, out] nodes Legendre Gauss Lobatto nodes + * @param [in, out] invWeights Legendre Gauss quadrature weights + * @param [in] invertWeights specifies if weights should be inverted + */ +void lglNodesWeights(const unsigned int polyDeg, Eigen::VectorXd& nodes, Eigen::VectorXd& invWeights, bool invertWeights = true); +/** + * @brief computes the Legendre polynomial and its derivative + * @param [in] polyDeg polynomial degree + * @param [in, out] leg Legendre polynomial + * @param [in, out] legDer Legendre polynomial derivative + * @param [in] x evaluation point + */ +void legendrePolynomialAndDerivative(const int polyDeg, double& leg, double& legDer, const double x); +/** + * @brief calculates the Gauss quadrature mass matrix for LGL interpolation polynomial on LG points + * @detail exact integration of polynomials up to order 2N - 1 + * @param [in] LGLnodes Legendre Gauss Lobatto nodes + * @param [in] nLGNodes number of Gauss quadrature nodes + */ +Eigen::MatrixXd gaussQuadratureMMatrix(const Eigen::VectorXd LGLnodes, const int nLGNodes); +/** + * @brief computation of barycentric weights for fast polynomial evaluation + * @param [in] polyDeg polynomial degree + * @param [in, out] baryWeights vector to store barycentric weights. Must already be initialized with ones! + */ +Eigen::VectorXd barycentricWeights(const unsigned int polyDeg, const Eigen::VectorXd nodes); +/** + * @brief computation of nodal (lagrange) polynomial derivative matrix + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +Eigen::MatrixXd derivativeMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes); +/** + * @brief calculates mass matrix via transformation to orthonormal Jacobi (modal) basis + * @detail exact integration for integrals of the form \int_E \ell_i(\xi) \ell_j(\xi) (1 - \xi)^\alpha (1 + \xi)^\beta d\xi + * @param [in] polyDeg polynomial degree + * @param [in] nodes polynomial interpolation nodes + */ +Eigen::MatrixXd invMMatrix(const unsigned int polyDeg, const Eigen::VectorXd nodes, const double alpha = 0.0, const double beta = 0.0); + +} // namespace dgtoolbox +} // namespace parts +} // namespace model +} // namespace cadet + +#endif // LIBCADET_DGTOOLBOX_HPP_ \ No newline at end of file