function [A,B] = formAb(f,p,t,nip)

% fe formulation A u = B

[ntriangles,n] = size(t);

A = sparse(nip,nip);
B = zeros(nip,1);

for i=1:ntriangles,

% Indices t and coordinates p of triangle

   t1 = t(i,1);
   t2 = t(i,2);
   t3 = t(i,3);
   p1 = p(t1,:);
   p2 = p(t2,:);
   p3 = p(t3,:);

% Area of triangle

   a = norm(p1-p2);
   b = norm(p2-p3);
   c = norm(p3-p1);
   s = (a+b+c)/2;
   trianglearea = sqrt(s*(s-a)*(s-b)*(s-c));

% Matrix and rhs entries
% using nodal quadrature rule

   grad1 = gradplate(p1,p2,p3);
   grad2 = gradplate(p2,p1,p3);
   grad3 = gradplate(p3,p2,p1);

   if (t1 <= nip)
      A(t1,t1) = A(t1,t1)+(grad1'*grad1)*trianglearea/2;
      B(t1) = B(t1)+feval(f,p1)*trianglearea/3;
      if (t2 <= nip)
         A(t1,t2) = A(t1,t2)+(grad1'*grad2)*trianglearea;
      end
   end
   if (t2 <= nip)
      A(t2,t2) = A(t2,t2)+(grad2'*grad2)*trianglearea/2;
      B(t2) = B(t2)+feval(f,p2)*trianglearea/3;
      if (t3  <= nip)
         A(t2,t3) = A(t2,t3)+(grad2'*grad3)*trianglearea;
      end
   end
   if (t3 <= nip)
      A(t3,t3) = A(t3,t3)+(grad3'*grad3)*trianglearea/2;
      B(t3) = B(t3)+feval(f,p3)*trianglearea/3;
      if (t1 <= nip)
         A(t1,t3) = A(t1,t3)+(grad1'*grad3)*trianglearea;
      end
   end


end

A = A + A';
