八節(jié)點(diǎn)等參元

1、?計(jì)算力學(xué)?課程大作業(yè)八節(jié)點(diǎn)等參單元平面有限元分析程序土木工程學(xué)院目錄?計(jì)算力學(xué)?課程大作業(yè)1 .概述通常情況下的有限元分析過程是運(yùn)用可視化分析軟件如ANSYS、SAP等進(jìn)行前處理和后處理,而中間的計(jì)算局部一般采用自己編制的程序來運(yùn)算.具有較強(qiáng)數(shù)值計(jì)算和處理能力的Fortran語言是傳統(tǒng)有限元計(jì)算的首選語言.隨著有限元技術(shù)的逐步成熟,它被應(yīng)用在越來越復(fù)雜的問題處理中,但在實(shí)際應(yīng)用中也暴露出一些問題.有時(shí)網(wǎng)格離散化的區(qū)域較大,而又限于研究精度的要求,使得劃分的網(wǎng)格數(shù)目極其龐大,結(jié)點(diǎn)數(shù)可多達(dá)數(shù)萬個(gè), 從而造成計(jì)算中要運(yùn)算的數(shù)據(jù)量巨大,程序運(yùn)行的時(shí)間較長的弊端,這就延長了問題解決的時(shí)間,使得求解效率

2、降低.由于運(yùn)行周期長,不利于程序的調(diào)試,特別是對于要計(jì)算多種運(yùn)行工況時(shí) 的情況；同時(shí)大數(shù)據(jù)量處理對計(jì)算機(jī)的內(nèi)存和CPU提出了更高的要求,而在實(shí)際應(yīng)用中,單靠計(jì)算機(jī)硬件水平的提升來解決問題的水平是有限的.因此,必須尋找新的編程語言.隨著有限元前后處理的不斷開展和完善 ,以及大型工程分析軟件對有限元接口的要求 , 有限元分析程序不應(yīng)只滿足解題功能,它還應(yīng)滿足軟件工程所要求的結(jié)構(gòu)化程序設(shè)計(jì)條件 , 能夠?qū)Υ鎯?chǔ)進(jìn)行動(dòng)態(tài)分配,以充分利用計(jì)算機(jī)資源,它還應(yīng)很容易地與其它軟件如 CAD的 實(shí)體造型,優(yōu)化設(shè)計(jì)等接口.現(xiàn)在可編寫工程應(yīng)用軟件的計(jì)算機(jī)語言較多 ,其中C語言是一 個(gè)較為優(yōu)秀的語言,很容易滿足現(xiàn)在有限

3、元分析程序編程的要求.C語言最初是為操作系統(tǒng)、 編譯器以及文字處理等編程而創(chuàng)造的.隨著不斷完善,它已應(yīng)用到其它領(lǐng)域,包括工程應(yīng)用軟件的編程.近年來,C語言已經(jīng)成為計(jì)算機(jī)領(lǐng)域最普及的一個(gè)編程語言,幾乎世界上所有的計(jì)算機(jī)都裝有C的編譯器,從PC機(jī)到巨型機(jī)到超巨型的并行機(jī),C與所有的硬件和操作系統(tǒng)聯(lián)系在一起.用C編寫的程序,可移植性極好,幾乎不用作多少修改,就可在任何一臺(tái)裝有 ANSI、C編譯器的計(jì)算機(jī)上運(yùn)行. C既是高級語言, 也是低級語言,也就是說,可用它作數(shù)值計(jì)算,也可用它對計(jì)算機(jī)存儲(chǔ)進(jìn)行操作.2 .編程思想本程序采用C語言編程,編制平面四邊形八節(jié)點(diǎn)等參元程序,用以求解平面結(jié)構(gòu)問題.程序采用二

