二、程序實現及說明
（一）、總體說明：
1、本程序是平面應力或平面應變問題的矩形單元的靜力分析程序，本主要用以較規則的平面結構的有限元分析。
2、本程序要求結構整體坐標水平向右為X軸正向，豎直向上為Y軸正向，矩形單元局部坐標和整體坐標一致。由于結構和矩形單元的規則性，這是很容易做到的。而且要求在編制單元結點號I，J，M，P的規則是從左到右、從下到上，逆時針。這一點也是很容易做到的。
(二)、標識符說明：
NO---算例標號。整型變量，輸入參數。
TL(30)—算例標題。字符型數組，輸入參數。
NJ---結點總數。整型變量，輸入參數。
NE---單元總數。整型變量，輸入參數。
NZ---約束結點數。整型變量，輸入參數。
NPJ---結點荷載數。整型變量，輸入參數。
IPS---問題類型碼。整型變量，輸入參數。
E---彈性模量。實型變量，輸入參數。
PR---泊松比。實型變量，輸入參數。
T---彈性體厚度。實型變量，輸入參數。
V---容重。實型變量，輸入參數。
LND(NE,4)---單元結點編號數組。LND(NE,1)、LND(NE,2)、LND(NE,3)、LND(NE,4)分別為單元結點的單元NE的四個結點編號。整型數組，輸入參數。
AA(NE)---單元a值數組。AA(I)為單元號為I的a值。實型數組，輸入參數。
BB(NE)---單元b值數組。BB(I)為單元號為I的b值。實型數組，輸入參數。
JZ(NZ,3)---支承結點數組。JZ(I,1)為第I個支承結點的結點編號。JZ(I,2)、JZ(I,2)分別是第I個支承結點在位移方向的約束信息。該信息為1時，表示有約束；為0時則無約束。整型數組，輸入參數。
PJ(NPJ,3)---結點荷載數組。PJ(I,1)表示第I個結點荷載作用結點的結點編號。PJ(I,2)、PJ(I,3)分別表示該結點沿方向作用的結點荷載數值。實型數組，輸入參數。
N---自由度。整型變量，N=NJ×2。
ND---半帶寬。整型變量。
NPJ0---存放結點荷載數NPJ信息，以保證當NPJ0=0時，亦可以按規定定義數組。
AK(N,ND)---整體剛度矩陣。
AKE(8,8)---單元剛度矩陣。
B(3,8)--- 位移-應力轉換矩陣(幾何矩陣)。
D(3,3)---彈性矩陣。
S(3,8)---應力矩陣，S=DB。
P(N)---結點荷載數組，存放結點荷載列向量；解方程后存結點位移。
DE(8)---單元的結點位移。
ST(NE,4)---單元結點應力。ST(I,1)、ST(I,2) 、ST(I,3) 、ST(I,4)分別表示單元I的第I，J，M，P號結點應力。
READ---子程序名，輸入數據。
MKE---子程序名，形成單元剛度矩陣。
MD---子程序名，計算彈性矩陣。
MB---子程序名，計算位移-應力轉換矩陣。
MF---子程序名，求結點荷載。
RKR---子程序名，引入支承條件，修改剛度方程。
SLOV---子程序名，求結點位移并輸出。
MADE---子程序名，計算應力并輸出。
(三)、數據輸入格式：
算例 
標題 
基本數據 結點總數 單元總數 約束結點數 結點荷載數 問題類型

材料性質 彈性模量 泊松比 單元厚度 容重

單元結點編號 按i，j，m，p結點順序依次填寫

單元a、b值 按單元號及a、b順序填寫

支承結點信息 按約束的結點號、X向信息、Y向信息依次填寫

結點荷載信息 按荷載作用結點號、X向荷載值、Y向荷載值依次填寫 

注意：該程序與三角形單元程序的對比
1、從所需輸入數據看。區別只在于三角形單元需要各結點的坐標值（X(NJ),Y(NJ)），而矩形單元需要的是單元的a、b值（AA(NE),BB(NE)）。
2、從程序的過程看。
（1）、一些數組的階數有所變化，這是由于三角形是三結點，而四邊形是四結點。區別在如下數組：單元結點編號數組、單元剛度矩陣、幾何矩陣、應力矩陣S、單元的結點位移、單元(結點)應力，三角形單元分別為LND(NE,3)、AKE(6,6) 、B(3,6) 、S(3,6) 、DE(6) 、ST(3)；矩形單元分別為LND(NE,3)、AKE(8,8)、 B(3,8)、 S(3,8)、 DE(8) 、ST(NE,4)。
（2）、由于矩形單元計算單元面積很方便，所以不需要專門的子程序來計算其單元面積；而三角形單元就不一樣，它需要用一個子程序來計算其單元面積（MA子程序）。
（3）、三角形單元是線性單元，從而位移-應力轉換矩陣(幾何矩陣)B是常量，又應力矩陣S=DB，因D是由材料性質決定，對同一材料而言D是常量，所以S也是常量（均對同一單元各點而言），所以同一單元上的各點應力相同；矩形單元是雙線性單元，位移-應力轉換矩陣(幾何矩陣)B對不同點而言是變量，從而應力矩陣S=DB對不同點也是變量，所以同一單元上各點的應力不同。這就使二者在形成B矩陣、S矩陣的程序方法上不一樣，從而導致二者生成單元剛度矩陣K 的程序方式也不一樣。