NeoPZ
pzelast3d.cpp
Go to the documentation of this file.
1 
6 #include "pzelast3d.h"
7 #include "pzbndcond.h"
8 #include "pzmatrix.h"
9 #include "pzfmatrix.h"
10 #include "pzerror.h"
11 #include "pzmanvector.h"
12 #include "pzaxestools.h"
13 #include <math.h>
14 #include <fstream>
15 
16 #define CODE3
17 
18 
19 STATE TPZElasticity3D::gTolerance = 1.e-11;
20 
21 TPZElasticity3D::TPZElasticity3D(int nummat, STATE E, STATE poisson, TPZVec<STATE> &force,
22  STATE preStressXX, STATE preStressYY, STATE preStressZZ) :
23 TPZRegisterClassId(&TPZElasticity3D::ClassId),
24  TPZMaterial(nummat),fFy(0.),fFrictionAngle(0.),fCohesion(0.),fPlasticPostProc(ENonePlasticProc)
25 {
26  this->fE = E;
27  this->fPoisson = poisson;
28 #ifdef PZDEBUG
29  if (force.NElements() != 3) PZError << __PRETTY_FUNCTION__ << " - error!" << std::endl;
30 #endif
31  int i;
32  this->fForce.Resize(3);
33  for(i = 0; i < 3; i++) this->fForce[i] = force[i];
34  //Default directions is {1,0,0}
35  this->fPostProcessDirection.Resize(3);
36  this->fPostProcessDirection.Fill(0.);
37  this->fPostProcessDirection[0] = 1.;
38  SetC();
39 
40  fPreStress.Resize(3);
41  fPreStress[0] = preStressXX;
42  fPreStress[1] = preStressYY;
43  fPreStress[2] = preStressZZ;
44 
45 }//method
46 
47 TPZElasticity3D::TPZElasticity3D(int nummat) :
48 TPZRegisterClassId(&TPZElasticity3D::ClassId),
49 TPZMaterial(nummat), fE(0.), fPoisson(0.),
50  fForce(3,0.),
51  fPostProcessDirection(3,0.), fFy(0.), fPreStress(3,0.),
52  fFrictionAngle(0.),fCohesion(0.),fPlasticPostProc(ENonePlasticProc)
53 {
54  SetC();
55 }
56 
57 TPZElasticity3D::TPZElasticity3D() : TPZRegisterClassId(&TPZElasticity3D::ClassId),
58 TPZMaterial(),fE(0.), fPoisson(0.), C1(-999.), C2(-999.), C3(-999.),
59  fForce(3,0.), fPostProcessDirection(3,0.), fFy(0.), fPreStress(3,0.),
60  fFrictionAngle(0.),fCohesion(0.),fPlasticPostProc(ENonePlasticProc)
61 {
62 }
63 
64 TPZElasticity3D::~TPZElasticity3D(){}
65 
66 TPZElasticity3D::TPZElasticity3D(const TPZElasticity3D &cp) :
67 TPZRegisterClassId(&TPZElasticity3D::ClassId),
68 TPZMaterial(cp), fE(cp.fE), fPoisson(cp.fPoisson),
69  fForce(cp.fForce),
70  fPostProcessDirection(cp.fPostProcessDirection), fFy(cp.fFy),
71  fPreStress(cp.fPreStress),
72  fFrictionAngle(cp.fFrictionAngle),fCohesion(cp.fCohesion),fPlasticPostProc(cp.fPlasticPostProc)
73 {
74  SetC();
75 }
76 
77 void TPZElasticity3D::Print(std::ostream & out){
78  out << "\nTPZElasticity3D material:\n";
79  out << "\tfE = " << this->fE << std::endl;
80  out << "\tfPoisson = " << this->fPoisson << std::endl;
81  out << "\tfForce = " << this->fForce << std::endl;
82 // DebugStop(); //TODO fFrictionAngle(0.),fCohesion(0.),fPlasticPostProc(
83  out << "\tBase class print\n";
84  TPZMaterial::Print(out);
85  out << "End of TPZElasticity3D::Print\n";
86 }
87 
88 void TPZElasticity3D::Contribute(TPZMaterialData &data,
89  REAL weight,
90  TPZFMatrix<STATE> &ek,
91  TPZFMatrix<STATE> &ef){
92 
93  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
94  if(shapetype == data.EVecShape){
95  ContributeVecShape(data,weight,ek,ef);
96  return;
97  }
98 
99 // TPZFMatrix<REAL> &dphi = data.dphix;
100  TPZFMatrix<REAL> &phi = data.phi;
101  TPZManVector<REAL,3> &x = data.x;
102 
103  TPZFNMatrix<666,REAL> dphi(3,data.dphix.Cols());
104 
105  TPZAxesTools<REAL>::Axes2XYZ(data.dphix,dphi,data.axes);
106  const int phr = phi.Rows();
107  TPZManVector<STATE,3> locForce(fForce);
108  if(this->fForcingFunction){
109  this->fForcingFunction->Execute(x,locForce);
110  }
111 #ifdef CODE0
112  TPZFNMatrix<9> Deriv(3,3);
113  const REAL E = this->fE;
114  const REAL nu = this->fPoisson;
115  const REAL C1 = E / (2.+ 2.*nu);
116  const REAL C2 = E * nu / (-1. + nu + 2.*nu*nu);
117  const REAL C3 = E * (nu - 1.) / (-1. + nu +2. * nu * nu);
118 
119  int in;
120  for(in = 0; in < phr; in++)
121  {
122  int kd;
123  for(kd = 0; kd < 3; kd++)
124  {
125  ef(in*3+kd, 0) += weight * ( locForce[kd] * phi(in,0) - fPreStress[kd] * dphi(kd,in) );
126  }//kd
127  REAL val;
128  for( int jn = 0; jn < phr; jn++ )
129  {
130  //Compute Deriv matrix
131  for(int ud = 0; ud < 3; ud++)
132  {
133  for(int vd = 0; vd < 3; vd++)
134  {
135  Deriv(vd,ud) = dphi(vd,in)*dphi(ud,jn);
136  }//ud
137  }//vd
138 
139  //First equation Dot[Sigma1, gradV1]
140  val = ( Deriv(1,1) + Deriv(2,2) ) * C1 + Deriv(0,0) * C3;
141  ek(in*3+0,jn*3+0) += weight * val;
142 
143  val = Deriv(1,0) * C1 - Deriv(0,1) * C2;
144  ek(in*3+0,jn*3+1) += weight * val;
145 
146  val = Deriv(2,0) * C1 - Deriv(0,2) * C2;
147  ek(in*3+0,jn*3+2) += weight * val;
148 
149  //Second equation Dot[Sigma2, gradV2]
150  val = Deriv(0,1) * C1 - Deriv(1,0) * C2;
151  ek(in*3+1,jn*3+0) += weight * val;
152 
153  val = ( Deriv(0,0) + Deriv(2,2) ) * C1 + Deriv(1,1) * C3;
154  ek(in*3+1,jn*3+1) += weight * val;
155 
156  val = Deriv(2,1) * C1 - Deriv(1,2) * C2;
157  ek(in*3+1,jn*3+2) += weight * val;
158 
159  //Third equation Dot[Sigma3, gradV3]
160  val = Deriv(0,2) * C1 - Deriv(2,0) * C2;
161  ek(in*3+2,jn*3+0) += weight * val;
162 
163  val = Deriv(1,2) * C1 - Deriv(2,1) * C2;
164  ek(in*3+2,jn*3+1) += weight * val;
165 
166  val = ( Deriv(0,0) + Deriv(1,1) ) * C1 + Deriv(2,2) * C3;
167  ek(in*3+2,jn*3+2) += weight * val;
168 
169  }//jn
170  }//in
171 #endif
172 #ifdef CODE1
173  REAL Deriv[3][3];
174 
175  for(int jn = 0; jn < phr; jn++)
176  {
177  int kd;
178  for(kd = 0; kd < 3; kd++)
179  {
180  ef(jn*3+kd, 0) += weight * ( locForce[kd] * phi(jn,0) - fPreStress[kd] * dphi(kd,jn) );
181  }//kd
182  for( int in = 0; in < phr; in++ )
183  {
184  //Compute Deriv matrix
185  for(int ud = 0; ud < 3; ud++)
186  {
187  for(int vd = 0; vd < 3; vd++)
188  {
189  Deriv[vd][ud] = dphi(vd,in)*dphi(ud,jn);
190  }//ud
191  }//vd
192 
193  //First equation Dot[Sigma1, gradV1]
194  STATE *ptr1 = &ek(in*3,jn*3);
195  /*ek(in*3+0,jn*3+0)*/ *ptr1++ += weight * (( Deriv[1][1] + Deriv[2][2] ) * C1 + Deriv[0][0] * C3);
196 
197  //Second equation Dot[Sigma2, gradV2]
198  /*ek(in*3+1,jn*3+0)*/ *ptr1++ += weight * (Deriv[0][1] * C1 - Deriv[1][0] * C2);
199 
200  //Third equation Dot[Sigma3, gradV3]
201  /*ek(in*3+2,jn*3+0)*/ *ptr1 += weight * (Deriv[0][2] * C1 - Deriv[2][0] * C2);
202 
203  STATE *ptr2 = &ek(in*3+0, jn*3+1);
204  /*ek(in*3+0,jn*3+1)*/ *ptr2++ += weight * (Deriv[1][0] * C1 - Deriv[0][1] * C2);
205 
206  /*ek(in*3+1,jn*3+1)*/ *ptr2++ += weight * (( Deriv[0][0] + Deriv[2][2] ) * C1 + Deriv[1][1] * C3);
207 
208  /*ek(in*3+2,jn*3+1)*/ *ptr2 += weight * (Deriv[1][2] * C1 - Deriv[2][1] * C2);
209 
210  STATE *ptr3 = &ek(in*3+0, jn*3+2);
211  /*ek(in*3+0,jn*3+2)*/ *ptr3++ += weight * (Deriv[2][0] * C1 - Deriv[0][2] * C2);
212 
213  /*ek(in*3+1,jn*3+2)*/ *ptr3++ += weight * (Deriv[2][1] * C1 - Deriv[1][2] * C2);
214 
215  /*ek(in*3+2,jn*3+2)*/ *ptr3 += weight * (( Deriv[0][0] + Deriv[1][1] ) * C1 + Deriv[2][2] * C3);
216 
217  }//jn
218  }//in
219 
220 #endif
221 #ifdef CODE2
222 
223  for(int jn = 0; jn < phr; jn++)
224  {
225  REAL dphij[3];
226  int kd;
227  for(kd = 0; kd < 3; kd++)
228  {
229  dphij[kd] = dphi(kd,jn);
230  ef(jn*3+kd, 0) += weight * ( locForce[kd] * phi(jn,0) - fPreStress[kd] * dphi(kd,jn) );
231  }//kd
232  for( int in = 0; in < phr; in++ )
233  {
234  REAL Deriv[3][3];
235  //Compute Deriv matrix
236  for(int ud = 0; ud < 3; ud++)
237  {
238  for(int vd = 0; vd < 3; vd++)
239  {
240  Deriv[vd][ud] = dphi(vd,in)*dphij[ud];
241  }//ud
242  }//vd
243 
244  //First equation Dot[Sigma1, gradV1]
245  STATE *ptr1 = &ek(in*3,jn*3);
246  STATE *ptr2 = &ek(in*3+0, jn*3+1);
247  STATE *ptr3 = &ek(in*3+0, jn*3+2);
248  /*ek(in*3+0,jn*3+0)*/ ptr1[0] += weight * (( Deriv[1][1] + Deriv[2][2] ) * C1 + Deriv[0][0] * C3);
249 
250  //Second equation Dot[Sigma2, gradV2]
251  /*ek(in*3+1,jn*3+0)*/ ptr1[1] += weight * (Deriv[0][1] * C1 - Deriv[1][0] * C2);
252 
253  //Third equation Dot[Sigma3, gradV3]
254  /*ek(in*3+2,jn*3+0)*/ ptr1[2] += weight * (Deriv[0][2] * C1 - Deriv[2][0] * C2);
255 
256  /*ek(in*3+0,jn*3+1)*/ ptr2[0] += weight * (Deriv[1][0] * C1 - Deriv[0][1] * C2);
257 
258  /*ek(in*3+1,jn*3+1)*/ ptr2[1] += weight * (( Deriv[0][0] + Deriv[2][2] ) * C1 + Deriv[1][1] * C3);
259 
260  /*ek(in*3+2,jn*3+1)*/ ptr2[2] += weight * (Deriv[1][2] * C1 - Deriv[2][1] * C2);
261 
262  /*ek(in*3+0,jn*3+2)*/ ptr3[0] += weight * (Deriv[2][0] * C1 - Deriv[0][2] * C2);
263 
264  /*ek(in*3+1,jn*3+2)*/ ptr3[1] += weight * (Deriv[2][1] * C1 - Deriv[1][2] * C2);
265 
266  /*ek(in*3+2,jn*3+2)*/ ptr3[2] += weight * (( Deriv[0][0] + Deriv[1][1] ) * C1 + Deriv[2][2] * C3);
267 
268  }//jn
269  }//in
270 
271 #endif
272 #ifdef CODE3
273 
274  for(int jn = 0; jn < phr; jn++)
275  {
276  REAL dphij[3];
277  int kd;
278  for(kd = 0; kd < 3; kd++)
279  {
280  dphij[kd] = dphi(kd,jn);
281  ef(jn*3+kd, 0) += weight * ( locForce[kd] * phi(jn,0) - fPreStress[kd] * dphi(kd,jn) );
282  }//kd
283 
284  const int stride = 3;
285  int phmax = (phr/stride)*stride;
286  int in=0;
287  for(in = 0; in < phmax; in+=stride )
288  {
289  REAL Deriv[3*stride][3];
290  //Compute Deriv matrix
291  for(int ud = 0; ud < 3; ud++)
292  {
293  for (int istr=0; istr<stride; istr++)
294  {
295  for(int vd = 0; vd < 3; vd++)
296  {
297  Deriv[vd+istr*stride][ud] = dphi(vd,in+istr)*dphij[ud];
298  }//ud
299  }
300  }//vd
301 
302  //First equation Dot[Sigma1, gradV1]
303  STATE *ptr1 = &ek(in*3,jn*3);
304  STATE *ptr2 = &ek(in*3+0, jn*3+1);
305  STATE *ptr3 = &ek(in*3+0, jn*3+2);
306  ptr1[0] += weight * (( Deriv[1][1] + Deriv[2][2] ) * C1 + Deriv[0][0] * C3);
307  ptr1[1] += weight * (Deriv[0][1] * C1 - Deriv[1][0] * C2);
308  ptr1[2] += weight * (Deriv[0][2] * C1 - Deriv[2][0] * C2);
309  ptr2[0] += weight * (Deriv[1][0] * C1 - Deriv[0][1] * C2);
310  ptr2[1] += weight * (( Deriv[0][0] + Deriv[2][2] ) * C1 + Deriv[1][1] * C3);
311  ptr2[2] += weight * (Deriv[1][2] * C1 - Deriv[2][1] * C2);
312  ptr3[0] += weight * (Deriv[2][0] * C1 - Deriv[0][2] * C2);
313  ptr3[1] += weight * (Deriv[2][1] * C1 - Deriv[1][2] * C2);
314  ptr3[2] += weight * (( Deriv[0][0] + Deriv[1][1] ) * C1 + Deriv[2][2] * C3);
315  const int add1 = 3;
316  ptr1[0+add1] += weight * (( Deriv[1+add1][1] + Deriv[2+add1][2] ) * C1 + Deriv[0+add1][0] * C3);
317  ptr1[1+add1] += weight * (Deriv[0+add1][1] * C1 - Deriv[1+add1][0] * C2);
318  ptr1[2+add1] += weight * (Deriv[0+add1][2] * C1 - Deriv[2+add1][0] * C2);
319  ptr2[0+add1] += weight * (Deriv[1+add1][0] * C1 - Deriv[0+add1][1] * C2);
320  ptr2[1+add1] += weight * (( Deriv[0+add1][0] + Deriv[2+add1][2] ) * C1 + Deriv[1+add1][1] * C3);
321  ptr2[2+add1] += weight * (Deriv[1+add1][2] * C1 - Deriv[2+add1][1] * C2);
322  ptr3[0+add1] += weight * (Deriv[2+add1][0] * C1 - Deriv[0+add1][2] * C2);
323  ptr3[1+add1] += weight * (Deriv[2+add1][1] * C1 - Deriv[1+add1][2] * C2);
324  ptr3[2+add1] += weight * (( Deriv[0+add1][0] + Deriv[1+add1][1] ) * C1 + Deriv[2+add1][2] * C3);
325  const int add2 = 6;
326  ptr1[0+add2] += weight * (( Deriv[1+add2][1] + Deriv[2+add2][2] ) * C1 + Deriv[0+add2][0] * C3);
327  ptr1[1+add2] += weight * (Deriv[0+add2][1] * C1 - Deriv[1+add2][0] * C2);
328  ptr1[2+add2] += weight * (Deriv[0+add2][2] * C1 - Deriv[2+add2][0] * C2);
329  ptr2[0+add2] += weight * (Deriv[1+add2][0] * C1 - Deriv[0+add2][1] * C2);
330  ptr2[1+add2] += weight * (( Deriv[0+add2][0] + Deriv[2+add2][2] ) * C1 + Deriv[1+add2][1] * C3);
331  ptr2[2+add2] += weight * (Deriv[1+add2][2] * C1 - Deriv[2+add2][1] * C2);
332  ptr3[0+add2] += weight * (Deriv[2+add2][0] * C1 - Deriv[0+add2][2] * C2);
333  ptr3[1+add2] += weight * (Deriv[2+add2][1] * C1 - Deriv[1+add2][2] * C2);
334  ptr3[2+add2] += weight * (( Deriv[0+add2][0] + Deriv[1+add2][1] ) * C1 + Deriv[2+add2][2] * C3);
335  }//jn
336  for(; in < phr; in++ )
337  {
338  REAL Deriv[3][3];
339  //Compute Deriv matrix
340  for(int ud = 0; ud < 3; ud++)
341  {
342  for(int vd = 0; vd < 3; vd++)
343  {
344  Deriv[vd][ud] = dphi(vd,in)*dphij[ud];
345  }//ud
346  }//vd
347 
348  //First equation Dot[Sigma1, gradV1]
349  STATE *ptr1 = &ek(in*3,jn*3);
350  STATE *ptr2 = &ek(in*3+0, jn*3+1);
351  STATE *ptr3 = &ek(in*3+0, jn*3+2);
352  /*ek(in*3+0,jn*3+0)*/ ptr1[0] += weight * (( Deriv[1][1] + Deriv[2][2] ) * C1 + Deriv[0][0] * C3);
353 
354  //Second equation Dot[Sigma2, gradV2]
355  /*ek(in*3+1,jn*3+0)*/ ptr1[1] += weight * (Deriv[0][1] * C1 - Deriv[1][0] * C2);
356 
357  //Third equation Dot[Sigma3, gradV3]
358  /*ek(in*3+2,jn*3+0)*/ ptr1[2] += weight * (Deriv[0][2] * C1 - Deriv[2][0] * C2);
359 
360  /*ek(in*3+0,jn*3+1)*/ ptr2[0] += weight * (Deriv[1][0] * C1 - Deriv[0][1] * C2);
361 
362  /*ek(in*3+1,jn*3+1)*/ ptr2[1] += weight * (( Deriv[0][0] + Deriv[2][2] ) * C1 + Deriv[1][1] * C3);
363 
364  /*ek(in*3+2,jn*3+1)*/ ptr2[2] += weight * (Deriv[1][2] * C1 - Deriv[2][1] * C2);
365 
366  /*ek(in*3+0,jn*3+2)*/ ptr3[0] += weight * (Deriv[2][0] * C1 - Deriv[0][2] * C2);
367 
368  /*ek(in*3+1,jn*3+2)*/ ptr3[1] += weight * (Deriv[2][1] * C1 - Deriv[1][2] * C2);
369 
370  /*ek(in*3+2,jn*3+2)*/ ptr3[2] += weight * (( Deriv[0][0] + Deriv[1][1] ) * C1 + Deriv[2][2] * C3);
371 
372  }//jn
373  }//in
374 
375 #endif
376 
377 #ifdef CODE4
378  static TPZFNMatrix<300,REAL> BMatrix(6,ek.Rows(),0.),DBMatrix(6,ek.Rows(),0.);
379  int64_t nphi = data.phi.Rows();
380  for (int64_t iph=0; iph<nphi; iph++) {
381  BMatrix(0,3*iph) = data.dphix(0,iph);
382  BMatrix(1,3*iph+1) = data.dphix(1,iph);
383  BMatrix(2,3*iph+2) = data.dphix(2,iph);
384  BMatrix(3,3*iph) = data.dphix(1,iph);
385  BMatrix(3,3*iph+1) = data.dphix(0,iph);
386  BMatrix(4,3*iph) = data.dphix(2,iph);
387  BMatrix(4,3*iph+2) = data.dphix(0,iph);
388  BMatrix(5,3*iph+1) = data.dphix(2,iph);
389  BMatrix(5,3*iph+2) = data.dphix(1,iph);
390  DBMatrix(0,3*iph) = data.dphix(0,iph)*(1.-fPoisson);
391  DBMatrix(0,3*iph+1) = data.dphix(1,iph)*fPoisson;
392  DBMatrix(0,3*iph+2) = data.dphix(2,iph)*fPoisson;
393  DBMatrix(1,3*iph) = data.dphix(0,iph)*fPoisson;
394  DBMatrix(1,3*iph+1) = data.dphix(1,iph)*(1.-fPoisson);
395  DBMatrix(1,3*iph+2) = data.dphix(2,iph)*fPoisson;
396  DBMatrix(2,3*iph) = data.dphix(0,iph)*fPoisson;
397  DBMatrix(2,3*iph+1) = data.dphix(1,iph)*fPoisson;
398  DBMatrix(2,3*iph+2) = data.dphix(2,iph)*(1.-fPoisson);
399  DBMatrix(3,3*iph) = data.dphix(1,iph)*(1.-2.*fPoisson)/2.;
400  DBMatrix(3,3*iph+1) = data.dphix(0,iph)*(1.-2.*fPoisson)/2.;
401  DBMatrix(4,3*iph) = data.dphix(2,iph)*(1.-2.*fPoisson)/2.;
402  DBMatrix(4,3*iph+2) = data.dphix(0,iph)*(1.-2.*fPoisson)/2.;
403  DBMatrix(5,3*iph+1) = data.dphix(2,iph)*(1.-2.*fPoisson)/2.;
404  DBMatrix(5,3*iph+2) = data.dphix(1,iph)*(1.-2.*fPoisson)/2.;
405  const REAL mult= fE/((1+fPoisson)*(1.-2.*fPoisson));
406  }
407 #endif
408 #if !defined CODE0 && !defined CODE1 && !defined CODE2 && !defined CODE3
409  DebugStop();
410 #endif
411 #ifdef PZDEBUG
412  if ( !ek.VerifySymmetry( 1.e-8 ) ) PZError << __PRETTY_FUNCTION__ << "\nERROR - NON SYMMETRIC MATRIX" << std::endl;
413 #endif
414 }//method
415 
416 
417 void TPZElasticity3D::ContributeVecShape(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
418 {
419  TPZFMatrix<REAL> & dphi = data.dphix;
420  TPZFMatrix<REAL> & phi = data.phi;
421 
422  int phc = phi.Cols();
423  int efc = ef.Cols();
424 
425  TPZManVector<STATE,3> locForce(fForce);
426  if(fForcingFunction)
427  {
428  TPZManVector<STATE> res(3);
429  fForcingFunction->Execute(data.x,res);
430  locForce[0] = res[0];
431  locForce[1] = res[1];
432  locForce[2] = res[2];
433  }
434 
435  REAL dvxdx, dvxdy, dvxdz;
436  REAL dvydx, dvydy, dvydz;
437  REAL dvzdx, dvzdy, dvzdz;
438 
439  REAL duxdx, duxdy, duxdz;
440  REAL duydx, duydy, duydz;
441  REAL duzdx, duzdy, duzdz;
442 
443  /*
444  * Plain strain materials values
445  */
446  REAL lambda = fE*fPoisson/((1.+fPoisson)*(1.-2.*fPoisson));
447  REAL mu = fE/(2.*(1.+fPoisson));
448 
449  for( int in = 0; in < phc; in++ )
450  {
451  //x
452  dvxdx = dphi(0,in);
453  dvxdy = dphi(1,in);
454  dvxdz = dphi(2,in);
455 
456  //y
457  dvydx = dphi(3,in);
458  dvydy = dphi(4,in);
459  dvydz = dphi(5,in);
460 
461  //z
462  dvzdx = dphi(6,in);
463  dvzdy = dphi(7,in);
464  dvzdz = dphi(8,in);
465 
466  for (int col = 0; col < efc; col++)
467  {
468  ef(in,col) += weight*( locForce[0] * phi(0, in)
469  + locForce[1] * phi(1, in)
470  + locForce[2] * phi(2, in)
471  - dvxdx * fPreStress[0]
472  - dvydy * fPreStress[1]
473  - dvzdz * fPreStress[2]);
474  }
475  for( int jn = 0; jn < phc; jn++ )
476  {
477  //x
478  duxdx = dphi(0,jn);
479  duxdy = dphi(1,jn);
480  duxdz = dphi(2,jn);
481 
482  //y
483  duydx = dphi(3,jn);
484  duydy = dphi(4,jn);
485  duydz = dphi(5,jn);
486 
487  //z
488  duzdx = dphi(6,jn);
489  duzdy = dphi(7,jn);
490  duzdz = dphi(8,jn);
491 
492  REAL eq1 = duydy*dvxdx*lambda + duzdz*dvxdx*lambda + duxdy*dvydx*mu +
493  duydx*dvydx*mu + duxdz*dvzdx*mu + duzdx*dvzdx*mu +
494  duxdx*dvxdx*(lambda + 2.*mu);
495 
496  REAL eq2 = duxdx*dvydy*lambda + duzdz*dvydy*lambda + duxdy*dvxdy*mu +
497  duydx*dvxdy*mu + duydz*dvzdy*mu + duzdy*dvzdy*mu +
498  duydy*dvydy*(lambda + 2.*mu);
499 
500  REAL eq3 = duxdx*dvzdz*lambda + duydy*dvzdz*lambda + duxdz*dvxdz*mu +
501  duzdx*dvxdz*mu + duydz*dvydz*mu + duzdy*dvydz*mu +
502  duzdz*dvzdz*(lambda + 2.*mu);
503 
504  ek(in,jn) += weight * (eq1 + eq2 + eq3);
505  }
506  }
507 }
508 
509 void TPZElasticity3D::ContributeVecShapeBC(TPZMaterialData & data, REAL weight,
510  TPZFMatrix<STATE> & ek, TPZFMatrix<STATE> & ef,TPZBndCond &bc)
511 {
512  TPZFMatrix<REAL> & phi = data.phi;
513 
514  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
515 
516  int phc = phi.Cols();
517  short in,jn;
518 
519  switch (bc.Type())
520  {
521  case 0:// Dirichlet condition
522  {
523  for(in = 0 ; in < phc; in++)
524  {
525  for (int il = 0; il < fNumLoadCases; il++)
526  {
527  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
528  ef(in,il) += weight * BIGNUMBER * ( v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in) + v2(2,il)*phi(2,in) );
529  }
530  for (jn = 0 ; jn < phc; jn++)
531  {
532  ek(in,jn) += weight * BIGNUMBER * ( phi(0,in)*phi(0,jn) + phi(1,in)*phi(1,jn) + phi(2,in)*phi(2,jn) );
533  }
534  }
535  break;
536  }
537  case 1:// Neumann condition
538  {
539  for (in = 0; in < phc; in++)
540  {
541  for (int il = 0; il <fNumLoadCases; il++)
542  {
543  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
544  ef(in,il) += weight * ( v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in) + v2(2,il)*phi(2,in) );
545  }
546  }
547  break;
548  }
549  case 2:// condicao mista
550  {
551  for(in = 0 ; in < phc; in++)
552  {
553  for (int il = 0; il <fNumLoadCases; il++)
554  {
555  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
556  ef(in,il)+= weight * ( v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in) + v2(2,il)*phi(2,in) );
557  }
558 
559  for (jn = 0; jn <phc; jn++)
560  {
561 
562  ek(in,jn) += bc.Val1()(0,0)*phi(0,in)*phi(0,jn)*weight
563 
564  + bc.Val1()(1,0)*phi(1,in)*phi(0,jn)*weight
565 
566  + bc.Val1()(2,0)*phi(2,in)*phi(0,jn)*weight
567 
568 
569  + bc.Val1()(0,1)*phi(0,in)*phi(1,jn)*weight
570 
571  + bc.Val1()(1,1)*phi(1,in)*phi(1,jn)*weight
572 
573  + bc.Val1()(2,1)*phi(2,in)*phi(1,jn)*weight
574 
575 
576  + bc.Val1()(0,2)*phi(0,in)*phi(2,jn)*weight
577 
578  + bc.Val1()(1,2)*phi(1,in)*phi(2,jn)*weight
579 
580  + bc.Val1()(2,2)*phi(2,in)*phi(2,jn)*weight;
581  }
582  break;
583  }
584  }
585  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
586  {
587  for(in = 0 ; in < phc; in++)
588  {
589  for (int il = 0; il < fNumLoadCases; il++)
590  {
591  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
592  for (jn = 0 ; jn < phc; jn++)
593  {
594  ek(in,jn) += weight * BIGNUMBER * ( v2(0,il)*phi(0,in)*phi(0,jn) + v2(1,il)*phi(1,in)*phi(1,jn) + v2(2,il)*phi(2,in)*phi(2,jn) );
595  }
596  }
597  }
598  break;
599  }
600  case 4: // stressField Neumann condition
601  {
602  DebugStop();//Nao implementado!!!
603  break;
604  }
605  default:
606  PZError << "TPZElastitity3D::ContributeBC error - Wrong boundary condition type" << std::endl;
607  }
608 }
609 
610 
611 
612 void TPZElasticity3D::ContributeBC(TPZMaterialData &data,
613  REAL weight,
614  TPZFMatrix<STATE> &ek,
615  TPZFMatrix<STATE> &ef,
616  TPZBndCond &bc){
617 
618  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
619  if(shapetype==data.EVecShape){
620  ContributeVecShapeBC(data,weight,ek, ef,bc);
621  return;
622  }
623 
624  TPZFMatrix<REAL> &phi = data.phi;
625 
626  const STATE BIGNUMBER = 1.e12;
627 
628  const int phr = phi.Rows();
629  int in,jn,idf,jdf;
630  TPZFNMatrix<9,STATE> val1loc(bc.Val1()),val2loc(bc.Val2());
631 
632  if (bc.HasForcingFunction()) {
633  TPZManVector<STATE,3> val2vec(3);
634 
635  bc.ForcingFunction()->Execute(data.x, val2vec, val1loc);
636  val2loc(0,0) = val2vec[0];
637  val2loc(1,0) = val2vec[1];
638  val2loc(2,0) = val2vec[2];
639  // we assume type 2 is displacement value is weakly imposed
640  if(bc.Type() == 2)
641  {
642  val1loc = bc.Val1();
643  for (int i=0; i<3; i++) {
644  val2loc(i,0) = 0.;
645  for (int j=0; j<3; j++) {
646  val2loc(i,0) += val1loc(i,j)*val2vec[j];
647  }
648  }
649  }
650  if(bc.Type() == 1)
651  {
652  for (int i=0; i<3; i++) {
653  val2loc(i,0) = 0.;
654  for (int j=0; j<3; j++) {
655  val2loc(i,0) += val1loc(i,j)*data.normal[j];
656  }
657  }
658  }
659  }
660 
661 
662  int numbersol = data.sol.size();
663  if (numbersol != 1) {
664  DebugStop();
665  }
666 
667  switch (bc.Type()) {
668  case 0: // Dirichlet condition
669  for(in = 0 ; in < phr; in++) {
670  ef(3*in+0,0) += BIGNUMBER * val2loc(0,0) * phi(in,0) * weight;
671  ef(3*in+1,0) += BIGNUMBER * val2loc(1,0) * phi(in,0) * weight;
672  ef(3*in+2,0) += BIGNUMBER * val2loc(2,0) * phi(in,0) * weight;
673 
674  for (jn = 0 ; jn < phr; jn++) {
675  ek(3*in+0,3*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
676  ek(3*in+1,3*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
677  ek(3*in+2,3*jn+2) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
678  }//jn
679  }//in
680  break;
681 
682  case 1: // Neumann condition
683  for(in = 0 ; in < phi.Rows(); in++) {
684  ef(3*in+0,0) += val2loc(0,0) * phi(in,0) * weight;
685  ef(3*in+1,0) += val2loc(1,0) * phi(in,0) * weight;
686  ef(3*in+2,0) += val2loc(2,0) * phi(in,0) * weight;
687  }//in
688  break;
689  case 2: // Mixed condition
690  for(in = 0 ; in < phi.Rows(); in++) {
691  ef(3*in+0,0) += val2loc(0,0) * phi(in,0) * weight;
692  ef(3*in+1,0) += val2loc(1,0) * phi(in,0) * weight;
693  ef(3*in+2,0) += val2loc(2,0) * phi(in,0) * weight;
694  for(jn=0; jn<phi.Rows(); jn++)
695  {
696  for(idf=0; idf<3; idf++) for(jdf=0; jdf<3; jdf++)
697  {
698  ek(3*in+idf,3*jn+jdf) += val1loc(idf,jdf)*weight*phi(in,0)*phi(jn,0);
699  }
700  }
701  }//in
702  break;
703  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
704  for(in = 0 ; in < phr; in++) {
705  for (jn = 0 ; jn < phr; jn++) {
706  ek(3*in+0,3*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * val2loc(0,0);
707  ek(3*in+1,3*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * val2loc(1,0);
708  ek(3*in+2,3*jn+2) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * val2loc(2,0);
709  }//jn
710  }//in
711  break;
712 
713  case 4: // stressField Neumann condition
714  for(in = 0; in < 3; in ++)
715  val2loc(in,0) = - ( val1loc(in,0) * data.normal[0] +
716  val1loc(in,1) * data.normal[1] +
717  val1loc(in,2) * data.normal[2] );
718  // The normal vector points towards the neighbour. The negative sign is there to
719  // reflect the outward normal vector.
720  for(in = 0 ; in < phi.Rows(); in++) {
721  ef(3*in+0,0) += val2loc(0,0) * phi(in,0) * weight;
722  ef(3*in+1,0) += val2loc(1,0) * phi(in,0) * weight;
723  ef(3*in+2,0) += val2loc(2,0) * phi(in,0) * weight;
724  //cout << "normal:" << data.normal[0] << ' ' << data.normal[1] << ' ' << data.normal[2] << endl;
725  //cout << "val2: " << val2loc(0,0) << ' ' << val2loc(1,0) << ' ' << val2loc(2,0) << endl;
726  }
727  break;
728  case 5: // Directional Dirichlet condition on x
729  for(in = 0 ; in < phr; in++) {
730  ef(3*in+0,0) += BIGNUMBER * val2loc(0,0) * phi(in,0) * weight;
731 
732  for (jn = 0 ; jn < phr; jn++) {
733  ek(3*in+0,3*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
734  }//jn
735  }//in
736  break;
737  case 6: // Directional Dirichlet condition on y
738  for(in = 0 ; in < phr; in++) {
739  ef(3*in+1,0) += BIGNUMBER * val2loc(1,0) * phi(in,0) * weight;
740 
741  for (jn = 0 ; jn < phr; jn++) {
742  ek(3*in+1,3*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
743  }//jn
744  }//in
745  break;
746  case 7: // Directional Dirichlet condition on z
747  for(in = 0 ; in < phr; in++) {
748  ef(3*in+2,0) += BIGNUMBER * val2loc(2,0) * phi(in,0) * weight;
749 
750  for (jn = 0 ; jn < phr; jn++) {
751  ek(3*in+2,3*jn+2) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
752  }//jn
753  }//in
754  break;
755  case 8: // Directional Dirichlet condition on xz
756  for(in = 0 ; in < phr; in++) {
757  ef(3*in+0,0) += BIGNUMBER * val2loc(0,0) * phi(in,0) * weight;
758  ef(3*in+2,0) += BIGNUMBER * val2loc(2,0) * phi(in,0) * weight;
759 
760  for (jn = 0 ; jn < phr; jn++) {
761  ek(3*in+0,3*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
762  ek(3*in+2,3*jn+2) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
763  }//jn
764  }//in
765  break;
766  default:
767  PZError << "TPZElastitity3D::ContributeBC error - Wrong boundary condition type" << std::endl;
768  }//switch
769 }//method
770 
771 int TPZElasticity3D::VariableIndex(const std::string &name) {
772  if(!strcmp("Displacement",name.c_str())) return TPZElasticity3D::EDisplacement;
773  if(!strcmp("state",name.c_str())) return TPZElasticity3D::EDisplacement;
774  if(!strcmp("DisplacementX",name.c_str())) return TPZElasticity3D::EDisplacementX;
775  if(!strcmp("DisplacementY",name.c_str())) return TPZElasticity3D::EDisplacementY;
776  if(!strcmp("DisplacementZ",name.c_str())) return TPZElasticity3D::EDisplacementZ;
777  if(!strcmp("PrincipalStress", name.c_str())) return TPZElasticity3D::EPrincipalStress;
778  if(!strcmp("PrincipalStrain", name.c_str())) return TPZElasticity3D::EPrincipalStrain;
779  if(!strcmp("VonMises", name.c_str())) return TPZElasticity3D::EVonMisesStress;
780  if(!strcmp("Stress", name.c_str())) return TPZElasticity3D::EStress;
781  if(!strcmp("Strain", name.c_str())) return TPZElasticity3D::EStrain;
782  if(!strcmp("Stress1", name.c_str())) return TPZElasticity3D::EStress1;
783  if(!strcmp("Strain1", name.c_str())) return TPZElasticity3D::EStrain1;
784  if(!strcmp("NormalStress",name.c_str())) return TPZElasticity3D::ENormalStress;
785  if(!strcmp("NormalStrain",name.c_str())) return TPZElasticity3D::ENormalStrain;
786  if(!strcmp("StressX",name.c_str())) return TPZElasticity3D::EStressX;
787  if(!strcmp("StressY",name.c_str())) return TPZElasticity3D::EStressY;
788  if(!strcmp("StressZ",name.c_str())) return TPZElasticity3D::EStressZ;
789  if(!strcmp("I1",name.c_str())) return TPZElasticity3D::EI1;
790  if(!strcmp("I2",name.c_str())) return TPZElasticity3D::EI2;
791  if(!strcmp("I3",name.c_str())) return TPZElasticity3D::EI3;
792  if(!strcmp("PlasticFunction",name.c_str())) return TPZElasticity3D::EPlasticFunction;
793 
794  // cout << "TPZElasticityMaterial::VariableIndex Error\n";
795  return TPZMaterial::VariableIndex(name);
796 }
797 
798 int TPZElasticity3D::NSolutionVariables(int var) {
799  switch(var) {
800  case TPZElasticity3D::EDisplacement :
801  case TPZElasticity3D::EPrincipalStress :
802  case TPZElasticity3D::EPrincipalStrain :
803  case TPZElasticity3D::EPrincipalDirection1 :
804  case TPZElasticity3D::EPrincipalDirection2 :
805  case TPZElasticity3D::EPrincipalDirection3 :
806  case TPZElasticity3D::EStress :
807  case TPZElasticity3D::EStrain :
808  case TPZElasticity3D::ENormalStress :
809  case TPZElasticity3D::ENormalStrain :
810  return 3;
811  case TPZElasticity3D::EDisplacementX :
812  case TPZElasticity3D::EDisplacementY :
813  case TPZElasticity3D::EDisplacementZ :
814  case TPZElasticity3D::EVonMisesStress :
815  case TPZElasticity3D::EStrain1 :
816  case TPZElasticity3D::EStress1 :
817  case TPZElasticity3D::EStressX :
818  case TPZElasticity3D::EStressY :
819  case TPZElasticity3D::EStressZ :
820  case TPZElasticity3D::EI1 :
821  case TPZElasticity3D::EI2 :
822  case TPZElasticity3D::EI3 :
823  case TPZElasticity3D::EPlasticFunction :
824  return 1;
825  default:
826  return TPZMaterial::NSolutionVariables(var);
827  }
828  return -1;
829 }
830 
831 void TPZElasticity3D::Solution(TPZVec<STATE> &Sol,TPZFMatrix<STATE> &DSol,TPZFMatrix<REAL> &axes,int var,TPZVec<STATE> &Solout) {
832 
833  TPZFNMatrix<9,STATE> DSolXY(3,3);
834  TPZAxesTools<STATE>::Axes2XYZ(DSol, DSolXY, axes);
835  if(var == TPZElasticity3D::EDisplacement) {
836  int i;
837  for(i = 0; i < 3; i++){
838  Solout[i] = Sol[i];
839  }//for
840  return;
841  }//TPZElasticity3D::EDisplacement
842 
843  if(var == TPZElasticity3D::EDisplacementX) {
844  // int i;
845  Solout[0] = Sol[0];
846  return;
847  }//TPZElasticity3D::EDisplacementX
848 
849  if(var == TPZElasticity3D::EDisplacementY) {
850  // int i;
851  Solout[0] = Sol[1];
852  return;
853  }//TPZElasticity3D::EDisplacementY
854 
855  if(var == TPZElasticity3D::EDisplacementZ) {
856  // int i;
857  Solout[0] = Sol[2];
858  return;
859  }//TPZElasticity3D::EDisplacementZ
860 
861  if(var == TPZElasticity3D::EPlasticFunction){
862  TPZFNMatrix<9,STATE> StressTensor(3,3);
863  this->ComputeStressTensor(StressTensor, DSolXY);
864  if(fPlasticPostProc == EVonMises) Solout[0] = this->VonMisesPlasticFunction(StressTensor);
865  else if(fPlasticPostProc == EMohrCoulomb) Solout[0] = this->MohrCoulombPlasticFunction(StressTensor);
866  else DebugStop();
867  return;
868  }
869 
870  if(var == TPZElasticity3D::EPrincipalStress) {
871  TPZFNMatrix<9,STATE> StressTensor(3,3);
872  this->ComputeStressTensor(StressTensor, DSolXY);
873  int64_t numiterations = 1000;
874  REAL tol = TPZElasticity3D::gTolerance;
875  TPZManVector<STATE,3> eigv;
876  bool result;
877  result = StressTensor.SolveEigenvaluesJacobi(numiterations, tol, &eigv);
878  for (int i=0; i<eigv.size(); i++) {
879  Solout[i] = eigv[i];
880  }
881 #ifdef PZDEBUG
882  if (result == false){
883  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
884  }
885 #endif
886  }//TPZElasticity3D::EPrincipalStress
887 
888  if(var == TPZElasticity3D::EStress1){
889  TPZFNMatrix<9,STATE> StressTensor(3,3);
890  TPZManVector<STATE, 3> PrincipalStress(3);
891  this->ComputeStressTensor(StressTensor, DSolXY);
892  int64_t numiterations = 1000;
893  REAL tol = TPZElasticity3D::gTolerance;
894  bool result;
895  result = StressTensor.SolveEigenvaluesJacobi(numiterations, tol, &PrincipalStress);
896  Solout[0] = PrincipalStress[0];
897 #ifdef PZDEBUG
898  if (result == false){
899  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
900  }
901 #endif
902  }//TPZElasticity3D::EStress1
903 
904 
905  if(var == TPZElasticity3D::EPrincipalStrain){
906  TPZFNMatrix<9,STATE> StrainTensor(3,3);
907  this->ComputeStrainTensor(StrainTensor, DSolXY);
908  int64_t numiterations = 1000;
909  REAL tol = TPZElasticity3D::gTolerance;
910  TPZManVector<STATE,3> eigv;
911  bool result;
912  result = StrainTensor.SolveEigenvaluesJacobi(numiterations, tol, &eigv);
913  for (int i=0; i<eigv.size(); i++) {
914  Solout[i] = eigv[i];
915  }
916 #ifdef PZDEBUG
917  if (result == false){
918  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
919  }
920 #endif
921  }//TPZElasticity3D::EPrincipalStrain
922 
923  if(var == TPZElasticity3D::EStrain1){
924  TPZFNMatrix<9,STATE> StrainTensor(3,3);
925  TPZManVector<STATE, 3> PrincipalStrain(3);
926  this->ComputeStrainTensor(StrainTensor, DSolXY);
927  int64_t numiterations = 1000;
928  REAL tol = TPZElasticity3D::gTolerance;
929  bool result;
930  result = StrainTensor.SolveEigenvaluesJacobi(numiterations, tol, &PrincipalStrain);
931  Solout[0] = PrincipalStrain[0];
932 #ifdef PZDEBUG
933  if (result == false){
934  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
935  }
936 #endif
937  }//TPZElasticity3D::EPrincipalStrain
938 
939  TPZManVector<STATE,3> eigvec;
940  if(var == TPZElasticity3D::EPrincipalDirection1){
941  this->PrincipalDirection(DSolXY, eigvec, 0);
942  for (int i=0; i<eigvec.size(); i++) {
943  Solout[i] = eigvec[i];
944  }
945  }//TPZElasticity3D::EPrincipalDirection1
946 
947  if(var == TPZElasticity3D::EPrincipalDirection2){
948  this->PrincipalDirection(DSolXY, eigvec, 1);
949  for (int i=0; i<eigvec.size(); i++) {
950  Solout[i] = eigvec[i];
951  }
952  }//TPZElasticity3D::EPrincipalDirection2
953 
954  if(var == TPZElasticity3D::EPrincipalDirection3){
955  this->PrincipalDirection(DSolXY, eigvec, 2);
956  for (int i=0; i<eigvec.size(); i++) {
957  Solout[i] = eigvec[i];
958  }
959  }//TPZElasticity3D::EPrincipalDirection3
960 
961 
962  if(var == TPZElasticity3D::EVonMisesStress){
963  TPZManVector<STATE,3> PrincipalStress(3);
964  TPZFNMatrix<9,STATE> StressTensor(3,3);
965  this->ComputeStressTensor(StressTensor, DSolXY);
966  int64_t numiterations = 1000;
967  REAL tol = TPZElasticity3D::gTolerance;
968  bool result;
969  result = StressTensor.SolveEigenvaluesJacobi(numiterations, tol, &PrincipalStress);
970 #ifdef PZDEBUG
971  if (result == false){
972  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
973  }
974 #endif
975 
976  Solout.Resize(1);
977  Solout[0] = ( PrincipalStress[0] - PrincipalStress[1] ) * ( PrincipalStress[0] - PrincipalStress[1] )
978  + ( PrincipalStress[1] - PrincipalStress[2] ) * ( PrincipalStress[1] - PrincipalStress[2] )
979  + ( PrincipalStress[2] - PrincipalStress[0] ) * ( PrincipalStress[2] - PrincipalStress[0] );
980  Solout[0] = Solout[0] / (2. * this->fFy * this->fFy);
981  }//TPZElasticity3D::EVonMisesStress
982 
983  if(var == TPZElasticity3D::EStress){
984  TPZFNMatrix<6,STATE> Stress(6,1);
985  this->ComputeStressVector(Stress, DSolXY);
986  TPZManVector<STATE,3> eigvec;
987  this->ApplyDirection(Stress, eigvec);
988  for (int i=0; i<eigvec.size(); i++) {
989  Solout[i] = eigvec[i];
990  }
991  return;
992  }//TPZElasticity3D::EStress
993 
994  if(var == TPZElasticity3D::EStrain){
995  TPZFNMatrix<6,STATE> Strain(6,1);
996  this->ComputeStrainVector(Strain, DSolXY);
997  TPZManVector<STATE,3> eigvec;
998  this->ApplyDirection(Strain, eigvec);
999  for (int i=0; i<eigvec.size(); i++) {
1000  Solout[i] = eigvec[i];
1001  }
1002  return;
1003  }//TPZElasticity3D::EStrain
1004 
1005  if(var == TPZElasticity3D::ENormalStress){
1006  TPZFNMatrix<6,STATE> Stress(6,1);
1007  this->ComputeStressVector(Stress, DSolXY);
1008  Solout[0] = Stress(0,0);
1009  Solout[1] = Stress(1,0);
1010  Solout[2] = Stress(2,0);
1011  return;
1012  }//TPZElasticity3D::ENormalStress
1013 
1014  if(var == TPZElasticity3D::ENormalStrain){
1015  TPZFNMatrix<6,STATE> Strain(6,1);
1016  this->ComputeStrainVector(Strain, DSolXY);
1017  Solout[0] = Strain(0,0);
1018  Solout[1] = Strain(1,0);
1019  Solout[2] = Strain(2,0);
1020  return;
1021  }//TPZElasticity3D::ENormalStrain
1022 
1023  if(var == TPZElasticity3D::EStressX){
1024  TPZFNMatrix<9,STATE> Stress(3,3);
1025  this->ComputeStressTensor(Stress, DSolXY);
1026  Solout[0] = Stress(0,0);
1027  return;
1028  }//TPZElasticity3D::EStressX
1029 
1030  if(var == TPZElasticity3D::EStressY){
1031  TPZFNMatrix<9,STATE> Stress(3,3);
1032  this->ComputeStressTensor(Stress, DSolXY);
1033  Solout[0] = Stress(1,1);
1034  return;
1035  }//TPZElasticity3D::EStressY
1036 
1037  if(var == TPZElasticity3D::EStressZ){
1038  TPZFNMatrix<9,STATE> Stress(3,3);
1039  this->ComputeStressTensor(Stress, DSolXY);
1040  Solout[0] = Stress(2,2);
1041  return;
1042  }//TPZElasticity3D::EStressZ
1043 
1044  if(var == TPZElasticity3D::EI1){
1045  TPZFNMatrix<9,STATE> Stress(3,3);
1046  this->ComputeStressTensor(Stress, DSolXY);
1047  Solout[0] = Stress(0,0) + Stress(1,1) + Stress(2,2);//Calculado no Mathematica
1048 
1049  return;
1050  }//I1
1051 
1052  if(var == TPZElasticity3D::EI2){
1053  TPZFNMatrix<9,STATE> Stress(3,3);
1054  this->ComputeStressTensor(Stress, DSolXY);
1055 
1056  Solout[0] = -(Stress(0,1)*Stress(1,0)) - Stress(0,2)*Stress(2,0) - Stress(1,2)*Stress(2,1) +
1057  Stress(1,1)*Stress(2,2) + Stress(0,0)*(Stress(1,1) + Stress(2,2));//Calculado no Mathematica
1058 
1059  return;
1060  }//I2
1061 
1062  if(var == TPZElasticity3D::EI3){
1063  TPZFNMatrix<9,STATE> Stress(3,3);
1064  this->ComputeStressTensor(Stress, DSolXY);
1065  Solout[0] = -(Stress(0,2)*Stress(1,1)*Stress(2,0)) + Stress(0,1)*Stress(1,2)*Stress(2,0) +
1066  Stress(0,2)*Stress(1,0)*Stress(2,1) - Stress(0,0)*Stress(1,2)*Stress(2,1) -
1067  Stress(0,1)*Stress(1,0)*Stress(2,2) + Stress(0,0)*Stress(1,1)*Stress(2,2);//Calculado no Mathematica
1068 
1069  return;
1070  }//I3
1071  TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
1072 }//Solution
1073 
1074 void TPZElasticity3D::Errors(TPZVec<REAL> &x,TPZVec<STATE> &u, TPZFMatrix<STATE> &dudaxes,
1075  TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux,
1076  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values){
1077  int i, j;
1078  TPZFNMatrix<9,STATE> dudx(3,3);
1079  TPZAxesTools<STATE>::Axes2XYZ(dudaxes, dudx, axes);
1081  STATE L2 = 0.;
1082  for(i = 0; i < 3; i++) L2 += (u[i] - u_exact[i]) * (u[i] - u_exact[i]);
1083 
1085  REAL SemiH1 = 0.;
1086  for(i = 0; i < 3; i++) for(j = 0; j < 3; j++) SemiH1 += (dudx(i,j) - du_exact(i,j)) * (dudx(i,j) - du_exact(i,j));
1087 
1088  TPZFNMatrix<9,STATE> diffdsol(3,3,0.), diffstress(3,3,0.);
1089  for (int i=0; i<3; i++) {
1090  for (int j=0; j<3; j++) {
1091  diffdsol(i,j) = dudx(i,j)-du_exact(i,j);
1092  }
1093  }
1094  ComputeStressTensor(diffstress, diffdsol);
1096  STATE energy = 0.;
1097  for (int i=0; i<3; i++) {
1098  for (int j=0; j<3; j++) {
1099  energy += diffdsol(i,j)*diffstress(i,j);
1100  }
1101  }
1102 
1103 
1104  //values[1] : eror em norma L2
1105  values[1] = L2;
1106 
1107  //values[2] : erro em semi norma H1
1108  values[2] = SemiH1;
1109 
1110  //values[0] : erro em norma H1 <=> norma Energia
1111  values[0] = energy;
1112 
1113 }
1114 
1115 void TPZElasticity3D::ComputeStrainTensor(TPZFMatrix<STATE> &Strain, TPZFMatrix<STATE> &DSol) const{
1116  Strain.Redim(3,3);
1117  Strain(0,0) = DSol(0,0);
1118  Strain(0,1) = 0.5 * ( DSol(1,0) + DSol(0,1) );
1119  Strain(0,2) = 0.5 * ( DSol(2,0) + DSol(0,2) );
1120  Strain(1,0) = Strain(0,1);
1121  Strain(1,1) = DSol(1,1);
1122  Strain(1,2) = 0.5 * ( DSol(2,1) + DSol(1,2) );
1123  Strain(2,0) = Strain(0,2);
1124  Strain(2,1) = Strain(1,2);
1125  Strain(2,2) = DSol(2,2);
1126 }
1127 
1128 void TPZElasticity3D::ComputeStressTensor(TPZFMatrix<STATE> &Stress, TPZMaterialData &data) const{
1129  ComputeStressTensor(Stress, data.dsol[0]);
1130 }
1131 
1132 void TPZElasticity3D::ComputeStressTensor(TPZFMatrix<STATE> &Stress, TPZFMatrix<STATE> &dsol) const{
1133  TPZFNMatrix<6,STATE> Vec(6,1);
1134  this->ComputeStressVector(Vec, dsol);
1135 
1136  Stress.Redim(3,3);
1137  Stress(0,0) = Vec(0,0);
1138  Stress(0,1) = Vec(3,0);
1139  Stress(0,2) = Vec(4,0);
1140  Stress(1,0) = Vec(3,0);
1141  Stress(1,1) = Vec(1,0);
1142  Stress(1,2) = Vec(5,0);
1143  Stress(2,0) = Vec(4,0);
1144  Stress(2,1) = Vec(5,0);
1145  Stress(2,2) = Vec(2,0);
1146 }
1147 
1148 void TPZElasticity3D::ComputeStrainVector(TPZFMatrix<STATE> &Strain, TPZFMatrix<STATE> &DSol) const{
1149  Strain.Redim(6,1);
1150  Strain(0,0) = DSol(0,0);
1151  Strain(1,0) = DSol(1,1);
1152  Strain(2,0) = DSol(2,2);
1153  Strain(3,0) = 0.5 * ( DSol(1,0) + DSol(0,1) );
1154  Strain(4,0) = 0.5 * ( DSol(2,0) + DSol(0,2) );
1155  Strain(5,0) = 0.5 * ( DSol(2,1) + DSol(1,2) );
1156 }
1157 
1158 void TPZElasticity3D::ComputeStressVector(TPZFMatrix<STATE> &Stress, TPZFMatrix<STATE> &DSol) const{
1159  REAL const1 = -1. + this->fPoisson;
1160  REAL const2 = -1. + this->fPoisson + 2. * fPoisson * fPoisson;
1161  const REAL E = this->fE;
1162  const REAL ni = this->fPoisson;
1163  Stress.Redim(6,1);
1164  Stress(0,0) = E * ( DSol(0,0) * const1 - ( DSol(1,1) + DSol(2,2) ) * ni ) / const2 + fPreStress[0];
1165  Stress(1,0) = E * ( DSol(1,1) * const1 - ( DSol(0,0) + DSol(2,2) ) * ni ) / const2 + fPreStress[1];
1166  Stress(2,0) = E * ( DSol(2,2) * const1 - ( DSol(0,0) + DSol(1,1) ) * ni ) / const2 + fPreStress[2];
1167 
1168  REAL const3 = 2. * ( 1. + ni );
1169  Stress(3,0) = E * ( DSol(1,0) + DSol(0,1) ) / const3;
1170  Stress(4,0) = E * ( DSol(2,0) + DSol(0,2) ) / const3;
1171  Stress(5,0) = E * ( DSol(2,1) + DSol(1,2) ) / const3;
1172 }
1173 
1174 void TPZElasticity3D::ApplyDirection(TPZFMatrix<STATE> &StrVec, TPZVec<STATE> &Out) const{
1175  Out.Resize(3);
1176  const TPZVec<REAL> &Dir = this->fPostProcessDirection;
1177  Out[0] = Dir[0] * StrVec(0,0) + Dir[1] * StrVec(3,0) + Dir[2] * StrVec(4,0);
1178  Out[1] = Dir[1] * StrVec(1,0) + Dir[0] * StrVec(3,0) + Dir[2] * StrVec(5,0);
1179  Out[2] = Dir[2] * StrVec(2,0) + Dir[0] * StrVec(4,0) + Dir[1] * StrVec(5,0);
1180 }
1181 
1182 void TPZElasticity3D::PrincipalDirection(TPZFMatrix<STATE> &DSol, TPZVec< STATE > &Solout, int direction) const{
1183 
1184  TPZFNMatrix<9,STATE> StrainTensor(3,3);
1185  TPZManVector<REAL,3> Eigenvalues;
1186  TPZFNMatrix<9,STATE> Eigenvectors(3,3);
1187 
1188  this->ComputeStrainTensor(StrainTensor, DSol);
1189  int64_t numiterations = 1000;
1190  REAL tol = TPZElasticity3D::gTolerance;
1191  bool result;
1192  result = StrainTensor.SolveEigensystemJacobi(numiterations, tol, Solout, Eigenvectors); //Solout is used to store Eigenvaleus, but its values will be replaced below
1193 #ifdef PZDEBUG
1194  if (result == false){
1195  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
1196  }
1197 #endif
1198  Solout.Resize(3);
1199  for(int i = 0; i < 3; i++) Solout[i] = Eigenvectors(direction,i);
1200 }
1201 
1202 void TPZElasticity3D::Invariants( TPZFMatrix<STATE> & A, STATE & I1, STATE & I2, STATE & I3) const
1203 {
1204  I1 = A(0,0)+A(1,1)+A(2,2);
1205  I2 = A(0,0)*A(1,1)-A(1,0)*A(0,1) + A(1,1)*A(2,2)-A(2,1)*A(1,2) + A(0,0)*A(2,2) - A(2,0)*A(0,2);
1206  I3 = A(0,0)*A(1,1)*A(2,2) + A(2,0)*A(0,1)*A(1,2) + A(1,0)*A(2,1)*A(0,2) - ( A(2,0)*A(1,1)*A(0,2) + A(1,0)*A(0,1)*A(2,2) + A(2,1)*A(1,2)*A(0,0) );
1207 }
1208 
1209 void TPZElasticity3D::StressDecomposition( TPZFMatrix<STATE> & A, TPZFMatrix<STATE> & Deviator, STATE & p) const
1210 {
1211  p = (A(0,0)+A(1,1)+A(2,2))/3.;
1212  Deviator = A;
1213  for(int i = 0; i < 3; i++) Deviator(i,i) -= p;
1214 }
1215 
1216 STATE TPZElasticity3D::VonMisesPlasticFunction( TPZFMatrix<STATE> & StressTensor) const
1217 {
1218  TPZFNMatrix<9,STATE> S(3,3);
1219  STATE p = 0.;
1220  this->StressDecomposition(StressTensor, S, p);
1221  STATE J1, J2, J3;
1222  this->Invariants(S, J1, J2, J3);
1223  STATE result = sqrt(3.*fabs(J2))-this->fFy;
1224  return result;
1225 }
1226 
1227 STATE TPZElasticity3D::MohrCoulombPlasticFunction( TPZFMatrix<STATE> & StressTensor) const
1228 {
1229  TPZFNMatrix<9,STATE> S(3,3);
1230  STATE p = 0.;
1231  this->StressDecomposition(StressTensor, S, p);
1232  STATE J1, J2, J3;
1233  this->Invariants(S, J1, J2, J3);
1234  STATE val = -3.*sqrt(3)*J3/(2.*pow(fabs(J2),3./2.));
1235  if(val < -1.) val = -1.;
1236  if(val > +1) val = +1.;
1237  const STATE theta = 1./3.*asin(val);
1238  const REAL c = fCohesion;
1239  const REAL phi = fFrictionAngle;
1240  STATE result = (cos(theta)-1./sqrt(3.)*sin(theta)*sin(phi))*sqrt(fabs(J2))+p*sin(phi)-c*cos(phi);
1241  return result;
1242 }
1243 
1245 void TPZElasticity3D::Write(TPZStream &buf, int withclassid) const
1246 {
1247  TPZMaterial::Write(buf,withclassid);
1248  buf.Write(&fE,1);
1249  buf.Write(&fForce[0],3);
1250  buf.Write(&fFy,1);
1251  DebugStop(); //TODO fFrictionAngle(0.),fCohesion(0.),fPlasticPostProc(
1252  buf.Write(&fPoisson,1);
1253  buf.Write(&fPostProcessDirection[0],3);
1254 
1255 }
1256 
1258 void TPZElasticity3D::Read(TPZStream &buf, void *context)
1259 {
1260  TPZMaterial::Read(buf,context);
1261  buf.Read(&fE,1);
1262  fForce.Resize(3);
1263  buf.Read(&fForce[0],3);
1264  buf.Read(&fFy,1);
1265  buf.Read(&fPoisson,1);
1266  DebugStop(); //TODO fFrictionAngle(0.),fCohesion(0.),fPlasticPostProc(
1267  fPostProcessDirection.Resize(3);
1268  buf.Read(&fPostProcessDirection[0],3);
1269  SetC();
1270 }
1271 
1272 int TPZElasticity3D::ClassId() const{
1273  return Hash("TPZElasticity3D") ^ TPZMaterial::ClassId() << 1;
1274 }
1275 
1276 void TPZElasticity3D::FillDataRequirements(TPZMaterialData &data){
1277 
1278  TPZMaterial::FillDataRequirements(data);
1279  data.fNeedsSol = false;
1280 }
1281 
1282 #ifndef BORLAND
1283 template class TPZRestoreClass<TPZElasticity3D>;
1284 #endif
TPZElasticity3D::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzelast3d.cpp:1245
TPZMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
Definition: TPZMaterial.cpp:200
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZMatrix::SolveEigenvaluesJacobi
virtual bool SolveEigenvaluesJacobi(int64_t &numiterations, REAL &tol, TPZVec< TVar > *Sort=0)
Transforms this matrix in a diagonal matrix, where the diagonal values are its eigenvalues. This method is efficient only for small matrices.
Definition: pzmatrix.cpp:1727
TPZElasticity3D::ComputeStrainVector
void ComputeStrainVector(TPZFMatrix< STATE > &Strain, TPZFMatrix< STATE > &DSol) const
Definition: pzelast3d.cpp:1148
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
TPZElasticity3D::fFy
REAL fFy
Yeilding stress for von mises post processing.
Definition: pzelast3d.h:273
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
TPZMaterial::fNumLoadCases
int fNumLoadCases
Defines the number of load cases generated by this material.
Definition: TPZMaterial.h:76
TPZElasticity3D::PrincipalDirection
void PrincipalDirection(TPZFMatrix< STATE > &DSol, TPZVec< STATE > &Solout, int direction) const
Definition: pzelast3d.cpp:1182
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZElasticity3D::ContributeVecShapeBC
virtual void ContributeVecShapeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzelast3d.cpp:509
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZElasticity3D::Solution
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
Post-processing method. Based on solution Sol and its derivatives DSol, it computes the post-processe...
Definition: pzelast3d.cpp:831
TPZAxesTools::Axes2XYZ
static void Axes2XYZ(const TPZFMatrix< TVar > &dudaxes, TPZFMatrix< TVar > &dudx, const TPZFMatrix< REAL > &axesv, bool colMajor=true)
Makes the basis transformation from axes basis to euclidian basis.
Definition: pzaxestools.h:61
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZMatrix::SolveEigensystemJacobi
virtual bool SolveEigensystemJacobi(int64_t &numiterations, REAL &tol, TPZVec< TVar > &Eigenvalues, TPZFMatrix< TVar > &Eigenvectors) const
Compute Eigenvalues and Eigenvectors of this matrix. This method is efficient only for small matrice...
Definition: pzmatrix.cpp:1583
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
pzerror.h
Defines PZError.
TPZElasticity3D
This class implements a 3D isotropic elasticity material.
Definition: pzelast3d.h:21
TPZElasticity3D::fPlasticPostProc
PLASTICPOSTPROC fPlasticPostProc
Plastic model for post-processing.
Definition: pzelast3d.h:279
TPZElasticity3D::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
It computes a contribution to the stiffness matrix and load vector at one integration point...
Definition: pzelast3d.cpp:88
TPZElasticity3D::ContributeVecShape
virtual void ContributeVecShape(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: pzelast3d.cpp:417
TPZMaterial::Write
void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: TPZMaterial.cpp:387
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZMaterialData::fShapeType
MShapeFunctionType fShapeType
Definition: pzmaterialdata.h:43
TPZElasticity3D::gTolerance
static STATE gTolerance
Definition: pzelast3d.h:304
TPZElasticity3D::ApplyDirection
void ApplyDirection(TPZFMatrix< STATE > &StrVec, TPZVec< STATE > &Out) const
Definition: pzelast3d.cpp:1174
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZElasticity3D::Invariants
void Invariants(TPZFMatrix< STATE > &A, STATE &I1, STATE &I2, STATE &I3) const
Definition: pzelast3d.cpp:1202
TPZElasticity3D::ComputeStressTensor
virtual void ComputeStressTensor(TPZFMatrix< STATE > &Stress, TPZMaterialData &data) const
Definition: pzelast3d.cpp:1128
TPZElasticity3D::StressDecomposition
void StressDecomposition(TPZFMatrix< STATE > &StressTensor, TPZFMatrix< STATE > &Deviator, STATE &p) const
Definition: pzelast3d.cpp:1209
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
TPZElasticity3D::fFrictionAngle
REAL fFrictionAngle
Mohr-Coulomb parameters.
Definition: pzelast3d.h:276
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
TPZElasticity3D::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: pzelast3d.cpp:1272
pzelast3d.h
Contains the TPZElasticity3D class which implements a 3D isotropic elasticity material.
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZMaterial::Print
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Definition: TPZMaterial.cpp:103
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...
TPZElasticity3D::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Implements Dirichlet and Neumann boundary conditions.
Definition: pzelast3d.cpp:612
TPZMaterial::ForcingFunction
TPZAutoPointer< TPZFunction< STATE > > & ForcingFunction()
Returns a procedure as source function for the material.
Definition: TPZMaterial.h:397
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
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
sin
sin
Definition: fadfunc.h:63
TPZElasticity3D::Print
virtual void Print(std::ostream &out) override
Print material report.
Definition: pzelast3d.cpp:77
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
TPZMaterial::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: TPZMaterial.cpp:402
TPZElasticity3D::fPostProcessDirection
TPZManVector< REAL, 3 > fPostProcessDirection
Direction to compute stress and strain.
Definition: pzelast3d.h:270
TPZElasticity3D::fPreStress
TPZManVector< REAL > fPreStress
Definition: pzelast3d.h:281
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZElasticity3D::~TPZElasticity3D
virtual ~TPZElasticity3D()
Default destructor.
Definition: pzelast3d.cpp:64
TPZElasticity3D::FillDataRequirements
void FillDataRequirements(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the Contribute method.
Definition: pzelast3d.cpp:1276
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TPZMaterial::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data)
Fill material data parameter with necessary requirements for the.
Definition: TPZMaterial.cpp:81
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
pzmanvector.h
Free store vector implementation.
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZElasticity3D::MohrCoulombPlasticFunction
STATE MohrCoulombPlasticFunction(TPZFMatrix< STATE > &StressTensor) const
Definition: pzelast3d.cpp:1227
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
TPZMaterialData::axes
TPZFNMatrix< 9, REAL > axes
axes indicating the directions of the derivatives of the shapefunctions
Definition: pzmaterialdata.h:63
test.res
string res
Definition: test.py:151
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZMaterial::HasForcingFunction
virtual int HasForcingFunction()
Directive that gives true if the material has a forcing function.
Definition: TPZMaterial.h:472
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZMaterial::gBigNumber
static REAL gBigNumber
Big number to penalization method, used for Dirichlet conditions.
Definition: TPZMaterial.h:83
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZElasticity3D::Read
virtual void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzelast3d.cpp:1258
TPZElasticity3D::ComputeStrainTensor
void ComputeStrainTensor(TPZFMatrix< STATE > &Strain, TPZFMatrix< STATE > &DSol) const
Definition: pzelast3d.cpp:1115
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZElasticity3D::ComputeStressVector
virtual void ComputeStressVector(TPZFMatrix< STATE > &Stress, TPZFMatrix< STATE > &DSol) const
Definition: pzelast3d.cpp:1158
TPZMaterial::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: TPZMaterial.cpp:382
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
idf
static int idf[6]
Definition: TPZCompMeshTools.cpp:135
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZMatrix::VerifySymmetry
virtual int VerifySymmetry(REAL tol=1.e-13) const
Checks if current matrix value is symmetric.
Definition: pzmatrix.cpp:1501
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZElasticity3D::fPoisson
STATE fPoisson
Poisson&#39;s ratio.
Definition: pzelast3d.h:261
asin
expr_ expr_ expr_ expr_ expr_ expr_ asin
Definition: tfadfunc.h:80
TPZVec::Fill
void Fill(const T &copy, const int64_t from=0, const int64_t numelem=-1)
Will fill the elements of the vector with a copy object.
Definition: pzvec.h:460
TPZElasticity3D::fForce
TPZManVector< STATE, 3 > fForce
External forces.
Definition: pzelast3d.h:267
pzaxestools.h
Contains declaration of the TPZAxesTools class which implements verifications over axes...
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZElasticity3D::fE
STATE fE
Young&#39;s modulus.
Definition: pzelast3d.h:258
TPZVec::NElements
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
TPZElasticity3D::Errors
virtual void Errors(TPZVec< REAL > &x, TPZVec< STATE > &u, TPZFMatrix< STATE > &dudx, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux, TPZVec< STATE > &u_exact, TPZFMatrix< STATE > &du_exact, TPZVec< REAL > &values) override
Evaluate error between approximate (FEM) and exact solutions.
Definition: pzelast3d.cpp:1074
rdt.values
def values
Definition: rdt.py:119
TPZElasticity3D::VonMisesPlasticFunction
STATE VonMisesPlasticFunction(TPZFMatrix< STATE > &StressTensor) const
Definition: pzelast3d.cpp:1216
TPZElasticity3D::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns index of post-processing variable.
Definition: pzelast3d.cpp:771
TPZElasticity3D::TPZElasticity3D
TPZElasticity3D()
Default constructor.
Definition: pzelast3d.cpp:57
cos
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZElasticity3D::NSolutionVariables
virtual int NSolutionVariables(int var) override
Number of data of variable var.
Definition: pzelast3d.cpp:798
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
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
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13