1 /// @file
2 /// Test creation, action, and destruction for mass matrix operator with multiple components
3 /// \test Test creation, action, and destruction for mass matrix operator with multiple components
4 #include "t502-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
22 CeedInit(argv[1], &ceed);
23
24 CeedVectorCreate(ceed, num_nodes_x, &x);
25 {
26 CeedScalar x_array[num_nodes_x];
27
28 for (CeedInt i = 0; i < num_nodes_x; i++) x_array[i] = (CeedScalar)i / (num_nodes_x - 1);
29 CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
30 }
31 CeedVectorCreate(ceed, 2 * num_nodes_u, &u);
32 CeedVectorCreate(ceed, 2 * num_nodes_u, &v);
33 CeedVectorCreate(ceed, num_elem * q, &q_data);
34
35 // Restrictions
36 for (CeedInt i = 0; i < num_elem; i++) {
37 ind_x[2 * i + 0] = i;
38 ind_x[2 * i + 1] = i + 1;
39 }
40 CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_nodes_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restriction_x);
41
42 for (CeedInt i = 0; i < num_elem; i++) {
43 for (CeedInt j = 0; j < p; j++) {
44 ind_u[p * i + j] = 2 * (i * (p - 1) + j);
45 }
46 }
47 CeedElemRestrictionCreate(ceed, num_elem, p, 2, 1, 2 * num_nodes_u, CEED_MEM_HOST, CEED_USE_POINTER, ind_u, &elem_restriction_u);
48
49 CeedInt strides_q_data[3] = {1, q, q};
50 CeedElemRestrictionCreateStrided(ceed, num_elem, q, 1, q * num_elem, strides_q_data, &elem_restriction_q_data);
51
52 // Bases
53 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, q, CEED_GAUSS, &basis_x);
54 CeedBasisCreateTensorH1Lagrange(ceed, 1, 2, p, q, CEED_GAUSS, &basis_u);
55
56 // QFunctions
57 CeedQFunctionCreateInterior(ceed, 1, setup, setup_loc, &qf_setup);
58 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
59 CeedQFunctionAddInput(qf_setup, "dx", 1 * 1, CEED_EVAL_GRAD);
60 CeedQFunctionAddOutput(qf_setup, "rho", 1, CEED_EVAL_NONE);
61
62 CeedQFunctionCreateInterior(ceed, 1, mass, mass_loc, &qf_mass);
63 CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE);
64 CeedQFunctionAddInput(qf_mass, "u", 2, CEED_EVAL_INTERP);
65 CeedQFunctionAddOutput(qf_mass, "v", 2, CEED_EVAL_INTERP);
66
67 // Operators
68 CeedOperatorCreate(ceed, qf_setup, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_setup);
69 CeedOperatorSetField(op_setup, "weight", CEED_ELEMRESTRICTION_NONE, basis_x, CEED_VECTOR_NONE);
70 CeedOperatorSetField(op_setup, "dx", elem_restriction_x, basis_x, CEED_VECTOR_ACTIVE);
71 CeedOperatorSetField(op_setup, "rho", elem_restriction_q_data, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE);
72
73 CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass);
74 CeedOperatorSetField(op_mass, "rho", elem_restriction_q_data, CEED_BASIS_NONE, q_data);
75 CeedOperatorSetField(op_mass, "u", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
76 CeedOperatorSetField(op_mass, "v", elem_restriction_u, basis_u, CEED_VECTOR_ACTIVE);
77
78 CeedOperatorApply(op_setup, x, q_data, CEED_REQUEST_IMMEDIATE);
79
80 {
81 CeedScalar *u_array;
82
83 CeedVectorGetArrayWrite(u, CEED_MEM_HOST, &u_array);
84 for (int i = 0; i < num_nodes_u; i++) {
85 u_array[2 * i] = 1.0;
86 u_array[2 * i + 1] = 2.0;
87 }
88 CeedVectorRestoreArray(u, &u_array);
89 }
90 CeedOperatorApply(op_mass, u, v, CEED_REQUEST_IMMEDIATE);
91
92 // Check output
93 {
94 const CeedScalar *v_array;
95 CeedScalar sum_1 = 0., sum_2 = 0.;
96
97 CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
98 for (CeedInt i = 0; i < num_nodes_u; i++) {
99 sum_1 += v_array[2 * i];
100 sum_2 += v_array[2 * i + 1];
101 }
102 CeedVectorRestoreArrayRead(v, &v_array);
103 if (fabs(sum_1 - 1.) > 1000. * CEED_EPSILON) printf("Computed Area: %f != True Area: 1.0\n", sum_1);
104 if (fabs(sum_2 - 2.) > 1000. * CEED_EPSILON) printf("Computed Area: %f != True Area: 2.0\n", sum_2);
105 }
106
107 CeedVectorDestroy(&x);
108 CeedVectorDestroy(&u);
109 CeedVectorDestroy(&v);
110 CeedVectorDestroy(&q_data);
111 CeedElemRestrictionDestroy(&elem_restriction_u);
112 CeedElemRestrictionDestroy(&elem_restriction_x);
113 CeedElemRestrictionDestroy(&elem_restriction_q_data);
114 CeedBasisDestroy(&basis_u);
115 CeedBasisDestroy(&basis_x);
116 CeedQFunctionDestroy(&qf_setup);
117 CeedQFunctionDestroy(&qf_mass);
118 CeedOperatorDestroy(&op_setup);
119 CeedOperatorDestroy(&op_mass);
120 CeedDestroy(&ceed);
121 return 0;
122 }
123