1 /// @file 2 /// Test creation and use of FDM element inverse 3 /// \test Test creation and use of FDM element inverse 4 #include "t540-operator.h" 5 6 #include <ceed.h> 7 #include <math.h> 8 #include <stdlib.h> 9 10 int main(int argc, char **argv) { 11 Ceed ceed; 12 CeedElemRestriction elem_restr_x_i, elem_restr_u_i, elem_restr_qd_i; 13 CeedBasis basis_x, basis_u; 14 CeedQFunction qf_setup_mass, qf_apply; 15 CeedOperator op_setup_mass, op_apply, op_inv; 16 CeedVector q_data_mass, X, U, V; 17 CeedInt num_elem = 1, P = 4, Q = 5, dim = 2; 18 CeedInt num_dofs = P * P, num_qpts = num_elem * Q * Q; 19 CeedScalar x[dim * num_elem * (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 } 31 CeedVectorCreate(ceed, dim * num_elem * (2 * 2), &X); 32 CeedVectorSetArray(X, CEED_MEM_HOST, CEED_USE_POINTER, x); 33 34 // Qdata Vector 35 CeedVectorCreate(ceed, num_qpts, &q_data_mass); 36 37 // Element Setup 38 39 // Restrictions 40 CeedInt strides_x[3] = {1, 2 * 2, 2 * 2 * dim}; 41 CeedElemRestrictionCreateStrided(ceed, num_elem, 2 * 2, dim, dim * num_elem * 2 * 2, strides_x, &elem_restr_x_i); 42 43 CeedInt strides_u[3] = {1, P * P, P * P}; 44 CeedElemRestrictionCreateStrided(ceed, num_elem, P * P, 1, num_dofs, strides_u, &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, &elem_restr_qd_i); 48 49 // Bases 50 CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, Q, CEED_GAUSS, &basis_x); 51 CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &basis_u); 52 53 // QFunction - setup mass 54 CeedQFunctionCreateInterior(ceed, 1, setup_mass, setup_mass_loc, &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, CEED_QFUNCTION_NONE, &op_setup_mass); 61 CeedOperatorSetField(op_setup_mass, "dx", elem_restr_x_i, basis_x, CEED_VECTOR_ACTIVE); 62 CeedOperatorSetField(op_setup_mass, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE); 63 CeedOperatorSetField(op_setup_mass, "qdata", elem_restr_qd_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); 64 65 // Apply Setup Operator 66 CeedOperatorApply(op_setup_mass, X, q_data_mass, CEED_REQUEST_IMMEDIATE); 67 68 // QFunction - apply 69 CeedQFunctionCreateInterior(ceed, 1, apply, apply_loc, &qf_apply); 70 CeedQFunctionAddInput(qf_apply, "u", 1, CEED_EVAL_INTERP); 71 CeedQFunctionAddInput(qf_apply, "mass qdata", 1, CEED_EVAL_NONE); 72 CeedQFunctionAddOutput(qf_apply, "v", 1, CEED_EVAL_INTERP); 73 74 // Operator - apply 75 CeedOperatorCreate(ceed, qf_apply, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_apply); 76 CeedOperatorSetField(op_apply, "u", elem_restr_u_i, basis_u, CEED_VECTOR_ACTIVE); 77 CeedOperatorSetField(op_apply, "mass qdata", elem_restr_qd_i, CEED_BASIS_COLLOCATED, q_data_mass); 78 CeedOperatorSetField(op_apply, "v", elem_restr_u_i, basis_u, CEED_VECTOR_ACTIVE); 79 80 // Apply original operator 81 CeedVectorCreate(ceed, num_dofs, &U); 82 CeedVectorSetValue(U, 1.0); 83 CeedVectorCreate(ceed, num_dofs, &V); 84 CeedVectorSetValue(V, 0.0); 85 CeedOperatorApply(op_apply, U, V, CEED_REQUEST_IMMEDIATE); 86 87 // Create FDM element inverse 88 CeedOperatorCreateFDMElementInverse(op_apply, &op_inv, CEED_REQUEST_IMMEDIATE); 89 90 // Apply FDM element inverse 91 CeedOperatorApply(op_inv, V, U, CEED_REQUEST_IMMEDIATE); 92 93 // Check output 94 CeedVectorGetArrayRead(U, CEED_MEM_HOST, &u); 95 for (int i = 0; i < num_dofs; i++) { 96 if (fabs(u[i] - 1.0) > 500. * CEED_EPSILON) printf("[%" CeedInt_FMT "] Error in inverse: %e - 1.0 = %e\n", i, u[i], u[i] - 1.); 97 } 98 CeedVectorRestoreArrayRead(U, &u); 99 100 // Cleanup 101 CeedQFunctionDestroy(&qf_setup_mass); 102 CeedQFunctionDestroy(&qf_apply); 103 CeedOperatorDestroy(&op_setup_mass); 104 CeedOperatorDestroy(&op_apply); 105 CeedOperatorDestroy(&op_inv); 106 CeedElemRestrictionDestroy(&elem_restr_u_i); 107 CeedElemRestrictionDestroy(&elem_restr_x_i); 108 CeedElemRestrictionDestroy(&elem_restr_qd_i); 109 CeedBasisDestroy(&basis_u); 110 CeedBasisDestroy(&basis_x); 111 CeedVectorDestroy(&X); 112 CeedVectorDestroy(&q_data_mass); 113 CeedVectorDestroy(&U); 114 CeedVectorDestroy(&V); 115 CeedDestroy(&ceed); 116 return 0; 117 } 118