diff --git a/backends/ref/ceed-ref-basis.c b/backends/ref/ceed-ref-basis.c index c4314dbcae..6c10f93af8 100644 --- a/backends/ref/ceed-ref-basis.c +++ b/backends/ref/ceed-ref-basis.c @@ -24,8 +24,6 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp)); CeedCallBackend(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); - CeedFESpace fe_space; - CeedCall(CeedBasisGetFESpace(basis, &fe_space)); CeedTensorContract contract; CeedCallBackend(CeedBasisGetTensorContract(basis, &contract)); const CeedInt add = (t_mode == CEED_TRANSPOSE); @@ -195,12 +193,13 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo switch (eval_mode) { // Interpolate to/from quadrature points case CEED_EVAL_INTERP: { - CeedInt qdim = (fe_space == CEED_FE_SPACE_H1) ? 1 : dim; - CeedInt P = num_nodes, Q = qdim * num_qpts; + CeedInt q_comp; + CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, &q_comp)); + CeedInt P = num_nodes, Q = q_comp * num_qpts; const CeedScalar *interp; CeedCallBackend(CeedBasisGetInterp(basis, &interp)); if (t_mode == CEED_TRANSPOSE) { - P = qdim * num_qpts; + P = q_comp * num_qpts; Q = num_nodes; } CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, interp, t_mode, add, u, v)); @@ -250,12 +249,13 @@ static int CeedBasisApply_Ref(CeedBasis basis, CeedInt num_elem, CeedTransposeMo } break; // Evaluate the curl to/from the quadrature points case CEED_EVAL_CURL: { - CeedInt cdim = (dim < 3) ? 1 : dim; - CeedInt P = num_nodes, Q = cdim * num_qpts; + CeedInt curl_comp; + CeedCallBackend(CeedBasisGetNumCurlComponents(basis, &curl_comp)); + CeedInt P = num_nodes, Q = curl_comp * num_qpts; const CeedScalar *curl; CeedCallBackend(CeedBasisGetCurl(basis, &curl)); if (t_mode == CEED_TRANSPOSE) { - P = cdim * num_qpts; + P = curl_comp * num_qpts; Q = num_nodes; } CeedCallBackend(CeedTensorContractApply(contract, num_comp, P, num_elem, Q, curl, t_mode, add, u, v)); diff --git a/include/ceed-impl.h b/include/ceed-impl.h index 6476930dfa..279297bd72 100644 --- a/include/ceed-impl.h +++ b/include/ceed-impl.h @@ -191,8 +191,8 @@ struct CeedBasis_private { CeedScalar *grad; /* row-major matrix of shape [dim * Q, P] matrix expressing derivatives of nodal basis functions at quadrature points */ CeedScalar *grad_1d; /* row-major matrix of shape [Q1d, P1d] matrix expressing derivatives of nodal basis functions at quadrature points */ 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 [cdim * Q, P], cdim = 1 if dim < 3 else dim, expressing the curl of basis functions at quadrature - points for H(curl) 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 */ }; diff --git a/include/ceed/backend.h b/include/ceed/backend.h index 17737059c8..f25ff6e8cf 100644 --- a/include/ceed/backend.h +++ b/include/ceed/backend.h @@ -280,4 +280,9 @@ CEED_EXTERN int CeedOperatorSetData(CeedOperator op, void *data); CEED_EXTERN int CeedOperatorReference(CeedOperator op); CEED_EXTERN int CeedOperatorSetSetupDone(CeedOperator op); +CEED_EXTERN int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n, + CeedInt kk); +CEED_EXTERN int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n, + CeedInt k, CeedInt row, CeedInt col); + #endif diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index 2f8be19bbc..6bffab4509 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -372,6 +372,8 @@ CEED_EXTERN int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed); CEED_EXTERN int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim); CEED_EXTERN int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo); CEED_EXTERN int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp); +CEED_EXTERN int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedInt *q_comp); +CEED_EXTERN int CeedBasisGetNumCurlComponents(CeedBasis basis, CeedInt *curl_comp); CEED_EXTERN int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P); CEED_EXTERN int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d); CEED_EXTERN int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q); @@ -391,10 +393,6 @@ CEED_EXTERN int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScala CEED_EXTERN int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n); CEED_EXTERN int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n); CEED_EXTERN int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *x, CeedScalar *lambda, CeedInt n); -CEED_EXTERN int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n, - CeedInt k, CeedInt row, CeedInt col); -CEED_EXTERN int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n, - CeedInt kk); /** Handle for the user provided CeedQFunction callback function diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index 5602bb7a2a..2cae78b22d 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -396,19 +396,19 @@ int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEva break; } } else { - CeedInt dim, num_comp, num_nodes, num_qpts; + CeedInt dim, num_comp, q_comp, curl_comp, num_nodes, num_qpts; CeedCall(CeedBasisGetDimension(basis, &dim)); CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, &q_comp)); + CeedCall(CeedBasisGetNumCurlComponents(basis, &curl_comp)); CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); - CeedFESpace fe_space; - CeedCall(CeedBasisGetFESpace(basis, &fe_space)); switch (eval_mode) { case CEED_EVAL_NONE: *flops = 0; break; case CEED_EVAL_INTERP: - *flops = num_nodes * num_qpts * num_comp * (fe_space == CEED_FE_SPACE_H1 ? 1 : dim); + *flops = num_nodes * num_qpts * num_comp * q_comp; break; case CEED_EVAL_GRAD: *flops = num_nodes * num_qpts * num_comp * dim; @@ -417,7 +417,7 @@ int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEva *flops = num_nodes * num_qpts * num_comp; break; case CEED_EVAL_CURL: - *flops = num_nodes * num_qpts * num_comp * (dim < 3 ? 1 : dim); + *flops = num_nodes * num_qpts * num_comp * curl_comp; break; case CEED_EVAL_WEIGHT: *flops = 0; @@ -429,32 +429,32 @@ int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEva } /** - @brief Get dimension for given CeedElemTopology + @brief Get CeedFESpace for a CeedBasis - @param[in] topo CeedElemTopology - @param[out] dim Variable to store dimension of topology + @param[in] basis CeedBasis + @param[out] fe_space Variable to store CeedFESpace @return An error code: 0 - success, otherwise - failure @ref Backend **/ -int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) { - *dim = (CeedInt)topo >> 16; +int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space) { + *fe_space = basis->fe_space; return CEED_ERROR_SUCCESS; } /** - @brief Get CeedFESpace for a CeedBasis + @brief Get dimension for given CeedElemTopology - @param[in] basis CeedBasis - @param[out] fe_space Variable to store CeedFESpace + @param[in] topo CeedElemTopology + @param[out] dim Variable to store dimension of topology @return An error code: 0 - success, otherwise - failure @ref Backend **/ -int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space) { - *fe_space = basis->fe_space; +int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) { + *dim = (CeedInt)topo >> 16; return CEED_ERROR_SUCCESS; } @@ -489,6 +489,66 @@ int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract) { return CEED_ERROR_SUCCESS; } +/** + @brief Return a reference implementation of matrix multiplication C = A B. + Note, this is a reference implementation for CPU CeedScalar pointers that is not intended for high performance. + + @param[in] ceed Ceed context for error handling + @param[in] mat_A Row-major matrix A + @param[in] mat_B Row-major matrix B + @param[out] mat_C Row-major output matrix C + @param[in] m Number of rows of C + @param[in] n Number of columns of C + @param[in] kk Number of columns of A/rows of B + + @return An error code: 0 - success, otherwise - failure + + @ref Utility +**/ +int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n, CeedInt kk) { + for (CeedInt i = 0; i < m; i++) { + for (CeedInt j = 0; j < n; j++) { + CeedScalar sum = 0; + for (CeedInt k = 0; k < kk; k++) sum += mat_A[k + i * kk] * mat_B[j + k * n]; + mat_C[j + i * n] = sum; + } + } + return CEED_ERROR_SUCCESS; +} + +/** + @brief Apply Householder Q matrix + + Compute A = Q A, where Q is mxm and A is mxn. + + @param[in,out] A Matrix to apply Householder Q to, in place + @param[in] Q Householder Q matrix + @param[in] tau Householder scaling factors + @param[in] t_mode Transpose mode for application + @param[in] m Number of rows in A + @param[in] n Number of columns in A + @param[in] k Number of elementary reflectors in Q, kQ_1d = Q_1d; (*basis)->P = CeedIntPow(P_1d, dim); (*basis)->Q = CeedIntPow(Q_1d, dim); - (*basis)->fe_space = 1; // H^1 space + (*basis)->fe_space = CEED_FE_SPACE_H1; CeedCall(CeedCalloc(Q_1d, &(*basis)->q_ref_1d)); CeedCall(CeedCalloc(Q_1d, &(*basis)->q_weight_1d)); if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d * sizeof(q_ref_1d[0])); @@ -745,7 +805,7 @@ int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedIn (*basis)->num_comp = num_comp; (*basis)->P = P; (*basis)->Q = Q; - (*basis)->fe_space = 1; // H^1 space + (*basis)->fe_space = CEED_FE_SPACE_H1; CeedCall(CeedCalloc(Q * dim, &(*basis)->q_ref_1d)); CeedCall(CeedCalloc(Q, &(*basis)->q_weight_1d)); if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); @@ -778,8 +838,8 @@ int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedIn **/ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { - CeedInt Q = num_qpts, P = num_nodes, dim = 0; - CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); + CeedInt Q = num_qpts, P = num_nodes, dim = 0, q_comp = 0; + if (!ceed->BasisCreateHdiv) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); @@ -814,6 +874,9 @@ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Ceed CeedCall(CeedCalloc(1, basis)); + CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); + q_comp = dim; + (*basis)->ceed = ceed; CeedCall(CeedReference(ceed)); (*basis)->ref_count = 1; @@ -823,14 +886,14 @@ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Ceed (*basis)->num_comp = num_comp; (*basis)->P = P; (*basis)->Q = Q; - (*basis)->fe_space = 2; // H(div) space + (*basis)->fe_space = CEED_FE_SPACE_HDIV; CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); - CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); + CeedCall(CeedMalloc(q_comp * Q * P, &(*basis)->interp)); CeedCall(CeedMalloc(Q * P, &(*basis)->div)); - if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); + if (interp) memcpy((*basis)->interp, interp, q_comp * Q * P * sizeof(interp[0])); if (div) memcpy((*basis)->div, div, Q * P * sizeof(div[0])); CeedCall(ceed->BasisCreateHdiv(topo, dim, P, Q, interp, div, q_ref, q_weight, *basis)); return CEED_ERROR_SUCCESS; @@ -845,8 +908,8 @@ int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Ceed @param[in] num_nodes Total number of nodes (dofs per element) @param[in] num_qpts Total number of quadrature points @param[in] interp Row-major (dim * num_qpts * num_nodes) matrix expressing the values of basis functions at quadrature points - @param[in] curl Row-major (cdim * num_qpts * num_nodes, cdim = 1 if dim < 3 else dim) matrix expressing curl of basis functions at quadrature -points + @param[in] curl Row-major (curl_comp * num_qpts * num_nodes, curl_comp = 1 if dim < 3 else dim) matrix expressing curl of basis functions at +quadrature points @param[in] q_ref Array of length num_qpts * dim holding the locations of quadrature points on the reference element @param[in] q_weight Array of length num_qpts holding the quadrature weights on the reference element @param[out] basis Address of the variable where the newly created CeedBasis will be stored. @@ -857,9 +920,8 @@ points **/ int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { - CeedInt Q = num_qpts, P = num_nodes, dim = 0, cdim; - CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); - cdim = (dim < 3) ? 1 : dim; + CeedInt Q = num_qpts, P = num_nodes, dim = 0, q_comp = 0, curl_comp = 0; + if (!ceed->BasisCreateHdiv) { Ceed delegate; CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); @@ -894,6 +956,10 @@ int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Cee CeedCall(CeedCalloc(1, basis)); + CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); + q_comp = dim; + curl_comp = (dim < 3) ? 1 : dim; + (*basis)->ceed = ceed; CeedCall(CeedReference(ceed)); (*basis)->ref_count = 1; @@ -903,15 +969,15 @@ int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, Cee (*basis)->num_comp = num_comp; (*basis)->P = P; (*basis)->Q = Q; - (*basis)->fe_space = 3; // H(curl) space + (*basis)->fe_space = CEED_FE_SPACE_HCURL; CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); - CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); - CeedCall(CeedMalloc(cdim * Q * P, &(*basis)->curl)); - if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); - if (curl) memcpy((*basis)->curl, curl, cdim * Q * P * sizeof(curl[0])); + CeedCall(CeedMalloc(q_comp * Q * P, &(*basis)->interp)); + CeedCall(CeedMalloc(curl_comp * Q * P, &(*basis)->curl)); + if (interp) memcpy((*basis)->interp, interp, q_comp * Q * P * sizeof(interp[0])); + if (curl) memcpy((*basis)->curl, curl, curl_comp * Q * P * sizeof(curl[0])); CeedCall(ceed->BasisCreateHcurl(topo, dim, P, Q, interp, curl, q_ref, q_weight, *basis)); return CEED_ERROR_SUCCESS; } @@ -1006,8 +1072,10 @@ int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) { @ref User **/ int CeedBasisView(CeedBasis basis, FILE *stream) { - CeedFESpace fe_space = basis->fe_space; CeedElemTopology topo = basis->topo; + CeedFESpace fe_space = basis->fe_space; + CeedInt q_comp = 0; + CeedCall(CeedBasisGetNumQuadratureComponents(basis, &q_comp)); // Print FE space and element topology of the basis if (basis->tensor_basis) { @@ -1026,7 +1094,7 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { } else { // non-tensor basis CeedCall(CeedScalarView("qref", "\t% 12.8f", 1, basis->Q * basis->dim, basis->q_ref_1d, stream)); CeedCall(CeedScalarView("qweight", "\t% 12.8f", 1, basis->Q, basis->q_weight_1d, stream)); - CeedCall(CeedScalarView("interp", "\t% 12.8f", (fe_space == CEED_FE_SPACE_H1 ? 1 : basis->dim) * basis->Q, basis->P, basis->interp, stream)); + CeedCall(CeedScalarView("interp", "\t% 12.8f", q_comp * basis->Q, basis->P, basis->interp, stream)); if (basis->grad) { CeedCall(CeedScalarView("grad", "\t% 12.8f", basis->dim * basis->Q, basis->P, basis->grad, stream)); } @@ -1034,7 +1102,9 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { CeedCall(CeedScalarView("div", "\t% 12.8f", basis->Q, basis->P, basis->div, stream)); } if (basis->curl) { - CeedCall(CeedScalarView("curl", "\t% 12.8f", (basis->dim < 3 ? 1 : basis->dim) * basis->Q, basis->P, basis->curl, stream)); + CeedInt curl_comp = 0; + CeedCall(CeedBasisGetNumQuadratureComponents(basis, &curl_comp)); + CeedCall(CeedScalarView("curl", "\t% 12.8f", curl_comp * basis->Q, basis->P, basis->curl, stream)); } } return CEED_ERROR_SUCCESS; @@ -1060,18 +1130,18 @@ int CeedBasisView(CeedBasis basis, FILE *stream) { @ref User **/ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { - CeedSize u_length = 0, v_length; - CeedInt dim, num_comp, num_nodes, num_qpts, qdim, cdim; - CeedFESpace fe_space; + CeedSize u_length = 0, v_length; + CeedInt dim, num_comp, q_comp, curl_comp, num_nodes, num_qpts; CeedCall(CeedBasisGetDimension(basis, &dim)); CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); + CeedCall(CeedBasisGetNumQuadratureComponents(basis, &q_comp)); + CeedCall(CeedBasisGetNumCurlComponents(basis, &curl_comp)); CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); CeedCall(CeedVectorGetLength(v, &v_length)); if (u) { CeedCall(CeedVectorGetLength(u, &u_length)); } - CeedCall(CeedBasisGetFESpace(basis, &fe_space)); if (!basis->Apply) { // LCOV_EXCL_START @@ -1090,14 +1160,10 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, // Check vector lengths to prevent out of bounds issues bool bad_dims = false; switch (eval_mode) { - case CEED_EVAL_WEIGHT: - bad_dims = v_length < num_elem * num_qpts; - break; case CEED_EVAL_NONE: case CEED_EVAL_INTERP: - qdim = (fe_space == CEED_FE_SPACE_H1) ? 1 : dim; - bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * qdim || v_length < num_elem * num_comp * num_nodes)) || - (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp * qdim || u_length < num_elem * num_comp * num_nodes))); + bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * q_comp || v_length < num_elem * num_comp * num_nodes)) || + (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp * q_comp || u_length < num_elem * num_comp * num_nodes))); break; case CEED_EVAL_GRAD: bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * dim || v_length < num_elem * num_comp * num_nodes)) || @@ -1108,9 +1174,12 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp || u_length < num_elem * num_comp * num_nodes))); break; case CEED_EVAL_CURL: - cdim = (dim < 3) ? 1 : dim; - bad_dims = ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * cdim || v_length < num_elem * num_comp * num_nodes)) || - (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp * cdim || u_length < num_elem * num_comp * num_nodes))); + bad_dims = + ((t_mode == CEED_TRANSPOSE && (u_length < num_elem * num_comp * num_qpts * curl_comp || v_length < num_elem * num_comp * num_nodes)) || + (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem * num_qpts * num_comp * curl_comp || u_length < num_elem * num_comp * num_nodes))); + break; + case CEED_EVAL_WEIGHT: + bad_dims = v_length < num_elem * num_qpts; break; } if (bad_dims) { @@ -1183,6 +1252,36 @@ int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) { return CEED_ERROR_SUCCESS; } +/** + @brief Get number of Q-vector components for given CeedBasis + + @param[in] basis CeedBasis + @param[out] q_comp Variable to store number of Q-vector components of basis + + @return An error code: 0 - success, otherwise - failure + + @ref Advanced +**/ +int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedInt *q_comp) { + *q_comp = (basis->fe_space == CEED_FE_SPACE_H1) ? 1 : basis->dim; + return CEED_ERROR_SUCCESS; +} + +/** + @brief Get number of curl vector components for given CeedBasis + + @param[in] basis CeedBasis + @param[out] curl_comp Variable to store number of curl vector components of basis + + @return An error code: 0 - success, otherwise - failure + + @ref Advanced +**/ +int CeedBasisGetNumCurlComponents(CeedBasis basis, CeedInt *curl_comp) { + *curl_comp = (basis->dim < 3) ? 1 : basis->dim; + return CEED_ERROR_SUCCESS; +} + /** @brief Get total number of nodes (in dim dimensions) of a CeedBasis @@ -1873,65 +1972,4 @@ CeedPragmaOptimizeOn } CeedPragmaOptimizeOn - /** - @brief Apply Householder Q matrix - - Compute A = Q A, where Q is mxm and A is mxn. - - @param[in,out] A Matrix to apply Householder Q to, in place - @param[in] Q Householder Q matrix - @param[in] tau Householder scaling factors - @param[in] t_mode Transpose mode for application - @param[in] m Number of rows in A - @param[in] n Number of columns in A - @param[in] k Number of elementary reflectors in Q, keval_mode; - CeedInt dim = 1, num_comp = 1, restr_num_comp = 1, size = qf_field->size, qdim = 1, cdim = 1; - CeedFESpace fe_space; + CeedInt dim = 1, num_comp = 1, q_comp = 1, curl_comp = 1, restr_num_comp = 1, size = qf_field->size; // Restriction if (r != CEED_ELEMRESTRICTION_NONE) { @@ -62,9 +61,8 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedEl } CeedCall(CeedBasisGetDimension(b, &dim)); CeedCall(CeedBasisGetNumComponents(b, &num_comp)); - CeedCall(CeedBasisGetFESpace(b, &fe_space)); - qdim = (fe_space == CEED_FE_SPACE_H1) ? 1 : dim; - cdim = (dim < 3) ? 1 : dim; + CeedCall(CeedBasisGetNumQuadratureComponents(b, &q_comp)); + CeedCall(CeedBasisGetNumCurlComponents(b, &curl_comp)); if (r != CEED_ELEMRESTRICTION_NONE && restr_num_comp != num_comp) { // LCOV_EXCL_START return CeedError(ceed, CEED_ERROR_DIMENSION, @@ -91,11 +89,11 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedEl } break; case CEED_EVAL_INTERP: - if (size != num_comp * qdim) { + if (size != num_comp * q_comp) { // LCOV_EXCL_START return CeedError(ceed, CEED_ERROR_DIMENSION, "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction/Basis has %" CeedInt_FMT " components", - qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * qdim); + qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * q_comp); // LCOV_EXCL_STOP } break; @@ -119,11 +117,11 @@ static int CeedOperatorCheckField(Ceed ceed, CeedQFunctionField qf_field, CeedEl } break; case CEED_EVAL_CURL: - if (size != num_comp * cdim) { + if (size != num_comp * curl_comp) { // LCOV_EXCL_START return CeedError(ceed, CEED_ERROR_DIMENSION, "Field '%s' of size %" CeedInt_FMT " and EvalMode %s: ElemRestriction/Basis has %" CeedInt_FMT " components", - qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * cdim); + qf_field->field_name, qf_field->size, CeedEvalModes[qf_field->eval_mode], num_comp * curl_comp); // LCOV_EXCL_STOP } break; diff --git a/interface/ceed-types.c b/interface/ceed-types.c index 4e7f2d8501..fa23a50af8 100644 --- a/interface/ceed-types.c +++ b/interface/ceed-types.c @@ -59,6 +59,7 @@ const char *const CeedContextFieldTypes[] = { }; const char *const CeedFESpaces[] = { - [CEED_FE_SPACE_H1] = "H^1 space", - [CEED_FE_SPACE_HDIV] = "H(div) space", + [CEED_FE_SPACE_H1] = "H^1 space", + [CEED_FE_SPACE_HDIV] = "H(div) space", + [CEED_FE_SPACE_HCURL] = "H(curl) space", }; diff --git a/julia/LibCEED.jl/src/generated/libceed_bindings.jl b/julia/LibCEED.jl/src/generated/libceed_bindings.jl index 887dd5968e..729b1401f2 100644 --- a/julia/LibCEED.jl/src/generated/libceed_bindings.jl +++ b/julia/LibCEED.jl/src/generated/libceed_bindings.jl @@ -389,6 +389,14 @@ function CeedBasisGetNumComponents(basis, num_comp) ccall((:CeedBasisGetNumComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, num_comp) end +function CeedBasisGetNumQuadratureComponents(basis, q_comp) + ccall((:CeedBasisGetNumQuadratureComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, q_comp) +end + +function CeedBasisGetNumCurlComponents(basis, curl_comp) + ccall((:CeedBasisGetNumCurlComponents, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, curl_comp) +end + function CeedBasisGetNumNodes(basis, P) ccall((:CeedBasisGetNumNodes, libceed), Cint, (CeedBasis, Ptr{CeedInt}), basis, P) end @@ -461,14 +469,6 @@ function CeedSimultaneousDiagonalization(ceed, mat_A, mat_B, x, lambda, n) ccall((:CeedSimultaneousDiagonalization, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt), ceed, mat_A, mat_B, x, lambda, n) end -function CeedHouseholderApplyQ(A, Q, tau, t_mode, m, n, k, row, col) - ccall((:CeedHouseholderApplyQ, libceed), Cint, (Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt), A, Q, tau, t_mode, m, n, k, row, col) -end - -function CeedMatrixMatrixMultiply(ceed, mat_A, mat_B, mat_C, m, n, kk) - ccall((:CeedMatrixMatrixMultiply, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt, CeedInt), ceed, mat_A, mat_B, mat_C, m, n, kk) -end - # typedef int ( * CeedQFunctionUser ) ( void * ctx , const CeedInt Q , const CeedScalar * const * in , CeedScalar * const * out ) const CeedQFunctionUser = Ptr{Cvoid} @@ -1338,6 +1338,14 @@ function CeedOperatorSetSetupDone(op) ccall((:CeedOperatorSetSetupDone, libceed), Cint, (CeedOperator,), op) end +function CeedMatrixMatrixMultiply(ceed, mat_A, mat_B, mat_C, m, n, kk) + ccall((:CeedMatrixMatrixMultiply, libceed), Cint, (Ceed, Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedInt, CeedInt, CeedInt), ceed, mat_A, mat_B, mat_C, m, n, kk) +end + +function CeedHouseholderApplyQ(A, Q, tau, t_mode, m, n, k, row, col) + ccall((:CeedHouseholderApplyQ, libceed), Cint, (Ptr{CeedScalar}, Ptr{CeedScalar}, Ptr{CeedScalar}, CeedTransposeMode, CeedInt, CeedInt, CeedInt, CeedInt, CeedInt), A, Q, tau, t_mode, m, n, k, row, col) +end + # Skipping MacroDefinition: CEED_EXTERN extern CEED_VISIBILITY ( default ) # Skipping MacroDefinition: CEED_QFUNCTION_HELPER CEED_QFUNCTION_ATTR static inline