三角形有限元法程序4等参元

弹性矩阵弹性矩阵0 D 01称10 0 0应变矩阵：B1i2Aci e单元刚度矩阵K BTF=FeK=Ke网格划分及数据文件准备整整体控制方程 Ka引入位移边界条 求解整体控制方

%单元结点%单元结点编%node1,node2,node3,gElt=[ ielt=size(gElt,1);%单元总数forie=1:ieltgEInfo(ie)1;各单元的材料%边界条件:type=1--fixeddisp[1--结点固定%globalgMatgNdtgEltgEInfo%材料 E Poissonratio, gMat=%结点 gNdt= e+000 e+000 [函数%0--%0--结点自由];type=2--point% y, *, gBct=[1.01.00.00.021.01.00.00.00.0-1.00.00.0DD0001称10 0 [函数 % 0[函数 %%%%输入数据 给定点的s、t%输出数据 % Ni？ y Ni Ni x [[函数%计算单元的刚度矩KeBT1=F,KeBT1=F,Jd HiHjFi,j j%%实 积分，计算 积分子函

[函数globalgEltgKAenm %单元的结点总iLR %单元的总自由度

for inod= viDof(2*il-1:2*il)=2*inod- end%for %%%ir为行局部编号,irg为对应整体

for 520002000 0110 0110 0002 0310134 4

%ic为列局部编号,icg为对应整体编forgKA(irg,icg)=

[函数global%计算并组集整体刚度forie=1:iLenEltifelseifAssembleEleStif(ie,Ke);%[函数globalgEltgBctgKAiBct=size(gBct,1);forib=1:iBctiNode=iType=if forif(gBct(ib,ind+1)==1) [函数globalgBctgKAgFAiBct=size(gBct,1);forib=1:iBctiNode=iType=ifiType==2 forind=1:2irow=(iNode-gFA(irow)=[函数functionfunctionglobalgNdtgElthold %变形图放大%绘制变形前view(2);axisequal;axisoff;axisview(2);axisequal;axisoff;axistight;－求解主程 %定义全局%定义全局globalgMatgNdtgEltgEInfogBctglobalgKAgFAglobal %直接输入%Initializetheresultdata

%形成LoopCalcGlobalStif;%计算并组集整体刚度LoopBoundaryLoad;%处理荷载边界条件LoopBoundaryDisp;%处理位移边界条件%求解forid=1:iNdt%在屏幕上绘制－练习题在三角形有限元程序的基础上，编制四边形等参单元的程序，进下面算例的求图示悬臂梁(8m*1m)，左端固定，