NeoPZ
TPZElast3Dnlinear.cpp
Go to the documentation of this file.
1 //
2 // TPZElast3Dnlinear.cpp
3 // PZ
4 //
5 // Created by Cesar Lucci on 22/10/13.
6 //
7 //
8 
9 #include "TPZElast3Dnlinear.h"
10 #include "pzbndcond.h"
11 
12 TPZElast3Dnlinear::TPZElast3Dnlinear() : TPZRegisterClassId(&TPZElast3Dnlinear::ClassId),
13 TPZElasticity3D()
14 {
15 
16 }
17 
18 TPZElast3Dnlinear::TPZElast3Dnlinear(int nummat, STATE E, STATE poisson, TPZVec<STATE> &force,
19  STATE preStressXX, STATE preStressYY, STATE preStressZZ) :
20 TPZRegisterClassId(&TPZElast3Dnlinear::ClassId),
21 TPZElasticity3D(nummat, E, poisson, force, preStressXX, preStressYY, preStressZZ)
22 {
23 
24 }
25 
26 TPZElast3Dnlinear::~TPZElast3Dnlinear()
27 {
28 
29 }
30 
31 void TPZElast3Dnlinear::Contribute(TPZMaterialData &data,
32  REAL weight,
33  TPZFMatrix<STATE> &ef)
34 {
35  DebugStop();//Not implemented!!!
36 }
37 
38 void TPZElast3Dnlinear::Contribute(TPZMaterialData &data,
39  REAL weight,
40  TPZFMatrix<STATE> &ek,
41  TPZFMatrix<STATE> &ef)
42 {
43 #ifdef PZDEBUG
44  if(ek.Rows() != ef.Rows())
45  {
46  std::cout << "\n\n" << "ek and ef should have same number of rows!!!\n";
47  std::cout << "See " << __PRETTY_FUNCTION__ << "\n\n";
48  DebugStop();
49  }
50 
51  if(ef.Cols() != this->fNumLoadCases)
52  {
53  std::cout << "\n\n" << "ef should have fNumLoadCases equals to its NCols()!!!\n";
54  std::cout << "See " << __PRETTY_FUNCTION__ << "\n\n";
55  DebugStop();
56  }
57 #endif
58 
59  ContributeVecShapeAux(data, weight, ek, ef);
60 }
61 
62 void TPZElast3Dnlinear::ContributeBC(TPZMaterialData &data,
63  REAL weight,
64  TPZFMatrix<STATE> &ek,
65  TPZFMatrix<STATE> &ef,
66  TPZBndCond &bc)
67 {
68 #ifdef PZDEBUG
69  if(ek.Rows() != ef.Rows())
70  {
71  std::cout << "\n\n" << "ek and ef should have same number of rows!!!\n";
72  std::cout << "See " << __PRETTY_FUNCTION__ << "\n\n";
73  DebugStop();
74  }
75 
76  if(ef.Cols() != this->fNumLoadCases)
77  {
78  std::cout << "\n\n" << "ef should have fNumLoadCases equals to its NCols()!!!\n";
79  std::cout << "See " << __PRETTY_FUNCTION__ << "\n\n";
80  DebugStop();
81  }
82 #endif
83 
84  ContributeVecShapeBCAux(data, weight, ek, ef, bc);
85 }
86 
87 int TPZElast3Dnlinear::ClassId() const{
88  return Hash("TPZElast3Dnlinear") ^ TPZElasticity3D::ClassId() << 1;
89 }
90 
91 void TPZElast3Dnlinear::FillDataRequirements(TPZMaterialData &data)
92 {
93  TPZMaterial::FillDataRequirements(data);
94  data.fNeedsSol = true;
95 }
96 
97 
98 //--------------------------------------------------------------------------------
99 
100 
101 void TPZElast3Dnlinear::ContributeVecShapeAux(TPZMaterialData &data,
102  REAL weight,
103  TPZFMatrix<STATE> &ek,
104  TPZFMatrix<STATE> &ef)
105 {
106  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
107  if(shapetype != data.EVecShape)
108  {
109  DebugStop();
110  }
111 
112  TPZFMatrix<REAL> & phi = data.phi;
113  TPZFMatrix<REAL> & dphi = data.dphix;
114  TPZFMatrix<STATE> & dsol = data.dsol[0];
115 
116  int phc = phi.Cols();
117  int efc = ef.Cols();
118 
119  if(fForcingFunction)
120  {
121  TPZManVector<STATE> res(3);
122  fForcingFunction->Execute(data.x,res);
123  fForce[0] = res[0];
124  fForce[1] = res[1];
125  fForce[2] = res[2];
126  }
127 
128  REAL dvxdx, dvxdy, dvxdz;
129  REAL dvydx, dvydy, dvydz;
130  REAL dvzdx, dvzdy, dvzdz;
131 
132  REAL duxdx, duxdy, duxdz;
133  REAL duydx, duydy, duydz;
134  REAL duzdx, duzdy, duzdz;
135 
136  /*
137  * Plain strain materials values
138  */
139  REAL lambda = fE*fPoisson/((1.+fPoisson)*(1.-2.*fPoisson));
140  REAL mu = fE/(2.*(1.+fPoisson));
141 
142  for( int in = 0; in < phc; in++ )
143  {
144  //x
145  dvxdx = dphi(0,in);
146  dvxdy = dphi(1,in);
147  dvxdz = dphi(2,in);
148 
149  //y
150  dvydx = dphi(3,in);
151  dvydy = dphi(4,in);
152  dvydz = dphi(5,in);
153 
154  //z
155  dvzdx = dphi(6,in);
156  dvzdy = dphi(7,in);
157  dvzdz = dphi(8,in);
158 
159  for(int col = 0; col < efc; col++)
160  {
161  ef(in,col) += weight*( fForce[0] * phi(0, in) +
162  fForce[1] * phi(1, in) +
163  fForce[2] * phi(2, in) -
164  dvxdx * fPreStress[0] -
165  dvydy * fPreStress[1] -
166  dvzdz * fPreStress[2] );
167 
168  //x
169  duxdx = dsol(0,0);
170  duxdy = dsol(1,0);
171  duxdz = dsol(2,0);
172 
173  //y
174  duydx = dsol(0,1);
175  duydy = dsol(1,1);
176  duydz = dsol(2,1);
177 
178  //z
179  duzdx = dsol(0,2);
180  duzdy = dsol(1,2);
181  duzdz = dsol(2,2);
182 
183  REAL eq1 = duydy*dvxdx*lambda + duzdz*dvxdx*lambda + duxdy*dvydx*mu +
184  duydx*dvydx*mu + duxdz*dvzdx*mu + duzdx*dvzdx*mu +
185  duxdx*dvxdx*(lambda + 2.*mu);
186 
187  REAL eq2 = duxdx*dvydy*lambda + duzdz*dvydy*lambda + duxdy*dvxdy*mu +
188  duydx*dvxdy*mu + duydz*dvzdy*mu + duzdy*dvzdy*mu +
189  duydy*dvydy*(lambda + 2.*mu);
190 
191  REAL eq3 = duxdx*dvzdz*lambda + duydy*dvzdz*lambda + duxdz*dvxdz*mu +
192  duzdx*dvxdz*mu + duydz*dvydz*mu + duzdy*dvydz*mu +
193  duzdz*dvzdz*(lambda + 2.*mu);
194 
195  ef(in,col) -= weight * (eq1 + eq2 + eq3);
196  }
197  for( int jn = 0; jn < phc; jn++ )
198  {
199  //x
200  duxdx = dphi(0,jn);
201  duxdy = dphi(1,jn);
202  duxdz = dphi(2,jn);
203 
204  //y
205  duydx = dphi(3,jn);
206  duydy = dphi(4,jn);
207  duydz = dphi(5,jn);
208 
209  //z
210  duzdx = dphi(6,jn);
211  duzdy = dphi(7,jn);
212  duzdz = dphi(8,jn);
213 
214  REAL eq1 = duydy*dvxdx*lambda + duzdz*dvxdx*lambda + duxdy*dvydx*mu +
215  duydx*dvydx*mu + duxdz*dvzdx*mu + duzdx*dvzdx*mu +
216  duxdx*dvxdx*(lambda + 2.*mu);
217 
218  REAL eq2 = duxdx*dvydy*lambda + duzdz*dvydy*lambda + duxdy*dvxdy*mu +
219  duydx*dvxdy*mu + duydz*dvzdy*mu + duzdy*dvzdy*mu +
220  duydy*dvydy*(lambda + 2.*mu);
221 
222  REAL eq3 = duxdx*dvzdz*lambda + duydy*dvzdz*lambda + duxdz*dvxdz*mu +
223  duzdx*dvxdz*mu + duydz*dvydz*mu + duzdy*dvydz*mu +
224  duzdz*dvzdz*(lambda + 2.*mu);
225 
226  ek(in,jn) += weight * (eq1 + eq2 + eq3);
227  }
228  }
229 }
230 
231 
232 void TPZElast3Dnlinear::ContributeVecShapeBCAux(TPZMaterialData &data,
233  REAL weight,
234  TPZFMatrix<STATE> &ek,
235  TPZFMatrix<STATE> &ef,
236  TPZBndCond &bc)
237 {
238  TPZMaterialData::MShapeFunctionType shapetype = data.fShapeType;
239  if(shapetype != data.EVecShape)
240  {
241  DebugStop();
242  }
243 
244  TPZFMatrix<REAL> & phi = data.phi;
245  TPZManVector<STATE,3> sol = data.sol[0];
246 
247  const REAL BIGNUMBER = TPZMaterial::gBigNumber;
248 
249  int phc = phi.Cols();
250  short in,jn;
251 
252  switch (bc.Type())
253  {
254  case 0:// Dirichlet condition
255  {
256  for(in = 0 ; in < phc; in++)
257  {
258  for(int il = 0; il < fNumLoadCases; il++)
259  {
260  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
261  ef(in,il) += weight * BIGNUMBER * ( v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in) + v2(2,il)*phi(2,in) );
262 
263  ef(in,il) -= weight * BIGNUMBER * ( sol[0]*phi(0,in) + sol[1]*phi(1,in) + sol[2]*phi(2,in) );
264  }
265  for(jn = 0 ; jn < phc; jn++)
266  {
267  ek(in,jn) += weight * BIGNUMBER * ( phi(0,jn)*phi(0,in) + phi(1,jn)*phi(1,in) + phi(2,jn)*phi(2,in) );
268  }
269  }
270  break;
271  }
272  case 1:// Neumann condition
273  {
274  for(in = 0; in < phc; in++)
275  {
276  for(int il = 0; il < fNumLoadCases; il++)
277  {
278  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
279  ef(in,il) += weight * ( v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in) + v2(2,il)*phi(2,in) );
280  }
281  //... continua no TPZPlaneFractCouplingMat
282  }
283  break;
284  }
285  case 2:// condicao mista
286  {
287  DebugStop();
288  for(in = 0 ; in < phc; in++)
289  {
290  for(int il = 0; il < fNumLoadCases; il++)
291  {
292  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
293  ef(in,il) += weight * ( v2(0,il)*phi(0,in) + v2(1,il)*phi(1,in) + v2(2,il)*phi(2,in) );
294 
295  ef(in,il) -= ( bc.Val1()(0,0)*phi(0,in)*sol[0]*weight
296 
297  + bc.Val1()(1,0)*phi(1,in)*sol[0]*weight
298 
299  + bc.Val1()(2,0)*phi(2,in)*sol[0]*weight
300 
301 
302  + bc.Val1()(0,1)*phi(0,in)*sol[1]*weight
303 
304  + bc.Val1()(1,1)*phi(1,in)*sol[1]*weight
305 
306  + bc.Val1()(2,1)*phi(2,in)*sol[1]*weight
307 
308 
309  + bc.Val1()(0,2)*phi(0,in)*sol[2]*weight
310 
311  + bc.Val1()(1,2)*phi(1,in)*sol[2]*weight
312 
313  + bc.Val1()(2,2)*phi(2,in)*sol[2]*weight );
314  }
315 
316  for(jn = 0; jn <phc; jn++)
317  {
318 
319  ek(in,jn) += bc.Val1()(0,0)*phi(0,in)*phi(0,jn)*weight
320 
321  + bc.Val1()(1,0)*phi(1,in)*phi(0,jn)*weight
322 
323  + bc.Val1()(2,0)*phi(2,in)*phi(0,jn)*weight
324 
325 
326  + bc.Val1()(0,1)*phi(0,in)*phi(1,jn)*weight
327 
328  + bc.Val1()(1,1)*phi(1,in)*phi(1,jn)*weight
329 
330  + bc.Val1()(2,1)*phi(2,in)*phi(1,jn)*weight
331 
332 
333  + bc.Val1()(0,2)*phi(0,in)*phi(2,jn)*weight
334 
335  + bc.Val1()(1,2)*phi(1,in)*phi(2,jn)*weight
336 
337  + bc.Val1()(2,2)*phi(2,in)*phi(2,jn)*weight;
338  }
339  }
340  break;
341  }
342  case 3: // Directional Null Dirichlet - displacement is set to null in the non-null vector component direction
343  {
344  for(in = 0 ; in < phc; in++)
345  {
346  for(int il = 0; il < fNumLoadCases; il++)
347  {
348  TPZFNMatrix<2,STATE> v2 = bc.Val2(il);
349  ef(in,il) -= weight * BIGNUMBER * ( v2(0,il)*phi(0,in)*sol[0] + v2(1,il)*phi(1,in)*sol[1] + v2(2,il)*phi(2,in)*sol[2] );
350  for(jn = 0 ; jn < phc; jn++)
351  {
352  ek(in,jn) += weight * BIGNUMBER * ( v2(0,il)*phi(0,in)*phi(0,jn) + v2(1,il)*phi(1,in)*phi(1,jn) + v2(2,il)*phi(2,in)*phi(2,jn) );
353  }
354  }
355  }
356  break;
357  }
358  case 4: // stressField Neumann condition
359  {
360  DebugStop();//Nao implementado!!!
361  break;
362  }
363  default:
364  PZError << "TPZElastitity3D::ContributeBC error - Wrong boundary condition type" << std::endl;
365  }
366 }
TPZFunction::Execute
virtual void Execute(const TPZVec< REAL > &x, TPZVec< TVar > &f, TPZFMatrix< TVar > &df)
Performs function computation.
Definition: pzfunction.h:38
TPZMaterial::fNumLoadCases
int fNumLoadCases
Defines the number of load cases generated by this material.
Definition: TPZMaterial.h:76
TPZMaterialData::x
TPZManVector< REAL, 3 > x
value of the coordinate at the integration point
Definition: pzmaterialdata.h:71
TPZManVector
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
bc
clarg::argBool bc("-bc", "binary checkpoints", false)
TPZElast3Dnlinear::Contribute
virtual void Contribute(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ef) override
It computes a contribution to the residual vector at one integration point.
Definition: TPZElast3Dnlinear.cpp:31
TPZBndCond::Type
int Type()
Definition: pzbndcond.h:249
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
TPZMaterialData::fShapeType
MShapeFunctionType fShapeType
Definition: pzmaterialdata.h:43
TPZElast3Dnlinear.h
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
TPZElasticity3D::ClassId
virtual int ClassId() const override
Define the class id associated with the class.
Definition: pzelast3d.cpp:1272
TPZMaterialData::dphix
TPZFNMatrix< 660, REAL > dphix
values of the derivative of the shape functions
Definition: pzmaterialdata.h:59
pzbndcond.h
Contains the TPZBndCond class which implements a boundary condition for TPZMaterial objects...
TPZElast3Dnlinear::ContributeVecShapeAux
void ContributeVecShapeAux(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef)
Definition: TPZElast3Dnlinear.cpp:101
TPZElasticity3D::fPreStress
TPZManVector< REAL > fPreStress
Definition: pzelast3d.h:281
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
TPZBndCond
This class defines the boundary condition for TPZMaterial objects.
Definition: pzbndcond.h:29
TPZElast3Dnlinear::ContributeBC
virtual void ContributeBC(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc) override
Implements Dirichlet and Neumann boundary conditions.
Definition: TPZElast3Dnlinear.cpp:62
TPZMatrix::Rows
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
test.res
string res
Definition: test.py:151
TPZElast3Dnlinear::FillDataRequirements
virtual void FillDataRequirements(TPZMaterialData &data) override
Fill material data parameter with necessary requirements for the Contribute method.
Definition: TPZElast3Dnlinear.cpp:91
TPZBndCond::Val1
TPZFMatrix< STATE > & Val1()
Definition: pzbndcond.h:253
TPZElast3Dnlinear::ContributeVecShapeBCAux
void ContributeVecShapeBCAux(TPZMaterialData &data, REAL weight, TPZFMatrix< STATE > &ek, TPZFMatrix< STATE > &ef, TPZBndCond &bc)
Definition: TPZElast3Dnlinear.cpp:232
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
TPZElast3Dnlinear::ClassId
int ClassId() const override
Define the class id associated with the class.
Definition: TPZElast3Dnlinear.cpp:87
TPZMaterial::fForcingFunction
TPZAutoPointer< TPZFunction< STATE > > fForcingFunction
Pointer to forcing function, it is the right member at differential equation.
Definition: TPZMaterial.h:47
TPZElasticity3D::fPoisson
STATE fPoisson
Poisson&#39;s ratio.
Definition: pzelast3d.h:261
TPZElasticity3D::fForce
TPZManVector< STATE, 3 > fForce
External forces.
Definition: pzelast3d.h:267
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
TPZElasticity3D::fE
STATE fE
Young&#39;s modulus.
Definition: pzelast3d.h:258
TPZMaterialData::sol
TPZSolVec sol
vector of the solutions at the integration point
Definition: pzmaterialdata.h:77
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
Generated on Sun Aug 16 2020 11:17:04 for NeoPZ by   doxygen 1.8.13