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 May 10, 2023
1 parent be301fa commit 8d0f4d9
Show file tree
Hide file tree
Showing 4 changed files with 193 additions and 8 deletions.
4 changes: 3 additions & 1 deletion include/ceed-impl.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,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 @@ -196,7 +197,8 @@ 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 */
CeedScalar *legendre_coeffs_1d; /* row-major matrix mapping from quadrature point values to Legendre 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 @@ -276,6 +276,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
194 changes: 187 additions & 7 deletions interface/ceed-basis.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -951,7 +951,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 @@ -978,15 +978,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 @@ -1360,9 +1360,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 @@ -1393,6 +1391,187 @@ 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, Q_1d = 0;

CeedCall(CeedBasisGetDimension(basis, &dim));
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
if (!basis->legendre_coeffs_1d) {
// Build matrix mapping from quadrature point values to Legendre coefficients
if (basis->is_tensor_basis) {
CeedScalar *tau, *C, *Q_inv;
const CeedScalar *q_ref;

// Build coefficient matrix
CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
CeedCall(CeedCalloc(Q_1d * Q_1d, &C));
CeedCall(CeedBasisGetQRef(basis, &q_ref));
for (CeedInt i = 0; i < Q_1d; i++) {
const CeedScalar x = q_ref[i];

C[i * Q_1d + 0] = 1.0;
C[i * Q_1d + 1] = x;
for (CeedInt j = 2; j < Q_1d; j++) {
C[i * Q_1d + j] = (2 * (j - 1) + 1) * x * C[i * Q_1d + j - 1] - (j - 1) * C[i * Q_1d + j - 2];
}
}

// Inverse of coefficient matrix
CeedCall(CeedCalloc(Q_1d * Q_1d, &basis->legendre_coeffs_1d));
CeedCall(CeedCalloc(Q_1d * Q_1d, &Q_inv));
CeedCall(CeedCalloc(Q_1d, &tau));
// -- QR Factorization, C = Q R
CeedCall(CeedQRFactorization(basis->ceed, C, tau, Q_1d, Q_1d));
// -- Cinv = Q^T Rinv
for (CeedInt i = 0; i < Q_1d; i++) Q_inv[i * Q_1d + i] = 1.0;
CeedCall(CeedHouseholderApplyQ(Q_inv, C, tau, CEED_TRANSPOSE, Q_1d, Q_1d, Q_1d, Q_1d, 1));
for (CeedInt j = 0; j < Q_1d; j++) { // Column j
basis->legendre_coeffs_1d[j + Q_1d * (Q_1d - 1)] = Q_inv[j + Q_1d * (Q_1d - 1)] / C[Q_1d * Q_1d - 1];
for (CeedInt i = Q_1d - 2; i >= 0; i--) { // Row i
basis->legendre_coeffs_1d[j + Q_1d * i] = Q_inv[j + Q_1d * i];
for (CeedInt k = i + 1; k < Q_1d; k++) {
basis->legendre_coeffs_1d[j + Q_1d * i] -= C[k + Q_1d * i] * basis->legendre_coeffs_1d[j + Q_1d * k];
}
basis->legendre_coeffs_1d[j + Q_1d * i] /= C[i + Q_1d * i];
}
}
// Cleanup
CeedCall(CeedFree(&Q_inv));
CeedCall(CeedFree(&tau));
CeedCall(CeedFree(&C));
} else {
// ToDo: Add non-tensor implementation
return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAtPoints not supported for non-tensor basis");
}
}

// Basis evaluation
CeedVector q = NULL;
const CeedScalar *x_ref_array = NULL;

CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_ref_array));
switch (t_mode) {
case CEED_NOTRANSPOSE: {
// Nodes to arbitrary points
CeedInt Q, Q_comp;
CeedScalar *v_array = NULL, *coeffs;

// -- Interpolate to quadrature points
CeedCall(CeedBasisGetNumQuadraturePoints(basis, &Q));
CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &Q_comp));
CeedCall(CeedVectorCreate(basis->ceed, Q * Q_comp, &q));
CeedCall(CeedBasisApply(basis, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, q));

// -- Map to Legendre coefficients
// A tensor contraction would be really handy here!

// -- Evaluate Legendre polynomial at arbitrary points
CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array));
for (CeedInt p = 0; p < num_points; p++) {
CeedInt P[dim * Q_1d];

// ---- Build Legendre polynomial values at points
for (CeedInt d = 0; d < dim; d++) {
P[d * Q_1d + 0] = 1.0;
P[d * Q_1d + 1] = x_ref_array[p * dim + d];
for (CeedInt j = 2; j < Q_1d; j++) {
P[d * Q_1d + j] = (2 * (j - 1) + 1) * x_ref_array[p * dim + d] * P[d * Q_1d + j - 1] - (j - 1) * P[d * Q_1d + j - 2];
}
}
// ---- Evaluate at points
for (CeedInt c = 0; c < num_q_comp; c++) v_array[p * num_q_comp + c] = 0.0;
for (CeedInt q = 0; q < Q; q++) {
CeedScalar f_x = 1.0;

for (CeedInt d = 0; d < dim; d++) {
CeedInt qq = (q / CeedIntPow(Q_1d, d)) % Q_1d;
f_x *= P[d * Q_1d + qq];
}
for (CeedInt c = 0; c < num_q_comp; c++) v_array[p * num_q_comp + c] += f_x * coeffs[q * num_q_comp + c];
}
}
CeedCall(CeedVectorRestoreArray(v, &v_array));
break;
}
case CEED_TRANSPOSE: {
// Arbitrary points to nodes
CeedScalar *u_array = NULL;
CeedInt Q, Q_comp;

// -- Arbitrary points to quadrature points
CeedCall(CeedBasisGetNumQuadraturePoints(basis, &Q));
CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &Q_comp));
CeedCall(CeedVectorCreate(basis->ceed, Q * Q_comp, &q));

// -- Transpose interpolation from quadrature points
CeedCall(CeedBasisApply(basis, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, q, u));
}
}
// -- Cleanup
CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_ref_array));
CeedCall(CeedVectorDestroy(&q));

return CEED_ERROR_SUCCESS;
}

/**
@brief Get Ceed associated with a CeedBasis
Expand Down Expand Up @@ -1702,6 +1881,7 @@ int CeedBasisDestroy(CeedBasis *basis) {
CeedCall(CeedFree(&(*basis)->grad_1d));
CeedCall(CeedFree(&(*basis)->div));
CeedCall(CeedFree(&(*basis)->curl));
CeedCall(CeedFree(&(*basis)->legendre_coeffs_1d));
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 @@ -819,6 +819,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

0 comments on commit 8d0f4d9

Please sign in to comment.