PressureDependMultiYield.cpp

00001 // $Revision: 1.36 $
00002 // $Date: 2006/08/04 18:29:48 $
00003 // $Source: /usr/local/cvs/OpenSees/SRC/material/nD/soil/PressureDependMultiYield.cpp,v $
00004                                                                         
00005 // Written: ZHY
00006 // Created: August 2000
00007 
00008 //
00009 // PressureDependMultiYield.cpp
00010 // -------------------
00011 //
00012 
00013 #include <math.h>
00014 #include <stdlib.h>
00015 #include <PressureDependMultiYield.h>
00016 #include <Information.h>
00017 #include <ID.h>
00018 #include <MaterialResponse.h>
00019 
00020 int PressureDependMultiYield::matCount=0;
00021 int* PressureDependMultiYield::loadStagex = 0;  //=0 if elastic; =1 if plastic
00022 int* PressureDependMultiYield::ndmx=0;  //num of dimensions (2 or 3)
00023 double* PressureDependMultiYield::rhox=0;
00024 double* PressureDependMultiYield::refShearModulusx=0;
00025 double* PressureDependMultiYield::refBulkModulusx=0;
00026 double* PressureDependMultiYield::frictionAnglex=0;
00027 double* PressureDependMultiYield::peakShearStrainx=0;
00028 double* PressureDependMultiYield::refPressurex=0;
00029 double* PressureDependMultiYield::cohesionx=0;
00030 double* PressureDependMultiYield::pressDependCoeffx=0;
00031 int*    PressureDependMultiYield::numOfSurfacesx=0;
00032 double* PressureDependMultiYield::phaseTransfAnglex=0; 
00033 double* PressureDependMultiYield::contractParam1x=0;
00034 double* PressureDependMultiYield::dilateParam1x=0;
00035 double* PressureDependMultiYield::dilateParam2x=0;
00036 double* PressureDependMultiYield::liquefyParam1x=0;
00037 double* PressureDependMultiYield::liquefyParam2x=0;
00038 double* PressureDependMultiYield::liquefyParam4x=0;
00039 double* PressureDependMultiYield::einitx=0;    //initial void ratio
00040 double* PressureDependMultiYield::volLimit1x=0;
00041 double* PressureDependMultiYield::volLimit2x=0;
00042 double* PressureDependMultiYield::volLimit3x=0;
00043 double* PressureDependMultiYield::residualPressx=0;
00044 double* PressureDependMultiYield::stressRatioPTx=0;
00045 double* PressureDependMultiYield::Hvx=0;
00046 double* PressureDependMultiYield::Pvx=0;
00047 
00048 double PressureDependMultiYield::pAtm = 101.;
00049 Matrix PressureDependMultiYield::theTangent(6,6);
00050 T2Vector PressureDependMultiYield::trialStrain;
00051 T2Vector PressureDependMultiYield::subStrainRate;
00052 Vector PressureDependMultiYield::workV6(6);
00053 T2Vector PressureDependMultiYield::workT2V;
00054 
00055 const   double pi = 3.14159265358979;
00056 
00057 PressureDependMultiYield::PressureDependMultiYield (int tag, int nd, 
00058                                                     double r, double refShearModul,
00059                                                     double refBulkModul, double frictionAng,
00060                                                     double peakShearStra, double refPress, 
00061                                                     double pressDependCoe,
00062                                                     double phaseTransformAng, 
00063                                                     double contractionParam1,
00064                                                     double dilationParam1,
00065                                                     double dilationParam2,
00066                                                     double liquefactionParam1,
00067                                                     double liquefactionParam2,
00068                                                     double liquefactionParam4,
00069                                                     int numberOfYieldSurf, 
00070                                                                 double * gredu,
00071                                                     double ei,
00072                                                     double volLim1, double volLim2, double volLim3,
00073                                                     double atm, double cohesi,
00074                                                         double hv, double pv)
00075  : NDMaterial(tag,ND_TAG_PressureDependMultiYield), currentStress(),
00076    trialStress(), currentStrain(), strainRate(),
00077    reversalStress(), PPZPivot(), PPZCenter(), 
00078    lockStress(), reversalStressCommitted(), 
00079    PPZPivotCommitted(), PPZCenterCommitted(),
00080    lockStressCommitted()
00081 {
00082   if (nd !=2 && nd !=3) {
00083     opserr << "FATAL:PressureDependMultiYield:: dimension error" << endln;
00084     opserr << "Dimension has to be 2 or 3, you give nd= " << nd << endln;
00085    exit(-1);
00086   }
00087   if (refShearModul <= 0) {
00088     opserr << "FATAL:PressureDependMultiYield:: refShearModulus <= 0" << endln;
00089    exit(-1);
00090   }
00091   if (refBulkModul <= 0) {
00092     opserr << "FATAL:PressureDependMultiYield:: refBulkModulus <= 0" << endln;
00093    exit(-1);
00094   }
00095   if (frictionAng <= 0.) {
00096     opserr << "FATAL:PressureDependMultiYield:: frictionAngle <= 0" << endln;
00097    exit(-1);
00098   }
00099   if (frictionAng >= 90.) {
00100     opserr << "FATAL:PressureDependMultiYield:: frictionAngle >= 90" << endln;
00101    exit(-1);
00102   }
00103   if (phaseTransformAng <= 0.) {
00104     opserr << "FATAL:PressureDependMultiYield:: phaseTransformAng <= 0" << endln;
00105    exit(-1);
00106   }
00107   if (phaseTransformAng > frictionAng) {
00108     opserr << "WARNING:PressureDependMultiYield:: phaseTransformAng > frictionAng" << endln;
00109     opserr << "Will set phaseTransformAng = frictionAng." <<endln;
00110     phaseTransformAng = frictionAng;
00111   }
00112   if (cohesi < 0) {
00113     opserr << "WARNING:PressureDependMultiYield:: cohesion < 0" << endln;
00114     opserr << "Will reset cohesion to zero." << endln;
00115     cohesi = 0.;
00116   } 
00117   if (peakShearStra <= 0) {
00118     opserr << "FATAL:PressureDependMultiYield:: peakShearStra <= 0" << endln;
00119    exit(-1);
00120   }
00121   if (refPress <= 0) {
00122     opserr << "FATAL:PressureDependMultiYield:: refPress <= 0" << endln;
00123    exit(-1);
00124   }
00125   if (pressDependCoe < 0) {
00126     opserr << "WARNING:PressureDependMultiYield:: pressDependCoe < 0" << endln;
00127     opserr << "Will reset pressDependCoe to zero." << endln;
00128     pressDependCoe = 0.;
00129   }
00130   if (numberOfYieldSurf <= 0) {
00131     opserr << "WARNING:PressureDependMultiYield:: numberOfSurfaces <= 0" << endln;
00132     opserr << "Will use 10 yield surfaces." << endln;
00133     numberOfYieldSurf = 10;
00134   }
00135   if (numberOfYieldSurf > 40) {
00136     opserr << "WARNING:PressureDependMultiYield::PressureDependMultiYield: numberOfSurfaces > 40" << endln;
00137     opserr << "Will use 40 yield surfaces." << endln;
00138     numberOfYieldSurf = 40;
00139   }
00140   if (volLim1 < 0) {
00141     opserr << "WARNING:PressureDependMultiYield:: volLim1 < 0" << endln;
00142     opserr << "Will reset volLimit to 0.8" << endln;
00143     volLim1 = 0.8;
00144   }
00145   if (r < 0) {
00146     opserr << "FATAL:PressureDependMultiYield:: rho <= 0" << endln;
00147    exit(-1);
00148   }
00149   if (ei < 0) {
00150     opserr << "FATAL:PressureDependMultiYield:: e <= 0" << endln;
00151    exit(-1);
00152   }
00153 
00154   if (matCount%20 == 0) {
00155      int * temp1 = loadStagex;
00156          int * temp2 = ndmx;
00157          double * temp3 = rhox;
00158      double * temp4 = refShearModulusx;
00159      double * temp5 = refBulkModulusx;
00160      double * temp6 = frictionAnglex;
00161      double * temp7 = peakShearStrainx;
00162      double * temp8 = refPressurex;
00163      double * temp9 = cohesionx;
00164      double * temp10 = pressDependCoeffx;
00165          int * temp11 = numOfSurfacesx;
00166      double * temp12 = residualPressx;
00167      double * temp13 = phaseTransfAnglex; 
00168      double * temp14 = contractParam1x;
00169      double * temp15 = dilateParam1x;
00170      double * temp16 = dilateParam2x;
00171      double * temp17 = liquefyParam1x;
00172      double * temp18 = liquefyParam2x;
00173      double * temp19 = liquefyParam4x;
00174      double * temp20 = einitx;    //initial void ratio
00175      double * temp21 = volLimit1x;
00176      double * temp22 = volLimit2x;
00177      double * temp23 = volLimit3x;
00178      double * temp24 = stressRatioPTx;
00179          double * temp25 = Hvx;
00180          double * temp26 = Pvx;
00181 
00182      loadStagex = new int[matCount+20];
00183      ndmx = new int[matCount+20];
00184      rhox = new double[matCount+20];
00185      refShearModulusx = new double[matCount+20];
00186      refBulkModulusx = new double[matCount+20];
00187      frictionAnglex = new double[matCount+20];
00188      peakShearStrainx = new double[matCount+20];
00189      refPressurex = new double[matCount+20];
00190          cohesionx = new double[matCount+20];
00191      pressDependCoeffx = new double[matCount+20];
00192      numOfSurfacesx = new int[matCount+20];
00193      residualPressx = new double[matCount+20];
00194      phaseTransfAnglex = new double[matCount+20]; 
00195      contractParam1x = new double[matCount+20];
00196      dilateParam1x = new double[matCount+20];
00197      dilateParam2x = new double[matCount+20];
00198      liquefyParam1x = new double[matCount+20];
00199      liquefyParam2x = new double[matCount+20];
00200      liquefyParam4x = new double[matCount+20];
00201      einitx = new double[matCount+20];    //initial void ratio
00202      volLimit1x = new double[matCount+20];
00203      volLimit2x = new double[matCount+20];
00204      volLimit3x = new double[matCount+20];
00205      stressRatioPTx = new double[matCount+20];
00206          Hvx = new double[matCount+20];
00207          Pvx = new double[matCount+20];
00208 
00209          for (int i=0; i<matCount; i++) {
00210          loadStagex[i] = temp1[i];
00211                  ndmx[i] = temp2[i];
00212          rhox[i] = temp3[i];
00213              refShearModulusx[i] = temp4[i];
00214          refBulkModulusx[i] = temp5[i];
00215          frictionAnglex[i] = temp6[i];
00216          peakShearStrainx[i] = temp7[i];
00217          refPressurex[i] = temp8[i];
00218          cohesionx[i] = temp9[i];
00219          pressDependCoeffx[i] = temp10[i];
00220          numOfSurfacesx[i] = temp11[i];
00221          residualPressx[i] = temp12[i];
00222          phaseTransfAnglex[i] = temp13[i]; 
00223          contractParam1x[i] = temp14[i];
00224          dilateParam1x[i] = temp15[i];
00225          dilateParam2x[i] = temp16[i];
00226          liquefyParam1x[i] = temp17[i];
00227          liquefyParam2x[i] = temp18[i];
00228          liquefyParam4x[i] = temp19[i];
00229          einitx[i] = temp20[i];    //initial void ratio
00230          volLimit1x[i] = temp21[i];
00231          volLimit2x[i] = temp22[i];
00232          volLimit3x[i] = temp23[i];
00233          stressRatioPTx[i] = temp24[i];
00234                  Hvx[i] = temp25[i];
00235                  Pvx[i] = temp26[i];
00236      }
00237      
00238          if (matCount > 0) {
00239              delete [] temp1; delete [] temp2; delete [] temp3; delete [] temp4; 
00240              delete [] temp5; delete [] temp6; delete [] temp7; delete [] temp8; 
00241              delete [] temp9; delete [] temp10; delete [] temp11; delete [] temp12; 
00242              delete [] temp13; delete [] temp14; delete [] temp15; delete [] temp16; 
00243              delete [] temp17; delete [] temp18; delete [] temp19; delete [] temp20; 
00244              delete [] temp21; delete [] temp22; delete [] temp23; delete [] temp24; 
00245          delete [] temp25; delete [] temp26; 
00246      }
00247   }
00248 
00249   ndmx[matCount] = nd;
00250   loadStagex[matCount] = 0;   //default
00251   refShearModulusx[matCount] = refShearModul;
00252   refBulkModulusx[matCount] = refBulkModul;
00253   frictionAnglex[matCount] = frictionAng;
00254   peakShearStrainx[matCount] = peakShearStra;
00255   refPressurex[matCount] = -refPress;  //compression is negative
00256   cohesionx[matCount] = cohesi;
00257   pressDependCoeffx[matCount] = pressDependCoe;
00258   numOfSurfacesx[matCount] = numberOfYieldSurf;
00259   rhox[matCount] = r;
00260   phaseTransfAnglex[matCount] = phaseTransformAng;
00261   contractParam1x[matCount] = contractionParam1;
00262   dilateParam1x[matCount] = dilationParam1;
00263   dilateParam2x[matCount] = dilationParam2;
00264   volLimit1x[matCount] = volLim1;
00265   volLimit2x[matCount] = volLim2;
00266   volLimit3x[matCount] = volLim3;
00267   liquefyParam1x[matCount] = liquefactionParam1;
00268   liquefyParam2x[matCount] = liquefactionParam2;
00269   liquefyParam4x[matCount] = liquefactionParam4;
00270   einitx[matCount] = ei;
00271   Hvx[matCount] = hv;
00272   Pvx[matCount] = pv;
00273 
00274   matN = matCount;
00275   matCount ++;
00276   pAtm = atm;
00277 
00278   int numOfSurfaces = numOfSurfacesx[matN];
00279   initPress = refPressurex[matN];
00280 
00281   e2p = committedActiveSurf = activeSurfaceNum = 0; 
00282   onPPZCommitted = onPPZ = -1 ; 
00283   PPZSizeCommitted = PPZSize = 0.;
00284   pressureDCommitted = pressureD = modulusFactor = 0.;
00285   cumuDilateStrainOctaCommitted    = cumuDilateStrainOcta = 0.;
00286   maxCumuDilateStrainOctaCommitted = maxCumuDilateStrainOcta = 0.;
00287   cumuTranslateStrainOctaCommitted = cumuTranslateStrainOcta = 0.;
00288   prePPZStrainOctaCommitted = prePPZStrainOcta = 0.;
00289   oppoPrePPZStrainOctaCommitted = oppoPrePPZStrainOcta = 0.;
00290   maxPress = 0.;
00291 
00292   theSurfaces = new MultiYieldSurface[numOfSurfaces+1]; //first surface not used
00293   committedSurfaces = new MultiYieldSurface[numOfSurfaces+1]; 
00294 
00295   setUpSurfaces(gredu);  // residualPress and stressRatioPT are calculated inside.
00296 }
00297    
00298 
00299 PressureDependMultiYield::PressureDependMultiYield () 
00300  : NDMaterial(0,ND_TAG_PressureDependMultiYield), 
00301    currentStress(), trialStress(), currentStrain(), 
00302   strainRate(), reversalStress(), PPZPivot(), 
00303   PPZCenter(), lockStress(), reversalStressCommitted(), 
00304   PPZPivotCommitted(), PPZCenterCommitted(),
00305   lockStressCommitted(), theSurfaces(0), committedSurfaces(0)
00306 {
00307   //does nothing
00308 }
00309 
00310 PressureDependMultiYield::PressureDependMultiYield (const PressureDependMultiYield & a)
00311  : NDMaterial(a.getTag(),ND_TAG_PressureDependMultiYield), 
00312    currentStress(a.currentStress), trialStress(a.trialStress), 
00313   currentStrain(a.currentStrain), strainRate(a.strainRate), 
00314    reversalStress(a.reversalStress), PPZPivot(a.PPZPivot), PPZCenter(a.PPZCenter), 
00315   lockStress(a.lockStress), reversalStressCommitted(a.reversalStressCommitted), 
00316   PPZPivotCommitted(a.PPZPivotCommitted), 
00317   PPZCenterCommitted(a.PPZCenterCommitted),
00318   lockStressCommitted(a.lockStressCommitted)
00319 {
00320   matN = a.matN;
00321 
00322   int numOfSurfaces = numOfSurfacesx[matN];
00323 
00324   e2p = a.e2p;
00325   strainPTOcta = a.strainPTOcta;
00326   modulusFactor = a.modulusFactor;
00327   activeSurfaceNum = a.activeSurfaceNum;
00328   committedActiveSurf = a.committedActiveSurf;
00329   pressureDCommitted     = a.pressureDCommitted;
00330   onPPZCommitted = a.onPPZCommitted; 
00331   PPZSizeCommitted      = a.PPZSizeCommitted;
00332   cumuDilateStrainOctaCommitted    = a.cumuDilateStrainOctaCommitted;
00333   maxCumuDilateStrainOctaCommitted = a.maxCumuDilateStrainOctaCommitted;
00334   cumuTranslateStrainOctaCommitted = a.cumuTranslateStrainOctaCommitted;
00335   prePPZStrainOctaCommitted        = a.prePPZStrainOctaCommitted;
00336   oppoPrePPZStrainOctaCommitted    = a.oppoPrePPZStrainOctaCommitted;
00337   pressureD     = a.pressureD;
00338   onPPZ = a.onPPZ; 
00339   PPZSize      = a.PPZSize;
00340   cumuDilateStrainOcta    = a.cumuDilateStrainOcta;
00341   maxCumuDilateStrainOcta = a.maxCumuDilateStrainOcta;
00342   cumuTranslateStrainOcta = a.cumuTranslateStrainOcta;
00343   prePPZStrainOcta        = a.prePPZStrainOcta;
00344   oppoPrePPZStrainOcta    = a.oppoPrePPZStrainOcta;
00345   initPress = a.initPress;
00346   maxPress = a.maxPress;
00347 
00348   theSurfaces = new MultiYieldSurface[numOfSurfaces+1];  //first surface not used
00349   committedSurfaces = new MultiYieldSurface[numOfSurfaces+1];  
00350   for(int i=1; i<=numOfSurfaces; i++) {
00351     committedSurfaces[i] = a.committedSurfaces[i];  
00352     theSurfaces[i] = a.theSurfaces[i];  
00353   }
00354 }
00355 
00356 PressureDependMultiYield::~PressureDependMultiYield ()
00357 {
00358   if (theSurfaces != 0) delete [] theSurfaces;
00359   if (committedSurfaces != 0) delete [] committedSurfaces;
00360 }
00361 
00362 void 
00363 PressureDependMultiYield::elast2Plast(void)
00364 {
00365   int loadStage = loadStagex[matN];
00366   int numOfSurfaces = numOfSurfacesx[matN];
00367 
00368   if (loadStage != 1 || e2p == 1) return;
00369   e2p = 1;
00370 
00371   if (currentStress.volume() > 0.) {
00372     //opserr << "WARNING:PressureDependMultiYield::elast2Plast(): material in tension." << endln;
00373     currentStress.setData(currentStress.deviator(),0);
00374   }
00375 
00376   // Active surface is 0, return
00377   if (currentStress.deviatorLength() == 0.) return;
00378 
00379   // Find active surface
00380   while (yieldFunc(currentStress, committedSurfaces, ++committedActiveSurf) > 0) {
00381     if (committedActiveSurf == numOfSurfaces) {
00382       //opserr <<"WARNING:PressureDependMultiYield::elast2Plast(): stress out of failure surface"<<endln;
00383       deviatorScaling(currentStress, committedSurfaces, numOfSurfaces);
00384       initSurfaceUpdate();
00385       return;
00386     }
00387   } 
00388 
00389   committedActiveSurf--;
00390   initSurfaceUpdate();
00391 }
00392 
00393  
00394 int 
00395 PressureDependMultiYield::setTrialStrain (const Vector &strain)
00396 {
00397   int ndm = ndmx[matN];
00398 
00399   if (ndm==3 && strain.Size()==6) 
00400     workV6 = strain;
00401   else if (ndm==2 && strain.Size()==3) {
00402     workV6[0] = strain[0];
00403     workV6[1] = strain[1];
00404     workV6[2] = 0.0;
00405     workV6[3] = strain[2];
00406     workV6[4] = 0.0;
00407     workV6[5] = 0.0;
00408   }
00409   else {
00410     opserr << "Fatal:PressureDependMultiYield:: Material dimension is: " << ndm << endln;
00411     opserr << "But strain vector size is: " << strain.Size() << endln;
00412    exit(-1);
00413   }
00414 
00415   //strainRate.setData(workV6-currentStrain.t2Vector(1),1);
00416   workV6 -= currentStrain.t2Vector(1);
00417   strainRate.setData(workV6, 1);
00418 
00419   return 0;
00420 }
00421 
00422 int 
00423 PressureDependMultiYield::setTrialStrain (const Vector &strain, const Vector &rate)
00424 {
00425   return setTrialStrain (strain);
00426 }
00427 
00428 int 
00429 PressureDependMultiYield::setTrialStrainIncr (const Vector &strain)
00430 {
00431   int ndm = ndmx[matN];
00432 
00433   if (ndm==3 && strain.Size()==6) 
00434     workV6 = strain;
00435   else if (ndm==2 && strain.Size()==3) {
00436     workV6[0] = strain[0];
00437     workV6[1] = strain[1];
00438     workV6[2] = 0.0;
00439     workV6[3] = strain[2];
00440     workV6[4] = 0.0;
00441     workV6[5] = 0.0;
00442   }
00443   else {
00444     opserr << "Fatal:PressureDependMultiYield:: Material dimension is: " << ndm << endln;
00445     opserr << "But strain vector size is: " << strain.Size() << endln;
00446    exit(-1);
00447   }
00448 
00449   strainRate.setData(workV6,1);
00450   return 0;
00451 }
00452 
00453 int 
00454 PressureDependMultiYield::setTrialStrainIncr (const Vector &strain, const Vector &rate)
00455 {
00456   return setTrialStrainIncr(strain);
00457 }
00458 
00459 const Matrix & 
00460 PressureDependMultiYield::getTangent (void)
00461 {
00462   int loadStage = loadStagex[matN];
00463   double refShearModulus = refShearModulusx[matN];
00464   double refBulkModulus = refBulkModulusx[matN];
00465   double pressDependCoeff = pressDependCoeffx[matN];
00466   double refPressure = refPressurex[matN];
00467   double residualPress = residualPressx[matN];
00468   int ndm = ndmx[matN];
00469 
00470   if (loadStage == 1 && e2p == 0) elast2Plast();
00471   if (loadStage==2 && initPress==refPressure) 
00472           initPress = currentStress.volume();
00473 
00474   if (loadStage==0 || loadStage==2) {  //linear elastic
00475         double factor;
00476         if (loadStage==0) factor = 1.0;
00477         else {
00478                 factor = (initPress-residualPress)/(refPressure-residualPress);
00479                 if (factor <= 1.e-10) factor = 1.e-10;
00480                 else factor = pow(factor, pressDependCoeff);
00481                 factor = (1.e-10>factor) ? 1.e-10 : factor;
00482         }
00483     for (int i=0;i<6;i++) 
00484       for (int j=0;j<6;j++) {
00485             theTangent(i,j) = 0.;
00486         if (i==j) theTangent(i,j) += refShearModulus*factor;
00487         if (i<3 && j<3 && i==j) theTangent(i,j) += refShearModulus*factor;
00488         if (i<3 && j<3) theTangent(i,j) += (refBulkModulus - 2.*refShearModulus/3.)*factor;
00489       }
00490   }
00491   else {
00492     double coeff1, coeff2, coeff3, coeff4;
00493     double factor = getModulusFactor(currentStress);
00494     double shearModulus = factor*refShearModulus;
00495     double bulkModulus = factor*refBulkModulus;         
00496 
00497         // volumetric plasticity
00498         if (Hvx[matN] != 0. && trialStress.volume()<=maxPress && strainRate.volume()<0.) {
00499           double tp = fabs(trialStress.volume() - residualPress);
00500       bulkModulus = (bulkModulus*Hvx[matN]*pow(tp,Pvx[matN]))/(bulkModulus+Hvx[matN]*pow(tp,Pvx[matN]));
00501         }
00502         
00503     if (loadStage!=0 && committedActiveSurf > 0) {
00504       getSurfaceNormal(currentStress, workT2V);
00505       workV6 = workT2V.deviator();
00506       double volume = workT2V.volume();
00507       double Ho = 9.*bulkModulus*volume*volume + 2.*shearModulus*(workV6 && workV6);
00508       double plastModul = factor*committedSurfaces[committedActiveSurf].modulus();
00509       coeff1 = 9.*bulkModulus*bulkModulus*volume*volume/(Ho+plastModul);
00510       coeff2 = 4.*shearModulus*shearModulus/(Ho+plastModul); 
00511       /* non-symmetric stiffness
00512       getSurfaceNormal(currentStress, workT2V);
00513       workV6 = workT2V.deviator();
00514       double qq = workT2V.volume();
00515                         double pp=getPlasticPotential(currentStress, workT2V); 
00516       double Ho = 9.*bulkModulus*pp*qq + 2.*shearModulus*(workV6 && workV6);
00517       double plastModul = factor*committedSurfaces[committedActiveSurf].modulus();
00518       coeff1 = 9.*bulkModulus*bulkModulus*pp*qq/(Ho+plastModul);
00519       coeff2 = 4.*shearModulus*shearModulus/(Ho+plastModul);
00520                         coeff3 = 6.*shearModulus*pp/(Ho+plastModul);
00521                         coeff4 = 6.*shearModulus*qq/(Ho+plastModul);*/
00522 
00523           }
00524     else {
00525       coeff1 = coeff2 = coeff3 = coeff4 = 0.;
00526       workV6.Zero();
00527     }
00528 
00529     for (int i=0;i<6;i++) 
00530       for (int j=0;j<6;j++) {
00531               theTangent(i,j) = - coeff2*workV6[i]*workV6[j];
00532         if (i==j) theTangent(i,j) += shearModulus;
00533         if (i<3 && j<3 && i==j) theTangent(i,j) += shearModulus;
00534               if (i<3 && j<3) theTangent(i,j) += (bulkModulus - 2.*shearModulus/3. - coeff1);
00535 /* non-symmetric stiffness
00536                                 if (i<3) theTangent(i,j) -= coeff3 * workV6[j];
00537                                 if (j<3) theTangent(i,j) -= coeff4 * workV6[i];*/
00538       }
00539   }  if (ndm==3) 
00540     return theTangent;
00541   else {
00542     static Matrix workM(3,3);
00543     workM(0,0) = theTangent(0,0);
00544     workM(0,1) = theTangent(0,1);
00545     workM(0,2) = 0.;
00546 
00547     workM(1,0) = theTangent(1,0);
00548     workM(1,1) = theTangent(1,1);
00549     workM(1,2) = 0.; 
00550 
00551     workM(2,0) = 0.; 
00552     workM(2,1) = 0.; 
00553     workM(2,2) = theTangent(3,3);
00554 
00555     /* non-symmetric stiffness
00556        workM(0,2) = theTangent(0,3);
00557        workM(1,2) = theTangent(1,3);
00558        workM(2,0) = theTangent(3,0);
00559        workM(2,1) = theTangent(3,1);*/
00560 
00561     return workM;
00562   }
00563 }
00564 
00565 const Matrix & 
00566 PressureDependMultiYield::getInitialTangent (void)
00567 {
00568   int loadStage = loadStagex[matN];
00569   double refShearModulus = refShearModulusx[matN];
00570   double refBulkModulus = refBulkModulusx[matN];
00571   double pressDependCoeff = pressDependCoeffx[matN];
00572   double refPressure = refPressurex[matN];
00573   double residualPress = residualPressx[matN];
00574   int ndm = ndmx[matN];
00575 
00576   if (loadStage==2 && initPress==refPressure) 
00577           initPress = currentStress.volume();
00578   double factor; 
00579   if (loadStage==0)
00580           factor = 1.;
00581   else if (loadStage==2) {
00582                 factor = (initPress-residualPress)/(refPressure-residualPress);
00583                 if (factor <= 1.e-10) factor = 1.e-10;
00584                 else factor = pow(factor, pressDependCoeff);
00585                 factor = (1.e-10>factor) ? 1.e-10 : factor;
00586   }
00587   else if (loadStage==1)
00588           factor = getModulusFactor(currentStress);
00589   
00590   for (int i=0;i<6;i++) 
00591     for (int j=0;j<6;j++) {
00592       theTangent(i,j) = 0.;
00593       if (i==j) theTangent(i,j) += refShearModulus*factor;
00594       if (i<3 && j<3 && i==j) theTangent(i,j) += refShearModulus*factor;
00595       if (i<3 && j<3) theTangent(i,j) += (refBulkModulus - 2.*refShearModulus/3.)*factor;
00596     }
00597 
00598   if (ndm==3) 
00599     return theTangent;
00600   else {
00601     static Matrix workM(3,3);
00602     workM(0,0) = theTangent(0,0);
00603     workM(0,1) = theTangent(0,1);
00604     workM(0,2) = 0.;
00605 
00606     workM(1,0) = theTangent(1,0);
00607     workM(1,1) = theTangent(1,1);
00608     workM(1,2) = 0.; 
00609 
00610     workM(2,0) = 0.; 
00611     workM(2,1) = 0.; 
00612     workM(2,2) = theTangent(3,3);
00613 
00614 /* non-symmetric stiffness
00615     workM(0,2) = theTangent(0,3);
00616     workM(1,2) = theTangent(1,3);
00617     workM(2,0) = theTangent(3,0);
00618     workM(2,1) = theTangent(3,1);*/
00619     return workM;
00620   }
00621 }
00622 
00623 const Vector &
00624 PressureDependMultiYield::getStress (void)
00625 {
00626   int loadStage = loadStagex[matN];
00627   int numOfSurfaces = numOfSurfacesx[matN];
00628   int ndm = ndmx[matN];
00629 
00630   int i, is;
00631   if (loadStage == 1 && e2p == 0) 
00632     elast2Plast();
00633 
00634   if (loadStage!=1) {  //linear elastic
00635     //trialStrain.setData(currentStrain.t2Vector() + strainRate.t2Vector());
00636     getTangent();
00637     workV6 = currentStress.t2Vector();
00638         workV6.addMatrixVector(1.0, theTangent, strainRate.t2Vector(1), 1.0);
00639     trialStress.setData(workV6);
00640   }
00641   else {
00642     for (i=1; i<=numOfSurfaces; i++) theSurfaces[i] = committedSurfaces[i];
00643     activeSurfaceNum = committedActiveSurf;
00644     pressureD = pressureDCommitted;
00645     reversalStress = reversalStressCommitted;
00646     onPPZ = onPPZCommitted;
00647     PPZSize = PPZSizeCommitted;
00648     cumuDilateStrainOcta = cumuDilateStrainOctaCommitted;
00649     maxCumuDilateStrainOcta = maxCumuDilateStrainOctaCommitted;
00650     cumuTranslateStrainOcta = cumuTranslateStrainOctaCommitted;
00651     prePPZStrainOcta = prePPZStrainOctaCommitted;
00652     oppoPrePPZStrainOcta = oppoPrePPZStrainOctaCommitted;
00653     PPZPivot = PPZPivotCommitted;
00654     PPZCenter = PPZCenterCommitted;
00655     lockStress = lockStressCommitted;
00656 
00657     subStrainRate = strainRate;
00658     setTrialStress(currentStress);
00659     if (activeSurfaceNum>0 && isLoadReversal(currentStress)) {
00660       updateInnerSurface();
00661       activeSurfaceNum = 0;
00662     }
00663 
00664     if (activeSurfaceNum==0 && !isCrossingNextSurface()) {
00665             workV6 = currentStrain.t2Vector();
00666             workV6.addVector(1.0, strainRate.t2Vector(), 1.0);
00667             trialStrain.setData(workV6);
00668                 }
00669                 else {
00670       int numSubIncre = setSubStrainRate();
00671 
00672       for (i=0; i<numSubIncre; i++) {
00673 //      trialStrain.setData(currentStrain.t2Vector() 
00674 //                           + subStrainRate.t2Vector()*(i+1));
00675               workV6 = currentStrain.t2Vector();
00676               workV6.addVector(1.0, subStrainRate.t2Vector(), (i+1));
00677               trialStrain.setData(workV6);
00678 
00679                 if (i==0)  {
00680                           setTrialStress(currentStress);
00681               is = isLoadReversal(currentStress);
00682                 }  
00683                 else {
00684                         workT2V.setData(trialStress.t2Vector());
00685                         setTrialStress(trialStress); 
00686             is = isLoadReversal(workT2V);
00687                 }
00688         if (activeSurfaceNum>0 && is) {
00689           updateInnerSurface();
00690           activeSurfaceNum = 0;
00691                 }
00692         if (activeSurfaceNum==0 && !isCrossingNextSurface()) continue;
00693         if (activeSurfaceNum==0) activeSurfaceNum++;
00694         int lock = stressCorrection(0);
00695         if(lock==0) updateActiveSurface();
00696                 //opserr<<i<<" "<<activeSurfaceNum<<" "<<is<<" "<<subStrainRate.t2Vector()[3]<<endln;
00697                 }
00698     }
00699   }
00700 
00701   if (ndm==3)
00702     return trialStress.t2Vector();
00703   else {
00704                 static Vector workV(3);
00705     workV[0] = trialStress.t2Vector()[0];
00706     workV[1] = trialStress.t2Vector()[1];
00707     workV[2] = trialStress.t2Vector()[3];
00708     return workV;
00709   }
00710 }
00711 
00712 const Vector & 
00713 PressureDependMultiYield::getStrain (void)
00714 {
00715   return getCommittedStrain();
00716 }
00717 
00718 int 
00719 PressureDependMultiYield::commitState (void)
00720 {
00721   int loadStage = loadStagex[matN];
00722   int numOfSurfaces = numOfSurfacesx[matN];
00723 
00724   currentStress = trialStress;
00725   //currentStrain = T2Vector(currentStrain.t2Vector() + strainRate.t2Vector());
00726   workV6 = currentStrain.t2Vector();
00727   workV6 += strainRate.t2Vector();
00728   currentStrain.setData(workV6);
00729 
00730   workV6.Zero();
00731   strainRate.setData(workV6);
00732 
00733   if (loadStage==1) {
00734     committedActiveSurf = activeSurfaceNum;
00735     for (int i=1; i<=numOfSurfaces; i++) committedSurfaces[i] = theSurfaces[i];
00736     pressureDCommitted = pressureD;
00737     reversalStressCommitted = reversalStress;
00738     onPPZCommitted = onPPZ;
00739     PPZSizeCommitted = PPZSize;
00740     cumuDilateStrainOctaCommitted = cumuDilateStrainOcta;
00741     maxCumuDilateStrainOctaCommitted = maxCumuDilateStrainOcta;
00742     cumuTranslateStrainOctaCommitted = cumuTranslateStrainOcta;
00743     prePPZStrainOctaCommitted = prePPZStrainOcta;
00744     oppoPrePPZStrainOctaCommitted = oppoPrePPZStrainOcta;
00745     PPZPivotCommitted = PPZPivot;
00746     PPZCenterCommitted = PPZCenter;
00747     lockStressCommitted = lockStress;
00748         if (currentStress.volume() < maxPress) maxPress = currentStress.volume();
00749   }
00750 
00751   return 0;
00752 }
00753 
00754 int 
00755 PressureDependMultiYield::revertToLastCommit (void)
00756 {
00757   return 0;
00758 }
00759 
00760 NDMaterial * 
00761 PressureDependMultiYield::getCopy (void)
00762 {
00763   PressureDependMultiYield * copy = new PressureDependMultiYield(*this);
00764   return copy;
00765 }
00766 
00767 NDMaterial * 
00768 PressureDependMultiYield::getCopy (const char *code)
00769 {
00770   if (strcmp(code,"PlaneStrain") == 0 ||
00771       strcmp(code,"ThreeDimensional") == 0) {
00772     PressureDependMultiYield * copy = new PressureDependMultiYield(*this);
00773     return copy;
00774   }
00775 
00776   return 0;
00777 }
00778 
00779 const char * 
00780 PressureDependMultiYield::getType (void) const
00781 {
00782   int ndm = ndmx[matN];
00783 
00784   return (ndm == 2) ? "PlaneStrain" : "ThreeDimensional";
00785 }
00786 
00787 int 
00788 PressureDependMultiYield::getOrder (void) const
00789 {
00790   int ndm = ndmx[matN];
00791 
00792   return (ndm == 2) ? 3 : 6;
00793 }
00794 
00795 int 
00796 PressureDependMultiYield::updateParameter(int responseID, Information &info)
00797 {
00798         if (responseID<10) 
00799                 loadStagex[matN] = responseID;
00800 
00801         else {
00802                 if (responseID==10) refShearModulusx[matN]=info.theDouble;
00803                 if (responseID==11) refBulkModulusx[matN]=info.theDouble;
00804         }
00805 
00806   return 0;
00807 }
00808 
00809 int 
00810 PressureDependMultiYield::sendSelf(int commitTag, Channel &theChannel)
00811 {
00812     int loadStage = loadStagex[matN];
00813     int ndm = ndmx[matN];
00814         double rho = rhox[matN];
00815     double residualPress = residualPressx[matN];
00816     int numOfSurfaces = numOfSurfacesx[matN];
00817     double refPressure = refPressurex[matN];
00818     double pressDependCoeff =pressDependCoeffx[matN];
00819     double refShearModulus = refShearModulusx[matN];
00820         double refBulkModulus = refBulkModulusx[matN];
00821     double frictionAngle = frictionAnglex[matN];
00822         double cohesion = cohesionx[matN];
00823     double peakShearStrain = peakShearStrainx[matN];
00824     double phaseTransfAngle = phaseTransfAnglex[matN];
00825         double stressRatioPT = stressRatioPTx[matN];
00826         double contractParam1 = contractParam1x[matN];
00827     double dilateParam1 = dilateParam1x[matN];
00828     double dilateParam2 = dilateParam2x[matN];
00829         double liquefyParam1 = liquefyParam1x[matN];
00830         double liquefyParam2 = liquefyParam2x[matN];
00831         double liquefyParam4 = liquefyParam4x[matN];
00832         double einit = einitx[matN];
00833         double volLimit1 = volLimit1x[matN];
00834         double volLimit2 = volLimit2x[matN];
00835         double volLimit3 = volLimit3x[matN];
00836 
00837   int i, res = 0;
00838 
00839   static ID idData(5);
00840   idData(0) = this->getTag();
00841   idData(1) = numOfSurfaces;
00842   idData(2) = loadStage;
00843   idData(3) = ndm;
00844   idData(4) = matN;
00845 
00846   res += theChannel.sendID(this->getDbTag(), commitTag, idData);
00847   if (res < 0) {
00848     opserr << "PressureDependMultiYield::sendSelf -- could not send ID\n";
00849     return res;
00850   }
00851 
00852   Vector data(70+numOfSurfaces*8);
00853   data(0) = rho;
00854   data(1) = einit;
00855   data(2) = refShearModulus;
00856   data(3) = refBulkModulus;
00857   data(4) = frictionAngle;
00858   data(5) = peakShearStrain;
00859   data(6) = refPressure;
00860   data(7) = cohesion;
00861   data(8) = pressDependCoeff;
00862   data(9) = phaseTransfAngle;
00863   data(10) = contractParam1;
00864   data(11) = dilateParam1;
00865   data(12) = dilateParam2;
00866   data(13) = volLimit1;
00867   data(14) = volLimit2;
00868   data(15) = volLimit3;
00869   data(16) = pAtm;
00870   data(17) = liquefyParam1;
00871   data(18) = liquefyParam2;
00872   data(19) = liquefyParam4;
00873   data(20) = residualPress;
00874   data(21) = stressRatioPT;   
00875   data(22) = e2p; 
00876   data(23) = committedActiveSurf;
00877   data(24) = strainPTOcta;
00878   data(25) = pressureDCommitted;
00879   data(26) = onPPZCommitted;
00880   data(27) = PPZSizeCommitted;
00881   data(28) = cumuDilateStrainOctaCommitted;
00882   data(29) = maxCumuDilateStrainOctaCommitted;
00883   data(30) = cumuTranslateStrainOctaCommitted;
00884   data(31) = prePPZStrainOctaCommitted;
00885   data(32) = oppoPrePPZStrainOctaCommitted;
00886   data(69) = initPress;
00887 
00888   workV6 = currentStress.t2Vector();
00889   for(i = 0; i < 6; i++) data(i+33) = workV6[i];
00890 
00891   workV6 = currentStrain.t2Vector();
00892   for(i = 0; i < 6; i++) data(i+39) = workV6[i];
00893           
00894   workV6 = PPZPivotCommitted.t2Vector();
00895   for(i = 0; i < 6; i++) data(i+45) = workV6[i];
00896 
00897   workV6 = PPZCenterCommitted.t2Vector();
00898   for(i = 0; i < 6; i++) data(i+51) = workV6[i];
00899 
00900   workV6 = lockStressCommitted.t2Vector();
00901   for(i = 0; i < 6; i++) data(i+57) = workV6[i];
00902         
00903   workV6 = reversalStressCommitted.t2Vector();
00904   for(i = 0; i < 6; i++) data(i+63) = workV6[i];
00905         
00906   for(i = 0; i < numOfSurfaces; i++) {
00907     int k = 70 + i*8;
00908     data(k) = committedSurfaces[i+1].size();
00909     data(k+1) = committedSurfaces[i+1].modulus();
00910     workV6 = committedSurfaces[i+1].center();
00911     data(k+2) = workV6(0);
00912     data(k+3) = workV6(1);
00913     data(k+4) = workV6(2);
00914     data(k+5) = workV6(3);
00915     data(k+6) = workV6(4);
00916     data(k+7) = workV6(5);
00917   }
00918 
00919   res += theChannel.sendVector(this->getDbTag(), commitTag, data);
00920   if (res < 0) {
00921     opserr << "PressureDependMultiYield::sendSelf -- could not send Vector\n";
00922     return res;
00923   }
00924 
00925   return res;
00926 }
00927 
00928 int 
00929 PressureDependMultiYield::recvSelf(int commitTag, Channel &theChannel, 
00930                                        FEM_ObjectBroker &theBroker)    
00931 {
00932   int i, res = 0;
00933 
00934   static ID idData(5);
00935   res += theChannel.recvID(this->getDbTag(), commitTag, idData);
00936   if (res < 0) {
00937     opserr << "PressureDependMultiYield::recvelf -- could not recv ID\n";
00938                             
00939     return res;
00940   }
00941 
00942   this->setTag((int)idData(0));
00943   int numOfSurfaces = idData(1);
00944   int loadStage = idData(2);
00945   int ndm = idData(3);
00946   matN = idData(4);
00947 
00948   Vector data(70+idData(1)*8);
00949   res += theChannel.recvVector(this->getDbTag(), commitTag, data);
00950   if (res < 0) {
00951     opserr << "PressureDependMultiYield::recvSelf -- could not recv Vector\n";
00952     return res;
00953   }
00954 
00955   double rho = data(0);
00956   double einit = data(1);
00957   double refShearModulus = data(2);
00958   double refBulkModulus = data(3);
00959   double frictionAngle = data(4);
00960   double peakShearStrain = data(5);
00961   double refPressure = data(6);
00962   double cohesion = data(7);
00963   double pressDependCoeff = data(8);
00964   double phaseTransfAngle = data(9);
00965   double contractParam1 = data(10);
00966   double dilateParam1 = data(11);
00967   double dilateParam2 = data(12);
00968   double volLimit1 = data(13);
00969   double volLimit2 = data(14);
00970   double volLimit3 = data(15);
00971   pAtm = data(16);
00972   double liquefyParam1 = data(17);
00973   double liquefyParam2 = data(18);
00974   double liquefyParam4 = data(19);
00975   double residualPress = data(20);
00976   double stressRatioPT = data(21);   
00977   e2p = data(22); 
00978   committedActiveSurf = data(23);
00979   strainPTOcta = data(24);
00980   pressureDCommitted = data(25);
00981   onPPZCommitted = data(26);
00982   PPZSizeCommitted = data(27);
00983   cumuDilateStrainOctaCommitted = data(28);
00984   maxCumuDilateStrainOctaCommitted = data(29);
00985   cumuTranslateStrainOctaCommitted = data(30);
00986   prePPZStrainOctaCommitted = data(31);
00987   oppoPrePPZStrainOctaCommitted = data(32);
00988   initPress = data(69);
00989 
00990   for(i = 0; i < 6; i++) workV6[i] = data(i+33);
00991   currentStress.setData(workV6);
00992 
00993   for(i = 0; i < 6; i++) workV6[i] = data(i+39);
00994   currentStrain.setData(workV6);
00995   
00996   for(i = 0; i < 6; i++) workV6[i] = data(i+45);
00997   PPZPivotCommitted.setData(workV6);
00998 
00999   for(i = 0; i < 6; i++) workV6[i] = data(i+51);
01000   PPZCenterCommitted.setData(workV6);
01001   
01002   for(i = 0; i < 6; i++) workV6[i] = data(i+57);
01003   lockStressCommitted.setData(workV6);
01004   
01005   for(i = 0; i < 6; i++) workV6[i] = data(i+63);
01006   reversalStressCommitted.setData(workV6);  if (committedSurfaces != 0) {
01007       delete [] committedSurfaces;
01008       delete [] theSurfaces;
01009   }
01010 
01011   theSurfaces = new MultiYieldSurface[numOfSurfaces+1]; //first surface not used
01012   committedSurfaces = new MultiYieldSurface[numOfSurfaces+1]; 
01013   for (i=1; i<=numOfSurfaces; i++) 
01014       committedSurfaces[i] = MultiYieldSurface();    
01015 
01016   for(i = 0; i < numOfSurfaces; i++) {
01017     int k = 70 + i*8;
01018     workV6(0) = data(k+2);
01019     workV6(1) = data(k+3);
01020     workV6(2) = data(k+4);
01021     workV6(3) = data(k+5);
01022     workV6(4) = data(k+6);
01023     workV6(5) = data(k+7);
01024     committedSurfaces[i+1].setData(workV6, data(k), data(k+1));
01025   }
01026 
01027     loadStagex[matN] = loadStage;
01028     ndmx[matN] = ndm;
01029         rhox[matN] = rho;
01030     residualPressx[matN] = residualPress;
01031     numOfSurfacesx[matN] = numOfSurfaces;
01032     refPressurex[matN] = refPressure;
01033     pressDependCoeffx[matN] =pressDependCoeff;
01034     refShearModulusx[matN] = refShearModulus;
01035         refBulkModulusx[matN] = refBulkModulus;
01036     frictionAnglex[matN] = frictionAngle;
01037         cohesionx[matN] = cohesion;
01038     peakShearStrainx[matN] = peakShearStrain;
01039     phaseTransfAnglex[matN] = phaseTransfAngle;
01040         stressRatioPTx[matN] = stressRatioPT;
01041         contractParam1x[matN] = contractParam1;
01042     dilateParam1x[matN] = dilateParam1;
01043     dilateParam2x[matN] = dilateParam2;
01044         liquefyParam1x[matN] = liquefyParam1;
01045         liquefyParam2x[matN] = liquefyParam2;
01046         liquefyParam4x[matN] = liquefyParam4;
01047         einitx[matN] = einit;
01048         volLimit1x[matN] = volLimit1;
01049         volLimit2x[matN] = volLimit2;
01050         volLimit3x[matN] = volLimit3;
01051 
01052   return res;
01053 }
01054 
01055 Response*
01056 PressureDependMultiYield::setResponse (const char **argv, int argc, Information &matInfo, OPS_Stream &output)
01057 {
01058   if (strcmp(argv[0],"stress") == 0 || strcmp(argv[0],"stresses") == 0)
01059                 return new MaterialResponse(this, 1, this->getCommittedStress());
01060 
01061   else if (strcmp(argv[0],"strain") == 0 || strcmp(argv[0],"strains") == 0)
01062                 return new MaterialResponse(this, 2, this->getCommittedStrain());
01063     
01064         else if (strcmp(argv[0],"tangent") == 0)
01065                 return new MaterialResponse(this, 3, this->getTangent());
01066     
01067         else if (strcmp(argv[0],"backbone") == 0) {
01068             int numOfSurfaces = numOfSurfacesx[matN];
01069         static Matrix curv(numOfSurfaces+1,(argc-1)*2);
01070                 for (int i=1; i<argc; i++)
01071                         curv(0,(i-1)*2) = atoi(argv[i]);
01072                 return new MaterialResponse(this, 4, curv);
01073   }
01074         else
01075                 return 0;
01076 }
01077 
01078 void 
01079 PressureDependMultiYield::getBackbone (Matrix & bb)
01080 {
01081   double residualPress = residualPressx[matN];
01082   double refPressure = refPressurex[matN];
01083   double pressDependCoeff =pressDependCoeffx[matN];
01084   double refShearModulus = refShearModulusx[matN];
01085   int numOfSurfaces = numOfSurfacesx[matN];
01086 
01087   double vol, conHeig, scale, factor, shearModulus, stress1, 
01088                      stress2, strain1, strain2, plastModulus, elast_plast, gre;
01089 
01090         for (int k=0; k<bb.noCols()/2; k++) {
01091                 vol = bb(0,k*2);
01092                 if (vol<=0.) {
01093                         opserr <<k<< "\nNDMaterial " <<this->getTag()
01094                           <<": invalid confinement for backbone recorder, " << vol << endln;
01095                         continue;
01096                 }
01097                 conHeig = vol + residualPress;
01098                 scale = -conHeig / (refPressure-residualPress);
01099                 factor = pow(scale, pressDependCoeff); 
01100                 shearModulus = factor*refShearModulus;
01101 
01102                 for (int i=1; i<=numOfSurfaces; i++) {
01103                         if (i==1) {
01104                                 stress2 = theSurfaces[i].size()*conHeig/sqrt(3.0);
01105                                 strain2 = stress2/shearModulus;
01106                                 bb(1,k*2) = strain2; bb(1,k*2+1) = shearModulus;
01107                         } else {
01108                                 stress1 = stress2; strain1 = strain2;
01109                                 plastModulus = factor*theSurfaces[i-1].modulus();
01110                                 elast_plast = 2*shearModulus*plastModulus/(2*shearModulus+plastModulus);
01111                                 stress2 = theSurfaces[i].size()*conHeig/sqrt(3.0);
01112                           strain2 = 2*(stress2-stress1)/elast_plast + strain1;
01113                                 gre = stress2/strain2;
01114         bb(i,k*2) = strain2; bb(i,k*2+1) = gre;
01115                         }
01116                 }
01117         }
01118 
01119 }
01120 
01121 int 
01122 PressureDependMultiYield::getResponse (int responseID, Information &matInfo)
01123 {
01124   switch (responseID) {
01125   case -1:
01126     return -1;
01127   case 1:
01128     if (matInfo.theVector != 0)
01129       *(matInfo.theVector) = getCommittedStress();
01130     return 0;
01131   case 2:
01132     if (matInfo.theVector != 0)
01133       *(matInfo.theVector) = getCommittedStrain();
01134     return 0;
01135   case 3:
01136     if (matInfo.theMatrix != 0)
01137       *(matInfo.theMatrix) = getTangent();
01138     return 0;
01139   case 4:
01140     if (matInfo.theMatrix != 0) 
01141       getBackbone(*(matInfo.theMatrix));
01142     return 0;
01143   default:
01144     return -1;
01145   }
01146 }
01147 
01148 void 
01149 PressureDependMultiYield::Print(OPS_Stream &s, int flag )
01150 {
01151   s << "PressureDependMultiYield" << endln;
01152 }
01153 
01154 const Vector & 
01155 PressureDependMultiYield::getCommittedStress (void)
01156 {
01157         int ndm = ndmx[matN];
01158         int numOfSurfaces = numOfSurfacesx[matN];
01159     double residualPress = residualPressx[matN];
01160 
01161         double scale = currentStress.deviatorRatio(residualPress)/committedSurfaces[numOfSurfaces].size();
01162         if (loadStagex[matN] != 1) scale = 0.;
01163   if (ndm==3) {
01164                 static Vector temp7(7);
01165                 workV6 = currentStress.t2Vector();
01166     temp7[0] = workV6[0];
01167     temp7[1] = workV6[1];
01168     temp7[2] = workV6[2];
01169     temp7[3] = workV6[3];
01170     temp7[4] = workV6[4];
01171     temp7[5] = workV6[5];
01172     temp7[6] = scale;
01173                 return temp7;
01174         }
01175 
01176   else {
01177     static Vector temp5(5);  
01178                 workV6 = currentStress.t2Vector();
01179     temp5[0] = workV6[0];
01180     temp5[1] = workV6[1];
01181     temp5[2] = workV6[2];
01182     temp5[3] = workV6[3];
01183     temp5[4] = scale;
01184     /*temp5[5] = committedActiveSurf;
01185         temp5[6] = stressRatioPTx[matN];
01186         temp5[7] = currentStress.deviatorRatio(residualPressx[matN]);
01187     temp5[8] = pressureDCommitted;
01188     temp5[9] = cumuDilateStrainOctaCommitted;
01189     temp5[10] = maxCumuDilateStrainOctaCommitted;
01190     temp5[11] = cumuTranslateStrainOctaCommitted;
01191     temp5[12] = onPPZCommitted;
01192     temp5[13] = PPZSizeCommitted;*/
01193     return temp5;
01194   }
01195 }
01196 
01197 const 
01198 Vector & PressureDependMultiYield::getCommittedStrain (void)
01199 {
01200         int ndm = ndmx[matN];
01201 
01202   if (ndm==3)
01203     return currentStrain.t2Vector(1);
01204   else {
01205                 static Vector workV(3);
01206                 workV6 = currentStrain.t2Vector(1);
01207     workV[0] = workV6[0];
01208     workV[1] = workV6[1];
01209     workV[2] = workV6[3];
01210     return workV;
01211   }
01212 }// NOTE: surfaces[0] is not used 
01213 
01214 
01215 void 
01216 PressureDependMultiYield::setUpSurfaces (double * gredu)
01217 { 
01218     double residualPress = residualPressx[matN];
01219     double refPressure = refPressurex[matN];
01220     double pressDependCoeff =pressDependCoeffx[matN];
01221     double refShearModulus = refShearModulusx[matN];
01222     int numOfSurfaces = numOfSurfacesx[matN];
01223     double frictionAngle = frictionAnglex[matN];
01224         double cohesion = cohesionx[matN];
01225     double peakShearStrain = peakShearStrainx[matN];
01226     double phaseTransfAngle = phaseTransfAnglex[matN];
01227         double stressRatioPT = stressRatioPTx[matN];
01228 
01229     double refStrain, peakShear, coneHeight;
01230     double stress1, stress2, strain1, strain2, size, elasto_plast_modul, plast_modul;
01231     double ratio1, ratio2;
01232   
01233         if (gredu==0) {
01234           double sinPhi = sin(frictionAngle * pi/180.);
01235     double Mnys = 6.*sinPhi/(3.-sinPhi);
01236     double sinPhiPT = sin(phaseTransfAngle * pi/180.);
01237     stressRatioPT = 6.*sinPhiPT/(3.-sinPhiPT);
01238                 // tao = cohesion * sqrt(8.0)/3.
01239     residualPress = 2 * cohesion / Mnys;
01240     // a small nonzero residualPress for numerical purpose only
01241     if (residualPress < 0.01) residualPress = 0.01; 
01242     coneHeight = - (refPressure - residualPress);
01243     peakShear = sqrt(2.) * coneHeight * Mnys / 3.; 
01244     refStrain = (peakShearStrain * peakShear) 
01245                                 / (refShearModulus * peakShearStrain - peakShear);
01246 
01247     double stressInc = peakShear / numOfSurfaces;
01248 
01249     for (int ii=1; ii<=numOfSurfaces; ii++){
01250       stress1 = ii * stressInc; 
01251       stress2 = stress1 + stressInc;
01252       ratio1 = 3. * stress1 / sqrt(2.) / coneHeight;
01253       ratio2 = 3. * stress2 / sqrt(2.) / coneHeight;
01254       strain1 = stress1 * refStrain / (refShearModulus * refStrain - stress1);
01255       strain2 = stress2 * refStrain / (refShearModulus * refStrain - stress2);
01256     
01257       if (ratio1 <= stressRatioPT && ratio2 >= stressRatioPT) {
01258         double ratio = (ratio2 - stressRatioPT)/(ratio2 - ratio1);
01259         strainPTOcta = strain2 - ratio * (strain2 - strain1);
01260                         }
01261 
01262       size = ratio1;
01263       elasto_plast_modul = 2.*(stress2 - stress1)/(strain2 - strain1);
01264       if ( (2.*refShearModulus - elasto_plast_modul) <= 0) 
01265         plast_modul = UP_LIMIT;
01266       else 
01267         plast_modul = (2.*refShearModulus * elasto_plast_modul)/
01268                             (2.*refShearModulus - elasto_plast_modul);
01269       if (plast_modul < 0) plast_modul = 0;
01270       if (plast_modul > UP_LIMIT) plast_modul = UP_LIMIT;
01271       if (ii==numOfSurfaces) plast_modul = 0;
01272       workV6.Zero();
01273           //opserr<<size<<endln;
01274       committedSurfaces[ii] = MultiYieldSurface(workV6,size,plast_modul);
01275                 }  // ii  
01276         } 
01277         else {  //user defined surfaces   
01278                 int ii = 2*(numOfSurfaces-1);
01279                 double tmax = refShearModulus*gredu[ii]*gredu[ii+1];
01280                 double Mnys = -(sqrt(3.) * tmax - 2.* cohesion) / refPressure;
01281     residualPress = 2 * cohesion / Mnys;
01282     if (residualPress < 0.01) residualPress = 0.01; 
01283     coneHeight = - (refPressure - residualPress);
01284 
01285     double sinPhi = 3*Mnys /(6+Mnys);
01286                 if (sinPhi<0. || sinPhi>1.) {
01287                         opserr <<"\nNDMaterial " <<this->getTag()<<": Invalid friction angle, please modify ref. pressure or G/Gmax curve."<<endln;
01288      exit(-1);
01289                 } 
01290 
01291                 frictionAngle = asin(sinPhi)*180/pi;
01292                 opserr << "\nNDMaterial " <<this->getTag()<<": Friction angle is "<<frictionAngle<<"\n"<<endln;
01293     if (phaseTransfAngle > frictionAngle) {
01294                         opserr << "\nNDMaterial " <<this->getTag()<<": phase Transformation Angle > friction Angle," 
01295                                    << "will set phase Transformation Angle = friction Angle.\n" <<endln;
01296                         phaseTransfAngle = frictionAngle;
01297                 }
01298                 double sinPhiPT = sin(phaseTransfAngle * pi/180.);
01299     stressRatioPT = 6.*sinPhiPT/(3.-sinPhiPT);
01300 
01301                 for (int i=1; i<numOfSurfaces; i++) {
01302                         int ii = 2*(i-1);
01303                         strain1 = gredu[ii]; 
01304       stress1 = refShearModulus*gredu[ii+1]*strain1; 
01305                         strain2 = gredu[ii+2]; 
01306       stress2 = refShearModulus*gredu[ii+3]*strain2; 
01307 
01308       ratio1 = sqrt(3.) * stress1 / coneHeight;
01309       ratio2 = sqrt(3.) * stress2 / coneHeight;  
01310       if (ratio1 <= stressRatioPT && ratio2 >= stressRatioPT) {
01311         double ratio = (ratio2 - stressRatioPT)/(ratio2 - ratio1);
01312                           // gamma_oct = sqrt(6.0)/3*gamma12
01313         strainPTOcta = sqrt(6.)/3 * (strain2 - ratio * (strain2 - strain1));
01314                         }
01315 
01316       size = ratio1;
01317       elasto_plast_modul = 2.*(stress2 - stress1)/(strain2 - strain1);
01318         
01319                         if ( (2.*refShearModulus - elasto_plast_modul) <= 0) 
01320                                         plast_modul = UP_LIMIT;
01321       else 
01322                                         plast_modul = (2.*refShearModulus * elasto_plast_modul)/
01323                         (2.*refShearModulus - elasto_plast_modul);
01324       if (plast_modul <= 0) {
01325                                 opserr << "\nNDMaterial " <<this->getTag()<<": Surface " << i 
01326                                            << " has plastic modulus < 0.\n Please modify G/Gmax curve.\n"<<endln;
01327        exit(-1);
01328       }
01329       if (plast_modul > UP_LIMIT) plast_modul = UP_LIMIT;
01330 
01331       workV6.Zero();
01332                         //opserr<<size<<" "<<i<<" "<<plast_modul<<" "<<gredu[ii]<<" "<<gredu[ii+1]<<endln;
01333       committedSurfaces[i] = MultiYieldSurface(workV6,size,plast_modul);
01334 
01335                         if (i==(numOfSurfaces-1)) {
01336                                 plast_modul = 0;
01337                                 size = ratio2;
01338                           //opserr<<size<<" "<<i+1<<" "<<plast_modul<<" "<<gredu[ii+2]<<" "<<gredu[ii+3]<<endln;
01339         committedSurfaces[i+1] = MultiYieldSurface(workV6,size,plast_modul);
01340                         }
01341                 }
01342   }  
01343 
01344   residualPressx[matN] = residualPress;
01345   frictionAnglex[matN] = frictionAngle;
01346   cohesionx[matN] = cohesion;
01347   phaseTransfAnglex[matN] = phaseTransfAngle;
01348   stressRatioPTx[matN] = stressRatioPT;
01349 }
01350 
01351 double 
01352 PressureDependMultiYield::yieldFunc(const T2Vector & stress, 
01353                                            const MultiYieldSurface * surfaces, 
01354                                            int surfaceNum)
01355 {
01356   double residualPress = residualPressx[matN];
01357 
01358   double coneHeight = stress.volume() - residualPress;
01359   //workV6 = stress.deviator() - surfaces[surfaceNum].center()*coneHeight;
01360   workV6 = stress.deviator();
01361   workV6.addVector(1.0, surfaces[surfaceNum].center(), -coneHeight);
01362 
01363   double sz = surfaces[surfaceNum].size()*coneHeight;
01364 
01365   return 3./2.*(workV6 && workV6) - sz * sz;
01366 }
01367 
01368 void 
01369 PressureDependMultiYield::deviatorScaling(T2Vector & stress, 
01370                                                const MultiYieldSurface * surfaces, 
01371                                                int surfaceNum)
01372 {
01373   double residualPress = residualPressx[matN];
01374   int numOfSurfaces = numOfSurfacesx[matN];
01375 
01376   double diff = yieldFunc(stress, surfaces, surfaceNum);
01377   double coneHeight = stress.volume() - residualPress;
01378 
01379   if ( surfaceNum < numOfSurfaces && diff < 0. ) {
01380     double sz = -surfaces[surfaceNum].size()*coneHeight;
01381     double deviaSz = sqrt(sz*sz + diff);
01382     static Vector devia(6);
01383     devia = stress.deviator(); 
01384     workV6 = devia;
01385     workV6.addVector(1.0, surfaces[surfaceNum].center(), -coneHeight);
01386     double coeff = (sz-deviaSz) / deviaSz;
01387     if (coeff < 1.e-13) coeff = 1.e-13;
01388     //devia += workV6 * coeff;
01389     devia.addVector(1.0, workV6, coeff);
01390     stress.setData(devia, stress.volume());
01391     deviatorScaling(stress, surfaces, surfaceNum);  // recursive call
01392   }
01393 
01394   if (surfaceNum==numOfSurfaces && fabs(diff) > LOW_LIMIT) {
01395     double sz = -surfaces[surfaceNum].size()*coneHeight;
01396     //workV6 = stress.deviator() * sz/sqrt(diff+sz*sz);
01397     workV6 = stress.deviator();
01398     workV6 *= sz/sqrt(diff+sz*sz);
01399     stress.setData(workV6, stress.volume());
01400   }
01401 }
01402 
01403 void 
01404 PressureDependMultiYield::initSurfaceUpdate(void)
01405 {
01406   double residualPress = residualPressx[matN];
01407   int numOfSurfaces = numOfSurfacesx[matN];
01408 
01409   if (committedActiveSurf == 0) return; 
01410   
01411   double coneHeight = - (currentStress.volume() - residualPress);
01412   static Vector devia(6); 
01413   devia = currentStress.deviator();
01414   double Ms = sqrt(3./2.*(devia && devia));
01415 
01416   if (committedActiveSurf < numOfSurfaces) { // failure surface can't move
01417     //workV6 = devia * (1. - committedSurfaces[committedActiveSurf].size()*coneHeight / Ms);
01418     workV6.addVector(0.0, devia, (1. - committedSurfaces[committedActiveSurf].size()*coneHeight / Ms));
01419          
01420     //workV6 = workV6 / -coneHeight;
01421     workV6 /= -coneHeight;
01422     committedSurfaces[committedActiveSurf].setCenter(workV6);
01423   }
01424 
01425   for (int i=1; i<committedActiveSurf; i++) {
01426     //workV6 = devia * (1. - committedSurfaces[i].size()*coneHeight / Ms);
01427     workV6.addVector(0.0, devia, (1. - committedSurfaces[i].size()*coneHeight / Ms));
01428     //workV6 = workV6 / -coneHeight;
01429     workV6 /= -coneHeight;
01430     committedSurfaces[i].setCenter(workV6); 
01431     theSurfaces[i] = committedSurfaces[i];
01432   }
01433   activeSurfaceNum = committedActiveSurf;
01434 }
01435 
01436 void 
01437 PressureDependMultiYield::initStrainUpdate(void)
01438 {
01439     double residualPress = residualPressx[matN];
01440     double refPressure = refPressurex[matN];
01441     double pressDependCoeff =pressDependCoeffx[matN];
01442     double refShearModulus = refShearModulusx[matN];
01443         double refBulkModulus = refBulkModulusx[matN];
01444     double stressRatioPT = stressRatioPTx[matN];
01445 
01446   // elastic strain state
01447   double stressRatio = currentStress.deviatorRatio(residualPress);
01448   double ratio = (-currentStress.volume()+residualPress)/(-refPressure+residualPress);
01449   ratio = pow(ratio, 1.-pressDependCoeff);
01450   modulusFactor = getModulusFactor(currentStress);
01451   double shearCoeff = 1./(2.*refShearModulus*modulusFactor);
01452   double bulkCoeff = 1./(3.*refBulkModulus*modulusFactor);
01453 
01454   //currentStrain = currentStress.deviator()*shearCoeff
01455   //              + currentStress.volume()*bulkCoeff;
01456 
01457   // modified fmk as discussed with z.yang
01458   workV6.addVector(0.0, currentStress.deviator(), shearCoeff);
01459   currentStrain.setData(workV6, currentStress.volume()*bulkCoeff);
01460         
01461   double octalStrain = currentStrain.octahedralShear(1);
01462   if (octalStrain <= LOW_LIMIT) octalStrain = LOW_LIMIT;
01463 
01464   // plastic strain state (scaled from elastic strain)
01465   double scale, PPZLimit;
01466   if (stressRatio >= stressRatioPT) {  //above PT
01467     onPPZ = 2;
01468     prePPZStrainOcta = strainPTOcta * ratio;
01469     PPZLimit = getPPZLimits(1,currentStress);
01470     scale = sqrt(prePPZStrainOcta+PPZLimit)/octalStrain;
01471   }
01472   else {  // below PT
01473     onPPZ = -1;
01474     prePPZStrainOcta = octalStrain;
01475     if (prePPZStrainOcta > strainPTOcta * ratio) prePPZStrainOcta=strainPTOcta*ratio;
01476     scale = sqrt(prePPZStrainOcta)/octalStrain;
01477   }
01478   //currentStrain.setData(currentStrain.deviator()*scale, currentStrain.volume());
01479   workV6.addVector(0.0, currentStrain.deviator(), scale);
01480   currentStrain.setData(workV6, currentStrain.volume());
01481   PPZPivot = currentStrain;
01482 }
01483 
01484 double 
01485 PressureDependMultiYield::getModulusFactor(T2Vector & stress)
01486 {
01487     double residualPress = residualPressx[matN];
01488     double refPressure = refPressurex[matN];
01489     double pressDependCoeff =pressDependCoeffx[matN];
01490 
01491   double conHeig = stress.volume() - residualPress;
01492   double scale = conHeig / (refPressure-residualPress);  
01493   scale = pow(scale, pressDependCoeff);
01494         
01495   return (1.e-10>scale) ? 1.e-10 : scale; 
01496 }
01497 
01498 void 
01499 PressureDependMultiYield::setTrialStress(T2Vector & stress)
01500 {
01501     double refShearModulus = refShearModulusx[matN];
01502         double refBulkModulus = refBulkModulusx[matN];
01503 
01504   modulusFactor = getModulusFactor(stress);
01505   //workV6 = stress.deviator() 
01506   //                 + subStrainRate.deviator()*2.*refShearModulus*modulusFactor;
01507   workV6 = stress.deviator();
01508   workV6.addVector(1.0, subStrainRate.deviator(), 2*refShearModulus*modulusFactor);
01509   
01510   double B = refBulkModulus*modulusFactor;
01511   
01512   if (Hvx[matN] != 0. && trialStress.volume()<=maxPress && subStrainRate.volume()<0.) {
01513      double tp = fabs(trialStress.volume() - residualPressx[matN]);
01514      B = (B*Hvx[matN]*pow(tp,Pvx[matN]))/(B+Hvx[matN]*pow(tp,Pvx[matN]));
01515   }
01516   
01517   double volume = stress.volume() + subStrainRate.volume()*3.*B;
01518 
01519   if (volume > 0.) volume = 0.;
01520   trialStress.setData(workV6, volume);
01521 }
01522 
01523 int 
01524 PressureDependMultiYield::setSubStrainRate(void)
01525 {
01526     double residualPress = residualPressx[matN];
01527     double refShearModulus = refShearModulusx[matN];
01528         int numOfSurfaces = numOfSurfacesx[matN];
01529 
01530   if (activeSurfaceNum==numOfSurfaces) return 1;
01531   if (strainRate.isZero()) return 0;  
01532 
01533   double elast_plast_modulus;
01534   double conHeig = -(currentStress.volume() - residualPress);
01535   double factor = getModulusFactor(currentStress);
01536   if (activeSurfaceNum==0) 
01537     elast_plast_modulus = 2*refShearModulus*factor;
01538   else {
01539     double plast_modulus = theSurfaces[activeSurfaceNum].modulus()*factor;
01540     elast_plast_modulus = 2*refShearModulus*factor*plast_modulus 
01541       / (2*refShearModulus*factor+plast_modulus);
01542   }
01543   //workV6 = strainRate.deviator()*elast_plast_modulus;
01544   workV6.addVector(0.0, strainRate.deviator(), elast_plast_modulus);
01545   workT2V.setData(workV6,0);
01546 
01547   double singleCross = theSurfaces[numOfSurfaces].size()*conHeig / numOfSurfaces;
01548   double totalCross = 3.*workT2V.octahedralShear() / sqrt(2.);
01549   int numOfSub = totalCross/singleCross + 1;
01550   if (numOfSub > numOfSurfaces) numOfSub = numOfSurfaces;
01551         
01552   int numOfSub1 = strainRate.octahedralShear(1) / 1.0e-4;
01553   int numOfSub2 = strainRate.volume() / 1.e-5;
01554   if (numOfSub1 > numOfSub) numOfSub = numOfSub1;
01555   if (numOfSub2 > numOfSub) numOfSub = numOfSub2;
01556 
01557   workV6.addVector(0.0, strainRate.t2Vector(), 1.0/numOfSub);
01558 
01559   subStrainRate.setData(workV6);
01560 
01561   return numOfSub;
01562 }
01563 
01564 void
01565 PressureDependMultiYield::getContactStress(T2Vector &contactStress)
01566 {
01567     double residualPress = residualPressx[matN];
01568 
01569   double conHeig = trialStress.volume() - residualPress;
01570   static Vector center(6);
01571   center = theSurfaces[activeSurfaceNum].center(); 
01572   //workV6 = trialStress.deviator() - center*conHeig;
01573   workV6 = trialStress.deviator();
01574   workV6.addVector(1.0, center, -conHeig);
01575   double Ms = sqrt(3./2.*(workV6 && workV6));
01576   //workV6 = workV6 * theSurfaces[activeSurfaceNum].size()*(-conHeig) / Ms + center*conHeig;
01577   workV6.addVector(theSurfaces[activeSurfaceNum].size()*(-conHeig) / Ms, center, conHeig);
01578   //return T2Vector(workV6,trialStress.volume()); 
01579   contactStress.setData(workV6,trialStress.volume()); 
01580 }
01581 
01582 int 
01583 PressureDependMultiYield::isLoadReversal(const T2Vector & stress)
01584 {
01585   if(activeSurfaceNum == 0) return 0;
01586 
01587   getSurfaceNormal(stress, workT2V);
01588 
01589   //if (((trialStress.t2Vector() - currentStress.t2Vector()) 
01590   //    && workT2V.t2Vector()) < 0) return 1;
01591   workV6 = trialStress.t2Vector();
01592   workV6 -= currentStress.t2Vector();
01593  
01594   if ((workV6 && workT2V.t2Vector()) < 0) return 1; 
01595 
01596   return 0;   
01597 }
01598 
01599 void
01600 PressureDependMultiYield::getSurfaceNormal(const T2Vector & stress, T2Vector &normal)
01601 {
01602     double residualPress = residualPressx[matN];
01603 
01604   double conHeig = stress.volume() - residualPress;
01605   workV6 = stress.deviator();
01606   static Vector center(6);
01607   center = theSurfaces[activeSurfaceNum].center(); 
01608   double sz = theSurfaces[activeSurfaceNum].size();
01609   double volume = conHeig*((center && center) - 2./3.*sz*sz) - (workV6 && center);
01610   //workT2V.setData((workV6-center*conHeig)*3., volume);
01611   workV6.addVector(1.0, center, -conHeig);
01612   workV6 *= 3.0;
01613   workT2V.setData(workV6, volume);
01614   
01615   normal.setData(workT2V.unitT2Vector());
01616 }
01617 
01618 double 
01619 PressureDependMultiYield::getPlasticPotential(const T2Vector & contactStress,
01620                                               const T2Vector & surfaceNormal)
01621 {
01622     double residualPress = residualPressx[matN];
01623     double stressRatioPT = stressRatioPTx[matN];
01624     int numOfSurfaces = numOfSurfacesx[matN];
01625         double contractParam1 = contractParam1x[matN];
01626     double dilateParam1 = dilateParam1x[matN];
01627     double dilateParam2 = dilateParam2x[matN];
01628 
01629   double plasticPotential, contractRule, unloadRule, dilateRule, shearLoading, temp;
01630 
01631   double contactRatio = contactStress.deviatorRatio(residualPress);
01632   temp = contactRatio/stressRatioPT;
01633   double factorPT = (temp*temp - 1)/(temp*temp + 1)/3.;
01634   double volume = contactStress.volume();
01635   contractRule = factorPT*contractParam1;
01636   if (contractRule > 0.) contractRule = -contractRule;
01637         if (contractRule<-5.0e4) contractRule = -5.0e4;
01638   temp = currentStress.volume() - pressureD;
01639   if (temp >= 0.) unloadRule = 0.;
01640   else {
01641     double conHeiD = pressureD-residualPress;
01642     double temp1 = sqrt(3./2.)*currentStress.deviatorLength() 
01643       + stressRatioPT*conHeiD;
01644     temp = temp1/(-temp);
01645     if (temp < theSurfaces[numOfSurfaces].size()) 
01646       temp = theSurfaces[numOfSurfaces].size();
01647     temp1 = (reversalStress.volume()-residualPress)/conHeiD;
01648     unloadRule = -sqrt(3./2.)*surfaceNormal.deviatorLength()*temp1/temp;
01649   }
01650           
01651   double currentRatio = currentStress.deviatorRatio(residualPress);
01652   double trialRatio = trialStress.deviatorRatio(residualPress);
01653   shearLoading = currentStress.deviator() && trialStress.deviator();
01654 
01655   if (factorPT < 0.) {  //below PT
01656     if (pressureD == 0.) 
01657       plasticPotential = contractRule;
01658     else if (trialStress.volume() >= pressureD) {
01659       pressureD = 0.;
01660       plasticPotential = contractRule;
01661     }
01662     else if (trialRatio > currentRatio && shearLoading >= 0.)
01663       plasticPotential = contractRule;
01664     else  plasticPotential = unloadRule;
01665   }
01666   
01667   else {  //above PT
01668     if (trialRatio > currentRatio && shearLoading >= 0.) {  //dilation
01669       if (pressureD == 0.) pressureD = currentStress.volume();
01670       reversalStress = currentStress;
01671       updatePPZ(contactStress);  
01672       if (onPPZ==-1 || onPPZ==1) return LOCK_VALUE; 
01673       if (isCriticalState(contactStress))  
01674               dilateRule = 0;
01675       else
01676         dilateRule = factorPT*dilateParam1*exp(dilateParam2*cumuDilateStrainOcta);
01677                         if (dilateRule>5.0e4) dilateRule = 5.0e4;
01678       return dilateRule;
01679     }
01680     else if (pressureD == 0.) plasticPotential = contractRule;
01681     else if (trialStress.volume() >= pressureD) {
01682       pressureD = 0.;
01683       plasticPotential = contractRule;
01684     }
01685     else plasticPotential = unloadRule;
01686   }
01687 
01688   if (onPPZ > 0) onPPZ = 0;
01689   if (onPPZ != -1) PPZTranslation(contactStress);  
01690   if (isCriticalState(contactStress)) return 0.;
01691   return plasticPotential;
01692 }
01693 
01694 int 
01695 PressureDependMultiYield::isCriticalState(const T2Vector & stress)
01696 {
01697         double einit = einitx[matN];
01698         double volLimit1 = volLimit1x[matN];
01699         double volLimit2 = volLimit2x[matN];
01700         double volLimit3 = volLimit3x[matN];
01701 
01702   double vol = trialStrain.volume()*3.0;
01703         double etria = einit + vol + vol*einit;
01704         vol = currentStrain.volume()*3.0;
01705         double ecurr = einit + vol + vol*einit;
01706  
01707         double ecr1, ecr2;
01708         if (volLimit3 != 0.) {
01709                 ecr1 = volLimit1 - volLimit2*pow(abs(-stress.volume()/pAtm), volLimit3);
01710           ecr2 = volLimit1 - volLimit2*pow(abs(-currentStress.volume()/pAtm), volLimit3);
01711         } else {
01712                 ecr1 = volLimit1 - volLimit2*log(abs(-stress.volume()/pAtm));
01713           ecr2 = volLimit1 - volLimit2*log(abs(-currentStress.volume()/pAtm));
01714   }
01715 
01716         if (ecurr < ecr2 && etria < ecr1) return 0;
01717         if (ecurr > ecr2 && etria > ecr1) return 0;     
01718         return 1;
01719 }
01720 
01721 void 
01722 PressureDependMultiYield::updatePPZ(const T2Vector & contactStress)
01723 {
01724   double liquefyParam1 = liquefyParam1x[matN];
01725   double residualPress = residualPressx[matN];
01726   double refPressure = refPressurex[matN];
01727   double pressDependCoeff =pressDependCoeffx[matN];
01728 
01729   // PPZ inactive if liquefyParam1==0.
01730   if (liquefyParam1==0.) {
01731     if (onPPZ==2) {
01732                   workT2V.setData(trialStrain.t2Vector() - PPZPivot.t2Vector());
01733       cumuDilateStrainOcta = workT2V.octahedralShear(1);
01734     }
01735     else if (onPPZ != 2) {
01736       onPPZ = 2;
01737       PPZPivot = trialStrain;
01738       cumuDilateStrainOcta = 0.;
01739     }
01740     return;
01741   }
01742 
01743   // dilation: calc. cumulated dilative strain
01744   if (onPPZ==2) {
01745     PPZPivot = trialStrain;
01746     workV6 = PPZPivot.t2Vector();
01747     workV6 -= PPZCenter.t2Vector();
01748     workT2V.setData(workV6);
01749     //cumuDilateStrainOcta = workT2V.octahedralShear(1) - PPZSize;
01750                 cumuDilateStrainOcta += subStrainRate.octahedralShear(1);
01751     if (cumuDilateStrainOcta > maxCumuDilateStrainOcta) 
01752       maxCumuDilateStrainOcta = cumuDilateStrainOcta;
01753     return;
01754   }
01755 
01756   // calc. PPZ size.
01757   double PPZLimit = getPPZLimits(1,contactStress);
01758         double TransLimit = getPPZLimits(2,contactStress);
01759         //if (PPZLimit==0.) return;
01760 
01761   if (onPPZ==-1 || onPPZ==0) {
01762     workV6 = trialStrain.t2Vector();
01763     workV6 -= PPZPivot.t2Vector();
01764     workT2V.setData(workV6);
01765                 double temp = workT2V.octahedralShear(1);
01766                 if (temp > cumuDilateStrainOcta) {
01767       double volume = -contactStress.volume();
01768       oppoPrePPZStrainOcta = prePPZStrainOcta;
01769       double ratio = (volume+residualPress)/(-refPressure+residualPress);
01770       ratio = pow(ratio, 1.-pressDependCoeff);
01771       prePPZStrainOcta = ratio * strainPTOcta;
01772       if (oppoPrePPZStrainOcta == 0.) oppoPrePPZStrainOcta = prePPZStrainOcta;
01773                 }
01774   }
01775   if (onPPZ > -1) PPZSize = PPZLimit 
01776     + (prePPZStrainOcta+oppoPrePPZStrainOcta+TransLimit+maxCumuDilateStrainOcta)/2.;
01777         else PPZSize = PPZLimit 
01778     + (prePPZStrainOcta+oppoPrePPZStrainOcta+maxCumuDilateStrainOcta)/2.;
01779 
01780   // calc. new PPZ center.
01781   if (onPPZ==0 || onPPZ==1) { 
01782     //workT2V.setData(PPZPivot.t2Vector() - PPZCenter.t2Vector());
01783     workV6 = PPZPivot.t2Vector();
01784     workV6 -= PPZCenter.t2Vector();
01785     workT2V.setData(workV6);
01786     
01787     double coeff = (PPZSize-cumuTranslateStrainOcta)/workT2V.octahedralShear(1);
01788     //PPZCenter.setData(PPZPivot.t2Vector() - workT2V.t2Vector()*coeff);
01789     workV6 = PPZPivot.t2Vector();
01790     workV6.addVector(1.0, workT2V.t2Vector(), -coeff);
01791     PPZCenter.setData(workV6);
01792   }
01793 
01794   //workT2V.setData(trialStrain.t2Vector() - PPZCenter.t2Vector());
01795   workV6 = trialStrain.t2Vector();
01796   workV6 -= PPZCenter.t2Vector();
01797   workT2V.setData(workV6);
01798   double temp = subStrainRate.t2Vector() && workV6;
01799 
01800   if (workT2V.octahedralShear(1) > PPZSize && temp > 0. || PPZLimit==0.) {  //outside PPZ
01801     workV6 = trialStrain.t2Vector();
01802     workV6 -= PPZPivot.t2Vector();
01803     workT2V.setData(workV6);
01804                 double temp1 = workT2V.octahedralShear(1);
01805                 if (temp1 > cumuDilateStrainOcta) {
01806       cumuDilateStrainOcta = 0.;
01807       if (PPZLimit == 0.) maxCumuDilateStrainOcta = 0.;
01808                 }
01809     onPPZ = 2;
01810     PPZPivot = trialStrain;
01811     cumuTranslateStrainOcta = 0.;
01812   }
01813   else {  //inside PPZ
01814     if (onPPZ == 0 || onPPZ == 1) PPZTranslation(contactStress);
01815     if (onPPZ == -1 || onPPZ == 0) lockStress = contactStress;
01816     if (onPPZ == 0) onPPZ = 1;
01817   }
01818 }
01819 
01820 void 
01821 PressureDependMultiYield::PPZTranslation(const T2Vector & contactStress)
01822 {
01823         double liquefyParam1 = liquefyParam1x[matN];
01824 
01825   if (liquefyParam1==0.) return;
01826 
01827   double PPZLimit = getPPZLimits(1,contactStress);
01828         if (PPZLimit==0.) return;
01829 
01830   double PPZTransLimit = getPPZLimits(2,contactStress);
01831 
01832   //workT2V.setData(trialStrain.deviator()-PPZPivot.deviator());
01833   workV6 = trialStrain.deviator();
01834   workV6 -= PPZPivot.deviator();
01835   workT2V.setData(workV6);
01836   
01837   double temp = workT2V.octahedralShear(1);
01838   if (cumuTranslateStrainOcta < temp) cumuTranslateStrainOcta = temp;
01839         if (maxCumuDilateStrainOcta == 0.) temp = PPZTransLimit;
01840         else temp = PPZTransLimit*cumuDilateStrainOcta/maxCumuDilateStrainOcta;
01841   if (cumuTranslateStrainOcta > temp) cumuTranslateStrainOcta = temp;
01842 }
01843 
01844 double 
01845 PressureDependMultiYield::getPPZLimits(int which, const T2Vector & contactStress)
01846 {
01847         double liquefyParam1 = liquefyParam1x[matN];
01848         double liquefyParam2 = liquefyParam2x[matN];
01849         double liquefyParam4 = liquefyParam4x[matN];
01850 
01851   double PPZLimit, temp;
01852   double volume = -contactStress.volume();
01853 
01854   if (volume >= liquefyParam1) PPZLimit = 0.;
01855   else {
01856     temp = volume*pi/liquefyParam1/2.;
01857                 // liquefyParam3 = 3.0 default
01858     PPZLimit = liquefyParam2 * pow(cos(temp), 3.);
01859   }
01860   
01861   if (which==1) 
01862     return PPZLimit;
01863   else if (which==2) 
01864     return liquefyParam4 * PPZLimit;
01865   else {
01866     opserr << "FATAL:PressureDependMultiYield::getPPZLimits: unknown argument value" << endln;
01867    exit(-1);
01868     return 0.0;
01869   }
01870 }
01871 
01872 double 
01873 PressureDependMultiYield::getLoadingFunc(const T2Vector & contactStress, 
01874                                          const T2Vector & surfaceNormal,
01875                                          double plasticPotential,
01876                                          int crossedSurface)
01877 {
01878     int numOfSurfaces = numOfSurfacesx[matN];
01879     double refShearModulus = refShearModulusx[matN];
01880         double refBulkModulus = refBulkModulusx[matN];
01881 
01882   double loadingFunc, limit;
01883   double modul = theSurfaces[activeSurfaceNum].modulus();
01884   double temp1 = 2. * refShearModulus * modulusFactor 
01885     * (surfaceNormal.deviator() && surfaceNormal.deviator()) ;
01886   double temp2 = 9. * refBulkModulus * modulusFactor 
01887     * surfaceNormal.volume() * plasticPotential ;
01888         
01889   //for the first crossing 
01890   double temp = temp1 + temp2 + modul * modulusFactor;
01891   if (activeSurfaceNum == numOfSurfaces) 
01892     limit = theSurfaces[activeSurfaceNum-1].modulus() * modulusFactor /2.;
01893   else limit = modul * modulusFactor /2.;
01894   if (temp < limit) {
01895     plasticPotential = (temp2 + limit - temp)/(9. * refBulkModulus * modulusFactor 
01896                                                * surfaceNormal.volume());
01897     temp = limit;
01898   }
01899   //loadingFunc = (surfaceNormal.t2Vector() 
01900   //                 && (trialStress.deviator()-contactStress.deviator()))/temp;
01901   workV6 = trialStress.deviator();
01902   workV6 -= contactStress.deviator();  
01903   loadingFunc = (surfaceNormal.t2Vector() && workV6) /temp;
01904   
01905   if (loadingFunc < 0.) loadingFunc = 0;
01906 
01907   //for more than one crossing 
01908   if(crossedSurface) {
01909     temp = (theSurfaces[activeSurfaceNum-1].modulus() - modul)
01910       /theSurfaces[activeSurfaceNum-1].modulus();
01911     loadingFunc *= temp;
01912   }
01913 
01914   return loadingFunc;
01915 }
01916 
01917 int 
01918 PressureDependMultiYield::stressCorrection(int crossedSurface)
01919 {
01920     double refShearModulus = refShearModulusx[matN];
01921         double refBulkModulus = refBulkModulusx[matN];
01922 
01923   static T2Vector contactStress;
01924   getContactStress(contactStress);
01925   static T2Vector surfNormal;
01926   getSurfaceNormal(contactStress, surfNormal);
01927   double plasticPotential = getPlasticPotential(contactStress,surfNormal);
01928   if (plasticPotential==LOCK_VALUE && (onPPZ == -1 || onPPZ == 1)) {
01929     trialStress = lockStress;
01930     return 1;
01931   }
01932 
01933         double tVolume = trialStress.volume();
01934   double loadingFunc = getLoadingFunc(contactStress, surfNormal, 
01935                                       plasticPotential, crossedSurface);
01936   double volume = tVolume - plasticPotential*3*refBulkModulus*modulusFactor*loadingFunc;
01937 
01938   //workV6 = trialStress.deviator() 
01939   //             - surfNormal.deviator()*2*refShearModulus*modulusFactor*loadingFunc;
01940   workV6 = trialStress.deviator();
01941         
01942   if (volume > 0. && volume != tVolume) {
01943                 double coeff = tVolume / (tVolume - volume);
01944     coeff *= -2*refShearModulus*modulusFactor*loadingFunc;
01945     workV6.addVector(1.0, surfNormal.deviator(), coeff);
01946                 volume = 0.;
01947   } else if (volume > 0.) {
01948                 volume = 0.;
01949         } else {
01950                 double coeff = -2*refShearModulus*modulusFactor*loadingFunc;
01951     workV6.addVector(1.0, surfNormal.deviator(), coeff);
01952         }
01953         /*
01954           if (volume>0.)volume = 0.;
01955                 double coeff = -2*refShearModulus*modulusFactor*loadingFunc;
01956     workV6.addVector(1.0, surfNormal.deviator(), coeff);
01957   */
01958         trialStress.setData(workV6, volume);
01959   deviatorScaling(trialStress, theSurfaces, activeSurfaceNum);
01960 
01961   if (isCrossingNextSurface()) {
01962     activeSurfaceNum++;
01963     return stressCorrection(1);  //recursive call
01964   }
01965   return 0;
01966 }
01967 
01968 void 
01969 PressureDependMultiYield::updateActiveSurface(void)
01970 {
01971     double residualPress = residualPressx[matN];
01972     int numOfSurfaces = numOfSurfacesx[matN];
01973 
01974   if (activeSurfaceNum == numOfSurfaces) return;
01975 
01976   double A, B, C, X;
01977   static Vector t1(6);
01978   static Vector t2(6);
01979   static Vector center(6);
01980   static Vector outcenter(6);
01981   double conHeig = trialStress.volume() - residualPress;
01982   center = theSurfaces[activeSurfaceNum].center();
01983   double size = theSurfaces[activeSurfaceNum].size();
01984   outcenter = theSurfaces[activeSurfaceNum+1].center();
01985   double outsize = theSurfaces[activeSurfaceNum+1].size();
01986 
01987   //t1 = trialStress.deviator() - center*conHeig;
01988   //t2 = (center - outcenter)*conHeig;
01989   t1 = trialStress.deviator();
01990   t1.addVector(1.0, center, -conHeig);
01991   t2 = center;
01992   t2 -= outcenter;
01993   t2 *=conHeig;
01994  
01995   A = t1 && t1;
01996   B = 2. * (t1 && t2);
01997   C = (t2 && t2) - 2./3.* outsize * outsize * conHeig * conHeig;
01998   X = secondOrderEqn(A,B,C,0);
01999   if ( fabs(X-1.) < LOW_LIMIT ) X = 1.;
02000   if (X < 1.) return;
02001 
02002   if (X < 1.){
02003     //t2 = trialStress.deviator() - outcenter*conHeig;
02004     t2 = trialStress.deviator();
02005     t2.addVector(1.0, outcenter, -conHeig);
02006     
02007     double xx1 = (t2 && t2) - 2./3.* outsize * outsize * conHeig * conHeig;
02008     double xx2 = (t1 && t1) - 2./3.* size * size * conHeig * conHeig;
02009     opserr << "FATAL:PressureDependMultiYield::updateActiveSurface(): error in Direction of surface motion." << endln; 
02010     opserr << "X-1= " << X-1 <<" A= "<<A<<" B= "<<B<<" C= "<<C <<" M= "<<activeSurfaceNum<<" low_limit="<<LOW_LIMIT<< endln;
02011     opserr << "diff1= "<<xx1 <<" diff2= "<<xx2 <<" p= "<<conHeig<<" size= "<<size<<" outs= "<<outsize<<endln; 
02012    exit(-1);
02013   }
02014 
02015   //workV6 = (t1 * X + center*conHeig) * (1. - size / outsize) 
02016   //         - (center - outcenter * size / outsize) * conHeig;
02017   
02018   workV6.addVector(0.0, t1, X);
02019   workV6.addVector(1.0, center, conHeig);
02020   workV6 *= (1.0 - size / outsize);
02021   t2 = center;
02022   t2.addVector(1.0, outcenter, -size/outsize);
02023   t2 *= conHeig;
02024   workV6 -= t2;
02025 
02026   workT2V.setData(workV6);
02027   if (workT2V.deviatorLength() < LOW_LIMIT) return;
02028 
02029   workV6 = workT2V.deviator();  
02030   A = conHeig * conHeig * (workV6 && workV6);
02031   B = 2 * conHeig * (t1 && workV6);
02032   if (fabs(B) < LOW_LIMIT) B = 0.; 
02033   C = (t1 && t1) - 2./3.* size * size * conHeig * conHeig;
02034   if ( fabs(C) < LOW_LIMIT || fabs(C)/(t1 && t1) < LOW_LIMIT ) return;
02035   if (B > 0. || C < 0.) {
02036     opserr << "FATAL:PressureDependMultiYield::updateActiveSurface(): error in surface motion.\n" 
02037          << "A= " <<A <<" B= " <<B <<" C= "<<C <<" (t1&&t1)= "<<(t1&&t1) <<endln; 
02038    exit(-1);
02039   }
02040   X = secondOrderEqn(A,B,C,1);  
02041 
02042   center.addVector(1.0, workV6, -X);
02043   theSurfaces[activeSurfaceNum].setCenter(center);
02044 }      
02045 
02046 void 
02047 PressureDependMultiYield::updateInnerSurface(void)
02048 {
02049     double residualPress = residualPressx[matN];
02050 
02051         if (activeSurfaceNum <= 1) return;
02052         static Vector devia(6);
02053         static Vector center(6);
02054 
02055         double conHeig = currentStress.volume() - residualPress;
02056         devia = currentStress.deviator();
02057         center = theSurfaces[activeSurfaceNum].center();
02058         double size = theSurfaces[activeSurfaceNum].size();
02059 
02060         for (int i=1; i<activeSurfaceNum; i++) {
02061                 workV6.addVector(0.0, center, conHeig);
02062                 workV6 -= devia;
02063                 workV6 *= theSurfaces[i].size()/size;
02064                 workV6 += devia;
02065 
02066                 //workV6 = devia - (devia - center*conHeig) * theSurfaces[i].size() / size;
02067                 
02068                 workV6 /= conHeig;
02069                 theSurfaces[i].setCenter(workV6);
02070         }
02071 }
02072 
02073 int 
02074 PressureDependMultiYield:: isCrossingNextSurface(void)
02075 {
02076     int numOfSurfaces = numOfSurfacesx[matN];
02077 
02078   if (activeSurfaceNum == numOfSurfaces) return 0;  
02079 
02080   if(yieldFunc(trialStress, theSurfaces, activeSurfaceNum+1) > 0) return 1;
02081   
02082   return 0;
02083 }
Generated on Mon Oct 23 15:05:15 2006 for OpenSees by 1.5.0
opensees-support @ berkeley.edu
Supported by the National Science Foundation