# Error and I can't follow up with new fixes

**URL:** <https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092>\
**Category:** Errors\
**Created:** [October 21, 2024, 3:52pm UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092 "2024-10-21T15:52:10Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![sohailahmedx](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sohailahmedx/32/7697_2.png) [@sohailahmedx](https://fenicsproject.discourse.group/u/sohailahmedx)\
**Post date:** [October 21, 2024, 3:52pm UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092/1 "2024-10-21T15:52:10Z")

</div>

I’m trying to learn python and get used to FenicsX, But I’m encountering a lot of issues, I’m trying to simulate 3d Wave Equation in a Steel Plate

```auto
import numpy as np
import dolfinx
from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import matplotlib.pyplot as plt

# --- 1. Material Properties ---
density = 7850 # kg/m^3
E = 210e9 # Young's modulus in Pascals

# --- 2. Mesh Generation ---
length = 1.0 # Length of the plate in meters
width = 0.5 # Width of the plate in meters
height = 0.01 # Thickness of the plate in meters

# Creating the mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 10, 2, cell_type=dolfinx.mesh.CellType.hexahedron)
mesh.geometry.x[:, 0] *= length # Scale in x-direction
mesh.geometry.x[:, 1] *= width # Scale in y-direction
mesh.geometry.x[:, 2] *= height # Scale in z-direction

# --- 3. Function Space ---
element = basix.ufl.element("Lagrange", basix.CellType.hexahedron, 1, shape=(3,))
V = dolfinx.fem.functionspace(mesh, element)

# --- 4. Boundary Conditions ---
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], length)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0
bc = fem.dirichletbc(u_bc, fem.locate_dofs_geometrical(V, boundary))

# --- 5. Time Parameters ---
T = 0.5 # Total time in seconds
num_steps = 100 # Number of time steps
dt = T / num_steps # Time step size

# --- 6. Wave Speed ---
c = np.sqrt(E / density)
c_squared = fem.Constant(mesh, PETSc.ScalarType(c**2).real)
dt_sq = fem.Constant(mesh, PETSc.ScalarType(dt**2).real)

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=PETSc.ScalarType)
    values[:, 0] = np.sin(np.pi * x[:, 0]) * np.sin(2 * np.pi * x[:, 1]) 
    return values

u_n1.interpolate(initial_condition)

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, 0.0])) 

# --- 10. Variational Formulation ---
dx = ufl.Measure("dx", domain=mesh)

# Define the variational forms
a = ufl.inner(u, v) * dx + dt_sq * c_squared * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = dt_sq * ufl.inner(f, v) * dx + ufl.inner(2 * u_n - u_n1, v) * dx

# --- 11. Time-Stepping Loop ---
u = fem.Function(V)
for n in range(num_steps):
    petsc_options = {"ksp_type": "cg", "pc_type": "ilu"} 
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
    problem.solve()

    if n % 10 == 0:
        with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"displacement_3d_{n}.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(u, n*dt)

    # Update solution 
    u_n1.x.array[:] = u_n.x.array
    u_n.x.array[:] = u.x.array

# --- 12. Plotting ---
# Get coordinates of the mesh nodes
x = mesh.geometry.x

# Reshape the displacement array (extract the x-component)
u_array = u.x.array.reshape(-1, 3)[:, 0]

# Get the number of nodes in x and y directions
num_nodes_x = 20  
num_nodes_y = 10  
u_array = u_array.reshape(num_nodes_y, num_nodes_x)

num_nodes_x = mesh.topology.index_map(0).size_local
num_nodes_y = mesh.topology.index_map(1).size_local

# Create the plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
X, Y = np.meshgrid(np.linspace(0, length, num_nodes_x), np.linspace(0, width, num_nodes_y))
ax.plot_surface(X, Y, u_array, cmap='viridis')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Displacement (m)')
ax.set_title('Displacement of the Steel Plate at Final Time')
plt.show()

```

And I get this error  
TypeError Traceback (most recent call last)  
Cell In[5], line 79  
77 for n in range(num\_steps):  
78 petsc\_options = {“ksp\_type”: “cg”, “pc\_type”: “ilu”}  
—\> 79 problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc\_options=petsc\_options)  
80 problem.solve()  
82 if n % 10 == 0:

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/petsc.py:782, in LinearProblem. **init** (self, a, L, bcs, u, petsc\_options, form\_compiler\_options, jit\_options)  
744 def **init** (  
745 self,  
746 a: ufl.Form,  
(…)  
752 jit\_options: typing.Optional[dict] = None,  
753 ):  
754 “”“Initialize solver for a linear variational problem.  
755  
756 Args:  
(…)  
780 “mumps”})  
781 “””  
 → 782 self.\_a = \_create\_form(  
783 a,  
784 dtype=PETSc.ScalarType,  
785 form\_compiler\_options=form\_compiler\_options,  
786 jit\_options=jit\_options,  
787 )  
788 self.\_A = create\_matrix(self.\_a)  
789 self.\_L = \_create\_form(  
790 L,  
791 dtype=PETSc.ScalarType,  
792 form\_compiler\_options=form\_compiler\_options,  
793 jit\_options=jit\_options,  
794 )

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/forms.py:337, in form(form, dtype, form\_compiler\_options, jit\_options, entity\_maps)  
334 else:  
335 return form  
 → 337 return \_create\_form(form)

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/forms.py:331, in form..\_create\_form(form)  
329 return None  
330 else:  
 → 331 return \_form(form)  
