NeoPZ
tpztetrahedron.cpp
Go to the documentation of this file.
1 
6 #include "tpztetrahedron.h"
7 
8 #include "pzmanvector.h"
9 #include "pzerror.h"
10 #include "pzreal.h"
11 #include "pzquad.h"
12 #include "pzeltype.h"
13 #include "tpztriangle.h"
14 
15 #ifdef _AUTODIFF
16 #include "fad.h"
17 #endif
18 
19 #include "pzlog.h"
20 
21 #ifdef LOG4CXX
22 static LoggerPtr logger(Logger::getLogger("pz.topology.pztetrahedron"));
23 #endif
24 
25 using namespace std;
26 
27 namespace pztopology {
28 
29  static int nhighdimsides[15] = {7,7,7,7,3,3,3,3,3,3,1,1,1,1,0};
30 
31  int TPZTetrahedron::FaceNodes[4][3] = { {0,1,2},{0,1,3},{1,2,3},{0,2,3} };
32 
33  int TPZTetrahedron::SideNodes[6][2] = { {0,1},{1,2},{2,0},{0,3},{1,3},{2,3} };
34 
35  int TPZTetrahedron::ShapeFaceId[4][3] = { {0,1,2},{0,1,3},{1,2,3},{0,2,3} };
36 
37  int TPZTetrahedron::fPermutations[24][15] =
38  {
39  {0,1,2,3,4,5,6,7,8,9,10,11,12,13,14},/*000*/
40  {0,1,3,2,4,8,7,6,5,9,11,10,12,13,14},/*001*/
41  {0,2,1,3,6,5,4,7,9,8,10,13,12,11,14},/*002*/
42  {0,2,3,1,6,9,7,4,5,8,13,10,12,11,14},/*003*/
43  {0,3,1,2,7,8,4,6,9,5,11,13,12,10,14},/*004*/
44  {0,3,2,1,7,9,6,4,8,5,13,11,12,10,14},/*005*/
45  {1,0,2,3,4,6,5,8,7,9,10,11,13,12,14},/*006*/
46  {1,0,3,2,4,7,8,5,6,9,11,10,13,12,14},/*007*/
47  {1,2,0,3,5,6,4,8,9,7,10,12,13,11,14},/*008*/
48  {1,2,3,0,5,9,8,4,6,7,12,10,13,11,14},/*009*/
49  {1,3,0,2,8,7,4,5,9,6,11,12,13,10,14},/*010*/
50  {1,3,2,0,8,9,5,4,7,6,12,11,13,10,14},/*011*/
51  {2,0,1,3,6,4,5,9,7,8,10,13,11,12,14},/*012*/
52  {2,0,3,1,6,7,9,5,4,8,13,10,11,12,14},/*013*/
53  {2,1,0,3,5,4,6,9,8,7,10,12,11,13,14},/*014*/
54  {2,1,3,0,5,8,9,6,4,7,12,10,11,13,14},/*015*/
55  {2,3,0,1,9,7,6,5,8,4,13,12,11,10,14},/*016*/
56  {2,3,1,0,9,8,5,6,7,4,12,13,11,10,14},/*017*/
57  {3,0,1,2,7,4,8,9,6,5,11,13,10,12,14},/*018*/
58  {3,0,2,1,7,6,9,8,4,5,13,11,10,12,14},/*019*/
59  {3,1,0,2,8,4,7,9,5,6,11,12,10,13,14},/*020*/
60  {3,1,2,0,8,5,9,7,4,6,12,11,10,13,14},/*021*/
61  {3,2,0,1,9,6,7,8,5,4,13,12,10,11,14},/*022*/
62  {3,2,1,0,9,5,8,7,6,4,12,13,10,11,14} /*023*/
63  };
64 
65  static int sidedimension[15] = {0,0,0,0,1,1,1,1,1,1,2,2,2,2,3};
66 
67 
68  static int FaceConnectLocId[4][7] = { {0,1,2,4,5,6,10},{0,1,3,4,8,7,11},
69  {1,2,3,5,9,8,12},{0,2,3,6,9,7,13} };
70 
71  static int highsides[15][7] = {
72  {4,6,7,10,11,13,14},
73  {4,5,8,10,11,12,14},
74  {5,6,9,10,12,13,14},
75  {7,8,9,11,12,13,14},
76  {10,11,14},
77  {10,12,14},
78  {10,13,14},
79  {11,13,14},
80  {11,12,14},
81  {12,13,14},
82  {14},
83  {14},
84  {14},
85  {14},
86  {-999}
87  };
88 
89  static int nsidenodes[15] =
90  {
91  1,1,1,1,
92  2,2,2,2,2,2,
93  3,3,3,3,
94  4};
95 
96  static REAL sidetosidetransforms[15][7][4][3] = {
97  {
98  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
99  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
100  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
101  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,-99}},
102  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,-99}},
103  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,-99}},
104  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,0}}
105  },
106  {
107  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
108  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
109  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
110  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,-99}},
111  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,-99}},
112  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,-99}},
113  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,0}}
114  },
115  {
116  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
117  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
118  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-1,-99,-99}},
119  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,-99}},
120  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,-99}},
121  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,0,-99}},
122  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,0}}
123  },
124  {
125  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
126  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
127  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{1,-99,-99}},
128  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,-99}},
129  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,-99}},
130  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,1,-99}},
131  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{0,0,1}}
132  },
133  {
134  {{0.5,0,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0,-99}},
135  {{0.5,0,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0,-99}},
136  {{0.5,0,0},{-99,-99,-99},{-99,-99,-99},{0.5,0,0}}
137  },
138  {
139  {{-0.5,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,-99}},
140  {{0.5,0,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0,-99}},
141  {{-0.5,0.5,0},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,0}}
142  },
143  {
144  {{0,-0.5,-99},{-99,-99,-99},{-99,-99,-99},{0,0.5,-99}},
145  {{-0.5,0,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0,-99}},
146  {{0,-0.5,0},{-99,-99,-99},{-99,-99,-99},{0,0.5,0}}
147  },
148  {
149  {{0,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0,0.5,-99}},
150  {{0,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0,0.5,-99}},
151  {{0,0,0.5},{-99,-99,-99},{-99,-99,-99},{0,0,0.5}}
152  },
153  {
154  {{-0.5,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,-99}},
155  {{0,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0,0.5,-99}},
156  {{-0.5,0,0.5},{-99,-99,-99},{-99,-99,-99},{0.5,0,0.5}}
157  },
158  {
159  {{-0.5,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,-99}},
160  {{-0.5,0.5,-99},{-99,-99,-99},{-99,-99,-99},{0.5,0.5,-99}},
161  {{0,-0.5,0.5},{-99,-99,-99},{-99,-99,-99},{0,0.5,0.5}}
162  },
163  {
164  {{1,0,0},{0,1,0},{-99,-99,-99},{0,0,0}}
165  },
166  {
167  {{1,0,0},{0,0,1},{-99,-99,-99},{0,0,0}}
168  },
169  {
170  {{-1,1,0},{-1,0,1},{-99,-99,-99},{1,0,0}}
171  },
172  {
173  {{0,1,0},{0,0,1},{-99,-99,-99},{0,0,0}}
174  },
175  {
176  {{-99,-99,-99},{-99,-99,-99},{-99,-99,-99},{-99,-99,-99}}
177  }
178  };
179 
180  static REAL MidSideNode[15][3] = {
181  /*00*/{.0,.0,.0},/*01*/{1.,.0,.0},/*02*/{0.,1.,.0},/*03*/{.0,0.,1.0},/*04*/{.5,.0,.0},
182  /*05*/{.5,.5,.0},/*06*/{0.,.5,.0},/*07*/{0.,0.,.5},/*08*/{.5,0.,0.5},/*09*/{.0,.5,.5},
183  /*10*/{1./3.,1./3., 0. } ,/*11*/{1./3., .0 ,1./3.},
184  /*12*/{1./3.,1./3.,1./3.} ,/*13*/{ 0. ,1./3.,1./3.},/*14*/{1./4.,1./4.,1./4.} };
185 
186  static REAL bTetra[45][3] = // direcao perpendicular ao lado
187  {
188  {0,0,-1}, {1,0,-1}, {0,1,-1}, {0,0,-1}, {0.5,0.5,-1}, {0,0,-1}, {0,0,-1},// face 0
189  {0,-1,0}, {1,-1,0}, {0,-1,1}, {0,-1,0}, {0.5,-1,0.5}, {0,-1,0}, {0,-1,0},// face 1
190  {1,0,0} , {0,1,0} , {0,0,1}, {1,1,0}, {0,0.5,0.5}, {0.5,0,0.5}, {1,1,1} ,// face 2
191  {-1,0,0}, {-1,1,0}, {-1,0,1}, {-1,0,0}, {-1,0.5,0.5}, {-1,0,0} , {-1,0,0},// face 3
192  //interior
193  //aresta
194  {1,0,0},{-1,1,0},{0,-1,0}, {0,0,1}, {-1,0,1}, {0,-1,1},
195  //faces
196  {-1,0,0}, {0,1,0},// face 0
197  {1,0,0}, {0,0,1},// face 1
198  {1,0,-1}, {-1,2,-1},// face 2
199  {0,0,1}, {0,1,0} ,// face 3
200  //interior
201  {1,0,0} ,
202  {0,1,0} ,
203  {0,0,1}
204  };
205  static REAL t1Tetra[45][3] =
206  {
207  {-1,0,0},{-1,0,0},{-1,0,0},{-1,0,0},{-1,0,0},{-1,0,0},{-1,0,0},//face 0
208  {1,0,0}, {1,0,0}, {1,0,0}, {1,0,0}, {1,0,0}, {1,0,0}, {1,0,0}, //face 1
209  {1/sqrt(2),0,-1/sqrt(2)},{1/sqrt(2),0,-1/sqrt(2)},{1/sqrt(2),0,-1/sqrt(2)},{1/sqrt(2),0,-1/sqrt(2)},{1/sqrt(2),0,-1/sqrt(2)},{1/sqrt(2),0,-1/sqrt(2)},{1/sqrt(2),0,-1/sqrt(2)},//face 2
210  {0,0,1} ,{0,0,1} ,{0,0,1} ,{0,0,1} ,{0,0,1} ,{0,0,1} ,{0,0,1} ,//face 3
211  //interior
212  //aresta
213  {0,-1,0},{1,1,0},{0,0,-1}, {0,-1,0}, {0,-1,0}, {1,1,1},
214  //faces
215  {0,1,0}, {1,0,0},// face 0
216  {0,0,1}, {-1,0,0},// face 1
217  {-1,2,-1}, {-1,0,1},// face 2
218  {0,1,0}, {0,0,-1} ,// face 3
219  //interior
220  {0,1,0} ,
221  {0,0,1} ,
222  {1,0,0}
223 
224  };
225  static REAL t2Tetra[45][3] =
226  {
227  {0,1,0}, {0,1,0}, {0,1,0}, {0,1,0}, {0,1,0}, {0,1,0}, {0,1,0}, // face 0
228  {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1},// face 1
229  {-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},{-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},{-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},{-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},{-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},{-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},{-1/sqrt(6),2/sqrt(6),-1/sqrt(6)},// face 2
230  {0,1,0},{0,1,0},{0,1,0},{0,1,0},{0,1,0},{0,1,0},{0,1,0},// face 3
231  //interior
232  //aresta
233  {0,0,-1},{0,0,-1},{-1,0,0}, {-1,0,0}, {1,1,1}, {-1,0,0},
234  //faces
235  {0,0,-1}, {0,0,-1},// face 0
236  {0,-1,0}, {0,-1,0},// face 1
237  {1,1,1}, {1,1,1},// face 2
238  {-1,0,0}, {-1,0,0} ,// face 3
239  //interior
240  {0,0,1} ,
241  {1,0,0} ,
242  {0,1,0}
243  };
244 
245  static int vectorsideorderTe [45] =
246  {
247  0,1,2,4,5,6,10, //face 0
248  0,1,3,4,8,7,11,//face 1
249  1,2,3,5,9,8,12,//face 2
250  0,2,3,6,9,7,13,//face 3
251  4,5,6,7,
252  8,9,
253  10,10,//tg face 0
254  11,11,//tg face 1
255  12,12,//tg face 2
256  13,13,//tg face 3
257  14,14,14
258  };
259 
260  static int bilinearounao [45] =
261  {
262  0,0,0,0,0,0,0,0,0,0,
263  0,0,0,0,0,0,0,0,0,0,
264  0,0,0,0,0,0,0,0,1,1,
265  1,1,1,1,1,1,1,1,1,1,
266  1,1,1,1,1
267  }; //P*k Pk
268 
269 // static int bilinearounao [45] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; //Pk Pk-1
270 
271  static int direcaoksioueta [45] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,2};
272 
273  template<class T>
274  inline void TPZTetrahedron::TShape(const TPZVec<T> &loc,TPZFMatrix<T> &phi,TPZFMatrix<T> &dphi) {
275  T qsi = loc[0], eta = loc[1] , zeta = loc[2];
276 
277  phi(0,0) = 1.0-qsi-eta-zeta;
278  phi(1,0) = qsi;
279  phi(2,0) = eta;
280  phi(3,0) = zeta;
281 
282  dphi(0,0) = -1.0;
283  dphi(1,0) = -1.0;
284  dphi(2,0) = -1.0;
285  dphi(0,1) = 1.0;
286  dphi(1,1) = 0.0;
287  dphi(2,1) = 0.0;
288  dphi(0,2) = 0.0;
289  dphi(1,2) = 1.0;
290  dphi(2,2) = 0.0;
291  dphi(0,3) = 0.0;
292  dphi(1,3) = 0.0;
293  dphi(2,3) = 1.0;
294 
295  }
296  template<class T>
297  void TPZTetrahedron::BlendFactorForSide(const int &side, const TPZVec<T> &xi, T &blendFactor,
298  TPZVec<T> &blendFactorDxi) {
299  blendFactorDxi.Resize(TPZTetrahedron::Dimension, (T) 0);
300  blendFactor = 0;
301  const REAL tol = pztopology::GetTolerance();
302  #ifdef PZDEBUG
303  std::ostringstream sout;
304  if(side < NCornerNodes || side >= NSides){
305  sout<<"The side\t"<<side<<"is invalid. Aborting..."<<std::endl;
306  }
307 
308  if(!pztopology::TPZTetrahedron::IsInParametricDomain(xi,tol)){
309  sout<<"The method BlendFactorForSide expects the point xi to correspond to coordinates of a point";
310  sout<<" inside the parametric domain. Aborting...";
311  }
312 
313  if(!sout.str().empty()){
314  PZError<<std::endl<<sout.str()<<std::endl;
315 #ifdef LOG4CXX
316  LOGPZ_FATAL(logger,sout.str().c_str());
317 #endif
318  DebugStop();
319  }
320  #endif
321  //if the point is singular, the blend factor and its derivatives should be zero
322  if(!CheckProjectionForSingularity(side,xi)){
323  std::cout<<"Side projection is not regular and it should have been checked earlier. Aborting.."<<std::endl;
324  DebugStop();
325  blendFactor = 0;
326  for(int i = 0; i < blendFactorDxi.size(); i++) blendFactorDxi[i] = 0;
327  return;
328  }
329 
330  TPZFNMatrix<4, T> phi(NCornerNodes, 1);
331  TPZFNMatrix<8, T> dphi(Dimension, NCornerNodes);
332  TPZTetrahedron::TShape(xi, phi, dphi);
333  for (int i = 0; i < TPZTetrahedron::NSideNodes(side); i++) {
334  const int currentNode = TPZTetrahedron::SideNodeLocId(side, i);
335  blendFactor += phi(currentNode, 0);
336  blendFactorDxi[0] += dphi(0, currentNode);
337  blendFactorDxi[1] += dphi(1, currentNode);
338  blendFactorDxi[2] += dphi(2, currentNode);
339  }
340  switch (side) {
341  case 0:
342  case 1:
343  case 2:
344  case 3:
345  blendFactorDxi[0] = 0;
346  blendFactorDxi[1] = 0;
347  blendFactorDxi[2] = 0;
348  blendFactor = 0;
349  return;
350  case 4:
351  case 5:
352  case 6:
353  case 7:
354  case 8:
355  case 9:
356  blendFactorDxi[0] *= 2 * blendFactor;
357  blendFactorDxi[1] *= 2 * blendFactor;
358  blendFactorDxi[2] *= 2 * blendFactor;
359  blendFactor *= blendFactor;
360  return;
361  case 10:
362  case 11:
363  case 12:
364  case 13:
365  blendFactorDxi[0] *= 3 * blendFactor * blendFactor;
366  blendFactorDxi[1] *= 3 * blendFactor * blendFactor;
367  blendFactorDxi[2] *= 3 * blendFactor * blendFactor;
368  blendFactor *= blendFactor * blendFactor;
369  return;
370  case 14:
371  blendFactorDxi[0] = 0;
372  blendFactorDxi[1] = 0;
373  blendFactorDxi[2] = 0;
374  blendFactor = 1;
375  return;
376  }
377  }
378 
379  int TPZTetrahedron::NBilinearSides()
380  {
381  DebugStop();
382  return 0;
383  }
384 
385  void TPZTetrahedron::LowerDimensionSides(int side,TPZStack<int> &smallsides)
386  {
387  smallsides.Resize(0);
388  int nsidecon = NContainedSides(side);
389  int is;
390  for(is=0; is<nsidecon-1; is++)
391  smallsides.Push(ContainedSideLocId(side,is));
392  }
393 
394  void TPZTetrahedron::LowerDimensionSides(int side,TPZStack<int> &smallsides, int DimTarget)
395  {
396  smallsides.Resize(0);
397  int nsidecon = NContainedSides(side);
398  for(int is = 0; is < nsidecon - 1; is++) {
399  if (SideDimension(ContainedSideLocId(side,is)) == DimTarget) smallsides.Push(ContainedSideLocId(side,is));
400  }
401  }
402 
403  void TPZTetrahedron::HigherDimensionSides(int side, TPZStack<int> &high)
404  {
405  if(side <0 || side >= NSides) {
406  PZError << "TPZTetrahedron::HigherDimensionSides side "<< side << endl;
407  }
408  int is;
409  for(is=0; is<nhighdimsides[side]; is++) high.Push(highsides[side][is]);
410 
411  }
412 
413  int TPZTetrahedron::NSideNodes(int side)
414  {
415  return nsidenodes[side];
416  }
417  //Tentando criar o metodo
418  int TPZTetrahedron::NumSides(int dimension) { if(dimension<0 || dimension> 3) {
419  PZError << "TPZTetrahedron::NumSides. Bad parameter i.\n";
420  return 0;
421  }
422  if(dimension==0) return 4;
423  if(dimension==1) return 6;
424  if(dimension==2) return 4;
425  if(dimension==3) return 1;
426  return -1;
427  }
428 
429  int TPZTetrahedron::SideNodeLocId(int side, int node)
430  {
431  if(side<4 && node == 0) return side;
432  if(side>=4 && side < 10 && node <2) return SideNodes[side-4][node];
433  if(side >= 10 && side < 14 && node <3) return FaceNodes[side-10][node];
434  if(side ==14 && node < 4) return node;
435  PZError << "TPZTetrahedron::SideNodeLocId inconsistent side or node " << side
436  << ' ' << node << endl;
437  DebugStop();
438  return -1;
439 
440  }
441 
442  void TPZTetrahedron::CenterPoint(int side, TPZVec<REAL> &center) {
443  if (center.size()!=Dimension) {
444  DebugStop();
445  }
446  int i;
447  for(i=0; i<Dimension; i++) {
448  center[i] = MidSideNode[side][i];
449  }
450  }
451 
452  int TPZTetrahedron::SideDimension(int side) {
453  if(side<0 || side >= NSides) {
454  PZError << "TPZTetrahedron::SideDimension side " << side << endl;
455  return -1;
456  }
457  return sidedimension[side];
458  }
459 
460  TPZTransform<> TPZTetrahedron::SideToSideTransform(int sidefrom, int sideto)
461  {
462  if(sidefrom <0 || sidefrom >= NSides || sideto <0 || sideto >= NSides) {
463  PZError << "TPZTetrahedron::HigherDimensionSides sidefrom "<< sidefrom <<
464  ' ' << sideto << endl;
465  return TPZTransform<>(0);
466  }
467  if(sidefrom == sideto) {
468  return TPZTransform<>(sidedimension[sidefrom]);
469  }
470  if(sidefrom == NSides-1) {
471  return TransformElementToSide(sideto);
472  }
473  int nhigh = nhighdimsides[sidefrom];
474  int is;
475  for(is=0; is<nhigh; is++) {
476  if(highsides[sidefrom][is] == sideto) {
477  int dfr = sidedimension[sidefrom];
478  int dto = sidedimension[sideto];
479  TPZTransform<> trans(dto,dfr);
480  int i,j;
481  for(i=0; i<dto; i++) {
482  for(j=0; j<dfr; j++) {
483  trans.Mult()(i,j) = sidetosidetransforms[sidefrom][is][j][i];
484  }
485  trans.Sum()(i,0) = sidetosidetransforms[sidefrom][is][3][i];
486  }
487  return trans;
488  }
489  }
490  PZError << "TPZTetrahedron::SideToSideTransform highside not found sidefrom "
491  << sidefrom << ' ' << sideto << endl;
492  return TPZTransform<>(0);
493  }
494 
495  TPZTransform<> TPZTetrahedron::TransformElementToSide(int side){
496 
497  if(side<0 || side>14){
498  PZError << "TPZTetrahedron::TransformElementToSide called with side error\n";
499  return TPZTransform<>(0,0);
500  }
501 
502  TPZTransform<> t(sidedimension[side],3);
503  t.Mult().Zero();
504  t.Sum().Zero();
505  switch(side){
506  case 0:
507  case 1:
508  case 2:
509  case 3:
510  return t;
511  case 4:
512  t.Mult()(0,0) = 2.0;
513  t.Mult()(0,1) = 1.0;
514  t.Mult()(0,2) = 1.0;
515 
516 
517  t.Sum()(0,0) = -1.0;
518  return t;
519  case 5:
520  t.Mult()(0,0) = -1.0;
521  t.Mult()(0,1) = 1.0;
522  return t;
523 
524  case 6:
525  t.Mult()(0,0) = -1.0;
526  t.Mult()(0,1) = -2.0;
527  t.Mult()(0,2) = -1.0;
528 
529  t.Sum()(0,0) = 1.0;
530  return t;
531  case 7:
532  t.Mult()(0,0) = 1.0;
533  t.Mult()(0,1) = 1.0;
534  t.Mult()(0,2) = 2.0;
535 
536  t.Sum()(0,0) = -1.0;
537  return t;
538  case 8:
539  t.Mult()(0,0) = -1.0;
540  t.Mult()(0,2) = 1.0;
541 
542  return t;
543 
544  case 9:
545  t.Mult()(0,1) = -1.0;
546  t.Mult()(0,2) = 1.0;
547 
548  return t;
549 
550  case 10:
551  t.Mult()(0,0) = 1.0;
552  t.Mult()(1,1) = 1.0;
553  return t;
554 
555  case 11:
556  t.Mult()(0,0) = 1.0;
557  t.Mult()(1,2) = 1.0;
558  return t;
559  case 12:
560 
561  t.Mult()(0,0) = -1.0/3.0;
562  t.Mult()(0,1) = 2.0/3.0;
563  t.Mult()(0,2) = -1.0/3.0;
564  t.Mult()(1,0) = -1.0/3.0;
565  t.Mult()(1,1) = -1.0/3.0;
566  t.Mult()(1,2) = 2.0/3.0;
567 
568  t.Sum()(0,0) = 1.0/3.0;
569  t.Sum()(1,0)= 1.0/3.0;
570 
571 
572  return t;
573  case 13:
574  t.Mult()(0,1) = 1.0;
575  t.Mult()(1,2) = 1.0;
576  return t;
577  case 14:
578  t.Mult()(0,0) = 1.0;
579  t.Mult()(1,1) = 1.0;
580  t.Mult()(2,2) = 1.0;
581  return t;
582  }
583  return TPZTransform<>(0,0);
584 
585 // switch(side){
586 // case 0:
587 // case 1:
588 // case 2:
589 // case 3:
590 // return t;
591 // case 4:
592 // t.Mult()(0,0) = 2.0;
593 // t.Mult()(0,1) = 1.0;
594 // t.Mult()(0,2) = 1.0;
595 //
596 //
597 // t.Sum()(0,0) = -1.0;
598 // return t;
599 // case 5:
600 // t.Mult()(0,0) = -1.0;
601 // t.Mult()(0,1) = 1.0;
602 // return t;
603 //
604 // case 6:
605 // t.Mult()(0,0) = -1.0;
606 // t.Mult()(0,1) = -2.0;
607 // t.Mult()(0,2) = -1.0;
608 //
609 // t.Sum()(0,0) = 1.0;
610 // return t;
611 // case 7:
612 // t.Mult()(0,0) = 1.0;
613 // t.Mult()(0,1) = 1.0;
614 // t.Mult()(0,2) = 2.0;
615 //
616 // t.Sum()(0,0) = -1.0;
617 // return t;
618 // case 8:
619 // t.Mult()(0,0) = -1.0;
620 // t.Mult()(0,2) = 1.0;
621 //
622 // return t;
623 //
624 // case 9:
625 // t.Mult()(0,1) = -1.0;
626 // t.Mult()(0,2) = 1.0;
627 //
628 // return t;
629 // case 10:
630 // t.Mult()(0,0) = 2.0;
631 // t.Mult()(1,1) = 2.0;
632 //
633 // t.Sum()(0,0) = -1.0;
634 // t.Sum()(1,0)= -1.0;
635 //
636 // return t;
637 // case 11:
638 // t.Mult()(0,0) = 2.0;
639 // t.Mult()(1,2) = 2.0;
640 //
641 // t.Sum()(0,0) = -1.0;
642 // t.Sum()(1,0)= -1.0;
643 //
644 // return t;
645 // case 12:
646 // t.Mult()(0,0) = -2.0/3.0;
647 // t.Mult()(0,1) = 4.0/3.0;
648 // t.Mult()(0,2) = -2.0/3.0;
649 // t.Mult()(1,0) = -2.0/3.0;
650 // t.Mult()(1,1) = -2.0/3.0;
651 // t.Mult()(1,2) = 4.0/3.0;
652 //
653 // t.Sum()(0,0) = -1.0/3.0;
654 // t.Sum()(1,0)= -1.0/3.0;
655 // return t;
656 // case 13:
657 // t.Mult()(0,1) = 2.0;
658 // t.Mult()(1,2) = 2.0;
659 // t.Sum()(0,0) = -1.0;
660 // t.Sum()(1,0)= -1.0;
661 //
662 // return t;
663 //
664 // case 14:
665 // t.Mult()(0,0) = 1.0;
666 // t.Mult()(1,1) = 1.0;
667 // t.Mult()(2,2) = 1.0;
668 // return t;
669 // }
670 // return TPZTransform<>(0,0);
671  }
672 
673  TPZTransform<> TPZTetrahedron::TransformSideToElement(int side){
674 
675  if(side<0 || side>14){
676  PZError << "TPZTetrahedron::TransformSideToElement side out range\n";
677  return TPZTransform<>(0,0);
678  }
679  TPZTransform<> t(3,sidedimension[side]);
680  t.Mult().Zero();
681  t.Sum().Zero();
682 
683  switch(side){
684  case 0:
685  return t;
686  case 1:
687  t.Sum()(0,0) = 1.0;
688  return t;
689  case 2:
690  t.Sum()(1,0) = 1.0;
691  return t;
692  case 3:
693  t.Sum()(2,0) = 1.0;
694  return t;
695  case 4:
696  t.Mult()(0,0) = 0.5;
697  t.Sum() (0,0) = 0.5;
698  return t;
699  case 5:
700  t.Mult()(0,0) = -0.5;
701  t.Mult()(1,0) = 0.5;
702  t.Sum() (0,0) = 0.5;
703  t.Sum() (1,0) = 0.5;
704  return t;
705  case 6:
706  t.Mult()(1,0) = 0.5; //estava -0.5
707  t.Sum() (1,0) = 0.5;
708  return t;
709  case 7:
710  t.Mult()(2,0) = 0.5;
711  t.Sum() (2,0) = 0.5;
712  return t;
713  case 8:
714  t.Mult()(0,0) = -0.5;
715  t.Mult()(2,0) = 0.5;
716  t.Sum() (0,0) = 0.5;
717  t.Sum() (2,0) = 0.5;
718  return t;
719  case 9:
720  t.Mult()(1,0) = -0.5;
721  t.Mult()(2,0) = 0.5;
722  t.Sum() (1,0) = 0.5;
723  t.Sum() (2,0) = 0.5;
724  return t;
725  case 10:
726  t.Mult()(0,0) = 1.0;
727  t.Mult()(1,1) = 1.0;
728  return t;
729  case 11:
730  t.Mult()(0,0) = 1.0;
731  t.Mult()(2,1) = 1.0;
732  return t;
733  case 12:
734 
735  t.Mult()(0,0) = -1.0;
736  t.Mult()(0,1) = -1.0;
737  t.Mult()(1,0) = 1.0;
738  t.Mult()(2,1) = 1.0;
739  t.Sum() (0,0) = 1.0;
740 
741 // t.Mult()(0,0) = 1.0;
742 // // t.Mult()(0,1) = -1.0;
743 // t.Mult()(1,1) = 1.0;
744 // t.Mult()(2,0) = -1.0;
745 // t.Mult()(2,1) = -1.0;
746 // t.Sum() (2,0) = 1.0;
747  return t;
748  case 13:
749  t.Mult()(1,0) = 1.0;
750  t.Mult()(2,1) = 1.0;
751  return t;
752  case 14:
753  t.Mult()(0,0) = 1.0;
754  t.Mult()(1,1) = 1.0;
755  t.Mult()(2,2) = 1.0;
756  return t;
757 
758  }
759  return TPZTransform<>(0,0);
760  }
761 
762  TPZIntPoints * TPZTetrahedron::CreateSideIntegrationRule(int side, int order) {
763  if(side<0 || side>15) {
764  PZError << "TPZTetrahedron::CreateSideIntegrationRule. bad side number.\n";
765  return 0;
766  }
767  if(side<4) return new TPZInt1Point(order); // sides 0 to 3 are points (vertices)
768  if(side<10) return new TPZInt1d(order); // sides 4 to 9 are lines
769  if(side<14) { // sides 10 to 13 are triangles
770  return new TPZIntTriang(order);
771  }
772  if(side==14) { // integration of the element
773  return new IntruleType(order);
774  }
775  return 0;
776  }
777 
778 
779  MElementType TPZTetrahedron::Type()
780  {
781  return ETetraedro;
782  }
783 
784  MElementType TPZTetrahedron::Type(int side)
785  {
786  switch(side) {
787  case 0:
788  case 1:
789  case 2:
790  case 3:
791  return EPoint;
792  case 4:
793  case 5:
794  case 6:
795  case 7:
796  case 8:
797  case 9:
798  return EOned;
799  case 10:
800  case 11:
801  case 12:
802  case 13:
803  return ETriangle;
804  case 14:
805  return ETetraedro;
806  default:
807  return ENoType;
808  }
809  }
810 
811 
812  int TPZTetrahedron::NumSides() {
813  return NSides;
814  }
815 
816 
817  int TPZTetrahedron::NContainedSides(int side) {
818  if(side<0) return -1;
819  if(side<4) return 1;//cantos : 0 a 3
820  if(side<10) return 3;//lados : 4 a 9
821  if(side<14) return 7;//faces : 10 a 13
822  if(side==14) return 15;//centro : 14
823  return -1;
824  }
825 
826  int TPZTetrahedron::ContainedSideLocId(int side, int node) {
827  if(side<0 || side>15) return -1;
828  if(side<4) {
829  if(node==0) return side;
830  } else
831  if(side<7) {//4,5,6
832  int s = side-4;//0,1,2
833  if(!node) return s;//0,1,2
834  if(node==1) return (s+1)%3;//1,2,0
835  if(node==2) return side;//4,5,6
836  } else
837  if(side<10) {//7,8,9
838  int s = side-7;//0,1,2
839  if(!node) return s;//0,1,2
840  if(node==1) return 3;//4,4,4
841  if(node==2) return side;//7,8,9
842  } else
843  if(side<14) {//10 a 13
844  int s = side-10;
845  if(node<7) return FaceConnectLocId[s][node];
846  } else
847  if(side==14 && node<15){
848  return node;
849  }
850  PZError << "TPZShapeTetra::ContainedSideLocId called for node = "
851  << node << " and side = " << side << "\n";
852  return -1;
853  }
854 
855  bool TPZTetrahedron::IsInParametricDomain(const TPZVec<REAL> &pt, REAL tol){
856  const REAL qsi = pt[0];
857  const REAL eta = pt[1];
858  const REAL zeta = pt[2];
859  if( (qsi < 0. - tol) || (qsi > 1. + tol) ||
860  (eta < 0. - tol) || (eta > 1. + tol) ||
861  (zeta < 0. -tol) || (zeta > 1. +tol) ||
862  (qsi+eta+zeta > 1.+tol) ) {
863  return false;
864  }
865  else{
866  return true;
867  }
868 
869  }//method
870 
871 
873  void TPZTetrahedron::RandomPoint(TPZVec<REAL> &pt)
874  {
875  REAL val = (REAL) rand() / (RAND_MAX);
876  pt[0] = val;
877  val = (1.-pt[0]) * (REAL) rand() / (RAND_MAX);
878  pt[1] = val;
879  val = (1.-pt[0]-pt[1]) * (REAL) rand() / (RAND_MAX);
880  pt[2] = val;
881  }
882 
883  template<class T>
884  bool TPZTetrahedron::CheckProjectionForSingularity(const int &side, const TPZVec<T> &xiInterior) {
885 
886  T zero = pztopology::GetTolerance();
887 
888  T qsi = xiInterior[0], eta = xiInterior[1], zeta = xiInterior[2];
889  bool regularmap = true;
890  switch(side)
891  {
892  case 0:
893  case 1:
894  case 2:
895  case 3:
896  break;
897  case 4://1D
898  if(fabs((T)(eta + zeta - 1.)) < zero) regularmap = false;
899  break;
900  case 5://1D
901  if(fabs((T)(qsi + eta)) < zero) regularmap = false;
902  break;
903 
904  case 6://1D
905  if(fabs((T)(qsi + zeta - 1.)) < zero) regularmap = false;
906  break;
907  case 7://1D
908  if(fabs((T)(qsi + eta - 1.)) < zero) regularmap = false;
909  break;
910 
911  case 8://1D
912  if(fabs((T)(qsi + zeta)) < zero) regularmap = false;
913  break;
914 
915  case 9://1D
916  if(fabs((T)(eta + zeta)) < zero) regularmap = false;
917  break;
918 
919  case 10://2D
920  if(fabs((T)(zeta - 1.)) < zero) regularmap = false;
921  break;
922 
923  case 11://2D
924  if(fabs((T)(eta - 1.)) < zero) regularmap = false;
925  break;
926 
927  case 12://2D
928  if(fabs((T)(qsi+eta+zeta)) < zero) regularmap = false;
929  break;
930 
931  case 13://2D
932  if(fabs((T)(qsi - 1.)) < zero) regularmap = false;
933  break;
934  case 14:
935  break;
936  }
937  if(side > 14)
938  {
939  cout << "Cant compute CheckProjectionForSingularity method in TPZTetrahedron class!\nParameter (SIDE) must be between 4 and 13!\nMethod Aborted!\n";
940  DebugStop();
941  }
942  return regularmap;
943  }
944 
945  template<class T>
946  void TPZTetrahedron::MapToSide(int side, TPZVec<T> &InternalPar, TPZVec<T> &SidePar, TPZFMatrix<T> &JacToSide) {
947 
948  T qsi = InternalPar[0], eta = InternalPar[1], zeta = InternalPar[2];
949  if(!CheckProjectionForSingularity(side,InternalPar)){
950  std::cout<<"Side projection is not regular and it should have been checked earlier. Aborting.."<<std::endl;
951  DebugStop();
952  }
953  switch(side)
954  {
955  case 0:
956  case 1:
957  case 2:
958  case 3:
959  {
960  SidePar.Resize(0); JacToSide.Resize(0,0);
961  break;
962  }
963  case 4://1D
964  SidePar.Resize(1); JacToSide.Resize(1,3);
965  {
966  T den = (-1 + eta + zeta);
967  SidePar[0] = -1 - (2*qsi)/den;
968  den *=den;
969  JacToSide(0,0) = -2/(-1 + eta + zeta);
970  JacToSide(0,1) = (2*qsi)/den;
971  JacToSide(0,2) = (2*qsi)/den;
972  }
973  break;
974 
975  case 5://1D
976  SidePar.Resize(1); JacToSide.Resize(1,3);
977  {
978  T den = eta + qsi;
979  SidePar[0] = -1 + (2*eta)/den;
980  den *= den;
981  JacToSide(0,0) = (-2*eta)/den;
982  JacToSide(0,1) = (2*qsi)/den;
983  JacToSide(0,2) = 0;
984  }
985  break;
986 
987  case 6://1D
988  SidePar.Resize(1); JacToSide.Resize(1,3);
989  {
990  T den = (-1 + qsi + zeta);
991  SidePar[0] = 1 + (2*eta)/den;
992 
993  JacToSide(0,0) = (-2*eta)/(den * den);
994  JacToSide(0,1) = 2/den;
995  JacToSide(0,2) = (-2*eta)/(den*den);
996 
997  }
998  break;
999 
1000  case 7://1D
1001  SidePar.Resize(1); JacToSide.Resize(1,3);
1002  {
1003  T den = (-1 + eta + qsi);
1004  SidePar[0] = -1 - (2*zeta)/den;
1005  JacToSide(0,0) = (2*zeta)/(den*den);
1006  JacToSide(0,1) = (2*zeta)/(den*den);
1007  JacToSide(0,2) = -2/den;
1008 
1009  }
1010  break;
1011 
1012  case 8://1D
1013  SidePar.Resize(1); JacToSide.Resize(1,3);
1014  {
1015  T den = (qsi + zeta);
1016  SidePar[0] = -1 + (2*zeta)/den;
1017  den *= den;
1018  JacToSide(0,0) = (-2*zeta)/den;
1019  JacToSide(0,1) = 0;
1020  JacToSide(0,2) = (2*qsi)/den;
1021  }
1022  break;
1023 
1024  case 9://1D
1025  SidePar.Resize(1); JacToSide.Resize(1,3);
1026  {
1027  T den = eta + zeta;
1028  SidePar[0] = -1 + (2*zeta)/den;
1029  den *= den;
1030  JacToSide(0,0) = 0;
1031  JacToSide(0,1) = (-2*zeta)/den;
1032  JacToSide(0,2) = (2*eta)/den;
1033  }
1034  break;
1035 
1036  case 10://2D
1037  SidePar.Resize(2); JacToSide.Resize(2,3);
1038  {
1039  SidePar[0] = qsi/(1 - zeta);
1040  SidePar[1] = eta/(1 - zeta);
1041  JacToSide(0,0) = 1/(1 - zeta);
1042  JacToSide(0,1) = 0;
1043  JacToSide(0,2) = qsi/((-1 + zeta)*(-1 + zeta));
1044  JacToSide(1,0) = 0;
1045  JacToSide(1,1) = 1/(1 - zeta);
1046  JacToSide(1,2) = eta/((-1 + zeta)*(-1 + zeta));
1047  }
1048  break;
1049 
1050  case 11://2D
1051  SidePar.Resize(2); JacToSide.Resize(2,3);
1052  {
1053  SidePar[0] = qsi/(1 - eta);
1054  SidePar[1] = zeta/(1 - eta);
1055  JacToSide(0,0) = 1/(1 - eta);
1056  JacToSide(0,1) = qsi/((-1 + eta)*(-1 + eta));
1057  JacToSide(0,2) = 0;
1058  JacToSide(1,0) = 0;
1059  JacToSide(1,1) = zeta/((-1 + eta)*(-1 + eta));
1060  JacToSide(1,2) = 1/(1 - eta);
1061  }
1062  break;
1063 
1064  case 12://2D
1065  SidePar.Resize(2); JacToSide.Resize(2,3);
1066  {
1067  SidePar[0] = eta/(eta + qsi + zeta);
1068  SidePar[1] = zeta/(eta + qsi + zeta);
1069  T den = ((eta + qsi + zeta)*(eta + qsi + zeta));
1070  JacToSide(0,0) = -(eta/den);
1071  JacToSide(0,1) = (qsi + zeta)/den;
1072  JacToSide(0,2) = -(eta/den);
1073  JacToSide(1,0) = -(zeta/den);
1074  JacToSide(1,1) = -(zeta/den);
1075  JacToSide(1,2) = (eta + qsi)/den;
1076  }
1077  break;
1078 
1079  case 13://2D
1080  SidePar.Resize(2); JacToSide.Resize(2,3);
1081  {
1082  SidePar[0] = eta/(1 - qsi);
1083  SidePar[1] = zeta/(1 - qsi);
1084  JacToSide(0,0) = eta/((-1 + qsi)*(-1 + qsi));
1085  JacToSide(0,1) = 1/(1 - qsi);
1086  JacToSide(0,2) = 0;
1087  JacToSide(1,0) = zeta/((-1 + qsi)*(-1 + qsi));
1088  JacToSide(1,1) = 0;
1089  JacToSide(1,2) = 1/(1 - qsi);
1090  }
1091  break;
1092  case 14:
1093  SidePar = InternalPar;
1094  JacToSide.Resize(3, 3);
1095  JacToSide.Identity();
1096  break;
1097  }
1098  if(side > 14)
1099  {
1100  cout << "Cant compute MapToSide method in TPZTetrahedron class!\nParameter (SIDE) must be between 4 and 13!\nMethod Aborted!\n";
1101  cout << "This should have been caught earlier in the execution, there is something wrong.\n";
1102  cout << "Check method TPZTetrahedron::CheckProjectionForSingularity<T>\n";
1103  DebugStop();
1104  }
1105  }
1106 
1107  void TPZTetrahedron::ParametricDomainNodeCoord(int node, TPZVec<REAL> &nodeCoord)
1108  {
1109  if(node > NCornerNodes)
1110  {
1111  DebugStop();
1112  }
1113  nodeCoord.Resize(Dimension, 0.);
1114  switch (node) {
1115  case (0):
1116  {
1117  nodeCoord[0] = 0.;
1118  nodeCoord[1] = 0.;
1119  nodeCoord[2] = 0.;
1120  break;
1121  }
1122  case (1):
1123  {
1124  nodeCoord[0] = 1.;
1125  nodeCoord[1] = 0.;
1126  nodeCoord[2] = 0.;
1127  break;
1128  }
1129  case (2):
1130  {
1131  nodeCoord[0] = 0.;
1132  nodeCoord[1] = 1.;
1133  nodeCoord[2] = 0.;
1134  break;
1135  }
1136  case (3):
1137  {
1138  nodeCoord[0] = 0.;
1139  nodeCoord[1] = 0.;
1140  nodeCoord[2] = 1.;
1141  break;
1142  }
1143  default:
1144  {
1145  DebugStop();
1146  break;
1147  }
1148  }
1149  }
1150 
1157  int TPZTetrahedron::GetTransformId(TPZVec<int64_t> &id)
1158  {
1159  LOGPZ_ERROR(logger,"Please implement me")
1160  return -1;
1161  }
1162 
1169  int TPZTetrahedron::GetTransformId(int side, TPZVec<int64_t> &id)
1170  {
1171  switch (side) {
1172  case 0:
1173  case 1:
1174  case 2:
1175  case 3:
1176 
1177  return 0;
1178  break;
1179  case 4:
1180  case 5:
1181  case 6:
1182  case 7:
1183  case 8:
1184  case 9:
1185 
1186  {
1187  int in1 = ContainedSideLocId(side,0);
1188  int in2 = ContainedSideLocId(side,1);
1189  return id[in1]<id[in2] ? 0 : 1;
1190  }
1191  break;
1192  case 10:
1193  case 11:
1194  case 12:
1195  case 13:
1196  {
1197  TPZManVector<int64_t,3> locid(3);
1198  int i;
1199  for(i=0; i<3; i++) locid[i] = id[ContainedSideLocId(side,i)];
1200  return pztopology::TPZTriangle::GetTransformId(locid);
1201  }
1202  break;
1203  case 14:
1204  {
1205  return 0;//that is not really true
1206  }
1207  default:
1208  DebugStop();
1209  }
1210  return -1;
1211  }
1212 
1220  void TPZTetrahedron::GetSideHDivPermutation(int transformationid, TPZVec<int> &permgather)
1221  {
1222  // Not complete
1223  DebugStop();
1224 
1225 //#ifdef PZDEBUG
1226 // if (SideDimension(side) != 2) {
1227 // DebugStop();
1228 // }
1229 //#endif
1230 // permgather.Resize(7);
1231 // TPZManVector<int64_t,7> locids(3);
1232 // for (int in=0; in<3; in++) {
1233 // locids[in] = id[SideNodeLocId(side, in)];
1234 // }
1235 // int transformid = pztopology::TPZTriangle::GetTransformId(locids);
1236 // pztopology::TPZTriangle::GetHDivGatherPermute(transformid,permgather);
1237  }
1238 
1239 
1240  void computedirectionsT3(int inicio, int fim, TPZFMatrix<REAL> &bvec, TPZFMatrix<REAL> &t1vec,
1241  TPZFMatrix<REAL> &t2vec, TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions);
1242 
1243  void computedirectionsT3(int inicio, int fim, TPZFMatrix<REAL> &bvec, TPZFMatrix<REAL> &t1vec,
1244  TPZFMatrix<REAL> &t2vec, TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions)
1245  {
1246  REAL detgrad = 0.0;
1247  TPZVec<REAL> u(3);
1248  TPZVec<REAL> v(3);
1249  TPZVec<REAL> uxv(3);// result
1250  int cont = 0;
1251 
1252  for (int ivet=inicio; ivet<=fim; ivet++)
1253  {
1254  for (int ilin=0; ilin<3; ilin++)
1255  {
1256  u[ilin] = t1vec(ilin,ivet);
1257  v[ilin] = t2vec(ilin,ivet);
1258  }
1259  TPZVec<REAL> e2(3);
1260  detgrad = 0.0;
1261  REAL normaX0xX1 = 0.0;
1262  //TPZNumeric::ProdVetorial(u,v,e2);
1263  e2[0] = u[1]*v[2]-u[2]*v[1];
1264  e2[1] = -(u[0]*v[2]-v[0]*u[2]);
1265  e2[2] = u[0]*v[1]-v[0]*u[1];
1266 
1267  // calc do v gradx*b
1268  TPZManVector<REAL,3> dxt1(3,0.),dxt2(3,0.),dxt3(3,0.),Vvec(3,0.);
1269  REAL be2 = 0.0, ne2 = 0.0;
1270  for(int i=0;i<3;i++)
1271  {
1272  ne2 += e2[i]*e2[i];
1273  }
1274  ne2 = sqrt(fabs(ne2));
1275  for (int il=0; il<3; il++)
1276  {
1277  for (int i = 0 ; i<3; i++)
1278  {
1279  dxt1[il] += gradx(il,i) * t1vec(i,ivet);
1280  dxt2[il] += gradx(il,i) * t2vec(i,ivet);
1281  dxt3[il] += gradx(il,i) * e2[i]/ne2;
1282  Vvec[il] += gradx(il,i) * bvec(i,ivet);
1283  }
1284  be2 += bvec(il,ivet)*e2[il]/ne2;
1285  }
1286  TPZManVector<REAL,3> normal(3,0.);
1287  //TPZNumeric::ProdVetorial(dxt1,dxt2,normal);
1288  normal[0] = dxt1[1]*dxt2[2]-dxt1[2]*dxt2[1];
1289  normal[1] = -(dxt1[0]*dxt2[2]-dxt2[0]*dxt1[2]);
1290  normal[2] = dxt1[0]*dxt2[1]-dxt2[0]*dxt1[1];
1291 
1292  for (int pos=0; pos<3; pos++)
1293  {
1294  detgrad += normal[pos]*dxt3[pos];//uxv[pos]*gradx.GetVal(pos, 2);
1295  normaX0xX1 += normal[pos]*normal[pos]; //uxv[pos]*uxv[pos];
1296  }
1297  TPZFMatrix<REAL> Wvec(3,1);
1298  detgrad = fabs(detgrad);
1299  normaX0xX1 = sqrt(normaX0xX1);
1300 
1301  for (int il=0; il<3; il++)
1302  {
1303  Wvec(il,0) = Vvec[il]*normaX0xX1/(detgrad*be2);
1304  directions(il,cont) = Wvec(il,0);
1305  }
1306  cont++;
1307  }
1308 
1309 
1310  }
1311 
1312 
1313  void TPZTetrahedron::ComputeDirections(int side, TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions, TPZVec<int> &sidevectors)
1314  {
1315  if(gradx.Cols()!=3)
1316  { std::cout << "Gradient dimensions are not compatible with this topology" << std::endl;
1317  DebugStop();
1318  }
1319  TPZFMatrix<REAL> bvec(3, 45);
1320  int numvec = bvec.Cols();
1321  TPZFMatrix<REAL> t1vec(3, numvec);
1322  TPZFMatrix<REAL> t2vec(3,numvec);
1323 
1324  directions.Redim(3, numvec);
1325 
1326  for (int lin = 0; lin<numvec ; lin++)
1327  {
1328  for(int col = 0;col<3;col++)
1329  {
1330  bvec.PutVal(col, lin, bTetra[lin][col]);
1331  t1vec.PutVal(col, lin, t1Tetra[lin][col]);
1332  t2vec.PutVal(col, lin, t2Tetra[lin][col]);
1333  }
1334  }
1335 
1336  // calcula os vetores
1337  switch (side) {
1338  case 10:
1339  {
1340  directions.Resize(3, 7);
1341  sidevectors.Resize(7);
1342  int inicio = 0, fim = 6;
1343  computedirectionsT3( inicio, fim, bvec, t1vec, t2vec, gradx, directions);
1344  int diff = fim-inicio+1;
1345  for (int ip = 0; ip < diff; ip++) {
1346  sidevectors[ip] = vectorsideorderTe[ip+inicio];
1347  }
1348  }
1349  break;
1350  case 11:
1351  {
1352  directions.Resize(3, 7);
1353  sidevectors.Resize(7);
1354  int inicio = 7, fim = 13;
1355  computedirectionsT3( inicio, fim, bvec, t1vec, t2vec, gradx, directions);
1356  int diff = fim-inicio+1;
1357  for (int ip = 0; ip < diff; ip++) {
1358  sidevectors[ip] = vectorsideorderTe[ip+inicio];
1359  }
1360  }
1361  break;
1362  case 12:
1363  {
1364  directions.Resize(3, 7);
1365  sidevectors.Resize(7);
1366  int inicio = 14, fim = 20;
1367  computedirectionsT3( inicio, fim, bvec, t1vec, t2vec, gradx, directions);
1368  int diff = fim-inicio+1;
1369  for (int ip = 0; ip < diff; ip++) {
1370  sidevectors[ip] = vectorsideorderTe[ip+inicio];
1371  }
1372  }
1373  break;
1374  case 13:
1375  {
1376  directions.Resize(3, 7);
1377  sidevectors.Resize(7);
1378  int inicio = 21, fim = 27;
1379  computedirectionsT3( inicio, fim, bvec, t1vec, t2vec, gradx, directions);
1380  int diff = fim-inicio+1;
1381  for (int ip = 0; ip < diff; ip++) {
1382  sidevectors[ip] = vectorsideorderTe[ip+inicio];
1383  }
1384  }
1385  break;
1386  case 14:
1387  {
1388  directions.Resize(3, 17);
1389  sidevectors.Resize(17);
1390  int inicio = 28, fim = 44;
1391  computedirectionsT3( inicio, fim, bvec, t1vec, t2vec, gradx, directions);
1392  int diff = fim-inicio+1;
1393  for (int ip = 0; ip < diff; ip++) {
1394  sidevectors[ip] = vectorsideorderTe[ip+inicio];
1395  }
1396  }
1397  break;
1398 
1399  default:
1400  break;
1401  }
1402 
1403 
1404  }
1405 
1406  template <class TVar>
1407  void TPZTetrahedron::ComputeHDivDirections(TPZFMatrix<TVar> &gradx, TPZFMatrix<TVar> &directions)
1408  {
1409  TVar detjac = TPZAxesTools<TVar>::ComputeDetjac(gradx);
1410 
1411  TPZManVector<TVar,3> v1(3),v2(3),v3(3),v1v2(3),v3v1(3),v2v3(3),vdiagxy(3),vi(3),vivdiagxy(3);
1412 
1413  for (int i=0; i<3; i++) {
1414  v1[i] = gradx(i,0);
1415  v2[i] = gradx(i,1);
1416  v3[i] = gradx(i,2);
1417  vdiagxy[i] = (gradx(i,0)-gradx(i,1));
1418  vi[i] = gradx(i,2)-gradx(i,0);//gradx(i,2)-0.5*(gradx(i,0)+gradx(i,1));
1419  }
1420 
1421 
1428  TVar Nv1v2 = 1.0;
1429  TVar Nv2v3 = 1.0;
1430  TVar Nv3v1 = 1.0;
1431  TVar Nvivdiagb = 1.0;
1432 
1433  {
1434  // the above constants are wrong
1435  for (int i=0; i<3; i++) {
1436  v1[i] /= detjac;
1437  v2[i] /= detjac;
1438  v3[i] /= detjac;
1439  }
1440  for (int i=0; i<3; i++)
1441  {
1442 
1443  //face 0
1444  directions(i,0) = -v3[i];
1445  directions(i,1) = (v1[i]-v3[i]);
1446  directions(i,2) = (v2[i]-v3[i]);
1447  directions(i,3) = (directions(i,0)+directions(i,1))/2.;
1448  directions(i,4) = (directions(i,1)+directions(i,2))/2.;
1449  directions(i,5) = (directions(i,0)+directions(i,2))/2.;
1450  directions(i,6) = (directions(i,3)+directions(i,4)+directions(i,5))/3.;
1451  //face 1
1452  directions(i,7) = -v2[i];
1453  directions(i,8) = (v1[i]-v2[i]);
1454  directions(i,9) = (v3[i]-v2[i]);
1455  directions(i,10) = (directions(i,7)+directions(i,8))/2.;
1456  directions(i,11) = (directions(i,8)+directions(i,9))/2.;
1457  directions(i,12) = (directions(i,7)+directions(i,9))/2.;
1458  directions(i,13) = (directions(i,10)+directions(i,11)+directions(i,12))/3.;
1459  //face 2 face diagonal
1460 
1461  directions(i,14) = v1[i];
1462  directions(i,15) = v2[i];
1463  directions(i,16) = v3[i];
1464  directions(i,17) = (directions(i,14)+directions(i,15))/2.;
1465  directions(i,18) = (directions(i,15)+directions(i,16))/2.;
1466  directions(i,19) = (directions(i,14)+directions(i,16))/2.;
1467  directions(i,20) = (directions(i,17)+directions(i,18)+directions(i,19))/3.;
1468  //face 3
1469  directions(i,21) = -v1[i];
1470  directions(i,22) = (v2[i]-v1[i]);
1471  directions(i,23) = (v3[i]-v1[i]);
1472  directions(i,24) = (directions(i,21)+directions(i,22))/2.;
1473  directions(i,25) = (directions(i,22)+directions(i,23))/2.;
1474  directions(i,26) = (directions(i,21)+directions(i,23))/2.;
1475  directions(i,27) = (directions(i,24)+directions(i,25)+directions(i,26))/3.;
1476 
1477  //arestas
1478  directions(i,28) = v1[i];
1479  directions(i,29) = (v2[i]-v1[i]);
1480  directions(i,30) = -v2[i];
1481  directions(i,31) = v3[i];
1482  directions(i,32) = (v3[i]-v1[i]);
1483  directions(i,33) = (v3[i]-v2[i]);
1484 
1485  //faces
1486  directions(i,34) = v1[i];
1487  directions(i,35) = v2[i];
1488  directions(i,36) = v1[i];
1489  directions(i,37) = v3[i];
1490  directions(i,38) = (v2[i]-v1[i]);
1491  directions(i,39) = (v3[i]-v1[i]);//v3[i]-0.5*(v1[i]+v2[i]);//
1492  directions(i,40) = v2[i];
1493  directions(i,41) = v3[i];
1494 
1495  directions(i,42) = v1[i];
1496  directions(i,43) = v2[i];
1497  directions(i,44) = v3[i];
1498 
1499  }
1500 
1501  }
1502 
1503  }
1504 
1505 
1506  void TPZTetrahedron::GetSideHDivDirections(TPZVec<int> &sides, TPZVec<int> &dir, TPZVec<int> &bilounao)
1507  {
1508  int nsides = NumSides()*3;
1509 
1510  sides.Resize(nsides);
1511  dir.Resize(nsides);
1512  bilounao.Resize(nsides);
1513 
1514  for (int is = 0; is<nsides; is++)
1515  {
1516  sides[is] = vectorsideorderTe[is];
1517  dir[is] = direcaoksioueta[is];
1518  bilounao[is] = bilinearounao[is];
1519  }
1520  }
1521 
1522  void TPZTetrahedron::GetSideHDivDirections(TPZVec<int> &sides, TPZVec<int> &dir, TPZVec<int> &bilounao, TPZVec<int> &sidevectors)
1523  {
1524  int nsides = NumSides()*3;
1525 
1526  sides.Resize(nsides);
1527  dir.Resize(nsides);
1528  bilounao.Resize(nsides);
1529 
1530  for (int is = 0; is<nsides; is++)
1531  {
1532  sides[is] = vectorsideorderTe[is];
1533  dir[is] = direcaoksioueta[is];
1534  bilounao[is] = bilinearounao[is];
1535  }
1536  for (int i=0; i<Dimension*NumSides(); i++) {
1537  sidevectors[i] = vectorsideorderTe[i];
1538  }
1539  }
1540 
1541  template <class TVar>
1542  void TPZTetrahedron::ComputeHCurlDirections(TPZFMatrix<TVar> &gradx, TPZFMatrix<TVar> &directions, const TPZVec<int> &transformationIds)
1543  {
1544  TPZManVector<TVar,3> v1(3),v2(3),v3(3);
1545 
1546  for (int i=0; i<3; i++) {
1547  v1[i] = gradx(i,0);
1548  v2[i] = gradx(i,1);
1549  v3[i] = gradx(i,2);
1550  }
1551  constexpr int nFaces{4},nEdges{6};
1552  //edges 4, 5 ,6,7, 8 , 9
1553  constexpr REAL edgeLength[nEdges]{1,M_SQRT2,1,1,M_SQRT2,M_SQRT2};
1554  constexpr REAL sqrt3 = 1.73205080756887729352744634150587237;
1555  //faces 10, 11, 12 , 13
1556  constexpr REAL faceArea[nFaces]{0.5,0.5,0.5*sqrt3,0.5};
1557  TPZManVector<REAL,nEdges> edgeSign(nEdges,0);
1558  for(auto iEdge = 0; iEdge < nEdges; iEdge++){
1559  edgeSign[iEdge] = transformationIds[iEdge] == 0 ? 1 : -1;
1560  }
1561  for (int i=0; i<3; i++)
1562  {
1563  //v^{e,a} constant vector fields associated with edge e and vertex a
1564  //they are defined in such a way that v^{e,a} is normal to the edge \hat{e}
1565  //adjacent to edge e by the vertex a. the tangential component is set to be 1 /edgeLength[e]
1566 
1567  directions(i,0) = (v1[i]) * edgeSign[0] / edgeLength[0];//edge 4 vertex 0
1568  directions(i,1) = (v1[i]+v2[i]+v3[i]) * edgeSign[0] / edgeLength[0];//edge 4 vertex 1
1569  directions(i,2) = (v2[i]*M_SQRT2) * edgeSign[1] / edgeLength[1];//edge 5 vertex 1
1570  directions(i,3) = (-v1[i]*M_SQRT2) * edgeSign[1] / edgeLength[1];//edge 5 vertex 2
1571  directions(i,4) = (v1[i]+v2[i]+v3[i]) * -1 * edgeSign[2] / edgeLength[2]; //edge 6 vertex 2
1572  directions(i,5) = (-v2[i]) * edgeSign[2] / edgeLength[2]; //edge 6 vertex 0
1573  directions(i,6) = (v3[i]) * edgeSign[3] / edgeLength[3]; //edge 7 vertex 0
1574  directions(i,7) = (v1[i]+v2[i]+v3[i]) * edgeSign[3] / edgeLength[3]; //edge 7 vertex 3
1575  directions(i,8) = v3[i] * M_SQRT2 * edgeSign[4] / edgeLength[4];//edge 8 vertex 1
1576  directions(i,9) = -v1[i] * M_SQRT2 * edgeSign[4] / edgeLength[4];//edge 8 vertex 3
1577  directions(i,10) = v3[i] * M_SQRT2 * edgeSign[5] / edgeLength[5];//edge 9 vertex 2
1578  directions(i,11) = -v2[i] * M_SQRT2 * edgeSign[5] / edgeLength[5];//edge 9 vertex 3
1579 
1580  //v^{e,T} constant vector fields associated with edge e and aligned with it
1581  directions(i,12) = (v1[i]) * edgeSign[0] / edgeLength[0];//edge 4
1582  directions(i,13) = ( (v2[i] - v1[i]) * M_SQRT1_2) * edgeSign[1] / edgeLength[1];//edge 5
1583  directions(i,14) = (-v2[i]) * edgeSign[2] / edgeLength[2];//edge 6
1584  directions(i,15) = (v3[i]) * edgeSign[3] / edgeLength[3];//edge 7
1585  directions(i,16) = (v3[i]-v1[i]) * M_SQRT1_2 * edgeSign[4] / edgeLength[4];//edge 8
1586  directions(i,17) = (v3[i]-v2[i]) * M_SQRT1_2 * edgeSign[5] / edgeLength[5];//edge 9
1587 
1588  //v^{F,e} constant vector fields associated with face F and edge e
1589  //they are defined in such a way that v^{F,e} is normal to the face \hat{F}
1590  //adjacent to face F by edge e
1591  directions(i,18) = v2[i] * edgeSign[4 - NCornerNodes] / faceArea[0];//face 10 edge 4
1592  directions(i,19) = -1 * (v1[i]+v2[i]+v3[i]) * M_SQRT1_2 * edgeSign[5 - NCornerNodes] / faceArea[0];//face 10 edge 5
1593  directions(i,20) = v1[i] * edgeSign[6 - NCornerNodes] / faceArea[0];//face 10 edge 6
1594  directions(i,21) = v3[i] * edgeSign[4 - NCornerNodes] / faceArea[1];//face 11 edge 4
1595  directions(i,22) = -1 * (v1[i]+v2[i]+v3[i]) * M_SQRT1_2 * edgeSign[8 - NCornerNodes] / faceArea[1];//face 11 edge 8
1596  directions(i,23) = -1 * v1[i] * edgeSign[7 - NCornerNodes] / faceArea[1];//face 11 edge 7
1597  directions(i,24) = v3[i] * sqrt3 * edgeSign[5 - NCornerNodes] / faceArea[2];//face 12 edge 5
1598  directions(i,25) = v1[i] * sqrt3 * M_SQRT1_2 * edgeSign[9 - NCornerNodes] / faceArea[2];//face 12 edge 9
1599  directions(i,26) = -1 * v2[i] * sqrt3 * edgeSign[8 - NCornerNodes] / faceArea[2];//face 12 edge 8
1600  directions(i,27) = -1 * v3[i] * edgeSign[6 - NCornerNodes] / faceArea[3];//face 13 edge 6
1601  directions(i,28) = -1* (v1[i]+v2[i]+v3[i]) * M_SQRT1_2 * edgeSign[9 - NCornerNodes] / faceArea[3];//face 13 edge 9
1602  directions(i,29) = -1 * v2[i] * edgeSign[7 - NCornerNodes] / faceArea[3];//face 13 edge 7
1603 
1604  //v^(F,T} vectors are calculated afterwards
1605 
1606  //v^{F,orth} vector associated with face F and normal to it
1607  directions(i,38) = -v3[i];//face 10
1608  directions(i,39) = -v2[i];//face 11
1609  directions(i,40) = (v1[i]+v2[i]+v3[i])/sqrt3;//face 12
1610  directions(i,41) = -v1[i];//face 13
1611 
1612  //v^{K,3}
1613  directions(i,42) = v1[i];
1614  directions(i,43) = v2[i];
1615  directions(i,44) = v3[i];
1616  }
1617  TPZManVector<REAL,2> vft1(2,0), vft2(2,0);
1618  constexpr auto firstVftVec = 30;
1619  //v^{F,T} orthonormal vectors associated with face F and tangent to it.
1620  for(auto iFace = 0; iFace < nFaces; iFace ++){
1621  TPZTriangle::ComputeHCurlFaceDirections(vft1,vft2,transformationIds[nEdges + iFace]);
1622  directions(0,firstVftVec+2*iFace) = 0;directions(1,firstVftVec+2*iFace) = 0;directions(2,firstVftVec+2*iFace) = 0;
1623  directions(0,firstVftVec+2*iFace+1) = 0;directions(1,firstVftVec+2*iFace+1) = 0;directions(2,firstVftVec+2*iFace+1) = 0;
1624  auto axes = TPZTetrahedron::TransformElementToSide(NCornerNodes+nEdges+iFace).Mult();
1625  axes.Transpose();
1626  for(auto x = 0; x < Dimension; x++){
1627  for(auto i = 0; i < 2; i++) {
1628  directions(x, firstVftVec + 2 * iFace) += axes(x,i) * vft1[i];
1629  directions(x, firstVftVec + 2 * iFace + 1) += axes(x,i) * vft2[i];
1630  }
1631  }
1632  }
1633  }
1634 
1635  int TPZTetrahedron::ClassId() const{
1636  return Hash("TPZTetrahedron");
1637  }
1638 
1639  void TPZTetrahedron::Read(TPZStream& buf, void* context) {
1640 
1641  }
1642 
1643  void TPZTetrahedron::Write(TPZStream& buf, int withclassid) const {
1644 
1645  }
1646 
1647 
1648 }
1649 
1650 /**********************************************************************************************************************
1651  * The following are explicit instantiation of member function template of this class, both with class T=REAL and its
1652  * respective FAD<REAL> version. In other to avoid potential errors, always declare the instantiation in the same order
1653  * in BOTH cases. @orlandini
1654  **********************************************************************************************************************/
1655 template bool pztopology::TPZTetrahedron::CheckProjectionForSingularity<REAL>(const int &side, const TPZVec<REAL> &xiInterior);
1656 
1657 template void pztopology::TPZTetrahedron::MapToSide<REAL>(int side, TPZVec<REAL> &InternalPar, TPZVec<REAL> &SidePar, TPZFMatrix<REAL> &JacToSide);
1658 
1659 template void pztopology::TPZTetrahedron::BlendFactorForSide<REAL>(const int &, const TPZVec<REAL> &, REAL &, TPZVec<REAL> &);
1660 
1661 template void pztopology::TPZTetrahedron::TShape<REAL>(const TPZVec<REAL> &loc,TPZFMatrix<REAL> &phi,TPZFMatrix<REAL> &dphi);
1662 
1663 template void pztopology::TPZTetrahedron::ComputeHDivDirections<REAL>(TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions);
1664 
1665 template void pztopology::TPZTetrahedron::ComputeHCurlDirections<REAL>(TPZFMatrix<REAL> &gradx, TPZFMatrix<REAL> &directions, const TPZVec<int> &transformationIds);
1666 #ifdef _AUTODIFF
1667 
1668 template bool pztopology::TPZTetrahedron::CheckProjectionForSingularity<Fad<REAL> >(const int &side, const TPZVec<Fad<REAL> > &xiInterior);
1669 
1670 template void pztopology::TPZTetrahedron::MapToSide<Fad<REAL> >(int side, TPZVec<Fad<REAL> > &InternalPar, TPZVec<Fad<REAL> > &SidePar, TPZFMatrix<Fad<REAL> > &JacToSide);
1671 
1672 template void pztopology::TPZTetrahedron::BlendFactorForSide<Fad<REAL>>(const int &, const TPZVec<Fad<REAL>> &, Fad<REAL> &,
1673  TPZVec<Fad<REAL>> &);
1674 template void pztopology::TPZTetrahedron::TShape<Fad<REAL>>(const TPZVec<Fad<REAL>> &loc,TPZFMatrix<Fad<REAL>> &phi,TPZFMatrix<Fad<REAL>> &dphi);
1675 
1676 template void pztopology::TPZTetrahedron::ComputeHDivDirections<Fad<REAL>>(TPZFMatrix<Fad<REAL>> &gradx, TPZFMatrix<Fad<REAL>> &directions);
1677 
1678 template void pztopology::TPZTetrahedron::ComputeHCurlDirections<Fad<REAL>>(TPZFMatrix<Fad<REAL>> &gradx, TPZFMatrix<Fad<REAL>> &directions, const TPZVec<int> &transformationIds);
1679 #endif
fabs
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
ETetraedro
Definition: pzeltype.h:58
pzquad.h
Contains the TPZInt1d, TPZIntTriang, TPZIntQuad, TPZIntCube3D, TPZIntTetra3D, TPZIntPyram3D and TPZIn...
pzlog.h
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.
TPZIntTriang
Handles the numerical integration for two-dimensional problems using triangular elements. Numerical Integration.
Definition: pzquad.h:102
ETriangle
Definition: pzeltype.h:56
pztopology::FaceConnectLocId
static int FaceConnectLocId[4][7]
Definition: tpztetrahedron.cpp:68
dimension
clarg::argInt dimension("-d", "Matrices dimension M x M", 1000)
ENoType
Definition: pzeltype.h:71
pztopology::sidedimension
static int sidedimension[15]
Definition: tpztetrahedron.cpp:65
pztopology::computedirectionsT3
void computedirectionsT3(int inicio, int fim, TPZFMatrix< REAL > &bvec, TPZFMatrix< REAL > &t1vec, TPZFMatrix< REAL > &t2vec, TPZFMatrix< REAL > &gradx, TPZFMatrix< REAL > &directions)
Definition: tpztetrahedron.cpp:1243
TPZTransform::Mult
const TPZFMatrix< T > & Mult() const
Definition: pztrnsform.h:64
TPZMatrix::Identity
virtual void Identity()
Converts the matrix in an identity matrix.
Definition: pzmatrix.cpp:199
pzerror.h
Defines PZError.
Fad
Definition: fad.h:54
pztopology::direcaoksioueta
static int direcaoksioueta[45]
Definition: tpztetrahedron.cpp:271
pzeltype.h
Defines enum MElementType and contains the implementation of MElementType_NNodes(...) functions.
TPZVec
This class implements a simple vector storage scheme for a templated class T. Utility.
Definition: pzgeopoint.h:19
val
REAL val(STATE &number)
Returns value of the variable.
Definition: pzartdiff.h:23
TPZAxesTools::ComputeDetjac
static TVar ComputeDetjac(TPZFMatrix< TVar > &gradx)
Definition: pzaxestools.h:111
TPZManVector::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object.
Definition: pzmanvector.h:426
tpztriangle.h
Contains the TPZTriangle class which defines the topology of a triangle.
TPZInt1d
Handles the numerical integration for one-dimensional problems. Numerical Integration.
Definition: pzquad.h:36
tpztetrahedron.h
Contains the TPZTetrahedron class which defines the topology of the tetrahedron element.
LOGPZ_FATAL
#define LOGPZ_FATAL(A, B)
Define log for fatal errors (cout)
Definition: pzlog.h:95
TPZIntPoints
Abstract class defining integration rules. Numerical Integration.
Definition: tpzintpoints.h:19
pztopology::highsides
static int highsides[15][7]
Definition: tpztetrahedron.cpp:71
TPZFMatrix::Zero
int Zero() override
Makes Zero all the elements.
Definition: pzfmatrix.h:651
pztopology::MidSideNode
static REAL MidSideNode[15][3]
Definition: tpztetrahedron.cpp:180
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
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
pztopology
Groups all classes defining the structure of the master element.
Definition: PrismExtend.cpp:15
TPZStack::Push
void Push(const T object)
Pushes a copy of the object on the stack.
Definition: pzstack.h:80
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
pztopology::t2Tetra
static REAL t2Tetra[45][3]
Definition: tpztetrahedron.cpp:225
pzmanvector.h
Free store vector implementation.
pztopology::TPZTriangle::GetTransformId
static int GetTransformId(TPZVec< int64_t > &id)
Method which identifies the transformation based on the IDs of the corner nodes.
Definition: tpztriangle.cpp:644
TPZTransform::Sum
const TPZFMatrix< T > & Sum() const
Definition: pztrnsform.h:66
sqrt
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ sqrt
Definition: tfadfunc.h:120
pztopology::t1Tetra
static REAL t1Tetra[45][3]
Definition: tpztetrahedron.cpp:205
pztopology::vectorsideorderTe
static int vectorsideorderTe[45]
Definition: tpztetrahedron.cpp:245
LOGPZ_ERROR
#define LOGPZ_ERROR(A, B)
Define log for errors (cout)
Definition: pzlog.h:93
TPZFMatrix::Redim
int Redim(const int64_t newRows, const int64_t newCols) override
Redimension a matrix and ZERO your elements.
Definition: pzfmatrix.h:616
pztopology::nsidenodes
static int nsidenodes[15]
Definition: tpztetrahedron.cpp:89
EPoint
Definition: pzeltype.h:54
Hash
int32_t Hash(std::string str)
Definition: TPZHash.cpp:10
TPZInt1Point
Integration rule for one point. Numerical Integration.
Definition: pzquad.h:478
pztopology::bilinearounao
static int bilinearounao[45]
Definition: tpztetrahedron.cpp:260
MElementType
MElementType
Define the element types.
Definition: pzeltype.h:52
pztopology::bTetra
static REAL bTetra[45][3]
Definition: tpztetrahedron.cpp:186
pztopology::sidetosidetransforms
static REAL sidetosidetransforms[15][7][4][3]
Definition: tpztetrahedron.cpp:96
TPZMatrix::Cols
int64_t Cols() const
Returns number of cols.
Definition: pzmatrix.h:809
pztopology::TPZTetrahedron::IsInParametricDomain
static bool IsInParametricDomain(const TPZVec< REAL > &pt, REAL tol=pztopology::gTolerance)
Verifies if the parametric point pt is in the element parametric domain.
Definition: tpztetrahedron.cpp:855
TPZFMatrix::Resize
int Resize(const int64_t newRows, const int64_t wCols) override
Redimension a matrix, but maintain your elements.
Definition: pzfmatrix.cpp:1016
TPZStream
Defines the interface for saving and reading data. Persistency.
Definition: TPZStream.h:50
pzreal.h
Contains the declaration of TPZFlopCounter class and TPZCounter struct.
TPZTransform
Implements an affine transformation between points in parameter space. Topology Utility.
Definition: pzmganalysis.h:14
EOned
Definition: pzeltype.h:55
TPZIntTetra3D
Handles the numerical integration for three-dimensional problems using tetraedra elements. Numerical Integration.
Definition: pzquad.h:322
pztopology::nhighdimsides
static int nhighdimsides[15]
Definition: tpztetrahedron.cpp:29
TPZFMatrix::PutVal
int PutVal(const int64_t row, const int64_t col, const TVar &value) override
Put values without bounds checking This method is faster than "Put" if DEBUG is defined.
Definition: pzfmatrix.h:548
TPZFMatrix::Transpose
void Transpose(TPZMatrix< TVar > *const T) const override
It makes *T the transpose of current matrix.
Definition: pzfmatrix.cpp:1073
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:11 for NeoPZ by   doxygen 1.8.13