SolidCrossSectionAnalysis

SolidCrossSectionAnalysis.py

# ------------------------------------------------------------------------
# The following Python code is implemented by Professor Terje Haukaas at
# the University of British Columbia in Vancouver, Canada. It is made
# freely available online at terje.civil.ubc.ca together with notes,
# examples, and additional Python code. Please be cautious when using
# this code; it may contain bugs and comes without warranty of any form.
# ------------------------------------------------------------------------

import matplotlib.pyplot as plt
import numpy as np
import terjesfunctions as fcns
from matplotlib.tri import Triangulation, UniformTriRefiner


# ------------------------------------------------------------------------
# FUNCTION THAT HELPS CREATE THE MESH
# ------------------------------------------------------------------------
def createRectangle(yLowerLeft, zLowerLeft, width, height):
    if width > height:
        intervals = round(width / float(height))
        forward = np.linspace(yLowerLeft, yLowerLeft+width, intervals+1)
        backward = np.linspace(yLowerLeft+width, yLowerLeft, intervals+1)
        yValues = np.concatenate((forward, backward))
        forward = np.linspace(zLowerLeft, zLowerLeft, intervals+1)
        backward = np.linspace(zLowerLeft+height, zLowerLeft+height, intervals+1)
        zValues = np.concatenate((forward, backward))
    else:
        intervals = round(height / float(width))
        up = np.linspace(yLowerLeft+width, yLowerLeft+width, intervals+1)
        down = np.linspace(yLowerLeft, yLowerLeft, intervals+1)
        yValues = np.concatenate((up, down))
        up = np.linspace(zLowerLeft, zLowerLeft+height, intervals+1)
        down = np.linspace(zLowerLeft+height, zLowerLeft, intervals+1)
        zValues = np.concatenate((up, down))
    return yValues, zValues


# ------------------------------------------------------------------------
# SET MESH REFINEMENT
# ------------------------------------------------------------------------

refinement = 2


# ------------------------------------------------------------------------
# INPUT, CONVEX AND NUMBERED COUNTERCLOCKWISE (The last one is used by the code)
# ------------------------------------------------------------------------

# L-SHAPE
t = 10.0
b = 100.0
h = 200.0
yBottom, zBottom = createRectangle(t/2, -t/2, b-t/2, t)
yMiddle, zMiddle = createRectangle(-t/2, -t/2, t, t)
yWeb, zWeb = createRectangle(-t/2, t/2, t, h+t/2)
y = [yBottom, yMiddle, yWeb]
z = [zBottom, zMiddle, zWeb]

# ------------------------------------------------------------------------

# I-BEAM
ttop = 15.0
tbottom = 15.0
tweb = 15.0
btop = 100.0
bbottom = 150.0
h = 250.0
yBottomLeft, zBottomLeft = createRectangle(-bbottom/2, -tbottom/2, bbottom/2-tweb/2, tbottom)
yBottomRight, zBottomRight = createRectangle(tweb/2, -tbottom/2, bbottom/2-tweb/2, tbottom)
yBottomMiddle, zBottomMiddle = createRectangle(-tweb/2, -tbottom/2, tweb, tbottom)
yWeb, zWeb = createRectangle(-tweb/2, tbottom/2, tweb, h-tbottom/2-ttop/2)
yTopLeft, zTopLeft = createRectangle(-btop/2, h-ttop/2, btop/2-tweb/2, ttop)
yTopRight, zTopRight = createRectangle(tweb/2, h-ttop/2, btop/2-tweb/2, ttop)
yTopMiddle, zTopMiddle = createRectangle(-tweb/2, h-ttop/2, tweb, ttop)
y = [yBottomLeft, yBottomMiddle, yBottomRight, yWeb, yTopLeft, yTopMiddle, yTopRight]
z = [zBottomLeft, zBottomMiddle, zBottomRight, zWeb, zTopLeft, zTopMiddle, zTopRight]

# ------------------------------------------------------------------------

