Poisson example

Click here to download this notebook

Example of approximations to the Poisson problem with different boundary conditions

The problem to solve is:

\[ \begin{align} -\Delta u &= f \;\;\;\;\; \text{in } \Omega \\ u &= g \;\;\;\;\; \text{in } \partial \Omega_{int} \\ \hat{n} \cdot \nabla u &= h \;\;\;\;\; \text{in } \partial \Omega_{ext} \\ \end{align} \]

For that we impose the weak version of it:

Find $u$ in $V = H^1(\Omega, f) \cap H^1_0(\Omega, \partial \Omega_{int})$ (that is with elements vanishing at $\partial \Omega_{int}$) such that:,

\[ \int_{\Omega} \nabla v \cdot \nabla u \; d\Omega - \int_{\Omega} v \; f \; d\Omega - \oint_{\partial \Omega_{ext}} v \; h \; d\Gamma = 0 \;\;\;\;\; \forall v \;\; \in V \]

If we obtain a solution $u$ and it is sufficiently smooth, we can integrate by parts the first term and obtain:

\[ \int_{\Omega} v \; (-\Delta u - f) \; d\Omega + \oint_{\partial \Omega_{ext}} v \; (\hat{n} \cdot \nabla u - h) \; d\Gamma = 0 \;\;\;\;\; \forall v \;\; \in V \]

Taking $v$ arbitrary but of compact support inside $Ω$ we see that $u$ must satisfy:

\[ -\Delta u = f \;\;\;\;\; \text{in } \Omega, \]

taking now $v$ arbitrary but different from zero in a neighborhood of the exterior boundary, we see that it must satisfy also the Neumann condition:

\[ \hat{n} \cdot \nabla u = h \;\;\;\;\; \text{in } \partial \Omega_{ext}. \]

The Dirichlet condition is automatic for this choice of functional space.

To solve the problem we shall use the a Julia package Gridap.jl. This example is a mixture of several examples in the Gridap tutorial.

using Gridap
using GridapGmsh
#import Gmsh: gmsh # NEVER use #using Gmsh
#mkdir("models")
#mkdir("images")

plot_s = true
if plot_s
    using GridapMakie, GLMakie #Para graficar 
    using FileIO #Gráficos y salidas
    GLMakie.activate!(inline=true) # For windows on the notebook itself. Comment out if you want the as pop-out plots.
end

ERROR: The following 1 direct dependency failed to precompile:

GLMakie 

Failed to precompile GLMakie [e9467ef8-e4e7-5192-8a1a-b1aee30e663a] to "/home/runner/.julia/compiled/v1.11/GLMakie/jl_DdTHt4".
┌ Warning:     OpenGL/GLFW wasn't loaded correctly or couldn't be initialized.
│     This likely means, you're on a headless server without having OpenGL support setup correctly.
│     Have a look at the troubleshooting section in the readme:
│     https://github.com/MakieOrg/Makie.jl/tree/master/GLMakie#troubleshooting-opengl.
└ @ GLMakie ~/.julia/packages/GLMakie/1dQSN/src/gl_backend.jl:4
ERROR: LoadError: InitError: Exception[GLFW.GLFWError(65550, "X11: The DISPLAY environment variable is missing"), ErrorException("glfwInit failed")]
Stacktrace:
  [1] __init__()
    @ GLFW ~/.julia/packages/GLFW/a10jC/src/GLFW.jl:69
  [2] run_module_init(mod::Module, i::Int64)
    @ Base ./loading.jl:1378
  [3] register_restored_modules(sv::Core.SimpleVector, pkg::Base.PkgId, path::String)
    @ Base ./loading.jl:1366
  [4] _include_from_serialized(pkg::Base.PkgId, path::String, ocachepath::String, depmods::Vector{Any}, ignore_native::Nothing; register::Bool)
    @ Base ./loading.jl:1254
  [5] _include_from_serialized (repeats 2 times)
    @ ./loading.jl:1210 [inlined]
  [6] _require_search_from_serialized(pkg::Base.PkgId, sourcepath::String, build_id::UInt128, stalecheck::Bool; reasons::Dict{String, Int64}, DEPOT_PATH::Vector{String})
    @ Base ./loading.jl:2057
  [7] _require(pkg::Base.PkgId, env::String)
    @ Base ./loading.jl:2527
  [8] __require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2388
  [9] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
 [10] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
 [11] _require_prelocked(uuidkey::Base.PkgId, env::String)
    @ Base ./loading.jl:2375
 [12] macro expansion
    @ ./loading.jl:2314 [inlined]
 [13] macro expansion
    @ ./lock.jl:273 [inlined]
 [14] __require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2271
 [15] #invoke_in_world#3
    @ ./essentials.jl:1089 [inlined]
 [16] invoke_in_world
    @ ./essentials.jl:1086 [inlined]
 [17] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:2260
 [18] top-level scope
    @ ~/.julia/packages/GLMakie/1dQSN/src/gl_backend.jl:2
 [19] include(mod::Module, _path::String)
    @ Base ./Base.jl:557
 [20] include(x::String)
    @ GLMakie ~/.julia/packages/GLMakie/1dQSN/src/GLMakie.jl:1
 [21] top-level scope
    @ ~/.julia/packages/GLMakie/1dQSN/src/GLMakie.jl:87
 [22] include
    @ ./Base.jl:557 [inlined]
 [23] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt128}}, source::Nothing)
    @ Base ./loading.jl:2881
 [24] top-level scope
    @ stdin:6
