(More notes from Professor Terje Haukaas at terje.civil.ubc.ca.)

Appendix A — Bilinear Material

import numpy as np

class bilinearMaterial():

    # -------------------------------------------------
    # Constructor
    # -------------------------------------------------
    def __init__(self, mat):

        # Material and geometry properties
        self.E = float(mat[1])      # Modulus of elasticity
        self.fy = float(mat[2])     # Yield stress
        self.alpha = float(mat[3])  # Second slope stiffness is alpha*E

        # Trial instantiation variables
        self.trialStrain = 0.0
        self.trialStress = 0.0
        self.trialBackStress = 0.0
        self.trialYielding = False

        # Committed instantiation variables
        self.committedStrain = 0.0
        self.committedStress = 0.0
        self.committedBackStress = 0.0
        self.committedYielding = False

        # Trial sensitivity history variables
        self.trialStrainDerivative = [0.0]
        self.trialStressDerivative = [0.0]
        self.trialBackStressDerivative = [0.0]

        # Committed sensitivity history variables
        self.committedStrainDerivative = [0.0]
        self.committedStressDerivative = [0.0]
        self.committedBackStressDerivative = [0.0]

        # Trial second-order sensitivity history variables
        self.trialSecondStrainDerivative = [0.0]
        self.trialSecondStressDerivative = [0.0]
        self.trialSecondBackStressDerivative = [0.0]

        # Committed second-order sensitivity history variables
        self.committedSecondStrainDerivative = [0.0]
        self.committedSecondStressDerivative = [0.0]
        self.committedSecondBackStressDerivative = [0.0]

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

        if parameter == 'E':
            self.E = value
        elif parameter == 'fy':
            self.fy = value
        elif parameter == 'alpha':
            self.alpha = value
        else:
            print("Cannot set parameter", parameter, "in material")

    # -------------------------------------------------
    # State determination
    # -------------------------------------------------
    def stateDetermination(self, eps):

        # This is an incremental material
        deps = eps[1]

        # Check for unloading
        unloading = False
        if (self.committedYielding == True and (self.committedStress - self.committedBackStress) * deps < 0):
            unloading = True

        # Check if the last state was elastic or if the strain increment implies unloading from yielding
        if self.committedYielding is False or unloading:

            # Strain increment that would cause yielding
            depsToYield = (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E

            # Keep elastic state handy, in case that becomes the conclusion
            self.trialStress = self.committedStress + self.E * deps
            Et = self.E
            self.trialYielding = False

            # Check if the strain increment causes yielding from an elastic state (initiation of unloading is elastic)
            if abs(deps) > abs(depsToYield) and not unloading:
                self.trialStress = self.trialStress + (self.alpha * self.E - self.E) * (deps - depsToYield)
                self.trialBackStress = self.committedBackStress + self.alpha * self.E * (deps - depsToYield)
                Et = self.alpha * self.E
                self.trialYielding = True

        else:

            # Continue plastic loading
            self.trialStress = self.committedStress + self.alpha * self.E * deps
            self.trialBackStress = self.committedBackStress + self.alpha * self.E * deps
            Et = self.alpha * self.E
            self.trialYielding = True

        return self.trialStress, Et

    # -------------------------------------------------
    # Equilibrium iterations to find strain for given stress
    # -------------------------------------------------
    def equilibrium(self, stressIncrement):

        # Set iteration controls
        maxiter = 100
        tol = 1e-5

        # Initialize incremental trial strain
        strainIncrement = 0.0

        # Iterate to achieve equilibrium
        for j in range(maxiter):

            self.trialStrain = self.committedStrain + strainIncrement

            trialStressIncrement, stiffness = self.stateDetermination(np.array([0.0, strainIncrement, 0.0]))
            stressIncrementResidual = stressIncrement - trialStressIncrement
            residualNorm = np.abs(stressIncrementResidual)

            if residualNorm < tol:
                break

            strainIncrement += stressIncrementResidual/stiffness

        if residualNorm >= tol:
            print("ERROR: Equilibrium iterations in bilinear material failed with residual", residualNorm)
            import sys
            sys.exit()

        return strainIncrement, stiffness

    # -------------------------------------------------
    # Initial stiffness
    # -------------------------------------------------
    def initialStiffness(self):

        return self.E

    # -------------------------------------------------
    # First-order sensitivity analysis
    # -------------------------------------------------
    def stateDerivative(self, eps, ddmParameter, ddmIndex, ddmIsHere, dKflag='none'):

        # Extend arrays in the first increment
        if ddmIndex > len(self.trialStressDerivative) - 1:
            self.trialStrainDerivative.append(0)
            self.trialStressDerivative.append(0)
            self.trialBackStressDerivative.append(0)
            self.committedStrainDerivative.append(0)
            self.committedStressDerivative.append(0)
            self.committedBackStressDerivative.append(0)

        # This is an incremental material
        deps = eps[1]

        # Check for unloading
        unloading = False
        if (self.committedYielding == True and (self.committedStress - self.committedBackStress) * deps < 0):
            unloading = True

        # Check if the last state was elastic or if the strain increment implies unloading from yielding
        if self.committedYielding is False or unloading:

            # Strain increment that would cause yielding
            depsToYield = (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E

            if ddmParameter == 'E' and ddmIsHere:
                depsToYieldDerivative = (self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
            elif ddmParameter == 'fy' and ddmIsHere:
                depsToYieldDerivative = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E
            elif ddmParameter == 'alpha' and ddmIsHere:
                depsToYieldDerivative = (self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E
            else:
                depsToYieldDerivative = (self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E

            # Keep elastic state handy, in case that becomes the conclusion
            if ddmParameter == 'E' and ddmIsHere:
                stressDerivative = self.committedStressDerivative[ddmIndex] + deps - self.E * self.committedStrainDerivative[ddmIndex]
            elif ddmParameter == 'fy' and ddmIsHere:
                stressDerivative = self.committedStressDerivative[ddmIndex] - self.E * self.committedStrainDerivative[ddmIndex]
            elif ddmParameter == 'alpha' and ddmIsHere:
                stressDerivative = self.committedStressDerivative[ddmIndex] - self.E * self.committedStrainDerivative[ddmIndex]
            else:
                stressDerivative = self.committedStressDerivative[ddmIndex] - self.E * self.committedStrainDerivative[ddmIndex]

            if ddmParameter == 'E' and ddmIsHere:
                dEt = 1.0
            elif ddmParameter == 'fy' and ddmIsHere:
                dEt = 0.0
            elif ddmParameter == 'alpha' and ddmIsHere:
                dEt = 0.0
            else:
                dEt = 0.0

            # Check if the strain increment causes yielding from an elastic state (initiation of unloading is elastic)
            if abs(deps) > abs(depsToYield) and not unloading:

                if ddmParameter == 'E' and ddmIsHere:
                    stressDerivative = stressDerivative + (self.alpha - 1) * (deps - depsToYield)  + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                elif ddmParameter == 'fy' and ddmIsHere:
                    stressDerivative = stressDerivative  + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                elif ddmParameter == 'alpha' and ddmIsHere:
                    stressDerivative = stressDerivative + self.E * (deps - depsToYield)  + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                else:
                    stressDerivative = stressDerivative  + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)

                if ddmParameter == 'E' and ddmIsHere:
                    dEt = self.alpha
                elif ddmParameter == 'fy' and ddmIsHere:
                    dEt = 0.0
                elif ddmParameter == 'alpha' and ddmIsHere:
                    dEt = self.E
                else:
                    dEt = 0.0
        else:

            # Continue plastic loading
            if ddmParameter == 'E' and ddmIsHere:
                stressDerivative = self.committedStressDerivative[ddmIndex] + self.alpha * deps + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'fy' and ddmIsHere:
                stressDerivative = self.committedStressDerivative[ddmIndex] + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'alpha' and ddmIsHere:
                stressDerivative = self.committedStressDerivative[ddmIndex] + self.E * deps + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex])
            else:
                stressDerivative = self.committedStressDerivative[ddmIndex] + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex])

            if ddmParameter == 'E' and ddmIsHere:
                dEt = self.alpha
            elif ddmParameter == 'fy' and ddmIsHere:
                dEt = 0.0
            elif ddmParameter == 'alpha' and ddmIsHere:
                dEt = self.E
            else:
                dEt = 0.0

        if dKflag == 'Initial':

            if ddmParameter == 'E' and ddmIsHere:
                dEt = 1.0
            elif ddmParameter == 'fy' and ddmIsHere:
                dEt = 0.0
            elif ddmParameter == 'alpha' and ddmIsHere:
                dEt = 0.0
            else:
                dEt = 0.0

        dKdu = 0

        return stressDerivative, dEt, dKdu

    # -------------------------------------------------
    # Second-order DDM sensitivity analysis
    # -------------------------------------------------
    def stateSecondDerivative(self, eps, secondOrderIndex, ddmParameter1, ddmIndex1, ddmIsHere1, ddmParameter2, ddmIndex2, ddmIsHere2, dKflag='none'):

        # Extend arrays in the first increment
        if secondOrderIndex > len(self.trialSecondStressDerivative) - 1:
            self.trialSecondStrainDerivative.append(0)
            self.trialSecondStressDerivative.append(0)
            self.trialSecondBackStressDerivative.append(0)
            self.committedSecondStrainDerivative.append(0)
            self.committedSecondStressDerivative.append(0)
            self.committedSecondBackStressDerivative.append(0)

        # This is an incremental material
        deps = eps[1]

        # Check for unloading
        unloading = False
        if (self.committedYielding == True and (self.committedStress - self.committedBackStress) * deps < 0):
            unloading = True

        # Check if the last state was elastic or if the strain increment implies unloading from yielding
        if self.committedYielding is False or unloading:

            # Strain increment that would cause yielding
            depsToYield = (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E
            if ddmParameter1 == 'E' and ddmIsHere1:
                depsToYieldDerivative1 = (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2 \
                                                  - (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E**2 \
                                                  + 2 * (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**3
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative =  (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                   - (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E ** 2
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E ** 2
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E ** 2

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                depsToYieldDerivative1 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (np.sign(deps) + self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                depsToYieldDerivative1 = (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E

            else:
                depsToYieldDerivative1 = (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E

            # Keep elastic state handy, in case that becomes the conclusion
            if ddmParameter1 == 'E' and ddmIsHere1:
                stressDerivative = self.committedStressDerivative[ddmIndex1] + deps - self.E * self.committedStrainDerivative[ddmIndex1]
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex2] \
                                             - self.committedStrainDerivative[ddmIndex1] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex2] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex2] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex2] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                stressDerivative = self.committedStressDerivative[ddmIndex1] - self.E * self.committedStrainDerivative[ddmIndex1]
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex1] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                stressDerivative = self.committedStressDerivative[ddmIndex1] - self.E * self.committedStrainDerivative[ddmIndex1]
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex1] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]

            else:
                stressDerivative = self.committedStressDerivative[ddmIndex1] - self.E * self.committedStrainDerivative[ddmIndex1]
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.committedStrainDerivative[ddmIndex1] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             - self.E * self.committedSecondStrainDerivative[secondOrderIndex]

            # Second-order derivative of stiffness
            ddEt = 0.0

            # Check if the strain increment causes yielding from an elastic state (initiation of unloading is elastic)
            if abs(deps) > abs(depsToYield) and not unloading:

                if ddmParameter1 == 'E' and ddmIsHere1:
                    stressDerivative = stressDerivative + (self.alpha - 1) * (deps - depsToYield) + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (deps - depsToYield) \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                elif ddmParameter1 == 'fy' and ddmIsHere1:
                    stressDerivative = stressDerivative  + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                elif ddmParameter1 == 'alpha' and ddmIsHere1:
                    stressDerivative = stressDerivative + self.E * (deps - depsToYield) + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (deps - depsToYield) \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        secondStressDerivative = secondStressDerivative \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                else:
                    stressDerivative = stressDerivative  + (self.alpha * self.E - self.E) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha - 1) * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        secondStressDerivative = secondStressDerivative \
                                                 + self.E * (-self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        secondStressDerivative = secondStressDerivative \
                                                 + (self.alpha * self.E - self.E) * (-self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                # Second-order derivative of stiffness
                ddEt = 0.0
                if ddmParameter1 == 'E' and ddmIsHere1:
                    if ddmParameter2 == 'alpha' and ddmIsHere2:
                        ddEt = 1.0
                elif ddmParameter1 == 'alpha' and ddmIsHere1:
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        ddEt = 1.0

        else:

            # Continue plastic loading
            if ddmParameter1 == 'E' and ddmIsHere1:
                stressDerivative = self.committedStressDerivative[ddmIndex1] + self.alpha * deps + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + deps \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                stressDerivative = self.committedStressDerivative[ddmIndex1] + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                stressDerivative = self.committedStressDerivative[ddmIndex1] + self.E * deps + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + deps \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex2]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
            else:
                stressDerivative = self.committedStressDerivative[ddmIndex1] + self.alpha * self.E * (-self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.E * (-self.committedStrainDerivative[ddmIndex1]) \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    secondStressDerivative = self.committedSecondStressDerivative[secondOrderIndex] \
                                             + self.alpha * self.E * (-self.committedSecondStrainDerivative[secondOrderIndex])

            # Second-order derivative of stiffness
            ddEt = 0.0
            if ddmParameter1 == 'E' and ddmIsHere1:
                if ddmParameter2 == 'alpha' and ddmIsHere2:
                    ddEt = 1.0
            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                if ddmParameter2 == 'E' and ddmIsHere2:
                    ddEt = 1.0

        if dKflag == 'Initial':
            ddEt = 0.0

        return secondStressDerivative, ddEt

    # -------------------------------------------------
    # Commit
    # -------------------------------------------------
    def commit(self):

        self.committedStrain = self.trialStrain
        self.committedStress = self.trialStress
        self.committedBackStress = self.trialBackStress
        self.committedYielding = self.trialYielding

        self.committedStrainDerivative[:] = self.trialStrainDerivative[:]
        self.committedStressDerivative[:] = self.trialStressDerivative[:]
        self.committedBackStressDerivative[:] = self.trialBackStressDerivative[:]

        self.committedSecondStrainDerivative[:] = self.trialSecondStrainDerivative[:]
        self.committedSecondStressDerivative[:] = self.trialSecondStressDerivative[:]
        self.committedSecondBackStressDerivative[:] = self.trialSecondBackStressDerivative[:]

    # -------------------------------------------------
    # Commit sensitivity history variables
    # -------------------------------------------------
    def commitSensitivity(self, eps, strainSensitivity, ddmParameter, ddmIndex, ddmIsHere):

        # This is an incremental material
        deps = eps[1]
        self.trialStrainDerivative[ddmIndex] = strainSensitivity

        # Check for unloading
        unloading = False
        if (self.committedYielding == True and (self.committedStress - self.committedBackStress) * deps < 0):
            unloading = True

        # Check if the last state was elastic or if the strain increment implies unloading from yielding
        if self.committedYielding is False or unloading:

            # Strain increment that would cause yielding
            depsToYield = (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E

            if ddmParameter == 'E' and ddmIsHere:
                depsToYieldDerivative = (self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E ** 2
            elif ddmParameter == 'fy' and ddmIsHere:
                depsToYieldDerivative = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E
            elif ddmParameter == 'alpha' and ddmIsHere:
                depsToYieldDerivative = (self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E
            else:
                depsToYieldDerivative = (self.committedBackStressDerivative[ddmIndex] - self.committedStressDerivative[ddmIndex]) / self.E

            # Keep elastic state handy, in case that becomes the conclusion
            if ddmParameter == 'E' and ddmIsHere:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + deps + self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'fy' and ddmIsHere:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'alpha' and ddmIsHere:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            else:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])

            # Check if the strain increment causes yielding from an elastic state (initiation of unloading is elastic)
            if abs(deps) > abs(depsToYield) and not unloading:

                if ddmParameter == 'E' and ddmIsHere:
                    self.trialStressDerivative[ddmIndex] = self.trialStressDerivative[ddmIndex] + (self.alpha - 1) * (deps - depsToYield) + (self.alpha * self.E - self.E) * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                elif ddmParameter == 'fy' and ddmIsHere:
                    self.trialStressDerivative[ddmIndex] = self.trialStressDerivative[ddmIndex] + (self.alpha * self.E - self.E) * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                elif ddmParameter == 'alpha' and ddmIsHere:
                    self.trialStressDerivative[ddmIndex] = self.trialStressDerivative[ddmIndex] + self.E * (deps - depsToYield) + (self.alpha * self.E - self.E) * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                else:
                    self.trialStressDerivative[ddmIndex] = self.trialStressDerivative[ddmIndex]  + (self.alpha * self.E - self.E) * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)

                if ddmParameter == 'E' and ddmIsHere:
                    self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.alpha * (deps - depsToYield) + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                elif ddmParameter == 'fy' and ddmIsHere:
                    self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                elif ddmParameter == 'alpha' and ddmIsHere:
                    self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.E * (deps - depsToYield) + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)
                else:
                    self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex] - depsToYieldDerivative)

        else:

            # Continue plastic loading
            if ddmParameter == 'E' and ddmIsHere:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.alpha * deps + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'fy' and ddmIsHere:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'alpha' and ddmIsHere:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.E * deps + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            else:
                self.trialStressDerivative[ddmIndex] = self.committedStressDerivative[ddmIndex] + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])

            if ddmParameter == 'E' and ddmIsHere:
                self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.alpha * deps + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'fy' and ddmIsHere:
                self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            elif ddmParameter == 'alpha' and ddmIsHere:
                self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.E * deps + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])
            else:
                self.trialBackStressDerivative[ddmIndex] = self.committedBackStressDerivative[ddmIndex] + self.alpha * self.E * (strainSensitivity - self.committedStrainDerivative[ddmIndex])

        return

    # -------------------------------------------------
    # Commit second-order sensitivity history variables
    # -------------------------------------------------
    def commitSecondSensitivity(self, eps, depsilon1, depsilon2, ddepsilon, secondOrderIndex, ddmParameter1, ddmIndex1, ddmIsHere1, ddmParameter2, ddmIndex2, ddmIsHere2):

        # This is an incremental material
        deps = eps[1]
        self.trialSecondStrainDerivative[secondOrderIndex] = ddepsilon

        # Check for unloading
        unloading = False
        if (self.committedYielding == True and (self.committedStress - self.committedBackStress) * deps < 0):
            unloading = True

        # Check if the last state was elastic or if the strain increment implies unloading from yielding
        if self.committedYielding is False or unloading:

            # Strain increment that would cause yielding
            depsToYield = (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E
            if ddmParameter1 == 'E' and ddmIsHere1:
                depsToYieldDerivative1 = (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2 \
                                                  - (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E**2 \
                                                  + 2 * (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**3
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative =  (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                   - (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E ** 2
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E ** 2
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E ** 2

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                depsToYieldDerivative1 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (np.sign(deps) + self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                depsToYieldDerivative1 = (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E

            else:
                depsToYieldDerivative1 = (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E
                if ddmParameter2 == 'E' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E - (np.sign(deps) * self.fy + self.committedBackStress - self.committedStress) / self.E**2
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E \
                                                  - (self.committedBackStressDerivative[ddmIndex1] - self.committedStressDerivative[ddmIndex1]) / self.E**2
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    depsToYieldDerivative2 = (np.sign(deps) + self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E
                else:
                    depsToYieldDerivative2 = (self.committedBackStressDerivative[ddmIndex2] - self.committedStressDerivative[ddmIndex2]) / self.E
                    secondDepsToYieldDerivative = (self.committedSecondBackStressDerivative[secondOrderIndex] - self.committedSecondStressDerivative[secondOrderIndex]) / self.E

            # Keep elastic state handy, in case that becomes the conclusion
            if ddmParameter1 == 'E' and ddmIsHere1:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + deps + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                                                                         + (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            else:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                                                                         + self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            # Check if the strain increment causes yielding from an elastic state (initiation of unloading is elastic)
            if abs(deps) > abs(depsToYield) and not unloading:

                if ddmParameter1 == 'E' and ddmIsHere1:
                    # self.trialStressDerivative[ddmIndex1] = self.trialStressDerivative[ddmIndex1] + (self.alpha - 1) * (deps - depsToYield) + (self.alpha * self.E - self.E) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha - 1) * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + (self.alpha - 1) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha - 1) * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (deps - depsToYield) \
                                                                             + (self.alpha - 1) * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha - 1) * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                elif ddmParameter1 == 'fy' and ddmIsHere1:
                    # self.trialStressDerivative[ddmIndex1] = self.trialStressDerivative[ddmIndex1] + (self.alpha * self.E - self.E) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha - 1) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                elif ddmParameter1 == 'alpha' and ddmIsHere1:
                    # self.trialStressDerivative[ddmIndex1] = self.trialStressDerivative[ddmIndex1] + self.E * (deps - depsToYield) + (self.alpha * self.E - self.E) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (deps - depsToYield) \
                                                                             + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + (self.alpha - 1) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                else:
                    # self.trialStressDerivative[ddmIndex1] = self.trialStressDerivative[ddmIndex1] + (self.alpha * self.E - self.E) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha - 1) * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondStressDerivative[secondOrderIndex] = self.trialSecondStressDerivative[secondOrderIndex] \
                                                                             + (self.alpha * self.E - self.E) * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                if ddmParameter1 == 'E' and ddmIsHere1:
                    # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.alpha * (deps - depsToYield) + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative1) \
                        + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative2) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + (deps - depsToYield) \
                        + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                elif ddmParameter1 == 'fy' and ddmIsHere1:
                    # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                elif ddmParameter1 == 'alpha' and ddmIsHere1:
                    # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.E * (deps - depsToYield) + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + (deps - depsToYield) \
                        + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2] - depsToYieldDerivative2) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)

                else:
                    # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1)
                    if ddmParameter2 == 'E' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'fy' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    elif ddmParameter2 == 'alpha' and ddmIsHere2:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1] - depsToYieldDerivative1) \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
                    else:
                        self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                        + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex] - secondDepsToYieldDerivative)
        else:

            # Continue plastic loading
            if ddmParameter1 == 'E' and ddmIsHere1:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.alpha * deps + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + deps \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.E * deps + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + deps \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            else:
                # self.trialStressDerivative[ddmIndex1] = self.committedStressDerivative[ddmIndex1] + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondStressDerivative[secondOrderIndex] = self.committedSecondStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            if ddmParameter1 == 'E' and ddmIsHere1:
                # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.alpha * deps + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + deps \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'fy' and ddmIsHere1:
                # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            elif ddmParameter1 == 'alpha' and ddmIsHere1:
                # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.E * deps + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + deps \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon2 - self.committedStrainDerivative[ddmIndex2]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

            else:
                # self.trialBackStressDerivative[ddmIndex1] = self.committedBackStressDerivative[ddmIndex1] + self.alpha * self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1])
                if ddmParameter2 == 'E' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'fy' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                elif ddmParameter2 == 'alpha' and ddmIsHere2:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.E * (depsilon1 - self.committedStrainDerivative[ddmIndex1]) \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])
                else:
                    self.trialSecondBackStressDerivative[secondOrderIndex] = self.committedSecondBackStressDerivative[secondOrderIndex] \
                    + self.alpha * self.E * (ddepsilon - self.committedSecondStrainDerivative[secondOrderIndex])

        return

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

        return self.committedStress