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