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