Hi,
I would like to load a rectangle with a traction boundary condition on a face normal to the x axis. The traction is applied on a circle in that face, and the traction follows a radially symmetric function, that is pre-computed for 10 points and then accessed using np.interp().
My current self-contained code is attached below. When I run the code I get the following warning:
/usr/lib/python3/dist-packages/ufl/core/terminal.py:60: UserWarning: Couldn't map 'f' to a float, returning ufl object without evaluation.
And the code freezes. Am I declaring the correct space in the function get_traction_yz()? If I’m applying the traction to a face, should the space be 2D?
If tried the following instead:
Trac_vals = fem.functionspace(domain, (“Lagrange”, 1, (domain.geometry.dim-1,)))
But get the following error:
File /usr/lib/petsc/lib/python3/dist-packages/dolfinx/fem/function.py:484 in interpolate
self._cpp_object.interpolate(np.asarray(u0(x), dtype=self.dtype), cells0) # type: ignoreRuntimeError: Interpolation data has the wrong shape/size.
Another option is not to use np.interp(), like this:
def traction_fun(x, g_n):
r = ufl.sqrt((x[1]-g_n[1])**2 + (x[2]-g_n[2])**2)
return ufl.conditional(ufl.gt(r, 2.0e-3), 2e4*(1.0 - (r-2.0e-3)**2/(4.0e-3-2.0e-3)**2), 2e4)
def get_traction_yz(domain, g_n, pre_rs, pre_vals):
Trac_vals = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim-1,)))
trac_vals = fem.Function(Trac_vals)
trac_vals.interpolate(lambda x: traction_fun(x, g_n))
return trac_vals
But I get this error:
File /usr/lib/python3/dist-packages/ufl/constantvalue.py:511 in as_ufl
raise ValueError(
ValueError: Invalid type conversion: [5.0e-05 5.0e-05 3.4e-05 ... 5.0e-05 3.4e-05 5.0e-05] can not be converted to any UFL type.
I would like to be able to use a function with ufl expressions and if possible something like what is coded now, using np.interp() for a pointwise defined function
Thank you,
Alex
import numpy as np
import time
from mpi4py import MPI
import ufl
from dolfinx import mesh, fem, io, default_scalar_type
import dolfinx.fem.petsc
import matplotlib.pyplot as plt
r_loading = 4.0e-3
def in_sphere(x, point, r_range):
return ((x[0]-point[0])**2 + (x[1]-point[1])**2 + (x[2]-point[2])**2) < r_range**2
def mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal):
dim = domain.topology.dim
fdim = dim - 1
domain.topology.create_connectivity(fdim, dim)
def left_facets(xq):
return in_sphere(xq, facet_node, r_loading)
def right_facets(xq):
return in_sphere(xq, grasp_node, r_loading)
facets_left = mesh.locate_entities(domain, fdim, left_facets)
facets_right = mesh.locate_entities(domain, fdim, right_facets)
facet_indices = np.concatenate([facets_left, facets_right])
facet_markers = np.concatenate([np.full_like(facets_left, 1, dtype=np.int32),
np.full_like(facets_right, 2, dtype=np.int32)])
# Sort by facet index
sort_perm = np.argsort(facet_indices)
facet_indices = facet_indices[sort_perm]
facet_markers = facet_markers[sort_perm]
mt = mesh.meshtags(domain, fdim, facet_indices, facet_markers)
return mt
def elasticity_parameters(E=1.0e6, nu=0.20):
"""Lamé parameters from E, nu."""
lambda_ = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))
mu = E * 0.5 / (1.0 + nu)
return lambda_, mu
def get_traction_yz(domain, g_n, pre_rs, pre_vals):
Trac_vals = fem.functionspace(domain, ("DG", 0))
trac_vals = fem.Function(Trac_vals)
trac_vals.interpolate(lambda x: np.interp(np.sqrt((x[1]-g_n[1])**2 + (x[2]-g_n[2])**2), pre_rs,pre_vals))
return trac_vals
def epsilon(u):
return ufl.sym(ufl.grad(u))
def sigma(u, lambda_, mu):
return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2.0 * mu * epsilon(u)
def build_primal_problem(domain, facet_tags, lambda_, mu, grasp_node, node_normal, lock_num, pre_rs, pre_vals, body_force=(0.0, 0.0, 0.0)):
dim = domain.topology.dim
fdim = dim - 1
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f_vec = fem.Constant(domain, default_scalar_type(body_force))
if lock_num == 1:
p_traction = get_traction_yz(domain, grasp_node, pre_rs, pre_vals)
# Traction direction: push along +x (change sign if needed)
T = fem.Constant(domain, -default_scalar_type((node_normal[0]*p_traction, node_normal[1], node_normal[2])))
a = ufl.inner(sigma(u, lambda_, mu), epsilon(v)) * ufl.dx
L = ufl.dot(f_vec, v) * ufl.dx + ufl.dot(T, v) * ds(2)
# Dirichlet BC on left jaw (tag 1)
left_facets = facet_tags.find(1)
domain.topology.create_connectivity(fdim, dim)
left_dofs = fem.locate_dofs_topological(V, fdim, left_facets)
u_D = np.array((0.0, 0.0, 0.0), dtype=default_scalar_type)
bc_left = fem.dirichletbc(u_D, left_dofs, V)
# Dirichlet BC to stop deflection normal to x (lock_num == 1) or y (lock_num == 2) axis
if lock_num == 1:
right_facets = facet_tags.find(2)
right_dofs1 = fem.locate_dofs_topological(V.sub(1), fdim, right_facets)
bc_1 = fem.dirichletbc(0.0, right_dofs1, V.sub(1))
right_dofs2 = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
return V, a, L, [bc_left, bc_1, bc_2]
if lock_num == 2:
right_facets = facet_tags.find(2)
right_dofs0 = fem.locate_dofs_topological(V.sub(0), fdim, right_facets)
bc_1 = fem.dirichletbc(0.0, right_dofs0, V.sub(0))
right_dofs2 = fem.locate_dofs_topological(V.sub(2), fdim, right_facets)
bc_2 = fem.dirichletbc(0.0, right_dofs2, V.sub(2))
return V, a, L, [bc_left, bc_1, bc_2]
return V, a, L, [bc_left]
def principal_stress(domain, uh, lambda_, mu):
#Compute principal stresses via invariants (analytical eigenvalue solution)
sigma_uh = sigma(uh, lambda_, mu)
q = ufl.tr(sigma_uh) / 3 # mean stress (I1/3)
B = sigma_uh - q * ufl.Identity(3) # deviatoric stress tensor
j = ufl.tr(B * B) / 2 # J2 = 1/2 tr(B^2)
b = ufl.tr(B * B * B) / 3 # J3 = 1/3 tr(B^3)
p = 2 / ufl.sqrt(3) * ufl.sqrt(j + 1e-14) # p = 2/sqrt(3)*sqrt(J2) (small eps**2 added)
r = 4 * b / (p**3) # r = 4 J3 / p^3
r = ufl.max_value(ufl.min_value(r,+1-1e-7),-1+1e-7) # clamp r to [-1+eps, 1-eps] to avoid acos domain error
phi = ufl.acos(r) / 3
lambda2 = q + p * ufl.cos(phi) # largest principal stress (max)
V_dg = fem.functionspace(domain, ("DG", 0))
pmax_expr = fem.Expression(lambda2, V_dg.element.interpolation_points)
pmax = fem.Function(V_dg, name="PrincipalStress")
pmax.interpolate(pmax_expr)
return pmax, V_dg
def array_dist_squared(x, p):
return (x[:,0]-p[0])**2 + (x[:,1]-p[1])**2 + (x[:,2]-p[2])**2
def fea_solve(domain, facet_node, grasp_node, grasp_node_normal, lock_num, pre_rs, pre_vals):
lambda_, mu = elasticity_parameters()
facet_tags = mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal)
Vp, a_p, L_p, bcs_p = build_primal_problem(domain, facet_tags, lambda_, mu, grasp_node, grasp_node_normal, lock_num, pre_rs, pre_vals)
problem = dolfinx.fem.petsc.LinearProblem(a_p, L_p, bcs=bcs_p,
petsc_options_prefix="linear_problem",
petsc_options={"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"})
uh = problem.solve()
# Principal stress field
pmax, V_dg = principal_stress(domain, uh, lambda_, mu)
# MPS maximum value not considering low E area
dof_coords = V_dg.tabulate_dof_coordinates()
indices = np.where(np.logical_and(array_dist_squared(dof_coords, facet_node) > 5.0e-3**2, \
array_dist_squared(dof_coords, grasp_node) > 5.0e-3**2))
max_pmax = np.sort(pmax.x.array[indices])[-1]
# Get displacement at node where load is applied
dof_coords = Vp.tabulate_dof_coordinates()
distances = np.linalg.norm(dof_coords - grasp_node, axis=1)
dof_index = np.argmin(distances)
u_node = np.sqrt(uh.x.array[3*dof_index]**2 + uh.x.array[3*dof_index+1]**2 + uh.x.array[3*dof_index+2]**2)
return domain, uh, pmax, max_pmax, u_node
# Part FEA
##################################################################
length = 20.0e-3
width = 10.0e-3
grasp_facet_node = np.array([0.0, width/2, width/2])
grasp_node2 = np.array([length, width/2, width/2])
grasp_node_normals = np.array([1.0, 0.0, 0.0 ])
cell_size = 2.0e-3
num_loops = 1
cell_sizes = np.empty((num_loops,1))
u_disps = np.empty((num_loops,1))
stresses = np.empty((num_loops,1))
lock_num = 1
def gen_pressure_profile(ri, ro, pmag, num_ps):
pre_r = np.linspace(0.0, ro, num=50, endpoint=True)
pre_val = np.zeros(pre_r.size)
for nn in range(pre_r.size):
if pre_r[nn] < ri:
pre_val[nn] = pmag
else:
pre_val[nn] = pmag*(1.0 - (pre_r[nn]-ri)**2/(ro-ri)**2)
return pre_r, pre_val
r_i = 2.0e-3
r_o = r_loading
p_mag = 1.0/(np.pi*r_loading**2)
pre_rs, pre_vals = gen_pressure_profile(r_i, r_o, p_mag, 10)
plt.figure()
plt.plot(pre_rs, pre_vals)
plt.show()
for n in range(num_loops):
n_l = int(length/cell_size)
n_w = int(width/cell_size )
domain = mesh.create_box(MPI.COMM_WORLD, ((0.0, 0.0, 0.0), (length,width,width)), (n_l,n_w,n_w), mesh.CellType.tetrahedron)
t_start = time.perf_counter()
domain_f, uh, pmax, max_stress, u_node = fea_solve(domain, grasp_facet_node, grasp_node2, grasp_node_normals, lock_num, pre_rs, pre_vals)
t_end = time.perf_counter()
print(f"FEA elapsed time: {t_end - t_start:.6f} seconds")
cell_sizes[n] = cell_size
u_disps[n] = u_node
stresses[n] = max_stress
cell_size = cell_size/2.0**0.25
with io.XDMFFile(domain_f.comm, "max_p_stress_part.xdmf", "w") as xdmf:
xdmf.write_mesh(domain_f)
uh.name = "u"
xdmf.write_function(uh)
pmax.name = "PrincipalStress"
xdmf.write_function(pmax)
print('\n')
print(cell_sizes)
print('\n')
print(u_disps)
print('\n')
print(stresses)