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