布料的模拟

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

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