4、維等帶寬存儲(chǔ)整體剛度矩陣,乘大數(shù)法引入約束,等帶寬高斯消去法求解位移.在有限元程序中,變量數(shù)據(jù)需賦值的可分為節(jié)點(diǎn)信息,單元信息,載荷信息等.對于一個(gè)節(jié)點(diǎn)來說,需以下信息：節(jié)點(diǎn)編號整型,節(jié)點(diǎn)坐標(biāo)實(shí)型,節(jié)點(diǎn)位移實(shí)型, 節(jié)點(diǎn)載荷實(shí)型,邊界條件實(shí)型和節(jié)點(diǎn)溫度實(shí)型等.同樣,對于一個(gè)單元來說, 需以下信息：單元的節(jié)點(diǎn)聯(lián)接信息整型,單元類型信息桁架、梁、板、殼等整型, 單元特性信息厚度、內(nèi)力矩等實(shí)型,材料信息彈性模量,泊松比等實(shí)型等.在FORTRAN程序中,以上這些變量混合在一起,很難識(shí)別,使程序的可讀性不好,如需要進(jìn)行單元網(wǎng)絡(luò)的自適應(yīng)劃分,節(jié)點(diǎn)及單元的修改將非常困難.在進(jìn)行C語言編譯過程中,采用結(jié)構(gòu)str