# ONE CELL
b = 500.0
h = 200.0
ttop = 10.0
tbottom = 10.0
tleft = 15.0
tright = 20.0
yBottom, zBottom = createRectangle(tleft/2, -tbottom/2, b-tleft/2-tright/2, tbottom)
yLowerLeft, zLowerLeft = createRectangle(-tleft/2, -tbottom/2, tleft, tbottom)
yLeft, zLeft = createRectangle(-tleft/2, tbottom/2, tleft, h-tbottom/2-ttop/2)
yUpperLeft, zUpperLeft = createRectangle(-tleft/2, h-ttop/2, tleft, ttop)
yTop, zTop = createRectangle(tleft/2, h-ttop/2, b-tleft/2-tright/2, ttop)
yUpperRight, zUpperRight = createRectangle(b-tright/2, h-ttop/2, tright, ttop)
yRight, zRight = createRectangle(b-tright/2, tbottom/2, tright, h-ttop/2-tbottom/2)
yLowerRight, zLowerRight = createRectangle(b-tright/2, -tbottom/2, tright, tbottom)
y = [yBottom, yLowerLeft, yLeft, yUpperLeft, yTop, yUpperRight, yRight, yLowerRight]
z = [zBottom, zLowerLeft, zLeft, zUpperLeft, zTop, zUpperRight, zRight, zLowerRight]

# ------------------------------------------------------------------------

# TRIANGLE
b = 100.0
h = 0.5*np.sqrt(3)*b
y = np.array([(-b/2, b/2, 0.0)])
z = np.array([(0.0, 0.0, h)])

# ------------------------------------------------------------------------

# RECTANGLE
b = 100.0
h = 200.0
yRect, zRect = createRectangle(0.0, 0.0, b, h)
y = [yRect]
z = [zRect]

# ------------------------------------------------------------------------

# CHANNEL PROFILE
t = 20.0
b = 80.0
h = 200.0
yBottomRight, zBottomRight = createRectangle(t/2, -t/2, b-t/2, t)
yBottomMiddle, zBottomMiddle = createRectangle(-t/2, -t/2, t, t)
yWeb, zWeb = createRectangle(-t/2, t/2, t, h-t/2-t/2)
yTopRight, zTopRight = createRectangle(t/2, h-t/2, b-t/2, t)
yTopMiddle, zTopMiddle = createRectangle(-t/2, h-t/2, t, t)
y = [yBottomMiddle, yBottomRight, yWeb, yTopMiddle, yTopRight]
z = [zBottomMiddle, zBottomRight, zWeb, zTopMiddle, zTopRight]


# ------------------------------------------------------------------------
# CHECK INPUT
# ------------------------------------------------------------------------

if len(y) != len(z):
    print("ERROR: Inconsistent number of shapes!")
for i in range(len(y)):
    if len(y[i]) != len(z[i]):
        print("ERROR: Inconsistent number of coordinates!")


# ------------------------------------------------------------------------
# CREATE TRIANGLES & FINITE ELEMENTS
# ------------------------------------------------------------------------
nodes = []
areas = []
localyCentroid = []
localzCentroid = []
connectivity = []
threeNodes = [0, 0, 0]
A = 0
sumyA = 0
sumzA = 0
sumySquaredA = 0
sumzSquaredA = 0
sumyzA = 0
for i in range(len(y)):

    # Create the triangles
    triang = Triangulation(y[i], z[i])
    refiner = UniformTriRefiner(triang)
    newtriang = refiner.refine_triangulation(subdiv=refinement)

    # Store nodes and connectivity
    for t in newtriang.triangles:

        # Triangle vertices
        n1 = t[0]
        n2 = t[1]
        n3 = t[2]

        # y-coordinates of vertices
        y1 = newtriang.x[n1]
        y2 = newtriang.x[n2]
        y3 = newtriang.x[n3]

        # y-coordinates of vertices
        z1 = newtriang.y[n1]
        z2 = newtriang.y[n2]
        z3 = newtriang.y[n3]

        # Centroid
        y0 = (y1 + y2 + y3) / 3.0
        z0 = (z1 + z2 + z3) / 3.0
        localyCentroid.append(y0)
        localzCentroid.append(z0)

        # Area
        dA = 0.5 * ((y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1))
        areas.append(dA)
        A += dA
        if dA < 0.0:
            print("ERROR: NODES ARE NUMBERED CLOCKWISE!")

        # First moment of the area about the origin
        sumyA += y0 * dA
        sumzA += z0 * dA

        # Prepare for calculation of the moment of inertia
        sumySquaredA += dA / 6.0 * (y1**2 + y2**2 + y3**2 + y1*y2 + y1*y3 + y2*y3)
        sumzSquaredA += dA / 6.0 * (z1**2 + z2**2 + z3**2 + z1*z2 + z1*z3 + z2*z3)

        # Linear shape functions gave that; constant coordinates value is simpler but less accurate:
        #sumySquaredA += y0 ** 2 * dA
        #sumzSquaredA += z0 ** 2 * dA

        # Prepare for calculation of the product of inertia, using constant coordinate values
        sumyzA += y0 * z0 * dA

        # Initialize the nodes array, if necessary
        if len(nodes) == 0:
            nodes.append([y1, z1])

        # Store nodal coordinates, if they exist
        node1Exists = False
        node2Exists = False
        node3Exists = False
        for j in range(len(nodes)):

            if abs((nodes[j])[0]-y1)<1e-9 and abs((nodes[j])[1]-z1)<1e-9:
                threeNodes[0] = int(j)
                node1Exists = True
            if abs((nodes[j])[0]-y2)<1e-9 and abs((nodes[j])[1]-z2)<1e-9:
                threeNodes[1] = int(j)
                node2Exists = True
            if abs((nodes[j])[0]-y3)<1e-9 and abs((nodes[j])[1]-z3)<1e-9:
                threeNodes[2] = int(j)
                node3Exists = True

        # Store NEW nodal coordinates
        if node1Exists == False:
            nodes.append([y1, z1])
            threeNodes[0] = int(len(nodes)-1)
        if node2Exists == False:
            nodes.append([y2, z2])
            threeNodes[1] = int(len(nodes)-1)
        if node3Exists == False:
            nodes.append([y3, z3])
            threeNodes[2] = int(len(nodes)-1)

        # Store the connectivity
        connectivity.append([threeNodes[0], threeNodes[1], threeNodes[2]])

