gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
HierarchicalBasisHcurlBrick.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 #include <stdexcept>
9 #include <algorithm>
10 #include "HierarchicalBasisHcurlBrick.h"
11 
12 HierarchicalBasisHcurlBrick::HierarchicalBasisHcurlBrick(int order)
13 {
14  _pb1 = order;
15  _pb2 = order;
16  _pb3 = order;
17  for(int i = 0; i < 12; i++) { _pOrderEdge[i] = order; }
18  for(int j = 0; j < 6; j++) {
19  _pOrderFace1[j] = order;
20  _pOrderFace2[j] = order;
21  }
22  _nvertex = 8;
23  _nedge = 12;
24  _nfaceQuad = 6;
25  _nfaceTri = 0;
26  _nVertexFunction = 0;
27  _nEdgeFunction = 12 * order + 12;
28  _nQuadFaceFunction = 12 * order * (order + 1);
29  _nTriFaceFunction = 0;
30  _nBubbleFunction = 3 * order * order * (order + 1);
31 }
32 
33 HierarchicalBasisHcurlBrick::~HierarchicalBasisHcurlBrick() {}
34 
35 unsigned int HierarchicalBasisHcurlBrick::getNumberOfOrientations() const
36 {
37  return 40320; // factorial 8
38 }
39 
40 double HierarchicalBasisHcurlBrick::_affineCoordinate(const int &j,
41  const double &u,
42  const double &v,
43  const double &w)
44 {
45  switch(j) {
46  case(1): return 0.5 * (1 + u);
47  case(2): return 0.5 * (1 - u);
48  case(3): return 0.5 * (1 + v);
49  case(4): return 0.5 * (1 - v);
50  case(5): return 0.5 * (1 + w);
51  case(6): return 0.5 * (1 - w);
52  default: throw std::runtime_error("j must be : 1<=j<=6");
53  }
54 }
55 
56 void HierarchicalBasisHcurlBrick::generateHcurlBasis(
57  double const &u, double const &v, double const &w,
58  std::vector<std::vector<double> > &edgeBasis,
59  std::vector<std::vector<double> > &faceBasis,
60  std::vector<std::vector<double> > &bubbleBasis)
61 {
62  std::vector<std::vector<double> > lobattoVector(3);
63  lobattoVector[0] = std::vector<double>(_pb1);
64  lobattoVector[1] = std::vector<double>(_pb2);
65  lobattoVector[2] = std::vector<double>(_pb3);
66  for(int it = 2; it <= _pb1 + 1; it++) {
67  lobattoVector[0][it - 2] = OrthogonalPoly::EvalLobatto(it, u);
68  }
69  for(int it = 2; it <= _pb2 + 1; it++) {
70  lobattoVector[1][it - 2] = OrthogonalPoly::EvalLobatto(it, v);
71  }
72  for(int it = 2; it <= _pb3 + 1; it++) {
73  lobattoVector[2][it - 2] = OrthogonalPoly::EvalLobatto(it, w);
74  }
75  std::vector<std::vector<double> > legendreVector(3);
76  legendreVector[0] = std::vector<double>(_pb1 + 1);
77  legendreVector[1] = std::vector<double>(_pb2 + 1);
78  legendreVector[2] = std::vector<double>(_pb3 + 1);
79  for(unsigned int k = 0; k < legendreVector[0].size(); k++) {
80  legendreVector[0][k] = OrthogonalPoly::EvalLegendre(k, u);
81  }
82  for(unsigned int k = 0; k < legendreVector[1].size(); k++) {
83  legendreVector[1][k] = OrthogonalPoly::EvalLegendre(k, v);
84  }
85  for(unsigned int k = 0; k < legendreVector[2].size(); k++) {
86  legendreVector[2][k] = OrthogonalPoly::EvalLegendre(k, w);
87  }
88  std::vector<double> lambda(6, 0);
89  lambda[0] = _affineCoordinate(1, u, v, w);
90  lambda[1] = _affineCoordinate(2, u, v, w);
91  lambda[2] = _affineCoordinate(3, u, v, w);
92  lambda[3] = _affineCoordinate(4, u, v, w);
93  lambda[4] = _affineCoordinate(5, u, v, w);
94  lambda[5] = _affineCoordinate(6, u, v, w);
95  std::vector<double> product(12, 0);
96  product[0] = lambda[3] * lambda[5];
97  product[1] = lambda[1] * lambda[5];
98  product[2] = lambda[1] * lambda[3];
99  product[3] = lambda[0] * lambda[5];
100  product[4] = lambda[3] * lambda[0];
101  product[5] = lambda[2] * lambda[5];
102  product[6] = lambda[2] * lambda[0];
103  product[7] = lambda[2] * lambda[1];
104  product[8] = lambda[3] * lambda[4];
105  product[9] = lambda[4] * lambda[1];
106  product[10] = lambda[4] * lambda[0];
107  product[11] = lambda[4] * lambda[2];
108  int indexEdgeBasis = 0;
109  for(int iEdge = 0; iEdge < _nedge; iEdge++) {
110  int uvw = 0;
111  std::vector<double> direction(3, 0);
112  switch(iEdge) {
113  case(0):
114  case(5):
115  case(8):
116  case(11):
117  uvw = 0;
118  direction[0] = 1;
119  break;
120  case(1):
121  case(3):
122  case(9):
123  case(10):
124  uvw = 1;
125  direction[1] = 1;
126  break;
127  case(2):
128  case(4):
129  case(6):;
130  case(7):
131  uvw = 2;
132  direction[2] = 1;
133  break;
134  }
135  for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] + 1;
136  indexEdgeFunc++) {
137  for(int j = 0; j < 3; j++) {
138  edgeBasis[indexEdgeBasis][j] =
139  direction[j] * legendreVector[uvw][indexEdgeFunc] * product[iEdge];
140  }
141  indexEdgeBasis++;
142  }
143  }
144 
145  // face functions:
146  int indexFaceFunction = 0;
147  for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
148  int indexLambda = 0;
149  std::vector<double> direction1(3, 0);
150  std::vector<double> direction2(3, 0);
151  int uvw1 = 0;
152  int uvw2 = 0;
153  switch(iFace) {
154  case(0):
155  indexLambda = 5;
156  uvw1 = 0;
157  uvw2 = 1;
158  direction1[0] = 1;
159  direction2[1] = 1;
160  break;
161  case(1):
162  indexLambda = 3;
163  uvw1 = 0;
164  uvw2 = 2;
165  direction1[0] = 1;
166  direction2[2] = 1;
167  break;
168  case(2):
169  indexLambda = 1;
170  uvw1 = 1;
171  uvw2 = 2;
172  direction1[1] = 1;
173  direction2[2] = 1;
174  break;
175  case(3):
176  indexLambda = 0;
177  uvw1 = 1;
178  uvw2 = 2;
179  direction1[1] = 1;
180  direction2[2] = 1;
181  break;
182  case(4):
183  indexLambda = 2;
184  uvw1 = 0;
185  uvw2 = 2;
186  direction1[0] = 1;
187  direction2[2] = 1;
188  break;
189  case(5):
190  indexLambda = 4;
191  uvw1 = 0;
192  uvw2 = 1;
193  direction1[0] = 1;
194  direction2[1] = 1;
195  break;
196  }
197  for(int index1 = 0; index1 < _pOrderFace1[iFace] + 1; index1++) {
198  for(int index2 = 0; index2 < _pOrderFace2[iFace]; index2++) {
199  for(int j = 0; j < 3; j++) {
200  faceBasis[indexFaceFunction][j] =
201  lambda[indexLambda] * legendreVector[uvw1][index1] *
202  lobattoVector[uvw2][index2] * direction1[j];
203  }
204 
205  indexFaceFunction++;
206  }
207  }
208  for(int index1 = 0; index1 < _pOrderFace1[iFace]; index1++) {
209  for(int index2 = 0; index2 < _pOrderFace2[iFace] + 1; index2++) {
210  for(int j = 0; j < 3; j++) {
211  faceBasis[indexFaceFunction][j] =
212  lambda[indexLambda] * lobattoVector[uvw1][index1] *
213  legendreVector[uvw2][index2] * direction2[j];
214  }
215 
216  indexFaceFunction++;
217  }
218  }
219  }
220  // bubble shape functions:
221  int indexBubbleBasis = 0;
222  for(int ipb1 = 0; ipb1 < _pb1 + 1; ipb1++) {
223  for(int ipb2 = 0; ipb2 < _pb2; ipb2++) {
224  for(int ipb3 = 0; ipb3 < _pb3; ipb3++) {
225  bubbleBasis[indexBubbleBasis][0] = legendreVector[0][ipb1] *
226  lobattoVector[1][ipb2] *
227  lobattoVector[2][ipb3];
228  bubbleBasis[indexBubbleBasis][1] = 0;
229  bubbleBasis[indexBubbleBasis][2] = 0;
230  indexBubbleBasis++;
231  }
232  }
233  }
234  for(int ipb1 = 0; ipb1 < _pb1; ipb1++) {
235  for(int ipb2 = 0; ipb2 < _pb2 + 1; ipb2++) {
236  for(int ipb3 = 0; ipb3 < _pb3; ipb3++) {
237  bubbleBasis[indexBubbleBasis][0] = 0;
238  bubbleBasis[indexBubbleBasis][1] = lobattoVector[0][ipb1] *
239  legendreVector[1][ipb2] *
240  lobattoVector[2][ipb3];
241  bubbleBasis[indexBubbleBasis][2] = 0;
242  indexBubbleBasis++;
243  }
244  }
245  }
246  for(int ipb1 = 0; ipb1 < _pb1; ipb1++) {
247  for(int ipb2 = 0; ipb2 < _pb2; ipb2++) {
248  for(int ipb3 = 0; ipb3 < _pb3 + 1; ipb3++) {
249  bubbleBasis[indexBubbleBasis][0] = 0;
250  bubbleBasis[indexBubbleBasis][1] = 0;
251  bubbleBasis[indexBubbleBasis][2] = lobattoVector[0][ipb1] *
252  lobattoVector[1][ipb2] *
253  legendreVector[2][ipb3];
254  indexBubbleBasis++;
255  }
256  }
257  }
258 }
259 
260 void HierarchicalBasisHcurlBrick::orientEdge(
261  int const &flagOrientation, int const &edgeNumber,
262  std::vector<std::vector<double> > &edgeFunctions,
263  const std::vector<std::vector<double> > &eTablePositiveFlag,
264  const std::vector<std::vector<double> > &eTableNegativeFlag)
265 {
266  if(flagOrientation == -1) {
267  int constant1 = 0;
268  int constant2 = 0;
269  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
270  constant2 = constant2 - 1;
271  constant1 = constant2 - _pOrderEdge[edgeNumber];
272  for(int k = constant1; k <= constant2; k++) {
273  edgeFunctions[k][0] = eTableNegativeFlag[k][0];
274  edgeFunctions[k][1] = eTableNegativeFlag[k][1];
275  edgeFunctions[k][2] = eTableNegativeFlag[k][2];
276  }
277  }
278  else {
279  int constant1 = 0;
280  int constant2 = 0;
281  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
282  constant2 = constant2 - 1;
283  constant1 = constant2 - _pOrderEdge[edgeNumber];
284  for(int k = constant1; k <= constant2; k++) {
285  edgeFunctions[k][0] = eTablePositiveFlag[k][0];
286  edgeFunctions[k][1] = eTablePositiveFlag[k][1];
287  edgeFunctions[k][2] = eTablePositiveFlag[k][2];
288  }
289  }
290 }
291 void HierarchicalBasisHcurlBrick::orientEdgeFunctionsForNegativeFlag(
292  std::vector<std::vector<double> > &edgeFunctions)
293 {
294  int constant1 = 0;
295  int constant2 = 0;
296  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
297  constant2 = 0;
298  constant2 = 0;
299  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
300  constant2 = constant2 - 1;
301  constant1 = constant2 - _pOrderEdge[edgeNumber];
302  for(int k = constant1; k <= constant2; k++) {
303  if((k - constant1) % 2 == 0) {
304  edgeFunctions[k][0] = edgeFunctions[k][0] * (-1);
305  edgeFunctions[k][1] = edgeFunctions[k][1] * (-1);
306  edgeFunctions[k][2] = edgeFunctions[k][2] * (-1);
307  }
308  }
309  }
310 }
311 void HierarchicalBasisHcurlBrick::orientOneFace(
312  double const &u, double const &v, double const &w, int const &flag1,
313  int const &flag2, int const &flag3, int const &faceNumber,
314  std::vector<std::vector<double> > &faceFunctions, std::string typeFunction)
315 {
316  if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
317  int iterator = 0;
318  for(int i = 0; i < faceNumber; i++) {
319  iterator += (_pOrderFace1[i] + 1) * _pOrderFace2[i] +
320  (_pOrderFace2[i] + 1) * _pOrderFace1[i];
321  }
322  if(flag3 == 1) {
323  for(int it1 = 0; it1 < _pOrderFace1[faceNumber] + 1; it1++) {
324  for(int it2 = 0; it2 < _pOrderFace2[faceNumber]; it2++) {
325  int impactFlag1 = 1;
326  int impactFlag2 = 1;
327  if(flag1 == -1 && it1 % 2 == 0) { impactFlag1 = -1; }
328  if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
329  faceFunctions[iterator][0] =
330  faceFunctions[iterator][0] * impactFlag1 * impactFlag2;
331  faceFunctions[iterator][1] =
332  faceFunctions[iterator][1] * impactFlag1 * impactFlag2;
333  faceFunctions[iterator][2] =
334  faceFunctions[iterator][2] * impactFlag1 * impactFlag2;
335  iterator++;
336  }
337  }
338  for(int it1 = 0; it1 < _pOrderFace1[faceNumber]; it1++) {
339  for(int it2 = 0; it2 < _pOrderFace2[faceNumber] + 1; it2++) {
340  int impactFlag1 = 1;
341  int impactFlag2 = 1;
342  if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
343  if(flag2 == -1 && it2 % 2 == 0) { impactFlag2 = -1; }
344  faceFunctions[iterator][0] =
345  faceFunctions[iterator][0] * impactFlag1 * impactFlag2;
346  faceFunctions[iterator][1] =
347  faceFunctions[iterator][1] * impactFlag1 * impactFlag2;
348  faceFunctions[iterator][2] =
349  faceFunctions[iterator][2] * impactFlag1 * impactFlag2;
350  iterator++;
351  }
352  }
353  }
354  else {
355  if(typeFunction == "HcurlLegendre") {
356  std::vector<double> uvw(3);
357  uvw[0] = u;
358  uvw[1] = v;
359  uvw[2] = w;
360  double lambda = 0;
361  int var1 = 0;
362  int var2 = 0;
363  switch(faceNumber) {
364  case(0):
365  lambda = _affineCoordinate(6, u, v, w);
366  var1 = 0;
367  var2 = 1;
368  break;
369  case(1):
370  lambda = _affineCoordinate(4, u, v, w);
371  var1 = 0;
372  var2 = 2;
373  break;
374  case(2):
375  lambda = _affineCoordinate(2, u, v, w);
376  var1 = 1;
377  var2 = 2;
378  break;
379  case(3):
380  lambda = _affineCoordinate(1, u, v, w);
381  var1 = 1;
382  var2 = 2;
383  break;
384  case(4):
385  lambda = _affineCoordinate(3, u, v, w);
386  var1 = 0;
387  var2 = 2;
388  break;
389  case(5):
390  lambda = _affineCoordinate(5, u, v, w);
391  var1 = 0;
392  var2 = 1;
393  break;
394  }
395  std::vector<double> lkVector1(_pOrderFace2[faceNumber]);
396  std::vector<double> lkVector2(_pOrderFace1[faceNumber]);
397  std::vector<double> legendreVector1(_pOrderFace1[faceNumber] + 1);
398  std::vector<double> legendreVector2(_pOrderFace2[faceNumber] + 1);
399  for(int it = 2; it <= _pOrderFace2[faceNumber] + 1; it++) {
400  lkVector1[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[var2]);
401  }
402  for(int it = 2; it <= _pOrderFace1[faceNumber] + 1; it++) {
403  lkVector2[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[var1]);
404  }
405  for(int it = 0; it <= _pOrderFace1[faceNumber]; it++) {
406  legendreVector1[it] = OrthogonalPoly::EvalLegendre(it, uvw[var1]);
407  }
408  for(int it = 0; it <= _pOrderFace2[faceNumber]; it++) {
409  legendreVector2[it] = OrthogonalPoly::EvalLegendre(it, uvw[var2]);
410  }
411  std::vector<double> direction1(3, 0);
412  direction1[var1] = 1;
413  std::vector<double> direction2(3, 0);
414  direction2[var2] = 1;
415  for(int it1 = 0; it1 < _pOrderFace2[faceNumber] + 1; it1++) {
416  for(int it2 = 0; it2 < _pOrderFace1[faceNumber]; it2++) {
417  int impactFlag1 = 1;
418  int impactFlag2 = 1;
419  if(flag2 == -1 && it1 % 2 == 0) { impactFlag1 = -1; }
420  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
421  for(int itVector = 0; itVector < 3; itVector++) {
422  faceFunctions[iterator][itVector] =
423  lambda * legendreVector2[it1] * lkVector2[it2] * impactFlag1 *
424  impactFlag2 * direction2[itVector];
425  }
426  iterator++;
427  }
428  }
429  for(int it1 = 0; it1 < _pOrderFace2[faceNumber]; it1++) {
430  for(int it2 = 0; it2 < _pOrderFace1[faceNumber] + 1; it2++) {
431  int impactFlag1 = 1;
432  int impactFlag2 = 1;
433  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
434  if(flag1 == -1 && it2 % 2 == 0) { impactFlag2 = -1; }
435  for(int itVector = 0; itVector < 3; itVector++) {
436  faceFunctions[iterator][itVector] =
437  lambda * legendreVector1[it2] * lkVector1[it1] * impactFlag1 *
438  impactFlag2 * direction1[itVector];
439  }
440  iterator++;
441  }
442  }
443  }
444  else if(typeFunction == "CurlHcurlLegendre") {
445  std::vector<double> vec1(3, 0);
446  std::vector<double> vec2(3, 0);
447  std::vector<double> uvw(3);
448  uvw[0] = u;
449  uvw[1] = v;
450  uvw[2] = w;
451  int uvw1 = 0;
452  int uvw2 = 0;
453  int il1 = 0;
454  int jl1 = 0;
455  int il2 = 0;
456  int jl2 = 0;
457  switch(faceNumber) {
458  case(0):
459  uvw1 = 0;
460  uvw2 = 1;
461  vec1[1] = -0.5;
462  vec1[2] = -_affineCoordinate(6, u, v, w);
463  il1 = 2;
464  jl1 = 1;
465  vec2[0] = 0.5;
466  vec2[2] = _affineCoordinate(6, u, v, w);
467  il2 = 2;
468  jl2 = 0;
469  break;
470  case(1):
471  uvw1 = 0;
472  uvw2 = 2;
473  vec1[1] = _affineCoordinate(4, u, v, w);
474  vec1[2] = 0.5;
475  il1 = 1;
476  jl1 = 2;
477  vec2[0] = -0.5;
478  vec2[1] = -_affineCoordinate(4, u, v, w);
479  il2 = 1;
480  jl2 = 0;
481  break;
482  case(2):
483  uvw1 = 1;
484  uvw2 = 2;
485  vec1[0] = -_affineCoordinate(2, u, v, w);
486  vec1[2] = -0.5;
487  il1 = 0;
488  jl1 = 2;
489  vec2[0] = _affineCoordinate(2, u, v, w);
490  vec2[1] = 0.5;
491  il2 = 0;
492  jl2 = 1;
493  break;
494  case(3):
495  uvw1 = 1;
496  uvw2 = 2;
497  vec1[0] = -_affineCoordinate(1, u, v, w);
498  vec1[2] = 0.5;
499  il1 = 0;
500  jl1 = 2;
501  vec2[0] = _affineCoordinate(1, u, v, w);
502  vec2[1] = -0.5;
503  il2 = 0;
504  jl2 = 1;
505  break;
506  case(4):
507  uvw1 = 0;
508  uvw2 = 2;
509  vec1[1] = _affineCoordinate(3, u, v, w);
510  vec1[2] = -0.5;
511  il1 = 1;
512  jl1 = 2;
513  vec2[0] = 0.5;
514  vec2[1] = -_affineCoordinate(3, u, v, w);
515  il2 = 1;
516  jl2 = 0;
517  break;
518  case(5):
519  uvw1 = 0;
520  uvw2 = 1;
521  vec1[1] = 0.5;
522  vec1[2] = -_affineCoordinate(5, u, v, w);
523  il1 = 2;
524  jl1 = 1;
525  vec2[0] = -0.5;
526  vec2[2] = _affineCoordinate(5, u, v, w);
527  il2 = 2;
528  jl2 = 0;
529  break;
530  }
531  std::vector<double> lkVector1(_pOrderFace2[faceNumber]);
532  std::vector<double> lkVector2(_pOrderFace1[faceNumber]);
533  std::vector<double> dlkVector1(_pOrderFace2[faceNumber]);
534  std::vector<double> dlkVector2(_pOrderFace1[faceNumber]);
535  std::vector<double> legendreVector1(_pOrderFace1[faceNumber] + 1);
536  std::vector<double> legendreVector2(_pOrderFace2[faceNumber] + 1);
537  for(int it = 2; it <= _pOrderFace2[faceNumber] + 1; it++) {
538  lkVector1[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[uvw2]);
539  dlkVector1[it - 2] = OrthogonalPoly::EvalDLobatto(it, uvw[uvw2]);
540  }
541  for(int it = 2; it <= _pOrderFace1[faceNumber] + 1; it++) {
542  lkVector2[it - 2] = OrthogonalPoly::EvalLobatto(it, uvw[uvw1]);
543  dlkVector2[it - 2] = OrthogonalPoly::EvalDLobatto(it, uvw[uvw1]);
544  }
545  for(int it = 0; it <= _pOrderFace1[faceNumber]; it++) {
546  legendreVector1[it] = OrthogonalPoly::EvalLegendre(it, uvw[uvw1]);
547  }
548  for(int it = 0; it <= _pOrderFace2[faceNumber]; it++) {
549  legendreVector2[it] = OrthogonalPoly::EvalLegendre(it, uvw[uvw2]);
550  }
551  for(int it1 = 0; it1 < _pOrderFace2[faceNumber] + 1; it1++) {
552  for(int it2 = 0; it2 < _pOrderFace1[faceNumber]; it2++) {
553  int impactFlag1 = 1;
554  int impactFlag2 = 1;
555  if(flag2 == -1 && it1 % 2 == 0) { impactFlag1 = -1; }
556  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
557  for(int j = 0; j < 3; j++) {
558  if(j == jl2) {
559  faceFunctions[iterator][j] = vec2[j] * lkVector2[it2] *
560  legendreVector2[it1] *
561  impactFlag1 * impactFlag2;
562  }
563  else if(j == il2) {
564  faceFunctions[iterator][j] = vec2[j] * dlkVector2[it2] *
565  legendreVector2[it1] *
566  impactFlag1 * impactFlag2;
567  }
568  else {
569  faceFunctions[iterator][j] = 0;
570  }
571  }
572  iterator++;
573  }
574  }
575  for(int it1 = 0; it1 < _pOrderFace2[faceNumber]; it1++) {
576  for(int it2 = 0; it2 < _pOrderFace1[faceNumber] + 1; it2++) {
577  int impactFlag1 = 1;
578  int impactFlag2 = 1;
579  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
580  if(flag1 == -1 && it2 % 2 == 0) { impactFlag2 = -1; }
581  for(int j = 0; j < 3; j++) {
582  if(j == jl1) {
583  faceFunctions[iterator][j] = vec1[j] * legendreVector1[it2] *
584  lkVector1[it1] * impactFlag1 *
585  impactFlag2;
586  }
587 
588  else if(j == il1) {
589  faceFunctions[iterator][j] = vec1[j] * legendreVector1[it2] *
590  dlkVector1[it1] * impactFlag1 *
591  impactFlag2;
592  }
593  else {
594  faceFunctions[iterator][j] = 0;
595  }
596  }
597  iterator++;
598  }
599  }
600  }
601  else {
602  throw std::runtime_error("unknown typeFunction");
603  }
604  }
605  }
606 }
607 void HierarchicalBasisHcurlBrick::orientFace(
608  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
609  const std::vector<std::vector<double> > &quadFaceFunctionsAllOrientation,
610  const std::vector<std::vector<double> > &triFaceFunctionsAllOrientation,
611  std::vector<std::vector<double> > &fTableCopy)
612 {
613  int iterator = 0;
614  for(int i = 0; i < faceNumber; i++) {
615  iterator += (_pOrderFace1[i] + 1) * _pOrderFace2[i] +
616  (_pOrderFace2[i] + 1) * _pOrderFace1[i];
617  }
618  int numFaceFunctions =
619  (_pOrderFace1[faceNumber] + 1) * _pOrderFace2[faceNumber] +
620  (_pOrderFace2[faceNumber] + 1) * _pOrderFace1[faceNumber];
621  int iOrientation = numberOrientationQuadFace(flag1, flag2, flag3);
622  int offset = iOrientation * _nQuadFaceFunction;
623  int offset2 = iterator + numFaceFunctions;
624  for(int i = iterator; i < offset2; i++) {
625  fTableCopy[i][0] = quadFaceFunctionsAllOrientation[i + offset][0];
626  fTableCopy[i][1] = quadFaceFunctionsAllOrientation[i + offset][1];
627  fTableCopy[i][2] = quadFaceFunctionsAllOrientation[i + offset][2];
628  }
629 }
630 
631 void HierarchicalBasisHcurlBrick::generateCurlBasis(
632  double const &u, double const &v, double const &w,
633  std::vector<std::vector<double> > &edgeBasis,
634  std::vector<std::vector<double> > &faceBasis,
635  std::vector<std::vector<double> > &bubbleBasis)
636 {
637  std::vector<std::vector<double> > lobattoVector(3);
638  lobattoVector[0] = std::vector<double>(_pb1);
639  lobattoVector[1] = std::vector<double>(_pb2);
640  lobattoVector[2] = std::vector<double>(_pb3);
641  std::vector<std::vector<double> > dlobattoVector(3);
642  dlobattoVector[0] = std::vector<double>(_pb1);
643  dlobattoVector[1] = std::vector<double>(_pb2);
644  dlobattoVector[2] = std::vector<double>(_pb3);
645  for(int it = 2; it <= _pb1 + 1; it++) {
646  lobattoVector[0][it - 2] = OrthogonalPoly::EvalLobatto(it, u);
647  dlobattoVector[0][it - 2] = OrthogonalPoly::EvalDLobatto(it, u);
648  }
649  for(int it = 2; it <= _pb2 + 1; it++) {
650  lobattoVector[1][it - 2] = OrthogonalPoly::EvalLobatto(it, v);
651  dlobattoVector[1][it - 2] = OrthogonalPoly::EvalDLobatto(it, v);
652  }
653  for(int it = 2; it <= _pb3 + 1; it++) {
654  lobattoVector[2][it - 2] = OrthogonalPoly::EvalLobatto(it, w);
655  dlobattoVector[2][it - 2] = OrthogonalPoly::EvalDLobatto(it, w);
656  }
657 
658  std::vector<std::vector<double> > legendreVector(3);
659  legendreVector[0] = std::vector<double>(_pb1 + 1);
660  legendreVector[1] = std::vector<double>(_pb2 + 1);
661  legendreVector[2] = std::vector<double>(_pb3 + 1);
662  for(unsigned int k = 0; k < legendreVector[0].size(); k++) {
663  legendreVector[0][k] = OrthogonalPoly::EvalLegendre(k, u);
664  }
665  for(unsigned int k = 0; k < legendreVector[1].size(); k++) {
666  legendreVector[1][k] = OrthogonalPoly::EvalLegendre(k, v);
667  }
668  for(unsigned int k = 0; k < legendreVector[2].size(); k++) {
669  legendreVector[2][k] = OrthogonalPoly::EvalLegendre(k, w);
670  }
671  std::vector<double> lambda(6, 0);
672  lambda[0] = _affineCoordinate(1, u, v, w);
673  lambda[1] = _affineCoordinate(2, u, v, w);
674  lambda[2] = _affineCoordinate(3, u, v, w);
675  lambda[3] = _affineCoordinate(4, u, v, w);
676  lambda[4] = _affineCoordinate(5, u, v, w);
677  lambda[5] = _affineCoordinate(6, u, v, w);
678 
679  std::vector<double> dlambda(6, 0);
680 
681  dlambda[0] = 0.5;
682  dlambda[1] = -0.5;
683  dlambda[2] = 0.5;
684  dlambda[3] = -0.5;
685  dlambda[4] = 0.5;
686  dlambda[5] = -0.5;
687  std::vector<std::vector<double> > curlProduct(12, std::vector<double>(3, 0));
688  curlProduct[0][1] = lambda[3] * dlambda[5];
689  curlProduct[0][2] = -dlambda[3] * lambda[5];
690  curlProduct[5][1] = lambda[2] * dlambda[5];
691  curlProduct[5][2] = -dlambda[2] * lambda[5];
692  curlProduct[8][1] = lambda[3] * dlambda[4];
693  curlProduct[8][2] = -dlambda[3] * lambda[4];
694  curlProduct[11][1] = lambda[2] * dlambda[4];
695  curlProduct[11][2] = -dlambda[2] * lambda[4];
696 
697  curlProduct[1][0] = -lambda[1] * dlambda[5];
698  curlProduct[1][2] = dlambda[1] * lambda[5];
699  curlProduct[3][0] = -lambda[0] * dlambda[5];
700  curlProduct[3][2] = dlambda[0] * lambda[5];
701  curlProduct[9][0] = -lambda[1] * dlambda[4];
702  curlProduct[9][2] = dlambda[1] * lambda[4];
703  curlProduct[10][0] = -lambda[0] * dlambda[4];
704  curlProduct[10][2] = dlambda[0] * lambda[4];
705 
706  curlProduct[2][0] = lambda[1] * dlambda[3];
707  curlProduct[2][1] = -dlambda[1] * lambda[3];
708  curlProduct[4][0] = lambda[0] * dlambda[3];
709  curlProduct[4][1] = -dlambda[0] * lambda[3];
710  curlProduct[6][0] = lambda[0] * dlambda[2];
711  curlProduct[6][1] = -dlambda[0] * lambda[2];
712  curlProduct[7][0] = lambda[1] * dlambda[2];
713  curlProduct[7][1] = -dlambda[1] * lambda[2];
714 
715  int indexEdgeBasis = 0;
716  for(int iEdge = 0; iEdge < _nedge; iEdge++) {
717  int uvw = 0;
718  switch(iEdge) {
719  case(0):
720  case(5):
721  case(8):
722  case(11): uvw = 0; break;
723  case(1):
724  case(3):
725  case(9):
726  case(10): uvw = 1; break;
727  case(2):
728  case(4):
729  case(6):;
730  case(7): uvw = 2; break;
731  }
732  for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] + 1;
733  indexEdgeFunc++) {
734  for(int j = 0; j < 3; j++) {
735  edgeBasis[indexEdgeBasis][j] =
736  curlProduct[iEdge][j] * legendreVector[uvw][indexEdgeFunc];
737  }
738  indexEdgeBasis++;
739  }
740  }
741  // face functions:
742  int indexFaceFunction = 0;
743  for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
744  std::vector<double> vec1(3, 0);
745  std::vector<double> vec2(3, 0);
746  int uvw1 = 0;
747  int uvw2 = 0;
748  int il1 = 0;
749  int jl1 = 0;
750  int il2 = 0;
751  int jl2 = 0;
752  switch(iFace) {
753  case(0):
754  uvw1 = 0;
755  uvw2 = 1;
756  vec1[1] = dlambda[5];
757  vec1[2] = -lambda[5];
758  il1 = 2;
759  jl1 = 1;
760  vec2[0] = -dlambda[5];
761  vec2[2] = lambda[5];
762  il2 = 2;
763  jl2 = 0;
764  break;
765  case(1):
766  uvw1 = 0;
767  uvw2 = 2;
768  vec1[1] = lambda[3];
769  vec1[2] = -dlambda[3];
770  il1 = 1;
771  jl1 = 2;
772  vec2[0] = dlambda[3];
773  vec2[1] = -lambda[3];
774  il2 = 1;
775  jl2 = 0;
776  break;
777  case(2):
778  uvw1 = 1;
779  uvw2 = 2;
780  vec1[0] = -lambda[1];
781  vec1[2] = dlambda[1];
782  il1 = 0;
783  jl1 = 2;
784  vec2[0] = lambda[1];
785  vec2[1] = -dlambda[1];
786  il2 = 0;
787  jl2 = 1;
788  break;
789  case(3):
790  uvw1 = 1;
791  uvw2 = 2;
792  vec1[0] = -lambda[0];
793  vec1[2] = dlambda[0];
794  il1 = 0;
795  jl1 = 2;
796  vec2[0] = lambda[0];
797  vec2[1] = -dlambda[0];
798  il2 = 0;
799  jl2 = 1;
800  break;
801  case(4):
802  uvw1 = 0;
803  uvw2 = 2;
804  vec1[1] = lambda[2];
805  vec1[2] = -dlambda[2];
806  il1 = 1;
807  jl1 = 2;
808  vec2[0] = dlambda[2];
809  vec2[1] = -lambda[2];
810  il2 = 1;
811  jl2 = 0;
812  break;
813  case(5):
814  uvw1 = 0;
815  uvw2 = 1;
816  vec1[1] = dlambda[4];
817  vec1[2] = -lambda[4];
818  il1 = 2;
819  jl1 = 1;
820  vec2[0] = -dlambda[4];
821  vec2[2] = lambda[4];
822  il2 = 2;
823  jl2 = 0;
824  break;
825  }
826  for(int index1 = 0; index1 < _pOrderFace1[iFace] + 1; index1++) {
827  for(int index2 = 0; index2 < _pOrderFace2[iFace]; index2++) {
828  for(int j = 0; j < 3; j++) {
829  if(j == jl1) {
830  faceBasis[indexFaceFunction][j] = vec1[j] *
831  legendreVector[uvw1][index1] *
832  lobattoVector[uvw2][index2];
833  }
834  else if(j == il1) {
835  faceBasis[indexFaceFunction][j] = vec1[j] *
836  legendreVector[uvw1][index1] *
837  dlobattoVector[uvw2][index2];
838  }
839  else {
840  faceBasis[indexFaceFunction][j] = 0;
841  }
842  }
843  indexFaceFunction++;
844  }
845  }
846  for(int index1 = 0; index1 < _pOrderFace1[iFace]; index1++) {
847  for(int index2 = 0; index2 < _pOrderFace2[iFace] + 1; index2++) {
848  for(int j = 0; j < 3; j++) {
849  if(j == jl2) {
850  faceBasis[indexFaceFunction][j] = vec2[j] *
851  lobattoVector[uvw1][index1] *
852  legendreVector[uvw2][index2];
853  }
854  else if(j == il2) {
855  faceBasis[indexFaceFunction][j] = vec2[j] *
856  dlobattoVector[uvw1][index1] *
857  legendreVector[uvw2][index2];
858  }
859  else {
860  faceBasis[indexFaceFunction][j] = 0;
861  }
862  }
863  indexFaceFunction++;
864  }
865  }
866  }
867  // bubble shape functions:
868  int indexBubbleBasis = 0;
869  for(int ipb1 = 0; ipb1 < _pb1 + 1; ipb1++) {
870  for(int ipb2 = 0; ipb2 < _pb2; ipb2++) {
871  for(int ipb3 = 0; ipb3 < _pb3; ipb3++) {
872  bubbleBasis[indexBubbleBasis][0] = 0;
873  bubbleBasis[indexBubbleBasis][1] = legendreVector[0][ipb1] *
874  lobattoVector[1][ipb2] *
875  dlobattoVector[2][ipb3];
876  bubbleBasis[indexBubbleBasis][2] = -legendreVector[0][ipb1] *
877  dlobattoVector[1][ipb2] *
878  lobattoVector[2][ipb3];
879  indexBubbleBasis++;
880  }
881  }
882  }
883  for(int ipb1 = 0; ipb1 < _pb1; ipb1++) {
884  for(int ipb2 = 0; ipb2 < _pb2 + 1; ipb2++) {
885  for(int ipb3 = 0; ipb3 < _pb3; ipb3++) {
886  bubbleBasis[indexBubbleBasis][0] = -lobattoVector[0][ipb1] *
887  legendreVector[1][ipb2] *
888  dlobattoVector[2][ipb3];
889  bubbleBasis[indexBubbleBasis][1] = 0;
890  bubbleBasis[indexBubbleBasis][2] = dlobattoVector[0][ipb1] *
891  legendreVector[1][ipb2] *
892  lobattoVector[2][ipb3];
893  indexBubbleBasis++;
894  }
895  }
896  }
897  for(int ipb1 = 0; ipb1 < _pb1; ipb1++) {
898  for(int ipb2 = 0; ipb2 < _pb2; ipb2++) {
899  for(int ipb3 = 0; ipb3 < _pb3 + 1; ipb3++) {
900  bubbleBasis[indexBubbleBasis][0] = lobattoVector[0][ipb1] *
901  dlobattoVector[1][ipb2] *
902  legendreVector[2][ipb3];
903  bubbleBasis[indexBubbleBasis][1] = -dlobattoVector[0][ipb1] *
904  lobattoVector[1][ipb2] *
905  legendreVector[2][ipb3];
906  bubbleBasis[indexBubbleBasis][2] = 0;
907 
908  indexBubbleBasis++;
909  }
910  }
911  }
912 }
913 
914 void HierarchicalBasisHcurlBrick::getKeysInfo(
915  std::vector<int> &functionTypeInfo, std::vector<int> &orderInfo)
916 {
917  int it = 0;
918  for(int numEdge = 0; numEdge < 12; numEdge++) {
919  for(int i = 0; i <= _pOrderEdge[numEdge]; i++) {
920  functionTypeInfo[it] = 1;
921  orderInfo[it] = i;
922  it++;
923  }
924  }
925  for(int iFace = 0; iFace < _nfaceQuad; iFace++) {
926  for(int index1 = 0; index1 <= _pOrderFace1[iFace]; index1++) {
927  for(int index2 = 2; index2 <= _pOrderFace2[iFace] + 1; index2++) {
928  functionTypeInfo[it] = 2;
929  orderInfo[it] = std::max(index1, index2);
930  it++;
931  }
932  }
933  for(int index1 = 2; index1 <= _pOrderFace1[iFace] + 1; index1++) {
934  for(int index2 = 0; index2 <= _pOrderFace2[iFace]; index2++) {
935  functionTypeInfo[it] = 2;
936  orderInfo[it] = std::max(index1, index2);
937  it++;
938  }
939  }
940  }
941  for(int ipb1 = 0; ipb1 < _pb1 + 1; ipb1++) {
942  for(int ipb2 = 2; ipb2 <= _pb2 + 1; ipb2++) {
943  for(int ipb3 = 2; ipb3 <= _pb3 + 1; ipb3++) {
944  functionTypeInfo[it] = 3;
945  orderInfo[it] = std::max(std::max(ipb1, ipb2), ipb3);
946  it++;
947  }
948  }
949  }
950  for(int ipb1 = 2; ipb1 <= _pb1 + 1; ipb1++) {
951  for(int ipb2 = 0; ipb2 < _pb2 + 1; ipb2++) {
952  for(int ipb3 = 2; ipb3 <= _pb3 + 1; ipb3++) {
953  functionTypeInfo[it] = 3;
954  orderInfo[it] = std::max(std::max(ipb1, ipb2), ipb3);
955  it++;
956  }
957  }
958  }
959  for(int ipb1 = 2; ipb1 <= _pb1 + 1; ipb1++) {
960  for(int ipb2 = 2; ipb2 <= _pb2 + 1; ipb2++) {
961  for(int ipb3 = 0; ipb3 < _pb3 + 1; ipb3++) {
962  functionTypeInfo[it] = 3;
963  orderInfo[it] = std::max(std::max(ipb1, ipb2), ipb3);
964  it++;
965  }
966  }
967  }
968 }
OrthogonalPoly::EvalLegendre
double EvalLegendre(int order, double x)
Definition: OrthogonalPoly.cpp:415
HierarchicalBasis::_nTriFaceFunction
int _nTriFaceFunction
Definition: HierarchicalBasis.h:27
HierarchicalBasisHcurlBrick::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< std::vector< double > > &faceFunctions, std::string typeFunction)
Definition: HierarchicalBasisHcurlBrick.cpp:311
SurfaceProjectorUtils::vec2
std::array< double, 2 > vec2
Definition: meshOctreeLibOL.cpp:27
HierarchicalBasisHcurlBrick::getKeysInfo
virtual void getKeysInfo(std::vector< int > &functionTypeInfo, std::vector< int > &orderInfo)
Definition: HierarchicalBasisHcurlBrick.cpp:914
HierarchicalBasisHcurlBrick.h
HierarchicalBasis::_nvertex
int _nvertex
Definition: HierarchicalBasis.h:20
HierarchicalBasis::_nBubbleFunction
int _nBubbleFunction
Definition: HierarchicalBasis.h:28
HierarchicalBasisHcurlBrick::generateHcurlBasis
virtual void generateHcurlBasis(double const &u, double const &v, double const &w, std::vector< std::vector< double > > &edgeBasis, std::vector< std::vector< double > > &faceBasis, std::vector< std::vector< double > > &bubbleBasis)
Definition: HierarchicalBasisHcurlBrick.cpp:56
HierarchicalBasisHcurlBrick::_pb2
int _pb2
Definition: HierarchicalBasisHcurlBrick.h:86
HierarchicalBasisHcurlBrick::generateCurlBasis
virtual void generateCurlBasis(double const &u, double const &v, double const &w, std::vector< std::vector< double > > &edgeBasis, std::vector< std::vector< double > > &faceBasis, std::vector< std::vector< double > > &bubbleBasis)
Definition: HierarchicalBasisHcurlBrick.cpp:631
OrthogonalPoly::EvalLobatto
double EvalLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:9
OrthogonalPoly::EvalDLobatto
double EvalDLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:127
HierarchicalBasisHcurlBrick::getNumberOfOrientations
virtual unsigned int getNumberOfOrientations() const
Definition: HierarchicalBasisHcurlBrick.cpp:35
HierarchicalBasisHcurlBrick::_pb3
int _pb3
Definition: HierarchicalBasisHcurlBrick.h:87
HierarchicalBasisHcurlBrick::orientFace
virtual void orientFace(int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, const std::vector< std::vector< double > > &quadFaceFunctionsAllOrientation, const std::vector< std::vector< double > > &triFaceFunctionsAllOrientation, std::vector< std::vector< double > > &fTableCopy)
Definition: HierarchicalBasisHcurlBrick.cpp:607
HierarchicalBasisHcurlBrick::~HierarchicalBasisHcurlBrick
virtual ~HierarchicalBasisHcurlBrick()
Definition: HierarchicalBasisHcurlBrick.cpp:33
HierarchicalBasisHcurlBrick::_pOrderFace2
int _pOrderFace2[6]
Definition: HierarchicalBasisHcurlBrick.h:92
HierarchicalBasisHcurlBrick::orientEdge
virtual void orientEdge(int const &flagOrientation, int const &edgeNumber, std::vector< std::vector< double > > &edgeBasis, const std::vector< std::vector< double > > &eTablePositiveFlag, const std::vector< std::vector< double > > &eTableNegativeFlag)
Definition: HierarchicalBasisHcurlBrick.cpp:260
HierarchicalBasisHcurlBrick::_pOrderEdge
int _pOrderEdge[12]
Definition: HierarchicalBasisHcurlBrick.h:88
HierarchicalBasisHcurlBrick::_pb1
int _pb1
Definition: HierarchicalBasisHcurlBrick.h:85
HierarchicalBasisHcurlBrick::_pOrderFace1
int _pOrderFace1[6]
Definition: HierarchicalBasisHcurlBrick.h:90
HierarchicalBasisHcurlBrick::orientEdgeFunctionsForNegativeFlag
virtual void orientEdgeFunctionsForNegativeFlag(std::vector< std::vector< double > > &edgeFunctions)
Definition: HierarchicalBasisHcurlBrick.cpp:291
HierarchicalBasis::_nQuadFaceFunction
int _nQuadFaceFunction
Definition: HierarchicalBasis.h:26
HierarchicalBasis::_nfaceTri
int _nfaceTri
Definition: HierarchicalBasis.h:23
HierarchicalBasis::_nfaceQuad
int _nfaceQuad
Definition: HierarchicalBasis.h:22
direction
static void direction(Vertex *v1, Vertex *v2, double d[3])
Definition: Geo.cpp:218
HierarchicalBasis::_nVertexFunction
int _nVertexFunction
Definition: HierarchicalBasis.h:24
HierarchicalBasisHcurlBrick::HierarchicalBasisHcurlBrick
HierarchicalBasisHcurlBrick(int order)
Definition: HierarchicalBasisHcurlBrick.cpp:12
HierarchicalBasis::_nEdgeFunction
int _nEdgeFunction
Definition: HierarchicalBasis.h:25
HierarchicalBasisHcurlBrick::_affineCoordinate
static double _affineCoordinate(const int &j, const double &u, const double &v, const double &w)
Definition: HierarchicalBasisHcurlBrick.cpp:40
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