PressureDependMultiYield02.cpp

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