function poi2D( ) % % -div(grad u) = f, in [0,1]x[0,1], % u = 0, on boundary % clear all, clc % triangulation g = [2 0 1 0 0 1 0; 2 1 1 0 1 1 0; 2 1 0 1 1 1 0; 2 0 0 1 0 1 0]'; [p,e,t] = initmesh(g,'hmax',0.05); figure(1); clf pdemesh(p,e,t) % assemble [A,b] = assemble(p,e,t,'f'); % solve U = A\b; % visualize figure(2); clf pdesurf(p,t,U) xlabel('x'), ylabel('y'), zlabel('U(x,y)') % subroutines ----------------------------------------------------------------- function z = f(x,y) z = 2*pi^2*sin(pi*x).*sin(pi*y); function [A,b] = assemble(p,e,t,f) Nt = size(t,2); Np = size(p,2); Ne = size(e,2); A = sparse(Np,Np); b = zeros(Np,1); for i = 1:Nt n = t(1:3,i); x = p(1,n); y = p(2,n); dx = [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; dy = [x(3)-x(2); x(1)-x(3); x(2)-x(1)]; area = 0.5*abs(x(2)*y(3)-y(2)*x(3)-x(1)*y(3)+y(1)*x(3)+x(1)*y(2)-y(1)*x(2)); A(n,n) = A(n,n) + (dx*dx'+dy*dy')/4/area; b(n) = b(n) + area/12*[2 1 1; 1 2 1; 1 1 2]*feval('f',x,y)'; end % BC for i = 1:Ne n = e(1,i); A(n,n) = 1e6; b(n) = 0; end