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