b = Eigen::VectorXd::Zero(V.rows()*2);
    Eigen::SparseMatrix<double> Dx, Dy, U, D, B;
    computeSurfaceGradientMatrix(Dx, Dy);
    igl::cat(2, Dx, Eigen::SparseMatrix<double>(-Dy), U);
    igl::cat(2, Dx, Dy, D);
    igl::cat(1, U, D, B);
    // Area diagonal
    Eigen::VectorXd areas;
    igl::doublearea(V, F, areas);
    Eigen::SparseMatrix<double> AA(areas.size()*2, areas.size()*2);
    for(int i=0; i<areas.size()*2; ++i){
      AA.insert(i, i) = std::sqrt(areas(i%areas.size()));
    }
    // Put it together
    A = (B.transpose()*AA)*B;