xref: /libCEED/tests/t504-operator.c (revision 91dfd1cde504855b06db4052ed3242027f67af19)
1 /// @file
2 /// Test viewing of mass matrix operator
3 /// \test Test viewing of mass matrix operator
4 #include <ceed.h>
5 #include <stdlib.h>
6 #include <math.h>
7 
8 #include "t500-operator.h"
9 
10 int main(int argc, char **argv) {
11   Ceed ceed;
12   CeedElemRestriction Erestrictx, Erestrictu, Erestrictui;
13   CeedBasis bx, bu;
14   CeedQFunction qf_setup, qf_mass;
15   CeedOperator op_setup, op_mass;
16   CeedVector qdata;
17   CeedInt nelem = 15, P = 5, Q = 8;
18   CeedInt Nx = nelem+1, Nu = nelem*(P-1)+1;
19   CeedInt indx[nelem*2], indu[nelem*P];
20 
21   CeedInit(argv[1], &ceed);
22 
23   for (CeedInt i=0; i<nelem; i++) {
24     indx[2*i+0] = i;
25     indx[2*i+1] = i+1;
26   }
27   // Restrictions
28   CeedElemRestrictionCreate(ceed, nelem, 2, 1, 1, Nx, CEED_MEM_HOST,
29                             CEED_USE_POINTER, indx, &Erestrictx);
30 
31 
32   for (CeedInt i=0; i<nelem; i++) {
33     for (CeedInt j=0; j<P; j++) {
34       indu[P*i+j] = 2*(i*(P-1) + j);
35     }
36   }
37   CeedElemRestrictionCreate(ceed, nelem, P, 2, 1, 2*Nu, CEED_MEM_HOST,
38                             CEED_USE_POINTER, indu, &Erestrictu);
39   CeedInt stridesu[3] = {1, Q, Q};
40   CeedElemRestrictionCreateStrided(ceed, nelem, Q, 1, Q*nelem, stridesu,
41                                    &Erestrictui);
42 
43   // Bases
44   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bx);
45   CeedBasisCreateTensorH1Lagrange(ceed, 1, 2, P, Q, CEED_GAUSS, &bu);
46 
47   // QFunctions
48   CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
49   CeedQFunctionAddInput(qf_setup, "_weight", 1, CEED_EVAL_WEIGHT);
50   CeedQFunctionAddInput(qf_setup, "dx", 1*1, CEED_EVAL_GRAD);
51   CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
52 
53   CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
54   CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE);
55   CeedQFunctionAddInput(qf_mass, "u", 2, CEED_EVAL_INTERP);
56   CeedQFunctionAddOutput(qf_mass, "v", 2, CEED_EVAL_INTERP);
57 
58   // Operators
59   CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
60                      &op_setup);
61 
62   CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
63                      &op_mass);
64 
65   CeedVectorCreate(ceed, nelem*Q, &qdata);
66 
67   CeedOperatorSetField(op_setup, "_weight", CEED_ELEMRESTRICTION_NONE, bx,
68                        CEED_VECTOR_NONE);
69   CeedOperatorSetField(op_setup, "dx", Erestrictx, bx, CEED_VECTOR_ACTIVE);
70   CeedOperatorSetField(op_setup, "rho", Erestrictui, CEED_BASIS_COLLOCATED,
71                        CEED_VECTOR_ACTIVE);
72 
73   CeedOperatorSetField(op_mass, "rho", Erestrictui, CEED_BASIS_COLLOCATED,
74                        qdata);
75   CeedOperatorSetField(op_mass, "u", Erestrictu, bu, CEED_VECTOR_ACTIVE);
76   CeedOperatorSetField(op_mass, "v", Erestrictu, bu, CEED_VECTOR_ACTIVE);
77 
78   CeedOperatorView(op_setup, stdout);
79   CeedOperatorView(op_mass, stdout);
80 
81   CeedQFunctionDestroy(&qf_setup);
82   CeedQFunctionDestroy(&qf_mass);
83   CeedOperatorDestroy(&op_setup);
84   CeedOperatorDestroy(&op_mass);
85   CeedElemRestrictionDestroy(&Erestrictu);
86   CeedElemRestrictionDestroy(&Erestrictx);
87   CeedElemRestrictionDestroy(&Erestrictui);
88 
89   CeedBasisDestroy(&bu);
90   CeedBasisDestroy(&bx);
91   CeedVectorDestroy(&qdata);
92   CeedDestroy(&ceed);
93   return 0;
94 }
95