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()
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.
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:
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)
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.
##-------------------- 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)
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).
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?)
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)
