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 
00031 
00032 #include <bio/comemoleculelink.h>
00033 #include <bio/comemoleculestissue.h>
00034 #include <bio/comemolecule.h>
00035 #include <math.h>
00036 
00040 
00041 COME_MoleculeLink::COME_MoleculeLink()
00042 {
00043         instanceCount++;
00044 }
00045 
00046 COME_MoleculeLink::COME_MoleculeLink( COME_Molecule *e1, COME_Molecule *e2 ){
00047 
00048         element1 = e1;
00049         element2 = e2;
00050 
00051         // set default values
00052         dampingConst = 0.0;
00053         hookeConst = 0.0;
00054         flowVolume1 = 0.0;
00055         flowVolume2 = 0.0;
00056         nominalDist = COME_Vector3D( e1->getPosition() - e2->getPosition() ).vpModule();
00057         viscousFraction = nominalDist / ( e1->getRadius() + e2->getRadius() );
00058         frictionConst = 0.0;
00059         instanceCount++;
00060         //printf( "Created link: %d\n", instanceCount );
00061 }
00062 
00063 COME_MoleculeLink::~COME_MoleculeLink(){
00064 
00065         instanceCount--;
00066         printf( "Destroied link: %d\n", instanceCount );
00067 }
00068 
00072 
00073 void
00074 COME_MoleculeLink::setDampingConst( double dampingN ){
00075 
00076         dampingConst = dampingN;
00077 }
00078 
00079 void
00080 COME_MoleculeLink::setHookeConst( double hookeN ){
00081 
00082         hookeConst = hookeN;
00083 }
00084 
00085 void
00086 COME_MoleculeLink::setNominalDist( double nominalN ){
00087 
00088         nominalDist = nominalN;
00089 }
00090 
00091 void
00092 COME_MoleculeLink::setViscousFraction( double fractionN ){
00093 
00094         viscousFraction = fractionN;
00095 }
00096 
00097 void
00098 COME_MoleculeLink::setFrictionConst( double frictionConstN ){
00099 
00100         frictionConst = frictionConstN;
00101 }
00102 
00103 void
00104 COME_MoleculeLink::setElement( int index, COME_Molecule* elementN ){
00105 
00106         if( index == 2 ){
00107                 element2 = elementN;
00108         } else {
00109                 element1 = elementN;
00110         }
00111 }
00112 
00113 void
00114 COME_MoleculeLink::setFlowVolume( COME_Molecule* element, double volumeNew ){
00115 
00116         if( element == element1 ){
00117                 flowVolume1 = volumeNew;
00118         } else if( element == element2 ) {
00119                 flowVolume2 = volumeNew;
00120         }
00121 }
00122 
00123 
00127 
00128 double
00129 COME_MoleculeLink::getDampingConst() const{
00130 
00131         return dampingConst;
00132 }
00133 
00134 double
00135 COME_MoleculeLink::getHookeConst() const{
00136 
00137         return hookeConst;
00138 }
00139 
00140 double
00141 COME_MoleculeLink::getNominalDist() const{
00142 
00143         return nominalDist;
00144 }
00145 
00146 
00147 double
00148 COME_MoleculeLink::getViscousFraction() const{
00149 
00150         return viscousFraction;
00151 }
00152 
00153 double
00154 COME_MoleculeLink::getFractureDist() const{
00155 
00156         return element1->getRadius() + element2->getRadius();
00157 }
00158 
00159 double
00160 COME_MoleculeLink::getFrictionConst() const{
00161 
00162         return frictionConst;
00163 }
00164 
00165 COME_Molecule*
00166 COME_MoleculeLink::getElement( int index ){
00167 
00168         if( index == 1 ){
00169                 return element2;
00170         }
00171         return element1;
00172 }
00173 
00174 COME_Molecule*
00175 COME_MoleculeLink::getOtherElement( COME_Molecule* thisElement ){
00176 
00177         if( thisElement == element1 ){
00178                 return element2;
00179         }
00180         return element1;
00181 }
00182 
00183 double
00184 COME_MoleculeLink::getFlowVolume( COME_Molecule* element ) const{
00185 
00186         if( element == element1 ){
00187                 return flowVolume2 - flowVolume1;
00188         } else if( element == element2 ) {
00189                 return flowVolume1 - flowVolume2;
00190         }
00191 
00192         return 0.0;
00193 
00194 }
00195 
00199 double
00200 COME_MoleculeLink::getLength() const {
00201 
00202         return element1->getPosition().vpDistance( element2->getPosition() );
00203 }
00204 
00205 
00209 vector<COME_Tetra>
00210 COME_MoleculeLink::getNeighbouringTetra() const {
00211 
00212         vector<COME_Tetra> tetras;
00213 
00214         list<COME_MoleculeLink*>::iterator iter1;
00215         int teste = 0;
00216         for (iter1 = element1->getConnectionList()->begin(); iter1 != element1->getConnectionList()->end(); iter1++, teste ++){
00217 
00218                 COME_Molecule* tetra2 = (*iter1)->getOtherElement(element1);
00219                 COME_Molecule* tetra3;
00220                 bool tetra2oK = false;
00221                 bool tetra3oK = false;
00222                 
00223                 if( tetra2 == element2 ){
00224                         continue;
00225                 }
00226 
00227                 list<COME_MoleculeLink*>::iterator iter2;
00228                 for (iter2 = tetra2->getConnectionList()->begin(); iter2 != tetra2->getConnectionList()->end(); iter2++){
00229                         
00230                         COME_Molecule* tetraTemp = (*iter2)->getOtherElement(tetra2);
00231                         if( tetraTemp == element2 ){
00232                                 tetra2oK = true;
00233                                 break;
00234                         }
00235                 }
00236 
00237                 if( tetra2oK ){
00238                         for (iter2 = tetra2->getConnectionList()->begin(); iter2 != tetra2->getConnectionList()->end(); iter2++){
00239                                 tetra3 = (*iter2)->getOtherElement(tetra2);
00240                                 if( ( tetra3 != element2 ) && ( tetra3 != element1 ) ){
00241                                         tetra3oK = true;
00242                                         break;
00243                                 }
00244                         }
00245 
00246                         if( tetra3oK ){
00247                                 COME_Tetra newTetra;
00248                                 newTetra.addMolecule( element1 ); 
00249                                 newTetra.addMolecule( tetra2 );
00250                                 newTetra.addMolecule( tetra3 );
00251                                 newTetra.addMolecule( element2 );
00252                                 tetras.push_back( newTetra );
00253                         }
00254                 }
00255 
00256         }
00257 
00258         return tetras;
00259 }
00260 
00264 void
00265 COME_MoleculeLink::recalculateSpring( int mode, double avgConnect ){
00266 
00267         if( mode == VOLUME_TETRAS ){
00268 
00269                 double totalVE = 0.0;
00270                 vector<COME_Tetra> tetras = getNeighbouringTetra();
00271                 for( int i= 0; i < tetras.size(); i++ ){
00272 
00273                         totalVE = tetras[i].getVolume()*tetras[i].getElasticity();
00274                 }
00275                 
00276                 // test when there is not any tetrahedron.
00277                 if( totalVE > 0.0 ){
00278                         printf( "old Hooke: %f\n", hookeConst );
00279                         hookeConst = totalVE * 1000 / ( getLength() * getLength() );
00280                         printf( "new Hooke: %f\n", hookeConst );
00281                 }
00282         } else if( mode == NUM_CONNECTIONS ){
00283 
00284                 int N1 = element1->getConnectionList()->size();
00285                 int N2 = element2->getConnectionList()->size();
00286 
00287                 int D1, D2;
00288 
00289                 if( N1 < 8 ){
00290                         D1 = N1;
00291                 } else if( ( N1 >= 9 ) && ( N1 < 12 ) ){
00292                         D1 = N1 - 1;
00293                 } else if( ( N1 >= 13 ) && ( N1 < 20 ) ){
00294                         D1 = N1 - 4;
00295                 } else if( N1 >= 20 ){
00296                         D1 = N1 /2;
00297                 }
00298 
00299                 if( N2 < 8 ){
00300                         D2 = N2;
00301                 } else if( ( N2 >= 9 ) && ( N2 < 12 ) ){
00302                         D2 = N2 - 1;
00303                 } else if( ( N2 >= 13 ) && ( N2 < 20 ) ){
00304                         D2 = N2 - 4;
00305                 } else if( N2 >= 20 ){
00306                         D2 = N2 /2;
00307                 }
00308 
00309                 double dd1 = (double)D1;
00310                 double dd2 = (double)D2;
00311 
00312                 double area1 = ( (2*element1->getRadius())*(2*element1->getRadius()) ) / (dd1/2);
00313                 double area2 = ( (2*element2->getRadius())*(2*element2->getRadius()) ) / (dd2/2);
00314                 
00315                 double hooke = ( ( ( element1->getMaterial()->getYoungsModulus()*area1 ) + ( element2->getMaterial()->getYoungsModulus()*area2 ) ) / 2 ) / getNominalDist();
00316 
00317                 hooke += (hooke / 3.0);
00318                 printf("old hooke: %f\n", hookeConst );
00319                 setHookeConst( hooke );
00320                 printf("new hooke: %f\n", hookeConst );
00321 
00322         } else if( mode == ANGLES ){
00323 
00324                 COME_Vector3D vThis( element2->getPosition() - element1->getPosition() );
00325 
00326                 list<COME_MoleculeLink*>::iterator iter1;
00327                 double discounts = 0.0;
00328                 double totalCos = 0.0;
00329                 for (iter1 = element1->getConnectionList()->begin(); iter1 != element1->getConnectionList()->end(); iter1++){
00330 
00331                         COME_Molecule* temp = (*iter1)->getOtherElement(element1);
00332                         if( temp != element2 ){
00333 
00334                                 
00335                                 double area1 = ( (2*element1->getRadius())*(2*element1->getRadius()) );
00336                                 double area2 = ( (2*temp->getRadius())*(2*temp->getRadius()) );
00337                 
00338                                 double hooke = ( ( ( element1->getMaterial()->getYoungsModulus()*area1 ) / (*iter1)->getNominalDist() ) + ( ( temp->getMaterial()->getYoungsModulus()*area2 ) / (*iter1)->getNominalDist() ) ) / 2;
00339 
00340                                 COME_Vector3D vCurr( temp->getPosition() - element1->getPosition() );
00341                                 if( vThis.getAngle( vCurr ) < (M_PI/2) ){
00342                                         double currCos = cos( vThis.getAngle( vCurr ) );
00343                                         if( vThis.vpModule() > vCurr.vpModule() ){
00344                                                 currCos = 1 / currCos;
00345                                         }
00346                                         totalCos += currCos * vCurr.vpModule();
00347                                         discounts += hooke * cos( vThis.getAngle( vCurr ) );
00348                                 }
00349 
00350                         } else {
00351                                 COME_Vector3D vCurr( element2->getPosition() - element1->getPosition() );
00352                                 totalCos += vCurr.vpModule();
00353                                 double area1 = ( (2*element1->getRadius())*(2*element1->getRadius()) );
00354                                 double area2 = ( (2*element2->getRadius())*(2*element2->getRadius()) );
00355                                 double hooke = ( ( ( element1->getMaterial()->getYoungsModulus()*area1 ) / getNominalDist() ) + ( ( element2->getMaterial()->getYoungsModulus()*area2 ) / getNominalDist() ) ) / 2;
00356                                 discounts += hooke;
00357                         }
00358                 }
00359 
00360                 COME_Vector3D vThis2( element1->getPosition() - element2->getPosition() );
00361 
00362                 for (iter1 = element2->getConnectionList()->begin(); iter1 != element2->getConnectionList()->end(); iter1++){
00363 
00364                         COME_Molecule* temp = (*iter1)->getOtherElement(element2);
00365                         if( temp != element1 ){
00366 
00367                                 
00368                                 double area1 = ( (2*element2->getRadius())*(2*element2->getRadius()) );
00369                                 double area2 = ( (2*temp->getRadius())*(2*temp->getRadius()) );
00370                 
00371                                 double hooke = ( ( ( element2->getMaterial()->getYoungsModulus()*area1 ) / (*iter1)->getNominalDist() ) + ( ( temp->getMaterial()->getYoungsModulus()*area2 ) / (*iter1)->getNominalDist() ) ) / 2;
00372 
00373                                 COME_Vector3D vCurr( temp->getPosition() - element2->getPosition() );
00374                                 if( vThis2.getAngle( vCurr ) < (M_PI/2) ){
00375                                         double currCos = cos( vThis2.getAngle( vCurr ) );
00376                                         if( vThis2.vpModule() > vCurr.vpModule() ){
00377                                                 currCos = 1 / currCos;
00378                                         }
00379                                         totalCos += currCos * vCurr.vpModule();
00380                                         discounts += hooke * cos( vThis2.getAngle( vCurr ) );
00381                                 }
00382 
00383                         } else {
00384                                 COME_Vector3D vCurr( element2->getPosition() - element1->getPosition() );
00385                                 totalCos += vCurr.vpModule();
00386                                 double area1 = ( (2*element1->getRadius())*(2*element1->getRadius()) );
00387                                 double area2 = ( (2*element2->getRadius())*(2*element2->getRadius()) );
00388                                 double hooke = ( ( ( element1->getMaterial()->getYoungsModulus()*area1 ) / getNominalDist() ) + ( ( element2->getMaterial()->getYoungsModulus()*area2 ) / getNominalDist() ) ) / 2;
00389                                 discounts += hooke;
00390                         }
00391                 }
00392 
00393                 double area1 = ( (2*element1->getRadius())*(2*element1->getRadius()) );
00394                 double area2 = ( (2*element2->getRadius())*(2*element2->getRadius()) );
00395                 double hooke = ( ( ( element1->getMaterial()->getYoungsModulus()*area1 ) / getNominalDist() ) + ( ( element2->getMaterial()->getYoungsModulus()*area2 ) / getNominalDist() ) ) / 2;
00396 
00397                 double newHooke =  getNominalDist() * ( discounts / totalCos );
00398                 printf( "newHooke: %f ; hooke: %f ; discounts: %f ; totalCos: %f \n", newHooke, hooke, discounts, totalCos );
00399 
00400                 setHookeConst( newHooke );
00401 
00402         } else if( mode == NUMERICALLY ){
00403 
00404                 double increment = avgConnect; // reuse parameter
00405                 setHookeConst( hookeConst * increment );
00406         }
00407 
00408 
00409         
00410 }
comemoleculelink.cpp
