Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

matrix assembly update #161

Merged
merged 4 commits into from
Dec 23, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 5 additions & 5 deletions GalerkinToolkitExamples/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,8 @@ authors = ["Francesc Verdugo <[email protected]>"]
version = "0.1.0"

[deps]
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
GalerkinToolkit = "5e3ba9c4-dd81-444d-b69a-0e7bd7bf60a4"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
Expand All @@ -20,14 +22,12 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

[compat]
PartitionedSolvers = "0.3"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]

[compat]
PartitionedSolvers = "0.3"

31 changes: 15 additions & 16 deletions GalerkinToolkitExamples/src/poisson.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,8 @@ using SparseArrays
using PartitionedArrays
using TimerOutputs
using WriteVTK
using AbstractTrees


function main(params_in)
params_default = default_params()
Expand All @@ -26,18 +28,15 @@ function main(params_in)
end
end

gradient(u) = x->ForwardDiff.gradient(u,x)
jacobian(u) = x->ForwardDiff.jacobian(u,x)
laplacian(u,x) = tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
laplacian(u) = x->laplacian(u,x)
# TODO hide this
laplacian(u::GT.AbstractQuantity,x::GT.AbstractQuantity) = GT.call(laplacian,u,x)

Δ(u) = GT.call(laplacian,u)
∇(u) = GT.call(gradient,u)
Δ(u,x) = Δ(u)(x)
∇(u,x) = ∇(u)(x)
const ∇ = ForwardDiff.gradient
const Δ = laplacian

mean(u,x) = 0.5*(u(x)[+]+u(x)[-])
jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x)
mean(u) = 0.5*(u[1]+u[2])
jump(u,n) = u[2]*n[2] + u[1]*n[1]

function main_automatic(params)
timer = params[:timer]
Expand All @@ -54,7 +53,7 @@ function main_automatic(params)

u = GT.analytical_field(params[:u],Ω)
f(x) = -Δ(u,x)
n_Γn = GT.unit_normal(Γn,Ω)
n_Γn = GT.unit_normal(mesh,D-1)
g(x) = n_Γn(x)⋅∇(u,x)

