布料的模拟

布料的模拟有多种方法,可以实现,下面分别介绍之

1. explicit mass spring simulation

 下面看下关键部分的实现代码

void Simulator::ExplicitMassSpring()
{
	for(int i=0;i<m_meshpool.size();i++)
	{
		for( int j = 0; j < m_meshpool[i]->m_vertices_number; ++j )
			m_meshpool[i]->m_external_force.block_vector(j) = EigenVector3(0, g_gravity, 0);
	}

	for(int i=0; i<m_meshpool.size(); i++)
	{
		VectorX force_current;
		computeForcesWithDamping(m_meshpool[i]->m_current_positions, m_meshpool[i]->m_current_velocities, force_current, i);
		//computeForces(m_meshpool[i]->m_current_positions, force_current, i);
		m_meshpool[i]->m_current_velocities	= m_meshpool[i]->m_current_velocities + m_timestep * m_meshpool[i]->m_identity_matrix * force_current;
		m_meshpool[i]->m_current_positions = m_meshpool[i]->m_current_positions + m_meshpool[i]->m_current_velocities * m_timestep;	
	}
	
}

void Simulator::computeForcesWithDamping(const VectorX& x, const VectorX& v, VectorX& force, unsigned int meshid)
{
    VectorX constraintForce;

    constraintForce.resize(m_meshpool[meshid]->m_system_dimension);
    constraintForce.setZero();

    // springs
    for (map<pair<int, int>, Constraint*>::iterator it = m_meshpool[meshid]->m_mapConstraint.begin(); it != m_meshpool[meshid]->m_mapConstraint.end(); ++it)
    {
        it->second->EvaluateForceWithDamping(x, v, constraintForce);
    }

    // external forces
    constraintForce += m_meshpool[meshid]->m_external_force;

    force = constraintForce;
}

其中,damping force  和spring force的计算如下:

void AttachmentConstraint::EvaluateForceWithDamping(const VectorX& x, const VectorX& v, VectorX& force)
{
	EigenVector3 x_i = m_fixd_point - x.block_vector(m_p0);
	EigenVector3 g_i;
	if( x_i.isZero() )
	{
		g_i = (*(m_stiffness))* x_i;
	}
	else
	{
		g_i = (*(m_stiffness))* x_i - (*(m_damping)) * v.block_vector(m_p0).dot( x_i.normalized() ) * x_i.normalized();
	}
    
    force.block_vector(m_p0) += g_i;

}

void SpringConstraint::EvaluateForceWithDamping(const VectorX& x, const VectorX& v, VectorX& force)
{
    EigenVector3 x_ij = x.block_vector(m_p2) - x.block_vector(m_p1);
    EigenVector3 g_ij = (*(m_stiffness))*(x_ij.norm()-m_rest_length)*x_ij.normalized();
    force.block_vector(m_p1) += g_ij;
    force.block_vector(m_p2) -= g_ij;

	EigenVector3 v_ij = v.block_vector(m_p2) - v.block_vector(m_p1);

	if( x_ij.isZero() )
	{
	}
	else
	{
		g_ij = (*(m_damping)) * v_ij.dot(x_ij.normalized()) * x_ij.normalized();
		force.block_vector(m_p1) += g_ij;
		force.block_vector(m_p2) -= g_ij;
	}
 }

视频地址：https://youtu.be/dMgJ1AEy4Ko

2.implicit mass spring simulation

这里注意Ki,i的求法：

关键代码如下：

void Simulator::ImplicitMassSpring()
{
	for( int i = 0; i < m_meshpool.size(); ++i )
	{
		VectorX force_current;
		//computeForcesWithDamping(m_meshpool[i]->m_current_positions, m_meshpool[i]->m_current_velocities, force_current, i);
		computeForces(m_meshpool[i]->m_current_positions, force_current, i);
		//PrintVectorInDebugWindow(L"force_current", force_current);
		SparseMatrix K;
		computeJacobianOfF(m_meshpool[i]->m_current_positions, K, i);
		//PrintMatrixInDebugWindow( L"K", K);
		ScalarType deltaTSquare = m_timestep * m_timestep;
		SparseMatrix A = ( m_meshpool[i]->m_mass_matrix - deltaTSquare * K);
		VectorX b = m_meshpool[i]->m_mass_matrix * m_meshpool[i]->m_current_velocities + m_timestep * force_current;
		//PrintVectorInDebugWindow( L"b", b);
		Eigen::ConjugateGradient<SparseMatrix> cg;
		cg.compute(A);
		VectorX nextVel = cg.solve(b);
		m_meshpool[i]->m_current_positions = m_meshpool[i]->m_current_positions + m_timestep * nextVel;
		m_meshpool[i]->m_current_velocities = nextVel;
	}
}

void Simulator::computeJacobianOfF(const VectorX& x, SparseMatrix& K, unsigned meshid)
{
	K.resize( m_meshpool[meshid]->m_system_dimension, m_meshpool[meshid]->m_system_dimension);
	std::vector<SparseMatrixTriplet> vecTriplet;
	vecTriplet.clear();
	for (map<pair<int, int>, Constraint*>::iterator it = m_meshpool[meshid]->m_mapConstraint.begin(); it != m_meshpool[meshid]->m_mapConstraint.end(); ++it)
    {
        it->second->EvaluateJacobianOfF(x, vecTriplet);
    }
	K.setFromTriplets(vecTriplet.begin(), vecTriplet.end());
}

void AttachmentConstraint::EvaluateJacobianOfF(const VectorX& x, std::vector<SparseMatrixTriplet>& jacobian_triplets)
{
	ScalarType ks = *(m_stiffness);
	jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p0, 3 * m_p0, -ks));
	jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p0 + 1, 3 * m_p0 + 1, -ks));
	jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p0 + 2, 3 * m_p0 + 2, -ks));

}

void SpringConstraint::EvaluateJacobianOfF(const VectorX& x, std::vector<SparseMatrixTriplet>& jacobian_triplets)
{
	EigenVector3 x_ij = x.block_vector(m_p2) - x.block_vector(m_p1);
	ScalarType l_ij = x_ij.norm();
	ScalarType l0 = m_rest_length;
	ScalarType ks = *(m_stiffness);

	EigenMatrix3 k = ks * (-EigenMatrix3::Identity() + l0 / l_ij * ( EigenMatrix3::Identity() - ( x_ij * x_ij.transpose()) / ( l_ij * l_ij)));

	for( int row = 0; row < 3; ++row)
	{
		for( int col = 0; col < 3; ++col)
		{
			ScalarType val = k(row, col);
			jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p1 + row, 3 * m_p1 + col, val));
			jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p1 + row, 3 * m_p2 + col, -val));
			jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p2 + row, 3 * m_p1 + col, -val));
			jacobian_triplets.push_back(SparseMatrixTriplet(3 * m_p2 + row, 3 * m_p2 + col, val));
		}
	}
}

视频地址：https://youtu.be/HGtPOt5Ng7A