function finelt

beta = 2E4;

% read triangles
file = fopen ( 'gitter.txt', 'r' );
npt = cell2mat ( textscan ( file, '%d,%d', 1 ) );
np = npt ( 1 ); nt = npt ( 2 );
TX = cell2mat ( textscan ( file, '%f', np, 'delimiter', ',' ) );
TY = cell2mat ( textscan ( file, '%f', np, 'delimiter', ',' ) );
TH = cell2mat ( textscan ( file, '%f', np, 'delimiter', ',' ) );
TT = cell2mat ( textscan ( file, '%d,%d,%d', nt ) );

% preallocate matrix data
nijv = 9 * nt;
LI = zeros ( nijv, 1 );
LJ = zeros ( nijv, 1 );
LV = zeros ( nijv, 1 );
MI = zeros ( nijv, 1 );
MJ = zeros ( nijv, 1 );
MV = zeros ( nijv, 1 );

% precompute some values
A = zeros ( nt, 1 ); % area of triangle
C = zeros ( nt, 3 ); % cosine of angle at node
S = zeros ( nt, 3 ); % sine of angle at node
H = zeros ( nt, 3 ); % slope of basis function
for t = 1 : nt
    T = TT ( t, : );
    for s = 1 : 3
        % vector pointing from node s to next node
        dx = TX ( T ( s ) ) - TX ( T ( n ( s ) ) );
        dy = TY ( T ( s ) ) - TY ( T ( n ( s ) ) );
        E ( s, : ) = [ dx dy ]; 
        L ( s ) = norm ( E ( s, : ) );
        D ( s, : ) = E ( s, : ) * 1 / L ( s );
    end
    for s = 1 : 3
        C ( t, s ) = D ( p ( s ), : ) * D ( s, : )';
        S ( t, s ) = det ( D ( [ ( p ( s ) ) s ], : ) );
    end
    for s = 1 : 3
        H ( t, s ) = 1 / ( L ( s ) * S ( t, n ( s ) ) );
    end
    A ( t ) = 0.5 * abs ( det ( E ( 1 : 2, : ) ) );
end

% assemble mass matrix
Mloc = [ 2 1 1; 1 2 1; 1 1 2 ];
k = 1;
for t = 1 : nt
    T = TT ( t, : );
    fac = A ( t ) / 12.;
    for i = 1 : 3
        for j = 1 : 3
            val = fac * Mloc ( i, j );
            MI ( k ) = T ( i ); MJ ( k ) = T ( j ); MV ( k ) = val; k = k + 1;
        end
    end
end
MI ( k : nijv ) = [];
MJ ( k : nijv ) = [];
MV ( k : nijv ) = [];
M = sparse ( MI, MJ, MV );

% assemble stiffness matrix
O = [ 0 3 2; 3 0 1; 2 1 0 ]; % O ( i, j ) is index of other vertex for i ~= j
k = 1;
for t = 1 : nt
    T = TT ( t, : );
    fac = A ( t );
    for i = 1 : 3
        val = fac * H ( t, i ) ^ 2;
        LI ( k ) = T ( i ); LJ ( k ) = T ( i ); LV ( k ) = val; k = k + 1;
        for j = i + 1 : 3 
            val = fac * H ( t, i ) * H ( t , j ) * C ( t, O ( i, j ) );
            LI ( k ) = T ( i ); LJ ( k ) = T ( j ); LV ( k ) = val; k = k + 1;
            LI ( k ) = T ( j ); LJ ( k ) = T ( i ); LV ( k ) = val; k = k + 1;
        end
    end
end
LI ( k : nijv ) = [];
LJ ( k : nijv ) = [];
LV ( k : nijv ) = [];
L = sparse ( LI, LJ, LV );

% solve
MTH = M * TH;
LPM = beta * L + M;
Z = LPM \ MTH;

% visualize
colormap ( winter );
h1 = subplot ( 221 );
trisurf ( TT, TX, TY, TH );
view ( 2 );
h2 = subplot ( 222 );
trisurf ( TT, TX, TY, Z );
view ( 2 );
axis ( [ h1 h2 ], 'square' );
subplot ( 223 );
trisurf ( TT, TX, TY, TH );
shading interp;
lighting phong;
light;
view ( -40, 60 );
subplot ( 224 );
trisurf ( TT, TX, TY, Z );
shading interp;
lighting phong;
light;
view ( -40, 60 );

% index of next and previous vertex
function n = n ( s )
n = mod ( s, 3 ) + 1;
function p = p ( s )
p = mod ( s + 1, 3 ) + 1;