# Count DOFs and allocate K and F
numDOFsPerNodeInStructure = 1
numDOFsPerNodeInElement = 1
numNodes = len(nodes)
totalNumDOFs = numDOFsPerNodeInStructure * numNodes
Ka = np.zeros((totalNumDOFs, totalNumDOFs))
Fa = np.zeros(totalNumDOFs)


# ------------------------------------------------------------------------
# MOMENTS OF INERTIA and PRINCIPAL AXES ORIENTATION
# ------------------------------------------------------------------------

# Centroid coordinates
yCentroid = sumyA / A
zCentroid = sumzA / A

# Moments of inertia
Iy = sumzSquaredA - 2.0 * zCentroid * sumzA + zCentroid**2 * A
Iz = sumySquaredA - 2.0 * yCentroid * sumyA + yCentroid**2 * A

# Product of inertia
Iyz = sumyzA + yCentroid * zCentroid * A - zCentroid * sumyA - yCentroid * sumzA

# Apply caution with the principal axes rotation: what if Iy=Iz
thetaPrincipal = 0.0
if np.abs(Iyz) < 1e-3:
    thetaPrincipal = 0.0
elif np.abs(Iy-Iz) < 1e-3:
    if np.abs(Iyz) < 1e-3:
        thetaPrincipal = 0.0
    else:
        thetaPrincipal = np.sign(Iyz) * 0.25 * np.pi
else:
    thetaPrincipal = 0.5 * np.arctan(2.0*Iyz/(Iy-Iz))

# Let counter-clockwise rotation of the system be positive
thetaPrincipal *= -1

# Transform the moments of inertia
IyRotated =  0.5*(Iy+Iz) + 0.5*np.sqrt((Iz-Iy)**2 + 4.0* Iyz**2)
IzRotated =  0.5*(Iy+Iz) - 0.5*np.sqrt((Iz-Iy)**2 + 4.0* Iyz**2)


# ------------------------------------------------------------------------
# ROTATE AXES IF Iyz HAS VALUE
# ------------------------------------------------------------------------

# Shift and rotate axis system to align with centroid and principal axes
if np.abs(thetaPrincipal) > 1e-3:
    for i in range(len(nodes)):
        yOld = (nodes[i])[0] - yCentroid
        zOld = (nodes[i])[1] - zCentroid
        (nodes[i])[0] = yOld * np.cos(thetaPrincipal) + zOld * np.sin(thetaPrincipal)
        (nodes[i])[1] = -yOld * np.sin(thetaPrincipal) + zOld * np.cos(thetaPrincipal)
    storeYcentroid = yCentroid
    storeZcentroid = zCentroid
    yCentroid = 0.0
    zCentroid = 0.0


# ------------------------------------------------------------------------
# ESTABLISH THE OMEGA DIAGRAM
# ------------------------------------------------------------------------

