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