地球物理中的有限单元法-第一类边界条件-三角剖分-二次插值 matlab编程实现

Li_Fu@csu

已于 2022-09-02 20:59:07 修改

阅读量606

点赞数 1

分类专栏：有限单元法文章标签： matlab

于 2022-09-02 20:57:46 首次发布

本文链接：https://blog.csdn.net/weixin_46973478/article/details/126670285

版权

有限单元法专栏收录该内容

3 篇文章 0 订阅

订阅专栏

懵懵懂懂，太难了

参考资料：《地球物理中的有限单元法》徐世浙

% ---------------------------------------------
% 主函数
% 第一类边界条件、二次插值、三角单元剖分
% 编写时间：2022-08-30
% ---------------------------------------------

clear all;clc;

%% 输入原始数据
ND = 45; % 节点总数
NE = 16; % 单元总数
ND1 = 24; % 第一类边界节点数

% 单元节点编号
fid1=fopen('单元节点编号2.dat','rt');
line1=fgets(fid1);
data1 = fscanf(fid1,'%e %e %e %e %e %e %e',[7 inf]);
fclose(fid1);
data1 = data1';
NEi = data1(:,1)'; % 单元序号
i1 = data1(:,2)';
i2 = data1(:,3)';
i3 = data1(:,4)';
i4 = data1(:,5)';
i5 = data1(:,6)';
i6 = data1(:,7)';
I6 = [i1;i2;i3;i4;i5;i6]; % 单元节点编号 [6,NE]

% 节点坐标
fid2 = fopen('节点坐标2.dat','rt');
line2 = fgets(fid2);
data2 = fscanf(fid2,'%e %e %e',[3 inf]);
fclose(fid2);
data2 = data2';
NDi = data2(:,1)'; % 节点序号
x = data2(:,2)';
y = data2(:,3)';
XY = [x;y]; % 节点坐标[2,ND]

% 计算节点坐标（存放三边中点节点坐标原始数据定义为0）
% 调用中点节点坐标计算函数
XY = XY_cal(XY);

% 第一类边界节点号和场值
fid3 = fopen('第一类边界节点号和场值2.dat','rt');
line3 = fgets(fid3);
data3 = fscanf(fid3,'%e %e %e',[3 inf]);
data3 = data3';
NBi = data3(:,1)'; % 第一类边界节点序号
NB1 = data3(:,2)'; % 第一类边界节点号 [1,ND1]
U1 = data3(:,3)'; % 第一类边界节点场值 [1,ND1]

[ID,NP] = FID(ND,NE,I6); % 形成ID
[GA] = UK2(NE,NP,I6,XY,ID); % 形成总体系数矩阵
[GA,U]=UB12(ND1,NB1,U1,ND,ID,GA); % 代入第一类边界条件
[B,IR] = BANDV(ND,GA,ID,U); % 解线性方程组

% ---------------------------------------------
% 程序说明
% 三角形中点节点坐标计算
% 输入三角形顶点节点编号及其坐标
% 变量说明
% XY 输入输出变量节点坐标[2,ND]
% NDi 节点序号
% ---------------------------------------------

function [XY] = XY_cal(XY)
XY(1,2) = (XY(1,1)+XY(1,3))/2;
XY(2,2) = (XY(2,1)+XY(2,3))/2;
XY(1,4) = (XY(1,3)+XY(1,5))/2;
XY(2,4) = (XY(2,3)+XY(2,5))/2;
XY(1,6) = (XY(1,1)+XY(1,11))/2;
XY(2,6) = (XY(2,1)+XY(2,11))/2;
XY(1,7) = (XY(1,3)+XY(1,11))/2;
XY(2,7) = (XY(2,3)+XY(2,11))/2;
XY(1,8) = (XY(1,3)+XY(1,13))/2;
XY(2,8) = (XY(2,3)+XY(2,13))/2;
XY(1,9) = (XY(1,5)+XY(1,13))/2;
XY(2,9) = (XY(2,5)+XY(2,13))/2;
XY(1,10) = (XY(1,5)+XY(1,15))/2;
XY(2,10) = (XY(2,5)+XY(2,15))/2;

XY(1,12) = (XY(1,11)+XY(1,13))/2;
XY(2,12) = (XY(2,11)+XY(2,13))/2;
XY(1,14) = (XY(1,13)+XY(1,15))/2;
XY(2,14) = (XY(2,13)+XY(2,15))/2;
XY(1,16) = (XY(1,11)+XY(1,21))/2;
XY(2,16) = (XY(2,11)+XY(2,21))/2;
XY(1,17) = (XY(1,13)+XY(1,21))/2;
XY(2,17) = (XY(2,13)+XY(2,21))/2;
XY(1,18) = (XY(1,13)+XY(1,23))/2;
XY(2,18) = (XY(2,13)+XY(2,23))/2;
XY(1,19) = (XY(1,15)+XY(1,23))/2;
XY(2,19) = (XY(2,15)+XY(2,23))/2;
XY(1,20) = (XY(1,15)+XY(1,25))/2;
XY(2,20) = (XY(2,15)+XY(2,25))/2;

