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