From 3b2de1f49242382f00cac6531aca5c6395e73bb2 Mon Sep 17 00:00:00 2001 From: Jeremy L Thompson Date: Mon, 12 Aug 2024 14:44:25 -0600 Subject: [PATCH] basis - add CeedBasisApplyAddAtPoints + default impl --- include/ceed-impl.h | 1 + include/ceed/ceed.h | 2 + interface/ceed-basis.c | 124 ++++++++++++++++++++++++++++++++++++----- interface/ceed.c | 1 + tests/t364-basis.c | 98 ++++++++++++++++++++++++++++++++ tests/t365-basis.c | 123 ++++++++++++++++++++++++++++++++++++++++ 6 files changed, 334 insertions(+), 15 deletions(-) create mode 100644 tests/t364-basis.c create mode 100644 tests/t365-basis.c diff --git a/include/ceed-impl.h b/include/ceed-impl.h index 889d8a5d31..0b76071a08 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -188,6 +188,7 @@ struct CeedBasis_private { int (*Apply)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector); int (*ApplyAdd)(CeedBasis, CeedInt, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector); int (*ApplyAtPoints)(CeedBasis, CeedInt, const CeedInt *, CeedTransposeMode, CeedEvalMode, CeedVector, CeedVector, CeedVector); + int (*ApplyAddAtPoints)(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 fb0d37a35d..d847525da8 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -308,6 +308,8 @@ CEED_EXTERN int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeM CEED_EXTERN int CeedBasisApplyAdd(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 CeedBasisApplyAddAtPoints(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); diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index afb62180b7..23e05781eb 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -1579,8 +1579,8 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, @param[in] basis `CeedBasis` to evaluate @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] t_mode @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes. - @ref CEED_NOTRANSPOSE is not valid for CeedBasisApplyAdd. + @param[in] t_mode @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes; + @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAdd()` @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, @ref CEED_EVAL_INTERP to use interpolated values, @ref CEED_EVAL_GRAD to use gradients, @@ -1603,7 +1603,7 @@ int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mod } /** - @brief Apply basis evaluation from nodes to arbitrary points + @brief Check input vector dimensions for CeedBasisApply[Add]AtPoints @param[in] basis `CeedBasis` to evaluate @param[in] num_elem The number of elements to apply the basis evaluation to; @@ -1620,11 +1620,10 @@ int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mod @return An error code: 0 - success, otherwise - failure - @ref User + @ref Developer **/ -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; +static int CeedBasisApplyAtPointsCheckDims(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, + CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 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; @@ -1677,16 +1676,48 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num // LCOV_EXCL_STOP } CeedCheck(has_good_dims, ceed, CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); + return CEED_ERROR_SUCCESS; +} - // Backend method - if (basis->ApplyAtPoints) { - CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); - return CEED_ERROR_SUCCESS; - } +/** + @brief Default implimentation to apply basis evaluation from nodes to arbitrary points + + @param[in] basis `CeedBasis` to evaluate + @param[in] apply_add Sum result into target vector or overwrite + @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 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, + @ref CEED_EVAL_GRAD to use gradients, + @ref CEED_EVAL_WEIGHT to use quadrature weights + @param[in] x_ref `CeedVector` holding reference coordinates of each point + @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE + @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP + + @return An error code: 0 - success, otherwise - failure + + @ref Developer +**/ +static int CeedBasisApplyAtPoints_Core(CeedBasis basis, bool apply_add, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, + CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { + CeedInt dim, num_comp, P_1d = 1, Q_1d = 1, total_num_points = num_points[0]; + Ceed ceed; + + CeedCall(CeedBasisGetCeed(basis, &ceed)); + CeedCall(CeedBasisGetDimension(basis, &dim)); + CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); + CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); + CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); // 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"); + { + bool is_tensor_basis; + + 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)); @@ -1858,13 +1889,76 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num CeedCall(CeedVectorRestoreArrayRead(u, &u_array)); // -- Interpolate transpose from Chebyshev coefficients - CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); + if (apply_add) CeedCall(CeedBasisApplyAdd(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); + else CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); break; } } return CEED_ERROR_SUCCESS; } +/** + @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; + the backend will specify the ordering in @ref CeedElemRestrictionCreate() + @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, + @ref CEED_EVAL_GRAD to use gradients, + @ref CEED_EVAL_WEIGHT to use quadrature weights + @param[in] x_ref `CeedVector` holding reference coordinates of each point + @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE + @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, + CeedVector x_ref, CeedVector u, CeedVector v) { + CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); + if (basis->ApplyAtPoints) { + CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); + } else { + CeedCall(CeedBasisApplyAtPoints_Core(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); + } + return CEED_ERROR_SUCCESS; +} + +/** + @brief Apply basis evaluation from nodes to arbitrary points and sum into target vector + + @param[in] basis `CeedBasis` to evaluate + @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 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_NOTRANSPOSE is not valid for `CeedBasisApplyAddAtPoints()` + @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, + @ref CEED_EVAL_GRAD to use gradients, + @ref CEED_EVAL_WEIGHT to use quadrature weights + @param[in] x_ref `CeedVector` holding reference coordinates of each point + @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE + @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedBasisApplyAddAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, + CeedVector x_ref, CeedVector u, CeedVector v) { + CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); + if (basis->ApplyAddAtPoints) { + CeedCall(basis->ApplyAddAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); + } else { + CeedCall(CeedBasisApplyAtPoints_Core(basis, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); + } + return CEED_ERROR_SUCCESS; +} + /** @brief Get the `Ceed` associated with a `CeedBasis` diff --git a/interface/ceed.c b/interface/ceed.c index 3a9be0c5dc..3bdd471454 100644 --- a/interface/ceed.c +++ b/interface/ceed.c @@ -954,6 +954,7 @@ int CeedInit(const char *resource, Ceed *ceed) { CEED_FTABLE_ENTRY(CeedBasis, Apply), CEED_FTABLE_ENTRY(CeedBasis, ApplyAdd), CEED_FTABLE_ENTRY(CeedBasis, ApplyAtPoints), + CEED_FTABLE_ENTRY(CeedBasis, ApplyAddAtPoints), CEED_FTABLE_ENTRY(CeedBasis, Destroy), CEED_FTABLE_ENTRY(CeedTensorContract, Apply), CEED_FTABLE_ENTRY(CeedTensorContract, Destroy), diff --git a/tests/t364-basis.c b/tests/t364-basis.c new file mode 100644 index 0000000000..6ab4058d30 --- /dev/null +++ b/tests/t364-basis.c @@ -0,0 +1,98 @@ +/// @file +/// Test polynomial interpolation transpose ApplyAdd from arbitrary points in 1D +/// \test Test polynomial interpolation transpose ApplyAdd from arbitrary points in 1D +#include +#include +#include + +#define ALEN(a) (sizeof(a) / sizeof((a)[0])) + +static CeedScalar Eval(CeedScalar x, CeedInt n, const CeedScalar *c) { + CeedScalar y = c[n - 1]; + for (CeedInt i = n - 2; i >= 0; i--) y = y * x + c[i]; + return y; +} + +int main(int argc, char **argv) { + Ceed ceed; + CeedVector x, x_nodes, x_points, x_point, u, v, u_point, v_point; + CeedBasis basis_x, basis_u; + const CeedInt p = 5, q = 5, num_points = 4; + const CeedScalar c[4] = {1, 2, 3, 4}; // 1 + 2x + 3x^2 + ... + + CeedInit(argv[1], &ceed); + + CeedVectorCreate(ceed, 2, &x); + CeedVectorCreate(ceed, p, &x_nodes); + CeedVectorCreate(ceed, num_points, &x_points); + CeedVectorCreate(ceed, 1, &x_point); + CeedVectorCreate(ceed, p, &u); + CeedVectorCreate(ceed, num_points, &v); + CeedVectorCreate(ceed, p, &u_point); + CeedVectorCreate(ceed, 1, &v_point); + CeedVectorSetValue(v_point, 1.0); + + // Get nodal coordinates + CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, p, CEED_GAUSS_LOBATTO, &basis_x); + { + CeedScalar x_array[2]; + + for (CeedInt i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1); + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes); + + // Set values of u at nodes + { + const CeedScalar *x_array; + CeedScalar u_array[p]; + + CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array); + for (CeedInt i = 0; i < p; i++) u_array[i] = Eval(x_array[i], ALEN(c), c); + CeedVectorRestoreArrayRead(x_nodes, &x_array); + CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); + } + + // Interpolate to arbitrary points + CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u); + { + CeedScalar x_array[4] = {-0.33, -0.65, 0.16, 0.99}; + + CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedBasisApplyAtPoints(basis_u, 1, &num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v); + + for (CeedInt i = 0; i < num_points; i++) { + 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, num_point, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point); + // Double it + CeedBasisApplyAddAtPoints(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] * 2.0 - fx) > 100. * CEED_EPSILON) printf("%f != %f = f(%f)\n", v_array[i] * 2.0, fx, x_array[i]); + CeedVectorRestoreArrayRead(u_point, &u_point_array); + CeedVectorRestoreArrayRead(x_points, &x_array); + CeedVectorRestoreArrayRead(u, &u_array); + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_nodes); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&x_point); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedVectorDestroy(&u_point); + CeedVectorDestroy(&v_point); + CeedBasisDestroy(&basis_x); + CeedBasisDestroy(&basis_u); + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t365-basis.c b/tests/t365-basis.c new file mode 100644 index 0000000000..74f93ce881 --- /dev/null +++ b/tests/t365-basis.c @@ -0,0 +1,123 @@ +/// @file +/// Test gradient transpose ApplyAdd in multiple dimensions at arbitrary points +/// \test Test gradient transpose ApplyAdd in multiple dimensions at arbitrary points +#include +#include +#include + +static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { + CeedScalar result = tanh(x[0] + 0.1); + if (dim > 1) result += atan(x[1] + 0.2); + if (dim > 2) result += exp(-(x[2] + 0.3) * (x[2] + 0.3)); + return result; +} + +static CeedScalar GetTolerance(CeedScalarType scalar_type, int dim) { + CeedScalar tol; + if (scalar_type == CEED_SCALAR_FP32) { + if (dim == 3) tol = 0.005; + else tol = 1.e-4; + } else { + tol = 1.e-11; + } + return tol; +} + +int main(int argc, char **argv) { + Ceed ceed; + + CeedInit(argv[1], &ceed); + + for (CeedInt dim = 1; dim <= 3; dim++) { + CeedVector x, x_nodes, x_points, u, u_points, v, ones; + CeedBasis basis_x, basis_u; + const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim); + CeedScalar sum_1 = 0, sum_2 = 0; + + CeedVectorCreate(ceed, x_dim * dim, &x); + CeedVectorCreate(ceed, p_dim * dim, &x_nodes); + CeedVectorCreate(ceed, num_points * dim, &x_points); + CeedVectorCreate(ceed, p_dim, &u); + CeedVectorCreate(ceed, num_points * dim, &u_points); + CeedVectorCreate(ceed, p_dim, &v); + CeedVectorCreate(ceed, num_points * dim, &ones); + + CeedVectorSetValue(ones, 1); + CeedVectorSetValue(v, 0); + + // Get nodal coordinates + CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x); + { + CeedScalar x_array[x_dim * dim]; + + for (CeedInt d = 0; d < dim; d++) { + for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, d + 1)) / CeedIntPow(2, d) ? 1 : -1; + } + CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes); + + // Set values of u at nodes + { + const CeedScalar *x_array; + CeedScalar u_array[p_dim]; + + CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array); + for (CeedInt i = 0; i < p_dim; i++) { + CeedScalar coord[dim]; + + for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i]; + u_array[i] = Eval(dim, coord); + } + CeedVectorRestoreArrayRead(x_nodes, &x_array); + CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); + } + + // Interpolate to arbitrary points + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u); + { + CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65}; + + CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); + } + + // Calculate G u at arbitrary points, G' * 1 at dofs + 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); + // Double it + CeedBasisApplyAddAtPoints(basis_u, 1, &num_points, CEED_TRANSPOSE, CEED_EVAL_GRAD, x_points, ones, v); + { + const CeedScalar *u_array, *v_array, *u_points_array; + + CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array); + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + CeedVectorGetArrayRead(u_points, CEED_MEM_HOST, &u_points_array); + for (CeedInt i = 0; i < p_dim; i++) sum_1 += v_array[i] * u_array[i]; + for (CeedInt i = 0; i < num_points * dim; i++) sum_2 += u_points_array[i]; + CeedVectorRestoreArrayRead(u, &u_array); + CeedVectorRestoreArrayRead(v, &v_array); + CeedVectorRestoreArrayRead(u_points, &u_points_array); + } + { + CeedScalarType scalar_type; + + CeedGetScalarType(&scalar_type); + + CeedScalar tol = GetTolerance(scalar_type, dim); + + if (fabs(sum_1 - 2.0 * sum_2) > tol) printf("[%" CeedInt_FMT "] %f != %f\n", dim, sum_1, 2.0 * sum_2); + } + + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_nodes); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&u); + CeedVectorDestroy(&u_points); + CeedVectorDestroy(&ones); + CeedVectorDestroy(&v); + CeedBasisDestroy(&basis_x); + CeedBasisDestroy(&basis_u); + } + CeedDestroy(&ceed); + return 0; +}