 |
TAPs 0.7.7.3
|
|
TAPs 0.7.7.3
|
00001 /****************************************************************************** 00002 TAPsModelStrand.cpp 00003 ******************************************************************************/ 00007 /****************************************************************************** 00008 SUKITTI PUNAK (07/10/2005) 00009 UPDATE (09/03/2010) 00010 ******************************************************************************/ 00011 #include "TAPsModelStrand.hpp" 00012 // Using Inclusion Model (i.e. definitions are included in declarations) 00013 // (this name.cpp is included in name.hpp) 00014 // Each friend is defined directly inside its declaration. 00015 00016 BEGIN_NAMESPACE_TAPs__OpenGL 00017 //============================================================================= 00018 #ifdef TAPs_ADVANCED_SIMULATION 00019 template <typename T> unsigned int ModelStrand<T>::TotalModels = 0; 00020 #endif//TAPs_ADVANCED_SIMULATION 00021 //----------------------------------------------------------------------------- 00022 // For Assisting Constructor(s) 00023 template <typename T> 00024 void ModelStrand<T>::AllocateLinkProps ( Vector3<T> &posOfVertex0 ) 00025 { 00026 //---------------------------------------------------------------- 00027 // Create the the link direction (from 0 to # of vertices-1) 00028 if ( ( m_prLinkDirection = new Vector3<T>[ m_iNoLinks ] ) == NULL ) { 00029 #ifdef TAPs_ENABLE_DEBUG 00030 std::cerr << "ERROR => ModelStrand Constructor:" 00031 << " Cannot allocate memory for strand links!" 00032 << std::endl; 00033 #endif 00034 delete this; 00035 return; 00036 } 00037 00038 #ifdef TAPs_USE_EXTERNAL_ODE_SOLVER 00039 //---------------------------------------------------------------- 00040 // Create Links 00041 if ( ( m_prVertex = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00042 #ifdef TAPs_ENABLE_DEBUG 00043 std::cerr << "ERROR => ModelStrand Constructor:" 00044 << " Cannot allocate memory for strand links!" 00045 << std::endl; 00046 #endif 00047 delete this; 00048 return; 00049 } 00050 //---------------------------------------------------------------- 00051 // Create Velocities 00052 if ( ( m_prVelocity = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00053 #ifdef TAPs_ENABLE_DEBUG 00054 std::cerr << "ERROR => ModelStrand Constructor:" 00055 << " Cannot allocate memory for strand velocities!" 00056 << std::endl; 00057 #endif 00058 delete this; 00059 return; 00060 } 00061 //---------------------------------------------------------------- 00062 // Create Forces 00063 if ( ( m_prForce = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00064 #ifdef TAPs_ENABLE_DEBUG 00065 std::cerr << "ERROR => ModelStrand Constructor:" 00066 << " Cannot allocate memory for strand forces!" 00067 << std::endl; 00068 #endif 00069 delete this; 00070 return; 00071 } 00072 00073 #else //TAPs_USE_EXTERNAL_ODE_SOLVER 00074 00075 //---------------------------------------------------------------- 00076 // Prepare counters 00077 00078 // For Euler Integration 00079 #ifdef TAPs_USE_INTERNAL_ODE_SOLVER_EULER 00080 int numPositionSets = 2; 00081 int numVelocitySets = 2; 00082 int numForceSets = 2; 00083 m_PositionCounter.SetValueHighest( numPositionSets ); 00084 m_VelocityCounter.SetValueHighest( numVelocitySets ); 00085 m_ForceCounter.SetValueHighest( numForceSets ); 00086 #endif//TAPs_USE_INTERNAL_ODE_SOLVER_EULER 00087 00088 //---------------------------------------------------------------- 00089 00090 //---------------------------------------------------------------- 00091 // Create Links 00092 if ( ( m_prVertexPositionRecord = new Vector3<T>*[numPositionSets] ) == NULL ) { 00093 #ifdef TAPs_ENABLE_DEBUG 00094 std::cerr << "ERROR => ModelStrand Constructor:" 00095 << " Cannot allocate memory for strand links!" 00096 << std::endl; 00097 #endif 00098 delete this; 00099 return; 00100 } 00101 for ( int i = 0; i < numPositionSets; ++i ) { 00102 if ( ( m_prVertexPositionRecord[i] = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00103 #ifdef TAPs_ENABLE_DEBUG 00104 std::cerr << "ERROR => ModelStrand Constructor:" 00105 << " Cannot allocate memory for strand links!" 00106 << std::endl; 00107 #endif 00108 delete this; 00109 return; 00110 } 00111 } 00112 m_prVertex = m_prVertexPositionRecord[0]; 00113 m_prVertexPrev = m_prVertex; 00114 //---------------------------------------------------------------- 00115 // Create Velocities 00116 if ( ( m_prVertexVelocityRecord = new Vector3<T>*[numVelocitySets] ) == NULL ) { 00117 #ifdef TAPs_ENABLE_DEBUG 00118 std::cerr << "ERROR => ModelStrand Constructor:" 00119 << " Cannot allocate memory for strand velocitys!" 00120 << std::endl; 00121 #endif 00122 delete this; 00123 return; 00124 } 00125 for ( int i = 0; i < numVelocitySets; ++i ) { 00126 if ( ( m_prVertexVelocityRecord[i] = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00127 #ifdef TAPs_ENABLE_DEBUG 00128 std::cerr << "ERROR => ModelStrand Constructor:" 00129 << " Cannot allocate memory for strand velocitys!" 00130 << std::endl; 00131 #endif 00132 delete this; 00133 return; 00134 } 00135 } 00136 m_prVelocity = m_prVertexVelocityRecord[0]; 00137 m_prVelocityPrev = m_prVelocity; 00138 //---------------------------------------------------------------- 00139 // Create Forces 00140 if ( ( m_prVertexForceRecord = new Vector3<T>*[numForceSets] ) == NULL ) { 00141 #ifdef TAPs_ENABLE_DEBUG 00142 std::cerr << "ERROR => ModelStrand Constructor:" 00143 << " Cannot allocate memory for strand forces!" 00144 << std::endl; 00145 #endif 00146 delete this; 00147 return; 00148 } 00149 for ( int i = 0; i < numForceSets; ++i ) { 00150 if ( ( m_prVertexForceRecord[i] = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00151 #ifdef TAPs_ENABLE_DEBUG 00152 std::cerr << "ERROR => ModelStrand Constructor:" 00153 << " Cannot allocate memory for strand forces!" 00154 << std::endl; 00155 #endif 00156 delete this; 00157 return; 00158 } 00159 } 00160 m_prForce = m_prVertexForceRecord[0]; 00161 m_prForcePrev = m_prForce; 00162 //---------------------------------------------------------------- 00163 00164 #endif//TAPs_USE_EXTERNAL_ODE_SOLVER 00165 00166 //---------------------------------------------------------------- 00167 // Setup the start shape of the strand 00168 SetupShape ( posOfVertex0 ); 00169 /* 00170 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00171 posOfVertex0[0] += m_tLinkLength; 00172 m_prVertex[i] = posOfVertex0; 00173 } 00174 //*/ 00175 //---------------------------------------------------------------- 00176 // Create Links 00177 if ( ( m_prOrigVertex = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00178 #ifdef TAPs_ENABLE_DEBUG 00179 std::cerr << "ERROR => ModelStrand Constructor:" 00180 << " Cannot allocate memory for strand links!" 00181 << std::endl; 00182 #endif 00183 delete this; 00184 return; 00185 } 00186 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00187 m_prOrigVertex[i] = m_prVertex[i]; 00188 } 00189 //---------------------------------------------------------------- 00190 // Create Prev Links 00191 if ( ( m_prVertexPrevPos = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00192 #ifdef TAPs_ENABLE_DEBUG 00193 std::cerr << "ERROR => ModelStrand Constructor:" 00194 << " Cannot allocate memory for strand links!" 00195 << std::endl; 00196 #endif 00197 delete this; 00198 return; 00199 } 00200 //---------------------------------------------------------------- 00201 // Create Temp Vertices Set 1 00202 if ( ( m_prTempVertex1 = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00203 #ifdef TAPs_ENABLE_DEBUG 00204 std::cerr << "ERROR => ModelStrand Constructor:" 00205 << " Cannot allocate memory for strand links (temp set 1)!" 00206 << std::endl; 00207 #endif 00208 delete this; 00209 return; 00210 } 00211 //---------------------------------------------------------------- 00212 // Create Temp Vertices Set 2 00213 if ( ( m_prTempVertex2 = new Vector3<T>[ m_iNoLinks + 1 ] ) == NULL ) { 00214 #ifdef TAPs_ENABLE_DEBUG 00215 std::cerr << "ERROR => ModelStrand Constructor:" 00216 << " Cannot allocate memory for strand links (temp set 1)!" 00217 << std::endl; 00218 #endif 00219 delete this; 00220 return; 00221 } 00222 //---------------------------------------------------------------- 00223 // Create Rotations 00224 if ( ( m_prDegRot = new T[ m_iNoLinks + 1 ] ) == NULL ) { 00225 #ifdef TAPs_ENABLE_DEBUG 00226 std::cerr << "ERROR => ModelStrand Constructor:" 00227 << " Cannot allocate memory for strand rotations!" 00228 << std::endl; 00229 #endif 00230 delete this; 00231 return; 00232 } 00233 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00234 m_prDegRot[i] = T(); 00235 } 00236 //---------------------------------------------------------------- 00237 // Create Fix Status Flags 00238 if ( ( m_prIsFixed = new bool[ m_iNoLinks + 1 ] ) == NULL ) { 00239 #ifdef TAPs_ENABLE_DEBUG 00240 std::cerr << "ERROR => ModelStrand Constructor:" 00241 << " Cannot allocate memory for strand fix flags!" 00242 << std::endl; 00243 #endif 00244 delete this; 00245 return; 00246 } 00247 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00248 m_prIsFixed[i] = false; 00249 } 00250 00251 /* 00252 //---------------------------------------------------------------- 00253 // Create Was Intersect Status Flags 00254 if ( ( m_prWasIntersect = new char[ m_iNoLinks ] ) == NULL ) { 00255 #ifdef TAPs_ENABLE_DEBUG 00256 std::cerr << "ERROR => ModelStrand Constructor:" 00257 << " Cannot allocate memory for strand 'was intersect' flags!" 00258 << std::endl; 00259 #endif 00260 delete this; 00261 return; 00262 } 00263 for ( int i = 0; i < m_iNoLinks; ++i ) { 00264 m_prWasIntersect[i] = 0; 00265 } 00266 //*/ 00267 00268 //---------------------------------------------------------------- 00269 // Create (Subdivision) links 00270 if ( ( m_prDispVertex = new Vector3<T>[ m_iNoLinks * 2 ] ) == NULL ) { 00271 #ifdef TAPs_ENABLE_DEBUG 00272 std::cerr << "ERROR => ModelStrand Constructor:" 00273 << " Cannot allocate memory for strand (display) links!" 00274 << std::endl; 00275 #endif 00276 delete this; 00277 return; 00278 } 00279 SetDisplayLinks(); 00280 //---------------------------------------------------------------- 00281 // Create Structure Forces 00282 if ( ( m_StructForce = new StructForce[ m_iNoLinks + 1 ] ) == NULL ) { 00283 #ifdef TAPs_ENABLE_DEBUG 00284 std::cerr << "ERROR => ModelStrand Constructor:" 00285 << " Cannot allocate memory for strand forces!" 00286 << std::endl; 00287 #endif 00288 delete this; 00289 return; 00290 } 00291 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00292 m_StructForce[i].gravity = Vector3<T>( 0, -Physics_SI::K::g*m_tPointWeight, 0 ); 00293 //m_StructForce[i].antiBending; 00294 } 00295 //---------------------------------------------------------------- 00296 // For simulation 00297 #ifdef TAPs_USE_EXTERNAL_ODE_SOLVER 00298 if ( ( m_prODESolver = new Simulation::ODESolverEuler<T>() ) == NULL ) { 00299 //if ( ( m_prODESolver = new Simulation::ODESolverMidpoint<T>() ) == NULL ) { 00300 //if ( ( m_prODESolver = new Simulation::ODESolverRungeKutta4<T>() ) == NULL ) { 00301 #ifdef TAPs_ENABLE_DEBUG 00302 std::cerr << "ERROR => ModelStrand Constructor:" 00303 << " Cannot allocate memory for strand ODE solver!" 00304 << std::endl; 00305 #endif 00306 delete this; 00307 return; 00308 } 00309 m_iStateSize = (m_iNoLinks + 1) * 6; 00310 m_prODESolver->SetSize( m_iStateSize ); 00311 x0.resize( m_iStateSize ); 00312 xEnd.resize( m_iStateSize ); 00313 StateToArray( &xEnd[0] ); 00314 #endif//TAPs_USE_EXTERNAL_ODE_SOLVER 00315 00316 //---------------------------------------------------------------- 00317 InitCurrentShape(); 00318 00319 // Save Previous Positions for roll back position 00320 SavePreviousPositions(); 00321 00322 BuildBVHTree(); 00323 CheckSelfIntersections(); 00324 00325 #ifdef TAPs_DEBUG_MODE 00326 std::cout << "Sphere Bounding Hierarchy for Strand with " 00327 << m_BVHTree->CountNodes() << " nodes" 00328 << " which " << m_BVHTree->CountLeafNodes() 00329 << " of them are leaf nodes." << "\n"; 00330 #endif//TAPs_DEBUG_MODE 00331 //---------------------------------------------------------------- 00332 00333 //=============================================================== 00334 #ifdef TAPs_STRAND_LOCK_SEGMENTS 00335 //--------------------------------------------------------------- 00336 // Initialization for lock segments 00337 //------------------------------------------- 00338 // The strand starts as one segment 00339 m_Segments.push_back( IC_Segment( 0, m_iNoLinks, false, true ) ); 00340 //------------------------------------------- 00341 // Create Lock Status Flags 00342 if ( ( m_prIsLocked = new bool[ m_iNoLinks + 1 ] ) == NULL ) { 00343 #ifdef TAPs_ENABLE_DEBUG 00344 std::cerr << "ERROR => ModelStrand Constructor:" 00345 << " Cannot allocate memory for strand lock flags!" 00346 << std::endl; 00347 #endif 00348 delete this; 00349 return; 00350 } 00351 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00352 m_prIsLocked[i] = false; 00353 } 00354 //--------------------------------------------------------------- 00355 #endif//TAPs_STRAND_LOCK_SEGMENTS 00356 //=============================================================== 00357 00358 //=============================================================== 00359 #ifdef TAPs_STRAND_KNOT_CONTROL 00360 //--------------------------------------------------------------- 00361 CreateKnotControl(); 00362 //--------------------------------------------------------------- 00363 #endif//TAPs_STRAND_KNOT_CONTROL 00364 //=============================================================== 00365 00366 //=============================================================== 00367 #ifdef TAPs_ADVANCED_SIMULATION 00368 //--------------------------------------------------------------- 00369 // Create Simulation Flags 00370 if ( ( m_SimFlags = new DS::SimulationFlags[ m_iNoLinks + 1 ] ) == NULL ) { 00371 #ifdef TAPs_ENABLE_DEBUG 00372 std::cerr << "ERROR => ModelStrand Constructor:" 00373 << " Cannot allocate memory for strand simulation flags!" 00374 << std::endl; 00375 #endif 00376 delete this; 00377 return; 00378 } 00379 //--------------------------------------------------------------- 00380 #endif//TAPs_ADVANCED_SIMULATION 00381 //=============================================================== 00382 00383 } //EOF: AllocateLinkProps ( Vector3<T> &posOfVertex0 ) 00384 //----------------------------------------------------------------------------- 00385 // Default constructor 00386 template <typename T> 00387 ModelStrand<T>::ModelStrand () 00388 : OpenGLModel<T>(), 00389 m_iNoLinks ( 128 ), 00390 m_tRadius ( 0.05 ), 00391 m_tLinkLength ( 0.1 ), 00392 m_tPointWeight ( 0.1 ), 00393 m_tKRatioOfStretchingAllowed ( 0.10 ), 00394 m_tKRatioOfCompressionAllowed ( 0.05 ), 00395 m_tKStretching ( 2.0 ), 00396 m_tKBending ( 0.5 ), 00397 m_tKTwisting ( 0.5 ), 00398 m_tKFriction ( 0.5 ), 00399 m_tKContact1 ( 0.5 ), 00400 m_tKContact2 ( 0.9 ), 00401 m_tKGravity ( 9.81 ), 00402 m_tMoveStepDistAllowed ( m_tRadius * 0.5 ) 00403 #ifdef TAPs_STRAND_ENABLE_DAMPING_LIMIT 00404 , uiCounterDamping( 0 ) 00405 , uiLimitCounterDamping( 100 ) 00406 #endif//TAPs_STRAND_ENABLE_DAMPING_LIMIT 00407 #ifdef TAPs_ADVANCED_SIMULATION 00408 , AdvSimID( TotalModels++ ) 00409 #endif//TAPs_ADVANCED_SIMULATION 00410 { 00411 #ifdef TAPs_USE_CUDA 00412 // Data Members for CUDA 00413 m_cudaID = 0; 00414 m_cudaVertexList = NULL; 00415 m_cudaPrevVertexList = NULL; 00416 00417 #ifdef TAPs_ADVANCED_SIMULATION 00418 m_cudaSimFlagsList = NULL; 00419 m_cudaPosConstraintList = NULL; 00420 #endif//TAPs_ADVANCED_SIMULATION 00421 #endif//TAPs_USE_CUDA 00422 00423 //---------------------------------------------------------------- 00424 AllocateLinkProps( Vector3<T>(0,0,0) ); 00425 //---------------------------------------------------------------- 00426 #ifdef TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 00427 CreateDataSpaceForKnotRecognitionByDowkerNotation(); 00428 #endif//TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 00429 //---------------------------------------------------------------- 00430 #ifdef TAPs_ENABLE_DEBUG 00431 std::cout << "Default ModelStrand<" << typeid(T).name() << "> with:\n" 00432 << " " << m_iNoLinks << " links with each link of length " 00433 << m_tLinkLength << " and weight " << m_tPointWeight << "\n"; 00434 #endif//TAPs_DEBUG_MODE 00435 00436 #ifdef TAPs_USE_CUDA 00437 CUDA_Initialize_All(); 00438 #endif//TAPs_USE_CUDA 00439 } 00440 //----------------------------------------------------------------------------- 00441 // Constructor 00442 template <typename T> 00443 ModelStrand<T>::ModelStrand ( int iNoLinks, T tLength, T tWeight, 00444 Vector3<T> & posOfVertex0, T tRadius ) 00445 00446 : OpenGLModel<T>(), 00447 m_iNoLinks ( iNoLinks ), 00448 m_tRadius ( tRadius ), 00449 m_tLinkLength ( tLength / iNoLinks ), 00450 m_tPointWeight ( tWeight / iNoLinks ), 00451 m_tKRatioOfStretchingAllowed ( 0.10 ), 00452 m_tKRatioOfCompressionAllowed ( 0.05 ), 00453 m_tKStretching ( 10.0 ), 00454 m_tKBending ( 0.02 ), 00455 m_tKTwisting ( 0.5 ), 00456 m_tKFriction ( 5 ), 00457 m_tKGravity ( 9.81 ), 00458 m_tKContact1 ( 0.5 ), 00459 m_tKContact2 ( 0.9 ), 00460 m_tMoveStepDistAllowed ( m_tRadius * 0.5 ) 00461 #ifdef TAPs_STRAND_ENABLE_DAMPING_LIMIT 00462 , uiCounterDamping( 0 ) 00463 , uiLimitCounterDamping( 100 ) 00464 #endif//TAPs_STRAND_ENABLE_DAMPING_LIMIT 00465 #ifdef TAPs_ADVANCED_SIMULATION 00466 , AdvSimID( TotalModels++ ) 00467 #endif//TAPs_ADVANCED_SIMULATION 00468 { 00469 #ifdef TAPs_USE_CUDA 00470 // Data Members for CUDA 00471 m_cudaID = 0; 00472 m_cudaVertexList = NULL; 00473 m_cudaPrevVertexList = NULL; 00474 00475 #ifdef TAPs_ADVANCED_SIMULATION 00476 m_cudaSimFlagsList = NULL; 00477 m_cudaPosConstraintList = NULL; 00478 #endif//TAPs_ADVANCED_SIMULATION 00479 #endif//TAPs_USE_CUDA 00480 00481 //---------------------------------------------------------------- 00482 AllocateLinkProps( posOfVertex0 ); 00483 //---------------------------------------------------------------- 00484 #ifdef TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 00485 CreateDataSpaceForKnotRecognitionByDowkerNotation(); 00486 #endif//TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 00487 //---------------------------------------------------------------- 00488 #ifdef TAPs_ENABLE_DEBUG 00489 std::cout << "ModelStrand<" << typeid(T).name() << "> with:\n" 00490 << " " << m_iNoLinks << " links with each link of length " 00491 << m_tLinkLength << " and weight " << m_tPointWeight << "\n"; 00492 #endif//TAPs_DEBUG_MODE 00493 00494 #ifdef TAPs_USE_CUDA 00495 CUDA_Initialize_All(); 00496 #endif//TAPs_USE_CUDA 00497 } 00498 //----------------------------------------------------------------------------- 00499 // Destructor 00500 template <typename T> 00501 ModelStrand<T>::~ModelStrand () 00502 { 00503 //---------------------------------------------------------------- 00504 if ( m_prLinkDirection ) { 00505 delete [] m_prLinkDirection; 00506 m_prLinkDirection = NULL; 00507 } 00508 //---------------------------------------------------------- 00509 if ( m_prOrigVertex ) { 00510 delete [] m_prOrigVertex; 00511 m_prOrigVertex = NULL; 00512 } 00513 00514 #ifdef TAPs_USE_EXTERNAL_ODE_SOLVER 00515 if ( m_prVertex ) { 00516 delete [] m_prVertex; 00517 m_prVertex = NULL; 00518 } 00519 if ( m_prVelocity ) { 00520 delete [] m_prVelocity; 00521 m_prVelocity = NULL; 00522 } 00523 if ( m_prForce ) { 00524 delete [] m_prForce; 00525 m_prForce = NULL; 00526 } 00527 #else //TAPs_USE_EXTERNAL_ODE_SOLVER 00528 if ( m_prVertexPositionRecord ) { 00529 for ( int i = 0; i < m_PositionCounter.GetValueHighest(); ++i ) { 00530 delete [] m_prVertexPositionRecord[i]; 00531 m_prVertexPositionRecord[i] = NULL; 00532 } 00533 delete [] m_prVertexPositionRecord; 00534 } 00535 if ( m_prVertexVelocityRecord ) { 00536 for ( int i = 0; i < m_VelocityCounter.GetValueHighest(); ++i ) { 00537 delete [] m_prVertexVelocityRecord[i]; 00538 m_prVertexVelocityRecord[i] = NULL; 00539 } 00540 delete [] m_prVertexVelocityRecord; 00541 } 00542 if ( m_prVertexForceRecord ) { 00543 for ( int i = 0; i < m_ForceCounter.GetValueHighest(); ++i ) { 00544 delete [] m_prVertexForceRecord[i]; 00545 m_prVertexForceRecord[i] = NULL; 00546 } 00547 delete [] m_prVertexForceRecord; 00548 } 00549 #endif//TAPs_USE_EXTERNAL_ODE_SOLVER 00550 00551 //---------------------------------------------------------- 00552 if ( m_prVertexPrevPos ) { 00553 delete [] m_prVertexPrevPos; 00554 m_prVertexPrevPos = NULL; 00555 } 00556 if ( m_prTempVertex1 ) { 00557 delete [] m_prTempVertex1; 00558 m_prTempVertex1 = NULL; 00559 } 00560 if ( m_prTempVertex2 ) { 00561 delete [] m_prTempVertex2; 00562 m_prTempVertex2 = NULL; 00563 } 00564 if ( m_prDegRot ) { 00565 delete [] m_prDegRot; 00566 m_prDegRot = NULL; 00567 } 00568 if ( m_prIsFixed ) { 00569 delete [] m_prIsFixed; 00570 m_prIsFixed = NULL; 00571 } 00572 00573 /* 00574 if ( m_prWasIntersect ) { 00575 delete [] m_prWasIntersect; 00576 m_prWasIntersect = NULL; 00577 } 00578 //*/ 00579 00580 if ( m_prDispVertex ) { 00581 delete [] m_prDispVertex; 00582 m_prDispVertex = NULL; 00583 } 00584 if ( m_StructForce ) { 00585 delete [] m_StructForce; 00586 m_StructForce = NULL; 00587 } 00588 //---------------------------------------------------------------- 00589 #ifdef TAPs_ENABLE_DEBUG 00590 std::cout << "ModelStrand<" << typeid(T).name() << "> Destructor\n"; 00591 #endif//TAPs_DEBUG_MODE 00592 00593 //---------------------------------------------------------------- 00594 #ifdef TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 00595 DeleteDataSpaceForKnotRecognitionByDowkerNotation(); 00596 #endif//TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 00597 //---------------------------------------------------------------- 00598 00599 //---------------------------------------------------------------- 00600 // For simulation 00601 #ifdef TAPs_USE_EXTERNAL_ODE_SOLVER 00602 delete m_prODESolver; 00603 #endif//TAPs_USE_EXTERNAL_ODE_SOLVER 00604 delete m_BVHTree; 00605 //---------------------------------------------------------------- 00606 00607 //=============================================================== 00608 #ifdef TAPs_STRAND_KNOT_CONTROL 00609 //--------------------------------------------------------------- 00610 DeleteKnotControl(); 00611 //--------------------------------------------------------------- 00612 #endif//TAPs_STRAND_KNOT_CONTROL 00613 //=============================================================== 00614 00615 //=============================================================== 00616 #ifdef TAPs_USE_CUDA 00617 //--------------------------------------------------------------- 00618 CUDA_Cleanup_All(); 00619 //--------------------------------------------------------------- 00620 #endif//TAPs_USE_CUDA 00621 //=============================================================== 00622 00623 //=============================================================== 00624 #ifdef TAPs_ADVANCED_SIMULATION 00625 //--------------------------------------------------------------- 00626 // Create Simulation Flags 00627 if ( m_SimFlags ) { 00628 delete [] m_SimFlags; 00629 m_SimFlags = NULL; 00630 } 00631 //--------------------------------------------------------------- 00632 #endif//TAPs_ADVANCED_SIMULATION 00633 //=============================================================== 00634 } 00635 //============================================================================= 00636 // Get/Set Fn(s) 00637 //----------------------------------------------------------------------------- 00638 template <typename T> 00639 inline void ModelStrand<T>::SetFixStatusOfPtNo ( int i, bool b ) 00640 { 00641 if ( i < 0 || m_iNoLinks < i ) return; 00642 00643 if ( m_prIsFixed[i] == b ) return; // m_prIsFixed[i] is already b status 00644 if ( b ) { // add its # from the set of fixed points 00645 m_prIsFixed[i] = true; 00646 m_setOfFixedPts.insert( i ); 00647 #ifdef TAPs_ADVANCED_SIMULATION 00648 //m_SimFlags[i].SetSimulationConstraints( TAPs::Enum::AddOn::FIXED ); 00649 #endif//TAPs_ADVANCED_SIMULATION 00650 } 00651 else { // remove its # from the set of fixed points 00652 m_prIsFixed[i] = false; 00653 m_setOfFixedPts.erase( i ); 00654 #ifdef TAPs_ADVANCED_SIMULATION 00655 //m_SimFlags[i].SetSimulationConstraints( TAPs::Enum::AddOn::CLEARED ); 00656 #endif//TAPs_ADVANCED_SIMULATION 00657 } 00658 } 00659 //----------------------------------------------------------------------------- 00660 template <typename T> 00661 inline void ModelStrand<T>::ToggleFixStatusOfPtNo ( int i ) 00662 { 00663 SetFixStatusOfPtNo( i, !GetFixStatusOfPtNo(i) ); 00664 } 00666 template <typename T> 00667 void ModelStrand<T>::FixSegment ( int start, int end ) 00668 { 00669 if ( start < 0 ) start = 0; 00670 if ( end > m_iNoLinks ) end = m_iNoLinks; 00671 while ( start <= end ) { 00672 SetFixStatusOfPtNo( start, true ); 00673 ++start; 00674 } 00675 } 00676 //----------------------------------------------------------------------------- 00678 template <typename T> 00679 void ModelStrand<T>::UnfixSegment ( int start, int end ) 00680 { 00681 if ( start < 0 ) start = 0; 00682 if ( end > m_iNoLinks ) end = m_iNoLinks; 00683 while ( start <= end ) { 00684 SetFixStatusOfPtNo( start, false ); 00685 ++start; 00686 } 00687 } 00688 //----------------------------------------------------------------------------- 00689 template <typename T> 00690 inline Vector3<T> const & ModelStrand<T>::RetPointPosition ( int i ) const 00691 { 00692 return m_prVertex[i]; 00693 } 00694 //----------------------------------------------------------------------------- 00695 template <typename T> 00696 inline Vector3<T> & ModelStrand<T>::RetPointPosition ( int i ) 00697 { 00698 return m_prVertex[i]; 00699 } 00700 //----------------------------------------------------------------------------- 00701 template <typename T> 00702 inline Vector3<T> ModelStrand<T>::GetPointPosition ( int i ) const 00703 { 00704 return m_prVertex[i]; 00705 } 00706 //----------------------------------------------------------------------------- 00707 template <typename T> 00708 inline void ModelStrand<T>::SetPointPosition ( 00709 int i, const Vector3<T> & position, 00710 //T tMoveDistanceLimit, 00711 MultiBoundingVolume<T> const * const pMBV, 00712 TransformationSupport<T> const * const pTransform 00713 ) 00714 { 00715 #ifdef TAPs_STRAND_ENABLE_DAMPING_LIMIT 00716 // Reset counter for damping (for stopping vibration) 00717 uiCounterDamping = 0; 00718 #endif//TAPs_STRAND_ENABLE_DAMPING_LIMIT 00719 00720 // Save Previous Positions for roll back position 00721 //SavePreviousPositions(); 00722 00723 //-------------------------------------------------------------------- 00724 // Adjust the time 00725 static TAPs::Simulation::SimClock<Real> gClock( 0, 1E10, 0, 0.01 ); 00726 static Real prevTime = glutGet( GLUT_ELAPSED_TIME ); 00727 static Real currTime; 00728 00729 //SetFixStatusOfPtNo( i, true ); 00730 Vector3<T> u = position - m_prVertex[i]; 00731 00732 //* 00733 //m_tMoveStepDistAllowed = m_tRadius * 0.2; 00734 //T moveAllowed = m_tRadius * 200; 00735 T moveDistance = u.Length(); 00736 u = u.Normalized() * m_tMoveStepDistAllowed; 00737 T moveStep = moveDistance / m_tMoveStepDistAllowed; 00738 //u = u.Normalized() * tMoveDistanceLimit; 00739 //T moveStep = moveDistance / tMoveDistanceLimit; 00740 // T moveStep = 1; 00741 00742 // Limit number of steps 00743 const unsigned int LIMIT_STEPS = 250; 00744 unsigned int currentStep = 0; 00745 //---------------------------------------------------------------- 00746 //* 00747 00748 while ( moveStep > Math<T>::ZERO && ++currentStep < LIMIT_STEPS ) { 00749 00750 // Save Previous Positions for roll back position 00751 //SavePreviousPositions(); 00752 00753 if ( moveStep > Math<T>::ONE ) { 00754 moveStep -= Math<T>::ONE; 00755 } 00756 else { 00757 u *= moveStep; 00758 moveStep = Math<T>::ZERO; 00759 } 00760 //------------------------------------------------------ 00761 m_prVertex[i] += u; 00762 //-------------------------------------------------------------------- 00763 00764 //-------------------------------------------------------------------- 00765 //currTime = glutGet( GLUT_ELAPSED_TIME ); 00766 //gClock.step = (currTime - prevTime) / static_cast<Real>( 1000.0 ); 00767 //prevTime = currTime; 00768 //AdvanceSimulation( gClock.current, gClock.current + gClock.step ); 00769 AdvanceSimulation( 0.0, 0.005 ); 00770 } 00771 } 00772 //============================================================================= 00773 // Helper Fn(s) 00774 //----------------------------------------------------------------------------- 00775 template <typename T> 00776 void ModelStrand<T>::SetDisplayLinks () 00777 { 00778 m_prDispVertex[0] = m_prVertex[0]; 00779 for ( int i = 1, p = 0; i < m_iNoLinks; ++i ) { 00780 m_prDispVertex[++p] = (m_prVertex[i-1] + 3*m_prVertex[i]) / 4; 00781 m_prDispVertex[++p] = (3*m_prVertex[i] + m_prVertex[i+1]) / 4; 00782 } 00783 m_prDispVertex[m_iNoLinks*2 - 1] = m_prVertex[m_iNoLinks]; 00784 00785 //---------------------------------------------------------------- 00786 // Current positions ==> Previous positions 00787 // for ( int i = 0; i <= m_iNoLinks; ++i ) { 00788 // m_prVertexPrevPos[i] = m_prVertex[i]; 00789 // } 00790 } 00791 //============================================================================= 00792 // For Simulation 00793 //----------------------------------------------------------------------------- 00794 //----------------------------------------------------------------------------- 00795 // DxDt 00796 template <typename T> 00797 bool ModelStrand<T>::DxDt ( 00798 T dt, // i/p: time step 00799 Simulation::VectorSet<T> &x, // i/p: array x = {pos, vel} 00800 Simulation::VectorSet<T> &xdot, // o/p: array xdot = {vel, accel} 00801 void *userData // o/p: array of user data 00802 ) 00803 { 00804 //------------------------------------------------------ 00805 // Convert userData ptr to ModelStrand ptr 00806 ModelStrand<T> *pThis = static_cast<ModelStrand<T> *>( userData ); 00807 assert( pThis ); 00808 /* 00809 //------------------------------------------------------ 00810 for ( int i = 0; i <= pThis->m_iNoLinks; ++i ) { 00811 if ( pThis->m_prIsFixed[i] ) { 00812 pThis->m_prVelocity[i][0] = 00813 pThis->m_prVelocity[i][1] = 00814 pThis->m_prVelocity[i][2] = 00815 pThis->m_prForce[i][0] = 00816 pThis->m_prForce[i][1] = 00817 pThis->m_prForce[i][2] = 0; 00818 } 00819 else { 00820 //pThis->m_prForce[i] = pThis->m_StructForce[i].gravity; 00821 //pThis->m_prVelocity[i][1] = pThis->m_prForce[i][1] * 00822 // pThis->m_prVertex[i][1] * 00823 // dt; 00824 } 00825 } 00826 //*/ 00827 //------------------------------------------------------ 00828 // Copy data from x array into state variables 00829 //pThis->ArrayToState( &x[0] ); 00830 //------------------------------------------------------ 00831 pThis->DdtStateToArray( &xdot[0], dt ); 00832 //pThis->ClearForces(); 00833 //------------------------------------------------------ 00834 00835 /* 00836 for ( i = 0; i <= pThis->m_iNoLinks; ++i ) { 00837 //if ( pThis->m_prVertex[i][1] > 2 * pThis->m_tRadius ) { 00838 if ( false ) { 00839 // pThis->m_prForce[i] += pThis->CalForceGravity( i ); 00840 } 00841 else { 00842 // pThis->m_prForce[i][1] = 0; 00843 // pThis->m_prVelocity[i][1] = 0; 00844 pThis->m_prVertex[i][1] = 0; 00845 } 00846 } 00847 //*/ 00848 00849 return true; 00850 } 00851 //----------------------------------------------------------------------------- 00852 // StateToArray 00853 template <typename T> 00854 void ModelStrand<T>::StateToArray ( T *dst ) 00855 { 00856 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00857 *(dst++) = m_prVertex[i][0]; 00858 *(dst++) = m_prVertex[i][1]; 00859 *(dst++) = m_prVertex[i][2]; 00860 *(dst++) = m_prVelocity[i][0]; 00861 *(dst++) = m_prVelocity[i][1]; 00862 *(dst++) = m_prVelocity[i][2]; 00863 } 00864 } 00865 //----------------------------------------------------------------------------- 00866 // ArrayToState 00867 template <typename T> 00868 void ModelStrand<T>::ArrayToState ( T *src ) 00869 { 00870 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00871 m_prVertex[i][0] = *(src++); 00872 m_prVertex[i][1] = *(src++); 00873 m_prVertex[i][2] = *(src++); 00874 m_prVelocity[i][0] = *(src++); 00875 m_prVelocity[i][1] = *(src++); 00876 m_prVelocity[i][2] = *(src++); 00877 } 00878 } 00879 //----------------------------------------------------------------------------- 00880 // DdtStateToArray 00881 template <typename T> 00882 void ModelStrand<T>::DdtStateToArray ( T *xdot, T dt ) 00883 { 00884 //------------------------------------------- 00885 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00886 /*if ( m_prIsFixed[i] ) { 00887 *(xdot++) = 00888 *(xdot++) = 00889 *(xdot++) = 00890 *(xdot++) = 00891 *(xdot++) = 00892 *(xdot++) = 0; 00893 } 00894 else {*/ 00895 //m_prForce[i] = m_StructForce[i].gravity; 00896 //m_prVelocity[i][1] = m_prForce[i][1] * m_prVertex[i][1] * dt; 00897 //---------------------------------------- 00898 *(xdot++) = m_prVelocity[i][0]; 00899 *(xdot++) = m_prVelocity[i][1]; 00900 *(xdot++) = m_prVelocity[i][2]; 00901 *(xdot++) = m_prForce[i][0] / m_tPointWeight; 00902 *(xdot++) = m_prForce[i][1] / m_tPointWeight; 00903 *(xdot++) = m_prForce[i][2] / m_tPointWeight; 00904 //} 00905 } 00906 } 00907 //* 00908 //----------------------------------------------------------------------------- 00909 // ClearForces 00910 template <typename T> 00911 void ModelStrand<T>::ClearForces () 00912 { 00913 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00914 m_prForce[i].SetXYZ( 0, 0, 0 ); 00915 } 00916 } 00917 //----------------------------------------------------------------------------- 00918 // ClearVelocities 00919 template <typename T> 00920 void ModelStrand<T>::ClearVelocities () 00921 { 00922 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00923 m_prVelocity[i].SetXYZ( 0, 0, 0 ); 00924 } 00925 } 00926 //*/ 00927 //----------------------------------------------------------------------------- 00928 // SimGravity 00929 template <typename T> 00930 void ModelStrand<T>::SimGravity ( T tCurrent, T tNext ) 00931 { 00932 static int beingAtLinkNo = 2; 00933 T LIMIT = 0.1; 00934 T dt = tNext - tCurrent; // in seconds 00935 if ( dt > LIMIT ) dt = LIMIT; 00936 00937 //* 00938 int count = 0; 00939 //---------------------------------------------------------------- 00940 for ( int i = beingAtLinkNo; i <= m_iNoLinks; ++i ) { 00941 if ( m_prIsFixed[i] == false && 00942 m_prVertex[i][1] > 6 * m_tRadius ) { 00943 ++count; 00944 } 00945 } 00946 //---------------------------------------------------------------- 00947 if ( count > 0 ) { 00948 for ( int i = beingAtLinkNo; i <= m_iNoLinks; ++i ) { 00949 if ( m_prIsFixed[i] == false ) { 00950 if ( m_prVertex[i][1] > 3 * m_tRadius ) { 00951 //m_prVertex[i][1] -= Physics_SI::K::g*m_tPointWeight * (tNext - tCurrent) 00952 // / static_cast<T>( 1000 ); 00953 m_prVertex[i][1] -= Physics_SI::K::g*m_tPointWeight * dt; 00954 } 00955 if ( m_prVertex[i][1] < m_tRadius ) { 00956 m_prVertex[i][1] = m_tRadius; 00957 } 00958 } 00959 } 00960 } 00961 //*/ 00962 } 00963 //----------------------------------------------------------------------------- 00964 // AdvanceSimulation 00965 template <typename T> 00966 void ModelStrand<T>::AdvanceSimulation ( T tCurrent, T tNext ) 00967 { 00968 //std::cout << "FILE: " << __FILE__ << " LINE: " << __LINE__ << " -- AdvanceSimulation( T tCurrent, T tNext )\n"; 00969 00970 #ifdef TAPs_STRAND_DEBUG 00971 T currLen = GetCurrentLength(); 00972 T ideaLen = GetLength(); 00973 std::cout << "Strand Length:" << " (current) " << currLen << " (idea) " << ideaLen << "\n"; 00974 std::cout << " Len diff: " << currLen - ideaLen << "\n"; 00975 #endif//TAPs_STRAND_DEBUG 00976 00977 #ifdef TAPs_STRAND_ENABLE_DAMPING_LIMIT 00978 if ( uiLimitCounterDamping < ++uiCounterDamping ) return; 00979 #endif//TAPs_STRAND_ENABLE_DAMPING_LIMIT 00980 00981 //int steps = 1; 00982 //T currentTime = tCurrent; 00983 //T timeStep = tNext / steps; 00984 //T nextTime = tCurrent + timeStep; 00985 //for ( int st = 0; st < steps; ++st ) 00986 { 00987 //------------------------------------------------------ 00988 //* 00989 //------------------------------------------------------ 00990 for ( int i = 0; i <= m_iNoLinks; ++i ) { 00991 if ( m_prIsFixed[i] ) { 00992 m_prVelocity[i][0] = 00993 m_prVelocity[i][1] = 00994 m_prVelocity[i][2] = 00995 m_prForce[i][0] = 00996 m_prForce[i][1] = 00997 m_prForce[i][2] = 0; 00998 } 00999 else { 01000 //pThis->m_prForce[i] = pThis->m_StructForce[i].gravity; 01001 //pThis->m_prVelocity[i][1] = pThis->m_prForce[i][1] * 01002 // pThis->m_prVertex[i][1] * 01003 // dt; 01004 m_prForce[i].SetXYZ( 0, 0, 0 ); 01005 } 01006 } 01007 //*/ 01008 01009 01010 //* 01011 // This is put in AdvanceSimulation Fn 01012 //------------------------------------------------------------------ 01013 // Reset Forces and Velocities 01014 // for ( int i = 0; i <= m_iNoLinks; ++i ) { 01015 // m_prForce[i].SetXYZ( 0, 0, 0 ); 01016 // //m_prVelocity[i].SetXYZ( 0, 0, 0 ); 01017 // } 01018 //------------------------------------------------------ 01019 // UpdateBVHTree(); 01020 // CheckSelfIntersections(); 01021 //------------------------------------------------------ 01022 //*/ 01023 01024 //=============================================================== 01025 //--------------------------------------------------------------- 01026 #ifdef TAPs_USE_EXTERNAL_ODE_SOLVER 01027 //--------------------------------------------------------------- 01028 01029 //------------------------------------------------- 01030 StateToArray( &x0[0] ); 01031 SimByForce( tCurrent, tNext ); 01032 // ODE Solver 01033 m_prODESolver->Solve( x0, xEnd, tCurrent, tNext, DxDt, this ); 01034 // Copy d/dt X(tNext) into state variables xEnd 01035 ArrayToState( &xEnd[0] ); 01036 //------------------------------------------------- 01037 01038 //--------------------------------------------------------------- 01039 #else //TAPs_USE_EXTERNAL_ODE_SOLVER 01040 //--------------------------------------------------------------- 01041 01042 01043 //------------------------------------------------- 01044 // Sim by Verlet Intetration using CUDA 01045 //------------------------------------------------- 01046 #ifdef TAPs_USE_CUDA 01047 01048 //PrintDebug_ConnectionList(); 01049 //PrintDebug(); 01050 01051 // Copy the (current) vertex positions of the object to memory for transfering to CUDA. 01052 CUDA_CopyVertexToMem(); 01053 #ifdef TAPs_ADVANCED_SIMULATION 01054 CUDA_CopySimFlagsToMem(); 01055 CUDA_CopyPosConstraintToMem(); 01056 #endif//TAPs_ADVANCED_SIMULATION 01057 //CUDA_CopyPrevVertexToMem(); 01058 01059 // DEBUG 01060 for ( int i = 0; i <= m_iNoLinks; i+=m_iNoLinks ) { 01061 int idx = i*4; 01062 printf( "m_cudaVertexList# %d : %g %g %g %g\n", i, m_cudaVertexList[idx], m_cudaVertexList[idx+1], m_cudaVertexList[idx+2], m_cudaVertexList[idx+3] ); 01063 printf( "m_cudaPrevVertexList# %d: %g %g %g %g\n", i, m_cudaPrevVertexList[idx], m_cudaPrevVertexList[idx+1], m_cudaPrevVertexList[idx+2], m_cudaPrevVertexList[idx+3] ); 01064 } 01065 01066 // Call the wrapper function for 01067 #ifdef TAPs_ADVANCED_SIMULATION 01068 TAPs::CUDA::Global__ModelStrand_AdvSim_ADVSIM( 01069 m_cudaID, 01070 m_iNoLinks+1, 01071 64, 01072 tCurrent, 01073 0.001*10, //GetPredefinedTimeStep(), //!< time step 01074 10, //GetNumSimSubSteps(), //!< number of sub-steps 01075 m_tPointWeight, 01076 GetKStretching(), 01077 GetKFriction(), 01078 GetLinkLength(), 01079 m_cudaVertexList, 01080 m_cudaPrevVertexList, 01081 m_cudaSimFlagsList, 01082 m_cudaPosConstraintList 01083 ); 01084 #else //TAPs_ADVANCED_SIMULATION 01085 TAPs::CUDA::Global__ModelStrand_AdvSim( 01086 m_cudaID, 01087 m_iNoLinks+1, 01088 64, 01089 tCurrent, 01090 0.001*10, //GetPredefinedTimeStep(), //!< time step 01091 10, //GetNumSimSubSteps(), //!< number of sub-steps 01092 m_tPointWeight, 01093 GetKStretching(), 01094 GetKFriction(), 01095 GetLinkLength(), 01096 m_cudaVertexList, 01097 m_cudaPrevVertexList 01098 //m_cudaHomeVertexList, 01099 //m_cudaVertexConnectionList, 01100 ); 01101 #endif//TAPs_ADVANCED_SIMULATION 01102 01103 // The new vertex positions are returned by CUDA to memory for previous vertex positions. 01104 // So the new vertex positions must be copied to the vertex positions of the object. 01105 CUDA_CopyMemToPrevVertex(); 01106 // Then the previous and current vertex pointers (that point to memory for transfering to CUDA) are swapped. 01107 // So that now the previous becomes current and the current becomes previous. 01108 CUDA_SwapBuffers(); 01109 01110 #ifdef TAPs_ADVANCED_SIMULATION 01111 // The new strand's simulation flags list is returned by CUDA. 01112 // It must be copied to the simulation flags list of the strand. 01113 CUDA_CopyMemToSimFlags(); 01114 CUDA_CopyMemToPosConstraint(); 01115 01116 // NEED TO CHANGE THE VALUES STORED IN m_pAdvSimData 01117 01118 #endif//TAPs_ADVANCED_SIMULATION 01119 01120 01121 GetBVHTree()->Update(); // commented out to rely on CDR to update the collided BVNodes 01122 01123 #else //TAPs_USE_CUDA 01124 //------------------------------------------------- 01125 // Sim by Verlet Intetration using CUDA 01126 //------------------------------------------------- 01127 01128 01129 //------------------------------------------------- 01130 // Sim by Euler Intetration 01131 //------------------------------------------------- 01132 #ifdef TAPs_USE_INTERNAL_ODE_SOLVER_EULER 01133 // Remark: here ' represents next value, not differentiation 01134 // a' = f/m 01135 // v' = v + a'*t = v + f*t/m 01136 // x' = x + v'*t 01137 01138 T t = tNext - tCurrent; 01139 //T t = 0.005; 01140 01141 // ERROR WHY??? 01142 //CycleBuffer(); // make next (a',v',x') become current and current become previous 01143 01144 //--------------------------------------- 01145 // Calculate F 01146 // For collision we use geometric constraints instead of (external) force constraints 01147 SimByForce( tCurrent, tNext ); // Calculate F_internal 01148 01149 //--------------------------------------- 01150 // Calculate V' and X' 01151 #ifdef TAPs_STRAND_LOCK_SEGMENTS 01152 // For segment# 0 01153 { 01154 unsigned int s = 0; 01155 if ( !m_Segments[s].isLocked ) { 01156 for ( unsigned int i = m_Segments[s].start; i <= m_Segments[s].end; ++i ) 01157 { 01158 if ( !m_prIsFixed[i] ) { 01159 m_prVelocity[i] = m_prVelocityPrev[i] + m_prForcePrev[i]*t/m_tPointWeight; 01160 m_prVertex[i] = m_prVertexPrev[i] + m_prVelocityPrev[i]*t; 01161 } 01162 } 01163 } 01164 } 01165 // For segment# > 0 01166 for ( unsigned int s = 1; s < m_Segments.size(); ++s ) { 01167 if ( !m_Segments[s].isLocked ) { 01168 for ( unsigned int i = m_Segments[s].start; i <= m_Segments[s].end; ++i ) 01169 #else //TAPs_STRAND_LOCK_SEGMENTS 01170 { 01171 { 01172 for ( int i = 0; i < m_iNoLinks+1; ++i ) 01173 #endif//TAPs_STRAND_LOCK_SEGMENTS 01174 { 01175 if ( !m_prIsFixed[i] ) { 01176 m_prVelocity[i] = m_prVelocityPrev[i] + m_prForcePrev[i]*t/m_tPointWeight; 01177 m_prVertex[i] = m_prVertexPrev[i] + m_prVelocityPrev[i]*t; 01178 } 01179 } 01180 } 01181 } 01182 01183 #endif//TAPs_USE_INTERNAL_ODE_SOLVER_EULER 01184 //------------------------------------------------- 01185 01186 #endif //TAPs_USE_CUDA 01187 01188 //--------------------------------------------------------------- 01189 #endif//TAPs_USE_EXTERNAL_ODE_SOLVER 01190 //--------------------------------------------------------------- 01191 //=============================================================== 01192 01193 UpdateBVHTree(); 01194 CheckSelfIntersections(); 01195 //------------------------------------------------------ 01196 } 01197 } 01198 //----------------------------------------------------------------------------- 01199 // AdvanceSimulation 01200 template <typename T> 01201 void ModelStrand<T>::AdvanceSimulation ( Simulation::SimClock<T> & simClock ) 01202 { 01203 std::cout << "FILE: " << __FILE__ << " LINE: " << __LINE__ << " -- Not Implemented Yet!\n"; 01204 } 01205 //----------------------------------------------------------------------------- 01206 01207 01208 //============================================================================= 01209 // Collision Detection 01210 //----------------------------------------------------------------------------- 01211 template <typename T> 01212 bool ModelStrand<T>::TestOverlapWith ( BVHTree<T> const * const that ) 01213 { 01214 m_svCollidedNode.clear(); 01215 m_svCollidedNodeThat.clear(); 01216 //-------------------------------------------------------------------- 01217 bool result = GetBVHTree()->TestOverlapWithTillLeafNodes( that ); 01218 if ( !result ) return result; 01219 //-------------------------------------------------------------------- 01220 BVHNode<T> * cylinder; 01221 BVHNode<T> * triangle; 01222 //-------------------------------------------------------------------- 01223 int size = static_cast<int>( GetBVHTree()->GetListOfCollidedNodes().size() ); 01224 int linkNo; 01225 for ( int i = 0; i < size; ++i ) { 01226 cylinder = GetBVHTree()->GetListOfCollidedNodes()[i]; 01227 triangle = GetBVHTree()->GetListOfCollidedNodesThat()[i]; 01228 linkNo = cylinder->GetID(); 01229 HEFace<T> * heFace = triangle->GetAPrimitiveHalfEdgeFace(); 01230 std::vector< HEVertex<T> * const > vertices = heFace->GetPtrsToVertices(); 01231 result = CGMath<T>::FindIntersectionCylinderTriangle( 01232 GetPointPosition( linkNo ), GetPointPosition( linkNo+1 ), GetRadius(), 01233 vertices[0]->GetPosition(), 01234 vertices[1]->GetPosition(), 01235 vertices[2]->GetPosition() ); 01236 01237 if ( result ) { 01238 m_svCollidedNode.push_back( GetBVHTree()->GetListOfCollidedNodes()[i] ); 01239 m_svCollidedNodeThat.push_back( GetBVHTree()->GetListOfCollidedNodesThat()[i] ); 01240 } 01241 } 01242 //-------------------------------------------------------------------- 01243 if ( static_cast<int>( m_svCollidedNode.size() ) > 0 ) return true; 01244 else return false; 01245 } 01246 01247 //----------------------------------------------------------------------------- 01248 // BuildBVHTree 01249 template <typename T> 01250 void ModelStrand<T>::BuildBVHTree () 01251 { 01252 m_BVHTree = NULL; 01253 int numOfNodes = -1; 01254 //---------------------------------------------------------------- 01255 Vector3<T> center; 01256 01257 #ifdef TAPs_DEBUG_MODE 01258 std::cout << "m_tRadius: " << m_tRadius << std::endl; 01259 std::cout << "m_tLinkLength: " << m_tLinkLength << std::endl; 01260 #endif//TAPs_DEBUG_MODE 01261 01262 T radius = m_tRadius + m_tLinkLength/2.0; 01263 01264 #ifdef TAPs_DEBUG_MODE 01265 std::cout << "radius: " << radius << std::endl; 01266 #endif//TAPs_DEBUG_MODE 01267 01268 //---------------------------------------------------------------- 01269 // Create Leaf Nodes 01270 BVHNode<T> **nodeList = new BVHNode<T> *[m_iNoLinks]; 01271 for ( int i = 0; i < m_iNoLinks; ++i ) { 01272 //nodeList[i] = new BVHNode<T>( TAPs::Enum::BVH_NODE_BINARY_SPHERE, ++numOfNodes, NULL, NULL ); 01273 //nodeList[i] = new BVHNode<T>( TAPs::Enum::BVH_NODE_BINARY_SPHERE, ++numOfNodes, NULL, NULL ); 01274 nodeList[i] = new BVHNodeLeafHalfEdgeAPrim<T>( 01275 TAPs::Enum::BVH_NODE_LEAF_SPHERE_HALFEDGE_APRIM, 01276 ++numOfNodes, NULL ); 01277 // nodeList[i] = new BVHNodeLeafHalfEdgePrims<T>( 01278 // TAPs::Enum::BVH_NODE_LEAF_SPHERE_HALFEDGE_PRIMS, 01279 // ++numOfNodes, NULL ); 01280 01281 nodeList[i]->SetRadius( radius ); 01282 nodeList[i]->SetCenter( (m_prVertex[i] + m_prVertex[i+1]) / 2.0 ); 01283 } 01284 //---------------------------------------------------------------- 01285 // Create Hierarchy Nodes 01286 m_BVHTree = BuildBVHTreeRecursively( numOfNodes+1, nodeList, m_iNoLinks ); 01287 01288 #ifdef TAPs_DEBUG_MODE 01289 std::cout << *m_BVHTree << std::endl; 01290 #endif//TAPs_DEBUG_MODE 01291 } 01292 //----------------------------------------------------------------------------- 01293 // BuildBVHTreeRecursively 01294 template <typename T> 01295 BVHTree<T> * ModelStrand<T>::BuildBVHTreeRecursively ( 01296 int startID, 01297 BVHNode<T> ** childNodeList, 01298 int numOfChildNodes ) 01299 { 01300 BVHNode<T> ** parentNodeList; 01301 if ( numOfChildNodes > 1 ) { 01302 int numOfParentNodes; 01303 if ( numOfChildNodes % 2 == 0 ) { 01304 numOfParentNodes = numOfChildNodes/2; 01305 parentNodeList = new BVHNode<T> *[ numOfParentNodes ]; 01306 for ( int i = 0, p = 0; i < numOfChildNodes; i+=2, ++p, ++startID ) { 01307 parentNodeList[p] = new BVHNode<T>( 01308 TAPs::Enum::BVH_NODE_BINARY_SPHERE, 01309 startID, NULL, NULL ); 01310 parentNodeList[p]->Child( 0, childNodeList[i] ); 01311 parentNodeList[p]->Child( 1, childNodeList[i+1] ); 01312 } 01313 } 01314 else { 01315 numOfParentNodes = numOfChildNodes/2 + 1; 01316 parentNodeList = new BVHNode<T> *[ numOfParentNodes ]; 01317 for ( int i = 0, p = 0; i < numOfChildNodes-1; i+=2, ++p, ++startID ) { 01318 parentNodeList[p] = new BVHNode<T>( 01319 TAPs::Enum::BVH_NODE_BINARY_SPHERE, 01320 startID, NULL, NULL ); 01321 parentNodeList[p]->Child( 0, childNodeList[i] ); 01322 parentNodeList[p]->Child( 1, childNodeList[i+1] ); 01323 } 01324 parentNodeList[numOfParentNodes-1] = childNodeList[numOfChildNodes-1]; 01325 } 01326 delete [] childNodeList; 01327 return BuildBVHTreeRecursively( startID, parentNodeList, numOfParentNodes ); 01328 } 01329 else { 01330 BVHTree<T> * tree = new BVHTree_BinarySphere<T>( GetTransform(), 01331 /*TAPs::Enum::BVH_TREE_BINARY_SPHERE,*/ startID, childNodeList[0] ); 01332 delete [] childNodeList; 01333 return tree; 01334 } 01335 } 01336 //----------------------------------------------------------------------------- 01337 // UpdateBVHTree 01338 template <typename T> 01339 void ModelStrand<T>::UpdateBVHTree () 01340 { 01341 if ( m_BVHTree->Root() ) { 01342 UpdateBVHNodeRecursively( m_BVHTree->Root() ); 01343 } 01344 } 01345 //----------------------------------------------------------------------------- 01346 // UpdateBVHNodeRecursively 01347 template <typename T> 01348 void ModelStrand<T>::UpdateBVHNodeRecursively ( BVHNode<T> * node ) 01349 { 01350 if ( node->Child(0) ) { 01351 UpdateBVHNodeRecursively( node->Child(0) ); 01352 } 01353 if ( node->Child(1) ) { 01354 UpdateBVHNodeRecursively( node->Child(1) ); 01355 } 01356 //---------------------------------------------------------------- 01357 if ( node->Child(0) == NULL && node->Child(1) == NULL ) { 01358 int id = node->GetID(); 01359 node->SetRadius( (m_prVertex[id] - m_prVertex[id+1]).Length()/2.0 + m_tRadius ); 01360 node->SetCenter( (m_prVertex[id] + m_prVertex[id+1]) / 2.0 ); 01361 } 01362 else { 01363 node->Update(); 01364 } 01365 } 01366 //----------------------------------------------------------------------------- 01367 // CheckSelfIntersections 01368 // Assume no intersections between immediate adjacent links 01369 // Assume a BVH node either has both children or none, therefore 01370 // checking either it has a left or right child is enough to determine 01371 // whether it is a leaf node or not. 01372 template <typename T> 01373 void ModelStrand<T>::CheckSelfIntersections () 01374 { 01375 #ifdef TAPs_CURRENT_DEBUG 01376 g_iDebugCounter01 = 0; 01377 #endif 01378 //---------------------------------------------------------------- 01379 //std::cout << "\n\nCheckSelfIntersections Result:\n"; 01380 //if ( !m_BVHTree->Root() ) return; 01381 //BVHNode<T> * leftNode = m_BVHTree->Root()->Child(0); 01382 //BVHNode<T> * rightNode = m_BVHTree->Root()->Child(1); 01383 //---------------------------------------------------------------- 01384 //CheckSelfIntersectionsNextLevelRecursively( leftNode, rightNode ); 01385 01386 /* 01387 for ( int i = 0; i < m_iNoLinks; ++i ) { 01388 if ( m_prWasIntersect[i] > 0 ) { 01389 --m_prWasIntersect[i]; 01390 } 01391 } 01392 //*/ 01393 01394 //---------------------------------------------------------------- 01395 CheckSelfIntersectionsNextLevelRecursively( m_BVHTree->Root()->Child(0), 01396 m_BVHTree->Root()->Child(1) ); 01397 //---------------------------------------------------------------- 01398 01399 #ifdef TAPs_CURRENT_DEBUG 01400 if ( g_iDebugCounter01 > 10 ) { 01401 sprintf( g_cstrDebugArray[32], "Possibly a knot is formed!!!" ); 01402 } 01403 else { 01404 sprintf( g_cstrDebugArray[32], "No knots" ); 01405 } 01406 // Clear the less of the string rows 01407 for ( ; g_iDebugCounter01 < 32; ++g_iDebugCounter01 ) { 01408 sprintf( g_cstrDebugArray[g_iDebugCounter01], "" ); 01409 } 01410 #endif 01411 } // EOF: CheckSelfIntersections () 01412 //----------------------------------------------------------------------------- 01413 // CheckSelfIntersectionsNextLevelRecursively 01414 // Assume no intersections between immediate adjacent links 01415 // Assume a BVH node either has both children or none, therefore 01416 // checking either it has a left or right child is enough to determine 01417 // whether it is a leaf node or not. 01418 template <typename T> 01419 void ModelStrand<T>::CheckSelfIntersectionsNextLevelRecursively ( 01420 BVHNode<T> * node1, 01421 BVHNode<T> * node2 ) 01422 { 01423 CheckSelfIntersectionsRecursively( node1, node2 ); 01424 if ( node1->Child(0) ) 01425 CheckSelfIntersectionsNextLevelRecursively( node1->Child(0), node1->Child(1) ); 01426 if ( node2->Child(0) ) 01427 CheckSelfIntersectionsNextLevelRecursively( node2->Child(0), node2->Child(1) ); 01428 } 01429 //----------------------------------------------------------------------------- 01430 // CheckSelfIntersectionsRecursively 01431 // Assume no intersections between immediate adjacent links 01432 // Use the fact that the ids of sphere bvh leaf nodes are the same as 01433 // the link id and 01434 template <typename T> 01435 T ModelStrand<T>::CheckSelfIntersectionsRecursively ( 01436 BVHNode<T> * node1, 01437 BVHNode<T> * node2 ) 01438 { 01439 static T LinkDiameter = static_cast<T>(2.0) * m_tRadius; 01440 static T epsilon = m_tRadius * static_cast<T>(0.05); // for collision response 01441 //---------------------------------------------------------------- 01442 //if ( !node1 || !node2 ) return Math<T>::MAX; 01443 if ( !node1->Child(0) && !node2->Child(0) && 01444 ( abs(node1->GetID() - node2->GetID()) == 1 ) ) 01445 { 01446 return DBL_MAX; 01447 } 01448 //---------------------------------------------------------------- 01449 T overlap = node1->TestOverlapSphereWithSphere_NoPrimitiveTests( node2 ); 01450 //std::cout << "Node#" << node1->GetID() << " intersects with " 01451 // << "Node#" << node2->GetID() << " by " << overlap << "\n"; 01452 //---------------------------------------------------------------- 01453 if ( overlap < DBL_MIN ) { 01454 if ( node2->Child(0) == NULL ) { //&& node2->RightChild() == NULL ) { 01455 if ( node1->Child(0) ) { 01456 overlap = CheckSelfIntersectionsRecursively( node1->Child(0), node2 ); 01457 } 01458 if ( node1->Child(1) ) { 01459 overlap = CheckSelfIntersectionsRecursively( node1->Child(1), node2 ); 01460 } 01461 } 01462 else { 01463 if ( node2->Child(0) ) { 01464 overlap = node1->TestOverlapSphereWithSphere_NoPrimitiveTests( node2->Child(0) ); 01465 if ( overlap < DBL_MIN ) { 01466 if ( node1->Child(0) ) { 01467 overlap = CheckSelfIntersectionsRecursively( node1->Child(0), node2->Child(0) ); 01468 } 01469 if ( node1->Child(1) ) { 01470 overlap = CheckSelfIntersectionsRecursively( node1->Child(1), node2->Child(0) ); 01471 } 01472 } 01473 } 01474 if ( node2->Child(1) ) { 01475 overlap = node1->TestOverlapSphereWithSphere_NoPrimitiveTests( node2->Child(1) ); 01476 if ( overlap < DBL_MIN ) { 01477 if ( node1->Child(0) ) { 01478 overlap = CheckSelfIntersectionsRecursively( node1->Child(0), node2->Child(1) ); 01479 } 01480 if ( node1->Child(1) ) { 01481 overlap = CheckSelfIntersectionsRecursively( node1->Child(1), node2->Child(1) ); 01482 } 01483 } 01484 } 01485 } 01486 //------------------------------------------------------------ 01487 if ( node1->Child(0) == NULL && //node1->Child(1) == NULL && 01488 node2->Child(0) == NULL ) //&& node2->Child(1) == NULL ) 01489 { 01490 T distance; 01491 Vector3<T> U; 01492 int i = node1->GetID(); 01493 int j = node2->GetID(); 01494 01495 /* 01496 //-------------------------------------------------- 01497 // For Knot Tying (in intersection) 01498 if ( m_iLinkKnotGrp[i] >= 0 && (m_iLinkKnotGrp[i] == m_iLinkKnotGrp[j]) ) { 01499 return 0; 01500 } 01501 //*/ 01502 01503 //-------------------------------------------------- 01504 CGMath<T>::FindClosestDistanceBetweenTwoLineSegments( 01505 m_prVertex[i], m_prVertex[i+1], // Link# i 01506 m_prVertex[j], m_prVertex[j+1], // Link# j 01507 distance, U 01508 ); 01509 //* 01510 //-------------------------------------------------- 01511 // If collide, push the nodes apart GEOMETRICALLY 01512 if ( distance < LinkDiameter ) { 01513 //m_prVertex 01514 //std::cout << "Node#" << node1->GetID() << " intersects with " 01515 // << "Node#" << node2->GetID() << " by " << overlap; 01516 //std::cout << " with CONTACT: " << distance << " < " << LinkDiameter << "\n"; 01517 U *= m_tRadius - distance/static_cast<T>(2.0) + epsilon; 01518 01519 /* 01520 T length = U.Length(); 01521 if ( length > m_tRadius*4 ) { 01522 U /= length * m_tRadius*4; 01523 } 01524 //*/ 01525 01526 //U /= 2; 01527 /* 01528 m_prVertex[i] -= U; 01529 m_prVertex[i+1] -= U; 01530 m_prVertex[j] += U; 01531 m_prVertex[j+1] += U; 01532 //*/ 01533 01534 01535 #ifdef TAPs_STRAND_LOCK_SEGMENTS 01536 { 01537 //if ( !m_prIsLocked[i] && !m_prIsFixed[i] ) m_prVertex[i] -= U; 01538 //if ( !m_prIsLocked[j] && !m_prIsFixed[j] ) m_prVertex[j] += U; 01539 //if ( !m_prIsLocked[i+1] && !m_prIsFixed[i+1] ) m_prVertex[i+1] -= U; 01540 //if ( !m_prIsLocked[j+1] && !m_prIsFixed[j+1] ) m_prVertex[j+1] += U; 01541 if ( !m_prIsLocked[i] && !m_prIsLocked[i+1] ) { 01542 /* 01543 // Limit the link lneght, to prevent the non-stop move of point(s). 01544 // Here the length is limited to be not over a length threshold. 01545 T length = ( m_prVertex[i] - m_prVertex[i+1] ).Length(); 01546 if ( length < 2 * m_tLinkLength ) { 01547 if ( !m_prIsFixed[i] ) m_prVertex[i] -= U; 01548 if ( !m_prIsFixed[i+1] ) m_prVertex[i+1] -= U; 01549 } 01550 //*/ 01551 if ( !m_prIsFixed[i] && !m_prIsFixed[i+1] ) { 01552 m_prVertex[i] -= U; 01553 m_prVertex[i+1] -= U; 01554 } 01555 } 01556 if ( !m_prIsLocked[j] && !m_prIsLocked[j+1] ) { 01557 /* 01558 // Limit the link lneght, to prevent the non-stop move of point(s). 01559 // Here the length is limited to be not over a length threshold. 01560 T length = ( m_prVertex[j] - m_prVertex[j+1] ).Length(); 01561 if ( length < 2 * m_tLinkLength ) { 01562 if ( !m_prIsFixed[j] ) m_prVertex[j] += U; 01563 if ( !m_prIsFixed[j+1] ) m_prVertex[j+1] += U; 01564 } 01565 //*/ 01566 if ( !m_prIsFixed[j] && !m_prIsFixed[j+1] ) { 01567 m_prVertex[j] += U; 01568 m_prVertex[j+1] += U; 01569 } 01570 } 01571 } 01572 #else //TAPs_STRAND_LOCK_SEGMENTS 01573 { 01574 if ( !m_prIsFixed[i] ) m_prVertex[i] -= U; 01575 if ( !m_prIsFixed[j] ) m_prVertex[j] += U; 01576 if ( !m_prIsFixed[i+1] ) m_prVertex[i+1] -= U; 01577 if ( !m_prIsFixed[j+1] ) m_prVertex[j+1] += U; 01578 } 01579 #endif//TAPs_STRAND_LOCK_SEGMENTS 01580 01581 //m_prVertex[i] -= U; 01582 //m_prVertex[j] += U; 01583 //m_prVertex[i+1] -= U; 01584 //m_prVertex[j+1] += U; 01585 } 01586 01587 //*/ 01588 //-------------------------------------------------- 01589 // If collide, push the nodes apart PHYSICALLY (BASED ON FORCE) 01590 //Vector3<T> dir = ( node2->GetCenter() - node1->GetCenter() ); 01591 //distance = ( node1->GetRadius() + node2->GetRadius() ) - dir.Length(); 01592 //if ( distance > 0 ) { 01593 //* 01594 01595 /* 01596 //distance *= 2; 01597 if ( distance < LinkDiameter ) { 01598 Vector3<T> dir = ( node2->GetCenter() - node1->GetCenter() ); 01599 dir.Normalized(); 01600 Vector3<T> contactForce = m_tKContact1 * exp(m_tKContact2 * distance) * dir; 01601 m_prForce[i] -= contactForce; 01602 m_prForce[i+1] -= contactForce; 01603 m_prForce[j] += contactForce; 01604 m_prForce[j+1] += contactForce; 01605 //----------------------------------------- 01606 // Set WasIntersect Flags 01607 m_prWasIntersect[i] = m_prWasIntersect[j] = 1; 01608 } 01609 //*/ 01610 } 01611 } 01612 //---------------------------------------------------------------- 01613 return overlap; 01614 } 01615 //----------------------------------------------------------------------------- 01616 01617 #ifdef TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 01618 //============================================================================= 01619 // START: Check Crossings for Knot Recognition by Dowker Notation 01620 //----------------------------------------------------------------------------- 01624 #ifdef TAPs_ALLOW_SETTING_PROJECTION_DIRECTION 01625 template <typename T> 01626 std::string ModelStrand<T>::GetKnotName ( Vector3<T> const & projectionDirection ) 01627 { 01628 #ifdef TAPs_STRAND_DEBUG 01629 special_debug_knotID = 01630 #endif//TAPs_STRAND_DEBUG 01631 CheckSelfCrossingForKnotRecognition( projectionDirection ); 01632 std::string knotName; 01633 DowkerNotationToKnotName( &knotName ); 01634 return knotName; 01635 } 01636 #else //TAPs_ALLOW_SETTING_PROJECTION_DIRECTION 01637 template <typename T> 01638 std::string ModelStrand<T>::GetKnotName () 01639 { 01640 #ifdef TAPs_STRAND_DEBUG 01641 special_debug_knotID = 01642 #endif//TAPs_STRAND_DEBUG 01643 CheckSelfCrossingForKnotRecognition(); 01644 std::string knotName; 01645 DowkerNotationToKnotName( &knotName ); 01646 return knotName; 01647 } 01648 #endif//TAPs_ALLOW_SETTING_PROJECTION_DIRECTION 01649 //----------------------------------------------------------------------------- 01650 01651 //----------------------------------------------------------------------------- 01652 01653 //----------------------------------------------------------------------------- 01655 template <typename T> 01656 void ModelStrand<T>::CheckSelfCrossingForKnotRecognition_zDir () 01657 { 01658 //---------------------------------------------------------------- 01659 // Reset 01660 for ( int i = 0; i < m_iMaxCrossingPairs; ++i ) m_iaDowkerNotation[i] = 0; 01661 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01662 m_crossings.clear(); 01663 01664 #ifdef DEBUG_TOO 01665 m_crossings_TOO.clear(); 01666 #endif//DEBUG_TOO 01667 01668 for ( int i = 0; i < m_iNoLinks; ++i ) m_iaCrossingCounters[i] = 0; 01669 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01670 for ( int i = 0; i < m_iNoLinks; ++i ) m_iaCrossingsForDowkerNotation[i] = 0; 01671 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01672 01673 //---------------------------------------------------------------- 01674 // Find Crossings 01675 CheckSelfCrossingsNextLevelRecursively_zDir( m_BVHTree->Root()->Child(0), 01676 m_BVHTree->Root()->Child(1) ); 01677 //---------------------------------------------------------------- 01678 // For Recording Dowker Notation 01679 std::vector<int> crosses; 01680 std::vector<int> crossesID; 01681 01682 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01683 01684 #ifdef DEBUG_TOO 01685 01686 std::list< IC_Crossing >::const_iterator it = m_crossings_TOO.begin(); 01687 while ( it != m_crossings_TOO.end() ) { 01688 01689 #else //DEBUG_TOO 01690 01691 std::list< IC_Crossing >::const_iterator it = m_crossings.begin(); 01692 while ( it != m_crossings.end() ) { 01693 01694 #endif//DEBUG_TOO 01695 01696 // crossesID records the node number that has another node cross over or under it. 01697 // crosses records the crossing identified by node number that crosses over the crossesID. 01698 // E.g. real cross# 1 2 3 4 5 6 01699 // crossesID: 30 31 32 74 77 79 01700 // crosses: -74 77 -79 30 -31 32 01701 // results in the Dowker notation below 01702 // Dowker#: (1,4) (3,6) (5,2) 01703 // which is a half knot 01704 crossesID.push_back( (*it).id ); 01705 crosses.push_back( (*it).number ); 01706 ++it; 01707 } 01708 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01709 for ( int i = 0; i < m_iNoLinks; ++i ) { 01710 if ( m_iaCrossingsForDowkerNotation[i] != 0 ) { 01711 // crossesID records the node number that has another node cross over or under it. 01712 // crosses records the crossing identified by node number that crosses over the crossesID. 01713 // E.g. real cross# 1 2 3 4 5 6 01714 // crossesID: 30 31 32 74 77 79 01715 // crosses: -74 77 -79 30 -31 32 01716 // results in the Dowker notation below 01717 // Dowker#: (1,4) (3,6) (5,2) 01718 // which is a half knot 01719 crossesID.push_back( i+1 ); 01720 crosses.push_back( m_iaCrossingsForDowkerNotation[i] ); 01721 } 01722 } 01723 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01724 01725 int numOfCrosses = static_cast<int>( crosses.size() ); 01726 if ( numOfCrosses == 0 ) { // Not a knot 01727 m_iKnotStartPt = m_iKnotEndPt = -1; 01728 m_iNumCrossingPairs = 0; 01729 } 01730 else if ( numOfCrosses % 2 == 0 ) { 01731 m_iNumCrossingPairs = numOfCrosses/2; 01732 assert( m_iNumCrossingPairs <= m_iMaxCrossingPairs ); 01733 int idx = 0; 01734 for ( int i = 0, odd = 1; i < numOfCrosses; i+=2, odd+=2 ) { 01735 for ( int j = 1, even = 2; j < numOfCrosses; j+=2, even+=2 ) { 01736 if ( crossesID[i] == abs( crosses[j] ) ) { 01737 //iaDowkerNotation[ idx++ ] = odd; 01738 if ( crosses[j] < 0 ) 01739 m_iaDowkerNotation[ idx++ ] = -even; 01740 else 01741 m_iaDowkerNotation[ idx++ ] = even; 01742 } 01743 } 01744 } 01745 01747 m_iKnotStartPt = crossesID[0]; 01748 m_iKnotEndPt = crossesID[numOfCrosses-1]; 01749 } 01750 else { // Not a knot 01751 m_iKnotStartPt = m_iKnotEndPt = -1; 01752 m_iNumCrossingPairs = numOfCrosses/2; 01753 } 01754 01755 #ifdef TAPs_STRAND_DEBUG 01756 01757 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01758 /* 01759 { 01760 std::cout << "NUMs: "; 01761 std::list< IC_Crossing >::const_iterator it = m_crossings.begin(); 01762 while ( it != m_crossings.end() ) { 01763 std::cout << "\t" << (*it).number; 01764 ++it; 01765 } 01766 std::cout << "\n"; 01767 std::cout << "IDs: "; 01768 it = m_crossings.begin(); 01769 while ( it != m_crossings.end() ) { 01770 std::cout << "\t" << (*it).id; 01771 ++it; 01772 } 01773 std::cout << "\n"; 01774 } 01775 //*/ 01776 01777 /* 01778 #ifdef DEBUG_TOO 01779 { 01780 std::cout << "NUMs: "; 01781 std::list< IC_Crossing >::const_iterator it = m_crossings_TOO.begin(); 01782 while ( it != m_crossings_TOO.end() ) { 01783 std::cout << "\t" << (*it).number; 01784 ++it; 01785 } 01786 std::cout << "\n"; 01787 std::cout << "IDs: "; 01788 it = m_crossings_TOO.begin(); 01789 while ( it != m_crossings_TOO.end() ) { 01790 std::cout << "\t" << (*it).id; 01791 ++it; 01792 } 01793 std::cout << "\n"; 01794 } 01795 #endif//DEBUG_TOO 01796 //*/ 01797 01798 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01799 01800 debug_crosses = crosses; 01801 debug_crossesID = crossesID; 01802 std::cout << Debug_Str(); 01803 #endif//TAPs_STRAND_DEBUG 01804 01805 } 01806 //----------------------------------------------------------------------------- 01807 01808 //----------------------------------------------------------------------------- 01810 template <typename T> 01811 void ModelStrand<T>::CheckSelfCrossingForKnotRecognition_yDir () 01812 { 01813 //---------------------------------------------------------------- 01814 // Reset 01815 for ( int i = 0; i < m_iMaxCrossingPairs; ++i ) m_iaDowkerNotation[i] = 0; 01816 for ( int i = 0; i < m_iNoLinks; ++i ) m_iaCrossingsForDowkerNotation[i] = 0; 01817 //---------------------------------------------------------------- 01818 // Find Crossings 01819 CheckSelfCrossingsNextLevelRecursively_yDir( m_BVHTree->Root()->Child(0), 01820 m_BVHTree->Root()->Child(1), 01821 m_iaCrossingsForDowkerNotation ); 01822 //---------------------------------------------------------------- 01823 // Record Dowker Notation 01824 std::vector<int> crosses; 01825 std::vector<int> crossesID; 01826 for ( int i = 0; i < m_iNoLinks; ++i ) { 01827 if ( m_iaCrossingsForDowkerNotation[i] != 0 ) { 01828 crosses.push_back( m_iaCrossingsForDowkerNotation[i] ); 01829 crossesID.push_back( i+1 ); 01830 } 01831 } 01832 int numOfCrosses = static_cast<int>( crosses.size() ); 01833 if ( numOfCrosses == 0 ) { // Not a knot 01834 m_iKnotStartPt = m_iKnotEndPt = -1; 01835 m_iNumCrossingPairs = 0; 01836 } 01837 else if ( numOfCrosses % 2 == 0 ) { 01838 m_iNumCrossingPairs = numOfCrosses/2; 01839 assert( m_iNumCrossingPairs <= m_iMaxCrossingPairs ); 01840 int idx = 0; 01841 for ( int i = 0, odd = 1; i < numOfCrosses; i+=2, odd+=2 ) { 01842 for ( int j = 1, even = 2; j < numOfCrosses; j+=2, even+=2 ) { 01843 if ( crossesID[i] == abs( crosses[j] ) ) { 01844 //iaDowkerNotation[ idx++ ] = odd; 01845 if ( crosses[j] < 0 ) 01846 m_iaDowkerNotation[ idx++ ] = -even; 01847 else 01848 m_iaDowkerNotation[ idx++ ] = even; 01849 } 01850 } 01851 } 01852 01854 m_iKnotStartPt = crossesID[0]; 01855 m_iKnotEndPt = crossesID[numOfCrosses-1]; 01856 } 01857 else { // Not a knot 01858 m_iKnotStartPt = m_iKnotEndPt = -1; 01859 m_iNumCrossingPairs = numOfCrosses/2; 01860 } 01861 } 01862 //----------------------------------------------------------------------------- 01863 01864 //----------------------------------------------------------------------------- 01866 template <typename T> 01867 void ModelStrand<T>::CheckSelfCrossingForKnotRecognition_xDir () 01868 { 01869 //---------------------------------------------------------------- 01870 // Reset 01871 for ( int i = 0; i < m_iMaxCrossingPairs; ++i ) m_iaDowkerNotation[i] = 0; 01872 for ( int i = 0; i < m_iNoLinks; ++i ) m_iaCrossingsForDowkerNotation[i] = 0; 01873 //---------------------------------------------------------------- 01874 // Find Crossings 01875 CheckSelfCrossingsNextLevelRecursively_xDir( m_BVHTree->Root()->Child(0), 01876 m_BVHTree->Root()->Child(1), 01877 m_iaCrossingsForDowkerNotation ); 01878 //---------------------------------------------------------------- 01879 // Record Dowker Notation 01880 std::vector<int> crosses; 01881 std::vector<int> crossesID; 01882 for ( int i = 0; i < m_iNoLinks; ++i ) { 01883 if ( m_iaCrossingsForDowkerNotation[i] != 0 ) { 01884 crosses.push_back( m_iaCrossingsForDowkerNotation[i] ); 01885 crossesID.push_back( i+1 ); 01886 } 01887 } 01888 int numOfCrosses = static_cast<int>( crosses.size() ); 01889 if ( numOfCrosses == 0 ) { // Not a knot 01890 m_iKnotStartPt = m_iKnotEndPt = -1; 01891 m_iNumCrossingPairs = 0; 01892 } 01893 else if ( numOfCrosses % 2 == 0 ) { 01894 m_iNumCrossingPairs = numOfCrosses/2; 01895 assert( m_iNumCrossingPairs <= m_iMaxCrossingPairs ); 01896 int idx = 0; 01897 for ( int i = 0, odd = 1; i < numOfCrosses; i+=2, odd+=2 ) { 01898 for ( int j = 1, even = 2; j < numOfCrosses; j+=2, even+=2 ) { 01899 if ( crossesID[i] == abs( crosses[j] ) ) { 01900 //iaDowkerNotation[ idx++ ] = odd; 01901 if ( crosses[j] < 0 ) 01902 m_iaDowkerNotation[ idx++ ] = -even; 01903 else 01904 m_iaDowkerNotation[ idx++ ] = even; 01905 } 01906 } 01907 } 01908 01910 m_iKnotStartPt = crossesID[0]; 01911 m_iKnotEndPt = crossesID[numOfCrosses-1]; 01912 } 01913 else { // Not a knot 01914 m_iKnotStartPt = m_iKnotEndPt = -1; 01915 m_iNumCrossingPairs = numOfCrosses/2; 01916 } 01917 } 01918 //----------------------------------------------------------------------------- 01919 01920 #ifdef TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 01921 //----------------------------------------------------------------------------- 01922 // CreateDataSpaceForKnotRecognitionByDowkerNotation 01923 template <typename T> 01924 bool ModelStrand<T>::CheckClumped ( int startPt, int endPt, T sphereRadius ) 01925 { 01926 // Return false if the clump is less then 6 links. 01927 if ( endPt - startPt < 5 ) { 01928 m_bIsClumped = false; 01929 return m_bIsClumped; 01930 } 01931 01932 #ifdef TAPs_STRAND_LOCK_SEGMENTS 01933 // Determine the sphere center 01934 Vector3<T> avgPt( 0,0,0 ); 01935 int count = 0; 01936 for ( int i = startPt; i <= endPt; ++i ) { 01937 if ( !m_prIsLocked[i] ) { // skip locked point 01938 avgPt += m_prVertex[i]; 01939 ++count; 01940 } 01941 } 01942 avgPt /= count; 01943 01944 // Return true if all points are within the sphere, else return false. 01945 for ( int i = startPt; i <= endPt; ++i ) { 01946 if ( !m_prIsLocked[i] ) { // skip locked point 01947 T length = ( avgPt - m_prVertex[i] ).Length(); 01948 if ( length > sphereRadius ) { 01949 m_bIsClumped = false; 01950 return m_bIsClumped; 01951 } 01952 } 01953 } 01954 #else //TAPs_STRAND_LOCK_SEGMENTS 01955 // Determine the sphere center 01956 Vector3<T> avgPt( m_prVertex[startPt] ); 01957 for ( int i = startPt+1; i <= endPt; ++i ) { 01958 avgPt += m_prVertex[i]; 01959 } 01960 avgPt /= ( endPt - startPt + 1 ); 01961 01962 // Return true if all points are within the sphere, else return false. 01963 for ( int i = startPt; i <= endPt; ++i ) { 01964 T length = ( avgPt - m_prVertex[i] ).Length(); 01965 if ( length > sphereRadius ) { 01966 m_bIsClumped = false; 01967 return m_bIsClumped; 01968 } 01969 } 01970 #endif//TAPs_STRAND_LOCK_SEGMENTS 01971 01972 m_bIsClumped = true; 01973 return m_bIsClumped; 01974 } 01975 #endif//TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 01976 //----------------------------------------------------------------------------- 01977 01978 //----------------------------------------------------------------------------- 01979 // CreateDataSpaceForKnotRecognitionByDowkerNotation 01980 template <typename T> 01981 void ModelStrand<T>::CreateDataSpaceForKnotRecognitionByDowkerNotation () 01982 { 01983 //m_ucNumKnots = 0; //! Init number of knots on the strand to zero 01984 m_iKnotStartPt = m_iKnotEndPt = -1; 01985 //---------------------------------------------------------------- 01986 m_iMaxCrossingPairs = m_iNoLinks/2; 01987 m_iNumCrossingPairs = 0; 01988 //---------------------------------------------------------------- 01989 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 01990 if ( ( m_iaCrossingCounters = new int[ m_iNoLinks ] ) == NULL ) { 01991 #ifdef TAPs_ENABLE_DEBUG 01992 std::cerr << "ERROR => ModelStrand Constructor:" 01993 << " Cannot allocate memory for knot recognition by Dowker notation!" 01994 << std::endl; 01995 #endif 01996 delete this; 01997 return; 01998 } 01999 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02000 if ( ( m_iaCrossingsForDowkerNotation = new int[ m_iNoLinks ] ) == NULL ) { 02001 #ifdef TAPs_ENABLE_DEBUG 02002 std::cerr << "ERROR => ModelStrand Constructor:" 02003 << " Cannot allocate memory for knot recognition by Dowker notation!" 02004 << std::endl; 02005 #endif 02006 delete this; 02007 return; 02008 } 02009 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02010 //---------------------------------------------------------------- 02011 if ( ( m_iaDowkerNotation = new int[ m_iMaxCrossingPairs ] ) == NULL ) { 02012 #ifdef TAPs_ENABLE_DEBUG 02013 std::cerr << "ERROR => ModelStrand Constructor:" 02014 << " Cannot allocate memory for knot recognition by Dowker notation!" 02015 << std::endl; 02016 #endif 02017 delete this; 02018 return; 02019 } 02020 } 02021 //----------------------------------------------------------------------------- 02022 // DeleteDataSpaceForKnotRecognitionByDowkerNotation 02023 template <typename T> 02024 void ModelStrand<T>::DeleteDataSpaceForKnotRecognitionByDowkerNotation () 02025 { 02026 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02027 if ( m_iaCrossingCounters ) { 02028 delete [] m_iaCrossingCounters; 02029 m_iaCrossingCounters = NULL; 02030 } 02031 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02032 if ( m_iaCrossingsForDowkerNotation ) { 02033 delete [] m_iaCrossingsForDowkerNotation; 02034 m_iaCrossingsForDowkerNotation = NULL; 02035 } 02036 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02037 02038 if ( m_iaDowkerNotation ) { 02039 delete [] m_iaDowkerNotation; 02040 m_iaDowkerNotation = NULL; 02041 } 02042 } 02043 //----------------------------------------------------------------------------- 02044 02045 //============================================================================= 02046 // START: Fns for CheckSelfCrossingForKnotRecognition_zDir function 02047 //----------------------------------------------------------------------------- 02048 // CheckSelfCrossingsNextLevelRecursively 02049 // Assume no intersections between immediate adjacent links 02050 // Assume a BVH node either has both children or none, therefore 02051 // checking either it has a left or right child is enough to determine 02052 // whether it is a leaf node or not. 02053 template <typename T> 02054 void ModelStrand<T>::CheckSelfCrossingsNextLevelRecursively_zDir ( 02055 BVHNode<T> * node1, 02056 BVHNode<T> * node2 02057 ) 02058 { 02059 CheckSelfCrossingsRecursively_zDir( node1, node2 ); 02060 if ( node1->Child(0) ) 02061 CheckSelfCrossingsNextLevelRecursively_zDir( node1->Child(0), node1->Child(1) ); 02062 if ( node2->Child(0) ) 02063 CheckSelfCrossingsNextLevelRecursively_zDir( node2->Child(0), node2->Child(1) ); 02064 } 02065 //----------------------------------------------------------------------------- 02066 // CheckSelfCrossingsRecursively 02067 // Assume no intersections between immediate adjacent links 02068 // Use the fact that the ids of sphere bvh leaf nodes are the same as 02069 // the link id and 02070 template <typename T> 02071 T ModelStrand<T>::CheckSelfCrossingsRecursively_zDir ( 02072 BVHNode<T> * node1, 02073 BVHNode<T> * node2 02074 ) 02075 { 02076 static T LinkDiameter = static_cast<T>(2.0) * m_tRadius; 02077 static T epsilon = m_tRadius * static_cast<T>(0.05); // for collision response 02078 //--------------------------------------------------------------- 02079 if ( !node1->Child(0) && !node2->Child(0) && 02080 ( abs(node1->GetID() - node2->GetID()) == 1 ) ) 02081 { 02082 return DBL_MAX; 02083 } 02084 //--------------------------------------------------------------- 02085 T overlap = TestOverlapCircle_zDir( node1, node2 ); 02086 //--------------------------------------------------------------- 02087 if ( overlap < DBL_MIN ) { 02088 if ( node2->Child(0) == NULL ) { //&& node2->RightChild() == NULL ) { 02089 if ( node1->Child(0) ) { 02090 overlap = CheckSelfCrossingsRecursively_zDir( node1->Child(0), node2 ); 02091 } 02092 if ( node1->Child(1) ) { 02093 overlap = CheckSelfCrossingsRecursively_zDir( node1->Child(1), node2 ); 02094 } 02095 } 02096 else { 02097 if ( node2->Child(0) ) { 02098 overlap = TestOverlapCircle_zDir( node1, node2->Child(0) ); 02099 if ( overlap < DBL_MIN ) { 02100 if ( node1->Child(0) ) { 02101 overlap = CheckSelfCrossingsRecursively_zDir( node1->Child(0), node2->Child(0) ); 02102 } 02103 if ( node1->Child(1) ) { 02104 overlap = CheckSelfCrossingsRecursively_zDir( node1->Child(1), node2->Child(0) ); 02105 } 02106 } 02107 } 02108 if ( node2->Child(1) ) { 02109 overlap = TestOverlapCircle_zDir( node1, node2->Child(1) ); 02110 if ( overlap < DBL_MIN ) { 02111 if ( node1->Child(0) ) { 02112 overlap = CheckSelfCrossingsRecursively_zDir( node1->Child(0), node2->Child(1) ); 02113 } 02114 if ( node1->Child(1) ) { 02115 overlap = CheckSelfCrossingsRecursively_zDir( node1->Child(1), node2->Child(1) ); 02116 } 02117 } 02118 } 02119 } 02120 //----------------------------------------------------------- 02121 if ( node1->Child(0) == NULL && //node1->Child(1) == NULL && 02122 node2->Child(0) == NULL ) //&& node2->Child(1) == NULL ) 02123 { 02124 CheckCrossings_zDir( node1, node2 ); 02125 } 02126 } 02127 //--------------------------------------------------------------- 02128 return overlap; 02129 } 02130 //----------------------------------------------------------------------------- 02131 template <typename T> 02132 T ModelStrand<T>::TestOverlapCircle_zDir ( 02133 BVHNode<T> * node1, BVHNode<T> * node2 ) 02134 { 02135 Vector3<T> C1 = node1->GetCenter(); 02136 Vector3<T> C2 = node2->GetCenter(); 02137 C1.SetZ( 0 ); 02138 C2.SetZ( 0 ); 02139 return (C1 - C2).Length() - ( node1->GetRadius() + node2->GetRadius() ); 02140 } 02141 //----------------------------------------------------------------------------- 02142 template <typename T> 02143 void ModelStrand<T>::CheckCrossings_zDir ( 02144 BVHNode<T> * node1, BVHNode<T> * node2 ) 02145 { 02146 int i = node1->GetID(); 02147 int j = node2->GetID(); 02148 02149 // Skip if a node is locked 02150 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02151 if ( m_prIsLocked[i] ) return; 02152 if ( m_prIsLocked[j] ) return; 02153 #endif//TAPs_STRAND_LOCK_SEGMENTS 02154 02155 // Link# i 02156 Vector3<T> p1 = m_prVertex[i]; 02157 Vector3<T> p2 = m_prVertex[i+1]; 02158 // Link# j 02159 Vector3<T> q1 = m_prVertex[j]; 02160 Vector3<T> q2 = m_prVertex[j+1]; 02161 02162 // Project in z-direction 02163 p1.SetZ( 0 ); 02164 p2.SetZ( 0 ); 02165 q1.SetZ( 0 ); 02166 q2.SetZ( 0 ); 02167 02168 bool isTouched; 02169 T pt, qt; // parameter for the 1st & 2nd line 02170 // where intersection point = p1 + pt(p2-p1) 02171 // where intersection point = q1 + qt(q2-q1) 02172 CGMath<T>::FindIntersectionLineSegmentLineSegment( 02173 p1, p2, q1, q2, // I/P 02174 pt, qt, // O/P 02175 isTouched // O/P 02176 ); 02177 if ( 0 <= pt && pt <= 1 && 0 <= qt && qt <= 1 ) { 02178 Vector3<T> P = m_prVertex[i] + pt*( m_prVertex[i+1] - m_prVertex[i] ); 02179 Vector3<T> Q = m_prVertex[j] + qt*( m_prVertex[j+1] - m_prVertex[j] ); 02180 02181 // i+1 and j+1 are for making the link# to be one-based instead of zero-based 02182 if ( P[2] > Q[2] ) { 02183 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02184 02185 #ifdef DEBUG_TOO 02186 int a = i*10 + m_iaCrossingCounters[i]++; 02187 int b = j*10 + m_iaCrossingCounters[j]++; 02188 InsertToList_TOO( -b, a ); 02189 InsertToList_TOO( a, b ); 02190 #endif//DEBUG_TOO 02191 02192 InsertToList( -(j), i ); 02193 InsertToList( +(i), j ); 02194 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02195 m_iaCrossingsForDowkerNotation[i] = -(j+1); // the link i is above the link j 02196 m_iaCrossingsForDowkerNotation[j] = +(i+1); // the link j is below the link i 02197 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02198 } 02199 else { 02200 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02201 02202 #ifdef DEBUG_TOO 02203 int a = i*10 + m_iaCrossingCounters[i]++; 02204 int b = j*10 + m_iaCrossingCounters[j]++; 02205 InsertToList_TOO( b, a ); 02206 InsertToList_TOO( -a, b ); 02207 #endif//DEBUG_TOO 02208 02209 InsertToList( +(j), i ); 02210 InsertToList( -(i), j ); 02211 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02212 m_iaCrossingsForDowkerNotation[i] = +(j+1); // the link j is above the link i 02213 m_iaCrossingsForDowkerNotation[j] = -(i+1); // the link i is below the link j 02214 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 02215 } 02216 } 02217 } 02218 //----------------------------------------------------------------------------- 02219 // END: Fns for CheckSelfCrossingForKnotRecognition_zDir function 02220 //============================================================================= 02221 02222 02223 02224 //============================================================================= 02225 // START: Fns for CheckSelfCrossingForKnotRecognition_yDir function 02226 //----------------------------------------------------------------------------- 02227 // CheckSelfCrossingsNextLevelRecursively 02228 // Assume no intersections between immediate adjacent links 02229 // Assume a BVH node either has both children or none, therefore 02230 // checking either it has a left or right child is enough to determine 02231 // whether it is a leaf node or not. 02232 template <typename T> 02233 void ModelStrand<T>::CheckSelfCrossingsNextLevelRecursively_yDir ( 02234 BVHNode<T> * node1, 02235 BVHNode<T> * node2 02236 ) 02237 { 02238 CheckSelfCrossingsRecursively_yDir( node1, node2 ); 02239 if ( node1->Child(0) ) 02240 CheckSelfCrossingsNextLevelRecursively_yDir( node1->Child(0), node1->Child(1) ); 02241 if ( node2->Child(0) ) 02242 CheckSelfCrossingsNextLevelRecursively_yDir( node2->Child(0), node2->Child(1) ); 02243 } 02244 //----------------------------------------------------------------------------- 02245 // CheckSelfCrossingsRecursively 02246 // Assume no intersections between immediate adjacent links 02247 // Use the fact that the ids of sphere bvh leaf nodes are the same as 02248 // the link id and 02249 template <typename T> 02250 T ModelStrand<T>::CheckSelfCrossingsRecursively_yDir ( 02251 BVHNode<T> * node1, 02252 BVHNode<T> * node2 02253 ) 02254 { 02255 static T LinkDiameter = static_cast<T>(2.0) * m_tRadius; 02256 static T epsilon = m_tRadius * static_cast<T>(0.05); // for collision response 02257 //--------------------------------------------------------------- 02258 if ( !node1->Child(0) && !node2->Child(0) && 02259 ( abs(node1->GetID() - node2->GetID()) == 1 ) ) 02260 { 02261 return DBL_MAX; 02262 } 02263 //--------------------------------------------------------------- 02264 T overlap = TestOverlapCircle_yDir( node1, node2 ); 02265 //--------------------------------------------------------------- 02266 if ( overlap < DBL_MIN ) { 02267 if ( node2->Child(0) == NULL ) { //&& node2->RightChild() == NULL ) { 02268 if ( node1->Child(0) ) { 02269 overlap = CheckSelfCrossingsRecursively_yDir( node1->Child(0), node2 ); 02270 } 02271 if ( node1->Child(1) ) { 02272 overlap = CheckSelfCrossingsRecursively_yDir( node1->Child(1), node2 ); 02273 } 02274 } 02275 else { 02276 if ( node2->Child(0) ) { 02277 overlap = TestOverlapCircle_yDir( node1, node2->Child(0) ); 02278 if ( overlap < DBL_MIN ) { 02279 if ( node1->Child(0) ) { 02280 overlap = CheckSelfCrossingsRecursively_yDir( node1->Child(0), node2->Child(0) ); 02281 } 02282 if ( node1->Child(1) ) { 02283 overlap = CheckSelfCrossingsRecursively_yDir( node1->Child(1), node2->Child(0) ); 02284 } 02285 } 02286 } 02287 if ( node2->Child(1) ) { 02288 overlap = TestOverlapCircle_yDir( node1, node2->Child(1) ); 02289 if ( overlap < DBL_MIN ) { 02290 if ( node1->Child(0) ) { 02291 overlap = CheckSelfCrossingsRecursively_yDir( node1->Child(0), node2->Child(1) ); 02292 } 02293 if ( node1->Child(1) ) { 02294 overlap = CheckSelfCrossingsRecursively_yDir( node1->Child(1), node2->Child(1) ); 02295 } 02296 } 02297 } 02298 } 02299 //----------------------------------------------------------- 02300 if ( node1->Child(0) == NULL && //node1->Child(1) == NULL && 02301 node2->Child(0) == NULL ) //&& node2->Child(1) == NULL ) 02302 { 02303 CheckCrossings_yDir( node1, node2 ); 02304 } 02305 } 02306 //--------------------------------------------------------------- 02307 return overlap; 02308 } 02309 //----------------------------------------------------------------------------- 02310 template <typename T> 02311 T ModelStrand<T>::TestOverlapCircle_yDir ( 02312 BVHNode<T> * node1, BVHNode<T> * node2 ) 02313 { 02314 Vector3<T> C1 = node1->GetCenter(); 02315 Vector3<T> C2 = node2->GetCenter(); 02316 C1.SetY( 0 ); 02317 C2.SetY( 0 ); 02318 return (C1 - C2).Length() - ( node1->GetRadius() + node2->GetRadius() ); 02319 } 02320 //----------------------------------------------------------------------------- 02321 template <typename T> 02322 void ModelStrand<T>::CheckCrossings_yDir ( 02323 BVHNode<T> * node1, BVHNode<T> * node2 ) 02324 { 02325 int i = node1->GetID(); 02326 int j = node2->GetID(); 02327 02328 // Skip if a node is locked 02329 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02330 if ( m_prIsLocked[i] ) return; 02331 if ( m_prIsLocked[j] ) return; 02332 #endif//TAPs_STRAND_LOCK_SEGMENTS 02333 02334 // Link# i 02335 Vector3<T> p1 = m_prVertex[i]; 02336 Vector3<T> p2 = m_prVertex[i+1]; 02337 // Link# j 02338 Vector3<T> q1 = m_prVertex[j]; 02339 Vector3<T> q2 = m_prVertex[j+1]; 02340 02341 // Project in z-direction 02342 p1.SetY( 0 ); 02343 p2.SetY( 0 ); 02344 q1.SetY( 0 ); 02345 q2.SetY( 0 ); 02346 02347 bool isTouched; 02348 T pt, qt; // parameter for the 1st & 2nd line 02349 // where intersection point = p1 + pt(p2-p1) 02350 // where intersection point = q1 + qt(q2-q1) 02351 CGMath<T>::FindIntersectionLineSegmentLineSegment( 02352 p1, p2, q1, q2, // I/P 02353 pt, qt, // O/P 02354 isTouched // O/P 02355 ); 02356 if ( 0 <= pt && pt <= 1 && 0 <= qt && qt <= 1 ) { 02357 Vector3<T> P = m_prVertex[i] + pt*( m_prVertex[i+1] - m_prVertex[i] ); 02358 Vector3<T> Q = m_prVertex[j] + qt*( m_prVertex[j+1] - m_prVertex[j] ); 02359 02360 // i+1 and j+1 are for making the link# to be one-based instead of zero-based 02361 if ( P[2] > Q[2] ) { 02362 m_iaCrossingsForDowkerNotation[i] = -(j+1); // the link i is above the link j 02363 m_iaCrossingsForDowkerNotation[j] = +(i+1); // the link j is below the link i 02364 } 02365 else { 02366 m_iaCrossingsForDowkerNotation[i] = +(j+1); // the link j is above the link i 02367 m_iaCrossingsForDowkerNotation[j] = -(i+1); // the link i is below the link j 02368 } 02369 } 02370 } 02371 //----------------------------------------------------------------------------- 02372 // END: Fns for CheckSelfCrossingForKnotRecognition_yDir function 02373 //============================================================================= 02374 02375 02376 02377 //============================================================================= 02378 // START: Fns for CheckSelfCrossingForKnotRecognition_xDir function 02379 //----------------------------------------------------------------------------- 02380 // CheckSelfCrossingsNextLevelRecursively 02381 // Assume no intersections between immediate adjacent links 02382 // Assume a BVH node either has both children or none, therefore 02383 // checking either it has a left or right child is enough to determine 02384 // whether it is a leaf node or not. 02385 template <typename T> 02386 void ModelStrand<T>::CheckSelfCrossingsNextLevelRecursively_xDir ( 02387 BVHNode<T> * node1, 02388 BVHNode<T> * node2 02389 ) 02390 { 02391 CheckSelfCrossingsRecursively_xDir( node1, node2 ); 02392 if ( node1->Child(0) ) 02393 CheckSelfCrossingsNextLevelRecursively_xDir( node1->Child(0), node1->Child(1) ); 02394 if ( node2->Child(0) ) 02395 CheckSelfCrossingsNextLevelRecursively_xDir( node2->Child(0), node2->Child(1) ); 02396 } 02397 //----------------------------------------------------------------------------- 02398 // CheckSelfCrossingsRecursively 02399 // Assume no intersections between immediate adjacent links 02400 // Use the fact that the ids of sphere bvh leaf nodes are the same as 02401 // the link id and 02402 template <typename T> 02403 T ModelStrand<T>::CheckSelfCrossingsRecursively_xDir ( 02404 BVHNode<T> * node1, 02405 BVHNode<T> * node2 02406 ) 02407 { 02408 static T LinkDiameter = static_cast<T>(2.0) * m_tRadius; 02409 static T epsilon = m_tRadius * static_cast<T>(0.05); // for collision response 02410 //--------------------------------------------------------------- 02411 if ( !node1->Child(0) && !node2->Child(0) && 02412 ( abs(node1->GetID() - node2->GetID()) == 1 ) ) 02413 { 02414 return DBL_MAX; 02415 } 02416 //--------------------------------------------------------------- 02417 T overlap = TestOverlapCircle_xDir( node1, node2 ); 02418 //--------------------------------------------------------------- 02419 if ( overlap < DBL_MIN ) { 02420 if ( node2->Child(0) == NULL ) { //&& node2->RightChild() == NULL ) { 02421 if ( node1->Child(0) ) { 02422 overlap = CheckSelfCrossingsRecursively_xDir( node1->Child(0), node2 ); 02423 } 02424 if ( node1->Child(1) ) { 02425 overlap = CheckSelfCrossingsRecursively_xDir( node1->Child(1), node2 ); 02426 } 02427 } 02428 else { 02429 if ( node2->Child(0) ) { 02430 overlap = TestOverlapCircle_xDir( node1, node2->Child(0) ); 02431 if ( overlap < DBL_MIN ) { 02432 if ( node1->Child(0) ) { 02433 overlap = CheckSelfCrossingsRecursively_xDir( node1->Child(0), node2->Child(0) ); 02434 } 02435 if ( node1->Child(1) ) { 02436 overlap = CheckSelfCrossingsRecursively_xDir( node1->Child(1), node2->Child(0) ); 02437 } 02438 } 02439 } 02440 if ( node2->Child(1) ) { 02441 overlap = TestOverlapCircle_xDir( node1, node2->Child(1) ); 02442 if ( overlap < DBL_MIN ) { 02443 if ( node1->Child(0) ) { 02444 overlap = CheckSelfCrossingsRecursively_xDir( node1->Child(0), node2->Child(1) ); 02445 } 02446 if ( node1->Child(1) ) { 02447 overlap = CheckSelfCrossingsRecursively_xDir( node1->Child(1), node2->Child(1) ); 02448 } 02449 } 02450 } 02451 } 02452 //----------------------------------------------------------- 02453 if ( node1->Child(0) == NULL && //node1->Child(1) == NULL && 02454 node2->Child(0) == NULL ) //&& node2->Child(1) == NULL ) 02455 { 02456 CheckCrossings_xDir( node1, node2 ); 02457 } 02458 } 02459 //--------------------------------------------------------------- 02460 return overlap; 02461 } 02462 //----------------------------------------------------------------------------- 02463 template <typename T> 02464 T ModelStrand<T>::TestOverlapCircle_xDir ( 02465 BVHNode<T> * node1, BVHNode<T> * node2 ) 02466 { 02467 Vector3<T> C1 = node1->GetCenter(); 02468 Vector3<T> C2 = node2->GetCenter(); 02469 C1.SetX( 0 ); 02470 C2.SetX( 0 ); 02471 return (C1 - C2).Length() - ( node1->GetRadius() + node2->GetRadius() ); 02472 } 02473 //----------------------------------------------------------------------------- 02474 template <typename T> 02475 void ModelStrand<T>::CheckCrossings_xDir ( 02476 BVHNode<T> * node1, BVHNode<T> * node2 ) 02477 { 02478 int i = node1->GetID(); 02479 int j = node2->GetID(); 02480 02481 // Skip if a node is locked 02482 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02483 if ( m_prIsLocked[i] ) return; 02484 if ( m_prIsLocked[j] ) return; 02485 #endif//TAPs_STRAND_LOCK_SEGMENTS 02486 02487 // Link# i 02488 Vector3<T> p1 = m_prVertex[i]; 02489 Vector3<T> p2 = m_prVertex[i+1]; 02490 // Link# j 02491 Vector3<T> q1 = m_prVertex[j]; 02492 Vector3<T> q2 = m_prVertex[j+1]; 02493 02494 // Project in z-direction 02495 p1.SetX( 0 ); 02496 p2.SetX( 0 ); 02497 q1.SetX( 0 ); 02498 q2.SetX( 0 ); 02499 02500 bool isTouched; 02501 T pt, qt; // parameter for the 1st & 2nd line 02502 // where intersection point = p1 + pt(p2-p1) 02503 // where intersection point = q1 + qt(q2-q1) 02504 CGMath<T>::FindIntersectionLineSegmentLineSegment( 02505 p1, p2, q1, q2, // I/P 02506 pt, qt, // O/P 02507 isTouched // O/P 02508 ); 02509 if ( 0 <= pt && pt <= 1 && 0 <= qt && qt <= 1 ) { 02510 Vector3<T> P = m_prVertex[i] + pt*( m_prVertex[i+1] - m_prVertex[i] ); 02511 Vector3<T> Q = m_prVertex[j] + qt*( m_prVertex[j+1] - m_prVertex[j] ); 02512 02513 // i+1 and j+1 are for making the link# to be one-based instead of zero-based 02514 if ( P[2] > Q[2] ) { 02515 m_iaCrossingsForDowkerNotation[i] = -(j+1); // the link i is above the link j 02516 m_iaCrossingsForDowkerNotation[j] = +(i+1); // the link j is below the link i 02517 } 02518 else { 02519 m_iaCrossingsForDowkerNotation[i] = +(j+1); // the link j is above the link i 02520 m_iaCrossingsForDowkerNotation[j] = -(i+1); // the link i is below the link j 02521 } 02522 } 02523 } 02524 //----------------------------------------------------------------------------- 02525 // END: Fns for CheckSelfCrossingForKnotRecognition_xDir function 02526 //============================================================================= 02527 02528 02529 02530 //----------------------------------------------------------------------------- 02531 template <typename T> 02532 int ModelStrand<T>::DowkerNotationToKnotName ( std::string * knotName ) 02533 { 02534 // Find the knot that match with the Dowker notation 02535 int knotID = Data::DatabaseDowkerNotation::FindKnotOfDowkerNotation( 02536 m_iNumCrossingPairs, m_iaDowkerNotation, knotName ); 02537 02538 // Include the Dowker notation in the return string 02539 if ( knotName ) { 02540 if ( knotID >= 0 ) { 02541 *knotName = Data::DatabaseDowkerNotation::GetStrEntry( knotID ); 02542 } 02543 else { 02544 *knotName = "n/a"; 02545 } 02546 } 02547 02548 return knotID; // return -1 if not found, otherwise the knotID 02549 } 02550 //----------------------------------------------------------------------------- 02551 // End Check Crossings for Knot Recognition by Dowker Notation 02552 //============================================================================= 02553 #endif//TAPs_STRAND_KNOT_RECOGNITION_BY_DOWKER_NOTATION 02554 02555 //----------------------------------------------------------------------------- 02556 // Reset 02557 template <typename T> 02558 void ModelStrand<T>::Reset () 02559 { 02560 for ( int i = 0; i < m_iNoLinks; ++i ) { 02561 //--------------------------------------- 02562 m_prIsFixed[i] = false; 02563 //--------------------------------------- 02564 //m_prWasIntersect[i] = 0; 02565 //--------------------------------------- 02566 m_prVelocity[i].SetXYZ( 0, 0, 0 ); 02567 m_prForce[i].SetXYZ( 0, 0, 0 ); 02568 //--------------------------------------- 02569 #ifdef TAPs_ADVANCED_SIMULATION 02570 m_SimFlags[i].ClearAllFlags(); 02571 #endif//TAPs_ADVANCED_SIMULATION 02572 } 02573 // m_vListOfFixedPoints.clear(); // Unused (i.e. doesn't exist!) 02574 m_setOfFixedPts.clear(); 02575 02576 #ifdef TAPs_USE_EXTERNAL_ODE_SOLVER 02577 for ( unsigned int i = 0; i < x0.size(); ++i ) { 02578 x0[i] = 0; 02579 xEnd[i] = 0; 02580 } 02581 #endif//TAPs_USE_EXTERNAL_ODE_SOLVER 02582 02583 //=============================================================== 02584 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02585 //--------------------------------------------------------------- 02586 ClearSegments(); 02587 //--------------------------------------------------------------- 02588 #endif//TAPs_STRAND_LOCK_SEGMENTS 02589 //=============================================================== 02590 02591 //=============================================================== 02592 #ifdef TAPs_STRAND_KNOT_CONTROL 02593 //--------------------------------------------------------------- 02594 ResetKnotControl(); 02595 //--------------------------------------------------------------- 02596 #endif//TAPs_STRAND_KNOT_CONTROL 02597 //=============================================================== 02598 02599 } // EOF: Reset () 02600 //----------------------------------------------------------------------------- 02601 // InitCurrentShape 02602 template <typename T> 02603 void ModelStrand<T>::InitCurrentShape () 02604 { 02605 static int m_iLoopSize = 16; 02606 } 02607 //----------------------------------------------------------------------------- 02608 // SetupShape Fn 02609 template <typename T> 02610 void ModelStrand<T>::SetupShape ( Vector3<T> startPosition ) 02611 { 02612 //for ( int i = 0; i <= m_iNoLinks; ++i ) { 02613 // startPosition[0] += m_tLinkLength; 02614 // m_prVertex[i] = startPosition; 02615 //} 02616 02617 /* 02618 // Case m_iNoLinks == 6 02619 m_prVertex[0] = Vector3<T>( 0,0,0 ); 02620 m_prVertex[1] = Vector3<T>( 1,0,0 ); 02621 m_prVertex[2] = Vector3<T>( 1,1,0 ); 02622 m_prVertex[3] = Vector3<T>( 1,1,1 ); 02623 m_prVertex[4] = Vector3<T>( 2,1,1 ); 02624 m_prVertex[5] = Vector3<T>( 2,0,1 ); 02625 m_prVertex[6] = Vector3<T>( 3,0,1 ); 02626 //*/ 02627 02628 //* 02629 //------------------------------------------- 02630 int i = 0; 02631 T posNev = 1; 02632 T countStep = 0.01; 02633 T count = countStep; 02634 T deg = 10; T step = -17; 02635 //T deg = 120; T step = -50; 02636 m_prVertex[i++] = startPosition; 02637 while ( i <= m_iNoLinks ) { 02638 Vector3<T> v( m_tLinkLength, 0, 0 ); 02639 T rad = deg * 3.14156 / 180.0; 02640 T c = cos(rad); 02641 T s = sin(rad); 02642 // y rotation matrix 02643 Matrix3x3<T> Ry( c, 0, -s, 02644 0, 1, 0, 02645 s, 0, c ); 02646 // x rotation matrix 02647 T cx = cos(count); 02648 T sx = sin(count); 02649 count += countStep; 02650 Matrix3x3<T> Rx( 1, 0, 0, 02651 0, cx, sx, 02652 0, -sx, cx ); 02653 v = Rx * Ry * v; 02654 m_prVertex[i] = m_prVertex[i-1] - v; 02655 ++i; 02656 deg += step; 02657 if ( deg <= -90 || 90 <= deg ) { 02658 step *= -1; 02659 } 02660 } 02661 //------------------------------------------- 02662 //*/ 02663 02664 /* 02665 //------------------------------------------- 02666 // Setup strand as a straight line 02667 int i = 0; 02668 m_prVertex[i++] = startPosition; 02669 while ( i <= m_iNoLinks ) { 02670 m_prVertex[i] = m_prVertex[i-1]; 02671 m_prVertex[i][0] += m_tLinkLength; 02672 ++i; 02673 } 02674 //------------------------------------------- 02675 //*/ 02676 } 02677 //----------------------------------------------------------------------------- 02678 // SavePreviousPositions Fn 02679 template <typename T> 02680 void ModelStrand<T>::SavePreviousPositions () 02681 { 02682 std::cout << __FILE__ << ":" << __LINE__ << ": SavePreviousPositions()\n"; 02683 for ( int i = 0; i <= m_iNoLinks; ++i ) { 02684 m_prVertexPrevPos[i] = m_prVertex[i]; 02685 } 02686 } 02687 //----------------------------------------------------------------------------- 02688 // RestorePreviousPositions Fn 02689 template <typename T> 02690 void ModelStrand<T>::RestorePreviousPositions () 02691 { 02692 for ( int i = 0; i <= m_iNoLinks; ++i ) { 02693 m_prVertex[i] = m_prVertexPrevPos[i]; 02694 } 02695 UpdateBVHTree(); 02696 CheckSelfIntersections(); 02697 SetDisplayLinks(); 02698 } 02699 //----------------------------------------------------------------------------- 02700 02701 //============================================================================= 02702 // Calculate Forces 02703 /* 02704 //----------------------------------------------------------------------------- 02705 // Stretching and Compression Force 02706 template <typename T> 02707 Vector3<T> ModelStrand<T>::CalForceStretchingAndCompression ( int pos ) 02708 { 02709 CalLinkDirection(); 02710 //--------------------------------------------------------------- 02711 Vector3<T> force( 0, 0, 0 ); 02712 if ( pos <= 0 ) { 02713 } 02714 else { 02715 force -= (GetLinkCurrentLength( pos-1 ) - m_tLinkLength) * GetLinkUnitVector( pos-1 ); 02716 } 02717 if ( pos >= m_iNoLinks ) { 02718 } 02719 else { 02720 force += (GetLinkCurrentLength( pos ) - m_tLinkLength) * GetLinkUnitVector( pos ); 02721 } 02722 //--------------------------------------------------------------- 02723 return force * m_tKStretching; 02724 } 02725 //----------------------------------------------------------------------------- 02726 // Bending Force 02727 template <typename T> 02728 Vector3<T> ModelStrand<T>::CalForceBending ( int pos ) 02729 { 02730 Vector<T> u_p1 = GetLinkCurrentLength( pos+1 ); 02731 Vector<T> u = GetLinkCurrentLength( pos ); 02732 Vector<T> u_n1 = GetLinkCurrentLength( pos-1 ); 02733 Vector<T> u_n2 = GetLinkCurrentLength( pos-2 ); 02734 //--------------------------------------------------------------- 02735 Vector3<T> force( 0, 0, 0 ); 02736 //--------------------------------------------------------------- 02737 return force; 02738 } 02739 //----------------------------------------------------------------------------- 02740 // Twisting Force 02741 template <typename T> 02742 Vector3<T> ModelStrand<T>::CalForceTwisting ( int pos ) 02743 { 02744 //--------------------------------------------------------------- 02745 Vector3<T> force( 0, 0, 0 ); 02746 //--------------------------------------------------------------- 02747 return force; 02748 } 02749 //----------------------------------------------------------------------------- 02750 // Friction Force 02751 template <typename T> 02752 Vector3<T> ModelStrand<T>::CalForceFriction ( int pos ) 02753 { 02754 //--------------------------------------------------------------- 02755 Vector3<T> force( 0, 0, 0 ); 02756 if ( pos <= 0 ) { 02757 } 02758 else { 02759 force += m_prVelocity[pos-1]; 02760 } 02761 if ( pos >= m_iNoLinks ) { 02762 } 02763 else { 02764 force -= m_prVelocity[pos]; 02765 } 02766 //--------------------------------------------------------------- 02767 return force * m_tKFriction; 02768 } 02769 //----------------------------------------------------------------------------- 02770 // Contact Force 02771 template <typename T> 02772 Vector3<T> ModelStrand<T>::CalForceContact ( int pos ) 02773 { 02774 //--------------------------------------------------------------- 02775 Vector3<T> force( 0, 0, 0 ); 02776 //--------------------------------------------------------------- 02777 return force; 02778 } 02779 //----------------------------------------------------------------------------- 02780 // Gravity Force 02781 template <typename T> 02782 Vector3<T> ModelStrand<T>::CalForceGravity ( int pos ) 02783 { 02784 //--------------------------------------------------------------- 02785 Vector3<T> force = - ( m_tGravity * CGMath<T>::VectorY ) * m_tPointWeight; 02786 //--------------------------------------------------------------- 02787 return force; 02788 } 02789 //----------------------------------------------------------------------------- 02790 // Get Link Vector 02791 template <typename T> 02792 Vector3<T> ModelStrand<T>::GetLinkUnitVector ( int linkNo ) 02793 { 02794 return (m_prVertex[linkNo+1] - m_prVertex[linkNo]).Normalized(); 02795 } 02796 //----------------------------------------------------------------------------- 02797 // Get Link Current Length 02798 template <typename T> 02799 T ModelStrand<T>::GetLinkCurrentLength ( int linkNo ) 02800 { 02801 return (m_prVertex[linkNo+1] - m_prVertex[linkNo]).Length(); 02802 } 02803 //*/ 02804 02805 //----------------------------------------------------------------------------- 02806 // SimByForce -- Calculate the new total internal force of each link 02807 template <typename T> 02808 void ModelStrand<T>::SimByForce ( T tCurrent, T tNext ) 02809 { 02810 T dt = tNext - tCurrent; // in milli seconds 02811 //------------------------------------------------------------------ 02812 // FOR TESTING/DEBUGGING 02813 //m_prIsFixed[0] = true; 02814 //m_prIsFixed[m_iNoLinks/2] = true; 02815 //m_prIsFixed[m_iNoLinks] = true; 02816 //------------------------------------------------------------------ 02817 T length; 02818 Vector3<T> force; 02819 bool useLinearSpring = false; 02820 02821 //* 02822 //---------------------------------------------------------------- 02823 // Add Stretching and Compression Force 02824 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02825 // For segment# 0 02826 { 02827 unsigned int s = 0; 02828 if ( !m_Segments[s].isLocked ) { 02829 for ( unsigned int i = m_Segments[s].start; i < m_Segments[s].end; ++i ) 02830 { 02831 m_prLinkDirection[i] = m_prVertex[i+1] - m_prVertex[i]; 02832 length = m_prLinkDirection[i].Length() - m_tLinkLength; 02833 02834 //------------------------------- 02835 // Calculate force 02836 if ( useLinearSpring ) { // Linear spring 02837 force = m_tKStretching * (length * m_prLinkDirection[i].Normalized()); 02838 } 02839 else { // Nonlinear spring 02840 T tKStretching = m_tKStretching; 02841 T absLength = fabs( length ); 02842 if ( absLength <= 0.50 ) { 02843 } 02844 else if ( absLength <= 0.75 ) { 02845 tKStretching *= 1.25; 02846 } 02847 else if ( absLength <= 0.90 ) { 02848 tKStretching *= 2.0; 02849 } 02850 else if ( absLength <= 1.50 ) { 02851 tKStretching *= 4.0; 02852 } 02853 else { 02854 tKStretching *= 8.0; 02855 } 02856 force = tKStretching * (length * m_prLinkDirection[i].Normalized()); 02857 } 02858 02859 //------------------------------- 02860 // Accumulate force 02861 if ( !m_prIsFixed[i] && !m_prIsFixed[i+1] ) { 02862 m_prForce[i] += force; 02863 m_prForce[i+1] -= force; 02864 } 02865 else if ( m_prIsFixed[i] && !m_prIsFixed[i+1] ) { 02866 m_prForce[i+1] -= force * 2; 02867 } 02868 else if ( !m_prIsFixed[i] && m_prIsFixed[i+1] ) { 02869 m_prForce[i] += force * 2; 02870 } 02871 } 02872 } 02873 } 02874 // For segment# > 0 02875 for ( unsigned int s = 1; s < m_Segments.size(); ++s ) { 02876 if ( !m_Segments[s].isLocked ) { 02877 for ( unsigned int i = m_Segments[s].start; i < m_Segments[s].end; ++i ) 02878 02879 #else //TAPs_STRAND_LOCK_SEGMENTS 02880 { 02881 { 02882 for ( int i = 0; i < m_iNoLinks; ++i ) 02883 #endif//TAPs_STRAND_LOCK_SEGMENTS 02884 { 02885 m_prLinkDirection[i] = m_prVertex[i+1] - m_prVertex[i]; 02886 length = m_prLinkDirection[i].Length() - m_tLinkLength; 02887 02888 //------------------------------- 02889 // Calculate force 02890 if ( useLinearSpring ) { // Linear spring 02891 force = m_tKStretching * (length * m_prLinkDirection[i].Normalized()); 02892 } 02893 else { // Nonlinear spring 02894 T tKStretching = m_tKStretching; 02895 T absLength = fabs( length ); 02896 if ( absLength <= 0.50 ) { 02897 } 02898 else if ( absLength <= 0.75 ) { 02899 tKStretching *= 1.25; 02900 } 02901 else if ( absLength <= 0.90 ) { 02902 tKStretching *= 2.0; 02903 } 02904 else if ( absLength <= 1.50 ) { 02905 tKStretching *= 4.0; 02906 } 02907 else { 02908 tKStretching *= 8.0; 02909 } 02910 force = tKStretching * (length * m_prLinkDirection[i].Normalized()); 02911 } 02912 02913 //------------------------------- 02914 // Accumulate force 02915 if ( !m_prIsFixed[i] && !m_prIsFixed[i+1] ) { 02916 m_prForce[i] += force; 02917 m_prForce[i+1] -= force; 02918 } 02919 else if ( m_prIsFixed[i] && !m_prIsFixed[i+1] ) { 02920 m_prForce[i+1] -= force * 2; 02921 } 02922 else if ( !m_prIsFixed[i] && m_prIsFixed[i+1] ) { 02923 m_prForce[i] += force * 2; 02924 } 02925 } 02926 } 02927 } 02928 //*/ 02929 02930 //* 02931 //---------------------------------------------------------------- 02932 // Add Bending Force 02933 Vector3<T> bendingForceA, bendingForceB; 02934 T bendingAngleA, bendingAngleB; 02935 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02936 // For segment# 0 02937 { 02938 unsigned int s = 0; 02939 if ( !m_Segments[s].isLocked ) { 02940 for ( unsigned int i = m_Segments[s].start+2; i <= m_Segments[s].end-2; ++i ) 02941 { 02942 if ( !m_prIsFixed[i] ) { 02943 bendingAngleA = Math<T>::PI - 02944 Math<T>::ACos( m_prLinkDirection[i-1] * (-m_prLinkDirection[i-2]) ); 02945 bendingAngleB = Math<T>::PI - 02946 Math<T>::ACos( m_prLinkDirection[i+1] * (-m_prLinkDirection[ i ]) ); 02947 bendingForceA = bendingAngleA * m_prLinkDirection[i-1] ^ (m_prLinkDirection[i-2]^m_prLinkDirection[i-1]); 02948 bendingForceB = bendingAngleB * m_prLinkDirection[ i ] ^ (m_prLinkDirection[ i ]^m_prLinkDirection[i+1]); 02949 force = m_tKBending * ( bendingForceA + bendingForceB ); 02950 02951 m_prForce[i] += force; 02952 } 02953 } 02954 } 02955 } 02956 // For segment# > 0 02957 for ( unsigned int s = 1; s < m_Segments.size(); ++s ) { 02958 if ( !m_Segments[s].isLocked ) { 02959 for ( unsigned int i = m_Segments[s].start+2; i <= m_Segments[s].end-2; ++i ) 02960 02961 #else //TAPs_STRAND_LOCK_SEGMENTS 02962 { 02963 { 02964 for ( int i = 2; i <= m_iNoLinks-2; ++i ) 02965 #endif//TAPs_STRAND_LOCK_SEGMENTS 02966 { 02967 if ( !m_prIsFixed[i] ) { 02968 bendingAngleA = Math<T>::PI - 02969 Math<T>::ACos( m_prLinkDirection[i-1] * (-m_prLinkDirection[i-2]) ); 02970 bendingAngleB = Math<T>::PI - 02971 Math<T>::ACos( m_prLinkDirection[i+1] * (-m_prLinkDirection[ i ]) ); 02972 bendingForceA = bendingAngleA * m_prLinkDirection[i-1] ^ (m_prLinkDirection[i-2]^m_prLinkDirection[i-1]); 02973 bendingForceB = bendingAngleB * m_prLinkDirection[ i ] ^ (m_prLinkDirection[ i ]^m_prLinkDirection[i+1]); 02974 force = m_tKBending * ( bendingForceA + bendingForceB ); 02975 02976 m_prForce[i] += force; 02977 } 02978 } 02979 } 02980 } 02981 //*/ 02982 02983 //* 02984 //---------------------------------------------------------------- 02985 // Add Friction Force and Gravity Force 02986 m_tGravity = m_tKGravity * dt; 02987 #ifdef TAPs_STRAND_LOCK_SEGMENTS 02988 // For segment# 0 02989 { 02990 unsigned int s = 0; 02991 if ( !m_Segments[s].isLocked ) { 02992 for ( unsigned int i = m_Segments[s].start; i <= m_Segments[s].end; ++i ) 02993 { 02994 if ( !m_prIsFixed[i] ) { 02995 //------------------------------------------------- 02996 // Gravity 02997 //if ( m_prVertex[i][1] > 2 * m_tRadius ) { 02998 // m_prForce[i] -= ( m_tGravity * CGMath<T>::VectorY ) * m_tPointWeight; 02999 //} 03000 //else { 03001 // m_prForce[i][1] = 0; 03002 // m_prVelocity[i][1] = 0; 03003 // m_prVertex[i][1] = 0; 03004 //} 03005 //------------------------------------------------- 03006 // Friction 03007 //Vector3<T> velocity = m_prForce[i] * dt / m_tPointWeight; 03008 //Vector3<T> velocity = m_prVertex[i] - m_prVertexPrevPos[i]; 03009 if ( i == 0 || i == m_iNoLinks ) { 03010 m_prForce[i] -= m_tKFriction * m_prVelocity[i] * 2.0; 03011 //m_prForce[i] -= m_tKFriction * m_prVelocity[i]; 03012 03013 //m_prForce[i] -= m_tKFriction * velocity * 0.5; 03014 } 03015 else { 03016 m_prForce[i] -= m_tKFriction * m_prVelocity[i]; 03017 //m_prForce[i] -= m_tKFriction * velocity; 03018 03019 //m_prForce[i] -= m_tKFriction * velocity; 03020 } 03021 } 03022 } 03023 } 03024 } 03025 // For segment# > 0 03026 for ( unsigned int s = 1; s < m_Segments.size(); ++s ) { 03027 if ( !m_Segments[s].isLocked ) { 03028 for ( unsigned int i = m_Segments[s].start; i <= m_Segments[s].end; ++i ) 03029 #else //TAPs_STRAND_LOCK_SEGMENTS 03030 { 03031 { 03032 for ( int i = 0; i <= m_iNoLinks; ++i ) 03033 #endif//TAPs_STRAND_LOCK_SEGMENTS 03034 { 03035 if ( !m_prIsFixed[i] ) { 03036 //------------------------------------------------- 03037 // Gravity 03038 //if ( m_prVertex[i][1] > 2 * m_tRadius ) { 03039 // m_prForce[i] -= ( m_tGravity * CGMath<T>::VectorY ) * m_tPointWeight; 03040 //} 03041 //else { 03042 // m_prForce[i][1] = 0; 03043 // m_prVelocity[i][1] = 0; 03044 // m_prVertex[i][1] = 0; 03045 //} 03046 //------------------------------------------------- 03047 // Friction 03048 //Vector3<T> velocity = m_prForce[i] * dt / m_tPointWeight; 03049 //Vector3<T> velocity = m_prVertex[i] - m_prVertexPrevPos[i]; 03050 if ( i == 0 || i == m_iNoLinks ) { 03051 m_prForce[i] -= m_tKFriction * m_prVelocity[i] * 2.0; 03052 //m_prForce[i] -= m_tKFriction * m_prVelocity[i]; 03053 03054 //m_prForce[i] -= m_tKFriction * velocity * 0.5; 03055 } 03056 else { 03057 m_prForce[i] -= m_tKFriction * m_prVelocity[i]; 03058 //m_prForce[i] -= m_tKFriction * velocity; 03059 03060 //m_prForce[i] -= m_tKFriction * velocity; 03061 } 03062 } 03063 } 03064 } 03065 } 03066 //*/ 03067 } 03068 //----------------------------------------------------------------------------- 03069 //============================================================================= 03070 03071 //----------------------------------------------------------------------------- 03072 // Collision Detection with a multi bounding object 03073 template <typename T> 03074 void ModelStrand<T>::CDAndRForStrandWithAnMBVObj ( 03075 int iLeader, int iFollower, 03076 MultiBoundingVolume<T> const * const pMBV, 03077 TransformationSupport<T> const * const pTransform = NULL 03078 ) 03079 { 03080 //--------------------------------------------------------------- 03081 //bool isThePointInsideTheBV; // the pt is inside or outside the bv 03082 Vector3<T> vTotalDistance; // total displacement vector 03083 Vector3<T> vDistance; // displacement vector 03084 TransformationSupport<T> transform; 03085 //--------------------------------------------------------------- 03086 std::vector< BoundingVolume<T> * >::const_iterator pos = 03087 pMBV->GetBoundingVolumeList().begin(); 03088 //--------------------------------------------------------------- 03089 while ( pos != pMBV->GetBoundingVolumeList().end() ) { 03090 assert( *pos != NULL ); 03091 //isThePointInsideTheBV = false; 03092 transform.MakeIdentity(); 03093 //------------------------------------------------------ 03094 if ( pTransform ) { 03095 transform.ReturnMatrixTransform() *= pTransform->GetMatrixTransform(); 03096 } 03097 transform.ReturnMatrixTransform() *= pMBV->GetTransform().GetMatrixTransform(); 03098 transform.ReturnMatrixTransform() *= (*pos)->GetTransform().GetMatrixTransform(); 03099 //------------------------------------------------------ 03100 // Transform each vertex of the strand to the bounding volume coordinate 03101 //for ( int i = 0; i <= gModelManager->strandObj->GetNumberOfStrandLinksWithoutNeedle(); ++i ) { 03102 if ( m_iLinkKnotGrp[i] < 0 ) // if this link doesn't belong to a knot group 03103 for ( int i = iFollower; i <= iFollower; ++i ) { 03104 vDistance.SetXYZ( 0, 0, 0 ); 03105 //------------------------------------------------------ 03106 // In the local coordinate of the bounding cylinder, 03107 // the cylinder axis is parallel to z-axis 03108 // where it is shifted by the cylinder center. 03109 // The cylinder has radius, height and center. 03110 // Its height is measured from -z and +z axis 03111 // where its center is situated between height/2 in -z and +z axis. 03112 Vector3<T> location = m_prVertex[i]; 03113 location -= transform.GetTranslation(); 03114 location = transform.GetMatrixRotation().GetInverse() * location; 03115 // Shift the location of point by the cylinder center 03116 // i.e. shift the cylinder to the origin (0,0,0) 03117 location -= (*pos)->GetCenter(); 03118 //------------------------------------------------------ 03119 // Now the cylinder is centered at the origin with its axis on the z-axis 03120 // where its height is from -height/2 to +height/2. 03121 // And the input point has been transformed into the cylinder coordinate. 03122 // The test can be performed below 03123 T halfHeight = (*pos)->GetHeight()/2.0; 03124 T ptRadius = sqrt( location[0]*location[0] + location[1]*location[1] ); 03125 //====================================================== 03126 // If the point is INSIDE the cylinder 03127 //------------------------------------------------------ 03128 // If the point is within the cylinder radius 03129 if ( ptRadius <= (*pos)->GetRadius() ) { 03130 // If the point is within the cylinder height 03131 if ( -halfHeight <= location[2] && location[2] <= halfHeight ) { 03132 // Find the shortest distance vector from the point to the surface 03133 // of the cylinder 03134 T difRadius = (*pos)->GetRadius() - ptRadius; 03135 T difHeight = halfHeight - fabs( location[2] ); 03136 if ( difRadius <= difHeight ) { 03137 // If the point is on the z-axis (degenerate case) 03138 if ( difRadius >= (*pos)->GetRadius() ) { 03139 vDistance.SetXYZ( 0, (*pos)->GetRadius(), 0 ); 03140 } 03141 // If the point is NOT on the z-axis 03142 else { 03143 Vector2<T> difDirection = Vector2<T>( location[0], location[1] ); 03144 difDirection.Normalized(); 03145 difDirection *= difRadius; 03146 vDistance.SetXYZ( difDirection[0], difDirection[1], 0 ); 03147 } 03148 } 03149 else { 03150 vDistance.SetXYZ( 0, 0, location[2] >= 0 ? difHeight : -difHeight ); 03151 } 03152 //------------------------- 03153 // Now rotate it back by the rotation matrix (from the transform matrix). 03154 // REMARK: Translation (from the transform matrix) is NOT applied back, 03155 // because vDistance represents a displacement vector from 03156 // the input point to the cylinder surface. 03157 vDistance = transform.GetMatrixRotation() * vDistance; 03158 //------------------------- 03159 //isThePointInsideTheBV = true; 03160 //------------------------- 03161 Vector3<T> vOffset = vDistance.GetUnit() * m_tRadius; 03162 vTotalDistance += vDistance + vOffset; 03163 //------------------------- 03164 } 03165 vTotalDistance = vTotalDistance + m_prVertex[i] - m_prVertex[iLeader]; 03166 m_prVertex[i] = m_prVertex[iLeader] + vTotalDistance.GetUnit() * m_tLinkLength; 03167 } 03168 //--------------------------------------------- 03169 } // Next Strand Point 03170 //------------------------------------------------------ 03171 ++pos; 03172 } 03173 } 03174 //----------------------------------------------------------------------------- 03175 03176 //============================================================================= 03177 #ifdef TAPs_STRAND_LOCK_SEGMENTS 03178 //----------------------------------------------------------------------------- 03179 template <typename T> 03180 int ModelStrand<T>::DivideAtPt ( unsigned int pt ) 03181 { 03182 /* 03183 // WORK IN VC++ 2003 (I.E., VC++ .NET) 03184 // Find the segment containing the point 03185 for ( unsigned int i = 0; i < static_cast<unsigned int>( m_Segments.size() ); ++i ) { 03186 // Found the segment containin the point 03187 if ( m_Segments[i].start < pt && pt < m_Segments[i].end ) { 03188 // Devide the old segment into two segments 03189 // The old one becomes the right segment and its ID is shifted by +1 03190 unsigned int newStart = m_Segments[i].start; 03191 // 03192 // &m_Segments[i] will not be converted into iterator in VC++ later version than 2003 03193 m_Segments.insert( &m_Segments[i], IC_Segment( newStart, pt, false, true ) ); 03194 03195 m_Segments[i+1].start = pt; 03196 03197 // Fix the divided point 03198 SetFixStatusOfPtNo( pt, true ); 03199 03200 return i; 03201 } 03202 } 03203 //*/ 03204 03205 // Find the segment containing the point 03206 std::vector< IC_Segment >::iterator pos = m_Segments.begin(); 03207 int i = 0; 03208 while ( pos != m_Segments.end() ) { 03209 // Found the segment containin the point 03210 if ( pos->start < pt && pt < pos->end ) { 03211 // Devide the old segment into two segments 03212 // The old one becomes the right segment and its ID is shifted by +1 03213 unsigned int newStart = pos->start; 03214 // 03215 m_Segments.insert( pos, IC_Segment( newStart, pt, false, true ) ); 03216 03217 // DOESN'T WORK!!! WHY??? 03218 //++pos; 03219 //pos->start = pt; 03220 03221 // USE THIS INSTEAD 03222 m_Segments[i+1].start = pt; 03223 03224 // Fix the divided point 03225 SetFixStatusOfPtNo( pt, true ); 03226 03227 return i; 03228 } 03229 ++i; 03230 ++pos; 03231 } 03232 03233 return -1; // cannot divide the strand at the pt 03234 } 03235 //----------------------------------------------------------------------------- 03236 template <typename T> 03237 int ModelStrand<T>::GetStartPtOfSegment ( unsigned int segment ) 03238 { 03239 if ( !ValidateSegment( segment ) ) return -1; 03240 return m_Segments[segment].start; 03241 } 03242 //----------------------------------------------------------------------------- 03243 template <typename T> 03244 int ModelStrand<T>::GetEndPtOfSegment ( unsigned int segment ) 03245 { 03246 if ( !ValidateSegment( segment ) ) return -1; 03247 return m_Segments[segment].end; 03248 } 03249 //----------------------------------------------------------------------------- 03250 template <typename T> 03251 bool ModelStrand<T>::GetSegmentLockStatus ( unsigned int segment ) 03252 { 03253 if ( !ValidateSegment( segment ) ) return false; 03254 return m_Segments[segment].isLocked; 03255 } 03256 //----------------------------------------------------------------------------- 03257 template <typename T> 03258 bool ModelStrand<T>::GetSegmentVisibleStatus ( unsigned int segment ) 03259 { 03260 if ( !ValidateSegment( segment ) ) return false; 03261 return m_Segments[segment].isVisible; 03262 } 03263 //----------------------------------------------------------------------------- 03264 template <typename T> 03265 void ModelStrand<T>::SetSegmentLockStatus ( unsigned int segment, bool status ) 03266 { 03267 if ( !ValidateSegment( segment ) ) return; 03268 m_Segments[segment].isLocked = status; 03269 03270 #ifdef TAPs_ADVANCED_SIMULATION 03271 if ( status ) { 03272 for ( unsigned int i = m_Segments[segment].start; i < m_Segments[segment].end; ++i ) { 03273 m_prIsLocked[i] = status; 03274 SetFixStatusOfPtNo( i, status ); 03275 m_SimFlags[i].SetSimulationConstraints( TAPs::Enum::AddOn::FIXED ); 03276 } 03277 } 03278 else { 03279 for ( unsigned int i = m_Segments[segment].start; i < m_Segments[segment].end; ++i ) { 03280 m_prIsLocked[i] = status; 03281 SetFixStatusOfPtNo( i, status ); 03282 m_SimFlags[i].ClearSimulationConstraints( TAPs::Enum::AddOn::FIXED ); 03283 } 03284 } 03285 #else //TAPs_ADVANCED_SIMULATION 03286 for ( unsigned int i = m_Segments[segment].start; i < m_Segments[segment].end; ++i ) { 03287 m_prIsLocked[i] = status; 03288 SetFixStatusOfPtNo( i, status ); 03289 } 03290 #endif//TAPs_ADVANCED_SIMULATION 03291 03292 /* 03293 if ( status == true ) { 03294 03295 // Tighten the first and the last link of the lock segment 03296 // I.e. make each one's length the rest length if its length is greater than the rest length 03297 03298 // The start of the segment 03299 T length = m_prVertex[m_Segments[segment].start].Length(); 03300 if ( length > m_tLinkLength ) { 03301 m_prVertex[m_Segments[segment].start] = 03302 (m_prVertex[m_Segments[segment].start] - m_prVertex[m_Segments[segment].start + 1]) * length * 2//m_tLinkLength/length 03303 + m_prVertex[m_Segments[segment].start + 1]; 03304 } 03305 03306 // The end of the segment 03307 length = m_prVertex[m_Segments[segment].end-1].Length(); 03308 if ( length > m_tLinkLength ) { 03309 m_prVertex[m_Segments[segment].end] = 03310 (m_prVertex[m_Segments[segment].end] - m_prVertex[m_Segments[segment].end - 1]) * m_tLinkLength/length 03311 + m_prVertex[m_Segments[segment].end - 1]; 03312 } 03313 03314 // DEBUG 03315 std::cout << "LOCK segment#" << segment << " [start,end]: [" << m_Segments[segment].start << "," << m_Segments[segment].end << "]\n"; 03316 } 03317 //*/ 03318 } 03319 //----------------------------------------------------------------------------- 03320 template <typename T> 03321 void ModelStrand<T>::SetSegmentVisibleStatus ( unsigned int segment, bool status ) 03322 { 03323 if ( !ValidateSegment( segment ) ) return; 03324 m_Segments[segment].isVisible = status; 03325 } 03326 //----------------------------------------------------------------------------- 03327 #endif//TAPs_STRAND_LOCK_SEGMENTS 03328 //============================================================================= 03329 03330 //============================================================================= 03331 #ifdef TAPs_STRAND_KNOT_CONTROL 03332 //----------------------------------------------------------------------------- 03333 template <typename T> 03334 void ModelStrand<T>::CreateKnotControl ( unsigned char numKnotAllowed ) 03335 { 03336 m_ucMaxNumKnotsAllowed = numKnotAllowed; 03337 m_piKnotID = new int[ m_ucMaxNumKnotsAllowed ]; 03338 ResetKnotControl(); 03339 //=============================================================== 03340 #ifdef TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03341 //--------------------------------------------------------------- 03342 m_WaitingTime = 5.0; // set waiting time 03343 03344 m_uiCounterThreshold = 20; // act as a timer for fixing the current knot 03345 03346 m_uiMinLengthThreshold = 0; // start at zero 03347 03348 m_uiLengthThreshold[0] = 12; // length threshold for half knot -- basic knot 03349 // extra length threshold for double knot will be added by m_uiLengthThreshold[0]/2 03350 m_uiLengthThreshold[1] = 7; // extra length threshold for double knot 03351 // i.e. for double knot the length threshold is 18 (of half knot) + 8 03352 03353 // Unused for now 03354 m_uiLengthThreshold[2] = 8; // length threshold for ... 03355 m_uiLengthThreshold[3] = 8; // length threshold for ... 03356 m_uiLengthThreshold[4] = 8; // length threshold for ... 03357 m_uiLengthThreshold[5] = 8; // length threshold for ... 03358 m_uiLengthThreshold[6] = 8; // length threshold for ... 03359 m_uiLengthThreshold[7] = 8; // length threshold for ... 03360 //--------------------------------------------------------------- 03361 #endif//TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03362 //=============================================================== 03363 } 03364 //----------------------------------------------------------------------------- 03365 template <typename T> 03366 void ModelStrand<T>::DeleteKnotControl () 03367 { 03368 if ( m_piKnotID ) { 03369 delete [] m_piKnotID; 03370 m_piKnotID = NULL; 03371 } 03372 } 03373 //----------------------------------------------------------------------------- 03374 template <typename T> 03375 void ModelStrand<T>::InitKnotControl ( unsigned char numKnotAllowed ) 03376 { 03377 DeleteKnotControl(); 03378 CreateKnotControl( numKnotAllowed ); 03379 } 03380 //----------------------------------------------------------------------------- 03381 template <typename T> 03382 void ModelStrand<T>::ResetKnotControl () 03383 { 03384 m_ucKnotCount = 0; 03385 for ( int i = 0; i < m_ucMaxNumKnotsAllowed; ++i ) { 03386 m_piKnotID[i] = -1; 03387 } 03388 //=============================================================== 03389 #ifdef TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03390 //--------------------------------------------------------------- 03391 m_uiCounterForFixingKnot = 0; // reset counter 03392 m_uiMinLengthThreshold = 0; 03393 ResetElapsedTime(); 03394 m_uiLengthAccumulateTarget = m_uiLengthThreshold[0]; // reset accumulated target length 03395 m_ucTrackingKnot = Data::DatabaseDowkerNotation::UNDEFINED; // reset tracker 03396 //--------------------------------------------------------------- 03397 #endif//TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03398 //=============================================================== 03399 } 03400 //----------------------------------------------------------------------------- 03401 template <typename T> 03402 std::string ModelStrand<T>::CheckSurgeonsKnot () 03403 { 03404 //=============================================================== 03405 #ifdef TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03406 //--------------------------------------------------------------- 03407 std::string str; 03408 static std::string name; 03409 03410 // Recognize the knot 03411 //* 03412 if ( m_ucTrackingKnot == Data::DatabaseDowkerNotation::UNDEFINED ) { 03413 03414 // Check if a new knot is being formed 03415 m_ucTrackingKnot = DowkerNotationToKnotName( &name ); 03416 03417 // For current update, CheckSelfCrossingForKnotRecognition fn should be before DowkerNotationToKnotName fn 03418 #ifdef TAPs_STRAND_DEBUG 03419 special_debug_knotID = 03420 #endif//TAPs_STRAND_DEBUG 03421 CheckSelfCrossingForKnotRecognition( GetProjectionDirection( m_prVertex[0], m_prVertex[m_iNoLinks/2], m_prVertex[m_iNoLinks] ) ); 03422 int start = GetKnotStartPt(); 03423 int end = GetKnotEndPt(); 03424 03425 // Recognize the knot -- Total Knots is One 03426 if ( m_ucKnotCount == 1 ) { 03427 str = GetNameOfKnot( 0, false ); // 2nd parameter: true with (-)/(+) for direction 03428 str += " Knot"; 03429 //str = "One Knot"; 03430 } 03431 03432 // Recognize the knot -- Total Knots is Two 03433 else if ( m_ucKnotCount == 2 ) { 03434 //----------------------------------- 03435 // Check Surgeon's Knot (NEG-then-POS or POS-then-NEG) 03436 if ( 03437 // neg double then pos half equals surgeon's knot 03438 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_NEG 03439 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03440 || 03441 // pos double then neg half equals surgeon's knot 03442 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_POS 03443 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03444 ) { 03445 str = "Surgeon's Knot"; 03446 } 03447 //----------------------------------- 03448 // Check Surgeon's Knot (NEG-then-NEG or POS-then-POS) 03449 else if ( 03450 // neg double then pos half equals surgeon's knot 03451 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_NEG 03452 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03453 || 03454 // pos double then neg half equals surgeon's knot 03455 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_POS 03456 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03457 ) { 03458 str = "Granny Double Knot"; 03459 } 03460 //----------------------------------- 03461 // Check Square Knot (NEG-then-POS or POS-then-NEG) 03462 else if ( 03463 // neg half then pos half equals square knot 03464 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_NEG 03465 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03466 || 03467 // pos half then neg half equals square knot 03468 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_POS 03469 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03470 ) { 03471 str = "Square Knot"; 03472 } 03473 //----------------------------------- 03474 // Check Square Knot (NEG-then-NEG or POS-then-POS) 03475 else if ( 03476 // neg half then pos half equals square knot 03477 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_NEG 03478 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03479 || 03480 // pos half then neg half equals square knot 03481 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_POS 03482 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03483 ) { 03484 str = "Granny Knot"; 03485 } 03486 //----------------------------------- 03487 // Unrecognized Knot Combination 03488 else { 03489 char s[16]; 03490 str = _itoa( m_ucKnotCount, s, 10 ); 03491 str += " knot combination: "; 03492 03493 //* 03494 #ifdef TAPs_STRAND_DEBUG 03495 for ( int i = 0; i < m_ucKnotCount; ++i ) { 03496 //str += itoa( m_piKnotID[i], s, 10 ); 03497 str += GetNameOfKnot( i, true ); // 2nd parameter: true with (-)/(+) for direction 03498 if ( i < m_ucKnotCount-1 ) str += ", "; 03499 } 03500 #endif//TAPs_STRAND_DEBUG 03501 //*/ 03502 str += GetStrNameForKnotCombination(); 03503 } 03504 } 03505 03506 // Recognize the knot -- Total Knots is Three 03507 else if ( m_ucKnotCount == 3 ) { 03508 //----------------------------------- 03509 // Check Surgeon's Knot (NEG-then-POS-then-NEG or POS-then-NEG-then-POS) 03510 if ( 03511 // neg double then pos half equals surgeon's knot 03512 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_NEG 03513 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS 03514 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03515 || 03516 // pos double then neg half equals surgeon's knot 03517 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_POS 03518 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG 03519 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03520 ) { 03521 str = "Square Surgeon's Knot"; 03522 } 03523 //----------------------------------- 03524 // Check Surgeon's Knot (NEG-then-NEG-then-NEG or POS-then-POS-then-POS) 03525 else if ( 03526 // neg double then pos half equals surgeon's knot 03527 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_NEG 03528 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG 03529 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03530 || 03531 // pos double then neg half equals surgeon's knot 03532 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::DOUBLE_POS 03533 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS 03534 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03535 ) { 03536 str = "Granny Granny Double Knot"; 03537 } 03538 //----------------------------------- 03539 // Check Square's Knot (NEG-then-POS-then-NEG or POS-then-NEG-then-POS) 03540 else if ( 03541 // neg double then pos half equals surgeon's knot 03542 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_NEG 03543 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS 03544 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03545 || 03546 // pos double then neg half equals surgeon's knot 03547 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_POS 03548 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG 03549 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03550 ) { 03551 str = "Square Square Knot"; 03552 } 03553 //----------------------------------- 03554 // Check Square's Knot (NEG-then-NEG-then-NEG or POS-then-POS-then-POS) 03555 else if ( 03556 // neg double then pos half equals surgeon's knot 03557 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_NEG 03558 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG 03559 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_NEG ) 03560 || 03561 // pos double then neg half equals surgeon's knot 03562 ( m_piKnotID[0] == Data::DatabaseDowkerNotation::HALF_POS 03563 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS 03564 && m_piKnotID[1] == Data::DatabaseDowkerNotation::HALF_POS ) 03565 ) { 03566 str = "Granny Granny Knot"; 03567 } 03568 //----------------------------------- 03569 // Unrecognized Knot Combination 03570 else { 03571 char s[16]; 03572 str = _itoa( m_ucKnotCount, s, 10 ); 03573 str += " knot combination: "; 03574 03575 //* 03576 #ifdef TAPs_STRAND_DEBUG 03577 for ( int i = 0; i < m_ucKnotCount; ++i ) { 03578 //str += itoa( m_piKnotID[i], s, 10 ); 03579 str += GetNameOfKnot( i, true ); // 2nd parameter: true with (-)/(+) for direction 03580 if ( i < m_ucKnotCount-1 ) str += ", "; 03581 } 03582 #endif//TAPs_STRAND_DEBUG 03583 //*/ 03584 str += GetStrNameForKnotCombination(); 03585 } 03586 } 03587 03588 // Recognize the knot -- Total Knots is greater than three 03589 else { 03590 char s[16]; 03591 str = _itoa( m_ucKnotCount, s, 10 ); 03592 str += " knot combination: "; 03593 03594 //* 03595 #ifdef TAPs_STRAND_DEBUG 03596 for ( int i = 0; i < m_ucKnotCount; ++i ) { 03597 //str += itoa( m_piKnotID[i], s, 10 ); 03598 str += GetNameOfKnot( i, true ); // 2nd parameter: true with (-)/(+) for direction 03599 if ( i < m_ucKnotCount-1 ) str += ", "; 03600 } 03601 #endif//TAPs_STRAND_DEBUG 03602 //*/ 03603 str += GetStrNameForKnotCombination(); 03604 } 03605 } 03606 //*/ 03607 03608 /* 03609 // DEBUG 03610 if ( m_ucTrackingKnot == Data::DatabaseDowkerNotation::UNDEFINED ) { 03611 // Check if a new knot is being formed 03612 m_ucTrackingKnot = DowkerNotationToKnotName( &name ); 03613 03614 // For current update, CheckSelfCrossingForKnotRecognition fn should be before DowkerNotationToKnotName fn 03615 #ifdef TAPs_STRAND_DEBUG 03616 special_debug_knotID = 03617 #endif//TAPs_STRAND_DEBUG 03618 Vector3<T> v; 03619 CheckSelfCrossingForKnotRecognition( GetProjectionDirection( StartPoint, MiddlePoint, EndPoint ) ); 03620 int start = GetKnotStartPt(); 03621 int end = GetKnotEndPt(); 03622 03623 char s[16]; 03624 str = itoa( m_ucKnotCount, s, 10 ); 03625 str += " knot combination: "; 03626 for ( int i = 0; i < m_ucKnotCount; ++i ) { 03627 str += itoa( m_piKnotID[i], s, 10 ); 03628 if ( i < m_ucKnotCount-1 ) str += ", "; 03629 } 03630 } 03631 //*/ 03632 03633 03634 // Track the current knot 03635 else { 03636 // For current update, CheckSelfCrossingForKnotRecognition fn should be before DowkerNotationToKnotName fn 03637 #ifdef TAPs_STRAND_DEBUG 03638 special_debug_knotID = 03639 #endif//TAPs_STRAND_DEBUG 03640 CheckSelfCrossingForKnotRecognition( GetProjectionDirection( m_prVertex[0], m_prVertex[m_iNoLinks/2], m_prVertex[m_iNoLinks] ) ); 03641 int start = GetKnotStartPt(); 03642 int end = GetKnotEndPt(); 03643 03644 //------------------------------------------------- 03645 // Change tracking knot -- for a more complex knot 03646 03647 //#ifdef TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 03648 //#else //TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 03649 // Current tracked knot is a half knot 03650 if ( m_ucTrackingKnot == Data::DatabaseDowkerNotation::HALF_NEG 03651 || m_ucTrackingKnot == Data::DatabaseDowkerNotation::HALF_POS ) 03652 { 03653 // Check the current knot 03654 std::string name2; 03655 int knot = DowkerNotationToKnotName( &name2 ); 03656 03657 // The current knot is not a half knot 03658 if ( knot != Data::DatabaseDowkerNotation::UNDEFINED ) { 03665 //if ( knot == Data::DatabaseDowkerNotation::HALF_NEG 03666 // || knot == Data::DatabaseDowkerNotation::HALF_POS ) 03667 //{ 03668 // // Change from half knot to half knot 03669 // str = name = name2; 03670 // m_ucTrackingKnot = knot; 03671 //} 03672 //else 03673 if ( knot == Data::DatabaseDowkerNotation::DOUBLE_NEG 03674 || knot == Data::DatabaseDowkerNotation::DOUBLE_POS ) 03675 { 03676 // Change from half knot to double knot 03677 str = name = name2; 03678 m_ucTrackingKnot = knot; 03679 // Add extra length threshold for the knot change 03680 m_uiLengthAccumulateTarget += m_uiLengthThreshold[0] / 2; 03681 } 03682 } 03683 } 03684 //#endif//TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 03685 03686 //------------------------------------------------- 03687 // A new knot is formed or being tracked 03688 if ( m_ucTrackingKnot != Data::DatabaseDowkerNotation::UNDEFINED ) { 03689 03690 /* 03691 #ifdef TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 03692 if ( IsClumped() ) { 03693 char val[16]; 03694 str = "Tracking Knot# "; 03695 str += itoa(m_ucKnotCount+1, val, 10); 03696 str += " (" + name + ")"; 03697 03698 str += " ID: "; 03699 str += itoa(m_ucTrackingKnot, val, 10); 03700 03701 int knotLength = end-start; 03702 03703 if ( m_ucKnotCount == 0 ) { 03704 if ( ++m_uiCounterForFixingKnot > m_uiCounterThreshold ) { 03705 // Divide strand into segments 03706 DivideAtPt( start ); 03707 int segmentID = DivideAtPt( end ); 03708 SetSegmentLockStatus( segmentID, true ); 03709 03710 // Record knot ID and increase knot count 03711 m_piKnotID[m_ucKnotCount++] = m_ucTrackingKnot; 03712 // Reset counter and set the new accumulated target length for next knot 03713 m_uiCounterForFixingKnot = 0; 03714 m_uiLengthAccumulateTarget = knotLength + m_uiLengthThreshold[0]; 03715 // Reset tracking knot 03716 m_ucTrackingKnot = Data::DatabaseDowkerNotation::UNDEFINED; 03717 str = ""; 03718 03719 // DEBUG 03720 //PrtSegments(); 03721 //std::cout << "Fix Knot#: " << m_ucKnotCount << "\n"; 03722 } 03723 } 03724 else if ( 0 < m_ucKnotCount && m_ucKnotCount < m_ucMaxNumKnotsAllowed ) { 03725 if ( ++m_uiCounterForFixingKnot > m_uiCounterThreshold ) { 03726 // Create left segment 03727 int segmentID1 = DivideAtPt( start ); 03728 SetSegmentLockStatus( segmentID1+1, true ); 03729 // Create right segment 03730 int segmentID2 = DivideAtPt( end ); 03731 SetSegmentLockStatus( segmentID2, true ); 03732 03733 // Record knot ID and increase knot count 03734 m_piKnotID[m_ucKnotCount++] = m_ucTrackingKnot; 03735 // Reset counter and set the new accumulated target length for next knot 03736 m_uiCounterForFixingKnot = 0; 03737 m_uiLengthAccumulateTarget = knotLength + m_uiLengthThreshold[m_ucKnotCount]; 03738 // Reset tracking knot 03739 m_ucTrackingKnot = Data::DatabaseDowkerNotation::UNDEFINED; 03740 str = ""; 03741 03742 // DEBUG 03743 //PrtSegments(); 03744 //std::cout << "Fix Knot#: " << m_ucKnotCount << "\n"; 03745 } 03746 } 03747 03748 return str; 03749 } 03750 #endif//TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 03751 //*/ 03752 03753 // Recovering from incorrect collision detection 03754 // Case when there was a knot then there is no crossing pairs 03755 if ( m_iNumCrossingPairs == 0 ) { 03756 m_ucTrackingKnot = Data::DatabaseDowkerNotation::UNDEFINED; 03757 str = "Recovering from an error!"; 03758 03759 m_uiCounterForFixingKnot = 0; 03760 ResetElapsedTime(); 03761 03762 // DEBUG 03763 //std::cout << "m_iNumCrossingPairs was ZERO!\n" << std::endl; 03764 } 03765 03766 char val[16]; 03767 str = "Tracking Knot# "; 03768 str += _itoa(m_ucKnotCount+1, val, 10); 03769 str += " (" + name + ")"; 03770 03771 str += " ID: "; 03772 str += _itoa(m_ucTrackingKnot, val, 10); 03773 03774 // DEBUG 03775 //std::cout << "end-start: " << end-start << " = " << end << " - " << start << "\n"; 03776 //std::cout << "TargetLength: " << m_uiLengthAccumulateTarget << "\n"; 03777 //std::cout << "CounterForFixingKnot: " << m_uiCounterForFixingKnot << "\n"; 03778 03779 //----------------------------------- 03780 // Skip invalid length 03781 if ( end < 0 || start < 0 ) { 03782 //++m_uiCounterForFixingKnot; 03783 } 03784 //----------------------------------- 03785 // Track the first knot 03786 else if ( m_ucKnotCount == 0 ) { 03787 int knotLength = end-start; 03788 if ( knotLength >= static_cast<int>( m_uiMinLengthThreshold ) || static_cast<int>( m_uiCounterForFixingKnot ) > m_uiCounterThreshold ) { 03789 if ( knotLength < static_cast<int>( m_uiLengthAccumulateTarget ) ) { 03790 // Reset counter and set new accumulated target length 03791 m_uiCounterForFixingKnot = 0; 03792 m_uiLengthAccumulateTarget = knotLength; 03793 ResetElapsedTime(); 03794 } 03795 else if ( ( knotLength <= static_cast<int>( m_uiLengthAccumulateTarget ) + 2 ) ) { 03796 if ( ++m_uiCounterForFixingKnot >= m_uiCounterThreshold ) { 03797 // Divide strand into segments 03798 DivideAtPt( start ); 03799 int segmentID = DivideAtPt( end ); 03800 SetSegmentLockStatus( segmentID, true ); 03801 03802 // Record knot ID and increase knot count 03803 m_piKnotID[m_ucKnotCount++] = m_ucTrackingKnot; 03804 // Reset counter and set the new accumulated target length for next knot 03805 m_uiCounterForFixingKnot = 0; 03806 ResetElapsedTime(); 03807 if ( knotLength >= static_cast<int>( m_uiMinLengthThreshold ) ) { 03808 m_uiMinLengthThreshold = knotLength; // new min length 03809 } 03810 m_uiLengthAccumulateTarget = m_uiMinLengthThreshold + m_uiLengthThreshold[m_ucKnotCount]; 03811 // Reset tracking knot 03812 m_ucTrackingKnot = Data::DatabaseDowkerNotation::UNDEFINED; 03813 str = ""; 03814 03815 // DEBUG 03816 //PrtSegments(); 03817 //std::cout << "Fix Knot#: " << m_ucKnotCount << "\n"; 03818 } 03819 } 03820 else { 03821 ++m_uiCounterForFixingKnot; 03822 AddElapsedTime(); 03823 } 03824 } 03825 } 03826 //----------------------------------- 03827 // Tracking the second, third, and so on 03828 else if ( 0 < m_ucKnotCount && m_ucKnotCount < m_ucMaxNumKnotsAllowed ) { 03829 int knotLength = end-start; 03830 if ( knotLength >= static_cast<int>( m_uiMinLengthThreshold ) || static_cast<int>( m_uiCounterForFixingKnot ) > m_uiCounterThreshold ) { 03831 if ( knotLength < static_cast<int>( m_uiLengthAccumulateTarget ) ) { 03832 // Reset counter and set new accumulated target length 03833 m_uiCounterForFixingKnot = 0; 03834 m_uiLengthAccumulateTarget = knotLength; 03835 ResetElapsedTime(); 03836 } 03837 else if ( ( knotLength <= static_cast<int>( m_uiLengthAccumulateTarget ) + 2 ) || 03838 03839 ( IsTimeLimitReached() && knotLength <= static_cast<int>( m_uiLengthAccumulateTarget ) + 4 ) ) { 03840 03841 //(m_uiCounterForFixingKnot >= 125) ) { 03842 03843 #ifdef TAPs_STRAND_DEBUG 03844 if ( IsTimeLimitReached() ) { 03845 std::cout << "Elapsed Time: " << m_AccTime << "\n"; 03846 } 03847 #endif//TAPs_STRAND_DEBUG 03848 03849 if ( ++m_uiCounterForFixingKnot >= m_uiCounterThreshold ) { 03850 // Create left segment 03851 int segmentID1 = DivideAtPt( start ); 03852 SetSegmentLockStatus( segmentID1+1, true ); 03853 // Create right segment 03854 int segmentID2 = DivideAtPt( end ); 03855 SetSegmentLockStatus( segmentID2, true ); 03856 03857 // Record knot ID and increase knot count 03858 m_piKnotID[m_ucKnotCount++] = m_ucTrackingKnot; 03859 // Reset counter and set the new accumulated target length for next knot 03860 m_uiCounterForFixingKnot = 0; 03861 ResetElapsedTime(); 03862 if ( knotLength >= static_cast<int>( m_uiMinLengthThreshold ) ) { 03863 m_uiMinLengthThreshold = knotLength; // new min length 03864 } 03865 m_uiLengthAccumulateTarget = m_uiMinLengthThreshold + m_uiLengthThreshold[m_ucKnotCount]; 03866 // Reset tracking knot 03867 m_ucTrackingKnot = Data::DatabaseDowkerNotation::UNDEFINED; 03868 str = ""; 03869 03870 // DEBUG 03871 //PrtSegments(); 03872 //std::cout << "Fix Knot#: " << m_ucKnotCount << "\n"; 03873 } 03874 } 03875 else { 03876 ++m_uiCounterForFixingKnot; 03877 AddElapsedTime(); 03878 } 03879 } 03880 } 03881 } 03882 //std::cout << "Acc Time: " << m_AccTime << "\n"; 03883 } 03884 ResetStartTime(); 03885 return str; 03886 03887 #else //TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03888 return ""; 03889 03890 //--------------------------------------------------------------- 03891 #endif//TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03892 //=============================================================== 03893 03894 } 03895 //----------------------------------------------------------------------------- 03896 template <typename T> 03897 std::string ModelStrand<T>::GetStrNameForKnotCombination () 03898 { 03899 std::string str; 03900 for ( int i = 1; i < m_ucKnotCount; ++i ) { 03901 if ( m_piKnotID[i-1] == 0 ) { 03902 if ( m_piKnotID[i] == 0 ) { 03903 str = "Granny " + str; 03904 } 03905 else if ( m_piKnotID[i] == 1 ) { 03906 str = "Square " + str; 03907 } 03908 } 03909 else if ( m_piKnotID[i-1] == 1 ) { 03910 if ( m_piKnotID[i] == 1 ) { 03911 str = "Granny " + str; 03912 } 03913 else if ( m_piKnotID[i] == 0 ) { 03914 str = "Square " + str; 03915 } 03916 } 03917 } 03918 return str + " Knot"; 03919 } 03920 03921 #ifdef TAPs_COLLECT_DATA 03922 //----------------------------------------------------------------------------- 03923 template <typename T> 03924 std::string ModelStrand<T>::GetNameOfKnot ( unsigned int i, bool withPosOrNegPostfix ) 03925 { 03926 if ( i >= m_ucKnotCount ) return "n\a"; 03927 if ( withPosOrNegPostfix ) { 03928 if ( m_piKnotID[i] % 2 == 0 ) { 03929 return Data::DatabaseDowkerNotation::GetKnotName( m_piKnotID[i]/2 ) + "(-)"; 03930 } 03931 else { 03932 return Data::DatabaseDowkerNotation::GetKnotName( m_piKnotID[i]/2 ) + "(+)"; 03933 } 03934 } 03935 return Data::DatabaseDowkerNotation::GetKnotName( m_piKnotID[i]/2 ); 03936 } 03937 //----------------------------------------------------------------------------- 03938 template <typename T> 03939 int ModelStrand<T>::GetIDOfKnot ( unsigned int i ) 03940 { 03941 if ( i >= m_ucKnotCount ) return -1; 03942 return m_piKnotID[i]; 03943 } 03944 //----------------------------------------------------------------------------- 03945 #endif//TAPs_COLLECT_DATA 03946 03947 03948 //=================================================================== 03949 #ifdef TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 03950 //------------------------------------------------------------------- 03951 template <typename T> 03952 bool ModelStrand<T>::SetKnotControlThreshold ( Threshold th, T val ) 03953 { 03954 if ( val < 0.0 ) return false; 03955 switch ( th ) { 03956 case THRESHOLD_FOR_FIXING_KNOT: 03957 m_uiCounterThreshold = static_cast< unsigned int >( val ); 03958 return true; 03959 case THRESHOLD_WAITING_TIME: 03960 m_WaitingTime = static_cast< double >( val ); 03961 return true; 03962 case LENGTH_THRESHOLD_0: 03963 m_uiLengthThreshold[0] = static_cast< unsigned int >( val ); 03964 return true; 03965 case LENGTH_THRESHOLD_1: 03966 m_uiLengthThreshold[1] = static_cast< unsigned int >( val ); 03967 return true; 03968 case LENGTH_THRESHOLD_2: 03969 m_uiLengthThreshold[2] = static_cast< unsigned int >( val ); 03970 return true; 03971 case LENGTH_THRESHOLD_3: 03972 m_uiLengthThreshold[3] = static_cast< unsigned int >( val ); 03973 return true; 03974 case LENGTH_THRESHOLD_4: 03975 m_uiLengthThreshold[4] = static_cast< unsigned int >( val ); 03976 return true; 03977 case LENGTH_THRESHOLD_5: 03978 m_uiLengthThreshold[5] = static_cast< unsigned int >( val ); 03979 return true; 03980 case LENGTH_THRESHOLD_6: 03981 m_uiLengthThreshold[6] = static_cast< unsigned int >( val ); 03982 return true; 03983 case LENGTH_THRESHOLD_7: 03984 m_uiLengthThreshold[7] = static_cast< unsigned int >( val ); 03985 return true; 03986 } 03987 return false; 03988 } 03989 //------------------------------------------------------------------- 03990 template <typename T> 03991 T ModelStrand<T>::GetKnotControlThreshold ( Threshold th ) const 03992 { 03993 switch ( th ) { 03994 case THRESHOLD_FOR_FIXING_KNOT: 03995 return static_cast< T >( m_uiCounterThreshold ); 03996 case THRESHOLD_WAITING_TIME: 03997 return static_cast< T >( m_WaitingTime ); 03998 case LENGTH_THRESHOLD_0: 03999 return static_cast< T >( m_uiLengthThreshold[0] ); 04000 case LENGTH_THRESHOLD_1: 04001 return static_cast< T >( m_uiLengthThreshold[1] ); 04002 case LENGTH_THRESHOLD_2: 04003 return static_cast< T >( m_uiLengthThreshold[2] ); 04004 case LENGTH_THRESHOLD_3: 04005 return static_cast< T >( m_uiLengthThreshold[3] ); 04006 case LENGTH_THRESHOLD_4: 04007 return static_cast< T >( m_uiLengthThreshold[4] ); 04008 case LENGTH_THRESHOLD_5: 04009 return static_cast< T >( m_uiLengthThreshold[5] ); 04010 case LENGTH_THRESHOLD_6: 04011 return static_cast< T >( m_uiLengthThreshold[6] ); 04012 case LENGTH_THRESHOLD_7: 04013 return static_cast< T >( m_uiLengthThreshold[7] ); 04014 } 04015 return -1; 04016 } 04017 //------------------------------------------------------------------- 04018 04019 04020 //=================================================================== 04021 #ifdef TAPs_ALLOW_SETTING_PROJECTION_DIRECTION 04022 //------------------------------------------------------------------- 04023 // Find projection direction from a set of three points 04024 template <typename T> 04025 Vector3<T> ModelStrand<T>::GetProjectionDirection ( 04026 Vector3<T> const & pt1, Vector3<T> const & pt2, Vector3<T> const & pt3 ) 04027 { 04028 return TAPs::CGMath<T>::NormalOfTriangle( pt1, pt2, pt3 ); 04029 } 04030 //------------------------------------------------------------------- 04031 // Check self-crossings 04032 template <typename T> 04033 int ModelStrand<T>::CheckSelfCrossingForKnotRecognition ( 04034 Vector3<T> const & projectionDirection ) 04035 { 04036 //----------------------------------------------------- 04037 // Reset 04038 for ( int i = 0; i < m_iMaxCrossingPairs; ++i ) m_iaDowkerNotation[i] = 0; 04039 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04040 m_crossings.clear(); 04041 04042 for ( int i = 0; i < m_iNoLinks; ++i ) m_iaCrossingCounters[i] = 0; 04043 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04044 for ( int i = 0; i < m_iNoLinks; ++i ) m_iaCrossingsForDowkerNotation[i] = 0; 04045 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04046 04047 //----------------------------------------------------- 04048 // Find Crossings 04049 CheckSelfCrossingsNextLevelRecursively( m_BVHTree->Root()->Child(0), 04050 m_BVHTree->Root()->Child(1), 04051 projectionDirection ); 04052 //----------------------------------------------------- 04053 // For Recording Dowker Notation 04054 std::vector<int> crosses; 04055 std::vector<int> crossesID; 04056 04057 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04058 std::list< IC_Crossing >::const_iterator it = m_crossings.begin(); 04059 while ( it != m_crossings.end() ) { 04060 // crossesID records the node number that has another node cross over or under it. 04061 // crosses records the crossing identified by node number that crosses over the crossesID. 04062 // E.g. real cross# 1 2 3 4 5 6 04063 // crossesID: 30 31 32 74 77 79 04064 // crosses: -74 77 -79 30 -31 32 04065 // results in the Dowker notation below 04066 // Dowker#: (1,4) (3,6) (5,2) 04067 // which is a half knot 04068 crossesID.push_back( (*it).id ); 04069 crosses.push_back( (*it).number ); 04070 ++it; 04071 } 04072 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04073 for ( int i = 0; i < m_iNoLinks; ++i ) { 04074 if ( m_iaCrossingsForDowkerNotation[i] != 0 ) { 04075 // crossesID records the node number that has another node cross over or under it. 04076 // crosses records the crossing identified by node number that crosses over the crossesID. 04077 // E.g. real cross# 1 2 3 4 5 6 04078 // crossesID: 30 31 32 74 77 79 04079 // crosses: -74 77 -79 30 -31 32 04080 // results in the Dowker notation below 04081 // Dowker#: (1,4) (3,6) (5,2) 04082 // which is a half knot 04083 crossesID.push_back( i+1 ); 04084 crosses.push_back( m_iaCrossingsForDowkerNotation[i] ); 04085 } 04086 } 04087 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04088 04089 int numOfCrosses = static_cast<int>( crosses.size() ); 04090 if ( numOfCrosses == 0 ) { // Not a knot 04091 m_iKnotStartPt = m_iKnotEndPt = -1; 04092 m_iNumCrossingPairs = 0; 04093 } 04094 else if ( numOfCrosses % 2 == 0 ) { 04095 m_iNumCrossingPairs = numOfCrosses/2; 04096 assert( m_iNumCrossingPairs <= m_iMaxCrossingPairs ); 04097 int idx = 0; 04098 for ( int i = 0, odd = 1; i < numOfCrosses; i+=2, odd+=2 ) { 04099 for ( int j = 1, even = 2; j < numOfCrosses; j+=2, even+=2 ) { 04100 if ( crossesID[i] == abs( crosses[j] ) ) { 04101 //iaDowkerNotation[ idx++ ] = odd; 04102 if ( crosses[j] < 0 ) 04103 m_iaDowkerNotation[ idx++ ] = -even; 04104 else 04105 m_iaDowkerNotation[ idx++ ] = even; 04106 } 04107 } 04108 } 04109 04111 m_iKnotStartPt = crossesID[0]; 04112 m_iKnotEndPt = crossesID[numOfCrosses-1]; 04113 } 04114 else { // Not a knot 04115 m_iKnotStartPt = m_iKnotEndPt = -1; 04116 m_iNumCrossingPairs = numOfCrosses/2; 04117 } 04118 04119 #ifdef TAPs_STRAND_DEBUG 04120 debug_crosses = crosses; 04121 debug_crossesID = crossesID; 04122 std::cout << Debug_Str(); 04123 #endif//TAPs_STRAND_DEBUG 04124 04125 // Find the knot that match with the Dowker notation 04126 int knotID = DowkerNotationToKnotName(); 04127 if ( knotID >= 0 ) { 04128 //std::cout << "zDir found Knot ID# " << knotID << "\n"; 04129 return knotID; 04130 } 04131 04132 //------------------------------------------------- 04133 #ifdef TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 04134 // Check Clumped 04135 if ( m_iNumCrossingPairs > 1 ) { 04136 if ( m_ucTrackingKnot != Data::DatabaseDowkerNotation::UNDEFINED ) { 04137 T sphereRadius = m_tLinkLength * 10.0; 04138 CheckClumped( m_iKnotStartPt, m_iKnotEndPt, sphereRadius ); 04139 if ( IsClumped() ) { 04140 #ifdef TAPs_STRAND_DEBUG 04141 std::cout << "CLUMPED: true (" << m_iKnotStartPt << "," << m_iKnotEndPt << ")" << "\n"; 04142 #endif//TAPs_STRAND_DEBUG 04143 04144 m_uiCounterForFixingKnot = m_uiCounterThreshold+1; 04145 04146 #ifdef TAPs_STRAND_DEBUG 04147 std::cout << "RET: " << m_ucTrackingKnot << "\n"; 04148 #endif//TAPs_STRAND_DEBUG 04149 04150 return m_ucTrackingKnot; 04151 } 04152 } 04153 else { 04154 //std::cout << "CLUMPED: false (" << m_iKnotStartPt << "," << m_iKnotEndPt << ")" << std::endl; 04155 } 04156 } 04157 #endif//TAPs_STRAND_KNOT_RECOGNITION_ADD_ONs 04158 //------------------------------------------------- 04159 04160 // return -1 if not found, otherwise the knotID 04161 //std::cout << "Not found Knot ID# " << knotID << "\n"; 04162 return knotID; 04163 } 04164 //------------------------------------------------------------------- 04165 04166 //================================================================= 04167 // START: Fns for CheckSelfCrossingForKnotRecognition function 04168 //----------------------------------------------------------------- 04169 // CheckSelfCrossingsNextLevelRecursively 04170 // Assume no intersections between immediate adjacent links 04171 // Assume a BVH node either has both children or none, therefore 04172 // checking either it has a left or right child is enough to determine 04173 // whether it is a leaf node or not. 04174 template <typename T> 04175 void ModelStrand<T>::CheckSelfCrossingsNextLevelRecursively ( 04176 BVHNode<T> * node1, BVHNode<T> * node2, 04177 Vector3<T> const & projectionDirection 04178 ) 04179 { 04180 CheckSelfCrossingsRecursively( node1, node2, projectionDirection ); 04181 if ( node1->Child(0) ) 04182 CheckSelfCrossingsNextLevelRecursively( node1->Child(0), node1->Child(1), projectionDirection ); 04183 if ( node2->Child(0) ) 04184 CheckSelfCrossingsNextLevelRecursively( node2->Child(0), node2->Child(1), projectionDirection ); 04185 } 04186 //----------------------------------------------------------------- 04187 // CheckSelfCrossingsRecursively 04188 // Assume no intersections between immediate adjacent links 04189 // Use the fact that the ids of sphere bvh leaf nodes are the same as 04190 // the link id and 04191 template <typename T> 04192 T ModelStrand<T>::CheckSelfCrossingsRecursively ( 04193 BVHNode<T> * node1, BVHNode<T> * node2, 04194 Vector3<T> const & projectionDirection 04195 ) 04196 { 04197 static T LinkDiameter = static_cast<T>(2.0) * m_tRadius; 04198 static T epsilon = m_tRadius * static_cast<T>(0.05); // for collision response 04199 //------------------------------------------------------------- 04200 if ( !node1->Child(0) && !node2->Child(0) && 04201 ( abs(node1->GetID() - node2->GetID()) == 1 ) ) 04202 { 04203 return DBL_MAX; 04204 } 04205 //------------------------------------------------------------- 04206 T overlap = TestOverlapCircle( node1, node2, projectionDirection ); 04207 //------------------------------------------------------------- 04208 if ( overlap < DBL_MIN ) { 04209 if ( node2->Child(0) == NULL ) { //&& node2->RightChild() == NULL ) { 04210 if ( node1->Child(0) ) { 04211 overlap = CheckSelfCrossingsRecursively( node1->Child(0), node2, projectionDirection ); 04212 } 04213 if ( node1->Child(1) ) { 04214 overlap = CheckSelfCrossingsRecursively( node1->Child(1), node2, projectionDirection ); 04215 } 04216 } 04217 else { 04218 if ( node2->Child(0) ) { 04219 overlap = TestOverlapCircle( node1, node2->Child(0), projectionDirection ); 04220 if ( overlap < DBL_MIN ) { 04221 if ( node1->Child(0) ) { 04222 overlap = CheckSelfCrossingsRecursively( node1->Child(0), node2->Child(0), projectionDirection ); 04223 } 04224 if ( node1->Child(1) ) { 04225 overlap = CheckSelfCrossingsRecursively( node1->Child(1), node2->Child(0), projectionDirection ); 04226 } 04227 } 04228 } 04229 if ( node2->Child(1) ) { 04230 overlap = TestOverlapCircle( node1, node2->Child(1), projectionDirection ); 04231 if ( overlap < DBL_MIN ) { 04232 if ( node1->Child(0) ) { 04233 overlap = CheckSelfCrossingsRecursively( node1->Child(0), node2->Child(1), projectionDirection ); 04234 } 04235 if ( node1->Child(1) ) { 04236 overlap = CheckSelfCrossingsRecursively( node1->Child(1), node2->Child(1), projectionDirection ); 04237 } 04238 } 04239 } 04240 } 04241 //--------------------------------------------------------- 04242 if ( node1->Child(0) == NULL && //node1->Child(1) == NULL && 04243 node2->Child(0) == NULL ) //&& node2->Child(1) == NULL ) 04244 { 04245 CheckCrossings( node1, node2, projectionDirection ); 04246 } 04247 } 04248 //------------------------------------------------------------- 04249 return overlap; 04250 } 04251 //----------------------------------------------------------------- 04252 template <typename T> 04253 T ModelStrand<T>::TestOverlapCircle ( 04254 BVHNode<T> * node1, BVHNode<T> * node2, 04255 Vector3<T> const & projectionDirection 04256 ) 04257 { 04258 Vector3<T> C1 = node1->GetCenter(); 04259 Vector3<T> C2 = node2->GetCenter(); 04260 04261 04262 C1.SetZ( 0 ); 04263 C2.SetZ( 0 ); 04264 return (C1 - C2).Length() - ( node1->GetRadius() + node2->GetRadius() ); 04265 } 04266 //----------------------------------------------------------------- 04267 template <typename T> 04268 void ModelStrand<T>::CheckCrossings ( 04269 BVHNode<T> * node1, BVHNode<T> * node2, 04270 Vector3<T> const & projectionDirection 04271 ) 04272 { 04273 int i = node1->GetID(); 04274 int j = node2->GetID(); 04275 04276 // Skip if a node is locked 04277 #ifdef TAPs_STRAND_LOCK_SEGMENTS 04278 if ( m_prIsLocked[i] ) return; 04279 if ( m_prIsLocked[j] ) return; 04280 #endif//TAPs_STRAND_LOCK_SEGMENTS 04281 04282 Matrix3x3<T> rotationMatrix = CGMath<T>::CreateRotationMatrix3x3FromVectorAtoVectorB( 04283 Vector3<T>( 0, 0, 1 ), projectionDirection ); 04284 04285 // Link# i 04286 Vector3<T> P1 = rotationMatrix * m_prVertex[i]; 04287 Vector3<T> P2 = rotationMatrix * m_prVertex[i+1]; 04288 // Link# j 04289 Vector3<T> Q1 = rotationMatrix * m_prVertex[j]; 04290 Vector3<T> Q2 = rotationMatrix * m_prVertex[j+1]; 04291 04292 // Project in z-direction 04293 Vector3<T> p1( P1[0], P1[1], 0 ); 04294 Vector3<T> p2( P2[0], P2[1], 0 ); 04295 Vector3<T> q1( Q1[0], Q1[1], 0 ); 04296 Vector3<T> q2( Q2[0], Q2[1], 0 ); 04297 04298 bool isTouched; 04299 T pt, qt; // parameter for the 1st & 2nd line 04300 // where intersection point = p1 + pt(p2-p1) 04301 // where intersection point = q1 + qt(q2-q1) 04302 CGMath<T>::FindIntersectionLineSegmentLineSegment( 04303 p1, p2, q1, q2, // I/P 04304 pt, qt, // O/P 04305 isTouched // O/P 04306 ); 04307 if ( 0 <= pt && pt <= 1 && 0 <= qt && qt <= 1 ) { 04308 Vector3<T> P = P1 + pt*( P2 - P1 ); 04309 Vector3<T> Q = Q1 + qt*( Q2 - Q1 ); 04310 04311 // i+1 and j+1 are for making the link# to be one-based instead of zero-based 04312 if ( P[2] > Q[2] ) { 04313 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04314 04315 #ifdef DEBUG_TOO 04316 int a = i*10 + m_iaCrossingCounters[i]++; 04317 int b = j*10 + m_iaCrossingCounters[j]++; 04318 InsertToList_TOO( -b, a ); 04319 InsertToList_TOO( a, b ); 04320 #endif//DEBUG_TOO 04321 04322 InsertToList( -(j), i ); 04323 InsertToList( +(i), j ); 04324 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04325 m_iaCrossingsForDowkerNotation[i] = -(j+1); // the link i is above the link j 04326 m_iaCrossingsForDowkerNotation[j] = +(i+1); // the link j is below the link i 04327 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04328 } 04329 else { 04330 #ifdef TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04331 04332 #ifdef DEBUG_TOO 04333 int a = i*10 + m_iaCrossingCounters[i]++; 04334 int b = j*10 + m_iaCrossingCounters[j]++; 04335 InsertToList_TOO( b, a ); 04336 InsertToList_TOO( -a, b ); 04337 #endif//DEBUG_TOO 04338 04339 InsertToList( +(j), i ); 04340 InsertToList( -(i), j ); 04341 #else //TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04342 m_iaCrossingsForDowkerNotation[i] = +(j+1); // the link j is above the link i 04343 m_iaCrossingsForDowkerNotation[j] = -(i+1); // the link i is below the link j 04344 #endif//TAPs_PROJECTION_ALLOWS_MULTIPLE_CONTACT_POINTS 04345 } 04346 } 04347 } 04348 //----------------------------------------------------------------- 04349 // END: Fns for CheckSelfCrossingForKnotRecognition function 04350 //================================================================= 04351 04352 //------------------------------------------------------------------- 04353 #endif//TAPs_ALLOW_SETTING_PROJECTION_DIRECTION 04354 //=================================================================== 04355 04356 04357 //------------------------------------------------------------------- 04358 #endif//TAPs_STRAND_KNOT_CONTROL_USE_DOWKER_NOTATION_VERSION_01 04359 //=================================================================== 04360 04361 04362 //----------------------------------------------------------------------------- 04363 #endif//TAPs_STRAND_KNOT_CONTROL 04364 //============================================================================= 04365 04366 04367 //============================================================================= 04368 #ifdef TAPs_USE_CUDA 04369 //----------------------------------------------------------------------------- 04370 template <typename T> 04371 bool ModelStrand<T>::CUDA_Initialize_All () 04372 { 04373 bool result = true; 04374 if ( !CUDA_Initialize_VertexList() ) result = false; 04375 if ( !CUDA_Initialize_PrevVertexList() ) result = false; 04376 04377 #ifdef TAPs_ADVANCED_SIMULATION 04378 if ( !CUDA_Initialize_SimFlagsList() ) result = false; 04379 if ( !CUDA_Initialize_PosConstraintList() ) result = false; 04380 #endif//TAPs_ADVANCED_SIMULATION 04381 04382 //std::cout << "m_cudaID: " << m_cudaID << "\n"; 04383 04384 #ifdef TAPs_ADVANCED_SIMULATION 04385 TAPs::CUDA::InitailizeDataForSutureModel_ADVSIM( 04386 m_cudaID, 04387 m_iNoLinks+1, 04388 m_cudaVertexList, 04389 m_cudaPrevVertexList, 04390 m_cudaSimFlagsList, 04391 m_cudaPosConstraintList 04392 ); 04393 #else //TAPs_ADVANCED_SIMULATION 04394 TAPs::CUDA::InitailizeDataForSutureModel( 04395 m_cudaID, 04396 m_iNoLinks+1, 04397 m_cudaVertexList, 04398 m_cudaPrevVertexList 04399 ); 04400 #endif//TAPs_ADVANCED_SIMULATION 04401 04402 //std::cout << "m_cudaID: " << m_cudaID << "\n"; 04403 04404 return result; 04405 } 04406 04407 template <typename T> 04408 bool ModelStrand<T>::CUDA_Initialize_VertexList () 04409 { 04410 CUDA_Cleanup_VertexList(); 04411 04412 std::cout << "m_cudaVertexList size: " << (m_iNoLinks+1) * 4 << "\n"; 04413 04414 unsigned int size = sizeof(float) * (m_iNoLinks+1) * 4; 04415 04416 // Vertex List 04417 m_cudaVertexList = (float *)malloc( size ); 04418 if ( m_cudaVertexList ) { 04419 for ( unsigned int i = 0, p = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04420 m_cudaVertexList[p++] = m_prVertex[i][0]; 04421 m_cudaVertexList[p++] = m_prVertex[i][1]; 04422 m_cudaVertexList[p++] = m_prVertex[i][2]; 04423 m_cudaVertexList[p++] = 1; 04424 } 04425 } 04426 else { 04427 std::cout << "ERROR: Could not allocate memory for (Suture Model) CUDA!" << std::endl; 04428 exit( -1 ); 04429 } 04430 04431 return true; 04432 } 04433 04434 template <typename T> 04435 bool ModelStrand<T>::CUDA_Initialize_PrevVertexList () 04436 { 04437 CUDA_Cleanup_PrevVertexList(); 04438 04439 std::cout << "m_cudaPrevVertexList size: " << (m_iNoLinks+1) * 4 << "\n"; 04440 04441 unsigned int size = sizeof(float) * (m_iNoLinks+1) * 4; 04442 04443 // Previous Vertex List 04444 m_cudaPrevVertexList = (float *)malloc( size ); 04445 if ( m_cudaPrevVertexList ) { 04446 CUDA_CopyPrevVertexToMem(); 04447 } 04448 else { 04449 std::cout << "ERROR: Could not allocate memory for (Previous Vertex List) CUDA!" << std::endl; 04450 exit( -1 ); 04451 } 04452 04453 return true; 04454 } 04455 04456 template <typename T> 04457 void ModelStrand<T>::CUDA_Cleanup_All () 04458 { 04459 CUDA_Cleanup_VertexList(); 04460 CUDA_Cleanup_PrevVertexList(); 04461 04462 #ifdef TAPs_ADVANCED_SIMULATION 04463 CUDA_Cleanup_SimFlagsList(); 04464 CUDA_Cleanup_PosConstraintList(); 04465 #endif//TAPs_ADVANCED_SIMULATION 04466 } 04467 04468 template <typename T> 04469 void ModelStrand<T>::CUDA_Cleanup_VertexList () 04470 { 04471 if ( m_cudaVertexList ) { 04472 free( m_cudaVertexList ); 04473 m_cudaVertexList = NULL; 04474 } 04475 } 04476 04477 template <typename T> 04478 void ModelStrand<T>::CUDA_Cleanup_PrevVertexList () 04479 { 04480 if ( m_cudaPrevVertexList ) { 04481 free( m_cudaPrevVertexList ); 04482 m_cudaPrevVertexList = NULL; 04483 } 04484 } 04485 04486 template <typename T> 04487 void ModelStrand<T>::CUDA_CopyVertexToMem () 04488 { 04489 if ( !m_cudaVertexList ) return; 04490 for ( unsigned int i = 0, p = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04491 m_cudaVertexList[p++] = m_prVertex[i][0]; 04492 m_cudaVertexList[p++] = m_prVertex[i][1]; 04493 m_cudaVertexList[p++] = m_prVertex[i][2]; 04494 //m_cudaVertexList[p++] = 1; 04495 if ( GetFixStatusOfPtNo( i ) ) m_cudaVertexList[p++] = 0; // the point is fixed 04496 else m_cudaVertexList[p++] = 1; // the point is not fixed 04497 } 04498 } 04499 04500 template <typename T> 04501 void ModelStrand<T>::CUDA_CopyPrevVertexToMem () 04502 { 04503 if ( !m_cudaPrevVertexList ) return; 04504 for ( unsigned int i = 0, p = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04505 //* 04506 m_cudaPrevVertexList[p++] = m_prVertexPrevPos[i][0]; 04507 m_cudaPrevVertexList[p++] = m_prVertexPrevPos[i][1]; 04508 m_cudaPrevVertexList[p++] = m_prVertexPrevPos[i][2]; 04509 m_cudaPrevVertexList[p++] = 1; 04510 //*/ 04511 04512 /* 04513 m_cudaPrevVertexList[p++] = m_prVertex[i][0]; 04514 m_cudaPrevVertexList[p++] = m_prVertex[i][1]; 04515 m_cudaPrevVertexList[p++] = m_prVertex[i][2]; 04516 m_cudaPrevVertexList[p++] = 1; 04517 //*/ 04518 } 04519 } 04520 04521 template <typename T> 04522 void ModelStrand<T>::CUDA_CopyMemToVertex () 04523 { 04524 if ( !m_cudaVertexList ) return; 04525 for ( unsigned int i = 0, p = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04526 m_prVertex[i].SetXYZ( 04527 m_cudaVertexList[p ], 04528 m_cudaVertexList[p+1], 04529 m_cudaVertexList[p+2] 04530 ); 04531 p += 4; 04532 } 04533 } 04534 04535 template <typename T> 04536 void ModelStrand<T>::CUDA_CopyMemToPrevVertex () 04537 { 04538 if ( !m_cudaPrevVertexList ) return; 04539 for ( unsigned int i = 0, p = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04540 m_prVertexPrevPos[i].SetXYZ( 04541 m_cudaPrevVertexList[p ], 04542 m_cudaPrevVertexList[p+1], 04543 m_cudaPrevVertexList[p+2] 04544 ); 04545 p += 4; 04546 } 04547 } 04548 04549 //----------------------------------------------------------------------------- 04550 #ifdef TAPs_ADVANCED_SIMULATION 04551 //----------------------------------------------------------------------------- 04552 template <typename T> 04553 bool ModelStrand<T>::CUDA_Initialize_SimFlagsList () 04554 { 04555 CUDA_Cleanup_SimFlagsList(); 04556 04557 std::cout << "m_cudaSimFlagsList size: " << (m_iNoLinks+1) << "\n"; 04558 04559 unsigned int size = sizeof(unsigned int) * (m_iNoLinks+1); 04560 04561 // Sim Flags List 04562 m_cudaSimFlagsList = (unsigned int *)malloc( size ); 04563 if ( m_cudaSimFlagsList ) { 04564 CUDA_CopySimFlagsToMem(); 04565 } 04566 else { 04567 std::cout << "ERROR: Could not allocate memory for (Suture Model) CUDA!" << std::endl; 04568 exit( -1 ); 04569 } 04570 04571 return true; 04572 } 04573 template <typename T> 04574 bool ModelStrand<T>::CUDA_Initialize_PosConstraintList () 04575 { 04576 CUDA_Cleanup_PosConstraintList(); 04577 04578 std::cout << "m_cudaPosConstraintList size: " << (m_iNoLinks+1) * 4 << "\n"; 04579 04580 unsigned int size = sizeof(float) * (m_iNoLinks+1) * 4; 04581 04582 // Position Constraint List 04583 m_cudaPosConstraintList = (float *)malloc( size ); 04584 if ( m_cudaPosConstraintList ) { 04585 CUDA_CopyPosConstraintToMem(); 04586 } 04587 else { 04588 std::cout << "ERROR: Could not allocate memory for (Suture Model) CUDA!" << std::endl; 04589 exit( -1 ); 04590 } 04591 04592 return true; 04593 } 04594 template <typename T> 04595 void ModelStrand<T>::CUDA_Cleanup_SimFlagsList () 04596 { 04597 if ( m_cudaSimFlagsList ) { 04598 free( m_cudaSimFlagsList ); 04599 m_cudaSimFlagsList = NULL; 04600 } 04601 } 04602 template <typename T> 04603 void ModelStrand<T>::CUDA_Cleanup_PosConstraintList () 04604 { 04605 if ( m_cudaPosConstraintList ) { 04606 free( m_cudaPosConstraintList ); 04607 m_cudaPosConstraintList = NULL; 04608 } 04609 } 04610 template <typename T> 04611 void ModelStrand<T>::CUDA_CopySimFlagsToMem () 04612 { 04613 if ( !m_cudaSimFlagsList ) return; 04614 04615 for ( unsigned int i = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04616 m_cudaSimFlagsList[i] = m_SimFlags[i].GetSimulationConstraints().FlagsToUnsignedInt(); 04617 } 04618 } 04619 template <typename T> 04620 void ModelStrand<T>::CUDA_CopyPosConstraintToMem () 04621 { 04622 if ( !m_cudaPosConstraintList ) return; 04623 04624 AdvSimSupport_DS<T> * pSimData = AdvSimSupport<T>::GetAdvSimData(); 04625 std::vector< AdvSimSupport_DS_InteractionSutureAndModel<T> > * pDataInteractionSutureAndModel = &( pSimData->GetInteractionSutureAndModel( AdvSimID ) ); 04626 for ( unsigned int k = 0; k < static_cast<unsigned int>( pDataInteractionSutureAndModel->size() ); ++k ) { 04627 unsigned int idx = 4 * (*pDataInteractionSutureAndModel)[k].VertexIDSuture; 04628 // Position constraint 04629 m_cudaPosConstraintList[idx++] = (*pDataInteractionSutureAndModel)[k].VertexModel->GetPosition()[0]; 04630 m_cudaPosConstraintList[idx++] = (*pDataInteractionSutureAndModel)[k].VertexModel->GetPosition()[1]; 04631 m_cudaPosConstraintList[idx++] = (*pDataInteractionSutureAndModel)[k].VertexModel->GetPosition()[2]; 04632 // Force ratio 04633 m_cudaPosConstraintList[idx++] = (*pDataInteractionSutureAndModel)[k].ForceRatio; 04634 } 04635 } 04636 template <typename T> 04637 void ModelStrand<T>::CUDA_CopyMemToSimFlags () 04638 { 04639 if ( !m_cudaSimFlagsList ) return; 04640 04641 for ( unsigned int i = 0; i <= static_cast<unsigned int>( m_iNoLinks ); ++i ) { 04642 m_SimFlags[i].GetSimulationConstraints() = m_cudaSimFlagsList[i]; 04643 } 04644 } 04645 template <typename T> 04646 void ModelStrand<T>::CUDA_CopyMemToPosConstraint () 04647 { 04648 if ( !m_cudaPosConstraintList ) return; 04649 04650 AdvSimSupport_DS<T> * pSimData = AdvSimSupport<T>::GetAdvSimData(); 04651 std::vector< AdvSimSupport_DS_InteractionSutureAndModel<T> > * pDataInteractionSutureAndModel = &( pSimData->GetInteractionSutureAndModel( AdvSimID ) ); 04652 for ( unsigned int k = 0; k < static_cast<unsigned int>( pDataInteractionSutureAndModel->size() ); ++k ) { 04653 int vertexID = (*pDataInteractionSutureAndModel)[k].VertexIDSuture; 04654 unsigned int idx = 4 * vertexID; 04655 T forceRatio = m_cudaPosConstraintList[idx+3]; 04656 04657 std::cout << "Suture Vertex Sim Flags# " << vertexID << ":\t" << m_SimFlags[vertexID] << " " << m_SimFlags[vertexID].GetSimulationConstraints() << "\n"; 04658 std::cout << "Constraint VertexID: " << vertexID << "" << "\tForceRatio: " << forceRatio << "\n"; 04659 04660 // SHIFTING OF SUTURE DUE TO SLIPDABLE CONSTRAINTS 04661 // Since AdvSimSupport_DS_InteractionSutureAndModel class contains 04662 //unsigned int VertexIDSuture; //!< vertex ID of suture 04663 //DS::SimulationFlags * SimFlagsSuture; //!< simulation flags of suture Vector3<T> * VertexSuture; //!< vertex of suture 04664 //Vector3<T> * VertexSuture; //!< a vertex of suture 04665 //HEVertex<T> * VertexModel; //!< vertex of model 04666 //T ForceRatio; //!< force ratio [0-1] (1 means suture's vertex completely dominates model's) 04667 // So VertexIDSuture, 7SimFlagsSuture, and &VertexSuture have to be updated! 04668 04669 // The shift values are encoded in forceRatio. 04670 04671 // Shift up 04672 if ( forceRatio > TAPs_Const_Check_Suture_ShiftUp ) { 04673 // Shift up 04674 if ( vertexID < GetNumberOfLinks() ) { 04675 //m_SimFlags[vertexID+1] = m_SimFlags[vertexID]; 04676 04677 /* 04678 std::cout << "m_SimFlags[vertexID+1]: " << m_SimFlags[vertexID+1] << "\n"; 04679 if ( m_SimFlags[vertexID+1].CheckSimulationConstraints( TAPs::Enum::AddOn::CLEARED ) ) { 04680 std::cout << "CLEARED\n"; 04681 } 04682 else { 04683 std::cout << "NOT CLEARED\n"; 04684 } 04685 //*/ 04686 04687 04688 if ( m_SimFlags[vertexID+1].AreSimulationFlagsCleared() ) { 04689 m_SimFlags[vertexID+1].SetSimulationConstraints( TAPs::Enum::AddOn::SLIDABLE ); 04690 m_SimFlags[vertexID].ClearSimulationConstraints( TAPs::Enum::AddOn::SLIDABLE ); 04691 ++((*pDataInteractionSutureAndModel)[k].VertexIDSuture); 04692 04693 //std::cout << "Shift up: " << "VertexID: " << vertexID 04694 // << " Prev #id: " << (*pDataInteractionSutureAndModel)[k].VertexIDSuture - 1 04695 // << " New #id: " << (*pDataInteractionSutureAndModel)[k].VertexIDSuture 04696 // << " sim flags of # " << vertexID << ": " << m_SimFlags[vertexID] 04697 // << " sim flags of # " << vertexID+1 << ": " << m_SimFlags[vertexID+1] << "\n"; 04698 04699 // Position constraint 04700 (*pDataInteractionSutureAndModel)[k].VertexModel->SetPosition( m_cudaPosConstraintList[idx], m_cudaPosConstraintList[idx+1], m_cudaPosConstraintList[idx+2] ); 04701 // Force ratio 04702 (*pDataInteractionSutureAndModel)[k].ForceRatio = forceRatio - TAPs_Const_Suture_ShiftUp; 04703 // Ptr to Suture Vertex 04704 (*pDataInteractionSutureAndModel)[k].VertexSuture = &m_prVertex[vertexID+1]; 04705 // Ptr to Suture SimFlag 04706 (*pDataInteractionSutureAndModel)[k].SimFlagsSuture = &m_SimFlags[vertexID+1]; 04707 } 04708 else { 04709 // Force ratio 04710 (*pDataInteractionSutureAndModel)[k].ForceRatio = forceRatio - TAPs_Const_Suture_ShiftUp; 04711 } 04712 } 04713 } 04714 // Shift down 04715 else if ( forceRatio > TAPs_Const_Check_Suture_ShiftDown ) { 04716 // Shift down 04717 if ( vertexID > 0 ) { 04718 //m_SimFlags[vertexID-1] = m_SimFlags[vertexID]; 04719 04720 if ( m_SimFlags[vertexID-1].AreSimulationFlagsCleared() ) { 04721 m_SimFlags[vertexID-1].SetSimulationConstraints( TAPs::Enum::AddOn::SLIDABLE ); 04722 m_SimFlags[vertexID].ClearSimulationConstraints( TAPs::Enum::AddOn::SLIDABLE ); 04723 --((*pDataInteractionSutureAndModel)[k].VertexIDSuture); 04724 04725 //std::cout << "Shift down: " << "VertexID: " << vertexID 04726 // << " Prev #id: " << (*pDataInteractionSutureAndModel)[k].VertexIDSuture + 1 04727 // << " New #id: " << (*pDataInteractionSutureAndModel)[k].VertexIDSuture 04728 // << " sim flags of # " << vertexID << ": " << m_SimFlags[vertexID] 04729 // << " sim flags of # " << vertexID-1 << ": " << m_SimFlags[vertexID-1] << "\n"; 04730 04731 // Position constraint 04732 (*pDataInteractionSutureAndModel)[k].VertexModel->SetPosition( m_cudaPosConstraintList[idx], m_cudaPosConstraintList[idx+1], m_cudaPosConstraintList[idx+2] ); 04733 // Force ratio 04734 (*pDataInteractionSutureAndModel)[k].ForceRatio = forceRatio - TAPs_Const_Suture_ShiftDown; 04735 // Ptr to Suture Vertex 04736 (*pDataInteractionSutureAndModel)[k].VertexSuture = &m_prVertex[vertexID-1]; 04737 // Ptr to Suture SimFlag 04738 (*pDataInteractionSutureAndModel)[k].SimFlagsSuture = &m_SimFlags[vertexID-1]; 04739 } 04740 else { 04741 // Force ratio 04742 (*pDataInteractionSutureAndModel)[k].ForceRatio = forceRatio - TAPs_Const_Suture_ShiftDown; 04743 } 04744 } 04745 } 04746 else // NO SHIFT 04747 { 04748 // Position constraint 04749 (*pDataInteractionSutureAndModel)[k].VertexModel->SetPosition( m_cudaPosConstraintList[idx], m_cudaPosConstraintList[idx+1], m_cudaPosConstraintList[idx+2] ); 04750 // Force ratio 04751 (*pDataInteractionSutureAndModel)[k].ForceRatio = forceRatio; 04752 } 04753 } 04754 } 04755 //----------------------------------------------------------------------------- 04756 #endif//TAPs_ADVANCED_SIMULATION 04757 //----------------------------------------------------------------------------- 04758 04759 template <typename T> 04760 void ModelStrand<T>::CUDA_SwapBuffers () 04761 { 04762 float * tmpPtr = m_cudaVertexList; 04763 m_cudaVertexList = m_cudaPrevVertexList; 04764 m_cudaPrevVertexList = tmpPtr; 04765 04766 Vector3<T> * tmpPtrV = m_prVertex; 04767 m_prVertex = m_prVertexPrevPos; 04768 m_prVertexPrevPos = tmpPtrV; 04769 } 04770 04771 //----------------------------------------------------------------------------- 04772 #endif//TAPs_USE_CUDA 04773 //============================================================================= 04774 04775 04776 //----------------------------------------------------------------------------- 04777 //============================================================================= 04778 END_NAMESPACE_TAPs__OpenGL 04779 //----------------------------------------------------------------------------- 04780 //345678901234567890123456789012345678901234567890123456789012345678901234567890 04781 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1.7.4 
