Read in quadratic mesh to represent CG2 elements

Hello everyone,

I am working with a sector of a cylinder mesh in 2D and I was wondering if reading a mesh with quadratic elements is possible. For example, I created a mesh in Gmsh (and converted it to XDMF using meshio) with quadratic triangles as shown in the figure.

The ideal case would be that I can read that mesh with dolfin and, when creating a quadratic continuous space the middle nodes were the same as the ones in the file. I tried to do that using the following code:

import dolfin as d
mesh = d.Mesh()
with d.XDMFFile('quadratic.xdmf') as infile:
    infile.read(mesh)
    
V = d.VectorFunctionSpace(mesh, 'CG', 2)
u = d.Function(V)

xdmf = d.XDMFFile(mesh.mpi_comm(), 'quadratic_fenics.xdmf')
xdmf.write_checkpoint(u, 'u', 0, d.XDMFFile.Encoding.HDF5, False)
xdmf.close()

but the result is this

I know I can load a linear mesh and FEniCS will interpolate the elements to create the middle nodes like this:

But I lose geometrical accuracy that will matter in the convergence of my problem.

I do not know if this is possible or not. I will appreciate your guidance.

Hi, high-order geometry is supported in dolfin-x only.

1 Like

To extend at @bleyerj’s reply:
Examples of how to read higher order Gmsh meshes into dolfin-X is covered in:
https://github.com/FEniCS/dolfinx/blob/master/python/demo/gmsh/demo_gmsh.py
and
https://github.com/FEniCS/dolfinx/blob/master/python/test/unit/io/test_xdmf_mesh.py#L98

2 Likes

Hello,
Thank you for your previous answers.
In my case I have a xdmf file with quadratic tetrahedron that I want to import in Dolfinx.
I defined a function of degree 2 for the displacements but I have problems when I use xdmf.write_function, I got the following error.

Please create a minimal example for your mesh.

The error above is due to the fact that DOLFINx is not able to to pull the degrees of freedom back to the reference element when interpolating the function you are writing to file to a CG-1 function space:

(You error should be easy to reproduce with the eval function as covered here: dolfinx/test_function.py at main · FEniCS/dolfinx · GitHub

An alternative workaround is to use the VTX writer, as it does not pull back values to the reference element; see for instance dolfinx/test_adios2.py at 114c268aab752b6669abc2718837116cec90dba6 · FEniCS/dolfinx · GitHub

Thank you very much for answering.
Regarding the mesh, it is a .inp mesh converted to .xdmf using meshio (the original mesh comes for a CT scan, elements are quadratic tetrahedron)

Regarding the function space I wrote

So I don’t completely understand what you said about the function interpolation to a CG-1 function space because I fixed a degree = 2. Thanks !

To be more specific, when you use the function XDMFFile.write_function it in turn calls compute_point_values, which computes an interpolation of your function onto the mesh nodes (10 nodes per tetrahedron).
This interpolation requires a pull_back, which causes the issue you describe above.

For me to be able to be of any further assistance, you need to supply a minimal code example which can reproduce your issue.

This is a minimal example where my issue appears.

##-------------------- Minimal example
import dolfinx
from mpi4py import MPI
import ufl
import dolfinx.io
import numpy as np

# Mesh construction
filename = "meshes/FFH6R_msh.xdmf"
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, "r" ) as file:
   mesh = file.read_mesh(name="Grid")

tdim=3 # problem in 3D
# Function space
element_res = ufl.FiniteElement("Lagrange", mesh.ufl_cell(), degree=2)
V_res = dolfinx.FunctionSpace(mesh, element_res)
res = dolfinx.Function(V_res)

# Values imported from the original file
fileresult = "FFH6R_B_590N_All_nodes_UVW.txt"
res_All_val = np.genfromtxt(fileresult, usecols=(3),  usemask=True)

res.vector.array = res_All_val

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "output/res1mini_x.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    res.name = "U (mm)"
    xdmf.write_function(res)

Thanks again

First of all, you need to either supply your xdmf/h5 file or an XML version of the xdmf file, such that one can redproduce your error message.

Secondly, I do not think you need the line:

to produce the error.

Another option would be to supply the .inp file and your meshio conversion script.

Hello,
I tried to use VTXwriter but I got 2 files (md.0 and data.0) that I cannot open with Paraview. Probably I did not use the code properly (please see below).

#-------------------- Minimal example

from mpi4py import MPI
import ufl
import dolfinx.io
import numpy as np
import os

Mesh construction

filename = “meshes/FFH6R_msh.xdmf”
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, filename, “r” ) as file:
mesh = file.read_mesh(name=“Grid”)

import dolfinx

Function space

