 |
TAPs 0.7.7.3
|
|
TAPs 0.7.7.3
|
00001 /****************************************************************************** 00002 TAPsVolPresPolygonalModel.cpp 00003 00004 A Polygonal Model that is deformable with volume preservation. 00005 00006 SUKITTI PUNAK (11/11/2004) 00007 UPDATE (11/15/2004) 00008 ******************************************************************************** 00009 +====================+ 00010 | DESCRIPTION / NOTE | 00011 +====================+ 00012 00013 The model is a polygonal mesh of V (number of vertices) masses, which each 00014 mass linked to its neighbors by massless springs of natural length greater than 00015 zeroes. 00016 00017 00018 CHANGE THE COMMENT TO FROM TRIGONAL TO POLYGONAL 00019 Since the model is triangulated mesh, the linkage between neighbors comprised 00020 of two kinds of linking springs: 00021 -- Structural springs linking vertices of a face f along its edges 00022 (i.e. metric deformation) 00023 e.g. 2-----1 2 00024 \ / /| 00025 \ / / | 00026 0 0--1 00027 -- Flexion springs linking diagonally between two faces 00028 (i.e. bending deformation) 00029 e.g. 2-----1 this example shows the link of vertex 2 and 3 00030 \ + / \ 00031 \ / + \ 00032 0-----3 00033 (-- No Shear springs (i.e. shearing deformation)) 00034 00035 (2)-----------(2)-----------(2) 00036 / \ /|\ / \ 00037 / \ / | \ / \ 00038 / \ / | \ / \ 00039 / \ / | \ / \ 00040 / \ / | \ / \ 00041 / \ / | \ / \ 00042 (2)-----------(1)-----|-----(1)-----------(2) 00043 / \ . / \ | / \ . / \ 00044 / \ . / \ | / \ . / \ 00045 / \ . / \ | / \ . / \ 00046 / \ / . \ | / . \ / \ 00047 / \ / . \ | / . \ / \ 00048 / \ / . \|/ . \ / \ 00049 (2)-----------(1)-----------(0)-----------(1)----------(2) 00050 \ / \ . /|\ . / \ / 00051 \ / \ . / | \ . / \ / 00052 \ / \ . / | \ . / \ / 00053 \ / . \ / | \ / . \ / 00054 \ / . \ / | \ / . \ / 00055 \ / . \ / | \ / . \ / 00056 (2)-----------(1)-----|-----(1)-----------(2) 00057 \ / \ | / \ / 00058 \ / \ | / \ / 00059 \ / \ | / \ / 00060 \ / \ | / \ / 00061 \ / \ | / \ / 00062 \ / \|/ \ / 00063 (2)------------(2)----------(2) 00064 00065 Data Structure: 00066 =============== 00067 Need edge list 00068 Each edge incidents to two faces, find cross list 00069 00070 00071 00072 In case of Quadrateral Mesh 00073 (0,0)---------(0,1)---------(0,2) 00074 | \ / | \ / | 00075 | \ / | \ / | 00076 | + | + | 00077 | / \ | / \ | 00078 | / \ | / \ | 00079 (1,0)---------(1,1)---------(1,2) 00080 | \ / | \ / | 00081 | \ / | \ / | 00082 | + | + | 00083 | / \ | / \ | 00084 | / \ | / \ | 00085 (2,0)---------(2,1)---------(2,2) 00086 00087 ******************************************************************************/ 00088 #include "TAPsVolPresPolygonalModel.hpp" 00089 // Using Inclusion Model (i.e. definitions are included in declarations) 00090 // (this name.cpp is included in name.hpp) 00091 // Each friend is defined directly inside its declaration. 00092 00093 BEGIN_NAMESPACE_TAPs__OpenGL 00094 //============================================================================= 00095 //----------------------------------------------------------------------------- 00096 // default constructor 00097 template <typename T> 00098 VolPresPolygonalModel<T>::VolPresPolygonalModel () 00099 : XPolygonalModel<T>(), m_tVolume(), m_vtCentroid(), 00100 m_iNoEdges( 0 ), m_prEdge( NULL ), 00101 // Spring 00102 m_pprParticleVertex( NULL ), 00103 m_pprSpringVertex( NULL ), 00104 m_bFlagFixed( NULL ), 00105 m_pvDeformVector( NULL ), 00106 m_iDeformedFaceNo( 0 ), 00107 m_pcVertexRings( NULL ), 00108 m_bFlagVolPresMode( true ), 00109 m_bFlagSpringForceMode( true ) 00110 00111 { 00112 #ifdef TAPs_ENABLE_DEBUG 00113 std::cout << "VolPresPolygonalModel<" << typeid(T).name() << "> Constructor\n"; 00114 #endif//TAPs_ENABLE_DEBUG 00115 } 00116 //----------------------------------------------------------------------------- 00117 // destructor 00118 template <typename T> 00119 VolPresPolygonalModel<T>::~VolPresPolygonalModel () 00120 { 00121 //---------------------------------------------------------------- 00122 // Delete All Edges 00123 delete [] m_prEdge; 00124 m_prEdge = NULL; 00125 00126 #ifdef TAPs_ENABLE_DEBUG 00127 std::cout << "VolPresPolygonalModel<" << typeid(T).name() << "> Destructor\n"; 00128 #endif//TAPs_ENABLE_DEBUG 00129 //---------------------------------------------------------------- 00130 // Delete All ParticleVertex 00131 for ( int i = 0; i < m_iNoVertices; ++i ) { 00132 delete m_pprParticleVertex[i]; 00133 } 00134 delete [] m_pprParticleVertex; 00135 //---------------------------------------------------------------- 00136 // Delete All SpringVertex 00137 for ( int i = 0; i < m_iNoEdges; ++i ) { 00138 delete m_pprSpringVertex[i]; 00139 } 00140 delete [] m_pprSpringVertex; 00141 //---------------------------------------------------------------- 00142 // Delete All Fixed Status 00143 delete [] m_bFlagFixed; 00144 //---------------------------------------------------------------- 00145 // Delete All Deform Vectors 00146 delete [] m_pvDeformVector; 00147 //---------------------------------------------------------------- 00148 // Delete Vertex Rings Info 00149 delete m_pcVertexRings; 00150 } 00151 //----------------------------------------------------------------------------- 00152 // Initialize 00153 template <typename T> 00154 void VolPresPolygonalModel<T>::Initialize () 00155 { 00156 //---------------------------------------------------------------- 00157 // Model Initialization 00158 DetermineAndSortRings(); 00159 CalAndSetNormals(); 00160 // DetermineAndSortRings(); 00161 ApplyMaterial(); 00162 DetermineEdgeList(); 00163 CalVolume2(); 00164 00165 #ifdef TAPs_ENABLE_DEBUG 00166 std::cout << "Volume: " << m_tVolume << std::endl; 00167 #endif//TAPs_ENABLE_DEBUG 00168 00169 if ( m_tVolume < Math<T>::ZERO ) { 00170 //if ( false ) { 00171 #ifdef TAPs_ENABLE_DEBUG 00172 std::cout << "Recalculate Volume\n"; 00173 #endif//TAPs_ENABLE_DEBUG 00174 00175 int vSize, vNo, k; 00176 T s, t; 00177 for ( int i = 0; i < GetNoFaces(); ++i ) { 00178 vSize = GetFaceList()[i].GetNoVertices(); 00179 for ( int j = 0; j < vSize/2; ++j ) { 00180 k = vSize - 1 - j; 00181 vNo = GetFaceList()[i].GetVertexNo( j ); 00182 GetFaceList()[i].SetVertexNo( j, GetFaceList()[i].GetVertexNo(k) ); 00183 GetFaceList()[i].SetVertexNo( k, vNo ); 00184 00185 s = GetFaceList()[i].GetTexCoordHalfNo( j*2 ); 00186 t = GetFaceList()[i].GetTexCoordHalfNo( j*2 + 1 ); 00187 GetFaceList()[i].SetTexCoordNo( j, 00188 GetFaceList()[i].GetTexCoordHalfNo( k*2 ), 00189 GetFaceList()[i].GetTexCoordHalfNo( k*2 + 1 ) ); 00190 GetFaceList()[i].SetTexCoordNo( k, s, t); 00191 00192 } 00193 /* 00194 if ( vSize % 2 == 1 ) { 00195 j = vSize/2 + 1; 00196 s = GetFaceList()[i].GetTexCoordHalfNo( j*2 ); 00197 t = GetFaceList()[i].GetTexCoordHalfNo( j*2 + 1 ); 00198 GetFaceList()[i].SetTexCoordNo( j, s, t ); 00199 } 00200 */ 00201 } 00202 DetermineAndSortRings(); 00203 CalAndSetNormals(); 00204 // DetermineAndSortRings(); 00205 ApplyMaterial(); 00206 DetermineEdgeList(); 00207 CalVolume2(); 00208 00209 #ifdef TAPs_ENABLE_DEBUG 00210 std::cout << "Volume: " << m_tVolume << std::endl; 00211 #endif//TAPs_ENABLE_DEBUG 00212 } 00213 00214 m_pcVertexRings = new VertexRings( m_pviVertexRing1List, m_iNoVertices ); 00215 //PreserveVolume( weight ); 00216 00217 //---------------------------------------------------------------- 00218 // Allocate Fix Status and Deform Vector for each Vertex 00219 delete [] m_bFlagFixed; 00220 m_bFlagFixed = new bool[m_iNoVertices]; 00221 for ( int i = 0; i < m_iNoVertices; ++i ) { 00222 m_bFlagFixed[i] = false; 00223 } 00224 //---------------------------------------------------------------- 00225 // Allocate Deform Vector for each Vertex 00226 delete [] m_pvDeformVector; 00227 m_pvDeformVector = new Vector3<T>[m_iNoVertices]; 00228 00229 // For Collision Detection 00230 CalBoundingAABB(); 00231 CalBoundingEllipsoid(); 00232 CalBoundingSphere(); 00233 } 00234 //----------------------------------------------------------------------------- 00235 // Find Edge List 00236 template <typename T> 00237 void VolPresPolygonalModel<T>::DetermineEdgeList () 00238 { 00239 std::vector<int> *edgeRecord; 00240 edgeRecord = new std::vector<int>[ m_iNoVertices ]; 00241 int iOrig, iDest; // original vertex# and destination vertex# of edge 00242 int iV; // temp vertex number 00243 int numberOfEdges = 0; 00244 //------------------------------------------------------------------------- 00245 // Get edges in a temp data 00246 for ( int f = 0; f < m_iNoFaces; ++f ) { 00247 iV = m_prFace[f].GetVertexNo( 0 ); 00248 //---------------------------------------------------------------------- 00249 for ( int v = 1; v < m_prFace[f].GetNoVertices(); ++v ) { 00250 iOrig = iV; 00251 iDest = m_prFace[f].GetVertexNo( v ); 00252 //---------------------------------------------------------------------- 00253 // Make the original vertex always smaller than the destination vertex 00254 if ( iOrig > iDest ) { 00255 int tmp = iOrig; 00256 iOrig = iDest; 00257 iDest = tmp; 00258 } 00259 iV = iDest; // keep the current vertex 00260 //---------------------------------------------------------------------- 00261 // Record the edge 00262 bool notStored = true; 00263 for ( int e = 0; e < static_cast<int>( edgeRecord[iOrig].size() ); ++e ) { 00264 if ( edgeRecord[iOrig][e] == iDest ) { 00265 notStored = false; 00266 break; 00267 } 00268 } 00269 //---------------------------------------------------------------------- 00270 if ( notStored ) { 00271 edgeRecord[iOrig].push_back( iDest ); 00272 ++numberOfEdges; 00273 } 00274 } 00275 } 00276 //------------------------------------------------------------------------- 00277 // Set edges 00278 int edgeNo = 0; 00279 m_iNoEdges = numberOfEdges; 00280 m_prEdge = new Edge[numberOfEdges]; 00281 for ( int v = 0; v < m_iNoVertices; ++v ) { 00282 for ( int e = 0; e < static_cast<int>( edgeRecord[v].size() ); ++e ) { 00283 m_prEdge[edgeNo++].SetOrigAndDest( v, edgeRecord[v][e] ); 00284 } 00285 } 00286 //------------------------------------------------------------------------- 00287 CreateParticleVerticesFromVertices(); 00288 CreateSpringVerticesFromEdges(); 00289 } 00290 //----------------------------------------------------------------------------- 00291 // Create Particles From Vertices 00292 template <typename T> 00293 void VolPresPolygonalModel<T>::CreateParticleVerticesFromVertices () 00294 { 00295 m_pprParticleVertex = new ParticleVertex<T>*[m_iNoVertices]; 00296 for ( int i = 0; i < m_iNoVertices; ++i ) { 00297 m_pprParticleVertex[i] = new ParticleVertex<T> ( 00298 2.0, // Mass 00299 m_prXVertex[i], // Vertex (including Position and Normal Vectors) 00300 Vector3<T>(), // Velocity 00301 Vector3<T>(), // Acceleration 00302 Vector3<T>() // Force 00303 ); 00304 } 00305 } 00306 //----------------------------------------------------------------------------- 00307 // Create Springs From Edges 00308 template <typename T> 00309 void VolPresPolygonalModel<T>::CreateSpringVerticesFromEdges () 00310 { 00311 m_pprSpringVertex = new SpringVertex<T>*[m_iNoEdges]; 00312 for ( int i = 0; i < m_iNoEdges; ++i ) { 00313 m_pprSpringVertex[i] = new SpringVertex<T> ( 00314 *m_pprParticleVertex[m_prEdge[i].GetOriginal()], // Vertex#1 00315 *m_pprParticleVertex[m_prEdge[i].GetDestination()], // Vertex#2 00316 //0.05, // K (spring constant) 00317 0.025, // K (spring constant) 00318 0.25 // D (damping constant) 00319 ); 00320 } 00321 } 00322 //----------------------------------------------------------------------------- 00323 // Calculate Volume 00324 // For Triangle Face Only 00325 // Treat the triangle as a tetrahedron and calculate its volume 00326 template <typename T> 00327 T VolPresPolygonalModel<T>::CalVolume () 00328 { 00329 Vector3<T> va, vb, vc; 00330 Vector3<T> vd; 00331 Vector3<T> vmom; 00332 T tv; 00333 T tVolume = 0.0; 00334 for ( int i = 0; i < m_iNoFaces; ++i ) { 00335 va.SetXYZ ( 00336 m_prXVertex[m_prFace[i].GetVertexNo(0)][0], 00337 m_prXVertex[m_prFace[i].GetVertexNo(0)][1], 00338 m_prXVertex[m_prFace[i].GetVertexNo(0)][2] 00339 ); 00340 vb.SetXYZ ( 00341 m_prXVertex[m_prFace[i].GetVertexNo(1)][0], 00342 m_prXVertex[m_prFace[i].GetVertexNo(1)][1], 00343 m_prXVertex[m_prFace[i].GetVertexNo(1)][2] 00344 ); 00345 vc.SetXYZ ( 00346 m_prXVertex[m_prFace[i].GetVertexNo(2)][0], 00347 m_prXVertex[m_prFace[i].GetVertexNo(2)][1], 00348 m_prXVertex[m_prFace[i].GetVertexNo(2)][2] 00349 ); 00350 00351 tVolume += tv = ( va * ( vb ^ vc ) ) / 6.0; 00352 vd = ( va + vb + vc ) / 4.0; 00353 vmom += vd * tv; 00354 } 00355 //m_tVolume = tVolume; 00356 m_vtCentroid = vmom / tVolume; 00357 return tVolume; 00358 } 00359 //----------------------------------------------------------------------------- 00360 // Calculate Volume Too 00361 // If not a triangle face find the center (average) point, then virtually crete 00362 // triangle fan from the center point for volume calculations 00363 // Using d3mom0 ( v, dv, wgt ) subroutine to calculate a volume contribution by 00364 // the triangle(s) with all dv zeroes (zero moment) 00365 template <typename T> 00366 T VolPresPolygonalModel<T>::CalVolume2 () 00367 { 00368 T wgt[4] = { 0, 0, 0, 0 }; 00369 T v[3][3]; 00370 T dv[3][3]; 00371 00373 dv[0][0] = dv[0][1] = dv[0][2] = 0; 00374 dv[1][0] = dv[1][1] = dv[1][2] = 0; 00375 dv[2][0] = dv[2][1] = dv[2][2] = 0; 00376 //*/ 00377 //------------------------------------------------------------------------- 00378 for ( int i = 0; i < m_iNoFaces; ++i ) { 00379 //---------------------------------------------------------------------- 00380 // Tri Face 00381 if ( m_prFace[i].GetNoVertices() == 3 ) { 00382 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][0]; 00383 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][1]; 00384 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][2]; 00385 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][0]; 00386 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][1]; 00387 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][2]; 00388 v[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][0]; 00389 v[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][1]; 00390 v[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][2]; 00391 00392 d3mom0 ( v, dv, wgt ); // call subroutine 00393 } 00394 //---------------------------------------------------------------------- 00395 // Not Tri Face 00396 else { 00397 //---------------------------------------------- 00398 // Find the avg (center) vertex of the face 00399 XVertex<T> avgXVertex; // default ctor --> position and normal are zeros 00400 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00401 avgXVertex.SetPosition( avgXVertex.GetPosition() 00402 + m_prXVertex[m_prFace[i].GetVertexNo(j)].GetPosition() ); 00403 avgXVertex.SetNormal( avgXVertex.GetNormal() 00404 + m_prXVertex[m_prFace[i].GetVertexNo(j)].GetNormal() ); 00405 } 00406 //---------------------------------------------- 00407 int k = m_prFace[i].GetNoVertices() - 1; 00408 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00409 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][0]; 00410 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][1]; 00411 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][2]; 00412 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][0]; 00413 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][1]; 00414 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][2]; 00415 v[2][0] = avgXVertex[0]; 00416 v[2][1] = avgXVertex[1]; 00417 v[2][2] = avgXVertex[2]; 00418 00419 d3mom0 ( v, dv, wgt ); // call subroutine 00420 k = j; 00421 } 00422 } 00423 } 00424 00425 m_tVolume = wgt[0]; 00426 return m_tVolume; 00427 } 00428 //----------------------------------------------------------------------------- 00429 // Preserve Volume using Vertex Normals 00430 // If not a triangle face find the center (average) point, then virtually crete 00431 // triangle fan from the center point for volume calculations 00432 // Using d3mom0 ( v, dv, wgt ) subroutine to calculate a volume contribution by 00433 // the triangle(s), dv zeroes mean the point is fixed, else the point will have 00434 // to be moved in order to preserve the volume of the model. 00435 template <typename T> 00436 void VolPresPolygonalModel<T>::PreserveVolumeByVertexNormals ( T wgt[4] ) 00437 { 00438 wgt[0] = wgt[1] = wgt[2] = wgt[3] = 0; 00439 T v[3][3]; 00440 T dv[3][3]; 00441 00442 for ( int i = 0; i < m_iNoFaces; ++i ) { 00443 // Tri Face 00444 if ( m_prFace[i].GetNoVertices() == 3 ) { 00445 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][0]; 00446 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][1]; 00447 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][2]; 00448 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][0]; 00449 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][1]; 00450 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][2]; 00451 v[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][0]; 00452 v[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][1]; 00453 v[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][2]; 00454 //if ( true ) { 00455 if ( !m_bFlagFixed[ m_prFace[i].GetVertexNo(0) ] ) { 00456 dv[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetX(); 00457 dv[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetY(); 00458 dv[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetZ(); 00459 } 00460 else { 00461 dv[0][0] = dv[0][1] = dv[0][2] = Math<T>::ZERO; 00462 } 00463 //if ( true ) { 00464 if ( !m_bFlagFixed[ m_prFace[i].GetVertexNo(1) ] ) { 00465 dv[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetX(); 00466 dv[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetY(); 00467 dv[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetZ(); 00468 } 00469 else { 00470 dv[1][0] = dv[1][1] = dv[1][2] = Math<T>::ZERO; 00471 } 00472 //if ( true ) { 00473 if ( !m_bFlagFixed[ m_prFace[i].GetVertexNo(2) ] ) { 00474 dv[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetX(); 00475 dv[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetY(); 00476 dv[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetZ(); 00477 } 00478 else { 00479 dv[2][0] = dv[2][1] = dv[2][2] = Math<T>::ZERO; 00480 } 00481 d3mom0 ( v, dv, wgt ); // call subroutine 00482 } 00483 // Not Tri Face 00484 else { 00485 //---------------------------------------------- 00486 // Find the avg (center) vertex of the face 00487 XVertex<T> avgXVertex; // default ctor --> position and normal are zeros 00488 bool isOneFixed = false; // if at least a vertex of the face is fixed, 00489 // then fix all other vertices of the faces 00490 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00491 if ( !isOneFixed ) { 00492 isOneFixed = m_bFlagFixed[m_prFace[i].GetVertexNo(j)]; 00493 } 00494 avgXVertex.SetPosition( avgXVertex.GetPosition() 00495 + m_prXVertex[m_prFace[i].GetVertexNo(j)].GetPosition() ); 00496 avgXVertex.SetNormal( avgXVertex.GetNormal() 00497 + m_prXVertex[m_prFace[i].GetVertexNo(j)].GetNormal() ); 00498 } 00499 avgXVertex.SetPosition( avgXVertex.GetPosition() / static_cast<T>(m_prFace[i].GetNoVertices()) ); 00500 avgXVertex.GetNormal().Normalized(); 00501 if ( !isOneFixed ) { 00502 //---------------------------------------------- 00503 int k = m_prFace[i].GetNoVertices() - 1; 00504 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00505 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][0]; 00506 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][1]; 00507 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][2]; 00508 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][0]; 00509 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][1]; 00510 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][2]; 00511 v[2][0] = avgXVertex[0]; 00512 v[2][1] = avgXVertex[1]; 00513 v[2][2] = avgXVertex[2]; 00514 00515 dv[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ].GetNormal().GetX(); 00516 dv[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ].GetNormal().GetY(); 00517 dv[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ].GetNormal().GetZ(); 00518 dv[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ].GetNormal().GetX(); 00519 dv[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ].GetNormal().GetY(); 00520 dv[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ].GetNormal().GetZ(); 00521 dv[2][0] = avgXVertex.GetNormal().GetX(); 00522 dv[2][1] = avgXVertex.GetNormal().GetY(); 00523 dv[2][2] = avgXVertex.GetNormal().GetZ(); 00524 00525 d3mom0 ( v, dv, wgt ); // call subroutine 00526 k = j; 00527 } 00528 } 00529 else { 00530 //---------------------------------------------- 00531 int k = m_prFace[i].GetNoVertices() - 1; 00532 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00533 m_bFlagFixed[ m_prFace[i].GetVertexNo(j) ] = true; 00534 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][0]; 00535 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][1]; 00536 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][2]; 00537 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][0]; 00538 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][1]; 00539 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][2]; 00540 v[2][0] = avgXVertex[0]; 00541 v[2][1] = avgXVertex[1]; 00542 v[2][2] = avgXVertex[2]; 00543 00544 dv[0][0] = dv[0][1] = dv[0][2] = 00545 dv[1][0] = dv[1][1] = dv[1][2] = 00546 dv[2][0] = dv[2][1] = dv[2][2] = Math<T>::ZERO; 00547 00548 d3mom0 ( v, dv, wgt ); // call subroutine 00549 k = j; 00550 } 00551 } 00552 } 00553 } 00554 00555 T diffVolume = wgt[0] - m_tVolume; 00556 if ( diffVolume != 0 ) { 00557 //if ( fabs( diffVolume ) > Math<T>::ROUND_ERROR ) { 00558 T sigma = SP_Maths::rootOfCubicPolynomial( wgt[3], wgt[2], wgt[1], diffVolume, 1 ); 00559 //std::cout << "\nsigma: " << sigma << "\n"; 00560 for ( int i = 0; i < m_iNoVertices; ++i ) { 00561 if ( !m_bFlagFixed[i] ) { 00562 m_prXVertex[i].SetPosition( m_prXVertex[i].GetPosition() + sigma * m_prXVertex[i].GetNormal() ); 00563 } 00564 } 00565 } 00566 00567 // Recalculate Normals 00568 CalAndSetFaceNormalsNotNormalized(); 00569 CalAndSetVertexNormals(); 00570 NormalizeFaceNormals(); 00571 00572 // FOR CHECKING 00573 //CalVolume2(); 00574 } 00575 //----------------------------------------------------------------------------- 00576 // Preserve Volume 00577 template <typename T> 00578 void VolPresPolygonalModel<T>::PreserveVolume ( T wgt[4] ) 00579 { 00580 //std::cout << "PreserveVolume Fn\n"; 00581 00582 static std::vector<int> deformedVertexList; 00583 00584 if ( deformedVertexList != m_pcVertexRings->GetVertexListOfRingZero() ) { 00585 //---------------------------------------------------------------- 00586 // Find Vertex Rings 00587 00588 // a good one 00589 //const int NoRings = 8; 00590 //T factor[NoRings+1] = { 0, -64, -1, 64, 8, 1, 0.125, 0.015625, 0.002 }; 00591 const int NoRings = 3; 00592 T factor[NoRings+1] = { 1, 1, 1 }; 00593 00594 //const int NoRings = 10; 00595 //T factor[NoRings+1] = { 0, -64, -1, 1024, 256, 64, 16, 4, 1, 0.25, 0.0625 }; 00596 00597 //const int NoRings = 8; 00598 //T factor[NoRings+1] = { 0, -100, -10, 100, 25, 6.25, 1.5625, 0.4, 0.1 }; 00599 00600 //const int NoRings = 5; 00601 //T factor[NoRings+1] = { 0, -1, -1, 1, 2, 10 };` 00602 00603 for ( int i = 1; i <= NoRings; ++i ) { 00604 m_pcVertexRings->IncreaseOneLayer(); 00605 //m_pcVertexRings->IncreaseOneLayer(0); 00606 for ( int j = m_pcVertexRings->GetRingLayers()[i]; j < m_pcVertexRings->GetRingLayers()[i+1]; ++j ) { 00607 m_pvDeformVector[ m_pcVertexRings->GetRingVertices()[j] ] 00608 = m_prXVertex[ m_pcVertexRings->GetRingVertices()[j] ].GetNormal() * factor[i]; 00609 } 00610 } 00611 } 00612 00613 //---------------------------------------------------------------- 00614 // Variables 00615 wgt[0] = wgt[1] = wgt[2] = wgt[3] = 0; 00616 T v[3][3]; 00617 T dv[3][3]; 00618 //---------------------------------------------------------------- 00619 // Set Fix Status of Vertices 00620 for ( int i = 0; i < m_iNoVertices; ++i ) { 00621 m_bFlagFixed[i] = true; 00622 } 00623 //---------------------------------------------------------------- 00624 // Set Unfix status for vertices in the rings, except ring#0 00625 //std::cout << *m_pcVertexRings << "\n"; 00626 for ( int i = 1; i < static_cast<int>( m_pcVertexRings->GetRingVertices().size() ); ++i ) { 00627 m_bFlagFixed[ m_pcVertexRings->GetRingVertices()[i] ] = false; 00628 } 00629 //---------------------------------------------------------------- 00630 // For FPS 00631 static int count = 0; 00632 static clock_t start = clock(); 00633 static clock_t finish; 00634 //start = clock(); 00635 00636 for ( int i = 0; i < m_iNoFaces; ++i ) { 00637 // Tri Face 00638 if ( m_prFace[i].GetNoVertices() == 3 ) { 00639 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][0]; 00640 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][1]; 00641 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ][2]; 00642 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][0]; 00643 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][1]; 00644 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ][2]; 00645 v[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][0]; 00646 v[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][1]; 00647 v[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ][2]; 00648 00649 //if ( true ) { 00650 if ( !m_bFlagFixed[ m_prFace[i].GetVertexNo(0) ] ) { 00651 //dv[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetX(); 00652 //dv[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetY(); 00653 //dv[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal().GetZ(); 00654 00655 //m_pvDeformVector[ m_prFace[i].GetVertexNo(0) ] = m_prXVertex[ m_prFace[i].GetVertexNo(0) ].GetNormal(); 00656 dv[0][0] = m_pvDeformVector[ m_prFace[i].GetVertexNo(0) ].GetX(); 00657 dv[0][1] = m_pvDeformVector[ m_prFace[i].GetVertexNo(0) ].GetY(); 00658 dv[0][2] = m_pvDeformVector[ m_prFace[i].GetVertexNo(0) ].GetZ(); 00659 00660 //dv[0][0] = dv[0][1] = dv[0][2] = 1/3; 00661 } 00662 else { 00663 dv[0][0] = dv[0][1] = dv[0][2] = 0; 00664 } 00665 //if ( true ) { 00666 if ( !m_bFlagFixed[ m_prFace[i].GetVertexNo(1) ] ) { 00667 //dv[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetX(); 00668 //dv[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetY(); 00669 //dv[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal().GetZ(); 00670 00671 //m_pvDeformVector[ m_prFace[i].GetVertexNo(1) ] = m_prXVertex[ m_prFace[i].GetVertexNo(1) ].GetNormal(); 00672 dv[1][0] = m_pvDeformVector[ m_prFace[i].GetVertexNo(1) ].GetX(); 00673 dv[1][1] = m_pvDeformVector[ m_prFace[i].GetVertexNo(1) ].GetY(); 00674 dv[1][2] = m_pvDeformVector[ m_prFace[i].GetVertexNo(1) ].GetZ(); 00675 00676 //dv[1][0] = dv[1][1] = dv[1][2] = 1/3; 00677 } 00678 else { 00679 dv[1][0] = dv[1][1] = dv[1][2] = 0; 00680 } 00681 //if ( true ) { 00682 if ( !m_bFlagFixed[ m_prFace[i].GetVertexNo(2) ] ) { 00683 //dv[2][0] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetX(); 00684 //dv[2][1] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetY(); 00685 //dv[2][2] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal().GetZ(); 00686 00687 //m_pvDeformVector[ m_prFace[i].GetVertexNo(2) ] = m_prXVertex[ m_prFace[i].GetVertexNo(2) ].GetNormal(); 00688 dv[2][0] = m_pvDeformVector[ m_prFace[i].GetVertexNo(2) ].GetX(); 00689 dv[2][1] = m_pvDeformVector[ m_prFace[i].GetVertexNo(2) ].GetY(); 00690 dv[2][2] = m_pvDeformVector[ m_prFace[i].GetVertexNo(2) ].GetZ(); 00691 00692 //dv[2][0] = dv[2][1] = dv[2][2] = 1/3; 00693 } 00694 else { 00695 dv[2][0] = dv[2][1] = dv[2][2] = 0; 00696 } 00697 d3mom0 ( v, dv, wgt ); // call subroutine 00698 } 00699 // Not Tri Face 00700 else { 00701 //---------------------------------------------- 00702 // Find the avg (center) vertex of the face 00703 XVertex<T> avgXVertex; // default ctor --> position and normal are zeros 00704 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00705 avgXVertex.SetPosition( avgXVertex.GetPosition() 00706 + m_prXVertex[m_prFace[i].GetVertexNo(j)].GetPosition() ); 00707 avgXVertex.SetNormal( avgXVertex.GetNormal() 00708 + m_prXVertex[m_prFace[i].GetVertexNo(j)].GetNormal() ); 00709 } 00710 //---------------------------------------------- 00711 int k = m_prFace[i].GetNoVertices() - 1; 00712 for ( int j = 0; j < m_prFace[i].GetNoVertices(); ++j ) { 00713 v[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][0]; 00714 v[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][1]; 00715 v[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ][2]; 00716 v[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][0]; 00717 v[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][1]; 00718 v[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ][2]; 00719 v[2][0] = avgXVertex[0]; 00720 v[2][1] = avgXVertex[1]; 00721 v[2][2] = avgXVertex[2]; 00722 00724 dv[0][0] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ].GetNormal().GetX(); 00725 dv[0][1] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ].GetNormal().GetY(); 00726 dv[0][2] = m_prXVertex[ m_prFace[i].GetVertexNo(k) ].GetNormal().GetZ(); 00727 dv[1][0] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ].GetNormal().GetX(); 00728 dv[1][1] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ].GetNormal().GetY(); 00729 dv[1][2] = m_prXVertex[ m_prFace[i].GetVertexNo(j) ].GetNormal().GetZ(); 00730 dv[2][0] = avgXVertex.GetNormal().GetX(); 00731 dv[2][1] = avgXVertex.GetNormal().GetY(); 00732 dv[2][2] = avgXVertex.GetNormal().GetZ(); 00733 //*/ 00734 00735 d3mom0 ( v, dv, wgt ); // call subroutine 00736 k = j; 00737 } 00738 } 00739 } 00740 00741 // For FPS 00742 finish = clock(); 00743 double duration = static_cast<double>(finish - start); 00744 int dFPS; 00745 //char fps[40]; 00746 if ( duration > 0.0 ) { 00747 dFPS = CLOCKS_PER_SEC / duration; 00748 //sprintf( fps, "FPS: %i", dFPS ); 00749 } 00750 else { 00751 dFPS = 1000; 00752 //sprintf( fps, "FPS: >1000" ); 00753 } 00754 //std::cout << "FPs: " << dFPS << "\n"; 00755 00756 std::cout << "Total of " << ++count << " VolumeCalculations takes " << duration << " msec.\n"; 00757 00758 //---------------------------------------------------------------- 00759 T diffVolume = wgt[0] - m_tVolume; 00760 if ( diffVolume != 0 ) { 00761 //if ( fabs( diffVolume ) > Math<T>::ROUND_ERROR ) { 00762 T sigma = SP_Maths::rootOfCubicPolynomial( wgt[3], wgt[2], wgt[1], diffVolume, 1 ); 00763 //std::cout << "\nsigma: " << sigma << "\n"; 00764 for ( int i = 0; i < m_iNoVertices; ++i ) { 00765 //if ( true ) { 00766 if ( !m_bFlagFixed[i] ) { 00767 //m_prXVertex[i].SetPosition( m_prXVertex[i].GetPosition() + (sigma * m_prXVertex[i].GetNormal()) ); 00768 m_prXVertex[i].SetPosition( m_prXVertex[i].GetPosition() + (sigma * m_pvDeformVector[i]) ); 00769 } 00770 } 00771 } 00772 //---------------------------------------------------------------- 00773 // Recalculate Normals 00774 CalAndSetFaceNormalsNotNormalized(); 00775 CalAndSetVertexNormals(); 00776 NormalizeFaceNormals(); 00777 //---------------------------------------------------------------- 00778 // FOR CHECKING 00779 //CalVolume2(); 00780 } 00781 //----------------------------------------------------------------------------- 00782 /* IN : coefficients v of degree 1 patch 00783 * and coefficients dv of the pertubation 00784 * OUT: weights wgt of moment 0: 00785 * vol = wgt[0]+wgt[1]*s+wgt[2]*s^2+wgt[3]*s^3 00786 * generated by Jorg Peters and Maple */ 00787 //#define PTS 3 00788 //#define XYZ 3 00789 //double d3mom0 (double v[PTS][XYZ], double 00790 // dv[PTS][XYZ], double wgt[4] ) 00791 template <typename T> 00792 inline void VolPresPolygonalModel<T>::d3mom0 ( T v[][3], T dv[][3], T wgt[4] ) 00793 //void VolPresPolygonalModel<T>::d3mom0 () 00794 { 00795 T t[184]; 00796 //------------------------------------------------------------------------- 00797 t[1] = v[1][2]; 00798 t[2] = v[2][0]; 00799 t[3] = t[1]*t[2]; 00800 t[4] = v[1][1]; 00801 t[6] = v[0][2]; 00802 t[7] = t[6]*t[2]; 00803 t[9] = v[1][0]; 00804 t[10] = t[6]*t[9]; 00805 t[11] = v[2][1]; 00806 t[13] = v[0][1]; 00807 t[15] = v[0][0]; 00808 t[16] = t[6]*t[15]; 00809 t[20] = t[1]*t[9]; 00810 t[23] = t[1]*t[15]; 00811 t[27] = v[2][2]; 00812 t[28] = t[27]*t[2]; 00813 t[30] = t[27]*t[9]; 00814 t[33] = t[27]*t[15]; 00815 t[37] = -t[3]*t[4]-t[7]*t[4]+t[10]*t[11]-t[10]*t[13]-t[16]*t[11]+t[7]*t[13]+t[16]*t[4]+t[20]*t[11]-t[20]*t[13]- 00816 t[23]*t[11]+t[3]*t[13]+t[23]*t[4]-t[28]*t[4]+t[30]*t[11]-t[30]*t[13]-t[33]*t[11]+t[28]*t[13]+t[33]*t[4]; 00817 wgt[0] += t[37]/6.0; 00818 //------------------------------------------------------------------------- 00819 t[38] = dv[1][0]; 00820 t[40] = dv[2][0]; 00821 t[43] = dv[1][1]; 00822 t[45] = dv[2][1]; 00823 t[47] = dv[0][1]; 00824 t[51] = dv[0][0]; 00825 t[56] = t[38]*t[11]-t[40]*t[4]+t[40]*t[13]+t[15]*t[43]+t[9]*t[45]-t[9]*t[47]-t[38]*t[13]-t[15]*t[45]-t[51]* 00826 t[11]+t[2]*t[47]-t[2]*t[43]+t[51]*t[4]; 00827 t[58] = dv[0][2]; 00828 t[65] = -t[2]*t[4]+t[9]*t[11]-t[9]*t[13]-t[15]*t[11]+t[2]*t[13]+t[15]*t[4]; 00829 t[68] = dv[1][2]; 00830 t[71] = dv[2][2]; 00831 wgt[1] += t[6]*t[56]/6.0 + t[58]*t[65]/6.0 + t[1]*t[56]/6.0 + t[68]*t[65]/6.0 + t[27]*t[56]/6.0 + 00832 t[71]*t[65]/6.0; 00833 //------------------------------------------------------------------------- 00834 t[74] = t[71]*t[15]; 00835 t[76] = t[68]*t[15]; 00836 t[78] = t[58]*t[40]; 00837 t[81] = t[71]*t[9]; 00838 t[83] = t[68]*t[38]; 00839 t[87] = t[71]*t[2]; 00840 t[89] = t[58]*t[15]; 00841 t[91] = t[68]*t[2]; 00842 t[93] = t[58]*t[2]; 00843 t[96] = t[74]*t[43]-t[76]*t[45]+t[78]*t[13]+t[76]*t[43]-t[81]*t[47]-t[83]*t[13]-t[78]*t[4]+t[83]*t[11]+t[87]* 00844 t[47]-t[89]*t[45]-t[91]*t[43]+t[93]*t[47]+t[81]*t[45]; 00845 t[98] = t[58]*t[38]; 00846 t[100] = t[71]*t[51]; 00847 t[102] = t[6]*t[38]; 00848 t[104] = t[27]*t[51]; 00849 t[106] = t[27]*t[38]; 00850 t[108] = t[27]*t[40]; 00851 t[111] = t[6]*t[51]; 00852 t[113] = t[1]*t[38]; 00853 t[115] = t[71]*t[40]; 00854 t[121] = t[98]*t[11]-t[100]*t[11]-t[102]*t[47]-t[104]*t[45]+t[106]*t[45]-t[108]*t[43]+t[108]*t[47]+t[111] 00855 *t[43]-t[113]*t[47]-t[115]*t[4]+t[91]*t[47]-t[87]*t[43]+t[89]*t[43]+t[113]*t[45]; 00856 t[124] = t[68]*t[40]; 00857 t[126] = t[68]*t[51]; 00858 t[130] = t[1]*t[51]; 00859 t[133] = t[6]*t[40]; 00860 t[136] = t[71]*t[38]; 00861 t[142] = -t[124]*t[4]+t[126]*t[4]+t[100]*t[4]-t[74]*t[45]-t[130]*t[45]+t[102]*t[45]+t[133]*t[47]-t[93]* 00862 t[43]-t[136]*t[13]+t[130]*t[43]-t[98]*t[13]-t[111]*t[45]-t[133]*t[43]; 00863 t[144] = t[58]*t[9]; 00864 t[147] = t[58]*t[51]; 00865 t[150] = t[1]*t[40]; 00866 t[154] = t[68]*t[9]; 00867 t[162] = t[144]*t[45]-t[144]*t[47]-t[147]*t[11]+t[147]*t[4]+t[150]*t[47]-t[150]*t[43]+t[124]*t[13]+t[154] 00868 *t[45]-t[154]*t[47]-t[126]*t[11]+t[104]*t[43]-t[106]*t[47]+t[136]*t[11]+t[115]*t[13]; 00869 wgt[2] += t[162]/6.0 + t[142]/6.0 + t[121]/6.0 + t[96]/6.0; 00870 //------------------------------------------------------------------------- 00871 t[183] = t[78]*t[47]+t[147]*t[43]+t[98]*t[45]-t[98]*t[47]-t[147]*t[45]-t[78]*t[43]+t[124]*t[47]+t[126]* 00872 t[43]+t[83]*t[45]-t[83]*t[47]-t[126]*t[45]-t[124]*t[43]+t[115]*t[47]+t[100]*t[43]+t[136]*t[45]-t[136]*t[47]-t[100]* 00873 t[45]-t[115]*t[43]; 00874 wgt[3] += t[183]/6.0; 00875 //------------------------------------------------------------------------- 00876 } 00877 //----------------------------------------------------------------------------- 00878 //----------------------------------------------------------------------------- 00879 // StepSimulation 00880 template <typename T> 00881 void VolPresPolygonalModel<T>::StepSimulation( T dt ) { 00882 static int stepNo = 0; 00883 00884 std::cout << "SimStep#" << stepNo++ << "\n"; 00885 00886 T divider = 10.0; 00887 00888 static int deformNo = 0; 00889 00890 // Deform Vertices 00891 const int size = 3; 00892 //if ( 1799 < m_iNoFaces ) { 00893 m_iDeformedFaceNo = 1; 00894 //} 00895 int deformedFaceNo = m_iDeformedFaceNo; 00896 int deformVertexList[size]; 00897 Vector3<T> deformVector; 00898 for ( int i = 0; i < size; ++i ) { 00899 deformVertexList[i] = m_prFace[deformedFaceNo].GetVertexNo( i ); 00900 } 00901 m_pcVertexRings->SetVertices( deformVertexList, size ); 00902 00903 if ( ++deformNo <= 75 ) { 00904 for ( int i = 0; i < size; ++i ) { 00905 DeformVertex( deformVertexList[i], Vector3<T>( m_prXVertex[deformVertexList[i]].GetPosition() * -0.01 ) ); 00906 } 00907 //m_pcVertexRings->SetExternalVertexRingInfo( m_pviVertexRing1List ); 00908 } 00909 else if ( ++deformNo <= 150 ) { 00910 for ( int i = 0; i < size; ++i ) { 00911 DeformVertex( deformVertexList[i], Vector3<T>( m_prXVertex[deformVertexList[i]].GetPosition() * +0.01 ) ); 00912 } 00913 //m_pcVertexRings->SetExternalVertexRingInfo( m_pviVertexRing1List ); 00914 } 00915 else { 00916 } 00917 00918 00919 //CalForces( Vector3<T>( (rand()%1000)/divider, (rand()%10)/divider, (rand()%10)/divider ) ); 00920 00921 if ( m_bFlagSpringForceMode ) { 00922 CalForces(); 00923 // Euler's Method for Integration for ODE Solver 00924 for ( int n = 0; n < m_iNoVertices; ++n ) { 00925 m_pprParticleVertex[n]->SetAcceleration( m_pprParticleVertex[n]->GetForce() / m_pprParticleVertex[n]->GetMass() ); 00926 m_pprParticleVertex[n]->SetVelocity( m_pprParticleVertex[n]->GetVelocity() 00927 + ( m_pprParticleVertex[n]->GetAcceleration() * dt ) ); 00928 m_pprParticleVertex[n]->SetPosition( m_pprParticleVertex[n]->GetPosition() 00929 + ( m_pprParticleVertex[n]->GetVelocity() * dt ) ); 00930 } 00931 } 00932 00933 //PreserveVolumeByVertexNormals ( weight ); 00934 if ( m_bFlagVolPresMode ) { 00935 if ( false ) { 00936 //if ( deformNo <= 25 ) { 00937 PreserveVolume( m_weight ); 00938 } 00939 else { 00940 PreserveVolumeByVertexNormals ( m_weight ); 00941 //PreserveVolume( weight ); 00942 } 00943 } 00944 } 00945 //----------------------------------------------------------------------------- 00946 // Calculate All Forces Acting on Particles 00947 template <typename T> 00948 void VolPresPolygonalModel<T>::CalForces( Vector3<T> vWind, T g ) 00949 { 00950 //CalParticleNormals( 0 ); // calculate particle normals use for finding force due to wind 00951 static int simCounter = 0; 00952 00953 T zero = 0, y = 0; 00954 for ( int n = 0; n < m_iNoVertices; ++n ) { 00955 if ( m_pprParticleVertex[n]->GetFixStatus() ) { // If the particle is fixed, skip it! 00956 m_pprParticleVertex[n]->SetForce( 0, 0, 0 ); 00957 continue; 00958 } 00959 // Clear and Calculate Force due to Acceleration of Gravity -------- 00960 // Fij = mass * g; where g is a vector 00961 y = m_pprParticleVertex[n]->GetMass() * (-g); 00962 00963 /* 00964 m_pprParticleVertex[n]->SetForce( 0, y, 0 ); 00965 */ 00966 00967 /* 00968 if ( simCounter++ < 10 && ( n == 0 || n == 4 ) ) { 00969 m_pprParticleVertex[n]->SetPosition( m_pprParticleVertex[n]->GetPosition() * 0.9999999999 ); 00970 } 00971 */ 00972 00973 // Add Force due to Viscous Damping -------------------------------- 00974 // Fij = -Kvd * Vij 00975 m_pprParticleVertex[n]->SetForce( m_pprParticleVertex[n]->GetForce() 00976 + 0.001/*m_tKvd (viscous damping)*/ * m_pprParticleVertex[n]->GetVelocity() ); 00977 00978 // Add Force due to viscous fluid moving (wind) -------------------- 00979 // Fij = Kvfm * [Nij * (vWind - Vij)] * Nij 00980 // where Nij is the unit normal on the surface at particle Pij 00981 m_pprParticleVertex[n]->SetForce( m_pprParticleVertex[n]->GetForce() 00982 + 0.1/*m_tKvfm (viscous fluid moving damping)*/ 00983 * ( m_pprParticleVertex[n]->GetNormal() * ( vWind - m_pprParticleVertex[n]->GetVelocity() ) ) 00984 * m_pprParticleVertex[n]->GetNormal() ); 00985 00986 // Dump 00987 //cout << "Mass: " << m_pcParticle[m][n]->m_tMass << " Force: " << m_pcParticle[m][n]->m_vtForce << "\n"; 00988 } 00989 00990 AddSpringForces(); // calculate and add spring forces to particle forces 00991 } 00992 //----------------------------------------------------------------------------- 00993 // Calculate and Add the Force of a Spring 00994 template <typename T> 00995 inline void VolPresPolygonalModel<T>::AddSpringForces() 00996 { 00997 for ( int i = 0; i < m_iNoEdges; ++i ) { 00998 01000 //------------------------------------------------------------ 01001 // Local Variables 01002 // spring current length vector 01003 Vector3<T> cL = ( m_pprSpringVertex[i]->GetParticleVertexOne().GetPosition() 01004 - m_pprSpringVertex[i]->GetParticleVertexTwo().GetPosition() ); 01005 // spring current length scalar 01006 T scl = cL.Length(); 01007 // velocity difference of the two particles linked by the spring 01008 Vector3<T> V = ( m_pprSpringVertex[i]->GetParticleVertexOne().GetVelocity() 01009 - m_pprSpringVertex[i]->GetParticleVertexTwo().GetVelocity() ); 01010 01011 // F1 = -{Ks(l - r) + Kd[(V1-V2)*L]/l}*L/l 01012 // F2 = -F1 01013 // where Ks = spring stiffness, Kd = spring damping, 01014 // V1 and V2 are velocity of particles linked by the spring 01015 // r = spring rest length 01016 // L = vector of difference of positions of particle#1 and #2 01017 // l = magnitude of L 01018 Vector3<T> f = -( 01019 ( m_pprSpringVertex[i]->GetConstantK() * ( scl - m_pprSpringVertex[i]->GetRestLengthL() ) ) 01020 + ( m_pprSpringVertex[i]->GetDampingD() * ((V*cL) / scl) ) ) 01021 * (cL / scl); 01023 //*/ 01024 01025 //std::cout << "Spring Force: " << f << std::endl; 01026 /* 01027 if ( !m_pprSpringVertex[i]->GetParticleVertexOne().GetFixStatus() ) 01028 m_pprSpringVertex[i]->GetParticleVertexOne().SetForce( f ); 01029 if ( !m_pprSpringVertex[i]->GetParticleVertexTwo().GetFixStatus() ) 01030 m_pprSpringVertex[i]->GetParticleVertexTwo().SetForce( -f ); 01031 */ 01033 if ( !m_pprSpringVertex[i]->GetParticleVertexOne().GetFixStatus() ) 01034 m_pprSpringVertex[i]->GetParticleVertexOne().SetForce( 01035 m_pprSpringVertex[i]->GetParticleVertexOne().GetForce() 01036 + m_pprSpringVertex[i]->GetForce() 01037 + f 01038 ); 01039 if ( !m_pprSpringVertex[i]->GetParticleVertexTwo().GetFixStatus() ) 01040 m_pprSpringVertex[i]->GetParticleVertexTwo().SetForce( 01041 m_pprSpringVertex[i]->GetParticleVertexTwo().GetForce() 01042 - m_pprSpringVertex[i]->GetForce() 01043 - f 01044 ); 01045 //*/ 01046 } 01047 } 01048 //----------------------------------------------------------------------------- 01049 // Defrom A Vertex # v by a deformVector 01050 template <typename T> 01051 inline void VolPresPolygonalModel<T>::DeformVertex ( int v, Vector3<T> & deformVector ) 01052 { 01053 m_prXVertex[v].SetPosition( m_prXVertex[v].GetPosition() + deformVector ); 01054 } 01055 //----------------------------------------------------------------------------- 01056 // Defrom Vertices referred by vList by a deformVector 01057 template <typename T> 01058 inline void VolPresPolygonalModel<T>::DeformVertices ( int vList[], int size, Vector3<T> & deformVector ) 01059 { 01060 for ( int i = 0; i < size; ++i ) { 01061 m_prXVertex[ vList[i] ].SetPosition( m_prXVertex[ vList[i] ].GetPosition() + deformVector ); 01062 } 01063 } 01064 //----------------------------------------------------------------------------- 01065 01066 //****************************************************************************** 01067 // Transformation by Model's translation, rotation, and scale 01068 //****************************************************************************** 01069 //----------------------------------------------------------------------------- 01070 template <typename T> 01071 void VolPresPolygonalModel<T>::ApplyAndResetTransform () 01072 { 01073 XPolygonalModel<T>::ApplyAndResetTransform(); 01074 //---------------------------------------------------------------- 01075 // Recalculate Normals 01076 CalAndSetFaceNormalsNotNormalized(); 01077 CalAndSetVertexNormals(); 01078 NormalizeFaceNormals(); 01079 // Calculate Volume 01080 CalVolume2(); 01081 } 01082 //----------------------------------------------------------------------------- 01083 //****************************************************************************** 01084 01085 //----------------------------------------------------------------------------- 01086 //============================================================================= 01087 END_NAMESPACE_TAPs__OpenGL 01088 //----------------------------------------------------------------------------- 01089 //345678901234567890123456789012345678901234567890123456789012345678901234567890 01090 //--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
1.7.4 
