Orbots

void computeGradients(const Vector3& pa, const Vector3& pb, const Vector3& pc, const Vector3& pd,
                      Vector3& grad_pa, Vector3& grad_pb, Vector3& grad_pc, Vector3& grad_pd) {
    // Step 1: Compute edge vectors using operator-
    Vector3 eca = pc - pa;
    Vector3 eda = pd - pa;
    Vector3 edb = pd - pb;
    Vector3 ecb = pc - pb;
    Vector3 pd_pc = pd - pc;
    double L = pd_pc.Norm();
    Vector3 ecd = (L < 1e-8) ? Vector3(0.0, 0.0, 0.0) : (pd_pc * (1.0 / L));