Writing XDMF file leads to inf & nan values

dolfinx
1

Hi,
I’m trying to solve an iterative FEM problem, but I am getting weird results exporting the results to ParaView. The problem needs to be solved for incremental loading, but for certain timesteps, the XDMF file appears to get inf & nan values. Here is an example result of a single step, where the part on the right looks OK, but there are regions with nan & inf values.

Screenshot from 2024-03-05 17-28-43
Screenshot from 2024-03-05 17-28-43886×641 143 KB

When I export to PDFs using pyvista, the result looks much better:
image
image703×572 90.7 KB

I’m still poorly understanding the whole FunctionSpace → Function → Interpolate process when I already have a mesh and the solution values, so I’m guessing I am doing something incorrect there.
Here is the code:

import os
import numpy as np
import ufl
import dolfinx
from dolfinx import mesh, fem, io
from dolfinx.io import XDMFFile
from dolfinx.plot import vtk_mesh
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
from petsc4py import PETSc
import pyvista as pv

print(f"dolfinx version: {dolfinx.__version__}")

# Read the mesh
mesh_file = 'meshes/L.xdmf'
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, mesh_file, "r") as xdmf_file:
    mesh_domain = xdmf_file.read_mesh(name="Grid")

tdim = mesh_domain.topology.dim - 1

# Vector function space for displacement field
element = ufl.VectorElement("CG", mesh_domain.ufl_cell(), 1)
V = fem.FunctionSpace(mesh_domain, element)

# Define trial and test functions
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Define material properties (Young's modulus, Poisson's ratio)
E = 210e9  # Example: 210 GPa for steel
nu = 0.3
lambda_ = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))

# Strain and stress
def epsilon(u):
    return ufl.sym(ufl.grad(u))

def sigma(u):
    return lambda_ * ufl.div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)

#### Dirichlet BC (0 displacement on top boundary)
# Define function for displacement boundary conditions
uD = fem.Function(V)
uD.interpolate(lambda x: np.zeros((3, x.shape[1])))     # zero displacements

def boundary_points(x):   # return top boundary points
    return np.isclose(x[1], max(x[1]))

boundary_dofs = fem.locate_dofs_geometrical(V, boundary_points)
bc = fem.dirichletbc(value=uD, dofs=boundary_dofs)

##### Neumann BC (stress/force on the right boundary)
def right_boundary(x):   # return right boundary points
    return np.isclose(x[0], max(x[0]))

load_dofs = mesh.locate_entities_boundary(mesh_domain, tdim, right_boundary)

bound_tag = mesh.MeshTags(load_dofs)

max_load = -6  # Maximum load magnitude
num_increments = 10  # Number of increments
increment = max_load / num_increments  # Load increment

g = fem.Constant(mesh_domain, PETSc.ScalarType((0, 0, 0)))

ds = ufl.Measure('ds', subdomain_data=bound_tag)

# Weak form
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
f = fem.Constant(mesh_domain, PETSc.ScalarType((0, 0, 0)))
L = ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

# Solve the problem for several steps, and save solutions
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
solutions = []

for step in range(1, num_increments + 1):
    # Update bc applied load value:
    g.value = PETSc.ScalarType((0, step * increment, 0))

    # Solve the problem
    uh = problem.solve()

    # Store the solution for analysis or post-processing
    solutions.append(uh.copy())

# Save results for ParaView
plot_space = fem.FunctionSpace(mesh_domain, ("CG", 1))
u_y = fem.Function(plot_space)

with XDMFFile(mesh_domain.comm, "linEla_debug/linOut.xdmf", "w") as write_file:
    write_file.write_mesh(mesh_domain)
    for step, uh in enumerate(solutions):
        # Plot result to xdmffile
        u_y.name = "$u_y$"
        u_y.interpolate(lambda x: uh.x.array[1::3])  # [1::3] selects y displacements from 3D vector
        write_file.write_function(u_y, step)

        # 2D y displacement pdf plots using pyvista instead
        u_topology, u_cell_types, u_geometry = vtk_mesh(V)
        u_grid = pv.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
        u_grid.point_data["u"] = uh.x.array[1::3]
        u_grid.set_active_scalars("u")
        u_plotter = pv.Plotter()
        u_plotter.add_mesh(u_grid, show_edges=True)
        u_plotter.view_xy()
        u_plotter.save_graphic(f'linEla_debug/linOutCont_y_{step}.pdf')

Many thanks!

2

Hi.

You did not supply the file for the mesh you are using. However, you can try to use fun.sub(i).collapse() to obtain the function component. In your test case, it would be uh.sub(1).collapse() to get the y component, instead of lambda x: uh.x.array[1::3]

BTW, welcome to the FEniCS community!

2 Likes
3

Hi Jesus, thank you for the reply!
Using fun.sub(i).collapse() works! I don’t understand why this is necessary, could you elaborate on that?

Also, a follow-up question if that is okay: in my real case I want to plot von Mises stresses as well, which I compute as a 1D numpy array. However, I get the same plotting behavior as in my first question. What is the correct way to interpolate this for plotting, since the numpy array can’t be collapsed like uh?

Here is the code snippet, let me know if I should provide the full code instead:

        sig = fem.Function(Vsig)
        sig.name = "$\sigma$"
        sig.interpolate(lambda x: von_mis_np)
        write_file.write_function(sig, step)
4

As stated here by @dokken for function spaces,sub(i) and collapse() gives you a function space over the ith element of the mixed element, which can be seen in your case as the “mixed element” between CG1 spaces (one per component). For an already defined fem.Function, you could do the same (see this link) to access its components if available. By using uh.x.array[1::3] you have to be aware of the dof ordering on FEniCS.

Regarding your second question, since you now have uh, you could obtain the stress by using an interpolation into a DG0 space. I would do the following

Tensorplot_space = fem.FunctionSpace(
mesh_domain, ("DG", 0, (3, 3)))
sig_h = fem.Function(Tensorplot_space)
sig = fem.Expression(ufl.as_ufl(sigma(uh)),
                 Tensorplot_space.element.interpolation_points())
sig_h.interpolate(sig)
sig_h.name = "$\sigma$"

write_file.write_function(sig_h, step)

Maybe there is a better (faster) way of doing that.

Cheers.

3 Likes
5

Thanks a lot, this is very useful.