G2 Element 10

G2Element10.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 kind.
# ------------------------------------------------------------------------
# Element 10: Rigid element with one hysteretic uniaxial material at each end

from G2MaterialBilinear import *
from G2MaterialPlasticity import *
from G2MaterialBoucWen import *

class element10():

    # -------------------------------------------------
    # Constructor
    # -------------------------------------------------
    def __init__(self, section, material, elno):

        self.type = "element10"

        # Store input
        self.no = elno
        self.EA = float(section[1])

        # Create materials
        if material[0] == 'Bilinear':
            self.hinge1 = bilinearMaterial(material)
            self.hinge2 = bilinearMaterial(material)
        elif material[0] == 'Plasticity':
            self.hinge1 = plasticityMaterial(material)
            self.hinge2 = plasticityMaterial(material)
        elif material[0] == 'BoucWen':
            self.hinge1 = boucWenMaterial(material)
            self.hinge2 = boucWenMaterial(material)
        else:
            print("Error: Wrong material given to element type 10")
            import sys
            sys.exit()

        # Store basic forces to ease getResponse
        self.FbTrial = []
        self.FbCommitted = []

    # -------------------------------------------------
    # Set parameter
    # -------------------------------------------------
    def setParameter(self, parameter, value):

        if parameter == 'EA':
            self.EA = value
        else:
            self.hinge1.setParameter(parameter, value)
            self.hinge2.setParameter(parameter, value)

    # -------------------------------------------------
    # Initialize
    # -------------------------------------------------
    def initialize(self, xyz):

        # Geometry restricted to 1,2 plane (x,y)
        dx = xyz[1,:] - xyz[0,:]
        L = np.sqrt(dx.dot(dx))

        # Check length
        if L < 1e-10:
            print('\n'"ERROR: Element number", self.no, "has zero length")
            import sys
            sys.exit()

    # -------------------------------------------------
    # State determination
    # -------------------------------------------------
    def state(self, xyz, ug, theLambda):

        # Direction cosines
        dx = xyz[1,:] - xyz[0,:]
        L = np.sqrt(dx.dot(dx))
        dx = dx / L

        # Transformation matrix from basic to global
        Tbg = np.array([[-dx[1]/L, dx[0]/L, 1.0, dx[1]/L, -dx[0]/L, 0.0],
                        [-dx[1]/L, dx[0]/L, 0.0, dx[1]/L, -dx[0]/L, 1.0],
                        [-dx[0],  -dx[1],   0.0, dx[0],    dx[1],   0.0]])

        # Basic deformations (rot-left, rot-right, elongation)
        ub = Tbg.dot(ug)

        # State determination for hinges
        moment1, stiffness1 = self.hinge1.state(ub[0,:])
        moment2, stiffness2 = self.hinge2.state(ub[1,:])

        # Axial portion
        axialKb = self.EA / L
        axialFb = axialKb * ub[2, 0]

        # Basic force and stiffness
        Fb = np.array([moment1, moment2, axialFb])
        Kb = np.diag([stiffness1, stiffness2, axialKb])

        # Store basic forces to avoid recomputing in getResponse
        self.FbTrial = Fb

        # Global restoring force and stiffness matrix
        Fg = np.transpose(Tbg).dot(Fb)
        Kg = np.einsum('ji,jk,kl->il', Tbg, Kb, Tbg)

        return Fg, Kg

    # -------------------------------------------------
    # DDM sensitivity analysis
    # -------------------------------------------------
    def stateDerivative(self, xyz, ug, theLambda, ddmParameter, ddmIndex, ddmIsHere, dKflag='none'):

        dx = xyz[1,:] - xyz[0,:]
        L = np.sqrt(dx.dot(dx))
        dx = dx / L

        Tbg = np.array([[-dx[1]/L, dx[0]/L, 1.0, dx[1]/L, -dx[0]/L, 0.0],
                        [-dx[1]/L, dx[0]/L, 0.0, dx[1]/L, -dx[0]/L, 1.0],
                        [-dx[0],  -dx[1],   0.0, dx[0],    dx[1],   0.0]])

        ub = Tbg.dot(ug)

        dmoment1, dstiffness1, dKdu1 = self.hinge1.stateDerivative(ub[0,:], ddmParameter, ddmIndex, ddmIsHere, dKflag)
        dmoment2, dstiffness2, dKdu2 = self.hinge2.stateDerivative(ub[1,:], ddmParameter, ddmIndex, ddmIsHere, dKflag)

        dFb = np.array([dmoment1, dmoment2, 0.0])
        dKb = np.diag([dstiffness1, dstiffness2, 0.0])
        dKbdub = np.array([[[dKdu1, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]],
                           [[0.0, 0.0, 0.0],[0.0, dKdu2, 0.0],[0.0, 0.0, 0.0]],
                           [[0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]]])

        dFg = (np.transpose(Tbg)).dot(dFb)
        dKg = np.einsum('ji,jk,kl->il', Tbg, dKb, Tbg)
        dKgdug = np.einsum('ki,knl,nm,lj->ijm', Tbg, dKbdub, Tbg, Tbg)

        return dFg, dKg, dKgdug

    # -------------------------------------------------
    # Commit
    # -------------------------------------------------
    def commit(self, xyz, u):

        self.FbCommitted = self.FbTrial
        self.hinge1.commit()
        self.hinge2.commit()

    # -------------------------------------------------
    # Commit sensitivity history variables
    # -------------------------------------------------
    def commitSensitivity(self, xyz, ug, ddmug, ddmParameter, ddmIndex, ddmIsHere):

        dx = xyz[1,:] - xyz[0,:]
        L = np.sqrt(dx.dot(dx))
        dx = dx / L

        Tbg = np.array([[-dx[1]/L, dx[0]/L, 1.0, dx[1]/L, -dx[0]/L, 0.0],
                        [-dx[1]/L, dx[0]/L, 0.0, dx[1]/L, -dx[0]/L, 1.0],
                        [-dx[0],  -dx[1],   0.0, dx[0],    dx[1],   0.0]])

        ub = Tbg.dot(ug)
        ddmub = Tbg.dot(ddmug)

        self.hinge1.commitSensitivity(ub[0,:], ddmub[0], ddmParameter, ddmIndex, ddmIsHere)
        self.hinge2.commitSensitivity(ub[1,:], ddmub[1], ddmParameter, ddmIndex, ddmIsHere)

    # -------------------------------------------------
    # Print response
    # -------------------------------------------------
    def getResponse(self, xyz):

        # Print element force
        print("Element %d (%s): No information " % (self.no, self.type))