From fc0f7cc68128f3536a834f19d72828f4c59a4439 Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Fri, 31 May 2024 11:23:11 -0600 Subject: [PATCH 1/2] basis - ApplyAtPoints should take number of elem --- backends/ref/ceed-ref-operator.c | 10 ++-- include/ceed-impl.h | 2 +- include/ceed/ceed.h | 34 ++++++------ interface/ceed-basis.c | 55 ++++++++++--------- .../src/generated/libceed_bindings.jl | 4 +- tests/t350-basis.c | 2 +- tests/t351-basis.c | 2 +- tests/t352-basis.c | 2 +- tests/t353-basis.c | 7 ++- tests/t354-basis.c | 7 ++- tests/t355-basis.c | 2 +- tests/t356-basis.c | 2 +- tests/t357-basis.c | 4 +- 13 files changed, 70 insertions(+), 63 deletions(-) diff --git a/backends/ref/ceed-ref-operator.c b/backends/ref/ceed-ref-operator.c index a0a6a29a2d..538accec17 100644 --- a/backends/ref/ceed-ref-operator.c +++ b/backends/ref/ceed-ref-operator.c @@ -654,7 +654,7 @@ static int CeedOperatorSetupFieldsAtPoints_Ref(CeedQFunction qf, CeedOperator op q_size = (CeedSize)max_num_points; CeedCallBackend(CeedVectorCreate(ceed, q_size, &q_vecs[i])); CeedCallBackend( - CeedBasisApplyAtPoints(basis, max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); + CeedBasisApplyAtPoints(basis, 1, &max_num_points, CEED_NOTRANSPOSE, CEED_EVAL_WEIGHT, CEED_VECTOR_NONE, CEED_VECTOR_NONE, q_vecs[i])); break; } // Initialize full arrays for E-vectors and Q-vectors @@ -767,7 +767,7 @@ static inline int CeedOperatorInputBasisAtPoints_Ref(CeedInt e, CeedInt num_poin CeedCallBackend(CeedVectorSetArray(impl->e_vecs_in[i], CEED_MEM_HOST, CEED_USE_POINTER, &e_data[i][(CeedSize)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])); + CeedBasisApplyAtPoints(basis, 1, &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 @@ -803,7 +803,7 @@ static inline int CeedOperatorOutputBasisAtPoints_Ref(CeedInt e, CeedInt num_poi 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])); + CeedBasisApplyAtPoints(basis, 1, &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: { @@ -1237,7 +1237,7 @@ static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, Ce case CEED_EVAL_DIV: case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); - CeedCallBackend(CeedBasisApplyAtPoints(basis, num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], + CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_NOTRANSPOSE, eval_mode, impl->point_coords_elem, impl->e_vecs_in[i], impl->q_vecs_in[i])); break; case CEED_EVAL_WEIGHT: @@ -1280,7 +1280,7 @@ static int CeedOperatorLinearAssembleAddDiagonalAtPoints_Ref(CeedOperator op, Ce case CEED_EVAL_DIV: case CEED_EVAL_CURL: CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); - CeedCallBackend(CeedBasisApplyAtPoints(basis, num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], + CeedCallBackend(CeedBasisApplyAtPoints(basis, 1, &num_points, CEED_TRANSPOSE, eval_mode, impl->point_coords_elem, impl->q_vecs_out[i], impl->e_vecs_out[i])); break; // LCOV_EXCL_START diff --git a/include/ceed-impl.h b/include/ceed-impl.h index 4fe0c67c9f..956b30b757 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -184,7 +184,7 @@ struct CeedElemRestriction_private { struct CeedBasis_private { Ceed ceed; int (*Apply)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector); - int (*ApplyAtPoints)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector, CeedVector); + int (*ApplyAtPoints)(CeedBasis, CeedInt, const CeedInt *, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector, CeedVector); int (*Destroy)(CeedBasis); int ref_count; bool is_tensor_basis; /* flag for tensor basis */ diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index 9ebb40534d..f0f6f65be7 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -289,23 +289,23 @@ CEED_EXTERN int CeedElemRestrictionDestroy(CeedElemRestriction *rstr); // \int_\Omega v^T f_0(u, \nabla u, qdata) + (\nabla v)^T f_1(u, \nabla u, qdata) // where gradients are with respect to the reference element. -CEED_EXTERN int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, - CeedBasis *basis); -CEED_EXTERN int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, - const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis *basis); -CEED_EXTERN int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, - const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); -CEED_EXTERN int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, - const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); -CEED_EXTERN int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, - const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); -CEED_EXTERN int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasis *basis_project); -CEED_EXTERN int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy); -CEED_EXTERN int CeedBasisView(CeedBasis basis, FILE *stream); -CEED_EXTERN int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v); -CEED_EXTERN int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, - CeedVector u, CeedVector v); -CEED_EXTERN int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed); +CEED_EXTERN int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, + CeedBasis *basis); +CEED_EXTERN int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, + const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis *basis); +CEED_EXTERN int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, + const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); +CEED_EXTERN int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, + const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); +CEED_EXTERN int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt nqpts, const CeedScalar *interp, + const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weights, CeedBasis *basis); +CEED_EXTERN int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasis *basis_project); +CEED_EXTERN int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy); +CEED_EXTERN int CeedBasisView(CeedBasis basis, FILE *stream); +CEED_EXTERN int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v); +CEED_EXTERN int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, + CeedVector x_ref, CeedVector u, CeedVector v); +CEED_EXTERN int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed); CEED_EXTERN Ceed CeedBasisReturnCeed(CeedBasis basis); CEED_EXTERN int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim); CEED_EXTERN int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo); diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index ed24d72c66..f64ac7bf50 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -1499,7 +1499,9 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, @brief Apply basis evaluation from nodes to arbitrary points @param[in] basis `CeedBasis` to evaluate - @param[in] num_points The number of points to apply the basis evaluation to + @param[in] num_elem The number of elements to apply the basis evaluation to; + the backend will specify the ordering in @ref CeedElemRestrictionCreate() + @param[in] num_points The number of points to apply the basis evaluation to in each element @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, @@ -1513,10 +1515,10 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, @ref User **/ -int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, - CeedVector v) { +int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, + CeedVector x_ref, CeedVector u, CeedVector v) { bool is_tensor_basis; - CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1; + CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1, total_num_points = 0; CeedSize x_length = 0, u_length = 0, v_length; Ceed ceed; @@ -1532,12 +1534,13 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod if (u != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(u, &u_length)); // Check compatibility of topological and geometrical dimensions + for (CeedInt i = 0; i < num_elem; i++) total_num_points += num_points[i]; CeedCheck((t_mode == CEED_TRANSPOSE && v_length % num_nodes == 0) || (t_mode == CEED_NOTRANSPOSE && u_length % num_nodes == 0) || (eval_mode == CEED_EVAL_WEIGHT), ceed, CEED_ERROR_DIMENSION, "Length of input/output vectors incompatible with basis dimensions and number of points"); // Check compatibility coordinates vector - CeedCheck((x_length >= num_points * dim) || (eval_mode == CEED_EVAL_WEIGHT), ceed, CEED_ERROR_DIMENSION, + CeedCheck((x_length >= total_num_points * dim) || (eval_mode == CEED_EVAL_WEIGHT), ceed, CEED_ERROR_DIMENSION, "Length of reference coordinate vector incompatible with basis dimension and number of points"); // Check CEED_EVAL_WEIGHT only on CEED_NOTRANSPOSE @@ -1548,15 +1551,16 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod bool has_good_dims = true; switch (eval_mode) { case CEED_EVAL_INTERP: - has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= num_points * num_q_comp || v_length >= num_nodes * num_comp)) || - (t_mode == CEED_NOTRANSPOSE && (v_length >= num_points * num_q_comp || u_length >= num_nodes * num_comp))); + has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= total_num_points * num_q_comp || v_length >= num_elem * num_nodes * num_comp)) || + (t_mode == CEED_NOTRANSPOSE && (v_length >= total_num_points * num_q_comp || u_length >= num_elem * num_nodes * num_comp))); break; case CEED_EVAL_GRAD: - has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= num_points * num_q_comp * dim || v_length >= num_nodes * num_comp)) || - (t_mode == CEED_NOTRANSPOSE && (v_length >= num_points * num_q_comp * dim || u_length >= num_nodes * num_comp))); + has_good_dims = + ((t_mode == CEED_TRANSPOSE && (u_length >= total_num_points * num_q_comp * dim || v_length >= num_elem * num_nodes * num_comp)) || + (t_mode == CEED_NOTRANSPOSE && (v_length >= total_num_points * num_q_comp * dim || u_length >= num_elem * num_nodes * num_comp))); break; case CEED_EVAL_WEIGHT: - has_good_dims = t_mode == CEED_NOTRANSPOSE && (v_length >= num_points); + has_good_dims = t_mode == CEED_NOTRANSPOSE && (v_length >= total_num_points); break; // LCOV_EXCL_START case CEED_EVAL_NONE: @@ -1569,13 +1573,14 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod // Backend method if (basis->ApplyAtPoints) { - CeedCall(basis->ApplyAtPoints(basis, num_points, t_mode, eval_mode, x_ref, u, v)); + CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); return CEED_ERROR_SUCCESS; } // Default implementation CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); CeedCheck(is_tensor_basis, ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points only supported for tensor product bases"); + CeedCheck(num_elem == 1, ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points only supported for a single element at a time"); if (eval_mode == CEED_EVAL_WEIGHT) { CeedCall(CeedVectorSetValue(v, 1.0)); return CEED_ERROR_SUCCESS; @@ -1652,18 +1657,18 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; // ---- Values at point - for (CeedInt p = 0; p < num_points; p++) { + for (CeedInt p = 0; p < total_num_points; p++) { CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; for (CeedInt d = 0; d < dim; d++) { // ------ Tensor contract with current Chebyshev polynomial values - CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * num_points + p], Q_1d, chebyshev_x)); + CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2])); pre /= Q_1d; post *= 1; } - for (CeedInt c = 0; c < num_comp; c++) v_array[c * num_points + p] = tmp[dim % 2][c]; + for (CeedInt c = 0; c < num_comp; c++) v_array[c * total_num_points + p] = tmp[dim % 2][c]; } break; } @@ -1671,21 +1676,21 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; // ---- Values at point - for (CeedInt p = 0; p < num_points; p++) { + for (CeedInt p = 0; p < total_num_points; p++) { // Dim**2 contractions, apply grad when pass == dim for (CeedInt pass = 0; pass < dim; pass++) { CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; for (CeedInt d = 0; d < dim; d++) { // ------ Tensor contract with current Chebyshev polynomial values - if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * num_points + p], Q_1d, chebyshev_x)); - else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * num_points + p], Q_1d, chebyshev_x)); + if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); + else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2])); pre /= Q_1d; post *= 1; } - for (CeedInt c = 0; c < num_comp; c++) v_array[(pass * num_comp + c) * num_points + p] = tmp[dim % 2][c]; + for (CeedInt c = 0; c < num_comp; c++) v_array[(pass * num_comp + c) * total_num_points + p] = tmp[dim % 2][c]; } } break; @@ -1715,13 +1720,13 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; // ---- Values at point - for (CeedInt p = 0; p < num_points; p++) { + for (CeedInt p = 0; p < total_num_points; p++) { CeedInt pre = num_comp * 1, post = 1; - for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[c * num_points + p]; + for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[c * total_num_points + p]; for (CeedInt d = 0; d < dim; d++) { // ------ Tensor contract with current Chebyshev polynomial values - CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * num_points + p], Q_1d, chebyshev_x)); + CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == (dim - 1), tmp[d % 2], d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); pre /= 1; @@ -1734,16 +1739,16 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; // ---- Values at point - for (CeedInt p = 0; p < num_points; p++) { + for (CeedInt p = 0; p < total_num_points; p++) { // Dim**2 contractions, apply grad when pass == dim for (CeedInt pass = 0; pass < dim; pass++) { CeedInt pre = num_comp * 1, post = 1; - for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[(pass * num_comp + c) * num_points + p]; + for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[(pass * num_comp + c) * total_num_points + p]; for (CeedInt d = 0; d < dim; d++) { // ------ Tensor contract with current Chebyshev polynomial values - if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * num_points + p], Q_1d, chebyshev_x)); - else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * num_points + p], Q_1d, chebyshev_x)); + if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); + else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, (p > 0 || (p == 0 && pass > 0)) && d == (dim - 1), tmp[d % 2], d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); diff --git a/julia/LibCEED.jl/src/generated/libceed_bindings.jl b/julia/LibCEED.jl/src/generated/libceed_bindings.jl index f814609b86..9cbf889dd9 100644 --- a/julia/LibCEED.jl/src/generated/libceed_bindings.jl +++ b/julia/LibCEED.jl/src/generated/libceed_bindings.jl @@ -436,8 +436,8 @@ function CeedBasisApply(basis, num_elem, t_mode, eval_mode, u, v) ccall((:CeedBasisApply, libceed), Cint, (CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector), basis, num_elem, t_mode, eval_mode, u, v) end -function CeedBasisApplyAtPoints(basis, num_points, t_mode, eval_mode, x_ref, u, v) - ccall((:CeedBasisApplyAtPoints, libceed), Cint, (CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector, CeedVector), basis, num_points, t_mode, eval_mode, x_ref, u, v) +function CeedBasisApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v) + ccall((:CeedBasisApplyAtPoints, libceed), Cint, (CeedBasis, CeedInt, Ptr{CeedInt}, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector, CeedVector), basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v) end function CeedBasisGetCeed(basis, ceed) diff --git a/tests/t350-basis.c b/tests/t350-basis.c index 54979bb9a0..becc0d98ea 100644 --- a/tests/t350-basis.c +++ b/tests/t350-basis.c @@ -56,7 +56,7 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); { const CeedScalar *x_array, *v_array; diff --git a/tests/t351-basis.c b/tests/t351-basis.c index 14b23730e1..84f59cc838 100644 --- a/tests/t351-basis.c +++ b/tests/t351-basis.c @@ -65,7 +65,7 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); { const CeedScalar *x_array, *v_array; diff --git a/tests/t352-basis.c b/tests/t352-basis.c index a4bf13d8b6..c2da0e2dd4 100644 --- a/tests/t352-basis.c +++ b/tests/t352-basis.c @@ -65,7 +65,7 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); { const CeedScalar *x_array, *v_array; diff --git a/tests/t353-basis.c b/tests/t353-basis.c index 22f80ddcdd..83fd16adb0 100644 --- a/tests/t353-basis.c +++ b/tests/t353-basis.c @@ -60,17 +60,18 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); for (CeedInt i = 0; i < num_points; i++) { - CeedScalar fx = 0.0; + const CeedInt num_point[1] = {1}; + CeedScalar fx = 0.0; const CeedScalar *x_array, *u_array, *v_array, *u_point_array; CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); CeedVectorSetValue(x_point, x_array[i]); - CeedBasisApplyAtPoints(basis_u, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point); + CeedBasisApplyAtPoints(basis_u, 1, num_point, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point); CeedVectorGetArrayRead(u_point, CEED_MEM_HOST, &u_point_array); for (CeedInt j = 0; j < p; j++) fx += u_array[j] * u_point_array[j]; if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) printf("%f != %f = f(%f)\n", v_array[i], fx, x_array[i]); diff --git a/tests/t354-basis.c b/tests/t354-basis.c index 85f0ac2293..e137f8c44f 100644 --- a/tests/t354-basis.c +++ b/tests/t354-basis.c @@ -69,10 +69,11 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); for (CeedInt i = 0; i < num_points; i++) { - CeedScalar fx = 0.0; + const CeedInt num_point[1] = {1}; + CeedScalar fx = 0.0; CeedScalar coord[dim]; const CeedScalar *x_array, *u_array, *v_array, *u_point_array; @@ -81,7 +82,7 @@ int main(int argc, char **argv) { CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * num_points + i]; CeedVectorSetArray(x_point, CEED_MEM_HOST, CEED_COPY_VALUES, coord); - CeedBasisApplyAtPoints(basis_u, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point); + CeedBasisApplyAtPoints(basis_u, 1, num_point, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point); CeedVectorGetArrayRead(u_point, CEED_MEM_HOST, &u_point_array); for (CeedInt j = 0; j < p_dim; j++) fx += u_array[j] * u_point_array[j]; if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) { diff --git a/tests/t355-basis.c b/tests/t355-basis.c index 7fd7906dcb..5b93764a7a 100644 --- a/tests/t355-basis.c +++ b/tests/t355-basis.c @@ -62,7 +62,7 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v); { const CeedScalar *x_array, *v_array; diff --git a/tests/t356-basis.c b/tests/t356-basis.c index 8eb3c57e7c..263cc43b66 100644 --- a/tests/t356-basis.c +++ b/tests/t356-basis.c @@ -75,7 +75,7 @@ int main(int argc, char **argv) { CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, v); { const CeedScalar *x_array, *v_array; diff --git a/tests/t357-basis.c b/tests/t357-basis.c index ecfa56476c..0f4e105a66 100644 --- a/tests/t357-basis.c +++ b/tests/t357-basis.c @@ -82,8 +82,8 @@ int main(int argc, char **argv) { } // Calculate G u at arbitrary points, G' * 1 at dofs - CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, u_points); - CeedBasisApplyAtPoints(basis_u, num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_GRAD, x_points, u, u_points); + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v); { const CeedScalar *u_array, *v_array, *u_points_array; From faed4840b5126336d09076d80e315366da92cc18 Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Thu, 6 Jun 2024 09:27:14 -0600 Subject: [PATCH 2/2] doc - fix docstring alignment Co-authored-by: James Wright --- interface/ceed-basis.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index f64ac7bf50..6acb3fdabc 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -1499,9 +1499,9 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, @brief Apply basis evaluation from nodes to arbitrary points @param[in] basis `CeedBasis` to evaluate - @param[in] num_elem The number of elements to apply the basis evaluation to; + @param[in] num_elem The number of elements to apply the basis evaluation to; the backend will specify the ordering in @ref CeedElemRestrictionCreate() - @param[in] num_points The number of points to apply the basis evaluation to in each element + @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values,