NeoPZ
pzelasAXImat.cpp
Go to the documentation of this file.
1 
6 #include "pzelasAXImat.h"
7 #include "pzelmat.h"
8 #include "pzbndcond.h"
9 #include "pzmatrix.h"
10 #include "pzfmatrix.h"
11 #include "pzerror.h"
12 #include <cmath>
13 #include <math.h>
14 #include "pzaxestools.h"
15 #include "pzlog.h"
16 
17 #ifdef LOG4CXX
18 static LoggerPtr logger(Logger::getLogger("pz.material.axisymetric"));
19 static LoggerPtr logdata(Logger::getLogger("pz.material.axisymetric.data"));
20 #endif
21 
22 #include <fstream>
23 using namespace std;
24 
25 TPZElasticityAxiMaterial::TPZElasticityAxiMaterial() :
26 TPZRegisterClassId(&TPZElasticityAxiMaterial::ClassId),
27 TPZDiscontinuousGalerkin(0),
28 fIntegral(0.),
29 fAlpha(1.e-5),
30 f_AxisR(3,0.),
31 f_AxisZ(3,0.),
32 f_Origin(3,0.),
33 fTemperatureFunction(0) {
34  f_AxisZ[1] = 1.;
35  f_AxisR[0] = 1.;
36  fDelTemperature = 0.;
37  fE = -1.; // Young modulus
38  fnu = -1.; // poisson coefficient
39  ff[0] = 0.; // X component of the body force
40  ff[1] = 0.; // Y component of the body force
41  ff[2] = 0.; // Z component of the body force - not used for this class
42  fEover1MinNu2 = -1.; //G = E/2(1-nu);
43  fEover21PlusNu = -1.;//E/(1-nu)
44 
45  f_c = 0.;
46  f_phi = 0.;
47  fSymmetric = 1.;
48  fPenaltyConstant = 1.;
49 }
50 
51 TPZElasticityAxiMaterial::TPZElasticityAxiMaterial(int num, REAL E, REAL nu, REAL fx, REAL fy) :
52 TPZRegisterClassId(&TPZElasticityAxiMaterial::ClassId),
53 TPZDiscontinuousGalerkin(num), fIntegral(0.), fAlpha(1.e-5), fDelTemperature(0.), f_AxisR(3,0.), f_AxisZ(3,0.),f_Origin(3,0.), fTemperatureFunction(0)
54 {
55 
56  f_AxisZ[1] = 1.;
57  f_AxisR[0] = 1.;
58 
59  fE = E; // Young modulus
60  fnu = nu; // poisson coefficient
61  ff[0] = fx; // X component of the body force
62  ff[1] = fy; // Y component of the body force
63  ff[2] = 0.; // Z component of the body force - not used for this class
64  fEover1MinNu2 = E/(1-fnu*fnu); //G = E/2(1-nu);
65  fEover21PlusNu = E/(2.*(1+fnu));//E/(1-nu)
66  f_c = 0.;
67  f_phi = 0.;
68  fSymmetric = 1.;
69  fPenaltyConstant = 1.;
70 }
71 
72 //--------------------------------------------------------------------------------------------------------------------------------------
73 TPZElasticityAxiMaterial::TPZElasticityAxiMaterial(int num, REAL E, REAL nu, REAL fx, REAL fy, REAL coefTheta, REAL coefAlpha) :
74 TPZRegisterClassId(&TPZElasticityAxiMaterial::ClassId),
75 TPZDiscontinuousGalerkin(num), fIntegral(0.), fAlpha(1.e-5), fDelTemperature(0.), f_AxisR(3,0.), f_AxisZ(3,0.),f_Origin(3,0.),
76 fTemperatureFunction(0)
77 {
78 
79  f_AxisZ[1] = 1.;
80  f_AxisR[0] = 1.;
81  fE = E; // Young modulus
82  fnu = nu; // poisson coefficient
83  ff[0] = fx; // X component of the body force
84  ff[1] = fy; // Y component of the body force
85  ff[2] = 0.; // Z component of the body force - not used for this class
86  fEover1MinNu2 = E/(1-fnu*fnu); //G = E/2(1-nu);
87  fEover21PlusNu = E/(2.*(1+fnu));//E/(1-nu)
88  f_c = 0.;
89  f_phi = 0.;
90  fSymmetric = coefTheta;
91  fPenaltyConstant = coefAlpha;
92 }
93 
94 TPZElasticityAxiMaterial::TPZElasticityAxiMaterial(const TPZElasticityAxiMaterial &copy) :
95 TPZRegisterClassId(&TPZElasticityAxiMaterial::ClassId),
96 TPZDiscontinuousGalerkin(copy), fIntegral(copy.fIntegral),f_phi(copy.f_phi),f_c(copy.f_c), fE(copy.fE),
97 fnu(copy.fnu), fAlpha(copy.fAlpha), fDelTemperature(copy.fDelTemperature), fEover21PlusNu(copy.fEover21PlusNu),
98 fEover1MinNu2(copy.fEover1MinNu2),f_AxisR(copy.f_AxisR),f_AxisZ(copy.f_AxisZ),
99 f_Origin(copy.f_Origin),fSymmetric(copy.fSymmetric),fPenaltyConstant(copy.fPenaltyConstant),fTemperatureFunction(copy.fTemperatureFunction)
100 {
101  ff[0] = copy.ff[0];
102  ff[1] = copy.ff[1];
103  ff[2] = copy.ff[2];
104 }
105 
106 
107 void TPZElasticityAxiMaterial::SetOrigin(TPZManVector<REAL> &Orig, TPZManVector<REAL> &AxisZ, TPZManVector<REAL> &AxisR)
108 {
109  if(Orig.size() == 3 && AxisZ.size() == 3 && AxisR.size() == 3)
110  {
111  TPZFNMatrix<6> Vecs(3,2,0.), VecsOrt(3,2,0.), JacVecsOrth(2,2,0.);
112  for(int i = 0; i < 3; i++)
113  {
114  Vecs.PutVal(i,0,AxisR[i]);
115  Vecs.PutVal(i,1,AxisZ[i]);
116  }
117  Vecs.GramSchmidt(VecsOrt,JacVecsOrth);
118 
119  for(int j = 0; j < 3; j++)
120  {
121  AxisR[j] = VecsOrt.GetVal(j,0);
122  AxisZ[j] = VecsOrt.GetVal(j,1);
123  }
124  f_Origin = Orig;
125  f_AxisZ = AxisZ;
126  f_AxisR = AxisR;
127  }
128  else
129  {
130  cout << "Invalid Origin and/or Axis vector on TPZElasticityAxiMaterial()!\n";
131  DebugStop();
132  }
133 }
134 
135 REAL TPZElasticityAxiMaterial::ComputeR(TPZVec<REAL> &x)
136 {
137  return (x[0] - f_Origin[0])*f_AxisR[0] + (x[1] - f_Origin[1])*f_AxisR[1] + (x[2] - f_Origin[2])*f_AxisR[2];
138 }
139 
140 TPZManVector<REAL> TPZElasticityAxiMaterial::GetAxisR()
141 {
142  return f_AxisR;
143 }
144 
145 TPZManVector<REAL> TPZElasticityAxiMaterial::GetAxisZ()
146 {
147  return f_AxisZ;
148 }
149 
150 TPZManVector<REAL> TPZElasticityAxiMaterial::GetOrigin()
151 {
152  return f_Origin;
153 }
154 
155 TPZElasticityAxiMaterial::~TPZElasticityAxiMaterial() {
156 }
157 
158 int TPZElasticityAxiMaterial::NStateVariables() const {
159  return 2;
160 }
161 
162 void TPZElasticityAxiMaterial::Print(std::ostream &out) {
163  TPZMaterial::Print(out);
164  out << "name of material : " << Name() << "\n";
165  out << "properties : \n";
166  out << "\tE = " << fE << endl;
167  out << "\tnu = " << fnu << endl;
168  out << "\tF = " << ff[0] << ' ' << ff[1] << endl;
169 }
170 
171 void TPZElasticityAxiMaterial::Contribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
172 {
173  TPZFMatrix<REAL> &dphi = data.dphix;
174  TPZFMatrix<REAL> &phi = data.phi;
175  TPZFMatrix<REAL> &axes = data.axes;
176 
177  int64_t phc,phr,dphc,dphr,efr,efc,ekr,ekc;
178  phc = phi.Cols();
179  phr = phi.Rows();
180  dphc = dphi.Cols();
181  dphr = dphi.Rows();
182  efr = ef.Rows();
183  efc = ef.Cols();
184  ekr = ek.Rows();
185  ekc = ek.Cols();
186  if(phc != 1 || dphr != 2 || phr != dphc || ekr != phr*2 || ekc != phr*2 || efr != phr*2 || efc != 1)
187  {
188  PZError << "\nTPZElasticityMaterial.contr, inconsistent input data : \n" <<
189  "phi.Cols() = " << phi.Cols() << " dphi.Cols() = " << dphi.Cols() <<
190  " phi.Rows = " << phi.Rows() << " dphi.Rows = " <<
191  dphi.Rows() << "\nek.Rows() = " << ek.Rows() << " ek.Cols() = "
192  << ek.Cols() <<
193  "\nef.Rows() = " << ef.Rows() << " ef.Cols() = "
194  << ef.Cols() << "\n";
195  return;
196  }
197  if(fForcingFunction)
198  {
199  TPZManVector<STATE> res(3);
200  fForcingFunction->Execute(data.x,res);
201  ff[0] = res[0];
202  ff[1] = res[1];
203  ff[2] = res[2];
204  }
205 
206  // R = Dot[{data.x - origin},{AxisR}] ***because AxisR is already normalized!
207  REAL R = (data.x[0] - f_Origin[0])*f_AxisR[0] + (data.x[1] - f_Origin[1])*f_AxisR[1] + (data.x[2] - f_Origin[2])*f_AxisR[2];
208  REAL Z = (data.x[0] - f_Origin[0])*f_AxisZ[0] + (data.x[1] - f_Origin[1])*f_AxisZ[1] + (data.x[2] - f_Origin[2])*f_AxisZ[2];
209 
210  int s = (R > 0)? 1:-1;
211  R = fabs(R);
212  if(R < 1.e-10) R = 1.e-10;
213 
214  REAL DelTemp(fDelTemperature);
215  if (fTemperatureFunction) {
216  TPZManVector<REAL,2> RZ(2);
217  RZ[0] = R;
218  RZ[1] = Z;
219  fTemperatureFunction(RZ,DelTemp);
220  }
224  // REAL nu1 = 1. - fnu;//(1-nu)
225  // REAL nu2 = (1.-2.*fnu)/2.;
226  // REAL F = fE/((1.+fnu)*(1.-2.*fnu));
227 
228  TPZFNMatrix<4> dphiRZi(2,1), dphiRZj(2,1);
229 
230  REAL axis0DOTr = 0., axis1DOTr = 0., axis0DOTz = 0., axis1DOTz = 0.;
231  for(int pos = 0; pos < 3; pos++)
232  {
233  axis0DOTr += axes.GetVal(0,pos) * f_AxisR[pos] * s;
234  axis1DOTr += axes.GetVal(1,pos) * f_AxisR[pos] * s;
235  axis0DOTz += axes.GetVal(0,pos) * f_AxisZ[pos];
236  axis1DOTz += axes.GetVal(1,pos) * f_AxisZ[pos];
237  }
238 
239  REAL R2PI = 2. * M_PI * R;
240 
241  REAL lambda = -((fE*fnu)/((1. + fnu)*(2.*fnu-1.)));
242  REAL mi = fE/(2.*(1. + fnu));
243 
244  REAL epsT = DelTemp*fAlpha;
245  REAL TensorThermico = (3.*lambda+2.*mi)*epsT;
246 
247  //criado para resolver o problema de Girkmann
248  fIntegral += R2PI*weight;
249 
250  for( int64_t in = 0; in < phr; in++ )
251  {
252 
253  //dphi_i/dr = dphi_i/axis0 <axes0,f_AxisR> + dphi_i/axis1 <axes1,f_AxisR>
254  dphiRZi.PutVal(0,0, dphi.GetVal(0,in)*axis0DOTr + dphi.GetVal(1,in)*axis1DOTr );
255 
256  //dphi_i/dz = dphi_i/axis0 <axes0,f_AxisZ> + dphi_i/axis1 <axes1,f_AxisZ>
257  dphiRZi.PutVal(1,0, dphi.GetVal(0,in)*axis0DOTz + dphi.GetVal(1,in)*axis1DOTz );
258 
259  ef(2*in, 0) += weight * R2PI * (ff[0] * phi(in,0) + TensorThermico*(phi(in,0)/R + dphiRZi(0,0))); // direcao x
260  ef(2*in+1, 0) += weight * R2PI * (ff[1] * phi(in,0) + TensorThermico*dphiRZi(1,0)); // direcao y
261 
262  for( int64_t jn = 0; jn < phr; jn++ )
263  {
264 
265  //dphi_j/dr = dphi_j/axis0 <axes0,f_AxisR> + dphi_j/axis1 <axes1,f_AxisR>
266  dphiRZj.PutVal(0,0, dphi.GetVal(0,jn)*axis0DOTr + dphi.GetVal(1,jn)*axis1DOTr );
267 
268  //dphi_j/dz = dphi_j/axis0 <axes0,f_AxisZ> + dphi_j/axis1 <axes1,f_AxisZ>
269  dphiRZj.PutVal(1,0, dphi.GetVal(0,jn)*axis0DOTz + dphi.GetVal(1,jn)*axis1DOTz );
270 
271  REAL term00 = dphiRZi(0,0) * (lambda + 2.*mi) * dphiRZj(0,0) +
272  dphiRZi(1,0) * mi * dphiRZj(1,0) +
273  phi(in,0) * lambda/R * dphiRZj(0,0) +
274  dphiRZi(0,0) * lambda/R * phi(jn,0) +
275  phi(in,0) * (lambda+2.*mi)/(R*R) * phi(jn,0);
276 
277  REAL term01 = dphiRZi(1,0) * mi * dphiRZj(0,0) +
278  dphiRZi(0,0) * lambda * dphiRZj(1,0) +
279  phi(in,0) * lambda/R * dphiRZj(1,0);
280 
281  REAL term10 = dphiRZi(1,0) * lambda * dphiRZj(0,0) +
282  dphiRZi(0,0) * mi * dphiRZj(1,0) +
283  dphiRZi(1,0) * lambda/R * phi(jn,0);
284 
285  REAL term11 = dphiRZi(0,0) * mi * dphiRZj(0,0) +
286  dphiRZi(1,0) * (lambda + 2.*mi) * dphiRZj(1,0);
287 
288  ek(2*in,2*jn) += weight * R2PI * term00;
289  ek(2*in,2*jn+1) += weight * R2PI * term01;
290  ek(2*in+1,2*jn) += weight * R2PI * term10;
291  ek(2*in+1,2*jn+1) += weight * R2PI * term11;
292  }
293  }
294 
295 #ifdef LOG4CXX
296  if(logdata->isDebugEnabled())
297  {
298  std::stringstream sout;
299  ek.Print("ek = ",sout,EMathematicaInput);
300  LOGPZ_DEBUG(logdata,sout.str())
301  }
302 #endif
303 
304 }
305 
306 void TPZElasticityAxiMaterial::ContributeBC(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef,TPZBndCond &bc)
307 {
308  TPZFMatrix<REAL> &phi = data.phi;
309 
310  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
311 
312  int64_t phr = phi.Rows();
313  short in,jn;
314  REAL v2[2];
315  v2[0] = bc.Val2()(0,0);
316  v2[1] = bc.Val2()(1,0);
317 
318  // R = Dot[{data.x - origin},{AxisR}] ***because AxisR is already normalized!
319  REAL R = (data.x[0] - f_Origin[0])*f_AxisR[0] + (data.x[1] - f_Origin[1])*f_AxisR[1] + (data.x[2] - f_Origin[2])*f_AxisR[2];
320 
321  int s = (R > 0) ? 1:-1;
322  R = fabs(R);
323  double R2PI = 2. * M_PI * R;
324 
325  static REAL accum1 = 0., accum2 = 0.;
326 
327  switch (bc.Type())
328  {
329  case 0 :// Dirichlet condition
330  {
331  for(in = 0 ; in < phr; in++)
332  {
333  ef(2*in,0) += BIGNUMBER * v2[0]* // x displacement
334  phi(in,0) * R2PI * weight; // forced v2 displacement
335 
336  ef(2*in+1,0) += BIGNUMBER * v2[1]* // x displacement
337  phi(in,0) * R2PI * weight; // forced v2 displacement
338 
339  for (jn = 0 ; jn < phi.Rows(); jn++)
340  {
341  ek(2*in,2*jn) += BIGNUMBER * phi(in,0) * phi(jn,0) * R2PI * weight;
342  ek(2*in+1,2*jn+1) += BIGNUMBER * phi(in,0) * phi(jn,0) * R2PI * weight;
343  }
344  }
345  }
346  break;
347 
348  case 1 :// Neumann condition
349  {
350  for(in = 0 ; in < phi.Rows(); in++)
351  { // componentes da tracao na direcao de v2
352  ef(2*in,0) += v2[0] * phi(in,0) * R2PI * weight; // tracao em x (ou pressao)
353  ef(2*in+1,0) += v2[1] * phi(in,0) * R2PI * weight; // tracao em y (ou pressao) , nula se n� h
354  } // ou deslocamento nulo v2 = 0
355  accum1 += v2[1] * R2PI *weight;
356 #ifdef LOG4CXX
357  if (logger->isDebugEnabled())
358  {
359  std::stringstream sout;
360  sout << "Accumulated force 1 " << accum1;
361  LOGPZ_DEBUG(logger,sout.str());
362  }
363 #endif
364  }
365  break;
366 
367  case 2 :// condicao mista
368  {
369  for(in = 0 ; in < phi.Rows(); in++)
370  {
371  ef(2*in, 0) += v2[0] * phi(in, 0) * R2PI * weight; // Neumann , Sigmaij
372  ef(2*in+1, 0) += v2[1] * phi(in, 0) * R2PI * weight; // Neumann
373 
374  for (jn = 0 ; jn < phi.Rows(); jn++)
375  {
376  ek(2*in,2*jn) += bc.Val1()(0,0) * phi(in,0) * phi(jn,0) * R2PI * weight; // peso de contorno => integral de contorno
377  ek(2*in+1,2*jn) += bc.Val1()(1,0) * phi(in,0) * phi(jn,0) * R2PI * weight;
378  ek(2*in+1,2*jn+1) += bc.Val1()(1,1) * phi(in,0) * phi(jn,0) * R2PI * weight;
379  ek(2*in,2*jn+1) += bc.Val1()(0,1) * phi(in,0) * phi(jn,0) * R2PI * weight;
380  }
381  } // este caso pode reproduzir o caso 0 quando o deslocamento
382  }
383  break;
384 
385  case 3 :// Neumann condition - Normal rotacionada 90 graus no sentido horário em relacao ao vetor axes (1D)
386  {
387  REAL Nxy[2],Nrz[2];
388  Nxy[0] = data.axes(0,1);
389  Nxy[1] = -data.axes(0,0);
390  Nrz[0] = s*(Nxy[0]*f_AxisR[0] + Nxy[1]*f_AxisR[1]);
391  Nrz[1] = Nxy[0]*f_AxisZ[0] + Nxy[1]*f_AxisZ[1];
392 
393  // #ifdef LOG4CXX
394  // {
395  // std::stringstream sout;
396  // sout << "CoordX: " << data.x << " R= " << R << " Nrz = " << Nrz[0] << "," << Nrz[1] << endl;
397  // LOGPZ_DEBUG(logger,sout.str());
398  // }
399  // #endif
400 
401  for(in = 0 ; in < phi.Rows(); in++)
402  { // componentes da tracao normal ao contorno
403  ef(2*in,0) += v2[0] * Nrz[0] * phi(in,0) * R2PI * weight; // tracao em x (ou pressao)
404  ef(2*in+1,0) += v2[0] * Nrz[1] * phi(in,0) * R2PI * weight; // tracao em y (ou pressao) , nula se n� h
405  } // ou deslocamento nulo v2 = 0
406  accum2 += v2[0] * Nrz[1] * R2PI *weight;
407 #ifdef LOG4CXX
408  if (logger->isDebugEnabled())
409  {
410  std::stringstream sout;
411  sout << "Accumulated force 2 " << accum2;
412  LOGPZ_DEBUG(logger,sout.str());
413  }
414 #endif
415  }
416  break;
417  } // �nulo introduzindo o BIGNUMBER pelos valores da condicao
418 } // 1 Val1 : a leitura �00 01 10 11
419 
420 //---------------------------- --------------------------------------
421 void TPZElasticityAxiMaterial::ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright,
422  REAL weight,
423  TPZFMatrix<STATE> &ek, TPZFMatrix<STATE> &ef){
424 
425  TPZFMatrix<REAL> &dphiL = dataleft.dphix;
426  TPZFMatrix<REAL> &dphiR = dataright.dphix;
427  TPZFMatrix<REAL> &phiL = dataleft.phi;
428  TPZFMatrix<REAL> &phiR = dataright.phi;
429  TPZFMatrix<REAL> &axesL = dataleft.axes;
430  TPZFMatrix<REAL> &axesR = dataright.axes;
431 
432  TPZManVector<REAL,3> &normal = data.normal;
433 
434  int &LeftPOrder=dataleft.p;
435  int &RightPOrder=dataright.p;
436 
437  REAL &faceSize=data.HSize;
438 
439 #ifdef LOG4CXX
440  if (logger->isDebugEnabled())
441  {
442  std::stringstream sout;
443  dataleft.phi.Print("phil",sout);
444  dataright.phi.Print("phir",sout);
445  LOGPZ_DEBUG(logger,sout.str())
446  }
447 #endif
448 
449 #ifdef LOG4CXX
450  if(logdata->isDebugEnabled())
451  {
452  std::stringstream sout;
453  sout << "Origin = { " << f_Origin << "};" << std::endl;
454  sout << "AxisR = { " << f_AxisR << "};" << std::endl;
455  sout << "AxisZ = { " << f_AxisZ << "};" << std::endl;
456  data.PrintMathematica(sout);
457  LOGPZ_DEBUG(logdata,sout.str())
458  }
459 #endif
460 
461  int64_t nrowl = phiL.Rows();
462  int64_t nrowr = phiR.Rows();
463  int64_t il,jl,ir,jr;
464 
466  REAL R = (data.x[0] - f_Origin[0])*f_AxisR[0] + (data.x[1] - f_Origin[1])*f_AxisR[1] + (data.x[2] - f_Origin[2])*f_AxisR[2];
467  int s = (R > 0)? 1:-1;
468  R = fabs(R);
469 
470 
471  TPZFNMatrix<16> dphiRZiL(2,1), dphiRZjL(2,1), dphiRZiR(2,1), dphiRZjR(2,1);
472 
473  double axis0DOTrL = 0., axis1DOTrL = 0., axis0DOTzL = 0., axis1DOTzL = 0.;
474  double axis0DOTrR = 0., axis1DOTrR = 0., axis0DOTzR = 0., axis1DOTzR = 0.;
475  for(int pos = 0; pos < 3; pos++)
476  {
477  axis0DOTrL += axesL.GetVal(0,pos) * f_AxisR[pos] * s;
478  axis1DOTrL += axesL.GetVal(1,pos) * f_AxisR[pos] * s;
479  axis0DOTzL += axesL.GetVal(0,pos) * f_AxisZ[pos];
480  axis1DOTzL += axesL.GetVal(1,pos) * f_AxisZ[pos];
481  axis0DOTrR += axesR.GetVal(0,pos) * f_AxisR[pos] * s;
482  axis1DOTrR += axesR.GetVal(1,pos) * f_AxisR[pos] * s;
483  axis0DOTzR += axesR.GetVal(0,pos) * f_AxisZ[pos];
484  axis1DOTzR += axesR.GetVal(1,pos) * f_AxisZ[pos];
485  }
486 
487  REAL R2PI = 2.0 * M_PI * R;
488 
489  REAL symmetry = fSymmetric; //thermo of symmetry ( -1: method symmetric; 1: method not symmetric)
490  REAL penalty = fPenaltyConstant*(0.5 * (LeftPOrder*LeftPOrder + RightPOrder*RightPOrder))/faceSize;
491 
492  REAL beta;
493  if (symmetry==1.0) {
494  beta=1.;
495  }else {
496  beta=1.0;
497  }
498 
499  double lambda = -((fE*fnu)/((1. + fnu)*(2.*fnu-1.)));
500  double mu = fE/(2.*(1. + fnu));
501 
502  int64_t numbersol = dataleft.dsol.size();
503  if (numbersol != 1) {
504  DebugStop();
505  }
506 
507 #ifdef LOG4CXX
508  if (logdata->isDebugEnabled())
509  {
510  std::stringstream sout;
511  sout.precision(16);
512  sout << "Penalty = " << fPenaltyConstant << ";" << std::endl;
513  sout << "R = " << R << ";" << std::endl;
514  sout << "fE = " << fE << ";" << std::endl;
515  sout << "fnu = " << fnu << ";" << std::endl;
516  sout << "Lambda = " << lambda << ";" << std::endl;
517  sout << "Mu = " << mu << ";" << std::endl;
518  sout << "weight = " << weight << ";" << std::endl;
519  data.dsol[0].Print("dsol = ",sout,EMathematicaInput);
520  LOGPZ_DEBUG(logdata,sout.str())
521  }
522 #endif
523 
524 #ifdef LOG4CXX
525  TPZFNMatrix<4> DSolLAxes(2,2), DSolRAxes(2,2);
526  DSolLAxes(0,0) = dataleft.dsol[0](0,0)*axis0DOTrL+dataleft.dsol[0](1,0)*axis1DOTrL;
527  DSolLAxes(1,0) = dataleft.dsol[0](0,0)*axis0DOTzL+dataleft.dsol[0](1,0)*axis1DOTzL;
528  DSolLAxes(0,1) = dataleft.dsol[0](0,1)*axis0DOTrL+dataleft.dsol[0](1,1)*axis1DOTrL;
529  DSolLAxes(1,1) = dataleft.dsol[0](0,1)*axis0DOTzL+dataleft.dsol[0](1,1)*axis1DOTzL;
530 
531  DSolRAxes(0,0) = dataright.dsol[0](0,0)*axis0DOTrR+dataright.dsol[0](1,0)*axis1DOTrR;
532  DSolRAxes(1,0) = dataright.dsol[0](0,0)*axis0DOTzR+dataright.dsol[0](1,0)*axis1DOTzR;
533  DSolRAxes(0,1) = dataright.dsol[0](0,1)*axis0DOTrR+dataright.dsol[0](1,1)*axis1DOTrR;
534  DSolRAxes(1,1) = dataright.dsol[0](0,1)*axis0DOTzR+dataright.dsol[0](1,1)*axis1DOTzR;
535 
536  TPZFNMatrix<9> DeformL(3,3,0.),DeformR(3,3,0.);
537  DeformL(0,0) = DSolLAxes(0,0);
538  DeformL(1,1) = DSolLAxes(1,1);
539  DeformL(0,1) = (DSolLAxes(0,1)+DSolLAxes(1,0))/2.;
540  DeformL(1,0) = DeformL(0,1);
541  DeformL(2,2) = dataleft.sol[0][0]/R;
542 
543  DeformR(0,0) = DSolRAxes(0,0);
544  DeformR(1,1) = DSolRAxes(1,1);
545  DeformR(0,1) = (DSolRAxes(0,1)+DSolRAxes(1,0))/2.;
546  DeformR(1,0) = DeformR(0,1);
547  DeformR(2,2) = dataright.sol[0][0]/R;
548 
549  TPZFNMatrix<9> TensorL(3,3,0.),TensorR(3,3,0.);
550  REAL TrDeformL, TrDeformR;
551  TrDeformL = DeformL(0,0)+DeformL(1,1)+DeformL(2,2);
552  TrDeformR = DeformR(0,0)+DeformR(1,1)+DeformR(2,2);
553  TensorL = REAL(2.*mu)*DeformL;
554  TensorR = REAL(2.*mu)*DeformR;
555  REAL normalCompL, normalCompR;
556  for (int i=0; i<3; i++) {
557  TensorL(i,i) += lambda*TrDeformL;
558  TensorR(i,i) += lambda*TrDeformR;
559  }
560  normalCompL = TensorL(1,0)*normal[0]+TensorL(1,1)*normal[1];
561  normalCompR = TensorR(1,0)*normal[0]+TensorR(1,1)*normal[1];
562  if (logdata->isDebugEnabled() && TrDeformL != 0.)
563  {
564  std::stringstream sout;
565  sout.precision(15);
566  TensorL.Print("TensorL = ",sout,EMathematicaInput);
567  TensorR.Print("TensorR = ",sout,EMathematicaInput);
568  sout << "NormalCompL = " << normalCompL << ";" << std::endl;
569  sout << "NormalCompR = " << normalCompR << ";" << std::endl;
570  LOGPZ_DEBUG(logdata,sout.str())
571  }
572 #endif
573 
574  //Calcule: Integrate { -[v].<sigma(u).n> + symmetry *[u].<sigma(v).n> + penalty*[u].[v] } ds
575 
576  // 1) Matrix Band: phi_I_Left, phi_J_Left
577 
578  for(il = 0; il< nrowl; il++ )
579  {
580  //dphiL_i/dr = dphi_i/axis0 <axes0,f_AxisR> + dphi_i/axis1 <axes1,f_AxisR>
581  dphiRZiL.PutVal(0,0, dphiL.GetVal(0,il)*axis0DOTrL + dphiL.GetVal(1,il)*axis1DOTrL );
582 
583  //dphiL_i/dz = dphi_i/axis0 <axes0,f_AxisZ> + dphi_i/axis1 <axes1,f_AxisZ>
584  dphiRZiL.PutVal(1,0, dphiL.GetVal(0,il)*axis0DOTzL + dphiL.GetVal(1,il)*axis1DOTzL );
585 
586 
587  for( jl = 0; jl < nrowl; jl++ )
588  {
589  //dphiL_j/dr = dphi_j/axis0 <axes0,f_AxisR> + dphi_j/axis1 <axes1,f_AxisR>
590  dphiRZjL.PutVal(0,0, dphiL.GetVal(0,jl)*axis0DOTrL + dphiL.GetVal(1,jl)*axis1DOTrL );
591 
592  //dphiL_j/dz = dphi_j/axis0 <axes0,f_AxisZ> + dphi_j/axis1 <axes1,f_AxisZ>
593  dphiRZjL.PutVal(1,0, dphiL.GetVal(0,jl)*axis0DOTzL + dphiL.GetVal(1,jl)*axis1DOTzL );
594 
595  double nr = normal[0];
596  double nz = normal[1];
597 
598  double term00 = -1.*(lambda*nr*phiL(il,0)*phiL(jl,0)/(2.*R) +
599  mu*nz*phiL(il,0)*dphiRZjL(1,0)/2. +
600  lambda*nr*phiL(il,0)*dphiRZjL(0,0)/2. +
601  mu*nr*phiL(il,0)*dphiRZjL(0,0));
602  term00 += symmetry*(lambda*nr*phiL(il,0)*phiL(jl,0)/(2.*R) +
603  mu*nz*phiL(jl,0)*dphiRZiL(1,0)/2. +
604  lambda*nr*phiL(jl,0)*dphiRZiL(0,0)/2. +
605  mu*nr*phiL(jl,0)*dphiRZiL(0,0));
606  term00 += beta*penalty*phiL(il,0)*phiL(jl,0);
607 
608 
609  double term01 = -1.*(lambda*nr*phiL(il,0)*dphiRZjL(1,0)/2. +
610  mu*nz*phiL(il,0)*dphiRZjL(0,0)/2.);
611  term01 += symmetry*(lambda*nz*phiL(il,0)*phiL(jl,0)/(2.*R) +
612  mu*nr*phiL(jl,0)*dphiRZiL(1,0)/2. +
613  lambda*nz*phiL(jl,0)*dphiRZiL(0,0)/2.);
614 
615 
616  double term10 = -1.*(lambda*nz*phiL(il,0)*phiL(jl,0)/(2.*R) +
617  mu*nr*phiL(il,0)*dphiRZjL(1,0)/2. +
618  lambda*nz*phiL(il,0)*dphiRZjL(0,0)/2.);
619  term10 += symmetry*(lambda*nr*phiL(jl,0)*dphiRZiL(1,0)/2. +
620  mu*nz*phiL(jl,0)*dphiRZiL(0,0)/2.);
621 
622 
623  double term11 = -1.*(lambda*nz*phiL(il,0)*dphiRZjL(1,0)/2. +
624  mu*nz*phiL(il,0)*dphiRZjL(1,0) +
625  mu*nr*phiL(il,0)*dphiRZjL(0,0)/2.);
626  term11 += symmetry*(lambda*nz*phiL(jl,0)*dphiRZiL(1,0)/2. +
627  mu*nz*phiL(jl,0)*dphiRZiL(1,0) +
628  mu*nr*phiL(jl,0)*dphiRZiL(0,0)/2.);
629  term11 +=beta*penalty*phiL(il,0)*phiL(jl,0);
630 
631 
632  ek(2*il, 2*jl) += weight*R2PI*term00;
633  ek(2*il, 2*jl+1) += weight*R2PI*term01;
634  ek(2*il+1, 2*jl) += weight*R2PI*term10;
635  ek(2*il+1, 2*jl+1) += weight*R2PI*term11;
636 
637  }
638  }
639 
640 
641  // 2) Matrix Band: phi_I_Left, phi_J_Right
642  for( il = 0; il < nrowl; il++ )
643  {
644  //dphiL_i/dr = dphi_i/axis0 <axes0,f_AxisR> + dphi_i/axis1 <axes1,f_AxisR>
645  dphiRZiL.PutVal(0,0, dphiL.GetVal(0,il)*axis0DOTrL + dphiL.GetVal(1,il)*axis1DOTrL );
646 
647  //dphiL_i/dz = dphi_i/axis0 <axes0,f_AxisZ> + dphi_i/axis1 <axes1,f_AxisZ>
648  dphiRZiL.PutVal(1,0, dphiL.GetVal(0,il)*axis0DOTzL + dphiL.GetVal(1,il)*axis1DOTzL );
649 
650  for(jr = 0; jr < nrowr; jr++ )
651  {
652  //dphiR_j/dr = dphi_j/axis0 <axes0,f_AxisR> + dphi_j/axis1 <axes1,f_AxisR>
653  dphiRZjR.PutVal(0,0, dphiR.GetVal(0,jr)*axis0DOTrR + dphiR.GetVal(1,jr)*axis1DOTrR );
654 
655  //dphiR_j/dz = dphi_j/axis0 <axes0,f_AxisZ> + dphi_j/axis1 <axes1,f_AxisZ>
656  dphiRZjR.PutVal(1,0, dphiR.GetVal(0,jr)*axis0DOTzR + dphiR.GetVal(1,jr)*axis1DOTzR );
657 
658 
659  double lambda = -((fE*fnu)/((1. + fnu)*(2.*fnu-1.)));
660  double mu = fE/(2.*(1. + fnu));
661  double nr = normal[0];
662  double nz = normal[1];
663 
664  double term02 = -1.*(lambda*nr*phiL(il,0)*phiR(jr,0)/(2.*R) +
665  mu*nz*phiL(il,0)*dphiRZjR(1,0)/2. +
666  lambda*nr*phiL(il,0)*dphiRZjR(0,0)/2. +
667  mu*nr*phiL(il,0)*dphiRZjR(0,0));
668  term02 += -1.*symmetry*(lambda*nr*phiL(il,0)*phiR(jr,0)/(2.*R) +
669  mu*nz*phiR(jr,0)*dphiRZiL(1,0)/2. +
670  lambda*nr*phiR(jr,0)*dphiRZiL(0,0)/2. +
671  mu*nr*phiR(jr,0)*dphiRZiL(0,0));
672  term02 += -1.*beta*penalty*phiL(il,0)*phiR(jr,0);
673 
674 
675  double term03 = -1.*(lambda*nr*phiL(il,0)*dphiRZjR(1,0)/2. +
676  mu*nz*phiL(il,0)*dphiRZjR(0,0)/2.);
677  term03 += -1.*symmetry*(lambda*nz*phiL(il,0)*phiR(jr,0)/(2.*R) +
678  mu*nr*phiR(jr,0)*dphiRZiL(1,0)/2. +
679  lambda*nz*phiR(jr,0)*dphiRZiL(0,0)/2.);
680 
681 
682  double term12 = -1.*(lambda*nz*phiL(il,0)*phiR(jr,0)/(2.*R) +
683  mu*nr*phiL(il,0)*dphiRZjR(1,0)/2. +
684  lambda*nz*phiL(il,0)*dphiRZjR(0,0)/2.);
685  term12 += -1.*symmetry*(lambda*nr*phiR(jr,0)*dphiRZiL(1,0)/2. +
686  mu*nz*phiR(jr,0)*dphiRZiL(0,0)/2.);
687 
688 
689  double term13 = -1.*(lambda*nz*phiL(il,0)*dphiRZjR(1,0)/2. +
690  mu*nz*phiL(il,0)*dphiRZjR(1,0) +
691  mu*nr*phiL(il,0)*dphiRZjR(0,0)/2.);
692  term13 += -1.*symmetry*(lambda*nz*phiR(jr,0)*dphiRZiL(1,0)/2. +
693  mu*nz*phiR(jr,0)*dphiRZiL(1,0) +
694  mu*nr*phiR(jr,0)*dphiRZiL(0,0)/2.);
695  term13 += -1.*beta*penalty*phiL(il,0)*phiR(jr,0);
696 
697  ek(2*il, 2*jr+2*nrowl) += weight*R2PI*term02;
698  ek(2*il, 2*jr+2*nrowl+1) += weight*R2PI*term03;
699  ek(2*il+1, 2*jr+2*nrowl) += weight*R2PI*term12;
700  ek(2*il+1, 2*jr+2*nrowl+1) += weight*R2PI*term13;
701 
702  }
703  }
704 
705  // 3) Matrix Band: phi_I_Right, phi_J_Left
706  for( ir = 0; ir < nrowr; ir++ )
707  {
708  //dphiR_i/dr = dphi_i/axis0 <axes0,f_AxisR> + dphi_i/axis1 <axes1,f_AxisR>
709  dphiRZiR.PutVal(0,0, dphiR.GetVal(0,ir)*axis0DOTrR + dphiR.GetVal(1,ir)*axis1DOTrR );
710 
711  //dphiR_i/dz = dphi_i/axis0 <axes0,f_AxisZ> + dphi_i/axis1 <axes1,f_AxisZ>
712  dphiRZiR.PutVal(1,0, dphiR.GetVal(0,ir)*axis0DOTzR + dphiR.GetVal(1,ir)*axis1DOTzR );
713 
714  for( jl = 0; jl < nrowl; jl++ )
715  {
716  //dphiL_j/dr = dphi_j/axis0 <axes0,f_AxisR> + dphi_j/axis1 <axes1,f_AxisR>
717  dphiRZjL.PutVal(0,0, dphiL.GetVal(0,jl)*axis0DOTrL + dphiL.GetVal(1,jl)*axis1DOTrL );
718 
719  //dphiL_j/dz = dphi_j/axis0 <axes0,f_AxisZ> + dphi_j/axis1 <axes1,f_AxisZ>
720  dphiRZjL.PutVal(1,0, dphiL.GetVal(0,jl)*axis0DOTzL + dphiL.GetVal(1,jl)*axis1DOTzL );
721 
722 
723  double lambda = -((fE*fnu)/((1. + fnu)*(2.*fnu-1.)));
724  double mu = fE/(2.*(1. + fnu));
725  double nr = normal[0];
726  double nz = normal[1];
727 
728  double term20 = (lambda*nr*phiR(ir,0)*phiL(jl,0)/(2.*R) +
729  mu*nz*phiR(ir,0)*dphiRZjL(1,0)/2. +
730  lambda*nr*phiR(ir,0)*dphiRZjL(0,0)/2. +
731  mu*nr*phiR(ir,0)*dphiRZjL(0,0));
732  term20 += symmetry*(lambda*nr*phiR(ir,0)*phiL(jl,0)/(2.*R) +
733  mu*nz*phiL(jl,0)*dphiRZiR(1,0)/2. +
734  lambda*nr*phiL(jl,0)*dphiRZiR(0,0)/2. +
735  mu*nr*phiL(jl,0)*dphiRZiR(0,0));
736  term20 += -1.*beta*penalty*phiR(ir,0)*phiL(jl,0);
737 
738 
739  double term21 = (lambda*nr*phiR(ir,0)*dphiRZjL(1,0)/2. +
740  mu*nz*phiR(ir,0)*dphiRZjL(0,0)/2.);
741  term21 += symmetry*(lambda*nz*phiR(ir,0)*phiL(jl,0)/(2.*R) +
742  mu*nr*phiL(jl,0)*dphiRZiR(1,0)/2. +
743  lambda*nz*phiL(jl,0)*dphiRZiR(0,0)/2.);
744 
745 
746  double term30 = (lambda*nz*phiR(ir,0)*phiL(jl,0)/(2.*R) +
747  mu*nr*phiR(ir,0)*dphiRZjL(1,0)/2. +
748  lambda*nz*phiR(ir,0)*dphiRZjL(0,0)/2.);
749  term30 += symmetry*(lambda*nr*phiL(jl,0)*dphiRZiR(1,0)/2. +
750  mu*nz*phiL(jl,0)*dphiRZiR(0,0)/2.);
751 
752 
753  double term31 = (lambda*nz*phiR(ir,0)*dphiRZjL(1,0)/2. +
754  mu*nz*phiR(ir,0)*dphiRZjL(1,0) +
755  mu*nr*phiR(ir,0)*dphiRZjL(0,0)/2.);
756  term31 += symmetry*(lambda*nz*phiL(jl,0)*dphiRZiR(1,0)/2. +
757  mu*nz*phiL(jl,0)*dphiRZiR(1,0) +
758  mu*nr*phiL(jl,0)*dphiRZiR(0,0)/2.);
759  term31 += -1.*beta*penalty*phiR(ir,0)*phiL(jl,0);
760 
761  ek(2*ir+2*nrowl, 2*jl) += weight*R2PI*term20;
762  ek(2*ir+2*nrowl, 2*jl+1) += weight*R2PI*term21;
763  ek(2*ir+2*nrowl+1, 2*jl) += weight*R2PI*term30;
764  ek(2*ir+2*nrowl+1, 2*jl+1) += weight*R2PI*term31;
765 
766  }
767  }
768 
769  // 4) Matrix Band: phi_I_Right, phi_J_Right
770 
771  for(ir = 0; ir < nrowr; ir++ )
772  {
773  //dphiR_i/dr = dphi_i/axis0 <axes0,f_AxisR> + dphi_i/axis1 <axes1,f_AxisR>
774  dphiRZiR.PutVal(0,0, dphiR.GetVal(0,ir)*axis0DOTrR + dphiR.GetVal(1,ir)*axis1DOTrR );
775 
776  //dphiR_i/dz = dphi_i/axis0 <axes0,f_AxisZ> + dphi_i/axis1 <axes1,f_AxisZ>
777  dphiRZiR.PutVal(1,0, dphiR.GetVal(0,ir)*axis0DOTzR + dphiR.GetVal(1,ir)*axis1DOTzR );
778 
779  for( jr = 0; jr < nrowr; jr++ )
780  {
781  //dphiR_j/dr = dphi_j/axis0 <axes0,f_AxisR> + dphi_j/axis1 <axes1,f_AxisR>
782  dphiRZjR.PutVal(0,0, dphiR.GetVal(0,jr)*axis0DOTrR + dphiR.GetVal(1,jr)*axis1DOTrR );
783 
784  //dphiR_j/dz = dphi_j/axis0 <axes0,f_AxisZ> + dphi_j/axis1 <axes1,f_AxisZ>
785  dphiRZjR.PutVal(1,0, dphiR.GetVal(0,jr)*axis0DOTzR + dphiR.GetVal(1,jr)*axis1DOTzR );
786 
787  double lambda = -((fE*fnu)/((1. + fnu)*(2.*fnu-1.)));
788  double mu = fE/(2.*(1. + fnu));
789  double nr = normal[0];
790  double nz = normal[1];
791 
792  double term22 = (lambda*nr*phiR(ir,0)*phiR(jr,0)/(2.*R) +
793  mu*nz*phiR(ir,0)*dphiRZjR(1,0)/2. +
794  lambda*nr*phiR(ir,0)*dphiRZjR(0,0)/2. +
795  mu*nr*phiR(ir,0)*dphiRZjR(0,0));
796  term22 += -1.*symmetry*(lambda*nr*phiR(ir,0)*phiR(jr,0)/(2.*R) +
797  mu*nz*phiR(jr,0)*dphiRZiR(1,0)/2. +
798  lambda*nr*phiR(jr,0)*dphiRZiR(0,0)/2. +
799  mu*nr*phiR(jr,0)*dphiRZiR(0,0));
800  term22 += beta*penalty*phiR(ir,0)*phiR(jr,0);
801 
802 
803  double term23 = (lambda*nr*phiR(ir,0)*dphiRZjR(1,0)/2. +
804  mu*nz*phiR(ir,0)*dphiRZjR(0,0)/2.);
805  term23 += -1.*symmetry*(lambda*nz*phiR(ir,0)*phiR(jr,0)/(2.*R) +
806  mu*nr*phiR(jr,0)*dphiRZiR(1,0)/2. +
807  lambda*nz*phiR(jr,0)*dphiRZiR(0,0)/2.);
808 
809 
810  double term32 = (lambda*nz*phiR(ir,0)*phiR(jr,0)/(2.*R) +
811  mu*nr*phiR(ir,0)*dphiRZjR(1,0)/2. +
812  lambda*nz*phiR(ir,0)*dphiRZjR(0,0)/2.);
813  term32 += -1.*symmetry*(lambda*nr*phiR(jr,0)*dphiRZiR(1,0)/2. +
814  mu*nz*phiR(jr,0)*dphiRZiR(0,0)/2.);
815 
816 
817  double term33 = (lambda*nz*phiR(ir,0)*dphiRZjR(1,0)/2. +
818  mu*nz*phiR(ir,0)*dphiRZjR(1,0) +
819  mu*nr*phiR(ir,0)*dphiRZjR(0,0)/2.);
820  term33 += -1.*symmetry*(lambda*nz*phiR(jr,0)*dphiRZiR(1,0)/2. +
821  mu*nz*phiR(jr,0)*dphiRZiR(1,0) +
822  mu*nr*phiR(jr,0)*dphiRZiR(0,0)/2.);
823  term33 += beta*penalty*phiR(ir,0)*phiR(jr,0);
824 
825 
826  ek(2*ir+2*nrowl, 2*jr+2*nrowl) += weight*R2PI*term22;
827  ek(2*ir+2*nrowl, 2*jr+2*nrowl+1) += weight*R2PI*term23;
828  ek(2*ir+2*nrowl+1, 2*jr+2*nrowl) += weight*R2PI*term32;
829  ek(2*ir+2*nrowl+1, 2*jr+2*nrowl+1) += weight*R2PI*term33;
830 
831  }
832  }
833 #ifdef LOG4CXX
834  if(logdata->isDebugEnabled())
835  {
836  std::stringstream sout;
837  ek.Print("ek = ",sout,EMathematicaInput);
838  LOGPZ_DEBUG(logdata,sout.str())
839  }
840 #endif
841 }
842 
843 //-------------------------------------------------------------------
844 
846 int TPZElasticityAxiMaterial::VariableIndex(const std::string &name)
847 {
848  if(!strcmp("Eigenvector1",name.c_str())) return 9;
849  if(!strcmp("Eigenvector2",name.c_str())) return 1;
850  if(!strcmp("Eigenvector3",name.c_str())) return 2;
851  if(!strcmp("Sigmarr",name.c_str())) return 3;
852  if(!strcmp("Sigmazz",name.c_str())) return 4;
853  if(!strcmp("Sigmatt",name.c_str())) return 5;
854  if(!strcmp("Taurz",name.c_str())) return 6;
855  if(!strcmp("displacement",name.c_str())) return 7;
856  if(!strcmp("MohrCoulomb",name.c_str())) return 8;
857 
858  return TPZMaterial::VariableIndex(name);
859  return -1;
860 }
861 
864 int TPZElasticityAxiMaterial::NSolutionVariables(int var)
865 {
866  switch(var) {
867  case 9:
868  return 3;
869  case 1:
870  return 3;
871  case 2:
872  return 3;
873  case 3:
874  return 1;
875  case 4:
876  return 1;
877  case 5:
878  return 1;
879  case 6:
880  return 1;
881  case 7:
882  return 3;
883  case 8:
884  return 1;
885  default:
886  return TPZMaterial::NSolutionVariables(var);
887  return 0;
888  }
889 }
890 
893 void TPZElasticityAxiMaterial::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout)
894 {
895  if(var == 0)
896  {
897  TPZMaterial::Solution(data,var,Solout);
898  return;
899  }
900  int64_t numbersol = data.dsol.size();
901  if (numbersol != 1) {
902  DebugStop();
903  }
904 
905  TPZFMatrix<REAL> &axes = data.axes;
906  TPZVec<STATE> &SolAxes = data.sol[0];
907  TPZFMatrix<STATE> &DSolAxes = data.dsol[0];
908 
909  // R = Dot[{data.x - origin},{AxisR}] ***because AxisR is already normalized!
910  REAL R = (data.x[0] - f_Origin[0])*f_AxisR[0] + (data.x[1] - f_Origin[1])*f_AxisR[1] + (data.x[2] - f_Origin[2])*f_AxisR[2];
911  REAL Z = (data.x[0] - f_Origin[0])*f_AxisZ[0] + (data.x[1] - f_Origin[1])*f_AxisZ[1] + (data.x[2] - f_Origin[2])*f_AxisZ[2];
912 
913  int s = (R > 0)? 1:-1;
914  R = fabs(R);
915  if(R < 1.e-10) R = 1.e-10;
916 
917  REAL DelTemp(fDelTemperature);
918 
919  if (fTemperatureFunction) {
920  TPZManVector<REAL,2> RZ(2);
921  RZ[0] = R;
922  RZ[1] = Z;
923  fTemperatureFunction(RZ,DelTemp);
924  }
925 
926 
927  double axis0DOTr = 0., axis1DOTr = 0., axis0DOTz = 0., axis1DOTz = 0.;
928  for(int pos = 0; pos < 3; pos++)
929  {
930  axis0DOTr += axes.GetVal(0,pos) * f_AxisR[pos] * s;
931  axis1DOTr += axes.GetVal(1,pos) * f_AxisR[pos] * s;
932  axis0DOTz += axes.GetVal(0,pos) * f_AxisZ[pos];
933  axis1DOTz += axes.GetVal(1,pos) * f_AxisZ[pos];
934  }
935  TPZFNMatrix<9> DSolrz(2,2,0.);
936  DSolrz.PutVal(0,0, DSolAxes(0,0)*axis0DOTr + DSolAxes(1,0)*axis1DOTr );
937  DSolrz.PutVal(0,1, DSolAxes(0,0)*axis0DOTz + DSolAxes(1,0)*axis1DOTz );
938  DSolrz.PutVal(1,0, DSolAxes(0,1)*axis0DOTr + DSolAxes(1,1)*axis1DOTr );
939  DSolrz.PutVal(1,1, DSolAxes(0,1)*axis0DOTz + DSolAxes(1,1)*axis1DOTz );
940 
941 
942 #ifdef LOG4CXX
943  if (logger->isDebugEnabled())
944  {
945  std::stringstream sout;
946  sout << "Point " << data.x << std::endl;
947  sout << "Solution " << data.sol << std::endl;
948  DSolrz.Print("Derivatives of the solution\n",sout);
949  sout << "Radius " << R << std::endl;
950  LOGPZ_DEBUG(logger,sout.str());
951  }
952 #endif
953 
954  //Infinitesimal Tensor
955  TPZFNMatrix<9> Einf(3,3,0.);
956  Einf.PutVal(0,0,DSolrz(0,0));
957  Einf.PutVal(0,1,0.5*(DSolrz(0,1) + DSolrz(1,0)));
958  Einf.PutVal(1,0,0.5*(DSolrz(0,1) + DSolrz(1,0)));
959  Einf.PutVal(1,1,DSolrz(1,1));
960  Einf.PutVal(2,2,data.sol[0][0]/R);
961 
962 #ifdef LOG4CXX
963  if (logger->isDebugEnabled())
964  {
965  std::stringstream sout;
966  Einf.Print("Deformation tensor",sout);
967  LOGPZ_DEBUG(logger,sout.str());
968  }
969 #endif
970 
971  double lambda = -((fE*fnu)/((1. + fnu)*(2.*fnu-1.)));
972  double mi = fE/(2.*(1. + fnu));
973  double trE = Einf(0,0) + Einf(1,1) + Einf(2,2);
974 
975 #ifdef LOG4CXX
976  if (logger->isDebugEnabled())
977  {
978  std::stringstream sout;
979  sout << "E = " << fE << " fnu = " << fnu << " mi = " << mi << " lambda = " << lambda << " trE = " << trE;
980  LOGPZ_DEBUG(logger,sout.str());
981  }
982 #endif
983 
984  //Stress Tensor
985  TPZFNMatrix<9> T(3,3,0.);
986  double cte = lambda*trE;
987  REAL epsT = DelTemp*fAlpha;
988  REAL TensorThermico = (3.*lambda+2.*mi)*epsT;
989 
990 
991  //T = lambda.tr(E).I + 2.mi.E
992  T.PutVal(0,0,cte + 2.*mi*Einf(0,0)-TensorThermico);
993  T.PutVal(0,1, 2.*mi*Einf(0,1));
994  T.PutVal(0,2, 2.*mi*Einf(0,2));
995  T.PutVal(1,0, 2.*mi*Einf(1,0));
996  T.PutVal(1,1,cte + 2.*mi*Einf(1,1)-TensorThermico);
997  T.PutVal(1,2, 2.*mi*Einf(1,2));
998  T.PutVal(2,0, 2.*mi*Einf(2,0));
999  T.PutVal(2,1, 2.*mi*Einf(2,1));
1000  T.PutVal(2,2,cte + 2.*mi*Einf(2,2)-TensorThermico);
1001 
1002 #ifdef LOG4CXX
1003  if (logger->isDebugEnabled())
1004  {
1005  std::stringstream sout;
1006  T.Print("Stress tensor",sout);
1007  LOGPZ_DEBUG(logger,sout.str());
1008  }
1009 #endif
1010 
1011  switch(var)
1012  {
1013  case 9: //Solout = 1stEigenvalue * {1stEigenvector}
1014  {
1015  int64_t NumIt = 1000;
1016  REAL tol = 1.E-5;
1017  TPZVec<REAL> EigValues(3,0.);
1018  TPZFNMatrix<9,REAL> EigVectors(3,3,0.);
1019  bool EigenWorks;
1020  EigenWorks = T.SolveEigensystemJacobi(NumIt, tol, EigValues, EigVectors);
1021  if(EigenWorks)
1022  {
1023  Solout.Resize(3);
1024  for(int i = 0; i < 3; i++) Solout[i] = EigValues[0] * EigVectors(0,i);
1025  }
1026  else cout << "TPZElasticityAxiMaterial::Solution Error -> case 0\n";
1027 #ifdef LOG4CXX
1028  if (logger->isDebugEnabled())
1029  {
1030  std::stringstream sout;
1031  sout << "First eigenvector " << std::endl;
1032  sout << Solout;
1033  LOGPZ_DEBUG(logger,sout.str());
1034  }
1035 #endif
1036  }
1037  break;
1038 
1039  case 1: //Solout = 2ndEigenvalue * {2ndEigenvector}
1040  {
1041  int64_t NumIt = 1000;
1042  REAL tol = 1.E-5;
1043  TPZVec<REAL> EigValues(3,0.);
1044  TPZFNMatrix<9,REAL> EigVectors(3,3,0.);
1045  bool EigenWorks;
1046  EigenWorks = T.SolveEigensystemJacobi(NumIt, tol, EigValues, EigVectors);
1047  if(EigenWorks)
1048  {
1049  Solout.Resize(3);
1050  for(int i = 0; i < 3; i++) Solout[i] = EigValues[1] * EigVectors(1,i);
1051  }
1052  else cout << "TPZElasticityAxiMaterial::Solution Error -> case 1\n";
1053 #ifdef LOG4CXX
1054  if (logger->isDebugEnabled())
1055  {
1056  std::stringstream sout;
1057  sout << "Second eigenvector " << std::endl;
1058  sout << Solout;
1059  LOGPZ_DEBUG(logger,sout.str());
1060  }
1061 #endif
1062  }
1063  break;
1064 
1065  case 2: //Solout = 3rdEigenvalue * {3rdEigenvector}
1066  {
1067  int64_t NumIt = 1000;
1068  REAL tol = 1.E-5;
1069  TPZVec<REAL> EigValues(3,0.);
1070  TPZFNMatrix<9,REAL> EigVectors(3,3,0.);
1071  bool EigenWorks;
1072  EigenWorks = T.SolveEigensystemJacobi(NumIt, tol, EigValues, EigVectors);
1073  if(EigenWorks)
1074  {
1075  Solout.Resize(3);
1076  for(int i = 0; i < 3; i++) Solout[i] = EigValues[2] * EigVectors(2,i);
1077  }
1078  else cout << "TPZElasticityAxiMaterial::Solution Error -> case 2\n";
1079 #ifdef LOG4CXX
1080  if (logger->isDebugEnabled())
1081  {
1082  std::stringstream sout;
1083  sout << "Third eigenvector " << std::endl;
1084  sout << Solout;
1085  LOGPZ_DEBUG(logger,sout.str());
1086  }
1087 #endif
1088  }
1089  break;
1090 
1091  case 3: //Solout = Sigma_r,r
1092  {
1093  Solout.Resize(1);
1094  Solout[0] = T(0,0);
1095 #ifdef LOG4CXX
1096  if (logger->isDebugEnabled())
1097  {
1098  std::stringstream sout;
1099  sout << "Sigma r " << T(0,0);
1100  LOGPZ_DEBUG(logger,sout.str());
1101  }
1102 #endif
1103  }
1104  break;
1105 
1106  case 4: //Solout = Sigma_z,z
1107  {
1108  Solout.Resize(1);
1109  Solout[0] = T(1,1);
1110 #ifdef LOG4CXX
1111  if (logger->isDebugEnabled())
1112  {
1113  std::stringstream sout;
1114  sout << "Sigma z " << T(1,1);
1115  LOGPZ_DEBUG(logger,sout.str());
1116  }
1117 #endif
1118  }
1119  break;
1120 
1121  case 5: //Solout = Sigma_theta,theta
1122  {
1123  Solout.Resize(1);
1124  Solout[0] = T(2,2);
1125 #ifdef LOG4CXX
1126  if (logger->isDebugEnabled())
1127  {
1128  std::stringstream sout;
1129  sout << "Sigma theta " << T(2,2);
1130  LOGPZ_DEBUG(logger,sout.str());
1131  }
1132 #endif
1133  }
1134  break;
1135 
1136  case 6: //Solout = Tau_r,z
1137  {
1138  Solout.Resize(1);
1139  Solout[0] = T(0,1);
1140 #ifdef LOG4CXX
1141  if (logger->isDebugEnabled())
1142  {
1143  std::stringstream sout;
1144  sout << "Sigma rz " << T(0,1);
1145  LOGPZ_DEBUG(logger,sout.str());
1146  }
1147 #endif
1148  }
1149  break;
1150 
1151  case 7: //Solout = Displacement
1152  {
1153  Solout.Resize(3);
1154  Solout[0] = SolAxes[0];
1155  Solout[1] = SolAxes[1];
1156  Solout[2] = SolAxes[2];
1157  }
1158  break;
1159 
1160  case 8: //MohrCoulomb plasticity criteria
1161  {
1162  Solout.Resize(1);
1163  REAL i1, i2, i3, j1, j2, j3;
1164 
1165  i1 = T(0,0) + T(1,1) + T(2,2);
1166  TPZFMatrix<REAL> T2(3,3,0.);
1167  T.Multiply(T,T2);
1168 
1169  i2 = 0.5*( i1*i1 - (T2(0,0) + T2(1,1) + T2(2,2)) );
1170 
1171  i3 = T(0,0)*T(1,1)*T(2,2) + T(0,1)*T(1,2)*T(2,0) + T(0,2)*T(1,0)*T(2,1) - T(0,2)*T(1,1)*T(2,0) - T(1,2)*T(2,1)*T(0,0) - T(2,2)*T(0,1)*T(1,0);
1172 
1173  j1 = 0.;
1174  j2 = 1./3.*(i1*i1 - 3.*i2);
1175  j3 = 1./27.*(2.*i1*i1*i1 - 9.*i1*i2 + 27.*i3);
1176 
1177  REAL cos3theta = 3.*sqrt(3.)/2. * j3/(pow(j2,(REAL)1.5));
1178 
1179  REAL theta = acos(cos3theta)/3.;
1180 
1181  Solout[0] = 1./3.*i1*sin(f_phi) + sqrt(i2)*sin(theta + M_PI/3.) + sqrt(i2/3.)*cos(theta + M_PI/3.)*sin(f_phi) - f_c*cos(f_phi);
1182 #ifdef LOG4CXX
1183  if (logger->isDebugEnabled())
1184  {
1185  std::stringstream sout;
1186  sout << "Criterio de Mohr Coulomb \ni1 " << i1 << " i2 " << i2 <<
1187  " i3 " << i3 << " \nj1 " << j1 << " j2 " << j2 << " j3 " << j3
1188  << " \nf_phi " << f_phi <<
1189  " f_c " << f_c << " theta " << theta << " MohrCoul " << Solout[0];
1190  LOGPZ_DEBUG(logger,sout.str());
1191  }
1192 #endif
1193  }
1194  break;
1195 
1196  default:
1197  {
1198  cout << "TPZElasticityAxiMaterial::Solution Error -> default\n";
1199  TPZMaterial::Solution(data.sol[0],data.dsol[0],data.axes,var,Solout);
1200  }
1201  break;
1202  }
1203 }
1204 
1205 void TPZElasticityAxiMaterial::Flux(TPZVec<REAL> &x, TPZVec<STATE> &Sol, TPZFMatrix<STATE> &DSol, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux)
1206 {
1207  if(fabs(axes(2,0)) >= 1.e-6 || fabs(axes(2,1)) >= 1.e-6)
1208  {
1209  cout << "TPZElasticityAxiMaterial::Flux only serves for xy configuration\n";
1210  axes.Print("axes");
1211  }
1212 }
1213 
1214 void TPZElasticityAxiMaterial::Errors(TPZVec<REAL> &x,TPZVec<STATE> &u, TPZFMatrix<STATE> &dudaxes, TPZFMatrix<REAL> &axes, TPZVec<STATE> &flux,
1215  TPZVec<STATE> &u_exact,TPZFMatrix<STATE> &du_exact,TPZVec<REAL> &values)
1216 {
1217  values[0] = 0.;
1218  TPZManVector<REAL> sigma(3,0.),sigma_exact(3,0.);
1219  REAL sigx,sigy,sigxy,gamma;
1220 
1221  TPZFNMatrix<9> du;
1222  TPZFMatrix<REAL> dudaxesReal = (TPZFMatrix<REAL> &)dudaxes;
1223  TPZAxesTools<REAL>::Axes2XYZ(dudaxesReal, du, axes);
1224 
1225  //tensoes aproximadas : uma forma
1226  gamma = du(1,0)+du(0,1);
1227  sigma[0] = fEover1MinNu2*(du(0,0)+fnu*du(1,1));
1228  sigma[1] = fEover1MinNu2*(fnu*du(0,0)+du(1,1));
1229  sigma[2] = fE*0.5/(1.+fnu)*gamma;
1230 
1231  TPZMaterialData mydata;
1232  mydata.sol[0] = u;
1233  mydata.dsol[0] = (TPZFNMatrix<9,STATE> &)du;
1234  mydata.axes = axes;
1235  mydata.x = x;
1236 
1237  //tensoes aproximadas : outra forma
1238  TPZVec<STATE> sol(1);
1239  Solution(mydata,5,sol);
1240  sigma[0] = sol[0];
1241  Solution(mydata,6,sol);
1242  sigma[1] = sol[0];
1243  Solution(mydata,8,sol);
1244  sigma[2] = sol[0];
1245 
1246  //exata
1247  gamma = du_exact(1,0)+du_exact(0,1);
1248  sigma_exact[0] = fEover1MinNu2*(du_exact(0,0)+fnu*du_exact(1,1));
1249  sigma_exact[1] = fEover1MinNu2*(fnu*du_exact(0,0)+du_exact(1,1));
1250  sigma_exact[2] = fE*0.5/(1.+fnu)*gamma;
1251  sigx = (sigma[0] - sigma_exact[0]);
1252  sigy = (sigma[1] - sigma_exact[1]);
1253  sigxy = (sigma[2] - sigma_exact[2]);
1254  //values[0] = calculo do erro estimado em norma Energia
1255  values[0] = fE*(sigx*sigx + sigy*sigy + 2*fnu*sigx*sigy)/(1-fnu*fnu);
1256  values[0] = (values[0] + .5*fE*sigxy*sigxy/(1+fnu));
1257 
1258  //values[1] : erro em norma L2 em tensoes
1259  //values[1] = sigx*sigx + sigy*sigy + sigxy*sigxy;
1260 
1261  //values[1] : erro em norma L2 em deslocamentos
1262  values[1] = pow((REAL)fabs(u[0] - u_exact[0]),(REAL)2.0)+pow((REAL)fabs(u[1] - u_exact[1]),(REAL)2.0); // It is important when REAL is long double, because fabs(...) returns long double then pow() must to return long double - Jorge
1263 
1264  //values[2] : erro estimado
1265  values[2] = 0.;
1266 }
1267 
1268 
1269 int TPZElasticityAxiMaterial::ClassId() const{
1270  return Hash("TPZElasticityAxiMaterial") ^ TPZDiscontinuousGalerkin::ClassId() << 1;
1271 }
1272 
1273 #ifndef BORLAND
1274 template class TPZRestoreClass<TPZElasticityAxiMaterial>;
1275 #endif
1276 
1277 void TPZElasticityAxiMaterial::Read(TPZStream &buf, void *context)
1278 {
1279  TPZMaterial::Read(buf,context);
1280  buf.Read(&fE,1);
1281  buf.Read(&fnu,1);
1282  buf.Read(&fEover21PlusNu,1);
1283  buf.Read(&fEover1MinNu2,1);
1284 
1285  buf.Read(ff,3);
1286 }
1287 
1288 void TPZElasticityAxiMaterial::Write(TPZStream &buf, int withclassid) const
1289 {
1290  TPZMaterial::Write(buf,withclassid);
1291  buf.Write(&fE,1);
1292  buf.Write(&fnu,1);
1293  buf.Write(&fEover21PlusNu,1);
1294  buf.Write(&fEover1MinNu2,1);
1295 
1296  buf.Write(ff,3);
1297 }
TPZDiscontinuousGalerkin
Defines the interface which material objects need to implement for discontinuous Galerkin formulation...
Definition: pzdiscgal.h:20
TPZElasticityAxiMaterial::fPenaltyConstant
REAL fPenaltyConstant
Penalty term.
Definition: pzelasAXImat.h:210
TPZElasticityAxiMaterial::Errors
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: pzelasAXImat.cpp:1214
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
TPZVec::Print
void Print(std::ostream &out=std::cout)
Prints the structural information of the vector object to the output stream. This method will not p...
Definition: pzvec.h:482
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
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
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
TPZElasticityAxiMaterial::f_AxisZ
TPZManVector< REAL > f_AxisZ
Revolution Axis.
Definition: pzelasAXImat.h:201
TPZElasticityAxiMaterial
Implements a two dimensional elastic material in plane stress or strain.
Definition: pzelasAXImat.h:21
TPZElasticityAxiMaterial::fAlpha
REAL fAlpha
Thermal expansion coeficient.
Definition: pzelasAXImat.h:183
TPZElasticityAxiMaterial::fDelTemperature
REAL fDelTemperature
Temperature difference.
Definition: pzelasAXImat.h:189
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZElasticityAxiMaterial::f_Origin
TPZManVector< REAL > f_Origin
Origin of AxisR and AxisZ.
Definition: pzelasAXImat.h:204
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
TPZElasticityAxiMaterial::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: pzelasAXImat.cpp:1205
TPZElasticityAxiMaterial::ContributeInterface
virtual void ContributeInterface(TPZMaterialData &data, TPZMaterialData &dataleft, TPZMaterialData &dataright, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
It computes a contribution to stiffness matrix and load vector at one integration point...
Definition: pzelasAXImat.cpp:421
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.
TPZElasticityAxiMaterial::ClassId
virtual int ClassId() const override
Unique identifier for serialization purposes.
Definition: pzelasAXImat.cpp:1269
TPZElasticityAxiMaterial::f_phi
REAL f_phi
Mohr Coulomb Plasticity Criteria Data.
Definition: pzelasAXImat.h:172
TPZElasticityAxiMaterial::SetOrigin
void SetOrigin(TPZManVector< REAL > &Orig, TPZManVector< REAL > &AxisZ, TPZManVector< REAL > &AxisR)
Set the origin of Revolution Axis ( ), the direction of Revolution Axis ( ), and the Radius vector (...
Definition: pzelasAXImat.cpp:107
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
TPZElasticityAxiMaterial::E
REAL E()
Returns the elasticity modulus E.
Definition: pzelasAXImat.h:142
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::p
int p
maximum polinomial order of the shape functions
Definition: pzmaterialdata.h:75
TPZElasticityAxiMaterial::TPZElasticityAxiMaterial
TPZElasticityAxiMaterial()
Default constructor.
Definition: pzelasAXImat.cpp:25
TPZElasticityAxiMaterial::Write
virtual void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzelasAXImat.cpp:1288
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZMaterial::Print
virtual void Print(std::ostream &out=std::cout)
Prints out the data associated with the material.
Definition: TPZMaterial.cpp:103
pzelmat.h
Contains declaration of TPZElementMatrix struct which associates an element matrix with the coeficien...
TPZElasticityAxiMaterial::GetAxisZ
TPZManVector< REAL > GetAxisZ()
Definition: pzelasAXImat.cpp:145
TPZElasticityAxiMaterial::fSymmetric
REAL fSymmetric
Symmetric.
Definition: pzelasAXImat.h:207
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
acos
expr_ expr_ expr_ expr_ acos
Definition: tfadfunc.h:75
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
TPZElasticityAxiMaterial::ComputeR
REAL ComputeR(TPZVec< REAL > &x)
Definition: pzelasAXImat.cpp:135
sin
sin
Definition: fadfunc.h:63
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
TPZElasticityAxiMaterial::ff
REAL ff[3]
Forcing vector.
Definition: pzelasAXImat.h:186
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZElasticityAxiMaterial::Print
virtual void Print(std::ostream &out=std::cout) override
Prints the material data.
Definition: pzelasAXImat.cpp:162
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
TPZElasticityAxiMaterial::GetAxisR
TPZManVector< REAL > GetAxisR()
Definition: pzelasAXImat.cpp:140
LOGPZ_DEBUG
#define LOGPZ_DEBUG(A, B)
Define log for debug info.
Definition: pzlog.h:87
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZElasticityAxiMaterial::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns the variable index associated with the name.
Definition: pzelasAXImat.cpp:846
TPZMaterialData::PrintMathematica
void PrintMathematica(std::ostream &out) const
Prints the data in a format suitable for Mathematica.
Definition: pzmaterialdata.cpp:193
TPZElasticityAxiMaterial::~TPZElasticityAxiMaterial
virtual ~TPZElasticityAxiMaterial()
Destructor.
Definition: pzelasAXImat.cpp:155
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
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
TPZElasticityAxiMaterial::Read
virtual void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzelasAXImat.cpp:1277
TPZElasticityAxiMaterial::NSolutionVariables
virtual int NSolutionVariables(int var) override
Returns the number of variables associated with the variable indexed by var. var is obtained by cal...
Definition: pzelasAXImat.cpp:864
TPZMaterialData::HSize
REAL HSize
measure of the size of the element
Definition: pzmaterialdata.h:83
TPZElasticityAxiMaterial::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: pzelasAXImat.cpp:893
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
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
pzmatrix.h
Contains TPZMatrix<TVar>class, root matrix class.
TPZElasticityAxiMaterial::Name
std::string Name() override
Returns the material name.
Definition: pzelasAXImat.h:96
TPZElasticityAxiMaterial::f_c
REAL f_c
Mohr Coulomb Plasticity Criteria Data.
Definition: pzelasAXImat.h:174
TPZElasticityAxiMaterial::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef) override
Calculates the element stiffness matrix.
Definition: pzelasAXImat.cpp:171
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZMaterial::NSolutionVariables
virtual int NSolutionVariables(int var)
Returns the number of variables associated with the variable indexed by var.
Definition: TPZMaterial.cpp:176
TPZElasticityAxiMaterial::fTemperatureFunction
void(* fTemperatureFunction)(const TPZVec< REAL > &rz, REAL &temperature)
Function which defines the temperature.
Definition: pzelasAXImat.h:213
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZElasticityAxiMaterial::GetOrigin
TPZManVector< REAL > GetOrigin()
Definition: pzelasAXImat.cpp:150
EMathematicaInput
Definition: pzmatrix.h:55
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
TPZElasticityAxiMaterial::fnu
REAL fnu
Poison coeficient.
Definition: pzelasAXImat.h:180
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
TPZElasticityAxiMaterial::NStateVariables
virtual int NStateVariables() const override
Returns the number of state variables associated with the material.
Definition: pzelasAXImat.cpp:158
rdt.values
def values
Definition: rdt.py:119
TPZElasticityAxiMaterial::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Applies the element boundary conditions.
Definition: pzelasAXImat.cpp:306
TPZElasticityAxiMaterial::fE
REAL fE
Elasticity modulus.
Definition: pzelasAXImat.h:177
pzelasAXImat.h
Contains the TPZElasticityAxiMaterial class which implements a two dimensional elastic material in pl...
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
TPZFMatrix::PutVal
int PutVal(const int64_t row, const int64_t col, const TVar &value) override
Put values without bounds checking This method is faster than "Put" if DEBUG is defined.
Definition: pzfmatrix.h:548
TPZFMatrix::GetVal
const TVar & GetVal(const int64_t row, const int64_t col) const override
Get values without bounds checking This method is faster than "Get" if DEBUG is defined.
Definition: pzfmatrix.h:566
TPZElasticityAxiMaterial::f_AxisR
TPZManVector< REAL > f_AxisR
Direction of Surface.
Definition: pzelasAXImat.h:198
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
TPZFNMatrix
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
TPZStream::Read
virtual void Read(bool &val)
Definition: TPZStream.cpp:91
Generated on Sun Aug 16 2020 11:17:03 for NeoPZ by   doxygen 1.8.13