1 /// @file
2 /// Test creation, action, and destruction for mass matrix operator
3 /// \test Test creation, action, and destruction for mass matrix operator
4 #include "t500-operator.h"
5
6 #include <ceed.h>
7 #include <math.h>
8 #include <stdio.h>
9 #include <stdlib.h>
10
main(int argc,char ** argv)11 int main(int argc, char **argv) {
12 Ceed ceed;
13 CeedElemRestriction elem_restriction_x, elem_restriction_u, elem_restriction_q_data;
14 CeedBasis basis_x, basis_u;
15 CeedQFunction qf_setup, qf_mass;
16 CeedOperator op_setup, op_mass;
17 CeedVector q_data, x, u, v;
18 CeedInt num_elem = 15, p = 5, q = 8;
19 CeedInt num_nodes_x = num_elem + 1, num_nodes_u = num_elem * (p - 1) + 1;
20 CeedInt ind_x[num_elem * 2], ind_u[num_elem * p];
21 CeedScalar x_array[num_nodes_x];
22
23 //! [Ceed Init]
24 CeedInit(argv[1], &ceed);
25 //! [Ceed Init]
26 for (CeedInt i = 0; i < num_nodes_x; i++) x_array[i] = (CeedScalar)i / (num_nodes_x - 1);
27 for (CeedInt i = 0; i < num_elem; i++) {
28 ind_x[2 * i + 0] = i;
29 ind_x[2 * i + 1] = i + 1;
30 }
31 //! [ElemRestr Create]
32 CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_nodes_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
33 //! [ElemRestr Create]
34
35 for (CeedInt i = 0; i < num_elem; i++) {
36 for (CeedInt j = 0; j < p; j++) {
37 ind_u[p * i + j] = i * (p - 1) + j;
38 }
39 }
40 //! [ElemRestrU Create]
41 CeedElemRestrictionCreate(ceed, num_elem, p, 1, 1, num_nodes_u, CEED_MEM_HOST, CEED_USE_POINTER, ind_u, &elem_restriction_u);
42 CeedInt strides_q_data[3] = {1, q, q};
43 CeedElemRestrictionCreateStrided(ceed, num_elem, q, 1, q * num_elem, strides_q_data, &elem_restriction_q_data);
44 //! [ElemRestrU Create]
45
46 //! [Basis Create]
47 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x);
48 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u);
49 //! [Basis Create]
50
51 //! [QFunction Create]
52 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
53 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
54 CeedQFunctionAddInput(qf_setup, "dx", 1, CEED_EVAL_GRAD);
55 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
56
57 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
58 CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE);
59 CeedQFunctionAddInput(qf_mass, "u", 1, CEED_EVAL_INTERP);
60 CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP);
61 //! [QFunction Create]
62
63 //! [Setup Create]
64 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
65 //! [Setup Create]
66
67 //! [Operator Create]
68 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass);
69 //! [Operator Create]
70
71 CeedVectorCreate(ceed, num_nodes_x, &x);
72 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_USE_POINTER, x_array);
73 CeedVectorCreate(ceed, num_elem * q, &q_data);
74
75 //! [Setup Set]
76 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
77 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
78 CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
79 //! [Setup Set]
80
81 //! [Operator Set]
82 CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
83 CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
84 CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
85 //! [Operator Set]
86
87 //! [Setup Apply]
88 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
89 //! [Setup Apply]
90
91 CeedVectorCreate(ceed, num_nodes_u, &u);
92 CeedVectorSetValue(u, 0.0);
93 CeedVectorCreate(ceed, num_nodes_u, &v);
94 //! [Operator Apply]
95 CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE);
96 //! [Operator Apply]
97
98 {
99 const CeedScalar *v_array;
100
101 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
102 for (CeedInt i = 0; i < num_nodes_u; i++) {
103 if (fabs(v_array[i]) > 1e-14) printf("[%" CeedInt_FMT "] v %g != 0.0\n", i, v_array[i]);
104 }
105 CeedVectorRestoreArrayRead(v, &v_array);
106 }
107
108 CeedVectorDestroy(&x);
109 CeedVectorDestroy(&u);
110 CeedVectorDestroy(&v);
111 CeedVectorDestroy(&q_data);
112 CeedElemRestrictionDestroy(&elem_restriction_x);
113 CeedElemRestrictionDestroy(&elem_restriction_u);
114 CeedElemRestrictionDestroy(&elem_restriction_q_data);
115 CeedBasisDestroy(&basis_x);
116 CeedBasisDestroy(&basis_u);
117 CeedQFunctionDestroy(&qf_setup);
118 CeedQFunctionDestroy(&qf_mass);
119 CeedOperatorDestroy(&op_setup);
120 CeedOperatorDestroy(&op_mass);
121 CeedDestroy(&ceed);
122 return 0;
123 }
124