xref: /libCEED/tests/t504-operator.c (revision ec3da8bcb94d9f0073544b37b5081a06981a86f7)
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   for (CeedInt i=0; i<nelem; i++) {
32     for (CeedInt j=0; j<P; j++) {
33       indu[P*i+j] = 2*(i*(P-1) + j);
34     }
35   }
36   CeedElemRestrictionCreate(ceed, nelem, P, 2, 1, 2*Nu, CEED_MEM_HOST,
37                             CEED_USE_POINTER, indu, &Erestrictu);
38   CeedInt stridesu[3] = {1, Q, Q};
39   CeedElemRestrictionCreateStrided(ceed, nelem, Q, 1, Q*nelem, stridesu,
40                                    &Erestrictui);
41 
42   // Bases
43   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, Q, CEED_GAUSS, &bx);
44   CeedBasisCreateTensorH1Lagrange(ceed, 1, 2, P, Q, CEED_GAUSS, &bu);
45 
46   // QFunctions
47   CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
48   CeedQFunctionAddInput(qf_setup, "_weight", 1, CEED_EVAL_WEIGHT);
49   CeedQFunctionAddInput(qf_setup, "dx", 1*1, CEED_EVAL_GRAD);
50   CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
51 
52   CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
53   CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE);
54   CeedQFunctionAddInput(qf_mass, "u", 2, CEED_EVAL_INTERP);
55   CeedQFunctionAddOutput(qf_mass, "v", 2, CEED_EVAL_INTERP);
56 
57   // Operators
58   CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
59                      &op_setup);
60 
61   CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
62                      &op_mass);
63 
64   CeedVectorCreate(ceed, nelem*Q, &qdata);
65 
66   CeedOperatorSetField(op_setup, "_weight", CEED_ELEMRESTRICTION_NONE, bx,
67                        CEED_VECTOR_NONE);
68   CeedOperatorSetField(op_setup, "dx", Erestrictx, bx, CEED_VECTOR_ACTIVE);
69   CeedOperatorSetField(op_setup, "rho", Erestrictui, CEED_BASIS_COLLOCATED,
70                        CEED_VECTOR_ACTIVE);
71 
72   CeedOperatorSetField(op_mass, "rho", Erestrictui, CEED_BASIS_COLLOCATED,
73                        qdata);
74   CeedOperatorSetField(op_mass, "u", Erestrictu, bu, CEED_VECTOR_ACTIVE);
75   CeedOperatorSetField(op_mass, "v", Erestrictu, bu, CEED_VECTOR_ACTIVE);
76 
77   CeedOperatorView(op_setup, stdout);
78   CeedOperatorView(op_mass, stdout);
79 
80   CeedQFunctionDestroy(&qf_setup);
81   CeedQFunctionDestroy(&qf_mass);
82   CeedOperatorDestroy(&op_setup);
83   CeedOperatorDestroy(&op_mass);
84   CeedElemRestrictionDestroy(&Erestrictu);
85   CeedElemRestrictionDestroy(&Erestrictx);
86   CeedElemRestrictionDestroy(&Erestrictui);
87 
88   CeedBasisDestroy(&bu);
89   CeedBasisDestroy(&bx);
90   CeedVectorDestroy(&qdata);
91   CeedDestroy(&ceed);
92   return 0;
93 }
94