function [A,f] = NewAssemble(TR,g)
% returns stiffness and load matrix for given triangulation class TR
% g is the function that is right hand side of Poisson equation.
N = size(TR.ConnectivityList,1);
M = size(TR.Points,1);
d = size(TR.Points,2);
% Initialize
A = zeros(M);
f = zeros(M,1);
% Phi derivatives of Lagrange basis functions wrt reference coords
Phi = [eye(d),-ones(d,1)]';
% Location of barycentre, for quadrature
lambda = ones(1,d+1)'/(d+1);
for i=1:N
    ind = TR.ConnectivityList(i,:);
    % X is a d x d+1 matrix with columns the vertices x_i
    X = TR.Points(ind,:)';
    % Jacobian of forward mapping
    % An equivalent code (2d case only): J = [X(:,1)-X(:,3),X(:,2)-X(:,3)];
    J = bsxfun(@minus, X(:,1:d),X(:,d+1));
    Jac = abs(det(J));
    % Gradient matrix, a d+1 x d matrix
    G = Phi / J;
    % element stiffness
    A1 = (G*G')*Jac/factorial(d);
    % element load by 1-point barycentre quadrature
    barycentre = X*lambda;
    f1 = lambda*g(barycentre)*Jac/factorial(d);
    % insert element stiffness/load into main matrix/vector
    A(ind,ind)=A(ind,ind)+A1;
    f(ind) = f(ind) + f1;
end