NeoPZ
pzshapecube.cpp
Go to the documentation of this file.
1 
6 #include "pzshapecube.h"
7 #include "pzshapequad.h"
8 #include "pzshapepoint.h"
9 #include "pzshapelinear.h"
10 #include "pzmanvector.h"
11 #include "pzerror.h"
12 #include "pzreal.h"
13 
14 using namespace std;
15 
16 namespace pzshape {
17 
18  REAL TPZShapeCube::gRibTrans3dCube1d[12][3] = {
19  {1., 0.,0.} , { 0.,1.,0.} , {-1., 0.,0.} ,
20  {0.,-1.,0.} , { 0.,0.,1.} , { 0., 0.,1.} ,
21  {0., 0.,1.} , { 0.,0.,1.} , { 1., 0.,0.} ,
22  {0., 1.,0.} , {-1.,0.,0.} , { 0.,-1.,0.}
23  };
24 
25  REAL TPZShapeCube::gFaceTrans3dCube2d[6][2][3] = {
26  { { 1.,0.,0.},{0.,1.,0.} },
27  { { 1.,0.,0.},{0.,0.,1.} },
28  { { 0.,1.,0.},{0.,0.,1.} },
29  { { 1.,0.,0.},{0.,0.,1.} },//{-1.,0.,0.},{0.,0.,1.}
30  { { 0.,1.,0.},{0.,0.,1.} },
31  { { 1.,0.,0.},{0.,1.,0.} }
32  };
33 
34  void TPZShapeCube::ShapeCorner(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
35 
36  REAL x[2],dx[2],y[2],dy[2],z[2],dz[2];
37  x[0] = (1.-pt[0])/2.;
38  x[1] = (1.+pt[0])/2.;
39  dx[0] = -0.5;
40  dx[1] = 0.5;
41  y[0] = (1.-pt[1])/2.;
42  y[1] = (1.+pt[1])/2.;
43  dy[0] = -0.5;
44  dy[1] = 0.5;
45  z[0] = (1.-pt[2])/2.;
46  z[1] = (1.+pt[2])/2.;
47  dz[0] = -0.5;
48  dz[1] = 0.5;
49 
50  phi(0,0) = x[0]*y[0]*z[0];
51  phi(1,0) = x[1]*y[0]*z[0];
52  phi(2,0) = x[1]*y[1]*z[0];
53  phi(3,0) = x[0]*y[1]*z[0];
54  phi(4,0) = x[0]*y[0]*z[1];
55  phi(5,0) = x[1]*y[0]*z[1];
56  phi(6,0) = x[1]*y[1]*z[1];
57  phi(7,0) = x[0]*y[1]*z[1];
58  dphi(0,0) = dx[0]*y[0]*z[0];
59  dphi(1,0) = x[0]*dy[0]*z[0];
60  dphi(2,0) = x[0]*y[0]*dz[0];
61  dphi(0,1) = dx[1]*y[0]*z[0];
62  dphi(1,1) = x[1]*dy[0]*z[0];
63  dphi(2,1) = x[1]*y[0]*dz[0];
64  dphi(0,2) = dx[1]*y[1]*z[0];
65  dphi(1,2) = x[1]*dy[1]*z[0];
66  dphi(2,2) = x[1]*y[1]*dz[0];
67  dphi(0,3) = dx[0]*y[1]*z[0];
68  dphi(1,3) = x[0]*dy[1]*z[0];
69  dphi(2,3) = x[0]*y[1]*dz[0];
70  dphi(0,4) = dx[0]*y[0]*z[1];
71  dphi(1,4) = x[0]*dy[0]*z[1];
72  dphi(2,4) = x[0]*y[0]*dz[1];
73  dphi(0,5) = dx[1]*y[0]*z[1];
74  dphi(1,5) = x[1]*dy[0]*z[1];
75  dphi(2,5) = x[1]*y[0]*dz[1];
76  dphi(0,6) = dx[1]*y[1]*z[1];
77  dphi(1,6) = x[1]*dy[1]*z[1];
78  dphi(2,6) = x[1]*y[1]*dz[1];
79  dphi(0,7) = dx[0]*y[1]*z[1];
80  dphi(1,7) = x[0]*dy[1]*z[1];
81  dphi(2,7) = x[0]*y[1]*dz[1];
82  }
83 
90  void TPZShapeCube::ShapeGenerating(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi)
91  {
92  int is;
93  // contribute the ribs
94  for(is=8; is<20; is++)
95  {
96  int is1,is2;
97  is1 = ContainedSideLocId(is,0);
98  is2 = ContainedSideLocId(is,1);
99  phi(is,0) = phi(is1,0)*phi(is2,0);
100  dphi(0,is) = dphi(0,is1)*phi(is2,0)+phi(is1,0)*dphi(0,is2);
101  dphi(1,is) = dphi(1,is1)*phi(is2,0)+phi(is1,0)*dphi(1,is2);
102  dphi(2,is) = dphi(2,is1)*phi(is2,0)+phi(is1,0)*dphi(2,is2);
103  }
104  // contribution of the faces
105  for(is=20; is<26; is++)
106  {
107  int is1,is2;
108  is1 = ContainedSideLocId(is,0);
109  is2 = ContainedSideLocId(is,2);
110  phi(is,0) = phi(is1,0)*phi(is2,0);
111  dphi(0,is) = dphi(0,is1)*phi(is2,0)+phi(is1,0)*dphi(0,is2);
112  dphi(1,is) = dphi(1,is1)*phi(is2,0)+phi(is1,0)*dphi(1,is2);
113  dphi(2,is) = dphi(2,is1)*phi(is2,0)+phi(is1,0)*dphi(2,is2);
114  }
115  // contribution of the volume
116  for(is=26; is<27; is++)
117  {
118  int is1,is2;
119  is1 = 0;
120  is2 = 6;
121  phi(is,0) = phi(is1,0)*phi(is2,0);
122  dphi(0,is) = dphi(0,is1)*phi(is2,0)+phi(is1,0)*dphi(0,is2);
123  dphi(1,is) = dphi(1,is1)*phi(is2,0)+phi(is1,0)*dphi(1,is2);
124  dphi(2,is) = dphi(2,is1)*phi(is2,0)+phi(is1,0)*dphi(2,is2);
125  }
126 
127  // Make the generating shape functions linear and unitary
128  for(is=8; is<27; is++)
129  {
130  TPZStack<int> highsides;
131  HigherDimensionSides(is,highsides);
132  int h, nh = highsides.NElements();
133  for(h=0; h<nh; h++)
134  {
135  int hs = highsides[h];
136  phi(is,0) += phi(hs,0);
137  dphi(0,is) += dphi(0,hs);
138  dphi(1,is) += dphi(1,hs);
139  dphi(2,is) += dphi(2,hs);
140  }
141  int dim = SideDimension(is);
142  int mult = (dim == 1) ? 4 : (dim == 2) ? 16 : (dim ==3) ? 64 : 0;
143  phi(is,0) *= mult;
144  dphi(0,is) *= mult;
145  dphi(1,is) *= mult;
146  dphi(2,is) *= mult;
147  }
148 
149  }
150 
151 
152  void TPZShapeCube::Shape(TPZVec<REAL> &pt, TPZVec<int64_t> &id, TPZVec<int> &order, TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
153  ShapeCorner(pt,phi,dphi);
154  bool linear = true;
155  int is,d;
156  for(is=NCornerNodes; is<NSides; is++) if(order[is-NCornerNodes] > 1) linear = false;
157  if(linear) return;
158 
159  TPZFNMatrix<100> phiblend(NSides,1),dphiblend(Dimension,NSides);
160  for(is=0; is<NCornerNodes; is++)
161  {
162  phiblend(is,0) = phi(is,0);
163  for(d=0; d<Dimension; d++)
164  {
165  dphiblend(d,is) = dphi(d,is);
166  }
167  }
168  ShapeGenerating(pt,phiblend,dphiblend);
169  int shape = 8;
170  //rib shapes
171  for (int rib = 0; rib < 12; rib++) {
172  REAL outval;
173  ProjectPoint3dCubeToRib(rib,pt,outval);
174  TPZManVector<REAL,1> outvalvec(1,outval);
175  TPZVec<int64_t> ids(2);
176  int id0,id1;
177  id0 = SideNodes[rib][0];
178  id1 = SideNodes[rib][1];
179  ids[0] = id[id0];
180  ids[1] = id[id1];
181  REAL store1[20],store2[60];
182  int nshape = order[rib]-1;//three orders : order in x , order in y and order in z
183  TPZFMatrix<REAL> phin(nshape,1,store1,20),dphin(3,nshape,store2,60);
184  phin.Zero();
185  dphin.Zero();
186  TPZShapeLinear::ShapeInternal(outvalvec,order[rib],phin,dphin,TPZShapeLinear::GetTransformId1d(ids));//ordin = ordem de um lado
187  TransformDerivativeFromRibToCube(rib,nshape,dphin);
188  for (int i = 0; i < nshape; i++) {
189  phi(shape,0) = phiblend(rib+8,0)*phin(i,0);
190  for(int xj=0;xj<3;xj++) {
191  dphi(xj,shape) = dphiblend(xj ,rib+8) * phin( i, 0) +
192  phiblend(rib+8, 0 ) * dphin(xj,i);
193  }
194  shape++;
195  }
196  }
197  //face shapes
198  for (int face = 0; face < 6; face++) {
199 
200  TPZVec<REAL> outval(2);
201  ProjectPoint3dCubeToFace(face,pt,outval);
202  REAL store1[20],store2[60];
203  int ord1,ord2;
204  ord1 = order[12+face];
205  ord2=ord1;
206  // FaceOrder(face,ord1,ord2);
207  if(ord1<2 || ord2<2) continue;
208  int ord = (ord1-1)*(ord2-1);
209  TPZFMatrix<REAL> phin(ord,1,store1,20),dphin(3,ord,store2,60);//ponto na face
210  phin.Zero();
211  dphin.Zero();
212  int ordin = (ord1 > ord2) ? ord1 : ord2;
213  ordin--;
214  TPZManVector<int64_t> ids(4);
215  //TPZVec<int> ids(4);
216  // int id0,id1;
217  int i;
218  for(i=0;i<4;i++) ids[i] = id[FaceNodes[face][i]];
219  //id0 = ShapeFaceId[face][0];//numero das shapes da face que compoem a shape atual
220  //id1 = ShapeFaceId[face][1];
221  TPZShapeQuad::ShapeInternal(outval,ord1-2,phin,dphin,TPZShapeQuad::GetTransformId2dQ(ids));//ordin = ordem de um lado
222  TransformDerivativeFromFaceToCube(face,ord,dphin);//ord = numero de shapes
223  for(i=0;i<ord;i++) {
224  phi(shape,0) = phiblend(face+20,0)*phin(i,0);
225  for(int xj=0;xj<3;xj++) {
226  dphi(xj,shape) = dphiblend(xj,face+20) * phin(i ,0) +
227  phiblend(face+20, 0)* dphin(xj,i);
228  }
229  shape++;
230  }
231  }
232  //volume shapes
233  REAL store1[20],store2[60];
234  int ordmin1 = (order[18]-1);
235  int ord = ordmin1*ordmin1*ordmin1;//(p-1)^3 : 0<=n1,n2,n3<=p-2
236  TPZFNMatrix<20, REAL> phin(ord,1),dphin(3,ord);
237  phin.Zero();
238  dphin.Zero();
239  ShapeInternal(pt,order[18],phin,dphin);
240  for(int i=0;i<ord;i++) {
241  phi(shape,0) = phiblend(26,0)*phin(i,0);
242  for(int xj=0;xj<3;xj++) {
243  dphi(xj,shape) = dphiblend(xj,26) * phin(i ,0) +
244  phiblend(26, 0)*dphin(xj,i);
245  }
246  shape++;
247  }
248  }
249 
250  void TPZShapeCube::SideShape(int side, TPZVec<REAL> &pt, TPZVec<int64_t> &id, TPZVec<int> &order, TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
251  if(side<0 || side>26) PZError << "TPZCompElC3d::SideShapeFunction. Bad paramenter side.\n";
252  else if(side==26) Shape(pt,id,order,phi,dphi);
253  else if(side<8) TPZShapePoint::Shape(pt,id,order,phi,dphi);
254  else if(side<20) {//8 a 19
255  TPZShapeLinear::Shape(pt,id,order,phi,dphi);
256  }
257  else if(side<26) {//faces
258  TPZShapeQuad::Shape(pt,id,order,phi,dphi);
259  }
260  else
261  {
262  Shape(pt,id,order,phi,dphi);
263  }
264  }
265 
266  void TPZShapeCube::ShapeOrder(TPZVec<int64_t> &id, TPZVec<int> &order, TPZGenMatrix<int> &shapeorders)//, TPZVec<int64_t> &sides
267  {
268  //DebugStop();
269 
270  int64_t nsides = TPZShapeCube::NSides;
271  int nshape = NCornerNodes;
272 
273  int linha = 0;
274  for (int side = 0; side < nsides; side++)
275  {
276 
277  nshape = 1;
278  if(side >= NCornerNodes) nshape = NConnectShapeF(side,order[side-NCornerNodes]);
279  int sideorder = 1;
280  if(side >= NCornerNodes) sideorder = order[side-NCornerNodes];
281 
282  TPZGenMatrix<int> locshapeorders(nshape,3);
283  SideShapeOrder(side, id, sideorder, locshapeorders);
284 
285  int nlin = locshapeorders.Rows();
286  int ncol = locshapeorders.Cols();
287 
288  for (int il = 0; il<nlin; il++)
289  {
290  for (int jc = 0; jc<ncol; jc++)
291  {
292  shapeorders(linha, jc) = locshapeorders(il, jc);
293  }
294  linha++;
295  }
296  }
297  }
298 
299 
300  void TPZShapeCube::SideShapeOrder(int side, TPZVec<int64_t> &id, int order, TPZGenMatrix<int> &shapeorders)
301  {
302  //DebugStop();
303  if (side<=7)
304  {
305  if (shapeorders.Rows() != 1)
306  {
307  DebugStop();
308  }
309  shapeorders(0,0) = 1;
310  shapeorders(0,1) = 0;
311  shapeorders(0,2) = 0;
312  }
313  else if (side == 26)
314  {
315  int nsei = (order-1)*(order-1)*(order-1);
316  if (shapeorders.Rows() != nsei) {
317  DebugStop();
318  }
319  int count = 0;
320  for (int i=2; i<order+1; i++) {
321  for (int j=2; j<order+1; j++) {
322  for (int k=2; k<order+1; k++) {
323  int a = i;
324  int b = j;
325  int c = k;
326  shapeorders(count,0) = a;
327  shapeorders(count,1) = b;
328  shapeorders(count,2) = c;
329  count++;
330 
331  }
332 
333  }
334  }
335  }
336  else if (side>7 && side<20)
337  {
338  int nshape = order-1;
339  if (shapeorders.Rows() != nshape)
340  {
341  DebugStop();
342  }
343  for (int ioy = 0; ioy < order-1; ioy++)
344  {
345  shapeorders(ioy,0) = ioy+2;
346  }
347  }
348  else
349  {
350 
351  if (shapeorders.Rows() != (order-1)*(order-1))
352  {
353  DebugStop();
354  }
355  TPZStack<int> lowersides;
356  LowerDimensionSides(side, lowersides);
357  lowersides.Push(side);
358 
359  //TPZVec<int> locsideorder(lowersides.size(),order);
360 
361  int nnodes = NSideNodes(side);
362 
363  TPZManVector<int64_t, 4> locid(nnodes);
364  for (int node=0; node<locid.size(); node++) {
365  locid[node] = id[ContainedSideLocId(side, node)];// SideNodeLocId( side, node);
366  }// sera que esta pegando os ids corretos mesmo?
367 
368  int nshape = (order-1)*(order-1);
369  TPZGenMatrix<int> locshapeorders(nshape,3);
370 
371 
372  TPZShapeQuad::SideShapeOrder(8,locid, order, locshapeorders);
373 
374  // temos que arrumar a saida de locshapeorders para adequar a orientacao dos vetores que geram
375  // a face do lado side
376 
377  // aqui o locshapeorders esta armazenado so para x e y
378  for (int il = 0; il<nshape; il++)
379  {
380  shapeorders(il, 0) = locshapeorders(il, 0);
381  shapeorders(il, 1) = locshapeorders(il, 1);
382  shapeorders(il, 2) = locshapeorders(il, 2);
383  }
384 
385  }
386 
387  }
388  void TPZShapeCube::ShapeInternal(TPZVec<REAL> &x, int order,TPZFMatrix<REAL> &phi,
389  TPZFMatrix<REAL> &dphi) {//,int cube_transformation_index
390  if((order-1) < 1) return;
391  int ord = order - 1;//fSideOrder[18]-1;
392  int nshape = ord*ord*ord;
393  phi.Resize(nshape,1);
394  dphi.Resize(3,nshape);
395  REAL store1[20],store2[20],store3[20],store4[20],store5[20],store6[20];
396  TPZFNMatrix<20, REAL> phi0(ord,1),phi1(ord,1),phi2(ord,1),
397  dphi0(1,ord),dphi1(1,ord),dphi2(1,ord);
398  TPZShapeLinear::fOrthogonal(x[0],ord,phi0,dphi0);
399  TPZShapeLinear::fOrthogonal(x[1],ord,phi1,dphi1);
400  TPZShapeLinear::fOrthogonal(x[2],ord,phi2,dphi2);
401  for (int i=0;i<ord;i++) {
402  for (int j=0;j<ord;j++) {
403  for (int k=0;k<ord;k++) {
404  int index = ord*(ord*i+j)+k;
405  phi(index,0) = phi0(i,0)* phi1(j,0)* phi2(k,0);
406  dphi(0,index) = dphi0(0,i)* phi1(j,0)* phi2(k,0);
407  dphi(1,index) = phi0(i,0)*dphi1(0,j)* phi2(k,0);
408  dphi(2,index) = phi0(i,0)* phi1(j,0)*dphi2(0,k);
409  }
410  }
411  }
412  }
413 
414 
415  void TPZShapeCube::ShapeInternal(int side, TPZVec<REAL> &x, int order,
416  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
417  if (side < 8 || side > 26) {
418  DebugStop();
419  }
420 
421  switch (side) {
422 
423  case 8:
424  case 9:
425  case 10:
426  case 11:
427  case 12:
428  case 13:
429  case 14:
430  case 15:
431  case 16:
432  case 17:
433  case 18:
434  case 19:
435  {
436  pzshape::TPZShapeLinear::ShapeInternal(x, order, phi, dphi);
437  }
438  break;
439  case 20:
440  case 21:
441  case 22:
442  case 23:
443  case 24:
444  case 25:
445  {
446  pzshape::TPZShapeQuad::ShapeInternal(x, order, phi, dphi);
447  }
448  break;
449  case 26:
450  {
451  ShapeInternal(x, order, phi, dphi);
452  }
453  break;
454  default:
455  std::cout << "Wrong side parameter side " << side << std::endl;
456  DebugStop();
457  break;
458  }
459 
460 
461 
462  }
463  void TPZShapeCube::TransformDerivativeFromRibToCube(int rib,int num,TPZFMatrix<REAL> &dphi) {
464  for (int j = 0;j<num;j++) {
465  dphi(2,j) = gRibTrans3dCube1d[rib][2]*dphi(0,j);
466  dphi(1,j) = gRibTrans3dCube1d[rib][1]*dphi(0,j);
467  dphi(0,j) = gRibTrans3dCube1d[rib][0]*dphi(0,j);
468  }
469  }
470 
471  void TPZShapeCube::ProjectPoint3dCubeToRib(int side, TPZVec<REAL> &in, REAL &outval) {
472  outval = gRibTrans3dCube1d[side][0]*in[0]+gRibTrans3dCube1d[side][1]*in[1]+gRibTrans3dCube1d[side][2]*in[2];
473  }
474 
475  void TPZShapeCube::ProjectPoint3dCubeToFace(int face, TPZVec<REAL> &in, TPZVec<REAL> &outval) {
476  outval[0] = gFaceTrans3dCube2d[face][0][0]*in[0]+gFaceTrans3dCube2d[face][0][1]*in[1]+gFaceTrans3dCube2d[face][0][2]*in[2];
477  outval[1] = gFaceTrans3dCube2d[face][1][0]*in[0]+gFaceTrans3dCube2d[face][1][1]*in[1]+gFaceTrans3dCube2d[face][1][2]*in[2];
478  }
479 
480  void TPZShapeCube::TransformDerivativeFromFaceToCube(int rib,int num,TPZFMatrix<REAL> &dphi) {
481 
482  for (int j = 0;j<num;j++) {
483  dphi(2,j) = gFaceTrans3dCube2d[rib][0][2]*dphi(0,j)+gFaceTrans3dCube2d[rib][1][2]*dphi(1,j);
484  REAL dphi1j = dphi(1,j);
485  dphi(1,j) = gFaceTrans3dCube2d[rib][0][1]*dphi(0,j)+gFaceTrans3dCube2d[rib][1][1]*dphi(1,j);
486  dphi(0,j) = gFaceTrans3dCube2d[rib][0][0]*dphi(0,j)+gFaceTrans3dCube2d[rib][1][0]*dphi1j;//dphi(1,j);
487  }
488  }
489 
490  void TPZShapeCube::ProjectPoint3dCubeSide(int side, TPZVec<REAL> &in, REAL &out) {
491 
492  out = gRibTrans3dCube1d[side][0]*in[0]+gRibTrans3dCube1d[side][1]*in[1]+gRibTrans3dCube1d[side][2]*in[2];
493  }
494 
495  void TPZShapeCube::ProjectPoint3dCubeFace(int face, TPZVec<REAL> &in, TPZVec<REAL> &out) {
496 
497  out[0] = gFaceTrans3dCube2d[face][0][0]*in[0]+gFaceTrans3dCube2d[face][0][1]*in[1]+gFaceTrans3dCube2d[face][0][2]*in[2];
498  out[1] = gFaceTrans3dCube2d[face][1][0]*in[0]+gFaceTrans3dCube2d[face][1][1]*in[1]+gFaceTrans3dCube2d[face][1][2]*in[2];
499  }
500 
501  int TPZShapeCube::NConnectShapeF(int side, int order){
502  if(side<8) return 1;//0 a 4
503  if(side<20) return (order-1);//6 a 14
504  if(side<26) {
505  return ((order-1)*(order-1));
506  }
507  if(side==26) {
508  return ((order-1)*(order-1)*(order-1));
509  }
510  PZError << "TPZShapeCube::NConnectShapeF, bad parameter side " << side << endl;
511  return 0;
512  }
513 
514  int TPZShapeCube::NShapeF(TPZVec<int> &order) {
515  int in,res = NCornerNodes;
516  for(in=NCornerNodes;in<NSides;in++) res += NConnectShapeF(in,order[in-NCornerNodes]);
517  return res;
518  }
519 
520 #ifdef _AUTODIFF
521 
522  void TPZShapeCube::ShapeCube(TPZVec<REAL> &point, TPZVec<int64_t> &id, TPZVec<int> &order, TPZVec<FADREAL> &phi)
523  {
524  const int ndim = 3;
525 
526  TPZVec<FADREAL> pt(3);
527  pt[0] = point[0];
528  pt[0].diff(0, ndim);
529 
530  pt[1] = point[1];
531  pt[1].diff(1, ndim);
532 
533  pt[2] = point[2];
534  pt[2].diff(2, ndim);
535 
536  ShapeCornerCube(pt,phi);
537 
538  int shape = 8;
539  //rib shapes
540  for (int rib = 0; rib < 12; rib++) {
541  FADREAL outval(ndim, 0.0);
542  ProjectPoint3dCubeToRib(rib,pt,outval);
543  TPZVec<int64_t> ids(2);
544  int id0,id1;
545  id0 = SideNodes[rib][0];
546  id1 = SideNodes[rib][1];
547  ids[0] = id[id0];
548  ids[1] = id[id1];
549  //REAL store1[20], store2[60];
550  int ordin = order[rib]-1;//three orders : order in x , order in y and order in z
551  //TPZFMatrix<REAL> phin(ordin,1,store1,20),dphin(3,ordin,store2,60);
552  //phin.Zero();
553  //dphin.Zero();
554  TPZVec<FADREAL> phin(20, FADREAL(ndim, 0.0)); //3d
555  TPZShapeLinear::ShapeInternal(outval,ordin,phin,TPZShapeLinear::GetTransformId1d(ids));//ordin = ordem de um lado
556  // TransformDerivativeFromRibToCube(rib,ordin,phin);
557  for (int i = 0; i < ordin; i++) {
558  //phi(shape,0) = phi(id0,0)*phi(id1,0)*phin(i,0);
559  phi[shape] = phi[id0] * phi[id1] * phin[i];
560  /*for(int xj=0;xj<3;xj++) {
561  dphi(xj,shape) = dphi(xj ,id0) * phi(id1, 0 ) * phin( i, 0) +
562  phi(id0, 0 ) * dphi(xj ,id1) * phin( i, 0) +
563  phi(id0, 0 ) * phi(id1, 0 ) * dphin(xj,i);
564  }*/ // implicitly done
565  shape++;
566  }
567 
568  }
569  //face shapes
570  for (int face = 0; face < 6; face++) {
571 
572  //TPZVec<REAL> outval(2);
573  TPZVec<FADREAL> outval(2, FADREAL(ndim, 0.0));
574  ProjectPoint3dCubeToFace(face,pt,outval);
575  // REAL store1[20],store2[60];
576  int ord1,ord2;
577  ord1 = order[12+face];
578  ord2=ord1;
579  // FaceOrder(face,ord1,ord2); // already commented in the non-FAD version
580  if(ord1<2 || ord2<2) continue;
581  int ord = (ord1-1)*(ord2-1);
582  //TPZFMatrix<REAL> phin(ord,1,store1,20),dphin(3,ord,store2,60);//ponto na face
583  TPZVec<FADREAL> phin(20, FADREAL(ndim, 0.0)); //3d
584  //phin.Zero();
585  //dphin.Zero();
586  int ordin = (ord1 > ord2) ? ord1 : ord2;
587  ordin--;
588  TPZManVector<int64_t> ids(4);
589  //TPZVec<int> ids(4);
590  int id0,id1,i;
591  for(i=0;i<4;i++) ids[i] = id[FaceNodes[face][i]];
592  id0 = ShapeFaceId[face][0];//numero das shapes da face que compoem a shape atual
593  id1 = ShapeFaceId[face][1];
594  TPZShapeQuad::Shape2dQuadInternal(outval,ord1-2,phin,TPZShapeQuad::GetTransformId2dQ(ids));//ordin = ordem de um lado
595  // TransformDerivativeFromFaceToCube(face,ord,phin);//ord = numero de shapes
596  for(i=0;i<ord;i++) {
597  // phi(shape,0) = phi(id0,0)*phi(id1,0)*phin(i,0);
598  phi[shape] = phi[id0] * phi[id1] * phin[i];
599  /* for(int xj=0;xj<3;xj++) {
600  dphi(xj,shape) = dphi(xj,id0)* phi(id1 , 0 )* phin(i ,0) +
601  phi(id0, 0)*dphi(xj ,id1)* phin(i ,0) +
602  phi(id0, 0)* phi(id1 , 0 )*dphin(xj,i); // implicitly done
603  }*/
604  shape++;
605  }
606  }
607 
608  //volume shapes
609  //REAL store1[20],store2[60];
610  int ordmin1 = (order[18]-1);
611  int ord = ordmin1*ordmin1*ordmin1;//(p-1)^3 : 0<=n1,n2,n3<=p-2
612  //TPZFMatrix<REAL> phin(ord,1,store1,20),dphin(3,ord,store2,60);
613  TPZVec<FADREAL> phin(20, FADREAL(ndim, 0.0)); //3d
614  //phin.Zero();
615  //dphin.Zero();
616  Shape3dCubeInternal(pt,ordmin1,phin);
617  for(int i=0;i<ord;i++) {
618  //phi(shape,0) = phi(0,0)*phi(6,0)*phin(i,0);
619  phi[shape] = phi[0] * phi[6] * phin[i];
620  shape++;
621  }
622  }
623 
624 
625  void TPZShapeCube::ShapeCornerCube(TPZVec<FADREAL> &pt, TPZVec<FADREAL> &phi)
626  {
627  FADREAL x[2], y[2], z[2];
628 
629  x[0] = (REAL(1.)-pt[0])/REAL(2.);
630  x[1] = (REAL(1.)+pt[0])/REAL(2.);
631  y[0] = (REAL(1.)-pt[1])/REAL(2.);
632  y[1] = (REAL(1.)+pt[1])/REAL(2.);
633  z[0] = (REAL(1.)-pt[2])/REAL(2.);
634  z[1] = (REAL(1.)+pt[2])/REAL(2.);
635 
636  phi[0] = x[0]*y[0]*z[0];
637  phi[1] = x[1]*y[0]*z[0];
638  phi[2] = x[1]*y[1]*z[0];
639  phi[3] = x[0]*y[1]*z[0];
640  phi[4] = x[0]*y[0]*z[1];
641  phi[5] = x[1]*y[0]*z[1];
642  phi[6] = x[1]*y[1]*z[1];
643  phi[7] = x[0]*y[1]*z[1];
644  }
645 
646  void TPZShapeCube::ProjectPoint3dCubeToRib(int side, TPZVec<FADREAL> &in, FADREAL &outval)
647  {
648  outval = gRibTrans3dCube1d[side][0]*in[0]+gRibTrans3dCube1d[side][1]*in[1]+gRibTrans3dCube1d[side][2]*in[2];
649  }
650 
651  void TPZShapeCube::ProjectPoint3dCubeToFace(int face, TPZVec<FADREAL> &in, TPZVec<FADREAL> &outval) {
652  outval[0] = gFaceTrans3dCube2d[face][0][0]*in[0]+gFaceTrans3dCube2d[face][0][1]*in[1]+gFaceTrans3dCube2d[face][0][2]*in[2];
653  outval[1] = gFaceTrans3dCube2d[face][1][0]*in[0]+gFaceTrans3dCube2d[face][1][1]*in[1]+gFaceTrans3dCube2d[face][1][2]*in[2];
654  }
655 
656  void TPZShapeCube::Shape3dCubeInternal(TPZVec<FADREAL> &x, int order,TPZVec<FADREAL> &phi)
657  {
658  const int ndim = 3;
659 
660  if(order < 1) return;
661  int ord = order;//fSideOrder[18]-1;
662  order = order*order*order;
663  phi.Resize(order, FADREAL(ndim, 0.0));
664  TPZVec<FADREAL> phi0(20, FADREAL(ndim, 0.0)),
665  phi1(20, FADREAL(ndim, 0.0)),
666  phi2(20, FADREAL(ndim, 0.0));
667  TPZShapeLinear::FADfOrthogonal(x[0],ord,phi0);
668  TPZShapeLinear::FADfOrthogonal(x[1],ord,phi1);
669  TPZShapeLinear::FADfOrthogonal(x[2],ord,phi2);
670  for (int i=0;i<ord;i++) {
671  for (int j=0;j<ord;j++) {
672  for (int k=0;k<ord;k++) {
673  int index = ord*(ord*i+j)+k;
674  //phi(index,0) = phi0(i,0)* phi1(j,0)* phi2(k,0);
675  phi[index] = phi0[i] * phi1[j] * phi2[k];
676  /*dphi(0,index) = dphi0(0,i)* phi1(j,0)* phi2(k,0);
677  dphi(1,index) = phi0(i,0)*dphi1(0,j)* phi2(k,0);
678  dphi(2,index) = phi0(i,0)* phi1(j,0)*dphi2(0,k);
679  */
680  }
681  }
682  }
683  }
684 
685 #endif
686 
687 };
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)
pzshapelinear.h
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
pzshapecube.h
Contains TPZShapeCube class which implements the shape functions of a hexaedral element.
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
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
dx
expr_ dx(i) *cos(expr_.val())
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
TPZGenMatrix::Cols
int64_t Cols() const
Definition: pzshtmat.h:47
TPZGenMatrix::Rows
int64_t Rows() const
Definition: pzshtmat.h:45
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
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...
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:09 for NeoPZ by   doxygen 1.8.13