for i in range(len(connectivity)):

    # Node numbers
    node1 = (connectivity[i])[0]
    node2 = (connectivity[i])[1]
    node3 = (connectivity[i])[2]

    # y-coordinates of vertices
    y1 = (nodes[node1])[0] - yCentroid
    y2 = (nodes[node2])[0] - yCentroid
    y3 = (nodes[node3])[0] - yCentroid

    # y-coordinates of vertices
    z1 = (nodes[node1])[1] - zCentroid
    z2 = (nodes[node2])[1] - zCentroid
    z3 = (nodes[node3])[1] - zCentroid

    # Coordinates for finite element analysis
    z23 = z2-z3
    z31 = z3-z1
    z12 = z1-z2
    y32 = y3-y2
    y13 = y1-y3
    y21 = y2-y1

    # Element stiffness matrix
    factor = 1.0/(4.0*areas[i])
    k11 = factor * (y32**2 + z23**2)
    k22 = factor * (y13**2 + z31**2)
    k33 = factor * (y21**2 + z12**2)
    k12 = factor * (y13*y32 + z23*z31)
    k13 = factor * (y21*y32 + z12*z23)
    k23 = factor * (y13*y21 + z12*z31)
    Kelement = np.array([(k11, k12, k13),
                         (k12, k22, k23),
                         (k13, k23, k33)])

    # Element load vector
    factor = 1.0/6.0
    F1 = factor * ((y1*y32 - z1*z23) + (y2*y32 - z2*z23) + (y3*y32 - z3*z23))
    F2 = factor * ((y1*y13 - z1*z31) + (y2*y13 - z2*z31) + (y3*y13 - z3*z31))
    F3 = factor * ((y1*y21 - z1*z12) + (y2*y21 - z2*z12) + (y3*y21 - z3*z12))
    Felement = np.array([-F1, -F2, -F3])

    # Element assembly
    nodeList = np.array([node1+1, node2+1, node3+1])
    Ka = fcns.assembleStiffnessMatrix(Ka, Kelement, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)
    Fa = fcns.assembleLoadVector(Fa, Felement, nodeList, numDOFsPerNodeInElement, numDOFsPerNodeInStructure)

# Lock the warping at one node and solve for the others
Ff = Fa[0:totalNumDOFs-1]
for i in range(totalNumDOFs-1):
    Ff[i] = Ff[i] - Ka[i, totalNumDOFs-1]
Kf = Ka[0:totalNumDOFs-1, 0:totalNumDOFs-1]
uf = np.linalg.solve(Kf, Ff)
ua = np.append(uf, 1.0)

# Determine shear centre coordinates
warpIntegral = 0
yscIntegral = 0
zscIntegral = 0
for i in range(len(connectivity)):

    # Node numbers
    node1 = (connectivity[i])[0]
    node2 = (connectivity[i])[1]
    node3 = (connectivity[i])[2]

    # y-coordinates of vertices
    y1 = (nodes[node1])[0] - yCentroid
    y2 = (nodes[node2])[0] - yCentroid
    y3 = (nodes[node3])[0] - yCentroid

    # y-coordinates of vertices
    z1 = (nodes[node1])[1] - zCentroid
    z2 = (nodes[node2])[1] - zCentroid
    z3 = (nodes[node3])[1] - zCentroid

    # Warping values
    u1 = ua[node1]
    u2 = ua[node2]
    u3 = ua[node3]

    # Accumulate integrals needed to get final omega diagram
    warpIntegral += (u1+u2+u3)/3.0 * areas[i]
    yscIntegral += areas[i]/12.0*(u1*(2*z1+z2+z3) + u2*(z1+2*z2+z3) + u3*(z1+z2+2*z3))
    zscIntegral += areas[i]/12.0*(u1*(2*y1+y2+y3) + u2*(y1+2*y2+y3) + u3*(y1+y2+2*y3))

# Normalizing constant and shear centre coordinates
normalizingConstant = -warpIntegral / A
if np.abs(thetaPrincipal) > 1e-3:
    ysc = -yscIntegral / IyRotated
    zsc = zscIntegral / IzRotated
else:
    ysc = -yscIntegral / Iy
    zsc = zscIntegral / Iz

# Final warping diagram
finalWarp = []
for i in range(totalNumDOFs):

    yValue = (nodes[i])[0]
    zValue = (nodes[i])[1]
    finalWarp.append(ua[i] + normalizingConstant + ysc*(zValue-zCentroid) - zsc*(yValue-yCentroid))

