function [ sysmat,sysvec ] = gen_system( ecken, kanten, knoten )
%GEN_SYSTEM Generate System matrix and vector, homogeneous Dirichlet.
%Note: The transformation from local to global matrices works
%for triangular Lagrange elements only, but can easily be expanded.
close all;
if (nargin==0)
     %Use the external triangulation algorithm.
     %[ecken,kanten,knoten]=triangulate_circ(30);
     %[ecken,kanten,knoten]=triangulate();
     gd=[4;0;0;1;1;0];sf='E1'; ns=[69;49];
     g=decsg(gd,sf,ns);
     [ecken,kanten,knoten]=initmesh(g);
end
pdemesh(ecken,kanten,knoten);
figure;
%waitforbuttonpress;
ecken=ecken'; kanten=kanten'; knoten=knoten';

%Compute system matrix and vector. For other elements, the number of base
%functions must be changed.
sysmat=sparse(size(ecken,1),size(ecken,1));
sysvec=zeros(size(ecken,1),1);
for i=1:size(knoten,1)
    x0=ecken(knoten(i,1),:)';x1=ecken(knoten(i,2),:)';x2=ecken(knoten(i,3),:)';
    [elmat,elvec]=compute_elmat_elvec(x0,x1,x2);
    for j=1:3
        for l=1:3
           sysmat(knoten(i,j),knoten(i,l))=sysmat(knoten(i,j),knoten(i,l))+elmat(j,l);
        end
        sysvec(knoten(i,j))=sysvec(knoten(i,j))+elvec(j);
    end
end

%Homogeneous Dirichlet boundary conditions.
%kanten=kanten(:,1:2);
for i=1:size(ecken,1)
    if (numel(find((kanten(:,1)==i|kanten(:,2)==i)&(kanten(:,6)==0|kanten(:,7)==0)))>0)
        sysmat(i,:)=0; sysmat(i,i)=1; sysvec(i)=0;
    end
%    if (numel(find(kanten=)>0)
%        sysmat(i,:)=0; sysmat(i,i)=1; sysvec(i)=0;
%    end    
end

%compute the solution
z=inv(sysmat)*sysvec;

%plot the result using voronoi.
x=ecken(:,1);
y=ecken(:,2);
N=100;
ti = (0:N)/N;
[XI,YI] = meshgrid(ti,ti);
YI=1-YI;
ymin=min(y);
ymax=max(y);
xmin=min(x);
xmax=max(x);
YI=(ymax-ymin)*YI+ymin;
XI=(xmax-xmin)*XI+xmin;
ZI = griddata(x,y,z,XI,YI,'cubic'); 
ZI = reshape(ZI,N+1,N+1);
imagesc(ZI); colorbar();
end

%Define the PDE.
function  y=f(x)
%This is the right hand side. x is in Omega.
y=4;
end
function Y=D(x)
%This must be negative definite.
Y=eye(2);
end
function y=ce(x)
%This is the constant term.
y=0;
end

%Computation of element matrices/vectors.
function [A, b, B, c]=compute_transformation_matrices(x0,x1,x2)
b=x0;
A=[x1-x0 x2-x0]; B=inv(A);
c=-B*b;
end
function [elmat elvec]=compute_elmat_elvec(x0,x1,x2)
%We integrate the reference triangle by evaluating at the center
%and multiplying with its size. Insert better integration formula here.
u=[1./3.;1./3.];
[A,b,B,c]=compute_transformation_matrices(x0,x1,x2);
[elmat,elvec]=integrand(u,A,b,B,c);
elmat=elmat/2.;
elvec=elvec/2.;
end
function [bsmat,bsvec]=integrand(x,A,b,B,c)
%compute integrand at x, x is on the reference triangle.
%Should share J,Y!!!!
J=refJakobimat(x);
Y=refformfunc(x);
ximg=A*x+b;
bsmat=(J*B*D(ximg)*B'*J'+ce(ximg)*Y*Y')*abs(det(A));
bsvec=f(ximg)*Y*abs(det(A));
end
function y=locformfunc(x,A,b,B,c)
y=refformfunc(B*x+c);
end

%And now: Linear Lagrange Triangular!
function A=refJakobimat(x)
A=zeros(3,2); A(1,1)=-1; A(1,2)=-1; A(2,1)= 1; A(2,2)= 0; A(3,1)= 0; A(3,2)= 1;
end
function y=refformfunc(x)
y=zeros(3,1); y(1)=1-x(1)-x(2); y(2)=x(1); y(3)=x(2);
end