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
Eval(CeedInt dim,const CeedScalar x[])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
EvalGrad(CeedInt dim,const CeedScalar x[])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
GetTolerance(CeedScalarType scalar_type,int dim)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
VerifyProjectedBasis(CeedBasis basis_project,CeedInt dim,CeedInt p_to_dim,CeedInt p_from_dim,CeedVector x_to,CeedVector x_from,CeedVector u_to,CeedVector u_from,CeedVector du_to)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
main(int argc,char ** argv)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 = 4, p_to = 5, q = 6, 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, num_qpts;
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, &num_qpts);
165 CeedBasisGetInterp(basis_from, &interp);
166 CeedBasisGetGrad(basis_from, &grad);
167 CeedBasisCreateH1(ceed, topo, num_comp, num_nodes, num_qpts, 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, &num_qpts);
173 CeedBasisGetInterp(basis_to, &interp);
174 CeedBasisGetGrad(basis_to, &grad);
175 CeedBasisCreateH1(ceed, topo, num_comp, num_nodes, num_qpts, 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 // LCOV_EXCL_START
229 printf("Mixed Topology Projection: interp[%" CeedInt_FMT "] expected %f, got %f\n", i, interp_proj[i], ((CeedScalar *)interp_proj_ref)[i]);
230 // LCOV_EXCL_STOP
231 }
232 }
233
234 for (CeedInt i = 0; i < 3 * 4 * 8; i++) {
235 if (fabs(grad_proj[i] - ((CeedScalar *)grad_proj_ref)[i]) > tol) {
236 // LCOV_EXCL_START
237 printf("Mixed Topology Projection: grad[%" CeedInt_FMT "] expected %f, got %f\n", i, grad_proj[i], ((CeedScalar *)grad_proj_ref)[i]);
238 // LCOV_EXCL_STOP
239 }
240 }
241
242 CeedBasisDestroy(&basis_face);
243 CeedBasisDestroy(&basis_cell_to_face);
244 CeedBasisDestroy(&basis_proj);
245 }
246 CeedDestroy(&ceed);
247 return 0;
248 }
249