Lines Matching +full:clang +full:- +full:format +full:- +full:16
1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
4 // SPDX-License-Identifier: BSD-2-Clause
8 #include <ceed-impl.h>
31 /// ----------------------------------------------------------------------------
33 /// ----------------------------------------------------------------------------
44 @return An error code: 0 - success, otherwise - failure
51 for (CeedInt i = 2; i < n; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2]; in CeedChebyshevPolynomialsAtPoint()
62 @return An error code: 0 - success, otherwise - failure
76 chebyshev_x[2] = 2 * x * chebyshev_x[1] - chebyshev_x[0]; in CeedChebyshevDerivativeAtPoint()
77 chebyshev_dx[i] = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[1] - chebyshev_dx[i - 2]; in CeedChebyshevDerivativeAtPoint()
85 …Computes \f$A = (I - b v v^T) A\f$, where \f$A\f$ is an \f$m \times n\f$ matrix indexed as `A[i*ro…
95 @return An error code: 0 - success, otherwise - failure
104 A[0 * row + j * col] -= b * w; in CeedHouseholderReflect()
105 for (CeedInt i = 1; i < m; i++) A[i * row + j * col] -= b * w * v[i]; in CeedHouseholderReflect()
118 …@param[in] t_mode @ref CEED_NOTRANSPOSE to rotate the basis counter-clockwise, which has the e…
125 @return An error code: 0 - success, otherwise - failure
142 A[i * stride_ik + j * stride_j] = c * tau1 - s * tau2; in CeedGivensRotation()
152 @param[in] fp_fmt Printing format
159 @return An error code: 0 - success, otherwise - failure
170 fprintf(stream, "%s %-10s", tabs, padded_name); in CeedScalarView()
174 …for (CeedInt j = 0; j < n; j++) fprintf(stream, fp_fmt, fabs(a[i * n + j]) > 1E-14 ? a[i * n + j] … in CeedScalarView()
186 @return An error code: 0 - success, otherwise - failure
200 @return An error code: 0 - success, otherwise - failure
222 @return An error code: 0 - success, otherwise - failure
239 // Check for matching tensor or non-tensor in CeedBasisCreateProjectionMatrices()
305 input_from[m] = (CeedScalar *)&grad_from_source[(m - 1) * Q * P_from]; in CeedBasisCreateProjectionMatrices()
306 output_project[m] = &((*grad_project)[(m - 1) * P_to * P_from]); in CeedBasisCreateProjectionMatrices()
341 @return An error code: 0 - success, otherwise - failure
394 @return An error code: 0 - success, otherwise - failure
472 @return An error code: 0 - success, otherwise - failure
481 // Inserting check because clang-tidy doesn't understand this cannot occur in CeedBasisApplyAtPoints_Core()
501 if (!basis->basis_chebyshev) { in CeedBasisApplyAtPoints_Core()
514 CeedCall(CeedVectorCreate(ceed, num_comp * CeedIntPow(Q_1d, dim), &basis->vec_chebyshev)); in CeedBasisApplyAtPoints_Core()
516 &basis->basis_chebyshev)); in CeedBasisApplyAtPoints_Core()
526 if (!basis->contract) { in CeedBasisApplyAtPoints_Core()
533 // Note - clang-tidy doesn't know basis_ref->contract must be valid here in CeedBasisApplyAtPoints_Core()
534 CeedCheck(basis_ref && basis_ref->contract, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, in CeedBasisApplyAtPoints_Core()
536 CeedCall(CeedTensorContractReferenceCopy(basis_ref->contract, &basis->contract)); in CeedBasisApplyAtPoints_Core()
548 // -- Interpolate to Chebyshev coefficients in CeedBasisApplyAtPoints_Core()
549 …CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, basis->v… in CeedBasisApplyAtPoints_Core()
551 // -- Evaluate Chebyshev polynomials at arbitrary points in CeedBasisApplyAtPoints_Core()
552 CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); in CeedBasisApplyAtPoints_Core()
559 // ---- Values at point in CeedBasisApplyAtPoints_Core()
561 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; in CeedBasisApplyAtPoints_Core()
564 // ------ Tensor contract with current Chebyshev polynomial values in CeedBasisApplyAtPoints_Core()
566 … CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, in CeedBasisApplyAtPoints_Core()
578 // ---- Values at point in CeedBasisApplyAtPoints_Core()
582 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; in CeedBasisApplyAtPoints_Core()
585 // ------ Tensor contract with current Chebyshev polynomial values in CeedBasisApplyAtPoints_Core()
588 … CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, in CeedBasisApplyAtPoints_Core()
602 CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs)); in CeedBasisApplyAtPoints_Core()
613 // -- Transpose of evaluation of Chebyshev polynomials at arbitrary points in CeedBasisApplyAtPoints_Core()
614 CeedCall(CeedVectorGetArrayWrite(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); in CeedBasisApplyAtPoints_Core()
622 // ---- Values at point in CeedBasisApplyAtPoints_Core()
628 // ------ Tensor contract with current Chebyshev polynomial values in CeedBasisApplyAtPoints_Core()
630 …CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 &… in CeedBasisApplyAtPoints_Core()
631 … d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); in CeedBasisApplyAtPoints_Core()
641 // ---- Values at point in CeedBasisApplyAtPoints_Core()
649 // ------ Tensor contract with current Chebyshev polynomial values in CeedBasisApplyAtPoints_Core()
652 … CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, in CeedBasisApplyAtPoints_Core()
653 … (p > 0 || (p == 0 && pass > 0)) && d == (dim - 1), tmp[d % 2], in CeedBasisApplyAtPoints_Core()
654 … d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); in CeedBasisApplyAtPoints_Core()
666 CeedCall(CeedVectorRestoreArray(basis->vec_chebyshev, &chebyshev_coeffs)); in CeedBasisApplyAtPoints_Core()
670 // -- Interpolate transpose from Chebyshev coefficients in CeedBasisApplyAtPoints_Core()
671 … (apply_add) CeedCall(CeedBasisApplyAdd(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTER… in CeedBasisApplyAtPoints_Core()
672 …else CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->v… in CeedBasisApplyAtPoints_Core()
681 /// ----------------------------------------------------------------------------
683 /// ----------------------------------------------------------------------------
688 …@brief Fallback to a reference implementation for a non tensor-product basis for \f$H^1\f$ discret…
698 …@param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal bas…
699 …@param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of no…
704 @return An error code: 0 - success, otherwise - failure
716 CeedCall(CeedReferenceCopy(delegate, &(basis)->obj.ceed)); in CeedBasisCreateH1Fallback()
718 CeedCall(delegate->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, basis)); in CeedBasisCreateH1Fallback()
727 …@param[out] collo_grad_1d Row-major (`Q_1d * Q_1d`) matrix expressing derivatives of basis functio…
729 @return An error code: 0 - success, otherwise - failure
760 …@param[out] chebyshev_interp_1d Row-major (`P_1d * Q_1d`) matrix interpolating from basis nodes to…
762 @return An error code: 0 - success, otherwise - failure
777 // -- Note: Clang-tidy needs this check in CeedBasisGetChebyshevInterp1D()
804 @return An error code: 0 - success, otherwise - failure
809 *is_tensor = basis->is_tensor_basis; in CeedBasisIsTensor()
819 @return An error code: 0 - success, otherwise - failure
824 if (basis->is_tensor_basis && (basis->Q_1d == basis->P_1d)) { in CeedBasisIsCollocated()
827 for (CeedInt i = 0; i < basis->P_1d; i++) { in CeedBasisIsCollocated()
828 …*is_collocated = *is_collocated && (fabs(basis->interp_1d[i + basis->P_1d * i] - 1.0) < 10 * CEED_… in CeedBasisIsCollocated()
829 for (CeedInt j = 0; j < basis->Q_1d; j++) { in CeedBasisIsCollocated()
830 …if (j != i) *is_collocated = *is_collocated && (fabs(basis->interp_1d[j + basis->P_1d * i]) < 10 *… in CeedBasisIsCollocated()
845 @return An error code: 0 - success, otherwise - failure
850 *(void **)data = basis->data; in CeedBasisGetData()
860 @return An error code: 0 - success, otherwise - failure
865 basis->data = data; in CeedBasisSetData()
874 @return An error code: 0 - success, otherwise - failure
884 @brief Get number of Q-vector components for given `CeedBasis`
891 @param[out] q_comp Variable to store number of Q-vector components of basis
893 @return An error code: 0 - success, otherwise - failure
942 …edBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Can only evaluate tensor-product bases at poin… in CeedBasisGetFlopsEstimate()
954 CeedInt tensor_flops = 0, pre = num_comp * CeedIntPow(P_1d, dim - 1), post = 1; in CeedBasisGetFlopsEstimate()
971 CeedInt chebyshev_flops = (Q_1d - 2) * 3 + 1, d_chebyshev_flops = (Q_1d - 2) * 8 + 1; in CeedBasisGetFlopsEstimate()
972 CeedInt point_tensor_flops = 0, pre = CeedIntPow(Q_1d, dim - 1), post = 1; in CeedBasisGetFlopsEstimate()
998 …dim * (2 * Q_1d * Q_1d + (t_mode == CEED_TRANSPOSE ? 2 : 3) * Q_1d) + (dim - 1) * (2 * chebyshev_f… in CeedBasisGetFlopsEstimate()
1002 …*flops += num_points * (is_gpu ? num_comp : 1) * dim * (d_chebyshev_flops + (dim - 1) * chebyshev_… in CeedBasisGetFlopsEstimate()
1074 @return An error code: 0 - success, otherwise - failure
1079 *fe_space = basis->fe_space; in CeedBasisGetFESpace()
1089 @return An error code: 0 - success, otherwise - failure
1094 *dim = (CeedInt)topo >> 16; in CeedBasisGetTopologyDimension()
1104 @return An error code: 0 - success, otherwise - failure
1109 *contract = basis->contract; in CeedBasisGetTensorContract()
1119 @return An error code: 0 - success, otherwise - failure
1124 basis->contract = contract; in CeedBasisSetTensorContract()
1135 @param[in] mat_A Row-major matrix `A`
1136 @param[in] mat_B Row-major matrix `B`
1137 @param[out] mat_C Row-major output matrix `C`
1142 @return An error code: 0 - success, otherwise - failure
1162 @param[in,out] mat Row-major matrix to be factorized in place
1167 @return An error code: 0 - success, otherwise - failure
1180 if (i >= m - 1) { // last row of matrix, no reflection needed in CeedQRFactorization()
1191 const CeedScalar R_ii = -copysign(norm, v[i]); in CeedQRFactorization()
1193 v[i] -= R_ii; in CeedQRFactorization()
1201 CeedHouseholderReflect(&mat[i * n + i + 1], &v[i], tau[i], m - i, n - i - 1, n, 1); in CeedQRFactorization()
1224 @return An error code: 0 - success, otherwise - failure
1234 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k - 1 - ii; in CeedHouseholderApplyQ()
1236 // Apply Householder reflector (I - tau v v^T) collo_grad_1d^T in CeedHouseholderApplyQ()
1237 CeedCall(CeedHouseholderReflect(&mat_A[i * row], &v[i], tau[i], m - i, n, row, col)); in CeedHouseholderApplyQ()
1247 @param[in] mat Row-major matrix to compute pseudoinverse of
1250 @param[out] mat_pinv Row-major pseudoinverse matrix
1252 @return An error code: 0 - success, otherwise - failure
1267 // -- Apply Q^T, I = Q^T * I in CeedMatrixPseudoinverse()
1270 // -- Apply R_inv, mat_pinv = R_inv * Q^T in CeedMatrixPseudoinverse()
1272 mat_pinv[j + m * (n - 1)] = I[j + m * (n - 1)] / mat_copy[n * n - 1]; in CeedMatrixPseudoinverse()
1273 for (CeedInt i = n - 2; i >= 0; i--) { // Row i in CeedMatrixPseudoinverse()
1275 …for (CeedInt k = i + 1; k < n; k++) mat_pinv[j + m * i] -= mat_copy[k + n * i] * mat_pinv[j + m * … in CeedMatrixPseudoinverse()
1291 @param[in,out] mat Row-major matrix to be factorized in place
1295 @return An error code: 0 - success, otherwise - failure
1301 // Check bounds for clang-tidy in CeedSymmetricSchurDecomposition()
1304 CeedScalar v[n - 1], tau[n - 1], mat_T[n * n]; in CeedSymmetricSchurDecomposition()
1313 for (CeedInt i = 0; i < n - 1; i++) { in CeedSymmetricSchurDecomposition()
1318 for (CeedInt j = i + 1; j < n - 1; j++) { in CeedSymmetricSchurDecomposition()
1322 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:n-1] in CeedSymmetricSchurDecomposition()
1323 const CeedScalar R_ii = -copysign(norm, v[i]); in CeedSymmetricSchurDecomposition()
1325 v[i] -= R_ii; in CeedSymmetricSchurDecomposition()
1329 tau[i] = i == n - 2 ? 2 : 2 * v[i] * v[i] / (v[i] * v[i] + sigma); in CeedSymmetricSchurDecomposition()
1330 for (CeedInt j = i + 1; j < n - 1; j++) v[j] /= v[i]; in CeedSymmetricSchurDecomposition()
1338 …CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, … in CeedSymmetricSchurDecomposition()
1339 …CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), 1, … in CeedSymmetricSchurDecomposition()
1344 for (CeedInt j = i + 1; j < n - 1; j++) { in CeedSymmetricSchurDecomposition()
1349 for (CeedInt i = n - 2; i >= 0; i--) { in CeedSymmetricSchurDecomposition()
1352 for (CeedInt j = i + 1; j < n - 1; j++) { in CeedSymmetricSchurDecomposition()
1356 …CeedHouseholderReflect(&mat[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1); in CeedSymmetricSchurDecomposition()
1368 for (CeedInt i = n - 2; i >= 0; i--) { in CeedSymmetricSchurDecomposition()
1372 for (CeedInt i = 0; i < n - q - 1; i++) { in CeedSymmetricSchurDecomposition()
1376 if (q == n - 1) break; // Finished reducing in CeedSymmetricSchurDecomposition()
1379 …CeedScalar t_nn = mat_T[(n - 1 - q) + n * (n - 1 - q)], t_nnm1 = mat_T[(n - 2 - q) + n * (n - 1 - … in CeedSymmetricSchurDecomposition()
1380 CeedScalar d = (mat_T[(n - 2 - q) + n * (n - 2 - q)] - t_nn) / 2; in CeedSymmetricSchurDecomposition()
1381 CeedScalar mu = t_nn - t_nnm1 * t_nnm1 / (d + copysign(sqrt(d * d + t_nnm1 * t_nnm1), d)); in CeedSymmetricSchurDecomposition()
1382 CeedScalar x = mat_T[p + n * p] - mu; in CeedSymmetricSchurDecomposition()
1385 for (CeedInt k = p; k < n - q - 1; k++) { in CeedSymmetricSchurDecomposition()
1391 const CeedScalar tau = -x / z; in CeedSymmetricSchurDecomposition()
1396 const CeedScalar tau = -z / x; in CeedSymmetricSchurDecomposition()
1411 if (k < n - q - 2) { in CeedSymmetricSchurDecomposition()
1436 @param[in] mat_A Row-major matrix to be factorized with eigenvalues
1437 @param[in] mat_B Row-major matrix to be factorized to identity
1438 @param[out] mat_X Row-major orthogonal matrix
1442 @return An error code: 0 - success, otherwise - failure
1459 for (CeedInt i = n - 1; i >= 0; i--) { in CeedSimultaneousDiagonalization()
1468 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T in CeedSimultaneousDiagonalization()
1469 // = D^-1/2 G^T A G D^-1/2 in CeedSimultaneousDiagonalization()
1470 // -- D = D^-1/2 in CeedSimultaneousDiagonalization()
1472 // -- G = G D^-1/2 in CeedSimultaneousDiagonalization()
1473 // -- C = D^-1/2 G^T in CeedSimultaneousDiagonalization()
1480 // -- X = (D^-1/2 G^T) A in CeedSimultaneousDiagonalization()
1482 // -- C = (D^-1/2 G^T A) (G D^-1/2) in CeedSimultaneousDiagonalization()
1489 for (CeedInt i = n - 1; i >= 0; i--) { in CeedSimultaneousDiagonalization()
1498 // Set X = (G D^1/2)^-T Q in CeedSimultaneousDiagonalization()
1499 // = G D^-1/2 Q in CeedSimultaneousDiagonalization()
1512 /// ----------------------------------------------------------------------------
1514 /// ----------------------------------------------------------------------------
1519 @brief Create a tensor-product basis for \f$H^1\f$ discretizations
1526 …@param[in] interp_1d Row-major (`Q_1d * P_1d`) matrix expressing the values of nodal basis func…
1527 …@param[in] grad_1d Row-major (`Q_1d * P_1d`) matrix expressing derivatives of nodal basis fun…
1528 …y of length `Q_1d` holding the locations of quadrature points on the 1D reference element `[-1, 1]`
1532 @return An error code: 0 - success, otherwise - failure
1538 if (!ceed->BasisCreateTensorH1) { in CeedBasisCreateTensorH1()
1556 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj)); in CeedBasisCreateTensorH1()
1557 (*basis)->is_tensor_basis = true; in CeedBasisCreateTensorH1()
1558 (*basis)->dim = dim; in CeedBasisCreateTensorH1()
1559 (*basis)->topo = topo; in CeedBasisCreateTensorH1()
1560 (*basis)->num_comp = num_comp; in CeedBasisCreateTensorH1()
1561 (*basis)->P_1d = P_1d; in CeedBasisCreateTensorH1()
1562 (*basis)->Q_1d = Q_1d; in CeedBasisCreateTensorH1()
1563 (*basis)->P = CeedIntPow(P_1d, dim); in CeedBasisCreateTensorH1()
1564 (*basis)->Q = CeedIntPow(Q_1d, dim); in CeedBasisCreateTensorH1()
1565 (*basis)->fe_space = CEED_FE_SPACE_H1; in CeedBasisCreateTensorH1()
1566 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_ref_1d)); in CeedBasisCreateTensorH1()
1567 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_weight_1d)); in CeedBasisCreateTensorH1()
1568 if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d * sizeof(q_ref_1d[0])); in CeedBasisCreateTensorH1()
1569 if (q_weight_1d) memcpy((*basis)->q_weight_1d, q_weight_1d, Q_1d * sizeof(q_weight_1d[0])); in CeedBasisCreateTensorH1()
1570 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->interp_1d)); in CeedBasisCreateTensorH1()
1571 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->grad_1d)); in CeedBasisCreateTensorH1()
1572 if (interp_1d) memcpy((*basis)->interp_1d, interp_1d, Q_1d * P_1d * sizeof(interp_1d[0])); in CeedBasisCreateTensorH1()
1573 if (grad_1d) memcpy((*basis)->grad_1d, grad_1d, Q_1d * P_1d * sizeof(grad_1d[0])); in CeedBasisCreateTensorH1()
1574 …CeedCall(ceed->BasisCreateTensorH1(dim, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, *ba… in CeedBasisCreateTensorH1()
1579 @brief Create a tensor-product \f$H^1\f$ Lagrange basis
1584 @param[in] P Number of Gauss-Lobatto nodes in one dimension.
1585 The polynomial degree of the resulting `Q_k` element is `k = P - 1`.
1590 @return An error code: 0 - success, otherwise - failure
1625 c3 = nodes[0] - q_ref_1d[i]; in CeedBasisCreateTensorH1Lagrange()
1630 c3 = nodes[j] - q_ref_1d[i]; in CeedBasisCreateTensorH1Lagrange()
1632 dx = nodes[j] - nodes[k]; in CeedBasisCreateTensorH1Lagrange()
1634 if (k == j - 1) { in CeedBasisCreateTensorH1Lagrange()
1635 grad_1d[i * P + j] = c1 * (interp_1d[i * P + k] - c4 * grad_1d[i * P + k]) / c2; in CeedBasisCreateTensorH1Lagrange()
1636 interp_1d[i * P + j] = -c1 * c4 * interp_1d[i * P + k] / c2; in CeedBasisCreateTensorH1Lagrange()
1638 grad_1d[i * P + k] = (c3 * grad_1d[i * P + k] - interp_1d[i * P + k]) / dx; in CeedBasisCreateTensorH1Lagrange()
1656 @brief Create a non tensor-product basis for \f$H^1\f$ discretizations
1663 …@param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal bas…
1664 …@param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of no…
1669 @return An error code: 0 - success, otherwise - failure
1677 if (!ceed->BasisCreateH1) { in CeedBasisCreateH1()
1694 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj)); in CeedBasisCreateH1()
1695 (*basis)->is_tensor_basis = false; in CeedBasisCreateH1()
1696 (*basis)->dim = dim; in CeedBasisCreateH1()
1697 (*basis)->topo = topo; in CeedBasisCreateH1()
1698 (*basis)->num_comp = num_comp; in CeedBasisCreateH1()
1699 (*basis)->P = P; in CeedBasisCreateH1()
1700 (*basis)->Q = Q; in CeedBasisCreateH1()
1701 (*basis)->fe_space = CEED_FE_SPACE_H1; in CeedBasisCreateH1()
1702 CeedCall(CeedCalloc(Q * dim, &(*basis)->q_ref_1d)); in CeedBasisCreateH1()
1703 CeedCall(CeedCalloc(Q, &(*basis)->q_weight_1d)); in CeedBasisCreateH1()
1704 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); in CeedBasisCreateH1()
1705 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); in CeedBasisCreateH1()
1706 CeedCall(CeedCalloc(Q * P, &(*basis)->interp)); in CeedBasisCreateH1()
1707 CeedCall(CeedCalloc(dim * Q * P, &(*basis)->grad)); in CeedBasisCreateH1()
1708 if (interp) memcpy((*basis)->interp, interp, Q * P * sizeof(interp[0])); in CeedBasisCreateH1()
1709 if (grad) memcpy((*basis)->grad, grad, dim * Q * P * sizeof(grad[0])); in CeedBasisCreateH1()
1710 CeedCall(ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, *basis)); in CeedBasisCreateH1()
1715 @brief Create a non tensor-product basis for \f$H(\mathrm{div})\f$ discretizations
1722 …@param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of bas…
1723 …@param[in] div Row-major (`num_qpts * num_nodes`) matrix expressing divergence of basis fun…
1728 @return An error code: 0 - success, otherwise - failure
1736 if (!ceed->BasisCreateHdiv) { in CeedBasisCreateHdiv()
1753 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj)); in CeedBasisCreateHdiv()
1754 (*basis)->is_tensor_basis = false; in CeedBasisCreateHdiv()
1755 (*basis)->dim = dim; in CeedBasisCreateHdiv()
1756 (*basis)->topo = topo; in CeedBasisCreateHdiv()
1757 (*basis)->num_comp = num_comp; in CeedBasisCreateHdiv()
1758 (*basis)->P = P; in CeedBasisCreateHdiv()
1759 (*basis)->Q = Q; in CeedBasisCreateHdiv()
1760 (*basis)->fe_space = CEED_FE_SPACE_HDIV; in CeedBasisCreateHdiv()
1761 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); in CeedBasisCreateHdiv()
1762 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); in CeedBasisCreateHdiv()
1763 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); in CeedBasisCreateHdiv()
1764 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); in CeedBasisCreateHdiv()
1765 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); in CeedBasisCreateHdiv()
1766 CeedCall(CeedMalloc(Q * P, &(*basis)->div)); in CeedBasisCreateHdiv()
1767 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); in CeedBasisCreateHdiv()
1768 if (div) memcpy((*basis)->div, div, Q * P * sizeof(div[0])); in CeedBasisCreateHdiv()
1769 CeedCall(ceed->BasisCreateHdiv(topo, dim, P, Q, interp, div, q_ref, q_weight, *basis)); in CeedBasisCreateHdiv()
1774 @brief Create a non tensor-product basis for \f$H(\mathrm{curl})\f$ discretizations
1781 …@param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of bas…
1782 …@param[in] curl Row-major (`curl_comp * num_qpts * num_nodes`, `curl_comp = 1` if `dim < 3` …
1787 @return An error code: 0 - success, otherwise - failure
1795 if (!ceed->BasisCreateHcurl) { in CeedBasisCreateHcurl()
1813 CeedCall(CeedObjectCreate(ceed, CeedBasisView_Object, CeedBasisDestroy_Object, &(*basis)->obj)); in CeedBasisCreateHcurl()
1814 (*basis)->is_tensor_basis = false; in CeedBasisCreateHcurl()
1815 (*basis)->dim = dim; in CeedBasisCreateHcurl()
1816 (*basis)->topo = topo; in CeedBasisCreateHcurl()
1817 (*basis)->num_comp = num_comp; in CeedBasisCreateHcurl()
1818 (*basis)->P = P; in CeedBasisCreateHcurl()
1819 (*basis)->Q = Q; in CeedBasisCreateHcurl()
1820 (*basis)->fe_space = CEED_FE_SPACE_HCURL; in CeedBasisCreateHcurl()
1821 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); in CeedBasisCreateHcurl()
1822 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); in CeedBasisCreateHcurl()
1823 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); in CeedBasisCreateHcurl()
1824 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); in CeedBasisCreateHcurl()
1825 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); in CeedBasisCreateHcurl()
1826 CeedCall(CeedMalloc(curl_comp * Q * P, &(*basis)->curl)); in CeedBasisCreateHcurl()
1827 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); in CeedBasisCreateHcurl()
1828 if (curl) memcpy((*basis)->curl, curl, curl_comp * Q * P * sizeof(curl[0])); in CeedBasisCreateHcurl()
1829 CeedCall(ceed->BasisCreateHcurl(topo, dim, P, Q, interp, curl, q_ref, q_weight, *basis)); in CeedBasisCreateHcurl()
1846 …Note: If either `basis_from` or `basis_to` are non-tensor, then `basis_project` will also be non-t…
1852 @return An error code: 0 - success, otherwise - failure
1904 …Note: If the value of `*basis_copy` passed into this function is non-`NULL`, then it is assumed th…
1910 @return An error code: 0 - success, otherwise - failure
1927 @return Error code: 0 - success, otherwise - failure
1942 @return Error code: 0 - success, otherwise - failure
1957 @return An error code: 0 - success, otherwise - failure
1983 …tream, "%s P: %" CeedInt_FMT "\n%s Q: %" CeedInt_FMT "\n", tabs, basis->P_1d, tabs, basis->Q_1d); in CeedBasisView()
1985 …intf(stream, "%s P: %" CeedInt_FMT "\n%s Q: %" CeedInt_FMT "\n", tabs, basis->P, tabs, basis->Q); in CeedBasisView()
1987 …CeedInt_FMT "\n%s field components: %" CeedInt_FMT "\n", tabs, basis->dim, tabs, basis->num_comp); in CeedBasisView()
2004 } else { // non-tensor basis in CeedBasisView()
2056 @return An error code: 0 - success, otherwise - failure
2062 …CeedCheck(basis->Apply, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not supp… in CeedBasisApply()
2063 CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v)); in CeedBasisApply()
2084 @return An error code: 0 - success, otherwise - failure
2091 …CeedCheck(basis->ApplyAdd, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not i… in CeedBasisApplyAdd()
2092 CeedCall(basis->ApplyAdd(basis, num_elem, t_mode, eval_mode, u, v)); in CeedBasisApplyAdd()
2112 @return An error code: 0 - success, otherwise - failure
2119 if (basis->ApplyAtPoints) { in CeedBasisApplyAtPoints()
2120 CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); in CeedBasisApplyAtPoints()
2143 @return An error code: 0 - success, otherwise - failure
2151 if (basis->ApplyAddAtPoints) { in CeedBasisApplyAddAtPoints()
2152 CeedCall(basis->ApplyAddAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); in CeedBasisApplyAddAtPoints()
2165 @return An error code: 0 - success, otherwise - failure
2191 @return An error code: 0 - success, otherwise - failure
2196 *dim = basis->dim; in CeedBasisGetDimension()
2206 @return An error code: 0 - success, otherwise - failure
2211 *topo = basis->topo; in CeedBasisGetTopology()
2221 @return An error code: 0 - success, otherwise - failure
2226 *num_comp = basis->num_comp; in CeedBasisGetNumComponents()
2236 @return An error code: 0 - success, otherwise - failure
2241 *P = basis->P; in CeedBasisGetNumNodes()
2251 @return An error code: 0 - success, otherwise - failure
2256 …CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply P_1… in CeedBasisGetNumNodes1D()
2257 *P_1d = basis->P_1d; in CeedBasisGetNumNodes1D()
2267 @return An error code: 0 - success, otherwise - failure
2272 *Q = basis->Q; in CeedBasisGetNumQuadraturePoints()
2282 @return An error code: 0 - success, otherwise - failure
2287 …CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply Q_1… in CeedBasisGetNumQuadraturePoints1D()
2288 *Q_1d = basis->Q_1d; in CeedBasisGetNumQuadraturePoints1D()
2298 @return An error code: 0 - success, otherwise - failure
2303 *q_ref = basis->q_ref_1d; in CeedBasisGetQRef()
2313 @return An error code: 0 - success, otherwise - failure
2318 *q_weight = basis->q_weight_1d; in CeedBasisGetQWeights()
2328 @return An error code: 0 - success, otherwise - failure
2333 if (!basis->interp && basis->is_tensor_basis) { in CeedBasisGetInterp()
2335 CeedCall(CeedMalloc(basis->Q * basis->P, &basis->interp)); in CeedBasisGetInterp()
2338 for (CeedInt i = 0; i < basis->Q * basis->P; i++) basis->interp[i] = 1.0; in CeedBasisGetInterp()
2341 for (CeedInt d = 0; d < basis->dim; d++) { in CeedBasisGetInterp()
2342 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) { in CeedBasisGetInterp()
2343 for (CeedInt node = 0; node < basis->P; node++) { in CeedBasisGetInterp()
2344 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; in CeedBasisGetInterp()
2345 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; in CeedBasisGetInterp()
2347 basis->interp[qpt * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p]; in CeedBasisGetInterp()
2352 *interp = basis->interp; in CeedBasisGetInterp()
2362 @return An error code: 0 - success, otherwise - failure
2371 *interp_1d = basis->interp_1d; in CeedBasisGetInterp1D()
2381 @return An error code: 0 - success, otherwise - failure
2386 if (!basis->grad && basis->is_tensor_basis) { in CeedBasisGetGrad()
2388 CeedCall(CeedMalloc(basis->dim * basis->Q * basis->P, &basis->grad)); in CeedBasisGetGrad()
2391 for (CeedInt i = 0; i < basis->dim * basis->Q * basis->P; i++) basis->grad[i] = 1.0; in CeedBasisGetGrad()
2394 for (CeedInt d = 0; d < basis->dim; d++) { in CeedBasisGetGrad()
2395 for (CeedInt i = 0; i < basis->dim; i++) { in CeedBasisGetGrad()
2396 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) { in CeedBasisGetGrad()
2397 for (CeedInt node = 0; node < basis->P; node++) { in CeedBasisGetGrad()
2398 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; in CeedBasisGetGrad()
2399 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; in CeedBasisGetGrad()
2401 …if (i == d) basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->grad_1d[q * basis->P_1… in CeedBasisGetGrad()
2402 …else basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p… in CeedBasisGetGrad()
2408 *grad = basis->grad; in CeedBasisGetGrad()
2418 @return An error code: 0 - success, otherwise - failure
2427 *grad_1d = basis->grad_1d; in CeedBasisGetGrad1D()
2437 @return An error code: 0 - success, otherwise - failure
2442 *div = basis->div; in CeedBasisGetDiv()
2452 @return An error code: 0 - success, otherwise - failure
2457 *curl = basis->curl; in CeedBasisGetCurl()
2466 @return An error code: 0 - success, otherwise - failure
2475 if ((*basis)->Destroy) CeedCall((*basis)->Destroy(*basis)); in CeedBasisDestroy()
2476 CeedCall(CeedTensorContractDestroy(&(*basis)->contract)); in CeedBasisDestroy()
2477 CeedCall(CeedFree(&(*basis)->q_ref_1d)); in CeedBasisDestroy()
2478 CeedCall(CeedFree(&(*basis)->q_weight_1d)); in CeedBasisDestroy()
2479 CeedCall(CeedFree(&(*basis)->interp)); in CeedBasisDestroy()
2480 CeedCall(CeedFree(&(*basis)->interp_1d)); in CeedBasisDestroy()
2481 CeedCall(CeedFree(&(*basis)->grad)); in CeedBasisDestroy()
2482 CeedCall(CeedFree(&(*basis)->grad_1d)); in CeedBasisDestroy()
2483 CeedCall(CeedFree(&(*basis)->div)); in CeedBasisDestroy()
2484 CeedCall(CeedFree(&(*basis)->curl)); in CeedBasisDestroy()
2485 CeedCall(CeedVectorDestroy(&(*basis)->vec_chebyshev)); in CeedBasisDestroy()
2486 CeedCall(CeedBasisDestroy(&(*basis)->basis_chebyshev)); in CeedBasisDestroy()
2487 CeedCall(CeedObjectDestroy_Private(&(*basis)->obj)); in CeedBasisDestroy()
2493 @brief Construct a Gauss-Legendre quadrature
2495 …@param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-1` exac…
2496 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]`
2499 @return An error code: 0 - success, otherwise - failure
2515 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); in CeedGaussQuadrature()
2520 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); in CeedGaussQuadrature()
2521 xi = xi - P2 / dP2; in CeedGaussQuadrature()
2527 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); in CeedGaussQuadrature()
2531 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); in CeedGaussQuadrature()
2532 xi = xi - P2 / dP2; in CeedGaussQuadrature()
2535 wi = 2.0 / ((1.0 - xi * xi) * dP2 * dP2); in CeedGaussQuadrature()
2537 q_weight_1d[Q - 1 - i] = wi; in CeedGaussQuadrature()
2538 q_ref_1d[i] = -xi; in CeedGaussQuadrature()
2539 q_ref_1d[Q - 1 - i] = xi; in CeedGaussQuadrature()
2545 @brief Construct a Gauss-Legendre-Lobatto quadrature
2547 …@param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-3` exac…
2548 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]`
2551 @return An error code: 0 - success, otherwise - failure
2561 wi = 2.0 / ((CeedScalar)(Q * (Q - 1))); in CeedLobattoQuadrature()
2564 q_weight_1d[Q - 1] = wi; in CeedLobattoQuadrature()
2566 q_ref_1d[0] = -1.0; in CeedLobattoQuadrature()
2567 q_ref_1d[Q - 1] = 1.0; in CeedLobattoQuadrature()
2569 for (CeedInt i = 1; i <= (Q - 1) / 2; i++) { in CeedLobattoQuadrature()
2571 xi = cos(PI * (CeedScalar)(i) / (CeedScalar)(Q - 1)); in CeedLobattoQuadrature()
2577 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); in CeedLobattoQuadrature()
2582 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); in CeedLobattoQuadrature()
2583 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi); in CeedLobattoQuadrature()
2584 xi = xi - dP2 / d2P2; in CeedLobattoQuadrature()
2590 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); in CeedLobattoQuadrature()
2594 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); in CeedLobattoQuadrature()
2595 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi); in CeedLobattoQuadrature()
2596 xi = xi - dP2 / d2P2; in CeedLobattoQuadrature()
2599 wi = 2.0 / (((CeedScalar)(Q * (Q - 1))) * P2 * P2); in CeedLobattoQuadrature()
2602 q_weight_1d[Q - 1 - i] = wi; in CeedLobattoQuadrature()
2604 q_ref_1d[i] = -xi; in CeedLobattoQuadrature()
2605 q_ref_1d[Q - 1 - i] = xi; in CeedLobattoQuadrature()