NeoPZ
pzshapequad.cpp
Go to the documentation of this file.
1 
6 #include "pzshapequad.h"
7 #include "pzshapelinear.h"
8 #include "pzshapepoint.h"
9 #include "pzmanvector.h"
10 #include "pzerror.h"
11 #include "pzreal.h"
12 #include "pzlog.h"
13 
14 #ifdef LOG4CXX
15 static LoggerPtr logger(Logger::getLogger("pz.shape.TPZShapeQuad"));
16 #endif
17 
18 using namespace std;
19 
20 namespace pzshape {
21 
23  REAL TPZShapeQuad::gTrans2dQ[8][2][2] = {//s* , t*
24  { { 1., 0.},{ 0., 1.} },
25  { { 0., 1.},{ 1., 0.} },
26  { { 0., 1.},{-1., 0.} },
27  { {-1., 0.},{ 0., 1.} },
28  { {-1., 0.},{ 0.,-1.} },//s* = -s t* = -t , etc
29  { { 0.,-1.},{-1., 0.} },
30  { { 0.,-1.},{ 1., 0.} },
31  { { 1., 0.},{ 0.,-1.} }
32  };
33 
34  REAL TPZShapeQuad::gRibTrans2dQ1d[4][2] = { {1.,0.},{0.,1.},{-1.,0.},{0.,-1.} };
35 
36  void TPZShapeQuad::ShapeCorner(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
37 
38  REAL x[2],dx[2],y[2],dy[2];
39  x[0] = (1.-pt[0])/2.;
40  x[1] = (1.+pt[0])/2.;
41  dx[0] = -0.5;
42  dx[1] = 0.5;
43  y[0] = (1.-pt[1])/2.;
44  y[1] = (1.+pt[1])/2.;
45  dy[0] = -0.5;
46  dy[1] = 0.5;
47  phi(0,0) = x[0]*y[0];
48  phi(1,0) = x[1]*y[0];
49  phi(2,0) = x[1]*y[1];
50  phi(3,0) = x[0]*y[1];
51  dphi(0,0) = dx[0]*y[0];
52  dphi(1,0) = x[0]*dy[0];
53  dphi(0,1) = dx[1]*y[0];
54  dphi(1,1) = x[1]*dy[0];
55  dphi(0,2) = dx[1]*y[1];
56  dphi(1,2) = x[1]*dy[1];
57  dphi(0,3) = dx[0]*y[1];
58  dphi(1,3) = x[0]*dy[1];
59  }
60 
61  /*
62  * Computes the generating shape functions for a quadrilateral element
63  * @param pt (input) point where the shape function is computed
64  * @param phi (input) value of the (4) shape functions
65  * @param dphi (input) value of the derivatives of the (4) shape functions holding the derivatives in a column
66  */
67  void TPZShapeQuad::ShapeGenerating(TPZVec<REAL> &pt, TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi)
68  {
69  int is;
70  for(is=4; is<8; is++)
71  {
72  phi(is,0) = phi(is%4,0)*phi((is+1)%4,0);
73  dphi(0,is) = dphi(0,is%4)*phi((is+1)%4,0)+phi(is%4,0)*dphi(0,(is+1)%4);
74  dphi(1,is) = dphi(1,is%4)*phi((is+1)%4,0)+phi(is%4,0)*dphi(1,(is+1)%4);
75  }
76  phi(8,0) = phi(0,0)*phi(2,0);
77  dphi(0,8) = dphi(0,0)*phi(2,0)+phi(0,0)*dphi(0,2);
78  dphi(1,8) = dphi(1,0)*phi(2,0)+phi(0,0)*dphi(1,2);
79 
80  // Make the generating shape functions linear and unitary
81  for(is=4; is<8; is++)
82  {
83  phi(is,0) += phi(8,0);
84  dphi(0,is) += dphi(0,8);
85  dphi(1,is) += dphi(1,8);
86  phi(is,0) *= 4.;
87  dphi(0,is) *= 4.;
88  dphi(1,is) *= 4.;
89  }
90  phi(8,0) *= 16.;
91  dphi(0,8) *= 16.;
92  dphi(1,8) *= 16.;
93  }
94 
97  ShapeCorner(pt,phi,dphi);
98  int is,d;
99  TPZFNMatrix<100> phiblend(NSides,1),dphiblend(Dimension,NSides);
100  for(is=0; is<NCornerNodes; is++)
101  {
102  phiblend(is,0) = phi(is,0);
103  for(d=0; d<Dimension; d++)
104  {
105  dphiblend(d,is) = dphi(d,is);
106  }
107  }
108  ShapeGenerating(pt,phiblend,dphiblend);
109  REAL out;
110  int shape = 4;
111  for (int rib = 0; rib < 4; rib++) {
112 
113  ProjectPoint2dQuadToRib(rib,pt,out);
114  TPZVec<int64_t> ids(2);
115  TPZManVector<REAL,1> outvec(1,out);
116  ids[0] = id[rib%4];
117  ids[1] = id[(rib+1)%4];
118  REAL store1[20],store2[40];
119  int ord2 = order[rib]-1;//two orders : order in x and order in y
120  TPZFMatrix<REAL> phin(ord2,1,store1,20),dphin(2,ord2,store2,40);
121  TPZShapeLinear::ShapeInternal(outvec,order[rib],phin,dphin,TPZShapeLinear::GetTransformId1d(ids));
122  TransformDerivativeFromRibToQuad(rib,ord2,dphin);
123  for (int i = 0; i < ord2; i++) {
124  phi(shape,0) = phiblend(rib+4,0)*phin(i,0);
125  for(int xj=0;xj<2;xj++) {
126  dphi(xj,shape) = dphiblend(xj,rib+4)*phin(i,0)+
127  phiblend(rib+4,0)*dphin(xj,i);
128  }
129  shape++;
130  }
131  }
132  REAL store1[20],store2[40];
133  int ord = (order[4]-1)*(order[4]-1);
134  TPZFMatrix<REAL> phin(ord,1,store1,20),dphin(2,ord,store2,40);
135  ShapeInternal(pt,order[4]-2,phin,dphin,GetTransformId2dQ(id));
136  for(int i=0;i<ord;i++) {//funcoes de interior s�o em numero ordem-1
137  phi(shape,0) = phiblend(8,0)*phin(i,0);
138  for(int xj=0;xj<2;xj++) {//x e y
139  dphi(xj,shape) = dphiblend(xj,8)*phin(i,0) +
140  phiblend(8,0)*dphin(xj,i);
141  }
142  shape++;
143  }
144  }
145 
146  void TPZShapeQuad::SideShape(int side,TPZVec<REAL> &pt, TPZVec<int64_t> &id, TPZVec<int> &order,
147  TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi) {
148  switch(side) {
149  case 0:
150  case 1:
151  case 2:
152  case 3:
153  TPZShapePoint::Shape(pt, id, order, phi, dphi);
154  break;
155  case 4:
156  case 5:
157  case 6:
158  case 7:
159  TPZShapeLinear::Shape(pt, id, order, phi, dphi);
160  break;
161  case 8:
162  Shape(pt, id, order, phi, dphi);
163  }
164  }
165 
166  void TPZShapeQuad::ShapeOrder(TPZVec<int64_t> &id, TPZVec<int> &order, TPZGenMatrix<int> &shapeorders)//, TPZVec<int64_t> &sides
167  {
168  int64_t nsides = TPZQuadrilateral::NSides;
169  // o que eh o vetor order?
170  // Eu suponho que em cada posicao tem a ordem de cada lado.
171  // Na shape ja esta associado a lados com dimensao maior que 1, order[0](lado 3) ...
172 
173  int nshape = NCornerNodes;
174 
175  for (int is = NCornerNodes; is < nsides; is++)
176  {
177  nshape += NConnectShapeF(is,order[is-NCornerNodes]);
178  }
179  // shapeorders ja vem com tamanho definido, entao a conta acima so serve para ver se as ordens estao corretas
180  if (nshape!=shapeorders.Rows())
181  {
182  // Algo errado
183  DebugStop();
184  }
185 
186  int linha = 0;
187  for (int side = 0; side < nsides; side++)
188  {
189 
190  nshape = 1;
191  if(side >= NCornerNodes) nshape = NConnectShapeF(side,order[side-NCornerNodes]);
192  int sideorder = 1;
193  if(side >= NCornerNodes) sideorder = order[side-NCornerNodes];
194 
195  TPZGenMatrix<int> locshapeorders(nshape,3);
196  SideShapeOrder(side, id, sideorder, locshapeorders);
197 
198  int nlin = locshapeorders.Rows();
199  int ncol = locshapeorders.Cols();
200 
201  for (int il = 0; il<nlin; il++)
202  {
203  for (int jc = 0; jc<ncol; jc++)
204  {
205  shapeorders(linha, jc) = locshapeorders(il, jc);
206  }
207  linha++;
208  }
209  }
210 
211  }
212 
213 
214  void TPZShapeQuad::SideShapeOrder(int side, TPZVec<int64_t> &id, int order, TPZGenMatrix<int> &shapeorders)
215  {
216 
217  if (side<=3)
218  {
219  if (shapeorders.Rows() != 1)
220  {
221  DebugStop();
222  }
223  shapeorders(0,0) = 1;
224  shapeorders(0,1) = 0;
225  shapeorders(0,2) = 0;
226  }
227  else if (side == 8)
228  {
229  if (shapeorders.Rows() != (order-1)*(order-1)) {
230  DebugStop();
231  }
232  int transid = GetTransformId(id);
233  int count = 0;
234  for (int i=2; i<order+1; i++) {
235  for (int j=2; j<order+1; j++) {
236  int a = i;
237  int b = j;
238  switch (transid)
239  {
240  case 0:
241  case 3:
242  case 4:
243  case 7:
244  break;
245  case 1:
246  case 2:
247  case 5:
248  case 6:
249  {
250  int c = a;
251  a = b;
252  b = c;
253  }
254  break;
255 
256  default:
257  DebugStop();
258  break;
259  }
260  shapeorders(count,0) = a;
261  shapeorders(count,1) = b;
262  count++;
263  }
264  }
265  }
266  else
267  {
268  int nshape = order-1;
269  if (shapeorders.Rows() != nshape)
270  {
271  DebugStop();
272  }
273  for (int ioy = 0; ioy < order-1; ioy++)
274  {
275  shapeorders(ioy,0) = ioy+2;
276  }
277  }
278 
279 
280  }
281  void TPZShapeQuad::ShapeInternal(TPZVec<REAL> &x, int order,
282  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
283 
284  if((order-2) < 0) return;
285  int ord1 = order - 1;
286  int numshape = (order-1)*(order-1);
287  if(numshape > phi.Rows() || phi.Cols() < 1) phi.Resize(numshape,1);
288  if(dphi.Rows() < 2 || dphi.Cols() < numshape) dphi.Resize(2,numshape);
289  TPZFNMatrix<20, REAL> phi0(ord1,1),phi1(ord1,1),dphi0(1,ord1),dphi1(1,ord1);
290  TPZShapeLinear::fOrthogonal(x[0],ord1,phi0,dphi0);
291  TPZShapeLinear::fOrthogonal(x[1],ord1,phi1,dphi1);
292  for (int i=0;i<ord1;i++) {
293  for (int j=0;j<ord1;j++) {
294  int index = i*ord1+j;
295  phi(index,0) = phi0(i,0)* phi1(j,0);
296  dphi(0,index) = dphi0(0,i)* phi1(j,0);
297  dphi(1,index) = phi0(i,0)*dphi1(0,j);
298  }
299  }
300 
301  }
302 
303  void TPZShapeQuad::ShapeInternal(TPZVec<REAL> &x, int order,
304  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi,int quad_transformation_index) {
305 
306  if(order < 0) return;
307  int ord1 = order+1;
308  int numshape = ord1*ord1;
309  TPZManVector<REAL> out(2);
310  TransformPoint2dQ(quad_transformation_index,x,out);
311 
312  if(numshape > phi.Rows() || phi.Cols() < 1) phi.Resize(numshape,1);
313  if(dphi.Rows() < 2 || dphi.Cols() < numshape) dphi.Resize(2,numshape);
314  REAL store1[20],store2[20],store3[20],store4[20];
315  TPZFMatrix<REAL> phi0(ord1,1,store1,20),phi1(ord1,1,store2,20),dphi0(1,ord1,store3,20),dphi1(1,ord1,store4,20);
316 
317  TPZShapeLinear::fOrthogonal(out[0],ord1,phi0,dphi0);
318  TPZShapeLinear::fOrthogonal(out[1],ord1,phi1,dphi1);
319  for (int i=0;i<ord1;i++) {
320  for (int j=0;j<ord1;j++) {
321  int index = i*ord1+j;
322  phi(index,0) = phi0(i,0)* phi1(j,0);
323  dphi(0,index) = dphi0(0,i)* phi1(j,0);
324  dphi(1,index) = phi0(i,0)*dphi1(0,j);
325  }
326  }
327  TransformDerivative2dQ(quad_transformation_index,numshape,dphi);
328  }
329 
330 
331 
332  void TPZShapeQuad::ShapeInternal(int side, TPZVec<REAL> &x, int order,
333  TPZFMatrix<REAL> &phi, TPZFMatrix<REAL> &dphi) {
334  if (side < 4 || side > 8) {
335  DebugStop();
336  }
337 
338  switch (side) {
339 
340  case 4:
341  case 5:
342  case 6:
343  case 7:
344  {
345  pzshape::TPZShapeLinear::ShapeInternal(x, order, phi, dphi);
346  }
347  break;
348  case 8:
349  {
350  ShapeInternal(x, order, phi, dphi);
351  }
352  break;
353  default:
354  LOGPZ_ERROR(logger,"Wrong side parameter")
355  DebugStop();
356  break;
357  }
358 
359 
360 
361  }
362 
363  void TPZShapeQuad::TransformDerivative2dQ(int transid, int num, TPZFMatrix<REAL> &in) {
364 
365  for(int i=0;i<num;i++) {
366  REAL aux[2];
367  aux[0] = in(0,i);
368  aux[1] = in(1,i);
369  in(0,i) = gTrans2dQ[transid][0][0]*aux[0]+gTrans2dQ[transid][1][0]*aux[1];
370  in(1,i) = gTrans2dQ[transid][0][1]*aux[0]+gTrans2dQ[transid][1][1]*aux[1];
371  }
372  }
373 
374  //transf. o ponto dentro da face quadrilateral
375  void TPZShapeQuad::TransformPoint2dQ(int transid, TPZVec<REAL> &in, TPZVec<REAL> &out) {
376 
377  out[0] = gTrans2dQ[transid][0][0]*in[0]+gTrans2dQ[transid][0][1]*in[1];//Cedric 23/02/99
378  out[1] = gTrans2dQ[transid][1][0]*in[0]+gTrans2dQ[transid][1][1]*in[1];//Cedric 23/02/99
379 
380  }
381 
382  TPZTransform<REAL> TPZShapeQuad::ParametricTransform(int trans_id){
383  TPZTransform<REAL> trans(2,2);
384  trans.Mult()(0,0) = gTrans2dQ[trans_id][0][0];
385  trans.Mult()(0,1) = gTrans2dQ[trans_id][0][1];
386  trans.Mult()(1,0) = gTrans2dQ[trans_id][1][0];
387  trans.Mult()(1,1) = gTrans2dQ[trans_id][1][1];
388  return trans;
389  }
390  void TPZShapeQuad::ProjectPoint2dQuadToRib(int rib, TPZVec<REAL> &in, REAL &out) {
391 
392  out = gRibTrans2dQ1d[rib][0]*in[0]+gRibTrans2dQ1d[rib][1]*in[1];
393  }
394 
395  int TPZShapeQuad::GetTransformId2dQ(TPZVec<int64_t> &id) {
396 
397  int id0,id1,minid;
398  id0 = (id[0] < id[1]) ? 0 : 1;
399  id1 = (id[2] < id[3]) ? 2 : 3;
400  minid = (id[id0] < id[id1]) ? id0 : id1;//minid : menor id local
401  id0 = (minid+1)%4;//id anterior local (sentido antihorario)
402  id1 = (minid+3)%4;//id posterior local (sentido horario)
403  minid = id[minid];//minid : menor id global
404 
405  if (id[id0] < id[id1]) {//antihorario
406 
407  if (minid == id[0]) return 0;
408  if (minid == id[1]) return 2;
409  if (minid == id[2]) return 4;
410  if (minid == id[3]) return 6;
411 
412  } else {//horario
413 
414  if (minid == id[0]) return 1;
415  if (minid == id[1]) return 3;
416  if (minid == id[2]) return 5;
417  if (minid == id[3]) return 7;
418  }
419  return 0;
420  }
421 
422 
423 
424  void TPZShapeQuad::TransformDerivativeFromRibToQuad(int rib,int num,TPZFMatrix<REAL> &dphi) {
425 
426  for (int j = 0;j<num;j++) {
427  dphi(1,j) = gRibTrans2dQ1d[rib][1]*dphi(0,j);
428  dphi(0,j) = gRibTrans2dQ1d[rib][0]*dphi(0,j);
429  }
430  }
431 
432  int TPZShapeQuad::NConnectShapeF(int side, int order) {
433  if(side<4) return 1;//0 a 4
434  // int s = side-4;//s = 0 a 14 ou side = 6 a 20
435 #ifdef PZDEBUG
436  {
437 // if(order <1) DebugStop();
438  }
439 #endif
440  if(side<8) return (order-1);//6 a 14
441  if(side==8) {
442  return ((order-1)*(order-1));
443  }
444  PZError << "TPZShapeQuad::NConnectShapeF, bad parameter side " << side << endl;
445  return 0;
446  }
447 
448  int TPZShapeQuad::NShapeF(TPZVec<int> &order) {
449 
450 
451  int in,res=NCornerNodes;
452  for(in=NCornerNodes;in<NSides;in++) {
453 
454  res += NConnectShapeF(in,order[in-NCornerNodes]);
455  }
456  return res;
457  }
458 
459 
460 
461 #ifdef _AUTODIFF
462 
463  void TPZShapeQuad::Shape2dQuadInternal(TPZVec<FADREAL> &x, int order,
464  TPZVec<FADREAL> &phi,int quad_transformation_index) {
465 
466  const int ndim = 3;
467  if(order < 0) return;
468  int ord1 = order+1;
469  int numshape = ord1*ord1;
470  TPZVec<FADREAL> out(2);
471  TransformPoint2dQ(quad_transformation_index,x,out);
472 
473  if(numshape > phi.NElements()/*Rows()*/ || phi[0].size()/*Cols()*/ < ndim) phi.Resize(numshape, FADREAL(ndim, 0.0));
474  //if(dphi.Rows() < 2 || dphi.Cols() < numshape) dphi.Resize(2,numshape);
475  //REAL store1[20],store2[20],store3[20],store4[20];
476  //TPZFMatrix<REAL> phi0(ord1,1,store1,20),phi1(ord1,1,store2,20),dphi0(1,ord1,store3,20),dphi1(1,ord1,store4,20);
477  TPZVec<FADREAL> phi0(20, FADREAL(ndim, 0.0)),
478  phi1(20, FADREAL(ndim, 0.0));
479 
480  TPZShapeLinear::FADfOrthogonal(out[0],ord1,phi0);
481  TPZShapeLinear::FADfOrthogonal(out[1],ord1,phi1);
482  for (int i=0;i<ord1;i++) {
483  for (int j=0;j<ord1;j++) {
484  int index = i*ord1+j;
485  //phi(index,0) = phi0(i,0)* phi1(j,0);
486  phi[index] = phi0[i] * phi1[j];
487  /*dphi(0,index) = dphi0(0,i)* phi1(j,0);
488  dphi(1,index) = phi0(i,0)*dphi1(0,j);*/
489  }
490  }
491  // TransformDerivative2dQ(quad_transformation_index,numshape,phi);
492  }
493 
494  void TPZShapeQuad::TransformPoint2dQ(int transid, TPZVec<FADREAL> &in, TPZVec<FADREAL> &out) {
495 
496  out[0] = gTrans2dQ[transid][0][0]*in[0]+gTrans2dQ[transid][0][1]*in[1];//Cedric 23/02/99
497  out[1] = gTrans2dQ[transid][1][0]*in[0]+gTrans2dQ[transid][1][1]*in[1];//Cedric 23/02/99
498  }
499 
500  /*
501  void TPZShapeQuad::TransformDerivative2dQ(int transid, int num, TPZVec<FADREAL> &in) {
502 
503  for(int i=0;i<num;i++) {
504  double aux[2];
505  aux[0] = in[i].d(0);
506  aux[1] = in[i].d(1);
507  in[i].fastAccessDx(0) = gTrans2dQ[transid][0][0]*aux[0]+gTrans2dQ[transid][1][0]*aux[1];
508  in[i].fastAccessDx(1) = gTrans2dQ[transid][0][1]*aux[0]+gTrans2dQ[transid][1][1]*aux[1];
509  }
510  }
511 
512  void TPZShapeQuad::TransformDerivativeFromRibToQuad(int rib,int num,TPZVec<FADREAL> &phi) {
513 
514  for (int j = 0;j<num;j++) {
515  //dphi(1,j) = gRibTrans2dQ1d[rib][1]*dphi(0,j);
516  //dphi(0,j) = gRibTrans2dQ1d[rib][0]*dphi(0,j);
517  phi[j].fastAccessDx(1) = gRibTrans2dQ1d[rib][1]*phi[j].d(0);
518  phi[j].fastAccessDx(0) = gRibTrans2dQ1d[rib][0]*phi[j].d(0);
519  }
520  }
521  */
522 
523 #endif
524 
525 };
526 
527 
528 
Contains definitions to LOGPZ_DEBUG, LOGPZ_INFO, LOGPZ_WARN, LOGPZ_ERROR and LOGPZ_FATAL, and the implementation of the inline InitializePZLOG(string) function using log4cxx library or not. It must to be called out of "#ifdef LOG4CXX" scope.
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
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
Defines PZError.
Contains TPZShapeLinear class which implements the shape functions of a linear one-dimensional elemen...
int64_t size() const
Returns the number of elements of the vector.
Definition: pzvec.h:196
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
expr_ dx(i) *cos(expr_.val())
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
Free store vector implementation.
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
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...
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
Implements generic class which holds a matrix of objects. Matrix.
Definition: pzshtmat.h:18
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
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
int64_t NElements() const
Returns the number of elements of the vector.
Definition: pzvec.h:190
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
Contains TPZShapeQuad class which implements the shape functions of a quadrilateral 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