Untitled

#3D principal stress determination
#Laoise Reynolds
#20400934

import math

# ----------------------------------------------------------------
# 1. Define the stress invariants based on the given tensor components.
# ----------------------------------------------------------------
# Given stress components:
s11 = 70.0
s22 = 150.0
s33 = 250.0
s12 = 100.0
s13 = -100.0
s23 = -150.0  

# Compute invariants:
I1 = s11 + s22 + s33
I2 = s11*s22 + s22*s33 + s33*s11 - (s12**2 + s13**2 + s23**2)
I3 = ( s11*(s22*s33 - s23**2) 
      - s12*(s12*s33 - s13*s23)
      + s13*(s12*s23 - s13*s22) )

print("Stress Invariants:")
print("I1 =", I1)
print("I2 =", I2)
print("I3 =", I3)

# ----------------------------------------------------------------
# 2. Define the polynomial and its derivative.
#    The characteristic polynomial is:
#       P(lambda) = lambda^3 - I1*lambda^2 + I2*lambda - I3 = 0.
# ----------------------------------------------------------------

# Coefficients of the cubic polynomial: [1, -I1, I2, -I3]
coeffs = [1.0, -I1, I2, -I3]

def poly(x, coeffs):
    """Evaluate polynomial with given coefficients at x.
    Coefficients are ordered from highest to constant term."""
    result = 0
    degree = len(coeffs) - 1
    for i, c in enumerate(coeffs):
        result += c * x**(degree - i)
    return result

def dpoly(x, coeffs):
    """Evaluate derivative of the polynomial at x."""
    result = 0
    degree = len(coeffs) - 1
    for i, c in enumerate(coeffs[:-1]):  # skip constant term
        power = degree - i
        result += power * c * x**(power - 1)
    return result

# ----------------------------------------------------------------
# 3. Newton-Raphson function for finding one root.
# ----------------------------------------------------------------
def newton_raphson_poly(initial_guess, coeffs, tol=1e-6, max_iter=100):
    x = initial_guess
    for i in range(max_iter):
        f_val = poly(x, coeffs)
        df_val = dpoly(x, coeffs)
        if df_val == 0:
            print("Zero derivative encountered!")
            break
        x_new = x - f_val / df_val
        if abs(x_new - x) < tol:
            return x_new, i+1
        x = x_new
    return x, max_iter

# ----------------------------------------------------------------
# 4. Deflation of the polynomial (synthetic division)
#    Given a root r, factor out (x - r) from the polynomial.
# ----------------------------------------------------------------
def deflate_poly(coeffs, root):
    n = len(coeffs)
    new_coeffs = [coeffs[0]]  # start with the leading coefficient
    for i in range(1, n):
        new_coeffs.append(coeffs[i] + root * new_coeffs[-1])
    # The last element is the remainder (should be ~0); remove it.
    return new_coeffs[:-1]

# ----------------------------------------------------------------
# 5. Find one root using Newton-Raphson.
#    Choose an initial guess; here we use 50. (You may try other values.)
# ----------------------------------------------------------------
initial_guess = 50.0
root1, iterations = newton_raphson_poly(initial_guess, coeffs)
print("\nNewton-Raphson on cubic:")
print("Found root 1: {:.6f} (in {} iterations)".format(root1, iterations))

# ----------------------------------------------------------------
# 6. Deflate the cubic to obtain a quadratic.
# ----------------------------------------------------------------
coeffs_quad = deflate_poly(coeffs, root1)
a, b, c = coeffs_quad  # Quadratic: a*x^2 + b*x + c = 0

# Solve the quadratic equation.
discriminant = b**2 - 4*a*c
if discriminant >= 0:
    root2 = (-b + math.sqrt(discriminant)) / (2*a)
    root3 = (-b - math.sqrt(discriminant)) / (2*a)
else:
    # In our case with a symmetric stress tensor the eigenvalues are real,
    # but we include the complex solution branch for completeness.
    root2 = (-b + complex(0, math.sqrt(-discriminant))) / (2*a)
    root3 = (-b - complex(0, math.sqrt(-discriminant))) / (2*a)

# ----------------------------------------------------------------
# 7. Display the principal stresses (eigenvalues).
# ----------------------------------------------------------------
# Gather roots and sort them (smallest first):
roots = [root1, root2, root3]
# For sorting, if any are complex, sort by real part.
roots.sort(key=lambda r: r.real if isinstance(r, complex) else r)

print("\nPrincipal Stresses (Eigenvalues):")
for i, r in enumerate(roots, start=1):
    if isinstance(r, complex):
        print("Principal Stress {}: {:.6f} + {:.6f}i".format(i, r.real, r.imag))
    else:
        print("Principal Stress {}: {:.6f}".format(i, r))