# Shear centre coordinates in original axis system
ysc = yCentroid + ysc
zsc = zCentroid + zsc


# ------------------------------------------------------------------------
# WARPING CONSTANT and ST. VENANT CONSTANT
# ------------------------------------------------------------------------

Cw = 0
J = 0
for i in range(len(connectivity)):

    # Node numbers
    node1 = (connectivity[i])[0]
    node2 = (connectivity[i])[1]
    node3 = (connectivity[i])[2]

    # y-coordinates of vertices
    y1 = (nodes[node1])[0] - yCentroid
    y2 = (nodes[node2])[0] - yCentroid
    y3 = (nodes[node3])[0] - yCentroid

    # y-coordinates of vertices
    z1 = (nodes[node1])[1] - zCentroid
    z2 = (nodes[node2])[1] - zCentroid
    z3 = (nodes[node3])[1] - zCentroid

    # Coordinates for finite element analysis
    z23 = z2-z3
    z31 = z3-z1
    z12 = z1-z2
    y32 = y3-y2
    y13 = y1-y3
    y21 = y2-y1

    # Warping values
    u1 = finalWarp[node1]
    u2 = finalWarp[node2]
    u3 = finalWarp[node3]

    # Warping constant
    dA = areas[i]
    Cw += dA / 6.0 * (u1**2 + u2**2 + u3**2 + u1*u2 + u1*u3 + u2*u3)

    # St. Venant constant
    Czeta1  = (u2 * y1 * y13 + u3 * y1 * y21 + u1 * y1 * y32 - u3 * z1 * z12 - u1 * z1 * z23 - u2 * z1 * z31)/(2*dA)
    Czeta2  = (u2 * y13 * y2 + u3 * y2 * y21 + u1 * y2 * y32 - u3 * z12 * z2 - u1 * z2 * z23 - u2 * z2 * z31)/(2*dA)
    Czeta3  = (u2 * y13 * y3 + u3 * y21 * y3 + u1 * y3 * y32 - u3 * z12 * z3 - u1 * z23 * z3 - u2 * z3 * z31)/(2*dA)
    Czeta12 = 2 * y1 * y2 + 2 * z1 * z2
    Czeta13 = 2 * y1 * y3 + 2 * z1 * z3
    Czeta23 = 2 * y2 * y3 + 2 * z2 * z3
    Czeta1s = y1**2 + z1**2
    Czeta2s = y2**2 + z2**2
    Czeta3s = y3**2 + z3**2
    J += (Czeta1+Czeta2+Czeta3)*dA/3.0 + (Czeta12+Czeta13+Czeta23)*dA/12.0 + (Czeta1s+Czeta2s+Czeta3s)*dA/6.0


# ------------------------------------------------------------------------
# ROTATE AXES BACK
# ------------------------------------------------------------------------

if np.abs(thetaPrincipal) > 1e-3:
    yCentroid = storeYcentroid
    zCentroid = storeZcentroid
    yscOff = ysc
    zscOff = zsc
    ysc = yscOff * np.cos(thetaPrincipal) - zscOff * np.sin(thetaPrincipal) + yCentroid
    zsc = yscOff * np.sin(thetaPrincipal) + zscOff * np.cos(thetaPrincipal) + zCentroid

    for i in range(len(nodes)):
        yOld = (nodes[i])[0]
        zOld = (nodes[i])[1]
        (nodes[i])[0] = yOld * np.cos(thetaPrincipal) - zOld * np.sin(thetaPrincipal) + yCentroid
        (nodes[i])[1] = yOld * np.sin(thetaPrincipal) + zOld * np.cos(thetaPrincipal) + zCentroid


# ------------------------------------------------------------------------
# PRINT RESULTS
# ------------------------------------------------------------------------

print('\n'"Area:", A)
print('\n'"Centroid y-coordinate: %.2f" % yCentroid)
print('\n'"Centroid z-coordinate: %.2f" % zCentroid)
print('\n'"Moment of inertia Iy: %.2f" % Iy)
print('\n'"Moment of inertia Iz: %.2f" % Iz)
if np.abs(thetaPrincipal) > 1e-3:
    print('\n'"Principal axis rotation (degrees): %.2f" % (thetaPrincipal/np.pi*180.0))
    print('\n'"Product of inertia Iyz: %.2f" % Iyz)
    print('\n'"New principal moment of inertia Iy: %.2f" % IyRotated)
    print('\n'"New principal moment of inertia Iz: %.2f" % IzRotated)
    print('\n'"Shear centre y-coordinate in original coordinate system: %.2f" % ysc)
    print('\n'"Shear centre z-coordinate in original coordinate system: %.2f" % zsc)
