Order for finite elements when using Gridap and Gmsh

[mikeliux]

I am experimenting with Gridap package for solving partial differential equations (PDEs) in Julia. I am creating my meshes with GMSH, using the API provided by the Gmsh package in Julia. I save them as *.msh files, and then import them with GridapGmsh package. Finally, I solve the PDE using Gridap package and export the results as *.vtk files.

I am using the following package versions in julia 1.11.3:

  [705231aa] Gmsh v0.3.1
⌃ [56d4f2e9] Gridap v0.18.8
  [3025c34a] GridapGmsh v0.7.2

In Gmsh, the function gmsh.model.mesh.setOrder(order) is said to “Change the order of the elements in the mesh of the current model to order” . So, using the API for julia, I could specify second order elements as follows:

#=
Sample code to create a quarter of circle 2D domain and mesh it with GMSH using the API in Julia.
Change the element_order 
=#
import Gmsh: gmsh
gmsh.initialize()
gmsh.model.add("Simple mesh with Gmsh API in Julia")

mesh_size = 0.4
element_order = 2 # Modify this to change order of the finite elements.

# Define points of a quarter of circle
p1 = gmsh.model.geo.addPoint(0, 0, 0, mesh_size)
p2 = gmsh.model.geo.addPoint(0, 1, 0, mesh_size)
p3 = gmsh.model.geo.addPoint(1, 0, 0, mesh_size)

# borders of the triangle
bottom = gmsh.model.geo.add_line(p1, p2) # bottom border
arc1 = gmsh.model.geo.addCircleArc(p2, p1, p3)  # cicle 1 (start=p2, center=p1, end=p3)
left = gmsh.model.geo.add_line(p3, p1) # left border

gmsh.model.geo.synchronize()

# Combine all borders into a curved loop
quarter_circle = gmsh.model.geo.add_curve_loop([bottom, arc1, left])
gmsh.model.geo.synchronize()

quarter_circle_surface = gmsh.model.geo.addPlaneSurface([quarter_circle])
gmsh.model.geo.synchronize()


# Define tags for the boundaries using Gmsh interface in Julia
gmsh.model.addPhysicalGroup(1, [arc1])
gmsh.model.setPhysicalName(1, 1, "arc1")

gmsh.model.addPhysicalGroup(1, [bottom])
gmsh.model.setPhysicalName(1, 2, "bottom")

gmsh.model.addPhysicalGroup(1, [left])
gmsh.model.setPhysicalName(1, 3, "left")

gmsh.model.addPhysicalGroup(2, [quarter_circle_surface])
gmsh.model.setPhysicalName(2, 4, "quarter_circle_surface")

# Create mesh (order = 1 by default)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)
gmsh.model.geo.synchronize()

# Change order of mesh elements
gmsh.model.mesh.setOrder(element_order)
gmsh.model.geo.synchronize()
gmsh.write("my_quartercircle_order"*string(element_order)*".msh")
elementOrder = 1

# Set mesh faces visible, mesh points invisible, set geometrical elements all visible with description
gmsh.option.setNumber("Mesh.SurfaceFaces", 1) # Make mesh faces visible
gmsh.option.setNumber("Mesh.Points", 1)      # Make mesh points visible
gmsh.option.setNumber("Mesh.NumSubEdges", 6)   # Subdivisions of mesh elements to use while rendering (does not affect model)
gmsh.option.setNumber("Geometry.PointNumbers", 1) # Show geometrical elements with descriptions
gmsh.option.setNumber("Geometry.CurveNumbers", 1) 
gmsh.option.setNumber("Geometry.SurfaceNumbers", 1)

gmsh.fltk.run()
gmsh.finalize()

As we increase the order of the element from 1 to 3, more points are added and the curved boundaries are more accurately represented.

Mesh examples in GMSH with different order of mesh elements

Mesh file “my_quartercircle_order1.msh”

Mesh file “my_quartercircle_order2.msh”

Mesh file “my_quartercircle_order3.msh”

Importing meshes with GridapGmsh
Notice that GridapGmsh can import meshes with elements of order 1 or 2, but not greater, so my_quartercicle_order3.msh cannot be imported.

Actually, as far as I understand, the input to gridap should be a mesh with element order 1, and importing my_quartercicle_order2.msh in the code below raises some issues.

Next, when I define the test functions, can assign an order for the polynomials of the finite elements discretization in the functions ReferenceFE and TestFESpace. Can I use polynomials of any degree here, regardless of the order of the mesh elements?

#=
Solve Heat a transfer equation in 2D : ∇²T = 0
=#

using Gridap
using GridapGmsh
using GridapMakie, CairoMakie
using FileIO
 
# Step 1: Load the GMSH mesh
mesh_elements_order = 1 # try setting this equal 2 for an unexpected behaviour.
order_str = string(mesh_elements_order)
mesh_file = "my_quartercicle"*order_str*".msh"

# Create a DiscreteModel from the GMSH mesh file
model = GmshDiscreteModel(mesh_file)
writevtk(model,"model_my_quartercicle_order"*order_str)


# Step 2: Set up a finite element space for H^1 functions
order = 2 
reffe = ReferenceFE(lagrangian,Float64,order)
V = TestFESpace(
    model,
    reffe;conformity=:H1,
    dirichlet_tags = ["bottom", "left"] 
    )

# Step 3: Define the weak form for the Laplace equation
Ω = Triangulation(model)
dΩ = Measure(Ω,2*order)

# Neumann boundary conditions
degree = order + 1
Γ_arc = BoundaryTriangulation(model, tags=["arc1"])
dΓ_arc = Measure(Γ_arc, degree)


