Hi,
I was able to run two elastic simulations using the original domain, and then a subdivided domain using the following code:
#Make new domain with all elements subdivided once
tdim = domain.topology.dim
domain.topology.create_entities(1) # edges
domain.topology.create_connectivity(1, tdim) # edge → cell
# IMPORTANT: refine returns a tuple
domain2,_,_ = mesh.refine(domain)
tdim = domain2.topology.dim
domain2.topology.create_entities(1) # edges
domain2.topology.create_connectivity(1, tdim) # edge → cell
The code seems to be working, but I am not sure why create entities and connectivity needs to be called the second time. If refine needs it, but the solver doesn’t, why do I need to create entities and connectivity the second time?
Thank you,
Alex
A self-contained running code is shown below:
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_parameters2(domain, grasp_node, 4.0e-3)
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 = 0.6e-3
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)
#Make new domain with all elements subdivided once
tdim = domain.topology.dim
domain.topology.create_entities(1) # edges
domain.topology.create_connectivity(1, tdim) # edge → cell
# IMPORTANT: refine returns a tuple
domain2,_,_ = mesh.refine(domain)
tdim = domain2.topology.dim
domain2.topology.create_entities(1) # edges
domain2.topology.create_connectivity(1, tdim) # edge → cell
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")
with io.XDMFFile(domain_f.comm, f"max_p_stress_part_no_sub.xdmf", "w") as xdmf:
xdmf.write_mesh(domain_f)
uh.name = "u"
xdmf.write_function(uh)
pmax.name = "PrincipalStress"
xdmf.write_function(pmax)
t_start = time.perf_counter()
domain_f, uh, pmax, max_stress2, _, _ = fea_solve(domain2, grasp_facet_node, grasp_node2, grasp_node_normals)
t_end = time.perf_counter()
print(f"FEA elapsed time: {t_end - t_start:.6f} seconds")
with io.XDMFFile(domain_f.comm, f"max_p_stress_part_one_sub.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((max_stress2-max_stress)/max_stress2)
print('\n')