• Click here to download this notebook

Example of solving the Poisson problem with different boundary conditions and manipulation of the solution for the capacitance calculation problem

The problem to be solved 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} \]

To do this, we will impose its weak form:

Find $u$ in $H^1(\Omega, f))$ (that is, with Dirichlet boundary conditions on $\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 H^1_0(\Omega) \]

If we obtain a $u$ satisfying this equation, and it is sufficiently smooth, then we can integrate the first term by parts 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 H^1_0(\Omega) \]

Taking $v$ arbitrary but with compact support, we see that $u$ must satisfy:

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

and taking $v$ arbitrary we also see that the Neumann condition must be satisfied,

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

The Dirichlet condition is automatic due to the choice of space.

We will then use the solution found for a capacitance problem.

To solve the problem, we will use the infrastructure of the Julia package Gridap.jl. This example is a compilation of several examples from the package's tutorial.

using Gridap
using FileIO #Gráficos y salidas
using GridapGmsh
#using gmsh #]add https://github.com/koehlerson/gmsh.jl.git
#NEVER USE using Gmsh!!!
#import Gmsh: gmsh

#plot_s = false #true
plot_s = true
if plot_s
    #using Pkg; Pkg.activate("./gridap_makie")
    #Pkg.add(Pkg.PackageSpec(;name="Makie", version="0.12"))
    using Makie
    using GridapMakie, GLMakie #Para graficar 
    GLMakie.activate!(inline=true) # For windows on the notebook itself. Comment out if you want the as pop-out plots.
end

ERROR: ArgumentError: Package Makie not found in current path.
- Run `import Pkg; Pkg.add("Makie")` to install the Makie package.

We will use meshes built with the gmsh library, through the script mesh_generator.jl. Note that in the script, the three boundaries are named: the external (rectangular) one as "ext" and the internal ones as "inner_circle" and "inner_square". Based on the script and following the gmsh tutorial, you can build other meshes as well. 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.

create_grid = false

#Grid resolution:
res = 1 # 1, 2 or 3

1

if create_grid

include("mesh_generator.jl")
grid_type = "rectangle_hole_square"

#Resoluciones:
res = 1

if res == 1
    lc = 1e-1
    lc_f = 0.25e-1
    name = grid_type * "_coarse"
elseif res == 2 
    lc = 5e-2
    lc_f = 2.5e-2
    name = grid_type * "_intermediate"
elseif res == 3
    lc = 1e-2
    lc_f = 0.25e-2
    name = grid_type * "_finner"
end
    
#Lados exteriores de la grilla rectangular
side_x = 2
side_y = 1

#Rectangulo interior
rec_base = 0.25  #Coordenada y de la base
rec_top = 0.75   #Coordenada y del lado superior
rec_left = 1.25  #Coordenada x del lado izquierdo
rec_right = 1.75 #Coordenada x del lado derecho

#Circulo interior
circ_center_x = 0.5  #Coordenada x del centro
circ_center_y = 0.5  #Coordenada y del centro
circ_radius = 0.25   #Radio



p = (name, side_x, side_y, circ_center_x, circ_center_y, circ_radius, rec_base, rec_top, rec_left, rec_right, lc, lc_f)

model = make_model(grid_type, p)

end

You can also upload the grid once is created and stored in the directory: models

if create_grid == false
    model = GmshDiscreteModel("models/rectangle_hole_square_coarse.msh")
end

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

Ω = Triangulation(model)

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

degree = 3
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
display(fig)

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

We can also plot the different boundaries:

boundary_tags = ["inner_circle", "inner_square", "ext"]

Γ = BoundaryTriangulation(model,tags=boundary_tags)
dΓ = Measure(Γ,degree)

if plot_s
    fig, ax = plot(Γ, linewidth=8)
    ax.aspect = AxisAspect(2)
    wireframe!(Γ, color=:black, linewidth=1)
    fig
end
display(fig)

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

We are going to choose two simple problems to calculate the capacitance matrix of a set of conductors. We will take the conductors as the two bodies: the circle and the square, and we will set constant potential conditions. The external boundary will be considered as infinity and we will always set the potential to zero there.