XY(1,22) = (XY(1,21)+XY(1,23))/2;
XY(2,22) = (XY(2,21)+XY(2,23))/2;
XY(1,24) = (XY(1,23)+XY(1,25))/2;
XY(2,24) = (XY(2,23)+XY(2,25))/2;
XY(1,26) = (XY(1,21)+XY(1,31))/2;
XY(2,26) = (XY(2,21)+XY(2,31))/2;
XY(1,27) = (XY(1,21)+XY(1,33))/2;
XY(2,27) = (XY(2,21)+XY(2,33))/2;
XY(1,28) = (XY(1,23)+XY(1,33))/2;
XY(2,28) = (XY(2,23)+XY(2,33))/2;
XY(1,29) = (XY(1,23)+XY(1,35))/2;
XY(2,29) = (XY(2,23)+XY(2,35))/2;
XY(1,30) = (XY(1,25)+XY(1,35))/2;
XY(2,30) = (XY(2,25)+XY(2,35))/2;

XY(1,32) = (XY(1,31)+XY(1,33))/2;
XY(2,32) = (XY(2,31)+XY(2,33))/2;
XY(1,34) = (XY(1,33)+XY(1,35))/2;
XY(2,34) = (XY(2,33)+XY(2,35))/2;
XY(1,36) = (XY(1,31)+XY(1,41))/2;
XY(2,36) = (XY(2,31)+XY(2,41))/2;
XY(1,37) = (XY(1,31)+XY(1,43))/2;
XY(2,37) = (XY(2,31)+XY(2,43))/2;
XY(1,38) = (XY(1,33)+XY(1,43))/2;
XY(2,38) = (XY(2,33)+XY(2,43))/2;
XY(1,39) = (XY(1,33)+XY(1,45))/2;
XY(2,39) = (XY(2,33)+XY(2,45))/2;
XY(1,40) = (XY(1,35)+XY(1,45))/2;
XY(2,40) = (XY(2,35)+XY(2,45))/2;

XY(1,42) = (XY(1,41)+XY(1,43))/2;
XY(2,42) = (XY(2,41)+XY(2,43))/2;
XY(1,44) = (XY(1,43)+XY(1,45))/2;
XY(2,44) = (XY(2,43)+XY(2,45))/2;

end

% ---------------------------------------------
% 程序说明
% 三角单元、二次插值
% 功能：在确定单元节点编号后，确定ID矩阵
% 变量说明：
% ND 节点总数
% NE 单元总数
% I6 存放单元节点编号，[6×NE]
% ID 输出参数，第i行对角元Kii在GA向量中的位置j 表示ID(i)=j
% NP 输出参数，GA 向量的长度
% ---------------------------------------------

function [ID,NP] = FID(ND,NE,I6)
ID = zeros(1,ND);
ID(1) = 1;
for I = 2:ND
NL = 0;
for L = 1:NE

for J = 1:6
if(I == I6(J,L))
Temp = 1; break;
else
Temp = 0;
end
end
switch Temp
case 1
L1 = 0;
for K = 1:6
LL = I-I6(K,L);
if(LL>L1)
L1 = LL;
end
end
if(L1>NL)
NL = L1;
end
case 0
continue;
end

end
ID(I) = ID(I-1)+NL+1;
end
NP = ID(ND);

end

% ---------------------------------------------
% 程序说明
% 三角单元、二次插值时计算单元系数矩阵 ke
% 变量说明：
% X 存放三角元顶点的x坐标，存放顺序x1、x2、x3
% Y 存放三角元顶点的y坐标，存放顺序y1、y2、y3
% KE 存放单元系数矩阵，[6×6]
% ---------------------------------------------

function KE = UKE2(X,Y)
KE = zeros(6,6);

A(1) = Y(2)-Y(3);
A(2) = Y(3)-Y(1);
A(3) = Y(1)-Y(2);
B(1) = X(3)-X(2);
B(2) = X(1)-X(3);
B(3) = X(2)-X(1);

S = 6*(A(1)*B(2)-A(2)*B(1));

