NeoPZ
TPZMatElasticity2D.cpp
Go to the documentation of this file.
1 //
2 // TPZMatElasticity2D.cpp
3 // PZ
4 //
5 // Created by Omar on 10/27/14.
6 //
7 //
8 
9 
10 #include <iostream>
11 #include <string>
12 #include "TPZMatElasticity2D.h"
13 #include "pzbndcond.h"
14 #include "pzaxestools.h"
15 #include "pzlog.h"
16 
17 #ifdef LOG4CXX
18 static LoggerPtr logger(Logger::getLogger("pz.elasticity"));
19 #endif
20 
21 
22 TPZMatElasticity2D::TPZMatElasticity2D()
23 : TPZRegisterClassId(&TPZMatElasticity2D::ClassId), TPZMaterial()
24 {
25  m_E = 0.;
26  m_nu = 0.;
27  m_lambda = 0.;
28  m_mu = 0.;
29  m_f.resize(2);
30  m_f[0]=0.;
31  m_f[1]=0.;
32  m_plane_stress = 1;
33  m_s0_xx = 0.0;
34  m_s0_xy = 0.0;
35  m_s0_yy = 0.0;
36  m_s0_zz = 0.0;
37 
38 }
39 
40 TPZMatElasticity2D::TPZMatElasticity2D(int matid)
41 : TPZRegisterClassId(&TPZMatElasticity2D::ClassId), TPZMaterial(matid)
42 {
43  m_E = 0.;
44  m_nu = 0.;
45  m_lambda = 0.;
46  m_mu = 0.;
47  m_f.resize(2);
48  m_f[0]=0.;
49  m_f[1]=0.;
50  m_plane_stress = 1;
51  m_s0_xx = 0.0;
52  m_s0_xy = 0.0;
53  m_s0_yy = 0.0;
54  m_s0_zz = 0.0;
55 }
56 
57 TPZMatElasticity2D::TPZMatElasticity2D(int matid, REAL E, REAL nu, REAL fx, REAL fy, int plainstress)
58 : TPZRegisterClassId(&TPZMatElasticity2D::ClassId), TPZMaterial(matid)
59 {
60  m_E = E;
61  m_nu = nu;
62  m_lambda = (E*nu)/((1+nu)*(1-2*nu));
63  m_mu = E/(2*(1+nu));
64  m_f.resize(2);
65  m_f[0]=fx;
66  m_f[1]=fy;
67  m_plane_stress = plainstress;
68  m_s0_xx = 0.0;
69  m_s0_xy = 0.0;
70  m_s0_yy = 0.0;
71  m_s0_zz = 0.0;
72 }
73 
74 TPZMatElasticity2D::~TPZMatElasticity2D()
75 {
76 }
77 
78 
79 TPZMatElasticity2D::TPZMatElasticity2D(const TPZMatElasticity2D &copy)
80 : TPZRegisterClassId(&TPZMatElasticity2D::ClassId), TPZMaterial(copy)
81 {
82  m_E = copy.m_E;
83  m_nu = copy.m_nu;
84  m_lambda = copy.m_lambda;
85  m_mu = copy.m_mu;
86  m_f.resize(copy.m_f.size());
87  for (int i = 0; i < copy.m_f.size(); i++) {
88  m_f[i] = copy.m_f[i];
89  }
90  m_plane_stress = copy.m_plane_stress;
91  m_s0_xx = copy.m_s0_xx;
92  m_s0_xy = copy.m_s0_xy;
93  m_s0_yy = copy.m_s0_yy;
94  m_s0_zz = copy.m_s0_zz;
95 }
96 
97 TPZMatElasticity2D & TPZMatElasticity2D::operator=(const TPZMatElasticity2D &copy)
98 {
99  TPZMaterial::operator = (copy);
100  m_E = copy.m_E;
101  m_nu = copy.m_nu;
102  m_lambda = copy.m_lambda;
103  m_mu = copy.m_mu;
104  m_s0_xx = copy.m_s0_xx;
105  m_s0_xy = copy.m_s0_xy;
106  m_s0_yy = copy.m_s0_yy;
107  m_s0_zz = copy.m_s0_zz;
108  m_f.resize(copy.m_f.size());
109  for (int i = 0; i < copy.m_f.size(); i++) {
110  m_f[i] = copy.m_f[i];
111  }
112  m_plane_stress = copy.m_plane_stress;
113  return *this;
114 }
115 
116 int TPZMatElasticity2D::NStateVariables() const {
117  return 2;
118 }
119 
120 
121 void TPZMatElasticity2D::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef) {
122 
123 
124  if (data.fShapeType == TPZMaterialData::EVecShape) {
125  ContributeVec(data,weight,ek,ef);
126  return;
127  }
128  // Getting weight functions
129  TPZFMatrix<REAL> &phi_u = data.phi;
130  TPZFMatrix<REAL> &dphi_u = data.dphix;
131  int n_phi_u = phi_u.Rows();
132  int first_u = 0;
133 
134  TPZFNMatrix<40,REAL> grad_phi_u(3,n_phi_u);
135  TPZAxesTools<REAL>::Axes2XYZ(data.dphix, grad_phi_u, data.axes);
136  TPZFMatrix<STATE> dsol_u = data.dsol[0];
137 
138  TPZManVector<STATE,3> grad_p(m_f);
139 
140  if(this->HasForcingFunction())
141  {
142  fForcingFunction->Execute(data.x,grad_p);
143  }
144 
145  REAL dvdx,dvdy,dudx,dudy;
146  REAL duxdx,duxdy,duydx,duydy;
147 
148  // Gradient for ux
149  duxdx = dsol_u(0,0)*data.axes(0,0)+dsol_u(1,0)*data.axes(1,0); // dux/dx
150  duxdy = dsol_u(0,0)*data.axes(0,1)+dsol_u(1,0)*data.axes(1,1); // dux/dy
151 
152  // Gradient for uy
153  duydx = dsol_u(0,1)*data.axes(0,0)+dsol_u(1,1)*data.axes(1,0); // duy/dx
154  duydy = dsol_u(0,1)*data.axes(0,1)+dsol_u(1,1)*data.axes(1,1); // duy/dy
155 
156  for(int iu = 0; iu < n_phi_u; iu++ )
157  {
158  dvdx = grad_phi_u(0,iu);
159  dvdy = grad_phi_u(1,iu);
160 
161  ef(2*iu + first_u) += weight * (grad_p[0] * phi_u(iu, 0) - (dvdx*m_s0_xx + dvdy*m_s0_xy) ); // x direction
162  ef(2*iu+1 + first_u) += weight * (grad_p[1] * phi_u(iu, 0) - (dvdx*m_s0_xy + dvdy*m_s0_yy) ); // y direction
163 
164  if (m_plane_stress == 1)
165  {
166  /* Plain stress state */
167  ef(2*iu + first_u) += weight*((4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dvdx*duxdx + (2*m_mu)*dvdy*duxdy);
168 
169  ef(2*iu + first_u) += weight*((2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dvdx*duydy + (2*m_mu)*dvdy*duydx);
170 
171  ef(2*iu+1 + first_u) += weight*((2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dvdy*duxdx + (2*m_mu)*dvdx*duxdy);
172 
173  ef(2*iu+1 + first_u) += weight*((4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dvdy*duydy + (2*m_mu)*dvdx*duydx);
174  }
175  else
176  {
177  /* Plain Strain State */
178  ef(2*iu + first_u) += weight* ((m_lambda + 2*m_mu)*dvdx*duxdx + (m_mu)*dvdy*(duxdy));
179 
180  ef(2*iu + first_u) += weight* (m_lambda*dvdx*duydy + (m_mu)*dvdy*(duydx));
181 
182  ef(2*iu+1 + first_u) += weight* (m_lambda*dvdy*duxdx + (m_mu)*dvdx*(duxdy));
183 
184  ef(2*iu+1 + first_u) += weight* ((m_lambda + 2*m_mu)*dvdy*duydy + (m_mu)*dvdx*(duydx));
185  }
186 
187  for(int ju = 0; ju < n_phi_u; ju++)
188  {
189 
190  dudx = grad_phi_u(0,ju);
191  dudy = grad_phi_u(1,ju);
192 
193  if (this->m_plane_stress == 1)
194  {
195  /* Plain stress state */
196  ek(2*iu + first_u, 2*ju + first_u) += weight*((4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dvdx*dudx + (m_mu)*dvdy*dudy);
197 
198  ek(2*iu + first_u, 2*ju+1 + first_u) += weight*((2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dvdx*dudy + (m_mu)*dvdy*dudx);
199 
200  ek(2*iu+1 + first_u, 2*ju + first_u) += weight*((2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dvdy*dudx + (m_mu)*dvdx*dudy);
201 
202  ek(2*iu+1 + first_u, 2*ju+1 + first_u) += weight*((4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dvdy*dudy + (m_mu)*dvdx*dudx);
203  }
204  else
205  {
206  /* Plain Strain State */
207  ek(2*iu + first_u,2*ju + first_u) += weight* ((m_lambda + 2*m_mu)*dvdx*dudx + (m_mu)*dvdy*dudy);
208 
209  ek(2*iu + first_u,2*ju+1 + first_u) += weight* (m_lambda*dvdx*dudy + (m_mu)*dvdy*dudx);
210 
211  ek(2*iu+1 + first_u,2*ju + first_u) += weight* (m_lambda*dvdy*dudx + (m_mu)*dvdx*dudy);
212 
213  ek(2*iu+1 + first_u,2*ju+1 + first_u) += weight* ((m_lambda + 2*m_mu)*dvdy*dudy + (m_mu)*dvdx*dudx);
214 
215  }
216  }
217  }
218 
219 }
220 
221 void TPZMatElasticity2D::ContributeVec(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef)
222 {
223  if (data.fShapeType != TPZMaterialData::EVecShape) {
224  DebugStop();
225  }
226 
227  // Getting weight functions
228  TPZFMatrix<REAL> &dphiU = data.dphix;
229  int phrU = dphiU.Cols();
230 
231  TPZFNMatrix<200,REAL> dudaxes(2,dphiU.Cols()), dvdaxes(2,dphiU.Cols()), dudx(3,phrU), dvdx(3,phrU);
232  for (int i=0; i<2; i++) {
233  for (int j=0; j<phrU; j++) {
234  dudaxes(i,j) = dphiU(i,j);
235  dvdaxes(i,j) = dphiU(2+i,j);
236  }
237  }
238  TPZAxesTools<REAL>::Axes2XYZ(dudaxes, dudx, data.axes);
239  TPZAxesTools<REAL>::Axes2XYZ(dvdaxes, dvdx, data.axes);
240 
241 
242  for(int iu = 0; iu < phrU; iu++ )
243  {
244  TPZFNMatrix<4,REAL> gradv(2,2);
245  gradv(0,0) = dudx(0,iu);
246  gradv(0,1) = dudx(1,iu);
247  gradv(1,0) = dvdx(0,iu);
248  gradv(1,1) = dvdx(1,iu);
249 
250  for(int ju = 0; ju < phrU; ju++)
251  {
252  TPZFNMatrix<4,REAL> gradu(2,2);
253  gradu(0,0) = dudx(0,ju);
254  gradu(0,1) = dudx(1,ju);
255  gradu(1,0) = dvdx(0,ju);
256  gradu(1,1) = dvdx(1,ju);
257 
258  if (this->m_plane_stress == 1)
259  {
260  /* Plain stress state
261  \sigma_x = E/(1-\nu\nu) (\epsilon_x + \nu \epsilon_y)
262  \sigma_x = \frac{4\mu(\lambda+\mu)}{\lambda+2\mu)}\epsilon_x + \frac{2\mu\lambda}{\lambda+2\mu} \epsilon_y
263  \sigma_y = E/(1-\nu\nu) (\epsilon_y + \nu \epsilon_x)
264  \sigma_y = \frac{4\mu(\lambda+\mu)}{\lambda+2\mu)}\epsilon_y + \frac{2\mu\lambda}{\lambda+2\mu} \epsilon_x
265  \tau_{xy} = \frac{E}{1+\nu} \epsilon_{xy}
266  \tau_{xy} = \frac{1}{2\mu} \epsilon_{xy}
267  */
268  TPZFNMatrix<4,REAL> sigma_u(2,2);
269  sigma_u(0,0) = m_E/(1-m_nu*m_nu) *(gradu(0,0)+m_nu*gradu(1,1));
270  sigma_u(1,1) = m_E/(1-m_nu*m_nu) *(gradu(1,1)+m_nu*gradu(0,0));
271  sigma_u(0,1) = m_E/(2.*(1+m_nu))*(gradu(0,1)+gradu(1,0));
272  sigma_u(1,0) = sigma_u(0,1);
273  ek(iu, ju) += weight*(sigma_u(0,0)*gradv(0,0)+sigma_u(1,1)*gradv(1,1)+
274  sigma_u(1,0)*gradv(1,0)+sigma_u(0,1)*gradv(0,1));
275  }
276  else
277  {
278  /* Plain Strain State */
279 
280  /*
281  \sigma_x = \frac{E}{(1+\nu)(1-2\nu)}((1-\nu)\epsilon_x + \nu\epsilon_y)
282  \sigma_x = (\lambda+2\mu)\epsilon_x + \lambda\epsilon_y
283  \sigma_y = \frac{E}{(1+\nu)(1-2\nu)}((1-\nu)\epsilon_y + \nu\epsilon_x)
284  \sigma_y = (\lambda+2\mu)\epsilon_y + \lambda\epsilon_x
285  \tau_{xy} = \frac{E}{1+\nu} \epsilon_{xy}
286  \tau_{xy} = \frac{1}{2\mu} \epsilon_{xy}
287  */
288  TPZFNMatrix<4,REAL> sigma_u(2,2);
289  sigma_u(0,0) = (m_lambda+2.*m_mu)*gradu(0,0)+m_lambda*gradu(1,1);
290  sigma_u(1,1) = (m_lambda+2.*m_mu)*gradu(1,1)+m_lambda*gradu(0,0);
291  sigma_u(0,1) = m_mu*(gradu(0,1)+gradu(1,0));
292  sigma_u(1,0) = sigma_u(0,1);
293  STATE energy = (sigma_u(0,0)*gradv(0,0)+sigma_u(1,1)*gradv(1,1)+
294  sigma_u(1,0)*gradv(1,0)+sigma_u(0,1)*gradv(0,1));
295  ek(iu, ju) += weight*energy;
296 
297  }
298  }
299  }
300  // ////////////////////////// Jacobian Matrix ///////////////////////////////////
301  this->ContributeVec(data,weight,ef);
302 
303 }
304 void TPZMatElasticity2D::ContributeVec(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ef)
305 {
306  // Getting weight functions
307  TPZFMatrix<REAL> &phi_u = data.phi;
308  int n_phi_u = phi_u.Rows()/2;
309  TPZFMatrix<REAL> &dphiU = data.dphix;
310 
311  TPZFNMatrix<200,REAL> dudaxes(2,dphiU.Cols()), dvdaxes(2,dphiU.Cols()), dudx(2,n_phi_u), dvdx(2,n_phi_u);
312  for (int i=0; i<2; i++) {
313  for (int j=0; j<n_phi_u; j++) {
314  dudaxes(i,j) = dphiU(i,j);
315  dvdaxes(i,j) = dphiU(2+i,j);
316  }
317  }
318  TPZAxesTools<REAL>::Axes2XYZ(dudaxes, dudx, data.axes);
319  TPZAxesTools<REAL>::Axes2XYZ(dvdaxes, dvdx, data.axes);
320 
321  TPZManVector<STATE,3> sol_u =data.sol[0];
322  TPZFNMatrix<4,STATE> dsol_xy(2,2), dsol_u = data.dsol[0];
323 
324  TPZAxesTools<STATE>::Axes2XYZ(dsol_u, dsol_xy, data.axes);
325 
326  TPZManVector<STATE,3> p(m_f);
327 // TPZFNMatrix<4,STATE> grad_p(2,2,0.0);
328 
329  if(this->HasForcingFunction())
330  {
331  fForcingFunction->Execute(data.x,p);
332  }
333 
334  for(int iu = 0; iu < n_phi_u; iu++ )
335  {
336  TPZFNMatrix<4,REAL> gradv(2,2);
337  gradv(0,0) = dudx(0,iu);
338  gradv(0,1) = dudx(1,iu);
339  gradv(1,0) = dvdx(0,iu);
340  gradv(1,1) = dvdx(1,iu);
341 
342  // Vector Force right hand term
343  ef(iu) += weight*(p[0]*phi_u(2*iu, 0) + p[1]*phi_u(2*iu+1,0)
344  -m_s0_xx*gradv(0,0)-m_s0_yy*gradv(1,1)-m_s0_xy*(gradv(0,1)+gradv(1,0)));
345 
346  if (m_plane_stress == 1)
347  {
348  /* Plain stress state */
349  TPZFNMatrix<4,REAL> sigma_u(2,2);
350  sigma_u(0,0) = m_E/(1-m_nu*m_nu) *(dsol_xy(0,0)+m_nu*dsol_xy(1,1));
351  sigma_u(1,1) = m_E/(1-m_nu*m_nu) *(dsol_xy(1,1)+m_nu*dsol_xy(0,0));
352  sigma_u(0,1) = m_E/(2.*(1+m_nu))*(dsol_xy(0,1)+dsol_xy(1,0));
353  sigma_u(1,0) = sigma_u(0,1);
354 
355  ef(iu) -= weight*(gradv(0,0)*sigma_u(0,0)+gradv(1,1)*sigma_u(1,1)+gradv(1,0)*sigma_u(1,0)+gradv(0,1)*sigma_u(0,1));
356  }
357  else
358  {
359  /* Plain Strain State */
360  TPZFNMatrix<4,REAL> sigma_u(2,2);
361  sigma_u(0,0) = (m_lambda+2.*m_mu)*dsol_xy(0,0)+m_lambda*dsol_xy(1,1);
362  sigma_u(1,1) = (m_lambda+2.*m_mu)*dsol_xy(1,1)+m_lambda*dsol_xy(0,0);
363  sigma_u(0,1) = m_mu*(dsol_xy(0,1)+dsol_xy(1,0));
364  sigma_u(1,0) = sigma_u(0,1);
365  ef(iu) -= weight*(gradv(0,0)*sigma_u(0,0)+gradv(1,1)*sigma_u(1,1)+gradv(1,0)*sigma_u(1,0)+gradv(0,1)*sigma_u(0,1));
366 
367  }
368  }
369 
370  // ////////////////////////// Residual Vector ///////////////////////////////////
371 
372 
373 }
374 
375 void TPZMatElasticity2D::Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ef) {
376 
377  if (data.fShapeType == TPZMaterialData::EVecShape) {
378  ContributeVec(data, weight, ef);
379  return;
380  }
381 
382  // Getting weight functions
383  TPZFMatrix<REAL> &phi_u = data.phi;
384  int n_phi_u = phi_u.Rows();
385  int first_u = 0;
386 
387  TPZFNMatrix<40,REAL> grad_phi_u(3,n_phi_u);
388  TPZAxesTools<REAL>::Axes2XYZ(data.dphix, grad_phi_u, data.axes);
389  TPZFMatrix<STATE> dsol_u = data.dsol[0];
390 
391  TPZManVector<STATE,3> grad_p(m_f);
392 
393  if(this->HasForcingFunction())
394  {
395  fForcingFunction->Execute(data.x,grad_p);
396  }
397 
398  REAL dvdx,dvdy;
399  REAL duxdx,duxdy,duydx,duydy;
400 
401  // Gradient for ux
402  duxdx = dsol_u(0,0)*data.axes(0,0)+dsol_u(1,0)*data.axes(1,0); // dux/dx
403  duxdy = dsol_u(0,0)*data.axes(0,1)+dsol_u(1,0)*data.axes(1,1); // dux/dy
404 
405  // Gradient for uy
406  duydx = dsol_u(0,1)*data.axes(0,0)+dsol_u(1,1)*data.axes(1,0); // duy/dx
407  duydy = dsol_u(0,1)*data.axes(0,1)+dsol_u(1,1)*data.axes(1,1); // duy/dy
408 
409 
410  for(int iu = 0; iu < n_phi_u; iu++ )
411  {
412  dvdx = grad_phi_u(0,iu);
413  dvdy = grad_phi_u(1,iu);
414 
415  ef(2*iu + first_u) += weight * (grad_p[0] * phi_u(iu, 0) - (dvdx*m_s0_xx + dvdy*m_s0_xy) ); // x direction
416  ef(2*iu+1 + first_u) += weight * (grad_p[1] * phi_u(iu, 0) - (dvdx*m_s0_xy + dvdy*m_s0_yy) ); // y direction
417 
418  if (m_plane_stress == 1)
419  {
420  /* Plain stress state */
421  ef(2*iu + first_u) += weight*((4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dvdx*duxdx + (2*m_mu)*dvdy*duxdy);
422 
423  ef(2*iu + first_u) += weight*((2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dvdx*duydy + (2*m_mu)*dvdy*duydx);
424 
425  ef(2*iu+1 + first_u) += weight*((2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dvdy*duxdx + (2*m_mu)*dvdx*duxdy);
426 
427  ef(2*iu+1 + first_u) += weight*((4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dvdy*duydy + (2*m_mu)*dvdx*duydx);
428  }
429  else
430  {
431  /* Plain Strain State */
432  ef(2*iu + first_u) += weight* ((m_lambda + 2*m_mu)*dvdx*duxdx + (m_mu)*dvdy*(duxdy));
433 
434  ef(2*iu + first_u) += weight* (m_lambda*dvdx*duydy + (m_mu)*dvdy*(duydx));
435 
436  ef(2*iu+1 + first_u) += weight* (m_lambda*dvdy*duxdx + (m_mu)*dvdx*(duxdy));
437 
438  ef(2*iu+1 + first_u) += weight* ((m_lambda + 2*m_mu)*dvdy*duydy + (m_mu)*dvdx*(duydx));
439  }
440  }
441 
442  // ////////////////////////// Residual Vector ///////////////////////////////////
443 
444 }
445 
447 void TPZMatElasticity2D::ComputeSigma(const TPZFMatrix<STATE> &dudx, TPZFMatrix<STATE> &sigma)
448 {
449 #ifdef PZDEBUG
450  if (dudx.Rows() < 2 || dudx.Cols() != 2 || sigma.Rows() != 2 || sigma.Cols() != 2) {
451  DebugStop();
452  }
453 #endif
454 
455  if (m_plane_stress == 1)
456  {
457  sigma(0,0) = (4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dudx.g(0,0)+(2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dudx.g(1,1);
458  sigma(0,1) = (m_mu)*(dudx.g(1,0)+dudx.g(0,1));
459  sigma(1,0) = sigma(0,1);
460  sigma(1,1) = (2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu))*dudx.g(0,0)+(4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu))*dudx.g(1,1);
461  /* Plain stress state */
462  }
463  else
464  {
465  sigma(0,0) = (m_lambda + 2*m_mu)*dudx.g(0,0)+m_lambda*dudx.g(1,1);
466  sigma(1,0) = m_mu*(dudx.g(1,0)+dudx.g(0,1));
467  sigma(0,1) = sigma(1,0);
468  sigma(1,1) = (m_lambda + 2*m_mu)*dudx.g(1,1)+m_lambda*dudx.g(0,0);
469  /* Plain Strain State */
470  }
471 
472 }
473 
474 
475 
476 void TPZMatElasticity2D::ContributeBC(TPZMaterialData &data,REAL weight, TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
477 {
478 
479  TPZFMatrix<REAL> &phiu = data.phi;
480  TPZManVector<STATE,3> sol_u = data.sol[0];
481  TPZFMatrix<STATE> dsol_u = data.dsol[0];
482 
483  REAL ux = sol_u[0];
484  REAL uy = sol_u[1];
485 
486  TPZFNMatrix<4,STATE> val1loc(bc.Val1()),val2loc(bc.Val2());
487 
488  if (bc.HasForcingFunction()) {
489  TPZManVector<STATE,2> val2vec(2);
490 
491  bc.ForcingFunction()->Execute(data.x, val2vec, val1loc);
492  val2loc(0,0) = val2vec[0];
493  val2loc(1,0) = val2vec[1];
494  // we assume type 2 is displacement value is weakly imposed
495  if(bc.Type() == 2)
496  {
497  val1loc = bc.Val1();
498  for (int i=0; i<2; i++) {
499  val2loc(i,0) = 0.;
500  for (int j=0; j<2; j++) {
501  val2loc(i,0) += val1loc(i,j)*val2vec[j];
502  }
503  }
504  }
505  if(bc.Type() == 1)
506  {
507  for (int i=0; i<2; i++) {
508  val2loc(i,0) = 0.;
509  for (int j=0; j<2; j++) {
510  val2loc(i,0) += val1loc(i,j)*data.normal[j];
511  }
512  }
513  }
514  }
515  int phru = phiu.Rows();
516  short in,jn;
517  STATE v2[3];
518  TPZFMatrix<STATE> &v1 = val1loc;
519  v2[0] = val2loc(0,0); // Ux displacement or Tnx
520  v2[1] = val2loc(1,0); // Uy displacement or Tny
521 
522  // Here each digit represent an individual boundary condition corresponding to each state variable.
523  // 0 means Dirichlet condition on x-y
524  // 1 means Neumann condition
525  // 7 means Dirichlet condition on x
526  // 8 means Dirichlet condition on y
527 
528  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
529  switch (bc.Type())
530  {
531  case 0 :
532  {
533  // Dirichlet condition for each state variable
534  // Elasticity Equation
535  for(in = 0 ; in < phru; in++)
536  {
537  // Contribution for load Vector
538  ef(2*in,0) += BIGNUMBER*(ux - v2[0])*phiu(in,0)*weight; // X displacement Value
539  ef(2*in+1,0) += BIGNUMBER*(uy - v2[1])*phiu(in,0)*weight; // y displacement Value
540 
541  for (jn = 0 ; jn < phru; jn++)
542  {
543  // Contribution for Stiffness Matrix
544  ek(2*in,2*jn) += BIGNUMBER*phiu(in,0)*phiu(jn,0)*weight; // X displacement
545  ek(2*in+1,2*jn+1) += BIGNUMBER*phiu(in,0)*phiu(jn,0)*weight; // Y displacement
546  }
547  }
548 
549  break;
550  }
551  case 1 :
552  {
553  // Neumann condition for each state variable
554  // Elasticity Equation
555  for(in = 0 ; in <phru; in++)
556  {
557  // Normal Tension Components on neumann boundary
558  ef(2*in,0) += -1.0*v2[0]*phiu(in,0)*weight; // Tnx
559  ef(2*in+1,0) += -1.0*v2[1]*phiu(in,0)*weight; // Tny
560  }
561  break;
562  }
563  case 2 :
564  {
565  // Mixed condition for each state variable no used here
566  // Elasticity Equation
567  TPZFNMatrix<2,STATE> res(2,1,0.);
568  for(int i=0; i<2; i++) for(int j=0; j<2; j++)
569  {
570  res(i,0) += v1(i,j)*data.sol[0][j];
571  }
572 
573  for(in = 0 ; in < phru; in++)
574  {
575  ef(2*in+0,0) += weight * (v2[0]-res(0,0)) * phiu(in,0);
576  ef(2*in+1,0) += weight * (v2[1]-res(1,0)) * phiu(in,0);
577 
578  for (jn = 0 ; jn < phru; jn++)
579  {
580  for(int idf=0; idf < this->Dimension(); idf++) for(int jdf=0; jdf< this->Dimension(); jdf++)
581  {
582  ek(2*in+idf,2*jn+jdf) += v1(idf,jdf)*phiu(in,0)*phiu(jn,0)*weight;
583  // Not Complete with val2? HERE! PHIL!!!!
584  // DebugStop();
585  }
586  }
587  }
588 
589  break;
590  }
591  case 3 :
592  {
593  // Null Dirichlet condition for each state variable
594  // Elasticity Equation
595  for(in = 0 ; in < phru; in++)
596  {
597  // Contribution for load Vector
598  ef(2*in,0) += BIGNUMBER*( v2[0])*phiu(in,0)*weight; // X displacement Value
599  ef(2*in+1,0) += BIGNUMBER*( v2[1])*phiu(in,0)*weight; // y displacement Value
600 
601  for (jn = 0 ; jn < phru; jn++)
602  {
603  // Contribution for Stiffness Matrix
604  ek(2*in,2*jn) += BIGNUMBER*phiu(in,0)*phiu(jn,0)*weight; // X displacement
605  ek(2*in+1,2*jn+1) += BIGNUMBER*phiu(in,0)*phiu(jn,0)*weight; // Y displacement
606  }
607  }
608 
609  break;
610  }
611  case 4 :
612  {
613  // Stress Field as Neumann condition for each state variable
614  // Elasticity Equation
615 
616  for(in = 0; in < this->Dimension(); in ++){
617  v2[in] = ( v1(in,0) * data.normal[0] + v1(in,1) * data.normal[1]);
618  }
619 
620  for(in = 0 ; in <phru; in++)
621  {
622  // Normal Tension Components on neumann boundary
623  ef(2*in,0) += -1.0*v2[0]*phiu(in,0)*weight; // Tnx
624  ef(2*in+1,0) += -1.0*v2[1]*phiu(in,0)*weight; // Tny
625  }
626 
627  break;
628  }
629  case 5 :
630  // Normal Pressure condition Pressure value Should be inserted in v2[0]
631  // Elasticity Equation
632  {
633  TPZFNMatrix<2,STATE> res(2,1,0.);
634  for(int i=0; i<2; i++) for(int j=0; j<2; j++)
635  {
636  res(i,0) += data.normal[i]*bc.Val1()(i,j)*data.sol[0][j]*data.normal[j];
637  }
638  for(int in = 0 ; in < phru; in++)
639  {
640  ef(2*in+0,0) += (v2[0]*data.normal[0]-res(0,0)) * phiu(in,0) * weight ;
641  ef(2*in+1,0) += (v2[0]*data.normal[1]-res(1,0)) * phiu(in,0) * weight ;
642  for(int jn=0; jn< phru; jn++)
643  {
644  for(int idf=0; idf < this->Dimension(); idf++) for(int jdf=0; jdf < this->Dimension(); jdf++)
645  {
646  ek(2*in+idf,2*jn+jdf) += v1(idf,jdf)*data.normal[idf]*data.normal[jdf]*phiu(in,0)*phiu(jn,0)*weight;
647  // Not Complete with val2? HERE! PHIL!!!!
648  // DebugStop();
649  }
650  }
651  }
652  }
653  break;
654  case 6 :
655  // Normal Pressure condition Pressure value Should be inserted in v2[0]
656  // Elasticity Equation
657  {
658  TPZFNMatrix<2,STATE> res(2,1,0.);
659  for(int i=0; i<2; i++) for(int j=0; j<2; j++)
660  {
661  res(i,0) += data.normal[i]*bc.Val1()(i,j)*data.sol[0][j]*data.normal[j];
662  }
663  for(int in = 0 ; in < phru; in++)
664  {
665  ef(2*in+0,0) += (v2[0]*data.normal[0]-res(0,0)) * phiu(in,0) * weight ;
666  ef(2*in+1,0) += (v2[0]*data.normal[1]-res(1,0)) * phiu(in,0) * weight ;
667  for(int jn=0; jn< phru; jn++)
668  {
669  for(int idf=0; idf < this->Dimension(); idf++) for(int jdf=0; jdf < this->Dimension(); jdf++)
670  {
671  ek(2*in+idf,2*jn+jdf) += v1(idf,jdf)*data.normal[idf]*data.normal[jdf]*phiu(in,0)*phiu(jn,0)*weight;
672  // Not Complete
673  // DebugStop();
674  }
675  }
676  }
677  }
678  break;
679  case 7 :
680  {
681  // Dirichlet condition for each state variable
682  // Elasticity Equation
683  for(in = 0 ; in < phru; in++)
684  {
685  // Contribution for load Vector
686  ef(2*in,0) += BIGNUMBER*(ux - v2[0])*phiu(in,0)*weight; // X displacement Value
687 
688  for (jn = 0 ; jn < phru; jn++)
689  {
690  // Contribution for Stiffness Matrix
691  ek(2*in,2*jn) += BIGNUMBER*phiu(in,0)*phiu(jn,0)*weight; // X displacement
692  }
693  }
694 
695  break;
696  }
697  case 8 :
698  {
699  // Dirichlet condition for uy
700  // Elasticity Equation
701  for(in = 0 ; in < phru; in++)
702  {
703  // Contribution for load Vector
704  ef(2*in+1,0) += BIGNUMBER*(uy - v2[1])*phiu(in,0)*weight; // y displacement Value
705 
706  for (jn = 0 ; jn < phru; jn++)
707  {
708  // Contribution for Stiffness Matrix
709  ek(2*in+1,2*jn+1) += BIGNUMBER*phiu(in,0)*phiu(jn,0)*weight; // Y displacement
710  }
711  }
712 
713  break;
714  }
715  default:
716  {
717  PZError << "TPZMatElasticity2D::ContributeBC error - Wrong boundary condition type" << std::endl;
718  DebugStop();
719  }
720  break;
721  }
722 
723 }
724 
725 void TPZMatElasticity2D::ContributeBC(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
726 {
727 
728  TPZFMatrix<REAL> &phiu = data.phi;
729  TPZManVector<STATE,3> sol_u = data.sol[0];
730  TPZFMatrix<STATE> dsol_u = data.dsol[0];
731 
732  REAL ux = sol_u[0];
733  REAL uy = sol_u[1];
734 
735  int phru = phiu.Rows();
736  short in;
737  STATE v2[3]; TPZFMatrix<STATE> &v1 = bc.Val1();
738  v2[0] = bc.Val2()(0,0); // Ux displacement or Tnx
739  v2[1] = bc.Val2()(1,0); // Uy displacement or Tny
740 
741  // Here each digit represent an individual boundary condition corresponding to each state variable.
742  // 0 means Dirichlet condition on x-y
743  // 1 means Neumann condition
744  // 7 means Dirichlet condition on x
745  // 8 means Dirichlet condition on y
746 
747  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
748  switch (bc.Type())
749  {
750  case 0 :
751  {
752  // Dirichlet condition for each state variable
753  // Elasticity Equation
754  for(in = 0 ; in < phru; in++)
755  {
756  // Contribution for load Vector
757  ef(2*in,0) += BIGNUMBER*(ux - v2[0])*phiu(in,0)*weight; // X displacement Value
758  ef(2*in+1,0) += BIGNUMBER*(uy - v2[1])*phiu(in,0)*weight; // y displacement Value
759 
760  }
761 
762  break;
763  }
764  case 1 :
765  {
766  // Neumann condition for each state variable
767  // Elasticity Equation
768  for(in = 0 ; in <phru; in++)
769  {
770  // Normal Tension Components on neumann boundary
771  ef(2*in,0) += -1.0*v2[0]*phiu(in,0)*weight; // Tnx
772  ef(2*in+1,0) += -1.0*v2[1]*phiu(in,0)*weight; // Tny
773  }
774  break;
775  }
776  case 2 :
777  {
778  // Mixed condition for each state variable no used here
779  // Elasticity Equation
780  TPZFNMatrix<2,STATE> res(2,1,0.);
781  for(int i=0; i<2; i++) for(int j=0; j<2; j++)
782  {
783  res(i,0) += bc.Val1()(i,j)*data.sol[0][j];
784  }
785 
786  for(in = 0 ; in < phru; in++)
787  {
788  ef(2*in+0,0) += weight * (v2[0]-res(0,0)) * phiu(in,0);
789  ef(2*in+1,0) += weight * (v2[1]-res(1,0)) * phiu(in,0);
790 
791  }
792 
793  break;
794  }
795  case 3 :
796  {
797  // Null Dirichlet condition for each state variable
798  // Elasticity Equation
799  for(in = 0 ; in < phru; in++)
800  {
801  // Contribution for load Vector
802  ef(2*in,0) += BIGNUMBER*(0.0 - v2[0])*phiu(in,0)*weight; // X displacement Value
803  ef(2*in+1,0) += BIGNUMBER*(0.0 - v2[1])*phiu(in,0)*weight; // y displacement Value
804 
805  }
806 
807  break;
808  }
809  case 4 :
810  {
811  // Stress Field as Neumann condition for each state variable
812  // Elasticity Equation
813 
814  for(in = 0; in < this->Dimension(); in ++){ v2[in] = ( v1(in,0) * data.normal[0] + v1(in,1) * data.normal[1]);}
815 
816  for(in = 0 ; in <phru; in++)
817  {
818  // Normal Tension Components on neumann boundary
819  ef(2*in,0) += -1.0*v2[0]*phiu(in,0)*weight; // Tnx
820  ef(2*in+1,0) += -1.0*v2[1]*phiu(in,0)*weight; // Tny
821  }
822 
823  break;
824  }
825  case 5 :
826  // Normal Pressure condition Pressure value Should be inserted in v2[0]
827  // Elasticity Equation
828  {
829  TPZFNMatrix<2,STATE> res(2,1,0.);
830  for(int i=0; i<2; i++) for(int j=0; j<2; j++)
831  {
832  res(i,0) += data.normal[i]*bc.Val1()(i,j)*data.sol[0][j]*data.normal[j];
833  }
834  for(int in = 0 ; in < phru; in++)
835  {
836  ef(2*in+0,0) += (v2[0]*data.normal[0]-res(0,0)) * phiu(in,0) * weight ;
837  ef(2*in+1,0) += (v2[0]*data.normal[1]-res(1,0)) * phiu(in,0) * weight ;
838  }
839  }
840  break;
841  case 6 :
842  // Normal Pressure condition Pressure value Should be inserted in v2[0]
843  // Elasticity Equation
844  {
845  TPZFNMatrix<2,STATE> res(2,1,0.);
846  for(int i=0; i<2; i++) for(int j=0; j<2; j++)
847  {
848  res(i,0) += data.normal[i]*bc.Val1()(i,j)*data.sol[0][j]*data.normal[j];
849  }
850  for(int in = 0 ; in < phru; in++)
851  {
852  ef(2*in+0,0) += (v2[0]*data.normal[0]-res(0,0)) * phiu(in,0) * weight ;
853  ef(2*in+1,0) += (v2[0]*data.normal[1]-res(1,0)) * phiu(in,0) * weight ;
854  }
855  }
856  break;
857  case 7 :
858  {
859  // Dirichlet condition for each state variable
860  // Elasticity Equation
861  for(in = 0 ; in < phru; in++)
862  {
863  // Contribution for load Vector
864  ef(2*in,0) += BIGNUMBER*(ux - v2[0])*phiu(in,0)*weight; // X displacement Value
865  }
866 
867  break;
868  }
869  case 8 :
870  {
871  // Dirichlet condition for each state variable
872  // Elasticity Equation
873  for(in = 0 ; in < phru; in++)
874  {
875  // Contribution for load Vector
876  ef(2*in+1,0) += BIGNUMBER*(uy - v2[1])*phiu(in,0)*weight; // y displacement Value
877  }
878 
879  break;
880  }
881  default:
882  {
883  PZError << "TPZMatElasticity2D::ContributeBC error - Wrong boundary condition type" << std::endl;
884  DebugStop();
885  }
886  break;
887  }
888 
889 }
890 
891 
892 void TPZMatElasticity2D::FillDataRequirements(TPZMaterialData &data)
893 {
894  data.SetAllRequirements(false);
895  data.fNeedsSol = true;
896  data.fNeedsNeighborSol = true;
897  data.fNeedsNeighborCenter = false;
898  data.fNeedsNormal = true;
899 }
900 
901 void TPZMatElasticity2D::FillBoundaryConditionDataRequirement(int type, TPZMaterialData &data){
902  data.SetAllRequirements(false);
903  data.fNeedsSol = true;
904  data.fNeedsNormal = true;
905 }
906 
907 
908 void TPZMatElasticity2D::Print(std::ostream &out)
909 {
910  out << "Material Name : " << Name() << "\n";
911  out << "Plane Problem (m_plane_stress = 0, for Plane Strain conditions) " << m_plane_stress << std::endl;
912  out << "Properties for elasticity: \n";
913  out << "\t Young modulus = " << m_E << std::endl;
914  out << "\t Poisson Ratio = " << m_nu << std::endl;
915  out << "\t First Lamé Parameter = " << m_lambda << std::endl;
916  out << "\t Second Lamé Parameter = " << m_mu << std::endl;
917  out << "\t Body force vector B {X-direction, Y-direction} = " << m_f[0] << ' ' << m_f[1] << std::endl;
918  out << "\t m_s0_xx = " << m_s0_xx << std::endl;
919  out << "\t m_s0_xy = " << m_s0_xy << std::endl;
920  out << "\t m_s0_yy = " << m_s0_yy << std::endl;
921  out << "\t m_s0_zz = " << m_s0_zz << std::endl;
922  out << "Class properties :";
923  TPZMaterial::Print(out);
924  out << "\n";
925 
926 }
927 
929 int TPZMatElasticity2D::VariableIndex(const std::string &name)
930 {
931  // Elasticity Variables
932  if(!strcmp("Displacement",name.c_str())) return 1;
933  if(!strcmp("SolidPressure",name.c_str())) return 2;
934  if(!strcmp("SigmaX",name.c_str())) return 3;
935  if(!strcmp("SigmaY",name.c_str())) return 4;
936  if(!strcmp("SigmaZ",name.c_str())) return 5;
937  if(!strcmp("TauXY",name.c_str())) return 6;
938  if(!strcmp("EpsX",name.c_str())) return 7;
939  if(!strcmp("EpsY",name.c_str())) return 8;
940  if(!strcmp("EpsZ",name.c_str())) return 9;
941  if(!strcmp("EpsXY",name.c_str())) return 10;
942 // PZError << "TPZMatElasticity2D::VariableIndex Error\n";
943 
944  return TPZMaterial::VariableIndex(name);
945 }
946 
950 void TPZMatElasticity2D::Write(TPZStream &buf, int withclassid) const
951 {
952  TPZMaterial::Write(buf,withclassid);
953  buf.Write(&m_E);
954  buf.Write(&m_nu);
955  buf.Write(&m_lambda);
956  buf.Write(&m_mu);
957  buf.Write( m_f);
958  buf.Write(&m_s0_xx);
959  buf.Write(&m_s0_xy);
960  buf.Write(&m_s0_yy);
961  buf.Write(&m_s0_zz);
962  buf.Write(&m_plane_stress);
963 
964 }
965 
969 void TPZMatElasticity2D::Read(TPZStream &buf, void *context)
970 {
971  TPZMaterial::Read(buf,context);
972  buf.Read(&m_E);
973  buf.Read(&m_nu);
974  buf.Read(&m_lambda);
975  buf.Read(&m_mu);
976  buf.Read( m_f);
977  buf.Read(&m_s0_xx);
978  buf.Read(&m_s0_xy);
979  buf.Read(&m_s0_yy);
980  buf.Read(&m_s0_zz);
981  buf.Read(&m_plane_stress);
982 
983 }
984 
985 int TPZMatElasticity2D::NSolutionVariables(int var){
986  if(var == 1) return 3;
987  if(var == 2) return 1;
988  if(var == 3) return 1;
989  if(var == 4) return 1;
990  if(var == 5) return 1;
991  if(var == 6) return 1;
992  if(var == 7) return 1;
993  if(var == 8) return 1;
994  if(var == 9) return 1;
995  if(var == 10) return 1;
996 
997  return TPZMaterial::NSolutionVariables(var);
998 }
999 
1000 // Calculate Secondary variables based on ux, uy, Pore pressure and their derivatives
1001 void TPZMatElasticity2D::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout){
1002 
1003  Solout.Resize(this->NSolutionVariables(var));
1004 
1005  TPZManVector<STATE,3> SolU, SolP;
1006  TPZFNMatrix <6,STATE> DSolU, DSolP;
1007  TPZFNMatrix <9> axesU, axesP;
1008 
1009  TPZVec<REAL> ptx(3);
1010  TPZVec<STATE> solExata(3);
1011  TPZFMatrix<STATE> flux(5,1);
1012 
1013  if (data.sol.size() != 1) {
1014  DebugStop();
1015  }
1016 
1017  SolU = data.sol[0];
1018  DSolU = data.dsol[0];
1019  axesU = data.axes;
1020 
1021 
1022  // Displacements
1023  if(var == 1 || var == 0){
1024  Solout[0] = SolU[0];
1025  Solout[1] = SolU[1];
1026  if(var==1) Solout[2] = 0.0;
1027  return;
1028  }
1029 
1030 
1031  REAL epsx;
1032  REAL epsy;
1033  REAL epsxy;
1034  REAL SigX;
1035  REAL SigY;
1036  REAL SigZ;
1037  REAL Tau, DSolxy[2][2];
1038  REAL divu;
1039 
1040  DSolxy[0][0] = DSolU(0,0)*axesU(0,0)+DSolU(1,0)*axesU(1,0); // dUx/dx
1041  DSolxy[1][0] = DSolU(0,0)*axesU(0,1)+DSolU(1,0)*axesU(1,1); // dUx/dy
1042 
1043  DSolxy[0][1] = DSolU(0,1)*axesU(0,0)+DSolU(1,1)*axesU(1,0); // dUy/dx
1044  DSolxy[1][1] = DSolU(0,1)*axesU(0,1)+DSolU(1,1)*axesU(1,1); // dUy/dy
1045 
1046  divu = DSolxy[0][0]+DSolxy[1][1]+0.0;
1047 
1048  epsx = DSolxy[0][0];// du/dx
1049  epsy = DSolxy[1][1];// dv/dy
1050  epsxy = 0.5*(DSolxy[1][0]+DSolxy[0][1]);
1051  REAL C11 = 4*(m_mu)*(m_lambda+m_mu)/(m_lambda+2*m_mu);
1052  REAL C22 = 2*(m_mu)*(m_lambda)/(m_lambda+2*m_mu);
1053 
1054  if (this->m_plane_stress)
1055  {
1056  SigX = C11*epsx+C22*epsy;
1057  SigY = C11*epsy+C22*epsx;
1058  SigZ = 0.0;
1059  Tau = 2.0*m_mu*epsxy;
1060  }
1061  else
1062  {
1063  SigX = ((m_lambda + 2*m_mu)*(epsx) + (m_lambda)*epsy);
1064  SigY = ((m_lambda + 2*m_mu)*(epsy) + (m_lambda)*epsx);
1065  SigZ = m_lambda*divu;
1066  Tau = 2.0*m_mu*epsxy;
1067  }
1068 
1069 
1070  // Hydrostatic stress
1071  if(var == 2)
1072  {
1073  Solout[0] = SigX+SigY+SigZ;
1074  return;
1075  }
1076 
1077  // Effective Stress x-direction
1078  if(var == 3) {
1079  Solout[0] = SigX + m_s0_xx;
1080  return;
1081  }
1082 
1083  // Effective Stress y-direction
1084  if(var == 4) {
1085  Solout[0] = SigY + m_s0_yy;
1086  return;
1087  }
1088 
1089  // Effective Stress y-direction
1090  if(var == 5) {
1091  Solout[0] = SigZ + m_s0_zz;
1092  return;
1093  }
1094 
1095  // Shear Stress
1096  if(var == 6) {
1097  Solout[0] = Tau + m_s0_xy;
1098  return;
1099  }
1100 
1101  // epsx
1102  if (var == 7) {
1103  Solout[0] = epsx;
1104  }
1105 
1106  // epsy
1107  if (var == 8) {
1108  Solout[0] = epsy;
1109  }
1110 
1111  // epsz
1112  if (var == 9) {
1113  if (m_plane_stress) {
1114  Solout[0] = -m_nu*(epsx+epsy);
1115  }
1116  else
1117  {
1118  Solout[0] = 0.;
1119  }
1120  }
1121 
1122  // epsxy
1123  if (var == 10) {
1124  Solout[0] = epsxy;
1125  }
1126 }
1127 
1128 int TPZMatElasticity2D::ClassId() const{
1129  return Hash("TPZMatElasticity2D") ^ TPZMaterial::ClassId() << 1;
1130 }
1131 
1136 void TPZMatElasticity2D::Errors(TPZVec<REAL> &x, TPZVec<STATE> &sol, TPZFMatrix<STATE> &dsol,
1137  TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux,
1138  TPZVec<STATE> &uexact, TPZFMatrix<STATE> &duexact,
1139  TPZVec<REAL> &val)
1140 {
1141  TPZFNMatrix<9,STATE> dudx(2,2), stress(2,2), stressexact(2,2);
1142  TPZAxesTools<STATE>::Axes2XYZ(dsol, dudx, axes);
1143  ComputeSigma(dudx, stress);
1144  ComputeSigma(duexact, stressexact);
1145  REAL L2 = 0.;
1146  L2 = (sol[0]-uexact[0])*(sol[0]-uexact[0])+(sol[1]-uexact[1])*(sol[1]-uexact[1]);
1147  REAL H1 = 0.;
1148  REAL energy = 0.;
1149  for (int i=0; i<2; i++) {
1150  for (int j=0; j<2; j++) {
1151  H1 += (dudx(i,j)-duexact(i,j))*(dudx(i,j)-duexact(i,j));
1152  energy += (stress(i,j)-stressexact(i,j))*(dudx(i,j)-duexact(i,j));
1153  }
1154  }
1155  val[0] = energy;
1156  val[1] = L2;
1157  val[2] = H1;
1158 }
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZMaterialData::normal
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
Definition: pzmaterialdata.h:69
TPZMatElasticity2D::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Definition: TPZMatElasticity2D.cpp:950
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.
TPZMatElasticity2D::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, 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...
Definition: TPZMatElasticity2D.cpp:476
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZMatElasticity2D::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: TPZMatElasticity2D.cpp:116
TPZMatElasticity2D::m_E
REAL m_E
Elasticity modulus.
Definition: TPZMatElasticity2D.h:40
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZMatElasticity2D::operator=
TPZMatElasticity2D & operator=(const TPZMatElasticity2D &copy)
Definition: TPZMatElasticity2D.cpp:97
TPZVec::resize
virtual void resize(const int64_t newsize)
Definition: pzvec.h:213
TPZMatElasticity2D::Print
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
Definition: TPZMatElasticity2D.cpp:908
TPZMatElasticity2D::Name
virtual std::string Name() override
Returns the name of the material.
Definition: TPZMatElasticity2D.h:95
TPZMatElasticity2D::m_plane_stress
int m_plane_stress
plain stress directive
Definition: TPZMatElasticity2D.h:59
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
TPZMaterial::operator=
TPZMaterial & operator=(const TPZMaterial &copy)
operator =
Definition: TPZMaterial.cpp:62
TPZMatElasticity2D::m_s0_xx
REAL m_s0_xx
Initial Stress.
Definition: TPZMatElasticity2D.h:52
TPZMaterial::VariableIndex
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Definition: TPZMaterial.cpp:141
TPZMaterialData::SetAllRequirements
void SetAllRequirements(bool set)
Set all flags at once.
Definition: pzmaterialdata.cpp:90
TPZMatElasticity2D::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: TPZMatElasticity2D.cpp:1128
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
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMatElasticity2D::Dimension
int Dimension() const override
Returns the integrable dimension of the material.
Definition: TPZMatElasticity2D.h:97
TPZMatElasticity2D::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the.
Definition: TPZMatElasticity2D.cpp:892
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
TPZMatElasticity2D::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMatElasticity2D.cpp:985
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
TPZMatElasticity2D::Errors
virtual void Errors(TPZVec< REAL > &x, TPZVec< STATE > &sol, TPZFMatrix< STATE > &dsol, TPZFMatrix< REAL > &axes, TPZVec< STATE > &flux, TPZVec< STATE > &uexact, TPZFMatrix< STATE > &duexact, TPZVec< REAL > &val) override
Computes the error due to the difference between the interpolated flux and the flux computed based o...
Definition: TPZMatElasticity2D.cpp:1136
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
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
TPZMatElasticity2D::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: TPZMatElasticity2D.cpp:1001
TPZMaterial::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: TPZMaterial.cpp:402
TPZMatElasticity2D::Read
void Read(TPZStream &buf, void *context) override
Definition: TPZMatElasticity2D.cpp:969
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZMatElasticity2D::ComputeSigma
void ComputeSigma(const TPZFMatrix< STATE > &dudx, TPZFMatrix< STATE > &sigma)
compute the stress tensor as a function of the solution gradient
Definition: TPZMatElasticity2D.cpp:447
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZMatElasticity2D::m_nu
REAL m_nu
Poison coeficient.
Definition: TPZMatElasticity2D.h:43
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
TPZMatElasticity2D::ContributeVec
void ContributeVec(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: TPZMatElasticity2D.cpp:221
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
TPZMatElasticity2D::m_lambda
REAL m_lambda
first Lame Parameter
Definition: TPZMatElasticity2D.h:46
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
TPZMaterial::ClassId
int ClassId() const override
Unique identifier for serialization purposes.
Definition: TPZMaterial.cpp:382
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
TPZMatElasticity2D::m_mu
REAL m_mu
Second Lame Parameter.
Definition: TPZMatElasticity2D.h:49
TPZMatElasticity2D::m_f
TPZManVector< STATE, 2 > m_f
Forcing vector.
Definition: TPZMatElasticity2D.h:37
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
pzaxestools.h
Contains declaration of the TPZAxesTools class which implements verifications over axes...
TPZMatElasticity2D::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: TPZMatElasticity2D.cpp:901
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZMatElasticity2D::VariableIndex
virtual int VariableIndex(const std::string &name) override
Definition: TPZMatElasticity2D.cpp:929
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZMatElasticity2D::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 to multip...
Definition: TPZMatElasticity2D.cpp:121
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZMatElasticity2D
Description Linear elastic equations.
Definition: TPZMatElasticity2D.h:32
TPZFMatrix::g
TVar & g(const int64_t row, const int64_t col) const
Definition: pzfmatrix.h:594
TPZMatElasticity2D.h
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:04 for NeoPZ by   doxygen 1.8.13