Order for finite elements when using Gridap and Gmsh

mikeliux · January 23, 2025, 2:20am

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

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.
Please confirm if I should only import elements of order 1 from Gmsh into Gridap
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?
Can Gridap handle curved geometries like Gmsh does?

stevengj · January 23, 2025, 12:46pm

I think you are confusing two things here:

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.
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 · January 23, 2025, 3:00pm

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.

mikeliux · January 23, 2025, 6:33pm

Thank you very much for your quick reply.
I think you point to a very important distinction: curver elements versus higher order trial functions.
Higher order meshes can be used to represent curved boundaries, and I thought a higher order mesh would always imply a higher order polynomial for the trial function. Apparently, this is a very common approach called Isoparametric Formulation. Here are two resources describing it:

This distinction, helps us answer Q1.

stevengj · January 23, 2025, 6:52pm

Usually there is no practical purpose in using a higher-mesh unless you are also using higher-order polynomial element functions, because the rate of convergence is limited by the slowest-converging term. ~~However, this is not mathematically prohibited.~~ Also, the degrees don’t necessarily match (sometimes the basis functions are somewhat higher degree than the mesh IIRC), and it’s quite reasonable to use higher-order basis functions with first-order meshes in cases where the interfaces are flat.

Correction: if you use first-order basis functions (affine/planar) with a second (or higher) order mesh (curved boundaries), you have a problem: the basis functions are no longer continuous across boundaries (as they would be with first-order basis functions and a first-order mesh).

mikeliux · January 23, 2025, 9:50pm

Thank @Domenico_Lahaye you for your quick and comprehensive reply!

Some brief comments on your answers (I hope to update this answer another day):

Q1: when I set elementOrder = 2, GMSH automatically draws curved boundaries as seen in the original question. I don’t know what happens when I use interpolationOrder = 1 and elementOrder = 2 . I will try it later and update this reply.

Q2: There is an open issue in GridapGmsh about higher order mesh compatibility
In the code provided in my question, when I set elementOrder = 2 and import that mesh using GridapGmsh, for some reason the mesh is imported, but then these things happen::

The program raises an error when trying to export the model to a *.vtk for Paraview.
It raises an error when trying to make a wireframe plot of the of domain \Omega in Julia.
If I plot the solution in Julia, I get a very strange shape (I may add it later)
I think this can be reproduced by others with the code in the question.

So I suppose the answer to Q2 is that this is not implemented.

Q3: By visualization, I mean if I can plot the markers at all the nodes in Julia.

Q4: Gridap sees the mesh. But will it “see” a mesh with curved elements?

Best regards.

mikeliux · January 25, 2025, 9:12pm

Below is the updated code where I tried to answer you question, by importing a mesh of order 2 from Gmsh and solving the PDE in Gridap.

Code snippet 0:

#=
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 = 2
reffe_order = 2
mesh_file = "my_quartercicle$mesh_elements_order.msh"
println(repeat("=", 55))
println("Running program with mesh of order $mesh_elements_order and Reffe order $reffe_order")
println(repeat("=", 55))
# Create a DiscreteModel from the GMSH mesh file
model = GmshDiscreteModel(mesh_file)
writevtk(model,"model_my_quartercicle_order $mesh_elements_order")


# Step 2: Set up a finite element space for H^1 functions
reffe = ReferenceFE(lagrangian, Float64, reffe_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*reffe_order)

# Neumann boundary conditions
degree = reffe_order * 2
Γ_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 $mesh_elements_order and Reffe of order $reffe_order");
uh = solve(solver,op)
κ = 5.0
heat_flux = -κ *  ∇(uh)
# Step 6: Write the solution to a VTK file
println("Writing solution to a Paraview vtk file");
writevtk(Ω, "solution_pde01_mesh_order_$mesh_elements_order _reffe_order_$reffe_order", cellfields = ["uh" => uh, "heat_flux" => heat_flux])
   
# Step 7: Plot the solution in Julia (optional)

# Plot the solution using GridapMakie
println("\nPlotting solution with GridapMakie")
figure1 = Figure()
axis = Axis(
    figure1[1, 1],
    xlabel="x", ylabel="y",
    title="Solution plotted with plot(Ω, uh)\nMesh file: $mesh_file,  Reffe order: $reffe_order",
    aspect=DataAspect())
plt = plot!(Ω, uh) # renders a strange plot when mesh_elements_order = 2
try
    # mesh_elements_order != 1 raises assertion error in GridapMakie at conversions.jl
    wireframe!(Ω, color=:black, linewidth=2)
catch ex
    println("$(typeof(ex)) found with mesh of order $mesh_elements_order using 'wireframe!'. Ignoring error") 
end
colorbar = Colorbar(figure1[1, 2], plt, label="Temperature (°C)")
save("fig1_solution_mesh_order_$mesh_elements_order reffe_order_$reffe_order.png", figure1)
figure1 # display


# Plot the mesh nodes using GridapMakie
println("\nPlotting meshg nodes with GridapMakie")
figure2 = Figure()
axis = Axis(
    figure2[1, 1],
    xlabel="x", ylabel="y",
    title="Mesh nodes plotted with scatter!(axis, Ω)\nMesh file: $mesh_file,  Reffe order: $reffe_order",
    aspect=DataAspect())
scatter!(axis, Ω)
save("fig2_mesh_points_mesh_order_$mesh_elements_order reffe_order_$reffe_order.png", figure2)
figure2


