1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed-impl.h> 9 #include <ceed.h> 10 #include <ceed/backend.h> 11 #include <assert.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 char *source_path_with_name = NULL; 39 CeedInt num_input_fields, num_output_fields; 40 CeedQFunctionField *input_fields, *output_fields; 41 42 // Check if NULL qf passed in 43 if (!qf) return CEED_ERROR_SUCCESS; 44 45 CeedDebug(CeedQFunctionReturnCeed(qf), "Creating fallback CeedQFunction\n"); 46 47 if (qf->source_path) { 48 size_t path_len = strlen(qf->source_path), name_len = strlen(qf->kernel_name); 49 50 CeedCall(CeedCalloc(path_len + name_len + 2, &source_path_with_name)); 51 memcpy(source_path_with_name, qf->source_path, path_len); 52 memcpy(&source_path_with_name[path_len], ":", 1); 53 memcpy(&source_path_with_name[path_len + 1], qf->kernel_name, name_len); 54 } else if (qf->user_source) { 55 CeedCall(CeedStringAllocCopy(qf->user_source, &source_path_with_name)); 56 } else { 57 CeedCall(CeedCalloc(1, &source_path_with_name)); 58 } 59 60 { 61 CeedInt vec_length; 62 CeedQFunctionUser f; 63 64 CeedCall(CeedQFunctionGetVectorLength(qf, &vec_length)); 65 CeedCall(CeedQFunctionGetUserFunction(qf, &f)); 66 CeedCall(CeedQFunctionCreateInterior(fallback_ceed, vec_length, f, source_path_with_name, qf_fallback)); 67 } 68 { 69 CeedQFunctionContext ctx; 70 71 CeedCall(CeedQFunctionGetContext(qf, &ctx)); 72 CeedCall(CeedQFunctionSetContext(*qf_fallback, ctx)); 73 CeedCall(CeedQFunctionContextDestroy(&ctx)); 74 } 75 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 76 for (CeedInt i = 0; i < num_input_fields; i++) { 77 const char *field_name; 78 CeedInt size; 79 CeedEvalMode eval_mode; 80 81 CeedCall(CeedQFunctionFieldGetData(input_fields[i], &field_name, &size, &eval_mode)); 82 CeedCall(CeedQFunctionAddInput(*qf_fallback, field_name, size, eval_mode)); 83 } 84 for (CeedInt i = 0; i < num_output_fields; i++) { 85 const char *field_name; 86 CeedInt size; 87 CeedEvalMode eval_mode; 88 89 CeedCall(CeedQFunctionFieldGetData(output_fields[i], &field_name, &size, &eval_mode)); 90 CeedCall(CeedQFunctionAddOutput(*qf_fallback, field_name, size, eval_mode)); 91 } 92 CeedCall(CeedFree(&source_path_with_name)); 93 return CEED_ERROR_SUCCESS; 94 } 95 96 /** 97 @brief Duplicate a `CeedOperator` with a reference `Ceed` to fallback for advanced `CeedOperator` functionality 98 99 @param[in,out] op `CeedOperator` to create fallback for 100 101 @return An error code: 0 - success, otherwise - failure 102 103 @ref Developer 104 **/ 105 static int CeedOperatorCreateFallback(CeedOperator op) { 106 bool is_composite; 107 Ceed ceed, ceed_fallback; 108 CeedOperator op_fallback; 109 110 // Check not already created 111 if (op->op_fallback) return CEED_ERROR_SUCCESS; 112 113 // Fallback Ceed 114 CeedCall(CeedOperatorGetCeed(op, &ceed)); 115 CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback)); 116 CeedCall(CeedDestroy(&ceed)); 117 if (!ceed_fallback) return CEED_ERROR_SUCCESS; 118 119 CeedDebug(CeedOperatorReturnCeed(op), "Creating fallback CeedOperator\n"); 120 121 // Clone Op 122 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 123 if (is_composite) { 124 CeedInt num_suboperators; 125 CeedOperator *sub_operators; 126 127 CeedCall(CeedOperatorCreateComposite(ceed_fallback, &op_fallback)); 128 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 129 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 130 for (CeedInt i = 0; i < num_suboperators; i++) { 131 CeedOperator op_sub_fallback; 132 133 CeedCall(CeedOperatorGetFallback(sub_operators[i], &op_sub_fallback)); 134 CeedCall(CeedOperatorCompositeAddSub(op_fallback, op_sub_fallback)); 135 } 136 } else { 137 bool is_at_points = false; 138 CeedInt num_input_fields, num_output_fields; 139 CeedQFunction qf_fallback = NULL, dqf_fallback = NULL, dqfT_fallback = NULL; 140 CeedOperatorField *input_fields, *output_fields; 141 142 CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->qf, &qf_fallback)); 143 CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqf, &dqf_fallback)); 144 CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqfT, &dqfT_fallback)); 145 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 146 if (is_at_points) { 147 CeedVector points; 148 CeedElemRestriction rstr_points; 149 150 CeedCall(CeedOperatorCreateAtPoints(ceed_fallback, qf_fallback, dqf_fallback, dqfT_fallback, &op_fallback)); 151 CeedCall(CeedOperatorAtPointsGetPoints(op, &rstr_points, &points)); 152 CeedCall(CeedOperatorAtPointsSetPoints(op_fallback, rstr_points, points)); 153 CeedCall(CeedVectorDestroy(&points)); 154 CeedCall(CeedElemRestrictionDestroy(&rstr_points)); 155 } else { 156 CeedCall(CeedOperatorCreate(ceed_fallback, qf_fallback, dqf_fallback, dqfT_fallback, &op_fallback)); 157 } 158 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 159 for (CeedInt i = 0; i < num_input_fields; i++) { 160 const char *field_name; 161 CeedVector vec; 162 CeedElemRestriction rstr; 163 CeedBasis basis; 164 165 CeedCall(CeedOperatorFieldGetData(input_fields[i], &field_name, &rstr, &basis, &vec)); 166 CeedCall(CeedOperatorSetField(op_fallback, field_name, rstr, basis, vec)); 167 CeedCall(CeedVectorDestroy(&vec)); 168 CeedCall(CeedElemRestrictionDestroy(&rstr)); 169 CeedCall(CeedBasisDestroy(&basis)); 170 } 171 for (CeedInt i = 0; i < num_output_fields; i++) { 172 const char *field_name; 173 CeedVector vec; 174 CeedElemRestriction rstr; 175 CeedBasis basis; 176 177 CeedCall(CeedOperatorFieldGetData(output_fields[i], &field_name, &rstr, &basis, &vec)); 178 CeedCall(CeedOperatorSetField(op_fallback, field_name, rstr, basis, vec)); 179 CeedCall(CeedVectorDestroy(&vec)); 180 CeedCall(CeedElemRestrictionDestroy(&rstr)); 181 CeedCall(CeedBasisDestroy(&basis)); 182 } 183 { 184 CeedQFunctionAssemblyData data; 185 186 CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); 187 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(data, &op_fallback->qf_assembled)); 188 } 189 // Cleanup 190 CeedCall(CeedQFunctionDestroy(&qf_fallback)); 191 CeedCall(CeedQFunctionDestroy(&dqf_fallback)); 192 CeedCall(CeedQFunctionDestroy(&dqfT_fallback)); 193 } 194 CeedCall(CeedOperatorSetName(op_fallback, op->name)); 195 CeedCall(CeedOperatorCheckReady(op_fallback)); 196 // Note: No ref-counting here so we don't get caught in a reference loop. 197 // The op holds the only reference to op_fallback and is responsible for deleting itself and op_fallback. 198 op->op_fallback = op_fallback; 199 op_fallback->op_fallback_parent = op; 200 CeedCall(CeedDestroy(&ceed_fallback)); 201 return CEED_ERROR_SUCCESS; 202 } 203 204 /** 205 @brief Core logic for assembling operator diagonal or point block diagonal 206 207 @param[in] op `CeedOperator` to assemble diagonal or point block diagonal 208 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 209 @param[in] is_point_block Boolean flag to assemble diagonal or point block diagonal 210 @param[out] assembled `CeedVector` to store assembled diagonal 211 212 @return An error code: 0 - success, otherwise - failure 213 214 @ref Developer 215 **/ 216 static inline int CeedOperatorLinearAssembleAddDiagonalSingle_Mesh(CeedOperator op, CeedRequest *request, const bool is_point_block, 217 CeedVector assembled) { 218 bool is_composite; 219 220 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 221 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 222 223 // Assemble QFunction 224 CeedInt layout_qf[3]; 225 const CeedScalar *assembled_qf_array; 226 CeedVector assembled_qf = NULL; 227 CeedElemRestriction assembled_elem_rstr = NULL; 228 229 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, request)); 230 CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, layout_qf)); 231 CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr)); 232 CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); 233 234 // Get assembly data 235 const CeedEvalMode **eval_modes_in, **eval_modes_out; 236 CeedInt num_active_bases_in, *num_eval_modes_in, num_active_bases_out, *num_eval_modes_out; 237 CeedSize **eval_mode_offsets_in, **eval_mode_offsets_out, num_output_components; 238 CeedBasis *active_bases_in, *active_bases_out; 239 CeedElemRestriction *active_elem_rstrs_in, *active_elem_rstrs_out; 240 CeedOperatorAssemblyData data; 241 242 CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); 243 CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases_in, &num_eval_modes_in, &eval_modes_in, &eval_mode_offsets_in, 244 &num_active_bases_out, &num_eval_modes_out, &eval_modes_out, &eval_mode_offsets_out, 245 &num_output_components)); 246 CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases_in, NULL, NULL, &active_bases_out, NULL)); 247 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, NULL, &active_elem_rstrs_in, NULL, &active_elem_rstrs_out)); 248 249 // Loop over all active bases (find matching input/output pairs) 250 for (CeedInt b = 0; b < CeedIntMin(num_active_bases_in, num_active_bases_out); b++) { 251 CeedInt b_in, b_out, num_elem, num_nodes, num_qpts, num_comp; 252 bool has_eval_none = false; 253 CeedScalar *elem_diag_array, *identity = NULL; 254 CeedVector elem_diag; 255 CeedElemRestriction diag_elem_rstr; 256 257 if (num_active_bases_in <= num_active_bases_out) { 258 b_in = b; 259 for (b_out = 0; b_out < num_active_bases_out; b_out++) { 260 if (active_bases_in[b_in] == active_bases_out[b_out]) { 261 break; 262 } 263 } 264 if (b_out == num_active_bases_out) { 265 continue; 266 } // No matching output basis found 267 } else { 268 b_out = b; 269 for (b_in = 0; b_in < num_active_bases_in; b_in++) { 270 if (active_bases_in[b_in] == active_bases_out[b_out]) { 271 break; 272 } 273 } 274 if (b_in == num_active_bases_in) { 275 continue; 276 } // No matching output basis found 277 } 278 CeedCheck(active_elem_rstrs_in[b_in] == active_elem_rstrs_out[b_out], CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 279 "Cannot assemble operator diagonal with different input and output active element restrictions"); 280 281 // Assemble point block diagonal restriction, if needed 282 if (is_point_block) { 283 CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstrs_in[b_in], &diag_elem_rstr)); 284 } else { 285 CeedCall(CeedElemRestrictionCreateUnsignedCopy(active_elem_rstrs_in[b_in], &diag_elem_rstr)); 286 } 287 288 // Create diagonal vector 289 CeedCall(CeedElemRestrictionCreateVector(diag_elem_rstr, NULL, &elem_diag)); 290 291 // Assemble element operator diagonals 292 CeedCall(CeedVectorSetValue(elem_diag, 0.0)); 293 CeedCall(CeedVectorGetArray(elem_diag, CEED_MEM_HOST, &elem_diag_array)); 294 CeedCall(CeedElemRestrictionGetNumElements(diag_elem_rstr, &num_elem)); 295 CeedCall(CeedBasisGetNumNodes(active_bases_in[b_in], &num_nodes)); 296 CeedCall(CeedBasisGetNumComponents(active_bases_in[b_in], &num_comp)); 297 if (active_bases_in[b_in] == CEED_BASIS_NONE) num_qpts = num_nodes; 298 else CeedCall(CeedBasisGetNumQuadraturePoints(active_bases_in[b_in], &num_qpts)); 299 300 // Construct identity matrix for basis if required 301 for (CeedInt i = 0; i < num_eval_modes_in[b_in]; i++) { 302 has_eval_none = has_eval_none || (eval_modes_in[b_in][i] == CEED_EVAL_NONE); 303 } 304 for (CeedInt i = 0; i < num_eval_modes_out[b_out]; i++) { 305 has_eval_none = has_eval_none || (eval_modes_out[b_out][i] == CEED_EVAL_NONE); 306 } 307 if (has_eval_none) { 308 CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); 309 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; 310 } 311 312 // Compute the diagonal of B^T D B 313 // Each element 314 for (CeedSize e = 0; e < num_elem; e++) { 315 // Each basis eval mode pair 316 CeedInt d_out = 0, q_comp_out; 317 CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; 318 319 for (CeedInt e_out = 0; e_out < num_eval_modes_out[b_out]; e_out++) { 320 CeedInt d_in = 0, q_comp_in; 321 const CeedScalar *B_t = NULL; 322 CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; 323 324 CeedCall(CeedOperatorGetBasisPointer(active_bases_out[b_out], eval_modes_out[b_out][e_out], identity, &B_t)); 325 CeedCall(CeedBasisGetNumQuadratureComponents(active_bases_out[b_out], eval_modes_out[b_out][e_out], &q_comp_out)); 326 if (q_comp_out > 1) { 327 if (e_out == 0 || eval_modes_out[b_out][e_out] != eval_mode_out_prev) d_out = 0; 328 else B_t = &B_t[(++d_out) * num_qpts * num_nodes]; 329 } 330 eval_mode_out_prev = eval_modes_out[b_out][e_out]; 331 332 for (CeedInt e_in = 0; e_in < num_eval_modes_in[b_in]; e_in++) { 333 const CeedScalar *B = NULL; 334 335 CeedCall(CeedOperatorGetBasisPointer(active_bases_in[b_in], eval_modes_in[b_in][e_in], identity, &B)); 336 CeedCall(CeedBasisGetNumQuadratureComponents(active_bases_in[b_in], eval_modes_in[b_in][e_in], &q_comp_in)); 337 if (q_comp_in > 1) { 338 if (e_in == 0 || eval_modes_in[b_in][e_in] != eval_mode_in_prev) d_in = 0; 339 else B = &B[(++d_in) * num_qpts * num_nodes]; 340 } 341 eval_mode_in_prev = eval_modes_in[b_in][e_in]; 342 343 // Each component 344 for (CeedInt c_out = 0; c_out < num_comp; c_out++) { 345 // Each qpt/node pair 346 for (CeedInt q = 0; q < num_qpts; q++) { 347 if (is_point_block) { 348 // Point Block Diagonal 349 for (CeedInt c_in = 0; c_in < num_comp; c_in++) { 350 const CeedSize c_offset = 351 (eval_mode_offsets_in[b_in][e_in] + c_in) * num_output_components + eval_mode_offsets_out[b_out][e_out] + c_out; 352 const CeedScalar qf_value = assembled_qf_array[q * layout_qf[0] + c_offset * layout_qf[1] + e * layout_qf[2]]; 353 354 for (CeedInt n = 0; n < num_nodes; n++) { 355 elem_diag_array[((e * num_comp + c_out) * num_comp + c_in) * num_nodes + n] += 356 B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n]; 357 } 358 } 359 } else { 360 // Diagonal Only 361 const CeedInt c_offset = 362 (eval_mode_offsets_in[b_in][e_in] + c_out) * num_output_components + eval_mode_offsets_out[b_out][e_out] + c_out; 363 const CeedScalar qf_value = assembled_qf_array[q * layout_qf[0] + c_offset * layout_qf[1] + e * layout_qf[2]]; 364 365 for (CeedInt n = 0; n < num_nodes; n++) { 366 elem_diag_array[(e * num_comp + c_out) * num_nodes + n] += B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n]; 367 } 368 } 369 } 370 } 371 } 372 } 373 } 374 CeedCall(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); 375 376 // Assemble local operator diagonal 377 CeedCall(CeedElemRestrictionApply(diag_elem_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); 378 379 // Cleanup 380 CeedCall(CeedElemRestrictionDestroy(&diag_elem_rstr)); 381 CeedCall(CeedVectorDestroy(&elem_diag)); 382 CeedCall(CeedFree(&identity)); 383 } 384 CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 385 CeedCall(CeedVectorDestroy(&assembled_qf)); 386 return CEED_ERROR_SUCCESS; 387 } 388 389 /** 390 @brief Core logic for assembling operator diagonal or point block diagonal 391 392 @param[in] op `CeedOperator` to assemble diagonal or point block diagonal 393 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 394 @param[in] is_point_block Boolean flag to assemble diagonal or point block diagonal 395 @param[out] assembled `CeedVector` to store assembled diagonal 396 397 @return An error code: 0 - success, otherwise - failure 398 399 @ref Developer 400 **/ 401 static inline int CeedOperatorLinearAssembleAddDiagonalSingle(CeedOperator op, CeedRequest *request, const bool is_point_block, 402 CeedVector assembled) { 403 bool is_at_points; 404 405 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 406 CeedCheck(!is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "AtPoints operator not supported"); 407 CeedCall(CeedOperatorLinearAssembleAddDiagonalSingle_Mesh(op, request, is_point_block, assembled)); 408 return CEED_ERROR_SUCCESS; 409 } 410 411 /** 412 @brief Core logic for assembling composite operator diagonal 413 414 @param[in] op `CeedOperator` to assemble point block diagonal 415 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 416 @param[in] is_point_block Boolean flag to assemble diagonal or point block diagonal 417 @param[out] assembled `CeedVector` to store assembled diagonal 418 419 @return An error code: 0 - success, otherwise - failure 420 421 @ref Developer 422 **/ 423 static inline int CeedOperatorLinearAssembleAddDiagonalComposite(CeedOperator op, CeedRequest *request, const bool is_point_block, 424 CeedVector assembled) { 425 CeedInt num_sub; 426 CeedOperator *suboperators; 427 428 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub)); 429 CeedCall(CeedOperatorCompositeGetSubList(op, &suboperators)); 430 for (CeedInt i = 0; i < num_sub; i++) { 431 if (is_point_block) { 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 CeedOperator`. 442 443 Users should generally use @ref 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 CeedOperatorAssembleSymbolicSingle(CeedOperator op, CeedInt offset, CeedInt *rows, CeedInt *cols) { 455 Ceed ceed; 456 bool is_composite; 457 CeedSize num_nodes_in, num_nodes_out, local_num_entries, count = 0; 458 CeedInt num_elem_in, elem_size_in, num_comp_in, layout_er_in[3]; 459 CeedInt num_elem_out, elem_size_out, num_comp_out, layout_er_out[3]; 460 CeedScalar *array; 461 const CeedScalar *elem_dof_a_in, *elem_dof_a_out; 462 CeedVector index_vec_in, index_vec_out, elem_dof_in, elem_dof_out; 463 CeedElemRestriction elem_rstr_in, elem_rstr_out, index_elem_rstr_in, index_elem_rstr_out; 464 465 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 466 CeedCall(CeedOperatorGetCeed(op, &ceed)); 467 CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 468 469 CeedCall(CeedOperatorGetActiveVectorLengths(op, &num_nodes_in, &num_nodes_out)); 470 CeedCall(CeedOperatorGetActiveElemRestrictions(op, &elem_rstr_in, &elem_rstr_out)); 471 CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_in, &num_elem_in)); 472 CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_in, &elem_size_in)); 473 CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_in, &num_comp_in)); 474 CeedCall(CeedElemRestrictionGetELayout(elem_rstr_in, layout_er_in)); 475 476 // Determine elem_dof relation for input 477 CeedCall(CeedVectorCreate(ceed, num_nodes_in, &index_vec_in)); 478 CeedCall(CeedVectorGetArrayWrite(index_vec_in, CEED_MEM_HOST, &array)); 479 for (CeedSize i = 0; i < num_nodes_in; i++) array[i] = i; 480 CeedCall(CeedVectorRestoreArray(index_vec_in, &array)); 481 CeedCall(CeedVectorCreate(ceed, num_elem_in * elem_size_in * num_comp_in, &elem_dof_in)); 482 CeedCall(CeedVectorSetValue(elem_dof_in, 0.0)); 483 CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr_in, &index_elem_rstr_in)); 484 CeedCall(CeedElemRestrictionApply(index_elem_rstr_in, CEED_NOTRANSPOSE, index_vec_in, elem_dof_in, CEED_REQUEST_IMMEDIATE)); 485 CeedCall(CeedVectorGetArrayRead(elem_dof_in, CEED_MEM_HOST, &elem_dof_a_in)); 486 CeedCall(CeedVectorDestroy(&index_vec_in)); 487 CeedCall(CeedElemRestrictionDestroy(&index_elem_rstr_in)); 488 489 if (elem_rstr_in != elem_rstr_out) { 490 CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_out, &num_elem_out)); 491 CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, 492 "Active input and output operator restrictions must have the same number of elements." 493 " Input has %" CeedInt_FMT " elements; output has %" CeedInt_FMT "elements.", 494 num_elem_in, num_elem_out); 495 CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_out, &elem_size_out)); 496 CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_out, &num_comp_out)); 497 CeedCall(CeedElemRestrictionGetELayout(elem_rstr_out, layout_er_out)); 498 499 // Determine elem_dof relation for output 500 CeedCall(CeedVectorCreate(ceed, num_nodes_out, &index_vec_out)); 501 CeedCall(CeedVectorGetArrayWrite(index_vec_out, CEED_MEM_HOST, &array)); 502 for (CeedSize i = 0; i < num_nodes_out; i++) array[i] = i; 503 CeedCall(CeedVectorRestoreArray(index_vec_out, &array)); 504 CeedCall(CeedVectorCreate(ceed, num_elem_out * elem_size_out * num_comp_out, &elem_dof_out)); 505 CeedCall(CeedVectorSetValue(elem_dof_out, 0.0)); 506 CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr_out, &index_elem_rstr_out)); 507 CeedCall(CeedElemRestrictionApply(index_elem_rstr_out, CEED_NOTRANSPOSE, index_vec_out, elem_dof_out, CEED_REQUEST_IMMEDIATE)); 508 CeedCall(CeedVectorGetArrayRead(elem_dof_out, CEED_MEM_HOST, &elem_dof_a_out)); 509 CeedCall(CeedVectorDestroy(&index_vec_out)); 510 CeedCall(CeedElemRestrictionDestroy(&index_elem_rstr_out)); 511 } else { 512 num_elem_out = num_elem_in; 513 elem_size_out = elem_size_in; 514 num_comp_out = num_comp_in; 515 layout_er_out[0] = layout_er_in[0]; 516 layout_er_out[1] = layout_er_in[1]; 517 layout_er_out[2] = layout_er_in[2]; 518 elem_dof_a_out = elem_dof_a_in; 519 } 520 local_num_entries = (CeedSize)elem_size_out * num_comp_out * elem_size_in * num_comp_in * num_elem_in; 521 522 // Determine i, j locations for element matrices 523 for (CeedInt e = 0; e < num_elem_in; e++) { 524 for (CeedInt comp_in = 0; comp_in < num_comp_in; comp_in++) { 525 for (CeedInt comp_out = 0; comp_out < num_comp_out; comp_out++) { 526 for (CeedInt i = 0; i < elem_size_out; i++) { 527 for (CeedInt j = 0; j < elem_size_in; j++) { 528 const CeedInt elem_dof_index_row = i * layout_er_out[0] + comp_out * layout_er_out[1] + e * layout_er_out[2]; 529 const CeedInt elem_dof_index_col = j * layout_er_in[0] + comp_in * layout_er_in[1] + e * layout_er_in[2]; 530 const CeedInt row = elem_dof_a_out[elem_dof_index_row]; 531 const CeedInt col = elem_dof_a_in[elem_dof_index_col]; 532 533 rows[offset + count] = row; 534 cols[offset + count] = col; 535 count++; 536 } 537 } 538 } 539 } 540 } 541 CeedCheck(count == local_num_entries, ceed, CEED_ERROR_MAJOR, "Error computing assembled entries"); 542 CeedCall(CeedVectorRestoreArrayRead(elem_dof_in, &elem_dof_a_in)); 543 CeedCall(CeedVectorDestroy(&elem_dof_in)); 544 if (elem_rstr_in != elem_rstr_out) { 545 CeedCall(CeedVectorRestoreArrayRead(elem_dof_out, &elem_dof_a_out)); 546 CeedCall(CeedVectorDestroy(&elem_dof_out)); 547 } 548 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in)); 549 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out)); 550 CeedCall(CeedDestroy(&ceed)); 551 return CEED_ERROR_SUCCESS; 552 } 553 554 /** 555 @brief Core logic to assemble `CeedQFunction` and store result internally. 556 557 Return copied references of stored data to the caller. 558 Caller is responsible for ownership and destruction of the copied references. 559 See also @ref CeedOperatorLinearAssembleQFunction(). 560 561 Note: If the value of `assembled` or `rstr` passed to this function are non-`NULL` , then it is assumed that they hold valid pointers. 562 These objects will be destroyed if `*assembled` or `*rstr` is the only reference to the object. 563 564 @param[in] op `CeedOperator` to assemble `CeedQFunction` 565 @param[in] use_parent Boolean flag to check for fallback parent implementation 566 @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points 567 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 568 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 569 570 @return An error code: 0 - success, otherwise - failure 571 572 @ref User 573 **/ 574 static int CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(CeedOperator op, bool use_parent, CeedVector *assembled, CeedElemRestriction *rstr, 575 CeedRequest *request) { 576 int (*LinearAssembleQFunctionUpdate)(CeedOperator, CeedVector, CeedElemRestriction, CeedRequest *) = NULL; 577 CeedOperator op_assemble = NULL; 578 CeedOperator op_fallback_parent = NULL; 579 580 CeedCall(CeedOperatorCheckReady(op)); 581 582 // Determine if fallback parent or operator has implementation 583 CeedCall(CeedOperatorGetFallbackParent(op, &op_fallback_parent)); 584 if (op_fallback_parent && use_parent && op_fallback_parent->LinearAssembleQFunctionUpdate) { 585 // -- Backend version for op fallback parent is faster, if it exists 586 CeedDebug(CeedOperatorReturnCeed(op), "Using fallback parent for CeedOperatorLinearAssembleQFunctionBuildOrUpdate\n"); 587 LinearAssembleQFunctionUpdate = op_fallback_parent->LinearAssembleQFunctionUpdate; 588 op_assemble = op_fallback_parent; 589 } else if (op->LinearAssembleQFunctionUpdate) { 590 // -- Backend version for op 591 LinearAssembleQFunctionUpdate = op->LinearAssembleQFunctionUpdate; 592 op_assemble = op; 593 } 594 595 // Assemble QFunction 596 if (LinearAssembleQFunctionUpdate) { 597 // Backend or fallback parent version 598 CeedQFunctionAssemblyData data; 599 bool data_is_setup; 600 CeedVector assembled_vec = NULL; 601 CeedElemRestriction assembled_rstr = NULL; 602 603 CeedCall(CeedOperatorGetQFunctionAssemblyData(op, &data)); 604 CeedCall(CeedQFunctionAssemblyDataIsSetup(data, &data_is_setup)); 605 if (data_is_setup) { 606 bool update_needed; 607 608 CeedCall(CeedQFunctionAssemblyDataGetObjects(data, &assembled_vec, &assembled_rstr)); 609 CeedCall(CeedQFunctionAssemblyDataIsUpdateNeeded(data, &update_needed)); 610 if (update_needed) CeedCall(LinearAssembleQFunctionUpdate(op_assemble, assembled_vec, assembled_rstr, request)); 611 } else { 612 CeedCall(CeedOperatorLinearAssembleQFunction(op_assemble, &assembled_vec, &assembled_rstr, request)); 613 CeedCall(CeedQFunctionAssemblyDataSetObjects(data, assembled_vec, assembled_rstr)); 614 } 615 CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(data, false)); 616 617 // Copy reference from internally held copy 618 CeedCall(CeedVectorReferenceCopy(assembled_vec, assembled)); 619 CeedCall(CeedElemRestrictionReferenceCopy(assembled_rstr, rstr)); 620 CeedCall(CeedVectorDestroy(&assembled_vec)); 621 CeedCall(CeedElemRestrictionDestroy(&assembled_rstr)); 622 } else { 623 // Operator fallback 624 CeedOperator op_fallback; 625 626 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleQFunctionBuildOrUpdate\n"); 627 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 628 if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op_fallback, assembled, rstr, request)); 629 else return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunctionUpdate"); 630 } 631 return CEED_ERROR_SUCCESS; 632 } 633 634 /** 635 @brief Assemble `CeedQFunction` and store result internally, but do not use fallback parent. 636 637 Return copied references of stored data to the caller. 638 Caller is responsible for ownership and destruction of the copied references. 639 See also @ref CeedOperatorLinearAssembleQFunction(). 640 641 Note: If the value of `assembled` or `rstr` passed to this function are non-`NULL` , then it is assumed that they hold valid pointers. 642 These objects will be destroyed if `*assembled` or `*rstr` is the only reference to the object. 643 644 @param[in] op `CeedOperator` to assemble `CeedQFunction` 645 @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points 646 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 647 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 648 649 @return An error code: 0 - success, otherwise - failure 650 651 @ref Developer 652 **/ 653 int CeedOperatorLinearAssembleQFunctionBuildOrUpdateFallback(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, 654 CeedRequest *request) { 655 return CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(op, false, assembled, rstr, request); 656 } 657 658 /** 659 @brief Assemble nonzero entries for non-composite `CeedOperator`. 660 661 Users should generally use @ref CeedOperatorLinearAssemble(). 662 663 @param[in] op `CeedOperator` to assemble 664 @param[in] offset Offset for number of entries 665 @param[out] values Values to assemble into matrix 666 667 @return An error code: 0 - success, otherwise - failure 668 669 @ref Developer 670 **/ 671 int CeedOperatorAssembleSingle(CeedOperator op, CeedInt offset, CeedVector values) { 672 bool is_composite, is_at_points; 673 674 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 675 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 676 677 // Early exit for empty operator 678 { 679 CeedInt num_elem = 0; 680 681 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 682 if (num_elem == 0) return CEED_ERROR_SUCCESS; 683 } 684 685 if (op->LinearAssembleSingle) { 686 // Backend version 687 CeedCall(op->LinearAssembleSingle(op, offset, values)); 688 return CEED_ERROR_SUCCESS; 689 } else { 690 // Operator fallback 691 CeedOperator op_fallback; 692 693 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorAssembleSingle\n"); 694 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 695 if (op_fallback) { 696 CeedCall(CeedOperatorAssembleSingle(op_fallback, offset, values)); 697 return CEED_ERROR_SUCCESS; 698 } 699 } 700 701 CeedCall(CeedOperatorIsAtPoints(op, &is_at_points)); 702 CeedCheck(!is_at_points, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 703 "Backend does not implement CeedOperatorLinearAssemble for AtPoints operator"); 704 705 // Assemble QFunction 706 CeedInt layout_qf[3]; 707 const CeedScalar *assembled_qf_array; 708 CeedVector assembled_qf = NULL; 709 CeedElemRestriction assembled_elem_rstr = NULL; 710 711 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, CEED_REQUEST_IMMEDIATE)); 712 CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, layout_qf)); 713 CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr)); 714 CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); 715 716 // Get assembly data 717 CeedInt num_elem_in, elem_size_in, num_comp_in, num_qpts_in; 718 CeedInt num_elem_out, elem_size_out, num_comp_out, num_qpts_out; 719 CeedSize local_num_entries, count = 0; 720 const CeedEvalMode **eval_modes_in, **eval_modes_out; 721 CeedInt num_active_bases_in, *num_eval_modes_in, num_active_bases_out, *num_eval_modes_out; 722 CeedBasis *active_bases_in, *active_bases_out, basis_in, basis_out; 723 const CeedScalar **B_mats_in, **B_mats_out, *B_mat_in, *B_mat_out; 724 CeedElemRestriction elem_rstr_in, elem_rstr_out; 725 CeedRestrictionType elem_rstr_type_in, elem_rstr_type_out; 726 const bool *elem_rstr_orients_in = NULL, *elem_rstr_orients_out = NULL; 727 const CeedInt8 *elem_rstr_curl_orients_in = NULL, *elem_rstr_curl_orients_out = NULL; 728 CeedOperatorAssemblyData data; 729 730 CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); 731 CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases_in, &num_eval_modes_in, &eval_modes_in, NULL, &num_active_bases_out, 732 &num_eval_modes_out, &eval_modes_out, NULL, NULL)); 733 734 CeedCheck(num_active_bases_in == 1 && num_active_bases_out == 1, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 735 "Cannot assemble operator with multiple active bases"); 736 CeedCheck(num_eval_modes_in[0] > 0 && num_eval_modes_out[0] > 0, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 737 "Cannot assemble operator without inputs/outputs"); 738 739 CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases_in, &B_mats_in, NULL, &active_bases_out, &B_mats_out)); 740 CeedCall(CeedOperatorGetActiveElemRestrictions(op, &elem_rstr_in, &elem_rstr_out)); 741 basis_in = active_bases_in[0]; 742 basis_out = active_bases_out[0]; 743 B_mat_in = B_mats_in[0]; 744 B_mat_out = B_mats_out[0]; 745 746 CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_in, &num_elem_in)); 747 CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_in, &elem_size_in)); 748 CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_in, &num_comp_in)); 749 if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in; 750 else CeedCall(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in)); 751 752 CeedCall(CeedElemRestrictionGetType(elem_rstr_in, &elem_rstr_type_in)); 753 if (elem_rstr_type_in == CEED_RESTRICTION_ORIENTED) { 754 CeedCall(CeedElemRestrictionGetOrientations(elem_rstr_in, CEED_MEM_HOST, &elem_rstr_orients_in)); 755 } else if (elem_rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 756 CeedCall(CeedElemRestrictionGetCurlOrientations(elem_rstr_in, CEED_MEM_HOST, &elem_rstr_curl_orients_in)); 757 } 758 759 if (elem_rstr_in != elem_rstr_out) { 760 CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_out, &num_elem_out)); 761 CeedCheck(num_elem_in == num_elem_out, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 762 "Active input and output operator restrictions must have the same number of elements." 763 " Input has %" CeedInt_FMT " elements; output has %" CeedInt_FMT "elements.", 764 num_elem_in, num_elem_out); 765 CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_out, &elem_size_out)); 766 CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_out, &num_comp_out)); 767 if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out; 768 else CeedCall(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out)); 769 CeedCheck(num_qpts_in == num_qpts_out, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 770 "Active input and output bases must have the same number of quadrature points." 771 " Input has %" CeedInt_FMT " points; output has %" CeedInt_FMT "points.", 772 num_qpts_in, num_qpts_out); 773 774 CeedCall(CeedElemRestrictionGetType(elem_rstr_out, &elem_rstr_type_out)); 775 if (elem_rstr_type_out == CEED_RESTRICTION_ORIENTED) { 776 CeedCall(CeedElemRestrictionGetOrientations(elem_rstr_out, CEED_MEM_HOST, &elem_rstr_orients_out)); 777 } else if (elem_rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 778 CeedCall(CeedElemRestrictionGetCurlOrientations(elem_rstr_out, CEED_MEM_HOST, &elem_rstr_curl_orients_out)); 779 } 780 } else { 781 num_elem_out = num_elem_in; 782 elem_size_out = elem_size_in; 783 num_comp_out = num_comp_in; 784 num_qpts_out = num_qpts_in; 785 786 elem_rstr_orients_out = elem_rstr_orients_in; 787 elem_rstr_curl_orients_out = elem_rstr_curl_orients_in; 788 } 789 local_num_entries = (CeedSize)elem_size_out * num_comp_out * elem_size_in * num_comp_in * num_elem_in; 790 791 // Loop over elements and put in data structure 792 // We store B_mat_in, B_mat_out, BTD, elem_mat in row-major order 793 CeedTensorContract contract; 794 CeedScalar *vals, *BTD_mat = NULL, *elem_mat = NULL, *elem_mat_b = NULL; 795 796 CeedCall(CeedBasisGetTensorContract(basis_in, &contract)); 797 CeedCall(CeedCalloc(elem_size_out * num_qpts_in * num_eval_modes_in[0], &BTD_mat)); 798 CeedCall(CeedCalloc(elem_size_out * elem_size_in, &elem_mat)); 799 if (elem_rstr_curl_orients_in || elem_rstr_curl_orients_out) CeedCall(CeedCalloc(elem_size_out * elem_size_in, &elem_mat_b)); 800 801 CeedCall(CeedVectorGetArray(values, CEED_MEM_HOST, &vals)); 802 for (CeedSize e = 0; e < num_elem_in; e++) { 803 for (CeedInt comp_in = 0; comp_in < num_comp_in; comp_in++) { 804 for (CeedInt comp_out = 0; comp_out < num_comp_out; comp_out++) { 805 // Compute B^T*D 806 for (CeedSize n = 0; n < elem_size_out; n++) { 807 for (CeedSize q = 0; q < num_qpts_in; q++) { 808 for (CeedInt e_in = 0; e_in < num_eval_modes_in[0]; e_in++) { 809 const CeedSize btd_index = n * (num_qpts_in * num_eval_modes_in[0]) + q * num_eval_modes_in[0] + e_in; 810 CeedScalar sum = 0.0; 811 812 for (CeedInt e_out = 0; e_out < num_eval_modes_out[0]; e_out++) { 813 const CeedSize b_out_index = (q * num_eval_modes_out[0] + e_out) * elem_size_out + n; 814 const CeedSize eval_mode_index = ((e_in * num_comp_in + comp_in) * num_eval_modes_out[0] + e_out) * num_comp_out + comp_out; 815 const CeedSize qf_index = q * layout_qf[0] + eval_mode_index * layout_qf[1] + e * layout_qf[2]; 816 817 sum += B_mat_out[b_out_index] * assembled_qf_array[qf_index]; 818 } 819 BTD_mat[btd_index] = sum; 820 } 821 } 822 } 823 824 // Form element matrix itself (for each block component) 825 if (contract) { 826 CeedCall(CeedTensorContractApply(contract, 1, num_qpts_in * num_eval_modes_in[0], elem_size_in, elem_size_out, BTD_mat, CEED_NOTRANSPOSE, 827 false, B_mat_in, elem_mat)); 828 } else { 829 Ceed ceed; 830 831 CeedCall(CeedOperatorGetCeed(op, &ceed)); 832 CeedCall(CeedMatrixMatrixMultiply(ceed, BTD_mat, B_mat_in, elem_mat, elem_size_out, elem_size_in, num_qpts_in * num_eval_modes_in[0])); 833 CeedCall(CeedDestroy(&ceed)); 834 } 835 836 // Transform the element matrix if required 837 if (elem_rstr_orients_out) { 838 const bool *elem_orients = &elem_rstr_orients_out[e * elem_size_out]; 839 840 for (CeedInt i = 0; i < elem_size_out; i++) { 841 const double orient = elem_orients[i] ? -1.0 : 1.0; 842 843 for (CeedInt j = 0; j < elem_size_in; j++) { 844 elem_mat[i * elem_size_in + j] *= orient; 845 } 846 } 847 } else if (elem_rstr_curl_orients_out) { 848 const CeedInt8 *elem_curl_orients = &elem_rstr_curl_orients_out[e * 3 * elem_size_out]; 849 850 // T^T*(B^T*D*B) 851 memcpy(elem_mat_b, elem_mat, elem_size_out * elem_size_in * sizeof(CeedScalar)); 852 for (CeedInt i = 0; i < elem_size_out; i++) { 853 for (CeedInt j = 0; j < elem_size_in; j++) { 854 elem_mat[i * elem_size_in + j] = elem_mat_b[i * elem_size_in + j] * elem_curl_orients[3 * i + 1] + 855 (i > 0 ? elem_mat_b[(i - 1) * elem_size_in + j] * elem_curl_orients[3 * i - 1] : 0.0) + 856 (i < elem_size_out - 1 ? elem_mat_b[(i + 1) * elem_size_in + j] * elem_curl_orients[3 * i + 3] : 0.0); 857 } 858 } 859 } 860 if (elem_rstr_orients_in) { 861 const bool *elem_orients = &elem_rstr_orients_in[e * elem_size_in]; 862 863 for (CeedInt i = 0; i < elem_size_out; i++) { 864 for (CeedInt j = 0; j < elem_size_in; j++) { 865 elem_mat[i * elem_size_in + j] *= elem_orients[j] ? -1.0 : 1.0; 866 } 867 } 868 } else if (elem_rstr_curl_orients_in) { 869 const CeedInt8 *elem_curl_orients = &elem_rstr_curl_orients_in[e * 3 * elem_size_in]; 870 871 // (B^T*D*B)*T 872 memcpy(elem_mat_b, elem_mat, elem_size_out * elem_size_in * sizeof(CeedScalar)); 873 for (CeedInt i = 0; i < elem_size_out; i++) { 874 for (CeedInt j = 0; j < elem_size_in; j++) { 875 elem_mat[i * elem_size_in + j] = elem_mat_b[i * elem_size_in + j] * elem_curl_orients[3 * j + 1] + 876 (j > 0 ? elem_mat_b[i * elem_size_in + j - 1] * elem_curl_orients[3 * j - 1] : 0.0) + 877 (j < elem_size_in - 1 ? elem_mat_b[i * elem_size_in + j + 1] * elem_curl_orients[3 * j + 3] : 0.0); 878 } 879 } 880 } 881 882 // Put element matrix in coordinate data structure 883 for (CeedInt i = 0; i < elem_size_out; i++) { 884 for (CeedInt j = 0; j < elem_size_in; j++) { 885 vals[offset + count] = elem_mat[i * elem_size_in + j]; 886 count++; 887 } 888 } 889 } 890 } 891 } 892 CeedCheck(count == local_num_entries, CeedOperatorReturnCeed(op), CEED_ERROR_MAJOR, "Error computing entries"); 893 CeedCall(CeedVectorRestoreArray(values, &vals)); 894 895 // Cleanup 896 CeedCall(CeedFree(&BTD_mat)); 897 CeedCall(CeedFree(&elem_mat)); 898 CeedCall(CeedFree(&elem_mat_b)); 899 if (elem_rstr_type_in == CEED_RESTRICTION_ORIENTED) { 900 CeedCall(CeedElemRestrictionRestoreOrientations(elem_rstr_in, &elem_rstr_orients_in)); 901 } else if (elem_rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { 902 CeedCall(CeedElemRestrictionRestoreCurlOrientations(elem_rstr_in, &elem_rstr_curl_orients_in)); 903 } 904 if (elem_rstr_in != elem_rstr_out) { 905 if (elem_rstr_type_out == CEED_RESTRICTION_ORIENTED) { 906 CeedCall(CeedElemRestrictionRestoreOrientations(elem_rstr_out, &elem_rstr_orients_out)); 907 } else if (elem_rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { 908 CeedCall(CeedElemRestrictionRestoreCurlOrientations(elem_rstr_out, &elem_rstr_curl_orients_out)); 909 } 910 } 911 CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); 912 CeedCall(CeedVectorDestroy(&assembled_qf)); 913 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in)); 914 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out)); 915 return CEED_ERROR_SUCCESS; 916 } 917 918 /** 919 @brief Count number of entries for assembled `CeedOperator` 920 921 @param[in] op `CeedOperator` to assemble 922 @param[out] num_entries Number of entries in assembled representation 923 924 @return An error code: 0 - success, otherwise - failure 925 926 @ref Utility 927 **/ 928 static int CeedOperatorAssemblyCountEntriesSingle(CeedOperator op, CeedSize *num_entries) { 929 bool is_composite; 930 CeedInt num_elem_in, elem_size_in, num_comp_in, num_elem_out, elem_size_out, num_comp_out; 931 CeedElemRestriction rstr_in, rstr_out; 932 933 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 934 CeedCheck(!is_composite, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); 935 936 CeedCall(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); 937 CeedCall(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in)); 938 CeedCall(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); 939 CeedCall(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in)); 940 if (rstr_in != rstr_out) { 941 CeedCall(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out)); 942 CeedCheck(num_elem_in == num_elem_out, CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, 943 "Active input and output operator restrictions must have the same number of elements." 944 " Input has %" CeedInt_FMT " elements; output has %" CeedInt_FMT "elements.", 945 num_elem_in, num_elem_out); 946 CeedCall(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); 947 CeedCall(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out)); 948 } else { 949 num_elem_out = num_elem_in; 950 elem_size_out = elem_size_in; 951 num_comp_out = num_comp_in; 952 } 953 CeedCall(CeedElemRestrictionDestroy(&rstr_in)); 954 CeedCall(CeedElemRestrictionDestroy(&rstr_out)); 955 *num_entries = (CeedSize)elem_size_in * num_comp_in * elem_size_out * num_comp_out * num_elem_in; 956 return CEED_ERROR_SUCCESS; 957 } 958 959 /** 960 @brief Count number of entries for assembled `CeedOperator` 961 962 @param[in] op `CeedOperator` to assemble 963 @param[out] num_entries Number of entries in assembled representation 964 965 @return An error code: 0 - success, otherwise - failure 966 967 @ref Utility 968 **/ 969 int CeedOperatorLinearAssembleGetNumEntries(CeedOperator op, CeedSize *num_entries) { 970 bool is_composite; 971 972 CeedCall(CeedOperatorCheckReady(op)); 973 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 974 975 if (is_composite) { 976 CeedInt num_suboperators; 977 CeedOperator *sub_operators; 978 979 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 980 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 981 982 *num_entries = 0; 983 for (CeedInt k = 0; k < num_suboperators; ++k) { 984 CeedSize single_entries; 985 986 CeedCall(CeedOperatorAssemblyCountEntriesSingle(sub_operators[k], &single_entries)); 987 *num_entries += single_entries; 988 } 989 } else { 990 CeedCall(CeedOperatorAssemblyCountEntriesSingle(op, num_entries)); 991 } 992 return CEED_ERROR_SUCCESS; 993 } 994 995 /** 996 @brief Common code for creating a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` 997 998 @param[in] op_fine Fine grid `CeedOperator` 999 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 1000 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 1001 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 1002 @param[in] basis_c_to_f `CeedBasis` for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction operators 1003 @param[out] op_coarse Coarse grid `CeedOperator` 1004 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 1005 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 1006 1007 @return An error code: 0 - success, otherwise - failure 1008 1009 @ref Developer 1010 **/ 1011 static int CeedOperatorMultigridLevelCreateSingle_Core(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, 1012 CeedBasis basis_coarse, CeedBasis basis_c_to_f, CeedOperator *op_coarse, 1013 CeedOperator *op_prolong, CeedOperator *op_restrict) { 1014 bool is_composite; 1015 Ceed ceed; 1016 CeedInt dim = 0, num_comp, num_input_fields, num_output_fields; 1017 CeedVector mult_vec = NULL; 1018 CeedElemRestriction rstr_p_mult_fine = NULL, rstr_fine = NULL; 1019 CeedOperatorField *input_fields, *output_fields; 1020 1021 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 1022 1023 // Check for composite operator 1024 CeedCall(CeedOperatorIsComposite(op_fine, &is_composite)); 1025 CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Automatic multigrid setup for composite operators not supported"); 1026 1027 // Coarse Grid 1028 { 1029 bool is_at_points; 1030 1031 CeedCall(CeedOperatorIsAtPoints(op_fine, &is_at_points)); 1032 if (is_at_points) { 1033 CeedVector point_coords; 1034 CeedElemRestriction rstr_points; 1035 1036 CeedCall(CeedOperatorCreateAtPoints(ceed, op_fine->qf, op_fine->dqf, op_fine->dqfT, op_coarse)); 1037 CeedCall(CeedOperatorAtPointsGetPoints(op_fine, &rstr_points, &point_coords)); 1038 CeedCall(CeedOperatorAtPointsSetPoints(*op_coarse, rstr_points, point_coords)); 1039 CeedCall(CeedVectorDestroy(&point_coords)); 1040 CeedCall(CeedElemRestrictionDestroy(&rstr_points)); 1041 } else { 1042 CeedCall(CeedOperatorCreate(ceed, op_fine->qf, op_fine->dqf, op_fine->dqfT, op_coarse)); 1043 } 1044 } 1045 CeedCall(CeedOperatorGetFields(op_fine, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); 1046 // -- Clone input fields 1047 for (CeedInt i = 0; i < num_input_fields; i++) { 1048 const char *field_name; 1049 CeedVector vec; 1050 CeedElemRestriction rstr = NULL; 1051 CeedBasis basis = NULL; 1052 1053 CeedCall(CeedOperatorFieldGetName(input_fields[i], &field_name)); 1054 CeedCall(CeedOperatorFieldGetVector(input_fields[i], &vec)); 1055 if (vec == CEED_VECTOR_ACTIVE) { 1056 CeedCall(CeedElemRestrictionReferenceCopy(rstr_coarse, &rstr)); 1057 CeedCall(CeedBasisReferenceCopy(basis_coarse, &basis)); 1058 if (!rstr_fine) CeedCall(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_fine)); 1059 } else { 1060 CeedCall(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr)); 1061 CeedCall(CeedOperatorFieldGetBasis(input_fields[i], &basis)); 1062 } 1063 if (dim == 0) CeedCall(CeedBasisGetDimension(basis, &dim)); 1064 CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec)); 1065 CeedCall(CeedVectorDestroy(&vec)); 1066 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1067 CeedCall(CeedBasisDestroy(&basis)); 1068 } 1069 // -- Clone output fields 1070 for (CeedInt i = 0; i < num_output_fields; i++) { 1071 const char *field_name; 1072 CeedVector vec; 1073 CeedElemRestriction rstr = NULL; 1074 CeedBasis basis = NULL; 1075 1076 CeedCall(CeedOperatorFieldGetName(output_fields[i], &field_name)); 1077 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 1078 if (vec == CEED_VECTOR_ACTIVE) { 1079 CeedCall(CeedElemRestrictionReferenceCopy(rstr_coarse, &rstr)); 1080 CeedCall(CeedBasisReferenceCopy(basis_coarse, &basis)); 1081 if (!rstr_fine) CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_fine)); 1082 } else { 1083 CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr)); 1084 CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 1085 } 1086 if (dim == 0) CeedCall(CeedBasisGetDimension(basis, &dim)); 1087 CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec)); 1088 CeedCall(CeedVectorDestroy(&vec)); 1089 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1090 CeedCall(CeedBasisDestroy(&basis)); 1091 } 1092 dim = dim ? dim : 1; 1093 // -- Clone QFunctionAssemblyData 1094 { 1095 CeedQFunctionAssemblyData fine_data; 1096 1097 CeedCall(CeedOperatorGetQFunctionAssemblyData(op_fine, &fine_data)); 1098 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(fine_data, &(*op_coarse)->qf_assembled)); 1099 } 1100 1101 // Multiplicity vector 1102 bool use_scalar_mult = true; 1103 1104 if (op_restrict || op_prolong) { 1105 CeedInt num_elem, num_comp, elem_size; 1106 CeedVector mult_l_vec, mult_e_vec; 1107 CeedRestrictionType rstr_type; 1108 CeedElemRestriction rstr_p_mult_full; 1109 1110 CeedCall(CeedElemRestrictionGetType(rstr_fine, &rstr_type)); 1111 CeedCheck(rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_UNSUPPORTED, 1112 "Element restrictions created with CeedElemRestrictionCreateCurlOriented are not supported"); 1113 CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector"); 1114 1115 // Create multiplicity multi-component l-vector 1116 CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_full)); 1117 CeedCall(CeedElemRestrictionCreateVector(rstr_fine, &mult_l_vec, &mult_e_vec)); 1118 CeedCall(CeedVectorSetValue(mult_e_vec, 0.0)); 1119 CeedCall(CeedElemRestrictionApply(rstr_p_mult_full, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE)); 1120 CeedCall(CeedVectorSetValue(mult_l_vec, 0.0)); 1121 CeedCall(CeedElemRestrictionApply(rstr_p_mult_full, CEED_TRANSPOSE, mult_e_vec, mult_l_vec, CEED_REQUEST_IMMEDIATE)); 1122 CeedCall(CeedVectorReciprocal(mult_l_vec)); 1123 1124 // Determine to use scalar multiplicity or not 1125 { 1126 const CeedInt p = pow(elem_size, 1.0 / dim); 1127 1128 use_scalar_mult = num_comp > 1 && (dim < 3 || num_comp - 1 > (3 * (pow(p, dim - 1) - pow(p, dim - 2)) + 1) / pow(p - 1, dim)); 1129 } 1130 1131 if (use_scalar_mult) { 1132 // Create multiplicity single component e-vector 1133 CeedCall(CeedElemRestrictionGetNumElements(rstr_p_mult_full, &num_elem)); 1134 CeedCall(CeedElemRestrictionGetNumComponents(rstr_p_mult_full, &num_comp)); 1135 CeedCall(CeedElemRestrictionGetElementSize(rstr_p_mult_full, &elem_size)); 1136 CeedCall(CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, 1, num_elem * elem_size, CEED_STRIDES_BACKEND, &rstr_p_mult_fine)); 1137 CeedCall(CeedElemRestrictionCreateVector(rstr_p_mult_fine, &mult_vec, NULL)); 1138 { 1139 CeedQFunction qf_to_scalar; 1140 CeedOperator op_to_scalar; 1141 1142 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Identity to scalar", &qf_to_scalar)); 1143 CeedCall(CeedQFunctionAddInput(qf_to_scalar, "input", num_comp, CEED_EVAL_NONE)); 1144 CeedCall(CeedQFunctionAddOutput(qf_to_scalar, "output", 1, CEED_EVAL_NONE)); 1145 1146 CeedCall(CeedOperatorCreate(ceed, qf_to_scalar, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_to_scalar)); 1147 CeedCall(CeedOperatorSetField(op_to_scalar, "input", rstr_p_mult_full, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1148 CeedCall(CeedOperatorSetField(op_to_scalar, "output", rstr_p_mult_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1149 1150 CeedCall(CeedOperatorApply(op_to_scalar, mult_l_vec, mult_vec, CEED_REQUEST_IMMEDIATE)); 1151 1152 // Clean-up 1153 CeedCall(CeedQFunctionDestroy(&qf_to_scalar)); 1154 CeedCall(CeedOperatorDestroy(&op_to_scalar)); 1155 } 1156 } else { 1157 mult_vec = NULL; 1158 CeedCall(CeedVectorReferenceCopy(mult_l_vec, &mult_vec)); 1159 rstr_p_mult_fine = NULL; 1160 CeedCall(CeedElemRestrictionReferenceCopy(rstr_p_mult_full, &rstr_p_mult_fine)); 1161 } 1162 // Clean-up 1163 CeedCall(CeedVectorDestroy(&mult_e_vec)); 1164 CeedCall(CeedVectorDestroy(&mult_l_vec)); 1165 CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_full)); 1166 } 1167 1168 // Clone name 1169 bool has_name = op_fine->name; 1170 size_t name_len = op_fine->name ? strlen(op_fine->name) : 0; 1171 CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name)); 1172 1173 // Check that coarse to fine basis is provided if prolong/restrict operators are requested 1174 CeedCheck(basis_c_to_f || (!op_restrict && !op_prolong), ceed, CEED_ERROR_INCOMPATIBLE, 1175 "Prolongation or restriction operator creation requires coarse-to-fine basis"); 1176 1177 // Restriction/Prolongation Operators 1178 CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp)); 1179 1180 // Restriction 1181 if (op_restrict) { 1182 CeedInt *num_comp_r_data; 1183 CeedQFunctionContext ctx_r; 1184 CeedQFunction qf_restrict; 1185 1186 CeedCall(CeedQFunctionCreateInteriorByName(ceed, use_scalar_mult ? "Scale (scalar)" : "Scale", &qf_restrict)); 1187 CeedCall(CeedCalloc(1, &num_comp_r_data)); 1188 num_comp_r_data[0] = num_comp; 1189 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r)); 1190 CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data)); 1191 CeedCall(CeedQFunctionSetContext(qf_restrict, ctx_r)); 1192 CeedCall(CeedQFunctionContextDestroy(&ctx_r)); 1193 CeedCall(CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE)); 1194 CeedCall(CeedQFunctionAddInput(qf_restrict, "scale", use_scalar_mult ? 1 : num_comp, CEED_EVAL_NONE)); 1195 CeedCall(CeedQFunctionAddOutput(qf_restrict, "output", num_comp, CEED_EVAL_INTERP)); 1196 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_restrict, num_comp)); 1197 1198 CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict)); 1199 CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1200 CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec)); 1201 CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); 1202 1203 // Set name 1204 char *restriction_name; 1205 1206 CeedCall(CeedCalloc(17 + name_len, &restriction_name)); 1207 sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); 1208 CeedCall(CeedOperatorSetName(*op_restrict, restriction_name)); 1209 CeedCall(CeedFree(&restriction_name)); 1210 1211 // Check 1212 CeedCall(CeedOperatorCheckReady(*op_restrict)); 1213 1214 // Cleanup 1215 CeedCall(CeedQFunctionDestroy(&qf_restrict)); 1216 } 1217 1218 // Prolongation 1219 if (op_prolong) { 1220 CeedInt *num_comp_p_data; 1221 CeedQFunctionContext ctx_p; 1222 CeedQFunction qf_prolong; 1223 1224 CeedCall(CeedQFunctionCreateInteriorByName(ceed, use_scalar_mult ? "Scale (scalar)" : "Scale", &qf_prolong)); 1225 CeedCall(CeedCalloc(1, &num_comp_p_data)); 1226 num_comp_p_data[0] = num_comp; 1227 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_p)); 1228 CeedCall(CeedQFunctionContextSetData(ctx_p, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_p_data), num_comp_p_data)); 1229 CeedCall(CeedQFunctionSetContext(qf_prolong, ctx_p)); 1230 CeedCall(CeedQFunctionContextDestroy(&ctx_p)); 1231 CeedCall(CeedQFunctionAddInput(qf_prolong, "input", num_comp, CEED_EVAL_INTERP)); 1232 CeedCall(CeedQFunctionAddInput(qf_prolong, "scale", use_scalar_mult ? 1 : num_comp, CEED_EVAL_NONE)); 1233 CeedCall(CeedQFunctionAddOutput(qf_prolong, "output", num_comp, CEED_EVAL_NONE)); 1234 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_prolong, num_comp)); 1235 1236 CeedCall(CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_prolong)); 1237 CeedCall(CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); 1238 CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec)); 1239 CeedCall(CeedOperatorSetField(*op_prolong, "output", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1240 1241 // Set name 1242 char *prolongation_name; 1243 1244 CeedCall(CeedCalloc(18 + name_len, &prolongation_name)); 1245 sprintf(prolongation_name, "prolongation%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); 1246 CeedCall(CeedOperatorSetName(*op_prolong, prolongation_name)); 1247 CeedCall(CeedFree(&prolongation_name)); 1248 1249 // Check 1250 CeedCall(CeedOperatorCheckReady(*op_prolong)); 1251 1252 // Cleanup 1253 CeedCall(CeedQFunctionDestroy(&qf_prolong)); 1254 } 1255 1256 // Check 1257 CeedCall(CeedOperatorCheckReady(*op_coarse)); 1258 1259 // Cleanup 1260 CeedCall(CeedDestroy(&ceed)); 1261 CeedCall(CeedVectorDestroy(&mult_vec)); 1262 CeedCall(CeedElemRestrictionDestroy(&rstr_fine)); 1263 CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_fine)); 1264 CeedCall(CeedBasisDestroy(&basis_c_to_f)); 1265 return CEED_ERROR_SUCCESS; 1266 } 1267 1268 /** 1269 @brief Build 1D mass matrix and Laplacian with perturbation 1270 1271 @param[in] interp_1d Interpolation matrix in one dimension 1272 @param[in] grad_1d Gradient matrix in one dimension 1273 @param[in] q_weight_1d Quadrature weights in one dimension 1274 @param[in] P_1d Number of basis nodes in one dimension 1275 @param[in] Q_1d Number of quadrature points in one dimension 1276 @param[in] dim Dimension of basis 1277 @param[out] mass Assembled mass matrix in one dimension 1278 @param[out] laplace Assembled perturbed Laplacian in one dimension 1279 1280 @return An error code: 0 - success, otherwise - failure 1281 1282 @ref Developer 1283 **/ 1284 CeedPragmaOptimizeOff 1285 static int CeedBuildMassLaplace(const CeedScalar *interp_1d, const CeedScalar *grad_1d, const CeedScalar *q_weight_1d, CeedInt P_1d, CeedInt Q_1d, 1286 CeedInt dim, CeedScalar *mass, CeedScalar *laplace) { 1287 for (CeedInt i = 0; i < P_1d; i++) { 1288 for (CeedInt j = 0; j < P_1d; j++) { 1289 CeedScalar sum = 0.0; 1290 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]; 1291 mass[i + j * P_1d] = sum; 1292 } 1293 } 1294 // -- Laplacian 1295 for (CeedInt i = 0; i < P_1d; i++) { 1296 for (CeedInt j = 0; j < P_1d; j++) { 1297 CeedScalar sum = 0.0; 1298 1299 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]; 1300 laplace[i + j * P_1d] = sum; 1301 } 1302 } 1303 CeedScalar perturbation = dim > 2 ? 1e-6 : 1e-4; 1304 for (CeedInt i = 0; i < P_1d; i++) laplace[i + P_1d * i] += perturbation; 1305 return CEED_ERROR_SUCCESS; 1306 } 1307 CeedPragmaOptimizeOn 1308 1309 /// @} 1310 1311 /// ---------------------------------------------------------------------------- 1312 /// CeedOperator Backend API 1313 /// ---------------------------------------------------------------------------- 1314 /// @addtogroup CeedOperatorBackend 1315 /// @{ 1316 1317 /** 1318 @brief Select correct basis matrix pointer based on @ref CeedEvalMode 1319 1320 @param[in] basis `CeedBasis` from which to get the basis matrix 1321 @param[in] eval_mode Current basis evaluation mode 1322 @param[in] identity Pointer to identity matrix 1323 @param[out] basis_ptr `CeedBasis` pointer to set 1324 1325 @ref Backend 1326 **/ 1327 int CeedOperatorGetBasisPointer(CeedBasis basis, CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar **basis_ptr) { 1328 switch (eval_mode) { 1329 case CEED_EVAL_NONE: 1330 *basis_ptr = identity; 1331 break; 1332 case CEED_EVAL_INTERP: 1333 CeedCall(CeedBasisGetInterp(basis, basis_ptr)); 1334 break; 1335 case CEED_EVAL_GRAD: 1336 CeedCall(CeedBasisGetGrad(basis, basis_ptr)); 1337 break; 1338 case CEED_EVAL_DIV: 1339 CeedCall(CeedBasisGetDiv(basis, basis_ptr)); 1340 break; 1341 case CEED_EVAL_CURL: 1342 CeedCall(CeedBasisGetCurl(basis, basis_ptr)); 1343 break; 1344 case CEED_EVAL_WEIGHT: 1345 break; // Caught by QF Assembly 1346 } 1347 assert(*basis_ptr != NULL); 1348 return CEED_ERROR_SUCCESS; 1349 } 1350 1351 /** 1352 @brief Create point block restriction for active `CeedOperatorField` 1353 1354 @param[in] rstr Original `CeedElemRestriction` for active field 1355 @param[out] point_block_rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1356 1357 @return An error code: 0 - success, otherwise - failure 1358 1359 @ref Backend 1360 **/ 1361 int CeedOperatorCreateActivePointBlockRestriction(CeedElemRestriction rstr, CeedElemRestriction *point_block_rstr) { 1362 Ceed ceed; 1363 CeedInt num_elem, num_comp, shift, elem_size, comp_stride, *point_block_offsets; 1364 CeedSize l_size; 1365 const CeedInt *offsets; 1366 1367 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1368 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 1369 1370 // Expand offsets 1371 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1372 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1373 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 1374 CeedCall(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 1375 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 1376 shift = num_comp; 1377 if (comp_stride != 1) shift *= num_comp; 1378 CeedCall(CeedCalloc(num_elem * elem_size, &point_block_offsets)); 1379 for (CeedInt i = 0; i < num_elem * elem_size; i++) { 1380 point_block_offsets[i] = offsets[i] * shift; 1381 } 1382 1383 // Create new restriction 1384 CeedCall(CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp * num_comp, 1, l_size * num_comp, CEED_MEM_HOST, CEED_OWN_POINTER, 1385 point_block_offsets, point_block_rstr)); 1386 1387 // Cleanup 1388 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 1389 CeedCall(CeedDestroy(&ceed)); 1390 return CEED_ERROR_SUCCESS; 1391 } 1392 1393 /** 1394 @brief Get `CeedQFunctionAssemblyData` 1395 1396 @param[in] op `CeedOperator` to assemble 1397 @param[out] data `CeedQFunctionAssemblyData` 1398 1399 @return An error code: 0 - success, otherwise - failure 1400 1401 @ref Backend 1402 **/ 1403 int CeedOperatorGetQFunctionAssemblyData(CeedOperator op, CeedQFunctionAssemblyData *data) { 1404 if (!op->qf_assembled) { 1405 CeedQFunctionAssemblyData data; 1406 1407 CeedCall(CeedQFunctionAssemblyDataCreate(CeedOperatorReturnCeed(op), &data)); 1408 op->qf_assembled = data; 1409 } 1410 *data = op->qf_assembled; 1411 return CEED_ERROR_SUCCESS; 1412 } 1413 1414 /** 1415 @brief Create object holding `CeedQFunction` assembly data for `CeedOperator` 1416 1417 @param[in] ceed `Ceed` object used to create the `CeedQFunctionAssemblyData` 1418 @param[out] data Address of the variable where the newly created `CeedQFunctionAssemblyData` will be stored 1419 1420 @return An error code: 0 - success, otherwise - failure 1421 1422 @ref Backend 1423 **/ 1424 int CeedQFunctionAssemblyDataCreate(Ceed ceed, CeedQFunctionAssemblyData *data) { 1425 CeedCall(CeedCalloc(1, data)); 1426 (*data)->ref_count = 1; 1427 CeedCall(CeedReferenceCopy(ceed, &(*data)->ceed)); 1428 return CEED_ERROR_SUCCESS; 1429 } 1430 1431 /** 1432 @brief Increment the reference counter for a `CeedQFunctionAssemblyData` 1433 1434 @param[in,out] data `CeedQFunctionAssemblyData` to increment the reference counter 1435 1436 @return An error code: 0 - success, otherwise - failure 1437 1438 @ref Backend 1439 **/ 1440 int CeedQFunctionAssemblyDataReference(CeedQFunctionAssemblyData data) { 1441 data->ref_count++; 1442 return CEED_ERROR_SUCCESS; 1443 } 1444 1445 /** 1446 @brief Set re-use of `CeedQFunctionAssemblyData` 1447 1448 @param[in,out] data `CeedQFunctionAssemblyData` to mark for reuse 1449 @param[in] reuse_data Boolean flag indicating data re-use 1450 1451 @return An error code: 0 - success, otherwise - failure 1452 1453 @ref Backend 1454 **/ 1455 int CeedQFunctionAssemblyDataSetReuse(CeedQFunctionAssemblyData data, bool reuse_data) { 1456 data->reuse_data = reuse_data; 1457 data->needs_data_update = true; 1458 return CEED_ERROR_SUCCESS; 1459 } 1460 1461 /** 1462 @brief Mark `CeedQFunctionAssemblyData` as stale 1463 1464 @param[in,out] data `CeedQFunctionAssemblyData` to mark as stale 1465 @param[in] needs_data_update Boolean flag indicating if update is needed or completed 1466 1467 @return An error code: 0 - success, otherwise - failure 1468 1469 @ref Backend 1470 **/ 1471 int CeedQFunctionAssemblyDataSetUpdateNeeded(CeedQFunctionAssemblyData data, bool needs_data_update) { 1472 data->needs_data_update = needs_data_update; 1473 return CEED_ERROR_SUCCESS; 1474 } 1475 1476 /** 1477 @brief Determine if `CeedQFunctionAssemblyData` needs update 1478 1479 @param[in] data `CeedQFunctionAssemblyData` to mark as stale 1480 @param[out] is_update_needed Boolean flag indicating if re-assembly is required 1481 1482 @return An error code: 0 - success, otherwise - failure 1483 1484 @ref Backend 1485 **/ 1486 int CeedQFunctionAssemblyDataIsUpdateNeeded(CeedQFunctionAssemblyData data, bool *is_update_needed) { 1487 *is_update_needed = !data->reuse_data || data->needs_data_update; 1488 return CEED_ERROR_SUCCESS; 1489 } 1490 1491 /** 1492 @brief Copy the pointer to a `CeedQFunctionAssemblyData`. 1493 1494 Both pointers should be destroyed with @ref CeedQFunctionAssemblyDataDestroy(). 1495 1496 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 `CeedQFunctionAssemblyData`. 1497 This `CeedQFunctionAssemblyData` will be destroyed if ` *data_copy` is the only reference to this `CeedQFunctionAssemblyData`. 1498 1499 @param[in] data `CeedQFunctionAssemblyData` to copy reference to 1500 @param[in,out] data_copy Variable to store copied reference 1501 1502 @return An error code: 0 - success, otherwise - failure 1503 1504 @ref Backend 1505 **/ 1506 int CeedQFunctionAssemblyDataReferenceCopy(CeedQFunctionAssemblyData data, CeedQFunctionAssemblyData *data_copy) { 1507 CeedCall(CeedQFunctionAssemblyDataReference(data)); 1508 CeedCall(CeedQFunctionAssemblyDataDestroy(data_copy)); 1509 *data_copy = data; 1510 return CEED_ERROR_SUCCESS; 1511 } 1512 1513 /** 1514 @brief Get setup status for internal objects for `CeedQFunctionAssemblyData` 1515 1516 @param[in] data `CeedQFunctionAssemblyData` to retrieve status 1517 @param[out] is_setup Boolean flag for setup status 1518 1519 @return An error code: 0 - success, otherwise - failure 1520 1521 @ref Backend 1522 **/ 1523 int CeedQFunctionAssemblyDataIsSetup(CeedQFunctionAssemblyData data, bool *is_setup) { 1524 *is_setup = data->is_setup; 1525 return CEED_ERROR_SUCCESS; 1526 } 1527 1528 /** 1529 @brief Set internal objects for `CeedQFunctionAssemblyData` 1530 1531 @param[in,out] data `CeedQFunctionAssemblyData` to set objects 1532 @param[in] vec `CeedVector` to store assembled `CeedQFunction` at quadrature points 1533 @param[in] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 1534 1535 @return An error code: 0 - success, otherwise - failure 1536 1537 @ref Backend 1538 **/ 1539 int CeedQFunctionAssemblyDataSetObjects(CeedQFunctionAssemblyData data, CeedVector vec, CeedElemRestriction rstr) { 1540 CeedCall(CeedVectorReferenceCopy(vec, &data->vec)); 1541 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &data->rstr)); 1542 1543 data->is_setup = true; 1544 return CEED_ERROR_SUCCESS; 1545 } 1546 1547 /** 1548 @brief Get internal objects for `CeedQFunctionAssemblyData` 1549 1550 @param[in,out] data `CeedQFunctionAssemblyData` to set objects 1551 @param[out] vec `CeedVector` to store assembled `CeedQFunction` at quadrature points 1552 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 1553 1554 @return An error code: 0 - success, otherwise - failure 1555 1556 @ref Backend 1557 **/ 1558 int CeedQFunctionAssemblyDataGetObjects(CeedQFunctionAssemblyData data, CeedVector *vec, CeedElemRestriction *rstr) { 1559 CeedCheck(data->is_setup, data->ceed, CEED_ERROR_INCOMPLETE, "Internal objects not set; must call CeedQFunctionAssemblyDataSetObjects first."); 1560 1561 CeedCall(CeedVectorReferenceCopy(data->vec, vec)); 1562 CeedCall(CeedElemRestrictionReferenceCopy(data->rstr, rstr)); 1563 return CEED_ERROR_SUCCESS; 1564 } 1565 1566 /** 1567 @brief Destroy `CeedQFunctionAssemblyData` 1568 1569 @param[in,out] data `CeedQFunctionAssemblyData` to destroy 1570 1571 @return An error code: 0 - success, otherwise - failure 1572 1573 @ref Backend 1574 **/ 1575 int CeedQFunctionAssemblyDataDestroy(CeedQFunctionAssemblyData *data) { 1576 if (!*data || --(*data)->ref_count > 0) { 1577 *data = NULL; 1578 return CEED_ERROR_SUCCESS; 1579 } 1580 CeedCall(CeedDestroy(&(*data)->ceed)); 1581 CeedCall(CeedVectorDestroy(&(*data)->vec)); 1582 CeedCall(CeedElemRestrictionDestroy(&(*data)->rstr)); 1583 1584 CeedCall(CeedFree(data)); 1585 return CEED_ERROR_SUCCESS; 1586 } 1587 1588 /** 1589 @brief Get `CeedOperatorAssemblyData` 1590 1591 @param[in] op `CeedOperator` to assemble 1592 @param[out] data `CeedOperatorAssemblyData` 1593 1594 @return An error code: 0 - success, otherwise - failure 1595 1596 @ref Backend 1597 **/ 1598 int CeedOperatorGetOperatorAssemblyData(CeedOperator op, CeedOperatorAssemblyData *data) { 1599 if (!op->op_assembled) { 1600 CeedOperatorAssemblyData data; 1601 1602 CeedCall(CeedOperatorAssemblyDataCreate(CeedOperatorReturnCeed(op), op, &data)); 1603 op->op_assembled = data; 1604 } 1605 *data = op->op_assembled; 1606 return CEED_ERROR_SUCCESS; 1607 } 1608 1609 /** 1610 @brief Create object holding `CeedOperator` assembly data. 1611 1612 The `CeedOperatorAssemblyData` holds an array with references to every active `CeedBasis` used in the `CeedOperator`. 1613 An array with references to the corresponding active `CeedElemRestriction` is also stored. 1614 For each active `CeedBasis, the `CeedOperatorAssemblyData` holds an array of all input and output @ref CeedEvalMode for this `CeedBasis`. 1615 The `CeedOperatorAssemblyData` holds an array of offsets for indexing into the assembled `CeedQFunction` arrays to the row representing each @ref CeedEvalMode. 1616 The number of input columns across all active bases for the assembled `CeedQFunction` is also stored. 1617 Lastly, the `CeedOperatorAssembly` data holds assembled matrices representing the full action of the `CeedBasis` for all @ref CeedEvalMode. 1618 1619 @param[in] ceed `Ceed` object used to create the `CeedOperatorAssemblyData` 1620 @param[in] op `CeedOperator` to be assembled 1621 @param[out] data Address of the variable where the newly created `CeedOperatorAssemblyData` will be stored 1622 1623 @return An error code: 0 - success, otherwise - failure 1624 1625 @ref Backend 1626 **/ 1627 int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) { 1628 CeedInt num_active_bases_in = 0, num_active_bases_out = 0, offset = 0; 1629 CeedInt num_input_fields, *num_eval_modes_in = NULL, num_output_fields, *num_eval_modes_out = NULL; 1630 CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL; 1631 CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL; 1632 CeedQFunctionField *qf_fields; 1633 CeedQFunction qf; 1634 CeedOperatorField *op_fields; 1635 bool is_composite; 1636 1637 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1638 CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Can only create CeedOperator assembly data for non-composite operators."); 1639 1640 // Allocate 1641 CeedCall(CeedCalloc(1, data)); 1642 CeedCall(CeedReferenceCopy(ceed, &(*data)->ceed)); 1643 1644 // Build OperatorAssembly data 1645 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1646 1647 // Determine active input basis 1648 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL)); 1649 CeedCall(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 1650 for (CeedInt i = 0; i < num_input_fields; i++) { 1651 CeedVector vec; 1652 1653 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1654 if (vec == CEED_VECTOR_ACTIVE) { 1655 CeedInt index = -1, num_comp, q_comp; 1656 CeedEvalMode eval_mode; 1657 CeedBasis basis_in = NULL; 1658 1659 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 1660 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1661 CeedCall(CeedBasisGetNumComponents(basis_in, &num_comp)); 1662 CeedCall(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1663 for (CeedInt i = 0; i < num_active_bases_in; i++) { 1664 if ((*data)->active_bases_in[i] == basis_in) index = i; 1665 } 1666 if (index == -1) { 1667 CeedElemRestriction elem_rstr_in; 1668 1669 index = num_active_bases_in; 1670 CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_bases_in)); 1671 (*data)->active_bases_in[num_active_bases_in] = NULL; 1672 CeedCall(CeedBasisReferenceCopy(basis_in, &(*data)->active_bases_in[num_active_bases_in])); 1673 CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_elem_rstrs_in)); 1674 (*data)->active_elem_rstrs_in[num_active_bases_in] = NULL; 1675 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_in)); 1676 CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_in, &(*data)->active_elem_rstrs_in[num_active_bases_in])); 1677 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in)); 1678 CeedCall(CeedRealloc(num_active_bases_in + 1, &num_eval_modes_in)); 1679 num_eval_modes_in[index] = 0; 1680 CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_modes_in)); 1681 eval_modes_in[index] = NULL; 1682 CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_mode_offsets_in)); 1683 eval_mode_offsets_in[index] = NULL; 1684 CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->assembled_bases_in)); 1685 (*data)->assembled_bases_in[index] = NULL; 1686 num_active_bases_in++; 1687 } 1688 if (eval_mode != CEED_EVAL_WEIGHT) { 1689 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1690 CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_modes_in[index])); 1691 CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_mode_offsets_in[index])); 1692 for (CeedInt d = 0; d < q_comp; d++) { 1693 eval_modes_in[index][num_eval_modes_in[index] + d] = eval_mode; 1694 eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset; 1695 offset += num_comp; 1696 } 1697 num_eval_modes_in[index] += q_comp; 1698 } 1699 CeedCall(CeedBasisDestroy(&basis_in)); 1700 } 1701 CeedCall(CeedVectorDestroy(&vec)); 1702 } 1703 1704 // Determine active output basis 1705 CeedCall(CeedQFunctionGetFields(qf, NULL, NULL, &num_output_fields, &qf_fields)); 1706 CeedCall(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 1707 offset = 0; 1708 for (CeedInt i = 0; i < num_output_fields; i++) { 1709 CeedVector vec; 1710 1711 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1712 if (vec == CEED_VECTOR_ACTIVE) { 1713 CeedInt index = -1, num_comp, q_comp; 1714 CeedEvalMode eval_mode; 1715 CeedBasis basis_out = NULL; 1716 1717 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 1718 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1719 CeedCall(CeedBasisGetNumComponents(basis_out, &num_comp)); 1720 CeedCall(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1721 for (CeedInt i = 0; i < num_active_bases_out; i++) { 1722 if ((*data)->active_bases_out[i] == basis_out) index = i; 1723 } 1724 if (index == -1) { 1725 CeedElemRestriction elem_rstr_out; 1726 1727 index = num_active_bases_out; 1728 CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_bases_out)); 1729 (*data)->active_bases_out[num_active_bases_out] = NULL; 1730 CeedCall(CeedBasisReferenceCopy(basis_out, &(*data)->active_bases_out[num_active_bases_out])); 1731 CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_elem_rstrs_out)); 1732 (*data)->active_elem_rstrs_out[num_active_bases_out] = NULL; 1733 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_out)); 1734 CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_out, &(*data)->active_elem_rstrs_out[num_active_bases_out])); 1735 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out)); 1736 CeedCall(CeedRealloc(num_active_bases_out + 1, &num_eval_modes_out)); 1737 num_eval_modes_out[index] = 0; 1738 CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_modes_out)); 1739 eval_modes_out[index] = NULL; 1740 CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_mode_offsets_out)); 1741 eval_mode_offsets_out[index] = NULL; 1742 CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->assembled_bases_out)); 1743 (*data)->assembled_bases_out[index] = NULL; 1744 num_active_bases_out++; 1745 } 1746 if (eval_mode != CEED_EVAL_WEIGHT) { 1747 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1748 CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_modes_out[index])); 1749 CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_mode_offsets_out[index])); 1750 for (CeedInt d = 0; d < q_comp; d++) { 1751 eval_modes_out[index][num_eval_modes_out[index] + d] = eval_mode; 1752 eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset; 1753 offset += num_comp; 1754 } 1755 num_eval_modes_out[index] += q_comp; 1756 } 1757 CeedCall(CeedBasisDestroy(&basis_out)); 1758 } 1759 CeedCall(CeedVectorDestroy(&vec)); 1760 } 1761 CeedCall(CeedQFunctionDestroy(&qf)); 1762 (*data)->num_active_bases_in = num_active_bases_in; 1763 (*data)->num_eval_modes_in = num_eval_modes_in; 1764 (*data)->eval_modes_in = eval_modes_in; 1765 (*data)->eval_mode_offsets_in = eval_mode_offsets_in; 1766 (*data)->num_active_bases_out = num_active_bases_out; 1767 (*data)->num_eval_modes_out = num_eval_modes_out; 1768 (*data)->eval_modes_out = eval_modes_out; 1769 (*data)->eval_mode_offsets_out = eval_mode_offsets_out; 1770 (*data)->num_output_components = offset; 1771 return CEED_ERROR_SUCCESS; 1772 } 1773 1774 /** 1775 @brief Get `CeedOperator` @ref CeedEvalMode for assembly. 1776 1777 Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. 1778 1779 @param[in] data `CeedOperatorAssemblyData` 1780 @param[out] num_active_bases_in Total number of active bases for input 1781 @param[out] num_eval_modes_in Pointer to hold array of numbers of input @ref CeedEvalMode, or `NULL`. 1782 `eval_modes_in[0]` holds an array of eval modes for the first active `CeedBasis`. 1783 @param[out] eval_modes_in Pointer to hold arrays of input @ref CeedEvalMode, or `NULL` 1784 @param[out] eval_mode_offsets_in Pointer to hold arrays of input offsets at each quadrature point 1785 @param[out] num_active_bases_out Total number of active bases for output 1786 @param[out] num_eval_modes_out Pointer to hold array of numbers of output @ref CeedEvalMode, or `NULL` 1787 @param[out] eval_modes_out Pointer to hold arrays of output @ref CeedEvalMode, or `NULL` 1788 @param[out] eval_mode_offsets_out Pointer to hold arrays of output offsets at each quadrature point 1789 @param[out] num_output_components The number of columns in the assembled `CeedQFunction` matrix for each quadrature point, including contributions of all active bases 1790 1791 @return An error code: 0 - success, otherwise - failure 1792 1793 @ref Backend 1794 **/ 1795 int CeedOperatorAssemblyDataGetEvalModes(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedInt **num_eval_modes_in, 1796 const CeedEvalMode ***eval_modes_in, CeedSize ***eval_mode_offsets_in, CeedInt *num_active_bases_out, 1797 CeedInt **num_eval_modes_out, const CeedEvalMode ***eval_modes_out, CeedSize ***eval_mode_offsets_out, 1798 CeedSize *num_output_components) { 1799 if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in; 1800 if (num_eval_modes_in) *num_eval_modes_in = data->num_eval_modes_in; 1801 if (eval_modes_in) *eval_modes_in = (const CeedEvalMode **)data->eval_modes_in; 1802 if (eval_mode_offsets_in) *eval_mode_offsets_in = data->eval_mode_offsets_in; 1803 if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out; 1804 if (num_eval_modes_out) *num_eval_modes_out = data->num_eval_modes_out; 1805 if (eval_modes_out) *eval_modes_out = (const CeedEvalMode **)data->eval_modes_out; 1806 if (eval_mode_offsets_out) *eval_mode_offsets_out = data->eval_mode_offsets_out; 1807 if (num_output_components) *num_output_components = data->num_output_components; 1808 return CEED_ERROR_SUCCESS; 1809 } 1810 1811 /** 1812 @brief Get `CeedOperator` `CeedBasis` data for assembly. 1813 1814 Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. 1815 1816 @param[in] data `CeedOperatorAssemblyData` 1817 @param[out] num_active_bases_in Number of active input bases, or `NULL` 1818 @param[out] active_bases_in Pointer to hold active input `CeedBasis`, or `NULL` 1819 @param[out] assembled_bases_in Pointer to hold assembled active input `B` , or `NULL` 1820 @param[out] num_active_bases_out Number of active output bases, or `NULL` 1821 @param[out] active_bases_out Pointer to hold active output `CeedBasis`, or `NULL` 1822 @param[out] assembled_bases_out Pointer to hold assembled active output `B` , or `NULL` 1823 1824 @return An error code: 0 - success, otherwise - failure 1825 1826 @ref Backend 1827 **/ 1828 int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedBasis **active_bases_in, 1829 const CeedScalar ***assembled_bases_in, CeedInt *num_active_bases_out, CeedBasis **active_bases_out, 1830 const CeedScalar ***assembled_bases_out) { 1831 // Assemble B_in, B_out if needed 1832 if (assembled_bases_in && !data->assembled_bases_in[0]) { 1833 CeedInt num_qpts; 1834 1835 if (data->active_bases_in[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[0], &num_qpts)); 1836 else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_in[0], &num_qpts)); 1837 for (CeedInt b = 0; b < data->num_active_bases_in; b++) { 1838 bool has_eval_none = false; 1839 CeedInt num_nodes; 1840 CeedScalar *B_in = NULL, *identity = NULL; 1841 1842 CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[b], &num_nodes)); 1843 CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_in[b], &B_in)); 1844 1845 for (CeedInt i = 0; i < data->num_eval_modes_in[b]; i++) { 1846 has_eval_none = has_eval_none || (data->eval_modes_in[b][i] == CEED_EVAL_NONE); 1847 } 1848 if (has_eval_none) { 1849 CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); 1850 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { 1851 identity[i * num_nodes + i] = 1.0; 1852 } 1853 } 1854 1855 for (CeedInt q = 0; q < num_qpts; q++) { 1856 for (CeedInt n = 0; n < num_nodes; n++) { 1857 CeedInt d_in = 0, q_comp_in; 1858 CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; 1859 1860 for (CeedInt e_in = 0; e_in < data->num_eval_modes_in[b]; e_in++) { 1861 const CeedInt qq = data->num_eval_modes_in[b] * q; 1862 const CeedScalar *B = NULL; 1863 1864 CeedCall(CeedOperatorGetBasisPointer(data->active_bases_in[b], data->eval_modes_in[b][e_in], identity, &B)); 1865 CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_in[b], data->eval_modes_in[b][e_in], &q_comp_in)); 1866 if (q_comp_in > 1) { 1867 if (e_in == 0 || data->eval_modes_in[b][e_in] != eval_mode_in_prev) d_in = 0; 1868 else B = &B[(++d_in) * num_qpts * num_nodes]; 1869 } 1870 eval_mode_in_prev = data->eval_modes_in[b][e_in]; 1871 B_in[(qq + e_in) * num_nodes + n] = B[q * num_nodes + n]; 1872 } 1873 } 1874 } 1875 if (identity) CeedCall(CeedFree(&identity)); 1876 data->assembled_bases_in[b] = B_in; 1877 } 1878 } 1879 1880 if (assembled_bases_out && !data->assembled_bases_out[0]) { 1881 CeedInt num_qpts; 1882 1883 if (data->active_bases_out[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[0], &num_qpts)); 1884 else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_out[0], &num_qpts)); 1885 for (CeedInt b = 0; b < data->num_active_bases_out; b++) { 1886 bool has_eval_none = false; 1887 CeedInt num_nodes; 1888 CeedScalar *B_out = NULL, *identity = NULL; 1889 1890 CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[b], &num_nodes)); 1891 CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_out[b], &B_out)); 1892 1893 for (CeedInt i = 0; i < data->num_eval_modes_out[b]; i++) { 1894 has_eval_none = has_eval_none || (data->eval_modes_out[b][i] == CEED_EVAL_NONE); 1895 } 1896 if (has_eval_none) { 1897 CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); 1898 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { 1899 identity[i * num_nodes + i] = 1.0; 1900 } 1901 } 1902 1903 for (CeedInt q = 0; q < num_qpts; q++) { 1904 for (CeedInt n = 0; n < num_nodes; n++) { 1905 CeedInt d_out = 0, q_comp_out; 1906 CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; 1907 1908 for (CeedInt e_out = 0; e_out < data->num_eval_modes_out[b]; e_out++) { 1909 const CeedInt qq = data->num_eval_modes_out[b] * q; 1910 const CeedScalar *B = NULL; 1911 1912 CeedCall(CeedOperatorGetBasisPointer(data->active_bases_out[b], data->eval_modes_out[b][e_out], identity, &B)); 1913 CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_out[b], data->eval_modes_out[b][e_out], &q_comp_out)); 1914 if (q_comp_out > 1) { 1915 if (e_out == 0 || data->eval_modes_out[b][e_out] != eval_mode_out_prev) d_out = 0; 1916 else B = &B[(++d_out) * num_qpts * num_nodes]; 1917 } 1918 eval_mode_out_prev = data->eval_modes_out[b][e_out]; 1919 B_out[(qq + e_out) * num_nodes + n] = B[q * num_nodes + n]; 1920 } 1921 } 1922 } 1923 if (identity) CeedCall(CeedFree(&identity)); 1924 data->assembled_bases_out[b] = B_out; 1925 } 1926 } 1927 1928 // Pass out assembled data 1929 if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in; 1930 if (active_bases_in) *active_bases_in = data->active_bases_in; 1931 if (assembled_bases_in) *assembled_bases_in = (const CeedScalar **)data->assembled_bases_in; 1932 if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out; 1933 if (active_bases_out) *active_bases_out = data->active_bases_out; 1934 if (assembled_bases_out) *assembled_bases_out = (const CeedScalar **)data->assembled_bases_out; 1935 return CEED_ERROR_SUCCESS; 1936 } 1937 1938 /** 1939 @brief Get `CeedOperator` `CeedBasis` data for assembly. 1940 1941 Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. 1942 1943 @param[in] data `CeedOperatorAssemblyData` 1944 @param[out] num_active_elem_rstrs_in Number of active input element restrictions, or `NULL` 1945 @param[out] active_elem_rstrs_in Pointer to hold active input `CeedElemRestriction`, or `NULL` 1946 @param[out] num_active_elem_rstrs_out Number of active output element restrictions, or `NULL` 1947 @param[out] active_elem_rstrs_out Pointer to hold active output `CeedElemRestriction`, or `NULL` 1948 1949 @return An error code: 0 - success, otherwise - failure 1950 1951 @ref Backend 1952 **/ 1953 int CeedOperatorAssemblyDataGetElemRestrictions(CeedOperatorAssemblyData data, CeedInt *num_active_elem_rstrs_in, 1954 CeedElemRestriction **active_elem_rstrs_in, CeedInt *num_active_elem_rstrs_out, 1955 CeedElemRestriction **active_elem_rstrs_out) { 1956 if (num_active_elem_rstrs_in) *num_active_elem_rstrs_in = data->num_active_bases_in; 1957 if (active_elem_rstrs_in) *active_elem_rstrs_in = data->active_elem_rstrs_in; 1958 if (num_active_elem_rstrs_out) *num_active_elem_rstrs_out = data->num_active_bases_out; 1959 if (active_elem_rstrs_out) *active_elem_rstrs_out = data->active_elem_rstrs_out; 1960 return CEED_ERROR_SUCCESS; 1961 } 1962 1963 /** 1964 @brief Destroy `CeedOperatorAssemblyData` 1965 1966 @param[in,out] data `CeedOperatorAssemblyData` to destroy 1967 1968 @return An error code: 0 - success, otherwise - failure 1969 1970 @ref Backend 1971 **/ 1972 int CeedOperatorAssemblyDataDestroy(CeedOperatorAssemblyData *data) { 1973 if (!*data) { 1974 *data = NULL; 1975 return CEED_ERROR_SUCCESS; 1976 } 1977 CeedCall(CeedDestroy(&(*data)->ceed)); 1978 for (CeedInt b = 0; b < (*data)->num_active_bases_in; b++) { 1979 CeedCall(CeedBasisDestroy(&(*data)->active_bases_in[b])); 1980 CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_in[b])); 1981 CeedCall(CeedFree(&(*data)->eval_modes_in[b])); 1982 CeedCall(CeedFree(&(*data)->eval_mode_offsets_in[b])); 1983 CeedCall(CeedFree(&(*data)->assembled_bases_in[b])); 1984 } 1985 for (CeedInt b = 0; b < (*data)->num_active_bases_out; b++) { 1986 CeedCall(CeedBasisDestroy(&(*data)->active_bases_out[b])); 1987 CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_out[b])); 1988 CeedCall(CeedFree(&(*data)->eval_modes_out[b])); 1989 CeedCall(CeedFree(&(*data)->eval_mode_offsets_out[b])); 1990 CeedCall(CeedFree(&(*data)->assembled_bases_out[b])); 1991 } 1992 CeedCall(CeedFree(&(*data)->active_bases_in)); 1993 CeedCall(CeedFree(&(*data)->active_bases_out)); 1994 CeedCall(CeedFree(&(*data)->active_elem_rstrs_in)); 1995 CeedCall(CeedFree(&(*data)->active_elem_rstrs_out)); 1996 CeedCall(CeedFree(&(*data)->num_eval_modes_in)); 1997 CeedCall(CeedFree(&(*data)->num_eval_modes_out)); 1998 CeedCall(CeedFree(&(*data)->eval_modes_in)); 1999 CeedCall(CeedFree(&(*data)->eval_modes_out)); 2000 CeedCall(CeedFree(&(*data)->eval_mode_offsets_in)); 2001 CeedCall(CeedFree(&(*data)->eval_mode_offsets_out)); 2002 CeedCall(CeedFree(&(*data)->assembled_bases_in)); 2003 CeedCall(CeedFree(&(*data)->assembled_bases_out)); 2004 2005 CeedCall(CeedFree(data)); 2006 return CEED_ERROR_SUCCESS; 2007 } 2008 2009 /** 2010 @brief Retrieve fallback `CeedOperator` with a reference `Ceed` for advanced `CeedOperator` functionality 2011 2012 @param[in] op `CeedOperator` to retrieve fallback for 2013 @param[out] op_fallback Fallback `CeedOperator` 2014 2015 @return An error code: 0 - success, otherwise - failure 2016 2017 @ref Backend 2018 **/ 2019 int CeedOperatorGetFallback(CeedOperator op, CeedOperator *op_fallback) { 2020 // Create if needed 2021 if (!op->op_fallback) CeedCall(CeedOperatorCreateFallback(op)); 2022 if (op->op_fallback) { 2023 bool is_debug; 2024 Ceed ceed; 2025 2026 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2027 CeedCall(CeedIsDebug(ceed, &is_debug)); 2028 if (is_debug) { 2029 Ceed ceed_fallback; 2030 const char *resource, *resource_fallback, *op_name; 2031 2032 CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback)); 2033 CeedCall(CeedGetResource(ceed, &resource)); 2034 CeedCall(CeedGetResource(ceed_fallback, &resource_fallback)); 2035 CeedCall(CeedOperatorGetName(op, &op_name)); 2036 2037 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "---------- CeedOperator Fallback ----------\n"); 2038 CeedDebug(ceed, "Falling back from Operator with backend %s at address %p to Operator with backend %s at address %p for CeedOperator \"%s\"\n", 2039 resource, op, resource_fallback, op->op_fallback, op_name); 2040 CeedCall(CeedDestroy(&ceed_fallback)); 2041 } 2042 CeedCall(CeedDestroy(&ceed)); 2043 } 2044 *op_fallback = op->op_fallback; 2045 return CEED_ERROR_SUCCESS; 2046 } 2047 2048 /** 2049 @brief Get the parent `CeedOperator` for a fallback `CeedOperator` 2050 2051 @param[in] op `CeedOperator` context 2052 @param[out] parent Variable to store parent `CeedOperator` context 2053 2054 @return An error code: 0 - success, otherwise - failure 2055 2056 @ref Backend 2057 **/ 2058 int CeedOperatorGetFallbackParent(CeedOperator op, CeedOperator *parent) { 2059 *parent = op->op_fallback_parent ? op->op_fallback_parent : NULL; 2060 return CEED_ERROR_SUCCESS; 2061 } 2062 2063 /** 2064 @brief Get the `Ceed` context of the parent `CeedOperator` for a fallback `CeedOperator` 2065 2066 @param[in] op `CeedOperator` context 2067 @param[out] parent Variable to store parent `Ceed` context 2068 2069 @return An error code: 0 - success, otherwise - failure 2070 2071 @ref Backend 2072 **/ 2073 int CeedOperatorGetFallbackParentCeed(CeedOperator op, Ceed *parent) { 2074 *parent = NULL; 2075 if (op->op_fallback_parent) CeedCall(CeedReferenceCopy(CeedOperatorReturnCeed(op->op_fallback_parent), parent)); 2076 else CeedCall(CeedReferenceCopy(CeedOperatorReturnCeed(op), parent)); 2077 return CEED_ERROR_SUCCESS; 2078 } 2079 2080 /// @} 2081 2082 /// ---------------------------------------------------------------------------- 2083 /// CeedOperator Public API 2084 /// ---------------------------------------------------------------------------- 2085 /// @addtogroup CeedOperatorUser 2086 /// @{ 2087 2088 /** 2089 @brief Assemble a linear `CeedQFunction` associated with a `CeedOperator`. 2090 2091 This returns a `CeedVector` containing a matrix at each quadrature point providing the action of the `CeedQFunction` associated with the `CeedOperator`. 2092 The vector `assembled` is of shape `[num_elements, num_input_fields, num_output_fields, num_quad_points]` and contains column-major matrices representing the action of the `CeedQFunction` for a corresponding quadrature point on an element. 2093 2094 Inputs and outputs are in the order provided by the user when adding `CeedOperator` fields. 2095 For example, a `CeedQFunction` with inputs `u` and `gradu` and outputs `gradv` and `v` , provided in that order, would result in an assembled `CeedQFunction` that consists of `(1 + dim) x (dim + 1)` matrices at each quadrature point acting on the input ` [u, du_0, du_1]` and producing the output `[dv_0, dv_1, v]`. 2096 2097 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2098 2099 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2100 @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points 2101 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 2102 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2103 2104 @return An error code: 0 - success, otherwise - failure 2105 2106 @ref User 2107 **/ 2108 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 2109 CeedCall(CeedOperatorCheckReady(op)); 2110 2111 if (op->LinearAssembleQFunction) { 2112 // Backend version 2113 CeedCall(op->LinearAssembleQFunction(op, assembled, rstr, request)); 2114 } else { 2115 // Operator fallback 2116 CeedOperator op_fallback; 2117 2118 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleQFunction\n"); 2119 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2120 if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunction(op_fallback, assembled, rstr, request)); 2121 else return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunction"); 2122 } 2123 return CEED_ERROR_SUCCESS; 2124 } 2125 2126 /** 2127 @brief Assemble `CeedQFunction` and store result internally. 2128 2129 Return copied references of stored data to the caller. 2130 Caller is responsible for ownership and destruction of the copied references. 2131 See also @ref CeedOperatorLinearAssembleQFunction(). 2132 2133 Note: If the value of `assembled` or `rstr` passed to this function are non-`NULL` , then it is assumed that they hold valid pointers. 2134 These objects will be destroyed if `*assembled` or `*rstr` is the only reference to the object. 2135 2136 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2137 @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points 2138 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 2139 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2140 2141 @return An error code: 0 - success, otherwise - failure 2142 2143 @ref User 2144 **/ 2145 int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 2146 return CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(op, true, assembled, rstr, request); 2147 } 2148 2149 /** 2150 @brief Assemble the diagonal of a square linear `CeedOperator` 2151 2152 This overwrites a `CeedVector` with the diagonal of a linear `CeedOperator`. 2153 2154 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2155 2156 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2157 2158 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2159 @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal 2160 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2161 2162 @return An error code: 0 - success, otherwise - failure 2163 2164 @ref User 2165 **/ 2166 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2167 bool is_composite; 2168 CeedSize input_size = 0, output_size = 0; 2169 2170 CeedCall(CeedOperatorCheckReady(op)); 2171 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2172 2173 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2174 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2175 2176 // Early exit for empty operator 2177 if (!is_composite) { 2178 CeedInt num_elem = 0; 2179 2180 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2181 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2182 } 2183 2184 if (op->LinearAssembleDiagonal) { 2185 // Backend version 2186 CeedCall(op->LinearAssembleDiagonal(op, assembled, request)); 2187 return CEED_ERROR_SUCCESS; 2188 } else if (op->LinearAssembleAddDiagonal) { 2189 // Backend version with zeroing first 2190 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2191 CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); 2192 return CEED_ERROR_SUCCESS; 2193 } else if (is_composite) { 2194 // Default to summing contributions of suboperators 2195 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2196 CeedCall(CeedOperatorLinearAssembleAddDiagonalComposite(op, request, false, assembled)); 2197 return CEED_ERROR_SUCCESS; 2198 } else { 2199 // Operator fallback 2200 CeedOperator op_fallback; 2201 2202 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleDiagonal\n"); 2203 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2204 if (op_fallback) { 2205 CeedCall(CeedOperatorLinearAssembleDiagonal(op_fallback, assembled, request)); 2206 return CEED_ERROR_SUCCESS; 2207 } 2208 } 2209 // Default interface implementation 2210 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2211 CeedCall(CeedOperatorLinearAssembleAddDiagonal(op, assembled, request)); 2212 return CEED_ERROR_SUCCESS; 2213 } 2214 2215 /** 2216 @brief Assemble the diagonal of a square linear `CeedOperator`. 2217 2218 This sums into a `CeedVector` the diagonal of a linear `CeedOperator`. 2219 2220 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2221 2222 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 2223 2224 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2225 @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal 2226 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2227 2228 @return An error code: 0 - success, otherwise - failure 2229 2230 @ref User 2231 **/ 2232 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2233 bool is_composite; 2234 CeedSize input_size = 0, output_size = 0; 2235 2236 CeedCall(CeedOperatorCheckReady(op)); 2237 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2238 2239 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2240 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2241 2242 // Early exit for empty operator 2243 if (!is_composite) { 2244 CeedInt num_elem = 0; 2245 2246 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2247 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2248 } 2249 2250 if (op->LinearAssembleAddDiagonal) { 2251 // Backend version 2252 CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); 2253 return CEED_ERROR_SUCCESS; 2254 } else if (is_composite) { 2255 // Default to summing contributions of suboperators 2256 CeedCall(CeedOperatorLinearAssembleAddDiagonalComposite(op, request, false, assembled)); 2257 return CEED_ERROR_SUCCESS; 2258 } else { 2259 // Operator fallback 2260 CeedOperator op_fallback; 2261 2262 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleAddDiagonal\n"); 2263 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2264 if (op_fallback) { 2265 CeedCall(CeedOperatorLinearAssembleAddDiagonal(op_fallback, assembled, request)); 2266 return CEED_ERROR_SUCCESS; 2267 } 2268 } 2269 // Default interface implementation 2270 CeedCall(CeedOperatorLinearAssembleAddDiagonalSingle(op, request, false, assembled)); 2271 return CEED_ERROR_SUCCESS; 2272 } 2273 2274 /** 2275 @brief Fully assemble the point-block diagonal pattern of a linear `CeedOperator`. 2276 2277 Expected to be used in conjunction with @ref CeedOperatorLinearAssemblePointBlockDiagonal(). 2278 2279 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 matrix in entry `(i, j)`. 2280 Note that the `(i, j)` pairs are unique. 2281 This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemblePointBlockDiagonal() provides the values in the same ordering. 2282 2283 This will generally be slow unless your operator is low-order. 2284 2285 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2286 2287 @param[in] op `CeedOperator` to assemble 2288 @param[out] num_entries Number of entries in coordinate nonzero pattern 2289 @param[out] rows Row number for each entry 2290 @param[out] cols Column number for each entry 2291 2292 @ref User 2293 **/ 2294 int CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { 2295 bool is_composite; 2296 CeedInt num_active_components, num_sub_operators; 2297 CeedOperator *sub_operators; 2298 2299 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2300 2301 CeedSize input_size = 0, output_size = 0; 2302 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2303 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2304 2305 if (is_composite) { 2306 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub_operators)); 2307 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2308 } else { 2309 sub_operators = &op; 2310 num_sub_operators = 1; 2311 } 2312 2313 // Verify operator can be assembled correctly 2314 { 2315 CeedOperatorAssemblyData data; 2316 CeedInt num_active_elem_rstrs, comp_stride; 2317 CeedElemRestriction *active_elem_rstrs; 2318 2319 // Get initial values to check against 2320 CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[0], &data)); 2321 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL)); 2322 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[0], &comp_stride)); 2323 CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[0], &num_active_components)); 2324 2325 // Verify that all active element restrictions have same component stride and number of components 2326 for (CeedInt k = 0; k < num_sub_operators; k++) { 2327 CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[k], &data)); 2328 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL)); 2329 for (CeedInt i = 0; i < num_active_elem_rstrs; i++) { 2330 CeedInt comp_stride_sub, num_active_components_sub; 2331 2332 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[i], &comp_stride_sub)); 2333 CeedCheck(comp_stride == comp_stride_sub, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, 2334 "Active element restrictions must have the same component stride: %d vs %d", comp_stride, comp_stride_sub); 2335 CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[i], &num_active_components_sub)); 2336 CeedCheck(num_active_components == num_active_components_sub, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, 2337 "All suboperators must have the same number of output components." 2338 " Previous: %" CeedInt_FMT " Current: %" CeedInt_FMT, 2339 num_active_components, num_active_components_sub); 2340 } 2341 } 2342 } 2343 *num_entries = input_size * num_active_components; 2344 CeedCall(CeedCalloc(*num_entries, rows)); 2345 CeedCall(CeedCalloc(*num_entries, cols)); 2346 2347 for (CeedInt o = 0; o < num_sub_operators; o++) { 2348 CeedElemRestriction active_elem_rstr, point_block_active_elem_rstr; 2349 CeedInt comp_stride, num_elem, elem_size; 2350 const CeedInt *offsets, *point_block_offsets; 2351 2352 CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[o], &active_elem_rstr)); 2353 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstr, &comp_stride)); 2354 CeedCall(CeedElemRestrictionGetNumElements(active_elem_rstr, &num_elem)); 2355 CeedCall(CeedElemRestrictionGetElementSize(active_elem_rstr, &elem_size)); 2356 CeedCall(CeedElemRestrictionGetOffsets(active_elem_rstr, CEED_MEM_HOST, &offsets)); 2357 2358 CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstr, &point_block_active_elem_rstr)); 2359 CeedCall(CeedElemRestrictionGetOffsets(point_block_active_elem_rstr, CEED_MEM_HOST, &point_block_offsets)); 2360 2361 for (CeedSize i = 0; i < num_elem * elem_size; i++) { 2362 for (CeedInt c_out = 0; c_out < num_active_components; c_out++) { 2363 for (CeedInt c_in = 0; c_in < num_active_components; c_in++) { 2364 (*rows)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_out * comp_stride; 2365 (*cols)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_in * comp_stride; 2366 } 2367 } 2368 } 2369 2370 CeedCall(CeedElemRestrictionRestoreOffsets(active_elem_rstr, &offsets)); 2371 CeedCall(CeedElemRestrictionRestoreOffsets(point_block_active_elem_rstr, &point_block_offsets)); 2372 CeedCall(CeedElemRestrictionDestroy(&active_elem_rstr)); 2373 CeedCall(CeedElemRestrictionDestroy(&point_block_active_elem_rstr)); 2374 } 2375 return CEED_ERROR_SUCCESS; 2376 } 2377 2378 /** 2379 @brief Assemble the point block diagonal of a square linear `CeedOperator`. 2380 2381 This overwrites a `CeedVector` with the point block diagonal of a linear `CeedOperator`. 2382 2383 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2384 2385 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2386 2387 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2388 @param[out] assembled `CeedVector` to store assembled `CeedOperator` point block diagonal, provided in row-major form with an `num_comp * num_comp` block at each node. 2389 The dimensions of this vector are derived from the active vector for the `CeedOperator`. 2390 The array has shape `[nodes, component out, component in]`. 2391 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2392 2393 @return An error code: 0 - success, otherwise - failure 2394 2395 @ref User 2396 **/ 2397 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2398 bool is_composite; 2399 CeedSize input_size = 0, output_size = 0; 2400 2401 CeedCall(CeedOperatorCheckReady(op)); 2402 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2403 2404 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2405 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2406 2407 // Early exit for empty operator 2408 if (!is_composite) { 2409 CeedInt num_elem = 0; 2410 2411 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2412 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2413 } 2414 2415 if (op->LinearAssemblePointBlockDiagonal) { 2416 // Backend version 2417 CeedCall(op->LinearAssemblePointBlockDiagonal(op, assembled, request)); 2418 return CEED_ERROR_SUCCESS; 2419 } else if (op->LinearAssembleAddPointBlockDiagonal) { 2420 // Backend version with zeroing first 2421 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2422 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); 2423 return CEED_ERROR_SUCCESS; 2424 } else { 2425 // Operator fallback 2426 CeedOperator op_fallback; 2427 2428 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssemblePointBlockDiagonal\n"); 2429 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2430 if (op_fallback) { 2431 CeedCall(CeedOperatorLinearAssemblePointBlockDiagonal(op_fallback, assembled, request)); 2432 return CEED_ERROR_SUCCESS; 2433 } 2434 } 2435 // Default interface implementation 2436 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2437 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); 2438 return CEED_ERROR_SUCCESS; 2439 } 2440 2441 /** 2442 @brief Assemble the point block diagonal of a square linear `CeedOperator`. 2443 2444 This sums into a `CeedVector` with the point block diagonal of a linear `CeedOperator`. 2445 2446 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2447 2448 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2449 2450 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2451 @param[out] assembled `CeedVector` to store assembled CeedOperator point block diagonal, provided in row-major form with an `num_comp * num_comp` block at each node. 2452 The dimensions of this vector are derived from the active vector for the `CeedOperator`. 2453 The array has shape `[nodes, component out, component in]`. 2454 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2455 2456 @return An error code: 0 - success, otherwise - failure 2457 2458 @ref User 2459 **/ 2460 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2461 bool is_composite; 2462 CeedSize input_size = 0, output_size = 0; 2463 2464 CeedCall(CeedOperatorCheckReady(op)); 2465 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2466 2467 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2468 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2469 2470 // Early exit for empty operator 2471 if (!is_composite) { 2472 CeedInt num_elem = 0; 2473 2474 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2475 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2476 } 2477 2478 if (op->LinearAssembleAddPointBlockDiagonal) { 2479 // Backend version 2480 CeedCall(op->LinearAssembleAddPointBlockDiagonal(op, assembled, request)); 2481 return CEED_ERROR_SUCCESS; 2482 } else { 2483 // Operator fallback 2484 CeedOperator op_fallback; 2485 2486 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleAddPointBlockDiagonal\n"); 2487 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2488 if (op_fallback) { 2489 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op_fallback, assembled, request)); 2490 return CEED_ERROR_SUCCESS; 2491 } 2492 } 2493 // Default interface implementation 2494 if (is_composite) { 2495 CeedCall(CeedOperatorLinearAssembleAddDiagonalComposite(op, request, true, assembled)); 2496 } else { 2497 CeedCall(CeedOperatorLinearAssembleAddDiagonalSingle(op, request, true, assembled)); 2498 } 2499 return CEED_ERROR_SUCCESS; 2500 } 2501 2502 /** 2503 @brief Fully assemble the nonzero pattern of a linear `CeedOperator`. 2504 2505 Expected to be used in conjunction with @ref CeedOperatorLinearAssemble(). 2506 2507 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 matrix in entry `(i, j)`. 2508 Note that the `(i, j)` pairs are not unique and may repeat. 2509 This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemble() provides the values in the same ordering. 2510 2511 This will generally be slow unless your operator is low-order. 2512 2513 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2514 2515 @param[in] op `CeedOperator` to assemble 2516 @param[out] num_entries Number of entries in coordinate nonzero pattern 2517 @param[out] rows Row number for each entry 2518 @param[out] cols Column number for each entry 2519 2520 @ref User 2521 **/ 2522 int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { 2523 bool is_composite; 2524 CeedInt num_suboperators, offset = 0; 2525 CeedSize single_entries; 2526 CeedOperator *sub_operators; 2527 2528 CeedCall(CeedOperatorCheckReady(op)); 2529 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2530 2531 if (op->LinearAssembleSymbolic) { 2532 // Backend version 2533 CeedCall(op->LinearAssembleSymbolic(op, num_entries, rows, cols)); 2534 return CEED_ERROR_SUCCESS; 2535 } else { 2536 // Operator fallback 2537 CeedOperator op_fallback; 2538 2539 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleSymbolic\n"); 2540 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2541 if (op_fallback) { 2542 CeedCall(CeedOperatorLinearAssembleSymbolic(op_fallback, num_entries, rows, cols)); 2543 return CEED_ERROR_SUCCESS; 2544 } 2545 } 2546 2547 // Default interface implementation 2548 2549 // Count entries and allocate rows, cols arrays 2550 CeedCall(CeedOperatorLinearAssembleGetNumEntries(op, num_entries)); 2551 CeedCall(CeedCalloc(*num_entries, rows)); 2552 CeedCall(CeedCalloc(*num_entries, cols)); 2553 2554 // Assemble nonzero locations 2555 if (is_composite) { 2556 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2557 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2558 for (CeedInt k = 0; k < num_suboperators; ++k) { 2559 CeedCall(CeedOperatorAssembleSymbolicSingle(sub_operators[k], offset, *rows, *cols)); 2560 CeedCall(CeedOperatorAssemblyCountEntriesSingle(sub_operators[k], &single_entries)); 2561 offset += single_entries; 2562 } 2563 } else { 2564 CeedCall(CeedOperatorAssembleSymbolicSingle(op, offset, *rows, *cols)); 2565 } 2566 return CEED_ERROR_SUCCESS; 2567 } 2568 2569 /** 2570 @brief Fully assemble the nonzero entries of a linear operator. 2571 2572 Expected to be used in conjunction with @ref CeedOperatorLinearAssembleSymbolic(). 2573 2574 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 matrix in entry `(i, j)`. 2575 Note that the `(i, j)` pairs are not unique and may repeat. 2576 This function returns the values of the nonzero entries to be added, their `(i, j)` locations are provided by @ref CeedOperatorLinearAssembleSymbolic(). 2577 2578 This will generally be slow unless your operator is low-order. 2579 2580 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2581 2582 @param[in] op `CeedOperator` to assemble 2583 @param[out] values Values to assemble into matrix 2584 2585 @ref User 2586 **/ 2587 int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { 2588 bool is_composite; 2589 CeedInt num_suboperators, offset = 0; 2590 CeedSize single_entries = 0; 2591 CeedOperator *sub_operators; 2592 2593 CeedCall(CeedOperatorCheckReady(op)); 2594 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2595 2596 // Early exit for empty operator 2597 if (!is_composite) { 2598 CeedInt num_elem = 0; 2599 2600 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2601 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2602 } 2603 2604 if (op->LinearAssemble) { 2605 // Backend version 2606 CeedCall(op->LinearAssemble(op, values)); 2607 return CEED_ERROR_SUCCESS; 2608 } else if (is_composite) { 2609 // Default to summing contributions of suboperators 2610 CeedCall(CeedVectorSetValue(values, 0.0)); 2611 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2612 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2613 for (CeedInt k = 0; k < num_suboperators; k++) { 2614 CeedCall(CeedOperatorAssembleSingle(sub_operators[k], offset, values)); 2615 CeedCall(CeedOperatorAssemblyCountEntriesSingle(sub_operators[k], &single_entries)); 2616 offset += single_entries; 2617 } 2618 return CEED_ERROR_SUCCESS; 2619 } else if (op->LinearAssembleSingle) { 2620 CeedCall(CeedVectorSetValue(values, 0.0)); 2621 CeedCall(CeedOperatorAssembleSingle(op, offset, values)); 2622 return CEED_ERROR_SUCCESS; 2623 } else { 2624 // Operator fallback 2625 CeedOperator op_fallback; 2626 2627 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssemble\n"); 2628 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2629 if (op_fallback) { 2630 CeedCall(CeedOperatorLinearAssemble(op_fallback, values)); 2631 return CEED_ERROR_SUCCESS; 2632 } 2633 } 2634 2635 // Default to interface version if non-composite and no fallback 2636 CeedCall(CeedVectorSetValue(values, 0.0)); 2637 CeedCall(CeedOperatorAssembleSingle(op, offset, values)); 2638 return CEED_ERROR_SUCCESS; 2639 } 2640 2641 /** 2642 @brief Get the multiplicity of nodes across sub-operators in a composite `CeedOperator`. 2643 2644 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2645 2646 @param[in] op Composite `CeedOperator` 2647 @param[in] num_skip_indices Number of sub-operators to skip 2648 @param[in] skip_indices Array of indices of sub-operators to skip 2649 @param[out] mult Vector to store multiplicity (of size `l_size` ) 2650 2651 @return An error code: 0 - success, otherwise - failure 2652 2653 @ref User 2654 **/ 2655 int CeedOperatorCompositeGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult) { 2656 Ceed ceed; 2657 CeedInt num_suboperators; 2658 CeedSize l_vec_len; 2659 CeedScalar *mult_array; 2660 CeedVector ones_l_vec; 2661 CeedElemRestriction elem_rstr, mult_elem_rstr; 2662 CeedOperator *sub_operators; 2663 2664 CeedCall(CeedOperatorCheckReady(op)); 2665 2666 // Zero mult vector 2667 CeedCall(CeedVectorSetValue(mult, 0.0)); 2668 2669 // Get suboperators 2670 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2671 if (num_suboperators == 0) return CEED_ERROR_SUCCESS; 2672 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2673 2674 // Work vector 2675 CeedCall(CeedVectorGetLength(mult, &l_vec_len)); 2676 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2677 CeedCall(CeedVectorCreate(ceed, l_vec_len, &ones_l_vec)); 2678 CeedCall(CeedDestroy(&ceed)); 2679 CeedCall(CeedVectorSetValue(ones_l_vec, 1.0)); 2680 CeedCall(CeedVectorGetArray(mult, CEED_MEM_HOST, &mult_array)); 2681 2682 // Compute multiplicity across suboperators 2683 for (CeedInt i = 0; i < num_suboperators; i++) { 2684 const CeedScalar *sub_mult_array; 2685 CeedVector sub_mult_l_vec, ones_e_vec; 2686 2687 // -- Check for suboperator to skip 2688 for (CeedInt j = 0; j < num_skip_indices; j++) { 2689 if (skip_indices[j] == i) continue; 2690 } 2691 2692 // -- Sub operator multiplicity 2693 CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[i], &elem_rstr)); 2694 CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr, &mult_elem_rstr)); 2695 CeedCall(CeedElemRestrictionDestroy(&elem_rstr)); 2696 CeedCall(CeedElemRestrictionCreateVector(mult_elem_rstr, &sub_mult_l_vec, &ones_e_vec)); 2697 CeedCall(CeedVectorSetValue(sub_mult_l_vec, 0.0)); 2698 CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_NOTRANSPOSE, ones_l_vec, ones_e_vec, CEED_REQUEST_IMMEDIATE)); 2699 CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_TRANSPOSE, ones_e_vec, sub_mult_l_vec, CEED_REQUEST_IMMEDIATE)); 2700 CeedCall(CeedVectorGetArrayRead(sub_mult_l_vec, CEED_MEM_HOST, &sub_mult_array)); 2701 // ---- Flag every node present in the current suboperator 2702 for (CeedSize j = 0; j < l_vec_len; j++) { 2703 if (sub_mult_array[j] > 0.0) mult_array[j] += 1.0; 2704 } 2705 CeedCall(CeedVectorRestoreArrayRead(sub_mult_l_vec, &sub_mult_array)); 2706 CeedCall(CeedVectorDestroy(&sub_mult_l_vec)); 2707 CeedCall(CeedVectorDestroy(&ones_e_vec)); 2708 CeedCall(CeedElemRestrictionDestroy(&mult_elem_rstr)); 2709 } 2710 CeedCall(CeedVectorRestoreArray(mult, &mult_array)); 2711 CeedCall(CeedVectorDestroy(&ones_l_vec)); 2712 return CEED_ERROR_SUCCESS; 2713 } 2714 2715 /** 2716 @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator`, creating the prolongation basis from the fine and coarse grid interpolation. 2717 2718 Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. 2719 2720 @param[in] op_fine Fine grid `CeedOperator` 2721 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 2722 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 2723 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 2724 @param[out] op_coarse Coarse grid `CeedOperator` 2725 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 2726 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 2727 2728 @return An error code: 0 - success, otherwise - failure 2729 2730 @ref User 2731 **/ 2732 int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2733 CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { 2734 CeedBasis basis_c_to_f = NULL; 2735 2736 CeedCall(CeedOperatorCheckReady(op_fine)); 2737 2738 // Build prolongation matrix, if required 2739 if (op_prolong || op_restrict) { 2740 CeedBasis basis_fine; 2741 2742 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2743 CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f)); 2744 CeedCall(CeedBasisDestroy(&basis_fine)); 2745 } 2746 2747 // Core code 2748 CeedCall(CeedOperatorMultigridLevelCreateSingle_Core(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, 2749 op_restrict)); 2750 return CEED_ERROR_SUCCESS; 2751 } 2752 2753 /** 2754 @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a tensor basis for the active basis. 2755 2756 Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. 2757 2758 @param[in] op_fine Fine grid `CeedOperator` 2759 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 2760 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 2761 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 2762 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator` 2763 @param[out] op_coarse Coarse grid `CeedOperator` 2764 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 2765 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 2766 2767 @return An error code: 0 - success, otherwise - failure 2768 2769 @ref User 2770 **/ 2771 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2772 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 2773 CeedOperator *op_restrict) { 2774 Ceed ceed; 2775 CeedInt Q_f, Q_c; 2776 CeedBasis basis_fine, basis_c_to_f = NULL; 2777 2778 CeedCall(CeedOperatorCheckReady(op_fine)); 2779 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 2780 2781 // Check for compatible quadrature spaces 2782 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2783 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); 2784 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); 2785 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, 2786 "Bases must have compatible quadrature spaces." 2787 " Fine grid: %" CeedInt_FMT " points, Coarse grid: %" CeedInt_FMT " points", 2788 Q_f, Q_c); 2789 2790 // Create coarse to fine basis, if required 2791 if (op_prolong || op_restrict) { 2792 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c; 2793 CeedScalar *q_ref, *q_weight, *grad; 2794 2795 // Check if interpolation matrix is provided 2796 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE, 2797 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix"); 2798 CeedCall(CeedBasisGetDimension(basis_fine, &dim)); 2799 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); 2800 CeedCall(CeedBasisGetNumNodes1D(basis_fine, &P_1d_f)); 2801 CeedCall(CeedBasisDestroy(&basis_fine)); 2802 CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); 2803 P_1d_c = dim == 1 ? num_nodes_c : dim == 2 ? sqrt(num_nodes_c) : cbrt(num_nodes_c); 2804 CeedCall(CeedCalloc(P_1d_f, &q_ref)); 2805 CeedCall(CeedCalloc(P_1d_f, &q_weight)); 2806 CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad)); 2807 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)); 2808 CeedCall(CeedFree(&q_ref)); 2809 CeedCall(CeedFree(&q_weight)); 2810 CeedCall(CeedFree(&grad)); 2811 } 2812 2813 // Core code 2814 CeedCall(CeedOperatorMultigridLevelCreateSingle_Core(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, 2815 op_restrict)); 2816 CeedCall(CeedDestroy(&ceed)); 2817 return CEED_ERROR_SUCCESS; 2818 } 2819 2820 /** 2821 @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a non-tensor basis for the active vector 2822 2823 Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. 2824 2825 @param[in] op_fine Fine grid `CeedOperator` 2826 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 2827 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 2828 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 2829 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator` 2830 @param[out] op_coarse Coarse grid `CeedOperator` 2831 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 2832 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 2833 2834 @return An error code: 0 - success, otherwise - failure 2835 2836 @ref User 2837 **/ 2838 int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2839 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 2840 CeedOperator *op_restrict) { 2841 Ceed ceed; 2842 CeedInt Q_f, Q_c; 2843 CeedBasis basis_fine, basis_c_to_f = NULL; 2844 2845 CeedCall(CeedOperatorCheckReady(op_fine)); 2846 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 2847 2848 // Check for compatible quadrature spaces 2849 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2850 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); 2851 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); 2852 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces"); 2853 2854 // Coarse to fine basis 2855 if (op_prolong || op_restrict) { 2856 CeedInt dim, num_comp, num_nodes_c, num_nodes_f; 2857 CeedScalar *q_ref, *q_weight, *grad; 2858 CeedElemTopology topo; 2859 2860 // Check if interpolation matrix is provided 2861 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE, 2862 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix"); 2863 CeedCall(CeedBasisGetTopology(basis_fine, &topo)); 2864 CeedCall(CeedBasisGetDimension(basis_fine, &dim)); 2865 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); 2866 CeedCall(CeedBasisGetNumNodes(basis_fine, &num_nodes_f)); 2867 CeedCall(CeedBasisDestroy(&basis_fine)); 2868 CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); 2869 CeedCall(CeedCalloc(num_nodes_f * dim, &q_ref)); 2870 CeedCall(CeedCalloc(num_nodes_f, &q_weight)); 2871 CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad)); 2872 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)); 2873 CeedCall(CeedFree(&q_ref)); 2874 CeedCall(CeedFree(&q_weight)); 2875 CeedCall(CeedFree(&grad)); 2876 } 2877 2878 // Core code 2879 CeedCall(CeedOperatorMultigridLevelCreateSingle_Core(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, 2880 op_restrict)); 2881 CeedCall(CeedDestroy(&ceed)); 2882 return CEED_ERROR_SUCCESS; 2883 } 2884 2885 /** 2886 @brief Build a FDM based approximate inverse for each element for a `CeedOperator`. 2887 2888 This returns a `CeedOperator` and `CeedVector` to apply a Fast Diagonalization Method based approximate inverse. 2889 This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, \f$M = V^T V, K = V^T S V\f$. 2890 The assembled `CeedQFunction` is used to modify the eigenvalues from simultaneous diagonalization and obtain an approximate inverse of the form \f$V^T \hat S V\f$. 2891 The `CeedOperator` must be linear and non-composite. 2892 The associated `CeedQFunction` must therefore also be linear. 2893 2894 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2895 2896 @param[in] op `CeedOperator` to create element inverses 2897 @param[out] fdm_inv `CeedOperator` to apply the action of a FDM based inverse for each element 2898 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2899 2900 @return An error code: 0 - success, otherwise - failure 2901 2902 @ref User 2903 **/ 2904 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, CeedRequest *request) { 2905 Ceed ceed, ceed_parent; 2906 bool interp = false, grad = false, is_tensor_basis = true; 2907 CeedInt num_input_fields, P_1d, Q_1d, num_nodes, num_qpts, dim, num_comp = 1, num_elem = 1; 2908 CeedScalar *mass, *laplace, *x, *fdm_interp, *lambda, *elem_avg; 2909 const CeedScalar *interp_1d, *grad_1d, *q_weight_1d; 2910 CeedVector q_data; 2911 CeedElemRestriction rstr = NULL, rstr_qd_i; 2912 CeedBasis basis = NULL, fdm_basis; 2913 CeedQFunctionContext ctx_fdm; 2914 CeedQFunctionField *qf_fields; 2915 CeedQFunction qf, qf_fdm; 2916 CeedOperatorField *op_fields; 2917 2918 CeedCall(CeedOperatorCheckReady(op)); 2919 2920 if (op->CreateFDMElementInverse) { 2921 // Backend version 2922 CeedCall(op->CreateFDMElementInverse(op, fdm_inv, request)); 2923 return CEED_ERROR_SUCCESS; 2924 } else { 2925 // Operator fallback 2926 CeedOperator op_fallback; 2927 2928 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorCreateFDMElementInverse\n"); 2929 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2930 if (op_fallback) { 2931 CeedCall(CeedOperatorCreateFDMElementInverse(op_fallback, fdm_inv, request)); 2932 return CEED_ERROR_SUCCESS; 2933 } 2934 } 2935 2936 // Default interface implementation 2937 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2938 CeedCall(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 2939 CeedCall(CeedOperatorGetQFunction(op, &qf)); 2940 2941 // Determine active input basis 2942 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_fields, NULL, NULL)); 2943 CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 2944 for (CeedInt i = 0; i < num_input_fields; i++) { 2945 CeedVector vec; 2946 2947 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 2948 if (vec == CEED_VECTOR_ACTIVE) { 2949 CeedEvalMode eval_mode; 2950 2951 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 2952 interp = interp || eval_mode == CEED_EVAL_INTERP; 2953 grad = grad || eval_mode == CEED_EVAL_GRAD; 2954 if (!basis) CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 2955 if (!rstr) CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr)); 2956 } 2957 CeedCall(CeedVectorDestroy(&vec)); 2958 } 2959 CeedCheck(basis, ceed, CEED_ERROR_BACKEND, "No active field set"); 2960 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 2961 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 2962 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 2963 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 2964 CeedCall(CeedBasisGetDimension(basis, &dim)); 2965 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 2966 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 2967 2968 // Build and diagonalize 1D Mass and Laplacian 2969 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 2970 CeedCheck(is_tensor_basis, ceed, CEED_ERROR_BACKEND, "FDMElementInverse only supported for tensor bases"); 2971 CeedCall(CeedCalloc(P_1d * P_1d, &mass)); 2972 CeedCall(CeedCalloc(P_1d * P_1d, &laplace)); 2973 CeedCall(CeedCalloc(P_1d * P_1d, &x)); 2974 CeedCall(CeedCalloc(P_1d * P_1d, &fdm_interp)); 2975 CeedCall(CeedCalloc(P_1d, &lambda)); 2976 // -- Build matrices 2977 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 2978 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 2979 CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d)); 2980 CeedCall(CeedBuildMassLaplace(interp_1d, grad_1d, q_weight_1d, P_1d, Q_1d, dim, mass, laplace)); 2981 2982 // -- Diagonalize 2983 CeedCall(CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d)); 2984 CeedCall(CeedFree(&mass)); 2985 CeedCall(CeedFree(&laplace)); 2986 for (CeedInt i = 0; i < P_1d; i++) { 2987 for (CeedInt j = 0; j < P_1d; j++) fdm_interp[i + j * P_1d] = x[j + i * P_1d]; 2988 } 2989 CeedCall(CeedFree(&x)); 2990 2991 { 2992 CeedInt layout[3], num_modes = (interp ? 1 : 0) + (grad ? dim : 0); 2993 CeedScalar max_norm = 0; 2994 const CeedScalar *assembled_array, *q_weight_array; 2995 CeedVector assembled = NULL, q_weight; 2996 CeedElemRestriction rstr_qf = NULL; 2997 2998 // Assemble QFunction 2999 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled, &rstr_qf, request)); 3000 CeedCall(CeedElemRestrictionGetELayout(rstr_qf, layout)); 3001 CeedCall(CeedElemRestrictionDestroy(&rstr_qf)); 3002 CeedCall(CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm)); 3003 3004 // Calculate element averages 3005 CeedCall(CeedVectorCreate(ceed_parent, num_qpts, &q_weight)); 3006 CeedCall(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_weight)); 3007 CeedCall(CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array)); 3008 CeedCall(CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array)); 3009 CeedCall(CeedCalloc(num_elem, &elem_avg)); 3010 const CeedScalar qf_value_bound = max_norm * 100 * CEED_EPSILON; 3011 3012 for (CeedInt e = 0; e < num_elem; e++) { 3013 CeedInt count = 0; 3014 3015 for (CeedInt q = 0; q < num_qpts; q++) { 3016 for (CeedInt i = 0; i < num_comp * num_comp * num_modes * num_modes; i++) { 3017 if (fabs(assembled_array[q * layout[0] + i * layout[1] + e * layout[2]]) > qf_value_bound) { 3018 elem_avg[e] += assembled_array[q * layout[0] + i * layout[1] + e * layout[2]] / q_weight_array[q]; 3019 count++; 3020 } 3021 } 3022 } 3023 if (count) { 3024 elem_avg[e] /= count; 3025 } else { 3026 elem_avg[e] = 1.0; 3027 } 3028 } 3029 CeedCall(CeedVectorRestoreArrayRead(assembled, &assembled_array)); 3030 CeedCall(CeedVectorDestroy(&assembled)); 3031 CeedCall(CeedVectorRestoreArrayRead(q_weight, &q_weight_array)); 3032 CeedCall(CeedVectorDestroy(&q_weight)); 3033 } 3034 3035 // Build FDM diagonal 3036 { 3037 CeedScalar *q_data_array, *fdm_diagonal; 3038 3039 CeedCall(CeedCalloc(num_comp * num_nodes, &fdm_diagonal)); 3040 const CeedScalar fdm_diagonal_bound = num_nodes * CEED_EPSILON; 3041 for (CeedInt c = 0; c < num_comp; c++) { 3042 for (CeedInt n = 0; n < num_nodes; n++) { 3043 if (interp) fdm_diagonal[c * num_nodes + n] = 1.0; 3044 if (grad) { 3045 for (CeedInt d = 0; d < dim; d++) { 3046 CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d; 3047 fdm_diagonal[c * num_nodes + n] += lambda[i]; 3048 } 3049 } 3050 if (fabs(fdm_diagonal[c * num_nodes + n]) < fdm_diagonal_bound) fdm_diagonal[c * num_nodes + n] = fdm_diagonal_bound; 3051 } 3052 } 3053 CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * num_nodes, &q_data)); 3054 CeedCall(CeedVectorSetValue(q_data, 0.0)); 3055 CeedCall(CeedVectorGetArrayWrite(q_data, CEED_MEM_HOST, &q_data_array)); 3056 for (CeedInt e = 0; e < num_elem; e++) { 3057 for (CeedInt c = 0; c < num_comp; c++) { 3058 for (CeedInt n = 0; n < num_nodes; n++) { 3059 q_data_array[(e * num_comp + c) * num_nodes + n] = 1. / (elem_avg[e] * fdm_diagonal[c * num_nodes + n]); 3060 } 3061 } 3062 } 3063 CeedCall(CeedFree(&elem_avg)); 3064 CeedCall(CeedFree(&fdm_diagonal)); 3065 CeedCall(CeedVectorRestoreArray(q_data, &q_data_array)); 3066 } 3067 3068 // Setup FDM operator 3069 // -- Basis 3070 { 3071 CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy; 3072 3073 CeedCall(CeedCalloc(P_1d * P_1d, &grad_dummy)); 3074 CeedCall(CeedCalloc(P_1d, &q_ref_dummy)); 3075 CeedCall(CeedCalloc(P_1d, &q_weight_dummy)); 3076 CeedCall(CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, fdm_interp, grad_dummy, q_ref_dummy, q_weight_dummy, &fdm_basis)); 3077 CeedCall(CeedFree(&fdm_interp)); 3078 CeedCall(CeedFree(&grad_dummy)); 3079 CeedCall(CeedFree(&q_ref_dummy)); 3080 CeedCall(CeedFree(&q_weight_dummy)); 3081 CeedCall(CeedFree(&lambda)); 3082 } 3083 3084 // -- Restriction 3085 { 3086 CeedInt strides[3] = {1, num_nodes, num_nodes * num_comp}; 3087 CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, num_nodes, num_comp, 3088 (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes, strides, &rstr_qd_i)); 3089 } 3090 3091 // -- QFunction 3092 CeedCall(CeedQFunctionCreateInteriorByName(ceed_parent, "Scale", &qf_fdm)); 3093 CeedCall(CeedQFunctionAddInput(qf_fdm, "input", num_comp, CEED_EVAL_INTERP)); 3094 CeedCall(CeedQFunctionAddInput(qf_fdm, "scale", num_comp, CEED_EVAL_NONE)); 3095 CeedCall(CeedQFunctionAddOutput(qf_fdm, "output", num_comp, CEED_EVAL_INTERP)); 3096 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_fdm, num_comp)); 3097 3098 // -- QFunction context 3099 { 3100 CeedInt *num_comp_data; 3101 3102 CeedCall(CeedCalloc(1, &num_comp_data)); 3103 num_comp_data[0] = num_comp; 3104 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_fdm)); 3105 CeedCall(CeedQFunctionContextSetData(ctx_fdm, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_data), num_comp_data)); 3106 } 3107 CeedCall(CeedQFunctionSetContext(qf_fdm, ctx_fdm)); 3108 CeedCall(CeedQFunctionContextDestroy(&ctx_fdm)); 3109 3110 // -- Operator 3111 CeedCall(CeedOperatorCreate(ceed_parent, qf_fdm, NULL, NULL, fdm_inv)); 3112 CeedCall(CeedOperatorSetField(*fdm_inv, "input", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); 3113 CeedCall(CeedOperatorSetField(*fdm_inv, "scale", rstr_qd_i, CEED_BASIS_NONE, q_data)); 3114 CeedCall(CeedOperatorSetField(*fdm_inv, "output", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); 3115 3116 // Cleanup 3117 CeedCall(CeedDestroy(&ceed)); 3118 CeedCall(CeedDestroy(&ceed_parent)); 3119 CeedCall(CeedVectorDestroy(&q_data)); 3120 CeedCall(CeedElemRestrictionDestroy(&rstr)); 3121 CeedCall(CeedElemRestrictionDestroy(&rstr_qd_i)); 3122 CeedCall(CeedBasisDestroy(&basis)); 3123 CeedCall(CeedBasisDestroy(&fdm_basis)); 3124 CeedCall(CeedQFunctionDestroy(&qf)); 3125 CeedCall(CeedQFunctionDestroy(&qf_fdm)); 3126 return CEED_ERROR_SUCCESS; 3127 } 3128 3129 /// @} 3130