NeoPZ
TPZSandlerExtended.cpp
Go to the documentation of this file.
1 //
2 // pzsandlerextPV.cpp
3 // PZ
4 //
5 // Created by Diogo Cecilio on 9/3/13.
6 //
7 //
8 
9 #include "TPZSandlerExtended.h"
10 #include "pzlog.h"
11 #include "pzreferredcompel.h"
12 #include "TPZHWTools.h"
13 #include "TPZConvergenceException.h"
14 #include "TPZInconsistentStateException.h"
15 
16 #ifdef LOG4CXX
17 static LoggerPtr logger(Logger::getLogger("plasticity.poroelastoplastic"));
18 #endif
19 
20 #ifdef LOG4CXX
21 static LoggerPtr loggerConvTest(Logger::getLogger("ConvTest"));
22 #endif
23 
24 TPZSandlerExtended::TPZSandlerExtended() : ftol(1e-8), fA(0), fB(0), fC(0), fD(0), fW(0), fK(0), fR(0), fG(0), fPhi(0), fN(0), fPsi(0), fE(0), fnu(0), fkappa_0(0) {
25 }
26 
27 TPZSandlerExtended::TPZSandlerExtended(const TPZSandlerExtended & copy) {
28  ftol = copy.ftol;
29  fA = copy.fA;
30  fB = copy.fB;
31  fC = copy.fC;
32  fD = copy.fD;
33  fW = copy.fW;
34  fK = copy.fK;
35  fR = copy.fR;
36  fG = copy.fG;
37  fPhi = copy.fPhi;
38  fN = copy.fN;
39  fPsi = copy.fPsi;
40  fE = copy.fE;
41  fnu = copy.fnu;
42  fkappa_0 = copy.fkappa_0;
43  fElasticResponse = copy.fElasticResponse;
44 }
45 
46 TPZSandlerExtended::TPZSandlerExtended(STATE A, STATE B, STATE C, STATE D, STATE K, STATE G, STATE W, STATE R, STATE Phi, STATE N, STATE Psi, STATE kappa_0) :
47 fA(A), fB(B), fC(C), fD(D), fW(W), fK(K), fR(R), fG(G), fPhi(Phi), fN(N), fPsi(Psi), fkappa_0(kappa_0) {
48  fE = (9. * fK * fG) / (3. * fK + fG);
49  fnu = ((3. * fK)-(2. * fG)) / (2 * (3. * fK + fG));
50  TPZElasticResponse ER;
51  ER.SetEngineeringData(fE, fnu);
52  fElasticResponse = ER;
53  ftol = 1.e-8;
54 }
55 
56 TPZSandlerExtended::~TPZSandlerExtended() {
57 
58 }
59 
60 template <class T>
61 T TPZSandlerExtended::F(const T x) const {
62  return (fA - fC * exp(x * fB) - fPhi * x);
63 }
64 
65 template <class T>
66 T TPZSandlerExtended::DF(const T x) const {
67  return -(exp(fB * x) * fB * fC) - fPhi;
68 }
69 
70 STATE TPZSandlerExtended::GetF(STATE x) const {
71  return F(x);
72 }
73 
74 template<class T>
75 T TPZSandlerExtended::X(const T L) const {
76  return (L - fR * F(L));
77 }
78 
79 STATE TPZSandlerExtended::GetX(STATE k) {
80  return X(k);
81 }
82 
83 void TPZSandlerExtended::SetUp(STATE A, STATE B, STATE C, STATE D, STATE K, STATE G, STATE W, STATE R, STATE Phi, STATE N, STATE Psi) {
84  fA = A;
85  fB = B;
86  fC = C;
87  fD = D;
88  fK = K;
89  fG = G;
90  fW = W;
91  fR = R;
92  fPhi = Phi;
93  fN = N;
94  fPsi = Psi;
95  fE = (9. * fK * fG) / (3. * fK + fG);
96  fnu = ((3. * fK)-(2. * fG)) / (2 * (3. * fK + fG));
97  TPZElasticResponse ER;
98  ER.SetEngineeringData(fE, fnu);
99  fElasticResponse = ER;
100 
101 }
102 
103 void TPZSandlerExtended::SetInitialDamage(STATE kappa_0) {
104  fkappa_0 = kappa_0;
105 }
106 
107 void TPZSandlerExtended::SetElasticResponse(const TPZElasticResponse &ER) {
108  fElasticResponse = ER;
109  fE = ER.E();
110  fnu = ER.Poisson();
111  fK = ER.K();
112  fG = ER.G();
113 }
114 
115 TPZElasticResponse TPZSandlerExtended::GetElasticResponse() const {
116  return fElasticResponse;
117 }
118 
119 void TPZSandlerExtended::Read(TPZStream& buf, void* context) { //ok
120  buf.Read(&ftol);
121  buf.Read(&fA);
122  buf.Read(&fB);
123  buf.Read(&fC);
124  buf.Read(&fD);
125  buf.Read(&fW);
126  buf.Read(&fK);
127  buf.Read(&fR);
128  buf.Read(&fG);
129  buf.Read(&fPhi);
130  buf.Read(&fN);
131  buf.Read(&fPsi);
132  buf.Read(&fE);
133  buf.Read(&fnu);
134  buf.Read(&fkappa_0);
135  fElasticResponse.Read(buf, context);
136 }
137 
138 void TPZSandlerExtended::Write(TPZStream& buf, int withclassid) const { //ok
139  buf.Write(&ftol);
140  buf.Write(&fA);
141  buf.Write(&fB);
142  buf.Write(&fC);
143  buf.Write(&fD);
144  buf.Write(&fW);
145  buf.Write(&fK);
146  buf.Write(&fR);
147  buf.Write(&fG);
148  buf.Write(&fPhi);
149  buf.Write(&fN);
150  buf.Write(&fPsi);
151  buf.Write(&fE);
152  buf.Write(&fnu);
153  buf.Write(&fkappa_0);
154  fElasticResponse.Write(buf, withclassid);
155 }
156 
157 TPZElasticResponse TPZSandlerExtended::GetElasticResponse() {
158  return fElasticResponse;
159 }
160 
161 STATE TPZSandlerExtended::GetR() {
162  return fR;
163 }
164 
165 template<class T>
166 T TPZSandlerExtended::EpsEqX(T X) const {
167  STATE CPer = CPerturbation();
168  return (fW * (exp(fD * X) - 1 + CPer * X));
169 }
170 
171 template<class T>
172 T TPZSandlerExtended::EpsEqk(const T k) const {
173  return EpsEqX(X(k));
174 }
175 
176 void TPZSandlerExtended::Firstk(STATE &epsp, STATE &k) const {
177  STATE f, df, kn1, kn, resnorm, diff;
178  int counter = 1;
179  resnorm = 1;
180  kn = epsp; //chute inicial
181  kn1 = kn;
182  while (resnorm > ftol && counter < 30) {
183 
184  f = EpsEqk(kn) - epsp;
185  df = fD * exp(fD * (kn - (fA - fC * exp(fB * kn)) * fR))*(1 + fB * fC * exp(fB * kn) * fR) * fW;
186  //df=fD*exp(fD*(kn - fR*(fA - fC*exp(fB*kn) - kn*fPhi)))*fW*(1 - fR*(-(fB*fC*exp(fB*kn)) - fPhi));
187  kn1 = kn - f / df;
188  diff = kn1 - kn;
189  resnorm = sqrt(diff * diff);
190  kn = kn1;
191  counter++;
192 
193  }
194  k = kn1;
195 }
196 
198 STATE TPZSandlerExtended::NormalToF1(STATE I1, STATE I1_ref) const {
199 
200 #ifdef PZDEBUG
201  if (I1 < I1_ref) { // normal function is constructed for I1 >= I1_ref
202  throw TPZInconsistentStateException("TPZSandlerExtended::NormalToF1: I1 < I1_ref.");
203  }
204 #endif
205 
206  STATE normal_f1 = (exp(fB*I1_ref)*fA*fB*fC - exp(2*fB*I1_ref)*fB*pow(fC,2) + I1 - I1_ref)/(exp(fB*I1_ref)*fB*fC);
207  return normal_f1;
208 }
209 
210 REAL TPZSandlerExtended::InitialDamage(const TPZVec<REAL> &stress_pv) const {
211 
212  TPZManVector<REAL,2> f(2);
213  YieldFunction(stress_pv, 0.0, f);
214 
215  bool Is_valid_stress_Q = fabs(f[0]) < ftol || f[0] < 0.0;
216 
217  if (Is_valid_stress_Q) {
218 
219  int n_iter = 30; // @TODO:: Variable to GUI.
220  REAL I1 = (stress_pv[0])+(stress_pv[1])+(stress_pv[2]);
221  REAL res, jac, dk, k;
222 
223  k = 0.0; // initial guess
224  bool stop_criterion_Q = false;
225  int i;
226  for (i = 0; i < n_iter; i++) {
227  res = I1 - X(k);
228  stop_criterion_Q = fabs(res) < ftol;
229  if (stop_criterion_Q) {
230  break;
231  }
232  jac = - 1.0 - fB * fC * fR * exp(fB*k);
233  dk = - res /jac;
234  k+=dk;
235  }
236 
237  if (!stop_criterion_Q) {
238  throw TPZConvergenceException(ftol, n_iter, res, i, "TPZSandlerExtended::InitialDamage:: Newton process did not converge in hydrostatic direction.");
239  }
240 
241  stop_criterion_Q = false;
242  REAL J2 = (1.0/3.0) * (stress_pv[0]*stress_pv[0] + stress_pv[1]*stress_pv[1] + stress_pv[2]*stress_pv[2] - stress_pv[1]*stress_pv[2] - stress_pv[0]*stress_pv[2] - stress_pv[0]*stress_pv[1]);
243 
244  bool pure_hydrostatic_state_Q = IsZero(J2) || J2 < 0.0;
245  if (!pure_hydrostatic_state_Q) {
246  k = X(k) + k; // guess from the outer part of the cap
247 
248  for (int i = 0; i < n_iter; i++) {
249 
250  res = -1 + pow(I1,2)/(pow(fA - exp(fB*k)*fC,2)*pow(fR,2)) +
251  J2/pow(fA - exp(fB*k)*fC,2) -
252  (2*I1*k)/(pow(fA - exp(fB*k)*fC,2)*pow(fR,2)) +
253  pow(k,2)/(pow(fA - exp(fB*k)*fC,2)*pow(fR,2));
254 
255  stop_criterion_Q = fabs(res) < ftol;
256  if (stop_criterion_Q) {
257  break;
258  }
259  jac = (2*(fA*(-I1 + k) + exp(fB*k)*fC*
260  (I1 + fB*pow(I1,2) + fB*pow(fR,2)*J2 - k - 2*fB*I1*k + fB*pow(k,2))))/
261  (pow(fA - exp(fB*k)*fC,3)*pow(fR,2));
262  dk = - res /jac;
263  k+=dk;
264  }
265 
266  if (!stop_criterion_Q) {
267  throw TPZConvergenceException(ftol, n_iter, res, i, "TPZSandlerExtended::InitialDamage:: Newton process did not converge in deviatoric direction.");
268  }
269 
270  }
271 
272  YieldFunction(stress_pv, k, f);
273  bool Is_valid_stress_on_cap_Q = fabs(f[1]) < ftol || f[1] < 0.0;
274 
275  if (!Is_valid_stress_on_cap_Q) {
276  throw TPZInconsistentStateException("TPZSandlerExtended::InitialDamage: Invalid stress state over cap.");
277  }
278  return k;
279 
280  }
281  else{
282  throw TPZInconsistentStateException("TPZSandlerExtended::InitialDamage: Invalid stress state over failure surface.");
283  }
284 
285  return -1;
286 
287 }
288 
289 template<class T>
290 T TPZSandlerExtended::ResLF2(const TPZVec<T> &trial_stress, T theta, T beta, T k, STATE kprev) const {
291  T trial_I1 = (trial_stress[0])+(trial_stress[1])+(trial_stress[2]);
292  T I1 = fR * F(k) * cos(theta) + k;
293  T delepsp = EpsEqk(k) - EpsEqk(kprev);
294  return (3. * fK * delepsp - (trial_I1 - I1));
295 }
296 
297 template<class T>
298 T TPZSandlerExtended::ResLF2IJ(const TPZVec<T> &sigtrialIJ, T theta, T k, STATE kprev) const {
299 
300  T I1tr = sigtrialIJ[0];
301  T I1 = fR * F(k) * cos(theta) + k;
302  T delepsp = EpsEqk(k) - EpsEqk(kprev);
303  return (3. * fK * delepsp - (I1tr - I1));
304 }
305 
307 
308 STATE TPZSandlerExtended::ResLF1(const TPZVec<STATE> &sigtrial, const TPZVec<STATE> &sigproj, const STATE k, const STATE kprev) const {
309  STATE I1trial = (sigtrial[0] + sigtrial[1] + sigtrial[2]);
310  STATE I1proj = (sigproj[0] + sigproj[1] + sigproj[2]);
311  STATE delepsp = EpsEqk(k) - EpsEqk(kprev);
312  return (3. * fK * delepsp - (I1trial - I1proj));
313 
314 }
315 
317 // the derivative are given in terms of k
318 
319 STATE TPZSandlerExtended::DResLF1(const TPZVec<STATE> &sigtrial, const TPZVec<STATE> &sigproj, const STATE k, const STATE kprev) const {
320  STATE expfBk = exp(fB * k);
321  STATE dreskk = 3. * exp(fD * (-((fA - expfBk * fC) * fR) + k)) * fD * fK *
322  (1. + expfBk * fB * fC * fR) * fW;
323  return dreskk;
324 }
325 
326 
328 
329 void TPZSandlerExtended::DResLF2(const TPZVec<STATE> &pt, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &dresl) const {
330  STATE expfBk = exp(fB * k);
331  STATE sintheta = sin(theta);
332  STATE costheta = cos(theta);
333 
334  STATE dreskk = 1. + costheta * (-(expfBk * fB * fC) - fPhi) * fR +
335  3. * exp(fD * (-((fA - expfBk * fC) * fR) + k)) * fD * fK *
336  (1. + expfBk * fB * fC * fR) * fW;
337 
338 
339  STATE dresktheta = -fR * (fA - fC * expfBk - k * fPhi) * sintheta;
340  dresl[0] = dresktheta;
341  dresl[1] = 0;
342  dresl[2] = dreskk;
343 
344 
345 }
346 
347 void TPZSandlerExtended::F1Cyl(STATE xi, STATE beta, TPZVec<STATE> &f1cyl) const {
348  STATE sqrt2 = M_SQRT2;
349  STATE sqrt3 = sqrt(3.);
350  STATE gamma = 0.5 * (1 + (1 - sin(3 * beta)) / fPsi + sin(3 * beta));
351  STATE I1 = xi*sqrt3;
352  STATE F1 = F(I1);
353  STATE sqrtj2 = F1 / gamma;
354  STATE rho = sqrt2*sqrtj2;
355  f1cyl[0] = xi;
356  f1cyl[1] = rho;
357  f1cyl[2] = beta;
358 
359 }
360 
361 void TPZSandlerExtended::SurfaceParamF1(TPZVec<STATE> &sigproj, STATE &xi, STATE &beta) const {
362  TPZManVector<STATE> sigHWCyl(3);
363  TPZHWTools::FromPrincipalToHWCyl(sigproj, sigHWCyl);
364  xi = sigHWCyl[0];
365  beta = sigHWCyl[2];
366 #ifdef PZDEBUG
367  STATE dist = DistF1(sigproj, xi, beta);
368  if (fabs(dist) > ftol) {
369  DebugStop();
370  }
371 #endif
372 }
373 
374 void TPZSandlerExtended::F2Cyl(STATE theta, STATE beta, STATE k, TPZVec<STATE> &f2cyl) const {
375  const STATE M_SQRT3 = sqrt(3.);
376  const STATE gamma = 0.5 * (1.0 + sin(3.0 * beta) + (1.0 - sin(3.0 * beta)) / fPsi);
377  const STATE Fk = F(k);
378  const STATE var = fR * Fk * cos(theta);
379  const STATE sig1_star = (k - var) / M_SQRT3; // The definition for theta was corrected.
380  const STATE sqrtj2 = Fk * sin(theta) / gamma;
381  const STATE rho = M_SQRT2*sqrtj2;
382  const STATE xi = sig1_star;
383  f2cyl[0] = xi;
384  f2cyl[1] = rho;
385  f2cyl[2] = beta;
386 
387 }
388 
389 void TPZSandlerExtended::SurfaceParamF2(const TPZVec<STATE> &sigproj, const STATE k, STATE &theta, STATE &beta) const {
390  TPZManVector<STATE> sigHWCyl(3);
391  TPZHWTools::FromPrincipalToHWCyl(sigproj, sigHWCyl);
392  // STATE xi,rho;
393  // xi=sigHWCyl[0];
394  // rho=sigHWCyl[1];
395  STATE I1 = sigHWCyl[0] * sqrt(3.);
396  beta = sigHWCyl[2];
397  STATE Fk = F(k);
398  STATE gamma = 0.5 * (1 + (1 - sin(3 * beta)) / fPsi + sin(3 * beta));
399  STATE costheta = (I1 - k) / (fR * Fk);
400  STATE sqrtj2 = sigHWCyl[1] / sqrt(2.);
401  STATE sintheta = sqrtj2 * gamma / (Fk - fN);
402  theta = atan2(sintheta, costheta);
403  //theta = acos(costheta);
404  //STATE theta2 = atan((rho*sin(beta))/xi);
405 #ifdef PZDEBUG
406  STATE err = 1. - sintheta * sintheta - costheta*costheta;
407  STATE dist = DistF2(sigproj, theta, beta, k);
408  if (fabs(dist) > ftol || err > ftol) {
409  DebugStop();
410  }
411 #endif
412 }
413 
414 STATE TPZSandlerExtended::DistF1(const TPZVec<STATE> &pt, STATE xi, STATE beta) const {
415  TPZManVector<STATE, 3> cyl(3);
416  F1Cyl(xi, beta, cyl);
417  TPZManVector<STATE, 3> cart(3);
418  TPZHWTools::FromHWCylToHWCart(cyl, cart);
419  TPZManVector<STATE, 3> carttrial(3);
420  TPZHWTools::FromPrincipalToHWCart(pt, carttrial);
421  return ((1. / (3. * fK))*(carttrial[0] - cart[0])*(carttrial[0] - cart[0]))
422  +(1. / (2. * fG))*((carttrial[1] - cart[1])*(carttrial[1] - cart[1])+(carttrial[2] - cart[2])*(carttrial[2] - cart[2]));
423 
424 }
425 
426 STATE TPZSandlerExtended::DistF2(const TPZVec<STATE> &pt, STATE theta, STATE beta, STATE k) const {
427  TPZManVector<STATE, 3> cart_trial(3); // cyl and cart trial are the same variable, it is renamed as cart_trial
428  TPZManVector<STATE, 3> cart(3);
429  F2Cyl(theta, beta, k, cart_trial);
430  TPZHWTools::FromHWCylToHWCart(cart_trial, cart);
431  TPZHWTools::FromPrincipalToHWCart(pt, cart_trial);
432  return ((1. / (3. * fK))*(cart_trial[0] - cart[0])*(cart_trial[0] - cart[0]))
433  +(1. / (2. * fG))*((cart_trial[1] - cart[1])*(cart_trial[1] - cart[1])+(cart_trial[2] - cart[2])*(cart_trial[2] - cart[2]));
434 
435 }
436 
437 STATE TPZSandlerExtended::DistF2IJ(const TPZVec<STATE> &sigtrialIJ, STATE theta, STATE k) const {
438  STATE I1 = sigtrialIJ[0];
439  STATE sqJ2 = sigtrialIJ[1];
440  STATE Fk;
441  Fk = F(k);
442  STATE y = (sqJ2 - Fk * sin(theta));
443  STATE x = 1. / (3 * fK)*(I1 - (k + Fk * fR * cos(theta)));
444  STATE res = x * x / (9. * fK) + y * y / (fG);
445  return res;
446 
447 }
448 
449 template<class T>
450 void TPZSandlerExtended::FromThetaKToSigIJ(const T &theta, const T &K, TPZVec<T> &sigIJ) const {
451  T Fk = F(K);
452  sigIJ[0] = K + Fk * fR * cos(theta);
453  sigIJ[1] = Fk * sin(theta);
454 }
455 
459 template<class T>
460 void TPZSandlerExtended::DDistF2IJ(TPZVec<T> &sigtrialIJ, T theta, T L, STATE LPrev, TPZVec<T> &ddistf2) const {
461  T I1 = sigtrialIJ[0];
462  T sqJ2 = sigtrialIJ[1];
463  T Fk;
464  Fk = F(L);
465  T y = (sqJ2 - Fk * sin(theta));
466  T x = (I1 - (L + Fk * fR * cos(theta)));
467  ddistf2[0] = T(2.) * x * Fk * fR * sin(theta) / T(9. * fK) - T(2.) * y * Fk * cos(theta);
468  ddistf2[1] = ResLF2IJ(sigtrialIJ, theta, L, LPrev);
469 #ifdef LOG4CXX
470  if (logger->isDebugEnabled()) {
471  std::stringstream sout;
472  sout << "x = " << x << " y = " << y << " theta = " << theta << " res = " << ddistf2[0];
473  LOGPZ_DEBUG(logger, sout.str())
474  }
475 #endif
476 }
477 
478 void TPZSandlerExtended::DDistFunc1(const TPZVec<STATE> &pt, STATE xi, STATE beta, TPZFMatrix<STATE> &ddistf1) const {
479  STATE sigstar1, sigstar2, sigstar3, DFf, Ff, I1, sb, cb, DGamma;
480  TPZManVector<STATE,3> ptcart(3);
481  sb = sin(beta);
482  cb = cos(beta);
483  STATE sin3b = sin(3 * beta);
484  STATE cos3b = cos(3 * beta);
485  TPZHWTools::FromPrincipalToHWCart(pt, ptcart);
486  sigstar1 = ptcart[0];
487  sigstar2 = ptcart[1];
488  sigstar3 = ptcart[2];
489  I1 = xi * sqrt(3);
490  Ff = F(I1);
491  const REAL gamma = (1 + sin3b + (1 - sin3b) / fPsi) / 2.;
492  const REAL gamma2 = gamma*gamma;
493  const REAL gamma3 = gamma*gamma2;
494  DFf = DF(I1);
495  DGamma = 0.5*(3 * cos3b - (3 * cos3b) / fPsi);
496  const REAL sqrt2 = M_SQRT2;
497  const REAL sqrt3 = sqrt(3.);
498 
499  ddistf1.Resize(2, 1);
500  ddistf1(0, 0) = (6 * sqrt3 * DFf * Ff * fK - 3 * sqrt3 * sqrt2 * DFf * fK * gamma * (cb * sigstar2 + sb * sigstar3) + 2 * fG * gamma2 * (-sigstar1 + xi)) / (3. * fG * fK * gamma2);
501  ddistf1(1, 0) = -((Ff * sqrt2 * (gamma2 * (-(sb * sigstar2) + cb * sigstar3) - DGamma * gamma * (cb * sigstar2 + sb * sigstar3) + DGamma * Ff * sqrt2)) / (fG * gamma3));
502 
503 }
504 
505 void TPZSandlerExtended::D2DistFunc1(const TPZVec<STATE> &pt, STATE xi, STATE beta, TPZFMatrix<STATE> &tangentf1) const {
506  STATE sig2, sig3, DFf, Gamma, Ff, I1, sb, cb, D2Ff, DGamma, D2Gamma, Gamma2, Gamma3, Sqrt2, Sqrt3;
507  TPZVec<STATE> ptcart(3);
508  sb = sin(beta);
509  cb = cos(beta);
510  STATE sin3b = sin(3 * beta);
511  STATE cos3b = cos(3 * beta);
512  TPZHWTools::FromPrincipalToHWCart(pt, ptcart);
513  // STATE sig1 = ptcart[0];
514  sig2 = ptcart[1];
515  sig3 = ptcart[2];
516  I1 = xi * sqrt(3);
517  Ff = F(I1);
518  Gamma = (1. + sin3b + (1. - sin3b) / fPsi) / 2.;
519  DFf = -DF(I1);
520  D2Ff = -(exp(fB * I1) * pow(fB, 2.) * fC);
521  DGamma = (3. * cos3b - (3. * cos3b) / fPsi) / 2.;
522  D2Gamma = (-9. * sin(3. * beta) + (9. * sin(3. * beta)) / fPsi) / 2.;
523  Gamma2 = Gamma*Gamma;
524  Gamma3 = Gamma*Gamma2;
525  // STATE D2Gamma2=D2Gamma*D2Gamma;
526  Sqrt2 = sqrt(2);
527  Sqrt3 = sqrt(3);
528 
529  tangentf1.Resize(2, 2);
530 
531 
532  tangentf1(0, 0) = (18 * (DFf * DFf + D2Ff * Ff) * fK + 2 * fG * Gamma2 -
533  9 * D2Ff * fK * Gamma * (cb * sig2 + sb * sig3) * Sqrt2) / (3. * fG * fK * Gamma2);
534 
535  tangentf1(0, 1) = -((Sqrt3 * Sqrt2 * DFf * (Gamma2 * (-(sb * sig2) + cb * sig3) - DGamma * Gamma * (cb * sig2 + sb * sig3) +
536  2 * DGamma * Ff * Sqrt2)) / (fG * Gamma3));
537 
538  tangentf1(1, 1) = -((Ff * (-6 * DGamma * DGamma * Ff - Gamma3 * (cb * sig2 + sb * sig3) * Sqrt2 -
539  Gamma2 * (2 * DGamma * (-(sb * sig2) + cb * sig3) + D2Gamma * (cb * sig2 + sb * sig3)) *
540  Sqrt2 + 2 * Gamma * (D2Gamma * Ff + DGamma * DGamma * (cb * sig2 + sb * sig3) * Sqrt2))) /
541  (fG * Gamma2 * Gamma2));
542 
543  tangentf1(1, 0) = tangentf1(0, 1);
544 
545 }
546 
547 
549 template<class T>
550 void TPZSandlerExtended::Res1(const TPZVec<T> &trial_stress, T i1, T beta, T k, T kprev, TPZVec<T> & residue_1) const{
551 
552  // In this implementation the definition for theta is given by the angle formed from -I1 axis to sqrt(J2) axis with origin on the damage variable kappa.
553 
554  residue_1.Resize(3, 1);
555  STATE CX0 = X_0();
556  STATE CPer = CPerturbation();
557  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
558  STATE CG = fE/(2.0*(1.0 + fnu));
559 
560  TPZManVector<REAL,3> rhw_sigma(3);
561  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
562 
563  residue_1[0] = (2*(i1 - sqrt(3)*rhw_sigma[0]))/(9.*CK) + (2*sqrt(2)*exp(fB*i1)*fB*fC*fPsi*cos(beta)*
564  (-2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*cos(beta) + rhw_sigma[1]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))
565  )/(CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2)) +
566  (2*sqrt(2)*exp(fB*i1)*fB*fC*fPsi*sin(beta)*
567  (-2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*sin(beta) + rhw_sigma[2]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))
568  )/(CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
569 
570 
571  residue_1[1] = (2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*(cos(beta)*(1 + fPsi + 8*(-1 + fPsi)*pow(sin(beta),3))*
572  (2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*sin(beta) -
573  rhw_sigma[2]*(1 + fPsi + (-1 + fPsi)*sin(3*beta))) +
574  (2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*cos(beta) -
575  rhw_sigma[1]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))*
576  (-3*(-1 + fPsi)*cos(beta)*cos(3*beta) - sin(beta)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))))/
577  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3));
578 
579  residue_1[2] = -i1 + 3*CK*fW*(-exp(-(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k))) +
580  exp(-(fD*(CX0 + fA*fR - exp(fB*kprev)*fC*fR - kprev))) - CPer*exp(fB*k)*fC*fR +
581  CPer*exp(fB*kprev)*fC*fR + CPer*(-k + kprev)) + sqrt(3)*rhw_sigma[0];
582 
583 }
584 
585 template<class T>
586 void TPZSandlerExtended::Res2(const TPZVec<T> &trial_stress, T theta, T beta, T k, T kprev, TPZVec<T> &residue_2) const {
587 
588  // In this implementation the definition for theta is given by the angle formed from -I1 axis to sqrt(J2) axis with origin on the damage variable kappa.
589 
590  residue_2.Resize(3, 1);
591  STATE CX0 = X_0();
592  STATE CPer = CPerturbation();
593  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
594  STATE CG = fE/(2.0*(1.0 + fnu));
595 
596  TPZManVector<REAL,3> rhw_sigma(3);
597  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
598 
599 
600 
601  residue_2[0] = (2.0*(fA - exp(fB*k)*fC)*(CG*fR*(k - sqrt(3.0)*rhw_sigma[0] - (fA - exp(fB*k)*fC)*fR*cos(theta))*
602  pow(1.0 + fPsi + (-1.0 + fPsi)*sin(3.0*beta),2.0)*sin(theta) -
603  9.0*sqrt(2.0)*CK*fPsi*cos(beta)*cos(theta)*
604  ((1.0 + fPsi)*rhw_sigma[1] + (-1.0 + fPsi)*rhw_sigma[1]*sin(3.0*beta) -
605  2.0*sqrt(2.0)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*sin(theta)) -
606  9.0*sqrt(2.0)*CK*fPsi*cos(theta)*sin(beta)*
607  ((1 + fPsi)*rhw_sigma[2] + (-1.0 + fPsi)*rhw_sigma[2]*sin(3.0*beta) -
608  2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*sin(beta)*sin(theta))))/
609  (9.0*CG*CK*pow(1.0 + fPsi + (-1.0 + fPsi)*sin(3.0*beta),2.0));
610 
611  residue_2[1] = (-2.0*sqrt(2.0)*(fA - exp(fB*k)*fC)*fPsi*sin(theta)*
612  ((-3.0*(-1.0 + fPsi)*cos(beta)*cos(3.0*beta) - sin(beta)*(1.0 + fPsi + (-1.0 + fPsi)*sin(3.0*beta)))*
613  (rhw_sigma[1]*(1.0 + fPsi + (-1.0 + fPsi)*sin(3.0*beta)) -
614  2.0*sqrt(2.0)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*sin(theta)) +
615  cos(beta)*(1.0 + fPsi + 8.0*(-1.0 + fPsi)*pow(sin(beta),3.0))*
616  (rhw_sigma[2]*(1.0 + fPsi + (-1.0 + fPsi)*sin(3.0*beta)) -
617  2.0*sqrt(2.0)*(fA - exp(fB*k)*fC)*fPsi*sin(beta)*sin(theta))))/
618  (CG*pow(1.0 + fPsi + (-1.0 + fPsi)*sin(3.0*beta),3.0));
619 
620 
621  residue_2[2] = -k + 3.0*CK*fW*(-exp(-(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k))) +
622  exp(-(fD*(CX0 + fA*fR - exp(fB*kprev)*fC*fR - kprev))) - CPer*exp(fB*k)*fC*fR +
623  CPer*exp(fB*kprev)*fC*fR - CPer*k + CPer*kprev) + sqrt(3.0)*rhw_sigma[0] +
624  (fA - exp(fB*k)*fC)*fR*cos(theta);
625 
626 }
627 
628 template<class T>
629 void TPZSandlerExtended::Res2CoVertex(const TPZVec<T> &trial_stress, T beta, T k, T kprev, TPZVec<T> & residue_covertex) const{
630 
631  // In this implementation the definition for theta is given by the angle formed from -I1 axis to sqrt(J2) axis with origin on the damage variable kappa.
632 
633  residue_covertex.Resize(2, 1);
634  STATE CX0 = X_0();
635  STATE CPer = CPerturbation();
636  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
637  STATE CG = fE/(2.0*(1.0 + fnu));
638  STATE theta = M_PI_2;
639 
640  TPZManVector<REAL,3> rhw_sigma(3);
641  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
642 
643  residue_covertex[0] = (2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*(cos(beta)*(1 + fPsi + 8*(-1 + fPsi)*pow(sin(beta),3))*
644  (2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*sin(beta) - rhw_sigma[2]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))
645  + (2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(beta) -
646  rhw_sigma[1]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))*
647  (-3*(-1 + fPsi)*cos(beta)*cos(3*beta) - sin(beta)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))))/
648  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3));
649 
650 
651  residue_covertex[1] = -k + 3*CK*fW*(-exp(-(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k))) +
652  exp(-(fD*(CX0 + fA*fR - exp(fB*kprev)*fC*fR - kprev))) - CPer*exp(fB*k)*fC*fR +
653  CPer*exp(fB*kprev)*fC*fR + CPer*(-k + kprev)) + sqrt(3)*rhw_sigma[0];
654 
655 }
656 
657 template<class T>
658 void TPZSandlerExtended::Res2Vertex(const TPZVec<T> &trial_stress, T k, T kprev, T & residue_vertex) const{
659 
660  STATE CX0 = X_0();
661  STATE CPer = CPerturbation();
662  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
663  STATE theta = 0.0;
664 
665  TPZManVector<REAL,3> rhw_sigma(3);
666  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
667 
668  residue_vertex = -k + 3*CK*fW*(-exp(-(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k))) +
669  exp(-(fD*(CX0 + fA*fR - exp(fB*kprev)*fC*fR - kprev))) - CPer*exp(fB*k)*fC*fR +
670  CPer*exp(fB*kprev)*fC*fR - CPer*k + CPer*kprev) + sqrt(3)*rhw_sigma[0] + (fA - exp(fB*k)*fC)*fR*cos(theta);
671 }
672 
674 void TPZSandlerExtended::Jacobianf1(const TPZVec<STATE> &trial_stress, STATE i1, STATE beta, STATE k, TPZFMatrix<STATE> &jacobianf1)const{
675 
676  jacobianf1.Resize(3, 3);
677  STATE CX0 = X_0();
678  STATE CPer = CPerturbation();
679  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
680  STATE CG = fE/(2.0*(1.0 + fnu));
681 
682  TPZManVector<REAL,3> rhw_sigma(3);
683  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
684 
685 
686  jacobianf1(0,0) = (2*(36*CK*exp(fB*i1)*pow(fB,2)*fC*(-fA + 2*exp(fB*i1)*fC)*pow(fPsi,2)*pow(cos(beta),2) +
687  36*CK*exp(fB*i1)*pow(fB,2)*fC*(-fA + 2*exp(fB*i1)*fC)*pow(fPsi,2)*
688  pow(sin(beta),2) + 9*sqrt(2)*CK*exp(fB*i1)*pow(fB,2)*fC*fPsi*rhw_sigma[1]*cos(beta)*
689  (1 + fPsi + (-1 + fPsi)*sin(3*beta)) +
690  9*sqrt(2)*CK*exp(fB*i1)*pow(fB,2)*fC*fPsi*rhw_sigma[2]*sin(beta)*
691  (1 + fPsi + (-1 + fPsi)*sin(3*beta)) + CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2)))/
692  (9.*CG*CK*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
693 
694 
695  jacobianf1(0,1) = -((sqrt(2)*exp(fB*i1)*fB*fC*fPsi*(-((3 + 2*fPsi + 3*pow(fPsi,2))*rhw_sigma[2]*cos(beta)) +
696  5*(-1 + pow(fPsi,2))*rhw_sigma[1]*cos(2*beta) + 24*sqrt(2)*fA*fPsi*cos(3*beta) -
697  24*sqrt(2)*exp(fB*i1)*fC*fPsi*cos(3*beta) - 24*sqrt(2)*fA*pow(fPsi,2)*cos(3*beta) +
698  24*sqrt(2)*exp(fB*i1)*fC*pow(fPsi,2)*cos(3*beta) - rhw_sigma[1]*cos(4*beta) +
699  pow(fPsi,2)*rhw_sigma[1]*cos(4*beta) + 2*rhw_sigma[2]*cos(5*beta) - 4*fPsi*rhw_sigma[2]*cos(5*beta) +
700  2*pow(fPsi,2)*rhw_sigma[2]*cos(5*beta) - rhw_sigma[2]*cos(7*beta) + 2*fPsi*rhw_sigma[2]*cos(7*beta) -
701  pow(fPsi,2)*rhw_sigma[2]*cos(7*beta) + 3*rhw_sigma[1]*sin(beta) + 2*fPsi*rhw_sigma[1]*sin(beta) +
702  3*pow(fPsi,2)*rhw_sigma[1]*sin(beta) + 5*rhw_sigma[2]*sin(2*beta) - 5*pow(fPsi,2)*rhw_sigma[2]*sin(2*beta) -
703  rhw_sigma[2]*sin(4*beta) + pow(fPsi,2)*rhw_sigma[2]*sin(4*beta) + 2*rhw_sigma[1]*sin(5*beta) -
704  4*fPsi*rhw_sigma[1]*sin(5*beta) + 2*pow(fPsi,2)*rhw_sigma[1]*sin(5*beta) + rhw_sigma[1]*sin(7*beta) -
705  2*fPsi*rhw_sigma[1]*sin(7*beta) + pow(fPsi,2)*rhw_sigma[1]*sin(7*beta)))/
706  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3)));
707 
708 
709  jacobianf1(0,2) = 0;
710 
711  jacobianf1(1,0) = -((sqrt(2)*exp(fB*i1)*fB*fC*fPsi*(-((3 + 2*fPsi + 3*pow(fPsi,2))*rhw_sigma[2]*cos(beta)) +
712  5*(-1 + pow(fPsi,2))*rhw_sigma[1]*cos(2*beta) + 24*sqrt(2)*fA*fPsi*cos(3*beta) -
713  24*sqrt(2)*exp(fB*i1)*fC*fPsi*cos(3*beta) - 24*sqrt(2)*fA*pow(fPsi,2)*cos(3*beta) +
714  24*sqrt(2)*exp(fB*i1)*fC*pow(fPsi,2)*cos(3*beta) - rhw_sigma[1]*cos(4*beta) +
715  pow(fPsi,2)*rhw_sigma[1]*cos(4*beta) + 2*rhw_sigma[2]*cos(5*beta) - 4*fPsi*rhw_sigma[2]*cos(5*beta) +
716  2*pow(fPsi,2)*rhw_sigma[2]*cos(5*beta) - rhw_sigma[2]*cos(7*beta) + 2*fPsi*rhw_sigma[2]*cos(7*beta) -
717  pow(fPsi,2)*rhw_sigma[2]*cos(7*beta) + 3*rhw_sigma[1]*sin(beta) + 2*fPsi*rhw_sigma[1]*sin(beta) +
718  3*pow(fPsi,2)*rhw_sigma[1]*sin(beta) + 5*rhw_sigma[2]*sin(2*beta) - 5*pow(fPsi,2)*rhw_sigma[2]*sin(2*beta) -
719  rhw_sigma[2]*sin(4*beta) + pow(fPsi,2)*rhw_sigma[2]*sin(4*beta) + 2*rhw_sigma[1]*sin(5*beta) -
720  4*fPsi*rhw_sigma[1]*sin(5*beta) + 2*pow(fPsi,2)*rhw_sigma[1]*sin(5*beta) + rhw_sigma[1]*sin(7*beta) -
721  2*fPsi*rhw_sigma[1]*sin(7*beta) + pow(fPsi,2)*rhw_sigma[1]*sin(7*beta)))/
722  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3)));
723 
724  jacobianf1(1,1) = ((fA - exp(fB*i1)*fC)*fPsi*(8*(fA - exp(fB*i1)*fC)*fPsi*pow(cos(beta),2)*
725  pow(1 + fPsi + 8*(-1 + fPsi)*pow(sin(beta),3),2) +
726  2*sqrt(2)*cos(beta)*(15 - 34*fPsi + 15*pow(fPsi,2) - 6*pow(-1 + fPsi,2)*cos(2*beta) +
727  6*pow(-1 + fPsi,2)*cos(4*beta) + 2*cos(6*beta) - 4*fPsi*cos(6*beta) +
728  2*pow(fPsi,2)*cos(6*beta) + 12*sin(beta) - 12*pow(fPsi,2)*sin(beta) - 13*sin(3*beta) +
729  13*pow(fPsi,2)*sin(3*beta))*(2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*cos(beta) -
730  rhw_sigma[1]*(1 + fPsi + (-1 + fPsi)*sin(3*beta))) +
731  8*(fA - exp(fB*i1)*fC)*fPsi*pow(3*(-1 + fPsi)*cos(beta)*cos(3*beta) +
732  sin(beta)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)),2) +
733  sqrt(2)*(2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*sin(beta) -
734  rhw_sigma[2]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)))*
735  ((-1 + pow(fPsi,2))*cos(2*beta) - 13*(-1 + pow(fPsi,2))*cos(4*beta) +
736  8*(3 - 7*fPsi + 3*pow(fPsi,2))*sin(beta) +
737  2*pow(-1 + fPsi,2)*(-4*sin(5*beta) + sin(7*beta)))))/
738  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),4));
739 
740  jacobianf1(1,2) = 0;
741 
742 
743  jacobianf1(2,0) = -1;
744 
745 
746  jacobianf1(2,1) = 0;
747 
748  jacobianf1(2,2) = -3*CK*(CPer + fD/exp(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k)))*(1 + exp(fB*k)*fB*fC*fR)*fW;
749 
750 
751 
752 
753 
754 
755 }
756 
757 
758 void TPZSandlerExtended::Jacobianf2(const TPZVec<STATE> &trial_stress, STATE theta, STATE beta, STATE k, TPZFMatrix<STATE> &jacobianf2)const {
759 
760 
761  jacobianf2.Resize(3, 3);
762  STATE CX0 = X_0();
763  STATE CPer = CPerturbation();
764  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
765  STATE CG = fE/(2.0*(1.0 + fnu));
766 
767  TPZManVector<REAL,3> rhw_sigma(3);
768  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
769 
770 
771  jacobianf2(0,0) = (2*(fA - exp(fB*k)*fC)*(108*CK*(fA - exp(fB*k)*fC)*pow(fPsi,2)*pow(cos(beta),2)*
772  pow(cos(theta),2) + 108*CK*(fA - exp(fB*k)*fC)*pow(fPsi,2)*pow(cos(theta),2)*
773  pow(sin(beta),2) - sqrt(3)*CG*fR*cos(theta)*
774  (3*rhw_sigma[0] + sqrt(3)*(-k + (fA - exp(fB*k)*fC)*fR*cos(theta)))*
775  pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2) +
776  3*CG*(fA - exp(fB*k)*fC)*pow(fR,2)*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2)*
777  pow(sin(theta),2) + 27*sqrt(2)*CK*fPsi*cos(beta)*sin(theta)*
778  (rhw_sigma[1]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)) -
779  2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*sin(theta)) +
780  27*sqrt(2)*CK*fPsi*sin(beta)*sin(theta)*
781  (rhw_sigma[2]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)) -
782  2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*sin(beta)*sin(theta))))/
783  (27.*CG*CK*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
784 
785 
786 
787  jacobianf2(0,1) = (sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(theta)*
788  (-((3 + 2*fPsi + 3*pow(fPsi,2))*rhw_sigma[2]*cos(beta)) + 5*(-1 + pow(fPsi,2))*rhw_sigma[1]*cos(2*beta) -
789  rhw_sigma[1]*cos(4*beta) + pow(fPsi,2)*rhw_sigma[1]*cos(4*beta) + 2*rhw_sigma[2]*cos(5*beta) -
790  4*fPsi*rhw_sigma[2]*cos(5*beta) + 2*pow(fPsi,2)*rhw_sigma[2]*cos(5*beta) - rhw_sigma[2]*cos(7*beta) +
791  2*fPsi*rhw_sigma[2]*cos(7*beta) - pow(fPsi,2)*rhw_sigma[2]*cos(7*beta) + 3*rhw_sigma[1]*sin(beta) +
792  2*fPsi*rhw_sigma[1]*sin(beta) + 3*pow(fPsi,2)*rhw_sigma[1]*sin(beta) + 5*rhw_sigma[2]*sin(2*beta) -
793  5*pow(fPsi,2)*rhw_sigma[2]*sin(2*beta) - rhw_sigma[2]*sin(4*beta) + pow(fPsi,2)*rhw_sigma[2]*sin(4*beta) +
794  2*rhw_sigma[1]*sin(5*beta) - 4*fPsi*rhw_sigma[1]*sin(5*beta) + 2*pow(fPsi,2)*rhw_sigma[1]*sin(5*beta) +
795  rhw_sigma[1]*sin(7*beta) - 2*fPsi*rhw_sigma[1]*sin(7*beta) + pow(fPsi,2)*rhw_sigma[1]*sin(7*beta) -
796  12*sqrt(2)*fA*fPsi*sin(3*beta - theta) + 12*sqrt(2)*exp(fB*k)*fC*fPsi*sin(3*beta - theta) +
797  12*sqrt(2)*fA*pow(fPsi,2)*sin(3*beta - theta) -
798  12*sqrt(2)*exp(fB*k)*fC*pow(fPsi,2)*sin(3*beta - theta) +
799  12*sqrt(2)*fA*fPsi*sin(3*beta + theta) - 12*sqrt(2)*exp(fB*k)*fC*fPsi*sin(3*beta + theta) -
800  12*sqrt(2)*fA*pow(fPsi,2)*sin(3*beta + theta) +
801  12*sqrt(2)*exp(fB*k)*fC*pow(fPsi,2)*sin(3*beta + theta)))/
802  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3));
803 
804 
805  jacobianf2(0,2) = (2*((9*sqrt(2)*exp(fB*k)*fB*fC*fPsi*rhw_sigma[1]*cos(beta)*cos(theta))/
806  (CG*(1 + fPsi + (-1 + fPsi)*sin(3*beta))) +
807  (9*sqrt(2)*exp(fB*k)*fB*fC*fPsi*rhw_sigma[2]*cos(theta)*sin(beta))/
808  (CG*(1 + fPsi + (-1 + fPsi)*sin(3*beta))) +
809  (sqrt(3)*exp(fB*k)*fB*fC*fR*rhw_sigma[0]*sin(theta))/CK +
810  ((fA - exp(fB*k)*fC)*fR*(1 + exp(fB*k)*fB*fC*fR*cos(theta))*sin(theta))/CK -
811  (exp(fB*k)*fB*fC*fR*(k - (fA - exp(fB*k)*fC)*fR*cos(theta))*sin(theta))/CK +
812  (36*exp(fB*k)*fB*fC*(-fA + exp(fB*k)*fC)*pow(fPsi,2)*pow(cos(beta),2)*sin(2*theta))/
813  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2)) +
814  (36*exp(fB*k)*fB*fC*(-fA + exp(fB*k)*fC)*pow(fPsi,2)*pow(sin(beta),2)*sin(2*theta))/
815  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2))))/9.0;
816 
817  jacobianf2(1,0) = (sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(theta)*
818  (-((3 + 2*fPsi + 3*pow(fPsi,2))*rhw_sigma[2]*cos(beta)) + 5*(-1 + pow(fPsi,2))*rhw_sigma[1]*cos(2*beta) -
819  rhw_sigma[1]*cos(4*beta) + pow(fPsi,2)*rhw_sigma[1]*cos(4*beta) + 2*rhw_sigma[2]*cos(5*beta) -
820  4*fPsi*rhw_sigma[2]*cos(5*beta) + 2*pow(fPsi,2)*rhw_sigma[2]*cos(5*beta) - rhw_sigma[2]*cos(7*beta) +
821  2*fPsi*rhw_sigma[2]*cos(7*beta) - pow(fPsi,2)*rhw_sigma[2]*cos(7*beta) + 3*rhw_sigma[1]*sin(beta) +
822  2*fPsi*rhw_sigma[1]*sin(beta) + 3*pow(fPsi,2)*rhw_sigma[1]*sin(beta) + 5*rhw_sigma[2]*sin(2*beta) -
823  5*pow(fPsi,2)*rhw_sigma[2]*sin(2*beta) - rhw_sigma[2]*sin(4*beta) + pow(fPsi,2)*rhw_sigma[2]*sin(4*beta) +
824  2*rhw_sigma[1]*sin(5*beta) - 4*fPsi*rhw_sigma[1]*sin(5*beta) + 2*pow(fPsi,2)*rhw_sigma[1]*sin(5*beta) +
825  rhw_sigma[1]*sin(7*beta) - 2*fPsi*rhw_sigma[1]*sin(7*beta) + pow(fPsi,2)*rhw_sigma[1]*sin(7*beta) -
826  12*sqrt(2)*fA*fPsi*sin(3*beta - theta) + 12*sqrt(2)*exp(fB*k)*fC*fPsi*sin(3*beta - theta) +
827  12*sqrt(2)*fA*pow(fPsi,2)*sin(3*beta - theta) -
828  12*sqrt(2)*exp(fB*k)*fC*pow(fPsi,2)*sin(3*beta - theta) +
829  12*sqrt(2)*fA*fPsi*sin(3*beta + theta) - 12*sqrt(2)*exp(fB*k)*fC*fPsi*sin(3*beta + theta) -
830  12*sqrt(2)*fA*pow(fPsi,2)*sin(3*beta + theta) +
831  12*sqrt(2)*exp(fB*k)*fC*pow(fPsi,2)*sin(3*beta + theta)))/
832  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3));
833 
834 
835  jacobianf2(1,1) = ((fA - exp(fB*k)*fC)*fPsi*sin(theta)*(8*(fA - exp(fB*k)*fC)*fPsi*
836  pow(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) + (1 + fPsi)*sin(beta),2)*sin(theta)\
837  + 8*(fA - exp(fB*k)*fC)*fPsi*pow(cos(beta),2)*
838  pow(1 + fPsi + 8*(-1 + fPsi)*pow(sin(beta),3),2)*sin(theta) -
839  2*sqrt(2)*(18*pow(-1 + fPsi,2)*cos(beta)*pow(cos(3*beta),2) +
840  6*(-1 + fPsi)*cos(3*beta)*sin(beta)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)) +
841  9*(-1 + fPsi)*cos(beta)*sin(3*beta)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)) -
842  cos(beta)*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2))*
843  (rhw_sigma[1]*(1 + fPsi + (-1 + fPsi)*sin(3*beta)) -
844  2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*sin(theta)) +
845  sqrt(2)*((-1 + pow(fPsi,2))*cos(2*beta) - 13*(-1 + pow(fPsi,2))*cos(4*beta) +
846  8*(3 - 7*fPsi + 3*pow(fPsi,2))*sin(beta) +
847  2*pow(-1 + fPsi,2)*(-4*sin(5*beta) + sin(7*beta)))*
848  (-((1 + fPsi)*rhw_sigma[2]) - 3*(-1 + fPsi)*rhw_sigma[2]*pow(cos(beta),2)*sin(beta) +
849  (-1 + fPsi)*rhw_sigma[2]*pow(sin(beta),3) +
850  2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*sin(beta)*sin(theta))))/
851  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),4));
852 
853 
854  jacobianf2(1,2) = -((sqrt(2)*exp(fB*k)*fB*fC*fPsi*sin(theta)*
855  (-((3 + 2*fPsi + 3*pow(fPsi,2))*rhw_sigma[2]*cos(beta)) +
856  5*(-1 + pow(fPsi,2))*rhw_sigma[1]*cos(2*beta) - rhw_sigma[1]*cos(4*beta) +
857  pow(fPsi,2)*rhw_sigma[1]*cos(4*beta) + 2*rhw_sigma[2]*cos(5*beta) - 4*fPsi*rhw_sigma[2]*cos(5*beta) +
858  2*pow(fPsi,2)*rhw_sigma[2]*cos(5*beta) - rhw_sigma[2]*cos(7*beta) + 2*fPsi*rhw_sigma[2]*cos(7*beta) -
859  pow(fPsi,2)*rhw_sigma[2]*cos(7*beta) + 3*rhw_sigma[1]*sin(beta) + 2*fPsi*rhw_sigma[1]*sin(beta) +
860  3*pow(fPsi,2)*rhw_sigma[1]*sin(beta) + 5*rhw_sigma[2]*sin(2*beta) - 5*pow(fPsi,2)*rhw_sigma[2]*sin(2*beta) -
861  rhw_sigma[2]*sin(4*beta) + pow(fPsi,2)*rhw_sigma[2]*sin(4*beta) + 2*rhw_sigma[1]*sin(5*beta) -
862  4*fPsi*rhw_sigma[1]*sin(5*beta) + 2*pow(fPsi,2)*rhw_sigma[1]*sin(5*beta) + rhw_sigma[1]*sin(7*beta) -
863  2*fPsi*rhw_sigma[1]*sin(7*beta) + pow(fPsi,2)*rhw_sigma[1]*sin(7*beta) -
864  12*sqrt(2)*fA*fPsi*sin(3*beta - theta) +
865  12*sqrt(2)*exp(fB*k)*fC*fPsi*sin(3*beta - theta) +
866  12*sqrt(2)*fA*pow(fPsi,2)*sin(3*beta - theta) -
867  12*sqrt(2)*exp(fB*k)*fC*pow(fPsi,2)*sin(3*beta - theta) +
868  12*sqrt(2)*fA*fPsi*sin(3*beta + theta) -
869  12*sqrt(2)*exp(fB*k)*fC*fPsi*sin(3*beta + theta) -
870  12*sqrt(2)*fA*pow(fPsi,2)*sin(3*beta + theta) +
871  12*sqrt(2)*exp(fB*k)*fC*pow(fPsi,2)*sin(3*beta + theta)))/
872  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),3)));
873 
874  jacobianf2(2,0) = -((fA - exp(fB*k)*fC)*fR*sin(theta));
875 
876  jacobianf2(2,1) = 0;
877 
878 
879  jacobianf2(2,2) = -1 - 3*CK*(CPer + fD/exp(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k)))*(1 + exp(fB*k)*fB*fC*fR)*
880  fW - exp(fB*k)*fB*fC*fR*cos(theta);
881 
882 
883 }
884 
886 void TPZSandlerExtended::Jacobianf2CoVertex(const TPZVec<STATE> &trial_stress, STATE beta, STATE k, TPZFMatrix<STATE> &jacobianf2_covertex)const{
887 
888  jacobianf2_covertex.Resize(2, 2);
889  STATE CX0 = X_0();
890  STATE CPer = CPerturbation();
891  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
892  STATE CG = fE/(2.0*(1.0 + fnu));
893  STATE theta = M_PI_2;
894 
895  TPZManVector<REAL,3> rhw_sigma(3);
896  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
897 
898  jacobianf2_covertex(0,0) = pow((2*sqrt(2)*(fA - exp(fB*k)*fC)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta))/
899  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) -
900  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)),2)/CG +
901  pow((2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi))/
902  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
903  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)),2)/CG +
904  ((rhw_sigma[1] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
905  ((-4*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*pow(3*cos(3*beta) - (3*cos(3*beta))/fPsi,2))/
906  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),3) -
907  (4*sqrt(2)*(fA - exp(fB*k)*fC)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta))/
908  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
909  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)) +
910  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(-9*sin(3*beta) + (9*sin(3*beta))/fPsi))/
911  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2)))/CG +
912  ((rhw_sigma[2] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
913  ((-4*sqrt(2)*(fA - exp(fB*k)*fC)*pow(3*cos(3*beta) - (3*cos(3*beta))/fPsi,2)*sin(beta))/
914  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),3) +
915  (4*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi))/
916  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
917  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)) +
918  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*(-9*sin(3*beta) + (9*sin(3*beta))/fPsi))/
919  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2)))/CG;
920 
921 
922  jacobianf2_covertex(0,1) = (2*sqrt(2)*exp(fB*k)*fB*fC*sin(beta)*((2*sqrt(2)*(fA - exp(fB*k)*fC)*
923  (3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta))/
924  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) -
925  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/
926  (CG*(1 + (1 - sin(3*beta))/fPsi + sin(3*beta))) +
927  ((rhw_sigma[1] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
928  ((-2*sqrt(2)*exp(fB*k)*fB*fC*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi))/
929  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) -
930  (2*sqrt(2)*exp(fB*k)*fB*fC*sin(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/CG +
931  (((-2*sqrt(2)*exp(fB*k)*fB*fC*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta))/
932  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
933  (2*sqrt(2)*exp(fB*k)*fB*fC*cos(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
934  (rhw_sigma[2] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta))/
935  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/CG +
936  (2*sqrt(2)*exp(fB*k)*fB*fC*cos(beta)*
937  ((2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi))/
938  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
939  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/
940  (CG*(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)));
941 
942  jacobianf2_covertex(1,0) = 0;
943 
944  jacobianf2_covertex(1,1) = -1 - 3*CK*(CPer*(1 + exp(fB*k)*fB*fC*fR) +
945  exp(fD*(-CX0 - (fA - exp(fB*k)*fC)*fR + k))*fD*(1 + exp(fB*k)*fB*fC*fR))*fW;
946 
947 }
948 
950 void TPZSandlerExtended::Jacobianf2Vertex(const TPZVec<STATE> &trial_stress, STATE k, STATE &jacobianf2_vertex)const{
951 
952  STATE CX0 = X_0();
953  STATE CPer = CPerturbation();
954  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
955  STATE theta = 0.0;
956 
957 // TPZManVector<REAL,3> rhw_sigma(3);
958 // TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
959 
960  jacobianf2_vertex = -1 - 3*CK*(CPer + fD/exp(fD*(CX0 + fA*fR - exp(fB*k)*fC*fR - k)))*(1 + exp(fB*k)*fB*fC*fR)*fW -
961  exp(fB*k)*fB*fC*fR*cos(theta);
962 }
963 
964 void TPZSandlerExtended::JacobianVertex(const TPZVec<STATE> &trial_stress, STATE k, STATE &jacobian_vertex) const {
965 
966  STATE theta = 0.0;
967  STATE CX0 = X_0();
968  STATE CPer = CPerturbation();
969  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
970  STATE CG = fE/(2.0*(1.0 + fnu));
971 
972  jacobian_vertex = -1 - 3*CK*(CPer*(1 + exp(fB*k)*fB*fC*fR) +
973  exp(fD*(-CX0 - (fA - exp(fB*k)*fC)*fR + k))*fD*(1 + exp(fB*k)*fB*fC*fR))*fW - exp(fB*k)*fB*fC*fR*cos(theta);
974 
975 }
976 
977 void TPZSandlerExtended::JacobianCoVertex(const TPZVec<STATE> &trial_stress, STATE beta, STATE k, TPZFMatrix<STATE> &jacobian_covertex) const {
978 
979  // In this implementation the definition for theta is given by the angle formed from -I1 axis to sqrt(J2) axis with origin on the damage variable kappa.
980  // Thus, covertex is located at theta = M_PI_2;
981  STATE theta = M_PI_2;
982 
983  jacobian_covertex.Resize(2, 2);
984  STATE CX0 = X_0();
985  STATE CPer = CPerturbation();
986  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
987  STATE CG = fE/(2.0*(1.0 + fnu));
988 
989  TPZManVector<REAL,3> rhw_sigma(3);
990  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_sigma);
991 
992  jacobian_covertex(0,0) = pow((2*sqrt(2)*(fA - exp(fB*k)*fC)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta)*sin(theta))/
993  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) -
994  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*sin(theta))/
995  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)),2)/CG +
996  pow((2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*
997  sin(theta))/pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
998  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*sin(theta))/
999  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)),2)/CG +
1000  ((rhw_sigma[1] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*sin(theta))/
1001  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
1002  ((-4*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*pow(3*cos(3*beta) - (3*cos(3*beta))/fPsi,2)*
1003  sin(theta))/pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),3) -
1004  (4*sqrt(2)*(fA - exp(fB*k)*fC)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta)*sin(theta))/
1005  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
1006  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*sin(theta))/
1007  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)) +
1008  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(-9*sin(3*beta) + (9*sin(3*beta))/fPsi)*
1009  sin(theta))/pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2)))/CG +
1010  ((rhw_sigma[2] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*sin(theta))/
1011  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
1012  ((-4*sqrt(2)*(fA - exp(fB*k)*fC)*pow(3*cos(3*beta) - (3*cos(3*beta))/fPsi,2)*sin(beta)*
1013  sin(theta))/pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),3) +
1014  (4*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(theta))/
1015  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
1016  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*sin(theta))/
1017  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)) +
1018  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*(-9*sin(3*beta) + (9*sin(3*beta))/fPsi)*
1019  sin(theta))/pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2)))/CG;
1020 
1021 
1022  jacobian_covertex(0,1) = (2*sqrt(2)*exp(fB*k)*fB*fC*sin(beta)*sin(theta)*
1023  ((2*sqrt(2)*(fA - exp(fB*k)*fC)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta)*sin(theta))/
1024  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) -
1025  (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*sin(theta))/
1026  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/(CG*(1 + (1 - sin(3*beta))/fPsi + sin(3*beta))) +
1027  ((rhw_sigma[1] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*sin(theta))/
1028  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))*
1029  ((-2*sqrt(2)*exp(fB*k)*fB*fC*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(theta))/
1030  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) -
1031  (2*sqrt(2)*exp(fB*k)*fB*fC*sin(beta)*sin(theta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))
1032  )/CG + (((-2*sqrt(2)*exp(fB*k)*fB*fC*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(beta)*
1033  sin(theta))/pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
1034  (2*sqrt(2)*exp(fB*k)*fB*fC*cos(beta)*sin(theta))/(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)))
1035  *(rhw_sigma[2] - (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*sin(theta))/
1036  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/CG +
1037  (2*sqrt(2)*exp(fB*k)*fB*fC*cos(beta)*sin(theta)*
1038  ((2*sqrt(2)*(fA - exp(fB*k)*fC)*cos(beta)*(3*cos(3*beta) - (3*cos(3*beta))/fPsi)*sin(theta))/
1039  pow(1 + (1 - sin(3*beta))/fPsi + sin(3*beta),2) +
1040  (2*sqrt(2)*(fA - exp(fB*k)*fC)*sin(beta)*sin(theta))/
1041  (1 + (1 - sin(3*beta))/fPsi + sin(3*beta))))/(CG*(1 + (1 - sin(3*beta))/fPsi + sin(3*beta)));
1042 
1043 
1044  jacobian_covertex(1,0) = 0;
1045 
1046  jacobian_covertex(1,1) = -1 - 3*CK*(CPer*(1 + exp(fB*k)*fB*fC*fR) +
1047  exp(fD*(-CX0 - (fA - exp(fB*k)*fC)*fR + k))*fD*(1 + exp(fB*k)*fB*fC*fR))*fW -
1048  exp(fB*k)*fB*fC*fR*cos(theta);
1049 
1050 }
1051 
1052 void TPZSandlerExtended::YieldFunction(const TPZVec<STATE> &sigma, STATE kprev, TPZVec<STATE> &yield) const {
1053 
1054  yield.resize(3);
1055  STATE II1, JJ2, ggamma, temp1, temp3, f1, f2, phi, sqrtj2, X, xi, rho, beta;
1056  TPZManVector<STATE, 3> cylstress(3);
1057  TPZHWTools::FromPrincipalToHWCyl(sigma, cylstress);
1058 
1059  // Zylinderkoordinaten
1060  xi = cylstress[0];
1061  rho = cylstress[1];
1062  beta = cylstress[2];
1063 
1064  II1 = sqrt(3.0)*xi;
1065  JJ2 = 0.5*rho*rho;
1066 
1067  if (IsZero(JJ2)) {
1068  JJ2 = 0.0;
1069  }
1070 
1071  sqrtj2 = sqrt(JJ2);
1072  ggamma = 0.5 * (1. + (1. - sin(3. * beta)) / fPsi + sin(3. * beta));
1073 
1074  temp1 = (kprev-II1) / (fR * F(kprev));
1075  temp3 = (ggamma * sqrtj2) / (F(kprev));
1076 
1077  f1 = sqrtj2 - F(II1)/ggamma;
1078  f2 = (temp1 * temp1 + temp3 * temp3 - 1);
1079 
1080  X = kprev - fR * F(kprev);
1081 
1082  // hardcoded
1083  if (II1 > kprev) {
1084  phi = f1;
1085  }else if (II1 > X || IsZero(II1-X) ) {
1086  phi = f2;
1087  }else{
1088  phi = 0.0;
1089  }
1090 
1091  yield[0] = f1;
1092  yield[1] = f2;
1093  yield[2] = phi;
1094 
1095 }
1096 
1097 void TPZSandlerExtended::Print(std::ostream& out) const {
1098  out << "TPZSandlerExtended\n";
1099  out << "A: " << fA << std::endl;
1100  out << "B: " << fB << std::endl;
1101  out << "C: " << fC << std::endl;
1102  out << "D: " << fD << std::endl;
1103  out << "R: " << fR << std::endl;
1104  out << "W: " << fW << std::endl;
1105  out << "k_0: " << fkappa_0 << std::endl;
1106 }
1107 
1108 
1109 void TPZSandlerExtended::Phi(TPZVec<REAL> sigma, STATE alpha, TPZVec<STATE> &phi)const {
1110 
1111  // TPZTensor<REAL>::TPZDecomposed DecompSig;
1112  // TPZTensor<STATE> sig;
1113  // TPZElasticResponse ER;
1114  // ER.Compute(eps,sig);
1115  // sig.EigenSystem(DecompSig);
1116  YieldFunction(sigma, alpha, phi);
1117 }
1118 
1119 std::map<int, int64_t> gF1Stat;
1120 std::map<int, int64_t> gF2Stat;
1121 std::vector<int64_t> gYield;
1122 
1123 STATE TPZSandlerExtended::NormalFunctionToF1(STATE & I1, STATE & k) const{
1124 
1125  STATE normal_to_failure = I1/(exp(fB*k)*fB*fC) + (sqrt(1 + exp(2*fB*k)*pow(fB,2)*pow(fC,2))*
1126  (-((fA - exp(fB*k)*fC + 1/sqrt(1 + exp(2*fB*k)*pow(fB,2)*pow(fC,2)))*k) +
1127  (fA - exp(fB*k)*fC)*((exp(fB*k)*fB*fC)/
1128  sqrt(1 + exp(2*fB*k)*pow(fB,2)*pow(fC,2)) + k)))/(exp(fB*k)*fB*fC);
1129  return normal_to_failure;
1130 
1131 }
1132 
1133 void TPZSandlerExtended::ProjectApex(const TPZVec<STATE> &sigmatrial, STATE kprev, TPZVec<STATE> &sigproj, STATE &kproj) const {
1134 
1135  REAL K = fElasticResponse.K();
1136  REAL xi_apex = Apex()/3.0;
1137 
1138  // Trial
1139  REAL ptr_np1 = 0.;
1140  for (int i = 0; i < 3; i++) {
1141  ptr_np1 += sigmatrial[i];
1142  }
1143  ptr_np1 /= 3.;
1144 
1145  REAL p_np1;
1146  REAL delta_eps_np1 = 0;
1147  REAL res = xi_apex - ptr_np1;
1148  REAL jac;
1149 
1150  bool stop_criterion;
1151  int n_iterations = 50; // @TODO : Define a numeric controls manager object and use it to obtain this information
1152  int i;
1153  for (i = 0; i < n_iterations; i++) {
1154  jac = K;
1155  delta_eps_np1 -= res / jac;
1156  p_np1 = ptr_np1 - K * delta_eps_np1;
1157  res = xi_apex - p_np1;
1158  stop_criterion = std::fabs(res) <= ftol; //IsZero(res);
1159  if (stop_criterion) {
1160  break;
1161  }
1162  }
1163 
1164 #ifdef PZDEBUG
1165  if (i == n_iterations) {
1166 #ifdef WIN32
1167  REAL tol = 1.e-10;
1168 #else
1169  REAL tol = 1.e-12;
1170 #endif
1171  throw TPZConvergenceException(tol, n_iterations, res, i, "TPZSandlerExtended::ProjectApex:: Newton process did not converge.");
1172  }
1173 #endif
1174 
1175  for (int i = 0; i < 3; i++) {
1176  sigproj[i] = p_np1;
1177  }
1178  kproj = kprev;
1179 }
1180 
1181 void TPZSandlerExtended::ProjectF1(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> & projected_stress, STATE &kproj) const {
1182 #ifdef LOG4CXX
1183  if (loggerConvTest->isDebugEnabled()) {
1184  std::stringstream outfile;
1185  outfile << "\n projection over F1 " << endl;
1186  LOGPZ_DEBUG(loggerConvTest, outfile.str());
1187  }
1188 #endif
1189 
1190  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1191  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1192  TPZHWTools::FromPrincipalToHWCyl(trial_stress, hw_space_xi_rho_beta);
1193  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_space_s1_s2_s3);
1194 
1195  STATE i1_guess = sqrt(3.0)*rhw_space_s1_s2_s3[0];
1196  STATE beta_guess = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
1197  STATE k_guess = sqrt(3.0)*rhw_space_s1_s2_s3[0];
1198 
1199  TPZFNMatrix<3, STATE> delta_par(3, 1, 0.), par(3, 1, 0.), residue(3, 1, 0.);
1200  par(0, 0) = i1_guess;
1201  par(1, 0) = beta_guess;
1202  par(2, 0) = k_guess;
1203 
1204  TPZFNMatrix<9, STATE> jac(3, 3);
1205  TPZFNMatrix<9, STATE> jac_inv(3,3);
1206 
1207  TPZManVector<STATE,3> residue_vec(3);
1208  STATE residue_norm;
1209  bool stop_criterion_res;
1210  int max_iterations = 50;
1211  int it;
1212  for (it = 0; it < max_iterations; it++) {
1213  // Computing the Residue vector for a Newton step
1214  Res1(trial_stress, par(0), par(1), par(2), kprev, residue_vec); // Residue
1215  for (int k = 0; k < 3; k++) residue(k, 0) = - 1.0 * residue_vec[k]; // Transfering to a Matrix object
1216  residue_norm = Norm(residue);
1217  stop_criterion_res = std::fabs(residue_norm) <= ftol;
1218  if (stop_criterion_res) {
1219  break;
1220  }
1221 
1222  Jacobianf1(trial_stress, par(0), par(1), par(2), jac); // Jacobian
1223  TPZHWTools::A3x3Inverse(jac, jac_inv);
1224  jac_inv.Multiply(residue, delta_par);
1225 
1226  par += delta_par;
1227  }
1228 
1229 #ifdef PZDEBUG
1230  if (it == max_iterations) {
1231  throw TPZConvergenceException(ftol, max_iterations, residue_norm, it, "TPZSandlerExtended::ProjectF1:: Newton process did not converge.");
1232  }
1233 #endif
1234 
1235  STATE xi, beta;
1236  xi = par(0)/sqrt(3.0);
1237  beta = par(1);
1238  kproj = par(2);
1239 
1240  TPZManVector<STATE, 3> f1cyl(3);
1241  F1Cyl(xi, beta, f1cyl);
1242  TPZHWTools::FromHWCylToPrincipal(f1cyl, projected_stress);
1243 
1244 }
1245 
1246 void TPZSandlerExtended::ProjectF2(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const {
1247 
1248  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1249  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1250  TPZHWTools::FromPrincipalToHWCyl(trial_stress, hw_space_xi_rho_beta);
1251  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_space_s1_s2_s3);
1252 
1253  STATE k_guess = kprev;
1254  STATE beta_guess = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
1255  STATE theta_guess = atan2(hw_space_xi_rho_beta[1],(k_guess/sqrt(3))-sqrt(3)*hw_space_xi_rho_beta[0]);
1256 
1257  bool cap_vertex_validity_Q = theta_guess < ftol;
1258  if (cap_vertex_validity_Q) {
1259 // std::cout << "Projecting on Cap Vertex " << std::endl;
1260  ProjectCapVertex(trial_stress, kprev, projected_stress, kproj);
1261  return;
1262  }
1263 
1264  if (theta_guess > M_PI_2) { // Restriction for theta
1265 // std::cout << "Reached restriction for theta guess = " << theta_guess*(180.0/M_PI) << std::endl;
1266 // std::cout << "Reached restriction for beta guess = " << beta_guess*(180.0/M_PI) << std::endl;
1267  theta_guess = M_PI_2;
1268  }
1269 
1270  TPZFNMatrix<3, STATE> delta_par(3, 1, 0.), par(3, 1, 0.), residue(3, 1, 0.);
1271  par(0, 0) = theta_guess;
1272  par(1, 0) = beta_guess;
1273  par(2, 0) = k_guess;
1274 
1275  TPZFNMatrix<9, STATE> jac(3, 3);
1276  TPZFNMatrix<9, STATE> jac_inv(3,3);
1277 
1278  TPZManVector<STATE,3> residue_vec(3);
1279  STATE residue_norm;
1280  bool stop_criterion_res;
1281  int max_iterations = 50;
1282  int it;
1283  for (it = 0; it < max_iterations; it++) {
1284  // Computing the Residue vector for a Newton step
1285  Res2(trial_stress, par(0), par(1), par(2), kprev, residue_vec); // Residue
1286  for (int k = 0; k < 3; k++) residue(k, 0) = - 1.0 * residue_vec[k]; // Transfering to a Matrix object
1287 
1288  residue_norm = Norm(residue);
1289  stop_criterion_res = std::fabs(residue_norm) <= ftol;
1290  if (stop_criterion_res) {
1291  break;
1292  }
1293 
1294  Jacobianf2(trial_stress, par(0), par(1), par(2), jac); // Jacobian
1295  TPZHWTools::A3x3Inverse(jac, jac_inv);
1296  jac_inv.Multiply(residue, delta_par);
1297  par += delta_par;
1298 
1299  if (par(0) > M_PI_2) { // Restriction for theta
1300 // std::cout << "Reached restriction for theta = " << par(0)*(180.0/M_PI) << std::endl;
1301 // std::cout << "Reached restriction for beta = " << par(1)*(180.0/M_PI) << std::endl;
1302 // std::cout << "Reached restriction for k = " << par(2) << std::endl;
1303  par(0) = M_PI_2;
1304  }
1305  }
1306 #ifdef PZDEBUG
1307  if (it == max_iterations) {
1308  throw TPZConvergenceException(ftol, max_iterations, residue_norm, it, "TPZSandlerExtended::ProjectF2:: Newton process did not converge.");
1309  }
1310 #endif
1311 
1312  STATE theta, beta;
1313  theta = par(0);
1314  beta = par(1);
1315  kproj = par(2);
1316 
1317  TPZManVector<STATE, 3> f2cyl(3);
1318  F2Cyl(theta, beta, kproj, f2cyl);
1319  TPZHWTools::FromHWCylToPrincipal(f2cyl, projected_stress);
1320 }
1321 
1322 void TPZSandlerExtended::ProjectCoVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const{
1323 
1324  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1325  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1326  TPZHWTools::FromPrincipalToHWCyl(trial_stress, hw_space_xi_rho_beta);
1327  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_space_s1_s2_s3);
1328 
1329  STATE k_guess = kprev;
1330  STATE beta_guess = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
1331  STATE theta_guess = M_PI_2;
1332 
1333  TPZFNMatrix<2, STATE> delta_par(2, 1, 0.), par(2, 1, 0.), residue(2, 1, 0.);
1334  par(0, 0) = beta_guess;
1335  par(1, 0) = k_guess;
1336 
1337  TPZFNMatrix<9, STATE> jac(2, 2);
1338  TPZFNMatrix<9, STATE> jac_inv(2,2);
1339 
1340  TPZManVector<STATE,3> residue_vec(2);
1341  STATE residue_norm;
1342  bool stop_criterion_res;
1343  int max_iterations = 50;
1344  int it;
1345  for (it = 0; it < max_iterations; it++) {
1346 
1347  // Computing the Residue vector for a Newton step
1348  Res2CoVertex(trial_stress, par(0), par(1), kprev, residue_vec); // Residue
1349  for (int k = 0; k < 2; k++) residue(k, 0) = - 1.0 * residue_vec[k]; // Transfering to a Matrix object
1350 
1351  residue_norm = Norm(residue);
1352  stop_criterion_res = std::fabs(residue_norm) <= ftol;
1353  if (stop_criterion_res) {
1354  break;
1355  }
1356 
1357  JacobianCoVertex(trial_stress, par(0), par(1), jac); // Jacobian
1358  TPZHWTools::A2x2Inverse(jac, jac_inv);
1359  jac_inv.Multiply(residue, delta_par);
1360  par += delta_par;
1361  }
1362 
1363 #ifdef PZDEBUG
1364  if (it == max_iterations) {
1365  throw TPZConvergenceException(ftol, max_iterations, residue_norm, it, "TPZSandlerExtended::ProjectCoVertex:: Newton process did not converge.");
1366  }
1367 #endif
1368 
1369  STATE theta, beta;
1370  theta = M_PI_2;
1371  beta = par(0);
1372  kproj = par(1);
1373 
1374  TPZManVector<STATE, 3> f2cyl(3);
1375  F2Cyl(theta, beta, kproj, f2cyl);
1376  TPZHWTools::FromHWCylToPrincipal(f2cyl, projected_stress);
1377 
1378 }
1379 
1380 void TPZSandlerExtended::ProjectVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const{
1381 
1382  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1383  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1384  TPZHWTools::FromPrincipalToHWCyl(trial_stress, hw_space_xi_rho_beta);
1385  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_space_s1_s2_s3);
1386 
1387  STATE k = kprev;
1388  STATE dk;
1389 
1390  STATE jac, jac_inv;
1391  STATE residue;
1392  bool stop_criterion_res;
1393  int max_iterations = 50;
1394  int it;
1395  for (it = 0; it < max_iterations; it++) {
1396  // Computing the Residue vector for a Newton step
1397  Res2Vertex(trial_stress, k, kprev, residue); // Residue
1398  stop_criterion_res = std::fabs(residue) <= ftol;
1399  if (stop_criterion_res) {
1400  break;
1401  }
1402 
1403  JacobianVertex(trial_stress, k, jac); // Jacobian
1404  jac_inv = 1.0/jac;
1405  dk = jac_inv * residue;
1406  k += dk;
1407  }
1408 #ifdef PZDEBUG
1409  if (it == max_iterations) {
1410  throw TPZConvergenceException(ftol, max_iterations, stop_criterion_res, it, "TPZSandlerExtended::ProjectVertex:: Newton process did not converge.");
1411  }
1412 #endif
1413 
1414  TPZManVector<STATE, 3> f2cyl(3);
1415  f2cyl[0] = (k - fR * F(k))/sqrt(3.0); // xi
1416  f2cyl[1] = 0.0; // rho
1417  f2cyl[2] = 0.0; // beta
1418  TPZHWTools::FromHWCylToPrincipal(f2cyl, projected_stress);
1419 
1420 }
1421 
1422 void TPZSandlerExtended::ProjectCapVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const{
1423 
1424  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1425  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1426  TPZHWTools::FromPrincipalToHWCyl(trial_stress, hw_space_xi_rho_beta);
1427  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_space_s1_s2_s3);
1428 
1429  STATE k = kprev, dk;
1430 
1431  TPZManVector<STATE,3> residue_vec(3);
1432  STATE residue, jac, jac_inv;
1433  bool stop_criterion_res;
1434  int max_iterations = 50;
1435  int it;
1436  for (it = 0; it < max_iterations; it++) {
1437  // Computing the Residue vector for a Newton step
1438  Res2Vertex(trial_stress, k, kprev, residue); // Residue
1439 
1440  stop_criterion_res = std::fabs(residue) <= ftol;
1441  if (stop_criterion_res) {
1442  break;
1443  }
1444 
1445  Jacobianf2Vertex(trial_stress, k, jac); // Jacobian
1446  jac_inv = 1.0 / jac;
1447  dk = - residue * jac_inv;
1448  k += dk;
1449  }
1450 #ifdef PZDEBUG
1451  if (it == max_iterations) {
1452  throw TPZConvergenceException(ftol, max_iterations, residue, it, "TPZSandlerExtended::ProjectCapVertex:: Newton process did not converge.");
1453  }
1454 #endif
1455 
1456 
1457  STATE theta, beta;
1458  theta = 0.0;
1459  beta = atan2(0.0, 0.0);
1460  kproj = k;
1461 
1462  TPZManVector<STATE, 3> f2cyl(3);
1463  F2Cyl(theta, beta, kproj, f2cyl); // The definition for theta was corrected.
1464  TPZHWTools::FromHWCylToPrincipal(f2cyl, projected_stress);
1465 
1466 }
1467 
1468 void TPZSandlerExtended::ProjectCapCoVertex(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj) const{
1469 
1470  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1471  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1472  TPZHWTools::FromPrincipalToHWCyl(trial_stress, hw_space_xi_rho_beta);
1473  TPZHWTools::FromPrincipalToHWCart(trial_stress, rhw_space_s1_s2_s3);
1474 
1475  STATE k_guess = kprev;
1476  STATE beta_guess = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
1477 
1478  TPZFNMatrix<2, STATE> delta_par(2, 1, 0.), par(2, 1, 0.), residue(2, 1, 0.);
1479  par(0, 0) = beta_guess;
1480  par(1, 0) = k_guess;
1481 
1482  TPZFNMatrix<9, STATE> jac(2, 2);
1483  TPZFNMatrix<9, STATE> jac_inv(2,2);
1484 
1485  TPZManVector<STATE,2> residue_vec(2);
1486  STATE residue_norm;
1487  bool stop_criterion_res;
1488  int max_terations = 50;
1489  int it;
1490  for (it = 0; it < max_terations; it++) {
1491  // Computing the Residue vector for a Newton step
1492  Res2CoVertex(trial_stress, par(0), par(1), kprev, residue_vec); // Residue
1493  for (int k = 0; k < 2; k++) residue(k, 0) = - 1.0 * residue_vec[k]; // Transfering to a Matrix object
1494 
1495  residue_norm = Norm(residue);
1496  stop_criterion_res = std::fabs(residue_norm) <= ftol;
1497  if (stop_criterion_res) {
1498  break;
1499  }
1500 
1501  Jacobianf2CoVertex(trial_stress, par(0), par(1), jac); // Jacobian
1502  TPZHWTools::A2x2Inverse(jac, jac_inv);
1503  jac_inv.Multiply(residue, delta_par);
1504  par += delta_par;
1505  }
1506 #ifdef PZDEBUG
1507  if (it == max_terations) {
1508  throw TPZConvergenceException(ftol, max_terations, residue_norm, it, "TPZSandlerExtended::ProjectCapCoVertex:: Newton process did not converge.");
1509  }
1510 #endif
1511 
1512  STATE theta, beta;
1513  theta = M_PI_2;
1514  beta = par(0);
1515  kproj = par(1);
1516 
1517  TPZManVector<STATE, 3> f2cyl(3);
1518  F2Cyl(theta, beta, kproj, f2cyl); // The definition for theta was corrected.
1519  TPZHWTools::FromHWCylToPrincipal(f2cyl, projected_stress);
1520 
1521 }
1522 
1523 void TPZSandlerExtended::ProjectRing(const TPZVec<STATE> &sigmatrial, STATE kprev, TPZVec<STATE> &sigproj, STATE &kproj) const {
1524 #ifdef LOG4CXX
1525  if (loggerConvTest->isDebugEnabled()) {
1526  std::stringstream outfile;
1527  outfile << "\n projection over Ring " << endl;
1528  LOGPZ_DEBUG(loggerConvTest, outfile.str());
1529  }
1530 #endif
1531  STATE beta = 0., distnew;
1532  STATE resnorm, disttheta;
1533  disttheta = 1.e8;
1534 
1535  for (STATE betaguess = 0; betaguess <= 2 * M_PI; betaguess += M_PI / 20.) {
1536  distnew = DistF2(sigmatrial, M_PI / 2, betaguess, kprev);
1537  if (fabs(distnew) < fabs(disttheta)) {
1538  // STATE theta = M_PI/2;
1539  beta = betaguess;
1540  disttheta = distnew;
1541  }
1542  }
1543  resnorm = 1;
1544  int64_t counter = 1;
1545  TPZFMatrix<STATE> xn1(3, 1, 0.), xn(3, 1, 0.), fxn(3, 1, 0.), diff(3, 1, 0.);
1546  TPZFNMatrix<3, STATE> sol(3, 1, 0.);
1547  xn(0, 0) = M_PI / 2;
1548  xn(1, 0) = beta;
1549  xn(2, 0) = kprev;
1550  while (resnorm > ftol && counter < 30) {
1551  TPZFNMatrix<9, STATE> jac(3, 3);
1552  Jacobianf2(sigmatrial, xn[0], xn[1], xn[2], jac);
1553  TPZManVector<STATE> fxnvec(3);
1554  Res2(sigmatrial, xn(0), xn(1), xn(2), kprev, fxnvec);
1555  for (int k = 0; k < 3; k++) fxn(k, 0) = fxnvec[k];
1556 
1557  for (int i = 0; i < 3; i++) {
1558  jac(i, 0) = 0.;
1559  jac(0, i) = 0.;
1560  }
1561  jac(0, 0) = 1.;
1562  fxn(0, 0) = 0.;
1563  sol = fxn;
1564  resnorm = Norm(sol);
1565  jac.Solve_LU(&sol);
1566 
1567  xn1(0, 0) = xn(0, 0);
1568  xn1(1, 0) = xn(1, 0) - sol(1, 0);
1569  xn1(2, 0) = xn(2, 0) - sol(2, 0);
1570 
1571 #ifdef LOG4CXX
1572  if (loggerConvTest->isDebugEnabled()) {
1573  std::stringstream outfile; //("convergencF1.txt");
1574  outfile << counter << " " << log(resnorm) << endl;
1575  //jac.Print(outfile);
1576  //outfile<< "\n xn " << " "<<fxnvec <<endl;
1577  //outfile<< "\n res " << " "<<fxnvec <<endl;
1578  LOGPZ_DEBUG(loggerConvTest, outfile.str());
1579  }
1580 #endif
1581  //diff=xn1-xn;
1582  //resnorm=Norm(diff);
1583  xn = xn1;
1584  counter++;
1585 
1586  }
1587 #ifdef PZDEBUG
1588  if (counter == 30) {
1589  throw TPZConvergenceException(ftol, 30, resnorm, counter, "TPZSandlerExtended::ProjectRing:: Newton process did not converge.");
1590  }
1591 #endif
1592  // cout<< "\n resnorm = "<<resnorm <<endl;
1593  // cout<< "\n counter = "<<xn1 <<endl;
1594  // cout<< "\n k = "<<xn1[2] <<endl;
1595  STATE thetasol, betasol, ksol;
1596 
1597  thetasol = xn1[0];
1598  betasol = xn1[1];
1599  ksol = xn1[2];
1600 
1601  TPZManVector<STATE, 3> f2cyl(3);
1602  F2Cyl(thetasol, betasol, ksol, f2cyl);
1603  TPZHWTools::FromHWCylToPrincipal(f2cyl, sigproj);
1604 
1605  kproj = ksol;
1606 
1607 }
1608 
1609 void TPZSandlerExtended::ProjectBetaConstF2(const TPZVec<STATE> &sigmatrial, STATE kprev, TPZVec<STATE> &sigproj, STATE &kproj) const {
1610  //#ifdef LOG4CXX
1611  // {
1612  // std::stringstream outfile;
1613  // outfile << "\n projection over F2 " <<endl;
1614  // LOGPZ_DEBUG(logger,outfile.str());
1615  //
1616  // }
1617  //#endif
1618 
1619  STATE theta = 0., beta = 0., distnew;
1620  STATE resnorm, disttheta;
1621  disttheta = 1.e8;
1622  STATE betaconst = 0;
1623  for (STATE thetaguess = 0; thetaguess <= M_PI; thetaguess += M_PI / 20.) {
1624  distnew = DistF2(sigmatrial, thetaguess, betaconst, kprev);
1625  if (fabs(distnew) < fabs(disttheta)) {
1626  theta = thetaguess;
1627  beta = betaconst;
1628  disttheta = distnew;
1629  }
1630  }
1631 
1632  resnorm = 1;
1633  int counter = 1;
1634  TPZFNMatrix<3, STATE> xn1(3, 1, 0.), xn(3, 1, 0.), sol(3, 1, 0.), fxn(3, 1, 0.), diff(3, 1, 0.);
1635  xn(0, 0) = theta;
1636  xn(1, 0) = beta;
1637  xn(2, 0) = kprev;
1638  while (resnorm > ftol && counter < 30) {
1639  TPZFNMatrix<9, STATE> jac(3, 3);
1640  Jacobianf2(sigmatrial, xn(0), xn(1), xn(2), jac);
1641  TPZManVector<STATE> fxnvec(3);
1642  Res2(sigmatrial, xn(0), xn(1), xn(2), kprev, fxnvec);
1643  for (int k = 0; k < 3; k++) fxn(k, 0) = fxnvec[k];
1644  for (int i = 0; i < 3; i++) {
1645  jac(i, 1) = 0.;
1646  jac(1, i) = 0.;
1647  }
1648  jac(1, 1) = 1.;
1649  fxn(1, 0) = 0.;
1650  sol = fxn;
1651  resnorm = Norm(sol);
1652  jac.Solve_LU(&sol);
1653 
1654  xn1(0, 0) = xn(0, 0) - sol(0, 0);
1655  xn1(1, 0) = xn(1, 0);
1656  xn1(2, 0) = xn(2, 0) - sol(2, 0);
1657 
1658  // diff=xn1-xn;
1659  // resnorm=Norm(diff);
1660  xn = xn1;
1661  counter++;
1662 
1663  }
1664  //if(counter == 30) cout << "resnorm = " << resnorm << std::endl;
1665  STATE thetasol, betasol, ksol;
1666 
1667 
1668 
1669 
1670  thetasol = xn1(0);
1671  betasol = xn1(1);
1672  ksol = xn1(2);
1673  kproj = ksol;
1674 
1675  TPZManVector<STATE, 3> f2cyl(3);
1676  F2Cyl(thetasol, betasol, ksol, f2cyl);
1677  TPZHWTools::FromHWCylToPrincipal(f2cyl, sigproj);
1678 }
1679 
1680 void TPZSandlerExtended::ComputeI1(TPZVec<STATE> stress, STATE &I1)const {
1681  STATE sig1, sig2, sig3;
1682  sig1 = stress[0];
1683  sig2 = stress[1];
1684  sig3 = stress[2];
1685  I1 = sig1 + sig2 + sig3;
1686 
1687 }
1688 
1689 void TPZSandlerExtended::ComputeJ2(TPZVec<STATE> stress, STATE &J2)const {
1690  STATE sig1, sig2, sig3;
1691  sig1 = stress[0];
1692  sig2 = stress[1];
1693  sig3 = stress[2];
1694  J2 = (2. * sig1 * sig2 + pow(sig1 + (-sig1 - sig2 - sig3) / 3., 2.) +
1695  pow(sig2 + (-sig1 - sig2 - sig3) / 3., 2.) + 2 * sig1 * sig3 + 2. * sig2 * sig3 +
1696  pow((-sig1 - sig2 - sig3) / 3. + sig3, 2.)) / 2.;
1697 }
1698 
1699 void TPZSandlerExtended::ApplyStrainComputeElasticStress(TPZVec<STATE> &strain, TPZVec<STATE> &stress)const {
1700  STATE sig1, sig2, sig3, s1, s2, s3;
1701  sig1 = strain[0];
1702  sig2 = strain[1];
1703  sig3 = strain[2];
1704 
1705  s1 = sig1 - (1. / 3.)*(sig1 + sig2 + sig3);
1706  s2 = sig2 - (1. / 3.)*(sig1 + sig2 + sig3);
1707  s3 = sig3 - (1. / 3.)*(sig1 + sig2 + sig3);
1708 
1709  stress[0] = s1 * (2 * fG) + fK * (sig1 + sig2 + sig3);
1710  stress[1] = s2 * (2 * fG) + fK * (sig1 + sig2 + sig3);
1711  stress[2] = s3 * (2 * fG) + fK * (sig1 + sig2 + sig3);
1712 }
1713 
1714 void TPZSandlerExtended::ApplyStressComputeElasticStrain(TPZVec<STATE> &stress, TPZVec<STATE> &strain)const {
1715  STATE sig1, sig2, sig3, s1, s2, s3;
1716  sig1 = stress[0];
1717  sig2 = stress[1];
1718  sig3 = stress[2];
1719 
1720  s1 = sig1 - (1. / 3.)*(sig1 + sig2 + sig3);
1721  s2 = sig2 - (1. / 3.)*(sig1 + sig2 + sig3);
1722  s3 = sig3 - (1. / 3.)*(sig1 + sig2 + sig3);
1723 
1724  strain[0] = s1 / (2. * fG)+(sig1 + sig2 + sig3) / (9. * fK);
1725  strain[1] = s2 / (2. * fG)+(sig1 + sig2 + sig3) / (9. * fK);
1726  strain[2] = s3 / (2. * fG)+(sig1 + sig2 + sig3) / (9. * fK);
1727 }
1728 
1735 void TPZSandlerExtended::ApplyStrainComputeSigma(TPZVec<STATE> &epst, TPZVec<STATE> &epsp, STATE & kprev, TPZVec<STATE> &epspnext, TPZVec<STATE> &stressnext, STATE & knext, TPZFMatrix<REAL> * tangent) const {
1736 
1737  bool require_tangent_Q = true;
1738  if (!tangent) {
1739  require_tangent_Q = false;
1740  }
1741 
1742 #ifdef PZDEBUG
1743  // Check for required dimensions of tangent
1744  if (!(tangent->Rows() == 6 && tangent->Cols() == 6)) {
1745  std::cerr << "Unable to compute the tangent operator. Required tangent array dimensions are 6x6." << std::endl;
1746  DebugStop();
1747  }
1748 #endif
1749 
1750  if (require_tangent_Q) {
1751  DebugStop(); // implemented this functionality.
1752  }
1753 
1754  STATE trial_I1, I1proj;
1755 
1756  TPZManVector<STATE, 3> trial_stress(3), yield(2), deltastress(3), delepsp(3), epsT(epst);
1757  epsT[0] -= epsp[0];
1758  epsT[1] -= epsp[1];
1759  epsT[2] -= epsp[2];
1760  ApplyStrainComputeElasticStress(trial_stress, epsT);
1761  YieldFunction(trial_stress, kprev, yield);
1762  ComputeI1(trial_stress, trial_I1);
1763 
1764  if ((yield[1] <= 0 && trial_I1 <= kprev) || (yield[0] <= 0 && trial_I1 > kprev)) {
1765  epspnext = epsp;
1766  knext = kprev;
1767  stressnext = trial_stress;
1768  //cout<<"\n elastic "<<endl;
1769  } else {
1770  //cout<<"\n plastic "<<endl;
1771  if (yield[1] > 0 && trial_I1 < kprev) {
1772  //cout<<"\n F2 "<<endl;
1773  ProjectF2(trial_stress, kprev, stressnext, knext);
1774  } else {
1775  //cout<<"\n F1 "<<endl;
1776  ProjectF1(trial_stress, kprev, stressnext, knext);
1777  ComputeI1(stressnext, I1proj);
1778  if (I1proj < knext) {
1779  //cout<<"\n Ring "<<endl;
1780  ProjectRing(trial_stress, kprev, stressnext, knext);
1781  }
1782  }
1783 
1784  for (int i = 0; i < 3; i++) {
1785  deltastress[i] = trial_stress[i] - stressnext[i];
1786  }
1787 
1788  ApplyStressComputeElasticStrain(deltastress, delepsp);
1789 
1790  for (int i = 0; i < 3; i++) {
1791  epspnext[i] = epsp[i] + delepsp[i];
1792  }
1793  }
1794 }
1795 
1796 void TPZSandlerExtended::ProjectSigma(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &sigproj, STATE &kproj, int &m_type, TPZFMatrix<REAL> * gradient) const {
1797 
1798  bool require_gradient_Q = (gradient != NULL);
1799 
1800 #ifdef PZDEBUG
1801  if (require_gradient_Q) {
1802  // Check for required dimensions of tangent
1803  if (!(gradient->Rows() == 3 && gradient->Cols() == 3)) {
1804  std::cerr << "Unable to compute the gradient operator. Required gradient array dimensions are 3x3." << std::endl;
1805  DebugStop();
1806  }
1807  }
1808 #endif
1809 
1810  //Firstk(epspv,k0);
1811  TPZManVector<STATE, 2> yield(2);
1812  STATE I1 = sigtrial[0] + sigtrial[1] + sigtrial[2];
1813  REAL J2 = (1.0/3.0) * (sigtrial[0]*sigtrial[0] + sigtrial[1]*sigtrial[1] + sigtrial[2]*sigtrial[2] - sigtrial[1]*sigtrial[2] - sigtrial[0]*sigtrial[2] - sigtrial[0]*sigtrial[1]);
1814 
1815  YieldFunction(sigtrial, kprev, yield);
1816 
1817  if (I1 < kprev) {
1818  if (yield[1] > 0.0 || IsZero(yield[1]) ) {
1819 
1820  m_type = 2; // cap behavior
1821 // std::cout << "Projecting on Cap " << std::endl;
1822  ProjectF2(sigtrial, kprev, sigproj, kproj);
1823 
1824  STATE proj_i1 = sigproj[0] + sigproj[1] + sigproj[2];
1825 #ifdef PZDEBUG
1826  if (proj_i1 > kproj) {
1827  std::ostringstream stringbuilder("TPZSandlerExtended::ProjectSigma: proj_i1 > kproj with proj_i1 = ");
1828  stringbuilder << proj_i1 << " and kproj = " << kproj << ".";
1829  throw TPZInconsistentStateException(stringbuilder.str());
1830  }
1831 #endif
1832  if (require_gradient_Q) {
1833  ComputeCapTangent(sigtrial, kprev, sigproj, kproj, gradient);
1834  }
1835  return;
1836 
1837  } else {
1838  m_type = 0; // Elastic behaviour
1839  sigproj = sigtrial;
1840  kproj = kprev;
1841  if (require_gradient_Q) {
1842  gradient->Identity();
1843  }
1844  return;
1845  }
1846  } else {
1847 
1848  if (yield[0] > 0.0 || IsZero(yield[0]) ) {
1849 
1850  bool failure_validity_Q = I1 > kprev || IsZero(I1 - kprev);
1851  if (failure_validity_Q) {
1852  STATE normal_to_f1_at_last_k = NormalToF1(I1, kprev);
1853  bool covertex_validity_Q = normal_to_f1_at_last_k < sqrt(J2) || IsZero(normal_to_f1_at_last_k - sqrt(J2));
1854  if (covertex_validity_Q) {
1855  m_type = 2; // cap behavior
1856  ProjectCapCoVertex(sigtrial, kprev, sigproj, kproj);
1857  if (require_gradient_Q) {
1858  ComputeCapCoVertexTangent(sigtrial, kprev, sigproj, kproj, gradient);
1859  }
1860  return;
1861  }else{
1862 
1863  REAL apex_i1 = Apex();
1864  bool inside_apex_region_Q = I1 > apex_i1 || IsZero(I1 - apex_i1);
1865  if (inside_apex_region_Q) {
1866  STATE normal_to_f1_at_appex = NormalToF1(I1, apex_i1);
1867  bool apex_validity_Q = normal_to_f1_at_appex > sqrt(J2) || IsZero(normal_to_f1_at_appex - sqrt(J2));
1868  if (apex_validity_Q) { // Tensile behavior
1869  m_type = -1;
1870 // std::cout << "Projecting on Apex " << std::endl;
1871  ProjectApex(sigtrial, kprev, sigproj, kproj);
1872  if (require_gradient_Q) {
1873  gradient->Identity(); // Because tangent here is zero matrix, it is preferred use the elastic one.
1874  }
1875  return;
1876  }
1877 
1878  m_type = 1; // failure behavior
1879 // std::cout << "Projecting on Failure " << std::endl;
1880  ProjectF1(sigtrial, kprev, sigproj, kproj);
1881  if (require_gradient_Q) {
1882  ComputeFailureTangent(sigtrial, kprev, sigproj, kproj, gradient);
1883  }
1884  return;
1885 
1886  }
1887  else{
1888  m_type = 1; // failure behavior
1889  ProjectF1(sigtrial, kprev, sigproj, kproj);
1890  if (require_gradient_Q) {
1891  ComputeFailureTangent(sigtrial, kprev, sigproj, kproj, gradient);
1892  }
1893  return;
1894  }
1895  }
1896  }
1897 
1898  DebugStop();
1899 
1900  } else {
1901  m_type = 0; // elastic behaviour
1902  sigproj = sigtrial;
1903  kproj = kprev;
1904  if (require_gradient_Q) {
1905  gradient->Identity();
1906  }
1907  return;
1908  }
1909  }
1910 }
1911 
1912 
1913 void TPZSandlerExtended::ComputeCapTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const {
1914 
1915  TPZFMatrix<STATE> Rot(3,3,0.0);
1916  TPZHWTools::GetRotMatrix(Rot);
1917 
1918 
1919  STATE fv = F(kproj);
1920  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
1921  STATE CG = fE/(2.0*(1.0 + fnu));
1922 
1923  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
1924  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
1925  TPZHWTools::FromPrincipalToHWCyl(projected_stress, hw_space_xi_rho_beta);
1926  TPZHWTools::FromPrincipalToHWCart(projected_stress, rhw_space_s1_s2_s3);
1927 
1928  STATE k = kproj;
1929  STATE i1v = sqrt(3)*hw_space_xi_rho_beta[0];
1930  STATE ratio = (k-i1v)/(fR*fv);
1931 
1932  bool cap_vertex_validity_Q = fabs(ratio - 1.0) <= ftol;
1933  if (cap_vertex_validity_Q) {
1934  ComputeCapVertexTangent(trial_stress, kprev, projected_stress, kproj, gradient);
1935  return;
1936  }
1937 
1938  STATE theta = acos(ratio);
1939  STATE beta = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
1940 
1941  // Computing the jacobian and the inverse
1942  TPZFMatrix<STATE> jac(3,3,0.0),jac_inv(3,3,0.0);
1943  Jacobianf2(trial_stress, theta, beta, k, jac); // Jacobian
1944  TPZHWTools::A3x3Inverse(jac, jac_inv); // Jacobian inverse
1945 
1946  TPZFMatrix<STATE> d_res_d_sig_trial(3,3,0.0);
1947  TPZFMatrix<STATE> d_sig_proj_d_theta_beta_kappa(3,3,0.0);
1948  TPZFMatrix<STATE> d_theta_beta_kappa_d_sig_trial(3,3,0.0);
1949  TPZFMatrix<STATE> d_sig_proj_d_sig_trial(3,3,0.0);
1950 
1951  // Derivative for the cap residue respect to trial stresses
1952  d_res_d_sig_trial(0,0) = (-2*(fA - exp(fB*k)*fC)*fR*sin(theta))/(3.*sqrt(3)*CK);
1953  d_res_d_sig_trial(0,1) = (-2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*cos(theta))/(CG*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1954  d_res_d_sig_trial(0,2) = (-2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(theta)*sin(beta))/(CG*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1955  d_res_d_sig_trial(1,0) = 0;
1956  d_res_d_sig_trial(1,1) = (2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) +
1957  (1 + fPsi)*sin(beta))*sin(theta))/(CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
1958  d_res_d_sig_trial(1,2) = (-2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*
1959  (1 + fPsi + 6*(-1 + fPsi)*sin(beta) - 2*(-1 + fPsi)*sin(3*beta))*sin(theta))/
1960  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
1961  d_res_d_sig_trial(2,0) = sqrt(3);
1962  d_res_d_sig_trial(2,1) = 0;
1963  d_res_d_sig_trial(2,2) = 0;
1964 
1965 
1966 
1967  // Derivative for the projected stresses respect to internal variables
1968 
1969  d_sig_proj_d_theta_beta_kappa(0,0) = ((fA - exp(fB*k)*fC)*(4*sqrt(3)*fPsi*cos(beta)*cos(theta) +
1970  fR*(1 + fPsi + (-1 + fPsi)*sin(3*beta))*sin(theta)))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1971  d_sig_proj_d_theta_beta_kappa(0,1) = (-4*(fA - exp(fB*k)*fC)*fPsi*(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) +
1972  (1 + fPsi)*sin(beta))*sin(theta))/(sqrt(3)*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
1973  d_sig_proj_d_theta_beta_kappa(0,2) = (1 + exp(fB*k)*fB*fC*fR*cos(theta) - (4*sqrt(3)*exp(fB*k)*fB*fC*fPsi*cos(beta)*sin(theta))/
1974  (1 + fPsi + (-1 + fPsi)*sin(3*beta)))/3.;
1975 
1976 
1977  d_sig_proj_d_theta_beta_kappa(1,0) = ((fA - exp(fB*k)*fC)*(-2*sqrt(3)*fPsi*cos(beta)*cos(theta) + 6*fPsi*cos(theta)*sin(beta) +
1978  fR*(1 + fPsi + (-1 + fPsi)*sin(3*beta))*sin(theta)))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1979  d_sig_proj_d_theta_beta_kappa(1,1) = (2*(fA - exp(fB*k)*fC)*fPsi*(3*(1 + fPsi)*cos(beta) + 2*sqrt(3)*(-1 + fPsi)*cos(2*beta) -
1980  sqrt(3)*cos(4*beta) + sqrt(3)*fPsi*cos(4*beta) + sqrt(3)*sin(beta) + sqrt(3)*fPsi*sin(beta) -
1981  6*sin(2*beta) + 6*fPsi*sin(2*beta) + 3*sin(4*beta) - 3*fPsi*sin(4*beta))*sin(theta))/
1982  (3.*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
1983  d_sig_proj_d_theta_beta_kappa(1,2) = (1 + fPsi + (-1 + fPsi)*sin(3*beta) + exp(fB*k)*fB*fC*fR*cos(theta)*
1984  (1 + fPsi + (-1 + fPsi)*sin(3*beta)) + 2*sqrt(3)*exp(fB*k)*fB*fC*fPsi*cos(beta)*sin(theta) -
1985  6*exp(fB*k)*fB*fC*fPsi*sin(beta)*sin(theta))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1986 
1987 
1988  d_sig_proj_d_theta_beta_kappa(2,0) = ((fA - exp(fB*k)*fC)*(-2*sqrt(3)*fPsi*cos(beta)*cos(theta) - 6*fPsi*cos(theta)*sin(beta) +
1989  fR*(1 + fPsi + (-1 + fPsi)*sin(3*beta))*sin(theta)))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1990  d_sig_proj_d_theta_beta_kappa(2,1) = (2*(fA - exp(fB*k)*fC)*fPsi*(-3*(1 + fPsi)*cos(beta) + 2*sqrt(3)*(-1 + fPsi)*cos(2*beta) -
1991  sqrt(3)*cos(4*beta) + sqrt(3)*fPsi*cos(4*beta) + sqrt(3)*sin(beta) + sqrt(3)*fPsi*sin(beta) +
1992  6*sin(2*beta) - 6*fPsi*sin(2*beta) - 3*sin(4*beta) + 3*fPsi*sin(4*beta))*sin(theta))/
1993  (3.*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
1994  d_sig_proj_d_theta_beta_kappa(2,2) = (1 + fPsi + (-1 + fPsi)*sin(3*beta) + exp(fB*k)*fB*fC*fR*cos(theta)*
1995  (1 + fPsi + (-1 + fPsi)*sin(3*beta)) + 2*sqrt(3)*exp(fB*k)*fB*fC*fPsi*cos(beta)*sin(theta) +
1996  6*exp(fB*k)*fB*fC*fPsi*sin(beta)*sin(theta))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
1997 
1998  // Composing the gradient
1999  jac_inv *= -1.0;
2000  jac_inv.Multiply(d_res_d_sig_trial, d_theta_beta_kappa_d_sig_trial);//(ok)
2001  d_sig_proj_d_theta_beta_kappa.Multiply(d_theta_beta_kappa_d_sig_trial, d_sig_proj_d_sig_trial);
2002  d_sig_proj_d_sig_trial.Multiply(Rot, *gradient);
2003 
2004 
2005 }
2006 
2008 void TPZSandlerExtended::ComputeCapVertexTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const{
2009 
2010  TPZFMatrix<STATE> Rot(3,3,0.0);
2011  TPZHWTools::GetRotMatrix(Rot);
2012 
2013  STATE fv = F(kproj);
2014  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
2015  STATE CG = fE/(2.0*(1.0 + fnu));
2016 
2017  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
2018  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
2019  TPZHWTools::FromPrincipalToHWCyl(projected_stress, hw_space_xi_rho_beta);
2020  TPZHWTools::FromPrincipalToHWCart(projected_stress, rhw_space_s1_s2_s3);
2021 
2022  STATE k = kproj;
2023  STATE theta = 0.0;
2024 
2025  // Computing the jacobian and the inverse
2026  STATE jac,jac_inv;
2027  Jacobianf2Vertex(trial_stress, k, jac); // Jacobian
2028  jac_inv = 1.0 / jac; // Jacobian inverse
2029 
2030  TPZFMatrix<STATE> d_res_d_sig_trial(1,3,0.0);
2031  TPZFMatrix<STATE> d_sig_proj_d_kappa(3,1,0.0);
2032  TPZFMatrix<STATE> d_kappa_d_sig_trial(1,3,0.0);
2033  TPZFMatrix<STATE> d_sig_proj_d_sig_trial(3,3,0.0);
2034 
2035  // Derivative for the Covertex residue respect to trial stresses (1x3)
2036  d_res_d_sig_trial(0,0) = sqrt(3);
2037  d_res_d_sig_trial(0,1) =0;
2038  d_res_d_sig_trial(0,2) =0;
2039 
2040  // Derivative for the projected stresses respect to internal variables beta and k (3x1)
2041  d_sig_proj_d_kappa(0,0) = (1 + exp(fB*k)*fB*fC*fR*cos(theta))/3.;
2042  d_sig_proj_d_kappa(1,0) = (1 + exp(fB*k)*fB*fC*fR*cos(theta))/3.;
2043  d_sig_proj_d_kappa(2,0) = (1 + exp(fB*k)*fB*fC*fR*cos(theta))/3.;
2044 
2045  // Composing the gradient
2046  jac_inv *= -1.0;
2047  d_kappa_d_sig_trial = jac_inv * d_res_d_sig_trial;
2048  d_sig_proj_d_kappa.Multiply(d_kappa_d_sig_trial, d_sig_proj_d_sig_trial);
2049  d_sig_proj_d_sig_trial.Multiply(Rot, *gradient);
2050 
2051 }
2052 
2054 void TPZSandlerExtended::ComputeCapCoVertexTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const{
2055 
2056  TPZFMatrix<STATE> Rot(3,3,0.0);
2057  TPZHWTools::GetRotMatrix(Rot);
2058 
2059  STATE fv = F(kproj);
2060  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
2061  STATE CG = fE/(2.0*(1.0 + fnu));
2062 
2063  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
2064  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);;
2065  TPZHWTools::FromPrincipalToHWCyl(projected_stress, hw_space_xi_rho_beta);
2066  TPZHWTools::FromPrincipalToHWCart(projected_stress, rhw_space_s1_s2_s3);
2067 
2068  STATE k = kproj;
2069  STATE i1v = sqrt(3)*hw_space_xi_rho_beta[0];
2070  STATE theta = M_PI_2;
2071  STATE beta = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
2072 
2073  // Computing the jacobian and the inverse
2074  TPZFNMatrix<4, STATE> jac(2,2,0.0),jac_inv(2,2,0.0);
2075  Jacobianf2CoVertex(trial_stress, beta, k, jac); // Jacobian
2076  TPZHWTools::A2x2Inverse(jac, jac_inv); // Jacobian inverse
2077 
2078  TPZFMatrix<STATE> d_res_d_sig_trial(2,3,0.0);
2079  TPZFMatrix<STATE> d_sig_proj_d_beta_kappa(3,2,0.0);
2080  TPZFMatrix<STATE> d_beta_kappa_d_sig_trial(2,3,0.0);
2081  TPZFMatrix<STATE> d_sig_proj_d_sig_trial(3,3,0.0);
2082 
2083  // Derivative for the Covertex residue respect to trial stresses (2x3)
2084  d_res_d_sig_trial(0,0) = 0;
2085  d_res_d_sig_trial(0,1) = (2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) +
2086  (1 + fPsi)*sin(beta)))/(CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2087  d_res_d_sig_trial(0,2) = (-2*sqrt(2)*(fA - exp(fB*k)*fC)*fPsi*cos(beta)*
2088  (1 + fPsi + 6*(-1 + fPsi)*sin(beta) - 2*(-1 + fPsi)*sin(3*beta)))/
2089  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2090 
2091  d_res_d_sig_trial(1,0) = sqrt(3);
2092  d_res_d_sig_trial(1,1) = 0;
2093  d_res_d_sig_trial(1,2) = 0;
2094 
2095  // Derivative for the projected stresses respect to internal variables beta and k (3x2)
2096 
2097  d_sig_proj_d_beta_kappa(0,0) = (-4*(fA - exp(fB*k)*fC)*fPsi*(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) +
2098  (1 + fPsi)*sin(beta)))/(sqrt(3)*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2099  d_sig_proj_d_beta_kappa(0,1) = (1.0/3.0) - (4*exp(fB*k)*fB*fC*fPsi*cos(beta))/
2100  (sqrt(3)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2101 
2102  d_sig_proj_d_beta_kappa(1,0) = (2*(fA - exp(fB*k)*fC)*fPsi*(3*(1 + fPsi)*cos(beta) + 2*sqrt(3)*(-1 + fPsi)*cos(2*beta) -
2103  sqrt(3)*cos(4*beta) + sqrt(3)*fPsi*cos(4*beta) + sqrt(3)*sin(beta) + sqrt(3)*fPsi*sin(beta) -
2104  6*sin(2*beta) + 6*fPsi*sin(2*beta) + 3*sin(4*beta) - 3*fPsi*sin(4*beta)))/
2105  (3.*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2106  d_sig_proj_d_beta_kappa(1,1) = (1 + fPsi + 2*sqrt(3)*exp(fB*k)*fB*fC*fPsi*cos(beta) - 6*exp(fB*k)*fB*fC*fPsi*sin(beta) -
2107  sin(3*beta) + fPsi*sin(3*beta))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2108 
2109 
2110  d_sig_proj_d_beta_kappa(2,0) = (2*(fA - exp(fB*k)*fC)*fPsi*(-3*(1 + fPsi)*cos(beta) + 2*sqrt(3)*(-1 + fPsi)*cos(2*beta) -
2111  sqrt(3)*cos(4*beta) + sqrt(3)*fPsi*cos(4*beta) + sqrt(3)*sin(beta) + sqrt(3)*fPsi*sin(beta) +
2112  6*sin(2*beta) - 6*fPsi*sin(2*beta) - 3*sin(4*beta) + 3*fPsi*sin(4*beta)))/
2113  (3.*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2114  d_sig_proj_d_beta_kappa(2,1) = (1 + fPsi + 2*sqrt(3)*exp(fB*k)*fB*fC*fPsi*cos(beta) + 6*exp(fB*k)*fB*fC*fPsi*sin(beta) -
2115  sin(3*beta) + fPsi*sin(3*beta))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2116 
2117 
2118  // Composing the gradient
2119  jac_inv *= -1.0;
2120  jac_inv.Multiply(d_res_d_sig_trial, d_beta_kappa_d_sig_trial);//(ok)
2121  d_sig_proj_d_beta_kappa.Multiply(d_beta_kappa_d_sig_trial, d_sig_proj_d_sig_trial);
2122  d_sig_proj_d_sig_trial.Multiply(Rot, *gradient);
2123 
2124 }
2125 
2127 void TPZSandlerExtended::ComputeFailureTangent(const TPZVec<STATE> &trial_stress, STATE kprev, TPZVec<STATE> &projected_stress, STATE &kproj, TPZFMatrix<REAL> * gradient) const{
2128 
2129  TPZFMatrix<STATE> Rot(3,3,0.0);
2130  TPZHWTools::GetRotMatrix(Rot);
2131 
2132 
2133  STATE fv = F(kproj);
2134  STATE CK = fE/(3.0*(1.0 - 2.0 *fnu));
2135  STATE CG = fE/(2.0*(1.0 + fnu));
2136 
2137  TPZManVector<STATE, 3> hw_space_xi_rho_beta(3);
2138  TPZManVector<STATE, 3> rhw_space_s1_s2_s3(3);
2139  TPZHWTools::FromPrincipalToHWCyl(projected_stress, hw_space_xi_rho_beta);
2140  TPZHWTools::FromPrincipalToHWCart(projected_stress, rhw_space_s1_s2_s3);
2141 
2142  STATE i1 = sqrt(3.0)*rhw_space_s1_s2_s3[0];
2143  STATE beta = atan2(rhw_space_s1_s2_s3[2],rhw_space_s1_s2_s3[1]);
2144  STATE k = kproj;
2145 
2146 
2147  // Computing the jacobian and the inverse
2148  TPZFMatrix<STATE> jac(3,3,0.0),jac_inv(3,3,0.0);
2149  Jacobianf1(trial_stress, i1, beta, k, jac); // Jacobian
2150  TPZHWTools::A3x3Inverse(jac, jac_inv); // Jacobian inverse
2151 
2152  TPZFMatrix<STATE> d_res_d_sig_trial(3,3,0.0);
2153  TPZFMatrix<STATE> d_sig_proj_d_i1_beta_kappa(3,3,0.0);
2154  TPZFMatrix<STATE> d_i1_beta_kappa_d_sig_trial(3,3,0.0);
2155  TPZFMatrix<STATE> d_sig_proj_d_sig_trial(3,3,0.0);
2156 
2157  // Derivative for the cap residue respect to trial stresses
2158  d_res_d_sig_trial(0,0) = -2/(3.*sqrt(3)*CK);
2159  d_res_d_sig_trial(0,1) = (2*sqrt(2)*exp(fB*i1)*fB*fC*fPsi*cos(beta))/(CG*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2160  d_res_d_sig_trial(0,2) = (2*sqrt(2)*exp(fB*i1)*fB*fC*fPsi*sin(beta))/(CG*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2161 
2162  d_res_d_sig_trial(1,0) = 0;
2163  d_res_d_sig_trial(1,1) = (2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) +
2164  (1 + fPsi)*sin(beta)))/(CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2165  d_res_d_sig_trial(1,2) = (-2*sqrt(2)*(fA - exp(fB*i1)*fC)*fPsi*cos(beta)*
2166  (1 + fPsi + 6*(-1 + fPsi)*sin(beta) - 2*(-1 + fPsi)*sin(3*beta)))/
2167  (CG*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2168 
2169  d_res_d_sig_trial(2,0) = sqrt(3);
2170  d_res_d_sig_trial(2,1) = 0;
2171  d_res_d_sig_trial(2,2) = 0;
2172 
2173  // Derivative for the projected stresses respect to internal variables
2174  d_sig_proj_d_i1_beta_kappa(0,0) = 1.0/3.0 - (4*exp(fB*i1)*fB*fC*fPsi*cos(beta))/(sqrt(3)*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2175  d_sig_proj_d_i1_beta_kappa(0,1) = (-4*(fA - exp(fB*i1)*fC)*fPsi*(2*(-1 + fPsi)*cos(2*beta) + (-1 + fPsi)*cos(4*beta) +
2176  (1 + fPsi)*sin(beta)))/(sqrt(3)*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2177  d_sig_proj_d_i1_beta_kappa(0,2) = 0;
2178 
2179  d_sig_proj_d_i1_beta_kappa(1,0) = (1 + fPsi + 2*sqrt(3)*exp(fB*i1)*fB*fC*fPsi*cos(beta) - 6*exp(fB*i1)*fB*fC*fPsi*sin(beta) -
2180  sin(3*beta) + fPsi*sin(3*beta))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2181  d_sig_proj_d_i1_beta_kappa(1,1) = (2*(fA - exp(fB*i1)*fC)*fPsi*(3*(1 + fPsi)*cos(beta) + 2*sqrt(3)*(-1 + fPsi)*cos(2*beta) -
2182  sqrt(3)*cos(4*beta) + sqrt(3)*fPsi*cos(4*beta) + sqrt(3)*sin(beta) + sqrt(3)*fPsi*sin(beta) -
2183  6*sin(2*beta) + 6*fPsi*sin(2*beta) + 3*sin(4*beta) - 3*fPsi*sin(4*beta)))/
2184  (3.*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2185  d_sig_proj_d_i1_beta_kappa(1,2) = 0;
2186 
2187  d_sig_proj_d_i1_beta_kappa(2,0) = (1 + fPsi + 2*sqrt(3)*exp(fB*i1)*fB*fC*fPsi*cos(beta) + 6*exp(fB*i1)*fB*fC*fPsi*sin(beta) -
2188  sin(3*beta) + fPsi*sin(3*beta))/(3.*(1 + fPsi + (-1 + fPsi)*sin(3*beta)));
2189  d_sig_proj_d_i1_beta_kappa(2,1) = (2*(fA - exp(fB*i1)*fC)*fPsi*(-3*(1 + fPsi)*cos(beta) + 2*sqrt(3)*(-1 + fPsi)*cos(2*beta) -
2190  sqrt(3)*cos(4*beta) + sqrt(3)*fPsi*cos(4*beta) + sqrt(3)*sin(beta) + sqrt(3)*fPsi*sin(beta) +
2191  6*sin(2*beta) - 6*fPsi*sin(2*beta) - 3*sin(4*beta) + 3*fPsi*sin(4*beta)))/
2192  (3.*pow(1 + fPsi + (-1 + fPsi)*sin(3*beta),2));
2193  d_sig_proj_d_i1_beta_kappa(2,2) = 0;
2194 
2195  // Composing the gradient
2196  jac_inv *= -1.0;
2197  jac_inv.Multiply(d_res_d_sig_trial, d_i1_beta_kappa_d_sig_trial);//(ok)
2198  d_sig_proj_d_i1_beta_kappa.Multiply(d_i1_beta_kappa_d_sig_trial, d_sig_proj_d_sig_trial);
2199  d_sig_proj_d_sig_trial.Multiply(Rot, *gradient);
2200 
2201 // jac.Print(std::cout);
2202 // jac_inv.Print(std::cout);
2203 // d_res_d_sig_trial.Print(std::cout);
2204 // d_i1_beta_kappa_d_sig_trial.Print(std::cout);
2205 // d_sig_proj_d_i1_beta_kappa.Print(std::cout);
2206 // d_sig_proj_d_sig_trial.Print(std::cout);
2207 // gradient->Print(std::cout);
2208 
2209 }
2210 
2211 void TPZSandlerExtended::ProjectSigmaDep(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &sigproj, STATE &kproj, TPZFMatrix<STATE> &GradSigma) const {
2212 
2213  DebugStop(); // Deprecated method.
2214 
2215  STATE I1;
2216  //Firstk(epspv,k0);
2217  TPZManVector<STATE, 2> yield(2);
2218  I1 = sigtrial[0] + sigtrial[1] + sigtrial[2];
2219 
2220  YieldFunction(sigtrial, kprev, yield);
2221  bool threeEigEqual = (fabs(sigtrial[0] - sigtrial[1]) < ftol && fabs(sigtrial[1] - sigtrial[2]) < ftol);
2222 
2223  // kprev corresponde a L do artigo
2224  if (I1 < kprev) {
2225  if (yield[1] > 0. && !threeEigEqual) {
2226 #ifdef LOG4CXX
2227  if (logger->isDebugEnabled()) {
2228  std::stringstream sout;
2229  sout << "Projecting on F2, distinct eigenvalues " << sigtrial;
2230  LOGPZ_DEBUG(logger, sout.str())
2231  }
2232 #endif
2233  ProjectF2(sigtrial, kprev, sigproj, kproj);
2234  // we can compute the tangent matrix
2235  TPZFNMatrix<9, STATE> dbetadsigtrial(3, 3), jacF2(3, 3), DF2cart(3, 3);
2236  STATE theta, beta;
2237  SurfaceParamF2(sigproj, kproj, theta, beta);
2238  GradF2SigmaTrial(sigtrial, theta, beta, kproj, kprev, dbetadsigtrial);
2239  Jacobianf2(sigtrial, theta, beta, kproj, jacF2);
2240  jacF2.Solve_LU(&dbetadsigtrial);
2241  DF2Cart(theta, beta, kproj, DF2cart);
2242  DF2cart.Multiply(dbetadsigtrial, GradSigma);
2243  GradSigma *= -1.;
2244  } else if (yield[1] > 0. && threeEigEqual) {
2245 #ifdef LOG4CXX
2246  if (logger->isDebugEnabled()) {
2247  std::stringstream sout;
2248  sout << "Projecting on F2, equal eigenvalues " << sigtrial;
2249  LOGPZ_DEBUG(logger, sout.str())
2250  }
2251 #endif
2252  ProjectBetaConstF2(sigtrial, kprev, sigproj, kproj);
2253  // we can compute the tangent matrix
2254  STATE theta, beta;
2255  // compute theta beta as a function of sigproj, kproj
2256  SurfaceParamF2(sigproj, kproj, theta, beta);
2257  // theta should be Pi
2258  // for hydrostatic stress beta doesn't mean anything
2259 #ifdef LOG4CXX
2260  if (logger->isDebugEnabled()) {
2261  std::stringstream sout;
2262  sout << "Surface parameters for sigproj = " << sigproj << " kproj " << kproj << " theta " << theta;
2263  LOGPZ_DEBUG(logger, sout.str())
2264  }
2265 #endif
2266  TPZManVector<STATE, 2> sigtrialIJ(2, 0.), sigprojIJ(2), ddistf2(2);
2267  sigtrialIJ[0] = sigtrial[0] + sigtrial[1] + sigtrial[2];
2268  DDistF2IJ(sigtrialIJ, theta, kproj, kprev, ddistf2);
2269 #ifdef LOG4CXX
2270  if (logger->isDebugEnabled()) {
2271  std::stringstream sout;
2272  sout << "Derivative of the distance function (should be zero) = " << ddistf2;
2273  LOGPZ_DEBUG(logger, sout.str())
2274  }
2275 #endif
2276  TPZFNMatrix<4, STATE> JacThetaK(2, 2), JacSigtrIJ(2, 2), JacSigprojThetaK(2, 2);
2277  {
2278  TPZManVector<TFad<2, STATE>, 2> sigtrialIJFAD(2), ddistf2fad(2);
2279  TFad<2, STATE> thetafad(theta, 0), kprojfad(kproj, 1);
2280  sigtrialIJFAD[0].val() = sigtrialIJ[0];
2281  sigtrialIJFAD[1].val() = sigtrialIJ[1];
2282  DDistF2IJ(sigtrialIJFAD, thetafad, kprojfad, kprev, ddistf2fad);
2283  JacThetaK(0, 0) = ddistf2fad[0].d(0);
2284  JacThetaK(0, 1) = ddistf2fad[0].d(1);
2285  JacThetaK(1, 0) = ddistf2fad[1].d(0);
2286  JacThetaK(1, 1) = ddistf2fad[1].d(1);
2287  }
2288  {
2289  TPZManVector<TFad<2, STATE>, 2> sigtrialIJFAD(2), ddistf2fad(2);
2290  TFad<2, STATE> thetafad(theta), kprojfad(kproj);
2291  sigtrialIJFAD[0].val() = sigtrialIJ[0];
2292  sigtrialIJFAD[0].fastAccessDx(0) = 1.;
2293  sigtrialIJFAD[1].val() = sigtrialIJ[1];
2294  sigtrialIJFAD[1].fastAccessDx(1) = 1.;
2295  DDistF2IJ(sigtrialIJFAD, thetafad, kprojfad, kprev, ddistf2fad);
2296  JacSigtrIJ(0, 0) = ddistf2fad[0].d(0);
2297  JacSigtrIJ(0, 1) = ddistf2fad[0].d(1);
2298  JacSigtrIJ(1, 0) = ddistf2fad[1].d(0);
2299  JacSigtrIJ(1, 1) = ddistf2fad[1].d(1);
2300  }
2301  FromThetaKToSigIJ(theta, kproj, sigprojIJ);
2302  {
2303  TPZManVector<TFad<2, STATE>, 2> sigprojIJFAD(2);
2304  TFad<2, STATE> thetafad(theta, 0), kprojfad(kproj, 1);
2305  FromThetaKToSigIJ(thetafad, kprojfad, sigprojIJFAD);
2306  JacSigprojThetaK(0, 0) = sigprojIJFAD[0].d(0);
2307  JacSigprojThetaK(0, 1) = sigprojIJFAD[0].d(1);
2308  JacSigprojThetaK(1, 0) = sigprojIJFAD[1].d(0);
2309  JacSigprojThetaK(1, 1) = sigprojIJFAD[1].d(1);
2310  }
2311 #ifdef LOG4CXX
2312  if (logger->isDebugEnabled()) {
2313  std::stringstream sout;
2314  sout << "Derivative of the distanceIJ " << ddistf2 << std::endl;
2315  JacThetaK.Print("Derivative of distanceIJ with respect to theta, K", sout);
2316  JacSigtrIJ.Print("Derivative of distanceIJ with respect to sigtialIJ", sout);
2317  sout << "SigmaProjected IJ " << sigprojIJ << std::endl;
2318  JacSigprojThetaK.Print("Derivative of sigproj with respect to theta K", sout);
2319  LOGPZ_DEBUG(logger, sout.str())
2320  }
2321 #endif
2322  std::list<int64_t> singular;
2323  JacThetaK.Solve_LU(&JacSigtrIJ, singular);
2324 #ifdef LOG4CXX
2325  if (logger->isDebugEnabled()) {
2326  std::stringstream sout;
2327  sout << "Negative of derivative of Theta,K with respect to sigtrIJ" << std::endl;
2328  JacSigtrIJ.Print("Derivative = ", sout);
2329  LOGPZ_DEBUG(logger, sout.str())
2330  }
2331 #endif
2332  TPZFNMatrix<4, STATE> dsigprojdsigtr(2, 2);
2333  JacSigprojThetaK.Multiply(JacSigtrIJ, dsigprojdsigtr);
2334  dsigprojdsigtr *= -1.;
2335 #ifdef LOG4CXX
2336  if (logger->isDebugEnabled()) {
2337  std::stringstream sout;
2338  dsigprojdsigtr.Print("Derivative of SigprojIJ with respect to SigtrialIJ", sout);
2339  LOGPZ_DEBUG(logger, sout.str())
2340  }
2341 #endif
2342  for (int i = 0; i < 3; i++) {
2343  for (int j = 0; j < 3; j++) {
2344  GradSigma(i, j) = dsigprojdsigtr(0, 0) / 3.;
2345  if (i == j) {
2346  GradSigma(i, j) += dsigprojdsigtr(1, 1)*2. / 3.;
2347  } else {
2348  GradSigma(i, j) -= dsigprojdsigtr(1, 1) / 3.;
2349  }
2350  }
2351  }
2352  } else {
2353  sigproj = sigtrial;
2354  kproj = kprev;
2355  GradSigma.Identity();
2356 
2357  }
2358  } else {
2359  if (yield[0] > 0.) {
2360 #ifdef LOG4CXX
2361  if (logger->isDebugEnabled()) {
2362  std::stringstream sout;
2363  sout << "Projecting on F1";
2364  LOGPZ_DEBUG(logger, sout.str())
2365  }
2366 #endif
2367  ProjectF1(sigtrial, kprev, sigproj, kproj);
2368 #ifdef PZDEBUG
2369  {
2370  TPZManVector<STATE, 2> yieldcopy(2);
2371  YieldFunction(sigproj, kproj, yieldcopy);
2372  if (fabs(yieldcopy[0]) > 1.e-3) {
2373  DebugStop();
2374  }
2375  }
2376 #endif
2377 
2378  I1 = 0.;
2379  for (int i = 0; i < 3; i++) {
2380  I1 += sigproj[i];
2381  }
2382  if (I1 < kproj) {
2383  ProjectRing(sigtrial, kprev, sigproj, kproj);
2384 
2385 
2386  // we can compute the tangent matrix
2387  TPZFNMatrix<9, STATE> dbetadsigtrial(3, 3), jacF2(3, 3), DF2cart(3, 3);
2388  STATE theta, beta;
2389  SurfaceParamF2(sigproj, kproj, theta, beta);
2390 #ifdef PZDEBUG
2391  if (fabs(theta) - M_PI_2 > 1.e-8) {
2392  DebugStop();
2393  }
2394 #endif
2395  GradF2SigmaTrial(sigtrial, theta, beta, kproj, kprev, dbetadsigtrial);
2396  for (int i = 0; i < 3; i++) dbetadsigtrial(0, i) = 0.;
2397  Jacobianf2(sigtrial, theta, beta, kproj, jacF2);
2398  for (int i = 0; i < 3; i++) {
2399  jacF2(i, 0) = 0.;
2400  jacF2(0, 1) = 0.;
2401  }
2402  jacF2(0, 0) = 1.;
2403  jacF2.Solve_LU(&dbetadsigtrial);
2404  DF2Cart(theta, beta, kproj, DF2cart);
2405  for (int i = 0; i < 3; i++) DF2cart(i, 0) = 0.;
2406  DF2cart.Multiply(dbetadsigtrial, GradSigma);
2407  GradSigma *= -1.;
2408 
2409  } else {
2410  // we can compute the tangent matrix
2411  TPZFNMatrix<9, STATE> dbetadsigtrial(2, 3), jacF1(2, 2), DF1cart(3, 2);
2412  STATE xi, beta;
2413  SurfaceParamF1(sigproj, xi, beta);
2414  GradF1SigmaTrial(sigtrial, xi, beta, dbetadsigtrial);
2415  D2DistFunc1(sigtrial, xi, beta, jacF1);
2416  jacF1.Solve_LU(&dbetadsigtrial);
2417  DF1Cart(xi, beta, DF1cart);
2418  DF1cart.Multiply(dbetadsigtrial, GradSigma);
2419  GradSigma *= -1.;
2420 
2421  }
2422  } else {
2423 #ifdef LOG4CXX
2424  {
2425  if (logger->isDebugEnabled()) {
2426  std::stringstream sout;
2427  sout << "Elastic Behaviour";
2428  LOGPZ_DEBUG(logger, sout.str())
2429  }
2430  }
2431 #endif
2432  // elastic behaviour
2433  sigproj = sigtrial;
2434  kproj = kprev;
2435  GradSigma.Identity();
2436  }
2437  }
2438 }
2439 
2440 
2441 #define Cos cos
2442 #define Sin sin
2443 #define Sqrt sqrt
2444 
2445 
2447 
2448 void TPZSandlerExtended::GradF1SigmaTrial(const TPZVec<STATE> &sigtrial, STATE xi, STATE beta, TPZFMatrix<STATE> &deriv) const {
2449  STATE s3beta = sin(3. * beta);
2450  STATE c3beta = cos(3. * beta);
2451  STATE Gamma = (1 + s3beta)+(1. - s3beta) / fPsi;
2452  STATE DGamma = 3. * c3beta * (1. - 1. / fPsi);
2453  STATE SQR3 = sqrt(3.);
2454  STATE FFI = fA - fPhi * SQR3 * xi - fC * exp(fB * SQR3 * xi);
2455  STATE NN = fN;
2456  deriv.Redim(2, 3);
2457  deriv(0, 0) = (-2 * (fG * Gamma - 6 * fK * (SQR3 * fA * fB - SQR3 * fB * FFI + SQR3 * fPhi - 3 * fB * fPhi * xi) * Cos(beta))) / (3. * SQR3 * fG * fK * Gamma);
2458  deriv(1, 0) = (4 * (FFI - NN)*(DGamma * Cos(beta) + Gamma * Sin(beta))) / (SQR3 * fG * Gamma * Gamma);
2459 
2460  deriv(0, 1) = (-2 * (SQR3 * fG * Gamma + 9 * fK * (fA * fB + fPhi - fB * (FFI + SQR3 * fPhi * xi)) * Cos(beta) -
2461  9 * fK * (SQR3 * fA * fB + SQR3 * fPhi - fB * (SQR3 * FFI + 3 * fPhi * xi)) * Sin(beta))) / (9. * fG * fK * Gamma);
2462  deriv(1, 1) = (-2 * (FFI - NN)*((SQR3 * DGamma + 3 * Gamma) * Cos(beta) + (-3 * DGamma + SQR3 * Gamma) * Sin(beta))) / (3. * fG * Gamma * Gamma);
2463 
2464  deriv(0, 2) = (-2 * (SQR3 * fG * Gamma + 9 * fK * (fA * fB + fPhi - fB * (FFI + SQR3 * fPhi * xi)) * Cos(beta) +
2465  9 * fK * (SQR3 * fA * fB + SQR3 * fPhi - fB * (SQR3 * FFI + 3 * fPhi * xi)) * Sin(beta))) / (9. * fG * fK * Gamma);
2466  deriv(1, 2) = (-2 * (FFI - NN)*((SQR3 * DGamma - 3 * Gamma) * Cos(beta) + (3 * DGamma + SQR3 * Gamma) * Sin(beta))) / (3. * fG * Gamma * Gamma);
2467 }
2468 
2470 
2471 void TPZSandlerExtended::GradF2SigmaTrial(const TPZVec<STATE> &sigtrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZFMatrix<STATE> &deriv) const {
2472  STATE s3beta = sin(3. * beta);
2473  STATE c3beta = cos(3. * beta);
2474  STATE Gamma = (1 + s3beta)+(1. - s3beta) / fPsi;
2475  STATE DGamma = 3. * c3beta * (1. - 1. / fPsi);
2476  STATE SQR3 = sqrt(3.);
2477  STATE FFK = fA - fPhi * k - fC * exp(fB * k);
2478  deriv.Redim(3, 3);
2479  deriv(0, 0) = (-4 * (FFK - fN) * Cos(beta) * Cos(theta)) / (SQR3 * fG * Gamma) + (2 * FFK * fR * Sin(theta)) / (9. * fK);
2480  deriv(1, 0) = (4 * (FFK - fN)*(DGamma * Cos(beta) + Gamma * Sin(beta)) * Sin(theta)) / (SQR3 * fG * Gamma * Gamma);
2481  deriv(2, 0) = -1.;
2482 
2483  deriv(0, 1) = (2 * (3 * SQR3 * fK * (FFK - fN) * Cos(beta) * Cos(theta) - 9 * fK * (FFK - fN) * Cos(theta) * Sin(beta) + FFK * fG * fR * Gamma * Sin(theta))) / (9. * fG * fK * Gamma);
2484  deriv(1, 1) = (-2 * (FFK - fN)*((SQR3 * DGamma + 3 * Gamma) * Cos(beta) + (-3 * DGamma + SQR3 * Gamma) * Sin(beta)) * Sin(theta)) / (3. * fG * Gamma * Gamma);
2485  deriv(2, 1) = -1.;
2486 
2487  deriv(0, 2) = (2 * (3 * SQR3 * fK * (FFK - fN) * Cos(beta) * Cos(theta) + 9 * fK * (FFK - fN) * Cos(theta) * Sin(beta) + FFK * fG * fR * Gamma * Sin(theta))) / (9. * fG * fK * Gamma);
2488  deriv(1, 2) = (-2 * (FFK - fN)*((SQR3 * DGamma - 3 * Gamma) * Cos(beta) + (3 * DGamma + SQR3 * Gamma) * Sin(beta)) * Sin(theta)) / (3. * fG * Gamma * Gamma);
2489  deriv(2, 2) = -1.;
2490 
2491 }
2492 
2493 void TPZSandlerExtended::TaylorCheckDistF1(const TPZVec<STATE> &sigmatrial, STATE xi, STATE beta, TPZVec<STATE> &xnorm,
2494  TPZVec<STATE> &errnorm) const {
2495  STATE deltaxi = 0.4;
2496  STATE deltabeta = 0.05;
2497  STATE dist0 = DistF1(sigmatrial, xi, beta);
2498  TPZFNMatrix<4, STATE> jac(2, 1);
2499  DDistFunc1(sigmatrial, xi, beta, jac);
2500  xnorm.resize(10);
2501  errnorm.resize(10);
2502  for (int i = 1; i <= 10; i++) {
2503  STATE diffxi = deltaxi * i / 10.;
2504  STATE xinext = xi + diffxi;
2505  STATE diffbeta = deltabeta * i / 10.;
2506  STATE betanext = beta + diffbeta;
2507  STATE distnext = DistF1(sigmatrial, xinext, betanext);
2508  STATE distguess = dist0 + jac(0, 0) * diffxi + jac(1, 0) * diffbeta;
2509  xnorm[i - 1] = sqrt(diffxi * diffxi + diffbeta * diffbeta);
2510  errnorm[i - 1] = fabs(distnext - distguess);
2511  }
2512 
2513 }
2514 
2515 void TPZSandlerExtended::TaylorCheckDDistF1(const TPZVec<STATE> &sigmatrial, STATE xi, STATE beta, TPZVec<STATE> &xnorm,
2516  TPZVec<STATE> &errnorm) const {
2517  STATE deltaxi = 0.4;
2518  STATE deltabeta = 0.05;
2519  TPZFNMatrix<2, STATE> res0(2, 1), resid(2, 1), residguess(2, 1), diff(2, 1);
2520  TPZFNMatrix<4, STATE> jac(2, 2);
2521  DDistFunc1(sigmatrial, xi, beta, res0);
2522  D2DistFunc1(sigmatrial, xi, beta, jac);
2523  xnorm.resize(10);
2524  errnorm.resize(10);
2525  for (int i = 1; i <= 10; i++) {
2526  STATE diffxi = deltaxi * i / 10.;
2527  STATE xinext = xi + diffxi;
2528  STATE diffbeta = deltabeta * i / 10.;
2529  diff(0) = diffxi;
2530  diff(1) = diffbeta;
2531  STATE betanext = beta + diffbeta;
2532  DDistFunc1(sigmatrial, xinext, betanext, resid);
2533  jac.Multiply(diff, residguess);
2534  residguess += res0;
2535  xnorm[i - 1] = Norm(diff);
2536  errnorm[i - 1] = Norm(resid - residguess);
2537  }
2538 
2539 }
2540 
2541 void TPZSandlerExtended::TaylorCheckDDistF1DSigtrial(const TPZVec<STATE> &sigmatrial, STATE xi, STATE beta, TPZVec<STATE> &xnorm,
2542  TPZVec<STATE> &errnorm) const {
2543  TPZManVector<STATE, 3> deltasigma(3, 0.3);
2544  deltasigma[1] *= -1.;
2545  deltasigma[2] *= 0.7;
2546 
2547  TPZFNMatrix<2, STATE> res0(2, 1), resid(2, 1), residguess(2, 1), diff(3, 1);
2548  TPZFNMatrix<6, STATE> jac(2, 3);
2549  DDistFunc1(sigmatrial, xi, beta, res0);
2550  GradF1SigmaTrial(sigmatrial, xi, beta, jac);
2551  xnorm.resize(10);
2552  errnorm.resize(10);
2553  for (int i = 1; i <= 10; i++) {
2554  for (int j = 0; j < 3; j++) diff(j) = deltasigma[j] * i / 10.;
2555  TPZManVector<STATE, 3> sigmanext(3);
2556  for (int j = 0; j < 3; j++) sigmanext[j] = sigmatrial[j] + diff(j);
2557  DDistFunc1(sigmanext, xi, beta, resid);
2558  jac.Multiply(diff, residguess);
2559  residguess += res0;
2560  xnorm[i - 1] = Norm(diff);
2561  errnorm[i - 1] = Norm(resid - residguess);
2562  }
2563 
2564 }
2565 
2566 void TPZSandlerExtended::ConvergenceRate(TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm, TPZVec<STATE> &convergence) {
2567  convergence.resize(xnorm.size() - 1);
2568  for (int i = 1; i < xnorm.size(); i++) {
2569  convergence[i - 1] = log(errnorm[i] / errnorm[i - 1]) / log(xnorm[i] / xnorm[i - 1]);
2570  }
2571 }
2572 
2573 void TPZSandlerExtended::TaylorCheckDistF2(const TPZVec<STATE> &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &xnorm,
2574  TPZVec<STATE> &errnorm) const {
2575  STATE deltatheta = 0.4;
2576  STATE deltabeta = 0.05;
2577  STATE deltak = 0.;
2578  STATE dist0 = DistF2(sigmatrial, theta, beta, k);
2579  TPZFNMatrix<4, STATE> jac(3, 1);
2580  TPZManVector<STATE> fxnvec(3);
2581  Res2(sigmatrial, theta, beta, k, kprev, fxnvec);
2582  for (int kk = 0; kk < 3; kk++) jac(kk, 0) = fxnvec[kk];
2583  // DDistFunc2(sigmatrial, theta, beta, k, kprev, jac);
2584  xnorm.resize(10);
2585  errnorm.resize(10);
2586  for (int i = 1; i <= 10; i++) {
2587  STATE difftheta = deltatheta * i / 10.;
2588  STATE thetanext = theta + difftheta;
2589  STATE diffbeta = deltabeta * i / 10.;
2590  STATE betanext = beta + diffbeta;
2591  STATE diffk = deltak * i / 10.;
2592  STATE knext = k + deltak;
2593  STATE distnext = DistF2(sigmatrial, thetanext, betanext, knext);
2594  STATE distguess = dist0 + jac(0, 0) * difftheta + jac(1, 0) * diffbeta + jac(2, 0) * diffk;
2595  xnorm[i - 1] = sqrt(difftheta * difftheta + diffbeta * diffbeta + diffk * diffk);
2596  errnorm[i - 1] = fabs(distnext - distguess);
2597  }
2598 
2599 }
2600 
2601 void TPZSandlerExtended::TaylorCheckDDistF2(const TPZVec<STATE> &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &xnorm,
2602  TPZVec<STATE> &errnorm) const {
2603  STATE deltatheta = 0.4;
2604  STATE deltabeta = 0.05;
2605  STATE deltak = 0.5;
2606  TPZFNMatrix<3, STATE> res0(3, 1), resid(3, 1), residguess(3, 1), diff(3, 1);
2607  TPZFNMatrix<9, STATE> jac(3, 3);
2608  TPZManVector<STATE> fxnvec(3);
2609  Res2(sigmatrial, theta, beta, k, kprev, fxnvec);
2610  for (int kk = 0; kk < 3; kk++) res0(kk, 0) = fxnvec[kk];
2611  // DDistFunc2(sigmatrial, theta, beta, k, kprev, res0);
2612  Jacobianf2(sigmatrial, theta, beta, k, jac);
2613  xnorm.resize(10);
2614  errnorm.resize(10);
2615  for (int i = 1; i <= 10; i++) {
2616  STATE difftheta = deltatheta * i / 10.;
2617  STATE thetanext = theta + difftheta;
2618  STATE diffbeta = deltabeta * i / 10.;
2619  STATE betanext = beta + diffbeta;
2620  STATE diffk = deltak * i / 10.;
2621  STATE knext = k + diffk;
2622  diff(0) = difftheta;
2623  diff(1) = diffbeta;
2624  diff(2) = diffk;
2625  TPZManVector<STATE> fxnvec(3);
2626  Res2(sigmatrial, thetanext, betanext, knext, kprev, fxnvec);
2627  for (int k = 0; k < 3; k++) resid(k, 0) = fxnvec[k];
2628  // DDistFunc2(sigmatrial, thetanext, betanext, knext, kprev, resid);
2629  jac.Multiply(diff, residguess);
2630  residguess += res0;
2631  xnorm[i - 1] = Norm(diff);
2632  errnorm[i - 1] = Norm(resid - residguess);
2633  }
2634 
2635 }
2636 
2638 
2639 void TPZSandlerExtended::TaylorCheckDDistF2DSigtrial(const TPZVec<STATE> &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec<STATE> &xnorm,
2640  TPZVec<STATE> &errnorm) const {
2641  TPZManVector<STATE, 3> deltasigma(3, 0.3);
2642  deltasigma[1] *= -1.;
2643  deltasigma[2] *= 0.7;
2644 
2645  TPZFNMatrix<3, STATE> res0(3, 1), resid(3, 1), residguess(3, 1), diff(3, 1);
2646  TPZFNMatrix<9, STATE> jac(3, 3);
2647  TPZManVector<STATE> fxnvec(3);
2648  Res2(sigmatrial, theta, beta, k, kprev, fxnvec);
2649  for (int kk = 0; kk < 3; kk++) res0(kk, 0) = fxnvec[kk];
2650  // DDistFunc2(sigmatrial, theta, beta, k, kprev, res0);
2651  GradF2SigmaTrial(sigmatrial, theta, beta, k, kprev, jac);
2652  xnorm.resize(10);
2653  errnorm.resize(10);
2654  for (int i = 1; i <= 10; i++) {
2655  for (int j = 0; j < 3; j++) diff(j) = deltasigma[j] * i / 10.;
2656  TPZManVector<STATE, 3> sigmanext(3);
2657  for (int j = 0; j < 3; j++) sigmanext[j] = sigmatrial[j] + diff(j);
2658  TPZManVector<STATE> fxnvec(3);
2659  Res2(sigmanext, theta, beta, k, kprev, fxnvec);
2660  for (int k = 0; k < 3; k++) resid(k, 0) = fxnvec[k];
2661  // DDistFunc2(sigmanext, theta, beta,k,kprev,resid);
2662  jac.Multiply(diff, residguess);
2663  residguess += res0;
2664  xnorm[i - 1] = Norm(diff);
2665  errnorm[i - 1] = Norm(resid - residguess);
2666  }
2667 
2668 }
2669 
2670 #include "pzvec_extras.h"
2671 
2672 void TPZSandlerExtended::CheckCoordinateTransformation(TPZVec<STATE> &cart) {
2673  TPZManVector<STATE, 3> HWCart(3), HWCyl(3), HWCart2(3), Cart2(3);
2674  TPZHWTools::FromPrincipalToHWCart(cart, HWCart);
2675  TPZHWTools::FromPrincipalToHWCyl(cart, HWCyl);
2676  TPZHWTools::FromHWCylToHWCart(HWCyl, HWCart2);
2677  TPZHWTools::FromHWCylToPrincipal(HWCyl, Cart2);
2678  REAL dist1 = dist(cart, Cart2);
2679  REAL dist2 = dist(HWCart, HWCart2);
2680  cout << __FUNCTION__ << " dist1 = " << dist1 << " dist2 = " << dist2 << endl;
2681 
2682 }
2683 
2685 
2686 void TPZSandlerExtended::DF1Cart(STATE xi, STATE beta, TPZFMatrix<STATE> &DF1) const {
2687  STATE s3beta = sin(3. * beta);
2688  STATE c3beta = cos(3. * beta);
2689  STATE Gamma = (1 + s3beta)+(1. - s3beta) / fPsi;
2690  STATE DGamma = 3. * c3beta * (1. - 1. / fPsi);
2691  STATE SQR3 = sqrt(3.);
2692  STATE FFI = fA - fPhi * SQR3 * xi - fC * exp(fB * SQR3 * xi);
2693  DF1.Resize(3, 2);
2694  DF1(0, 0) = 1 / SQR3 + (4 * (-(SQR3 * fPhi) - SQR3 * fB * (fA - FFI - SQR3 * fPhi * xi)) * Cos(beta)) /
2695  (SQR3 * Gamma);
2696  DF1(1, 0) = 1 / SQR3 + (4 * (-(SQR3 * fPhi) -
2697  SQR3 * fB * (fA - FFI - SQR3 * fPhi * xi)) * Sin(beta - M_PI / 6.)) / (SQR3 * Gamma);
2698  DF1(2, 0) =
2699  1 / SQR3 - (4 * (-(SQR3 * fPhi) - SQR3 * fB * (fA - FFI - SQR3 * fPhi * xi)) *
2700  Sin(beta + M_PI / 6.)) / (SQR3 * Gamma);
2701 
2702  DF1(0, 1) = (-4 * DGamma * (FFI - fN) * Cos(beta)) / (SQR3 * Gamma * Gamma) -
2703  (4 * (FFI - fN) * Sin(beta)) / (SQR3 * Gamma);
2704  DF1(1, 1) = (4 * (FFI - fN) * Cos(beta - M_PI / 6.)) / (SQR3 * Gamma) -
2705  (4 * DGamma * (FFI - fN) * Sin(beta - M_PI / 6.)) / (SQR3 * Gamma * Gamma);
2706 
2707  DF1(2, 1) = (-4 * (FFI - fN) * Cos(beta + M_PI / 6.)) / (SQR3 * Gamma) +
2708  (4 * DGamma * (FFI - fN) * Sin(beta + M_PI / 6.)) / (SQR3 * Gamma * Gamma);
2709 }
2710 
2712 
2713 void TPZSandlerExtended::DF2Cart(STATE theta, STATE beta, STATE k, TPZFMatrix<STATE> &DF2) const {
2714  STATE s3beta = sin(3. * beta);
2715  STATE c3beta = cos(3. * beta);
2716  STATE Gamma = (1 + s3beta)+(1. - s3beta) / fPsi;
2717  STATE DGamma = 3. * c3beta * (1. - 1. / fPsi);
2718  STATE SQR3 = sqrt(3.);
2719  STATE FFK = fA - fPhi * k - fC * exp(fB * k);
2720  DF2.Resize(3, 3);
2721  DF2(0, 0) = (4 * (FFK - fN) * Cos(beta) * Cos(theta)) / (SQR3 * Gamma) - (FFK * fR * Sin(theta)) / 3.;
2722  DF2(1, 0) = (4 * (FFK - fN) * Cos(theta) * Sin(beta - M_PI / 6.)) / (SQR3 * Gamma) - (FFK * fR * Sin(theta)) / 3.;
2723  DF2(2, 0) = (4 * (-FFK + fN) * Cos(theta) * Sin(beta + M_PI / 6.)) / (SQR3 * Gamma) - (FFK * fR * Sin(theta)) / 3.;
2724 
2725  DF2(0, 1) = (-4 * (FFK - fN)*(DGamma * Cos(beta) + Gamma * Sin(beta)) * Sin(theta)) / (SQR3 * Gamma * Gamma);
2726  DF2(1, 1) = (4 * (FFK - fN)*(Gamma * Cos(beta - M_PI / 6.) - DGamma * Sin(beta - M_PI / 6.)) * Sin(theta)) / (SQR3 * Gamma * Gamma);
2727  DF2(2, 1) = (-4 * (FFK - fN)*(Gamma * Cos(beta + M_PI / 6.) - DGamma * Sin(beta + M_PI / 6.)) * Sin(theta)) / (SQR3 * Gamma * Gamma);
2728 
2729  DF2(0, 2) = (Gamma - (fA * fB + fPhi - fB * (FFK + fPhi * k))*
2730  (fR * Gamma * Cos(theta) + 4 * SQR3 * Cos(beta) * Sin(theta))) / (3. * Gamma);
2731  DF2(1, 2) = (Gamma - (fA * fB + fPhi - fB * (FFK + fPhi * k))*
2732  (fR * Gamma * Cos(theta) + 4 * SQR3 * Sin(beta - M_PI / 6.) * Sin(theta))) / (3. * Gamma);
2733  DF2(2, 2) = (Gamma - (fA * fB + fPhi - fB * (FFK + fPhi * k))*
2734  (fR * Gamma * Cos(theta) - 4 * SQR3 * Sin(beta + M_PI / 6.) * Sin(theta))) / (3. * Gamma);
2735 }
2736 
2737 void TPZSandlerExtended::TaylorCheckDF1Cart(STATE xi, STATE beta, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2738  STATE deltaxi = 0.4;
2739  STATE deltabeta = 0.05;
2740  TPZFNMatrix<3, STATE> res0(3, 1), resid(3, 1), residguess(3, 1), diff(2, 1);
2741  TPZFNMatrix<9, STATE> jac(3, 2);
2742  TPZManVector<STATE> sigHWCyl(3), sigCart(3);
2743  F1Cyl(xi, beta, sigHWCyl);
2744  TPZHWTools::FromHWCylToPrincipal(sigHWCyl, sigCart);
2745  for (int i = 0; i < 3; i++) {
2746  res0(i) = sigCart[i];
2747  }
2748  DF1Cart(xi, beta, jac);
2749  xnorm.resize(10);
2750  errnorm.resize(10);
2751  for (int i = 1; i <= 10; i++) {
2752  STATE diffxi = deltaxi * i / 10.;
2753  STATE xinext = xi + diffxi;
2754  STATE diffbeta = deltabeta * i / 10.;
2755  diff(0) = diffxi;
2756  diff(1) = diffbeta;
2757  STATE betanext = beta + diffbeta;
2758  F1Cyl(xinext, betanext, sigHWCyl);
2759  TPZHWTools::FromHWCylToPrincipal(sigHWCyl, sigCart);
2760  for (int ii = 0; ii < 3; ii++) {
2761  resid(ii) = sigCart[ii];
2762  }
2763  jac.Multiply(diff, residguess);
2764  residguess += res0;
2765  xnorm[i - 1] = Norm(diff);
2766  errnorm[i - 1] = Norm(resid - residguess);
2767  }
2768 }
2769 
2770 void TPZSandlerExtended::TaylorCheckDF2Cart(STATE theta, STATE beta, STATE k, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2771  STATE deltatheta = 0.4;
2772  STATE deltabeta = 0.05;
2773  STATE deltak = 0.5;
2774  TPZFNMatrix<3, STATE> res0(3, 1), resid(3, 1), residguess(3, 1), diff(3, 1);
2775  TPZFNMatrix<9, STATE> jac(3, 3);
2776  DF2Cart(theta, beta, k, jac);
2777  TPZManVector<STATE> sigHWCyl(3), sigCart(3);
2778  F2Cyl(theta, beta, k, sigHWCyl);
2779  TPZHWTools::FromHWCylToPrincipal(sigHWCyl, sigCart);
2780  for (int i = 0; i < 3; i++) {
2781  res0(i) = sigCart[i];
2782  }
2783  xnorm.resize(10);
2784  errnorm.resize(10);
2785  for (int i = 1; i <= 10; i++) {
2786  STATE difftheta = deltatheta * i / 10.;
2787  STATE thetanext = theta + difftheta;
2788  STATE diffbeta = deltabeta * i / 10.;
2789  STATE betanext = beta + diffbeta;
2790  STATE diffk = deltak * i / 10.;
2791  STATE knext = k + diffk;
2792  diff(0) = difftheta;
2793  diff(1) = diffbeta;
2794  diff(2) = diffk;
2795  F2Cyl(thetanext, betanext, knext, sigHWCyl);
2796  TPZHWTools::FromHWCylToPrincipal(sigHWCyl, sigCart);
2797  for (int ii = 0; ii < 3; ii++) {
2798  resid(ii) = sigCart[ii];
2799  }
2800  jac.Multiply(diff, residguess);
2801  residguess += res0;
2802  xnorm[i - 1] = Norm(diff);
2803  errnorm[i - 1] = Norm(resid - residguess);
2804  }
2805 
2806 }
2807 
2808 void TPZSandlerExtended::TaylorCheckProjectSigma(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2809  TPZManVector<STATE, 3> deltasigma(3, -0.000012), sigproj(3);
2810  deltasigma[1] = -4.e-6;
2811  deltasigma[2] = -4.e-6;
2812  // TPZManVector<STATE,3> deltasigma(3,-0.01), sigproj(3);
2813  // deltasigma[1] *= 3;
2814  // deltasigma[2] *= -2;
2815 
2816  TPZFNMatrix<3, STATE> res0(3, 1), diff(3, 1), resid(3, 1), residguess(3, 1);
2817  TPZFNMatrix<9, STATE> jac(3, 3);
2818  STATE kproj;
2819  int m_type;
2820  ProjectSigmaDep(sigtrial, kprev, sigproj, kproj, jac);
2821  for (int j = 0; j < 3; j++) res0(j) = sigproj[j];
2822  xnorm.resize(10);
2823  errnorm.resize(10);
2824  for (int i = 1; i <= 10; i++) {
2825  TPZManVector<STATE, 3> diffsigma(3), nextsigma(3);
2826  for (int j = 0; j < 3; j++) {
2827  diffsigma[j] = deltasigma[j] * i / 10.;
2828  nextsigma[j] = sigtrial[j] + diffsigma[j];
2829  diff(j) = diffsigma[j];
2830  }
2831  jac.Multiply(diff, residguess);
2832  residguess += res0;
2833  ProjectSigma(nextsigma, kprev, sigproj, kproj, m_type);
2834  for (int j = 0; j < 3; j++) resid(j) = sigproj[j];
2835  xnorm[i - 1] = Norm(diff);
2836  errnorm[i - 1] = Norm(resid - residguess);
2837  }
2838 }
2839 
2841 
2842 void TPZSandlerExtended::TaylorCheckParamF1Sigtrial(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2843  TPZManVector<STATE, 3> deltasigma(3, -0.01), sigproj(3);
2844  deltasigma[1] *= 3.;
2845  deltasigma[2] *= -2.;
2846  TPZFNMatrix<3, STATE> res0(2, 1), diff(3, 1), resid(2, 1), residguess(2, 1);
2847  TPZFNMatrix<9, STATE> jac(2, 3), jacF1(2, 2), gradF1(2, 3);
2848  STATE kproj;
2849  ProjectF1(sigtrial, kprev, sigproj, kproj);
2850  STATE xi, beta;
2851  SurfaceParamF1(sigproj, xi, beta);
2852  res0(0) = xi;
2853  res0(1) = beta;
2854  D2DistFunc1(sigtrial, xi, beta, jacF1);
2855  GradF1SigmaTrial(sigtrial, xi, beta, gradF1);
2856  //gradF1.Transpose();
2857  jacF1.Solve_LU(&gradF1);
2858  jac = -1. * gradF1;
2859  xnorm.resize(10);
2860  errnorm.resize(10);
2861  for (int i = 1; i <= 10; i++) {
2862  TPZManVector<STATE, 3> diffsigma(3), nextsigma(3);
2863  for (int j = 0; j < 3; j++) {
2864  diffsigma[j] = deltasigma[j] * i / 10.;
2865  nextsigma[j] = sigtrial[j] + diffsigma[j];
2866  diff(j) = diffsigma[j];
2867  }
2868  jac.Multiply(diff, residguess);
2869  residguess += res0;
2870  ProjectF1(nextsigma, kprev, sigproj, kproj);
2871  SurfaceParamF1(sigproj, xi, beta);
2872  resid(0) = xi;
2873  resid(1) = beta;
2874  xnorm[i - 1] = Norm(diff);
2875  errnorm[i - 1] = Norm(resid - residguess);
2876  }
2877 
2878 }
2879 
2880 void TPZSandlerExtended::TaylorCheckProjectF1(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2881  TPZManVector<STATE, 3> deltasigma(3, -0.01), sigproj(3);
2882  deltasigma[1] *= 3.;
2883  deltasigma[2] *= -2.;
2884  TPZFNMatrix<3, STATE> res0(3, 1), diff(3, 1), resid(3, 1), residguess(3, 1);
2885  TPZFNMatrix<9, STATE> jac(3, 3), jacF1(2, 2), gradF1(2, 3), DF1cart(2, 3), GradSigma(3, 3);
2886  STATE kproj;
2887  ProjectF1(sigtrial, kprev, sigproj, kproj);
2888  STATE xi, beta;
2889  SurfaceParamF1(sigproj, xi, beta);
2890  res0(0) = sigproj[0];
2891  ;
2892  res0(1) = sigproj[1];
2893  res0(2) = sigproj[2];
2894  D2DistFunc1(sigtrial, xi, beta, jacF1);
2895  GradF1SigmaTrial(sigtrial, xi, beta, gradF1);
2896  //gradF1.Transpose();
2897  jacF1.Solve_LU(&gradF1);
2898  DF1Cart(xi, beta, DF1cart);
2899  DF1cart.Multiply(gradF1, GradSigma);
2900  GradSigma *= -1.;
2901 
2902  jac = GradSigma;
2903  xnorm.resize(10);
2904  errnorm.resize(10);
2905  for (int i = 1; i <= 10; i++) {
2906  TPZManVector<STATE, 3> diffsigma(3), nextsigma(3);
2907  for (int j = 0; j < 3; j++) {
2908  diffsigma[j] = deltasigma[j] * i / 10.;
2909  nextsigma[j] = sigtrial[j] + diffsigma[j];
2910  diff(j) = diffsigma[j];
2911  }
2912  jac.Multiply(diff, residguess);
2913  residguess += res0;
2914  ProjectF1(nextsigma, kprev, sigproj, kproj);
2915  resid(0) = sigproj[0];
2916  resid(1) = sigproj[1];
2917  resid(2) = sigproj[2];
2918  xnorm[i - 1] = Norm(diff);
2919  errnorm[i - 1] = Norm(resid - residguess);
2920  }
2921 
2922 }
2923 
2924 void TPZSandlerExtended::TaylorCheckProjectF2(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2925  // TPZManVector<STATE,3> deltasigma(3,-0.01), sigproj(3);
2926  // deltasigma[1] *= 3.;
2927  // deltasigma[2] *= -2.;
2928  TPZManVector<STATE, 3> deltasigma(3, -0.000012), sigproj(3);
2929  deltasigma[1] = -4.e-6;
2930  deltasigma[2] = -4.e-6;
2931  TPZFNMatrix<3, STATE> res0(3, 1), diff(3, 1), resid(3, 1), residguess(3, 1);
2932  TPZFNMatrix<9, STATE> jac(3, 3), jacF2(3, 3), gradF2(3, 3), DF2cart(3, 3), GradSigma(3, 3);
2933  STATE kproj;
2934  ProjectF2(sigtrial, kprev, sigproj, kproj);
2935  STATE theta, beta;
2936  SurfaceParamF2(sigproj, kproj, theta, beta);
2937  res0(0) = sigproj[0];
2938  res0(1) = sigproj[1];
2939  res0(2) = sigproj[2];
2940  Jacobianf2(sigtrial, theta, beta, kproj, jacF2);
2941 
2942  TFad<3, STATE> thetafad(theta, 0), betafad(beta, 1), kprojfad(kproj, 2);
2943  TPZManVector<TFad<3, STATE>, 3> sigtrialfad(3), ddistf2(3);
2944  for (int m = 0; m < 3; m++) {
2945  sigtrialfad[m] = sigtrial[m];
2946  }
2947  //DDistFunc2(sigtrialfad, thetafad, betafad, kprojfad, kprev, ddistf2);
2948 
2949  TPZFNMatrix<9, STATE> diffjac(3, 3);
2950  for (int m = 0; m < 3; m++) {
2951  for (int n = 0; n < 3; n++) {
2952  diffjac(m, n) = jacF2(m, n) - ddistf2[m].fastAccessDx(n);
2953  jacF2(m, n) = ddistf2[m].fastAccessDx(n);
2954  }
2955  }
2956 
2957 
2958  GradF2SigmaTrial(sigtrial, theta, beta, kproj, kprev, gradF2);
2959  //gradF1.Transpose();
2960  jacF2.Solve_LU(&gradF2);
2961  DF2Cart(theta, beta, kproj, DF2cart);
2962  DF2cart.Multiply(gradF2, GradSigma);
2963  GradSigma *= -1.;
2964 
2965  jac = GradSigma;
2966  xnorm.resize(10);
2967  errnorm.resize(10);
2968  for (int i = 1; i <= 10; i++) {
2969  TPZManVector<STATE, 3> diffsigma(3), nextsigma(3);
2970  for (int j = 0; j < 3; j++) {
2971  diffsigma[j] = deltasigma[j] * i / 10.;
2972  nextsigma[j] = sigtrial[j] + diffsigma[j];
2973  diff(j) = diffsigma[j];
2974  }
2975  jac.Multiply(diff, residguess);
2976  residguess += res0;
2977  ProjectF2(nextsigma, kprev, sigproj, kproj);
2978  resid(0) = sigproj[0];
2979  resid(1) = sigproj[1];
2980  resid(2) = sigproj[2];
2981  xnorm[i - 1] = Norm(diff);
2982  errnorm[i - 1] = Norm(resid - residguess);
2983  }
2984 
2985 }
2986 
2987 void TPZSandlerExtended::TaylorCheckDtbkDsigtrial(const TPZVec<STATE> &sigtrial, STATE kprev, TPZVec<STATE> &xnorm, TPZVec<STATE> &errnorm) const {
2988  // TPZManVector<STATE,3> deltasigma(3,-0.01), sigproj(3);
2989  // deltasigma[1] *= 3.;
2990  // deltasigma[2] *= -2.;
2991  TPZManVector<STATE, 3> deltasigma(3, -0.000012), sigproj(3);
2992  deltasigma[1] = -4.e-6;
2993  deltasigma[2] = -4.e-6;
2994  TPZFNMatrix<3, STATE> res0(3, 1), diff(3, 1), resid(3, 1), residguess(3, 1);
2995  TPZFNMatrix<9, STATE> jac(3, 3), jacF2(3, 3), gradF2(3, 3), GradTBK(3, 3);
2996  STATE kproj;
2997  ProjectF2(sigtrial, kprev, sigproj, kproj);
2998  STATE theta, beta;
2999  SurfaceParamF2(sigproj, kproj, theta, beta);
3000  res0(0) = theta;
3001  res0(1) = beta;
3002  res0(2) = kproj;
3003  Jacobianf2(sigtrial, theta, beta, kproj, jacF2);
3004  TFad<3, STATE> thetafad(theta, 0), betafad(beta, 1), kprojfad(kproj, 2);
3005  TPZManVector<TFad<3, STATE>, 3> sigtrialfad(3), ddistf2(3);
3006  for (int m = 0; m < 3; m++) {
3007  sigtrialfad[m] = sigtrial[m];
3008  }
3009  //DDistFunc2(sigtrialfad, thetafad, betafad, kprojfad, kprev, ddistf2);
3010 
3011  TPZFNMatrix<9, STATE> diffjac(3, 3);
3012  for (int m = 0; m < 3; m++) {
3013  for (int n = 0; n < 3; n++) {
3014  diffjac(m, n) = jacF2(m, n) - ddistf2[m].fastAccessDx(n);
3015  jacF2(m, n) = ddistf2[m].fastAccessDx(n);
3016  }
3017  }
3018  diffjac.Print("DiffMatrix");
3019 
3020  GradF2SigmaTrial(sigtrial, theta, beta, kproj, kprev, gradF2);
3021  //gradF1.Transpose();
3022  TPZFMatrix<STATE> jacF2temp(jacF2), gradF2temp(gradF2);
3023  jacF2temp.Solve_LU(&gradF2temp);
3024  GradTBK = gradF2temp;
3025  GradTBK *= -1.;
3026 
3027  jac = GradTBK;
3028 
3029  TPZFNMatrix<3, STATE> tbk(3, 1, 0.), rhs;
3030  tbk(1, 0) = 0.1;
3031  rhs = tbk;
3032  GradTBK.Solve_LU(&rhs);
3033  for (int i = 0; i < 3; i++) {
3034  deltasigma[i] = rhs(i, 0);
3035  }
3036 
3037 
3038  xnorm.resize(10);
3039  errnorm.resize(10);
3040  TPZFNMatrix<30, STATE> erros(3, 10);
3041  for (int i = 1; i <= 10; i++) {
3042  TPZManVector<STATE, 3> diffsigma(3), nextsigma(3);
3043  TPZFNMatrix<3, STATE> difftbk(3, 1);
3044  for (int j = 0; j < 3; j++) {
3045  diffsigma[j] = deltasigma[j] * i / 10.;
3046  difftbk(j) = tbk[j] * i / 10.;
3047  nextsigma[j] = sigtrial[j] + diffsigma[j];
3048  diff(j) = diffsigma[j];
3049  }
3050  jac.Multiply(diff, residguess);
3051  residguess += res0;
3052  // TPZFNMatrix<3,STATE> multsig(3,0),multbk(3,0);
3053  // jacF2.Multiply(difftbk, multbk);
3054  // gradF2.Multiply(diff, multsig);
3055  // residguess = multbk;
3056  // residguess = multsig;
3057  // TPZFNMatrix<3,STATE> fxn(3,1);
3058  // DDistFunc2(nextsigma, theta+difftbk[0],beta+difftbk[1],kproj+difftbk[2],kprev,fxn);
3059  // TPZManVector<STATE> fxnvec(3);
3060  // DDistFunc2<STATE>(sigtrial,theta+difftbk[0],beta+difftbk[1],kproj+difftbk[2],kprev,fxnvec);
3061  // for(int k=0; k<3; k++) fxn(k,0) = fxnvec[k];
3062  // DDistFunc2(sigtrial, theta+difftbk[0],beta+difftbk[1],kproj+difftbk[2],kprev,fxn);
3063  // DDistFunc2(nextsigma, theta,beta,kproj,kprev,fxn);
3064  ProjectF2(nextsigma, kprev, sigproj, kproj);
3065  STATE theta, beta;
3066  SurfaceParamF2(sigproj, kproj, theta, beta);
3067  resid(0) = theta;
3068  resid(1) = beta;
3069  resid(2) = kproj;
3070  xnorm[i - 1] = Norm(diff);
3071  errnorm[i - 1] = Norm(resid - residguess);
3072  for (int a = 0; a < 3; a++) erros(a, i - 1) = fabs(resid[a] - residguess[a]);
3073  }
3074  erros.Print(cout);
3075 
3076 }
3077 
3078 void TPZSandlerExtended::MCormicRanchSand(TPZSandlerExtended &mat)//em ksi
3079 {
3080  STATE E = 100, nu = 0.25, A = 0.25, B = 0.67, C = 0.18, D = 0.67, R = 2.5, W = 0.066, N = 0., phi = 0, psi = 1.0;
3081  STATE G = E / (2. * (1. + nu));
3082  STATE K = E / (3. * (1. - 2 * nu));
3083  mat.fA = A;
3084  mat.fB = B;
3085  mat.fC = C;
3086  mat.fD = D;
3087  mat.fK = K;
3088  mat.fG = G;
3089  mat.fW = W;
3090  mat.fR = R;
3091  mat.fPhi = phi;
3092  mat.fN = N;
3093  mat.fPsi = psi;
3094  mat.fE = E;
3095  mat.fnu = nu;
3096  TPZElasticResponse ER;
3097  ER.SetEngineeringData(E, nu);
3098  mat.fElasticResponse = ER;
3099 
3100 }
3101 
3102 void TPZSandlerExtended::ReservoirSandstone(TPZSandlerExtended &mat)//em ksi
3103 {
3104  STATE E = 1305, nu = 0.25, A = 2.61, B = 0.169, C = 2.57, D = 0.05069, R = 1.5, W = 0.0908, N = 0., phi = 0, psi = 1.0;
3105  STATE G = E / (2. * (1. + nu));
3106  STATE K = E / (3. * (1. - 2 * nu));
3107  mat.fA = A;
3108  mat.fB = B;
3109  mat.fC = C;
3110  mat.fD = D;
3111  mat.fK = K;
3112  mat.fG = G;
3113  mat.fW = W;
3114  mat.fR = R;
3115  mat.fPhi = phi;
3116  mat.fN = N;
3117  mat.fPsi = psi;
3118  mat.fE = E;
3119  mat.fnu = nu;
3120  TPZElasticResponse ER;
3121  ER.SetEngineeringData(E, nu);
3122  mat.fElasticResponse = ER;
3123 
3124 
3125 }
3126 
3127 void TPZSandlerExtended::SalemLimestone(TPZSandlerExtended &mat)// em MPa
3128 {
3129  STATE E = 22547., nu = 0.2524, A = 689.2,
3130  B = 3.94e-4, C = 675.2, D = 1.47e-3, R = 28, W = 0.08, N = 6., phi = 0, psi = 1.0;
3131  STATE G = E / (2. * (1. + nu));
3132  STATE K = E / (3. * (1. - 2 * nu));
3133  mat.fA = A;
3134  mat.fB = B;
3135  mat.fC = C;
3136  mat.fD = D;
3137  mat.fK = K;
3138  mat.fG = G;
3139  mat.fW = W;
3140  mat.fR = R;
3141  mat.fPhi = phi;
3142  mat.fN = N;
3143  mat.fPsi = psi;
3144  mat.fE = E;
3145  mat.fnu = nu;
3146  TPZElasticResponse ER;
3147  ER.SetEngineeringData(E, nu);
3148  mat.fElasticResponse = ER;
3149 }
3150 
3151 void TPZSandlerExtended::PreSMat(TPZSandlerExtended &mat)// em MPa
3152 {
3153 
3154  STATE E = 29269, nu = 0.203, A = 116.67,
3155  B = 0.0036895, C = 111.48, D = 0.018768, R = 0.91969, W = 0.006605, N = 0., phi = 0, psi = 1.0;
3156  STATE G = E / (2. * (1. + nu));
3157  STATE K = E / (3. * (1. - 2 * nu));
3158  mat.fA = A;
3159  mat.fB = B;
3160  mat.fC = C;
3161  mat.fD = D;
3162  mat.fK = K;
3163  mat.fG = G;
3164  mat.fW = W;
3165  mat.fR = R;
3166  mat.fPhi = phi;
3167  mat.fN = N;
3168  mat.fPsi = psi;
3169  mat.fE = E;
3170  mat.fnu = nu;
3171  TPZElasticResponse ER;
3172  ER.SetEngineeringData(E, nu);
3173  mat.fElasticResponse = ER;
3174 }
3175 
3176 int TPZSandlerExtended::ClassId() const {
3177  return Hash("TPZSandlerExtended");
3178 }
TPZSandlerExtended::DF2Cart
void DF2Cart(STATE theta, STATE beta, STATE k, TPZFMatrix< STATE > &DF1) const
Compute the derivative of the stress (principal s;tresses) as a function of xi and beta...
Definition: TPZSandlerExtended.cpp:2713
TPZVec::Print
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
atan2
TPZFlopCounter atan2(const TPZFlopCounter &val1, const TPZFlopCounter &val2)
Returns the arc tangent in radians and increments the counter of the Arc Tangent. ATAN2 returns the ...
Definition: pzreal.h:544
TPZSandlerExtended::TaylorCheckDistF1
void TaylorCheckDistF1(const TPZVec< STATE > &sigmatrial, STATE xi, STATE beta, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2493
TPZSandlerExtended::YieldFunction
virtual void YieldFunction(const TPZVec< STATE > &sigma, STATE kprev, TPZVec< STATE > &yield) const override
Definition: TPZSandlerExtended.cpp:1052
TPZSandlerExtended::GetElasticResponse
TPZElasticResponse GetElasticResponse()
Definition: TPZSandlerExtended.cpp:157
TPZSandlerExtended::ProjectVertex
void ProjectVertex(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1380
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZSandlerExtended::Res2CoVertex
void Res2CoVertex(const TPZVec< T > &trial_stress, T beta, T k, T kprev, TPZVec< T > &residue_covertex) const
Compute the derivative of the distance function to the covertex cap function and the result of covert...
Definition: TPZSandlerExtended.cpp:629
TPZSandlerExtended::F2Cyl
void F2Cyl(const STATE theta, const STATE beta, const STATE k, TPZVec< STATE > &f2cyl) const
Definition: TPZSandlerExtended.cpp:374
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
TPZSandlerExtended::~TPZSandlerExtended
virtual ~TPZSandlerExtended()
Desctructor.
Definition: TPZSandlerExtended.cpp:56
TPZMatrix::Solve_LU
int Solve_LU(TPZFMatrix< TVar > *B, std::list< int64_t > &singular)
Solves the linear system using LU method .
Definition: pzmatrix.h:900
pzlog.h
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
TPZHWTools::FromPrincipalToHWCyl
static void FromPrincipalToHWCyl(const TPZVec< REAL > &PrincipalCoords, TPZVec< REAL > &HWCyl)
Definition: TPZHWTools.h:39
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
TPZSandlerExtended::ProjectF2
void ProjectF2(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1246
TPZInconsistentStateException.h
TPZSandlerExtended::ResLF2IJ
T ResLF2IJ(const TPZVec< T > &sigtrIJ, T theta, T k, STATE kprev) const
Compute the residual of the equation which defines the update of the damage variable.
Definition: TPZSandlerExtended.cpp:298
TPZSandlerExtended::fElasticResponse
TPZElasticResponse fElasticResponse
Definition: TPZSandlerExtended.h:404
TPZSandlerExtended::ApplyStressComputeElasticStrain
void ApplyStressComputeElasticStrain(TPZVec< STATE > &stress, TPZVec< STATE > &strain) const
Definition: TPZSandlerExtended.cpp:1714
TPZSandlerExtended::TaylorCheckParamF1Sigtrial
void TaylorCheckParamF1Sigtrial(const TPZVec< STATE > &sigtrial, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
verifies the validity of dxi/dsigtrial and dbeta/dsigtrial
Definition: TPZSandlerExtended.cpp:2842
Sin
#define Sin
Definition: TPZSandlerExtended.cpp:2442
TPZSandlerExtended::ProjectCapCoVertex
void ProjectCapCoVertex(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1468
TPZVec::resize
virtual void resize(const int64_t newsize)
Definition: pzvec.h:213
TPZSandlerExtended::Print
virtual void Print(std::ostream &out) const override
Definition: TPZSandlerExtended.cpp:1097
TPZSandlerExtended.h
TPZSandlerExtended::TPZSandlerExtended
TPZSandlerExtended()
Empty constructor.
Definition: TPZSandlerExtended.cpp:24
TPZSandlerExtended::TaylorCheckProjectSigma
void TaylorCheckProjectSigma(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Compute the approximation rate for the derivative of the projected stresses respect to trial stresses...
Definition: TPZSandlerExtended.cpp:2808
TPZSandlerExtended::TaylorCheckDF2Cart
void TaylorCheckDF2Cart(STATE theta, STATE beta, STATE k, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2770
TPZSandlerExtended::Jacobianf1
void Jacobianf1(const TPZVec< STATE > &trial_stress, STATE i1, STATE beta, STATE k, TPZFMatrix< STATE > &jacobianf1) const
Compute the jacobian function of the f1 (failure) distance as a function of i1, beta and k...
Definition: TPZSandlerExtended.cpp:674
TPZElasticResponse::Read
void Read(TPZStream &buf, void *context) override
Definition: TPZElasticResponse.cpp:38
TPZSandlerExtended::NormalFunctionToF1
STATE NormalFunctionToF1(STATE &I1, STATE &k) const
Definition: TPZSandlerExtended.cpp:1123
TPZHWTools.h
TPZSandlerExtended::TaylorCheckDDistF1
void TaylorCheckDDistF1(const TPZVec< STATE > &sigmatrial, STATE xi, STATE beta, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2515
TPZSandlerExtended::DistF2
STATE DistF2(const TPZVec< STATE > &pt, const STATE theta, const STATE beta, const STATE k) const
Compute the distance of sigtrial to the point on the cap.
Definition: TPZSandlerExtended.cpp:426
TPZMatrix::Identity
virtual void Identity()
Converts the matrix in an identity matrix.
Definition: pzmatrix.cpp:199
TPZSandlerExtended::X
T X(const T k) const
The function which defines the plastic surface.
Definition: TPZSandlerExtended.cpp:75
TPZSandlerExtended::ComputeCapTangent
void ComputeCapTangent(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj, TPZFMatrix< REAL > *gradient) const
Compute the derivative of the projected stresses respect to trial stresses (tangent) over the cap...
Definition: TPZSandlerExtended.cpp:1913
TPZSandlerExtended::CheckCoordinateTransformation
static void CheckCoordinateTransformation(TPZVec< STATE > &cart)
Definition: TPZSandlerExtended.cpp:2672
TPZSandlerExtended::FromThetaKToSigIJ
void FromThetaKToSigIJ(const T &theta, const T &K, TPZVec< T > &sigIJ) const
Definition: TPZSandlerExtended.cpp:450
TPZSandlerExtended::Write
void Write(TPZStream &buf, int withclassid) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: TPZSandlerExtended.cpp:138
TPZSandlerExtended::GradF1SigmaTrial
void GradF1SigmaTrial(const TPZVec< STATE > &sigtrial, STATE xi, STATE beta, TPZFMatrix< STATE > &deriv) const
Compute the derivative of the residual with respect to sigtrial.
Definition: TPZSandlerExtended.cpp:2448
TPZSandlerExtended::SetElasticResponse
void SetElasticResponse(const TPZElasticResponse &ER)
Definition: TPZSandlerExtended.cpp:107
TPZSandlerExtended::SetInitialDamage
void SetInitialDamage(STATE kappa_0)
Definition: TPZSandlerExtended.cpp:103
TPZSandlerExtended::TaylorCheckProjectF2
void TaylorCheckProjectF2(const TPZVec< STATE > &sigtrial, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2924
TPZSandlerExtended::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: TPZSandlerExtended.cpp:3176
TPZSandlerExtended::NormalToF1
STATE NormalToF1(STATE I1, STATE I1_ref) const
Compute the normal function to the failure surface based on a reference point (I1_ref,f1(I1_ref))
Definition: TPZSandlerExtended.cpp:198
TPZSandlerExtended::ProjectRing
void ProjectRing(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1523
TPZSandlerExtended::ApplyStrainComputeSigma
void ApplyStrainComputeSigma(TPZVec< STATE > &epst, TPZVec< STATE > &epsp, STATE &kprev, TPZVec< STATE > &epspnext, TPZVec< STATE > &stressnext, STATE &knext, TPZFMatrix< REAL > *tangent=NULL) const
Definition: TPZSandlerExtended.cpp:1735
TPZSandlerExtended::ConvergenceRate
static void ConvergenceRate(TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm, TPZVec< STATE > &convergence)
Definition: TPZSandlerExtended.cpp:2566
TPZSandlerExtended::TaylorCheckDDistF2DSigtrial
void TaylorCheckDDistF2DSigtrial(const TPZVec< STATE > &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
teste da derivada D(ResF2)/D(sigtrial)
Definition: TPZSandlerExtended.cpp:2639
TPZSandlerExtended::Firstk
void Firstk(STATE &epsp, STATE &k) const
Compute k as a function of epsp using Newton iterations.
Definition: TPZSandlerExtended.cpp:176
TPZSandlerExtended::ProjectSigma
void ProjectSigma(const TPZVec< STATE > &sigmatrial, STATE kprev, TPZVec< STATE > &sigmaproj, STATE &kproj, int &m_type, TPZFMatrix< REAL > *gradient=NULL) const
Definition: TPZSandlerExtended.cpp:1796
TPZSandlerExtended::ResLF2
T ResLF2(const TPZVec< T > &pt, T theta, T beta, T k, STATE kprev) const
Compute the residual of the equation which defines the update of the damage variable.
Definition: TPZSandlerExtended.cpp:290
TPZSandlerExtended::ApplyStrainComputeElasticStress
void ApplyStrainComputeElasticStress(TPZVec< STATE > &strain, TPZVec< STATE > &stress) const
Definition: TPZSandlerExtended.cpp:1699
TPZSandlerExtended::ProjectCapVertex
void ProjectCapVertex(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1422
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZSandlerExtended::ProjectCoVertex
void ProjectCoVertex(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1322
acos
expr_ expr_ expr_ expr_ acos
Definition: tfadfunc.h:75
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZSandlerExtended::SalemLimestone
static void SalemLimestone(TPZSandlerExtended &mat)
Definition: TPZSandlerExtended.cpp:3127
TPZHWTools::FromHWCylToHWCart
static void FromHWCylToHWCart(const TPZVec< REAL > &HWCylCoords, TPZVec< REAL > &HWCart)
Definition: TPZHWTools.h:28
sin
sin
Definition: fadfunc.h:63
TPZSandlerExtended::GradF2SigmaTrial
void GradF2SigmaTrial(const TPZVec< STATE > &sigtrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZFMatrix< STATE > &deriv) const
Compute the derivative of the F2 residual with respecto do sigtrial.
Definition: TPZSandlerExtended.cpp:2471
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
test.f
f
Definition: test.py:287
TPZSandlerExtended::TaylorCheckDDistF1DSigtrial
void TaylorCheckDDistF1DSigtrial(const TPZVec< STATE > &sigmatrial, STATE xi, STATE beta, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2541
TPZSandlerExtended::TaylorCheckDistF2
void TaylorCheckDistF2(const TPZVec< STATE > &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2573
TPZSandlerExtended::ProjectSigmaDep
void ProjectSigmaDep(const TPZVec< STATE > &sigmatrial, STATE kprev, TPZVec< STATE > &sigmaproj, STATE &kproj, TPZFMatrix< STATE > &GradSigma) const
Definition: TPZSandlerExtended.cpp:2211
TPZSandlerExtended::SurfaceParamF1
void SurfaceParamF1(TPZVec< STATE > &sigproj, STATE &xi, STATE &beta) const
Definition: TPZSandlerExtended.cpp:361
TPZSandlerExtended::DResLF2
void DResLF2(const TPZVec< STATE > &pt, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec< STATE > &dresl) const
Compute the derivative of the equation which determines the evolution of k.
Definition: TPZSandlerExtended.cpp:329
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZHWTools::A2x2Inverse
static void A2x2Inverse(TPZFMatrix< REAL > &A, TPZFMatrix< REAL > &Ainv)
Definition: TPZHWTools.h:78
TPZSandlerExtended::Jacobianf2CoVertex
void Jacobianf2CoVertex(const TPZVec< STATE > &trial_stress, STATE beta, STATE k, TPZFMatrix< STATE > &jacobianf2_covertex) const
Compute the jacobian function of the f2 (cap) distance as a function of beta and k.
Definition: TPZSandlerExtended.cpp:886
TPZHWTools::GetRotMatrix
static void GetRotMatrix(TPZFMatrix< REAL > &Rot)
Computes the rotation matrix.
Definition: TPZHWTools.h:113
Norm
TVar Norm(const TPZFMatrix< TVar > &A)
Returns the norm of the matrix A.
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZSandlerExtended::Phi
void Phi(TPZVec< REAL > sigma, STATE alpha, TPZVec< STATE > &phi) const
Definition: TPZSandlerExtended.cpp:1109
Cos
#define Cos
Definition: TPZSandlerExtended.cpp:2441
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZSandlerExtended::ReservoirSandstone
static void ReservoirSandstone(TPZSandlerExtended &mat)
Definition: TPZSandlerExtended.cpp:3102
TFad
Definition: tfad.h:64
TPZSandlerExtended::TaylorCheckDDistF2
void TaylorCheckDDistF2(const TPZVec< STATE > &sigmatrial, STATE theta, STATE beta, STATE k, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2601
TPZHWTools::A3x3Inverse
static void A3x3Inverse(TPZFMatrix< REAL > &A, TPZFMatrix< REAL > &Ainv)
Definition: TPZHWTools.h:93
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
TPZSandlerExtended::JacobianVertex
void JacobianVertex(const TPZVec< STATE > &trial_stress, STATE k, STATE &jacobian_vertex) const
Compute the jacobian of the distance function to the cap vertex function and the result of Vertex res...
Definition: TPZSandlerExtended.cpp:964
TPZSandlerExtended::DistF1
STATE DistF1(const TPZVec< STATE > &pt, const STATE xi, const STATE beta) const
Compute the distance of sigtrial to the point on the yield surface.
Definition: TPZSandlerExtended.cpp:414
TPZSandlerExtended::DResLF1
STATE DResLF1(const TPZVec< STATE > &sigtrial, const TPZVec< STATE > &sigproj, const STATE k, const STATE kprev) const
Compute the derivative of the equation which determines the evolution of k.
Definition: TPZSandlerExtended.cpp:319
test.res
string res
Definition: test.py:151
TPZSandlerExtended::ProjectApex
void ProjectApex(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1133
TPZSandlerExtended::DDistFunc1
void DDistFunc1(const TPZVec< STATE > &pt, STATE xi, STATE beta, TPZFMatrix< STATE > &ddistf1) const
Compute the derivative of the distance function to the yield surface as a function of xi and beta...
Definition: TPZSandlerExtended.cpp:478
TPZSandlerExtended::ComputeJ2
void ComputeJ2(TPZVec< STATE > stress, STATE &J2) const
Definition: TPZSandlerExtended.cpp:1689
TPZSandlerExtended::JacobianCoVertex
void JacobianCoVertex(const TPZVec< STATE > &trial_stress, STATE beta, STATE k, TPZFMatrix< STATE > &jacobian_covertex) const
Compute the jacobian of the distance function to the cap covertex function and the result of Covertex...
Definition: TPZSandlerExtended.cpp:977
TPZSandlerExtended::ComputeFailureTangent
void ComputeFailureTangent(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj, TPZFMatrix< REAL > *gradient) const
Compute the derivative of the projected stresses respect to trial stresses (tangent) over the failure...
Definition: TPZSandlerExtended.cpp:2127
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZHWTools::FromHWCylToPrincipal
static void FromHWCylToPrincipal(const TPZVec< REAL > &HWCylCoords, TPZVec< REAL > &PrincipalCoords)
Definition: TPZHWTools.h:67
TPZConvergenceException.h
TPZSandlerExtended::DF1Cart
void DF1Cart(STATE xi, STATE beta, TPZFMatrix< STATE > &DF1) const
Compute the derivative of the stress (principal s;tresses) as a function of xi and beta...
Definition: TPZSandlerExtended.cpp:2686
TPZSandlerExtended::EpsEqX
T EpsEqX(T X) const
compute the damage variable as a function of the X function
Definition: TPZSandlerExtended.cpp:166
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZSandlerExtended::PreSMat
static void PreSMat(TPZSandlerExtended &mat)
Definition: TPZSandlerExtended.cpp:3151
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
TPZSandlerExtended::Jacobianf2Vertex
void Jacobianf2Vertex(const TPZVec< STATE > &trial_stress, STATE k, STATE &jacobianf2_vertex) const
Compute the jacobian function of the vertex on f2 (cap) distance as a function of k...
Definition: TPZSandlerExtended.cpp:950
gF2Stat
std::map< int, int64_t > gF2Stat
Definition: TPZSandlerExtended.cpp:1120
TPZSandlerExtended::ProjectF1
void ProjectF1(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1181
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZHWTools::FromPrincipalToHWCart
static void FromPrincipalToHWCart(const TPZVec< REAL > &PrincipalCoords, TPZVec< REAL > &HWCart)
Definition: TPZHWTools.h:52
gF1Stat
std::map< int, int64_t > gF1Stat
Definition: TPZSandlerExtended.cpp:1119
CG
int CG(Matrix &A, Vector &x, const Vector &b, Preconditioner &M, Vector *residual, int64_t &max_iter, Real &tol, const int FromCurrent)
CG solves the symmetric positive definite linear system using the Conjugate Gradient method...
Definition: cg.h:31
TPZSandlerExtended::Res1
void Res1(const TPZVec< T > &trial_stress, T i1, T beta, T k, T kprev, TPZVec< T > &residue_1) const
Compute the derivative of the distance function to the failure function and the result of Residue 1 (...
Definition: TPZSandlerExtended.cpp:550
TPZSandlerExtended::TaylorCheckProjectF1
void TaylorCheckProjectF1(const TPZVec< STATE > &sigtrial, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2880
dist
REAL dist(TPZVec< T1 > &vec1, TPZVec< T1 > &vec2)
Definition: pzvec_extras.h:124
log
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ log
Definition: tfadfunc.h:130
TPZSandlerExtended::ResLF1
STATE ResLF1(const TPZVec< STATE > &sigtrial, const TPZVec< STATE > &sigproj, const STATE k, const STATE kprev) const
Compute the residual of the equation which defines the update of the damage variable.
Definition: TPZSandlerExtended.cpp:308
TPZSandlerExtended::DistF2IJ
STATE DistF2IJ(const TPZVec< STATE > &sigtrialIJ, STATE theta, STATE k) const
Compute the distance considering the sigtrial is given as a funcion of I1, sqJ2.
Definition: TPZSandlerExtended.cpp:437
TPZElasticResponse::SetEngineeringData
void SetEngineeringData(REAL Eyoung, REAL Poisson)
Definition: TPZElasticResponse.cpp:76
TPZSandlerExtended::D2DistFunc1
void D2DistFunc1(const TPZVec< STATE > &pt, STATE xi, STATE beta, TPZFMatrix< STATE > &d2distf1) const
Compute the second derivative of the distance as a function of xi and beta.
Definition: TPZSandlerExtended.cpp:505
TPZSandlerExtended::CPerturbation
STATE CPerturbation() const
Definition: TPZSandlerExtended.h:391
gYield
std::vector< int64_t > gYield
Definition: TPZSandlerExtended.cpp:1121
TPZSandlerExtended::DDistF2IJ
void DDistF2IJ(TPZVec< T > &sigtrialIJ, T theta, T L, STATE Lprev, TPZVec< T > &ddistf2) const
Compute the value of the equation which determines the orthogonality of the projection.
Definition: TPZSandlerExtended.cpp:460
TPZSandlerExtended::MCormicRanchSand
static void MCormicRanchSand(TPZSandlerExtended &mat)
Definition: TPZSandlerExtended.cpp:3078
TPZSandlerExtended::Res2Vertex
void Res2Vertex(const TPZVec< T > &trial_stress, T k, T kprev, T &residue_vertex) const
Compute the derivative of the distance function to the vertex cap function and the result of vertex R...
Definition: TPZSandlerExtended.cpp:658
TPZMatrix::Multiply
virtual void Multiply(const TPZFMatrix< TVar > &A, TPZFMatrix< TVar > &res, int opt=0) const
It mutiplies itself by TPZMatrix<TVar>A putting the result in res.
Definition: pzmatrix.cpp:1916
TPZSandlerExtended::Jacobianf2
void Jacobianf2(const TPZVec< STATE > &trial_stress, STATE theta, STATE beta, STATE k, TPZFMatrix< STATE > &jacobianf2) const
Compute the jacobian function of the f2 (cap) distance as a function of theta, beta and k...
Definition: TPZSandlerExtended.cpp:758
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
exp
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ exp
Definition: tfadfunc.h:125
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZSandlerExtended::GetF
STATE GetF(STATE x) const
Definition: TPZSandlerExtended.cpp:70
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZSandlerExtended::SetUp
void SetUp(STATE A, STATE B, STATE C, STATE D, STATE K, STATE G, STATE W, STATE R, STATE Phi, STATE N, STATE Psi)
Definition: TPZSandlerExtended.cpp:83
TPZSandlerExtended::TaylorCheckDtbkDsigtrial
void TaylorCheckDtbkDsigtrial(const TPZVec< STATE > &sigtrial, STATE kprev, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
verify D(theta,beta,k)/D(sigtrial)
Definition: TPZSandlerExtended.cpp:2987
pzreferredcompel.h
Contains declaration of TPZReferredCompEl class which generates computational elements.
TPZSandlerExtended::Read
void Read(TPZStream &buf, void *context) override
read objects from the stream
Definition: TPZSandlerExtended.cpp:119
TPZSandlerExtended::TaylorCheckDF1Cart
void TaylorCheckDF1Cart(STATE xi, STATE beta, TPZVec< STATE > &xnorm, TPZVec< STATE > &errnorm) const
Definition: TPZSandlerExtended.cpp:2737
pzvec_extras.h
Extra utilities for TPZVec. Implementations of the saxpy, sscal, sdot, intercept, max and min functio...
TPZSandlerExtended::ComputeI1
void ComputeI1(TPZVec< STATE > stress, STATE &I1) const
Definition: TPZSandlerExtended.cpp:1680
m
clarg::argString m("-m", "input matrix file name (text format)", "matrix.txt")
TPZSandlerExtended::ComputeCapVertexTangent
void ComputeCapVertexTangent(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj, TPZFMatrix< REAL > *gradient) const
Compute the derivative of the projected stresses respect to trial stresses (tangent) over the cap...
Definition: TPZSandlerExtended.cpp:2008
cos
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TPZSandlerExtended::SurfaceParamF2
void SurfaceParamF2(const TPZVec< STATE > &sigproj, const STATE k, STATE &theta, STATE &beta) const
Definition: TPZSandlerExtended.cpp:389
TPZSandlerExtended::Res2
void Res2(const TPZVec< T > &trial_stress, T theta, T beta, T k, T kprev, TPZVec< T > &residue_2) const
Compute the derivative of the distance function to the cap function and the result of Residue 2 (Cap)...
Definition: TPZSandlerExtended.cpp:586
TPZSandlerExtended::ComputeCapCoVertexTangent
void ComputeCapCoVertexTangent(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj, TPZFMatrix< REAL > *gradient) const
Compute the derivative of the projected stresses respect to trial stresses (tangent) over the cap...
Definition: TPZSandlerExtended.cpp:2054
TPZSandlerExtended::ProjectBetaConstF2
void ProjectBetaConstF2(const TPZVec< STATE > &trial_stress, STATE kprev, TPZVec< STATE > &projected_stress, STATE &kproj) const
Definition: TPZSandlerExtended.cpp:1609
TPZSandlerExtended::DF
T DF(const T x) const
Definition: TPZSandlerExtended.cpp:66
TPZElasticResponse::Write
void Write(TPZStream &buf, int withclassid) const override
Definition: TPZElasticResponse.cpp:32
TPZSandlerExtended::F1Cyl
void F1Cyl(STATE xi, STATE beta, TPZVec< STATE > &f1cyl) const
Compute the point on F1 in HW Cylindrical coordinates.
Definition: TPZSandlerExtended.cpp:347
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZSandlerExtended::F
T F(const T x) const
Definition: TPZSandlerExtended.cpp:61
TPZSandlerExtended::EpsEqk
T EpsEqk(const T k) const
Compute the damage variable as a function of the position of the cap k.
Definition: TPZSandlerExtended.cpp:172
Generated on Sun Aug 16 2020 11:17:02 for NeoPZ by   doxygen 1.8.13