during initialization of module GLFW
in expression starting at /home/runner/.julia/packages/GLMakie/1dQSN/src/gl_backend.jl:1
in expression starting at /home/runner/.julia/packages/GLMakie/1dQSN/src/GLMakie.jl:1
in expression starting at stdin:

We will use a mesh previously built with the gmsh library. In the models directory you will find a script with the .geo extension (rectangle_hole.geo), which was used to build the example. Based on this script, and following the gmsh tutorial, you can build different meshes. Other libraries can also be used to build meshes. These are imported into the Gridap infrastructure, and with them the triangulation to be used is constructed. Note that in the script, the two boundaries are named: the external (rectangular) one as "ext" and the internal (circle) one as "int".

#model = GmshDiscreteModel("models/rectangle_hole_fine.msh")
model = GmshDiscreteModel("models/rectangle_hole_coarse.msh")
#model = GmshDiscreteModel("models/rectangle_hole_finer.msh")

ERROR: Msh file not found: models/rectangle_hole_coarse.msh

Ω = Triangulation(model)

ERROR: UndefVarError: `model` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

degree = 2
dΩ = Measure(Ω,degree)

ERROR: UndefVarError: `Ω` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

if plot_s
    fig, ax = plot(Ω)
    ax.aspect = AxisAspect(2)
    wireframe!(Ω, color=:black, linewidth=1)
    scatter!(Ω, marker=:star8, markersize=4, color=:blue)
    fig
end

ERROR: UndefVarError: `Ω` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

full_dirichlet = false
int_dirichlet = false
#full_dirichlet = true
int_dirichlet = true

true

Once we have the mesh, we begin to define the finite elements that we will use. In this case, we will use Lagrangian elements of order 1 that will satisfy a Dirichlet condition on the region $\partial \Omega_{int}$. When the mesh is built, this region has been marked as the interior boundary of the rectangle with the tag "int".

order = 2
reffe = ReferenceFE(lagrangian,Float64,order)
if int_dirichlet
    dirichlet_tags="int, " 
elseif full_dirichlet
    dirichlet_tags=["int","ext"]
end
V = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags = dirichlet_tags)

ERROR: UndefVarError: `model` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

In this example, we are going to test our code using a solution given by the function $u_e(x,y) = 4((x-x_0)^2 - (y-y_0)^3) - 5y$.

Thus, $\Delta u_e = -f = 8 - 24(y-y0)$ and the Neumann condition will be given by $h = \hat{n} \cdot \nabla u_e$.

x0 = 0.5
y0 = 0.5
ue(x) = 4*((x[1]-x0)^2 - (x[2]-y0)^3) - 5.0*x[2] # "exact solution"

if plot_s
    fig, ax, plt = plot(Ω, ue, shading=false)
    ax.aspect = AxisAspect(2)
    Colorbar(fig[1,2], plt)
    fig
end

ERROR: UndefVarError: `Ω` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

# internal Dirichlet boundary condition
U = TrialFESpace(V,ue)
Γ₁ = BoundaryTriangulation(model,tags=dirichlet_tags)
n₁ = get_normal_vector(Γ₁)

