Skip to content

Commit

Permalink
basis - add CeedBasisApplyAtPoints (wip)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jeremylt
jeremylt committed Jun 30, 2023
1 parent 62e4d60 commit 1ad0850
Show file tree
Hide file tree
Showing 6 changed files with 315 additions and 14 deletions.
5 changes: 4 additions & 1 deletion include/ceed-impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 */
Expand All @@ -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 {
Expand Down
2 changes: 2 additions & 0 deletions include/ceed/ceed.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
238 changes: 225 additions & 13 deletions interface/ceed-basis.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
Expand Down Expand Up @@ -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");
Expand All @@ -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) {
Expand Down Expand Up @@ -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");

Expand Down Expand Up @@ -1386,6 +1383,219 @@ 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));

/*
// QR Factorization, interp_1d = Q R
CeedCall(CeedBasisGetCeed(basis, &ceed));
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 (CeedInt i = 0; i < Q_1d; i++) { // Row i
collo_grad_1d[Q_1d * i] = grad_1d[P_1d * i] / interp_1d[0];
for (CeedInt j = 1; j < P_1d; j++) { // Column j
collo_grad_1d[j + Q_1d * i] = grad_1d[j + P_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;
}
// Apply Q^T, collo_grad_1d = collo_grad_1d Q^T
CeedCall(CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, P_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, CeedIntPow(Q_1d, basis->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));
}

// 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];

// ---- Clear output array
for (CeedInt i = 0; i < num_points * num_comp; i++) v_array[i] = 0.0;
// ---- Values at point
for (CeedInt p = 0; p < num_points; p++) {
for (CeedInt d = 0; d < dim; d++) {
CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
// ------ 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
{
const CeedScalar *in = d == 0 ? chebyshev_coeffs : tmp[d % 2];
CeedScalar *out = d == dim - 1 ? &v_array[p * num_comp] : tmp[(d + 1) % 2];

for (CeedInt q = 0; q < pre * post; q++) out[q] = (CeedScalar)0.0;
for (CeedInt a = 0; a < pre; a++) {
for (CeedInt b = 0; b < Q_1d; b++) {
for (CeedInt c = 0; c < post; c++) out[a * post + c] += chebyshev_x[b] * in[(a * Q_1d + b) * post + c];
}
}
}
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
Expand Down Expand Up @@ -1695,6 +1905,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;
Expand Down
1 change: 1 addition & 0 deletions interface/ceed.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down
1 change: 1 addition & 0 deletions tests/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading

0 comments on commit 1ad0850

Please sign in to comment.