NeoPZ
TPZMatElastoPlastic2D_impl.h
Go to the documentation of this file.
1 //
2 // TPZMatElastoPlastic2D_impl.h
3 // pz
4 //
5 // Created by Omar Durán on 9/8/18.
6 //
7 
8 #include "TPZMatElastoPlastic2D.h"
9 #include "TPZMatElastoPlastic_impl.h"
10 
11 #include "pzbndcond.h"
12 #include "TPZBndCondWithMem.h"
13 
14 #ifndef WIN32
15 #include <fenv.h>//NAN DETECTOR
16 #endif
17 
18 #ifdef LOG4CXX
19 #include "pzlog.h"
20 static LoggerPtr elastoplastic2dLogger(Logger::getLogger("material.pzElastoPlastic2D"));
21 #endif
22 
23 
24 template <class T, class TMEM>
25 TPZMatElastoPlastic2D<T,TMEM>::TPZMatElastoPlastic2D() : TPZMatElastoPlastic<T,TMEM>()
26 {
27  fPlaneStrain = true;
28 }
29 
30 template <class T, class TMEM>
31 TPZMatElastoPlastic2D<T,TMEM>::TPZMatElastoPlastic2D(int id , int PlaneStrainOrPlaneStress) : TPZMatElastoPlastic<T,TMEM>(id)
32 {
33  fPlaneStrain = PlaneStrainOrPlaneStress;
34 }
35 
36 template <class T, class TMEM>
37 TPZMatElastoPlastic2D<T,TMEM>::TPZMatElastoPlastic2D(const TPZMatElastoPlastic2D<T,TMEM> &mat) : TPZMatElastoPlastic<T,TMEM>(mat), fPlaneStrain(mat.fPlaneStrain)
38 {
39 #ifdef LOG4CXX
40  if(elastoplastic2dLogger->isDebugEnabled())
41  {
42  std::stringstream sout;
43  sout << ">>> TPZMatElastoPlastic2D<T,TMEM>() copy constructor called ***";
44  LOGPZ_INFO(elastoplastic2dLogger,sout.str().c_str());
45  }
46 #endif
47 }
48 
49 template <class T, class TMEM>
50 TPZMatElastoPlastic2D<T,TMEM>::~TPZMatElastoPlastic2D()
51 {
52 
53 }
54 
55 template <class T, class TMEM>
56 void TPZMatElastoPlastic2D<T,TMEM>::ApplyDeltaStrain(TPZMaterialData & data, TPZFMatrix<REAL> & DeltaStrain,TPZFMatrix<REAL> & Stress)
57 {
58 
59 
60 #ifdef PZDEBUG
61  if (DeltaStrain.Rows() != 6) {
62  DebugStop();
63  }
64 #endif
65 
66  if (!fPlaneStrain) //
67  {//
68  DebugStop();//PlaneStress
69  }
70 
71  TPZMatElastoPlastic<T,TMEM>::ApplyDeltaStrain(data,DeltaStrain,Stress);//
72 
73 }
74 
75 
76 template <class T, class TMEM>
77 void TPZMatElastoPlastic2D<T,TMEM>::ApplyDeltaStrainComputeDep(TPZMaterialData & data, TPZFMatrix<REAL> & DeltaStrain,TPZFMatrix<REAL> & Stress, TPZFMatrix<REAL> & Dep)
78 {
79 #ifdef PZDEBUG
80  if (DeltaStrain.Rows() != 6) {
81  DebugStop();
82  }
83 #endif
84  if (!fPlaneStrain) //
85  {//
86  DebugStop();//PlaneStress
87  }
88  TPZMatElastoPlastic<T,TMEM>::ApplyDeltaStrainComputeDep(data,DeltaStrain,Stress,Dep);//
89 }
90 
91 template <class T, class TMEM>
92 void TPZMatElastoPlastic2D<T, TMEM>::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<REAL> &ek, TPZFMatrix<REAL> &ef) {
93 
94  TPZFMatrix<REAL> &dphi = data.dphix, dphiXY;
95  TPZFMatrix<REAL> &phi = data.phi;
96  TPZFMatrix<REAL> &axes = data.axes, axesT;
97 
98  axes.Transpose(&axesT);
99  axesT.Multiply(dphi, dphiXY);
100 
101  const int phr = phi.Rows();
102 
103  TPZFNMatrix<4> Deriv(2, 2);
104  TPZFNMatrix<36> Dep(6, 6,0.0);
105  TPZFNMatrix<6> DeltaStrain(6, 1);
106  TPZFNMatrix<6> Stress(6, 1);
107  int ptindex = data.intGlobPtIndex;
108 
109  if (TPZMatWithMem<TMEM>::fUpdateMem && data.sol.size() > 0) {
110  // Loop over the solutions if update memory is true
111  TPZSolVec locsol(data.sol);
112  TPZGradSolVec locdsol(data.dsol);
113  int numsol = locsol.size();
114 
115  for (int is = 0; is < numsol; is++) {
116  data.sol[0] = locsol[is];
117  data.dsol[0] = locdsol[is];
118 
119  this->ComputeDeltaStrainVector(data, DeltaStrain);
120  this->ApplyDeltaStrainComputeDep(data, DeltaStrain, Stress, Dep);
121  }
122  } else {
123  this->ComputeDeltaStrainVector(data, DeltaStrain);
124  this->ApplyDeltaStrainComputeDep(data, DeltaStrain, Stress, Dep);
125  }
126 
127 #ifdef MACOS
128  feclearexcept(FE_ALL_EXCEPT);
129  if (fetestexcept(/*FE_DIVBYZERO*/ FE_ALL_EXCEPT)) {
130  std::cout << "division by zero reported\n";
131  DebugStop();
132  }
133 #endif
134 
135 #ifdef LOG4CXX
136  if (elastoplastic2dLogger->isDebugEnabled()) {
137  std::stringstream sout;
138  sout << ">>> TPZMatElastoPlastic<T,TMEM>::Contribute ***";
139  sout << "\nIntegration Local Point index = " << data.intGlobPtIndex;
140  sout << "\nIntegration Global Point index = " << data.intGlobPtIndex;
141  sout << "\ndata.axes = " << data.axes;
142  sout << "\nDep " << endl;
143  sout << Dep(_XX_, _XX_) << "\t" << Dep(_XX_, _YY_) << "\t" << Dep(_XX_, _XY_) << "\n";
144  sout << Dep(_YY_, _XX_) << "\t" << Dep(_YY_, _YY_) << "\t" << Dep(_YY_, _XY_) << "\n";
145  sout << Dep(_XY_, _XX_) << "\t" << Dep(_XY_, _YY_) << "\t" << Dep(_XY_, _XY_) << "\n";
146 
147  sout << "\nStress " << endl;
148  sout << Stress(_XX_, 0) << "\t" << Stress(_YY_, 0) << "\t" << Stress(_XY_, 0) << "\n";
149 
150  sout << "\nDELTA STRAIN " << endl;
151  sout << DeltaStrain(0, 0) << "\t" << DeltaStrain(1, 0) << "\t" << DeltaStrain(2, 0) << "\n";
152  sout << "data.phi" << data.phi;
153 
154  LOGPZ_DEBUG(elastoplastic2dLogger, sout.str().c_str());
155  }
156 #endif
157  ptindex = 0;
158  int nstate = NStateVariables();
159  REAL val;
160 
161  TPZManVector<STATE, 2> ForceLoc(nstate,0.0);
162  if (this->fForcingFunction) {
163  this->fForcingFunction->Execute(data.x, ForceLoc);
164  }
165 
166  int in;
167  for (in = 0; in < phr; in++) {
168 
169  val = ForceLoc[0] * phi(in, 0);
170  val -= Stress(_XX_, 0) * dphiXY(0, in);
171  val -= Stress(_XY_, 0) * dphiXY(1, in);
172  ef(in * nstate + 0, 0) += weight * val;
173 
174  val = ForceLoc[1] * phi(in, 0);
175  val -= Stress(_XY_, 0) * dphiXY(0, in);
176  val -= Stress(_YY_, 0) * dphiXY(1, in);
177  ef(in * nstate + 1, 0) += weight * val;
178 
179  for (int jn = 0; jn < phr; jn++) {
180  for (int ud = 0; ud < 2; ud++) {
181  for (int vd = 0; vd < 2; vd++) {
182  Deriv(vd, ud) = dphiXY(vd, in) * dphiXY(ud, jn);
183  }
184  }
185 
186  val = 2. * Dep(_XX_, _XX_) * Deriv(0, 0); //dvdx*dudx
187  val += Dep(_XX_, _XY_) * Deriv(0, 1); //dvdx*dudy
188  val += 2. * Dep(_XY_, _XX_) * Deriv(1, 0); //dvdy*dudx
189  val += Dep(_XY_, _XY_) * Deriv(1, 1); //dvdy*dudy
190  val *= 0.5;
191  ek(in * nstate + 0, jn * nstate + 0) += weight * val;
192 
193  val = Dep(_XX_, _XY_) * Deriv(0, 0);
194  val += 2. * Dep(_XX_, _YY_) * Deriv(0, 1);
195  val += Dep(_XY_, _XY_) * Deriv(1, 0);
196  val += 2. * Dep(_XY_, _YY_) * Deriv(1, 1);
197  val *= 0.5;
198  ek(in * nstate + 0, jn * nstate + 1) += weight * val;
199 
200  val = 2. * Dep(_XY_, _XX_) * Deriv(0, 0);
201  val += Dep(_XY_, _XY_) * Deriv(0, 1);
202  val += 2. * Dep(_YY_, _XX_) * Deriv(1, 0);
203  val += Dep(_YY_, _XY_) * Deriv(1, 1);
204  val *= 0.5;
205  ek(in * nstate + 1, jn * nstate + 0) += weight * val;
206 
207  val = Dep(_XY_, _XY_) * Deriv(0, 0);
208  val += 2. * Dep(_XY_, _YY_) * Deriv(0, 1);
209  val += Dep(_YY_, _XY_) * Deriv(1, 0);
210  val += 2. * Dep(_YY_, _YY_) * Deriv(1, 1);
211  val *= 0.5;
212  ek(in * nstate + 1, jn * nstate + 1) += weight * val;
213  }
214  }
215 
216 #ifdef LOG4CXX
217  if (elastoplastic2dLogger->isDebugEnabled()) {
218  std::stringstream sout;
219  sout << "<<< TPZMatElastoPlastic2D<T,TMEM>::Contribute ***";
220  sout << " Resultant rhs vector:\n" << ef;
221  sout << " Resultant stiff vector:\n" << ek;
222  LOGPZ_DEBUG(elastoplastic2dLogger, sout.str().c_str());
223  }
224 #endif
225 }
226 
227 template <class T, class TMEM>
228 void TPZMatElastoPlastic2D<T, TMEM>::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<REAL> &ef) {
229  TPZFMatrix<REAL> &dphi = data.dphix;
230  TPZFMatrix<REAL> &phi = data.phi;
231  TPZFMatrix<REAL> &axes = data.axes;
232  TPZFNMatrix<9, REAL> axesT;
233  TPZFNMatrix<50, REAL> dphiXY;
234 
235  axes.Transpose(&axesT);
236  axesT.Multiply(dphi, dphiXY);
237 
238  const int phr = phi.Rows();
239 
240  //TPZFNMatrix<36> Deriv(6, 6);
241  TPZFNMatrix<6> DeltaStrain(6, 1);
242  TPZFNMatrix<6> Stress(6, 1);
243  int ptindex = data.intGlobPtIndex;
244 
245  if (TPZMatWithMem<TMEM>::fUpdateMem && data.sol.size() > 0) {
246 
247  TPZSolVec locsol(data.sol);
248  TPZGradSolVec locdsol(data.dsol);
249  int numsol = locsol.size();
250 
251  for (int is = 0; is < numsol; is++) {
252  data.sol[0] = locsol[is];
253  data.dsol[0] = locdsol[is];
254 
255  this->ComputeDeltaStrainVector(data, DeltaStrain);
256  this->ApplyDeltaStrain(data, DeltaStrain, Stress);
257  }
258  } else {
259  this->ComputeDeltaStrainVector(data, DeltaStrain);
260  this->ApplyDeltaStrain(data, DeltaStrain, Stress);
261  }
262 #ifdef MACOS
263  feclearexcept(FE_ALL_EXCEPT);
264  if (fetestexcept(/*FE_DIVBYZERO*/ FE_ALL_EXCEPT)) {
265  std::cout << "division by zero reported\n";
266  DebugStop();
267  }
268 #endif
269 
270 #ifdef LOG4CXX
271  if (elastoplastic2dLogger->isDebugEnabled()) {
272  std::stringstream sout;
273  sout << ">>> TPZMatElastoPlastic<T,TMEM>::Contribute ***";
274  sout << "\nIntegration Local Point index = " << data.intGlobPtIndex;
275  sout << "\nIntegration Global Point index = " << data.intGlobPtIndex;
276  sout << "\ndata.axes = " << data.axes;
277  sout << "\nStress " << endl;
278  sout << Stress(_XX_, 0) << "\t" << Stress(_YY_, 0) << "\t" << Stress(_XY_, 0) << "\n";
279  sout << "\nDELTA STRAIN " << endl;
280  sout << DeltaStrain(0, 0) << "\t" << DeltaStrain(1, 0) << "\t" << DeltaStrain(2, 0) << "\n";
281  sout << "data.phi" << data.phi<< std::endl;
282  TPZMatWithMem<TMEM>::MemItem(ptindex).Print(sout);
283 
284  LOGPZ_DEBUG(elastoplastic2dLogger, sout.str().c_str());
285  }
286 #endif
287  int nstate = NStateVariables();
288  REAL val;
289 
290  TPZManVector<STATE, 2> ForceLoc(nstate,0.0);
291  if (this->fForcingFunction) {
292  this->fForcingFunction->Execute(data.x, ForceLoc);
293  }
294 
295  int in;
296  for (in = 0; in < phr; in++) {
297  val = ForceLoc[0] * phi(in, 0);
298  val -= Stress(_XX_, 0) * dphiXY(0, in);
299  val -= Stress(_XY_, 0) * dphiXY(1, in);
300  ef(in * nstate + 0, 0) += weight * val;
301 
302  val = ForceLoc[1] * phi(in, 0);
303  val -= Stress(_XY_, 0) * dphiXY(0, in);
304  val -= Stress(_YY_, 0) * dphiXY(1, in);
305  ef(in * nstate + 1, 0) += weight * val;
306  }
307 
308 
309 #ifdef LOG4CXX
310  if (elastoplastic2dLogger->isDebugEnabled()) {
311  std::stringstream sout;
312  sout << "<<< TPZMatElastoPlastic2D<T,TMEM>::Contribute ***";
313  sout << " Resultant rhs vector:\n" << ef;
314  LOGPZ_DEBUG(elastoplastic2dLogger, sout.str().c_str());
315  }
316 #endif
317 
318 }
319 
320 
321 template <class T, class TMEM>
322 void TPZMatElastoPlastic2D<T,TMEM>::FillBoundaryConditionDataRequirement(int type,TPZMaterialData &data)
323 {
324 
325  data.fNeedsSol = true;
326  if (type == 4 || type ==5 || type == 6) {
327  data.fNeedsNormal = true;
328  }
329  else {
330  data.fNeedsNormal = false;
331  }
332 }
333 
334 template <class T, class TMEM>
335 void TPZMatElastoPlastic2D<T, TMEM>::ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix<REAL> &ek, TPZFMatrix<REAL> &ef, TPZBndCond &bc) {
336  int nstate = NStateVariables();
337  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
338 
339  TPZBndCondWithMem<TMEM> & bc_with_memory = dynamic_cast<TPZBndCondWithMem<TMEM> &>(bc);
340 
342  int gp_index = data.intGlobPtIndex;
343  TPZFMatrix<REAL> &phi = data.phi;
344  TPZManVector<STATE,3> delta_u = data.sol[0];
345  TPZManVector<STATE,3> u_n(nstate,0.0);
346  TPZManVector<STATE,3> u(bc_with_memory.MemItem(gp_index).m_u);
347  for (int i = 0; i < nstate; i++) {
348  u_n[i] = delta_u[i] + u[i];
349  }
350 
351  if(TPZMatWithMem<TMEM>::fUpdateMem)
352  {
353  bc_with_memory.MemItem(gp_index).m_u = u_n;
354  }
355 
356  const int phr = phi.Rows();
357  int in, jn, idf, jdf;
358  REAL v2[2];
359  v2[0] = bc.Val2()(0, 0);
360  v2[1] = bc.Val2()(1, 0);
361 
362  TPZFMatrix<REAL> &v1 = bc.Val1();
363  switch (bc.Type()) {
364  case 0: // Dirichlet condition
365  for (in = 0; in < phr; in++) {
366  ef(nstate * in + 0, 0) += BIGNUMBER * (v2[0] - u_n[0]) * phi(in, 0) * weight;
367  ef(nstate * in + 1, 0) += BIGNUMBER * (v2[1] - u_n[1]) * phi(in, 0) * weight;
368 
369  for (jn = 0; jn < phr; jn++) {
370  ek(nstate * in + 0, nstate * jn + 0) += BIGNUMBER * phi(in, 0) * phi(jn, 0) * weight;
371  ek(nstate * in + 1, nstate * jn + 1) += BIGNUMBER * phi(in, 0) * phi(jn, 0) * weight;
372 
373  }//jn
374  }//in
375  break;
376 
377  case 1: // Neumann condition
378  for (in = 0; in < phi.Rows(); in++) {
379  ef(nstate * in + 0, 0) += v2[0] * phi(in, 0) * weight;
380  ef(nstate * in + 1, 0) += v2[1] * phi(in, 0) * weight;
381  }
382  break;
383 
384  case 2: // Mixed condition
385  {
386  TPZFNMatrix<2, STATE> res(2, 1, 0.);
387  for (int i = 0; i < 2; i++) {
388  for (int j = 0; j < 2; j++) {
389  res(i, 0) += bc.Val1()(i, j) * u_n[j];
390  }
391  }
392 
393  for (in = 0; in < phi.Rows(); in++) {
394  ef(nstate * in + 0, 0) += (v2[0] - res(0, 0)) * phi(in, 0) * weight;
395  ef(nstate * in + 1, 0) += (v2[1] - res(1, 0)) * phi(in, 0) * weight;
396  for (jn = 0; jn < phi.Rows(); jn++) {
397  for (idf = 0; idf < 2; idf++) {
398  for (jdf = 0; jdf < 2; jdf++) {
399  ek(nstate * in + idf, nstate * jn + jdf) += bc.Val1()(idf, jdf) * phi(in, 0) * phi(jn, 0) * weight;
400  //BUG FALTA COLOCAR VAL2
401  //DebugStop();
402  }
403  }
404  }
405  }//in
406  }
407  break;
408 
409  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
410  for (in = 0; in < phr; in++) {
411  ef(nstate * in + 0, 0) += BIGNUMBER * (0. - u_n[0]) * v2[0] * phi(in, 0) * weight;
412  ef(nstate * in + 1, 0) += BIGNUMBER * (0. - u_n[1]) * v2[1] * phi(in, 0) * weight;
413  for (jn = 0; jn < phr; jn++) {
414  ek(nstate * in + 0, nstate * jn + 0) += BIGNUMBER * phi(in, 0) * phi(jn, 0) * weight * v2[0];
415  ek(nstate * in + 1, nstate * jn + 1) += BIGNUMBER * phi(in, 0) * phi(jn, 0) * weight * v2[1];
416  }//jn
417  }//in
418  break;
419 
420  case 4: // stressField Neumann condition
421  v2[0] = v1(0, 0) * data.normal[0] + v1(0, 1) * data.normal[1];
422  v2[1] = v1(1, 0) * data.normal[0] + v1(1, 1) * data.normal[1];
423  // The normal vector points towards the neighbor. The negative sign is there to
424  // reflect the outward normal vector.
425  for (in = 0; in < phi.Rows(); in++) {
426  ef(nstate * in + 0, 0) += v2[0] * phi(in, 0) * weight;
427  ef(nstate * in + 1, 0) += v2[1] * phi(in, 0) * weight;
428  }
429  break;
430  case 5://PRESSAO DEVE SER POSTA NA POSICAO 0 DO VETOR v2
431  {
432  TPZFNMatrix<2, STATE> res(2, 1, 0.);
433  for (int i = 0; i < 2; i++) {
434  for (int j = 0; j < 2; j++) {
435  res(i, 0) += data.normal[i] * bc.Val1()(i, j) * u_n[j] * data.normal[j];
436  }
437  }
438  for (in = 0; in < phi.Rows(); in++) {
439  ef(nstate * in + 0, 0) += (v2[0] * data.normal[0] - res(0, 0)) * phi(in, 0) * weight;
440  ef(nstate * in + 1, 0) += (v2[0] * data.normal[1] - res(1, 0)) * phi(in, 0) * weight;
441  for (jn = 0; jn < phi.Rows(); jn++) {
442  for (idf = 0; idf < 2; idf++) {
443  for (jdf = 0; jdf < 2; jdf++) {
444  ek(nstate * in + idf, nstate * jn + jdf) += bc.Val1()(idf, jdf) * data.normal[idf] * data.normal[jdf] * phi(in, 0) * phi(jn, 0) * weight;
445  // BUG FALTA COLOCAR VAL2
446  // DebugStop();
447  }
448  }
449  }
450 
451  }
452  }
453  break;
454 
455  case 6:
456  {
457 
458  REAL v[1];
459  v[0] = bc.Val2()(0,0); // Tn normal component of normal traction (T)
460 
461  REAL tn = v[0];
462  TPZManVector<REAL,3> n = data.normal;
465  for(in = 0 ; in < phi.Rows(); in++)
466  {
468  ef(nstate*in+0,0) += weight * tn * n[0] * phi(in,0); // Tnx
469  ef(nstate*in+1,0) += weight * tn * n[1] * phi(in,0); // Tny
470  }
471 
472  }
473  break;
474 
475  case 7 : {
476 
477  REAL v_null[2];
478  v_null[0] = bc.Val1()(0, 0);
479  v_null[1] = bc.Val1()(1, 1);
480 
481  for (in = 0; in < phr; in++) {
482  ef(nstate * in + 0, 0) += BIGNUMBER * (v2[0] - u_n[0]) * v_null[0] * phi(in, 0) * weight;
483  ef(nstate * in + 1, 0) += BIGNUMBER * (v2[1] - u_n[1]) * v_null[1] * phi(in, 0) * weight;
484  for (jn = 0; jn < phr; jn++) {
485  ek(nstate * in + 0, nstate * jn + 0) += BIGNUMBER * phi(in, 0) * phi(jn, 0) * weight * v_null[0];
486  ek(nstate * in + 1, nstate * jn + 1) += BIGNUMBER * phi(in, 0) * phi(jn, 0) * weight * v_null[1];
487  }//jn
488  }//in
489  break;
490 
491  }
492 
493  default:
494 #ifdef LOG4CXX
495  if (elastoplastic2dLogger->isDebugEnabled()) {
496  std::stringstream sout;
497  sout << "<<< TPZMatElastoPlastic2D<T,TMEM>::ContributeBC *** WRONG BOUNDARY CONDITION TYPE = " << bc.Type();
498  LOGPZ_ERROR(elastoplastic2dLogger, sout.str().c_str());
499  }
500 #endif
501  PZError << "TPZMatElastoPlastic2D::ContributeBC error - Wrong boundary condition type" << std::endl;
502  }
503 }
504 
505 template <class T, class TMEM>
506 void TPZMatElastoPlastic2D<T, TMEM>::ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ef, TPZBndCond &bc) {
507  TPZFMatrix<STATE> ek_fake(ef.Rows(),ef.Rows());
508  int nstate = NStateVariables();
509  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
510 
511  TPZBndCondWithMem<TMEM> * bc_with_memory = dynamic_cast<TPZBndCondWithMem<TMEM> *>(&bc);
512  if(!bc_with_memory){
513  PZError << "TPZMatElastoPlastic2D::ContributeBC error - Wrong boundary objected: expected TPZBndCondWithMem<TMEM>" << std::endl;
514  DebugStop();
515  }
516 
518  int gp_index = data.intGlobPtIndex;
519  TPZFMatrix<REAL> &phi = data.phi;
520  TPZManVector<STATE,3> delta_u = data.sol[0];
521  TPZManVector<STATE,3> u_n(nstate,0.0);
522  TPZManVector<STATE,3> u(bc_with_memory->MemItem(gp_index).m_u);
523  for (int i = 0; i < nstate; i++) {
524  u_n[i] = delta_u[i] + u[i];
525  }
526 
527  if(TPZMatWithMem<TMEM>::fUpdateMem)
528  {
529  bc_with_memory->MemItem(gp_index).m_u = u_n;
530  }
531 
532  const int phr = phi.Rows();
533  int in, jn, idf, jdf;
534  REAL v2[2];
535  v2[0] = bc.Val2()(0, 0);
536  v2[1] = bc.Val2()(1, 0);
537 
538  TPZFMatrix<REAL> &v1 = bc.Val1();
539  switch (bc.Type()) {
540  case 0: // Dirichlet condition
541  for (in = 0; in < phr; in++) {
542  ef(nstate * in + 0, 0) += BIGNUMBER * (v2[0] - u_n[0]) * phi(in, 0) * weight;
543  ef(nstate * in + 1, 0) += BIGNUMBER * (v2[1] - u_n[1]) * phi(in, 0) * weight;
544 
545  }//in
546  break;
547 
548  case 1: // Neumann condition
549  for (in = 0; in < phi.Rows(); in++) {
550  ef(nstate * in + 0, 0) += v2[0] * phi(in, 0) * weight;
551  ef(nstate * in + 1, 0) += v2[1] * phi(in, 0) * weight;
552  }
553  break;
554 
555  case 2: // Mixed condition
556  {
557  TPZFNMatrix<2, STATE> res(2, 1, 0.);
558  for (int i = 0; i < 2; i++) {
559  for (int j = 0; j < 2; j++) {
560  res(i, 0) += bc.Val1()(i, j) * u_n[j];
561  }
562  }
563 
564  for (in = 0; in < phi.Rows(); in++) {
565  ef(nstate * in + 0, 0) += (v2[0] - res(0, 0)) * phi(in, 0) * weight;
566  ef(nstate * in + 1, 0) += (v2[1] - res(1, 0)) * phi(in, 0) * weight;
567  }//in
568  }
569  break;
570 
571  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
572  for (in = 0; in < phr; in++) {
573  ef(nstate * in + 0, 0) += BIGNUMBER * (0. - u_n[0]) * v2[0] * phi(in, 0) * weight;
574  ef(nstate * in + 1, 0) += BIGNUMBER * (0. - u_n[1]) * v2[1] * phi(in, 0) * weight;
575  }//in
576  break;
577 
578  case 4: // stressField Neumann condition
579  v2[0] = v1(0, 0) * data.normal[0] + v1(0, 1) * data.normal[1];
580  v2[1] = v1(1, 0) * data.normal[0] + v1(1, 1) * data.normal[1];
581  // The normal vector points towards the neighbor. The negative sign is there to
582  // reflect the outward normal vector.
583  for (in = 0; in < phi.Rows(); in++) {
584  ef(nstate * in + 0, 0) += v2[0] * phi(in, 0) * weight;
585  ef(nstate * in + 1, 0) += v2[1] * phi(in, 0) * weight;
586  }
587  break;
588  case 5://PRESSAO DEVE SER POSTA NA POSICAO 0 DO VETOR v2
589  {
590  TPZFNMatrix<2, STATE> res(2, 1, 0.);
591  for (int i = 0; i < 2; i++) {
592  for (int j = 0; j < 2; j++) {
593  res(i, 0) += data.normal[i] * bc.Val1()(i, j) * u_n[j] * data.normal[j];
594  }
595  }
596  for (in = 0; in < phi.Rows(); in++) {
597  ef(nstate * in + 0, 0) += (v2[0] * data.normal[0] - res(0, 0)) * phi(in, 0) * weight;
598  ef(nstate * in + 1, 0) += (v2[0] * data.normal[1] - res(1, 0)) * phi(in, 0) * weight;
599  }
600  }
601  break;
602 
603  case 6:
604  {
605 
606  REAL v[1];
607  v[0] = bc.Val2()(0,0); // Tn normal component of normal traction (T)
608 
609  REAL tn = v[0];
610  TPZManVector<REAL,3> n = data.normal;
613  for(in = 0 ; in < phi.Rows(); in++)
614  {
616  ef(nstate*in+0,0) += weight * tn * n[0] * phi(in,0); // Tnx
617  ef(nstate*in+1,0) += weight * tn * n[1] * phi(in,0); // Tny
618  }
619 
620  }
621  break;
622 
623  case 7 : {
624 
625  REAL v_null[2];
626  v_null[0] = bc.Val1()(0, 0);
627  v_null[1] = bc.Val1()(1, 1);
628 
629  for (in = 0; in < phr; in++) {
630  ef(nstate * in + 0, 0) += BIGNUMBER * (v2[0] - u_n[0]) * v_null[0] * phi(in, 0) * weight;
631  ef(nstate * in + 1, 0) += BIGNUMBER * (v2[1] - u_n[1]) * v_null[1] * phi(in, 0) * weight;
632  }//in
633  break;
634 
635  }
636 
637  default:
638 #ifdef LOG4CXX
639  if (elastoplastic2dLogger->isDebugEnabled()) {
640  std::stringstream sout;
641  sout << "<<< TPZMatElastoPlastic2D<T,TMEM>::ContributeBC *** WRONG BOUNDARY CONDITION TYPE = " << bc.Type();
642  LOGPZ_ERROR(elastoplastic2dLogger, sout.str().c_str());
643  }
644 #endif
645  PZError << "TPZMatElastoPlastic2D::ContributeBC error - Wrong boundary condition type" << std::endl;
646  }
647 }
648 
649 
650 
651 template <class T, class TMEM>
652 void TPZMatElastoPlastic2D<T,TMEM>::Solution(TPZMaterialData &data, int var, TPZVec<REAL> &Solout)
653 {
654 
655  TPZMaterialData datalocal(data);
656  datalocal.sol[0].Resize(3,0.);
657  datalocal.dsol[0].Resize(3,3);
658  datalocal.dsol[0](2,0) = 0.;
659  datalocal.dsol[0](2,1) = 0.;
660  datalocal.dsol[0](2,2) = 0.;
661  datalocal.dsol[0](0,2) = 0.;
662  datalocal.dsol[0](1,2) = 0.;
663  TPZMatElastoPlastic<T,TMEM>::Solution(datalocal,var,Solout);
664 
665 }
666 
667 template <class T, class TMEM>
668 void TPZMatElastoPlastic2D<T, TMEM>::ComputeDeltaStrainVector(TPZMaterialData & data, TPZFMatrix<REAL> &DeltaStrain) {
669  TPZFNMatrix<9> DSolXYZ(3, 3, 0.);
670  data.axes.Multiply(data.dsol[0], DSolXYZ, 1/*transpose*/);
671  if (DeltaStrain.Rows() != 6) {
672  DebugStop();
673  }
674 
675  DeltaStrain(_XX_, 0) = DSolXYZ(0, 0);
676  DeltaStrain(_YY_, 0) = DSolXYZ(1, 1);
677  DeltaStrain(_XY_, 0) = 0.5 * (DSolXYZ(1, 0) + DSolXYZ(0, 1));
678  DeltaStrain(_XZ_, 0) = 0.;
679  DeltaStrain(_YZ_, 0) = 0.;
680  DeltaStrain(_ZZ_, 0) = 0.;
681 
682 }
683 
684 
685 template <class T, class TMEM>
686 TPZMaterial * TPZMatElastoPlastic2D<T,TMEM>::NewMaterial()
687 {
688  return new TPZMatElastoPlastic2D<T,TMEM>(*this);
689 }
690 
691 #include "TPZSandlerExtended.h"
692 #include "TPZPlasticStepPV.h"
693 #include "TPZYCMohrCoulombPV.h"
694 
695 template <class T, class TMEM>
696 std::string TPZMatElastoPlastic2D<T,TMEM>::Name()
697 {
698  return "TPZMatElastoPlastic<T,TMEM>";
699 }
700 
701 template <class T, class TMEM>
702 void TPZMatElastoPlastic2D<T,TMEM>::Write(TPZStream &buf, int withclassid) const{
703  TPZMatElastoPlastic<T,TMEM>::Write(buf,withclassid);
704  int classid = ClassId();
705  buf.Write(&classid);
706 }
707 
708 template <class T, class TMEM>
709 void TPZMatElastoPlastic2D<T,TMEM>::Read(TPZStream &buf, void *context)
710 {
711  TPZMatElastoPlastic<T,TMEM>::Read(buf,context);
712  int classid;
713  buf.Read(&classid);
714  if (classid != ClassId()) {
715  DebugStop();
716  }
717 }
718 
719 
720 template <class T, class TMEM>
721 void TPZMatElastoPlastic2D<T,TMEM>::Print(std::ostream &out, const int memory)
722 {
723  out << __PRETTY_FUNCTION__ << std::endl;
724  TPZMatElastoPlastic<T,TMEM>::Print(out,memory);
725 }
726 
727 template <class T, class TMEM>
728 void TPZMatElastoPlastic2D<T,TMEM>::Print(std::ostream &out)
729 {
730  out << __PRETTY_FUNCTION__ << std::endl;
731  out << "Plane strain " << fPlaneStrain << std::endl;
732  TPZMatElastoPlastic<T,TMEM>::Print(out);
733 }
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
_XZ_
#define _XZ_
Definition: TPZTensor.h:29
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
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.
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZMatElastoPlastic2D::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
TPZMatElastoPlastic::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Definition: TPZMatElastoPlastic_impl.h:1175
TPZSandlerExtended.h
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZYCMohrCoulombPV.h
TPZMatElastoPlastic2D::ComputeDeltaStrainVector
virtual void ComputeDeltaStrainVector(TPZMaterialData &data, TPZFMatrix< REAL > &DeltaStrain)
Definition: TPZMatElastoPlastic2D_impl.h:668
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZMatElastoPlastic2D::FillBoundaryConditionDataRequirement
virtual void FillBoundaryConditionDataRequirement(int type, TPZMaterialData &data) override
Definition: TPZMatElastoPlastic2D_impl.h:322
TPZMatElastoPlastic::ApplyDeltaStrain
void ApplyDeltaStrain(TPZMaterialData &data, TPZFMatrix< REAL > &DeltaStrain, TPZFMatrix< REAL > &Stress)
Definition: TPZMatElastoPlastic_impl.h:1059
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMatElastoPlastic::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< REAL > &Solout) override
Definition: TPZMatElastoPlastic_impl.h:243
TPZMatElastoPlastic2D::Read
virtual void Read(TPZStream &buf, void *context) override
Definition: TPZMatElastoPlastic2D_impl.h:709
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
TPZMatWithMem
Implements an abstract class implementing the memory features.
Definition: TPZMatWithMem.h:23
TPZPlasticStepPV.h
TPZBndCond::Val2
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZMaterialData::phi
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
Definition: pzmaterialdata.h:55
TPZMatElastoPlastic2D.h
TPZMatElastoPlastic::Read
virtual void Read(TPZStream &buf, void *context) override
Definition: TPZMatElastoPlastic_impl.h:1185
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZMaterial
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZMatElastoPlastic2D::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Definition: TPZMatElastoPlastic2D_impl.h:702
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZMatElastoPlastic::Print
virtual void Print(std::ostream &out, const int memory)
Definition: TPZMatElastoPlastic_impl.h:145
LOGPZ_INFO
#define LOGPZ_INFO(A, B)
Define log for informations.
Definition: pzlog.h:89
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
_XX_
#define _XX_
Definition: TPZTensor.h:27
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
_YZ_
#define _YZ_
Definition: TPZTensor.h:31
_XY_
#define _XY_
Definition: TPZTensor.h:28
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZMatElastoPlastic2D::ClassId
virtual int ClassId() const override
Definition: TPZMatElastoPlastic2D.h:166
TPZBndCondWithMem.h
TPZMatElastoPlastic2D::ApplyDeltaStrainComputeDep
virtual void ApplyDeltaStrainComputeDep(TPZMaterialData &data, TPZFMatrix< REAL > &DeltaStrain, TPZFMatrix< REAL > &Stress, TPZFMatrix< REAL > &Dep)
Definition: TPZMatElastoPlastic2D_impl.h:77
TPZMaterialData::axes
TPZFNMatrix< 9, REAL > axes
axes indicating the directions of the derivatives of the shapefunctions
Definition: pzmaterialdata.h:63
TPZMatElastoPlastic2D::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< REAL > &Solout) override
Definition: TPZMatElastoPlastic2D_impl.h:652
test.res
string res
Definition: test.py:151
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZMaterialData::intGlobPtIndex
int intGlobPtIndex
global point index
Definition: pzmaterialdata.h:115
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
_YY_
#define _YY_
Definition: TPZTensor.h:30
TPZMaterial::gBigNumber
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
TPZMatElastoPlastic2D::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< REAL > &ef) override
Definition: TPZMatElastoPlastic2D_impl.h:228
_ZZ_
#define _ZZ_
Definition: TPZTensor.h:32
TPZMatElastoPlastic2D::ApplyDeltaStrain
virtual void ApplyDeltaStrain(TPZMaterialData &data, TPZFMatrix< REAL > &DeltaStrain, TPZFMatrix< REAL > &Stress)
Definition: TPZMatElastoPlastic2D_impl.h:56
TPZMatElastoPlastic2D::Name
virtual std::string Name() override
Definition: TPZMatElastoPlastic2D_impl.h:696
idf
static int idf[6]
Definition: TPZCompMeshTools.cpp:135
TPZMatElastoPlastic2D::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: TPZMatElastoPlastic2D_impl.h:728
TPZMatElastoPlastic2D::NStateVariables
virtual int NStateVariables() const override
Definition: TPZMatElastoPlastic2D.h:54
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
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZMatElastoPlastic_impl.h
TPZMatWithMem::MemItem
virtual TMEM & MemItem(const int i) const
Definition: TPZMatWithMem.h:172
TPZMatElastoPlastic::ApplyDeltaStrainComputeDep
void ApplyDeltaStrainComputeDep(TPZMaterialData &data, TPZFMatrix< REAL > &DeltaStrain, TPZFMatrix< REAL > &Stress, TPZFMatrix< REAL > &Dep)
Definition: TPZMatElastoPlastic_impl.h:1011
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZFMatrix::Transpose
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
TPZBndCondWithMem::MemItem
virtual TMEM & MemItem(const int i) const
Definition: TPZBndCondWithMem_impl.h:18
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
TPZMatElastoPlastic2D::NewMaterial
virtual TPZMaterial * NewMaterial() override
Definition: TPZMatElastoPlastic2D_impl.h:686
Generated on Sun Aug 16 2020 11:17:02 for NeoPZ by   doxygen 1.8.13