5、uct使每個(gè)節(jié)點(diǎn)信息存儲(chǔ)在一個(gè)結(jié)構(gòu)體數(shù)組中, 提升程序的可讀性, 使數(shù) 據(jù)結(jié)構(gòu)更趨于合理.2.1. 八節(jié)點(diǎn)矩形單元介紹八節(jié)點(diǎn)矩形單元編號如錯(cuò)誤！未找到引用源.所示八節(jié)點(diǎn)矩形單元的位移函數(shù)為：2u12X3y4X2v910X11 y12X其形函數(shù)為2225xy6y7X y 8xy22213xy14 y15X y 16xyUN1U1N2U2N3U3N4U4N5U5山見N7U7NgUgv NmN2v2N3v3N4v4N5v5N6v6N7v7N8v8式(2.3)和式(2.4)中 Ni Ni(,并且采用等參變換,那么單元的坐標(biāo)變換式可取為XN1X1N2X2N3X3N4X4N5X5N7X7Ng%YMyN2y

6、2N3y3N4y4N5y5N6y6N7y7N8y8(2.(1)(2.(2)(2.(3)(2.(4)(2.(5)單元應(yīng)變矩陣為(2.6)xy式(2.6) 一般簡寫為其中B的子塊矩陣為Nix由于Ni是、的函數(shù),在NiNi從而有因此,單元應(yīng)力矩陣為單元?jiǎng)偠染仃嚍?2.7)BiBi中的NiNi yNi(2.8)y要根據(jù)復(fù)合函數(shù)來求導(dǎo),即NiNixNixNi(2.9)NixNiyNiNi(2.10)(2.11)D B hdxdy(2.12)其中積分采用三點(diǎn)高斯積分,1111f ( , )di jf( i, j):.nipWf( , )i(2.13)nipn2 (高斯積分點(diǎn)的總數(shù))j或Wi是加權(quán)系數(shù),j是

7、單元內(nèi)的坐標(biāo).對2 0.0, 28.0/9.0 ,于三點(diǎn)高斯積分,高斯積分點(diǎn)的位置:3706, 3 5.09.0.單元等效節(jié)點(diǎn)荷載結(jié)構(gòu)剛度矩陣結(jié)構(gòu)結(jié)點(diǎn)荷載列陣TP ds(2.14)KePe(2.15)(2.16)注意,對于式(2.15)和式(2.16)中的理解不是簡單的疊加而是集成o總剛平衡方程從式(2.17)求出(2.17)(2.18)將 回代入式(2.7)和式(2.11),得到2.2. 有限元分析的模塊組織一個(gè)典型的有限元分析過程主要包括以下幾個(gè)步驟：1讀輸入數(shù)據(jù),定義節(jié)點(diǎn)及單元數(shù)組.2由邊界條件計(jì)算方程個(gè)數(shù),賦值荷載列陣.3讀入在帶狀存儲(chǔ)的總剛度矩陣中單元和載荷信息.4定義總剛度陣數(shù)組.

8、5組裝總剛度陣.6解方程組.其流程圖可見 錯(cuò)誤!未找到引用源.3.1.主要變量說明:主要程序變量說明wide分析模型的寬hight分析模型的高wsize寬度方向網(wǎng)格劃分尺寸hsize高度方向網(wǎng)格劃分尺寸npoin節(jié)點(diǎn)總數(shù)nelem單元總數(shù)struct node500節(jié)點(diǎn)結(jié)構(gòu)體變量struct element 500單元結(jié)構(gòu)體變量ifpre500節(jié)點(diǎn)約束信息posgp3高斯積分點(diǎn)weigp3權(quán)系數(shù)gpcod29高斯積分點(diǎn)總體坐標(biāo)bmatx316單元變形矩陣dmatx33單元彈性矩陣xjacm22雅可比矩陣xjaci22雅可比矩陣的逆矩陣djacb雅可比矩陣行列式的值estif136單元?jiǎng)偠染仃噑

9、hape8單元形函數(shù)deriv28形函數(shù)對局部坐標(biāo)的導(dǎo)數(shù)cartd28形函數(shù)對整體坐標(biāo)的導(dǎo)數(shù)A30000總體剛度矩陣eload16等效節(jié)點(diǎn)荷載gpcod29高斯積分點(diǎn)的總體坐標(biāo)主要函數(shù)說明void meshing ()網(wǎng)格劃分void coordinate()節(jié)點(diǎn)整體坐標(biāo)計(jì)算void input()信息輸入void stif()形成單元?jiǎng)偠染仃噕oid sfr2 ()計(jì)算形函數(shù)對當(dāng)前積分點(diǎn)及形函數(shù)對局部坐標(biāo)的到數(shù)值void jacobi2()形成雅可比矩陣void modps()形成彈性矩陣void bmatps()形成應(yīng)變矩陣void load()計(jì)算等效節(jié)點(diǎn)荷載void asem()形成整

10、體剛度矩陣和整體荷載列陣void solve()求解整體方程void stre()計(jì)算單元應(yīng)力任何一個(gè)有限元程序處理過程,都可以劃分為三個(gè)階段:1前處理：讀入數(shù)據(jù)和生成數(shù)據(jù).2處理器：有限元的矩陣計(jì)算、組集和求解.3后處理：輸出節(jié)點(diǎn)位移、單元應(yīng)變等.為了更清楚地描述程序,我們給出了程序的流程圖.5.程序應(yīng)用與ANSYS分析的比擬為了驗(yàn)證程序的正確性, 我們?nèi)×艘粋€(gè)簡單框架結(jié)構(gòu),分別用自編程序和 ANSYS進(jìn)行 分析和比擬.5.1. 問題說明錯(cuò)誤！未找到引用源. 所示一簡支梁,高 3m,長18m,承受均布荷載 10N/m2, 10 E 2 10 Pa ,0.167 ,取t 1m,作為平面應(yīng)力問題

11、.18 m14111111用1箱用11“口儲(chǔ)“11 °N/m：圖5-1由于結(jié)構(gòu)和荷載對稱,只對右邊一半進(jìn)行有限單元計(jì)算,如錯(cuò)誤！未找到引用源.所示,而在y軸各節(jié)點(diǎn)布置水平連桿支座.5.2. ANSYS分析結(jié)果2ANSYS分析中,米用plane82單兀,在板單兀上邊施加均布何載 10N /m ,在y軸上 的各結(jié)點(diǎn)約束UX方向的自由度,在板單元右下角的結(jié)點(diǎn)約束 UX自由度,結(jié)點(diǎn)布置、網(wǎng)格 劃分方案如 錯(cuò)誤！未找到引用源. 所示.ANSYSFEB IS 202112:51!： 15圖5-3應(yīng)力云圖.錯(cuò)誤!未找到引用源.和錯(cuò)誤！未找到引用源.分別給出了結(jié)構(gòu)的位移圖和ANSYS圖5-4位移圖圖

12、5-5各單元 x圖X從錯(cuò)誤!未找到引用源. 和錯(cuò)誤!未找到引用源. 可以看到,跨中的位移和正應(yīng)力最大, 與簡支梁承受均布荷載,跨中撓度和正應(yīng)力最大的力學(xué)常識(shí)相符合,可以初步認(rèn)為模型是正確的.錯(cuò)誤!未找到引用源.給出了 x 0.75m的截面上的正應(yīng)力X和錯(cuò)誤！未找到引用源.給出了 y 1.5m處各橫截面Y方向位移,其中X的單位為Pa, y的單位為mo表格1考查點(diǎn)的y/m0正應(yīng)力 X表格2考查點(diǎn)的x/m0Y方向位移(10-6)05.3. 自編程序分析結(jié)果用自編程序進(jìn)行分析時(shí),采用了整個(gè)結(jié)構(gòu)模型進(jìn)行計(jì)算,其他條件不變,因編程水平有限,在后處理階段沒有給出節(jié)點(diǎn)處的位移與應(yīng)力,而只能得到高斯積分點(diǎn)處的位

13、移與應(yīng)力,得到積分點(diǎn)處的應(yīng)力后, 根據(jù)數(shù)值分析相關(guān)知識(shí)采用了插值進(jìn)行處理,從而得到x 0.75m的截面上的正應(yīng)力X和y1.5m處各橫截面Y方向位移,分別列出表格如下：表格3考查點(diǎn)的y/m0正應(yīng)力 X93表格4考查點(diǎn)的x/m0Y方向位移(10-6)05.4. 結(jié)果分析比擬為了更直觀的比擬自編程序和ANSYS的計(jì)算結(jié)果,將 錯(cuò)誤！未找到引用源.找到引用源.的數(shù)據(jù)繪入錯(cuò)誤！未找到引用源.,將錯(cuò)誤！未找到引用源.源.的數(shù)據(jù)繪入 錯(cuò)誤！未找到引用源.圖5-6應(yīng)力比擬- - stress stress(ansys)(Ed)監(jiān)a>號0D- displacement* displacement(ansy

14、s)圖5-7位移比擬自編程序所得結(jié)果與 ANSYS分析結(jié)果進(jìn)行比擬發(fā)現(xiàn),應(yīng)力偏大而位移偏小.分析其中的原因,一方面是編程水平有限,自編程序有很多不完善之處,有很多因素沒有考慮 溫度、變形、非線性等,所得結(jié)果可信度不是很高,只能得到應(yīng)力以及位移大概的分布情況；另一方面是在后處理階段,在對高斯積分點(diǎn)結(jié)果進(jìn)行處理時(shí),誤差難以防止,還有一方面的原因是在進(jìn)行求解時(shí)保存數(shù)據(jù)精度不夠,造成誤差較大,并且進(jìn)行數(shù)值運(yùn)算時(shí)可能會(huì)造成大 數(shù)吃小數(shù),從而引起結(jié)果的差異.參考文獻(xiàn)1(美)史密斯(Smith,.)著；王感等譯.有限單元法編程(第三版)M.北京：電子工業(yè)出版社,20032王勖成.有限單元法M.北京：清華大學(xué)

15、出版社,20033宋建輝,涂志剛.MATLAB 語言及其在有限元編程中的應(yīng)用J.湛江師范學(xué)院學(xué)報(bào).(24):101-1054鄭大玉,連宏玉,周躍發(fā).有限元編程與 C語言J.黑龍江商學(xué)院學(xué)報(bào).(13):23-285王偉,劉德富.有限元編程中應(yīng)用面向?qū)ο缶幊碳夹g(shù)的探討J.三峽大學(xué)學(xué)報(bào).(23):124-1286汪文立,呂士俊二級C語言程序設(shè)計(jì)教程M.北京：中國水利水電出版社,20067趙翔龍.FORTRAN 90學(xué)習(xí)教程M.北京：北京大學(xué)出版社,2002附錄程序源代碼#include<>#include<>/*The gemotry of the model*/float

