21 REAL TPZShapePiram::gRibTrans3dPiram1d[8][3] = {
22 { 1., 0.,0.} , { 0., 1.,0.} ,
23 {-1., 0.,0.} , { 0.,-1.,0.} ,
24 { .5, .5,1.} , {-.5,.5 ,1.} ,
25 {-.5,-.5,1.} , { .5,-.5,1.}
30 REAL TPZShapePiram::gFaceTrans3dPiram2d[5][2][3] = {
31 { { 1., 0., 0.},{ 0., 1.,0.} },
32 { { 1.,-.5,-.5},{ 0., 1.,1.} },
33 { { .5, 1.,-.5},{-1., 0.,1.} },
34 { { 1., .5,-.5},{ 0.,-1.,1.} },
35 { {-.5, 1.,-.5},{1.,0.,1.} }
38 REAL TPZShapePiram::gFaceSum3dPiram2d[5][2] = { {.0,0.},{-.5,0.},{-.5,0.},{-.5,0.},{-.5,0.} };
51 for(is=NCornerNodes; is<NCornerNodes+8; is++)
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);
66 is1 = ShapeFaceId[0][0];
67 is2 = ShapeFaceId[0][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);
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);
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);
99 for(is=NCornerNodes; is<NCornerNodes+4; is++)
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);
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++)
112 phi(is,0) *= sidescale[is];
113 dphi(0,is) *= sidescale[is];
114 dphi(1,is) *= sidescale[is];
115 dphi(2,is) *= sidescale[is];
124 CornerShape(pt,phi,dphi);
127 for(is=NCornerNodes; is<NSides; is++){
128 if(order[is-NCornerNodes] > 1){
135 for(is=0; is<NCornerNodes; is++)
137 phiblend(is,0) = phi(is,0);
138 for(d=0; d<Dimension; d++)
140 dphiblend(d,is) = dphi(d,is);
143 ShapeGenerating(pt,phiblend,dphiblend);
147 for (
int rib = 0; rib < 8; rib++) {
148 if (order[rib] <2 )
continue;
150 ProjectPoint3dPiramToRib(rib,pt,outval);
154 id0 = SideNodes[rib][0];
155 id1 = SideNodes[rib][1];
158 REAL store1[20],store2[60];
159 int nshape = order[rib]-1;
163 TPZShapeLinear::ShapeInternal(outvalvec,order[rib],phin,dphin,TPZShapeLinear::GetTransformId1d(ids));
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);
175 for (
int face = 0; face < 5; face++) {
176 if (order[face+8] < 2)
continue;
177 if(face>0 && order[face+8]<=2)
continue;
179 ProjectPoint3dPiramToFace(face,pt,outval);
181 ord1 = order[8+face];
184 if(face && ord1<3)
continue;
186 if(!face) nshape = (ord1-1)*(ord1-1);
187 else nshape = (ord1-2)*(ord1-1)/2;
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]];
201 transid = TPZShapeQuad::GetTransformId2dQ(ids);
202 TPZShapeQuad::ShapeInternal(outval,ord1-2,phin,dphin,transid);
205 transid = TPZShapeTriang::GetTransformId2dT(ids);
206 outval[0] = (outval[0]+1.)/2.;
207 outval[1] = (outval[1]+1.)/2.;
208 TPZShapeTriang::ShapeInternal(outval,ord1-2,phin,dphin,transid);
209 int c = dphin.
Cols();
210 for(
int i=0;i<c;i++) {
216 TransformDerivativeFromFaceToPiram(face,nshape,dphin);
217 for(i=0;i<nshape;i++) {
218 phi(shape,0) = phiblend(face+13,0)*phin(i,0);
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);
226 if(order[13]<3)
return;
229 for(i=0;i<order[13]-2;i++) {
230 nshape += (i+1)*(i+2) / 2;
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);
247 if(
fabs(pt[0])<1.e-10 &&
fabs(pt[1])<1.e-10 && pt[2]==1.) {
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;
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;
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;
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;
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;
312 if(side<0 || side>18)
PZError <<
"TPZCompElPi3d::SideShapeFunction. Bad paramenter side.\n";
313 else if(side==18)
Shape(point,
id,order,phi,dphi);
317 }
else if(side == 13) {
328 int64_t nsides = TPZShapePiram::NSides;
332 for (
int side = 0; side < nsides; side++)
336 if(side >= NCornerNodes) nshape = NConnectShapeF(side,order[side-NCornerNodes]);
338 if(side >= NCornerNodes) sideorder = order[side-NCornerNodes];
341 SideShapeOrder(side,
id, sideorder, locshapeorders);
343 int nlin = locshapeorders.
Rows();
344 int ncol = locshapeorders.
Cols();
346 for (
int il = 0; il<nlin; il++)
348 for (
int jc = 0; jc<ncol; jc++)
350 shapeorders(linha, jc) = locshapeorders(il, jc);
363 if (shapeorders.
Rows() != 1)
367 shapeorders(0,0) = 1;
368 shapeorders(0,1) = 0;
369 shapeorders(0,2) = 0;
371 else if (side>=5 && side<=12)
373 int nshape = order-1;
374 if (shapeorders.
Rows() != nshape)
378 for (
int ioy = 0; ioy < order-1; ioy++)
380 shapeorders(ioy,0) = ioy+2;
386 if (shapeorders.
Rows() != (order-1)*(order-1))
391 LowerDimensionSides(side, lowersides);
392 lowersides.
Push(side);
394 int nnodes = NSideNodes(side);
397 for (
int node=0; node<locid.
size(); node++) {
398 locid[node] =
id[ContainedSideLocId(side, node)];
401 int nshape = (order-1)*(order-1);
405 TPZShapeQuad::SideShapeOrder(8,locid, order, locshapeorders);
411 for (
int il = 0; il<nshape; il++)
413 shapeorders(il, 0) = locshapeorders(il, 0);
414 shapeorders(il, 1) = locshapeorders(il, 1);
415 shapeorders(il, 2) = locshapeorders(il, 2);
418 else if (side >= 14 && side <= 17)
421 int nshape = (order-2)*(order-1)/2;
422 if (shapeorders.
Rows() != nshape)
427 LowerDimensionSides(side, lowersides);
428 lowersides.
Push(side);
430 int nnodes = NSideNodes(side);
433 for (
int node=0; node<locid.
size(); node++) {
434 locid[node] =
id[ContainedSideLocId(side, node)];
440 TPZShapeTriang::SideShapeOrder(6,locid, order, locshapeorders);
442 for (
int il = 0; il<nshape; il++)
444 shapeorders(il, 0) = locshapeorders(il, 0);
445 shapeorders(il, 1) = locshapeorders(il, 1);
446 shapeorders(il, 2) = locshapeorders(il, 2);
453 for(
int i=0;i<order - 2;i++) {
454 totsum += (i+1)*(i+2) / 2;
458 if (shapeorders.
Rows() != nshape) {
463 for (
int i=0;i<ord;i++) {
464 for (
int j=0;j<ord;j++) {
465 for (
int k=0;k<ord;k++) {
469 shapeorders(count,0) = 3 + soma;
470 shapeorders(count,1) = 3 + soma;
471 shapeorders(count,2) = 3 + soma;
484 if(order < 3)
return;
487 dphi0(1,ord),dphi1(1,ord),dphi2(1,ord);
488 TPZShapeLinear::fOrthogonal(x[0],ord,phi0,dphi0);
489 TPZShapeLinear::fOrthogonal(x[1],ord,phi1,dphi1);
490 TPZShapeLinear::fOrthogonal(2.*x[2]-1.,ord,phi2,dphi2);
492 for (
int i=0;i<ord;i++) {
493 for (
int j=0;j<ord;j++) {
494 for (
int k=0;k<ord;k++) {
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);
508 void TPZShapePiram::ShapeInternal(
int side,
TPZVec<REAL> &x,
int order,
510 if (side < 5 || side > 18) {
544 ShapeInternal(x, order, phi, dphi);
548 std::cout << __PRETTY_FUNCTION__ <<
" Wrong side parameter side " << side << std::endl;
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);
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;
572 void TPZShapePiram::TransformDerivativeFace3dPiram(
int transid,
int face,
int num,
TPZFMatrix<REAL> &in) {
573 if (!face) TPZShapeQuad::TransformDerivative2dQ(transid,num,in);
574 else TPZShapeTriang::TransformDerivative2dT(transid,num,in);
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];
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];
590 if (!face) TPZShapeQuad::TransformPoint2dQ(transid,in,out);
591 else TPZShapeTriang::TransformPoint2dT(transid,in,out);
594 int TPZShapePiram::NConnectShapeF(
int side,
int order) {
597 if(side<13)
return (order-1);
599 return ((order-1)*(order-1));
602 return ((order-2)*(order-1)/2);
606 for(
int i=1;i<order-1;i++) {
612 PZError <<
"TPZShapePiram::NConnectShapeF, bad parameter side " << side << endl;
617 int in,
res=NCornerNodes;
618 for(in=NCornerNodes;in<NSides;in++) res += NConnectShapeF(in,order[in-NCornerNodes]);
623 int TPZShapePiram::ClassId()
const{
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
Implements a vector class which allows to use external storage provided by the user. Utility.
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
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
int Zero() override
Makes Zero all the elements.
int64_t size() const
Returns the number of elements of the vector.
void Push(const T object)
Pushes a copy of the object on the stack.
#define DebugStop()
Returns a message to user put a breakpoint in.
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)
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.
int32_t Hash(std::string str)
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
Returns number of cols.
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...
#define PZError
Defines the output device to error messages and the DebugStop() function.