Hi all,
I’m trying to convert some FEniCS code to the current FEniCSx (based on dolfinx). I have a NumPy array representing spatially varying Young’s modulus values on a 2D mesh, and I want to assign it to a Functiondefined on a CG1 function space.
Here is what I used previously:
E_vals = np.load("E_field.npy")
assert E_vals.shape == (Ny+1, Nx+1)
ordering = df.dof_to_vertex_map(V_E)
E_func.vector()[:] = E_vals.flatten(order='C')[ordering]
E_func.vector().apply("insert")
Here is my current Fenicsx implementation that is giving wrong result:
from time import time
import numpy as np
import matplotlib.pyplot as plt
from mpi4py import MPI
from ufl import sym, grad, Identity, tr, inner, Measure, TestFunction, TrialFunction, sin, pi, SpatialCoordinate
from dolfinx import fem, io
import dolfinx.fem.petsc
from dolfinx.mesh import create_rectangle, CellType
from dolfinx.io import XDMFFile
from dolfinx import fem
import ufl
import pyvista
from dolfinx import plot
import matplotlib.pyplot as plt
from ufl import sym, grad, tr, Identity
from dolfinx.fem import Expression
import os
# ------------------------------
# Domain and mesh parameters (1×1 physical domain)
# ------------------------------
def epsilon(v):
return sym(grad(v))
def sigma(v,lmbda,dim,mu):
return lmbda * tr(epsilon(v)) * Identity(dim) + 2 * mu * epsilon(v)
def left(x):
return np.isclose(x[0], 0.0)
def right(x,length):
return np.isclose(x[0], length)
def bottom(x):
return np.isclose(x[1], 0.0)
def top(x,height):
return np.isclose(x[1], height)
def displ(u_sol,domain,dim,direction,Nx = 64):
if direction == 'x' :
i = 0
elif direction == 'y' :
i = 1
else:
return('wrong direction')
V = fem.functionspace(domain, ('CG',1,(dim,)))
u_interpolated = fem.Function(V,name='Interpolated Solution')
u_interpolated.interpolate(u_sol)
array = u_interpolated.x.array[i::2] # Displacements along x or y
dof_coord = V.tabulate_dof_coordinates()[:, :2]
N=Nx
indices = dof_coord / (1./N)
input_indices = np.round(indices).astype(np.int32)
y_sort = np.argsort(input_indices[:, 1])
sorted_coords = dof_coord[y_sort]
sorted_values = array[y_sort]
# Check that x-axis is sorted by increasing x coordinate (per y coord)
for i in range(N+1):
sub_indices = np.flatnonzero(input_indices[y_sort, 1]==i)
sub_sort = np.argsort(input_indices[y_sort[sub_indices], 0])
sub_coords = sorted_coords[sub_indices]
sub_vals = sorted_values[sub_indices]
sorted_coords[sub_indices] = sub_coords[sub_sort]
sorted_values[sub_indices] = sub_vals[sub_sort]
sorted_values = sorted_values.reshape(N, N)
return sorted_values, sorted_coords
def interp_scalar_with_coords(f_sol, domain, Nx=64):
"""
Take a scalar CG1 Function f_sol on 'domain' and return
- sorted_values: an (Nx x Nx) array of nodal values, ordered by y then x
- sorted_coords: the corresponding (Nx*Nx x 2) array of node coords
matching your displ() output.
"""
# 1) Build the same CG1 scalar space
V0 = fem.functionspace(domain, ("CG", 1))
f_i = fem.Function(V0, name=f_sol.name)
f_i.interpolate(f_sol)
# 2) Raw values + coordinates
vals = f_i.x.array # length = (Nx+1)*(Nx+1)
dof_coord = V0.tabulate_dof_coordinates()[:, :2]
# 3) Map coords → integer grid indices
indices = np.round(dof_coord / (1.0 / Nx)).astype(int) # shape=(n_dofs,2)
# 4) First sort by y, then by x within each row
# (exactly like your displ logic)
y_sort = np.argsort(indices[:, 1])
sorted_vals = vals[y_sort].copy()
sorted_coords= dof_coord[y_sort].copy()
sorted_idx = indices[y_sort]
for yy in range(Nx+1):
mask = np.where(sorted_idx[:,1] == yy)[0]
x_order = np.argsort(sorted_idx[mask, 0])
sorted_vals[mask] = sorted_vals[mask][x_order]
sorted_coords[mask] = sorted_coords[mask][x_order]
# 5) reshape into a 2D grid
grid_vals = sorted_vals.reshape(Nx, Nx) # note Nx+1 nodes per side
# if you really want 64×64, use reshape(Nx, Nx) on a 64×64 mesh of 63×63 elements.
return grid_vals, sorted_coords
def compute_fem(file) :
samples = np.load(file)[:2].shape[0]
print(samples)
solutions_disp = np.zeros((samples,2,64,64)) # shape 256,2,65,65
solutions_stress = np.zeros((samples, 3, 64, 64))
solutions_strain = np.zeros((samples, 3, 64, 64))
for i in range(samples):
length, height = 1.0, 1.0 # Physical domain: 1 x 1
Nx, Ny = 63, 63 # Mesh resolution
domain = create_rectangle(
MPI.COMM_WORLD,
[np.array([0.0, 0.0]), np.array([length, height])],
[Nx, Ny],
cell_type=CellType.quadrilateral,
)
dim = domain.topology.dim
degree = 2
#V = fem.functionspace(domain, ("P", degree, (dim,)))
V = fem.functionspace(domain, ('CG',1,(dim,)))
u_sol = fem.Function(V, name="Displacement")
# ------------------------------
# Load the GRF for Young's modulus E(x)
# ------------------------------
#E_vals = np.load(f"{file}")[:,:,i] # shape: (64, 64,4) JE CHANGE POUR ADAPTER AU NOUVEAU FORMAT (256,64,64)
E_vals = np.load(f"{file}")[i,:,:] # shape: (64, 64,4)
nx_grf, ny_grf = E_vals.shape
# Create grid coordinates corresponding to the GRF domain.
x_grf = np.linspace(0, length, nx_grf)
y_grf = np.linspace(0, height, ny_grf)
dx_grf = x_grf[1] - x_grf[0]
dy_grf = y_grf[1] - y_grf[0]
# ------------------------------
# Project GRF values onto the mesh
# ------------------------------
V_E = fem.functionspace(domain, ("CG", 1))
E_func = fem.Function(V_E, name="Youngs_modulus")
tdim = domain.topology.dim
domain.topology.create_connectivity(0, tdim) # vertex → cell
domain.topology.create_connectivity(tdim, 0) # cell → vertex
num_vertices = domain.topology.index_map(0).size_local
vertex_indices = np.arange(num_vertices, dtype=np.int32)
vertex_to_dof = fem.locate_dofs_topological(V_E, 0, vertex_indices)
flat = E_vals.flatten(order="C")[vertex_to_dof]
E_func.x.array[:] = flat
E_func.x.scatter_forward()
# ------------------------------
# Material parameters (spatially varying)
# ------------------------------
nu = fem.Constant(domain, 0.2) # Poisson's ratio
lmbda = (E_func * nu) / ((1 + nu) * (1 - (2 * nu)))
mu = (E_func) / (2 * (1 + nu))
# ------------------------------
# Variational problem setup
# ------------------------------
u = TrialFunction(V)
v = TestFunction(V)
rho = 2e-3
g = 9.81
f = fem.Constant(domain, np.array([0, -rho * g]))
dx_measure = Measure("dx", domain=domain)
a = inner(sigma(u,lmbda,dim,mu), epsilon(v)) * dx_measure
L = inner(f, v) * dx_measure
# ------------------------------
# Boundary conditions (Fixed on Bottom, Left, and Right)
# ------------------------------
left_dofs = fem.locate_dofs_geometrical(V, left)
right_dofs = fem.locate_dofs_geometrical(V, lambda x: right(x, length))
bottom_dofs = fem.locate_dofs_geometrical(V, bottom)
bcs = [
fem.dirichletbc(np.zeros((dim,), dtype=np.float64), left_dofs, V),
fem.dirichletbc(np.zeros((dim,), dtype=np.float64), right_dofs, V),
fem.dirichletbc(np.zeros((dim,), dtype=np.float64), bottom_dofs, V),
]
top_dofs = fem.locate_dofs_geometrical(V,lambda x: top(x, height))
# Define spatial coordinate
x = SpatialCoordinate(domain)
# Define cosine load function in UFL form
F0 = 10 # Load magnitude
cosine_load = F0 * sin(pi * x[0]) # f_y(x) = F0 * sin(pi * x), zero at x=0, x=1, max at x=0.5
# Apply as Neumann boundary condition (traction force) to the y-component
ds = Measure("ds", domain=domain)
L += inner(cosine_load, v[1]) * ds # Fixed: Apply force only to the y-component
# The fix ensures shape compatibility in UFL expressions
# ------------------------------
# Solve the linear elasticity problem
# ------------------------------
problem = fem.petsc.LinearProblem(
a, L, u=u_sol, bcs=bcs,
petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
)
problem.solve()
# ------------------------------
# Output the solution
# ------------------------------
vtk = io.VTKFile(domain.comm, "displacement_field.pvd", "w")
vtk.write_function(u_sol)
vtk.close()
ux_array, sorted_coords_x = displ(u_sol,domain,dim,'x',Nx = 64)
uy_array, sorted_coords_y = displ(u_sol,domain,dim,'y',Nx = 64)
ux_array = ux_array[::-1, :]
uy_array = uy_array[::-1, :]
solutions_disp[i,0,:,:] = ux_array
solutions_disp[i,1,:,:] = uy_array
dx = Measure("dx", domain=domain)
eps = sym(grad(u_sol))
sigma_expr = lmbda*tr(eps)*Identity(dim) + 2*mu*eps
W = fem.functionspace(domain, ("CG", 1, (dim, dim)))
V1 = fem.functionspace(domain, ("CG", 1))
expr_xx = Expression(sigma_expr[0,0], V1.element.interpolation_points())
expr_yy = Expression(sigma_expr[1,1], V1.element.interpolation_points())
expr_xy = Expression(sigma_expr[0,1], V1.element.interpolation_points())
sig_xx = fem.Function(V1, name="sigma_xx")
sig_yy = fem.Function(V1, name="sigma_yy")
sig_xy = fem.Function(V1, name="sigma_xy")
sig_xx.interpolate(expr_xx)
sig_yy.interpolate(expr_yy)
sig_xy.interpolate(expr_xy)
# write them all into one PVD
vtk_s = io.VTKFile(domain.comm, "stress_components.pvd", "w")
vtk_s.write_function(sig_xx)
vtk_s.write_function(sig_yy)
vtk_s.write_function(sig_xy)
vtk_s.close()
# CG1 projection space
p = TrialFunction(V1)
q = TestFunction(V1)
a_proj = inner(p, q) * dx
# components: (xx,yy,xy)
for comp, (i0,i1) in enumerate([(0,0), (1,1), (0,1)]):
expr = sigma_expr[i0, i1]
L_proj = inner(expr, q) * dx
proj = fem.petsc.LinearProblem(
a_proj, L_proj, bcs=[],
petsc_options={"ksp_type":"preonly","pc_type":"lu"}
)
sigma_cg1 = proj.solve() # CG1 nodal Function
# reorder to 64×64 array
grid_vals, _ = interp_scalar_with_coords(sigma_cg1, domain, Nx=64)
solutions_stress[i, comp, :, :] = grid_vals
for comp, (i0,i1) in enumerate([(0,0), (1,1), (0,1)]):
expr_ep = eps[i0, i1]
L_proj = inner(expr_ep, q) * dx
proj = fem.petsc.LinearProblem(
a_proj, L_proj, bcs=[],
petsc_options={"ksp_type":"preonly","pc_type":"lu"}
)
eps_cg1 = proj.solve() # CG1 nodal Function
# reorder to 64×64 array
grid_vals, _ = interp_scalar_with_coords(eps_cg1, domain, Nx=64)
solutions_strain[i, comp, :, :] = grid_vals
if i >= samples - 1 :
N = 64
X_ = sorted_coords_y[:, 0].reshape(N, N)
Y_ = sorted_coords_y[:, 1].reshape(N, N)
fig = plt.figure()
ax = plt.gca()
im1 = ax.imshow(uy_array, cmap="viridis")
fig.colorbar(im1)
return solutions_disp, solutions_stress[:,:,::-1,:] , solutions_strain[:,:,::-1,:]
Thank you in advance for your help!