Skip to content
Snippets Groups Projects
Commit 950657fd authored by Florian Atteneder's avatar Florian Atteneder
Browse files

wip: add hybrid mesh

parent 839a6a95
No related branches found
No related tags found
No related merge requests found
...@@ -43,6 +43,7 @@ module dg1d ...@@ -43,6 +43,7 @@ module dg1d
include("spectralelement.jl") include("spectralelement.jl")
include("subvolumeelement.jl") include("subvolumeelement.jl")
include("mesh.jl") include("mesh.jl")
include("hybridmesh.jl")
# custom evolution # custom evolution
include("evolve.jl") include("evolve.jl")
......
export compute_rhs_weak_form!, export compute_rhs_weak_form!,
compute_rhs_strong_form! compute_rhs_strong_form!,
compute_rhs_fv_form!
""" """
...@@ -55,3 +56,36 @@ function compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::AbstractMesh) ...@@ -55,3 +56,36 @@ function compute_rhs_strong_form!(rhs, f, s, nf_lhs, nf_rhs, mesh::AbstractMesh)
rhs .+= s rhs .+= s
return return
end end
"""
compute_rhs_weak_form!(rhs, sv_flag, f, mesh::HybridMesh)
compute_rhs_weak_form!(rhs, sv_flag, f, s, mesh::HybridMesh)
"""
function compute_rhs_fv_form!(rhs, sv_flag, f, hybridmesh::HybridMesh)
@unpack mesh, sv_element = hybridmesh
@unpack invjac = mesh
shape = layout(mesh)
mat_rhs = reshape(view(rhs,:), shape)
mat_f = reshape(view(f,:), shape)
Npts, K = shape
for (k#=cell nr=#, flag) in enumerate(sv_flag)
flag == 0 && continue
for i = 1:Npts-1
mat_rhs[i,k] = - (mat_f[i+1,k] - mat_f[i,k]) / sv_element.w[i] * invjac[k]
end
end
return
end
function compute_rhs_fv_form!(rhs, sv_flag, f, s, hybridmesh::HybridMesh)
compute_rhs_weak_form!(rhs, tci_flag, f, hybridmesh)
shape = layout(mesh)
mat_rhs = reshape(view(rhs,:), shape)
mat_s = reshape(view(s, :), shape)
for (k#=cell nr=#, flag) in enumerate(sv_flag)
flag == 0 && continue
mat_rhs[:,k] .+= mat_s[:,k]
end
rhs .+= s
return
end
struct HybridMesh <: AbstractMesh
mesh::Mesh
sv_element::SubVolumeElement
function HybridMesh(mesh::Mesh)
sv_element = SubVolumeElement(mesh.element)
return new(mesh, sv_element)
end
end
function Base.getproperty(hmesh::HybridMesh, p::Symbol)
if p in propertynames(hmesh.mesh)
return getproperty(hmesh.mesh, p)
else
getfield(hmesh.mesh, p)
end
end
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment