COME - Biomodel

00001 /***************************************************************************
00002  *   Copyright (C) 2002 by Anderson Maciel                                 *
00003  *   andi.maciel@gmail.com                                                 *
00004  *                                                                         *
00005  *   This program is free software; you can redistribute it and/or modify  *
00006  *   it under the terms of the GNU General Public License as published by  *
00007  *   the Free Software Foundation; either version 2 of the License, or     *
00008  *   (at your option) any later version.                                   *
00009  *                                                                         *
00010  *   This program is distributed in the hope that it will be useful,       *
00011  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
00012  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
00013  *   GNU General Public License for more details.                          *
00014  *                                                                         *
00015  *   You should have received a copy of the GNU General Public License     *
00016  *   along with this program; if not, write to the                         *
00017  *   Free Software Foundation, Inc.,                                       *
00018  *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
00019  **************************************************************************/
00030 
00031 #include <algebra/comeface.h>
00032 #include <algebra/comemesh.h>
00033 #include <algebra/comematrix.h>
00034 #include <bio/comebiostructure.h>
00035 #include <math.h>
00036 #include <bio/comemolecule.h>
00037 
00041 COME_Face :: COME_Face(){
00042 
00043         color[0] = 1;
00044         color[1] = .7;
00045         color[2] = .3;
00046         parent = NULL;
00047         nearestMolecule = NULL;
00048         collisionDetectable = true;
00049 }
00050 
00051 /*COME_Face :: COME_Face( const COME_Face& pFace ){
00052 
00053         normal = pFace.getNormal();
00054         normalGlobalPosition = pFace.getNormalGlobalPosition();
00055         nearestMolecule = pFace.getNearestMolecule();
00056         neighbourFaces = pFace.getNeighbourFaces();
00057         color[0] = pFace.getColor()[0];
00058         color[1] = pFace.getColor()[1];
00059         color[2] = pFace.getColor()[2];
00060         color[3] = pFace.getColor()[3];
00061         parent = pFace.getParent();
00062 }*/
00063 
00064 COME_Face :: COME_Face( vector<int> verticesN, COME_Mesh* parentN ){
00065 
00066         vertices = verticesN;
00067         parent  = parentN;
00068         updateNormal();
00069         nearestMolecule = NULL;
00070         color[0] = 1;
00071         color[1] = .7;
00072         color[2] = .3;  
00073         collisionDetectable = true;
00074 }
00075 
00076 COME_Face :: COME_Face( vector<int> verticesN ){
00077 
00078         vertices = verticesN;
00079         if( parent ){
00080                 updateNormal();
00081         }
00082         nearestMolecule = NULL;
00083         color[0] = 1;
00084         color[1] = .7;
00085         color[2] = .3;  
00086         collisionDetectable = true;
00087 }
00088 
00089 COME_Face :: ~COME_Face(){
00090 
00091         
00092 }
00093 
00097 void            
00098 COME_Face :: setVertex( vector<int> verticesN ) {
00099         
00100         vertices = verticesN;
00101 }
00102 
00103 void
00104 COME_Face :: setNormal( COME_Vector3D normalN ){
00105 
00106         normal = normalN;
00107 }
00108 
00109 void
00110 COME_Face::setNearestMolecule( COME_Molecule *nearest ){
00111 
00112         nearestMolecule = nearest;
00113 }
00114 
00115 void 
00116 COME_Face::setColor(double* c){
00117         
00118         color[0] = c[0];
00119         color[1] = c[1];
00120         color[2] = c[2];
00121         color[3] = c[3];
00122 }
00123 
00124 void
00125 COME_Face::setCollisionDetectable( bool yesno ){
00126         
00127         collisionDetectable = yesno;
00128 }
00129         
00130 
00134 
00135 int
00136 COME_Face :: getIndexVertex( int iPosN ) const{
00137 
00138         return vertices[ iPosN ];
00139 }
00140 
00141 vector<int>
00142 COME_Face:: getIndexVertices(){
00143         return vertices;
00144 }
00145 int     
00146 COME_Face :: getNumberVertices() const{
00147         
00148         return vertices.size();
00149 }
00150 
00151 bool
00152 COME_Face::getCollisionDetectable(){
00153         
00154         return collisionDetectable;
00155 }
00156 
00157 COME_Vertex3D   
00158 COME_Face :: getVertex( int iPosN ) const{
00159 
00160         if( parent)
00161                 return ((COME_Mesh*)parent)->getAVertex( vertices[ iPosN ] );
00162         printf( "COME_Face::getVertex(int): A parent mesh is not set for a face trying to return a vertex.\n" );
00163         exit(0);
00164         
00165 }
00166 
00167 COME_Vertex3D*
00168 COME_Face :: getClosestVertexPt( COME_Point3D currentPoint ){
00169 
00170         double distance = INFINITE;
00171         int closest = 0;
00172         vector<COME_Vertex3D*> globalVertices = getVerticesGlobalPositionPt();
00173         for( int i = 0; i < globalVertices.size(); i++ ){
00174                 
00175                 if( distance > globalVertices[i]->getPoint3D().vpDistance( currentPoint ) ){
00176                         distance = globalVertices[i]->getPoint3D().vpDistance( currentPoint );
00177                         closest = i;
00178                 }
00179         }
00180         return globalVertices[closest];
00181 }
00182 
00183 double
00184 COME_Face :: geoSolidAngle( COME_Point3D point2Check ){
00185 
00186         double area = 0.0, ang, s, l1, l2;
00187         COME_Vector3D p1, p2, r1, a, b, n1, n2, plane;
00188         
00189         vector<COME_Vertex3D*> globalVertices = getVerticesGlobalPositionPt();
00190         plane = getNormalGlobalPositionPt()->getPoint3D();
00191         p2 = getVertexPt(2)->getPoint3D();
00192         p1 = getVertexPt(0)->getPoint3D();
00193         a = p2 - p1;
00194         
00195         for( int i = 0; i < vertices.size(); i++   ){
00196                 
00197                 r1 = p1 - point2Check;
00198                 p2 = globalVertices[(i+1)%vertices.size()]->getPoint3D();
00199                 b = p2 - p1;
00200                 n1 = a * r1;
00201                 n2 = r1 * b;
00202                 
00203                 l1 = n1.vpModule();
00204                 l2 = n2.vpModule();
00205                 s = n1.vpDotProduct(n2) / (l1 * l2 );
00206                 ang = acos( max( -1.0, min( 1.0, s ) ) );
00207                 s = ( b * a ).vpDotProduct( plane );
00208                 area += s > 0.0 ? M_PI - ang : M_PI + ang;
00209                 
00210                 a = b * -1.0;
00211                 p1 = p2;
00212         }
00213         
00214         area -= M_PI*(vertices.size()-2);
00215         
00216         return (plane.vpDotProduct(r1)> 0.0) ? -area : area;
00217 }
00218 
00219 COME_Vertex3D*  
00220 COME_Face :: getVertexPt( int iPosN ) {
00221 
00222         if( parent)
00223                 return ((COME_Mesh*)parent)->getAVertexPt( vertices[ iPosN ] );
00224         printf( "COME_Face::getVertexPt(int): A parent mesh is not set for a face trying to return a vertex.\n" );
00225         exit(0);
00226 }
00227 
00228 
00229 COME_Vertex3D   
00230 COME_Face :: getVertexGlobalPosition( int iPosN ) const{
00231         return ((COME_Mesh*)parent)->getAVertexGlobalPosition( vertices[ iPosN ] );
00232 }
00233 
00234 COME_Vertex3D*  
00235 COME_Face :: getVertexGlobalPositionPt( int iPosN ){
00236         return ((COME_Mesh*)parent)->getAVertexGlobalPositionPt( vertices[ iPosN ] );
00237 }
00238 
00239 COME_Vector3D
00240 COME_Face :: getNormal() const{
00241 
00242         return normal;
00243 }
00244 
00245 COME_Vector3D
00246 COME_Face :: getNormalGlobalPosition() const {
00247 
00248         return normalGlobalPosition;
00249 }
00250 
00251 COME_Vector3D*
00252 COME_Face :: getNormalGlobalPositionPt(){
00253 
00254         return &normalGlobalPosition;
00255 }
00256 
00257 const double*
00258 COME_Face::getColor() const {
00259         return color;
00260 }
00261 
00265 void
00266 COME_Face::updateNormal(){
00267 
00268         COME_Vector3D v1; v1.setVector3D( ((COME_Mesh*)parent)->getAVertex(vertices[1]).getPoint3D() - ((COME_Mesh*)parent)->getAVertex(vertices[0]).getPoint3D() );
00269         COME_Vector3D v2; v2.setVector3D( ((COME_Mesh*)parent)->getAVertex(vertices[2]).getPoint3D() - ((COME_Mesh*)parent)->getAVertex(vertices[0]).getPoint3D() );
00270         COME_Vector3D v3 = v1 * v2;
00271         v3.vpNormalize();
00272         normal = COME_Vector3D( v3 );
00273 
00274 }
00275 
00276 void
00277 COME_Face::updateNormalGlobalPosition(){
00278         
00279         vector<COME_Vertex3D*> globalVertices = getVerticesGlobalPositionPt();
00280         COME_Vector3D v1; v1.setVector3D(globalVertices[1]->getPoint3D() - globalVertices[0]->getPoint3D() );
00281         COME_Vector3D v2; v2.setVector3D(globalVertices[2]->getPoint3D() - globalVertices[0]->getPoint3D() );
00282         COME_Vector3D v3 = (v1 * v2);
00283         v3.vpNormalize();
00284         normalGlobalPosition = v3;
00285 }
00286 
00287 
00288 vector<COME_Vertex3D>   
00289 COME_Face :: getVertices() const{
00290 
00291         vector<COME_Vertex3D> verticesN;
00292         for ( int i = 0; i < vertices.size(); i++ ){
00293                 
00294                 verticesN.push_back( getVertex( i ));           
00295         
00296         }
00297 
00298         return verticesN;
00299 }
00300 
00301 
00302 vector<COME_Vertex3D>   
00303 COME_Face :: getVerticesGlobalPosition() const{
00304 
00305         vector<COME_Vertex3D> verticesN;
00306         for ( int i = 0; i < vertices.size(); i++ ){
00307                 
00308                 verticesN.push_back( getVertexGlobalPosition( i ));             
00309         
00310         }
00311 
00312         return verticesN;
00313 }
00314 
00315 vector<COME_Vertex3D*>  
00316 COME_Face :: getVerticesGlobalPositionPt(){
00317 
00318         vector<COME_Vertex3D*> verticesN;
00319         for ( int i = 0; i < vertices.size(); i++ ){
00320                 
00321                 verticesN.push_back( getVertexGlobalPositionPt( i ));           
00322         
00323         }
00324 
00325         return verticesN;
00326 }
00327 
00328 COME_Molecule*
00329 COME_Face::getNearestMolecule() const {
00330 
00331         return nearestMolecule;
00332 }
00333 
00339 double
00340 COME_Face::getElasticity(){
00341 
00342         return ( nearestMolecule->getMaterial()->getYoungsModulus() );
00343 }
00344 
00345 
00349 COME_Point3D
00350 COME_Face::getCenter(){
00351 
00352         COME_Vector3D center;
00353         vector<COME_Vertex3D> verticesFace  = getVertices();
00354 
00355         for ( int i = 0; i < verticesFace.size(); i++ ){
00356                 
00357                 center.setX( center.getX() + verticesFace[ i ].getX() );
00358                 center.setY( center.getY() + verticesFace[ i ].getY() );
00359                 center.setZ( center.getZ() + verticesFace[ i ].getZ() );
00360         }
00361 
00362         center /= 3.0;
00363 
00364         return center.getPoint3D();
00365 
00366 }
00367 
00368 COME_Vertex3D
00369 COME_Face::getCenterGlobalPosition(){
00370 
00371         COME_Vertex3D center;
00372         vector<COME_Vertex3D> verticesFace  = getVerticesGlobalPosition();
00373 
00374         for ( int i = 0; i < verticesFace.size(); i++ ){
00375                 
00376                 center.x = center.x + verticesFace[ i ].x;
00377                 center.y = center.y + verticesFace[ i ].y;
00378                 center.z = center.z + verticesFace[ i ].z;
00379         }
00380 
00381         center.x /= 3.0;
00382         center.y /= 3.0;
00383         center.z /= 3.0;
00384 
00385         return center;
00386 
00387 }
00388 
00389 COME_Vertex3D*
00390 COME_Face::getCenterGlobalPositionPt(){
00391 
00392         COME_Vertex3D *center = new COME_Vertex3D();
00393         vector<COME_Vertex3D*> verticesFace  = getVerticesGlobalPositionPt();
00394 
00395         for ( int i = 0; i < verticesFace.size(); i++ ){
00396                 
00397                 center->x = center->x + verticesFace[ i ]->x;
00398                 center->y = center->y + verticesFace[ i ]->y;
00399                 center->z = center->z + verticesFace[ i ]->z;
00400         }
00401 
00402         center->x /= 3.0;
00403         center->y /= 3.0;
00404         center->z /= 3.0;
00405 
00406         return center;
00407 
00408 }
00409 
00414 double
00415 COME_Face :: getArea(){
00416 
00417         vector<COME_Vertex3D> verticesFace  = getVertices();
00418         double area = ( ( COME_Vector3D( ( verticesFace[1] - verticesFace[0] ) + ( verticesFace[2] - verticesFace[0] ) ).vpModule() / 2.0 ) * ( COME_Vector3D( verticesFace[2] - verticesFace[1] ).vpModule() ) ) / 2.0;
00419         return area;
00420 
00421 }
00422 
00426 double
00427 COME_Face :: triangleArea( COME_Point3D A, COME_Point3D B, COME_Point3D C ){
00428 
00429         COME_Vector3D bNormal = C - A; bNormal.vpNormalize();
00430         COME_Vector3D cNormal = B - A; cNormal.vpNormalize();
00431         
00432         double area = 0.5* COME_Vector3D( C - A ).vpModule() * COME_Vector3D( B - A ).vpModule() * sin( acos( bNormal.vpDotProduct( cNormal ) ) );
00433         //printf( "area: %f\n", area );
00434         return area;
00435 }
00436 
00440 double
00441 COME_Face :: getStress(){
00442 
00443         vector<COME_Vertex3D> verticesFace  = getVertices();
00444         double stress = ( verticesFace[0].getStress() + verticesFace[1].getStress() + verticesFace[2].getStress() ) / 3.0;
00445         return stress;
00446 
00447 }
00448 
00454 double                          
00455 COME_Face :: distancePointFace( COME_Point3D pPos ){
00456 
00457         COME_Vector3D p_a = pPos - COME_Vector3D( getVertexPt( 0 )->x, getVertexPt( 0 )->y, getVertexPt( 0 )->z );
00458         double distPointPlan = fabs( normal.vpDotProduct( p_a ) );
00459 
00460         COME_Point3D pPosPr = pPos - ( normal * distPointPlan );
00461         
00462         if( !isInside2D( pPosPr ) ){
00463                 
00464                 COME_Point3D pA, pB, pC;
00465                 pA = COME_Vector3D( getVertexPt( 0 )->x, getVertexPt( 0 )->y, getVertexPt( 0 )->z );
00466                 pB = COME_Vector3D( getVertexPt( 1 )->x, getVertexPt( 1 )->y, getVertexPt( 1 )->z );
00467                 pC = COME_Vector3D( getVertexPt( 2 )->x, getVertexPt( 2 )->y, getVertexPt( 2 )->z );
00468                 
00469                 COME_Vector3D a_b = pB  - pA;
00470                 a_b.vpNormalize();
00471                 COME_Vector3D b_c = pC  - pB;
00472                 b_c.vpNormalize();
00473                 COME_Vector3D c_a = pA  - pC;
00474                 c_a.vpNormalize();
00475                 
00477                 COME_Point3D pOnEdgeAB = pA + ( a_b * ( (COME_Vector3D( pPosPr - pA )).vpDotProduct(a_b) )  );
00478                 COME_Point3D pOnEdgeBC = pB + ( b_c * ( (COME_Vector3D( pPosPr - pB )).vpDotProduct(b_c) )  );
00479                 COME_Point3D pOnEdgeCA = pC + ( c_a * ( (COME_Vector3D( pPosPr - pC )).vpDotProduct(c_a) )  );
00480                 
00482                 if( ( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() + (COME_Vector3D(pB-pOnEdgeAB)).vpModule() ) > (COME_Vector3D(pA-pB)).vpModule() ){
00483                         if( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() < (COME_Vector3D(pB-pOnEdgeAB)).vpModule() )
00484                                 pOnEdgeAB = pA;
00485                         else pOnEdgeAB = pB;
00486                 }
00487                 if( ( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() + (COME_Vector3D(pC-pOnEdgeBC)).vpModule() ) > (COME_Vector3D(pB-pC)).vpModule() ){
00488                         if( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() < (COME_Vector3D(pC-pOnEdgeBC)).vpModule() )
00489                                 pOnEdgeBC = pB;
00490                         else pOnEdgeBC = pC;
00491                 }
00492                 if( ( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() + (COME_Vector3D(pA-pOnEdgeCA)).vpModule() ) > (COME_Vector3D(pC-pA)).vpModule() ){
00493                         if( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() < (COME_Vector3D(pA-pOnEdgeCA)).vpModule() )
00494                                 pOnEdgeCA = pC;
00495                         else pOnEdgeCA = pA;
00496                 }
00497                 
00498                 double distProj = pPos.vpDistance( pOnEdgeAB );
00499                 if( distProj > pPos.vpDistance( pOnEdgeBC ) ){
00500                         distProj = pPos.vpDistance( pOnEdgeBC );
00501                         if( distProj > pPos.vpDistance( pOnEdgeCA ) ){
00502                                 distProj = pPos.vpDistance( pOnEdgeCA );
00503                         }
00504                 }
00505                 return distProj;
00506         }
00507         
00508         return distPointPlan;
00509 }
00510 
00511 
00515 void
00516 COME_Face :: getClosestEdge( COME_Point3D &pA, COME_Point3D &pB, COME_Point3D p ){
00517 
00518         pA = COME_Point3D( getVertexPt( 0 )->x, getVertexPt( 0 )->y, getVertexPt( 0 )->z );
00519         pB = COME_Point3D( getVertexPt( 1 )->x, getVertexPt( 1 )->y, getVertexPt( 1 )->z );
00520         COME_Point3D pC = COME_Point3D( getVertexPt( 2 )->x, getVertexPt( 2 )->y, getVertexPt( 2 )->z );
00521         
00522         if( p.vpDistance( pA ) < p.vpDistance( pB ) ){
00523                 if( p.vpDistance( pB ) < p.vpDistance( pC ) ){
00524                         return;
00525                 } else {
00526                         pB = pC;
00527                         return;
00528                 }
00529         } else if ( p.vpDistance( pB ) < p.vpDistance( pC ) ){
00530                 if( p.vpDistance( pC ) < p.vpDistance( pA ) ){
00531                         pA = pC;
00532                         return;
00533                 } else {
00534                         return;
00535                 }
00536         }       
00537         
00538 }
00539 
00543 bool
00544 COME_Face :: isInside2D( COME_Point3D pPos ){
00545 
00546         //find pA
00547         COME_Point3D pA = COME_Point3D( getVertexPt( 0 )->x, getVertexPt( 0 )->y, getVertexPt( 0 )->z );
00548         COME_Point3D pB = COME_Point3D( getVertexPt( 1 )->x, getVertexPt( 1 )->y, getVertexPt( 1 )->z );
00549         COME_Point3D pC = COME_Point3D( getVertexPt( 2 )->x, getVertexPt( 2 )->y, getVertexPt( 2 )->z );
00550         
00551         // build vectors
00552         COME_Vector3D c_b = pC  - pB;
00553         COME_Vector3D p_b = pPos  - pB;
00554         COME_Vector3D a_b = pA  - pB;
00555         COME_Vector3D c_a = pC  - pA;
00556         COME_Vector3D p_a = pPos  - pA;
00557         COME_Vector3D b_a = pB  - pA;
00558         
00559         if( ( c_b * p_b ).vpDotProduct( c_b * a_b ) < 0 ){
00560                 return false;
00561         } else if( ( c_a * p_a ).vpDotProduct( c_a * b_a ) < 0 ){
00562                 return false;
00563         } else if( ( b_a * p_a ).vpDotProduct( b_a * c_a ) < 0 ){
00564                 return false;
00565         } 
00566         return true;
00567 }
00568 
00575 double                          
00576 COME_Face :: distancePointFaceGlobal( COME_Point3D pPos ){
00577 
00578         COME_Vector3D p_a = pPos - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00579         double distPointPlan = normalGlobalPosition.vpDotProduct( p_a );
00580 
00581         COME_Point3D pPosPr = pPos - ( normalGlobalPosition * distPointPlan );
00582         
00583         distPointPlan = fabs( distPointPlan );
00584         
00585         if( !isInside2DGlobal( pPosPr ) ){
00586                 
00587                 COME_Point3D pA, pB, pC;
00588                 pA = COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00589                 pB = COME_Vector3D( getVertexGlobalPositionPt( 1 )->x, getVertexGlobalPositionPt( 1 )->y, getVertexGlobalPositionPt( 1 )->z );
00590                 pC = COME_Vector3D( getVertexGlobalPositionPt( 2 )->x, getVertexGlobalPositionPt( 2 )->y, getVertexGlobalPositionPt( 2 )->z );
00591                 
00592                 COME_Vector3D a_b = pB  - pA;
00593                 a_b.vpNormalize();
00594                 COME_Vector3D b_c = pC  - pB;
00595                 b_c.vpNormalize();
00596                 COME_Vector3D c_a = pA  - pC;
00597                 c_a.vpNormalize();
00598                 
00600                 COME_Point3D pOnEdgeAB = pA + ( a_b * ( (COME_Vector3D( pPosPr - pA )).vpDotProduct(a_b) )  );
00601                 COME_Point3D pOnEdgeBC = pB + ( b_c * ( (COME_Vector3D( pPosPr - pB )).vpDotProduct(b_c) )  );
00602                 COME_Point3D pOnEdgeCA = pC + ( c_a * ( (COME_Vector3D( pPosPr - pC )).vpDotProduct(c_a) )  );
00603                 
00605                 double EPSILON = 0.0000000001;
00606                 if( ( ( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() + EPSILON ) > (COME_Vector3D(pA-pB)).vpModule() )
00607                 ||  ( ( (COME_Vector3D(pB-pOnEdgeAB)).vpModule() + EPSILON ) > (COME_Vector3D(pA-pB)).vpModule() ) ){
00608                         if( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() < (COME_Vector3D(pB-pOnEdgeAB)).vpModule() )
00609                                 pOnEdgeAB = pA;
00610                         else pOnEdgeAB = pB;
00611                 }
00612                 if( ( ( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() + EPSILON ) > (COME_Vector3D(pB-pC)).vpModule() )
00613                 ||  ( ( (COME_Vector3D(pC-pOnEdgeBC)).vpModule() + EPSILON ) > (COME_Vector3D(pB-pC)).vpModule() ) ){
00614                         if( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() < (COME_Vector3D(pC-pOnEdgeBC)).vpModule() )
00615                                 pOnEdgeBC = pB;
00616                         else pOnEdgeBC = pC;
00617                 }
00618                 if( ( ( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() + EPSILON ) > (COME_Vector3D(pC-pA)).vpModule() )
00619                 ||  ( ( (COME_Vector3D(pA-pOnEdgeCA)).vpModule() + EPSILON ) > (COME_Vector3D(pC-pA)).vpModule() ) ){
00620                         if( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() < (COME_Vector3D(pA-pOnEdgeCA)).vpModule() )
00621                                 pOnEdgeCA = pC;
00622                         else pOnEdgeCA = pA;
00623                 }
00624         
00625                 double distProj = pPos.vpDistance( pOnEdgeAB );
00626                 if( distProj > pPos.vpDistance( pOnEdgeBC ) ){
00627                         distProj = pPos.vpDistance( pOnEdgeBC );
00628                 }
00629                 if( distProj > pPos.vpDistance( pOnEdgeCA ) ){
00630                         distProj = pPos.vpDistance( pOnEdgeCA );
00631                 }
00632                 return distProj;
00633         }
00634         return distPointPlan;
00635 }
00636 
00643 double                          
00644 COME_Face :: signedDistancePointFaceGlobal( COME_Point3D pPos ){
00645 
00646         COME_Vector3D p_a = pPos - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00647         double distPointPlan = normalGlobalPosition.vpDotProduct( p_a );
00648         double signal = 1.0;
00649         if( distPointPlan < 0 ){
00650                 signal = -1;
00651         }
00652 
00653         COME_Point3D pPosPr = pPos - ( normalGlobalPosition * distPointPlan );
00654         
00655         distPointPlan = fabs( distPointPlan );
00656         
00657         if( !isInside2DGlobal( pPosPr ) ){
00658                 
00659                 COME_Point3D pA, pB, pC;
00660                 pA = COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00661                 pB = COME_Vector3D( getVertexGlobalPositionPt( 1 )->x, getVertexGlobalPositionPt( 1 )->y, getVertexGlobalPositionPt( 1 )->z );
00662                 pC = COME_Vector3D( getVertexGlobalPositionPt( 2 )->x, getVertexGlobalPositionPt( 2 )->y, getVertexGlobalPositionPt( 2 )->z );
00663                 
00664                 COME_Vector3D a_b = pB  - pA;
00665                 a_b.vpNormalize();
00666                 COME_Vector3D b_c = pC  - pB;
00667                 b_c.vpNormalize();
00668                 COME_Vector3D c_a = pA  - pC;
00669                 c_a.vpNormalize();
00670                 
00672                 COME_Point3D pOnEdgeAB = pA + ( a_b * ( (COME_Vector3D( pPosPr - pA )).vpDotProduct(a_b) )  );
00673                 COME_Point3D pOnEdgeBC = pB + ( b_c * ( (COME_Vector3D( pPosPr - pB )).vpDotProduct(b_c) )  );
00674                 COME_Point3D pOnEdgeCA = pC + ( c_a * ( (COME_Vector3D( pPosPr - pC )).vpDotProduct(c_a) )  );
00675                 
00677                 double EPSILON = 0.0000000001;
00678                 if( ( ( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() + EPSILON ) > (COME_Vector3D(pA-pB)).vpModule() )
00679                 ||  ( ( (COME_Vector3D(pB-pOnEdgeAB)).vpModule() + EPSILON ) > (COME_Vector3D(pA-pB)).vpModule() ) ){
00680                         if( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() < (COME_Vector3D(pB-pOnEdgeAB)).vpModule() )
00681                                 pOnEdgeAB = pA;
00682                         else pOnEdgeAB = pB;
00683                 }
00684                 if( ( ( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() + EPSILON ) > (COME_Vector3D(pB-pC)).vpModule() )
00685                 ||  ( ( (COME_Vector3D(pC-pOnEdgeBC)).vpModule() + EPSILON ) > (COME_Vector3D(pB-pC)).vpModule() ) ){
00686                         if( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() < (COME_Vector3D(pC-pOnEdgeBC)).vpModule() )
00687                                 pOnEdgeBC = pB;
00688                         else pOnEdgeBC = pC;
00689                 }
00690                 if( ( ( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() + EPSILON ) > (COME_Vector3D(pC-pA)).vpModule() )
00691                 ||  ( ( (COME_Vector3D(pA-pOnEdgeCA)).vpModule() + EPSILON ) > (COME_Vector3D(pC-pA)).vpModule() ) ){
00692                         if( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() < (COME_Vector3D(pA-pOnEdgeCA)).vpModule() )
00693                                 pOnEdgeCA = pC;
00694                         else pOnEdgeCA = pA;
00695                 }
00696         
00697                 double distProj = pPos.vpDistance( pOnEdgeAB );
00698                 if( distProj > pPos.vpDistance( pOnEdgeBC ) ){
00699                         distProj = pPos.vpDistance( pOnEdgeBC );
00700                 }
00701                 if( distProj > pPos.vpDistance( pOnEdgeCA ) ){
00702                         distProj = pPos.vpDistance( pOnEdgeCA );
00703                 }
00704                 return distProj * signal;
00705         }
00706         return distPointPlan * signal;
00707 }
00708 
00709 
00714 /*void
00715 COME_Face :: getClosestEdgeGlobal( COME_Point3D &pA, COME_Point3D &pB, COME_Point3D p ){
00716 
00717         pA = COME_Point3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00718         pB = COME_Point3D( getVertexGlobalPositionPt( 1 )->x, getVertexGlobalPositionPt( 1 )->y, getVertexGlobalPositionPt( 1 )->z );
00719         COME_Point3D pC = COME_Point3D( getVertexGlobalPositionPt( 2 )->x, getVertexGlobalPositionPt( 2 )->y, getVertexGlobalPositionPt( 2 )->z );
00720         
00721         if( p.vpDistance( pA ) < p.vpDistance( pB ) ){
00722                 if( p.vpDistance( pB ) < p.vpDistance( pC ) ){
00723                         return;
00724                 } else {
00725                         pB = pC;
00726                         return;
00727                 }
00728         } else if ( p.vpDistance( pB ) < p.vpDistance( pC ) ){
00729                 if( p.vpDistance( pC ) < p.vpDistance( pA ) ){
00730                         pA = pC;
00731                         return;
00732                 } else {
00733                         return;
00734                 }
00735         }       
00736         
00737 }*/
00738 
00743 bool
00744 COME_Face :: isInside2DGlobal( COME_Point3D pPos ){
00745 
00746         //find pA
00747         COME_Point3D pA = COME_Point3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00748         COME_Point3D pB = COME_Point3D( getVertexGlobalPositionPt( 1 )->x, getVertexGlobalPositionPt( 1 )->y, getVertexGlobalPositionPt( 1 )->z );
00749         COME_Point3D pC = COME_Point3D( getVertexGlobalPositionPt( 2 )->x, getVertexGlobalPositionPt( 2 )->y, getVertexGlobalPositionPt( 2 )->z );
00750 
00751         // build vectors
00752         COME_Vector3D c_b = pC  - pB;
00753         COME_Vector3D p_b = pPos  - pB;
00754         COME_Vector3D a_b = pA  - pB;
00755         COME_Vector3D c_a = pC  - pA;
00756         COME_Vector3D p_a = pPos  - pA;
00757         COME_Vector3D b_a = pB  - pA;
00758         
00759         if( ( c_b * p_b ).vpDotProduct( c_b * a_b ) < 0 ){
00760                 return false;
00761         } else if( ( c_a * p_a ).vpDotProduct( c_a * b_a ) < 0 ){
00762                 return false;
00763         } else if( ( b_a * p_a ).vpDotProduct( b_a * c_a ) < 0 ){
00764                 return false;
00765         } 
00766         return true;
00767 }
00768 
00772 double                          
00773 COME_Face :: distancePointPlaneGlobalPosition( COME_Point3D pPos ){
00774 
00775         COME_Vector3D p_a = pPos - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00776         double dist = normalGlobalPosition.vpDotProduct( p_a );
00777 
00778         return fabs(dist);
00779 
00780 }
00781 
00785 double                          
00786 COME_Face :: distanceVertexPlaneGlobalPosition( COME_Vertex3D vertex ){
00787 
00788         COME_Vector3D p_a = COME_Vector3D( vertex.x, vertex.y, vertex.z ) - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00789         double dist = normalGlobalPosition.vpDotProduct( p_a );
00790 
00791         return fabs(dist);
00792 
00793 }
00794 
00795 // Returns the vector having the given vertex as origin and its projection
00796 // on the face as end point.
00797 /*COME_Vector3D
00798 COME_Face :: vectorVertexFaceGlobalPosition( COME_Vertex3D vertex ){
00799 
00800         //COME_Vector3D direction = (vertex - getCenterGlobalPosition());       
00801         //return (normalGlobalPosition.vpNormalize().vpDotProduct(direction));
00802         
00803         vector<COME_Vertex3D*> verticesFace = getVerticesGlobalPositionPt(); 
00804 
00805         COME_Vector3D unitAB;
00806         unitAB.setVector3D( *(verticesFace[1]) - *(verticesFace[0]) );
00807         unitAB.vpNormalize();
00808 
00809         COME_Vector3D AP;
00810         AP.setVector3D( vertex - *(verticesFace[0]) );
00811 
00812         COME_Vector3D AB;
00813         AB.setVector3D( *(verticesFace[1]) - *(verticesFace[0]) );
00814 
00815         COME_Vector3D w1 =  unitAB * AP.vpModule() * AP.vpDotProduct( AB );                  
00816 
00817         COME_Vector3D unitAC;
00818         unitAC.setVector3D( *(verticesFace[2]) - *(verticesFace[0]) );
00819         unitAC.vpNormalize();
00820         
00821         COME_Vector3D AC;
00822         AC.setVector3D( *(verticesFace[2]) - *(verticesFace[0]) );
00823                              
00824         COME_Vector3D w2 =  unitAC * AP.vpModule() * AP.vpDotProduct( AC );     
00825 
00826         COME_Point3D pointPlane = w1 + w2 + vertex;
00827 
00828         
00829         return COME_Vector3D( pointPlane - vertex );
00830 
00831 }*/
00832 
00838 double                          
00839 COME_Face :: distanceFaceFace( COME_Face face ){
00840 
00841         vector<double> vDistance;
00842         double distance;
00843 
00844         vector<COME_Vertex3D> realVertices = getVertices();
00845         for ( int i = 0; i < realVertices.size(); i++ ){
00846         
00847                 if ( face.isInside( realVertices[ i ] )) {
00848                 
00849                         vDistance.push_back( face.distancePointFace( realVertices[ i ] ));
00850                 
00851                 }       
00852         }
00853 
00854         if ( vDistance.size() != 0 ){
00855         
00856                 distance = 0;
00857                 for ( int j = 0; j < vDistance.size(); j++ ){
00858 
00859                         if ( vDistance[ j ] > distance ){
00860 
00861                                 distance = vDistance[ j ];      
00862                         }
00863                 } 
00864         }
00865         else {
00866                 distance = 0;
00867         }
00868         return distance;
00869 }
00870 
00871 double
00872 COME_Face :: distanceFaceFaceGlobalPosition( COME_Face* faceN ){
00873 
00874         vector<double> vDistance;
00875         double distance;
00876         
00877         vector<COME_Vertex3D*> vertices = getVerticesGlobalPositionPt();
00878         
00879         //Distance is negative if vertex inside the mesh of the face!!
00880         for ( int i = 0; i < vertices.size(); i++ ){
00881                 
00882                 vDistance.push_back( faceN->distanceVertexPlaneGlobalPosition( *(vertices[ i ]) ));
00883         }
00884         
00885         distance = 0;
00886         for ( int j = 0; j < vDistance.size(); j++ ){
00887                 
00888                 if ( fabs( vDistance[ j ] ) > fabs( distance ) ){
00889                         distance = vDistance[ j ];
00890                         }
00891         } 
00892         
00893         //printf ("DIST=%f\n", distance);
00894         return distance;
00895 }
00896 
00897 
00904 bool
00905 COME_Face :: isInsideGlobalPosition( COME_Point3D currentPoint ){
00906         
00907         COME_Vector3D direction;
00908         
00909         direction.setVector3D( currentPoint - getCenterGlobalPosition());
00910         direction.vpNormalize();
00911         
00912         COME_Vector3D globalNormal =  getNormalGlobalPosition();
00913                 
00914         if( globalNormal.vpDotProduct( direction ) >= 0 ) return false;
00915         
00916         vector <int> neighbourFaces = getNeighbourFaces();
00917         for( int j = 0; j < neighbourFaces.size(); j++ ){
00918         
00919                 direction.setVector3D( currentPoint - ((COME_Mesh*)(parent))->getAFace(neighbourFaces[j]).getCenterGlobalPosition());
00920                 direction.vpNormalize();
00921                 
00922                 globalNormal =  ((COME_Mesh*)(parent))->getAFace(neighbourFaces[j]).getNormalGlobalPosition();
00923                 
00924                 if( globalNormal.vpDotProduct( direction ) >= 0 ) return false;
00925         }
00926         
00927         return true;
00928 }
00929 
00937 bool
00938 COME_Face::isInsideGlobalGOOD( COME_Point3D pPos ){
00939 
00940         COME_Vector3D direction;
00941         direction.setVector3D( pPos - getVertexGlobalPositionPt( 0 )->getVertexAsPoint3D() );
00942         direction.vpNormalize();
00943         if( normalGlobalPosition.vpDotProduct( direction ) >= 0 ) return false;
00944         
00945         COME_Vector3D p_a = pPos - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
00946         double distPointPlan = normalGlobalPosition.vpDotProduct( p_a );
00947 
00948         COME_Point3D pPosPr = pPos - ( normalGlobalPosition * distPointPlan );
00949         
00950         if( !isInside2DGlobal( pPosPr ) ) return false;
00951 
00952         return true;
00953 }
00954 
00955 
00961 bool
00962 COME_Face :: isInside( COME_Point3D currentPoint ){
00963 
00964         COME_Vector3D direction;
00965         vector<COME_Vertex3D> vert = getVertices();
00966         direction.setVector3D( currentPoint - vert[ 0 ] );
00967 
00968         direction.vpNormalize();
00969         if( getNormal().vpDotProduct( direction ) < 0 ){
00970 
00971                 return true; //inside
00972         } else {
00973                 
00974                 return false;
00975         }
00976 
00977 }
00978 
00979 
00983 bool
00984 COME_Face::containsVertexIndex( int v ){
00985 
00986         for( int i = 0; i < vertices.size(); i++ ){
00987 
00988                 if( vertices[i] == v ){
00989 
00990                         return true;
00991                 }
00992         }
00993         return false;
00994 }
00995 
01001 void
01002 COME_Face::replaceIndex( int v1, int v2 ){
01003 
01004         for( int i = 0; i < vertices.size(); i++ ){
01005 
01006                 if( vertices[i] == v1 ){
01007 
01008                         vertices[i] = v2;
01009                 }
01010         }
01011 }
01012 
01017 void
01018 COME_Face::update(){
01019 
01020         vector<COME_Vertex3D> vert = getVertices();
01021         COME_Vector3D displacement = nearestMolecule->getLastDisplacement();
01022         for ( int i = 0; i < vert.size(); i++ ){
01023 
01024                 vert[i] += displacement;
01025                 ((COME_Mesh*)parent)->setAVertex( vert[ i ], vertices[ i ] );
01026         }
01027 
01028 }
01029 
01030 void
01031 COME_Face:: setNeighbourFaces( vector<int> faces){
01032         neighbourFaces = faces;
01033 }
01034 
01035 vector <int>
01036 COME_Face:: getNeighbourFaces() const {
01037         return neighbourFaces;
01038 }
01039 
01040 void
01041 COME_Face:: addNeighbour( int faceIndex){
01042         neighbourFaces.push_back(faceIndex);
01043 }
01044 
01045 
01046 bool
01047 COME_Face::isCompletelyInside(COME_Face* face2){
01048         
01049         vector<COME_Vertex3D*> vertices = face2->getVerticesGlobalPositionPt();
01050         int in=0;
01051         for (int i = 0; i < vertices.size(); i++ ){
01052                 //cout << "Vertex to point (x,y,z) = " << vertices[i].getPoint3D().x << "," << vertices[i].getPoint3D().y << "," << vertices[i].getPoint3D().z << "\n";
01053                 if (!isInsideGlobalPosition(vertices[i]->getVertexAsPoint3D())){
01054                         return false;
01055                 }               
01056         }
01057         return true;
01058 }
01059 
01060 
01064 bool
01065 COME_Face::intersectionRayTriangle( COME_Vector3D vRay, COME_Point3D &intersection ){
01066 
01067         COME_Vector3D p_a = COME_Point3D() - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
01068         double distPointPlan = normalGlobalPosition.vpDotProduct( p_a );
01069 
01070         vRay.vpNormalize();
01071         
01072         double cosinus = vRay.vpDotProduct( -normalGlobalPosition );
01073         if( cosinus <= 0.0 ){
01074                 intersection = COME_Point3D();
01075                 return false;
01076         }
01077         
01078         double distOPlanRay = distPointPlan / cosinus;
01079         
01080         intersection = vRay * distOPlanRay;
01081         
01082         if( !isInside2DGlobal( intersection ) ){
01083                 intersection = COME_Point3D();
01084                 return false;
01085         }
01086         
01087         return true;
01088 }
01089 
01090 
01091 COME_Point3D
01092 COME_Face::getClosestPointOnFaceToPoint( COME_Point3D pPos ){
01093 
01094         COME_Vector3D p_a = pPos - COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
01095         double distPointPlan = normalGlobalPosition.vpDotProduct( p_a );
01096 
01097         COME_Point3D pPosPr = pPos - ( normalGlobalPosition * distPointPlan );
01098         
01099         if( !isInside2DGlobal( pPosPr ) ){
01100                 
01101                 COME_Point3D pA, pB, pC;
01102                 pA = COME_Vector3D( getVertexGlobalPositionPt( 0 )->x, getVertexGlobalPositionPt( 0 )->y, getVertexGlobalPositionPt( 0 )->z );
01103                 pB = COME_Vector3D( getVertexGlobalPositionPt( 1 )->x, getVertexGlobalPositionPt( 1 )->y, getVertexGlobalPositionPt( 1 )->z );
01104                 pC = COME_Vector3D( getVertexGlobalPositionPt( 2 )->x, getVertexGlobalPositionPt( 2 )->y, getVertexGlobalPositionPt( 2 )->z );
01105                 
01106                 COME_Vector3D a_b = pB  - pA;
01107                 a_b.vpNormalize();
01108                 COME_Vector3D b_c = pC  - pB;
01109                 b_c.vpNormalize();
01110                 COME_Vector3D c_a = pA  - pC;
01111                 c_a.vpNormalize();
01112                 
01114                 COME_Point3D pOnEdgeAB = pA + ( a_b * ( (COME_Vector3D( pPosPr - pA )).vpDotProduct(a_b) )  );
01115                 COME_Point3D pOnEdgeBC = pB + ( b_c * ( (COME_Vector3D( pPosPr - pB )).vpDotProduct(b_c) )  );
01116                 COME_Point3D pOnEdgeCA = pC + ( c_a * ( (COME_Vector3D( pPosPr - pC )).vpDotProduct(c_a) )  );
01117                 
01119                 double EPSILON = 0.0000000001;
01120                 if( ( ( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() + EPSILON ) > (COME_Vector3D(pA-pB)).vpModule() )
01121                 ||  ( ( (COME_Vector3D(pB-pOnEdgeAB)).vpModule() + EPSILON ) > (COME_Vector3D(pA-pB)).vpModule() ) ){
01122                         if( ( COME_Vector3D(pA-pOnEdgeAB)).vpModule() < (COME_Vector3D(pB-pOnEdgeAB)).vpModule() )
01123                                 pOnEdgeAB = pA;
01124                         else pOnEdgeAB = pB;
01125                 }
01126                 if( ( ( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() + EPSILON ) > (COME_Vector3D(pB-pC)).vpModule() )
01127                 ||  ( ( (COME_Vector3D(pC-pOnEdgeBC)).vpModule() + EPSILON ) > (COME_Vector3D(pB-pC)).vpModule() ) ){
01128                         if( ( COME_Vector3D(pB-pOnEdgeBC)).vpModule() < (COME_Vector3D(pC-pOnEdgeBC)).vpModule() )
01129                                 pOnEdgeBC = pB;
01130                         else pOnEdgeBC = pC;
01131                 }
01132                 if( ( ( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() + EPSILON ) > (COME_Vector3D(pC-pA)).vpModule() )
01133                 ||  ( ( (COME_Vector3D(pA-pOnEdgeCA)).vpModule() + EPSILON ) > (COME_Vector3D(pC-pA)).vpModule() ) ){
01134                         if( ( COME_Vector3D(pC-pOnEdgeCA)).vpModule() < (COME_Vector3D(pA-pOnEdgeCA)).vpModule() )
01135                                 pOnEdgeCA = pC;
01136                         else pOnEdgeCA = pA;
01137                 }
01138                 
01139                 double distProj = pPos.vpDistance( pOnEdgeAB );
01140                 COME_Point3D pPosOnEdge = pOnEdgeAB;
01141                 if( distProj > pPos.vpDistance( pOnEdgeBC ) ){
01142                         distProj = pPos.vpDistance( pOnEdgeBC );
01143                         pPosOnEdge = pOnEdgeBC;
01144                 }
01145                 if( distProj > pPos.vpDistance( pOnEdgeCA ) ){
01146                         distProj = pPos.vpDistance( pOnEdgeCA );
01147                         pPosOnEdge = pOnEdgeCA;
01148                 }
01149                 return pPosOnEdge;
01150         }
01151         
01152         return pPosPr;
01153 }
01154 
01155 COME_Face *
01156 COME_Face::getNeighborClosestTo( COME_Vertex3D *vPos ){
01157 
01158         COME_Face * fClosest = this;
01159         double distMin = this->distancePointFaceGlobal(  vPos->getVertexAsPoint3D() );
01160         
01161         for( int iNF = 0; iNF < neighbourFaces.size(); iNF++ ){
01162         
01163                 if( ((COME_Mesh*)(parent))->getAFacePt( neighbourFaces[iNF] )->distancePointFaceGlobal( vPos->getVertexAsPoint3D() ) < distMin ){
01164                         distMin = ((COME_Mesh*)(parent))->getAFacePt( neighbourFaces[iNF] )->distancePointFaceGlobal( vPos->getVertexAsPoint3D() );
01165                         fClosest = ((COME_Mesh*)(parent))->getAFacePt( neighbourFaces[iNF] );
01166                 }
01167                 //printf( "%p dist Neighbor: %f min: %f index neighbor:%d \n", this, ((COME_Mesh*)(parent))->getAFacePt( neighbourFaces[iNF] )->distancePointFaceGlobal( vPos->getVertexAsPoint3D() ), distMin, iNF );
01168         }
01169         return fClosest;
01170 }
01171 
01172 COME_Vector3D
01173 COME_Face::getBlendedVelocity( COME_Point3D pPos ){
01174 
01175         COME_Vector3D velocBlended;
01176         
01177         COME_Vertex3D *A =  getVertexGlobalPositionPt( 0 );
01178         COME_Vertex3D *B =  getVertexGlobalPositionPt( 1 );
01179         COME_Vertex3D *C =  getVertexGlobalPositionPt( 2 );
01180         
01181         double m1 = COME_Face::triangleArea( B->getVertexAsPoint3D(), C->getVertexAsPoint3D(), pPos );
01182         double m2 = COME_Face::triangleArea( C->getVertexAsPoint3D(), A->getVertexAsPoint3D(), pPos );
01183         double m3 = COME_Face::triangleArea( A->getVertexAsPoint3D(), B->getVertexAsPoint3D(), pPos );
01184         
01185         double sumMasses = m1 + m2 + m3;
01186         
01187         double w1 = m1 / sumMasses;
01188         double w2 = m2 / sumMasses;
01189         double w3 = m3 / sumMasses;
01190         
01191         velocBlended = ( A->getBlendedVelocity() * w1) + ( B->getBlendedVelocity() * w2) + ( C->getBlendedVelocity() * w3);
01192         
01193         return velocBlended;
01194 }
01195 
01196         /* // old and imprecise version
01197         vector<COME_Vertex3D*> verticesFace = getVerticesGlobalPositionPt(); 
01198 
01199         COME_Vector3D unitAB;
01200         unitAB.setVector3D( *(verticesFace[1]) - *(verticesFace[0]) );
01201         unitAB.vpNormalize();
01202         
01203         COME_Vector3D unitAC;
01204         unitAC.setVector3D( *(verticesFace[2]) - *(verticesFace[0]) );
01205         unitAC.vpNormalize();
01206 
01207         COME_Vector3D AP;
01208         AP.setVector3D( pPos - *(verticesFace[0]) );
01209 
01210         COME_Vector3D w1 =  unitAB * AP.vpDotProduct( unitAB );
01211         COME_Vector3D w2 =  (normalGlobalPosition * w1) * AP.vpDotProduct( unitAC );
01212 
01213         COME_Point3D pointPlane = *(verticesFace[0]) + w1 + w2;
01214         
01215         return pointPlane;*/
comeface.cpp