16、mesh2;float wide,hight;float wsize,hsize,young,poiss,thick;int i,j,k,knelem,npoin,elem500,ielem;float coord21,props31,lnods8500,ifpre1;float posgp3,weigp3,estif136,elcod28;int npoin,nelem,kevab,igaus,jgaus,kgasp,ngash,djacb;float gpcod29,bmatx316,dmatx33;float shape8,deriv28;float xjacm22,xjaci22,el

17、cod28;float cartd28;float bmatx316,dmatx33;float nload1;float press32,pgash2,dgash2;struct node int nodenu;/*The number of the node*/float cor2;/*The coordinate of the node*/float displacement2;/*The displacement of the node*/float load2; /*The load of the node*/float boundary2; node500;struct int e

18、lementnu;/*The number of element*/int localnu8;/*Local number*/int localcorx8;int localcory8; /*Local coordinate of node*/float qx,qy,t;/*The stress and strain*/void meshing(float wide,float hight,float mesh2) printf("Please input the meshing sizen");scanf("%f%f",&wsize,&

19、hsize);getchar();mesh0=wide/wsize;mesh1=hight/hsize;/*The global coordinate of node*/void coordinate()int i,j=1;float x,y;for(i=1;i<=2*mesh1+1;i+)x=0;while(x<=wide)if(i%2!=0) nodej.cor1+=wsize/2;nodej.cor2+=(i-1)*hsize;j+; else nodej.cor1+=wsize;nodej.cor2+=(i-1)*hsize;j+; main() int top500,bo