332 elif isinstance(form, collections.abc.Iterable):  
333 return list(map(lambda sub\_form: \_create\_form(sub\_form), form))

File /usr/local/dolfinx-complex/lib/python3.12/dist-packages/dolfinx/fem/forms.py:306, in form..\_form(form)  
303 else:  
304 \_entity\_maps = {msh.\_cpp\_object: emap for (msh, emap) in entity\_maps.items()}  
 → 306 f = ftype(  
307 module.ffi.cast(“uintptr\_t”, module.ffi.addressof(ufcx\_form)),  
308 V,  
309 coeffs,  
310 constants,  
311 subdomains,  
312 \_entity\_maps,  
313 mesh,  
314 )  
315 return Form(f, ufcx\_form, code, module)

TypeError: **init** (): incompatible function arguments. The following argument types are supported:  
1. **init** (self, spaces: collections.abc.Sequence[dolfinx.cpp.fem.FunctionSpace\_float64], integrals: collections.abc.Mapping[dolfinx.cpp.fem.IntegralType, collections.abc.Sequence[tuple[int, int, ndarray[dtype=int32, shape=(_), order=‘C’, writable=False], ndarray[dtype=int8, shape=(_), order=‘C’, writable=False]]]], coefficients: collections.abc.Sequence[dolfinx.cpp.fem.Function\_complex128], constants: collections.abc.Sequence[dolfinx.cpp.fem.Constant\_complex128], need\_permutation\_data: bool, entity\_maps: collections.abc.Mapping[dolfinx.cpp.mesh.Mesh\_float64, ndarray[dtype=int32, shape=(_), order=‘C’, writable=False]], mesh: dolfinx.cpp.mesh.Mesh\_float64 | None) → None  
2. **init** (self, form: int, spaces: collections.abc.Sequence[dolfinx.cpp.fem.FunctionSpace\_float64], coefficients: collections.abc.Sequence[dolfinx.cpp.fem.Function\_complex128], constants: collections.abc.Sequence[dolfinx.cpp.fem.Constant\_complex128], subdomains: collections.abc.Mapping[dolfinx.cpp.fem.IntegralType, collections.abc.Sequence[tuple[int, ndarray[dtype=int32, shape=(_), order=‘C’, writable=False]]]], entity\_maps: collections.abc.Mapping[dolfinx.cpp.mesh.Mesh\_float64, ndarray[dtype=int32, shape=(\*), order=‘C’, writable=False]], mesh: dolfinx.cpp.mesh.Mesh\_float64 | None) → None

Invoked with types: dolfinx.cpp.fem.Form\_complex128, \_CDataBase, list, list, list, dict, dict, dolfinx.cpp.mesh.Mesh\_float64

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 21, 2024, 7:18pm UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092/2 "2024-10-21T19:18:28Z")

</div>

I would change all dtype definitions to `dolfinx.default_scalar_type`, i.e.

```python
c_squared = fem.Constant(mesh, dolfinx.default_scalar_type(c**2))
dt_sq = fem.Constant(mesh, dolfinx.default_scalar_type(dt**2))

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=dolfinx.default_scalar_type)
    values[:, 0] = np.sin(np.pi * x[:, 0]) * np.sin(2 * np.pi * x[:, 1]) 
    return values

u_n1.interpolate(initial_condition)

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, 0.0], dtype=dolfinx.default_scalar_type)) 

```

that will bring you all the way down to postprocessing  
where

