GalerkinToolkit

This package aims at providing a fully-fledged finite-element toolbox in pure Julia, with support for different computing systems from laptops to supercomputers and GPUs.

NB. This package is work in progress; a proof-of-concept API is already available (for CPUs). The package is not production ready at this point. Planned performance and documentation improvements are needed.

This package follows a new approach to implement finite-element methods based on the lessons learned in the Gridap project.

Since July 2024, this package is being developed with support from the Netherlands eScience Center under grant ID NLESC.SS.2023.008.

• STABLE — Documentation for the most recently tagged version.
• LATEST — Documentation for the in-development version.

• You can open a new discussion to ask questions here.
• If you have found a bug, open an issue here. Do not forget to include a (minimal) reproducer.

This package is under active development and there are several ways to contribute:

• by enhancing the documentation (e.g., fixing typos, enhancing doc strings, adding examples).
• by addressing one of the issues waiting for help.
• by adding more tests to increase the code coverage.
• by extending the current functionality. In this case, open a discussion here to coordinate with the package maintainers before proposing significant changes.

Discuss with the package authors before working on any non-trivial contribution.

The following code solves a Laplace PDE with Dirichlet boundary conditions.

import GalerkinToolkit as GT
import ForwardDiff
using LinearAlgebra
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 = GT.lagrange_space(Ω,k;dirichlet_boundary=Γd)
uhd = GT.dirichlet_field(Float64,V)
u = GT.analytical_field(sum,Ω)
GT.interpolate_dirichlet!(u,uhd)
dΩ = GT.measure(Ω,2*k)
gradient(u) = x->ForwardDiff.gradient(u,x)
∇(u,x) = GT.call(gradient,u)(x)
a(u,v) = GT.∫( x->∇(u,x)⋅∇(v,x), dΩ)
l(v) = 0
x,A,b = GT.linear_problem(uhd,a,l)
x .= A\b
uh = GT.solution_field(uhd,x)
GT.vtk_plot("results",Ω) do plt
    GT.plot!(plt,u;label="u")
    GT.plot!(plt,uh;label="uh")
end

This code solves the same boundary value problem, but using an auxiliary field of Lagrange multipliers to impose Dirichlet boundary conditions.

import GalerkinToolkit as GT
import ForwardDiff
using LinearAlgebra
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 = GT.lagrange_space(Ω,k)
Q = GT.lagrange_space(Γd,k-1; conformity=:L2)
VxQ = V × Q
dΓd = GT.measure(Γd,2*k)
gradient(u) = x->ForwardDiff.gradient(u,x)
∇(u,x) = GT.call(gradient,u)(x)
a((u,p),(v,q)) =
    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)) = GT.∫(x->u(x)*q(x), dΓd)
x,A,b = GT.linear_problem(Float64,VxQ,a,l)
x .= A\b
uh,qh = GT.solution_field(VxQ,x)
GT.vtk_plot("results",Ω) do plt
    GT.plot!(plt,u;label="u")
    GT.plot!(plt,uh;label="uh")
end