function G = extractG(VList, FList)

nV = size(VList,1);
nF = size(FList,1);

%G = zeros(nF*4, nV*2);

nz = nF*36;
rowi = zeros(nz,1);
coli = zeros(nz,1);
vali = zeros(nz,1);

rowId = 1; idx=1;
for i=1:nF
    fvids = FList(i,:);
    a = VList(fvids(1),1:3)';
    b = VList(fvids(2),1:3)';
    c = VList(fvids(3),1:3)';
    d = VList(fvids(4),1:3)';
    restX = [b-a c-a d-a];
    invX = inv(restX);
    
    for cj=1:3
        for ci=1:3                                          % 3*(vid-1)+ci-1+1
           rowi(idx,1) = rowId+3*(cj-1)+ci-1; coli(idx,1) = 3*fvids(1)+ci-3; vali(idx,1) = -invX(1,cj)-invX(2,cj)-invX(3,cj); idx = idx + 1;%G(rowId+2*(ci-1)+cj-1,2*fvids(1)+ci-2) = -invX(1,cj)-invX(2,cj);
           rowi(idx,1) = rowId+3*(cj-1)+ci-1; coli(idx,1) = 3*fvids(2)+ci-3; vali(idx,1) = invX(1,cj); idx = idx + 1;%G(rowId+2*(ci-1)+cj-1,2*fvids(2)+ci-2) = invX(1,cj);
           rowi(idx,1) = rowId+3*(cj-1)+ci-1; coli(idx,1) = 3*fvids(3)+ci-3; vali(idx,1) = invX(2,cj); idx = idx + 1;%G(rowId+2*(ci-1)+cj-1,2*fvids(3)+ci-2) = invX(2,cj);
           rowi(idx,1) = rowId+3*(cj-1)+ci-1; coli(idx,1) = 3*fvids(4)+ci-3; vali(idx,1) = invX(3,cj); idx = idx + 1;
        end
    end
    rowId = rowId + 9;
    
end

G = sparse(rowi, coli, vali, 9*nF, 3*nV);


end

