Two dimensional examples

A second order equation

Consider the following equation:

\[ \left\{ \begin{array}{l} -\Delta u + u = f \textrm{ in } \Omega \\ u = 0 \textrm{ on } \Omega, \end{array} \right.\]

with \( \Omega = (0, 1) \times (0, 1)\).

using SymFEL
using SymPy
using LinearAlgebra
using SparseArrays
import gmsh
using WriteVTK

x, h = symbols("x h")

discretization parameters

We use the mesh square.msh (in gmsh format).

This mesh is formed by quadratic quad elements.

gmsh.jl module should be installed manually following the instructions on gmsh web page

gmsh.initialize()
gmsh.open("square.msh")
nodes = gmsh.model.mesh.getNodes()
nodes_label = convert(Array{Int64}, nodes[1])
nodes_N = length(nodes_label)
nodes_coordinate = reshape(nodes[2], (3, nodes_N))
nodes_boundary = convert(Array{Int64}, unique(gmsh.model.mesh.getNodesByElementType(1)[1]))
nodes_boundary_N = length(nodes_boundary)

elements = gmsh.model.mesh.getElements()

elements_bound_label = convert(Array{Int64}, elements[2][1])
elements_bound_N = length(elements_bound_label)

elements_int_label = convert(Array{Int64}, elements[2][2])
elements_int_N = length(elements_int_label)

elements_bound = reshape(convert(Array{Int64}, elements[3][1]), (3, elements_bound_N))
elements_int = reshape(convert(Array{Int64}, elements[3][2]), (9, elements_int_N))
#gmsh.fltk.run()
gmsh.finalize()

Info    : Reading 'square.msh'...
Info    : 9 entities
Info    : 1681 nodes
Info    : 480 elements
Info    : Done reading 'square.msh'

elementary matrices - P2 x P2

elem_Mxy = SymFEL.get_square_lagrange_em((2, 2), (0, 0), (0, 0))
elem_Kxy = SymFEL.get_square_lagrange_em((2, 2), (1, 0), (1, 0)) +
    SymFEL.get_square_lagrange_em((2, 2), (0, 1), (0, 1))

dx = norm(nodes_coordinate[:, elements_bound[1,1]] - nodes_coordinate[:, elements_bound[2,1]])

elem_Kxy_dx = convert(Matrix{Float64}, elem_Kxy.subs(h, dx))
elem_Mxy_dx = convert(Matrix{Float64}, elem_Mxy.subs(h, dx));

global matrices

K = SymFEL.assemble_squaremesh_FE_matrix(elem_Kxy_dx, elements_int, order1=2, order2=2)
M = SymFEL.assemble_squaremesh_FE_matrix(elem_Mxy_dx, elements_int, order1=2, order2=2)

right hand side

We consider \(f(x, y) = (2\pi^2 + 1) \sin(\pi x) \sin(\pi y)\).

This corresponds to the exact solution \(u(x,y) = \sin(\pi x) \sin(\pi y)\).

f = (2*pi^2 + 1) * sin.(pi * nodes_coordinate[1,:]) .* sin.(pi * nodes_coordinate[2,:])

F = M * f

A = K + M

boundary condition

tgv = 1e30
A[nodes_boundary, nodes_boundary] += tgv * sparse(Matrix{Float64}(I, nodes_boundary_N, nodes_boundary_N))

solve the linear system

u = A \ F

compute and display the error

u_exact = sin.(pi*nodes_coordinate[1,:]) .* sin.(pi*nodes_coordinate[2,:])
err = u - u_exact
println("L2 error : ", sqrt(err' * M * err))
println("H1 error : ", sqrt(err' * K * err))

L2 error : 4.0733705143585556e-7
H1 error : 1.0411926041117035e-5

export to vtk

cells = [MeshCell(VTKCellTypes.VTK_QUADRATIC_QUAD, elements_int[1:8, i]) for i = 1:elements_int_N]
points_x = nodes_coordinate[1, :]
points_y = nodes_coordinate[2, :]
vtkfile = vtk_grid("ex3-output", points_x, points_y, cells)

vtkfile["my_point_data", VTKPointData()] = u
outfiles = vtk_save(vtkfile)

1-element Vector{String}:
 "ex3-output.vtu"

affichage de la solution

L'affichage est obtenu à l'aide de Paraview à partir du fichier "ex3-output.vtu"

A fourth order equation

The corresponding program can be downloaded here.

Consider the following equation:

\[ \left\{ \begin{array}{l} \Delta^2 u = f \textrm{ in } \Omega \\ u = \Delta u = 0 \textrm{ on } \Omega, \end{array} \right.\]

with \( \Omega = (0, 1) \times (0, 1)\).

using Revise
using SymFEL
using SymPy
using LinearAlgebra
using SparseArrays
import gmsh
using WriteVTK

x, h = symbols("x h")

discretization parameters

We use the mesh simple-square.msh (in gmsh format).

This mesh is formed by square elements.

we use gmsh module for read the mesh

gmsh.initialize()
gmsh.open("square-simple.msh")
nodes = gmsh.model.mesh.getNodes()
nodes_label = convert(Array{Int64}, nodes[1])
nodes_N = length(nodes_label)
nodes_coordinate = reshape(nodes[2], (3, nodes_N))
nodes_boundary = convert(Array{Int64}, unique(gmsh.model.mesh.getNodesByElementType(1)[1]))
nodes_boundary_N = length(nodes_boundary)

elements = gmsh.model.mesh.getElements()

elements_bound_label = convert(Array{Int64}, elements[2][1])
elements_bound_N = length(elements_bound_label)

elements_int_label = convert(Array{Int64}, elements[2][2])
elements_int_N = length(elements_int_label)

elements_bound = reshape(convert(Array{Int64}, elements[3][1]), (2, elements_bound_N))
elements_int = reshape(convert(Array{Int64}, elements[3][2]), (4, elements_int_N))
#gmsh.fltk.run()
gmsh.finalize()

Info    : Reading 'square-simple.msh'...
Info    : 9 entities
Info    : 441 nodes
Info    : 480 elements
Info    : Done reading 'square-simple.msh'

elementary matrices - Hermite P3 x Hermite P3

elem_M = SymFEL.get_square_hermite_em((3, 3), (0, 0), (0, 0))
elem_K = SymFEL.get_square_hermite_em((3, 3), (2, 0), (2, 0)) +
    SymFEL.get_square_hermite_em((3, 3), (2, 0), (0, 2)) +
    SymFEL.get_square_hermite_em((3, 3), (0, 2), (2, 0)) +
    SymFEL.get_square_hermite_em((3, 3), (0, 2), (0, 2))

dx = norm(nodes_coordinate[:, elements_bound[1,1]] - nodes_coordinate[:, elements_bound[2,1]])
elem_K_dx = convert(Matrix{Float64}, elem_K.subs(h, dx))
elem_M_dx = convert(Matrix{Float64}, elem_M.subs(h, dx));

global matrices

K = SymFEL.assemble_squaremesh_FE_matrix(elem_K_dx, elements_int,
                                           order1=1, order2=1,
                                           dof1=4, dof2=4)
M = SymFEL.assemble_squaremesh_FE_matrix(elem_M_dx, elements_int,
                                           order1=1, order2 = 1,
                                           dof1=4, dof2=4)

Right hand side

f = zeros(Float64, 4*nodes_N)
f[4*((1:nodes_N) .- 1) .+ 1] = 4*pi^4*sin.(pi * nodes_coordinate[1,:]) .* sin.(pi * nodes_coordinate[2,:])
f[4*((1:nodes_N) .- 1) .+ 2] = 4*pi^5*cos.(pi * nodes_coordinate[1,:]) .* sin.(pi * nodes_coordinate[2,:])
f[4*((1:nodes_N) .- 1) .+ 3] = 4*pi^5*sin.(pi * nodes_coordinate[1,:]) .* cos.(pi * nodes_coordinate[2,:])
f[4*((1:nodes_N) .- 1) .+ 4] = 4*pi^6*cos.(pi * nodes_coordinate[1,:]) .* cos.(pi * nodes_coordinate[2,:])

F = M * f

F[4*(nodes_boundary.-1).+1] = zeros(nodes_boundary_N)

boundary condition

A = copy(K)
tgv = 1e30
A[4*(nodes_boundary.-1).+1, 4*(nodes_boundary.-1).+1] += tgv * sparse(Matrix{Float64}(I, nodes_boundary_N, nodes_boundary_N))

linear system and explicit solution

u = A \ F
u_exact = zeros(Float64, 4*nodes_N)
u_exact[4*((1:nodes_N) .- 1) .+ 1] = sin.(pi * nodes_coordinate[1,:]) .* sin.(pi * nodes_coordinate[2,:])
u_exact[4*((1:nodes_N) .- 1) .+ 2] = pi * cos.(pi * nodes_coordinate[1,:]) .* sin.(pi * nodes_coordinate[2,:])
u_exact[4*((1:nodes_N) .- 1) .+ 3] = pi * sin.(pi * nodes_coordinate[1,:]) .* cos.(pi * nodes_coordinate[2,:])
u_exact[4*((1:nodes_N) .- 1) .+ 4] = pi^2 * cos.(pi * nodes_coordinate[1,:]) .* cos.(pi * nodes_coordinate[2,:])

error

err = u - u_exact

println("L2 error : ", sqrt(err' * M * err))
println("H2 error : ", sqrt(err' * K * err))

L2 error : 1.34151525591616e-6
H2 error : 0.013034241233629965

export to vtk

cells = [MeshCell(VTKCellTypes.VTK_QUAD, elements_int[1:4, i]) for i = 1:elements_int_N]


points_x = nodes_coordinate[1, :]
points_y = nodes_coordinate[2, :]
vtkfile = vtk_grid("ex4-output", points_x, points_y, cells)

vtkfile["u", VTKPointData()] = u[4*((1:nodes_N) .- 1) .+ 1]
vtkfile["ux", VTKPointData()] = u[4*((1:nodes_N) .- 1) .+ 2]
vtkfile["uy", VTKPointData()] = u[4*((1:nodes_N) .- 1) .+ 3]
vtkfile["uxy", VTKPointData()] = u[4*((1:nodes_N) .- 1) .+ 4]
outfiles = vtk_save(vtkfile)