function [tau1,lambda] = adjointCurl(gridNodes,e,u,dIdB)

% [tau1,lambda] = adjointCurl(gridNodes,e,u,dIdB)
%          applies gradient of the curl of the scalar potential u
%          defined at the nodes gridNodes of the triangulation e
%          with respect to the gridNodes and with respect to u,
%          respectively, to the vector dIdB, and stores the results
%          to tau1 and lambda.


n = max(size(u));
ne = max(size(e));

tau1 = zeros(size(gridNodes,1)*size(gridNodes,2),1);
lambda = zeros(n,1);

for i=1:ne
    ei = e(i,:);
    
    x1 = gridNodes(:,ei(1));
    x2 = gridNodes(:,ei(2));
    x3 = gridNodes(:,ei(3));
    
    R = [x2-x1 x3-x1];
    invR = inv(R);
    detR = det(R);
    graddetR = [x2(2)-x3(2) x3(2)-x1(2) x1(2)-x2(2);x3(1)-x2(1) x1(1)-x3(1) x2(1)-x1(1)]';
    Rtilde = [x3(2)-x1(2) x1(1)-x3(1);x1(2)-x2(2) x2(1)-x1(1)];
    
    dinvRdx11 = (1/power(detR,2)) * ( detR*[0 1;0 -1] - graddetR(1,1)*Rtilde );
    dinvRdx12 = (1/power(detR,2)) * ( detR*[-1 0;1 0] - graddetR(1,2)*Rtilde );
    dinvRdx21 = (1/power(detR,2)) * ( detR*[0 0;0 1] - graddetR(2,1)*Rtilde );
    dinvRdx22 = (1/power(detR,2)) * ( detR*[0 0;-1 0] - graddetR(2,2)*Rtilde );
    dinvRdx31 = (1/power(detR,2)) * ( detR*[0 -1;0 0] - graddetR(3,1)*Rtilde );
    dinvRdx32 = (1/power(detR,2)) * ( detR*[1 0;0 0] - graddetR(3,2)*Rtilde );
    
    
    u1 = u(ei(1));
    u2 = u(ei(2));
    u3 = u(ei(3));
    
    tau1(2*ei(1)-1) = tau1(2*ei(1)-1) + ((u2-u1)*dinvRdx11(1,2)+(u3-u1)*dinvRdx11(2,2))*dIdB(2*i-1) - ((u2-u1)*dinvRdx11(1,1)+(u3-u1)*dinvRdx11(2,1))*dIdB(2*i);
    tau1(2*ei(1)) = tau1(2*ei(1)) + ((u2-u1)*dinvRdx12(1,2)+(u3-u1)*dinvRdx12(2,2))*dIdB(2*i-1) - ((u2-u1)*dinvRdx12(1,1)+(u3-u1)*dinvRdx12(2,1))*dIdB(2*i);
    tau1(2*ei(2)-1) = tau1(2*ei(2)-1) + ((u2-u1)*dinvRdx21(1,2)+(u3-u1)*dinvRdx21(2,2))*dIdB(2*i-1) - ((u2-u1)*dinvRdx21(1,1)+(u3-u1)*dinvRdx21(2,1))*dIdB(2*i);
    tau1(2*ei(2)) = tau1(2*ei(2)) + ((u2-u1)*dinvRdx22(1,2)+(u3-u1)*dinvRdx22(2,2))*dIdB(2*i-1) - ((u2-u1)*dinvRdx22(1,1)+(u3-u1)*dinvRdx22(2,1))*dIdB(2*i);
    tau1(2*ei(3)-1) = tau1(2*ei(3)-1) + ((u2-u1)*dinvRdx31(1,2)+(u3-u1)*dinvRdx31(2,2))*dIdB(2*i-1) - ((u2-u1)*dinvRdx31(1,1)+(u3-u1)*dinvRdx31(2,1))*dIdB(2*i);
    tau1(2*ei(3)) = tau1(2*ei(3)) + ((u2-u1)*dinvRdx32(1,2)+(u3-u1)*dinvRdx32(2,2))*dIdB(2*i-1) - ((u2-u1)*dinvRdx32(1,1)+(u3-u1)*dinvRdx32(2,1))*dIdB(2*i);
    
    lambda(ei(1)) = lambda(ei(1)) + (-invR(1,2)-invR(2,2))*dIdB(2*i-1) + (invR(1,1)+invR(2,1))*dIdB(2*i);
    lambda(ei(2)) = lambda(ei(2)) + invR(1,2)*dIdB(2*i-1) - invR(1,1)*dIdB(2*i);
    lambda(ei(3)) = lambda(ei(3)) + invR(2,2)*dIdB(2*i-1) - invR(2,1)*dIdB(2*i);
    
%    B(1,i) = (u2-u1)*invR(1,2) + (u3-u1)*invR(2,2);
%    B(2,i) = - ( (u2-u1)*invR(1,1) + (u3-u1)*invR(2,1) );
end
