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