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