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