Contents

%%%%%%%%%%%%%%%%%%%%
%%%%%% FEM 2D %%%%%%
%%%%%%%%%%%%%%%%%%%%

% Solve: -laplace(u) = f    on unit square
% u = g on boundary
% *NEGATIVE LAPLACIAN**  plug in -f if you just want laplace(u) = f
clear all;

% PROCEDURE FOR MESHING SQUARE REGION AND ASSIGNING ELEMENTS
%_______________________
%| \        | \        |
%|   \   7  |   \  8   |
%|  5  \    |  6  \    |
%________\__|_______\_ |  % gridline = 1 in this case
%| \        | \        |    one gridline for both x and y components
%|   \   3  |   \   4  |
%| 1   \    |  2  \    |
%|_______\__|_______\__|

%%%%%% ERRORS (u = sin(x+y))%%%%%%%
% elements = 8 -- error = 0.5721  %
% elements = 98 -- error = 0.2000 %
% elements = 800 -- error =0.1004 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%% CONTROL CENTER %%%%%
ep2 = 10; % size of square
g = @(x,y) sin(x)+cos(y); % boundary constant
f = @(x,y) sin(x)+cos(y); % the function f
gridlines = 10; % decide how you'd like to divide up the square
resolution = 50; % resolution of plotted solution (e.g. 10 --> plot 10x10 grid)
%%%%%%%%%%%%%%%%%%%%%%%%%%

ep1 = 0;
elements = 2; % dont change this
for i = 1:gridlines
    elements = elements + 2*(2*i + 1); % make elements from grid #
end
elements
local = 3; % number of linear basis functions
Gnodes = (sqrt(elements/2)+1)^2; % total global nodes

% Local linear basis elements
N{1} = @(x,y) 1-x-y;
N{2} = @(x,y) x;
N{3} = @(x,y) y;

disp('ID, IEN, LM...')
% Assemble ID
root = sqrt(Gnodes);
ID = 1:Gnodes;
ID = reshape(ID,root,root);
ID = ID';
ID(1,:) = 0; ID(root,:) = 0; ID(:,1) = 0; ID(:,root) = 0;
ID = reshape(ID,1,Gnodes);
pp = size(find(ID>0),2);
C = 1:pp;
ID(find(ID>0)) = C;

% Record interior and boundary global function sets (use in assembling f)
I = find(ID>0); % Global interior function numbers
B = find(ID == 0); % Global boundary function numbers

% Assemble IEN
% first need to make a template matrix which has the global nodes in order
% just like how theyd look on paper.
% The IEN is assembled by looking through the template and recording which
% entries the local basis function is nonzero on.
template = Gnodes:-1:1;
template = reshape(template,root,root);
template = template';
if mod(root,2) == 1
    for i = 1:(root-1)/2
        temp = template(:,i);
        template(:,i) = template(:,root+1-i);
        template(:,root+1-i) = temp;
    end
else
    for i = 1:root/2
        temp = template(:,i);
        template(:,i) = template(:,(root+1-i));
        template(:,(root+1-i)) = temp;
    end
end

% First row numbers for IEN
row1 = [];
for i = root:-1:1
    row1 = [row1 template(i,2:end)];
    row1 = [row1 template(i,1:end-1)];
end
row1(1:root-1) = [];
row1(end:-1:end-root+2)=[];
% second row
row2 = [];
for i = root:-1:1
   row2 = [row2 template(i,1:end-1)];
   row2 = [row2 template(i,2:end)];
end
row2(1:root-1) = [];
row2(end:-1:end-root+2)=[];
% third row
row3 = [];
for i = 1:root-1
   row3 = [row3 template(root-i,1:root-1)];
   row3 = [row3 template(root-i+1,2:end)];
end

% IEN
IEN = [row1; row2; row3;];

% LM
LM = ID(IEN);

% Integral involved here:
% integral (Ni_x*Nj_x + Ni_y*Nj_y)dydx (0<x<1, 0<y<1-x)
% quad2d does not work for constant functions for some strange reason, so i
% will just hardcode in all the integrals...
k1 = zeros(local,local);
k1(1,1)=1;k1(1,2)=-.5;k1(1,3)=-.5;k1(2,1) = -.5;k1(2,2)=.5;
k1(2,3)=0;k1(3,1)=-.5;k1(3,2)=0;k1(3,3)=.5;

% % Compile translation functions
% set of lower triangles (element numbers)
lower = [];
for i = 1:2*(gridlines+1):elements-gridlines-2
    lower = [lower i:i+gridlines];
end
% % set of upper triangles (translation functions different for each
% % triangle)
upper = 1:elements;
upper(lower) = [];

