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