function [A,b] = fem2dPlaneStress(P,T,ID,IN,lambdamu,f,g)

N = size(P,2);
m = size(T,2);

Atilda = sparse(2*N,2*N);
btilda = zeros(2*N,1);
for l=1:m
    idx = [T(:,l);N+T(:,l)];
    Pl = P(:,T(:,l));
    
    Al = localStiffnessMatrix2dPlaneStress(lambdamu(:,l),Pl);
    Atilda(idx,idx) = Atilda(idx,idx) + Al;
%     spy(Atilda)
%     pause
    
    bl = localLoadVector2dPlaneStress(f(:,l),Pl);
    btilda(idx) = btilda(idx)+bl;
end

mN = size(IN,2);
for l=1:mN
    idx = [IN(:,l);N+IN(:,l)];
    Pl = P(:,IN(:,l));
    
    cl = localBoundaryLoadVector2dPlaneStress(g(:,l),Pl);
    btilda(idx) = btilda(idx)+cl;
%     btilda
%     pause
end

nonDirichletIdx = setdiff(1:N,ID);
nonDirichletIdx = [nonDirichletIdx,N+nonDirichletIdx];
A = Atilda(nonDirichletIdx,nonDirichletIdx);
b = btilda(nonDirichletIdx);
