How to plot ratios of stress tensor components?

dolfinx
2

You can for instance use:

  1. Locate which cell your point is in
  2. Pull it back to the reference element
  3. Use dolfinx.fem.Expression to evaluate the expression P[0,0]/P[1,1] at the point of interest.

Most of the methodology is shown in:

and can thus be added to the demo above as:

def compute_expression_value(expr, points):
    # Determine what process owns a point and what cells it lies within
    mesh = V.mesh
    collision_data = dolfinx.cpp.geometry.determine_point_ownership(
        mesh._cpp_object, points, 1e-6
    )
    owning_points = collision_data.dest_points
    cells = collision_data.dest_cells

    # Pull owning points back to reference cell
    mesh_nodes = mesh.geometry.x
    cmap = mesh.geometry.cmap
    ref_x = np.zeros((len(cells), mesh.geometry.dim), dtype=mesh.geometry.x.dtype)
    for i, (point, cell) in enumerate(zip(owning_points, cells)):
        geom_dofs = mesh.geometry.dofmap[cell]
        ref_x[i] = cmap.pull_back(point.reshape(-1, 3), mesh_nodes[geom_dofs])

    import math

    expr_size = math.prod(expr.ufl_shape)
    if len(cells) > 0:
        # NOTE: Expression lives on only this communicator rank
        compiled_expr = dolfinx.fem.Expression(expr, ref_x, comm=MPI.COMM_SELF)
        values = compiled_expr.eval(mesh, np.asarray(cells, dtype=np.int32))
        # Extract values for each (cell, point) pair
        stripped_values = np.diagonal(
            values.reshape(
                len(cells),
                ref_x.shape[0],
                expr_size,
            ),
            axis1=0,
            axis2=1,
        ).T.reshape(len(cells), *expr.ufl_shape)
        return stripped_values
    else:
        return np.empty((0, *expr.ufl_shape))


# Example to illustrate that this works for any expression of any shape
# is commented out here
# P_rate = ufl.as_tensor([[x[0], x[1]], [x[2], -x[0]]])

P_rate = P[0, 0] / P[1, 1]
x = ufl.SpatialCoordinate(domain)
values = compute_expression_value(
    P_rate, np.array([[0.5, 0.2, 0.3], [1, 0, 0.1], [0.764, 0.33, 0.23]])
)
print(f"P_rate: {values}")