Mehanny.cpp

00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** Developed by:                                                      **
00015 **   Frank McKenna (fmckenna@ce.berkeley.edu)                         **
00016 **   Gregory L. Fenves (fenves@ce.berkeley.edu)                       **
00017 **   Filip C. Filippou (filippou@ce.berkeley.edu)                     **
00018 **                                                                    **
00019 ** ****************************************************************** */
00020                                                                         
00021 // $Revision: 1.1 $
00022 // $Date: 2004/09/01 03:54:28 $
00023 // $Source: /usr/local/cvs/OpenSees/SRC/damage/Mehanny.cpp,v $
00024 
00025 // Written: Arash Altoontash, Gregory Deierlein 
00026 // Created: 10/02
00027 // Revision: AA
00028 //
00029 // Description: This file contains the class implementation for Mehanny 
00030 // damage model. Mehanny damage model calculates the damage index based on
00031 // primary half cycle and the summation of follower half cycles.
00032 //
00033 
00034 #include <Mehanny.h>
00035 #include <DamageResponse.h>
00036 
00037 
00038 #define DEBG   0
00039 
00040 Mehanny::Mehanny(int tag, double alpha, double beta, double gamma,
00041                  double ultimatePosValue , double ultimateNegValue, double abstol, double reltol, 
00042                  double posmodifier, double negmodifier)
00043   :DamageModel(tag,DMG_TAG_Mehanny),
00044    Alpha(alpha) , Beta(beta) , Gamma(gamma) , UltimatePosValue(ultimatePosValue) , UltimateNegValue(ultimateNegValue),
00045    PosModifier(posmodifier), NegModifier(negmodifier),AbsTol(abstol), RelTol(reltol)
00046 {
00047   if ( UltimatePosValue<=0 || Alpha<0 || Beta<0 || Gamma<0 )
00048     opserr << "CumulativePeak::CumulativePeak : Incorrect arguments for the damage model";
00049   
00050   if ( UltimateNegValue == 0.0 ) { 
00051     UltimateNegValue = UltimatePosValue; 
00052   } else { 
00053     UltimateNegValue = fabs ( UltimateNegValue ) ; 
00054   }
00055   
00056   if ( AbsTol < 0.0 ) AbsTol = 1.0;
00057   if ( RelTol < 0.0 ) RelTol = 1.0;
00058   if ( PosModifier < 0.0 ) PosModifier = 1.0;
00059   if ( NegModifier < 0.0 ) NegModifier = 1.0;
00060   
00061   this->revertToStart();
00062 
00063   /*
00064   if ( DEBG ==1 )
00065     {
00066       // Open an output file for debugging
00067       char FileName[20];                                                // debugging
00068       sprintf(FileName, "Mehanny%d.out", tag);
00069       OutputFile = fopen ( FileName , "w" );    // debugging
00070       fprintf( OutputFile , "\t dp\t\t pcyc\t\t pFHC\t\t pPHC\t\t ncyc\t\t nFHC\t\t nPHC\t\t dmg \n" ) ;
00071     }
00072   */
00073 }
00074 
00075 Mehanny::Mehanny()
00076 :DamageModel(0,DMG_TAG_Mehanny)
00077 {
00078   // Does nothing
00079 }
00080 
00081 Mehanny::~Mehanny()
00082 {
00083   /*
00084   if ( DEBG == 1 )
00085     {
00086       fclose( OutputFile );                                             // debugging
00087     }
00088   */
00089 }
00090 
00091 
00092 int Mehanny::setTrial (Vector trialVector )
00093 {
00094   if ( trialVector.Size() != 3 ) {
00095     opserr << "WARNING: Mehanny::setTrial Wrong vector size for trial data" << endln;
00096     return -1;
00097         }
00098   
00099   double TrialDefo = trialVector(0);
00100   double TrialForce = trialVector(1);
00101   double TrialKU = trialVector(2);
00102   
00103   double TrialScalar = 0.0;
00104   
00105   // calculate the plastic deformation once the unloading stiffness is defined as non zero
00106   if ( TrialKU != 0.0 )
00107     {
00108       TrialScalar = TrialDefo - TrialForce/TrialKU;
00109     } else {
00110       
00111       // use the elastic deformation instead
00112       TrialScalar = TrialDefo;
00113     }
00114   
00115   return this->processData( TrialScalar );
00116 }
00117 
00118 
00119 double
00120 Mehanny::getDamage (void)
00121 {
00122   double PosDamage , NegDamage , OveralDamage;
00123   
00124   PosDamage = ( pow(TrialPosPHC,Alpha) + pow(TrialSumPosFHC,Beta) ) / ( pow(UltimatePosValue,Alpha) + pow(TrialSumPosFHC,Beta) ) ;
00125   
00126   NegDamage = ( pow(fabs(TrialNegPHC),Alpha) + pow(fabs(TrialSumNegFHC),Beta) ) / ( pow(fabs(UltimateNegValue),Alpha) + pow(fabs(TrialSumNegFHC),Beta) ) ;
00127   
00128   OveralDamage = pow( ( pow(PosDamage,Gamma) + pow(NegDamage,Gamma) ) , 1/Gamma );
00129   
00130   if ( OveralDamage < CommDamage ) OveralDamage = CommDamage;
00131   return OveralDamage;
00132 }
00133 
00134 
00135 double Mehanny::getPosDamage (void)
00136 {
00137   double PosDamage , NegDamage , OveralDamage;
00138   
00139   PosDamage = ( pow(TrialPosPHC,Alpha) + pow(TrialSumPosFHC,Beta) ) / ( pow(UltimatePosValue,Alpha) + pow(TrialSumPosFHC,Beta) ) ;
00140   
00141   NegDamage = ( pow(fabs(TrialNegPHC),Alpha) + pow(fabs(TrialSumNegFHC),Beta) ) / ( pow(fabs(UltimateNegValue),Alpha) + pow(fabs(TrialSumNegFHC),Beta) ) ;;
00142 
00143   OveralDamage = pow( ( 1.0 * pow(PosDamage,Gamma) + NegModifier * pow(NegDamage,Gamma) ) , 1/Gamma );
00144   
00145   return OveralDamage;
00146 }
00147 
00148 
00149 double Mehanny::getNegDamage (void)
00150 {
00151   double PosDamage , NegDamage , OveralDamage;
00152   
00153   PosDamage = ( pow(TrialPosPHC,Alpha) + pow(TrialSumPosFHC,Beta) ) / ( pow(UltimatePosValue,Alpha) + pow(TrialSumPosFHC,Beta) ) ;
00154   
00155   NegDamage = ( pow(fabs(TrialNegPHC),Alpha) + pow(fabs(TrialSumNegFHC),Beta) ) / ( pow(fabs(UltimateNegValue),Alpha) + pow(fabs(TrialSumNegFHC),Beta) ) ;;
00156   
00157   OveralDamage = pow( ( PosModifier * pow(PosDamage,Gamma) + 1.0 * pow(NegDamage,Gamma) ) , 1/Gamma );
00158   
00159   return OveralDamage;
00160 }
00161 
00162     
00163 int
00164 Mehanny::commitState (void)
00165 {
00166   /*
00167   if ( DEBG ==1 )
00168     {           
00169       fprintf( OutputFile , "\t %f\t %f\t %f\t %f\t %f\t %f\t %f\t %f \n",      TrialPlasticDefo, TrialPosCycle, TrialSumPosFHC, TrialPosPHC, TrialNegCycle, TrialSumNegFHC, TrialNegPHC, this->getDamage() ) ;
00170     }
00171   */
00172   LCommPlasticDefo = CommPlasticDefo;
00173   LCommDefoIncr = CommDefoIncr;
00174   LCommTempPDefo = CommTempPDefo;
00175   LCommPosCycle = CommPosCycle;
00176   LCommNegCycle = CommNegCycle; 
00177   LCommSumPosFHC = CommSumPosFHC;
00178   LCommPosPHC = CommPosPHC;
00179   LCommSumNegFHC = CommSumNegFHC;
00180   LCommNegPHC = CommNegPHC;
00181   LCommDamage = CommDamage;
00182   
00183   CommPlasticDefo = TrialPlasticDefo;
00184   CommDefoIncr = TrialDefoIncr;
00185   CommTempPDefo = TrialTempPDefo;
00186   CommPosCycle = TrialPosCycle;
00187   CommNegCycle = TrialNegCycle; 
00188   CommSumPosFHC = TrialSumPosFHC;
00189   CommPosPHC = TrialPosPHC;
00190   CommSumNegFHC = TrialSumNegFHC;
00191   CommNegPHC = TrialNegPHC;
00192   CommDamage = TrialDamage;
00193   
00194   return 0;
00195 }
00196 
00197 int
00198 Mehanny::revertToLastCommit (void)
00199 {
00200   CommPlasticDefo = LCommPlasticDefo;
00201   CommDefoIncr = LCommDefoIncr;
00202   CommTempPDefo = LCommTempPDefo;
00203   CommPosCycle = LCommPosCycle;
00204   CommNegCycle = LCommNegCycle; 
00205   CommSumPosFHC = LCommSumPosFHC;
00206   CommPosPHC = LCommPosPHC;
00207   CommSumNegFHC = LCommSumNegFHC;
00208   CommNegPHC = LCommNegPHC;
00209   CommDamage = LCommDamage;
00210   
00211   return 0;
00212 }
00213 
00214 int
00215 Mehanny::revertToStart (void)
00216 {
00217   CommPlasticDefo = LCommPlasticDefo = 0.0;
00218   CommDefoIncr = LCommDefoIncr = 0.0;
00219   CommTempPDefo = LCommTempPDefo = 0.0;
00220   CommPosCycle = LCommPosCycle = 0.0;
00221   CommNegCycle = LCommNegCycle = 0.0;   
00222   CommSumPosFHC = LCommSumPosFHC = 0.0;
00223   CommPosPHC = LCommPosPHC = 0.0;
00224   CommSumNegFHC = LCommSumNegFHC = 0.0;
00225   CommNegPHC = LCommNegPHC = 0.0;
00226   CommDamage = LCommDamage = 0.0;
00227   
00228   return 0;
00229 }
00230 
00231 DamageModel*
00232 Mehanny::getCopy (void)
00233 {
00234   Mehanny *theCopy = new Mehanny(this->getTag(), Alpha , Beta , Gamma , UltimatePosValue , UltimateNegValue, AbsTol, RelTol, PosModifier, NegModifier );
00235   
00236   theCopy->TrialPlasticDefo = TrialPlasticDefo;
00237   theCopy->TrialDefoIncr = TrialDefoIncr;
00238   theCopy->TrialTempPDefo = TrialTempPDefo;
00239   theCopy->TrialPosCycle = TrialPosCycle;
00240   theCopy->TrialNegCycle = TrialNegCycle;       
00241   theCopy->TrialSumPosFHC = TrialSumPosFHC;
00242   theCopy->TrialPosPHC = TrialPosPHC;
00243   theCopy->TrialSumNegFHC = TrialSumNegFHC;
00244   theCopy->TrialNegPHC = TrialNegPHC;
00245   theCopy->TrialDamage = TrialDamage;
00246   
00247   // Commited state
00248   theCopy->CommPlasticDefo = CommPlasticDefo;
00249   theCopy->CommDefoIncr= CommDefoIncr;
00250   theCopy->CommTempPDefo = CommTempPDefo;
00251   theCopy->CommPosCycle = CommPosCycle;
00252   theCopy->CommNegCycle = CommNegCycle; 
00253   theCopy->CommSumPosFHC = CommSumPosFHC;
00254   theCopy->CommPosPHC = CommPosPHC;
00255   theCopy->CommSumNegFHC = CommSumNegFHC;
00256   theCopy->CommNegPHC = CommNegPHC;
00257   theCopy->CommDamage = CommDamage;
00258   
00259   // Last commit
00260   theCopy->LCommPlasticDefo = LCommPlasticDefo;
00261   theCopy->LCommDefoIncr = LCommDefoIncr;
00262   theCopy->LCommTempPDefo = LCommTempPDefo;
00263   theCopy->LCommPosCycle = LCommPosCycle;
00264   theCopy->LCommNegCycle = LCommNegCycle;       
00265   theCopy->LCommSumPosFHC = LCommSumPosFHC ;
00266   theCopy->LCommPosPHC = LCommPosPHC;
00267   theCopy->LCommSumNegFHC = LCommSumNegFHC;
00268   theCopy->LCommNegPHC = LCommNegPHC;
00269   theCopy->LCommDamage = LCommDamage;
00270   
00271   return theCopy;
00272 }
00273 
00274 int
00275 Mehanny::setVariable(const char *argv)
00276 {
00277   return -1;
00278 }
00279 
00280 int
00281 Mehanny::getVariable(int variableID, double &info)
00282 {
00283   return -1;
00284 }
00285 
00286 int
00287 Mehanny::setParameter(char **argv, int argc, Information &eleInformation)
00288 {
00289   return -1;
00290 }
00291 
00292 int
00293 Mehanny::updateParameter(int responseID, Information &eleInformation)
00294 {
00295   return -1;
00296 }
00297 
00298 Response*
00299 Mehanny::setResponse(char **argv, int argc, Information &info)
00300 {
00301   //
00302   // we compare argv[0] for known response types for the Truss
00303   //
00304   
00305   if ( strcmp(argv[0],"damage") == 0 || strcmp(argv[0],"damageindex") == 0 )
00306     return new DamageResponse( this , 1 , 0.0 );
00307   
00308   else if (strcmp(argv[0],"Value") == 0 || strcmp(argv[0],"defo") == 0 || strcmp(argv[0],
00309                                                                                  "deformation") == 0 )
00310     return new DamageResponse( this , 2 , 0.0 );
00311   
00312   else if (strcmp(argv[0],"trial") == 0 || strcmp(argv[0],"trialinfo") == 0 )
00313     return new DamageResponse( this , 3 , Vector(4) );
00314   
00315   else 
00316     return 0;
00317   
00318 }
00319 
00320 int 
00321 Mehanny::getResponse(int responseID, Information &info)
00322 {
00323   switch (responseID) {
00324   case -1:
00325     return -1;
00326     
00327   case 1:
00328     return info.setDouble( this->getDamage() );
00329     
00330   case 2:
00331     return info.setDouble( TrialPlasticDefo );
00332     
00333   case 3:
00334     if(info.theVector!=0)
00335       {
00336         (*(info.theVector))(0) = TrialPosPHC;
00337         (*(info.theVector))(1) = TrialSumPosFHC;
00338         (*(info.theVector))(2) = TrialNegPHC;
00339         (*(info.theVector))(3) = TrialSumNegFHC;
00340       }
00341     return 0;
00342     
00343   default:
00344     return -1;
00345   }
00346 }
00347 
00348 
00349 int
00350 Mehanny::sendSelf(int commitTag, Channel &theChannel)
00351 {
00352   return 0;
00353 }
00354 
00355 
00356 int
00357 Mehanny::recvSelf(int commitTag, Channel &theChannel,
00358                   FEM_ObjectBroker &theBroker)
00359 {
00360   return 0;
00361 }
00362 
00363 
00364 void
00365 Mehanny::Print(OPS_Stream &s, int flag )
00366 {
00367   s << "CumulativePeak tag: " << this->getTag() << endln;
00368   s << "  Alpha: " << Alpha << " Beta: " << Beta << "  Gamma: " << Gamma << endln;
00369   s << " UltimatePosValue: " << UltimatePosValue << " UltimateNegValue: " << UltimateNegValue << endln;
00370 }
00371 
00372 
00373 // Perivate functions
00374 int Mehanny::processData (double PDefo)
00375 {
00376   TrialPlasticDefo = PDefo;
00377   TrialDefoIncr = PDefo - CommPlasticDefo;
00378   TrialTempPDefo = CommTempPDefo;
00379   TrialPosCycle = CommPosCycle;
00380   TrialNegCycle = CommNegCycle;
00381   TrialSumPosFHC = CommSumPosFHC;
00382   TrialPosPHC = CommPosPHC;
00383   TrialSumNegFHC = CommSumNegFHC;
00384   TrialNegPHC = CommNegPHC;
00385   TrialDamage = CommDamage;
00386   
00387   if ( TrialDefoIncr != 0.0 ) 
00388     {
00389       // For a non-zero step
00390       if ( ( (TrialDefoIncr >= AbsTol) && (TrialDefoIncr >= RelTol*TrialPosPHC) ) ||
00391            ( (TrialDefoIncr+TrialTempPDefo) >= AbsTol && (TrialDefoIncr+TrialTempPDefo) >= RelTol*TrialPosPHC )  ||
00392            ( (TrialDefoIncr <= -AbsTol) && (TrialDefoIncr >= -RelTol*TrialPosPHC) ) ||
00393            ( (TrialDefoIncr+TrialTempPDefo) <= -AbsTol && (TrialDefoIncr+TrialTempPDefo) <= -RelTol*TrialPosPHC )       )
00394         {
00395           // in case the plastic deformation increament is significant enough
00396           // or the current increment is larger than the tolerence
00397           
00398           if ( TrialPosCycle == 0.0 && TrialNegCycle == 0.0 )
00399             {
00400               // For a brand new cycle
00401               if ( TrialDefoIncr > 0.0 )
00402                 { 
00403                   TrialPosCycle = TrialDefoIncr;
00404                 } 
00405               else
00406                 { 
00407                   TrialNegCycle = TrialDefoIncr;
00408                 }
00409             }
00410           else if ( TrialPosCycle > 0.0 && TrialNegCycle == 0.0 )
00411             {
00412               // Check the status for a positive half cycle
00413               // Determine the status of this step, if a new cycle has started
00414               // or still on the last cycle
00415               //
00416               if ( TrialDefoIncr + TrialTempPDefo >= 0.0 )
00417                 {
00418                   // just add the temporarily saved and recent cycle to the current positive cycle
00419                   TrialPosCycle = TrialPosCycle + TrialDefoIncr + TrialTempPDefo;                       
00420                 }
00421               else
00422                 {
00423                   // end this positive half cycle and initiate a new negative half cycle
00424                   TrialPosCycle = 0.0;
00425                   TrialNegCycle = TrialDefoIncr + TrialTempPDefo;       
00426                 }
00427             }
00428           else if ( TrialPosCycle == 0.0 && TrialNegCycle < 0.0 )
00429             {
00430               // Check the status for a negative half cycle
00431               // Determine the status of this step, if a new cycle has started
00432               // or still on the last cycle
00433               //
00434               
00435               if ( TrialDefoIncr+TrialTempPDefo <= 0.0 )
00436                 {
00437                   // just add the temporarily saved and recent cycle to the current negative cycle
00438                   TrialNegCycle = TrialNegCycle + TrialDefoIncr + TrialTempPDefo;                       
00439                 } else
00440                   {
00441                     // end this negative half cycle and initiate a new positive half cycle
00442                     TrialNegCycle = 0.0;
00443                     TrialPosCycle = TrialDefoIncr + TrialTempPDefo;     
00444                   }
00445             }
00446           else 
00447             {
00448               // This is an internal error case
00449               opserr << "Mehanny::processData :Error, Can not detect a half cycle" << endln;
00450               return -1;
00451             }
00452           
00453           // reset the temporary plastic deformation
00454           TrialTempPDefo = 0.0;
00455         } else 
00456           {
00457             // for very small increments, the increment is only added to the temporary plastic deformation
00458             TrialTempPDefo = TrialTempPDefo + TrialDefoIncr;
00459             
00460           }
00461       
00462       // Process the current step
00463       // now detect the peaks
00464       if ( TrialPosCycle > 0.0 && TrialNegCycle == 0.0 )
00465         {
00466           // deal with a positive half cycle
00467           if ( TrialPosCycle > TrialPosPHC)
00468             {
00469               TrialPosPHC = TrialPosCycle;
00470             }
00471           else 
00472             {
00473               TrialSumPosFHC = TrialSumPosFHC - CommPosCycle + TrialPosCycle;
00474             }
00475         }
00476       else if ( TrialPosCycle == 0.0 && TrialNegCycle < 0.0 )
00477         {
00478           // deal with a negative half cycle
00479           if ( TrialNegCycle < TrialNegPHC)
00480             {
00481               TrialNegPHC = TrialNegCycle;
00482             }
00483           else 
00484             {
00485               TrialSumNegFHC = TrialSumNegFHC - CommNegCycle + TrialNegCycle;
00486             }
00487         }
00488     }
00489   
00490   return 0;
00491 }
00492 
00493 
00494 int Mehanny::setInputResponse ( Element *elem , const char **argv , int argc, int ndof )
00495 {
00496   return -1;
00497 }
Generated on Mon Oct 23 15:05:00 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation