Utility function to help implement boundary conditions by specifying 
    them for and applying them to the elastic thickness grid

def BC_Rigidity(self):
    """
    Utility function to help implement boundary conditions by specifying 
    them for and applying them to the elastic thickness grid
    """

    #########################################
    # FLEXURAL RIGIDITY BOUNDARY CONDITIONS #
    #########################################
    # West
    if self.BC_W == 'Periodic':
      self.BC_Rigidity_W = 'periodic'
    elif (self.BC_W == np.array(['0Displacement0Slope', '0Moment0Shear', '0Slope0Shear'])).any():
      self.BC_Rigidity_W = '0 curvature'
    elif self.BC_W == 'Mirror':
      self.BC_Rigidity_W = 'mirror symmetry'
    else:
      sys.exit("Invalid Te B.C. case")
    # East
    if self.BC_E == 'Periodic':
      self.BC_Rigidity_E = 'periodic'
    elif (self.BC_E == np.array(['0Displacement0Slope', '0Moment0Shear', '0Slope0Shear'])).any():
      self.BC_Rigidity_E = '0 curvature'
    elif self.BC_E == 'Mirror':
      self.BC_Rigidity_E = 'mirror symmetry'
    else:
      sys.exit("Invalid Te B.C. case")
    
    #############
    # PAD ARRAY #
    #############
    if np.isscalar(self.Te):
      self.D *= np.ones(self.qs.shape) # And leave Te as a scalar for checks
    else:
      self.Te_unpadded = self.Te.copy()
    # F2D keeps this inside the "else" and handles this differently, 
    # largely because it has different ways of computing the flexural
    # response with variable Te. We'll keep everything simpler here and 
    # just pad this array so it can be sent through the same process
    # to create the coefficient arrays.
    self.D = np.hstack([np.nan, self.D, np.nan])

    ###############################################################
    # APPLY FLEXURAL RIGIDITY BOUNDARY CONDITIONS TO PADDED ARRAY #
    ###############################################################
    if self.BC_Rigidity_W == "0 curvature":
      self.D[0] = 2*self.D[1] - self.D[2]
    if self.BC_Rigidity_E == "0 curvature":
      self.D[-1] = 2*self.D[-2] - self.D[-3]
    if self.BC_Rigidity_W == "mirror symmetry":
      self.D[0] = self.D[2]
    if self.BC_Rigidity_E == "mirror symmetry":
      self.D[-1] = self.D[-3]
    if self.BC_Rigidity_W == "periodic":
      self.D[0] = self.D[-2]
    if self.BC_Rigidity_E == "periodic":
      self.D[-1] = self.D[-3]