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