capacity_cs = false # potencial 1 en el círculo y potencial 0 en el cuadrado.
capacity_sc = false # potencial 0 en el círculo y potencial 1 en el cuadrado.
capacity_cs = true
#capacity_sc = 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 = 1
reffe = ReferenceFE(lagrangian,Float64,order)

dirichlet_tags= ["inner_circle", "inner_square","ext"] 

V = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags = dirichlet_tags)
#V = TestFESpace(model,reffe;conformity=:L2,dirichlet_tags = dirichlet_tags) #no funciona la inversión.

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

We now assign the boundary values:

# internal Dirichlet boundary condition
g(x) = 1.0 # esta puede ser una función de x (vector posición)
if capacity_cs
    U = TrialFESpace(V,[g 0.0 0.0])
elseif capacity_sc
    U = TrialFESpace(V,[0.0 g 0.0])
end

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

Next we define the weak problem in abstract form:

f(x) = 0 # in this case we are taking the source to vanish, but in principle we could put a charge distribution inside.

a(u,v) = ∫( ∇(v)⋅∇(u) )*dΩ  # in a(u,v) goes all the terms linear in the unknown, u. 

b(v) = ∫(v*f )*dΩ # here everything else.

b (generic function with 1 method)

From this point on the package Gridap.jl generates a system of the form: $Ax=b$ and solves it for $u$ in the finite element version given. It first computes the matrices:

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

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

Then it solves with some of the possible known strategies. Here we give two of them, the first is just an L-U decomposition (Lower triangular Unitary). The second it the / operator of Julia, which in general would use the most powerful method for the case.

ls = LUSolver()
lb = BackslashSolver() # x = A \ b
solver = LinearFESolver(lb)

LinearFESolver()

And we solve it. We call the solution $u_h$.

uh = solve(solver,op)
#uh = solve(op)

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

Now we can plot the found approximation and also save it in a format that can be given to Paraview or VisIt for ploting and processing.

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

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

if capacity_cs
    writevtk(Ω,"images/solucion_cs_$res",cellfields=["uh_cs_$res" => uh])
    writevtk(Ω,"images/grad_cs_$res",cellfields=["grad_uh_cs_$res" => ∇(uh)])
elseif capacity_sc
    writevtk(Ω,"images/solucion_sc_$res",cellfields=["uh_sc_$res" => uh])
    writevtk(Ω,"images/grad_cs_$res",cellfields=["grad_uh_cs_$res" => ∇(uh)])
end

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

Once we obtain the approximations we can compute some physical quantities of interest. Like, for instance, the charges on each of the conductors. Here we compute the charge in the exterior and interior. Note the sign change for the interior, this is due to the fact that the package takes alwais the normal pointing to the exterior of the mesh.

Recall the charge is defined as:

\begin{equation} Qi = \int{\partial \Omegai} \sigma dS = \frac{1}{4\pi}\int{\partial \Omega_i} E \cdot \; dS \end{equation}

While the capacity matrix as:

\begin{equation} Qi = C{ij}V^j \end{equation}

Γ_ext = BoundaryTriangulation(model,tags="ext")
dΓ_ext = Measure(Γ_ext,degree)
nb_ext = get_normal_vector(Γ_ext) # this is the normal to the boundary.
Q_ext = -sum(∫((nb_ext ⋅ ∇(uh)))*dΓ_ext)/4/π

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

Γ_square = BoundaryTriangulation(model,tags="inner_square")
dΓ_square = Measure(Γ_square,degree)
nb_square = get_normal_vector(Γ_square)
Q_square = sum(∫((nb_square ⋅ ∇(uh)))*dΓ_square)/4/π

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

Γ_circle = BoundaryTriangulation(model,tags="inner_circle")
dΓ_circle = Measure(Γ_circle,degree)
nb_circle = get_normal_vector(Γ_circle)
Q_circle = sum(∫((nb_circle ⋅ ∇(uh)))*dΓ_circle)/4/π

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

Q_ext - Q_circle - Q_square # this number should vanish if there are no sources (charges) present.

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

We show an image with the program VisIt: