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