20、ttom500;/*The number of top and bottom element*/ void input();void stif();void jacobi2();void load();void asem();void solve( );void stre();npoin=8;for(i=1;i<=8;i+)lnodsi1=i;for(i=1;i<=3;i+)scanf("%f",&propsi1);printf("The EX is %enThe uo is %8fnThe bt is %8f",props11,pr

21、ops21,props31);getch();printf("Please input the gemotry of the modeln");scanf("%f%f",&wide,&hight);getchar();printf("The wide is %,the hight is %",wide,hight);getchar();meshing(wide,hight,mesh);printf("The wide size is %f,the hight size is %fn",mesh0,m

22、esh1);input();getchar();nelem=mesh0*mesh1;for(i=1;i<=nelem;i+) /*The element number*/elementi.elementnu=i;for(i=1;i<mesh0;i+)npoin+=5;for(i=1;i<mesh1;i+)npoin+=3*mesh0+2;printf("%d",npoin);for(i=1;i<=npoin;i+)nodei.nodenu=i;for(i=1;i<=nelem;i+)for(j=1;j<=8;j+)elementi.loc

23、alnuj=j;for (i=1;i<=mesh0;i+)bottomi=i;j=1;for(i=mesh0*(mesh1-1)+1;i<=nelem;i+)topj+=i;for(i=1;i<j；i+)printf("the top number is %dn",topi);printf("%dn",element1.localnu8);coordinate();stif();jacobi2();load();asem();solve( );stre();getchar();/*data input*/void input() int

24、 nnode=8;int notreadchar;/*input element node number*/printf("input element node numbern");for(ielem=1;ielem<=nelem;ielem+)for(i=1;i<=nnode;i+)scanf("%d",&lnodsiielem);/*Input restriction message*/printf("double-digit is symbol of the restriction,1 representates st

25、able,2 representates gived displacementn");i=1;while(i) scanf("%d",&i);if(i!=0)scanf("%d",&j);scanf("%d",&ifprej);else break;/*Element stiffness matrix*/void stif() int kevab,nstre,ievab,istre,lnode,jstre,jevab;float kgasp,exisp,etasp, dvolu;float btdbm

26、,dbmat3;void sfr();/*Gauss point*/posgp1=-sqrt;posgp2=0;posgp3=sqrt;/*weight coefficient*/weigp1=;weigp2=;weigp3=;/*number of stress*/nstre=3;/*form elastic matrix*/for(ielem=1;ielem<=nelem;ielem+) printf("The number of element is %dn",ielem); for(i=1;i<=8;i+) lnode=lnodsiielem;for(j

27、=1;j<=2;j+)coord皿lnode=nodelnode.corj;elcodji=coordjlnode;young=props11;poiss=props21;thick=props31;modps(young,poiss);/*The initial value of k*/kevab=0;for(i=1;i<=16;i+)for(j=1;j<=i;j+) kevab+;estifkevab=;)/*The gauss point shape function and deriv*/kgasp=0;for(igaus=1;igaus<=3;igaus+)f

