I’m experiencing an issue with my FEM analysis using Dolfinx (version 0.10.0) and PETSc where the computed displacement field is identically zero—even though nonzero forcing is applied (gravity and an attachment force). I suspect that either the applied Dirichlet BC on the floor is over-constraining the domain or there is an assembly/forcing error in the variational forms.
Details:
- I apply a gravitational body force and an additional attachment force (as a Neumann term on facets marked with a specific marker).
- The floor boundary is enforced by a Dirichlet condition (zero displacement on all DOFs located via a geometric query on the minimum y‑coordinate).
- I have verified that the forcing terms are assembled and the boundary markers are set correctly.
- The assembled system is solved using PETSc’s GMRES/KSP with an ILU preconditioner.
- Upon solution, the displacement function is set using the underlying vector’s array property.
- For example, my code and solving is as follows:
import petsc4py
petsc4py.init()
import os
import json
import numpy as np
import dolfinx
from dolfinx import mesh, fem
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import basix
from basix.ufl import blocked_element
from dolfinx.io import XDMFFile
import h5py
import requests
from db import get_db_connection, get_latest_iteration
import dolfinx.fem.petsc
from dolfinx.fem.petsc import set_bc
# Define SI units: (all coordinates in meters)
BOUNDING_BOX = {"x": [-0.05, 0.05], "y": [0.0, 0.13], "z": [-0.13, 0.13]}
ATTACHMENT_POSITION = np.array([0.0, 0.13, -0.11]) # in m
TOLERANCE = 0.005 # 5 mm in m
# URL for the Mesh Service
MESH_SERVICE_URL = "http://mesh_service:8003"
SHARED_ASSETS_DIR = "/app/src/assets/"
class FEMAnalysis:
def __init__(self, updated_density=None):
os.environ["OMP_NUM_THREADS"] = "1"
self.xdmf_path, self.h5_path = self._fetch_mesh_paths()
self.mesh, self.disp_fun = self._load_mesh()
self.V = self._setup_function_space()
# Set up the boundary markers for Dirichlet BCs.
self.boundary_markers = self._setup_boundaries()
self.updated_density = updated_density
print(f"📌 Using Mesh:\n XDMF: {self.xdmf_path}\n H5: {self.h5_path}", flush=True)
def _fetch_mesh_paths(self):
print("🔵 Fetching latest mesh iteration from Mesh Service...", flush=True)
response = requests.get(f"{MESH_SERVICE_URL}/latest_mesh_iteration")
latest_iteration = response.json().get("iteration")
response = requests.get(f"{MESH_SERVICE_URL}/get_mesh_info/{latest_iteration}")
mesh_data = response.json()
xdmf_path = mesh_data.get("xdmf_path")
h5_path = mesh_data.get("h5_path")
if not xdmf_path or not h5_path:
raise RuntimeError("❌ Mesh paths missing!")
print(f"✅ Using mesh files: {xdmf_path}, {h5_path}", flush=True)
return xdmf_path, h5_path
def _load_mesh(self):
print(f"🔵 Loading mesh from {self.xdmf_path}...", flush=True)
with XDMFFile(MPI.COMM_WORLD, self.xdmf_path, "r") as xdmf:
loaded_mesh = xdmf.read_mesh(name="mesh")
# Convert coordinates from mm to m (if needed; here assume they are in mm)
print("🔵 Converting mesh coordinates from mm to m.", flush=True)
loaded_mesh.geometry.x[:] /= 1000.0
print(f"🔵 Loading displacement function from {self.h5_path}...", flush=True)
basix_element = basix.create_element(
family=basix.ElementFamily.P,
celltype=basix.CellType.tetrahedron,
degree=1,
lagrange_variant=basix.LagrangeVariant.equispaced
)
ufl_element = basix.ufl.wrap_element(basix_element)
vector_element = blocked_element(ufl_element, (3,))
V_temp = fem.functionspace(loaded_mesh, vector_element)
disp_fun = fem.Function(V_temp)
comm = MPI.COMM_WORLD
if comm.rank == 0:
with h5py.File(self.h5_path, "r") as h5_file:
func_group = h5_file.get("Function")
f_group = func_group.get("f")
dataset_keys = list(f_group.keys())
disp_array = np.array(f_group[dataset_keys[0]]).flatten()
else:
disp_array = None
disp_array = comm.bcast(disp_array, root=0)
if disp_array.shape[0] != len(disp_fun.x.array):
raise RuntimeError("❌ Mismatch in displacement data and function space DOFs!")
disp_fun.x.array[:] = disp_array.reshape(-1)
print("✅ Displacement function loaded.", flush=True)
return loaded_mesh, disp_fun
def _setup_function_space(self):
print("🔵 Setting up function space...", flush=True)
basix_element = basix.create_element(
family=basix.ElementFamily.P,
celltype=basix.CellType.tetrahedron,
degree=1,
lagrange_variant=basix.LagrangeVariant.equispaced
)
ufl_element = basix.ufl.wrap_element(basix_element)
self.ulf_element = ufl_element
V = fem.functionspace(self.mesh, ufl_element)
print(f"✅ Function space created with {V.dofmap.index_map.size_global} DOFs.", flush=True)
return V
def _setup_boundaries(self):
print("🔵 Setting up boundary conditions...", flush=True)
# Locate the floor facets by finding facets where y is minimal.
floor_facets = dolfinx.mesh.locate_entities_boundary(
self.mesh, self.mesh.topology.dim - 1,
lambda x: np.isclose(x[1], np.min(self.mesh.geometry.x[:, 1]), atol=1e-6)
)
print(f"👉 Floor boundary: {len(floor_facets)} facets found.", flush=True)
markers = np.full(len(floor_facets), 1, dtype=np.int32)
meshtags = dolfinx.mesh.meshtags(self.mesh, self.mesh.topology.dim - 1, floor_facets, markers)
return meshtags
def _apply_forces(self, u, v):
print("🔵 Applying forces to FEM model...", flush=True)
rho = 7850 # kg/m^3
g = 9.81 # m/s²
gravity_force = fem.Constant(self.mesh, PETSc.ScalarType(ufl.as_vector([0, -rho * g, 0])))
body_force = ufl.dot(gravity_force, v) * ufl.dx
# Apply attachment force as a Neumann term.
attachment_point = ATTACHMENT_POSITION # already in m
tolerance = TOLERANCE # 5 mm in m
print("🔍 Checking attachment point proximity...", flush=True)
all_points = self.mesh.geometry.x
close_points = [p for p in all_points if np.linalg.norm(p - attachment_point) < tolerance]
if len(close_points) == 0:
raise ValueError("❌ No mesh points found near the attachment point!")
print(f"✅ Found {len(close_points)} points near the attachment.", flush=True)
attachment_dofs = fem.locate_dofs_geometrical(
self.V, lambda x: np.linalg.norm(x.T - attachment_point, axis=1) < tolerance
)
# Here we add the attachment force as a Neumann term.
if len(attachment_dofs) > 0:
force_vector = np.array([0, -800, 0], dtype=np.float64)
force_constant = fem.Constant(self.mesh, PETSc.ScalarType(ufl.as_vector(force_vector)))
# Assume that facets marked with 2 represent the attachment region.
ds = ufl.Measure("ds", domain=self.mesh, subdomain_data=self.boundary_markers)
F_attach = ufl.dot(force_constant, v) * ds(2)
print(f"✅ Applying attachment force {force_vector} over marker 2.", flush=True)
else:
print("⚠️ No attachment DOFs found. Skipping attachment force.", flush=True)
F_attach = 0
print("✅ Forces successfully applied.", flush=True)
return body_force + F_attach
def run_analysis(self, updated_density=None):
try:
iteration = get_latest_iteration() + 1
print(f"📌 Running FEM analysis for iteration {iteration}.", flush=True)
# Check that the mesh has cells.
num_cells = self.mesh.topology.index_map(self.mesh.topology.dim).size_global
if num_cells == 0:
raise ValueError("❌ No cells in mesh! Mesh generation failed.")
nodes = self.mesh.geometry.x
print(f"[DEBUG] Mesh has {nodes.shape[0]} nodes.")
print(f"[DEBUG] Mesh Y-range: min={np.min(nodes[:,1])}, max={np.max(nodes[:,1])}")
print(f"[DEBUG] Material properties: E={210e6}, nu={0.3}, mu={210e6/(2*(1.3))}, lambda={210e6*0.3/(1.3*0.4)}", flush=True)
# Create vector-valued function space.
vector_element = blocked_element(self.ulf_element, (3,))
self.V = fem.functionspace(self.mesh, vector_element)
# Optionally (re)create scalar function space for density.
P1 = fem.functionspace(self.mesh, ("CG", 1))
density_function = fem.Function(P1)
# Set default material properties.
E_base, nu = 210e6, 0.3
mu_base = E_base / (2.0 * (1.0 + nu))
lmbda_base = E_base * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
if updated_density is not None:
density_array = np.array(updated_density, dtype=np.float64)
if density_array.shape[0] != self.V.dofmap.index_map.size_global:
raise RuntimeError("❌ Density field size mismatch!")
density_function.x.array[:] = density_array
E_scaled = density_function * E_base
mu_scaled = E_scaled / (2.0 * (1.0 + nu))
lmbda_scaled = E_scaled * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
else:
mu_scaled = mu_base
lmbda_scaled = lmbda_base
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u):
return lmbda_scaled * ufl.tr(epsilon(u)) * ufl.Identity(3) + 2 * mu_scaled * epsilon(u)
U_trial = ufl.TrialFunction(self.V)
V_test = ufl.TestFunction(self.V)
a = ufl.inner(sigma(U_trial), epsilon(V_test)) * ufl.dx
L = self._apply_forces(U_trial, V_test)
print("🔵 Assembling the system...", flush=True)
# Assemble the matrix and vector.
A = dolfinx.fem.petsc.assemble_matrix(fem.form(a))
A.assemble()
b = dolfinx.fem.petsc.assemble_vector(fem.form(L))
# Apply the floor (Dirichlet) boundary condition.
zero_vector = fem.Constant(self.mesh, np.array([0.0, 0.0, 0.0], dtype=np.float64))
floor_dofs = fem.locate_dofs_geometrical(
self.V,
lambda x: np.isclose(x[1], np.min(self.mesh.geometry.x[:, 1]), atol=1e-6)
)
bc = fem.dirichletbc(zero_vector, floor_dofs, self.V)
# dolfinx.fem.petsc.set_bc(b, [bc])
if not hasattr(b, "array_w"):
b.array_w = b.array # assign the writable array property to array_w
dolfinx.fem.petsc.set_bc(b, [bc])
# Solve the system.
# Solve the system.
u_sol = fem.Function(self.V)
ksp = PETSc.KSP().create(MPI.COMM_WORLD)
ksp.setOperators(A)
ksp.setType(PETSc.KSP.Type.GMRES)
ksp.getPC().setType(PETSc.PC.Type.ILU)
x = A.createVecRight()
ksp.solve(b, x)
# Instead of using x.array, we use get_local()
sol_array = x.array.copy()
u_sol.x.array[:] = sol_array
print("✅ FEM Solver succeeded!", flush=True)
# Write results.
os.makedirs(SHARED_ASSETS_DIR, exist_ok=True)
xdmf_filename = os.path.join(SHARED_ASSETS_DIR, f"fem_results_{iteration}.xdmf")
h5_filename = os.path.join(SHARED_ASSETS_DIR, f"fem_results_{iteration}.h5")
json_filename = os.path.join(SHARED_ASSETS_DIR, f"fem_results_{iteration}.json")
with XDMFFile(MPI.COMM_WORLD, xdmf_filename, "w") as xdmf:
xdmf.write_mesh(self.mesh)
xdmf.write_function(u_sol, t=0.0)
displacement_data = {"iteration": iteration, "displacement": u_sol.x.array.tolist()}
with open(json_filename, "w") as f:
json.dump(displacement_data, f, indent=4)
print(f"✅ FEM results stored (iteration {iteration}).", flush=True)
return {"iteration": iteration, "xdmf_file": xdmf_filename, "h5_file": h5_filename, "json_file": json_filename}
except Exception as e:
print(f"❌ FEM Analysis Error: {e}", flush=True)
raise RuntimeError(e)
After solving, I output the displacement array to a JSON file, but it contains only zeros. I have also verified that the forces (gravity and attachment) are nonzero and that the boundary markers (both for the floor and the attachment) are set as expected.
here are the docker logs:
🔵 Fetching latest mesh iteration from Mesh Service...
🔵 Fetching latest mesh iteration from Mesh Service...
✅ Using mesh files: /app/src/assets/convex_hull_1_mesh.xdmf, /app/src/assets/convex_hull_1_disp.h5
✅ Using mesh files: /app/src/assets/convex_hull_1_mesh.xdmf, /app/src/assets/convex_hull_1_disp.h5
🔵 Loading mesh from /app/src/assets/convex_hull_1_mesh.xdmf...
🔵 Converting mesh coordinates from mm to m.
🔵 Loading mesh from /app/src/assets/convex_hull_1_mesh.xdmf...
🔵 Loading displacement function from /app/src/assets/convex_hull_1_disp.h5...
🔵 Converting mesh coordinates from mm to m.
🔵 Loading displacement function from /app/src/assets/convex_hull_1_disp.h5...
✅ Displacement function loaded.
✅ Displacement function loaded.
🔵 Setting up function space...
🔵 Setting up function space...
✅ Function space created with 126 DOFs.
✅ Function space created with 126 DOFs.
🔵 Setting up boundary conditions...
🔵 Setting up boundary conditions...
👉 Floor boundary: 32 facets found.
👉 Floor boundary: 32 facets found.
📌 Using Mesh:
XDMF: /app/src/assets/convex_hull_1_mesh.xdmf
H5: /app/src/assets/convex_hull_1_disp.h5
📌 Using Mesh:
XDMF: /app/src/assets/convex_hull_1_mesh.xdmf
H5: /app/src/assets/convex_hull_1_disp.h5
📌 Running FEM analysis for iteration 1.
[DEBUG] Mesh has 126 nodes.
[DEBUG] Mesh Y-range: min=0.1275, max=0.13
[DEBUG] Material properties: E=210000000.0, nu=0.3, mu=80769230.76923077, lambda=121153846.15384614
🔵 Applying forces to FEM model...
📌 Running FEM analysis for iteration 1.
🔍 Checking attachment point proximity...
[DEBUG] Mesh has 126 nodes.
[DEBUG] Mesh Y-range: min=0.1275, max=0.13
[DEBUG] Material properties: E=210000000.0, nu=0.3, mu=80769230.76923077, lambda=121153846.15384614
✅ Found 126 points near the attachment.
🔵 Applying forces to FEM model...
✅ Applying attachment force [ 0. -800. 0.] over marker 2.
🔍 Checking attachment point proximity...
✅ Forces successfully applied.
🔵 Assembling the system...
✅ Found 126 points near the attachment.
✅ Applying attachment force [ 0. -800. 0.] over marker 2.
✅ Forces successfully applied.
🔵 Assembling the system...
✅ FEM Solver succeeded!
✅ FEM Solver succeeded!
✅ FEM results stored (iteration 1).
✅ FEM results stored (iteration 1).
✅ FEM results stored in database under iteration 1.
INFO: 172.26.0.1:59682 - "GET /run_fem HTTP/1.1" 200 OK
Questions & Request:
- Could the zero displacement result be due to over-constraining with the floor Dirichlet BC? How can I verify whether these BCs are too restrictive?
- Is it possible that there is an issue with the assembly of the forcing term? For example, are there known pitfalls in assembling Neumann terms for the attachment in Dolfinx that might lead to an effectively zero right-hand side?
- Are there any recommended debugging strategies or additional printing/debugging options in Dolfinx/PETSc that could help diagnose why the computed displacements are zero?