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 CeedInt qf_size_in, qf_size_out, Q, num_elem, num_input_fields, num_output_fields; 386 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 387 CeedQFunctionField *qf_input_fields, *qf_output_fields; 388 CeedQFunction qf; 389 CeedOperatorField *op_input_fields, *op_output_fields; 390 CeedOperator_Ref *impl; 391 392 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 393 CeedCallBackend(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 394 CeedCallBackend(CeedOperatorGetData(op, &impl)); 395 qf_size_in = impl->qf_size_in; 396 qf_size_out = impl->qf_size_out; 397 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 398 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 399 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 400 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 401 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 402 403 // Setup 404 CeedCallBackend(CeedOperatorSetup_Ref(op)); 405 406 // Check for restriction only operator 407 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 408 409 // Input Evecs and Restriction 410 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 411 412 // Count number of active input fields 413 if (qf_size_in == 0) { 414 for (CeedInt i = 0; i < num_input_fields; i++) { 415 CeedInt field_size; 416 CeedVector vec; 417 418 // Get input vector 419 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 420 // Check if active input 421 if (vec == CEED_VECTOR_ACTIVE) { 422 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 423 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 424 qf_size_in += field_size; 425 } 426 } 427 CeedCheck(qf_size_in > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 428 impl->qf_size_in = qf_size_in; 429 } 430 431 // Count number of active output fields 432 if (qf_size_out == 0) { 433 for (CeedInt i = 0; i < num_output_fields; i++) { 434 CeedInt field_size; 435 CeedVector vec; 436 437 // Get output vector 438 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 439 // Check if active output 440 if (vec == CEED_VECTOR_ACTIVE) { 441 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 442 qf_size_out += field_size; 443 } 444 } 445 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 446 impl->qf_size_out = qf_size_out; 447 } 448 449 // Build objects if needed 450 if (build_objects) { 451 const CeedSize l_size = (CeedSize)num_elem * Q * qf_size_in * qf_size_out; 452 CeedInt strides[3] = {1, Q, qf_size_in * qf_size_out * Q}; /* *NOPAD* */ 453 454 // Create output restriction 455 CeedCallBackend( 456 CeedElemRestrictionCreateStrided(ceed_parent, num_elem, Q, qf_size_in * qf_size_out, qf_size_in * qf_size_out * num_elem * Q, strides, rstr)); 457 // Create assembled vector 458 CeedCallBackend(CeedVectorCreate(ceed_parent, l_size, assembled)); 459 } 460 // Clear output vector 461 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 462 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 463 464 // Loop through elements 465 for (CeedInt e = 0; e < num_elem; e++) { 466 // Input basis apply 467 CeedCallBackend(CeedOperatorInputBasis_Ref(e, Q, qf_input_fields, op_input_fields, num_input_fields, true, e_data_full, impl)); 468 469 // Assemble QFunction 470 471 for (CeedInt i = 0; i < num_input_fields; i++) { 472 CeedInt field_size; 473 CeedVector vec; 474 475 // Set Inputs 476 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 477 if (vec != CEED_VECTOR_ACTIVE) continue; 478 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 479 for (CeedInt field = 0; field < field_size; field++) { 480 // Set current portion of input to 1.0 481 { 482 CeedScalar *array; 483 484 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 485 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 1.0; 486 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 487 } 488 489 if (!impl->is_identity_qf) { 490 // Set Outputs 491 for (CeedInt out = 0; out < num_output_fields; out++) { 492 CeedVector vec; 493 494 // Get output vector 495 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 496 // Check if active output 497 if (vec == CEED_VECTOR_ACTIVE) { 498 CeedInt field_size; 499 500 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 501 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 502 assembled_array += field_size * Q; // Advance the pointer by the size of the output 503 } 504 } 505 // Apply QFunction 506 CeedCallBackend(CeedQFunctionApply(qf, Q, impl->q_vecs_in, impl->q_vecs_out)); 507 } else { 508 CeedInt field_size; 509 const CeedScalar *array; 510 511 // Copy Identity Outputs 512 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 513 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array)); 514 for (CeedInt j = 0; j < field_size * Q; j++) assembled_array[j] = array[j]; 515 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array)); 516 assembled_array += field_size * Q; 517 } 518 // Reset input to 0.0 519 { 520 CeedScalar *array; 521 522 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 523 for (CeedInt j = 0; j < Q; j++) array[field * Q + j] = 0.0; 524 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 525 } 526 } 527 } 528 } 529 530 // Un-set output Qvecs to prevent accidental overwrite of Assembled 531 if (!impl->is_identity_qf) { 532 for (CeedInt out = 0; out < num_output_fields; out++) { 533 CeedVector vec; 534 535 // Get output vector 536 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 537 // Check if active output 538 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 539 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 540 } 541 } 542 } 543 544 // Restore input arrays 545 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 546 547 // Restore output 548 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 549 return CEED_ERROR_SUCCESS; 550 } 551 552 //------------------------------------------------------------------------------ 553 // Assemble Linear QFunction 554 //------------------------------------------------------------------------------ 555 static int CeedOperatorLinearAssembleQFunction_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 556 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, true, assembled, rstr, request); 557 } 558 559 //------------------------------------------------------------------------------ 560 // Update Assembled Linear QFunction 561 //------------------------------------------------------------------------------ 562 static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, CeedRequest *request) { 563 return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); 564 } 565 566 //------------------------------------------------------------------------------ 567 // Setup Input/Output Fields 568 //------------------------------------------------------------------------------ 569 static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedVector *e_vecs_full, CeedVector *e_vecs, 570 CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, CeedInt Q) { 571 Ceed ceed; 572 CeedSize e_size, q_size; 573 CeedInt max_num_points, num_comp, size, P; 574 CeedQFunctionField *qf_fields; 575 CeedOperatorField *op_fields; 576 577 { 578 Ceed ceed_parent; 579 580 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 581 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 582 if (ceed_parent) ceed = ceed_parent; 583 } 584 if (is_input) { 585 CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 586 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 587 } else { 588 CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 589 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); 590 } 591 592 // Get max number of points 593 { 594 CeedInt dim; 595 CeedElemRestriction rstr_points = NULL; 596 CeedOperator_Ref *impl; 597 598 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 599 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 600 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); 601 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 602 CeedCallBackend(CeedOperatorGetData(op, &impl)); 603 if (is_input) { 604 CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); 605 CeedCallBackend(CeedVectorSetValue(impl->point_coords_elem, 0.0)); 606 } 607 } 608 609 // Loop over fields 610 for (CeedInt i = 0; i < num_fields; i++) { 611 CeedEvalMode eval_mode; 612 CeedBasis basis; 613 614 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 615 if (eval_mode != CEED_EVAL_WEIGHT) { 616 CeedElemRestriction elem_rstr; 617 618 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr)); 619 CeedCallBackend(CeedElemRestrictionCreateVector(elem_rstr, NULL, &e_vecs_full[i + start_e])); 620 CeedCallBackend(CeedVectorSetValue(e_vecs_full[i + start_e], 0.0)); 621 } 622 623 switch (eval_mode) { 624 case CEED_EVAL_NONE: { 625 CeedVector vec; 626 627 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 628 e_size = (CeedSize)max_num_points * size; 629 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 630 CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &vec)); 631 if (vec == CEED_VECTOR_ACTIVE || !is_input) { 632 CeedCallBackend(CeedVectorReferenceCopy(e_vecs[i], &q_vecs[i])); 633 } else { 634 q_size = (CeedSize)max_num_points * size; 635 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 636 } 637 break; 638 } 639 case CEED_EVAL_INTERP: 640 case CEED_EVAL_GRAD: 641 case CEED_EVAL_DIV: 642 case CEED_EVAL_CURL: 643 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 644 CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); 645 CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); 646 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 647 e_size = (CeedSize)P * num_comp; 648 CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); 649 q_size = (CeedSize)max_num_points * size; 650 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 651 break; 652 case CEED_EVAL_WEIGHT: // Only on input fields 653 CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 654 q_size = (CeedSize)max_num_points; 655 CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); 656 CeedCallBackend( 657 CeedBasisApplyAtPoints(basis, 1, &max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); 658 break; 659 } 660 // Initialize full arrays for E-vectors and Q-vectors 661 if (e_vecs[i]) CeedCallBackend(CeedVectorSetValue(e_vecs[i], 0.0)); 662 if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorSetValue(q_vecs[i], 0.0)); 663 } 664 return CEED_ERROR_SUCCESS; 665 } 666 667 //------------------------------------------------------------------------------ 668 // Setup Operator 669 //------------------------------------------------------------------------------ 670 static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { 671 bool is_setup_done; 672 CeedInt Q, num_input_fields, num_output_fields; 673 CeedQFunctionField *qf_input_fields, *qf_output_fields; 674 CeedQFunction qf; 675 CeedOperatorField *op_input_fields, *op_output_fields; 676 CeedOperator_Ref *impl; 677 678 CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); 679 if (is_setup_done) return CEED_ERROR_SUCCESS; 680 681 CeedCallBackend(CeedOperatorGetData(op, &impl)); 682 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 683 CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); 684 CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); 685 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 686 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 687 688 // Allocate 689 CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); 690 691 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); 692 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); 693 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); 694 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); 695 CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); 696 697 impl->num_inputs = num_input_fields; 698 impl->num_outputs = num_output_fields; 699 700 // Set up infield and outfield pointer arrays 701 // Infields 702 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, num_input_fields, Q)); 703 // Outfields 704 CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, num_input_fields, 705 num_output_fields, Q)); 706 707 // Identity QFunctions 708 if (impl->is_identity_qf) { 709 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); 710 CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->e_vecs_out[0])); 711 } 712 713 CeedCallBackend(CeedOperatorSetSetupDone(op)); 714 return CEED_ERROR_SUCCESS; 715 } 716 717 //------------------------------------------------------------------------------ 718 // Input Basis Action 719 //------------------------------------------------------------------------------ 720 static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, 721 CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, 722 CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], 723 CeedOperator_Ref *impl, CeedRequest *request) { 724 for (CeedInt i = 0; i < num_input_fields; i++) { 725 bool is_active_input = false; 726 CeedInt elem_size, size, num_comp; 727 CeedRestrictionType rstr_type; 728 CeedEvalMode eval_mode; 729 CeedVector vec; 730 CeedElemRestriction elem_rstr; 731 CeedBasis basis; 732 733 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 734 // Skip active input 735 is_active_input = vec == CEED_VECTOR_ACTIVE; 736 if (skip_active && is_active_input) continue; 737 738 // Get elem_size, eval_mode, size 739 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 740 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 741 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 742 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 743 // Restrict block active input 744 if (is_active_input) { 745 if (rstr_type == CEED_RESTRICTION_POINTS) { 746 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 747 } else { 748 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); 749 } 750 } 751 // Basis action 752 switch (eval_mode) { 753 case CEED_EVAL_NONE: 754 if (!is_active_input) { 755 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); 756 } 757 break; 758 // Note - these basis eval modes require FEM fields 759 case CEED_EVAL_INTERP: 760 case CEED_EVAL_GRAD: 761 case CEED_EVAL_DIV: 762 case CEED_EVAL_CURL: 763 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 764 if (!is_active_input) { 765 CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); 766 CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); 767 CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][(CeedSize)e * elem_size * num_comp])); 768 } 769 CeedCallBackend( 770 CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); 771 break; 772 case CEED_EVAL_WEIGHT: 773 break; // No action 774 } 775 } 776 return CEED_ERROR_SUCCESS; 777 } 778 779 //------------------------------------------------------------------------------ 780 // Output Basis Action 781 //------------------------------------------------------------------------------ 782 static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, 783 CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, 784 CeedOperator op, CeedVector out_vec, CeedVector point_coords_elem, CeedOperator_Ref *impl, 785 CeedRequest *request) { 786 for (CeedInt i = 0; i < num_output_fields; i++) { 787 CeedRestrictionType rstr_type; 788 CeedEvalMode eval_mode; 789 CeedVector vec; 790 CeedElemRestriction elem_rstr; 791 CeedBasis basis; 792 793 // Get elem_size, eval_mode, size 794 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 795 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 796 // Basis action 797 switch (eval_mode) { 798 case CEED_EVAL_NONE: 799 break; // No action 800 case CEED_EVAL_INTERP: 801 case CEED_EVAL_GRAD: 802 case CEED_EVAL_DIV: 803 case CEED_EVAL_CURL: 804 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 805 CeedCallBackend( 806 CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); 807 break; 808 // LCOV_EXCL_START 809 case CEED_EVAL_WEIGHT: { 810 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 811 // LCOV_EXCL_STOP 812 } 813 } 814 // Restrict output block 815 // Get output vector 816 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 817 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 818 if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; 819 // Restrict 820 if (rstr_type == CEED_RESTRICTION_POINTS) { 821 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 822 } else { 823 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); 824 } 825 } 826 return CEED_ERROR_SUCCESS; 827 } 828 829 //------------------------------------------------------------------------------ 830 // Operator Apply 831 //------------------------------------------------------------------------------ 832 static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { 833 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; 834 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 835 CeedVector point_coords = NULL; 836 CeedElemRestriction rstr_points = NULL; 837 CeedQFunctionField *qf_input_fields, *qf_output_fields; 838 CeedQFunction qf; 839 CeedOperatorField *op_input_fields, *op_output_fields; 840 CeedOperator_Ref *impl; 841 842 CeedCallBackend(CeedOperatorGetData(op, &impl)); 843 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 844 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 845 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 846 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 847 848 // Setup 849 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 850 851 // Point coordinates 852 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 853 854 // Input Evecs and Restriction 855 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 856 857 // Loop through elements 858 for (CeedInt e = 0; e < num_elem; e++) { 859 CeedInt num_points; 860 861 // Setup points for element 862 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 863 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 864 865 // Input basis apply 866 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 867 impl->point_coords_elem, false, e_data, impl, request)); 868 869 // Q function 870 if (!impl->is_identity_qf) { 871 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 872 } 873 874 // Output basis apply and restriction 875 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 876 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 877 878 num_points_offset += num_points; 879 } 880 881 // Restore input arrays 882 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 883 884 // Cleanup point coordinates 885 CeedCallBackend(CeedVectorDestroy(&point_coords)); 886 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 887 return CEED_ERROR_SUCCESS; 888 } 889 890 //------------------------------------------------------------------------------ 891 // Core code for assembling linear QFunction 892 //------------------------------------------------------------------------------ 893 static inline int CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(CeedOperator op, bool build_objects, CeedVector *assembled, 894 CeedElemRestriction *rstr, CeedRequest *request) { 895 Ceed ceed; 896 CeedInt qf_size_in, qf_size_out, max_num_points, num_elem, num_input_fields, num_output_fields, num_points_offset = 0; 897 CeedScalar *assembled_array, *e_data_full[2 * CEED_FIELD_MAX] = {NULL}; 898 CeedVector point_coords = NULL; 899 CeedQFunctionField *qf_input_fields, *qf_output_fields; 900 CeedQFunction qf; 901 CeedOperatorField *op_input_fields, *op_output_fields; 902 CeedOperator_Ref *impl; 903 CeedElemRestriction rstr_points = NULL; 904 905 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 906 CeedCallBackend(CeedOperatorGetData(op, &impl)); 907 qf_size_in = impl->qf_size_in; 908 qf_size_out = impl->qf_size_out; 909 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 910 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 911 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 912 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 913 914 // Setup 915 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 916 917 // Check for restriction only operator 918 CeedCheck(!impl->is_identity_rstr_op, ceed, CEED_ERROR_BACKEND, "Assembling restriction only operators is not supported"); 919 920 // Point coordinates 921 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 922 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); 923 924 // Input Evecs and Restriction 925 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data_full, impl, request)); 926 927 // Count number of active input fields 928 if (qf_size_in == 0) { 929 for (CeedInt i = 0; i < num_input_fields; i++) { 930 CeedInt field_size; 931 CeedVector vec; 932 933 // Get input vector 934 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 935 // Check if active input 936 if (vec == CEED_VECTOR_ACTIVE) { 937 // Check that all active inputs are nodal fields 938 { 939 CeedElemRestriction elem_rstr; 940 bool is_at_points = false; 941 942 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 943 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 944 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 945 } 946 // Get size of active input 947 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 948 qf_size_in += field_size; 949 } 950 } 951 CeedCheck(qf_size_in, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 952 impl->qf_size_in = qf_size_in; 953 } 954 955 // Count number of active output fields 956 if (qf_size_out == 0) { 957 for (CeedInt i = 0; i < num_output_fields; i++) { 958 CeedInt field_size; 959 CeedVector vec; 960 961 // Get output vector 962 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 963 // Check if active output 964 if (vec == CEED_VECTOR_ACTIVE) { 965 // Check that all active inputs are nodal fields 966 { 967 CeedElemRestriction elem_rstr; 968 bool is_at_points = false; 969 970 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 971 CeedCallBackend(CeedElemRestrictionIsAtPoints(elem_rstr, &is_at_points)); 972 CeedCheck(!is_at_points, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction with active input at points"); 973 } 974 // Get size of active output 975 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[i], &field_size)); 976 CeedCallBackend(CeedVectorSetValue(impl->q_vecs_in[i], 0.0)); 977 qf_size_out += field_size; 978 } 979 } 980 CeedCheck(qf_size_out > 0, ceed, CEED_ERROR_BACKEND, "Cannot assemble QFunction without active inputs and outputs"); 981 impl->qf_size_out = qf_size_out; 982 } 983 984 // Build objects if needed 985 if (build_objects) { 986 CeedInt num_points_total; 987 const CeedInt *offsets; 988 989 CeedCallBackend(CeedElemRestrictionGetNumPoints(rstr_points, &num_points_total)); 990 991 // Create output restriction (at points) 992 CeedCallBackend(CeedElemRestrictionGetOffsets(rstr_points, CEED_MEM_HOST, &offsets)); 993 CeedCallBackend(CeedElemRestrictionCreateAtPoints(ceed, num_elem, num_points_total, qf_size_in * qf_size_out, 994 qf_size_in * qf_size_out * num_points_total, CEED_MEM_HOST, CEED_COPY_VALUES, offsets, rstr)); 995 CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr_points, &offsets)); 996 997 // Create assembled vector 998 CeedCallBackend(CeedElemRestrictionCreateVector(*rstr, assembled, NULL)); 999 } 1000 // Clear output vector 1001 CeedCallBackend(CeedVectorSetValue(*assembled, 0.0)); 1002 CeedCallBackend(CeedVectorGetArray(*assembled, CEED_MEM_HOST, &assembled_array)); 1003 1004 // Loop through elements 1005 for (CeedInt e = 0; e < num_elem; e++) { 1006 CeedInt num_points; 1007 1008 // Setup points for element 1009 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1010 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1011 1012 // Input basis apply 1013 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, NULL, 1014 impl->point_coords_elem, true, e_data_full, impl, request)); 1015 1016 // Assemble QFunction 1017 for (CeedInt i = 0; i < num_input_fields; i++) { 1018 CeedInt field_size; 1019 CeedVector vec; 1020 1021 // Get input vector 1022 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1023 // Check if active input 1024 if (vec != CEED_VECTOR_ACTIVE) continue; 1025 // Get size of active input 1026 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &field_size)); 1027 for (CeedInt field = 0; field < field_size; field++) { 1028 // Set current portion of input to 1.0 1029 { 1030 CeedScalar *array; 1031 1032 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 1033 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 1.0; 1034 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 1035 } 1036 1037 if (!impl->is_identity_qf) { 1038 // Set Outputs 1039 for (CeedInt out = 0; out < num_output_fields; out++) { 1040 CeedVector vec; 1041 CeedInt field_size; 1042 1043 // Get output vector 1044 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1045 // Check if active output 1046 if (vec == CEED_VECTOR_ACTIVE) { 1047 CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[out], CEED_MEM_HOST, CEED_USE_POINTER, assembled_array)); 1048 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[out], &field_size)); 1049 assembled_array += field_size * num_points; // Advance the pointer by the size of the output 1050 } 1051 } 1052 // Apply QFunction 1053 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1054 } else { 1055 const CeedScalar *array; 1056 CeedInt field_size; 1057 1058 // Copy Identity Outputs 1059 CeedCallBackend(CeedQFunctionFieldGetSize(qf_output_fields[0], &field_size)); 1060 CeedCallBackend(CeedVectorGetArrayRead(impl->q_vecs_out[0], CEED_MEM_HOST, &array)); 1061 for (CeedInt j = 0; j < field_size * num_points; j++) assembled_array[j] = array[j]; 1062 CeedCallBackend(CeedVectorRestoreArrayRead(impl->q_vecs_out[0], &array)); 1063 assembled_array += field_size * num_points; 1064 } 1065 // Reset input to 0.0 1066 { 1067 CeedScalar *array; 1068 1069 CeedCallBackend(CeedVectorGetArray(impl->q_vecs_in[i], CEED_MEM_HOST, &array)); 1070 for (CeedInt j = 0; j < num_points; j++) array[field * num_points + j] = 0.0; 1071 CeedCallBackend(CeedVectorRestoreArray(impl->q_vecs_in[i], &array)); 1072 } 1073 } 1074 } 1075 num_points_offset += num_points; 1076 } 1077 1078 // Un-set output Qvecs to prevent accidental overwrite of Assembled 1079 if (!impl->is_identity_qf) { 1080 for (CeedInt out = 0; out < num_output_fields; out++) { 1081 CeedVector vec; 1082 1083 // Get output vector 1084 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[out], &vec)); 1085 // Check if active output 1086 if (vec == CEED_VECTOR_ACTIVE && num_elem > 0) { 1087 CeedCallBackend(CeedVectorTakeArray(impl->q_vecs_out[out], CEED_MEM_HOST, NULL)); 1088 } 1089 } 1090 } 1091 1092 // Restore input arrays 1093 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data_full, impl)); 1094 1095 // Restore output 1096 CeedCallBackend(CeedVectorRestoreArray(*assembled, &assembled_array)); 1097 1098 // Cleanup 1099 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1100 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1101 return CEED_ERROR_SUCCESS; 1102 } 1103 1104 //------------------------------------------------------------------------------ 1105 // Assemble Linear QFunction 1106 //------------------------------------------------------------------------------ 1107 static int CeedOperatorLinearAssembleQFunctionAtPoints_Ref(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1108 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, true, assembled, rstr, request); 1109 } 1110 1111 //------------------------------------------------------------------------------ 1112 // Update Assembled Linear QFunction 1113 //------------------------------------------------------------------------------ 1114 static int CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref(CeedOperator op, CeedVector assembled, CeedElemRestriction rstr, 1115 CeedRequest *request) { 1116 return CeedOperatorLinearAssembleQFunctionAtPointsCore_Ref(op, false, &assembled, &rstr, request); 1117 } 1118 1119 //------------------------------------------------------------------------------ 1120 // Assemble Operator Diagonal AtPoints 1121 //------------------------------------------------------------------------------ 1122 static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1123 bool is_active_at_points = true; 1124 CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem, elem_size_active = 1, num_comp_active = 1; 1125 CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; 1126 Ceed ceed; 1127 CeedVector point_coords = NULL, in_vec, out_vec; 1128 CeedElemRestriction rstr_points = NULL; 1129 CeedQFunctionField *qf_input_fields, *qf_output_fields; 1130 CeedQFunction qf; 1131 CeedOperatorField *op_input_fields, *op_output_fields; 1132 CeedOperator_Ref *impl; 1133 1134 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1135 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 1136 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 1137 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 1138 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 1139 1140 // Setup 1141 CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); 1142 1143 // Ceed 1144 { 1145 Ceed ceed_parent; 1146 1147 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1148 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 1149 if (ceed_parent) ceed = ceed_parent; 1150 } 1151 1152 // Point coordinates 1153 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); 1154 1155 // Input and output vectors 1156 { 1157 CeedSize input_size, output_size; 1158 1159 CeedCallBackend(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1160 CeedCallBackend(CeedVectorCreate(ceed, input_size, &in_vec)); 1161 CeedCallBackend(CeedVectorCreate(ceed, output_size, &out_vec)); 1162 CeedCallBackend(CeedVectorSetValue(out_vec, 0.0)); 1163 } 1164 1165 // Input Evecs and Restriction 1166 CeedCallBackend(CeedOperatorSetupInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, NULL, true, e_data, impl, request)); 1167 1168 // Check if active field is at points 1169 for (CeedInt i = 0; i < num_input_fields; i++) { 1170 CeedRestrictionType rstr_type; 1171 CeedVector vec; 1172 CeedElemRestriction elem_rstr; 1173 1174 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1175 // Skip non-active input 1176 if (vec != CEED_VECTOR_ACTIVE) continue; 1177 1178 // Get active restriction type 1179 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1180 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1181 CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp_active)); 1182 is_active_at_points = rstr_type == CEED_RESTRICTION_POINTS; 1183 if (!is_active_at_points) CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size_active)); 1184 } 1185 1186 // Loop through elements 1187 for (CeedInt e = 0; e < num_elem; e++) { 1188 CeedInt num_points, e_vec_size = 0; 1189 1190 // Setup points for element 1191 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); 1192 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); 1193 1194 // Input basis apply for non-active bases 1195 CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, 1196 impl->point_coords_elem, true, e_data, impl, request)); 1197 1198 // Loop over points on element 1199 e_vec_size = (is_active_at_points ? num_points : elem_size_active) * num_comp_active; 1200 for (CeedInt s = 0; s < e_vec_size; s++) { 1201 for (CeedInt i = 0; i < num_input_fields; i++) { 1202 bool is_active_input = false; 1203 CeedInt size; 1204 CeedRestrictionType rstr_type; 1205 CeedEvalMode eval_mode; 1206 CeedVector vec; 1207 CeedElemRestriction elem_rstr; 1208 CeedBasis basis; 1209 1210 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 1211 // Skip non-active input 1212 is_active_input = vec == CEED_VECTOR_ACTIVE; 1213 if (!is_active_input) continue; 1214 1215 // Get elem_size, eval_mode, size 1216 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); 1217 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1218 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 1219 CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); 1220 // Update unit vector 1221 { 1222 CeedScalar *array; 1223 1224 if (s == 0) CeedCallBackend(CeedVectorSetValue(impl->e_vecs_in[i], 0.0)); 1225 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_in[i], CEED_MEM_HOST, &array)); 1226 array[s] = 1.0; 1227 if (s > 0) array[s - 1] = 0.0; 1228 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_in[i], &array)); 1229 } 1230 // Basis action 1231 switch (eval_mode) { 1232 case CEED_EVAL_NONE: 1233 break; 1234 // Note - these basis eval modes require FEM fields 1235 case CEED_EVAL_INTERP: 1236 case CEED_EVAL_GRAD: 1237 case CEED_EVAL_DIV: 1238 case CEED_EVAL_CURL: 1239 CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); 1240 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], 1241 impl->q_vecs_in[i])); 1242 break; 1243 case CEED_EVAL_WEIGHT: 1244 break; // No action 1245 } 1246 } 1247 1248 // -- Q function 1249 if (!impl->is_identity_qf) { 1250 CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); 1251 } 1252 1253 // -- Output basis apply and restriction 1254 CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, 1255 num_output_fields, op, out_vec, impl->point_coords_elem, impl, request)); 1256 1257 // -- Grab diagonal value 1258 for (CeedInt i = 0; i < num_output_fields; i++) { 1259 bool is_active_output = false; 1260 CeedRestrictionType rstr_type; 1261 CeedEvalMode eval_mode; 1262 CeedVector vec; 1263 CeedElemRestriction elem_rstr; 1264 CeedBasis basis; 1265 1266 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 1267 // ---- Skip non-active output 1268 is_active_output = vec == CEED_VECTOR_ACTIVE; 1269 if (!is_active_output) continue; 1270 1271 // ---- Get elem_size, eval_mode, size 1272 CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); 1273 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 1274 // ---- Basis action 1275 switch (eval_mode) { 1276 case CEED_EVAL_NONE: 1277 break; // No action 1278 case CEED_EVAL_INTERP: 1279 case CEED_EVAL_GRAD: 1280 case CEED_EVAL_DIV: 1281 case CEED_EVAL_CURL: 1282 CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); 1283 CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], 1284 impl->e_vecs_out[i])); 1285 break; 1286 // LCOV_EXCL_START 1287 case CEED_EVAL_WEIGHT: { 1288 return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); 1289 // LCOV_EXCL_STOP 1290 } 1291 } 1292 // ---- Update output vector 1293 { 1294 CeedScalar *array, current_value = 0.0; 1295 1296 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, &array)); 1297 current_value = array[s]; 1298 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &array)); 1299 CeedCallBackend(CeedVectorSetValue(impl->e_vecs_out[i], 0.0)); 1300 CeedCallBackend(CeedVectorGetArray(impl->e_vecs_out[i], CEED_MEM_HOST, &array)); 1301 array[s] = current_value; 1302 CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &array)); 1303 } 1304 // ---- Restrict output block 1305 CeedCallBackend(CeedElemRestrictionGetType(elem_rstr, &rstr_type)); 1306 if (rstr_type == CEED_RESTRICTION_POINTS) { 1307 CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], assembled, request)); 1308 } else { 1309 CeedCallBackend(CeedElemRestrictionApplyBlock(elem_rstr, e, CEED_TRANSPOSE, impl->e_vecs_out[i], assembled, request)); 1310 } 1311 } 1312 } 1313 num_points_offset += num_points; 1314 } 1315 1316 // Restore input arrays 1317 CeedCallBackend(CeedOperatorRestoreInputs_Ref(num_input_fields, qf_input_fields, op_input_fields, true, e_data, impl)); 1318 1319 // Cleanup 1320 CeedCallBackend(CeedVectorDestroy(&in_vec)); 1321 CeedCallBackend(CeedVectorDestroy(&out_vec)); 1322 CeedCallBackend(CeedVectorDestroy(&point_coords)); 1323 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 1324 return CEED_ERROR_SUCCESS; 1325 } 1326 1327 //------------------------------------------------------------------------------ 1328 // Operator Destroy 1329 //------------------------------------------------------------------------------ 1330 static int CeedOperatorDestroy_Ref(CeedOperator op) { 1331 CeedOperator_Ref *impl; 1332 1333 CeedCallBackend(CeedOperatorGetData(op, &impl)); 1334 for (CeedInt i = 0; i < impl->num_inputs + impl->num_outputs; i++) { 1335 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_full[i])); 1336 } 1337 CeedCallBackend(CeedFree(&impl->e_vecs_full)); 1338 CeedCallBackend(CeedFree(&impl->input_states)); 1339 1340 for (CeedInt i = 0; i < impl->num_inputs; i++) { 1341 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_in[i])); 1342 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_in[i])); 1343 } 1344 CeedCallBackend(CeedFree(&impl->e_vecs_in)); 1345 CeedCallBackend(CeedFree(&impl->q_vecs_in)); 1346 1347 for (CeedInt i = 0; i < impl->num_outputs; i++) { 1348 CeedCallBackend(CeedVectorDestroy(&impl->e_vecs_out[i])); 1349 CeedCallBackend(CeedVectorDestroy(&impl->q_vecs_out[i])); 1350 } 1351 CeedCallBackend(CeedFree(&impl->e_vecs_out)); 1352 CeedCallBackend(CeedFree(&impl->q_vecs_out)); 1353 CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); 1354 1355 CeedCallBackend(CeedFree(&impl)); 1356 return CEED_ERROR_SUCCESS; 1357 } 1358 1359 //------------------------------------------------------------------------------ 1360 // Operator Create 1361 //------------------------------------------------------------------------------ 1362 int CeedOperatorCreate_Ref(CeedOperator op) { 1363 Ceed ceed; 1364 CeedOperator_Ref *impl; 1365 1366 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1367 CeedCallBackend(CeedCalloc(1, &impl)); 1368 CeedCallBackend(CeedOperatorSetData(op, impl)); 1369 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunction_Ref)); 1370 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionUpdate_Ref)); 1371 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Ref)); 1372 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1373 return CEED_ERROR_SUCCESS; 1374 } 1375 1376 //------------------------------------------------------------------------------ 1377 // Operator Create At Points 1378 //------------------------------------------------------------------------------ 1379 int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { 1380 Ceed ceed; 1381 CeedOperator_Ref *impl; 1382 1383 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 1384 CeedCallBackend(CeedCalloc(1, &impl)); 1385 CeedCallBackend(CeedOperatorSetData(op, impl)); 1386 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunction", CeedOperatorLinearAssembleQFunctionAtPoints_Ref)); 1387 CeedCallBackend( 1388 CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleQFunctionUpdate", CeedOperatorLinearAssembleQFunctionAtPointsUpdate_Ref)); 1389 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "LinearAssembleAddDiagonal", CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref)); 1390 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); 1391 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); 1392 return CEED_ERROR_SUCCESS; 1393 } 1394 1395 //------------------------------------------------------------------------------ 1396