xref: /libCEED/tests/t319-basis.c (revision db2becc9f302fe8eb3a32ace50ce3f3a5d42e6c4) 
1 /// @file
2 /// Test projection interp and grad in multiple dimensions
3 /// \test Test projection interp and grad in multiple dimensions
4 #include "t319-basis.h"
5 #include <ceed.h>
6 #include <math.h>
7 #include <stdio.h>
8 
9 static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) {
10   CeedScalar result = (x[0] + 0.1) * (x[0] + 0.1);
11   if (dim > 1) result += (x[1] + 0.2) * (x[1] + 0.2);
12   if (dim > 2) result += -(x[2] + 0.3) * (x[2] + 0.3);
13   return result;
14 }
15 
16 static CeedScalar EvalGrad(CeedInt dim, const CeedScalar x[]) {
17   switch (dim) {
18     case 0:
19       return 2 * x[0] + 0.2;
20     case 1:
21       return 2 * x[1] + 0.4;
22     default:
23       return -2 * x[2] - 0.6;
24   }
25 }
26 
27 static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) {
28   CeedScalar tol;
29   if (scalar_type == CEED_SCALAR_FP32) {
30     if (dim == 3) tol = 1.e-4;
31     else tol = 1.e-5;
32   } else {
33     tol = 1.e-11;
34   }
35   return tol;
36 }
37 
38 static void VerifyProjectedBasis(CeedBasis basis_project, CeedInt dim, CeedInt p_to_dim, CeedInt p_from_dim, CeedVector x_to, CeedVector x_from,
39                                  CeedVector u_to, CeedVector u_from, CeedVector du_to) {
40   CeedScalar tol;
41 
42   {
43     CeedScalarType scalar_type;
44 
45     CeedGetScalarType(&scalar_type);
46     tol = GetTolerance(scalar_type, dim);
47   }
48 
49   // Setup coarse solution
50   {
51     const CeedScalar *x_array;
52     CeedScalar        u_array[p_from_dim];
53 
54     CeedVectorGetArrayRead(x_from, CEED_MEM_HOST, &x_array);
55     for (CeedInt i = 0; i < p_from_dim; i++) {
56       CeedScalar coord[dim];
57       for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[p_from_dim * d + i];
58       u_array[i] = Eval(dim, coord);
59     }
60     CeedVectorRestoreArrayRead(x_from, &x_array);
61     CeedVectorSetArray(u_from, CEED_MEM_HOST, CEED_COPY_VALUES, u_array);
62   }
63 
64   // Project to fine basis
65   CeedBasisApply(basis_project, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u_from, u_to);
66 
67   // Check solution
68   {
69     const CeedScalar *x_array, *u_array;
70 
71     CeedVectorGetArrayRead(x_to, CEED_MEM_HOST, &x_array);
72     CeedVectorGetArrayRead(u_to, CEED_MEM_HOST, &u_array);
73     for (CeedInt i = 0; i < p_to_dim; i++) {
74       CeedScalar coord[dim];
75       for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_to_dim + i];
76       const CeedScalar u = Eval(dim, coord);
77       if (fabs(u - u_array[i]) > tol) printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f\n", dim, i, u_array[i], u);
78     }
79     CeedVectorRestoreArrayRead(x_to, &x_array);
80     CeedVectorRestoreArrayRead(u_to, &u_array);
81   }
82 
83   // Project and take gradient
84   CeedBasisApply(basis_project, 1, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, u_from, du_to);
85 
86   // Check solution
87   {
88     const CeedScalar *x_array, *du_array;
89 
90     CeedVectorGetArrayRead(x_to, CEED_MEM_HOST, &x_array);
91     CeedVectorGetArrayRead(du_to, CEED_MEM_HOST, &du_array);
92     for (CeedInt i = 0; i < p_to_dim; i++) {
93       CeedScalar coord[dim];
94 
95       for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[p_to_dim * d + i];
96       for (CeedInt d = 0; d < dim; d++) {
97         const CeedScalar du = EvalGrad(d, coord);
98 
99         if (fabs(du - du_array[p_to_dim * d + i]) > tol) {
100           // LCOV_EXCL_START
101           printf("[%" CeedInt_FMT ", %" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f\n", dim, i, d, du_array[p_to_dim * (dim - 1 - d) + i], du);
102           // LCOV_EXCL_STOP
103         }
104       }
105     }
106     CeedVectorRestoreArrayRead(x_to, &x_array);
107     CeedVectorRestoreArrayRead(du_to, &du_array);
108   }
109 }
110 
111 int main(int argc, char **argv) {
112   Ceed ceed;
113 
114   CeedInit(argv[1], &ceed);
115 
116   for (CeedInt dim = 1; dim <= 3; dim++) {
117     CeedVector x_corners, x_from, x_to, u_from, u_to, du_to;
118     CeedBasis  basis_x, basis_from, basis_to, basis_project;
119     CeedInt    p_from = 5, p_to = 6, q = 7, x_dim = CeedIntPow(2, dim), p_from_dim = CeedIntPow(p_from, dim), p_to_dim = CeedIntPow(p_to, dim);
120 
121     CeedVectorCreate(ceed, x_dim * dim, &x_corners);
122     {
123       CeedScalar x_array[x_dim * dim];
124 
125       for (CeedInt d = 0; d < dim; d++) {
126         for (CeedInt i = 0; i < x_dim; i++) x_array[x_dim * d + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1;
127       }
128       CeedVectorSetArray(x_corners, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
129     }
130     CeedVectorCreate(ceed, p_from_dim * dim, &x_from);
131     CeedVectorCreate(ceed, p_to_dim * dim, &x_to);
132     CeedVectorCreate(ceed, p_from_dim, &u_from);
133     CeedVectorSetValue(u_from, 0);
134     CeedVectorCreate(ceed, p_to_dim, &u_to);
135     CeedVectorSetValue(u_to, 0);
136     CeedVectorCreate(ceed, p_to_dim * dim, &du_to);
137     CeedVectorSetValue(du_to, 0);
138 
139     // Get nodal coordinates
140     CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p_from, CEED_GAUSS_LOBATTO, &basis_x);
141     CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_corners, x_from);
142     CeedBasisDestroy(&basis_x);
143     CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p_to, CEED_GAUSS_LOBATTO, &basis_x);
144     CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_corners, x_to);
145     CeedBasisDestroy(&basis_x);
146 
147     // Create U and projection bases
148     CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_from, q, CEED_GAUSS, &basis_from);
149     CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p_to, q, CEED_GAUSS, &basis_to);
150     CeedBasisCreateProjection(basis_from, basis_to, &basis_project);
151 
152     VerifyProjectedBasis(basis_project, dim, p_to_dim, p_from_dim, x_to, x_from, u_to, u_from, du_to);
153 
154     // Create non-tensor bases
155     CeedBasis basis_from_nontensor, basis_to_nontensor;
156     {
157       CeedElemTopology  topo;
158       CeedInt           num_comp, num_nodes, nqpts;
159       const CeedScalar *interp, *grad;
160 
161       CeedBasisGetTopology(basis_from, &topo);
162       CeedBasisGetNumComponents(basis_from, &num_comp);
163       CeedBasisGetNumNodes(basis_from, &num_nodes);
164       CeedBasisGetNumQuadraturePoints(basis_from, &nqpts);
165       CeedBasisGetInterp(basis_from, &interp);
166       CeedBasisGetGrad(basis_from, &grad);
167       CeedBasisCreateH1(ceed, topo, num_comp, num_nodes, nqpts, interp, grad, NULL, NULL, &basis_from_nontensor);
168 
169       CeedBasisGetTopology(basis_to, &topo);
170       CeedBasisGetNumComponents(basis_to, &num_comp);
171       CeedBasisGetNumNodes(basis_to, &num_nodes);
172       CeedBasisGetNumQuadraturePoints(basis_to, &nqpts);
173       CeedBasisGetInterp(basis_to, &interp);
174       CeedBasisGetGrad(basis_to, &grad);
175       CeedBasisCreateH1(ceed, topo, num_comp, num_nodes, nqpts, interp, grad, NULL, NULL, &basis_to_nontensor);
176     }
177 
178     // Test projection on non-tensor bases
179     CeedBasisDestroy(&basis_project);
180     CeedBasisCreateProjection(basis_from_nontensor, basis_to_nontensor, &basis_project);
181     VerifyProjectedBasis(basis_project, dim, p_to_dim, p_from_dim, x_to, x_from, u_to, u_from, du_to);
182 
183     // Test projection from non-tensor to tensor
184     CeedBasisDestroy(&basis_project);
185     CeedBasisCreateProjection(basis_from_nontensor, basis_to, &basis_project);
186     VerifyProjectedBasis(basis_project, dim, p_to_dim, p_from_dim, x_to, x_from, u_to, u_from, du_to);
187 
188     // Test projection from tensor to non-tensor
189     CeedBasisDestroy(&basis_project);
190     CeedBasisCreateProjection(basis_from, basis_to_nontensor, &basis_project);
191     VerifyProjectedBasis(basis_project, dim, p_to_dim, p_from_dim, x_to, x_from, u_to, u_from, du_to);
192 
193     CeedVectorDestroy(&x_corners);
194     CeedVectorDestroy(&x_from);
195     CeedVectorDestroy(&x_to);
196     CeedVectorDestroy(&u_from);
197     CeedVectorDestroy(&u_to);
198     CeedVectorDestroy(&du_to);
199     CeedBasisDestroy(&basis_from);
200     CeedBasisDestroy(&basis_from_nontensor);
201     CeedBasisDestroy(&basis_to);
202     CeedBasisDestroy(&basis_to_nontensor);
203     CeedBasisDestroy(&basis_project);
204   }
205 
206   // Test projection between basis of different topological dimension
207   {
208     CeedInt   face_dim = 2, P_1D = 2;
209     CeedBasis basis_face, basis_cell_to_face, basis_proj;
210 
211     CeedScalar       *q_ref = NULL, *q_weights = NULL;
212     const CeedScalar *grad, *interp;
213     CeedInt           P, Q;
214     GetCellToFaceTabulation(CEED_GAUSS, &P, &Q, &interp, &grad);
215 
216     CeedBasisCreateTensorH1Lagrange(ceed, face_dim, 1, 2, P_1D, CEED_GAUSS, &basis_face);
217     CeedBasisCreateH1(ceed, CEED_TOPOLOGY_HEX, 1, P, Q, (CeedScalar *)interp, (CeedScalar *)grad, q_ref, q_weights, &basis_cell_to_face);
218     CeedBasisCreateProjection(basis_cell_to_face, basis_face, &basis_proj);
219     const CeedScalar *interp_proj, *grad_proj, *interp_proj_ref, *grad_proj_ref;
220 
221     GetCellToFaceTabulation(CEED_GAUSS_LOBATTO, NULL, NULL, &interp_proj_ref, &grad_proj_ref);
222     CeedBasisGetInterp(basis_proj, &interp_proj);
223     CeedBasisGetGrad(basis_proj, &grad_proj);
224     CeedScalar tol = 100 * CEED_EPSILON;
225 
226     for (CeedInt i = 0; i < 4 * 8; i++) {
227       if (fabs(interp_proj[i] - ((CeedScalar *)interp_proj_ref)[i]) > tol)
228         printf("Mixed Topology Projection: interp[%" CeedInt_FMT "] expected %f, got %f\n", i, interp_proj[i], ((CeedScalar *)interp_proj_ref)[i]);
229     }
230 
231     for (CeedInt i = 0; i < 3 * 4 * 8; i++) {
232       if (fabs(grad_proj[i] - ((CeedScalar *)grad_proj_ref)[i]) > tol)
233         printf("Mixed Topology Projection: grad[%" CeedInt_FMT "] expected %f, got %f\n", i, grad_proj[i], ((CeedScalar *)grad_proj_ref)[i]);
234     }
235 
236     CeedBasisDestroy(&basis_face);
237     CeedBasisDestroy(&basis_cell_to_face);
238     CeedBasisDestroy(&basis_proj);
239   }
240   CeedDestroy(&ceed);
241   return 0;
242 }
243