I found an example using triangular elements here. I then went on to modify the way it generates the mesh to replace the triangular elements with rectangular ones, but I don't know how to integrate over them. Here is my version of it:
%3.6 femcode.m
% [p,t,b] = squaregrid(m,n) % create grid of N=mn nodes to be listed in p
% generate mesh of T=2(m-1)(n-1) right triangles in unit square
m=6; n=5; % includes boundary nodes, mesh spacing 1/(m-1) and 1/(n-1)
[x,y]=ndgrid((0:m-1)/(m-1),(0:n-1)/(n-1)); % matlab forms x and y lists
p=[x(:),y(:)]; % N by 2 matrix listing x,y coordinates of all N=mn nodes
nelems=(m-1)*(n-1);
t=zeros(nelems,3);
for e=1:nelems
t(e) = e + floor(e/6);
t(e, 2) = e + floor(e/6) + 1;
t(e, 3) = e + floor(e/6) + 6;
t(e, 4) = e + floor(e/6) + 7;
end
% final t lists 4 node numbers of all rectangles in T by 4 matrix
b=[1:m,m+1:m:m*n,2*m:m:m*n,m*n-m+2:m*n-1]; % bottom, left, right, top
% b = numbers of all 2m+2n **boundary nodes** preparing for U(b)=0
% [K,F] = assemble(p,t) % K and F for any mesh of rectangles: linear phi's
N=size(p,1);T=nelems; % number of nodes, number of rectangles
% p lists x,y coordinates of N nodes, t lists rectangles by 4 node numbers
K=sparse(N,N); % zero matrix in sparse format: zeros(N) would be "dense"
F=zeros(N,1); % load vector F to hold integrals of phi's times load f(x,y)
for e=1:T % integration over one rectangular element at a time
nodes=t(e,:); % row of t = node numbers of the 3 corners of triangle e
Pe=[ones(4,1),p(nodes,:),p(nodes,1).*p(nodes,2)]; % 4 by 4 matrix with rows=[1 x y xy]
Area=abs(det(Pe)); % area of triangle e = half of parallelogram area
C=inv(Pe); % columns of C are coeffs in a+bx+cy to give phi=1,0,0 at nodes
% now compute 3 by 3 Ke and 3 by 1 Fe for element e
grad=C(2:3,:);Ke=Area*grad'*grad; % element matrix from slopes b,c in grad
Fe=Area/3; % integral of phi over triangle is volume of pyramid: f(x,y)=1
% multiply Fe by f at centroid for load f(x,y): one-point quadrature!
% centroid would be mean(p(nodes,:)) = average of 3 node coordinates
K(nodes,nodes)=K(nodes,nodes)+Ke; % add Ke to 9 entries of global K
F(nodes)=F(nodes)+Fe; % add Fe to 3 components of load vector F
end % all T element matrices and vectors now assembled into K and F
% [Kb,Fb] = dirichlet(K,F,b) % assembled K was singular! K*ones(N,1)=0
% Implement Dirichlet boundary conditions U(b)=0 at nodes in list b
K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F
K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K
Kb=K; Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb
% Solving for the vector U will produce U(b)=0 at boundary nodes
U=Kb\Fb % The FEM approximation is U_1 phi_1 + ... + U_N phi_N
% Plot the FEM approximation U(x,y) with values U_1 to U_N at the nodes
% surf(p(:,1),p(:,2),0*p(:,1),U,'edgecolor','k','facecolor','interp');
% view(2),axis equal,colorbar
I started to edit the section of the code which used to integrate over triangular elements but I'm not sure how to proceed or if it even is done in a similar way for rectangular elements.
UPDATE
So I've tried to incorporate what Dohyun suggested with my limited understanding and here's what I've got now:
m=12; n=12; % includes boundary nodes, mesh spacing 1/(m-1) and 1/(n-1)
[x,y]=ndgrid((0:m-1)/(m-1),(0:n-1)/(n-1)); % matlab forms x and y lists
p=[x(:),y(:)]; % N by 2 matrix listing x,y coordinates of all N=mn nodes
nelems=(m-1)*(n-1);
t=zeros(nelems,4);
a=0;
for e=1:nelems
t(e) = e + a;
t(e, 2) = e + a + 1;
t(e, 3) = e + a + m + 1;
t(e, 4) = e + a + m;
a = floor(e/n);
end
% final t lists 4 node numbers of all rectangles in T by 4 matrix
b=[1:m,m+1:m:m*n,2*m:m:m*n,m*n-m+2:m*n-1]; % bottom, left, right, top
% b = numbers of all 2m+2n **boundary nodes** preparing for U(b)=0
N=size(p,1);T=nelems; % number of nodes, number of rectangles
% p lists x,y coordinates of N nodes, t lists rectangles by 4 node numbers
K=sparse(N,N); % zero matrix in sparse format: zeros(N) would be "dense"
F=zeros(N,1); % load vector F to hold integrals of phi's times load f(x,y)
for e=1:T % integration over one rectangular element at a time
nodes=t(e,:); % row of t = node numbers of the 3 corners of triangle e
Pe=[ones(4,1),p(nodes,:),p(nodes,1).*p(nodes,2)]; % 4 by 4 matrix with rows=[1 x y xy]
Area=abs(det(Pe(1:3,1:3))); % area of triangle e = half of parallelogram area
C=inv(Pe); % columns of C are coeffs in a+bx+cy to give phi=1,0,0 at nodes
% now compute 3 by 3 Ke and 3 by 1 Fe for element e
% grad=C(2:3,:);
% constantKe=Area*grad'*grad; % element matrix from slopes b,c in grad
for i=1:4
for j=1:4
syms x y
Kn = int(int( ...
C(i,2)*C(j,2)+ ...
(C(i,2)*C(j,4)+C(i,4)*C(j,2))*y + ...
C(i,4)*C(j,4)*y^2 + ...
C(i,3)*C(j,3) + ...
(C(i,4)*C(j,3)+C(i,3)*C(j,4))*x + ...
C(i,4)*C(j,4)*x^2 ...
, x, Pe(1, 2), Pe(2, 2)), y, Pe(1, 3), Pe(3, 3));
K(nodes(i),nodes(j)) = K(nodes(i),nodes(j)) + Kn;
end
end
Fe=Area / 3; % integral of phi over triangle is volume of pyramid: f(x,y)=1
% multiply Fe by f at centroid for load f(x,y): one-point quadrature!
% centroid would be mean(p(nodes,:)) = average of 3 node coordinates
% K(nodes,nodes)=K(nodes,nodes)+Ke; % add Ke to 9 entries of global K
F(nodes)=F(nodes)+Fe; % add Fe to 4 components of load vector F
end % all T element matrices and vectors now assembled into K and F
% [Kb,Fb] = dirichlet(K,F,b) % assembled K was singular! K*ones(N,1)=0
% Implement Dirichlet boundary conditions U(b)=0 at nodes in list b
K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F
K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K
Kb=K; Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb
% Solving for the vector U will produce U(b)=0 at boundary nodes
U=Kb\Fb; % The FEM approximation is U_1 phi_1 + ... + U_N phi_N
U2=reshape(U',m,n);
% Plot the FEM approximation U(x,y) with values U_1 to U_N at the nodes
surf(U2)
I think I can use definite integrals instead of numerical ones, but the result doesn't match that of the original program.
For P1 triangular element case, gradient happens to be constant. Hence integration is simply
area*grad'*grad.However, for bilinear case, inner product of gradient is 2nd order polynomial. Hence you need to use numerical integration.
So, inside of simple multiplication, you need another loop which computes basis value at quadrature points.
Also, your formula for area is wrong. Search for the affine mapping.
I have repository on github which implements poisson equation from 1D to 3D with arbitrary order polynomial. If you are interested, come visit https://github.com/dohyun64/fem_dohyun