fem2dt2.m

使用有限元方法, 2D FEM计算热扩散方程. 三角形分割.

clear all;
clc;

n = 4;
h = 1./ n;
nsq = n*n;
A = zeros(nsq,nsq);  
B = zeros(nsq,1);

for j=1:1:n                                 % Outer loop over y
    yj_1 = (j-1)*h;
    yj = j*h;
    for i=1:1:n                             % Inner loop over x
        xi_1 = (i-1)*h;
        xi = i*h;
        
        
        % Handle two triangles within each square
        iUpperRight = (j-1)*n + i;          % Indices of nodes in
        iUpperLeft  = iUpperRight - 1;      % each triangle
        iLowerRight = iUpperRight - n;
        iLowerLeft  = iLowerRight - 1; 
        area = 0.5*h^2;                      % Area of each triangle
        
        % Handle two triangles within each square
        if i > 1 & j > 1           % Lower left node is a free node
            k = iLowerLeft;
            A(k,k) = A(k,k) + area*(1/h)^2;   % Contribution from right triangle
            A(k,k+1) = A(k,k+1) - area*(1/h)^2; 
            A(k,k) = A(k,k) + area*(1/h)^2;   % Contribution from left triangle
            A(k,k+n) = A(k,k+n) - area*(1/h)^2;
        end
        if j > 1                   % Lower right node is a free node
            k = iLowerRight;
            A(k,k) = A(k,k) + area*((1/h)^2 + (1/h)^2);
            ti = area*((1/h)^2 + (1/h)^2);
            if i > 1, 
                A(k,k-1) = A(k,k-1) - area*(1/h)^2; 
            end;
            A(k,k+n) = A(k,k+n) - area*(1/h)^2;
        end
        if i > 1                   % Upper left node is a free node
            k = iUpperLeft;
            A(k,k) = A(k,k) + area*((1/h)^2 + (1/h)^2);
            if j > 1, 
                A(k,k-n) = A(k,k-n) - area*(1/h)^2; 
            end;
            A(k,k+1) = A(k,k+1) - area*(1/h)^2;
        end
        
        % Upper right node is always a free node
        k = iUpperRight;
        A(k,k) = A(k,k) + area*(1/h)^2;      % Contribution from right triangle
        if j > 1
            A(k,k-n) = A(k,k-n) - area*(1/h)^2; 
        end        
        A(k,k) = A(k,k) + area*(1/h)^2;      % Contribution from left triangle
        if i > 1 
            A(k,k-1) = A(k,k-1) - area*(1/h)^2; 
        end 
        
        % Right-hand side vector.  Use second order quadrature formula
        % involving function values at midpoints of each side.     
        MidPt_R1 = [0.5*(xi+xi_1), yj_1];      % Midpoints of sides of right triangle
        MidPt_R2 = [xi, 0.5*(yj+yj_1)];
        MidPt_R3 = [0.5*(xi+xi_1), 0.5*(yj+yj_1)];       
      
        MidPt_L1 = [xi_1, 0.5*(yj+yj_1)];      % Midpoints of sides of left triangle
        MidPt_L2 = [0.5*(xi+xi_1), yj];                    
        MidPt_L3 = MidPt_R3;
      
        % Evaluate f at midpoints.
        xt = MidPt_R1(1);     yt = MidPt_R1(2);
        f_R1 = xt*(1-xt) + yt*(2-yt);
      
        xt = MidPt_R2(1);     yt = MidPt_R2(2);
        f_R2 = xt*(1-xt) + yt*(2-yt);
      
        xt = MidPt_R3(1);     yt = MidPt_R3(2);
        f_R3 = xt*(1-xt) + yt*(2-yt);
      
        xt = MidPt_L1(1);     yt = MidPt_L1(2);
        f_L1 = xt*(1-xt) + yt*(2-yt);
      
        xt = MidPt_L2(1);     yt = MidPt_L2(2);
        f_L2 = xt*(1-xt) + yt*(2-yt);
      
        xt = MidPt_L3(1);     yt = MidPt_L3(2);
        f_L3 = xt*(1-xt) + yt*(2-yt);
            
        % detlaX, deltaY: Translate triangle to have lower left corner at (0,0) 
        MidPt_R1 = MidPt_R1 - [xi_1, yj_1];
        MidPt_R2 = MidPt_R2 - [xi_1, yj_1];
        MidPt_R3 = MidPt_R3 - [xi_1, yj_1];
      
        MidPt_L1 = MidPt_L1 - [xi_1, yj_1];
        MidPt_L2 = MidPt_L2 - [xi_1, yj_1];
        MidPt_L3 = MidPt_L3 - [xi_1, yj_1];
        
        if i >=2  & j >= 2                          % Lower left node is a free node
            k = iLowerLeft;
            Phi_R1 = 1 - MidPt_R1(1) / h;           % 1 - dX/h, Contribution from right triangle
            Phi_R2 = 1 - MidPt_R2(1) / h;
            Phi_R3 = 1 - MidPt_R3(1) / h;
            B(k) = B(k) + (area/3)*(f_R1*Phi_R1 + f_R2*Phi_R2 + f_R3*Phi_R3);       %??
      
            Phi_L1 = 1 - MidPt_L1(2) / h;           % Contribution from left triangle
            Phi_L2 = 1 - MidPt_L2(2) / h;
            Phi_L3 = 1 - MidPt_L3(2) / h;
            B(k) = B(k) + (area/3)*(f_L1*Phi_L1 + f_L2*Phi_L2 + f_L3*Phi_L3);
        end
        
        if j >= 2                   % Lower right node is a free node
            k = iLowerRight;
            Phi_R1 = (1/h)*MidPt_R1(1) - (1/h)*MidPt_R1(2);     % Phi_R = (dXr- dYr) / h
            Phi_R2 = (1/h)*MidPt_R2(1) - (1/h)*MidPt_R2(2);
            Phi_R3 = (1/h)*MidPt_R3(1) - (1/h)*MidPt_R3(2);
           
            B(k) = B(k) + (area/3)*(f_R1*Phi_R1 + f_R2*Phi_R2 + f_R3*Phi_R3);
        end
        
        if i >= 2                   % Upper left node is a free node
            k = iUpperLeft;
            Phi_L1 = -(1/h)*MidPt_L1(1) + (1/h)*MidPt_L1(2);    % Phi_L =(dYl - dXl) / h
            Phi_L2 = -(1/h)*MidPt_L2(1) + (1/h)*MidPt_L2(2);
            Phi_L3 = -(1/h)*MidPt_L3(1) + (1/h)*MidPt_L3(2);
            B(k) = B(k) + (area/3)*(f_L1*Phi_L1 + f_L2*Phi_L2 + f_L3*Phi_L3);
        end
    
        % Upper right node is always a free node, 
        k = iUpperRight;                                  
        Phi_R1 = (1/h)*MidPt_R1(2);                              % Phi_R = dYr / h, Contribution from right triangle 
        Phi_R2 = (1/h)*MidPt_R2(2);
        Phi_R3 = (1/h)*MidPt_R3(2);
        B(k) = B(k) + (area/3)*(f_R1*Phi_R1 + f_R2*Phi_R2 + f_R3*Phi_R3);
        
        Phi_L1 = (1/h)*MidPt_L1(1);            % Contribution from left triangle
        Phi_L2 = (1/h)*MidPt_L2(1);
        Phi_L3 = (1/h)*MidPt_L3(1);
        B(k) = B(k) + (area/3)*(f_L1*Phi_L1 + f_L2*Phi_L2 + f_L3*Phi_L3);

        % Extra boundary term.  Use 2-point Gauss quadrature here.
        % du(1, y)/dx = y*(1-0.5*y)
        if i==n
            Pt1 = 0.5*(yj+yj_1) - h/sqrt(3);    % Gauss quadrature points
            Pt2 = 0.5*(yj+yj_1) + h/sqrt(3);
            
            usubn1 = -Pt1 * (1 - 0.5*Pt1);          % u_n at quadrature points
            usubn2 = -Pt2 * (1 - 0.5*Pt2);
            Pt1 = Pt1 - yj_1;                   % Translate interval
            Pt2 = Pt2 - yj_1;

            k = iUpperRight;                % Upper right basis function
            Phi1 = (1/h)*Pt1;
            Phi2 = (1/h)*Pt2;
            B(k) = B(k) + 0.5*h*(usubn1*Phi1 + usubn2*Phi2);
            
            if j >= 2              
                k = iLowerRight;            % Lower right basis function
                Phi1 = 1 - (1/h)*Pt1;
                Phi2 = 1 - (1/h)*Pt2;
                B(k) = B(k) + 0.5*h*(usubn1*Phi1 + usubn2*Phi2);
            end
        end     
    end
end

%  Solve linear system.
u = A \ B;

u2d = reshape(u, n, n);   % B = reshape(A,m,n) returns the m-by-n matrix B whose elements are taken column-wise from A
surf(u2d);
xlabel('x');
ylabel('y');
zlabel('u_FEM');

u_true = zeros(n, n);
for j=1:1:n
    yj = j*h;
    for i=1:1:n
        xi = i*h;
        u_true(i,j) = xi*(1-xi)*yj*(1-.5*yj);
  end
end

figure(2);
surf(u_true);
xlabel('x');
ylabel('y');
zlabel('u_True');