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
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
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.3415152766051647e-6
H2 error : 0.013034241233630735
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)
affichage de la solution