NeoPZ
TPZYCSandlerDimaggio.h
Go to the documentation of this file.
1 // $Id: TPZYCSandlerDimaggio.h,v 1.11 2009-06-29 22:54:01 erick Exp $
2 
3 #ifndef TPZYCSANDLERDIMAGGIO_H
4 #define TPZYCSANDLERDIMAGGIO_H
5 
6 #include "TPZTensor.h"
7 #include "pzfmatrix.h"
8 #include "TPZElasticResponse.h"
9 #include "pzlog.h"
10 
11 #ifndef CHECKCONV
12 #define CHECKCONV
13 #include "checkconv.h"
14 #endif
15 
16 #include "fadType.h"
17 #include "TPZPlasticCriterion.h"
18 
19 #ifdef LOG4CXX
20 #include "pzlog.h"
21 
22 static LoggerPtr loggerSM(Logger::getLogger("material.plasticity.SM"));
23 
24 #endif
25 
31 class TPZYCSandlerDimaggio : public TPZPlasticCriterion {
32 public:
33 
34  enum {
35  NYield = 2
36  };
37 
38  virtual int ClassId() const override;
39 
40  TPZYCSandlerDimaggio() : fA(0.), fB(0.), fC(0.), fD(0.), fW(0.), fR(0.), fIsonCap(false) {
41  }
42 
43  TPZYCSandlerDimaggio(const TPZYCSandlerDimaggio & source) {
44  fA = source.fA;
45  fB = source.fB;
46  fC = source.fC;
47  fD = source.fD;
48  fW = source.fW;
49  fR = source.fR;
50  fIsonCap = source.fIsonCap;
51  }
52 
53  TPZYCSandlerDimaggio & operator=(const TPZYCSandlerDimaggio & source) {
54  fA = source.fA;
55  fB = source.fB;
56  fC = source.fC;
57  fD = source.fD;
58  fW = source.fW;
59  fR = source.fR;
60  fIsonCap = source.fIsonCap;
61  return *this;
62  }
63 
64  void Write(TPZStream& buf, int withclassid) const override {
65  buf.Write(&fA);
66  buf.Write(&fB);
67  buf.Write(&fC);
68  buf.Write(&fD);
69  buf.Write(&fW);
70  buf.Write(&fR);
71  }
72 
73  void Read(TPZStream& buf, void* context) override {
74  buf.Read(&fA);
75  buf.Read(&fB);
76  buf.Read(&fC);
77  buf.Read(&fD);
78  buf.Read(&fW);
79  buf.Read(&fR);
80  }
81 
82  virtual ~TPZYCSandlerDimaggio() {
83 
84  }
85 
86  const char * Name() const {
87  return "TPZYCSandlerDimaggio";
88  }
89 
90  void Print(std::ostream & out) const override {
91  out << "\n" << this->Name();
92  out << "\n fA = " << fA;
93  out << "\n fB = " << fB;
94  out << "\n fC = " << fC;
95  out << "\n fD = " << fD;
96  out << "\n fR = " << fR;
97  out << "\n fW = " << fW;
98  out << "\n IsonCap " << fIsonCap;
99  }
100 
101  int GetForceYield() {
102  return 0; // nothing to be done in this yield criterium
103  }
104 
105  void SetForceYield(const int forceYield) {
106  // nothing to be done in this yield criterium
107  }
108 
116  void SetYieldStatusMode(const TPZTensor<REAL> & sigma, const REAL & A) {
117  // nothing to be done in this yield criterium
118  }
119 
127  template < class T>
128  void Compute(const TPZTensor<T> & sigma, const T & A, TPZVec<T> &res, int checkForcedYield) const;
129 
137  template <class T>
138  void N(const TPZTensor<T> & sigma, const T & A, TPZVec<TPZTensor<T> > & Ndir, int checkForcedYield) const;
139 
147  template <class T>
148  void H(const TPZTensor<T> & sigma, const T & A, TPZVec<T> & h, int checkForcedYield) const;
149 
150  inline void SetUp(REAL A, REAL B, REAL C, REAL D, REAL R, REAL W) {
151  fA = A;
152  fB = B;
153  fC = C;
154  fD = D;
155  fR = R;
156  fW = W;
157  fIsonCap = false;
158  }
159 
163  template <class T>
164  void AlphaMultiplier(const T &A, T &multiplier) const {
165  multiplier = T(1.);
166  }
167 
168  inline REAL InitialDamage() {
169  return 0.;
170  }
171 
182  void InitialGuess(const TPZElasticResponse &ER, REAL epsp, TPZTensor<REAL> &sigtrial, REAL &epspproj,
183  TPZVec<REAL> &delgamma, TPZTensor<REAL> &sigproj);
184 
188  REAL FZero() const {
189  return log(0.9 * fA / fC) / fB;
190  }
191 
195  REAL LMax() const {
196  return FZero();
197  }
198 
202  template<class T>
203  void EpspFromL(const T &L, T &epsp) const;
204 
208  template<class T>
209  void DEpspDL(const T &L, T& depspdL) const;
210 
214  template<class T>
215  void D2EpspDL2(const T &L, T& d2epspdL2) const;
216 
217  void YieldFunction(const TPZVec<STATE>& sigma, STATE kprev, TPZVec<STATE>& yield) const override{
218  TPZTensor<STATE> sigmaTensor;
219  sigmaTensor.XX() = sigma[0];
220  sigmaTensor.YY() = sigma[1];
221  sigmaTensor.ZZ() = sigma[2];
222  Compute(sigmaTensor, kprev, yield, 0);
223  }
224 
225  virtual int GetNYield() const override {
226  return as_integer(NYield);
227  }
228 
236  template <class T>
237  void SolveL(const T & X, T & L, REAL relTol = 1.e-10) const;
238 
239 protected:
245  template <class T>
246  void ComputeDL(const T &L, const T &A, T &DL) const;
251  template <class T>
252  void LInitialGuess(const T & X, T & L) const;
253 
254 
255 public:
259  template <class T>
260  void ComputeF(const T & L, T & F) const;
261 
265  template <class T>
266  void ComputedF(const T & L, T & dF) const;
267 
268 protected:
272  void ComputeD2F(REAL L, REAL &D2F) const;
273 
277  REAL Distance(const TPZElasticResponse &ER, REAL L, TPZVec<REAL> &sigtrialIJ) const;
278 
282  REAL DDistance(const TPZElasticResponse &ER, REAL L, TPZVec<REAL> &sigtrialIJ) const;
283 
287  REAL D2Distance(const TPZElasticResponse &ER, REAL L, TPZVec<REAL> &sigtrialIJ) const;
288 
289 public:
295  template <class T>
296  void ComputeX(const T & A, T & X) const;
297 
298 protected:
299 
303  REAL ComputeEpsp(const REAL L) const {
304  REAL FL;
305  ComputeF(L, FL);
306  REAL X = L - fR*FL;
307  REAL eps = fW * (exp(fD * X) - 1);
308  return eps;
309  }
310 
314  REAL FuncEpsp(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &sigtrialIJ) {
315  REAL X, L;
316  ComputeX(epsp + delepsp, X);
317  LInitialGuess(X, L);
318  SolveL(X, L);
319  REAL F;
320  ComputeF(L, F);
321  REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
322  REAL fac1 = 3. * K*delepsp;
323  REAL fac2 = (sigtrialIJ[0]-(L + F * fR * cos(theta)));
324 
325  REAL funcepsp = fac1 - fac2;
326  return funcepsp;
327  }
328 
332  REAL FuncEpspUsingL(const TPZElasticResponse &ER, REAL theta, REAL epspini, REAL L, TPZVec<REAL> &sigtrialIJ) {
333  REAL epsp = ComputeEpsp(L);
334  REAL delepsp = epsp - epspini;
335  REAL F;
336  ComputeF(L, F);
337  REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
338  REAL fac1 = 3. * K*delepsp;
339  REAL fac2 = (sigtrialIJ[0]-(L + F * fR * cos(theta)));
340 
341  REAL funcepsp = fac1 - fac2;
342  return funcepsp;
343  }
344 
348  void DFuncEpspUsingL(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec<REAL> &sigtrialIJ, TPZVec<REAL> &result) const {
349  REAL DepspDL, Dtheta, DerivL;
350  DEpspDL(L, DepspDL);
351  REAL F, DF;
352  ComputeF(L, F);
353  ComputedF(L, DF);
354  Dtheta = -F * fR * sin(theta);
355  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
356  DerivL = 3. * K * DepspDL + (1. + DF * fR * cos(theta));
357  result[0] = Dtheta;
358  result[1] = DerivL;
359 
360  }
361 
365  void DFuncEpsp(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &result) const {
366  REAL X, L, DL, F, DF, Dtheta, Depsp;
367  ComputeX(epsp + delepsp, X);
368  LInitialGuess(X, L);
369  SolveL(X, L);
370  ComputeDL(L, epsp + delepsp, DL);
371  ComputeF(L, F);
372  ComputedF(L, DF);
373  DF *= DL;
374  Dtheta = -F * fR * sin(theta);
375  REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
376  Depsp = 3. * K + (DL + DF * fR * cos(theta));
377  result[0] = Dtheta;
378  result[1] = Depsp;
379  }
380 
384  REAL FuncEpspL(const TPZElasticResponse &ER, REAL theta, REAL L, REAL deltaL, TPZVec<REAL> &sigtrialIJ) {
385  REAL DepspDL;
386  DEpspDL(L, DepspDL);
387  REAL F;
388  ComputeF(L, F);
389  REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
390  REAL fac1 = 3. * K * deltaL*DepspDL;
391  REAL fac2 = (sigtrialIJ[0]-(L + F * fR * cos(theta)));
392  REAL funcepsp = fac1 - fac2;
393 #ifdef LOG4CXX
394  {
395  std::stringstream sout;
396  sout << "funcepsp " << funcepsp << " theta " << theta << " L " << L << " deltaL " << deltaL;
397  LOGPZ_DEBUG(loggerSM, sout.str())
398  }
399 #endif
400 #ifdef PZDEBUG2
401  REAL epsp, epspini;
402  EpspFromL(L, epsp);
403  EpspFromL(L - deltaL, epspini);
404  // REAL funcepsp2 = FuncEpsp(ER, theta, epspini, epsp-epspini, sigtrialIJ);
405  // REAL diff = 3.*K*((deltaL*DepspDL)-(epsp-epspini));
406  // REAL diff2 = funcepsp-funcepsp2;
407 #endif
408  return funcepsp;
409  }
410 
414  void DFuncEpspL(const TPZElasticResponse &ER, REAL theta, REAL L, REAL deltaL, TPZVec<REAL> &sigtrialIJ, TPZVec<REAL> &result) const {
415  REAL DepspDL, D2epspDL2, Dtheta, DerivL;
416  DEpspDL(L, DepspDL);
417  D2EpspDL2(L, D2epspDL2);
418  REAL F, DF;
419  ComputeF(L, F);
420  ComputedF(L, DF);
421  Dtheta = -F * fR * sin(theta);
422  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
423  DerivL = 3. * K * DepspDL + 3. * K * deltaL * D2epspDL2 + (1. + DF * fR * cos(theta));
424  result[0] = Dtheta;
425  result[1] = DerivL;
426 
427  }
428 
432  REAL FuncThetaL(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec<REAL> &sigtrialIJ) const {
433  REAL I1 = sigtrialIJ[0];
434  REAL sqJ2 = sigtrialIJ[1];
435  REAL F;
436  ComputeF(L, F);
437  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
438  const REAL G = ER.Mu();
439  REAL y = 9. * K * (sqJ2 - F * sin(theta));
440  REAL x = G * fR * (I1 - (L + F * fR * cos(theta)));
441  REAL res = theta - atan2(y, x);
442  return res;
443  }
444 
448  REAL FuncTheta(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &sigtrialIJ) const {
449  REAL X, L;
450  ComputeX(epsp + delepsp, X);
451  SolveL(X, L);
452  REAL res = FuncThetaL(ER, theta, L, sigtrialIJ);
453  return res;
454  }
455 
459  REAL FuncTheta2L(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec<REAL> &sigtrialIJ) const {
460  REAL I1 = sigtrialIJ[0];
461  REAL sqJ2 = sigtrialIJ[1];
462  REAL F;
463  ComputeF(L, F);
464  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
465  const REAL G = ER.Mu();
466  REAL y = (sqJ2 - F * sin(theta));
467  REAL x = G * fR / (9. * K)*(I1 - (L + F * fR * cos(theta)));
468  REAL res = x * sin(theta) - y * cos(theta);
469 #ifdef LOG4CXX
470  {
471  std::stringstream sout;
472  sout << "x = " << x << " y = " << y << " theta = " << theta << " res = " << res;
473  LOGPZ_DEBUG(loggerSM, sout.str())
474  }
475 #endif
476  return res;
477  }
478 
482  REAL FuncTheta2(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &sigtrialIJ) const {
483  REAL X, L;
484  ComputeX(epsp + delepsp, X);
485  SolveL(X, L);
486  return FuncTheta2L(ER, theta, L, sigtrialIJ);
487  }
488 
492  REAL DistThetaL(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec<REAL> &sigtrialIJ) const {
493  REAL I1 = sigtrialIJ[0];
494  REAL sqJ2 = sigtrialIJ[1];
495  REAL F;
496  ComputeF(L, F);
497  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
498  const REAL G = ER.Mu();
499  REAL y = (sqJ2 - F * sin(theta));
500  REAL x = (I1 - (L + F * fR * cos(theta)));
501  REAL dist = G * x * x / 2. + 9. * K * y * y / 2.;
502  return dist;
503  }
504 
508  REAL DistTheta(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &sigtrialIJ) const {
509  REAL X, L;
510  ComputeX(epsp + delepsp, X);
511  SolveL(X, L);
512  return DistThetaL(ER, theta, L, sigtrialIJ);
513  }
514 
518  void DFuncTheta(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &sigtrialIJ, TPZVec<REAL> &result) const {
519  REAL I1 = sigtrialIJ[0];
520  REAL sqJ2 = sigtrialIJ[1];
521  REAL X, L, DL, F, DF, Dtheta, Depsp;
522  ComputeX(epsp + delepsp, X);
523  LInitialGuess(X, L);
524  SolveL(X, L);
525  ComputeDL(L, epsp + delepsp, DL);
526  ComputeF(L, F);
527  ComputedF(L, DF);
528  DF *= DL;
529  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
530  const REAL G = ER.Mu();
531  REAL x = G * fR * (I1 - (L + F * fR * cos(theta)));
532  REAL y = 9. * K * (sqJ2 - F * sin(theta));
533  REAL dxepsp = -G * fR * (DL + DF * fR * cos(theta));
534  REAL dyepsp = -9. * K * DF * sin(theta);
535  REAL dxtheta = G * fR * F * fR * sin(theta);
536  REAL dytheta = -9. * K * F * cos(theta);
537  REAL denom = x * x + y*y;
538  Dtheta = 1 - (-y * dxtheta + x * dytheta) / denom;
539  Depsp = -(-y * dxepsp + x * dyepsp) / denom;
540  result[0] = Dtheta;
541  result[1] = Depsp;
542  }
543 
547  void DFuncTheta2(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec<REAL> &sigtrialIJ, TPZVec<REAL> &result) const {
548  REAL I1 = sigtrialIJ[0];
549  REAL sqJ2 = sigtrialIJ[1];
550  REAL X, L, DL, F, DF, Dtheta, Depsp;
551  ComputeX(epsp + delepsp, X);
552  LInitialGuess(X, L);
553  SolveL(X, L);
554  ComputeDL(L, epsp + delepsp, DL);
555  ComputeF(L, F);
556  ComputedF(L, DF);
557  DF *= DL;
558  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
559  const REAL G = ER.Mu();
560  REAL x = G * fR / (9. * K)*(I1 - (L + F * fR * cos(theta)));
561  REAL y = (sqJ2 - F * sin(theta));
562  REAL dxepsp = -G * fR / (9. * K)*(DL + DF * fR * cos(theta));
563  REAL dyepsp = -DF * sin(theta);
564  REAL dxtheta = G * fR / (9. * K) * F * fR * sin(theta);
565  REAL dytheta = -F * cos(theta);
566  Dtheta = dxtheta * sin(theta) - dytheta * cos(theta) + x * cos(theta) + y * sin(theta);
567  Depsp = dxepsp * sin(theta) - dyepsp * cos(theta);
568  result[0] = Dtheta;
569  result[1] = Depsp;
570  }
571 
575  void DFuncTheta2L(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec<REAL> &sigtrialIJ, TPZVec<REAL> &result) const {
576  REAL I1 = sigtrialIJ[0];
577  REAL sqJ2 = sigtrialIJ[1];
578  REAL DL, F, DF, Dtheta, Depsp;
579  DL = 1.;
580  ComputeF(L, F);
581  ComputedF(L, DF);
582  const REAL K = ER.Lambda() + 2. * ER.Mu() / 3.;
583  const REAL G = ER.Mu();
584  REAL x = G * fR / (9. * K)*(I1 - (L + F * fR * cos(theta)));
585  REAL y = (sqJ2 - F * sin(theta));
586  REAL dxepsp = -G * fR / (9. * K)*(DL + DF * fR * cos(theta));
587  REAL dyepsp = -DF * sin(theta);
588  REAL dxtheta = G * fR / (9. * K) * F * fR * sin(theta);
589  REAL dytheta = -F * cos(theta);
590  Dtheta = dxtheta * sin(theta) - dytheta * cos(theta) + x * cos(theta) + y * sin(theta);
591  Depsp = dxepsp * sin(theta) - dyepsp * cos(theta);
592  result[0] = Dtheta;
593  result[1] = Depsp;
594  }
595 
596  void UpdateSigtrialIJ(const TPZElasticResponse &ER, REAL epsp, REAL theta, TPZVec<REAL> &sigtrialIJ) {
597  REAL X, L;
598  ComputeX(epsp, X);
599  LInitialGuess(X, L);
600  SolveL(X, L);
601  UpdateSigtrialIJL(ER, L, theta, sigtrialIJ);
602  }
603 
604  void UpdateSigtrialIJL(const TPZElasticResponse &ER, REAL L, REAL theta, TPZVec<REAL> &sigtrialIJ) {
605  REAL F;
606  ComputeF(L, F);
607  sigtrialIJ[0] = L + F * fR * cos(theta);
608  sigtrialIJ[1] = F * sin(theta);
609  }
610 
611 protected:
613  void NewtonF1(const TPZElasticResponse &ER, REAL &L, TPZVec<REAL> &sigtrialIJ);
614 
615  void NewtonF2(const TPZElasticResponse &ER, REAL &epsp, TPZVec<REAL> &sigtrialIJ) {
616  REAL restheta, resdelepsp, disttheta;
617  REAL theta = 0.;
618  REAL delepsp = 0.;
619  disttheta = DistTheta(ER, theta, epsp, delepsp, sigtrialIJ);
620  // Look for a best guess for theta
621  for (REAL thetaguess = 0.; thetaguess <= M_PI; thetaguess += M_PI / 20.) {
622  REAL distnew = DistTheta(ER, thetaguess, epsp, delepsp, sigtrialIJ);
623  if (fabs(distnew) < fabs(disttheta)) {
624  theta = thetaguess;
625  disttheta = distnew;
626  }
627  }
628  REAL Xini, Lini, Fini;
629  ComputeX(epsp, Xini);
630  LInitialGuess(Xini, Lini);
631  SolveL(Xini, Lini);
632  ComputeF(Lini, Fini);
633  bool secondquadrant = false;
634  if (sigtrialIJ[0] < Lini) {
635  secondquadrant = true;
636  }
637  if (theta < M_PI / 2 && secondquadrant) {
638  theta = M_PI / 2 + M_PI / 20.;
639  }
640  if (Lini < -5.) {
641  REAL arcs = sigtrialIJ[1] / Fini;
642  if (arcs < 1.) {
643  theta = asin(arcs);
644  if (secondquadrant) {
645  theta = M_PI - theta;
646  }
647  }
648  REAL Lnew = sigtrialIJ[0] - Fini * fR * cos(theta);
649  REAL epsnew = ComputeEpsp(Lnew);
650  delepsp = epsnew - epsp;
651  }
652 
653  // perform Newton iterations
654  restheta = FuncTheta2(ER, theta, epsp, delepsp, sigtrialIJ);
655  resdelepsp = FuncEpsp(ER, theta, epsp, delepsp, sigtrialIJ);
656  REAL error = sqrt(restheta * restheta + resdelepsp * resdelepsp);
657  int64_t count = 0;
658  while ((fabs(restheta) > 1.e-10 || fabs(resdelepsp) > 1.e-10) && count < 100) {
659  REAL errprev = error;
660  TPZFNMatrix<4, REAL> tangent(2, 2);
661  TPZFNMatrix<2, REAL> resmat(2, 1);
662  TPZManVector<REAL, 2> tantheta(2, 0.), tandelepsp(2, 0.);
663  DFuncEpsp(ER, theta, epsp, delepsp, tandelepsp);
664  DFuncTheta2(ER, theta, epsp, delepsp, sigtrialIJ, tantheta);
665  for (int i = 0; i < 2; i++) {
666  tangent(1, i) = tandelepsp[i];
667  tangent(0, i) = tantheta[i];
668  }
669  resmat(0, 0) = restheta;
670  resmat(1, 0) = resdelepsp;
671 #ifdef LOG4CXX
672  {
673  std::stringstream sout;
674  tangent.Print("tangent matrix", sout);
675  resmat.Print("residual", sout);
676  LOGPZ_DEBUG(loggerSM, sout.str())
677  }
678 #endif
679  std::list<int64_t> singular;
680  tangent.SolveDirect(resmat, ELU, singular);
681  REAL scale = 1.;
682  if (epsp + delepsp - resmat(1, 0) < -fW) {
683  scale = 0.9999 * (fW + epsp + delepsp) / resmat(1, 0);
684  resmat *= scale;
685  }
686  REAL thetaprev = theta;
687  REAL delepspprev = delepsp;
688  theta = thetaprev - resmat(0, 0);
689  delepsp = delepspprev - resmat(1, 0);
690  restheta = FuncTheta2(ER, theta, epsp, delepsp, sigtrialIJ);
691  resdelepsp = FuncEpsp(ER, theta, epsp, delepsp, sigtrialIJ);
692  error = sqrt(restheta * restheta + resdelepsp * resdelepsp);
693  int iline = 0;
694  while (error > errprev && iline < 5) {
695  resmat *= 0.5;
696  theta = thetaprev - resmat(0, 0);
697  delepsp = delepspprev - resmat(1, 0);
698  restheta = FuncTheta2(ER, theta, epsp, delepsp, sigtrialIJ);
699  resdelepsp = FuncEpsp(ER, theta, epsp, delepsp, sigtrialIJ);
700  error = sqrt(restheta * restheta + resdelepsp * resdelepsp);
701  iline++;
702  }
703  count++;
704  }
705 #ifdef LOG4CXX
706  if (count > 10) {
707  std::stringstream sout;
708  sout << "iteration count " << count;
709  LOGPZ_DEBUG(loggerSM, sout.str())
710  }
711 #endif
712 
713 #ifdef LOG4CXX
714  {
715  std::stringstream sout;
716  sout << "sigtrialIJ input " << sigtrialIJ << " epsp input " << epsp << "\ndelepsp = " << delepsp
717  << " theta = " << theta;
718  LOGPZ_DEBUG(loggerSM, sout.str())
719  }
720 #endif
721  epsp += delepsp;
722  UpdateSigtrialIJ(ER, epsp, theta, sigtrialIJ);
723  }
724 
725  void NewtonF3(const TPZElasticResponse &ER, REAL &epsp, TPZVec<REAL> &sigtrialIJ) {
726  REAL restheta, resdelepsp, disttheta;
727  REAL theta = 0.;
728  REAL epspini = epsp;
729  REAL X, L, F;
730  ComputeX(epsp, X);
731  LInitialGuess(X, L);
732  SolveL(X, L);
733  ComputeF(L, F);
734  REAL Lini;
735  Lini = L;
736  disttheta = DistThetaL(ER, theta, L, sigtrialIJ);
737  // Look for a best guess for theta
738  for (REAL thetaguess = 0.; thetaguess <= M_PI; thetaguess += M_PI / 20.) {
739  REAL distnew = DistThetaL(ER, thetaguess, L, sigtrialIJ);
740  if (fabs(distnew) < fabs(disttheta)) {
741  theta = thetaguess;
742  disttheta = distnew;
743  }
744  }
745  bool secondquadrant = false;
746  if (sigtrialIJ[0] < L) {
747  secondquadrant = true;
748  }
749  if (theta < M_PI / 2 && secondquadrant) {
750  theta = M_PI / 2 + M_PI / 20.;
751  }
752  if (L < -5.) {
753  REAL arcs = sigtrialIJ[1] / F;
754  if (arcs < 1.) {
755  theta = asin(arcs);
756  if (secondquadrant) {
757  theta = M_PI - theta;
758  }
759  }
760  L = sigtrialIJ[0] - F * fR * cos(theta);
761  }
762 
763  // perform Newton iterations
764  restheta = FuncTheta2L(ER, theta, L, sigtrialIJ);
765  resdelepsp = FuncEpspUsingL(ER, theta, epspini, L, sigtrialIJ);
766  REAL error = sqrt(restheta * restheta + resdelepsp * resdelepsp);
767  int64_t count = 0;
768  while ((fabs(restheta) > 1.e-13 || fabs(resdelepsp) > 1.e-13) && count < 100) {
769  REAL errprev = error;
770  TPZFNMatrix<4, REAL> tangent(2, 2);
771  TPZFNMatrix<2, REAL> resmat(2, 1);
772  TPZManVector<REAL, 2> tantheta(2, 0.), tandelepsp(2, 0.);
773  DFuncEpspUsingL(ER, theta, L, sigtrialIJ, tandelepsp);
774  DFuncTheta2L(ER, theta, L, sigtrialIJ, tantheta);
775  for (int i = 0; i < 2; i++) {
776  tangent(1, i) = tandelepsp[i];
777  tangent(0, i) = tantheta[i];
778  }
779  resmat(0, 0) = restheta;
780  resmat(1, 0) = resdelepsp;
781 #ifdef LOG4CXX
782  {
783  std::stringstream sout;
784  tangent.Print("tangent matrix", sout);
785  resmat.Print("residual", sout);
786  LOGPZ_DEBUG(loggerSM, sout.str())
787  }
788 #endif
789  std::list<int64_t> singular;
790  tangent.SolveDirect(resmat, ELU, singular);
791  REAL thetaprev = theta;
792  REAL Lprev = L;
793  theta = thetaprev - resmat(0, 0);
794  L = Lprev - resmat(1, 0);
795  restheta = FuncTheta2L(ER, theta, L, sigtrialIJ);
796  resdelepsp = FuncEpspUsingL(ER, theta, epspini, L, sigtrialIJ);
797  error = sqrt(restheta * restheta + resdelepsp * resdelepsp);
798  int iline = 0;
799  while (error > errprev && iline < 5) {
800  resmat *= 0.5;
801  theta = thetaprev - resmat(0, 0);
802  L = Lprev - resmat(1, 0);
803  restheta = FuncTheta2L(ER, theta, L, sigtrialIJ);
804  resdelepsp = FuncEpspUsingL(ER, theta, epspini, L, sigtrialIJ);
805  error = sqrt(restheta * restheta + resdelepsp * resdelepsp);
806  iline++;
807  }
808  count++;
809  }
810 #ifdef LOG4CXX
811  if (count > 10) {
812  std::stringstream sout;
813  sout << "iteration count " << count;
814  LOGPZ_DEBUG(loggerSM, sout.str())
815  }
816 #endif
817 
818  epsp = ComputeEpsp(L);
819 #ifdef LOG4CXX
820  {
821  std::stringstream sout;
822  sout << "sigtrialIJ input " << sigtrialIJ << " epsp input " << epspini << "\ndelepsp = " << epsp - epspini
823  << " theta = " << theta << " Linitial " << Lini << " L " << L;
824  LOGPZ_DEBUG(loggerSM, sout.str())
825  }
826 #endif
827  UpdateSigtrialIJL(ER, L, theta, sigtrialIJ);
828  }
829 
830 protected:
831 
832  void NewtonF2L(const TPZElasticResponse &ER, REAL &L, TPZVec<REAL> &sigtrialIJ) {
833  REAL restheta, resdeltaL, disttheta;
834  REAL theta = 0.;
835  disttheta = DistThetaL(ER, theta, L, sigtrialIJ);
836  // Look for a best guess for theta
837  for (REAL thetaguess = 0.; thetaguess <= M_PI; thetaguess += M_PI / 20.) {
838  REAL distnew = DistThetaL(ER, thetaguess, L, sigtrialIJ);
839  if (fabs(distnew) < fabs(disttheta)) {
840  theta = thetaguess;
841  disttheta = distnew;
842  }
843  }
844  REAL Lini(L), Fini;
845  ComputeF(Lini, Fini);
846  bool secondquadrant = false;
847  if (sigtrialIJ[0] < Lini) {
848  secondquadrant = true;
849  }
850  if (theta < M_PI / 2 && secondquadrant) {
851  theta = M_PI / 2 + M_PI / 20.;
852  }
853  REAL deltaL = 0.;
854  // if (Lini < -5. || sigtrialIJ[0] < -5.) {
855  if (Lini * fB * log(fC) < -5. && sigtrialIJ[0] < 0.) {
856  REAL arcs = sigtrialIJ[1] / Fini;
857  if (arcs < 1.) {
858  theta = asin(arcs);
859  if (secondquadrant) {
860  theta = M_PI - theta;
861  }
862  }
863  L = sigtrialIJ[0] - Fini * fR * cos(theta);
864  deltaL = L - Lini;
865  }
866 
867  // perform Newton iterations
868  restheta = FuncTheta2L(ER, theta, L, sigtrialIJ);
869  resdeltaL = FuncEpspL(ER, theta, L, deltaL, sigtrialIJ);
870  REAL error = sqrt(restheta * restheta + resdeltaL * resdeltaL);
871  int64_t count = 0;
872  REAL diagTheta = 1., diagL = 1.;
873  REAL maxres = max(fabs(restheta / diagTheta), fabs(resdeltaL / diagL));
874  REAL maxresprev = maxres + 1.;
875  while ((maxres > 1.e-10 || maxres < maxresprev) && count < 150) {
876  maxresprev = maxres;
877  REAL errprev = error;
878  TPZFNMatrix<4, REAL> tangent(2, 2);
879  TPZFNMatrix<2, REAL> resmat(2, 1);
880  TPZManVector<REAL, 2> tantheta(2, 0.), tandeltaL(2, 0.);
881  DFuncEpspL(ER, theta, L, deltaL, sigtrialIJ, tandeltaL);
882  DFuncTheta2L(ER, theta, L, sigtrialIJ, tantheta);
883  for (int i = 0; i < 2; i++) {
884  tangent(1, i) = tandeltaL[i];
885  tangent(0, i) = tantheta[i];
886  }
887  diagTheta = tangent(0, 0);
888  diagL = tangent(1, 1);
889  resmat(0, 0) = restheta;
890  resmat(1, 0) = resdeltaL;
891 #ifdef LOG4CXX
892  {
893  std::stringstream sout;
894  tangent.Print("tangent matrix", sout);
895  resmat.Print("residual", sout);
896  LOGPZ_DEBUG(loggerSM, sout.str())
897  }
898 #endif
899  std::list<int64_t> singular;
900  tangent.SolveDirect(resmat, ELU, singular);
901  REAL thetaprev = theta;
902  REAL deltaLprev = deltaL;
903  theta = thetaprev - resmat(0, 0);
904  deltaL = deltaLprev - resmat(1, 0);
905  L = Lini + deltaL;
906  restheta = FuncTheta2L(ER, theta, L, sigtrialIJ);
907  resdeltaL = FuncEpspL(ER, theta, L, deltaL, sigtrialIJ);
908  error = sqrt(restheta * restheta + resdeltaL * resdeltaL);
909  int iline = 0;
910  while (error > errprev && iline < 5 && maxres > 1.e-10) {
911  resmat *= 0.5;
912  theta = thetaprev - resmat(0, 0);
913  deltaL = deltaLprev - resmat(1, 0);
914  L = Lini + deltaL;
915  restheta = FuncTheta2L(ER, theta, L, sigtrialIJ);
916  resdeltaL = FuncEpspL(ER, theta, L, deltaL, sigtrialIJ);
917  error = sqrt(restheta * restheta + resdeltaL * resdeltaL);
918  iline++;
919  }
920  maxres = max(fabs(restheta / diagTheta), fabs(resdeltaL / diagL));
921  count++;
922  if (count > 9 && maxres < 1.e-10) {
923  break;
924  }
925  }
926 #ifdef LOG4CXX
927  if (count > 10) {
928  std::stringstream sout;
929  sout << "iteration count " << count;
930  LOGPZ_DEBUG(loggerSM, sout.str())
931  }
932 #endif
933 
934 #ifdef LOG4CXX
935  {
936  std::stringstream sout;
937  sout << "sigtrialIJ input " << sigtrialIJ << " L input " << Lini << "\ndeltaL = " << deltaL
938  << " theta = " << theta;
939  LOGPZ_DEBUG(loggerSM, sout.str())
940  }
941 #endif
942  UpdateSigtrialIJL(ER, L, theta, sigtrialIJ);
943  }
944 public:
945 
946 
952  REAL fA, fB, fC;
953 
960  REAL fD, fW;
961 
967  REAL fR;
968 
972  bool fIsonCap;
973 
974 
975 public:
977 
981  inline int NumCases();
982 
983  static TPZTensor<REAL> gRefTension;
987  inline void LoadState(TPZFMatrix<REAL> &state);
988 
989  inline void ComputeTangent(TPZFMatrix<REAL> &tangent, TPZVec<REAL> &, int icase);
990 
991  inline void Residual(TPZFMatrix<REAL> &res, int icase);
992 
993  static void CheckConv();
994 
996 
998 
999  static void TestSolveL();
1001 
1002  static void McCormicRanchSand(TPZYCSandlerDimaggio & material);
1003 public:
1004 
1005 };
1006 
1007 template <class T>
1008 inline void TPZYCSandlerDimaggio::Compute(const TPZTensor<T> & sigma, const T & A, TPZVec<T> &res, int checkForcedYield) const {
1009  // the thermo force A in this case is assumed to be the
1010  // plastic volumetric strain itself. In fact it is not,
1011  // but the resultant derivatives are correct for practical purposes.
1012 
1013  // The following line evaluates L, the value of the first
1014  // invariant of stresses at the intersection of the
1015  // shear and hardening cap yield criteria / plastic potential.
1016  // It is first evaluated as REAL type to avoid unnecessary
1017  // derivatives evaluation.
1018 
1019 #ifdef LOG4CXX_PLASTICITY
1020  {
1021  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1022  std::stringstream sout;
1023  sout << ">>> TPZYCSandlerDimaggio::Compute *** - Plastic Potential / Yield - associative";
1024  LOGPZ_INFO(logger, sout.str().c_str());
1025  }
1026 #endif
1027 
1028  REAL L_REAL, X_REAL;
1029  ComputeX((REAL) TPZExtractVal::val(A), X_REAL);
1030  LInitialGuess(X_REAL, L_REAL);
1031  SolveL(X_REAL, L_REAL);
1032 
1033  T I1 = sigma.I1();
1034  T J2 = sigma.J2();
1035 
1036  // f1 - Modified Drucker-Prager as shear Yield Criterium
1037  T FI1;
1038  ComputeF(I1, FI1);
1039  if (fabs((REAL) TPZExtractVal::val(J2)) < 1.e-6) {
1040  res[0] = -FI1; // avoiding nan derivatives
1041  } else {
1042  res[0] = sqrt(J2) - FI1;
1043  }
1044 
1045  // f2 - ellipsoidal hardening/softening cap
1046  T FL, L(L_REAL), X;
1047  ComputeX(A, X);
1048  SolveL(X, L); // evaluating the derivatives of L
1049  ComputeF(L, FL);
1050 
1051  if (fabs((REAL) TPZExtractVal::val(FL)) < 0.00001) {
1052 #ifdef LOG4CXX_PLASTICITY
1053  {
1054  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1055  std::stringstream sout;
1056  sout << "*** TPZYCSandlerDimaggio::ComputePlasticPotential ***";
1057  sout << "\nDivision by F=" << TPZExtractVal::val(L) << " at f2 - ellipsoidal hardening/softening cap";
1058  LOGPZ_WARN(logger, sout.str().c_str());
1059  }
1060 #endif
1061  }
1062 
1063  T Temp1((L - I1) / (FL * T(fR)));
1064  Temp1 *= Temp1;
1065  T Temp2 = J2 / FL / FL;
1066 
1067  res[1] = Temp1 + Temp2 - T(1.);
1068 
1069  return;
1070 }
1071 
1072 template <class T>
1073 inline void TPZYCSandlerDimaggio::N(const TPZTensor<T> & sigma, const T & A, TPZVec<TPZTensor<T> > & Ndir, int checkForcedYield) const {
1074 
1075  // the thermo force A in this case is assumed to be the
1076  // plastic volumetric strain itself. In fact it is not,
1077  // but the resultant derivatives are correct for practical purposes.
1078 
1079  //The following line evaluates L, the value of the first
1080  // invariant of stresses at the intersection of the
1081  // shear and hardening cap yield criteria / plastic potential.
1082  // It is first evaluated as REAL type to avoid unnecessary
1083  // derivatives evaluation.
1084 
1085 
1086 #ifdef LOG4CXX_PLASTICITY
1087  {
1088  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1089  std::stringstream sout;
1090  sout << ">>> TPZYCSandlerDimaggio::N *** - Plastification direction - associative";
1091  LOGPZ_INFO(logger, sout.str().c_str());
1092  }
1093 #endif
1094 
1095  REAL ResTol = 1.e-6;
1096 
1097  REAL L_REAL, X_REAL;
1098  ComputeX((REAL) TPZExtractVal::val(A), X_REAL);
1099  LInitialGuess(X_REAL, L_REAL);
1100  SolveL(X_REAL, L_REAL, ResTol);
1101 
1102  T I1 = sigma.I1();
1103  T J2 = sigma.J2();
1104  T SQRTJ2;
1105  if (TPZExtractVal::val(J2) > 1.e-12) {
1106  SQRTJ2 = sqrt(J2);
1107  } else {
1108  SQRTJ2 = 1.e-6;
1109  }
1110 
1111  {
1112  //f1 - Modified Drucker-Prager as shear Yield Criterium / Plastic Potential
1113  REAL fz = FZero();
1114  T Temp1(0.);
1115  if (TPZExtractVal::val(I1) < fz) {
1116  Temp1 = I1 * T(fB);
1117  Temp1 = exp(Temp1) * T(fB * fC);
1118 
1119  if ((REAL) TPZExtractVal::val(SQRTJ2) < 1.e-6) // just for robustness. f1 shouldn't be reached when J2 = 0.
1120  {
1121 #ifdef LOG4CXX_PLASTICITY
1122  {
1123  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1124  std::stringstream sout;
1125  sout << "*** TPZYCSandlerDimaggio::N *** - SQRT(J2) = " << TPZExtractVal::val(SQRTJ2) << " < 1.e-6 causes error in 0-th yield function. Imposing J2 = 1.e-6 instead";
1126  LOGPZ_WARN(logger, sout.str().c_str());
1127  }
1128 #endif
1129  SQRTJ2 = T(1.e-6);
1130  }
1131  }
1132  Temp1 = Temp1 - I1 / SQRTJ2 / T(6.);
1133  T Temp2 = T(1.) / SQRTJ2;
1134  T Temp3 = Temp2 / T(2.);
1135 
1136  Ndir[0].XX() = Temp1 + sigma.XX() * Temp3;
1137  Ndir[0].YY() = Temp1 + sigma.YY() * Temp3;
1138  Ndir[0].ZZ() = Temp1 + sigma.ZZ() * Temp3;
1139  // Ndir[0].YZ() = sigma.YZ() * Temp2;
1140  // Ndir[0].XZ() = sigma.XZ() * Temp2;
1141  // Ndir[0].XY() = sigma.XY() * Temp2;
1142  Ndir[0].YZ() = sigma.YZ() * Temp3;
1143  Ndir[0].XZ() = sigma.XZ() * Temp3;
1144  Ndir[0].XY() = sigma.XY() * Temp3;
1145  }
1146 
1147  {//f2 - ellipsoidal hardening/softening cap
1148 
1149  T FL, X, L(L_REAL * 1. - ResTol); // guaranteeing that the function will be evaluated
1150  ComputeX(A, X);
1151  SolveL(X, L, ResTol); // evaluating the derivatives of L
1152 
1153  ComputeF(L, FL);
1154  // the radius of the ellips needs to be positive
1155  // this should be taken care of by the computation of L which is limited by LMax()
1156  if (TPZExtractVal::val(FL) <= 0.) {
1157  DebugStop();
1158  }
1159  T FL2 = FL * FL;
1160  T FL3 = FL2; // / T(2.);
1161 
1162  T Temp = (I1 - L) / T(fR * fR) - I1 / T(6.);
1163  Temp = Temp / FL2 * T(2.);
1164 
1165 #ifdef LOG4CXX_PLASTICITY
1166  {
1167  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1168  std::stringstream sout;
1169  sout << "*** TPZYCSandlerDimaggio::N *** X = " << X
1170  << "\n L = " << L << " L_REAL = " << L_REAL
1171  << "\n FL = " << FL
1172  << "\n Temp = " << Temp;
1173  LOGPZ_DEBUG(logger, sout.str().c_str());
1174  }
1175 #endif
1176 
1177  Ndir[1].XX() = Temp + sigma.XX() / FL2;
1178  Ndir[1].YY() = Temp + sigma.YY() / FL2;
1179  Ndir[1].ZZ() = Temp + sigma.ZZ() / FL2;
1180  Ndir[1].YZ() = sigma.YZ() / FL3;
1181  Ndir[1].XZ() = sigma.XZ() / FL3;
1182  Ndir[1].XY() = sigma.XY() / FL3;
1183  }
1184 
1185 #ifdef LOG4CXX
1186  {
1187  LoggerPtr logger(Logger::getLogger("pz.plasticity.SandlerDimaggio.main"));
1188  std::stringstream sout;
1189  sout << "<< TPZYCSandlerDimaggio::N *** \n sigma = \n" << sigma
1190  << "\nI1 = " << I1
1191  << "\nJ2 = " << J2
1192  << "\nSQRTJ2 = " << SQRTJ2
1193  << "\nNdir = \n" << Ndir;
1194  //LOGPZ_DEBUG(logger,sout.str().c_str());
1195  }
1196 #endif
1197 
1198  return;
1199 }
1200 
1201 template <class T>
1202 inline void TPZYCSandlerDimaggio::H(const TPZTensor<T> & sigma, const T & A, TPZVec<T> & h, int checkForcedYield) const {
1203 
1204  // the thermo force A in this case is assumed to be the
1205  // plastic volumetric strain itself. In fact it is not,
1206  // but the resultant derivatives are correct for practical purposes.
1207 
1208  //The following line evaluates L, the value of the first
1209  // invariant of stresses at the intersection of the
1210  // shear and hardening cap yield criteria / plastic potential.
1211  // It is first evaluated as REAL type to avoid unnecessary
1212  // derivatives evaluation.
1213 
1214 #ifdef LOG4CXX_PLASTICITY
1215  {
1216  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1217  std::stringstream sout;
1218  sout << ">>> TPZYCSandlerDimaggio::H *** - Hardening modulus";
1219  LOGPZ_INFO(logger, sout.str().c_str());
1220  }
1221 #endif
1222 
1223  REAL L_REAL, X_REAL;
1224  ComputeX((REAL) TPZExtractVal::val(A), X_REAL);
1225  LInitialGuess(X_REAL, L_REAL);
1226  SolveL(X_REAL, L_REAL);
1227 
1228  T I1 = sigma.I1();
1229 
1230  {//f1 - Modified Drucker-Prager as shear Yield Criterium / Plastic Potential
1231 
1232  h[0] = exp(I1 * T(fB)) * T(3. * fB * fC);
1233  // h > 0 because plastic deformation is dilatant
1234  }
1235 
1236  {//f2 - ellipsoidal hardening/softening cap
1237  T FL, X, L(L_REAL);
1238  ComputeX(A, X);
1239  SolveL(X, L); // evaluating the derivatives of L
1240  ComputeF(L, FL);
1241  T FL2 = FL * FL;
1242 
1243  h[1] = (I1 - L) / FL2 * T(6. / fR / fR);
1244  // h <=0 because plastic deformation is compactant
1245  }
1246 
1247  return;
1248 
1249 }
1250 
1251 template <class T>
1252 inline void TPZYCSandlerDimaggio::SolveL(const T & X, T & L, REAL relTol) const {
1253  T F, dF, res, dRes;
1254 
1255  ComputeF(L, F);
1256  res = F * T(fR) + X - L;
1257 
1258  int i = 0; // ensuring evaluating at least once the Newton method
1259  // so that the evaluation with FAD Type and converged L evaluates the
1260  // derivatives
1261  while (
1262  fabs((REAL) TPZExtractVal::val(res)) > relTol || i < 1) {
1263  i++;
1264 
1265  ComputedF(L, dF);
1266  dRes = dF * T(fR) - T(1.);
1267 
1268  L -= res / dRes;
1269 
1270  ComputeF(L, F);
1271  res = F * T(fR) + X - L;
1272  }
1273  if (TPZExtractVal::val(L) > LMax()) {
1274  L = T(LMax());
1275  }
1276 
1277 }
1278 
1279 template <class T>
1280 inline void TPZYCSandlerDimaggio::LInitialGuess(const T & X, T & L) const {
1281  T FAprox;
1282  ComputeF(X / T(2.), FAprox);
1283  L = X + FAprox * T(fR);
1284 }
1285 
1291 template <class T>
1292 inline void TPZYCSandlerDimaggio::ComputeDL(const T &L, const T &A, T &DL) const {
1293  REAL LMx = LMax();
1294  if (TPZExtractVal::val(L) >= LMx) {
1295  DL = T(0.);
1296  } else {
1297  DL = T(1.) / (T(fD)*(A + T(fW))*(T(1.) + T(fR * fB * fC) * exp(T(fB) * L)));
1298  }
1299 }
1300 
1301 template <class T>
1302 inline void TPZYCSandlerDimaggio::ComputeF(const T & L, T & F) const {
1303  F = T(fA) - exp(L * T(fB)) * T(fC);
1304 }
1305 
1306 template <class T>
1307 inline void TPZYCSandlerDimaggio::ComputedF(const T & L, T & dF) const {
1308  dF = exp(L * T(fB)) * T(-fC * fB);
1309 }
1310 
1311 inline void TPZYCSandlerDimaggio::ComputeD2F(const REAL L, REAL & d2F) const {
1312  d2F = exp(L * fB) * (-fC * fB * fB);
1313 }
1314 
1315 template <class T>
1316 inline void TPZYCSandlerDimaggio::ComputeX(const T & A, T & X) const {
1317  REAL ep = -0.999999 * fW;
1318 
1319  if (TPZExtractVal::val(A) < ep) {
1320  T dXdep = T(fW / (ep + fW) / fD);
1321  X = T(log(ep / fW + 1.) / fD);
1322  X = X + dXdep * (A - T(ep));
1323 #ifdef LOG4CXX_PLASTICITY
1324  {
1325  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1326  std::stringstream sout;
1327  sout << "*** TPZYCSandlerDimaggio::ComputeX *** ##### Changing hardening equation to adjusted linear - Excessive volumetric plastic strain - Please check aftwerwards if alpha = epsP.I1() to verify if results are consistent! ######";
1328  LOGPZ_WARN(logger, sout.str().c_str());
1329  }
1330 #endif
1331  } else {
1332  X = A / T(fW) + T(1.);
1333  if (TPZExtractVal::val(X) <= 0.) {
1334  DebugStop();
1335  }
1336  X = log(X) / T(fD);
1337  }
1338 
1339 }
1340 
1342 
1343 inline int TPZYCSandlerDimaggio::NumCases() {
1344  return 4;
1345 }
1346 
1347 inline void TPZYCSandlerDimaggio::LoadState(TPZFMatrix<REAL> &state) {
1348 
1349  int i;
1350  for (i = 0; i < 6; i++) gRefTension.fData[i] = state(i, 0);
1351 
1352 #ifdef LOG4CXX_PLASTICITY
1353  {
1354  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1355 
1356  std::stringstream sout;
1357  sout << ">>> LoadState *** Tension " << state;
1358  LOGPZ_DEBUG(logger, sout.str().c_str());
1359  }
1360 #endif
1361 }
1362 
1363 inline void TPZYCSandlerDimaggio::ComputeTangent(TPZFMatrix<REAL> &tangent, TPZVec<REAL> &, int icase) {
1364 
1365  const int nVars = 6;
1366  typedef TFad<nVars, REAL> TFAD;
1367 
1368  int i, j;
1369  TPZVec< TPZTensor < REAL > > N_Dir(2);
1370  REAL A = -0.05; // example epsVP value to be used with checkconv
1371  TPZTensor < TFAD > Sigma_FAD;
1372  TFAD A_FAD(A);
1373  TPZVec< TPZTensor < TFAD > > N_Dir_FAD(2);
1374 
1375  switch (icase) {
1376  case 0:
1377  case 1:
1378  //Compute N
1379  tangent.Redim(1, nVars);
1380  N(gRefTension, A, N_Dir, 0);
1381  for (i = 0; i < nVars; i++)
1382  tangent(0, i) = N_Dir[icase].fData[i];
1383  break;
1384  case 2:
1385  case 3:
1386  //Compute derivatives of N
1387  tangent.Redim(nVars, nVars);
1388  gRefTension.CopyTo(Sigma_FAD);
1389  for (i = 0; i < nVars; i++)
1390  Sigma_FAD.fData[i].diff(i, nVars);
1391  N(Sigma_FAD, A_FAD, N_Dir_FAD, 0);
1392  for (i = 0; i < nVars; i++)
1393  for (j = 0; j < nVars; j++)
1394  tangent(i, j) = N_Dir_FAD[icase - 2].fData[i].dx(j);
1395  break;
1396 
1397  }
1398 #ifdef LOG4CXX_PLASTICITY
1399  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1400  std::stringstream sout;
1401  sout << ">>> ComputeTangent *** " << tangent;
1402  LOGPZ_DEBUG(logger, sout.str().c_str());
1403 #endif
1404 }
1405 
1406 inline void TPZYCSandlerDimaggio::Residual(TPZFMatrix<REAL> &res, int icase) {
1407 
1408  int i;
1409  const int nVars = 6;
1410  // typedef TFad<nVars,REAL> TFAD;
1411 
1412  TPZVec< REAL > PlasticPot(2);
1413  REAL A = -0.05; // example epsVP value to be used with checkconv
1414  TPZVec< TPZTensor<REAL> > N_Dir(2);
1415  switch (icase) {
1416  case 0:
1417  case 1:
1418  //Compute PlasticPotential
1419  res.Redim(1, 1);
1420  Compute(gRefTension, A, PlasticPot, 0);
1421  res(0, 0) = PlasticPot[icase];
1422  break;
1423  case 2:
1424  case 3:
1425  //Compute Ndir
1426  res.Redim(nVars, 1);
1427  N(gRefTension, A, N_Dir, 0);
1428  for (i = 0; i < nVars; i++)
1429  res(i, 0) = N_Dir[icase - 2].fData[i];
1430  break;
1431  }
1432 
1433 #ifdef LOG4CXX_PLASTICITY
1434  LoggerPtr logger(Logger::getLogger("plasticity.SandlerDimaggio"));
1435  std::stringstream sout;
1436  sout << ">>> Residual *** " << res;
1437  LOGPZ_DEBUG(logger, sout.str().c_str());
1438 #endif
1439 
1440 }
1441 
1442 inline void TPZYCSandlerDimaggio::CheckConv() {
1443 
1444  const int nVars = 6;
1445 
1446  // Creating the Sandler Dimaggio obejct
1447  TPZYCSandlerDimaggio YCSandlerDimaggio;
1448 
1449  TPZYCSandlerDimaggio::McCormicRanchSand(YCSandlerDimaggio);
1450 
1451  REAL multipl = 1.;
1452 
1453  TPZFNMatrix<nVars> sigma(nVars, 1), Range(nVars, 1);
1454  sigma(_XX_, 0) = -0.17 * multipl;
1455  sigma(_YY_, 0) = -0.13 * multipl;
1456  sigma(_ZZ_, 0) = -0.11 * multipl;
1457  sigma(_XY_, 0) = -0.7 * multipl;
1458  sigma(_XZ_, 0) = -0.5 * multipl;
1459  sigma(_ZZ_, 0) = -0.3 * multipl;
1460 
1461  Range = sigma * (1. / 19.);
1462  TPZVec< REAL > Coefs(1, 1.);
1463  CheckConvergence(YCSandlerDimaggio, sigma, Range, Coefs);
1464 
1465 }
1466 
1468 
1470 
1471 inline void TPZYCSandlerDimaggio::TestSolveL() {
1472 #ifdef LOG4CXX_PLASTICITY
1473  LoggerPtr loggerSandlerDimaggio(Logger::getLogger("plasticity.SandlerDimaggio"));
1474  {
1475  std::stringstream sout;
1476  sout << "<<< TPZYCSandlerDimaggio::TestSolveL ***";
1477  LOGPZ_INFO(loggerSandlerDimaggio, sout.str().c_str());
1478  }
1479 #endif
1480 
1481  const int oneVar = 1;
1482  typedef TFad<oneVar, REAL> TFAD_ONE;
1483 
1484  // Creating the Sandler Dimaggio obejct
1485  TPZYCSandlerDimaggio YCSandlerDimaggio;
1486 
1487  TPZYCSandlerDimaggio::McCormicRanchSand(YCSandlerDimaggio);
1488 
1489  // Verifying if the SolveL routines are working
1490  // It needs a correct evaluation of both F and dF
1491 
1492  //Supposing a plastic volumetric strain
1493  REAL X, L, dF, epsVP = -.05;
1494  TFAD_ONE L_FAD, F_FAD;
1495  YCSandlerDimaggio.ComputeX(epsVP, X);
1496 
1497  YCSandlerDimaggio.LInitialGuess(X, L);
1498  L_FAD = L;
1499  L_FAD.diff(0, 0);
1500 
1501  YCSandlerDimaggio.ComputeF(L_FAD, F_FAD);
1502  YCSandlerDimaggio.ComputedF((REAL) TPZExtractVal::val(L_FAD), dF);
1503 
1504 #ifdef LOG4CXX_PLASTICITY
1505  {
1506  std::stringstream sout;
1507  sout << "*** TPZYCSandlerDimaggio::TestSolveL ***";
1508  sout << "\nDerivative evaluated with FAD (dF/dL)FAD = " << F_FAD.dx(0);
1509  sout << "\nDerivative evaluated explicitly (dF/dL) = " << dF;
1510  sout << "\nProposed L = " << L;
1511  LOGPZ_DEBUG(loggerSandlerDimaggio, sout.str().c_str());
1512  }
1513 #endif
1514 
1515  YCSandlerDimaggio.SolveL(X, L);
1516 #ifdef LOG4CXX_PLASTICITY
1517  {
1518  std::stringstream sout;
1519  sout << "*** TPZYCSandlerDimaggio::TestSolveL ***";
1520  sout << "\nSolved L = " << L;
1521  LOGPZ_DEBUG(loggerSandlerDimaggio, sout.str().c_str());
1522  }
1523 #endif
1524 
1525  //REAL multipl = 1.; // testing the shear yield criterium (f1)
1526  REAL multipl = 10.; // testing the elipsoidal compressive cap (f2)
1527 
1528  // Checking if NDir:(1,1,1,0,0,0) equals H
1529  // verifying if the TFAD derivatives of N equal the Ndir vector.
1530  const int sixVars = 6;
1531  typedef TFad<sixVars, REAL> TFAD_SIX;
1532 
1533  TPZTensor<REAL> sigma;
1534  TPZVec<TPZTensor<REAL> > Ndir(2);
1535  TPZVec<REAL> h(2);
1536  TPZTensor<TFAD_SIX> sigma_FAD;
1537  TFAD_SIX epsVP_FAD(epsVP);
1538  TPZVec< TFAD_SIX > PlasticPot(2);
1539 
1540  sigma.XX() = -0.17 * multipl;
1541  sigma.YY() = -0.13 * multipl;
1542  sigma.ZZ() = -0.11 * multipl;
1543  sigma.XY() = -0.7 * multipl;
1544  sigma.YZ() = -0.5 * multipl;
1545  sigma.XZ() = -0.3 * multipl;
1546 
1547  sigma_FAD.XX() = sigma.XX();
1548  sigma_FAD.XX().diff(0, sixVars);
1549 
1550  sigma_FAD.YY() = sigma.YY();
1551  sigma_FAD.YY().diff(3, sixVars);
1552 
1553  sigma_FAD.ZZ() = sigma.ZZ();
1554  sigma_FAD.ZZ().diff(5, sixVars);
1555 
1556  sigma_FAD.XY() = sigma.XY();
1557  sigma_FAD.XY().diff(1, sixVars);
1558 
1559  sigma_FAD.YZ() = sigma.YZ();
1560  sigma_FAD.YZ().diff(4, sixVars);
1561 
1562  sigma_FAD.XZ() = sigma.XZ();
1563  sigma_FAD.XZ().diff(2, sixVars);
1564 
1565  YCSandlerDimaggio.N(sigma, epsVP, Ndir, 0);
1566  YCSandlerDimaggio.H(sigma, epsVP, h, 0);
1567  YCSandlerDimaggio.Compute(sigma_FAD, epsVP_FAD, PlasticPot, 0);
1568 
1569 #ifdef LOG4CXX_PLASTICITY
1570  {
1571  std::stringstream sout;
1572  sout << "*** TPZYCSandlerDimaggio::TestSolveL ***";
1573  sout << "\nVerifying if H equals depsVP/dSigma = Ndir:(1 1 1 0 0 0)";
1574 
1575  for (int i = 0; i < NYield; i++) {
1576 
1577  sout << "\n" << i << "-th Yield Criterium";
1578  sout << "\nepsVP calculated from Ndir = " << Ndir[i].XX() + Ndir[i].YY() + Ndir[i].ZZ();
1579  sout << "\nepsVP calculated from H = " << h[i];
1580  sout << "\nVerifying if dPlasticPot/dSigma equals Ndir";
1581  sout << "\nNdir evaluated explicitly with function N():" << Ndir[i];
1582  sout << "\nN evaluated through FAD evaluations: dPlasticPot/dsigma:" << PlasticPot[i];
1583 
1584  }
1585 
1586  LOGPZ_DEBUG(loggerSandlerDimaggio, sout.str().c_str());
1587  }
1588 #endif
1589 
1590 }
1592 
1593 inline void TPZYCSandlerDimaggio::McCormicRanchSand(TPZYCSandlerDimaggio & material) {
1594 #ifdef LOG4CXX_PLASTICITY
1595  LoggerPtr loggerSandlerDimaggio(Logger::getLogger("plasticity.SandlerDimaggio"));
1596  {
1597  std::stringstream sout;
1598  sout << ">>> TPZYCSandlerDimaggio::McCormicRanchSand ***";
1599  LOGPZ_INFO(loggerSandlerDimaggio, sout.str().c_str());
1600  }
1601 #endif
1602  // setup with data from McCormic Ranch Sand
1603  // Dimaggio, Frank L. Sandler, Ivan S. Material model for granular soils
1604  // J. of the Eng. Mech. Div. vol. 97 n0 EM3
1605  // pp 935-949,1971
1606  // OBS: stresses in ksi
1607  REAL A = 0.25,
1608  B = 0.67,
1609  C = 0.18,
1610  D = 0.67,
1611  R = 2.5,
1612  W = 0.066;
1613 
1614  material.SetUp(A, B, C, D, R, W);
1615 }
1616 
1627 inline void TPZYCSandlerDimaggio::InitialGuess(const TPZElasticResponse &ER, REAL epsp, TPZTensor<REAL> &sigtrial, REAL &epspproj,
1628  TPZVec<REAL> &delgamma, TPZTensor<REAL> &sigproj) {
1629  TPZManVector<REAL, 2> yield(2, 0.);
1630  Compute(sigtrial, epsp, yield, 0);
1631  if (yield[0] <= 0. && yield[1] <= 0.) {
1632  epspproj = epsp;
1633  sigproj = sigtrial;
1634  delgamma.Fill(0.);
1635 #ifdef LOG4CXX
1636  if (loggerSM->isDebugEnabled()) {
1637  std::stringstream sout;
1638  sout << "Deformation elastic yield = " << yield;
1639  sout << "delgamma condition " << (delgamma[0] > 0.) << " " << (delgamma[1] > 0.) << std::endl;
1640  LOGPZ_DEBUG(loggerSM, sout.str())
1641  }
1642 #endif
1643  return;
1644  }
1645  TPZTensor<REAL> S;
1646  sigtrial.S(S);
1647  REAL J2 = S.J2();
1648  REAL sqJ2 = sqrt(fabs(J2));
1649  REAL I1 = sigtrial.I1();
1650  if (yield[1] > 0.) {
1651  // project the point on surface F2
1652  TPZManVector<REAL, 2> sigtrialIJ(2);
1653  sigtrialIJ[0] = I1;
1654  sigtrialIJ[1] = sqJ2;
1655  REAL epspprev = epsp;
1656  NewtonF3(ER, epsp, sigtrialIJ);
1657  sigtrial.Adjust(sigtrialIJ, sigproj);
1658  epspproj = epsp;
1659  TPZManVector<TPZTensor<REAL>, 2> Ndir(2);
1660  this->N(sigproj, epspproj, Ndir, 1);
1661  TPZTensor<REAL> sigPlast(sigtrial), epsPlast;
1662  sigPlast.Add(sigproj, -1.);
1663  ER.ComputeStrain(sigPlast, epsPlast);
1664  REAL scale = epsPlast.Norm() / Ndir[1].Norm();
1665  for (int i = 0; i < 6; i++) {
1666  REAL diff = fabs(scale * Ndir[1][i] - epsPlast[i]);
1667  if (diff > 1.e-6) {
1668  epsp = epspprev;
1669  sigtrialIJ[0] = I1;
1670  sigtrialIJ[1] = sqJ2;
1671 
1672  NewtonF3(ER, epsp, sigtrialIJ);
1673 
1674  // DebugStop();
1675  }
1676  }
1677  delgamma[0] = 0.;
1678  delgamma[1] = scale;
1679  }
1680  Compute(sigproj, epspproj, yield, 0);
1681 #ifdef LOG4CXX
1682  if (loggerSM->isDebugEnabled()) {
1683  std::stringstream sout;
1684  sout << "After projecting the point yield = " << yield;
1685  sout << "\ndelgamma = " << delgamma;
1686  LOGPZ_DEBUG(loggerSM, sout.str())
1687  }
1688 #endif
1689  if (yield[0] > 0.) {
1690  ;
1691  // DebugStop();
1692  }
1693  if (yield[1] > 1.e-6) {
1694  DebugStop();
1695  }
1696  // falta calcular delgamma
1697 }
1698 
1702 template<class T>
1703 void TPZYCSandlerDimaggio::EpspFromL(const T &L, T &epsp) const {
1704  T FL;
1705  ComputeF(L, FL);
1706  T fac1 = T(fD)*(L - FL * T(fR));
1707  epsp = (exp(fac1) - T(1.)) * T(fW);
1708 }
1709 
1713 template<class T>
1714 void TPZYCSandlerDimaggio::DEpspDL(const T &L, T& depspdL) const {
1715  T FL;
1716  ComputeF(L, FL);
1717  T fac1 = T(fD)*(L - FL * T(fR));
1718  T fac2 = (T(1.) + T(fB * fC * fR) * exp(T(fB) * L));
1719  depspdL = T(fD * fW) * exp(fac1) * fac2;
1720 }
1721 
1725 template<class T>
1726 void TPZYCSandlerDimaggio::D2EpspDL2(const T &L, T& d2epspdL2) const {
1727  T FL;
1728  ComputeF(L, FL);
1729  T fac1 = T(fD)*(L - FL * T(fR));
1730  T BCR = T(fB * fC * fR) * exp(T(fB) * L);
1731  T fac2 = T(1.) + BCR;
1732  T fac3 = T(fB) * BCR + T(fD) * fac2*fac2;
1733  d2epspdL2 = T(fD * fW) * exp(fac1) * fac3;
1734 
1735 }
1736 
1738 
1739 inline void TPZYCSandlerDimaggio::NewtonF1(const TPZElasticResponse &ER, REAL &L, TPZVec<REAL> &sigtrialIJ) {
1740 
1741  REAL resultL = L;
1742  TPZManVector<STATE, 2> sigProj(sigtrialIJ);
1743  REAL F;
1744  ComputeF(sigProj[0], F);
1745  if (F > 0.) {
1746  resultL = sigProj[0];
1747  }
1748  REAL ddist = DDistance(ER, resultL, sigProj);
1749  while (ddist < 0.) {
1750  resultL += 1.;
1751  ddist = DDistance(ER, resultL, sigProj);
1752  }
1753  REAL ddistnext = fabs(ddist) + 1;
1754  REAL correct = 1.;
1755  int count = 0;
1756  while ((count < 20) && ((fabs(ddist) > 1.e-10 && fabs(correct) > 1.e-14) || (fabs(ddistnext) > fabs(ddist)))) {
1757  ddist = ddistnext;
1758  REAL d2dist = D2Distance(ER, resultL, sigProj);
1759  correct = ddist / d2dist;
1760 
1761  resultL -= correct;
1762  if (fabs(resultL) > 1.) {
1763  correct /= resultL;
1764  }
1765 
1766  ddistnext = DDistance(ER, resultL, sigProj);
1767 
1768  count++;
1769  }
1770  sigProj[0] = resultL;
1771  ComputeF(resultL, F);
1772  sigProj[1] = F;
1773  // compute the increase in epsp
1774  // compute L such that depsp/dL (delta L) = delta epsp
1775  REAL K = ER.K();
1776  STATE depspdl;
1777  DEpspDL(resultL, depspdl);
1778  STATE residueL;
1779 #ifdef PZDEBUG_KEEP
1780  TPZFMatrix<STATE> table(2, 200, 0.);
1781  int64_t count = 0;
1782  for (resultL = -1.; resultL < 10.; resultL += 0.1) {
1783  residueL = 3. * K * depspdl * (resultL - L)-(sigtrialIJ[0] - sigProj[0]);
1784  table(0, count) = resultL;
1785  table(1, count) = residueL;
1786  count++;
1787  }
1788  table.Resize(2, count);
1789  {
1790  std::ofstream fout("resfunc.nb");
1791  table.Print("resfunc", fout, EMathematicaInput);
1792  }
1793 #endif
1794  resultL = L;
1795  DEpspDL(resultL, depspdl);
1796  residueL = 3. * K * depspdl * (resultL - L)-(sigtrialIJ[0] - sigProj[0]);
1797  while (residueL < 0.) {
1798  resultL += 1.;
1799  DEpspDL(resultL, depspdl);
1800  residueL = 3. * K * depspdl * (resultL - L)-(sigtrialIJ[0] - sigProj[0]);
1801  }
1802  while (fabs(residueL) > 1.e-10) {
1803  STATE d2epspdl2;
1804  D2EpspDL2(resultL, d2epspdl2);
1805  STATE dres = 3. * K * depspdl + 3. * K * (resultL - L) * d2epspdl2;
1806  resultL -= residueL / dres;
1807  DEpspDL(resultL, depspdl);
1808  residueL = 3. * K * depspdl * (resultL - L)-(sigtrialIJ[0] - sigProj[0]);
1809  }
1810  sigtrialIJ = sigProj;
1811  L = resultL;
1812 }
1813 
1817 inline REAL TPZYCSandlerDimaggio::Distance(const TPZElasticResponse &ER, REAL L, TPZVec<REAL> &sigtrialIJ) const {
1818  REAL I1 = sigtrialIJ[0];
1819  REAL sqj2 = sigtrialIJ[1];
1820  REAL x = L;
1821  REAL y;
1822  ComputeF(L, y);
1823  REAL delx = x - I1;
1824  REAL dely = y - sqj2;
1825  REAL dist = ER.G() * delx * delx / 2. + dely * dely * 9. * ER.K() / 2.;
1826  return dist;
1827 }
1828 
1832 inline REAL TPZYCSandlerDimaggio::DDistance(const TPZElasticResponse &ER, REAL L, TPZVec<REAL> &sigtrialIJ) const {
1833  REAL I1 = sigtrialIJ[0];
1834  REAL sqj2 = sigtrialIJ[1];
1835  REAL x = L;
1836  REAL y, dy;
1837  ComputeF(L, y);
1838  ComputedF(L, dy);
1839  REAL delx = x - I1;
1840  REAL dely = y - sqj2;
1841  REAL ddely = dy;
1842  REAL ddelx = 1;
1843  REAL ddist = ER.G() * ddelx * delx + ddely * dely * 9. * ER.K();
1844  return ddist;
1845 }
1846 
1850 inline REAL TPZYCSandlerDimaggio::D2Distance(const TPZElasticResponse &ER, REAL L, TPZVec<REAL> &sigtrialIJ) const {
1851  // REAL I1 = sigtrialIJ[0];
1852  REAL sqj2 = sigtrialIJ[1];
1853  // REAL x = L;
1854  REAL y, dy, d2y;
1855  ComputeF(L, y);
1856  ComputedF(L, dy);
1857  ComputeD2F(L, d2y);
1858  // REAL delx = x-I1;
1859  REAL dely = y - sqj2;
1860  REAL ddely = dy;
1861  REAL ddelx = 1;
1862  REAL d2dist = ER.G() * ddelx * ddelx + ddely * ddely * 9. * ER.K() + d2y * dely * 9. * ER.K();
1863  return d2dist;
1864 
1865 }
1866 
1867 
1868 
1869 #endif //TPZYCSANDLERDIMAGGIO_H
TPZYCSandlerDimaggio::GetNYield
virtual int GetNYield() const override
Definition: TPZYCSandlerDimaggio.h:225
_XZ_
#define _XZ_
Definition: TPZTensor.h:29
TPZYCSandlerDimaggio::LInitialGuess
void LInitialGuess(const T &X, T &L) const
Definition: TPZYCSandlerDimaggio.h:1280
TPZYCSandlerDimaggio::LoadState
void LoadState(TPZFMatrix< REAL > &state)
Definition: TPZYCSandlerDimaggio.h:1347
TPZYCSandlerDimaggio::SetForceYield
void SetForceYield(const int forceYield)
Definition: TPZYCSandlerDimaggio.h:105
atan2
TPZFlopCounter atan2(const TPZFlopCounter &val1, const TPZFlopCounter &val2)
Returns the arc tangent in radians and increments the counter of the Arc Tangent. ATAN2 returns the ...
Definition: pzreal.h:544
CheckConvergence
void CheckConvergence(TConv &obj, TPZFMatrix< STATE > &state, TPZFMatrix< STATE > &range, TPZVec< REAL > &coefs)
Implements a general procedure to check whether the class TConv implements a consistente tangent matr...
Definition: checkconv.h:23
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZTensor::CopyTo
void CopyTo(TPZTensor< T1 > &target) const
Definition: TPZTensor.h:819
TPZTensor::Add
void Add(const TPZTensor< T1 > &tensor, const T2 &constant)
Definition: TPZTensor.h:757
TPZTensor::Norm
T Norm() const
Definition: TPZTensor.h:854
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
TPZYCSandlerDimaggio::ComputeX
void ComputeX(const T &A, T &X) const
Definition: TPZYCSandlerDimaggio.h:1316
TPZTensor::fData
TPZManVector< T, 6 > fData
Definition: TPZTensor.h:652
TPZTensor::I1
T I1() const
Definition: TPZTensor.h:903
TPZYCSandlerDimaggio::FuncEpspUsingL
REAL FuncEpspUsingL(const TPZElasticResponse &ER, REAL theta, REAL epspini, REAL L, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:332
TPZYCSandlerDimaggio::DEpspDL
void DEpspDL(const T &L, T &depspdL) const
Definition: TPZYCSandlerDimaggio.h:1714
TPZYCSandlerDimaggio::AlphaMultiplier
void AlphaMultiplier(const T &A, T &multiplier) const
Definition: TPZYCSandlerDimaggio.h:164
TPZYCSandlerDimaggio::SetUp
void SetUp(REAL A, REAL B, REAL C, REAL D, REAL R, REAL W)
Definition: TPZYCSandlerDimaggio.h:150
TPZYCSandlerDimaggio::FuncThetaL
REAL FuncThetaL(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec< REAL > &sigtrialIJ) const
Definition: TPZYCSandlerDimaggio.h:432
TPZYCSandlerDimaggio::DFuncEpsp
void DFuncEpsp(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &result) const
Definition: TPZYCSandlerDimaggio.h:365
TPZPlasticCriterion.h
as_integer
std::underlying_type< Enumeration >::type as_integer(const Enumeration value)
Definition: pzreal.h:37
TPZYCSandlerDimaggio::DFuncEpspUsingL
void DFuncEpspUsingL(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec< REAL > &sigtrialIJ, TPZVec< REAL > &result) const
Definition: TPZYCSandlerDimaggio.h:348
TPZYCSandlerDimaggio::FuncEpspL
REAL FuncEpspL(const TPZElasticResponse &ER, REAL theta, REAL L, REAL deltaL, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:384
TPZTensor::YY
T & YY() const
Definition: TPZTensor.h:578
TPZYCSandlerDimaggio::UpdateSigtrialIJ
void UpdateSigtrialIJ(const TPZElasticResponse &ER, REAL epsp, REAL theta, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:596
TPZYCSandlerDimaggio::ComputeTangent
void ComputeTangent(TPZFMatrix< REAL > &tangent, TPZVec< REAL > &, int icase)
Definition: TPZYCSandlerDimaggio.h:1363
TPZYCSandlerDimaggio::NewtonF3
void NewtonF3(const TPZElasticResponse &ER, REAL &epsp, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:725
TPZElasticResponse.h
h
clarg::argBool h("-h", "help message", false)
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
TPZYCSandlerDimaggio::NewtonF2
void NewtonF2(const TPZElasticResponse &ER, REAL &epsp, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:615
LOGPZ_WARN
#define LOGPZ_WARN(A, B)
Define log for warnings.
Definition: pzlog.h:91
TPZYCSandlerDimaggio::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: TPZYCSandlerDimaggio.cpp:4
TPZYCSandlerDimaggio::UpdateSigtrialIJL
void UpdateSigtrialIJL(const TPZElasticResponse &ER, REAL L, REAL theta, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:604
TPZYCSandlerDimaggio::McCormicRanchSand
static void McCormicRanchSand(TPZYCSandlerDimaggio &material)
Definition: TPZYCSandlerDimaggio.h:1593
TPZYCSandlerDimaggio::TPZYCSandlerDimaggio
TPZYCSandlerDimaggio(const TPZYCSandlerDimaggio &source)
Definition: TPZYCSandlerDimaggio.h:43
TPZYCSandlerDimaggio::ComputeEpsp
REAL ComputeEpsp(const REAL L) const
Definition: TPZYCSandlerDimaggio.h:303
TPZYCSandlerDimaggio::EpspFromL
void EpspFromL(const T &L, T &epsp) const
Definition: TPZYCSandlerDimaggio.h:1703
TPZYCSandlerDimaggio::ComputeF
void ComputeF(const T &L, T &F) const
Definition: TPZYCSandlerDimaggio.h:1302
TPZYCSandlerDimaggio::Read
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: TPZYCSandlerDimaggio.h:73
TPZYCSandlerDimaggio::D2EpspDL2
void D2EpspDL2(const T &L, T &d2epspdL2) const
Definition: TPZYCSandlerDimaggio.h:1726
error
void error(char *string)
Definition: testShape.cc:7
TPZYCSandlerDimaggio::DFuncTheta
void DFuncTheta(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &sigtrialIJ, TPZVec< REAL > &result) const
Definition: TPZYCSandlerDimaggio.h:518
TPZYCSandlerDimaggio::Name
const char * Name() const
Definition: TPZYCSandlerDimaggio.h:86
TPZYCSandlerDimaggio::Distance
REAL Distance(const TPZElasticResponse &ER, REAL L, TPZVec< REAL > &sigtrialIJ) const
compute the distance between the trial point and the point on the F1 curve
Definition: TPZYCSandlerDimaggio.h:1817
TPZTensor::YZ
T & YZ() const
Definition: TPZTensor.h:582
sin
sin
Definition: fadfunc.h:63
LOGPZ_INFO
#define LOGPZ_INFO(A, B)
Define log for informations.
Definition: pzlog.h:89
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZElasticResponse::ComputeStrain
void ComputeStrain(const TPZTensor< T > &sigma, TPZTensor< T > &epsilon) const
Definition: TPZElasticResponse.h:112
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
_XX_
#define _XX_
Definition: TPZTensor.h:27
TPZYCSandlerDimaggio::SetYieldStatusMode
void SetYieldStatusMode(const TPZTensor< REAL > &sigma, const REAL &A)
Definition: TPZYCSandlerDimaggio.h:116
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TFad
Definition: tfad.h:64
ELU
Definition: pzmatrix.h:52
_XY_
#define _XY_
Definition: TPZTensor.h:28
TPZYCSandlerDimaggio::operator=
TPZYCSandlerDimaggio & operator=(const TPZYCSandlerDimaggio &source)
Definition: TPZYCSandlerDimaggio.h:53
TPZYCSandlerDimaggio::InitialGuess
void InitialGuess(const TPZElasticResponse &ER, REAL epsp, TPZTensor< REAL > &sigtrial, REAL &epspproj, TPZVec< REAL > &delgamma, TPZTensor< REAL > &sigproj)
Definition: TPZYCSandlerDimaggio.h:1627
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZYCSandlerDimaggio::DFuncTheta2L
void DFuncTheta2L(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec< REAL > &sigtrialIJ, TPZVec< REAL > &result) const
Definition: TPZYCSandlerDimaggio.h:575
test.res
string res
Definition: test.py:151
TPZTensor::J2
T J2() const
Definition: TPZTensor.h:927
TPZTensor::Adjust
void Adjust(TPZVec< T > &sigIJ, TPZTensor< T > &result) const
adjust the tensor to the given values of I1 and sqj2
Definition: TPZTensor.h:490
TPZExtractVal::val
static REAL val(const int number)
Definition: pzextractval.h:21
TPZYCSandlerDimaggio::DFuncTheta2
void DFuncTheta2(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &sigtrialIJ, TPZVec< REAL > &result) const
Definition: TPZYCSandlerDimaggio.h:547
_YY_
#define _YY_
Definition: TPZTensor.h:30
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZYCSandlerDimaggio::YieldFunction
void YieldFunction(const TPZVec< STATE > &sigma, STATE kprev, TPZVec< STATE > &yield) const override
Definition: TPZYCSandlerDimaggio.h:217
TPZYCSandlerDimaggio::ComputeDL
void ComputeDL(const T &L, const T &A, T &DL) const
Definition: TPZYCSandlerDimaggio.h:1292
TPZYCSandlerDimaggio::NewtonF2L
void NewtonF2L(const TPZElasticResponse &ER, REAL &L, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:832
TPZYCSandlerDimaggio::FuncTheta2
REAL FuncTheta2(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &sigtrialIJ) const
Definition: TPZYCSandlerDimaggio.h:482
_ZZ_
#define _ZZ_
Definition: TPZTensor.h:32
TPZYCSandlerDimaggio::FuncEpsp
REAL FuncEpsp(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &sigtrialIJ)
Definition: TPZYCSandlerDimaggio.h:314
TPZYCSandlerDimaggio::Compute
void Compute(const TPZTensor< T > &sigma, const T &A, TPZVec< T > &res, int checkForcedYield) const
Calculo do criterio de plastificacao.
Definition: TPZYCSandlerDimaggio.h:1008
TPZYCSandlerDimaggio::Print
void Print(std::ostream &out) const override
Definition: TPZYCSandlerDimaggio.h:90
FL
#define FL(A)
Function for dynamic cast of the material based on map A (second data)
Definition: pzflowcmesh.cpp:147
dist
REAL dist(TPZVec< T1 > &vec1, TPZVec< T1 > &vec2)
Definition: pzvec_extras.h:124
log
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ log
Definition: tfadfunc.h:130
TPZYCSandlerDimaggio::gRefTension
static TPZTensor< REAL > gRefTension
Definition: TPZYCSandlerDimaggio.h:983
TPZYCSandlerDimaggio::FuncTheta
REAL FuncTheta(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &sigtrialIJ) const
Definition: TPZYCSandlerDimaggio.h:448
TPZYCSandlerDimaggio::ComputedF
void ComputedF(const T &L, T &dF) const
Definition: TPZYCSandlerDimaggio.h:1307
TPZTensor::XX
T & XX() const
Definition: TPZTensor.h:566
TPZYCSandlerDimaggio::ComputeD2F
void ComputeD2F(REAL L, REAL &D2F) const
Definition: TPZYCSandlerDimaggio.h:1311
TPZYCSandlerDimaggio::DFuncEpspL
void DFuncEpspL(const TPZElasticResponse &ER, REAL theta, REAL L, REAL deltaL, TPZVec< REAL > &sigtrialIJ, TPZVec< REAL > &result) const
Definition: TPZYCSandlerDimaggio.h:414
TPZYCSandlerDimaggio::N
void N(const TPZTensor< T > &sigma, const T &A, TPZVec< TPZTensor< T > > &Ndir, int checkForcedYield) const
Definition: TPZYCSandlerDimaggio.h:1073
TPZYCSandlerDimaggio::D2Distance
REAL D2Distance(const TPZElasticResponse &ER, REAL L, TPZVec< REAL > &sigtrialIJ) const
compute the second derivative of the distance with respect to L
Definition: TPZYCSandlerDimaggio.h:1850
TPZYCSandlerDimaggio::DDistance
REAL DDistance(const TPZElasticResponse &ER, REAL L, TPZVec< REAL > &sigtrialIJ) const
compute the derivative of the distance with respect to L
Definition: TPZYCSandlerDimaggio.h:1832
TPZTensor::XY
T & XY() const
Definition: TPZTensor.h:570
EMathematicaInput
Definition: pzmatrix.h:55
asin
expr_ expr_ expr_ expr_ expr_ expr_ asin
Definition: tfadfunc.h:80
TPZVec::Fill
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
TPZTensor::ZZ
T & ZZ() const
Definition: TPZTensor.h:586
exp
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ exp
Definition: tfadfunc.h:125
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZYCSandlerDimaggio::Write
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: TPZYCSandlerDimaggio.h:64
TPZYCSandlerDimaggio::DistTheta
REAL DistTheta(const TPZElasticResponse &ER, REAL theta, REAL epsp, REAL delepsp, TPZVec< REAL > &sigtrialIJ) const
Definition: TPZYCSandlerDimaggio.h:508
TPZYCSandlerDimaggio::NewtonF1
void NewtonF1(const TPZElasticResponse &ER, REAL &L, TPZVec< REAL > &sigtrialIJ)
Projeto o ponto sobre a superficie F1, atualiza o L.
Definition: TPZYCSandlerDimaggio.h:1739
TPZTensor.h
TPZTensor::XZ
T & XZ() const
Definition: TPZTensor.h:574
checkconv.h
Contains the implementation of the CheckConvergence function.
cos
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TPZYCSandlerDimaggio::DistThetaL
REAL DistThetaL(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec< REAL > &sigtrialIJ) const
Definition: TPZYCSandlerDimaggio.h:492
TPZYCSandlerDimaggio::SolveL
void SolveL(const T &X, T &L, REAL relTol=1.e-10) const
Definition: TPZYCSandlerDimaggio.h:1252
TPZYCSandlerDimaggio::Residual
void Residual(TPZFMatrix< REAL > &res, int icase)
Definition: TPZYCSandlerDimaggio.h:1406
TPZYCSandlerDimaggio::FuncTheta2L
REAL FuncTheta2L(const TPZElasticResponse &ER, REAL theta, REAL L, TPZVec< REAL > &sigtrialIJ) const
Definition: TPZYCSandlerDimaggio.h:459
TPZTensor::S
void S(TPZTensor< T > &s) const
Definition: TPZTensor.h:894
TPZMatrix::SolveDirect
virtual int SolveDirect(TPZFMatrix< TVar > &F, const DecomposeType dt, std::list< int64_t > &singular)
Solves the linear system using Direct methods.
Definition: pzmatrix.cpp:710
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
TPZYCSandlerDimaggio::H
void H(const TPZTensor< T > &sigma, const T &A, TPZVec< T > &h, int checkForcedYield) const
Definition: TPZYCSandlerDimaggio.h:1202
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13