xref: /libCEED/tests/t509-operator.c (revision 2b730f8b5a9c809740a0b3b302db43a719c636b1)
10b454692Sjeremylt /// @file
20b454692Sjeremylt /// Test creation, action, and destruction for identity operator
30b454692Sjeremylt /// \test Test creation, action, and destruction for identity operator
40b454692Sjeremylt #include <ceed.h>
50b454692Sjeremylt #include <math.h>
6*2b730f8bSJeremy L Thompson #include <stdlib.h>
70b454692Sjeremylt 
80b454692Sjeremylt int main(int argc, char **argv) {
90b454692Sjeremylt   Ceed                ceed;
100b454692Sjeremylt   CeedElemRestriction elem_restr_u, elem_restr_u_i;
110b454692Sjeremylt   CeedBasis           basis_u;
120b454692Sjeremylt   CeedQFunction       qf_identity;
130b454692Sjeremylt   CeedOperator        op_identity;
140b454692Sjeremylt   CeedVector          U, V;
150b454692Sjeremylt   const CeedScalar   *hv;
160b454692Sjeremylt   CeedInt             num_elem = 15, P = 5, Q = 8;
170b454692Sjeremylt   CeedInt             elem_size = P, num_nodes = num_elem * (P - 1) + 1;
180b454692Sjeremylt   CeedInt             ind_u[num_elem * P];
190b454692Sjeremylt 
200b454692Sjeremylt   CeedInit(argv[1], &ceed);
210b454692Sjeremylt 
220b454692Sjeremylt   // Restrictions
230b454692Sjeremylt   for (CeedInt i = 0; i < num_elem; i++) {
240b454692Sjeremylt     for (CeedInt j = 0; j < P; j++) {
250b454692Sjeremylt       ind_u[P * i + j] = i * (P - 1) + j;
260b454692Sjeremylt     }
270b454692Sjeremylt   }
28*2b730f8bSJeremy L Thompson   CeedElemRestrictionCreate(ceed, num_elem, elem_size, 1, 1, num_nodes, CEED_MEM_HOST, CEED_USE_POINTER, ind_u, &elem_restr_u);
290b454692Sjeremylt   CeedInt strides_u_i[3] = {1, P, P};
30*2b730f8bSJeremy L Thompson   CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, 1, elem_size * num_elem, strides_u_i, &elem_restr_u_i);
310b454692Sjeremylt 
320b454692Sjeremylt   // Bases
330b454692Sjeremylt   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, Q, CEED_GAUSS, &basis_u);
340b454692Sjeremylt 
350b454692Sjeremylt   // QFunction
36*2b730f8bSJeremy L Thompson   CeedQFunctionCreateIdentity(ceed, 1, CEED_EVAL_NONE, CEED_EVAL_NONE, &qf_identity);
370b454692Sjeremylt 
380b454692Sjeremylt   // Operators
39*2b730f8bSJeremy L Thompson   CeedOperatorCreate(ceed, qf_identity, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_identity);
40*2b730f8bSJeremy L Thompson   CeedOperatorSetField(op_identity, "input", elem_restr_u, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
41*2b730f8bSJeremy L Thompson   CeedOperatorSetField(op_identity, "output", elem_restr_u_i, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
420b454692Sjeremylt   CeedOperatorSetNumQuadraturePoints(op_identity, elem_size);
430b454692Sjeremylt 
440b454692Sjeremylt   CeedVectorCreate(ceed, num_nodes, &U);
450b454692Sjeremylt   CeedVectorSetValue(U, 3.0);
460b454692Sjeremylt   CeedVectorCreate(ceed, elem_size * num_elem, &V);
470b454692Sjeremylt   CeedOperatorApply(op_identity, U, V, CEED_REQUEST_IMMEDIATE);
480b454692Sjeremylt 
490b454692Sjeremylt   // Check output
500b454692Sjeremylt   CeedVectorGetArrayRead(V, CEED_MEM_HOST, &hv);
51*2b730f8bSJeremy L Thompson   for (CeedInt i = 0; i < elem_size * num_elem; i++) {
52*2b730f8bSJeremy L Thompson     if (fabs(hv[i] - 3.) > 1e-14) printf("[%" CeedInt_FMT "] Computed Value: %f != True Value: 1.0\n", i, hv[i]);
53*2b730f8bSJeremy L Thompson   }
540b454692Sjeremylt   CeedVectorRestoreArrayRead(V, &hv);
550b454692Sjeremylt 
560b454692Sjeremylt   CeedQFunctionDestroy(&qf_identity);
570b454692Sjeremylt   CeedOperatorDestroy(&op_identity);
580b454692Sjeremylt   CeedElemRestrictionDestroy(&elem_restr_u);
590b454692Sjeremylt   CeedElemRestrictionDestroy(&elem_restr_u_i);
600b454692Sjeremylt   CeedBasisDestroy(&basis_u);
610b454692Sjeremylt   CeedVectorDestroy(&U);
620b454692Sjeremylt   CeedVectorDestroy(&V);
630b454692Sjeremylt   CeedDestroy(&ceed);
640b454692Sjeremylt   return 0;
650b454692Sjeremylt }
66