element_res = ufl.FiniteElement(“Lagrange”, mesh.ufl_cell(), degree=2)
V_res = dolfinx.FunctionSpace(mesh, element_res)
res = dolfinx.Function(V_res)
res.interpolate(lambda x: (1*np.ones_like(x[0])))

filename_VTX = os.path.join(“output/”, “v.bp”)
f = dolfinx.cpp.io.VTXWriter(MPI.COMM_WORLD, filename_VTX, [res._cpp_object])

Set two cells to 0

for c in [0, 1]:
res.x.array[V_res.dofmap.cell_dofs(c)] = 1

res.x.scatter_forward()

Save twice and update geometry

for t in [0.1, 1]:
mesh.geometry.x[:, :2] += 0.1
f.write(t)

f.close()

Other possiblity to use XDMF but does not work with XDMFFile.write_function and the pull back

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, “output/res_cube1.xdmf”, “w”) as xdmf:

xdmf.write_mesh(mesh)

res.name = “U (mm)”

xdmf.write_function(res)

To open the files made with VTXWriter you must open the folder (ending with .BP) and not an individual file in the folder when using paraview

Thank you very much, it seems to work very well !

I have a second question regarding the capability of DOLFINx to handle p-version of finite element analysis (increasing Legendre polynomial shape functions).
Are we able to use FeniCs for 3D linear elasticity problem and Legendre polynomial shape functions of order p>4 considering the complexity to build the stiffness matrix in these conditions.
Thank you in advance.

Hello, we have experimental support for Legendre polynomial shape functions in FEniCSx.

You should be able to create these using the following snippet:

from ufl import FiniteElement
from dolfinx.fem import FunctionSpace

mesh = ...
element = FiniteElement("DG", mesh.ufl_cell(), 4, variant="legendre")
space = FunctionSpace(mesh, element)

Note that the element must be “DG” as these polynomials are not defined using point evaluations at the endpoints of the interval. (Unless I’ve misunderstood what you mean by Legendre polynomials?)

edit: @dokken’s correction. Also I’ve opened a PR in Basix to allow "Legendre" with a cap to work too as I imagine this will be a common typo (Convert variant names to lowercase to allow eg "GLL" to be used by mscroggs · Pull Request #392 · FEniCS/basix · GitHub)

1 Like

Note that it should be non-capitalized legendre, i.e.:

from ufl import FiniteElement
from dolfinx.fem import FunctionSpace

mesh = ...
element = FiniteElement("DG", mesh.ufl_cell(), 4, variant="legendre")
space = FunctionSpace(mesh, element)
1 Like

Hello,
I moved to the recent version of dolfinx (0.6.0) where I am trying to use shape functions based on Legendre polynomials of order p>4.
However it seems it is not possible anymore to use these polynomials (with ufl.FiniteElement - Form language — Unified Form Language (UFL) 2021.1.0 documentation), could you tell me more about that (if it is correct or not) ?
Thanks !

You would need to use the basix.create_element and the basix.ufl_wrapper to use this element with DOLFINx v0.6.0.
Consider the following minimal working example

import dolfinx
from mpi4py import MPI
import basix
import basix.ufl_wrapper
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
el = basix.create_element(basix.ElementFamily.P, basix.cell.string_to_type(mesh.ufl_cell().cellname()), 4, basix.LagrangeVariant.legendre, True)
element = basix.ufl_wrapper.BasixElement(el)
space = dolfinx.fem.FunctionSpace(mesh, element)

Here

el = basix.create_element(basix.ElementFamily.P, basix.cell.string_to_type(mesh.ufl_cell().cellname()), 4, basix.LagrangeVariant.legendre, True)

is equivalent to

where the last boolean tells DOLFINx that the element should be discontinuous.

Thank you very much, it works.
Otherwise I am not sure it is what I expected.
I have a 1D bar (in traction):
comm = MPI.COMM_WORLD
Nb_elem= 200
from dolfinx.mesh import create_interval
mesh = create_interval(comm, Nb_elem, [0,Lx])
And I want to implement a p-FE version it means that for each element I will get many dofs but 2 of them will be associated to the end points of the element
Defining the element as discontinuous make my current way of imposing boundary conditions (clamp at x=0 and traction at x=Lx not working :
element_u = basix.ufl_wrapper.BasixElement(element)
V_u = FunctionSpace(mesh, element_u)

Boundary definitions dofs

dofs_u_left = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], 0.0))
dofs_u_right = locate_dofs_geometrical(V_u, lambda x: np.isclose(x[0], Lx))

Indeed I get
RuntimeError: Cannot evaluate dof coordinates - this element does not have pointwise evaluation.

Please make a minimal reproducible example, as I’m missing several crucial definitions here to recreate the error message.