%% Clean up
clear
close all

%% Load data
% Load mesh
fid=fopen('gitter.grid','r');
num=fscanf(fid,'%i',2);
np=num(1);
nt=num(2);
ptmp=fscanf(fid,'%f',[2,np]);
ttmp=fscanf(fid,'%f',[3,nt]);
fclose(fid);
% Change for Matlab
px=ptmp(1,:)';
py=ptmp(2,:)';
triangles=ttmp'+1;
clear ptmp ttmp
% Load right-hand side
load f.dat

%  np = 4;
%  nt = 2;
%  px = [0;1;0;1];
%  py = [0;0;1;1];
%  triangles = [1,2,3;2,3,4];
%  f = [1,2,3,4]';

%% Plot
% Plot mesh
subplot(1,2,1);
triplot(triangles,px,py)
axis square
% Plot height data
subplot(1,2,2);
trisurf(triangles,px,py,f)
colormap winter;
shading interp;
lighting phong;
light;
view ( -40, 80 );
axis square

%% mass and stiffness matrix assembly
beta = 1e5;
I = zeros( 9 * nt, 1 );
J = I;
VM = I;
VL = I;
MT = [2,1,1;1,2,1;1,1,2]/12;
counter = 1;
for t = 1:nt
    % compute triangle area
    points = [px(triangles(t,:)),py(triangles(t,:))]';
    sides = [points(:,3)-points(:,2), points(:,1)-points(:,3), points(:,2)-points(:,1)];
    T = abs(det ( [ sides(:,1), sides(:,2) ] )/2);
    % local assembly
    for nodei = 1:3
        for nodej = 1:3
            I(counter) = triangles(t,nodei);
            J(counter) = triangles(t,nodej);
            VM(counter) = MT(nodei,nodej) * T;
            VL(counter) = sides(:,nodei)'*sides(:,nodej)/(4*T);
            counter = counter + 1;
        end
    end
end
M = sparse(I,J,VM);
L = sparse(I,J,VL);

u = (M+beta*L)\(M*f);
% Plot height data
subplot(1,3,3);
trisurf(triangles,px,py,u)
colormap winter;
shading interp;
lighting phong;
light;
view ( -40, 80 );
axis square