# Order of element on a 1D sutructure **URL:** \ **Category:** General\ **Created:** [October 11, 2024, 11:55am UTC](https://fenicsproject.discourse.group/t/order-of-element-on-a-1d-sutructure/16024 "2024-10-11T11:55:02Z")\ **Posts on this page:** 2\ **Page:** 1
**Author:** ![pierremrtt](https://avatars.discourse-cdn.com/v4/letter/p/e0b2c6/32.png) [@pierremrtt](https://fenicsproject.discourse.group/u/pierremrtt)\ **Post date:** [October 11, 2024, 11:55am UTC](https://fenicsproject.discourse.group/t/order-of-element-on-a-1d-sutructure/16024/1 "2024-10-11T11:55:02Z")
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 : ```auto 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 : ```auto --------------------------------------------------------------------------- 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 ---
**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\ **Post date:** [October 13, 2024, 5:47am UTC](https://fenicsproject.discourse.group/t/order-of-element-on-a-1d-sutructure/16024/2 "2024-10-13T05:47:49Z")
Running with the current `main` I get ```auto Traceback (most recent call last): File "/file.py", line 94, in 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.