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 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 CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec)); 1064 CeedCall(CeedVectorDestroy(&vec)); 1065 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1066 CeedCall(CeedBasisDestroy(&basis)); 1067 } 1068 // -- Clone output fields 1069 for (CeedInt i = 0; i < num_output_fields; i++) { 1070 const char *field_name; 1071 CeedVector vec; 1072 CeedElemRestriction rstr = NULL; 1073 CeedBasis basis = NULL; 1074 1075 CeedCall(CeedOperatorFieldGetName(output_fields[i], &field_name)); 1076 CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); 1077 if (vec == CEED_VECTOR_ACTIVE) { 1078 CeedCall(CeedElemRestrictionReferenceCopy(rstr_coarse, &rstr)); 1079 CeedCall(CeedBasisReferenceCopy(basis_coarse, &basis)); 1080 if (!rstr_fine) CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_fine)); 1081 } else { 1082 CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr)); 1083 CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis)); 1084 } 1085 CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec)); 1086 CeedCall(CeedVectorDestroy(&vec)); 1087 CeedCall(CeedElemRestrictionDestroy(&rstr)); 1088 CeedCall(CeedBasisDestroy(&basis)); 1089 } 1090 // -- Clone QFunctionAssemblyData 1091 { 1092 CeedQFunctionAssemblyData fine_data; 1093 1094 CeedCall(CeedOperatorGetQFunctionAssemblyData(op_fine, &fine_data)); 1095 CeedCall(CeedQFunctionAssemblyDataReferenceCopy(fine_data, &(*op_coarse)->qf_assembled)); 1096 } 1097 1098 // Multiplicity vector 1099 if (op_restrict || op_prolong) { 1100 CeedInt num_elem, num_comp, elem_size; 1101 CeedVector mult_l_vec, mult_e_vec; 1102 CeedRestrictionType rstr_type; 1103 CeedElemRestriction rstr_p_mult_full; 1104 1105 CeedCall(CeedElemRestrictionGetType(rstr_fine, &rstr_type)); 1106 CeedCheck(rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_UNSUPPORTED, 1107 "Element restrictions created with CeedElemRestrictionCreateCurlOriented are not supported"); 1108 CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector"); 1109 1110 // Create multiplicity multi-component l-vector 1111 CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_full)); 1112 CeedCall(CeedElemRestrictionCreateVector(rstr_fine, &mult_l_vec, &mult_e_vec)); 1113 CeedCall(CeedVectorSetValue(mult_e_vec, 0.0)); 1114 CeedCall(CeedElemRestrictionApply(rstr_p_mult_full, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE)); 1115 CeedCall(CeedVectorSetValue(mult_l_vec, 0.0)); 1116 CeedCall(CeedElemRestrictionApply(rstr_p_mult_full, CEED_TRANSPOSE, mult_e_vec, mult_l_vec, CEED_REQUEST_IMMEDIATE)); 1117 CeedCall(CeedVectorReciprocal(mult_l_vec)); 1118 1119 // Create multiplicity single component e-vector 1120 CeedCall(CeedElemRestrictionGetNumElements(rstr_p_mult_full, &num_elem)); 1121 CeedCall(CeedElemRestrictionGetNumComponents(rstr_p_mult_full, &num_comp)); 1122 CeedCall(CeedElemRestrictionGetElementSize(rstr_p_mult_full, &elem_size)); 1123 CeedCall(CeedElemRestrictionCreateStrided(ceed, num_elem, elem_size, 1, num_elem * elem_size, CEED_STRIDES_BACKEND, &rstr_p_mult_fine)); 1124 CeedCall(CeedElemRestrictionCreateVector(rstr_p_mult_fine, &mult_vec, NULL)); 1125 { 1126 CeedQFunction qf_to_scalar; 1127 CeedOperator op_to_scalar; 1128 1129 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Identity to scalar", &qf_to_scalar)); 1130 CeedCall(CeedQFunctionAddInput(qf_to_scalar, "input", num_comp, CEED_EVAL_NONE)); 1131 CeedCall(CeedQFunctionAddOutput(qf_to_scalar, "output", 1, CEED_EVAL_NONE)); 1132 1133 CeedCall(CeedOperatorCreate(ceed, qf_to_scalar, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_to_scalar)); 1134 CeedCall(CeedOperatorSetField(op_to_scalar, "input", rstr_p_mult_full, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1135 CeedCall(CeedOperatorSetField(op_to_scalar, "output", rstr_p_mult_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1136 1137 CeedCall(CeedOperatorApply(op_to_scalar, mult_l_vec, mult_vec, CEED_REQUEST_IMMEDIATE)); 1138 1139 // Clean-up 1140 CeedCall(CeedQFunctionDestroy(&qf_to_scalar)); 1141 CeedCall(CeedOperatorDestroy(&op_to_scalar)); 1142 } 1143 // Clean-up 1144 CeedCall(CeedVectorDestroy(&mult_e_vec)); 1145 CeedCall(CeedVectorDestroy(&mult_l_vec)); 1146 CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_full)); 1147 } 1148 1149 // Clone name 1150 bool has_name = op_fine->name; 1151 size_t name_len = op_fine->name ? strlen(op_fine->name) : 0; 1152 CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name)); 1153 1154 // Check that coarse to fine basis is provided if prolong/restrict operators are requested 1155 CeedCheck(basis_c_to_f || (!op_restrict && !op_prolong), ceed, CEED_ERROR_INCOMPATIBLE, 1156 "Prolongation or restriction operator creation requires coarse-to-fine basis"); 1157 1158 // Restriction/Prolongation Operators 1159 CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp)); 1160 1161 // Restriction 1162 if (op_restrict) { 1163 CeedInt *num_comp_r_data; 1164 CeedQFunctionContext ctx_r; 1165 CeedQFunction qf_restrict; 1166 1167 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale (scalar)", &qf_restrict)); 1168 CeedCall(CeedCalloc(1, &num_comp_r_data)); 1169 num_comp_r_data[0] = num_comp; 1170 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r)); 1171 CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data)); 1172 CeedCall(CeedQFunctionSetContext(qf_restrict, ctx_r)); 1173 CeedCall(CeedQFunctionContextDestroy(&ctx_r)); 1174 CeedCall(CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE)); 1175 CeedCall(CeedQFunctionAddInput(qf_restrict, "scale", 1, CEED_EVAL_NONE)); 1176 CeedCall(CeedQFunctionAddOutput(qf_restrict, "output", num_comp, CEED_EVAL_INTERP)); 1177 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_restrict, num_comp)); 1178 1179 CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict)); 1180 CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1181 CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec)); 1182 CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); 1183 1184 // Set name 1185 char *restriction_name; 1186 1187 CeedCall(CeedCalloc(17 + name_len, &restriction_name)); 1188 sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); 1189 CeedCall(CeedOperatorSetName(*op_restrict, restriction_name)); 1190 CeedCall(CeedFree(&restriction_name)); 1191 1192 // Check 1193 CeedCall(CeedOperatorCheckReady(*op_restrict)); 1194 1195 // Cleanup 1196 CeedCall(CeedQFunctionDestroy(&qf_restrict)); 1197 } 1198 1199 // Prolongation 1200 if (op_prolong) { 1201 CeedInt *num_comp_p_data; 1202 CeedQFunctionContext ctx_p; 1203 CeedQFunction qf_prolong; 1204 1205 CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale (scalar)", &qf_prolong)); 1206 CeedCall(CeedCalloc(1, &num_comp_p_data)); 1207 num_comp_p_data[0] = num_comp; 1208 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_p)); 1209 CeedCall(CeedQFunctionContextSetData(ctx_p, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_p_data), num_comp_p_data)); 1210 CeedCall(CeedQFunctionSetContext(qf_prolong, ctx_p)); 1211 CeedCall(CeedQFunctionContextDestroy(&ctx_p)); 1212 CeedCall(CeedQFunctionAddInput(qf_prolong, "input", num_comp, CEED_EVAL_INTERP)); 1213 CeedCall(CeedQFunctionAddInput(qf_prolong, "scale", 1, CEED_EVAL_NONE)); 1214 CeedCall(CeedQFunctionAddOutput(qf_prolong, "output", num_comp, CEED_EVAL_NONE)); 1215 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_prolong, num_comp)); 1216 1217 CeedCall(CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_prolong)); 1218 CeedCall(CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); 1219 CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec)); 1220 CeedCall(CeedOperatorSetField(*op_prolong, "output", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); 1221 1222 // Set name 1223 char *prolongation_name; 1224 1225 CeedCall(CeedCalloc(18 + name_len, &prolongation_name)); 1226 sprintf(prolongation_name, "prolongation%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); 1227 CeedCall(CeedOperatorSetName(*op_prolong, prolongation_name)); 1228 CeedCall(CeedFree(&prolongation_name)); 1229 1230 // Check 1231 CeedCall(CeedOperatorCheckReady(*op_prolong)); 1232 1233 // Cleanup 1234 CeedCall(CeedQFunctionDestroy(&qf_prolong)); 1235 } 1236 1237 // Check 1238 CeedCall(CeedOperatorCheckReady(*op_coarse)); 1239 1240 // Cleanup 1241 CeedCall(CeedDestroy(&ceed)); 1242 CeedCall(CeedVectorDestroy(&mult_vec)); 1243 CeedCall(CeedElemRestrictionDestroy(&rstr_fine)); 1244 CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_fine)); 1245 CeedCall(CeedBasisDestroy(&basis_c_to_f)); 1246 return CEED_ERROR_SUCCESS; 1247 } 1248 1249 /** 1250 @brief Build 1D mass matrix and Laplacian with perturbation 1251 1252 @param[in] interp_1d Interpolation matrix in one dimension 1253 @param[in] grad_1d Gradient matrix in one dimension 1254 @param[in] q_weight_1d Quadrature weights in one dimension 1255 @param[in] P_1d Number of basis nodes in one dimension 1256 @param[in] Q_1d Number of quadrature points in one dimension 1257 @param[in] dim Dimension of basis 1258 @param[out] mass Assembled mass matrix in one dimension 1259 @param[out] laplace Assembled perturbed Laplacian in one dimension 1260 1261 @return An error code: 0 - success, otherwise - failure 1262 1263 @ref Developer 1264 **/ 1265 CeedPragmaOptimizeOff 1266 static int CeedBuildMassLaplace(const CeedScalar *interp_1d, const CeedScalar *grad_1d, const CeedScalar *q_weight_1d, CeedInt P_1d, CeedInt Q_1d, 1267 CeedInt dim, CeedScalar *mass, CeedScalar *laplace) { 1268 for (CeedInt i = 0; i < P_1d; i++) { 1269 for (CeedInt j = 0; j < P_1d; j++) { 1270 CeedScalar sum = 0.0; 1271 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]; 1272 mass[i + j * P_1d] = sum; 1273 } 1274 } 1275 // -- Laplacian 1276 for (CeedInt i = 0; i < P_1d; i++) { 1277 for (CeedInt j = 0; j < P_1d; j++) { 1278 CeedScalar sum = 0.0; 1279 1280 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]; 1281 laplace[i + j * P_1d] = sum; 1282 } 1283 } 1284 CeedScalar perturbation = dim > 2 ? 1e-6 : 1e-4; 1285 for (CeedInt i = 0; i < P_1d; i++) laplace[i + P_1d * i] += perturbation; 1286 return CEED_ERROR_SUCCESS; 1287 } 1288 CeedPragmaOptimizeOn 1289 1290 /// @} 1291 1292 /// ---------------------------------------------------------------------------- 1293 /// CeedOperator Backend API 1294 /// ---------------------------------------------------------------------------- 1295 /// @addtogroup CeedOperatorBackend 1296 /// @{ 1297 1298 /** 1299 @brief Select correct basis matrix pointer based on @ref CeedEvalMode 1300 1301 @param[in] basis `CeedBasis` from which to get the basis matrix 1302 @param[in] eval_mode Current basis evaluation mode 1303 @param[in] identity Pointer to identity matrix 1304 @param[out] basis_ptr `CeedBasis` pointer to set 1305 1306 @ref Backend 1307 **/ 1308 int CeedOperatorGetBasisPointer(CeedBasis basis, CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar **basis_ptr) { 1309 switch (eval_mode) { 1310 case CEED_EVAL_NONE: 1311 *basis_ptr = identity; 1312 break; 1313 case CEED_EVAL_INTERP: 1314 CeedCall(CeedBasisGetInterp(basis, basis_ptr)); 1315 break; 1316 case CEED_EVAL_GRAD: 1317 CeedCall(CeedBasisGetGrad(basis, basis_ptr)); 1318 break; 1319 case CEED_EVAL_DIV: 1320 CeedCall(CeedBasisGetDiv(basis, basis_ptr)); 1321 break; 1322 case CEED_EVAL_CURL: 1323 CeedCall(CeedBasisGetCurl(basis, basis_ptr)); 1324 break; 1325 case CEED_EVAL_WEIGHT: 1326 break; // Caught by QF Assembly 1327 } 1328 assert(*basis_ptr != NULL); 1329 return CEED_ERROR_SUCCESS; 1330 } 1331 1332 /** 1333 @brief Create point block restriction for active `CeedOperatorField` 1334 1335 @param[in] rstr Original `CeedElemRestriction` for active field 1336 @param[out] point_block_rstr Address of the variable where the newly created `CeedElemRestriction` will be stored 1337 1338 @return An error code: 0 - success, otherwise - failure 1339 1340 @ref Backend 1341 **/ 1342 int CeedOperatorCreateActivePointBlockRestriction(CeedElemRestriction rstr, CeedElemRestriction *point_block_rstr) { 1343 Ceed ceed; 1344 CeedInt num_elem, num_comp, shift, elem_size, comp_stride, *point_block_offsets; 1345 CeedSize l_size; 1346 const CeedInt *offsets; 1347 1348 CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); 1349 CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); 1350 1351 // Expand offsets 1352 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 1353 CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 1354 CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 1355 CeedCall(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 1356 CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 1357 shift = num_comp; 1358 if (comp_stride != 1) shift *= num_comp; 1359 CeedCall(CeedCalloc(num_elem * elem_size, &point_block_offsets)); 1360 for (CeedInt i = 0; i < num_elem * elem_size; i++) { 1361 point_block_offsets[i] = offsets[i] * shift; 1362 } 1363 1364 // Create new restriction 1365 CeedCall(CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp * num_comp, 1, l_size * num_comp, CEED_MEM_HOST, CEED_OWN_POINTER, 1366 point_block_offsets, point_block_rstr)); 1367 1368 // Cleanup 1369 CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); 1370 CeedCall(CeedDestroy(&ceed)); 1371 return CEED_ERROR_SUCCESS; 1372 } 1373 1374 /** 1375 @brief Get `CeedQFunctionAssemblyData` 1376 1377 @param[in] op `CeedOperator` to assemble 1378 @param[out] data `CeedQFunctionAssemblyData` 1379 1380 @return An error code: 0 - success, otherwise - failure 1381 1382 @ref Backend 1383 **/ 1384 int CeedOperatorGetQFunctionAssemblyData(CeedOperator op, CeedQFunctionAssemblyData *data) { 1385 if (!op->qf_assembled) { 1386 CeedQFunctionAssemblyData data; 1387 1388 CeedCall(CeedQFunctionAssemblyDataCreate(CeedOperatorReturnCeed(op), &data)); 1389 op->qf_assembled = data; 1390 } 1391 *data = op->qf_assembled; 1392 return CEED_ERROR_SUCCESS; 1393 } 1394 1395 /** 1396 @brief Create object holding `CeedQFunction` assembly data for `CeedOperator` 1397 1398 @param[in] ceed `Ceed` object used to create the `CeedQFunctionAssemblyData` 1399 @param[out] data Address of the variable where the newly created `CeedQFunctionAssemblyData` will be stored 1400 1401 @return An error code: 0 - success, otherwise - failure 1402 1403 @ref Backend 1404 **/ 1405 int CeedQFunctionAssemblyDataCreate(Ceed ceed, CeedQFunctionAssemblyData *data) { 1406 CeedCall(CeedCalloc(1, data)); 1407 (*data)->ref_count = 1; 1408 CeedCall(CeedReferenceCopy(ceed, &(*data)->ceed)); 1409 return CEED_ERROR_SUCCESS; 1410 } 1411 1412 /** 1413 @brief Increment the reference counter for a `CeedQFunctionAssemblyData` 1414 1415 @param[in,out] data `CeedQFunctionAssemblyData` to increment the reference counter 1416 1417 @return An error code: 0 - success, otherwise - failure 1418 1419 @ref Backend 1420 **/ 1421 int CeedQFunctionAssemblyDataReference(CeedQFunctionAssemblyData data) { 1422 data->ref_count++; 1423 return CEED_ERROR_SUCCESS; 1424 } 1425 1426 /** 1427 @brief Set re-use of `CeedQFunctionAssemblyData` 1428 1429 @param[in,out] data `CeedQFunctionAssemblyData` to mark for reuse 1430 @param[in] reuse_data Boolean flag indicating data re-use 1431 1432 @return An error code: 0 - success, otherwise - failure 1433 1434 @ref Backend 1435 **/ 1436 int CeedQFunctionAssemblyDataSetReuse(CeedQFunctionAssemblyData data, bool reuse_data) { 1437 data->reuse_data = reuse_data; 1438 data->needs_data_update = true; 1439 return CEED_ERROR_SUCCESS; 1440 } 1441 1442 /** 1443 @brief Mark `CeedQFunctionAssemblyData` as stale 1444 1445 @param[in,out] data `CeedQFunctionAssemblyData` to mark as stale 1446 @param[in] needs_data_update Boolean flag indicating if update is needed or completed 1447 1448 @return An error code: 0 - success, otherwise - failure 1449 1450 @ref Backend 1451 **/ 1452 int CeedQFunctionAssemblyDataSetUpdateNeeded(CeedQFunctionAssemblyData data, bool needs_data_update) { 1453 data->needs_data_update = needs_data_update; 1454 return CEED_ERROR_SUCCESS; 1455 } 1456 1457 /** 1458 @brief Determine if `CeedQFunctionAssemblyData` needs update 1459 1460 @param[in] data `CeedQFunctionAssemblyData` to mark as stale 1461 @param[out] is_update_needed Boolean flag indicating if re-assembly is required 1462 1463 @return An error code: 0 - success, otherwise - failure 1464 1465 @ref Backend 1466 **/ 1467 int CeedQFunctionAssemblyDataIsUpdateNeeded(CeedQFunctionAssemblyData data, bool *is_update_needed) { 1468 *is_update_needed = !data->reuse_data || data->needs_data_update; 1469 return CEED_ERROR_SUCCESS; 1470 } 1471 1472 /** 1473 @brief Copy the pointer to a `CeedQFunctionAssemblyData`. 1474 1475 Both pointers should be destroyed with @ref CeedQFunctionAssemblyDataDestroy(). 1476 1477 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`. 1478 This `CeedQFunctionAssemblyData` will be destroyed if ` *data_copy` is the only reference to this `CeedQFunctionAssemblyData`. 1479 1480 @param[in] data `CeedQFunctionAssemblyData` to copy reference to 1481 @param[in,out] data_copy Variable to store copied reference 1482 1483 @return An error code: 0 - success, otherwise - failure 1484 1485 @ref Backend 1486 **/ 1487 int CeedQFunctionAssemblyDataReferenceCopy(CeedQFunctionAssemblyData data, CeedQFunctionAssemblyData *data_copy) { 1488 CeedCall(CeedQFunctionAssemblyDataReference(data)); 1489 CeedCall(CeedQFunctionAssemblyDataDestroy(data_copy)); 1490 *data_copy = data; 1491 return CEED_ERROR_SUCCESS; 1492 } 1493 1494 /** 1495 @brief Get setup status for internal objects for `CeedQFunctionAssemblyData` 1496 1497 @param[in] data `CeedQFunctionAssemblyData` to retrieve status 1498 @param[out] is_setup Boolean flag for setup status 1499 1500 @return An error code: 0 - success, otherwise - failure 1501 1502 @ref Backend 1503 **/ 1504 int CeedQFunctionAssemblyDataIsSetup(CeedQFunctionAssemblyData data, bool *is_setup) { 1505 *is_setup = data->is_setup; 1506 return CEED_ERROR_SUCCESS; 1507 } 1508 1509 /** 1510 @brief Set internal objects for `CeedQFunctionAssemblyData` 1511 1512 @param[in,out] data `CeedQFunctionAssemblyData` to set objects 1513 @param[in] vec `CeedVector` to store assembled `CeedQFunction` at quadrature points 1514 @param[in] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 1515 1516 @return An error code: 0 - success, otherwise - failure 1517 1518 @ref Backend 1519 **/ 1520 int CeedQFunctionAssemblyDataSetObjects(CeedQFunctionAssemblyData data, CeedVector vec, CeedElemRestriction rstr) { 1521 CeedCall(CeedVectorReferenceCopy(vec, &data->vec)); 1522 CeedCall(CeedElemRestrictionReferenceCopy(rstr, &data->rstr)); 1523 1524 data->is_setup = true; 1525 return CEED_ERROR_SUCCESS; 1526 } 1527 1528 /** 1529 @brief Get internal objects for `CeedQFunctionAssemblyData` 1530 1531 @param[in,out] data `CeedQFunctionAssemblyData` to set objects 1532 @param[out] vec `CeedVector` to store assembled `CeedQFunction` at quadrature points 1533 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 1534 1535 @return An error code: 0 - success, otherwise - failure 1536 1537 @ref Backend 1538 **/ 1539 int CeedQFunctionAssemblyDataGetObjects(CeedQFunctionAssemblyData data, CeedVector *vec, CeedElemRestriction *rstr) { 1540 CeedCheck(data->is_setup, data->ceed, CEED_ERROR_INCOMPLETE, "Internal objects not set; must call CeedQFunctionAssemblyDataSetObjects first."); 1541 1542 CeedCall(CeedVectorReferenceCopy(data->vec, vec)); 1543 CeedCall(CeedElemRestrictionReferenceCopy(data->rstr, rstr)); 1544 return CEED_ERROR_SUCCESS; 1545 } 1546 1547 /** 1548 @brief Destroy `CeedQFunctionAssemblyData` 1549 1550 @param[in,out] data `CeedQFunctionAssemblyData` to destroy 1551 1552 @return An error code: 0 - success, otherwise - failure 1553 1554 @ref Backend 1555 **/ 1556 int CeedQFunctionAssemblyDataDestroy(CeedQFunctionAssemblyData *data) { 1557 if (!*data || --(*data)->ref_count > 0) { 1558 *data = NULL; 1559 return CEED_ERROR_SUCCESS; 1560 } 1561 CeedCall(CeedDestroy(&(*data)->ceed)); 1562 CeedCall(CeedVectorDestroy(&(*data)->vec)); 1563 CeedCall(CeedElemRestrictionDestroy(&(*data)->rstr)); 1564 1565 CeedCall(CeedFree(data)); 1566 return CEED_ERROR_SUCCESS; 1567 } 1568 1569 /** 1570 @brief Get `CeedOperatorAssemblyData` 1571 1572 @param[in] op `CeedOperator` to assemble 1573 @param[out] data `CeedOperatorAssemblyData` 1574 1575 @return An error code: 0 - success, otherwise - failure 1576 1577 @ref Backend 1578 **/ 1579 int CeedOperatorGetOperatorAssemblyData(CeedOperator op, CeedOperatorAssemblyData *data) { 1580 if (!op->op_assembled) { 1581 CeedOperatorAssemblyData data; 1582 1583 CeedCall(CeedOperatorAssemblyDataCreate(CeedOperatorReturnCeed(op), op, &data)); 1584 op->op_assembled = data; 1585 } 1586 *data = op->op_assembled; 1587 return CEED_ERROR_SUCCESS; 1588 } 1589 1590 /** 1591 @brief Create object holding `CeedOperator` assembly data. 1592 1593 The `CeedOperatorAssemblyData` holds an array with references to every active `CeedBasis` used in the `CeedOperator`. 1594 An array with references to the corresponding active `CeedElemRestriction` is also stored. 1595 For each active `CeedBasis, the `CeedOperatorAssemblyData` holds an array of all input and output @ref CeedEvalMode for this `CeedBasis`. 1596 The `CeedOperatorAssemblyData` holds an array of offsets for indexing into the assembled `CeedQFunction` arrays to the row representing each @ref CeedEvalMode. 1597 The number of input columns across all active bases for the assembled `CeedQFunction` is also stored. 1598 Lastly, the `CeedOperatorAssembly` data holds assembled matrices representing the full action of the `CeedBasis` for all @ref CeedEvalMode. 1599 1600 @param[in] ceed `Ceed` object used to create the `CeedOperatorAssemblyData` 1601 @param[in] op `CeedOperator` to be assembled 1602 @param[out] data Address of the variable where the newly created `CeedOperatorAssemblyData` will be stored 1603 1604 @return An error code: 0 - success, otherwise - failure 1605 1606 @ref Backend 1607 **/ 1608 int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) { 1609 CeedInt num_active_bases_in = 0, num_active_bases_out = 0, offset = 0; 1610 CeedInt num_input_fields, *num_eval_modes_in = NULL, num_output_fields, *num_eval_modes_out = NULL; 1611 CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL; 1612 CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL; 1613 CeedQFunctionField *qf_fields; 1614 CeedQFunction qf; 1615 CeedOperatorField *op_fields; 1616 bool is_composite; 1617 1618 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 1619 CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Can only create CeedOperator assembly data for non-composite operators."); 1620 1621 // Allocate 1622 CeedCall(CeedCalloc(1, data)); 1623 CeedCall(CeedReferenceCopy(ceed, &(*data)->ceed)); 1624 1625 // Build OperatorAssembly data 1626 CeedCall(CeedOperatorGetQFunction(op, &qf)); 1627 1628 // Determine active input basis 1629 CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL)); 1630 CeedCall(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); 1631 for (CeedInt i = 0; i < num_input_fields; i++) { 1632 CeedVector vec; 1633 1634 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1635 if (vec == CEED_VECTOR_ACTIVE) { 1636 CeedInt index = -1, num_comp, q_comp; 1637 CeedEvalMode eval_mode; 1638 CeedBasis basis_in = NULL; 1639 1640 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); 1641 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1642 CeedCall(CeedBasisGetNumComponents(basis_in, &num_comp)); 1643 CeedCall(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); 1644 for (CeedInt i = 0; i < num_active_bases_in; i++) { 1645 if ((*data)->active_bases_in[i] == basis_in) index = i; 1646 } 1647 if (index == -1) { 1648 CeedElemRestriction elem_rstr_in; 1649 1650 index = num_active_bases_in; 1651 CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_bases_in)); 1652 (*data)->active_bases_in[num_active_bases_in] = NULL; 1653 CeedCall(CeedBasisReferenceCopy(basis_in, &(*data)->active_bases_in[num_active_bases_in])); 1654 CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_elem_rstrs_in)); 1655 (*data)->active_elem_rstrs_in[num_active_bases_in] = NULL; 1656 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_in)); 1657 CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_in, &(*data)->active_elem_rstrs_in[num_active_bases_in])); 1658 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_in)); 1659 CeedCall(CeedRealloc(num_active_bases_in + 1, &num_eval_modes_in)); 1660 num_eval_modes_in[index] = 0; 1661 CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_modes_in)); 1662 eval_modes_in[index] = NULL; 1663 CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_mode_offsets_in)); 1664 eval_mode_offsets_in[index] = NULL; 1665 CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->assembled_bases_in)); 1666 (*data)->assembled_bases_in[index] = NULL; 1667 num_active_bases_in++; 1668 } 1669 if (eval_mode != CEED_EVAL_WEIGHT) { 1670 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1671 CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_modes_in[index])); 1672 CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_mode_offsets_in[index])); 1673 for (CeedInt d = 0; d < q_comp; d++) { 1674 eval_modes_in[index][num_eval_modes_in[index] + d] = eval_mode; 1675 eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset; 1676 offset += num_comp; 1677 } 1678 num_eval_modes_in[index] += q_comp; 1679 } 1680 CeedCall(CeedBasisDestroy(&basis_in)); 1681 } 1682 CeedCall(CeedVectorDestroy(&vec)); 1683 } 1684 1685 // Determine active output basis 1686 CeedCall(CeedQFunctionGetFields(qf, NULL, NULL, &num_output_fields, &qf_fields)); 1687 CeedCall(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); 1688 offset = 0; 1689 for (CeedInt i = 0; i < num_output_fields; i++) { 1690 CeedVector vec; 1691 1692 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 1693 if (vec == CEED_VECTOR_ACTIVE) { 1694 CeedInt index = -1, num_comp, q_comp; 1695 CeedEvalMode eval_mode; 1696 CeedBasis basis_out = NULL; 1697 1698 CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); 1699 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 1700 CeedCall(CeedBasisGetNumComponents(basis_out, &num_comp)); 1701 CeedCall(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); 1702 for (CeedInt i = 0; i < num_active_bases_out; i++) { 1703 if ((*data)->active_bases_out[i] == basis_out) index = i; 1704 } 1705 if (index == -1) { 1706 CeedElemRestriction elem_rstr_out; 1707 1708 index = num_active_bases_out; 1709 CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_bases_out)); 1710 (*data)->active_bases_out[num_active_bases_out] = NULL; 1711 CeedCall(CeedBasisReferenceCopy(basis_out, &(*data)->active_bases_out[num_active_bases_out])); 1712 CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_elem_rstrs_out)); 1713 (*data)->active_elem_rstrs_out[num_active_bases_out] = NULL; 1714 CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_out)); 1715 CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_out, &(*data)->active_elem_rstrs_out[num_active_bases_out])); 1716 CeedCall(CeedElemRestrictionDestroy(&elem_rstr_out)); 1717 CeedCall(CeedRealloc(num_active_bases_out + 1, &num_eval_modes_out)); 1718 num_eval_modes_out[index] = 0; 1719 CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_modes_out)); 1720 eval_modes_out[index] = NULL; 1721 CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_mode_offsets_out)); 1722 eval_mode_offsets_out[index] = NULL; 1723 CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->assembled_bases_out)); 1724 (*data)->assembled_bases_out[index] = NULL; 1725 num_active_bases_out++; 1726 } 1727 if (eval_mode != CEED_EVAL_WEIGHT) { 1728 // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly 1729 CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_modes_out[index])); 1730 CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_mode_offsets_out[index])); 1731 for (CeedInt d = 0; d < q_comp; d++) { 1732 eval_modes_out[index][num_eval_modes_out[index] + d] = eval_mode; 1733 eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset; 1734 offset += num_comp; 1735 } 1736 num_eval_modes_out[index] += q_comp; 1737 } 1738 CeedCall(CeedBasisDestroy(&basis_out)); 1739 } 1740 CeedCall(CeedVectorDestroy(&vec)); 1741 } 1742 CeedCall(CeedQFunctionDestroy(&qf)); 1743 (*data)->num_active_bases_in = num_active_bases_in; 1744 (*data)->num_eval_modes_in = num_eval_modes_in; 1745 (*data)->eval_modes_in = eval_modes_in; 1746 (*data)->eval_mode_offsets_in = eval_mode_offsets_in; 1747 (*data)->num_active_bases_out = num_active_bases_out; 1748 (*data)->num_eval_modes_out = num_eval_modes_out; 1749 (*data)->eval_modes_out = eval_modes_out; 1750 (*data)->eval_mode_offsets_out = eval_mode_offsets_out; 1751 (*data)->num_output_components = offset; 1752 return CEED_ERROR_SUCCESS; 1753 } 1754 1755 /** 1756 @brief Get `CeedOperator` @ref CeedEvalMode for assembly. 1757 1758 Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. 1759 1760 @param[in] data `CeedOperatorAssemblyData` 1761 @param[out] num_active_bases_in Total number of active bases for input 1762 @param[out] num_eval_modes_in Pointer to hold array of numbers of input @ref CeedEvalMode, or `NULL`. 1763 `eval_modes_in[0]` holds an array of eval modes for the first active `CeedBasis`. 1764 @param[out] eval_modes_in Pointer to hold arrays of input @ref CeedEvalMode, or `NULL` 1765 @param[out] eval_mode_offsets_in Pointer to hold arrays of input offsets at each quadrature point 1766 @param[out] num_active_bases_out Total number of active bases for output 1767 @param[out] num_eval_modes_out Pointer to hold array of numbers of output @ref CeedEvalMode, or `NULL` 1768 @param[out] eval_modes_out Pointer to hold arrays of output @ref CeedEvalMode, or `NULL` 1769 @param[out] eval_mode_offsets_out Pointer to hold arrays of output offsets at each quadrature point 1770 @param[out] num_output_components The number of columns in the assembled `CeedQFunction` matrix for each quadrature point, including contributions of all active bases 1771 1772 @return An error code: 0 - success, otherwise - failure 1773 1774 @ref Backend 1775 **/ 1776 int CeedOperatorAssemblyDataGetEvalModes(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedInt **num_eval_modes_in, 1777 const CeedEvalMode ***eval_modes_in, CeedSize ***eval_mode_offsets_in, CeedInt *num_active_bases_out, 1778 CeedInt **num_eval_modes_out, const CeedEvalMode ***eval_modes_out, CeedSize ***eval_mode_offsets_out, 1779 CeedSize *num_output_components) { 1780 if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in; 1781 if (num_eval_modes_in) *num_eval_modes_in = data->num_eval_modes_in; 1782 if (eval_modes_in) *eval_modes_in = (const CeedEvalMode **)data->eval_modes_in; 1783 if (eval_mode_offsets_in) *eval_mode_offsets_in = data->eval_mode_offsets_in; 1784 if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out; 1785 if (num_eval_modes_out) *num_eval_modes_out = data->num_eval_modes_out; 1786 if (eval_modes_out) *eval_modes_out = (const CeedEvalMode **)data->eval_modes_out; 1787 if (eval_mode_offsets_out) *eval_mode_offsets_out = data->eval_mode_offsets_out; 1788 if (num_output_components) *num_output_components = data->num_output_components; 1789 return CEED_ERROR_SUCCESS; 1790 } 1791 1792 /** 1793 @brief Get `CeedOperator` `CeedBasis` data for assembly. 1794 1795 Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. 1796 1797 @param[in] data `CeedOperatorAssemblyData` 1798 @param[out] num_active_bases_in Number of active input bases, or `NULL` 1799 @param[out] active_bases_in Pointer to hold active input `CeedBasis`, or `NULL` 1800 @param[out] assembled_bases_in Pointer to hold assembled active input `B` , or `NULL` 1801 @param[out] num_active_bases_out Number of active output bases, or `NULL` 1802 @param[out] active_bases_out Pointer to hold active output `CeedBasis`, or `NULL` 1803 @param[out] assembled_bases_out Pointer to hold assembled active output `B` , or `NULL` 1804 1805 @return An error code: 0 - success, otherwise - failure 1806 1807 @ref Backend 1808 **/ 1809 int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedBasis **active_bases_in, 1810 const CeedScalar ***assembled_bases_in, CeedInt *num_active_bases_out, CeedBasis **active_bases_out, 1811 const CeedScalar ***assembled_bases_out) { 1812 // Assemble B_in, B_out if needed 1813 if (assembled_bases_in && !data->assembled_bases_in[0]) { 1814 CeedInt num_qpts; 1815 1816 if (data->active_bases_in[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[0], &num_qpts)); 1817 else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_in[0], &num_qpts)); 1818 for (CeedInt b = 0; b < data->num_active_bases_in; b++) { 1819 bool has_eval_none = false; 1820 CeedInt num_nodes; 1821 CeedScalar *B_in = NULL, *identity = NULL; 1822 1823 CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[b], &num_nodes)); 1824 CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_in[b], &B_in)); 1825 1826 for (CeedInt i = 0; i < data->num_eval_modes_in[b]; i++) { 1827 has_eval_none = has_eval_none || (data->eval_modes_in[b][i] == CEED_EVAL_NONE); 1828 } 1829 if (has_eval_none) { 1830 CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); 1831 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { 1832 identity[i * num_nodes + i] = 1.0; 1833 } 1834 } 1835 1836 for (CeedInt q = 0; q < num_qpts; q++) { 1837 for (CeedInt n = 0; n < num_nodes; n++) { 1838 CeedInt d_in = 0, q_comp_in; 1839 CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; 1840 1841 for (CeedInt e_in = 0; e_in < data->num_eval_modes_in[b]; e_in++) { 1842 const CeedInt qq = data->num_eval_modes_in[b] * q; 1843 const CeedScalar *B = NULL; 1844 1845 CeedCall(CeedOperatorGetBasisPointer(data->active_bases_in[b], data->eval_modes_in[b][e_in], identity, &B)); 1846 CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_in[b], data->eval_modes_in[b][e_in], &q_comp_in)); 1847 if (q_comp_in > 1) { 1848 if (e_in == 0 || data->eval_modes_in[b][e_in] != eval_mode_in_prev) d_in = 0; 1849 else B = &B[(++d_in) * num_qpts * num_nodes]; 1850 } 1851 eval_mode_in_prev = data->eval_modes_in[b][e_in]; 1852 B_in[(qq + e_in) * num_nodes + n] = B[q * num_nodes + n]; 1853 } 1854 } 1855 } 1856 if (identity) CeedCall(CeedFree(&identity)); 1857 data->assembled_bases_in[b] = B_in; 1858 } 1859 } 1860 1861 if (assembled_bases_out && !data->assembled_bases_out[0]) { 1862 CeedInt num_qpts; 1863 1864 if (data->active_bases_out[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[0], &num_qpts)); 1865 else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_out[0], &num_qpts)); 1866 for (CeedInt b = 0; b < data->num_active_bases_out; b++) { 1867 bool has_eval_none = false; 1868 CeedInt num_nodes; 1869 CeedScalar *B_out = NULL, *identity = NULL; 1870 1871 CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[b], &num_nodes)); 1872 CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_out[b], &B_out)); 1873 1874 for (CeedInt i = 0; i < data->num_eval_modes_out[b]; i++) { 1875 has_eval_none = has_eval_none || (data->eval_modes_out[b][i] == CEED_EVAL_NONE); 1876 } 1877 if (has_eval_none) { 1878 CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); 1879 for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { 1880 identity[i * num_nodes + i] = 1.0; 1881 } 1882 } 1883 1884 for (CeedInt q = 0; q < num_qpts; q++) { 1885 for (CeedInt n = 0; n < num_nodes; n++) { 1886 CeedInt d_out = 0, q_comp_out; 1887 CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; 1888 1889 for (CeedInt e_out = 0; e_out < data->num_eval_modes_out[b]; e_out++) { 1890 const CeedInt qq = data->num_eval_modes_out[b] * q; 1891 const CeedScalar *B = NULL; 1892 1893 CeedCall(CeedOperatorGetBasisPointer(data->active_bases_out[b], data->eval_modes_out[b][e_out], identity, &B)); 1894 CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_out[b], data->eval_modes_out[b][e_out], &q_comp_out)); 1895 if (q_comp_out > 1) { 1896 if (e_out == 0 || data->eval_modes_out[b][e_out] != eval_mode_out_prev) d_out = 0; 1897 else B = &B[(++d_out) * num_qpts * num_nodes]; 1898 } 1899 eval_mode_out_prev = data->eval_modes_out[b][e_out]; 1900 B_out[(qq + e_out) * num_nodes + n] = B[q * num_nodes + n]; 1901 } 1902 } 1903 } 1904 if (identity) CeedCall(CeedFree(&identity)); 1905 data->assembled_bases_out[b] = B_out; 1906 } 1907 } 1908 1909 // Pass out assembled data 1910 if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in; 1911 if (active_bases_in) *active_bases_in = data->active_bases_in; 1912 if (assembled_bases_in) *assembled_bases_in = (const CeedScalar **)data->assembled_bases_in; 1913 if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out; 1914 if (active_bases_out) *active_bases_out = data->active_bases_out; 1915 if (assembled_bases_out) *assembled_bases_out = (const CeedScalar **)data->assembled_bases_out; 1916 return CEED_ERROR_SUCCESS; 1917 } 1918 1919 /** 1920 @brief Get `CeedOperator` `CeedBasis` data for assembly. 1921 1922 Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. 1923 1924 @param[in] data `CeedOperatorAssemblyData` 1925 @param[out] num_active_elem_rstrs_in Number of active input element restrictions, or `NULL` 1926 @param[out] active_elem_rstrs_in Pointer to hold active input `CeedElemRestriction`, or `NULL` 1927 @param[out] num_active_elem_rstrs_out Number of active output element restrictions, or `NULL` 1928 @param[out] active_elem_rstrs_out Pointer to hold active output `CeedElemRestriction`, or `NULL` 1929 1930 @return An error code: 0 - success, otherwise - failure 1931 1932 @ref Backend 1933 **/ 1934 int CeedOperatorAssemblyDataGetElemRestrictions(CeedOperatorAssemblyData data, CeedInt *num_active_elem_rstrs_in, 1935 CeedElemRestriction **active_elem_rstrs_in, CeedInt *num_active_elem_rstrs_out, 1936 CeedElemRestriction **active_elem_rstrs_out) { 1937 if (num_active_elem_rstrs_in) *num_active_elem_rstrs_in = data->num_active_bases_in; 1938 if (active_elem_rstrs_in) *active_elem_rstrs_in = data->active_elem_rstrs_in; 1939 if (num_active_elem_rstrs_out) *num_active_elem_rstrs_out = data->num_active_bases_out; 1940 if (active_elem_rstrs_out) *active_elem_rstrs_out = data->active_elem_rstrs_out; 1941 return CEED_ERROR_SUCCESS; 1942 } 1943 1944 /** 1945 @brief Destroy `CeedOperatorAssemblyData` 1946 1947 @param[in,out] data `CeedOperatorAssemblyData` to destroy 1948 1949 @return An error code: 0 - success, otherwise - failure 1950 1951 @ref Backend 1952 **/ 1953 int CeedOperatorAssemblyDataDestroy(CeedOperatorAssemblyData *data) { 1954 if (!*data) { 1955 *data = NULL; 1956 return CEED_ERROR_SUCCESS; 1957 } 1958 CeedCall(CeedDestroy(&(*data)->ceed)); 1959 for (CeedInt b = 0; b < (*data)->num_active_bases_in; b++) { 1960 CeedCall(CeedBasisDestroy(&(*data)->active_bases_in[b])); 1961 CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_in[b])); 1962 CeedCall(CeedFree(&(*data)->eval_modes_in[b])); 1963 CeedCall(CeedFree(&(*data)->eval_mode_offsets_in[b])); 1964 CeedCall(CeedFree(&(*data)->assembled_bases_in[b])); 1965 } 1966 for (CeedInt b = 0; b < (*data)->num_active_bases_out; b++) { 1967 CeedCall(CeedBasisDestroy(&(*data)->active_bases_out[b])); 1968 CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_out[b])); 1969 CeedCall(CeedFree(&(*data)->eval_modes_out[b])); 1970 CeedCall(CeedFree(&(*data)->eval_mode_offsets_out[b])); 1971 CeedCall(CeedFree(&(*data)->assembled_bases_out[b])); 1972 } 1973 CeedCall(CeedFree(&(*data)->active_bases_in)); 1974 CeedCall(CeedFree(&(*data)->active_bases_out)); 1975 CeedCall(CeedFree(&(*data)->active_elem_rstrs_in)); 1976 CeedCall(CeedFree(&(*data)->active_elem_rstrs_out)); 1977 CeedCall(CeedFree(&(*data)->num_eval_modes_in)); 1978 CeedCall(CeedFree(&(*data)->num_eval_modes_out)); 1979 CeedCall(CeedFree(&(*data)->eval_modes_in)); 1980 CeedCall(CeedFree(&(*data)->eval_modes_out)); 1981 CeedCall(CeedFree(&(*data)->eval_mode_offsets_in)); 1982 CeedCall(CeedFree(&(*data)->eval_mode_offsets_out)); 1983 CeedCall(CeedFree(&(*data)->assembled_bases_in)); 1984 CeedCall(CeedFree(&(*data)->assembled_bases_out)); 1985 1986 CeedCall(CeedFree(data)); 1987 return CEED_ERROR_SUCCESS; 1988 } 1989 1990 /** 1991 @brief Retrieve fallback `CeedOperator` with a reference `Ceed` for advanced `CeedOperator` functionality 1992 1993 @param[in] op `CeedOperator` to retrieve fallback for 1994 @param[out] op_fallback Fallback `CeedOperator` 1995 1996 @return An error code: 0 - success, otherwise - failure 1997 1998 @ref Backend 1999 **/ 2000 int CeedOperatorGetFallback(CeedOperator op, CeedOperator *op_fallback) { 2001 // Create if needed 2002 if (!op->op_fallback) CeedCall(CeedOperatorCreateFallback(op)); 2003 if (op->op_fallback) { 2004 bool is_debug; 2005 Ceed ceed; 2006 2007 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2008 CeedCall(CeedIsDebug(ceed, &is_debug)); 2009 if (is_debug) { 2010 Ceed ceed_fallback; 2011 const char *resource, *resource_fallback, *op_name; 2012 2013 CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback)); 2014 CeedCall(CeedGetResource(ceed, &resource)); 2015 CeedCall(CeedGetResource(ceed_fallback, &resource_fallback)); 2016 CeedCall(CeedOperatorGetName(op, &op_name)); 2017 2018 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "---------- CeedOperator Fallback ----------\n"); 2019 CeedDebug(ceed, "Falling back from Operator with backend %s at address %p to Operator with backend %s at address %p for CeedOperator \"%s\"\n", 2020 resource, op, resource_fallback, op->op_fallback, op_name); 2021 CeedCall(CeedDestroy(&ceed_fallback)); 2022 } 2023 CeedCall(CeedDestroy(&ceed)); 2024 } 2025 *op_fallback = op->op_fallback; 2026 return CEED_ERROR_SUCCESS; 2027 } 2028 2029 /** 2030 @brief Get the parent `CeedOperator` for a fallback `CeedOperator` 2031 2032 @param[in] op `CeedOperator` context 2033 @param[out] parent Variable to store parent `CeedOperator` context 2034 2035 @return An error code: 0 - success, otherwise - failure 2036 2037 @ref Backend 2038 **/ 2039 int CeedOperatorGetFallbackParent(CeedOperator op, CeedOperator *parent) { 2040 *parent = op->op_fallback_parent ? op->op_fallback_parent : NULL; 2041 return CEED_ERROR_SUCCESS; 2042 } 2043 2044 /** 2045 @brief Get the `Ceed` context of the parent `CeedOperator` for a fallback `CeedOperator` 2046 2047 @param[in] op `CeedOperator` context 2048 @param[out] parent Variable to store parent `Ceed` context 2049 2050 @return An error code: 0 - success, otherwise - failure 2051 2052 @ref Backend 2053 **/ 2054 int CeedOperatorGetFallbackParentCeed(CeedOperator op, Ceed *parent) { 2055 *parent = NULL; 2056 if (op->op_fallback_parent) CeedCall(CeedReferenceCopy(CeedOperatorReturnCeed(op->op_fallback_parent), parent)); 2057 else CeedCall(CeedReferenceCopy(CeedOperatorReturnCeed(op), parent)); 2058 return CEED_ERROR_SUCCESS; 2059 } 2060 2061 /// @} 2062 2063 /// ---------------------------------------------------------------------------- 2064 /// CeedOperator Public API 2065 /// ---------------------------------------------------------------------------- 2066 /// @addtogroup CeedOperatorUser 2067 /// @{ 2068 2069 /** 2070 @brief Assemble a linear `CeedQFunction` associated with a `CeedOperator`. 2071 2072 This returns a `CeedVector` containing a matrix at each quadrature point providing the action of the `CeedQFunction` associated with the `CeedOperator`. 2073 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. 2074 2075 Inputs and outputs are in the order provided by the user when adding `CeedOperator` fields. 2076 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]`. 2077 2078 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2079 2080 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2081 @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points 2082 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 2083 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2084 2085 @return An error code: 0 - success, otherwise - failure 2086 2087 @ref User 2088 **/ 2089 int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 2090 CeedCall(CeedOperatorCheckReady(op)); 2091 2092 if (op->LinearAssembleQFunction) { 2093 // Backend version 2094 CeedCall(op->LinearAssembleQFunction(op, assembled, rstr, request)); 2095 } else { 2096 // Operator fallback 2097 CeedOperator op_fallback; 2098 2099 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleQFunction\n"); 2100 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2101 if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunction(op_fallback, assembled, rstr, request)); 2102 else return CeedError(CeedOperatorReturnCeed(op), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunction"); 2103 } 2104 return CEED_ERROR_SUCCESS; 2105 } 2106 2107 /** 2108 @brief Assemble `CeedQFunction` and store result internally. 2109 2110 Return copied references of stored data to the caller. 2111 Caller is responsible for ownership and destruction of the copied references. 2112 See also @ref CeedOperatorLinearAssembleQFunction(). 2113 2114 Note: If the value of `assembled` or `rstr` passed to this function are non-`NULL` , then it is assumed that they hold valid pointers. 2115 These objects will be destroyed if `*assembled` or `*rstr` is the only reference to the object. 2116 2117 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2118 @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points 2119 @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` 2120 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2121 2122 @return An error code: 0 - success, otherwise - failure 2123 2124 @ref User 2125 **/ 2126 int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { 2127 return CeedOperatorLinearAssembleQFunctionBuildOrUpdate_Core(op, true, assembled, rstr, request); 2128 } 2129 2130 /** 2131 @brief Assemble the diagonal of a square linear `CeedOperator` 2132 2133 This overwrites a `CeedVector` with the diagonal of a linear `CeedOperator`. 2134 2135 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2136 2137 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2138 2139 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2140 @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal 2141 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2142 2143 @return An error code: 0 - success, otherwise - failure 2144 2145 @ref User 2146 **/ 2147 int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2148 bool is_composite; 2149 CeedSize input_size = 0, output_size = 0; 2150 2151 CeedCall(CeedOperatorCheckReady(op)); 2152 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2153 2154 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2155 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2156 2157 // Early exit for empty operator 2158 if (!is_composite) { 2159 CeedInt num_elem = 0; 2160 2161 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2162 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2163 } 2164 2165 if (op->LinearAssembleDiagonal) { 2166 // Backend version 2167 CeedCall(op->LinearAssembleDiagonal(op, assembled, request)); 2168 return CEED_ERROR_SUCCESS; 2169 } else if (op->LinearAssembleAddDiagonal) { 2170 // Backend version with zeroing first 2171 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2172 CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); 2173 return CEED_ERROR_SUCCESS; 2174 } else if (is_composite) { 2175 // Default to summing contributions of suboperators 2176 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2177 CeedCall(CeedOperatorLinearAssembleAddDiagonalComposite(op, request, false, assembled)); 2178 return CEED_ERROR_SUCCESS; 2179 } else { 2180 // Operator fallback 2181 CeedOperator op_fallback; 2182 2183 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleDiagonal\n"); 2184 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2185 if (op_fallback) { 2186 CeedCall(CeedOperatorLinearAssembleDiagonal(op_fallback, assembled, request)); 2187 return CEED_ERROR_SUCCESS; 2188 } 2189 } 2190 // Default interface implementation 2191 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2192 CeedCall(CeedOperatorLinearAssembleAddDiagonal(op, assembled, request)); 2193 return CEED_ERROR_SUCCESS; 2194 } 2195 2196 /** 2197 @brief Assemble the diagonal of a square linear `CeedOperator`. 2198 2199 This sums into a `CeedVector` the diagonal of a linear `CeedOperator`. 2200 2201 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2202 2203 Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. 2204 2205 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2206 @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal 2207 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2208 2209 @return An error code: 0 - success, otherwise - failure 2210 2211 @ref User 2212 **/ 2213 int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2214 bool is_composite; 2215 CeedSize input_size = 0, output_size = 0; 2216 2217 CeedCall(CeedOperatorCheckReady(op)); 2218 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2219 2220 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2221 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2222 2223 // Early exit for empty operator 2224 if (!is_composite) { 2225 CeedInt num_elem = 0; 2226 2227 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2228 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2229 } 2230 2231 if (op->LinearAssembleAddDiagonal) { 2232 // Backend version 2233 CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); 2234 return CEED_ERROR_SUCCESS; 2235 } else if (is_composite) { 2236 // Default to summing contributions of suboperators 2237 CeedCall(CeedOperatorLinearAssembleAddDiagonalComposite(op, request, false, assembled)); 2238 return CEED_ERROR_SUCCESS; 2239 } else { 2240 // Operator fallback 2241 CeedOperator op_fallback; 2242 2243 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleAddDiagonal\n"); 2244 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2245 if (op_fallback) { 2246 CeedCall(CeedOperatorLinearAssembleAddDiagonal(op_fallback, assembled, request)); 2247 return CEED_ERROR_SUCCESS; 2248 } 2249 } 2250 // Default interface implementation 2251 CeedCall(CeedOperatorLinearAssembleAddDiagonalSingle(op, request, false, assembled)); 2252 return CEED_ERROR_SUCCESS; 2253 } 2254 2255 /** 2256 @brief Fully assemble the point-block diagonal pattern of a linear `CeedOperator`. 2257 2258 Expected to be used in conjunction with @ref CeedOperatorLinearAssemblePointBlockDiagonal(). 2259 2260 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)`. 2261 Note that the `(i, j)` pairs are unique. 2262 This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemblePointBlockDiagonal() provides the values in the same ordering. 2263 2264 This will generally be slow unless your operator is low-order. 2265 2266 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2267 2268 @param[in] op `CeedOperator` to assemble 2269 @param[out] num_entries Number of entries in coordinate nonzero pattern 2270 @param[out] rows Row number for each entry 2271 @param[out] cols Column number for each entry 2272 2273 @ref User 2274 **/ 2275 int CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { 2276 bool is_composite; 2277 CeedInt num_active_components, num_sub_operators; 2278 CeedOperator *sub_operators; 2279 2280 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2281 2282 CeedSize input_size = 0, output_size = 0; 2283 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2284 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2285 2286 if (is_composite) { 2287 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_sub_operators)); 2288 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2289 } else { 2290 sub_operators = &op; 2291 num_sub_operators = 1; 2292 } 2293 2294 // Verify operator can be assembled correctly 2295 { 2296 CeedOperatorAssemblyData data; 2297 CeedInt num_active_elem_rstrs, comp_stride; 2298 CeedElemRestriction *active_elem_rstrs; 2299 2300 // Get initial values to check against 2301 CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[0], &data)); 2302 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL)); 2303 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[0], &comp_stride)); 2304 CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[0], &num_active_components)); 2305 2306 // Verify that all active element restrictions have same component stride and number of components 2307 for (CeedInt k = 0; k < num_sub_operators; k++) { 2308 CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[k], &data)); 2309 CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL)); 2310 for (CeedInt i = 0; i < num_active_elem_rstrs; i++) { 2311 CeedInt comp_stride_sub, num_active_components_sub; 2312 2313 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[i], &comp_stride_sub)); 2314 CeedCheck(comp_stride == comp_stride_sub, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, 2315 "Active element restrictions must have the same component stride: %d vs %d", comp_stride, comp_stride_sub); 2316 CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[i], &num_active_components_sub)); 2317 CeedCheck(num_active_components == num_active_components_sub, CeedOperatorReturnCeed(op), CEED_ERROR_INCOMPATIBLE, 2318 "All suboperators must have the same number of output components." 2319 " Previous: %" CeedInt_FMT " Current: %" CeedInt_FMT, 2320 num_active_components, num_active_components_sub); 2321 } 2322 } 2323 } 2324 *num_entries = input_size * num_active_components; 2325 CeedCall(CeedCalloc(*num_entries, rows)); 2326 CeedCall(CeedCalloc(*num_entries, cols)); 2327 2328 for (CeedInt o = 0; o < num_sub_operators; o++) { 2329 CeedElemRestriction active_elem_rstr, point_block_active_elem_rstr; 2330 CeedInt comp_stride, num_elem, elem_size; 2331 const CeedInt *offsets, *point_block_offsets; 2332 2333 CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[o], &active_elem_rstr)); 2334 CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstr, &comp_stride)); 2335 CeedCall(CeedElemRestrictionGetNumElements(active_elem_rstr, &num_elem)); 2336 CeedCall(CeedElemRestrictionGetElementSize(active_elem_rstr, &elem_size)); 2337 CeedCall(CeedElemRestrictionGetOffsets(active_elem_rstr, CEED_MEM_HOST, &offsets)); 2338 2339 CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstr, &point_block_active_elem_rstr)); 2340 CeedCall(CeedElemRestrictionGetOffsets(point_block_active_elem_rstr, CEED_MEM_HOST, &point_block_offsets)); 2341 2342 for (CeedSize i = 0; i < num_elem * elem_size; i++) { 2343 for (CeedInt c_out = 0; c_out < num_active_components; c_out++) { 2344 for (CeedInt c_in = 0; c_in < num_active_components; c_in++) { 2345 (*rows)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_out * comp_stride; 2346 (*cols)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_in * comp_stride; 2347 } 2348 } 2349 } 2350 2351 CeedCall(CeedElemRestrictionRestoreOffsets(active_elem_rstr, &offsets)); 2352 CeedCall(CeedElemRestrictionRestoreOffsets(point_block_active_elem_rstr, &point_block_offsets)); 2353 CeedCall(CeedElemRestrictionDestroy(&active_elem_rstr)); 2354 CeedCall(CeedElemRestrictionDestroy(&point_block_active_elem_rstr)); 2355 } 2356 return CEED_ERROR_SUCCESS; 2357 } 2358 2359 /** 2360 @brief Assemble the point block diagonal of a square linear `CeedOperator`. 2361 2362 This overwrites a `CeedVector` with the point block diagonal of a linear `CeedOperator`. 2363 2364 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2365 2366 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2367 2368 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2369 @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. 2370 The dimensions of this vector are derived from the active vector for the `CeedOperator`. 2371 The array has shape `[nodes, component out, component in]`. 2372 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2373 2374 @return An error code: 0 - success, otherwise - failure 2375 2376 @ref User 2377 **/ 2378 int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2379 bool is_composite; 2380 CeedSize input_size = 0, output_size = 0; 2381 2382 CeedCall(CeedOperatorCheckReady(op)); 2383 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2384 2385 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2386 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2387 2388 // Early exit for empty operator 2389 if (!is_composite) { 2390 CeedInt num_elem = 0; 2391 2392 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2393 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2394 } 2395 2396 if (op->LinearAssemblePointBlockDiagonal) { 2397 // Backend version 2398 CeedCall(op->LinearAssemblePointBlockDiagonal(op, assembled, request)); 2399 return CEED_ERROR_SUCCESS; 2400 } else if (op->LinearAssembleAddPointBlockDiagonal) { 2401 // Backend version with zeroing first 2402 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2403 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); 2404 return CEED_ERROR_SUCCESS; 2405 } else { 2406 // Operator fallback 2407 CeedOperator op_fallback; 2408 2409 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssemblePointBlockDiagonal\n"); 2410 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2411 if (op_fallback) { 2412 CeedCall(CeedOperatorLinearAssemblePointBlockDiagonal(op_fallback, assembled, request)); 2413 return CEED_ERROR_SUCCESS; 2414 } 2415 } 2416 // Default interface implementation 2417 CeedCall(CeedVectorSetValue(assembled, 0.0)); 2418 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); 2419 return CEED_ERROR_SUCCESS; 2420 } 2421 2422 /** 2423 @brief Assemble the point block diagonal of a square linear `CeedOperator`. 2424 2425 This sums into a `CeedVector` with the point block diagonal of a linear `CeedOperator`. 2426 2427 Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. 2428 2429 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2430 2431 @param[in] op `CeedOperator` to assemble `CeedQFunction` 2432 @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. 2433 The dimensions of this vector are derived from the active vector for the `CeedOperator`. 2434 The array has shape `[nodes, component out, component in]`. 2435 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2436 2437 @return An error code: 0 - success, otherwise - failure 2438 2439 @ref User 2440 **/ 2441 int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { 2442 bool is_composite; 2443 CeedSize input_size = 0, output_size = 0; 2444 2445 CeedCall(CeedOperatorCheckReady(op)); 2446 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2447 2448 CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); 2449 CeedCheck(input_size == output_size, CeedOperatorReturnCeed(op), CEED_ERROR_DIMENSION, "Operator must be square"); 2450 2451 // Early exit for empty operator 2452 if (!is_composite) { 2453 CeedInt num_elem = 0; 2454 2455 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2456 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2457 } 2458 2459 if (op->LinearAssembleAddPointBlockDiagonal) { 2460 // Backend version 2461 CeedCall(op->LinearAssembleAddPointBlockDiagonal(op, assembled, request)); 2462 return CEED_ERROR_SUCCESS; 2463 } else { 2464 // Operator fallback 2465 CeedOperator op_fallback; 2466 2467 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleAddPointBlockDiagonal\n"); 2468 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2469 if (op_fallback) { 2470 CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op_fallback, assembled, request)); 2471 return CEED_ERROR_SUCCESS; 2472 } 2473 } 2474 // Default interface implementation 2475 if (is_composite) { 2476 CeedCall(CeedOperatorLinearAssembleAddDiagonalComposite(op, request, true, assembled)); 2477 } else { 2478 CeedCall(CeedOperatorLinearAssembleAddDiagonalSingle(op, request, true, assembled)); 2479 } 2480 return CEED_ERROR_SUCCESS; 2481 } 2482 2483 /** 2484 @brief Fully assemble the nonzero pattern of a linear `CeedOperator`. 2485 2486 Expected to be used in conjunction with @ref CeedOperatorLinearAssemble(). 2487 2488 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)`. 2489 Note that the `(i, j)` pairs are not unique and may repeat. 2490 This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemble() provides the values in the same ordering. 2491 2492 This will generally be slow unless your operator is low-order. 2493 2494 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2495 2496 @param[in] op `CeedOperator` to assemble 2497 @param[out] num_entries Number of entries in coordinate nonzero pattern 2498 @param[out] rows Row number for each entry 2499 @param[out] cols Column number for each entry 2500 2501 @ref User 2502 **/ 2503 int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { 2504 bool is_composite; 2505 CeedInt num_suboperators, offset = 0; 2506 CeedSize single_entries; 2507 CeedOperator *sub_operators; 2508 2509 CeedCall(CeedOperatorCheckReady(op)); 2510 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2511 2512 if (op->LinearAssembleSymbolic) { 2513 // Backend version 2514 CeedCall(op->LinearAssembleSymbolic(op, num_entries, rows, cols)); 2515 return CEED_ERROR_SUCCESS; 2516 } else { 2517 // Operator fallback 2518 CeedOperator op_fallback; 2519 2520 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssembleSymbolic\n"); 2521 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2522 if (op_fallback) { 2523 CeedCall(CeedOperatorLinearAssembleSymbolic(op_fallback, num_entries, rows, cols)); 2524 return CEED_ERROR_SUCCESS; 2525 } 2526 } 2527 2528 // Default interface implementation 2529 2530 // Count entries and allocate rows, cols arrays 2531 CeedCall(CeedOperatorLinearAssembleGetNumEntries(op, num_entries)); 2532 CeedCall(CeedCalloc(*num_entries, rows)); 2533 CeedCall(CeedCalloc(*num_entries, cols)); 2534 2535 // Assemble nonzero locations 2536 if (is_composite) { 2537 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2538 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2539 for (CeedInt k = 0; k < num_suboperators; ++k) { 2540 CeedCall(CeedOperatorAssembleSymbolicSingle(sub_operators[k], offset, *rows, *cols)); 2541 CeedCall(CeedOperatorAssemblyCountEntriesSingle(sub_operators[k], &single_entries)); 2542 offset += single_entries; 2543 } 2544 } else { 2545 CeedCall(CeedOperatorAssembleSymbolicSingle(op, offset, *rows, *cols)); 2546 } 2547 return CEED_ERROR_SUCCESS; 2548 } 2549 2550 /** 2551 @brief Fully assemble the nonzero entries of a linear operator. 2552 2553 Expected to be used in conjunction with @ref CeedOperatorLinearAssembleSymbolic(). 2554 2555 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)`. 2556 Note that the `(i, j)` pairs are not unique and may repeat. 2557 This function returns the values of the nonzero entries to be added, their `(i, j)` locations are provided by @ref CeedOperatorLinearAssembleSymbolic(). 2558 2559 This will generally be slow unless your operator is low-order. 2560 2561 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2562 2563 @param[in] op `CeedOperator` to assemble 2564 @param[out] values Values to assemble into matrix 2565 2566 @ref User 2567 **/ 2568 int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { 2569 bool is_composite; 2570 CeedInt num_suboperators, offset = 0; 2571 CeedSize single_entries = 0; 2572 CeedOperator *sub_operators; 2573 2574 CeedCall(CeedOperatorCheckReady(op)); 2575 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 2576 2577 // Early exit for empty operator 2578 if (!is_composite) { 2579 CeedInt num_elem = 0; 2580 2581 CeedCall(CeedOperatorGetNumElements(op, &num_elem)); 2582 if (num_elem == 0) return CEED_ERROR_SUCCESS; 2583 } 2584 2585 if (op->LinearAssemble) { 2586 // Backend version 2587 CeedCall(op->LinearAssemble(op, values)); 2588 return CEED_ERROR_SUCCESS; 2589 } else if (is_composite) { 2590 // Default to summing contributions of suboperators 2591 CeedCall(CeedVectorSetValue(values, 0.0)); 2592 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2593 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2594 for (CeedInt k = 0; k < num_suboperators; k++) { 2595 CeedCall(CeedOperatorAssembleSingle(sub_operators[k], offset, values)); 2596 CeedCall(CeedOperatorAssemblyCountEntriesSingle(sub_operators[k], &single_entries)); 2597 offset += single_entries; 2598 } 2599 return CEED_ERROR_SUCCESS; 2600 } else if (op->LinearAssembleSingle) { 2601 CeedCall(CeedVectorSetValue(values, 0.0)); 2602 CeedCall(CeedOperatorAssembleSingle(op, offset, values)); 2603 return CEED_ERROR_SUCCESS; 2604 } else { 2605 // Operator fallback 2606 CeedOperator op_fallback; 2607 2608 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorLinearAssemble\n"); 2609 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2610 if (op_fallback) { 2611 CeedCall(CeedOperatorLinearAssemble(op_fallback, values)); 2612 return CEED_ERROR_SUCCESS; 2613 } 2614 } 2615 2616 // Default to interface version if non-composite and no fallback 2617 CeedCall(CeedVectorSetValue(values, 0.0)); 2618 CeedCall(CeedOperatorAssembleSingle(op, offset, values)); 2619 return CEED_ERROR_SUCCESS; 2620 } 2621 2622 /** 2623 @brief Get the multiplicity of nodes across sub-operators in a composite `CeedOperator`. 2624 2625 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2626 2627 @param[in] op Composite `CeedOperator` 2628 @param[in] num_skip_indices Number of sub-operators to skip 2629 @param[in] skip_indices Array of indices of sub-operators to skip 2630 @param[out] mult Vector to store multiplicity (of size `l_size` ) 2631 2632 @return An error code: 0 - success, otherwise - failure 2633 2634 @ref User 2635 **/ 2636 int CeedOperatorCompositeGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult) { 2637 Ceed ceed; 2638 CeedInt num_suboperators; 2639 CeedSize l_vec_len; 2640 CeedScalar *mult_array; 2641 CeedVector ones_l_vec; 2642 CeedElemRestriction elem_rstr, mult_elem_rstr; 2643 CeedOperator *sub_operators; 2644 2645 CeedCall(CeedOperatorCheckReady(op)); 2646 2647 // Zero mult vector 2648 CeedCall(CeedVectorSetValue(mult, 0.0)); 2649 2650 // Get suboperators 2651 CeedCall(CeedOperatorCompositeGetNumSub(op, &num_suboperators)); 2652 if (num_suboperators == 0) return CEED_ERROR_SUCCESS; 2653 CeedCall(CeedOperatorCompositeGetSubList(op, &sub_operators)); 2654 2655 // Work vector 2656 CeedCall(CeedVectorGetLength(mult, &l_vec_len)); 2657 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2658 CeedCall(CeedVectorCreate(ceed, l_vec_len, &ones_l_vec)); 2659 CeedCall(CeedDestroy(&ceed)); 2660 CeedCall(CeedVectorSetValue(ones_l_vec, 1.0)); 2661 CeedCall(CeedVectorGetArray(mult, CEED_MEM_HOST, &mult_array)); 2662 2663 // Compute multiplicity across suboperators 2664 for (CeedInt i = 0; i < num_suboperators; i++) { 2665 const CeedScalar *sub_mult_array; 2666 CeedVector sub_mult_l_vec, ones_e_vec; 2667 2668 // -- Check for suboperator to skip 2669 for (CeedInt j = 0; j < num_skip_indices; j++) { 2670 if (skip_indices[j] == i) continue; 2671 } 2672 2673 // -- Sub operator multiplicity 2674 CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[i], &elem_rstr)); 2675 CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr, &mult_elem_rstr)); 2676 CeedCall(CeedElemRestrictionDestroy(&elem_rstr)); 2677 CeedCall(CeedElemRestrictionCreateVector(mult_elem_rstr, &sub_mult_l_vec, &ones_e_vec)); 2678 CeedCall(CeedVectorSetValue(sub_mult_l_vec, 0.0)); 2679 CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_NOTRANSPOSE, ones_l_vec, ones_e_vec, CEED_REQUEST_IMMEDIATE)); 2680 CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_TRANSPOSE, ones_e_vec, sub_mult_l_vec, CEED_REQUEST_IMMEDIATE)); 2681 CeedCall(CeedVectorGetArrayRead(sub_mult_l_vec, CEED_MEM_HOST, &sub_mult_array)); 2682 // ---- Flag every node present in the current suboperator 2683 for (CeedSize j = 0; j < l_vec_len; j++) { 2684 if (sub_mult_array[j] > 0.0) mult_array[j] += 1.0; 2685 } 2686 CeedCall(CeedVectorRestoreArrayRead(sub_mult_l_vec, &sub_mult_array)); 2687 CeedCall(CeedVectorDestroy(&sub_mult_l_vec)); 2688 CeedCall(CeedVectorDestroy(&ones_e_vec)); 2689 CeedCall(CeedElemRestrictionDestroy(&mult_elem_rstr)); 2690 } 2691 CeedCall(CeedVectorRestoreArray(mult, &mult_array)); 2692 CeedCall(CeedVectorDestroy(&ones_l_vec)); 2693 return CEED_ERROR_SUCCESS; 2694 } 2695 2696 /** 2697 @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator`, creating the prolongation basis from the fine and coarse grid interpolation. 2698 2699 Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. 2700 2701 @param[in] op_fine Fine grid `CeedOperator` 2702 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 2703 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 2704 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 2705 @param[out] op_coarse Coarse grid `CeedOperator` 2706 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 2707 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 2708 2709 @return An error code: 0 - success, otherwise - failure 2710 2711 @ref User 2712 **/ 2713 int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2714 CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { 2715 CeedBasis basis_c_to_f = NULL; 2716 2717 CeedCall(CeedOperatorCheckReady(op_fine)); 2718 2719 // Build prolongation matrix, if required 2720 if (op_prolong || op_restrict) { 2721 CeedBasis basis_fine; 2722 2723 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2724 CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f)); 2725 CeedCall(CeedBasisDestroy(&basis_fine)); 2726 } 2727 2728 // Core code 2729 CeedCall(CeedOperatorMultigridLevelCreateSingle_Core(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, 2730 op_restrict)); 2731 return CEED_ERROR_SUCCESS; 2732 } 2733 2734 /** 2735 @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a tensor basis for the active basis. 2736 2737 Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. 2738 2739 @param[in] op_fine Fine grid `CeedOperator` 2740 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 2741 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 2742 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 2743 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator` 2744 @param[out] op_coarse Coarse grid `CeedOperator` 2745 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 2746 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 2747 2748 @return An error code: 0 - success, otherwise - failure 2749 2750 @ref User 2751 **/ 2752 int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2753 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 2754 CeedOperator *op_restrict) { 2755 Ceed ceed; 2756 CeedInt Q_f, Q_c; 2757 CeedBasis basis_fine, basis_c_to_f = NULL; 2758 2759 CeedCall(CeedOperatorCheckReady(op_fine)); 2760 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 2761 2762 // Check for compatible quadrature spaces 2763 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2764 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); 2765 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); 2766 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, 2767 "Bases must have compatible quadrature spaces." 2768 " Fine grid: %" CeedInt_FMT " points, Coarse grid: %" CeedInt_FMT " points", 2769 Q_f, Q_c); 2770 2771 // Create coarse to fine basis, if required 2772 if (op_prolong || op_restrict) { 2773 CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c; 2774 CeedScalar *q_ref, *q_weight, *grad; 2775 2776 // Check if interpolation matrix is provided 2777 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE, 2778 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix"); 2779 CeedCall(CeedBasisGetDimension(basis_fine, &dim)); 2780 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); 2781 CeedCall(CeedBasisGetNumNodes1D(basis_fine, &P_1d_f)); 2782 CeedCall(CeedBasisDestroy(&basis_fine)); 2783 CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); 2784 P_1d_c = dim == 1 ? num_nodes_c : dim == 2 ? sqrt(num_nodes_c) : cbrt(num_nodes_c); 2785 CeedCall(CeedCalloc(P_1d_f, &q_ref)); 2786 CeedCall(CeedCalloc(P_1d_f, &q_weight)); 2787 CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad)); 2788 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)); 2789 CeedCall(CeedFree(&q_ref)); 2790 CeedCall(CeedFree(&q_weight)); 2791 CeedCall(CeedFree(&grad)); 2792 } 2793 2794 // Core code 2795 CeedCall(CeedOperatorMultigridLevelCreateSingle_Core(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, 2796 op_restrict)); 2797 CeedCall(CeedDestroy(&ceed)); 2798 return CEED_ERROR_SUCCESS; 2799 } 2800 2801 /** 2802 @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a non-tensor basis for the active vector 2803 2804 Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. 2805 2806 @param[in] op_fine Fine grid `CeedOperator` 2807 @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` 2808 @param[in] rstr_coarse Coarse grid `CeedElemRestriction` 2809 @param[in] basis_coarse Coarse grid active vector `CeedBasis` 2810 @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator` 2811 @param[out] op_coarse Coarse grid `CeedOperator` 2812 @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` 2813 @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` 2814 2815 @return An error code: 0 - success, otherwise - failure 2816 2817 @ref User 2818 **/ 2819 int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, 2820 const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, 2821 CeedOperator *op_restrict) { 2822 Ceed ceed; 2823 CeedInt Q_f, Q_c; 2824 CeedBasis basis_fine, basis_c_to_f = NULL; 2825 2826 CeedCall(CeedOperatorCheckReady(op_fine)); 2827 CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); 2828 2829 // Check for compatible quadrature spaces 2830 CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); 2831 CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); 2832 CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); 2833 CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces"); 2834 2835 // Coarse to fine basis 2836 if (op_prolong || op_restrict) { 2837 CeedInt dim, num_comp, num_nodes_c, num_nodes_f; 2838 CeedScalar *q_ref, *q_weight, *grad; 2839 CeedElemTopology topo; 2840 2841 // Check if interpolation matrix is provided 2842 CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE, 2843 "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix"); 2844 CeedCall(CeedBasisGetTopology(basis_fine, &topo)); 2845 CeedCall(CeedBasisGetDimension(basis_fine, &dim)); 2846 CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); 2847 CeedCall(CeedBasisGetNumNodes(basis_fine, &num_nodes_f)); 2848 CeedCall(CeedBasisDestroy(&basis_fine)); 2849 CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); 2850 CeedCall(CeedCalloc(num_nodes_f * dim, &q_ref)); 2851 CeedCall(CeedCalloc(num_nodes_f, &q_weight)); 2852 CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad)); 2853 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)); 2854 CeedCall(CeedFree(&q_ref)); 2855 CeedCall(CeedFree(&q_weight)); 2856 CeedCall(CeedFree(&grad)); 2857 } 2858 2859 // Core code 2860 CeedCall(CeedOperatorMultigridLevelCreateSingle_Core(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, 2861 op_restrict)); 2862 CeedCall(CeedDestroy(&ceed)); 2863 return CEED_ERROR_SUCCESS; 2864 } 2865 2866 /** 2867 @brief Build a FDM based approximate inverse for each element for a `CeedOperator`. 2868 2869 This returns a `CeedOperator` and `CeedVector` to apply a Fast Diagonalization Method based approximate inverse. 2870 This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, \f$M = V^T V, K = V^T S V\f$. 2871 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$. 2872 The `CeedOperator` must be linear and non-composite. 2873 The associated `CeedQFunction` must therefore also be linear. 2874 2875 Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. 2876 2877 @param[in] op `CeedOperator` to create element inverses 2878 @param[out] fdm_inv `CeedOperator` to apply the action of a FDM based inverse for each element 2879 @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE 2880 2881 @return An error code: 0 - success, otherwise - failure 2882 2883 @ref User 2884 **/ 2885 int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, CeedRequest *request) { 2886 Ceed ceed, ceed_parent; 2887 bool interp = false, grad = false, is_tensor_basis = true; 2888 CeedInt num_input_fields, P_1d, Q_1d, num_nodes, num_qpts, dim, num_comp = 1, num_elem = 1; 2889 CeedScalar *mass, *laplace, *x, *fdm_interp, *lambda, *elem_avg; 2890 const CeedScalar *interp_1d, *grad_1d, *q_weight_1d; 2891 CeedVector q_data; 2892 CeedElemRestriction rstr = NULL, rstr_qd_i; 2893 CeedBasis basis = NULL, fdm_basis; 2894 CeedQFunctionContext ctx_fdm; 2895 CeedQFunctionField *qf_fields; 2896 CeedQFunction qf, qf_fdm; 2897 CeedOperatorField *op_fields; 2898 2899 CeedCall(CeedOperatorCheckReady(op)); 2900 2901 if (op->CreateFDMElementInverse) { 2902 // Backend version 2903 CeedCall(op->CreateFDMElementInverse(op, fdm_inv, request)); 2904 return CEED_ERROR_SUCCESS; 2905 } else { 2906 // Operator fallback 2907 CeedOperator op_fallback; 2908 2909 CeedDebug(CeedOperatorReturnCeed(op), "\nFalling back for CeedOperatorCreateFDMElementInverse\n"); 2910 CeedCall(CeedOperatorGetFallback(op, &op_fallback)); 2911 if (op_fallback) { 2912 CeedCall(CeedOperatorCreateFDMElementInverse(op_fallback, fdm_inv, request)); 2913 return CEED_ERROR_SUCCESS; 2914 } 2915 } 2916 2917 // Default interface implementation 2918 CeedCall(CeedOperatorGetCeed(op, &ceed)); 2919 CeedCall(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); 2920 CeedCall(CeedOperatorGetQFunction(op, &qf)); 2921 2922 // Determine active input basis 2923 CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_fields, NULL, NULL)); 2924 CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); 2925 for (CeedInt i = 0; i < num_input_fields; i++) { 2926 CeedVector vec; 2927 2928 CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); 2929 if (vec == CEED_VECTOR_ACTIVE) { 2930 CeedEvalMode eval_mode; 2931 2932 CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); 2933 interp = interp || eval_mode == CEED_EVAL_INTERP; 2934 grad = grad || eval_mode == CEED_EVAL_GRAD; 2935 if (!basis) CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis)); 2936 if (!rstr) CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr)); 2937 } 2938 CeedCall(CeedVectorDestroy(&vec)); 2939 } 2940 CeedCheck(basis, ceed, CEED_ERROR_BACKEND, "No active field set"); 2941 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 2942 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 2943 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 2944 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 2945 CeedCall(CeedBasisGetDimension(basis, &dim)); 2946 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 2947 CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 2948 2949 // Build and diagonalize 1D Mass and Laplacian 2950 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 2951 CeedCheck(is_tensor_basis, ceed, CEED_ERROR_BACKEND, "FDMElementInverse only supported for tensor bases"); 2952 CeedCall(CeedCalloc(P_1d * P_1d, &mass)); 2953 CeedCall(CeedCalloc(P_1d * P_1d, &laplace)); 2954 CeedCall(CeedCalloc(P_1d * P_1d, &x)); 2955 CeedCall(CeedCalloc(P_1d * P_1d, &fdm_interp)); 2956 CeedCall(CeedCalloc(P_1d, &lambda)); 2957 // -- Build matrices 2958 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 2959 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 2960 CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d)); 2961 CeedCall(CeedBuildMassLaplace(interp_1d, grad_1d, q_weight_1d, P_1d, Q_1d, dim, mass, laplace)); 2962 2963 // -- Diagonalize 2964 CeedCall(CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d)); 2965 CeedCall(CeedFree(&mass)); 2966 CeedCall(CeedFree(&laplace)); 2967 for (CeedInt i = 0; i < P_1d; i++) { 2968 for (CeedInt j = 0; j < P_1d; j++) fdm_interp[i + j * P_1d] = x[j + i * P_1d]; 2969 } 2970 CeedCall(CeedFree(&x)); 2971 2972 { 2973 CeedInt layout[3], num_modes = (interp ? 1 : 0) + (grad ? dim : 0); 2974 CeedScalar max_norm = 0; 2975 const CeedScalar *assembled_array, *q_weight_array; 2976 CeedVector assembled = NULL, q_weight; 2977 CeedElemRestriction rstr_qf = NULL; 2978 2979 // Assemble QFunction 2980 CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled, &rstr_qf, request)); 2981 CeedCall(CeedElemRestrictionGetELayout(rstr_qf, layout)); 2982 CeedCall(CeedElemRestrictionDestroy(&rstr_qf)); 2983 CeedCall(CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm)); 2984 2985 // Calculate element averages 2986 CeedCall(CeedVectorCreate(ceed_parent, num_qpts, &q_weight)); 2987 CeedCall(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_weight)); 2988 CeedCall(CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array)); 2989 CeedCall(CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array)); 2990 CeedCall(CeedCalloc(num_elem, &elem_avg)); 2991 const CeedScalar qf_value_bound = max_norm * 100 * CEED_EPSILON; 2992 2993 for (CeedInt e = 0; e < num_elem; e++) { 2994 CeedInt count = 0; 2995 2996 for (CeedInt q = 0; q < num_qpts; q++) { 2997 for (CeedInt i = 0; i < num_comp * num_comp * num_modes * num_modes; i++) { 2998 if (fabs(assembled_array[q * layout[0] + i * layout[1] + e * layout[2]]) > qf_value_bound) { 2999 elem_avg[e] += assembled_array[q * layout[0] + i * layout[1] + e * layout[2]] / q_weight_array[q]; 3000 count++; 3001 } 3002 } 3003 } 3004 if (count) { 3005 elem_avg[e] /= count; 3006 } else { 3007 elem_avg[e] = 1.0; 3008 } 3009 } 3010 CeedCall(CeedVectorRestoreArrayRead(assembled, &assembled_array)); 3011 CeedCall(CeedVectorDestroy(&assembled)); 3012 CeedCall(CeedVectorRestoreArrayRead(q_weight, &q_weight_array)); 3013 CeedCall(CeedVectorDestroy(&q_weight)); 3014 } 3015 3016 // Build FDM diagonal 3017 { 3018 CeedScalar *q_data_array, *fdm_diagonal; 3019 3020 CeedCall(CeedCalloc(num_comp * num_nodes, &fdm_diagonal)); 3021 const CeedScalar fdm_diagonal_bound = num_nodes * CEED_EPSILON; 3022 for (CeedInt c = 0; c < num_comp; c++) { 3023 for (CeedInt n = 0; n < num_nodes; n++) { 3024 if (interp) fdm_diagonal[c * num_nodes + n] = 1.0; 3025 if (grad) { 3026 for (CeedInt d = 0; d < dim; d++) { 3027 CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d; 3028 fdm_diagonal[c * num_nodes + n] += lambda[i]; 3029 } 3030 } 3031 if (fabs(fdm_diagonal[c * num_nodes + n]) < fdm_diagonal_bound) fdm_diagonal[c * num_nodes + n] = fdm_diagonal_bound; 3032 } 3033 } 3034 CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * num_nodes, &q_data)); 3035 CeedCall(CeedVectorSetValue(q_data, 0.0)); 3036 CeedCall(CeedVectorGetArrayWrite(q_data, CEED_MEM_HOST, &q_data_array)); 3037 for (CeedInt e = 0; e < num_elem; e++) { 3038 for (CeedInt c = 0; c < num_comp; c++) { 3039 for (CeedInt n = 0; n < num_nodes; n++) { 3040 q_data_array[(e * num_comp + c) * num_nodes + n] = 1. / (elem_avg[e] * fdm_diagonal[c * num_nodes + n]); 3041 } 3042 } 3043 } 3044 CeedCall(CeedFree(&elem_avg)); 3045 CeedCall(CeedFree(&fdm_diagonal)); 3046 CeedCall(CeedVectorRestoreArray(q_data, &q_data_array)); 3047 } 3048 3049 // Setup FDM operator 3050 // -- Basis 3051 { 3052 CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy; 3053 3054 CeedCall(CeedCalloc(P_1d * P_1d, &grad_dummy)); 3055 CeedCall(CeedCalloc(P_1d, &q_ref_dummy)); 3056 CeedCall(CeedCalloc(P_1d, &q_weight_dummy)); 3057 CeedCall(CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, fdm_interp, grad_dummy, q_ref_dummy, q_weight_dummy, &fdm_basis)); 3058 CeedCall(CeedFree(&fdm_interp)); 3059 CeedCall(CeedFree(&grad_dummy)); 3060 CeedCall(CeedFree(&q_ref_dummy)); 3061 CeedCall(CeedFree(&q_weight_dummy)); 3062 CeedCall(CeedFree(&lambda)); 3063 } 3064 3065 // -- Restriction 3066 { 3067 CeedInt strides[3] = {1, num_nodes, num_nodes * num_comp}; 3068 CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, num_nodes, num_comp, 3069 (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes, strides, &rstr_qd_i)); 3070 } 3071 3072 // -- QFunction 3073 CeedCall(CeedQFunctionCreateInteriorByName(ceed_parent, "Scale", &qf_fdm)); 3074 CeedCall(CeedQFunctionAddInput(qf_fdm, "input", num_comp, CEED_EVAL_INTERP)); 3075 CeedCall(CeedQFunctionAddInput(qf_fdm, "scale", num_comp, CEED_EVAL_NONE)); 3076 CeedCall(CeedQFunctionAddOutput(qf_fdm, "output", num_comp, CEED_EVAL_INTERP)); 3077 CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_fdm, num_comp)); 3078 3079 // -- QFunction context 3080 { 3081 CeedInt *num_comp_data; 3082 3083 CeedCall(CeedCalloc(1, &num_comp_data)); 3084 num_comp_data[0] = num_comp; 3085 CeedCall(CeedQFunctionContextCreate(ceed, &ctx_fdm)); 3086 CeedCall(CeedQFunctionContextSetData(ctx_fdm, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_data), num_comp_data)); 3087 } 3088 CeedCall(CeedQFunctionSetContext(qf_fdm, ctx_fdm)); 3089 CeedCall(CeedQFunctionContextDestroy(&ctx_fdm)); 3090 3091 // -- Operator 3092 CeedCall(CeedOperatorCreate(ceed_parent, qf_fdm, NULL, NULL, fdm_inv)); 3093 CeedCall(CeedOperatorSetField(*fdm_inv, "input", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); 3094 CeedCall(CeedOperatorSetField(*fdm_inv, "scale", rstr_qd_i, CEED_BASIS_NONE, q_data)); 3095 CeedCall(CeedOperatorSetField(*fdm_inv, "output", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); 3096 3097 // Cleanup 3098 CeedCall(CeedDestroy(&ceed)); 3099 CeedCall(CeedDestroy(&ceed_parent)); 3100 CeedCall(CeedVectorDestroy(&q_data)); 3101 CeedCall(CeedElemRestrictionDestroy(&rstr)); 3102 CeedCall(CeedElemRestrictionDestroy(&rstr_qd_i)); 3103 CeedCall(CeedBasisDestroy(&basis)); 3104 CeedCall(CeedBasisDestroy(&fdm_basis)); 3105 CeedCall(CeedQFunctionDestroy(&qf)); 3106 CeedCall(CeedQFunctionDestroy(&qf_fdm)); 3107 return CEED_ERROR_SUCCESS; 3108 } 3109 3110 /// @} 3111