Structural-Mechanical
/
FEM-Course-Matlab
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
FEM-Course-Matlab/18. 材料非线性问题matlab有限元编程/qiexiangangdu_0/qiexiangangdu/Groundoverload_2111238.m

660 lines
21 KiB
Mathematica
Raw Permalink Normal View History

init project
2024-01-28 16:46:36 +00:00
function Groundoverload_2111238(file_in)
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ó<EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ε<EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ھֲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>µı<EFBFBD><EFBFBD>κ<EFBFBD>Ӧ<EFBFBD><EFBFBD>
% ȫ<EFBFBD>ֱ<EFBFBD><EFBFBD><EFBFBD>
% gNode ------------- <EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gElement ---------- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gMaterial --------- <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gBC1 -------------- Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gDF --------------- <EFBFBD>ֲ<EFBFBD><EFBFBD><EFBFBD>
% gK ---------------- <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
% gDelta ------------ <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gNodeStress ------- <EFBFBD>ڵ<EFBFBD>Ӧ<EFBFBD><EFBFBD>
% gElementStress ---- <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD>
% gElementStrain ---- <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD>
global gDelta
if nargin < 1
file_in = 'Groundoverload.dat' ;
end
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if exist( file_in ) == 0
fprintf( '<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD> %s <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n', file_in )
fprintf( '<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֹ\n' )
return ;
end
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
[path_str,name_str,ext_str] = fileparts( file_in ) ;
ext_str_out = '.mat' ;
file_out = fullfile( path_str, [name_str, ext_str_out] ) ;
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
if exist( file_out ) ~= 0
answer = input( sprintf( '<EFBFBD>ļ<EFBFBD> %s <EFBFBD>Ѿ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڣ<EFBFBD><EFBFBD>Ƿ񸲸<EFBFBD>? ( Yes / [No] ): ', file_out ), 's' ) ;
if length( answer ) == 1
answer = 'no' ;
end
answer = lower( answer ) ;
if answer ~= 'y' | answer ~= 'yes'
fprintf( '<EFBFBD><EFBFBD>ʹ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>򱸷<EFBFBD><EFBFBD><EFBFBD><EFBFBD>е<EFBFBD><EFBFBD>ļ<EFBFBD>\n' ) ;
fprintf( '<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֹ\n' ) ;
return ;
end
end
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD>Ͳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>⣬<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
FemModel( file_in ) ; % <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD>
SolveModel ; % <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD>
SaveResults( file_out ) ; % <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
DisplayModel;
PlotDisplacement( 2 );
PlotDisplacementContour( 2, 10, 'r' );
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
fprintf( '<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD><EFBFBD><EFBFBD>ֱλ<EFBFBD><EFBFBD> %f m\n', full( gDelta(170,1) ) ) ;
fprintf( '<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD> %s <EFBFBD><EFBFBD>\n', file_out ) ;
return ;
function FemModel(filename)
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD>
% <EFBFBD>ú<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ݣ<EFBFBD>
% gNode ------- <EFBFBD>ڵ㶨<EFBFBD><EFBFBD>
% gElement ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gMaterial --- <EFBFBD><EFBFBD><EFBFBD>϶<EFBFBD><EFBFBD>壬<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ģ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ľ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ŀ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ծ<EFBFBD>
% gBC1 -------- Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% gDF --------- <EFBFBD>ֲ<EFBFBD><EFBFBD><EFBFBD>
global gNode gElement gMaterial gBC1 gDF
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD>
fid = fopen( filename, 'r' ) ;
% <EFBFBD><EFBFBD>ȡ<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
node_number = fscanf( fid, '%d', 1 ) ;
gNode = zeros( node_number, 2 ) ;
for i=1:node_number
dummy = fscanf( fid, '%d', 1 ) ;
gNode( i, : ) = fscanf( fid, '%f', [1, 2] ) ;
end
% <EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
element_number = fscanf( fid, '%d', 1 ) ;
gElement = zeros( element_number, 4 ) ;
for i=1:element_number
dummy = fscanf( fid, '%d', 1 ) ;
gElement( i, : ) = fscanf( fid, '%d', [1, 4] ) ;
end
% <EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϣ
load_number = fscanf( fid, '%d', 1 ) ;
gDF = zeros( load_number, 4 ) ;
for i=1:load_number
dummy = fscanf( fid, '%d', 1 ) ;
gDF( i, : ) = fscanf( fid, '%f', [1, 4] ) ;
end
% <EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϣ
material_number = fscanf( fid, '%d', 1 ) ;
gMaterial = zeros( material_number, 4 ) ;
for i=1:material_number
dummy = fscanf( fid, '%d', 1 ) ;
gMaterial( i, : ) = fscanf( fid, '%f', [1,4] ) ;
end
% <EFBFBD><EFBFBD>ȡ<EFBFBD>߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
bc1_number = fscanf( fid, '%d', 1 ) ;
gBC1 = zeros( bc1_number, 3 ) ;
for i=1:bc1_number
gBC1( i, 1 ) = fscanf( fid, '%d', 1 ) ;
gBC1( i, 2 ) = fscanf( fid, '%d', 1 ) ;
gBC1( i, 3 ) = fscanf( fid, '%f', 1 ) ;
end
% <EFBFBD>ر<EFBFBD><EFBFBD>ļ<EFBFBD>
fclose( fid ) ;
return
function SolveModel
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD>
% <EFBFBD>ú<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD>ͣ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% 1. <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD>󣬼<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
% 2. <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>ĵ<EFBFBD>Ч<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% 3. <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>޸<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͽڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% 4. <EFBFBD><EFBFBD><EFBFBD>ⷽ<EFBFBD><EFBFBD><EFBFBD>飬<EFBFBD>õ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% 5. <EFBFBD><EFBFBD><EFBFBD>㵥ԪӦ<EFBFBD><EFBFBD><EFBFBD>ͽڵ<EFBFBD>Ӧ<EFBFBD><EFBFBD>
global gNode gElement gMaterial gBC1 gK gDelta gNodeStress gElementStress gDF gElementStrain gFE
%% step 1. <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͽڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
[node_number,dummy] = size( gNode ) ;
gK = sparse( node_number * 2, node_number * 2 ) ; %ϵͳ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
gFE = sparse( node_number * 2, 1 ) ; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
f = sparse( node_number * 2, 1 ) ; %<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>Ч<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% step 2. <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʼȫ<EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
[element_number, dunmmy] = size( gElement ) ;
gElementStrain = zeros( element_number, 3) ; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
gElementStress = zeros( element_number, 3) ; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
%% step 3. <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>泬<EFBFBD>ز<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ĵ<EFBFBD>Ч<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
[df_number,dummy] = size( gDF ) ;
for idf = 1:1:df_number
enf = EquivalentNodeForce( gDF(idf,1), gDF(idf,2), gDF(idf,3), gDF(idf,4) )*20;% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>10<EFBFBD>Ǹı<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ĵ<EFBFBD>С
i = gElement( gDF(idf,1), 1 ) ;
j = gElement( gDF(idf,1), 2 ) ;
m = gElement( gDF(idf,1), 3 ) ;
f( (i-1)*2+1 : (i-1)*2+2 ) = f( (i-1)*2+1 : (i-1)*2+2 ) + enf( 1:2 ) ;
f( (j-1)*2+1 : (j-1)*2+2 ) = f( (j-1)*2+1 : (j-1)*2+2 ) + enf( 3:4 ) ;
f( (m-1)*2+1 : (m-1)*2+2 ) = f( (m-1)*2+1 : (m-1)*2+2 ) + enf( 5:6 ) ;
end
%% step 5 <EFBFBD><EFBFBD><EFBFBD>߸նȷ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
gDelta=zeros(node_number * 2,1); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> ȡ<EFBFBD><EFBFBD>ֵdelta0=0
js=0; %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͳ<EFBFBD><EFBFBD>ָ<EFBFBD><EFBFBD>
while true
gK=zeros( node_number * 2, node_number * 2 ); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
gFE=zeros(node_number * 2,1); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ч<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for ie=1:1:element_number
delta = NodeDe( ie,gDelta); %<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ȡ<EFBFBD><EFBFBD>Ԫ<EFBFBD>Ľڵ<EFBFBD>λ<EFBFBD><EFBFBD>
eps = MatrixB( ie ) * delta; %<EFBFBD><EFBFBD><EFBFBD>η<EFBFBD><EFBFBD>̣<EFBFBD><EFBFBD><EFBFBD>ʽ2<EFBFBD><EFBFBD>epsilon0
sigma0 = unlinerD(ie,eps)* (eps-gElementStrain(ie,:)')+gElementStress(ie,:)';%g<EFBFBD><EFBFBD>ʽ3<EFBFBD><EFBFBD>sigmma0
h = gMaterial( gElement(ie, 4), 3 ) ;
B = MatrixB( ie );
area = ElementArea( ie );
kt = transpose(B)*unlinerD(ie,eps)*B*h*abs(area) ; %<EFBFBD><EFBFBD>ʽ4
AssembleStiffnessMatrix( ie, kt ) ;
for ii = 1:1:3
gElementStrain(ie,ii) = eps(ii);
gElementStress(ie,ii) = sigma0(ii);
end
FE = elementforce( ie ,gElementStress) ;% <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
AssembleFE( ie, FE ) ; % <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>gFE
end
Phi=( f-gFE );
%% step 4. <EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>޸ĸնȾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͽڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ó˴<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
[bc_number,dummy] = size( gBC1 ) ;
for ibc=1:1:bc_number
n = gBC1(ibc, 1 ) ;
d = gBC1(ibc, 2 ) ;
m = (n-1)*2 + d ;
Phi(m) = gBC1(ibc, 3)* gK(m,m) * 1e20 ;
gK(m,m) = gK(m,m) * 1e20 ;
end
dd = gK \ Phi; %<EFBFBD><EFBFBD>ʽ5 λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
gDelta=dd+gDelta;
conv=norm(Phi)/norm(f);
fprintf('<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%d <EFBFBD><EFBFBD>ƽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>/<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>%e \n',js,conv)
if js>100 || conv<1e-8 %<EFBFBD>ж<EFBFBD><EFBFBD>Ƿ<EFBFBD><EFBFBD>ﵽ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ҫ<EFBFBD><EFBFBD>
break
else
js=js+1;
end
end
%% step 6. <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD>Ӧ<EFBFBD><EFBFBD>(<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ƽڵ<EFBFBD><EFBFBD><EFBFBD>Ȩƽ<EFBFBD><EFBFBD>)
gNodeStress = zeros( node_number, 6 ) ;
for i=1:node_number
S = zeros( 1, 3 ) ;
A = 0 ;
for ie=1:1:element_number
for k=1:1:3
if i == gElement( ie, k )
area= ElementArea( ie ) ;
S = S + gElementStress(ie,1:3 ) * area ;
A = A + area ;
break ;
end
end
end
gNodeStress(i,1:3) = S / A ;
gNodeStress(i,6) = 0.5*sqrt( (gNodeStress(i,1)-gNodeStress(i,2))^2 + 4*gNodeStress(i,3)^2 ) ;
gNodeStress(i,4) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) + gNodeStress(i,6) ;
gNodeStress(i,5) = 0.5*(gNodeStress(i,1)+gNodeStress(i,2)) - gNodeStress(i,6) ;
end
return
function B = MatrixB( ie )
% <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD><EFBFBD>Ӧ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>B
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% B ---- <EFBFBD><EFBFBD>ԪӦ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
global gNode gElement
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai = xj*ym - xm*yj ;
aj = xm*yi - xi*ym ;
am = xi*yj - xj*yi ;
bi = yj - ym ;
bj = ym - yi ;
bm = yi - yj ;
ci = -(xj-xm) ;
cj = -(xm-xi) ;
cm = -(xi-xj) ;
area = abs((ai+aj+am)/2) ;
B = [bi 0 bj 0 bm 0
0 ci 0 cj 0 cm
ci bi cj bj cm bm] ;
B = B/2/area ;
return
function area = ElementArea( ie )
% <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% area ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
global gNode gElement
xi = gNode( gElement( ie, 1 ), 1 ) ;
yi = gNode( gElement( ie, 1 ), 2 ) ;
xj = gNode( gElement( ie, 2 ), 1 ) ;
yj = gNode( gElement( ie, 2 ), 2 ) ;
xm = gNode( gElement( ie, 3 ), 1 ) ;
ym = gNode( gElement( ie, 3 ), 2 ) ;
ai = xj*ym - xm*yj ;
aj = xm*yi - xi*ym ;
am = xi*yj - xj*yi ;
area = abs((ai+aj+am)/2) ;
return
function k = StiffnessMatrix( ie )
% <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% k ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
global gElement gMaterial
k = zeros( 6, 6 ) ;
h = gMaterial( gElement(ie, 4), 3 ) ;
B = MatrixB( ie );
area = ElementArea( ie );
E = gMaterial( gElement(ie, 4), 1 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
D = [ 1-mu mu 0
mu 1-mu 0
0 0 (1-2*mu)/2] ;
D = D*E/(1-2*mu)/(1+mu) ;
k = transpose(B)*D*B*h*abs(area) ;
return
function FE = elementforce( ie ,gElementStress)
% <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% FE ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
global gElement gMaterial
FE = zeros( 6, 6 ) ;
thickness = gMaterial( gElement(ie, 4), 3 ) ;
area = ElementArea(ie);
B = MatrixB( ie );
sigma0 = [gElementStress(ie,1);gElementStress(ie,2);gElementStress(ie,3)];
FE = transpose(B)*sigma0*thickness*abs(area) ;
return
function FEE = AssembleFE(ie, FE)
global gElement gFE
for i = 1:1:3
gFE(gElement(ie,i)*2-1)=gFE(gElement(ie,i)*2-1)+FE(i*2-1);
gFE(gElement(ie,i)*2)=gFE(gElement(ie,i)*2)+FE(i*2);
end
FEE=0;
return
function D = unlinerD (ie,eps)
%<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ե<EFBFBD><EFBFBD><EFBFBD>D<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie ---- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% D ---- D<EFBFBD><EFBFBD><EFBFBD><EFBFBD>
global gElement gMaterial
E = gMaterial( gElement(ie, 4), 1 ) ;
mu = gMaterial( gElement(ie, 4), 2 ) ;
D = [ 1-mu mu 0
mu 1-mu 0
0 0 (1-2*mu)/2] ;
epsx = eps(1);
epsy = eps(2);
D = D*E*(1-100*epsx^2-100*epsy^2)/(1-2*mu)/(1+mu) ;
return
function delta = NodeDe( ie ,gDelta)
% <EFBFBD><EFBFBD><EFBFBD>㵥Ԫ<EFBFBD>ڵ<EFBFBD>λ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
global gElement
delta = zeros( 6, 1 ) ;
for j=1:1:3
delta( 2*j-1 ) = gDelta( 2*gElement( ie, j )-1 ) ;
delta( 2*j ) = gDelta( 2*gElement( ie, j ) ) ;
end
return
function AssembleStiffnessMatrix( ie, k )
% <EFBFBD>ѵ<EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD>󼯳ɵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie --- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% k --- <EFBFBD><EFBFBD>Ԫ<EFBFBD>նȾ<EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% <EFBFBD><EFBFBD>
global gElement gK
for i=1:1:3
for j=1:1:3
for p=1:1:2
for q=1:1:2
m = (i-1)*2+p ;
n = (j-1)*2+q ;
M = (gElement(ie,i)-1)*2+p ;
N = (gElement(ie,j)-1)*2+q ;
gK(M,N) = gK(M,N) + k(m,n) ;
end
end
end
end
return
function enf = EquivalentNodeForce( ie, aa, bb, ps )
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Էֲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>صĵ<EFBFBD>Ч<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>:
% ie ----- <EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>
% aa ----- <EFBFBD>յ<EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
% bb ----- <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD>
% ps ----- <EFBFBD><EFBFBD><EFBFBD>طֲ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ:
% enf ----- <EFBFBD><EFBFBD>Ч<EFBFBD>ڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
global gElement gNode gMaterial
enf = zeros(6,1) ;
h = gMaterial( gElement( ie, 4 ), 3 ) ;
xj = gNode( aa, 1 ) ;
yj = gNode( aa, 2 ) ;
xi = gNode( bb, 1 ) ;
yi = gNode( bb, 2 ) ;
f1 = h*ps*(yi-yj)/2 ;
f2 = h*ps*(xj-xi)/2 ;
f3 = h*ps*(yi-yj)/2 ;
f4 = h*ps*(xj-xi)/2 ;
% <EFBFBD>ֲ<EFBFBD>ǰ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>еĽڵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
switch ie
case 446
enf = [f1; f2; f3; f4; 0; 0] ;
case {492, 449, 456, 458, 460, 462}
enf = [0; 0; f1; f2; f3; f4] ;
case 437
enf = [f3; f4; 0; 0; f1; f2] ;
end
return
function SaveResults( file_out )
% <EFBFBD><EFBFBD>ʾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% file_out --- <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ļ<EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>
global gNode gElement gMaterial gBC1 gDelta gNodeStress gElementStress gK
save( file_out, 'gNode', 'gElement', 'gMaterial', 'gBC1', 'gDelta', 'gNodeStress', 'gElementStress', 'gK' ) ;
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DisplayModel
% <EFBFBD><EFBFBD>ͼ<EFBFBD>η<EFBFBD>ʽ<EFBFBD><EFBFBD>ʾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>
global gNode gElement gMaterial gBC1
figure ;
axis equal ;
axis off ;
set( gcf, 'NumberTitle', 'off' ) ;
set( gcf, 'Name', '<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD><EFBFBD>' ) ;
% <EFBFBD><EFBFBD><EFBFBD>ݲ<EFBFBD>ͬ<EFBFBD>IJ<EFBFBD><EFBFBD>ϣ<EFBFBD><EFBFBD><EFBFBD>ʾ<EFBFBD><EFBFBD>Ԫ<EFBFBD><EFBFBD>ɫ
[element_number, dummy] = size( gElement ) ;
material_color = [ 'r','g','b','c','m','y','w','k'] ;
for i=1:element_number
x1 = gNode( gElement( i, 1 ), 1 ) ;
x2 = gNode( gElement( i, 2 ), 1 ) ;
x3 = gNode( gElement( i, 3 ), 1 ) ;
x4 = gNode( gElement( i, 4 ), 1 ) ;
y1 = gNode( gElement( i, 1 ), 2 ) ;
y2 = gNode( gElement( i, 2 ), 2 ) ;
y3 = gNode( gElement( i, 3 ), 2 ) ;
y4 = gNode( gElement( i, 4 ), 2 ) ;
color_index = mod( 1, length( material_color ) ) ;
if color_index == 0
color_index = length( material_color ) ;
end
patch( [x1;x2;x3], [y1;y2;y3], material_color( color_index ) ) ;
end
% <EFBFBD><EFBFBD>ʾ<EFBFBD>߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
DisplayBC( 'blue' ) ;
return
function DisplayBC( color )
% <EFBFBD><EFBFBD>ͼ<EFBFBD>η<EFBFBD>ʽ<EFBFBD><EFBFBD>ʾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ԫģ<EFBFBD>͵ı߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% color ---- <EFBFBD>߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɫ
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>
global gNode gBC1
% ȷ<EFBFBD><EFBFBD><EFBFBD>߽<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ĵ<EFBFBD>С
xmin = min( gNode(:,1) ) ;
xmax = max( gNode(:,1) ) ;
factor = ( xmax - xmin ) / 25 ;
[bc1_number,dummy] = size( gBC1 ) ;
dBCSize = factor ;
for i=1:bc1_number
if( gBC1( i, 2 ) == 1 ) % x<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>
x0 = gNode( gBC1( i, 1 ), 1 ) ;
y0 = gNode( gBC1( i, 1 ), 2 ) ;
x1 = x0 - dBCSize ;
y1 = y0 + dBCSize/2 ;
x2 = x1 ;
y2 = y0 - dBCSize/2 ;
hLine = line( [x0 x1 x2 x0], [y0 y1 y2 y0] ) ;
set( hLine, 'Color', color ) ;
xCenter = x1 - dBCSize/6 ;
yCenter = y0 + dBCSize/4 ;
radius = dBCSize/6 ;
theta=0:pi/6:2*pi ;
x = radius * cos( theta ) ;
y = radius * sin( theta ) ;
hLine = line( x+xCenter, y+yCenter ) ;
set( hLine, 'Color', color ) ;
hLine = line( x+xCenter, y+yCenter-dBCSize/2 ) ;
set( hLine, 'Color', color ) ;
x0 = x0 - dBCSize - dBCSize/3 ;
y0 = y0 + dBCSize/2 ;
x1 = x0 ;
y1 = y0 - dBCSize ;
hLine = line( [x0, x1], [y0, y1] ) ;
set( hLine, 'Color', color ) ;
x = [x0 x0-dBCSize/6] ;
y = [y0 y0-dBCSize/6] ;
hLine = line( x, y ) ;
set( hLine, 'Color', color ) ;
for j=1:1:4
hLine = line( x, y - dBCSize/4*j );
set( hLine, 'Color', color ) ;
end
else % y<EFBFBD><EFBFBD><EFBFBD><EFBFBD>Լ<EFBFBD><EFBFBD>
x0 = gNode( gBC1( i, 1 ), 1 ) ;
y0 = gNode( gBC1( i, 1 ), 2 ) ;
x1 = x0 - dBCSize/2 ;
y1 = y0 - dBCSize ;
x2 = x1 + dBCSize ;
y2 = y1 ;
hLine = line( [x0 x1 x2 x0], [y0 y1 y2 y0] ) ;
set( hLine, 'Color', color ) ;
xCenter = x0 - dBCSize/4 ;
yCenter = y1 - dBCSize/6 ;
radius = dBCSize/6 ;
theta=0:pi/6:2*pi ;
x = radius * cos( theta ) ;
y = radius * sin( theta ) ;
hLine = line( x+xCenter, y+yCenter ) ;
set( hLine, 'Color', color ) ;
hLine = line( x+xCenter+dBCSize/2, y+yCenter ) ;
set( hLine, 'Color', color ) ;
hLine = line( [x1, x1+dBCSize], [y1-dBCSize/3, y1-dBCSize/3] ) ;
set( hLine, 'Color', color ) ;
x = [x1 x1-dBCSize/6] ;
y = [y1-dBCSize/3 y1-dBCSize/2] ;
hLine = line( x, y ) ;
set( hLine, 'Color', color ) ;
for j=1:1:4
hLine = line( x+dBCSize/4*j, y );
set( hLine, 'Color', color ) ;
end
end
end
return
function PlotDisplacement( iDisp )
% <EFBFBD><EFBFBD>ʾλ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͼ
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% iDisp --- λ<EFBFBD>Ʒ<EFBFBD><EFBFBD><EFBFBD>ָʾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
% 1 -- x<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
% 2 -- y<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>
global gNode gElement gDelta
switch iDisp
case 1
title = ' x <EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>' ;
case 2
title = ' y <EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>' ;
otherwise
fprintf( 'λ<EFBFBD>Ʒ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>\n' ) ;
return ;
end
figure ;
axis equal ;
axis off ;
set( gcf, 'NumberTitle', 'off' ) ;
set( gcf, 'Name', title ) ;
dispMin = min( gDelta( iDisp:2:length(gDelta) ) ) ;
dispMax = max( gDelta( iDisp:2:length(gDelta) ) ) ;
caxis( [dispMin, dispMax] ) ;
colormap( 'jet' ) ;
[element_number, dummy] = size( gElement ) ;
for ie=1:1:element_number
x = [ gNode( gElement( ie, 1 ), 1 ) ;
gNode( gElement( ie, 2 ), 1 ) ;
gNode( gElement( ie, 3 ), 1 ) ] ;
y = [ gNode( gElement( ie, 1 ), 2 ) ;
gNode( gElement( ie, 2 ), 2 ) ;
gNode( gElement( ie, 3 ), 2 ) ] ;
c = [ gDelta( ( gElement( ie, 1 ) - 1 ) * 2 + iDisp ) ;
gDelta( ( gElement( ie, 2 ) - 1 ) * 2 + iDisp ) ;
gDelta( ( gElement( ie, 3 ) - 1 ) * 2 + iDisp )] ;
set( patch( x, y, c ), 'EdgeColor', 'interp' ) ;
end
yTick = dispMin:(dispMax-dispMin)/10:dispMax ;
Label = cell( 1, length(yTick) );
for i=1:length(yTick)
Label{i} = sprintf( '%.2e', yTick(i) ) ;
end
set( colorbar( 'vert' ), 'YTick', yTick, 'YTickLabelMode', 'Manual', 'YTickLabel', Label ) ;
PlotDisplacementContour( iDisp, 10, 'white' ) ;
return
function PlotDisplacementContour( iDisp, nContour, color )
% <EFBFBD><EFBFBD>ʾλ<EFBFBD>Ƶ<EFBFBD>ֵ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% iDisp ----- λ<EFBFBD>Ʒ<EFBFBD><EFBFBD><EFBFBD>ָʾ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ
% 1 -- x<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
% 2 -- y<EFBFBD><EFBFBD><EFBFBD><EFBFBD>λ<EFBFBD><EFBFBD>
% nContour -- <EFBFBD><EFBFBD>ֵ<EFBFBD>ߵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
% color ---- <EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɫ
% <EFBFBD><EFBFBD><EFBFBD><EFBFBD>ֵ<EFBFBD><EFBFBD>
% <EFBFBD><EFBFBD>
global gNode gElement gDelta
[element_number, dummy] = size( gElement ) ;
[node_number, dummy] = size( gNode ) ;
dispMin = min( gDelta( iDisp:2:length(gDelta) ) ) ;
dispMax = max( gDelta( iDisp:2:length(gDelta) ) ) ;
dispDelta = (dispMax-dispMin)/( nContour+1 ) ;
dispContour = dispMin+dispDelta : dispDelta : dispMax - dispDelta ;
for ie=1:1:element_number
x = [ gNode( gElement( ie, 1 ), 1 ) ;
gNode( gElement( ie, 2 ), 1 ) ;
gNode( gElement( ie, 3 ), 1 ) ] ;
y = [ gNode( gElement( ie, 1 ), 2 ) ;
gNode( gElement( ie, 2 ), 2 ) ;
gNode( gElement( ie, 3 ), 2 )] ;
s = [ gDelta( ( gElement( ie, 1 ) - 1 ) * 2 + iDisp ) ;
gDelta( ( gElement( ie, 2 ) - 1 ) * 2 + iDisp ) ;
gDelta( ( gElement( ie, 3 ) - 1 ) * 2 + iDisp ) ] ;
for is = 1:1:nContour
[smax, ismax] = max( s ) ;
[smin, ismin] = min( s ) ;
if dispContour(is) > smax || dispContour(is) < smin
continue ;
end
x1 = x(ismin) + ( dispContour(is)- smin ) / (smax-smin) * ( x(ismax) - x(ismin) ) ;
y1 = y(ismin) + ( dispContour(is)- smin ) / (smax-smin) * ( y(ismax) - y(ismin) ) ;
for ismed=1:1:4
if ismed ~= ismax && ismed ~= ismin
break ;
end
end
if dispContour(is) < s( ismed )
x2 = x(ismin) + (dispContour(is)-smin)/(s(ismed)-smin)*(x(ismed)-x(ismin)) ;
y2 = y(ismin) + (dispContour(is)-smin)/(s(ismed)-smin)*(y(ismed)-y(ismin)) ;
else
x2 = x(ismed) + (dispContour(is)-s(ismed))/(smax-s(ismed))*(x(ismax)-x(ismed)) ;
y2 = y(ismed) + (dispContour(is)-s(ismed))/(smax-s(ismed))*(y(ismax)-y(ismed)) ;
end
set( line( [x1;x2], [y1;y2] ), 'color', color ) ;
end
end
return