Export functionalities of Dolfinx 0.9 and troubleshooting with the interpolate function

dolfinx
1

Hello everyone,
I want to explore the export functionalities of dolfinx 0.9. Specifically, if higher order finite element functions can directly be exported without having to do any interpolations.

For this I constructed the following minimal working example:

from mpi4py import MPI
import dolfinx
import dolfinx.mesh
import dolfinx.fem
import dolfinx.io
import basix
import numpy as np
from ufl import *

# Create a mesh
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

# Define P2 function space
# Ve = basix.ufl.element("Lagrange", mesh.basix_cell(), degree=2)
V = dolfinx.fem.functionspace(mesh, ("CG", 2))

# Create and name a function
u = dolfinx.fem.Function(V)
u.name = "solution"
def f(mod, x):
    return sin(pi * x[0]) * sin(pi * x[1])

u.interpolate(f)
expr = dolfinx.fem.Expression(u, V.element.interpolation_points())

# Write mesh and function to XDMF
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "solution.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(u)

I get the following error :

Traceback (most recent call last):
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/dolfinx/fem/function.py", line 459, in interpolate
    _interpolate(u0)
    ~~~~~~~~~~~~^^^^
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/functools.py", line 934, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
           ~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/dolfinx/fem/function.py", line 440, in _interpolate
    self._cpp_object.interpolate(u0, cells0, cells1)  # type: ignore
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. interpolate(self, f: ndarray[dtype=complex128, writable=False, shape=(*), order='C'], cells: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    2. interpolate(self, f: ndarray[dtype=complex128, writable=False, shape=(*, *), order='C'], cells: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    3. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells0: ndarray[dtype=int32, writable=False, shape=(*), order='C'], cells1: ndarray[dtype=int32, writable=False, shape=(*), order='C']) -> None
    4. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells: ndarray[dtype=int32, writable=False, shape=(*), order='C'], interpolation_data: dolfinx.cpp.geometry.PointOwnershipData_float64) -> None
    5. interpolate(self, e0: dolfinx.cpp.fem.Expression_complex128, cells0: ndarray[dtype=int32, writable=False, order='C'], cells1: ndarray[dtype=int32, writable=False, order='C']) -> None

Invoked with types: dolfinx.cpp.fem.Function_complex128, function, ndarray, ndarray

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "src/new_dolfinx/P2Fields.py", line 25, in <module>
    u.interpolate(f)
    ~~~~~~~~~~~~~^^^
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/dolfinx/fem/function.py", line 466, in interpolate
    self._cpp_object.interpolate(np.asarray(u0(x), dtype=self.dtype), cells0)  # type: ignore
                                            ~~^^^
  File "/Users/anant/PhD/felics_simulations/src/new_dolfinx/P2Fields.py", line 24, in f
    return ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
           ~~~~~~~^^^^^^^^^^^^^^^
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/ufl/operators.py", line 628, in sin
    return _mathfunction(f, Sin)
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/ufl/operators.py", line 597, in _mathfunction
    f = as_ufl(f)
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/ufl/constantvalue.py", line 517, in as_ufl
    raise ValueError(
        f"Invalid type conversion: {expression} can not be converted to any UFL type."
    )
ValueError: Invalid type conversion: [2.82743339e+00 3.14159265e+00 3.14159265e+00 ... 1.57079633e-01
 1.57079633e-01 7.51366423e-34] can not be converted to any UFL type.

If I simply comment out :

def f(mod, x):
    return sin(pi * x[0]) * sin(pi * x[1])

u.interpolate(f)
expr = dolfinx.fem.Expression(u, V.element.interpolation_points())

and run the remaining code, then I get the error :

Traceback (most recent call last):
  File "src/new_dolfinx/P2Fields.py", line 37, in <module>
    xdmf.write_function(u)
    ~~~~~~~~~~~~~~~~~~~^^^
  File "/opt/anaconda3/envs/envt_name/lib/python3.13/site-packages/dolfinx/io/utils.py", line 252, in write_function
    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?

Would be highly grateful for any leads. Thanks and regards! :slight_smile:

2

The XDMF format isn’t particularly extensible for high order elements. Consider the demo dolfinx/python/demo/demo_interpolation-io.py at main · FEniCS/dolfinx · GitHub and VTXWriter in particular.

1 Like
3

Thank you very much for your response, @nate .

Do you know what I am doing wrong with the interpolate function?

4

I don’t see anything wrong with your code. It’s just that when you call

    xdmf.write_function(u)

XDMF files do not support high order (i.e. p > 1) Lagrange elements.

You can try interpolating into a lower order space (if that fits your needs) or use VTXWriter.