NeoPZ
pzshapepiram.cpp
Go to the documentation of this file.
1 
6 #include "pzshapepiram.h"
7 #include "pzshapequad.h"
8 #include "pzshapetriang.h"
9 #include "pzshapelinear.h"
10 #include "pzshapepoint.h"
11 #include "pzmanvector.h"
12 #include "pzerror.h"
13 #include "pzreal.h"
14 
15 using namespace std;
16 
18 namespace pzshape {
19 
20  /*Projection of the point within a piramide to a rib*/
21  REAL TPZShapePiram::gRibTrans3dPiram1d[8][3] = {//par�metros de arestas
22  { 1., 0.,0.} , { 0., 1.,0.} ,//percorre o sentido
23  {-1., 0.,0.} , { 0.,-1.,0.} ,//da aresta segundo
24  { .5, .5,1.} , {-.5,.5 ,1.} ,//SideNodes[8][2]
25  {-.5,-.5,1.} , { .5,-.5,1.}
26  };
27  //n�o tem vetor associado -> OK!
28 
29  //Projection of the point within a piramide to a face
30  REAL TPZShapePiram::gFaceTrans3dPiram2d[5][2][3] = {//par�metros de faces
31  { { 1., 0., 0.},{ 0., 1.,0.} },//0 1 2 3 : percorre os eixos segundo
32  { { 1.,-.5,-.5},{ 0., 1.,1.} },//0 1 4 : FaceNodes[5][4]
33  { { .5, 1.,-.5},{-1., 0.,1.} },//1 2 4
34  { { 1., .5,-.5},{ 0.,-1.,1.} },//3 2 4 ; original-> 2 3 4 : {-1., .5,-.5},{ 0.,-1.,1.}
35  { {-.5, 1.,-.5},{1.,0.,1.} } //0 3 4 ; original-> 3 0 4 : {-.5,-1.,-.5},{ 1., 0.,1.}
36  };
37 
38  REAL TPZShapePiram::gFaceSum3dPiram2d[5][2] = { {.0,0.},{-.5,0.},{-.5,0.},{-.5,0.},{-.5,0.} };//{ {.0,0.},{-.5,0.},{-.5,0.},{-.5,0.},{-.5,0.} };//original
39 
40 
47  void TPZShapePiram::ShapeGenerating(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi)
48  {
49  int is;
50  // contribute the ribs
51  for(is=NCornerNodes; is<NCornerNodes+8; is++)
52  {
53  int is1,is2;
54  is1 = ContainedSideLocId(is,0);
55  is2 = ContainedSideLocId(is,1);
56  phi(is,0) = phi(is1,0)*phi(is2,0);
57  dphi(0,is) = dphi(0,is1)*phi(is2,0)+phi(is1,0)*dphi(0,is2);
58  dphi(1,is) = dphi(1,is1)*phi(is2,0)+phi(is1,0)*dphi(1,is2);
59  dphi(2,is) = dphi(2,is1)*phi(is2,0)+phi(is1,0)*dphi(2,is2);
60  }
61  // contribution of the faces
62  // quadrilateral face
63  is = 13;
64  {
65  int is1,is2;
66  is1 = ShapeFaceId[0][0];// SideConnectLocId(is,0);
67  is2 = ShapeFaceId[0][2];// SideConnectLocId(is,2);
68  phi(is,0) = phi(is1,0)*phi(is2,0);
69  dphi(0,is) = dphi(0,is1)*phi(is2,0)+phi(is1,0)*dphi(0,is2);
70  dphi(1,is) = dphi(1,is1)*phi(is2,0)+phi(is1,0)*dphi(1,is2);
71  dphi(2,is) = dphi(2,is1)*phi(is2,0)+phi(is1,0)*dphi(2,is2);
72  }
73  is++;
74  // triangular faces
75  for(;is<18; is++)
76  {
77  int is1,is2,is3;
78  is1 = ContainedSideLocId(is,0);
79  is2 = ContainedSideLocId(is,1);
80  is3 = ContainedSideLocId(is,2);
81  phi(is,0) = phi(is1,0)*phi(is2,0)*phi(is3,0);
82  dphi(0,is) = dphi(0,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(0,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(0,is3);
83  dphi(1,is) = dphi(1,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(1,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(1,is3);
84  dphi(2,is) = dphi(2,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(2,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(2,is3);
85 
86  }
87  {
88  int is1 = 0;
89  int is2 = 2;
90  int is3 = 4;
91  phi(is,0) = phi(is1,0)*phi(is2,0)*phi(is3,0);
92  dphi(0,is) = dphi(0,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(0,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(0,is3);
93  dphi(1,is) = dphi(1,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(1,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(1,is3);
94  dphi(2,is) = dphi(2,is1)*phi(is2,0)*phi(is3,0)+phi(is1,0)*dphi(2,is2)*phi(is3,0)+phi(is1,0)*phi(is2,0)*dphi(2,is3);
95  }
96 
97  // Make the generating shape functions linear and unitary
98  // contribute the ribs
99  for(is=NCornerNodes; is<NCornerNodes+4; is++)
100  {
101  int isface = 13;
102  phi(is,0) += phi(isface,0);
103  dphi(0,is) += dphi(0,isface);
104  dphi(1,is) += dphi(1,isface);
105  dphi(2,is) += dphi(2,isface);
106  }
107  // scaling the shapefunctions
108  {
109  REAL sidescale[] = {1.,1.,1.,1.,1.,4.,4.,4.,4.,4.,4.,4.,4.,16.,27.,27.,27.,27.,64.};
110  for(is=5; is<NSides; is++)
111  {
112  phi(is,0) *= sidescale[is];
113  dphi(0,is) *= sidescale[is];
114  dphi(1,is) *= sidescale[is];
115  dphi(2,is) *= sidescale[is];
116  }
117  }
118 
119 
120  }
121 
123 
124  CornerShape(pt,phi,dphi);
125  bool linear = true;
126  int is,d;
127  for(is=NCornerNodes; is<NSides; is++){
128  if(order[is-NCornerNodes] > 1){
129  linear = false;
130  }
131  }
132  if(linear) return;
133 
134  TPZFNMatrix<100> phiblend(NSides,1),dphiblend(Dimension,NSides);
135  for(is=0; is<NCornerNodes; is++)
136  {
137  phiblend(is,0) = phi(is,0);
138  for(d=0; d<Dimension; d++)
139  {
140  dphiblend(d,is) = dphi(d,is);
141  }
142  }
143  ShapeGenerating(pt,phiblend,dphiblend);
144  // if(order[13]<2) return;//order tem as ordens dos lados do elemento
145  int shape = 5;
146  //rib shapes
147  for (int rib = 0; rib < 8; rib++) {//todas as arestas
148  if (order[rib] <2 ) continue;
149  REAL outval;
150  ProjectPoint3dPiramToRib(rib,pt,outval);
151  TPZVec<int64_t> ids(2);
152  TPZManVector<REAL,1> outvalvec(1,outval);
153  int id0,id1;
154  id0 = SideNodes[rib][0];
155  id1 = SideNodes[rib][1];
156  ids[0] = id[id0];
157  ids[1] = id[id1];
158  REAL store1[20],store2[60];
159  int nshape = order[rib]-1;//three orders : order in x , order in y and order in z
160  TPZFMatrix<REAL> phin(nshape,1,store1,20),dphin(3,nshape,store2,60);
161  phin.Zero();
162  dphin.Zero();
163  TPZShapeLinear::ShapeInternal(outvalvec,order[rib],phin,dphin,TPZShapeLinear::GetTransformId1d(ids));//ordin = ordem de um lado
164  TransformDerivativeFromRibToPiram(rib,nshape,dphin);
165  for (int i = 0; i < nshape; i++) {
166  phi(shape,0) = phiblend(rib+5,0)*phin(i,0);
167  for(int xj=0;xj<3;xj++) {
168  dphi(xj,shape) = dphiblend(xj ,rib+5) * phin( i, 0) +
169  phiblend(rib+5, 0 ) * dphin(xj,i);
170  }
171  shape++;
172  }
173  }
174 
175  for (int face = 0; face < 5; face++) {
176  if (order[face+8] < 2) continue;
177  if(face>0 && order[face+8]<=2) continue;//s� a face 13 tem shape associada com ordem p=2
178  TPZManVector<REAL,2> outval(2);
179  ProjectPoint3dPiramToFace(face,pt,outval);
180  int ord1;//,ord2;
181  ord1 = order[8+face];
182  //ord2 = ord1;
183  //FaceOrder(face,ord1,ord2);//ordem da face
184  if(face && ord1<3) continue;//uma face com ordem < 3 n�o tem shape associada
185  int nshape;
186  if(!face) nshape = (ord1-1)*(ord1-1);//face quadrada
187  else nshape = (ord1-2)*(ord1-1)/2;//face triangular
188  TPZFNMatrix<60> phin(nshape,1),dphin(3,nshape);//ponto na face
189  phin.Zero();
190  dphin.Zero();
191  TPZManVector<int64_t> ids(4);
192  // int id0,id1,id2;
193  int i;
194  if(!face) for(i=0;i<4;i++) ids[i] = id[FaceNodes[face][i]];
195  else for(i=0;i<3;i++) ids[i] = id[FaceNodes[face][i]];
196  // id0 = ShapeFaceId[face][0];//indice das shapes da face que compoem a shape atual
197  // id1 = ShapeFaceId[face][1];//equivale a FaceIdsCube(face,ids,id,id0,id1);
198  // id2 = ShapeFaceId[face][2];//if(face == 0) id3 = ShapeFaceId[face][3];
199  int transid;
200  if(!face) {
201  transid = TPZShapeQuad::GetTransformId2dQ(ids);
202  TPZShapeQuad::ShapeInternal(outval,ord1-2,phin,dphin,transid);//ordin = ordem de um lado
203  } else {
204  ids.Resize(3);
205  transid = TPZShapeTriang::GetTransformId2dT(ids);
206  outval[0] = (outval[0]+1.)/2.;//devido a corre��o na fun��o
207  outval[1] = (outval[1]+1.)/2.;//Shape2dTriangleInternal(..) : correto aqui
208  TPZShapeTriang::ShapeInternal(outval,ord1-2,phin,dphin,transid);//ordin = ordem de um lado
209  int c = dphin.Cols();//isto da (p-2)(p-1)/2 ; ord1 = p ; correto aqui
210  for(int i=0;i<c;i++) {
211  dphin(0,i) /= 2.;//correcao da derivada OK! aqui
212  dphin(1,i) /= 2.;
213  dphin(2,i) /= 2.;
214  }
215  }
216  TransformDerivativeFromFaceToPiram(face,nshape,dphin);//ordin = numero de shapes
217  for(i=0;i<nshape;i++) {
218  phi(shape,0) = phiblend(face+13,0)*phin(i,0);//face quadril�teral
219  for(int xj=0;xj<3;xj++) {
220  dphi(xj,shape) = dphiblend(xj,face+13)* phin(i ,0) +
221  phiblend(face+13, 0)*dphin(xj,i);
222  }
223  shape++;
224  }
225  }
226  if(order[13]<3) return;//n�o h� ordens para cantos
227  //volume shapes
228  int nshape=0,i;
229  for(i=0;i<order[13]-2;i++) {
230  nshape += (i+1)*(i+2) / 2;
231  }
232  TPZFNMatrix<40> phin(nshape,1),dphin(3,nshape);
233  phin.Zero();
234  dphin.Zero();
235  ShapeInternal(pt,order[13],phin,dphin);
236  for(i=0;i<nshape;i++) {
237  phi(shape,0) = phiblend(NSides-1,0)*phin(i,0);
238  for(int xj=0;xj<3;xj++) {
239  dphi(xj,shape) = dphiblend(xj,NSides-1) * phin(i ,0) +
240  phiblend(NSides-1, 0) * dphin(xj,i);
241  }
242  shape++;
243  }
244  }
245 
246  void TPZShapePiram::ShapeCorner(const TPZVec<REAL> &pt,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
247  if(fabs(pt[0])<1.e-10 && fabs(pt[1])<1.e-10 && pt[2]==1.) {
248  //para testes com transformaçoes geometricas-->>Que é o que faz o RefPattern!!
249  //(0,0,1) nunca é um ponto de integração
250  phi(0,0) = 0.;
251  phi(1,0) = 0.;
252  phi(2,0) = 0.;
253  phi(3,0) = 0.;
254  phi(4,0) = 1.;
255  dphi(0,0) = -0.25;
256  dphi(1,0) = -0.25;
257  dphi(2,0) = -0.25;
258  dphi(0,1) = 0.25;
259  dphi(1,1) = -0.25;
260  dphi(2,1) = -0.25;
261  dphi(0,2) = 0.25;
262  dphi(1,2) = 0.25;
263  dphi(2,2) = -0.25;
264  dphi(0,3) = -0.25;
265  dphi(1,3) = 0.25;
266  dphi(2,3) = -0.25;
267  dphi(0,4) = 0;
268  dphi(1,4) = 0;
269  dphi(2,4) = 1.;
270 
271 
272 
273  return;
274  }
275 
276  REAL T0xz = .5*(1.-pt[2]-pt[0]) / (1.-pt[2]);
277  REAL T0yz = .5*(1.-pt[2]-pt[1]) / (1.-pt[2]);
278  REAL T1xz = .5*(1.-pt[2]+pt[0]) / (1.-pt[2]);
279  REAL T1yz = .5*(1.-pt[2]+pt[1]) / (1.-pt[2]);
280  REAL lmez = (1.-pt[2]);
281  phi(0,0) = T0xz*T0yz*lmez;
282  phi(1,0) = T1xz*T0yz*lmez;
283  phi(2,0) = T1xz*T1yz*lmez;
284  phi(3,0) = T0xz*T1yz*lmez;
285  phi(4,0) = pt[2];
286  REAL lmexmez = 1.-pt[0]-pt[2];
287  REAL lmeymez = 1.-pt[1]-pt[2];
288  REAL lmaxmez = 1.+pt[0]-pt[2];
289  REAL lmaymez = 1.+pt[1]-pt[2];
290  dphi(0,0) = -.25*lmeymez / lmez;
291  dphi(1,0) = -.25*lmexmez / lmez;
292  dphi(2,0) = -.25*(lmeymez+lmexmez-lmexmez*lmeymez/lmez) / lmez;
293 
294  dphi(0,1) = .25*lmeymez / lmez;
295  dphi(1,1) = -.25*lmaxmez / lmez;
296  dphi(2,1) = -.25*(lmeymez+lmaxmez-lmaxmez*lmeymez/lmez) / lmez;
297 
298  dphi(0,2) = .25*lmaymez / lmez;
299  dphi(1,2) = .25*lmaxmez / lmez;
300  dphi(2,2) = -.25*(lmaymez+lmaxmez-lmaxmez*lmaymez/lmez) / lmez;
301 
302  dphi(0,3) = -.25*lmaymez / lmez;
303  dphi(1,3) = .25*lmexmez / lmez;
304  dphi(2,3) = -.25*(lmaymez+lmexmez-lmexmez*lmaymez/lmez) / lmez;
305 
306  dphi(0,4) = 0.0;
307  dphi(1,4) = 0.0;
308  dphi(2,4) = 1.0;
309  }
310  void TPZShapePiram::SideShape(int side, TPZVec<REAL> &point, TPZVec<int64_t> &id, TPZVec<int> &order, TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
311 
312  if(side<0 || side>18) PZError << "TPZCompElPi3d::SideShapeFunction. Bad paramenter side.\n";
313  else if(side==18) Shape(point,id,order,phi,dphi);
314  else if(side<5) TPZShapePoint::Shape(point,id,order,phi,dphi);
315  else if(side<13) {//5 a 12
316  TPZShapeLinear::Shape(point,id,order,phi,dphi);
317  } else if(side == 13) {
318  TPZShapeQuad::Shape(point,id,order,phi,dphi);
319  } else if(side<18) {//faces 13 a 17
320  TPZShapeTriang::Shape(point,id,order,phi,dphi);
321  }
322 
323  }
324 
325  void TPZShapePiram::ShapeOrder(TPZVec<int64_t> &id, TPZVec<int> &order, TPZGenMatrix<int> &shapeorders)//, TPZVec<int64_t> &sides
326  {
327  //DebugStop();
328  int64_t nsides = TPZShapePiram::NSides;
329  int nshape;
330 
331  int linha = 0;
332  for (int side = 0; side < nsides; side++)
333  {
334 
335  nshape = 1;
336  if(side >= NCornerNodes) nshape = NConnectShapeF(side,order[side-NCornerNodes]);
337  int sideorder = 1;
338  if(side >= NCornerNodes) sideorder = order[side-NCornerNodes];
339 
340  TPZGenMatrix<int> locshapeorders(nshape,3);
341  SideShapeOrder(side, id, sideorder, locshapeorders);
342 
343  int nlin = locshapeorders.Rows();
344  int ncol = locshapeorders.Cols();
345 
346  for (int il = 0; il<nlin; il++)
347  {
348  for (int jc = 0; jc<ncol; jc++)
349  {
350  shapeorders(linha, jc) = locshapeorders(il, jc);
351  }
352  linha++;
353  }
354  }
355  }
356 
357 
358  void TPZShapePiram::SideShapeOrder(int side, TPZVec<int64_t> &id, int order, TPZGenMatrix<int> &shapeorders)
359  {
360  //DebugStop();
361  if (side<=4)
362  {
363  if (shapeorders.Rows() != 1)
364  {
365  DebugStop();
366  }
367  shapeorders(0,0) = 1;
368  shapeorders(0,1) = 0;
369  shapeorders(0,2) = 0;
370  }
371  else if (side>=5 && side<=12)
372  {
373  int nshape = order-1;
374  if (shapeorders.Rows() != nshape)
375  {
376  DebugStop();
377  }
378  for (int ioy = 0; ioy < order-1; ioy++)
379  {
380  shapeorders(ioy,0) = ioy+2;
381  }
382  }
383  else if (side == 13)
384  {
385  // quadrilatero
386  if (shapeorders.Rows() != (order-1)*(order-1))
387  {
388  DebugStop();
389  }
390  TPZStack<int> lowersides;
391  LowerDimensionSides(side, lowersides);
392  lowersides.Push(side);
393 
394  int nnodes = NSideNodes(side);
395 
396  TPZManVector<int64_t, 4> locid(nnodes);
397  for (int node=0; node<locid.size(); node++) {
398  locid[node] = id[ContainedSideLocId(side, node)];// SideNodeLocId( side, node);
399  }
400 
401  int nshape = (order-1)*(order-1);
402  TPZGenMatrix<int> locshapeorders(nshape,3);
403 
404 
405  TPZShapeQuad::SideShapeOrder(8,locid, order, locshapeorders);
406 
407  // temos que arrumar a saida de locshapeorders para adequar a orientacao dos vetores que geram
408  // a face do lado side
409 
410  // aqui o locshapeorders esta armazenado so para x e y
411  for (int il = 0; il<nshape; il++)
412  {
413  shapeorders(il, 0) = locshapeorders(il, 0);
414  shapeorders(il, 1) = locshapeorders(il, 1);
415  shapeorders(il, 2) = locshapeorders(il, 2);
416  }
417  }
418  else if (side >= 14 && side <= 17)
419  {
420  //triangulos
421  int nshape = (order-2)*(order-1)/2;
422  if (shapeorders.Rows() != nshape)
423  {
424  DebugStop();
425  }
426  TPZStack<int> lowersides;
427  LowerDimensionSides(side, lowersides);
428  lowersides.Push(side);
429 
430  int nnodes = NSideNodes(side);
431 
432  TPZManVector<int64_t, 4> locid(nnodes);
433  for (int node=0; node<locid.size(); node++) {
434  locid[node] = id[ContainedSideLocId(side, node)];
435  }
436 
437  TPZGenMatrix<int> locshapeorders(nshape,3);
438 
439 
440  TPZShapeTriang::SideShapeOrder(6,locid, order, locshapeorders);
441 
442  for (int il = 0; il<nshape; il++)
443  {
444  shapeorders(il, 0) = locshapeorders(il, 0);
445  shapeorders(il, 1) = locshapeorders(il, 1);
446  shapeorders(il, 2) = locshapeorders(il, 2);
447  }
448  }
449  else
450  {
451  // interno
452  int totsum = 0;
453  for(int i=0;i<order - 2;i++) {
454  totsum += (i+1)*(i+2) / 2;
455  }
456  int nshape = totsum;
457 
458  if (shapeorders.Rows() != nshape) {
459  DebugStop();
460  }
461  int count = 0;
462  int ord = order-2;
463  for (int i=0;i<ord;i++) {
464  for (int j=0;j<ord;j++) {
465  for (int k=0;k<ord;k++) {
466  int soma = i+j+k;
467  if( i+j+k < ord ) // Duvida
468  {
469  shapeorders(count,0) = 3 + soma;
470  shapeorders(count,1) = 3 + soma;
471  shapeorders(count,2) = 3 + soma;
472  count++;
473  }
474  }
475  }
476  }
477  }
478 
479  }
480 
481  void TPZShapePiram::ShapeInternal(TPZVec<REAL> &x, int order,TPZFMatrix<REAL> &phi,
482  TPZFMatrix<REAL> &dphi) {
483  // calculate the values of the function and derivatives of the product of the orthogonal functions
484  if(order < 3) return;
485  int ord = order-2;
486  TPZFNMatrix<20, REAL> phi0(ord,1),phi1(ord,1),phi2(ord,1),
487  dphi0(1,ord),dphi1(1,ord),dphi2(1,ord);
488  TPZShapeLinear::fOrthogonal(x[0],ord,phi0,dphi0);//f e df -1<=x0<=1
489  TPZShapeLinear::fOrthogonal(x[1],ord,phi1,dphi1);//g e dg -1<=x1<=1
490  TPZShapeLinear::fOrthogonal(2.*x[2]-1.,ord,phi2,dphi2);//h e dh 0<=x3<=1 -> -1<=2*x2-1<=1
491  int index = 0; // integration point internal to pyramid
492  for (int i=0;i<ord;i++) {
493  for (int j=0;j<ord;j++) {
494  for (int k=0;k<ord;k++) {
495  if( i+j+k < ord ) {
496  //int index = ord*(ord*i+j)+k; //i,j,k � o grau das fun��es ortogonais
497  phi(index,0) = phi0(i,0)* phi1(j,0)* phi2(k,0);
498  dphi(0,index) = dphi0(0,i)* phi1(j,0)* phi2(k,0);
499  dphi(1,index) = phi0(i,0)*dphi1(0,j)* phi2(k,0);
500  dphi(2,index) = 2.*phi0(i,0)* phi1(j,0)*dphi2(0,k);
501  index++;
502  }
503  }
504  }
505  }
506  }
507 
508  void TPZShapePiram::ShapeInternal(int side, TPZVec<REAL> &x, int order,
509  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi){
510  if (side < 5 || side > 18) {
511  DebugStop();
512  }
513 
514  switch (side) {
515 
516  case 5:
517  case 6:
518  case 7:
519  case 8:
520  case 9:
521  case 10:
522  case 11:
523  case 12:
524  {
525  pzshape::TPZShapeLinear::ShapeInternal(x, order, phi, dphi);
526  }
527  break;
528  case 13:
529  {
530  pzshape::TPZShapeQuad::ShapeInternal(x, order, phi, dphi);
531  }
532  break;
533  case 14:
534  case 15:
535  case 16:
536  case 17:
537  {
538  pzshape::TPZShapeTriang::ShapeInternal(x, order, phi, dphi);
539  }
540  break;
541 
542  case 18:
543  {
544  ShapeInternal(x, order, phi, dphi);
545  }
546  break;
547  default:
548  std::cout << __PRETTY_FUNCTION__ << " Wrong side parameter side " << side << std::endl;
549  DebugStop();
550  break;
551  }
552 
553 
554  }
555  void TPZShapePiram::TransformDerivativeFromRibToPiram(int rib,int num,TPZFMatrix<REAL> &dphi) {
556  for (int j = 0;j<num;j++) {
557  dphi(2,j) = gRibTrans3dPiram1d[rib][2]*dphi(0,j);
558  dphi(1,j) = gRibTrans3dPiram1d[rib][1]*dphi(0,j);
559  dphi(0,j) = gRibTrans3dPiram1d[rib][0]*dphi(0,j);
560  }
561  }
562 
563  void TPZShapePiram::TransformDerivativeFromFaceToPiram(int face,int num,TPZFMatrix<REAL> &dphi) {
564  for (int j = 0;j<num;j++) {
565  dphi(2,j) = gFaceTrans3dPiram2d[face][0][2]*dphi(0,j)+gFaceTrans3dPiram2d[face][1][2]*dphi(1,j);
566  REAL dphi1j = dphi(1,j);
567  dphi(1,j) = gFaceTrans3dPiram2d[face][0][1]*dphi(0,j)+gFaceTrans3dPiram2d[face][1][1]*dphi(1,j);
568  dphi(0,j) = gFaceTrans3dPiram2d[face][0][0]*dphi(0,j)+gFaceTrans3dPiram2d[face][1][0]*dphi1j;//dphi(1,j);
569  }
570  }
571  //transforma a derivada no ponto dentro da face
572  void TPZShapePiram::TransformDerivativeFace3dPiram(int transid, int face, int num, TPZFMatrix<REAL> &in) {
573  if (!face) TPZShapeQuad::TransformDerivative2dQ(transid,num,in);//face 13
574  else TPZShapeTriang::TransformDerivative2dT(transid,num,in);//outras
575  }
576 
577  //projeta o ponto do interior para o lado
578  void TPZShapePiram::ProjectPoint3dPiramToRib(int rib, TPZVec<REAL> &in, REAL &outval) {
579  outval = gRibTrans3dPiram1d[rib][0]*in[0]+gRibTrans3dPiram1d[rib][1]*in[1]+gRibTrans3dPiram1d[rib][2]*in[2];
580  }
581 
582  //projeta o ponto do interior para a face
583  void TPZShapePiram::ProjectPoint3dPiramToFace(int face, TPZVec<REAL> &in, TPZVec<REAL> &outval) {
584  outval[0] = gFaceTrans3dPiram2d[face][0][0]*in[0]+gFaceTrans3dPiram2d[face][0][1]*in[1]+gFaceTrans3dPiram2d[face][0][2]*in[2]+gFaceSum3dPiram2d[face][0];
585  outval[1] = gFaceTrans3dPiram2d[face][1][0]*in[0]+gFaceTrans3dPiram2d[face][1][1]*in[1]+gFaceTrans3dPiram2d[face][1][2]*in[2]+gFaceSum3dPiram2d[face][1];
586  }
587 
588  //transforma o ponto dentro da face
589  void TPZShapePiram::TransformPoint3dPiramFace(int transid, int face, TPZVec<REAL> &in, TPZVec<REAL> &out) {
590  if (!face) TPZShapeQuad::TransformPoint2dQ(transid,in,out);//face zero ou 13
591  else TPZShapeTriang::TransformPoint2dT(transid,in,out);//outras 14 a 17
592  }
593 
594  int TPZShapePiram::NConnectShapeF(int side, int order) {
595  if(side<5) return 1;//0 a 4
596  // int s = side-5;//s = 0 a 14 ou side = 6 a 20
597  if(side<13) return (order-1);//6 a 14
598  if(side==13) {
599  return ((order-1)*(order-1));
600  }
601  if(side<18) {//16,17,18
602  return ((order-2)*(order-1)/2);
603  }
604  if(side==18) {
605  int totsum = 0,sum;
606  for(int i=1;i<order-1;i++) {
607  sum = i*(i+1) / 2;
608  totsum += sum;
609  }
610  return totsum;
611  }
612  PZError << "TPZShapePiram::NConnectShapeF, bad parameter side " << side << endl;
613  return 0;
614  }
615 
616  int TPZShapePiram::NShapeF(TPZVec<int> &order) {
617  int in,res=NCornerNodes;
618  for(in=NCornerNodes;in<NSides;in++) res += NConnectShapeF(in,order[in-NCornerNodes]);
619  return res;
620  }
621 
622 
623  int TPZShapePiram::ClassId() const{
624  return Hash("TPZShapePiram") ^ pztopology::TPZPyramid::ClassId() << 1;
625  }
626 
627 };
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
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
static void ShapeInternal(TPZVec< REAL > &x, int ord, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, int transformation_index)
Computes the values of the orthogonal shapefunctions before multiplying them by the corner shapefunct...
groups all classes dedicated to the computation of shape functions
Definition: pzshapeextend.h:16
Defines PZError.
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
Free store vector implementation.
int ClassId() const override
Define the class id associated with the class.
void Shape(TPZVec< REAL > &pt, TPZVec< int > orders, TPZVec< TPZTransform< REAL > > &transvec, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi)
string res
Definition: test.py:151
Contains TPZShapePoint class which implements the shape function associated with a point...
static void ShapeInternal(TPZVec< REAL > &x, int order, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, int quad_transformation_index)
Compute the internal functions of the quadrilateral shape function at a point.
Implements generic class which holds a matrix of objects. Matrix.
Definition: pzshtmat.h:18
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
static void ShapeInternal(TPZVec< REAL > &x, int order, TPZFMatrix< REAL > &phi, TPZFMatrix< REAL > &dphi, int triangle_transformation_index)
Compute the internal functions of the triangle shape function at a point.
int64_t Cols() const
Definition: pzshtmat.h:47
int64_t Rows() const
Definition: pzshtmat.h:45
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
Contains TPZShapePiram class which implements the shape functions of a pyramid element.
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
Contains TPZShapeQuad class which implements the shape functions of a quadrilateral element...
Contains TPZShapeTriang class which implements the shape functions of a triangular element...
Non abstract class which implements full matrices with preallocated storage with (N+1) entries...
Definition: pzfmatrix.h:716
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15