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