gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
HierarchicalBasisH1Tria.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 "HierarchicalBasisH1Tria.h"
13 
14 HierarchicalBasisH1Tria::HierarchicalBasisH1Tria(int pf, int pe0, int pe1,
15  int pe2)
16 {
17  _nvertex = 3;
18  _nedge = 3;
19  _nfaceTri = 1;
20  _nfaceQuad = 0;
21  _nVertexFunction = 3;
22  _nEdgeFunction = pe0 + pe1 + pe2 - 3;
23  _nQuadFaceFunction = 0;
24  _nTriFaceFunction = (pf - 1) * (pf - 2) / 2;
25  _nBubbleFunction = 0;
26  _pf = pf;
27 
28  if(pe0 > pf || pe2 > pf || pe1 > pf) {
29  throw std::runtime_error("pe0, pe1 and pe2 must be <=pf");
30  }
31  _pOrderEdge[0] = pe0;
32  _pOrderEdge[1] = pe1;
33  _pOrderEdge[2] = pe2;
34 }
35 HierarchicalBasisH1Tria::HierarchicalBasisH1Tria(int order)
36 {
37  _nvertex = 3;
38  _nedge = 3;
39  _nfaceTri = 1;
40  _nfaceQuad = 0;
41  _nVertexFunction = 3;
42  _nEdgeFunction = 3 * order - 3;
43  _nQuadFaceFunction = 0;
44  _nTriFaceFunction = (order - 1) * (order - 2) / 2;
45  _nBubbleFunction = 0;
46  _pf = order;
47  _pOrderEdge[0] = order;
48  _pOrderEdge[1] = order;
49  _pOrderEdge[2] = order;
50 }
51 HierarchicalBasisH1Tria::~HierarchicalBasisH1Tria() {}
52 
53 unsigned int HierarchicalBasisH1Tria::getNumberOfOrientations() const
54 {
55  return 6; // factorial 3
56 }
57 
58 double HierarchicalBasisH1Tria::_affineCoordinate(int const &j, double const &u,
59  double const &v)
60 {
61  switch(j) {
62  case(1): return 0.5 * (1 + v);
63  case(2): return -0.5 * (u + v);
64  case(3): return 0.5 * (1 + u);
65  default: throw std::runtime_error("j must be : 1<=j<=3");
66  }
67 }
68 
69 void HierarchicalBasisH1Tria::generateBasis(double const &u, double const &v,
70  double const &w,
71  std::vector<double> &vertexBasis,
72  std::vector<double> &edgeBasis,
73  std::vector<double> &faceBasis,
74  std::vector<double> &bubbleBasis)
75 {
76  //***
77  // to map onto the reference domain of gmsh:
78  double uc = 2 * u - 1;
79  double vc = 2 * v - 1;
80  //*****
81  double lambda1 = _affineCoordinate(1, uc, vc);
82  double lambda2 = _affineCoordinate(2, uc, vc);
83  double lambda3 = _affineCoordinate(3, uc, vc);
84  double product32 = lambda2 * lambda3;
85  double subtraction32 = lambda3 - lambda2;
86  double product13 = lambda3 * lambda1;
87  double subtraction13 = lambda1 - lambda3;
88  double product21 = lambda2 * lambda1;
89  double subtraction21 = lambda2 - lambda1;
90  double product = lambda1 * lambda2 * lambda3;
91  // vertex shape functions:
92  vertexBasis[0] = lambda2;
93  vertexBasis[1] = lambda3;
94  vertexBasis[2] = lambda1;
95  // edge 0 shape functions and a part of face functions :
96  int iterator2 = 0;
97  for(int k = 0; k <= _pOrderEdge[0] - 2; k++) {
98  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction32);
99  edgeBasis[k] = product32 * kernel;
100  int iterator = 0;
101  while(iterator <= _pf - 3 - k) {
102  faceBasis[iterator2 + iterator] = product * kernel;
103  iterator++;
104  }
105  iterator2 = iterator2 + _pf - 2 - k;
106  }
107  for(int k = _pOrderEdge[0] - 1; k <= _pf - 3; k++) {
108  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction32);
109  int iterator = 0;
110  while(iterator <= _pf - 3 - k) {
111  faceBasis[iterator2 + iterator] = product * kernel;
112  iterator++;
113  }
114  iterator2 = iterator2 + _pf - 2 - k;
115  }
116  // edge 1 shape functions :
117  for(int k = 0; k <= _pOrderEdge[1] - 2; k++) {
118  edgeBasis[_pOrderEdge[0] + k - 1] =
119  product13 * OrthogonalPoly::EvalKernelFunction(k, subtraction13);
120  }
121  // edge 2 shape functions and a part of face functions :
122  for(int k = 0; k <= _pOrderEdge[2] - 2; k++) {
123  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction21);
124  edgeBasis[k + _pOrderEdge[0] + _pOrderEdge[1] - 2] = product21 * kernel;
125  int iterator2 = k;
126  int iterator1 = 0;
127  int iterator3 = _pf - 2;
128  while(iterator1 <= _pf - 3 - k) {
129  faceBasis[iterator2] = faceBasis[iterator2] * kernel;
130  iterator2 = iterator2 + iterator3;
131  iterator1++;
132  iterator3--;
133  }
134  }
135  for(int k = _pOrderEdge[2] - 1; k <= _pf - 3; k++) {
136  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction21);
137  int iterator2 = k;
138  int iterator1 = 0;
139  int iterator3 = _pf - 2;
140  while(iterator1 <= _pf - 3 - k) {
141  faceBasis[iterator2] = faceBasis[iterator2] * kernel;
142  iterator2 = iterator2 + iterator3;
143  iterator1++;
144  iterator3--;
145  }
146  }
147 }
148 
149 void HierarchicalBasisH1Tria::generateGradientBasis(
150  double const &u, double const &v, double const &w,
151  std::vector<std::vector<double> > &gradientVertex,
152  std::vector<std::vector<double> > &gradientEdge,
153  std::vector<std::vector<double> > &gradientFace,
154  std::vector<std::vector<double> > &gradientBubble)
155 {
156  // to map onto the reference domain of gmsh:
157  // ****
158  double uc = 2 * u - 1;
159  double vc = 2 * v - 1;
160  double jacob = 2; // jacobian=((2,0),(0,2))
161  //****
162  double dlambda1U = 0;
163  double dlambda1V = 0.5;
164  double dlambda2U = -0.5;
165  double dlambda2V = -0.5;
166  double dlambda3U = 0.5;
167  double dlambda3V = 0;
168  double lambda1 = _affineCoordinate(1, uc, vc);
169  double lambda2 = _affineCoordinate(2, uc, vc);
170  double lambda3 = _affineCoordinate(3, uc, vc);
171  double product32 = lambda2 * lambda3;
172  double subtraction32 = lambda3 - lambda2;
173  double product13 = lambda3 * lambda1;
174  double subtraction13 = lambda1 - lambda3;
175  double product21 = lambda2 * lambda1;
176  double subtraction21 = lambda2 - lambda1;
177  double product = lambda1 * lambda2 * lambda3;
178  // vertex gradient functions:
179  gradientVertex[0][0] = jacob * dlambda2U;
180  gradientVertex[0][1] = jacob * dlambda2V;
181  gradientVertex[1][0] = jacob * dlambda3U;
182  gradientVertex[1][1] = jacob * dlambda3V;
183  gradientVertex[2][0] = jacob * dlambda1U;
184  gradientVertex[2][1] = jacob * dlambda1V;
185  std::vector<double> tablIntermU(_nTriFaceFunction);
186  std::vector<double> tablIntermV(_nTriFaceFunction);
187  // edge 0 gradient and a part of face functions gradient :
188  int iterator2 = 0;
189  for(int k = 0; k <= _pOrderEdge[0] - 2; k++) {
190  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction32);
191  double dKernel = OrthogonalPoly::EvalDKernelFunction(k, subtraction32);
192  gradientEdge[k][0] =
193  (gradientVertex[1][0] * lambda2 + gradientVertex[0][0] * lambda3) *
194  kernel +
195  product32 * (gradientVertex[1][0] - gradientVertex[0][0]) * dKernel;
196  gradientEdge[k][1] =
197  (gradientVertex[1][1] * lambda2 + gradientVertex[0][1] * lambda3) *
198  kernel +
199  product32 * (gradientVertex[1][1] - gradientVertex[0][1]) * dKernel;
200  int iterator = 0;
201  while(iterator <= _pf - 3 - k) {
202  gradientFace[iterator2 + iterator][0] =
203  (gradientVertex[1][0] * product21 + gradientVertex[0][0] * product13 +
204  gradientVertex[2][0] * product32) *
205  kernel +
206  product * (gradientVertex[1][0] - gradientVertex[0][0]) * dKernel;
207  gradientFace[iterator2 + iterator][1] =
208  (gradientVertex[1][1] * product21 + gradientVertex[0][1] * product13 +
209  gradientVertex[2][1] * product32) *
210  kernel +
211  product * (gradientVertex[1][1] - gradientVertex[0][1]) * dKernel;
212  tablIntermU[iterator2 + iterator] =
213  product * (gradientVertex[0][0] - gradientVertex[2][0]) * kernel;
214  tablIntermV[iterator2 + iterator] =
215  product * (gradientVertex[0][1] - gradientVertex[2][1]) * kernel;
216  iterator++;
217  }
218 
219  iterator2 = iterator2 + _pf - 2 - k;
220  }
221  for(int k = _pOrderEdge[0] - 1; k <= _pf - 3; k++) {
222  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction32);
223  double dKernel = OrthogonalPoly::EvalDKernelFunction(k, subtraction32);
224  int iterator = 0;
225  while(iterator <= _pf - 3 - k) {
226  gradientFace[iterator2 + iterator][0] =
227  (gradientVertex[1][0] * product21 + gradientVertex[0][0] * product13 +
228  gradientVertex[2][0] * product32) *
229  kernel +
230  product * (gradientVertex[1][0] - gradientVertex[0][0]) * dKernel;
231  gradientFace[iterator2 + iterator][1] =
232  (gradientVertex[1][1] * product21 + gradientVertex[0][1] * product13 +
233  gradientVertex[2][1] * product32) *
234  kernel +
235  product * (gradientVertex[1][1] - gradientVertex[0][1]) * dKernel;
236  tablIntermU[iterator2 + iterator] =
237  product * (gradientVertex[0][0] - gradientVertex[2][0]) * kernel;
238  tablIntermV[iterator2 + iterator] =
239  product * (gradientVertex[0][1] - gradientVertex[2][1]) * kernel;
240  iterator++;
241  }
242 
243  iterator2 = iterator2 + _pf - 2 - k;
244  }
245  // edge 1 shape functions gradient :
246  for(int k = 0; k <= _pOrderEdge[1] - 2; k++) {
247  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction13);
248  double dKernel = OrthogonalPoly::EvalDKernelFunction(k, subtraction13);
249  gradientEdge[_pOrderEdge[0] + k - 1][0] =
250  (lambda3 * gradientVertex[2][0] + lambda1 * gradientVertex[1][0]) *
251  kernel +
252  product13 * (gradientVertex[2][0] - gradientVertex[1][0]) * dKernel;
253  gradientEdge[_pOrderEdge[0] + k - 1][1] =
254  (lambda3 * gradientVertex[2][1] + lambda1 * gradientVertex[1][1]) *
255  kernel +
256  product13 * (gradientVertex[2][1] - gradientVertex[1][1]) * dKernel;
257  }
258  // edge 2 gradient and a part of face functions gradient :
259  for(int k = 0; k <= _pOrderEdge[2] - 2; k++) {
260  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction21);
261  double dKernel = OrthogonalPoly::EvalDKernelFunction(k, subtraction21);
262  gradientEdge[k + _pOrderEdge[0] + _pOrderEdge[1] - 2][0] =
263  (lambda2 * gradientVertex[2][0] + lambda1 * gradientVertex[0][0]) *
264  kernel +
265  product21 * (gradientVertex[0][0] - gradientVertex[2][0]) * dKernel;
266  gradientEdge[k + _pOrderEdge[0] + _pOrderEdge[1] - 2][1] =
267  (lambda2 * gradientVertex[2][1] + lambda1 * gradientVertex[0][1]) *
268  kernel +
269  product21 * (gradientVertex[0][1] - gradientVertex[2][1]) * dKernel;
270  int iterator2 = k;
271  int iterator1 = 0;
272  int iterator3 = _pf - 2;
273  while(iterator1 <= _pf - 3 - k) {
274  gradientFace[iterator2][0] =
275  gradientFace[iterator2][0] * kernel + tablIntermU[iterator2] * dKernel;
276  gradientFace[iterator2][1] =
277  gradientFace[iterator2][1] * kernel + tablIntermV[iterator2] * dKernel;
278  iterator2 = iterator2 + iterator3;
279  iterator1++;
280  iterator3--;
281  }
282  }
283  for(int k = _pOrderEdge[2] - 1; k <= _pf - 3; k++) {
284  double kernel = OrthogonalPoly::EvalKernelFunction(k, subtraction21);
285  double dKernel = OrthogonalPoly::EvalDKernelFunction(k, subtraction21);
286  int iterator2 = k;
287  int iterator1 = 0;
288  int iterator3 = _pf - 2;
289  while(iterator1 <= _pf - 3 - k) {
290  gradientFace[iterator2][0] =
291  gradientFace[iterator2][0] * kernel + tablIntermU[iterator2] * dKernel;
292  gradientFace[iterator2][1] =
293  gradientFace[iterator2][1] * kernel + tablIntermV[iterator2] * dKernel;
294  iterator2 = iterator2 + iterator3;
295  iterator1++;
296  iterator3--;
297  }
298  }
299 }
300 
301 void HierarchicalBasisH1Tria::orientEdge(
302  int const &flagOrientation, int const &edgeNumber,
303  std::vector<double> &edgeFunctions,
304  const std::vector<double> &eTablePositiveFlag,
305  const std::vector<double> &eTableNegativeFlag)
306 {
307  if(flagOrientation == -1) {
308  int constant1 = 0;
309  int constant2 = 0;
310  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
311  constant2 = constant2 - 1;
312  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
313  for(int k = constant1; k <= constant2; k++) {
314  edgeFunctions[k] = eTableNegativeFlag[k];
315  }
316  }
317  else {
318  int constant1 = 0;
319  int constant2 = 0;
320  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
321  constant2 = constant2 - 1;
322  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
323  for(int k = constant1; k <= constant2; k++) {
324  edgeFunctions[k] = eTablePositiveFlag[k];
325  }
326  }
327 }
328 
329 void HierarchicalBasisH1Tria::orientEdge(
330  int const &flagOrientation, int const &edgeNumber,
331  std::vector<std::vector<double> > &edgeFunctions,
332  const std::vector<std::vector<double> > &eTablePositiveFlag,
333  const std::vector<std::vector<double> > &eTableNegativeFlag)
334 {
335  if(flagOrientation == -1) {
336  int constant1 = 0;
337  int constant2 = 0;
338  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
339  constant2 = constant2 - 1;
340  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
341  for(int k = constant1; k <= constant2; k++) {
342  edgeFunctions[k][0] = eTableNegativeFlag[k][0];
343  edgeFunctions[k][1] = eTableNegativeFlag[k][1];
344  }
345  }
346  else {
347  int constant1 = 0;
348  int constant2 = 0;
349  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
350  constant2 = constant2 - 1;
351  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
352  for(int k = constant1; k <= constant2; k++) {
353  edgeFunctions[k][0] = eTablePositiveFlag[k][0];
354  edgeFunctions[k][1] = eTablePositiveFlag[k][1];
355  }
356  }
357 }
358 void HierarchicalBasisH1Tria::orientEdgeFunctionsForNegativeFlag(
359  std::vector<double> &edgeFunctions)
360 {
361  int constant1 = 0;
362  int constant2 = 0;
363  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
364  constant2 = 0;
365  constant2 = 0;
366  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
367  constant2 = constant2 - 1;
368  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
369  for(int k = constant1; k <= constant2; k++) {
370  if((k - constant1) % 2 != 0) {
371  edgeFunctions[k] = edgeFunctions[k] * (-1);
372  }
373  }
374  }
375 }
376 
377 void HierarchicalBasisH1Tria::orientEdgeFunctionsForNegativeFlag(
378  std::vector<std::vector<double> > &edgeFunctions)
379 {
380  int constant1 = 0;
381  int constant2 = 0;
382  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
383  constant2 = 0;
384  constant2 = 0;
385  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] - 1; }
386  constant2 = constant2 - 1;
387  constant1 = constant2 - _pOrderEdge[edgeNumber] + 2;
388  for(int k = constant1; k <= constant2; k++) {
389  if((k - constant1) % 2 != 0) {
390  edgeFunctions[k][0] = edgeFunctions[k][0] * (-1);
391  edgeFunctions[k][1] = edgeFunctions[k][1] * (-1);
392  edgeFunctions[k][2] = edgeFunctions[k][2] * (-1);
393  }
394  }
395  }
396 }
397 void HierarchicalBasisH1Tria::orientOneFace(double const &u, double const &v,
398  double const &w, int const &flag1,
399  int const &flag2, int const &flag3,
400  int const &faceNumber,
401  std::vector<double> &faceBasis)
402 {
403  if(!(flag1 == 0 && flag2 == 1)) {
404  // to map onto the reference domain of gmsh:
405  double uc = 2 * u - 1;
406  double vc = 2 * v - 1;
407  //*****
408  int iterator = 0;
409  std::vector<double> lambda(3);
410 
411  lambda[0] = _affineCoordinate(2, uc, vc);
412  lambda[1] = _affineCoordinate(3, uc, vc);
413  lambda[2] = _affineCoordinate(1, uc, vc);
414 
415  double product = lambda[0] * lambda[1] * lambda[2];
416  if(flag1 == 1 && flag2 == -1) {
417  double copy = lambda[0];
418  lambda[0] = lambda[1];
419  lambda[1] = copy;
420  }
421  else if(flag1 == 0 && flag2 == -1) {
422  double copy = lambda[2];
423  lambda[2] = lambda[1];
424  lambda[1] = copy;
425  }
426  else if(flag1 == 2 && flag2 == -1) {
427  double copy = lambda[2];
428  lambda[2] = lambda[0];
429  lambda[0] = copy;
430  }
431  else if(flag1 == 1 && flag2 == 1) {
432  double copy = lambda[0];
433  lambda[0] = lambda[1];
434  lambda[1] = lambda[2];
435  lambda[2] = copy;
436  }
437  else if(flag1 == 2 && flag2 == 1) {
438  double copy = lambda[0];
439  lambda[0] = lambda[2];
440  lambda[2] = lambda[1];
441  lambda[1] = copy;
442  }
443  double subs1 = lambda[1] - lambda[0];
444  double subs2 = lambda[0] - lambda[2];
445  std::vector<double> phiSubs2(_pf - 2);
446  for(int it = 0; it < _pf - 2; it++) {
447  phiSubs2[it] = OrthogonalPoly::EvalKernelFunction(it, subs2);
448  }
449  for(int n1 = 0; n1 < _pf - 2; n1++) {
450  double phiSubs1 = OrthogonalPoly::EvalKernelFunction(n1, subs1);
451  for(int n2 = 0; n2 < _pf - 2 - n1; n2++) {
452  faceBasis[iterator] = product * phiSubs1 * phiSubs2[n2];
453  iterator++;
454  }
455  }
456  }
457 }
458 void HierarchicalBasisH1Tria::orientOneFace(
459  double const &u, double const &v, double const &w, int const &flag1,
460  int const &flag2, int const &flag3, int const &faceNumber,
461  std::vector<std::vector<double> > &gradientFace, std::string typeFunction)
462 {
463  if(!(flag1 == 0 && flag2 == 1)) {
464  // to map onto the reference domain of gmsh:
465  double uc = 2 * u - 1;
466  double vc = 2 * v - 1;
467  //*****
468  int iterator = 0;
469  std::vector<double> lambda(3);
470  std::vector<std::vector<double> > dlambda(3, std::vector<double>(2, 0));
471  std::vector<double> dProduct(2); // gradient of (lambdaA*lambdaB*lambdaC)
472 
473  lambda[0] = _affineCoordinate(2, uc, vc);
474  lambda[1] = _affineCoordinate(3, uc, vc);
475  lambda[2] = _affineCoordinate(1, uc, vc);
476  dlambda[0][0] = -1; //* jacobian
477  dlambda[0][1] = -1; //* jacobian
478  dlambda[1][0] = 1; //* jacobian
479  dlambda[2][1] = 1; //* jacobian
480  double pl3l1 = lambda[1] * lambda[2];
481  dProduct[0] = lambda[2] * lambda[0] - pl3l1;
482  dProduct[1] = lambda[0] * lambda[1] - pl3l1;
483  double product = lambda[0] * lambda[1] * lambda[2];
484  if(flag1 == 1 && flag2 == -1) {
485  double copy = lambda[0];
486  lambda[0] = lambda[1];
487  lambda[1] = copy;
488  std::vector<double> dcopy = dlambda[0];
489  dlambda[0] = dlambda[1];
490  dlambda[1] = dcopy;
491  }
492  else if(flag1 == 0 && flag2 == -1) {
493  double copy = lambda[2];
494  lambda[2] = lambda[1];
495  lambda[1] = copy;
496  std::vector<double> dcopy = dlambda[2];
497  dlambda[2] = dlambda[1];
498  dlambda[1] = dcopy;
499  }
500  else if(flag1 == 2 && flag2 == -1) {
501  double copy = lambda[2];
502  lambda[2] = lambda[0];
503  lambda[0] = copy;
504  std::vector<double> dcopy = dlambda[2];
505  dlambda[2] = dlambda[0];
506  dlambda[0] = dcopy;
507  }
508  else if(flag1 == 1 && flag2 == 1) {
509  double copy = lambda[0];
510  lambda[0] = lambda[1];
511  lambda[1] = lambda[2];
512  lambda[2] = copy;
513  std::vector<double> dcopy = dlambda[0];
514  dlambda[0] = dlambda[1];
515  dlambda[1] = dlambda[2];
516  dlambda[2] = dcopy;
517  }
518  else if(flag1 == 2 && flag2 == 1) {
519  double copy = lambda[0];
520  lambda[0] = lambda[2];
521  lambda[2] = lambda[1];
522  lambda[1] = copy;
523  std::vector<double> dcopy = dlambda[0];
524  dlambda[0] = dlambda[2];
525  dlambda[2] = dlambda[1];
526  dlambda[1] = dcopy;
527  }
528  double subsBA = lambda[1] - lambda[0];
529  double subsAC = lambda[0] - lambda[2];
530  std::vector<double> dsubsBA(2);
531  std::vector<double> dsubsAC(2);
532  for(int i = 0; i < 2; i++) {
533  dsubsBA[i] = dlambda[1][i] - dlambda[0][i];
534  dsubsAC[i] = dlambda[0][i] - dlambda[2][i];
535  }
536  std::vector<double> phiSubsAC(_pf - 2);
537  std::vector<double> dphiSubsAC(_pf - 2);
538  for(int it = 0; it < _pf - 2; it++) {
539  phiSubsAC[it] = OrthogonalPoly::EvalKernelFunction(it, subsAC);
540  dphiSubsAC[it] = OrthogonalPoly::EvalDKernelFunction(it, subsAC);
541  }
542  for(int n1 = 0; n1 < _pf - 2; n1++) {
543  double phiBA = OrthogonalPoly::EvalKernelFunction(n1, subsBA);
544  double dphiBA = OrthogonalPoly::EvalDKernelFunction(n1, subsBA);
545  for(int n2 = 0; n2 < _pf - 2 - n1; n2++) {
546  for(int i = 0; i < 2; i++) {
547  gradientFace[iterator][i] =
548  dProduct[i] * phiBA * phiSubsAC[n2] +
549  product * dphiBA * dsubsBA[i] * phiSubsAC[n2] +
550  product * phiBA * dsubsAC[i] * dphiSubsAC[n2];
551  }
552  iterator++;
553  }
554  }
555  }
556 }
557 void HierarchicalBasisH1Tria::orientFace(
558  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
559  const std::vector<double> &quadFaceFunctionsAllOrientation,
560  const std::vector<double> &triFaceFunctionsAllOrientation,
561  std::vector<double> &fTableCopy)
562 {
563  int iOrientation = numberOrientationTriFace(flag1, flag2);
564  int offset = iOrientation * _nTriFaceFunction;
565  for(int i = 0; i < _nTriFaceFunction; i++) {
566  fTableCopy[i] = triFaceFunctionsAllOrientation[i + offset];
567  }
568 }
569 void HierarchicalBasisH1Tria::orientFace(
570  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
571  const std::vector<std::vector<double> > &quadFaceFunctionsAllOrientation,
572  const std::vector<std::vector<double> > &triFaceFunctionsAllOrientation,
573  std::vector<std::vector<double> > &fTableCopy)
574 {
575  int iOrientation = numberOrientationTriFace(flag1, flag2);
576  int offset = iOrientation * _nTriFaceFunction;
577  for(int i = 0; i < _nTriFaceFunction; i++) {
578  fTableCopy[i][0] = triFaceFunctionsAllOrientation[i + offset][0];
579  fTableCopy[i][1] = triFaceFunctionsAllOrientation[i + offset][1];
580  fTableCopy[i][2] = triFaceFunctionsAllOrientation[i + offset][2];
581  }
582 }
583 
584 void HierarchicalBasisH1Tria::getKeysInfo(std::vector<int> &functionTypeInfo,
585  std::vector<int> &orderInfo)
586 {
587  functionTypeInfo[0] = 0;
588  functionTypeInfo[1] = 0;
589  functionTypeInfo[2] = 0;
590  orderInfo[0] = 1;
591  orderInfo[1] = 1;
592  orderInfo[2] = 1;
593  int it = 3;
594  for(int numEdge = 0; numEdge < 3; numEdge++) {
595  for(int i = 2; i <= _pOrderEdge[numEdge]; i++) {
596  functionTypeInfo[it] = 1;
597  orderInfo[it] = i;
598  it++;
599  }
600  }
601  for(int n1 = 1; n1 < _pf - 1; n1++) {
602  for(int n2 = 1; n2 <= _pf - 1 - n1; n2++) {
603  functionTypeInfo[it] = 2;
604  orderInfo[it] = n1 + n2 + 1;
605  it++;
606  }
607  }
608 }
HierarchicalBasisH1Tria::_pf
int _pf
Definition: HierarchicalBasisH1Tria.h:92
HierarchicalBasis::_nTriFaceFunction
int _nTriFaceFunction
Definition: HierarchicalBasis.h:27
OrthogonalPoly::EvalKernelFunction
double EvalKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:238
HierarchicalBasisH1Tria::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: HierarchicalBasisH1Tria.cpp:149
HierarchicalBasisH1Tria::orientEdgeFunctionsForNegativeFlag
virtual void orientEdgeFunctionsForNegativeFlag(std::vector< double > &edgeFunctions)
Definition: HierarchicalBasisH1Tria.cpp:358
HierarchicalBasisH1Tria::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: HierarchicalBasisH1Tria.cpp:69
HierarchicalBasisH1Tria::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: HierarchicalBasisH1Tria.cpp:301
HierarchicalBasisH1Tria::getNumberOfOrientations
virtual unsigned int getNumberOfOrientations() const
Definition: HierarchicalBasisH1Tria.cpp:53
HierarchicalBasis::_nvertex
int _nvertex
Definition: HierarchicalBasis.h:20
HierarchicalBasis::_nBubbleFunction
int _nBubbleFunction
Definition: HierarchicalBasis.h:28
HierarchicalBasisH1Tria::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: HierarchicalBasisH1Tria.cpp:397
HierarchicalBasisH1Tria::_pOrderEdge
int _pOrderEdge[3]
Definition: HierarchicalBasisH1Tria.h:93
HierarchicalBasisH1Tria::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: HierarchicalBasisH1Tria.cpp:557
pe2
const double pe2
Definition: GaussQuadratureQuad.cpp:76
HierarchicalBasisH1Tria::HierarchicalBasisH1Tria
HierarchicalBasisH1Tria(int pf, int pe0, int pe1, int pe2)
Definition: HierarchicalBasisH1Tria.cpp:14
HierarchicalBasisH1Tria::~HierarchicalBasisH1Tria
virtual ~HierarchicalBasisH1Tria()
Definition: HierarchicalBasisH1Tria.cpp:51
HierarchicalBasisH1Tria.h
OrthogonalPoly::EvalDKernelFunction
double EvalDKernelFunction(int order, double x)
Definition: OrthogonalPoly.cpp:329
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
HierarchicalBasis::_nQuadFaceFunction
int _nQuadFaceFunction
Definition: HierarchicalBasis.h:26
HierarchicalBasis::_nfaceTri
int _nfaceTri
Definition: HierarchicalBasis.h:23
HierarchicalBasis::_nfaceQuad
int _nfaceQuad
Definition: HierarchicalBasis.h:22
HierarchicalBasis::_nVertexFunction
int _nVertexFunction
Definition: HierarchicalBasis.h:24
HierarchicalBasisH1Tria::_affineCoordinate
static double _affineCoordinate(int const &j, double const &u, double const &v)
Definition: HierarchicalBasisH1Tria.cpp:58
HierarchicalBasis::numberOrientationTriFace
int numberOrientationTriFace(int const &flag1, int const &flag2)
Definition: HierarchicalBasis.h:131
HierarchicalBasis::_nEdgeFunction
int _nEdgeFunction
Definition: HierarchicalBasis.h:25
HierarchicalBasis::_nedge
int _nedge
Definition: HierarchicalBasis.h:21
HierarchicalBasisH1Tria::getKeysInfo
virtual void getKeysInfo(std::vector< int > &functionTypeInfo, std::vector< int > &orderInfo)
Definition: HierarchicalBasisH1Tria.cpp:584