NeoPZ
pzviscoelastic.cpp
Go to the documentation of this file.
1 
6 #include "pzviscoelastic.h"
7 
8 TPZViscoelastic::TPZViscoelastic() :
9 TPZRegisterClassId(&TPZViscoelastic::ClassId),
10 TPZMatWithMem<TPZFMatrix<STATE>, TPZElasticity3D>(), fAlpha(-1), fDeltaT(-1), fLambdaV(-1.), fmuV(-1.)
11 {
12 
13 }
14 
15 TPZViscoelastic::TPZViscoelastic(int id) :
16 TPZRegisterClassId(&TPZViscoelastic::ClassId),
17 TPZMatWithMem<TPZFMatrix<STATE>, TPZElasticity3D>(id), fAlpha(-1), fDeltaT(-1), fLambdaV(-1.), fmuV(-1.)
18 {
19 
20 }
21 
22 TPZViscoelastic::TPZViscoelastic(int id,STATE ElaE,STATE poissonE, STATE lambdaV, STATE muV, STATE alpha, STATE deltaT, TPZVec <STATE> &force): TPZRegisterClassId(&TPZViscoelastic::ClassId),
23 TPZMatWithMem<TPZFMatrix<STATE>, TPZElasticity3D>(id), fAlpha(alpha), fDeltaT(deltaT), fLambdaV(lambdaV), fmuV(muV)
24 
25 {
26  STATE lambdaE = (poissonE * ElaE)/((1+poissonE)*(1-2*poissonE));
27  STATE muE = ElaE/(2*(1+poissonE));
28  STATE lambdaVE = lambdaE-(fDeltaT*lambdaV)/(1+fAlpha*fDeltaT);
29  STATE muVE = muE -(fDeltaT*muV)/(1+fAlpha*fDeltaT);
30  STATE ElaVE = muVE*(3*lambdaVE+2*muVE)/(lambdaVE+muVE);
31  STATE PoissonVE = lambdaVE/(2*(lambdaVE+muVE));
32  if (lambdaVE < 0 || muVE < 0)
33  {
34  PZError << "lambdaVE and muVE must be positive. Check your constants values\n";
35  DebugStop();
36  }
37  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::SetMaterialDataHook(ElaVE, PoissonVE); // Creating an elastic material with the viscoelastic properties
38  SetForce(force);
39  //SetC(); already set in SetMaterialDataHook from elastic material
40 }
41 
42 void TPZViscoelastic::SetMaterialDataHooke(STATE ElaE, STATE poissonE, STATE ElaV, STATE poissonV, STATE alpha, STATE deltaT, TPZVec <STATE> &force)
43 {
44  fAlpha = alpha;
45  fDeltaT = deltaT;
46  STATE lambdaE = (poissonE * ElaE)/((1+poissonE)*(1-2*poissonE));
47  STATE muE = ElaE/(2*(1+poissonE));
48  STATE lambdaV = (poissonV * ElaV)/((1+poissonV)*(1-2*poissonV));
49  STATE muV = ElaV/(2*(1+poissonV));
50  fLambdaV = lambdaV;
51  fmuV = muV;
52  STATE lambdaVE = lambdaE-(fDeltaT*lambdaV)/(1+fAlpha*fDeltaT);
53  STATE muVE = muE -(fDeltaT*muV)/(1+fAlpha*fDeltaT);
54  STATE ElaVE = muVE*(3*lambdaVE+2*muVE)/(lambdaVE+muVE);
55  STATE PoissonVE = lambdaVE/(2*(lambdaVE+muVE));
56  if (lambdaVE < 0 || muVE < 0)
57  {
58  PZError << "lambdaVE and muVE must be positive. Check your constants values\n";
59  DebugStop();
60  }
61  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::SetMaterialDataHook(ElaVE, PoissonVE); // Creating an elastic material with the viscoelastic properties
62  SetForce(force);
63  //SetC(); already set in SetMaterialDataHook from elastic material
64 }
65 
66 void TPZViscoelastic::Contribute(TPZMaterialData &data,REAL weight,TPZFMatrix<STATE> &ek,TPZFMatrix<STATE> &ef)
67 {
68  int nstate = this->NStateVariables();
69  int in;
70  STATE val;
71  if(fUpdateMem != 0)
72  {
73  UpdateQsi(data);
74  }
75  TPZFMatrix<REAL> dphi = data.dphix;
76  TPZFMatrix<REAL> phi = data.phi;
77  const int phr = phi.Rows();
78  int index = data.intGlobPtIndex;
79 
80  TPZFNMatrix<6,STATE> qsi(6,1);
81  int rows = this->MemItem(index).Rows();
82  if (rows < 6 || rows > 6)
83  {
84  DebugStop(); //deve inicializar o qsi pelo SetDefaultMemory
85  }
86  else
87  {
88  qsi = this->MemItem(index);
89  }
90 
91  for(in = 0; in < phr; in++)
92  {
93  //in: test function index. First equation
94  val = 0.;
95  val -= qsi(_XX_,0) * dphi(0,in); // |
96  val -= qsi(_XY_,0) * dphi(1,in); // fk
97  val -= qsi(_XZ_,0) * dphi(2,in); // |
98  val *= 1./(1+fAlpha*fDeltaT);
99  ef(in*nstate+0,0) += weight * val;
100 
101  //Second equation: fb and fk
102  val = 0.;
103  val -= qsi(_XY_,0) * dphi(0,in); // |
104  val -= qsi(_YY_,0) * dphi(1,in); // fk
105  val -= qsi(_YZ_,0) * dphi(2,in); // |
106  val *= 1./(1+fAlpha*fDeltaT);
107  ef(in*nstate+1,0) += weight * val;
108 
109  //third equation: fb and fk
110  val = 0.;
111  val -= qsi(_XZ_,0) * dphi(0,in); // |
112  val -= qsi(_YZ_,0) * dphi(1,in); // fk
113  val -= qsi(_ZZ_,0) * dphi(2,in); // |
114  val *= 1./(1+fAlpha*fDeltaT);
115  ef(in*nstate+2,0) += weight * val;
116  }
117  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::Contribute(data,weight,ek,ef);
118 }
119 
120 void TPZViscoelastic::UpdateQsi(TPZMaterialData &data)
121 {
122  TPZFNMatrix<6,STATE> qsi;
123  TPZFNMatrix<6,STATE> qsin1(6,1,0.);
124  TPZFNMatrix<9,STATE> DSolXYZ(3,3,0.);
125  TPZFNMatrix<6,STATE> Strain(6,1);
126  int index = data.intGlobPtIndex;
127  int numbersol = data.dsol.size();
128  if (numbersol != 1) {
129  DebugStop();
130  }
131 
132  DSolXYZ = data.dsol[0];
133  qsi = this->MemItem(index);
134 
135  Strain.Redim(6,1);
136  Strain(_XX_,0) = DSolXYZ(0,0);
137  Strain(_YY_,0) = DSolXYZ(1,1);
138  Strain(_ZZ_,0) = DSolXYZ(2,2);
139  Strain(_XY_,0) = 0.5 * ( DSolXYZ(1,0) + DSolXYZ(0,1) );
140  Strain(_XZ_,0) = 0.5 * ( DSolXYZ(2,0) + DSolXYZ(0,2) );
141  Strain(_YZ_,0) = 0.5 * ( DSolXYZ(2,1) + DSolXYZ(1,2) );
142 
143  STATE tr;
144  tr = Strain(_XX_,0)+Strain(_YY_,0)+Strain(_ZZ_,0);
145 
146  qsin1(_XX_,0) = (fDeltaT*(-(tr)*fLambdaV - 2*Strain(_XX_,0)*fmuV) + qsi(_XX_,0))/(1 + fAlpha*fDeltaT);
147  qsin1(_YY_,0) = (fDeltaT*(-(tr)*fLambdaV - 2*Strain(_YY_,0)*fmuV) + qsi(_YY_,0))/(1 + fAlpha*fDeltaT);
148  qsin1(_ZZ_,0) = (fDeltaT*(-(tr)*fLambdaV - 2*Strain(_ZZ_,0)*fmuV) + qsi(_ZZ_,0))/(1 + fAlpha*fDeltaT);
149  qsin1(_XY_,0) = (-2*fDeltaT*Strain(_XY_,0)*fmuV + qsi(_XY_,0))/(1 + fAlpha*fDeltaT);
150  qsin1(_XZ_,0) = (-2*fDeltaT*Strain(_XZ_,0)*fmuV + qsi(_XZ_,0))/(1 + fAlpha*fDeltaT);
151  qsin1(_YZ_,0) = (-2*fDeltaT*Strain(_YZ_,0)*fmuV + qsi(_YZ_,0))/(1 + fAlpha*fDeltaT);
152 
153  MemItem(index) = qsin1;
154 }
155 
156 
157 int TPZViscoelastic::VariableIndex(const std::string &name)
158 {
159  if(!strcmp("Displacement",name.c_str())) return TPZElasticity3D::EDisplacement;
160  if(!strcmp("state",name.c_str())) return TPZElasticity3D::EDisplacement;
161  if(!strcmp("DisplacementX",name.c_str())) return TPZElasticity3D::EDisplacementX;
162  if(!strcmp("DisplacementY",name.c_str())) return TPZElasticity3D::EDisplacementY;
163  if(!strcmp("DisplacementZ",name.c_str())) return TPZElasticity3D::EDisplacementZ;
164  if(!strcmp("PrincipalStrain", name.c_str())) return TPZElasticity3D::EPrincipalStrain;
165  if(!strcmp("ViscoStressX",name.c_str())) return TPZViscoelastic::EViscoStressX;
166  if(!strcmp("ViscoStressY",name.c_str())) return TPZViscoelastic::EViscoStressY;
167  if(!strcmp("ViscoStressZ",name.c_str())) return TPZViscoelastic::EViscoStressZ;
168  return -1;
169 }
170 
171 int TPZViscoelastic::NSolutionVariables(int var)
172 {
173  if(var == TPZElasticity3D::EDisplacement) return 3;
174  if(var == TPZElasticity3D::EDisplacementX) return 1;
175  if(var == TPZElasticity3D::EDisplacementY) return 1;
176  if(var == TPZElasticity3D::EDisplacementZ) return 1;
177  if(var == TPZElasticity3D::EPrincipalStrain) return 3;
178  if(var == TPZViscoelastic::EViscoStressX) return 1;
179  if(var == TPZViscoelastic::EViscoStressY) return 1;
180  if(var == TPZViscoelastic::EViscoStressZ) return 1;
181  PZError << "TPZViscoelastic::NSolutionVariables Error\n";
182  return -1;
183 }
184 
185 void TPZViscoelastic::ComputeStressTensor(TPZFMatrix<STATE> &Stress, TPZMaterialData &data) const
186 {
187  int numbersol = data.dsol.size();
188  if (numbersol != 1) {
189  DebugStop();
190  }
191 
192  TPZFNMatrix<9,STATE> Dsol = data.dsol[0];
193  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::ComputeStressTensor(Stress,Dsol);
194  TPZFNMatrix<6,STATE> qsi;
195  const int index = data.intGlobPtIndex;
196  qsi = *MemItem(index);
197  const REAL denominador = 1 + fAlpha*fDeltaT;
198  qsi(_XX_,0)*= 1/(1+denominador);
199  qsi(_XY_,0)*= 1/(1+denominador);
200  qsi(_XZ_,0)*= 1/(1+denominador);
201  qsi(_YY_,0)*= 1/(1+denominador);
202  qsi(_YZ_,0)*= 1/(1+denominador);
203  qsi(_ZZ_,0)*= 1/(1+denominador);
204 
205  Stress(0,0) += qsi(_XX_,0);
206  Stress(1,1) += qsi(_YY_,0);
207  Stress(2,2) += qsi(_ZZ_,0);
208  Stress(0,1) += qsi(_XY_,0);
209  Stress(0,2) += qsi(_XZ_,0);
210  Stress(1,0) += qsi(_XY_,0);
211  Stress(1,2) += qsi(_YZ_,0);
212  Stress(2,0) += qsi(_XZ_,0);
213  Stress(2,1) += qsi(_YZ_,0);
214 }
215 
216 void TPZViscoelastic::Solution(TPZMaterialData &data, int var, TPZVec<STATE> &Solout)
217 {
218  int numbersol = data.sol.size();
219  if (numbersol != 1) {
220  DebugStop();
221  }
222 
223  if(var == TPZElasticity3D::EDisplacement){
224  int i;
225  for(i = 0; i < 3; i++){
226  TPZVec<STATE> Sol(data.sol[0]);
227  Solout[i] = Sol[i];
228  }//for
229  return;
230  }//TPZElasticity3D::EDisplacement
231 
232  if(var == TPZElasticity3D::EDisplacementX){
233  // int i;
234  TPZVec<STATE> Sol(data.sol[0]);
235  Solout[0] = Sol[0];
236  return;
237  }//TPZElasticity3D::EDisplacementX
238 
239  if(var == TPZElasticity3D::EDisplacementY){
240  // int i;
241  TPZVec<STATE> Sol(data.sol[0]);
242  Solout[0] = Sol[1];
243  return;
244  }//TPZElasticity3D::EDisplacementY
245 
246  if(var == TPZElasticity3D::EDisplacementZ){
247  // int i;
248  TPZVec<STATE> Sol(data.sol[0]);
249  Solout[0] = Sol[2];
250  return;
251  }//TPZElasticity3D::EDisplacementZ
252 
253 
254  if(var == TPZElasticity3D::EPrincipalStrain){
255  TPZFNMatrix<9,STATE> StrainTensor(3,3);
256  TPZFMatrix<STATE> DSol(data.dsol[0]);
257  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::ComputeStrainTensor(StrainTensor, DSol);
258  int64_t numiterations = 1000;
259  REAL tol = TPZElasticity3D::gTolerance;
260  TPZManVector<STATE,3> eigv(3,0.);
261  bool result;
262  result = StrainTensor.SolveEigenvaluesJacobi(numiterations, tol, &eigv);
263  for (int i=0; i<3; i++) {
264  Solout[i] = eigv[i];
265  }
266 #ifdef PZDEBUG
267  if (result == false){
268  PZError << __PRETTY_FUNCTION__ << " - ERROR! - result = false - numiterations = " << numiterations << " - tol = " << tol << std::endl;
269  }
270 #endif
271  }//TPZElasticity3D::EPrincipalStrain
272 
273 
274  if(var == TPZViscoelastic::EViscoStressX){
275  TPZFNMatrix<9,STATE> Stress(3,3);
276  this->ComputeStressTensor(Stress, data);
277  Solout[0] = Stress(0,0);
278  return;
279  }
280  if(var == TPZViscoelastic::EViscoStressY){
281  TPZFMatrix<STATE> Stress(3,3);
282  this->ComputeStressTensor(Stress, data);
283  Solout[0] = Stress(1,1);
284  return;
285  }
286  if(var == TPZViscoelastic::EViscoStressZ){
287  TPZFMatrix<STATE> Stress(3,3);
288  this->ComputeStressTensor(Stress, data);
289  Solout[0] = Stress(2,2);
290  return;
291  }
292 }
293 
294 void TPZViscoelastic::FillDataRequirements(TPZMaterialData &data){
295 
296  TPZMaterial::FillDataRequirements(data);
297  data.fNeedsSol = true;
298 }
299 
301 void TPZViscoelastic::Write(TPZStream &buf, int withclassid) const
302 {
303  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::Write(buf,withclassid);
304  buf.Write(&fAlpha, 1);
305  buf.Write(&fDeltaT, 1);
306  buf.Write(&fLambdaV, 1);
307  buf.Write(&fmuV, 1);
308 }
309 
311 void TPZViscoelastic::Read(TPZStream &buf, void *context)
312 {
313  TPZMatWithMem<TPZFMatrix<STATE>,TPZElasticity3D>::Read(buf,context);
314  buf.Read(&fAlpha, 1);
315  buf.Read(&fDeltaT, 1);
316  buf.Read(&fLambdaV, 1);
317  buf.Read(&fmuV, 1);
318 
319 }
320 
321 int TPZViscoelastic::ClassId() const{
322  return Hash("TPZViscoelastic") ^ TPZMatWithMem<TPZFMatrix<STATE>, TPZElasticity3D>::ClassId() << 1;
323 }
324 
325 #ifndef BORLAND
326 template class TPZRestoreClass<TPZViscoelastic>;
327 #endif
_XZ_
#define _XZ_
Definition: TPZTensor.h:29
TPZMatrix::SolveEigenvaluesJacobi
virtual bool SolveEigenvaluesJacobi(int64_t &numiterations, REAL &tol, TPZVec< TVar > *Sort=0)
Transforms this matrix in a diagonal matrix, where the diagonal values are its eigenvalues. This method is efficient only for small matrices.
Definition: pzmatrix.cpp:1727
TPZViscoelastic::ComputeStressTensor
virtual void ComputeStressTensor(TPZFMatrix< STATE > &Stress, TPZMaterialData &data) const override
Definition: pzviscoelastic.cpp:185
TPZElasticity3D
This class implements a 3D isotropic elasticity material.
Definition: pzelast3d.h:21
TPZMaterialData::dsol
TPZGradSolVec dsol
vector of the derivatives of the solution at the integration point
Definition: pzmaterialdata.h:79
TPZViscoelastic::Write
void Write(TPZStream &buf, int withclassid) const override
Saves the element data to a stream.
Definition: pzviscoelastic.cpp:301
TPZElasticity3D::gTolerance
static STATE gTolerance
Definition: pzelast3d.h:304
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZMatWithMem
Implements an abstract class implementing the memory features.
Definition: TPZMatWithMem.h:23
TPZMaterialData::phi
TPZFNMatrix< 220, REAL > phi
vector of shapefunctions (format is dependent on the value of shapetype)
Definition: pzmaterialdata.h:55
TPZViscoelastic::Read
void Read(TPZStream &buf, void *context) override
Reads the element data from a stream.
Definition: pzviscoelastic.cpp:311
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
TPZVec::size
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
TPZMatWithMem< TPZFMatrix< STATE >, TPZElasticity3D >::fUpdateMem
bool fUpdateMem
Flag to indicate whether the memory data are to be updated in an assemble loop.
Definition: TPZMatWithMem.h:115
TPZElasticity3D::SetForce
void SetForce(TPZVec< STATE > force)
Definition: pzelast3d.h:223
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
TPZStream::Write
virtual void Write(const bool val)
Definition: TPZStream.cpp:8
TPZMaterial::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data)
Fill material data parameter with necessary requirements for the.
Definition: TPZMaterial.cpp:81
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
_XX_
#define _XX_
Definition: TPZTensor.h:27
_YZ_
#define _YZ_
Definition: TPZTensor.h:31
_XY_
#define _XY_
Definition: TPZTensor.h:28
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
TPZElasticity3D::NStateVariables
int NStateVariables() const override
Number of state variables.
Definition: pzelast3d.h:73
TPZMaterialData::intGlobPtIndex
int intGlobPtIndex
global point index
Definition: pzmaterialdata.h:115
_YY_
#define _YY_
Definition: TPZTensor.h:30
TPZFMatrix
Full matrix class. Matrix.
Definition: pzfmatrix.h:32
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZElasticity3D::ComputeStrainTensor
void ComputeStrainTensor(TPZFMatrix< STATE > &Strain, TPZFMatrix< STATE > &DSol) const
Definition: pzelast3d.cpp:1115
_ZZ_
#define _ZZ_
Definition: TPZTensor.h:32
TPZViscoelastic::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: pzviscoelastic.cpp:321
TPZElasticity3D::SetMaterialDataHook
void SetMaterialDataHook(REAL Ela, REAL poisson)
Definition: pzelast3d.h:202
TPZViscoelastic::fDeltaT
STATE fDeltaT
the stability coeficient alpha
Definition: pzviscoelastic.h:112
TPZViscoelastic::fAlpha
STATE fAlpha
the stability coeficient alpha
Definition: pzviscoelastic.h:107
TPZViscoelastic::TPZViscoelastic
TPZViscoelastic()
Empty constructor.
Definition: pzviscoelastic.cpp:8
TPZViscoelastic::NSolutionVariables
virtual int NSolutionVariables(int var) override
Number of data of variable var.
Definition: pzviscoelastic.cpp:171
TPZViscoelastic
This class implements an isotropic viscoelasticity material.
Definition: pzviscoelastic.h:36
TPZViscoelastic::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the Contribute method.
Definition: pzviscoelastic.cpp:294
pzviscoelastic.h
Contains the TPZViscoelastic class which implements an isotropic viscoelasticity material.
TPZViscoelastic::fLambdaV
STATE fLambdaV
viscoelasticity coeficients
Definition: pzviscoelastic.h:117
TPZViscoelastic::UpdateQsi
void UpdateQsi(TPZMaterialData &data)
Definition: pzviscoelastic.cpp:120
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
TPZViscoelastic::VariableIndex
virtual int VariableIndex(const std::string &name) override
Returns index of post-processing variable.
Definition: pzviscoelastic.cpp:157
TPZViscoelastic::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< REAL > &ek, TPZFMatrix< REAL > &ef)
Definition: pzviscoelastic.cpp:36
TPZMatWithMem< TPZFMatrix< STATE >, TPZElasticity3D >::MemItem
virtual TPZFMatrix< STATE > & MemItem(const int i) const
Definition: TPZMatWithMem.h:172
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
TPZViscoelastic::SetMaterialDataHooke
void SetMaterialDataHooke(STATE ElaE, STATE poissonE, STATE ElaV, STATE poissonV, STATE alpha, STATE deltaT, TPZVec< STATE > &force)
Set material Data with hooke constants.
Definition: pzviscoelastic.cpp:42
TPZRestoreClass
Implements an interface to register a class id and a restore function. Persistence.
Definition: TPZSavable.h:150
TPZViscoelastic::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: pzviscoelastic.cpp:216
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:01 for NeoPZ by   doxygen 1.8.13