function P1FE_poisson(red)

%creation of triangulation
c4n = [0,0;1,0;1,1;0,1];
n4e = [1,2,3;1,3,4];
Db = [1,2;3,4;4,1];
Nb = [2,3];

for i = 1 : red
    [c4n,n4e,Db,Nb] = red_refine(c4n,n4e,Db,Nb);
end

%plot triangulation
figure(1); clf;
trimesh(n4e,c4n(:,1),c4n(:,2));

%required constants
[nC,d] = size(c4n);
nE = size(n4e,1);

DiriNodes = unique(Db);
freeNodes = setdiff( 1:nC, DiriNodes );

%initialization of vectors for sparse matrix A
ctr_max = nE * (d + 1)^2;
I = zeros(ctr_max,1);
J = zeros(ctr_max,1);
X = zeros(ctr_max,1);
ctr = 0;

%right hand side and solution vector
b = zeros(nC,1);
x = zeros(nC,1);

%assembly of stiffness matrix A and right hand side b
for j = 1 : nE
    grads_T = [ones(1,d+1);c4n(n4e(j,:),:)'] \ [zeros(1,d);eye(d)];
    area_T = det([ones(1,d+1);c4n(n4e(j,:),:)']) / factorial(d);
    mp_T = sum(c4n(n4e(j,:),:)) / (d+1);
    for m = 1 : d + 1
        b(n4e(j,m)) = b(n4e(j,m)) + (1/3) * area_T * f(mp_T);
        for n = 1 : d + 1
            ctr = ctr + 1;
            I(ctr) = n4e(j,m);
            J(ctr) = n4e(j,n);
            X(ctr) = area_T * grads_T(m,:) * grads_T(n,:)';
        end
    end
end
A = sparse(I,J,X,nC,nC);

%neumann boundary condition
for j = 1 : size(Nb,1)
    length_S = norm(c4n(Nb(j,1),:) - c4n(Nb(j,2),:));
    mp_S = (c4n(Nb(j,1),:) + c4n(Nb(j,2),:))/2;
    b(Nb(j,1)) = b(Nb(j,1)) + (1/2) * length_S * g(mp_S);
    b(Nb(j,2)) = b(Nb(j,2)) + (1/2) * length_S * g(mp_S);
end

%dirichlet boundary condition;
for j = 1 : size(DiriNodes,1)
    x(DiriNodes(j)) = u_D(c4n(DiriNodes(j),:));
end
b = b - A * x;

%computation of the solution for non-dirichlet-boundary nodes
x(freeNodes) = A(freeNodes,freeNodes) \ b(freeNodes);

%plot solution
figure(2); clf;
trisurf(n4e,c4n(:,1),c4n(:,2),x);


%function for neumann boundary condition
function val = g(x)
val = 2 * pi;

%function for dirichlet boundary condition
function val = u_D(x)
val = sin(2*pi*x(1)) * cos(2 * pi * x(2));

%function for right hand side
function val = f(x)
val = 5;
