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