NeoPZ
TPZPlasticityTest.h
Go to the documentation of this file.
1 // $Id: TPZPlasticityTest.h,v 1.32 2010-08-20 20:56:45 diogo Exp $
2 // This file implements the interfaces for the
3 // most general plasticity load tests
4 
5 #include <iostream>
6 
7 using namespace std;
8 
9 #include "pzlog.h"
10 
11 #ifdef LOG4CXX // LOG4CXX may be defined alone or with LOG4CXX_PLASTICITY. The latter shall not be used alone.
12 #include <log4cxx/logger.h>
13 #include <log4cxx/basicconfigurator.h>
14 #include <log4cxx/propertyconfigurator.h>
15 
16 static LoggerPtr plasticIntegrLogger(Logger::getLogger("plasticity.plasticIntegr"));
17 #endif
18 
19 
20 #ifdef LOG4CXX_PLASTICITY
21 static LoggerPtr testLogger(Logger::getLogger("plasticity.test"));
22 #endif
23 
24 
25 
26 #ifdef LOG4CXX_PLASTICITY
27 static LoggerPtr MaterialPoint(Logger::getLogger("MaterialPointTest"));
28 #endif
29 
30 #include "pzvec.h"
31 #include "pzfmatrix.h"
32 #include "TPZTensor.h"
33 
34 #include "TPZLadeKim.h"
35 #include "TPZSandlerDimaggio.h"
36 
37 #include "TPZYCDruckerPrager.h"
38 #include "TPZThermoForceA.h"
39 #include "TPZElasticResponse.h"
40 
41 
42 #include "TPZYCMohrCoulomb.h"
44 
45 #include "TPZYCWillamWarnke.h"
46 
47 #include "TPZYCVonMises.h"
48 #include "TPZVonMises.h"
49 #include "TPZDruckerPrager.h"
50 
51 
52 
54 {
55 public:
56 
57  static void InitializeLOG();
58 
59  template <class T>
60  static void ReciprocityTest(T & plasticModel, TPZTensor<REAL> strain1);
61 
62  template <class T>
63  static void StressTest(T & plasticModel, const char * filename, REAL stressMultiplier=1);
64 
65  template <class T>
66  static void StrainTest(T & plasticModel, const char * filename, REAL strainMultiplier=1);
67 
68  static void LoadTest(const char * filename);
69 
70  template <class T>
71  static int CreatePlasticModel(T * ( & plasticModel), const char * line);
72 
73  template <class T>
74  static void GlobalCheckConv(T & plasticModel, TPZTensor<REAL> & strain, REAL maxDeltaStrain = 0.01);
75 
76  static void LadeKimTriaxialLooseSand();
77 
78 
79  template <class T>
80  static void MultiDirectionsMaterialPointTest(T & plasticModel, REAL dirMult);
81 
82  static void DruckerPragerTest();
83 
84  static void MohrCoulombTest();
85 
86  static void ModifiedMohrCoulombTest();
87 
88  static void WillamWarnkeTest();
89 
90  static void VonMisesTest();
91 
92  static void UndocumentedTest2();
93 
94  static void UndocumentedTest3();
95 
96  static void UndocumentedTest4();
97 
99 
100  static void LKFineSilicaLoadTest();
101 
102  static void LKIsotropicCompression();
103 
104  static void LKKoCompressionLoadTest();
105 
106  static void LKLoadingTest();
107 
108  static void DruckerIsotropicCompression();
109 
110  static void DruckerTest();
111 
112  static void LKBiaxialTest();
113 
114  //static void MaterialPointTests();
116 
120  int NumCases()
121  {
122  return 9;
123  }
124 
126 
127  // TPZMatElastoPlastic<TPZMaterial> mate;
129 
135  {
136  int i;
137  for(i=0; i<6; i++) gRefTension.fData[i] = state(i,0);
138  }
139 
140  void ComputeTangent(TPZFMatrix<REAL> &tangent, TPZVec<REAL> &coefs, int icase)
141  {
142  switch(icase)
143  {
144  case 0:
145  {
146  TPZTensor<REAL> grad,EpsT,DiagonalStress;
147  TPZFNMatrix<6*6> Dep(6,6,0.);
148  gPlasticModel.ApplyStrainComputeDep(gRefTension,DiagonalStress,Dep);
149  tangent.Redim(1,1);
150  REAL norm = Norm(Dep);
151  tangent(0,0) = norm;
152 
153  break;
154  }
155 
156  }
157  }
158 
159  void Residual(TPZFMatrix<REAL> &res,int icase)
160  {
161 
162  res.Redim(1,1);
163  TPZTensor<REAL> grad,DiagonalStress;
164  TPZFNMatrix<36> Dep(6,6,0.);
165  gPlasticModel.ApplyStrainComputeDep(gRefTension,DiagonalStress,Dep);
166  res.Redim(1,1);
167  REAL norm = DiagonalStress.Norm();
168  switch(icase)
169  {
170  case 0:
171  {
172  res(0,0) = norm;
173  break;
174  }
175 
176  }
177 
178  }
179 
180 
181 
182 
183  static void RotateMatrix(TPZFMatrix<REAL> &Mat, REAL thetaRad,int rotateaboutaxes);
184 
185  //static void RotationMatrix(TPZFMatrix<REAL> &R, double thetaRad, int axis);
186 
187  template <class T>
188  static void PlasticIntegratorCheck(T mat);
189 
190  static void VerifyIntegrationAtPoint(TPZVec< TPZTensor<REAL> > vectensor);
191  // static void DruckerTest();
192 
193 
194 };
195 
197 {
198 
199  ofstream outfiletxt("LadeKimTriaxialLooseSand.txt");
200  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
201  TPZFNMatrix<6*6> Dep(6,6,0.);
202 
203 
204  TPZLadeKim LK;
206  LK.ApplyLoad(stress,strain);
207  deltastress.XX() = -0.005;//Confining stress = 17.07 Psi
208  deltastress.XY() = 0.;
209  deltastress.XZ() = 0.;
210  deltastress.YY() = -0.0003;//Confining stress = 17.07 Psi
211  deltastress.YZ() = 0.;
212  deltastress.ZZ() = -0.0003;//Confining stress = 17.07 Psi
213  stress=deltastress;
214 
215  LK.ApplyLoad(stress,strain);
216 
217 
218  int steps=60;
219  for(int i=0;i<steps;i++)
220  {
221 
222  LK.ApplyLoad(stress,strain);
223  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()-stress.ZZ()) << "\n";
224  stress+=deltastress;
225  }
226 
227 
228 
229 
230 }
231 //inline void MaterialPointTests()
232 //{
233 //
234 // cout << "\nChoose Plasticity test:";
235 // cout << "\n0 - Isotropic compression ";
236 // cout << "\n1 - Biaxial Tests ";
237 // cout << "\n2 - Uniaxial traction ";
238 // cout << "\n";
239 // int choice;
240 // cin >> choice;
241 //
242 // switch(choice)
243 // {
244 // case(0):
245 // cout << "\n Choose the Plastic model tou need to run Isotropic compression: ";
246 // cout << "\n0 - Lade - Kim ";
247 // cout << "\n1 - Sandler Dimaggio ";
248 // cout << "\n2 - Drucker Prager ";
249 // cout << "\n";
250 // int choice2;
251 // cin >> choice2;
252 // switch(choice2)
253 // {
254 // case(0):
255 // LKIsotropicCompression();
256 // break;
257 // case(1):
258 // SandlerDimaggioIsotropicCompression();
259 // break;
260 // case(2):
261 // DruckerIsotropicCompression();
262 // break;
263 // }
264 //
265 //
266 // break;
267 // case(1):
268 // LKBiaxialTest();
269 // break;
270 // case(2):
271 // cout << "NOT IMPLEMENTED YET";
272 // // LadeKim_ReversalTest();
273 // break;
274 // default:
275 // cout << "Unknown Test Type. Exiting...";
276 // }
277 //
278 //}
279 
280 
281 
283 {
284  ofstream outfiletxt5("SDMcCormicRanchSand.txt");
285  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
286  TPZFNMatrix<6*6> Dep(6,6,0.);
287 
288 
290 
291  cout << "\n Put the value of strain you want to add in each step of your loat test: (sugg.0.0001) ";
292  REAL straininput;
293  cin >> straininput;
294  deltastrain.XX() = -straininput;
295  deltastrain.XY() = 0.;
296  deltastrain.XZ() = 0.;
297  deltastrain.YY() = -straininput;
298  deltastrain.YZ() = 0.;
299  deltastrain.ZZ() = -straininput;
300  strain=deltastrain;
301  cout << "Choose the material pareameters you want to set to SandlerDimaggio Test :";
302  cout << "\n0 - McCormicRanchSandMod";
303  cout << "\n1 - McCormicRanchSandMod2 ";
304  cout << "\n2 - UncDeepSandRes ";
305  cout << "\n3 - UncDeepSandResPSI";
306  cout << "\n4 - UncDeepSandResMPa";
307  cout << "\n5 - Put the material parameters you want ";
308  int choice;
309  cin >> choice;
310 
311  cout << "\n Put the numbers of steps you want: (sugg. 20)";
312  int length;
313  cin >> length;
314 
315  switch (choice) {
316  case(0):
317  {
319  std::ofstream outfiletxt("TPZSandlerDimaggioMcCormicRanchSandMod(SD).txt");
320 
321  for(int step=0;step<length;step++)
322  {
323  cout << "\nstep "<< step;
324  SD.ApplyStrainComputeDep(strain, stress,Dep);
325  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
326  strain += deltastrain;
327  }
328  break;
329  }
330  case(1):
331  {
333  std::ofstream outfiletxt("TPZSandlerDimaggioMcCormicRanchSandMod2.txt");
334  for(int step=0;step<length;step++)
335  {
336  cout << "\nstep "<< step;
337  SD.ApplyStrainComputeDep(strain, stress,Dep);
338  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
339  strain += deltastrain;
340 
341  }
342  break;
343  }
344  case(2):
345  {
347  std::ofstream outfiletxt("TPZLadeKimUncDeepSandRes.txt");
348  for(int step=0;step<length;step++)
349  {
350  cout << "\nstep "<< step;
351  SD.ApplyStrainComputeDep(strain, stress,Dep);
352  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
353  strain += deltastrain;
354 
355  }
356  break;
357  }
358  case(3):
359  {
361  std::ofstream outfiletxt("TPZLadeKimUncDeepSandResPSI.txt");
362  for(int step=0;step<length;step++)
363  {
364  cout << "\nstep "<< step;
365  SD.ApplyStrainComputeDep(strain, stress,Dep);
366  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
367  strain += deltastrain;
368 
369  }
370  break;
371  }
372  case(4):
373  {
375  std::ofstream outfiletxt("TPZLadeKimUncDeepSandResMPa.txt");
376  for(int step=0;step<length;step++)
377  {
378  cout << "\nstep "<< step;
379  SD.ApplyStrainComputeDep(strain, stress,Dep);
380  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
381  strain += deltastrain;
382 
383  }
384  break;
385 
386  break;
387 
388  }
389  case(5):
390  {
391 
392  // REAL E = 9000,
393  // poisson = 0.25;
394  //
395  // material.fER.SetUp(E, poisson);
396  //
397  // REAL A = 18,
398  // B = 0.0245,
399  // C = 17.7,
400  // D = 0.00735,
401  // R = 1.5,
402  // W = 0.0908;
403  cout<< "\n Young Modulus 9000.";
404  REAL E;
405  cin >> E;
406 
407  cout<< "\n poisson 0.25 ";
408  REAL poisson;
409  cin >> poisson;
410 
411 
412  SD.fER.SetUp(E, poisson);
413 
414  cout<< "\n A (sugg.18)";
415  REAL A;
416  cin >> A;
417 
418  cout << "\n B (sugg. 0.0245) ";
419  REAL B;
420  cin >> B;
421 
422  cout << "\n C (sugg. 17.7) ";
423  REAL C;
424  cin >> C;
425 
426  cout << "\n D (sugg.0.00735) ";
427  REAL D;
428  cin >> D;
429 
430  cout << "\n R (sugg. 1.5) ";
431  REAL R;
432  cin >> R;
433 
434  cout << "\n W (sugg. 0.0908) ";
435  REAL W;
436  cin >> W;
437 
438  SD.fYC.SetUp(A, B, C, D, R, W);
439 
440  std::ofstream outfiletxt("SandlerDimaggioYOURMODEL.txt");
441  for(int step=0;step<length;step++)
442  {
443  cout << "\nstep "<< step;
444  SD.ApplyStrainComputeDep(strain, stress,Dep);
445  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
446  strain += deltastrain;
447 
448  }
449  break;
450  }
451  default:
452  {
453  cout << "Unknown Test Type. Exiting...";
454  break;
455  }
456  }
457 
458 
459 
460 }
461 
462 inline void LKFineSilicaLoadTest()//
463 {
464  ofstream outfiletxt1("FineSilica.txt");
465  TPZLadeKim LK;
467  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
468  TPZFNMatrix<6*6> Dep(6,6,0.);
469  deltastress.XX() = -4.;
470  deltastress.XY() = 0.;
471  deltastress.XZ() = 0.;
472  deltastress.YY() = -4.;
473  deltastress.YZ() = 0.;
474  deltastress.ZZ() = -4.;
475  stress = deltastress;
476 
477 
478  int length =30;
479  for(int step=0;step<length;step++)
480  {
481  cout << "\nstep "<< step;
482  if(step == 69)deltastress *=-1.;
483  LK.ApplyLoad(stress,strain);
484  outfiletxt1 << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
485  stress += deltastress;
486  cout << "strain = "<<strain <<"\n";
487  cout << "sigma = "<< stress <<"\t "<< "I1 = "<< stress.I1() <<"\n";
488 
489  }
490 
491 
492 
493 
494 }
495 
496 
497 
499 {
500 
501  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
502  TPZFNMatrix<6*6> Dep(6,6,0.);
503 
504  deltastress.XX() = -147.;
505  deltastress.XY() = 0.;
506  deltastress.XZ() = 0.;
507  deltastress.YY() = -147.;
508  deltastress.YZ() = 0.;
509  deltastress.ZZ() = -147.;
510  stress=deltastress;
511 
512 
513  TPZLadeKim LK2;
514  // TPZLadeKim::FineSilicaSandPaperIII(LK2);
515  // std::ofstream outfiletxt("TPZLadeKim_IsotropicCompression_FineSilicaSand.txt");
517  std::ofstream outfiletxt("TPZLadeKim_IsotropicCompression_PlaneConcrete.txt");
518 
519  int length =120;
520  for(int i=0;i<length;i++)
521  {
522  // if(i==9)deltastress*=-1;
523  cout << "\nstep = "<< i << endl;
524  LK2.ApplyLoad(stress,strain);
525  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX())/14.7 << "\n";
526  stress+=deltastress;
527  }
528 
529 }
530 
531 
532 //inline void LKIsotropicCompression()
533 //{
534 //
535 // TPZTensor<REAL> stress, strain, deltastress, deltastrain;
536 // TPZFNMatrix<6*6> Dep(6,6,0.);
537 //
538 // cout << "\n Put the value of strain you want to add in each step of your loat test: (sugg. 0.0001) ";
539 // REAL straininput = 0.00001;
540 // // cin >> straininput;
541 // /*
542 // deltastrain.XX() = -straininput;
543 // deltastrain.XY() = 0.;
544 // deltastrain.XZ() = 0.;
545 // deltastrain.YY() = -straininput;
546 // deltastrain.YZ() = 0.;
547 // deltastrain.ZZ() = -straininput;
548 // strain=deltastrain;
549 // */
550 // deltastrain.XX() = -0.0001;
551 // deltastrain.XY() = 0.;
552 // deltastrain.XZ() = 0.;
553 // deltastrain.YY() = -0.00005;
554 // deltastrain.YZ() = 0.;
555 // deltastrain.ZZ() = -0.;
556 // strain=deltastrain;
557 //
558 // cout << "Choose the material pareameters you want to set to Lade Kim Test :";
559 // cout << "\n0 - Plain Concrete ";
560 // cout << "\n1 - Loose Sacramento River Sand ";
561 // cout << "\n2 - Dense Sacramento River Sand ";
562 // cout << "\n3 - Fine Silica Sand";
563 // cout << "\n4 - Put the material parameters you want";
564 // int choice = 0;
565 // //cin >> choice;
566 //
567 // cout << "\n Put the numbers of steps you want:(sugg. 20)";
568 // int length =120;
569 // //cin >> length;
570 //
571 // TPZLadeKim LK2;
572 // switch (choice) {
573 // case(0):
574 // {
575 // TPZLadeKim::PlainConcrete(LK2);
576 // std::ofstream outfiletxt("TPZLadeKim::PlainConcrete.txt");
577 //
578 // for(int step=0;step<length;step++)
579 // {
580 // cout << "\nstep "<< step;
581 // LK2.ApplyStrainComputeDep(strain, stress,Dep);
582 // outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
583 // strain += deltastrain;
584 // }
585 // break;
586 // }
587 // case(1):
588 // {
589 // TPZLadeKim::LooseSacrRiverSand(LK2);
590 // std::ofstream outfiletxt("TPZLadeKim::LooseSacrRiverSand.txt");
591 // for(int step=0;step<length;step++)
592 // {
593 // cout << "\nstep "<< step;
594 // LK2.ApplyStrainComputeDep(strain, stress,Dep);
595 // outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
596 // strain += deltastrain;
597 //
598 // }
599 // break;
600 // }
601 // case(2):
602 // {
603 // TPZLadeKim::DenseSacrRiverSand(LK2);
604 // std::ofstream outfiletxt("TPZLadeKim::DenseSacrRiverSand.txt");
605 // for(int step=0;step<length;step++)
606 // {
607 // cout << "\nstep "<< step;
608 // LK2.ApplyStrainComputeDep(strain, stress,Dep);
609 // outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
610 // strain += deltastrain;
611 //
612 // }
613 // break;
614 // }
615 // case(3):
616 // {
617 // TPZLadeKim::FineSilicaSand(LK2);
618 // std::ofstream outfiletxt("TPZLadeKim::FineSilicaSand.txt");
619 // for(int step=0;step<length;step++)
620 // {
621 // cout << "\nstep "<< step;
622 // LK2.ApplyStrainComputeDep(strain, stress,Dep);
623 // outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
624 // strain += deltastrain;
625 //
626 // }
627 // break;
628 // }
629 // case(4):
630 // {
631 // cout << "\n poisson (sugg. 0.18)";
632 // REAL poisson;// = 0.18;
633 // cin>>poisson;
634 //
635 // cout << "\n M (sugg. 361800.)";
636 // REAL M;// = 361800.;
637 // cin >> M;
638 //
639 // cout << "\nlambda (sugg. 0.)";
640 // REAL lambda;// = 0.;
641 // cin >> lambda;
642 //
643 // cout << "\n a (sugg. 28.5)";
644 // REAL a;// = 28.5;
645 // cin >> a;
646 //
647 // cout << "\n m (sugg. 1.113)";
648 // REAL m;// = 1.113;
649 // cin >> m;
650 //
651 // cout << "\n neta1 (sugg. 159800.)";
652 // REAL neta1;// = 159800.;
653 // cin >> neta1;
654 //
655 // cout << "\n ksi2 (sugg. -2.92)";
656 // REAL ksi2; // = -2.92;
657 // cin >> ksi2;
658 //
659 // cout << "\n mu (sugg. 5.06)";
660 // REAL mu;// = 5.06;
661 // cin >> mu;
662 //
663 // cout << "\n C (sugg. 0.712E-12)";
664 // REAL C;// = 0.712E-12;
665 // cin >> C;
666 //
667 // cout << "\n p (sugg. 3.8)";
668 // REAL p;// = 3.8;
669 // cin >> p;
670 //
671 // cout <<"\n h (sugg. 1.990) ";
672 // REAL h;// = 1.990;
673 // cin >> h;
674 //
675 // cout << "\n alpha (sugg. 0.75) ";
676 // REAL alpha;// = 0.75;
677 // cin >> alpha;
678 //
679 // cout << "\n pa (sugg. 14.7) ";
680 // REAL pa;// = 14.7;
681 // cin >> pa;
682 //
683 // REAL restol;
684 // cout << "\n Tolerance (sugg. 0.0001) ";
685 // cin >> restol;
686 //
687 // LK2.fResTol = restol;
688 //
689 // LK2.SetUp(poisson, M, lambda,
690 // a, m, neta1,
691 // ksi2, mu,
692 // C, p,
693 // h, alpha,
694 // pa);
695 // std::ofstream outfiletxt("TPZLadeKim::YOURMODEL.txt");
696 // for(int step=0;step<length;step++)
697 // {
698 // cout << "\nstep "<< step;
699 // LK2.ApplyStrainComputeDep(strain, stress,Dep);
700 // outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
701 // strain += deltastrain;
702 //
703 // }
704 // break;
705 // }
706 // default:
707 // {
708 // cout << "Unknown Test Type. Exiting...";
709 // break;
710 // }
711 // }
712 //
713 //}
714 
715 
716 
718 {
719  ofstream outfiletxt1("LKKoCompressionLoadTest.txt");
720  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
721  TPZFNMatrix<6*6> Dep(6,6,0.);
722 
723  deltastress.XX()=-1.;
724  deltastress.YY()=-0.5;
725  deltastress.ZZ()=-0.5;
726  stress=deltastress;
727 
728 
729  TPZLadeKim LK2;
731 
732  int length2 =20;
733  for(int step=0;step<length2;step++)
734  {
735  cout << "\nstep "<< step;
736  LK2.ApplyLoad(stress,strain);
737  LK2.ApplyStrainComputeDep(strain, stress, Dep);
738  outfiletxt1 << fabs(strain.I1()) << " " << fabs(stress.XX()/stress.ZZ()) << "\n";
739  stress += deltastress;
740  cout << "strain = "<<strain <<"\n";
741  cout << "sigma = "<< stress <<"\t "<< "I1 = "<< stress.I1() <<"\n";
742 
743  }
744 
745 }
746 
747 
748 
749 
750 inline void LKLoadingTest()
751 {
752 
753  ofstream outfiletxt1("FineSilica.txt");
754  ofstream outfiletxt2("PlainConcretee1.txt");
755  ofstream outfiletxt3("PlainConcretee2.txt");
756  ofstream outfiletxt4("PlainConcretee3.txt");
757  ofstream outfiletxt5("SDMcCormicRanchSand.txt");
758  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
759  TPZFNMatrix<6*6> Dep(6,6,0.);
760 
761 
762 
763 
764  deltastress.XX() = -0.001;
765  deltastress.XY() = 0.;
766  deltastress.XZ() = 0.;
767  deltastress.YY() = -0.001;
768  deltastress.YZ() = 0.;
769  deltastress.ZZ() = -0.001;
770  stress = deltastress;
771 
772  TPZLadeKim LK;
773 
775  // LK.ApplyLoad(stress,deltastrain);
776 
777  deltastress.XX() = -4.;
778  deltastress.XY() = 0.;
779  deltastress.XZ() = 0.;
780  deltastress.YY() = -4.;
781  deltastress.YZ() = 0.;
782  deltastress.ZZ() = -4.;
783  stress = deltastress;
784 
785 
786  // deltastrain.XX() = -0.0001;
787  // deltastrain.XY() = 0.;
788  // deltastrain.XZ() = 0.;
789  // deltastrain.YY() = -0.0001;
790  // deltastrain.YZ() = 0.;
791  // deltastrain.ZZ() = -0.0001;
792  // strain=deltastrain;
793 
794  int length =72;
795  for(int step=0;step<length;step++)
796  {
797  cout << "\nstep "<< step;
798  if(step == 69)deltastress *=-1.;
799  LK.ApplyLoad(stress,strain);
800  // if(step == 16)deltastrain *=-1.;
801  // LK.ApplyStrainComputeDep(strain, stress,Dep);
802  if(step==0)
803  {
804 
805  //outfiletxt1 << 0. << " " << 0. << "\n";
806 
807  }
808 
809  else
810  {
811 
812  outfiletxt1 << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
813 
814  }
815 
816  // deltastress.Multiply(1.1, 1.);
817  stress += deltastress;
818  // strain +=deltastrain;
819  cout << "strain = "<<strain <<"\n";
820  cout << "sigma = "<< stress <<"\t "<< "I1 = "<< stress.I1() <<"\n";
821 
822  }
823 
824 
825  deltastrain.XX() = -0.0001;
826  deltastrain.XY() = 0.;
827  deltastrain.XZ() = 0.;
828  deltastrain.YY() = -0.;
829  deltastrain.YZ() = 0.;
830  deltastrain.ZZ() = -0.;
831  strain=deltastrain;
832  TPZLadeKim LK2;
834 
835  int length2 =30;
836  for(int step=0;step<length2;step++)
837  {
838  cout << "\nstep "<< step;
839  // if(step == 10 || step == 16|| step == 40 || step==51 || step==80)deltastrain *=-1.;
840  LK2.ApplyStrainComputeDep(strain, stress,Dep);
841  outfiletxt2 << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
842  outfiletxt3 << fabs(strain.YY()) << " " << fabs(stress.XX()) << "\n";
843  outfiletxt4 << fabs(strain.ZZ()) << " " << fabs(stress.XX()) << "\n";
844  strain += deltastrain;
845  cout << "strain = "<<strain <<"\n";
846  cout << "sigma = "<< stress <<"\t "<< "I1 = "<< stress.I1() <<"\n";
847 
848  }
849 
852 
853  deltastrain.XX() = -0.005;
854  deltastrain.XY() = 0.;
855  deltastrain.XZ() = 0.;
856  deltastrain.YY() = -0.;
857  deltastrain.YZ() = 0.;
858  deltastrain.ZZ() = -0.;
859  strain=deltastrain;
860 
861  int length3 =23;
862  for(int step=0;step<length3;step++)
863  {
864  cout << "\nstep "<< step;
865  if(step == 14 )deltastrain *=-1.;
866  SD.ApplyStrainComputeDep(strain, stress,Dep);
867  outfiletxt5 << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
868  strain += deltastrain;
869  cout << "strain = "<<strain <<"\n";
870  cout << "sigma = "<< stress <<"\t "<< "I1 = "<< stress.I1() <<"\n";
871 
872  }
873 
874 
875 }
876 
877 
878 
879 //int main()
880 //{
881 //
882 //
883 //
884 //
885 // PressureCilinder();
886 // return 0;
887 //}
888 
889 
891 {
892  TPZDruckerPrager DP;
893  ofstream outfiletxt("DruckerPragerIsotropicCompression.txt");
894  cout << "\n Put the value of strain you want to add in each step of your loat test: (sugg. 0.0001)";
895  REAL straininput;
896  cin >> straininput;
897  TPZFNMatrix<6*6> Dep(6,6,0.);
898  TPZTensor<REAL> deltastrain,strain,stress,deltastress;
899  deltastrain.XX() = -straininput;
900  deltastrain.XY() = 0.;
901  deltastrain.XZ() = 0.;
902  deltastrain.YY() = -straininput;
903  deltastrain.YZ() = 0.;
904  deltastrain.ZZ() = -straininput;
905  strain=deltastrain;
906 
907  cout << "\n4 - Put the material parameters you want";
908 
909  cout << "\n Young modulus (sugg. 20000.)";
910  REAL E;
911  cin >> E;
912 
913  cout << "\n Poisson (sugg. 0.2)";
914  REAL poisson;
915  cin >> poisson;
916 
917  int mcfit;
918  cout << "\n choose 0 for Iner Morh-Coulomb fit or 1 for outer Morh-Coulomb Fit (sugg. 0) ";
919  cin >> mcfit;
920 
921  // if(mcfit!= 0 || mcfit!= 1)
922  // {
923  // cout << "\n wrong choice in Morh-Coulomb fit tipe 0 or 1";
924  // return;
925  // }
926 
927  REAL phi;
928  cout << "\n Type the internal frictional angle in degrees(sugg. 20.)";
929  cin >> phi;
930 
931  REAL c;
932  cout << "\n Type the material coesion (sugg. 9.)";
933  cin >> c;
934 
935  REAL h;
936  cout << "\n Type the material hardening modulus (sugg. 1000.)";
937  cin >> h;
938 
939  DP.fYC.SetUp(phi/180. *M_PI ,mcfit);
940  DP.fTFA.SetUp(c,h);
941  DP.fER.SetUp(E,poisson);
942 
943  int length =100;
944  for(int step=0;step<length;step++)
945  {
946  cout << "\nstep "<< step;
947  DP.ApplyStrainComputeDep(strain, stress, Dep);
948  strain += deltastrain;
949  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
950  }
951 
952 }
953 
954 inline void DruckerBiaxialTest()
955 {
956  TPZDruckerPrager DP;
957  ofstream outfiletxt("DruckerBiaxialTest.txt");
958  cout << "\n Put the value of strain you want to add in each step of your loat test: (sugg. 0.0001)";
959  REAL straininput;
960  cin >> straininput;
961  TPZFNMatrix<6*6> Dep(6,6,0.);
962  TPZTensor<REAL> deltastrain,strain,stress,deltastress;
963  deltastrain.XX() = -straininput;
964  deltastrain.XY() = 0.;
965  deltastrain.XZ() = 0.;
966  deltastrain.YY() = -straininput/0.52;
967  deltastrain.YZ() = 0.;
968  deltastrain.ZZ() = 0.;
969  strain=deltastrain;
970 
971  cout << "\n4 - Put the material parameters you want";
972 
973  cout << "\n Young modulus (sugg. 20000.)";
974  REAL E;
975  cin >> E;
976 
977  cout << "\n Poisson (sugg. 0.2)";
978  REAL poisson;
979  cin >> poisson;
980 
981  int mcfit;
982  cout << "\n choose 0 for Iner Morh-Coulomb fit or 1 for outer Morh-Coulomb Fit (sugg. 0) ";
983  cin >> mcfit;
984 
985  // if(mcfit!= 0 || mcfit!= 1)
986  // {
987  // cout << "\n wrong choice in Morh-Coulomb fit tipe 0 or 1";
988  // return;
989  // }
990 
991  REAL phi;
992  cout << "\n Type the internal frictional angle in degrees(sugg. 20.)";
993  cin >> phi;
994 
995  REAL c;
996  cout << "\n Type the material coesion (sugg. 9.)";
997  cin >> c;
998 
999  REAL h;
1000  cout << "\n Type the material hardening modulus (sugg. 1000.)";
1001  cin >> h;
1002 
1003  DP.fYC.SetUp(phi/180. *M_PI ,mcfit);
1004  DP.fTFA.SetUp(c,h);
1005  DP.fER.SetUp(E,poisson);
1006 
1007  int length =100;
1008  for(int step=0;step<length;step++)
1009  {
1010  cout << "\nstep "<< step;
1011  DP.ApplyStrainComputeDep(strain, stress, Dep);
1012  strain += deltastrain;
1013  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
1014  }
1015 
1016 }
1017 
1018 
1019 inline void LKBiaxialTest()
1020 {
1021  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
1022  TPZFNMatrix<6*6> Dep(6,6,0.);
1023 
1024  cout << "\n Put the value of stress you want to put sx ";
1025  REAL stressinputx;
1026  cin >> stressinputx;
1027 
1028  cout << "\n Put the value of stress you want to put sy ";
1029  REAL stressinputy;
1030  cin >> stressinputy;
1031 
1032  cout << "\n Put the value of stress you want to put sz ";
1033  REAL stressinputz;
1034  cin >> stressinputz;
1035 
1036  deltastress.XX() = -stressinputx;
1037  deltastress.XY() = 0.;
1038  deltastress.XZ() = 0.;
1039  deltastress.YY() = -stressinputy;
1040  deltastress.YZ() = 0.;
1041  deltastress.ZZ() = -stressinputz;
1042  stress=deltastress;
1043 
1044  cout << "Choose the material pareameters you want to set to Lade Kim Test :";
1045  cout << "\n0 - Plain Concrete ";
1046  cout << "\n1 - Loose Sacramento River Sand ";
1047  cout << "\n2 - Dense Sacramento River Sand ";
1048  cout << "\n3 - Fine Silica Sand";
1049  cout << "\n4 - Put the material parameters you want";
1050  int choice;
1051  cin >> choice;
1052 
1053  cout << "\n Put the numbers of steps you want: (sugg. 20)";
1054  int length;
1055  cin >> length;
1056 
1057  TPZLadeKim LK2;
1058  switch (choice) {
1059  case(0):
1060  {
1062  std::ofstream outfiletxt1("BiaxialXTPZLadeKimPlainConcrete.txt");
1063  std::ofstream outfiletxt2("BiaxialYTPZLadeKimPlainConcrete.txt");
1064  std::ofstream outfiletxt3("BiaxialZTPZLadeKimPlainConcrete.txt");
1065  LK2.ApplyLoad(stress, strain);
1066  for(int step=0;step<length;step++)
1067  {
1068  cout << "\nstep "<< step;
1069  LK2.ApplyLoad(stress, strain);
1070  outfiletxt1 << -strain.XX() << " " << fabs(stress.XX()/14.7) << "\n";
1071  outfiletxt2 << -strain.YY() << " " << fabs(stress.XX()/14.7) << "\n";
1072  outfiletxt3 << -strain.ZZ() << " " << fabs(stress.XX()/14.7) << "\n";
1073  stress += deltastress;
1074  }
1075  break;
1076  }
1077  case(1):
1078  {
1080  std::ofstream outfiletxt1("BiaxialXTPZLadeKimLooseSacrRiverSand.txt");
1081  std::ofstream outfiletxt2("BiaxialYTPZLadeKimLooseSacrRiverSand.txt");
1082  std::ofstream outfiletxt3("BiaxialZTPZLadeKimLooseSacrRiverSand.txt");
1083  LK2.ApplyLoad(stress, strain);
1084  for(int step=0;step<length;step++)
1085  {
1086  cout << "\nstep "<< step;
1087  LK2.ApplyLoad(stress, strain);
1088  outfiletxt1 << strain.XX() << " " << fabs(stress.XX()) << "\n";
1089  outfiletxt2 << strain.YY() << " " << fabs(stress.XX()) << "\n";
1090  outfiletxt3 << strain.ZZ() << " " << fabs(stress.XX()) << "\n";
1091  stress += deltastress;
1092  }
1093  break;
1094  }
1095  case(2):
1096  {
1098  std::ofstream outfiletxt("TPZLadeKim::DenseSacrRiverSand.txt");
1099  LK2.ApplyLoad(stress, strain);
1100  for(int step=0;step<length;step++)
1101  {
1102  cout << "\nstep "<< step;
1103  LK2.ApplyStrainComputeDep(strain, stress,Dep);
1104  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
1105  outfiletxt << fabs(strain.YY()) << " " << fabs(stress.XX()) << "\n";
1106  outfiletxt << fabs(strain.ZZ()) << " " << fabs(stress.XX()) << "\n";
1107  strain += deltastrain;
1108 
1109  }
1110  break;
1111  }
1112  case(3):
1113  {
1115  std::ofstream outfiletxt("TPZLadeKim::FineSilicaSand.txt");
1116  LK2.ApplyLoad(stress, strain);
1117  for(int step=0;step<length;step++)
1118  {
1119  cout << "\nstep "<< step;
1120  LK2.ApplyStrainComputeDep(strain, stress,Dep);
1121  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
1122  outfiletxt << fabs(strain.YY()) << " " << fabs(stress.XX()) << "\n";
1123  outfiletxt << fabs(strain.ZZ()) << " " << fabs(stress.XX()) << "\n";
1124  strain += deltastrain;
1125 
1126  }
1127  break;
1128  }
1129  case(4):
1130  {
1131  cout << "\n poisson (sugg. 0.18)";
1132  REAL poisson;// = 0.18;
1133  cin>>poisson;
1134 
1135  cout << "\n M (sugg. 361800.)";
1136  REAL M;// = 361800.;
1137  cin >> M;
1138 
1139  cout << "\nlambda (sugg. 0.)";
1140  REAL lambda;// = 0.;
1141  cin >> lambda;
1142 
1143  cout << "\n a (sugg. 28.5)";
1144  REAL a;// = 28.5;
1145  cin >> a;
1146 
1147  cout << "\n m (sugg. 1.113)";
1148  REAL m;// = 1.113;
1149  cin >> m;
1150 
1151  cout << "\n neta1 (sugg. 159800.)";
1152  REAL neta1;// = 159800.;
1153  cin >> neta1;
1154 
1155  cout << "\n ksi2 (sugg. -2.92)";
1156  REAL ksi2; // = -2.92;
1157  cin >> ksi2;
1158 
1159  cout << "\n mu (sugg. 5.06)";
1160  REAL mu;// = 5.06;
1161  cin >> mu;
1162 
1163  cout << "\n C (sugg. 0.712E-12)";
1164  REAL C;// = 0.712E-12;
1165  cin >> C;
1166 
1167  cout << "\n p (sugg. 3.8)";
1168  REAL p;// = 3.8;
1169  cin >> p;
1170 
1171  cout <<"\n h (sugg. 1.990) ";
1172  REAL h;// = 1.990;
1173  cin >> h;
1174 
1175  cout << "\n alpha (sugg. 0.75) ";
1176  REAL alpha;// = 0.75;
1177  cin >> alpha;
1178 
1179  cout << "\n pa (sugg. 14.7) ";
1180  REAL pa;// = 14.7;
1181  cin >> pa;
1182 
1183  REAL restol;
1184  cout << "\n Tolerance (sugg. 0.0001) ";
1185  cin >> restol;
1186 
1187  LK2.SetResidualTolerance(restol);
1188 
1189  LK2.SetUp(poisson, M, lambda,
1190  a, m, neta1,
1191  ksi2, mu,
1192  C, p,
1193  h, alpha,
1194  pa);
1195  std::ofstream outfiletxt("TPZLadeKim::YOURMODEL.txt");
1196  for(int step=0;step<length;step++)
1197  {
1198  cout << "\nstep "<< step;
1199  LK2.ApplyStrainComputeDep(strain, stress,Dep);
1200  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
1201  outfiletxt << fabs(strain.YY()) << " " << fabs(stress.XX()) << "\n";
1202  outfiletxt << fabs(strain.ZZ()) << " " << fabs(stress.XX()) << "\n";
1203  strain += deltastrain;
1204 
1205  }
1206  break;
1207  }
1208  default:
1209  {
1210  cout << "Unknown Test Type. Exiting...";
1211  break;
1212  }
1213  }
1214 
1215 
1216 
1217 
1218 }
1219 
1220 
1221 //inline void MaterialPointTests()
1222 //{
1223 //
1224 // cout << "\nChoose Plasticity test:";
1225 // cout << "\n0 - Isotropic compression ";
1226 // cout << "\n1 - Biaxial Tests ";
1227 // cout << "\n2 - Uniaxial traction ";
1228 // cout << "\n";
1229 // int choice;
1230 // cin >> choice;
1231 //
1232 // switch(choice)
1233 // {
1234 // case(0):
1235 // cout << "\n Choose the Plastic model tou need to run Isotropic compression: ";
1236 // cout << "\n0 - Lade - Kim ";
1237 // cout << "\n1 - Sandler Dimaggio ";
1238 // cout << "\n2 - Drucker Prager ";
1239 // cout << "\n";
1240 // int choice2;
1241 // cin >> choice2;
1242 // switch(choice2)
1243 // {
1244 // case(0):
1245 // TPZPlasticTest::LKIsotropicCompression();
1246 // break;
1247 // case(1):
1248 // TPZPlasticTest::SandlerDimaggioIsotropicCompression();
1249 // break;
1250 // case(2):
1251 // TPZPlasticTest::DruckerIsotropicCompression();
1252 // break;
1253 // }
1254 //
1255 //
1256 // break;
1257 // case(1):
1258 // TPZPlasticTest::LKBiaxialTest();
1259 // break;
1260 // case(2):
1261 // cout << "NOT IMPLEMENTED YET";
1262 // // LadeKim_ReversalTest();
1263 // break;
1264 // default:
1265 // cout << "Unknown Test Type. Exiting...";
1266 // }
1267 //
1268 //}
1269 
1270 
1271 
1273 {
1274 #ifdef LOG4CXX
1275 
1276  std::string path;
1277  std::string configfile;
1278 #ifdef HAVE_CONFIG_H
1279  path = PLASTICITYSOURCEDIR;
1280  path += "/src/";
1281  cout << path.c_str() << endl;
1282  cout.flush();
1283 #else
1284  path = "";
1285 #endif
1286  configfile = path;
1287 
1288  configfile += "log4cxx.cfg";
1289  log4cxx::PropertyConfigurator::configure(configfile.c_str());
1290 
1291  std::stringstream sout;
1292  sout << __PRETTY_FUNCTION__ << "\nLOG4CXX configured.\n"
1293  << "LOG4CXX config file:" << configfile;
1294 
1295  LOGPZ_INFO(plasticIntegrLogger,sout.str().c_str());
1296 
1297 #ifdef LOG4CXX_PLASTICITY
1298  LOGPZ_INFO(testLogger,sout.str().c_str());
1299 #endif
1300 
1301 #endif
1302 }
1303 
1304 
1305 template <class T>
1306 inline void TPZPlasticTest::ReciprocityTest(T & plasticModel, TPZTensor<REAL> strain1)
1307 {
1308 
1310 
1311  T plasticModelCopy = plasticModel;
1312 
1313  TPZTensor<REAL> stress1, stress2, strain2;
1314  TPZFNMatrix < 6*6 > tangent (6,6,0.);
1315 
1316  TPZVec<REAL> phi(2,0.), phi2(2,0.);
1317 
1318  plasticModel.Phi(strain1,phi);
1319 
1320 #ifdef LOG4CXX_PLASTICITY
1321  {
1322  std::stringstream sout;
1323  sout << __PRETTY_FUNCTION__
1324  << "\nImposing strain = \n" << strain1
1325  << "\nPhi before plastic loop = \n" << phi ;
1326  LOGPZ_INFO(testLogger,sout.str().c_str());
1327 
1328  }
1329 #endif
1330 
1331  // applying the strain to the element and processing plastic loop
1332 
1333  plasticModel.ApplyStrainComputeDep(strain1, stress1, tangent);
1334 
1335  plasticModel.Phi(strain1,phi);
1336 
1337 #ifdef LOG4CXX_PLASTICITY
1338  {
1339  std::stringstream sout;
1340  sout << __PRETTY_FUNCTION__
1341  << "\nphi after plastic loop =\n" << phi
1342  << "\ncausing stress = " << stress1
1343  << "\nNow applying the opposite operation...";
1344  LOGPZ_INFO(testLogger,sout.str().c_str());
1345 
1346  }
1347 #endif
1348 
1349  stress1.CopyTo(stress2);
1350 
1351  plasticModelCopy.ApplyLoad(stress2, strain2);
1352 
1353  plasticModelCopy.Phi(strain2,phi2);
1354 
1355 #ifdef LOG4CXX_PLASTICITY
1356  {
1357  std::stringstream sout;
1358  sout << __PRETTY_FUNCTION__
1359  << "\nphi after plastic loop =\n" << phi2
1360  << "\ncausing strain = " << strain2
1361  << "\nEnd of test";
1362  LOGPZ_INFO(testLogger,sout.str().c_str());
1363 
1364  }
1365 #endif
1366 
1367 }
1368 
1369 template <class T>
1370 inline void TPZPlasticTest::StressTest(T & plasticModel, const char * filename, REAL stressMultiplier)
1371 {
1372 
1373  // Load test: in this test the material is submitted to a stress loading
1374  // path read from a file in the format below.
1375 
1377 
1378 #ifdef LOG4CXX_PLASTICITY
1379  {
1380  std::stringstream sout;
1381  sout << __PRETTY_FUNCTION__
1382  << "\nApply Load Test\nPlease enter file name containing load path (sigmaxx xy xz yy yz zz) for all iterations:\n"
1383  << "\nReading file named \"" << filename << "\"\n";
1384  LOGPZ_INFO(testLogger,sout.str().c_str());
1385  }
1386 #endif
1387 
1388  const int linelen = 1024;
1389  char outfilename[linelen], line[linelen];
1390 
1391  ifstream file;
1392  ofstream outFile;
1393 
1394  file.open(filename);
1395 
1396  if(!file.is_open())
1397  {
1398  cout << "\nFile not open.\nExiting...\n";
1399  return;
1400  }
1401 
1402  file.getline(line, linelen);
1403  strncpy(outfilename, filename, 120);
1404  strcpy(outfilename+strlen(outfilename), ".out");
1405 
1406 #ifdef LOG4CXX_PLASTICITY
1407  {
1408  std::stringstream sout;
1409  sout << __PRETTY_FUNCTION__
1410  << "\nReading file named \"" << filename << "\"\n"
1411  << "\nTest data description: \"" << line << "\"\n"
1412  << "\nOutput file named \"" << outfilename << "\"\n";
1413  LOGPZ_INFO(testLogger,sout.str().c_str());
1414  }
1415 #endif
1416 
1417  outFile.open(outfilename);
1418  outFile << line << endl;
1419 
1420  TPZManVector<TPZTensor<REAL>, 2> loadPath, totalStrainPath;
1421  int i, j=0;
1422 
1423  while(file.good())
1424  {
1425  file.getline(line, linelen);
1426  stringstream strLine(line);
1427  loadPath.Resize(j+1);
1428  for(i = 0; i < 6; i++)
1429  {
1430  strLine >> loadPath[j].fData[i];
1431  }
1432 
1433  j++;
1434  }
1435  int nsteps = loadPath.NElements();
1436  totalStrainPath.Resize(nsteps);
1437 
1438  for(i = 0; i < nsteps; i++)
1439  {
1440  TPZTensor<REAL> tempLoad(loadPath[i]);
1441  for(j=0;j<6;j++)tempLoad.fData[j]*=stressMultiplier;
1442  if(i > 1)totalStrainPath[i] = totalStrainPath[i-1];
1443 
1444  plasticModel.ApplyLoad(tempLoad, totalStrainPath[i]);
1445 
1446  std::stringstream outputLine;
1447  outputLine << "step " << i << ", sigma: " << loadPath[i] ;
1448  outputLine << "\nfN:/n" << plasticModel.GetState() << endl;
1449 
1450  outFile << outputLine.str();
1451 #ifdef LOG4CXX_PLASTICITY
1452  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1453 #endif
1454  outFile.flush();
1455  }
1456 
1457  file.close();
1458  outFile.close();
1459 
1460 }
1461 
1462 
1463 template <class T>
1464 inline void TPZPlasticTest::StrainTest(T & plasticModel, const char * filename, REAL strainMultiplier)
1465 {
1466 
1467  // Strain test: in this test the material is submitted to a strain
1468  // path according to that read from a file in the format below.
1469 
1471 
1472 #ifdef LOG4CXX_PLASTICITY
1473  {
1474  std::stringstream sout;
1475  sout << __PRETTY_FUNCTION__
1476  << "\nApply Strain Test\nPlease enter file name containing strain path (sigmaxx xy xz yy yz zz) for all iterations:\n"
1477  << "\nReading file named \"" << filename << "\"\n";
1478  LOGPZ_INFO(testLogger,sout.str().c_str());
1479  }
1480 #endif
1481 
1482  const int linelen = 1024;
1483  char outfilename[linelen], line[linelen];
1484 
1485  ifstream file;
1486  ofstream outFile;
1487 
1488  file.open(filename);
1489 
1490  if(!file.is_open())
1491  {
1492  cout << "\nFile not open.\nExiting...\n";
1493  return;
1494  }
1495 
1496  file.getline(line, linelen);
1497  strncpy(outfilename, filename, 120);
1498  strcpy(outfilename+strlen(outfilename), ".out");
1499 
1500 #ifdef LOG4CXX_PLASTICITY
1501  {
1502  std::stringstream sout;
1503  sout << __PRETTY_FUNCTION__
1504  << "\nReading file named \"" << filename << "\"\n"
1505  << "\nTest data description: \"" << line << "\"\n"
1506  << "\nOutput file named \"" << outfilename << "\"\n";
1507  LOGPZ_INFO(testLogger,sout.str().c_str());
1508  }
1509 #endif
1510 
1511  outFile.open(outfilename);
1512  outFile << line << endl;
1513 
1514  TPZManVector<TPZTensor<REAL>, 2> strainPath, stressPath;
1515  TPZFNMatrix < 6*6 > tangent (6,6,0.); // satisfying interface only
1516  int i, j=0;
1517 
1518  while(file.good())
1519  {
1520  file.getline(line, linelen);
1521  stringstream strLine(line);
1522  strainPath.Resize(j+1);
1523  for(i = 0; i < 6; i++)
1524  {
1525  strLine >> strainPath[j].fData[i];
1526  }
1527 
1528  j++;
1529  }
1530  int nsteps = strainPath.NElements() - 1;// the last line loaded always return null path
1531  stressPath.Resize(nsteps);
1532 
1533  for(i = 0; i < nsteps; i++)
1534  {
1535  TPZTensor<REAL> tempStrain(strainPath[i]);
1536  for(j=0;j<6;j++)tempStrain.fData[j]*=strainMultiplier;
1537  if(i > 1)stressPath[i] = stressPath[i-1];
1538 
1539  plasticModel.ApplyStrain(tempStrain);
1540  plasticModel.ApplyStrainComputeDep(tempStrain, stressPath[i], tangent);
1541 
1542  TPZPlasticState<REAL> state = plasticModel.GetState();
1543  TPZTensor<REAL> epsp(state.EpsP());
1544  REAL alpha=state.VolHardening();
1545 
1546 
1547  // TPZTensor<REAL> tempStrain(strainPath[i]);
1548  // for(j=0;j<6;j++)tempStrain.fData[j]*=strainMultiplier;
1549  // if(i > 1)stressPath[i] = stressPath[i-1];
1550  // plasticModel.ApplyStrain(tempStrain);
1551  // plasticModel.Sigma(tempStrain, stressPath[i], tangent);
1552  // TPZTensor<REAL> epsp;
1553  // plasticModel.GetPlasticStrain(epsp);
1554  // REAL alpha = plasticModel.GetAlpha();
1555 
1556 
1557 
1558  std::stringstream outputLine;
1559  outputLine << "step " << i << ", imposedEps: " << strainPath[i]
1560  << ", sigma: " << stressPath[i]
1561  << ", epsP= " << epsp
1562  << ", alpha= " << alpha << endl;
1563 
1564  outFile << outputLine.str();
1565 #ifdef LOG4CXX_PLASTICITY
1566  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1567 #endif
1568  outFile.flush();
1569  }
1570 
1571  file.close();
1572  outFile.close();
1573 
1574 }
1575 
1577 
1578 inline void TPZPlasticTest::LoadTest(const char * filename)
1579 {
1580 
1581  // Load test: in this test the material is submitted to a loading path
1582  // read from a file in the format below. The keywords strain and stress
1583  // switches to the proper loading conditions.
1584 
1586 
1587 #ifdef LOG4CXX_PLASTICITY
1588  {
1589  std::stringstream sout;
1590  sout << __PRETTY_FUNCTION__
1591  << "\nApply Load Test\nPlease enter file name containing load path (strain/stress sigmaxx xy xz yy yz zz) for all iterations:\n"
1592  << "\nReading file named \"" << filename << "\"\n";
1593  LOGPZ_INFO(testLogger,sout.str().c_str());
1594  }
1595 #endif
1596 
1597  const int linelen = 1024;
1598  char outfilename[linelen], line[linelen];
1599 
1600  ifstream file;
1601  ofstream outFile;
1602 
1603  file.open(filename);
1604 
1605  if(!file.is_open())
1606  {
1607  cout << "\nFile not open.\nExiting...\n";
1608  return;
1609  }
1610 
1612  file.getline(line, linelen);
1613  strncpy(outfilename, filename, 120);
1615  strcpy(outfilename+strlen(outfilename), ".txt");
1616 
1617 #ifdef LOG4CXX_PLASTICITY
1618  {
1619  std::stringstream sout;
1620  sout << __PRETTY_FUNCTION__
1621  << "\nReading file named \"" << filename << "\"\n"
1622  << "\nTest data description: \"" << line << "\"\n"
1623  << "\nOutput file named \"" << outfilename << "\"\n";
1624  LOGPZ_INFO(testLogger,sout.str().c_str());
1625  }
1626 #endif
1627 
1628  outFile.open(outfilename);
1629  outFile << line << endl;
1630 
1631  //reading second line description (CVS label data for instance)
1632  file.getline(line, linelen);
1633  outFile << line << endl;
1634 
1635  //reading plasticity model
1636  file.getline(line, linelen);
1637  TPZPlasticBase * pPlasticModel;
1638 
1639  if(!CreatePlasticModel(pPlasticModel, line))
1640  {
1641  std::stringstream sout;
1642  sout << "Could not create plastic model named " << line;
1643 #ifdef LOG4CXX_PLASTICITY
1644  {
1645  LOGPZ_INFO(testLogger,sout.str().c_str());
1646  }
1647 #endif
1648 
1649  outFile << endl << sout.str();
1650  outFile.close();
1651  file.close();
1652  return;
1653  }
1654 
1655 
1656  TPZManVector<TPZTensor<REAL>, 2> loadPath, totalStrainPath;
1657  TPZTensor<REAL> stress, tempStress, strain, tempStrain, epsp, strainP;
1658  //TPZFNMatrix < 6*6 > tangent (6,6,0.); // satisfying interface only
1659  REAL stressMultiplier = 1., strainMultiplier = 1., integrTol;
1660  int i, j = 0;
1661 
1662  while(file.good())
1663  {
1664  file.getline(line, linelen);
1665  if(!strncmp(line, "strainLoad",10))
1666  {
1667  stringstream strLine(line + 10);
1668  for(i = 0; i < 6; i++)
1669  {
1670  strLine >> strain.fData[i];
1671  }
1672  strain.CopyTo(tempStrain);
1673  for(i = 0; i < 6; i++)tempStrain.fData[i]*=strainMultiplier;
1674 
1675  pPlasticModel->ApplyStrainComputeSigma(tempStrain, tempStress);
1676 
1677  tempStress.CopyTo(stress);
1678  for(i = 0; i < 6; i++)stress.fData[i]/=stressMultiplier;
1679 
1680  pPlasticModel->GetState().m_eps_p.CopyTo(strainP);
1681  for(i = 0; i < 6; i++)strainP.fData[i]/=strainMultiplier;
1682 
1683  std::stringstream outputLine;
1684 
1685  // outputLine << "stress step " << j
1686  // << ", strain: " << strain.XX()
1687  // << ", stress: " << stress.XX()
1688  // << ", epsP = " << strainP
1689  // << ", alpha = " << pPlasticModel->GetState().m_hardening
1690  // << ", integrationSteps = " << pPlasticModel->IntegrationSteps() << endl;
1691  cout
1692  << " " << strain.XX()
1693  << " " << stress.XX() << endl;
1694 
1695 
1696  outputLine
1697  << " " << strain.XX()
1698  << " " << stress.XX() << endl;
1699 
1700  outFile << outputLine.str();
1701  outFile.flush();
1702 #ifdef LOG4CXX_PLASTICITY
1703  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1704 #endif
1705  j++;
1706  }
1707  if(!strncmp(line, "stressLoad",10))
1708  {
1709  stringstream strLine(line + 10);
1710  for(i = 0; i < 6; i++)
1711  {
1712  strLine >> stress.fData[i];
1713  }
1714  stress.CopyTo(tempStress);
1715  for(i = 0; i < 6; i++)tempStress.fData[i]*=stressMultiplier;
1716 
1717  pPlasticModel->ApplyLoad(tempStress, tempStrain);
1718 
1719  int intSteps = 0;
1720  intSteps = pPlasticModel->IntegrationSteps();
1721 
1722  //evaluating the converged sigma
1723  pPlasticModel->ApplyStrainComputeSigma(tempStrain, tempStress);
1724 
1725  tempStrain.CopyTo(strain);
1726  for(i = 0; i < 6; i++)strain.fData[i]/=strainMultiplier;
1727 
1728  pPlasticModel->GetState().m_eps_p.CopyTo(strainP);
1729  for(i = 0; i < 6; i++)strainP.fData[i]/=strainMultiplier;
1730 
1731  tempStress.CopyTo(stress);
1732  for(i = 0; i < 6; i++)stress.fData[i]/=stressMultiplier;
1733 
1734 
1735  std::stringstream outputLine;
1736 
1737 
1738  // outputLine << "stress step " << j
1739  // << ", strain: " << strain.XX()
1740  // << ", stress: " << stress.XX()
1741  // << ", epsP = " << strainP
1742  // << ", alpha = " << pPlasticModel->GetState().m_hardening
1743  // << ", integrationSteps = " << intSteps << endl;
1744 
1745  cout
1746  << " " << strain.XX()
1747  << " " << stress.XX() << endl;
1748 
1749 
1750  outputLine
1751  << " " << strain.XX()
1752  << " " << stress.XX() << endl;
1753 
1754  outFile << outputLine.str();
1755  outFile.flush();
1756 #ifdef LOG4CXX_PLASTICITY
1757  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1758 #endif
1759  j++;
1760  }
1761  if(!strncmp(line, "stressMult",10))
1762  {
1763  stringstream strLine(line + 10);
1764  strLine >> stressMultiplier;
1765 
1766  std::stringstream outputLine;
1767  outputLine << "stressMult = " << stressMultiplier << endl;
1768 
1769  outFile << outputLine.str();
1770  outFile.flush();
1771 #ifdef LOG4CXX_PLASTICITY
1772  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1773 #endif
1774  }
1775  if(!strncmp(line, "strainMult",10))
1776  {
1777  stringstream strLine(line + 10);
1778  strLine >> strainMultiplier;
1779 
1780  std::stringstream outputLine;
1781  outputLine << "strainMult = " << strainMultiplier << endl;
1782 
1783  outFile << outputLine.str();
1784  outFile.flush();
1785 #ifdef LOG4CXX_PLASTICITY
1786  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1787 #endif
1788  }
1789  if(!strncmp(line, "resetMater",10))
1790  {
1791  TPZPlasticState<REAL> nullState;
1792 
1793  pPlasticModel->SetState(nullState);
1794 
1795  std::stringstream outputLine;
1796  outputLine << "ResetMaterial" << endl;
1797 
1798  outFile << outputLine.str();
1799  outFile.flush();
1800 #ifdef LOG4CXX_PLASTICITY
1801  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1802 #endif
1803  }
1804  if(!strncmp(line, "integrTol ",10))
1805  {
1806  stringstream strLine(line + 10);
1807  strLine >> integrTol;
1808  pPlasticModel->SetIntegrTol(integrTol);
1809 
1810  std::stringstream outputLine;
1811  outputLine << "integrTol " << integrTol << endl;
1812 
1813  outFile << outputLine.str();
1814  outFile.flush();
1815 #ifdef LOG4CXX_PLASTICITY
1816  LOGPZ_INFO(testLogger,outputLine.str().c_str());
1817 #endif
1818  }
1819  }
1820 
1821  file.close();
1822  outFile.close();
1823 
1824  delete pPlasticModel;
1825 
1826 }
1827 
1828 template <class T>
1829 inline int TPZPlasticTest::CreatePlasticModel(T * ( & pPlasticModel), const char * line)
1830 {
1831  if(!strncmp(line,"LadeKim.PlainConcrete",21))
1832  {
1833  TPZLadeKim * pLK = new TPZLadeKim();
1835  pPlasticModel = pLK;
1836  return 1;
1837  }
1838  if(!strncmp(line,"LadeKim.LooseSacrRiverSand",26))
1839  {
1840  TPZLadeKim * pLK = new TPZLadeKim();
1842  pPlasticModel = pLK;
1843  return 1;
1844  }
1845  if(!strncmp(line,"LadeKim.DenseSacrRiverSand",26))
1846  {
1847  TPZLadeKim * pLK = new TPZLadeKim();
1849  pPlasticModel = pLK;
1850  return 1;
1851  }
1852  if(!strncmp(line,"LadeKim.FineSilicaSand",22))
1853  {
1854  TPZLadeKim * pLK = new TPZLadeKim();
1856  pPlasticModel = pLK;
1857  return 1;
1858  }
1859  if(!strncmp(line,"SandlerDimaggio.McCormicRanchSandMod",36))
1860  {
1863  pPlasticModel = pSD;
1864  return 1;
1865  }
1866  if(!strncmp(line,"SandlerDimaggio.McCormicRanchSand",33))
1867  {
1870  pPlasticModel = pSD;
1871  return 1;
1872  }
1873  if(!strncmp(line,"DruckerPrager.DummyConcrete",27))
1874  {
1876  TPZDruckerPrager * pDP = new TPZDruckerPrager();
1877  pDP->fYC.SetUp(/*phi=30*/asin(0.5),/*innerMCFit*/0);
1878  pDP->fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ 9.2376, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1000.);
1879  pDP->fER.SetUp(/*young*/ 20000, /*poisson*/ 0.2);
1880  pPlasticModel = pDP;
1881  return 1;
1882  }
1883  if(!strncmp(line,"MohrCoulomb.DummyConcrete",25))
1884  {
1886  TPZDruckerPrager * pDP = new TPZDruckerPrager();
1887  pDP->fYC.SetUp(/*phi=30*/asin(0.5),/*innerMCFit*/0);
1888  pDP->fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ 9.2376, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1000.);
1889  pDP->fER.SetUp(/*young*/ 20000, /*poisson*/ 0.2);
1890  pPlasticModel = pDP;
1891  return 1;
1892  }
1893 
1894 
1895  return 0;
1896 }
1897 
1898 template <class T>
1899 inline void TPZPlasticTest::GlobalCheckConv(T & plasticModel, TPZTensor<REAL> & strain, REAL maxDeltaStrain)
1900 {
1901 
1902  // GlobalCheckConv Test: tests if the elastoplastic evaluated stiffness matrix
1903  // is really the jacobian of the stress tensor with respect to the total
1904  // strain.
1905  // The imported plasticModel is submitted to the imposed strain. The evaluated stress
1906  // and jacobian matrix are kept in memory.
1907  // Sucessive copies of this plastic model are subjected to further strains and
1908  // the check conv method is evaluated.
1909 
1910  const int nRetries = 30; // must be at least 2
1911  const int nVars = 6;
1912  TPZFNMatrix<nVars*nVars> Dep(nVars,nVars), tempMatrix(nVars,nVars);
1913  TPZTensor<REAL> stress, tempStress;
1914  //The below matrices hold the values of the stress vector (column wise) due to the
1915  //strain component changes.
1916  TPZManVector<TPZFNMatrix<nVars*nVars>, nRetries> FcnHistory(nRetries);
1917  // The below strain tensors hold the changes in strain components.
1918  TPZManVector<TPZTensor<REAL>, nRetries> StrainHistory(nRetries);
1919 
1920  plasticModel.ApplyStrainComputeDep(strain, stress, Dep);
1921 
1922 
1923  //TPZPlasticTest::InitializeLOG();
1924 
1925 #ifdef LOG4CXX
1926  {
1927  std::stringstream sout;
1928  sout << __PRETTY_FUNCTION__
1929  << "\nGlobalCheckConv Test with " << nRetries << " samplings\n";
1930  sout << "\n Stress = " << stress
1931  << "\n Dep=\n" << Dep;
1932  LOGPZ_INFO(plasticIntegrLogger,sout.str().c_str());
1933  }
1934 #endif
1935 
1936  int i,j,k;
1937 
1938 
1939  //Evaluating the plastic integration for all cases.
1940  for(k = 0; k < nRetries; k++)
1941  {
1942  FcnHistory[k].Resize(nVars,nVars);
1943  for(i = 0; i < nVars; i++)
1944  {
1945  T plasticModelCopy(plasticModel);
1946  REAL random = (rand()&1000)/1000.;
1947  REAL deltaStrain = random * maxDeltaStrain;
1948  TPZTensor<REAL> tempStrain;
1949  tempStrain.fData[i] = deltaStrain;
1950  StrainHistory[k].fData[i] = deltaStrain;
1951  tempStrain.Add(strain,1.);
1952  plasticModelCopy.ApplyStrainComputeDep(tempStrain, tempStress, tempMatrix);
1953  for(j = 0; j < nVars; j++)
1954  FcnHistory[k](j,i) = tempStress.fData[j]
1955  - stress.fData[j]
1956  - /*Dep(j,i)*/tempMatrix(j,i) * deltaStrain;
1957 
1958  /*
1959  for(j = 0; j < nVars; j++)
1960  if(i==j)
1961  {
1962  FcnHistory[k](j,i) = pow(-strain.fData[j],1.5) +
1963  1.5 * pow(-strain.fData[j],0.5) * (-deltaStrain)
1964  - pow(-tempStrain.fData[j],1.5);
1965  }else{
1966  FcnHistory[k](j,i) = pow(-strain.fData[j],1.5) +
1967  1.5 * pow(-strain.fData[j],0.5) * 0.
1968  - pow(-tempStrain.fData[j],1.5);
1969  }
1970  */
1971 
1972 #ifdef LOG4CXX
1973  {
1974  std::stringstream sout;
1975  sout << "\nCase " << k
1976  << "\nStrain= " << tempStrain
1977  << "\nStress= " << tempStress
1978  << "\ntempMatrix= " << tempMatrix;
1979  LOGPZ_INFO(plasticIntegrLogger,sout.str().c_str());
1980  }
1981 #endif
1982  }
1983 
1984  }
1985 
1986  // Evaluating the checkConv routines
1987 
1988  std::stringstream output;
1989  for(k = 1; k < nRetries; k++)
1990  {
1991  output << "\n\n Result of checkconv with steps " << k-1 << " and " << k << "\n";
1992  output << "\nStrainHistory[" << k-1 << "]=" << StrainHistory[k-1];
1993  output << "\nFcnHistory[" << k-1 << "]=" << FcnHistory[k-1];
1994  output << "\nStrainHistory[" << k << "]=" << StrainHistory[k];
1995  output << "\nFcnHistory[" << k << "]=" << FcnHistory[k];
1996  output << "\nCheckConv:\n" ;
1997  for(j = 0; j < nVars; j++)
1998  {
1999  for(i = 0; i < nVars; i++)
2000  {
2001  REAL ykm1 = FcnHistory[k-1](j,i);
2002  REAL yk = FcnHistory[k ](j,i);
2003  REAL xkm1 = StrainHistory[k-1].fData[i];
2004  REAL xk = StrainHistory[k ].fData[i];
2005 
2006  output.width(12);
2007  if(ykm1 * yk < 1.e-24)
2008  {
2009  if(ykm1 * yk < -1.e-24)
2010  {
2011  output << "+/-y?"; // They should be of the same sign
2012  }else
2013  {
2014  output << "Exact"; // at least one of the evaluations was exact
2015  }
2016 
2017  }else
2018  {
2019  if(xkm1 * xk < 1.e-24)
2020  {
2021  output << "+/-x?";
2022  }else
2023  {
2024  output << log(ykm1/yk) / log(xkm1 / xk);
2025  }
2026  }
2027  }
2028  output << "\n";
2029  }
2030  output << "\n" ;
2031  }
2032 
2033 #ifdef LOG4CXX
2034  {
2035  // LOGPZ_INFO(testLogger,output.str().c_str());
2036  LOGPZ_INFO(plasticIntegrLogger,output.str().c_str());
2037  }
2038 #endif
2039 
2040 }
2041 
2043 {
2044 
2046  ofstream outfile("comparamathematica1.nb");
2047 
2048  int choice;
2049  choice =0;
2050  cout << "\nChoose Load Case:";
2051  cout << "\n0 - Strain Step";
2052  cout << "\n1 - Stress Step";
2053  cout << "\n";
2054 
2055  //cin >> choice;
2056  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
2057 
2058  //deltastress.XX() = -1.;
2059  // deltastress.YY() = 0.;
2060  // deltastress.ZZ()= 0.;
2061  // deltastress.XY()= 0.;
2062  // deltastress.XZ()= 0.;
2063  // deltastress.YZ()= 0.;
2064 
2065  deltastress.XX() = -10.;
2066  deltastress.ZZ() = 1.;
2067  stress = deltastress;
2068 
2069  /*
2070  deltastrain.XX() = -0.0001;
2071  // deltastrain.XY() = 0.0001;
2072  // deltastrain.XZ() = -0.0001;
2073  // deltastrain.YZ() = -0.0001;
2074  deltastrain.ZZ() = 0.00002;
2075  strain=deltastrain;
2076  */
2077 
2079  TPZDruckerPrager Pstep;
2080  REAL pi = M_PI;
2081  /*innerMCFit = 0*/
2082  /*OuterMCFit = 1*/
2083  Pstep.fYC.SetUp(/*phi=20*/ 20./180. * pi ,/*MCFit*/1);
2084  REAL coesao = 9.2376;
2085  Pstep.fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ coesao, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1000.);
2086  Pstep.fER.SetUp(/*young*/ 20000., /*poisson*/ 0.2);
2087 
2088  int length = 300;
2089  for(int step=0;step<length;step++)
2090  {
2091  cout << "\nstep "<< step;
2092 
2093 
2094  if( step==100 || step==200)
2095  {
2096  deltastrain *= -1.;
2097  deltastress *= -1.;
2098  }
2099 
2100  stress += deltastress;
2101  Pstep.ApplyLoad(stress,strain);
2102 
2103 
2104  // strain += deltastrain;
2105  // Pstep.ApplyStrainComputeSigma(strain,stress);
2106 
2107  if(step==0)outfile <<"points={";
2108 
2109  if(step!=length && step!=length-1)
2110  {
2111  outfile <<"{"<<-strain.XX()<< ","<<-stress.XX()<<" }, ";
2112  //outfile <<"{"<<-stress.J2()<< ","<<-stress.I1()<<" }, ";
2113  //outfile <<"{"<<stress.ZZ()<< ","<<stress.XX()<<" }, ";
2114  }
2115 
2116  if(step==length-2)outfile <<"{"<<-strain.XX()<< ","<<-stress.XX()<<"}";
2117  //if(step==length-2)outfile <<"{"<<-stress.J2()<< ","<<-stress.I1()<<"}";
2118  //if(step==length-2)outfile <<"{"<<stress.ZZ()<< ","<<stress.XX()<<"}";
2119 
2120  if(step==length-1)
2121  {
2122  //do = SolAprox = ListPlot[points, PlotRange -> All, AxesLabel -> {epsilonx, sigmax}, PlotLabel -> DRUCKER*PRAGER*OUTERMCFIT, AxesOrigin -> {0, 0},PlotStyle -> {Thick, Blue, Dashed}, Joined -> True, FillingStyle -> Directive[Opacity[0.5], Orange]]
2123  outfile << " };\n";
2124  outfile<<"SolAprox = ListPlot[points, PlotRange-> All,AxesLabel->{epsilon,sigma}, PlotLabel->DRUCKER PRAGER,AxesOrigin->{0,0},Filling -> Axis , FillingStyle -> Directive[Opacity[0.5], Orange]]";
2125  }
2126 
2127  }
2128 
2129 }
2130 
2132 {
2133 // int choice =1;
2134  //TPZPlasticTest::InitializeLOG();
2135  ofstream outfiletxt("mohrcoulomb.txt");
2136 
2137  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
2138 
2139  deltastress.XX() = -13.5;
2140  deltastress.XY() = -0.01;
2141  deltastress.XZ() = 0.;
2142  deltastress.YY() = 0.;
2143  deltastress.YZ() = 0.;
2144  deltastress.ZZ() = -0.01;
2145  stress = deltastress;
2146 
2147  // TPZFNMatrix<6*6> Dep(6,6,0.);
2148  // deltastrain.XX() = -0.0002;
2149  // deltastrain.XY() = 0.;
2150  // deltastrain.XZ() = 0.;
2151  // deltastrain.YY() = -0.00000001;
2152  // deltastrain.YZ() = 0.;
2153  // deltastrain.ZZ() = -0.00000001;
2154  // strain=deltastrain;
2155 
2157  TPZMohrCoulomb Pstep;
2158  Pstep.fYC.SetUp(/*phi=20*/ 20./180. * M_PI);
2159  REAL coesao = 9.2376;
2160  Pstep.fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ coesao, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1000.);
2161  Pstep.fER.SetUp(/*young*/ 20000., /*poisson*/ 0.);
2162 
2163  int length =1;
2164  for(int step=0;step<length;step++)
2165  {
2166  cout << "\nstep "<< step;
2167 
2168  //
2169  // if(step == 30 || step== 30 || step == 60)
2170  // {
2171  // deltastrain *= -1.;
2172  // deltastress *= -1.;
2173  // }
2174 
2175  stress += deltastress;
2176  // strain+=deltastrain;
2177  Pstep.ApplyLoad(stress,strain);
2178  // Pstep.ApplyStrainComputeDep(strain,stress, Dep);
2179  //TPZTensor<REAL> eigenval,dSigma1,dSigma2,dSigma3;
2180  //stress.Eigenvalue(eigenval,dSigma1,dSigma2,dSigma3);
2181  cout<< "\nstress " << stress << endl;
2182  cout<< "\nstrain " << strain << endl;
2183  Pstep.Print(cout);
2184  TPZVec<REAL> phis(1);
2185  Pstep.Phi(strain, phis);
2186  cout << "\nphis " << phis << endl;
2187  //cout<< "\nEigen " << eigenval << endl;
2188 
2189  outfiletxt << fabs(strain.XX()) << " " << fabs(stress.XX()) << "\n";
2190 
2191  }
2192 
2193 
2194 
2195 }
2196 
2197 
2199 {
2200 
2202  ofstream outfile("mohr");
2203 
2204  int choice;
2205  choice =1;
2206  cout << "\nChoose Load Case:";
2207  cout << "\n0 - Strain Step";
2208  cout << "\n1 - Stress Step";
2209  cout << "\n";
2210 
2211  //cin >> choice;
2212  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
2213 
2214  deltastrain.XX() = -0.000005;
2215  deltastrain.YY() = -0.000005;
2216  deltastress.XX()= -10.;
2217  deltastress.YY()= -20.;
2218 
2220  TPZModifiedMohrCoulomb Pstep;
2221  REAL pi = M_PI;
2222  Pstep.fYC.SetUp(/*phi=20*/ 20./180. * pi );
2223  REAL coesao = 9.2376;
2224  Pstep.fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ coesao, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1000.);
2225  Pstep.fER.SetUp(/*young*/ 20000., /*poisson*/ 0.);
2226 
2227 
2228  int length = 1;
2229  for(int step=0;step<length;step++)
2230  {
2231 
2232  if(choice == 0)
2233  {
2234  strain += deltastrain;
2235  Pstep.ApplyStrainComputeSigma(strain,stress);
2236  }
2237 
2238  if(choice == 1)
2239  {
2240  stress += deltastress;
2241  Pstep.ApplyLoad(stress, strain);
2242  }
2243 
2244 #ifdef LOG4CXX
2245  {
2246  std::stringstream sout;
2247  sout << "\nMOHRCOULOMB\n";
2248  sout<<"\n";
2249  sout << "STRAIN" <<std::endl;
2250  sout<<strain<<std::endl;
2251  sout << "STRESS"<<std::endl;
2252  sout<<stress<<std::endl;
2253 
2254  sout<<"STRESS-I1"<<std::endl;
2255  sout<<stress.I1()<<std::endl;
2256  sout<<"STRESS-I2"<<std::endl;
2257  sout<<stress.I2()<<std::endl;
2258  sout<<"STRESS-I3"<<std::endl;
2259  sout<<stress.I3()<<std::endl;
2260 
2261  sout<<"STRESS-J2"<<std::endl;
2262  sout<<stress.J2()<<std::endl;
2263  sout<<"STRESS-J3"<<std::endl;
2264  sout<<stress.J3()<<std::endl;
2265 
2266 
2267  // LOGPZ_INFO(MaterialPoint,sout.str());
2268  }
2269 #endif
2270 
2271  }
2272 }
2273 
2275 {
2276 // int choice =1;
2277  //TPZPlasticTest::InitializeLOG();
2278  ofstream outfiletxt("TPZYCWillamWarnke.txt");
2279 
2280  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
2281  TPZFNMatrix<6*6> Dep(6,6,0.);
2282 
2283  deltastress.XX() = -1.;
2284  deltastress.YY()= 0.;
2285  deltastress.ZZ()= 0.;
2286  deltastress.XY() = 0.;
2287  deltastress.XZ() = 0.;
2288  deltastress.YZ()= 0.;
2289  stress = deltastress;
2290 
2291  // deltastrain.XX() = -0.0001;
2292  // deltastrain.XY() = 0.0000001;
2293  // deltastrain.XZ() = -0.0000001;
2294  // deltastrain.YY() = -0.0000001;
2295  // deltastrain.YZ() = -0.0000001;
2296  // deltastrain.ZZ() = 0.0000002;
2297  // strain=deltastrain;
2298 
2300  TPZWillamWarnke WW;
2301  WW.fYC.SetUp(1.,1.,20.);
2302  REAL coesao = 9.2376;
2303  WW.fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ coesao, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1000.);
2304  WW.fER.SetUp(/*young*/ 20000., /*poisson*/ 0.2);
2305 
2306 
2307  int length = 40;
2308  for(int step=0;step<length;step++)
2309  {
2310  cout << "\nstep "<< step;
2311 
2312 
2313  // if(step == 50 || step==100 || step == 150)
2314  // {
2315  // deltastrain *= -1.;
2316  // deltastress *= -1.;
2317  // }
2318 
2319  stress += deltastress;
2320  // strain+=deltastrain;
2321  WW.ApplyLoad(stress,strain);
2322 #ifdef LOG4CXX
2323  {
2324  std::stringstream sout;
2325  sout << "\n *********************************** STEP LOAD IN WILLANWERNAKE = "<< step << "*********************************** \n";
2326  LOGPZ_INFO(plasticIntegrLogger,sout.str());
2327  }
2328 #endif
2329  // WW.ApplyStrainComputeDep(strain,stress,Dep);
2330  TPZTensor<REAL> eigenval,dSigma1,dSigma2,dSigma3;
2331  stress.Eigenvalue(eigenval,dSigma1,dSigma2,dSigma3);
2332  cout<< "\nstress " << stress << endl;
2333  cout<< "\nEigen " << eigenval << endl;
2334 
2335  outfiletxt << strain.XX() << " " << stress.XX() << "\n";
2336 
2337  }
2338 
2339 
2340 
2341 
2342 
2343 
2344 }
2345 
2346 
2347 
2348 #include <math.h>
2349 inline void RotationMatrix(TPZFMatrix<REAL> &R, REAL thetaRad, int axis)
2350 {
2351  R.Resize(3,3);
2352 
2353 
2354  switch (axis)
2355  {
2356 
2357  case 0://ROTATE ABOUT X
2358 
2359  R.Put(0,0,1.);
2360  R.Put(1,1,cos(thetaRad));R.Put(1,2,sin(thetaRad));
2361  R.Put(2,1,-sin(thetaRad));R.Put(2,2,cos(thetaRad));
2362 
2363  break;
2364 
2365  case 1://ROTATE ABOUT Y
2366 
2367  R.Put(1,1,1.);
2368  R.Put(0,0,cos(thetaRad));R.Put(0,2,sin(thetaRad));
2369  R.Put(2,0,-sin(thetaRad));R.Put(2,2,cos(thetaRad));
2370 
2371  break;
2372 
2373  case 2://ROTATE ABOUT Z
2374 
2375  R.Put(0,0,cos(thetaRad)); R.Put(0, 1,sin(thetaRad));
2376  R.Put(1,0,-sin(thetaRad)); R.Put(1, 1,cos(thetaRad));
2377  R.Put(2,2,1.);
2378 
2379  break;
2380 
2381  case 3:
2382 
2383 
2384  R.Put(0,0,(1./3.)*(1+2*cos(thetaRad)));R.Put( 0,1, (1./3.) * (1-cos(thetaRad) - sqrt(3)*sin(thetaRad)) );R.Put( 0,2, (1./3.) * (1-cos(thetaRad) + sqrt(3)*sin(thetaRad)) );
2385  R.Put( 1,0, (1./3.) * (1-cos(thetaRad) + sqrt(3)*sin(thetaRad)) );R.Put(1,1,(1./3.)*(1+2*cos(thetaRad)));R.Put( 1,2, (1./3.) * (1-cos(thetaRad) - sqrt(3)*sin(thetaRad) ));
2386  R.Put( 2,0, (1./3.) * ( 1-cos(thetaRad) - sqrt(3)*sin(thetaRad) ) );R.Put( 2,1, (1./3.) * ( 1-cos(thetaRad) + sqrt(3)*sin(thetaRad) ) );R.Put(2,2,(1./3.)*(1+2*cos(thetaRad)));
2387 
2388  break;
2389 
2390  default:
2391 
2392  std::cout << " NON SPECIFIED AXIS "<<std::endl;
2393  break;
2394  }
2395 
2396 }
2397 
2398 
2399 template <class T>
2400 inline void MultiDirectionsMaterialPointTest(T & plasticModel, REAL dirMult)
2401 {
2402 
2403  std::ifstream input("../SnubDodecahedron.txt");
2404 
2405  int sizedirs;
2406  input >>sizedirs;
2407  TPZFMatrix<REAL> directions(sizedirs,3,0.);
2408 
2409  TPZTensor<REAL> DiagonalStress,epst,epsp,DeltaDiagonalStress;
2410  int nyield = plasticModel.NYield;
2411 
2412  TPZVec<REAL> funcs(nyield);
2413  TPZFNMatrix<6*6> Dep(6,6,0.);
2414  REAL coordxx,coordyy,coordzz;
2415  for(int i=0;i<sizedirs;i++)
2416  {
2417 
2418  input >> coordxx >> coordyy >> coordzz;
2419  DiagonalStress.XX()=coordxx;
2420  DiagonalStress.YY()=coordyy;
2421  DiagonalStress.ZZ()=coordzz;
2422  DiagonalStress*=dirMult;
2423 
2424  cout << "\n DiagonalStress" << DiagonalStress << endl;
2425  plasticModel.ApplyLoad(DiagonalStress,epst);
2426 
2427  /*
2428 
2429  TPZPlasticState<REAL> stateN = plasticModel.GetState();
2430  TPZPlasticState<REAL> stateN1;
2431  plasticModel.FindPointAtYield(epst,stateN1);
2432 
2433  cout << "\n stateN " << stateN << endl;
2434  cout << "\n stateN1 " << stateN1 << endl;
2435  TPZTensor<REAL> epsTatYield = stateN1.EpsT();
2436 
2437  plasticModel.ApplyStrainComputeDep(epsTatYield,DiagonalStress,Dep);
2438  cout << "\n DiagonalStressAtYield " << DiagonalStress << endl;
2439  */
2440 
2441  bool Plastifica = false;
2442 
2443 #ifdef LOG4CXX
2444  {
2445 
2446  std::stringstream sout;
2447  //sout << " \n Dep = " << Dep << endl;
2448  sout << " \n DIRECTION Number = " << i << "\n " <<endl;
2449  sout << " \n DiagonalStress = " << DiagonalStress << "\n " <<endl;
2450  LOGPZ_INFO(plasticIntegrLogger,sout.str());
2451  }
2452 #endif
2453  int count =0;
2454 
2455  do{
2456 
2457  //epst*=1.1;
2458  DiagonalStress*=1.1;
2459  cout << " \ncount "<< count << endl;
2460  cout << "\n DiagonalStress " << DiagonalStress << endl;
2461  cout << "\n epsT " << epst << endl;
2462  cout << "\n funcs " << funcs << endl;
2463  plasticModel.Phi(epst, funcs);
2464  //plasticModel.ApplyStrainComputeDep(epst,DiagonalStress,Dep);
2465  plasticModel.ApplyLoad(DiagonalStress,epst);
2466 
2467  // Dep.VerifySymmetry();
2468 #ifdef LOG4CXX
2469  {
2470 
2471  std::stringstream sout;
2472 
2473 
2474  sout << " \n\n While loop number = " << count << endl;
2475  sout << " \n\n funcs = " << funcs <<"\n"<< endl;
2476  sout << " \n\n DiagonalStressInsideWhile = " << DiagonalStress <<endl;
2477  sout<< "\nvector"<<count<<" = {" << DiagonalStress.XX()<<","<<DiagonalStress.YY()<<","<<DiagonalStress.ZZ()<< "};" << endl;
2478  //sout << " \n fTFA = " << plasticModel.fTFA.Compute( plasticModel.GetState().VolHardening() ) <<endl;
2479  sout << " \n Alpha() = " << plasticModel.GetState().VolHardening() << endl;
2480  sout << " \n epst = " << plasticModel.GetState().EpsT() <<endl;
2481  sout << " \n epsP = " << plasticModel.GetState().EpsP() <<endl;
2482  //sout << "\n Dep = " << Dep << endl;
2483 
2484  LOGPZ_INFO(plasticIntegrLogger,sout.str());
2485  }
2486 #endif
2487  if(funcs[0] < -1000.)break;
2488 
2489  for(int j = 0;j<nyield;j++)
2490  {
2491  if(funcs[j]>=0.)
2492  {
2493  Plastifica = true;
2494  }
2495 
2496  }
2497 
2498  count++;
2499 
2500  }while(Plastifica==false);
2501 
2502 #ifdef LOG4CXX
2503  {
2504 
2505  std::stringstream sout;
2506  sout << " ######## saindo do while =####### "<< endl;
2507  LOGPZ_INFO(plasticIntegrLogger,sout.str());
2508  }
2509 #endif
2510 
2511  }
2512 
2513 }
2514 
2515 
2516 //Metodo que aplica a tensao e calcula a deformacao de um estado de tensoes um pouco fora da zona elastica.
2517 //Este metodo tambem rotaciona este estado de tensoes de forma a cobrir toda a superficie do plano em questao oriunda das constantes materiais escolhida pelo usuario.
2518 //A variacao na direcao do eixo hidrostatico pode ser modificada alterando as constantes materiais no caso de materiais sensiveis a pressao hidrostatica.
2519 //Neste metodo tambem apos a aplicacao da tensao ligeiramente fora da zona elastica, e obtida a correspondente deformacao total, e chamado o metodo GlobalCheckConv,
2520 //que verificase a matriz de rigidez elastoplastica e realmente a jacobiana dSigma/dEpsT. Neste metodo tambem o modelo plastico importado e submetido a deformacao imposta. A tensao
2521 //e a matriz jacobiana calculadas sao mantidas na memoria. Sucessivas copias deste modelo de plastico sao submetidos a mais deformacao e a verificação de convergencia do método é avaliada.
2522 template <class T>
2524 {
2525 
2526  std::ifstream input("../SnubDodecahedron.txt");
2527  int sizedirs;
2528  input >>sizedirs;
2529  TPZFMatrix<REAL> directions(sizedirs,3,0.);
2530  TPZTensor<REAL> DiagonalStress,epst,epsp,stress;
2531  TPZFNMatrix<6*6> Dep(6,6,0.);
2532  int nyield = mat.NYield;
2533  TPZVec<REAL> funcs(nyield);
2534  REAL pa = 1.;
2535  bool Plastifica = false;
2536  bool MustStop = false;
2537  for(int i=0;i<sizedirs;i++)
2538  {
2539 
2540  REAL coordxx,coordyy,coordzz;
2541  input >> coordxx >> coordyy >> coordzz;
2542 
2543  DiagonalStress.XX()=pa*coordxx;
2544  DiagonalStress.YY()=pa*coordyy;
2545  DiagonalStress.ZZ()=pa*coordzz;
2546  T plasticModelCopy(mat);
2547  plasticModelCopy.ApplyLoad(DiagonalStress,epst);
2548 
2549 
2550  int count =0;
2551  while(( Plastifica == false && MustStop == false) || count < 30 ){
2552 
2553 
2554  plasticModelCopy.ApplyStrainComputeDep(epst,stress,Dep);
2555 
2556 
2557  cout << "\n count " <<count <<endl;
2558  cout << "\n EPST = "<< epst << endl;
2559  cout << "\n EPSP = "<< plasticModelCopy.GetState().EpsP() << endl;
2560  cout << "\n STRESS = "<< stress << endl;
2561 
2562  epst*=1.1;
2563 
2564 
2565 
2566  for(int j = 0;j<nyield;j++)
2567  {
2568  if(funcs[j]>=0.)
2569  {
2570  Plastifica = true;
2571 
2572  }
2573  if(funcs[j] < -1000.)
2574  {
2575  MustStop=true;
2576  }
2577 
2578  }
2579  count++;
2580 
2581 
2582  }
2583 
2584 #ifdef LOG4CXX
2585  {
2586 
2587  std::stringstream sout;
2588  sout << " \n count = "<< count << endl;
2589  sout << " \n PlasticIntegratorCheck = "<< endl;
2590  sout << " \n funcs = " << funcs << endl;
2591  sout << " \n DiagonalStress = " << stress <<endl;
2592  sout << " \n Alpha() = " << plasticModelCopy.GetState().VolHardening() << endl;
2593  sout << " \n epst = " << plasticModelCopy.GetState().EpsT() <<endl;
2594  sout << " \n epsP = " << plasticModelCopy.GetState().EpsP() <<endl;
2595  sout << "\n Dep = " << Dep << endl;
2596  LOGPZ_INFO(plasticIntegrLogger,sout.str());
2597 
2598  }
2599 #endif
2600  GlobalCheckConv(plasticModelCopy, epst, 0.0000001);
2601 
2602  }
2603 
2604 }
2605 
2606 #include "TPZDruckerPrager.h"
2607 
2609 {
2610 
2611  ofstream outfiletxt1("e1dp052NewThermoAModulusCorrected.txt");
2612  ofstream outfiletxt2("e2dp052NewThermoAModulusCorrected.txt");
2613  ofstream outfiletxt3("e3dp052NewThermoAModulusCorrected.txt");
2614  ofstream outfiletxt4("Voldp052NewThermoAModulusCorrected.txt");
2615  TPZTensor<REAL> stress, strain, deltastress, deltastrain;
2616  //
2617  // deltastress.XX() = -0.5;
2618  // deltastress.XY() = -0.001;
2619  // deltastress.XZ() = -0.001;
2620  // deltastress.YY() = -0.001;
2621  // deltastress.YZ() = -0.001;
2622  // deltastress.ZZ() = -0.001;
2623 
2624  deltastress.XX() = -60.;
2625  deltastress.XY() = -0.;
2626  deltastress.XZ() = -0.;
2627  deltastress.YY() = -0.;
2628  deltastress.YZ() = -0.;
2629  deltastress.ZZ() = -3.18;
2630 
2631 
2632 
2633  // deltastrain.XX() = -0.00001;
2634  // deltastrain.XY() = -0.0000001;
2635  // deltastrain.XZ() = -0.0000003;
2636  // deltastrain.YY() = -0.0000015;
2637  // deltastrain.YZ() = -0.00000015;
2638  // deltastrain.ZZ() = -0.0000017;
2639  // strain=deltastrain;
2640 
2641  // TPZSandlerDimaggio SD;
2642  // TPZSandlerDimaggio::McCormicRanchSand(SD);
2643 
2644  TPZLadeKim LK;
2645  //TPZLadeKim::FineSilicaSand(LK);
2646  //TPZLadeKim::DenseSacrRiverSand(LK);
2648 
2650  TPZDruckerPrager DP;
2651  REAL pi = M_PI;
2652  /*innerMCFit = 0*/
2653  /*OuterMCFit = 1*/
2654  DP.fYC.SetUp(/*phi=20*/ 20./180. * pi ,/*MCFit*/0);
2655  REAL coesao = 1351;//PSI//9.2376;
2656  DP.fTFA.SetUp(/*yield- coesao inicial correspondeno a fck igual 32 Mpa */ coesao, /*k Modulo de hardening da coesao equivante 1 Mpa a cada 0.1% de deformacao */1.);
2657  DP.fER.SetUp(/*young*/5318389., /*poisson*/ 0.18);
2658 
2659  // LK.Print(cout);
2660  DP.ApplyLoad(stress, strain);
2661 
2662  // deltastress.XX() = -0.8;
2663  // deltastress.XY() = -0.001;
2664  // deltastress.XZ() = -0.001;
2665  // deltastress.YY() = -0.001;
2666  // deltastress.YZ() = -0.001;
2667  // deltastress.ZZ() = -0.001;
2668  // deltastress.XX() = -1000.;
2669  // deltastress.XY() = -0.001;
2670  // deltastress.XZ() = -0.003;
2671  // deltastress.YY() = -0.15;
2672  // deltastress.YZ() = -0.0015;
2673  // deltastress.ZZ() = -0.17;
2674  stress = deltastress;
2675 
2676  int length =100;
2677  for(int step=0;step<length;step++)
2678  {
2679  /*if(step == 40 || step == 80){
2680 
2681  deltastress*=-1.;
2682 
2683  }*/
2684  cout << "\nstep "<< step;
2685  DP.ApplyLoad(stress,strain);
2686 // REAL pa = stress.I1()/3.;
2687  outfiletxt1 << strain.XX() << " " << fabs(stress.XX()) << "\n";
2688  outfiletxt2 << strain.YY() << " " << fabs(stress.XX()) << "\n";
2689  outfiletxt3 << strain.ZZ() << " " << fabs(stress.XX()) << "\n";
2690  outfiletxt4 << strain.I1() << " " << fabs(stress.I1()) << "\n";
2691  stress += deltastress;
2692  cout << "strain = "<<strain <<"\n";
2693  cout << "sigma = "<< stress <<"\t "<< "I1 = "<< stress.I1() <<"\n";
2694 
2695  }
2696 
2697 
2698 }
2699 
2700 
2702 {
2703 
2704  ofstream outfiletxt1("Mises25.txt");
2705  ofstream outfiletxt2("Mises50.txt");
2706  TPZTensor<REAL> stress, strain, deltastress, deltastrain,stress2,strain2;
2707 
2708  deltastress.XX() = 250.;
2709  deltastress.XY() = 0.0001;
2710  deltastress.XY() = 0.0001;
2711  deltastress.YY() = 0.0001;
2712  deltastress.YZ() = 0.0001;
2713  deltastress.ZZ() = 0.0001;
2714  stress = deltastress;
2715  stress2= deltastress;
2716 
2717 
2718  TPZVonMises VM1;/*CA25*/
2719  VM1.fTFA.SetUp(250.,2500./*(300/0.12)*/);
2720  VM1.fER.SetUp(/*young*/ 210000., /*poisson*/ 0.3);
2721 
2722  TPZVonMises VM2;/*CA50*/
2723  VM2.fTFA.SetUp(500.,7812.5/*(625/0.08)*/);
2724  VM2.fER.SetUp(/*young*/ 210000., /*poisson*/ 0.3);
2725 
2726  int length = 2;
2727  for(int step=0;step<length;step++)
2728  {
2729  cout << "\nstep "<< step;
2730 
2731  if(step==350|| step ==700)
2732  {
2733  deltastrain *= -1.;
2734  deltastress *= -1.;
2735  }
2736 
2737  stress += deltastress;
2738  VM1.ApplyLoad(stress,strain);
2739 
2740  outfiletxt1 << strain.XX() << " " << stress.XX() << "\n";
2741 
2742 
2743  }
2744  /*
2745  int length2 = 2000;
2746  for(int step=0;step<length2;step++)
2747  {
2748  cout << "\nstep "<< step;
2749 
2750  if(step==600|| step ==1200)
2751  {
2752  deltastrain *= -1.;
2753  deltastress *= -1.;
2754  }
2755 
2756  stress2 += deltastress;
2757  VM2.ApplyLoad(stress2,strain2);
2758  outfiletxt2 << strain2.XX() << " " << stress2.XX() << "\n";
2759 
2760 
2761  }
2762  */
2763 }
2764 
2765 
2766 
2768 {
2769  TPZTensor<REAL> epsT, epsTNp1, deltaEpsT;
2770  TPZTensor<REAL> epsP;
2771  REAL alpha;
2772 
2773  epsT.fData[_XX_] = -0.000279312;
2774  epsT.fData[_XY_] = -0.0033376;
2775  epsT.fData[_XZ_] = 3.43561e-10;
2776  epsT.fData[_YY_] = 0.0180517;
2777  epsT.fData[_YZ_] = 3.96378e-10;
2778  epsT.fData[_ZZ_] = 0.0318554;
2779 
2780  epsP.fData[_XX_] = -0.00170072;
2781  epsP.fData[_XY_] = -0.00130612;
2782  epsP.fData[_XZ_] = 1.4072e-10;
2783  epsP.fData[_YY_] = 0.0140323;
2784  epsP.fData[_YZ_] = 1.48982e-10;
2785  epsP.fData[_ZZ_] = 0.028329;
2786 
2787  alpha = 50.4337;
2788 
2789  epsTNp1.fData[_XX_] = -0.00142113;
2790  epsTNp1.fData[_XY_] = -0.00528525;
2791  epsTNp1.fData[_XZ_] = 9.33585e-10;
2792  epsTNp1.fData[_YY_] = 0.0184717;
2793  epsTNp1.fData[_YZ_] = 9.64833e-10;
2794  epsTNp1.fData[_ZZ_] = 0.0318554;
2795 
2796  TPZPlasticState<REAL> state, state2;
2797 
2798  state.m_eps_t = epsT;
2799  state.m_eps_p = epsP;
2800  state.m_hardening = alpha;
2801 
2802  TPZLadeKim LK;
2805  LK.SetState(state);
2806 
2807  LK.SetIntegrTol(0.0001);
2808 
2809  //LK.Print(cout);
2810 
2811  TPZTensor<REAL> sigma;
2812  TPZFNMatrix<6*6> Dep(6,6,0.);
2813 
2814  LK.ApplyStrainComputeDep(epsTNp1, sigma,Dep);
2815  //cout << " \n DEP \n";
2816  //cout << Dep << endl;
2817 
2818  // LK.ApplyStrain(epsTNp1);
2819  // LK.Print(cout);
2820 }
2821 
2822 //Verifica a se a derivada N = dphi/dsigma e consistente
2824 {
2825  typedef TFad<6, REAL> TFAD;
2826 
2827  REAL A;
2828  TFAD A_FAD;
2829 
2830  TPZTensor<REAL> sigma;
2831  sigma.XX() = -0.1;
2832  sigma.XY() = -0.001;
2833  sigma.XZ() = -0.003;
2834  sigma.YY() = -0.15;
2835  sigma.YZ() = -0.0015;
2836  sigma.ZZ() = -0.17;
2837 
2838  TPZTensor<TFAD> sigma_FAD;
2839  sigma.CopyTo(sigma_FAD);
2840 
2841  int i;
2842  for(i = 0; i < 6; i++)sigma_FAD.fData[i].diff(i,6);
2843 
2844  TPZYCSandlerDimaggio YCSD;
2846 
2847  TPZVec<TFAD> res_FAD(2);
2848  TPZVec<TPZTensor<REAL> > NDir(2);
2849 
2850  YCSD.N(sigma, A, NDir,0);
2851  YCSD.Compute(sigma_FAD, A_FAD, res_FAD,0);
2852 
2853  cout << "\nNDir = " << NDir;
2854  cout << "\nres_FAD = " << res_FAD;
2855 
2856 
2857 }
2858 
2860 {
2861  typedef TFad<6, REAL> TFAD;
2862 
2863  //REAL A;
2864  TFAD A_FAD;
2865 
2866  TPZTensor<REAL> sigma, strain;
2867  strain.XX() = -0.1;
2868  strain.XY() = -0.001;
2869  strain.XZ() = -0.003;
2870  strain.YY() = -0.15;
2871  strain.YZ() = -0.0015;
2872  strain.ZZ() = -0.17;
2873 
2874  strain.XX() = -1;
2875  strain.XY() = 0;
2876  strain.XZ() = 0;
2877  strain.YY() = 0.5;
2878  strain.YZ() = 0;
2879  strain.ZZ() = 0.5;
2880 
2881  REAL multipl;
2882 
2883  cin >> multipl;
2884 
2885  strain *= multipl;
2886 
2887  TPZFNMatrix<6*6> Dep(6,6,0.);
2888 
2891  SD.SetIntegrTol(1.e-0);
2892 
2893  SD.ApplyStrainComputeDep(strain,sigma,Dep);
2894 
2895  cout << "\n sigma: " << sigma;
2896 
2897  cout << "\n Dep: " << Dep;
2898 
2899  cout << "\n PlasticState: " << SD.GetState();
2900 
2901  strain *= 2.;
2902 
2903  SD.ApplyStrainComputeDep(strain,sigma,Dep);
2904 
2905  cout << "\n sigma: " << sigma;
2906 
2907  cout << "\n Dep: " << Dep;
2908 
2909  cout << "\n PlasticState: " << SD.GetState();
2910 
2911  strain.XX() -= 0.00001;
2912 
2913  SD.ApplyStrainComputeDep(strain,sigma,Dep);
2914 
2915  cout << "\n sigma: " << sigma;
2916 
2917  cout << "\n Dep: " << Dep;
2918 
2919  cout << "\n PlasticState: " << SD.GetState();
2920 
2921 }
void RotationMatrix(TPZFMatrix< REAL > &R, REAL thetaRad, int axis)
#define _XZ_
Definition: TPZTensor.h:29
virtual void ApplyLoad(const TPZTensor< REAL > &sigma, TPZTensor< REAL > &epsTotal) override
Definition: TPZLadeKim.h:187
Classe que efetua avanco de um passo de plastificacao utilizando o metodo de Newton.
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 LoadState(TPZFMatrix< REAL > &state)
void CopyTo(TPZTensor< T1 > &target) const
Definition: TPZTensor.h:819
void LKLoadingTest()
static void UncDeepSandRes(TPZSandlerDimaggio &material)
static void DruckerPragerTest()
void Add(const TPZTensor< T1 > &tensor, const T2 &constant)
Definition: TPZTensor.h:757
filename
Definition: stats.py:82
YC_t fYC
Object which represents the yield criterium.
T Norm() const
Definition: TPZTensor.h:854
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.
static void UndocumentedTest4()
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
static void UncDeepSandResPSI(TPZSandlerDimaggio &material)
void LKBiaxialTest()
static void DenseSacrRiverSand(TPZLadeKim &material)
Definition: TPZLadeKim.h:383
TPZTensor< T > m_eps_t
Tensors representing the total and plastic strain states.
TPZManVector< T, 6 > fData
Definition: TPZTensor.h:652
virtual TPZPlasticState< REAL > GetState() const =0
void DruckerBiaxialTest()
static void McCormicRanchSandMod2(TPZSandlerDimaggio &material)
T I1() const
Definition: TPZTensor.h:903
virtual void SetState(const TPZPlasticState< REAL > &state)=0
Templated vector implementation.
T m_hardening
Plastic volumetric hardeing variable.
static void UndocumentedTest3()
void LKFineSilicaLoadTest()
T & YY() const
Definition: TPZTensor.h:578
static void McCormicRanchSand(TPZSandlerDimaggio &material)
clarg::argBool h("-h", "help message", false)
clarg::argString input("-if", "input file", "cube1.txt")
static void MohrCoulombTest()
static void StrainTest(T &plasticModel, const char *filename, REAL strainMultiplier=1)
TPZTensor< REAL > gRefTension
static void McCormicRanchSandMod(TPZSandlerDimaggio &material)
void DruckerIsotropicCompression()
static void McCormicRanchSand(TPZYCSandlerDimaggio &material)
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
virtual void Print(std::ostream &out) const override
void SetUp(REAL poisson, REAL M, REAL lambda, REAL a, REAL m, REAL neta1, REAL ksi2, REAL mu, REAL C, REAL p, REAL h, REAL alpha, REAL pa)
Definition: TPZLadeKim.h:85
static void VonMisesTest()
virtual void ApplyLoad(const TPZTensor< REAL > &sigma, TPZTensor< REAL > &epsTotal)=0
static void ModifiedMohrCoulombTest()
sin
Definition: fadfunc.h:63
T & YZ() const
Definition: TPZTensor.h:582
#define LOGPZ_INFO(A, B)
Define log for informations.
Definition: pzlog.h:89
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
Contains TPZMatrixclass which implements full matrix (using column major representation).
void LKKoCompressionLoadTest()
#define _XX_
Definition: TPZTensor.h:27
TPZVonMises gPlasticModel
Definition: tfad.h:64
#define _YZ_
Definition: TPZTensor.h:31
void Residual(TPZFMatrix< REAL > &res, int icase)
#define _XY_
Definition: TPZTensor.h:28
static void LooseSantaMonicaBeachSand(TPZLadeKim &material)
Definition: TPZLadeKim.h:467
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
static void UndocumentedTest2()
virtual int Put(const int64_t row, const int64_t col, const TVar &value)
Put values with bounds checking if DEBUG variable is defined.
Definition: pzmatrix.h:836
string res
Definition: test.py:151
T J2() const
Definition: TPZTensor.h:927
void SetUp(REAL poisson, REAL E, REAL A, REAL B, REAL C, REAL R, REAL D, REAL W)
static void UncDeepSandResMPa(TPZSandlerDimaggio &material)
void SetUp(REAL &cohesion, REAL &phi, REAL &hardening, REAL &young, REAL &poisson)
#define _YY_
Definition: TPZTensor.h:30
void SetUp(REAL young, REAL poisson, REAL fangle, REAL coesion, REAL hardeningModulus, int InnerOuter=0)
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZTensor< T > m_eps_p
void LKIsotropicCompression()
static void FineSilicaSand(TPZLadeKim &material)
Definition: TPZLadeKim.h:411
static void WillamWarnkeTest()
static void StressTest(T &plasticModel, const char *filename, REAL stressMultiplier=1)
virtual int IntegrationSteps() const
static int CreatePlasticModel(T *(&plasticModel), const char *line)
#define _ZZ_
Definition: TPZTensor.h:32
void SandlerDimaggioIsotropicCompression()
static void ReciprocityTest(T &plasticModel, TPZTensor< REAL > strain1)
void Eigenvalue(TPZTensor< T > &eigenval, TPZTensor< T > &dSigma1, TPZTensor< T > &dSigma2, TPZTensor< T > &dSigma3) const
Definition: TPZTensor.h:1744
void Compute(const TPZTensor< T > &sigma, const T &A, TPZVec< T > &res, int checkForcedYield) const
Calculo do criterio de plastificacao.
T I2() const
Definition: TPZTensor.h:908
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
static void DruckerTest()
void MultiDirectionsMaterialPointTest(T &plasticModel, REAL dirMult)
static void PlasticIntegratorCheck(T mat)
const T & VolHardening() const
virtual void ApplyStrainComputeSigma(const TPZTensor< REAL > &epsTotal, TPZTensor< REAL > &sigma, TPZFMatrix< REAL > *tangent=NULL)=0
T & XX() const
Definition: TPZTensor.h:566
void ComputeTangent(TPZFMatrix< REAL > &tangent, TPZVec< REAL > &coefs, int icase)
T I3() const
Definition: TPZTensor.h:918
void N(const TPZTensor< T > &sigma, const T &A, TPZVec< TPZTensor< T > > &Ndir, int checkForcedYield) const
virtual void ApplyStrainComputeDep(const TPZTensor< REAL > &epsTotal, TPZTensor< REAL > &sigma, TPZFMatrix< REAL > &Dep) override
Definition: TPZLadeKim.h:214
virtual void SetUp(const TPZTensor< REAL > &epsTotal) override
Definition: TPZVonMises.h:53
const TPZTensor< T > & EpsP() const
T & XY() const
Definition: TPZTensor.h:570
expr_ expr_ expr_ expr_ expr_ expr_ asin
Definition: tfadfunc.h:80
T & ZZ() const
Definition: TPZTensor.h:586
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
static void GlobalCheckConv(T &plasticModel, TPZTensor< REAL > &strain, REAL maxDeltaStrain=0.01)
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
static void LooseSacrRiverSand(TPZLadeKim &material)
Definition: TPZLadeKim.h:355
clarg::argString m("-m", "input matrix file name (text format)", "matrix.txt")
virtual void ApplyStrainComputeDep(const TPZTensor< REAL > &epsTotal, TPZTensor< REAL > &sigma, TPZFMatrix< REAL > &Dep) override
static void InitializeLOG()
T & XZ() const
Definition: TPZTensor.h:574
static void LoadTest(const char *filename)
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
virtual void SetState(const TPZPlasticState< REAL > &state) override
Definition: TPZLadeKim.h:162
virtual TPZPlasticState< REAL > GetState() const override
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
static void PlainConcrete(TPZLadeKim &material)
Definition: TPZLadeKim.h:299
T J3() const
Definition: TPZTensor.h:946
void LadeKimTriaxialLooseSand()