gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
HierarchicalBasisHcurlTetra.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 <stdexcept>
12 #include "HierarchicalBasisHcurlTetra.h"
13 
14 HierarchicalBasisHcurlTetra::HierarchicalBasisHcurlTetra(int order)
15 {
16  _nvertex = 4;
17  _nedge = 6;
18  _nfaceTri = 4;
19  _nfaceQuad = 0;
20  _nVertexFunction = 0;
21  _nEdgeFunction = 6 * order + 6;
22  _nQuadFaceFunction = 0;
23  if(order == 0) {
24  _nTriFaceFunction = 0;
25  _nBubbleFunction = 0;
26  }
27  else {
28  _nTriFaceFunction = 12 * (order - 1) + 4 * (order - 2) * (order - 1);
29  _nBubbleFunction = (order - 1) * (order - 2) * (order - 3) / 2 +
30  2 * (order - 2) * (order - 1);
31  }
32  _pb = order;
33  for(int i = 0; i < 4; i++) { _pOrderFace[i] = order; }
34  for(int i = 0; i < 6; i++) { _pOrderEdge[i] = order; }
35 }
36 
37 HierarchicalBasisHcurlTetra::~HierarchicalBasisHcurlTetra() {}
38 
39 unsigned int HierarchicalBasisHcurlTetra::getNumberOfOrientations() const
40 {
41  return 24; // factorial 4
42 }
43 
44 double HierarchicalBasisHcurlTetra::_affineCoordinate(const int &j,
45  const double &u,
46  const double &v,
47  const double &w)
48 {
49  switch(j) {
50  case(1): return 0.5 * (1 + v);
51  case(2): return -0.5 * (1 + u + v + w);
52  case(3): return 0.5 * (1 + u);
53  case(4): return 0.5 * (1 + w);
54  default: throw std::runtime_error("j must be : 1<=j<=4");
55  }
56 }
57 
58 double HierarchicalBasisHcurlTetra::dotProduct(const std::vector<double> &u,
59  const std::vector<double> &v)
60 {
61  return u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
62 }
63 
64 void HierarchicalBasisHcurlTetra::generateHcurlBasis(
65  double const &u, double const &v, double const &w,
66  std::vector<std::vector<double> > &edgeBasis,
67  std::vector<std::vector<double> > &faceBasis,
68  std::vector<std::vector<double> > &bubbleBasis)
69 {
70  //***
71  // to map onto the reference domain of gmsh:
72  double uc = 2 * u - 1;
73  double vc = 2 * v - 1;
74  double wc = 2 * w - 1;
75  double jacob = 2;
76  //*****
77  // compute all needed terms
78  double lambda1 = _affineCoordinate(1, uc, vc, wc);
79  double lambda2 = _affineCoordinate(2, uc, vc, wc);
80  double lambda3 = _affineCoordinate(3, uc, vc, wc);
81  double lambda4 = _affineCoordinate(4, uc, vc, wc);
82  std::vector<double> n1 = std::vector<double>(3, 0);
83  n1[1] = 1;
84  std::vector<double> n2 = std::vector<double>(3, 0);
85  n2[0] = -1. / sqrt(3);
86  n2[1] = -1. / sqrt(3);
87  n2[2] = -1. / sqrt(3);
88  std::vector<double> n3 = std::vector<double>(3, 0);
89  n3[0] = 1;
90  std::vector<double> n4 = std::vector<double>(3, 0);
91  n4[2] = 1;
92  std::vector<double> t1 = std::vector<double>(3, 0);
93  t1[0] = 1;
94  std::vector<double> t2 = std::vector<double>(3, 0);
95  t2[0] = -1;
96  t2[1] = 1;
97  std::vector<double> t3 = std::vector<double>(3, 0);
98  t3[1] = -1;
99  std::vector<double> t4 = std::vector<double>(3, 0);
100  t4[2] = 1;
101  std::vector<double> t5 = std::vector<double>(3, 0);
102  t5[0] = -1;
103  t5[2] = 1;
104  std::vector<double> t6 = std::vector<double>(3, 0);
105  t6[1] = -1;
106  t6[2] = 1;
107  // Whitney functions:
108  std::vector<std::vector<double> > psie_0(6, std::vector<double>(3, 0));
109  std::vector<std::vector<double> > psie_1(6, std::vector<double>(3, 0));
110  for(int i = 0; i < 3; i++) {
111  psie_0[0][i] = lambda3 * n2[i] / dotProduct(n2, t1) +
112  lambda2 * n3[i] / dotProduct(n3, t1);
113  psie_0[1][i] = lambda1 * n3[i] / dotProduct(n3, t2) +
114  lambda3 * n1[i] / dotProduct(n1, t2);
115  psie_0[2][i] = lambda2 * n1[i] / dotProduct(n1, t3) +
116  lambda1 * n2[i] / dotProduct(n2, t3);
117  psie_0[3][i] = lambda4 * n2[i] / dotProduct(n2, t4) +
118  lambda2 * n4[i] / dotProduct(n4, t4);
119  psie_0[4][i] = lambda4 * n1[i] / dotProduct(n1, t6) +
120  lambda1 * n4[i] / dotProduct(n4, t6);
121  psie_0[5][i] = lambda4 * n3[i] / dotProduct(n3, t5) +
122  lambda3 * n4[i] / dotProduct(n4, t5);
123 
124  psie_1[0][i] = lambda3 * n2[i] / dotProduct(n2, t1) -
125  lambda2 * n3[i] / dotProduct(n3, t1);
126  psie_1[1][i] = lambda1 * n3[i] / dotProduct(n3, t2) -
127  lambda3 * n1[i] / dotProduct(n1, t2);
128  psie_1[2][i] = lambda2 * n1[i] / dotProduct(n1, t3) -
129  lambda1 * n2[i] / dotProduct(n2, t3);
130  psie_1[3][i] = lambda4 * n2[i] / dotProduct(n2, t4) -
131  lambda2 * n4[i] / dotProduct(n4, t4);
132  psie_1[4][i] = lambda4 * n1[i] / dotProduct(n1, t6) -
133  lambda1 * n4[i] / dotProduct(n4, t6);
134  psie_1[5][i] = lambda4 * n3[i] / dotProduct(n3, t5) -
135  lambda3 * n4[i] / dotProduct(n4, t5);
136  }
137 
138  std::vector<double> sub(9, 0);
139  sub[0] = lambda3 - lambda2;
140  sub[1] = lambda1 - lambda3;
141  sub[2] = lambda2 - lambda1;
142  sub[3] = lambda4 - lambda2;
143  sub[4] = lambda4 - lambda1;
144  sub[5] = lambda4 - lambda3;
145  sub[6] = lambda2 - lambda4;
146  sub[7] = -lambda2 + lambda1;
147  sub[8] = lambda3 - lambda4;
148  std::vector<std::vector<double> > legendreVector(9);
149  legendreVector[0] = std::vector<double>(std::max(
150  std::max(std::max(_pOrderEdge[0], _pOrderFace[0] - 1), _pOrderFace[1] - 1),
151  0));
152  legendreVector[1] = std::vector<double>(std::max(
153  std::max(std::max(_pOrderEdge[1], _pOrderFace[0] - 1), _pOrderFace[3] - 1),
154  0));
155  legendreVector[2] = std::vector<double>(std::max(
156  std::max(std::max(_pOrderEdge[2], _pOrderFace[0] - 1), _pOrderFace[2] - 1),
157  0));
158  legendreVector[3] = std::vector<double>(_pOrderEdge[3], 0);
159  legendreVector[4] = std::vector<double>(std::max(
160  std::max(std::max(_pOrderEdge[4], _pOrderFace[2] - 1), _pOrderFace[3] - 1),
161  0));
162  legendreVector[5] = std::vector<double>(
163  std::max(std::max(_pOrderEdge[5], _pOrderFace[1] - 1), 0));
164  legendreVector[6] = std::vector<double>(
165  std::max(std::max(_pOrderFace[1] - 1, _pOrderFace[2] - 1), 0));
166  legendreVector[7] = std::vector<double>(std::max(_pOrderFace[2] - 1, 0));
167  legendreVector[8] = std::vector<double>(std::max(_pOrderFace[3] - 1, 0));
168  for(int i = 0; i < 9; i++) {
169  for(unsigned int j = 0; j < legendreVector[i].size(); j++) {
170  legendreVector[i][j] = OrthogonalPoly::EvalLegendre(j, sub[i]);
171  }
172  }
173  // compute edge Basis
174  int edgeIt = 0;
175  for(int i = 0; i < _nedge; i++) {
176  for(int j = 0; j < 3; j++) { edgeBasis[edgeIt][j] = jacob * psie_0[i][j]; }
177  edgeIt++;
178  if(_pOrderEdge[i] >= 1) {
179  for(int j = 0; j < 3; j++) {
180  edgeBasis[edgeIt][j] = jacob * psie_1[i][j];
181  }
182  edgeIt++;
183  for(int iedge = 2; iedge <= _pOrderEdge[i]; iedge++) {
184  for(int j = 0; j < 3; j++) {
185  edgeBasis[edgeIt][j] =
186  jacob * ((2 * float(iedge) - 1) / float(iedge) *
187  legendreVector[i][iedge - 1] * psie_1[i][j] -
188  (float(iedge) - 1) / float(iedge) *
189  legendreVector[i][iedge - 2] * psie_0[i][j]);
190  }
191  edgeIt++;
192  }
193  }
194  }
195 
196  // edge-based face functions , Genuine face functions
197  // and Face-based interior functions
198  int faceIt = 0;
199  int bubbleIt = 0;
200  for(int nFace = 0; nFace < _nfaceTri; nFace++) {
201  int indexVector1 = 0;
202  int indexVector2 = 0;
203  std::vector<double> vecTangent1(3, 0);
204  std::vector<double> vecTangent2(3, 0);
205  std::vector<double> niT(3, 0);
206  double faceProduct = 0;
207  switch(nFace) {
208  case(0):
209  for(int i = 0; i < 3; i++) {
210  double product = 0;
211  std::vector<double> nD(3, 0);
212  int index = 0;
213  switch(i) {
214  case(0):
215  product = lambda3 * lambda2;
216  nD[1] = 0.5;
217  index = 0;
218  break;
219  case(1):
220  product = lambda1 * lambda3;
221  nD[0] = -0.5;
222  nD[1] = -0.5;
223  nD[2] = -0.5;
224  index = 1;
225  break;
226  case(2):
227  product = lambda1 * lambda2;
228  nD[0] = 0.5;
229  index = 2;
230  break;
231  }
232  for(int i1 = 2; i1 <= _pOrderFace[0]; i1++) {
233  for(int j = 0; j < 3; j++) {
234  faceBasis[faceIt][j] =
235  jacob * product * legendreVector[index][i1 - 2] * nD[j];
236  }
237  faceIt++;
238  }
239  }
240 
241  indexVector1 = 0;
242  vecTangent1[0] = 0.5;
243  indexVector2 = 2;
244  vecTangent2[1] = 0.5;
245  niT[2] = 0.5;
246  faceProduct = lambda1 * lambda2 * lambda3;
247  break;
248  case(1):
249  for(int i = 0; i < 3; i++) {
250  double product = 0;
251  std::vector<double> nD(3, 0);
252  int index = 0;
253  switch(i) {
254  case(0):
255  product = lambda3 * lambda2;
256  nD[2] = 0.5;
257  index = 0;
258  break;
259  case(1):
260  product = lambda4 * lambda3;
261  nD[0] = -0.5;
262  nD[1] = -0.5;
263  nD[2] = -0.5;
264  index = 5;
265  break;
266  case(2):
267  product = lambda4 * lambda2;
268  nD[0] = 0.5;
269  index = 6;
270  break;
271  }
272  for(int i1 = 2; i1 <= _pOrderFace[1]; i1++) {
273  for(int j = 0; j < 3; j++) {
274  faceBasis[faceIt][j] =
275  jacob * product * legendreVector[index][i1 - 2] * nD[j];
276  }
277  faceIt++;
278  }
279  }
280  indexVector1 = 0;
281  vecTangent1[0] = 0.5;
282  indexVector2 = 6;
283  vecTangent2[2] = 0.5;
284  niT[1] = 0.5;
285  faceProduct = lambda4 * lambda2 * lambda3;
286  break;
287  case(2):
288  for(int i = 0; i < 3; i++) {
289  double product = 0;
290  std::vector<double> nD(3, 0);
291  int index = 0;
292  switch(i) {
293  case(0):
294  product = lambda1 * lambda2;
295  nD[2] = 0.5;
296  index = 7;
297  break;
298  case(1):
299  product = lambda4 * lambda1;
300  nD[0] = -0.5;
301  nD[1] = -0.5;
302  nD[2] = -0.5;
303  index = 4;
304  break;
305  case(2):
306  product = lambda4 * lambda2;
307  nD[1] = 0.5;
308  index = 6;
309  break;
310  }
311  for(int i1 = 2; i1 <= _pOrderFace[2]; i1++) {
312  for(int j = 0; j < 3; j++) {
313  faceBasis[faceIt][j] =
314  jacob * product * legendreVector[index][i1 - 2] * nD[j];
315  }
316  faceIt++;
317  }
318  }
319  indexVector1 = 7;
320  vecTangent1[1] = 0.5;
321  indexVector2 = 6;
322  vecTangent2[2] = 0.5;
323  niT[0] = 0.5;
324  faceProduct = lambda1 * lambda4 * lambda2;
325  break;
326  case(3):
327  for(int i = 0; i < 3; i++) {
328  double product = 0;
329  std::vector<double> nD(3, 0);
330  int index = 0;
331  switch(i) {
332  case(0):
333  product = lambda1 * lambda3;
334  nD[2] = 0.5;
335  index = 1;
336  break;
337  case(1):
338  product = lambda4 * lambda1;
339  nD[0] = 0.5;
340  index = 4;
341  break;
342  case(2):
343  product = lambda4 * lambda3;
344  nD[1] = 0.5;
345  index = 8;
346  break;
347  }
348  for(int i1 = 2; i1 <= _pOrderFace[3]; i1++) {
349  for(int j = 0; j < 3; j++) {
350  faceBasis[faceIt][j] =
351  jacob * product * legendreVector[index][i1 - 2] * nD[j];
352  }
353  faceIt++;
354  }
355  }
356  indexVector1 = 1;
357  vecTangent1[1] = 0.5;
358  indexVector2 = 8;
359  vecTangent2[2] = 0.5;
360  niT[0] = -0.5;
361  niT[1] = -0.5;
362  niT[2] = -0.5;
363  faceProduct = lambda1 * lambda3 * lambda4;
364  break;
365  }
366  std::vector<double> copy(
367  int((_pOrderFace[nFace] - 2) * (_pOrderFace[nFace] - 1) / 2), 0);
368  int itCopy = 0;
369  for(int n1 = 0; n1 < _pOrderFace[nFace] - 2; n1++) {
370  for(int n2 = 0; n2 < _pOrderFace[nFace] - 2 - n1; n2++) {
371  copy[itCopy] = jacob * faceProduct * legendreVector[indexVector1][n1] *
372  legendreVector[indexVector2][n2];
373  for(int i = 0; i < 3; i++) {
374  faceBasis[faceIt][i] = copy[itCopy] * vecTangent1[i];
375  }
376  itCopy++;
377  faceIt++;
378  }
379  }
380  itCopy = 0;
381  for(int n1 = 0; n1 < _pOrderFace[nFace] - 2; n1++) {
382  for(int n2 = 0; n2 < _pOrderFace[nFace] - 2 - n1; n2++) {
383  for(int i = 0; i < 3; i++) {
384  faceBasis[faceIt][i] = copy[itCopy] * vecTangent2[i];
385  }
386  faceIt++;
387  itCopy++;
388  }
389  }
390 
391  for(int n1 = 0; n1 < _pb - 2; n1++) {
392  for(int n2 = 0; n2 < _pb - 2 - n1; n2++) {
393  for(int i = 0; i < 3; i++) {
394  bubbleBasis[bubbleIt][i] = jacob * faceProduct *
395  legendreVector[indexVector1][n1] *
396  legendreVector[indexVector2][n2] * niT[i];
397  }
398  bubbleIt++;
399  }
400  }
401  }
402  if(_pb > 3) {
403  int itCopy = 0;
404  std::vector<double> copyBubbleTerm((_pb - 3) * (_pb - 2) * (_pb - 1) / 6.,
405  0);
406  for(int n1 = 0; n1 < _pb - 3; n1++) {
407  double phi1 = OrthogonalPoly::EvalKernelFunction(n1, sub[7]);
408  for(int n2 = 0; n2 < _pb - 3 - n1; n2++) {
409  double phi2 = OrthogonalPoly::EvalKernelFunction(n2, sub[0]);
410  for(int n3 = 0; n3 < _pb - 3 - n2 - n1; n3++) {
411  copyBubbleTerm[itCopy] =
412  jacob * lambda1 * lambda2 * lambda3 * lambda4 * phi1 * phi2 *
413  OrthogonalPoly::EvalKernelFunction(n3, sub[3]);
414  bubbleBasis[bubbleIt][0] = copyBubbleTerm[itCopy];
415  itCopy++;
416  bubbleIt++;
417  }
418  }
419  }
420  itCopy = 0;
421  for(int n1 = 0; n1 < _pb - 3; n1++) {
422  for(int n2 = 0; n2 < _pb - 3 - n1; n2++) {
423  for(int n3 = 0; n3 < _pb - 3 - n2 - n1; n3++) {
424  bubbleBasis[bubbleIt][1] = copyBubbleTerm[itCopy];
425  itCopy++;
426  bubbleIt++;
427  }
428  }
429  }
430  itCopy = 0;
431  for(int n1 = 0; n1 < _pb - 3; n1++) {
432  for(int n2 = 0; n2 < _pb - 3 - n1; n2++) {
433  for(int n3 = 0; n3 < _pb - 3 - n2 - n1; n3++) {
434  bubbleBasis[bubbleIt][2] = copyBubbleTerm[itCopy];
435  itCopy++;
436  bubbleIt++;
437  }
438  }
439  }
440  }
441 }
442 
443 void HierarchicalBasisHcurlTetra::orientEdge(
444  int const &flagOrientation, int const &edgeNumber,
445  std::vector<std::vector<double> > &edgeFunctions,
446  const std::vector<std::vector<double> > &eTablePositiveFlag,
447  const std::vector<std::vector<double> > &eTableNegativeFlag)
448 {
449  if(flagOrientation == -1) {
450  int constant1 = 0;
451  int constant2 = 0;
452  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
453  constant2 = constant2 - 1;
454  constant1 = constant2 - _pOrderEdge[edgeNumber];
455  for(int k = constant1; k <= constant2; k++) {
456  edgeFunctions[k][0] = eTableNegativeFlag[k][0];
457  edgeFunctions[k][1] = eTableNegativeFlag[k][1];
458  edgeFunctions[k][2] = eTableNegativeFlag[k][2];
459  }
460  }
461  else {
462  int constant1 = 0;
463  int constant2 = 0;
464  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
465  constant2 = constant2 - 1;
466  constant1 = constant2 - _pOrderEdge[edgeNumber];
467  for(int k = constant1; k <= constant2; k++) {
468  edgeFunctions[k][0] = eTablePositiveFlag[k][0];
469  edgeFunctions[k][1] = eTablePositiveFlag[k][1];
470  edgeFunctions[k][2] = eTablePositiveFlag[k][2];
471  }
472  }
473 }
474 void HierarchicalBasisHcurlTetra::orientEdgeFunctionsForNegativeFlag(
475  std::vector<std::vector<double> > &edgeFunctions)
476 {
477  int constant1 = 0;
478  int constant2 = 0;
479  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
480  constant2 = 0;
481  constant2 = 0;
482  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
483  constant2 = constant2 - 1;
484  constant1 = constant2 - _pOrderEdge[edgeNumber];
485  for(int k = constant1; k <= constant2; k++) {
486  if((k - constant1) % 2 == 0) {
487  edgeFunctions[k][0] = edgeFunctions[k][0] * (-1);
488  edgeFunctions[k][1] = edgeFunctions[k][1] * (-1);
489  edgeFunctions[k][2] = edgeFunctions[k][2] * (-1);
490  }
491  }
492  }
493 }
494 void HierarchicalBasisHcurlTetra::orientOneFace(
495  double const &u, double const &v, double const &w, int const &flag1,
496  int const &flag2, int const &flag3, int const &faceNumber,
497  std::vector<std::vector<double> > &faceFunctions, std::string typeFunction)
498 {
499  if(!(flag1 == 0 && flag2 == 1)) {
500  if(typeFunction == "HcurlLegendre") {
501  // orient Edge-based interior functions
502  // to map onto the reference domain of gmsh:
503  double uc = 2 * u - 1;
504  double vc = 2 * v - 1;
505  double wc = 2 * w - 1;
506  double jacob = 2;
507  //*****
508 
509  int faceIt = 0;
510  for(int k = 0; k < faceNumber; k++) {
511  faceIt = faceIt + 3 * (_pOrderFace[k] - 1) +
512  (_pOrderFace[k] - 2) * (_pOrderFace[k] - 1);
513  }
514  std::vector<double> lambda(3);
515  std::vector<std::vector<double> > dlambda(3, std::vector<double>(3, 0));
516  switch(faceNumber) {
517  case(0):
518  lambda[0] = _affineCoordinate(2, uc, vc, wc);
519  lambda[1] = _affineCoordinate(3, uc, vc, wc);
520  lambda[2] = _affineCoordinate(1, uc, vc, wc);
521  dlambda[0][0] = -0.5;
522  dlambda[0][1] = -0.5;
523  dlambda[0][2] = -0.5;
524  dlambda[1][0] = 0.5;
525  dlambda[2][1] = 0.5;
526  break;
527  case(1):
528  lambda[0] = _affineCoordinate(2, uc, vc, wc);
529  lambda[1] = _affineCoordinate(3, uc, vc, wc);
530  lambda[2] = _affineCoordinate(4, uc, vc, wc);
531  dlambda[0][0] = -0.5;
532  dlambda[0][1] = -0.5;
533  dlambda[0][2] = -0.5;
534  dlambda[1][0] = 0.5;
535  dlambda[2][2] = 0.5;
536  break;
537  case(2):
538  lambda[0] = _affineCoordinate(2, uc, vc, wc);
539  lambda[1] = _affineCoordinate(1, uc, vc, wc);
540  lambda[2] = _affineCoordinate(4, uc, vc, wc);
541  dlambda[0][0] = -0.5;
542  dlambda[0][1] = -0.5;
543  dlambda[0][2] = -0.5;
544  dlambda[1][1] = 0.5;
545  dlambda[2][2] = 0.5;
546  break;
547  case(3):
548  lambda[0] = _affineCoordinate(3, uc, vc, wc);
549  lambda[1] = _affineCoordinate(1, uc, vc, wc);
550  lambda[2] = _affineCoordinate(4, uc, vc, wc);
551  dlambda[0][0] = 0.5;
552  dlambda[1][1] = 0.5;
553  dlambda[2][2] = 0.5;
554  break;
555  }
556  if(flag1 == 1 && flag2 == -1) {
557  double copy = lambda[0];
558  lambda[0] = lambda[1];
559  lambda[1] = copy;
560  std::vector<double> dcopy = dlambda[0];
561  dlambda[0] = dlambda[1];
562  dlambda[1] = dcopy;
563  }
564  else if(flag1 == 0 && flag2 == -1) {
565  double copy = lambda[2];
566  lambda[2] = lambda[1];
567  lambda[1] = copy;
568  std::vector<double> dcopy = dlambda[2];
569  dlambda[2] = dlambda[1];
570  dlambda[1] = dcopy;
571  }
572  else if(flag1 == 2 && flag2 == -1) {
573  double copy = lambda[2];
574  lambda[2] = lambda[0];
575  lambda[0] = copy;
576  std::vector<double> dcopy = dlambda[2];
577  dlambda[2] = dlambda[0];
578  dlambda[0] = dcopy;
579  }
580  else if(flag1 == 1 && flag2 == 1) {
581  double copy = lambda[0];
582  lambda[0] = lambda[1];
583  lambda[1] = lambda[2];
584  lambda[2] = copy;
585  std::vector<double> dcopy = dlambda[0];
586  dlambda[0] = dlambda[1];
587  dlambda[1] = dlambda[2];
588  dlambda[2] = dcopy;
589  }
590  else if(flag1 == 2 && flag2 == 1) {
591  double copy = lambda[0];
592  lambda[0] = lambda[2];
593  lambda[2] = lambda[1];
594  lambda[1] = copy;
595  std::vector<double> dcopy = dlambda[0];
596  dlambda[0] = dlambda[2];
597  dlambda[2] = dlambda[1];
598  dlambda[1] = dcopy;
599  }
600  std::vector<double> sub(3);
601  sub[0] = lambda[1] - lambda[0];
602  sub[1] = lambda[2] - lambda[1];
603  sub[2] = lambda[0] - lambda[2];
604  std::vector<double> n1 = std::vector<double>(3, 0);
605  n1[0] = dlambda[2][0];
606  n1[1] = dlambda[2][1];
607  n1[2] = dlambda[2][2];
608  std::vector<double> n2 = std::vector<double>(3, 0);
609  n2[0] = dlambda[0][0];
610  n2[1] = dlambda[0][1];
611  n2[2] = dlambda[0][2];
612  std::vector<double> n3 = std::vector<double>(3, 0);
613  n3[0] = dlambda[1][0];
614  n3[1] = dlambda[1][1];
615  n3[2] = dlambda[1][2];
616  // edge-based face functions
617  for(int i = 0; i < 3; i++) {
618  double product2 = 0;
619  std::vector<double> *normal(nullptr);
620  switch(i) {
621  case(0):
622  product2 = lambda[1] * lambda[0];
623  normal = &n1;
624  break;
625  case(1):
626  product2 = lambda[1] * lambda[2];
627  normal = &n2;
628  break;
629  case(2):
630  product2 = lambda[2] * lambda[0];
631  normal = &n3;
632  break;
633  }
634  for(int i1 = 2; i1 <= _pOrderFace[faceNumber]; i1++) {
635  faceFunctions[faceIt][0] =
636  jacob * product2 * (*normal)[0] *
637  OrthogonalPoly::EvalLegendre(i1 - 2, sub[i]);
638  faceFunctions[faceIt][1] =
639  jacob * product2 * (*normal)[1] *
640  OrthogonalPoly::EvalLegendre(i1 - 2, sub[i]);
641  faceFunctions[faceIt][2] =
642  jacob * product2 * (*normal)[2] *
643  OrthogonalPoly::EvalLegendre(i1 - 2, sub[i]);
644  faceIt++;
645  }
646  }
647  // Genuine face functions
648  double sub1 = sub[0];
649  double sub2 = sub[2];
650  std::vector<double> LSub2(_pOrderFace[faceNumber] - 2);
651  for(int it = 0; it < _pOrderFace[faceNumber] - 2; it++) {
652  LSub2[it] = OrthogonalPoly::EvalLegendre(it, sub2);
653  }
654  double product = lambda[0] * lambda[1] * lambda[2];
655 
656  std::vector<double> copy(int((_pOrderFace[faceNumber] - 1) *
657  (_pOrderFace[faceNumber] - 2) / 2.));
658  int it1 = 0;
659  for(int n1 = 0; n1 < _pOrderFace[faceNumber] - 2; n1++) {
660  double LSub1 = OrthogonalPoly::EvalLegendre(n1, sub1);
661  for(int n2 = 0; n2 < _pOrderFace[faceNumber] - 2 - n1; n2++) {
662  copy[it1] = jacob * product * LSub1 * LSub2[n2];
663  for(int i = 0; i < 3; i++) {
664  faceFunctions[faceIt][i] = copy[it1] * dlambda[1][i];
665  }
666  it1++;
667  faceIt++;
668  }
669  }
670  it1 = 0;
671  for(int n1 = 0; n1 < _pOrderFace[faceNumber] - 2; n1++) {
672  for(int n2 = 0; n2 < _pOrderFace[faceNumber] - 2 - n1; n2++) {
673  for(int i = 0; i < 3; i++) {
674  faceFunctions[faceIt][i] = copy[it1] * dlambda[2][i];
675  }
676  faceIt++;
677  it1++;
678  }
679  }
680  }
681 
682  else if("CurlHcurlLegendre" == typeFunction) {
683  // to map onto the reference domain of gmsh:
684  double uc = 2 * u - 1;
685  double vc = 2 * v - 1;
686  double wc = 2 * w - 1;
687  double jacob = 4; // 8*0.5
688  //*****
689  int faceIt = 0;
690  for(int k = 0; k < faceNumber; k++) {
691  faceIt = faceIt + 3 * (_pOrderFace[k] - 1) +
692  (_pOrderFace[k] - 2) * (_pOrderFace[k] - 1);
693  }
694  std::vector<double> lambda(3);
695  std::vector<std::vector<double> > dlambda(3, std::vector<double>(3, 0));
696  switch(faceNumber) {
697  case(0):
698  lambda[0] = _affineCoordinate(2, uc, vc, wc);
699  lambda[1] = _affineCoordinate(3, uc, vc, wc);
700  lambda[2] = _affineCoordinate(1, uc, vc, wc);
701  dlambda[0][0] = -0.5;
702  dlambda[0][1] = -0.5;
703  dlambda[0][2] = -0.5;
704  dlambda[1][0] = 0.5;
705  dlambda[2][1] = 0.5;
706  break;
707  case(1):
708  lambda[0] = _affineCoordinate(2, uc, vc, wc);
709  lambda[1] = _affineCoordinate(3, uc, vc, wc);
710  lambda[2] = _affineCoordinate(4, uc, vc, wc);
711  dlambda[0][0] = -0.5;
712  dlambda[0][1] = -0.5;
713  dlambda[0][2] = -0.5;
714  dlambda[1][0] = 0.5;
715  dlambda[2][2] = 0.5;
716  break;
717  case(2):
718  lambda[0] = _affineCoordinate(2, uc, vc, wc);
719  lambda[1] = _affineCoordinate(1, uc, vc, wc);
720  lambda[2] = _affineCoordinate(4, uc, vc, wc);
721  dlambda[0][0] = -0.5;
722  dlambda[0][1] = -0.5;
723  dlambda[0][2] = -0.5;
724  dlambda[1][1] = 0.5;
725  dlambda[2][2] = 0.5;
726  break;
727  case(3):
728  lambda[0] = _affineCoordinate(3, uc, vc, wc);
729  lambda[1] = _affineCoordinate(1, uc, vc, wc);
730  lambda[2] = _affineCoordinate(4, uc, vc, wc);
731  dlambda[0][0] = 0.5;
732  dlambda[1][1] = 0.5;
733  dlambda[2][2] = 0.5;
734  break;
735  }
736  if(flag1 == 1 && flag2 == -1) {
737  double copy = lambda[0];
738  lambda[0] = lambda[1];
739  lambda[1] = copy;
740  std::vector<double> dcopy = dlambda[0];
741  dlambda[0] = dlambda[1];
742  dlambda[1] = dcopy;
743  }
744  else if(flag1 == 0 && flag2 == -1) {
745  double copy = lambda[2];
746  lambda[2] = lambda[1];
747  lambda[1] = copy;
748  std::vector<double> dcopy = dlambda[2];
749  dlambda[2] = dlambda[1];
750  dlambda[1] = dcopy;
751  }
752  else if(flag1 == 2 && flag2 == -1) {
753  double copy = lambda[2];
754  lambda[2] = lambda[0];
755  lambda[0] = copy;
756  std::vector<double> dcopy = dlambda[2];
757  dlambda[2] = dlambda[0];
758  dlambda[0] = dcopy;
759  }
760  else if(flag1 == 1 && flag2 == 1) {
761  double copy = lambda[0];
762  lambda[0] = lambda[1];
763  lambda[1] = lambda[2];
764  lambda[2] = copy;
765  std::vector<double> dcopy = dlambda[0];
766  dlambda[0] = dlambda[1];
767  dlambda[1] = dlambda[2];
768  dlambda[2] = dcopy;
769  }
770  else if(flag1 == 2 && flag2 == 1) {
771  double copy = lambda[0];
772  lambda[0] = lambda[2];
773  lambda[2] = lambda[1];
774  lambda[1] = copy;
775  std::vector<double> dcopy = dlambda[0];
776  dlambda[0] = dlambda[2];
777  dlambda[2] = dlambda[1];
778  dlambda[1] = dcopy;
779  }
780  std::vector<double> sub(3);
781  sub[0] = lambda[1] - lambda[0];
782  sub[1] = lambda[2] - lambda[1];
783  sub[2] = lambda[0] - lambda[2];
784  std::vector<double> n1 = std::vector<double>(3, 0);
785  n1[0] = dlambda[2][0];
786  n1[1] = dlambda[2][1];
787  n1[2] = dlambda[2][2];
788  std::vector<double> n2 = std::vector<double>(3, 0);
789  n2[0] = dlambda[0][0];
790  n2[1] = dlambda[0][1];
791  n2[2] = dlambda[0][2];
792  std::vector<double> n3 = std::vector<double>(3, 0);
793  n3[0] = dlambda[1][0];
794  n3[1] = dlambda[1][1];
795  n3[2] = dlambda[1][2];
796 
797  std::vector<std::vector<double> > dsub(3, std::vector<double>(3, 0));
798  for(int p = 0; p < 3; p++) {
799  dsub[0][p] = dlambda[1][p] - dlambda[0][p];
800  dsub[1][p] = dlambda[2][p] - dlambda[1][p];
801  dsub[2][p] = dlambda[0][p] - dlambda[2][p];
802  }
803  // edge-based face functions
804  double dlambda23U = dlambda[0][0] * lambda[1] + dlambda[1][0] * lambda[0];
805  double dlambda23V = dlambda[0][1] * lambda[1] + dlambda[1][1] * lambda[0];
806  double dlambda23W = dlambda[0][2] * lambda[1] + dlambda[1][2] * lambda[0];
807  double prod32 = lambda[0] * lambda[1];
808  for(int i1 = 2; i1 <= _pOrderFace[faceNumber]; i1++) {
809  double dphiU =
810  dlambda23U * OrthogonalPoly::EvalLegendre(i1 - 2, sub[0]) +
811  prod32 * dsub[0][0] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[0]);
812  double dphiV =
813  dlambda23V * OrthogonalPoly::EvalLegendre(i1 - 2, sub[0]) +
814  prod32 * dsub[0][1] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[0]);
815  double dphiW =
816  dlambda23W * OrthogonalPoly::EvalLegendre(i1 - 2, sub[0]) +
817  prod32 * dsub[0][2] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[0]);
818  faceFunctions[faceIt][0] = jacob * (n1[2] * dphiV - n1[1] * dphiW);
819  faceFunctions[faceIt][1] = jacob * (n1[0] * dphiW - n1[2] * dphiU);
820  faceFunctions[faceIt][2] = jacob * (n1[1] * dphiU - n1[0] * dphiV);
821  faceIt++;
822  }
823  double dlambda13U = dlambda[2][0] * lambda[1] + dlambda[1][0] * lambda[2];
824  double dlambda13V = dlambda[2][1] * lambda[1] + dlambda[1][1] * lambda[2];
825  double dlambda13W = dlambda[2][2] * lambda[1] + dlambda[1][2] * lambda[2];
826  double prod13 = lambda[2] * lambda[1];
827  for(int i1 = 2; i1 <= _pOrderFace[faceNumber]; i1++) {
828  double dphiU =
829  dlambda13U * OrthogonalPoly::EvalLegendre(i1 - 2, sub[1]) +
830  prod13 * dsub[1][0] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[1]);
831  double dphiV =
832  dlambda13V * OrthogonalPoly::EvalLegendre(i1 - 2, sub[1]) +
833  prod13 * dsub[1][1] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[1]);
834  double dphiW =
835  dlambda13W * OrthogonalPoly::EvalLegendre(i1 - 2, sub[1]) +
836  prod13 * dsub[1][2] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[1]);
837 
838  faceFunctions[faceIt][0] = jacob * (n2[2] * dphiV - n2[1] * dphiW);
839  faceFunctions[faceIt][1] = jacob * (n2[0] * dphiW - n2[2] * dphiU);
840 
841  faceFunctions[faceIt][2] = jacob * (n2[1] * dphiU - n2[0] * dphiV);
842 
843  faceIt++;
844  }
845  double dlambda12U = dlambda[2][0] * lambda[0] + dlambda[0][0] * lambda[2];
846  double dlambda12V = dlambda[2][1] * lambda[0] + dlambda[0][1] * lambda[2];
847  double dlambda12W = dlambda[2][2] * lambda[0] + dlambda[0][2] * lambda[2];
848  double prod12 = lambda[0] * lambda[2];
849  for(int i1 = 2; i1 <= _pOrderFace[faceNumber]; i1++) {
850  double dphiU =
851  dlambda12U * OrthogonalPoly::EvalLegendre(i1 - 2, sub[2]) +
852  prod12 * dsub[2][0] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[2]);
853  double dphiV =
854  dlambda12V * OrthogonalPoly::EvalLegendre(i1 - 2, sub[2]) +
855  prod12 * dsub[2][1] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[2]);
856  double dphiW =
857  dlambda12W * OrthogonalPoly::EvalLegendre(i1 - 2, sub[2]) +
858  prod12 * dsub[2][2] * OrthogonalPoly::EvalDLegendre(i1 - 2, sub[2]);
859 
860  faceFunctions[faceIt][0] = jacob * (n3[2] * dphiV - n3[1] * dphiW);
861  faceFunctions[faceIt][1] = jacob * (n3[0] * dphiW - n3[2] * dphiU);
862  faceFunctions[faceIt][2] = jacob * (n3[1] * dphiU - n3[0] * dphiV);
863  faceIt++;
864  }
865  // Genuine face functions
866  double subBA = sub[0];
867  double subAC = sub[2];
868  std::vector<double> dsubBA(3);
869  std::vector<double> dsubAC(3);
870  for(int i = 0; i < 3; i++) {
871  dsubBA[i] = dsub[0][i];
872  dsubAC[i] = dsub[2][i];
873  }
874  std::vector<double> LSubAC(_pOrderFace[faceNumber] - 2);
875  std::vector<double> dLSubAC(_pOrderFace[faceNumber] - 2);
876  for(int it = 0; it < _pOrderFace[faceNumber] - 2; it++) {
877  LSubAC[it] = OrthogonalPoly::EvalLegendre(it, subAC);
878  dLSubAC[it] = OrthogonalPoly::EvalDLegendre(it, subAC);
879  }
880  std::vector<double> LSubBA(_pOrderFace[faceNumber] - 2);
881  std::vector<double> dLSubBA(_pOrderFace[faceNumber] - 2);
882  for(int it = 0; it < _pOrderFace[faceNumber] - 2; it++) {
883  LSubBA[it] = OrthogonalPoly::EvalLegendre(it, subBA);
884  dLSubBA[it] = OrthogonalPoly::EvalDLegendre(it, subBA);
885  }
886  std::vector<double> dProduct(3, 0);
887  for(int i = 0; i < 3; i++) {
888  dProduct[i] = dlambda[0][i] * lambda[1] * lambda[2] +
889  dlambda[1][i] * lambda[0] * lambda[2] +
890  dlambda[2][i] * lambda[1] * lambda[0];
891  }
892  double product = lambda[0] * lambda[1] * lambda[2];
893  std::vector<std::vector<double> > copy(
894  int((_pOrderFace[faceNumber] - 2) * (_pOrderFace[faceNumber] - 1) / 2),
895  std::vector<double>(3, 0));
896  int itCopy = 0;
897  for(int n1 = 0; n1 < _pOrderFace[faceNumber] - 2; n1++) {
898  for(int n2 = 0; n2 < _pOrderFace[faceNumber] - 2 - n1; n2++) {
899  std::vector<double> gradFace(3, 0);
900 
901  for(int i = 0; i < 3; i++) {
902  gradFace[i] = dProduct[i] * LSubAC[n2] * LSubBA[n1] +
903  product * dsubBA[i] * LSubAC[n2] * dLSubBA[n1] +
904  product * dsubAC[i] * dLSubAC[n2] * LSubBA[n1];
905  copy[itCopy][i] = gradFace[i];
906  }
907  itCopy++;
908  curlFunction(jacob, dlambda[1], gradFace, faceFunctions[faceIt]);
909  faceIt++;
910  }
911  }
912  itCopy = 0;
913  for(int n1 = 0; n1 < _pOrderFace[faceNumber] - 2; n1++) {
914  for(int n2 = 0; n2 < _pOrderFace[faceNumber] - 2 - n1; n2++) {
915  curlFunction(jacob, dlambda[2], copy[itCopy], faceFunctions[faceIt]);
916  faceIt++;
917  itCopy++;
918  }
919  }
920  }
921  else {
922  throw std::runtime_error("unknown typeFunction");
923  }
924  }
925 }
926 void HierarchicalBasisHcurlTetra::orientFace(
927  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
928  const std::vector<std::vector<double> > &quadFaceFunctionsAllOrientation,
929  const std::vector<std::vector<double> > &triFaceFunctionsAllOrientation,
930  std::vector<std::vector<double> > &fTableCopy)
931 {
932  int iterator = 0;
933  for(int i = 0; i < faceNumber; i++) {
934  iterator +=
935  3 * (_pOrderFace[i] - 1) + (_pOrderFace[i] - 2) * (_pOrderFace[i] - 1);
936  }
937  int numFaceFunctions =
938  3 * (_pOrderFace[faceNumber] - 1) +
939  (_pOrderFace[faceNumber] - 2) * (_pOrderFace[faceNumber] - 1);
940  int iOrientation = numberOrientationTriFace(flag1, flag2);
941  int offset = iOrientation * _nTriFaceFunction;
942  int offset2 = iterator + numFaceFunctions;
943  for(int i = iterator; i < offset2; i++) {
944  fTableCopy[i][0] = triFaceFunctionsAllOrientation[i + offset][0];
945  fTableCopy[i][1] = triFaceFunctionsAllOrientation[i + offset][1];
946  fTableCopy[i][2] = triFaceFunctionsAllOrientation[i + offset][2];
947  }
948 }
949 
950 inline void HierarchicalBasisHcurlTetra::curlFunction(
951  const double &a, const std::vector<double> &nD,
952  const std::vector<double> &grad, std::vector<double> &result)
953 {
954  result[0] = a * (nD[2] * grad[1] - nD[1] * grad[2]);
955  result[1] = a * (nD[0] * grad[2] - nD[2] * grad[0]);
956  result[2] = a * (nD[1] * grad[0] - nD[0] * grad[1]);
957 }
958 inline void HierarchicalBasisHcurlTetra::gradient(
959  const double &lambda1, const double &lambda2,
960  const std::vector<double> &dlambda1, const std::vector<double> &dlambda2,
961  std::vector<double> &result)
962 {
963  for(int i = 0; i < 3; i++) {
964  result[i] = lambda1 * dlambda2[i] + lambda2 * dlambda1[i];
965  }
966 }
967 
968 void HierarchicalBasisHcurlTetra::generateCurlBasis(
969  double const &u, double const &v, double const &w,
970  std::vector<std::vector<double> > &edgeBasis,
971  std::vector<std::vector<double> > &faceBasis,
972  std::vector<std::vector<double> > &bubbleBasis)
973 {
974  //***
975  // to map onto the reference domain of gmsh:
976  double uc = 2 * u - 1;
977  double vc = 2 * v - 1;
978  double wc = 2 * w - 1;
979  double det = 8;
980  double jacob = 0.5;
981  //*****
982  // compute all needed terms
983  double lambda1 = _affineCoordinate(1, uc, vc, wc);
984  double lambda2 = _affineCoordinate(2, uc, vc, wc);
985  double lambda3 = _affineCoordinate(3, uc, vc, wc);
986  double lambda4 = _affineCoordinate(4, uc, vc, wc);
987  std::vector<double> dlambda1(3, 0);
988  std::vector<double> dlambda2(3, -0.5);
989  std::vector<double> dlambda3(3, 0);
990  std::vector<double> dlambda4(3, 0);
991  dlambda1[1] = 0.5;
992  dlambda3[0] = 0.5;
993  dlambda4[2] = 0.5;
994  std::vector<double> n1 = std::vector<double>(3, 0);
995  n1[1] = 1;
996  std::vector<double> n2 = std::vector<double>(3, 0);
997  n2[0] = -1. / sqrt(3);
998  n2[1] = -1. / sqrt(3);
999  n2[2] = -1. / sqrt(3);
1000  std::vector<double> n3 = std::vector<double>(3, 0);
1001  n3[0] = 1;
1002  std::vector<double> n4 = std::vector<double>(3, 0);
1003  n4[2] = 1;
1004  std::vector<double> t1 = std::vector<double>(3, 0);
1005  t1[0] = 1;
1006  std::vector<double> t2 = std::vector<double>(3, 0);
1007  t2[0] = -1;
1008  t2[1] = 1;
1009  std::vector<double> t3 = std::vector<double>(3, 0);
1010  t3[1] = -1;
1011  std::vector<double> t4 = std::vector<double>(3, 0);
1012  t4[2] = 1;
1013  std::vector<double> t5 = std::vector<double>(3, 0);
1014  t5[0] = -1;
1015  t5[2] = 1;
1016  std::vector<double> t6 = std::vector<double>(3, 0);
1017  t6[1] = -1;
1018  t6[2] = 1;
1019  std::vector<double> sub(9, 0);
1020  sub[0] = lambda3 - lambda2;
1021  sub[1] = lambda1 - lambda3;
1022  sub[2] = lambda2 - lambda1;
1023  sub[3] = lambda4 - lambda2;
1024  sub[4] = lambda4 - lambda1;
1025  sub[5] = lambda4 - lambda3;
1026  sub[6] = lambda2 - lambda4;
1027  sub[7] = -lambda2 + lambda1;
1028  sub[8] = lambda3 - lambda4;
1029 
1030  std::vector<std::vector<double> > dsubtraction(9, std::vector<double>(3, 0));
1031  for(int k = 0; k < 3; k++) {
1032  dsubtraction[0][k] = dlambda3[k] - dlambda2[k];
1033  dsubtraction[1][k] = dlambda1[k] - dlambda3[k];
1034  dsubtraction[2][k] = dlambda2[k] - dlambda1[k];
1035  dsubtraction[3][k] = dlambda4[k] - dlambda2[k];
1036  dsubtraction[4][k] = dlambda4[k] - dlambda1[k];
1037  dsubtraction[5][k] = dlambda4[k] - dlambda3[k];
1038  dsubtraction[6][k] = dlambda2[k] - dlambda4[k];
1039  dsubtraction[7][k] = -dlambda2[k] + dlambda1[k];
1040  dsubtraction[8][k] = dlambda3[k] - dlambda4[k];
1041  }
1042  std::vector<std::vector<double> > legendreVector(9);
1043  std::vector<std::vector<double> > dlegendreVector(9);
1044  legendreVector[0] = std::vector<double>(std::max(
1045  std::max(std::max(_pOrderEdge[0], _pOrderFace[0] - 1), _pOrderFace[1] - 1),
1046  0));
1047  legendreVector[1] = std::vector<double>(std::max(
1048  std::max(std::max(_pOrderEdge[1], _pOrderFace[0] - 1), _pOrderFace[3] - 1),
1049  0));
1050  legendreVector[2] = std::vector<double>(std::max(
1051  std::max(std::max(_pOrderEdge[2], _pOrderFace[0] - 1), _pOrderFace[2] - 1),
1052  0));
1053  legendreVector[3] = std::vector<double>(std::max(_pOrderEdge[3], 0));
1054  legendreVector[4] = std::vector<double>(std::max(
1055  std::max(std::max(_pOrderEdge[4], _pOrderFace[2] - 1), _pOrderFace[3] - 1),
1056  0));
1057  legendreVector[5] = std::vector<double>(
1058  std::max(std::max(_pOrderEdge[5], _pOrderFace[1] - 1), 0));
1059  legendreVector[6] = std::vector<double>(
1060  std::max(std::max(_pOrderFace[1] - 1, _pOrderFace[2] - 1), 0));
1061  legendreVector[7] = std::vector<double>(std::max(_pOrderFace[2] - 1, 0));
1062  legendreVector[8] = std::vector<double>(std::max(_pOrderFace[3] - 1, 0));
1063 
1064  dlegendreVector[0] = std::vector<double>(std::max(
1065  std::max(std::max(_pOrderEdge[0], _pOrderFace[0] - 1), _pOrderFace[1] - 1),
1066  0));
1067  dlegendreVector[1] = std::vector<double>(std::max(
1068  std::max(std::max(_pOrderEdge[1], _pOrderFace[0] - 1), _pOrderFace[3] - 1),
1069  0));
1070  dlegendreVector[2] = std::vector<double>(std::max(
1071  std::max(std::max(_pOrderEdge[2], _pOrderFace[0] - 1), _pOrderFace[2] - 1),
1072  0));
1073  dlegendreVector[3] = std::vector<double>(std::max(_pOrderEdge[3], 0));
1074  dlegendreVector[4] = std::vector<double>(std::max(
1075  std::max(std::max(_pOrderEdge[4], _pOrderFace[2] - 1), _pOrderFace[3] - 1),
1076  0));
1077  dlegendreVector[5] = std::vector<double>(
1078  std::max(std::max(_pOrderEdge[5], _pOrderFace[1] - 1), 0));
1079  dlegendreVector[6] = std::vector<double>(
1080  std::max(std::max(_pOrderFace[1] - 1, _pOrderFace[2] - 1), 0));
1081  dlegendreVector[7] = std::vector<double>(std::max(_pOrderFace[2] - 1, 0));
1082  dlegendreVector[8] = std::vector<double>(std::max(_pOrderFace[3] - 1, 0));
1083  for(int i = 0; i < 9; i++) {
1084  for(unsigned int j = 0; j < legendreVector[i].size(); j++) {
1085  legendreVector[i][j] = OrthogonalPoly::EvalLegendre(j, sub[i]);
1086  dlegendreVector[i][j] = OrthogonalPoly::EvalDLegendre(j, sub[i]);
1087  }
1088  }
1089 
1090  std::vector<std::vector<double> > dfaceProduct(_nfaceTri,
1091  std::vector<double>(3));
1092  for(int i = 0; i < 3; i++) {
1093  dfaceProduct[0][i] = lambda3 * lambda2 * dlambda1[i] +
1094  lambda3 * dlambda2[i] * lambda1 +
1095  dlambda3[i] * lambda2 * lambda1;
1096  dfaceProduct[1][i] = lambda3 * lambda2 * dlambda4[i] +
1097  lambda3 * dlambda2[i] * lambda4 +
1098  dlambda3[i] * lambda2 * lambda4;
1099  dfaceProduct[2][i] = lambda1 * lambda2 * dlambda4[i] +
1100  lambda1 * dlambda2[i] * lambda4 +
1101  dlambda1[i] * lambda2 * lambda4;
1102  dfaceProduct[3][i] = lambda1 * lambda3 * dlambda4[i] +
1103  lambda1 * dlambda3[i] * lambda4 +
1104  dlambda1[i] * lambda3 * lambda4;
1105  }
1106  // Whitney functions:
1107  std::vector<std::vector<double> > psie_0(6, std::vector<double>(3, 0));
1108  std::vector<std::vector<double> > psie_1(6, std::vector<double>(3, 0));
1109  for(int i = 0; i < 3; i++) {
1110  psie_0[0][i] = lambda3 * n2[i] / dotProduct(n2, t1) +
1111  lambda2 * n3[i] / dotProduct(n3, t1);
1112  psie_0[1][i] = lambda1 * n3[i] / dotProduct(n3, t2) +
1113  lambda3 * n1[i] / dotProduct(n1, t2);
1114  psie_0[2][i] = lambda2 * n1[i] / dotProduct(n1, t3) +
1115  lambda1 * n2[i] / dotProduct(n2, t3);
1116  psie_0[3][i] = lambda4 * n2[i] / dotProduct(n2, t4) +
1117  lambda2 * n4[i] / dotProduct(n4, t4);
1118  psie_0[4][i] = lambda4 * n1[i] / dotProduct(n1, t6) +
1119  lambda1 * n4[i] / dotProduct(n4, t6);
1120  psie_0[5][i] = lambda4 * n3[i] / dotProduct(n3, t5) +
1121  lambda3 * n4[i] / dotProduct(n4, t5);
1122 
1123  psie_1[0][i] = lambda3 * n2[i] / dotProduct(n2, t1) -
1124  lambda2 * n3[i] / dotProduct(n3, t1);
1125  psie_1[1][i] = lambda1 * n3[i] / dotProduct(n3, t2) -
1126  lambda3 * n1[i] / dotProduct(n1, t2);
1127  psie_1[2][i] = lambda2 * n1[i] / dotProduct(n1, t3) -
1128  lambda1 * n2[i] / dotProduct(n2, t3);
1129  psie_1[3][i] = lambda4 * n2[i] / dotProduct(n2, t4) -
1130  lambda2 * n4[i] / dotProduct(n4, t4);
1131  psie_1[4][i] = lambda4 * n1[i] / dotProduct(n1, t6) -
1132  lambda1 * n4[i] / dotProduct(n4, t6);
1133  psie_1[5][i] = lambda4 * n3[i] / dotProduct(n3, t5) -
1134  lambda3 * n4[i] / dotProduct(n4, t5);
1135  }
1136  // curl of withney functions
1137  std::vector<std::vector<double> > curlpsie_0(6, std::vector<double>(3, 0));
1138  curlpsie_0[0][1] = -1;
1139  curlpsie_0[0][2] = 1;
1140  curlpsie_0[1][2] = 1;
1141  curlpsie_0[2][0] = -1;
1142  curlpsie_0[2][2] = 1;
1143  curlpsie_0[3][0] = -1;
1144  curlpsie_0[3][1] = 1;
1145  curlpsie_0[4][0] = 1;
1146  curlpsie_0[5][1] = -1;
1147  // edge functions
1148  int edgeIt = 0;
1149  for(int i = 0; i < _nedge; i++) {
1150  for(int j = 0; j < 3; j++) {
1151  edgeBasis[edgeIt][j] = det * jacob * curlpsie_0[i][j];
1152  }
1153  edgeIt++;
1154  if(_pOrderEdge[i] >= 1) {
1155  for(int j = 0; j < 3; j++) { edgeBasis[edgeIt][j] = 0; }
1156  edgeIt++;
1157  }
1158  for(int iedge = 2; iedge <= _pOrderEdge[i]; iedge++) {
1159  edgeBasis[edgeIt][0] =
1160  det * jacob *
1161  ((2 * float(iedge) - 1) / float(iedge) *
1162  (dsubtraction[i][1] * dlegendreVector[i][iedge - 1] * psie_1[i][2] -
1163  dsubtraction[i][2] * dlegendreVector[i][iedge - 1] * psie_1[i][1]) -
1164  (float(iedge) - 1) / float(iedge) *
1165  (curlpsie_0[i][0] * legendreVector[i][iedge - 2] +
1166  dsubtraction[i][1] * dlegendreVector[i][iedge - 2] * psie_0[i][2] -
1167  dsubtraction[i][2] * dlegendreVector[i][iedge - 2] * psie_0[i][1]));
1168  edgeBasis[edgeIt][1] =
1169  det * jacob *
1170  ((2 * float(iedge) - 1) / float(iedge) *
1171  (dsubtraction[i][2] * dlegendreVector[i][iedge - 1] * psie_1[i][0] -
1172  dsubtraction[i][0] * dlegendreVector[i][iedge - 1] * psie_1[i][2]) -
1173  (float(iedge) - 1) / float(iedge) *
1174  (curlpsie_0[i][1] * legendreVector[i][iedge - 2] +
1175  dsubtraction[i][2] * dlegendreVector[i][iedge - 2] * psie_0[i][0] -
1176  dsubtraction[i][0] * dlegendreVector[i][iedge - 2] * psie_0[i][2]));
1177  edgeBasis[edgeIt][2] =
1178  det * jacob *
1179  ((2 * float(iedge) - 1) / float(iedge) *
1180  (dsubtraction[i][0] * dlegendreVector[i][iedge - 1] * psie_1[i][1] -
1181  dsubtraction[i][1] * dlegendreVector[i][iedge - 1] * psie_1[i][0]) -
1182  (float(iedge) - 1) / float(iedge) *
1183  (curlpsie_0[i][2] * legendreVector[i][iedge - 2] +
1184  dsubtraction[i][0] * dlegendreVector[i][iedge - 2] * psie_0[i][1] -
1185  dsubtraction[i][1] * dlegendreVector[i][iedge - 2] * psie_0[i][0]));
1186  edgeIt++;
1187  }
1188  }
1189  // edge-based face functions , Genuine face functions
1190  // and Face-based interior functions
1191  int faceIt = 0;
1192  int bubbleIt = 0;
1193  for(int nFace = 0; nFace < _nfaceTri; nFace++) {
1194  int indexVector1 = 0;
1195  int indexVector2 = 0;
1196  std::vector<double> vecTangent1(3, 0);
1197  std::vector<double> vecTangent2(3, 0);
1198  std::vector<double> niT(3, 0);
1199  double faceProduct = 0;
1200  switch(nFace) {
1201  case(0):
1202  indexVector1 = 0;
1203  vecTangent1[0] = 0.5;
1204  indexVector2 = 2;
1205  vecTangent2[1] = 0.5;
1206  niT[2] = 0.5;
1207  faceProduct = lambda1 * lambda2 * lambda3;
1208  for(int i = 0; i < 3; i++) {
1209  double product = 0;
1210  std::vector<double> dproduct(3, 0);
1211  std::vector<double> nD(3, 0);
1212  switch(i) {
1213  case(0):
1214  product = lambda3 * lambda2;
1215  gradient(lambda3, lambda2, dlambda3, dlambda2, dproduct);
1216  nD[1] = 0.5;
1217  break;
1218  case(1):
1219  product = lambda1 * lambda3;
1220  nD[0] = -0.5;
1221  nD[1] = -0.5;
1222  nD[2] = -0.5;
1223  gradient(lambda3, lambda1, dlambda3, dlambda1, dproduct);
1224  break;
1225  case(2):
1226  product = lambda1 * lambda2;
1227  gradient(lambda1, lambda2, dlambda1, dlambda2, dproduct);
1228  nD[0] = 0.5;
1229  break;
1230  }
1231  for(int i1 = 2; i1 <= _pOrderFace[0]; i1++) {
1232  std::vector<double> grad(3, 0);
1233  std::vector<double> gradLegendre(3, 0);
1234  gradLegendre[0] = dsubtraction[i][0] * dlegendreVector[i][i1 - 2];
1235  gradLegendre[1] = dsubtraction[i][1] * dlegendreVector[i][i1 - 2];
1236  gradLegendre[2] = dsubtraction[i][2] * dlegendreVector[i][i1 - 2];
1237 
1238  gradient(product, legendreVector[i][i1 - 2], dproduct, gradLegendre,
1239  grad);
1240  curlFunction(jacob * det, nD, grad, faceBasis[faceIt]);
1241 
1242  faceIt++;
1243  }
1244  }
1245  break;
1246  case(1):
1247  indexVector1 = 0;
1248  vecTangent1[0] = 0.5;
1249  indexVector2 = 6;
1250  vecTangent2[2] = 0.5;
1251  niT[1] = 0.5;
1252  faceProduct = lambda4 * lambda2 * lambda3;
1253 
1254  for(int i = 0; i < 3; i++) {
1255  double product = 0;
1256  std::vector<double> dproduct(3, 0);
1257  std::vector<double> nD(3, 0);
1258  int index = 0;
1259  switch(i) {
1260  case(0):
1261  product = lambda3 * lambda2;
1262  gradient(lambda3, lambda2, dlambda3, dlambda2, dproduct);
1263  nD[2] = 0.5;
1264  index = 0;
1265  break;
1266  case(1):
1267  product = lambda4 * lambda3;
1268  nD[0] = -0.5;
1269  nD[1] = -0.5;
1270  nD[2] = -0.5;
1271  gradient(lambda4, lambda3, dlambda4, dlambda3, dproduct);
1272  index = 5;
1273  break;
1274  case(2):
1275  product = lambda4 * lambda2;
1276  gradient(lambda4, lambda2, dlambda4, dlambda2, dproduct);
1277  nD[0] = 0.5;
1278  index = 6;
1279  break;
1280  }
1281  for(int i1 = 2; i1 <= _pOrderFace[1]; i1++) {
1282  std::vector<double> grad(3, 0);
1283  std::vector<double> gradLegendre(3, 0);
1284  gradLegendre[0] =
1285  dsubtraction[index][0] * dlegendreVector[index][i1 - 2];
1286  gradLegendre[1] =
1287  dsubtraction[index][1] * dlegendreVector[index][i1 - 2];
1288  gradLegendre[2] =
1289  dsubtraction[index][2] * dlegendreVector[index][i1 - 2];
1290  gradient(product, legendreVector[index][i1 - 2], dproduct,
1291  gradLegendre, grad);
1292  curlFunction(jacob * det, nD, grad, faceBasis[faceIt]);
1293  faceIt++;
1294  }
1295  }
1296  break;
1297  case(2):
1298 
1299  indexVector1 = 7;
1300  vecTangent1[1] = 0.5;
1301  indexVector2 = 6;
1302  vecTangent2[2] = 0.5;
1303  niT[0] = 0.5;
1304  faceProduct = lambda1 * lambda2 * lambda4;
1305  for(int i = 0; i < 3; i++) {
1306  double product = 0;
1307  std::vector<double> nD(3, 0);
1308  std::vector<double> dproduct(3, 0);
1309  int index = 0;
1310  switch(i) {
1311  case(0):
1312  product = lambda1 * lambda2;
1313  gradient(lambda1, lambda2, dlambda1, dlambda2, dproduct);
1314  nD[2] = 0.5;
1315  index = 7;
1316  break;
1317  case(1):
1318  product = lambda4 * lambda1;
1319  gradient(lambda4, lambda1, dlambda4, dlambda1, dproduct);
1320  nD[0] = -0.5;
1321  nD[1] = -0.5;
1322  nD[2] = -0.5;
1323  index = 4;
1324  break;
1325  case(2):
1326  product = lambda4 * lambda2;
1327  gradient(lambda4, lambda2, dlambda4, dlambda2, dproduct);
1328  nD[1] = 0.5;
1329  index = 6;
1330  break;
1331  }
1332  for(int i1 = 2; i1 <= _pOrderFace[2]; i1++) {
1333  std::vector<double> grad(3, 0);
1334  std::vector<double> gradLegendre(3, 0);
1335  gradLegendre[0] =
1336  dsubtraction[index][0] * dlegendreVector[index][i1 - 2];
1337  gradLegendre[1] =
1338  dsubtraction[index][1] * dlegendreVector[index][i1 - 2];
1339  gradLegendre[2] =
1340  dsubtraction[index][2] * dlegendreVector[index][i1 - 2];
1341  gradient(product, legendreVector[index][i1 - 2], dproduct,
1342  gradLegendre, grad);
1343  curlFunction(jacob * det, nD, grad, faceBasis[faceIt]);
1344  faceIt++;
1345  }
1346  }
1347  break;
1348  case(3):
1349  indexVector1 = 1;
1350  vecTangent1[1] = 0.5;
1351  indexVector2 = 8;
1352  vecTangent2[2] = 0.5;
1353  niT[0] = -0.5;
1354  niT[1] = -0.5;
1355  niT[2] = -0.5;
1356  faceProduct = lambda1 * lambda4 * lambda3;
1357  for(int i = 0; i < 3; i++) {
1358  std::vector<double> dproduct(3, 0);
1359  double product = 0;
1360  std::vector<double> nD(3, 0);
1361  int index = 0;
1362  switch(i) {
1363  case(0):
1364  product = lambda1 * lambda3;
1365  gradient(lambda3, lambda1, dlambda3, dlambda1, dproduct);
1366  nD[2] = 0.5;
1367  index = 1;
1368  break;
1369  case(1):
1370  product = lambda4 * lambda1;
1371  gradient(lambda4, lambda1, dlambda4, dlambda1, dproduct);
1372  nD[0] = 0.5;
1373  index = 4;
1374  break;
1375  case(2):
1376  product = lambda4 * lambda3;
1377  gradient(lambda4, lambda3, dlambda4, dlambda3, dproduct);
1378  nD[1] = 0.5;
1379  index = 8;
1380  break;
1381  }
1382  for(int i1 = 2; i1 <= _pOrderFace[3]; i1++) {
1383  std::vector<double> grad(3, 0);
1384  std::vector<double> gradLegendre(3, 0);
1385  gradLegendre[0] =
1386  dsubtraction[index][0] * dlegendreVector[index][i1 - 2];
1387  gradLegendre[1] =
1388  dsubtraction[index][1] * dlegendreVector[index][i1 - 2];
1389  gradLegendre[2] =
1390  dsubtraction[index][2] * dlegendreVector[index][i1 - 2];
1391  gradient(product, legendreVector[index][i1 - 2], dproduct,
1392  gradLegendre, grad);
1393  curlFunction(jacob * det, nD, grad, faceBasis[faceIt]);
1394  faceIt++;
1395  }
1396  }
1397  break;
1398  }
1399 
1400  std::vector<std::vector<double> > copy(
1401  int((_pOrderFace[nFace] - 2) * (_pOrderFace[nFace] - 1) / 2),
1402  std::vector<double>(3, 0));
1403  int itCopy = 0;
1404  for(int n1 = 0; n1 < _pOrderFace[nFace] - 2; n1++) {
1405  for(int n2 = 0; n2 < _pOrderFace[nFace] - 2 - n1; n2++) {
1406  std::vector<double> gradFace(3, 0);
1407  for(int i = 0; i < 3; i++) {
1408  gradFace[i] =
1409  dfaceProduct[nFace][i] * legendreVector[indexVector1][n1] *
1410  legendreVector[indexVector2][n2] +
1411  faceProduct * dlegendreVector[indexVector1][n1] *
1412  dsubtraction[indexVector1][i] * legendreVector[indexVector2][n2] +
1413  faceProduct * legendreVector[indexVector1][n1] *
1414  dsubtraction[indexVector2][i] * dlegendreVector[indexVector2][n2];
1415  copy[itCopy][i] = gradFace[i];
1416  }
1417  itCopy++;
1418  curlFunction(jacob * det, vecTangent1, gradFace, faceBasis[faceIt]);
1419 
1420  faceIt++;
1421  }
1422  }
1423  itCopy = 0;
1424  for(int n1 = 0; n1 < _pOrderFace[nFace] - 2; n1++) {
1425  for(int n2 = 0; n2 < _pOrderFace[nFace] - 2 - n1; n2++) {
1426  curlFunction(jacob * det, vecTangent2, copy[itCopy], faceBasis[faceIt]);
1427  faceIt++;
1428  itCopy++;
1429  }
1430  }
1431  for(int n1 = 0; n1 < _pb - 2; n1++) {
1432  for(int n2 = 0; n2 < _pb - 2 - n1; n2++) {
1433  std::vector<double> gradFace(3, 0);
1434  for(int i = 0; i < 3; i++) {
1435  gradFace[i] =
1436  dfaceProduct[nFace][i] * legendreVector[indexVector1][n1] *
1437  legendreVector[indexVector2][n2] +
1438  faceProduct * dlegendreVector[indexVector1][n1] *
1439  dsubtraction[indexVector1][i] * legendreVector[indexVector2][n2] +
1440  faceProduct * legendreVector[indexVector1][n1] *
1441  dsubtraction[indexVector2][i] * dlegendreVector[indexVector2][n2];
1442  }
1443  curlFunction(jacob * det, niT, gradFace, bubbleBasis[bubbleIt]);
1444  bubbleIt++;
1445  }
1446  }
1447  }
1448  if(_pb > 3) {
1449  double l1l2l3l4 = lambda1 * lambda2 * lambda3 * lambda4;
1450  std::vector<double> dl1l2l3l4(3, 0);
1451  for(int i = 0; i < 3; i++) {
1452  dl1l2l3l4[i] = lambda1 * lambda2 * lambda3 * dlambda4[i] +
1453  lambda1 * lambda2 * dlambda3[i] * lambda4 +
1454  lambda1 * dlambda2[i] * lambda3 * lambda4 +
1455  dlambda1[i] * lambda2 * lambda3 * lambda4;
1456  }
1457 
1458  for(int n1 = 0; n1 < _pb - 3; n1++) {
1459  double phi1 = OrthogonalPoly::EvalKernelFunction(n1, sub[7]);
1460  double dphi1 = OrthogonalPoly::EvalDKernelFunction(n1, sub[7]);
1461  for(int n2 = 0; n2 < _pb - 3 - n1; n2++) {
1462  double phi2 = OrthogonalPoly::EvalKernelFunction(n2, sub[0]);
1463  double dphi2 = OrthogonalPoly::EvalDKernelFunction(n2, sub[0]);
1464  for(int n3 = 0; n3 < _pb - 3 - n2 - n1; n3++) {
1465  bubbleBasis[bubbleIt][0] = 0;
1466  bubbleBasis[bubbleIt][1] =
1467  det * jacob *
1468  (dl1l2l3l4[2] * phi1 * phi2 *
1469  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1470  l1l2l3l4 * phi2 * dsubtraction[7][2] * dphi1 *
1471  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1472  l1l2l3l4 * dphi2 * dsubtraction[0][2] * phi1 *
1473  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1474  l1l2l3l4 * phi2 * dsubtraction[3][2] * phi1 *
1475  OrthogonalPoly::EvalDKernelFunction(n3, sub[3]));
1476  bubbleBasis[bubbleIt][2] =
1477  -det * jacob *
1478  (dl1l2l3l4[1] * phi1 * phi2 *
1479  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1480  l1l2l3l4 * phi2 * dsubtraction[7][1] * dphi1 *
1481  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1482  l1l2l3l4 * dphi2 * dsubtraction[0][1] * phi1 *
1483  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1484  l1l2l3l4 * phi2 * dsubtraction[3][1] * phi1 *
1485  OrthogonalPoly::EvalDKernelFunction(n3, sub[3]));
1486 
1487  bubbleIt++;
1488  }
1489  }
1490  }
1491 
1492  for(int n1 = 0; n1 < _pb - 3; n1++) {
1493  double phi1 = OrthogonalPoly::EvalKernelFunction(n1, sub[7]);
1494  double dphi1 = OrthogonalPoly::EvalDKernelFunction(n1, sub[7]);
1495  for(int n2 = 0; n2 < _pb - 3 - n1; n2++) {
1496  double phi2 = OrthogonalPoly::EvalKernelFunction(n2, sub[0]);
1497  double dphi2 = OrthogonalPoly::EvalDKernelFunction(n2, sub[0]);
1498  for(int n3 = 0; n3 < _pb - 3 - n2 - n1; n3++) {
1499  bubbleBasis[bubbleIt][0] =
1500  -det * jacob *
1501  (dl1l2l3l4[2] * phi1 * phi2 *
1502  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1503  l1l2l3l4 * phi2 * dsubtraction[7][2] * dphi1 *
1504  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1505  l1l2l3l4 * dphi2 * dsubtraction[0][2] * phi1 *
1506  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1507  l1l2l3l4 * phi2 * dsubtraction[3][2] * phi1 *
1508  OrthogonalPoly::EvalDKernelFunction(n3, sub[3]));
1509  bubbleBasis[bubbleIt][1] = 0;
1510  bubbleBasis[bubbleIt][2] =
1511  det * jacob *
1512  (dl1l2l3l4[0] * phi1 * phi2 *
1513  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1514  l1l2l3l4 * phi2 * dsubtraction[7][0] * dphi1 *
1515  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1516  l1l2l3l4 * dphi2 * dsubtraction[0][0] * phi1 *
1517  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1518  l1l2l3l4 * phi2 * dsubtraction[3][0] * phi1 *
1519  OrthogonalPoly::EvalDKernelFunction(n3, sub[3]));
1520 
1521  bubbleIt++;
1522  }
1523  }
1524  }
1525 
1526  for(int n1 = 0; n1 < _pb - 3; n1++) {
1527  double phi1 = OrthogonalPoly::EvalKernelFunction(n1, sub[7]);
1528  double dphi1 = OrthogonalPoly::EvalDKernelFunction(n1, sub[7]);
1529  for(int n2 = 0; n2 < _pb - 3 - n1; n2++) {
1530  double phi2 = OrthogonalPoly::EvalKernelFunction(n2, sub[0]);
1531  double dphi2 = OrthogonalPoly::EvalDKernelFunction(n2, sub[0]);
1532  for(int n3 = 0; n3 < _pb - 3 - n2 - n1; n3++) {
1533  bubbleBasis[bubbleIt][0] =
1534  det * jacob *
1535  (dl1l2l3l4[1] * phi1 * phi2 *
1536  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1537  l1l2l3l4 * phi2 * dsubtraction[7][1] * dphi1 *
1538  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1539  l1l2l3l4 * dphi2 * dsubtraction[0][1] * phi1 *
1540  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1541  l1l2l3l4 * phi2 * dsubtraction[3][1] * phi1 *
1542  OrthogonalPoly::EvalDKernelFunction(n3, sub[3]));
1543  bubbleBasis[bubbleIt][1] =
1544  -det * jacob *
1545  (dl1l2l3l4[0] * phi1 * phi2 *
1546  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1547  l1l2l3l4 * phi2 * dsubtraction[7][0] * dphi1 *
1548  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1549  l1l2l3l4 * dphi2 * dsubtraction[0][0] * phi1 *
1550  OrthogonalPoly::EvalKernelFunction(n3, sub[3]) +
1551  l1l2l3l4 * phi2 * dsubtraction[3][0] * phi1 *
1552  OrthogonalPoly::EvalDKernelFunction(n3, sub[3]));
1553  bubbleBasis[bubbleIt][2] = 0;
1554 
1555  bubbleIt++;
1556  }
1557  }
1558  }
1559  }
1560 }
1561 
1562 void HierarchicalBasisHcurlTetra::getKeysInfo(
1563  std::vector<int> &functionTypeInfo, std::vector<int> &orderInfo)
1564 {
1565  int it = 0;
1566  for(int numEdge = 0; numEdge < 6; numEdge++) {
1567  for(int i = 0; i <= _pOrderEdge[numEdge]; i++) {
1568  functionTypeInfo[it] = 1;
1569  orderInfo[it] = i;
1570  it++;
1571  }
1572  }
1573  for(int numFace = 0; numFace < 4; numFace++) {
1574  for(int i = 0; i < 3; i++) {
1575  for(int i1 = 2; i1 <= _pOrderFace[numFace]; i1++) {
1576  functionTypeInfo[it] = 2;
1577  orderInfo[it] = i1;
1578  it++;
1579  }
1580  }
1581  for(int n1 = 1; n1 < _pOrderFace[numFace] - 1; n1++) {
1582  for(int n2 = 1; n2 <= _pOrderFace[numFace] - 1 - n1; n2++) {
1583  functionTypeInfo[it] = 2;
1584  orderInfo[it] = n1 + n2 + 1;
1585  it++;
1586  }
1587  }
1588  for(int n1 = 1; n1 < _pOrderFace[numFace] - 1; n1++) {
1589  for(int n2 = 1; n2 <= _pOrderFace[numFace] - 1 - n1; n2++) {
1590  functionTypeInfo[it] = 2;
1591  orderInfo[it] = n1 + n2 + 1;
1592  it++;
1593  }
1594  }
1595  }
1596  for(int numFace = 0; numFace < 4; numFace++) {
1597  for(int n1 = 1; n1 < _pb - 1; n1++) {
1598  for(int n2 = 1; n2 <= _pb - 1 - n1; n2++) {
1599  functionTypeInfo[it] = 3;
1600  orderInfo[it] = n1 + n2 + 1;
1601  it++;
1602  }
1603  }
1604  }
1605  for(int i = 0; i < 3; i++) {
1606  for(int n1 = 1; n1 < _pb - 2; n1++) {
1607  for(int n2 = 1; n2 <= _pb - 2 - n1; n2++) {
1608  for(int n3 = 1; n3 <= _pb - 1 - n2 - n1; n3++) {
1609  functionTypeInfo[it] = 3;
1610  orderInfo[it] = n1 + n2 + n3 + 1;
1611  it++;
1612  }
1613  }
1614  }
1615  }
1616 }
OrthogonalPoly::EvalDLegendre
double EvalDLegendre(int order, double x)
Definition: OrthogonalPoly.cpp:465
OrthogonalPoly::EvalLegendre
double EvalLegendre(int order, double x)
Definition: OrthogonalPoly.cpp:415
HierarchicalBasis::_nTriFaceFunction
int _nTriFaceFunction
Definition: HierarchicalBasis.h:27
HierarchicalBasisHcurlTetra::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: HierarchicalBasisHcurlTetra.cpp:494
OrthogonalPoly::EvalKernelFunction
double EvalKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:238
HierarchicalBasisHcurlTetra::_pOrderFace
int _pOrderFace[4]
Definition: HierarchicalBasisHcurlTetra.h:95
HierarchicalBasisHcurlTetra::_affineCoordinate
static double _affineCoordinate(const int &j, const double &u, const double &v, const double &w)
Definition: HierarchicalBasisHcurlTetra.cpp:44
HierarchicalBasisHcurlTetra::_pb
int _pb
Definition: HierarchicalBasisHcurlTetra.h:92
HierarchicalBasisHcurlTetra::HierarchicalBasisHcurlTetra
HierarchicalBasisHcurlTetra(int order)
Definition: HierarchicalBasisHcurlTetra.cpp:14
HierarchicalBasis::_nvertex
int _nvertex
Definition: HierarchicalBasis.h:20
HierarchicalBasis::_nBubbleFunction
int _nBubbleFunction
Definition: HierarchicalBasis.h:28
HierarchicalBasisHcurlTetra::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: HierarchicalBasisHcurlTetra.cpp:64
HierarchicalBasisHcurlTetra::orientEdgeFunctionsForNegativeFlag
virtual void orientEdgeFunctionsForNegativeFlag(std::vector< std::vector< double > > &edgeFunctions)
Definition: HierarchicalBasisHcurlTetra.cpp:474
HierarchicalBasisHcurlTetra::dotProduct
static double dotProduct(const std::vector< double > &u, const std::vector< double > &v)
Definition: HierarchicalBasisHcurlTetra.cpp:58
HierarchicalBasisHcurlTetra::getKeysInfo
virtual void getKeysInfo(std::vector< int > &functionTypeInfo, std::vector< int > &orderInfo)
Definition: HierarchicalBasisHcurlTetra.cpp:1562
HierarchicalBasisHcurlTetra::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: HierarchicalBasisHcurlTetra.cpp:443
HierarchicalBasisHcurlTetra::gradient
void gradient(const double &lambda1, const double &lambda2, const std::vector< double > &dlambda1, const std::vector< double > &dlambda2, std::vector< double > &result)
Definition: HierarchicalBasisHcurlTetra.cpp:958
HierarchicalBasisHcurlTetra::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: HierarchicalBasisHcurlTetra.cpp:968
HierarchicalBasisHcurlTetra::getNumberOfOrientations
virtual unsigned int getNumberOfOrientations() const
Definition: HierarchicalBasisHcurlTetra.cpp:39
HierarchicalBasisHcurlTetra::_pOrderEdge
int _pOrderEdge[6]
Definition: HierarchicalBasisHcurlTetra.h:93
HierarchicalBasisHcurlTetra::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: HierarchicalBasisHcurlTetra.cpp:926
HierarchicalBasisHcurlTetra::~HierarchicalBasisHcurlTetra
virtual ~HierarchicalBasisHcurlTetra()
Definition: HierarchicalBasisHcurlTetra.cpp:37
HierarchicalBasisHcurlTetra.h
OrthogonalPoly::EvalDKernelFunction
double EvalDKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:329
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
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
HierarchicalBasisHcurlTetra::curlFunction
void curlFunction(const double &a, const std::vector< double > &nD, const std::vector< double > &grad, std::vector< double > &result)
Definition: HierarchicalBasisHcurlTetra.cpp:950
HierarchicalBasis::numberOrientationTriFace
int numberOrientationTriFace(int const &flag1, int const &flag2)
Definition: HierarchicalBasis.h:131
HierarchicalBasis::_nEdgeFunction
int _nEdgeFunction
Definition: HierarchicalBasis.h:25
HierarchicalBasis::_nedge
int _nedge
Definition: HierarchicalBasis.h:21