gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
HierarchicalBasisH1Pri.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 "HierarchicalBasisH1Pri.h"
13 
14 HierarchicalBasisH1Pri::HierarchicalBasisH1Pri(int order)
15 {
16  _nvertex = 6;
17  _nedge = 9;
18  if(order < 3) { _nfaceTri = 0; }
19  else {
20  _nfaceTri = 2;
21  }
22  _nfaceQuad = 3;
23  _nVertexFunction = 6;
24  _nEdgeFunction = 9 * order - 9;
25  _nQuadFaceFunction = 3 * (order - 1) * (order - 1);
26  _nTriFaceFunction = (order - 2) * (order - 1);
27  _nBubbleFunction = (order - 1) * (order - 2) * (order - 1) / 2;
28  _pb1 = order;
29  _pb2 = order;
30  for(int i = 0; i < 3; i++) {
31  _pOrderQuadFace1[i] = order;
32  _pOrderQuadFace2[i] = order;
33  }
34  for(int i = 0; i < 2; i++) { _pOrderTriFace[i] = order; }
35  for(int i = 0; i < 9; i++) { _pOrderEdge[i] = order; }
36 }
37 
38 HierarchicalBasisH1Pri::~HierarchicalBasisH1Pri() {}
39 
40 unsigned int HierarchicalBasisH1Pri::getNumberOfOrientations() const
41 {
42  return 720; // factorial 6
43 }
44 
45 double HierarchicalBasisH1Pri::_affineCoordinate(const int &j, const double &u,
46  const double &v,
47  const double &w)
48 {
49  switch(j) {
50  case(1): return 0.5 * (1 + v);
51  case(2): return -0.5 * (u + v);
52  case(3): return 0.5 * (1 + u);
53  case(4): return 0.5 * (1 + w);
54  case(5): return 0.5 * (1 - w);
55  default: throw std::runtime_error("j must be : 1<=j<=5");
56  }
57 }
58 void HierarchicalBasisH1Pri::generateBasis(double const &u, double const &v,
59  double const &w,
60  std::vector<double> &vertexBasis,
61  std::vector<double> &edgeBasis,
62  std::vector<double> &faceBasis,
63  std::vector<double> &bubbleBasis)
64 {
65  //***
66  // to map onto the reference domain of gmsh:
67  double uc = 2 * u - 1;
68  double vc = 2 * v - 1;
69  double wc = w;
70  //*****
71  double lambda1 = _affineCoordinate(1, uc, vc, wc);
72  double lambda2 = _affineCoordinate(2, uc, vc, wc);
73  double lambda3 = _affineCoordinate(3, uc, vc, wc);
74  double lambda4 = _affineCoordinate(4, uc, vc, wc);
75  double lambda5 = _affineCoordinate(5, uc, vc, wc);
76  // vertex shape functions:
77  vertexBasis[0] = lambda2 * lambda5;
78  vertexBasis[1] = lambda3 * lambda5;
79  vertexBasis[2] = lambda1 * lambda5;
80  vertexBasis[3] = lambda2 * lambda4;
81  vertexBasis[4] = lambda4 * lambda3;
82  vertexBasis[5] = lambda1 * lambda4;
83  // compute the terms to assemble
84  std::vector<double> product(_nedge);
85  product[0] = vertexBasis[0] * lambda3;
86  product[1] = vertexBasis[0] * lambda1;
87  product[2] = vertexBasis[0] * lambda4;
88  product[3] = vertexBasis[1] * lambda1;
89  product[4] = vertexBasis[1] * lambda4;
90  product[5] = vertexBasis[2] * lambda4;
91  product[6] = vertexBasis[3] * lambda3;
92  product[7] = vertexBasis[3] * lambda1;
93  product[8] = vertexBasis[4] * lambda1;
94  std::vector<double> substraction(5);
95  substraction[0] = lambda3 - lambda2;
96  substraction[1] = lambda1 - lambda2;
97  substraction[2] = lambda4 - lambda5;
98  substraction[3] = lambda1 - lambda3;
99  substraction[4] = lambda2 - lambda1;
100  std::vector<std::vector<double> > phi(5);
101  phi[0] = std::vector<double>(std::max(
102  std::max(std::max(std::max(std::max(_pOrderEdge[0] - 1, _pOrderEdge[6] - 1),
103  _pOrderQuadFace1[0] - 1),
104  _pOrderTriFace[0] - 2),
105  _pOrderTriFace[1] - 2),
106  _pb1 - 2));
107  phi[1] = std::vector<double>(std::max(
108  std::max(_pOrderEdge[1] - 1, _pOrderEdge[7] - 1), _pOrderQuadFace1[1] - 1));
109  phi[2] = std::vector<double>(std::max(
110  std::max(std::max(std::max(std::max(_pOrderEdge[2] - 1, _pOrderEdge[4] - 1),
111  _pOrderEdge[5] - 1),
112  _pOrderQuadFace2[0] - 1),
113  _pOrderQuadFace2[1] - 1),
114  _pOrderQuadFace2[2] - 1));
115 
116  phi[3] = std::vector<double>(std::max(
117  std::max(_pOrderEdge[8] - 1, _pOrderEdge[3] - 1), _pOrderQuadFace1[2] - 1));
118  phi[4] = std::vector<double>(std::max(
119  std::max(std::max(_pOrderTriFace[1] - 2, _pb1 - 2), _pOrderTriFace[0] - 2),
120  0));
121  for(int i = 0; i < 5; i++) {
122  for(unsigned int j = 0; j < phi[i].size(); j++) {
123  phi[i][j] = OrthogonalPoly::EvalKernelFunction(j, substraction[i]);
124  }
125  }
126  // edge shape functions:
127  int indexEdgeBasis = 0;
128  for(int iEdge = 0; iEdge < _nedge; iEdge++) {
129  int numPhi = 0;
130  switch(iEdge) {
131  case(0):
132  case(6): numPhi = 0; break;
133  case(1):
134  case(7): numPhi = 1; break;
135  case(2):
136  case(4):
137  case(5): numPhi = 2; break;
138  case(3):
139  case(8): numPhi = 3; break;
140  }
141  for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] - 1;
142  indexEdgeFunc++) {
143  edgeBasis[indexEdgeBasis] = product[iEdge] * phi[numPhi][indexEdgeFunc];
144 
145  indexEdgeBasis++;
146  }
147  }
148  // face shape functions:
149  int indexFaceFunction = 0;
150  for(int iFace = 0; iFace < _nfaceQuad + _nfaceTri; iFace++) {
151  int indexPhi1 = 0;
152  int indexPhi2 = 0;
153  double lambdaProduct = 0;
154  switch(iFace) {
155  case(3):
156  indexPhi1 = 0;
157  indexPhi2 = 4;
158  lambdaProduct = product[0] * lambda1;
159  break;
160  case(4):
161  indexPhi1 = 0;
162  indexPhi2 = 4;
163  lambdaProduct = product[7] * lambda3;
164  break;
165  case(0):
166  indexPhi1 = 0;
167  indexPhi2 = 2;
168  lambdaProduct = product[0] * lambda4;
169  break;
170  case(1):
171  indexPhi1 = 1;
172  indexPhi2 = 2;
173  lambdaProduct = product[7] * lambda5;
174  break;
175  case(2):
176  indexPhi1 = 3;
177  indexPhi2 = 2;
178  lambdaProduct = product[8] * lambda5;
179  break;
180  }
181  if(iFace < 3) {
182  for(int n1 = 0; n1 < _pOrderQuadFace1[iFace] - 1; n1++) {
183  for(int n2 = 0; n2 < _pOrderQuadFace2[iFace] - 1; n2++) {
184  faceBasis[indexFaceFunction] =
185  lambdaProduct * phi[indexPhi1][n1] * phi[indexPhi2][n2];
186  indexFaceFunction++;
187  }
188  }
189  }
190  else {
191  for(int n1 = 0; n1 < _pOrderTriFace[iFace - 3] - 2; n1++) {
192  for(int n2 = 0; n2 < _pOrderTriFace[iFace - 3] - 2 - n1; n2++) {
193  faceBasis[indexFaceFunction] =
194  lambdaProduct * phi[indexPhi1][n1] * phi[indexPhi2][n2];
195  indexFaceFunction++;
196  }
197  }
198  }
199  }
200  // bubble shape functions:
201  int indexBubbleBasis = 0;
202  double lambdaProductBubble = lambda1 * lambda2 * lambda3;
203  for(int n1 = 0; n1 < _pb1 - 2; n1++) {
204  for(int n2 = 0; n2 < _pb1 - 2 - n1; n2++) {
205  for(int n3 = 2; n3 <= _pb2; n3++) {
206  bubbleBasis[indexBubbleBasis] = lambdaProductBubble * phi[0][n1] *
207  phi[4][n2] *
208  OrthogonalPoly::EvalLobatto(n3, w);
209  indexBubbleBasis++;
210  }
211  }
212  }
213 }
214 void HierarchicalBasisH1Pri::generateGradientBasis(
215  double const &u, double const &v, double const &w,
216  std::vector<std::vector<double> > &gradientVertex,
217  std::vector<std::vector<double> > &gradientEdge,
218  std::vector<std::vector<double> > &gradientFace,
219  std::vector<std::vector<double> > &gradientBubble)
220 {
221  // to map onto the reference domain of gmsh:
222  // ****
223  double uc = 2 * u - 1;
224  double vc = 2 * v - 1;
225  double wc = w;
226  // jacobian=((2,0,0),(0,2,0),(0,0,1))
227  //*******
228  double lambda1 = _affineCoordinate(1, uc, vc, wc);
229  double lambda2 = _affineCoordinate(2, uc, vc, wc);
230  double lambda3 = _affineCoordinate(3, uc, vc, wc);
231  double lambda4 = _affineCoordinate(4, uc, vc, wc);
232  double lambda5 = _affineCoordinate(5, uc, vc, wc);
233  // gradient:
234  std::vector<double> dlambda1(3, 0);
235  std::vector<double> dlambda2(3, 0);
236  std::vector<double> dlambda3(3, 0);
237  std::vector<double> dlambda4(3, 0);
238  std::vector<double> dlambda5(3, 0);
239  dlambda1[1] = 1; // * jacob
240  dlambda2[0] = -1;
241  dlambda2[1] = -1;
242  dlambda3[0] = 1;
243  dlambda4[2] = 0.5;
244  dlambda5[2] = -0.5;
245  for(int i = 0; i < 3; i++) {
246  gradientVertex[0][i] = lambda5 * dlambda2[i] + dlambda5[i] * lambda2;
247  gradientVertex[1][i] = lambda5 * dlambda3[i] + dlambda5[i] * lambda3;
248  gradientVertex[2][i] = lambda5 * dlambda1[i] + dlambda5[i] * lambda1;
249  gradientVertex[3][i] = lambda2 * dlambda4[i] + dlambda2[i] * lambda4;
250  gradientVertex[4][i] = lambda4 * dlambda3[i] + dlambda4[i] * lambda3;
251  gradientVertex[5][i] = lambda1 * dlambda4[i] + dlambda1[i] * lambda4;
252  }
253  // compute the terms to assemble
254  std::vector<double> substraction(5);
255  substraction[0] = lambda3 - lambda2;
256  substraction[1] = lambda1 - lambda2;
257  substraction[2] = lambda4 - lambda5;
258  substraction[3] = lambda1 - lambda3;
259  substraction[4] = lambda2 - lambda1;
260  std::vector<std::vector<double> > dsubstraction(
261  5, std::vector<double>(3, 0)); // = dsubstraction*jacob
262  for(int i = 0; i < 3; i++) {
263  dsubstraction[0][i] = dlambda3[i] - dlambda2[i];
264  dsubstraction[1][i] = dlambda1[i] - dlambda2[i];
265  dsubstraction[2][i] = dlambda4[i] - dlambda5[i];
266  dsubstraction[3][i] = dlambda1[i] - dlambda3[i];
267  dsubstraction[4][i] = dlambda2[i] - dlambda1[i];
268  }
269  std::vector<std::vector<double> > phi(5);
270  std::vector<std::vector<double> > dphi(5);
271  phi[0] = std::vector<double>(std::max(
272  std::max(std::max(std::max(std::max(_pOrderEdge[0] - 1, _pOrderEdge[6] - 1),
273  _pOrderQuadFace1[0] - 1),
274  _pOrderTriFace[0] - 2),
275  _pOrderTriFace[1] - 2),
276  _pb1 - 2));
277  phi[1] = std::vector<double>(std::max(
278  std::max(_pOrderEdge[1] - 1, _pOrderEdge[7] - 1), _pOrderQuadFace1[1] - 1));
279  phi[2] = std::vector<double>(std::max(
280  std::max(std::max(std::max(std::max(_pOrderEdge[2] - 1, _pOrderEdge[4] - 1),
281  _pOrderEdge[5] - 1),
282  _pOrderQuadFace2[0] - 1),
283  _pOrderQuadFace2[1] - 1),
284  _pOrderQuadFace2[2] - 1));
285 
286  phi[3] = std::vector<double>(std::max(
287  std::max(_pOrderEdge[8] - 1, _pOrderEdge[3] - 1), _pOrderQuadFace1[2] - 1));
288  phi[4] = std::vector<double>(std::max(
289  std::max(std::max(_pOrderTriFace[1] - 2, _pb1 - 2), _pOrderTriFace[0] - 2),
290  0));
291  dphi[0] = std::vector<double>(std::max(
292  std::max(std::max(std::max(std::max(_pOrderEdge[0] - 1, _pOrderEdge[6] - 1),
293  _pOrderQuadFace1[0] - 1),
294  _pOrderTriFace[0] - 2),
295  _pOrderTriFace[1] - 2),
296  _pb1 - 2));
297  dphi[1] = std::vector<double>(std::max(
298  std::max(_pOrderEdge[1] - 1, _pOrderEdge[7] - 1), _pOrderQuadFace1[1] - 1));
299  dphi[2] = std::vector<double>(std::max(
300  std::max(std::max(std::max(std::max(_pOrderEdge[2] - 1, _pOrderEdge[4] - 1),
301  _pOrderEdge[5] - 1),
302  _pOrderQuadFace2[0] - 1),
303  _pOrderQuadFace2[1] - 1),
304  _pOrderQuadFace2[2] - 1));
305 
306  dphi[3] = std::vector<double>(std::max(
307  std::max(_pOrderEdge[8] - 1, _pOrderEdge[3] - 1), _pOrderQuadFace1[2] - 1));
308  dphi[4] = std::vector<double>(std::max(
309  std::max(std::max(_pOrderTriFace[1] - 2, _pb1 - 2), _pOrderTriFace[0] - 2),
310  0));
311  for(int i = 0; i < 5; i++) {
312  for(unsigned int j = 0; j < phi[i].size(); j++) {
313  phi[i][j] = OrthogonalPoly::EvalKernelFunction(j, substraction[i]);
314  dphi[i][j] = OrthogonalPoly::EvalDKernelFunction(j, substraction[i]);
315  }
316  }
317  std::vector<double> product1(6);
318  product1[0] = lambda2 * lambda5;
319  product1[1] = lambda3 * lambda5;
320  product1[2] = lambda1 * lambda5;
321  product1[3] = lambda2 * lambda4;
322  product1[4] = lambda4 * lambda3;
323  product1[5] = lambda1 * lambda4;
324  std::vector<double> productEdgeTerm(_nedge);
325  productEdgeTerm[0] = product1[0] * lambda3;
326  productEdgeTerm[1] = product1[0] * lambda1;
327  productEdgeTerm[2] = product1[0] * lambda4;
328  productEdgeTerm[3] = product1[1] * lambda1;
329  productEdgeTerm[4] = product1[1] * lambda4;
330  productEdgeTerm[5] = product1[2] * lambda4;
331  productEdgeTerm[6] = product1[3] * lambda3;
332  productEdgeTerm[7] = product1[3] * lambda1;
333  productEdgeTerm[8] = product1[4] * lambda1;
334  std::vector<std::vector<double> > dproductEdgeTerm(_nedge,
335  std::vector<double>(3));
336  for(int i = 0; i < 3; i++) {
337  dproductEdgeTerm[0][i] =
338  gradientVertex[0][i] * lambda3 + dlambda3[i] * product1[0];
339  dproductEdgeTerm[1][i] =
340  gradientVertex[0][i] * lambda1 + dlambda1[i] * product1[0];
341  dproductEdgeTerm[2][i] =
342  gradientVertex[0][i] * lambda4 + dlambda4[i] * product1[0];
343  dproductEdgeTerm[3][i] =
344  gradientVertex[1][i] * lambda1 + dlambda1[i] * product1[1];
345  dproductEdgeTerm[4][i] =
346  gradientVertex[1][i] * lambda4 + dlambda4[i] * product1[1];
347  dproductEdgeTerm[5][i] =
348  gradientVertex[2][i] * lambda4 + dlambda4[i] * product1[2];
349  dproductEdgeTerm[6][i] =
350  gradientVertex[3][i] * lambda3 + dlambda3[i] * product1[3];
351  dproductEdgeTerm[7][i] =
352  gradientVertex[3][i] * lambda1 + dlambda1[i] * product1[3];
353  dproductEdgeTerm[8][i] =
354  gradientVertex[4][i] * lambda1 + dlambda1[i] * product1[4];
355  }
356 
357  // edge shape functions:
358  int indexEdgeBasis = 0;
359  for(int iEdge = 0; iEdge < _nedge; iEdge++) {
360  int numPhi = 0;
361  switch(iEdge) {
362  case(0):
363  case(6): numPhi = 0; break;
364  case(1):
365  case(7): numPhi = 1; break;
366  case(2):
367  case(4):
368  case(5): numPhi = 2; break;
369  case(3):
370  case(8): numPhi = 3; break;
371  }
372 
373  for(int indexEdgeFunc = 0; indexEdgeFunc < _pOrderEdge[iEdge] - 1;
374  indexEdgeFunc++) {
375  for(int i = 0; i < 3; i++) {
376  gradientEdge[indexEdgeBasis][i] =
377  dproductEdgeTerm[iEdge][i] * phi[numPhi][indexEdgeFunc] +
378  dsubstraction[numPhi][i] * productEdgeTerm[iEdge] *
379  dphi[numPhi][indexEdgeFunc];
380  }
381 
382  indexEdgeBasis++;
383  }
384  }
385  // face shape functions:
386  int indexFaceFunction = 0;
387  for(int iFace = 0; iFace < _nfaceQuad + _nfaceTri; iFace++) {
388  int indexPhi1 = 0;
389  int indexPhi2 = 0;
390  double productFaceTerm = 0;
391  std::vector<double> dproductFaceTerm(3, 0);
392  switch(iFace) {
393  case(3):
394  indexPhi1 = 0;
395  indexPhi2 = 4;
396  productFaceTerm = productEdgeTerm[0] * lambda1;
397  for(int i = 0; i < 3; i++) {
398  dproductFaceTerm[i] =
399  dproductEdgeTerm[0][i] * lambda1 + productEdgeTerm[0] * dlambda1[i];
400  }
401  break;
402  case(4):
403  indexPhi1 = 0;
404  indexPhi2 = 4;
405  productFaceTerm = productEdgeTerm[7] * lambda3;
406  for(int i = 0; i < 3; i++) {
407  dproductFaceTerm[i] =
408  dproductEdgeTerm[7][i] * lambda3 + productEdgeTerm[7] * dlambda3[i];
409  }
410  break;
411  case(0):
412  indexPhi1 = 0;
413  indexPhi2 = 2;
414  productFaceTerm = productEdgeTerm[0] * lambda4;
415  for(int i = 0; i < 3; i++) {
416  dproductFaceTerm[i] =
417  dproductEdgeTerm[0][i] * lambda4 + productEdgeTerm[0] * dlambda4[i];
418  }
419  break;
420  case(1):
421  indexPhi1 = 1;
422  indexPhi2 = 2;
423  productFaceTerm = productEdgeTerm[7] * lambda5;
424  for(int i = 0; i < 3; i++) {
425  dproductFaceTerm[i] =
426  dproductEdgeTerm[7][i] * lambda5 + productEdgeTerm[7] * dlambda5[i];
427  }
428  break;
429  case(2):
430  indexPhi1 = 3;
431  indexPhi2 = 2;
432  productFaceTerm = productEdgeTerm[8] * lambda5;
433  for(int i = 0; i < 3; i++) {
434  dproductFaceTerm[i] =
435  dproductEdgeTerm[8][i] * lambda5 + productEdgeTerm[8] * dlambda5[i];
436  }
437  break;
438  }
439  if(iFace < 3) {
440  for(int n1 = 0; n1 < _pOrderQuadFace1[iFace] - 1; n1++) {
441  for(int n2 = 0; n2 < _pOrderQuadFace2[iFace] - 1; n2++) {
442  for(int i = 0; i < 3; i++) {
443  gradientFace[indexFaceFunction][i] =
444  dproductFaceTerm[i] * phi[indexPhi1][n1] * phi[indexPhi2][n2] +
445  productFaceTerm * dsubstraction[indexPhi1][i] *
446  dphi[indexPhi1][n1] * phi[indexPhi2][n2] +
447  productFaceTerm * dsubstraction[indexPhi2][i] *
448  phi[indexPhi1][n1] * dphi[indexPhi2][n2];
449  }
450  indexFaceFunction++;
451  }
452  }
453  }
454  else {
455  for(int n1 = 0; n1 < _pOrderTriFace[iFace - 3] - 2; n1++) {
456  for(int n2 = 0; n2 < _pOrderTriFace[iFace - 3] - 2 - n1; n2++) {
457  for(int i = 0; i < 3; i++) {
458  gradientFace[indexFaceFunction][i] =
459  dproductFaceTerm[i] * phi[indexPhi1][n1] * phi[indexPhi2][n2] +
460  productFaceTerm * dsubstraction[indexPhi1][i] *
461  dphi[indexPhi1][n1] * phi[indexPhi2][n2] +
462  productFaceTerm * dsubstraction[indexPhi2][i] *
463  phi[indexPhi1][n1] * dphi[indexPhi2][n2];
464  }
465 
466  indexFaceFunction++;
467  }
468  }
469  }
470  }
471  // bubble shape functions:
472  int indexBubbleBasis = 0;
473  double lambdaProductBubble = lambda1 * lambda2 * lambda3;
474  std::vector<double> dlambdaProductBubble(3);
475  for(int i = 0; i < 3; i++) {
476  dlambdaProductBubble[i] = dlambda1[i] * lambda2 * lambda3 +
477  lambda1 * dlambda2[i] * lambda3 +
478  lambda1 * lambda2 * dlambda3[i];
479  }
480 
481  for(int n1 = 0; n1 < _pb1 - 2; n1++) {
482  for(int n2 = 0; n2 < _pb1 - 2 - n1; n2++) {
483  for(int n3 = 2; n3 <= _pb2; n3++) {
484  gradientBubble[indexBubbleBasis][0] =
485  (dlambdaProductBubble[0] * phi[0][n1] * phi[4][n2] +
486  lambdaProductBubble * dsubstraction[0][0] * dphi[0][n1] *
487  phi[4][n2] +
488  lambdaProductBubble * dsubstraction[4][0] * phi[0][n1] *
489  dphi[4][n2]) *
490  OrthogonalPoly::EvalLobatto(n3, w);
491 
492  gradientBubble[indexBubbleBasis][1] =
493  (dlambdaProductBubble[1] * phi[0][n1] * phi[4][n2] +
494  lambdaProductBubble * dsubstraction[0][1] * dphi[0][n1] *
495  phi[4][n2] +
496  lambdaProductBubble * dsubstraction[4][1] * phi[0][n1] *
497  dphi[4][n2]) *
498  OrthogonalPoly::EvalLobatto(n3, w);
499 
500  gradientBubble[indexBubbleBasis][2] =
501  (dlambdaProductBubble[2] * phi[0][n1] * phi[4][n2] +
502  lambdaProductBubble * dsubstraction[0][2] * dphi[0][n1] *
503  phi[4][n2] +
504  lambdaProductBubble * dsubstraction[4][2] * phi[0][n1] *
505  dphi[4][n2]) *
506  OrthogonalPoly::EvalLobatto(n3, w) +
507  OrthogonalPoly::EvalDLobatto(n3, w) * lambdaProductBubble *
508  phi[0][n1] * phi[4][n2];
509 
510  indexBubbleBasis++;
511  }
512  }
513  }
514 }
515 
516 void HierarchicalBasisH1Pri::orientEdge(
517  int const &flagOrientation, int const &edgeNumber,
518  std::vector<double> &edgeFunctions,
519  const std::vector<double> &eTablePositiveFlag,
520  const std::vector<double> &eTableNegativeFlag)
521 {
522  if(flagOrientation == -1) {
523  int constant1 = 0;
524  int constant2 = 0;
525  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
526  constant2 = constant2 - 1;
527  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
528  for(int k = constant1; k <= constant2; k++) {
529  edgeFunctions[k] = eTableNegativeFlag[k];
530  }
531  }
532  else {
533  int constant1 = 0;
534  int constant2 = 0;
535  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
536  constant2 = constant2 - 1;
537  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
538  for(int k = constant1; k <= constant2; k++) {
539  edgeFunctions[k] = eTablePositiveFlag[k];
540  }
541  }
542 }
543 void HierarchicalBasisH1Pri::orientEdge(
544  int const &flagOrientation, int const &edgeNumber,
545  std::vector<std::vector<double> > &edgeFunctions,
546  const std::vector<std::vector<double> > &eTablePositiveFlag,
547  const std::vector<std::vector<double> > &eTableNegativeFlag)
548 {
549  if(flagOrientation == -1) {
550  int constant1 = 0;
551  int constant2 = 0;
552  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
553  constant2 = constant2 - 1;
554  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
555  for(int k = constant1; k <= constant2; k++) {
556  edgeFunctions[k][0] = eTableNegativeFlag[k][0];
557  edgeFunctions[k][1] = eTableNegativeFlag[k][1];
558  edgeFunctions[k][2] = eTableNegativeFlag[k][2];
559  }
560  }
561  else {
562  int constant1 = 0;
563  int constant2 = 0;
564  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
565  constant2 = constant2 - 1;
566  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
567  for(int k = constant1; k <= constant2; k++) {
568  edgeFunctions[k][0] = eTablePositiveFlag[k][0];
569  edgeFunctions[k][1] = eTablePositiveFlag[k][1];
570  edgeFunctions[k][2] = eTablePositiveFlag[k][2];
571  }
572  }
573 }
574 
575 void HierarchicalBasisH1Pri::orientEdgeFunctionsForNegativeFlag(
576  std::vector<double> &edgeFunctions)
577 {
578  int constant1 = 0;
579  int constant2 = 0;
580  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
581  constant2 = 0;
582  constant2 = 0;
583  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
584  constant2 = constant2 - 1;
585  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
586  for(int k = constant1; k <= constant2; k++) {
587  if((k - constant1) % 2 != 0) {
588  edgeFunctions[k] = edgeFunctions[k] * (-1);
589  }
590  }
591  }
592 }
593 void HierarchicalBasisH1Pri::orientEdgeFunctionsForNegativeFlag(
594  std::vector<std::vector<double> > &edgeFunctions)
595 {
596  int constant1 = 0;
597  int constant2 = 0;
598  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
599  constant2 = 0;
600  constant2 = 0;
601  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
602  constant2 = constant2 - 1;
603  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
604  for(int k = constant1; k <= constant2; k++) {
605  if((k - constant1) % 2 != 0) {
606  edgeFunctions[k][0] = edgeFunctions[k][0] * (-1);
607  edgeFunctions[k][1] = edgeFunctions[k][1] * (-1);
608  edgeFunctions[k][2] = edgeFunctions[k][2] * (-1);
609  }
610  }
611  }
612 }
613 
614 void HierarchicalBasisH1Pri::orientOneFace(double const &u, double const &v,
615  double const &w, int const &flag1,
616  int const &flag2, int const &flag3,
617  int const &faceNumber,
618  std::vector<double> &faceBasis)
619 {
620  if(faceNumber < 3) {
621  if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
622  int iterator = 0;
623  for(int i = 0; i < faceNumber; i++) {
624  iterator += (_pOrderQuadFace1[i] - 1) * (_pOrderQuadFace2[i] - 1);
625  }
626  if(flag3 == 1) {
627  for(int it1 = 2; it1 <= _pOrderQuadFace1[faceNumber]; it1++) {
628  for(int it2 = 2; it2 <= _pOrderQuadFace2[faceNumber]; it2++) {
629  int impactFlag1 = 1;
630  int impactFlag2 = 1;
631  if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
632  if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
633  faceBasis[iterator] =
634  faceBasis[iterator] * impactFlag1 * impactFlag2;
635  iterator++;
636  }
637  }
638  }
639  else {
640  // to map onto the reference domain of gmsh:
641  double uc = 2 * u - 1;
642  double vc = 2 * v - 1;
643  double wc = w;
644  //*****
645  double lambdaProduct = 0;
646  double var1 = 0;
647  double var2 = 0;
648  switch(faceNumber) {
649  case(0): {
650  double lambda2 = _affineCoordinate(2, uc, vc, wc);
651  double lambda3 = _affineCoordinate(3, uc, vc, wc);
652  double lambda4 = _affineCoordinate(4, uc, vc, wc);
653  double lambda5 = _affineCoordinate(5, uc, vc, wc);
654  lambdaProduct = lambda2 * lambda3 * lambda4 * lambda5;
655  var1 = lambda3 - lambda2;
656  var2 = lambda4 - lambda5;
657 
658  } break;
659  case(1): {
660  double lambda1 = _affineCoordinate(1, uc, vc, wc);
661  double lambda2 = _affineCoordinate(2, uc, vc, wc);
662  double lambda4 = _affineCoordinate(4, uc, vc, wc);
663  double lambda5 = _affineCoordinate(5, uc, vc, wc);
664  lambdaProduct = lambda2 * lambda1 * lambda4 * lambda5;
665  var1 = lambda1 - lambda2;
666  var2 = lambda4 - lambda5;
667  } break;
668  case(2): {
669  double lambda1 = _affineCoordinate(1, uc, vc, wc);
670  double lambda3 = _affineCoordinate(3, uc, vc, wc);
671  double lambda4 = _affineCoordinate(4, uc, vc, wc);
672  double lambda5 = _affineCoordinate(5, uc, vc, wc);
673  lambdaProduct = lambda1 * lambda3 * lambda4 * lambda5;
674  var1 = lambda1 - lambda3;
675  var2 = lambda4 - lambda5;
676  } break;
677  }
678  std::vector<double> phi1(_pOrderQuadFace1[faceNumber] - 1);
679  std::vector<double> phi2(_pOrderQuadFace2[faceNumber] - 1);
680  for(int it = 0; it <= _pOrderQuadFace1[faceNumber] - 2; it++) {
681  phi1[it] = OrthogonalPoly::EvalKernelFunction(it, var1);
682  }
683  for(int it = 0; it <= _pOrderQuadFace2[faceNumber] - 2; it++) {
684  phi2[it] = OrthogonalPoly::EvalKernelFunction(it, var2);
685  }
686 
687  for(int it1 = 0; it1 <= _pOrderQuadFace2[faceNumber] - 2; it1++) {
688  for(int it2 = 0; it2 <= _pOrderQuadFace1[faceNumber] - 2; it2++) {
689  int impactFlag1 = 1;
690  int impactFlag2 = 1;
691  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
692  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
693  faceBasis[iterator] =
694  lambdaProduct * phi1[it2] * phi2[it1] * impactFlag1 * impactFlag2;
695  iterator++;
696  }
697  }
698  }
699  }
700  }
701  else {
702  int constant = faceNumber - 3;
703  if(!(flag1 == 0 && flag2 == 1)) {
704  // to map onto the reference domain of gmsh:
705  double uc = 2 * u - 1;
706  double vc = 2 * v - 1;
707  double wc = w;
708  //*****
709  int iterator = (_pOrderQuadFace2[0] - 1) * (_pOrderQuadFace1[0] - 1) +
710  (_pOrderQuadFace2[1] - 1) * (_pOrderQuadFace1[1] - 1) +
711  (_pOrderQuadFace2[2] - 1) * (_pOrderQuadFace1[2] - 1);
712  for(int i = 0; i < constant; i++) {
713  iterator += int((_pOrderTriFace[i] - 1) * (_pOrderTriFace[i] - 2) / 2);
714  }
715  std::vector<double> lambda(3);
716  double lambdai = 0;
717  switch(faceNumber) {
718  case(3):
719  lambda[0] = _affineCoordinate(2, uc, vc, wc);
720  lambda[1] = _affineCoordinate(3, uc, vc, wc);
721  lambda[2] = _affineCoordinate(1, uc, vc, wc);
722  lambdai = _affineCoordinate(5, uc, vc, wc);
723  break;
724  case(4):
725  lambda[0] = _affineCoordinate(2, uc, vc, wc);
726  lambda[1] = _affineCoordinate(3, uc, vc, wc);
727  lambda[2] = _affineCoordinate(1, uc, vc, wc);
728  lambdai = _affineCoordinate(4, uc, vc, wc);
729  break;
730  }
731  double product = lambda[0] * lambda[1] * lambda[2] * lambdai;
732  if(flag1 == 1 && flag2 == -1) {
733  double copy = lambda[0];
734  lambda[0] = lambda[1];
735  lambda[1] = copy;
736  }
737  else if(flag1 == 0 && flag2 == -1) {
738  double copy = lambda[2];
739  lambda[2] = lambda[1];
740  lambda[1] = copy;
741  }
742  else if(flag1 == 2 && flag2 == -1) {
743  double copy = lambda[2];
744  lambda[2] = lambda[0];
745  lambda[0] = copy;
746  }
747  else if(flag1 == 1 && flag2 == 1) {
748  double copy = lambda[0];
749  lambda[0] = lambda[1];
750  lambda[1] = lambda[2];
751  lambda[2] = copy;
752  }
753  else if(flag1 == 2 && flag2 == 1) {
754  double copy = lambda[0];
755  lambda[0] = lambda[2];
756  lambda[2] = lambda[1];
757  lambda[1] = copy;
758  }
759  double subs1 = lambda[1] - lambda[0];
760  double subs2 = lambda[0] - lambda[2];
761  std::vector<double> phiSubs2(_pOrderTriFace[constant] - 2);
762  for(int it = 0; it < _pOrderTriFace[constant] - 2; it++) {
763  phiSubs2[it] = OrthogonalPoly::EvalKernelFunction(it, subs2);
764  }
765  for(int n1 = 0; n1 < _pOrderTriFace[constant] - 2; n1++) {
766  double phiSubs1 = OrthogonalPoly::EvalKernelFunction(n1, subs1);
767  for(int n2 = 0; n2 < _pOrderTriFace[constant] - 2 - n1; n2++) {
768  faceBasis[iterator] = product * phiSubs1 * phiSubs2[n2];
769  iterator++;
770  }
771  }
772  }
773  }
774 }
775 void HierarchicalBasisH1Pri::orientOneFace(
776  double const &u, double const &v, double const &w, int const &flag1,
777  int const &flag2, int const &flag3, int const &faceNumber,
778  std::vector<std::vector<double> > &gradientFace, std::string typeFunction)
779 {
780  if(faceNumber < 3) {
781  if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
782  int iterator = 0;
783  for(int i = 0; i < faceNumber; i++) {
784  iterator += (_pOrderQuadFace1[i] - 1) * (_pOrderQuadFace2[i] - 1);
785  }
786  if(flag3 == 1) {
787  for(int it1 = 2; it1 <= _pOrderQuadFace1[faceNumber]; it1++) {
788  for(int it2 = 2; it2 <= _pOrderQuadFace2[faceNumber]; it2++) {
789  int impactFlag1 = 1;
790  int impactFlag2 = 1;
791  if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
792  if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
793  gradientFace[iterator][0] =
794  gradientFace[iterator][0] * impactFlag1 * impactFlag2;
795  gradientFace[iterator][1] =
796  gradientFace[iterator][1] * impactFlag1 * impactFlag2;
797  gradientFace[iterator][2] =
798  gradientFace[iterator][2] * impactFlag1 * impactFlag2;
799  iterator++;
800  }
801  }
802  }
803  else {
804  // to map onto the reference domain of gmsh:
805  double uc = 2 * u - 1;
806  double vc = 2 * v - 1;
807  double wc = w;
808  //*****
809  double lambdaProduct = 0;
810  std::vector<double> dlambdaProduct(3, 0);
811  double var1 = 0;
812  std::vector<double> dvar1(3, 0);
813  double var2 = 0;
814  std::vector<double> dvar2(3, 0);
815  switch(faceNumber) {
816  case(0): {
817  double lambda2 = _affineCoordinate(2, uc, vc, wc);
818  double lambda3 = _affineCoordinate(3, uc, vc, wc);
819  double lambda4 = _affineCoordinate(4, uc, vc, wc);
820  double lambda5 = _affineCoordinate(5, uc, vc, wc);
821  double lambda45 = lambda4 * lambda5;
822  double lambda23 = lambda2 * lambda3;
823  lambdaProduct = lambda23 * lambda45;
824  var1 = lambda3 - lambda2;
825  var2 = lambda4 - lambda5;
826  dlambdaProduct[0] = -lambda45 * var1;
827  dlambdaProduct[1] = -lambda45 * lambda3;
828  dlambdaProduct[2] = -0.5 * lambda23 * var2;
829  dvar1[0] = 2;
830  dvar1[1] = 1;
831  dvar2[2] = 1;
832 
833  } break;
834  case(1): {
835  double lambda1 = _affineCoordinate(1, uc, vc, wc);
836  double lambda2 = _affineCoordinate(2, uc, vc, wc);
837  double lambda4 = _affineCoordinate(4, uc, vc, wc);
838  double lambda5 = _affineCoordinate(5, uc, vc, wc);
839  var1 = lambda1 - lambda2;
840  var2 = lambda4 - lambda5;
841  double lambda45 = lambda4 * lambda5;
842  double lambda21 = lambda2 * lambda1;
843  lambdaProduct = lambda21 * lambda45;
844  dlambdaProduct[0] = -lambda45 * lambda1;
845  dlambdaProduct[1] = lambda45 * (lambda2 - lambda1);
846  dlambdaProduct[2] = -0.5 * lambda21 * var2;
847  dvar1[0] = 1;
848  dvar1[1] = 2;
849  dvar2[2] = 1;
850 
851  } break;
852  case(2): {
853  double lambda1 = _affineCoordinate(1, uc, vc, wc);
854  double lambda3 = _affineCoordinate(3, uc, vc, wc);
855  double lambda4 = _affineCoordinate(4, uc, vc, wc);
856  double lambda5 = _affineCoordinate(5, uc, vc, wc);
857  var1 = lambda1 - lambda3;
858  var2 = lambda4 - lambda5;
859  double lambda45 = lambda4 * lambda5;
860  double lambda13 = lambda1 * lambda3;
861  lambdaProduct = lambda13 * lambda45;
862  dlambdaProduct[0] = lambda45 * lambda1;
863  dlambdaProduct[1] = lambda45 * lambda3;
864  dlambdaProduct[2] = -0.5 * lambda13 * var2;
865  dvar1[0] = -1;
866  dvar1[1] = 1;
867  dvar2[2] = 1;
868 
869  } break;
870  }
871  std::vector<double> phi1(_pOrderQuadFace1[faceNumber] - 1);
872  std::vector<double> phi2(_pOrderQuadFace2[faceNumber] - 1);
873  std::vector<double> dphi1(_pOrderQuadFace1[faceNumber] - 1);
874  std::vector<double> dphi2(_pOrderQuadFace2[faceNumber] - 1);
875  for(int it = 0; it <= _pOrderQuadFace1[faceNumber] - 2; it++) {
876  phi1[it] = OrthogonalPoly::EvalKernelFunction(it, var1);
877  dphi1[it] = OrthogonalPoly::EvalDKernelFunction(it, var1);
878  }
879  for(int it = 0; it <= _pOrderQuadFace2[faceNumber] - 2; it++) {
880  phi2[it] = OrthogonalPoly::EvalKernelFunction(it, var2);
881  dphi2[it] = OrthogonalPoly::EvalDKernelFunction(it, var2);
882  }
883 
884  for(int it1 = 0; it1 <= _pOrderQuadFace2[faceNumber] - 2; it1++) {
885  for(int it2 = 0; it2 <= _pOrderQuadFace1[faceNumber] - 2; it2++) {
886  int impactFlag1 = 1;
887  int impactFlag2 = 1;
888  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
889  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
890  for(int i = 0; i < 3; i++) {
891  gradientFace[iterator][i] =
892  impactFlag1 * impactFlag2 *
893  (dlambdaProduct[i] * phi1[it2] * phi2[it1] +
894  lambdaProduct * dvar1[i] * dphi1[it2] * phi2[it1] +
895  dvar2[i] * lambdaProduct * phi1[it2] * dphi2[it1]);
896  }
897  iterator++;
898  }
899  }
900  }
901  }
902  }
903  else {
904  int constant = faceNumber - 3;
905  if(!(flag1 == 0 && flag2 == 1)) {
906  // to map onto the reference domain of gmsh:
907  double uc = 2 * u - 1;
908  double vc = 2 * v - 1;
909  double wc = w;
910  //*****
911  int iterator = (_pOrderQuadFace2[0] - 1) * (_pOrderQuadFace1[0] - 1) +
912  (_pOrderQuadFace2[1] - 1) * (_pOrderQuadFace1[1] - 1) +
913  (_pOrderQuadFace2[2] - 1) * (_pOrderQuadFace1[2] - 1);
914  for(int i = 0; i < constant; i++) {
915  iterator += int((_pOrderTriFace[i] - 1) * (_pOrderTriFace[i] - 2) / 2);
916  }
917  std::vector<double> lambda(3);
918  std::vector<std::vector<double> > dlambda(3, std::vector<double>(3, 0));
919  double lambdai = 0;
920  std::vector<double> dlambdai(3, 0);
921  switch(faceNumber) {
922  case(3):
923  lambda[0] = _affineCoordinate(2, uc, vc, wc);
924  lambda[1] = _affineCoordinate(3, uc, vc, wc);
925  lambda[2] = _affineCoordinate(1, uc, vc, wc);
926  dlambda[0][0] = -1;
927  dlambda[0][1] = -1;
928  dlambda[1][0] = 1;
929  dlambda[2][1] = 1;
930  lambdai = _affineCoordinate(5, uc, vc, wc);
931  dlambdai[2] = -0.5;
932  break;
933  case(4):
934  lambda[0] = _affineCoordinate(2, uc, vc, wc);
935  lambda[1] = _affineCoordinate(3, uc, vc, wc);
936  lambda[2] = _affineCoordinate(1, uc, vc, wc);
937  dlambda[0][0] = -1;
938  dlambda[0][1] = -1;
939  dlambda[1][0] = 1;
940  dlambda[2][1] = 1;
941  lambdai = _affineCoordinate(4, uc, vc, wc);
942  dlambdai[2] = 0.5;
943  break;
944  }
945  if(flag1 == 1 && flag2 == -1) {
946  double copy = lambda[0];
947  lambda[0] = lambda[1];
948  lambda[1] = copy;
949  std::vector<double> dcopy = dlambda[0];
950  dlambda[0] = dlambda[1];
951  dlambda[1] = dcopy;
952  }
953  else if(flag1 == 0 && flag2 == -1) {
954  double copy = lambda[2];
955  lambda[2] = lambda[1];
956  lambda[1] = copy;
957  std::vector<double> dcopy = dlambda[2];
958  dlambda[2] = dlambda[1];
959  dlambda[1] = dcopy;
960  }
961  else if(flag1 == 2 && flag2 == -1) {
962  double copy = lambda[2];
963  lambda[2] = lambda[0];
964  lambda[0] = copy;
965  std::vector<double> dcopy = dlambda[2];
966  dlambda[2] = dlambda[0];
967  dlambda[0] = dcopy;
968  }
969  else if(flag1 == 1 && flag2 == 1) {
970  double copy = lambda[0];
971  lambda[0] = lambda[1];
972  lambda[1] = lambda[2];
973  lambda[2] = copy;
974  std::vector<double> dcopy = dlambda[0];
975  dlambda[0] = dlambda[1];
976  dlambda[1] = dlambda[2];
977  dlambda[2] = dcopy;
978  }
979  else if(flag1 == 2 && flag2 == 1) {
980  double copy = lambda[0];
981  lambda[0] = lambda[2];
982  lambda[2] = lambda[1];
983  lambda[1] = copy;
984  std::vector<double> dcopy = dlambda[0];
985  dlambda[0] = dlambda[2];
986  dlambda[2] = dlambda[1];
987  dlambda[1] = dcopy;
988  }
989  double subs1 = lambda[1] - lambda[0];
990  double subs2 = lambda[0] - lambda[2];
991  std::vector<double> dsubs1(3, 0);
992  std::vector<double> dsubs2(3, 0);
993  for(int i = 0; i < 3; i++) {
994  dsubs1[i] = dlambda[1][i] - dlambda[0][i];
995  dsubs2[i] = dlambda[0][i] - dlambda[2][i];
996  }
997  double product = lambda[0] * lambda[1] * lambda[2] * lambdai;
998  std::vector<double> dproduct(3, 0);
999  for(int i = 0; i < 3; i++) {
1000  dproduct[i] = dlambda[0][i] * lambda[1] * lambda[2] * lambdai +
1001  lambda[0] * dlambda[1][i] * lambda[2] * lambdai +
1002  lambda[0] * lambda[1] * dlambda[2][i] * lambdai +
1003  lambda[0] * lambda[1] * lambda[2] * dlambdai[i];
1004  }
1005  std::vector<double> phiSubs2(_pOrderTriFace[constant] - 2);
1006  std::vector<double> dphiSubs2(_pOrderTriFace[constant] - 2);
1007  for(int it = 0; it < _pOrderTriFace[constant] - 2; it++) {
1008  phiSubs2[it] = OrthogonalPoly::EvalKernelFunction(it, subs2);
1009  dphiSubs2[it] = OrthogonalPoly::EvalDKernelFunction(it, subs2);
1010  }
1011  for(int n1 = 0; n1 < _pOrderTriFace[constant] - 2; n1++) {
1012  double phiSubs1 = OrthogonalPoly::EvalKernelFunction(n1, subs1);
1013  double dphiSubs1 = OrthogonalPoly::EvalDKernelFunction(n1, subs1);
1014  for(int n2 = 0; n2 < _pOrderTriFace[constant] - 2 - n1; n2++) {
1015  for(int i = 0; i < 3; i++) {
1016  gradientFace[iterator][i] =
1017  dproduct[i] * phiSubs1 * phiSubs2[n2] +
1018  product * dphiSubs1 * dsubs1[i] * phiSubs2[n2] +
1019  product * phiSubs1 * dsubs2[i] * dphiSubs2[n2];
1020  }
1021 
1022  iterator++;
1023  }
1024  }
1025  }
1026  }
1027 }
1028 
1029 void HierarchicalBasisH1Pri::orientFace(
1030  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
1031  const std::vector<double> &quadFaceFunctionsAllOrientation,
1032  const std::vector<double> &triFaceFunctionsAllOrientation,
1033  std::vector<double> &fTableCopy)
1034 {
1035  if(faceNumber < 3) {
1036  int iterator = 0;
1037  for(int i = 0; i < faceNumber; i++) {
1038  iterator += (_pOrderQuadFace1[i] - 1) * (_pOrderQuadFace2[i] - 1);
1039  }
1040  int numFaceFunctions =
1041  (_pOrderQuadFace1[faceNumber] - 1) * (_pOrderQuadFace2[faceNumber] - 1);
1042  int iOrientation = numberOrientationQuadFace(flag1, flag2, flag3);
1043  int offset = iOrientation * _nQuadFaceFunction;
1044  int offset2 = iterator + numFaceFunctions;
1045  for(int i = iterator; i < offset2; i++) {
1046  fTableCopy[i] = quadFaceFunctionsAllOrientation[i + offset];
1047  }
1048  }
1049  else {
1050  int iterator = _nQuadFaceFunction;
1051  int numface = faceNumber - 3;
1052  for(int i = 0; i < numface; i++) {
1053  iterator += int((_pOrderTriFace[i] - 1) * (_pOrderTriFace[i] - 2) / 2);
1054  }
1055  int numFaceFunctions =
1056  int((_pOrderTriFace[numface] - 1) * (_pOrderTriFace[numface] - 2) / 2);
1057  int iOrientation = numberOrientationTriFace(flag1, flag2);
1058  int offset = iOrientation * _nTriFaceFunction - _nQuadFaceFunction;
1059  int offset2 = iterator + numFaceFunctions;
1060  for(int i = iterator; i < offset2; i++) {
1061  fTableCopy[i] = triFaceFunctionsAllOrientation[i + offset];
1062  }
1063  }
1064 }
1065 
1066 void HierarchicalBasisH1Pri::orientFace(
1067  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
1068  const std::vector<std::vector<double> > &quadFaceFunctionsAllOrientation,
1069  const std::vector<std::vector<double> > &triFaceFunctionsAllOrientation,
1070  std::vector<std::vector<double> > &fTableCopy)
1071 {
1072  if(faceNumber < 3) {
1073  int iterator = 0;
1074  for(int i = 0; i < faceNumber; i++) {
1075  iterator += (_pOrderQuadFace1[i] - 1) * (_pOrderQuadFace2[i] - 1);
1076  }
1077  int numFaceFunctions =
1078  (_pOrderQuadFace1[faceNumber] - 1) * (_pOrderQuadFace2[faceNumber] - 1);
1079  int iOrientation = numberOrientationQuadFace(flag1, flag2, flag3);
1080  int offset = iOrientation * _nQuadFaceFunction;
1081  int offset2 = iterator + numFaceFunctions;
1082  for(int i = iterator; i < offset2; i++) {
1083  fTableCopy[i][0] = quadFaceFunctionsAllOrientation[i + offset][0];
1084  fTableCopy[i][1] = quadFaceFunctionsAllOrientation[i + offset][1];
1085  fTableCopy[i][2] = quadFaceFunctionsAllOrientation[i + offset][2];
1086  }
1087  }
1088  else {
1089  int iterator = _nQuadFaceFunction;
1090  int numface = faceNumber - 3;
1091  for(int i = 0; i < numface; i++) {
1092  iterator += int((_pOrderTriFace[i] - 1) * (_pOrderTriFace[i] - 2) / 2);
1093  }
1094  int numFaceFunctions =
1095  int((_pOrderTriFace[numface] - 1) * (_pOrderTriFace[numface] - 2) / 2);
1096  int iOrientation = numberOrientationTriFace(flag1, flag2);
1097  int offset = iOrientation * _nTriFaceFunction - _nQuadFaceFunction;
1098  int offset2 = iterator + numFaceFunctions;
1099  for(int i = iterator; i < offset2; i++) {
1100  fTableCopy[i][0] = triFaceFunctionsAllOrientation[i + offset][0];
1101  fTableCopy[i][1] = triFaceFunctionsAllOrientation[i + offset][1];
1102  fTableCopy[i][2] = triFaceFunctionsAllOrientation[i + offset][2];
1103  }
1104  }
1105 }
1106 
1107 void HierarchicalBasisH1Pri::getKeysInfo(std::vector<int> &functionTypeInfo,
1108  std::vector<int> &orderInfo)
1109 {
1110  for(int i = 0; i < 6; i++) {
1111  functionTypeInfo[i] = 0;
1112  orderInfo[i] = 1;
1113  }
1114  int it = 6;
1115  for(int numEdge = 0; numEdge < 9; numEdge++) {
1116  for(int i = 2; i <= _pOrderEdge[numEdge]; i++) {
1117  functionTypeInfo[it] = 1;
1118  orderInfo[it] = i;
1119  it++;
1120  }
1121  }
1122  for(int iFace = 0; iFace < _nfaceQuad + _nfaceTri; iFace++) {
1123  if(iFace < 3) {
1124  for(int n1 = 2; n1 <= _pOrderQuadFace1[iFace]; n1++) {
1125  for(int n2 = 2; n2 <= _pOrderQuadFace2[iFace]; n2++) {
1126  functionTypeInfo[it] = 2;
1127  orderInfo[it] = std::max(n1, n2);
1128  it++;
1129  }
1130  }
1131  }
1132  else {
1133  for(int n1 = 1; n1 < _pOrderTriFace[iFace - 3] - 1; n1++) {
1134  for(int n2 = 1; n2 <= _pOrderTriFace[iFace - 3] - 1 - n1; n2++) {
1135  functionTypeInfo[it] = 2;
1136  orderInfo[it] = n1 + n2 + 1;
1137  it++;
1138  }
1139  }
1140  }
1141  }
1142  for(int n1 = 1; n1 < _pb1 - 1; n1++) {
1143  for(int n2 = 1; n2 <= _pb1 - 1 - n1; n2++) {
1144  for(int n3 = 2; n3 <= _pb2; n3++) {
1145  functionTypeInfo[it] = 3;
1146  orderInfo[it] = std::max(n1 + n2 + 1, n3);
1147  it++;
1148  }
1149  }
1150  }
1151 }
HierarchicalBasisH1Pri::_pOrderTriFace
int _pOrderTriFace[2]
Definition: HierarchicalBasisH1Pri.h:118
HierarchicalBasisH1Pri::generateBasis
virtual void generateBasis(double const &u, double const &v, double const &w, std::vector< double > &vertexBasis, std::vector< double > &edgeBasis, std::vector< double > &faceBasis, std::vector< double > &bubbleBasis)
Definition: HierarchicalBasisH1Pri.cpp:58
HierarchicalBasis::_nTriFaceFunction
int _nTriFaceFunction
Definition: HierarchicalBasis.h:27
HierarchicalBasisH1Pri::_pOrderQuadFace2
int _pOrderQuadFace2[3]
Definition: HierarchicalBasisH1Pri.h:115
OrthogonalPoly::EvalKernelFunction
double EvalKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:238
HierarchicalBasisH1Pri::_affineCoordinate
static double _affineCoordinate(const int &j, const double &u, const double &v, const double &w)
Definition: HierarchicalBasisH1Pri.cpp:45
HierarchicalBasisH1Pri::generateGradientBasis
void generateGradientBasis(double const &u, double const &v, double const &w, std::vector< std::vector< double > > &gradientVertex, std::vector< std::vector< double > > &gradientEdge, std::vector< std::vector< double > > &gradientFace, std::vector< std::vector< double > > &gradientBubble)
Definition: HierarchicalBasisH1Pri.cpp:214
HierarchicalBasisH1Pri::~HierarchicalBasisH1Pri
virtual ~HierarchicalBasisH1Pri()
Definition: HierarchicalBasisH1Pri.cpp:38
HierarchicalBasisH1Pri::_pOrderEdge
int _pOrderEdge[9]
Definition: HierarchicalBasisH1Pri.h:110
HierarchicalBasis::_nvertex
int _nvertex
Definition: HierarchicalBasis.h:20
HierarchicalBasis::_nBubbleFunction
int _nBubbleFunction
Definition: HierarchicalBasis.h:28
HierarchicalBasisH1Pri::orientFace
virtual void orientFace(int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, const std::vector< double > &quadFaceFunctionsAllOrientation, const std::vector< double > &triFaceFunctionsAllOrientation, std::vector< double > &fTableCopy)
Definition: HierarchicalBasisH1Pri.cpp:1029
HierarchicalBasisH1Pri::orientEdge
virtual void orientEdge(int const &flagOrientation, int const &edgeNumber, std::vector< double > &edgeFunctions, const std::vector< double > &eTablePositiveFlag, const std::vector< double > &eTableNegativeFlag)
Definition: HierarchicalBasisH1Pri.cpp:516
OrthogonalPoly::EvalLobatto
double EvalLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:9
OrthogonalPoly::EvalDLobatto
double EvalDLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:127
HierarchicalBasisH1Pri.h
HierarchicalBasisH1Pri::HierarchicalBasisH1Pri
HierarchicalBasisH1Pri(int order)
Definition: HierarchicalBasisH1Pri.cpp:14
HierarchicalBasisH1Pri::_pOrderQuadFace1
int _pOrderQuadFace1[3]
Definition: HierarchicalBasisH1Pri.h:112
HierarchicalBasisH1Pri::orientOneFace
virtual void orientOneFace(double const &u, double const &v, double const &w, int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, std::vector< double > &faceBasis)
Definition: HierarchicalBasisH1Pri.cpp:614
HierarchicalBasisH1Pri::_pb2
int _pb2
Definition: HierarchicalBasisH1Pri.h:109
HierarchicalBasisH1Pri::getKeysInfo
virtual void getKeysInfo(std::vector< int > &functionTypeInfo, std::vector< int > &orderInfo)
Definition: HierarchicalBasisH1Pri.cpp:1107
HierarchicalBasisH1Pri::_pb1
int _pb1
Definition: HierarchicalBasisH1Pri.h:108
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
HierarchicalBasisH1Pri::orientEdgeFunctionsForNegativeFlag
virtual void orientEdgeFunctionsForNegativeFlag(std::vector< double > &edgeFunctions)
Definition: HierarchicalBasisH1Pri.cpp:575
HierarchicalBasis::numberOrientationTriFace
int numberOrientationTriFace(int const &flag1, int const &flag2)
Definition: HierarchicalBasis.h:131
HierarchicalBasis::_nEdgeFunction
int _nEdgeFunction
Definition: HierarchicalBasis.h:25
HierarchicalBasisH1Pri::getNumberOfOrientations
virtual unsigned int getNumberOfOrientations() const
Definition: HierarchicalBasisH1Pri.cpp:40
HierarchicalBasis::_nedge
int _nedge
Definition: HierarchicalBasis.h:21
HierarchicalBasis::numberOrientationQuadFace
int numberOrientationQuadFace(int const &flag1, int const &flag2, int const &flag3)
Definition: HierarchicalBasis.h:105