Order of element on a 1D sutructure

General
1

Hi community,

I try to create a simple case of a 1D structure where the Helmholtz equation is solved with bc. It works perfectly when I use a space function with Lagrangian polynomial of degree 1, but it doesn’t work with quadratic element, here is the code and the output error :

import numpy as np
import matplotlib.pyplot as plt

from mpi4py import MPI
from petsc4py import PETSc
from ufl import TrialFunction, TestFunction, dx, grad, inner, div, ds, as_vector, Measure, conj
from basix.ufl import element

import dolfinx
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, functionspace, assemble_vector
from dolfinx.fem.petsc import LinearProblem
from dolfinx.fem import petsc
from dolfinx.mesh import create_interval, locate_entities_boundary, meshtags
from dolfinx.io import XDMFFile

# Problem parameters
Lf = 1.0  # Domain length (1 meter)
rho_0 = 1.21  # Volumic masse of air (kg/m^3)
c = 343.0  # Speed of sound in the air (m/s)
omega = 2 * np.pi * 1000  # Angular frequency
k = omega / c  # Wave number
Z_n = rho_0 * c  # Acoustic impedance

# 1D mesh
nb_cells = 100
mesh = create_interval(MPI.COMM_WORLD, nb_cells, [0, Lf])
cells_mesh = mesh.ufl_cell()

# Space function
deg = 2
#Ve = element("Lagrange", mesh.ufl_cell(), deg)
V = functionspace(mesh, ("Lagrange", deg))

# Test and Trial functions
p = TrialFunction(V)
v = TestFunction(V)

# Variational formulation
a = inner(grad(p), grad(v)) * dx - k**2 * inner(p, v) * dx

# Detect the boundaries
def left_boundary(x):
    return np.isclose(x[0], 0)

def right_boundary(x):
    return np.isclose(x[0], Lf)

# Find the facets
facets_left = locate_entities_boundary(mesh, mesh.topology.dim - 1, left_boundary)
facets_right = locate_entities_boundary(mesh, mesh.topology.dim - 1, right_boundary)

# Create a MeshTag for the facets
facet_tags = meshtags(mesh, mesh.topology.dim - 1,
                      np.concatenate([facets_left, facets_right]),
                      np.array([0] * len(facets_left) + [1] * len(facets_right)))

# Define a Measure, surface integral on boundaries
ds_boundary = Measure("ds", domain=mesh, subdomain_data=facet_tags)

a += inner((1j * omega * rho_0 / Z_n) * p,  v) * ds_boundary(1)  

# Source term
f = omega**2 * rho_0 * 1.0
L = inner(f, v) * ds_boundary(0)

# Solving the linear system
problem = LinearProblem(a, L, petsc_options={"ksp_type": "preonly",
                                            "pc_type": "lu",
                                            "pc_factor_mat_solver_type": "mumps"})
                                                           
p_solution = problem.solve()

u_f = p_solution.x.array

x = np.linspace(0, Lf, 1000)
k = omega/c
p_ana = 1j*rho_0*c*omega*1.0*np.exp(1j*k*x)

plt.plot(x, p_ana.real, label = 'Analytical')
plt.plot(np.linspace(0,Lf, (nb_cells*deg) + 1), u_f.real, label = 'Numerical')
plt.legend()
plt.show()



# Sauvegarde et visualisation
if True:
    with XDMFFile(MPI.COMM_WORLD, "solution_helmholtz.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        #xdmf.write_function(p_solution) #This line is changed
        u1 = Function(V)
        u1.interpolate(p_solution)
        # Sauvegarde de la fonction projetée
        xdmf.write_function(u1)

The error :

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/usr/local/lib/python3.10/dist-packages/dolfinx/fem/function.py in interpolate(self, u0, cells0, cells1)
    458             # u is a Function or Expression (or pointer to one)
--> 459             _interpolate(u0)
    460         except TypeError:

3 frames
TypeError: interpolate(): incompatible function arguments. The following argument types are supported:
    1. interpolate(self, f: ndarray[dtype=complex128, shape=(*), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
    2. interpolate(self, f: ndarray[dtype=complex128, shape=(*, *), order='C', writable=False], cells: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
    3. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells0: ndarray[dtype=int32, shape=(*), order='C', writable=False], cells1: ndarray[dtype=int32, shape=(*), order='C', writable=False]) -> None
    4. interpolate(self, u: dolfinx.cpp.fem.Function_complex128, cells: ndarray[dtype=int32, shape=(*), order='C', writable=False], interpolation_data: dolfinx.cpp.geometry.PointOwnershipData_float64) -> None
    5. interpolate(self, e0: dolfinx.cpp.fem.Expression_complex128, cells0: ndarray[dtype=int32, order='C', writable=False], cells1: ndarray[dtype=int32, order='C', writable=False]) -> None

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

During handling of the above exception, another exception occurred:

AssertionError                            Traceback (most recent call last)
/usr/local/lib/python3.10/dist-packages/dolfinx/fem/function.py in interpolate(self, u0, cells0, cells1)
    460         except TypeError:
    461             # u0 is callable
--> 462             assert callable(u0)
    463             x = _cpp.fem.interpolation_coords(
    464                 self._V.element, self._V.mesh.geometry._cpp_object, cells0

AssertionError:

Thanks a lot

2

Running with the current main I get

Traceback (most recent call last):
  File "/file.py", line 94, in <module>
    xdmf.write_function(u1)
  File "lib/python3.12/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?

which may be more informative in telling you what to do to fix the error.