KE(1,1) = 3*(A(1)*B(1)+A(1)*B(1))/S;
KE(2,1) = -(A(1)*B(2)+A(1)*B(2))/S;
KE(2,2) = 3*(A(2)*B(2)+A(2)*B(2))/S;
KE(3,1) = -(A(1)*B(3)+A(1)*B(3))/S;
KE(3,2) = -(A(2)*B(3)+A(2)*B(3))/S;
KE(3,3) = 3*(A(3)*B(3)+A(3)*B(3))/S;
KE(4,1) = 0;
KE(4,2) = 4*(A(2)*A(3)+B(2)*B(3))/S;
KE(4,3) = KE(4,2);
KE(4,4) = 8*(A(2)*A(2)-A(1)*A(3)+B(2)*B(2)-B(1)*B(3))/S;
KE(5,1) = 4*(A(1)*A(3)+B(1)*B(3))/S;
KE(5,2) = 0;
KE(5,3) = KE(5,1);
KE(5,4) = 8*(A(1)*A(2)+B(1)*B(2))/S;
KE(5,5) = 8*(A(3)*A(3)-A(1)*A(2)+B(3)*B(3)-B(1)*B(2));
KE(6,1) = 4*(A(1)*A(2)+B(1)*B(2))/S;
KE(6,2) = KE(6,1);
KE(6,3) = 0;
KE(6,4) = 8*(A(1)*A(3)+B(1)*B(3))/S;
KE(6,5) = 8*(A(2)*A(3)+B(2)*B(3))/S;
KE(6,6) = 8*(A(1)*A(1)-A(2)*A(3)+B(1)*B(1)-B(2)*B(3))/S;

end

% ---------------------------------------------
% 程序说明
% 将总体系数矩阵以变带宽方式集成在一维向量GA中
% 变量说明：
% ND 节点总数
% XY 存放节点x、y坐标，[2×ND]
% NE 单元总数
% I6 存放单元节点编号，[6×NE]
% ID 存放第i行对角元Kii在GA向量中的位置j 表示ID(i)=j
% NP GA向量的长度
% GA 输出参数，存放总体系数矩阵
% ---------------------------------------------

function [GA] = UK2(NE,NP,I6,XY,ID)
GA = zeros(1,NP);
X = zeros(1,3);
Y = zeros(1,3);

for L = 1:NE
for J = 1:3
I = I6(J,L);
X(J) = XY(1,I);
Y(J) = XY(2,I);
end
KE = UKE2(X,Y);
for J = 1:6
NJ = I6(J,L);
for K = 1:J
NK = I6(K,L);
if(NJ<NK)
M = ID(NK)-NK+NJ;
else
M = ID(NJ)-NJ+NK;
end
GA(M) = GA(M)+KE(J,K);
end
end
end

end

% ------------------------------------------------------------------
% 程序说明
% 在 GA 和右侧列向量中代入第一类边界条件
% 变量说明：
% ND1 第一类边界条件节点数
% NB1 存放第一类边界条件的节点号，[1×ND1]
% U1 存放第一类边界节点的场值，[1×ND1]
% ID 存放第i行对角元Kii在GA向量中的位置j 表示ID(i)=j
% NP GA向量的长度
% GA 输入输出参数，存放总体系数矩阵，输出时，第一类边界条件已代入
% U 输出参数，输出加入第一类边界条件后 Ku = 0 的右侧列向量
% ------------------------------------------------------------------

function [GA,U]=UB12(ND1,NB1,U1,ND,ID,GA)
U = zeros(1,ND);
for I = 1:ND1
J = NB1(I);
K = ID(J);
GA(K) = GA(K)*10^10;
U(J) = GA(K)*U1(I);
end

end

% ------------------------------------------------------------------
% 程序说明
% 解对称变带宽线性方程组
% 方程组的系数已存放向量A中
% 线性方程组第i行对角元在A中位置j，事先由ID确定
% 变量说明：
% N 线性方程组的阶数
% NP 向量A的长度
% A 存在总体系数矩阵，[1,NP]
% ID 存放第i行对角元Kii在GA向量中的位置j 表示ID(i)=j
% B 输入输出参数，开始存放方程组的右侧向量，工作结束时，存放解向量
% IR 输出参数，标志，IR=0 时表示程序正常结束，否则表示系数矩阵奇异
% ------------------------------------------------------------------
function [B,IR] = BANDV(N,A,ID,B)
for I = 1:N;
I0 = ID(I)-I;
if(I~=1)
MI = ID(I-1)-I0+1;
for J = MI:I
J0 = ID(J)-J;
MJ = 1;
if(J>1)
MJ = ID(J-1)-J0+1;
end
MIJ = MI;
if(MJ>MI)
MIJ = MJ;
end
IJ = I0+J;
JM1 = J-1;
% if(MIJ>JM1)
% continue;
% end
for K = MIJ:JM1
IK = I0+K;
KK = ID(K);
JK = J0+K;
A(IJ) = A(IJ)-A(IK)*A(KK)*A(JK);
end
if(J~=I)
JJ = ID(J);
A(IJ) = A(IJ)/A(JJ);
B(I) = B(I)-A(IJ)*A(JJ)*B(J);
else
break;
end
end
end

II = I0+I;
if(A(II)==0)
IR = 1;return;
else
B(I) = B(I)/A(II);
end
end