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")
    # North
    if self.BC_N == 'Periodic':
      self.BC_Rigidity_N = 'periodic'
    elif (self.BC_N == np.array(['0Displacement0Slope', '0Moment0Shear', '0Slope0Shear'])).any():
      self.BC_Rigidity_N = '0 curvature'
    elif self.BC_N == 'Mirror':
      self.BC_Rigidity_N = 'mirror symmetry'
    else:
      sys.exit("Invalid Te B.C. case")
    # South
    if self.BC_S == 'Periodic':
      self.BC_Rigidity_S = 'periodic'
    elif (self.BC_S == np.array(['0Displacement0Slope', '0Moment0Shear', '0Slope0Shear'])).any():
      self.BC_Rigidity_S = '0 curvature'
    elif self.BC_S == 'Mirror':
      self.BC_Rigidity_S = '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()
      self.Te = np.hstack(( np.nan*np.zeros((self.Te.shape[0], 1)), self.Te, np.nan*np.zeros((self.Te.shape[0], 1)) ))
      self.Te = np.vstack(( np.nan*np.zeros(self.Te.shape[1]), self.Te, np.nan*np.zeros(self.Te.shape[1]) ))
      self.D = np.hstack(( np.nan*np.zeros((self.D.shape[0], 1)), self.D, np.nan*np.zeros((self.D.shape[0], 1)) ))
      self.D = np.vstack(( np.nan*np.zeros(self.D.shape[1]), self.D, np.nan*np.zeros(self.D.shape[1]) ))

    ###############################################################
    # 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_N == "0 curvature":
      self.D[0,:] = 2*self.D[1,:] - self.D[2,:]
    if self.BC_Rigidity_S == "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_N == "mirror symmetry":
      self.D[0,:] = self.D[2,:] # Yes, will work on corners -- double-reflection
    if self.BC_Rigidity_S == "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]
    if self.BC_Rigidity_N == "periodic":
      self.D[0,:] = self.D[-2,:]
    if self.BC_Rigidity_S == "periodic":
      self.D[-1,:] = self.D[-3,:]