else:
    print('\n'"Shear centre y-coordinate: %.2f" % ysc)
    print('\n'"Shear centre z-coordinate: %.2f" % zsc)
print('\n'"St. Venant torsion constant J: %.2f" % J)
print('\n'"Warping torsion constant Cw: %.2f" % Cw)

# ------------------------------------------------------------------------
# PLOT
# ------------------------------------------------------------------------

# Determine warping bounds
maxWarp = np.amax(finalWarp)
minWarp = np.amin(finalWarp)

# Plot the edges of the triangles
alreadyPlotted = []
for i in range(len(connectivity)):

    # Node numbers
    node1 = (connectivity[i])[0]
    node2 = (connectivity[i])[1]
    node3 = (connectivity[i])[2]

    # y-coordinates of vertices
    y1 = (nodes[node1])[0]
    y2 = (nodes[node2])[0]
    y3 = (nodes[node3])[0]

    # y-coordinates of vertices
    z1 = (nodes[node1])[1]
    z2 = (nodes[node2])[1]
    z3 = (nodes[node3])[1]

    # Warping value at nodes
    warp1 = finalWarp[node1]
    warp2 = finalWarp[node2]
    warp3 = finalWarp[node3]

    # Plot first edge
    if [node1, node2] not in alreadyPlotted:
        warpAverage = (warp1+warp2)/2.0
        colorCode = (warpAverage - minWarp) / (maxWarp - minWarp)
        plt.plot(np.array([y1, y2]), np.array([z1, z2]), color=(colorCode, 0, 1-colorCode), zorder=1)
        alreadyPlotted.append([node1, node2])

    # Plot second edge
    if [node1, node3] not in alreadyPlotted:
        warpAverage = (warp1+warp3)/2.0
        colorCode = (warpAverage - minWarp) / (maxWarp - minWarp)
        plt.plot(np.array([y1, y3]), np.array([z1, z3]), color=(colorCode, 0, 1-colorCode), zorder=1)
        alreadyPlotted.append([node1, node3])

    # Plot third edge
    if [node2, node3] not in alreadyPlotted:
        warpAverage = (warp2+warp3)/2.0
        colorCode = (warpAverage - minWarp) / (maxWarp - minWarp)
        plt.plot(np.array([y2, y3]), np.array([z2, z3]), color=(colorCode, 0, 1-colorCode), zorder=1)
        alreadyPlotted.append([node2, node3])

# Create the plot
max = 0
min = 0
for i in range(len(y)):
    theMax = np.max(y[i])
    theMin = np.min(y[i])
    if theMax > max:
        max = theMax
    if theMin < min:
        min = theMin
for i in range(len(z)):
    theMax = np.max(z[i])
    theMin = np.min(z[i])
    if theMax > max:
        max = theMax
    if theMin < min:
        min = theMin
yMin = min - (max-min)*0.4
border = (max - min) * 0.1
plt.plot(np.array([yCentroid, max+border]), np.array([zCentroid, zCentroid+(np.abs(max)+border-yCentroid)*np.tan(thetaPrincipal)]), 'k-',linewidth=1.0)
plt.plot(np.array([yCentroid, yMin-border]), np.array([zCentroid, zCentroid-(np.abs(yMin)+border+yCentroid)*np.tan(thetaPrincipal)]), 'k-',linewidth=1.0)
plt.plot(np.array([yCentroid, yCentroid-(np.abs(max)+border-zCentroid)*np.tan(thetaPrincipal)]), np.array([zCentroid, max+border]), 'k-',linewidth=1.0)
plt.plot(np.array([yCentroid, yCentroid+(np.abs(min)+border+zCentroid)*np.tan(thetaPrincipal)]), np.array([zCentroid, min-border]), 'k-',linewidth=1.0)
plt.scatter([yCentroid], [zCentroid], s=100, c='k', marker='o', zorder=2)
plt.scatter([ysc], [zsc], s=100, c='k', marker='o', zorder=2)
plt.axis([yMin-border, max+border, min-border, max+border])
plt.title("Centroid + Shear Centre + Warping")
plt.show()