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