## Appendix B Listing from Author'S Program
nu = 0.2
#Intensityoftheload q_0=10_000.00#N
#Flexuralrigidityofaplet D=(E*h*h*3)/(12*(1-nu*2))
#---------------------------------------------------------------- #Helperfunctions
#----------------------------------------------------------------  defa(s):  ifs==1:  returna_1 elifs==2:  returna_2 defx(s):  defwrapper(x_1,x_2):  ifs==1:  returnx_1 elifs==2:  returnx_2  defgamma(k,s):  return(k*np.pi)/a(s) defdelta(k,s):  return((2*k-1)*np.pi)/(2*a(s)) defkappa(k,p,s):  ifp==1:  returngamma(k,s) elifp==2:  returndelta(k,s) defT(k,p,s):  defwrapper(x_1,x_2):  returnnp.cos(kappa(k,p,s)*x(s)(x_1,x_2))  returnwrapper defcalculate(fs,solution,x_1,x_2):  coeffs=np.array([f(x_1,x_2)forfinfs])  returnnp.einsum('ijk,i->jk',coeffs,solution) defparticular_solution(C,f):  defwrapper(x_1,x_2):  returnC*f(x_1,x_2) returnwrapper deffinal(general_solution,particular_solution):  result=[]  forfingeneral_solution:  
result.append(f) result.append(particular_solution) returnresult
#----------------------------------------------------------------
#Relationsknownfromtheoryofthinisotropicplates
#----------------------------------------------------------------  
defphi_1(w): defwrapper(x_1,x_2): returngrad(w,0)(x_1,x_2) returnwrapper defphi_2(w): defwrapper(x_1,x_2): returngrad(w,1)(x_1,x_2) returnwrapper defM_11(w): defwrapper(x_1,x_2): return-D*(grad(grad(w,0),0)(x_1,x_2)+ nu*grad(grad(w,1),1)(x_1,x_2)) returnwrapper defM_22(w): defwrapper(x_1,x_2): return-D*(grad(grad(w,1),1)(x_1,x_2)+ nu*grad(grad(w,0),0)(x_1,x_2)) returnwrapper defM_12(w): defwrapper(x_1,x_2): return-D*(1-nu)*grad(grad(w,0),1)(x_1,x_2) returnwrapper defQ_1(w): defwrapper(x_1,x_2): return-D*(grad(grad(grad(w,0),0),0)(x_1,x_2)+ grad(grad(grad(w,0),1),1)(x_1,x_2)) returnwrapper defQ_2(w): defwrapper(x_1,x_2): return-D*(grad(grad(grad(w,0),0),1)(x_1,x_2)+ grad(grad(grad(w,1),1),1)(x_1,x_2)) returnwrapper defV_1(w): defwrapper(x_1,x_2): return-D*(grad(grad(grad(w,0),0),0)(x_1,x_2)+ (2-nu)*grad(grad(grad(grad(w,0),1),1)(x_1,x_2)) returnwrapper defV_2(w): defwrapper(x_1,x_2):return-D*(grad(grad(grad((w,1),1),1)(x_1,x_2)+ +  (2-nu)*grad(grad(grad((w,0),0),1)(x_1,x_2)) return wrapper
#--------------------------------------------------------------------------------------------------------------------------------
#Forcefunctions
#----------------------------------------------------------------  
defforce_functions(m,n,p,q): defwrapper(x_1,x_2): returnT(m,p,1)(x_1,x_2)*T(n,q,2)(x_1,x_2) returnwrapper
#----------------------------------------------------------------
#Shapefunctions
#----------------------------------------------------------------
#Basefunctionsofthesolution defB(k,p,s,ni): defwrapper(x_1,x_2): ifni=1: return(np.cosh(kappa(k,p,3-s)*x(s)(x_1,x_2))) elifni--2: return((x(s)(x_1,x_2)/a(s))* (np.sinh(kappa(k,p,3-s)*x(s)(x_1,x_2)))) returnwrapper
#Shapefunctions defW(k,p,s,ni): defwrapper(x_1,x_2): returnB(k,p,s,ni)(x_1,x_2)*T(k,p,3-s)(x_1,x_2) returnwrapper  
defshape_functions(K): result=[] forkinrange(1,K+1): forpinrange(1,3): forsinrange(1,3): forniinrange(1,3): result.append(W(k,p,s,ni)) returnresult
#----------------------------------------------------------------
#Buildmodel
#----------------------------------------------------------------  
W_g=shape_functions(_K) W_p=force_functions(1,1,2,2)
#----------------------------------------------------------------
#----------------------------------------------------------------
#Calculatederivativesofshapefunction U_g=[phi_1(f)forfinW_g] V_g=[phi_2(f)forfinW_g] X_g=[M_11(f)forfinW_g] Y_g=[M_22(f)forfinW_g] Z_g=[M_12(f)forfinW_g] B_g=[Q_1(f)forfinW_g]  
n = 2*_X
Initialize matrix blocks B_1 = np.zeros((m, n)) B_2 = np.zeros((m, n)) B_3 = np.zeros((m, n)) B_4 = np.zeros((m, n))
Fill matrix blocks with values for boundary conditions for i in range(m): B_1[i] = W[i](b_11, b_12) B_2[i] = U[i](b_11, b_12) B_3[i] = W[i](b_21, b_22) B_4[i] = V[i](b_21, b_22)