nad.m

利用matlab实现电磁场数值方法里面的有限差分与有限元法

function [A, B, C, num, bound] = NaD(gx, gy, NOE)  %NaD函数用来划分网格gx与gy分别为横纵方向的网格数,NOE为三角元总数
nx = gx + 1;
ny = gy + 1; %nx,ny分别为横纵方向的结点数
num(ny, nx) = 0; %用来存储节点号的数组
if(ny > nx)%用来存储边界节点的数组
    bound(4, ny) = 0; 
else
    bound(4, nx) = 0;
end
%将各个结点的编号存入num数组中
for j = 1 : nx
    for i = 1 : ny
        if (j == 1)
            num(i, j) = i; %第一列的结点编号
        else
            num(i, j) = num(i, j - 1) + ny; %非第一列的节点编号
        end
    end
end

%将边界结点的编号存入bound数组中，其中1，2，3，4行分别对应为上边界下边界左边界和右边界
for j = 1 : nx
    bound(1, j) = num(1, j); %上边界
    bound(2, j) = num(ny, j);%下边界    
end

for i = 2 : ny - 1
    bound(3, i) = num(i, 1); %左边界
    bound(4, i) = num(i, nx);%右边界
end

A(1) = 1; B(1) = 2; C(1) = ny + 2; A(2) = 1; B(2) = ny + 2; C(2) = ny + 1;
epc = 2 * gy;
%计算奇数号三角元的顶点编号
for e = 3 :2 : NOE
    for e1 = 3 : 2 :e
        if (mod(A(e1 - 2), ny) == ny - 1)
            A(e1) = A(e1 - 2) + 2;
            B(e1) = B(e1 - 2) + 2;
            C(e1) = C(e1 - 2) + 2;
        else
            A(e1) = A(e1 - 2) + 1;
            B(e1) = B(e1 - 2) + 1;
            C(e1) = C(e1 - 2) + 1;
        end
    end
end
%计算偶数号三角元的顶点编号
for e = 4 :2 : NOE
    for e1 = 4 : 2 :e
        if (mod(A(e1 - 2), ny) == ny - 1)
            A(e1) = A(e1 - 2) + 2;
            B(e1) = B(e1 - 2) + 2;
            C(e1) = C(e1 - 2) + 2;
        else
            A(e1) = A(e1 - 2) + 1;
            B(e1) = B(e1 - 2) + 1;
            C(e1) = C(e1 - 2) + 1;
        end
    end
end
end