// Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. // // SPDX-License-Identifier: BSD-2-Clause // // This file is part of CEED: http://github.com/ceed #include #include #include #include #include #include #include #include /// @file /// Implementation of CeedOperator preconditioning interfaces /// ---------------------------------------------------------------------------- /// CeedOperator Library Internal Preconditioning Functions /// ---------------------------------------------------------------------------- /// @addtogroup CeedOperatorDeveloper /// @{ /** @brief Duplicate a `CeedQFunction` with a reference `Ceed` to fallback for advanced `CeedOperator` functionality @param[in] fallback_ceed `Ceed` on which to create fallback `CeedQFunction` @param[in] qf `CeedQFunction` to create fallback for @param[out] qf_fallback Fallback `CeedQFunction` @return An error code: 0 - success, otherwise - failure @ref Developer **/ static int CeedQFunctionCreateFallback(Ceed fallback_ceed, CeedQFunction qf, CeedQFunction *qf_fallback) { char *source_path_with_name = NULL; CeedInt num_input_fields, num_output_fields; Ceed ceed; CeedQFunctionField *input_fields, *output_fields; // Check if NULL qf passed in if (!qf) return CEED_ERROR_SUCCESS; CeedCall(CeedQFunctionGetCeed(qf, &ceed)); CeedDebug256(ceed, 1, "---------- CeedOperator Fallback ----------\n"); CeedDebug(ceed, "Creating fallback CeedQFunction\n"); if (qf->source_path) { size_t path_len = strlen(qf->source_path), name_len = strlen(qf->kernel_name); CeedCall(CeedCalloc(path_len + name_len + 2, &source_path_with_name)); memcpy(source_path_with_name, qf->source_path, path_len); memcpy(&source_path_with_name[path_len], ":", 1); memcpy(&source_path_with_name[path_len + 1], qf->kernel_name, name_len); } else { CeedCall(CeedCalloc(1, &source_path_with_name)); } { CeedInt vec_length; CeedQFunctionUser f; CeedCall(CeedQFunctionGetVectorLength(qf, &vec_length)); CeedCall(CeedQFunctionGetUserFunction(qf, &f)); CeedCall(CeedQFunctionCreateInterior(fallback_ceed, vec_length, f, source_path_with_name, qf_fallback)); } { CeedQFunctionContext ctx; CeedCall(CeedQFunctionGetContext(qf, &ctx)); CeedCall(CeedQFunctionSetContext(*qf_fallback, ctx)); } CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); for (CeedInt i = 0; i < num_input_fields; i++) { const char *field_name; CeedInt size; CeedEvalMode eval_mode; CeedCall(CeedQFunctionFieldGetData(input_fields[i], &field_name, &size, &eval_mode)); CeedCall(CeedQFunctionAddInput(*qf_fallback, field_name, size, eval_mode)); } for (CeedInt i = 0; i < num_output_fields; i++) { const char *field_name; CeedInt size; CeedEvalMode eval_mode; CeedCall(CeedQFunctionFieldGetData(output_fields[i], &field_name, &size, &eval_mode)); CeedCall(CeedQFunctionAddOutput(*qf_fallback, field_name, size, eval_mode)); } CeedCall(CeedFree(&source_path_with_name)); return CEED_ERROR_SUCCESS; } /** @brief Duplicate a `CeedOperator` with a reference `Ceed` to fallback for advanced `CeedOperator` functionality @param[in,out] op `CeedOperator` to create fallback for @return An error code: 0 - success, otherwise - failure @ref Developer **/ static int CeedOperatorCreateFallback(CeedOperator op) { bool is_composite; Ceed ceed, ceed_fallback; CeedOperator op_fallback; // Check not already created if (op->op_fallback) return CEED_ERROR_SUCCESS; // Fallback Ceed CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback)); if (!ceed_fallback) return CEED_ERROR_SUCCESS; CeedDebug256(ceed, 1, "---------- CeedOperator Fallback ----------\n"); CeedDebug(ceed, "Creating fallback CeedOperator\n"); // Clone Op CeedCall(CeedOperatorIsComposite(op, &is_composite)); if (is_composite) { CeedInt num_suboperators; CeedOperator *sub_operators; CeedCall(CeedCompositeOperatorCreate(ceed_fallback, &op_fallback)); CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); for (CeedInt i = 0; i < num_suboperators; i++) { CeedOperator op_sub_fallback; CeedCall(CeedOperatorGetFallback(sub_operators[i], &op_sub_fallback)); CeedCall(CeedCompositeOperatorAddSub(op_fallback, op_sub_fallback)); } } else { CeedInt num_input_fields, num_output_fields; CeedQFunction qf_fallback = NULL, dqf_fallback = NULL, dqfT_fallback = NULL; CeedOperatorField *input_fields, *output_fields; CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->qf, &qf_fallback)); CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqf, &dqf_fallback)); CeedCall(CeedQFunctionCreateFallback(ceed_fallback, op->dqfT, &dqfT_fallback)); CeedCall(CeedOperatorCreate(ceed_fallback, qf_fallback, dqf_fallback, dqfT_fallback, &op_fallback)); CeedCall(CeedOperatorGetFields(op, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); for (CeedInt i = 0; i < num_input_fields; i++) { const char *field_name; CeedVector vec; CeedElemRestriction rstr; CeedBasis basis; CeedCall(CeedOperatorFieldGetData(input_fields[i], &field_name, &rstr, &basis, &vec)); CeedCall(CeedOperatorSetField(op_fallback, field_name, rstr, basis, vec)); } for (CeedInt i = 0; i < num_output_fields; i++) { const char *field_name; CeedVector vec; CeedElemRestriction rstr; CeedBasis basis; CeedCall(CeedOperatorFieldGetData(output_fields[i], &field_name, &rstr, &basis, &vec)); CeedCall(CeedOperatorSetField(op_fallback, field_name, rstr, basis, vec)); } CeedCall(CeedQFunctionAssemblyDataReferenceCopy(op->qf_assembled, &op_fallback->qf_assembled)); // Cleanup CeedCall(CeedQFunctionDestroy(&qf_fallback)); CeedCall(CeedQFunctionDestroy(&dqf_fallback)); CeedCall(CeedQFunctionDestroy(&dqfT_fallback)); } CeedCall(CeedOperatorSetName(op_fallback, op->name)); CeedCall(CeedOperatorCheckReady(op_fallback)); // Note: No ref-counting here so we don't get caught in a reference loop. // The op holds the only reference to op_fallback and is responsible for deleting itself and op_fallback. op->op_fallback = op_fallback; op_fallback->op_fallback_parent = op; return CEED_ERROR_SUCCESS; } /** @brief Core logic for assembling operator diagonal or point block diagonal @param[in] op `CeedOperator` to assemble point block diagonal @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @param[in] is_point_block Boolean flag to assemble diagonal or point block diagonal @param[out] assembled `CeedVector` to store assembled diagonal @return An error code: 0 - success, otherwise - failure @ref Developer **/ static inline int CeedSingleOperatorAssembleAddDiagonal_Core(CeedOperator op, CeedRequest *request, const bool is_point_block, CeedVector assembled) { Ceed ceed; bool is_composite; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); // Assemble QFunction CeedInt layout_qf[3]; const CeedScalar *assembled_qf_array; CeedVector assembled_qf = NULL; CeedElemRestriction assembled_elem_rstr = NULL; CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, request)); CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, layout_qf)); CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr)); CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); // Get assembly data const CeedEvalMode **eval_modes_in, **eval_modes_out; CeedInt num_active_bases_in, *num_eval_modes_in, num_active_bases_out, *num_eval_modes_out; CeedSize **eval_mode_offsets_in, **eval_mode_offsets_out, num_output_components; CeedBasis *active_bases_in, *active_bases_out; CeedElemRestriction *active_elem_rstrs_in, *active_elem_rstrs_out; CeedOperatorAssemblyData data; CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases_in, &num_eval_modes_in, &eval_modes_in, &eval_mode_offsets_in, &num_active_bases_out, &num_eval_modes_out, &eval_modes_out, &eval_mode_offsets_out, &num_output_components)); CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases_in, NULL, NULL, &active_bases_out, NULL)); CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, NULL, &active_elem_rstrs_in, NULL, &active_elem_rstrs_out)); // Loop over all active bases (find matching input/output pairs) for (CeedInt b = 0; b < CeedIntMin(num_active_bases_in, num_active_bases_out); b++) { CeedInt b_in, b_out, num_elem, num_nodes, num_qpts, num_comp; bool has_eval_none = false; CeedScalar *elem_diag_array, *identity = NULL; CeedVector elem_diag; CeedElemRestriction diag_elem_rstr; if (num_active_bases_in <= num_active_bases_out) { b_in = b; for (b_out = 0; b_out < num_active_bases_out; b_out++) { if (active_bases_in[b_in] == active_bases_out[b_out]) { break; } } if (b_out == num_active_bases_out) { continue; } // No matching output basis found } else { b_out = b; for (b_in = 0; b_in < num_active_bases_in; b_in++) { if (active_bases_in[b_in] == active_bases_out[b_out]) { break; } } if (b_in == num_active_bases_in) { continue; } // No matching output basis found } CeedCheck(active_elem_rstrs_in[b_in] == active_elem_rstrs_out[b_out], ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator diagonal with different input and output active element restrictions"); // Assemble point block diagonal restriction, if needed if (is_point_block) { CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstrs_in[b_in], &diag_elem_rstr)); } else { CeedCall(CeedElemRestrictionCreateUnsignedCopy(active_elem_rstrs_in[b_in], &diag_elem_rstr)); } // Create diagonal vector CeedCall(CeedElemRestrictionCreateVector(diag_elem_rstr, NULL, &elem_diag)); // Assemble element operator diagonals CeedCall(CeedVectorSetValue(elem_diag, 0.0)); CeedCall(CeedVectorGetArray(elem_diag, CEED_MEM_HOST, &elem_diag_array)); CeedCall(CeedElemRestrictionGetNumElements(diag_elem_rstr, &num_elem)); CeedCall(CeedBasisGetNumNodes(active_bases_in[b_in], &num_nodes)); CeedCall(CeedBasisGetNumComponents(active_bases_in[b_in], &num_comp)); if (active_bases_in[b_in] == CEED_BASIS_NONE) num_qpts = num_nodes; else CeedCall(CeedBasisGetNumQuadraturePoints(active_bases_in[b_in], &num_qpts)); // Construct identity matrix for basis if required for (CeedInt i = 0; i < num_eval_modes_in[b_in]; i++) { has_eval_none = has_eval_none || (eval_modes_in[b_in][i] == CEED_EVAL_NONE); } for (CeedInt i = 0; i < num_eval_modes_out[b_out]; i++) { has_eval_none = has_eval_none || (eval_modes_out[b_out][i] == CEED_EVAL_NONE); } if (has_eval_none) { CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) identity[i * num_nodes + i] = 1.0; } // Compute the diagonal of B^T D B // Each element for (CeedSize e = 0; e < num_elem; e++) { // Each basis eval mode pair CeedInt d_out = 0, q_comp_out; CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; for (CeedInt e_out = 0; e_out < num_eval_modes_out[b_out]; e_out++) { CeedInt d_in = 0, q_comp_in; const CeedScalar *B_t = NULL; CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; CeedCall(CeedOperatorGetBasisPointer(active_bases_out[b_out], eval_modes_out[b_out][e_out], identity, &B_t)); CeedCall(CeedBasisGetNumQuadratureComponents(active_bases_out[b_out], eval_modes_out[b_out][e_out], &q_comp_out)); if (q_comp_out > 1) { if (e_out == 0 || eval_modes_out[b_out][e_out] != eval_mode_out_prev) d_out = 0; else B_t = &B_t[(++d_out) * num_qpts * num_nodes]; } eval_mode_out_prev = eval_modes_out[b_out][e_out]; for (CeedInt e_in = 0; e_in < num_eval_modes_in[b_in]; e_in++) { const CeedScalar *B = NULL; CeedCall(CeedOperatorGetBasisPointer(active_bases_in[b_in], eval_modes_in[b_in][e_in], identity, &B)); CeedCall(CeedBasisGetNumQuadratureComponents(active_bases_in[b_in], eval_modes_in[b_in][e_in], &q_comp_in)); if (q_comp_in > 1) { if (e_in == 0 || eval_modes_in[b_in][e_in] != eval_mode_in_prev) d_in = 0; else B = &B[(++d_in) * num_qpts * num_nodes]; } eval_mode_in_prev = eval_modes_in[b_in][e_in]; // Each component for (CeedInt c_out = 0; c_out < num_comp; c_out++) { // Each qpt/node pair for (CeedInt q = 0; q < num_qpts; q++) { if (is_point_block) { // Point Block Diagonal for (CeedInt c_in = 0; c_in < num_comp; c_in++) { const CeedSize c_offset = (eval_mode_offsets_in[b_in][e_in] + c_in) * num_output_components + eval_mode_offsets_out[b_out][e_out] + c_out; const CeedScalar qf_value = assembled_qf_array[q * layout_qf[0] + c_offset * layout_qf[1] + e * layout_qf[2]]; for (CeedInt n = 0; n < num_nodes; n++) { elem_diag_array[((e * num_comp + c_out) * num_comp + c_in) * num_nodes + n] += B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n]; } } } else { // Diagonal Only const CeedInt c_offset = (eval_mode_offsets_in[b_in][e_in] + c_out) * num_output_components + eval_mode_offsets_out[b_out][e_out] + c_out; const CeedScalar qf_value = assembled_qf_array[q * layout_qf[0] + c_offset * layout_qf[1] + e * layout_qf[2]]; for (CeedInt n = 0; n < num_nodes; n++) { elem_diag_array[(e * num_comp + c_out) * num_nodes + n] += B_t[q * num_nodes + n] * qf_value * B[q * num_nodes + n]; } } } } } } } CeedCall(CeedVectorRestoreArray(elem_diag, &elem_diag_array)); // Assemble local operator diagonal CeedCall(CeedElemRestrictionApply(diag_elem_rstr, CEED_TRANSPOSE, elem_diag, assembled, request)); // Cleanup CeedCall(CeedElemRestrictionDestroy(&diag_elem_rstr)); CeedCall(CeedVectorDestroy(&elem_diag)); CeedCall(CeedFree(&identity)); } CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); CeedCall(CeedVectorDestroy(&assembled_qf)); return CEED_ERROR_SUCCESS; } /** @brief Core logic for assembling composite operator diagonal @param[in] op `CeedOperator` to assemble point block diagonal @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @param[in] is_point_block Boolean flag to assemble diagonal or point block diagonal @param[out] assembled `CeedVector` to store assembled diagonal @return An error code: 0 - success, otherwise - failure @ref Developer **/ static inline int CeedCompositeOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedRequest *request, const bool is_point_block, CeedVector assembled) { CeedInt num_sub; CeedOperator *suboperators; CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub)); CeedCall(CeedCompositeOperatorGetSubList(op, &suboperators)); for (CeedInt i = 0; i < num_sub; i++) { if (is_point_block) { CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(suboperators[i], assembled, request)); } else { CeedCall(CeedOperatorLinearAssembleAddDiagonal(suboperators[i], assembled, request)); } } return CEED_ERROR_SUCCESS; } /** @brief Build nonzero pattern for non-composite CeedOperator`. Users should generally use @ref CeedOperatorLinearAssembleSymbolic(). @param[in] op `CeedOperator` to assemble nonzero pattern @param[in] offset Offset for number of entries @param[out] rows Row number for each entry @param[out] cols Column number for each entry @return An error code: 0 - success, otherwise - failure @ref Developer **/ static int CeedSingleOperatorAssembleSymbolic(CeedOperator op, CeedInt offset, CeedInt *rows, CeedInt *cols) { Ceed ceed; bool is_composite; CeedSize num_nodes_in, num_nodes_out, local_num_entries, count = 0; CeedInt num_elem_in, elem_size_in, num_comp_in, layout_er_in[3]; CeedInt num_elem_out, elem_size_out, num_comp_out, layout_er_out[3]; CeedScalar *array; const CeedScalar *elem_dof_a_in, *elem_dof_a_out; CeedVector index_vec_in, index_vec_out, elem_dof_in, elem_dof_out; CeedElemRestriction elem_rstr_in, elem_rstr_out, index_elem_rstr_in, index_elem_rstr_out; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); CeedCall(CeedOperatorGetActiveVectorLengths(op, &num_nodes_in, &num_nodes_out)); CeedCall(CeedOperatorGetActiveElemRestrictions(op, &elem_rstr_in, &elem_rstr_out)); CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_in, &num_elem_in)); CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_in, &elem_size_in)); CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_in, &num_comp_in)); CeedCall(CeedElemRestrictionGetELayout(elem_rstr_in, layout_er_in)); // Determine elem_dof relation for input CeedCall(CeedVectorCreate(ceed, num_nodes_in, &index_vec_in)); CeedCall(CeedVectorGetArrayWrite(index_vec_in, CEED_MEM_HOST, &array)); for (CeedInt i = 0; i < num_nodes_in; i++) array[i] = i; CeedCall(CeedVectorRestoreArray(index_vec_in, &array)); CeedCall(CeedVectorCreate(ceed, num_elem_in * elem_size_in * num_comp_in, &elem_dof_in)); CeedCall(CeedVectorSetValue(elem_dof_in, 0.0)); CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr_in, &index_elem_rstr_in)); CeedCall(CeedElemRestrictionApply(index_elem_rstr_in, CEED_NOTRANSPOSE, index_vec_in, elem_dof_in, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedVectorGetArrayRead(elem_dof_in, CEED_MEM_HOST, &elem_dof_a_in)); CeedCall(CeedVectorDestroy(&index_vec_in)); CeedCall(CeedElemRestrictionDestroy(&index_elem_rstr_in)); if (elem_rstr_in != elem_rstr_out) { CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_out, &num_elem_out)); CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, "Active input and output operator restrictions must have the same number of elements"); CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_out, &elem_size_out)); CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_out, &num_comp_out)); CeedCall(CeedElemRestrictionGetELayout(elem_rstr_out, layout_er_out)); // Determine elem_dof relation for output CeedCall(CeedVectorCreate(ceed, num_nodes_out, &index_vec_out)); CeedCall(CeedVectorGetArrayWrite(index_vec_out, CEED_MEM_HOST, &array)); for (CeedInt i = 0; i < num_nodes_out; i++) array[i] = i; CeedCall(CeedVectorRestoreArray(index_vec_out, &array)); CeedCall(CeedVectorCreate(ceed, num_elem_out * elem_size_out * num_comp_out, &elem_dof_out)); CeedCall(CeedVectorSetValue(elem_dof_out, 0.0)); CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr_out, &index_elem_rstr_out)); CeedCall(CeedElemRestrictionApply(index_elem_rstr_out, CEED_NOTRANSPOSE, index_vec_out, elem_dof_out, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedVectorGetArrayRead(elem_dof_out, CEED_MEM_HOST, &elem_dof_a_out)); CeedCall(CeedVectorDestroy(&index_vec_out)); CeedCall(CeedElemRestrictionDestroy(&index_elem_rstr_out)); } else { num_elem_out = num_elem_in; elem_size_out = elem_size_in; num_comp_out = num_comp_in; layout_er_out[0] = layout_er_in[0]; layout_er_out[1] = layout_er_in[1]; layout_er_out[2] = layout_er_in[2]; elem_dof_a_out = elem_dof_a_in; } local_num_entries = elem_size_out * num_comp_out * elem_size_in * num_comp_in * num_elem_in; // Determine i, j locations for element matrices for (CeedInt e = 0; e < num_elem_in; e++) { for (CeedInt comp_in = 0; comp_in < num_comp_in; comp_in++) { for (CeedInt comp_out = 0; comp_out < num_comp_out; comp_out++) { for (CeedInt i = 0; i < elem_size_out; i++) { for (CeedInt j = 0; j < elem_size_in; j++) { const CeedInt elem_dof_index_row = i * layout_er_out[0] + comp_out * layout_er_out[1] + e * layout_er_out[2]; const CeedInt elem_dof_index_col = j * layout_er_in[0] + comp_in * layout_er_in[1] + e * layout_er_in[2]; const CeedInt row = elem_dof_a_out[elem_dof_index_row]; const CeedInt col = elem_dof_a_in[elem_dof_index_col]; rows[offset + count] = row; cols[offset + count] = col; count++; } } } } } CeedCheck(count == local_num_entries, ceed, CEED_ERROR_MAJOR, "Error computing assembled entries"); CeedCall(CeedVectorRestoreArrayRead(elem_dof_in, &elem_dof_a_in)); CeedCall(CeedVectorDestroy(&elem_dof_in)); if (elem_rstr_in != elem_rstr_out) { CeedCall(CeedVectorRestoreArrayRead(elem_dof_out, &elem_dof_a_out)); CeedCall(CeedVectorDestroy(&elem_dof_out)); } return CEED_ERROR_SUCCESS; } /** @brief Assemble nonzero entries for non-composite `CeedOperator`. Users should generally use @ref CeedOperatorLinearAssemble(). @param[in] op `CeedOperator` to assemble @param[in] offset Offset for number of entries @param[out] values Values to assemble into matrix @return An error code: 0 - success, otherwise - failure @ref Developer **/ static int CeedSingleOperatorAssemble(CeedOperator op, CeedInt offset, CeedVector values) { Ceed ceed; bool is_composite; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); // Early exit for empty operator { CeedInt num_elem = 0; CeedCall(CeedOperatorGetNumElements(op, &num_elem)); if (num_elem == 0) return CEED_ERROR_SUCCESS; } if (op->LinearAssembleSingle) { // Backend version CeedCall(op->LinearAssembleSingle(op, offset, values)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedSingleOperatorAssemble(op_fallback, offset, values)); return CEED_ERROR_SUCCESS; } } // Assemble QFunction CeedInt layout_qf[3]; const CeedScalar *assembled_qf_array; CeedVector assembled_qf = NULL; CeedElemRestriction assembled_elem_rstr = NULL; CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled_qf, &assembled_elem_rstr, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedElemRestrictionGetELayout(assembled_elem_rstr, layout_qf)); CeedCall(CeedElemRestrictionDestroy(&assembled_elem_rstr)); CeedCall(CeedVectorGetArrayRead(assembled_qf, CEED_MEM_HOST, &assembled_qf_array)); // Get assembly data CeedInt num_elem_in, elem_size_in, num_comp_in, num_qpts_in; CeedInt num_elem_out, elem_size_out, num_comp_out, num_qpts_out; CeedSize local_num_entries, count = 0; const CeedEvalMode **eval_modes_in, **eval_modes_out; CeedInt num_active_bases_in, *num_eval_modes_in, num_active_bases_out, *num_eval_modes_out; CeedBasis *active_bases_in, *active_bases_out, basis_in, basis_out; const CeedScalar **B_mats_in, **B_mats_out, *B_mat_in, *B_mat_out; CeedElemRestriction elem_rstr_in, elem_rstr_out; CeedRestrictionType elem_rstr_type_in, elem_rstr_type_out; const bool *elem_rstr_orients_in = NULL, *elem_rstr_orients_out = NULL; const CeedInt8 *elem_rstr_curl_orients_in = NULL, *elem_rstr_curl_orients_out = NULL; CeedOperatorAssemblyData data; CeedCall(CeedOperatorGetOperatorAssemblyData(op, &data)); CeedCall(CeedOperatorAssemblyDataGetEvalModes(data, &num_active_bases_in, &num_eval_modes_in, &eval_modes_in, NULL, &num_active_bases_out, &num_eval_modes_out, &eval_modes_out, NULL, NULL)); CeedCheck(num_active_bases_in == num_active_bases_out && num_active_bases_in == 1, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator with multiple active bases"); CeedCheck(num_eval_modes_in[0] > 0 && num_eval_modes_out[0] > 0, ceed, CEED_ERROR_UNSUPPORTED, "Cannot assemble operator without inputs/outputs"); CeedCall(CeedOperatorAssemblyDataGetBases(data, NULL, &active_bases_in, &B_mats_in, NULL, &active_bases_out, &B_mats_out)); CeedCall(CeedOperatorGetActiveElemRestrictions(op, &elem_rstr_in, &elem_rstr_out)); basis_in = active_bases_in[0]; basis_out = active_bases_out[0]; B_mat_in = B_mats_in[0]; B_mat_out = B_mats_out[0]; CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_in, &num_elem_in)); CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_in, &elem_size_in)); CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_in, &num_comp_in)); if (basis_in == CEED_BASIS_NONE) num_qpts_in = elem_size_in; else CeedCall(CeedBasisGetNumQuadraturePoints(basis_in, &num_qpts_in)); CeedCall(CeedElemRestrictionGetType(elem_rstr_in, &elem_rstr_type_in)); if (elem_rstr_type_in == CEED_RESTRICTION_ORIENTED) { CeedCall(CeedElemRestrictionGetOrientations(elem_rstr_in, CEED_MEM_HOST, &elem_rstr_orients_in)); } else if (elem_rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { CeedCall(CeedElemRestrictionGetCurlOrientations(elem_rstr_in, CEED_MEM_HOST, &elem_rstr_curl_orients_in)); } if (elem_rstr_in != elem_rstr_out) { CeedCall(CeedElemRestrictionGetNumElements(elem_rstr_out, &num_elem_out)); CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, "Active input and output operator restrictions must have the same number of elements"); CeedCall(CeedElemRestrictionGetElementSize(elem_rstr_out, &elem_size_out)); CeedCall(CeedElemRestrictionGetNumComponents(elem_rstr_out, &num_comp_out)); if (basis_out == CEED_BASIS_NONE) num_qpts_out = elem_size_out; else CeedCall(CeedBasisGetNumQuadraturePoints(basis_out, &num_qpts_out)); CeedCheck(num_qpts_in == num_qpts_out, ceed, CEED_ERROR_UNSUPPORTED, "Active input and output bases must have the same number of quadrature points"); CeedCall(CeedElemRestrictionGetType(elem_rstr_out, &elem_rstr_type_out)); if (elem_rstr_type_out == CEED_RESTRICTION_ORIENTED) { CeedCall(CeedElemRestrictionGetOrientations(elem_rstr_out, CEED_MEM_HOST, &elem_rstr_orients_out)); } else if (elem_rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { CeedCall(CeedElemRestrictionGetCurlOrientations(elem_rstr_out, CEED_MEM_HOST, &elem_rstr_curl_orients_out)); } } else { num_elem_out = num_elem_in; elem_size_out = elem_size_in; num_comp_out = num_comp_in; num_qpts_out = num_qpts_in; elem_rstr_orients_out = elem_rstr_orients_in; elem_rstr_curl_orients_out = elem_rstr_curl_orients_in; } local_num_entries = elem_size_out * num_comp_out * elem_size_in * num_comp_in * num_elem_in; // Loop over elements and put in data structure // We store B_mat_in, B_mat_out, BTD, elem_mat in row-major order CeedTensorContract contract; CeedScalar *vals, *BTD_mat = NULL, *elem_mat = NULL, *elem_mat_b = NULL; CeedCall(CeedBasisGetTensorContract(basis_in, &contract)); CeedCall(CeedCalloc(elem_size_out * num_qpts_in * num_eval_modes_in[0], &BTD_mat)); CeedCall(CeedCalloc(elem_size_out * elem_size_in, &elem_mat)); if (elem_rstr_curl_orients_in || elem_rstr_curl_orients_out) CeedCall(CeedCalloc(elem_size_out * elem_size_in, &elem_mat_b)); CeedCall(CeedVectorGetArray(values, CEED_MEM_HOST, &vals)); for (CeedSize e = 0; e < num_elem_in; e++) { for (CeedInt comp_in = 0; comp_in < num_comp_in; comp_in++) { for (CeedInt comp_out = 0; comp_out < num_comp_out; comp_out++) { // Compute B^T*D for (CeedSize n = 0; n < elem_size_out; n++) { for (CeedSize q = 0; q < num_qpts_in; q++) { for (CeedInt e_in = 0; e_in < num_eval_modes_in[0]; e_in++) { const CeedSize btd_index = n * (num_qpts_in * num_eval_modes_in[0]) + q * num_eval_modes_in[0] + e_in; CeedScalar sum = 0.0; for (CeedInt e_out = 0; e_out < num_eval_modes_out[0]; e_out++) { const CeedSize b_out_index = (q * num_eval_modes_out[0] + e_out) * elem_size_out + n; const CeedSize eval_mode_index = ((e_in * num_comp_in + comp_in) * num_eval_modes_out[0] + e_out) * num_comp_out + comp_out; const CeedSize qf_index = q * layout_qf[0] + eval_mode_index * layout_qf[1] + e * layout_qf[2]; sum += B_mat_out[b_out_index] * assembled_qf_array[qf_index]; } BTD_mat[btd_index] = sum; } } } // Form element matrix itself (for each block component) if (contract) { CeedCall(CeedTensorContractApply(contract, 1, num_qpts_in * num_eval_modes_in[0], elem_size_in, elem_size_out, BTD_mat, CEED_NOTRANSPOSE, false, B_mat_in, elem_mat)); } else { CeedCall(CeedMatrixMatrixMultiply(ceed, BTD_mat, B_mat_in, elem_mat, elem_size_out, elem_size_in, num_qpts_in * num_eval_modes_in[0])); } // Transform the element matrix if required if (elem_rstr_orients_out) { const bool *elem_orients = &elem_rstr_orients_out[e * elem_size_out]; for (CeedInt i = 0; i < elem_size_out; i++) { const double orient = elem_orients[i] ? -1.0 : 1.0; for (CeedInt j = 0; j < elem_size_in; j++) { elem_mat[i * elem_size_in + j] *= orient; } } } else if (elem_rstr_curl_orients_out) { const CeedInt8 *elem_curl_orients = &elem_rstr_curl_orients_out[e * 3 * elem_size_out]; // T^T*(B^T*D*B) memcpy(elem_mat_b, elem_mat, elem_size_out * elem_size_in * sizeof(CeedScalar)); for (CeedInt i = 0; i < elem_size_out; i++) { for (CeedInt j = 0; j < elem_size_in; j++) { elem_mat[i * elem_size_in + j] = elem_mat_b[i * elem_size_in + j] * elem_curl_orients[3 * i + 1] + (i > 0 ? elem_mat_b[(i - 1) * elem_size_in + j] * elem_curl_orients[3 * i - 1] : 0.0) + (i < elem_size_out - 1 ? elem_mat_b[(i + 1) * elem_size_in + j] * elem_curl_orients[3 * i + 3] : 0.0); } } } if (elem_rstr_orients_in) { const bool *elem_orients = &elem_rstr_orients_in[e * elem_size_in]; for (CeedInt i = 0; i < elem_size_out; i++) { for (CeedInt j = 0; j < elem_size_in; j++) { elem_mat[i * elem_size_in + j] *= elem_orients[j] ? -1.0 : 1.0; } } } else if (elem_rstr_curl_orients_in) { const CeedInt8 *elem_curl_orients = &elem_rstr_curl_orients_in[e * 3 * elem_size_in]; // (B^T*D*B)*T memcpy(elem_mat_b, elem_mat, elem_size_out * elem_size_in * sizeof(CeedScalar)); for (CeedInt i = 0; i < elem_size_out; i++) { for (CeedInt j = 0; j < elem_size_in; j++) { elem_mat[i * elem_size_in + j] = elem_mat_b[i * elem_size_in + j] * elem_curl_orients[3 * j + 1] + (j > 0 ? elem_mat_b[i * elem_size_in + j - 1] * elem_curl_orients[3 * j - 1] : 0.0) + (j < elem_size_in - 1 ? elem_mat_b[i * elem_size_in + j + 1] * elem_curl_orients[3 * j + 3] : 0.0); } } } // Put element matrix in coordinate data structure for (CeedInt i = 0; i < elem_size_out; i++) { for (CeedInt j = 0; j < elem_size_in; j++) { vals[offset + count] = elem_mat[i * elem_size_in + j]; count++; } } } } } CeedCheck(count == local_num_entries, ceed, CEED_ERROR_MAJOR, "Error computing entries"); CeedCall(CeedVectorRestoreArray(values, &vals)); // Cleanup CeedCall(CeedFree(&BTD_mat)); CeedCall(CeedFree(&elem_mat)); CeedCall(CeedFree(&elem_mat_b)); if (elem_rstr_type_in == CEED_RESTRICTION_ORIENTED) { CeedCall(CeedElemRestrictionRestoreOrientations(elem_rstr_in, &elem_rstr_orients_in)); } else if (elem_rstr_type_in == CEED_RESTRICTION_CURL_ORIENTED) { CeedCall(CeedElemRestrictionRestoreCurlOrientations(elem_rstr_in, &elem_rstr_curl_orients_in)); } if (elem_rstr_in != elem_rstr_out) { if (elem_rstr_type_out == CEED_RESTRICTION_ORIENTED) { CeedCall(CeedElemRestrictionRestoreOrientations(elem_rstr_out, &elem_rstr_orients_out)); } else if (elem_rstr_type_out == CEED_RESTRICTION_CURL_ORIENTED) { CeedCall(CeedElemRestrictionRestoreCurlOrientations(elem_rstr_out, &elem_rstr_curl_orients_out)); } } CeedCall(CeedVectorRestoreArrayRead(assembled_qf, &assembled_qf_array)); CeedCall(CeedVectorDestroy(&assembled_qf)); return CEED_ERROR_SUCCESS; } /** @brief Count number of entries for assembled `CeedOperator` @param[in] op `CeedOperator` to assemble @param[out] num_entries Number of entries in assembled representation @return An error code: 0 - success, otherwise - failure @ref Utility **/ static int CeedSingleOperatorAssemblyCountEntries(CeedOperator op, CeedSize *num_entries) { bool is_composite; CeedInt num_elem_in, elem_size_in, num_comp_in, num_elem_out, elem_size_out, num_comp_out; Ceed ceed; CeedElemRestriction rstr_in, rstr_out; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Composite operator not supported"); CeedCall(CeedOperatorGetActiveElemRestrictions(op, &rstr_in, &rstr_out)); CeedCall(CeedElemRestrictionGetNumElements(rstr_in, &num_elem_in)); CeedCall(CeedElemRestrictionGetElementSize(rstr_in, &elem_size_in)); CeedCall(CeedElemRestrictionGetNumComponents(rstr_in, &num_comp_in)); if (rstr_in != rstr_out) { CeedCall(CeedElemRestrictionGetNumElements(rstr_out, &num_elem_out)); CeedCheck(num_elem_in == num_elem_out, ceed, CEED_ERROR_UNSUPPORTED, "Active input and output operator restrictions must have the same number of elements"); CeedCall(CeedElemRestrictionGetElementSize(rstr_out, &elem_size_out)); CeedCall(CeedElemRestrictionGetNumComponents(rstr_out, &num_comp_out)); } else { num_elem_out = num_elem_in; elem_size_out = elem_size_in; num_comp_out = num_comp_in; } *num_entries = (CeedSize)elem_size_in * num_comp_in * elem_size_out * num_comp_out * num_elem_in; return CEED_ERROR_SUCCESS; } /** @brief Common code for creating a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` @param[in] op_fine Fine grid `CeedOperator` @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` @param[in] rstr_coarse Coarse grid `CeedElemRestriction` @param[in] basis_coarse Coarse grid active vector `CeedBasis` @param[in] basis_c_to_f `CeedBasis` for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction operators @param[out] op_coarse Coarse grid `CeedOperator` @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` @return An error code: 0 - success, otherwise - failure @ref Developer **/ static int CeedSingleOperatorMultigridLevel(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, CeedBasis basis_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { bool is_composite; Ceed ceed; CeedInt num_comp, num_input_fields, num_output_fields; CeedVector mult_vec = NULL; CeedElemRestriction rstr_p_mult_fine = NULL, rstr_fine = NULL; CeedOperatorField *input_fields, *output_fields; CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); // Check for composite operator CeedCall(CeedOperatorIsComposite(op_fine, &is_composite)); CeedCheck(!is_composite, ceed, CEED_ERROR_UNSUPPORTED, "Automatic multigrid setup for composite operators not supported"); // Coarse Grid CeedCall(CeedOperatorCreate(ceed, op_fine->qf, op_fine->dqf, op_fine->dqfT, op_coarse)); CeedCall(CeedOperatorGetFields(op_fine, &num_input_fields, &input_fields, &num_output_fields, &output_fields)); // -- Clone input fields for (CeedInt i = 0; i < num_input_fields; i++) { const char *field_name; CeedVector vec; CeedElemRestriction rstr; CeedBasis basis; CeedCall(CeedOperatorFieldGetName(input_fields[i], &field_name)); CeedCall(CeedOperatorFieldGetVector(input_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { rstr = rstr_coarse; basis = basis_coarse; CeedCall(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr_fine)); } else { CeedCall(CeedOperatorFieldGetElemRestriction(input_fields[i], &rstr)); CeedCall(CeedOperatorFieldGetBasis(input_fields[i], &basis)); } CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec)); } // -- Clone output fields for (CeedInt i = 0; i < num_output_fields; i++) { const char *field_name; CeedVector vec; CeedElemRestriction rstr; CeedBasis basis; CeedCall(CeedOperatorFieldGetName(output_fields[i], &field_name)); CeedCall(CeedOperatorFieldGetVector(output_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { rstr = rstr_coarse; basis = basis_coarse; CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr_fine)); } else { CeedCall(CeedOperatorFieldGetElemRestriction(output_fields[i], &rstr)); CeedCall(CeedOperatorFieldGetBasis(output_fields[i], &basis)); } CeedCall(CeedOperatorSetField(*op_coarse, field_name, rstr, basis, vec)); } // -- Clone QFunctionAssemblyData CeedCall(CeedQFunctionAssemblyDataReferenceCopy(op_fine->qf_assembled, &(*op_coarse)->qf_assembled)); // Multiplicity vector if (op_restrict || op_prolong) { CeedVector mult_e_vec; CeedRestrictionType rstr_type; CeedCall(CeedElemRestrictionGetType(rstr_fine, &rstr_type)); CeedCheck(rstr_type != CEED_RESTRICTION_CURL_ORIENTED, ceed, CEED_ERROR_UNSUPPORTED, "Element restrictions created with CeedElemRestrictionCreateCurlOriented are not supported"); CeedCheck(p_mult_fine, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires fine grid multiplicity vector"); CeedCall(CeedElemRestrictionCreateUnsignedCopy(rstr_fine, &rstr_p_mult_fine)); CeedCall(CeedElemRestrictionCreateVector(rstr_fine, &mult_vec, &mult_e_vec)); CeedCall(CeedVectorSetValue(mult_e_vec, 0.0)); CeedCall(CeedElemRestrictionApply(rstr_p_mult_fine, CEED_NOTRANSPOSE, p_mult_fine, mult_e_vec, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedVectorSetValue(mult_vec, 0.0)); CeedCall(CeedElemRestrictionApply(rstr_p_mult_fine, CEED_TRANSPOSE, mult_e_vec, mult_vec, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedVectorDestroy(&mult_e_vec)); CeedCall(CeedVectorReciprocal(mult_vec)); } // Clone name bool has_name = op_fine->name; size_t name_len = op_fine->name ? strlen(op_fine->name) : 0; CeedCall(CeedOperatorSetName(*op_coarse, op_fine->name)); // Check that coarse to fine basis is provided if prolong/restrict operators are requested CeedCheck(basis_c_to_f || (!op_restrict && !op_prolong), ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires coarse-to-fine basis"); // Restriction/Prolongation Operators CeedCall(CeedBasisGetNumComponents(basis_coarse, &num_comp)); // Restriction if (op_restrict) { CeedInt *num_comp_r_data; CeedQFunctionContext ctx_r; CeedQFunction qf_restrict; CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_restrict)); CeedCall(CeedCalloc(1, &num_comp_r_data)); num_comp_r_data[0] = num_comp; CeedCall(CeedQFunctionContextCreate(ceed, &ctx_r)); CeedCall(CeedQFunctionContextSetData(ctx_r, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_r_data), num_comp_r_data)); CeedCall(CeedQFunctionSetContext(qf_restrict, ctx_r)); CeedCall(CeedQFunctionContextDestroy(&ctx_r)); CeedCall(CeedQFunctionAddInput(qf_restrict, "input", num_comp, CEED_EVAL_NONE)); CeedCall(CeedQFunctionAddInput(qf_restrict, "scale", num_comp, CEED_EVAL_NONE)); CeedCall(CeedQFunctionAddOutput(qf_restrict, "output", num_comp, CEED_EVAL_INTERP)); CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_restrict, num_comp)); CeedCall(CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_restrict)); CeedCall(CeedOperatorSetField(*op_restrict, "input", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); CeedCall(CeedOperatorSetField(*op_restrict, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec)); CeedCall(CeedOperatorSetField(*op_restrict, "output", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); // Set name char *restriction_name; CeedCall(CeedCalloc(17 + name_len, &restriction_name)); sprintf(restriction_name, "restriction%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); CeedCall(CeedOperatorSetName(*op_restrict, restriction_name)); CeedCall(CeedFree(&restriction_name)); // Check CeedCall(CeedOperatorCheckReady(*op_restrict)); // Cleanup CeedCall(CeedQFunctionDestroy(&qf_restrict)); } // Prolongation if (op_prolong) { CeedInt *num_comp_p_data; CeedQFunctionContext ctx_p; CeedQFunction qf_prolong; CeedCall(CeedQFunctionCreateInteriorByName(ceed, "Scale", &qf_prolong)); CeedCall(CeedCalloc(1, &num_comp_p_data)); num_comp_p_data[0] = num_comp; CeedCall(CeedQFunctionContextCreate(ceed, &ctx_p)); CeedCall(CeedQFunctionContextSetData(ctx_p, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_p_data), num_comp_p_data)); CeedCall(CeedQFunctionSetContext(qf_prolong, ctx_p)); CeedCall(CeedQFunctionContextDestroy(&ctx_p)); CeedCall(CeedQFunctionAddInput(qf_prolong, "input", num_comp, CEED_EVAL_INTERP)); CeedCall(CeedQFunctionAddInput(qf_prolong, "scale", num_comp, CEED_EVAL_NONE)); CeedCall(CeedQFunctionAddOutput(qf_prolong, "output", num_comp, CEED_EVAL_NONE)); CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_prolong, num_comp)); CeedCall(CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, op_prolong)); CeedCall(CeedOperatorSetField(*op_prolong, "input", rstr_coarse, basis_c_to_f, CEED_VECTOR_ACTIVE)); CeedCall(CeedOperatorSetField(*op_prolong, "scale", rstr_p_mult_fine, CEED_BASIS_NONE, mult_vec)); CeedCall(CeedOperatorSetField(*op_prolong, "output", rstr_fine, CEED_BASIS_NONE, CEED_VECTOR_ACTIVE)); // Set name char *prolongation_name; CeedCall(CeedCalloc(18 + name_len, &prolongation_name)); sprintf(prolongation_name, "prolongation%s%s", has_name ? " for " : "", has_name ? op_fine->name : ""); CeedCall(CeedOperatorSetName(*op_prolong, prolongation_name)); CeedCall(CeedFree(&prolongation_name)); // Check CeedCall(CeedOperatorCheckReady(*op_prolong)); // Cleanup CeedCall(CeedQFunctionDestroy(&qf_prolong)); } // Check CeedCall(CeedOperatorCheckReady(*op_coarse)); // Cleanup CeedCall(CeedVectorDestroy(&mult_vec)); CeedCall(CeedElemRestrictionDestroy(&rstr_p_mult_fine)); CeedCall(CeedBasisDestroy(&basis_c_to_f)); return CEED_ERROR_SUCCESS; } /** @brief Build 1D mass matrix and Laplacian with perturbation @param[in] interp_1d Interpolation matrix in one dimension @param[in] grad_1d Gradient matrix in one dimension @param[in] q_weight_1d Quadrature weights in one dimension @param[in] P_1d Number of basis nodes in one dimension @param[in] Q_1d Number of quadrature points in one dimension @param[in] dim Dimension of basis @param[out] mass Assembled mass matrix in one dimension @param[out] laplace Assembled perturbed Laplacian in one dimension @return An error code: 0 - success, otherwise - failure @ref Developer **/ CeedPragmaOptimizeOff static int CeedBuildMassLaplace(const CeedScalar *interp_1d, const CeedScalar *grad_1d, const CeedScalar *q_weight_1d, CeedInt P_1d, CeedInt Q_1d, CeedInt dim, CeedScalar *mass, CeedScalar *laplace) { for (CeedInt i = 0; i < P_1d; i++) { for (CeedInt j = 0; j < P_1d; j++) { CeedScalar sum = 0.0; 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]; mass[i + j * P_1d] = sum; } } // -- Laplacian for (CeedInt i = 0; i < P_1d; i++) { for (CeedInt j = 0; j < P_1d; j++) { CeedScalar sum = 0.0; 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]; laplace[i + j * P_1d] = sum; } } CeedScalar perturbation = dim > 2 ? 1e-6 : 1e-4; for (CeedInt i = 0; i < P_1d; i++) laplace[i + P_1d * i] += perturbation; return CEED_ERROR_SUCCESS; } CeedPragmaOptimizeOn /// @} /// ---------------------------------------------------------------------------- /// CeedOperator Backend API /// ---------------------------------------------------------------------------- /// @addtogroup CeedOperatorBackend /// @{ /** @brief Select correct basis matrix pointer based on @ref CeedEvalMode @param[in] basis `CeedBasis` from which to get the basis matrix @param[in] eval_mode Current basis evaluation mode @param[in] identity Pointer to identity matrix @param[out] basis_ptr `CeedBasis` pointer to set @ref Backend **/ int CeedOperatorGetBasisPointer(CeedBasis basis, CeedEvalMode eval_mode, const CeedScalar *identity, const CeedScalar **basis_ptr) { switch (eval_mode) { case CEED_EVAL_NONE: *basis_ptr = identity; break; case CEED_EVAL_INTERP: CeedCall(CeedBasisGetInterp(basis, basis_ptr)); break; case CEED_EVAL_GRAD: CeedCall(CeedBasisGetGrad(basis, basis_ptr)); break; case CEED_EVAL_DIV: CeedCall(CeedBasisGetDiv(basis, basis_ptr)); break; case CEED_EVAL_CURL: CeedCall(CeedBasisGetCurl(basis, basis_ptr)); break; case CEED_EVAL_WEIGHT: break; // Caught by QF Assembly } assert(*basis_ptr != NULL); return CEED_ERROR_SUCCESS; } /** @brief Create point block restriction for active `CeedOperatorField` @param[in] rstr Original `CeedElemRestriction` for active field @param[out] point_block_rstr Address of the variable where the newly created `CeedElemRestriction` will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorCreateActivePointBlockRestriction(CeedElemRestriction rstr, CeedElemRestriction *point_block_rstr) { Ceed ceed; CeedInt num_elem, num_comp, shift, elem_size, comp_stride, *point_block_offsets; CeedSize l_size; const CeedInt *offsets; CeedCall(CeedElemRestrictionGetCeed(rstr, &ceed)); CeedCall(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); // Expand offsets CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); CeedCall(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); CeedCall(CeedElemRestrictionGetElementSize(rstr, &elem_size)); CeedCall(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); shift = num_comp; if (comp_stride != 1) shift *= num_comp; CeedCall(CeedCalloc(num_elem * elem_size, &point_block_offsets)); for (CeedInt i = 0; i < num_elem * elem_size; i++) { point_block_offsets[i] = offsets[i] * shift; } // Create new restriction CeedCall(CeedElemRestrictionCreate(ceed, num_elem, elem_size, num_comp * num_comp, 1, l_size * num_comp, CEED_MEM_HOST, CEED_OWN_POINTER, point_block_offsets, point_block_rstr)); // Cleanup CeedCall(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); return CEED_ERROR_SUCCESS; } /** @brief Create object holding `CeedQFunction` assembly data for `CeedOperator` @param[in] ceed `Ceed` object used to create the `CeedQFunctionAssemblyData` @param[out] data Address of the variable where the newly created `CeedQFunctionAssemblyData` will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataCreate(Ceed ceed, CeedQFunctionAssemblyData *data) { CeedCall(CeedCalloc(1, data)); (*data)->ref_count = 1; (*data)->ceed = ceed; CeedCall(CeedReference(ceed)); return CEED_ERROR_SUCCESS; } /** @brief Increment the reference counter for a `CeedQFunctionAssemblyData` @param[in,out] data `CeedQFunctionAssemblyData` to increment the reference counter @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataReference(CeedQFunctionAssemblyData data) { data->ref_count++; return CEED_ERROR_SUCCESS; } /** @brief Set re-use of `CeedQFunctionAssemblyData` @param[in,out] data `CeedQFunctionAssemblyData` to mark for reuse @param[in] reuse_data Boolean flag indicating data re-use @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataSetReuse(CeedQFunctionAssemblyData data, bool reuse_data) { data->reuse_data = reuse_data; data->needs_data_update = true; return CEED_ERROR_SUCCESS; } /** @brief Mark `CeedQFunctionAssemblyData` as stale @param[in,out] data `CeedQFunctionAssemblyData` to mark as stale @param[in] needs_data_update Boolean flag indicating if update is needed or completed @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataSetUpdateNeeded(CeedQFunctionAssemblyData data, bool needs_data_update) { data->needs_data_update = needs_data_update; return CEED_ERROR_SUCCESS; } /** @brief Determine if `CeedQFunctionAssemblyData` needs update @param[in] data `CeedQFunctionAssemblyData` to mark as stale @param[out] is_update_needed Boolean flag indicating if re-assembly is required @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataIsUpdateNeeded(CeedQFunctionAssemblyData data, bool *is_update_needed) { *is_update_needed = !data->reuse_data || data->needs_data_update; return CEED_ERROR_SUCCESS; } /** @brief Copy the pointer to a `CeedQFunctionAssemblyData`. Both pointers should be destroyed with @ref CeedQFunctionAssemblyDataDestroy(). 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`. This `CeedQFunctionAssemblyData` will be destroyed if ` *data_copy` is the only reference to this `CeedQFunctionAssemblyData`. @param[in] data `CeedQFunctionAssemblyData` to copy reference to @param[in,out] data_copy Variable to store copied reference @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataReferenceCopy(CeedQFunctionAssemblyData data, CeedQFunctionAssemblyData *data_copy) { CeedCall(CeedQFunctionAssemblyDataReference(data)); CeedCall(CeedQFunctionAssemblyDataDestroy(data_copy)); *data_copy = data; return CEED_ERROR_SUCCESS; } /** @brief Get setup status for internal objects for `CeedQFunctionAssemblyData` @param[in] data `CeedQFunctionAssemblyData` to retrieve status @param[out] is_setup Boolean flag for setup status @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataIsSetup(CeedQFunctionAssemblyData data, bool *is_setup) { *is_setup = data->is_setup; return CEED_ERROR_SUCCESS; } /** @brief Set internal objects for `CeedQFunctionAssemblyData` @param[in,out] data `CeedQFunctionAssemblyData` to set objects @param[in] vec `CeedVector` to store assembled `CeedQFunction` at quadrature points @param[in] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataSetObjects(CeedQFunctionAssemblyData data, CeedVector vec, CeedElemRestriction rstr) { CeedCall(CeedVectorReferenceCopy(vec, &data->vec)); CeedCall(CeedElemRestrictionReferenceCopy(rstr, &data->rstr)); data->is_setup = true; return CEED_ERROR_SUCCESS; } /** @brief Get internal objects for `CeedQFunctionAssemblyData` @param[in,out] data `CeedQFunctionAssemblyData` to set objects @param[out] vec `CeedVector` to store assembled `CeedQFunction` at quadrature points @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataGetObjects(CeedQFunctionAssemblyData data, CeedVector *vec, CeedElemRestriction *rstr) { CeedCheck(data->is_setup, data->ceed, CEED_ERROR_INCOMPLETE, "Internal objects not set; must call CeedQFunctionAssemblyDataSetObjects first."); CeedCall(CeedVectorReferenceCopy(data->vec, vec)); CeedCall(CeedElemRestrictionReferenceCopy(data->rstr, rstr)); return CEED_ERROR_SUCCESS; } /** @brief Destroy `CeedQFunctionAssemblyData` @param[in,out] data `CeedQFunctionAssemblyData` to destroy @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedQFunctionAssemblyDataDestroy(CeedQFunctionAssemblyData *data) { if (!*data || --(*data)->ref_count > 0) { *data = NULL; return CEED_ERROR_SUCCESS; } CeedCall(CeedDestroy(&(*data)->ceed)); CeedCall(CeedVectorDestroy(&(*data)->vec)); CeedCall(CeedElemRestrictionDestroy(&(*data)->rstr)); CeedCall(CeedFree(data)); return CEED_ERROR_SUCCESS; } /** @brief Get `CeedOperatorAssemblyData` @param[in] op `CeedOperator` to assemble @param[out] data `CeedQFunctionAssemblyData` @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorGetOperatorAssemblyData(CeedOperator op, CeedOperatorAssemblyData *data) { if (!op->op_assembled) { CeedOperatorAssemblyData data; CeedCall(CeedOperatorAssemblyDataCreate(op->ceed, op, &data)); op->op_assembled = data; } *data = op->op_assembled; return CEED_ERROR_SUCCESS; } /** @brief Create object holding `CeedOperator` assembly data. The `CeedOperatorAssemblyData` holds an array with references to every active `CeedBasis` used in the `CeedOperator`. An array with references to the corresponding active `CeedElemRestriction` is also stored. For each active `CeedBasis, the `CeedOperatorAssemblyData` holds an array of all input and output @ref CeedEvalMode for this `CeedBasis`. The `CeedOperatorAssemblyData` holds an array of offsets for indexing into the assembled `CeedQFunction` arrays to the row representing each @ref CeedEvalMode. The number of input columns across all active bases for the assembled `CeedQFunction` is also stored. Lastly, the `CeedOperatorAssembly` data holds assembled matrices representing the full action of the `CeedBasis` for all @ref CeedEvalMode. @param[in] ceed `Ceed` object used to create the `CeedOperatorAssemblyData` @param[in] op `CeedOperator` to be assembled @param[out] data Address of the variable where the newly created `CeedOperatorAssemblyData` will be stored @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorAssemblyDataCreate(Ceed ceed, CeedOperator op, CeedOperatorAssemblyData *data) { CeedInt num_active_bases_in = 0, num_active_bases_out = 0, offset = 0; CeedInt num_input_fields, *num_eval_modes_in = NULL, num_output_fields, *num_eval_modes_out = NULL; CeedSize **eval_mode_offsets_in = NULL, **eval_mode_offsets_out = NULL; CeedEvalMode **eval_modes_in = NULL, **eval_modes_out = NULL; CeedQFunctionField *qf_fields; CeedQFunction qf; CeedOperatorField *op_fields; bool is_composite; CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCheck(!is_composite, ceed, CEED_ERROR_INCOMPATIBLE, "Can only create CeedOperator assembly data for non-composite operators."); // Allocate CeedCall(CeedCalloc(1, data)); (*data)->ceed = ceed; CeedCall(CeedReference(ceed)); // Build OperatorAssembly data CeedCall(CeedOperatorGetQFunction(op, &qf)); // Determine active input basis CeedCall(CeedQFunctionGetFields(qf, &num_input_fields, &qf_fields, NULL, NULL)); CeedCall(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); for (CeedInt i = 0; i < num_input_fields; i++) { CeedVector vec; CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { CeedInt index = -1, num_comp, q_comp; CeedEvalMode eval_mode; CeedBasis basis_in = NULL; CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_in)); CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); CeedCall(CeedBasisGetNumComponents(basis_in, &num_comp)); CeedCall(CeedBasisGetNumQuadratureComponents(basis_in, eval_mode, &q_comp)); for (CeedInt i = 0; i < num_active_bases_in; i++) { if ((*data)->active_bases_in[i] == basis_in) index = i; } if (index == -1) { CeedElemRestriction elem_rstr_in; index = num_active_bases_in; CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_bases_in)); (*data)->active_bases_in[num_active_bases_in] = NULL; CeedCall(CeedBasisReferenceCopy(basis_in, &(*data)->active_bases_in[num_active_bases_in])); CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->active_elem_rstrs_in)); (*data)->active_elem_rstrs_in[num_active_bases_in] = NULL; CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_in)); CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_in, &(*data)->active_elem_rstrs_in[num_active_bases_in])); CeedCall(CeedRealloc(num_active_bases_in + 1, &num_eval_modes_in)); num_eval_modes_in[index] = 0; CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_modes_in)); eval_modes_in[index] = NULL; CeedCall(CeedRealloc(num_active_bases_in + 1, &eval_mode_offsets_in)); eval_mode_offsets_in[index] = NULL; CeedCall(CeedRealloc(num_active_bases_in + 1, &(*data)->assembled_bases_in)); (*data)->assembled_bases_in[index] = NULL; num_active_bases_in++; } if (eval_mode != CEED_EVAL_WEIGHT) { // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_modes_in[index])); CeedCall(CeedRealloc(num_eval_modes_in[index] + q_comp, &eval_mode_offsets_in[index])); for (CeedInt d = 0; d < q_comp; d++) { eval_modes_in[index][num_eval_modes_in[index] + d] = eval_mode; eval_mode_offsets_in[index][num_eval_modes_in[index] + d] = offset; offset += num_comp; } num_eval_modes_in[index] += q_comp; } } } // Determine active output basis CeedCall(CeedQFunctionGetFields(qf, NULL, NULL, &num_output_fields, &qf_fields)); CeedCall(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); offset = 0; for (CeedInt i = 0; i < num_output_fields; i++) { CeedVector vec; CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { CeedInt index = -1, num_comp, q_comp; CeedEvalMode eval_mode; CeedBasis basis_out = NULL; CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis_out)); CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); CeedCall(CeedBasisGetNumComponents(basis_out, &num_comp)); CeedCall(CeedBasisGetNumQuadratureComponents(basis_out, eval_mode, &q_comp)); for (CeedInt i = 0; i < num_active_bases_out; i++) { if ((*data)->active_bases_out[i] == basis_out) index = i; } if (index == -1) { CeedElemRestriction elem_rstr_out; index = num_active_bases_out; CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_bases_out)); (*data)->active_bases_out[num_active_bases_out] = NULL; CeedCall(CeedBasisReferenceCopy(basis_out, &(*data)->active_bases_out[num_active_bases_out])); CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->active_elem_rstrs_out)); (*data)->active_elem_rstrs_out[num_active_bases_out] = NULL; CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &elem_rstr_out)); CeedCall(CeedElemRestrictionReferenceCopy(elem_rstr_out, &(*data)->active_elem_rstrs_out[num_active_bases_out])); CeedCall(CeedRealloc(num_active_bases_out + 1, &num_eval_modes_out)); num_eval_modes_out[index] = 0; CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_modes_out)); eval_modes_out[index] = NULL; CeedCall(CeedRealloc(num_active_bases_out + 1, &eval_mode_offsets_out)); eval_mode_offsets_out[index] = NULL; CeedCall(CeedRealloc(num_active_bases_out + 1, &(*data)->assembled_bases_out)); (*data)->assembled_bases_out[index] = NULL; num_active_bases_out++; } if (eval_mode != CEED_EVAL_WEIGHT) { // q_comp = 1 if CEED_EVAL_NONE, CEED_EVAL_WEIGHT caught by QF Assembly CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_modes_out[index])); CeedCall(CeedRealloc(num_eval_modes_out[index] + q_comp, &eval_mode_offsets_out[index])); for (CeedInt d = 0; d < q_comp; d++) { eval_modes_out[index][num_eval_modes_out[index] + d] = eval_mode; eval_mode_offsets_out[index][num_eval_modes_out[index] + d] = offset; offset += num_comp; } num_eval_modes_out[index] += q_comp; } } } (*data)->num_active_bases_in = num_active_bases_in; (*data)->num_eval_modes_in = num_eval_modes_in; (*data)->eval_modes_in = eval_modes_in; (*data)->eval_mode_offsets_in = eval_mode_offsets_in; (*data)->num_active_bases_out = num_active_bases_out; (*data)->num_eval_modes_out = num_eval_modes_out; (*data)->eval_modes_out = eval_modes_out; (*data)->eval_mode_offsets_out = eval_mode_offsets_out; (*data)->num_output_components = offset; return CEED_ERROR_SUCCESS; } /** @brief Get `CeedOperator` @ref CeedEvalMode for assembly. Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. @param[in] data `CeedOperatorAssemblyData` @param[out] num_active_bases_in Total number of active bases for input @param[out] num_eval_modes_in Pointer to hold array of numbers of input @ref CeedEvalMode, or `NULL`. `eval_modes_in[0]` holds an array of eval modes for the first active `CeedBasis`. @param[out] eval_modes_in Pointer to hold arrays of input @ref CeedEvalMode, or `NULL` @param[out] eval_mode_offsets_in Pointer to hold arrays of input offsets at each quadrature point @param[out] num_active_bases_out Total number of active bases for output @param[out] num_eval_modes_out Pointer to hold array of numbers of output @ref CeedEvalMode, or `NULL` @param[out] eval_modes_out Pointer to hold arrays of output @ref CeedEvalMode, or `NULL` @param[out] eval_mode_offsets_out Pointer to hold arrays of output offsets at each quadrature point @param[out] num_output_components The number of columns in the assembled `CeedQFunction` matrix for each quadrature point, including contributions of all active bases @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorAssemblyDataGetEvalModes(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedInt **num_eval_modes_in, const CeedEvalMode ***eval_modes_in, CeedSize ***eval_mode_offsets_in, CeedInt *num_active_bases_out, CeedInt **num_eval_modes_out, const CeedEvalMode ***eval_modes_out, CeedSize ***eval_mode_offsets_out, CeedSize *num_output_components) { if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in; if (num_eval_modes_in) *num_eval_modes_in = data->num_eval_modes_in; if (eval_modes_in) *eval_modes_in = (const CeedEvalMode **)data->eval_modes_in; if (eval_mode_offsets_in) *eval_mode_offsets_in = data->eval_mode_offsets_in; if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out; if (num_eval_modes_out) *num_eval_modes_out = data->num_eval_modes_out; if (eval_modes_out) *eval_modes_out = (const CeedEvalMode **)data->eval_modes_out; if (eval_mode_offsets_out) *eval_mode_offsets_out = data->eval_mode_offsets_out; if (num_output_components) *num_output_components = data->num_output_components; return CEED_ERROR_SUCCESS; } /** @brief Get `CeedOperator` `CeedBasis` data for assembly. Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. @param[in] data `CeedOperatorAssemblyData` @param[out] num_active_bases_in Number of active input bases, or `NULL` @param[out] active_bases_in Pointer to hold active input `CeedBasis`, or `NULL` @param[out] assembled_bases_in Pointer to hold assembled active input `B` , or `NULL` @param[out] num_active_bases_out Number of active output bases, or `NULL` @param[out] active_bases_out Pointer to hold active output `CeedBasis`, or `NULL` @param[out] assembled_bases_out Pointer to hold assembled active output `B` , or `NULL` @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorAssemblyDataGetBases(CeedOperatorAssemblyData data, CeedInt *num_active_bases_in, CeedBasis **active_bases_in, const CeedScalar ***assembled_bases_in, CeedInt *num_active_bases_out, CeedBasis **active_bases_out, const CeedScalar ***assembled_bases_out) { // Assemble B_in, B_out if needed if (assembled_bases_in && !data->assembled_bases_in[0]) { CeedInt num_qpts; if (data->active_bases_in[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[0], &num_qpts)); else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_in[0], &num_qpts)); for (CeedInt b = 0; b < data->num_active_bases_in; b++) { bool has_eval_none = false; CeedInt num_nodes; CeedScalar *B_in = NULL, *identity = NULL; CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_in[b], &num_nodes)); CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_in[b], &B_in)); for (CeedInt i = 0; i < data->num_eval_modes_in[b]; i++) { has_eval_none = has_eval_none || (data->eval_modes_in[b][i] == CEED_EVAL_NONE); } if (has_eval_none) { CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { identity[i * num_nodes + i] = 1.0; } } for (CeedInt q = 0; q < num_qpts; q++) { for (CeedInt n = 0; n < num_nodes; n++) { CeedInt d_in = 0, q_comp_in; CeedEvalMode eval_mode_in_prev = CEED_EVAL_NONE; for (CeedInt e_in = 0; e_in < data->num_eval_modes_in[b]; e_in++) { const CeedInt qq = data->num_eval_modes_in[b] * q; const CeedScalar *B = NULL; CeedCall(CeedOperatorGetBasisPointer(data->active_bases_in[b], data->eval_modes_in[b][e_in], identity, &B)); CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_in[b], data->eval_modes_in[b][e_in], &q_comp_in)); if (q_comp_in > 1) { if (e_in == 0 || data->eval_modes_in[b][e_in] != eval_mode_in_prev) d_in = 0; else B = &B[(++d_in) * num_qpts * num_nodes]; } eval_mode_in_prev = data->eval_modes_in[b][e_in]; B_in[(qq + e_in) * num_nodes + n] = B[q * num_nodes + n]; } } } if (identity) CeedCall(CeedFree(&identity)); data->assembled_bases_in[b] = B_in; } } if (assembled_bases_out && !data->assembled_bases_out[0]) { CeedInt num_qpts; if (data->active_bases_out[0] == CEED_BASIS_NONE) CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[0], &num_qpts)); else CeedCall(CeedBasisGetNumQuadraturePoints(data->active_bases_out[0], &num_qpts)); for (CeedInt b = 0; b < data->num_active_bases_out; b++) { bool has_eval_none = false; CeedInt num_nodes; CeedScalar *B_out = NULL, *identity = NULL; CeedCall(CeedElemRestrictionGetElementSize(data->active_elem_rstrs_out[b], &num_nodes)); CeedCall(CeedCalloc(num_qpts * num_nodes * data->num_eval_modes_out[b], &B_out)); for (CeedInt i = 0; i < data->num_eval_modes_out[b]; i++) { has_eval_none = has_eval_none || (data->eval_modes_out[b][i] == CEED_EVAL_NONE); } if (has_eval_none) { CeedCall(CeedCalloc(num_qpts * num_nodes, &identity)); for (CeedInt i = 0; i < (num_nodes < num_qpts ? num_nodes : num_qpts); i++) { identity[i * num_nodes + i] = 1.0; } } for (CeedInt q = 0; q < num_qpts; q++) { for (CeedInt n = 0; n < num_nodes; n++) { CeedInt d_out = 0, q_comp_out; CeedEvalMode eval_mode_out_prev = CEED_EVAL_NONE; for (CeedInt e_out = 0; e_out < data->num_eval_modes_out[b]; e_out++) { const CeedInt qq = data->num_eval_modes_out[b] * q; const CeedScalar *B = NULL; CeedCall(CeedOperatorGetBasisPointer(data->active_bases_out[b], data->eval_modes_out[b][e_out], identity, &B)); CeedCall(CeedBasisGetNumQuadratureComponents(data->active_bases_out[b], data->eval_modes_out[b][e_out], &q_comp_out)); if (q_comp_out > 1) { if (e_out == 0 || data->eval_modes_out[b][e_out] != eval_mode_out_prev) d_out = 0; else B = &B[(++d_out) * num_qpts * num_nodes]; } eval_mode_out_prev = data->eval_modes_out[b][e_out]; B_out[(qq + e_out) * num_nodes + n] = B[q * num_nodes + n]; } } } if (identity) CeedCall(CeedFree(&identity)); data->assembled_bases_out[b] = B_out; } } // Pass out assembled data if (num_active_bases_in) *num_active_bases_in = data->num_active_bases_in; if (active_bases_in) *active_bases_in = data->active_bases_in; if (assembled_bases_in) *assembled_bases_in = (const CeedScalar **)data->assembled_bases_in; if (num_active_bases_out) *num_active_bases_out = data->num_active_bases_out; if (active_bases_out) *active_bases_out = data->active_bases_out; if (assembled_bases_out) *assembled_bases_out = (const CeedScalar **)data->assembled_bases_out; return CEED_ERROR_SUCCESS; } /** @brief Get `CeedOperator` `CeedBasis` data for assembly. Note: See @ref CeedOperatorAssemblyDataCreate() for a full description of the data stored in this object. @param[in] data `CeedOperatorAssemblyData` @param[out] num_active_elem_rstrs_in Number of active input element restrictions, or `NULL` @param[out] active_elem_rstrs_in Pointer to hold active input `CeedElemRestriction`, or `NULL` @param[out] num_active_elem_rstrs_out Number of active output element restrictions, or `NULL` @param[out] active_elem_rstrs_out Pointer to hold active output `CeedElemRestriction`, or `NULL` @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorAssemblyDataGetElemRestrictions(CeedOperatorAssemblyData data, CeedInt *num_active_elem_rstrs_in, CeedElemRestriction **active_elem_rstrs_in, CeedInt *num_active_elem_rstrs_out, CeedElemRestriction **active_elem_rstrs_out) { if (num_active_elem_rstrs_in) *num_active_elem_rstrs_in = data->num_active_bases_in; if (active_elem_rstrs_in) *active_elem_rstrs_in = data->active_elem_rstrs_in; if (num_active_elem_rstrs_out) *num_active_elem_rstrs_out = data->num_active_bases_out; if (active_elem_rstrs_out) *active_elem_rstrs_out = data->active_elem_rstrs_out; return CEED_ERROR_SUCCESS; } /** @brief Destroy `CeedOperatorAssemblyData` @param[in,out] data `CeedOperatorAssemblyData` to destroy @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorAssemblyDataDestroy(CeedOperatorAssemblyData *data) { if (!*data) { *data = NULL; return CEED_ERROR_SUCCESS; } CeedCall(CeedDestroy(&(*data)->ceed)); for (CeedInt b = 0; b < (*data)->num_active_bases_in; b++) { CeedCall(CeedBasisDestroy(&(*data)->active_bases_in[b])); CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_in[b])); CeedCall(CeedFree(&(*data)->eval_modes_in[b])); CeedCall(CeedFree(&(*data)->eval_mode_offsets_in[b])); CeedCall(CeedFree(&(*data)->assembled_bases_in[b])); } for (CeedInt b = 0; b < (*data)->num_active_bases_out; b++) { CeedCall(CeedBasisDestroy(&(*data)->active_bases_out[b])); CeedCall(CeedElemRestrictionDestroy(&(*data)->active_elem_rstrs_out[b])); CeedCall(CeedFree(&(*data)->eval_modes_out[b])); CeedCall(CeedFree(&(*data)->eval_mode_offsets_out[b])); CeedCall(CeedFree(&(*data)->assembled_bases_out[b])); } CeedCall(CeedFree(&(*data)->active_bases_in)); CeedCall(CeedFree(&(*data)->active_bases_out)); CeedCall(CeedFree(&(*data)->active_elem_rstrs_in)); CeedCall(CeedFree(&(*data)->active_elem_rstrs_out)); CeedCall(CeedFree(&(*data)->num_eval_modes_in)); CeedCall(CeedFree(&(*data)->num_eval_modes_out)); CeedCall(CeedFree(&(*data)->eval_modes_in)); CeedCall(CeedFree(&(*data)->eval_modes_out)); CeedCall(CeedFree(&(*data)->eval_mode_offsets_in)); CeedCall(CeedFree(&(*data)->eval_mode_offsets_out)); CeedCall(CeedFree(&(*data)->assembled_bases_in)); CeedCall(CeedFree(&(*data)->assembled_bases_out)); CeedCall(CeedFree(data)); return CEED_ERROR_SUCCESS; } /** @brief Retrieve fallback `CeedOperator` with a reference `Ceed` for advanced `CeedOperator` functionality @param[in] op `CeedOperator` to retrieve fallback for @param[out] op_fallback Fallback `CeedOperator` @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorGetFallback(CeedOperator op, CeedOperator *op_fallback) { // Create if needed if (!op->op_fallback) CeedCall(CeedOperatorCreateFallback(op)); if (op->op_fallback) { bool is_debug; Ceed ceed; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedIsDebug(ceed, &is_debug)); if (is_debug) { Ceed ceed_fallback; const char *resource, *resource_fallback; CeedCall(CeedGetOperatorFallbackCeed(ceed, &ceed_fallback)); CeedCall(CeedGetResource(ceed, &resource)); CeedCall(CeedGetResource(ceed_fallback, &resource_fallback)); CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "---------- CeedOperator Fallback ----------\n"); CeedDebug(ceed, "Falling back from %s operator at address %p to %s operator at address %p\n", resource, op, resource_fallback, op->op_fallback); } } *op_fallback = op->op_fallback; return CEED_ERROR_SUCCESS; } /** @brief Get the parent `CeedOperator` for a fallback `CeedOperator` @param[in] op `CeedOperator` context @param[out] parent Variable to store parent `CeedOperator` context @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorGetFallbackParent(CeedOperator op, CeedOperator *parent) { *parent = op->op_fallback_parent ? op->op_fallback_parent : NULL; return CEED_ERROR_SUCCESS; } /** @brief Get the `Ceed` context of the parent `CeedOperator` for a fallback `CeedOperator` @param[in] op `CeedOperator` context @param[out] parent Variable to store parent `Ceed` context @return An error code: 0 - success, otherwise - failure @ref Backend **/ int CeedOperatorGetFallbackParentCeed(CeedOperator op, Ceed *parent) { *parent = op->op_fallback_parent ? op->op_fallback_parent->ceed : op->ceed; return CEED_ERROR_SUCCESS; } /// @} /// ---------------------------------------------------------------------------- /// CeedOperator Public API /// ---------------------------------------------------------------------------- /// @addtogroup CeedOperatorUser /// @{ /** @brief Assemble a linear `CeedQFunction` associated with a `CeedOperator`. This returns a `CeedVector` containing a matrix at each quadrature point providing the action of the `CeedQFunction` associated with the `CeedOperator`. 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. Inputs and outputs are in the order provided by the user when adding `CeedOperator` fields. 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]`. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble `CeedQFunction` @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorLinearAssembleQFunction(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { CeedCall(CeedOperatorCheckReady(op)); if (op->LinearAssembleQFunction) { // Backend version CeedCall(op->LinearAssembleQFunction(op, assembled, rstr, request)); } else { // Operator fallback Ceed ceed; CeedOperator op_fallback; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunction(op_fallback, assembled, rstr, request)); else return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunction"); } return CEED_ERROR_SUCCESS; } /** @brief Assemble `CeedQFunction` and store result internally. Return copied references of stored data to the caller. Caller is responsible for ownership and destruction of the copied references. See also @ref CeedOperatorLinearAssembleQFunction(). Note: If the value of `assembled` or `rstr` passed to this function are non-`NULL` , then it is assumed that they hold valid pointers. These objects will be destroyed if `*assembled` or `*rstr` is the only reference to the object. @param[in] op `CeedOperator` to assemble `CeedQFunction` @param[out] assembled `CeedVector` to store assembled `CeedQFunction` at quadrature points @param[out] rstr `CeedElemRestriction` for `CeedVector` containing assembled `CeedQFunction` @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorLinearAssembleQFunctionBuildOrUpdate(CeedOperator op, CeedVector *assembled, CeedElemRestriction *rstr, CeedRequest *request) { int (*LinearAssembleQFunctionUpdate)(CeedOperator, CeedVector, CeedElemRestriction, CeedRequest *) = NULL; CeedOperator op_assemble = NULL; CeedOperator op_fallback_parent = NULL; CeedCall(CeedOperatorCheckReady(op)); // Determine if fallback parent or operator has implementation CeedCall(CeedOperatorGetFallbackParent(op, &op_fallback_parent)); if (op_fallback_parent && op_fallback_parent->LinearAssembleQFunctionUpdate) { // -- Backend version for op fallback parent is faster, if it exists LinearAssembleQFunctionUpdate = op_fallback_parent->LinearAssembleQFunctionUpdate; op_assemble = op_fallback_parent; } else if (op->LinearAssembleQFunctionUpdate) { // -- Backend version for op LinearAssembleQFunctionUpdate = op->LinearAssembleQFunctionUpdate; op_assemble = op; } // Assemble QFunction if (LinearAssembleQFunctionUpdate) { // Backend or fallback parent version bool qf_assembled_is_setup; CeedVector assembled_vec = NULL; CeedElemRestriction assembled_rstr = NULL; CeedCall(CeedQFunctionAssemblyDataIsSetup(op->qf_assembled, &qf_assembled_is_setup)); if (qf_assembled_is_setup) { bool update_needed; CeedCall(CeedQFunctionAssemblyDataGetObjects(op->qf_assembled, &assembled_vec, &assembled_rstr)); CeedCall(CeedQFunctionAssemblyDataIsUpdateNeeded(op->qf_assembled, &update_needed)); if (update_needed) CeedCall(LinearAssembleQFunctionUpdate(op_assemble, assembled_vec, assembled_rstr, request)); } else { CeedCall(CeedOperatorLinearAssembleQFunction(op_assemble, &assembled_vec, &assembled_rstr, request)); CeedCall(CeedQFunctionAssemblyDataSetObjects(op->qf_assembled, assembled_vec, assembled_rstr)); } CeedCall(CeedQFunctionAssemblyDataSetUpdateNeeded(op->qf_assembled, false)); // Copy reference from internally held copy CeedCall(CeedVectorReferenceCopy(assembled_vec, assembled)); CeedCall(CeedElemRestrictionReferenceCopy(assembled_rstr, rstr)); CeedCall(CeedVectorDestroy(&assembled_vec)); CeedCall(CeedElemRestrictionDestroy(&assembled_rstr)); } else { // Operator fallback Ceed ceed; CeedOperator op_fallback; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op_fallback, assembled, rstr, request)); else return CeedError(ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support CeedOperatorLinearAssembleQFunctionUpdate"); } return CEED_ERROR_SUCCESS; } /** @brief Assemble the diagonal of a square linear `CeedOperator` This overwrites a `CeedVector` with the diagonal of a linear `CeedOperator`. Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble `CeedQFunction` @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorLinearAssembleDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { bool is_composite; CeedSize input_size = 0, output_size = 0; Ceed ceed; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square"); // Early exit for empty operator if (!is_composite) { CeedInt num_elem = 0; CeedCall(CeedOperatorGetNumElements(op, &num_elem)); if (num_elem == 0) return CEED_ERROR_SUCCESS; } if (op->LinearAssembleDiagonal) { // Backend version CeedCall(op->LinearAssembleDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } else if (op->LinearAssembleAddDiagonal) { // Backend version with zeroing first CeedCall(CeedVectorSetValue(assembled, 0.0)); CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorLinearAssembleDiagonal(op_fallback, assembled, request)); return CEED_ERROR_SUCCESS; } } // Default interface implementation CeedCall(CeedVectorSetValue(assembled, 0.0)); CeedCall(CeedOperatorLinearAssembleAddDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } /** @brief Assemble the diagonal of a square linear `CeedOperator`. This sums into a `CeedVector` the diagonal of a linear `CeedOperator`. Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. Note: Calling this function asserts that setup is complete and sets the CeedOperator as immutable. @param[in] op `CeedOperator` to assemble `CeedQFunction` @param[out] assembled `CeedVector` to store assembled `CeedOperator` diagonal @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorLinearAssembleAddDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { bool is_composite; CeedSize input_size = 0, output_size = 0; Ceed ceed; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square"); // Early exit for empty operator if (!is_composite) { CeedInt num_elem = 0; CeedCall(CeedOperatorGetNumElements(op, &num_elem)); if (num_elem == 0) return CEED_ERROR_SUCCESS; } if (op->LinearAssembleAddDiagonal) { // Backend version CeedCall(op->LinearAssembleAddDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorLinearAssembleAddDiagonal(op_fallback, assembled, request)); return CEED_ERROR_SUCCESS; } } // Default interface implementation if (is_composite) { CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, false, assembled)); } else { CeedCall(CeedSingleOperatorAssembleAddDiagonal_Core(op, request, false, assembled)); } return CEED_ERROR_SUCCESS; } /** @brief Fully assemble the point-block diagonal pattern of a linear `CeedOperator`. Expected to be used in conjunction with @ref CeedOperatorLinearAssemblePointBlockDiagonal(). 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)`. Note that the `(i, j)` pairs are unique. This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemblePointBlockDiagonal() provides the values in the same ordering. This will generally be slow unless your operator is low-order. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble @param[out] num_entries Number of entries in coordinate nonzero pattern @param[out] rows Row number for each entry @param[out] cols Column number for each entry @ref User **/ int CeedOperatorLinearAssemblePointBlockDiagonalSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { Ceed ceed; bool is_composite; CeedInt num_active_components, num_sub_operators; CeedOperator *sub_operators; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedSize input_size = 0, output_size = 0; CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square"); if (is_composite) { CeedCall(CeedCompositeOperatorGetNumSub(op, &num_sub_operators)); CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); } else { sub_operators = &op; num_sub_operators = 1; } // Verify operator can be assembled correctly { CeedOperatorAssemblyData data; CeedInt num_active_elem_rstrs, comp_stride; CeedElemRestriction *active_elem_rstrs; // Get initial values to check against CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[0], &data)); CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL)); CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[0], &comp_stride)); CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[0], &num_active_components)); // Verify that all active element restrictions have same component stride and number of components for (CeedInt k = 0; k < num_sub_operators; k++) { CeedCall(CeedOperatorGetOperatorAssemblyData(sub_operators[k], &data)); CeedCall(CeedOperatorAssemblyDataGetElemRestrictions(data, &num_active_elem_rstrs, &active_elem_rstrs, NULL, NULL)); for (CeedInt i = 0; i < num_active_elem_rstrs; i++) { CeedInt comp_stride_sub, num_active_components_sub; CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstrs[i], &comp_stride_sub)); CeedCheck(comp_stride == comp_stride_sub, ceed, CEED_ERROR_DIMENSION, "Active element restrictions must have the same component stride: %d vs %d", comp_stride, comp_stride_sub); CeedCall(CeedElemRestrictionGetNumComponents(active_elem_rstrs[i], &num_active_components_sub)); CeedCheck(num_active_components == num_active_components_sub, ceed, CEED_ERROR_INCOMPATIBLE, "All suboperators must have the same number of output components"); } } } *num_entries = input_size * num_active_components; CeedCall(CeedCalloc(*num_entries, rows)); CeedCall(CeedCalloc(*num_entries, cols)); for (CeedInt o = 0; o < num_sub_operators; o++) { CeedElemRestriction active_elem_rstr, point_block_active_elem_rstr; CeedInt comp_stride, num_elem, elem_size; const CeedInt *offsets, *point_block_offsets; CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[o], &active_elem_rstr)); CeedCall(CeedElemRestrictionGetCompStride(active_elem_rstr, &comp_stride)); CeedCall(CeedElemRestrictionGetNumElements(active_elem_rstr, &num_elem)); CeedCall(CeedElemRestrictionGetElementSize(active_elem_rstr, &elem_size)); CeedCall(CeedElemRestrictionGetOffsets(active_elem_rstr, CEED_MEM_HOST, &offsets)); CeedCall(CeedOperatorCreateActivePointBlockRestriction(active_elem_rstr, &point_block_active_elem_rstr)); CeedCall(CeedElemRestrictionGetOffsets(point_block_active_elem_rstr, CEED_MEM_HOST, &point_block_offsets)); for (CeedSize i = 0; i < num_elem * elem_size; i++) { for (CeedInt c_out = 0; c_out < num_active_components; c_out++) { for (CeedInt c_in = 0; c_in < num_active_components; c_in++) { (*rows)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_out * comp_stride; (*cols)[point_block_offsets[i] + c_out * num_active_components + c_in] = offsets[i] + c_in * comp_stride; } } } CeedCall(CeedElemRestrictionRestoreOffsets(active_elem_rstr, &offsets)); CeedCall(CeedElemRestrictionRestoreOffsets(point_block_active_elem_rstr, &point_block_offsets)); CeedCall(CeedElemRestrictionDestroy(&point_block_active_elem_rstr)); } return CEED_ERROR_SUCCESS; } /** @brief Assemble the point block diagonal of a square linear `CeedOperator`. This overwrites a `CeedVector` with the point block diagonal of a linear `CeedOperator`. Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble `CeedQFunction` @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. The dimensions of this vector are derived from the active vector for the `CeedOperator`. The array has shape `[nodes, component out, component in]`. @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorLinearAssemblePointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { bool is_composite; CeedSize input_size = 0, output_size = 0; Ceed ceed; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square"); // Early exit for empty operator if (!is_composite) { CeedInt num_elem = 0; CeedCall(CeedOperatorGetNumElements(op, &num_elem)); if (num_elem == 0) return CEED_ERROR_SUCCESS; } if (op->LinearAssemblePointBlockDiagonal) { // Backend version CeedCall(op->LinearAssemblePointBlockDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } else if (op->LinearAssembleAddPointBlockDiagonal) { // Backend version with zeroing first CeedCall(CeedVectorSetValue(assembled, 0.0)); CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorLinearAssemblePointBlockDiagonal(op_fallback, assembled, request)); return CEED_ERROR_SUCCESS; } } // Default interface implementation CeedCall(CeedVectorSetValue(assembled, 0.0)); CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } /** @brief Assemble the point block diagonal of a square linear `CeedOperator`. This sums into a `CeedVector` with the point block diagonal of a linear `CeedOperator`. Note: Currently only non-composite `CeedOperator` with a single field and composite `CeedOperator` with single field sub-operators are supported. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble `CeedQFunction` @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. The dimensions of this vector are derived from the active vector for the `CeedOperator`. The array has shape `[nodes, component out, component in]`. @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorLinearAssembleAddPointBlockDiagonal(CeedOperator op, CeedVector assembled, CeedRequest *request) { bool is_composite; CeedSize input_size = 0, output_size = 0; Ceed ceed; CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); CeedCall(CeedOperatorGetActiveVectorLengths(op, &input_size, &output_size)); CeedCheck(input_size == output_size, ceed, CEED_ERROR_DIMENSION, "Operator must be square"); // Early exit for empty operator if (!is_composite) { CeedInt num_elem = 0; CeedCall(CeedOperatorGetNumElements(op, &num_elem)); if (num_elem == 0) return CEED_ERROR_SUCCESS; } if (op->LinearAssembleAddPointBlockDiagonal) { // Backend version CeedCall(op->LinearAssembleAddPointBlockDiagonal(op, assembled, request)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorLinearAssembleAddPointBlockDiagonal(op_fallback, assembled, request)); return CEED_ERROR_SUCCESS; } } // Default interface implementation if (is_composite) { CeedCall(CeedCompositeOperatorLinearAssembleAddDiagonal(op, request, true, assembled)); } else { CeedCall(CeedSingleOperatorAssembleAddDiagonal_Core(op, request, true, assembled)); } return CEED_ERROR_SUCCESS; } /** @brief Fully assemble the nonzero pattern of a linear `CeedOperator`. Expected to be used in conjunction with @ref CeedOperatorLinearAssemble(). 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)`. Note that the `(i, j)` pairs are not unique and may repeat. This function returns the number of entries and their `(i, j)` locations, while @ref CeedOperatorLinearAssemble() provides the values in the same ordering. This will generally be slow unless your operator is low-order. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble @param[out] num_entries Number of entries in coordinate nonzero pattern @param[out] rows Row number for each entry @param[out] cols Column number for each entry @ref User **/ int CeedOperatorLinearAssembleSymbolic(CeedOperator op, CeedSize *num_entries, CeedInt **rows, CeedInt **cols) { bool is_composite; CeedInt num_suboperators, offset = 0; CeedSize single_entries; CeedOperator *sub_operators; CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); if (op->LinearAssembleSymbolic) { // Backend version CeedCall(op->LinearAssembleSymbolic(op, num_entries, rows, cols)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorLinearAssembleSymbolic(op_fallback, num_entries, rows, cols)); return CEED_ERROR_SUCCESS; } } // Default interface implementation // Count entries and allocate rows, cols arrays *num_entries = 0; if (is_composite) { CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); for (CeedInt k = 0; k < num_suboperators; ++k) { CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries)); *num_entries += single_entries; } } else { CeedCall(CeedSingleOperatorAssemblyCountEntries(op, &single_entries)); *num_entries += single_entries; } CeedCall(CeedCalloc(*num_entries, rows)); CeedCall(CeedCalloc(*num_entries, cols)); // Assemble nonzero locations if (is_composite) { CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); for (CeedInt k = 0; k < num_suboperators; ++k) { CeedCall(CeedSingleOperatorAssembleSymbolic(sub_operators[k], offset, *rows, *cols)); CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries)); offset += single_entries; } } else { CeedCall(CeedSingleOperatorAssembleSymbolic(op, offset, *rows, *cols)); } return CEED_ERROR_SUCCESS; } /** @brief Fully assemble the nonzero entries of a linear operator. Expected to be used in conjunction with @ref CeedOperatorLinearAssembleSymbolic(). 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)`. Note that the `(i, j)` pairs are not unique and may repeat. This function returns the values of the nonzero entries to be added, their `(i, j)` locations are provided by @ref CeedOperatorLinearAssembleSymbolic(). This will generally be slow unless your operator is low-order. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to assemble @param[out] values Values to assemble into matrix @ref User **/ int CeedOperatorLinearAssemble(CeedOperator op, CeedVector values) { bool is_composite; CeedInt num_suboperators, offset = 0; CeedSize single_entries = 0; CeedOperator *sub_operators; CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorIsComposite(op, &is_composite)); // Early exit for empty operator if (!is_composite) { CeedInt num_elem = 0; CeedCall(CeedOperatorGetNumElements(op, &num_elem)); if (num_elem == 0) return CEED_ERROR_SUCCESS; } if (op->LinearAssemble) { // Backend version CeedCall(op->LinearAssemble(op, values)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorLinearAssemble(op_fallback, values)); return CEED_ERROR_SUCCESS; } } // Default interface implementation CeedCall(CeedVectorSetValue(values, 0.0)); if (is_composite) { CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); for (CeedInt k = 0; k < num_suboperators; k++) { CeedCall(CeedSingleOperatorAssemble(sub_operators[k], offset, values)); CeedCall(CeedSingleOperatorAssemblyCountEntries(sub_operators[k], &single_entries)); offset += single_entries; } } else { CeedCall(CeedSingleOperatorAssemble(op, offset, values)); } return CEED_ERROR_SUCCESS; } /** @brief Get the multiplicity of nodes across sub-operators in a composite `CeedOperator`. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op Composite `CeedOperator` @param[in] num_skip_indices Number of sub-operators to skip @param[in] skip_indices Array of indices of sub-operators to skip @param[out] mult Vector to store multiplicity (of size `l_size` ) @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedCompositeOperatorGetMultiplicity(CeedOperator op, CeedInt num_skip_indices, CeedInt *skip_indices, CeedVector mult) { Ceed ceed; CeedInt num_suboperators; CeedSize l_vec_len; CeedScalar *mult_array; CeedVector ones_l_vec; CeedElemRestriction elem_rstr, mult_elem_rstr; CeedOperator *sub_operators; CeedCall(CeedOperatorCheckReady(op)); CeedCall(CeedOperatorGetCeed(op, &ceed)); // Zero mult vector CeedCall(CeedVectorSetValue(mult, 0.0)); // Get suboperators CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); if (num_suboperators == 0) return CEED_ERROR_SUCCESS; // Work vector CeedCall(CeedVectorGetLength(mult, &l_vec_len)); CeedCall(CeedVectorCreate(ceed, l_vec_len, &ones_l_vec)); CeedCall(CeedVectorSetValue(ones_l_vec, 1.0)); CeedCall(CeedVectorGetArray(mult, CEED_MEM_HOST, &mult_array)); // Compute multiplicity across suboperators for (CeedInt i = 0; i < num_suboperators; i++) { const CeedScalar *sub_mult_array; CeedVector sub_mult_l_vec, ones_e_vec; // -- Check for suboperator to skip for (CeedInt j = 0; j < num_skip_indices; j++) { if (skip_indices[j] == i) continue; } // -- Sub operator multiplicity CeedCall(CeedOperatorGetActiveElemRestriction(sub_operators[i], &elem_rstr)); CeedCall(CeedElemRestrictionCreateUnorientedCopy(elem_rstr, &mult_elem_rstr)); CeedCall(CeedElemRestrictionCreateVector(mult_elem_rstr, &sub_mult_l_vec, &ones_e_vec)); CeedCall(CeedVectorSetValue(sub_mult_l_vec, 0.0)); CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_NOTRANSPOSE, ones_l_vec, ones_e_vec, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedElemRestrictionApply(mult_elem_rstr, CEED_TRANSPOSE, ones_e_vec, sub_mult_l_vec, CEED_REQUEST_IMMEDIATE)); CeedCall(CeedVectorGetArrayRead(sub_mult_l_vec, CEED_MEM_HOST, &sub_mult_array)); // ---- Flag every node present in the current suboperator for (CeedInt j = 0; j < l_vec_len; j++) { if (sub_mult_array[j] > 0.0) mult_array[j] += 1.0; } CeedCall(CeedVectorRestoreArrayRead(sub_mult_l_vec, &sub_mult_array)); CeedCall(CeedVectorDestroy(&sub_mult_l_vec)); CeedCall(CeedVectorDestroy(&ones_e_vec)); CeedCall(CeedElemRestrictionDestroy(&mult_elem_rstr)); } CeedCall(CeedVectorRestoreArray(mult, &mult_array)); CeedCall(CeedVectorDestroy(&ones_l_vec)); return CEED_ERROR_SUCCESS; } /** @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator`, creating the prolongation basis from the fine and coarse grid interpolation. Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. @param[in] op_fine Fine grid `CeedOperator` @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` @param[in] rstr_coarse Coarse grid `CeedElemRestriction` @param[in] basis_coarse Coarse grid active vector `CeedBasis` @param[out] op_coarse Coarse grid `CeedOperator` @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorMultigridLevelCreate(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { CeedBasis basis_c_to_f = NULL; CeedCall(CeedOperatorCheckReady(op_fine)); // Build prolongation matrix, if required if (op_prolong || op_restrict) { CeedBasis basis_fine; CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); CeedCall(CeedBasisCreateProjection(basis_coarse, basis_fine, &basis_c_to_f)); } // Core code CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict)); return CEED_ERROR_SUCCESS; } /** @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a tensor basis for the active basis. Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. @param[in] op_fine Fine grid `CeedOperator` @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` @param[in] rstr_coarse Coarse grid `CeedElemRestriction` @param[in] basis_coarse Coarse grid active vector `CeedBasis` @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator` @param[out] op_coarse Coarse grid `CeedOperator` @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorMultigridLevelCreateTensorH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { Ceed ceed; CeedInt Q_f, Q_c; CeedBasis basis_fine, basis_c_to_f = NULL; CeedCall(CeedOperatorCheckReady(op_fine)); CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); // Check for compatible quadrature spaces CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces"); // Create coarse to fine basis, if required if (op_prolong || op_restrict) { CeedInt dim, num_comp, num_nodes_c, P_1d_f, P_1d_c; CeedScalar *q_ref, *q_weight, *grad; // Check if interpolation matrix is provided CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix"); CeedCall(CeedBasisGetDimension(basis_fine, &dim)); CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); CeedCall(CeedBasisGetNumNodes1D(basis_fine, &P_1d_f)); CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); P_1d_c = dim == 1 ? num_nodes_c : dim == 2 ? sqrt(num_nodes_c) : cbrt(num_nodes_c); CeedCall(CeedCalloc(P_1d_f, &q_ref)); CeedCall(CeedCalloc(P_1d_f, &q_weight)); CeedCall(CeedCalloc(P_1d_f * P_1d_c * dim, &grad)); 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)); CeedCall(CeedFree(&q_ref)); CeedCall(CeedFree(&q_weight)); CeedCall(CeedFree(&grad)); } // Core code CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict)); return CEED_ERROR_SUCCESS; } /** @brief Create a multigrid coarse `CeedOperator` and level transfer `CeedOperator` for a `CeedOperator` with a non-tensor basis for the active vector Note: Calling this function asserts that setup is complete and sets all four `CeedOperator` as immutable. @param[in] op_fine Fine grid `CeedOperator` @param[in] p_mult_fine L-vector multiplicity in parallel gather/scatter, or `NULL` if not creating prolongation/restriction `CeedOperator` @param[in] rstr_coarse Coarse grid `CeedElemRestriction` @param[in] basis_coarse Coarse grid active vector `CeedBasis` @param[in] interp_c_to_f Matrix for coarse to fine interpolation, or `NULL` if not creating prolongation/restriction `CeedOperator` @param[out] op_coarse Coarse grid `CeedOperator` @param[out] op_prolong Coarse to fine `CeedOperator`, or `NULL` @param[out] op_restrict Fine to coarse `CeedOperator`, or `NULL` @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorMultigridLevelCreateH1(CeedOperator op_fine, CeedVector p_mult_fine, CeedElemRestriction rstr_coarse, CeedBasis basis_coarse, const CeedScalar *interp_c_to_f, CeedOperator *op_coarse, CeedOperator *op_prolong, CeedOperator *op_restrict) { Ceed ceed; CeedInt Q_f, Q_c; CeedBasis basis_fine, basis_c_to_f = NULL; CeedCall(CeedOperatorCheckReady(op_fine)); CeedCall(CeedOperatorGetCeed(op_fine, &ceed)); // Check for compatible quadrature spaces CeedCall(CeedOperatorGetActiveBasis(op_fine, &basis_fine)); CeedCall(CeedBasisGetNumQuadraturePoints(basis_fine, &Q_f)); CeedCall(CeedBasisGetNumQuadraturePoints(basis_coarse, &Q_c)); CeedCheck(Q_f == Q_c, ceed, CEED_ERROR_DIMENSION, "Bases must have compatible quadrature spaces"); // Coarse to fine basis if (op_prolong || op_restrict) { CeedInt dim, num_comp, num_nodes_c, num_nodes_f; CeedScalar *q_ref, *q_weight, *grad; CeedElemTopology topo; // Check if interpolation matrix is provided CeedCheck(interp_c_to_f, ceed, CEED_ERROR_INCOMPATIBLE, "Prolongation or restriction operator creation requires coarse-to-fine interpolation matrix"); CeedCall(CeedBasisGetTopology(basis_fine, &topo)); CeedCall(CeedBasisGetDimension(basis_fine, &dim)); CeedCall(CeedBasisGetNumComponents(basis_fine, &num_comp)); CeedCall(CeedBasisGetNumNodes(basis_fine, &num_nodes_f)); CeedCall(CeedElemRestrictionGetElementSize(rstr_coarse, &num_nodes_c)); CeedCall(CeedCalloc(num_nodes_f * dim, &q_ref)); CeedCall(CeedCalloc(num_nodes_f, &q_weight)); CeedCall(CeedCalloc(num_nodes_f * num_nodes_c * dim, &grad)); 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)); CeedCall(CeedFree(&q_ref)); CeedCall(CeedFree(&q_weight)); CeedCall(CeedFree(&grad)); } // Core code CeedCall(CeedSingleOperatorMultigridLevel(op_fine, p_mult_fine, rstr_coarse, basis_coarse, basis_c_to_f, op_coarse, op_prolong, op_restrict)); return CEED_ERROR_SUCCESS; } /** @brief Build a FDM based approximate inverse for each element for a `CeedOperator`. This returns a `CeedOperator` and `CeedVector` to apply a Fast Diagonalization Method based approximate inverse. This function obtains the simultaneous diagonalization for the 1D mass and Laplacian operators, \f$M = V^T V, K = V^T S V\f$. 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$. The `CeedOperator` must be linear and non-composite. The associated `CeedQFunction` must therefore also be linear. Note: Calling this function asserts that setup is complete and sets the `CeedOperator` as immutable. @param[in] op `CeedOperator` to create element inverses @param[out] fdm_inv `CeedOperator` to apply the action of a FDM based inverse for each element @param[in] request Address of @ref CeedRequest for non-blocking completion, else @ref CEED_REQUEST_IMMEDIATE @return An error code: 0 - success, otherwise - failure @ref User **/ int CeedOperatorCreateFDMElementInverse(CeedOperator op, CeedOperator *fdm_inv, CeedRequest *request) { Ceed ceed, ceed_parent; bool interp = false, grad = false, is_tensor_basis = true; CeedInt num_input_fields, P_1d, Q_1d, num_nodes, num_qpts, dim, num_comp = 1, num_elem = 1; CeedSize l_size = 1; CeedScalar *mass, *laplace, *x, *fdm_interp, *lambda, *elem_avg; const CeedScalar *interp_1d, *grad_1d, *q_weight_1d; CeedVector q_data; CeedElemRestriction rstr = NULL, rstr_qd_i; CeedBasis basis = NULL, fdm_basis; CeedQFunctionContext ctx_fdm; CeedQFunctionField *qf_fields; CeedQFunction qf, qf_fdm; CeedOperatorField *op_fields; CeedCall(CeedOperatorCheckReady(op)); if (op->CreateFDMElementInverse) { // Backend version CeedCall(op->CreateFDMElementInverse(op, fdm_inv, request)); return CEED_ERROR_SUCCESS; } else { // Operator fallback CeedOperator op_fallback; CeedCall(CeedOperatorGetFallback(op, &op_fallback)); if (op_fallback) { CeedCall(CeedOperatorCreateFDMElementInverse(op_fallback, fdm_inv, request)); return CEED_ERROR_SUCCESS; } } // Default interface implementation CeedCall(CeedOperatorGetCeed(op, &ceed)); CeedCall(CeedOperatorGetFallbackParentCeed(op, &ceed_parent)); CeedCall(CeedOperatorGetQFunction(op, &qf)); // Determine active input basis CeedCall(CeedOperatorGetFields(op, &num_input_fields, &op_fields, NULL, NULL)); CeedCall(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); for (CeedInt i = 0; i < num_input_fields; i++) { CeedVector vec; CeedCall(CeedOperatorFieldGetVector(op_fields[i], &vec)); if (vec == CEED_VECTOR_ACTIVE) { CeedEvalMode eval_mode; CeedCall(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); interp = interp || eval_mode == CEED_EVAL_INTERP; grad = grad || eval_mode == CEED_EVAL_GRAD; CeedCall(CeedOperatorFieldGetBasis(op_fields[i], &basis)); CeedCall(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr)); } } CeedCheck(basis, ceed, CEED_ERROR_BACKEND, "No active field set"); CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); CeedCall(CeedBasisGetDimension(basis, &dim)); CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); CeedCall(CeedElemRestrictionGetNumElements(rstr, &num_elem)); CeedCall(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); // Build and diagonalize 1D Mass and Laplacian CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); CeedCheck(is_tensor_basis, ceed, CEED_ERROR_BACKEND, "FDMElementInverse only supported for tensor bases"); CeedCall(CeedCalloc(P_1d * P_1d, &mass)); CeedCall(CeedCalloc(P_1d * P_1d, &laplace)); CeedCall(CeedCalloc(P_1d * P_1d, &x)); CeedCall(CeedCalloc(P_1d * P_1d, &fdm_interp)); CeedCall(CeedCalloc(P_1d, &lambda)); // -- Build matrices CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d)); CeedCall(CeedBuildMassLaplace(interp_1d, grad_1d, q_weight_1d, P_1d, Q_1d, dim, mass, laplace)); // -- Diagonalize CeedCall(CeedSimultaneousDiagonalization(ceed, laplace, mass, x, lambda, P_1d)); CeedCall(CeedFree(&mass)); CeedCall(CeedFree(&laplace)); for (CeedInt i = 0; i < P_1d; i++) { for (CeedInt j = 0; j < P_1d; j++) fdm_interp[i + j * P_1d] = x[j + i * P_1d]; } CeedCall(CeedFree(&x)); { CeedInt layout[3], num_modes = (interp ? 1 : 0) + (grad ? dim : 0); CeedScalar max_norm = 0; const CeedScalar *assembled_array, *q_weight_array; CeedVector assembled = NULL, q_weight; CeedElemRestriction rstr_qf = NULL; // Assemble QFunction CeedCall(CeedOperatorLinearAssembleQFunctionBuildOrUpdate(op, &assembled, &rstr_qf, request)); CeedCall(CeedElemRestrictionGetELayout(rstr_qf, layout)); CeedCall(CeedElemRestrictionDestroy(&rstr_qf)); CeedCall(CeedVectorNorm(assembled, CEED_NORM_MAX, &max_norm)); // Calculate element averages CeedCall(CeedVectorCreate(ceed_parent, num_qpts, &q_weight)); CeedCall(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, q_weight)); CeedCall(CeedVectorGetArrayRead(assembled, CEED_MEM_HOST, &assembled_array)); CeedCall(CeedVectorGetArrayRead(q_weight, CEED_MEM_HOST, &q_weight_array)); CeedCall(CeedCalloc(num_elem, &elem_avg)); const CeedScalar qf_value_bound = max_norm * 100 * CEED_EPSILON; for (CeedInt e = 0; e < num_elem; e++) { CeedInt count = 0; for (CeedInt q = 0; q < num_qpts; q++) { for (CeedInt i = 0; i < num_comp * num_comp * num_modes * num_modes; i++) { if (fabs(assembled_array[q * layout[0] + i * layout[1] + e * layout[2]]) > qf_value_bound) { elem_avg[e] += assembled_array[q * layout[0] + i * layout[1] + e * layout[2]] / q_weight_array[q]; count++; } } } if (count) { elem_avg[e] /= count; } else { elem_avg[e] = 1.0; } } CeedCall(CeedVectorRestoreArrayRead(assembled, &assembled_array)); CeedCall(CeedVectorDestroy(&assembled)); CeedCall(CeedVectorRestoreArrayRead(q_weight, &q_weight_array)); CeedCall(CeedVectorDestroy(&q_weight)); } // Build FDM diagonal { CeedScalar *q_data_array, *fdm_diagonal; CeedCall(CeedCalloc(num_comp * num_nodes, &fdm_diagonal)); const CeedScalar fdm_diagonal_bound = num_nodes * CEED_EPSILON; for (CeedInt c = 0; c < num_comp; c++) { for (CeedInt n = 0; n < num_nodes; n++) { if (interp) fdm_diagonal[c * num_nodes + n] = 1.0; if (grad) { for (CeedInt d = 0; d < dim; d++) { CeedInt i = (n / CeedIntPow(P_1d, d)) % P_1d; fdm_diagonal[c * num_nodes + n] += lambda[i]; } } if (fabs(fdm_diagonal[c * num_nodes + n]) < fdm_diagonal_bound) fdm_diagonal[c * num_nodes + n] = fdm_diagonal_bound; } } CeedCall(CeedVectorCreate(ceed_parent, num_elem * num_comp * num_nodes, &q_data)); CeedCall(CeedVectorSetValue(q_data, 0.0)); CeedCall(CeedVectorGetArrayWrite(q_data, CEED_MEM_HOST, &q_data_array)); for (CeedInt e = 0; e < num_elem; e++) { for (CeedInt c = 0; c < num_comp; c++) { for (CeedInt n = 0; n < num_nodes; n++) q_data_array[(e * num_comp + c) * num_nodes + n] = 1. / (elem_avg[e] * fdm_diagonal[c * num_nodes + n]); } } CeedCall(CeedFree(&elem_avg)); CeedCall(CeedFree(&fdm_diagonal)); CeedCall(CeedVectorRestoreArray(q_data, &q_data_array)); } // Setup FDM operator // -- Basis { CeedScalar *grad_dummy, *q_ref_dummy, *q_weight_dummy; CeedCall(CeedCalloc(P_1d * P_1d, &grad_dummy)); CeedCall(CeedCalloc(P_1d, &q_ref_dummy)); CeedCall(CeedCalloc(P_1d, &q_weight_dummy)); CeedCall(CeedBasisCreateTensorH1(ceed_parent, dim, num_comp, P_1d, P_1d, fdm_interp, grad_dummy, q_ref_dummy, q_weight_dummy, &fdm_basis)); CeedCall(CeedFree(&fdm_interp)); CeedCall(CeedFree(&grad_dummy)); CeedCall(CeedFree(&q_ref_dummy)); CeedCall(CeedFree(&q_weight_dummy)); CeedCall(CeedFree(&lambda)); } // -- Restriction { CeedInt strides[3] = {1, num_nodes, num_nodes * num_comp}; CeedCall(CeedElemRestrictionCreateStrided(ceed_parent, num_elem, num_nodes, num_comp, num_elem * num_comp * num_nodes, strides, &rstr_qd_i)); } // -- QFunction CeedCall(CeedQFunctionCreateInteriorByName(ceed_parent, "Scale", &qf_fdm)); CeedCall(CeedQFunctionAddInput(qf_fdm, "input", num_comp, CEED_EVAL_INTERP)); CeedCall(CeedQFunctionAddInput(qf_fdm, "scale", num_comp, CEED_EVAL_NONE)); CeedCall(CeedQFunctionAddOutput(qf_fdm, "output", num_comp, CEED_EVAL_INTERP)); CeedCall(CeedQFunctionSetUserFlopsEstimate(qf_fdm, num_comp)); // -- QFunction context { CeedInt *num_comp_data; CeedCall(CeedCalloc(1, &num_comp_data)); num_comp_data[0] = num_comp; CeedCall(CeedQFunctionContextCreate(ceed, &ctx_fdm)); CeedCall(CeedQFunctionContextSetData(ctx_fdm, CEED_MEM_HOST, CEED_OWN_POINTER, sizeof(*num_comp_data), num_comp_data)); } CeedCall(CeedQFunctionSetContext(qf_fdm, ctx_fdm)); CeedCall(CeedQFunctionContextDestroy(&ctx_fdm)); // -- Operator CeedCall(CeedOperatorCreate(ceed_parent, qf_fdm, NULL, NULL, fdm_inv)); CeedCall(CeedOperatorSetField(*fdm_inv, "input", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); CeedCall(CeedOperatorSetField(*fdm_inv, "scale", rstr_qd_i, CEED_BASIS_NONE, q_data)); CeedCall(CeedOperatorSetField(*fdm_inv, "output", rstr, fdm_basis, CEED_VECTOR_ACTIVE)); // Cleanup CeedCall(CeedVectorDestroy(&q_data)); CeedCall(CeedBasisDestroy(&fdm_basis)); CeedCall(CeedElemRestrictionDestroy(&rstr_qd_i)); CeedCall(CeedQFunctionDestroy(&qf_fdm)); return CEED_ERROR_SUCCESS; } /// @}