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  **************************************************************************/
00020 
00022 //
00023 //  PROJECT.....: CO-ME
00024 //  RESPONSIBLE.: 
00025 //
00026 //  AUTHOR......: Anderson Maciel
00027 //  DATE........: August/05/2002
00028 //  DESCRIPTION.: Class definition.
00029 //
00031 
00032 
00033 #include <bio/comemoleculestissue.h>
00034 #include <bio/comemoleculescartilage.h>
00035 #include <bio/comebiostructure.h>
00036 #include <bio/comematerial.h>
00037 #include <math.h>
00038 #include <algebra/comemesh.h>
00039 #include <algebra/comevertex3d.h>
00040 
00041 //#include <LinAlg/linalg.h>
00042 
00043 //#include <conio.h>
00044 
00045 //using namespace LinAlg;
00046 
00047 //#include <MarchingCubes/marching_cubes.h>
00048 //#include <MarchingCubes/blob.h>
00049 
00053 
00054 COME_MoleculesTissue::COME_MoleculesTissue()
00055 {
00056         radiusAverage = 0.0;
00057         mass = 0.0;
00058         accu = 0.0;
00059         averageNumConnectionsPerMolecule = 0;
00060         dymv = NULL;
00061         dyma = NULL;
00062         dytv = NULL;
00063         dyta = NULL;
00064         ytp = NULL;
00065         ytv = NULL;
00066         gpuForces = NULL;
00067         parent = NULL;
00068 }
00069 
00070 COME_MoleculesTissue::~COME_MoleculesTissue()
00071 {
00072         if( dymv ) delete [] dymv;
00073         if( dyma ) delete [] dyma;
00074         if( dytv ) delete [] dytv;
00075         if( dyta ) delete [] dyta;
00076         if( ytp ) delete [] ytp;
00077         if( ytv ) delete [] ytv;
00078         if( gpuForces ) delete [] gpuForces;
00079         
00080         list<COME_Molecule*>::iterator iter;
00081         for (iter = shape.begin(); iter != shape.end(); iter++){
00082                 delete (*iter);
00083         }
00084 }
00085 
00086 void
00087 COME_MoleculesTissue::initForSimulation(){
00088 
00089         dymv = new COME_Vector3D[shape.size()];
00090         dyma = new COME_Vector3D[shape.size()];
00091         
00092         dytv = new COME_Vector3D[shape.size()];
00093         dyta = new COME_Vector3D[shape.size()];
00094         
00095         ytp = new COME_Vector3D[shape.size()];
00096         ytv = new COME_Vector3D[shape.size()];
00097         
00098         gpuForces = new COME_Vector3D[shape.size()];
00099 }
00100 
00101 list<COME_Molecule*>*
00102 COME_MoleculesTissue::getShape(){
00103 
00104         return &shape;
00105 }
00106 
00107 void
00108 COME_MoleculesTissue::setRadiusAverage( double radiusN ){
00109 
00110         radiusAverage = radiusN;
00111 }
00112 
00113 double
00114 COME_MoleculesTissue::getRadiusAverage(){
00115 
00116         return radiusAverage;
00117 }
00118 
00119 int
00120 COME_MoleculesTissue::getAverageNumConnections(){
00121 
00122         return averageNumConnectionsPerMolecule;
00123 }
00124 
00129 bool
00130 COME_MoleculesTissue::isTooClose( COME_Point3D currentPoint, double distance ){
00131 
00132         list<COME_Molecule*>::iterator iter;
00133         for (iter = shape.begin(); iter != shape.end(); iter++){
00134                 
00135                 if( (*iter)->getPosition().vpDistance( currentPoint ) < distance ){
00136                         return true;
00137                 } 
00138         }
00139         return false;
00140 }
00141 
00147 void
00148 COME_MoleculesTissue::setMoleculesMaterial( COME_Material* newMaterial ){
00149 
00150         list<COME_Molecule*>::iterator iter;
00151         for (iter = shape.begin(); iter != shape.end(); iter++){
00152                 
00153                 (*iter)->setMaterial( newMaterial );
00154         }
00155 }
00156 
00157 void
00158 COME_MoleculesTissue::setMoleculesDensity( double densityN ){
00159 
00160         list<COME_Molecule*>::iterator iter;
00161         for (iter = shape.begin(); iter != shape.end(); iter++){
00162                 
00163                 (*iter)->getMaterial()->setDensity( densityN );
00164         }
00165 }
00166 
00167 void
00168 COME_MoleculesTissue::setMoleculesElasticity( double elasticityN ){
00169 
00170         list<COME_Molecule*>::iterator iter;
00171         for (iter = shape.begin(); iter != shape.end(); iter++){
00172                 
00173                 (*iter)->getMaterial()->setYoungsModulus( elasticityN );
00174         }
00175 }
00176 
00177 void
00178 COME_MoleculesTissue::setMoleculesFrictionConst( double frictionN ){
00179 
00180         list<COME_Molecule*>::iterator iter;
00181         for (iter = shape.begin(); iter != shape.end(); iter++){
00182                 
00183                 (*iter)->setFrictionConst( frictionN );
00184         }
00185 }
00186 
00189 COME_Material*
00190 COME_MoleculesTissue::getMoleculesMaterial(){
00191 
00192         return shape.front()->getMaterial();
00193 }
00194 
00195 double
00196 COME_MoleculesTissue::getDensity(){
00197 
00198         return shape.front()->getMaterial()->getDensity();
00199 }
00200 
00201 double
00202 COME_MoleculesTissue::getElasticity(){
00203 
00204         return shape.front()->getMaterial()->getYoungsModulus();
00205 }
00206 
00207 double
00208 COME_MoleculesTissue::getFrictionConst(){
00209 
00210         return shape.front()->getFrictionConst();
00211 }
00212 
00213 double
00214 COME_MoleculesTissue::getTopStress(){
00215 
00216         double top = 0.0;
00217         list<COME_Molecule*>::iterator iter;
00218         for (iter = shape.begin(); iter != shape.end(); iter++ ){
00219                 
00220                 if( top < (*iter)->getStress() )
00221                         top=(*iter)->getStress();
00222         }
00223         return top;
00224 }
00225 
00226 COME_Force
00227 COME_MoleculesTissue::getGPUForce( int i ){
00228 
00229         return gpuForces[i];
00230 
00231 }
00233 
00234 
00239 void
00240 COME_MoleculesTissue::addMolecule( COME_Molecule *moleculeN ){
00241         
00242         // test for all existent molecules if distance to moleculeN is smaller than its fracture distance
00243         bool canAdd = true;
00244         list<COME_Molecule*>::iterator iter;
00245         for (iter = shape.begin(); iter != shape.end(); iter++){
00246                 // if it is, create a connection between them
00247                 if( ((COME_Point3D)((*iter)->getPosition())).vpDistance( moleculeN->getPosition() ) < ( (*iter)->getRadius() + moleculeN->getRadius() ) ){
00248                         
00249                         if( ((COME_BioStructure*)parent)->getGlobalHooke() == 0.0 ){
00250                                 moleculeN->createLink( (*iter) );
00251                         } else {
00252                                 moleculeN->createLink( (*iter), ((COME_BioStructure*)parent)->getGlobalHooke() );
00253                         }
00254                         canAdd = true;
00255                 }
00256                 
00257         }
00258 
00259         if( canAdd ){
00260                 radiusAverage = ( radiusAverage*shape.size() + moleculeN->getRadius() ) / ( shape.size() + 1 );
00261                 moleculeN->setIndex( shape.size() );
00262                 shape.push_back( moleculeN );
00263                 mass += moleculeN->getMass();
00264         }
00265 }
00266 
00272 void
00273 COME_MoleculesTissue::addMolecule( COME_Molecule *moleculeN, double minDist ){
00274         
00275         
00276         list<COME_Molecule*>::iterator iter;
00277 
00278         // test if this new molecule is too close of any other in the list (test minDist)
00279         for (iter = shape.begin(); iter != shape.end(); iter++){
00280                 // if it is, don't add it
00281                 if( ((COME_Point3D)((*iter)->getPosition())).vpDistance( moleculeN->getPosition() ) >= ( minDist ) ){
00282                         // do nothing
00283                 } else {
00284                         return;
00285                 }               
00286         }
00287 
00288         // test for all existent molecules if distance to moleculeN is smaller than its fracture distance
00289         for (iter = shape.begin(); iter != shape.end(); iter++){
00290                 // if it is, create a connection between them
00291                 if( ((COME_Point3D)((*iter)->getPosition())).vpDistance( moleculeN->getPosition() ) < ( (*iter)->getRadius() + moleculeN->getRadius() ) ) {
00292                         moleculeN->createLink( (*iter) );
00293                 }
00294                 
00295         }
00296 
00297         radiusAverage = ( radiusAverage*shape.size() + moleculeN->getRadius() ) / ( shape.size() + 1 );
00298         moleculeN->setIndex( shape.size() );
00299         shape.push_back( moleculeN );
00300         mass += moleculeN->getMass();
00301 }
00302 
00303 bool
00304 COME_MoleculesTissue::updateEuler( double timestep, double simClock, string fileNameToSaveDeformation ){
00305 
00306         //printf( "ENTROU\n");
00307         simTime = simClock;
00308         if( !dymv ){
00309                 initForSimulation();
00310         }
00311         int i = 0;
00312         list<COME_Molecule*>::iterator iter;
00313         for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00314         
00315                 ytp[i] = (*iter)->getGlobalPosition();
00316                 ytv[i] = (*iter)->getVelocity();
00317                 //printf( "pos: %f vel: %f \n", (*iter)->getGlobalPosition().getY(), (*iter)->getVelocity().getY() );
00318         }
00319         
00320         i = 0;
00321         for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00322                 
00323                 if( !(*iter)->isFixed() ){
00325                         (*iter)->derivs( simClock, timestep, ytp, ytv, dytv, dyta, i ); 
00326                         
00328                         (*iter)->setPositionMinusOne( (*iter)->getGlobalPosition() );
00329                         COME_Point3D newpos( (*iter)->getGlobalPosition() + dytv[i]*timestep );
00330                         (*iter)->setPosition( newpos );
00331                         (*iter)->setVelocity( dytv[i] );
00332                 }
00333         }
00334         
00336          if( ( COME::flagCollisionTreatment == NEIGHBORS ) || ( COME::flagCollisionTreatment == SPHERICAL_SLIDING ) ){
00337                 i = 0;
00338                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00339                         
00340                         if( !(*iter)->isFixed() ){
00341                                 
00342                                 vector<COME_Vertex3D*> *vertices = (*iter)->getBuois();
00343                                 int iv;
00344                                 
00345                                 COME_Vector3D dispPosition;
00346                                 COME_Vector3D dispVelocity;
00347                                 for( iv = 0; iv < vertices->size(); iv++ ){
00348                                         
00349                                         if( dispPosition.vpModule() < ( (*vertices)[iv]->getCollisionDispAvg().vpModule()  ) )
00350                                                 dispPosition = ( (*vertices)[iv]->getCollisionDispAvg()  );
00351                                         if( dispVelocity.vpModule() < ( (*vertices)[iv]->getVelocityDispAvg().vpModule() ) )
00352                                                 dispVelocity = ( (*vertices)[iv]->getVelocityDispAvg() );
00353                                 }
00354                                 COME_Vector3D newVeloc = dytv[i] + dispVelocity;
00355                                 dytv[i] = newVeloc;
00356                                 
00357                                 COME_Point3D newpos( (*iter)->getGlobalPosition() + dispPosition );
00358                                 
00360                                 (*iter)->setPosition( newpos );
00361                                 (*iter)->setVelocity( dytv[i] );
00362                         }
00363                 }
00364         }
00365         
00367         accu += timestep;
00368         if( COME::flagExportToPrecalculatedFile ){
00369                 if( ( fileNameToSaveDeformation != "" ) && ( accu >= COME::flagAnimationResolution ) ){
00370                         //printf( "CALL accu: %f ts: \n", accu, timestep );
00371                         saveDeformationFile( UPDATE );
00372                         accu = 0.0;
00373                 }
00374         }
00375                 
00376         makeAllLocalFrames();
00377 
00378         return true;
00379 }
00380 
00387 bool
00388 COME_MoleculesTissue::update( double timestep, double simClock, string fileNameToSaveDeformation ){
00389 
00390         int i;
00391         double xh , hh, h6;
00392         list<COME_Molecule*>::iterator iter;
00393         
00394         simTime = simClock;
00395         
00396         //if( deformation.size() <= 0 ){
00397 
00398                 if( !dymv ){
00399                         initForSimulation();
00400                 }
00401                 
00402                 #ifdef GPU_FORCES
00403                 // Call CG function to compute current forces on each molecule
00404                 // in this case, derivs never calls calculateForce
00405                 gpuForces = NULL; //TODAS AS FORÇAS
00406                 #endif
00407 
00408                 hh = timestep * 0.5;
00409                 h6 = timestep / 6.0;
00410                 xh = simClock + hh;
00411                 
00412                 // First step //////////////////////////////////////////////////////
00413                 i = 0;
00414                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00415                         ytp[i] = ( (*iter)->getVelocity() * hh ) + (*iter)->getGlobalPosition();
00416                         ytv[i] = ( (*iter)->getAcceleration() * hh ) + (*iter)->getVelocity();
00417                 }
00418 
00419                 // Second step //////////////////////////////////////////////////////
00420                 i = 0;
00421                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00422                         (*iter)->derivs( xh, hh, ytp, ytv, dytv, dyta, i ); 
00423                 }
00424 
00425                 i = 0;
00426                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00427                         ytp[i] = ( dytv[i] * hh ) + (*iter)->getGlobalPosition();
00428                         ytv[i] = ( dyta[i] * hh ) + (*iter)->getVelocity();
00429                 }
00430 
00431                 // Third step //////////////////////////////////////////////////////
00432                 i = 0;
00433                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00434                         (*iter)->derivs( xh, hh, ytp, ytv, dymv, dyma, i ); 
00435                 }
00436 
00437                 i = 0;
00438                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00439                         ytp[i] = ( dymv[i] * timestep ) + (*iter)->getGlobalPosition();
00440                         ytv[i] = ( dyma[i] * timestep ) + (*iter)->getVelocity();
00441                         
00442                         dymv[i] += dytv[i];
00443                         dyma[i] += dyta[i];
00444                 }
00445 
00446                 // Fourth step //////////////////////////////////////////////////////
00447                 i = 0;
00448                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00449                         (*iter)->derivs( simClock + timestep, timestep, ytp, ytv, dytv, dyta, i ); 
00450                 }
00451 
00452                 i = 0;
00453                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00454                         
00455                         ytp[i] = ( ( ( (*iter)->getVelocity() + dytv[i] + ( dymv[i] * 2.0 ) ) * h6 ) + (*iter)->getGlobalPosition() );
00456                         ytv[i] = ( ( ( (*iter)->getAcceleration() + dyta[i] + ( dyma[i] * 2.0 ) ) * h6 ) + (*iter)->getVelocity() );
00457                         
00458                         if( ( timestep < 0 ) || ( ( ytv[i].vpModule() * timestep ) < ( (*iter)->getRadius() / 20.0 ) ) ){
00459                                 // do nothing
00460                         } else {
00461                                 //return false;
00462                         }
00463                 }
00464 
00465                 i = 0;
00466                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00467                         
00468                         if( !(*iter)->isFixed() ){
00469                                 (*iter)->setPositionMinusOne( (*iter)->getGlobalPosition() );
00470                                 (*iter)->setPosition( ytp[i] );
00471                                 //(*iter)->setAcceleration( ( ytv[i] - (*iter)->getVelocity() ) / timestep ); // a = ( V- Vo ) / deltaT (Don't need to do it because it's done in derivs.)
00472                                 (*iter)->setVelocity( ytv[i] );
00473                                 //(*iter)->makeLocalFrame();
00474                         }
00475                 }
00476                 
00478                 if( ( COME::flagCollisionTreatment == NEIGHBORS ) || ( COME::flagCollisionTreatment == SPHERICAL_SLIDING ) ){
00479                         i = 0;
00480                         for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00481                                 
00482                                 if( !(*iter)->isFixed() ){
00483                                         
00484                                         vector<COME_Vertex3D*> *vertices = (*iter)->getBuois();
00485                                         int iv;
00486                                         
00487                                         COME_Vector3D dispPosition;
00488                                         COME_Vector3D dispVelocity;
00489                                         for( iv = 0; iv < vertices->size(); iv++ ){
00490                                                 
00491                                                 if( dispPosition.vpModule() < ( (*vertices)[iv]->getCollisionDispAvg().vpModule()  ) )
00492                                                         dispPosition = ( (*vertices)[iv]->getCollisionDispAvg()  );
00493                                                 if( dispVelocity.vpModule() < ( (*vertices)[iv]->getVelocityDispAvg().vpModule() ) )
00494                                                         dispVelocity = ( (*vertices)[iv]->getVelocityDispAvg() );
00495                                         }
00496                                         COME_Vector3D newVeloc = ytv[i] + dispVelocity;
00497                                         ytv[i] = newVeloc;
00498                                         
00499                                         COME_Point3D newpos( (*iter)->getGlobalPosition() + dispPosition );
00500                                         
00502                                         (*iter)->setPosition( newpos );
00503                                         (*iter)->setVelocity( ytv[i] );
00504                                 }
00505                         }
00506                 }
00507                 
00508 
00510                 accu += timestep;
00511                 if( flagExportToPrecalculatedFile ){
00512                         if( ( fileNameToSaveDeformation != "" ) && ( accu >= COME::flagAnimationResolution ) ){
00513                                 saveDeformationFile( UPDATE );
00514                                 accu = 0.0;
00515                         }
00516                 }
00517 
00518         /*} else {
00520                 accu += timestep;
00521                 
00522 
00523                 unsigned int t = (unsigned int)(simClock / COME::flagAnimationResolution);
00524                 int i = 0;
00525                 for (iter = shape.begin(); iter != shape.end(); iter++, i++){
00526                         
00527                         (*iter)->setPositionMinusOne( (*iter)->getGlobalPosition() );
00528                         (*iter)->setPosition( deformation[t][i] );
00529                         (*iter)->setStress( deformationStresses[t][i] );
00530                         (*iter)->setStrain( deformationStrains[t][i] );
00531                         //(*iter)->makeLocalFrame();
00532                 }
00533                 if( accu >= COME::flagAnimationResolution ){
00534                         accu = 0.0;
00535                 }
00536 
00537         }*/
00538         makeAllLocalFrames();
00539         return true;
00540 
00541 }
00542 
00543 
00549 void
00550 COME_MoleculesTissue::getEnvelop( COME_Point3D& mins, COME_Point3D& maxs ){
00551 
00552         if( shape.size() == 0 ){
00553                 mins = COME_Point3D();
00554                 maxs = COME_Point3D();
00555                 return;
00556         }
00557 
00558         mins = maxs = shape.front()->getGlobalPosition();
00559 
00560         list<COME_Molecule*>::iterator iter;
00561         for (iter = shape.begin(); iter != shape.end(); iter++){
00562                 
00563                 COME_Point3D minsM = (*iter)->getGlobalPosition();
00564                 COME_Point3D maxsM = (*iter)->getGlobalPosition();
00565 
00566                 if( minsM.getX() < mins.getX() ){
00567                         mins.setX( minsM.getX() );
00568                 }
00569                 if( minsM.getY() < mins.getY() ){
00570                         mins.setY( minsM.getY() );
00571                 }
00572                 if( minsM.getZ() < mins.getZ() ){
00573                         mins.setZ( minsM.getZ() );
00574                 }
00575 
00576                 if( maxsM.getX() > maxs.getX() ){
00577                         maxs.setX( maxsM.getX() );
00578                 }
00579                 if( maxsM.getY() > maxs.getY() ){
00580                         maxs.setY( maxsM.getY() );
00581                 }
00582                 if( maxsM.getZ() > maxs.getZ() ){
00583                         maxs.setZ( maxsM.getZ() );
00584                 }
00585         }
00586 
00587         // if spheres size change this will not work perfectly
00588         mins -= shape.front()->getRadius();
00589         maxs += shape.back()->getRadius();
00590 
00591 }
00592 
00596 COME_Molecule*
00597 COME_MoleculesTissue::nearestTo( COME_Point3D point ){
00598 
00599         list<COME_Molecule*>::iterator nearest = shape.begin();
00600         double minDist = point.vpDistance( (*nearest)->getPosition() );
00601 
00602         list<COME_Molecule*>::iterator iter;    
00603         for (iter = shape.begin(); iter != shape.end(); iter++){
00604 
00605                 if( minDist > point.vpDistance( (*iter)->getPosition() ) ){
00606 
00607                         minDist = point.vpDistance( (*iter)->getPosition() );
00608                         nearest = iter;
00609                 }
00610         }
00611 
00612         return (*nearest);
00613 }
00614 
00618 double
00619 COME_MoleculesTissue::getMass(){
00620 
00621         return mass;
00622 }
00623 
00628 COME_Point3D
00629 COME_MoleculesTissue::getMassCenter(){
00630 
00631         double x, y, z, mass;
00632         x = y = z = mass = 0.0;
00633 
00634         list<COME_Molecule*>::iterator iter;    
00635         for (iter = shape.begin(); iter != shape.end(); iter++){
00636 
00637                 x += (*iter)->getMass() * (*iter)->getPosition().getX();
00638                 y += (*iter)->getMass() * (*iter)->getPosition().getY();
00639                 z += (*iter)->getMass() * (*iter)->getPosition().getZ();
00640 
00641                 mass += (*iter)->getMass();
00642         }
00643 
00644         COME_Point3D massCenter( x/mass, y/mass, z/mass );
00645         return massCenter;
00646 }
00647 
00652 void
00653 COME_MoleculesTissue::recalculateAllSpringConstants( int mode, double currTime ){
00654 
00655         updateAverageNumConnections();
00656         list<COME_Molecule*>::iterator iter;    
00657         for (iter = shape.begin(); iter != shape.end(); iter++){
00658 
00659                 double increment = 1.0;
00660                 if( mode == NUMERICALLY ){
00661                         if( chain.currentLength > 0.0 ){
00662                                 
00663                                 chain.updateLengths( currTime );
00664 
00665                                 chain.targetLength = chain.nominalLength + ( ( (chain.force*100) * chain.nominalLength ) / ( 100000.0 * 1.0 ) );
00666 
00667                                 if( chain.currentLength > chain.targetLength ){
00668 
00669                                         if( chain.currentLength > chain.lastLength ){
00670 
00671                                                 // reduce k
00672                                                 increment = 1.001;
00673 
00674                                         } else {
00675 
00676                                                 // do nothing
00677                                                 increment = 1.0;
00678 
00679                                         }
00680 
00681                                 } else if( chain.currentLength < chain.targetLength ){
00682 
00683                                         if( chain.currentLength > chain.lastLength ){
00684 
00685                                                 // do nothing
00686                                                 increment = 1.0;
00687 
00688                                         } else {
00689 
00690                                                 // increase k
00691                                                 increment = 0.999;
00692 
00693                                         }
00694                                 }
00695                         }
00696                 }
00697 
00698                 list<COME_MoleculeLink*>::iterator iterLink;
00699                 for (iterLink = (*iter)->getConnectionList()->begin(); iterLink != (*iter)->getConnectionList()->end(); iterLink++ ){
00700 
00701                         if( mode == NUMERICALLY ){
00702                                 (*iterLink)->recalculateSpring( NUMERICALLY, increment );
00703                         } else {
00704                                 (*iterLink)->recalculateSpring( mode, averageNumConnectionsPerMolecule );
00705                         }
00706                 }
00707         }
00708 }
00709 
00711 void
00712 COME_MoleculesTissue::setSurfaceMolecules(){
00713 
00714         list<COME_Molecule*>::iterator iter;    
00715         for (iter = shape.begin(); iter != shape.end(); iter++){
00716                 
00717                 if( (*iter)->getConnectionList()->size() <= 17 ){
00718                         (*iter)->setOnSurface( true );
00719                 }
00720         }
00721 }
00722 
00727 void
00728 COME_MoleculesTissue::updateAverageNumConnections(){
00729 
00730         unsigned int total = 0;
00731 
00732         list<COME_Molecule*>::iterator iter;    
00733         for (iter = shape.begin(); iter != shape.end(); iter++){
00734                 
00735                 total += (*iter)->getConnectionList()->size();
00736         }
00737 
00738         if( shape.size() > 0 )
00739                 averageNumConnectionsPerMolecule = total / shape.size();
00740         else
00741                 averageNumConnectionsPerMolecule = 0;
00742 }
00743 
00748 void
00749 COME_MoleculesTissue::loadDeformationFile( string fileName ){
00750 
00751         FILE *file;
00752 
00753         unsigned int numberOfMolecules = getShape()->size();
00754 
00755         if( file = fopen( fileName.c_str(), "rb" ) ){
00756 
00757                 while( !feof( file ) ){
00758 
00759                         COME_Point3D *localDeformation =  new COME_Point3D [numberOfMolecules];
00760                         double *localDeformationStresses = new double [numberOfMolecules];
00761                         double *localDeformationStrains = new double [numberOfMolecules];
00762                         
00763                         for( int j = 0; j < numberOfMolecules; j++ ){
00764 
00765                                 double array[5];
00766                                 fread( array, sizeof(double), 5, file );
00767 
00768                                 localDeformation[j].setXYZ( array[0], array[1], array[2] );
00769                                 localDeformationStresses[j] = array[3];
00770                                 localDeformationStrains[j] = array[4];
00771                         }
00772                         deformation.push_back( localDeformation );
00773                         deformationStresses.push_back( localDeformationStresses );
00774                         deformationStrains.push_back( localDeformationStrains );
00775                 }               
00776 
00777                 fclose(file);
00778         }
00779 
00780 }
00781 
00785 /* old code for molecules
00786         /*list<COME_Molecule*>::iterator iter;
00787         for (iter = shape.begin(); iter != shape.end(); iter++ ){
00788                 
00789                 double array[6];
00790                 COME_Point3D localPoint = (*iter)->getGlobalPosition(); 
00791                 array[0] = simTime;
00792                 array[1] = localPoint.getX();
00793                 array[2] = localPoint.getY();
00794                 array[3] = localPoint.getZ();
00795                 array[4] = (*iter)->getStress();
00796                 array[5] = (*iter)->getStrain();
00797                 fwrite( array, sizeof(double), 6, file );
00798         }
00799 */
00801 /*
00802 void
00803 COME_MoleculesTissue::saveDeformationFile( int fileMode ){
00804         
00805         FILE *file;
00806         string strMode;
00807 
00808         if( fileMode == NEW ){
00809                 strMode = "wb";// create new
00810                 fclose( fopen( ((COME_MoleculesCartilage*)parent)->getDeformationFile().c_str(), strMode.c_str() ) );
00811         } else {
00812                 strMode = "a+b"; // update at the end
00813         
00814                 if( file = fopen( ((COME_MoleculesCartilage*)parent)->getDeformationFile().c_str(), strMode.c_str() ) ){
00815         
00816                         vector<COME_Vertex3D> *verts = ((COME_MoleculesCartilage*)parent)->getSurface()->getVerticesGlobalPt();
00817                         //vector<COME_Vertex3D> *vertsOriginal = ((COME_MoleculesCartilage*)parent)->getSurface()->getVerticesPt();
00818                         for ( int i = 0; i < verts->size(); i++ ){
00819 
00820                                 COME_Vertex3D vertex = verts->at(i);
00821                                 double array[6];
00822                                 array[0] = simTime;
00823                                 array[1] = vertex.x;
00824                                 array[2] = vertex.y;
00825                                 array[3] = vertex.z;
00826                                 array[4] = vertex.getStress();
00827                                 array[5] = vertex.getStrain();
00828                                 fwrite( array, sizeof(double), 6, file );
00829                         }
00830                         fclose(file);
00831                 }
00832         }
00833 }
00834 */
00835 
00836 void
00837 COME_MoleculesTissue::saveDeformationFile( int fileMode ){
00838         
00839         FILE *file;
00840         string strMode;
00841 
00842         if( fileMode == NEW ){
00843                 strMode = "wb";// create new
00844                 
00845                 fclose( fopen( ((COME_MoleculesCartilage*)parent)->getDeformationFile().c_str(), strMode.c_str() ) );
00846                 
00847                 string f_namei = ((COME_MoleculesCartilage*)parent)->getDeformationFile();
00848                 string tempStri( f_namei, 0, f_namei.size()-10 );
00849                 f_namei = tempStri + "stress.dat";
00850                 fclose( fopen( f_namei.c_str(), strMode.c_str() ) );
00851 
00852                 framesOnFile = 0;
00853 
00854         } else {
00855                 strMode = "a+b"; // update at the end
00856         
00857                 // save vertices
00858                 if( file = fopen( ((COME_MoleculesCartilage*)parent)->getDeformationFile().c_str(), strMode.c_str() ) ){
00859         
00860                         vector<COME_Vertex3D> *verts = ((COME_MoleculesCartilage*)parent)->getSurface()->getVerticesGlobalPt();
00861                         for ( int i = 0; i < verts->size(); i++ ){
00862 
00863                                 COME_Vertex3D vertex = verts->at(i);
00864                                 double array[3];
00865                                 array[0] = vertex.x;
00866                                 array[1] = vertex.y;
00867                                 array[2] = vertex.z;
00868                                 fwrite( array, sizeof(double), 3, file );
00869                         }
00870                         fclose(file);
00871                         framesOnFile++;
00872                 }
00873 
00874                 // save stress
00875                 string f_name = ((COME_MoleculesCartilage*)parent)->getDeformationFile();
00876                 string tempStr( f_name, 0, f_name.size()-10 );
00877                 f_name = tempStr + "stress.dat";
00878 
00879                 if( file = fopen( f_name.c_str(), strMode.c_str() ) ){
00880         
00881                         vector<COME_Vertex3D> *verts = ((COME_MoleculesCartilage*)parent)->getSurface()->getVerticesGlobalPt();
00882                         for ( int i = 0; i < verts->size(); i++ ){
00883 
00884                                 COME_Vertex3D vertex = verts->at(i);
00885                                 double array[1];
00886                                 array[0] = vertex.getStress();
00887                                 fwrite( array, sizeof(double), 1, file );
00888                         }
00889                         fclose(file);
00890                 }
00891         }
00892 }
00893 
00894 // visco //
00895 void
00896 COME_MoleculesTissue::flow( double timestep ){
00897 
00898         list<COME_Molecule*>::iterator iter;
00899 
00900         //if( (*iter)->getPermeability() != 0.0 ){
00901         
00902         for (iter = shape.begin(); iter != shape.end(); iter++ ){
00903         
00904                 double volume = (*iter)->getLiquidFraction() * ( 4.0 / 3.0 ) * M_PI * pow( (*iter)->getNominalRadius(), 3 ) * (*iter)->getLiquidLevel();
00905                 
00906                 double totalFlowingVolume = (*iter)->getStress() * volume * (*iter)->getPermeability() / timestep;
00907 
00908                 // Update radius of the current molecule from new calculated volume
00909                 //printf( "avant: %f ", (*iter)->getRadius() );
00910                 //(*iter)->setRadius( pow( ( 3.0 * ( (*iter)->getVolume() - totalFlowingVolume ) ) / ( 4.0 * M_PI ), 1.0/3.0 ) );
00911                 //printf( "apres: %f \n", (*iter)->getRadius() );
00912                 //getch();
00913 
00914 
00915                 //if( (*iter)->isOnSurface() ){
00916                 //      totalFlowingVolume = 0.8 * totalFlowingVolume;
00917                 //}
00918 
00919                 list<COME_MoleculeLink*>::iterator iterO;
00920                 double totalStress = 0.0;
00921                 double *stressArray = new double [(*iter)->getConnectionList()->size()];
00922                 int iStress;
00923                 for (iterO = (*iter)->getConnectionList()->begin(), iStress = 0; iterO != (*iter)->getConnectionList()->end(); iterO++, iStress++ ){
00924                         
00925                         COME_Molecule *temp;
00926                         if( (*iterO)->getElement(0) == (*iter) ){
00927                                 temp = (*iterO)->getElement(1);
00928                         } else {
00929                                 temp = (*iterO)->getElement(0);
00930                         }
00931 
00932                         totalStress += temp->getStress();
00933                         stressArray[iStress] = temp->getStress();
00934                 }
00935 
00936                 for (iterO = (*iter)->getConnectionList()->begin(), iStress = (*iter)->getConnectionList()->size()-1 ; iterO != (*iter)->getConnectionList()->end(); iterO++, iStress-- ){
00937                         
00938                         COME_Molecule *temp;
00939                         if( (*iterO)->getElement(0) == (*iter) ){
00940                                 temp = (*iterO)->getElement(1);
00941                         } else {
00942                                 temp = (*iterO)->getElement(0);
00943                         }
00944 
00945                         // 
00946                         if( totalStress ){
00947                                 (*iterO)->setFlowVolume( (*iter), ( stressArray[iStress] / totalStress ) * totalFlowingVolume );
00948                         }
00949 
00950                 }
00951                                                 
00952         }
00953 
00954         for (iter = shape.begin(); iter != shape.end(); iter++ ){
00955 
00956                 list<COME_MoleculeLink*>::iterator iterO;
00957                 double totalFlowingVolume = 0.0;
00958                 for (iterO = (*iter)->getConnectionList()->begin(); iterO != (*iter)->getConnectionList()->end(); iterO++ ){
00959                         
00960                         totalFlowingVolume += (*iterO)->getFlowVolume(*iter);
00961                 }
00962 
00963                 if( (*iter)->isOnSurface() ){
00964                         
00965                         totalFlowingVolume = 1.2 * totalFlowingVolume;
00966                 }
00967 
00968                 (*iter)->setRadius( pow( ( 3.0 * ( (*iter)->getVolume() + totalFlowingVolume ) ) / ( 4.0 * M_PI ), 1.0/3.0 ) );
00969 
00970                 for (iterO = (*iter)->getConnectionList()->begin(); iterO != (*iter)->getConnectionList()->end(); iterO++ ){
00971                 
00972                         COME_Molecule *temp;
00973                         if( (*iterO)->getElement(0) == (*iter) ){
00974                                 temp = (*iterO)->getElement(1);
00975                         } else {
00976                                 temp = (*iterO)->getElement(0);
00977                         }
00978                         
00979                         double incDist = fabs( ( (*iterO)->getViscousFraction() * ( temp->getRadius() + (*iter)->getRadius() ) ) - (*iterO)->getNominalDist() );
00980 
00981                         if( (*iterO)->getElement(1)->getPosition().vpDistance( (*iterO)->getElement(2)->getPosition() ) > (*iterO)->getNominalDist() ){
00982                                 (*iterO)->setNominalDist( (*iterO)->getNominalDist() - incDist );
00983                         } else {
00984                                 (*iterO)->setNominalDist( (*iterO)->getNominalDist() + incDist );
00985                         }
00986                 }
00987 
00988         }
00989 }
00990 
00991 
00996 void
00997 COME_MoleculesTissue :: transform( COME_Matrix M ){
00998 
00999         list<COME_Molecule*>::iterator iter;
01000         for( iter = shape.begin(); iter != shape.end(); iter++ ){
01001                 
01002                 if( (*iter)->isFixed() ){
01003 
01004                         (*iter)->setPosition( M * (*iter)->getPosition() );
01005 
01006                 }
01007         }
01008 }
01009 
01015 void
01016 COME_MoleculesTissue :: transform( COME_Matrix M, COME_Matrix MP ){
01017 
01018         list<COME_Molecule*>::iterator iter;
01019         for( iter = shape.begin(); iter != shape.end(); iter++ ){
01020                 
01021                 if( ( (*iter)->isFixed() ) && ( (*iter)->getDescription() == "origin" ) ){
01022 
01023                         (*iter)->setPosition( MP * (*iter)->getPosition() );
01024 
01025                 } else if( ( (*iter)->isFixed() ) && ( (*iter)->getDescription() == "free" ) ){
01026 
01027                         (*iter)->setPosition( M * (*iter)->getPosition() );
01028                 }
01029         }
01030 }
01031 
01036 void
01037 COME_MoleculesTissue :: scale( double factorx, double factory, double factorz ){
01038 
01039         // Find central point
01040         double accX = 0.0, accY = 0.0, accZ = 0.0;
01041         int n = 0;
01042         list<COME_Molecule*>::iterator iter;
01043         for (iter = shape.begin(); iter != shape.end(); iter++, n++ ){
01044                 accX+=(*iter)->getPosition().getX();
01045                 accY+=(*iter)->getPosition().getY();
01046                 accZ+=(*iter)->getPosition().getZ();
01047         }
01048         accX/=(double)n; accY/=(double)n; accZ/=(double)n;
01049         
01050         // Translate to origin
01051         translate( -accX, -accY, -accZ );
01052         
01053         //Apply scaling
01054         for (iter = shape.begin(); iter != shape.end(); iter++){
01055 
01056                 COME_Point3D newPos( (*iter)->getPosition().getX() * factorx, (*iter)->getPosition().getY() * factory, (*iter)->getPosition().getZ() * factorz );
01057                 (*iter)->setPosition( newPos );
01058                 (*iter)->setRadius( (*iter)->getRadius() * ( ( factorx + factory + factorz ) / 3.0 ) );
01059         }
01060         
01061         // Translate back to position
01062         translate( accX, accY, accZ );
01063 
01064 }
01065 
01069 void
01070 COME_MoleculesTissue :: translate( double dx, double dy, double dz ){
01071 
01072         list<COME_Molecule*>::iterator iter;
01073         for (iter = shape.begin(); iter != shape.end(); iter++ ){
01074         
01075                 COME_Point3D newPos( (*iter)->getPosition().getX() + dx, (*iter)->getPosition().getY() + dy, (*iter)->getPosition().getZ() + dz );
01076                 (*iter)->setPosition( newPos );
01077         }
01078 
01079 }
01080 
01086 void
01087 COME_MoleculesTissue :: rotate( double rx, double ry, double rz ){
01088 
01089         // Find central point
01090         double accX = 0.0, accY = 0.0, accZ = 0.0;
01091         int n = 0;
01092         list<COME_Molecule*>::iterator iter;
01093         for (iter = shape.begin(); iter != shape.end(); iter++, n++ ){
01094                 accX+=(*iter)->getPosition().getX();
01095                 accY+=(*iter)->getPosition().getY();
01096                 accZ+=(*iter)->getPosition().getZ();
01097         }
01098         accX/=(double)n; accY/=(double)n; accZ/=(double)n;
01099         
01100         // Translate to origin
01101         translate( -accX, -accY, -accZ );
01102         
01103         // Apply rotations
01104         for (iter = shape.begin(); iter != shape.end(); iter++){
01105 
01106                 COME_Point3D newPos;
01108                 newPos = COME_Point3D (
01109                         (*iter)->getPosition().getX(),
01110                         (*iter)->getPosition().getY() * cos( rx ) + (*iter)->getPosition().getZ() * sin( rx ),
01111                         (*iter)->getPosition().getZ() * cos( rx ) - (*iter)->getPosition().getY() * sin( rx )
01112                 );
01113                 (*iter)->setPosition( newPos );
01115                 newPos = COME_Point3D (
01116                         (*iter)->getPosition().getX() * cos( ry ) - (*iter)->getPosition().getZ() * sin( ry ),
01117                         (*iter)->getPosition().getY(),
01118                         (*iter)->getPosition().getX() * sin( ry ) + (*iter)->getPosition().getZ() * cos( ry )
01119                 );
01120                 (*iter)->setPosition( newPos );
01122                 newPos = COME_Point3D (
01123                         (*iter)->getPosition().getX() * cos( rz ) + (*iter)->getPosition().getY() * sin( rz ),
01124                         (*iter)->getPosition().getY() * cos( rz ) - (*iter)->getPosition().getX() * sin( rz ),
01125                         (*iter)->getPosition().getZ()
01126                 );
01127                 (*iter)->setPosition( newPos );
01128         }
01129         
01130         // Translate back to position
01131         translate( accX, accY, accZ );
01132 
01133 }
01134 
01135 void
01136 COME_MoleculesTissue::scaleNominals( double factor ){
01137 
01138         list<COME_Molecule*>::iterator iter;
01139         for (iter = shape.begin(); iter != shape.end(); iter++ ){
01140                 list<COME_MoleculeLink*>::iterator iterO;
01141                 for (iterO = (*iter)->getConnectionList()->begin(); iterO != (*iter)->getConnectionList()->end(); iterO++ ){
01142 
01143                         (*iterO)->setNominalDist( COME_Vector3D( (*iterO)->getElement(0)->getPosition() - (*iterO)->getElement(1)->getPosition() ).vpModule() * factor );
01144                 }
01145         }
01146 }
01147 
01151 void
01152 COME_MoleculesTissue::makeAllLocalFrames(){
01153         
01154         if (COME::flagSkinning ){
01155                 //printf( "Making all local frames\n" );
01156                 list<COME_Molecule*>::iterator iter;
01157                 for (iter = shape.begin(); iter != shape.end(); iter++){
01158                         (*iter)->makeLocalFrame();
01159                 }
01160         }
01161 }
01162 
01166 void
01167 COME_MoleculesTissue::checkInitialPositions(){
01168 
01169         list<COME_Molecule*>::iterator iter;
01170         for (iter = shape.begin(); iter != shape.end(); iter++){
01171                 (*iter)->checkPosition();
01172         }
01173 }
01174 
01178 void
01179 COME_MoleculesTissue::resetInitialPositions(){
01180 
01181         list<COME_Molecule*>::iterator iter;
01182         for (iter = shape.begin(); iter != shape.end(); iter++){
01183                 (*iter)->resetPosition();
01184         }
01185 }
comemoleculestissue.cpp