# When mesh_elements_order = 2, evaluate! does not work at some nodes.
# Get all nodes in mesh Ω, evaluate uh, and track which return an error
println("\nEvaluating function and at nodes and creating a plot of valid/invalid nodes")
valid_coords = []
invalid_coords = []
Ω_node_coords = Gridap.Geometry.get_node_coordinates(Ω)
cache = Gridap.Fields.return_cache(uh, Ω_node_coords)
for coord in Ω_node_coords
    try
        val = evaluate!(cache, uh, coord)
        push!(valid_coords, coord)
    catch ex
        push!(invalid_coords, coord)        
        println("evaluate! return an $(typeof(ex)) with message: $(ex.msg)")
    end
end

as_matrix(v::Vector) = reshape([col[i] for col in v, i in 1:2], :, 2)
valid_coords = as_matrix(valid_coords)
invalid_coords = as_matrix(invalid_coords)

# Plot valid and invalid node coordinates
figure3 = Figure()
axis = Axis(
    figure3[1, 1],
    xlabel="x", ylabel="y",
    title="Result of evaluate!(cache, uh, coord)\nMesh file:$mesh_file,  Reffe order: $reffe_order",
    aspect=DataAspect())
scatter!(axis, valid_coords[:,1], valid_coords[:,2],color=:blue, label="valid nodes ($(size(valid_coords,1)))")
scatter!(axis, invalid_coords[:,1], invalid_coords[:,2],color=:red, label="invalid nodes ($(size(invalid_coords,1)))")
axislegend(axis)
save("fig3_mesh_points_classified_mesh_order_$mesh_elements_order reffe_order_$reffe_order.png", figure3)
figure3

This is the result of running the code snippet
0:
Code snippet 1:

=======================================================
Running program with mesh of order 2 and Reffe order 2
=======================================================
Info    : Reading 'my_quartercicle2.msh'...
Info    : 7 entities
Info    : 39 nodes
Info    : 24 elements
Info    : Done reading 'my_quartercicle2.msh'
Number of cells in Ω: 14
Solving PDE with mesh elements of order 2 and Reffe of order 2
Writing solution to a Paraview vtk file

Plotting solution with GridapMakie
AssertionError found with mesh of order 2 using 'wireframe!'. Ignoring error

Plotting meshg nodes with GridapMakie

Evaluating function and at nodes and creating a plot of valid/invalid nodes
evaluate! return an AssertionError with message: Point (0.0, 0.8333333333328944) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.1950903225681118, 0.980785280293434) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.2440169358561554, 0.2440169358555612) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.3436652596429693, 0.7207388309057189) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.7349294964951394, 0.3489659475084919) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.1220084679280777, 0.1220084679277806) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.1220084679280777, 0.2886751345940318) is not inside any active cell
evaluate! return an AssertionError with message: Point (0.2743320108720906, 0.3808075327917608) is not inside any active cell

To summarize, this happens when I import a mesh with triangular elements of second order from Gmsh into Gridap (mesh_elements_order in the given code):

GridapGmsh imports the mesh (see code snippet 0)
GridapGmsh can export the solution to vtk, and it contains the additional nodes (see Figure 0).
I was unable to plot the curved boundary with curved elements, neither using Paraview nor using Gridap (see Figure 0 - nodes along the arc segments are connected by straight lines).
GridapMakie raises an error trying to draw the wireframe of the domain \Omega (See code snippet 1)
GridapMakie makes an abnormal plot of the solution (see Figure 1)
GridapMakie can plot the node coordinates with a scatter plot (see Figure 2).
Gridap is unable to evaluate the solution at some nodes (se Figure 3).

Domenico_Lahaye · January 26, 2025, 9:18pm

Thank you for this update. Much appreciated.

I am not sure at all.

I wonder, however, whether the curved elements that GMSH generates can be visualized prior to any PDE being solved. I never verified whether the GMSH GUI properly renders the curved boundaries. I am not sure whether this a good starting point for further investigations.

I do not know whether GMSH marks curved boundary in any way in the *.msh file and whether one could learn from how Gridap (GridapGmsh) reads these elements.

mikeliux · January 27, 2025, 2:03am

Yes, GMSH does render curved elements on curved boundaries when using elements of order > 1. Please see the images with the green meshes which I posted in the original question. These were created prior to solving the PDE.

Domenico_Lahaye · January 27, 2025, 8:22am

Wonderful to have this cleared up.

The figures in green do not involve any Julia coding at all.

Can the same curved elements be viewed after reading and writing the mesh in Julia? This question decomposes into two subquestions:

Q1: Can the curved elements be identified (separated from non-curved elements) in a loop over all elements? (alternatively, how does GMSH label curved elements?)

Q2: What syntax should be used to write curved elements to a VTK file. I currently fail to see functionality for iso-parametric elements in the documentation of writeVTK.jl (see Unstructured grid formats · WriteVTK.jl). Do I overlook something?

Formulated as above, the issue in independent of how the PDE on the iso-parametric mesh is solved.

Topic		Replies	Views
GridapGmsh on gmsh.jl object General Usage fem , gridap , gmsh	2	621	January 25, 2025
GridapGmsh error loading .msh file when line and surface elements are used General Usage gridap , gmsh	0	329	February 2, 2023
Dirichlet boundary condition in gmsh New to Julia question , gmsh	4	379	January 8, 2025
How to create JSON files for finite element meshes? Numerics question	0	631	April 8, 2021
GMSH Mesh Expert to Julia using koehlerson/gmsh.jl General Usage question , package , gmsh	4	1255	March 30, 2021

Order for finite elements when using Gridap and Gmsh

Related topics