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