function [x,A,dof4e] = DGFEM(c4n,n4e,n4ed,e4ed,eta)
m = size(n4e,1); n = size(e4ed,1);
A = zeros(3,3,m); b = zeros(m,3);
dof4e = reshape(1:3*m,m,3);
for j = 1 : m
    area = abs(det([ones(1,3);c4n(n4e(j,:),:)'])/2);
    grads = [ones(1,3);c4n(n4e(j,:),:)']\[zeros(1,2);eye(2)];
    A(:,:,j) = area*(grads*grads');
    b(j,:) = area/3;
end
% edge contributions
Ac = zeros(6,4,n);
Ap = zeros(4,4,n);
dof4ed = ones(n,6);
for j = 1:n
    % local edge patch enumeration
    nodes = zeros(6,1);
    nodes(1) = find(n4e(e4ed(j,1),:)==n4ed(j,1));
    nodes(2) = find(n4e(e4ed(j,1),:)==n4ed(j,2));
    nodes(3) = find(n4e(e4ed(j,1),:)~=n4ed(j,1) & n4e(e4ed(j,1),:)~=n4ed(j,2));
    if e4ed(j,2)~=0 % interior edge
        nodes(4) = find(n4e(e4ed(j,2),:)==n4ed(j,1));
        nodes(5) = find(n4e(e4ed(j,2),:)==n4ed(j,2));
        nodes(6) = find(n4e(e4ed(j,2),:)~=n4ed(j,1) & n4e(e4ed(j,2),:)~=n4ed(j,2));
    end
    % consistency
    t = c4n(n4ed(j,2),:)-c4n(n4ed(j,1),:);% not normalized!
    [~,indx] = find(n4e(e4ed(j,1),:)==n4ed(j,1)); %figure out orientation
    nE = (-1)^(n4e(e4ed(j,1),mod(indx,3)+1)~=n4ed(j,2)).*[t(2),-t(1)]; %normal pointing from T_1 to T_2
    grads1 = [ones(1,3);c4n(n4e(e4ed(j,1),nodes(1:3)),:)']\[zeros(1,2);eye(2)]; % grad on T_1
    if e4ed(j,2)~=0
        grads2 = [ones(1,3);c4n(n4e(e4ed(j,2),nodes(4:6)),:)']\[zeros(1,2);eye(2)]; % grad on T_2
        Ac(:,:,j) = 1/2*[sum(grads1.*nE,2);sum(grads2.*nE,2)].*[1 1 -1 -1]./2;
    else
        Ac(1:3,1:2,j) = sum(grads1.*nE,2).*[1 1]./2;
    end
    % penalty
    p = eta.*[1/3 1/6; 1/6 1/3];
    Ap(1:2,1:2,j) = p;
    if e4ed(j,2)~=0
        Ap(3:4,1:2,j) = -p;
        Ap(1:2,3:4,j) = -p;
        Ap(3:4,3:4,j) = p;
    end
    % edge dofs indices
    dof4ed(j,1:3) = dof4e(e4ed(j,1),nodes(1:3));
    if e4ed(j,2)~=0
        dof4ed(j,4:6) = dof4e(e4ed(j,2),nodes(4:6));
    end
end
% assemble volume
J = dof4e(:,cumsum(ones(3,3)))'; K = dof4e(:,cumsum(ones(3,3))')';
A = sparse(J(:),K(:),A(:),3*m,3*m);
b = accumarray(dof4e(:),b(:),[3*m 1]);
% assemble consistency
J = dof4ed(:,cumsum(ones(6,4)))';
dof4ed = dof4ed(:,[1 2 4 5]);
K = dof4ed(:,cumsum(ones(4,6))')';
Ac = sparse(J(:),K(:),Ac(:),3*m,3*m);
% assemble penalty
J = dof4ed(:,cumsum(ones(4,4)))'; K = dof4ed(:,cumsum(ones(4,4))')';
Ap = sparse(J(:),K(:),Ap(:),3*m,3*m);
A = A - Ac - Ac' + Ap;
x = A\b;
end