fem.m

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

function v = FEM(a, gx, gy) %gx gy分别为横纵方向的网格数
NOE = 2 * gx * gy; %NOE为三角元总数
[A, B, C, num, bound] = NaD(gx, gy, NOE);%网格划分

nx = gx + 1; %横边上的节点数
ny = gy + 1; %纵边上的节点数
NON = nx * ny; %总节点数
%计算各点在矩阵中的坐标再添加到X和Y数组中
for n = 1 : nx * ny
    for j = 1 : nx
        for i = 1 : ny
            if (num(i, j) == n)
                X(n) = j;
                Y(n) = i; %输入节点坐标
            end
        end
    end
end
K = zeros(NON, NON);%系数矩阵
for e = 1 : NOE
    I = A(e);
    J = B(e);
    M = C(e);
    xi = (X(I) - 1) * (2 * a / gx);
    xj = (X(J) - 1) * (2 * a / gx);
    xm = (X(M) - 1) * (2 * a / gx);
    yi = - (Y(I) - 1) * (a / gy);
    yj = - (Y(J) - 1) * (a / gy);
    ym = - (Y(M) - 1) * (a / gy);
    bi = yj - ym;
    bj = ym - yi;
    bm = yi - yj;
    ci = xm - xj;
    cj = xi - xm;
    cm = xj - xi;
    s = 0.5 * (bi * cj - bj * ci);
    K(I, I) = K(I, I) + (bi * bi + ci * ci) / (4 * s);
    K(J, J) = K(J, J) + (bj * bj + cj * cj) / (4 * s);
    K(M, M) = K(M, M) + (bm * bm + cm * cm) / (4 * s);
    K(I, J) = K(I, J) + (bi * bj + ci * cj) / (4 * s);
    K(J, I) = K(I, J);
    K(I, M) = K(I, M) + (bi * bm + ci * cm) / (4 * s);
    K(M, I) = K(I, M);
    K(J, M) = K(J, M) + (bj * bm + cj * cm) / (4 * s);
    K(M, J) = K(J, M);
end

P = zeros(NON, 1);%方程右边的矩阵
%强制边界条件，修改系数矩阵和方程右边的矩阵
for n = 1 : nx
    ubn = bound(1, n);
    P = P - 10 * K(: , ubn);
end
%针对上下边界进行修改
for n = 1 : nx
    ubn = bound(1, n);
    dbn = bound(2, n);%ubn和dbn分别为上下边界点的编号
    P(ubn, 1) = 10;
    P(dbn, 1) = 0;
    K(ubn, :) = zeros(1, NON);
    K(:, ubn) = zeros(NON, 1);
    K(ubn, ubn) = 1;
    K(dbn, :) = zeros(1, NON);
    K(:, dbn) = zeros(NON, 1);
    K(dbn, dbn) = 1;
end
%针对左右边界进行修改
for n = 2 : ny - 1
    lbn = bound(3, n);
    rbn = bound(4, n);%lbn和rbn分别为左右边界点的编号
    P(lbn, 1) = 0;
    P(rbn, 1) = 0;
    K(lbn, :) = zeros(1, NON);
    K(:, lbn) = zeros(NON, 1);
    K(lbn, lbn) = 1;
    K(rbn, :) = zeros(1, NON);
    K(:, rbn) = zeros(NON, 1);
    K(rbn, rbn) = 1;
end
v1 = zeros(NON, 1);
v1 = inv(K) * P;
v = zeros(ny, nx);
m = 0;
for j = 1 : nx
    for i = 1 : ny
        m = m + 1;
        v(i, j) = v1(m, 1);
    end
end
contour(v, 30);
end