Hi,
I wrote the following simple code to solve for linear elasticity and to compute the maximum principal stress. The displacement elements are linear and the computed stresses are interpolated to a DG 0 space.
-
My first question is: from a quick look of the code, is there anything I can do to make it simpler/easier to read?
-
I am using a wrapper for the CGAL library for meshing, and there are always 1 or 2 low quality nodes generated. How could I compute the mesh quality and find out what elements not to consider when I am looking for the maximum principal stress:
max_stress_per_N = np.max(pmax.x.array)
-
I tried increasing the order of the Lagrange elements by 1, and using DG 1 for the stress, but the program crashes. Is there a reason for this?
-
When I look for the displacement I use the line below, is there another way of obtaining this if you know the node, but for an arbitrary element degree?
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)
This is the full code:
import numpy as np
#import pygalmesh
import meshio
import time
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import basix.ufl as bu
from dolfinx import mesh, fem, io, default_scalar_type
from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, set_bc
r_loading = 3.0e-3
def make_dolfinx_mesh_from_pygalmesh(vol_mesh):
points = vol_mesh.points
tet_cells = None
for c in vol_mesh.cells:
if c.type in ["tetra", "tetra10", "tetrahedron"]:
tet_cells = c.data
break
if tet_cells is None:
raise RuntimeError("No tetrahedral cells found in the volume mesh.")
# Vector Lagrange P1 element for mesh
finiteElement = bu.element("Lagrange", "tetrahedron", 1, shape=(3,))
domain = mesh.create_mesh(MPI.COMM_WORLD, cells=tet_cells, x=points, e=finiteElement)
return domain
def mark_boundaries(domain, facet_node,grasp_node,grasp_node_normal, eps=1e-3):
dim = domain.topology.dim
fdim = dim - 1
domain.topology.create_connectivity(fdim, dim)
def left_facets(xq):
return np.sqrt((xq[0]-facet_node[0])**2 + (xq[1]-facet_node[1])**2 + (xq[2]-facet_node[2])**2) < r_loading
def right_facets(xq):
return np.sqrt((xq[0]-grasp_node[0])**2 + (xq[1]-grasp_node[1])**2 + (xq[2]-grasp_node[2])**2) < 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.35):
"""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 setup_function_space(domain):
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim,)))
return V
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, -1000.0*9.8)):
dim = domain.topology.dim
fdim = dim - 1
ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags)
V = setup_function_space(domain)
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 (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 solve_linear_system(domain, V, a, L, bcs):
a_form = fem.form(a)
L_form = fem.form(L)
A = assemble_matrix(a_form, bcs=bcs)
A.assemble()
b = assemble_vector(L_form)
apply_lifting(b, [a_form], bcs=[bcs])
b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, bcs)
uh = fem.Function(V)
solver = PETSc.KSP().create(domain.comm)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
solver.setOperators(A)
solver.solve(b, uh.x.petsc_vec)
uh.x.scatter_forward()
return uh
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-12) # p = 2/sqrt(3)*sqrt(J2) (small epsilon added)
r = 4 * b / (p**3) # r = 4 J3 / p^3
r = ufl.max_value(ufl.min_value(r, +1-1e-6), -1+1e-6) # 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
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)
uh = solve_linear_system(domain, Vp, a_p, L_p, bcs_p)
# Principal stress
pmax = principal_stress(domain, uh, lambda_, mu)
# 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, u_node, pmax
# Part dependent mesh for FEA
##################################################################
# min_facet_angle = 25.0
# max_radius_surface_delaunay_ball = 0.0006
# max_facet_distance = 0.0006
# max_circumradius_edge_ratio = 3.0
# t_start = time.perf_counter()
# vol_mesh = pygalmesh.generate_volume_mesh_from_surface_mesh("bracket_ver2.STL",
# min_facet_angle=min_facet_angle,
# max_radius_surface_delaunay_ball=max_radius_surface_delaunay_ball,
# max_facet_distance=max_facet_distance,
# max_circumradius_edge_ratio=max_circumradius_edge_ratio,
# odt=True, lloyd=True,
# verbose=False,
# )
# t_end = time.perf_counter()
# print(f"Meshing elapsed time: {t_end - t_start:.6f} seconds")
# vol_mesh.write("mesh_braket.vtk")
domain = make_dolfinx_mesh_from_pygalmesh(meshio.read("mesh_braket.vtk"))
grasp_facet_node = [0.02, 0.018, 0.008]
grasp_node2 = [0.0, 0.018, 0.008]
grasp_node_normals = [-1.0, 0.0, 0.0]
t_start = time.perf_counter()
domain_f, uh, u_at_node, pmax = 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")
print('\n')
with io.XDMFFile(domain_f.comm, "max_p_stress_braket.xdmf", "w") as xdmf:
xdmf.write_mesh(domain_f)
uh.name = "u"
xdmf.write_function(uh)
pmax.name = "PrincipalStress"
xdmf.write_function(pmax)
max_stress = np.max(pmax.x.array)
print(u_at_node)
print(max_stress)
The mesh .vtk file can be found here:
Thank you,
Alex