width = (ep2-ep1)/sqrt(elements/2); %width of each triangle
% x = horizontal global coordinate, y = vertical global coordinate
% z = horizontal local coordinate, w = vertical local coordinate
% x,y,z,w translations on lower triangle elements
disp('Translation functions...')
j = 1; p = 1;
for e = lower
    xofz{e} = @(z) width*z + (j-1)*width; % x(z) local to global
    zofx{e} = @(z) z/width - (j-1); % z(x) global to local
    yofw{e} = @(z) width*z + (p-1)*width; % y(w) local to global
    wofy{e} = @(z) z/width - (p-1); % w(y) global to local
    % counters to keep track of constants
    j = j+1;
    if j == gridlines + 2
        j = 1;
        p = p+1;
    end
end
% x,y,z,w translations on upper triangle elements
j = 1; p = 1;
for e = upper
    xofz{e} = @(z) width*(1-z) + (j-1)*width;
    zofx{e} = @(z) 1-z/width + (j-1);
    yofw{e} = @(z) width*(1-z) + (p-1)*width;
    wofy{e} = @(z) 1-z/width + (p-1);
    j = j+1;
    if j == gridlines + 2
        j = 1;
        p = p+1;
    end
end

% Multiplier for global to local integral
% Here, the jacobian for (f,N) integral is x'(z)*y'(w). Both x' and y' are
% equal to width, but change signs on some elements. Since they both share
% the same sign in each case, the jacobian is always just width^2
% (pretty convenient)
jacobian = width*width;

%%%%% SERENDIPITY %%%%%%
% It turns out that, through 1D linear change of variables, chain rule, and
% the jacobian, the global integral of delNi*delNj is EXACTLY the same as
% the local integral of the local functions. No multipliers or anything
% needed - local integration is all that's required (for K at least)

%%%%%% First: Boundary function g needs to be converted into aN1+bN2+cN3 ( a linear
% combination of local basis functions). Here we find the coefficients
% corresponding to each element and store them in a matrix g(i,e) to be
% used in the assembly of local f
% Get all the necessary points in order:
xvec = linspace(ep1,ep2, gridlines+2);
yvec = linspace(ep1,ep2, gridlines+2);

% x,y coordinates of A at each element
% (A is the corner of each triangle)
disp('Coefficients of boundary function...')
newtemplate = 1:root;
temp1 = [];
temp2 = [];
for i = 1:gridlines+1
     temp1 = [temp1 newtemplate(1:end-1)];
     temp1 = [temp1 newtemplate(2:end)];
     temp2 = [temp2 newtemplate(i)*ones(1,gridlines+1) ...
          newtemplate(i+1)*ones(1,gridlines+1)];
end
A1(:,1) = xvec(temp1)';
A1(:,2) = yvec(temp2)';
% Start again but assemble coordinates of B (the right vertex of each
% triangular element)
temp1 = [];
temp2 = [];
for i = 1:gridlines+1
     temp1 = [temp1 newtemplate(2:end)];
     temp1 = [temp1 newtemplate(1:end-1)];
     temp2 = [temp2 newtemplate(i)*ones(1,gridlines+1) ...
          newtemplate(i+1)*ones(1,gridlines+1)];
end
B1(:,1) = xvec(temp1)';
B1(:,2) = yvec(temp2)';
% Start again but assemble coordinates of C (the upper vertex of each
% triangular element)
temp1 = [];
temp2 = [];
for i = 1:gridlines+1
     temp1 = [temp1 newtemplate(1:end-1)];
     temp1 = [temp1 newtemplate(2:end)];
     temp2 = [temp2 newtemplate(i+1)*ones(1,gridlines+1) ...
          newtemplate(i)*ones(1,gridlines+1)];
end
C1(:,1) = xvec(temp1)';
C1(:,2) = yvec(temp2)';

Nmat = [1 0 0; -1 1 0 ; -1 0 1]; % Matrix of coefficients of local basis functions
gtild = zeros(local,elements); % This will be the matrix containing the essential coefficients
% These coefficients describe a local version of the boundary function on
% each element.

for e = 1:elements
    gatA = g(A1(e,1),A1(e,2)); % g evaluated at corner of element
    gatB = g(B1(e,1),B1(e,2)); % g evaluated at right vertex
    gatC = g(C1(e,1),C1(e,2)); % g evaluated at upper vertex

    vecb = [gatA;gatB-gatA;gatC-gatA]; % These are the coefficients of a function
    % that interpolates g at the vertices of each element
    % i.e. g(interpolated on element e) = gatA + (gatB-gatA)*x + (gatC-gatA)*y

    gtild(:,e) = Nmat\vecb;
    % this is solving the matrix equaton Nmat*gtild = vecb. This
    % provides the needed coefficients that describe the
    % linear combination of local basis functions that match the boundary
    % function g at the vertices of each element.
    % e.g. g(on element 1) = gtild(1,1)*N1 + gtild(2,1)*N2 + gtild(3,1)*N3
end

