1 /// @file 2 /// Test creation and use of FDM element inverse 3 /// \test Test creation and use of FDM element inverse 4 #include <ceed.h> 5 #include <stdlib.h> 6 #include <math.h> 7 #include "t540-operator.h" 8 9 int main(int argc, char **argv) { 10 Ceed ceed; 11 CeedTransposeMode lmode = CEED_NOTRANSPOSE; 12 CeedElemRestriction Erestrictxi, Erestrictui, Erestrictqi; 13 CeedBasis bx, bu; 14 CeedQFunction qf_setup_mass, qf_apply; 15 CeedOperator op_setup_mass, op_apply, op_inv; 16 CeedVector qdata_mass, X, U, V; 17 CeedInt nelem = 1, P = 4, Q = 5, dim = 2; 18 CeedInt ndofs = P*P, nqpts = nelem*Q*Q; 19 CeedScalar x[dim*nelem*(2*2)]; 20 const CeedScalar *u; 21 22 CeedInit(argv[1], &ceed); 23 24 // DoF Coordinates 25 for (CeedInt i=0; i<2; i++) 26 for (CeedInt j=0; j<2; j++) { 27 x[i+j*2+0*4] = i; 28 x[i+j*2+1*4] = j; 29 } 30 CeedVectorCreate(ceed, dim*nelem*(2*2), &X); 31 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 32 33 // Qdata Vector 34 CeedVectorCreate(ceed, nqpts, &qdata_mass); 35 36 // Element Setup 37 38 // Restrictions 39 CeedElemRestrictionCreateIdentity(ceed, lmode, nelem, 2*2, nelem*2*2, dim, 40 &Erestrictxi); 41 42 CeedElemRestrictionCreateIdentity(ceed, lmode, nelem, P*P, ndofs, 1, 43 &Erestrictui); 44 45 CeedElemRestrictionCreateIdentity(ceed, lmode, nelem, Q*Q, nqpts, 1, 46 &Erestrictqi); 47 48 // Bases 49 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &bx); 50 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &bu); 51 52 // QFunction - setup mass 53 CeedQFunctionCreateInterior(ceed, 1, setup_mass, setup_mass_loc, 54 &qf_setup_mass); 55 CeedQFunctionAddInput(qf_setup_mass, "dx", dim*dim, CEED_EVAL_GRAD); 56 CeedQFunctionAddInput(qf_setup_mass, "_weight", 1, CEED_EVAL_WEIGHT); 57 CeedQFunctionAddOutput(qf_setup_mass, "qdata", 1, CEED_EVAL_NONE); 58 59 // Operator - setup mass 60 CeedOperatorCreate(ceed, qf_setup_mass, CEED_QFUNCTION_NONE, 61 CEED_QFUNCTION_NONE, &op_setup_mass); 62 CeedOperatorSetField(op_setup_mass, "dx", Erestrictxi, bx, 63 CEED_VECTOR_ACTIVE); 64 CeedOperatorSetField(op_setup_mass, "_weight", Erestrictxi, bx, 65 CEED_VECTOR_NONE); 66 CeedOperatorSetField(op_setup_mass, "qdata", Erestrictqi, 67 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 68 69 // Apply Setup Operator 70 CeedOperatorApply(op_setup_mass, X, qdata_mass, CEED_REQUEST_IMMEDIATE); 71 72 // QFunction - apply 73 CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply); 74 CeedQFunctionAddInput(qf_apply, "u", 1, CEED_EVAL_INTERP); 75 CeedQFunctionAddInput(qf_apply, "qdata_mass", 1, CEED_EVAL_NONE); 76 CeedQFunctionAddOutput(qf_apply, "v", 1, CEED_EVAL_INTERP); 77 78 // Operator - apply 79 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, 80 &op_apply); 81 CeedOperatorSetField(op_apply, "u", Erestrictui, bu, CEED_VECTOR_ACTIVE); 82 CeedOperatorSetField(op_apply, "qdata_mass", Erestrictqi, 83 CEED_BASIS_COLLOCATED, qdata_mass); 84 CeedOperatorSetField(op_apply, "v", Erestrictui, bu, CEED_VECTOR_ACTIVE); 85 86 // Apply original operator 87 CeedVectorCreate(ceed, ndofs, &U); 88 CeedVectorSetValue(U, 1.0); 89 CeedVectorCreate(ceed, ndofs, &V); 90 CeedVectorSetValue(V, 0.0); 91 CeedOperatorApply(op_apply, U, V, CEED_REQUEST_IMMEDIATE); 92 93 // Create FDM element inverse 94 CeedOperatorCreateFDMElementInverse(op_apply, &op_inv, CEED_REQUEST_IMMEDIATE); 95 96 // Apply FDM element inverse 97 CeedOperatorApply(op_inv, V, U, CEED_REQUEST_IMMEDIATE); 98 99 // Check output 100 CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); 101 for (int i=0; i<ndofs; i++) 102 if (fabs(u[i] - 1.0) > 1E-14) 103 // LCOV_EXCL_START 104 printf("[%d] Error in inverse: %f != 1.0\n", i, u[i]); 105 // LCOV_EXCL_STOP 106 CeedVectorRestoreArrayRead(U, &u); 107 108 // Cleanup 109 CeedQFunctionDestroy(&qf_setup_mass); 110 CeedQFunctionDestroy(&qf_apply); 111 CeedOperatorDestroy(&op_setup_mass); 112 CeedOperatorDestroy(&op_apply); 113 CeedOperatorDestroy(&op_inv); 114 CeedElemRestrictionDestroy(&Erestrictui); 115 CeedElemRestrictionDestroy(&Erestrictxi); 116 CeedElemRestrictionDestroy(&Erestrictqi); 117 CeedBasisDestroy(&bu); 118 CeedBasisDestroy(&bx); 119 CeedVectorDestroy(&X); 120 CeedVectorDestroy(&qdata_mass); 121 CeedVectorDestroy(&U); 122 CeedVectorDestroy(&V); 123 CeedDestroy(&ceed); 124 return 0; 125 } 126