# Bilinear form for Laplace operator
a(u, v) = ∫( ∇(v) ⋅ ∇(u) )dΩ
# Right-hand side (homogeneous)
h_arc(x) = -10.0
l(v) = ∫( v*h_arc )dΓ_arc

println("Number of cells in Ω: ", num_cells(Ω))

# Step 4: Define the Dirichlet boundary condition values in a verctor
T_bottom(x) = 40.0
T_left(x) = 20.0
U = TrialFESpace(V, [T_left, T_bottom])

# Step 5: Solve the system
op = AffineFEOperator(a, l, U, V)
ls = LUSolver()
solver = LinearFESolver(ls)
println("Solving PDE with mesh elements of order "*order_str*"and Reffe of order "*string(order));
uh = solve(solver,op)
κ = 5.0
heat_flux = -κ *  ∇(uh)
# Step 6: Write the solution to a VTK file (works only with meshes of order 1)
if mesh_elements_order==1
    println("Writing solution to a Paraview vtk file");
    writevtk(Ω, "solution_pde01", cellfields = ["uh" => uh, "heat_flux" => heat_flux])
end

# Step 7: Plot the solution in Julia (optional)

# Plot the values
figure = Figure()
axis = Axis(
    figure[1, 1],
    xlabel="x", ylabel="y",
    title="Mesh file: $mesh_file,  Reffe order: $order",
    aspect=DataAspect())
plt = plot!(Ω, uh)
if mesh_elements_order ==1
    # order 2 raises an error
    wireframe!(Ω, color=:black, linewidth=2)
end
colorbar = Colorbar(figure[1, 2], plt, label="Temperature (°C)")
save("fig_solution_mesh_order"*order_str*"reffe_order"*string(order)*".png", figure)
figure # display

The number of nodes changes when the model is imported with GmshDiscreteModel:

Element order = 1 in GMSH

Info    : Reading 'my_quartercicle1.msh'...  
Info    : 7 entities
Info    : 13 nodes
Info    : 24 elements
Info    : Done reading 'my_quartercicle1.msh'

Element order = 2 in GMSH

Info    : Reading 'my_quartercicle2.msh'...
Info    : 7 entities
Info    : 39 nodes
Info    : 24 elements
Info    : Done reading 'my_quartercicle2.msh'

For the rest of this question, refer to this mesh of order 1: ‘my_quartercicle1.msh’

Suppose I have created a mesh with elements of order 1 and saved it as “my_quartercicle_order1.msh”. Now I import this mesh into Julia using the function GmshDiscreteModel from the package GridapGmsh.

So, using elements of order 1 in the mesh, I can change the order of the ReferenceFE and solve the same PDE. As I increase the order of the polynomials in Gridap, the accuracy of the results improves as expected eventhough I start with the same “my_quartercicle_order1.msh” file.

I can plot the mesh using ReferenceFE elements of order 1 and 2, and it does not show curved edges.

Questions

1. If I understand correctly, the order of the FEM elements defines the order of the polynomials, so I am not sure what happens when I import a mesh with elements of order 1 from Gmsh and then use ReferenceFE of order 2.
2. Please confirm if I should only import elements of order 1 from Gmsh into Gridap
3. I suppose using ReferenceFE of order 2 creates more nodes in the mesh for gridap. How can I plot the inner nodes in my FEM discretization when I set the ReferenceFE order to 2?
4. Can Gridap handle curved geometries like Gmsh does?

[stevengj]

@mikeliux

If I understand correctly, the order of the FEM elements defines the order of the polynomials, so I am not sure what happens when I import a mesh with elements of order 1 from Gmsh and then use ReferenceFE of order 2.

I think you are confusing two things here:

1. The order of the mesh (what you specify in Gmsh): this indicates whether your geometric edges/faces are straight (order 1) or curved (higher order). The latter are useful to conform the mesh to curved boundaries.
2. The order of your FEM basis functions (what you specify in Gridap): this is the degree of the polynomials used to describe/interpolate unknown functions within each geometric element. This is often much higher than the order of your mesh.

These two are related only in the sense that if you want to solve PDE with high-order accuracy (i.e. convergence/errors scale as a high power of the mesh resolution), and you have some curved interface where the solution has a boundary condition or discontinuity, then you may need both a high-order mesh (in order to accurately describe the interface) and high-order basis functions (in order to accurately describe the solution).

[Domenico_Lahaye]

Many thanks for a wonderful post!

Q1: In Gridap, the order of FEM elements defines the polynomial degree of the approximation or the second item in. the list of @stevengj. One could verify the size of the linear system the that Gridap generates to verify this claim. Gridap has its own mesh generator for rectangular geometries. On these simple geometries, Gridap is able to generate second order approximations. I therefore imagine that Gridap has functionality to add second order information (i.e., midpoints of edges) to some extend. I do not know how complete this functionality this. I do not know either if one is able to override this functionality with second information provided by the mesh file (generated by GMSH) on input (in case one wishes to do so).

GMSH provides separate options for second order basis functions (elementOrder) and for representation of curved boundaries (interpolationOrder). See e.g. the GMSH tutorial x6 tutorials/julia/x6.jl · gmsh_4_13_1 · gmsh / gmsh · GitLab

Q2: Not sure. What happens if you generate a mesh with elementOrder = 2 as input for Gridap?

Q3: Yes, elementOrder = 2 generates a longer list of nodes compared to elementOrder = 1. It suffices to look into the *.msh files generated. Not sure what you means by visualization. I imagine that we could explore this by looking into the VTK files that Gridap generates.

Q4: Gridap obviously sees only the mesh, not the geometry. See also Q1.