1 /// @file 2 /// Test assembly of mass and Poisson operator QFunction 3 /// \test Test assembly of mass and Poisson operator QFunction 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 #include "t532-operator.h" 8 9 int main(int argc, char **argv) { 10 Ceed ceed; 11 CeedElemRestriction Erestrictx, Erestrictu, 12 Erestrictxi, Erestrictui, 13 Erestrictqi, Erestrictlini; 14 CeedBasis bx, bu; 15 CeedQFunction qf_setup_mass, qf_setup_diff, qf_apply, qf_apply_lin; 16 CeedOperator op_setup_mass, op_setup_diff, op_apply, op_apply_lin; 17 CeedVector qdata_mass, qdata_diff, X, A, u, v; 18 CeedInt nelem = 6, P = 3, Q = 4, dim = 2; 19 CeedInt nx = 3, ny = 2; 20 CeedInt ndofs = (nx*2+1)*(ny*2+1), nqpts = nelem*Q*Q; 21 CeedInt indx[nelem*P*P]; 22 CeedScalar x[dim*ndofs]; 23 24 CeedInit(argv[1], &ceed); 25 26 // DoF Coordinates 27 for (CeedInt i=0; i<nx*2+1; i++) 28 for (CeedInt j=0; j<ny*2+1; j++) { 29 x[i+j*(nx*2+1)+0*ndofs] = (CeedScalar) i / (2*nx); 30 x[i+j*(nx*2+1)+1*ndofs] = (CeedScalar) j / (2*ny); 31 } 32 CeedVectorCreate(ceed, dim*ndofs, &X); 33 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 34 35 // Qdata Vectors 36 CeedVectorCreate(ceed, nqpts, &qdata_mass); 37 CeedVectorCreate(ceed, nqpts*dim*(dim+1)/2, &qdata_diff); 38 39 // Element Setup 40 for (CeedInt i=0; i<nelem; i++) { 41 CeedInt col, row, offset; 42 col = i % nx; 43 row = i / nx; 44 offset = col*(P-1) + row*(nx*2+1)*(P-1); 45 for (CeedInt j=0; j<P; j++) 46 for (CeedInt k=0; k<P; k++) 47 indx[P*(P*i+k)+j] = offset + k*(nx*2+1) + j; 48 } 49 50 // Restrictions 51 CeedElemRestrictionCreate(ceed, nelem, P*P, ndofs, dim, CEED_MEM_HOST, 52 CEED_USE_POINTER, indx, &Erestrictx); 53 CeedElemRestrictionCreateIdentity(ceed, nelem, P*P, nelem*P*P, dim, 54 &Erestrictxi); 55 56 CeedElemRestrictionCreate(ceed, nelem, P*P, ndofs, 1, CEED_MEM_HOST, 57 CEED_USE_POINTER, indx, &Erestrictu); 58 CeedElemRestrictionCreateIdentity(ceed, nelem, Q*Q, nqpts, 1, &Erestrictui); 59 60 CeedElemRestrictionCreateIdentity(ceed, nelem, Q*Q, nqpts, dim*(dim+1)/2, 61 &Erestrictqi); 62 63 // Bases 64 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, P, Q, CEED_GAUSS, &bx); 65 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &bu); 66 67 // QFunction - setup mass 68 CeedQFunctionCreateInterior(ceed, 1, setup_mass, setup_mass_loc, 69 &qf_setup_mass); 70 CeedQFunctionAddInput(qf_setup_mass, "dx", dim*dim, CEED_EVAL_GRAD); 71 CeedQFunctionAddInput(qf_setup_mass, "_weight", 1, CEED_EVAL_WEIGHT); 72 CeedQFunctionAddOutput(qf_setup_mass, "qdata", 1, CEED_EVAL_NONE); 73 74 // Operator - setup mass 75 CeedOperatorCreate(ceed, qf_setup_mass, NULL, NULL, &op_setup_mass); 76 CeedOperatorSetField(op_setup_mass, "dx", Erestrictx, CEED_NOTRANSPOSE, bx, 77 CEED_VECTOR_ACTIVE); 78 CeedOperatorSetField(op_setup_mass, "_weight", Erestrictxi, CEED_NOTRANSPOSE, 79 bx, CEED_VECTOR_NONE); 80 CeedOperatorSetField(op_setup_mass, "qdata", Erestrictui, CEED_NOTRANSPOSE, 81 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 82 83 // QFunction - setup diff 84 CeedQFunctionCreateInterior(ceed, 1, setup_diff, setup_diff_loc, 85 &qf_setup_diff); 86 CeedQFunctionAddInput(qf_setup_diff, "dx", dim*dim, CEED_EVAL_GRAD); 87 CeedQFunctionAddInput(qf_setup_diff, "_weight", 1, CEED_EVAL_WEIGHT); 88 CeedQFunctionAddOutput(qf_setup_diff, "qdata", dim*(dim+1)/2, CEED_EVAL_NONE); 89 90 // Operator - setup diff 91 CeedOperatorCreate(ceed, qf_setup_diff, NULL, NULL, &op_setup_diff); 92 CeedOperatorSetField(op_setup_diff, "dx", Erestrictx, CEED_NOTRANSPOSE, bx, 93 CEED_VECTOR_ACTIVE); 94 CeedOperatorSetField(op_setup_diff, "_weight", Erestrictxi, CEED_NOTRANSPOSE, 95 bx, CEED_VECTOR_NONE); 96 CeedOperatorSetField(op_setup_diff, "qdata", Erestrictqi, CEED_NOTRANSPOSE, 97 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 98 99 // Apply Setup Operators 100 CeedOperatorApply(op_setup_mass, X, qdata_mass, CEED_REQUEST_IMMEDIATE); 101 CeedOperatorApply(op_setup_diff, X, qdata_diff, CEED_REQUEST_IMMEDIATE); 102 103 // QFunction - apply 104 CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply); 105 CeedQFunctionAddInput(qf_apply, "du", dim, CEED_EVAL_GRAD); 106 CeedQFunctionAddInput(qf_apply, "qdata_mass", 1, CEED_EVAL_NONE); 107 CeedQFunctionAddInput(qf_apply, "qdata_diff", dim*(dim+1)/2, CEED_EVAL_NONE); 108 CeedQFunctionAddInput(qf_apply, "u", 1, CEED_EVAL_INTERP); 109 CeedQFunctionAddOutput(qf_apply, "v", 1, CEED_EVAL_INTERP); 110 CeedQFunctionAddOutput(qf_apply, "dv", dim, CEED_EVAL_GRAD); 111 112 // Operator - apply 113 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply); 114 CeedOperatorSetField(op_apply, "du", Erestrictu, CEED_NOTRANSPOSE, bu, 115 CEED_VECTOR_ACTIVE); 116 CeedOperatorSetField(op_apply, "qdata_mass", Erestrictui, CEED_NOTRANSPOSE, 117 CEED_BASIS_COLLOCATED, qdata_mass); 118 CeedOperatorSetField(op_apply, "qdata_diff", Erestrictqi, CEED_NOTRANSPOSE, 119 CEED_BASIS_COLLOCATED, qdata_diff); 120 CeedOperatorSetField(op_apply, "u", Erestrictu, CEED_NOTRANSPOSE, bu, 121 CEED_VECTOR_ACTIVE); 122 CeedOperatorSetField(op_apply, "v", Erestrictu, CEED_NOTRANSPOSE, bu, 123 CEED_VECTOR_ACTIVE); 124 CeedOperatorSetField(op_apply, "dv", Erestrictu, CEED_NOTRANSPOSE, bu, 125 CEED_VECTOR_ACTIVE); 126 127 // Apply original operator 128 CeedVectorCreate(ceed, ndofs, &u); 129 CeedVectorSetValue(u, 1.0); 130 CeedVectorCreate(ceed, ndofs, &v); 131 CeedVectorSetValue(v, 0.0); 132 CeedOperatorApply(op_apply, u, v, CEED_REQUEST_IMMEDIATE); 133 134 // Check output 135 CeedScalar area = 0.0; 136 const CeedScalar *vv; 137 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &vv); 138 for (CeedInt i=0; i<ndofs; i++) 139 area += vv[i]; 140 CeedVectorRestoreArrayRead(v, &vv); 141 if (fabs(area - 1.0) > 1E-14) 142 // LCOV_EXCL_START 143 printf("Error: True operator computed area = %f != 1.0\n", area); 144 // LCOV_EXCL_STOP 145 146 // Assemble QFunction 147 CeedOperatorAssembleLinearQFunction(op_apply, &A, &Erestrictlini, 148 CEED_REQUEST_IMMEDIATE); 149 150 // QFunction - apply assembled 151 CeedQFunctionCreateInterior(ceed, 1, apply_lin, apply_lin_loc, &qf_apply_lin); 152 CeedQFunctionAddInput(qf_apply_lin, "du", dim, CEED_EVAL_GRAD); 153 CeedQFunctionAddInput(qf_apply_lin, "qdata", (dim+1)*(dim+1), CEED_EVAL_NONE); 154 CeedQFunctionAddInput(qf_apply_lin, "u", 1, CEED_EVAL_INTERP); 155 CeedQFunctionAddOutput(qf_apply_lin, "v", 1, CEED_EVAL_INTERP); 156 CeedQFunctionAddOutput(qf_apply_lin, "dv", dim, CEED_EVAL_GRAD); 157 158 // Operator - apply assembled 159 CeedOperatorCreate(ceed, qf_apply_lin, NULL, NULL, &op_apply_lin); 160 CeedOperatorSetField(op_apply_lin, "du", Erestrictu, CEED_NOTRANSPOSE, bu, 161 CEED_VECTOR_ACTIVE); 162 CeedOperatorSetField(op_apply_lin, "qdata", Erestrictlini, CEED_NOTRANSPOSE, 163 CEED_BASIS_COLLOCATED, A); 164 CeedOperatorSetField(op_apply_lin, "u", Erestrictu, CEED_NOTRANSPOSE, bu, 165 CEED_VECTOR_ACTIVE); 166 CeedOperatorSetField(op_apply_lin, "v", Erestrictu, CEED_NOTRANSPOSE, bu, 167 CEED_VECTOR_ACTIVE); 168 CeedOperatorSetField(op_apply_lin, "dv", Erestrictu, CEED_NOTRANSPOSE, bu, 169 CEED_VECTOR_ACTIVE); 170 171 // Apply assembled QFunction operator 172 CeedVectorSetValue(v, 0.0); 173 CeedOperatorApply(op_apply_lin, u, v, CEED_REQUEST_IMMEDIATE); 174 175 // Check output 176 area = 0.0; 177 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &vv); 178 for (CeedInt i=0; i<ndofs; i++) 179 area += vv[i]; 180 CeedVectorRestoreArrayRead(v, &vv); 181 if (fabs(area - 1.0) > 1E-14) 182 // LCOV_EXCL_START 183 printf("Error: Assembled operator computed area = %f != 1.0\n", area); 184 // LCOV_EXCL_STOP 185 186 // Cleanup 187 CeedQFunctionDestroy(&qf_setup_mass); 188 CeedQFunctionDestroy(&qf_setup_diff); 189 CeedQFunctionDestroy(&qf_apply); 190 CeedQFunctionDestroy(&qf_apply_lin); 191 CeedOperatorDestroy(&op_setup_mass); 192 CeedOperatorDestroy(&op_setup_diff); 193 CeedOperatorDestroy(&op_apply); 194 CeedOperatorDestroy(&op_apply_lin); 195 CeedElemRestrictionDestroy(&Erestrictu); 196 CeedElemRestrictionDestroy(&Erestrictx); 197 CeedElemRestrictionDestroy(&Erestrictui); 198 CeedElemRestrictionDestroy(&Erestrictxi); 199 CeedElemRestrictionDestroy(&Erestrictqi); 200 CeedElemRestrictionDestroy(&Erestrictlini); 201 CeedBasisDestroy(&bu); 202 CeedBasisDestroy(&bx); 203 CeedVectorDestroy(&X); 204 CeedVectorDestroy(&A); 205 CeedVectorDestroy(&qdata_mass); 206 CeedVectorDestroy(&qdata_diff); 207 CeedVectorDestroy(&u); 208 CeedVectorDestroy(&v); 209 CeedDestroy(&ceed); 210 return 0; 211 } 212