-
Notifications
You must be signed in to change notification settings - Fork 46
/
t17-basis.c
60 lines (52 loc) · 1.89 KB
/
t17-basis.c
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
// Test grad in multiple dimensions
#include <ceed.h>
#include <math.h>
static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) {
CeedScalar result = tanh(x[0] + 0.1);
if (dim > 1) result += atan(x[1] + 0.2);
if (dim > 2) result += exp(-(x[2] + 0.3)*(x[2] + 0.3));
return result;
}
int main(int argc, char **argv) {
Ceed ceed;
CeedInit(argv[1], &ceed);
for (CeedInt dim=1; dim<=3; dim++) {
CeedBasis bxl, bug;
CeedInt P = 8, Q = 10, Pdim = CeedPowInt(P, dim), Qdim = CeedPowInt(Q, dim),
Xdim = CeedPowInt(2, dim);
CeedScalar x[Xdim*dim], ones[dim*Qdim], gtposeones[Pdim];
CeedScalar xq[Pdim*dim], uq[dim*Qdim], u[Pdim], sum1 = 0, sum2 = 0;
for (CeedInt i=0; i<dim*Qdim; i++) ones[i] = 1;
for (CeedInt d=0; d<dim; d++) {
for (CeedInt i=0; i<Xdim; i++) {
x[d*Xdim + i] = (i % CeedPowInt(2, dim-d)) / CeedPowInt(2, dim-d-1) ? 1 : -1;
}
}
// Get function values at quadrature points
CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, P, CEED_GAUSS_LOBATTO, &bxl);
CeedBasisApply(bxl, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, xq);
for (CeedInt i=0; i<Pdim; i++) {
CeedScalar xx[dim];
for (CeedInt d=0; d<dim; d++) xx[d] = xq[d*Pdim + i%Pdim];
u[i] = Eval(dim, xx);
}
// Calculate G u at quadrature points, G' * 1 at dofs
CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &bug);
CeedBasisApply(bug, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u, uq);
CeedBasisApply(bug, CEED_TRANSPOSE, CEED_EVAL_GRAD, ones, gtposeones);
// Check if 1' * G * u = u' * (G' * 1)
for (CeedInt i=0; i<Pdim; i++) {
sum1 += gtposeones[i]*u[i];
}
for (CeedInt i=0; i<dim*Qdim; i++) {
sum2 += uq[i];
}
if (fabs(sum1 - sum2) > 1e-10) {
printf("[%d] %f != %f\n", dim, sum1, sum2);
}
CeedBasisDestroy(&bxl);
CeedBasisDestroy(&bug);
}
CeedDestroy(&ceed);
return 0;
}