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