gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
MPrism.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 #include "MPrism.h"
7 #include "Numeric.h"
8 #include "BasisFactory.h"
9 #include "Context.h"
10 #include "pointsGenerators.h"
11 #include "nodalBasis.h"
12 
13 #if defined(HAVE_MESH)
14 #include "qualityMeasures.h"
15 #endif
16 
17 std::map<int, IndicesReversed> MPrismN::_order2indicesReversedPri;
18 
19 void MPrism::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
20  SVector3 *n)
21 {
22  MVertex *v0 = _v[edges_prism(num, 0)];
23  MVertex *v1 = _v[edges_prism(num, 1)];
24  x[0] = v0->x();
25  y[0] = v0->y();
26  z[0] = v0->z();
27  x[1] = v1->x();
28  y[1] = v1->y();
29  z[1] = v1->z();
30  // just one of the potential normals - did not bother computing the normal of
31  // one of the faces - don't use MElement::_getEdgeRep as it uses MFace, which
32  // is slow
33  double nn[3];
34  normal2points(x[0], y[0], z[0], x[1], y[1], z[1], nn);
35  n[0] = n[1] = SVector3(nn[0], nn[1], nn[2]);
36 }
37 
38 int MPrism::getVolumeSign()
39 {
40  double mat[3][3];
41  mat[0][0] = _v[1]->x() - _v[0]->x();
42  mat[0][1] = _v[2]->x() - _v[0]->x();
43  mat[0][2] = _v[3]->x() - _v[0]->x();
44  mat[1][0] = _v[1]->y() - _v[0]->y();
45  mat[1][1] = _v[2]->y() - _v[0]->y();
46  mat[1][2] = _v[3]->y() - _v[0]->y();
47  mat[2][0] = _v[1]->z() - _v[0]->z();
48  mat[2][1] = _v[2]->z() - _v[0]->z();
49  mat[2][2] = _v[3]->z() - _v[0]->z();
50  double d = det3x3(mat);
51  if(d < 0.)
52  return -1;
53  else if(d > 0.)
54  return 1;
55  else
56  return 0;
57 }
58 
59 void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
60 {
61  *npts = getNGQPriPts(pOrder);
62  *pts = getGQPriPts(pOrder);
63 }
64 
65 double MPrism::getInnerRadius()
66 {
67  double dist[3], k = 0.;
68  int triEdges[3] = {0, 1, 3};
69  for(int i = 0; i < 3; i++) {
70  MEdge e = getEdge(triEdges[i]);
71  dist[i] = e.getVertex(0)->distance(e.getVertex(1));
72  k += 0.5 * dist[i];
73  }
74  double radTri = sqrt(k * (k - dist[0]) * (k - dist[1]) * (k - dist[2])) / k;
75  double radVert = 0.5 * getVertex(0)->distance(getVertex(3));
76  return std::min(radTri, radVert);
77 }
78 
79 bool MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign,
80  int &rot) const
81 {
82  for(ithFace = 0; ithFace < 5; ithFace++) {
83  if(_getFaceInfo(getFace(ithFace), face, sign, rot)) return true;
84  }
85  Msg::Error("Could not get face information for prism %d", getNum());
86  return false;
87 }
88 
89 int MPrism::numCommonNodesInDualGraph(const MElement *const other) const
90 {
91  switch(other->getType()) {
92  case TYPE_PNT: return 1;
93  case TYPE_LIN: return 2;
94  case TYPE_QUA: return 4;
95  case TYPE_HEX: return 4;
96  default: return 3;
97  }
98 }
99 
100 static void _myGetEdgeRep(MPrism *pri, int num, double *x, double *y, double *z,
101  SVector3 *n, int numSubEdges)
102 {
103  // const int numSubEdges = CTX::instance()->mesh.numSubEdges;
104  static double pp[6][3] = {{0, 0, -1}, {1, 0, -1}, {0, 1, -1},
105  {0, 0, 1}, {1, 0, 1}, {0, 1, 1}};
106  static int ed[9][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {1, 4},
107  {2, 5}, {3, 4}, {3, 5}, {4, 5}};
108  int iEdge = num / numSubEdges;
109  int iSubEdge = num % numSubEdges;
110 
111  int iVertex1 = ed[iEdge][0];
112  int iVertex2 = ed[iEdge][1];
113  double t1 = (double)iSubEdge / (double)numSubEdges;
114  double u1 = pp[iVertex1][0] * (1. - t1) + pp[iVertex2][0] * t1;
115  double v1 = pp[iVertex1][1] * (1. - t1) + pp[iVertex2][1] * t1;
116  double w1 = pp[iVertex1][2] * (1. - t1) + pp[iVertex2][2] * t1;
117 
118  double t2 = (double)(iSubEdge + 1) / (double)numSubEdges;
119  double u2 = pp[iVertex1][0] * (1. - t2) + pp[iVertex2][0] * t2;
120  double v2 = pp[iVertex1][1] * (1. - t2) + pp[iVertex2][1] * t2;
121  double w2 = pp[iVertex1][2] * (1. - t2) + pp[iVertex2][2] * t2;
122 
123  SPoint3 pnt1, pnt2;
124  pri->pnt(u1, v1, w1, pnt1);
125  pri->pnt(u2, v2, w2, pnt2);
126  x[0] = pnt1.x();
127  x[1] = pnt2.x();
128  y[0] = pnt1.y();
129  y[1] = pnt2.y();
130  z[0] = pnt1.z();
131  z[1] = pnt2.z();
132 
133  // not great, but better than nothing
134  // static const int f[6] = {0, 0, 0, 1, 2, 3};
135  n[0] = n[1] = 1;
136 }
137 
138 void MPrism15::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
139  SVector3 *n)
140 {
141  if(curved)
142  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
143  else
144  MPrism::getEdgeRep(false, num, x, y, z, n);
145 }
146 
147 void MPrism18::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
148  SVector3 *n)
149 {
150  if(curved)
151  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
152  else
153  MPrism::getEdgeRep(false, num, x, y, z, n);
154 }
155 
156 void MPrismN::getEdgeRep(bool curved, int num, double *x, double *y, double *z,
157  SVector3 *n)
158 {
159  if(curved)
160  _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
161  else
162  MPrism::getEdgeRep(false, num, x, y, z, n);
163 }
164 
165 int MPrism15::getNumEdgesRep(bool curved)
166 {
167  return curved ? 9 * CTX::instance()->mesh.numSubEdges : 9;
168 }
169 
170 int MPrism18::getNumEdgesRep(bool curved)
171 {
172  return curved ? 9 * CTX::instance()->mesh.numSubEdges : 9;
173 }
174 
175 int MPrismN::getNumEdgesRep(bool curved)
176 {
177  return curved ? 9 * CTX::instance()->mesh.numSubEdges : 9;
178 }
179 
180 static void _myGetFaceRep(MPrism *pri, int num, double *x, double *y, double *z,
181  SVector3 *n, int numSubEdges)
182 {
183  static double pp[6][3] = {{0, 0, -1}, {1, 0, -1}, {0, 1, -1},
184  {0, 0, 1}, {1, 0, 1}, {0, 1, 1}};
185 
186  int iFace = num / (numSubEdges * numSubEdges);
187  int iSubFace = num % (numSubEdges * numSubEdges);
188 
189  if(iFace > 1) {
190  iFace = num / (2 * numSubEdges * numSubEdges) + 1;
191  iSubFace = num % (2 * numSubEdges * numSubEdges);
192  }
193 
194  int iVertex1 = pri->faces_prism(iFace, 0);
195  int iVertex2 = pri->faces_prism(iFace, 1);
196  int iVertex3 = pri->faces_prism(iFace, 2);
197  int iVertex4 = pri->faces_prism(iFace, 3);
198 
199  SPoint3 pnt1, pnt2, pnt3;
200  // double J1[3][3], J2[3][3], J3[3][3];
201 
202  /*
203  0
204  0 1
205  0 1 2
206  0 1 2 3
207  0 1 2 3 4
208  0 1 2 3 4 5
209  */
210 
211  // on each layer, we have (numSubEdges) * 2 triangles
212  // ix and iy are the coordinates of the sub-quadrangle
213 
214  if(iFace > 1) {
215  int io = iSubFace % 2;
216  int ix = (iSubFace / 2) / numSubEdges;
217  int iy = (iSubFace / 2) % numSubEdges;
218 
219  const double d = 2. / numSubEdges;
220  double ox = -1. + d * ix;
221  double oy = -1. + d * iy;
222 
223  if(io == 0) {
224  double U1 = pp[iVertex1][0] * (1. - ox) * (1 - oy) * .25 +
225  pp[iVertex2][0] * (1. + ox) * (1 - oy) * .25 +
226  pp[iVertex3][0] * (1. + ox) * (1 + oy) * .25 +
227  pp[iVertex4][0] * (1. - ox) * (1 + oy) * .25;
228  double V1 = pp[iVertex1][1] * (1. - ox) * (1 - oy) * .25 +
229  pp[iVertex2][1] * (1. + ox) * (1 - oy) * .25 +
230  pp[iVertex3][1] * (1. + ox) * (1 + oy) * .25 +
231  pp[iVertex4][1] * (1. - ox) * (1 + oy) * .25;
232  double W1 = pp[iVertex1][2] * (1. - ox) * (1 - oy) * .25 +
233  pp[iVertex2][2] * (1. + ox) * (1 - oy) * .25 +
234  pp[iVertex3][2] * (1. + ox) * (1 + oy) * .25 +
235  pp[iVertex4][2] * (1. - ox) * (1 + oy) * .25;
236 
237  ox += d;
238 
239  double U2 = pp[iVertex1][0] * (1. - ox) * (1 - oy) * .25 +
240  pp[iVertex2][0] * (1. + ox) * (1 - oy) * .25 +
241  pp[iVertex3][0] * (1. + ox) * (1 + oy) * .25 +
242  pp[iVertex4][0] * (1. - ox) * (1 + oy) * .25;
243  double V2 = pp[iVertex1][1] * (1. - ox) * (1 - oy) * .25 +
244  pp[iVertex2][1] * (1. + ox) * (1 - oy) * .25 +
245  pp[iVertex3][1] * (1. + ox) * (1 + oy) * .25 +
246  pp[iVertex4][1] * (1. - ox) * (1 + oy) * .25;
247  double W2 = pp[iVertex1][2] * (1. - ox) * (1 - oy) * .25 +
248  pp[iVertex2][2] * (1. + ox) * (1 - oy) * .25 +
249  pp[iVertex3][2] * (1. + ox) * (1 + oy) * .25 +
250  pp[iVertex4][2] * (1. - ox) * (1 + oy) * .25;
251 
252  oy += d;
253 
254  double U3 = pp[iVertex1][0] * (1. - ox) * (1 - oy) * .25 +
255  pp[iVertex2][0] * (1. + ox) * (1 - oy) * .25 +
256  pp[iVertex3][0] * (1. + ox) * (1 + oy) * .25 +
257  pp[iVertex4][0] * (1. - ox) * (1 + oy) * .25;
258  double V3 = pp[iVertex1][1] * (1. - ox) * (1 - oy) * .25 +
259  pp[iVertex2][1] * (1. + ox) * (1 - oy) * .25 +
260  pp[iVertex3][1] * (1. + ox) * (1 + oy) * .25 +
261  pp[iVertex4][1] * (1. - ox) * (1 + oy) * .25;
262  double W3 = pp[iVertex1][2] * (1. - ox) * (1 - oy) * .25 +
263  pp[iVertex2][2] * (1. + ox) * (1 - oy) * .25 +
264  pp[iVertex3][2] * (1. + ox) * (1 + oy) * .25 +
265  pp[iVertex4][2] * (1. - ox) * (1 + oy) * .25;
266 
267  pri->pnt(U1, V1, W1, pnt1);
268  pri->pnt(U2, V2, W2, pnt2);
269  pri->pnt(U3, V3, W3, pnt3);
270  }
271  else {
272  double U1 = pp[iVertex1][0] * (1. - ox) * (1 - oy) * .25 +
273  pp[iVertex2][0] * (1. + ox) * (1 - oy) * .25 +
274  pp[iVertex3][0] * (1. + ox) * (1 + oy) * .25 +
275  pp[iVertex4][0] * (1. - ox) * (1 + oy) * .25;
276  double V1 = pp[iVertex1][1] * (1. - ox) * (1 - oy) * .25 +
277  pp[iVertex2][1] * (1. + ox) * (1 - oy) * .25 +
278  pp[iVertex3][1] * (1. + ox) * (1 + oy) * .25 +
279  pp[iVertex4][1] * (1. - ox) * (1 + oy) * .25;
280  double W1 = pp[iVertex1][2] * (1. - ox) * (1 - oy) * .25 +
281  pp[iVertex2][2] * (1. + ox) * (1 - oy) * .25 +
282  pp[iVertex3][2] * (1. + ox) * (1 + oy) * .25 +
283  pp[iVertex4][2] * (1. - ox) * (1 + oy) * .25;
284 
285  ox += d;
286  oy += d;
287 
288  double U2 = pp[iVertex1][0] * (1. - ox) * (1 - oy) * .25 +
289  pp[iVertex2][0] * (1. + ox) * (1 - oy) * .25 +
290  pp[iVertex3][0] * (1. + ox) * (1 + oy) * .25 +
291  pp[iVertex4][0] * (1. - ox) * (1 + oy) * .25;
292  double V2 = pp[iVertex1][1] * (1. - ox) * (1 - oy) * .25 +
293  pp[iVertex2][1] * (1. + ox) * (1 - oy) * .25 +
294  pp[iVertex3][1] * (1. + ox) * (1 + oy) * .25 +
295  pp[iVertex4][1] * (1. - ox) * (1 + oy) * .25;
296  double W2 = pp[iVertex1][2] * (1. - ox) * (1 - oy) * .25 +
297  pp[iVertex2][2] * (1. + ox) * (1 - oy) * .25 +
298  pp[iVertex3][2] * (1. + ox) * (1 + oy) * .25 +
299  pp[iVertex4][2] * (1. - ox) * (1 + oy) * .25;
300 
301  ox -= d;
302 
303  double U3 = pp[iVertex1][0] * (1. - ox) * (1 - oy) * .25 +
304  pp[iVertex2][0] * (1. + ox) * (1 - oy) * .25 +
305  pp[iVertex3][0] * (1. + ox) * (1 + oy) * .25 +
306  pp[iVertex4][0] * (1. - ox) * (1 + oy) * .25;
307  double V3 = pp[iVertex1][1] * (1. - ox) * (1 - oy) * .25 +
308  pp[iVertex2][1] * (1. + ox) * (1 - oy) * .25 +
309  pp[iVertex3][1] * (1. + ox) * (1 + oy) * .25 +
310  pp[iVertex4][1] * (1. - ox) * (1 + oy) * .25;
311  double W3 = pp[iVertex1][2] * (1. - ox) * (1 - oy) * .25 +
312  pp[iVertex2][2] * (1. + ox) * (1 - oy) * .25 +
313  pp[iVertex3][2] * (1. + ox) * (1 + oy) * .25 +
314  pp[iVertex4][2] * (1. - ox) * (1 + oy) * .25;
315 
316  pri->pnt(U1, V1, W1, pnt1);
317  pri->pnt(U2, V2, W2, pnt2);
318  pri->pnt(U3, V3, W3, pnt3);
319  }
320  }
321  else {
322  int ix = 0, iy = 0;
323  int nbt = 0;
324  for(int i = 0; i < numSubEdges; i++) {
325  int nbl = (numSubEdges - i - 1) * 2 + 1;
326  nbt += nbl;
327  if(nbt > iSubFace) {
328  iy = i;
329  ix = nbl - (nbt - iSubFace);
330  break;
331  }
332  }
333 
334  const double d = 1. / numSubEdges;
335 
336  double u1, v1, u2, v2, u3, v3;
337  if(ix % 2 == 0) {
338  u1 = ix / 2 * d;
339  v1 = iy * d;
340  u2 = (ix / 2 + 1) * d;
341  v2 = iy * d;
342  u3 = ix / 2 * d;
343  v3 = (iy + 1) * d;
344  }
345  else {
346  u1 = (ix / 2 + 1) * d;
347  v1 = iy * d;
348  u2 = (ix / 2 + 1) * d;
349  v2 = (iy + 1) * d;
350  u3 = ix / 2 * d;
351  v3 = (iy + 1) * d;
352  }
353 
354  double U1 = pp[iVertex1][0] * (1. - u1 - v1) + pp[iVertex2][0] * u1 +
355  pp[iVertex3][0] * v1;
356  double U2 = pp[iVertex1][0] * (1. - u2 - v2) + pp[iVertex2][0] * u2 +
357  pp[iVertex3][0] * v2;
358  double U3 = pp[iVertex1][0] * (1. - u3 - v3) + pp[iVertex2][0] * u3 +
359  pp[iVertex3][0] * v3;
360 
361  double V1 = pp[iVertex1][1] * (1. - u1 - v1) + pp[iVertex2][1] * u1 +
362  pp[iVertex3][1] * v1;
363  double V2 = pp[iVertex1][1] * (1. - u2 - v2) + pp[iVertex2][1] * u2 +
364  pp[iVertex3][1] * v2;
365  double V3 = pp[iVertex1][1] * (1. - u3 - v3) + pp[iVertex2][1] * u3 +
366  pp[iVertex3][1] * v3;
367 
368  double W1 = pp[iVertex1][2] * (1. - u1 - v1) + pp[iVertex2][2] * u1 +
369  pp[iVertex3][2] * v1;
370  double W2 = pp[iVertex1][2] * (1. - u2 - v2) + pp[iVertex2][2] * u2 +
371  pp[iVertex3][2] * v2;
372  double W3 = pp[iVertex1][2] * (1. - u3 - v3) + pp[iVertex2][2] * u3 +
373  pp[iVertex3][2] * v3;
374 
375  pri->pnt(U1, V1, W1, pnt1);
376  pri->pnt(U2, V2, W2, pnt2);
377  pri->pnt(U3, V3, W3, pnt3);
378  }
379 
380  x[0] = pnt1.x();
381  x[1] = pnt2.x();
382  x[2] = pnt3.x();
383  y[0] = pnt1.y();
384  y[1] = pnt2.y();
385  y[2] = pnt3.y();
386  z[0] = pnt1.z();
387  z[1] = pnt2.z();
388  z[2] = pnt3.z();
389 
390  SVector3 d1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
391  SVector3 d2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
392  n[0] = crossprod(d1, d2);
393  n[0].normalize();
394  n[1] = n[0];
395  n[2] = n[0];
396 }
397 
398 void MPrism::getFaceRep(bool curved, int num, double *x, double *y, double *z,
399  SVector3 *n)
400 {
401 #if defined(HAVE_VISUDEV)
402  static const int fquad[12][4] = {{0, 1, 4, 3}, {1, 4, 3, 0}, {4, 3, 0, 1},
403  {3, 0, 1, 4}, {1, 2, 5, 4}, {2, 5, 4, 1},
404  {5, 4, 1, 2}, {4, 1, 2, 5}, {2, 0, 3, 5},
405  {0, 3, 5, 2}, {3, 5, 2, 0}, {5, 2, 0, 3}};
406  if(CTX::instance()->heavyVisu) {
407  if(CTX::instance()->mesh.numSubEdges > 1) {
408  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
409  return;
410  }
411  if(num > 1) {
412  int i = num - 2;
413  _getFaceRepQuad(getVertex(fquad[i][0]), getVertex(fquad[i][1]),
414  getVertex(fquad[i][2]), getVertex(fquad[i][3]), x, y, z,
415  n);
416  return;
417  }
418  }
419 #endif
420  static const int f[8][3] = {{0, 2, 1}, {3, 4, 5}, {0, 1, 4}, {0, 4, 3},
421  {0, 3, 5}, {0, 5, 2}, {1, 2, 5}, {1, 5, 4}};
422  _getFaceRep(_v[f[num][0]], _v[f[num][1]], _v[f[num][2]], x, y, z, n);
423 }
424 
425 void MPrism15::getFaceRep(bool curved, int num, double *x, double *y, double *z,
426  SVector3 *n)
427 {
428  if(curved)
429  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
430  else
431  MPrism::getFaceRep(false, num, x, y, z, n);
432 }
433 
434 void MPrism18::getFaceRep(bool curved, int num, double *x, double *y, double *z,
435  SVector3 *n)
436 {
437  if(curved)
438  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
439  else
440  MPrism::getFaceRep(false, num, x, y, z, n);
441 }
442 
443 void MPrismN::getFaceRep(bool curved, int num, double *x, double *y, double *z,
444  SVector3 *n)
445 {
446  if(curved)
447  _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
448  else
449  MPrism::getFaceRep(false, num, x, y, z, n);
450 }
451 
452 int MPrism::getNumFacesRep(bool curved)
453 {
454 #if defined(HAVE_VISUDEV)
455  if(CTX::instance()->heavyVisu) {
456  if(CTX::instance()->mesh.numSubEdges == 1) return 14;
457  return 8 * std::pow(CTX::instance()->mesh.numSubEdges, 2);
458  }
459  if(CTX::instance()->heavyVisu) return 14;
460 #endif
461  return 8;
462 }
463 
464 int MPrism15::getNumFacesRep(bool curved)
465 {
466  return curved ? 8 * std::pow(CTX::instance()->mesh.numSubEdges, 2) :
467  MPrism::getNumFacesRep(curved);
468 }
469 
470 int MPrism18::getNumFacesRep(bool curved)
471 {
472  return curved ? 8 * std::pow(CTX::instance()->mesh.numSubEdges, 2) :
473  MPrism::getNumFacesRep(curved);
474 }
475 
476 int MPrismN::getNumFacesRep(bool curved)
477 {
478  return curved ? 8 * std::pow(CTX::instance()->mesh.numSubEdges, 2) :
479  MPrism::getNumFacesRep(curved);
480 }
481 
482 static void _addEdgeNodes(int num, bool reverse, int order,
483  const std::vector<MVertex *> &vs, int &ind,
484  std::vector<MVertex *> &v)
485 {
486  const int nNode = order - 1, startNode = num * nNode,
487  endNode = startNode + nNode - 1;
488 
489  if(reverse)
490  for(int i = endNode; i >= startNode; i--, ind++) v[ind] = vs[i];
491  else
492  for(int i = startNode; i <= endNode; i++, ind++) v[ind] = vs[i];
493 }
494 
495 static void _addFaceNodes(int num, int order, const std::vector<MVertex *> &vs,
496  int &ind, std::vector<MVertex *> &v)
497 {
498  const int nNodeEd = order - 1, nNodeTri = (order - 2) * (order - 1) / 2;
499 
500  int startNode, endNode;
501  if(num < 2) {
502  startNode = 9 * nNodeEd + num * nNodeTri;
503  endNode = startNode + nNodeTri;
504  }
505  else {
506  const int nNodeQuad = (order - 1) * (order - 1);
507  startNode = 9 * nNodeEd + 2 * nNodeTri + (num - 2) * nNodeQuad;
508  endNode = startNode + nNodeQuad;
509  }
510 
511  for(int i = startNode; i < endNode; i++, ind++) v[ind] = vs[i];
512 }
513 
514 // To be tested
515 void MPrismN::getFaceVertices(const int num, std::vector<MVertex *> &v) const
516 {
517  // FIXME serendipity case
518  static const int edge[5][4] = {
519  {1, 3, 0, -1}, {6, 8, 7, -1}, {0, 4, 6, 2}, {2, 7, 5, 1}, {3, 5, 8, 4}};
520  static const bool reverse[5][4] = {{false, true, true, false},
521  {false, false, true, false},
522  {false, false, true, true},
523  {false, false, true, true},
524  {false, false, true, true}};
525 
526  int nNodeTotal, nEdge;
527  if(num < 2) { // Triangular face
528  nNodeTotal = (_order + 1) * (_order + 2) / 2;
529  nEdge = 3;
530  }
531  else { // Quad face
532  nNodeTotal = (_order + 1) * (_order + 1);
533  nEdge = 4;
534  }
535 
536  v.resize(nNodeTotal);
537 
538  // Add corner nodes (there are nEdge corner nodes)
539  MPrism::_getFaceVertices(num, v);
540  int ind = nEdge;
541 
542  // Add edge nodes
543  for(int iE = 0; iE < nEdge; iE++)
544  _addEdgeNodes(edge[num][iE], reverse[num][iE], _order, _vs, ind, v);
545 
546  // Add face nodes
547  _addFaceNodes(num, _order, _vs, ind, v);
548 }
549 
550 double MPrism::gammaShapeMeasure()
551 {
552 #if defined(HAVE_MESH)
553  return qmPrism::minNCJ(this);
554 #else
555  return 0.;
556 #endif
557 }
558 
559 void _getIndicesReversedPri(int order, IndicesReversed &indices)
560 {
561  fullMatrix<double> ref = gmshGenerateMonomialsPrism(order);
562 
563  indices.resize(ref.size1());
564  for(int i = 0; i < ref.size1(); ++i) {
565  const double u = ref(i, 0);
566  const double v = ref(i, 1);
567  const double w = ref(i, 2);
568  for(int j = 0; j < ref.size1(); ++j) {
569  if(u == ref(j, 1) && v == ref(j, 0) && w == ref(j, 2)) {
570  indices[i] = j;
571  break;
572  }
573  }
574  }
575 }
576 
577 void MPrismN::reverse()
578 {
579  auto it = _order2indicesReversedPri.find(_order);
580  if(it == _order2indicesReversedPri.end()) {
581  IndicesReversed indices;
582  _getIndicesReversedPri(_order, indices);
583  _order2indicesReversedPri[_order] = indices;
584  it = _order2indicesReversedPri.find(_order);
585  }
586 
587  IndicesReversed &indices = it->second;
588 
589  // copy vertices
590  std::vector<MVertex *> oldv(6 + _vs.size());
591  std::copy(_v, _v + 6, oldv.begin());
592  std::copy(_vs.begin(), _vs.end(), oldv.begin() + 6);
593 
594  // reverse
595  for(int i = 0; i < 6; ++i) { _v[i] = oldv[indices[i]]; }
596  for(std::size_t i = 0; i < _vs.size(); ++i) { _vs[i] = oldv[indices[6 + i]]; }
597 }
598 
599 void MPrismN::getNode(int num, double &u, double &v, double &w) const
600 {
601  const fullMatrix<double> &p = getFunctionSpace()->points;
602  u = p(num, 0);
603  v = p(num, 1);
604  w = p(num, 2);
605 }
MElement::getNum
virtual std::size_t getNum() const
Definition: MElement.h:68
MPrism::_v
MVertex * _v[6]
Definition: MPrism.h:36
crossprod
SVector3 crossprod(const SVector3 &a, const SVector3 &b)
Definition: SVector3.h:150
MPrismN::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:443
qualityMeasures.h
MPrism::faces_prism
static int faces_prism(const int face, const int vert)
Definition: MPrism.h:188
MPrism::getEdge
virtual MEdge getEdge(int num) const
Definition: MPrism.h:75
MEdge
Definition: MEdge.h:14
MPrismN::getNumFacesRep
virtual int getNumFacesRep(bool curved)
Definition: MPrism.cpp:476
nodalBasis::points
fullMatrix< double > points
Definition: nodalBasis.h:16
_getIndicesReversedPri
void _getIndicesReversedPri(int order, IndicesReversed &indices)
Definition: MPrism.cpp:559
TYPE_LIN
#define TYPE_LIN
Definition: GmshDefines.h:65
MPrismN::reverse
virtual void reverse()
Definition: MPrism.cpp:577
MPrism::getInnerRadius
virtual double getInnerRadius()
Definition: MPrism.cpp:65
MVertex
Definition: MVertex.h:24
nodalBasis.h
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
MVertex::z
double z() const
Definition: MVertex.h:62
det3x3
double det3x3(double mat[3][3])
Definition: Numeric.cpp:126
SPoint3
Definition: SPoint3.h:14
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
TYPE_PNT
#define TYPE_PNT
Definition: GmshDefines.h:64
MPrism::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:19
SVector3
Definition: SVector3.h:16
MPrism::gammaShapeMeasure
virtual double gammaShapeMeasure()
Definition: MPrism.cpp:550
MPrism
Definition: MPrism.h:34
qmPrism::minNCJ
static double minNCJ(const MPrism *el)
Definition: qualityMeasures.cpp:769
w1
const double w1
Definition: GaussQuadratureHex.cpp:18
MPrismN::_order
const char _order
Definition: MPrism.h:541
CTX::instance
static CTX * instance()
Definition: Context.cpp:122
MPrism::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:398
SPoint3::x
double x(void) const
Definition: SPoint3.h:125
fullMatrix< double >
MFace
Definition: MFace.h:20
getNGQPriPts
int getNGQPriPts(int order, bool forceTensorRule=false)
Definition: GaussQuadraturePri.cpp:47
MEdge::getVertex
MVertex * getVertex(std::size_t i) const
Definition: MEdge.h:39
MPrism::_getFaceVertices
void _getFaceVertices(const int num, std::vector< MVertex * > &v) const
Definition: MPrism.h:42
contextMeshOptions::numSubEdges
int numSubEdges
Definition: Context.h:85
MPrism::getVertex
virtual MVertex * getVertex(int num)
Definition: MPrism.h:71
MElement::getType
virtual int getType() const =0
_myGetEdgeRep
static void _myGetEdgeRep(MPrism *pri, int num, double *x, double *y, double *z, SVector3 *n, int numSubEdges)
Definition: MPrism.cpp:100
MPrism::getFaceInfo
virtual bool getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) const
Definition: MPrism.cpp:79
Numeric.h
MPrism::edges_prism
static int edges_prism(const int edge, const int vert)
Definition: MPrism.h:182
MPrism18::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:147
MElement::_getFaceInfo
static bool _getFaceInfo(const MFace &face, const MFace &other, int &sign, int &rot)
Definition: MElement.cpp:66
SPoint3::y
double y(void) const
Definition: SPoint3.h:127
CTX::mesh
contextMeshOptions mesh
Definition: Context.h:313
MElement
Definition: MElement.h:30
MPrismN::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:156
MPrism::getVolumeSign
virtual int getVolumeSign()
Definition: MPrism.cpp:38
MPrismN::getFaceVertices
virtual void getFaceVertices(const int num, std::vector< MVertex * > &v) const
Definition: MPrism.cpp:515
MPrism::getIntegrationPoints
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
Definition: MPrism.cpp:59
MPrism15::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:425
MPrism18::getNumEdgesRep
virtual int getNumEdgesRep(bool curved)
Definition: MPrism.cpp:170
MPrismN::getNode
virtual void getNode(int num, double &u, double &v, double &w) const
Definition: MPrism.cpp:599
MPrism15::getEdgeRep
virtual void getEdgeRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:138
MElement::pnt
virtual void pnt(double u, double v, double w, SPoint3 &p) const
Definition: MElement.cpp:1072
_myGetFaceRep
static void _myGetFaceRep(MPrism *pri, int num, double *x, double *y, double *z, SVector3 *n, int numSubEdges)
Definition: MPrism.cpp:180
MElement::_getFaceRep
void _getFaceRep(MVertex *v0, MVertex *v1, MVertex *v2, double *x, double *y, double *z, SVector3 *n)
Definition: MElement.cpp:146
TYPE_QUA
#define TYPE_QUA
Definition: GmshDefines.h:67
MPrismN::getNumEdgesRep
virtual int getNumEdgesRep(bool curved)
Definition: MPrism.cpp:175
_addFaceNodes
static void _addFaceNodes(int num, int order, const std::vector< MVertex * > &vs, int &ind, std::vector< MVertex * > &v)
Definition: MPrism.cpp:495
_addEdgeNodes
static void _addEdgeNodes(int num, bool reverse, int order, const std::vector< MVertex * > &vs, int &ind, std::vector< MVertex * > &v)
Definition: MPrism.cpp:482
Context.h
IntPt
Definition: GaussIntegration.h:12
fullMatrix::size1
int size1() const
Definition: fullMatrix.h:274
z
const double z
Definition: GaussQuadratureQuad.cpp:56
MPrismN::_order2indicesReversedPri
static std::map< int, IndicesReversed > _order2indicesReversedPri
Definition: MPrism.h:537
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
MPrism15::getNumFacesRep
virtual int getNumFacesRep(bool curved)
Definition: MPrism.cpp:464
TYPE_HEX
#define TYPE_HEX
Definition: GmshDefines.h:71
MPrism.h
SPoint3::z
double z(void) const
Definition: SPoint3.h:129
MPrismN::_vs
std::vector< MVertex * > _vs
Definition: MPrism.h:540
MPrism::getNumFacesRep
virtual int getNumFacesRep(bool curved)
Definition: MPrism.cpp:452
MPrism18::getNumFacesRep
virtual int getNumFacesRep(bool curved)
Definition: MPrism.cpp:470
MVertex::y
double y() const
Definition: MVertex.h:61
normal2points
void normal2points(double x0, double y0, double z0, double x1, double y1, double z1, double n[3])
Definition: Numeric.cpp:85
MPrism18::getFaceRep
virtual void getFaceRep(bool curved, int num, double *x, double *y, double *z, SVector3 *n)
Definition: MPrism.cpp:434
pointsGenerators.h
MElement::getFunctionSpace
virtual const nodalBasis * getFunctionSpace(int order=-1, bool serendip=false) const
Definition: MElement.cpp:666
MPrism15::getNumEdgesRep
virtual int getNumEdgesRep(bool curved)
Definition: MPrism.cpp:165
MVertex::distance
double distance(MVertex *const v)
Definition: MVertex.h:105
gmshGenerateMonomialsPrism
fullMatrix< double > gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
Definition: pointsGenerators.cpp:402
getGQPriPts
IntPt * getGQPriPts(int order, bool forceTensorRule=false)
Definition: GaussQuadraturePri.cpp:12
MPrism::getFace
virtual MFace getFace(int num) const
Definition: MPrism.h:94
MVertex::x
double x() const
Definition: MVertex.h:60
SVector3::normalize
double normalize()
Definition: SVector3.h:38
IndicesReversed
std::vector< int > IndicesReversed
Definition: MHexahedron.h:605
BasisFactory.h
MPrism::numCommonNodesInDualGraph
virtual int numCommonNodesInDualGraph(const MElement *const other) const
Definition: MPrism.cpp:89