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