1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 #include <ceed.h>
9 #include <ceed/backend.h>
10 #include <stdbool.h>
11 #include <stddef.h>
12 #include <stdint.h>
13
14 #include "ceed-ref.h"
15
16 //------------------------------------------------------------------------------
17 // Setup Input/Output Fields
18 //------------------------------------------------------------------------------
CeedOperatorSetupFields_Ref(CeedQFunction qf,CeedOperator op,bool is_input,bool * skip_rstr,CeedInt * e_data_out_indices,bool * apply_add_basis,CeedVector * e_vecs_full,CeedVector * e_vecs,CeedVector * q_vecs,CeedInt start_e,CeedInt num_fields,CeedInt Q)19 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool is_input, bool *skip_rstr, CeedInt *e_data_out_indices,
20 bool *apply_add_basis, CeedVector *e_vecs_full, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e,
21 CeedInt num_fields, CeedInt Q) {
22 Ceed ceed;
23 CeedSize e_size, q_size;
24 CeedInt num_comp, size, P;
25 CeedQFunctionField *qf_fields;
26 CeedOperatorField *op_fields;
27
28 {
29 Ceed ceed_parent;
30
31 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
32 CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
33 CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed));
34 CeedCallBackend(CeedDestroy(&ceed_parent));
35 }
36 if (is_input) {
37 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
38 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
39 } else {
40 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
41 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
42 }
43
44 // Loop over fields
45 for (CeedInt i = 0; i < num_fields; i++) {
46 CeedEvalMode eval_mode;
47 CeedElemRestriction elem_rstr;
48 CeedBasis basis;
49
50 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
51 if (eval_mode != CEED_EVAL_WEIGHT) {
52 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr));
53 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e]));
54 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
55 }
56
57 switch (eval_mode) {
58 case CEED_EVAL_NONE:
59 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
60 q_size = (CeedSize)Q * size;
61 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
62 break;
63 case CEED_EVAL_INTERP:
64 case CEED_EVAL_GRAD:
65 case CEED_EVAL_DIV:
66 case CEED_EVAL_CURL:
67 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
68 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
69 CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
70 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
71 e_size = (CeedSize)P * num_comp;
72 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i]));
73 q_size = (CeedSize)Q * size;
74 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
75 CeedCallBackend(CeedBasisDestroy(&basis));
76 break;
77 case CEED_EVAL_WEIGHT: // Only on input fields
78 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
79 q_size = (CeedSize)Q;
80 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
81 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i]));
82 CeedCallBackend(CeedBasisDestroy(&basis));
83 break;
84 }
85 }
86 // Drop duplicate restrictions
87 if (is_input) {
88 for (CeedInt i = 0; i < num_fields; i++) {
89 CeedVector vec_i;
90 CeedElemRestriction rstr_i;
91
92 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i));
93 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i));
94 for (CeedInt j = i + 1; j < num_fields; j++) {
95 CeedVector vec_j;
96 CeedElemRestriction rstr_j;
97
98 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j));
99 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j));
100 if (vec_i == vec_j && rstr_i == rstr_j) {
101 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j]));
102 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e]));
103 skip_rstr[j] = true;
104 }
105 CeedCallBackend(CeedVectorDestroy(&vec_j));
106 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
107 }
108 CeedCallBackend(CeedVectorDestroy(&vec_i));
109 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
110 }
111 } else {
112 for (CeedInt i = num_fields - 1; i >= 0; i--) {
113 CeedVector vec_i;
114 CeedElemRestriction rstr_i;
115
116 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i));
117 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i));
118 for (CeedInt j = i - 1; j >= 0; j--) {
119 CeedVector vec_j;
120 CeedElemRestriction rstr_j;
121
122 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j));
123 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j));
124 if (vec_i == vec_j && rstr_i == rstr_j) {
125 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j]));
126 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e]));
127 skip_rstr[j] = true;
128 apply_add_basis[i] = true;
129 e_data_out_indices[j] = i;
130 }
131 CeedCallBackend(CeedVectorDestroy(&vec_j));
132 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
133 }
134 CeedCallBackend(CeedVectorDestroy(&vec_i));
135 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
136 }
137 }
138 CeedCallBackend(CeedDestroy(&ceed));
139 return CEED_ERROR_SUCCESS;
140 }
141
142 //------------------------------------------------------------------------------
143 // Setup Operator
144 //------------------------------------------------------------------------------/*
CeedOperatorSetup_Ref(CeedOperator op)145 static int CeedOperatorSetup_Ref(CeedOperator op) {
146 bool is_setup_done;
147 CeedInt Q, num_input_fields, num_output_fields;
148 CeedQFunctionField *qf_input_fields, *qf_output_fields;
149 CeedQFunction qf;
150 CeedOperatorField *op_input_fields, *op_output_fields;
151 CeedOperator_Ref *impl;
152
153 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done));
154 if (is_setup_done) return CEED_ERROR_SUCCESS;
155
156 CeedCallBackend(CeedOperatorGetData(op, &impl));
157 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
158 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
159 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf));
160 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
161 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
162
163 // Allocate
164 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full));
165
166 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in));
167 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_out));
168 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_data_out_indices));
169 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->apply_add_basis_out));
170 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states));
171 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in));
172 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out));
173 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in));
174 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out));
175
176 impl->num_inputs = num_input_fields;
177 impl->num_outputs = num_output_fields;
178
179 // Set up infield and outfield e_vecs and q_vecs
180 // Infields
181 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, true, impl->skip_rstr_in, NULL, NULL, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0,
182 num_input_fields, Q));
183 // Outfields
184 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, false, impl->skip_rstr_out, impl->e_data_out_indices, impl->apply_add_basis_out,
185 impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q));
186
187 // Identity QFunctions
188 if (impl->is_identity_qf) {
189 CeedEvalMode in_mode, out_mode;
190 CeedQFunctionField *in_fields, *out_fields;
191
192 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields));
193 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode));
194 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode));
195
196 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) {
197 impl->is_identity_rstr_op = true;
198 } else {
199 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0]));
200 }
201 }
202
203 CeedCallBackend(CeedOperatorSetSetupDone(op));
204 CeedCallBackend(CeedQFunctionDestroy(&qf));
205 return CEED_ERROR_SUCCESS;
206 }
207
208 //------------------------------------------------------------------------------
209 // Setup Operator Inputs
210 //------------------------------------------------------------------------------
CeedOperatorSetupInputs_Ref(CeedInt num_input_fields,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,CeedVector in_vec,const bool skip_active,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Ref * impl,CeedRequest * request)211 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
212 CeedVector in_vec, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX],
213 CeedOperator_Ref *impl, CeedRequest *request) {
214 for (CeedInt i = 0; i < num_input_fields; i++) {
215 bool is_active;
216 uint64_t state;
217 CeedEvalMode eval_mode;
218 CeedVector vec;
219
220 // Get input vector
221 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
222 is_active = vec == CEED_VECTOR_ACTIVE;
223 if (is_active) {
224 if (skip_active) continue;
225 else vec = in_vec;
226 }
227
228 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
229 // Restrict and Evec
230 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip
231 } else {
232 // Restrict
233 CeedCallBackend(CeedVectorGetState(vec, &state));
234 // Skip restriction if input is unchanged
235 if ((state != impl->input_states[i] || vec == in_vec) && !impl->skip_rstr_in[i]) {
236 CeedElemRestriction elem_rstr;
237
238 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
239 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request));
240 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
241 }
242 impl->input_states[i] = state;
243 // Get evec
244 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i]));
245 }
246 if (!is_active) CeedCallBackend(CeedVectorDestroy(&vec));
247 }
248 return CEED_ERROR_SUCCESS;
249 }
250
251 //------------------------------------------------------------------------------
252 // Input Basis Action
253 //------------------------------------------------------------------------------
CeedOperatorInputBasis_Ref(CeedInt e,CeedInt Q,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,CeedInt num_input_fields,const bool skip_active,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Ref * impl)254 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
255 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX],
256 CeedOperator_Ref *impl) {
257 for (CeedInt i = 0; i < num_input_fields; i++) {
258 CeedInt elem_size, size, num_comp;
259 CeedEvalMode eval_mode;
260 CeedElemRestriction elem_rstr;
261 CeedBasis basis;
262
263 // Skip active input
264 if (skip_active) {
265 bool is_active;
266 CeedVector vec;
267
268 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
269 is_active = vec == CEED_VECTOR_ACTIVE;
270 CeedCallBackend(CeedVectorDestroy(&vec));
271 if (is_active) continue;
272 }
273 // Get elem_size, eval_mode, size
274 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
275 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
276 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
277 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
278 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size));
279 // Basis action
280 switch (eval_mode) {
281 case CEED_EVAL_NONE:
282 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * Q * size]));
283 break;
284 case CEED_EVAL_INTERP:
285 case CEED_EVAL_GRAD:
286 case CEED_EVAL_DIV:
287 case CEED_EVAL_CURL:
288 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
289 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
290 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * elem_size * num_comp]));
291 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i]));
292 CeedCallBackend(CeedBasisDestroy(&basis));
293 break;
294 case CEED_EVAL_WEIGHT:
295 break; // No action
296 }
297 }
298 return CEED_ERROR_SUCCESS;
299 }
300
301 //------------------------------------------------------------------------------
302 // Output Basis Action
303 //------------------------------------------------------------------------------
CeedOperatorOutputBasis_Ref(CeedInt e,CeedInt Q,CeedQFunctionField * qf_output_fields,CeedOperatorField * op_output_fields,CeedInt num_input_fields,CeedInt num_output_fields,bool * apply_add_basis,CeedOperator op,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Ref * impl)304 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields,
305 CeedInt num_input_fields, CeedInt num_output_fields, bool *apply_add_basis, CeedOperator op,
306 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) {
307 for (CeedInt i = 0; i < num_output_fields; i++) {
308 CeedInt elem_size, num_comp;
309 CeedEvalMode eval_mode;
310 CeedElemRestriction elem_rstr;
311 CeedBasis basis;
312
313 // Get elem_size, eval_mode
314 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
315 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
316 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
317 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
318 // Basis action
319 switch (eval_mode) {
320 case CEED_EVAL_NONE:
321 break; // No action
322 case CEED_EVAL_INTERP:
323 case CEED_EVAL_GRAD:
324 case CEED_EVAL_DIV:
325 case CEED_EVAL_CURL:
326 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
327 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
328 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER,
329 &e_data_full[i + num_input_fields][(CeedSize)e * elem_size * num_comp]));
330 if (apply_add_basis[i]) {
331 CeedCallBackend(CeedBasisApplyAdd(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
332 } else {
333 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i]));
334 }
335 CeedCallBackend(CeedBasisDestroy(&basis));
336 break;
337 // LCOV_EXCL_START
338 case CEED_EVAL_WEIGHT: {
339 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
340 // LCOV_EXCL_STOP
341 }
342 }
343 }
344 return CEED_ERROR_SUCCESS;
345 }
346
347 //------------------------------------------------------------------------------
348 // Restore Input Vectors
349 //------------------------------------------------------------------------------
CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,const bool skip_active,CeedScalar * e_data_full[2* CEED_FIELD_MAX],CeedOperator_Ref * impl)350 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields,
351 const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) {
352 for (CeedInt i = 0; i < num_input_fields; i++) {
353 CeedEvalMode eval_mode;
354
355 // Skip active inputs
356 if (skip_active) {
357 bool is_active;
358 CeedVector vec;
359
360 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
361 is_active = vec == CEED_VECTOR_ACTIVE;
362 CeedCallBackend(CeedVectorDestroy(&vec));
363 if (is_active) continue;
364 }
365 // Restore input
366 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
367 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip
368 } else {
369 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i]));
370 }
371 }
372 return CEED_ERROR_SUCCESS;
373 }
374
375 //------------------------------------------------------------------------------
376 // Operator Apply
377 //------------------------------------------------------------------------------
CeedOperatorApplyAdd_Ref(CeedOperator op,CeedVector in_vec,CeedVector out_vec,CeedRequest * request)378 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) {
379 CeedInt Q, num_elem, num_input_fields, num_output_fields, size;
380 CeedEvalMode eval_mode;
381 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {NULL};
382 CeedQFunctionField *qf_input_fields, *qf_output_fields;
383 CeedQFunction qf;
384 CeedOperatorField *op_input_fields, *op_output_fields;
385 CeedOperator_Ref *impl;
386
387 // Setup
388 CeedCallBackend(CeedOperatorSetup_Ref(op));
389
390 CeedCallBackend(CeedOperatorGetData(op, &impl));
391 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
392
393 // Restriction only operator
394 if (impl->is_identity_rstr_op) {
395 CeedElemRestriction elem_rstr;
396
397 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_rstr));
398 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request));
399 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
400 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_rstr));
401 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request));
402 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
403 return CEED_ERROR_SUCCESS;
404 }
405
406 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
407 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
408 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
409 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
410
411 // Input Evecs and Restriction
412 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request));
413
414 // Output Evecs
415 for (CeedInt i = num_output_fields - 1; i >= 0; i--) {
416 if (impl->skip_rstr_out[i]) {
417 e_data_full[i + num_input_fields] = e_data_full[impl->e_data_out_indices[i] + num_input_fields];
418 } else {
419 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields]));
420 }
421 }
422
423 // Loop through elements
424 for (CeedInt e = 0; e < num_elem; e++) {
425 // Output pointers
426 for (CeedInt i = 0; i < num_output_fields; i++) {
427 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
428 if (eval_mode == CEED_EVAL_NONE) {
429 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size));
430 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER,
431 &e_data_full[i + num_input_fields][(CeedSize)e * Q * size]));
432 }
433 }
434
435 // Input basis apply
436 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, false, e_data_full, impl));
437
438 // Q function
439 if (!impl->is_identity_qf) {
440 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out));
441 }
442
443 // Output basis apply
444 CeedCallBackend(CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, num_input_fields, num_output_fields,
445 impl->apply_add_basis_out, op, e_data_full, impl));
446 }
447
448 // Output restriction
449 for (CeedInt i = 0; i < num_output_fields; i++) {
450 bool is_active;
451 CeedVector vec;
452 CeedElemRestriction elem_rstr;
453
454 if (impl->skip_rstr_out[i]) continue;
455 // Restore Evec
456 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields]));
457 // Get output vector
458 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
459 // Active
460 is_active = vec == CEED_VECTOR_ACTIVE;
461 if (is_active) vec = out_vec;
462 // Restrict
463 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
464 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec, request));
465 if (!is_active) CeedCallBackend(CeedVectorDestroy(&vec));
466 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
467 }
468
469 // Restore input arrays
470 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl));
471 CeedCallBackend(CeedQFunctionDestroy(&qf));
472 return CEED_ERROR_SUCCESS;
473 }
474
475 //------------------------------------------------------------------------------
476 // Core code for assembling linear QFunction
477 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op,bool build_objects,CeedVector * assembled,CeedElemRestriction * rstr,CeedRequest * request)478 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr,
479 CeedRequest *request) {
480 Ceed ceed_parent;
481 CeedInt qf_size_in, qf_size_out, Q, num_elem, num_input_fields, num_output_fields;
482 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL};
483 CeedQFunctionField *qf_input_fields, *qf_output_fields;
484 CeedQFunction qf;
485 CeedOperatorField *op_input_fields, *op_output_fields;
486 CeedOperator_Ref *impl;
487
488 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent));
489 CeedCallBackend(CeedOperatorGetData(op, &impl));
490 qf_size_in = impl->qf_size_in;
491 qf_size_out = impl->qf_size_out;
492 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
493 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
494 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
495 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
496 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
497
498 // Setup
499 CeedCallBackend(CeedOperatorSetup_Ref(op));
500
501 // Check for restriction only operator
502 CeedCheck(!impl->is_identity_rstr_op, CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported");
503
504 // Input Evecs and Restriction
505 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request));
506
507 // Count number of active input fields
508 if (qf_size_in == 0) {
509 for (CeedInt i = 0; i < num_input_fields; i++) {
510 CeedInt field_size;
511 CeedVector vec;
512
513 // Get input vector
514 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
515 // Check if active input
516 if (vec == CEED_VECTOR_ACTIVE) {
517 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size));
518 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
519 qf_size_in += field_size;
520 }
521 CeedCallBackend(CeedVectorDestroy(&vec));
522 }
523 CeedCheck(qf_size_in > 0, CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
524 impl->qf_size_in = qf_size_in;
525 }
526
527 // Count number of active output fields
528 if (qf_size_out == 0) {
529 for (CeedInt i = 0; i < num_output_fields; i++) {
530 CeedInt field_size;
531 CeedVector vec;
532
533 // Get output vector
534 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
535 // Check if active output
536 if (vec == CEED_VECTOR_ACTIVE) {
537 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size));
538 qf_size_out += field_size;
539 }
540 CeedCallBackend(CeedVectorDestroy(&vec));
541 }
542 CeedCheck(qf_size_out > 0, CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
543 impl->qf_size_out = qf_size_out;
544 }
545
546 // Build objects if needed
547 if (build_objects) {
548 const CeedSize l_size = (CeedSize)num_elem * Q * qf_size_in * qf_size_out;
549 CeedInt strides[3] = {1, Q, qf_size_in * qf_size_out * Q}; /* *NOPAD* */
550
551 // Create output restriction
552 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, qf_size_in * qf_size_out,
553 (CeedSize)qf_size_in * (CeedSize)qf_size_out * (CeedSize)num_elem * (CeedSize)Q, strides, rstr));
554 // Create assembled vector
555 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled));
556 }
557 // Clear output vector
558 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0));
559 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array));
560
561 // Loop through elements
562 for (CeedInt e = 0; e < num_elem; e++) {
563 // Input basis apply
564 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl));
565
566 // Assemble QFunction
567
568 for (CeedInt i = 0; i < num_input_fields; i++) {
569 bool is_active;
570 CeedInt field_size;
571 CeedVector vec;
572
573 // Set Inputs
574 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
575 is_active = vec == CEED_VECTOR_ACTIVE;
576 CeedCallBackend(CeedVectorDestroy(&vec));
577 if (!is_active) continue;
578 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size));
579 for (CeedInt field = 0; field < field_size; field++) {
580 // Set current portion of input to 1.0
581 {
582 CeedScalar *array;
583
584 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array));
585 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 1.0;
586 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array));
587 }
588
589 if (!impl->is_identity_qf) {
590 // Set Outputs
591 for (CeedInt out = 0; out < num_output_fields; out++) {
592 CeedVector vec;
593
594 // Get output vector
595 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
596 // Check if active output
597 if (vec == CEED_VECTOR_ACTIVE) {
598 CeedInt field_size;
599
600 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array));
601 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size));
602 assembled_array += field_size * Q; // Advance the pointer by the size of the output
603 }
604 CeedCallBackend(CeedVectorDestroy(&vec));
605 }
606 // Apply QFunction
607 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out));
608 } else {
609 CeedInt field_size;
610 const CeedScalar *array;
611
612 // Copy Identity Outputs
613 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size));
614 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array));
615 for (CeedInt j = 0; j < field_size * Q; j++) assembled_array[j] = array[j];
616 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array));
617 assembled_array += field_size * Q;
618 }
619 // Reset input to 0.0
620 {
621 CeedScalar *array;
622
623 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array));
624 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 0.0;
625 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array));
626 }
627 }
628 }
629 }
630
631 // Un-set output Qvecs to prevent accidental overwrite of Assembled
632 if (!impl->is_identity_qf) {
633 for (CeedInt out = 0; out < num_output_fields; out++) {
634 CeedVector vec;
635
636 // Get output vector
637 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
638 // Check if active output
639 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) {
640 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL));
641 }
642 CeedCallBackend(CeedVectorDestroy(&vec));
643 }
644 }
645
646 // Restore input arrays
647 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl));
648
649 // Restore output
650 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array));
651 CeedCallBackend(CeedDestroy(&ceed_parent));
652 CeedCallBackend(CeedQFunctionDestroy(&qf));
653 return CEED_ERROR_SUCCESS;
654 }
655
656 //------------------------------------------------------------------------------
657 // Assemble Linear QFunction
658 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op,CeedVector * assembled,CeedElemRestriction * rstr,CeedRequest * request)659 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
660 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request);
661 }
662
663 //------------------------------------------------------------------------------
664 // Update Assembled Linear QFunction
665 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op,CeedVector assembled,CeedElemRestriction rstr,CeedRequest * request)666 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) {
667 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request);
668 }
669
670 //------------------------------------------------------------------------------
671 // Setup Input/Output Fields
672 //------------------------------------------------------------------------------
CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf,CeedOperator op,bool is_input,bool * skip_rstr,bool * apply_add_basis,CeedVector * e_vecs_full,CeedVector * e_vecs,CeedVector * q_vecs,CeedInt start_e,CeedInt num_fields,CeedInt Q)673 static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, bool *skip_rstr, bool *apply_add_basis,
674 CeedVector *e_vecs_full, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields,
675 CeedInt Q) {
676 Ceed ceed;
677 CeedSize e_size, q_size;
678 CeedInt max_num_points, num_comp, size, P;
679 CeedQFunctionField *qf_fields;
680 CeedOperatorField *op_fields;
681
682 {
683 Ceed ceed_parent;
684
685 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
686 CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
687 CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed));
688 CeedCallBackend(CeedDestroy(&ceed_parent));
689 }
690 if (is_input) {
691 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL));
692 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL));
693 } else {
694 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields));
695 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields));
696 }
697
698 // Get max number of points
699 {
700 CeedInt dim;
701 CeedElemRestriction rstr_points = NULL;
702 CeedOperator_Ref *impl;
703
704 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL));
705 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points));
706 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim));
707 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
708 CeedCallBackend(CeedOperatorGetData(op, &impl));
709 if (is_input) {
710 CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem));
711 CeedCallBackend(CeedVectorSetValue(impl->point_coords_elem, 0.0));
712 }
713 }
714
715 // Loop over fields
716 for (CeedInt i = 0; i < num_fields; i++) {
717 CeedEvalMode eval_mode;
718 CeedBasis basis;
719
720 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode));
721 if (eval_mode != CEED_EVAL_WEIGHT) {
722 CeedElemRestriction elem_rstr;
723
724 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr));
725 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e]));
726 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
727 CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0));
728 }
729
730 switch (eval_mode) {
731 case CEED_EVAL_NONE: {
732 CeedVector vec;
733
734 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
735 e_size = (CeedSize)max_num_points * size;
736 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i]));
737 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec));
738 if (vec == CEED_VECTOR_ACTIVE || !is_input) {
739 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i]));
740 } else {
741 q_size = (CeedSize)max_num_points * size;
742 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
743 }
744 CeedCallBackend(CeedVectorDestroy(&vec));
745 break;
746 }
747 case CEED_EVAL_INTERP:
748 case CEED_EVAL_GRAD:
749 case CEED_EVAL_DIV:
750 case CEED_EVAL_CURL:
751 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
752 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size));
753 CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
754 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
755 e_size = (CeedSize)P * num_comp;
756 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i]));
757 q_size = (CeedSize)max_num_points * size;
758 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
759 CeedCallBackend(CeedBasisDestroy(&basis));
760 break;
761 case CEED_EVAL_WEIGHT: // Only on input fields
762 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis));
763 q_size = (CeedSize)max_num_points;
764 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i]));
765 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE,
766 q_vecs[i]));
767 CeedCallBackend(CeedBasisDestroy(&basis));
768 break;
769 }
770 // Initialize full arrays for E-vectors and Q-vectors
771 if (e_vecs[i]) CeedCallBackend(CeedVectorSetValue(e_vecs[i], 0.0));
772 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0));
773 }
774 // Drop duplicate restrictions
775 if (is_input) {
776 for (CeedInt i = 0; i < num_fields; i++) {
777 CeedVector vec_i;
778 CeedElemRestriction rstr_i;
779
780 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i));
781 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i));
782 for (CeedInt j = i + 1; j < num_fields; j++) {
783 CeedVector vec_j;
784 CeedElemRestriction rstr_j;
785
786 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j));
787 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j));
788 if (vec_i == vec_j && rstr_i == rstr_j) {
789 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j]));
790 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e]));
791 skip_rstr[j] = true;
792 }
793 CeedCallBackend(CeedVectorDestroy(&vec_j));
794 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
795 }
796 CeedCallBackend(CeedVectorDestroy(&vec_i));
797 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
798 }
799 } else {
800 for (CeedInt i = num_fields - 1; i >= 0; i--) {
801 CeedVector vec_i;
802 CeedElemRestriction rstr_i;
803
804 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec_i));
805 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr_i));
806 for (CeedInt j = i - 1; j >= 0; j--) {
807 CeedVector vec_j;
808 CeedElemRestriction rstr_j;
809
810 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[j], &vec_j));
811 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[j], &rstr_j));
812 if (vec_i == vec_j && rstr_i == rstr_j) {
813 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &e_vecs[j]));
814 CeedCallBackend(CeedVectorReferenceCopy(e_vecs_full[i + start_e], &e_vecs_full[j + start_e]));
815 skip_rstr[j] = true;
816 apply_add_basis[i] = true;
817 }
818 CeedCallBackend(CeedVectorDestroy(&vec_j));
819 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_j));
820 }
821 CeedCallBackend(CeedVectorDestroy(&vec_i));
822 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_i));
823 }
824 }
825 CeedCallBackend(CeedDestroy(&ceed));
826 return CEED_ERROR_SUCCESS;
827 }
828
829 //------------------------------------------------------------------------------
830 // Setup Operator
831 //------------------------------------------------------------------------------
CeedOperatorSetupAtPoints_Ref(CeedOperator op)832 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) {
833 bool is_setup_done;
834 CeedInt Q, num_input_fields, num_output_fields;
835 CeedQFunctionField *qf_input_fields, *qf_output_fields;
836 CeedQFunction qf;
837 CeedOperatorField *op_input_fields, *op_output_fields;
838 CeedOperator_Ref *impl;
839
840 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done));
841 if (is_setup_done) return CEED_ERROR_SUCCESS;
842
843 CeedCallBackend(CeedOperatorGetData(op, &impl));
844 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
845 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q));
846 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf));
847 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
848 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
849
850 // Allocate
851 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full));
852
853 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_in));
854 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->skip_rstr_out));
855 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->apply_add_basis_out));
856 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states));
857 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in));
858 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out));
859 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in));
860 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out));
861
862 impl->num_inputs = num_input_fields;
863 impl->num_outputs = num_output_fields;
864
865 // Set up infield and outfield pointer arrays
866 // Infields
867 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->skip_rstr_in, NULL, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0,
868 num_input_fields, Q));
869 // Outfields
870 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->skip_rstr_out, impl->apply_add_basis_out, impl->e_vecs_full,
871 impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q));
872
873 // Identity QFunctions
874 if (impl->is_identity_qf) {
875 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0]));
876 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->e_vecs_out[0]));
877 }
878
879 CeedCallBackend(CeedOperatorSetSetupDone(op));
880 CeedCallBackend(CeedQFunctionDestroy(&qf));
881 return CEED_ERROR_SUCCESS;
882 }
883
884 //------------------------------------------------------------------------------
885 // Input Basis Action
886 //------------------------------------------------------------------------------
CeedOperatorInputBasisAtPoints_Ref(CeedInt e,CeedInt num_points_offset,CeedInt num_points,CeedQFunctionField * qf_input_fields,CeedOperatorField * op_input_fields,CeedInt num_input_fields,CeedVector in_vec,CeedVector point_coords_elem,bool skip_active,bool skip_passive,CeedScalar * e_data[2* CEED_FIELD_MAX],CeedOperator_Ref * impl,CeedRequest * request)887 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields,
888 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec,
889 CeedVector point_coords_elem, bool skip_active, bool skip_passive,
890 CeedScalar *e_data[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl, CeedRequest *request) {
891 for (CeedInt i = 0; i < num_input_fields; i++) {
892 bool is_active;
893 CeedInt elem_size, size, num_comp;
894 CeedRestrictionType rstr_type;
895 CeedEvalMode eval_mode;
896 CeedVector vec;
897 CeedElemRestriction elem_rstr;
898 CeedBasis basis;
899
900 // Skip active input
901 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
902 is_active = vec == CEED_VECTOR_ACTIVE;
903 CeedCallBackend(CeedVectorDestroy(&vec));
904 if (skip_active && is_active) continue;
905 if (skip_passive && !is_active) continue;
906
907 // Get elem_size, eval_mode, size
908 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
909 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
910 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode));
911 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size));
912 // Restrict block active input
913 // When skipping passive inputs, we're doing assembly and should not restrict
914 if (is_active && !impl->skip_rstr_in[i] && !skip_passive) {
915 if (rstr_type == CEED_RESTRICTION_POINTS) {
916 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request));
917 } else {
918 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request));
919 }
920 }
921 // Basis action
922 switch (eval_mode) {
923 case CEED_EVAL_NONE:
924 if (!is_active) {
925 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size]));
926 }
927 break;
928 // Note - these basis eval modes require FEM fields
929 case CEED_EVAL_INTERP:
930 case CEED_EVAL_GRAD:
931 case CEED_EVAL_DIV:
932 case CEED_EVAL_CURL:
933 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis));
934 if (!is_active) {
935 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
936 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
937 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][(CeedSize)e * elem_size * num_comp]));
938 }
939 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i],
940 impl->q_vecs_in[i]));
941 CeedCallBackend(CeedBasisDestroy(&basis));
942 break;
943 case CEED_EVAL_WEIGHT:
944 break; // No action
945 }
946 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
947 }
948 return CEED_ERROR_SUCCESS;
949 }
950
951 //------------------------------------------------------------------------------
952 // Output Basis Action
953 //------------------------------------------------------------------------------
CeedOperatorOutputBasisAtPoints_Ref(CeedInt e,CeedInt num_points_offset,CeedInt num_points,CeedQFunctionField * qf_output_fields,CeedOperatorField * op_output_fields,CeedInt num_input_fields,CeedInt num_output_fields,bool * apply_add_basis,bool * skip_rstr,CeedOperator op,CeedVector out_vec,CeedVector point_coords_elem,bool skip_passive,CeedOperator_Ref * impl,CeedRequest * request)954 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields,
955 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields,
956 bool *apply_add_basis, bool *skip_rstr, CeedOperator op, CeedVector out_vec,
957 CeedVector point_coords_elem, bool skip_passive, CeedOperator_Ref *impl, CeedRequest *request) {
958 for (CeedInt i = 0; i < num_output_fields; i++) {
959 bool is_active;
960 CeedRestrictionType rstr_type;
961 CeedEvalMode eval_mode;
962 CeedVector vec;
963 CeedElemRestriction elem_rstr;
964 CeedBasis basis;
965
966 // Skip active input
967 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
968 is_active = vec == CEED_VECTOR_ACTIVE;
969 CeedCallBackend(CeedVectorDestroy(&vec));
970 if (skip_passive && !is_active) continue;
971
972 // Get elem_size, eval_mode, size
973 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
974 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode));
975 // Basis action
976 switch (eval_mode) {
977 case CEED_EVAL_NONE:
978 break; // No action
979 case CEED_EVAL_INTERP:
980 case CEED_EVAL_GRAD:
981 case CEED_EVAL_DIV:
982 case CEED_EVAL_CURL:
983 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis));
984 if (apply_add_basis[i]) {
985 CeedCallBackend(CeedBasisApplyAddAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i],
986 impl->e_vecs_out[i]));
987 } else {
988 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i],
989 impl->e_vecs_out[i]));
990 }
991 CeedCallBackend(CeedBasisDestroy(&basis));
992 break;
993 // LCOV_EXCL_START
994 case CEED_EVAL_WEIGHT: {
995 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode");
996 // LCOV_EXCL_STOP
997 }
998 }
999 // Restrict output block
1000 // When skipping passive outputs, we're doing assembly and should not restrict
1001 if (skip_rstr[i] || skip_passive) {
1002 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1003 continue;
1004 }
1005
1006 // Get output vector
1007 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1008 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
1009 if (is_active) vec = out_vec;
1010 // Restrict
1011 if (rstr_type == CEED_RESTRICTION_POINTS) {
1012 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request));
1013 } else {
1014 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request));
1015 }
1016 if (!is_active) CeedCallBackend(CeedVectorDestroy(&vec));
1017 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1018 }
1019 return CEED_ERROR_SUCCESS;
1020 }
1021
1022 //------------------------------------------------------------------------------
1023 // Operator Apply
1024 //------------------------------------------------------------------------------
CeedOperatorApplyAddAtPoints_Ref(CeedOperator op,CeedVector in_vec,CeedVector out_vec,CeedRequest * request)1025 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) {
1026 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem;
1027 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0};
1028 CeedVector point_coords = NULL;
1029 CeedElemRestriction rstr_points = NULL;
1030 CeedQFunctionField *qf_input_fields, *qf_output_fields;
1031 CeedQFunction qf;
1032 CeedOperatorField *op_input_fields, *op_output_fields;
1033 CeedOperator_Ref *impl;
1034
1035 CeedCallBackend(CeedOperatorGetData(op, &impl));
1036 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
1037 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1038 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1039 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1040
1041 // Setup
1042 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op));
1043
1044 // Point coordinates
1045 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords));
1046
1047 // Input Evecs and Restriction
1048 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request));
1049
1050 // Loop through elements
1051 for (CeedInt e = 0; e < num_elem; e++) {
1052 CeedInt num_points;
1053
1054 // Setup points for element
1055 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request));
1056 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points));
1057 if (num_points < 1) continue;
1058
1059 // Input basis apply
1060 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec,
1061 impl->point_coords_elem, false, false, e_data, impl, request));
1062
1063 // Q function
1064 if (!impl->is_identity_qf) {
1065 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out));
1066 }
1067
1068 // Output basis apply and restriction
1069 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields,
1070 num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec,
1071 impl->point_coords_elem, false, impl, request));
1072
1073 num_points_offset += num_points;
1074 }
1075
1076 // Restore input arrays
1077 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl));
1078
1079 // Cleanup point coordinates
1080 CeedCallBackend(CeedVectorDestroy(&point_coords));
1081 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1082 CeedCallBackend(CeedQFunctionDestroy(&qf));
1083 return CEED_ERROR_SUCCESS;
1084 }
1085
1086 //------------------------------------------------------------------------------
1087 // Core code for assembling linear QFunction
1088 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(CeedOperator op,bool build_objects,CeedVector * assembled,CeedElemRestriction * rstr,CeedRequest * request)1089 static inline int CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled,
1090 CeedElemRestriction *rstr, CeedRequest *request) {
1091 Ceed ceed;
1092 CeedInt qf_size_in, qf_size_out, max_num_points, num_elem, num_input_fields, num_output_fields, num_points_offset = 0;
1093 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL};
1094 CeedVector point_coords = NULL;
1095 CeedQFunctionField *qf_input_fields, *qf_output_fields;
1096 CeedQFunction qf;
1097 CeedOperatorField *op_input_fields, *op_output_fields;
1098 CeedOperator_Ref *impl;
1099 CeedElemRestriction rstr_points = NULL;
1100
1101 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1102 CeedCallBackend(CeedOperatorGetData(op, &impl));
1103 qf_size_in = impl->qf_size_in;
1104 qf_size_out = impl->qf_size_out;
1105 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1106 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
1107 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1108 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1109
1110 // Setup
1111 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op));
1112
1113 // Check for restriction only operator
1114 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported");
1115
1116 // Point coordinates
1117 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords));
1118 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points));
1119
1120 // Input Evecs and Restriction
1121 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request));
1122
1123 // Count number of active input fields
1124 if (qf_size_in == 0) {
1125 for (CeedInt i = 0; i < num_input_fields; i++) {
1126 CeedInt field_size;
1127 CeedVector vec;
1128
1129 // Get input vector
1130 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1131 // Check if active input
1132 if (vec == CEED_VECTOR_ACTIVE) {
1133 // Check that all active inputs are nodal fields
1134 {
1135 CeedElemRestriction elem_rstr;
1136 bool is_at_points = false;
1137
1138 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
1139 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points));
1140 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1141 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points");
1142 }
1143 // Get size of active input
1144 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size));
1145 qf_size_in += field_size;
1146 }
1147 CeedCallBackend(CeedVectorDestroy(&vec));
1148 }
1149 CeedCheck(qf_size_in, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
1150 impl->qf_size_in = qf_size_in;
1151 }
1152
1153 // Count number of active output fields
1154 if (qf_size_out == 0) {
1155 for (CeedInt i = 0; i < num_output_fields; i++) {
1156 CeedInt field_size;
1157 CeedVector vec;
1158
1159 // Get output vector
1160 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec));
1161 // Check if active output
1162 if (vec == CEED_VECTOR_ACTIVE) {
1163 // Check that all active inputs are nodal fields
1164 {
1165 CeedElemRestriction elem_rstr;
1166 bool is_at_points = false;
1167
1168 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr));
1169 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points));
1170 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1171 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points");
1172 }
1173 // Get size of active output
1174 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size));
1175 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0));
1176 qf_size_out += field_size;
1177 }
1178 CeedCallBackend(CeedVectorDestroy(&vec));
1179 }
1180 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs");
1181 impl->qf_size_out = qf_size_out;
1182 }
1183
1184 // Build objects if needed
1185 if (build_objects) {
1186 CeedInt num_points_total;
1187 const CeedInt *offsets;
1188
1189 CeedCallBackend(CeedElemRestrictionGetNumPoints(rstr_points, &num_points_total));
1190
1191 // Create output restriction (at points)
1192 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr_points, CEED_MEM_HOST, &offsets));
1193 CeedCallBackend(CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points_total, qf_size_in * qf_size_out,
1194 qf_size_in * qf_size_out * num_points_total, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, rstr));
1195 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr_points, &offsets));
1196
1197 // Create assembled vector
1198 CeedCallBackend(CeedElemRestrictionCreateVector(*rstr, assembled, NULL));
1199 }
1200 // Clear output vector
1201 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0));
1202 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array));
1203
1204 // Loop through elements
1205 for (CeedInt e = 0; e < num_elem; e++) {
1206 CeedInt num_points;
1207
1208 // Setup points for element
1209 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request));
1210 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points));
1211 if (num_points < 1) continue;
1212
1213 // Input basis apply
1214 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, NULL,
1215 impl->point_coords_elem, true, false, e_data_full, impl, request));
1216
1217 // Assemble QFunction
1218 for (CeedInt i = 0; i < num_input_fields; i++) {
1219 bool is_active;
1220 CeedInt field_size;
1221 CeedVector vec;
1222
1223 // Get input vector
1224 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1225 is_active = vec == CEED_VECTOR_ACTIVE;
1226 CeedCallBackend(CeedVectorDestroy(&vec));
1227 // Check if active input
1228 if (!is_active) continue;
1229 // Get size of active input
1230 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size));
1231 for (CeedInt field = 0; field < field_size; field++) {
1232 // Set current portion of input to 1.0
1233 {
1234 CeedScalar *array;
1235
1236 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array));
1237 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 1.0;
1238 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array));
1239 }
1240
1241 if (!impl->is_identity_qf) {
1242 // Set Outputs
1243 for (CeedInt out = 0; out < num_output_fields; out++) {
1244 CeedVector vec;
1245 CeedInt field_size;
1246
1247 // Get output vector
1248 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
1249 // Check if active output
1250 if (vec == CEED_VECTOR_ACTIVE) {
1251 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array));
1252 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size));
1253 assembled_array += field_size * num_points; // Advance the pointer by the size of the output
1254 }
1255 CeedCallBackend(CeedVectorDestroy(&vec));
1256 }
1257 // Apply QFunction
1258 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out));
1259 } else {
1260 const CeedScalar *array;
1261 CeedInt field_size;
1262
1263 // Copy Identity Outputs
1264 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size));
1265 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array));
1266 for (CeedInt j = 0; j < field_size * num_points; j++) assembled_array[j] = array[j];
1267 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array));
1268 assembled_array += field_size * num_points;
1269 }
1270 // Reset input to 0.0
1271 {
1272 CeedScalar *array;
1273
1274 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array));
1275 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 0.0;
1276 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array));
1277 }
1278 }
1279 }
1280 num_points_offset += num_points;
1281 }
1282
1283 // Un-set output Qvecs to prevent accidental overwrite of Assembled
1284 if (!impl->is_identity_qf) {
1285 for (CeedInt out = 0; out < num_output_fields; out++) {
1286 CeedVector vec;
1287
1288 // Get output vector
1289 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec));
1290 // Check if active output
1291 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) {
1292 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL));
1293 }
1294 CeedCallBackend(CeedVectorDestroy(&vec));
1295 }
1296 }
1297
1298 // Restore input arrays
1299 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl));
1300
1301 // Restore output
1302 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array));
1303
1304 // Cleanup
1305 CeedCallBackend(CeedDestroy(&ceed));
1306 CeedCallBackend(CeedVectorDestroy(&point_coords));
1307 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1308 CeedCallBackend(CeedQFunctionDestroy(&qf));
1309 return CEED_ERROR_SUCCESS;
1310 }
1311
1312 //------------------------------------------------------------------------------
1313 // Assemble Linear QFunction
1314 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionAtPoints_Ref(CeedOperator op,CeedVector * assembled,CeedElemRestriction * rstr,CeedRequest * request)1315 static int CeedOperatorLinearAssembleQFunctionAtPoints_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) {
1316 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, true, assembled, rstr, request);
1317 }
1318
1319 //------------------------------------------------------------------------------
1320 // Update Assembled Linear QFunction
1321 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op,CeedVector assembled,CeedElemRestriction rstr,CeedRequest * request)1322 static int CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr,
1323 CeedRequest *request) {
1324 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, false, &assembled, &rstr, request);
1325 }
1326
1327 //------------------------------------------------------------------------------
1328 // Assemble Operator Diagonal AtPoints
1329 //------------------------------------------------------------------------------
CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op,CeedVector assembled,CeedRequest * request)1330 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, CeedVector assembled, CeedRequest *request) {
1331 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1;
1332 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0};
1333 Ceed ceed;
1334 CeedVector point_coords = NULL, in_vec, out_vec;
1335 CeedElemRestriction rstr_points = NULL;
1336 CeedQFunctionField *qf_input_fields, *qf_output_fields;
1337 CeedQFunction qf;
1338 CeedOperatorField *op_input_fields, *op_output_fields;
1339 CeedOperator_Ref *impl;
1340
1341 CeedCallBackend(CeedOperatorGetData(op, &impl));
1342 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
1343 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1344 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1345 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1346
1347 // Setup
1348 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op));
1349
1350 // Ceed
1351 {
1352 Ceed ceed_parent;
1353
1354 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1355 CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
1356 CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed));
1357 CeedCallBackend(CeedDestroy(&ceed_parent));
1358 }
1359
1360 // Point coordinates
1361 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords));
1362
1363 // Input and output vectors
1364 {
1365 CeedSize input_size, output_size;
1366
1367 CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
1368 CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec));
1369 CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec));
1370 CeedCallBackend(CeedVectorSetValue(out_vec, 0.0));
1371 }
1372
1373 // Clear input Evecs
1374 for (CeedInt i = 0; i < num_input_fields; i++) {
1375 bool is_active;
1376 CeedVector vec;
1377
1378 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1379 is_active = vec == CEED_VECTOR_ACTIVE;
1380 CeedCallBackend(CeedVectorDestroy(&vec));
1381 if (!is_active || impl->skip_rstr_in[i]) continue;
1382 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0));
1383 }
1384
1385 // Input Evecs and Restriction
1386 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request));
1387
1388 // Loop through elements
1389 for (CeedInt e = 0; e < num_elem; e++) {
1390 CeedInt num_points, e_vec_size = 0;
1391
1392 // Setup points for element
1393 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request));
1394 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points));
1395 if (num_points < 1) continue;
1396
1397 // Input basis apply for non-active bases
1398 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec,
1399 impl->point_coords_elem, true, false, e_data, impl, request));
1400
1401 // Loop over points on element
1402 for (CeedInt i = 0; i < num_input_fields; i++) {
1403 bool is_active_at_points = true, is_active;
1404 CeedInt elem_size_active = 1;
1405 CeedRestrictionType rstr_type;
1406 CeedVector vec;
1407 CeedElemRestriction elem_rstr;
1408
1409 // -- Skip non-active input
1410 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1411 is_active = vec == CEED_VECTOR_ACTIVE;
1412 CeedCallBackend(CeedVectorDestroy(&vec));
1413 if (!is_active || impl->skip_rstr_in[i]) continue;
1414
1415 // -- Get active restriction type
1416 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
1417 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1418 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS;
1419 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active));
1420 else elem_size_active = num_points;
1421 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active));
1422 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1423
1424 e_vec_size = elem_size_active * num_comp_active;
1425 for (CeedInt s = 0; s < e_vec_size; s++) {
1426 // -- Update unit vector
1427 {
1428 CeedScalar *array;
1429
1430 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array));
1431 array[s] = 1.0;
1432 if (s > 0) array[s - 1] = 0.0;
1433 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
1434 }
1435 // Input basis apply for active bases
1436 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields,
1437 in_vec, impl->point_coords_elem, false, true, e_data, impl, request));
1438
1439 // -- Q function
1440 if (!impl->is_identity_qf) {
1441 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out));
1442 }
1443
1444 // -- Output basis apply and restriction
1445 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields,
1446 num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec,
1447 impl->point_coords_elem, true, impl, request));
1448
1449 // -- Grab diagonal value
1450 for (CeedInt j = 0; j < num_output_fields; j++) {
1451 bool is_active;
1452 CeedInt elem_size = 0;
1453 CeedRestrictionType rstr_type;
1454 CeedVector vec;
1455 CeedElemRestriction elem_rstr;
1456
1457 // ---- Skip non-active output
1458 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec));
1459 is_active = vec == CEED_VECTOR_ACTIVE;
1460 CeedCallBackend(CeedVectorDestroy(&vec));
1461 if (!is_active || impl->skip_rstr_out[j]) continue;
1462
1463 // ---- Check if elem size matches
1464 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr));
1465 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1466 if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) {
1467 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1468 continue;
1469 }
1470 if (rstr_type == CEED_RESTRICTION_POINTS) {
1471 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, e, &elem_size));
1472 } else {
1473 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
1474 }
1475 {
1476 CeedInt num_comp = 0;
1477
1478 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1479 if (e_vec_size != num_comp * elem_size) {
1480 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1481 continue;
1482 }
1483 }
1484 // ---- Update output vector
1485 {
1486 CeedScalar *array, current_value = 0.0;
1487
1488 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &array));
1489 current_value = array[s];
1490 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &array));
1491 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_out[j], 0.0));
1492 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[j], CEED_MEM_HOST, &array));
1493 array[s] = current_value;
1494 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[j], &array));
1495 }
1496 // ---- Restrict output block
1497 if (rstr_type == CEED_RESTRICTION_POINTS) {
1498 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[j], assembled, request));
1499 } else {
1500 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[j], assembled, request));
1501 }
1502 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1503 }
1504 // -- Reset unit vector
1505 if (s == e_vec_size - 1) {
1506 CeedScalar *array;
1507
1508 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array));
1509 array[s] = 0.0;
1510 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
1511 }
1512 }
1513 }
1514 num_points_offset += num_points;
1515 }
1516
1517 // Restore input arrays
1518 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl));
1519
1520 // Cleanup
1521 CeedCallBackend(CeedDestroy(&ceed));
1522 CeedCallBackend(CeedVectorDestroy(&in_vec));
1523 CeedCallBackend(CeedVectorDestroy(&out_vec));
1524 CeedCallBackend(CeedVectorDestroy(&point_coords));
1525 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1526 CeedCallBackend(CeedQFunctionDestroy(&qf));
1527 return CEED_ERROR_SUCCESS;
1528 }
1529
1530 //------------------------------------------------------------------------------
1531 // Assemble Operator AtPoints
1532 //------------------------------------------------------------------------------
CeedOperatorAssembleSingleAtPoints_Ref(CeedOperator op,CeedInt offset,CeedVector values)1533 static int CeedOperatorAssembleSingleAtPoints_Ref(CeedOperator op, CeedInt offset, CeedVector values) {
1534 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, num_comp_active = 1;
1535 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}, *assembled;
1536 Ceed ceed;
1537 CeedVector point_coords = NULL, in_vec, out_vec;
1538 CeedElemRestriction rstr_points = NULL;
1539 CeedQFunctionField *qf_input_fields, *qf_output_fields;
1540 CeedQFunction qf;
1541 CeedOperatorField *op_input_fields, *op_output_fields;
1542 CeedOperator_Ref *impl;
1543
1544 CeedCallBackend(CeedOperatorGetData(op, &impl));
1545 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem));
1546 CeedCallBackend(CeedOperatorGetQFunction(op, &qf));
1547 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields));
1548 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields));
1549
1550 // Setup
1551 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op));
1552
1553 // Ceed
1554 {
1555 Ceed ceed_parent;
1556
1557 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1558 CeedCallBackend(CeedGetParent(ceed, &ceed_parent));
1559 CeedCallBackend(CeedReferenceCopy(ceed_parent, &ceed));
1560 CeedCallBackend(CeedDestroy(&ceed_parent));
1561 }
1562
1563 // Point coordinates
1564 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords));
1565
1566 // Input and output vectors
1567 {
1568 CeedSize input_size, output_size;
1569
1570 CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size));
1571 CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec));
1572 CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec));
1573 CeedCallBackend(CeedVectorSetValue(out_vec, 0.0));
1574 }
1575
1576 // Assembled array
1577 CeedCallBackend(CeedVectorGetArray(values, CEED_MEM_HOST, &assembled));
1578
1579 // Clear input Evecs
1580 for (CeedInt i = 0; i < num_input_fields; i++) {
1581 bool is_active;
1582 CeedVector vec;
1583
1584 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1585 is_active = vec == CEED_VECTOR_ACTIVE;
1586 CeedCallBackend(CeedVectorDestroy(&vec));
1587 if (!is_active || impl->skip_rstr_in[i]) continue;
1588 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0));
1589 }
1590
1591 // Input Evecs and Restriction
1592 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, CEED_REQUEST_IMMEDIATE));
1593
1594 // Loop through elements
1595 for (CeedInt e = 0; e < num_elem; e++) {
1596 CeedInt num_points, e_vec_size = 0;
1597
1598 // Setup points for element
1599 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem,
1600 CEED_REQUEST_IMMEDIATE));
1601 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points));
1602 if (num_points < 1) continue;
1603
1604 // Input basis apply for non-active bases
1605 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec,
1606 impl->point_coords_elem, true, false, e_data, impl, CEED_REQUEST_IMMEDIATE));
1607
1608 // Loop over points on element
1609 for (CeedInt i = 0; i < num_input_fields; i++) {
1610 bool is_active_at_points = true, is_active;
1611 CeedInt elem_size_active = 1;
1612 CeedRestrictionType rstr_type;
1613 CeedVector vec;
1614 CeedElemRestriction elem_rstr;
1615
1616 // -- Skip non-active input
1617 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec));
1618 is_active = vec == CEED_VECTOR_ACTIVE;
1619 CeedCallBackend(CeedVectorDestroy(&vec));
1620 if (!is_active || impl->skip_rstr_in[i]) continue;
1621
1622 // -- Get active restriction type
1623 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr));
1624 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1625 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS;
1626 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active));
1627 else elem_size_active = num_points;
1628 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active));
1629 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1630
1631 e_vec_size = elem_size_active * num_comp_active;
1632 for (CeedInt s = 0; s < e_vec_size; s++) {
1633 const CeedInt comp_in = s / elem_size_active;
1634 const CeedInt node_in = s % elem_size_active;
1635
1636 // -- Update unit vector
1637 {
1638 CeedScalar *array;
1639
1640 if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0));
1641 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array));
1642 array[s] = 1.0;
1643 if (s > 0) array[s - 1] = 0.0;
1644 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
1645 }
1646 // Input basis apply for active bases
1647 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields,
1648 in_vec, impl->point_coords_elem, false, true, e_data, impl, CEED_REQUEST_IMMEDIATE));
1649
1650 // -- Q function
1651 if (!impl->is_identity_qf) {
1652 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out));
1653 }
1654
1655 // -- Output basis apply and restriction
1656 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields,
1657 num_output_fields, impl->apply_add_basis_out, impl->skip_rstr_out, op, out_vec,
1658 impl->point_coords_elem, true, impl, CEED_REQUEST_IMMEDIATE));
1659
1660 // -- Build element matrix
1661 for (CeedInt j = 0; j < num_output_fields; j++) {
1662 bool is_active;
1663 CeedInt elem_size = 0;
1664 CeedRestrictionType rstr_type;
1665 CeedVector vec;
1666 CeedElemRestriction elem_rstr;
1667
1668 // ---- Skip non-active output
1669 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[j], &vec));
1670 is_active = vec == CEED_VECTOR_ACTIVE;
1671 CeedCallBackend(CeedVectorDestroy(&vec));
1672 if (!is_active || impl->skip_rstr_out[j]) continue;
1673
1674 // ---- Check if elem size matches
1675 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[j], &elem_rstr));
1676 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type));
1677 if (is_active_at_points && rstr_type != CEED_RESTRICTION_POINTS) {
1678 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1679 continue;
1680 }
1681 if (rstr_type == CEED_RESTRICTION_POINTS) {
1682 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(elem_rstr, e, &elem_size));
1683 } else {
1684 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size));
1685 }
1686 {
1687 CeedInt num_comp = 0;
1688
1689 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp));
1690 CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr));
1691 if (e_vec_size != num_comp * elem_size) continue;
1692 }
1693 // ---- Copy output
1694 {
1695 const CeedScalar *output;
1696
1697 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_out[j], CEED_MEM_HOST, &output));
1698 for (CeedInt k = 0; k < e_vec_size; k++) {
1699 const CeedInt comp_out = k / elem_size_active;
1700 const CeedInt node_out = k % elem_size_active;
1701
1702 assembled[offset + e * e_vec_size * e_vec_size + (comp_in * num_comp_active + comp_out) * elem_size_active * elem_size_active +
1703 node_out * elem_size_active + node_in] = output[k];
1704 }
1705 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_out[j], &output));
1706 }
1707 }
1708 // -- Reset unit vector
1709 if (s == e_vec_size - 1) {
1710 CeedScalar *array;
1711
1712 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array));
1713 array[s] = 0.0;
1714 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array));
1715 }
1716 }
1717 }
1718 num_points_offset += num_points;
1719 }
1720
1721 // Restore input arrays
1722 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl));
1723
1724 // Restore assembled values
1725 CeedCallBackend(CeedVectorRestoreArray(values, &assembled));
1726
1727 // Cleanup
1728 CeedCallBackend(CeedDestroy(&ceed));
1729 CeedCallBackend(CeedVectorDestroy(&in_vec));
1730 CeedCallBackend(CeedVectorDestroy(&out_vec));
1731 CeedCallBackend(CeedVectorDestroy(&point_coords));
1732 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points));
1733 CeedCallBackend(CeedQFunctionDestroy(&qf));
1734 return CEED_ERROR_SUCCESS;
1735 }
1736
1737 //------------------------------------------------------------------------------
1738 // Operator Destroy
1739 //------------------------------------------------------------------------------
CeedOperatorDestroy_Ref(CeedOperator op)1740 static int CeedOperatorDestroy_Ref(CeedOperator op) {
1741 CeedOperator_Ref *impl;
1742
1743 CeedCallBackend(CeedOperatorGetData(op, &impl));
1744 CeedCallBackend(CeedFree(&impl->skip_rstr_in));
1745 CeedCallBackend(CeedFree(&impl->skip_rstr_out));
1746 CeedCallBackend(CeedFree(&impl->e_data_out_indices));
1747 CeedCallBackend(CeedFree(&impl->apply_add_basis_out));
1748 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) {
1749 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i]));
1750 }
1751 CeedCallBackend(CeedFree(&impl->e_vecs_full));
1752 CeedCallBackend(CeedFree(&impl->input_states));
1753
1754 for (CeedInt i = 0; i < impl->num_inputs; i++) {
1755 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i]));
1756 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i]));
1757 }
1758 CeedCallBackend(CeedFree(&impl->e_vecs_in));
1759 CeedCallBackend(CeedFree(&impl->q_vecs_in));
1760
1761 for (CeedInt i = 0; i < impl->num_outputs; i++) {
1762 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i]));
1763 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i]));
1764 }
1765 CeedCallBackend(CeedFree(&impl->e_vecs_out));
1766 CeedCallBackend(CeedFree(&impl->q_vecs_out));
1767 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem));
1768
1769 CeedCallBackend(CeedFree(&impl));
1770 return CEED_ERROR_SUCCESS;
1771 }
1772
1773 //------------------------------------------------------------------------------
1774 // Operator Create
1775 //------------------------------------------------------------------------------
CeedOperatorCreate_Ref(CeedOperator op)1776 int CeedOperatorCreate_Ref(CeedOperator op) {
1777 Ceed ceed;
1778 CeedOperator_Ref *impl;
1779
1780 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1781 CeedCallBackend(CeedCalloc(1, &impl));
1782 CeedCallBackend(CeedOperatorSetData(op, impl));
1783 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref));
1784 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref));
1785 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref));
1786 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref));
1787 CeedCallBackend(CeedDestroy(&ceed));
1788 return CEED_ERROR_SUCCESS;
1789 }
1790
1791 //------------------------------------------------------------------------------
1792 // Operator Create At Points
1793 //------------------------------------------------------------------------------
CeedOperatorCreateAtPoints_Ref(CeedOperator op)1794 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) {
1795 Ceed ceed;
1796 CeedOperator_Ref *impl;
1797
1798 CeedCallBackend(CeedOperatorGetCeed(op, &ceed));
1799 CeedCallBackend(CeedCalloc(1, &impl));
1800 CeedCallBackend(CeedOperatorSetData(op, impl));
1801 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Ref));
1802 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate",
1803 CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref));
1804 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref));
1805 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleSingle", CeedOperatorAssembleSingleAtPoints_Ref));
1806 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref));
1807 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref));
1808 CeedCallBackend(CeedDestroy(&ceed));
1809 return CEED_ERROR_SUCCESS;
1810 }
1811
1812 //------------------------------------------------------------------------------
1813