interpolation_degree = params[:interpolation_degree]
Expand All @@ -68,7 +67,7 @@ function main_automatic(params)
GT.label_interior_faces!(mesh;physical_name="__INTERIOR_FACES__")
Λ = GT.skeleton(mesh;physical_names=["__INTERIOR_FACES__"])
dΛ = GT.measure(Λ,integration_degree)
n_Λ = GT.unit_normal(Λ,Ω)
n_Λ = GT.unit_normal(mesh,D-1)
h_Λ = GT.face_diameter_field(Λ)
else
conformity = :default
Expand All @@ -79,7 +78,7 @@ function main_automatic(params)
uhd = GT.dirichlet_field(Float64,V)
GT.interpolate_dirichlet!(u,uhd)
else
n_Γd = GT.unit_normal(Γd,Ω)
n_Γd = GT.unit_normal(mesh,D-1)
h_Γd = GT.face_diameter_field(Γd)
dΓd = GT.measure(Γd,integration_degree)
V = GT.lagrange_space(Ω,interpolation_degree;conformity)
Expand All @@ -98,7 +97,7 @@ function main_automatic(params)
end
if params[:discretization_method] === :interior_penalty
r += ∫( x->
(γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-jump(v,n_Λ,x)⋅mean(∇(u),x)-mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ)
(γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-jump(v(x),n_Λ(x))⋅mean(∇(u,x))-mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ)
end
r
end
Expand Down Expand Up @@ -202,7 +201,7 @@ function main_hand_written(params)
dΩ = GT.measure(Ω,integration_degree)

u = params[:u]
f(x) = -laplacian(u,x)
f(x) = -Δ(u,x)

# Set interpolation space -> Lagrange polynomial basis function, the element type and geometry
V = GT.lagrange_space(Ω,interpolation_degree;dirichlet_boundary=Γd)
Expand Down Expand Up @@ -564,8 +563,8 @@ function default_params()
params[:dirichlet_tags] = ["boundary"]
params[:neumann_tags] = String[]
params[:dirichlet_method] = :strong
params[:integration_degree] = 1
params[:interpolation_degree] = 2*params[:integration_degree]
params[:interpolation_degree] = 1
params[:integration_degree] = 2*params[:interpolation_degree]
params[:solver] = PS.LinearAlgebra_lu
params[:example_path] = joinpath(mkpath(joinpath(@__DIR__,"..","output")),"example_001")
params[:export_vtu] = true
Expand Down
4 changes: 2 additions & 2 deletions GalerkinToolkitExamples/test/poisson_tests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,7 @@ for discretization_method in (:continuous_galerkin,)
params[:dirichlet_tags] = ["1-face-1","1-face-2","1-face-3","1-face-4"]
params[:discretization_method] = discretization_method
params[:dirichlet_method] = dirichlet_method
params[:integration_degree] = 2
params[:interpolation_degree] = 2
results = Poisson.main(params)
@test results[:error_l2_norm] < tol
@test results[:error_h1_norm] < tol
Expand All @@ -35,7 +35,7 @@ for discretization_method in (:interior_penalty,:continuous_galerkin)
params[:neumann_tags] = ["1-face-2"]
params[:discretization_method] = discretization_method
params[:dirichlet_method] = dirichlet_method
params[:integration_degree] = 2
params[:interpolation_degree] = 2
results = Poisson.main(params)
@test results[:error_l2_norm] < tol
@test results[:error_h1_norm] < tol
Expand Down
35 changes: 18 additions & 17 deletions docs/src/examples/methods.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,21 +11,22 @@ import FileIO # hide
domain = (0,1,0,1)
cells = (10,10)
mesh = GT.cartesian_mesh(domain,cells)
D = GT.num_dims(mesh)
GT.label_boundary_faces!(mesh;physical_name="Γd")
Ω = GT.interior(mesh)
Γd = GT.boundary(mesh;physical_names=["Γd"])
n_Γd = GT.unit_normal(Γd,Ω)
n_Γd = GT.unit_normal(mesh,D-1)
u = GT.analytical_field(x->sum(x),Ω)
#TODO these two definitions should be enough
#∇ = ForwardDiff.gradient
#Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x))
gradient(u) = x->ForwardDiff.gradient(u,x)
jacobian(u) = x->ForwardDiff.jacobian(u,x)
laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
Δ(u) = GT.call(laplacian,u)
∇(u) = GT.call(gradient,u)
Δ(u,x) = Δ(u)(x)
∇(u,x) = ∇(u)(x)
∇ = ForwardDiff.gradient
Δ = (u,x) -> tr(ForwardDiff.jacobian(y->∇(u,y),x))
#gradient(u) = x->ForwardDiff.gradient(u,x)
#jacobian(u) = x->ForwardDiff.jacobian(u,x)
#laplacian(u) = x-> tr(ForwardDiff.jacobian(y->ForwardDiff.gradient(u,y),x))
#Δ(u) = GT.call(laplacian,u)
#∇(u) = GT.call(gradient,u)
#Δ(u,x) = Δ(u)(x)
#∇(u,x) = ∇(u)(x)
tol = 1e-10
order = 2
γ = order*(order+1)/10
Expand Down Expand Up @@ -54,7 +55,7 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ)))
# ## Nitche Method

V = GT.lagrange_space(Ω,order)
n_Γd = GT.unit_normal(Γd,Ω)
n_Γd = GT.unit_normal(mesh,D-1)
h_Γd = GT.face_diameter_field(Γd)
dΓd = GT.measure(Γd,2*order)
a_Γd = (u,v) -> GT.∫(
Expand Down Expand Up @@ -90,18 +91,18 @@ el2 = sqrt(sum(GT.∫( x->abs2(eh(x)), dΩ)))

# ## Interior Penalty

mean(u,x) = 0.5*(u(x)[+]+u(x)[-])
jump(u,n,x) = u(x)[+]*n[+](x) + u(x)[-]*n[-](x)
mean(u) = 0.5*(u[1]+u[2])
jump(u,n) = u[2]*n[2] + u[1]*n[1]
GT.label_interior_faces!(mesh;physical_name="Λ")
Λ = GT.skeleton(mesh;physical_names=["Λ"])
dΛ = GT.measure(Λ,2*order)
n_Λ = GT.unit_normal(Λ,Ω)
n_Λ = GT.unit_normal(mesh,D-1)
h_Λ = GT.face_diameter_field(Λ)
V = GT.lagrange_space(Ω,order;conformity=:L2)
a_Λ = (u,v) -> GT.∫(
x-> (γ/h_Λ(x))*jump(v,n_Λ,x)⋅jump(u,n_Λ,x)-
jump(v,n_Λ,x)⋅mean(∇(u),x)-
mean(∇(v),x)⋅jump(u,n_Λ,x), dΛ)
x-> (γ/h_Λ(x))*jump(v(x),n_Λ(x))⋅jump(u(x),n_Λ(x))-
jump(v(x),n_Λ(x))⋅mean(∇(u,x))-
mean(∇(v,x))⋅jump(u(x),n_Λ(x)), dΛ)
a_ip = (u,v) -> a_Ω(u,v) + a_Γd(u,v) + a_Λ(u,v)
l_ip = l_ni
p = GT.linear_problem(Float64,V,a_ip,l_ip)
Expand Down
Loading
Loading