function [x,A] = CFEM(c4n,n4e,n4Db)
n = size(c4n,1); m = size(n4e,1); d = size(c4n,2);
A = zeros(d+1,d+1,m); b = zeros(m,d+1); x = zeros(n,1);
for j = 1 : m
    area = abs(det([ones(1,d+1);c4n(n4e(j,:),:)'])/prod(2:d));
    grads = [ones(1,d+1);c4n(n4e(j,:),:)']\[zeros(1,d);eye(d)];
    A(:,:,j) = area*(grads*grads');
    b(j,:) = ones(d+1,1).*area/(d+1);
end
J = n4e(:,cumsum(ones(d+1,d+1)))'; K = n4e(:,cumsum(ones(d+1,d+1))')';
A = sparse(J(:),K(:),A(:),n,n);
b = accumarray(n4e(:),b(:),[n 1]);
dof = setdiff(1:n,n4Db(:));
x(dof) = A(dof,dof)\b(dof);
end