1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
|
// Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
#include <math.h>
#include "GaussIntegration.h"
#include "GaussLegendre1D.h"
static void brickToTet(double xi, double eta, double zeta, double *r, double *s,
double *t, double *J)
{
double r1, rs1;
*r = 0.5e0 * (1.0e0 + xi);
r1 = 1.0e0 - (*r);
*s = 0.5e0 * (1.0e0 + eta) * r1;
rs1 = 1.0e0 - (*r) - (*s);
*t = 0.5e0 * (1.0e0 + zeta) * rs1;
*J = 0.125e0 * r1 * rs1;
}
void quadToTri(double xi, double eta, double *r, double *s, double *J)
{
double r1;
*r = 0.5e0 * (1.0e0 + xi);
r1 = 1.0e0 - (*r);
*s = 0.5e0 * (1.0e0 + eta) * r1;
*J = 0.25e0 * r1;
}
double quadToTriJac(double, double eta) { return 0.125e0 * (1.0e0 - eta); }
int GaussLegendreTet(int n1, int n2, int n3, IntPt *pts)
{
/* get degenerate n1Xn2Xn3 Gauss-Legendre scheme to integrate over a tet */
int i, j, k, index = 0;
double *pt1, *pt2, *pt3, *wt1, *wt2, *wt3, dJ;
gmshGaussLegendre1D(n1, &pt1, &wt1);
gmshGaussLegendre1D(n2, &pt2, &wt2);
gmshGaussLegendre1D(n3, &pt3, &wt3);
for(i = 0; i < n1; i++) {
for(j = 0; j < n2; j++) {
for(k = 0; k < n3; k++) {
brickToTet(pt1[i], pt2[j], pt3[k], &(pts[index].pt[0]),
&(pts[index].pt[1]), &(pts[index].pt[2]), &dJ);
pts[index++].weight = dJ * wt1[i] * wt2[j] * wt3[k];
}
}
}
return index;
}
int GaussLegendreTri(int n1, int n2, IntPt *pts)
{
/* get degenerate n1Xn2 Gauss-Legendre scheme to integrate over a tri */
int i, j, index = 0;
double *pt1, *pt2, *wt1, *wt2, dJ;
// const double two = 2.0000000000000000;
gmshGaussLegendre1D(n1, &pt1, &wt1);
gmshGaussLegendre1D(n2, &pt2, &wt2);
for(i = 0; i < n1; i++) {
for(j = 0; j < n2; j++) {
quadToTri(pt1[i], pt2[j], &(pts[index].pt[0]), &(pts[index].pt[1]), &dJ);
pts[index].pt[2] = 0;
pts[index++].weight = dJ * wt1[i] * wt2[j];
}
}
return index;
}
|
