NeoPZ
pzelasmat.cpp
Go to the documentation of this file.
1 
6 #include "pzelasmat.h"
7 #include "pzelmat.h"
8 #include "pzbndcond.h"
9 #include "pzmatrix.h"
10 #include "pzfmatrix.h"
11 #include "pzerror.h"
12 #include "pzaxestools.h"
13 
14 #include <math.h>
15 
16 #include "pzlog.h"
17 #ifdef LOG4CXX
18 static LoggerPtr logdata(Logger::getLogger("pz.material.elasticity.data"));
19 #endif
20 
21 #include <fstream>
22 using namespace std;
23 
24 TPZElasticityMaterial::TPZElasticityMaterial() :
25 TPZRegisterClassId(&TPZElasticityMaterial::ClassId),
26 TPZDiscontinuousGalerkin(0), ff(3,0.) {
27  fE_def = -1.; // Young modulus
28  fnu_def = -1.; // poisson coefficient
29  ff[0] = 0.; // X component of the body force
30  ff[1] = 0.; // Y component of the body force
31  ff[2] = 0.; // Z component of the body force - not used for this class
32 
33 
34  //Added by Cesar 2001/03/16
35  fPreStressXX = 0.; //Prestress in the x direction
36  fPreStressYY = 0.; //Prestress in the y direction
37  fPreStressXY = 0.; //Prestress in the z direction
38  fPreStressZZ = 0.; //Prestress in the z direction
39  fPlaneStress = 0;
40 
41  // Added by Philippe 2012
42  fPostProcIndex = 0;
43 }
44 
45 TPZElasticityMaterial::TPZElasticityMaterial(int id) : TPZRegisterClassId(&TPZElasticityMaterial::ClassId),
46 TPZDiscontinuousGalerkin(id), ff(3,0.) {
47  fE_def = -1.; // Young modulus
48  fnu_def = -1.; // poisson coefficient
49  ff[0] = 0.; // X component of the body force
50  ff[1] = 0.; // Y component of the body force
51  ff[2] = 0.; // Z component of the body force - not used for this class
52 
53 
54  //Added by Cesar 2001/03/16
55  fPreStressXX = 0.; //Prestress in the x direction
56  fPreStressYY = 0.; //Prestress in the y direction
57  fPreStressXY = 0.; //Prestress in the z direction
58  fPreStressZZ = 0.; //Prestress in the z direction
59  fPlaneStress = 0;
60 
61  // Added by Philippe 2012
62  fPostProcIndex = 0;
63 }
64 
65 TPZElasticityMaterial::TPZElasticityMaterial(int num, REAL E, REAL nu, REAL fx, REAL fy, int plainstress) :
66 TPZRegisterClassId(&TPZElasticityMaterial::ClassId),
67 TPZDiscontinuousGalerkin(num), ff(3,0.) {
68 
69  fE_def = E; // Young modulus
70  fnu_def = nu; // poisson coefficient
71  ff[0] = fx; // X component of the body force
72  ff[1] = fy; // Y component of the body force
73  ff[2] = 0.; // Z component of the body force - not used for this class
74 
75  //Added by Cesar 2001/03/16
76  fPreStressXX = 0.; //Prestress in the x direction
77  fPreStressYY = 0.; //Prestress in the y direction
78  fPreStressXY = 0.; //Prestress in the z direction
79  fPreStressZZ = 0.; //Prestress in the z direction
80  fPlaneStress = plainstress;
81  // Added by Philippe 2012
82  fPostProcIndex = 0;
83 }
84 
85 TPZElasticityMaterial::~TPZElasticityMaterial() {
86 }
87 
88 int TPZElasticityMaterial::NStateVariables() const {
89  return 2;
90 }
91 
92 void TPZElasticityMaterial::Print(std::ostream &out) {
93  out << "name of material : " << Name() << "\n";
94  out << "properties : \n";
95  if(fElasticity)
96  {
97  out << "Elasticity coeficients are determined by a function\n";
98  }
99  else
100  {
101  out << "\tE = " << fE_def << endl;
102  out << "\tnu = " << fnu_def << endl;
103  }
104  out << "\tF = " << ff[0] << ' ' << ff[1] << endl;
105  out << "\t PreStress: \n"
106  << "Sigma xx = \t" << fPreStressXX << "\t"
107  << "Sigma yy = \t" << fPreStressYY << "\t"
108  << "Sigma xy = \t" << fPreStressXY << "Sigma zz = \t" << fPreStressZZ << endl;
109 }
110 
111 //Added by Cesar 2001/03/16
112 void TPZElasticityMaterial::SetPreStress(REAL Sigxx, REAL Sigyy, REAL Sigxy, REAL Sigzz){
113  fPreStressXX = Sigxx;
114  fPreStressYY = Sigyy;
115  fPreStressXY = Sigxy;
116  fPreStressZZ = Sigzz;
117 }
118 
119 
120 void TPZElasticityMaterial::Contribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef) {
121 
122  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
123  if(shapetype==data.EVecShape){
124  ContributeVecShape(data,weight,ek, ef);
125  return;
126  }
127 
128  TPZFMatrix<REAL> &dphi = data.dphix;
129  TPZFMatrix<REAL> &phi = data.phi;
130  TPZFMatrix<REAL> &axes=data.axes;
131 
132  int phc,phr,dphc,dphr,efr,efc,ekr,ekc;
133  phc = phi.Cols();
134  phr = phi.Rows();
135  dphc = dphi.Cols();
136  dphr = dphi.Rows();
137  efr = ef.Rows();
138  efc = ef.Cols();
139  ekr = ek.Rows();
140  ekc = ek.Cols();
141  if(phc != 1 || dphr != 2 || phr != dphc){
142  PZError << "\nTPZElasticityMaterial.contr, inconsistent input data : \n" <<
143  "phi.Cols() = " << phi.Cols() << " dphi.Cols() = " << dphi.Cols() <<
144  " phi.Rows = " << phi.Rows() << " dphi.Rows = " <<
145  dphi.Rows() << "\nek.Rows() = " << ek.Rows() << " ek.Cols() = "
146  << ek.Cols() <<
147  "\nef.Rows() = " << ef.Rows() << " ef.Cols() = "
148  << ef.Cols() << "\n";
149  return;
150  // PZError.show();
151  }
152  TPZManVector<STATE,3> floc(ff);
153  if(fForcingFunction) { // phi(in, 0) : node in associated forcing function
154  TPZManVector<STATE,3> res(3,0.);
155  fForcingFunction->Execute(data.x,res);
156  floc[0] = res[0];
157  floc[1] = res[1];
158  floc[2] = res[2];
159  }
160 
161  REAL E(fE_def), nu(fnu_def);
162 
163  if (fElasticity) {
164  TPZManVector<STATE,2> result(2);
165  TPZFNMatrix<4,STATE> Dres(0,0);
166  fElasticity->Execute(data.x, result, Dres);
167  E = result[0];
168  nu = result[1];
169  }
170 
171  REAL Eover1MinNu2 = E/(1-nu*nu);
172  REAL Eover21PlusNu = E/(2.*(1+nu));
173 
174  TPZFNMatrix<4,STATE> du(2,2);
175  /*
176  * Plain strain materials values
177  */
178  REAL nu1 = 1. - nu;//(1-nu)
179  REAL nu2 = (1.-2.*nu)/2.;
180  REAL F = E/((1.+nu)*(1.-2.*nu));
181 
182  for( int in = 0; in < phr; in++ ) {
183  du(0,0) = dphi(0,in)*axes(0,0)+dphi(1,in)*axes(1,0);//dvx
184  du(1,0) = dphi(0,in)*axes(0,1)+dphi(1,in)*axes(1,1);//dvy
185 
186  for (int col = 0; col < efc; col++)
187  {
188  ef(2*in, col) += weight * (floc[0]*phi(in,0) - du(0,0)*fPreStressXX - du(1,0)*fPreStressXY); // direcao x
189  ef(2*in+1, col) += weight * (floc[1]*phi(in,0) - du(0,0)*fPreStressXY - du(1,0)*fPreStressYY);// direcao y <<<----
190  }
191  for( int jn = 0; jn < phr; jn++ ) {
192  du(0,1) = dphi(0,jn)*axes(0,0)+dphi(1,jn)*axes(1,0);//dux
193  du(1,1) = dphi(0,jn)*axes(0,1)+dphi(1,jn)*axes(1,1);//duy
194 
195 
196  if (fPlaneStress != 1){
197  /* Plane Strain State */
198  ek(2*in,2*jn) += weight * (
199  nu1 * du(0,0)*du(0,1)+ nu2 * du(1,0)*du(1,1)
200  ) * F;
201 
202  ek(2*in,2*jn+1) += weight * (
203  nu*du(0,0)*du(1,1)+ nu2*du(1,0)*du(0,1)
204  ) * F;
205 
206  ek(2*in+1,2*jn) += weight * (
207  nu*du(1,0)*du(0,1)+ nu2*du(0,0)*du(1,1)
208  ) * F;
209 
210  ek(2*in+1,2*jn+1) += weight * (
211  nu1*du(1,0)*du(1,1)+ nu2*du(0,0)*du(0,1)
212  ) * F;
213  }
214  else{
215  /* Plain stress state */
216  ek(2*in,2*jn) += weight * (
217  Eover1MinNu2 * du(0,0)*du(0,1)+ Eover21PlusNu * du(1,0)*du(1,1)
218  );
219 
220  ek(2*in,2*jn+1) += weight * (
221  Eover1MinNu2*nu*du(0,0)*du(1,1)+ Eover21PlusNu*du(1,0)*du(0,1)
222  );
223 
224  ek(2*in+1,2*jn) += weight * (
225  Eover1MinNu2*nu*du(1,0)*du(0,1)+ Eover21PlusNu*du(0,0)*du(1,1)
226  );
227 
228  ek(2*in+1,2*jn+1) += weight * (
229  Eover1MinNu2*du(1,0)*du(1,1)+ Eover21PlusNu*du(0,0)*du(0,1)
230  );
231  }
232  }
233  }
234 
235 //#ifdef LOG4CXX
236 // if(logdata->isDebugEnabled())
237 // {
238 // std::stringstream sout;
239 // ek.Print("ek_elastmat = ",sout,EMathematicaInput);
240 // ef.Print("ef_elastmat = ",sout,EMathematicaInput);
241 // LOGPZ_DEBUG(logdata,sout.str())
242 // }
243 //#endif
244 
245 }
246 
247 void TPZElasticityMaterial::Contribute(TPZVec<TPZMaterialData> &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef) {
248 
249  TPZMaterialData::MShapeFunctionType shapetype = data[0].fShapeType;
250  if(shapetype==data[0].EVecShape){
251  DebugStop();
252  return;
253  }
254  TPZFMatrix<REAL> &dphi = data[0].dphix;
255  TPZFMatrix<REAL> &phi = data[0].phi;
256  TPZFMatrix<REAL> &axes=data[0].axes;
257 
258  int phc,phr,dphc,dphr,efr,efc,ekr,ekc;
259  phc = phi.Cols();
260  phr = phi.Rows();
261  dphc = dphi.Cols();
262  dphr = dphi.Rows();
263  efr = ef.Rows();
264  efc = ef.Cols();
265  ekr = ek.Rows();
266  ekc = ek.Cols();
267  if(phc != 1 || dphr != 2 || phr != dphc || 2*phr != ekr){
268  PZError << "\nTPZElasticityMaterial.contr, inconsistent input data : \n" <<
269  "phi.Cols() = " << phi.Cols() << " dphi.Cols() = " << dphi.Cols() <<
270  " phi.Rows = " << phi.Rows() << " dphi.Rows = " <<
271  dphi.Rows() << "\nek.Rows() = " << ek.Rows() << " ek.Cols() = "
272  << ek.Cols() <<
273  "\nef.Rows() = " << ef.Rows() << " ef.Cols() = "
274  << ef.Cols() << "\n";
275  DebugStop();
276  return;
277  // PZError.show();
278  }
279  if(fForcingFunction) { // phi(in, 0) : node in associated forcing function
280  TPZManVector<STATE,3> res(3);
281  fForcingFunction->Execute(data[0].x,res);
282  ff[0] = res[0];
283  ff[1] = res[1];
284  ff[2] = res[2];
285  }
286 
287  REAL E(fE_def), nu(fnu_def);
288 
289  if (fElasticity) {
290  TPZManVector<STATE,2> result(2);
291  TPZFNMatrix<4,STATE> Dres(0,0);
292  fElasticity->Execute(data[0].x, result, Dres);
293  E = result[0];
294  nu = result[1];
295  }
296 
297  REAL Eover1MinNu2 = E/(1-nu*nu);
298  REAL Eover21PlusNu = E/(2.*(1+nu));
299 
300 
301  TPZFNMatrix<4,STATE> du(2,2);
302  /*
303  * Plane strain materials values
304  */
305  REAL nu1 = 1. - nu;//(1-nu)
306  REAL nu2 = (1.-2.*nu)/2.;
307  REAL F = E/((1.+nu)*(1.-2.*nu));
308  STATE epsx, epsy,epsxy,epsz;
309  TPZFNMatrix<9,STATE> DSolxy(2,2);
310  // dudx - dudy
311  DSolxy(0,0) = data[0].dsol[0](0,0)*axes(0,0)+data[0].dsol[0](1,0)*axes(1,0);
312  DSolxy(1,0) = data[0].dsol[0](0,0)*axes(0,1)+data[0].dsol[0](1,0)*axes(1,1);
313  // dvdx - dvdy
314  DSolxy(0,1) = data[0].dsol[0](0,1)*axes(0,0)+data[0].dsol[0](1,1)*axes(1,0);
315  DSolxy(1,1) = data[0].dsol[0](0,1)*axes(0,1)+data[0].dsol[0](1,1)*axes(1,1);
316  epsx = DSolxy(0,0);// du/dx
317  epsy = DSolxy(1,1);// dv/dy
318  epsxy = 0.5*(DSolxy(1,0)+DSolxy(0,1));
319  epsz = data[1].sol[0][0];
320  STATE SigX = E/((1.-2.*nu)*(1.+nu))*((1.-nu)*epsx+nu*(epsy+epsz))+fPreStressXX;
321  STATE SigY = E/((1.-2.*nu)*(1.+nu))*(nu*epsx+(1.-nu)*(epsy+epsz))+fPreStressYY;
322  REAL lambda = GetLambda(E,nu);
323  REAL mu = GetMU(E,nu);
324 
325  STATE SigZ = lambda*(epsx+epsy+epsz)+2.*mu*epsz;
326  STATE TauXY = 2*mu*epsxy+fPreStressXY;
327 
328 
329  for( int in = 0; in < phr; in++ ) {
330  du(0,0) = dphi(0,in)*axes(0,0)+dphi(1,in)*axes(1,0);//dvx
331  du(1,0) = dphi(0,in)*axes(0,1)+dphi(1,in)*axes(1,1);//dvy
332 
333 
334  for (int col = 0; col < efc; col++)
335  {
336  ef(2*in, col) += weight * (ff[0]*phi(in,0) - du(0,0)*(SigX) - du(1,0)*(TauXY)); // direcao x
337  ef(2*in+1, col) += weight * (ff[1]*phi(in,0) - du(0,0)*(TauXY) - du(1,0)*(SigY)); // direcao y <<<----
338  }
339  for( int jn = 0; jn < phr; jn++ ) {
340  du(0,1) = dphi(0,jn)*axes(0,0)+dphi(1,jn)*axes(1,0);//dux
341  du(1,1) = dphi(0,jn)*axes(0,1)+dphi(1,jn)*axes(1,1);//duy
342 
343 
344  if (fPlaneStress != 1){
345  /* Plane Strain State */
346  ek(2*in,2*jn) += weight * (
347  nu1 * du(0,0)*du(0,1)+ nu2 * du(1,0)*du(1,1)
348  ) * F;
349 
350  ek(2*in,2*jn+1) += weight * (
351  nu*du(0,0)*du(1,1)+ nu2*du(1,0)*du(0,1)
352  ) * F;
353 
354  ek(2*in+1,2*jn) += weight * (
355  nu*du(1,0)*du(0,1)+ nu2*du(0,0)*du(1,1)
356  ) * F;
357 
358  ek(2*in+1,2*jn+1) += weight * (
359  nu1*du(1,0)*du(1,1)+ nu2*du(0,0)*du(0,1)
360  ) * F;
361  }
362  else{
363  /* Plain stress state */
364  ek(2*in,2*jn) += weight * (
365  Eover1MinNu2 * du(0,0)*du(0,1)+ Eover21PlusNu * du(1,0)*du(1,1)
366  );
367 
368  ek(2*in,2*jn+1) += weight * (
369  Eover1MinNu2*nu*du(0,0)*du(1,1)+ Eover21PlusNu*du(1,0)*du(0,1)
370  );
371 
372  ek(2*in+1,2*jn) += weight * (
373  Eover1MinNu2*nu*du(1,0)*du(0,1)+ Eover21PlusNu*du(0,0)*du(1,1)
374  );
375 
376  ek(2*in+1,2*jn+1) += weight * (
377  Eover1MinNu2*du(1,0)*du(1,1)+ Eover21PlusNu*du(0,0)*du(0,1)
378  );
379  }
380  }
381  const STATE C2 = E * nu / (-1. + nu + 2.*nu*nu);
382  const int nstate = 2;
383 
384 // ek(in*nstate,phr*nstate) += weight*(-C2*dphiXY(0,in));
385 // ek(in*nstate+1,phr*nstate) += weight*(-C2*dphiXY(1,in));
386 //
387 // ek(phr*nstate,in*nstate) += weight*(-C2*dphiXY(0,in));
388 // ek(phr*nstate,in*nstate+1) += weight*(-C2*dphiXY(1,in));
389  ek(in*nstate,phr*nstate) += weight*(-C2*du(0,0));
390  ek(in*nstate+1,phr*nstate) += weight*(-C2*du(1,0));
391 
392  ek(phr*nstate,in*nstate) += weight*(-C2*du(0,0));
393  ek(phr*nstate,in*nstate+1) += weight*(-C2*du(1,0));
394 
395  }
396  ek(phr*2,phr*2) += weight*(lambda+2.*mu);
397  for (int col = 0; col < efc; col++)
398  {
399  ef(2*phr, col) += weight * (-SigZ); // direcao z
400  }
401 
402 
403  //#ifdef LOG4CXX
404  // if(logdata->isDebugEnabled())
405  // {
406  // std::stringstream sout;
407  // ek.Print("ek_elastmat = ",sout,EMathematicaInput);
408  // ef.Print("ef_elastmat = ",sout,EMathematicaInput);
409  // LOGPZ_DEBUG(logdata,sout.str())
410  // }
411  //#endif
412 
413 }
414 
415 
416 
417 void TPZElasticityMaterial::FillDataRequirements(TPZMaterialData &data)
418 {
419  data.fNeedsSol = true;
420  data.fNeedsNormal = false;
421 
422 }
423 
424 void TPZElasticityMaterial::FillBoundaryConditionDataRequirement(int type,TPZMaterialData &data)
425 {
426  data.fNeedsSol = false;
427  data.fNeedsNormal = false;
428  if (type == 4 || type == 5 || type == 6) {
429  data.fNeedsNormal = true;
430  }
431 }
432 
433 void TPZElasticityMaterial::ContributeVecShape(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
434 {
435  TPZFMatrix<REAL> &dphi = data.dphix;
436  TPZFMatrix<REAL> &phi = data.phi;
437  TPZFMatrix<REAL> &axes=data.axes;
438 
439  int phc,phr,dphc,dphr,efr,efc,ekr,ekc;
440  phc = phi.Cols();
441  phr = phi.Rows();
442  dphc = dphi.Cols();
443  dphr = dphi.Rows();
444  efr = ef.Rows();
445  efc = ef.Cols();
446  ekr = ek.Rows();
447  ekc = ek.Cols();
448 
449  if(fForcingFunction) { // phi(in, 0) : node in associated forcing function
450  TPZManVector<STATE> res(3);
451  fForcingFunction->Execute(data.x,res);
452  ff[0] = res[0];
453  ff[1] = res[1];
454  ff[2] = res[2];
455  }
456 
457  REAL E(fE_def), nu(fnu_def);
458 
459  if (fElasticity) {
460  TPZManVector<STATE,2> result(2);
461  TPZFNMatrix<4,STATE> Dres(0,0);
462  fElasticity->Execute(data.x, result, Dres);
463  E = result[0];
464  nu = result[1];
465  }
466 
467  REAL Eover1MinNu2 = E/(1-nu*nu);
468  REAL Eover21PlusNu = E/(2.*(1+nu));
469 
470 
471  TPZFNMatrix<4,STATE> dphix_i(2,1),dphiy_i(2,1), dphix_j(2,1), dphiy_j(2,1);
472  /*
473  * Plain strain materials values
474  */
475  REAL nu1 = 1 - nu;//(1-nu)
476  REAL nu2 = (1-2*nu)/2;
477  REAL F = E/((1+nu)*(1-2*nu));
478 
479  for( int in = 0; in < phc; in++ )
480  {
481  dphix_i(0,0) = dphi(0,in)*axes(0,0)+dphi(1,in)*axes(1,0);
482  dphix_i(1,0) = dphi(0,in)*axes(0,1)+dphi(1,in)*axes(1,1);
483  dphiy_i(0,0) = dphi(2,in)*axes(0,0)+dphi(3,in)*axes(1,0);
484  dphiy_i(1,0) = dphi(2,in)*axes(0,1)+dphi(3,in)*axes(1,1);
485 
486  for (int col = 0; col < efc; col++)
487  {
488  ef(in,col) += weight*( ff[0] * phi(0, in)- dphix_i(0,0)*fPreStressXX - dphix_i(1,0)*fPreStressXY
489  + ff[1] * phi(1, in)- dphiy_i(0,0)*fPreStressYY - dphiy_i(1,0)*fPreStressXY);
490  }
491  for( int jn = 0; jn < phc; jn++ ) {
492 
493  dphix_j(0,0) = dphi(0,jn)*axes(0,0)+dphi(1,jn)*axes(1,0);
494  dphix_j(1,0) = dphi(0,jn)*axes(0,1)+dphi(1,jn)*axes(1,1);
495  dphiy_j(0,0) = dphi(2,jn)*axes(0,0)+dphi(3,jn)*axes(1,0);
496  dphiy_j(1,0) = dphi(2,jn)*axes(0,1)+dphi(3,jn)*axes(1,1);
497 
498 
499  if (fPlaneStress != 1){
500  /* Plane Strain State */
501  ek(in,jn) += weight*(nu1*dphix_i(0,0)*dphix_j(0,0) + nu2*dphix_i(1,0)*dphix_j(1,0) +
502 
503  nu*dphix_i(0,0)*dphiy_j(1,0) + nu2*dphix_i(1,0)*dphiy_j(0,0) +
504 
505  nu*dphiy_i(1,0)*dphix_j(0,0) + nu2*dphiy_i(0,0)*dphix_j(1,0) +
506 
507  nu1*dphiy_i(1,0)*dphiy_j(1,0) + nu2*dphiy_i(0,0)*dphiy_j(0,0))*F;
508  }
509  else{
510  /* Plain stress state */
511 
512  ek(in,jn) += weight*(Eover1MinNu2*dphix_i(0,0)*dphix_j(0,0) + Eover21PlusNu*dphix_i(1,0)*dphix_j(1,0) +
513 
514  Eover1MinNu2*dphix_i(0,0)*dphiy_j(1,0) + Eover21PlusNu*dphix_i(1,0)*dphiy_j(0,0) +
515 
516  Eover1MinNu2*dphiy_i(1,0)*dphix_j(0,0) + Eover21PlusNu*dphiy_i(0,0)*dphix_j(1,0) +
517 
518  Eover1MinNu2*dphiy_i(1,0)*dphiy_j(1,0) + Eover21PlusNu*dphiy_i(0,0)*dphiy_j(0,0));
519  }
520  }
521  }
522 }
523 
524 void TPZElasticityMaterial::ContributeBC(TPZVec<TPZMaterialData> &datavec, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, TPZBndCond &bc) {
525  bc.Contribute(datavec[0],weight,ek,ef);
526 }
527 
528 void TPZElasticityMaterial::ContributeBC(TPZMaterialData &data,REAL weight,
529  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
530 
531 
532  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
533  if(shapetype==data.EVecShape){
534  ContributeVecShapeBC(data,weight,ek, ef,bc);
535  return;
536  }
537 
538  TPZFMatrix<REAL> &phi = data.phi;
539  int dim = Dimension();
540 
541  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
542 
543  int phr = phi.Rows();
544  short in,jn;
545 
546  if (ef.Cols() != bc.NumLoadCases()) {
547  DebugStop();
548  }
549 
550 // In general when the problem is needed to stablish any convention for ContributeBC implementations
551 
552 // REAL v2[2];
553 // v2[0] = bc.Val2()(0,0);
554 // v2[1] = bc.Val2()(1,0);
555  int nstate = NStateVariables();
556 
557  TPZFMatrix<STATE> &v1 = bc.Val1();
558 
559 
560  switch (bc.Type()) {
561  case 0 : // Dirichlet condition
562  {
563  for(in = 0 ; in < phr; in++) {
564  for (int il = 0; il<NumLoadCases(); il++)
565  {
566  REAL v2[2];
567  v2[0] = bc.Val2(il)(0,0);
568  v2[1] = bc.Val2(il)(1,0);
569  ef(2*in,il) += BIGNUMBER * v2[0] * phi(in,0) * weight; // forced v2 displacement
570  ef(2*in+1,il) += BIGNUMBER * v2[1] * phi(in,0) * weight; // forced v2 displacement
571  }
572  for (jn = 0 ; jn < phi.Rows(); jn++)
573  {
574  ek(2*in,2*jn) += BIGNUMBER * phi(in,0) *phi(jn,0) * weight;
575  ek(2*in+1,2*jn+1) += BIGNUMBER * phi(in,0) *phi(jn,0) * weight;
576  }
577  }
578  }
579  break;
580 
581  case 1 : // Neumann condition
582  {
583  for (in = 0; in < phr; in++)
584  {
585  for (int il = 0; il <fNumLoadCases; il++)
586  {
587  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
588  ef(2*in,il) += v2(0,0) * phi(in,0) * weight; // force in x direction
589  ef(2*in+1,il) += v2(1,0) * phi(in,0) * weight; // force in y direction
590  }
591  }
592  }
593  break;
594 
595  case 2 : // Mixed Condition
596  {
597  for(in = 0 ; in < phi.Rows(); in++)
598  {
599  for (int il = 0; il <fNumLoadCases; il++)
600  {
601  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
602  ef(2*in,il) += v2(0,0) * phi(in,0) * weight; // force in x direction
603  ef(2*in+1,il) += v2(1,0) * phi(in,0) * weight; // forced in y direction
604  }
605 
606  for (jn = 0 ; jn < phi.Rows(); jn++) {
607  ek(2*in,2*jn) += bc.Val1()(0,0) * phi(in,0) * phi(jn,0) * weight; // peso de contorno => integral de contorno
608  ek(2*in+1,2*jn) += bc.Val1()(1,0) * phi(in,0) * phi(jn,0) * weight;
609  ek(2*in+1,2*jn+1) += bc.Val1()(1,1) * phi(in,0) * phi(jn,0) * weight;
610  ek(2*in,2*jn+1) += bc.Val1()(0,1) * phi(in,0) * phi(jn,0) * weight;
611  }
612  } // este caso pode reproduzir o caso 0 quando o deslocamento
613 
614  break;
615  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
616  for(in = 0 ; in < phr; in++) {
617 // ef(nstate*in+0,0) += BIGNUMBER * (0. - data.sol[0][0]) * v2[0] * phi(in,0) * weight;
618 // ef(nstate*in+1,0) += BIGNUMBER * (0. - data.sol[0][1]) * v2[1] * phi(in,0) * weight;
619  for (jn = 0 ; jn < phr; jn++) {
620  ek(nstate*in+0,nstate*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * bc.Val2()(0,0);
621  ek(nstate*in+1,nstate*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * bc.Val2()(1,0);
622  }//jn
623  }//in
624  break;
625 
626 
627  case 4: // stressField Neumann condition
628  {
629  REAL v2[2];
630  for(in = 0; in < dim; in ++)
631  {
632  v2[in] = ( v1(in,0) * data.normal[0] +
633  v1(in,1) * data.normal[1]);
634  }
635  // The normal vector points towards the neighbour. The negative sign is there to
636  // reflect the outward normal vector.
637  for(in = 0 ; in < phi.Rows(); in++) {
638  ef(nstate*in+0,0) += v2[0] * phi(in,0) * weight;
639  ef(nstate*in+1,0) += v2[1] * phi(in,0) * weight;
640  // cout << "normal:" << data.normal[0] << ' ' << data.normal[1] << ' ' << data.normal[2] << endl;
641  // cout << "val2: " << v2[0] << endl;
642  }
643  }
644  break;
645 
646  case 5://PRESSAO DEVE SER POSTA NA POSICAO 0 DO VETOR v2
647  {
648  TPZFNMatrix<2,STATE> res(2,1,0.);
649  for(in = 0 ; in < phi.Rows(); in++)
650  {
651  for (int il=0; il<NumLoadCases(); il++)
652  {
653  ef(nstate*in+0,0) += (bc.Val2(il)(0,0)*data.normal[0]) * phi(in,0) * weight ;
654  ef(nstate*in+1,0) += (bc.Val2(il)(0,0)*data.normal[1]) * phi(in,0) * weight ;
655  }
656  for(jn=0; jn<phi.Rows(); jn++)
657  {
658  for(int idf=0; idf<2; idf++) for(int jdf=0; jdf<2; jdf++)
659  {
660  ek(nstate*in+idf,nstate*jn+jdf) += bc.Val1()(idf,jdf)*data.normal[idf]*data.normal[jdf]*phi(in,0)*phi(jn,0)*weight;
661  //BUG FALTA COLOCAR VAL2
662  // DebugStop();
663  }
664  }
665 
666  }
667  }
668  break;
669 
670  case 6://PRESSAO DEVE SER POSTA NA POSICAO 0 DO VETOR v2
671  {
672  TPZFNMatrix<2,STATE> res(2,1,0.);
673  for(in = 0 ; in < phi.Rows(); in++)
674  {
675  for (int il=0; il<NumLoadCases(); il++)
676  {
677  ef(nstate*in+0,0) += (bc.Val2(il)(0,0)*data.normal[0]) * phi(in,0) * weight ;
678  ef(nstate*in+1,0) += (bc.Val2(il)(0,0)*data.normal[1]) * phi(in,0) * weight ;
679  }
680  for(jn=0; jn<phi.Rows(); jn++)
681  {
682  for(int idf=0; idf<2; idf++) for(int jdf=0; jdf<2; jdf++)
683  {
684  ek(nstate*in+idf,nstate*jn+jdf) += bc.Val1()(idf,jdf)*phi(in,0)*phi(jn,0)*weight;
685  //BUG FALTA COLOCAR VAL2
686  // DebugStop();
687  }
688  }
689 
690  }
691 
692  }
693  break;
694 
695  } // �nulo introduzindo o BIGNUMBER pelos valores da condi�o
696  } // 1 Val1 : a leitura �00 01 10 11
697 }
698 
699 
700 void TPZElasticityMaterial::ContributeVecShapeBC(TPZMaterialData &data,REAL weight,
701  TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc) {
702 
703  TPZFMatrix<REAL> &phi = data.phi;
704 
705  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
706 
707  int phc = phi.Cols();
708  short in,jn;
709 
710  switch (bc.Type()) {
711  case 0 : // Dirichlet condition
712  for(in = 0 ; in < phc; in++) {
713  for (int il = 0; il <fNumLoadCases; il++)
714  {
715  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
716 
717  ef(in,il) += weight*BIGNUMBER*(v2(0,il)*phi(0,in) + v2(1,il) * phi(1,in));
718  }
719  for (jn = 0 ; jn < phc; jn++) {
720 
721  ek(in,jn) += weight*BIGNUMBER*(phi(0,in)*phi(0,jn) + phi(1,in)*phi(1,jn));
722  }
723  }
724  break;
725 
726  case 1 : // Neumann condition
727  for (in = 0; in < phc; in++)
728  {
729  for (int il = 0; il <fNumLoadCases; il++)
730  {
731  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
732  ef(in,il)+= weight*(v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in));
733  }
734  }
735  break;
736 
737  case 2 : // condicao mista
738  for(in = 0 ; in < phc; in++)
739  {
740  for (int il = 0; il <fNumLoadCases; il++)
741  {
742  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
743  ef(in,il) += weight * (v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in));
744  }
745 
746  for (jn = 0; jn <phc; jn++) {
747 
748  ek(in,jn) += bc.Val1()(0,0)*phi(0,in)*phi(0,jn)*weight
749 
750  + bc.Val1()(1,0)*phi(1,in)*phi(0,jn)*weight
751 
752  + bc.Val1()(0,1)*phi(0,in)*phi(1,jn)*weight
753 
754  + bc.Val1()(1,1)*phi(1,in)*phi(1,jn)*weight;
755  }
756  }// este caso pode reproduzir o caso 0 quando o deslocamento
757  } // eh nulo introduzindo o BIGNUMBER pelos valores da condicao
758 }
759 
760 
762 int TPZElasticityMaterial::VariableIndex(const std::string &name){
763 
764 
765  /*
766  if(!strcmp("Displacement", name.c_str())) return TPZElasticityMaterial::EDisplacement;
767  if(!strcmp("DisplacementX", name.c_str())) return TPZElasticityMaterial::EDisplacementX;
768  if(!strcmp("DisplacementY", name.c_str())) return TPZElasticityMaterial::EDisplacementY;
769  if(!strcmp("DisplacementZ", name.c_str())) return TPZElasticityMaterial::EDisplacementZ;
770  if(!strcmp("NormalStress", name.c_str())) return TPZElasticityMaterial::ENormalStress;
771  if(!strcmp("ShearStress", name.c_str())) return TPZElasticityMaterial::EShearStress;
772  if(!strcmp("NormalStrain", name.c_str())) return TPZElasticityMaterial::ENormalStrain;
773  if(!strcmp("ShearStrain", name.c_str())) return TPZElasticityMaterial::EShearStrain;
774  if(!strcmp("PrincipalStress", name.c_str())) return TPZElasticityMaterial::EPrincipalStress;
775  if(!strcmp("Stress1", name.c_str())) return TPZElasticityMaterial::EStress1;
776  if(!strcmp("PrincipalStrain", name.c_str())) return TPZElasticityMaterial::EPrincipalStrain;
777  if(!strcmp("Strain1", name.c_str())) return TPZElasticityMaterial::EStrain1;
778  if(!strcmp("PrincipalStressDirection1",name.c_str())) return TPZElasticityMaterial::EPrincipalStressDirection1;
779  if(!strcmp("PrincipalStressDirection2",name.c_str())) return TPZElasticityMaterial::EPrincipalStressDirection2;
780  if(!strcmp("PrincipalStressDirection3",name.c_str())) return TPZElasticityMaterial::EPrincipalStressDirection3;
781  if(!strcmp("I1Stress", name.c_str())) return TPZElasticityMaterial::EI1Stress;
782  if(!strcmp("J2Stress", name.c_str())) return TPZElasticityMaterial::EJ2Stress;
783  if(!strcmp("I1J2Stress", name.c_str())) return TPZElasticityMaterial::EI1J2Stress;
784  if(!strcmp("DirStress", name.c_str())) return TPZElasticityMaterial::EDirStress;
785  if(!strcmp("DirStrain", name.c_str())) return TPZElasticityMaterial::EDirStrain;
786  if(!strcmp("VolElasticStrain", name.c_str())) return TPZElasticityMaterial::EVolElasticStrain;
787  if(!strcmp("VolPlasticStrain", name.c_str())) return TPZElasticityMaterial::EVolPlasticStrain;
788  if(!strcmp("VolTotalStrain", name.c_str())) return TPZElasticityMaterial::EVolTotalStrain;
789  if(!strcmp("VolTEPStrain", name.c_str())) return TPZElasticityMaterial::EVolTEPStrain;
790  if(!strcmp("Alpha", name.c_str())) return TPZElasticityMaterial::EAlpha;
791  if(!strcmp("PlasticSteps", name.c_str())) return TPZElasticityMaterial::EPlasticSteps;
792  if(!strcmp("YieldSurface", name.c_str())) return TPZElasticityMaterial::EYield;
793  if(!strcmp("TotalPlasticStrain", name.c_str())) return TPZElasticityMaterial::ENormalPlasticStrain;
794  if(!strcmp("EMisesStress", name.c_str())) return TPZElasticityMaterial::EMisesStress;
795  PZError << "TPZMatElastoPlastic::VariableIndex Error\n";
796  return -1;
797  */
798 
799  if(!strcmp("displacement",name.c_str())) return 9;
800  if(!strcmp("Displacement",name.c_str())) return 9;
801  if(!strcmp("DisplacementMem",name.c_str())) return 9;
802  if(!strcmp("Pressure",name.c_str())) return 1;
803  if(!strcmp("MaxStress",name.c_str())) return 2;
804  if(!strcmp("PrincipalStress1",name.c_str())) return 3;
805  if(!strcmp("PrincipalStress2",name.c_str())) return 4;
806  if(!strcmp("SigmaX",name.c_str())) return 5;
807  if(!strcmp("SigmaY",name.c_str())) return 6;
808  if(!strcmp("TauXY",name.c_str())) return 8;//Cedric
809  if(!strcmp("Strain",name.c_str())) return 11;//Philippe
810  if(!strcmp("SigmaZ",name.c_str())) return 12;//Philippe
811 
812  if(!strcmp("sig_x",name.c_str())) return 5;
813  if(!strcmp("sig_y",name.c_str())) return 6;
814  if(!strcmp("tau_xy",name.c_str())) return 8;//Cedric
815  if(!strcmp("Displacement6",name.c_str())) return 7;
816  if(!strcmp("Stress",name.c_str())) return 10;
817  if(!strcmp("Flux",name.c_str())) return 10;
818  if(!strcmp("J2",name.c_str())) return 20;
819  if(!strcmp("I1",name.c_str())) return 21;
820  if(!strcmp("J2Stress",name.c_str())) return 20;
821  if(!strcmp("I1Stress",name.c_str())) return 21;
822  if(!strcmp("Alpha",name.c_str())) return 22;
823  if(!strcmp("PlasticSqJ2",name.c_str())) return 22;
824  if(!strcmp("PlasticSqJ2El",name.c_str())) return 22;
825  if(!strcmp("YieldSurface",name.c_str())) return 27;
826  if(!strcmp("NormalStress",name.c_str())) return 23;
827  if(!strcmp("ShearStress",name.c_str())) return 24;
828  if(!strcmp("NormalStrain",name.c_str())) return 25;
829  if(!strcmp("ShearStrain",name.c_str())) return 26;
830  if(!strcmp("Young_Modulus",name.c_str())) return 28;
831  if(!strcmp("Poisson",name.c_str())) return 29;
832 
833 
834 
835  // cout << "TPZElasticityMaterial::VariableIndex Error\n";
836  return TPZMaterial::VariableIndex(name);
837 }
838 
840 int TPZElasticityMaterial::NSolutionVariables(int var){
841 
842  switch(var) {
843  case 0:
844  return 2;
845  case 1:
846  case 2:
847  return 1;
848  case 3:
849  case 4:
850  return 2;
851  case 5:
852  case 6:
853  case 8:
854  return 1;
855  case 7:
856  return 6;
857  case 9:
858  return 3;
859  case 10 : //Stress Tensor
860  return 3;
861  case 11 : //Strain Tensor
862  return 3;
863  // SigZ
864  case 12:
865  return 1;
866  case 20:
867  return 1;
868  case 21:
869  return 1;
870  case 22:
871  return 1;
872  case 23:
873  case 24:
874  case 25:
875  case 26:
876  case 27:
877  return 3;
878  case 28:
879  case 29:
880  return 1;
881  default:
882  return TPZMaterial::NSolutionVariables(var);
883  }
884 }
885 
886 void TPZElasticityMaterial::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout)
887 {
888  int numbersol = data.dsol.size();
889  int ipos = 0;
890  if (fPostProcIndex < numbersol) {
891  ipos = fPostProcIndex;
892  }
893 
894  REAL E(fE_def), nu(fnu_def);
895 
896  if (fElasticity) {
897  TPZManVector<STATE, 2> result(2);
898  TPZFNMatrix<4, STATE> Dres(0, 0);
899  fElasticity->Execute(data.x, result, Dres);
900  E = result[0];
901  nu = result[1];
902  }
903 
904  if(var == 28)
905  {
906  Solout[0] = E;
907  return ;
908  }
909  if(var == 29)
910  {
911  Solout[0] = nu;
912  return;
913  }
914 
915 
916  REAL Eover1MinNu2 = E/(1-nu*nu);
917  REAL Eover21PlusNu = E/(2.*(1+nu));
918 
919 
920  TPZVec<STATE> &Sol = data.sol[ipos];
921  TPZFMatrix<STATE> &DSol = data.dsol[ipos];
922  TPZFMatrix<REAL> &axes = data.axes;
923  TPZFNMatrix<4,STATE> DSolxy(2,2);
924 
925  REAL epsx;
926  REAL epsy;
927  REAL epsxy;
928  REAL epsz = 0.;
929  REAL SigX;
930  REAL SigY;
931  REAL SigZ;
932  REAL TauXY,aux,Sig1,Sig2,angle;
933 
934  // dudx - dudy
935  DSolxy(0,0) = DSol(0,0)*axes(0,0)+DSol(1,0)*axes(1,0);
936  DSolxy(1,0) = DSol(0,0)*axes(0,1)+DSol(1,0)*axes(1,1);
937  // dvdx - dvdy
938  DSolxy(0,1) = DSol(0,1)*axes(0,0)+DSol(1,1)*axes(1,0);
939  DSolxy(1,1) = DSol(0,1)*axes(0,1)+DSol(1,1)*axes(1,1);
940 
941  epsx = DSolxy(0,0);// du/dx
942  epsy = DSolxy(1,1);// dv/dy
943  epsxy = 0.5*(DSolxy(1,0)+DSolxy(0,1));
944 
945  REAL lambda = GetLambda(E,nu);
946  REAL mu = GetMU(E,nu);
947  if (this->fPlaneStress == 1) {
948  epsz = -lambda*(epsx+epsy)/(lambda+2.*mu);
949  }
950  else {
951  epsz = 0.;
952  }
953  TauXY = 2*mu*epsxy+fPreStressXY;
954 #ifdef PZDEBUG
955  REAL tol = 0.;
956  ZeroTolerance(tol);
957  tol *= 100;
958  if(fabs(TauXY) > 1.) tol *= fabs(TauXY);
959  REAL TauXY2 = E*epsxy/(1.+nu)+fPreStressXY;
960 #ifdef REALfloat
961  if (fabs(TauXY-TauXY2) > tol) {
962  DebugStop();
963  }
964 #else
965  if (fabs(TauXY-TauXY2) > tol) {
966  DebugStop();
967  }
968 #endif
969 #endif
970  if (this->fPlaneStress == 1){
971  SigX = Eover1MinNu2*(epsx+nu*epsy)+fPreStressXX;
972  SigY = Eover1MinNu2*(nu*epsx+epsy)+fPreStressYY;
973  SigZ = fPreStressZZ;
974  }
975  else
976  {
977  SigX = E/((1.-2.*nu)*(1.+nu))*((1.-nu)*epsx+nu*epsy)+fPreStressXX;
978  SigY = E/((1.-2.*nu)*(1.+nu))*(nu*epsx+(1.-nu)*epsy)+fPreStressYY;
979  SigZ = fPreStressZZ+lambda*(epsx+epsy);
980  }
981 
982  switch(var) {
983  case 0:
984  //numvar = 2;
985  Solout[0] = Sol[0];
986  Solout[1] = Sol[1];
987  break;
988  case 7:
989  //numvar = 6;
990  Solout[0] = Sol[0];
991  Solout[1] = Sol[1];
992  Solout[2] = 0.;
993  Solout[3] = 0.;
994  Solout[4] = 0.;
995  Solout[5] = 0.;
996  break;
997  case 1:
998  case 2:
999  case 3:
1000  case 4:
1001  case 5:
1002  case 6:
1003  case 8:
1004  case 10:
1005 
1006  //numvar = 1;
1007  Solout[0] = SigX+SigY+SigZ;
1008  // Pressure variable
1009  if(var == 1) {
1010  Solout[0] = SigX+SigY+SigZ;
1011  return;
1012  }
1013  // TauXY variable
1014  if(var == 8) {
1015  Solout[0] = TauXY;
1016  return;
1017  }
1018  if(var ==5) {
1019  Solout[0] = SigX;
1020  return;
1021  }
1022  if(var == 6) {
1023  Solout[0] = SigY;
1024  return;
1025  }
1026  aux = sqrt(0.25*(SigX-SigY)*(SigX-SigY)
1027  +(TauXY)*(TauXY));
1028  // Philippe 13/5/99
1029  // if(abs(Tau) < 1.e-10 && abs(SigY-SigX) < 1.e-10) angle = 0.;
1030  if(fabs(TauXY) < 1.e-10 && fabs(SigY-SigX) < 1.e-10) angle = 0.;
1031  else angle = atan2(2*TauXY,SigY-SigX)/2.;
1032  Sig1 = 0.5*(SigX+SigY)+aux;
1033  Sig2 = 0.5*(SigX+SigY)-aux;
1034  if(var == 3 ){
1035  //numvar = 2;
1036  Solout[0] = Sig1*cos(angle);
1037  Solout[1] = Sig1*sin(angle);
1038  return;
1039  }
1040  if(var == 4 ) {
1041  //numvar = 2;
1042  Solout[0] = -Sig2*sin(angle);
1043  Solout[1] = Sig2*cos(angle);
1044  return;
1045  }
1046  if(var == 2) {
1047  REAL sigmax;
1048  sigmax = (fabs(Sig1) < fabs(Sig2))? fabs(Sig2) : fabs(Sig1);
1049  Solout[0] = sigmax;
1050  return;
1051  }
1052  if (var ==10)
1053  {
1054  Solout[0] = SigX;
1055  Solout[1] = SigY;
1056  Solout[2] = TauXY;
1057  return;
1058  }
1059  cout << "Very critical error TPZElasticityMaterial::Solution\n";
1060  exit(-1);
1061  // Solout[0] /= 0.;
1062  break;
1063  case 9:
1064  Solout[0] = Sol[0];
1065  Solout[1] = Sol[1];
1066  Solout[2] = 0.;
1067  break;
1068  case 11:
1069  Solout[0] = epsx;
1070  Solout[1] = epsy;
1071  Solout[2] = epsxy;
1072  break;
1073  case 12:
1074  Solout[0] = SigZ;
1075  break;
1076 
1077  case 20:
1078  {
1079 
1080  REAL J2 = (pow(SigX + SigY,2) - (3*(-pow(SigX,2) - pow(SigY,2) + pow(SigX + SigY,2) - 2*pow(TauXY,2)))/2.)/2.;
1081 
1082  Solout[0]=J2;
1083  break;
1084  }
1085  case 21:
1086  {
1087  REAL I1 = SigX+SigY;
1088  Solout[0]=I1;
1089  break;
1090  }
1091  case 22:
1092  Solout[0] = 0.;
1093  break;
1094  case 23:
1095  // normal stress
1096  Solout[0] = SigX;
1097  Solout[1] = SigY;
1098  Solout[2] = SigZ;
1099  break;
1100  case 24:
1101  // shear stress
1102  Solout[0] = TauXY;
1103  Solout[1] = 0.;
1104  Solout[2] = 0.;
1105  break;
1106  case 25:
1107  Solout[0] = epsx;
1108  Solout[1] = epsy;
1109  Solout[2] = epsz;
1110  break;
1111  case 26:
1112  Solout[0] = epsxy;
1113  Solout[1] = 0.;
1114  Solout[2] = 0.;
1115  break;
1116  case 27:
1117  Solout[0] = 0.;
1118  Solout[1] = 0.;
1119  Solout[2] = 0.;
1120  break;
1121  default:
1122  TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
1123  break;
1124  }
1125 }
1126 
1127 
1128 void TPZElasticityMaterial::Flux(TPZVec<REAL> &x, TPZVec<STATE> &Sol, TPZFMatrix<STATE> &DSol, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux) {
1129  if(fabs(axes(0,2)) >= 1.e-6 || fabs(axes(1,2)) >= 1.e-6) {
1130  cout << "TPZElasticityMaterial::Flux only serves for xy configuration\n";
1131  axes.Print("axes");
1132  }
1133 }
1134 
1135 void TPZElasticityMaterial::Errors(TPZVec<REAL> &x,TPZVec<STATE> &u,
1136  TPZFMatrix<STATE> &dudx, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux,
1137  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) {
1138  values[0] = 0.;
1139  TPZManVector<REAL,3> sigma(3,0.),sigma_exact(3,0.);
1140  REAL sigx,sigy,sigxy;
1141  TPZFNMatrix<6,STATE> du(3,dudx.Cols());
1142  TPZAxesTools<STATE>::Axes2XYZ(dudx,du,axes);
1143 // du(0,0) = dudx(0,0)*axes(0,0)+dudx(1,0)*axes(1,0);
1144 // du(1,0) = dudx(0,0)*axes(0,1)+dudx(1,0)*axes(1,1);
1145 // du(0,1) = dudx(0,1)*axes(0,0)+dudx(1,1)*axes(1,0);
1146 // du(1,1) = dudx(0,1)*axes(0,1)+dudx(1,1)*axes(1,1);
1147 
1148  REAL E(fE_def), nu(fnu_def);
1149 
1150  if (fElasticity) {
1151  TPZManVector<STATE,2> result(2);
1152  TPZFNMatrix<4,STATE> Dres(0,0);
1153  fElasticity->Execute(x, result, Dres);
1154  E = result[0];
1155  nu = result[1];
1156  }
1157 
1158  REAL Eover1MinNu2 = E/(1-nu*nu);
1159  REAL Eover21PlusNu = E/(2.*(1+nu));
1160 
1161  STATE epsx,epsy,epsxy, epsz;
1162  epsx = du(0,0);// du/dx
1163  epsy = du(1,1);// dv/dy
1164  epsxy = 0.5*(du(1,0)+du(0,1));
1165  REAL lambda = GetLambda(E,nu);
1166  REAL mu = GetMU(E,nu);
1167  if (fPlaneStress) {
1168  epsz = -lambda*(epsx+epsy)/(lambda+2.*mu);
1169  }
1170  else
1171  {
1172  epsz = 0.;
1173  }
1174  // epsz = data[1].sol[0][0];
1175  STATE SigX = lambda*(epsx+epsy+epsz)+2.*mu*epsx + fPreStressXX;
1176  STATE SigY = lambda*(epsx+epsy+epsz)+2.*mu*epsy + fPreStressYY;
1177 
1178  STATE SigZ = lambda*(epsx+epsy+epsz)+2.*mu*epsz;
1179  STATE TauXY = 2*mu*epsxy+fPreStressXY;
1180 
1181  //tensoes aproximadas : uma forma
1182  sigma[0] = SigX;
1183  sigma[1] = SigY;
1184  sigma[2] = TauXY;
1185 
1186  //exata
1187  STATE epsx_exact,epsy_exact,epsxy_exact, epsz_exact;
1188  epsx_exact = du_exact(0,0);// du/dx
1189  epsy_exact = du_exact(1,1);// dv/dy
1190  epsxy_exact = 0.5*(du_exact(1,0)+du_exact(0,1));
1191  if (fPlaneStress) {
1192  epsz_exact = -lambda*(epsx+epsy)/(lambda+2.*mu);
1193  }
1194  else
1195  {
1196  epsz_exact = 0.;
1197  }
1198  // epsz = data[1].sol[0][0];
1199  SigX = lambda*(epsx_exact+epsy_exact+epsz_exact)+2.*mu*epsx_exact + fPreStressXX;
1200  SigY = lambda*(epsx_exact+epsy_exact+epsz_exact)+2.*mu*epsy_exact + fPreStressYY;
1201 
1202  SigZ = lambda*(epsx_exact+epsy_exact+epsz_exact)+2.*mu*epsz_exact;
1203  TauXY = 2*mu*epsxy_exact+fPreStressXY;
1204 
1205  sigma_exact[0] = SigX;
1206  sigma_exact[1] = SigY;
1207  sigma_exact[2] = TauXY;
1208  sigx = (sigma[0] - sigma_exact[0]);
1209  sigy = (sigma[1] - sigma_exact[1]);
1210  sigxy = (sigma[2] - sigma_exact[2]);
1211 
1212  //values[0] = calculo do erro estimado em norma Energia
1213  values[0] = (sigx*(epsx-epsx_exact)+sigy*(epsy-epsy_exact)+2.*sigxy*(epsxy-epsxy_exact));
1214 
1215  //values[3] = calculo da energia da solucao exata
1216  values[3] = (SigX*(epsx_exact)+SigY*(epsy_exact)+2.*TauXY*(epsxy_exact));
1217 
1218  //values[4] : erro em norma L2 em tensoes
1219  values[4] = sigx*sigx + sigy*sigy + 2.*sigxy*sigxy;
1220 
1221  //values[5] : erro em norma L2 em sig_xx
1222  values[5] = sigx*sigx;
1223 
1224  //values[1] : erro em norma L2 em deslocamentos
1225  values[1] = (u[0] - u_exact[0])*(u[0] - u_exact[0])+(u[1] - u_exact[1])*(u[1] - u_exact[1]);
1226 
1227  //values[2] : erro estimado na norma H1
1228  REAL SemiH1 =0.;
1229  for(int i = 0; i < 2; i++) for(int j = 0; j < 2; j++) SemiH1 += (du(i,j) - du_exact(i,j)) * (du(i,j) - du_exact(i,j));
1230  values[2] = values[1] + SemiH1;
1231 }
1232 
1233 
1234 TPZElasticityMaterial::TPZElasticityMaterial(const TPZElasticityMaterial &copy) :
1235 TPZRegisterClassId(&TPZElasticityMaterial::ClassId),
1236 TPZDiscontinuousGalerkin(copy),
1237 fE_def(copy.fE_def),
1238 fnu_def(copy.fnu_def),
1239 fElasticity(copy.fElasticity),
1240 fPreStressXX(copy.fPreStressXX),
1241 fPreStressYY(copy.fPreStressYY),
1242 fPreStressXY(copy.fPreStressXY),
1243 fPreStressZZ(copy.fPreStressZZ)
1244 {
1245  ff = copy.ff;
1246  fPlaneStress = copy.fPlaneStress;
1247  // Added by Philippe 2012
1248  fPostProcIndex = copy.fPostProcIndex;
1249 
1250 }
1251 
1252 
1253 int TPZElasticityMaterial::ClassId() const{
1254  return Hash("TPZElasticityMaterial") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
1255 }
1256 
1257 #ifndef BORLAND
1258 template class TPZRestoreClass<TPZElasticityMaterial>;
1259 #endif
1260 
1261 void TPZElasticityMaterial::Read(TPZStream &buf, void *context)
1262 {
1263  TPZMaterial::Read(buf,context);
1264  buf.Read(&fE_def,1);
1265  buf.Read(&fnu_def,1);
1266  buf.Read(&fPreStressXX,1);
1267  buf.Read(&fPreStressYY,1);
1268  buf.Read(&fPreStressXY,1);
1269  buf.Read(&fPreStressZZ,1);
1270 
1271  buf.Read(&ff[0],3);
1272  buf.Read(&fPlaneStress,1);
1273  buf.Read(&fPostProcIndex);
1274 
1275 }
1276 
1277 void TPZElasticityMaterial::Write(TPZStream &buf, int withclassid) const
1278 {
1279  TPZMaterial::Write(buf,withclassid);
1280  buf.Write(&fE_def,1);
1281  buf.Write(&fnu_def,1);
1282  buf.Write(&fPreStressXX,1);
1283  buf.Write(&fPreStressYY,1);
1284  buf.Write(&fPreStressXY,1);
1285  buf.Write(&fPreStressZZ,1);
1286 
1287  buf.Write(&ff[0],3);
1288  buf.Write(&fPlaneStress,1);
1289  buf.Write(&fPostProcIndex);
1290 
1291 }
1292 
TPZElasticityMaterial::Print
virtual void Print(std::ostream &out=std::cout) override
Print the material data.
Definition: pzelasmat.cpp:92
TPZDiscontinuousGalerkin
Defines the interface which material objects need to implement for discontinuous Galerkin formulation...
Definition: pzdiscgal.h:20
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
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
TPZElasticityMaterial::fPlaneStress
int fPlaneStress
Uses plain stress.
Definition: pzelasmat.h:270
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
TPZElasticityMaterial::GetMU
STATE GetMU(REAL E, REAL nu) const
Definition: pzelasmat.h:186
TPZDiscontinuousGalerkin::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzdiscgal.cpp:110
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
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
TPZElasticityMaterial::fPreStressYY
REAL fPreStressYY
Pre Stress Tensor - Sigma YY.
Definition: pzelasmat.h:261
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZElasticityMaterial::Flux
virtual void Flux(TPZVec< REAL > &x, TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux) override
Computes the value of the flux function to be used by ZZ error estimator.
Definition: pzelasmat.cpp:1128
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
TPZElasticityMaterial::Solution
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout) override
Returns the solution associated with the var index based on the finite element approximation.
Definition: pzelasmat.cpp:886
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
TPZElasticityMaterial::fPreStressZZ
REAL fPreStressZZ
Pre Stress Tensor - Sigma ZZ.
Definition: pzelasmat.h:267
TPZElasticityMaterial::TPZElasticityMaterial
TPZElasticityMaterial()
Default constructor.
Definition: pzelasmat.cpp:24
TPZElasticityMaterial::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
Computes the error due to the difference between the interpolated flux and the flux computed based o...
Definition: pzelasmat.cpp:1135
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.
TPZElasticityMaterial::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the.
Definition: pzelasmat.cpp:417
TPZElasticityMaterial::fPreStressXY
REAL fPreStressXY
Pre Stress Tensor - Sigma XY.
Definition: pzelasmat.h:264
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
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
TPZElasticityMaterial::SetPreStress
void SetPreStress(REAL Sigxx, REAL Sigyy, REAL Sigxy, REAL Sigzz)
Set PresStress Tensor.
Definition: pzelasmat.cpp:112
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
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZElasticityMaterial::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: pzelasmat.cpp:840
TPZMaterial::fPostProcIndex
int fPostProcIndex
indicates which solution should be used for post processing
Definition: TPZMaterial.h:79
pzelmat.h
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZElasticityMaterial::ContributeVecShapeBC
void ContributeVecShapeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: pzelasmat.cpp:700
TPZBndCond::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: pzbndcond.cpp:142
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
sin
sin
Definition: fadfunc.h:63
TPZElasticityMaterial::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzelasmat.cpp:1253
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
TPZElasticityMaterial::Read
virtual void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzelasmat.cpp:1261
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
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
pzelasmat.h
Contains the TPZElasticityMaterial class which implements a two dimensional elastic material in plane...
TPZElasticityMaterial::ContributeVecShape
void ContributeVecShape(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Calculates the element stiffness matrix - simulate compaction as aditional variable.
Definition: pzelasmat.cpp:433
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
TPZElasticityMaterial::fElasticity
TPZAutoPointer< TPZFunction< STATE > > fElasticity
Definition: pzelasmat.h:246
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
TPZElasticityMaterial::E
REAL E()
Returns the elasticity modulus E.
Definition: pzelasmat.h:220
TPZElasticityMaterial::ff
TPZManVector< STATE, 3 > ff
Forcing vector.
Definition: pzelasmat.h:249
TPZElasticityMaterial::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzelasmat.cpp:1277
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
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZElasticityMaterial::Name
std::string Name() override
Returns the material name.
Definition: pzelasmat.h:90
angle
static REAL angle
Angle in radians to test.
Definition: pzsubcmesh.cpp:53
TPZElasticityMaterial::GetLambda
STATE GetLambda(REAL E, REAL nu) const
Definition: pzelasmat.h:180
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZElasticityMaterial::FillBoundaryConditionDataRequirement
virtual void FillBoundaryConditionDataRequirement(int type, TPZMaterialData &data) override
This method defines which parameters need to be initialized in order to compute the contribution of t...
Definition: pzelasmat.cpp:424
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
TPZElasticityMaterial::~TPZElasticityMaterial
virtual ~TPZElasticityMaterial()
Default destructor.
Definition: pzelasmat.cpp:85
TPZElasticityMaterial::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Calculates the element stiffness matrix.
Definition: pzelasmat.cpp:120
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZElasticityMaterial
This class implements a two dimensional elastic material in plane stress or strain.
Definition: pzelasmat.h:19
TPZMaterial::NumLoadCases
int NumLoadCases()
returns the number of load cases for this material object
Definition: TPZMaterial.h:186
TPZElasticityMaterial::Dimension
int Dimension() const override
Returns the model dimension.
Definition: pzelasmat.h:81
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
TPZMatrix::Print
virtual void Print(std::ostream &out) const
Definition: pzmatrix.h:253
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZElasticityMaterial::fPreStressXX
REAL fPreStressXX
Pre Stress Tensor - Sigma XX.
Definition: pzelasmat.h:258
rdt.values
def values
Definition: rdt.py:119
TPZElasticityMaterial::fnu_def
REAL fnu_def
Poison coeficient.
Definition: pzelasmat.h:243
ZeroTolerance
REAL ZeroTolerance()
Returns the tolerance to Zero value. Actually: .
Definition: pzreal.h:633
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
TPZElasticityMaterial::fE_def
REAL fE_def
Elasticity modulus.
Definition: pzelasmat.h:240
TPZElasticityMaterial::ContributeBC
virtual void ContributeBC(TPZVec< TPZMaterialData > &datavec, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
It computes a contribution to the stiffness matrix and load vector at one BC integration point to mul...
Definition: pzelasmat.cpp:524
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
TPZElasticityMaterial::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: pzelasmat.cpp:88
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
TPZElasticityMaterial::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns the variable index associated with the name.
Definition: pzelasmat.cpp:762
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13