NeoPZ
pzmathyperelastic.cpp
Go to the documentation of this file.
1 
6 #include "pzmathyperelastic.h"
7 #include "pzelmat.h"
8 #include "pzbndcond.h"
9 #include "pzmatrix.h"
10 #include "pzfmatrix.h"
11 #include "pzerror.h"
12 #include <math.h>
13 
14 #include <cmath>
15 using namespace std;
16 
17 TPZMatHyperElastic::TPZMatHyperElastic(int nummat,STATE e,STATE mu,STATE nu,
18  STATE lambda,STATE coef1,STATE coef2,STATE coef3) :
20 TPZMaterial(nummat)
21 {
22 
23  fE = e;
24  fMu = mu;
25  fNu = .5*fE/(1.+fMu);
26  fLambda = 2.*fNu*fMu/(1.-2.*fMu);
27  if(coef1 != -1.0) fCoef1 = coef1;
28  else fCoef1 = .25*fLambda;//=a
29  if(coef2 != -1.0) fCoef2 = coef2;
30  else fCoef2 = -(.5*fLambda+fNu);//=b
31  if(coef3 != -1.0) fCoef3 = coef3;
32  else fCoef3 = .5*fNu;//=c
33  fE1[0] = 2.;
34  fE1[1] = 2.;
35  fE1[2] = 2.;
36  fE5[0] = 2.;
37  fE5[1] = 2.;
38  fE5[2] = 2.;
39  fE9[0] = 2.;
40  fE9[1] = 2.;
41  fE9[2] = 2.;
42 }
43 
45 }
46 
48  return 3;//3 deslocamentos
49 }
50 
51 void TPZMatHyperElastic::Print(std::ostream &out) {
52  out << "name of material : " << Name() << "\n";
53  out << "properties : \n";
54  TPZMaterial::Print(out);
55 }
56 
58  TPZFMatrix<REAL> &dphi = data.dphix;
59  TPZFMatrix<REAL> &phi = data.phi;
60  TPZManVector<REAL,3> &x = data.x;
61  int numbersol = data.dsol.size();
62  if (numbersol != 1) {
63  DebugStop();
64  }
65  TPZFMatrix<STATE> &dsol=data.dsol[0];
66 
67  if(fForcingFunction) {
69  fForcingFunction->Execute(x,res);
70  fXf[0] = res[0];
71  fXf[1] = res[1];
72  fXf[2] = res[2];
73  }
74 #ifdef _AUTODIFF
76  ComputeEnergy(fLambda,fNu,dsol,U);
77  int nshape = phi.Rows();
78  int ish,jsh, i,j;
79  for(ish=0; ish<nshape; ish++) {
80  for(i=0; i<3; i++) {
81  ef(ish*3+i) -= (U.val().d(i)*dphi(0,ish)+U.val().d(3+i)*dphi(1,ish)+U.val().d(6+i)*dphi(2,ish))*weight;
82  for(jsh=0; jsh<nshape; jsh++) {
83  for(j=0; j<3; j++) {
84  ek(ish*3+i,jsh*3+j) += (
85  U.d(i).d(j )*dphi(0,ish)*dphi(0,jsh)+U.d(i+3).d(j )*dphi(1,ish)*dphi(0,jsh)+U.d(i+6).d(j )*dphi(2,ish)*dphi(0,jsh)+
86  U.d(i).d(j+3)*dphi(0,ish)*dphi(1,jsh)+U.d(i+3).d(j+3)*dphi(1,ish)*dphi(1,jsh)+U.d(i+6).d(j+3)*dphi(2,ish)*dphi(1,jsh)+
87  U.d(i).d(j+6)*dphi(0,ish)*dphi(2,jsh)+U.d(i+3).d(j+6)*dphi(1,ish)*dphi(2,jsh)+U.d(i+6).d(j+6)*dphi(2,ish)*dphi(2,jsh)
88  )*weight;
89  }
90  }
91  }
92  }
93 /* for(int ii=0;ii<3;ii++)
94  for(int jj=0;jj<3;jj++) {
95  fK2[ii][jj] = (STATE)0.;
96  fK3[ii][jj] = (STATE)0.;
97  fK4[ii][jj] = (STATE)0.;
98  fK6[ii][jj] = (STATE)0.;
99  fK7[ii][jj] = (STATE)0.;
100  fK8[ii][jj] = (STATE)0.;
101  }
102 */
103 #else
104  int i;
105  STATE global[3][3][9];
106  STATE ux,uy,uz,vx,vy,vz,wx,wy,wz;
107  ux = dsol(0,0);
108  uy = dsol(1,0);
109  uz = dsol(2,0);
110 
111  vx = dsol(0,1);
112  vy = dsol(1,1);
113  vz = dsol(2,1);
114 
115  wx = dsol(0,2);
116  wy = dsol(1,2);
117  wz = dsol(2,2);
118 
119  for(i=0; i<3; i++) fK2[i][i] = 0.;
120  fK2[0][0] = 0.;
121  fK2[0][1] = wz+1.;
122  fK2[0][2] = -vz;
123  fK2[1][0] = -wz-1.;
124  fK2[1][1] = 0.;
125  fK2[1][2] = uz;
126  fK2[2][0] = vz;
127  fK2[2][1] = -uz;
128  fK2[2][2] = 0.;
129 
130  for(i=0; i<3; i++) fK3[i][i] = 0.;
131  fK3[0][1] = -wy;
132  fK3[0][2] = vy+1.;
133  fK3[1][0] = wy;
134  fK3[1][2] = -uy;
135  fK3[2][0] = -vy-1.;
136  fK3[2][1] = uy;
137 
138  for(i=0; i<3; i++) fK4[i][i] = 0.;
139  fK4[0][1] = -wz-1.;
140  fK4[0][2] = vz;
141  fK4[1][0] = wz+1.;
142  fK4[1][2] = -uz;
143  fK4[2][0] = -vz;
144  fK4[2][1] = uz;
145 
146  for(i=0; i<3; i++) fK6[i][i] = 0.;
147  fK6[0][1] = wx;
148  fK6[0][2] = -vx;
149  fK6[1][0] = -wx;
150  fK6[1][2] = ux+1.;
151  fK6[2][0] = vx;
152  fK6[2][1] = -ux-1.;
153 
154  fK7[0][0] = fK7[1][1] = fK7[2][2] = 0.;
155  fK7[0][1] = wy;
156  fK7[0][2] = -vy-1.;
157  fK7[1][0] = -wy;
158  fK7[1][2] = uy;
159  fK7[2][0] = vy+1.;
160  fK7[2][1] = -uy;
161 
162  for(i=0; i<3; i++) fK8[i][i] = 0.;
163  fK8[0][1] = -wx;
164  fK8[0][2] = vx;
165  fK8[1][0] = wx;
166  fK8[1][2] = -ux-1.;
167  fK8[2][0] = -vx;
168  fK8[2][1] = ux+1.;
169 
170  fGradDetF[0][0] = (vy+1.)*(wz+1.) - wy*vz;//dux
171  fGradDetF[1][0] = wx*vz - vx*(wz+1.); //duy
172  fGradDetF[2][0] = vx*wy - wx*(vy+1.); //duz
173  fGradDetF[0][1] = wy*uz - uy*(wz+1.); //dvx
174  fGradDetF[1][1] = (ux+1.)*(wz+1.) - wx*uz;//dvy
175  fGradDetF[2][1] = wx*uy - (ux+1.)*wy; //dvz
176  fGradDetF[0][2] = uy*vz - (vy+1.)*uz; //dwx
177  fGradDetF[1][2] = vx*uz - (ux+1.)*vz; //dwy
178  fGradDetF[2][2] = (ux+1.)*(vy+1.) - vx*uy;//dwz
179 
180  fGradtrC[0][0] = 2.*(ux+1.);//dux
181  fGradtrC[1][0] = 2.*uy; //duy | dux dvx dwx |
182  fGradtrC[2][0] = 2.*uz; //duz : | duy dvy dwy |
183  fGradtrC[0][1] = 2.*vx; //dvx | duz dvz dwz |
184  fGradtrC[1][1] = 2.*(vy+1.);//dvy
185  fGradtrC[2][1] = 2.*vz; //dvz
186  fGradtrC[0][2] = 2.*wx; //dwx
187  fGradtrC[1][2] = 2.*wy; //dwy
188  fGradtrC[2][2] = 2.*(wz+1.);//dwz
189 
190  //foi transposta da original
191  fL1[0][0] = fGradDetF[0][0]*fGradDetF[0][0];//dux*dux
192  fL1[0][1] = fGradDetF[0][0]*fGradDetF[0][1];//dux*dvx
193  fL1[0][2] = fGradDetF[0][0]*fGradDetF[0][2];//dux*dwx
194  fL1[1][0] = fGradDetF[0][1]*fGradDetF[0][0];//dvx*dux
195  fL1[1][1] = fGradDetF[0][1]*fGradDetF[0][1];//dvx*dvx
196  fL1[1][2] = fGradDetF[0][1]*fGradDetF[0][2];//dvx*dwx
197  fL1[2][0] = fGradDetF[0][2]*fGradDetF[0][0];//dwx*dux
198  fL1[2][1] = fGradDetF[0][2]*fGradDetF[0][1];//dwx*dvx
199  fL1[2][2] = fGradDetF[0][2]*fGradDetF[0][2];//dwx*dwx
200 
201  fL2[0][0] = fGradDetF[0][0]*fGradDetF[1][0];//duy*dux
202  fL2[0][1] = fGradDetF[0][0]*fGradDetF[1][1];//dvy*dux
203  fL2[0][2] = fGradDetF[0][0]*fGradDetF[1][2];//dwy*dux
204  fL2[1][0] = fGradDetF[0][1]*fGradDetF[1][0];//duy*dvx
205  fL2[1][1] = fGradDetF[0][1]*fGradDetF[1][1];//dvy*dvx
206  fL2[1][2] = fGradDetF[0][1]*fGradDetF[1][2];//dwy*dvx
207  fL2[2][0] = fGradDetF[0][2]*fGradDetF[1][0];//duy*dwx
208  fL2[2][1] = fGradDetF[0][2]*fGradDetF[1][1];//dvy*dwx
209  fL2[2][2] = fGradDetF[0][2]*fGradDetF[1][2];//dwy*dwx
210 
211  fL3[0][0] = fGradDetF[0][0]*fGradDetF[2][0];//duz*dux
212  fL3[0][1] = fGradDetF[0][0]*fGradDetF[2][1];//dvz*dux
213  fL3[0][2] = fGradDetF[0][0]*fGradDetF[2][2];//dwz*dux
214  fL3[1][0] = fGradDetF[0][1]*fGradDetF[2][0];//duz*dvx
215  fL3[1][1] = fGradDetF[0][1]*fGradDetF[2][1];//dvz*dvx
216  fL3[1][2] = fGradDetF[0][1]*fGradDetF[2][2];//duz*dwx
217  fL3[2][0] = fGradDetF[0][2]*fGradDetF[2][0];//dvz*dwx
218  fL3[2][1] = fGradDetF[0][2]*fGradDetF[2][1];//dwz*dvx
219  fL3[2][2] = fGradDetF[0][2]*fGradDetF[2][2];//dwz*dwx
220 
221  fL4[0][0] = fGradDetF[1][0]*fGradDetF[0][0];//dux*duy
222  fL4[0][1] = fGradDetF[1][0]*fGradDetF[0][1];//dvx*duy
223  fL4[0][2] = fGradDetF[1][0]*fGradDetF[0][2];//dwx*duy
224  fL4[1][0] = fGradDetF[1][1]*fGradDetF[0][0];//dux*dvy
225  fL4[1][1] = fGradDetF[1][1]*fGradDetF[0][1];//dvx*dvy
226  fL4[1][2] = fGradDetF[1][1]*fGradDetF[0][2];//dwx*dvy
227  fL4[2][0] = fGradDetF[1][2]*fGradDetF[0][0];//dux*dwy
228  fL4[2][1] = fGradDetF[1][2]*fGradDetF[0][1];//dvx*dwy
229  fL4[2][2] = fGradDetF[1][2]*fGradDetF[0][2];//dwx*dwy
230 
231  fL5[0][0] = fGradDetF[1][0]*fGradDetF[1][0];//duy*duy
232  fL5[0][1] = fGradDetF[1][0]*fGradDetF[1][1];//dvy*duy
233  fL5[0][2] = fGradDetF[1][0]*fGradDetF[1][2];//dwy*duy
234  fL5[1][0] = fGradDetF[1][1]*fGradDetF[1][0];//duy*dvy
235  fL5[1][1] = fGradDetF[1][1]*fGradDetF[1][1];//dvy*dvy
236  fL5[1][2] = fGradDetF[1][1]*fGradDetF[1][2];//dwy*dvy
237  fL5[2][0] = fGradDetF[1][2]*fGradDetF[1][0];//duy*dwy
238  fL5[2][1] = fGradDetF[1][2]*fGradDetF[1][1];//dvy*dwy
239  fL5[2][2] = fGradDetF[1][2]*fGradDetF[1][2];//dwy*dwy
240 
241  fL6[0][0] = fGradDetF[1][0]*fGradDetF[2][0];//duz*duy
242  fL6[0][1] = fGradDetF[1][0]*fGradDetF[2][1];//dvz*duy
243  fL6[0][2] = fGradDetF[1][0]*fGradDetF[2][2];//dwz*duy
244  fL6[1][0] = fGradDetF[1][1]*fGradDetF[2][0];//duz*dvy
245  fL6[1][1] = fGradDetF[1][1]*fGradDetF[2][1];//dvz*dvy
246  fL6[1][2] = fGradDetF[1][1]*fGradDetF[2][2];//dwz*dvy
247  fL6[2][0] = fGradDetF[1][2]*fGradDetF[2][0];//duz*dwy
248  fL6[2][1] = fGradDetF[1][2]*fGradDetF[2][1];//dvz*dwy
249  fL6[2][2] = fGradDetF[1][2]*fGradDetF[2][2];//dwz*dwy
250 
251  fL7[0][0] = fGradDetF[2][0]*fGradDetF[0][0];//dux*duz
252  fL7[0][1] = fGradDetF[2][0]*fGradDetF[0][1];//dvx*duz
253  fL7[0][2] = fGradDetF[2][0]*fGradDetF[0][2];//dwx*duz
254  fL7[1][0] = fGradDetF[2][1]*fGradDetF[0][0];//dux*dvz
255  fL7[1][1] = fGradDetF[2][1]*fGradDetF[0][1];//dvx*dvz
256  fL7[1][2] = fGradDetF[2][1]*fGradDetF[0][2];//dwx*dvz
257  fL7[2][0] = fGradDetF[2][2]*fGradDetF[0][0];//dux*dwz
258  fL7[2][1] = fGradDetF[2][2]*fGradDetF[0][1];//dvx*dwz
259  fL7[2][2] = fGradDetF[2][2]*fGradDetF[0][2];//dwx*dwz
260 
261  fL8[0][0] = fGradDetF[2][0]*fGradDetF[1][0];//duy*duz
262  fL8[0][1] = fGradDetF[2][0]*fGradDetF[1][1];//dvy*duz
263  fL8[0][2] = fGradDetF[2][0]*fGradDetF[1][2];//dwy*duz
264  fL8[1][0] = fGradDetF[2][1]*fGradDetF[1][0];//duy*dvz
265  fL8[1][1] = fGradDetF[2][1]*fGradDetF[1][1];//dvy*dvz
266  fL8[1][2] = fGradDetF[2][1]*fGradDetF[1][2];//dwy*dvz
267  fL8[2][0] = fGradDetF[2][2]*fGradDetF[1][0];//duy*dwz
268  fL8[2][1] = fGradDetF[2][2]*fGradDetF[1][1];//dvy*dwz
269  fL8[2][2] = fGradDetF[2][2]*fGradDetF[1][2];//dwy*dwz
270 
271  fL9[0][0] = fGradDetF[2][0]*fGradDetF[2][0];//duz*duz
272  fL9[0][1] = fGradDetF[2][0]*fGradDetF[2][1];//dvz*duz
273  fL9[0][2] = fGradDetF[2][0]*fGradDetF[2][2];//dwz*duz
274  fL9[1][0] = fGradDetF[2][1]*fGradDetF[2][0];//duz*dvz
275  fL9[1][1] = fGradDetF[2][1]*fGradDetF[2][1];//dvz*dvz
276  fL9[1][2] = fGradDetF[2][1]*fGradDetF[2][2];//dwz*dvz
277  fL9[2][0] = fGradDetF[2][2]*fGradDetF[2][0];//duz*dwz
278  fL9[2][1] = fGradDetF[2][2]*fGradDetF[2][1];//dvz*dwz
279  fL9[2][2] = fGradDetF[2][2]*fGradDetF[2][2];//dwz*dwz
280 
281 
282  STATE detF = (ux+1.)*(vy+1.)*(wz+1.) + vx*wy*uz + wx*uy*vz - wx*(vy+1.)*uz - (ux+1.)*wy*vz - vx*uy*(wz+1.);
283 
284  if(detF < 0) {
285  cout << "\nDeterminante negativo!\n";
286  }
287 
288  STATE fac1 = 2.*fCoef1 - fCoef2/detF/detF;//A
289  STATE fac2 = 2.*fCoef1*detF + fCoef2/detF;//B
290  STATE c = fCoef3;//C
291 
292  int j;
293  for(i=0; i<3; i++) {
294  for(j=0; j<3; j++) {
295  global[i][j][0] = fac1*fL1[i][j];
296  global[i][j][1] = fac1*fL2[i][j]+fac2*fK2[i][j];
297  global[i][j][2] = fac1*fL3[i][j]+fac2*fK3[i][j];
298  global[i][j][3] = fac1*fL4[i][j]+fac2*fK4[i][j];
299  global[i][j][4] = fac1*fL5[i][j];
300  global[i][j][5] = fac1*fL6[i][j]+fac2*fK6[i][j];
301  global[i][j][6] = fac1*fL7[i][j]+fac2*fK7[i][j];
302  global[i][j][7] = fac1*fL8[i][j]+fac2*fK8[i][j];
303  global[i][j][8] = fac1*fL9[i][j];
304  }
305  global[i][i][0] += c*fE1[i];
306  global[i][i][4] += c*fE5[i];
307  global[i][i][8] += c*fE9[i];
308  }
309 
310  //AT�AQUI AS TRES MATRIZES A INTEGRAR S� fK , fL E fE , SENDO fE CONSTANTE
311  STATE gradJx[3],gradJy[3],gradJz[3];
312  STATE gradtrCx[3],gradtrCy[3],gradtrCz[3];
313  int k,l;
314  int nshape = phi.Rows();
315  for(i=0; i<3; i++) {
316  gradJx[i] = fGradDetF[0][i];//dux , dvx , dwx
317  gradJy[i] = fGradDetF[1][i];//duy , dvy , dwy
318  gradJz[i] = fGradDetF[2][i];//duz , dvz , dwz
319  gradtrCx[i] = fGradtrC[0][i];//dux , dvx , dwx
320  gradtrCy[i] = fGradtrC[1][i];//duy , dvy , dwy
321  gradtrCz[i] = fGradtrC[2][i];//duz , dvz , dwz
322  }
323 
324  STATE *efptr = &ef(0,0);
325  REAL *dphiptr = &dphi(0,0);
326  int nrowek = ek.Rows();
327  STATE *ekptr = &ek(0,0);
328  for(k=0; k<3; k++) {
329  STATE kval[3] = {(fac2*gradJx[k]+c*gradtrCx[k]),(fac2*gradJy[k]+c*gradtrCy[k]),(fac2*gradJz[k]+c*gradtrCz[k])};
330  for(i=0; i<nshape; i++) {
331  efptr[i*3+k] += -weight*(dphiptr[3*i]*kval[0]+
332  dphiptr[3*i+1]*kval[1]+
333  dphiptr[3*i+2]*kval[2]);//G = grad W
334  }
335  }
336  for(k=0; k<3; k++) {
337  for(l=0; l<3; l++) {
338  for(i=0; i<nshape; i++) {
339  for(j=0; j<nshape; j++) {
340  ekptr[i*3+k+nrowek*(j*3+l)] += weight*(
341  dphiptr[3*i]*dphiptr[3*j]*global[k][l][0] + dphiptr[3*i]*dphiptr[1+3*j]*global[k][l][1] +
342  dphiptr[3*i]*dphiptr[2+3*j]*global[k][l][2] + dphiptr[1+3*i]*dphiptr[3*j]*global[k][l][3] +
343  dphiptr[1+3*i]*dphiptr[1+3*j]*global[k][l][4] + dphiptr[1+3*i]*dphiptr[2+3*j]*global[k][l][5] +
344  dphiptr[2+3*i]*dphiptr[3*j]*global[k][l][6] + dphiptr[2+3*i]*dphiptr[1+3*j]*global[k][l][7] +
345  dphiptr[2+3*i]*dphiptr[2+3*j]*global[k][l][8] );
346  }
347  }
348  }
349  }
350 
351 #endif
352 }
353 
355 int TPZMatHyperElastic::VariableIndex(const std::string &name){
356  if(!strcmp("Displacement6",name.c_str())) return 1;
357  if(!strcmp("displacement",name.c_str())) return 2;
358  if(!strcmp("Solution",name.c_str())) return 2;
359  if(!strcmp("Derivate",name.c_str())) return 3;
360  if(!strcmp("VonMises",name.c_str())) return 4;
361  if(!strcmp("POrder",name.c_str())) return 10;
362  return TPZMaterial::VariableIndex(name);
363 }
364 
366 
367  if(var == 1) return 6;
368  if(var == 2) return 3;
369  if(var == 3) return 9;
370  if(var == 4) return 1;
371  if(var == 10) return 1;
373 }
374 
376 
377  if(var == 1) Solout.Resize(6,0.);
378  if(var == 2) Solout.Resize(3,0.);
379  if(var == 1|| var == 2) {
380  Solout[0] = Sol[0];//function
381  Solout[1] = Sol[1];//function
382  Solout[2] = Sol[2];//function
383  } else
384  if(var == 3) {
385  Solout.Resize(9);
386  int k=0;
387  for(int i=0;i<3;i++) {
388  Solout[k++] = DSol(0,i);//derivate
389  Solout[k++] = DSol(1,i);//derivate
390  Solout[k++] = DSol(2,i);//derivate
391  }
392  }
393  else if(var == 10) {
394  Solout.Resize(1);
395  Solout[0] = 1.;
396  }
397  else if(var == 4) {
398 
399  TPZFMatrix<STATE> F(DSol),Ft(3,3,0.0),I(3,3,0.),S(3,3,0.);
400  F(0,0)+=1;
401  F(1,1)+=1;
402  F(2,2)+=1;//a diagonal e' sempre a mesma
403  F.Transpose(&Ft);
404  TPZFMatrix<STATE> B = F*Ft;
405 
406  STATE ux,uy,uz,vx,vy,vz,wx,wy,wz;
407  ux = DSol(0,0);
408  uy = DSol(1,0);
409  uz = DSol(2,0);
410 
411  vx = DSol(0,1);
412  vy = DSol(1,1);
413  vz = DSol(2,1);
414 
415  wx = DSol(0,2);
416  wy = DSol(1,2);
417  wz = DSol(2,2);
418 
419  STATE J = (ux+1.)*(vy+1.)*(wz+1.) + vx*wy*uz + wx*uy*vz - wx*(vy+1.)*uz - (ux+1.)*wy*vz - vx*uy*(wz+1.);
420 
421  I(0,0)=1.0;
422  I(1,1)=1.0;
423  I(2,2)=1.0;
424  TPZFMatrix<STATE> sigmaF = (STATE(fNu/J)*B+STATE((fLambda*STATE(0.5)*(J*J-1.0)-fNu)/J)*I);
425  STATE trsigma = sigmaF(0,0)+ sigmaF(1,1)+ sigmaF(2,2);
426  S = sigmaF - STATE(trsigma/3.0)*I;
427  int i,j;
428  STATE J2=0.0;
429  for(i=0;i<3;i++) for(j=0;j<3;j++) J2 += S(i,j)* S(i,j);
430  Solout[0] = sqrt(3.0*J2);
431 
432 
433  }
434  else TPZMaterial::Solution(Sol,DSol,axes,var,Solout);
435 }
436 
438  //Flux(TPZVec<REAL> &x, TPZVec<REAL> &Sol, TPZFMatrix<REAL> &DSol, TPZFMatrix<REAL> &axes, TPZVec<REAL> &flux)
439 }
440 
442  TPZFMatrix<STATE> &dudx, TPZFMatrix<REAL> &axes, TPZVec<STATE> &/*flux*/,
443  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values) {
444 
445  //TPZVec<REAL> sol(1),dsol(3);
446  TPZManVector<STATE> sol(3),dsol(9);
447  Solution(u,dudx,axes,2,sol);
448  Solution(u,dudx,axes,3,dsol);
449  //values[1] : erro em norma L2
450  values[1] = pow(sol[0] - u_exact[0],(STATE)2.0);
451  values[1] += pow(sol[1] - u_exact[1],(STATE)2.0);
452  values[1] += pow(sol[2] - u_exact[2],(STATE)2.0);
453  //values[2] : erro em semi norma H1
454  int k=0;
455  values[2] = 0.;
456  for(int i=0;i<3;i++) {
457  values[2] += pow(dsol[k++] - du_exact(0,i),(STATE)2.0);
458  values[2] += pow(dsol[k++] - du_exact(1,i),(STATE)2.0);
459  values[2] += pow(dsol[k++] - du_exact(2,i),(STATE)2.0);
460  }
461  //values[0] : erro em norma H1 <=> norma Energia
462  values[0] = values[1]+values[2];
463 }
464 
466  REAL weight,
467  TPZFMatrix<STATE> &ek,
468  TPZFMatrix<STATE> &ef,
469  TPZBndCond &bc){
470  TPZFMatrix<REAL> &phi = data.phi;
471 
472  const STATE BIGNUMBER = 1.e12;
473 
474  const int phr = phi.Rows();
475  int in,jn,idf,jdf;
476  STATE v2[3];
477  v2[0] = bc.Val2()(0,0);
478  v2[1] = bc.Val2()(1,0);
479  v2[2] = bc.Val2()(2,0);
480  TPZFMatrix<STATE> &v1 = bc.Val1();
481 
482  int numbersol = data.sol.size();
483  if (numbersol != 1) {
484  DebugStop();
485  }
486 
487  switch (bc.Type()) {
488  case 0: // Dirichlet condition
489  for(in = 0 ; in < phr; in++) {
490  ef(3*in+0,0) += BIGNUMBER * v2[0] * phi(in,0) * weight;
491  ef(3*in+1,0) += BIGNUMBER * v2[1] * phi(in,0) * weight;
492  ef(3*in+2,0) += BIGNUMBER * v2[2] * phi(in,0) * weight;
493 
494  for (jn = 0 ; jn < phr; jn++) {
495  ek(3*in+0,3*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
496  ek(3*in+1,3*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
497  ek(3*in+2,3*jn+2) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight;
498  }//jn
499  }//in
500  break;
501 
502  case 1: // Neumann condition
503  for(in = 0 ; in < phi.Rows(); in++) {
504  ef(3*in+0,0) += v2[0] * phi(in,0) * weight;
505  ef(3*in+1,0) += v2[1] * phi(in,0) * weight;
506  ef(3*in+2,0) += v2[2] * phi(in,0) * weight;
507  }//in
508  break;
509  case 2: // Mixed condition
510  for(in = 0 ; in < phi.Rows(); in++) {
511  ef(3*in+0,0) += v2[0] * phi(in,0) * weight;
512  ef(3*in+1,0) += v2[1] * phi(in,0) * weight;
513  ef(3*in+2,0) += v2[2] * phi(in,0) * weight;
514  for(jn=0; jn<phi.Rows(); jn++)
515  {
516  for(idf=0; idf<3; idf++) for(jdf=0; jdf<3; jdf++)
517  {
518  ek(3*in+idf,3*jn+jdf) += bc.Val1()(idf,jdf);
519  }
520  }
521  }//in
522  break;
523  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
524  for(in = 0 ; in < phr; in++) {
525  ef(3*in+0,0) += BIGNUMBER * (0. - data.sol[0][0]) * v2[0] * phi(in,0) * weight;
526  ef(3*in+1,0) += BIGNUMBER * (0. - data.sol[0][1]) * v2[1] * phi(in,0) * weight;
527  ef(3*in+2,0) += BIGNUMBER * (0. - data.sol[0][2]) * v2[2] * phi(in,0) * weight;
528  for (jn = 0 ; jn < phr; jn++) {
529  ek(3*in+0,3*jn+0) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * v2[0];
530  ek(3*in+1,3*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * v2[1];
531  ek(3*in+2,3*jn+2) += BIGNUMBER * phi(in,0) * phi(jn,0) * weight * v2[2];
532  }//jn
533  }//in
534  break;
535 
536  case 4: // stressField Neumann condition
537  for(in = 0; in < 3; in ++)
538  v2[in] = - ( v1(in,0) * data.normal[0] +
539  v1(in,1) * data.normal[1] +
540  v1(in,2) * data.normal[2] );
541  // The normal vector points towards the neighbour. The negative sign is there to
542  // reflect the outward normal vector.
543  for(in = 0 ; in < phi.Rows(); in++) {
544  ef(3*in+0,0) += v2[0] * phi(in,0) * weight;
545  ef(3*in+1,0) += v2[1] * phi(in,0) * weight;
546  ef(3*in+2,0) += v2[2] * phi(in,0) * weight;
547  //cout << "normal:" << data.normal[0] << ' ' << data.normal[1] << ' ' << data.normal[2] << endl;
548  //cout << "val2: " << v2[0] << ' ' << v2[1] << ' ' << v2[2] << endl;
549  }
550  break;
551  default:
552  PZError << "TPZElastitity3D::ContributeBC error - Wrong boundary condition type" << std::endl;
553  }//switch
554 
555 }//method
556 
557 //
558 //void TPZMatHyperElastic::ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef, TPZBndCond &bc){
559 //
560 // TPZFMatrix<REAL> &phi = data.phi;
561 // int numbersol = data.sol.size();
562 // if (numbersol != 1) {
563 // DebugStop();
564 // }
565 // TPZVec<STATE> &sol=data.sol[0];
566 //
567 // if(bc.Material() != this){
568 // PZError << "TPZMatHyperElastic.ContributeBC : this material don't exists \n";
569 // }
570 //
571 // if(bc.Type() < 0 && bc.Type() > 2){
572 // PZError << "ContributeBC.aplybc, unknown boundary condition type : "<<bc.Type() << endl;
573 // }
574 //
575 // int ndof = NStateVariables();
576 // int nnod = ek.Rows()/ndof;
577 // int r = ndof;
578 //
579 // int idf,jdf,in,jn;
580 // switch(bc.Type()){
581 // case 0:
582 // for(in=0 ; in<nnod ; ++in){
583 // for(idf = 0;idf<r;idf++) {
584 // (ef)(in*r+idf,0) += gBigNumber*phi(in,0)*(bc.Val2()(idf,0)-sol[idf])*weight;
585 // }
586 // for(jn=0 ; jn<nnod ; ++jn) {
587 // for(idf = 0;idf<r;idf++) {
588 // ek(in*r+idf,jn*r+idf) += gBigNumber*phi(in,0)*phi(jn,0)*weight;
589 // }
590 // }
591 // }
592 // break;
593 //
594 // case 1:
595 // for(in=0 ; in<nnod ; ++in){
596 // for(idf = 0;idf<r;idf++) {
597 // //(ef)(in*r+idf,0) += weight*phi(in,0)*(bc.Val2()(idf,0)-sol[idf]);
598 // (ef)(in*r+idf,0) += weight*phi(in,0)*(bc.Val2()(idf,0));
599 // }
600 // }
601 // break;
602 //
603 // case 2:
604 // for(in=0 ; in<nnod ; ++in){
605 // for(idf = 0;idf<r;idf++) {
606 // for (jdf=0; jdf<r; jdf++){
607 // (ef)(in*r+idf,0) += phi(in,0)*bc.Val1()(idf,jdf)*(bc.Val2()(jdf,0)-sol[jdf])*weight;
608 // }
609 // for(jn=0 ; jn<nnod ; ++jn) {
610 // for(idf = 0;idf<r;idf++) {
611 // for(jdf = 0;jdf<r;jdf++) {
612 // ek(in*r+idf,jn*r+jdf) += bc.Val1()(idf,jdf)*phi(in,0)*phi(jn,0)*weight;
613 // }
614 // }
615 // }
616 // }
617 //
618 // }
619 //
620 // }//fim switch
621 //}
622 
623 #ifdef _AUTODIFF
624 
626 inline int ith(const int i, const int j)
627 {
628  return i*3+j;
629 }
630 
631 void TPZMatHyperElastic::ContributeEnergy(TPZVec<REAL> &x,
633  FADFADREAL &U, REAL weight)
634 {
635  DebugStop();
636  /*
637  FADFADREAL J, TrC; // J = det(F); TrC = Trace(C)
638  FADFADREAL DiagF0(dsol[ 0]);
639  DiagF0.val().val() += 1.;// element [0][0]
640 
641  FADFADREAL DiagF1(dsol[ 3+1]);
642  DiagF1.val().val() += 1.;// element [0][0]
643 
644  FADFADREAL DiagF2(dsol[2*3+2]);
645  DiagF2.val().val() += 1.;// element [0][0]
646  TrC = DiagF0*DiagF0;
647  TrC += dsol[ith(1,0)] * dsol[ith(1,0)];
648  TrC += dsol[ith(2,0)] * dsol[ith(2,0)];
649 
650  TrC += dsol[ith(0,1)] * dsol[ith(0,1)];
651  TrC += DiagF1*DiagF1;
652  TrC += dsol[ith(2,1)] * dsol[ith(2,1)];
653 
654  TrC += dsol[ith(0,2)] * dsol[ith(0,2)];
655  TrC += dsol[ith(1,2)] * dsol[ith(1,2)];
656  TrC += DiagF2*DiagF2;
657  // cout << "TrC\n" << TrC << endl;
658  J = DiagF0 * DiagF1 * DiagF2;
659  J += dsol[ith(0,1)] * dsol[ith(1,2)] * dsol[ith(2,0)];
660  J += dsol[ith(0,2)] * dsol[ith(1,0)] * dsol[ith(2,1)];
661  J -= dsol[ith(0,2)] * DiagF1 * dsol[ith(2,0)];
662  J -= dsol[ith(0,1)] * dsol[ith(1,0)] * DiagF2;
663  J -= DiagF0 * dsol[ith(1,2)] * dsol[ith(2,1)]; // J = det(F)
664 
665  // cout << "J\n " << J << endl;
666  U += (J*J - FADREAL(1.)) * FADREAL(weight*fLambda/4.);
667  U -= log( J ) * FADREAL(weight*(fLambda/2.+fNu));
668  U += (TrC - FADREAL(3.)) * FADREAL(weight*fNu/2.);
669  */
670 }
671 
672 void TPZMatHyperElastic::ContributeBCEnergy(TPZVec<REAL> & x,
673  TPZVec<FADFADREAL> & sol, FADFADREAL &U,
674  REAL weight, TPZBndCond &bc)
675 {
676  DebugStop();
677  /*
678  if(bc.Material() != this){
679  PZError << "TPZMatHyperElastic.ContributeBC : this material doesn't exist \n";
680  }
681 
682  if(bc.Type() < 0 && bc.Type() > 2){
683  PZError << "ContributeBC.aplybc, unknown boundary condition type : "<<bc.Type() << endl;
684  }
685 
686  TPZVec<FADFADREAL> solMinBC(3), BCsolMinBC(3);
687  solMinBC[0] = sol[0] - FADREAL(bc.Val2()(0,0));
688  solMinBC[1] = sol[1] - FADREAL(bc.Val2()(1,0));
689  solMinBC[2] = sol[2] - FADREAL(bc.Val2()(2,0));
690 
691  switch(bc.Type()){
692  case 0:// Dirichlet condition
693  // U += 1/2* Big * weight * Integral((u - u0)^2 dOmega)
694  U += ( (solMinBC[0] * solMinBC[0]) +
695  (solMinBC[1] * solMinBC[1]) +
696  (solMinBC[2] * solMinBC[2]) ) *
697  FADREAL(gBigNumber * weight / 2.);
698  break;
699  case 1:// Neumann condition
700  // U -= weight * Integral([g].u dOmega)
701  U -= ( sol[0] * FADREAL(bc.Val2()(0,0) ) +
702  sol[1] * FADREAL(bc.Val2()(1,0) ) +
703  sol[2] * FADREAL(bc.Val2()(2,0) ) ) *
704  FADREAL(weight);
705  break;
706  case 2:// condi�o mista
707  // U += 1/2 * weight * Integral(<(u-u0), [g].(u-u0)> dOmega)
708  BCsolMinBC[0] = solMinBC[0] * FADREAL(bc.Val1()(0,0)) +
709  solMinBC[1] * FADREAL(bc.Val1()(0,1)) +
710  solMinBC[2] * FADREAL(bc.Val1()(0,2));
711  BCsolMinBC[1] = solMinBC[0] * FADREAL(bc.Val1()(1,0)) +
712  solMinBC[1] * FADREAL(bc.Val1()(1,1)) +
713  solMinBC[2] * FADREAL(bc.Val1()(1,2));
714  BCsolMinBC[2] = solMinBC[0] * FADREAL(bc.Val1()(2,0)) +
715  solMinBC[1] * FADREAL(bc.Val1()(2,1)) +
716  solMinBC[2] * FADREAL(bc.Val1()(2,2));
717  U += ( solMinBC[0] * BCsolMinBC[0] +
718  solMinBC[1] * BCsolMinBC[1] +
719  solMinBC[2] * BCsolMinBC[2] ) *
720  FADREAL(weight / 2.);
721  break;
722  }
723  */
724 }
725 
726 void TPZMatHyperElastic::ComputeEnergy(STATE lambda, STATE mu, TPZFMatrix<STATE> &dsol, TFad<9,TFad<9,STATE> > &energy) {
727  DebugStop();
728  /*
729  TFad<9,TFad<9,STATE> > tensor[3][3],J,TrC;
730  int i,j;
731  for(i=0; i<3; i++) {
732  for(j=0; j<3; j++) {
733  tensor[i][j].val().val() = dsol(j,i);
734  tensor[i][j].fastAccessDx(j*3+i).val() = 1.;
735  tensor[i][j].val().fastAccessDx(j*3+i) = 1.;
736  }
737  tensor[i][i].val().val() += 1.;
738  }
739  TrC = tensor[0][0]*tensor[0][0]+tensor[0][1]*tensor[0][1]+tensor[0][2]*tensor[0][2]+tensor[1][0]*tensor[1][0]+tensor[1][1]*tensor[1][1]+tensor[1][2]*tensor[1][2]+
740  tensor[2][0]*tensor[2][0]+tensor[2][1]*tensor[2][1]+tensor[2][2]*tensor[2][2];
741 
742  J = tensor[0][0] * tensor[1][1] * tensor[2][2] +
743  tensor[0][1] * tensor[1][2] * tensor[2][0] +
744  tensor[0][2] * tensor[1][0] * tensor[2][1] -
745  tensor[0][2] * tensor[1][1] * tensor[2][0] -
746  tensor[0][1] * tensor[1][0] * tensor[2][2] -
747  tensor[0][0] * tensor[1][2] * tensor[2][1]; // J = det(F)
748 
749  energy = (J*J - TFad<9,STATE>(1.)) * TFad<9,STATE>(lambda/4.) -
750  log( J ) * TFad<9,STATE>((lambda/2.+mu)) +
751  (TrC - TFad<9,STATE>(3.)) * TFad<9,STATE>(mu/2.);
752  */
753 }
754 
755 #endif
756 
758  return Hash("TPZMatHyperElastic") ^ TPZMaterial::ClassId() << 1;
759 }
virtual void Solution(TPZMaterialData &data, int var, TPZVec< STATE > &Solout)
Returns the solution associated with the var index based on the finite element approximation.
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZManVector< REAL, 3 > normal
normal to the element at the integration point
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
clarg::argBool bc("-bc", "binary checkpoints", false)
Implements a hyper elasticity material.
int Type()
Definition: pzbndcond.h:249
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.
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
Defines PZError.
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
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...
virtual int VariableIndex(const std::string &name) override
TPZFMatrix< STATE > & Val2(int loadcase=0)
Definition: pzbndcond.h:255
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
std::string Name() override
Returns the name of the material.
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
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...
This abstract class defines the behaviour which each derived class needs to implement.
Definition: TPZMaterial.h:39
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
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
Contains TPZMatrixclass which implements full matrix (using column major representation).
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
Contains the TPZMatHyperElastic class which implements a hyper elasticity material.
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
Definition: tfad.h:64
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
const T & d(int i) const
Definition: tfad.h:96
virtual void Print(std::ostream &out) override
Prints out the data associated with the material.
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
string res
Definition: test.py:151
TPZMatHyperElastic(int nummat, STATE e, STATE mu, STATE nu=-1., STATE lambda=-1., STATE coef1=-1., STATE coef2=-1., STATE coef3=-1.)
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
Contains TPZMatrix<TVar>class, root matrix class.
virtual void Solution(TPZVec< STATE > &Sol, TPZFMatrix< STATE > &DSol, TPZFMatrix< REAL > &axes, int var, TPZVec< STATE > &Solout) override
int ClassId() const override
Unique identifier for serialization purposes.
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
static int idf[6]
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
virtual int ClassId() const override
Define the class id associated with the class.
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var.
def values
Definition: rdt.py:119
TPZSolVec sol
vector of the solutions at the integration point
const T & val() const
Definition: tfad.h:89
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...
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.