28、or(jgaus=1;jgaus<=3;jgaus+)kgasp=kgasp+1;exisp=posgpigaus;etasp=posgpjgaus;printf("The number of gauss point is %dn",kgasp);sfr2(exisp,etasp);jacob2(ielem,djacb,kgasp,elcod);dvolu=djacb*weigpigaus*weigpjgaus*thick;bmatps()/*The initial value of elastic matrix D*/kevab=0;for(ievab=1;ieva

29、b<=16;ievab+)for(istre=1;istre<=nstre;istre+)dbmatistre=;for(jstre=1;jstre<=nstre;jstre+)dbmatistre=dbmatistre+bmatxjstreievab*dmatxjstreistre;)for(jevab=1;jevab<=ievab;ievab+) kevab=kevab+1;btdbm=;for(istre=1;istre<=nstre;istre+)btdbm=btdbm+dbmatistre*bmatxistrejevab;estifkevab=estif

30、kevab+btdbm*dvolu;)/*float sfr2(float s,float t)float s2,t2,ss,tt,st,stt,sst,st2;/*Shape function*/*these variables are to simplify the formula */s2=s*;t2=t*;ss=s*s;tt=t*t;st=s*t;stt=s*t*t;sst=s*s*t;st2=s*t*;/*shape function*/shape1=+st+ss+tt-sst-stt)/;shape2=+sst)/;shape3=+ss+tt-sst+stt)/;shape4=+s

31、-tt-stt)/;shape5=+st+ss+tt+sst+stt)/;shape6=+t-ss-sst)/;shape7=+ss+tt+sst-stt)/;shape8=+stt)/;/* differential coefficient to local coordinate*/deriv11=(t+s2-st2-tt)/;deriv12=-s+st;deriv13=(-t+s2-st2+tt)/;deriv14=/;deriv15=(t+s2+st2+tt)/;deriv16=-s-st;deriv17=(-t+s2+st2-tt)/;deriv18=+tt)/;deriv21=(s+

32、t2-ss-st2)/;deriv22=+ss)/;deriv23=(-s+t2-ss+st2)/;deriv24=-t-st;deriv25=(s+t2+ss+st2)/;deriv26=/;deriv27=(-s+t2+ss-st2)/;deriv28=-t+st; */*Jacobi matrix*/void jacobi2() /* xjacm22:jacobi matrix;xjaci22:j-1;elcod28:local coordinate*/float cartd28,gpcod29;int i,j,k;for(i=1;i<=2;i+) gpcodikgasp=;for