> [@sohailahmedx](#):
>
> ```auto
> # Get the number of nodes in x and y directions
> num_nodes_x = 20  
> num_nodes_y = 10  
> u_array = u_array.reshape(num_nodes_y, num_nodes_x)
> 
> ```

this doesn’t make sense to me.

---

<div class="post-metadata">

**Author:** ![sohailahmedx](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sohailahmedx/32/7697_2.png) [@sohailahmedx](https://fenicsproject.discourse.group/u/sohailahmedx)\
**Post date:** [October 22, 2024, 7:04am UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092/3 "2024-10-22T07:04:08Z")

</div>

For the displacement data (u\_array), and reorganizes it into a grid (num\_nodes\_y rows by num\_nodes\_x columns) that matches the structure of the mesh and especially for the plotting. Please correct me if I’m wrong

---

<div class="post-metadata">

**Author:** ![sohailahmedx](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sohailahmedx/32/7697_2.png) [@sohailahmedx](https://fenicsproject.discourse.group/u/sohailahmedx)\
**Post date:** [October 22, 2024, 8:41pm UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092/6 "2024-10-22T20:41:10Z")

</div>

I have made a few changes but, There’s no displacement values

```auto
import numpy as np
import dolfinx
from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import matplotlib.pyplot as plt

# --- 1. Material Properties ---
density = 7850 # kg/m^3
E = 210e9 # Young's modulus in Pascals

# --- 2. Mesh Generation ---
length = 1.0 # Length of the plate in meters
width = 0.5 # Width of the plate in meters
height = 0.01 # Thickness of the plate in meters

# Creating the mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 10, 2, cell_type=dolfinx.mesh.CellType.tetrahedron)
mesh.geometry.x[:, 0] *= length # Scale in x-direction
mesh.geometry.x[:, 1] *= width # Scale in y-direction
mesh.geometry.x[:, 2] *= height # Scale in z-direction

# Create edge entities (dimension 1)
mesh.topology.create_entities(1)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "my_mesh.xdmf", "w") as file:
    file.write_mesh(mesh)

# --- 3. Function Space ---
element = basix.ufl.element("Lagrange", basix.CellType.tetrahedron, 1, shape=(3,))
V = dolfinx.fem.functionspace(mesh, element)

# --- 4. Boundary Conditions ---
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], length)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0
bc = fem.dirichletbc(u_bc, fem.locate_dofs_geometrical(V, boundary))

# --- 5. Time Parameters ---
T = 1.0 # Total time in seconds
num_steps = 100 # Number of time steps
dt = T / num_steps # Time step size

# --- 6. Wave Speed ---
c = np.sqrt(E / density)
c_squared = fem.Constant(mesh, dolfinx.default_scalar_type(c**2))
dt_sq = fem.Constant(mesh, dolfinx.default_scalar_type(dt**2))

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=dolfinx.default_scalar_type)
    center_x = length / 2
    center_y = width / 2
    sigma = 0.1  
    amplitude = 0.01  
    values[:, 0] = amplitude * np.exp(-((x[:, 0] - center_x)**2 + (x[:, 1] - center_y)**2) / (2 * sigma**2))
    return values

u_n1.interpolate(initial_condition)

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, -9.81 * density], dtype=dolfinx.default_scalar_type)) 

# --- 10. Variational Formulation ---
dx = ufl.Measure("dx", domain=mesh)

a = ufl.inner(u, v) * dx + dt_sq * c_squared * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = dt_sq * ufl.inner(f, v) * dx + ufl.inner(2 * u_n - u_n1, v) * dx

# --- 11. Time-Stepping Loop ---
u = fem.Function(V)
for n in range(num_steps):
    petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}

    # Initialize LinearProblem.
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
    
    problem.solve()

    # Check displacement values every few steps to debug.
    if n % 10 == 0:
        print(f"Step {n}, Displacement max: {np.max(u.x.array)}, min: {np.min(u.x.array)}")

        with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"displacement_3d_{n}.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(u, n*dt)

    # Update solution 
    u_n1.x.array[:] = u_n.x.array
    u_n.x.array[:] = u.x.array

# --- 12. Plotting ---
x = mesh.geometry.x
u_array = u.x.array.reshape(-1, 3)[:, 0]

num_nodes_x = mesh.topology.index_map(0).size_local  
total_nodes = u_array.size  
num_nodes_y = total_nodes // num_nodes_x  

if total_nodes % num_nodes_x != 0:
    raise ValueError("Total nodes are not divisible by number of nodes in x direction.")

X, Y = np.meshgrid(np.linspace(0, length, num_nodes_x), np.linspace(0, width, num_nodes_y))
Z_array = u_array.reshape(num_nodes_y, num_nodes_x)  

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z_array, cmap='viridis')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Displacement (m)')
ax.set_title('Displacement of the Steel Plate at Final Time')
plt.show()

```

```auto
Step 0, Displacement max: (nan+nanj), min: (nan+nanj)
Step 10, Displacement max: (nan+nanj), min: (nan+nanj)
Step 20, Displacement max: (nan+nanj), min: (nan+nanj)
Step 30, Displacement max: (nan+nanj), min: (nan+nanj)
Step 40, Displacement max: (nan+nanj), min: (nan+nanj)
Step 50, Displacement max: (nan+nanj), min: (nan+nanj)
Step 60, Displacement max: (nan+nanj), min: (nan+nanj)
Step 70, Displacement max: (nan+nanj), min: (nan+nanj)
Step 80, Displacement max: (nan+nanj), min: (nan+nanj)
Step 90, Displacement max: (nan+nanj), min: (nan+nanj)

```

Any help?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 22, 2024, 9:42pm UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092/7 "2024-10-22T21:42:21Z")

</div>

> [@sohailahmedx](#):
>
> ```auto
> petsc_options = {"ksp_type": "preonly", "pc_type": "lu"}
> 
> # Initialize LinearProblem.
> problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
>     
> problem.solve()
> 
> ```

Try adding  
` "pc_factor_mat_solver_type": "mumps"`  
and  
`"ksp_error_if_not_converged": True`  
to `petsc_options` get some output from the petsc ksp solver (and use a more stable LU solver

---

<div class="post-metadata">

**Author:** ![sohailahmedx](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sohailahmedx/32/7697_2.png) [@sohailahmedx](https://fenicsproject.discourse.group/u/sohailahmedx)\
**Post date:** [October 22, 2024, 10:53pm UTC](https://fenicsproject.discourse.group/t/error-and-i-cant-follow-up-with-new-fixes/16092/8 "2024-10-22T22:53:17Z")

</div>

Tried that

```auto
import numpy as np
import dolfinx
from dolfinx import fem, mesh
import ufl
from mpi4py import MPI
import dolfinx.fem.petsc
from petsc4py import PETSc
import basix
import matplotlib.pyplot as plt

# --- 1. Material Properties ---
density = 7850 # kg/m^3
E = 210e9 # Young's modulus in Pascals

# --- 2. Mesh Generation ---
length = 1.0 # Length of the plate in meters
width = 0.5 # Width of the plate in meters
height = 0.01 # Thickness of the plate in meters

# Creating the mesh
mesh = dolfinx.mesh.create_unit_cube(MPI.COMM_WORLD, 20, 10, 2, cell_type=dolfinx.mesh.CellType.tetrahedron)
mesh.geometry.x[:, 0] *= length # Scale in x-direction
mesh.geometry.x[:, 1] *= width # Scale in y-direction
mesh.geometry.x[:, 2] *= height # Scale in z-direction

# Create edge entities (dimension 1)
mesh.topology.create_entities(1)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "my_mesh.xdmf", "w") as file:
    file.write_mesh(mesh)

# --- 3. Function Space ---
element = basix.ufl.element("Lagrange", basix.CellType.tetrahedron, 1, shape=(3,))
V = dolfinx.fem.functionspace(mesh, element)

# --- 4. Boundary Conditions ---
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], length)

u_bc = fem.Function(V)
u_bc.x.array[:] = 0.0 # Ensure real values for boundary conditions
bc = fem.dirichletbc(u_bc, fem.locate_dofs_geometrical(V, boundary))

print(f"Boundary condition values: {u_bc.x.array}")

# --- 5. Time Parameters ---
T = 1.0 # Total time
num_steps = 700 # Increase the number of steps for smaller time step size
dt = T / num_steps # Smaller time step

# --- 6. Wave Speed ---
c = np.sqrt(E / density)
c_squared = fem.Constant(mesh, dolfinx.default_scalar_type(c**2))
dt_sq = fem.Constant(mesh, dolfinx.default_scalar_type(dt**2))
c = np.sqrt(E / density)
print(f"Wave speed: {c}")

# --- 7. Trial and Test Functions ---
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# --- 8. Initial Conditions ---
u_n = fem.Function(V)
u_n1 = fem.Function(V)

def initial_condition(x):
    values = np.zeros((x.shape[0], 3), dtype=PETSc.ScalarType)
    center_x = length / 2
    center_y = width / 2
    sigma = 0.1
    amplitude = 0.01
    values[:, 2] = amplitude * np.exp(-((x[:, 0] - center_x)**2 + (x[:, 1] - center_y)**2) / (2 * sigma**2)) # Z-direction only
    return values

u_n.interpolate(initial_condition)
u_n1.interpolate(initial_condition)

print(f"Initial max displacement (u_n): {np.max(u_n.x.array)}, min displacement (u_n): {np.min(u_n.x.array)}")

# --- 9. External Force ---
f = fem.Constant(mesh, np.array([0.0, 0.0, -9.81 * density], dtype=dolfinx.default_scalar_type)) 

# --- 10. Variational Formulation ---
dx = ufl.Measure("dx", domain=mesh)

a = ufl.inner(u, v) * dx + dt_sq * c_squared * ufl.inner(ufl.grad(u), ufl.grad(v)) * dx
L = ufl.inner(2 * u_n - u_n1, v) * dx # Temporarily remove the external force for debugging

print(f"L form: {L}")

# --- 11. Time-Stepping Loop ---
u = fem.Function(V)
for n in range(num_steps):
    # Define PETSc solver options with MUMPS and error reporting
    petsc_options = {
        "ksp_type": "preonly", # Direct solver
        "pc_type": "lu", # LU decomposition
        "pc_factor_mat_solver_type": "mumps", # Use MUMPS solver
        "ksp_error_if_not_converged": True # Raise error if the solver doesn't converge
    }

    # Initialize the LinearProblem.
    problem = dolfinx.fem.petsc.LinearProblem(a, L, u=u, bcs=[bc], petsc_options=petsc_options)
    
    # Solve the linear problem
    problem.solve()

    # Check displacement values every few steps to debug.
    if n % 10 == 0:
        print(f"Step {n}, Displacement max: {np.max(u.x.array)}, min: {np.min(u.x.array)}")

        # Write displacement field to file
        with dolfinx.io.XDMFFile(MPI.COMM_WORLD, f"displacement_3d_{n}.xdmf", "w") as file:
            file.write_mesh(mesh)
            file.write_function(u, n * dt)

    # Update solution for the next time step
    u_n.x.array[:] = np.real(u_n.x.array)
    u_n1.x.array[:] = np.real(u_n1.x.array)

print(f"Initial max displacement (u_n): {np.max(u_n.x.array)}, min displacement (u_n): {np.min(u_n.x.array)}")
print(f"Initial max displacement (u_n1): {np.max(u_n1.x.array)}, min displacement (u_n1): {np.min(u_n1.x.array)}")

# --- 12. Plotting ---
x = mesh.geometry.x
u_array = u.x.array.reshape(-1, 3)[:, 0]

num_nodes_x = mesh.topology.index_map(0).size_local  
total_nodes = u_array.size  
num_nodes_y = total_nodes // num_nodes_x  

print(f"Step {n}, Max displacement: {np.max(u.x.array)}, Min displacement: {np.min(u.x.array)}")

if total_nodes % num_nodes_x != 0:
    raise ValueError("Total nodes are not divisible by number of nodes in x direction.")

X, Y = np.meshgrid(np.linspace(0, length, num_nodes_x), np.linspace(0, width, num_nodes_y))
Z_array = u_array.reshape(num_nodes_y, num_nodes_x)  

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, Z_array, cmap='viridis')
ax.set_xlabel('X coordinate (m)')
ax.set_ylabel('Y coordinate (m)')
ax.set_zlabel('Displacement (m)')
ax.set_title('Displacement of the Steel Plate at Final Time')
plt.show()

```

Gives this and there’s nothing in the plot

 ![Screenshot 2024-10-23 010052](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/9ea022343fecec3cede7a0441dbe0f4b38ca9f76.png)

```auto
Boundary condition values: [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]
Wave speed: 5172.194153034851
Initial max displacement (u_n): (4.093e-320+nanj), min displacement (u_n): (4.093e-320+nanj)
L form: { conj(((v_0) : (({ A | A_{i_{58}} = 2 * f[i_{58}] }) + ({ A | A_{i_{59}} = -1 * f[i_{59}] })))) } * dx(<Mesh #19>[everywhere], {})
Step 0, Displacement max: (nan+nanj), min: (nan+nanj)
Step 10, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 20, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 30, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 40, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 50, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 60, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 70, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 80, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 90, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 100, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 110, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 120, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 130, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 140, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 150, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 160, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 170, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 180, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 190, Displacement max: (1.3635697704911794e+290+0j), min: 0j
Step 200, Displacement max: (1.3635697704911794e+290+0j), min: 0j

```