Code:

void constrainedForwardEulerStep(ConstrainedSystem *system, double dt) {    forwardEulerStep(system, dt);    projectPositions(system);    projectVelocities(system);}

Code:

void projectPositions(ConstrainedSystem *system) {    int n = system->getDOFs(),        m = system->getConstraints();    static VectorXd c(m);system->getConstraintValues(c);     static MatrixXd J(m,n);system->getConstraintJacobian(J);    static MatrixXd Mi(n, n);system->getInverseInertia(Mi);    MatrixXd A = J*Mi*J.transpose();    VectorXd lambda = solve(A, -c);     static VectorXd dx(n);    dx = Mi*J.transpose()*lambda;        static VectorXd x(n), v(n);    system->getState(x, v);    system->setState(x + dx, v);}

Code:

void projectVelocities(ConstrainedSystem *system) {    int n = system->getDOFs(),         m = system->getConstraints();    static VectorXd x(n), v(n);    system->getState(x, v);    v = v * 0.995; //global dampening        static MatrixXd J(m,n);    system->getConstraintJacobian(J);        static MatrixXd Mi(n,n);    system->getInverseInertia(Mi);        MatrixXd A = J*Mi*J.transpose();        VectorXd b = -J*v;    VectorXd lambda = solve(A, b);        static VectorXd dv(n);    dv = Mi*J.transpose()*lambda;        system->setState(x, v + dv);}

Code:

 void FastProjection(ParticleSystem* system, float dt) {   //allocate memory for X0,Xt,Vt, dx, C, J, Mi;     //get initial values    system->getInverseInertia(Mi);    system->getState(Xt,Vt);   //outer solver iteration loop   for(int j=0; j<MAX_ITERATIONS;++j ) {            X0 = Xt+dt*Vt; //take initial unconstrained step      system->setState(X0,Vt);      for (int i=0; i<constraints.size(); ++i ) {         const Constraint& c = constraints[i];         float strain(0);         do {            //calculate strain            Vector2d delta = X0[c.p1]-X0[c.p0];     strain = delta.mag() / c.rest_length;                              //obtain lambda by solving linear system             system->getConstraintJacobian(J);            system->getConstraintValues(C);            Jt=J.transpose();                         MatrixXd A = J*Mi*Jt; //from eq. 7             VectorXd b = C;             VectorXd lambda = solve(A, b);             static VectorXd dx(n); //from eq. 5             dx = Mi*Jt*lambda;                         X0 = X0 + dx;              system->setState(X0,Vt);        } while(strain > STRAIN_THRESHOLD);         } //end for all constraints     }//end for all solver iterations      V = (X0-Xt)/dt;   system->setState(X0,V);}