gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
HierarchicalBasisHcurlQuad.cpp
Go to the documentation of this file.
1 // Gmsh - Copyright (C) 1997-2022 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5 //
6 // Contributed by Ismail Badia.
7 
8 #include <stdexcept>
9 #include "HierarchicalBasisHcurlQuad.h"
10 
11 HierarchicalBasisHcurlQuad::HierarchicalBasisHcurlQuad(int order)
12 
13 {
14  _nvertex = 4;
15  _nedge = 4;
16  _nfaceQuad = 1;
17  _nfaceTri = 0;
18  _nVertexFunction = 0;
19  _nEdgeFunction = 4 * order + 4;
20  _nQuadFaceFunction = 2 * order * (order + 1);
21  _nTriFaceFunction = 0;
22  _nBubbleFunction = 0;
23  _pf1 = order;
24  _pf2 = order;
25  _pOrderEdge[0] = order;
26  _pOrderEdge[1] = order;
27  _pOrderEdge[2] = order;
28  _pOrderEdge[3] = order;
29 }
30 
31 HierarchicalBasisHcurlQuad::~HierarchicalBasisHcurlQuad() {}
32 
33 unsigned int HierarchicalBasisHcurlQuad::getNumberOfOrientations() const
34 {
35  return 24; // factorial 4
36 }
37 
38 double HierarchicalBasisHcurlQuad::_affineCoordinate(int const &j,
39  double const &u,
40  double const &v)
41 {
42  switch(j) {
43  case(1): return 0.5 * (1 + u);
44  case(2): return 0.5 * (1 - u);
45  case(3): return 0.5 * (1 + v);
46  case(4): return 0.5 * (1 - v);
47  default: throw std::runtime_error("j must be : 1<=j<=4");
48  }
49 }
50 
51 void HierarchicalBasisHcurlQuad::generateHcurlBasis(
52  double const &u, double const &v, double const &w,
53  std::vector<std::vector<double> > &edgeBasis,
54  std::vector<std::vector<double> > &faceBasis,
55  std::vector<std::vector<double> > &bubbleBasis)
56 {
57  double lambda1 = _affineCoordinate(1, u, v);
58  double lambda2 = _affineCoordinate(2, u, v);
59  double lambda3 = _affineCoordinate(3, u, v);
60  double lambda4 = _affineCoordinate(4, u, v);
61  std::vector<std::vector<double> > legendreVector(2);
62  legendreVector[0] = std::vector<double>(_pf1 + 1);
63  legendreVector[1] = std::vector<double>(_pf2 + 1);
64  for(unsigned int k = 0; k < legendreVector[0].size(); k++) {
65  legendreVector[0][k] = OrthogonalPoly::EvalLegendre(k, u);
66  }
67  for(unsigned int k = 0; k < legendreVector[1].size(); k++) {
68  legendreVector[1][k] = OrthogonalPoly::EvalLegendre(k, v);
69  }
70  int edgeIt = 0;
71  for(int i = 0; i < _nedge; i++) {
72  int uv = 0;
73  double lambda = 0;
74  std::vector<double> direction(3, 0);
75  switch(i) {
76  case(0):
77  lambda = lambda4;
78  uv = 0;
79  direction[0] = 1;
80  break;
81  case(1):
82  lambda = lambda1;
83  uv = 1;
84  direction[1] = 1;
85  break;
86  case(2):
87  lambda = lambda3;
88  uv = 0;
89  direction[0] = 1;
90  break;
91  case(3):
92  lambda = lambda2;
93  uv = 1;
94  direction[1] = 1;
95  break;
96  }
97  for(int iedge = 0; iedge < _pOrderEdge[i] + 1; iedge++) {
98  for(int j = 0; j < 3; j++) {
99  edgeBasis[edgeIt][j] =
100  lambda * legendreVector[uv][iedge] * direction[j];
101  }
102  edgeIt++;
103  }
104  }
105  int faceIt = 0;
106  for(int n1 = 0; n1 <= _pf1; n1++) {
107  for(int n2 = 2; n2 <= _pf2 + 1; n2++) {
108  faceBasis[faceIt][0] =
109  legendreVector[0][n1] * OrthogonalPoly::EvalLobatto(n2, v);
110  faceBasis[faceIt][1] = 0;
111  faceBasis[faceIt][2] = 0;
112  faceIt++;
113  }
114  }
115  for(int n1 = 2; n1 <= _pf1 + 1; n1++) {
116  for(int n2 = 0; n2 <= _pf2; n2++) {
117  faceBasis[faceIt][0] = 0;
118  faceBasis[faceIt][1] =
119  legendreVector[1][n2] * OrthogonalPoly::EvalLobatto(n1, u);
120  faceBasis[faceIt][2] = 0;
121  faceIt++;
122  }
123  }
124 }
125 
126 void HierarchicalBasisHcurlQuad::orientEdge(
127  int const &flagOrientation, int const &edgeNumber,
128  std::vector<std::vector<double> > &edgeFunctions,
129  const std::vector<std::vector<double> > &eTablePositiveFlag,
130  const std::vector<std::vector<double> > &eTableNegativeFlag)
131 {
132  if(flagOrientation == -1) {
133  int constant1 = 0;
134  int constant2 = 0;
135  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
136  constant2 = constant2 - 1;
137  constant1 = constant2 - _pOrderEdge[edgeNumber];
138  for(int k = constant1; k <= constant2; k++) {
139  edgeFunctions[k][0] = eTableNegativeFlag[k][0];
140  edgeFunctions[k][1] = eTableNegativeFlag[k][1];
141  edgeFunctions[k][2] = eTableNegativeFlag[k][2];
142  }
143  }
144  else {
145  int constant1 = 0;
146  int constant2 = 0;
147  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
148  constant2 = constant2 - 1;
149  constant1 = constant2 - _pOrderEdge[edgeNumber];
150  for(int k = constant1; k <= constant2; k++) {
151  edgeFunctions[k][0] = eTablePositiveFlag[k][0];
152  edgeFunctions[k][1] = eTablePositiveFlag[k][1];
153  edgeFunctions[k][2] = eTablePositiveFlag[k][2];
154  }
155  }
156 }
157 void HierarchicalBasisHcurlQuad::orientEdgeFunctionsForNegativeFlag(
158  std::vector<std::vector<double> > &edgeFunctions)
159 {
160  int constant1 = 0;
161  int constant2 = 0;
162  for(int edgeNumber = 0; edgeNumber < _nedge; edgeNumber++) {
163  constant2 = 0;
164  constant2 = 0;
165  for(int i = 0; i <= edgeNumber; i++) { constant2 += _pOrderEdge[i] + 1; }
166  constant2 = constant2 - 1;
167  constant1 = constant2 - _pOrderEdge[edgeNumber];
168  for(int k = constant1; k <= constant2; k++) {
169  if((k - constant1) % 2 == 0) {
170  edgeFunctions[k][0] = edgeFunctions[k][0] * (-1);
171  edgeFunctions[k][1] = edgeFunctions[k][1] * (-1);
172  edgeFunctions[k][2] = edgeFunctions[k][2] * (-1);
173  }
174  }
175  }
176 }
177 
178 void HierarchicalBasisHcurlQuad::generateCurlBasis(
179  double const &u, double const &v, double const &w,
180  std::vector<std::vector<double> > &edgeBasis,
181  std::vector<std::vector<double> > &faceBasis,
182  std::vector<std::vector<double> > &bubbleBasis)
183 {
184  double dlambda1 = 0.5;
185  double dlambda2 = -0.5;
186  double dlambda3 = 0.5;
187  double dlambda4 = -0.5;
188  std::vector<std::vector<double> > legendreVector(2);
189  legendreVector[0] = std::vector<double>(_pf1 + 1);
190  legendreVector[1] = std::vector<double>(_pf2 + 1);
191  for(unsigned int k = 0; k < legendreVector[0].size(); k++) {
192  legendreVector[0][k] = OrthogonalPoly::EvalLegendre(k, u);
193  }
194  for(unsigned int k = 0; k < legendreVector[1].size(); k++) {
195  legendreVector[1][k] = OrthogonalPoly::EvalLegendre(k, v);
196  }
197  int edgeIt = 0;
198  for(int i = 0; i < _nedge; i++) {
199  int uv = 0;
200  double dlambda = 0;
201  std::vector<double> direction(3, 0);
202  direction[2] = 1;
203  switch(i) {
204  case(0):
205  dlambda = -dlambda4;
206  uv = 0;
207  break;
208  case(1):
209  dlambda = dlambda1;
210  uv = 1;
211  break;
212  case(2):
213  dlambda = -dlambda3;
214  uv = 0;
215  break;
216  case(3):
217  dlambda = dlambda2;
218  uv = 1;
219  break;
220  }
221  for(int iedge = 0; iedge < _pOrderEdge[i] + 1; iedge++) {
222  for(int j = 0; j < 3; j++) {
223  edgeBasis[edgeIt][j] =
224  dlambda * legendreVector[uv][iedge] * direction[j];
225  }
226  edgeIt++;
227  }
228  }
229  int faceIt = 0;
230  for(int n1 = 0; n1 <= _pf1; n1++) {
231  for(int n2 = 2; n2 <= _pf2 + 1; n2++) {
232  faceBasis[faceIt][0] = 0;
233  faceBasis[faceIt][1] = 0;
234  faceBasis[faceIt][2] =
235  -legendreVector[0][n1] * OrthogonalPoly::EvalDLobatto(n2, v);
236  faceIt++;
237  }
238  }
239  for(int n1 = 2; n1 <= _pf1 + 1; n1++) {
240  for(int n2 = 0; n2 <= _pf2; n2++) {
241  faceBasis[faceIt][0] = 0;
242  faceBasis[faceIt][1] = 0;
243  faceBasis[faceIt][2] =
244  legendreVector[1][n2] * OrthogonalPoly::EvalDLobatto(n1, u);
245  faceIt++;
246  }
247  }
248 }
249 void HierarchicalBasisHcurlQuad::orientOneFace(
250  double const &u, double const &v, double const &w, int const &flag1,
251  int const &flag2, int const &flag3, int const &faceNumber,
252  std::vector<std::vector<double> > &faceFunctions, std::string typeFunction)
253 {
254  if(!(flag1 == 1 && flag2 == 1 && flag3 == 1)) {
255  if(flag3 == 1) {
256  int iterator = 0;
257  for(int it1 = 0; it1 <= _pf1; it1++) {
258  for(int it2 = 2; it2 <= _pf2 + 1; it2++) {
259  int impactFlag1 = 1;
260  int impactFlag2 = 1;
261  if(flag1 == -1 && it1 % 2 == 0) { impactFlag1 = -1; }
262  if(flag2 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
263  faceFunctions[iterator][0] =
264  faceFunctions[iterator][0] * impactFlag1 * impactFlag2;
265  faceFunctions[iterator][1] =
266  faceFunctions[iterator][1] * impactFlag1 * impactFlag2;
267  faceFunctions[iterator][2] =
268  faceFunctions[iterator][2] * impactFlag1 * impactFlag2;
269  iterator++;
270  }
271  }
272  for(int it1 = 2; it1 <= _pf1 + 1; it1++) {
273  for(int it2 = 0; it2 <= _pf2; it2++) {
274  int impactFlag1 = 1;
275  int impactFlag2 = 1;
276  if(flag1 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
277  if(flag2 == -1 && it2 % 2 == 0) { impactFlag2 = -1; }
278  faceFunctions[iterator][0] =
279  faceFunctions[iterator][0] * impactFlag1 * impactFlag2;
280  faceFunctions[iterator][1] =
281  faceFunctions[iterator][1] * impactFlag1 * impactFlag2;
282  faceFunctions[iterator][2] =
283  faceFunctions[iterator][2] * impactFlag1 * impactFlag2;
284  iterator++;
285  }
286  }
287  }
288  else {
289  if(typeFunction == "HcurlLegendre") {
290  std::vector<std::vector<double> > legendreVector(2);
291  legendreVector[0] = std::vector<double>(_pf1 + 1);
292  legendreVector[1] = std::vector<double>(_pf2 + 1);
293  for(unsigned int k = 0; k < legendreVector[0].size(); k++) {
294  legendreVector[0][k] = OrthogonalPoly::EvalLegendre(k, u);
295  }
296  for(unsigned int k = 0; k < legendreVector[1].size(); k++) {
297  legendreVector[1][k] = OrthogonalPoly::EvalLegendre(k, v);
298  }
299  int iterator = 0;
300 
301  for(int it1 = 0; it1 <= _pf2; it1++) {
302  for(int it2 = 2; it2 <= _pf1 + 1; it2++) {
303  int impactFlag1 = 1;
304  int impactFlag2 = 1;
305  if(flag2 == -1 && it1 % 2 == 0) { impactFlag1 = -1; }
306  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
307  faceFunctions[iterator][0] = 0;
308  faceFunctions[iterator][1] = legendreVector[1][it1] *
309  OrthogonalPoly::EvalLobatto(it2, u) *
310  impactFlag1 * impactFlag2;
311  faceFunctions[iterator][2] = 0;
312  iterator++;
313  }
314  }
315  for(int it1 = 2; it1 <= _pf2 + 1; it1++) {
316  for(int it2 = 0; it2 <= _pf1; it2++) {
317  int impactFlag1 = 1;
318  int impactFlag2 = 1;
319  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
320  if(flag1 == -1 && it2 % 2 == 0) { impactFlag2 = -1; }
321  faceFunctions[iterator][0] = legendreVector[0][it2] *
322  OrthogonalPoly::EvalLobatto(it1, v) *
323  impactFlag1 * impactFlag2;
324  faceFunctions[iterator][1] = 0;
325  faceFunctions[iterator][2] = 0;
326  iterator++;
327  }
328  }
329  }
330  else if(typeFunction == "CurlHcurlLegendre") {
331  std::vector<std::vector<double> > legendreVector(2);
332  legendreVector[0] = std::vector<double>(_pf1 + 1);
333  legendreVector[1] = std::vector<double>(_pf2 + 1);
334  for(unsigned int k = 0; k < legendreVector[0].size(); k++) {
335  legendreVector[0][k] = OrthogonalPoly::EvalLegendre(k, u);
336  }
337  for(unsigned int k = 0; k < legendreVector[1].size(); k++) {
338  legendreVector[1][k] = OrthogonalPoly::EvalLegendre(k, v);
339  }
340  int iterator = 0;
341  for(int it1 = 0; it1 <= _pf2; it1++) {
342  for(int it2 = 2; it2 <= _pf1 + 1; it2++) {
343  int impactFlag1 = 1;
344  int impactFlag2 = 1;
345  if(flag2 == -1 && it1 % 2 == 0) { impactFlag1 = -1; }
346  if(flag1 == -1 && it2 % 2 != 0) { impactFlag2 = -1; }
347  faceFunctions[iterator][0] = 0;
348  faceFunctions[iterator][1] = 0;
349  faceFunctions[iterator][2] = legendreVector[1][it1] *
350  OrthogonalPoly::EvalDLobatto(it2, u) *
351  impactFlag1 * impactFlag2;
352  iterator++;
353  }
354  }
355  for(int it1 = 2; it1 <= _pf2 + 1; it1++) {
356  for(int it2 = 0; it2 <= _pf1; it2++) {
357  int impactFlag1 = 1;
358  int impactFlag2 = 1;
359  if(flag2 == -1 && it1 % 2 != 0) { impactFlag1 = -1; }
360  if(flag1 == -1 && it2 % 2 == 0) { impactFlag2 = -1; }
361  faceFunctions[iterator][0] = 0;
362  faceFunctions[iterator][1] = 0;
363  faceFunctions[iterator][2] = -legendreVector[0][it2] *
364  OrthogonalPoly::EvalDLobatto(it1, v) *
365  impactFlag1 * impactFlag2;
366  iterator++;
367  }
368  }
369  }
370  else {
371  throw std::runtime_error("unknown typeFunction");
372  }
373  }
374  }
375 }
376 void HierarchicalBasisHcurlQuad::orientFace(
377  int const &flag1, int const &flag2, int const &flag3, int const &faceNumber,
378  const std::vector<std::vector<double> > &quadFaceFunctionsAllOrientation,
379  const std::vector<std::vector<double> > &triFaceFunctionsAllOrientation,
380  std::vector<std::vector<double> > &fTableCopy)
381 {
382  int iOrientation = numberOrientationQuadFace(flag1, flag2, flag3);
383  int offset = iOrientation * _nQuadFaceFunction;
384  for(int i = 0; i < _nQuadFaceFunction; i++) {
385  fTableCopy[i][0] = quadFaceFunctionsAllOrientation[i + offset][0];
386  fTableCopy[i][1] = quadFaceFunctionsAllOrientation[i + offset][1];
387  fTableCopy[i][2] = quadFaceFunctionsAllOrientation[i + offset][2];
388  }
389 }
390 
391 void HierarchicalBasisHcurlQuad::getKeysInfo(std::vector<int> &functionTypeInfo,
392  std::vector<int> &orderInfo)
393 {
394  int it = 0;
395  for(int numEdge = 0; numEdge < 4; numEdge++) {
396  for(int i = 0; i <= _pOrderEdge[numEdge]; i++) {
397  functionTypeInfo[it] = 1;
398  orderInfo[it] = i;
399  it++;
400  }
401  }
402  for(int n1 = 0; n1 <= _pf1; n1++) {
403  for(int n2 = 2; n2 <= _pf2 + 1; n2++) {
404  functionTypeInfo[it] = 2;
405  orderInfo[it] = std::max(n1, n2);
406  it++;
407  }
408  }
409  for(int n1 = 2; n1 <= _pf1 + 1; n1++) {
410  for(int n2 = 0; n2 <= _pf2; n2++) {
411  functionTypeInfo[it] = 2;
412  orderInfo[it] = std::max(n1, n2);
413  it++;
414  }
415  }
416 }
OrthogonalPoly::EvalLegendre
double EvalLegendre(int order, double x)
Definition: OrthogonalPoly.cpp:415
HierarchicalBasis::_nTriFaceFunction
int _nTriFaceFunction
Definition: HierarchicalBasis.h:27
HierarchicalBasisHcurlQuad::HierarchicalBasisHcurlQuad
HierarchicalBasisHcurlQuad(int order)
Definition: HierarchicalBasisHcurlQuad.cpp:11
HierarchicalBasisHcurlQuad::orientEdgeFunctionsForNegativeFlag
virtual void orientEdgeFunctionsForNegativeFlag(std::vector< std::vector< double > > &edgeFunctions)
Definition: HierarchicalBasisHcurlQuad.cpp:157
HierarchicalBasisHcurlQuad::orientEdge
virtual void orientEdge(int const &flagOrientation, int const &edgeNumber, std::vector< std::vector< double > > &edgeBasis, const std::vector< std::vector< double > > &eTablePositiveFlag, const std::vector< std::vector< double > > &eTableNegativeFlag)
Definition: HierarchicalBasisHcurlQuad.cpp:126
HierarchicalBasis::_nvertex
int _nvertex
Definition: HierarchicalBasis.h:20
HierarchicalBasis::_nBubbleFunction
int _nBubbleFunction
Definition: HierarchicalBasis.h:28
HierarchicalBasisHcurlQuad::orientFace
virtual void orientFace(int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, const std::vector< std::vector< double > > &quadFaceFunctionsAllOrientation, const std::vector< std::vector< double > > &triFaceFunctionsAllOrientation, std::vector< std::vector< double > > &fTableCopy)
Definition: HierarchicalBasisHcurlQuad.cpp:376
OrthogonalPoly::EvalLobatto
double EvalLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:9
OrthogonalPoly::EvalDLobatto
double EvalDLobatto(int order, double x)
Definition: OrthogonalPoly.cpp:127
HierarchicalBasisHcurlQuad::_pOrderEdge
int _pOrderEdge[4]
Definition: HierarchicalBasisHcurlQuad.h:78
HierarchicalBasisHcurlQuad::~HierarchicalBasisHcurlQuad
virtual ~HierarchicalBasisHcurlQuad()
Definition: HierarchicalBasisHcurlQuad.cpp:31
HierarchicalBasisHcurlQuad::getKeysInfo
virtual void getKeysInfo(std::vector< int > &functionTypeInfo, std::vector< int > &orderInfo)
Definition: HierarchicalBasisHcurlQuad.cpp:391
HierarchicalBasisHcurlQuad.h
HierarchicalBasisHcurlQuad::generateCurlBasis
virtual void generateCurlBasis(double const &u, double const &v, double const &w, std::vector< std::vector< double > > &edgeBasis, std::vector< std::vector< double > > &faceBasis, std::vector< std::vector< double > > &bubbleBasis)
Definition: HierarchicalBasisHcurlQuad.cpp:178
HierarchicalBasisHcurlQuad::generateHcurlBasis
virtual void generateHcurlBasis(double const &u, double const &v, double const &w, std::vector< std::vector< double > > &edgeBasis, std::vector< std::vector< double > > &faceBasis, std::vector< std::vector< double > > &bubbleBasis)
Definition: HierarchicalBasisHcurlQuad.cpp:51
HierarchicalBasis::_nQuadFaceFunction
int _nQuadFaceFunction
Definition: HierarchicalBasis.h:26
HierarchicalBasis::_nfaceTri
int _nfaceTri
Definition: HierarchicalBasis.h:23
HierarchicalBasis::_nfaceQuad
int _nfaceQuad
Definition: HierarchicalBasis.h:22
direction
static void direction(Vertex *v1, Vertex *v2, double d[3])
Definition: Geo.cpp:218
HierarchicalBasis::_nVertexFunction
int _nVertexFunction
Definition: HierarchicalBasis.h:24
HierarchicalBasisHcurlQuad::getNumberOfOrientations
virtual unsigned int getNumberOfOrientations() const
Definition: HierarchicalBasisHcurlQuad.cpp:33
HierarchicalBasisHcurlQuad::orientOneFace
virtual void orientOneFace(double const &u, double const &v, double const &w, int const &flag1, int const &flag2, int const &flag3, int const &faceNumber, std::vector< std::vector< double > > &faceFunctions, std::string typeFunction)
Definition: HierarchicalBasisHcurlQuad.cpp:249
HierarchicalBasisHcurlQuad::_pf2
int _pf2
Definition: HierarchicalBasisHcurlQuad.h:77
HierarchicalBasisHcurlQuad::_pf1
int _pf1
Definition: HierarchicalBasisHcurlQuad.h:76
HierarchicalBasisHcurlQuad::_affineCoordinate
static double _affineCoordinate(int const &j, double const &u, double const &v)
Definition: HierarchicalBasisHcurlQuad.cpp:38
HierarchicalBasis::_nEdgeFunction
int _nEdgeFunction
Definition: HierarchicalBasis.h:25
HierarchicalBasis::_nedge
int _nedge
Definition: HierarchicalBasis.h:21
HierarchicalBasis::numberOrientationQuadFace
int numberOrientationQuadFace(int const &flag1, int const &flag2, int const &flag3)
Definition: HierarchicalBasis.h:105