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