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