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