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