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