if plot_s
    fig, ax , plt  = plot(Γ₁,ue, colormap=:heat, linewidth=10)
    ax.aspect = AxisAspect(2)
    Colorbar(fig[1,2], plt)
    fig
end

ERROR: UndefVarError: `V` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

if int_dirichlet
    neumann_tags = "ext"
    #neumanntags = "int"
    Γ = BoundaryTriangulation(model,tags=neumann_tags)
    dΓ = Measure(Γ,degree)
    nb = get_normal_vector(Γ)
end

ERROR: UndefVarError: `model` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

To make sure every thing is OK, we plot the values of $\hat{n}\cdot \nabla u_e$ at the external boundary.

if int_dirichlet && plot_s
    fig, ax , plt = plot(Γ, (nb ⋅ ∇(ue)), colormap=:algae, linewidth=10)
    ax.aspect = AxisAspect(2)
    Colorbar(fig[1,2], plt)
    fig
end

ERROR: UndefVarError: `nb` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

We now define the weak problem in abstract form:

f(x) = -8.0 + 24*(x[2] - y0) # source
#h(x) =  #external Neumann bc.
a(u,v) = ∫( ∇(v)⋅∇(u) )*dΩ  #  in a(u,v) goes all the terms linear in the unknown, u. 
if full_dirichlet
    b(v) = ∫(v*f )*dΩ # 
elseif int_dirichlet
    b(v) = ∫(v*f )*dΩ + ∫(v*(nb ⋅ ∇(ue)))*dΓ #  here everything else.
end

b (generic function with 1 method)

From here the package Gridap.jl generates a equation system of the form $Ax=b$ and solves it for the finite element version of $u$.

op = AffineFEOperator(a,b,U,V)

ERROR: UndefVarError: `U` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

ls = LUSolver()
solver = LinearFESolver(ls)

LinearFESolver()

uh = solve(solver,op)

ERROR: UndefVarError: `op` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

if plot_s
    fig, ax, plt = plot(Ω, uh, shading=NoShading)
    ax.aspect = AxisAspect(2)
    Colorbar(fig[1,2], plt)
    fig
end

ERROR: UndefVarError: `NoShading` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name may be made accessible by importing Makie in the current active module Main

if full_dirichlet
    writevtk(Ω,"images/solución_dir",cellfields=["uh_dir" => uh])
elseif int_dirichlet
    writevtk(Ω,"images/solución_newmann",cellfields=["uh_neu" => uh])
end

ERROR: UndefVarError: `uh` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

We now are going to validate our approximation by comparing it with the exact solution. For that we introduce several tools:

e = ue - uh

ERROR: UndefVarError: `uh` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

if plot_s
    fig, ax, plt = plot(Ω, e
        , shading=NoShading
        )
    ax.aspect = AxisAspect(2)
    Colorbar(fig[1,2], plt)
    fig
end
fig

ERROR: UndefVarError: `NoShading` not defined in `Main`
Suggestion: check for spelling errors or missing imports.
Hint: a global variable of this name may be made accessible by importing Makie in the current active module Main

if full_dirichlet
    writevtk(Ω,"images/error_dir",cellfields=["e_dir" => e])
elseif int_dirichlet
    writevtk(Ω,"images/error_newmann",cellfields=["e_neu" => e])
end

ERROR: UndefVarError: `e` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

Next we compute the $L^2$ and $H^1$ norm for the error.

el2 = sqrt(sum( ∫( e*e )*dΩ ))
println("l2 error = ",el2)
eh1 = sqrt(sum( ∫( e*e + ∇(e)⋅∇(e) )*dΩ ))
uh1 = sqrt(sum( ∫( uh*uh + ∇(uh)⋅∇(uh) )*dΩ ))
println("h1 error = ",eh1/uh1)

ERROR: UndefVarError: `e` not defined in `Main`
Suggestion: check for spelling errors or missing imports.

if plot_s
    fig, ax, plt = plot(Ω, ∇(e)⋅∇(e), shading=NoShading)
    ax.aspect = AxisAspect(2)
    Colorbar(fig[1,2], plt)
    fig
end

ERROR: UndefVarError: `e` not defined in `Main`
Suggestion: check for spelling errors or missing imports.