33、(j=1;j<=8;j+)gpcodikgasp=gpcodikgasp+elcodij*shapej;for(i=1;i<=2;i+)for(j=1;j<=2;j+) xjacmij=;for(k=1;k<=8;k+)xjacmij=xjacmij+derivik*elcodjk;/*Caculate the value of Jacobi matrix*/djacb=xjacm11*xjacm22-xjacm12*xjacm21;printf("The det of jacobi matrix is %fn",djacb);if(djacb<

34、;=1e-6)printf("The jacobi det of %d is less or equal 0n",ielem);/*The element of j-1*/xjaci11=xjacm22/djacb;xjaci22=xjacm11/djacb;xjaci12=-xjacm12/djacb;xjaci21=-xjacm21/djacb;/*The deriv of shape function to global coordinate*/for(i=1;i<=2;i+)for(k=1;k<=8;k+)cartdik=;for(j=1;j<=2

35、;j+)cartdik=cartdik+xjaciij*derivjk;)/* Form elastic matris */*float modps() float bmatx316,dmatx33;float constant;int i,j,y,p;y=props11;p=props21;for(i=1;i<=3;i+)for(j=1;j<=3;j+)dmatxij=;/*If Ntype=1,it is plane stress question*/constant=y/*p);dmatx11=constant;dmatx22=constant;dmatx12=constan

36、t*p;dmatx21=constant*p;dmatx33=constant*/;)*/*Form strain matrix*/*void bmatps() float cartd28,shape8,deriv28,bmatx316,dmatx33;int npoin,nelem,ngash,mgash;float gpcod29;/*The initial value ofB*/for(i=1;i<=16;i+) for(j=1;j<=3;j+) bmatxji=;/*Form B*/ngash=0;for(i=1;i<=8;i+)(mgash=ngash+1;ngas

37、h=mgash+1;bmatx1mgash=cartd1i;bmatx1ngash=;bmatx2mgash=;bmatx2ngash=cartd2i;bmatx3mgash=cartd2i;bmatx3ngash=cartd1i; */*equal load of node*/void load()(float nload500;float elcod23,eload16,posgp3,weigp3;float press32,pgash2,dgash2,noprs3;int iedge,nedge,kount;float sgmar,tau,s,dvolu,pxcom,pycom,n,m,

38、t;int lnode,number,nloca; /*The number of load side*/*element load*/*The position of gauss point*/printf("input the number of load side*/n");scanf("%d",&number);posgp1=-sqrt;posgp2=;posgp3=sqrt;/*The weight of gauss point*/weigp1=;weigp2=;weigp3=;/*The initial load of element

39、*/for(i=1;i<=nelem;i+)nloadi=;for(i=1;i<=number;i+)/*The initial value of equal node load*/for(i=1;i<=16;i+)eloadi=0;scanf("%d",&nedge);for(iedge=1;iedge<=nedge;iedge+)for(i=1;i<=3;i+)scanf("%f%f",&sgmar,&tau);/*if sgmar=0,input node load q(1,2,3)and t(

40、1,2,3)*/if(fabs(sgmar)<1e-8)&&(fabs(tau)<1e-8)for(j=1;j<=3;j+)for(i=1;i<=2;i+)scanf("%f",&pressji);else for(i=1;i<=3;i+)pressi1=sgmar;pressi2=tau;/*The coordinate of the load node*/for(i=1;i<=3;i+) /*input load node*/printf("input node load numbern")

41、;scanf("%d",&noprsi);lnode=noprsi;for(j=1;j<=2;j+)coordjlnode=nodelnode.corj;elcodji=coordjlnode;t=;for(igaus=1;igaus<=3;igaus+) s=posgpigaus;sfr(s,t);for(j=1;j<=2;j+) pgashj=;dgashj=;/*current point q,t and the value of deriv*/ pgashj=pgashj+pressij*shapei;dgashj=dgashj+elcod

42、ji*deriv1i;dvolu=weigpigaus;pxcom=dgash1*pgash2-dgash2*pgash1;pycom=dgash1*pgash1+dgash2*pgash2;for(i=1;i<=8;i+) nloca=lnodsiielem;if(nloca=noprs1) j=i+2;break;j=i+2;kount=0;for(k=i;k<=j;k+) kount=kount+1;n=(k-1)*2+1;m=n+1;/*if the local node numner >8,then it is 1*/if(k>8) n=1;m=2;/*sum

43、 to get the equal node load*/eloadn=eloadn+shapekount*pxcom*dvolu;eloadm=eloadm+shapekount*pycom*dvolu;)print("The element load matrx is/n");for(i=1;i<=16;i+)printf("%fn",eloadi);)/*To form global stiffness matrix and load matrix*/void asem()float v1,a1;float lnods81,ifpre1,nl

