1 // Copyright (c) 2017-2024, 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 //------------------------------------------------------------------------------ 19 static int CeedOperatorSetupFields_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs_full, CeedVector *e_vecs, 20 CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 21 Ceed ceed; 22 CeedSize e_size, q_size; 23 CeedInt num_comp, size, P; 24 CeedQFunctionField *qf_fields; 25 CeedOperatorField *op_fields; 26 27 { 28 Ceed ceed_parent; 29 30 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 31 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 32 if (ceed_parent) ceed = ceed_parent; 33 } 34 if (is_input) { 35 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 36 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 37 } else { 38 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 39 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 40 } 41 42 // Loop over fields 43 for (CeedInt i = 0; i < num_fields; i++) { 44 CeedEvalMode eval_mode; 45 CeedElemRestriction elem_rstr; 46 CeedBasis basis; 47 48 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 49 if (eval_mode != CEED_EVAL_WEIGHT) { 50 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 51 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 52 } 53 54 switch (eval_mode) { 55 case CEED_EVAL_NONE: 56 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 57 q_size = (CeedSize)Q * size; 58 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 59 break; 60 case CEED_EVAL_INTERP: 61 case CEED_EVAL_GRAD: 62 case CEED_EVAL_DIV: 63 case CEED_EVAL_CURL: 64 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 65 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 66 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 67 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 68 e_size = (CeedSize)P * num_comp; 69 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 70 q_size = (CeedSize)Q * size; 71 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 72 break; 73 case CEED_EVAL_WEIGHT: // Only on input fields 74 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 75 q_size = (CeedSize)Q; 76 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 77 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_vecs[i])); 78 break; 79 } 80 } 81 return CEED_ERROR_SUCCESS; 82 } 83 84 //------------------------------------------------------------------------------ 85 // Setup Operator 86 //------------------------------------------------------------------------------/* 87 static int CeedOperatorSetup_Ref(CeedOperator op) { 88 bool is_setup_done; 89 CeedInt Q, num_input_fields, num_output_fields; 90 CeedQFunctionField *qf_input_fields, *qf_output_fields; 91 CeedQFunction qf; 92 CeedOperatorField *op_input_fields, *op_output_fields; 93 CeedOperator_Ref *impl; 94 95 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 96 if (is_setup_done) return CEED_ERROR_SUCCESS; 97 98 CeedCallBackend(CeedOperatorGetData(op, &impl)); 99 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 100 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 101 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 102 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 103 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 104 105 // Allocate 106 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 107 108 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 109 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 110 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 111 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 112 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 113 114 impl->num_inputs = num_input_fields; 115 impl->num_outputs = num_output_fields; 116 117 // Set up infield and outfield e_vecs and q_vecs 118 // Infields 119 CeedCallBackend(CeedOperatorSetupFields_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 120 // Outfields 121 CeedCallBackend( 122 CeedOperatorSetupFields_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, num_output_fields, Q)); 123 124 // Identity QFunctions 125 if (impl->is_identity_qf) { 126 CeedEvalMode in_mode, out_mode; 127 CeedQFunctionField *in_fields, *out_fields; 128 129 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); 130 CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); 131 CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); 132 133 if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { 134 impl->is_identity_rstr_op = true; 135 } else { 136 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 137 } 138 } 139 140 CeedCallBackend(CeedOperatorSetSetupDone(op)); 141 return CEED_ERROR_SUCCESS; 142 } 143 144 //------------------------------------------------------------------------------ 145 // Setup Operator Inputs 146 //------------------------------------------------------------------------------ 147 static inline int CeedOperatorSetupInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 148 CeedVector in_vec, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 149 CeedOperator_Ref *impl, CeedRequest *request) { 150 for (CeedInt i = 0; i < num_input_fields; i++) { 151 uint64_t state; 152 CeedEvalMode eval_mode; 153 CeedVector vec; 154 CeedElemRestriction elem_rstr; 155 156 // Get input vector 157 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 158 if (vec == CEED_VECTOR_ACTIVE) { 159 if (skip_active) continue; 160 else vec = in_vec; 161 } 162 163 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 164 // Restrict and Evec 165 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 166 } else { 167 // Restrict 168 CeedCallBackend(CeedVectorGetState(vec, &state)); 169 // Skip restriction if input is unchanged 170 if (state != impl->input_states[i] || vec == in_vec) { 171 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 172 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); 173 impl->input_states[i] = state; 174 } 175 // Get evec 176 CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data_full[i])); 177 } 178 } 179 return CEED_ERROR_SUCCESS; 180 } 181 182 //------------------------------------------------------------------------------ 183 // Input Basis Action 184 //------------------------------------------------------------------------------ 185 static inline int CeedOperatorInputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 186 CeedInt num_input_fields, const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], 187 CeedOperator_Ref *impl) { 188 for (CeedInt i = 0; i < num_input_fields; i++) { 189 CeedInt elem_size, size, num_comp; 190 CeedEvalMode eval_mode; 191 CeedElemRestriction elem_rstr; 192 CeedBasis basis; 193 194 // Skip active input 195 if (skip_active) { 196 CeedVector vec; 197 198 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 199 if (vec == CEED_VECTOR_ACTIVE) continue; 200 } 201 // Get elem_size, eval_mode, size 202 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 203 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 204 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 205 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 206 // Basis action 207 switch (eval_mode) { 208 case CEED_EVAL_NONE: 209 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * Q * size])); 210 break; 211 case CEED_EVAL_INTERP: 212 case CEED_EVAL_GRAD: 213 case CEED_EVAL_DIV: 214 case CEED_EVAL_CURL: 215 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 216 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 217 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i][(CeedSize)e * elem_size * num_comp])); 218 CeedCallBackend(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, eval_mode, impl->e_vecs_in[i], impl->q_vecs_in[i])); 219 break; 220 case CEED_EVAL_WEIGHT: 221 break; // No action 222 } 223 } 224 return CEED_ERROR_SUCCESS; 225 } 226 227 //------------------------------------------------------------------------------ 228 // Output Basis Action 229 //------------------------------------------------------------------------------ 230 static inline int CeedOperatorOutputBasis_Ref(CeedInt e, CeedInt Q, CeedQFunctionField *qf_output_fields, CeedOperatorField *op_output_fields, 231 CeedInt num_input_fields, CeedInt num_output_fields, CeedOperator op, 232 CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 233 for (CeedInt i = 0; i < num_output_fields; i++) { 234 CeedInt elem_size, num_comp; 235 CeedEvalMode eval_mode; 236 CeedElemRestriction elem_rstr; 237 CeedBasis basis; 238 239 // Get elem_size, eval_mode 240 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 241 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 242 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 243 // Basis action 244 switch (eval_mode) { 245 case CEED_EVAL_NONE: 246 break; // No action 247 case CEED_EVAL_INTERP: 248 case CEED_EVAL_GRAD: 249 case CEED_EVAL_DIV: 250 case CEED_EVAL_CURL: 251 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 252 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 253 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, 254 &e_data_full[i + num_input_fields][(CeedSize)e * elem_size * num_comp])); 255 CeedCallBackend(CeedBasisApply(basis, 1, CEED_TRANSPOSE, eval_mode, impl->q_vecs_out[i], impl->e_vecs_out[i])); 256 break; 257 // LCOV_EXCL_START 258 case CEED_EVAL_WEIGHT: { 259 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 260 // LCOV_EXCL_STOP 261 } 262 } 263 } 264 return CEED_ERROR_SUCCESS; 265 } 266 267 //------------------------------------------------------------------------------ 268 // Restore Input Vectors 269 //------------------------------------------------------------------------------ 270 static inline int CeedOperatorRestoreInputs_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, CeedOperatorField *op_input_fields, 271 const bool skip_active, CeedScalar *e_data_full[2 * CEED_FIELD_MAX], CeedOperator_Ref *impl) { 272 for (CeedInt i = 0; i < num_input_fields; i++) { 273 CeedEvalMode eval_mode; 274 275 // Skip active inputs 276 if (skip_active) { 277 CeedVector vec; 278 279 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 280 if (vec == CEED_VECTOR_ACTIVE) continue; 281 } 282 // Restore input 283 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 284 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 285 } else { 286 CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data_full[i])); 287 } 288 } 289 return CEED_ERROR_SUCCESS; 290 } 291 292 //------------------------------------------------------------------------------ 293 // Operator Apply 294 //------------------------------------------------------------------------------ 295 static int CeedOperatorApplyAdd_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 296 CeedInt Q, num_elem, num_input_fields, num_output_fields, size; 297 CeedEvalMode eval_mode; 298 CeedScalar *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 299 CeedQFunctionField *qf_input_fields, *qf_output_fields; 300 CeedQFunction qf; 301 CeedOperatorField *op_input_fields, *op_output_fields; 302 CeedOperator_Ref *impl; 303 304 CeedCallBackend(CeedOperatorGetData(op, &impl)); 305 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 306 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 307 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 308 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 309 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 310 311 // Setup 312 CeedCallBackend(CeedOperatorSetup_Ref(op)); 313 314 // Restriction only operator 315 if (impl->is_identity_rstr_op) { 316 CeedElemRestriction elem_rstr; 317 318 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[0], &elem_rstr)); 319 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_full[0], request)); 320 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[0], &elem_rstr)); 321 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[0], out_vec, request)); 322 return CEED_ERROR_SUCCESS; 323 } 324 325 // Input Evecs and Restriction 326 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, false, e_data_full, impl, request)); 327 328 // Output Evecs 329 for (CeedInt i = 0; i < num_output_fields; i++) { 330 CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_full[i + impl->num_inputs], CEED_MEM_HOST, &e_data_full[i + num_input_fields])); 331 } 332 333 // Loop through elements 334 for (CeedInt e = 0; e < num_elem; e++) { 335 // Output pointers 336 for (CeedInt i = 0; i < num_output_fields; i++) { 337 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 338 if (eval_mode == CEED_EVAL_NONE) { 339 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 340 CeedCallBackend( 341 CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data_full[i + num_input_fields][(CeedSize)e * Q * size])); 342 } 343 } 344 345 // Input basis apply 346 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, false, e_data_full, impl)); 347 348 // Q function 349 if (!impl->is_identity_qf) { 350 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 351 } 352 353 // Output basis apply 354 CeedCallBackend( 355 CeedOperatorOutputBasis_Ref(e, Q, qf_output_fields, op_output_fields, num_input_fields, num_output_fields, op, e_data_full, impl)); 356 } 357 358 // Output restriction 359 for (CeedInt i = 0; i < num_output_fields; i++) { 360 CeedVector vec; 361 CeedElemRestriction elem_rstr; 362 363 // Restore Evec 364 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_full[i + impl->num_inputs], &e_data_full[i + num_input_fields])); 365 // Get output vector 366 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 367 // Active 368 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 369 // Restrict 370 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 371 CeedCallBackend(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, impl->e_vecs_full[i + impl->num_inputs], vec, request)); 372 } 373 374 // Restore input arrays 375 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, false, e_data_full, impl)); 376 return CEED_ERROR_SUCCESS; 377 } 378 379 //------------------------------------------------------------------------------ 380 // Core code for assembling linear QFunction 381 //------------------------------------------------------------------------------ 382 static inline int CeedOperatorLinearAssembleQFunctionCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, CeedElemRestriction *rstr, 383 CeedRequest *request) { 384 Ceed ceed, ceed_parent; 385 CeedSize q_size; 386 CeedInt num_active_in, num_active_out, Q, num_elem, num_input_fields, num_output_fields, size; 387 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 388 CeedVector *active_in; 389 CeedQFunctionField *qf_input_fields, *qf_output_fields; 390 CeedQFunction qf; 391 CeedOperatorField *op_input_fields, *op_output_fields; 392 CeedOperator_Ref *impl; 393 394 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 395 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 396 CeedCallBackend(CeedOperatorGetData(op, &impl)); 397 active_in = impl->qf_active_in; 398 num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 399 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 400 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 401 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 402 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 403 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 404 405 // Setup 406 CeedCallBackend(CeedOperatorSetup_Ref(op)); 407 408 // Check for restriction only operator 409 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 410 411 // Input Evecs and Restriction 412 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 413 414 // Count number of active input fields 415 if (!num_active_in) { 416 for (CeedInt i = 0; i < num_input_fields; i++) { 417 CeedScalar *q_vec_array; 418 CeedVector vec; 419 420 // Get input vector 421 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 422 // Check if active input 423 if (vec == CEED_VECTOR_ACTIVE) { 424 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 425 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 426 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &q_vec_array)); 427 CeedCallBackend(CeedRealloc(num_active_in + size, &active_in)); 428 for (CeedInt field = 0; field < size; field++) { 429 q_size = (CeedSize)Q; 430 CeedCallBackend(CeedVectorCreate(ceed_parent, q_size, &active_in[num_active_in + field])); 431 CeedCallBackend(CeedVectorSetArray(active_in[num_active_in + field], CEED_MEM_HOST, CEED_USE_POINTER, &q_vec_array[field * Q])); 432 } 433 num_active_in += size; 434 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array)); 435 } 436 } 437 impl->num_active_in = num_active_in; 438 impl->qf_active_in = active_in; 439 } 440 441 // Count number of active output fields 442 if (!num_active_out) { 443 for (CeedInt i = 0; i < num_output_fields; i++) { 444 CeedVector vec; 445 446 // Get output vector 447 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 448 // Check if active output 449 if (vec == CEED_VECTOR_ACTIVE) { 450 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &size)); 451 num_active_out += size; 452 } 453 } 454 impl->num_active_out = num_active_out; 455 } 456 457 // Check sizes 458 CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 459 460 // Build objects if needed 461 if (build_objects) { 462 const CeedSize l_size = (CeedSize)num_elem * Q * num_active_in * num_active_out; 463 CeedInt strides[3] = {1, Q, num_active_in * num_active_out * Q}; /* *NOPAD* */ 464 465 // Create output restriction 466 CeedCallBackend(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, num_active_in * num_active_out, 467 num_active_in * num_active_out * num_elem * Q, strides, rstr)); 468 // Create assembled vector 469 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 470 } 471 // Clear output vector 472 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 473 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 474 475 // Loop through elements 476 for (CeedInt e = 0; e < num_elem; e++) { 477 // Input basis apply 478 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl)); 479 480 // Assemble QFunction 481 for (CeedInt in = 0; in < num_active_in; in++) { 482 // Set Inputs 483 CeedCallBackend(CeedVectorSetValue(active_in[in], 1.0)); 484 if (num_active_in > 1) { 485 CeedCallBackend(CeedVectorSetValue(active_in[(in + num_active_in - 1) % num_active_in], 0.0)); 486 } 487 if (!impl->is_identity_qf) { 488 // Set Outputs 489 for (CeedInt out = 0; out < num_output_fields; out++) { 490 CeedVector vec; 491 492 // Get output vector 493 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 494 // Check if active output 495 if (vec == CEED_VECTOR_ACTIVE) { 496 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 497 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &size)); 498 assembled_array += size * Q; // Advance the pointer by the size of the output 499 } 500 } 501 // Apply QFunction 502 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 503 } else { 504 const CeedScalar *q_vec_array; 505 506 // Copy Identity Outputs 507 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &size)); 508 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &q_vec_array)); 509 for (CeedInt i = 0; i < size * Q; i++) assembled_array[i] = q_vec_array[i]; 510 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &q_vec_array)); 511 assembled_array += size * Q; 512 } 513 } 514 } 515 516 // Un-set output Qvecs to prevent accidental overwrite of Assembled 517 if (!impl->is_identity_qf) { 518 for (CeedInt out = 0; out < num_output_fields; out++) { 519 CeedVector vec; 520 521 // Get output vector 522 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 523 // Check if active output 524 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 525 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 526 } 527 } 528 } 529 530 // Restore input arrays 531 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 532 533 // Restore output 534 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 535 return CEED_ERROR_SUCCESS; 536 } 537 538 //------------------------------------------------------------------------------ 539 // Assemble Linear QFunction 540 //------------------------------------------------------------------------------ 541 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 542 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request); 543 } 544 545 //------------------------------------------------------------------------------ 546 // Update Assembled Linear QFunction 547 //------------------------------------------------------------------------------ 548 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 549 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); 550 } 551 552 //------------------------------------------------------------------------------ 553 // Setup Input/Output Fields 554 //------------------------------------------------------------------------------ 555 static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs_full, CeedVector *e_vecs, 556 CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 557 Ceed ceed; 558 CeedSize e_size, q_size; 559 CeedInt max_num_points, num_comp, size, P; 560 CeedQFunctionField *qf_fields; 561 CeedOperatorField *op_fields; 562 563 { 564 Ceed ceed_parent; 565 566 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 567 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 568 if (ceed_parent) ceed = ceed_parent; 569 } 570 if (is_input) { 571 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 572 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 573 } else { 574 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 575 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 576 } 577 578 // Get max number of points 579 { 580 CeedInt dim; 581 CeedElemRestriction rstr_points = NULL; 582 CeedOperator_Ref *impl; 583 584 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 585 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 586 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); 587 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 588 CeedCallBackend(CeedOperatorGetData(op, &impl)); 589 if (is_input) { 590 CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); 591 CeedCallBackend(CeedVectorSetValue(impl->point_coords_elem, 0.0)); 592 } 593 } 594 595 // Loop over fields 596 for (CeedInt i = 0; i < num_fields; i++) { 597 CeedEvalMode eval_mode; 598 CeedBasis basis; 599 600 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 601 if (eval_mode != CEED_EVAL_WEIGHT) { 602 CeedElemRestriction elem_rstr; 603 604 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 605 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 606 CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); 607 } 608 609 switch (eval_mode) { 610 case CEED_EVAL_NONE: { 611 CeedVector vec; 612 613 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 614 e_size = (CeedSize)max_num_points * size; 615 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 616 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 617 if (vec == CEED_VECTOR_ACTIVE || !is_input) { 618 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i])); 619 } else { 620 q_size = (CeedSize)max_num_points * size; 621 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 622 } 623 break; 624 } 625 case CEED_EVAL_INTERP: 626 case CEED_EVAL_GRAD: 627 case CEED_EVAL_DIV: 628 case CEED_EVAL_CURL: 629 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 630 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 631 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 632 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 633 e_size = (CeedSize)P * num_comp; 634 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 635 q_size = (CeedSize)max_num_points * size; 636 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 637 break; 638 case CEED_EVAL_WEIGHT: // Only on input fields 639 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 640 q_size = (CeedSize)max_num_points; 641 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 642 CeedCallBackend( 643 CeedBasisApplyAtPoints(basis, max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 644 break; 645 } 646 // Initialize full arrays for E-vectors and Q-vectors 647 if (e_vecs[i]) CeedCallBackend(CeedVectorSetValue(e_vecs[i], 0.0)); 648 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0)); 649 } 650 return CEED_ERROR_SUCCESS; 651 } 652 653 //------------------------------------------------------------------------------ 654 // Setup Operator 655 //------------------------------------------------------------------------------ 656 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { 657 bool is_setup_done; 658 CeedInt Q, num_input_fields, num_output_fields; 659 CeedQFunctionField *qf_input_fields, *qf_output_fields; 660 CeedQFunction qf; 661 CeedOperatorField *op_input_fields, *op_output_fields; 662 CeedOperator_Ref *impl; 663 664 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 665 if (is_setup_done) return CEED_ERROR_SUCCESS; 666 667 CeedCallBackend(CeedOperatorGetData(op, &impl)); 668 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 669 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 670 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 671 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 672 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 673 674 // Allocate 675 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 676 677 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 678 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 679 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 680 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 681 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 682 683 impl->num_inputs = num_input_fields; 684 impl->num_outputs = num_output_fields; 685 686 // Set up infield and outfield pointer arrays 687 // Infields 688 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 689 // Outfields 690 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, 691 num_output_fields, Q)); 692 693 // Identity QFunctions 694 if (impl->is_identity_qf) { 695 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 696 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->e_vecs_out[0])); 697 } 698 699 CeedCallBackend(CeedOperatorSetSetupDone(op)); 700 return CEED_ERROR_SUCCESS; 701 } 702 703 //------------------------------------------------------------------------------ 704 // Input Basis Action 705 //------------------------------------------------------------------------------ 706 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, 707 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, 708 CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 709 CeedOperator_Ref *impl, CeedRequest *request) { 710 for (CeedInt i = 0; i < num_input_fields; i++) { 711 bool is_active_input = false; 712 CeedInt elem_size, size, num_comp; 713 CeedRestrictionType rstr_type; 714 CeedEvalMode eval_mode; 715 CeedVector vec; 716 CeedElemRestriction elem_rstr; 717 CeedBasis basis; 718 719 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 720 // Skip active input 721 is_active_input = vec == CEED_VECTOR_ACTIVE; 722 if (skip_active && is_active_input) continue; 723 724 // Get elem_size, eval_mode, size 725 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 726 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 727 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 728 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 729 // Restrict block active input 730 if (is_active_input) { 731 if (rstr_type == CEED_RESTRICTION_POINTS) { 732 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 733 } else { 734 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 735 } 736 } 737 // Basis action 738 switch (eval_mode) { 739 case CEED_EVAL_NONE: 740 if (!is_active_input) { 741 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); 742 } 743 break; 744 // Note - these basis eval modes require FEM fields 745 case CEED_EVAL_INTERP: 746 case CEED_EVAL_GRAD: 747 case CEED_EVAL_DIV: 748 case CEED_EVAL_CURL: 749 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 750 if (!is_active_input) { 751 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 752 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 753 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][(CeedSize)e * elem_size * num_comp])); 754 } 755 CeedCallBackend( 756 CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); 757 break; 758 case CEED_EVAL_WEIGHT: 759 break; // No action 760 } 761 } 762 return CEED_ERROR_SUCCESS; 763 } 764 765 //------------------------------------------------------------------------------ 766 // Output Basis Action 767 //------------------------------------------------------------------------------ 768 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, 769 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, 770 CeedOperator op, CeedVector out_vec, CeedVector point_coords_elem, CeedOperator_Ref *impl, 771 CeedRequest *request) { 772 for (CeedInt i = 0; i < num_output_fields; i++) { 773 CeedRestrictionType rstr_type; 774 CeedEvalMode eval_mode; 775 CeedVector vec; 776 CeedElemRestriction elem_rstr; 777 CeedBasis basis; 778 779 // Get elem_size, eval_mode, size 780 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 781 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 782 // Basis action 783 switch (eval_mode) { 784 case CEED_EVAL_NONE: 785 break; // No action 786 case CEED_EVAL_INTERP: 787 case CEED_EVAL_GRAD: 788 case CEED_EVAL_DIV: 789 case CEED_EVAL_CURL: 790 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 791 CeedCallBackend( 792 CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); 793 break; 794 // LCOV_EXCL_START 795 case CEED_EVAL_WEIGHT: { 796 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 797 // LCOV_EXCL_STOP 798 } 799 } 800 // Restrict output block 801 // Get output vector 802 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 803 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 804 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 805 // Restrict 806 if (rstr_type == CEED_RESTRICTION_POINTS) { 807 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 808 } else { 809 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 810 } 811 } 812 return CEED_ERROR_SUCCESS; 813 } 814 815 //------------------------------------------------------------------------------ 816 // Operator Apply 817 //------------------------------------------------------------------------------ 818 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 819 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; 820 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 821 CeedVector point_coords = NULL; 822 CeedElemRestriction rstr_points = NULL; 823 CeedQFunctionField *qf_input_fields, *qf_output_fields; 824 CeedQFunction qf; 825 CeedOperatorField *op_input_fields, *op_output_fields; 826 CeedOperator_Ref *impl; 827 828 CeedCallBackend(CeedOperatorGetData(op, &impl)); 829 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 830 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 831 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 832 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 833 834 // Setup 835 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 836 837 // Point coordinates 838 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 839 840 // Input Evecs and Restriction 841 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 842 843 // Loop through elements 844 for (CeedInt e = 0; e < num_elem; e++) { 845 CeedInt num_points; 846 847 // Setup points for element 848 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 849 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 850 851 // Input basis apply 852 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 853 impl->point_coords_elem, false, e_data, impl, request)); 854 855 // Q function 856 if (!impl->is_identity_qf) { 857 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 858 } 859 860 // Output basis apply and restriction 861 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 862 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 863 864 num_points_offset += num_points; 865 } 866 867 // Restore input arrays 868 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 869 870 // Cleanup point coordinates 871 CeedCallBackend(CeedVectorDestroy(&point_coords)); 872 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 873 return CEED_ERROR_SUCCESS; 874 } 875 876 //------------------------------------------------------------------------------ 877 // Core code for assembling linear QFunction 878 //------------------------------------------------------------------------------ 879 static inline int CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, 880 CeedElemRestriction *rstr, CeedRequest *request) { 881 Ceed ceed; 882 CeedSize q_size; 883 CeedInt num_active_in, num_active_out, max_num_points, num_elem, num_input_fields, num_output_fields, num_points_offset = 0; 884 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 885 CeedVector *active_in, point_coords = NULL; 886 CeedQFunctionField *qf_input_fields, *qf_output_fields; 887 CeedQFunction qf; 888 CeedOperatorField *op_input_fields, *op_output_fields; 889 CeedOperator_Ref *impl; 890 CeedElemRestriction rstr_points = NULL; 891 892 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 893 CeedCallBackend(CeedOperatorGetData(op, &impl)); 894 active_in = impl->qf_active_in; 895 num_active_in = impl->num_active_in, num_active_out = impl->num_active_out; 896 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 897 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 898 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 899 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 900 901 // Setup 902 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 903 904 // Check for restriction only operator 905 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 906 907 // Point coordinates 908 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 909 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 910 911 // Input Evecs and Restriction 912 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 913 914 // Count number of active input fields 915 if (!num_active_in) { 916 for (CeedInt i = 0; i < num_input_fields; i++) { 917 CeedScalar *q_vec_array; 918 CeedInt field_size; 919 CeedVector vec; 920 921 // Get input vector 922 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 923 // Check if active input 924 if (vec == CEED_VECTOR_ACTIVE) { 925 // Check that all active inputs are nodal fields 926 { 927 CeedElemRestriction elem_rstr; 928 bool is_at_points = false; 929 930 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 931 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 932 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 933 } 934 // Get size of active input 935 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 936 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 937 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &q_vec_array)); 938 CeedCallBackend(CeedRealloc(num_active_in + field_size, &active_in)); 939 for (CeedInt field = 0; field < field_size; field++) { 940 q_size = (CeedSize)max_num_points; 941 CeedCallBackend(CeedVectorCreate(ceed, q_size, &active_in[num_active_in + field])); 942 CeedCallBackend(CeedVectorSetArray(active_in[num_active_in + field], CEED_MEM_HOST, CEED_USE_POINTER, &q_vec_array[field * q_size])); 943 } 944 num_active_in += field_size; 945 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &q_vec_array)); 946 } 947 } 948 impl->num_active_in = num_active_in; 949 impl->qf_active_in = active_in; 950 } 951 952 // Count number of active output fields 953 if (!num_active_out) { 954 for (CeedInt i = 0; i < num_output_fields; i++) { 955 CeedVector vec; 956 CeedInt field_size; 957 958 // Get output vector 959 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 960 // Check if active output 961 if (vec == CEED_VECTOR_ACTIVE) { 962 // Check that all active inputs are nodal fields 963 { 964 CeedElemRestriction elem_rstr; 965 bool is_at_points = false; 966 967 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 968 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 969 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 970 } 971 // Get size of active output 972 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 973 num_active_out += field_size; 974 } 975 } 976 impl->num_active_out = num_active_out; 977 } 978 979 // Check sizes 980 CeedCheck(num_active_in > 0 && num_active_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 981 982 // Build objects if needed 983 if (build_objects) { 984 CeedInt num_points_total; 985 const CeedInt *offsets; 986 987 CeedCallBackend(CeedElemRestrictionGetNumPoints(rstr_points, &num_points_total)); 988 989 // Create output restriction (at points) 990 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr_points, CEED_MEM_HOST, &offsets)); 991 CeedCallBackend(CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points_total, num_active_in * num_active_out, 992 num_active_in * num_active_out * num_points_total, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, 993 rstr)); 994 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr_points, &offsets)); 995 996 // Create assembled vector 997 CeedCallBackend(CeedElemRestrictionCreateVector(*rstr, assembled, NULL)); 998 } 999 // Clear output vector 1000 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 1001 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 1002 1003 // Loop through elements 1004 for (CeedInt e = 0; e < num_elem; e++) { 1005 CeedInt num_points; 1006 1007 // Setup points for element 1008 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1009 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1010 1011 // Input basis apply 1012 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, NULL, 1013 impl->point_coords_elem, true, e_data_full, impl, request)); 1014 1015 // Assemble QFunction 1016 for (CeedInt in = 0; in < num_active_in; in++) { 1017 // Set Inputs 1018 CeedCallBackend(CeedVectorSetValue(active_in[in], 1.0)); 1019 if (num_active_in > 1) { 1020 CeedCallBackend(CeedVectorSetValue(active_in[(in + num_active_in - 1) % num_active_in], 0.0)); 1021 } 1022 if (!impl->is_identity_qf) { 1023 // Set Outputs 1024 for (CeedInt out = 0; out < num_output_fields; out++) { 1025 CeedVector vec; 1026 CeedInt field_size; 1027 1028 // Get output vector 1029 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1030 // Check if active output 1031 if (vec == CEED_VECTOR_ACTIVE) { 1032 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 1033 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 1034 assembled_array += field_size * num_points; // Advance the pointer by the size of the output 1035 } 1036 } 1037 // Apply QFunction 1038 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1039 } else { 1040 const CeedScalar *q_vec_array; 1041 CeedInt field_size; 1042 1043 // Copy Identity Outputs 1044 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 1045 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &q_vec_array)); 1046 for (CeedInt i = 0; i < field_size * num_points; i++) assembled_array[i] = q_vec_array[i]; 1047 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &q_vec_array)); 1048 assembled_array += field_size * num_points; 1049 } 1050 } 1051 num_points_offset += num_points; 1052 } 1053 1054 // Un-set output Qvecs to prevent accidental overwrite of Assembled 1055 if (!impl->is_identity_qf) { 1056 for (CeedInt out = 0; out < num_output_fields; out++) { 1057 CeedVector vec; 1058 1059 // Get output vector 1060 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1061 // Check if active output 1062 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 1063 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 1064 } 1065 } 1066 } 1067 1068 // Restore input arrays 1069 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 1070 1071 // Restore output 1072 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 1073 1074 // Cleanup 1075 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1076 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1077 return CEED_ERROR_SUCCESS; 1078 } 1079 1080 //------------------------------------------------------------------------------ 1081 // Assemble Linear QFunction 1082 //------------------------------------------------------------------------------ 1083 static int CeedOperatorLinearAssembleQFunctionAtPoints_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1084 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, true, assembled, rstr, request); 1085 } 1086 1087 //------------------------------------------------------------------------------ 1088 // Update Assembled Linear QFunction 1089 //------------------------------------------------------------------------------ 1090 static int CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, 1091 CeedRequest *request) { 1092 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, false, &assembled, &rstr, request); 1093 } 1094 1095 //------------------------------------------------------------------------------ 1096 // Assemble Operator Diagonal AtPoints 1097 //------------------------------------------------------------------------------ 1098 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1099 bool is_active_at_points = true; 1100 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, elem_size_active = 1, num_comp_active = 1; 1101 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 1102 Ceed ceed; 1103 CeedVector point_coords = NULL, in_vec, out_vec; 1104 CeedElemRestriction rstr_points = NULL; 1105 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1106 CeedQFunction qf; 1107 CeedOperatorField *op_input_fields, *op_output_fields; 1108 CeedOperator_Ref *impl; 1109 1110 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1111 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1112 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1113 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1114 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1115 1116 // Setup 1117 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 1118 1119 // Ceed 1120 { 1121 Ceed ceed_parent; 1122 1123 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1124 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 1125 if (ceed_parent) ceed = ceed_parent; 1126 } 1127 1128 // Point coordinates 1129 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1130 1131 // Input and output vectors 1132 { 1133 CeedSize input_size, output_size; 1134 1135 CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1136 CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec)); 1137 CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec)); 1138 CeedCallBackend(CeedVectorSetValue(out_vec, 0.0)); 1139 } 1140 1141 // Input Evecs and Restriction 1142 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 1143 1144 // Check if active field is at points 1145 for (CeedInt i = 0; i < num_input_fields; i++) { 1146 CeedRestrictionType rstr_type; 1147 CeedVector vec; 1148 CeedElemRestriction elem_rstr; 1149 1150 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1151 // Skip non-active input 1152 if (vec != CEED_VECTOR_ACTIVE) continue; 1153 1154 // Get active restriction type 1155 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1156 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1157 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); 1158 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; 1159 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); 1160 } 1161 1162 // Loop through elements 1163 for (CeedInt e = 0; e < num_elem; e++) { 1164 CeedInt num_points, e_vec_size = 0; 1165 1166 // Setup points for element 1167 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1168 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1169 1170 // Input basis apply for non-active bases 1171 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 1172 impl->point_coords_elem, true, e_data, impl, request)); 1173 1174 // Loop over points on element 1175 e_vec_size = (is_active_at_points ? num_points : elem_size_active) * num_comp_active; 1176 for (CeedInt s = 0; s < e_vec_size; s++) { 1177 for (CeedInt i = 0; i < num_input_fields; i++) { 1178 bool is_active_input = false; 1179 CeedInt size; 1180 CeedRestrictionType rstr_type; 1181 CeedEvalMode eval_mode; 1182 CeedVector vec; 1183 CeedElemRestriction elem_rstr; 1184 CeedBasis basis; 1185 1186 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1187 // Skip non-active input 1188 is_active_input = vec == CEED_VECTOR_ACTIVE; 1189 if (!is_active_input) continue; 1190 1191 // Get elem_size, eval_mode, size 1192 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1193 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1194 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1195 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 1196 // Update unit vector 1197 { 1198 CeedScalar *array; 1199 1200 if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0)); 1201 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array)); 1202 array[s] = 1.0; 1203 if (s > 0) array[s - 1] = 0.0; 1204 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array)); 1205 } 1206 // Basis action 1207 switch (eval_mode) { 1208 case CEED_EVAL_NONE: 1209 break; 1210 // Note - these basis eval modes require FEM fields 1211 case CEED_EVAL_INTERP: 1212 case CEED_EVAL_GRAD: 1213 case CEED_EVAL_DIV: 1214 case CEED_EVAL_CURL: 1215 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1216 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], 1217 impl->q_vecs_in[i])); 1218 break; 1219 case CEED_EVAL_WEIGHT: 1220 break; // No action 1221 } 1222 } 1223 1224 // -- Q function 1225 if (!impl->is_identity_qf) { 1226 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1227 } 1228 1229 // -- Output basis apply and restriction 1230 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 1231 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 1232 1233 // -- Grab diagonal value 1234 for (CeedInt i = 0; i < num_output_fields; i++) { 1235 bool is_active_output = false; 1236 CeedRestrictionType rstr_type; 1237 CeedEvalMode eval_mode; 1238 CeedVector vec; 1239 CeedElemRestriction elem_rstr; 1240 CeedBasis basis; 1241 1242 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1243 // ---- Skip non-active output 1244 is_active_output = vec == CEED_VECTOR_ACTIVE; 1245 if (!is_active_output) continue; 1246 1247 // ---- Get elem_size, eval_mode, size 1248 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1249 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1250 // ---- Basis action 1251 switch (eval_mode) { 1252 case CEED_EVAL_NONE: 1253 break; // No action 1254 case CEED_EVAL_INTERP: 1255 case CEED_EVAL_GRAD: 1256 case CEED_EVAL_DIV: 1257 case CEED_EVAL_CURL: 1258 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1259 CeedCallBackend(CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], 1260 impl->e_vecs_out[i])); 1261 break; 1262 // LCOV_EXCL_START 1263 case CEED_EVAL_WEIGHT: { 1264 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1265 // LCOV_EXCL_STOP 1266 } 1267 } 1268 // ---- Update output vector 1269 { 1270 CeedScalar *array, current_value = 0.0; 1271 1272 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, &array)); 1273 current_value = array[s]; 1274 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &array)); 1275 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_out[i], 0.0)); 1276 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, &array)); 1277 array[s] = current_value; 1278 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &array)); 1279 } 1280 // ---- Restrict output block 1281 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1282 if (rstr_type == CEED_RESTRICTION_POINTS) { 1283 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], assembled, request)); 1284 } else { 1285 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], assembled, request)); 1286 } 1287 } 1288 } 1289 num_points_offset += num_points; 1290 } 1291 1292 // Restore input arrays 1293 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 1294 1295 // Cleanup 1296 CeedCallBackend(CeedVectorDestroy(&in_vec)); 1297 CeedCallBackend(CeedVectorDestroy(&out_vec)); 1298 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1299 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1300 return CEED_ERROR_SUCCESS; 1301 } 1302 1303 //------------------------------------------------------------------------------ 1304 // Operator Destroy 1305 //------------------------------------------------------------------------------ 1306 static int CeedOperatorDestroy_Ref(CeedOperator op) { 1307 CeedOperator_Ref *impl; 1308 1309 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1310 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 1311 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 1312 } 1313 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 1314 CeedCallBackend(CeedFree(&impl->input_states)); 1315 1316 for (CeedInt i = 0; i < impl->num_inputs; i++) { 1317 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 1318 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 1319 } 1320 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 1321 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 1322 1323 for (CeedInt i = 0; i < impl->num_outputs; i++) { 1324 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 1325 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 1326 } 1327 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 1328 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 1329 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 1330 1331 // QFunction assembly 1332 for (CeedInt i = 0; i < impl->num_active_in; i++) { 1333 CeedCallBackend(CeedVectorDestroy(&impl->qf_active_in[i])); 1334 } 1335 CeedCallBackend(CeedFree(&impl->qf_active_in)); 1336 1337 CeedCallBackend(CeedFree(&impl)); 1338 return CEED_ERROR_SUCCESS; 1339 } 1340 1341 //------------------------------------------------------------------------------ 1342 // Operator Create 1343 //------------------------------------------------------------------------------ 1344 int CeedOperatorCreate_Ref(CeedOperator op) { 1345 Ceed ceed; 1346 CeedOperator_Ref *impl; 1347 1348 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1349 CeedCallBackend(CeedCalloc(1, &impl)); 1350 CeedCallBackend(CeedOperatorSetData(op, impl)); 1351 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 1352 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 1353 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 1354 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1355 return CEED_ERROR_SUCCESS; 1356 } 1357 1358 //------------------------------------------------------------------------------ 1359 // Operator Create At Points 1360 //------------------------------------------------------------------------------ 1361 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { 1362 Ceed ceed; 1363 CeedOperator_Ref *impl; 1364 1365 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1366 CeedCallBackend(CeedCalloc(1, &impl)); 1367 CeedCallBackend(CeedOperatorSetData(op, impl)); 1368 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Ref)); 1369 CeedCallBackend( 1370 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref)); 1371 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref)); 1372 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); 1373 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1374 return CEED_ERROR_SUCCESS; 1375 } 1376 1377 //------------------------------------------------------------------------------ 1378