How to compute maximum stresses in elstic solver outside of sphere and not in nodes belonging to low quality elements

General
1

Hi,

I would like to compute the highest maximum principal stress outside of a sphere (specified by function in_sphere()) and not in nodes belonging to low quality elements. I have a function that can compute the element quality given the 4 tetrahedron nodes, belonging to a given node. How could I go about adding this as an output to the fea_solve() function? What I am not very clear on is how the DG0 element stress values are projected to the nodes; is this like an average of the elements around the node?

The short version of the code is posted below.

The mesh can be found here: Dropbox

Thank you,
Alex

import numpy as np
import pygalmesh
import meshio
import time


from mpi4py import MPI

import ufl
import basix.ufl as bu
from dolfinx import mesh, fem, io, default_scalar_type



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 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_parameters2(domain, point, r_range, E1=10.0e4, E2=1.0e4, nu=0.35):
    """Lamé parameters from E, nu."""
    
    param1 = nu  / ((1.0 + nu) * (1.0 - 2.0 * nu))
    param2 = 0.5 /  (1.0 + nu)
    
    lambda1_ = E1 * param1
    mu1      = E1 * param2
    
    lambda2_ = E2 * param1
    mu2      = E2 * param2
    
    Lambda_ = fem.functionspace(domain, ("DG", 0))
    lambda_ = fem.Function(Lambda_)
    lambda_.interpolate(lambda x: np.where(in_sphere(x, point, r_range), lambda2_, lambda1_))
    
    Mu = fem.functionspace(domain, ("DG", 0))
    mu = fem.Function(Mu)
    mu.interpolate(lambda x: np.where(in_sphere(x, point, r_range), mu2, mu1))
      
    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, -4000.0*9.8)):
    """
    BCs:
      - left side (tag 1): fixed
      - right side (tag 2): traction in +x direction
    """
     
    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-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_parameters2(domain, facet_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 = 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
    pmax = principal_stress(domain, uh, lambda_, mu)
    
    return domain, uh, 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,
#     seed=0, 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, 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_per_N = np.max(pmax.x.array)
#max_stress_per_N = np.sort(pmax.x.array)[-4]


print(max_stress_per_N)
2

Both element quality and principal stress are not defined for nodes, so I would recommend analyzing the solution from values in the element centers instead. This is well defined for both mesh quality and the principal stress and would involve looking at the DG0-dof of the computed principal stress function for a given element. Is there a reason you need the values in mesh nodes specifically?

How to interpret the DG0 stress values in the nodes is not well defined, since the value is computed from the solution gradient, which is not defined in nodes. If you really need a value of the gradient, or a value computed from the gradient, you could use a recovery operator like in a posteriori error estimation. This does not give you the actual value of the gradient of the computed solution or the actual solution, however. Clement interpolation is one such recovery operator and is very similar to your suggestion of taking the average from elements around the node.

3

Hi, thank you for the reply. Regarding not considering stresses inside a sphere centered at a point, I was able to do this:

# MPS 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(array_dist_squared(dof_coords, facet_node) > 4.5e-3**2)
max_pmax   = np.sort(pmax.x.array[indices])[-1]

Passing the DG 0 function space (V_dg) from the principal_stress() function, and using this function:

def array_dist_squared(x, p):
    
    return (x[:,0]-p[0])**2 + (x[:,1]-p[1])**2 + (x[:,2]-p[2])**2

Is this correct?

For not considering low quality elements, one way to do it is the following:

  1. For every element in V_dg, find the indices for dof_coords belonging to each element
  2. Compute the quality of the element using its dofs
  3. If the quality is low enough, go to indices, and remove the dof_coord indices belonging to the element

How can I do this?

One quality metric I can use is the radius ratio metric for tets:

import numpy as np


# Face areas
def tri_area(p1, p2, p3):
    return 0.5 * np.linalg.norm(np.cross(p2-p1, p3-p1))


def tetrahedron_quality(nodes):
    """
    Compute the quality of a tetrahedron using radius ratio.
    
    Quality = 3 * (inradius / circumradius)
    Perfect regular tetrahedron has quality = 1.0.
    
    Parameters:
        nodes: array-like of shape (4, 3) - 4 vertices [x, y, z]
    
    Returns:
        float: Quality value (1.0 or higher), where 1.0 is perfect.
    """
    A, B, C, D = [np.array(p, dtype=float) for p in nodes]
    

    # Volume
    V = np.abs(np.dot(B-A, np.cross(C-A, D-A))) / 6.0
    
    if V < 1e-12:
        return 0.0
    
    
    
    # Inradius
    surf_total = tri_area(A,B,C) + tri_area(A,C,D) + tri_area(A,B,D) + tri_area(B,C,D)
    inradius   = 3 * V / surf_total
    
    
        
    # Circumradius using edge lengths  
    a  = np.linalg.norm(A - B)
    b  = np.linalg.norm(A - C)
    c  = np.linalg.norm(A - D)
    A2 = np.linalg.norm(C - D)
    B2 = np.linalg.norm(B - D)
    C2 = np.linalg.norm(B - C)
    
    v1 =  a*A2 + b*B2 + c*C2
    v2 =  a*A2 + b*B2 - c*C2
    v3 =  a*A2 - b*B2 + c*C2
    v4 = -a*A2 + b*B2 + c*C2
    
    circumradius = np.sqrt(v1*v2*v3*v4)/(24*V)
    
    if circumradius < 1e-12:
        return 0.0
    
    
    return 3 * inradius / circumradius





# Regular tetrahedron: quality should be 1
nodes = [[ 1,  0,  1/np.sqrt(2)],
         [-1,  0,  1/np.sqrt(2)],
         [ 0,  1, -1/np.sqrt(2)],
         [ 0, -1, -1/np.sqrt(2)]]
       
q = tetrahedron_quality(nodes)
print(f"Tetrahedron Quality: {q:.6f}")

Thank you,
Alex

4

I think a difficulty here is that dofs of the DG0 space are located in the cell midpoints, so figuring out the mesh quality related to each dof is not directly easy. However, the dofs of the DG0 space are sorted in the same ordering as the mesh cells, so you could compute mesh quality and distance from the sphere center using cell vertices accessed like this:

N = 4
mesh = dfx.mesh.create_unit_square(MPI.COMM_WORLD, N, N)

cells = mesh.geometry.dofmap
cell_vertices = mesh.geometry.x[cells] # shape (num_cells, num_verts_in_cell, 3)
cell_midpoints = np.mean(cell_vertices, axis=1) # shape (num_cells, 3)

DG0 = dfx.fem.functionspace(mesh, ("DG", 0))
DG0_dof_locations = DG0.tabulate_dof_coordinates() # shape (num_cells, 3)

print(f"{np.allclose(cell_midpoints, DG0_dof_locations) = }")
# np.allclose(cell_midpoints, DG0_dof_locations) = True

So you could do the following:

  1. Compute principal stress and store in a DG0 function.
  2. Compute mesh quality using cell_vertices.
  3. Compute distance from sphere center using cell_midpoints.
  4. Find minimum of principal stress while filtering with results of 2. and 3.

Hope this is helpful.