44、oad1,maxa1;float ncodf2,estif136,eload16;float dlarge,x;int ipoin,kmin,lnode,khigh,kdofn,nn,nnm,nwk,iwk,in;int inodi,inode,icodf,idofn,ievab,icoln,jnode,lnodj,jdofn,jevab,jrow,ldis,lnodi;maxa1=1;for(ipoin=1;ipoin<=npoin;ipoin+) kmin=npoin+1;for(ielem=1;ielem<=nelem;ielem+)for(i=1;i<=8;i+) l

45、node=lnodsiielem;if(lnode=ipoin)break;)for(i=1;i<=8;i+) lnode=lnodsiielem;if(lnode>=kmin)break;)kmin=lnode;khigh=(ipoin-1)*2+j+1;for(j=1;j<=2;j+) kdofn=(ipoin-1)*2+j+1;maxakdofn=maxakdofn-1+khigh+j;)nn=npoin*2;nnm=nn+1;nwk=maxannm-1;printf("The sum is %dn",nwk);for(iwk=1;iwk<=n

46、wk;iwk+)aiwk=;for(in=1;in<=nn;nn+)vin=;dlarge=+15;for(ielem=1;ielem<=nelem;ielem+) printf("The number of element is %dn",ielem);for(inode=1;inode<=8;inode+) inodi=lnodsinodeielem;icodf=ifprelnodi;ncodf1=icodf/10;ncodf2=icodf-ncodf1*10;for(idofn=1;idofn<=2;idofn+) ievab=(inode-1

47、)*2+idofn;icoln=(lnodi-1)*2+idofn;for(jnode=1;jnode<=8;jnode+) lnodj=lnodsjnodeielem;for(jdofn=1;jdofn<=2;jdofn+)jevab=(jnode-1)*2+jdofn;jrow=(lnodj-1)*2+jnode;if(jrow>icoln)break;if(jevab>ievab)kevab=jevab*(jevab-1)/2+ievab;else kevab=ievab*(ievab-1)/2+jevab;ldis=maxaicoln+(icoln-jrow);

48、aldis=aldis+estifkevab;vicoln=vicoln+eloadievab;if(ncodfidofn=0) break;kevab=ievab*(ievab+1)/2;ldis=maxaicoln;aldis=aldis+dlarge*estifkevab; printf("Enter gived displacementn");for(i=1;i<=npoin;i+) icodf=ifprei;ncodf1=icodf/10;ncodf2=icodf-ncodf1*10;for(j=1;j<=2;j+) if(ncodfj=2) scan

49、f("%f",&x); icoln=(i-1)*2+j; vicoln=aldis*x; /*To solve the equation*/void solve( ) float v1,a1,maxa1;int npoin,nelem,ntype,nmats;int l,k,nn,nnm,nwk,n,ku,kl,kh,kn;float ic,b,c,klt,ki,nd,kk,m1,m2,kul;nn=npoin*2;nnm=nn+1;for(n=1;n<=nn;n+) kn=maxan;ku=maxan+1-1;kh=ku-kl;if(kh<0)if(a

50、kn<=0)printf("*Stop run! Coefficient matrix is not illegal!*Non+positive pivot for equation %d,pivot=%fn",n,akn);break;)elseif(kh=O)k=n;b=;for(kk=kl;kk<=ku;kk+)k=k-1;ki=maxak;c=akk/aki;b=b+c*akk;akk=c;akn=akn-b;)if(kh>0)k=n-kh;ic=0;klt=ku;for(j=1 ;j<=kh;j+)ic=ic+1;klt=klt-1;ki=maxak;nd=maxak+1-ki-1;if(nd<=0)k=k+1;c=0;for(l=1 ;l<=kk;l+)m1=ki+l;m2=klt+l;c=c+am1*am2;)printf("Reduce right-hand-side load vectorn");for(n=1;n<=nn;n+) kl=maxan+1;ku=maxan+1-1;kul=ku-kl;if(kul>=0) k=n;c=;for(kk=kl;kk<=ku;kk+) k=k-1;c=c+akk*vk;vn=vn-c;)printf(&quo