gmsh-TingyuanDoc  0.1
An Open-Source Timing-driven Analytical Mixed-size FPGA Placer
shapeFunctions.h
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 #ifndef SHAPE_FUNCTIONS_H
7 #define SHAPE_FUNCTIONS_H
8 
9 #include "Numeric.h"
10 #include "GmshMessage.h"
11 
12 class element {
13 protected:
14  bool _ownData;
15  double *_x, *_y, *_z;
16 
17 public:
18  element(double *x, double *y, double *z, int numNodes = 0)
19  {
20  if(!numNodes) {
21  _ownData = false;
22  _x = x;
23  _y = y;
24  _z = z;
25  }
26  else {
27  _ownData = true;
28  _x = new double[numNodes];
29  _y = new double[numNodes];
30  _z = new double[numNodes];
31  for(int i = 0; i < numNodes; i++) {
32  _x[i] = x[i];
33  _y[i] = y[i];
34  _z[i] = z[i];
35  }
36  }
37  }
38  virtual ~element()
39  {
40  if(_ownData) {
41  delete[] _x;
42  delete[] _y;
43  delete[] _z;
44  }
45  }
46  virtual void getXYZ(int num, double &x, double &y, double &z)
47  {
48  if(num < 0 || num >= getNumNodes()) return;
49  x = _x[num];
50  y = _y[num];
51  z = _z[num];
52  }
53  double getTolerance() const;
54  virtual int getDimension() = 0;
55  virtual int getNumNodes() = 0;
56  virtual void getNode(int num, double &u, double &v, double &w) = 0;
57  virtual int getNumEdges() = 0;
58  virtual void getEdge(int num, int &start, int &end) = 0;
59  virtual int getNumGaussPoints() = 0;
60  virtual void getGaussPoint(int num, double &u, double &v, double &w,
61  double &weight) = 0;
62  virtual void getShapeFunction(int num, double u, double v, double w,
63  double &s) = 0;
64  virtual void getGradShapeFunction(int num, double u, double v, double w,
65  double s[3]) = 0;
66  double getJacobian(double u, double v, double w, double jac[3][3])
67  {
68  jac[0][0] = jac[0][1] = jac[0][2] = 0.;
69  jac[1][0] = jac[1][1] = jac[1][2] = 0.;
70  jac[2][0] = jac[2][1] = jac[2][2] = 0.;
71  double s[3];
72  switch(getDimension()) {
73  case 3:
74  for(int i = 0; i < getNumNodes(); i++) {
75  getGradShapeFunction(i, u, v, w, s);
76  jac[0][0] += _x[i] * s[0];
77  jac[0][1] += _y[i] * s[0];
78  jac[0][2] += _z[i] * s[0];
79  jac[1][0] += _x[i] * s[1];
80  jac[1][1] += _y[i] * s[1];
81  jac[1][2] += _z[i] * s[1];
82  jac[2][0] += _x[i] * s[2];
83  jac[2][1] += _y[i] * s[2];
84  jac[2][2] += _z[i] * s[2];
85  }
86  return std::abs(
87  jac[0][0] * jac[1][1] * jac[2][2] + jac[0][2] * jac[1][0] * jac[2][1] +
88  jac[0][1] * jac[1][2] * jac[2][0] - jac[0][2] * jac[1][1] * jac[2][0] -
89  jac[0][0] * jac[1][2] * jac[2][1] - jac[0][1] * jac[1][0] * jac[2][2]);
90  case 2:
91  for(int i = 0; i < getNumNodes(); i++) {
92  getGradShapeFunction(i, u, v, w, s);
93  jac[0][0] += _x[i] * s[0];
94  jac[0][1] += _y[i] * s[0];
95  jac[0][2] += _z[i] * s[0];
96  jac[1][0] += _x[i] * s[1];
97  jac[1][1] += _y[i] * s[1];
98  jac[1][2] += _z[i] * s[1];
99  }
100  {
101  double a[3], b[3], c[3];
102  a[0] = _x[1] - _x[0];
103  a[1] = _y[1] - _y[0];
104  a[2] = _z[1] - _z[0];
105  b[0] = _x[2] - _x[0];
106  b[1] = _y[2] - _y[0];
107  b[2] = _z[2] - _z[0];
108  prodve(a, b, c);
109  jac[2][0] = c[0];
110  jac[2][1] = c[1];
111  jac[2][2] = c[2];
112  }
113  return std::sqrt(
114  std::pow(jac[0][0] * jac[1][1] - jac[0][1] * jac[1][0], 2) +
115  std::pow(jac[0][2] * jac[1][0] - jac[0][0] * jac[1][2], 2) +
116  std::pow(jac[0][1] * jac[1][2] - jac[0][2] * jac[1][1], 2));
117  case 1:
118  for(int i = 0; i < getNumNodes(); i++) {
119  getGradShapeFunction(i, u, v, w, s);
120  jac[0][0] += _x[i] * s[0];
121  jac[0][1] += _y[i] * s[0];
122  jac[0][2] += _z[i] * s[0];
123  }
124  {
125  double a[3], b[3], c[3];
126  a[0] = _x[1] - _x[0];
127  a[1] = _y[1] - _y[0];
128  a[2] = _z[1] - _z[0];
129  if((std::abs(a[0]) >= std::abs(a[1]) &&
130  std::abs(a[0]) >= std::abs(a[2])) ||
131  (std::abs(a[1]) >= std::abs(a[0]) &&
132  std::abs(a[1]) >= std::abs(a[2]))) {
133  b[0] = a[1];
134  b[1] = -a[0];
135  b[2] = 0.;
136  }
137  else {
138  b[0] = 0.;
139  b[1] = a[2];
140  b[2] = -a[1];
141  }
142  prodve(a, b, c);
143  jac[1][0] = b[0];
144  jac[1][1] = b[1];
145  jac[1][2] = b[2];
146  jac[2][0] = c[0];
147  jac[2][1] = c[1];
148  jac[2][2] = c[2];
149  }
150  return std::sqrt(std::pow(jac[0][0], 2) + std::pow(jac[0][1], 2) +
151  std::pow(jac[0][2], 2));
152  default: jac[0][0] = jac[1][1] = jac[2][2] = 1.; return 1.;
153  }
154  }
155  double interpolate(double val[], double u, double v, double w, int stride = 1)
156  {
157  double sum = 0;
158  int j = 0;
159  for(int i = 0; i < getNumNodes(); i++) {
160  double s;
161  getShapeFunction(i, u, v, w, s);
162  sum += val[j] * s;
163  j += stride;
164  }
165  return sum;
166  }
167  void interpolateGrad(double val[], double u, double v, double w, double f[3],
168  int stride = 1, double invjac[3][3] = nullptr)
169  {
170  double dfdu[3] = {0., 0., 0.};
171  int j = 0;
172  for(int i = 0; i < getNumNodes(); i++) {
173  double s[3];
174  getGradShapeFunction(i, u, v, w, s);
175  dfdu[0] += val[j] * s[0];
176  dfdu[1] += val[j] * s[1];
177  dfdu[2] += val[j] * s[2];
178  j += stride;
179  }
180  if(invjac) { matvec(invjac, dfdu, f); }
181  else {
182  double jac[3][3], inv[3][3];
183  getJacobian(u, v, w, jac);
184  inv3x3(jac, inv);
185  matvec(inv, dfdu, f);
186  }
187  }
188  void interpolateCurl(double val[], double u, double v, double w, double f[3],
189  int stride = 3)
190  {
191  double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
192  getJacobian(u, v, w, jac);
193  inv3x3(jac, inv);
194  interpolateGrad(&val[0], u, v, w, fx, stride, inv);
195  interpolateGrad(&val[1], u, v, w, fy, stride, inv);
196  interpolateGrad(&val[2], u, v, w, fz, stride, inv);
197  f[0] = fz[1] - fy[2];
198  f[1] = -(fz[0] - fx[2]);
199  f[2] = fy[0] - fx[1];
200  }
201  double interpolateDiv(double val[], double u, double v, double w,
202  int stride = 3)
203  {
204  double fx[3], fy[3], fz[3], jac[3][3], inv[3][3];
205  getJacobian(u, v, w, jac);
206  inv3x3(jac, inv);
207  interpolateGrad(&val[0], u, v, w, fx, stride, inv);
208  interpolateGrad(&val[1], u, v, w, fy, stride, inv);
209  interpolateGrad(&val[2], u, v, w, fz, stride, inv);
210  return fx[0] + fy[1] + fz[2];
211  }
212  double integrate(double val[], int stride = 1)
213  {
214  double sum = 0;
215  for(int i = 0; i < getNumGaussPoints(); i++) {
216  double u, v, w, weight, jac[3][3];
217  getGaussPoint(i, u, v, w, weight);
218  double det = getJacobian(u, v, w, jac);
219  double d = interpolate(val, u, v, w, stride);
220  sum += d * weight * det;
221  }
222  return sum;
223  }
224  double integrateLevelsetPositive(double val[])
225  {
226  // FIXME: explain + generalize this
227  double ones[8] = {1., 1., 1., 1., 1., 1., 1., 1.};
228  double area = integrate(ones);
229  double sum = 0, sumabs = 0.;
230  for(int i = 0; i < getNumNodes(); i++) {
231  sum += val[i];
232  sumabs += std::abs(val[i]);
233  }
234  double res = 0.;
235  if(sumabs) res = area * (1 + sum / sumabs) * 0.5;
236  return res;
237  }
238  virtual double integrateCirculation(double val[])
239  {
240  Msg::Error("integrateCirculation not available for this element");
241  return 0.;
242  }
243  virtual double integrateFlux(double val[])
244  {
245  Msg::Error("integrateFlux not available for this element");
246  return 0.;
247  }
248  virtual void xyz2uvw(double xyz[3], double uvw[3])
249  {
250  // general newton routine for the nonlinear case (more efficient
251  // routines are implemented for simplices, where the basis
252  // functions are linear)
253  uvw[0] = uvw[1] = uvw[2] = 0.0;
254 
255  int iter = 1, maxiter = 20;
256  double error = 1., tol = 1.e-6;
257 
258  while(error > tol && iter < maxiter) {
259  double jac[3][3];
260  if(!getJacobian(uvw[0], uvw[1], uvw[2], jac)) break;
261 
262  double xn = 0., yn = 0., zn = 0.;
263  for(int i = 0; i < getNumNodes(); i++) {
264  double s;
265  getShapeFunction(i, uvw[0], uvw[1], uvw[2], s);
266  xn += _x[i] * s;
267  yn += _y[i] * s;
268  zn += _z[i] * s;
269  }
270  double inv[3][3];
271  inv3x3(jac, inv);
272 
273  double un = uvw[0] + inv[0][0] * (xyz[0] - xn) +
274  inv[1][0] * (xyz[1] - yn) + inv[2][0] * (xyz[2] - zn);
275  double vn = uvw[1] + inv[0][1] * (xyz[0] - xn) +
276  inv[1][1] * (xyz[1] - yn) + inv[2][1] * (xyz[2] - zn);
277  double wn = uvw[2] + inv[0][2] * (xyz[0] - xn) +
278  inv[1][2] * (xyz[1] - yn) + inv[2][2] * (xyz[2] - zn);
279 
280  error = std::sqrt(std::pow(un - uvw[0], 2) + std::pow(vn - uvw[1], 2) +
281  std::pow(wn - uvw[2], 2));
282  uvw[0] = un;
283  uvw[1] = vn;
284  uvw[2] = wn;
285  iter++;
286  }
287  // if(error > tol) Msg::Warning("Newton did not converge in xyz2uvw") ;
288  }
289  virtual int isInside(double u, double v, double w) = 0;
290  double maxEdgeLength()
291  {
292  double max = 0.;
293  for(int i = 0; i < getNumEdges(); i++) {
294  int n1, n2;
295  getEdge(i, n1, n2);
296  double d =
297  std::sqrt(std::pow(_x[n1] - _x[n2], 2) + std::pow(_y[n1] - _y[n2], 2) +
298  std::pow(_z[n1] - _z[n2], 2));
299  if(d > max) max = d;
300  }
301  return max;
302  }
303 };
304 
305 class point : public element {
306 public:
307  point(double *x, double *y, double *z, int numNodes = 0)
308  : element(x, y, z, numNodes)
309  {
310  }
311  inline int getDimension() { return 0; }
312  inline int getNumNodes() { return 1; }
313  void getNode(int num, double &u, double &v, double &w) { u = v = w = 0.; }
314  inline int getNumEdges() { return 0; }
315  void getEdge(int num, int &start, int &end) { start = end = 0; }
316  inline int getNumGaussPoints() { return 1; }
317  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
318  {
319  u = v = w = 0.;
320  weight = 1.;
321  }
322  void getShapeFunction(int num, double u, double v, double w, double &s)
323  {
324  switch(num) {
325  case 0: s = 1.; break;
326  default: s = 0.; break;
327  }
328  }
329  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
330  {
331  s[0] = s[1] = s[2] = 0.;
332  }
333  void xyz2uvw(double xyz[3], double uvw[3]) { uvw[0] = uvw[1] = uvw[2] = 0.; }
334  int isInside(double u, double v, double w)
335  {
336  double TOL = getTolerance();
337  if(std::abs(u) > TOL || std::abs(v) > TOL || std::abs(w) > TOL) return 0;
338  return 1;
339  }
340 };
341 
342 class line : public element {
343 public:
344  line(double *x, double *y, double *z, int numNodes = 0)
345  : element(x, y, z, numNodes)
346  {
347  }
348  inline int getDimension() { return 1; }
349  inline int getNumNodes() { return 2; }
350  void getNode(int num, double &u, double &v, double &w)
351  {
352  v = w = 0.;
353  switch(num) {
354  case 0: u = -1.; break;
355  case 1: u = 1.; break;
356  default: u = 0.; break;
357  }
358  }
359  inline int getNumEdges() { return 1; }
360  void getEdge(int num, int &start, int &end)
361  {
362  start = 0;
363  end = 1;
364  }
365  inline int getNumGaussPoints() { return 1; }
366  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
367  {
368  if(num < 0 || num > 0) return;
369  u = v = w = 0.;
370  weight = 2.;
371  }
372  void getShapeFunction(int num, double u, double v, double w, double &s)
373  {
374  switch(num) {
375  case 0: s = 0.5 * (1. - u); break;
376  case 1: s = 0.5 * (1. + u); break;
377  default: s = 0.; break;
378  }
379  }
380  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
381  {
382  switch(num) {
383  case 0:
384  s[0] = -0.5;
385  s[1] = 0.;
386  s[2] = 0.;
387  break;
388  case 1:
389  s[0] = 0.5;
390  s[1] = 0.;
391  s[2] = 0.;
392  break;
393  default: s[0] = s[1] = s[2] = 0.; break;
394  }
395  }
396  double integrateCirculation(double val[])
397  {
398  double t[3] = {_x[1] - _x[0], _y[1] - _y[0], _z[1] - _z[0]};
399  norme(t);
400  double v[3];
401  for(int i = 0; i < 3; i++) v[i] = integrate(&val[i], 3);
402  return prosca(t, v);
403  }
404  int isInside(double u, double v, double w)
405  {
406  double TOL = getTolerance();
407  if(u < -(1. + TOL) || u > (1. + TOL) || std::abs(v) > TOL ||
408  std::abs(w) > TOL)
409  return 0;
410  return 1;
411  }
412 };
413 
414 class triangle : public element {
415 public:
416  triangle(double *x, double *y, double *z, int numNodes = 0)
417  : element(x, y, z, numNodes)
418  {
419  }
420  inline int getDimension() { return 2; }
421  inline int getNumNodes() { return 3; }
422  void getNode(int num, double &u, double &v, double &w)
423  {
424  w = 0.;
425  switch(num) {
426  case 0:
427  u = 0.;
428  v = 0.;
429  break;
430  case 1:
431  u = 1.;
432  v = 0.;
433  break;
434  case 2:
435  u = 0.;
436  v = 1.;
437  break;
438  default:
439  u = 0.;
440  v = 0.;
441  break;
442  }
443  }
444  inline int getNumEdges() { return 3; }
445  void getEdge(int num, int &start, int &end)
446  {
447  switch(num) {
448  case 0:
449  start = 0;
450  end = 1;
451  break;
452  case 1:
453  start = 1;
454  end = 2;
455  break;
456  case 2:
457  start = 2;
458  end = 0;
459  break;
460  default: start = end = 0; break;
461  }
462  }
463  inline int getNumGaussPoints() { return /* 3 */ 1; }
464  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
465  {
466  /*
467  static double u3[3] =
468  {0.16666666666666,0.66666666666666,0.16666666666666}; static double v3[3]
469  = {0.16666666666666,0.16666666666666,0.66666666666666}; static double
470  p3[3] = {0.16666666666666,0.16666666666666,0.16666666666666}; if(num < 0
471  || num > 2) return; u = u3[num]; v = v3[num]; w = 0.; weight = p3[num];
472  */
473  if(num < 0 || num > 0) return;
474  u = v = 0.333333333333333;
475  w = 0.;
476  weight = 0.5;
477  }
478  void getShapeFunction(int num, double u, double v, double w, double &s)
479  {
480  switch(num) {
481  case 0: s = 1. - u - v; break;
482  case 1: s = u; break;
483  case 2: s = v; break;
484  default: s = 0.; break;
485  }
486  }
487  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
488  {
489  switch(num) {
490  case 0:
491  s[0] = -1.;
492  s[1] = -1.;
493  s[2] = 0.;
494  break;
495  case 1:
496  s[0] = 1.;
497  s[1] = 0.;
498  s[2] = 0.;
499  break;
500  case 2:
501  s[0] = 0.;
502  s[1] = 1.;
503  s[2] = 0.;
504  break;
505  default: s[0] = s[1] = s[2] = 0.; break;
506  }
507  }
508  double integrateFlux(double val[])
509  {
510  double t1[3] = {_x[1] - _x[0], _y[1] - _y[0], _z[1] - _z[0]};
511  double t2[3] = {_x[2] - _x[0], _y[2] - _y[0], _z[2] - _z[0]};
512  double n[3];
513  prodve(t1, t2, n);
514  norme(n);
515  double v[3];
516  for(int i = 0; i < 3; i++) v[i] = integrate(&val[i], 3);
517  return prosca(n, v);
518  }
519  void xyz2uvw(double xyz[3], double uvw[3])
520  {
521  const double O[3] = {_x[0], _y[0], _z[0]};
522 
523  const double d[3] = {xyz[0] - O[0], xyz[1] - O[1], xyz[2] - O[2]};
524  const double d1[3] = {_x[1] - O[0], _y[1] - O[1], _z[1] - O[2]};
525  const double d2[3] = {_x[2] - O[0], _y[2] - O[1], _z[2] - O[2]};
526 
527  const double Jxy = d1[0] * d2[1] - d1[1] * d2[0];
528  const double Jxz = d1[0] * d2[2] - d1[2] * d2[0];
529  const double Jyz = d1[1] * d2[2] - d1[2] * d2[1];
530 
531  if((std::abs(Jxy) > std::abs(Jxz)) && (std::abs(Jxy) > std::abs(Jyz))) {
532  uvw[0] = (d[0] * d2[1] - d[1] * d2[0]) / Jxy;
533  uvw[1] = (d[1] * d1[0] - d[0] * d1[1]) / Jxy;
534  }
535  else if(std::abs(Jxz) > std::abs(Jyz)) {
536  uvw[0] = (d[0] * d2[2] - d[2] * d2[0]) / Jxz;
537  uvw[1] = (d[2] * d1[0] - d[0] * d1[2]) / Jxz;
538  }
539  else {
540  uvw[0] = (d[1] * d2[2] - d[2] * d2[1]) / Jyz;
541  uvw[1] = (d[2] * d1[1] - d[1] * d1[2]) / Jyz;
542  }
543  uvw[2] = 0.0;
544  }
545  int isInside(double u, double v, double w)
546  {
547  double TOL = getTolerance();
548  if(u < -TOL || v < -TOL || u > ((1. + TOL) - v) || std::abs(w) > TOL)
549  return 0;
550  return 1;
551  }
552 };
553 
554 class quadrangle : public element {
555 public:
556  quadrangle(double *x, double *y, double *z, int numNodes = 0)
557  : element(x, y, z, numNodes)
558  {
559  }
560  inline int getDimension() { return 2; }
561  inline int getNumNodes() { return 4; }
562  void getNode(int num, double &u, double &v, double &w)
563  {
564  w = 0.;
565  switch(num) {
566  case 0:
567  u = -1.;
568  v = -1.;
569  break;
570  case 1:
571  u = 1.;
572  v = -1.;
573  break;
574  case 2:
575  u = 1.;
576  v = 1.;
577  break;
578  case 3:
579  u = -1.;
580  v = 1.;
581  break;
582  default:
583  u = 0.;
584  v = 0.;
585  break;
586  }
587  }
588  inline int getNumEdges() { return 4; }
589  void getEdge(int num, int &start, int &end)
590  {
591  switch(num) {
592  case 0:
593  start = 0;
594  end = 1;
595  break;
596  case 1:
597  start = 1;
598  end = 2;
599  break;
600  case 2:
601  start = 2;
602  end = 3;
603  break;
604  case 3:
605  start = 3;
606  end = 0;
607  break;
608  default: start = end = 0; break;
609  }
610  }
611  inline int getNumGaussPoints() { return 4; }
612  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
613  {
614  static double u4[4] = {0.577350269189, -0.577350269189, 0.577350269189,
615  -0.577350269189};
616  static double v4[4] = {0.577350269189, 0.577350269189, -0.577350269189,
617  -0.577350269189};
618  static double p4[4] = {1., 1., 1., 1.};
619  if(num < 0 || num > 3) return;
620  u = u4[num];
621  v = v4[num];
622  w = 0.;
623  weight = p4[num];
624  }
625  void getShapeFunction(int num, double u, double v, double w, double &s)
626  {
627  switch(num) {
628  case 0: s = 0.25 * (1. - u) * (1. - v); break;
629  case 1: s = 0.25 * (1. + u) * (1. - v); break;
630  case 2: s = 0.25 * (1. + u) * (1. + v); break;
631  case 3: s = 0.25 * (1. - u) * (1. + v); break;
632  default: s = 0.; break;
633  }
634  }
635  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
636  {
637  switch(num) {
638  case 0:
639  s[0] = -0.25 * (1. - v);
640  s[1] = -0.25 * (1. - u);
641  s[2] = 0.;
642  break;
643  case 1:
644  s[0] = 0.25 * (1. - v);
645  s[1] = -0.25 * (1. + u);
646  s[2] = 0.;
647  break;
648  case 2:
649  s[0] = 0.25 * (1. + v);
650  s[1] = 0.25 * (1. + u);
651  s[2] = 0.;
652  break;
653  case 3:
654  s[0] = -0.25 * (1. + v);
655  s[1] = 0.25 * (1. - u);
656  s[2] = 0.;
657  break;
658  default: s[0] = s[1] = s[2] = 0.; break;
659  }
660  }
661  double integrateFlux(double val[])
662  {
663  double t1[3] = {_x[1] - _x[0], _y[1] - _y[0], _z[1] - _z[0]};
664  double t2[3] = {_x[2] - _x[0], _y[2] - _y[0], _z[2] - _z[0]};
665  double n[3];
666  prodve(t1, t2, n);
667  norme(n);
668  double v[3];
669  for(int i = 0; i < 3; i++) v[i] = integrate(&val[i], 3);
670  return prosca(n, v);
671  }
672  int isInside(double u, double v, double w)
673  {
674  double TOL = getTolerance();
675  if(u < -(1. + TOL) || v < -(1. + TOL) || u > (1. + TOL) || v > (1. + TOL) ||
676  std::abs(w) > TOL)
677  return 0;
678  return 1;
679  }
680 };
681 
682 class tetrahedron : public element {
683 public:
684  tetrahedron(double *x, double *y, double *z, int numNodes = 0)
685  : element(x, y, z, numNodes)
686  {
687  }
688  inline int getDimension() { return 3; }
689  inline int getNumNodes() { return 4; }
690  void getNode(int num, double &u, double &v, double &w)
691  {
692  switch(num) {
693  case 0:
694  u = 0.;
695  v = 0.;
696  w = 0.;
697  break;
698  case 1:
699  u = 1.;
700  v = 0.;
701  w = 0.;
702  break;
703  case 2:
704  u = 0.;
705  v = 1.;
706  w = 0.;
707  break;
708  case 3:
709  u = 0.;
710  v = 0.;
711  w = 1.;
712  break;
713  default:
714  u = 0.;
715  v = 0.;
716  w = 0.;
717  break;
718  }
719  }
720  inline int getNumEdges() { return 6; }
721  void getEdge(int num, int &start, int &end)
722  {
723  switch(num) {
724  case 0:
725  start = 0;
726  end = 1;
727  break;
728  case 1:
729  start = 1;
730  end = 2;
731  break;
732  case 2:
733  start = 2;
734  end = 0;
735  break;
736  case 3:
737  start = 3;
738  end = 0;
739  break;
740  case 4:
741  start = 3;
742  end = 2;
743  break;
744  case 5:
745  start = 3;
746  end = 1;
747  break;
748  default: start = end = 0; break;
749  }
750  }
751  inline int getNumGaussPoints() { return 4; }
752  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
753  {
754  static double u4[4] = {0.138196601125, 0.138196601125, 0.138196601125,
755  0.585410196625};
756  static double v4[4] = {0.138196601125, 0.138196601125, 0.585410196625,
757  0.138196601125};
758  static double w4[4] = {0.138196601125, 0.585410196625, 0.138196601125,
759  0.138196601125};
760  static double p4[4] = {0.0416666666667, 0.0416666666667, 0.0416666666667,
761  0.0416666666667};
762  if(num < 0 || num > 3) return;
763  u = u4[num];
764  v = v4[num];
765  w = w4[num];
766  weight = p4[num];
767  }
768  void getShapeFunction(int num, double u, double v, double w, double &s)
769  {
770  switch(num) {
771  case 0: s = 1. - u - v - w; break;
772  case 1: s = u; break;
773  case 2: s = v; break;
774  case 3: s = w; break;
775  default: s = 0.; break;
776  }
777  }
778  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
779  {
780  switch(num) {
781  case 0:
782  s[0] = -1.;
783  s[1] = -1.;
784  s[2] = -1.;
785  break;
786  case 1:
787  s[0] = 1.;
788  s[1] = 0.;
789  s[2] = 0.;
790  break;
791  case 2:
792  s[0] = 0.;
793  s[1] = 1.;
794  s[2] = 0.;
795  break;
796  case 3:
797  s[0] = 0.;
798  s[1] = 0.;
799  s[2] = 1.;
800  break;
801  default: s[0] = s[1] = s[2] = 0.; break;
802  }
803  }
804  void xyz2uvw(double xyz[3], double uvw[3])
805  {
806  double mat[3][3], b[3];
807  mat[0][0] = _x[1] - _x[0];
808  mat[0][1] = _x[2] - _x[0];
809  mat[0][2] = _x[3] - _x[0];
810  mat[1][0] = _y[1] - _y[0];
811  mat[1][1] = _y[2] - _y[0];
812  mat[1][2] = _y[3] - _y[0];
813  mat[2][0] = _z[1] - _z[0];
814  mat[2][1] = _z[2] - _z[0];
815  mat[2][2] = _z[3] - _z[0];
816  b[0] = xyz[0] - _x[0];
817  b[1] = xyz[1] - _y[0];
818  b[2] = xyz[2] - _z[0];
819  double det;
820  sys3x3(mat, b, uvw, &det);
821  }
822  int isInside(double u, double v, double w)
823  {
824  double TOL = getTolerance();
825  if(u < (-TOL) || v < (-TOL) || w < (-TOL) || u > ((1. + TOL) - v - w))
826  return 0;
827  return 1;
828  }
829 };
830 
831 class hexahedron : public element {
832 public:
833  hexahedron(double *x, double *y, double *z, int numNodes = 0)
834  : element(x, y, z, numNodes)
835  {
836  }
837  inline int getDimension() { return 3; }
838  inline int getNumNodes() { return 8; }
839  void getNode(int num, double &u, double &v, double &w)
840  {
841  switch(num) {
842  case 0:
843  u = -1.;
844  v = -1.;
845  w = -1.;
846  break;
847  case 1:
848  u = 1.;
849  v = -1.;
850  w = -1.;
851  break;
852  case 2:
853  u = 1.;
854  v = 1.;
855  w = -1.;
856  break;
857  case 3:
858  u = -1.;
859  v = 1.;
860  w = -1.;
861  break;
862  case 4:
863  u = -1.;
864  v = -1.;
865  w = 1.;
866  break;
867  case 5:
868  u = 1.;
869  v = -1.;
870  w = 1.;
871  break;
872  case 6:
873  u = 1.;
874  v = 1.;
875  w = 1.;
876  break;
877  case 7:
878  u = -1.;
879  v = 1.;
880  w = 1.;
881  break;
882  default:
883  u = 0.;
884  v = 0.;
885  w = 0.;
886  break;
887  }
888  }
889  inline int getNumEdges() { return 12; }
890  void getEdge(int num, int &start, int &end)
891  {
892  switch(num) {
893  case 0:
894  start = 0;
895  end = 1;
896  break;
897  case 1:
898  start = 0;
899  end = 3;
900  break;
901  case 2:
902  start = 0;
903  end = 4;
904  break;
905  case 3:
906  start = 1;
907  end = 2;
908  break;
909  case 4:
910  start = 1;
911  end = 5;
912  break;
913  case 5:
914  start = 2;
915  end = 3;
916  break;
917  case 6:
918  start = 2;
919  end = 6;
920  break;
921  case 7:
922  start = 3;
923  end = 7;
924  break;
925  case 8:
926  start = 4;
927  end = 5;
928  break;
929  case 9:
930  start = 4;
931  end = 7;
932  break;
933  case 10:
934  start = 5;
935  end = 6;
936  break;
937  case 11:
938  start = 6;
939  end = 7;
940  break;
941  default: start = end = 0; break;
942  }
943  }
944  inline int getNumGaussPoints() { return 6; }
945  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
946  {
947  static double u6[6] = {0.40824826, 0.40824826, -0.40824826,
948  -0.40824826, -0.81649658, 0.81649658};
949  static double v6[6] = {0.70710678, -0.70710678, 0.70710678,
950  -0.70710678, 0., 0.};
951  static double w6[6] = {-0.57735027, -0.57735027, 0.57735027,
952  0.57735027, -0.57735027, 0.57735027};
953  static double p6[6] = {1.3333333333, 1.3333333333, 1.3333333333,
954  1.3333333333, 1.3333333333, 1.3333333333};
955  if(num < 0 || num > 5) return;
956  u = u6[num];
957  v = v6[num];
958  w = w6[num];
959  weight = p6[num];
960  }
961  void getShapeFunction(int num, double u, double v, double w, double &s)
962  {
963  switch(num) {
964  case 0: s = (1. - u) * (1. - v) * (1. - w) * 0.125; break;
965  case 1: s = (1. + u) * (1. - v) * (1. - w) * 0.125; break;
966  case 2: s = (1. + u) * (1. + v) * (1. - w) * 0.125; break;
967  case 3: s = (1. - u) * (1. + v) * (1. - w) * 0.125; break;
968  case 4: s = (1. - u) * (1. - v) * (1. + w) * 0.125; break;
969  case 5: s = (1. + u) * (1. - v) * (1. + w) * 0.125; break;
970  case 6: s = (1. + u) * (1. + v) * (1. + w) * 0.125; break;
971  case 7: s = (1. - u) * (1. + v) * (1. + w) * 0.125; break;
972  default: s = 0.; break;
973  }
974  }
975  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
976  {
977  switch(num) {
978  case 0:
979  s[0] = -0.125 * (1. - v) * (1. - w);
980  s[1] = -0.125 * (1. - u) * (1. - w);
981  s[2] = -0.125 * (1. - u) * (1. - v);
982  break;
983  case 1:
984  s[0] = 0.125 * (1. - v) * (1. - w);
985  s[1] = -0.125 * (1. + u) * (1. - w);
986  s[2] = -0.125 * (1. + u) * (1. - v);
987  break;
988  case 2:
989  s[0] = 0.125 * (1. + v) * (1. - w);
990  s[1] = 0.125 * (1. + u) * (1. - w);
991  s[2] = -0.125 * (1. + u) * (1. + v);
992  break;
993  case 3:
994  s[0] = -0.125 * (1. + v) * (1. - w);
995  s[1] = 0.125 * (1. - u) * (1. - w);
996  s[2] = -0.125 * (1. - u) * (1. + v);
997  break;
998  case 4:
999  s[0] = -0.125 * (1. - v) * (1. + w);
1000  s[1] = -0.125 * (1. - u) * (1. + w);
1001  s[2] = 0.125 * (1. - u) * (1. - v);
1002  break;
1003  case 5:
1004  s[0] = 0.125 * (1. - v) * (1. + w);
1005  s[1] = -0.125 * (1. + u) * (1. + w);
1006  s[2] = 0.125 * (1. + u) * (1. - v);
1007  break;
1008  case 6:
1009  s[0] = 0.125 * (1. + v) * (1. + w);
1010  s[1] = 0.125 * (1. + u) * (1. + w);
1011  s[2] = 0.125 * (1. + u) * (1. + v);
1012  break;
1013  case 7:
1014  s[0] = -0.125 * (1. + v) * (1. + w);
1015  s[1] = 0.125 * (1. - u) * (1. + w);
1016  s[2] = 0.125 * (1. - u) * (1. + v);
1017  break;
1018  default: s[0] = s[1] = s[2] = 0.; break;
1019  }
1020  }
1021  int isInside(double u, double v, double w)
1022  {
1023  double TOL = getTolerance();
1024  if(u < -(1. + TOL) || v < -(1. + TOL) || w < -(1. + TOL) ||
1025  u > (1. + TOL) || v > (1. + TOL) || w > (1. + TOL))
1026  return 0;
1027  return 1;
1028  }
1029 };
1030 
1031 class prism : public element {
1032 public:
1033  prism(double *x, double *y, double *z, int numNodes = 0)
1034  : element(x, y, z, numNodes)
1035  {
1036  }
1037  inline int getDimension() { return 3; }
1038  inline int getNumNodes() { return 6; }
1039  void getNode(int num, double &u, double &v, double &w)
1040  {
1041  switch(num) {
1042  case 0:
1043  u = 0.;
1044  v = 0.;
1045  w = -1.;
1046  break;
1047  case 1:
1048  u = 1.;
1049  v = 0.;
1050  w = -1.;
1051  break;
1052  case 2:
1053  u = 0.;
1054  v = 1.;
1055  w = -1.;
1056  break;
1057  case 3:
1058  u = 0.;
1059  v = 0.;
1060  w = 1.;
1061  break;
1062  case 4:
1063  u = 1.;
1064  v = 0.;
1065  w = 1.;
1066  break;
1067  case 5:
1068  u = 0.;
1069  v = 1.;
1070  w = 1.;
1071  break;
1072  default:
1073  u = 0.;
1074  v = 0.;
1075  w = 0.;
1076  break;
1077  }
1078  }
1079  inline int getNumEdges() { return 9; }
1080  void getEdge(int num, int &start, int &end)
1081  {
1082  switch(num) {
1083  case 0:
1084  start = 0;
1085  end = 1;
1086  break;
1087  case 1:
1088  start = 0;
1089  end = 2;
1090  break;
1091  case 2:
1092  start = 0;
1093  end = 3;
1094  break;
1095  case 3:
1096  start = 1;
1097  end = 2;
1098  break;
1099  case 4:
1100  start = 1;
1101  end = 4;
1102  break;
1103  case 5:
1104  start = 2;
1105  end = 5;
1106  break;
1107  case 6:
1108  start = 3;
1109  end = 4;
1110  break;
1111  case 7:
1112  start = 3;
1113  end = 5;
1114  break;
1115  case 8:
1116  start = 4;
1117  end = 5;
1118  break;
1119  default: start = end = 0; break;
1120  }
1121  }
1122  inline int getNumGaussPoints() { return 6; }
1123  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
1124  {
1125  static double u6[6] = {0.166666666666666, 0.333333333333333,
1126  0.166666666666666, 0.166666666666666,
1127  0.333333333333333, 0.166666666666666};
1128  static double v6[6] = {0.166666666666666, 0.166666666666666,
1129  0.333333333333333, 0.166666666666666,
1130  0.166666666666666, 0.333333333333333};
1131  static double w6[6] = {-0.577350269189, -0.577350269189, -0.577350269189,
1132  0.577350269189, 0.577350269189, 0.577350269189};
1133  static double p6[6] = {
1134  0.166666666666666, 0.166666666666666, 0.166666666666666,
1135  0.166666666666666, 0.166666666666666, 0.166666666666666,
1136  };
1137  if(num < 0 || num > 5) return;
1138  u = u6[num];
1139  v = v6[num];
1140  w = w6[num];
1141  weight = p6[num];
1142  }
1143  void getShapeFunction(int num, double u, double v, double w, double &s)
1144  {
1145  switch(num) {
1146  case 0: s = (1. - u - v) * (1. - w) * 0.5; break;
1147  case 1: s = u * (1. - w) * 0.5; break;
1148  case 2: s = v * (1. - w) * 0.5; break;
1149  case 3: s = (1. - u - v) * (1. + w) * 0.5; break;
1150  case 4: s = u * (1. + w) * 0.5; break;
1151  case 5: s = v * (1. + w) * 0.5; break;
1152  default: s = 0.; break;
1153  }
1154  }
1155  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
1156  {
1157  switch(num) {
1158  case 0:
1159  s[0] = -0.5 * (1. - w);
1160  s[1] = -0.5 * (1. - w);
1161  s[2] = -0.5 * (1. - u - v);
1162  break;
1163  case 1:
1164  s[0] = 0.5 * (1. - w);
1165  s[1] = 0.;
1166  s[2] = -0.5 * u;
1167  break;
1168  case 2:
1169  s[0] = 0.;
1170  s[1] = 0.5 * (1. - w);
1171  s[2] = -0.5 * v;
1172  break;
1173  case 3:
1174  s[0] = -0.5 * (1. + w);
1175  s[1] = -0.5 * (1. + w);
1176  s[2] = 0.5 * (1. - u - v);
1177  break;
1178  case 4:
1179  s[0] = 0.5 * (1. + w);
1180  s[1] = 0.;
1181  s[2] = 0.5 * u;
1182  break;
1183  case 5:
1184  s[0] = 0.;
1185  s[1] = 0.5 * (1. + w);
1186  s[2] = 0.5 * v;
1187  break;
1188  default: s[0] = s[1] = s[2] = 0.; break;
1189  }
1190  }
1191  int isInside(double u, double v, double w)
1192  {
1193  double TOL = getTolerance();
1194  if(w > (1. + TOL) || w < -(1. + TOL) || u < (-TOL) || v < (-TOL) ||
1195  u > ((1. + TOL) - v))
1196  return 0;
1197  return 1;
1198  }
1199 };
1200 
1201 class pyramid : public element {
1202 public:
1203  pyramid(double *x, double *y, double *z, int numNodes = 0)
1204  : element(x, y, z, numNodes)
1205  {
1206  }
1207  inline int getDimension() { return 3; }
1208  inline int getNumNodes() { return 5; }
1209  void getNode(int num, double &u, double &v, double &w)
1210  {
1211  switch(num) {
1212  case 0:
1213  u = -1.;
1214  v = -1.;
1215  w = 0.;
1216  break;
1217  case 1:
1218  u = 1.;
1219  v = -1.;
1220  w = 0.;
1221  break;
1222  case 2:
1223  u = 1.;
1224  v = 1.;
1225  w = 0.;
1226  break;
1227  case 3:
1228  u = -1.;
1229  v = 1.;
1230  w = 0.;
1231  break;
1232  case 4:
1233  u = 0.;
1234  v = 0.;
1235  w = 1.;
1236  break;
1237  default:
1238  u = 0.;
1239  v = 0.;
1240  w = 0.;
1241  break;
1242  }
1243  }
1244  inline int getNumEdges() { return 8; }
1245  void getEdge(int num, int &start, int &end)
1246  {
1247  switch(num) {
1248  case 0:
1249  start = 0;
1250  end = 1;
1251  break;
1252  case 1:
1253  start = 0;
1254  end = 3;
1255  break;
1256  case 2:
1257  start = 0;
1258  end = 4;
1259  break;
1260  case 3:
1261  start = 1;
1262  end = 2;
1263  break;
1264  case 4:
1265  start = 1;
1266  end = 4;
1267  break;
1268  case 5:
1269  start = 2;
1270  end = 3;
1271  break;
1272  case 6:
1273  start = 2;
1274  end = 4;
1275  break;
1276  case 7:
1277  start = 3;
1278  end = 4;
1279  break;
1280  default: start = end = 0; break;
1281  }
1282  }
1283  inline int getNumGaussPoints() { return 8; }
1284  void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
1285  {
1286  static double u8[8] = {0.2631840555694285, -0.2631840555694285,
1287  0.2631840555694285, -0.2631840555694285,
1288  0.5066163033492386, -0.5066163033492386,
1289  0.5066163033492386, -0.5066163033492386};
1290  static double v8[8] = {0.2631840555694285, 0.2631840555694285,
1291  -0.2631840555694285, -0.2631840555694285,
1292  0.5066163033492386, 0.5066163033492386,
1293  -0.5066163033492386, -0.5066163033492386};
1294  static double w8[8] = {0.544151844011225, 0.544151844011225,
1295  0.544151844011225, 0.544151844011225,
1296  0.122514822655441, 0.122514822655441,
1297  0.122514822655441, 0.122514822655441};
1298  static double p8[8] = {0.100785882079825, 0.100785882079825,
1299  0.100785882079825, 0.100785882079825,
1300  0.232547451253508, 0.232547451253508,
1301  0.232547451253508, 0.232547451253508};
1302  if(num < 0 || num > 7) return;
1303  u = u8[num];
1304  v = v8[num];
1305  w = w8[num];
1306  weight = p8[num];
1307  }
1308  void getShapeFunction(int num, double u, double v, double w, double &s)
1309  {
1310  double r;
1311 
1312  if(w != 1. && num != 4)
1313  r = u * v * w / (1. - w);
1314  else
1315  r = 0.;
1316  switch(num) {
1317  case 0: s = 0.25 * ((1. - u) * (1. - v) - w + r); break;
1318  case 1: s = 0.25 * ((1. + u) * (1. - v) - w - r); break;
1319  case 2: s = 0.25 * ((1. + u) * (1. + v) - w + r); break;
1320  case 3: s = 0.25 * ((1. - u) * (1. + v) - w - r); break;
1321  case 4: s = w; break;
1322  default: s = 0.; break;
1323  }
1324  }
1325  void getGradShapeFunction(int num, double u, double v, double w, double s[3])
1326  {
1327  if(w == 1.) {
1328  switch(num) {
1329  case 0:
1330  s[0] = -0.25;
1331  s[1] = -0.25;
1332  s[2] = -0.25;
1333  break;
1334  case 1:
1335  s[0] = 0.25;
1336  s[1] = -0.25;
1337  s[2] = -0.25;
1338  break;
1339  case 2:
1340  s[0] = 0.25;
1341  s[1] = 0.25;
1342  s[2] = -0.25;
1343  break;
1344  case 3:
1345  s[0] = -0.25;
1346  s[1] = 0.25;
1347  s[2] = -0.25;
1348  break;
1349  case 4:
1350  s[0] = 0.;
1351  s[1] = 0.;
1352  s[2] = 1.;
1353  break;
1354  default: s[0] = s[1] = s[2] = 0.; break;
1355  }
1356  }
1357  else {
1358  switch(num) {
1359  case 0:
1360  s[0] = 0.25 * (-(1. - v) + v * w / (1. - w));
1361  s[1] = 0.25 * (-(1. - u) + u * w / (1. - w));
1362  s[2] =
1363  0.25 * (-1. + u * v / (1. - w) + u * v * w / std::pow(1. - w, 2));
1364  break;
1365  case 1:
1366  s[0] = 0.25 * ((1. - v) - v * w / (1. - w));
1367  s[1] = 0.25 * (-(1. + u) - u * w / (1. - w));
1368  s[2] =
1369  0.25 * (-1. - u * v / (1. - w) - u * v * w / std::pow(1. - w, 2));
1370  break;
1371  case 2:
1372  s[0] = 0.25 * ((1. + v) + v * w / (1. - w));
1373  s[1] = 0.25 * ((1. + u) + u * w / (1. - w));
1374  s[2] =
1375  0.25 * (-1. + u * v / (1. - w) + u * v * w / std::pow(1. - w, 2));
1376  break;
1377  case 3:
1378  s[0] = 0.25 * (-(1. + v) - v * w / (1. - w));
1379  s[1] = 0.25 * ((1. - u) - u * w / (1. - w));
1380  s[2] =
1381  0.25 * (-1. - u * v / (1. - w) - u * v * w / std::pow(1. - w, 2));
1382  break;
1383  case 4:
1384  s[0] = 0.;
1385  s[1] = 0.;
1386  s[2] = 1.;
1387  break;
1388  default: s[0] = s[1] = s[2] = 0.; break;
1389  }
1390  }
1391  }
1392  int isInside(double u, double v, double w)
1393  {
1394  double TOL = getTolerance();
1395  if(u < (w - (1. + TOL)) || u > ((1. + TOL) - w) || v < (w - (1. + TOL)) ||
1396  v > ((1. + TOL) - w) || w < (-TOL) || w > (1. + TOL))
1397  return 0;
1398  return 1;
1399  }
1400 };
1401 
1402 class elementFactory {
1403 public:
1404  element *create(int numNodes, int dimension, double *x, double *y, double *z,
1405  bool copy = false)
1406  {
1407  switch(dimension) {
1408  case 0: return new point(x, y, z, copy ? numNodes : 0);
1409  case 1: return new line(x, y, z, copy ? numNodes : 0);
1410  case 2:
1411  if(numNodes == 4)
1412  return new quadrangle(x, y, z, copy ? numNodes : 0);
1413  else
1414  return new triangle(x, y, z, copy ? numNodes : 0);
1415  case 3:
1416  if(numNodes == 8)
1417  return new hexahedron(x, y, z, copy ? numNodes : 0);
1418  else if(numNodes == 6)
1419  return new prism(x, y, z, copy ? numNodes : 0);
1420  else if(numNodes == 5)
1421  return new pyramid(x, y, z, copy ? numNodes : 0);
1422  else
1423  return new tetrahedron(x, y, z, copy ? numNodes : 0);
1424  default: Msg::Error("Unknown type of element in factory"); return nullptr;
1425  }
1426  }
1427 };
1428 
1429 #undef SQU
1430 
1431 #endif
line::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:372
element::getXYZ
virtual void getXYZ(int num, double &x, double &y, double &z)
Definition: shapeFunctions.h:46
element::isInside
virtual int isInside(double u, double v, double w)=0
quadrangle::getDimension
int getDimension()
Definition: shapeFunctions.h:560
inv3x3
double inv3x3(double mat[3][3], double inv[3][3])
Definition: Numeric.cpp:232
triangle
Definition: shapeFunctions.h:414
line::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:360
pyramid::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:1208
element::_z
double * _z
Definition: shapeFunctions.h:15
element::maxEdgeLength
double maxEdgeLength()
Definition: shapeFunctions.h:290
point::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:313
element::integrate
double integrate(double val[], int stride=1)
Definition: shapeFunctions.h:212
triangle::getDimension
int getDimension()
Definition: shapeFunctions.h:420
element::element
element(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:18
triangle::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:487
triangle::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:463
quadrangle::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:672
quadrangle::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:635
hexahedron::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:839
quadrangle::integrateFlux
double integrateFlux(double val[])
Definition: shapeFunctions.h:661
prism::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:1155
point::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:314
pyramid::pyramid
pyramid(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:1203
element::getNumGaussPoints
virtual int getNumGaussPoints()=0
c
static double c(int i, int j, fullMatrix< double > &CA, const std::vector< SPoint3 > &P, const std::vector< SPoint3 > &Q)
Definition: discreteFrechetDistance.cpp:15
hexahedron::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:889
point::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:316
Msg::Error
static void Error(const char *fmt,...)
Definition: GmshMessage.cpp:482
triangle::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:421
element::_x
double * _x
Definition: shapeFunctions.h:15
point::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:315
prism::getDimension
int getDimension()
Definition: shapeFunctions.h:1037
LegendrePolynomials::f
void f(int n, double u, double *val)
Definition: orthogonalBasis.cpp:77
hexahedron::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:838
hexahedron::getDimension
int getDimension()
Definition: shapeFunctions.h:837
element::getEdge
virtual void getEdge(int num, int &start, int &end)=0
line::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:404
point::point
point(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:307
quadrangle::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:612
line::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:380
prism::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:1079
tetrahedron::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:689
GmshMessage.h
point::xyz2uvw
void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:333
dimension
int dimension
Definition: GModelIO_TOCHNOG.cpp:18
sys3x3
int sys3x3(double mat[3][3], double b[3], double res[3], double *det)
Definition: Numeric.cpp:147
element::interpolateDiv
double interpolateDiv(double val[], double u, double v, double w, int stride=3)
Definition: shapeFunctions.h:201
prism::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:1191
element::_y
double * _y
Definition: shapeFunctions.h:15
element::getGradShapeFunction
virtual void getGradShapeFunction(int num, double u, double v, double w, double s[3])=0
pyramid::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:1245
hexahedron::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:890
hexahedron::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:945
quadrangle::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:589
elementFactory
Definition: shapeFunctions.h:1402
tetrahedron::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:768
line::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:366
element::getDimension
virtual int getDimension()=0
element::getTolerance
double getTolerance() const
Definition: shapeFunctions.cpp:9
triangle::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:478
prism::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:1122
tetrahedron::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:751
prism::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:1038
hexahedron::hexahedron
hexahedron(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:833
tetrahedron::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:778
quadrangle::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:611
hexahedron::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:944
point::getDimension
int getDimension()
Definition: shapeFunctions.h:311
element::getJacobian
double getJacobian(double u, double v, double w, double jac[3][3])
Definition: shapeFunctions.h:66
tetrahedron::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:721
prism
Definition: shapeFunctions.h:1031
Numeric.h
tetrahedron::getDimension
int getDimension()
Definition: shapeFunctions.h:688
triangle::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:422
point::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:329
line::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:350
line::line
line(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:344
quadrangle::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:562
line::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:359
prism::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:1080
prosca
double prosca(double const a[3], double const b[3])
Definition: Numeric.h:112
hexahedron::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:1021
prism::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:1143
hexahedron::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:975
pyramid::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:1244
triangle::integrateFlux
double integrateFlux(double val[])
Definition: shapeFunctions.h:508
element::~element
virtual ~element()
Definition: shapeFunctions.h:38
pyramid::getGradShapeFunction
void getGradShapeFunction(int num, double u, double v, double w, double s[3])
Definition: shapeFunctions.h:1325
prism::prism
prism(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:1033
elementFactory::create
element * create(int numNodes, int dimension, double *x, double *y, double *z, bool copy=false)
Definition: shapeFunctions.h:1404
pyramid::getDimension
int getDimension()
Definition: shapeFunctions.h:1207
prodve
void prodve(double a[3], double b[3], double c[3])
Definition: Numeric.h:105
tetrahedron::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:690
point::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:312
line::integrateCirculation
double integrateCirculation(double val[])
Definition: shapeFunctions.h:396
element
Definition: shapeFunctions.h:12
element::interpolateGrad
void interpolateGrad(double val[], double u, double v, double w, double f[3], int stride=1, double invjac[3][3]=nullptr)
Definition: shapeFunctions.h:167
quadrangle::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:561
element::getNumNodes
virtual int getNumNodes()=0
tetrahedron::tetrahedron
tetrahedron(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:684
triangle::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:464
pyramid::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:1284
triangle::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:545
tetrahedron::xyz2uvw
void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:804
matvec
void matvec(double mat[3][3], double vec[3], double res[3])
Definition: Numeric.cpp:47
line::getDimension
int getDimension()
Definition: shapeFunctions.h:348
prism::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:1123
element::xyz2uvw
virtual void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:248
z
const double z
Definition: GaussQuadratureQuad.cpp:56
pyramid::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:1283
element::getNumEdges
virtual int getNumEdges()=0
O
#define O
Definition: DefaultOptions.h:22
picojson::copy
void copy(const std::string &s, Iter oi)
Definition: picojson.h:510
tetrahedron::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:752
element::integrateCirculation
virtual double integrateCirculation(double val[])
Definition: shapeFunctions.h:238
tetrahedron
Definition: shapeFunctions.h:682
line::getNumGaussPoints
int getNumGaussPoints()
Definition: shapeFunctions.h:365
pyramid::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:1308
quadrangle::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:625
element::getNode
virtual void getNode(int num, double &u, double &v, double &w)=0
triangle::getEdge
void getEdge(int num, int &start, int &end)
Definition: shapeFunctions.h:445
triangle::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:444
line::getNumNodes
int getNumNodes()
Definition: shapeFunctions.h:349
quadrangle
Definition: shapeFunctions.h:554
element::integrateFlux
virtual double integrateFlux(double val[])
Definition: shapeFunctions.h:243
prism::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:1039
pyramid
Definition: shapeFunctions.h:1201
element::getShapeFunction
virtual void getShapeFunction(int num, double u, double v, double w, double &s)=0
line
Definition: shapeFunctions.h:342
element::integrateLevelsetPositive
double integrateLevelsetPositive(double val[])
Definition: shapeFunctions.h:224
element::interpolate
double interpolate(double val[], double u, double v, double w, int stride=1)
Definition: shapeFunctions.h:155
norme
double norme(double a[3])
Definition: Numeric.h:123
quadrangle::quadrangle
quadrangle(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:556
quadrangle::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:588
tetrahedron::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:822
point::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:322
pyramid::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:1392
pyramid::getNode
void getNode(int num, double &u, double &v, double &w)
Definition: shapeFunctions.h:1209
element::getGaussPoint
virtual void getGaussPoint(int num, double &u, double &v, double &w, double &weight)=0
element::interpolateCurl
void interpolateCurl(double val[], double u, double v, double w, double f[3], int stride=3)
Definition: shapeFunctions.h:188
hexahedron::getShapeFunction
void getShapeFunction(int num, double u, double v, double w, double &s)
Definition: shapeFunctions.h:961
triangle::xyz2uvw
void xyz2uvw(double xyz[3], double uvw[3])
Definition: shapeFunctions.h:519
triangle::triangle
triangle(double *x, double *y, double *z, int numNodes=0)
Definition: shapeFunctions.h:416
hexahedron
Definition: shapeFunctions.h:831
point::isInside
int isInside(double u, double v, double w)
Definition: shapeFunctions.h:334
point::getGaussPoint
void getGaussPoint(int num, double &u, double &v, double &w, double &weight)
Definition: shapeFunctions.h:317
tetrahedron::getNumEdges
int getNumEdges()
Definition: shapeFunctions.h:720
element::_ownData
bool _ownData
Definition: shapeFunctions.h:14