gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
HierarchicalBasisH1Brick.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributed by Ismail Badia.
7 
8 // Reference : "Higher-Order Finite Element Methods"; Pavel Solin, Karel
9 // Segeth, Ivo Dolezel, Chapman and Hall/CRC; Edition : Har/Cdr (2003).
10 
11 #include <algorithm>
12 #include <stdexcept>
13 #include "HierarchicalBasisH1Brick.h"
14 
15 HierarchicalBasisH1Brick::HierarchicalBasisH1Brick(int order)
16 {
17  _pb1 = order;
18  _pb2 = order;
19  _pb3 = order;
20  for(int i = 0; i < 12; i++) { _pOrderEdge[i] = order; }
21  for(int j = 0; j < 6; j++) {
22  _pOrderFace1[j] = order;
23  _pOrderFace2[j] = order;
24  }
25  _nvertex = 8;
26  _nedge = 12;
27  _nfaceQuad = 6;
28  _nfaceTri = 0;
29  _nVertexFunction = 8;
30  _nEdgeFunction = 12 * order - 12;
31  _nQuadFaceFunction = 6 * (order - 1) * (order - 1);
32  _nTriFaceFunction = 0;
33  _nBubbleFunction = (order - 1) * (order - 1) * (order - 1);
34 }
35 
36 HierarchicalBasisH1Brick::~HierarchicalBasisH1Brick() {}
37 
38 unsigned int HierarchicalBasisH1Brick::getNumberOfOrientations() const
39 {
40  return 40320; // factorial 8
41 }
42 
43 double HierarchicalBasisH1Brick::_affineCoordinate(const int &j,
44  const double &u,
45  const double &v,
46  const double &w)
47 {
48  switch(j) {
49  case(1): return 0.5 * (1 + u);
50  case(2): return 0.5 * (1 - u);
51  case(3): return 0.5 * (1 + v);
52  case(4): return 0.5 * (1 - v);
53  case(5): return 0.5 * (1 + w);
54  case(6): return 0.5 * (1 - w);
55  default: throw std::runtime_error("j must be : 1<=j<=6");
56  }
57 }
58 inline void HierarchicalBasisH1Brick::_someProduct(double const &u,
59  double const &v,
60  double const &w,
61  std::vector<double> &product,
62  std::vector<double> &lambda)
63 {
64  lambda[0] = _affineCoordinate(1, u, v, w);
65  lambda[1] = _affineCoordinate(2, u, v, w);
66  lambda[2] = _affineCoordinate(3, u, v, w);
67  lambda[3] = _affineCoordinate(4, u, v, w);
68  lambda[4] = _affineCoordinate(5, u, v, w);
69  lambda[5] = _affineCoordinate(6, u, v, w);
70  product[0] = lambda[3] * lambda[5];
71  product[1] = lambda[1] * lambda[5];
72  product[2] = lambda[1] * lambda[3];
73  product[3] = lambda[0] * lambda[5];
74  product[4] = lambda[3] * lambda[0];
75  product[5] = lambda[2] * lambda[5];
76  product[6] = lambda[2] * lambda[0];
77  product[7] = lambda[2] * lambda[1];
78  product[8] = lambda[3] * lambda[4];
79  product[9] = lambda[4] * lambda[1];
80  product[10] = lambda[4] * lambda[0];
81  product[11] = lambda[4] * lambda[2];
82 }
83 void HierarchicalBasisH1Brick::generateBasis(double const &u, double const &v,
84  double const &w,
85  std::vector<double> &vertexBasis,
86  std::vector<double> &edgeBasis,
87  std::vector<double> &faceBasis,
88  std::vector<double> &bubbleBasis)
89 {
90  std::vector<double> product(12, 0);
91  std::vector<double> lambda(6, 0);
92  HierarchicalBasisH1Brick::_someProduct(u, v, w, product, lambda);
93  // vertex shape functions:
94  vertexBasis[0] = lambda[1] * product[0];
95  vertexBasis[1] = lambda[0] * product[0];
96  vertexBasis[2] = lambda[0] * product[5];
97  vertexBasis[3] = lambda[1] * product[5];
98  vertexBasis[4] = lambda[1] * product[8];
99  vertexBasis[5] = lambda[0] * product[8];
100  vertexBasis[6] = lambda[0] * product[11];
101  vertexBasis[7] = lambda[1] * product[11];
102  std::vector<double> lkVectorU(_pb1 - 1);
103  std::vector<double> lkVectorV(_pb2 - 1);
104  std::vector<double> lkVectorW(_pb3 - 1);
105  for(int it = 2; it <= _pb1; it++) {
106  lkVectorU[it - 2] = OrthogonalPoly::EvalLobatto(it, u);
107  }
108  for(int it = 2; it <= _pb2; it++) {
109  lkVectorV[it - 2] = OrthogonalPoly::EvalLobatto(it, v);
110  }
111  for(int it = 2; it <= _pb3; it++) {
112  lkVectorW[it - 2] = OrthogonalPoly::EvalLobatto(it, w);
113  }
114  // edge shape functions:
115  int indexEdgeBasis = 0;
116  std::vector<double> *vectorTarget1(nullptr);
117  for(int iEdge = 0; iEdge < _nedge; iEdge++) {
118  switch(iEdge) {
119  case(0):
120  case(5):
121  case(8):
122  case(11): vectorTarget1 = &lkVectorU; break;
123  case(1):
124  case(3):
125  case(9):
126  case(10): vectorTarget1 = &lkVectorV; break;
127  case(2):
128  case(4):
129  case(6):
130  case(7): vectorTarget1 = &lkVectorW; break;
131  }
132  for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] - 1;
133  indexEdgeFunc++) {
134  edgeBasis[indexEdgeBasis] =
135  (*vectorTarget1)[indexEdgeFunc] * product[iEdge];
136  indexEdgeBasis++;
137  }
138  }
139  // face shape functions:
140  int indexFaceFunction = 0;
141  std::vector<double> *vectorTarget2(nullptr);
142  for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
143  int indexLambda;
144  switch(iFace) {
145  case(0):
146  indexLambda = 5;
147  vectorTarget1 = &lkVectorU;
148  vectorTarget2 = &lkVectorV;
149  break;
150  case(1):
151  indexLambda = 3;
152  vectorTarget1 = &lkVectorU;
153  vectorTarget2 = &lkVectorW;
154  break;
155  case(2):
156  indexLambda = 1;
157  vectorTarget1 = &lkVectorV;
158  vectorTarget2 = &lkVectorW;
159  break;
160  case(3):
161  indexLambda = 0;
162  vectorTarget1 = &lkVectorV;
163  vectorTarget2 = &lkVectorW;
164  break;
165  case(4):
166  indexLambda = 2;
167  vectorTarget1 = &lkVectorU;
168  vectorTarget2 = &lkVectorW;
169  break;
170  case(5):
171  indexLambda = 4;
172  vectorTarget1 = &lkVectorU;
173  vectorTarget2 = &lkVectorV;
174  break;
175  }
176  for(int index1 = 0; index1 < _pOrderFace1[iFace] - 1; index1++) {
177  for(int index2 = 0; index2 < _pOrderFace2[iFace] - 1; index2++) {
178  faceBasis[indexFaceFunction] = lambda[indexLambda] *
179  (*vectorTarget1)[index1] *
180  (*vectorTarget2)[index2];
181  indexFaceFunction++;
182  }
183  }
184  }
185  // bubble shape functions:
186  int indexBubbleBasis = 0;
187  for(int ipb1 = 0; ipb1 < _pb1 - 1; ipb1++) {
188  for(int ipb2 = 0; ipb2 < _pb2 - 1; ipb2++) {
189  for(int ipb3 = 0; ipb3 < _pb3 - 1; ipb3++) {
190  bubbleBasis[indexBubbleBasis] =
191  lkVectorU[ipb1] * lkVectorV[ipb2] * lkVectorW[ipb3];
192  indexBubbleBasis++;
193  }
194  }
195  }
196 }
197 
198 void HierarchicalBasisH1Brick::_someProductGrad(
199  double const &u, double const &v, double const &w,
200  std::vector<double> &product,
201  std::vector<std::vector<double> > &gradientProduct,
202  std::vector<double> &lambda,
203  std::vector<std::vector<double> > &gradientLambda)
204 {
205  lambda[0] = _affineCoordinate(1, u, v, w);
206  lambda[1] = _affineCoordinate(2, u, v, w);
207  lambda[2] = _affineCoordinate(3, u, v, w);
208  lambda[3] = _affineCoordinate(4, u, v, w);
209  lambda[4] = _affineCoordinate(5, u, v, w);
210  lambda[5] = _affineCoordinate(6, u, v, w);
211  gradientLambda[0][0] = 0.5;
212  gradientLambda[1][0] = -0.5;
213  gradientLambda[2][1] = 0.5;
214  gradientLambda[3][1] = -0.5;
215  gradientLambda[4][2] = 0.5;
216  gradientLambda[5][2] = -0.5;
217  product[0] = lambda[3] * lambda[5];
218  product[1] = lambda[1] * lambda[5];
219  product[2] = lambda[1] * lambda[3];
220  product[3] = lambda[0] * lambda[5];
221  product[4] = lambda[3] * lambda[0];
222  product[5] = lambda[2] * lambda[5];
223  product[6] = lambda[2] * lambda[0];
224  product[7] = lambda[2] * lambda[1];
225  product[8] = lambda[3] * lambda[4];
226  product[9] = lambda[4] * lambda[1];
227  product[10] = lambda[4] * lambda[0];
228  product[11] = lambda[4] * lambda[2];
229  gradientProduct[0][1] = -0.5 * lambda[5];
230  gradientProduct[0][2] = -0.5 * lambda[3];
231  gradientProduct[1][0] = -0.5 * lambda[5];
232  gradientProduct[1][2] = -0.5 * lambda[1];
233  gradientProduct[2][0] = -0.5 * lambda[3];
234  gradientProduct[2][1] = -0.5 * lambda[1];
235  gradientProduct[3][0] = 0.5 * lambda[5];
236  gradientProduct[3][2] = -0.5 * lambda[0];
237  gradientProduct[4][0] = 0.5 * lambda[3];
238  gradientProduct[4][1] = -0.5 * lambda[0];
239  gradientProduct[5][1] = 0.5 * lambda[5];
240  gradientProduct[5][2] = -0.5 * lambda[2];
241  gradientProduct[6][0] = 0.5 * lambda[2];
242  gradientProduct[6][1] = 0.5 * lambda[0];
243  gradientProduct[7][0] = -0.5 * lambda[2];
244  gradientProduct[7][1] = 0.5 * lambda[1];
245  gradientProduct[8][1] = -0.5 * lambda[4];
246  gradientProduct[8][2] = 0.5 * lambda[3];
247  gradientProduct[9][0] = -0.5 * lambda[4];
248  gradientProduct[9][2] = 0.5 * lambda[1];
249  gradientProduct[10][0] = 0.5 * lambda[4];
250  gradientProduct[10][2] = 0.5 * lambda[0];
251  gradientProduct[11][1] = 0.5 * lambda[4];
252  gradientProduct[11][2] = 0.5 * lambda[2];
253 }
254 
255 void HierarchicalBasisH1Brick::generateGradientBasis(
256  double const &u, double const &v, double const &w,
257  std::vector<std::vector<double> > &gradientVertex,
258  std::vector<std::vector<double> > &gradientEdge,
259  std::vector<std::vector<double> > &gradientFace,
260  std::vector<std::vector<double> > &gradientBubble)
261 {
262  std::vector<double> product(12, 0);
263  std::vector<std::vector<double> > gradientProduct(12,
264  std::vector<double>(3, 0));
265  std::vector<double> lambda(6, 0);
266  std::vector<std::vector<double> > gradientLambda(6,
267  std::vector<double>(3, 0));
268  HierarchicalBasisH1Brick::_someProductGrad(u, v, w, product, gradientProduct,
269  lambda, gradientLambda);
270  // vertex gradient:
271  for(int it = 0; it < 3; it++) {
272  gradientVertex[0][it] =
273  gradientLambda[1][it] * product[0] + gradientProduct[0][it] * lambda[1];
274  gradientVertex[1][it] =
275  gradientLambda[0][it] * product[0] + gradientProduct[0][it] * lambda[0];
276  gradientVertex[2][it] =
277  gradientLambda[0][it] * product[5] + gradientProduct[5][it] * lambda[0];
278  gradientVertex[3][it] =
279  gradientLambda[1][it] * product[5] + gradientProduct[5][it] * lambda[1];
280  gradientVertex[4][it] =
281  gradientLambda[1][it] * product[8] + gradientProduct[8][it] * lambda[1];
282  gradientVertex[5][it] =
283  gradientLambda[0][it] * product[8] + gradientProduct[8][it] * lambda[0];
284  gradientVertex[6][it] =
285  gradientLambda[0][it] * product[11] + gradientProduct[11][it] * lambda[0];
286  gradientVertex[7][it] =
287  gradientLambda[1][it] * product[11] + gradientProduct[11][it] * lambda[1];
288  }
289  std::vector<double> lkVectorU(_pb1 - 1);
290  std::vector<double> lkVectorV(_pb2 - 1);
291  std::vector<double> lkVectorW(_pb3 - 1);
292  std::vector<std::vector<double> > dlkVectorU(_pb1 - 1,
293  std::vector<double>(3, 0.));
294  std::vector<std::vector<double> > dlkVectorV(_pb2 - 1,
295  std::vector<double>(3, 0.));
296  std::vector<std::vector<double> > dlkVectorW(_pb3 - 1,
297  std::vector<double>(3, 0.));
298  for(int it = 2; it <= _pb1; it++) {
299  lkVectorU[it - 2] = OrthogonalPoly::EvalLobatto(it, u);
300  dlkVectorU[it - 2][0] = OrthogonalPoly::EvalDLobatto(it, u);
301  }
302  for(int it = 2; it <= _pb2; it++) {
303  lkVectorV[it - 2] = OrthogonalPoly::EvalLobatto(it, v);
304  dlkVectorV[it - 2][1] = OrthogonalPoly::EvalDLobatto(it, v);
305  }
306  for(int it = 2; it <= _pb3; it++) {
307  lkVectorW[it - 2] = OrthogonalPoly::EvalLobatto(it, w);
308  dlkVectorW[it - 2][2] = OrthogonalPoly::EvalDLobatto(it, w);
309  }
310  // edge gradient:
311  int indexEdgeBasis = 0;
312  std::vector<double> *vectorTarget1(nullptr);
313  std::vector<std::vector<double> > *dvectorTarget1(nullptr);
314  for(int iEdge = 0; iEdge < _nedge; iEdge++) {
315  switch(iEdge) {
316  case(0):
317  case(5):
318  case(8):
319  case(11):
320  vectorTarget1 = &lkVectorU;
321  dvectorTarget1 = &dlkVectorU;
322  break;
323  case(1):
324  case(3):
325  case(9):
326  case(10):
327  vectorTarget1 = &lkVectorV;
328  dvectorTarget1 = &dlkVectorV;
329  break;
330  case(2):
331  case(4):
332  case(6):
333  case(7):
334  vectorTarget1 = &lkVectorW;
335  dvectorTarget1 = &dlkVectorW;
336  break;
337  }
338  for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] - 1;
339  indexEdgeFunc++) {
340  for(int it = 0; it < 3; it++) {
341  gradientEdge[indexEdgeBasis][it] =
342  (*dvectorTarget1)[indexEdgeFunc][it] * product[iEdge] +
343  (*vectorTarget1)[indexEdgeFunc] * gradientProduct[iEdge][it];
344  }
345  indexEdgeBasis++;
346  }
347  }
348  // face gradient:
349  int indexFaceFunction = 0;
350  std::vector<double> *vectorTarget2(nullptr);
351  std::vector<std::vector<double> > *dvectorTarget2(nullptr);
352  for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
353  int indexLambda = 0;
354  switch(iFace) {
355  case(0):
356  indexLambda = 5;
357  vectorTarget1 = &lkVectorU;
358  vectorTarget2 = &lkVectorV;
359  dvectorTarget1 = &dlkVectorU;
360  dvectorTarget2 = &dlkVectorV;
361  break;
362  case(1):
363  indexLambda = 3;
364  vectorTarget1 = &lkVectorU;
365  vectorTarget2 = &lkVectorW;
366  dvectorTarget1 = &dlkVectorU;
367  dvectorTarget2 = &dlkVectorW;
368  break;
369  case(2):
370  indexLambda = 1;
371  vectorTarget1 = &lkVectorV;
372  vectorTarget2 = &lkVectorW;
373  dvectorTarget1 = &dlkVectorV;
374  dvectorTarget2 = &dlkVectorW;
375  break;
376  case(3):
377  indexLambda = 0;
378  vectorTarget1 = &lkVectorV;
379  vectorTarget2 = &lkVectorW;
380  dvectorTarget1 = &dlkVectorV;
381  dvectorTarget2 = &dlkVectorW;
382  break;
383  case(4):
384  indexLambda = 2;
385  vectorTarget1 = &lkVectorU;
386  vectorTarget2 = &lkVectorW;
387  dvectorTarget1 = &dlkVectorU;
388  dvectorTarget2 = &dlkVectorW;
389  break;
390  case(5):
391  indexLambda = 4;
392  vectorTarget1 = &lkVectorU;
393  vectorTarget2 = &lkVectorV;
394  dvectorTarget1 = &dlkVectorU;
395  dvectorTarget2 = &dlkVectorV;
396  break;
397  }
398  for(int index1 = 0; index1 < _pOrderFace1[iFace] - 1; index1++) {
399  for(int index2 = 0; index2 < _pOrderFace2[iFace] - 1; index2++) {
400  for(int it = 0; it < 3; it++) {
401  gradientFace[indexFaceFunction][it] =
402  gradientLambda[indexLambda][it] * (*vectorTarget1)[index1] *
403  (*vectorTarget2)[index2] +
404  lambda[indexLambda] * (*dvectorTarget1)[index1][it] *
405  (*vectorTarget2)[index2] +
406  lambda[indexLambda] * (*vectorTarget1)[index1] *
407  (*dvectorTarget2)[index2][it];
408  }
409  indexFaceFunction++;
410  }
411  }
412  }
413  // bubble shape functions:
414  int indexBubbleBasis = 0;
415  for(int ipb1 = 0; ipb1 < _pb1 - 1; ipb1++) {
416  for(int ipb2 = 0; ipb2 < _pb2 - 1; ipb2++) {
417  for(int ipb3 = 0; ipb3 < _pb3 - 1; ipb3++) {
418  gradientBubble[indexBubbleBasis][0] =
419  dlkVectorU[ipb1][0] * lkVectorV[ipb2] * lkVectorW[ipb3];
420  gradientBubble[indexBubbleBasis][1] =
421  lkVectorU[ipb1] * dlkVectorV[ipb2][1] * lkVectorW[ipb3];
422  gradientBubble[indexBubbleBasis][2] =
423  lkVectorU[ipb1] * lkVectorV[ipb2] * dlkVectorW[ipb3][2];
424  indexBubbleBasis++;
425  }
426  }
427  }
428 }
429 
430 void HierarchicalBasisH1Brick::orientEdge(
431  int const &flagOrientation, int const &edgeNumber,
432  std::vector<double> &edgeFunctions,
433  const std::vector<double> &eTablePositiveFlag,
434  const std::vector<double> &eTableNegativeFlag)
435 {
436  if(flagOrientation == -1) {
437  int constant1 = 0;
438  int constant2 = 0;
439  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
440  constant2 = constant2 - 1;
441  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
442  for(int k = constant1; k <= constant2; k++) {
443  edgeFunctions[k] = eTableNegativeFlag[k];
444  }
445  }
446  else {
447  int constant1 = 0;
448  int constant2 = 0;
449  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
450  constant2 = constant2 - 1;
451  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
452  for(int k = constant1; k <= constant2; k++) {
453  edgeFunctions[k] = eTablePositiveFlag[k];
454  }
455  }
456 }
457 
458 void HierarchicalBasisH1Brick::orientEdge(
459  int const &flagOrientation, int const &edgeNumber,
460  std::vector<std::vector<double> > &edgeFunctions,
461  const std::vector<std::vector<double> > &eTablePositiveFlag,
462  const std::vector<std::vector<double> > &eTableNegativeFlag)
463 {
464  if(flagOrientation == -1) {
465  int constant1 = 0;
466  int constant2 = 0;
467  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
468  constant2 = constant2 - 1;
469  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
470  for(int k = constant1; k <= constant2; k++) {
471  edgeFunctions[k][0] = eTableNegativeFlag[k][0];
472  edgeFunctions[k][1] = eTableNegativeFlag[k][1];
473  edgeFunctions[k][2] = eTableNegativeFlag[k][2];
474  }
475  }
476  else {
477  int constant1 = 0;
478  int constant2 = 0;
479  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
480  constant2 = constant2 - 1;
481  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
482  for(int k = constant1; k <= constant2; k++) {
483  edgeFunctions[k][0] = eTablePositiveFlag[k][0];
484  edgeFunctions[k][1] = eTablePositiveFlag[k][1];
485  edgeFunctions[k][2] = eTablePositiveFlag[k][2];
486  }
487  }
488 }
489 
490 void HierarchicalBasisH1Brick::orientEdgeFunctionsForNegativeFlag(
491  std::vector<double> &edgeFunctions)
492 {
493  int constant1 = 0;
494  int constant2 = 0;
495  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
496  constant2 = 0;
497  constant2 = 0;
498  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
499  constant2 = constant2 - 1;
500  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
501  for(int k = constant1; k <= constant2; k++) {
502  if((k - constant1) % 2 != 0) {
503  edgeFunctions[k] = edgeFunctions[k] * (-1);
504  }
505  }
506  }
507 }
508 
509 void HierarchicalBasisH1Brick::orientEdgeFunctionsForNegativeFlag(
510  std::vector<std::vector<double> > &edgeFunctions)
511 {
512  int constant1 = 0;
513  int constant2 = 0;
514  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
515  constant2 = 0;
516  constant2 = 0;
517  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
518  constant2 = constant2 - 1;
519  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
520  for(int k = constant1; k <= constant2; k++) {
521  if((k - constant1) % 2 != 0) {
522  edgeFunctions[k][0] = edgeFunctions[k][0] * (-1);
523  edgeFunctions[k][1] = edgeFunctions[k][1] * (-1);
524  edgeFunctions[k][2] = edgeFunctions[k][2] * (-1);
525  }
526  }
527  }
528 }
529 
530 void HierarchicalBasisH1Brick::orientOneFace(double const &u, double const &v,
531  double const &w, int const &flag1,
532  int const &flag2, int const &flag3,
533  int const &faceNumber,
534  std::vector<double> &faceBasis)
535 {
536  if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
537  int iterator = 0;
538  for(int i = 0; i < faceNumber; i++) {
539  iterator += (_pOrderFace1[i] - 1) * (_pOrderFace2[i] - 1);
540  }
541  if(flag3 == 1) {
542  for(int it1 = 2; it1 <= _pOrderFace1[faceNumber]; it1++) {
543  for(int it2 = 2; it2 <= _pOrderFace2[faceNumber]; it2++) {
544  int impactFlag1 = 1;
545  int impactFlag2 = 1;
546  if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
547  if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
548  faceBasis[iterator] = faceBasis[iterator] * impactFlag1 * impactFlag2;
549  iterator++;
550  }
551  }
552  }
553  else {
554  double lambda = 0;
555  double var1 = 0;
556  double var2 = 0;
557  switch(faceNumber) {
558  case(0):
559  lambda = _affineCoordinate(6, u, v, w);
560  var1 = u;
561  var2 = v;
562  break;
563  case(1):
564  lambda = _affineCoordinate(4, u, v, w);
565  var1 = u;
566  var2 = w;
567  break;
568  case(2):
569  lambda = _affineCoordinate(2, u, v, w);
570  var1 = v;
571  var2 = w;
572  break;
573  case(3):
574  lambda = _affineCoordinate(1, u, v, w);
575  var1 = v;
576  var2 = w;
577  break;
578  case(4):
579  lambda = _affineCoordinate(3, u, v, w);
580  var1 = u;
581  var2 = w;
582  break;
583  case(5):
584  lambda = _affineCoordinate(5, u, v, w);
585  var1 = u;
586  var2 = v;
587  break;
588  }
589  std::vector<double> lkVector1(_pOrderFace1[faceNumber] - 1);
590  std::vector<double> lkVector2(_pOrderFace2[faceNumber] - 1);
591  for(int it = 2; it <= _pOrderFace1[faceNumber]; it++) {
592  lkVector1[it - 2] = OrthogonalPoly::EvalLobatto(it, var1);
593  }
594  for(int it = 2; it <= _pOrderFace2[faceNumber]; it++) {
595  lkVector2[it - 2] = OrthogonalPoly::EvalLobatto(it, var2);
596  }
597 
598  for(int it1 = 2; it1 <= _pOrderFace2[faceNumber]; it1++) {
599  for(int it2 = 2; it2 <= _pOrderFace1[faceNumber]; it2++) {
600  int impactFlag1 = 1;
601  int impactFlag2 = 1;
602  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
603  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
604  faceBasis[iterator] = lambda * lkVector1[it2 - 2] *
605  lkVector2[it1 - 2] * impactFlag1 * impactFlag2;
606  iterator++;
607  }
608  }
609  }
610  }
611 }
612 void HierarchicalBasisH1Brick::orientOneFace(
613  double const &u, double const &v, double const &w, int const &flag1,
614  int const &flag2, int const &flag3, int const &faceNumber,
615  std::vector<std::vector<double> > &gradientFace, std::string typeFunction)
616 {
617  if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
618  int iterator = 0;
619  for(int i = 0; i < faceNumber; i++) {
620  iterator += (_pOrderFace1[i] - 1) * (_pOrderFace2[i] - 1);
621  }
622  if(flag3 == 1) {
623  for(int it1 = 2; it1 <= _pOrderFace1[faceNumber]; it1++) {
624  for(int it2 = 2; it2 <= _pOrderFace2[faceNumber]; it2++) {
625  int impactFlag1 = 1;
626  int impactFlag2 = 1;
627  if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
628  if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
629  gradientFace[iterator][0] =
630  gradientFace[iterator][0] * impactFlag1 * impactFlag2;
631  gradientFace[iterator][1] =
632  gradientFace[iterator][1] * impactFlag1 * impactFlag2;
633  gradientFace[iterator][2] =
634  gradientFace[iterator][2] * impactFlag1 * impactFlag2;
635  iterator++;
636  }
637  }
638  }
639  else {
640  std::vector<double> uvw(3);
641  uvw[0] = u;
642  uvw[1] = v;
643  uvw[2] = w;
644  double lambda = 0;
645  int var1 = 0;
646  int var2 = 0;
647  std::vector<double> gradientLambda(3, 0);
648  switch(faceNumber) {
649  case(0):
650  lambda = _affineCoordinate(6, u, v, w);
651  var1 = 0;
652  var2 = 1;
653  gradientLambda[2] = -0.5;
654  break;
655  case(1):
656  lambda = _affineCoordinate(4, u, v, w);
657  var1 = 0;
658  var2 = 2;
659  gradientLambda[1] = -0.5;
660  break;
661  case(2):
662  lambda = _affineCoordinate(2, u, v, w);
663  var1 = 1;
664  var2 = 2;
665  gradientLambda[0] = -0.5;
666  break;
667  case(3):
668  lambda = _affineCoordinate(1, u, v, w);
669  var1 = 1;
670  var2 = 2;
671  gradientLambda[0] = 0.5;
672  break;
673  case(4):
674  lambda = _affineCoordinate(3, u, v, w);
675  var1 = 0;
676  var2 = 2;
677  gradientLambda[1] = 0.5;
678  break;
679  case(5):
680  lambda = _affineCoordinate(5, u, v, w);
681  var1 = 0;
682  var2 = 1;
683  gradientLambda[2] = 0.5;
684  break;
685  }
686  std::vector<double> lkVector1(_pOrderFace1[faceNumber] - 1);
687  std::vector<double> lkVector2(_pOrderFace2[faceNumber] - 1);
688  std::vector<std::vector<double> > dlkVector1(_pOrderFace1[faceNumber] - 1,
689  std::vector<double>(3, 0.));
690  std::vector<std::vector<double> > dlkVector2(_pOrderFace2[faceNumber] - 1,
691  std::vector<double>(3, 0.));
692  for(int it = 2; it <= _pOrderFace1[faceNumber]; it++) {
693  lkVector1[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[var1]);
694  dlkVector1[it - 2][var1] = OrthogonalPoly::EvalDLobatto(it, uvw[var1]);
695  }
696  for(int it = 2; it <= _pOrderFace2[faceNumber]; it++) {
697  lkVector2[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[var2]);
698  dlkVector2[it - 2][var2] = OrthogonalPoly::EvalDLobatto(it, uvw[var2]);
699  }
700  for(int it1 = 2; it1 <= _pOrderFace2[faceNumber]; it1++) {
701  for(int it2 = 2; it2 <= _pOrderFace1[faceNumber]; it2++) {
702  int impactFlag1 = 1;
703  int impactFlag2 = 1;
704  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
705  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
706  for(int itVector = 0; itVector < 3; itVector++) {
707  gradientFace[iterator][itVector] =
708  (gradientLambda[itVector] * lkVector1[it2 - 2] *
709  lkVector2[it1 - 2] +
710  lambda * dlkVector1[it2 - 2][itVector] * lkVector2[it1 - 2] +
711  lambda * lkVector1[it2 - 2] * dlkVector2[it1 - 2][itVector]) *
712  impactFlag1 * impactFlag2;
713  }
714  iterator++;
715  }
716  }
717  }
718  }
719 }
720 
721 void HierarchicalBasisH1Brick::orientFace(
722  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
723  const std::vector<double> &quadFaceFunctionsAllOrientation,
724  const std::vector<double> &triFaceFunctionsAllOrientation,
725  std::vector<double> &fTableCopy)
726 {
727  int iterator = 0;
728  for(int i = 0; i < faceNumber; i++) {
729  iterator += (_pOrderFace1[i] - 1) * (_pOrderFace2[i] - 1);
730  }
731  int numFaceFunctions =
732  (_pOrderFace1[faceNumber] - 1) * (_pOrderFace2[faceNumber] - 1);
733  int iOrientation = numberOrientationQuadFace(flag1, flag2, flag3);
734  int offset = iOrientation * _nQuadFaceFunction;
735  int offset2 = iterator + numFaceFunctions;
736  for(int i = iterator; i < offset2; i++) {
737  fTableCopy[i] = quadFaceFunctionsAllOrientation[i + offset];
738  }
739 }
740 void HierarchicalBasisH1Brick::orientFace(
741  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
742  const std::vector<std::vector<double> > &quadFaceFunctionsAllOrientation,
743  const std::vector<std::vector<double> > &triFaceFunctionsAllOrientation,
744  std::vector<std::vector<double> > &fTableCopy)
745 {
746  int iterator = 0;
747  for(int i = 0; i < faceNumber; i++) {
748  iterator += (_pOrderFace1[i] - 1) * (_pOrderFace2[i] - 1);
749  }
750  int numFaceFunctions =
751  (_pOrderFace1[faceNumber] - 1) * (_pOrderFace2[faceNumber] - 1);
752  int iOrientation = numberOrientationQuadFace(flag1, flag2, flag3);
753  int offset = iOrientation * _nQuadFaceFunction;
754  int offset2 = iterator + numFaceFunctions;
755  for(int i = iterator; i < offset2; i++) {
756  fTableCopy[i][0] = quadFaceFunctionsAllOrientation[i + offset][0];
757  fTableCopy[i][1] = quadFaceFunctionsAllOrientation[i + offset][1];
758  fTableCopy[i][2] = quadFaceFunctionsAllOrientation[i + offset][2];
759  }
760 }
761 
762 void HierarchicalBasisH1Brick::getKeysInfo(std::vector<int> &functionTypeInfo,
763  std::vector<int> &orderInfo)
764 {
765  for(int i = 0; i < 8; i++) {
766  functionTypeInfo[i] = 0;
767  orderInfo[i] = 1;
768  }
769  int it = 8;
770  for(int numEdge = 0; numEdge < 12; numEdge++) {
771  for(int i = 2; i <= _pOrderEdge[numEdge]; i++) {
772  functionTypeInfo[it] = 1;
773  orderInfo[it] = i;
774  it++;
775  }
776  }
777  for(int numFace = 0; numFace < 6; numFace++) {
778  for(int n1 = 2; n1 <= _pOrderFace1[numFace]; n1++) {
779  for(int n2 = 2; n2 <= _pOrderFace2[numFace]; n2++) {
780  functionTypeInfo[it] = 2;
781  orderInfo[it] = std::max(n1, n2);
782  it++;
783  }
784  }
785  }
786  for(int ipb1 = 2; ipb1 <= _pb1; ipb1++) {
787  for(int ipb2 = 2; ipb2 <= _pb2; ipb2++) {
788  for(int ipb3 = 2; ipb3 <= _pb3; ipb3++) {
789  functionTypeInfo[it] = 3;
790  orderInfo[it] = std::max(std::max(ipb1, ipb2), ipb3);
791  it++;
792  }
793  }
794  }
795 }
HierarchicalBasisH1Brick::orientEdge
virtual void orientEdge(int const &flagOrientation, int const &edgeNumber, std::vector< double > &edgeFunctions, const std::vector< double > &eTablePositiveFlag, const std::vector< double > &eTableNegativeFlag)
Definition: HierarchicalBasisH1Brick.cpp:430
HierarchicalBasis::_nTriFaceFunction
int _nTriFaceFunction
Definition: HierarchicalBasis.h:27
HierarchicalBasisH1Brick::getKeysInfo
virtual void getKeysInfo(std::vector< int > &functionTypeInfo, std::vector< int > &orderInfo)
Definition: HierarchicalBasisH1Brick.cpp:762
HierarchicalBasisH1Brick::orientFace
virtual void orientFace(int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, const std::vector< double > &quadFaceFunctionsAllOrientation, const std::vector< double > &triFaceFunctionsAllOrientation, std::vector< double > &fTableCopy)
Definition: HierarchicalBasisH1Brick.cpp:721
HierarchicalBasisH1Brick::HierarchicalBasisH1Brick
HierarchicalBasisH1Brick(int order)
Definition: HierarchicalBasisH1Brick.cpp:15
HierarchicalBasis::_nvertex
int _nvertex
Definition: HierarchicalBasis.h:20
HierarchicalBasis::_nBubbleFunction
int _nBubbleFunction
Definition: HierarchicalBasis.h:28
HierarchicalBasisH1Brick::_pb3
int _pb3
Definition: HierarchicalBasisH1Brick.h:103
HierarchicalBasisH1Brick::_someProduct
void _someProduct(double const &u, double const &v, double const &w, std::vector< double > &product, std::vector< double > &lambda)
Definition: HierarchicalBasisH1Brick.cpp:58
HierarchicalBasisH1Brick::generateGradientBasis
void generateGradientBasis(double const &u, double const &v, double const &w, std::vector< std::vector< double > > &gradientVertex, std::vector< std::vector< double > > &gradientEdge, std::vector< std::vector< double > > &gradientFace, std::vector< std::vector< double > > &gradientBubble)
Definition: HierarchicalBasisH1Brick.cpp:255
HierarchicalBasisH1Brick.h
HierarchicalBasisH1Brick::_pOrderFace1
int _pOrderFace1[6]
Definition: HierarchicalBasisH1Brick.h:106
OrthogonalPoly::EvalLobatto
double EvalLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:9
OrthogonalPoly::EvalDLobatto
double EvalDLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:127
HierarchicalBasisH1Brick::_pOrderEdge
int _pOrderEdge[12]
Definition: HierarchicalBasisH1Brick.h:104
HierarchicalBasisH1Brick::_pb2
int _pb2
Definition: HierarchicalBasisH1Brick.h:102
HierarchicalBasisH1Brick::getNumberOfOrientations
virtual unsigned int getNumberOfOrientations() const
Definition: HierarchicalBasisH1Brick.cpp:38
HierarchicalBasisH1Brick::_someProductGrad
void _someProductGrad(double const &u, double const &v, double const &w, std::vector< double > &product, std::vector< std::vector< double > > &gradientProduct, std::vector< double > &lambda, std::vector< std::vector< double > > &gradientLambda)
Definition: HierarchicalBasisH1Brick.cpp:198
HierarchicalBasisH1Brick::_affineCoordinate
static double _affineCoordinate(const int &j, const double &u, const double &v, const double &w)
Definition: HierarchicalBasisH1Brick.cpp:43
HierarchicalBasisH1Brick::generateBasis
virtual void generateBasis(double const &u, double const &v, double const &w, std::vector< double > &vertexBasis, std::vector< double > &edgeBasis, std::vector< double > &faceBasis, std::vector< double > &bubbleBasis)
Definition: HierarchicalBasisH1Brick.cpp:83
HierarchicalBasisH1Brick::_pOrderFace2
int _pOrderFace2[6]
Definition: HierarchicalBasisH1Brick.h:108
HierarchicalBasisH1Brick::orientEdgeFunctionsForNegativeFlag
virtual void orientEdgeFunctionsForNegativeFlag(std::vector< double > &edgeFunctions)
Definition: HierarchicalBasisH1Brick.cpp:490
HierarchicalBasisH1Brick::_pb1
int _pb1
Definition: HierarchicalBasisH1Brick.h:101
HierarchicalBasis::_nQuadFaceFunction
int _nQuadFaceFunction
Definition: HierarchicalBasis.h:26
HierarchicalBasis::_nfaceTri
int _nfaceTri
Definition: HierarchicalBasis.h:23
HierarchicalBasis::_nfaceQuad
int _nfaceQuad
Definition: HierarchicalBasis.h:22
HierarchicalBasis::_nVertexFunction
int _nVertexFunction
Definition: HierarchicalBasis.h:24
HierarchicalBasisH1Brick::orientOneFace
virtual void orientOneFace(double const &u, double const &v, double const &w, int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, std::vector< double > &faceBasis)
Definition: HierarchicalBasisH1Brick.cpp:530
HierarchicalBasis::_nEdgeFunction
int _nEdgeFunction
Definition: HierarchicalBasis.h:25
HierarchicalBasisH1Brick::~HierarchicalBasisH1Brick
virtual ~HierarchicalBasisH1Brick()
Definition: HierarchicalBasisH1Brick.cpp:36
HierarchicalBasis::_nedge
int _nedge
Definition: HierarchicalBasis.h:21
HierarchicalBasis::numberOrientationQuadFace
int numberOrientationQuadFace(int const &flag1, int const &flag2, int const &flag3)
Definition: HierarchicalBasis.h:105