This is the code fragment

PositionBasedCosseratRods.cpp

bool PositionBasedCosseratRods::solve_StretchShearConstraint(
	const glm::vec3& p0, float invMass0,
	const glm::vec3& p1, float invMass1,
	const glm::quat& q0, float invMassq0,
	const glm::vec3& stretchingAndShearingKs,
	const float restLength,
	glm::vec3& corr0, glm::vec3&  corr1, glm::quat&  corrq0)
{
	glm::vec3 d3;	//third director d3 = q0 * e_3 * q0_conjugate
	d3[0] = 2.0 * (q0.x * q0.z + q0.w * q0.y);
	d3[1] = 2.0 * (q0.y * q0.z - q0.w * q0.x);
	d3[2] = q0.w * q0.w - q0.x * q0.x - q0.y* q0.y + q0.z * q0.z;
	//d3 = glm::normalize(d3);

	glm::vec3 gamma = (p1 - p0) / restLength - d3;
	gamma /= (invMass1 + invMass0) / restLength + invMassq0 * 4.0*restLength + 1.0e-6;
	for (int i = 0; i<3; i++) gamma[i] *= stretchingAndShearingKs[i];

	corr0 = invMass0 * gamma;
	corr1 = -invMass1 * gamma;
	//quat(w,x,y,z)
	glm::quat q_e_3_bar(q0.z, -q0.y, q0.x, -q0.w);//compute q*e_3.conjugate (cheaper than quaternion product)
	corrq0 = glm::quat(0.0, gamma.x, gamma.y, gamma.z) * q_e_3_bar;
	corrq0 *=  invMassq0 * restLength;

	return true;
}
bool PositionBasedCosseratRods::solve_BendTwistConstraint(
	const glm::quat& q0, float invMassq0,
	const glm::quat& q1, float invMassq1,
	const glm::vec3& bendingAndTwistingKs,
	const glm::quat& restDarbouxVector,
	glm::quat& corrq0, glm::quat& corrq1)
{
	glm::quat omega = glm::conjugate(q0) * q1;   //darboux vector

	glm::quat omega_plus;
	omega_plus = omega + restDarbouxVector;     //delta Omega with -Omega_0
	omega = omega - restDarbouxVector;          //delta Omega with + omega_0

	float omega_plus_norm = glm::dot(omega_plus, omega_plus);
	float omega_norm = glm::dot(omega, omega);

	if (omega_norm>omega_plus_norm) 
		omega = omega_plus;

	for (int i = 0; i < 3; i++) omega[i] *= bendingAndTwistingKs[i] / (invMassq0 + invMassq1 + 1.0e-6);
	omega.w = 0.0;    //discrete Darboux vector does not have vanishing scalar part

	corrq0 = q1 * omega;
	corrq1 = q0 * omega;
	corrq0 *= invMassq0;
	corrq1 *= -invMassq1;
	return true;
}

TimeController.cpp

void TimeStepController::positionConstraintProjection(SimulationModel &model)
{
	unsigned int iter = 0;

	// init constraint groups if necessary
	model.initConstraintGroups();

	SimulationModel::ConstraintVector &constraints = model.getConstraints();
	SimulationModel::ConstraintGroupVector &groups = model.getConstraintGroups();

	/* Solving the constraints in a sequential order from one end of the rod to the other one can lead to instablilities
	*  bilaterial interleaving order can overcome this problem
	*/
	while (iter < m_maxIter)
	{
		for (unsigned int group = 0; group < groups.size(); group++)
		{
			if (group % 2 == 0)
			{
				const int groupSize = (int)groups[group].size();
				{
					for (int i = 0; i < groupSize; i++)
					{
						const unsigned int constraintIndex = groups[group][i];

						constraints[constraintIndex]->updateConstraint(model);
						constraints[constraintIndex]->solvePositionConstraint(model);
					}
				}
			}
			else
			{
				const int groupSize = (int)groups[group].size();
				{
					for (int i = groupSize-1; i >= 0; i--)
					{
						const unsigned int constraintIndex = groups[group][i];

						constraints[constraintIndex]->updateConstraint(model);
						constraints[constraintIndex]->solvePositionConstraint(model);
					}
				}
			}

		}

		iter++;
	}
}

TimeIntegration.cpp

void TimeIntegration::semiImplicitEuler(
  const float h,
	const float mass,
	glm::vec3 &position,
	glm::vec3 &velocity,
	const glm::vec3 &acceleration)
{
	if (mass != 0.0)
	{
		velocity += acceleration * h;
		position += velocity * h;
	}
}

void TimeIntegration::semiImplicitEulerRotation(
	const float h,
	const float mass,
	const glm::mat3 &invInertiaW,
	glm::quat &rotation,
	glm::vec3 &angularVelocity,
	const glm::vec3 &torque)
{
	if (mass != 0.0)
	{
		// simple form without nutation effect
		angularVelocity += h * invInertiaW * torque;

		glm::quat angVelQ(0.0, angularVelocity[0], angularVelocity[1], angularVelocity[2]);
		rotation += (float)(h * 0.5) * angVelQ * rotation;
		glm::normalize(rotation);
	}
}

void TimeIntegration::velocityUpdateFirstOrder(
	const float h,
	const float mass,
	const glm::vec3 &position,
	const glm::vec3 &oldPosition,
	glm::vec3 &velocity)
{
	if (mass != 0.0)
		velocity = (float)(1.0 / h) * (position - oldPosition);
}

void TimeIntegration::angularVelocityUpdateFirstOrder(
	const float h,
	const float mass,
	const glm::quat &rotation,
	const glm::quat &oldRotation,
	glm::vec3 &angularVelocity)
{
	if (mass != 0.0)
	{
		const glm::quat relRot = (rotation * glm::conjugate(oldRotation));
		angularVelocity = glm::vec3(relRot.x, relRot.y, relRot.z) *(float)(2.0 / h);
	}
}

Position and Orientation Cosserat Rods