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