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

正在学习有限单元法，大家有需要的可以看看

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

在线求地球物理的大佬指点

% ---------------------------------------------
% 主函数
% 第一类边界条件、线性插值，三角单元剖分
% 编写时间：2022-08-30
% ---------------------------------------------
clear all;clc;

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

% 单元节点编号
fid1=fopen('单元节点编号1.dat','rt');
line1=fgets(fid1);
data1 = fscanf(fid1,'%e %e %e %e',[4 inf]);
fclose(fid1);
data1 = data1';
NEi = data1(:,1)';          % 单元序号
i = data1(:,2)';
j = data1(:,3)';
m = data1(:,4)';
I3 = [i;j;m];               % 单元节点编号 [3,NE]

% 节点坐标
fid2 = fopen('节点坐标1.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]

% 第一类边界节点号和场值
fid3 = fopen('第一类边界节点号和场值1.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]

%%

IW = MBW(NE,I3);                            % 计算半带宽
SK = UK1(ND,NE,IW,I3,XY);                   % 形成总体系数矩阵
[SK,U] = UB1(ND1,NB1,U1,ND,IW,SK);          % 代入第一类边界条件
[P,IE] = LDLT(SK,ND,IW,U);                  % 解线性方程组
 

% ---------------------------------------------
% 程序说明
% 功能：三角单元、线性插值时计算单元系数矩阵 ke
% 变量说明： 
% X 存放节点 x 坐标，存放顺序 xi xj xm
% Y 存放节点 y 坐标，存放顺序 yi yj ym
% KE 输出参数，存放单元系数矩阵
% ---------------------------------------------

function [KE] = UKE1(X,Y)
KE = zeros(3,3);
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 = 2*(A(1)*B(2)-A(2)*B(1));
for i = 1:3
    for j = 1:i
        KE(i,j) = (A(i)*A(j)+B(i)*B(j))/S;
    end
end

end

% ------------------------------------------------------
% 程序说明
% 三角单元、线性插值时用定带宽储存的方法集成总体矩阵
% 变量说明
% ND 节点总数
% NE 单元总数
% IW 半带宽
% I3 存放单元节点编号，[3×NE]
% XY 存放节点的xy坐标，存放顺序：x1，y1，……，xND，yND;[2×ND]
% SK 输出参数，存放定带宽存储的总体系数矩阵，[ND×IW]
% ------------------------------------------------------

function [SK] = UK1(ND,NE,IW,I3,XY)
SK = zeros(ND,IW);
for L=1:NE
    X = zeros(1,3);
    Y = zeros(1,3);
    
    for J1 = 1:3
       I = I3(J1,L);
       X(J1) = XY(1,I);
       Y(J1) = XY(2,I);           
    end
    
    KE = UKE1(X,Y);
    
    for J2 = 1:3
        NJ = I3(J2,L);
        for K = 1:J2
            NK = I3(K,L);
            if(NJ<NK)     
                NJ = NJ-NK+IW;
                SK(NK,NJ) = SK(NK,NJ)+KE(J2,K);
                NJ = NJ+NK-IW;
            else
                NK = NK-NJ+IW;
                SK(NJ,NK) = SK(NJ,NK)+KE(J2,K);
            end           
        end
    end
    
end

end

% ------------------------------------------------------------------------
% 程序说明
% 在定带宽储存的总体系数矩阵和右侧列向量上加上第一类边界条件
% 变量说明：
% ND1 第一类边界条件的节点数
% NB1 存放第一类边界条件的节点号
% U1 存放第一类边界点的场值
% IW 半带宽
% SK 输出参数，定带宽存放总体系数矩阵，输出时，第一类边界条件已代入 [ND×IW]
% U  输出参数，输出加入第一类边界条件后的 Ku=0 式的右端列向量 [ND×1]
% ------------------------------------------------------------------------

function [SK,U] = UB1(ND1,NB1,U1,ND,IW,SK)
U = zeros(ND,1);

for i = 1:ND1
    j = NB1(i);
    SK(j,IW) = SK(j,IW)*10^10;
    U(j) = SK(j,IW)*U1(i);
end

end
 

% ---------------------------------------------
% 程序说明
% 计算总体系数矩阵的半带宽
% 变量说明：
% NE 单元总数
% I3 存放单元节点编号，[3×NE]
% IW 输出参数，半带宽
% ---------------------------------------------

function [IW] = MBW(NE,I3)
IW = 0;
for i = 1:NE
    M = max([abs(I3(1,i)-I3(2,i)),abs(I3(2,i)-I3(3,i)),abs(I3(3,i)-I3(1,i))]);
    if M+1>IW
        IW = M+1;
    end
end

end

% ------------------------------------------------------------------------
% 程序说明
% 解定带宽储存的对称带型线性方程组
% 参数说明：
% A 存放对称带型方程组的系数矩阵的下三角部分,[N×IW]
% N 方程组的阶数
% IW 对称带型方程组的半带宽
% P 输入输出参数，[N×1]
%   开始存放方程组的右侧列向量，工作结束时，存放解向量
% IE 输出参数，标志，IE=0 表示程序正常结束，IE=1 表示系数矩阵奇异
% ------------------------------------------------------------------------

function [P,IE] = LDLT(A,N,IW,P)
for I = 1:N
    
    if(I<=IW)
        IT = 1;
    else
        IT = I-IW+1;
    end
    
    K = I-1;
    
    if(I==1)
        if(A(I,IW)==0)
            IE = 1;
            return;
        end
        continue;
    end
    
    for L = IT:K
        IL = L+IW-I;
        B = A(I,IL);
        A(I,IL) = B/A(L,IW);
        P(I) = P(I)-A(I,IL)*P(L);
        MI = L+1;
        for J = MI:I
            IJ = J+IW-I;
            JL = L+IW-J;
            A(I,IJ) = A(I,IJ)-A(J,JL)*B;
        end
    end    
    if(A(I,IW)==0)
        IE = 1;
        return;
    end
end

for J = 1:N
    if(J<=IW)
        IT=N;
    else
        IT = N-J+IW;
    end
    
    I = N-J+1;
    P(I) = P(I)/A(I,IW);
    
    if(J ~= 1)
        K = I+1;
        for MJ = K:IT
            IJ = I-MJ+IW;
            P(I) = P(I)-P(MJ)*A(MJ,IJ);
        end
    end
    continue
end
IE = 0;

end