NeoPZ
TPZMohrCoulombNeto.h
Go to the documentation of this file.
1 /*
2  * TPZMohrCoulomb.h
3  * FEMPZ
4  *
5  * Created by Diogo Cecilio on 5/4/10.
6  * Copyright 2010 __MyCompanyName__. All rights reserved.
7  *
8  */
9 
10 
11 /* Generated by Together */// $Id: TPZMohrCoulomb.h,v 1.2 2010-06-11 22:12:14 diogo Exp $
12 
13 #ifndef TPZMOHRCOULOMBNETO_H
14 #define TPZMOHRCOULOMBNETO_H
15 
16 #include "pzlog.h"
17 #include "TPZTensor.h"
18 #include "pzvec_extras.h"
19 #include "TPZPlasticState.h"
20 
21 #ifdef LOG4CXX
22 static LoggerPtr loggerMohrCoulomb(Logger::getLogger("pz.plasticity.mohrcoulombneto"));
23 #endif
24 
25 
26 
28 {
29 
30  REAL fYoung;
31  REAL fPoisson;
32  REAL fPhi;
33  REAL fPsi;
34  REAL coesion;
35 
36 public:
37 
40  {
42  {
43 
44  }
45 
47  {
48 
49  }
50 
52  {
53  fEpsPlastic = copy.fEpsPlastic;
55  return *this;
56  }
57 
58  void Print(std::ostream &out) const
59  {
60  out << "Plastic Deformation tensor ";
61  fEpsPlastic.Print(out);
62  out << "Acumulated plastic deformation " << fEpsPlasticBar << std::endl;
63  }
64 
67 
70  };
71 
72 
74  // we can only expect a consistent tangent matrix if the decision tree remains the same
76  {
77  TComputeSequence() : fWhichPlane(ENoPlane), fGamma(0)
78  {
79 
80  }
81 
82  TComputeSequence(const TComputeSequence &copy) : fWhichPlane(copy.fWhichPlane), fGamma(copy.fGamma)
83  {
84 
85  }
86 
88  {
89  fWhichPlane = copy.fWhichPlane;
90  fGamma = copy.fGamma;
91  return *this;
92  }
93 
94  enum MPlane {ENoPlane, EElastic, EMainPlane, ERightEdge, ELeftEdge, EHydroStatic };
95 
97 
99 
100 
101  };
102 
103 protected:
104 
107 
108 
109 public:
110 
111  TPZMohrCoulombNeto() : fYoung(20000.), fPoisson(0.), fPhi(M_PI/9.), fPsi(M_PI/9.)
112  {
113  this->coesion = 9.35;
114  }
115  REAL Lambda()
116  {
117  return fPoisson*fYoung/(1.+fPoisson)*(1.-2.*fPoisson);
118  }
119  REAL Mu()
120  {
121  return fYoung/(2.*(1.+fPoisson));
122  }
123  REAL G()
124  {
125  return fYoung/(2.*(1+fPoisson));
126  }
127  REAL K()
128  {
129  return fYoung/(3.*(1.-2.*fPoisson));
130  }
131 
132  void Print(std::ostream &out) const
133  {
134  out << "TPZMohrCoulombNeto\n";
135  fState.Print(out);
136  }
137 
139  template<class T>
140  void PlasticityFunction(T epsp, T &sigmay, T &H) const
141  {
142  //sigmay = T(15.)+(T(2571.43)-T(2.95238e6)*epsp)*(T(-0.0035)+epsp);
143  // H = T(12904.8)-T(5.90476e6)*epsp;
144  //sigmay=20.;
145  // H=0.;
146  sigmay = T(15.)+(T(2571.43))*(T(-0.0035)) + T(12904.8)*epsp;
147  H = T(12904.8);
148  }
149 
151  template<class T>
152  void PieceWise(T epbar, T &m, T & fx)const
153  {
154 
155  ifstream file("curvadehardening.txt");
156  TPZFMatrix<REAL> mat;
157  int sz;
158  file >>sz;
159  mat.Resize(sz,2);
160 
161  for(int i=0;i<sz-1;i++)
162  {
163 
164  file >> mat(i,0);
165  file >> mat(i,1);
166  }
167 
168  //cout << "\n mat = "<< mat <<endl;
169  for(int i=0;i<sz-1;i++)
170  {
171  if(epbar >= mat(i,0) && epbar <= mat(i+1,0))
172  {
173  REAL x0 = mat(i,0);
174  REAL y0 = mat(i,1);
175  REAL x = mat(i+1,0);
176  REAL y = mat(i+1,1);
177  m = (y - y0)/(x - x0);
178  fx=m*(epbar-x0)+y0;
179  // return;
180  }
181 
182  }
183  }
184 
185 
186  template<class T>
188  {
189  T trdeform = deform.I1();
190  TPZTensor<T> result;
191  result.Identity();
192  result *= (Lambda()*trdeform);
193  result.Add(deform,2.*Mu());
194  return result;
195  }
196 
197  template<class T>
199  {
200  TPZTensor<T> epslocal(epstotal);
201  for(int i=0; i<6; i++)
202  {
203  epslocal[i] -= T(fState.fEpsPlastic[i]);
204  }
205  TPZTensor<T> sigma;
206  sigma = SigmaElast(epslocal);
207  typename TPZTensor<T>::TPZDecomposed sigma_trial;
208  sigma.EigenSystem(sigma_trial);
209 
210 #ifdef LOG4CXX
211  if (loggerMohrCoulomb->isDebugEnabled()) {
212  std::stringstream sout;
213  sout << "Input stress tensor ";
214  sigma.Print(sout);
215  sout << "Input tensor in decomposed form\n";
216  sigma_trial.Print(sout);
217  LOGPZ_DEBUG(loggerMohrCoulomb, sout.str())
218  }
219 #endif
220  return sigma_trial;
221  }
222 
224  {
225  typedef TFad<6,REAL> fadtype;
226  TPZTensor<fadtype> epstotalFAD, sigmaElastFAD, sigmaFAD;
227  for (int i=0; i<6; i++) {
228  epstotalFAD[i].val() = epstotal[i];
229  epstotalFAD[i].fastAccessDx(i) = 1.;
230  }
231  sigmaElastFAD = SigmaElast(epstotalFAD);
232 
233  switch (memory.fWhichPlane) {
235  DebugStop();
236  break;
238  sigmaFAD = sigmaElastFAD;
239  break;
243  {
244  TPZTensor<fadtype>::TPZDecomposed sigma_trial = SigmaTrial(epstotalFAD);
245  TPZTensor<fadtype>::TPZDecomposed sigma_projected;
246  TComputeSequence locmem(memory);
247  switch (memory.fWhichPlane) {
249  ReturnMapPlane<fadtype>(sigma_trial, sigma_projected, locmem);
250  break;
252  ReturnMapLeftEdge<fadtype>(sigma_trial, sigma_projected, locmem);
253  break;
255  ReturnMapRightEdge<fadtype>(sigma_trial, sigma_projected, locmem);
256  break;
257  default:
258  DebugStop();
259  break;
260  }
261  sigmaFAD = TPZTensor<fadtype>(sigma_projected);
262  }
263  default:
264  break;
265  }
266  for (int i=0; i<6; i++) {
267  sigma[i] = sigmaFAD[i].val();
268  for (int j=0; j<6; j++) {
269  tangent(i,j) = sigmaFAD[i].fastAccessDx(j);
270  }
271  }
272  }
273 
275  {
276  const REAL cosphi = cos(fPhi);
277  TPZTensor<REAL> sigma;
278  switch (memory.fWhichPlane) {
280  {
281  TPZTensor<REAL> epslocal(epstotal);
282  epslocal -= fState.fEpsPlastic;
283  sigma = SigmaElast(epslocal);
284  //sigma = SigmaElast(epstotal);
285  }
286  break;
290  {
291  TPZTensor<REAL>::TPZDecomposed sigma_trial = SigmaTrial(epstotal);
292  TPZTensor<REAL>::TPZDecomposed sigma_projected;
293  TComputeSequence locmem(memory);
294  switch (memory.fWhichPlane) {
296  ReturnMapPlane<REAL>(sigma_trial, sigma_projected, locmem);
297  fState.fEpsPlasticBar+=(locmem.fGamma[0]*2.*cosphi);
298  break;
300  ReturnMapLeftEdge<REAL>(sigma_trial, sigma_projected, locmem);
301  fState.fEpsPlasticBar+=(locmem.fGamma[0]+locmem.fGamma[1])*2.*cosphi;
302  break;
304  ReturnMapRightEdge<REAL>(sigma_trial, sigma_projected, locmem);
305  fState.fEpsPlasticBar+=(locmem.fGamma[0]+locmem.fGamma[1])*2.*cosphi;
306  break;
307  default:
308  DebugStop();
309  break;
310  }
311  sigma = TPZTensor<REAL>(sigma_projected);
312  break;
313  }
314  default:
315  {
316  break;
317  }
318  }
319 
320 
321  REAL tempval;
322  TPZTensor<REAL> StressDeviatoric,P,I,epsplastic(epstotal),epselastic;
323 
324  P.XX()=1; I.XX()=1;
325  P.YY()=1; I.YY()=1;
326  P.ZZ()=1; I.ZZ()=1;
327 
328 
329  tempval=(sigma.I1()/3.)*(1./(3.*K()));
330  P.Multiply(tempval,1);
331 
332 
333  cout << " \n P = "<< P <<endl;
334  sigma.S(StressDeviatoric);
335 
336  StressDeviatoric*=1./(2.*G());
337  cout << " \n S = "<< StressDeviatoric <<endl;
338 
339  StressDeviatoric.Add(P,1);
340  epselastic=StressDeviatoric;
341 
342  epsplastic-=epselastic;
343 
344 
345  fState.fEpsPlastic = epsplastic;
346 
347 
348 
349  }
350 
351  template<class T>
353  {
354  typename TPZTensor<T>::TPZDecomposed sigma_trial = SigmaTrial(epstotal);
355  TComputeSequence memory;
356  T phi = PhiPlane<T>(sigma_trial);
357  if (TPZExtractVal::val(phi) <= 0.) {
359  memory.fGamma.Resize(0);
360  sigma = TPZTensor<T>(sigma_trial);
361  //state.m_eps_t = epstotal;
362  return memory;
363  }
364  typename TPZTensor<T>::TPZDecomposed sigma_projected;
365  memory.fGamma.Resize(1);
366  memory.fGamma[0] = 0.;
367  if (ReturnMapPlane<T>(sigma_trial, sigma_projected, memory)) {
368  sigma = TPZTensor<T>(sigma_projected);
370  }
371  else {
372  memory.fGamma.Resize(2);
373  memory.fGamma[0] = 0.;
374  memory.fGamma[1] = 0.;
375 
376 
377  const REAL sinpsi = sin(fPsi);
378  TPZManVector<T,3> &eigenvalues = sigma_trial.fEigenvalues;
379  REAL val = (1-sinpsi)*TPZExtractVal::val(eigenvalues[0])-2.*TPZExtractVal::val(eigenvalues[2])+(1+sinpsi)*TPZExtractVal::val(eigenvalues[1]);
380  if (val > 0.) {
381  ReturnMapRightEdge<T>(sigma_trial, sigma_projected, memory);
383  }
384  else {
385  ReturnMapLeftEdge<T>(sigma_trial, sigma_projected, memory);
387  }
388 #ifdef LOG4CXX
389  {
390  std::stringstream sout;
391  sout << "After the map to the edge, sigma_projected :\n";
392  sigma_projected.Print(sout);
393  LOGPZ_DEBUG(loggerMohrCoulomb, sout.str())
394  }
395 #endif
396  sigma = TPZTensor<T>(sigma_projected);
397  }
398  CommitDeformation(epstotal,memory);
399  return memory;
400  }
401 
402  template<class T>
403  T PhiPlane(typename TPZTensor<T>::TPZDecomposed &sigma) const
404  {
405  const REAL sinphi = sin(fPhi);
406  const REAL cosphi = cos(fPhi);
407  T sigmay,H;
408  PlasticityFunction(T(fState.fEpsPlasticBar),sigmay, H);
409  return sigma.fEigenvalues[0]-sigma.fEigenvalues[2]+(sigma.fEigenvalues[0]+sigma.fEigenvalues[2])*sinphi-2.*sigmay*cosphi;
410  }
411 
412  template<class T>
413  bool ReturnMapPlane(const typename TPZTensor<T>::TPZDecomposed &sigma_trial, typename TPZTensor<T>::TPZDecomposed &sigma_projected,
414  TComputeSequence &memory)
415  {
416  sigma_projected = sigma_trial;
417  TPZManVector<T,3> &eigenvalues = sigma_projected.fEigenvalues;
418 // TPZManVector<TPZTensor<T>,3> &eigenvectors = sigma_projected.fEigenvectors;
419  const REAL sinphi = sin(fPhi);
420  const REAL sinpsi = sin(fPsi);
421  const REAL cosphi = cos(fPhi);
422  const REAL sinphi2 = sinphi*sinphi;
423  const REAL cosphi2 = 1.-sinphi2;
424  const REAL constA = 4.* G() *(1.+ sinphi*sinpsi/3.) + 4.*K() * sinphi*sinpsi;
425  T sigmay,H;
426  T epsbar = T(fState.fEpsPlasticBar+memory.fGamma[0]*2.*cosphi);//diogo aqui
427  PlasticityFunction(epsbar,sigmay, H);
428  T phi = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*sinphi-2.*sigmay*cosphi;
429  T gamma = memory.fGamma[0];
430  REAL phival = TPZExtractVal::val(phi);
431  REAL tolerance = 1.e-8;
432  do {
433  T denom = -constA- T(4.*cosphi2)*H;
434 // T d = T(-4.*G()*(1.+sinphi*sinpsi/3.)-4.*K()*sinphi*sinpsi)-T(4.*cosphi2)*H;
435  T deriv_gamma = -phi/denom;
436  gamma += deriv_gamma;
437  epsbar = T(fState.fEpsPlasticBar)+gamma*T(2.*cosphi);
438  PlasticityFunction(epsbar, sigmay, H);
439  if (TPZExtractVal::val(H) < 0.) {
440  DebugStop();
441  }
442  phi = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*sinphi-2.*sigmay*cosphi-constA*gamma;
443  phival = TPZExtractVal::val(phi);
444 
445  } while (abs(phival) > tolerance);
446 
447 
448  memory.fGamma[0] = TPZExtractVal::val(gamma);
449  eigenvalues[0] -= T(2.*G()*(1+sinpsi/3.)+2.*K()*sinpsi)*gamma;
450  eigenvalues[1] += T((4.*G()/3. - K()*2.)*sinpsi)*gamma;
451  eigenvalues[2] += T(2.*G()*(1-sinpsi/3.)-2.*K()*sinpsi)*gamma;
452 #ifdef PZDEBUG
453  phi = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*sinphi-2.*sigmay*cosphi;
454 #endif
455  return (TPZExtractVal::val(eigenvalues[0])>TPZExtractVal::val(eigenvalues[1]) && TPZExtractVal::val(eigenvalues[1]) > TPZExtractVal::val(eigenvalues[2]));
456  }
457 
458  template<class T>
459  bool ReturnMapLeftEdge(const typename TPZTensor<T>::TPZDecomposed &sigma_trial, typename TPZTensor<T>::TPZDecomposed &sigma_projected,
460  TComputeSequence &memory)
461  {
462 
463  sigma_projected = sigma_trial;
464  TPZManVector<T,3> &eigenvalues = sigma_projected.fEigenvalues;
465 // TPZManVector<TPZTensor<T>,3> &eigenvectors = sigma_projected.fEigenvectors;
466  const REAL sinphi = sin(fPhi);
467  const REAL sinpsi = sin(fPsi);
468  const REAL cosphi = cos(fPhi);
469  const REAL sinphi2 = sinphi*sinphi;
470  const REAL cosphi2 = 1.-sinphi2;
471  TPZManVector<T,2> gamma(2,0.),phi(2,0.),sigma_bar(2,0.),ab(2,0.);
472  gamma[0] = memory.fGamma[0];
473  gamma[1] = memory.fGamma[1];
474  TPZManVector<REAL,2> phival(2,0.);
475  TPZFNMatrix<4,T> d(2,2,0.), dinverse(2,2,0.);
476  sigma_bar[0] = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*T(sinphi);
477  sigma_bar[1] = eigenvalues[1]-eigenvalues[2]+(eigenvalues[1]+eigenvalues[2])*T(sinphi);
478  T sigmay,H;
479  T epsbar = T(fState.fEpsPlasticBar) + (gamma[0]+gamma[1])*T(2.*cosphi);//diogo aqui
480  PlasticityFunction(epsbar,sigmay, H);
481  phi[0] = sigma_bar[0] - T(2.*cosphi)*sigmay;
482  phi[1] = sigma_bar[1] - T(2.*cosphi)*sigmay;
483  ab[0] = T(4.*G()*(1+sinphi*sinpsi/3.)+4.*K()*sinphi*sinpsi);
484  ab[1] = T(2.*G()*(1.-sinphi-sinpsi-sinphi*sinpsi/3.)+4.*K()*sinphi*sinpsi);
485  T residual =1;
486  REAL tolerance = 1.e-8;
487  do {
488  d(0,0) = -ab[0]-T(4.*cosphi2)*H;
489  d(1,0) = -ab[1]-T(4.*cosphi2)*H;
490  d(0,1) = -ab[1]-T(4.*cosphi2)*H;
491  d(1,1) = -ab[0]-T(4.*cosphi2)*H;
492  T detd = d(0,0)*d(1,1)-d(0,1)*d(1,0);
493  dinverse(0,0) = d(1,1)/detd;
494  dinverse(1,0) = -d(1,0)/detd;
495  dinverse(0,1) = -d(0,1)/detd;
496  dinverse(1,1) = d(0,0)/detd;
497  gamma[0] -= (dinverse(0,0)*phi[0]+dinverse(0,1)*phi[1]);
498  gamma[1] -= (dinverse(1,0)*phi[0]+dinverse(1,1)*phi[1]);
499  //T epsbar = T(fState.fEpsPlasticBar)+(gamma[0]+gamma[1])*T(2.*cosphi); diogo aqui
500  epsbar = T(fState.fEpsPlasticBar)+(gamma[0]+gamma[1])*T(2.*cosphi);
501  PlasticityFunction(epsbar, sigmay, H);
502  phi[0] = sigma_bar[0] - ab[0]*gamma[0] - ab[1]*gamma[1] - T(2.*cosphi)*sigmay;
503  phi[1] = sigma_bar[1] - ab[1]*gamma[0] - ab[0]*gamma[0] - T(2.*cosphi)*sigmay;
504  phival[0] = TPZExtractVal::val(phi[0]);
505  phival[1] = TPZExtractVal::val(phi[1]);
506  residual=(fabs(phival[0])+fabs(phival[1]))/sigmay;//aqui diogo
507  //} while (abs(phival[0]) > tolerance || abs(phival[1]) > tolerance);//aqui diogo
508  }while (residual>tolerance);//aqui diogo
509 // eigenvalues[0] -= T(2.*G()*(1+sinpsi/3.)+2.*K()*sinpsi)*gamma;
510 // eigenvalues[1] += T((4.*G()/3. - K()*2.)*sinpsi)*gamma;
511 // eigenvalues[2] += T(2.*G()*(1-sinpsi/3.)-2.*K()*sinpsi)*gamma;
512 
513  memory.fGamma[0] = TPZExtractVal::val(gamma[0]);
514  memory.fGamma[1] = TPZExtractVal::val(gamma[1]);
515  eigenvalues[0] -= T(2.*G()*(1+sinpsi/3.)+2.*K()*sinpsi)*gamma[0]+T((4.*G()/3.-2.*K())*sinpsi)*gamma[1];
516  eigenvalues[1] += T((4.*G()/3.- K()*2.)*sinpsi)*gamma[0]-T(2.*G()*(1.+sinpsi/3.)+2.*K()*sinpsi)*gamma[1];
517  eigenvalues[2] -= T(2.*G()*(1-sinpsi/3.)-2.*K()*sinpsi)*(gamma[0]+gamma[1]);
518  return (TPZExtractVal::val(eigenvalues[0])>TPZExtractVal::val(eigenvalues[1]) && TPZExtractVal::val(eigenvalues[1]) > TPZExtractVal::val(eigenvalues[2]));
519  }
520 
521  template<class T>
522  bool ReturnMapRightEdge(const typename TPZTensor<T>::TPZDecomposed &sigma_trial, typename TPZTensor<T>::TPZDecomposed &sigma_projected,
523  TComputeSequence &memory)
524  {
525  sigma_projected = sigma_trial;
526  TPZManVector<T,3> &eigenvalues = sigma_projected.fEigenvalues;
527 // TPZManVector<TPZTensor<T>,3> &eigenvectors = sigma_projected.fEigenvectors;
528  const REAL sinphi = sin(fPhi);
529  const REAL sinpsi = sin(fPsi);
530  const REAL cosphi = cos(fPhi);
531  const REAL sinphi2 = sinphi*sinphi;
532  const REAL cosphi2 = 1.-sinphi2;
533  const REAL KV = K();
534  const REAL GV = G();
535  TPZManVector<T,2> gamma(2,0.),phi(2,0.),sigma_bar(2,0.),ab(2,0.);
536  gamma[0] = memory.fGamma[0];
537  gamma[1] = memory.fGamma[1];
538  TPZManVector<REAL,2> phival(2,0.);
539  TPZFNMatrix<4,T> d(2,2,0.), dinverse(2,2,0.);
540  sigma_bar[0] = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*T(sinphi);
541  sigma_bar[1] = eigenvalues[0]-eigenvalues[1]+(eigenvalues[0]+eigenvalues[1])*T(sinphi);
542  T sigmay,H;
543  T epsbar = T(fState.fEpsPlasticBar)+(gamma[0]+gamma[1])*T(2.*cosphi);
544  PlasticityFunction(epsbar,sigmay, H);
545  phi[0] = sigma_bar[0] - T(2.*cosphi)*sigmay;
546  phi[1] = sigma_bar[1] - T(2.*cosphi)*sigmay;
547  ab[0] = T(4.*GV*(1+sinphi*sinpsi/3.)+4.*KV*sinphi*sinpsi);
548  ab[1] = T(2.*GV*(1.+sinphi+sinpsi-sinphi*sinpsi/3.)+4.*KV*sinphi*sinpsi);
549 #ifdef LOG4CXX
550  {
551  std::stringstream sout;
552  sout << "phi = " << phi << std::endl;
553  LOGPZ_DEBUG(loggerMohrCoulomb, sout.str())
554  }
555 #endif
556 
557 //#ifdef PZDEBUG
558 // gamma[0] = 0.;
559 // gamma[1] = 1.;
560 // T v[3];
561 // v[0] = -T(2.*GV*(1+sinpsi/3.)+2.*KV*sinpsi)*(gamma[0]+gamma[1]);
562 // v[1] = T((4.*GV/3.- KV*2.)*sinpsi)*gamma[0]+T(2.*GV*(1.-sinpsi/3.)-2.*KV*sinpsi)*gamma[1];
563 // v[2] = T(2.*GV*(1-sinpsi/3.)-2.*KV*sinpsi)*gamma[0]+T((4.*GV/3.-2.*KV)*sinpsi)*gamma[1];
564 // eigenvalues[0] = sigma_trial.fEigenvalues[0]+v[0];
565 // eigenvalues[1] = sigma_trial.fEigenvalues[1]+v[1];
566 // eigenvalues[2] = sigma_trial.fEigenvalues[2]+v[2];
567 //
568 // T test1 = (v[0]-v[2])+(v[0]+v[2])*sinphi;
569 // T test2 = (v[0]-v[1])+(v[0]+v[1])*sinphi;
570 //
571 // sigma_bar[0] = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*T(sinphi);
572 // sigma_bar[1] = eigenvalues[0]-eigenvalues[1]+(eigenvalues[0]+eigenvalues[1])*T(sinphi);
573 // T A = phi[0] - sigma_bar[0] + T(2.*cosphi)*sigmay;
574 // T B = phi[1] - sigma_bar[1] + T(2.*cosphi)*sigmay;
575 //
576 //#endif
577  REAL tolerance = 1.e-8;
578  int iter = 0;
579  T residual =1;
580  do {
581 #ifdef LOG4CXX
582  {
583  std::stringstream sout;
584  sout << "epsbar = " << epsbar << std::endl;
585  sout << "sigmay = " << sigmay << std::endl;
586  sout << "H = " << H << std::endl;
587  LOGPZ_DEBUG(loggerMohrCoulomb, sout.str())
588  }
589 #endif
590  d(0,0) = -ab[0]-T(4.*cosphi2)*H;
591  d(1,0) = -ab[1]-T(4.*cosphi2)*H;
592  d(0,1) = -ab[1]-T(4.*cosphi2)*H;
593  d(1,1) = -ab[0]-T(4.*cosphi2)*H;
594  T detd = d(0,0)*d(1,1)-d(0,1)*d(1,0);
595  dinverse(0,0) = d(1,1)/detd;
596  dinverse(1,0) = -d(1,0)/detd;
597  dinverse(0,1) = -d(0,1)/detd;
598  dinverse(1,1) = d(0,0)/detd;
599  gamma[0] -= (dinverse(0,0)*phi[0]+dinverse(0,1)*phi[1]);
600  gamma[1] -= (dinverse(1,0)*phi[0]+dinverse(1,1)*phi[1]);
601  epsbar = T(fState.fEpsPlasticBar)+(gamma[0]+gamma[1])*T(2.*cosphi);
602  PlasticityFunction(epsbar, sigmay, H);
603  if (TPZExtractVal::val(H) < 0.) {
604  DebugStop();
605  }
606  iter++;
607  phi[0] = sigma_bar[0] - ab[0]*gamma[0] - ab[1]*gamma[1] - T(2.*cosphi)*sigmay;
608  phi[1] = sigma_bar[1] - ab[1]*gamma[0] - ab[0]*gamma[1] - T(2.*cosphi)*sigmay;
609  phival[0] = TPZExtractVal::val(phi[0]);
610  phival[1] = TPZExtractVal::val(phi[1]);
611 #ifdef LOG4CXX
612  {
613  std::stringstream sout;
614  sout << "iter = " << iter << " phi = " << phival << std::endl;
615  LOGPZ_DEBUG(loggerMohrCoulomb, sout.str())
616  }
617 #endif
618  residual=(fabs(phival[0])+fabs(phival[1]))/sigmay;//aqui diogo
619  cout << "\n residula = "<< endl;
620  //} while (abs(phival[0]) > tolerance || abs(phival[1]) > tolerance);//aqui diogo
621  }while (residual>tolerance);//aqui diogo
622 
623 // eigenvalues[0] -= T(2.*GV*(1+sinpsi/3.)+2.*KV*sinpsi)*gamma;
624 // eigenvalues[1] += T((4.*GV/3. - KV*2.)*sinpsi)*gamma;
625 // eigenvalues[2] += T(2.*GV*(1-sinpsi/3.)-2.*KV*sinpsi)*gamma;
626  // epsbar = T(state.m_hardening)+(gamma[0]+gamma[1])*T(2.*cosphi);
627 
628 
629  memory.fGamma[0] = TPZExtractVal::val(gamma[0]);
630  memory.fGamma[1] = TPZExtractVal::val(gamma[1]);
631 #ifdef LOG4CXX
632  {
633  std::stringstream sout;
634  sout << "gamma = " << gamma << std::endl;
635  sout << "phival = " << phival << std::endl;
636  sout << "ab = " << ab << std::endl;
637  sout << "sigma_bar = " << sigma_bar << std::endl;
638  d.Print("Jacobian",sout);
639  dinverse.Print("Inverse Jacobian",sout);
640  sout << "epsbar = " << epsbar << std::endl;
641  LOGPZ_DEBUG(loggerMohrCoulomb, sout.str())
642  }
643 #endif
644  eigenvalues[0] -= T(2.*GV*(1+sinpsi/3.)+2.*KV*sinpsi)*(gamma[0]+gamma[1]);
645  eigenvalues[1] += T((4.*GV/3.- KV*2.)*sinpsi)*gamma[0]+T(2.*GV*(1.-sinpsi/3.)-2.*KV*sinpsi)*gamma[1];
646  eigenvalues[2] += T(2.*GV*(1-sinpsi/3.)-2.*KV*sinpsi)*gamma[0]+T((4.*GV/3.-2.*KV)*sinpsi)*gamma[1];
647 /*#ifdef PZDEBUG
648  sigma_bar[0] = eigenvalues[0]-eigenvalues[2]+(eigenvalues[0]+eigenvalues[2])*T(sinphi);
649  sigma_bar[1] = eigenvalues[0]-eigenvalues[1]+(eigenvalues[0]+eigenvalues[1])*T(sinphi);
650  //epsbar = T(state.m_hardening)+(gamma[0]+gamma[1])*T(2.*cosphi);
651  PlasticityFunction(epsbar, sigmay, H);
652 
653  phi[0] = sigma_bar[0] - T(2.*cosphi)*sigmay;
654  phi[1] = sigma_bar[1] - T(2.*cosphi)*sigmay;
655 #endif*/
656  return (TPZExtractVal::val(eigenvalues[0])>TPZExtractVal::val(eigenvalues[1]) && TPZExtractVal::val(eigenvalues[1]) > TPZExtractVal::val(eigenvalues[2]));
657  }
658 };
659 
660 
661 #endif //TPZMohrCoulomb_H
void CommitDeformation(TPZTensor< REAL > &epstotal, TComputeSequence &memory)
void Print(std::ostream &out=std::cout)
Prints the structural information of the vector object to the output stream. This method will not p...
Definition: pzvec.h:482
bool ReturnMapRightEdge(const typename TPZTensor< T >::TPZDecomposed &sigma_trial, typename TPZTensor< T >::TPZDecomposed &sigma_projected, TComputeSequence &memory)
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
void Add(const TPZTensor< T1 > &tensor, const T2 &constant)
Definition: TPZTensor.h:757
void ComputeSigmaTangent(TPZTensor< REAL > &epstotal, TPZTensor< REAL > &sigma, TPZFNMatrix< 36, REAL > &tangent, const TComputeSequence &memory)
void Multiply(const T1 &multipl, const T2 &constant)
Definition: TPZTensor.h:766
void Print(std::ostream &out) const
Definition: TPZTensor.h:71
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.
TFad< 3, REAL > fadtype
const double tolerance
Tolerance value (Is zero)
TPZTensor< T >::TPZDecomposed SigmaTrial(const TPZTensor< T > &epstotal)
T I1() const
Definition: TPZTensor.h:903
T & YY() const
Definition: TPZTensor.h:578
void PieceWise(T epbar, T &m, T &fx) const
a piecewise linear hardening function
void Print(std::ostream &out) const
TComputeSequence(const TComputeSequence &copy)
TinyFad< 8, T > abs(const TinyFad< 8, T > &in)
Definition: tinyfadeight.h:846
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
void Print(std::ostream &out) const
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
void PlasticityFunction(T epsp, T &sigmay, T &H) const
the hardening function and its derivative
bool ReturnMapLeftEdge(const typename TPZTensor< T >::TPZDecomposed &sigma_trial, typename TPZTensor< T >::TPZDecomposed &sigma_projected, TComputeSequence &memory)
bool ReturnMapPlane(const typename TPZTensor< T >::TPZDecomposed &sigma_trial, typename TPZTensor< T >::TPZDecomposed &sigma_projected, TComputeSequence &memory)
sin
Definition: fadfunc.h:63
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
Internal structure to represent the plastic memory (plastic deformation and damage) ...
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
void EigenSystem(TPZDecomposed &eigensystem) const
Definition: TPZTensor.h:1265
Definition: tfad.h:64
TComputeSequence ComputeSigma(TPZTensor< T > &epstotal, TPZTensor< T > &sigma)
REAL fEpsPlasticBar
accumulated damage
void Identity()
Definition: TPZTensor.h:705
static REAL val(const int number)
Definition: pzextractval.h:21
TPZTensor< T > SigmaElast(const TPZTensor< T > &deform)
T PhiPlane(typename TPZTensor< T >::TPZDecomposed &sigma) const
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
void Print(std::ostream &out) const
Definition: TPZTensor.h:1847
TPZMohrCoulombNeto::TPlasticState fState
information of the plastic state of the material point
TPZTensor< REAL > fEpsPlastic
plastic deformation tensor
T & XX() const
Definition: TPZTensor.h:566
structure which contains the decision tree of the return map
T & ZZ() const
Definition: TPZTensor.h:586
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZManVector< T, 3 > fEigenvalues
Definition: TPZTensor.h:48
Extra utilities for TPZVec. Implementations of the saxpy, sscal, sdot, intercept, max and min functio...
clarg::argString m("-m", "input matrix file name (text format)", "matrix.txt")
TPlasticState(const TPlasticState &copy)
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TComputeSequence & operator=(const TComputeSequence &copy)
TPlasticState & operator=(const TPlasticState &copy)
void S(TPZTensor< T > &s) const
Definition: TPZTensor.h:894
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716