From 3080dfe79048fef85a00e8623ae49d18c749ac6c Mon Sep 17 00:00:00 2001
From: Florian Atteneder <florian.atteneder@uni-jena.de>
Date: Tue, 13 Aug 2024 12:11:11 +0200
Subject: [PATCH] wip
---
src/mesh.jl | 18 ++++++++++--------
src/spectralelement.jl | 22 +++++++++++-----------
2 files changed, 21 insertions(+), 19 deletions(-)
diff --git a/src/mesh.jl b/src/mesh.jl
index 75ed9541..e527633a 100644
--- a/src/mesh.jl
+++ b/src/mesh.jl
@@ -30,7 +30,8 @@ const Mesh2d{Element} = Mesh{Element,2,4}
#######################################################################
-function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(true,), scheme="DG",
+function Mesh1d(; N=5, K=10, range=[-1.0,1.0], periodic=(true,),
+ basis="lgl_lumped", dg_kind="nodal_legendre", scheme="DG",
filter_parameters=parameters(:Mesh)["filter_parameters"])
@toggled_assert N > 0
@@ -45,7 +46,7 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
# TODO Remove this after updating parameter names in SpectralElement
quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
ηc, sh, α = filter_parameters
- SpectralElement(N, Symbol(quadrature_method); ηc, sh, α)
+ SpectralElement(N, Symbol(quadrature_method), Symbol(dg_kind); ηc, sh, α)
elseif scheme == "FV"
FVElement()
else
@@ -98,9 +99,8 @@ function Mesh1d(; N=5, K=10, range=[-1.0,1.0], basis="lgl_lumped", periodic=(tru
end
-function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
- basis="lgl_lumped", periodic=(true,true),
- scheme="DG",
+function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0], periodic=(true,true),
+ basis="lgl_lumped", dg_kind="nodal_legendre", scheme="DG",
filter_parameters=parameters(:Mesh)["filter_parameters"])
@toggled_assert length(N) == 2
@@ -124,8 +124,8 @@ function Mesh2d(; N=[5,5], K=[4,4], range=[-1.0,1.0, -1.0,1.0],
# TODO Remove this after updating parameter names in SpectralElement
quadrature_method = Dict("lgl" => :LGL, "glgl" => :GLGL, "lgl_lumped" => :LGL_lumped)[basis]
ηc, sh, α = filter_parameters
- element_x = SpectralElement(Nx, Symbol(quadrature_method); ηc, sh, α)
- element_y = SpectralElement(Ny, Symbol(quadrature_method); ηc, sh, α)
+ element_x = SpectralElement(Nx, Symbol(quadrature_method), Symbol(dg_kind); ηc, sh, α)
+ element_y = SpectralElement(Ny, Symbol(quadrature_method), Symbol(dg_kind); ηc, sh, α)
elseif scheme == "FV"
TODO()
element_x = FVElement()
@@ -235,12 +235,13 @@ function Mesh(; N=parameters(:Mesh)["n"],
K=parameters(:Mesh)["k"],
range=parameters(:Mesh)["range"],
basis=parameters(:Mesh)["basis"],
+ dg_kind=parameters(:Mesh)["kind"],
periodic=parameters(:Mesh)["periodic"],
scheme=parameters(:Mesh)["scheme"],
filter_parameters=parameters(:Mesh)["filter_parameters"])
dims = Int(length(range)/2)
MType = dims == 1 ? Mesh1d : Mesh2d
- return MType(; N, K, range, basis, periodic, scheme, filter_parameters)
+ return MType(; N, K, range, basis, dg_kind, periodic, scheme, filter_parameters)
end
@@ -269,6 +270,7 @@ extends(m::Mesh,cell::Cell) = m.boxes[cell.index].extends
function min_grid_spacing(m::Mesh{SpectralElement})
+ # TODO Use nodal_points(m)
@unpack z = m.element
ws = widths(m.boxes[1])
minw = minimum(ws)
diff --git a/src/spectralelement.jl b/src/spectralelement.jl
index 30000da3..3a259856 100644
--- a/src/spectralelement.jl
+++ b/src/spectralelement.jl
@@ -36,35 +36,33 @@ struct SpectralElement
if quadr === :LGL
qaudr_z, qaudr_w = LGL.rule(Npts); qaudr_v = -1.0 # unused here
purge_zeros!(qaudr_z)
- z = quadr_z
- D = first_derivative_matrix(z)
S = LGL.stiffness_matrix(Npts)
invM = LGL.inverse_mass_matrix_exact(Npts)
elseif quadr === :LGL_lumped
quadr_z, quadr_w = LGL.rule(Npts); quadr_v = -1.0 # unused here
purge_zeros!(quadr_z)
- z = quadr_z
- D = first_derivative_matrix(z)
S = LGL.stiffness_matrix(Npts)
invM = LGL.inverse_mass_matrix_lumped(Npts)
elseif quadr === :GLGL
quadr_z, quadr_w, quadr_v = GLGL.rule(Npts)
purge_zeros!(quadr_z)
- z = quadr_z
- D = first_derivative_matrix(z)
- D2 = second_derivative_matrix(z)
- M = GLGL.mass_matrix(z, w, v, D)
- S = GLGL.stiffness_matrix(z, w, v, D, D2)
+ quadr_D = first_derivative_matrix(quadr_z)
+ quadr_D2 = second_derivative_matrix(quadr_z)
+ M = GLGL.mass_matrix(quadr_z, quadr_w, quadr_v, quadr_D)
+ S = GLGL.stiffness_matrix(quadr_z, quadr_w, quadr_v, quadr_D, quadr_D2)
invM = inv(M)
else
error("Invalid quadrature rule specified for kind $kind: $quadr. Use one of :LGL_lumped, :LGL, :GLGL.")
end
+ z = quadr_z
+ D = first_derivative_matrix(z)
+ quadr_D = D
l_lhs = zeros(Float64, N+1)
l_rhs = zeros(Float64, N+1)
l_lhs[1] = 1.0
l_rhs[end] = 1.0
- quadr_D = D
V = vandermonde_matrix_legendre(z)
+ # In preparation for Bernstein basis
# elseif kind === :nodal_bernstein
# if quadr === :LGL || quadr === :LGL_lumped
# qaudr_z, qaudr_w = LGL.rule(Npts)
@@ -82,9 +80,11 @@ struct SpectralElement
# S = Bernstein.stiffness_matrix(Npts)
# invM = inv(M)
# D = Bernstein.first_derivative_matrix(Npts)
- # D2 = Bernstein.second_derivative_matrix(Npts)
+ # TODO("quadr_D")
+ # TODO("quadr_w")
# V = Bernstein.vandermonde_matrix(z)
# I = stack(Bernstein.bernstein_polynomial.(quadr_z, i-1, N) for i in 1:Npts)
+ # quadr_w = transpose(I) * quadr_w
# quadr_D = D * I
# l_lhs = zeros(Float64, N+1)
# l_rhs = zeros(Float64, N+1)
--
GitLab