NeoPZ
TPZYCWillamWarnke.h
Go to the documentation of this file.
1 
5 #ifndef TPZYCWILLAMWARNKE_H
6 #define TPZYCWILLAMWARNKE_H
7 
8 #include "TPZTensor.h"
9 #include "pzfmatrix.h"
10 
11 #include "TPZSavable.h"
12 #include "TPZPlasticCriterion.h"
13 
19 class TPZYCWillamWarnke: public TPZPlasticCriterion {
20 
21 public:
22 
23  TPZYCWillamWarnke():fa(1.),fb(1.),fcConcrete(20.),fcoesion(9.2376),fa0(0.1025),fa1(-0.8403),fa2(-0.0910),fb0(0.1025),fb1(-0.4507),fb2(-0.1018){}
24 
25  TPZYCWillamWarnke(const TPZYCWillamWarnke & source)
26  {
27  fa = source.fa;
28  fb = source.fb;
29  fcConcrete = source.fcConcrete;
30  fa0 = source.fa0;
31  fa1 = source.fa1;
32  fa2 = source.fa2;
33  fb0 = source.fb0;
34  fb1 = source.fb1;
35  fb2 = source.fb2;
36  }
37 
38  enum {NYield = 1};
39 
40  const char * Name() const
41  {
42  return "TPZYCDruckerPrager";
43  }
44 
45  void Print(std::ostream & out) const override
46  {
47  out << Name();
48  }
49 
53  void SetUp(const REAL a, const REAL b,const REAL fConcrete)
54  {
55  fa = a;
56  fb = b;
57  fcConcrete = fConcrete;
58 
59  }
60 
61  int GetForceYield()
62  {
63  return 0; //nothing to be done in this yield criterium
64  }
65 
66  void SetForceYield(const int forceYield)
67  {
68  //nothing to be done in this yield criterium
69  }
70 
78  void SetYieldStatusMode(const TPZTensor<REAL> & sigma, const REAL & A)
79  {
80  // nothing to be done in this yield criterium
81  }
82 
90  template < class T>
91  void Compute(const TPZTensor<T> & sigma, const T & A, TPZVec<T> &res, int checkForcedYield = 0) const;
92 
100  template <class T>
101  void N(const TPZTensor<T> & sigma,const T & A, TPZVec<TPZTensor<T> > & Ndir, int checkForcedYield = 0) const;
102 
110  template <class T>
111  void H(const TPZTensor<T> & sigma,const T & A, TPZVec<T> & h, int checkForcedYield = 0) const;
112 
116  template <class T>
117  void AlphaMultiplier(const T &A, T &multiplier) const
118  {
119  multiplier = T(1.);
120  }
121 
122 // void Write(TPZStream &buf, int withclassid) const override;
123 // void Read(TPZStream &buf, void *context) override;
124  public:
125 virtual int ClassId() const override;
126 
127 
128  TPZTensor<REAL> gRefTension;
129 
130  int NumCases()
131  {
132  return 1;
133  }
137  void LoadState(TPZFMatrix<REAL> &state)
138  {
139  int i;
140  for(i=0; i<6; i++) gRefTension[i] = state(i,0);
141  }
142 
143  void ComputeTangent(TPZFMatrix<REAL> &tangent, TPZVec<REAL> &coefs, int icase)
144  {
145  switch(icase)
146  {
147 
148  case 0:
149  {
150  TPZVec<TPZTensor<REAL> > Ndir(1);
151  REAL yield = 1.e6;
152  this->N<REAL>(gRefTension,yield, Ndir, 0);
153  tangent.Redim(1,6);
154  for(int i=0; i<6; i++)
155  {
156  tangent(0,i) = Ndir[0][i];
157  }
158  break;
159  }
160 
161  }
162  }
163 
164  void Residual(TPZFMatrix<REAL> &res,int icase)
165  {
166 
167  res.Redim(1,1);
168  TPZTensor<REAL> grad;
169 
170  switch(icase)
171  {
172  case 0:
173  {
174  TPZVec<REAL> phi(1);
175  REAL yield = 1.e6;
176  this->Compute(gRefTension,yield,phi,0);
177  res(0,0) = phi[0];
178  break;
179  }
180  }
181 
182  }
183 
184  void Write(TPZStream &out, int withclassid = 0) const override
185  {
186  out.Write(&fPhi);
187  out.Write(&fa);
188  out.Write(&fb);
189  out.Write(&fcConcrete);
190  out.Write(&fcoesion);
191  out.Write(&fa0);
192  out.Write(&fa1);
193  out.Write(&fa2);
194  out.Write(&fb0);
195  out.Write(&fb1);
196  out.Write(&fb2);
197  }
198 
199  void Read(TPZStream &input, void *context = 0) override
200  {
201  input.Read(&fPhi);
202  input.Read(&fa);
203  input.Read(&fb);
204  input.Read(&fcConcrete);
205  input.Read(&fcoesion);
206  input.Read(&fa0);
207  input.Read(&fa1);
208  input.Read(&fa2);
209  input.Read(&fb0);
210  input.Read(&fb1);
211  input.Read(&fb2);
212  }
213 
214  void YieldFunction(const TPZVec<STATE>& sigma, STATE kprev, TPZVec<STATE>& yield) const override{
215  TPZTensor<STATE> sigmaTensor;
216  sigmaTensor.XX() = sigma[0];
217  sigmaTensor.YY() = sigma[1];
218  sigmaTensor.ZZ() = sigma[2];
219  Compute(sigmaTensor, kprev, yield, 0);
220  }
221 
222  virtual int GetNYield() const override{
223  return as_integer(NYield);
224  }
225 
226 public:
227  REAL fPhi;
228  REAL fa, fb, fcConcrete,fcoesion,fa0,fa1,fa2,fb0,fb1,fb2;
229 
230 protected:
231 
232 
234 };
235 
236 
237 template <class T>
238 void TPZYCWillamWarnke::Compute(const TPZTensor<T> & sigma,const T & A, TPZVec<T> &res, int checkForcedYield) const
239 {
240 
241  TPZTensor<T> gradlode,dJ2,dJ3;
242  T lode,I1,J2,J3;
243  TPZTensor<T> eigenval,dSigma1,dSigma2,dSigma3;
244  sigma.Eigenvalue(eigenval, dSigma1, dSigma2, dSigma3);
245 
246  I1 = sigma.I1();
247  J2 = sigma.J2();
248  J3 = sigma.J3();
249  sigma.dJ2(dJ2);
250  sigma.dJ3(dJ3);
251 /*
252  T I13=I1/T(3.);
253  T rhoc = -(1./(2.*fb2))*(fb1+ sqrt((fb1*fb1)-T(4.)*T(fb2)*(T(fb0)-I13)));
254  T rhot = -(1./(2.*fa2))*(fa1+ sqrt((fa1*fa1)-T(4.)*T(fa2)*(T(fa0)-I13)));
255 
256 // std::cout << " rhoc " << std::endl;
257 // std::cout << rhoc<<std::endl;
258 //
259 // std::cout << " rhot " << std::endl;
260 // std::cout << rhot <<std::endl;
261 
262  sigma.Lodeangle(gradlode,lode);
263 
264  T s = (T(2.)*rhoc)*((rhoc*rhoc)-(rhot*rhot));
265  T u = T(4.)*(rhoc*rhoc-rhot*rhot)*cos(lode)*cos(lode)+T(5.)*rhot*rhot-T(4.)*rhot*rhoc;
266  T t = rhoc*(T(2.)*rhot-rhoc)*sqrt(u);
267  T v = T(4.)*(rhoc*rhoc)*cos(lode)*cos(lode)+(rhoc-T(2.)*rhot)*(rhoc-T(2.)*rhot);
268 // std::cout << " s " << std::endl;
269 // std::cout << s <<std::endl;
270 // std::cout << " u " << std::endl;
271 // std::cout << u <<std::endl;
272 // std::cout << " v " << std::endl;
273 // std::cout << v <<std::endl;
274 // std::cout << " t " << std::endl;
275 // std::cout << t <<std::endl;
276  res[0] = (s+t)/v;
277  */
278 // std::cout << " PHI " << std::endl;
279 // std::cout << res[0] <<std::endl;
280 
281  //-1 + (a*I1(x))/fc + ((b - c*Cos(3*theta(x)))*Sqrt(J2(x)))/fc
282 
283 
284 
285 
286  sigma.Lodeangle(gradlode,lode);
287 // if(fabs( TPZExtractVal::val(J2) ) < 1.e-6)J2 += 1.e-6;
288 // T sqrtJ2 = sqrt(J2);
289 // T lode1 =-asin( ( T( 3.) * sqrt( T( 3.) ) * J3 ) /( T( 2.) * sqrt(J2*J2*J2) ) )/T( 3.);
290 // cout << " sigma " << endl;
291 // cout << sigma << endl;
292 // cout << " gradlode " << endl;
293 // cout << gradlode << endl;
294 // cout << " lode " << endl;
295 // cout << lode << endl;
296  T temp1;
297  temp1 = ((I1*T(fa))/(T(fcConcrete)));
298 // cout << "temp1" <<endl;
299 // cout << temp1<<endl;
300  T cos3lode = cos(T(3.)*lode);
301  T temp4 = ( T(fb) - ( A * cos3lode ) );
302 // cout << "fb" <<endl;
303 // cout << fb<<endl;
304 
305 // cout << "A" <<endl;
306 // cout << A <<endl;
307 
308 // cout << "cos3lode" <<endl;
309 // cout << cos3lode <<endl;
310 
311 // cout << "temp4" <<endl;
312 // cout << temp4<<endl;
313  T temp2 = ( ( temp4 ) * sqrt(J2) )/(T(fcConcrete));
314  T resp = (temp1 + temp2) - 1.;
315 // std::cout << " PHI " << std::endl;
316 // std::cout << resp <<std::endl;
317 
318  res[0] = resp;
319 
320 
321 }
322 
323 template <class T>
324 void TPZYCWillamWarnke::N(const TPZTensor<T> & sigma,const T & A, TPZVec<TPZTensor<T> > & Ndir, int checkForcedYield) const
325 {
326  TPZTensor<T> gradlode,dJ2,dJ3,dI1,S;
327  T lode,I1,J2,J3;
328 
329  I1 = sigma.I1();
330  J2 = sigma.J2();
331  J3 = sigma.J3();
332  sigma.dJ2(dJ2);
333  sigma.dJ3(dJ3);
334  sigma.S(S);
335  dI1.Identity();
336  sigma.Lodeangle(gradlode,lode);
337 
338 /*
339  T I13=I1/T(3.);
340 
341  T rhoc = -(1./(2.*fb2))*(fb1+ sqrt((fb1*fb1)-T(4.)*T(fb2)*(T(fb0)-I13)));
342  T rhot = -(1./(2.*fa2))*(fa1+sqrt((fa1*fa1)-T(4.)*T(fa2)*(T(fa0)-I13)));
343 
344 
345  T rhotL = T(1.)/(fa1+T(2.)*T(fa2)*rhot);
346  T rhocL = T(1.)/(fb1+T(2.)*T(fb2)*rhoc);
347 
348  T s = T(2.)*rhoc*(rhoc*rhoc-rhot*rhot);
349  T u = T(4.)*(rhoc*rhoc-rhot*rhot)*cos(lode)*cos(lode)+T(5.)*rhot*rhot-T(4.)*rhot*rhoc;
350  T t = rhoc*(T(2.)*rhot-rhoc)*sqrt(u);
351  T v = T(4.)*(rhoc*rhoc)*cos(lode)*cos(lode)+(rhoc-T(2.)*rhot)*(rhoc-T(2.)*rhot);
352 
353  T rho = sqrt(T(2.)*J2);
354  T temp1 = -(rhotL/(T(3.)*v))*(T(4.)*rho*(rhoc-T(2.)*rhot)+ T(8.)* rho * rhot * cos(lode) * cos(lode) - ( T(4.)*rhoc*rhot *cos(lode) - T(2.)*rhoc*sqrt(u)+rhoc*(T(2.)*rhot-rhoc)*(T(1.)/sqrt(u))*(T(4.)*rhot*cos(lode)*cos(lode)-T(5.)*rhot+T(2.)*rhoc)));
355  T temp21 = T(2.)*rho*(rhoc-T(2.)*rhot)+ T(8.) * rho * rhoc * cos(lode) * cos(lode);//rhocL/(T(3.)*v))
356  T temp22 = -( T(2.)*(T(3.)*rhoc*rhoc-rhot*rhot)*cos(lode)+T(2.)*(rhot-rhoc)*sqrt(u)+T(2.)*rhoc*(T(2.)*rhot-rhoc)*(T(1.)/sqrt(u))*((T(2.)*rhoc*cos(lode)*cos(lode))-rhot));
357  T temp2 = (rhocL/(T(3.)*v))*(temp21 + temp22);
358 // T temp2 = (rhocL/(T(3.)*v))*(T(2.)*rho*(rhoc-T(2.)*rhot)+ T(8.) * rho * rhoc * cos(lode) * cos(lode) - ( T(2.)*(T(3.)*rhoc*rhoc-rhot*rhot)*cos(lode)+T(2.)*(rhot-rhoc)*sqrt(u)+T(2.)*rhoc*(T(2.)*rhot-rhoc)*(T(1.)/sqrt(u))*(T(2.)*rhoc*cos(lode)*cos(lode)-rhot));
359  //A0 = GRAD f/I1
360  T A0 = temp1+temp2;
361  T temp3 = (sqrt(T(3.))*(rhoc*rhoc-rhot*rhot))/(v*(T(4.)*cos(lode)*cos(lode)-1)*sqrt(J2*J2*J2));
362  T tt = ((T(2.)*(T(1.)/sqrt(u))))* cos(lode);
363  T ttt = T(4.)*rho*cos(lode);
364  T temp4 = ( ttt - rhoc * ( T(1.) + tt * ( T(2.) * rhot - rhoc ) ) );
365  //A2 = GRAD f/J3
366  T A2 = temp3*temp4;
367  //A1 = GRAD f/J2
368  T A1 = -((T(3.)*J3)/(T(2.)*J2))*A2;
369 
370 
371  dI1.Multiply(A0, T(1.));
372  S.Multiply(A1, T(1.));
373  dJ3.Multiply(A2, T(1.));
374 //
375 // cout << "dI1*A0" <<endl;
376 // cout << dI1 << endl;
377 //
378 // cout << "dJ2*A1" <<endl;
379 // cout << dJ2 << endl;
380 //
381 // cout << "dJ3*A2" <<endl;
382 // cout << dJ3 << endl;
383 
384  dI1.Add(dJ2, T(1));
385  dI1.Add(dJ3, T(1));
386 
387  Ndir[0] = dI1;
388 */
389 
390 /*
391  T sqrtJ2 = sqrt(J2);
392  T lode1 =-asin( ( T( 3.) * sqrt( T( 3.) ) * J3 ) /( T( 2.) * sqrt(J2*J2*J2) ) )/T( 3.);
393  T temp11= -( T(9.*sqrt(3.)) * J3) / ( T(4.)*sqrt(J2*J2*J2*J2*J2) );
394  TPZTensor<T> resultdj2(dJ2);
395  resultdj2.Multiply(temp11,T(1.));
396  T temp2 = ( T(3.)*sqrt( T( 3. ) ) )/ ( T( 2. ) * sqrt( J2 * J2 * J2 ) );
397  TPZTensor<T> resultdj3(dJ3);
398  resultdj3.Multiply(temp2,T(1.));
399  T temp33 = ( T( 27. ) * ( J3 * J3 ) )/( T( 4. ) * J2 * J2 * J2);
400  T temp3 = ( T( 3. ) * sqrt( T( 1. )- temp33 ));
401  TPZTensor<T> RESULT(resultdj2);
402  RESULT.Add(resultdj3, T(1.));
403  RESULT.Multiply( T(1.) / temp3 ,T(-1.));
404 
405 */
406 
407 // (a*DI1)/fc + (DJ2*(b - c*Cos(3*theta(x) ) ) )/(2.*fc*Sqrt(J2(x))) +(3*c*Dtheta*Sqrt(J2(x))*Sin(3*theta(x)))/fc
408  T temp111 = T(fa)/T(fcConcrete);
409  dI1.Identity();
410  dI1.Multiply(temp111,T(1.));
411  T temp222 = ( T(fb) - (A * cos( T(3.) * lode )) )/(T(2.) * T(fcConcrete) * sqrt(J2));
412  dJ2.Multiply(temp222, T(1.));
413 
414  T temp333 = ( T(3.) * A *sqrt(J2) * sin(T(3.) * lode) )/T(fcConcrete);
415  gradlode.Multiply(temp333, T(1.));
416  dI1.Add(dJ2,T(1));
417  dI1.Add(gradlode,T(1));
418  Ndir[0] = dI1;
419 
420 
421 }
422 
423 template <class T>
424 void TPZYCWillamWarnke::H(const TPZTensor<T> & sigma,const T & A, TPZVec<T> & h, int checkForcedYield) const
425 {
426 
427  h[0] = 1.;
428 }
429 
430 
431 #endif//TPZYCWillamWarnke
TPZSavable.h
Contains declaration of the TPZSavable class which defines the interface to save and restore objects ...
TPZYCWillamWarnke::LoadState
void LoadState(TPZFMatrix< REAL > &state)
Definition: TPZYCWillamWarnke.h:137
TPZTensor::Add
void Add(const TPZTensor< T1 > &tensor, const T2 &constant)
Definition: TPZTensor.h:757
TPZYCWillamWarnke::Name
const char * Name() const
Definition: TPZYCWillamWarnke.h:40
TPZTensor::Multiply
void Multiply(const T1 &multipl, const T2 &constant)
Definition: TPZTensor.h:766
TPZYCWillamWarnke::Write
void Write(TPZStream &out, int withclassid=0) const override
Writes this object to the TPZStream buffer. Include the classid if withclassid = true.
Definition: TPZYCWillamWarnke.h:184
TPZTensor::I1
T I1() const
Definition: TPZTensor.h:903
TPZYCWillamWarnke::AlphaMultiplier
void AlphaMultiplier(const T &A, T &multiplier) const
Definition: TPZYCWillamWarnke.h:117
TPZPlasticCriterion.h
as_integer
std::underlying_type< Enumeration >::type as_integer(const Enumeration value)
Definition: pzreal.h:37
TPZTensor::YY
T & YY() const
Definition: TPZTensor.h:578
h
clarg::argBool h("-h", "help message", false)
input
clarg::argString input("-if", "input file", "cube1.txt")
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
TPZYCWillamWarnke::Compute
void Compute(const TPZTensor< T > &sigma, const T &A, TPZVec< T > &res, int checkForcedYield=0) const
Definition: TPZYCWillamWarnke.h:238
TPZYCWillamWarnke::YieldFunction
void YieldFunction(const TPZVec< STATE > &sigma, STATE kprev, TPZVec< STATE > &yield) const override
Definition: TPZYCWillamWarnke.h:214
TPZYCWillamWarnke::N
void N(const TPZTensor< T > &sigma, const T &A, TPZVec< TPZTensor< T > > &Ndir, int checkForcedYield=0) const
Definition: TPZYCWillamWarnke.h:324
sin
sin
Definition: fadfunc.h:63
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZYCWillamWarnke
Implementa a plastificacao do criterio de Von Mises.
Definition: TPZYCWillamWarnke.h:19
TPZYCWillamWarnke::gRefTension
TPZTensor< REAL > gRefTension
Definition: TPZYCWillamWarnke.h:128
pzfmatrix.h
Contains TPZMatrixclass which implements full matrix (using column major representation).
TPZTensor::dJ3
void dJ3(TPZTensor< T > &deriv) const
Definition: TPZTensor.h:956
TPZTensor::Lodeangle
void Lodeangle(TPZTensor< T > &GradLode, T &Lode) const
Definition: TPZTensor.h:1675
TPZTensor::dJ2
void dJ2(TPZTensor< T > &Tangent) const
Definition: TPZTensor.h:935
TPZYCWillamWarnke::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: TPZYCWillamWarnke.cpp:12
TPZYCWillamWarnke::SetForceYield
void SetForceYield(const int forceYield)
Definition: TPZYCWillamWarnke.h:66
TPZYCWillamWarnke::Read
void Read(TPZStream &input, void *context=0) override
read objects from the stream
Definition: TPZYCWillamWarnke.h:199
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
test.res
string res
Definition: test.py:151
TPZTensor::J2
T J2() const
Definition: TPZTensor.h:927
TPZTensor::Identity
void Identity()
Definition: TPZTensor.h:705
TPZYCWillamWarnke::Print
void Print(std::ostream &out) const override
Definition: TPZYCWillamWarnke.h:45
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
TPZYCWillamWarnke::GetNYield
virtual int GetNYield() const override
Definition: TPZYCWillamWarnke.h:222
TPZTensor::Eigenvalue
void Eigenvalue(TPZTensor< T > &eigenval, TPZTensor< T > &dSigma1, TPZTensor< T > &dSigma2, TPZTensor< T > &dSigma3) const
Definition: TPZTensor.h:1744
TPZTensor::XX
T & XX() const
Definition: TPZTensor.h:566
TPZTensor::ZZ
T & ZZ() const
Definition: TPZTensor.h:586
TPZYCWillamWarnke::SetUp
void SetUp(const REAL a, const REAL b, const REAL fConcrete)
Definition: TPZYCWillamWarnke.h:53
TPZYCWillamWarnke::SetYieldStatusMode
void SetYieldStatusMode(const TPZTensor< REAL > &sigma, const REAL &A)
Definition: TPZYCWillamWarnke.h:78
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZTensor.h
TPZYCWillamWarnke::Residual
void Residual(TPZFMatrix< REAL > &res, int icase)
Definition: TPZYCWillamWarnke.h:164
cos
TPZFlopCounter cos(const TPZFlopCounter &orig)
Returns the cosine in radians and increments the counter of the Cosine.
Definition: pzreal.h:514
TPZYCWillamWarnke::H
void H(const TPZTensor< T > &sigma, const T &A, TPZVec< T > &h, int checkForcedYield=0) const
Definition: TPZYCWillamWarnke.h:424
TPZYCWillamWarnke::ComputeTangent
void ComputeTangent(TPZFMatrix< REAL > &tangent, TPZVec< REAL > &coefs, int icase)
Definition: TPZYCWillamWarnke.h:143
TPZYCWillamWarnke::TPZYCWillamWarnke
TPZYCWillamWarnke(const TPZYCWillamWarnke &source)
Definition: TPZYCWillamWarnke.h:25
TPZTensor::S
void S(TPZTensor< T > &s) const
Definition: TPZTensor.h:894
TPZTensor::J3
T J3() const
Definition: TPZTensor.h:946
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