function poisson1d(red)
d=1;
[c4n,n4e,Db,Nb] = triang_cube(1);
for red = 1 : red
    [c4n,n4e,Db,Nb] = red_refine(c4n,n4e,Db,Nb);
end
nC = size(c4n,1); nE = size(n4e,1);
dNodes = unique(Db); fNodes = setdiff(1:nC,dNodes);

%% P1 Poisson
A = sparse(nC,nC);
u = zeros(nC,1); b = zeros(nC,1);
for j = 1 : nE
    X_T = [ones(1,d+1);c4n(n4e(j,:),:)'];
    grads_T = X_T\[zeros(1,d);eye(d)];
    vol_T = det(X_T)/factorial(d);
    mp_T = sum(c4n(n4e(j,:),:),1)/(d+1);
    for m = 1:d+1
        b(n4e(j,m)) = b(n4e(j,m))+(1/(d+1))*vol_T*f(mp_T);
        for n = 1:d+1
            A(n4e(j,m),n4e(j,n)) = A(n4e(j,m),n4e(j,n)) + ...
                vol_T * grads_T(m,:) * grads_T(n,:)';
        end
    end
end
u(fNodes) = A(fNodes,fNodes)\b(fNodes);
close all;
fplot(@(x) x.*(1-x)./2,[0 1],'g'); hold on;
plot(c4n(n4e),u(n4e),'g');


%% P1-P1
A = sparse(nC,nC);
B = sparse(nC,nC);
g = zeros(nC,1);
for j = 1: nE
    X_T = [ones(1,d+1);c4n(n4e(j,:),:)'];
    grads_T = X_T\[zeros(1,d);eye(d)];
    vol_T = det(X_T)/factorial(d);
    mp_T = sum(c4n(n4e(j,:),:),1)/(d+1);
    m_loc = (ones(d+1,d+1)+eye(d+1,d+1))/((d+1)*(d+2));
    for m = 1:d+1
         g(n4e(j,m)) = g(n4e(j,m))+(1/(d+1))*vol_T*f(mp_T);
        for n = 1:d+1
            A(n4e(j,m),n4e(j,n)) = A(n4e(j,m),n4e(j,n)) + ...
                vol_T*m_loc(m,n);
            B(n4e(j,m),n4e(j,n)) = B(n4e(j,m),n4e(j,n)) + ...
                vol_T * grads_T(m,:) * (1/2);
        end
    end
end
M = [A,B;B',zeros(nC,nC)];
b = [zeros(nC,1);-g];
pu = M\b;
u = pu(nC+(1:nC));
plot(c4n(n4e),u(n4e),'b');
maxosc  = max(u(1:(nC-1))-u(2:nC))

%% P1-P0
BB = sparse(nE,nC);
gg = zeros(nE,1);
AA = A;
for j = 1: nE
    X_T = [ones(1,d+1);c4n(n4e(j,:),:)'];
    grads_T = X_T\[zeros(1,d);eye(d)];
    vol_T = det(X_T)/factorial(d);
    mp_T = sum(c4n(n4e(j,:),:),1)/(d+1);
    gg(j) = vol_T*f(mp_T);
    for m = 1:d+1
          BB(j,n4e(j,m)) = BB(j,n4e(j,m)) + vol_T * grads_T(m,:);
    end
end
MM = [AA,BB';BB,zeros(nE,nE)];
bb = [zeros(nC,1);-gg];
ppuu = MM\bb;
uu = ppuu(nC+(1:nE));
stairs([0;c4n(n4e(:,1));1],[0;uu;0]);
%vv = [0;uu(1:(nE-1))];
%plot(c4n(n4e),[vv,uu(nC+(1:nE))],'g');
end

function val = f(x)
    val = 1; 
end
