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 <assert.h> 9 #include <ceed-impl.h> 10 #include <ceed/backend.h> 11 #include <ceed/ceed.h> 12 #include <math.h> 13 #include <stdbool.h> 14 #include <stdio.h> 15 #include <string.h> 16 17 /// @file 18 /// Implementation of CeedOperator preconditioning interfaces 19 20 /// ---------------------------------------------------------------------------- 21 /// CeedOperator Library Internal Preconditioning Functions 22 /// ---------------------------------------------------------------------------- 23 /// @addtogroup CeedOperatorDeveloper 24 /// @{ 25 26 /** 27 @brief Duplicate a CeedQFunction with a reference Ceed to fallback for advanced CeedOperator functionality 28 29 @param[in] fallback_ceed Ceed on which to create fallback CeedQFunction 30 @param[in] qf CeedQFunction to create fallback for 31 @param[out] qf_fallback fallback CeedQFunction 32 33 @return An error code: 0 - success, otherwise - failure 34 35 @ref Developer 36 **/ 37 static int CeedQFunctionCreateFallback(Ceed fallback_ceed, CeedQFunction qf, CeedQFunction *qf_fallback) { 38 // Check if NULL qf passed in 39 if (!qf) return CEED_ERROR_SUCCESS; 40 41 CeedDebug256(qf->ceed, 1, "---------- CeedOperator Fallback ----------\n"); 42 CeedDebug(qf->ceed, "Creating fallback CeedQFunction\n"); 43 44 char *source_path_with_name = ""; 45 if (qf->source_path) { 46 size_t path_len = strlen(qf->source_path), name_len = strlen(qf->kernel_name); 47 CeedCall(CeedCalloc(path_len + name_len + 2, &source_path_with_name)); 48 memcpy(source_path_with_name, qf->source_path, path_len); 49 memcpy(&source_path_with_name[path_len], ":", 1); 50 memcpy(&source_path_with_name[path_len + 1], qf->kernel_name, name_len); 51 } else { 52 CeedCall(CeedCalloc(1, &source_path_with_name)); 53 } 54 55 CeedCall(CeedQFunctionCreateInterior(fallback_ceed, qf->vec_length, qf->function, source_path_with_name, qf_fallback)); 56 { 57 CeedQFunctionContext ctx; 58 59 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 60 CeedCall(CeedQFunctionSetContext(*qf_fallback, ctx)); 61 } 62 for (CeedInt i = 0; i < qf->num_input_fields; i++) { 63 CeedCall(CeedQFunctionAddInput(*qf_fallback, qf->input_fields[i]->field_name, qf->input_fields[i]->size, qf->input_fields[i]->eval_mode)); 64 } 65 for (CeedInt i = 0; i < qf->num_output_fields; i++) { 66 CeedCall(CeedQFunctionAddOutput(*qf_fallback, qf->output_fields[i]->field_name, qf->output_fields[i]->size, qf->output_fields[i]->eval_mode)); 67 } 68 CeedCall(CeedFree(&source_path_with_name)); 69 70 return CEED_ERROR_SUCCESS; 71 } 72 73 /** 74 @brief Duplicate a CeedOperator with a reference Ceed to fallback for advanced CeedOperator functionality 75 76 @param[in,out] op CeedOperator to create fallback for 77 78 @return An error code: 0 - success, otherwise - failure 79 80 @ref Developer 81 **/ 82 static int CeedOperatorCreateFallback(CeedOperator op) { 83 bool is_composite; 84 Ceed ceed_fallback; 85 86 // Check not already created 87 if (op->op_fallback) return CEED_ERROR_SUCCESS; 88 89 // Fallback Ceed 90 CeedCall(CeedGetOperatorFallbackCeed(op->ceed, &ceed_fallback)); 91 if (!ceed_fallback) return CEED_ERROR_SUCCESS; 92 93 CeedDebug256(op->ceed, 1, "---------- CeedOperator Fallback ----------\n"); 94 CeedDebug(op->ceed, "Creating fallback CeedOperator\n"); 95 96 // Clone Op 97 CeedOperator op_fallback; 98 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 99 if (is_composite) { 100 CeedInt num_suboperators; 101 CeedOperator *sub_operators; 102 103 CeedCall(CeedCompositeOperatorCreate(ceed_fallback, &op_fallback)); 104 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 105 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 106 for (CeedInt i = 0; i < num_suboperators; i++) { 107 CeedOperator op_sub_fallback; 108 109 CeedCall(CeedOperatorGetFallback(sub_operators[i], &op_sub_fallback)); 110 CeedCall(CeedCompositeOperatorAddSub(op_fallback, op_sub_fallback)); 111 } 112 } else { 113 CeedQFunction qf_fallback = NULL, dqf_fallback = NULL, dqfT_fallback = NULL; 114 CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->qf, &qf_fallback)); 115 CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqf, &dqf_fallback)); 116 CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqfT, &dqfT_fallback)); 117 CeedCall(CeedOperatorCreate(ceed_fallback, qf_fallback, dqf_fallback, dqfT_fallback, &op_fallback)); 118 for (CeedInt i = 0; i < op->qf->num_input_fields; i++) { 119 CeedCall(CeedOperatorSetField(op_fallback, op->input_fields[i]->field_name, op->input_fields[i]->elem_rstr, op->input_fields[i]->basis, 120 op->input_fields[i]->vec)); 121 } 122 for (CeedInt i = 0; i < op->qf->num_output_fields; i++) { 123 CeedCall(CeedOperatorSetField(op_fallback, op->output_fields[i]->field_name, op->output_fields[i]->elem_rstr, op->output_fields[i]->basis, 124 op->output_fields[i]->vec)); 125 } 126 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(op->qf_assembled, &op_fallback->qf_assembled)); 127 if (op_fallback->num_qpts == 0) { 128 CeedCall(CeedOperatorSetNumQuadraturePoints(op_fallback, op->num_qpts)); 129 } 130 // Cleanup 131 CeedCall(CeedQFunctionDestroy(&qf_fallback)); 132 CeedCall(CeedQFunctionDestroy(&dqf_fallback)); 133 CeedCall(CeedQFunctionDestroy(&dqfT_fallback)); 134 } 135 CeedCall(CeedOperatorSetName(op_fallback, op->name)); 136 CeedCall(CeedOperatorCheckReady(op_fallback)); 137 op->op_fallback = op_fallback; 138 139 return CEED_ERROR_SUCCESS; 140 } 141 142 /** 143 @brief Retrieve fallback CeedOperator with a reference Ceed for advanced CeedOperator functionality 144 145 @param[in] op CeedOperator to retrieve fallback for 146 @param[out] op_fallback Fallback CeedOperator 147 148 @return An error code: 0 - success, otherwise - failure 149 150 @ref Developer 151 **/ 152 int CeedOperatorGetFallback(CeedOperator op, CeedOperator *op_fallback) { 153 // Create if needed 154 if (!op->op_fallback) { 155 CeedCall(CeedOperatorCreateFallback(op)); 156 } 157 if (op->op_fallback) { 158 bool is_debug; 159 160 CeedCall(CeedIsDebug(op->ceed, &is_debug)); 161 if (is_debug) { 162 Ceed ceed, ceed_fallback; 163 const char *resource, *resource_fallback; 164 165 CeedCall(CeedOperatorGetCeed(op, &ceed)); 166 CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback)); 167 CeedCall(CeedGetResource(ceed, &resource)); 168 CeedCall(CeedGetResource(ceed_fallback, &resource_fallback)); 169 170 CeedDebug256(ceed, 1, "---------- CeedOperator Fallback ----------\n"); 171 CeedDebug(ceed, "Falling back from %s operator at address %ld to %s operator at address %ld\n", resource, op, resource_fallback, 172 op->op_fallback); 173 } 174 } 175 *op_fallback = op->op_fallback; 176 177 return CEED_ERROR_SUCCESS; 178 } 179 180 /** 181 @brief Select correct basis matrix pointer based on CeedEvalMode 182 183 @param[in] eval_mode Current basis evaluation mode 184 @param[in] identity Pointer to identity matrix 185 @param[in] interp Pointer to interpolation matrix 186 @param[in] grad Pointer to gradient matrix 187 @param[out] basis_ptr Basis pointer to set 188 189 @ref Developer 190 **/ 191 static inline void CeedOperatorGetBasisPointer(CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar *interp, const CeedScalar *grad, 192 const CeedScalar **basis_ptr) { 193 switch (eval_mode) { 194 case CEED_EVAL_NONE: 195 *basis_ptr = identity; 196 break; 197 case CEED_EVAL_INTERP: 198 *basis_ptr = interp; 199 break; 200 case CEED_EVAL_GRAD: 201 *basis_ptr = grad; 202 break; 203 case CEED_EVAL_WEIGHT: 204 case CEED_EVAL_DIV: 205 case CEED_EVAL_CURL: 206 break; // Caught by QF Assembly 207 } 208 assert(*basis_ptr != NULL); 209 } 210 211 /** 212 @brief Create point block restriction for active operator field 213 214 @param[in] rstr Original CeedElemRestriction for active field 215 @param[out] pointblock_rstr Address of the variable where the newly created CeedElemRestriction will be stored 216 217 @return An error code: 0 - success, otherwise - failure 218 219 @ref Developer 220 **/ 221 static int CeedOperatorCreateActivePointBlockRestriction(CeedElemRestriction rstr, CeedElemRestriction *pointblock_rstr) { 222 Ceed ceed; 223 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 224 const CeedInt *offsets; 225 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 226 227 // Expand offsets 228 CeedInt num_elem, num_comp, elem_size, comp_stride, *pointblock_offsets; 229 CeedSize l_size; 230 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 231 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 232 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 233 CeedCall(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 234 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 235 CeedInt shift = num_comp; 236 if (comp_stride != 1) shift *= num_comp; 237 CeedCall(CeedCalloc(num_elem * elem_size, &pointblock_offsets)); 238 for (CeedInt i = 0; i < num_elem * elem_size; i++) { 239 pointblock_offsets[i] = offsets[i] * shift; 240 } 241 242 // Create new restriction 243 CeedCall(CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp * num_comp, 1, l_size * num_comp, CEED_MEM_HOST, CEED_OWN_POINTER, 244 pointblock_offsets, pointblock_rstr)); 245 246 // Cleanup 247 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 248 249 return CEED_ERROR_SUCCESS; 250 } 251 252 /** 253 @brief Core logic for assembling operator diagonal or point block diagonal 254 255 @param[in] op CeedOperator to assemble point block diagonal 256 @param[in] request Address of CeedRequest for non-blocking completion, else CEED_REQUEST_IMMEDIATE 257 @param[in] is_pointblock Boolean flag to assemble diagonal or point block diagonal 258 @param[out] assembled CeedVector to store assembled diagonal 259 260 @return An error code: 0 - success, otherwise - failure 261 262 @ref Developer 263 **/ 264 static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, CeedRequest *request, const bool is_pointblock, CeedVector assembled) { 265 Ceed ceed; 266 CeedCall(CeedOperatorGetCeed(op, &ceed)); 267 268 // Assemble QFunction 269 CeedQFunction qf; 270 const CeedScalar *assembled_qf_array; 271 CeedVector assembled_qf; 272 CeedElemRestriction assembled_elem_rstr; 273 CeedInt num_input_fields, num_output_fields; 274 CeedInt layout[3]; 275 276 CeedCall(CeedOperatorGetQFunction(op, &qf)); 277 CeedCall(CeedQFunctionGetNumArgs(qf, &num_input_fields, &num_output_fields)); 278 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, request)); 279 CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, &layout)); 280 CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr)); 281 CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); 282 283 // Get assembly data 284 CeedOperatorAssemblyData data; 285 const CeedEvalMode **eval_modes_in, **eval_modes_out; 286 CeedInt *num_eval_modes_in, *num_eval_modes_out, num_active_bases; 287 CeedSize **eval_mode_offsets_in, **eval_mode_offsets_out, num_output_components; 288 CeedBasis *active_bases; 289 CeedElemRestriction *active_elem_rstrs; 290 291 CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); 292 CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases, &num_eval_modes_in, &eval_modes_in, &eval_mode_offsets_in, 293 &num_eval_modes_out, &eval_modes_out, &eval_mode_offsets_out, &num_output_components)); 294 CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases, NULL, NULL)); 295 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, NULL, &active_elem_rstrs)); 296 297 // Loop over all active bases 298 for (CeedInt b = 0; b < num_active_bases; b++) { 299 // Assemble point block diagonal restriction, if needed 300 CeedElemRestriction diag_elem_rstr = active_elem_rstrs[b]; 301 302 if (is_pointblock) { 303 CeedElemRestriction point_block_elem_rstr; 304 305 CeedCall(CeedOperatorCreateActivePointBlockRestriction(diag_elem_rstr, &point_block_elem_rstr)); 306 diag_elem_rstr = point_block_elem_rstr; 307 } 308 309 // Create diagonal vector 310 CeedVector elem_diag; 311 CeedCall(CeedElemRestrictionCreateVector(diag_elem_rstr, NULL, &elem_diag)); 312 313 // Assemble element operator diagonals 314 CeedScalar *elem_diag_array; 315 CeedInt num_elem, num_nodes, num_qpts, num_components; 316 317 CeedCall(CeedVectorSetValue(elem_diag, 0.0)); 318 CeedCall(CeedVectorGetArray(elem_diag, CEED_MEM_HOST, &elem_diag_array)); 319 CeedCall(CeedElemRestrictionGetNumElements(diag_elem_rstr, &num_elem)); 320 CeedCall(CeedBasisGetNumNodes(active_bases[b], &num_nodes)); 321 CeedCall(CeedBasisGetNumComponents(active_bases[b], &num_components)); 322 CeedCall(CeedBasisGetNumQuadraturePoints(active_bases[b], &num_qpts)); 323 324 // Basis matrices 325 const CeedScalar *interp, *grad; 326 CeedScalar *identity = NULL; 327 bool has_eval_none = false; 328 for (CeedInt i = 0; i < num_eval_modes_in[b]; i++) { 329 has_eval_none = has_eval_none || (eval_modes_in[b][i] == CEED_EVAL_NONE); 330 } 331 for (CeedInt i = 0; i < num_eval_modes_out[b]; i++) { 332 has_eval_none = has_eval_none || (eval_modes_out[b][i] == CEED_EVAL_NONE); 333 } 334 if (has_eval_none) { 335 CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); 336 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; 337 } 338 CeedCall(CeedBasisGetInterp(active_bases[b], &interp)); 339 CeedCall(CeedBasisGetGrad(active_bases[b], &grad)); 340 // Compute the diagonal of B^T D B 341 // Each element 342 for (CeedInt e = 0; e < num_elem; e++) { 343 CeedInt d_out = -1; 344 // Each basis eval mode pair 345 for (CeedInt e_out = 0; e_out < num_eval_modes_out[b]; e_out++) { 346 const CeedScalar *B_t = NULL; 347 if (eval_modes_out[b][e_out] == CEED_EVAL_GRAD) d_out += 1; 348 CeedOperatorGetBasisPointer(eval_modes_out[b][e_out], identity, interp, &grad[d_out * num_qpts * num_nodes], &B_t); 349 CeedInt d_in = -1; 350 for (CeedInt e_in = 0; e_in < num_eval_modes_in[b]; e_in++) { 351 const CeedScalar *B = NULL; 352 if (eval_modes_in[b][e_in] == CEED_EVAL_GRAD) d_in += 1; 353 CeedOperatorGetBasisPointer(eval_modes_in[b][e_in], identity, interp, &grad[d_in * num_qpts * num_nodes], &B); 354 // Each component 355 for (CeedInt c_out = 0; c_out < num_components; c_out++) { 356 // Each qpt/node pair 357 for (CeedInt q = 0; q < num_qpts; q++) { 358 if (is_pointblock) { 359 // Point Block Diagonal 360 for (CeedInt c_in = 0; c_in < num_components; c_in++) { 361 const CeedInt c_offset = (eval_mode_offsets_in[b][e_in] + c_in) * num_output_components + eval_mode_offsets_out[b][e_out] + c_out; 362 const CeedScalar qf_value = assembled_qf_array[q * layout[0] + c_offset * layout[1] + e * layout[2]]; 363 for (CeedInt n = 0; n < num_nodes; n++) { 364 elem_diag_array[((e * num_components + c_out) * num_components + c_in) * num_nodes + n] += 365 B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n]; 366 } 367 } 368 } else { 369 // Diagonal Only 370 const CeedInt c_offset = (eval_mode_offsets_in[b][e_in] + c_out) * num_output_components + eval_mode_offsets_out[b][e_out] + c_out; 371 const CeedScalar qf_value = assembled_qf_array[q * layout[0] + c_offset * layout[1] + e * layout[2]]; 372 for (CeedInt n = 0; n < num_nodes; n++) { 373 elem_diag_array[(e * num_components + c_out) * num_nodes + n] += B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n]; 374 } 375 } 376 } 377 } 378 } 379 } 380 } 381 CeedCall(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); 382 383 // Assemble local operator diagonal 384 CeedCall(CeedElemRestrictionApply(diag_elem_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); 385 386 // Cleanup 387 if (is_pointblock) CeedCall(CeedElemRestrictionDestroy(&diag_elem_rstr)); 388 CeedCall(CeedVectorDestroy(&elem_diag)); 389 CeedCall(CeedFree(&identity)); 390 } 391 CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 392 CeedCall(CeedVectorDestroy(&assembled_qf)); 393 394 return CEED_ERROR_SUCCESS; 395 } 396 397 /** 398 @brief Core logic for assembling composite operator diagonal 399 400 @param[in] op CeedOperator to assemble point block diagonal 401 @param[in] request Address of CeedRequest for non-blocking completion, else CEED_REQUEST_IMMEDIATE 402 @param[in] is_pointblock Boolean flag to assemble diagonal or point block diagonal 403 @param[out] assembled CeedVector to store assembled diagonal 404 405 @return An error code: 0 - success, otherwise - failure 406 407 @ref Developer 408 **/ 409 static inline int CeedCompositeOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedRequest *request, const bool is_pointblock, 410 CeedVector assembled) { 411 CeedInt num_sub; 412 CeedOperator *suboperators; 413 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); 414 CeedCall(CeedCompositeOperatorGetSubList(op, &suboperators)); 415 for (CeedInt i = 0; i < num_sub; i++) { 416 if (is_pointblock) { 417 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(suboperators[i], assembled, request)); 418 } else { 419 CeedCall(CeedOperatorLinearAssembleAddDiagonal(suboperators[i], assembled, request)); 420 } 421 } 422 return CEED_ERROR_SUCCESS; 423 } 424 425 /** 426 @brief Build nonzero pattern for non-composite operator 427 428 Users should generally use CeedOperatorLinearAssembleSymbolic() 429 430 @param[in] op CeedOperator to assemble nonzero pattern 431 @param[in] offset Offset for number of entries 432 @param[out] rows Row number for each entry 433 @param[out] cols Column number for each entry 434 435 @return An error code: 0 - success, otherwise - failure 436 437 @ref Developer 438 **/ 439 static int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, CeedInt *rows, CeedInt *cols) { 440 Ceed ceed; 441 bool is_composite; 442 CeedCall(CeedOperatorGetCeed(op, &ceed)); 443 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 444 445 if (is_composite) { 446 // LCOV_EXCL_START 447 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 448 // LCOV_EXCL_STOP 449 } 450 451 CeedSize num_nodes; 452 CeedCall(CeedOperatorGetActiveVectorLengths(op, &num_nodes, NULL)); 453 CeedElemRestriction rstr_in; 454 CeedCall(CeedOperatorGetActiveElemRestriction(op, &rstr_in)); 455 CeedInt num_elem, elem_size, num_comp; 456 CeedCall(CeedElemRestrictionGetNumElements(rstr_in, &num_elem)); 457 CeedCall(CeedElemRestrictionGetElementSize(rstr_in, &elem_size)); 458 CeedCall(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp)); 459 CeedInt layout_er[3]; 460 CeedCall(CeedElemRestrictionGetELayout(rstr_in, &layout_er)); 461 462 CeedInt local_num_entries = elem_size * num_comp * elem_size * num_comp * num_elem; 463 464 // Determine elem_dof relation 465 CeedVector index_vec; 466 CeedCall(CeedVectorCreate(ceed, num_nodes, &index_vec)); 467 CeedScalar *array; 468 CeedCall(CeedVectorGetArrayWrite(index_vec, CEED_MEM_HOST, &array)); 469 for (CeedInt i = 0; i < num_nodes; i++) array[i] = i; 470 CeedCall(CeedVectorRestoreArray(index_vec, &array)); 471 CeedVector elem_dof; 472 CeedCall(CeedVectorCreate(ceed, num_elem * elem_size * num_comp, &elem_dof)); 473 CeedCall(CeedVectorSetValue(elem_dof, 0.0)); 474 CeedCall(CeedElemRestrictionApply(rstr_in, CEED_NOTRANSPOSE, index_vec, elem_dof, CEED_REQUEST_IMMEDIATE)); 475 const CeedScalar *elem_dof_a; 476 CeedCall(CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a)); 477 CeedCall(CeedVectorDestroy(&index_vec)); 478 479 // Determine i, j locations for element matrices 480 CeedInt count = 0; 481 for (CeedInt e = 0; e < num_elem; e++) { 482 for (CeedInt comp_in = 0; comp_in < num_comp; comp_in++) { 483 for (CeedInt comp_out = 0; comp_out < num_comp; comp_out++) { 484 for (CeedInt i = 0; i < elem_size; i++) { 485 for (CeedInt j = 0; j < elem_size; j++) { 486 const CeedInt elem_dof_index_row = i * layout_er[0] + (comp_out)*layout_er[1] + e * layout_er[2]; 487 const CeedInt elem_dof_index_col = j * layout_er[0] + comp_in * layout_er[1] + e * layout_er[2]; 488 489 const CeedInt row = elem_dof_a[elem_dof_index_row]; 490 const CeedInt col = elem_dof_a[elem_dof_index_col]; 491 492 rows[offset + count] = row; 493 cols[offset + count] = col; 494 count++; 495 } 496 } 497 } 498 } 499 } 500 if (count != local_num_entries) { 501 // LCOV_EXCL_START 502 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing assembled entries"); 503 // LCOV_EXCL_STOP 504 } 505 CeedCall(CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a)); 506 CeedCall(CeedVectorDestroy(&elem_dof)); 507 508 return CEED_ERROR_SUCCESS; 509 } 510 511 /** 512 @brief Assemble nonzero entries for non-composite operator 513 514 Users should generally use CeedOperatorLinearAssemble() 515 516 @param[in] op CeedOperator to assemble 517 @param[in] offset Offset for number of entries 518 @param[out] values Values to assemble into matrix 519 520 @return An error code: 0 - success, otherwise - failure 521 522 @ref Developer 523 **/ 524 static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVector values) { 525 Ceed ceed; 526 bool is_composite; 527 CeedCall(CeedOperatorGetCeed(op, &ceed)); 528 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 529 530 if (is_composite) { 531 // LCOV_EXCL_START 532 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 533 // LCOV_EXCL_STOP 534 } 535 536 // Early exit for empty operator 537 { 538 CeedInt num_elem = 0; 539 540 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 541 if (num_elem == 0) return CEED_ERROR_SUCCESS; 542 } 543 544 if (op->LinearAssembleSingle) { 545 // Backend version 546 CeedCall(op->LinearAssembleSingle(op, offset, values)); 547 return CEED_ERROR_SUCCESS; 548 } else { 549 // Operator fallback 550 CeedOperator op_fallback; 551 552 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 553 if (op_fallback) { 554 CeedCall(CeedSingleOperatorAssemble(op_fallback, offset, values)); 555 return CEED_ERROR_SUCCESS; 556 } 557 } 558 559 // Assemble QFunction 560 CeedQFunction qf; 561 CeedCall(CeedOperatorGetQFunction(op, &qf)); 562 CeedVector assembled_qf; 563 CeedElemRestriction rstr_q; 564 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &rstr_q, CEED_REQUEST_IMMEDIATE)); 565 CeedSize qf_length; 566 CeedCall(CeedVectorGetLength(assembled_qf, &qf_length)); 567 568 CeedInt num_input_fields, num_output_fields; 569 CeedOperatorField *input_fields; 570 CeedOperatorField *output_fields; 571 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 572 573 // Get assembly data 574 CeedOperatorAssemblyData data; 575 CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); 576 const CeedEvalMode **eval_modes_in, **eval_modes_out; 577 CeedInt *num_eval_modes_in, *num_eval_modes_out, num_active_bases; 578 CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases, &num_eval_modes_in, &eval_modes_in, NULL, &num_eval_modes_out, 579 &eval_modes_out, NULL, NULL)); 580 CeedBasis *bases; 581 CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &bases, NULL, NULL)); 582 CeedBasis basis_in = bases[0]; 583 584 if (num_active_bases > 1) { 585 // LCOV_EXCL_START 586 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator with multiple active bases"); 587 // LCOV_EXCL_STOP 588 } 589 if (num_eval_modes_in[0] == 0 || num_eval_modes_out[0] == 0) { 590 // LCOV_EXCL_START 591 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator with out inputs/outputs"); 592 // LCOV_EXCL_STOP 593 } 594 595 CeedElemRestriction active_rstr; 596 CeedInt num_elem, elem_size, num_qpts, num_comp; 597 CeedCall(CeedOperatorGetActiveElemRestriction(op, &active_rstr)); 598 CeedCall(CeedElemRestrictionGetNumElements(active_rstr, &num_elem)); 599 CeedCall(CeedElemRestrictionGetElementSize(active_rstr, &elem_size)); 600 CeedCall(CeedElemRestrictionGetNumComponents(active_rstr, &num_comp)); 601 CeedCall(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts)); 602 603 CeedInt local_num_entries = elem_size * num_comp * elem_size * num_comp * num_elem; 604 605 // loop over elements and put in data structure 606 const CeedScalar *interp_in, *grad_in; 607 CeedCall(CeedBasisGetInterp(basis_in, &interp_in)); 608 CeedCall(CeedBasisGetGrad(basis_in, &grad_in)); 609 610 const CeedScalar *assembled_qf_array; 611 CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); 612 613 CeedInt layout_qf[3]; 614 CeedCall(CeedElemRestrictionGetELayout(rstr_q, &layout_qf)); 615 CeedCall(CeedElemRestrictionDestroy(&rstr_q)); 616 617 // we store B_mat_in, B_mat_out, BTD, elem_mat in row-major order 618 const CeedScalar **B_mats_in, **B_mats_out; 619 CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, NULL, &B_mats_in, &B_mats_out)); 620 const CeedScalar *B_mat_in = B_mats_in[0], *B_mat_out = B_mats_out[0]; 621 CeedScalar BTD_mat[elem_size * num_qpts * num_eval_modes_in[0]]; 622 CeedScalar elem_mat[elem_size * elem_size]; 623 CeedInt count = 0; 624 CeedScalar *vals; 625 CeedCall(CeedVectorGetArrayWrite(values, CEED_MEM_HOST, &vals)); 626 for (CeedInt e = 0; e < num_elem; e++) { 627 for (CeedInt comp_in = 0; comp_in < num_comp; comp_in++) { 628 for (CeedInt comp_out = 0; comp_out < num_comp; comp_out++) { 629 // Compute B^T*D 630 for (CeedInt n = 0; n < elem_size; n++) { 631 for (CeedInt q = 0; q < num_qpts; q++) { 632 for (CeedInt e_in = 0; e_in < num_eval_modes_in[0]; e_in++) { 633 const CeedInt btd_index = n * (num_qpts * num_eval_modes_in[0]) + (num_eval_modes_in[0] * q + e_in); 634 CeedScalar sum = 0.0; 635 for (CeedInt e_out = 0; e_out < num_eval_modes_out[0]; e_out++) { 636 const CeedInt b_out_index = (num_eval_modes_out[0] * q + e_out) * elem_size + n; 637 const CeedInt eval_mode_index = ((e_in * num_comp + comp_in) * num_eval_modes_out[0] + e_out) * num_comp + comp_out; 638 const CeedInt qf_index = q * layout_qf[0] + eval_mode_index * layout_qf[1] + e * layout_qf[2]; 639 sum += B_mat_out[b_out_index] * assembled_qf_array[qf_index]; 640 } 641 BTD_mat[btd_index] = sum; 642 } 643 } 644 } 645 // form element matrix itself (for each block component) 646 CeedCall(CeedMatrixMatrixMultiply(ceed, BTD_mat, B_mat_in, elem_mat, elem_size, elem_size, num_qpts * num_eval_modes_in[0])); 647 648 // put element matrix in coordinate data structure 649 for (CeedInt i = 0; i < elem_size; i++) { 650 for (CeedInt j = 0; j < elem_size; j++) { 651 vals[offset + count] = elem_mat[i * elem_size + j]; 652 count++; 653 } 654 } 655 } 656 } 657 } 658 if (count != local_num_entries) { 659 // LCOV_EXCL_START 660 return CeedError(ceed, CEED_ERROR_MAJOR, "Error computing entries"); 661 // LCOV_EXCL_STOP 662 } 663 CeedCall(CeedVectorRestoreArray(values, &vals)); 664 665 CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 666 CeedCall(CeedVectorDestroy(&assembled_qf)); 667 668 return CEED_ERROR_SUCCESS; 669 } 670 671 /** 672 @brief Count number of entries for assembled CeedOperator 673 674 @param[in] op CeedOperator to assemble 675 @param[out] num_entries Number of entries in assembled representation 676 677 @return An error code: 0 - success, otherwise - failure 678 679 @ref Utility 680 **/ 681 static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedInt *num_entries) { 682 bool is_composite; 683 CeedElemRestriction rstr; 684 CeedInt num_elem, elem_size, num_comp; 685 686 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 687 if (is_composite) { 688 // LCOV_EXCL_START 689 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 690 // LCOV_EXCL_STOP 691 } 692 CeedCall(CeedOperatorGetActiveElemRestriction(op, &rstr)); 693 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 694 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 695 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 696 *num_entries = elem_size * num_comp * elem_size * num_comp * num_elem; 697 698 return CEED_ERROR_SUCCESS; 699 } 700 701 /** 702 @brief Common code for creating a multigrid coarse operator and level transfer operators for a CeedOperator 703 704 @param[in] op_fine Fine grid operator 705 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 706 @param[in] rstr_coarse Coarse grid restriction 707 @param[in] basis_coarse Coarse grid active vector basis 708 @param[in] basis_c_to_f Basis for coarse to fine interpolation 709 @param[out] op_coarse Coarse grid operator 710 @param[out] op_prolong Coarse to fine operator 711 @param[out] op_restrict Fine to coarse operator 712 713 @return An error code: 0 - success, otherwise - failure 714 715 @ref Developer 716 **/ 717 static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 718 CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { 719 Ceed ceed; 720 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 721 722 // Check for composite operator 723 bool is_composite; 724 CeedCall(CeedOperatorIsComposite(op_fine, &is_composite)); 725 if (is_composite) { 726 // LCOV_EXCL_START 727 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Automatic multigrid setup for composite operators not supported"); 728 // LCOV_EXCL_STOP 729 } 730 731 // Coarse Grid 732 CeedCall(CeedOperatorCreate(ceed, op_fine->qf, op_fine->dqf, op_fine->dqfT, op_coarse)); 733 CeedElemRestriction rstr_fine = NULL; 734 // -- Clone input fields 735 for (CeedInt i = 0; i < op_fine->qf->num_input_fields; i++) { 736 if (op_fine->input_fields[i]->vec == CEED_VECTOR_ACTIVE) { 737 rstr_fine = op_fine->input_fields[i]->elem_rstr; 738 CeedCall(CeedOperatorSetField(*op_coarse, op_fine->input_fields[i]->field_name, rstr_coarse, basis_coarse, CEED_VECTOR_ACTIVE)); 739 } else { 740 CeedCall(CeedOperatorSetField(*op_coarse, op_fine->input_fields[i]->field_name, op_fine->input_fields[i]->elem_rstr, 741 op_fine->input_fields[i]->basis, op_fine->input_fields[i]->vec)); 742 } 743 } 744 // -- Clone output fields 745 for (CeedInt i = 0; i < op_fine->qf->num_output_fields; i++) { 746 if (op_fine->output_fields[i]->vec == CEED_VECTOR_ACTIVE) { 747 CeedCall(CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name, rstr_coarse, basis_coarse, CEED_VECTOR_ACTIVE)); 748 } else { 749 CeedCall(CeedOperatorSetField(*op_coarse, op_fine->output_fields[i]->field_name, op_fine->output_fields[i]->elem_rstr, 750 op_fine->output_fields[i]->basis, op_fine->output_fields[i]->vec)); 751 } 752 } 753 // -- Clone QFunctionAssemblyData 754 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(op_fine->qf_assembled, &(*op_coarse)->qf_assembled)); 755 756 // Multiplicity vector 757 CeedVector mult_vec, mult_e_vec; 758 CeedCall(CeedElemRestrictionCreateVector(rstr_fine, &mult_vec, &mult_e_vec)); 759 CeedCall(CeedVectorSetValue(mult_e_vec, 0.0)); 760 CeedCall(CeedElemRestrictionApply(rstr_fine, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE)); 761 CeedCall(CeedVectorSetValue(mult_vec, 0.0)); 762 CeedCall(CeedElemRestrictionApply(rstr_fine, CEED_TRANSPOSE, mult_e_vec, mult_vec, CEED_REQUEST_IMMEDIATE)); 763 CeedCall(CeedVectorDestroy(&mult_e_vec)); 764 CeedCall(CeedVectorReciprocal(mult_vec)); 765 766 // Restriction 767 CeedInt num_comp; 768 CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp)); 769 CeedQFunction qf_restrict; 770 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_restrict)); 771 CeedInt *num_comp_r_data; 772 CeedCall(CeedCalloc(1, &num_comp_r_data)); 773 num_comp_r_data[0] = num_comp; 774 CeedQFunctionContext ctx_r; 775 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r)); 776 CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data)); 777 CeedCall(CeedQFunctionSetContext(qf_restrict, ctx_r)); 778 CeedCall(CeedQFunctionContextDestroy(&ctx_r)); 779 CeedCall(CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE)); 780 CeedCall(CeedQFunctionAddInput(qf_restrict, "scale", num_comp, CEED_EVAL_NONE)); 781 CeedCall(CeedQFunctionAddOutput(qf_restrict, "output", num_comp, CEED_EVAL_INTERP)); 782 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_restrict, num_comp)); 783 784 CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict)); 785 CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE)); 786 CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_fine, CEED_BASIS_COLLOCATED, mult_vec)); 787 CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); 788 789 // Prolongation 790 CeedQFunction qf_prolong; 791 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_prolong)); 792 CeedInt *num_comp_p_data; 793 CeedCall(CeedCalloc(1, &num_comp_p_data)); 794 num_comp_p_data[0] = num_comp; 795 CeedQFunctionContext ctx_p; 796 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_p)); 797 CeedCall(CeedQFunctionContextSetData(ctx_p, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_p_data), num_comp_p_data)); 798 CeedCall(CeedQFunctionSetContext(qf_prolong, ctx_p)); 799 CeedCall(CeedQFunctionContextDestroy(&ctx_p)); 800 CeedCall(CeedQFunctionAddInput(qf_prolong, "input", num_comp, CEED_EVAL_INTERP)); 801 CeedCall(CeedQFunctionAddInput(qf_prolong, "scale", num_comp, CEED_EVAL_NONE)); 802 CeedCall(CeedQFunctionAddOutput(qf_prolong, "output", num_comp, CEED_EVAL_NONE)); 803 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_prolong, num_comp)); 804 805 CeedCall(CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_prolong)); 806 CeedCall(CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); 807 CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_fine, CEED_BASIS_COLLOCATED, mult_vec)); 808 CeedCall(CeedOperatorSetField(*op_prolong, "output", rstr_fine, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE)); 809 810 // Clone name 811 bool has_name = op_fine->name; 812 size_t name_len = op_fine->name ? strlen(op_fine->name) : 0; 813 CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name)); 814 { 815 char *prolongation_name; 816 CeedCall(CeedCalloc(18 + name_len, &prolongation_name)); 817 sprintf(prolongation_name, "prolongation%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); 818 CeedCall(CeedOperatorSetName(*op_prolong, prolongation_name)); 819 CeedCall(CeedFree(&prolongation_name)); 820 } 821 { 822 char *restriction_name; 823 CeedCall(CeedCalloc(17 + name_len, &restriction_name)); 824 sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); 825 CeedCall(CeedOperatorSetName(*op_restrict, restriction_name)); 826 CeedCall(CeedFree(&restriction_name)); 827 } 828 829 // Check 830 CeedCall(CeedOperatorCheckReady(*op_coarse)); 831 CeedCall(CeedOperatorCheckReady(*op_prolong)); 832 CeedCall(CeedOperatorCheckReady(*op_restrict)); 833 834 // Cleanup 835 CeedCall(CeedVectorDestroy(&mult_vec)); 836 CeedCall(CeedBasisDestroy(&basis_c_to_f)); 837 CeedCall(CeedQFunctionDestroy(&qf_restrict)); 838 CeedCall(CeedQFunctionDestroy(&qf_prolong)); 839 840 return CEED_ERROR_SUCCESS; 841 } 842 843 /** 844 @brief Build 1D mass matrix and Laplacian with perturbation 845 846 @param[in] interp_1d Interpolation matrix in one dimension 847 @param[in] grad_1d Gradient matrix in one dimension 848 @param[in] q_weight_1d Quadrature weights in one dimension 849 @param[in] P_1d Number of basis nodes in one dimension 850 @param[in] Q_1d Number of quadrature points in one dimension 851 @param[in] dim Dimension of basis 852 @param[out] mass Assembled mass matrix in one dimension 853 @param[out] laplace Assembled perturbed Laplacian in one dimension 854 855 @return An error code: 0 - success, otherwise - failure 856 857 @ref Developer 858 **/ 859 CeedPragmaOptimizeOff static int CeedBuildMassLaplace(const CeedScalar *interp_1d, const CeedScalar *grad_1d, const CeedScalar *q_weight_1d, 860 CeedInt P_1d, CeedInt Q_1d, CeedInt dim, CeedScalar *mass, CeedScalar *laplace) { 861 for (CeedInt i = 0; i < P_1d; i++) { 862 for (CeedInt j = 0; j < P_1d; j++) { 863 CeedScalar sum = 0.0; 864 for (CeedInt k = 0; k < Q_1d; k++) sum += interp_1d[k * P_1d + i] * q_weight_1d[k] * interp_1d[k * P_1d + j]; 865 mass[i + j * P_1d] = sum; 866 } 867 } 868 // -- Laplacian 869 for (CeedInt i = 0; i < P_1d; i++) { 870 for (CeedInt j = 0; j < P_1d; j++) { 871 CeedScalar sum = 0.0; 872 for (CeedInt k = 0; k < Q_1d; k++) sum += grad_1d[k * P_1d + i] * q_weight_1d[k] * grad_1d[k * P_1d + j]; 873 laplace[i + j * P_1d] = sum; 874 } 875 } 876 CeedScalar perturbation = dim > 2 ? 1e-6 : 1e-4; 877 for (CeedInt i = 0; i < P_1d; i++) laplace[i + P_1d * i] += perturbation; 878 return CEED_ERROR_SUCCESS; 879 } 880 CeedPragmaOptimizeOn; 881 882 /// @} 883 884 /// ---------------------------------------------------------------------------- 885 /// CeedOperator Backend API 886 /// ---------------------------------------------------------------------------- 887 /// @addtogroup CeedOperatorBackend 888 /// @{ 889 890 /** 891 @brief Create object holding CeedQFunction assembly data for CeedOperator 892 893 @param[in] ceed A Ceed object where the CeedQFunctionAssemblyData will be created 894 @param[out] data Address of the variable where the newly created CeedQFunctionAssemblyData will be stored 895 896 @return An error code: 0 - success, otherwise - failure 897 898 @ref Backend 899 **/ 900 int CeedQFunctionAssemblyDataCreate(Ceed ceed, CeedQFunctionAssemblyData *data) { 901 CeedCall(CeedCalloc(1, data)); 902 (*data)->ref_count = 1; 903 (*data)->ceed = ceed; 904 CeedCall(CeedReference(ceed)); 905 906 return CEED_ERROR_SUCCESS; 907 } 908 909 /** 910 @brief Increment the reference counter for a CeedQFunctionAssemblyData 911 912 @param[in,out] data CeedQFunctionAssemblyData to increment the reference counter 913 914 @return An error code: 0 - success, otherwise - failure 915 916 @ref Backend 917 **/ 918 int CeedQFunctionAssemblyDataReference(CeedQFunctionAssemblyData data) { 919 data->ref_count++; 920 return CEED_ERROR_SUCCESS; 921 } 922 923 /** 924 @brief Set re-use of CeedQFunctionAssemblyData 925 926 @param[in,out] data CeedQFunctionAssemblyData to mark for reuse 927 @param[in] reuse_data Boolean flag indicating data re-use 928 929 @return An error code: 0 - success, otherwise - failure 930 931 @ref Backend 932 **/ 933 int CeedQFunctionAssemblyDataSetReuse(CeedQFunctionAssemblyData data, bool reuse_data) { 934 data->reuse_data = reuse_data; 935 data->needs_data_update = true; 936 return CEED_ERROR_SUCCESS; 937 } 938 939 /** 940 @brief Mark QFunctionAssemblyData as stale 941 942 @param[in,out] data CeedQFunctionAssemblyData to mark as stale 943 @param[in] needs_data_update Boolean flag indicating if update is needed or completed 944 945 @return An error code: 0 - success, otherwise - failure 946 947 @ref Backend 948 **/ 949 int CeedQFunctionAssemblyDataSetUpdateNeeded(CeedQFunctionAssemblyData data, bool needs_data_update) { 950 data->needs_data_update = needs_data_update; 951 return CEED_ERROR_SUCCESS; 952 } 953 954 /** 955 @brief Determine if QFunctionAssemblyData needs update 956 957 @param[in] data CeedQFunctionAssemblyData to mark as stale 958 @param[out] is_update_needed Boolean flag indicating if re-assembly is required 959 960 @return An error code: 0 - success, otherwise - failure 961 962 @ref Backend 963 **/ 964 int CeedQFunctionAssemblyDataIsUpdateNeeded(CeedQFunctionAssemblyData data, bool *is_update_needed) { 965 *is_update_needed = !data->reuse_data || data->needs_data_update; 966 return CEED_ERROR_SUCCESS; 967 } 968 969 /** 970 @brief Copy the pointer to a CeedQFunctionAssemblyData. 971 Both pointers should be destroyed with `CeedCeedQFunctionAssemblyDataDestroy()`. 972 Note: If `*data_copy` is non-NULL, then it is assumed that `*data_copy` is a pointer to a CeedQFunctionAssemblyData. 973 This CeedQFunctionAssemblyData will be destroyed if `*data_copy` is the only reference to this CeedQFunctionAssemblyData. 974 975 @param[in] data CeedQFunctionAssemblyData to copy reference to 976 @param[in,out] data_copy Variable to store copied reference 977 978 @return An error code: 0 - success, otherwise - failure 979 980 @ref Backend 981 **/ 982 int CeedQFunctionAssemblyDataReferenceCopy(CeedQFunctionAssemblyData data, CeedQFunctionAssemblyData *data_copy) { 983 CeedCall(CeedQFunctionAssemblyDataReference(data)); 984 CeedCall(CeedQFunctionAssemblyDataDestroy(data_copy)); 985 *data_copy = data; 986 return CEED_ERROR_SUCCESS; 987 } 988 989 /** 990 @brief Get setup status for internal objects for CeedQFunctionAssemblyData 991 992 @param[in] data CeedQFunctionAssemblyData to retrieve status 993 @param[out] is_setup Boolean flag for setup status 994 995 @return An error code: 0 - success, otherwise - failure 996 997 @ref Backend 998 **/ 999 int CeedQFunctionAssemblyDataIsSetup(CeedQFunctionAssemblyData data, bool *is_setup) { 1000 *is_setup = data->is_setup; 1001 return CEED_ERROR_SUCCESS; 1002 } 1003 1004 /** 1005 @brief Set internal objects for CeedQFunctionAssemblyData 1006 1007 @param[in,out] data CeedQFunctionAssemblyData to set objects 1008 @param[in] vec CeedVector to store assembled CeedQFunction at quadrature points 1009 @param[in] rstr CeedElemRestriction for CeedVector containing assembled CeedQFunction 1010 1011 @return An error code: 0 - success, otherwise - failure 1012 1013 @ref Backend 1014 **/ 1015 int CeedQFunctionAssemblyDataSetObjects(CeedQFunctionAssemblyData data, CeedVector vec, CeedElemRestriction rstr) { 1016 CeedCall(CeedVectorReferenceCopy(vec, &data->vec)); 1017 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &data->rstr)); 1018 1019 data->is_setup = true; 1020 return CEED_ERROR_SUCCESS; 1021 } 1022 1023 int CeedQFunctionAssemblyDataGetObjects(CeedQFunctionAssemblyData data, CeedVector *vec, CeedElemRestriction *rstr) { 1024 if (!data->is_setup) { 1025 // LCOV_EXCL_START 1026 return CeedError(data->ceed, CEED_ERROR_INCOMPLETE, "Internal objects not set; must call CeedQFunctionAssemblyDataSetObjects first."); 1027 // LCOV_EXCL_STOP 1028 } 1029 1030 CeedCall(CeedVectorReferenceCopy(data->vec, vec)); 1031 CeedCall(CeedElemRestrictionReferenceCopy(data->rstr, rstr)); 1032 1033 return CEED_ERROR_SUCCESS; 1034 } 1035 1036 /** 1037 @brief Destroy CeedQFunctionAssemblyData 1038 1039 @param[in,out] data CeedQFunctionAssemblyData to destroy 1040 1041 @return An error code: 0 - success, otherwise - failure 1042 1043 @ref Backend 1044 **/ 1045 int CeedQFunctionAssemblyDataDestroy(CeedQFunctionAssemblyData *data) { 1046 if (!*data || --(*data)->ref_count > 0) return CEED_ERROR_SUCCESS; 1047 1048 CeedCall(CeedDestroy(&(*data)->ceed)); 1049 CeedCall(CeedVectorDestroy(&(*data)->vec)); 1050 CeedCall(CeedElemRestrictionDestroy(&(*data)->rstr)); 1051 1052 CeedCall(CeedFree(data)); 1053 return CEED_ERROR_SUCCESS; 1054 } 1055 1056 /** 1057 @brief Get CeedOperatorAssemblyData 1058 1059 @param[in] op CeedOperator to assemble 1060 @param[out] data CeedQFunctionAssemblyData 1061 1062 @return An error code: 0 - success, otherwise - failure 1063 1064 @ref Backend 1065 **/ 1066 int CeedOperatorGetOperatorAssemblyData(CeedOperator op, CeedOperatorAssemblyData *data) { 1067 if (!op->op_assembled) { 1068 CeedOperatorAssemblyData data; 1069 1070 CeedCall(CeedOperatorAssemblyDataCreate(op->ceed, op, &data)); 1071 op->op_assembled = data; 1072 } 1073 *data = op->op_assembled; 1074 1075 return CEED_ERROR_SUCCESS; 1076 } 1077 1078 /** 1079 @brief Create object holding CeedOperator assembly data 1080 1081 @param[in] ceed Ceed object where the CeedOperatorAssemblyData will be created 1082 @param[in] op CeedOperator to be assembled 1083 @param[out] data Address of the variable where the newly created CeedOperatorAssemblyData will be stored 1084 1085 @return An error code: 0 - success, otherwise - failure 1086 1087 @ref Backend 1088 **/ 1089 int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) { 1090 CeedInt num_active_bases = 0; 1091 1092 // Allocate 1093 CeedCall(CeedCalloc(1, data)); 1094 (*data)->ceed = ceed; 1095 CeedCall(CeedReference(ceed)); 1096 1097 // Build OperatorAssembly data 1098 CeedQFunction qf; 1099 CeedQFunctionField *qf_fields; 1100 CeedOperatorField *op_fields; 1101 CeedInt num_input_fields; 1102 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1103 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL)); 1104 CeedCall(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 1105 1106 // Determine active input basis 1107 CeedInt *num_eval_modes_in = NULL, *num_eval_modes_out = NULL, offset = 0; 1108 CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL; 1109 CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL; 1110 1111 for (CeedInt i = 0; i < num_input_fields; i++) { 1112 CeedVector vec; 1113 1114 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1115 if (vec == CEED_VECTOR_ACTIVE) { 1116 CeedInt index = -1, dim = 1, num_components; 1117 CeedBasis basis_in = NULL; 1118 CeedEvalMode eval_mode; 1119 1120 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 1121 CeedCall(CeedBasisGetDimension(basis_in, &dim)); 1122 CeedCall(CeedBasisGetNumComponents(basis_in, &num_components)); 1123 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1124 for (CeedInt i = 0; i < num_active_bases; i++) { 1125 if ((*data)->active_bases[i] == basis_in) index = i; 1126 } 1127 if (index == -1) { 1128 CeedElemRestriction elem_rstr_in; 1129 1130 index = num_active_bases; 1131 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->active_bases)); 1132 (*data)->active_bases[num_active_bases] = NULL; 1133 CeedCall(CeedBasisReferenceCopy(basis_in, &(*data)->active_bases[num_active_bases])); 1134 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->active_elem_rstrs)); 1135 (*data)->active_elem_rstrs[num_active_bases] = NULL; 1136 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_in)); 1137 CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_in, &(*data)->active_elem_rstrs[num_active_bases])); 1138 CeedCall(CeedRealloc(num_active_bases + 1, &num_eval_modes_in)); 1139 CeedCall(CeedRealloc(num_active_bases + 1, &num_eval_modes_out)); 1140 num_eval_modes_in[index] = 0; 1141 num_eval_modes_out[index] = 0; 1142 CeedCall(CeedRealloc(num_active_bases + 1, &eval_modes_in)); 1143 CeedCall(CeedRealloc(num_active_bases + 1, &eval_modes_out)); 1144 eval_modes_in[index] = NULL; 1145 eval_modes_out[index] = NULL; 1146 CeedCall(CeedRealloc(num_active_bases + 1, &eval_mode_offsets_in)); 1147 CeedCall(CeedRealloc(num_active_bases + 1, &eval_mode_offsets_out)); 1148 eval_mode_offsets_in[index] = NULL; 1149 eval_mode_offsets_out[index] = NULL; 1150 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->assembled_bases_in)); 1151 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->assembled_bases_out)); 1152 (*data)->assembled_bases_in[index] = NULL; 1153 (*data)->assembled_bases_out[index] = NULL; 1154 num_active_bases++; 1155 } 1156 switch (eval_mode) { 1157 case CEED_EVAL_NONE: 1158 case CEED_EVAL_INTERP: 1159 CeedCall(CeedRealloc(num_eval_modes_in[index] + 1, &eval_modes_in[index])); 1160 CeedCall(CeedRealloc(num_eval_modes_in[index] + 1, &eval_mode_offsets_in[index])); 1161 eval_modes_in[index][num_eval_modes_in[index]] = eval_mode; 1162 eval_mode_offsets_in[index][num_eval_modes_in[index]] = offset; 1163 offset += num_components; 1164 num_eval_modes_in[index] += 1; 1165 break; 1166 case CEED_EVAL_GRAD: 1167 CeedCall(CeedRealloc(num_eval_modes_in[index] + dim, &eval_modes_in[index])); 1168 CeedCall(CeedRealloc(num_eval_modes_in[index] + dim, &eval_mode_offsets_in[index])); 1169 for (CeedInt d = 0; d < dim; d++) { 1170 eval_modes_in[index][num_eval_modes_in[index] + d] = eval_mode; 1171 eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset; 1172 offset += num_components; 1173 } 1174 num_eval_modes_in[index] += dim; 1175 break; 1176 case CEED_EVAL_WEIGHT: 1177 case CEED_EVAL_DIV: 1178 case CEED_EVAL_CURL: 1179 break; // Caught by QF Assembly 1180 } 1181 } 1182 } 1183 (*data)->num_eval_modes_in = num_eval_modes_in; 1184 (*data)->eval_modes_in = eval_modes_in; 1185 (*data)->eval_mode_offsets_in = eval_mode_offsets_in; 1186 1187 // Determine active output basis 1188 CeedInt num_output_fields; 1189 1190 CeedCall(CeedQFunctionGetFields(qf, NULL, NULL, &num_output_fields, &qf_fields)); 1191 CeedCall(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 1192 offset = 0; 1193 for (CeedInt i = 0; i < num_output_fields; i++) { 1194 CeedVector vec; 1195 1196 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1197 if (vec == CEED_VECTOR_ACTIVE) { 1198 CeedInt index = -1, dim = 1, num_components; 1199 CeedBasis basis_out = NULL; 1200 CeedEvalMode eval_mode; 1201 1202 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 1203 CeedCall(CeedBasisGetDimension(basis_out, &dim)); 1204 CeedCall(CeedBasisGetNumComponents(basis_out, &num_components)); 1205 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1206 for (CeedInt i = 0; i < num_active_bases; i++) { 1207 if ((*data)->active_bases[i] == basis_out) index = i; 1208 } 1209 if (index == -1) { 1210 CeedElemRestriction elem_rstr_out; 1211 1212 index = num_active_bases; 1213 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->active_bases)); 1214 (*data)->active_bases[num_active_bases] = NULL; 1215 CeedCall(CeedBasisReferenceCopy(basis_out, &(*data)->active_bases[num_active_bases])); 1216 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->active_elem_rstrs)); 1217 (*data)->active_elem_rstrs[num_active_bases] = NULL; 1218 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_out)); 1219 CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_out, &(*data)->active_elem_rstrs[num_active_bases])); 1220 CeedCall(CeedRealloc(num_active_bases + 1, &num_eval_modes_in)); 1221 CeedCall(CeedRealloc(num_active_bases + 1, &num_eval_modes_out)); 1222 num_eval_modes_in[index] = 0; 1223 num_eval_modes_out[index] = 0; 1224 CeedCall(CeedRealloc(num_active_bases + 1, &eval_modes_in)); 1225 CeedCall(CeedRealloc(num_active_bases + 1, &eval_modes_out)); 1226 eval_modes_in[index] = NULL; 1227 eval_modes_out[index] = NULL; 1228 CeedCall(CeedRealloc(num_active_bases + 1, &eval_mode_offsets_in)); 1229 CeedCall(CeedRealloc(num_active_bases + 1, &eval_mode_offsets_out)); 1230 eval_mode_offsets_in[index] = NULL; 1231 eval_mode_offsets_out[index] = NULL; 1232 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->assembled_bases_in)); 1233 CeedCall(CeedRealloc(num_active_bases + 1, &(*data)->assembled_bases_out)); 1234 (*data)->assembled_bases_in[index] = NULL; 1235 (*data)->assembled_bases_out[index] = NULL; 1236 num_active_bases++; 1237 } 1238 switch (eval_mode) { 1239 case CEED_EVAL_NONE: 1240 case CEED_EVAL_INTERP: 1241 CeedCall(CeedRealloc(num_eval_modes_out[index] + 1, &eval_modes_out[index])); 1242 CeedCall(CeedRealloc(num_eval_modes_out[index] + 1, &eval_mode_offsets_out[index])); 1243 eval_modes_out[index][num_eval_modes_out[index]] = eval_mode; 1244 eval_mode_offsets_out[index][num_eval_modes_out[index]] = offset; 1245 offset += num_components; 1246 num_eval_modes_out[index] += 1; 1247 break; 1248 case CEED_EVAL_GRAD: 1249 CeedCall(CeedRealloc(num_eval_modes_out[index] + dim, &eval_modes_out[index])); 1250 CeedCall(CeedRealloc(num_eval_modes_out[index] + dim, &eval_mode_offsets_out[index])); 1251 for (CeedInt d = 0; d < dim; d++) { 1252 eval_modes_out[index][num_eval_modes_out[index] + d] = eval_mode; 1253 eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset; 1254 offset += num_components; 1255 } 1256 num_eval_modes_out[index] += dim; 1257 break; 1258 case CEED_EVAL_WEIGHT: 1259 case CEED_EVAL_DIV: 1260 case CEED_EVAL_CURL: 1261 break; // Caught by QF Assembly 1262 } 1263 } 1264 } 1265 (*data)->num_output_components = offset; 1266 (*data)->num_eval_modes_out = num_eval_modes_out; 1267 (*data)->eval_modes_out = eval_modes_out; 1268 (*data)->eval_mode_offsets_out = eval_mode_offsets_out; 1269 (*data)->num_active_bases = num_active_bases; 1270 1271 return CEED_ERROR_SUCCESS; 1272 } 1273 1274 /** 1275 @brief Get CeedOperator CeedEvalModes for assembly 1276 1277 @param[in] data CeedOperatorAssemblyData 1278 @param[out] num_active_bases Number of active bases 1279 @param[out] num_eval_mode_in Pointer to hold array of numbers of input CeedEvalModes, or NULL 1280 @param[out] eval_mode_in Pointer to hold arrays of input CeedEvalModes, or NULL 1281 @param[out] eval_mode_offsets_in Pointer to hold arrays of input offsets at each quadrature point 1282 @param[out] num_eval_mode_out Pointer to hold array of numbers of output CeedEvalModes, or NULL 1283 @param[out] eval_mode_out Pointer to hold arrays of output CeedEvalModes, or NULL 1284 @param[out] eval_mode_offsets_out Pointer to hold arrays of output offsets at each quadrature point 1285 @param[out] num_output_components Total number of output components at each quadrature point 1286 1287 @return An error code: 0 - success, otherwise - failure 1288 1289 @ref Backend 1290 **/ 1291 int CeedOperatorAssemblyDataGetEvalModes(CeedOperatorAssemblyData data, CeedInt *num_active_bases, CeedInt **num_eval_modes_in, 1292 const CeedEvalMode ***eval_modes_in, CeedSize ***eval_mode_offsets_in, CeedInt **num_eval_modes_out, 1293 const CeedEvalMode ***eval_modes_out, CeedSize ***eval_mode_offsets_out, CeedSize *num_output_components) { 1294 if (num_active_bases) *num_active_bases = data->num_active_bases; 1295 if (num_eval_modes_in) *num_eval_modes_in = data->num_eval_modes_in; 1296 if (eval_modes_in) *eval_modes_in = (const CeedEvalMode **)data->eval_modes_in; 1297 if (eval_mode_offsets_in) *eval_mode_offsets_in = data->eval_mode_offsets_in; 1298 if (num_eval_modes_out) *num_eval_modes_out = data->num_eval_modes_out; 1299 if (eval_modes_out) *eval_modes_out = (const CeedEvalMode **)data->eval_modes_out; 1300 if (eval_mode_offsets_out) *eval_mode_offsets_out = data->eval_mode_offsets_out; 1301 if (num_output_components) *num_output_components = data->num_output_components; 1302 1303 return CEED_ERROR_SUCCESS; 1304 } 1305 1306 /** 1307 @brief Get CeedOperator CeedBasis data for assembly 1308 1309 @param[in] data CeedOperatorAssemblyData 1310 @param[out] num_active_bases Number of active bases, or NULL 1311 @param[out] active_bases Pointer to hold active CeedBasis, or NULL 1312 @param[out] assembled_bases_in Pointer to hold assembled active input B, or NULL 1313 @param[out] assembled_bases_out Pointer to hold assembled active output B, or NULL 1314 1315 @return An error code: 0 - success, otherwise - failure 1316 1317 @ref Backend 1318 **/ 1319 int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num_active_bases, CeedBasis **active_bases, 1320 const CeedScalar ***assembled_bases_in, const CeedScalar ***assembled_bases_out) { 1321 // Assemble B_in, B_out if needed 1322 if (assembled_bases_in && !data->assembled_bases_in[0]) { 1323 CeedInt num_qpts; 1324 1325 CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases[0], &num_qpts)); 1326 for (CeedInt b = 0; b < data->num_active_bases; b++) { 1327 CeedInt elem_size; 1328 CeedScalar *B_in = NULL, *identity = NULL; 1329 const CeedScalar *interp_in, *grad_in; 1330 bool has_eval_none = false; 1331 1332 CeedCall(CeedBasisGetNumNodes(data->active_bases[b], &elem_size)); 1333 CeedCall(CeedCalloc(num_qpts * elem_size * data->num_eval_modes_in[b], &B_in)); 1334 1335 for (CeedInt i = 0; i < data->num_eval_modes_in[b]; i++) { 1336 has_eval_none = has_eval_none || (data->eval_modes_in[b][i] == CEED_EVAL_NONE); 1337 } 1338 if (has_eval_none) { 1339 CeedCall(CeedCalloc(num_qpts * elem_size, &identity)); 1340 for (CeedInt i = 0; i < (elem_size < num_qpts ? elem_size : num_qpts); i++) { 1341 identity[i * elem_size + i] = 1.0; 1342 } 1343 } 1344 CeedCall(CeedBasisGetInterp(data->active_bases[b], &interp_in)); 1345 CeedCall(CeedBasisGetGrad(data->active_bases[b], &grad_in)); 1346 1347 for (CeedInt q = 0; q < num_qpts; q++) { 1348 for (CeedInt n = 0; n < elem_size; n++) { 1349 CeedInt d_in = -1; 1350 for (CeedInt e_in = 0; e_in < data->num_eval_modes_in[b]; e_in++) { 1351 const CeedInt qq = data->num_eval_modes_in[b] * q; 1352 const CeedScalar *B = NULL; 1353 1354 if (data->eval_modes_in[b][e_in] == CEED_EVAL_GRAD) d_in++; 1355 CeedOperatorGetBasisPointer(data->eval_modes_in[b][e_in], identity, interp_in, &grad_in[d_in * num_qpts * elem_size], &B); 1356 B_in[(qq + e_in) * elem_size + n] = B[q * elem_size + n]; 1357 } 1358 } 1359 } 1360 if (identity) CeedCall(CeedFree(identity)); 1361 data->assembled_bases_in[b] = B_in; 1362 } 1363 } 1364 1365 if (assembled_bases_out && !data->assembled_bases_out[0]) { 1366 CeedInt num_qpts; 1367 1368 CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases[0], &num_qpts)); 1369 for (CeedInt b = 0; b < data->num_active_bases; b++) { 1370 CeedInt elem_size; 1371 const CeedScalar *interp_out, *grad_out; 1372 bool has_eval_none = false; 1373 CeedScalar *B_out = NULL, *identity = NULL; 1374 1375 CeedCall(CeedBasisGetNumNodes(data->active_bases[b], &elem_size)); 1376 CeedCall(CeedCalloc(num_qpts * elem_size * data->num_eval_modes_out[b], &B_out)); 1377 1378 for (CeedInt i = 0; i < data->num_eval_modes_out[b]; i++) { 1379 has_eval_none = has_eval_none || (data->eval_modes_out[b][i] == CEED_EVAL_NONE); 1380 } 1381 if (has_eval_none) { 1382 CeedCall(CeedCalloc(num_qpts * elem_size, &identity)); 1383 for (CeedInt i = 0; i < (elem_size < num_qpts ? elem_size : num_qpts); i++) { 1384 identity[i * elem_size + i] = 1.0; 1385 } 1386 } 1387 CeedCall(CeedBasisGetInterp(data->active_bases[b], &interp_out)); 1388 CeedCall(CeedBasisGetGrad(data->active_bases[b], &grad_out)); 1389 1390 for (CeedInt q = 0; q < num_qpts; q++) { 1391 for (CeedInt n = 0; n < elem_size; n++) { 1392 CeedInt d_out = -1; 1393 for (CeedInt e_out = 0; e_out < data->num_eval_modes_out[b]; e_out++) { 1394 const CeedInt qq = data->num_eval_modes_out[b] * q; 1395 const CeedScalar *B = NULL; 1396 1397 if (data->eval_modes_out[b][e_out] == CEED_EVAL_GRAD) d_out++; 1398 CeedOperatorGetBasisPointer(data->eval_modes_out[b][e_out], identity, interp_out, &grad_out[d_out * num_qpts * elem_size], &B); 1399 B_out[(qq + e_out) * elem_size + n] = B[q * elem_size + n]; 1400 } 1401 } 1402 } 1403 if (identity) CeedCall(CeedFree(identity)); 1404 data->assembled_bases_out[b] = B_out; 1405 } 1406 } 1407 1408 // Pass out assembled data 1409 if (active_bases) *active_bases = data->active_bases; 1410 if (assembled_bases_in) *assembled_bases_in = (const CeedScalar **)data->assembled_bases_in; 1411 if (assembled_bases_out) *assembled_bases_out = (const CeedScalar **)data->assembled_bases_out; 1412 1413 return CEED_ERROR_SUCCESS; 1414 } 1415 1416 /** 1417 @brief Get CeedOperator CeedBasis data for assembly 1418 1419 @param[in] data CeedOperatorAssemblyData 1420 @param[out] num_active_elem_rstrs Number of active element restrictions, or NULL 1421 @param[out] active_elem_rstrs Pointer to hold active CeedElemRestrictions, or NULL 1422 1423 @return An error code: 0 - success, otherwise - failure 1424 1425 @ref Backend 1426 **/ 1427 int CeedOperatorAssemblyDataGetElemRestrictions(CeedOperatorAssemblyData data, CeedInt *num_active_elem_rstrs, 1428 CeedElemRestriction **active_elem_rstrs) { 1429 if (num_active_elem_rstrs) *num_active_elem_rstrs = data->num_active_bases; 1430 if (active_elem_rstrs) *active_elem_rstrs = data->active_elem_rstrs; 1431 1432 return CEED_ERROR_SUCCESS; 1433 } 1434 1435 /** 1436 @brief Destroy CeedOperatorAssemblyData 1437 1438 @param[in,out] data CeedOperatorAssemblyData to destroy 1439 1440 @return An error code: 0 - success, otherwise - failure 1441 1442 @ref Backend 1443 **/ 1444 int CeedOperatorAssemblyDataDestroy(CeedOperatorAssemblyData *data) { 1445 if (!*data) return CEED_ERROR_SUCCESS; 1446 1447 CeedCall(CeedDestroy(&(*data)->ceed)); 1448 for (CeedInt b = 0; b < (*data)->num_active_bases; b++) { 1449 CeedCall(CeedBasisDestroy(&(*data)->active_bases[b])); 1450 CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs[b])); 1451 CeedCall(CeedFree(&(*data)->eval_modes_in[b])); 1452 CeedCall(CeedFree(&(*data)->eval_modes_out[b])); 1453 CeedCall(CeedFree(&(*data)->eval_mode_offsets_in[b])); 1454 CeedCall(CeedFree(&(*data)->eval_mode_offsets_out[b])); 1455 CeedCall(CeedFree(&(*data)->assembled_bases_in[b])); 1456 CeedCall(CeedFree(&(*data)->assembled_bases_out[b])); 1457 } 1458 CeedCall(CeedFree(&(*data)->active_bases)); 1459 CeedCall(CeedFree(&(*data)->active_elem_rstrs)); 1460 CeedCall(CeedFree(&(*data)->num_eval_modes_in)); 1461 CeedCall(CeedFree(&(*data)->num_eval_modes_out)); 1462 CeedCall(CeedFree(&(*data)->eval_modes_in)); 1463 CeedCall(CeedFree(&(*data)->eval_modes_out)); 1464 CeedCall(CeedFree(&(*data)->eval_mode_offsets_in)); 1465 CeedCall(CeedFree(&(*data)->eval_mode_offsets_out)); 1466 CeedCall(CeedFree(&(*data)->assembled_bases_in)); 1467 CeedCall(CeedFree(&(*data)->assembled_bases_out)); 1468 1469 CeedCall(CeedFree(data)); 1470 return CEED_ERROR_SUCCESS; 1471 } 1472 1473 /// @} 1474 1475 /// ---------------------------------------------------------------------------- 1476 /// CeedOperator Public API 1477 /// ---------------------------------------------------------------------------- 1478 /// @addtogroup CeedOperatorUser 1479 /// @{ 1480 1481 /** 1482 @brief Assemble a linear CeedQFunction associated with a CeedOperator 1483 1484 This returns a CeedVector containing a matrix at each quadrature point providing the action of the CeedQFunction associated with the CeedOperator. 1485 The vector 'assembled' is of shape [num_elements, num_input_fields, num_output_fields, num_quad_points] and contains column-major matrices 1486 representing the action of the CeedQFunction for a corresponding quadrature point on an element. Inputs and outputs are in the order provided by the 1487 user when adding CeedOperator fields. For example, a CeedQFunction with inputs 'u' and 'gradu' and outputs 'gradv' and 'v', provided in that order, 1488 would result in an assembled QFunction that consists of (1 + dim) x (dim + 1) matrices at each quadrature point acting on the input [u, du_0, du_1] 1489 and producing the output [dv_0, dv_1, v]. 1490 1491 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1492 1493 @param[in] op CeedOperator to assemble CeedQFunction 1494 @param[out] assembled CeedVector to store assembled CeedQFunction at quadrature points 1495 @param[out] rstr CeedElemRestriction for CeedVector containing assembled CeedQFunction 1496 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1497 1498 @return An error code: 0 - success, otherwise - failure 1499 1500 @ref User 1501 **/ 1502 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1503 CeedCall(CeedOperatorCheckReady(op)); 1504 1505 if (op->LinearAssembleQFunction) { 1506 // Backend version 1507 CeedCall(op->LinearAssembleQFunction(op, assembled, rstr, request)); 1508 } else { 1509 // Operator fallback 1510 CeedOperator op_fallback; 1511 1512 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1513 if (op_fallback) { 1514 CeedCall(CeedOperatorLinearAssembleQFunction(op_fallback, assembled, rstr, request)); 1515 } else { 1516 // LCOV_EXCL_START 1517 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunction"); 1518 // LCOV_EXCL_STOP 1519 } 1520 } 1521 return CEED_ERROR_SUCCESS; 1522 } 1523 1524 /** 1525 @brief Assemble CeedQFunction and store result internally. 1526 Return copied references of stored data to the caller. 1527 Caller is responsible for ownership and destruction of the copied references. 1528 See also @ref CeedOperatorLinearAssembleQFunction 1529 1530 @param[in] op CeedOperator to assemble CeedQFunction 1531 @param[out] assembled CeedVector to store assembled CeedQFunction at quadrature points 1532 @param[out] rstr CeedElemRestriction for CeedVector containing assembledCeedQFunction 1533 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1534 1535 @return An error code: 0 - success, otherwise - failure 1536 1537 @ref User 1538 **/ 1539 int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 1540 CeedCall(CeedOperatorCheckReady(op)); 1541 1542 if (op->LinearAssembleQFunctionUpdate) { 1543 // Backend version 1544 bool qf_assembled_is_setup; 1545 CeedVector assembled_vec = NULL; 1546 CeedElemRestriction assembled_rstr = NULL; 1547 1548 CeedCall(CeedQFunctionAssemblyDataIsSetup(op->qf_assembled, &qf_assembled_is_setup)); 1549 if (qf_assembled_is_setup) { 1550 bool update_needed; 1551 1552 CeedCall(CeedQFunctionAssemblyDataGetObjects(op->qf_assembled, &assembled_vec, &assembled_rstr)); 1553 CeedCall(CeedQFunctionAssemblyDataIsUpdateNeeded(op->qf_assembled, &update_needed)); 1554 if (update_needed) { 1555 CeedCall(op->LinearAssembleQFunctionUpdate(op, assembled_vec, assembled_rstr, request)); 1556 } 1557 } else { 1558 CeedCall(op->LinearAssembleQFunction(op, &assembled_vec, &assembled_rstr, request)); 1559 CeedCall(CeedQFunctionAssemblyDataSetObjects(op->qf_assembled, assembled_vec, assembled_rstr)); 1560 } 1561 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(op->qf_assembled, false)); 1562 1563 // Copy reference from internally held copy 1564 *assembled = NULL; 1565 *rstr = NULL; 1566 CeedCall(CeedVectorReferenceCopy(assembled_vec, assembled)); 1567 CeedCall(CeedVectorDestroy(&assembled_vec)); 1568 CeedCall(CeedElemRestrictionReferenceCopy(assembled_rstr, rstr)); 1569 CeedCall(CeedElemRestrictionDestroy(&assembled_rstr)); 1570 } else { 1571 // Operator fallback 1572 CeedOperator op_fallback; 1573 1574 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1575 if (op_fallback) { 1576 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op_fallback, assembled, rstr, request)); 1577 } else { 1578 // LCOV_EXCL_START 1579 return CeedError(op->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunctionUpdate"); 1580 // LCOV_EXCL_STOP 1581 } 1582 } 1583 1584 return CEED_ERROR_SUCCESS; 1585 } 1586 1587 /** 1588 @brief Assemble the diagonal of a square linear CeedOperator 1589 1590 This overwrites a CeedVector with the diagonal of a linear CeedOperator. 1591 1592 Note: Currently only non-composite CeedOperators with a single field and composite CeedOperators with single field sub-operators are supported. 1593 1594 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1595 1596 @param[in] op CeedOperator to assemble CeedQFunction 1597 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1598 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1599 1600 @return An error code: 0 - success, otherwise - failure 1601 1602 @ref User 1603 **/ 1604 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1605 bool is_composite; 1606 CeedCall(CeedOperatorCheckReady(op)); 1607 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1608 1609 CeedSize input_size = 0, output_size = 0; 1610 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1611 if (input_size != output_size) { 1612 // LCOV_EXCL_START 1613 return CeedError(op->ceed, CEED_ERROR_DIMENSION, "Operator must be square"); 1614 // LCOV_EXCL_STOP 1615 } 1616 1617 // Early exit for empty operator 1618 if (!is_composite) { 1619 CeedInt num_elem = 0; 1620 1621 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1622 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1623 } 1624 1625 if (op->LinearAssembleDiagonal) { 1626 // Backend version 1627 CeedCall(op->LinearAssembleDiagonal(op, assembled, request)); 1628 return CEED_ERROR_SUCCESS; 1629 } else if (op->LinearAssembleAddDiagonal) { 1630 // Backend version with zeroing first 1631 CeedCall(CeedVectorSetValue(assembled, 0.0)); 1632 CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); 1633 return CEED_ERROR_SUCCESS; 1634 } else { 1635 // Operator fallback 1636 CeedOperator op_fallback; 1637 1638 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1639 if (op_fallback) { 1640 CeedCall(CeedOperatorLinearAssembleDiagonal(op_fallback, assembled, request)); 1641 return CEED_ERROR_SUCCESS; 1642 } 1643 } 1644 // Default interface implementation 1645 CeedCall(CeedVectorSetValue(assembled, 0.0)); 1646 CeedCall(CeedOperatorLinearAssembleAddDiagonal(op, assembled, request)); 1647 1648 return CEED_ERROR_SUCCESS; 1649 } 1650 1651 /** 1652 @brief Assemble the diagonal of a square linear CeedOperator 1653 1654 This sums into a CeedVector the diagonal of a linear CeedOperator. 1655 1656 Note: Currently only non-composite CeedOperators with a single field and composite CeedOperators with single field sub-operators are supported. 1657 1658 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1659 1660 @param[in] op CeedOperator to assemble CeedQFunction 1661 @param[out] assembled CeedVector to store assembled CeedOperator diagonal 1662 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1663 1664 @return An error code: 0 - success, otherwise - failure 1665 1666 @ref User 1667 **/ 1668 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1669 bool is_composite; 1670 CeedCall(CeedOperatorCheckReady(op)); 1671 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1672 1673 CeedSize input_size = 0, output_size = 0; 1674 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1675 if (input_size != output_size) { 1676 // LCOV_EXCL_START 1677 return CeedError(op->ceed, CEED_ERROR_DIMENSION, "Operator must be square"); 1678 // LCOV_EXCL_STOP 1679 } 1680 1681 // Early exit for empty operator 1682 if (!is_composite) { 1683 CeedInt num_elem = 0; 1684 1685 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1686 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1687 } 1688 1689 if (op->LinearAssembleAddDiagonal) { 1690 // Backend version 1691 CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); 1692 return CEED_ERROR_SUCCESS; 1693 } else { 1694 // Operator fallback 1695 CeedOperator op_fallback; 1696 1697 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1698 if (op_fallback) { 1699 CeedCall(CeedOperatorLinearAssembleAddDiagonal(op_fallback, assembled, request)); 1700 return CEED_ERROR_SUCCESS; 1701 } 1702 } 1703 // Default interface implementation 1704 if (is_composite) { 1705 CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, false, assembled)); 1706 } else { 1707 CeedCall(CeedSingleOperatorAssembleAddDiagonal_Core(op, request, false, assembled)); 1708 } 1709 1710 return CEED_ERROR_SUCCESS; 1711 } 1712 1713 /** 1714 @brief Assemble the point block diagonal of a square linear CeedOperator 1715 1716 This overwrites a CeedVector with the point block diagonal of a linear CeedOperator. 1717 1718 Note: Currently only non-composite CeedOperators with a single field and composite CeedOperators with single field sub-operators are supported. 1719 1720 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1721 1722 @param[in] op CeedOperator to assemble CeedQFunction 1723 @param[out] assembled CeedVector to store assembled CeedOperator point block diagonal, provided in row-major form with an @a num_comp * @a num_comp 1724 block at each node. The dimensions of this vector are derived from the active vector for the CeedOperator. The array has shape [nodes, component out, 1725 component in]. 1726 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1727 1728 @return An error code: 0 - success, otherwise - failure 1729 1730 @ref User 1731 **/ 1732 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1733 bool is_composite; 1734 CeedCall(CeedOperatorCheckReady(op)); 1735 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1736 1737 CeedSize input_size = 0, output_size = 0; 1738 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1739 if (input_size != output_size) { 1740 // LCOV_EXCL_START 1741 return CeedError(op->ceed, CEED_ERROR_DIMENSION, "Operator must be square"); 1742 // LCOV_EXCL_STOP 1743 } 1744 1745 // Early exit for empty operator 1746 if (!is_composite) { 1747 CeedInt num_elem = 0; 1748 1749 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1750 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1751 } 1752 1753 if (op->LinearAssemblePointBlockDiagonal) { 1754 // Backend version 1755 CeedCall(op->LinearAssemblePointBlockDiagonal(op, assembled, request)); 1756 return CEED_ERROR_SUCCESS; 1757 } else if (op->LinearAssembleAddPointBlockDiagonal) { 1758 // Backend version with zeroing first 1759 CeedCall(CeedVectorSetValue(assembled, 0.0)); 1760 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); 1761 return CEED_ERROR_SUCCESS; 1762 } else { 1763 // Operator fallback 1764 CeedOperator op_fallback; 1765 1766 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1767 if (op_fallback) { 1768 CeedCall(CeedOperatorLinearAssemblePointBlockDiagonal(op_fallback, assembled, request)); 1769 return CEED_ERROR_SUCCESS; 1770 } 1771 } 1772 // Default interface implementation 1773 CeedCall(CeedVectorSetValue(assembled, 0.0)); 1774 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); 1775 1776 return CEED_ERROR_SUCCESS; 1777 } 1778 1779 /** 1780 @brief Assemble the point block diagonal of a square linear CeedOperator 1781 1782 This sums into a CeedVector with the point block diagonal of a linear CeedOperator. 1783 1784 Note: Currently only non-composite CeedOperators with a single field and composite CeedOperators with single field sub-operators are supported. 1785 1786 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1787 1788 @param[in] op CeedOperator to assemble CeedQFunction 1789 @param[out] assembled CeedVector to store assembled CeedOperator point block diagonal, provided in row-major form with an @a num_comp * @a num_comp 1790 block at each node. The dimensions of this vector are derived from the active vector for the CeedOperator. The array has shape [nodes, component out, 1791 component in]. 1792 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 1793 1794 @return An error code: 0 - success, otherwise - failure 1795 1796 @ref User 1797 **/ 1798 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 1799 bool is_composite; 1800 CeedCall(CeedOperatorCheckReady(op)); 1801 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1802 1803 CeedSize input_size = 0, output_size = 0; 1804 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 1805 if (input_size != output_size) { 1806 // LCOV_EXCL_START 1807 return CeedError(op->ceed, CEED_ERROR_DIMENSION, "Operator must be square"); 1808 // LCOV_EXCL_STOP 1809 } 1810 1811 // Early exit for empty operator 1812 if (!is_composite) { 1813 CeedInt num_elem = 0; 1814 1815 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1816 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1817 } 1818 1819 if (op->LinearAssembleAddPointBlockDiagonal) { 1820 // Backend version 1821 CeedCall(op->LinearAssembleAddPointBlockDiagonal(op, assembled, request)); 1822 return CEED_ERROR_SUCCESS; 1823 } else { 1824 // Operator fallback 1825 CeedOperator op_fallback; 1826 1827 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1828 if (op_fallback) { 1829 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op_fallback, assembled, request)); 1830 return CEED_ERROR_SUCCESS; 1831 } 1832 } 1833 // Default interface implementation 1834 if (is_composite) { 1835 CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, true, assembled)); 1836 } else { 1837 CeedCall(CeedSingleOperatorAssembleAddDiagonal_Core(op, request, true, assembled)); 1838 } 1839 1840 return CEED_ERROR_SUCCESS; 1841 } 1842 1843 /** 1844 @brief Fully assemble the nonzero pattern of a linear operator. 1845 1846 Expected to be used in conjunction with CeedOperatorLinearAssemble(). 1847 1848 The assembly routines use coordinate format, with num_entries tuples of the form (i, j, value) which indicate that value should be added to the 1849 matrix in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. This function returns the number of entries and their (i, j) 1850 locations, while CeedOperatorLinearAssemble() provides the values in the same ordering. 1851 1852 This will generally be slow unless your operator is low-order. 1853 1854 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1855 1856 @param[in] op CeedOperator to assemble 1857 @param[out] num_entries Number of entries in coordinate nonzero pattern 1858 @param[out] rows Row number for each entry 1859 @param[out] cols Column number for each entry 1860 1861 @ref User 1862 **/ 1863 int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { 1864 CeedInt num_suboperators, single_entries; 1865 CeedOperator *sub_operators; 1866 bool is_composite; 1867 CeedCall(CeedOperatorCheckReady(op)); 1868 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1869 1870 if (op->LinearAssembleSymbolic) { 1871 // Backend version 1872 CeedCall(op->LinearAssembleSymbolic(op, num_entries, rows, cols)); 1873 return CEED_ERROR_SUCCESS; 1874 } else { 1875 // Operator fallback 1876 CeedOperator op_fallback; 1877 1878 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1879 if (op_fallback) { 1880 CeedCall(CeedOperatorLinearAssembleSymbolic(op_fallback, num_entries, rows, cols)); 1881 return CEED_ERROR_SUCCESS; 1882 } 1883 } 1884 1885 // Default interface implementation 1886 1887 // count entries and allocate rows, cols arrays 1888 *num_entries = 0; 1889 if (is_composite) { 1890 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1891 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1892 for (CeedInt k = 0; k < num_suboperators; ++k) { 1893 CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries)); 1894 *num_entries += single_entries; 1895 } 1896 } else { 1897 CeedCall(CeedSingleOperatorAssemblyCountEntries(op, &single_entries)); 1898 *num_entries += single_entries; 1899 } 1900 CeedCall(CeedCalloc(*num_entries, rows)); 1901 CeedCall(CeedCalloc(*num_entries, cols)); 1902 1903 // assemble nonzero locations 1904 CeedInt offset = 0; 1905 if (is_composite) { 1906 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1907 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1908 for (CeedInt k = 0; k < num_suboperators; ++k) { 1909 CeedCall(CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows, *cols)); 1910 CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries)); 1911 offset += single_entries; 1912 } 1913 } else { 1914 CeedCall(CeedSingleOperatorAssembleSymbolic(op, offset, *rows, *cols)); 1915 } 1916 1917 return CEED_ERROR_SUCCESS; 1918 } 1919 1920 /** 1921 @brief Fully assemble the nonzero entries of a linear operator. 1922 1923 Expected to be used in conjunction with CeedOperatorLinearAssembleSymbolic(). 1924 1925 The assembly routines use coordinate format, with num_entries tuples of the form (i, j, value) which indicate that value should be added to the 1926 matrix in entry (i, j). Note that the (i, j) pairs are not unique and may repeat. This function returns the values of the nonzero entries to be added, 1927 their (i, j) locations are provided by CeedOperatorLinearAssembleSymbolic() 1928 1929 This will generally be slow unless your operator is low-order. 1930 1931 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1932 1933 @param[in] op CeedOperator to assemble 1934 @param[out] values Values to assemble into matrix 1935 1936 @ref User 1937 **/ 1938 int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { 1939 CeedInt num_suboperators, single_entries = 0; 1940 CeedOperator *sub_operators; 1941 bool is_composite; 1942 CeedCall(CeedOperatorCheckReady(op)); 1943 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1944 1945 // Early exit for empty operator 1946 if (!is_composite) { 1947 CeedInt num_elem = 0; 1948 1949 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 1950 if (num_elem == 0) return CEED_ERROR_SUCCESS; 1951 } 1952 1953 if (op->LinearAssemble) { 1954 // Backend version 1955 CeedCall(op->LinearAssemble(op, values)); 1956 return CEED_ERROR_SUCCESS; 1957 } else { 1958 // Operator fallback 1959 CeedOperator op_fallback; 1960 1961 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 1962 if (op_fallback) { 1963 CeedCall(CeedOperatorLinearAssemble(op_fallback, values)); 1964 return CEED_ERROR_SUCCESS; 1965 } 1966 } 1967 1968 // Default interface implementation 1969 CeedInt offset = 0; 1970 if (is_composite) { 1971 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 1972 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 1973 for (CeedInt k = 0; k < num_suboperators; k++) { 1974 CeedCall(CeedSingleOperatorAssemble(sub_operators[k], offset, values)); 1975 CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries)); 1976 offset += single_entries; 1977 } 1978 } else { 1979 CeedCall(CeedSingleOperatorAssemble(op, offset, values)); 1980 } 1981 1982 return CEED_ERROR_SUCCESS; 1983 } 1984 1985 /** 1986 @brief Get the multiplicity of nodes across suboperators in a composite CeedOperator 1987 1988 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 1989 1990 @param[in] op Composite CeedOperator 1991 @param[in] num_skip_indices Number of suboperators to skip 1992 @param[in] skip_indices Array of indices of suboperators to skip 1993 @param[out] mult Vector to store multiplicity (of size l_size) 1994 1995 @return An error code: 0 - success, otherwise - failure 1996 1997 @ref User 1998 **/ 1999 int CeedCompositeOperatorGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult) { 2000 CeedCall(CeedOperatorCheckReady(op)); 2001 2002 Ceed ceed; 2003 CeedInt num_suboperators; 2004 CeedSize l_vec_len; 2005 CeedScalar *mult_array; 2006 CeedVector ones_l_vec; 2007 CeedElemRestriction elem_rstr; 2008 CeedOperator *sub_operators; 2009 2010 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2011 2012 // Zero mult vector 2013 CeedCall(CeedVectorSetValue(mult, 0.0)); 2014 2015 // Get suboperators 2016 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 2017 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 2018 if (num_suboperators == 0) return CEED_ERROR_SUCCESS; 2019 2020 // Work vector 2021 CeedCall(CeedVectorGetLength(mult, &l_vec_len)); 2022 CeedCall(CeedVectorCreate(ceed, l_vec_len, &ones_l_vec)); 2023 CeedCall(CeedVectorSetValue(ones_l_vec, 1.0)); 2024 CeedCall(CeedVectorGetArray(mult, CEED_MEM_HOST, &mult_array)); 2025 2026 // Compute multiplicity across suboperators 2027 for (CeedInt i = 0; i < num_suboperators; i++) { 2028 const CeedScalar *sub_mult_array; 2029 CeedVector sub_mult_l_vec, ones_e_vec; 2030 2031 // -- Check for suboperator to skip 2032 for (CeedInt j = 0; j < num_skip_indices; j++) { 2033 if (skip_indices[j] == i) continue; 2034 } 2035 2036 // -- Sub operator multiplicity 2037 CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[i], &elem_rstr)); 2038 CeedCall(CeedElemRestrictionCreateVector(elem_rstr, &sub_mult_l_vec, &ones_e_vec)); 2039 CeedCall(CeedVectorSetValue(sub_mult_l_vec, 0.0)); 2040 CeedCall(CeedElemRestrictionApply(elem_rstr, CEED_NOTRANSPOSE, ones_l_vec, ones_e_vec, CEED_REQUEST_IMMEDIATE)); 2041 CeedCall(CeedElemRestrictionApply(elem_rstr, CEED_TRANSPOSE, ones_e_vec, sub_mult_l_vec, CEED_REQUEST_IMMEDIATE)); 2042 CeedCall(CeedVectorGetArrayRead(sub_mult_l_vec, CEED_MEM_HOST, &sub_mult_array)); 2043 // ---- Flag every node present in the current suboperator 2044 for (CeedInt j = 0; j < l_vec_len; j++) { 2045 if (sub_mult_array[j] > 0.0) mult_array[j] += 1.0; 2046 } 2047 CeedCall(CeedVectorRestoreArrayRead(sub_mult_l_vec, &sub_mult_array)); 2048 CeedCall(CeedVectorDestroy(&sub_mult_l_vec)); 2049 CeedCall(CeedVectorDestroy(&ones_e_vec)); 2050 } 2051 CeedCall(CeedVectorRestoreArray(mult, &mult_array)); 2052 CeedCall(CeedVectorDestroy(&ones_l_vec)); 2053 2054 return CEED_ERROR_SUCCESS; 2055 } 2056 2057 /** 2058 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator, creating the prolongation basis from the fine and coarse 2059 grid interpolation 2060 2061 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable. 2062 2063 @param[in] op_fine Fine grid operator 2064 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 2065 @param[in] rstr_coarse Coarse grid restriction 2066 @param[in] basis_coarse Coarse grid active vector basis 2067 @param[out] op_coarse Coarse grid operator 2068 @param[out] op_prolong Coarse to fine operator 2069 @param[out] op_restrict Fine to coarse operator 2070 2071 @return An error code: 0 - success, otherwise - failure 2072 2073 @ref User 2074 **/ 2075 int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2076 CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { 2077 CeedCall(CeedOperatorCheckReady(op_fine)); 2078 2079 // Build prolongation matrix 2080 CeedBasis basis_fine, basis_c_to_f; 2081 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2082 CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f)); 2083 2084 // Core code 2085 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict)); 2086 2087 return CEED_ERROR_SUCCESS; 2088 } 2089 2090 /** 2091 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator with a tensor basis for the active basis 2092 2093 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable. 2094 2095 @param[in] op_fine Fine grid operator 2096 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 2097 @param[in] rstr_coarse Coarse grid restriction 2098 @param[in] basis_coarse Coarse grid active vector basis 2099 @param[in] interp_c_to_f Matrix for coarse to fine interpolation 2100 @param[out] op_coarse Coarse grid operator 2101 @param[out] op_prolong Coarse to fine operator 2102 @param[out] op_restrict Fine to coarse operator 2103 2104 @return An error code: 0 - success, otherwise - failure 2105 2106 @ref User 2107 **/ 2108 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2109 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 2110 CeedOperator *op_restrict) { 2111 CeedCall(CeedOperatorCheckReady(op_fine)); 2112 Ceed ceed; 2113 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 2114 2115 // Check for compatible quadrature spaces 2116 CeedBasis basis_fine; 2117 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2118 CeedInt Q_f, Q_c; 2119 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); 2120 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); 2121 if (Q_f != Q_c) { 2122 // LCOV_EXCL_START 2123 return CeedError(ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces"); 2124 // LCOV_EXCL_STOP 2125 } 2126 2127 // Coarse to fine basis 2128 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c; 2129 CeedCall(CeedBasisGetDimension(basis_fine, &dim)); 2130 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); 2131 CeedCall(CeedBasisGetNumNodes1D(basis_fine, &P_1d_f)); 2132 CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); 2133 P_1d_c = dim == 1 ? num_nodes_c : dim == 2 ? sqrt(num_nodes_c) : cbrt(num_nodes_c); 2134 CeedScalar *q_ref, *q_weight, *grad; 2135 CeedCall(CeedCalloc(P_1d_f, &q_ref)); 2136 CeedCall(CeedCalloc(P_1d_f, &q_weight)); 2137 CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad)); 2138 CeedBasis basis_c_to_f; 2139 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_c, P_1d_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f)); 2140 CeedCall(CeedFree(&q_ref)); 2141 CeedCall(CeedFree(&q_weight)); 2142 CeedCall(CeedFree(&grad)); 2143 2144 // Core code 2145 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict)); 2146 return CEED_ERROR_SUCCESS; 2147 } 2148 2149 /** 2150 @brief Create a multigrid coarse operator and level transfer operators for a CeedOperator with a non-tensor basis for the active vector 2151 2152 Note: Calling this function asserts that setup is complete and sets all four CeedOperators as immutable. 2153 2154 @param[in] op_fine Fine grid operator 2155 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter 2156 @param[in] rstr_coarse Coarse grid restriction 2157 @param[in] basis_coarse Coarse grid active vector basis 2158 @param[in] interp_c_to_f Matrix for coarse to fine interpolation 2159 @param[out] op_coarse Coarse grid operator 2160 @param[out] op_prolong Coarse to fine operator 2161 @param[out] op_restrict Fine to coarse operator 2162 2163 @return An error code: 0 - success, otherwise - failure 2164 2165 @ref User 2166 **/ 2167 int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2168 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 2169 CeedOperator *op_restrict) { 2170 CeedCall(CeedOperatorCheckReady(op_fine)); 2171 Ceed ceed; 2172 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 2173 2174 // Check for compatible quadrature spaces 2175 CeedBasis basis_fine; 2176 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2177 CeedInt Q_f, Q_c; 2178 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); 2179 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); 2180 if (Q_f != Q_c) { 2181 // LCOV_EXCL_START 2182 return CeedError(ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces"); 2183 // LCOV_EXCL_STOP 2184 } 2185 2186 // Coarse to fine basis 2187 CeedElemTopology topo; 2188 CeedCall(CeedBasisGetTopology(basis_fine, &topo)); 2189 CeedInt dim, num_comp, num_nodes_c, num_nodes_f; 2190 CeedCall(CeedBasisGetDimension(basis_fine, &dim)); 2191 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); 2192 CeedCall(CeedBasisGetNumNodes(basis_fine, &num_nodes_f)); 2193 CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); 2194 CeedScalar *q_ref, *q_weight, *grad; 2195 CeedCall(CeedCalloc(num_nodes_f * dim, &q_ref)); 2196 CeedCall(CeedCalloc(num_nodes_f, &q_weight)); 2197 CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad)); 2198 CeedBasis basis_c_to_f; 2199 CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_c, num_nodes_f, interp_c_to_f, grad, q_ref, q_weight, &basis_c_to_f)); 2200 CeedCall(CeedFree(&q_ref)); 2201 CeedCall(CeedFree(&q_weight)); 2202 CeedCall(CeedFree(&grad)); 2203 2204 // Core code 2205 CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict)); 2206 return CEED_ERROR_SUCCESS; 2207 } 2208 2209 /** 2210 @brief Build a FDM based approximate inverse for each element for a CeedOperator 2211 2212 This returns a CeedOperator and CeedVector to apply a Fast Diagonalization Method based approximate inverse. 2213 This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, M = V^T V, K = V^T S V. 2214 The assembled QFunction is used to modify the eigenvalues from simultaneous diagonalization and obtain an approximate inverse of the form V^T 2215 S^hat V. The CeedOperator must be linear and non-composite. The associated CeedQFunction must therefore also be linear. 2216 2217 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 2218 2219 @param[in] op CeedOperator to create element inverses 2220 @param[out] fdm_inv CeedOperator to apply the action of a FDM based inverse for each element 2221 @param[in] request Address of CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2222 2223 @return An error code: 0 - success, otherwise - failure 2224 2225 @ref User 2226 **/ 2227 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, CeedRequest *request) { 2228 CeedCall(CeedOperatorCheckReady(op)); 2229 2230 if (op->CreateFDMElementInverse) { 2231 // Backend version 2232 CeedCall(op->CreateFDMElementInverse(op, fdm_inv, request)); 2233 return CEED_ERROR_SUCCESS; 2234 } else { 2235 // Operator fallback 2236 CeedOperator op_fallback; 2237 2238 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2239 if (op_fallback) { 2240 CeedCall(CeedOperatorCreateFDMElementInverse(op_fallback, fdm_inv, request)); 2241 return CEED_ERROR_SUCCESS; 2242 } 2243 } 2244 2245 // Default interface implementation 2246 Ceed ceed, ceed_parent; 2247 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2248 CeedCall(CeedGetOperatorFallbackParentCeed(ceed, &ceed_parent)); 2249 ceed_parent = ceed_parent ? ceed_parent : ceed; 2250 CeedQFunction qf; 2251 CeedCall(CeedOperatorGetQFunction(op, &qf)); 2252 2253 // Determine active input basis 2254 bool interp = false, grad = false; 2255 CeedBasis basis = NULL; 2256 CeedElemRestriction rstr = NULL; 2257 CeedOperatorField *op_fields; 2258 CeedQFunctionField *qf_fields; 2259 CeedInt num_input_fields; 2260 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_fields, NULL, NULL)); 2261 CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 2262 for (CeedInt i = 0; i < num_input_fields; i++) { 2263 CeedVector vec; 2264 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 2265 if (vec == CEED_VECTOR_ACTIVE) { 2266 CeedEvalMode eval_mode; 2267 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 2268 interp = interp || eval_mode == CEED_EVAL_INTERP; 2269 grad = grad || eval_mode == CEED_EVAL_GRAD; 2270 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 2271 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr)); 2272 } 2273 } 2274 if (!basis) { 2275 // LCOV_EXCL_START 2276 return CeedError(ceed, CEED_ERROR_BACKEND, "No active field set"); 2277 // LCOV_EXCL_STOP 2278 } 2279 CeedSize l_size = 1; 2280 CeedInt P_1d, Q_1d, elem_size, num_qpts, dim, num_comp = 1, num_elem = 1; 2281 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 2282 CeedCall(CeedBasisGetNumNodes(basis, &elem_size)); 2283 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 2284 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 2285 CeedCall(CeedBasisGetDimension(basis, &dim)); 2286 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 2287 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 2288 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 2289 2290 // Build and diagonalize 1D Mass and Laplacian 2291 bool tensor_basis; 2292 CeedCall(CeedBasisIsTensor(basis, &tensor_basis)); 2293 if (!tensor_basis) { 2294 // LCOV_EXCL_START 2295 return CeedError(ceed, CEED_ERROR_BACKEND, "FDMElementInverse only supported for tensor bases"); 2296 // LCOV_EXCL_STOP 2297 } 2298 CeedScalar *mass, *laplace, *x, *fdm_interp, *lambda; 2299 CeedCall(CeedCalloc(P_1d * P_1d, &mass)); 2300 CeedCall(CeedCalloc(P_1d * P_1d, &laplace)); 2301 CeedCall(CeedCalloc(P_1d * P_1d, &x)); 2302 CeedCall(CeedCalloc(P_1d * P_1d, &fdm_interp)); 2303 CeedCall(CeedCalloc(P_1d, &lambda)); 2304 // -- Build matrices 2305 const CeedScalar *interp_1d, *grad_1d, *q_weight_1d; 2306 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 2307 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 2308 CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d)); 2309 CeedCall(CeedBuildMassLaplace(interp_1d, grad_1d, q_weight_1d, P_1d, Q_1d, dim, mass, laplace)); 2310 2311 // -- Diagonalize 2312 CeedCall(CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d)); 2313 CeedCall(CeedFree(&mass)); 2314 CeedCall(CeedFree(&laplace)); 2315 for (CeedInt i = 0; i < P_1d; i++) { 2316 for (CeedInt j = 0; j < P_1d; j++) fdm_interp[i + j * P_1d] = x[j + i * P_1d]; 2317 } 2318 CeedCall(CeedFree(&x)); 2319 2320 // Assemble QFunction 2321 CeedVector assembled; 2322 CeedElemRestriction rstr_qf; 2323 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled, &rstr_qf, request)); 2324 CeedInt layout[3]; 2325 CeedCall(CeedElemRestrictionGetELayout(rstr_qf, &layout)); 2326 CeedCall(CeedElemRestrictionDestroy(&rstr_qf)); 2327 CeedScalar max_norm = 0; 2328 CeedCall(CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm)); 2329 2330 // Calculate element averages 2331 CeedInt num_modes = (interp ? 1 : 0) + (grad ? dim : 0); 2332 CeedScalar *elem_avg; 2333 const CeedScalar *assembled_array, *q_weight_array; 2334 CeedVector q_weight; 2335 CeedCall(CeedVectorCreate(ceed_parent, num_qpts, &q_weight)); 2336 CeedCall(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_weight)); 2337 CeedCall(CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array)); 2338 CeedCall(CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array)); 2339 CeedCall(CeedCalloc(num_elem, &elem_avg)); 2340 const CeedScalar qf_value_bound = max_norm * 100 * CEED_EPSILON; 2341 for (CeedInt e = 0; e < num_elem; e++) { 2342 CeedInt count = 0; 2343 for (CeedInt q = 0; q < num_qpts; q++) { 2344 for (CeedInt i = 0; i < num_comp * num_comp * num_modes * num_modes; i++) { 2345 if (fabs(assembled_array[q * layout[0] + i * layout[1] + e * layout[2]]) > qf_value_bound) { 2346 elem_avg[e] += assembled_array[q * layout[0] + i * layout[1] + e * layout[2]] / q_weight_array[q]; 2347 count++; 2348 } 2349 } 2350 } 2351 if (count) { 2352 elem_avg[e] /= count; 2353 } else { 2354 elem_avg[e] = 1.0; 2355 } 2356 } 2357 CeedCall(CeedVectorRestoreArrayRead(assembled, &assembled_array)); 2358 CeedCall(CeedVectorDestroy(&assembled)); 2359 CeedCall(CeedVectorRestoreArrayRead(q_weight, &q_weight_array)); 2360 CeedCall(CeedVectorDestroy(&q_weight)); 2361 2362 // Build FDM diagonal 2363 CeedVector q_data; 2364 CeedScalar *q_data_array, *fdm_diagonal; 2365 CeedCall(CeedCalloc(num_comp * elem_size, &fdm_diagonal)); 2366 const CeedScalar fdm_diagonal_bound = elem_size * CEED_EPSILON; 2367 for (CeedInt c = 0; c < num_comp; c++) { 2368 for (CeedInt n = 0; n < elem_size; n++) { 2369 if (interp) fdm_diagonal[c * elem_size + n] = 1.0; 2370 if (grad) { 2371 for (CeedInt d = 0; d < dim; d++) { 2372 CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d; 2373 fdm_diagonal[c * elem_size + n] += lambda[i]; 2374 } 2375 } 2376 if (fabs(fdm_diagonal[c * elem_size + n]) < fdm_diagonal_bound) fdm_diagonal[c * elem_size + n] = fdm_diagonal_bound; 2377 } 2378 } 2379 CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * elem_size, &q_data)); 2380 CeedCall(CeedVectorSetValue(q_data, 0.0)); 2381 CeedCall(CeedVectorGetArrayWrite(q_data, CEED_MEM_HOST, &q_data_array)); 2382 for (CeedInt e = 0; e < num_elem; e++) { 2383 for (CeedInt c = 0; c < num_comp; c++) { 2384 for (CeedInt n = 0; n < elem_size; n++) q_data_array[(e * num_comp + c) * elem_size + n] = 1. / (elem_avg[e] * fdm_diagonal[c * elem_size + n]); 2385 } 2386 } 2387 CeedCall(CeedFree(&elem_avg)); 2388 CeedCall(CeedFree(&fdm_diagonal)); 2389 CeedCall(CeedVectorRestoreArray(q_data, &q_data_array)); 2390 2391 // Setup FDM operator 2392 // -- Basis 2393 CeedBasis fdm_basis; 2394 CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy; 2395 CeedCall(CeedCalloc(P_1d * P_1d, &grad_dummy)); 2396 CeedCall(CeedCalloc(P_1d, &q_ref_dummy)); 2397 CeedCall(CeedCalloc(P_1d, &q_weight_dummy)); 2398 CeedCall(CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, fdm_interp, grad_dummy, q_ref_dummy, q_weight_dummy, &fdm_basis)); 2399 CeedCall(CeedFree(&fdm_interp)); 2400 CeedCall(CeedFree(&grad_dummy)); 2401 CeedCall(CeedFree(&q_ref_dummy)); 2402 CeedCall(CeedFree(&q_weight_dummy)); 2403 CeedCall(CeedFree(&lambda)); 2404 2405 // -- Restriction 2406 CeedElemRestriction rstr_qd_i; 2407 CeedInt strides[3] = {1, elem_size, elem_size * num_comp}; 2408 CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, elem_size, num_comp, num_elem * num_comp * elem_size, strides, &rstr_qd_i)); 2409 // -- QFunction 2410 CeedQFunction qf_fdm; 2411 CeedCall(CeedQFunctionCreateInteriorByName(ceed_parent, "Scale", &qf_fdm)); 2412 CeedCall(CeedQFunctionAddInput(qf_fdm, "input", num_comp, CEED_EVAL_INTERP)); 2413 CeedCall(CeedQFunctionAddInput(qf_fdm, "scale", num_comp, CEED_EVAL_NONE)); 2414 CeedCall(CeedQFunctionAddOutput(qf_fdm, "output", num_comp, CEED_EVAL_INTERP)); 2415 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_fdm, num_comp)); 2416 // -- QFunction context 2417 CeedInt *num_comp_data; 2418 CeedCall(CeedCalloc(1, &num_comp_data)); 2419 num_comp_data[0] = num_comp; 2420 CeedQFunctionContext ctx_fdm; 2421 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_fdm)); 2422 CeedCall(CeedQFunctionContextSetData(ctx_fdm, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_data), num_comp_data)); 2423 CeedCall(CeedQFunctionSetContext(qf_fdm, ctx_fdm)); 2424 CeedCall(CeedQFunctionContextDestroy(&ctx_fdm)); 2425 // -- Operator 2426 CeedCall(CeedOperatorCreate(ceed_parent, qf_fdm, NULL, NULL, fdm_inv)); 2427 CeedCall(CeedOperatorSetField(*fdm_inv, "input", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); 2428 CeedCall(CeedOperatorSetField(*fdm_inv, "scale", rstr_qd_i, CEED_BASIS_COLLOCATED, q_data)); 2429 CeedCall(CeedOperatorSetField(*fdm_inv, "output", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); 2430 2431 // Cleanup 2432 CeedCall(CeedVectorDestroy(&q_data)); 2433 CeedCall(CeedBasisDestroy(&fdm_basis)); 2434 CeedCall(CeedElemRestrictionDestroy(&rstr_qd_i)); 2435 CeedCall(CeedQFunctionDestroy(&qf_fdm)); 2436 2437 return CEED_ERROR_SUCCESS; 2438 } 2439 2440 /// @} 2441