NeoPZ
tpztriangle.cpp
Go to the documentation of this file.
1 
6 #include "tpztriangle.h"
7 #include "pzquad.h"
8 
9 //#include "pzshapetriang.h"
10 
11 #ifdef _AUTODIFF
12 #include "fad.h"
13 #endif
14 
15 #include "pzlog.h"
16 #include <cmath>
17 
18 #ifdef LOG4CXX
19 static LoggerPtr logger(Logger::getLogger("pz.topology.pztriangle"));
20 #endif
21 
22 using namespace std;
23 
24 namespace pztopology {
25 
26  template<class T>
27  inline void TPZTriangle::TShape(const TPZVec<T> &loc,TPZFMatrix<T> &phi,TPZFMatrix<T> &dphi) {
28  T qsi = loc[0], eta = loc[1];
29  phi(0,0) = 1.0-qsi-eta;
30  phi(1,0) = qsi;
31  phi(2,0) = eta;
32  qsi *= (T)0.; // Keep this for FAD calculation
33  dphi(0,0) = dphi(1,0) = -1.0 + qsi;
34  dphi(0,1) = dphi(1,2) = 1.0 + qsi;
35  dphi(1,1) = dphi(0,2) = qsi;
36  }
37 
38  int TPZTriangle::SideNodes[3][2] = { {0,1},{1,2},{2,0} };
39  int TPZTriangle::FaceNodes[1][3] = { {0,1,2} };
40 
41 
43  REAL TPZTriangle::gTrans2dT[6][2][2] = {//s* , t*
44  { { 1., 0.},{ 0., 1.} },
45  { { 0., 1.},{ 1., 0.} },
46  { { 0., 1.},{-1.,-1.} },//s* = t t* = -s-t-1 , etc
47  { {-1.,-1.},{ 0., 1.} },
48  { {-1.,-1.},{ 1., 0.} },
49  { { 1., 0.},{-1.,-1.} }
50  };
51 
52  template<class T>
53  void TPZTriangle::BlendFactorForSide(const int &side, const TPZVec<T> &xi, T &blendFactor,
54  TPZVec<T> &blendFactorDxi){
55  const REAL tol = pztopology::GetTolerance();
56  blendFactorDxi.Resize(TPZTriangle::Dimension, (T) 0);
57 #ifdef PZDEBUG
58  std::ostringstream sout;
59  if(side < NCornerNodes || side >= NSides){
60  sout<<"The side\t"<<side<<"is invalid. Aborting..."<<std::endl;
61  }
62 
64  sout<<"The method BlendFactorForSide expects the point xi to correspond to coordinates of a point";
65  sout<<" inside the parametric domain. Aborting...";
66  }
67 
68  if(!CheckProjectionForSingularity(side,xi)){
69  sout<<"The projection of xi "<<xi[0]<<" "<<xi[1]<<" to side "<<side<<" is singular."<<std::endl;
70  sout<<"This should have been caught by MapToSide method. Aborting..."<<std::endl;
71  }
72  if(!sout.str().empty()){
73  PZError<<std::endl<<sout.str()<<std::endl;
74 #ifdef LOG4CXX
75  LOGPZ_FATAL(logger,sout.str().c_str());
76 #endif
77  DebugStop();
78  }
79 #endif
80  //if the point is singular, the blend factor and its derivatives should be zero
81  if(!CheckProjectionForSingularity(side,xi)){
82  std::cout<<"Side projection is not regular and it should have been checked earlier. Aborting.."<<std::endl;
83  DebugStop();
84  blendFactor = 0;
85  for(int i = 0; i < blendFactorDxi.size(); i++) blendFactorDxi[i] = 0;
86  return;
87  }
88 
89  TPZFNMatrix<4,T> phi(NCornerNodes,1);
90  TPZFNMatrix<8,T> dphi(Dimension,NCornerNodes);
91  TPZTriangle::TShape(xi,phi,dphi);
92  int i = -1;
93  switch(side){
94  case 0:
95  case 1:
96  case 2:
97  blendFactor = 0;
98  return;
99  case 3:
100  i = 0;
101  break;
102  case 4:
103  i = 1;
104  break;
105  case 5:
106  i = 2;
107  break;
108  case 6:
109  blendFactor = 1;
110  return;
111  }
112  blendFactor = phi(i,0) + phi((i+1)%NCornerNodes,0);
113  blendFactor *= blendFactor;
114  blendFactorDxi[0] = 2 * ( phi(i,0) + phi((i+1)%NCornerNodes,0) ) * ( dphi(0,i) + dphi(0,(i+1)%NCornerNodes) );
115  blendFactorDxi[1] = 2 * ( phi(i,0) + phi((i+1)%NCornerNodes,0) ) * ( dphi(1,i) + dphi(1,(i+1)%NCornerNodes) );
116 
117  }
118 
119  static int sidedimension[7] = {0,0,0,1,1,1,2};
120 
121  static int nhighdimsides[7] = {3,3,3,1,1,1,0};
122 
123  static int highsides[7][3] = {
124  {3,5,6},
125  {3,4,6},
126  {4,5,6},
127  {6},
128  {6},
129  {6},
130  {-999}
131  };
132 
133  static REAL sidetosidetransforms[7][3][4][3] = {
134  {
135  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
136  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
137  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,-99}}
138  },
139  {
140  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
141  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
142  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,-99}}
143  },
144  {
145  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
146  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
147  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,-99}}
148  },
149  {
150  {{0.5,0,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0,-99}}
151  },
152  {
153  {{-0.5,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,-99}}
154  },
155  {
156  {{0,-0.5,-99},{-99,-99,-99},{-99,-99,-99},{0,0.5,-99}}
157  },
158  {
159  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-99,-99,-99}}
160  }
161  };
162 
163  static REAL MidSideNode[7][3] = {
164  /*00*/{.0,0.},/*01*/{1.0,.0},/*02*/{0.,1.0},
165  /*03*/{.5,0.},/*04*/{0.5,.5},/*05*/{0.,0.5},
166  /*06*/{ 1./3.,1./3.} };
167 
168  static int nsidenodes[7] = {1,1,1,2,2,2,3};
169 
170  static REAL bTriang[14][2] =
171  {
172  {0,-1},//0
173  {1,-1},
174  {0,-1},
175  {1,0},
176  {0,1},
177  {0.5,0.5},//5
178  {-1,1},
179  {-1,0},
180  {-1,0},//8
181  {1,0},
182  {-0.5,0.5},
183  {0,-1},
184  {1,0},
185  {0,1}
186 
187  };
188 
189  static REAL tTriang[14][2] =
190  {
191  {-1,0},
192  {-1,0},
193  {-1,0},
194  {1,-1},
195  {1,-1},
196  {1,-1},
197  {0,1},
198  {0,1},
199  {0,1},//
200  {0,-1},
201  {1,1},
202  {-1,0},
203  {0,-1},
204  {1,0}
205  };
206 
207  static int vectorsideorder [14] = {0,1,3,1,2,4,2,0,5,3,4,5,6,6};
208 
209  //static int bilinearounao [14] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0};//Pk Pk-1
210  static int bilinearounao [14] = {0,0,0,0,0,0,0,0,0,1,1,1,1,1};//P*k Pk
211 
212  static int direcaoksioueta [14] = {0,0,0,0,0,0,0,0,0,0,0,0,0,1};
213 
214  int TPZTriangle::fPermutations [6][7] =
215  {
216  {0,1,2,3,4,5,6}, // id 0
217  {0,2,1,5,4,3,6}, // id 1
218  {1,2,0,4,5,3,6}, // id 2
219  {1,0,2,3,5,4,6}, // id 3
220  {2,0,1,5,3,4,6}, // id 4
221  {2,1,0,4,3,5,6} // id 5
222  };
223 
224  REAL TPZTriangle::fTangentVectors [12][2] =
225  {
226  {2,0}, // id 0
227  {0,2}, // id 0
228  {0,2}, // id 1
229  {2,0}, // id 1
230  {0,2}, // id 2
231  {-2,-2}, //id 2
232  {-2,-2},// id 3
233  {0,2},// id 3
234  {-2,-2},// id 4
235  {2,0}, //id 4
236  {2,0}, //id 5
237  {-2,-2}, //id 5
238  };
239 
240  int TPZTriangle::NBilinearSides()
241  {
242  DebugStop();
243  return 5;
244  }
245 
246 
247  int TPZTriangle::NSideNodes(int side)
248  {
249  return nsidenodes[side];
250  }
251 
252  bool TPZTriangle::IsInParametricDomain(const TPZVec<REAL> &pt, REAL tol){
253  const REAL qsi = pt[0];
254  const REAL eta = pt[1];
255  if( ( qsi <= 1. + tol ) && ( qsi >= 0. - tol ) &&
256  ( eta <= 1. + tol ) && ( eta >= 0. - tol ) &&
257  ( eta <= 1. - qsi + tol ) ){
258  return true;
259  }
260  else{
261  return false;
262  }
263  }//method
264 
265  template<class T>
266  bool TPZTriangle::CheckProjectionForSingularity(const int &side, const TPZVec<T> &xiInterior) {
267 
268  double zero = pztopology::GetTolerance();
269  T qsi = xiInterior[0]; T eta = xiInterior[1];
270 
271  switch(side)
272  {
273  case 0:
274  case 1:
275  case 2:
276  return true;
277  case 3:
278  if(fabs((T)(eta - 1.)) < zero) return false;
279  case 4:
280  if((T)(qsi+eta) < (T)zero) return false;
281  case 5:
282  if(fabs((T)(qsi - 1.)) < zero) return false;
283  case 6: return true;
284  }
285  if(side > 6)
286  {
287  cout << "Cant compute CheckProjectionForSingularity method in TPZTriangle class!\nParameter (SIDE) must be 3, 4 or 5!\nMethod Aborted!\n";
288  DebugStop();
289  }
290  return true;
291  }
292 
293  template<class T>
294  void TPZTriangle::MapToSide(int side, TPZVec<T> &InternalPar, TPZVec<T> &SidePar, TPZFMatrix<T> &JacToSide) {
295 
296  T qsi = InternalPar[0]; T eta = InternalPar[1];
297  SidePar.Resize(1); JacToSide.Resize(1,2);
298 
299  if(!CheckProjectionForSingularity(side,InternalPar)){
300  std::cout<<"Side projection is not regular and it should have been checked earlier. Aborting.."<<std::endl;
301  DebugStop();
302  }
303 
304  switch(side)
305  {
306  case 0:
307  case 1:
308  case 2:
309  SidePar.Resize(0); JacToSide.Resize(0,0);
310  break;
311  case 3:
312  SidePar[0] = 2.*qsi/(1.-eta) - 1.;
313  JacToSide(0,0) = 2./(1.-eta); JacToSide(0,1) = 2.*qsi/((1.-eta)*(1.-eta));
314  break;
315 
316  case 4:
317  SidePar[0] = 1. - 2.*qsi/(qsi + eta);
318  JacToSide(0,0) = -2.*eta/((qsi+eta)*(qsi+eta)); JacToSide(0,1) = 2.*qsi/((qsi+eta)*(qsi+eta));
319  break;
320 
321  case 5:
322  SidePar[0] = 1. - 2.*eta/(1.-qsi);
323  JacToSide(0,0) = -2.*eta/((1.-qsi)*(1.-qsi)); JacToSide(0,1) = -2./(1.-qsi);
324  break;
325  case 6:
326  SidePar = InternalPar;
327  JacToSide.Resize(2, 2);
328  JacToSide.Identity();
329  }
330  }
331 
332  void TPZTriangle::ParametricDomainNodeCoord(int node, TPZVec<REAL> &nodeCoord)
333  {
334  if(node > NCornerNodes)
335  {
336  DebugStop();
337  }
338  nodeCoord.Resize(Dimension, 0.);
339  switch (node) {
340  case (0):
341  {
342  nodeCoord[0] = 0.;
343  nodeCoord[1] = 0.;
344  break;
345  }
346  case (1):
347  {
348  nodeCoord[0] = 1.;
349  nodeCoord[1] = 0.;
350  break;
351  }
352  case (2):
353  {
354  nodeCoord[0] = 0.;
355  nodeCoord[1] = 1.;
356  break;
357  }
358  default:
359  {
360  DebugStop();
361  break;
362  }
363  }
364  }
365 
367  void TPZTriangle::RandomPoint(TPZVec<REAL> &pt)
368  {
369  REAL val = (REAL) rand() / (RAND_MAX);
370  pt[0] = val;
371  val = (1.-pt[0]) * (REAL) rand() / (RAND_MAX);
372  pt[1] = val;
373  }
374 
375 
376  TPZTransform<> TPZTriangle::SideToSideTransform(int sidefrom, int sideto)
377  {
378  if(sidefrom <0 || sidefrom >= NSides || sideto <0 || sideto >= NSides) {
379  PZError << "TPZTriangle::HigherDimensionSides sidefrom "<< sidefrom <<
380  ' ' << sideto << endl;
381  return TPZTransform<>(0);
382  }
383  if(sidefrom == sideto) {
384  return TPZTransform<>(sidedimension[sidefrom]);
385  }
386  if(sidefrom == NSides-1) {
387  return TransformElementToSide(sideto);
388  }
389  if (sideto == NSides-1) {
390  return TransformSideToElement(sidefrom);
391  }
392 
393  int nhigh = nhighdimsides[sidefrom];
394  int is;
395  for(is=0; is<nhigh; is++) {
396  if(highsides[sidefrom][is] == sideto) {
397  int dfr = sidedimension[sidefrom];
398  int dto = sidedimension[sideto];
399  TPZTransform<> trans(dto,dfr);
400  int i,j;
401  for(i=0; i<dto; i++) {
402  for(j=0; j<dfr; j++) {
403  trans.Mult()(i,j) = sidetosidetransforms[sidefrom][is][j][i];
404  }
405  trans.Sum()(i,0) = sidetosidetransforms[sidefrom][is][3][i];
406  }
407  return trans;
408  }
409  }
410  PZError << "TPZTriangle::SideToSideTransform highside not found sidefrom "
411  << sidefrom << ' ' << sideto << endl;
412  return TPZTransform<>(0);
413  }
414 
415  int TPZTriangle::SideNodeLocId(int side, int node)
416  {
417  if(side<3 && node == 0) return side;
418  if(side>=3 && side<6 && node <2) return (side-3+node) %3;
419  if(side==6 && node <3) return node;
420  PZError << "TPZTriangle::SideNodeLocId inconsistent side or node " << side
421  << ' ' << node << endl;
422  return -1;
423  }
424 
425  void TPZTriangle::LowerDimensionSides(int side,TPZStack<int> &smallsides)
426  {
427  smallsides.Resize(0);
428  int nsidecon = NContainedSides(side);
429  int is;
430  for(is=0; is<nsidecon-1; is++)
431  smallsides.Push(ContainedSideLocId(side,is));
432  }
433 
434  void TPZTriangle::LowerDimensionSides(int side,TPZStack<int> &smallsides, int DimTarget)
435  {
436  smallsides.Resize(0);
437  int nsidecon = NContainedSides(side);
438  for(int is = 0; is < nsidecon - 1; is++) {
439  if (SideDimension(ContainedSideLocId(side,is)) == DimTarget) smallsides.Push(ContainedSideLocId(side,is));
440  }
441  }
442 
443  void TPZTriangle::HigherDimensionSides(int side, TPZStack<int> &high)
444  {
445  if(side <0 || side >= NSides) {
446  PZError << "TPZTriangle::HigherDimensionSides side "<< side << endl;
447  }
448  int is;
449  for(is=0; is<nhighdimsides[side]; is++) high.Push(highsides[side][is]);
450 
451  }
452 
453  //Tentando criar o metodo
454  int TPZTriangle::NumSides(int dimension) {
455  if(dimension<0 || dimension> 2) {
456  PZError << "TPZTriangle::NumSides. Bad parameter i.\n";
457  return 0;
458  }
459  if(dimension==0) return 3;
460  if(dimension==1) return 3;
461  if(dimension==2) return 1;
462  return -1;
463  }
464  int TPZTriangle::SideDimension(int side) {
465  if(side<0 || side >= NSides) {
466  PZError << "TPZTriangle::SideDimension side " << side << endl;
467  return -1;
468  }
469  return sidedimension[side];
470  }
471 
472  void TPZTriangle::CenterPoint(int side, TPZVec<REAL> &center) {
473  if (center.size()!=Dimension) {
474  DebugStop();
475  }
476  int i;
477  for(i=0; i<Dimension; i++) {
478  center[i] = MidSideNode[side][i];
479  }
480  }
481 
482 
483  TPZTransform<> TPZTriangle::TransformElementToSide(int side) {
484  if(side<0 || side>6){
485  PZError << "TPZTriangle::TransformElementToSide called with side error\n";
486  return TPZTransform<>(0,0);
487  }
488 
489  TPZTransform<> t(sidedimension[side],2);//t(dimto,2)
490  t.Mult().Zero();
491  t.Sum().Zero();
492 
493  switch(side){
494  case 0:
495  case 1:
496  case 2:
497  return t;
498  case 3:
499  t.Mult()(0,0) = 2.0;//par. var.
500  t.Mult()(0,1) = 1.0;//par. var.
501  t.Sum()(0,0) = -1.0;
502  return t;
503  case 4:
504  t.Mult()(0,0) = -1.0;
505  t.Mult()(0,1) = 1.0;
506  return t;
507  case 5:
508  t.Mult()(0,0) = -1.0;
509  t.Mult()(0,1) = -2.0;
510  t.Sum()(0,0) = 1.0;
511  return t;
512  case 6:
513  t.Mult()(0,0) = 1.0;
514  t.Mult()(1,1) = 1.0;
515 // t.Sum()(0,0) = -1;
516 // t.Sum()(1,0) = -1;
517  return t;
518  }
519  return TPZTransform<>(0,0);
520  }
521 
522 
523  TPZTransform<> TPZTriangle::TransformSideToElement(int side){
524 
525  if(side<0 || side>6){
526  PZError << "TPZTriangle::TransformSideToElement side out range\n";
527  return TPZTransform<>(0,0);
528  }
529  TPZTransform<> t(2,sidedimension[side]);
530  t.Mult().Zero();
531  t.Sum().Zero();
532 
533  switch(side){
534  case 0:
535  return t;
536  case 1:
537  t.Sum()(0,0) = 1.0;
538  return t;
539  case 2:
540  t.Sum()(1,0) = 1.0;
541  return t;
542  case 3:
543  t.Mult()(0,0) = 0.5;
544  t.Sum()(0,0) = 0.5;
545  return t;
546  case 4:
547  t.Mult()(0,0) = -0.5;
548  t.Mult()(1,0) = 0.5;
549  t.Sum() (0,0) = 0.5;
550  t.Sum() (1,0) = 0.5;
551  return t;
552  case 5:
553  t.Mult()(1,0) = -0.5;
554  t.Sum() (1,0) = 0.5;
555  return t;
556  case 6:
557  t.Mult()(0,0) = 1.0;
558  t.Mult()(1,1) = 1.0;
559  return t;
560  }
561  return TPZTransform<>(0,0);
562  }
563 
564  TPZIntPoints * TPZTriangle::CreateSideIntegrationRule(int side, int order){
565  if(side < 0 || side>6) {
566  PZError << "TPZTriangle::CreateSideIntegrationRule wrong side = " << side << ".\n";
567  return 0;
568  }
569  if(side<3) return new TPZInt1Point(order); // sides 0 to 2 are vertices (corners)
570  if(side<6) return new TPZInt1d(order); // sides 3 to 5 are lines
571  if(side==6)return new IntruleType(order); // integration of the element
572  return 0;
573  }
574 
575 
576  MElementType TPZTriangle::Type()
577  {
578  return ETriangle;
579  }
580 
581  MElementType TPZTriangle::Type(int side)
582  {
583  switch(side) {
584  case 0:
585  case 1:
586  case 2:
587  return EPoint;
588  case 3:
589  case 4:
590  case 5:
591  return EOned;
592  case 6:
593  return ETriangle;
594  default:
595  return ENoType;
596  }
597  }
598 
599 
600  int TPZTriangle::NumSides() {
601  return NSides;
602  }
603 
604 
605 
606  int TPZTriangle::NContainedSides(int side) {
607  if(side<0 || side>6) {
608  PZError << "TPZShapeTriang::NContainedSides. Bad parameter side = " << side << ".\n";
609  return 0;
610  }
611  if(side<3) return 1;
612  if(side<6) return 3;
613  return 7;
614  }
615 
617  // side eh o lado do elemento, c eh o noh do lado
618  int TPZTriangle::ContainedSideLocId(int side, int c) {
619  switch(side) {
620  case 0:
621  case 1:
622  case 2:
623  return side;
624  case 3:
625  case 4:
626  case 5:
627  if(!c) return side-3;
628  if(c==1) return (side-2)%3;
629  if(c==2) return side;
630  case 6:
631  return c;
632  default:
633  PZError << "TPZShapeTriang::ContainedSideLocId, connect = " << c << ".\n";
634  return -1;
635  }
636  }
637 
644  int TPZTriangle::GetTransformId(TPZVec<int64_t> &id)
645  {
646  int id0, id1, minid;
647  id0 = (id[0] < id[1]) ? 0 : 1;
648  minid = (id[2] < id[id0]) ? 2 : id0;
649  id0 = (minid + 1) % 3;
650  id1 = (minid + 2) % 3;
651 
652  if (id[id0] < id[id1]) {//antihorario
653 
654  if (minid == 0) return 0;
655  if (minid == 1) return 2;
656  if (minid == 2) return 4;
657 
658  }
659  else {//horario
660 
661  if (minid == 0) return 1;
662  if (minid == 1) return 3;
663  if (minid == 2) return 5;
664  }
665  return 0;
666 
667  }
674  int TPZTriangle::GetTransformId(int side, TPZVec<int64_t> &id)
675  {
676  switch (side) {
677  case 0:
678  case 1:
679  case 2:
680  return 0;
681  break;
682  case 3:
683  case 4:
684  case 5:
685  {
686  int in1 = ContainedSideLocId(side,0);
687  int in2 = ContainedSideLocId(side,1);
688  return id[in1]<id[in2] ? 0 : 1;
689  }
690  break;
691  case 6:
692  {
693  return GetTransformId(id);
694  }
695  break;
696  default:
697  break;
698  }
699  LOGPZ_ERROR(logger,"Wrong side parameter")
700  return -1;
701  }
702 
706 void TPZTriangle::GetHDivGatherPermute(int transformid, TPZVec<int> &permute)
707 {
708 #ifdef PZDEBUG
709  if (permute.size() != 7 || transformid >= 6 || transformid < 0) {
710  DebugStop();
711  }
712 #endif
713  for (int i=0; i<7; i++) {
714  permute[i] = fPermutations[transformid][i];
715  }
716  return;
717  int dir = 1;
718  if (transformid%2 ==1) dir = -1;
719  int runsmall = transformid/2;
720  int runlarge = runsmall+3;
721  if (dir == -1) {
722  runlarge += dir;
723  if (runlarge < 3) {
724  runlarge += 3;
725  }
726  }
727  permute[6] = 6;
728  for (int is=0; is<3; is++) {
729  permute[is] = runsmall;
730  permute[is+3] = runlarge;
731  runsmall += dir;
732  runlarge += dir;
733  if (dir == 1 && runsmall > 2) {
734  runsmall -= 3;
735  runlarge -= 3;
736  }
737  else if(dir == -1 && runsmall < 0)
738  {
739  runsmall += 3;
740  }
741  if (dir == -1 && runlarge < 3) {
742  runlarge += 3;
743  }
744  }
745 }
746 
747 
755  void TPZTriangle::GetSideHDivPermutation(int transformationid, TPZVec<int> &permgather)
756  {
757  permgather.Resize(7);
758 
759  for (int i=0; i<7; i++)
760  {
761  permgather[i] = fPermutations[transformationid][i];
762  }
763 
764  /*
765  switch (transformationid)
766  {
767  case 0:
768  {
769  for (int i=0; i<7; i++)
770  {
771  permgather[i] = fPermutations[0][i];
772  }
773  }
774  break;
775  case 1:
776  {
777  for (int i=0; i<7; i++)
778  {
779  permgather[i] = fPermutations[1][i];
780  }
781  }
782  break;
783  case 2:
784  {
785  for (int i=0; i<7; i++)
786  {
787  permgather[i] = fPermutations[2][i];
788  }
789  }
790  break;
791  case 3:
792  {
793  for (int i=0; i<7; i++)
794  {
795  permgather[i] = fPermutations[3][i];
796  }
797  }
798  break;
799  case 4:
800  {
801  for (int i=0; i<7; i++)
802  {
803  permgather[i] = fPermutations[4][i];
804  }
805  }
806  break;
807  case 5:
808  {
809  for (int i=0; i<7; i++)
810  {
811  permgather[i] = fPermutations[5][i];
812  }
813  }
814  break;
815  default:
816  DebugStop();
817  break;
818  }
819  */
820 
821  /*switch (side) {
822  case 0:
823  case 1:
824  case 2:
825  permgather[0] = 0;
826  break;
827  case 3:
828  case 4:
829  case 5:
830  {
831  int in1 = ContainedSideLocId(side,0);
832  int in2 = ContainedSideLocId(side,1);
833  if(in1<in2)
834  {
835  permgather[0] = 0;
836  permgather[1] = 1;
837  permgather[2] = 2;
838  }
839  else
840  {
841  permgather[0] = 1;
842  permgather[1] = 0;
843  permgather[2] = 2;
844  }
845  }
846  break;
847  case 6:
848  {
849  int i;
850  int tid = pzshape::TPZShapeTriang::GetTransformId2dT(id);
851  if(tid%2 == 0)
852  {
853  switch (tid)
854  {
855  case 0:
856  for(i=0; i<7; i++) permgather[i] = i;
857  break;
858  case 2:
859  for(i=0; i<3; i++) permgather[i] = (i+1)%3;
860  for(i=4; i<6; i++) permgather[i] = 3+(i+1)%3;
861  permgather[6] = 6;
862  break;
863  case 4:
864  for(i=0; i<3; i++) permgather[i] = (i+2)%3;
865  for(i=4; i<6; i++) permgather[i] = 3+(i+2)%3;
866  permgather[6] = 6;
867  break;
868  }
869  }
870  else
871  {
872  TPZManVector<int,3> invid(3);
873  invid[0] = 0;
874  invid[1] = 2;
875  invid[2] = 1;
876  switch (tid) {
877  case 1:
878  for(i=0; i<3; i++) permgather[i] = invid[i];
879  for(i=4; i<6; i++) permgather[i] = 3+invid[(i+0)%3];
880  permgather[6] = 6;
881  break;
882  case 3:
883  for(i=0; i<3; i++) permgather[i] = invid[(i+1)%3];
884  for(i=4; i<6; i++) permgather[i] = 3+invid[(i+1)%3];
885  permgather[6] = 6;
886  break;
887  case 5:
888  for(i=0; i<3; i++) permgather[i] = invid[(i+2)%3];
889  for(i=4; i<6; i++) permgather[i] = 3+invid[(i+2)%3];
890  permgather[6] = 6;
891  break;
892  default:
893  break;
894  }
895  }
896  }
897  break;
898  default:
899  break;
900  }
901  LOGPZ_ERROR(logger,"Wrong side parameter")*/
902  }
903 
904 
905 
906 // void computedirectionst(int inicio, int fim, TPZFMatrix<REAL> &bvec, TPZFMatrix<REAL> &t1vec,
907 // TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions);
908 // void computedirectionst(int inicio, int fim, TPZFMatrix<REAL> &bvec, TPZFMatrix<REAL> &t1vec,
909 // TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions)
910 // {
911 // REAL detgrad = 0.0;
912 // TPZVec<REAL> u(3);
913 // TPZVec<REAL> v(3);
914 // TPZVec<REAL> uxv(3);// result
915 //
916 // for (int ilin=0; ilin<3; ilin++)
917 // {
918 // u[ilin] = gradx(ilin, 0);
919 // v[ilin] = gradx(ilin, 1);
920 // }
921 //
922 // //TPZNumeric::ProdVetorial(u,v,uxv);
923 // uxv[0] = u[1]*v[2]-u[2]*v[1];
924 // uxv[1] = -(u[0]*v[2]-v[0]*u[2]);
925 // uxv[2] = u[0]*v[1]-v[0]*u[1];
926 //
927 // for (int pos=0; pos<3; pos++)
928 // {
929 // detgrad += uxv[pos]*uxv[pos];
930 // }
931 // detgrad = sqrt(fabs(detgrad));
932 //
933 // int cont = 0;
934 //
935 // for (int ivet=inicio; ivet<=fim; ivet++)
936 // {
937 // TPZFMatrix<REAL> Wvec(3,1);
938 // TPZVec<REAL> uxvtmp(3);
939 // REAL acumng = 0.0;
940 // // calc do g gradx*t
941 // TPZManVector<REAL,3> gvec(3,0.),Vvec(3,0.);
942 // REAL gvecnorm;
943 // for (int il=0; il<3; il++)
944 // {
945 // for (int i = 0 ; i<2; i++)
946 // {
947 // gvec[il] += gradx(il,i) * t1vec(i,ivet);
948 // Vvec[il] += gradx(il,i) * bvec(i,ivet);
949 // }
950 // u[il] = gvec[il];
951 // acumng += gvec[il]*gvec[il];
952 // }
953 // gvecnorm = sqrt(acumng);
954 //
955 // for (int il=0; il<3; il++)
956 // {
957 // Wvec(il,0) = Vvec[il]*gvecnorm/detgrad;
958 // directions(il,cont) = Wvec(il,0);
959 // }
960 // cont++;
961 // }
962 //
963 //
964 // }
965 
966 // void TPZTriangle::ComputeDirections(int side, TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions, TPZVec<int> &sidevectors)
967 // {
968 //
969 // if(gradx.Cols()!=2)
970 // {
971 // DebugStop();
972 // }
973 //
974 // TPZFMatrix<REAL> bvec(2,14);
975 // TPZFMatrix<REAL> t1vec(2,14);
976 //
977 // bvec.Redim(2, 14);
978 // t1vec.Redim(2, 14);
979 //
980 // for (int lin = 0; lin<14; lin++)
981 // {
982 // for(int col = 0;col<2;col++)
983 // {
984 // bvec.PutVal(col, lin, bTriang[lin][col]);
985 // t1vec.PutVal(col, lin, tTriang[lin][col]);
986 // }
987 // }
988 // // calcula os vetores
989 // //int numvec = bvec.Cols();
990 //
991 // switch (side)
992 // {
993 // case 0:
994 // {
995 //
996 // }
997 // break;
998 // case 1:
999 // {
1000 //
1001 // }
1002 // break;
1003 // case 2:
1004 // {
1005 //
1006 // }
1007 // break;
1008 // case 3:
1009 // {
1010 // directions.Redim(3, 3);
1011 // sidevectors.Resize(3);
1012 // int inumvec = 0, fnumvec = 2;
1013 // computedirectionst(inumvec, fnumvec, bvec, t1vec, gradx, directions);
1014 // for (int ip = 0; ip < 3; ip++) {
1015 // sidevectors[ip] = vectorsideorder[ip+inumvec];
1016 // }
1017 //
1018 // }
1019 // break;
1020 // case 4:
1021 // {
1022 // directions.Redim(3, 3);
1023 // sidevectors.Resize(3);
1024 // int inumvec = 3, fnumvec = 5;
1025 // computedirectionst(inumvec, fnumvec, bvec, t1vec, gradx, directions);
1026 // for (int ip = 0; ip < 3; ip++) {
1027 // sidevectors[ip] = vectorsideorder[ip+inumvec];
1028 // }
1029 // }
1030 // break;
1031 // case 5:
1032 // {
1033 // directions.Redim(3, 3);
1034 // sidevectors.Resize(3);
1035 // int inumvec = 6, fnumvec = 8;
1036 // computedirectionst(inumvec, fnumvec, bvec, t1vec, gradx, directions);
1037 // for (int ip = 0; ip < 3; ip++) {
1038 // sidevectors[ip] = vectorsideorder[ip+inumvec];
1039 // }
1040 // }
1041 // break;
1042 // case 6:
1043 // {
1044 // directions.Redim(3, 5);
1045 // sidevectors.Resize(5);
1046 // int inumvec = 9, fnumvec = 13;
1047 // computedirectionst(inumvec, fnumvec, bvec, t1vec, gradx, directions);
1048 // for (int ip = 0; ip < 5; ip++) {
1049 // sidevectors[ip] = vectorsideorder[ip+inumvec];
1050 // }
1051 // }
1052 // break;
1053 // default:
1054 // DebugStop();
1055 // break;
1056 // }
1057 //
1058 // }
1059 
1060  template <class TVar>
1061  void TPZTriangle::ComputeHDivDirections(TPZFMatrix<TVar> &gradx, TPZFMatrix<TVar> &directions)
1062  {
1063  TVar detjac = TPZAxesTools<TVar>::ComputeDetjac(gradx);
1064 
1065  TPZManVector<TVar, 3> v1(3),v2(3), vdiag(3);
1066  for (int i=0; i<3; i++) {
1067  v1[i] = gradx(i,0);
1068  v2[i] = gradx(i,1);
1069  vdiag[i] = (gradx(i,0)-gradx(i,1));
1070  }
1071 
1072  TVar Nv1 = TPZNumeric::Norm(v1);
1073  TVar Nv2 = TPZNumeric::Norm(v2);
1074  TVar Nvdiag = TPZNumeric::Norm(vdiag);
1075 
1076 
1082  TPZManVector<TVar,3> NormalScales(3,1.);
1083 
1084 
1085  {
1086  NormalScales[0] = 2./Nv1;
1087  NormalScales[1] = 2./Nvdiag;
1088  NormalScales[2] = 2./Nv2;
1089  }
1090 
1091  for (int i=0; i<3; i++) {
1092  v1[i] /= detjac;
1093  v2[i] /= detjac;
1094 
1095  }
1096 
1097  for (int i=0; i<3; i++)
1098  {
1099  directions(i,0) = -v2[i]*Nv1*NormalScales[0];
1100  directions(i,1) = (v1[i]-v2[i])*Nv1*NormalScales[0];
1101  directions(i,2) = (directions(i,0)+directions(i,1))/2.;
1102  directions(i,3) = v1[i]*Nvdiag*NormalScales[1];
1103  directions(i,4) = v2[i]*Nvdiag*NormalScales[1];
1104  directions(i,5) = (directions(i,3)+directions(i,4))/2.;
1105  directions(i,6) = (v2[i]-v1[i])*Nv2*NormalScales[2];
1106  directions(i,7) = -v1[i]*Nv2*NormalScales[2];
1107  directions(i,8) = (directions(i,6)+directions(i,7))/2.;
1108 
1109  directions(i,9) = v1[i]*Nv1*NormalScales[0];
1110  directions(i,10) = (v2[i]-v1[i])*Nvdiag*NormalScales[1];
1111  directions(i,11) = -v2[i]*Nv2*NormalScales[2];
1112 
1113  directions(i,12) = v1[i]*Nv2*NormalScales[0];
1114  directions(i,13) = v2[i]*Nv1*NormalScales[1];
1115  }
1116  }
1117 
1118  template <class TVar>
1119  void TPZTriangle::ComputeHCurlDirections(TPZFMatrix<TVar> &gradx, TPZFMatrix<TVar> &directions, const TPZVec<int> &transformationIds)
1120  {
1121  const auto dim = gradx.Rows();
1122  TPZManVector<TVar, 3> v1(dim),v2(dim);
1123  for (int i=0; i<dim; i++) {
1124  v1[i] = gradx(i,0);
1125  v2[i] = gradx(i,1);
1126  }
1127  constexpr auto nEdges{3};
1128  constexpr REAL faceArea{TPZTriangle::RefElVolume()};
1129  constexpr REAL edgeLength[nEdges]{1,M_SQRT2,1};
1130  const int facePermute = transformationIds[nEdges];
1131  TPZManVector<REAL,nEdges> edgeSign(nEdges,0);
1132  for(auto iEdge = 0; iEdge < nEdges; iEdge++){
1133  edgeSign[iEdge] = transformationIds[iEdge] == 0 ? 1 : -1;
1134  }
1135  for (int i=0; i<dim; i++)
1136  {
1137  //v^{e,a} constant vector fields associated with edge e and vertex a
1138  //they are defined in such a way that v^{e,a} is normal to the edge \hat{e}
1139  //adjacent to edge e by the vertex a. the tangential component is set to be 1 /edgeLength[e]
1140  directions(i,0) = (v1[i]) * edgeSign[0] / edgeLength[0];
1141  directions(i,1) = (v1[i]+v2[i]) * edgeSign[0] / edgeLength[0];
1142  directions(i,2) = (v2[i]*M_SQRT2) * edgeSign[1] / edgeLength[1];
1143  directions(i,3) = (-v1[i]*M_SQRT2) * edgeSign[1] / edgeLength[1];
1144  directions(i,4) = (v1[i] + v2[i]) * -1 * edgeSign[2] / edgeLength[2];
1145  directions(i,5) = (-v2[i]) * edgeSign[2] / edgeLength[2];
1146 
1147  //v^{e,T} constant vector fields associated with edge e and aligned with it
1148  directions(i,6) = (v1[i]) * edgeSign[0] / edgeLength[0];
1149  directions(i,7) = ( (v2[i] - v1[i]) * M_SQRT1_2) * edgeSign[1] / edgeLength[1];
1150  directions(i,8) = (-v2[i]) * edgeSign[2] / edgeLength[2];
1151 
1152  //v^{F,e} constant vector fields associated with face F and edge e
1153  //they are defined in such a way that v^{F,e} is normal to the face \hat{F}
1154  //adjacent to face F by edge e
1155  directions(i,9) = v2[i] * edgeSign[0] /faceArea;//* edgeLength[0];
1156  directions(i,10) = ((-v2[i]-v1[i]) * M_SQRT1_2) * edgeSign[1] /faceArea;//* edgeLength[1];
1157  directions(i,11) = v1[i] * edgeSign[2] /faceArea;//* edgeLength[0];
1158 
1159  //v^{F,T} orthonormal vectors associated with face F and tangent to it.
1160  directions(i,12) = fTangentVectors[2*facePermute][i];
1161  directions(i,13) = fTangentVectors[2*facePermute + 1][i];
1162  }
1163  }
1164 
1165  template <class TVar>
1166  void TPZTriangle::ComputeHCurlFaceDirections(TPZVec<TVar> &v1, TPZVec<TVar> &v2, int transformationId)
1167  {
1168  for (auto i=0; i< Dimension; i++){
1169  //v^{F,T} orthonormal vectors associated with face F and tangent to it.
1170  v1[i] = fTangentVectors[2*transformationId][i];
1171  v2[i] = fTangentVectors[2*transformationId + 1][i];
1172  }//for
1173  }
1174 
1175  void TPZTriangle::GetSideHDivDirections(TPZVec<int> &sides, TPZVec<int> &dir, TPZVec<int> &bilounao)
1176  {
1177  int nsides = NumSides()*2;
1178 
1179  sides.Resize(nsides);
1180  dir.Resize(nsides);
1181  bilounao.Resize(nsides);
1182 
1183  for (int is = 0; is<nsides; is++)
1184  {
1185  sides[is] = vectorsideorder[is];
1186  dir[is] = direcaoksioueta[is];
1187  bilounao[is] = bilinearounao[is];
1188  }
1189  }
1190 
1191  void TPZTriangle::GetSideHDivDirections(TPZVec<int> &sides, TPZVec<int> &dir, TPZVec<int> &bilounao, TPZVec<int> &sidevectors)
1192  {
1193  int nsides = NumSides()*2;
1194 
1195  sides.Resize(nsides);
1196  dir.Resize(nsides);
1197  bilounao.Resize(nsides);
1198 
1199  for (int is = 0; is<nsides; is++)
1200  {
1201  sides[is] = vectorsideorder[is];
1202  dir[is] = direcaoksioueta[is];
1203  bilounao[is] = bilinearounao[is];
1204  }
1205 
1206  for (int i=0; i<Dimension*NumSides(); i++) {
1207  sidevectors[i] = vectorsideorder[i];
1208  }
1209  }
1210 
1211  int TPZTriangle::ClassId() const{
1212  return Hash("TPZTriangle");
1213  }
1214 
1215  void TPZTriangle::Read(TPZStream& buf, void* context) {
1216 
1217  }
1218 
1219  void TPZTriangle::Write(TPZStream& buf, int withclassid) const {
1220 
1221  }
1222 }
1223 
1224 /**********************************************************************************************************************
1225  * The following are explicit instantiation of member function template of this class, both with class T=REAL and its
1226  * respective FAD<REAL> version. In other to avoid potential errors, always declare the instantiation in the same order
1227  * in BOTH cases. @orlandini
1228  **********************************************************************************************************************/
1229 
1230 template bool pztopology::TPZTriangle::CheckProjectionForSingularity<REAL>(const int &side, const TPZVec<REAL> &xiInterior);
1231 
1232 template void pztopology::TPZTriangle::MapToSide<REAL>(int side, TPZVec<REAL> &InternalPar, TPZVec<REAL> &SidePar, TPZFMatrix<REAL> &JacToSide);
1233 
1234 template void pztopology::TPZTriangle::BlendFactorForSide<REAL>(const int &, const TPZVec<REAL> &, REAL &, TPZVec<REAL> &);
1235 
1236 template void pztopology::TPZTriangle::TShape<REAL>(const TPZVec<REAL> &loc,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi);
1237 
1238 template void pztopology::TPZTriangle::ComputeHDivDirections<REAL>(TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions);
1239 
1240 template void pztopology::TPZTriangle::ComputeHCurlDirections<REAL>(TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions, const TPZVec<int> &transformationIds);
1241 
1242 template void pztopology::TPZTriangle::ComputeHCurlFaceDirections<REAL>(TPZVec<REAL> &v1, TPZVec<REAL> &v2, int transformationId);
1243 #ifdef _AUTODIFF
1244 
1245 template bool pztopology::TPZTriangle::CheckProjectionForSingularity<Fad<REAL> >(const int &side, const TPZVec<Fad<REAL> > &xiInterior);
1246 
1247 template void pztopology::TPZTriangle::MapToSide<Fad<REAL> >(int side, TPZVec<Fad<REAL> > &InternalPar, TPZVec<Fad<REAL> > &SidePar, TPZFMatrix<Fad<REAL> > &JacToSide);
1248 
1249 template void pztopology::TPZTriangle::BlendFactorForSide<Fad<REAL>>(const int &, const TPZVec<Fad<REAL>> &, Fad<REAL> &,
1250  TPZVec<Fad<REAL>> &);
1251 template void pztopology::TPZTriangle::TShape<Fad<REAL>>(const TPZVec<Fad<REAL>> &loc,TPZFMatrix<Fad<REAL>> &phi,TPZFMatrix<Fad<REAL>> &dphi);
1252 
1253 template void pztopology::TPZTriangle::ComputeHDivDirections<Fad<REAL>>(TPZFMatrix<Fad<REAL>> &gradx, TPZFMatrix<Fad<REAL>> &directions);
1254 
1255 template void pztopology::TPZTriangle::ComputeHCurlDirections<Fad<REAL>>(TPZFMatrix<Fad<REAL>> &gradx, TPZFMatrix<Fad<REAL>> &directions, const TPZVec<int> &transformationIds);
1256 
1257 template void pztopology::TPZTriangle::ComputeHCurlFaceDirections<Fad<REAL>>(TPZVec<Fad<REAL>> &v1, TPZVec<Fad<REAL>> &v2, int transformationId);
1258 #endif
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
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
static int direcaoksioueta[14]
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.
Handles the numerical integration for two-dimensional problems using triangular elements. Numerical Integration.
Definition: pzquad.h:102
Implements a vector class which allows to use external storage provided by the user. Utility.
Definition: pzquad.h:16
static bool IsInParametricDomain(const TPZVec< REAL > &pt, REAL tol=pztopology::gTolerance)
Verifies if the parametric point pt is in the element parametric domain.
clarg::argInt dimension("-d", "Matrices dimension M x M", 1000)
static int bilinearounao[14]
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
virtual void Identity()
Converts the matrix in an identity matrix.
Definition: pzmatrix.cpp:199
Definition: fad.h:54
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
static TVar ComputeDetjac(TPZFMatrix< TVar > &gradx)
Definition: pzaxestools.h:111
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
Contains the TPZTriangle class which defines the topology of a triangle.
static REAL tTriang[14][2]
Handles the numerical integration for one-dimensional problems. Numerical Integration.
Definition: pzquad.h:36
#define LOGPZ_FATAL(A, B)
Define log for fatal errors (cout)
Definition: pzlog.h:95
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
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
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
static const double tol
Definition: pzgeoprism.cpp:23
Groups all classes defining the structure of the master element.
Definition: PrismExtend.cpp:15
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
int64_t Rows() const
Returns number of rows.
Definition: pzmatrix.h:803
const TPZFMatrix< T > & Sum() const
Definition: pztrnsform.h:66
static int sidedimension[7]
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
static REAL MidSideNode[7][3]
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
static int vectorsideorder[14]
static REAL bTriang[14][2]
static Tvar Norm(const TPZVec< Tvar > &vetor)
Returns the L2-norm of the vector.
Definition: pznumeric.cpp:16
Integration rule for one point. Numerical Integration.
Definition: pzquad.h:478
MElementType
Define the element types.
Definition: pzeltype.h:52
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
static int highsides[7][3]
Definition: pzeltype.h:55
static int nsidenodes[7]
static REAL sidetosidetransforms[7][3][4][3]
static int nhighdimsides[7]
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