Triangle2D3Node有限元基础教程曾攀算例4.3.2更新内容

主函数，加入了单元应变以及单元应力的求解过程，对前述文件

http://t.csdnimg.cn/9VSZPhttp://t.csdnimg.cn/9VSZP稍作修改

%%-----主程序
clc
clear all;
 %---物理参数------------------
E=1;              %弹性模量,单位Pa
t=1;              %单元厚度,单位m
% a=1;                 %尺寸a，单位m
% q=10e6;              %均布力，单位N/m2
NU=0.25;             %泊松比
%-------------------------------
x1=0;
y1=0;
x2=1;
y2=0;
x3=1;
y3=1;
x4=0;
y4=1;
Triangle2D3Node=Triangle2D3Node_Func;
[k1,B1,D1]=Triangle2D3Node.Stiffness(E,NU,t,x1,y1,x2,y2,x4,y4,1)
[k2,B2,D2]=Triangle2D3Node.Stiffness(E,NU,t,x2,y2,x3,y3,x4,y4,1)
% 两单元刚度矩阵相同
KK=zeros(8);%（8*8）初始KK==0提供位置
KK=Triangle2D3Node.Assembly(KK,k1,1,2,4); %杆1
KK=Triangle2D3Node.Assembly(KK,k2,2,3,4)  %加杆2
% KK=Triangle2D3Node.Assembly(KK,k1,3,4,2)  %或者用这个节点顺序与k1一样
k=KK([3 4 5 6 8],[3 4 5 6 8])
p=[-1; 0; 1; 0 ;0]  %5*1  力边界条件
u=k\p   %由位移边界条件解
q=[0,0,u(1),u(2),u(3),u(4),0,u(5)]' % 全部节点位移
P=KK*q  %节点力
fp=[0 0 p(1) p(2) p(3) p(4) 0 p(5)]'
F=P-fp  %支反力
% 单元位移
q1=[q(1);q(2);q(3);q(4);q(7);q(8)]
q2=[q(3);q(4);q(5);q(6);q(7);q(8)]
%单元应变
epusilon1=B1*q1
epusilon2=B2*q2
%单元应力
sigama1=D1*B1*q1
sigama2=D2*B2*q2

功能函数：

%% 单元刚度矩阵计算函数 （输入节点、单元信息，输出单刚、弹性矩阵、应变矩阵）
function Triangle2D3Node=Triangle2D3Node_Func
Triang