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.