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