From fd8bae8a685332de00e3c140b0b0e6a72d96d203 Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Thu, 9 Nov 2023 15:26:45 -0700 Subject: [PATCH] frame in OperatorApplyAtPoints in backend, utterly untested --- backends/ref/ceed-ref-operator.c | 482 +++++++++++++++++++++++++++++++ backends/ref/ceed-ref.c | 1 + backends/ref/ceed-ref.h | 3 +- 3 files changed, 485 insertions(+), 1 deletion(-) diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index 0b92d6fb4a..e4c5cff9fc 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -547,6 +547,473 @@ static int CeedOperatorLinearAssembleQFunctionUpdate_Ref(CeedOperator op, CeedVe return CeedOperatorLinearAssembleQFunctionCore_Ref(op, false, &assembled, &rstr, request); } +//------------------------------------------------------------------------------ +// Setup Input/Output Fields +//------------------------------------------------------------------------------ +static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op, bool is_input, CeedElemRestriction *block_rstr, + CeedVector *e_vecs_full, CeedVector *e_vecs, CeedVector *q_vecs, CeedInt start_e, CeedInt num_fields, + CeedInt Q) { + Ceed ceed; + CeedSize e_size, q_size; + CeedInt max_num_points, num_comp, size, P; + CeedQFunctionField *qf_fields; + CeedOperatorField *op_fields; + + CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); + if (is_input) { + CeedCallBackend(CeedOperatorGetFields(op, NULL, &op_fields, NULL, NULL)); + CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_fields, NULL, NULL)); + } else { + CeedCallBackend(CeedOperatorGetFields(op, NULL, NULL, NULL, &op_fields)); + CeedCallBackend(CeedQFunctionGetFields(qf, NULL, NULL, NULL, &qf_fields)); + } + + // Get max number of points + { + CeedInt dim; + CeedElemRestriction rstr_points; + CeedOperator_Ref *impl; + + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); + CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr_points, &max_num_points)); + CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr_points, &dim)); + CeedCallBackend(CeedOperatorGetData(op, &impl)); + CeedCallBackend(CeedVectorCreate(ceed, dim * max_num_points, &impl->point_coords_elem)); + } + + // Loop over fields + for (CeedInt i = 0; i < num_fields; i++) { + CeedEvalMode eval_mode; + CeedBasis basis; + + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_fields[i], &eval_mode)); + if (eval_mode != CEED_EVAL_WEIGHT) { + Ceed ceed_rstr; + CeedSize l_size; + CeedInt num_elem, elem_size, comp_stride; + CeedRestrictionType rstr_type; + CeedElemRestriction rstr; + + CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); + CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); + if (rstr_type == CEED_RESTRICTION_POINTS) { + // -- Only vector needed if field is at points + CeedCallBackend(CeedElemRestrictionReferenceCopy(rstr, &block_rstr[i + start_e])); + CeedCallBackend(CeedVectorCreate(ceed, num_comp * max_num_points, &e_vecs_full[i + start_e])); + } else { + // -- FEM fields need blocked restriction + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_fields[i], &rstr)); + CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed_rstr)); + CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem)); + CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size)); + CeedCallBackend(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); + CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); + + switch (rstr_type) { + case CEED_RESTRICTION_STANDARD: { + const CeedInt *offsets = NULL; + + CeedCallBackend(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); + CeedCallBackend(CeedElemRestrictionCreateBlocked(ceed_rstr, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, CEED_MEM_HOST, + CEED_COPY_VALUES, offsets, &block_rstr[i + start_e])); + CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); + } break; + case CEED_RESTRICTION_ORIENTED: { + const bool *orients = NULL; + const CeedInt *offsets = NULL; + + CeedCallBackend(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); + CeedCallBackend(CeedElemRestrictionGetOrientations(rstr, CEED_MEM_HOST, &orients)); + CeedCallBackend(CeedElemRestrictionCreateBlockedOriented(ceed_rstr, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, + CEED_MEM_HOST, CEED_COPY_VALUES, offsets, orients, &block_rstr[i + start_e])); + CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); + CeedCallBackend(CeedElemRestrictionRestoreOrientations(rstr, &orients)); + } break; + case CEED_RESTRICTION_CURL_ORIENTED: { + const CeedInt8 *curl_orients = NULL; + const CeedInt *offsets = NULL; + + CeedCallBackend(CeedElemRestrictionGetOffsets(rstr, CEED_MEM_HOST, &offsets)); + CeedCallBackend(CeedElemRestrictionGetCurlOrientations(rstr, CEED_MEM_HOST, &curl_orients)); + CeedCallBackend(CeedElemRestrictionCreateBlockedCurlOriented(ceed_rstr, num_elem, elem_size, block_size, num_comp, comp_stride, l_size, + CEED_MEM_HOST, CEED_COPY_VALUES, offsets, curl_orients, + &block_rstr[i + start_e])); + CeedCallBackend(CeedElemRestrictionRestoreOffsets(rstr, &offsets)); + CeedCallBackend(CeedElemRestrictionRestoreCurlOrientations(rstr, &curl_orients)); + } break; + case CEED_RESTRICTION_STRIDED: { + CeedInt strides[3]; + + CeedCallBackend(CeedElemRestrictionGetStrides(rstr, &strides)); + CeedCallBackend(CeedElemRestrictionCreateBlockedStrided(ceed_rstr, num_elem, elem_size, block_size, num_comp, l_size, strides, + &block_rstr[i + start_e])); + } break; + case CEED_RESTRICTION_POINTS: + // Empty case - won't occur + break; + } + CeedCallBackend(CeedElemRestrictionCreateVector(block_rstr[i + start_e], NULL, &e_vecs_full[i + start_e])); + } + } + + switch (eval_mode) { + case CEED_EVAL_NONE: + CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); + e_size = (CeedSize)max_num_points * size; + CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); + q_size = (CeedSize)max_num_points * size; + CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); + break; + case CEED_EVAL_INTERP: + case CEED_EVAL_GRAD: + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + CeedCallBackend(CeedOperatorFieldGetBasis(op_fields[i], &basis)); + CeedCallBackend(CeedQFunctionFieldGetSize(qf_fields[i], &size)); + CeedCallBackend(CeedBasisGetNumNodes(basis, &P)); + CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); + e_size = (CeedSize)P * num_comp; + CeedCallBackend(CeedVectorCreate(ceed, e_size, &e_vecs[i])); + q_size = (CeedSize)max_num_points * size; + CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); + break; + case CEED_EVAL_WEIGHT: + // Not used at this time + break; + } + if (is_input && e_vecs[i]) { + CeedCallBackend(CeedVectorSetArray(e_vecs[i], CEED_MEM_HOST, CEED_COPY_VALUES, NULL)); + } + } + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ +// Setup Operator +//------------------------------------------------------------------------------ +static int CeedOperatorSetupAtPoints_Ref(CeedOperator op) { + bool is_setup_done; + Ceed ceed; + Ceed_Ref *ceed_impl; + CeedInt Q, num_input_fields, num_output_fields; + CeedQFunctionField *qf_input_fields, *qf_output_fields; + CeedQFunction qf; + CeedOperatorField *op_input_fields, *op_output_fields; + CeedOperator_Ref *impl; + + CeedCallBackend(CeedOperatorIsSetupDone(op, &is_setup_done)); + if (is_setup_done) return CEED_ERROR_SUCCESS; + + CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); + CeedCallBackend(CeedGetData(ceed, &ceed_impl)); + CeedCallBackend(CeedOperatorGetData(op, &impl)); + CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); + CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); + CeedCallBackend(CeedQFunctionIsIdentity(qf, &impl->is_identity_qf)); + CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); + CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); + + // Allocate + CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->block_rstr)); + CeedCallBackend(CeedCalloc(num_input_fields + num_output_fields, &impl->e_vecs_full)); + + CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->input_states)); + CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_in)); + CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->e_vecs_out)); + CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_in)); + CeedCallBackend(CeedCalloc(CEED_FIELD_MAX, &impl->q_vecs_out)); + + impl->num_inputs = num_input_fields; + impl->num_outputs = num_output_fields; + + // Set up infield and outfield pointer arrays + // Infields + CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, true, impl->block_rstr, impl->e_vecs_full, impl->e_vecs_in, impl->q_vecs_in, 0, + num_input_fields, Q)); + // Outfields + CeedCallBackend(CeedOperatorSetupFieldsAtPoints_Ref(qf, op, false, impl->block_rstr, impl->e_vecs_full, impl->e_vecs_out, impl->q_vecs_out, + num_input_fields, num_output_fields, Q)); + + // Identity QFunctions + if (impl->is_identity_qf) { + CeedEvalMode in_mode, out_mode; + CeedQFunctionField *in_fields, *out_fields; + + CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &in_fields, NULL, &out_fields)); + CeedCallBackend(CeedQFunctionFieldGetEvalMode(in_fields[0], &in_mode)); + CeedCallBackend(CeedQFunctionFieldGetEvalMode(out_fields[0], &out_mode)); + + if (in_mode == CEED_EVAL_NONE && out_mode == CEED_EVAL_NONE) { + impl->is_identity_rstr_op = true; + } else { + CeedCallBackend(CeedVectorReferenceCopy(impl->q_vecs_in[0], &impl->q_vecs_out[0])); + } + } + + CeedCallBackend(CeedOperatorSetSetupDone(op)); + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ +// Setup Input Fields +//------------------------------------------------------------------------------ +static inline int CeedOperatorSetupInputsAtPoints_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, + CeedOperatorField *op_input_fields, CeedVector in_vec, CeedScalar *e_data[2 * CEED_FIELD_MAX], + CeedOperator_Ref *impl, CeedRequest *request) { + for (CeedInt i = 0; i < num_input_fields; i++) { + uint64_t state; + CeedEvalMode eval_mode; + CeedVector vec; + + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); + if (eval_mode == CEED_EVAL_WEIGHT) { // Skip + } else { + // Get input vector + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + if (vec != CEED_VECTOR_ACTIVE) { + // Restrict + CeedCallBackend(CeedVectorGetState(vec, &state)); + if (state != impl->input_states[i]) { + CeedCallBackend(CeedElemRestrictionApply(impl->block_rstr[i], CEED_NOTRANSPOSE, vec, impl->e_vecs_full[i], request)); + impl->input_states[i] = state; + } + // Get evec + CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_full[i], CEED_MEM_HOST, (const CeedScalar **)&e_data[i])); + } else { + // Set Qvec for CEED_EVAL_NONE + if (eval_mode == CEED_EVAL_NONE) { + CeedCallBackend(CeedVectorGetArrayRead(impl->e_vecs_in[i], CEED_MEM_HOST, (const CeedScalar **)&e_data[i])); + CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, e_data[i])); + CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_in[i], (const CeedScalar **)&e_data[i])); + } + } + } + } + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ +// Input Basis Action +//------------------------------------------------------------------------------ +static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_input_fields, + CeedOperatorField *op_input_fields, CeedInt num_input_fields, CeedVector in_vec, + CeedVector point_coords_elem, bool skip_active, CeedScalar *e_data[2 * CEED_FIELD_MAX], + CeedOperator_Ref *impl, CeedRequest *request) { + for (CeedInt i = 0; i < num_input_fields; i++) { + bool is_active_input = false; + CeedInt elem_size, size, num_comp; + CeedEvalMode eval_mode; + CeedVector vec; + CeedRestrictionType rstr_type; + CeedElemRestriction elem_rstr; + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + // Skip active input + is_active_input = vec == CEED_VECTOR_ACTIVE; + if (skip_active && is_active_input) continue; + + // Get elem_size, eval_mode, size + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); + CeedCallBackend(CeedElemRestrictionGetType(impl->block_rstr[i], &rstr_type)); + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); + CeedCallBackend(CeedQFunctionFieldGetSize(qf_input_fields[i], &size)); + // Restrict block active input + if (is_active_input) { + if (elem_type == CEED_RESTRICTION_POINTS) { + CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(impl->block_rstr[i], e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)) + } else { + CeedCallBackend(CeedElemRestrictionApplyBlock(impl->block_rstr[i], e, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[i], request)); + } + } + // Basis action + switch (eval_mode) { + case CEED_EVAL_NONE: + if (!is_active_input) { + CeedCallBackend(CeedVectorSetArray(impl->q_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][num_points_offset * size])); + } + break; + // Note - these basis eval modes require FEM fields + case CEED_EVAL_INTERP: + case CEED_EVAL_GRAD: + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); + if (!is_active_input) { + CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); + CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][e * elem_size * num_comp])); + } + CeedCallBackend( + CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); + break; + case CEED_EVAL_WEIGHT: + break; // No action + } + } + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ +// Output Basis Action +//------------------------------------------------------------------------------ +static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_points_offset, CeedInt num_points, CeedQFunctionField *qf_output_fields, + CeedOperatorField *op_output_fields, CeedInt num_input_fields, CeedInt num_output_fields, + CeedOperator op, CeedVector out_vec, CeedVector point_coords_elem, CeedOperator_Ref *impl, + CeedRequest *request) { + for (CeedInt i = 0; i < num_output_fields; i++) { + CeedEvalMode eval_mode; + CeedVector vec; + CeedRestrictionType rstr_type; + CeedElemRestriction elem_rstr; + CeedBasis basis; + + // Get elem_size, eval_mode, size + CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); + // Basis action + switch (eval_mode) { + case CEED_EVAL_NONE: + break; // No action + case CEED_EVAL_INTERP: + case CEED_EVAL_GRAD: + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); + CeedCallBackend( + CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); + break; + // LCOV_EXCL_START + case CEED_EVAL_WEIGHT: { + Ceed ceed; + CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); + return CeedError(ceed, CEED_ERROR_BACKEND, "CEED_EVAL_WEIGHT cannot be an output evaluation mode"); + // LCOV_EXCL_STOP + } + } + // Restrict output block + // Get output vector + CeedCallBackend(CeedElemRestrictionGetType(impl->block_rstr[i + impl->num_inputs], &rstr_type)); + CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); + if (vec == CEED_VECTOR_ACTIVE) vec = out_vec; + // Restrict + if (elem_type == CEED_RESTRICTION_POINTS) { + CeedCallBackend( + CeedElemRestrictionApplyAtPointsInElement(impl->block_rstr[i + impl->num_inputs], e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)) + } else { + CeedCallBackend(CeedElemRestrictionApplyBlock(impl->block_rstr[i + impl->num_inputs], e, CEED_TRANSPOSE, impl->e_vecs_out[i], vec, request)); + } + } + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ +// Restore Input Vectors +//------------------------------------------------------------------------------ +static inline int CeedOperatorRestoreInputsAtPoints_Ref(CeedInt num_input_fields, CeedQFunctionField *qf_input_fields, + CeedOperatorField *op_input_fields, CeedScalar *e_data[2 * CEED_FIELD_MAX], + CeedOperator_Ref *impl) { + for (CeedInt i = 0; i < num_input_fields; i++) { + CeedEvalMode eval_mode; + CeedVector vec; + + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); + CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); + if (eval_mode != CEED_EVAL_WEIGHT && vec != CEED_VECTOR_ACTIVE) { + CeedCallBackend(CeedVectorRestoreArrayRead(impl->e_vecs_full[i], (const CeedScalar **)&e_data[i])); + } + } + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ +// Operator Apply +//------------------------------------------------------------------------------ +static int CeedOperatorApplyAddAtPoints_Ref(CeedOperator op, CeedVector in_vec, CeedVector out_vec, CeedRequest *request) { + Ceed ceed; + Ceed_Ref *ceed_impl; + CeedInt num_points_offset = 0, num_input_fields, num_output_fields, num_elem; + CeedEvalMode eval_mode; + CeedScalar *e_data[2 * CEED_FIELD_MAX] = {0}; + CeedVector point_coords = NULL; + CeedElemRestriction rstr_points = NULL; + CeedQFunctionField *qf_input_fields, *qf_output_fields; + CeedQFunction qf; + CeedOperatorField *op_input_fields, *op_output_fields; + CeedOperator_Ref *impl; + + CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); + CeedCallBackend(CeedGetData(ceed, &ceed_impl)); + CeedCallBackend(CeedOperatorGetData(op, &impl)); + CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); + CeedCallBackend(CeedOperatorGetNumQuadraturePoints(op, &Q)); + CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); + CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); + CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); + + // Setup + CeedCallBackend(CeedOperatorSetupAtPoints_Ref(op)); + + // Restriction only operator + if (impl->is_identity_rstr_op) { + for (CeedInt b = 0; b < num_blocks; b++) { + CeedCallBackend(CeedElemRestrictionApplyBlock(impl->block_rstr[0], b, CEED_NOTRANSPOSE, in_vec, impl->e_vecs_in[0], request)); + CeedCallBackend(CeedElemRestrictionApplyBlock(impl->block_rstr[1], b, CEED_TRANSPOSE, impl->e_vecs_in[0], out_vec, request)); + } + return CEED_ERROR_SUCCESS; + } + + // Point coordinates + CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, &point_coords)); + + // Input Evecs and Restriction + CeedCallBackend(CeedOperatorSetupInputsAtPoints_Ref(num_input_fields, qf_input_fields, op_input_fields, in_vec, e_data, impl, request)); + + // Output Lvecs, Evecs, and Qvecs + for (CeedInt i = 0; i < num_output_fields; i++) { + // Set Qvec if needed + CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); + if (eval_mode == CEED_EVAL_NONE) { + // Set qvec to single block evec + CeedCallBackend(CeedVectorGetArrayWrite(impl->e_vecs_out[i], CEED_MEM_HOST, &e_data[i + num_input_fields])); + CeedCallBackend(CeedVectorSetArray(impl->q_vecs_out[i], CEED_MEM_HOST, CEED_USE_POINTER, e_data[i + num_input_fields])); + CeedCallBackend(CeedVectorRestoreArray(impl->e_vecs_out[i], &e_data[i + num_input_fields])); + } + } + + // Loop through elements + for (CeedInt e = 0; e < num_elem; e++) { + CeedInt num_points; + + // Setup points for element + CeedCallBackend(CeedElemRestrictionApplyAtPointsInElement(rstr_points, e, CEED_NOTRANSPOSE, point_coords, impl->point_coords_elem, request)); + CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points)); + + // Input basis apply + CeedCallBackend(CeedOperatorInputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_input_fields, op_input_fields, num_input_fields, in_vec, + false, e_data, impl, request)); + + // Q function + if (!impl->is_identity_qf) { + CeedCallBackend(CeedQFunctionApply(qf, num_points, impl->q_vecs_in, impl->q_vecs_out)); + } + + // Output basis apply and restriction + CeedCallBackend(CeedOperatorOutputBasisAtPoints_Ref(e, num_points_offset, num_points, qf_output_fields, op_output_fields, num_input_fields, + num_output_fields, op, out_vec, impl, request)); + + num_points_offset += num_points; + } + + // Restore input arrays + CeedCallBackend(CeedOperatorRestoreInputsAtPoints_Ref(num_input_fields, qf_input_fields, op_input_fields, e_data, impl)); + + // Cleanup point coordinates + CeedCallBackend(CeedVectorDestroy(&point_coords)); + CeedCallBackend(CeedElemRestriction(&rstr_points)); + return CEED_ERROR_SUCCESS; +} + //------------------------------------------------------------------------------ // Operator Destroy //------------------------------------------------------------------------------ @@ -573,6 +1040,7 @@ static int CeedOperatorDestroy_Ref(CeedOperator op) { } CeedCallBackend(CeedFree(&impl->e_vecs_out)); CeedCallBackend(CeedFree(&impl->q_vecs_out)); + CeedCallBackend(CeedVectorDestroy(&impl->point_coords_elem)); // QFunction assembly for (CeedInt i = 0; i < impl->num_active_in; i++) { @@ -602,3 +1070,17 @@ int CeedOperatorCreate_Ref(CeedOperator op) { } //------------------------------------------------------------------------------ +// Operator Create At Points +//------------------------------------------------------------------------------ +int CeedOperatorCreateAtPoints_Ref(CeedOperator op) { + Ceed ceed; + CeedOperator_Ref *impl; + + CeedCallBackend(CeedCalloc(1, &impl)); + CeedCallBackend(CeedOperatorSetData(op, impl)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAddAtPoints_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Ref)); + return CEED_ERROR_SUCCESS; +} + +//------------------------------------------------------------------------------ diff --git a/backends/ref/ceed-ref.c b/backends/ref/ceed-ref.c index 11e4137261..088ef3a2b2 100644 --- a/backends/ref/ceed-ref.c +++ b/backends/ref/ceed-ref.c @@ -31,6 +31,7 @@ static int CeedInit_Ref(const char *resource, Ceed ceed) { CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "QFunctionCreate", CeedQFunctionCreate_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "QFunctionContextCreate", CeedQFunctionContextCreate_Ref)); CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "OperatorCreate", CeedOperatorCreate_Ref)); + CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "OperatorCreateAtPoints", CeedOperatorCreateAtPoints_Ref)); return CEED_ERROR_SUCCESS; } diff --git a/backends/ref/ceed-ref.h b/backends/ref/ceed-ref.h index 36e581211a..bfd5445771 100644 --- a/backends/ref/ceed-ref.h +++ b/backends/ref/ceed-ref.h @@ -56,7 +56,7 @@ typedef struct { CeedVector *q_vecs_out; /* Single element output Q-vectors */ CeedInt num_inputs, num_outputs; CeedInt num_active_in, num_active_out; - CeedVector *qf_active_in; + CeedVector *qf_active_in, point_coords_elem; } CeedOperator_Ref; CEED_INTERN int CeedVectorCreate_Ref(CeedSize n, CeedVector vec); @@ -80,5 +80,6 @@ CEED_INTERN int CeedQFunctionCreate_Ref(CeedQFunction qf); CEED_INTERN int CeedQFunctionContextCreate_Ref(CeedQFunctionContext ctx); CEED_INTERN int CeedOperatorCreate_Ref(CeedOperator op); +CEED_INTERN int CeedOperatorCreateAtPoints_Ref(CeedOperator op); #endif // CEED_REF_H