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