[EDITED BASED ON COMMENTS] I having some trouble updating some code from legacy Fenics to DOLFINx. In the code a temperature field is defined using a Karhunen-Loeve transformation via UserExpression. From reading other posts I understand that UserExpression is no longer supported in DOLFINx but I still haven’t found a way to reformulate the expression.
Here is an example of the problem using a mean temperature field.
from mpi4py import MPI
import numpy as np
import dolfinx
import ufl
import matplotlib.pyplot as plt
from petsc4py.PETSc import ScalarType
#WORKING FEA CODE
x_0 = 0
x_l = 10
Num_elements = 100
Bar_mesh = dolfinx.mesh.create_interval(comm=MPI.COMM_WORLD, nx=Num_elements, points=. [x_0, x_l])
V = dolfinx.fem.FunctionSpace(Bar_mesh, ("Lagrange", 1))
Bar_boundary_elements = dolfinx.mesh.locate_entities_boundary(Bar_mesh, dim=0, marker=lambda x:
np.logical_or(np.isclose(x[0], x_0), np.isclose(x[0], x_l)))
boundary_dofs = dolfinx.fem.locate_dofs_topological(V=V, entity_dim=0, entities=Bar_boundary_elements)
bc = dolfinx.fem.dirichletbc(value=ScalarType(0), dofs=boundary_dofs, V=V)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(Bar_mesh)
a = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
T = .5 * 10 * (x[0] - 5) ** 2 + 200
L = 0.0032 * v.dx(0) * T * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Solution = problem.solve()
x_loc = np.linspace(x_0, x_l, Num_elements+1)
t = .5 * 10 * (x_loc - 5) ** 2 + 200
fig1 = plt.figure()
ax1 = fig1.add_subplot(1, 1, 1)
ax1.plot(x_loc, Solution.x.array.real, label='FE Solution')
ax1.legend(loc='lower left')
ax1.set_xlabel("Bar Location", fontsize=14)
ax1.set_ylabel("Displacement", fontsize=14)
ax11 = ax1.twinx()
ax11.plot(x_loc, t, color="red", label='Temperature Field')
ax11.legend(loc='lower right')
ax11.set_ylabel("Temperature", color="red", fontsize=14)
plt.show()
However for my problem I want to incorporate a random temperature field using a Karhunen-Loeve transformation, in place of T = .5 * 10 * (x[0] - 5) ** 2 + 200. In the legacy version of Fenics this was done as follows using a UserExpression.
fkls = np.zeros(self.nmesh)
for j in range(0, self.ndim):
fkls += np.sqrt(self.wlam[j]) * self.vhat[:, j] * self.yhat[isim, j]
class MyExpression0(UserExpression):
def __init__(self, b, nmesh, **kwargs):
super().__init__(**kwargs)
self.A = b
self.nmesh = nmesh
def eval(self, value, x):
lh = (xL - x0) / 2.
i1 = np.min([np.floor(x[0] / deltax).astype(int), self.nmesh - 1])
value[0] = .5 * 10 * (x[0] - lh) ** 2 + 200 + self.A * fkls[i1]
def value_shape(self):
return ()
T = MyExpression0(Acoeff, self.nmesh, degree=2)
However since UserExpression is depreciated in DOLFINx I’m unsure how to make the change. Any help would be appreciated!
ACCEPTED SOLUTION - with help from nate
#CREAT RANDOM TEMPURATURE FIELD
def inputUncertaintyLHS(coords, nmesh, ndim, Npt):
x0 = coords[0]
xL = coords[1]
x = np.linspace(x0, xL, nmesh)
k0 = .25
R = np.exp(-k0 * x ** 2)
fcov = np.zeros([nmesh, nmesh])
flim = np.array([-3, 3]) # upper and lower bound for scaling Chebychev roots
# to underlying input uncertainty
for j in range(0, nmesh):
for j1 in range(j, nmesh):
ixi = j1 - j
fcov[j, j1] = R[ixi]
fcov = fcov + np.transpose(fcov) - np.eye(nmesh) * np.diag(fcov)
(wlam, vhat) = np.linalg.eig(fcov)
wlam = np.real(wlam) # these should be real b/c fcov is SPD matrix but due to numerical
vhat = np.real(vhat) # errors there are very small imaginary values that i'm eliminating here
mfac = np.zeros(nmesh)
wsum = np.sum(wlam)
for j in range(0, nmesh):
mfac[j] = np.sum(wlam[0:j]) / wsum
pts = np.random.rand(Npt, ndim)
yhat = (flim[1] - flim[0]) * pts + flim[0]
mux = 0.0
sigx = 1.0
yhat = ((yhat - mux) / sigx).astype(np.float32)
vhat = vhat[:, 0:ndim]
return (yhat, wlam, vhat)
Acoeff = 32
Num_dim = 4
Num_pt = 1
(yhat, wlam, vhat) = inputUncertaintyLHS([x_0, x_l], Num_elements, Num_dim, Num_pt)
fkls = np.zeros(Num_elements)
for n in range(0, 4):
fkls += np.sqrt(wlam[n]) * vhat[:, n] * yhat[0, n]
def my_expression(x):
lh = (x_l - x_0) / 2.
deltax = (x_l - x_0) / Num_elements
i1 = np.minimum(np.floor(x[0] / deltax).astype(int), Num_elements - 1).astype(int)
return .5 * 10 * (x[0] - lh) ** 2 + 200 + Acoeff * fkls[i1]
T_new= dolfinx.fem.Function(V)
T_new.interpolate(my_expression)
L = 0.0032 * v.dx(0) * T_new * ufl.dx
problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
Solution = problem.solve()
x_loc = np.linspace(x_0, x_l, Num_elements+1)
fig1 = plt.figure()
ax1 = fig1.add_subplot(1, 1, 1)
ax1.plot(x_loc, Solution.x.array.real, label='FE Solution')
ax1.legend(loc='lower left')
ax1.set_xlabel("Bar Location", fontsize=14)
ax1.set_ylabel("Displacement", fontsize=14)
plt.show()