disp('Local f...')
% Assemble local f
wmax = @(z) 1-z; % hypotenuse of local triangle - use for local integral
fe = zeros(local,elements); % here im saying fe is a matrix
% written convention f_i^e ---> program: fe(i,e) i = local basis #, e = element #
for e = 1:elements
    % put 'e' here if youd like to track progress
    floc = @(z,w) f(xofz{e}(z),yofw{e}(w)); % local f function for integration
    for i = 1:local
        integrand = @(z,w) floc(z,w).*N{i}(z,w);
        % adding in (f,N)
        fe(i,e) = fe(i,e) + jacobian*quad2d(integrand,0,1,0,wmax);
        %%%% ADD BILINEAR BOUNDARY/INTERIOR PART TO LOCAL F %%%%%
        % Looking at the global weak form of this whole problem, it is
        % apparent that the only bilinear parts added to local f are ones
        % that consist of a(N_i,N_j) where i and j strictly belong to
        % DIFFERENT sets (e.g. i in Interior nodes, j in Boundary nodes - we
        % check this by referencing IEN and the sets I and B).
        % The following code emulates this process.
        j = [1 2 3];
        j(i) = [];% other entries of IEN column e BESIDES row i to be checked
        a = IEN(i,e); % the entry in question
        b = IEN(j(1),e); % other entry of IEN column e
        c = IEN(j(2),e); % other entry of IEN column e
        ainI = ismember(a,I); % is IEN(i,e) a member of I?
        if ainI == 1 % IEN(i,e) is a member of I
            binB = ismember(b,B);
            cinB = ismember(c,B); %check if other entries are in B
            if binB == 1
                fe(i,e) = fe(i,e) - gtild(i,e)*k1(i,j(1)); % add integral to f if IEN row
                % i is in I and row j(1) is in B
            end
            if cinB == 1
               fe(i,e) = fe(i,e) - gtild(i,e)*k1(i,j(2)); % add integral to f if IEN row
                % i is in I and row j(2) is in B
            end
        else % this is assuming IEN(i,e) is in B
            binI = ismember(b,I);
            cinI = ismember(c,I); %check if other entries are in I
            if binI == 1
                fe(i,e) = fe(i,e) - gtild(i,e)*k1(i,j(1)); % add integral to f if IEN row
                % i is in B and row j(1) is in I
            end
            if cinI == 1
               fe(i,e) = fe(i,e) - gtild(i,e)*k1(i,j(2)); % add integral to f if IEN row
                % i is in B and row j(2) is in I
            end
        end
    end
end

disp('Global K and F...')
% Assemble Global K and F
o = length(find(ID>0)); % number of nonzero entries of ID matrix - # of rows
% of K and F
K = zeros(o,o);
F = zeros(o,1);

% Master loop
for e = 1:elements
    for i = 1:local
        if LM(i,e) ~= 0
           F(LM(i,e)) = F(LM(i,e)) + fe(i,e);
           for j = 1:local
               if LM(j,e) ~= 0
                   K(LM(i,e),LM(j,e)) = K(LM(i,e),LM(j,e)) + k1(i,j);
               end
           end
        end
    end
end

disp('Solving Kd = F...')
% Solution to Kd = F
d = K\F;

% Assemble true coefficients of u
d_new = zeros(Gnodes,1);
points = combvec(xvec,yvec);
for i = B
   d_new(i) = g(points(1,i),points(2,i)); % pretty simple - boundary coefficients of
   % the solution u should just be g evaluated at that specific point
end
d_new(I) = d; % add in calculated interior coefficients

disp('Creating FEM solution points...')

Create matrix of solution points to plot

horizontal and vertical mesh to evaluate and plot solution

xmesh = linspace(ep1,ep2,resolution);
ymesh = linspace(ep1,ep2,resolution);
FEM = zeros(resolution);
for y = 1:resolution
    for x = 1:resolution
        xpoint = xmesh(x);
        ypoint = ymesh(y);  % the points to be evaluated at
        E = TheElement(gridlines,xpoint,ypoint,ep1,ep2); % corresponding element
        for i = 1:local
            globalnum = IEN(i,E); % index for d_new coefficient (global function #)
            FEM(x,y) = FEM(x,y) ...
                + d_new(globalnum)*N{i}(zofx{E}(xpoint),wofy{E}(ypoint));
        end
    end
end

disp('Plotting...')
[X, Y] = meshgrid(xmesh,ymesh);
EXACT = sin(X)+cos(Y);
EXACT = reshape(EXACT, resolution, resolution);
figure(1)
surf(X,Y, EXACT); rotate3d on;
title('Exact Solution');

figure(2)
FEM = FEM';
surf(X,Y, FEM); rotate3d on;
title('FEM Solution');

FEM = reshape(FEM, resolution^2,1);
EXACT = reshape(EXACT, resolution^2,1);
error = norm(FEM-EXACT, 'inf')

elements =

   242

ID, IEN, LM...
Translation functions...
Coefficients of boundary function...
Local f...
Global K and F...
Solving Kd = F...
Creating FEM solution points...
Plotting...

error =

    0.6933

Published with MATLAB® 7.12