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