NeoPZ
TPZPlacaOrthotropic.cpp
Go to the documentation of this file.
1 
5 #include "pzintel.h"
6 #include "TPZPlacaOrthotropic.h"
7 #include "pzgnode.h"
8 #include "pzgeoel.h"
9 #include "pzquad.h"
10 //#include "pztempmat.h"
11 #include "TPZMaterial.h"
12 
13 using namespace std;
14 
16 
17  fGeoEl = gel;
18  fH = zmax-zmin;
19  fZMin = zmin;
20  fZMax = zmax;
21  fTensorVar = -1;
22  fIntel = dynamic_cast<TPZInterpolatedElement *>(gel->Reference());
23  if(fIntel) {
24  fTensorVar = fIntel->Material()->VariableIndex("Tensor");
25  }
26 }
27 
29 
30  fGeoEl = 0;
31  fH = 0.;
32  fZMin = 0.;
33  fZMax = 0.;
34  fIntel = 0;
35  fTensorVar = -1;
36 }
37 
39  if(fTensorVar == -1) {
40  if(!fIntel) fIntel = dynamic_cast<TPZInterpolatedElement *>(fGeoEl->Reference());
41  if(!fIntel) return;
42  fTensorVar = fIntel->Material()->VariableIndex("Tensor");
43  }
44  if(fTensorVar == -1) return;
45  // REAL ksi = -1+2.*(z-fZMin)/(fZMax-fZMin);
46  // if(ksi < -1. || ksi > 1.) return;
47  // TPZManVector<REAL,3> co(3,0.);
48  // co[2] = ksi;
49  TPZManVector<STATE> tensor(9);
50  fIntel->Solution(ksi,fTensorVar,tensor);
51  int i;
52  for(i=0; i<9; i++) {//original
53  T(i%3,i/3) = tensor[i];
54  }
55  // for(i=0; i<3; i++) {//�sim�rico
56  // for(j=0; j<3; j++) {
57  // T(i,j) = tensor[3*i+j];
58  // }
59  // }
60 }
61 
62 REAL TPZPlacaOrthotropic::Moment(REAL zref, TPZVec<REAL> &normal, TPZVec<REAL> &direction){
63  TPZInt1d rule(8);
64  int np = rule.NPoints();
65  TPZManVector<REAL,3> pos(1,0.), ksi(3,0.);
66  TPZFNMatrix<9> tensor(3,3,0.);
67  int ip;
68  REAL moment = 0.;
69  for(ip = 0; ip<np; ip++) {
70  REAL weight;
71  rule.Point(ip,pos,weight);
72  ksi[2] = pos[0];
73  Tensor(ksi,tensor);
74  REAL z = fZMin + (fZMax-fZMin)*(pos[0]+1.)/2.;//original
75  //REAL f = (fZMin + fZMax)/2.0;//n�
76  //REAL z = fH*pos[0]/2. + f;//n�
77  REAL tension = 0.;
78  int n,d;
79  for(n=0; n<3; n++) {
80  for(d=0; d<3; d++) {
81  tension += normal[n]*tensor(n,d)*direction[d];
82  }
83  }
84  moment += weight*tension*(z-zref);//original
85  //moment += weight*tension*z;//n�
86  }
87  moment *= fH/2.;
88  return moment;
89 
90 }
91 
93  TPZInt1d rule(8);
94  int np = rule.NPoints();
95  TPZManVector<REAL,3> pos(1,0.), ksi(3,0.);
96  TPZFNMatrix<9> tensor(3,3,0.);
97  int ip;
98  REAL force = 0.;
99  for(ip = 0; ip<np; ip++) {
100  REAL weight;
101  rule.Point(ip,pos,weight);
102  ksi[2] = pos[0];
103  Tensor(ksi,tensor);
104  REAL tension = 0.;
105  int n,d;
106  for(n=0; n<3; n++) {
107  for(d=0; d<3; d++) {
108  tension += normal[n]*tensor(n,d)*direction[d];
109  }
110  }
111  force += weight*tension;
112  }
113  force *= fH/2.;
114  return force;
115 }
116 
118 
119  fIntel->Print();
120  cout << "Plate height : " ;
121  cout << this->Height();//this �o objeto placa
122 
123 }
124 
125 //se fosse TPZMulticamadaOrthotropic::Print, o this seria multcam
126 //se a func� fosse declarada como void Print(){...}, n� haveria This.
127 //caso chamemos a func� Print com o objeto placa2, ao inv� de placa->Print, o This �o objeto placa2.
128 
129 
131  fIntel = dynamic_cast<TPZInterpolatedElement *>(fGeoEl->Reference());
132  if(fIntel) {
133  fTensorVar = fIntel->Material()->VariableIndex("Tensor");
134  }
135 
136 }
137 
138 REAL TensionNorm(TPZFMatrix<REAL> &tension,int dimrow,int dimcol);
139 void TPZPlacaOrthotropic::PrintTensors(std::ostream &out,TPZFMatrix<REAL> &tensorin,TPZFMatrix<REAL> &tensorout) {
140 
141  TPZInt1d rule(8);
142  int np = rule.NPoints();
143  TPZFMatrix<REAL> TensorOut(np,9),DiffTension(np,9);
144  TPZManVector<REAL,3> pos(1,0.), ksi(3,0.), x(3,0.);
145  TPZFNMatrix<9> tensor(3,3,0.);
146  int ip,i,j;
147  for(ip = 0; ip<np; ip++) {
148  REAL weight;
149  rule.Point(ip,pos,weight);
150  ksi[2] = pos[0];
151  fGeoEl->X(ksi,x);
152  Tensor(ksi,tensor);
153  out << "Tensor at pos " << pos[0] << " x " << x << endl;
154  for(i=0; i<3; i++) for(j=0; j<3; j++) tensorout(ip,3*i+j) = tensor(i,j);
155  for(j=0; j<9; j++) DiffTension(ip,j) = tensorout(ip,j) - tensorin(ip,j);
156  }
157  out << " tensor(n) - tensor(n-1) = " << endl;
158  for(i=0; i<np; i++){
159  for(j=0; j<9; j++){
160  out << DiffTension(i,j) << " ";
161  }
162  out << endl;
163  }
164  REAL normat = TensionNorm(DiffTension,np,9);
165  out << "Norma of tensor out - tensor in = " << normat << endl;
166 }
167 
168 void TPZPlacaOrthotropic::PrintTensors(std::ostream &out) {
169  TPZInt1d rule(8);
170  int np = rule.NPoints();
171  TPZManVector<REAL,3> pos(1,0.), ksi(3,0.), x(3,0.);
172  TPZFNMatrix<9> tensor(3,3,0.);
173  int ip;
174  REAL weight;
175  TPZStack<REAL> point;
176  point.Push(-1.0);
177  for(ip = 0; ip<np; ip++){
178  rule.Point(ip,pos,weight);
179  point.Push(pos[0]);
180  }
181  point.Push(1.0);
182  np += 2;
183  for(ip = 0; ip<np; ip++) {
184 
185  ksi[2] = point[ip];
186  fGeoEl->X(ksi,x);
187  Tensor(ksi,tensor);
188  out << "Tensor at pos " << point[ip] << " x " << x << " -> ";
189  int i,j;
190  for(i=0; i<3; i++) {
191  for(j=0; j<3; j++){
192  REAL tension = tensor(i,j);
193  if(fabs(tension) < 1.e-10) out << 0 << " ";
194  else out << tension << " ";
195  }
196  //out << endl;
197  }
198  out << endl;
199  REAL normat = TensionNorm(tensor,3,3);
200  // out << "Norma do tensor " << normat << endl;
201  cout << "Norma do tensor " << normat << endl;
202  }
203 }
204 
206 REAL TensionNorm(TPZFMatrix<REAL> &tension,int dimrow,int dimcol) {
207 
208  int i,j;
209  REAL val = 0.0;
210  for(i=0;i<dimrow;i++) for(j=0;j<dimcol;j++) val += tension(i,j)*tension(i,j);
211  return sqrt(val);
212 }
213 
214 REAL TPZPlacaOrthotropic::GradMoment(REAL zref, TPZVec<REAL> &graddir, TPZVec<REAL> &normal, TPZVec<REAL> &direction){
215  TPZInt1d rule(8);
216  int np = rule.NPoints();
217  TPZManVector<REAL,3> pos(1,0.), ksi1(3,0.), ksi2(3,0.), x1(3,0.), x2(3,0.);
218  TPZFNMatrix<9> tensor1(3,3,0.), tensor2(3,3,0.), tensordif(3,3,0.);
219  ksi1[0] = -graddir[0];
220  ksi1[1] = -graddir[1];
221  ksi2[0] = graddir[0];
222  ksi2[1] = graddir[1];
223  fGeoEl->X(ksi1,x1);
224  fGeoEl->X(ksi2,x2);
225  REAL dist = sqrt((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]));
226  int ip;
227  REAL moment = 0.;
228  for(ip = 0; ip<np; ip++) {
229  REAL weight;
230  rule.Point(ip,pos,weight);
231  ksi1[2] = pos[0];
232  ksi2[2] = pos[0];
233  Tensor(ksi1,tensor1);
234  Tensor(ksi2,tensor2);
235  tensordif = tensor2;
236  tensordif -= tensor1;
237  REAL z = fZMin + (fZMax-fZMin)*(pos[0]+1.)/2.;// = (z+f)
238  REAL tension = 0.;
239  int n,d;
240  for(n=0; n<3; n++) {
241  for(d=0; d<3; d++) {
242  tension += normal[n]*tensordif(n,d)*direction[d];
243  }
244  }
245  moment += weight*tension*(z-zref);//original
246  //moment += weight*tension*z;//n�
247  }
248  moment *= fH/(2.*dist);
249  return moment;
250 
251 }
252 
254  TPZInt1d rule(8);
255  int np = rule.NPoints();
256  TPZManVector<REAL,3> pos(1,0.), ksi1(3,0.), ksi2(3,0.), x1(3,0.), x2(3,0.);
257  TPZFNMatrix<9> tensor1(3,3,0.), tensor2(3,3,0.), tensordif(3,3,0.);
258  ksi1[0] = -graddir[0];
259  ksi1[1] = -graddir[1];
260  ksi2[0] = graddir[0];
261  ksi2[1] = graddir[1];
262  fGeoEl->X(ksi1,x1);
263  fGeoEl->X(ksi2,x2);
264  REAL dist = sqrt((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]));
265  int ip;
266  REAL force = 0.;
267  for(ip = 0; ip<np; ip++) {
268  REAL weight;
269  rule.Point(ip,pos,weight);
270  ksi1[2] = pos[0];
271  ksi2[2] = pos[0];
272  Tensor(ksi1,tensor1);
273  Tensor(ksi2,tensor2);
274  tensordif = tensor2;
275  tensordif -= tensor1;
276  REAL tension = 0.;
277  int n,d;
278  for(n=0; n<3; n++) {
279  for(d=0; d<3; d++) {
280  tension += normal[n]*tensordif(n,d)*direction[d];
281  }
282  }
283  force += weight*tension;
284  }
285  force *= fH/(2.*dist);
286  return force;
287 }
288 
290  if(fTensorVar == -1) {
291  if(!fIntel) fIntel = dynamic_cast<TPZInterpolatedElement *>(fGeoEl->Reference());
292  if(!fIntel) return;
293  fTensorVar = fIntel->Material()->VariableIndex("Tensor");
294  }
295  if(fTensorVar == -1) return;
296  // REAL ksi = -1+2.*(z-fZMin)/(fZMax-fZMin);
297  // if(ksi < -1. || ksi > 1.) return;
298  // TPZManVector<REAL,3> co(3,0.);
299  // co[2] = ksi;
300  TPZManVector<REAL,3> pos(1,0.), ksi1(3,0.), ksi2(3,0.), x1(3,0.), x2(3,0.);
301  TPZManVector<STATE> tensor1(0,0.), tensor2(9,0.);
302  ksi1[0] = ksi[0]-graddir[0];
303  ksi1[1] = ksi[1]-graddir[1];
304  ksi1[2] = ksi[2]-graddir[2];
305  ksi2[0] = ksi[0]+graddir[0];
306  ksi2[1] = ksi[1]+graddir[1];
307  ksi2[2] = ksi[2]+graddir[2];
308  fGeoEl->X(ksi1,x1);
309  fGeoEl->X(ksi2,x2);
310  REAL dist = sqrt((x1[0]-x2[0])*(x1[0]-x2[0])+(x1[1]-x2[1])*(x1[1]-x2[1]));
311  TPZManVector<REAL> tensor(9);
312  fIntel->Solution(ksi1,fTensorVar,tensor1);
313  fIntel->Solution(ksi2,fTensorVar,tensor2);
314  int i;
315  for(i=0; i<9; i++) {//�sim�rico OK
316  T(i%3,i/3) = (tensor2[i]-tensor1[i])/dist;
317  }
318  // for(i=0; i<3; i++) {//n� sim�rico
319  // for(j=0; j<3; j++) {
320  // T(i,j) = (tensor2[3*i+j]-tensor1[3*i+j])/dist;
321  // }
322  // }
323 }
virtual void Point(int ip, TPZVec< REAL > &pos, REAL &w) const
Returns i-th point at master element and related weight.
Definition: pzquad.cpp:101
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
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
Contains declaration of TPZGeoNode class which defines a geometrical node.
TPZPlacaOrthotropic()
Default constructor.
virtual int VariableIndex(const std::string &name)
Returns the variable index associated with the name.
REAL Force(TPZVec< REAL > &normal, TPZVec< REAL > &direction)
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
void PrintTensors(std::ostream &out)
Handles the numerical integration for one-dimensional problems. Numerical Integration.
Definition: pzquad.h:36
virtual TPZMaterial * Material() const
Identify the material object associated with the element.
Definition: pzcompel.cpp:959
void Tensor(TPZVec< REAL > &ksi, TPZFMatrix< REAL > &T)
Returns the tensions tensor of the shell.
void GradTensor(TPZVec< REAL > &graddir, TPZVec< REAL > &ksi, TPZFMatrix< REAL > &gradtensor)
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
Defines the behaviour of all geometric elements. GeometryTPZGeoEl is the common denominator for all g...
Definition: pzgeoel.h:43
REAL GradMoment(REAL zref, TPZVec< REAL > &graddir, TPZVec< REAL > &normal, TPZVec< REAL > &direction)
REAL GradForce(TPZVec< REAL > &graddir, TPZVec< REAL > &normal, TPZVec< REAL > &direction)
TPZCompEl * Reference() const
Return a pointer to the element referenced by the geometric element.
Definition: pzgeoel.h:750
virtual int NPoints() const
Returns number of points for the cubature rule related.
Definition: pzquad.cpp:95
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
Contains the TPZPlacaOrthotropic class.
REAL TensionNorm(TPZFMatrix< REAL > &tension, int dimrow, int dimcol)
Returns norm of the tension.
REAL dist(TPZVec< T1 > &vec1, TPZVec< T1 > &vec2)
Definition: pzvec_extras.h:124
REAL Height(TPZGeoMesh *gmesh)
Definition: substruct.cpp:1245
Contains declaration of TPZInterpolatedElement class which implements computational element of the in...
REAL Moment(REAL zref, TPZVec< REAL > &normal, TPZVec< REAL > &direction)
Implements computational element based on an interpolation space. Computational Element.
Definition: pzintel.h:27