Hello everyone,
Except some little projects with Fenics in my past I still trying to learn and work with FENCIS, whenever I can. Currently I tried to understand higher-order elements in FENICS better and found it at best confusing.
To illustrate my point I have developed this small WORKING example:
from mpi4py import MPI
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.fem import functionspace, Function
from dolfinx.io import XDMFFile
from basix.ufl import element
import numpy as np
element_type = "Lagrange"
element_dof = 1 # <-- This is what I will change later
mesh = create_unit_square(MPI.COMM_WORLD, 96, 96, CellType.triangle)
elem = element(element_type, mesh.basix_cell(), element_dof)
V = functionspace(mesh, elem)
u = Function(V)
# Zero u
u.x.array[:] = 0.0
# Interpolate some random function
u.interpolate(lambda x: 0.63 + 0.02 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()
from pathlib import Path
output = Path("/home/fenics/shared/") # this is important for jupyter notebooks
# Output file
file = XDMFFile(MPI.COMM_WORLD, str(output /f "test-{element_dof}.xdmf"), "w")
file.write_mesh(mesh)
file.write_function(u, 0)
file.close()
From my little understanding changing element_dof to something other then 1 (i.e. 2, 3, etc) should work with the same code, because more “interpolation points” are added, BUT doing so results in the following error:
RuntimeError Traceback (most recent call last)
/fenics_tests.ipynb Cell 10 line 3
29 file = XDMFFile(MPI.COMM_WORLD, str(output / f "test-{element_dof}.xdmf"), "w")
30 file.write_mesh(mesh)
---> 31 file.write_function(u, 0)
32 file.close()
File /usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/utils.py:235, in XDMFFile.write_function(self, u, t, mesh_xpath)
219 def write_function(self, u: Function, t: float = 0.0, mesh_xpath="/Xdmf/Domain/Grid[@GridType='Uniform'][1]"):
220 """Write function to file for a given time.
221
222 Note:
ref='/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/io/utils.py:0'>0</a>;32m (...)
233
234 """
--> 235 super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?
I studied issues like:
But still do I not know what to change.
Besides a workaround in code, which is degree independent, it would by nice to get an understanding how FENICS handles higher-order elements so I can better understand why this problem occur.
I would also appreciate if someone could show my how to debug such an issue, because debugging FENICS related issues, still seems like a mystery to me.
Thanks for your time!