diff --git a/include/ceed-impl.h b/include/ceed-impl.h index 8aaf2ef7b6..9d4dfdfec4 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -175,6 +175,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 (*Destroy)(CeedBasis); int ref_count; bool is_tensor_basis; /* flag for tensor basis */ @@ -197,7 +198,9 @@ struct CeedBasis_private { CeedScalar *div; /* row-major matrix of shape [Q, P] expressing the divergence of basis functions at quadrature points for H(div) discretizations */ CeedScalar *curl; /* row-major matrix of shape [curl_dim * Q, P], curl_dim = 1 if dim < 3 else dim, expressing the curl of basis functions at quadrature points for H(curl) discretizations */ - void *data; /* place for the backend to store any data */ + CeedVector vec_chebyshev; + CeedBasis basis_chebyshev; /* basis interpolating from nodes to Chebyshev polynomial coefficients */ + void *data; /* place for the backend to store any data */ }; struct CeedTensorContract_private { diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index 60c1b8715d..96878e5096 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -277,6 +277,8 @@ CEED_EXTERN int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_ 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 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 c11a70b60d..87684e25dc 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -258,7 +258,6 @@ static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis bas @ref Backend **/ int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { - int i, j, k; Ceed ceed; CeedInt P_1d = (basis)->P_1d, Q_1d = (basis)->Q_1d; CeedScalar *interp_1d, *grad_1d, *tau; @@ -274,15 +273,15 @@ int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { CeedCall(CeedQRFactorization(ceed, interp_1d, tau, Q_1d, P_1d)); // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure. - // Apply Rinv, collo_grad_1d = grad_1d Rinv - for (i = 0; i < Q_1d; i++) { // Row i + // Apply R_inv, collo_grad_1d = grad_1d R_inv + for (CeedInt i = 0; i < Q_1d; i++) { // Row i collo_grad_1d[Q_1d * i] = grad_1d[P_1d * i] / interp_1d[0]; - for (j = 1; j < P_1d; j++) { // Column j + for (CeedInt j = 1; j < P_1d; j++) { // Column j collo_grad_1d[j + Q_1d * i] = grad_1d[j + P_1d * i]; - for (k = 0; k < j; k++) collo_grad_1d[j + Q_1d * i] -= interp_1d[j + P_1d * k] * collo_grad_1d[k + Q_1d * i]; + for (CeedInt k = 0; k < j; k++) collo_grad_1d[j + Q_1d * i] -= interp_1d[j + P_1d * k] * collo_grad_1d[k + Q_1d * i]; collo_grad_1d[j + Q_1d * i] /= interp_1d[j + P_1d * j]; } - for (j = P_1d; j < Q_1d; j++) collo_grad_1d[j + Q_1d * i] = 0; + for (CeedInt j = P_1d; j < Q_1d; j++) collo_grad_1d[j + Q_1d * i] = 0; } // Apply Q^T, collo_grad_1d = collo_grad_1d Q^T @@ -950,7 +949,7 @@ int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_ **/ int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, CeedBasis *basis) { // Allocate - int ierr = CEED_ERROR_SUCCESS, i, j, k; + int ierr = CEED_ERROR_SUCCESS; CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d; CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "Basis dimension must be a positive value"); @@ -977,15 +976,15 @@ int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, Ce // Build B, D matrix // Fornberg, 1998 - for (i = 0; i < Q; i++) { + for (CeedInt i = 0; i < Q; i++) { c1 = 1.0; c3 = nodes[0] - q_ref_1d[i]; interp_1d[i * P + 0] = 1.0; - for (j = 1; j < P; j++) { + for (CeedInt j = 1; j < P; j++) { c2 = 1.0; c4 = c3; c3 = nodes[j] - q_ref_1d[i]; - for (k = 0; k < j; k++) { + for (CeedInt k = 0; k < j; k++) { dx = nodes[j] - nodes[k]; c2 *= dx; if (k == j - 1) { @@ -1353,9 +1352,7 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); CeedCall(CeedVectorGetLength(v, &v_length)); - if (u) { - CeedCall(CeedVectorGetLength(u, &u_length)); - } + if (u) CeedCall(CeedVectorGetLength(u, &u_length)); CeedCheck(basis->Apply, basis->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support BasisApply"); @@ -1386,6 +1383,198 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, return CEED_ERROR_SUCCESS; } +/** + @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] 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 + @param[in] x_ref CeedVector holding reference coordinates of each point + @param[in] u Input CeedVector, of length `num_nodes * num_comp` for `CEED_NOTRANSPOSE` + @param[out] v Output CeedVector, of length `num_points * num_q_comp` for `CEED_NOTRANSPOSE` with `CEED_EVAL_INTERP` + + @return An error code: 0 - success, otherwise - failure + + @ref User +**/ +int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, + CeedVector v) { + CeedSize x_length = 0, u_length = 0, v_length; + CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1; + + CeedCall(CeedBasisGetDimension(basis, &dim)); + CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); + CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); + CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &num_q_comp)); + CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); + CeedCall(CeedVectorGetLength(x_ref, &x_length)); + CeedCall(CeedVectorGetLength(v, &v_length)); + CeedCall(CeedVectorGetLength(u, &u_length)); + + // Check compatibility of topological and geometrical dimensions + CeedCheck((t_mode == CEED_TRANSPOSE && v_length % num_nodes == 0) || (t_mode == CEED_NOTRANSPOSE && u_length % num_nodes == 0), basis->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, basis->ceed, CEED_ERROR_DIMENSION, + "Length of reference coordinate vector incompatible with basis dimension and number of points"); + + // Check vector lengths to prevent out of bounds issues + bool good_dims = false; + switch (eval_mode) { + case CEED_EVAL_INTERP: + 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))); + break; + case CEED_EVAL_GRAD: + case CEED_EVAL_NONE: + case CEED_EVAL_WEIGHT: + case CEED_EVAL_DIV: + case CEED_EVAL_CURL: + // LCOV_EXCL_START + return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points not supported for %s", CeedEvalModes[eval_mode]); + // LCOV_EXCL_STOP + } + CeedCheck(good_dims, basis->ceed, CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); + + // Backend method + if (basis->ApplyAtPoints) { + CeedCall(basis->ApplyAtPoints(basis, num_points, t_mode, eval_mode, x_ref, u, v)); + return CEED_ERROR_SUCCESS; + } + + // Default implementation + CeedCheck(basis->is_tensor_basis, basis->ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points only supported for tensor product bases"); + if (!basis->basis_chebyshev) { + // Build matrix mapping from quadrature point values to Chebyshev coefficients + CeedScalar *tau, *C, *I, *chebyshev_coeffs_1d; + const CeedScalar *q_ref_1d; + + // Build coefficient matrix + // -- Note: Clang-tidy needs this check because it does not understand the is_tensor_basis check above + CeedCheck(P_1d > 0 && Q_1d > 0, basis->ceed, CEED_ERROR_INCOMPATIBLE, "Basis dimensions are malformed"); + CeedCall(CeedCalloc(Q_1d * Q_1d, &C)); + CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); + for (CeedInt i = 0; i < Q_1d; i++) { + const CeedScalar x = q_ref_1d[i]; + + C[i * Q_1d + 0] = 1.0; + C[i * Q_1d + 1] = 2 * x; + for (CeedInt j = 2; j < Q_1d; j++) C[i * Q_1d + j] = 2 * x * C[i * Q_1d + j - 1] - C[i * Q_1d + j - 2]; + } + + // Inverse of coefficient matrix + CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_coeffs_1d)); + CeedCall(CeedCalloc(Q_1d * Q_1d, &I)); + CeedCall(CeedCalloc(Q_1d, &tau)); + // -- QR Factorization, C = Q R + CeedCall(CeedQRFactorization(basis->ceed, C, tau, Q_1d, Q_1d)); + // -- chebyshev_coeffs_1d = R_inv Q^T + for (CeedInt i = 0; i < Q_1d; i++) I[i * Q_1d + i] = 1.0; + // ---- Apply R_inv, chebyshev_coeffs_1d = I R_inv + for (CeedInt i = 0; i < Q_1d; i++) { // Row i + chebyshev_coeffs_1d[Q_1d * i] = I[Q_1d * i] / C[0]; + for (CeedInt j = 1; j < Q_1d; j++) { // Column j + chebyshev_coeffs_1d[j + Q_1d * i] = I[j + Q_1d * i]; + for (CeedInt k = 0; k < j; k++) chebyshev_coeffs_1d[j + Q_1d * i] -= C[j + Q_1d * k] * chebyshev_coeffs_1d[k + Q_1d * i]; + chebyshev_coeffs_1d[j + Q_1d * i] /= C[j + Q_1d * j]; + } + } + // ---- Apply Q^T, chebyshev_coeffs_1d = R_inv Q^T + CeedCall(CeedHouseholderApplyQ(chebyshev_coeffs_1d, C, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, Q_1d, 1, Q_1d)); + + // Build basis mapping from nodes to Chebyshev coefficients + CeedScalar *chebyshev_interp_1d, *chebyshev_grad_1d, *chebyshev_q_weight_1d; + const CeedScalar *interp_1d; + + CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_interp_1d)); + CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_grad_1d)); + CeedCall(CeedCalloc(Q_1d, &chebyshev_q_weight_1d)); + CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); + CeedCall(CeedMatrixMatrixMultiply(basis->ceed, chebyshev_coeffs_1d, interp_1d, chebyshev_interp_1d, Q_1d, P_1d, Q_1d)); + + CeedCall(CeedVectorCreate(basis->ceed, num_comp * CeedIntPow(Q_1d, dim), &basis->vec_chebyshev)); + CeedCall(CeedBasisCreateTensorH1(basis->ceed, dim, num_comp, Q_1d, Q_1d, chebyshev_interp_1d, chebyshev_grad_1d, q_ref_1d, chebyshev_q_weight_1d, + &basis->basis_chebyshev)); + + // Cleanup + CeedCall(CeedFree(&C)); + CeedCall(CeedFree(&chebyshev_coeffs_1d)); + CeedCall(CeedFree(&I)); + CeedCall(CeedFree(&tau)); + CeedCall(CeedFree(&chebyshev_interp_1d)); + CeedCall(CeedFree(&chebyshev_grad_1d)); + CeedCall(CeedFree(&chebyshev_q_weight_1d)); + } + + // Create TensorContract object if needed, such as a basis from the GPU backends + if (!basis->contract) { + Ceed ceed_ref; + CeedBasis basis_ref; + + CeedCall(CeedInit("/cpu/self", &ceed_ref)); + // Only need matching tensor contraction dimensions, any type of basis will work + CeedCall(CeedBasisCreateTensorH1Lagrange(ceed_ref, dim, num_comp, Q_1d, Q_1d, CEED_GAUSS, &basis_ref)); + CeedCall(CeedTensorContractReference(basis_ref->contract)); + basis->contract = basis_ref->contract; + CeedCall(CeedBasisDestroy(&basis_ref)); + CeedCall(CeedDestroy(&ceed_ref)); + } + + // Basis evaluation + switch (t_mode) { + case CEED_NOTRANSPOSE: { + // Nodes to arbitrary points + CeedScalar *v_array; + const CeedScalar *chebyshev_coeffs, *x_array_read; + + // -- Interpolate to Chebyshev coefficients + CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, basis->vec_chebyshev)); + + // -- Evaluate Chebyshev polynomials at arbitrary points + CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); + CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read)); + CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array)); + { + CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; + + // ---- Values at point + for (CeedInt p = 0; p < num_points; p++) { + CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; + + for (CeedInt d = dim - 1; d >= 0; d--) { + // ------ Compute Chebyshev polynomial values + { + const CeedScalar x = x_array_read[p * dim + d]; + + chebyshev_x[0] = 1.0; + chebyshev_x[1] = 2 * x; + for (CeedInt j = 2; j < Q_1d; j++) chebyshev_x[j] = 2 * x * chebyshev_x[j - 1] - chebyshev_x[j - 2]; + } + // ------ Tensor contract + CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, + d == (dim - 1) ? chebyshev_coeffs : tmp[d % 2], d == 0 ? &v_array[p * num_comp] : tmp[(d + 1) % 2])); + pre /= Q_1d; + post *= 1; + } + } + } + CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs)); + CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read)); + CeedCall(CeedVectorRestoreArray(v, &v_array)); + break; + } + case CEED_TRANSPOSE: + return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, "CEED_TRANSPOSE unsupported for arbitrary basis point evaluation"); + } + + return CEED_ERROR_SUCCESS; +} + /** @brief Get Ceed associated with a CeedBasis @@ -1695,6 +1884,8 @@ int CeedBasisDestroy(CeedBasis *basis) { CeedCall(CeedFree(&(*basis)->grad_1d)); CeedCall(CeedFree(&(*basis)->div)); CeedCall(CeedFree(&(*basis)->curl)); + CeedCall(CeedVectorDestroy(&(*basis)->vec_chebyshev)); + CeedCall(CeedBasisDestroy(&(*basis)->basis_chebyshev)); CeedCall(CeedDestroy(&(*basis)->ceed)); CeedCall(CeedFree(basis)); return CEED_ERROR_SUCCESS; diff --git a/interface/ceed.c b/interface/ceed.c index 537e4e970f..dbe837fab1 100644 --- a/interface/ceed.c +++ b/interface/ceed.c @@ -841,6 +841,7 @@ int CeedInit(const char *resource, Ceed *ceed) { CEED_FTABLE_ENTRY(CeedElemRestriction, GetOffsets), CEED_FTABLE_ENTRY(CeedElemRestriction, Destroy), CEED_FTABLE_ENTRY(CeedBasis, Apply), + CEED_FTABLE_ENTRY(CeedBasis, ApplyAtPoints), CEED_FTABLE_ENTRY(CeedBasis, Destroy), CEED_FTABLE_ENTRY(CeedTensorContract, Apply), CEED_FTABLE_ENTRY(CeedTensorContract, Destroy), diff --git a/tests/README.md b/tests/README.md index a7a5600a7b..0d4fcec581 100644 --- a/tests/README.md +++ b/tests/README.md @@ -15,6 +15,7 @@ The tests are organized by API object, and some tests are further organized, as 3.2. CeedBasis simplex basis tests 3.3. CeedBasis non-tensor H(div) basis tests 3.4. CeedBasis non-tensor H(curl) basis tests + 3.5. CeedBasis evaluation at arbitrary points tests 4. CeedQFunction Tests 4.0. CeedQFunction user code tests 4.1. CeedQFunction gallery code tests diff --git a/tests/t311-basis.c b/tests/t311-basis.c index 489513747b..93929648a4 100644 --- a/tests/t311-basis.c +++ b/tests/t311-basis.c @@ -32,7 +32,7 @@ int main(int argc, char **argv) { { CeedScalar x_array[2]; - for (int i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1); + for (CeedInt i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1); CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array); } diff --git a/tests/t313-basis.c b/tests/t313-basis.c index f5f2fa6d87..f834c22df9 100644 --- a/tests/t313-basis.c +++ b/tests/t313-basis.c @@ -18,6 +18,7 @@ int main(int argc, char **argv) { Ceed ceed; CeedInit(argv[1], &ceed); + for (CeedInt dim = 1; dim <= 3; dim++) { CeedVector x, x_q, u, u_q; CeedBasis basis_x_lobatto, basis_u_lobatto, basis_x_gauss, basis_u_gauss; @@ -77,8 +78,8 @@ int main(int argc, char **argv) { CeedScalar fx = Eval(dim, coord); if (fabs(u_array[i] - fx) > 1E-4) { // LCOV_EXCL_START - printf("[%" CeedInt_FMT "] %f != %f=f(%f", dim, u_array[i], fx, coord[0]); - for (CeedInt d = 1; d < dim; d++) printf(",%f", coord[d]); + printf("[%" CeedInt_FMT "] %f != %f = f(%f", dim, u_array[i], fx, coord[0]); + for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); puts(")"); // LCOV_EXCL_STOP } diff --git a/tests/t350-basis.c b/tests/t350-basis.c new file mode 100644 index 0000000000..1b02c0b99d --- /dev/null +++ b/tests/t350-basis.c @@ -0,0 +1,83 @@ +/// @file +/// Test polynomial interpolation to arbirtary points in 1D +/// \test Test polynomial interpolation to 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, u, v; + 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, p, &u); + CeedVectorCreate(ceed, num_points, &v); + + // 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, .16, 0.99}; + + 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); + + { + const CeedScalar *x_array, *v_array; + + CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < num_points; i++) { + CeedScalar fx = Eval(x_array[i], ALEN(c), c); + if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) printf("%f != %f = f(%f)\n", v_array[i], fx, x_array[i]); + } + CeedVectorRestoreArrayRead(x_points, &x_array); + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_nodes); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedBasisDestroy(&basis_x); + CeedBasisDestroy(&basis_u); + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t351-basis.c b/tests/t351-basis.c new file mode 100644 index 0000000000..40092916ed --- /dev/null +++ b/tests/t351-basis.c @@ -0,0 +1,103 @@ +/// @file +/// Test polynomial interpolation to arbirtary points in multiple dimensions +/// \test Test polynomial interpolation to arbitrary points in multiple dimensions +#include +#include +#include + +static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) { + CeedScalar result = 1, center = 0.1; + for (CeedInt d = 0; d < dim; d++) { + result *= tanh(x[d] - center); + center += 0.1; + } + return result; +} + +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, v; + CeedBasis basis_x, basis_u; + const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim); + + 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, &v); + + // 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, dim - d)) / CeedIntPow(2, dim - d - 1) ? 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, .16, 0.99, -0.65, .16, 0.99, -0.33, .16, 0.99, -0.33, -0.65}; + + 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); + + { + const CeedScalar *x_array, *v_array; + + CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < num_points; i++) { + CeedScalar coord[dim]; + + for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d + i * dim]; + const CeedScalar fx = Eval(dim, coord); + if (fabs(v_array[i] - fx) > 1E-4) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT "] %f != %f = f(%f", dim, v_array[i], fx, coord[0]); + for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); + puts(")"); + // LCOV_EXCL_STOP + } + } + CeedVectorRestoreArrayRead(x_points, &x_array); + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_nodes); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedBasisDestroy(&basis_x); + CeedBasisDestroy(&basis_u); + } + + CeedDestroy(&ceed); + return 0; +} diff --git a/tests/t352-basis.c b/tests/t352-basis.c new file mode 100644 index 0000000000..378341cba2 --- /dev/null +++ b/tests/t352-basis.c @@ -0,0 +1,105 @@ +/// @file +/// Test polynomial interpolation to arbirtary points with multiple components in multiple dimensions +/// \test Test polynomial interpolation to arbitrary points with multiple components in multiple dimensions +#include +#include +#include + +static CeedScalar Eval(CeedInt dim, CeedScalar scale, const CeedScalar x[]) { + CeedScalar result = 1, center = 0.1; + for (CeedInt d = 0; d < dim; d++) { + result *= tanh(x[d] - center); + center += 0.1; + } + return scale * result; +} + +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, v; + CeedBasis basis_x, basis_u; + const CeedInt p = 9, q = 9, num_comp = 3, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim); + + CeedVectorCreate(ceed, x_dim * dim, &x); + CeedVectorCreate(ceed, p_dim * dim, &x_nodes); + CeedVectorCreate(ceed, num_points * dim, &x_points); + CeedVectorCreate(ceed, num_comp * p_dim, &u); + CeedVectorCreate(ceed, num_comp * num_points, &v); + + // 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, dim - d)) / CeedIntPow(2, dim - d - 1) ? 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[num_comp * 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]; + for (CeedInt c = 0; c < num_comp; c++) u_array[i + c * p_dim] = Eval(dim, c, coord); + } + CeedVectorRestoreArrayRead(x_nodes, &x_array); + CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array); + } + + // Interpolate to arbitrary points + CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, p, q, CEED_GAUSS, &basis_u); + { + CeedScalar x_array[12] = {-0.33, -0.65, .16, 0.99, -0.65, .16, 0.99, -0.33, .16, 0.99, -0.33, -0.65}; + + 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); + + { + const CeedScalar *x_array, *v_array; + + CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array); + CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array); + for (CeedInt i = 0; i < num_points; i++) { + CeedScalar coord[dim]; + + for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d + i * dim]; + for (CeedInt c = 0; c < num_comp; c++) { + CeedScalar fx = Eval(dim, c, coord); + if (fabs(v_array[c + i * num_comp] - fx) > 1E-4) { + // LCOV_EXCL_START + printf("[%" CeedInt_FMT ", %" CeedInt_FMT "] %f != %f = f(%f", dim, c, v_array[c + i * num_comp], fx, coord[0]); + for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]); + puts(")"); + // LCOV_EXCL_STOP + } + } + } + CeedVectorRestoreArrayRead(x_points, &x_array); + CeedVectorRestoreArrayRead(v, &v_array); + } + + CeedVectorDestroy(&x); + CeedVectorDestroy(&x_nodes); + CeedVectorDestroy(&x_points); + CeedVectorDestroy(&u); + CeedVectorDestroy(&v); + CeedBasisDestroy(&basis_x); + CeedBasisDestroy(&basis_u); + } + + CeedDestroy(&ceed); + return 0; +}