Running your code with the additional line:
print(dolfinx.__version__, dolfinx.git_commit_hash) yields:
FEA elapsed time: 0.013041 seconds
FEA elapsed time: 0.009913 seconds
FEA elapsed time: 0.022576 seconds
FEA elapsed time: 0.033547 seconds
FEA elapsed time: 0.082800 seconds
FEA elapsed time: 0.129380 seconds
FEA elapsed time: 0.373714 seconds
FEA elapsed time: 0.700876 seconds
FEA elapsed time: 2.039325 seconds
FEA elapsed time: 4.543628 seconds
FEA elapsed time: 12.998073 seconds
0.11.0.post0 fefdb2201b80a8f59527de2d461b9056906661d8
[[0.002 ]
[0.00168179]
[0.00141421]
[0.00118921]
[0.001 ]
[0.0008409 ]
[0.00070711]
[0.0005946 ]
[0.0005 ]
[0.00042045]
[0.00035355]]
[[0.00145655]
[0.00147336]
[0.00089138]
[0.00265791]
[0.00206985]
[0.00304512]
[0.00301411]
[0.00394732]
[0.00352334]
[0.00390616]
[0.00401537]]
[[1841.58226109]
[2034.78366816]
[1075.04116855]
[2240.10755091]
[1554.49175714]
[2822.55699444]
[2325.71849015]
[3492.19236626]
[2919.02979709]
[3469.49339551]
[3544.64612005]]
code:
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
r_loading = 2.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=10.0e4, 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 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, node_normal,
p_traction=1.0/(np.pi*r_loading**2), 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))
# Traction direction: push along +x (change sign if needed)
T = fem.Constant(domain, -p_traction * default_scalar_type((node_normal[0], 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)
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):
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_normal)
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))[0]
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, dof_coords[dof_index], 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 = 11
cell_sizes = np.empty((num_loops,1))
u_disps = np.empty((num_loops,1))
stresses = np.empty((num_loops,1))
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, point, u_node = fea_solve(domain, grasp_facet_node, grasp_node2, grasp_node_normals)
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, f"max_p_stress_part_iter_{n+1}.xdmf", "w") as xdmf:
xdmf.write_mesh(domain_f)
uh.name = "u"
xdmf.write_function(uh)
pmax.name = "PrincipalStress"
xdmf.write_function(pmax)
print(dolfinx.__version__, dolfinx.git_commit_hash)
print('\n')
print(cell_sizes)
print('\n')
print(u_disps)
print('\n')
print(stresses)
Running on 0.10.0 yields the results that you are seeing:
FEA elapsed time: 0.706027 seconds
FEA elapsed time: 0.010824 seconds
FEA elapsed time: 0.023024 seconds
FEA elapsed time: 0.034646 seconds
FEA elapsed time: 0.084435 seconds
FEA elapsed time: 0.129224 seconds
FEA elapsed time: 0.367234 seconds
FEA elapsed time: 0.691599 seconds
FEA elapsed time: 2.003611 seconds
FEA elapsed time: 4.567333 seconds
FEA elapsed time: 13.119260 seconds
0.10.0 b8375fa9fb932cbfce43005eeb4a3d52293f428b
[[0.002 ]
[0.00168179]
[0.00141421]
[0.00118921]
[0.001 ]
[0.0008409 ]
[0.00070711]
[0.0005946 ]
[0.0005 ]
[0.00042045]
[0.00035355]]
[[0.00145655]
[0.00147336]
[0.00849693]
[0.01369466]
[0.01771337]
[0.0261431 ]
[0.03383505]
[0.0467272 ]
[0.06021858]
[0.07440669]
[0.0972786 ]]
[[ 1841.58226109]
[ 2034.78366816]
[11929.89750367]
[12856.75539342]
[15471.79283942]
[26687.19345534]
[28523.12707249]
[42557.26616307]
[54022.69178219]
[69590.04903437]
[92535.84861181]]
while 0.9 is back to what 0.10 had
FEA elapsed time: 0.170569 seconds
FEA elapsed time: 0.010710 seconds
FEA elapsed time: 0.022907 seconds
FEA elapsed time: 0.033624 seconds
FEA elapsed time: 0.083407 seconds
FEA elapsed time: 0.129651 seconds
FEA elapsed time: 0.367752 seconds
FEA elapsed time: 0.694417 seconds
FEA elapsed time: 2.005706 seconds
FEA elapsed time: 4.524475 seconds
FEA elapsed time: 13.015335 seconds
0.9.0 4116cca8cf91f1a7e877d38039519349b3e08620
[[0.002 ]
[0.00168179]
[0.00141421]
[0.00118921]
[0.001 ]
[0.0008409 ]
[0.00070711]
[0.0005946 ]
[0.0005 ]
[0.00042045]
[0.00035355]]
[[0.00145655]
[0.00147336]
[0.00089138]
[0.00265791]
[0.00206985]
[0.00304512]
[0.00301411]
[0.00394732]
[0.00352334]
[0.00390616]
[0.00401537]]
[[1841.58226109]
[2034.78366816]
[1075.04116855]
[2240.10755091]
[1554.49175714]
[2822.55699444]
[2325.71849015]
[3492.19236626]
[2919.02979709]
[3469.49339551]
[3544.64612005]]
which to me indicates that what you are seeing is this bug:
as the facets marked with 2 in your code includes interior nodes, i.e. the following would trigger:
domain.topology.create_connectivity(domain.topology.dim-1, domain.topology.dim)
f_to_c = domain.topology.connectivity(domain.topology.dim-1, domain.topology.dim)
for facet in facet_tags.find(2):
cells = f_to_c.links(facet)
assert len(cells) == 1, "Facet should be on the boundary"
this has been corrected in the 0.11 release.