Skip to content

Commit

Permalink
gk -> GT
Browse files Browse the repository at this point in the history
  • Loading branch information
fverdugo committed Jul 6, 2024
1 parent ae8c105 commit de3a4cd
Show file tree
Hide file tree
Showing 14 changed files with 1,190 additions and 1,190 deletions.
70 changes: 35 additions & 35 deletions GalerkinToolkitExamples/src/poisson.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@ This implements a Poisson solver with several methods
"""
module Poisson

import GalerkinToolkit as gk
import GalerkinToolkit as GT
using GalerkinToolkit:
import PartitionedSolvers as ps
import ForwardDiff
Expand All @@ -29,8 +29,8 @@ 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) = gk.call(laplacian,u)
(u) = gk.call(gradient,u)
Δ(u) = GT.call(laplacian,u)
(u) = GT.call(gradient,u)
Δ(u,x) = Δ(u)(x)
(u,x) = (u)(x)

Expand All @@ -42,17 +42,17 @@ function main_automatic(params)
results = Dict{Symbol,Any}()

mesh = params[:mesh]
D = gk.num_dims(mesh)
Ω = gk.interior(mesh,physical_names=params[:domain_tags])
Γd = gk.boundary(mesh;physical_names=params[:dirichlet_tags])
Γn = gk.boundary(mesh;physical_names=params[:neumann_tags])
D = GT.num_dims(mesh)
Ω = GT.interior(mesh,physical_names=params[:domain_tags])
Γd = GT.boundary(mesh;physical_names=params[:dirichlet_tags])
Γn = GT.boundary(mesh;physical_names=params[:neumann_tags])
integration_degree = params[:integration_degree]
= gk.measure(Ω,integration_degree)
dΓn = gk.measure(Γn,integration_degree)
= GT.measure(Ω,integration_degree)
dΓn = GT.measure(Γn,integration_degree)

u = gk.analytical_field(params[:u],Ω)
u = GT.analytical_field(params[:u],Ω)
f(x) = -Δ(u,x)
n_Γn = gk.unit_normal(Γn,Ω)
n_Γn = GT.unit_normal(Γn,Ω)
g(x) = n_Γn(x)(u,x)

interpolation_degree = params[:interpolation_degree]
Expand All @@ -63,31 +63,31 @@ function main_automatic(params)

if params[:discretization_method] !== :continuous_galerkin
conformity = :L2
gk.label_interior_faces!(mesh;physical_name="__INTERIOR_FACES__")
Λ = gk.skeleton(mesh;physical_names=["__INTERIOR_FACES__"])
= gk.measure(Λ,integration_degree)
n_Λ = gk.unit_normal(Λ,Ω)
h_Λ = gk.face_diameter_field(Λ)
GT.label_interior_faces!(mesh;physical_name="__INTERIOR_FACES__")
Λ = GT.skeleton(mesh;physical_names=["__INTERIOR_FACES__"])
= GT.measure(Λ,integration_degree)
n_Λ = GT.unit_normal(Λ,Ω)
h_Λ = GT.face_diameter_field(Λ)
else
conformity = :default
end

if params[:dirichlet_method] === :strong
V = gk.lagrange_space(Ω,interpolation_degree;conformity,dirichlet_boundary=Γd)
uh = gk.zero_field(Float64,V)
gk.interpolate_dirichlet!(u,uh)
V = GT.lagrange_space(Ω,interpolation_degree;conformity,dirichlet_boundary=Γd)
uh = GT.zero_field(Float64,V)
GT.interpolate_dirichlet!(u,uh)
else
n_Γd = gk.unit_normal(Γd,Ω)
h_Γd = gk.face_diameter_field(Γd)
dΓd = gk.measure(Γd,integration_degree)
V = gk.lagrange_space(Ω,interpolation_degree;conformity)
uh = gk.zero_field(Float64,V)
n_Γd = GT.unit_normal(Γd,Ω)
h_Γd = GT.face_diameter_field(Γd)
dΓd = GT.measure(Γd,integration_degree)
V = GT.lagrange_space(Ω,interpolation_degree;conformity)
uh = GT.zero_field(Float64,V)
end

if params[:dirichlet_method] === :multipliers
Q = gk.lagrange_space(Γd,interpolation_degree-1;conformity=:L2)
Q = GT.lagrange_space(Γd,interpolation_degree-1;conformity=:L2)
VxQ = V × Q
uh_qh = gk.zero_field(Float64,VxQ)
uh_qh = GT.zero_field(Float64,VxQ)
uh, qh = uh_qh
end

Expand Down Expand Up @@ -132,7 +132,7 @@ function main_automatic(params)

@timeit timer "assembly" begin
# TODO give a hint when dirichlet BCS are homogeneous or not present
x,A,b = gk.linear_problem(Uh,A,L)
x,A,b = GT.linear_problem(Uh,A,L)
end

@timeit timer "solver" begin
Expand All @@ -150,14 +150,14 @@ function main_automatic(params)
end

@timeit timer "vtk" if params[:export_vtu]
gk.vtk_plot(params[:example_path]*"",Ω;refinement=4) do plt
gk.plot!(plt,u;label="u")
gk.plot!(plt,f;label="f")
gk.plot!(plt,uh;label="uh")
GT.vtk_plot(params[:example_path]*"",Ω;refinement=4) do plt
GT.plot!(plt,u;label="u")
GT.plot!(plt,f;label="f")
GT.plot!(plt,uh;label="uh")
end
gk.vtk_plot(params[:example_path]*"_Γn",Γn;refinement=4) do plt
gk.plot!(plt,n_Γn;label="n")
gk.plot!(plt,g;label="g")
GT.vtk_plot(params[:example_path]*"_Γn",Γn;refinement=4) do plt
GT.plot!(plt,n_Γn;label="n")
GT.plot!(plt,g;label="g")
end
end

Expand Down Expand Up @@ -193,7 +193,7 @@ function default_params()
params[:discretization_method] = :continuous_galerkin
params[:verbose] = true
params[:timer] = TimerOutput()
params[:mesh] = gk.cartesian_mesh((0,1,0,1),(10,10))
params[:mesh] = GT.cartesian_mesh((0,1,0,1),(10,10))
params[:u] = (x) -> sum(x)
params[:domain_tags] = ["interior"]
params[:dirichlet_tags] = ["boundary"]
Expand Down
4 changes: 2 additions & 2 deletions GalerkinToolkitExamples/test/poisson_tests.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
module PoissonTests

import GalerkinToolkit as gk
import GalerkinToolkit as GT
using GalerkinToolkitExamples: Poisson
using Test

Expand All @@ -9,7 +9,7 @@ tol = 1.0e-10
# TODO do not use a large value of n here.
# The code now is VERY slow.
n = 2
mesh = gk.cartesian_mesh((0,2,0,2),(n,n))
mesh = GT.cartesian_mesh((0,2,0,2),(n,n))
for discretization_method in (:interior_penalty,:continuous_galerkin)
for dirichlet_method in (:multipliers,:nitsche,:strong)
params = Dict{Symbol,Any}()
Expand Down
66 changes: 33 additions & 33 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -47,28 +47,28 @@ Discuss with the package authors before working on any non-trivial contribution.
The following code solves a Laplace PDE with Dirichlet boundary conditions.

```julia
import GalerkinToolkit as gk
import GalerkinToolkit as GT
import ForwardDiff
using LinearAlgebra
mesh = gk.mesh_from_gmsh("assets/demo.msh")
gk.label_boundary_faces!(mesh;physical_name="boundary")
Ω = gk.interior(mesh)
Γd = gk.boundary(mesh;physical_names=["boundary"])
mesh = GT.mesh_from_gmsh("assets/demo.msh")
GT.label_boundary_faces!(mesh;physical_name="boundary")
Ω = GT.interior(mesh)
Γd = GT.boundary(mesh;physical_names=["boundary"])
k = 2
V = gk.lagrange_space(Ω,k;dirichlet_boundary=Γd)
uh = gk.zero_field(Float64,V)
u = gk.analytical_field(sum,Ω)
gk.interpolate_dirichlet!(u,uh)
= gk.measure(Ω,2*k)
V = GT.lagrange_space(Ω,k;dirichlet_boundary=Γd)
uh = GT.zero_field(Float64,V)
u = GT.analytical_field(sum,Ω)
GT.interpolate_dirichlet!(u,uh)
= GT.measure(Ω,2*k)
gradient(u) = x->ForwardDiff.gradient(u,x)
(u,x) = gk.call(gradient,u)(x)
a(u,v) = gk.( x->(u,x)(v,x), dΩ)
(u,x) = GT.call(gradient,u)(x)
a(u,v) = GT.( x->(u,x)(v,x), dΩ)
l(v) = 0
x,A,b = gk.linear_problem(uh,a,l)
x,A,b = GT.linear_problem(uh,a,l)
x .= A\b
gk.vtk_plot("results",Ω) do plt
gk.plot!(plt,u;label="u")
gk.plot!(plt,uh;label="uh")
GT.vtk_plot("results",Ω) do plt
GT.plot!(plt,u;label="u")
GT.plot!(plt,uh;label="uh")
end
```

Expand All @@ -78,33 +78,33 @@ This code solves the same boundary value problem, but using an auxiliary field o
multipliers to impose Dirichlet boundary conditions.

```julia
import GalerkinToolkit as gk
import GalerkinToolkit as GT
import ForwardDiff
using LinearAlgebra
mesh = gk.mesh_from_gmsh("assets/demo.msh")
gk.label_boundary_faces!(mesh;physical_name="boundary")
Ω = gk.interior(mesh)
Γd = gk.boundary(mesh;physical_names=["boundary"])
mesh = GT.mesh_from_gmsh("assets/demo.msh")
GT.label_boundary_faces!(mesh;physical_name="boundary")
Ω = GT.interior(mesh)
Γd = GT.boundary(mesh;physical_names=["boundary"])
k = 2
V = gk.lagrange_space(Ω,k)
Q = gk.lagrange_space(Γd,k-1; conformity=:L2)
V = GT.lagrange_space(Ω,k)
Q = GT.lagrange_space(Γd,k-1; conformity=:L2)
VxQ = V × Q
dΓd = gk.measure(Γd,2*k)
dΓd = GT.measure(Γd,2*k)
gradient(u) = x->ForwardDiff.gradient(u,x)
(u,x) = gk.call(gradient,u)(x)
(u,x) = GT.call(gradient,u)(x)
a((u,p),(v,q)) =
gk.( x->(u,x)(v,x), dΩ) +
gk.(x->
GT.( x->(u,x)(v,x), dΩ) +
GT.(x->
(u(x)+p(x))*(v(x)+q(x))
-u(x)*v(x)-p(x)*q(x), dΓd)
l((v,q)) = gk.(x->u(x)*q(x), dΓd)
uh_qh = gk.zero_field(Float64,VxQ)
x,A,b = gk.linear_problem(uh_qh,a,l)
l((v,q)) = GT.(x->u(x)*q(x), dΓd)
uh_qh = GT.zero_field(Float64,VxQ)
x,A,b = GT.linear_problem(uh_qh,a,l)
x .= A\b
uh,qh = uh_qh
gk.vtk_plot("results",Ω) do plt
gk.plot!(plt,u;label="u")
gk.plot!(plt,uh;label="uh")
GT.vtk_plot("results",Ω) do plt
GT.plot!(plt,u;label="u")
GT.plot!(plt,uh;label="uh")
end
```

2 changes: 1 addition & 1 deletion src/GalerkinToolkit.jl
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
module GalerkinToolkit
const gk = GalerkinToolkit
const GT = GalerkinToolkit

using WriteVTK
using FastGaussQuadrature
Expand Down
Loading

0 comments on commit de3a4cd

Please sign in to comment.