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 = 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