How to calculate the principal stretches along the deformed neutral axis of a deformed beam?

General
1

Suppose if have a unit cube. I apply the deformation U1 = Expression(("0.0", "cos(theta) - sin(theta)", "sin(theta) + cos(theta)"), theta = pi/3, degree=0) on the cube. Using Paraview, I would like to visualize and graph the principal stretches (\lambda_x, \lambda_y, and \lambda_z) along the deformed neutral axis of the deformed cube. By neutral axis, I mean a line passing through the center of the cube parallel to its length.

My actual problem deals with pure bending of a rod. I use a Neo-Hookean strain energy density function with a uniform pressure applied on the top of the beam, keeping the middle cross section of the beam fixed. I need to calculate the principal stretches along the deformed neutral axis of the bent beam.

Here is the MWE:

from dolfin import *     

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}
              
# Create mesh and define function space
mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 1.0)
                
U1 = Expression(("0.0",
                "cos(theta) - sin(theta)",
                "sin(theta) + cos(theta)"),
                scale = 50.0, y0 = 0.5, z0 = 0.5, theta = pi/3, degree=0)
                
disp = project(U1, V)

File("displacement.pvd")<< disp

How do I find the principal stretches along the deformed neutral axis so that I can visualize and graph them using Paraview? If I am using a deformation gradient tensor, then can we get the principal stretches along the deformed neutral axis from the elastic tensor?