1 // Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed-impl.h> 9 #include <ceed.h> 10 #include <ceed/backend.h> 11 #include <math.h> 12 #include <stdbool.h> 13 #include <stdio.h> 14 #include <string.h> 15 16 /// @file 17 /// Implementation of CeedBasis interfaces 18 19 /// @cond DOXYGEN_SKIP 20 static struct CeedBasis_private ceed_basis_none; 21 /// @endcond 22 23 /// @addtogroup CeedBasisUser 24 /// @{ 25 26 /// Argument for @ref CeedOperatorSetField() indicating that the field does not require a `CeedBasis` 27 const CeedBasis CEED_BASIS_NONE = &ceed_basis_none; 28 29 /// @} 30 31 /// ---------------------------------------------------------------------------- 32 /// CeedBasis Library Internal Functions 33 /// ---------------------------------------------------------------------------- 34 /// @addtogroup CeedBasisDeveloper 35 /// @{ 36 37 /** 38 @brief Compute Chebyshev polynomial values at a point 39 40 @param[in] x Coordinate to evaluate Chebyshev polynomials at 41 @param[in] n Number of Chebyshev polynomials to evaluate, `n >= 2` 42 @param[out] chebyshev_x Array of Chebyshev polynomial values 43 44 @return An error code: 0 - success, otherwise - failure 45 46 @ref Developer 47 **/ 48 static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_x) { 49 chebyshev_x[0] = 1.0; 50 chebyshev_x[1] = 2 * x; 51 for (CeedInt i = 2; i < n; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2]; 52 return CEED_ERROR_SUCCESS; 53 } 54 55 /** 56 @brief Compute values of the derivative of Chebyshev polynomials at a point 57 58 @param[in] x Coordinate to evaluate derivative of Chebyshev polynomials at 59 @param[in] n Number of Chebyshev polynomials to evaluate, `n >= 2` 60 @param[out] chebyshev_dx Array of Chebyshev polynomial derivative values 61 62 @return An error code: 0 - success, otherwise - failure 63 64 @ref Developer 65 **/ 66 static int CeedChebyshevDerivativeAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_dx) { 67 CeedScalar chebyshev_x[3]; 68 69 chebyshev_x[1] = 1.0; 70 chebyshev_x[2] = 2 * x; 71 chebyshev_dx[0] = 0.0; 72 chebyshev_dx[1] = 2.0; 73 for (CeedInt i = 2; i < n; i++) { 74 chebyshev_x[0] = chebyshev_x[1]; 75 chebyshev_x[1] = chebyshev_x[2]; 76 chebyshev_x[2] = 2 * x * chebyshev_x[1] - chebyshev_x[0]; 77 chebyshev_dx[i] = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[1] - chebyshev_dx[i - 2]; 78 } 79 return CEED_ERROR_SUCCESS; 80 } 81 82 /** 83 @brief Compute Householder reflection. 84 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*row + j*col]`. 86 87 @param[in,out] A Matrix to apply Householder reflection to, in place 88 @param[in] v Householder vector 89 @param[in] b Scaling factor 90 @param[in] m Number of rows in `A` 91 @param[in] n Number of columns in `A` 92 @param[in] row Row stride 93 @param[in] col Col stride 94 95 @return An error code: 0 - success, otherwise - failure 96 97 @ref Developer 98 **/ 99 static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v, CeedScalar b, CeedInt m, CeedInt n, CeedInt row, CeedInt col) { 100 for (CeedInt j = 0; j < n; j++) { 101 CeedScalar w = A[0 * row + j * col]; 102 103 for (CeedInt i = 1; i < m; i++) w += v[i] * A[i * row + j * col]; 104 A[0 * row + j * col] -= b * w; 105 for (CeedInt i = 1; i < m; i++) A[i * row + j * col] -= b * w * v[i]; 106 } 107 return CEED_ERROR_SUCCESS; 108 } 109 110 /** 111 @brief Compute Givens rotation 112 113 Computes \f$A = G A\f$ (or \f$G^T A\f$ in transpose mode), where \f$A\f$ is an \f$m \times n\f$ matrix indexed as `A[i*n + j*m]`. 114 115 @param[in,out] A Row major matrix to apply Givens rotation to, in place 116 @param[in] c Cosine factor 117 @param[in] s Sine factor 118 @param[in] t_mode @ref CEED_NOTRANSPOSE to rotate the basis counter-clockwise, which has the effect of rotating columns of `A` clockwise; 119 @ref CEED_TRANSPOSE for the opposite rotation 120 @param[in] i First row/column to apply rotation 121 @param[in] k Second row/column to apply rotation 122 @param[in] m Number of rows in `A` 123 @param[in] n Number of columns in `A` 124 125 @return An error code: 0 - success, otherwise - failure 126 127 @ref Developer 128 **/ 129 static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s, CeedTransposeMode t_mode, CeedInt i, CeedInt k, CeedInt m, CeedInt n) { 130 CeedInt stride_j = 1, stride_ik = m, num_its = n; 131 132 if (t_mode == CEED_NOTRANSPOSE) { 133 stride_j = n; 134 stride_ik = 1; 135 num_its = m; 136 } 137 138 // Apply rotation 139 for (CeedInt j = 0; j < num_its; j++) { 140 CeedScalar tau1 = A[i * stride_ik + j * stride_j], tau2 = A[k * stride_ik + j * stride_j]; 141 142 A[i * stride_ik + j * stride_j] = c * tau1 - s * tau2; 143 A[k * stride_ik + j * stride_j] = s * tau1 + c * tau2; 144 } 145 return CEED_ERROR_SUCCESS; 146 } 147 148 /** 149 @brief View an array stored in a `CeedBasis` 150 151 @param[in] name Name of array 152 @param[in] fp_fmt Printing format 153 @param[in] m Number of rows in array 154 @param[in] n Number of columns in array 155 @param[in] a Array to be viewed 156 @param[in] stream Stream to view to, e.g., `stdout` 157 158 @return An error code: 0 - success, otherwise - failure 159 160 @ref Developer 161 **/ 162 static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m, CeedInt n, const CeedScalar *a, FILE *stream) { 163 if (m > 1) { 164 fprintf(stream, " %s:\n", name); 165 } else { 166 char padded_name[12]; 167 168 snprintf(padded_name, 11, "%s:", name); 169 fprintf(stream, " %-10s", padded_name); 170 } 171 for (CeedInt i = 0; i < m; i++) { 172 if (m > 1) fprintf(stream, " [%" CeedInt_FMT "]", i); 173 for (CeedInt j = 0; j < n; j++) fprintf(stream, fp_fmt, fabs(a[i * n + j]) > 1E-14 ? a[i * n + j] : 0); 174 fputs("\n", stream); 175 } 176 return CEED_ERROR_SUCCESS; 177 } 178 179 /** 180 @brief Create the interpolation and gradient matrices for projection from the nodes of `basis_from` to the nodes of `basis_to`. 181 182 The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization. 183 The gradient is given by `grad_project = interp_to^+ * grad_from`, and is only computed for \f$H^1\f$ spaces otherwise it should not be used. 184 185 Note: `basis_from` and `basis_to` must have compatible quadrature spaces. 186 187 @param[in] basis_from `CeedBasis` to project from 188 @param[in] basis_to `CeedBasis` to project to 189 @param[out] interp_project Address of the variable where the newly created interpolation matrix will be stored 190 @param[out] grad_project Address of the variable where the newly created gradient matrix will be stored 191 192 @return An error code: 0 - success, otherwise - failure 193 194 @ref Developer 195 **/ 196 static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis basis_to, CeedScalar **interp_project, CeedScalar **grad_project) { 197 bool are_both_tensor; 198 CeedInt Q, Q_to, Q_from, P_to, P_from; 199 200 // Check for compatible quadrature spaces 201 CeedCall(CeedBasisGetNumQuadraturePoints(basis_to, &Q_to)); 202 CeedCall(CeedBasisGetNumQuadraturePoints(basis_from, &Q_from)); 203 CeedCheck(Q_to == Q_from, CeedBasisReturnCeed(basis_to), CEED_ERROR_DIMENSION, 204 "Bases must have compatible quadrature spaces." 205 " 'basis_from' has %" CeedInt_FMT " points and 'basis_to' has %" CeedInt_FMT, 206 Q_from, Q_to); 207 Q = Q_to; 208 209 // Check for matching tensor or non-tensor 210 { 211 bool is_tensor_to, is_tensor_from; 212 213 CeedCall(CeedBasisIsTensor(basis_to, &is_tensor_to)); 214 CeedCall(CeedBasisIsTensor(basis_from, &is_tensor_from)); 215 are_both_tensor = is_tensor_to && is_tensor_from; 216 } 217 if (are_both_tensor) { 218 CeedCall(CeedBasisGetNumNodes1D(basis_to, &P_to)); 219 CeedCall(CeedBasisGetNumNodes1D(basis_from, &P_from)); 220 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis_from, &Q)); 221 } else { 222 CeedCall(CeedBasisGetNumNodes(basis_to, &P_to)); 223 CeedCall(CeedBasisGetNumNodes(basis_from, &P_from)); 224 } 225 226 // Check for matching FE space 227 CeedFESpace fe_space_to, fe_space_from; 228 229 CeedCall(CeedBasisGetFESpace(basis_to, &fe_space_to)); 230 CeedCall(CeedBasisGetFESpace(basis_from, &fe_space_from)); 231 CeedCheck(fe_space_to == fe_space_from, CeedBasisReturnCeed(basis_to), CEED_ERROR_MINOR, 232 "Bases must both be the same FE space type." 233 " 'basis_from' is a %s and 'basis_to' is a %s", 234 CeedFESpaces[fe_space_from], CeedFESpaces[fe_space_to]); 235 236 // Get source matrices 237 CeedInt dim, q_comp = 1; 238 CeedScalar *interp_to_inv, *interp_from; 239 const CeedScalar *interp_to_source = NULL, *interp_from_source = NULL, *grad_from_source = NULL; 240 241 CeedCall(CeedBasisGetDimension(basis_from, &dim)); 242 if (are_both_tensor) { 243 CeedCall(CeedBasisGetInterp1D(basis_to, &interp_to_source)); 244 CeedCall(CeedBasisGetInterp1D(basis_from, &interp_from_source)); 245 } else { 246 CeedCall(CeedBasisGetNumQuadratureComponents(basis_from, CEED_EVAL_INTERP, &q_comp)); 247 CeedCall(CeedBasisGetInterp(basis_to, &interp_to_source)); 248 CeedCall(CeedBasisGetInterp(basis_from, &interp_from_source)); 249 } 250 CeedCall(CeedMalloc(Q * P_from * q_comp, &interp_from)); 251 CeedCall(CeedCalloc(P_to * P_from, interp_project)); 252 253 // `grad_project = interp_to^+ * grad_from` is computed for the H^1 space case so the 254 // projection basis will have a gradient operation (allocated even if not H^1 for the 255 // basis construction later on) 256 if (fe_space_to == CEED_FE_SPACE_H1) { 257 if (are_both_tensor) { 258 CeedCall(CeedBasisGetGrad1D(basis_from, &grad_from_source)); 259 } else { 260 CeedCall(CeedBasisGetGrad(basis_from, &grad_from_source)); 261 } 262 } 263 CeedCall(CeedCalloc(P_to * P_from * (are_both_tensor ? 1 : dim), grad_project)); 264 265 // Compute interp_to^+, pseudoinverse of interp_to 266 CeedCall(CeedCalloc(Q * q_comp * P_to, &interp_to_inv)); 267 CeedCall(CeedMatrixPseudoinverse(CeedBasisReturnCeed(basis_to), interp_to_source, Q * q_comp, P_to, interp_to_inv)); 268 // Build matrices 269 CeedInt num_matrices = 1 + (fe_space_to == CEED_FE_SPACE_H1) * (are_both_tensor ? 1 : dim); 270 CeedScalar *input_from[num_matrices], *output_project[num_matrices]; 271 272 input_from[0] = (CeedScalar *)interp_from_source; 273 output_project[0] = *interp_project; 274 for (CeedInt m = 1; m < num_matrices; m++) { 275 input_from[m] = (CeedScalar *)&grad_from_source[(m - 1) * Q * P_from]; 276 output_project[m] = &((*grad_project)[(m - 1) * P_to * P_from]); 277 } 278 for (CeedInt m = 0; m < num_matrices; m++) { 279 // output_project = interp_to^+ * interp_from 280 memcpy(interp_from, input_from[m], Q * P_from * q_comp * sizeof(input_from[m][0])); 281 CeedCall(CeedMatrixMatrixMultiply(CeedBasisReturnCeed(basis_to), interp_to_inv, input_from[m], output_project[m], P_to, P_from, Q * q_comp)); 282 // Round zero to machine precision 283 for (CeedInt i = 0; i < P_to * P_from; i++) { 284 if (fabs(output_project[m][i]) < 10 * CEED_EPSILON) output_project[m][i] = 0.0; 285 } 286 } 287 288 // Cleanup 289 CeedCall(CeedFree(&interp_to_inv)); 290 CeedCall(CeedFree(&interp_from)); 291 return CEED_ERROR_SUCCESS; 292 } 293 294 /** 295 @brief Check input vector dimensions for CeedBasisApply[Add]AtPoints 296 297 @param[in] basis `CeedBasis` to evaluate 298 @param[in] num_elem The number of elements to apply the basis evaluation to; 299 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 300 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 301 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 302 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes 303 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 304 @ref CEED_EVAL_GRAD to use gradients, 305 @ref CEED_EVAL_WEIGHT to use quadrature weights 306 @param[in] x_ref `CeedVector` holding reference coordinates of each point 307 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 308 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 309 310 @return An error code: 0 - success, otherwise - failure 311 312 @ref Developer 313 **/ 314 static int CeedBasisApplyAtPointsCheckDims(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 315 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 316 CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1, total_num_points = 0; 317 CeedSize x_length = 0, u_length = 0, v_length; 318 319 CeedCall(CeedBasisGetDimension(basis, &dim)); 320 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 321 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 322 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 323 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &num_q_comp)); 324 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 325 CeedCall(CeedVectorGetLength(v, &v_length)); 326 if (x_ref != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(x_ref, &x_length)); 327 if (u != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(u, &u_length)); 328 329 // Check compatibility coordinates vector 330 for (CeedInt i = 0; i < num_elem; i++) total_num_points += num_points[i]; 331 CeedCheck((x_length >= (CeedSize)total_num_points * (CeedSize)dim) || (eval_mode == CEED_EVAL_WEIGHT), CeedBasisReturnCeed(basis), 332 CEED_ERROR_DIMENSION, 333 "Length of reference coordinate vector incompatible with basis dimension and number of points." 334 " Found reference coordinate vector of length %" CeedSize_FMT ", not of length %" CeedSize_FMT ".", 335 x_length, (CeedSize)total_num_points * (CeedSize)dim); 336 337 // Check CEED_EVAL_WEIGHT only on CEED_NOTRANSPOSE 338 CeedCheck(eval_mode != CEED_EVAL_WEIGHT || t_mode == CEED_NOTRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 339 "CEED_EVAL_WEIGHT only supported with CEED_NOTRANSPOSE"); 340 341 // Check vector lengths to prevent out of bounds issues 342 bool has_good_dims = true; 343 switch (eval_mode) { 344 case CEED_EVAL_INTERP: 345 has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp || 346 v_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)) || 347 (t_mode == CEED_NOTRANSPOSE && (v_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp || 348 u_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp))); 349 break; 350 case CEED_EVAL_GRAD: 351 has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp * (CeedSize)dim || 352 v_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)) || 353 (t_mode == CEED_NOTRANSPOSE && (v_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp * (CeedSize)dim || 354 u_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp))); 355 break; 356 case CEED_EVAL_WEIGHT: 357 has_good_dims = t_mode == CEED_NOTRANSPOSE && (v_length >= total_num_points); 358 break; 359 // LCOV_EXCL_START 360 case CEED_EVAL_NONE: 361 case CEED_EVAL_DIV: 362 case CEED_EVAL_CURL: 363 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points not supported for %s", 364 CeedEvalModes[eval_mode]); 365 // LCOV_EXCL_STOP 366 } 367 CeedCheck(has_good_dims, CeedBasisReturnCeed(basis), CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); 368 return CEED_ERROR_SUCCESS; 369 } 370 371 /** 372 @brief Default implimentation to apply basis evaluation from nodes to arbitrary points 373 374 @param[in] basis `CeedBasis` to evaluate 375 @param[in] apply_add Sum result into target vector or overwrite 376 @param[in] num_elem The number of elements to apply the basis evaluation to; 377 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 378 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 379 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 380 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes 381 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 382 @ref CEED_EVAL_GRAD to use gradients, 383 @ref CEED_EVAL_WEIGHT to use quadrature weights 384 @param[in] x_ref `CeedVector` holding reference coordinates of each point 385 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 386 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 387 388 @return An error code: 0 - success, otherwise - failure 389 390 @ref Developer 391 **/ 392 static int CeedBasisApplyAtPoints_Core(CeedBasis basis, bool apply_add, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 393 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 394 CeedInt dim, num_comp, P_1d = 1, Q_1d = 1, total_num_points = num_points[0]; 395 396 CeedCall(CeedBasisGetDimension(basis, &dim)); 397 // Inserting check because clang-tidy doesn't understand this cannot occur 398 CeedCheck(dim > 0, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Malformed CeedBasis, dim > 0 is required"); 399 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 400 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 401 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 402 403 // Default implementation 404 { 405 bool is_tensor_basis; 406 407 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 408 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 409 "Evaluation at arbitrary points only supported for tensor product bases"); 410 } 411 CeedCheck(num_elem == 1, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 412 "Evaluation at arbitrary points only supported for a single element at a time"); 413 if (eval_mode == CEED_EVAL_WEIGHT) { 414 CeedCall(CeedVectorSetValue(v, 1.0)); 415 return CEED_ERROR_SUCCESS; 416 } 417 if (!basis->basis_chebyshev) { 418 // Build basis mapping from nodes to Chebyshev coefficients 419 CeedScalar *chebyshev_interp_1d, *chebyshev_grad_1d, *chebyshev_q_weight_1d; 420 const CeedScalar *q_ref_1d; 421 Ceed ceed; 422 423 CeedCall(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 424 CeedCall(CeedCalloc(P_1d * Q_1d, &chebyshev_grad_1d)); 425 CeedCall(CeedCalloc(Q_1d, &chebyshev_q_weight_1d)); 426 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); 427 CeedCall(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 428 429 CeedCall(CeedBasisGetCeed(basis, &ceed)); 430 CeedCall(CeedVectorCreate(ceed, num_comp * CeedIntPow(Q_1d, dim), &basis->vec_chebyshev)); 431 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d, Q_1d, chebyshev_interp_1d, chebyshev_grad_1d, q_ref_1d, chebyshev_q_weight_1d, 432 &basis->basis_chebyshev)); 433 434 // Cleanup 435 CeedCall(CeedFree(&chebyshev_interp_1d)); 436 CeedCall(CeedFree(&chebyshev_grad_1d)); 437 CeedCall(CeedFree(&chebyshev_q_weight_1d)); 438 CeedCall(CeedDestroy(&ceed)); 439 } 440 441 // Create TensorContract object if needed, such as a basis from the GPU backends 442 if (!basis->contract) { 443 Ceed ceed_ref; 444 CeedBasis basis_ref = NULL; 445 446 CeedCall(CeedInit("/cpu/self", &ceed_ref)); 447 // Only need matching tensor contraction dimensions, any type of basis will work 448 CeedCall(CeedBasisCreateTensorH1Lagrange(ceed_ref, dim, num_comp, P_1d, Q_1d, CEED_GAUSS, &basis_ref)); 449 // Note - clang-tidy doesn't know basis_ref->contract must be valid here 450 CeedCheck(basis_ref && basis_ref->contract, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 451 "Reference CPU ceed failed to create a tensor contraction object"); 452 CeedCall(CeedTensorContractReferenceCopy(basis_ref->contract, &basis->contract)); 453 CeedCall(CeedBasisDestroy(&basis_ref)); 454 CeedCall(CeedDestroy(&ceed_ref)); 455 } 456 457 // Basis evaluation 458 switch (t_mode) { 459 case CEED_NOTRANSPOSE: { 460 // Nodes to arbitrary points 461 CeedScalar *v_array; 462 const CeedScalar *chebyshev_coeffs, *x_array_read; 463 464 // -- Interpolate to Chebyshev coefficients 465 CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, basis->vec_chebyshev)); 466 467 // -- Evaluate Chebyshev polynomials at arbitrary points 468 CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); 469 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read)); 470 CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array)); 471 switch (eval_mode) { 472 case CEED_EVAL_INTERP: { 473 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 474 475 // ---- Values at point 476 for (CeedInt p = 0; p < total_num_points; p++) { 477 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; 478 479 for (CeedInt d = 0; d < dim; d++) { 480 // ------ Tensor contract with current Chebyshev polynomial values 481 CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 482 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, 483 d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2])); 484 pre /= Q_1d; 485 post *= 1; 486 } 487 for (CeedInt c = 0; c < num_comp; c++) v_array[c * total_num_points + p] = tmp[dim % 2][c]; 488 } 489 break; 490 } 491 case CEED_EVAL_GRAD: { 492 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 493 494 // ---- Values at point 495 for (CeedInt p = 0; p < total_num_points; p++) { 496 // Dim**2 contractions, apply grad when pass == dim 497 for (CeedInt pass = 0; pass < dim; pass++) { 498 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; 499 500 for (CeedInt d = 0; d < dim; d++) { 501 // ------ Tensor contract with current Chebyshev polynomial values 502 if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 503 else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 504 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, 505 d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2])); 506 pre /= Q_1d; 507 post *= 1; 508 } 509 for (CeedInt c = 0; c < num_comp; c++) v_array[(pass * num_comp + c) * total_num_points + p] = tmp[dim % 2][c]; 510 } 511 } 512 break; 513 } 514 default: 515 // Nothing to do, excluded above 516 break; 517 } 518 CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs)); 519 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read)); 520 CeedCall(CeedVectorRestoreArray(v, &v_array)); 521 break; 522 } 523 case CEED_TRANSPOSE: { 524 // Note: No switch on e_mode here because only CEED_EVAL_INTERP is supported at this time 525 // Arbitrary points to nodes 526 CeedScalar *chebyshev_coeffs; 527 const CeedScalar *u_array, *x_array_read; 528 529 // -- Transpose of evaluation of Chebyshev polynomials at arbitrary points 530 CeedCall(CeedVectorGetArrayWrite(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); 531 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read)); 532 CeedCall(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array)); 533 534 switch (eval_mode) { 535 case CEED_EVAL_INTERP: { 536 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 537 538 // ---- Values at point 539 for (CeedInt p = 0; p < total_num_points; p++) { 540 CeedInt pre = num_comp * 1, post = 1; 541 542 for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[c * total_num_points + p]; 543 for (CeedInt d = 0; d < dim; d++) { 544 // ------ Tensor contract with current Chebyshev polynomial values 545 CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 546 CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == (dim - 1), tmp[d % 2], 547 d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); 548 pre /= 1; 549 post *= Q_1d; 550 } 551 } 552 break; 553 } 554 case CEED_EVAL_GRAD: { 555 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 556 557 // ---- Values at point 558 for (CeedInt p = 0; p < total_num_points; p++) { 559 // Dim**2 contractions, apply grad when pass == dim 560 for (CeedInt pass = 0; pass < dim; pass++) { 561 CeedInt pre = num_comp * 1, post = 1; 562 563 for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[(pass * num_comp + c) * total_num_points + p]; 564 for (CeedInt d = 0; d < dim; d++) { 565 // ------ Tensor contract with current Chebyshev polynomial values 566 if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 567 else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 568 CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, 569 (p > 0 || (p == 0 && pass > 0)) && d == (dim - 1), tmp[d % 2], 570 d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); 571 pre /= 1; 572 post *= Q_1d; 573 } 574 } 575 } 576 break; 577 } 578 default: 579 // Nothing to do, excluded above 580 break; 581 } 582 CeedCall(CeedVectorRestoreArray(basis->vec_chebyshev, &chebyshev_coeffs)); 583 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read)); 584 CeedCall(CeedVectorRestoreArrayRead(u, &u_array)); 585 586 // -- Interpolate transpose from Chebyshev coefficients 587 if (apply_add) CeedCall(CeedBasisApplyAdd(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); 588 else CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); 589 break; 590 } 591 } 592 return CEED_ERROR_SUCCESS; 593 } 594 595 /// @} 596 597 /// ---------------------------------------------------------------------------- 598 /// Ceed Backend API 599 /// ---------------------------------------------------------------------------- 600 /// @addtogroup CeedBasisBackend 601 /// @{ 602 603 /** 604 @brief Return collocated gradient matrix 605 606 @param[in] basis `CeedBasis` 607 @param[out] collo_grad_1d Row-major (`Q_1d * Q_1d`) matrix expressing derivatives of basis functions at quadrature points 608 609 @return An error code: 0 - success, otherwise - failure 610 611 @ref Backend 612 **/ 613 int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { 614 Ceed ceed; 615 CeedInt P_1d, Q_1d; 616 CeedScalar *interp_1d_pinv; 617 const CeedScalar *grad_1d, *interp_1d; 618 619 // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure. 620 CeedCall(CeedBasisGetCeed(basis, &ceed)); 621 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 622 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 623 624 // Compute interp_1d^+, pseudoinverse of interp_1d 625 CeedCall(CeedCalloc(P_1d * Q_1d, &interp_1d_pinv)); 626 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 627 CeedCall(CeedMatrixPseudoinverse(ceed, interp_1d, Q_1d, P_1d, interp_1d_pinv)); 628 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 629 CeedCall(CeedMatrixMatrixMultiply(ceed, grad_1d, (const CeedScalar *)interp_1d_pinv, collo_grad_1d, Q_1d, Q_1d, P_1d)); 630 631 CeedCall(CeedFree(&interp_1d_pinv)); 632 CeedCall(CeedDestroy(&ceed)); 633 return CEED_ERROR_SUCCESS; 634 } 635 636 /** 637 @brief Return 1D interpolation matrix to Chebyshev polynomial coefficients on quadrature space 638 639 @param[in] basis `CeedBasis` 640 @param[out] chebyshev_interp_1d Row-major (`P_1d * Q_1d`) matrix interpolating from basis nodes to Chebyshev polynomial coefficients 641 642 @return An error code: 0 - success, otherwise - failure 643 644 @ref Backend 645 **/ 646 int CeedBasisGetChebyshevInterp1D(CeedBasis basis, CeedScalar *chebyshev_interp_1d) { 647 CeedInt P_1d, Q_1d; 648 CeedScalar *C, *chebyshev_coeffs_1d_inv; 649 const CeedScalar *interp_1d, *q_ref_1d; 650 Ceed ceed; 651 652 CeedCall(CeedBasisGetCeed(basis, &ceed)); 653 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 654 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 655 656 // Build coefficient matrix 657 // -- Note: Clang-tidy needs this check 658 CeedCheck((P_1d > 0) && (Q_1d > 0), ceed, CEED_ERROR_INCOMPATIBLE, "CeedBasis dimensions are malformed"); 659 CeedCall(CeedCalloc(Q_1d * Q_1d, &C)); 660 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); 661 for (CeedInt i = 0; i < Q_1d; i++) CeedCall(CeedChebyshevPolynomialsAtPoint(q_ref_1d[i], Q_1d, &C[i * Q_1d])); 662 663 // Compute C^+, pseudoinverse of coefficient matrix 664 CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_coeffs_1d_inv)); 665 CeedCall(CeedMatrixPseudoinverse(ceed, C, Q_1d, Q_1d, chebyshev_coeffs_1d_inv)); 666 667 // Build mapping from nodes to Chebyshev coefficients 668 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 669 CeedCall(CeedMatrixMatrixMultiply(ceed, chebyshev_coeffs_1d_inv, interp_1d, chebyshev_interp_1d, Q_1d, P_1d, Q_1d)); 670 671 // Cleanup 672 CeedCall(CeedFree(&C)); 673 CeedCall(CeedFree(&chebyshev_coeffs_1d_inv)); 674 CeedCall(CeedDestroy(&ceed)); 675 return CEED_ERROR_SUCCESS; 676 } 677 678 /** 679 @brief Get tensor status for given `CeedBasis` 680 681 @param[in] basis `CeedBasis` 682 @param[out] is_tensor Variable to store tensor status 683 684 @return An error code: 0 - success, otherwise - failure 685 686 @ref Backend 687 **/ 688 int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor) { 689 *is_tensor = basis->is_tensor_basis; 690 return CEED_ERROR_SUCCESS; 691 } 692 693 /** 694 @brief Get backend data of a `CeedBasis` 695 696 @param[in] basis `CeedBasis` 697 @param[out] data Variable to store data 698 699 @return An error code: 0 - success, otherwise - failure 700 701 @ref Backend 702 **/ 703 int CeedBasisGetData(CeedBasis basis, void *data) { 704 *(void **)data = basis->data; 705 return CEED_ERROR_SUCCESS; 706 } 707 708 /** 709 @brief Set backend data of a `CeedBasis` 710 711 @param[in,out] basis `CeedBasis` 712 @param[in] data Data to set 713 714 @return An error code: 0 - success, otherwise - failure 715 716 @ref Backend 717 **/ 718 int CeedBasisSetData(CeedBasis basis, void *data) { 719 basis->data = data; 720 return CEED_ERROR_SUCCESS; 721 } 722 723 /** 724 @brief Increment the reference counter for a `CeedBasis` 725 726 @param[in,out] basis `CeedBasis` to increment the reference counter 727 728 @return An error code: 0 - success, otherwise - failure 729 730 @ref Backend 731 **/ 732 int CeedBasisReference(CeedBasis basis) { 733 basis->ref_count++; 734 return CEED_ERROR_SUCCESS; 735 } 736 737 /** 738 @brief Get number of Q-vector components for given `CeedBasis` 739 740 @param[in] basis `CeedBasis` 741 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 742 @ref CEED_EVAL_GRAD to use gradients, 743 @ref CEED_EVAL_DIV to use divergence, 744 @ref CEED_EVAL_CURL to use curl 745 @param[out] q_comp Variable to store number of Q-vector components of basis 746 747 @return An error code: 0 - success, otherwise - failure 748 749 @ref Backend 750 **/ 751 int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedEvalMode eval_mode, CeedInt *q_comp) { 752 CeedInt dim; 753 754 CeedCall(CeedBasisGetDimension(basis, &dim)); 755 switch (eval_mode) { 756 case CEED_EVAL_INTERP: { 757 CeedFESpace fe_space; 758 759 CeedCall(CeedBasisGetFESpace(basis, &fe_space)); 760 *q_comp = (fe_space == CEED_FE_SPACE_H1) ? 1 : dim; 761 } break; 762 case CEED_EVAL_GRAD: 763 *q_comp = dim; 764 break; 765 case CEED_EVAL_DIV: 766 *q_comp = 1; 767 break; 768 case CEED_EVAL_CURL: 769 *q_comp = (dim < 3) ? 1 : dim; 770 break; 771 case CEED_EVAL_NONE: 772 case CEED_EVAL_WEIGHT: 773 *q_comp = 1; 774 break; 775 } 776 return CEED_ERROR_SUCCESS; 777 } 778 779 /** 780 @brief Estimate number of FLOPs required to apply `CeedBasis` in `t_mode` and `eval_mode` 781 782 @param[in] basis `CeedBasis` to estimate FLOPs for 783 @param[in] t_mode Apply basis or transpose 784 @param[in] eval_mode @ref CeedEvalMode 785 @param[in] is_at_points Evaluate the basis at points or quadrature points 786 @param[in] num_points Number of points basis is evaluated at 787 @param[out] flops Address of variable to hold FLOPs estimate 788 789 @ref Backend 790 **/ 791 int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEvalMode eval_mode, bool is_at_points, CeedInt num_points, 792 CeedSize *flops) { 793 bool is_tensor; 794 795 CeedCall(CeedBasisIsTensor(basis, &is_tensor)); 796 CeedCheck(!is_at_points || is_tensor, CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Can only evaluate tensor-product bases at points"); 797 if (is_tensor) { 798 CeedInt dim, num_comp, P_1d, Q_1d; 799 800 CeedCall(CeedBasisGetDimension(basis, &dim)); 801 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 802 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 803 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 804 if (t_mode == CEED_TRANSPOSE) { 805 P_1d = Q_1d; 806 Q_1d = P_1d; 807 } 808 CeedInt tensor_flops = 0, pre = num_comp * CeedIntPow(P_1d, dim - 1), post = 1; 809 810 for (CeedInt d = 0; d < dim; d++) { 811 tensor_flops += 2 * pre * P_1d * post * Q_1d; 812 pre /= P_1d; 813 post *= Q_1d; 814 } 815 if (is_at_points) { 816 CeedInt chebyshev_flops = (Q_1d - 2) * 3 + 1, d_chebyshev_flops = (Q_1d - 2) * 8 + 1; 817 CeedInt point_tensor_flops = 0, pre = CeedIntPow(Q_1d, dim - 1), post = 1; 818 819 for (CeedInt d = 0; d < dim; d++) { 820 point_tensor_flops += 2 * pre * Q_1d * post * 1; 821 pre /= P_1d; 822 post *= Q_1d; 823 } 824 825 switch (eval_mode) { 826 case CEED_EVAL_NONE: 827 *flops = 0; 828 break; 829 case CEED_EVAL_INTERP: 830 *flops = tensor_flops + num_points * (dim * chebyshev_flops + point_tensor_flops + (t_mode == CEED_TRANSPOSE ? CeedIntPow(Q_1d, dim) : 0)); 831 break; 832 case CEED_EVAL_GRAD: 833 *flops = tensor_flops + num_points * (dim * (d_chebyshev_flops + (dim - 1) * chebyshev_flops + point_tensor_flops + 834 (t_mode == CEED_TRANSPOSE ? CeedIntPow(Q_1d, dim) : 0))); 835 break; 836 case CEED_EVAL_DIV: 837 case CEED_EVAL_CURL: { 838 // LCOV_EXCL_START 839 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Tensor basis evaluation for %s not supported", 840 CeedEvalModes[eval_mode]); 841 break; 842 // LCOV_EXCL_STOP 843 } 844 case CEED_EVAL_WEIGHT: 845 *flops = num_points; 846 break; 847 } 848 } else { 849 switch (eval_mode) { 850 case CEED_EVAL_NONE: 851 *flops = 0; 852 break; 853 case CEED_EVAL_INTERP: 854 *flops = tensor_flops; 855 break; 856 case CEED_EVAL_GRAD: 857 *flops = tensor_flops * 2; 858 break; 859 case CEED_EVAL_DIV: 860 case CEED_EVAL_CURL: { 861 // LCOV_EXCL_START 862 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Tensor basis evaluation for %s not supported", 863 CeedEvalModes[eval_mode]); 864 break; 865 // LCOV_EXCL_STOP 866 } 867 case CEED_EVAL_WEIGHT: 868 *flops = dim * CeedIntPow(Q_1d, dim); 869 break; 870 } 871 } 872 } else { 873 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; 874 875 CeedCall(CeedBasisGetDimension(basis, &dim)); 876 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 877 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 878 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 879 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 880 switch (eval_mode) { 881 case CEED_EVAL_NONE: 882 *flops = 0; 883 break; 884 case CEED_EVAL_INTERP: 885 case CEED_EVAL_GRAD: 886 case CEED_EVAL_DIV: 887 case CEED_EVAL_CURL: 888 *flops = num_nodes * num_qpts * num_comp * q_comp; 889 break; 890 case CEED_EVAL_WEIGHT: 891 *flops = 0; 892 break; 893 } 894 } 895 return CEED_ERROR_SUCCESS; 896 } 897 898 /** 899 @brief Get `CeedFESpace` for a `CeedBasis` 900 901 @param[in] basis `CeedBasis` 902 @param[out] fe_space Variable to store `CeedFESpace` 903 904 @return An error code: 0 - success, otherwise - failure 905 906 @ref Backend 907 **/ 908 int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space) { 909 *fe_space = basis->fe_space; 910 return CEED_ERROR_SUCCESS; 911 } 912 913 /** 914 @brief Get dimension for given `CeedElemTopology` 915 916 @param[in] topo `CeedElemTopology` 917 @param[out] dim Variable to store dimension of topology 918 919 @return An error code: 0 - success, otherwise - failure 920 921 @ref Backend 922 **/ 923 int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) { 924 *dim = (CeedInt)topo >> 16; 925 return CEED_ERROR_SUCCESS; 926 } 927 928 /** 929 @brief Get `CeedTensorContract` of a `CeedBasis` 930 931 @param[in] basis `CeedBasis` 932 @param[out] contract Variable to store `CeedTensorContract` 933 934 @return An error code: 0 - success, otherwise - failure 935 936 @ref Backend 937 **/ 938 int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) { 939 *contract = basis->contract; 940 return CEED_ERROR_SUCCESS; 941 } 942 943 /** 944 @brief Set `CeedTensorContract` of a `CeedBasis` 945 946 @param[in,out] basis `CeedBasis` 947 @param[in] contract `CeedTensorContract` to set 948 949 @return An error code: 0 - success, otherwise - failure 950 951 @ref Backend 952 **/ 953 int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract) { 954 basis->contract = contract; 955 CeedCall(CeedTensorContractReference(contract)); 956 return CEED_ERROR_SUCCESS; 957 } 958 959 /** 960 @brief Return a reference implementation of matrix multiplication \f$C = A B\f$. 961 962 Note: This is a reference implementation for CPU `CeedScalar` pointers that is not intended for high performance. 963 964 @param[in] ceed `Ceed` context for error handling 965 @param[in] mat_A Row-major matrix `A` 966 @param[in] mat_B Row-major matrix `B` 967 @param[out] mat_C Row-major output matrix `C` 968 @param[in] m Number of rows of `C` 969 @param[in] n Number of columns of `C` 970 @param[in] kk Number of columns of `A`/rows of `B` 971 972 @return An error code: 0 - success, otherwise - failure 973 974 @ref Utility 975 **/ 976 int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n, CeedInt kk) { 977 for (CeedInt i = 0; i < m; i++) { 978 for (CeedInt j = 0; j < n; j++) { 979 CeedScalar sum = 0; 980 981 for (CeedInt k = 0; k < kk; k++) sum += mat_A[k + i * kk] * mat_B[j + k * n]; 982 mat_C[j + i * n] = sum; 983 } 984 } 985 return CEED_ERROR_SUCCESS; 986 } 987 988 /** 989 @brief Return QR Factorization of a matrix 990 991 @param[in] ceed `Ceed` context for error handling 992 @param[in,out] mat Row-major matrix to be factorized in place 993 @param[in,out] tau Vector of length `m` of scaling factors 994 @param[in] m Number of rows 995 @param[in] n Number of columns 996 997 @return An error code: 0 - success, otherwise - failure 998 999 @ref Utility 1000 **/ 1001 int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n) { 1002 CeedScalar v[m]; 1003 1004 // Check matrix shape 1005 CeedCheck(n <= m, ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute QR factorization with n > m"); 1006 1007 for (CeedInt i = 0; i < n; i++) { 1008 CeedScalar sigma = 0.0; 1009 1010 if (i >= m - 1) { // last row of matrix, no reflection needed 1011 tau[i] = 0.; 1012 break; 1013 } 1014 // Calculate Householder vector, magnitude 1015 v[i] = mat[i + n * i]; 1016 for (CeedInt j = i + 1; j < m; j++) { 1017 v[j] = mat[i + n * j]; 1018 sigma += v[j] * v[j]; 1019 } 1020 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:m] 1021 const CeedScalar R_ii = -copysign(norm, v[i]); 1022 1023 v[i] -= R_ii; 1024 // norm of v[i:m] after modification above and scaling below 1025 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1026 // tau = 2 / (norm*norm) 1027 tau[i] = 2 * v[i] * v[i] / (v[i] * v[i] + sigma); 1028 for (CeedInt j = i + 1; j < m; j++) v[j] /= v[i]; 1029 1030 // Apply Householder reflector to lower right panel 1031 CeedHouseholderReflect(&mat[i * n + i + 1], &v[i], tau[i], m - i, n - i - 1, n, 1); 1032 // Save v 1033 mat[i + n * i] = R_ii; 1034 for (CeedInt j = i + 1; j < m; j++) mat[i + n * j] = v[j]; 1035 } 1036 return CEED_ERROR_SUCCESS; 1037 } 1038 1039 /** 1040 @brief Apply Householder Q matrix 1041 1042 Compute `mat_A = mat_Q mat_A`, where `mat_Q` is \f$m \times m\f$ and `mat_A` is \f$m \times n\f$. 1043 1044 @param[in,out] mat_A Matrix to apply Householder Q to, in place 1045 @param[in] mat_Q Householder Q matrix 1046 @param[in] tau Householder scaling factors 1047 @param[in] t_mode Transpose mode for application 1048 @param[in] m Number of rows in `A` 1049 @param[in] n Number of columns in `A` 1050 @param[in] k Number of elementary reflectors in Q, `k < m` 1051 @param[in] row Row stride in `A` 1052 @param[in] col Col stride in `A` 1053 1054 @return An error code: 0 - success, otherwise - failure 1055 1056 @ref Utility 1057 **/ 1058 int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n, 1059 CeedInt k, CeedInt row, CeedInt col) { 1060 CeedScalar *v; 1061 1062 CeedCall(CeedMalloc(m, &v)); 1063 for (CeedInt ii = 0; ii < k; ii++) { 1064 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k - 1 - ii; 1065 for (CeedInt j = i + 1; j < m; j++) v[j] = mat_Q[j * k + i]; 1066 // Apply Householder reflector (I - tau v v^T) collo_grad_1d^T 1067 CeedCall(CeedHouseholderReflect(&mat_A[i * row], &v[i], tau[i], m - i, n, row, col)); 1068 } 1069 CeedCall(CeedFree(&v)); 1070 return CEED_ERROR_SUCCESS; 1071 } 1072 1073 /** 1074 @brief Return pseudoinverse of a matrix 1075 1076 @param[in] ceed Ceed context for error handling 1077 @param[in] mat Row-major matrix to compute pseudoinverse of 1078 @param[in] m Number of rows 1079 @param[in] n Number of columns 1080 @param[out] mat_pinv Row-major pseudoinverse matrix 1081 1082 @return An error code: 0 - success, otherwise - failure 1083 1084 @ref Utility 1085 **/ 1086 int CeedMatrixPseudoinverse(Ceed ceed, const CeedScalar *mat, CeedInt m, CeedInt n, CeedScalar *mat_pinv) { 1087 CeedScalar *tau, *I, *mat_copy; 1088 1089 CeedCall(CeedCalloc(m, &tau)); 1090 CeedCall(CeedCalloc(m * m, &I)); 1091 CeedCall(CeedCalloc(m * n, &mat_copy)); 1092 memcpy(mat_copy, mat, m * n * sizeof mat[0]); 1093 1094 // QR Factorization, mat = Q R 1095 CeedCall(CeedQRFactorization(ceed, mat_copy, tau, m, n)); 1096 1097 // -- Apply Q^T, I = Q^T * I 1098 for (CeedInt i = 0; i < m; i++) I[i * m + i] = 1.0; 1099 CeedCall(CeedHouseholderApplyQ(I, mat_copy, tau, CEED_TRANSPOSE, m, m, n, m, 1)); 1100 // -- Apply R_inv, mat_pinv = R_inv * Q^T 1101 for (CeedInt j = 0; j < m; j++) { // Column j 1102 mat_pinv[j + m * (n - 1)] = I[j + m * (n - 1)] / mat_copy[n * n - 1]; 1103 for (CeedInt i = n - 2; i >= 0; i--) { // Row i 1104 mat_pinv[j + m * i] = I[j + m * i]; 1105 for (CeedInt k = i + 1; k < n; k++) mat_pinv[j + m * i] -= mat_copy[k + n * i] * mat_pinv[j + m * k]; 1106 mat_pinv[j + m * i] /= mat_copy[i + n * i]; 1107 } 1108 } 1109 1110 // Cleanup 1111 CeedCall(CeedFree(&I)); 1112 CeedCall(CeedFree(&tau)); 1113 CeedCall(CeedFree(&mat_copy)); 1114 return CEED_ERROR_SUCCESS; 1115 } 1116 1117 /** 1118 @brief Return symmetric Schur decomposition of the symmetric matrix mat via symmetric QR factorization 1119 1120 @param[in] ceed `Ceed` context for error handling 1121 @param[in,out] mat Row-major matrix to be factorized in place 1122 @param[out] lambda Vector of length n of eigenvalues 1123 @param[in] n Number of rows/columns 1124 1125 @return An error code: 0 - success, otherwise - failure 1126 1127 @ref Utility 1128 **/ 1129 CeedPragmaOptimizeOff 1130 int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n) { 1131 // Check bounds for clang-tidy 1132 CeedCheck(n > 1, ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute symmetric Schur decomposition of scalars"); 1133 1134 CeedScalar v[n - 1], tau[n - 1], mat_T[n * n]; 1135 1136 // Copy mat to mat_T and set mat to I 1137 memcpy(mat_T, mat, n * n * sizeof(mat[0])); 1138 for (CeedInt i = 0; i < n; i++) { 1139 for (CeedInt j = 0; j < n; j++) mat[j + n * i] = (i == j) ? 1 : 0; 1140 } 1141 1142 // Reduce to tridiagonal 1143 for (CeedInt i = 0; i < n - 1; i++) { 1144 // Calculate Householder vector, magnitude 1145 CeedScalar sigma = 0.0; 1146 1147 v[i] = mat_T[i + n * (i + 1)]; 1148 for (CeedInt j = i + 1; j < n - 1; j++) { 1149 v[j] = mat_T[i + n * (j + 1)]; 1150 sigma += v[j] * v[j]; 1151 } 1152 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:n-1] 1153 const CeedScalar R_ii = -copysign(norm, v[i]); 1154 1155 v[i] -= R_ii; 1156 // norm of v[i:m] after modification above and scaling below 1157 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1158 // tau = 2 / (norm*norm) 1159 tau[i] = i == n - 2 ? 2 : 2 * v[i] * v[i] / (v[i] * v[i] + sigma); 1160 for (CeedInt j = i + 1; j < n - 1; j++) v[j] /= v[i]; 1161 1162 // Update sub and super diagonal 1163 for (CeedInt j = i + 2; j < n; j++) { 1164 mat_T[i + n * j] = 0; 1165 mat_T[j + n * i] = 0; 1166 } 1167 // Apply symmetric Householder reflector to lower right panel 1168 CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1); 1169 CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), 1, n); 1170 1171 // Save v 1172 mat_T[i + n * (i + 1)] = R_ii; 1173 mat_T[(i + 1) + n * i] = R_ii; 1174 for (CeedInt j = i + 1; j < n - 1; j++) { 1175 mat_T[i + n * (j + 1)] = v[j]; 1176 } 1177 } 1178 // Backwards accumulation of Q 1179 for (CeedInt i = n - 2; i >= 0; i--) { 1180 if (tau[i] > 0.0) { 1181 v[i] = 1; 1182 for (CeedInt j = i + 1; j < n - 1; j++) { 1183 v[j] = mat_T[i + n * (j + 1)]; 1184 mat_T[i + n * (j + 1)] = 0; 1185 } 1186 CeedHouseholderReflect(&mat[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1); 1187 } 1188 } 1189 1190 // Reduce sub and super diagonal 1191 CeedInt p = 0, q = 0, itr = 0, max_itr = n * n * n * n; 1192 CeedScalar tol = CEED_EPSILON; 1193 1194 while (itr < max_itr) { 1195 // Update p, q, size of reduced portions of diagonal 1196 p = 0; 1197 q = 0; 1198 for (CeedInt i = n - 2; i >= 0; i--) { 1199 if (fabs(mat_T[i + n * (i + 1)]) < tol) q += 1; 1200 else break; 1201 } 1202 for (CeedInt i = 0; i < n - q - 1; i++) { 1203 if (fabs(mat_T[i + n * (i + 1)]) < tol) p += 1; 1204 else break; 1205 } 1206 if (q == n - 1) break; // Finished reducing 1207 1208 // Reduce tridiagonal portion 1209 CeedScalar t_nn = mat_T[(n - 1 - q) + n * (n - 1 - q)], t_nnm1 = mat_T[(n - 2 - q) + n * (n - 1 - q)]; 1210 CeedScalar d = (mat_T[(n - 2 - q) + n * (n - 2 - q)] - t_nn) / 2; 1211 CeedScalar mu = t_nn - t_nnm1 * t_nnm1 / (d + copysign(sqrt(d * d + t_nnm1 * t_nnm1), d)); 1212 CeedScalar x = mat_T[p + n * p] - mu; 1213 CeedScalar z = mat_T[p + n * (p + 1)]; 1214 1215 for (CeedInt k = p; k < n - q - 1; k++) { 1216 // Compute Givens rotation 1217 CeedScalar c = 1, s = 0; 1218 1219 if (fabs(z) > tol) { 1220 if (fabs(z) > fabs(x)) { 1221 const CeedScalar tau = -x / z; 1222 1223 s = 1 / sqrt(1 + tau * tau); 1224 c = s * tau; 1225 } else { 1226 const CeedScalar tau = -z / x; 1227 1228 c = 1 / sqrt(1 + tau * tau); 1229 s = c * tau; 1230 } 1231 } 1232 1233 // Apply Givens rotation to T 1234 CeedGivensRotation(mat_T, c, s, CEED_NOTRANSPOSE, k, k + 1, n, n); 1235 CeedGivensRotation(mat_T, c, s, CEED_TRANSPOSE, k, k + 1, n, n); 1236 1237 // Apply Givens rotation to Q 1238 CeedGivensRotation(mat, c, s, CEED_NOTRANSPOSE, k, k + 1, n, n); 1239 1240 // Update x, z 1241 if (k < n - q - 2) { 1242 x = mat_T[k + n * (k + 1)]; 1243 z = mat_T[k + n * (k + 2)]; 1244 } 1245 } 1246 itr++; 1247 } 1248 1249 // Save eigenvalues 1250 for (CeedInt i = 0; i < n; i++) lambda[i] = mat_T[i + n * i]; 1251 1252 // Check convergence 1253 CeedCheck(itr < max_itr || q > n, ceed, CEED_ERROR_MINOR, "Symmetric QR failed to converge"); 1254 return CEED_ERROR_SUCCESS; 1255 } 1256 CeedPragmaOptimizeOn 1257 1258 /** 1259 @brief Return Simultaneous Diagonalization of two matrices. 1260 1261 This solves the generalized eigenvalue problem `A x = lambda B x`, where `A` and `B` are symmetric and `B` is positive definite. 1262 We generate the matrix `X` and vector `Lambda` such that `X^T A X = Lambda` and `X^T B X = I`. 1263 This is equivalent to the LAPACK routine 'sygv' with `TYPE = 1`. 1264 1265 @param[in] ceed `Ceed` context for error handling 1266 @param[in] mat_A Row-major matrix to be factorized with eigenvalues 1267 @param[in] mat_B Row-major matrix to be factorized to identity 1268 @param[out] mat_X Row-major orthogonal matrix 1269 @param[out] lambda Vector of length `n` of generalized eigenvalues 1270 @param[in] n Number of rows/columns 1271 1272 @return An error code: 0 - success, otherwise - failure 1273 1274 @ref Utility 1275 **/ 1276 CeedPragmaOptimizeOff 1277 int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *mat_X, CeedScalar *lambda, CeedInt n) { 1278 CeedScalar *mat_C, *mat_G, *vec_D; 1279 1280 CeedCall(CeedCalloc(n * n, &mat_C)); 1281 CeedCall(CeedCalloc(n * n, &mat_G)); 1282 CeedCall(CeedCalloc(n, &vec_D)); 1283 1284 // Compute B = G D G^T 1285 memcpy(mat_G, mat_B, n * n * sizeof(mat_B[0])); 1286 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n)); 1287 1288 // Sort eigenvalues 1289 for (CeedInt i = n - 1; i >= 0; i--) { 1290 for (CeedInt j = 0; j < i; j++) { 1291 if (fabs(vec_D[j]) > fabs(vec_D[j + 1])) { 1292 CeedScalarSwap(vec_D[j], vec_D[j + 1]); 1293 for (CeedInt k = 0; k < n; k++) CeedScalarSwap(mat_G[k * n + j], mat_G[k * n + j + 1]); 1294 } 1295 } 1296 } 1297 1298 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T 1299 // = D^-1/2 G^T A G D^-1/2 1300 // -- D = D^-1/2 1301 for (CeedInt i = 0; i < n; i++) vec_D[i] = 1. / sqrt(vec_D[i]); 1302 // -- G = G D^-1/2 1303 // -- C = D^-1/2 G^T 1304 for (CeedInt i = 0; i < n; i++) { 1305 for (CeedInt j = 0; j < n; j++) { 1306 mat_G[i * n + j] *= vec_D[j]; 1307 mat_C[j * n + i] = mat_G[i * n + j]; 1308 } 1309 } 1310 // -- X = (D^-1/2 G^T) A 1311 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_C, (const CeedScalar *)mat_A, mat_X, n, n, n)); 1312 // -- C = (D^-1/2 G^T A) (G D^-1/2) 1313 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_X, (const CeedScalar *)mat_G, mat_C, n, n, n)); 1314 1315 // Compute Q^T C Q = lambda 1316 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_C, lambda, n)); 1317 1318 // Sort eigenvalues 1319 for (CeedInt i = n - 1; i >= 0; i--) { 1320 for (CeedInt j = 0; j < i; j++) { 1321 if (fabs(lambda[j]) > fabs(lambda[j + 1])) { 1322 CeedScalarSwap(lambda[j], lambda[j + 1]); 1323 for (CeedInt k = 0; k < n; k++) CeedScalarSwap(mat_C[k * n + j], mat_C[k * n + j + 1]); 1324 } 1325 } 1326 } 1327 1328 // Set X = (G D^1/2)^-T Q 1329 // = G D^-1/2 Q 1330 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_G, (const CeedScalar *)mat_C, mat_X, n, n, n)); 1331 1332 // Cleanup 1333 CeedCall(CeedFree(&mat_C)); 1334 CeedCall(CeedFree(&mat_G)); 1335 CeedCall(CeedFree(&vec_D)); 1336 return CEED_ERROR_SUCCESS; 1337 } 1338 CeedPragmaOptimizeOn 1339 1340 /// @} 1341 1342 /// ---------------------------------------------------------------------------- 1343 /// CeedBasis Public API 1344 /// ---------------------------------------------------------------------------- 1345 /// @addtogroup CeedBasisUser 1346 /// @{ 1347 1348 /** 1349 @brief Create a tensor-product basis for \f$H^1\f$ discretizations 1350 1351 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1352 @param[in] dim Topological dimension 1353 @param[in] num_comp Number of field components (1 for scalar fields) 1354 @param[in] P_1d Number of nodes in one dimension 1355 @param[in] Q_1d Number of quadrature points in one dimension 1356 @param[in] interp_1d Row-major (`Q_1d * P_1d`) matrix expressing the values of nodal basis functions at quadrature points 1357 @param[in] grad_1d Row-major (`Q_1d * P_1d`) matrix expressing derivatives of nodal basis functions at quadrature points 1358 @param[in] q_ref_1d Array of length `Q_1d` holding the locations of quadrature points on the 1D reference element `[-1, 1]` 1359 @param[in] q_weight_1d Array of length `Q_1d` holding the quadrature weights on the reference element 1360 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1361 1362 @return An error code: 0 - success, otherwise - failure 1363 1364 @ref User 1365 **/ 1366 int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, 1367 const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis *basis) { 1368 if (!ceed->BasisCreateTensorH1) { 1369 Ceed delegate; 1370 1371 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1372 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateTensorH1"); 1373 CeedCall(CeedBasisCreateTensorH1(delegate, dim, num_comp, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, basis)); 1374 CeedCall(CeedDestroy(&delegate)); 1375 return CEED_ERROR_SUCCESS; 1376 } 1377 1378 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis dimension must be a positive value"); 1379 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1380 CeedCheck(P_1d > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1381 CeedCheck(Q_1d > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1382 1383 CeedElemTopology topo = dim == 1 ? CEED_TOPOLOGY_LINE : dim == 2 ? CEED_TOPOLOGY_QUAD : CEED_TOPOLOGY_HEX; 1384 1385 CeedCall(CeedCalloc(1, basis)); 1386 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1387 (*basis)->ref_count = 1; 1388 (*basis)->is_tensor_basis = true; 1389 (*basis)->dim = dim; 1390 (*basis)->topo = topo; 1391 (*basis)->num_comp = num_comp; 1392 (*basis)->P_1d = P_1d; 1393 (*basis)->Q_1d = Q_1d; 1394 (*basis)->P = CeedIntPow(P_1d, dim); 1395 (*basis)->Q = CeedIntPow(Q_1d, dim); 1396 (*basis)->fe_space = CEED_FE_SPACE_H1; 1397 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_ref_1d)); 1398 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_weight_1d)); 1399 if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d * sizeof(q_ref_1d[0])); 1400 if (q_weight_1d) memcpy((*basis)->q_weight_1d, q_weight_1d, Q_1d * sizeof(q_weight_1d[0])); 1401 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->interp_1d)); 1402 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->grad_1d)); 1403 if (interp_1d) memcpy((*basis)->interp_1d, interp_1d, Q_1d * P_1d * sizeof(interp_1d[0])); 1404 if (grad_1d) memcpy((*basis)->grad_1d, grad_1d, Q_1d * P_1d * sizeof(grad_1d[0])); 1405 CeedCall(ceed->BasisCreateTensorH1(dim, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, *basis)); 1406 return CEED_ERROR_SUCCESS; 1407 } 1408 1409 /** 1410 @brief Create a tensor-product \f$H^1\f$ Lagrange basis 1411 1412 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1413 @param[in] dim Topological dimension of element 1414 @param[in] num_comp Number of field components (1 for scalar fields) 1415 @param[in] P Number of Gauss-Lobatto nodes in one dimension. 1416 The polynomial degree of the resulting `Q_k` element is `k = P - 1`. 1417 @param[in] Q Number of quadrature points in one dimension. 1418 @param[in] quad_mode Distribution of the `Q` quadrature points (affects order of accuracy for the quadrature) 1419 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1420 1421 @return An error code: 0 - success, otherwise - failure 1422 1423 @ref User 1424 **/ 1425 int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, CeedBasis *basis) { 1426 // Allocate 1427 int ierr = CEED_ERROR_SUCCESS; 1428 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d; 1429 1430 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis dimension must be a positive value"); 1431 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1432 CeedCheck(P > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1433 CeedCheck(Q > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1434 1435 // Get Nodes and Weights 1436 CeedCall(CeedCalloc(P * Q, &interp_1d)); 1437 CeedCall(CeedCalloc(P * Q, &grad_1d)); 1438 CeedCall(CeedCalloc(P, &nodes)); 1439 CeedCall(CeedCalloc(Q, &q_ref_1d)); 1440 CeedCall(CeedCalloc(Q, &q_weight_1d)); 1441 if (CeedLobattoQuadrature(P, nodes, NULL) != CEED_ERROR_SUCCESS) goto cleanup; 1442 switch (quad_mode) { 1443 case CEED_GAUSS: 1444 ierr = CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); 1445 break; 1446 case CEED_GAUSS_LOBATTO: 1447 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); 1448 break; 1449 } 1450 if (ierr != CEED_ERROR_SUCCESS) goto cleanup; 1451 1452 // Build B, D matrix 1453 // Fornberg, 1998 1454 for (CeedInt i = 0; i < Q; i++) { 1455 c1 = 1.0; 1456 c3 = nodes[0] - q_ref_1d[i]; 1457 interp_1d[i * P + 0] = 1.0; 1458 for (CeedInt j = 1; j < P; j++) { 1459 c2 = 1.0; 1460 c4 = c3; 1461 c3 = nodes[j] - q_ref_1d[i]; 1462 for (CeedInt k = 0; k < j; k++) { 1463 dx = nodes[j] - nodes[k]; 1464 c2 *= dx; 1465 if (k == j - 1) { 1466 grad_1d[i * P + j] = c1 * (interp_1d[i * P + k] - c4 * grad_1d[i * P + k]) / c2; 1467 interp_1d[i * P + j] = -c1 * c4 * interp_1d[i * P + k] / c2; 1468 } 1469 grad_1d[i * P + k] = (c3 * grad_1d[i * P + k] - interp_1d[i * P + k]) / dx; 1470 interp_1d[i * P + k] = c3 * interp_1d[i * P + k] / dx; 1471 } 1472 c1 = c2; 1473 } 1474 } 1475 // Pass to CeedBasisCreateTensorH1 1476 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, interp_1d, grad_1d, q_ref_1d, q_weight_1d, basis)); 1477 cleanup: 1478 CeedCall(CeedFree(&interp_1d)); 1479 CeedCall(CeedFree(&grad_1d)); 1480 CeedCall(CeedFree(&nodes)); 1481 CeedCall(CeedFree(&q_ref_1d)); 1482 CeedCall(CeedFree(&q_weight_1d)); 1483 return CEED_ERROR_SUCCESS; 1484 } 1485 1486 /** 1487 @brief Create a non tensor-product basis for \f$H^1\f$ discretizations 1488 1489 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1490 @param[in] topo Topology of element, e.g. hypercube, simplex, etc 1491 @param[in] num_comp Number of field components (1 for scalar fields) 1492 @param[in] num_nodes Total number of nodes 1493 @param[in] num_qpts Total number of quadrature points 1494 @param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal basis functions at quadrature points 1495 @param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of nodal basis functions at quadrature points 1496 @param[in] q_ref Array of length `num_qpts` * dim holding the locations of quadrature points on the reference element 1497 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 1498 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1499 1500 @return An error code: 0 - success, otherwise - failure 1501 1502 @ref User 1503 **/ 1504 int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 1505 const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { 1506 CeedInt P = num_nodes, Q = num_qpts, dim = 0; 1507 1508 if (!ceed->BasisCreateH1) { 1509 Ceed delegate; 1510 1511 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1512 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateH1"); 1513 CeedCall(CeedBasisCreateH1(delegate, topo, num_comp, num_nodes, num_qpts, interp, grad, q_ref, q_weight, basis)); 1514 CeedCall(CeedDestroy(&delegate)); 1515 return CEED_ERROR_SUCCESS; 1516 } 1517 1518 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1519 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1520 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1521 1522 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 1523 1524 CeedCall(CeedCalloc(1, basis)); 1525 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1526 (*basis)->ref_count = 1; 1527 (*basis)->is_tensor_basis = false; 1528 (*basis)->dim = dim; 1529 (*basis)->topo = topo; 1530 (*basis)->num_comp = num_comp; 1531 (*basis)->P = P; 1532 (*basis)->Q = Q; 1533 (*basis)->fe_space = CEED_FE_SPACE_H1; 1534 CeedCall(CeedCalloc(Q * dim, &(*basis)->q_ref_1d)); 1535 CeedCall(CeedCalloc(Q, &(*basis)->q_weight_1d)); 1536 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); 1537 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); 1538 CeedCall(CeedCalloc(Q * P, &(*basis)->interp)); 1539 CeedCall(CeedCalloc(dim * Q * P, &(*basis)->grad)); 1540 if (interp) memcpy((*basis)->interp, interp, Q * P * sizeof(interp[0])); 1541 if (grad) memcpy((*basis)->grad, grad, dim * Q * P * sizeof(grad[0])); 1542 CeedCall(ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, *basis)); 1543 return CEED_ERROR_SUCCESS; 1544 } 1545 1546 /** 1547 @brief Create a non tensor-product basis for \f$H(\mathrm{div})\f$ discretizations 1548 1549 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1550 @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below 1551 @param[in] num_comp Number of components (usually 1 for vectors in H(div) bases) 1552 @param[in] num_nodes Total number of nodes (DoFs per element) 1553 @param[in] num_qpts Total number of quadrature points 1554 @param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of basis functions at quadrature points 1555 @param[in] div Row-major (`num_qpts * num_nodes`) matrix expressing divergence of basis functions at quadrature points 1556 @param[in] q_ref Array of length `num_qpts` * dim holding the locations of quadrature points on the reference element 1557 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 1558 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1559 1560 @return An error code: 0 - success, otherwise - failure 1561 1562 @ref User 1563 **/ 1564 int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 1565 const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { 1566 CeedInt Q = num_qpts, P = num_nodes, dim = 0; 1567 1568 if (!ceed->BasisCreateHdiv) { 1569 Ceed delegate; 1570 1571 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1572 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHdiv"); 1573 CeedCall(CeedBasisCreateHdiv(delegate, topo, num_comp, num_nodes, num_qpts, interp, div, q_ref, q_weight, basis)); 1574 CeedCall(CeedDestroy(&delegate)); 1575 return CEED_ERROR_SUCCESS; 1576 } 1577 1578 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1579 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1580 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1581 1582 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 1583 1584 CeedCall(CeedCalloc(1, basis)); 1585 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1586 (*basis)->ref_count = 1; 1587 (*basis)->is_tensor_basis = false; 1588 (*basis)->dim = dim; 1589 (*basis)->topo = topo; 1590 (*basis)->num_comp = num_comp; 1591 (*basis)->P = P; 1592 (*basis)->Q = Q; 1593 (*basis)->fe_space = CEED_FE_SPACE_HDIV; 1594 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); 1595 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); 1596 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); 1597 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); 1598 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); 1599 CeedCall(CeedMalloc(Q * P, &(*basis)->div)); 1600 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); 1601 if (div) memcpy((*basis)->div, div, Q * P * sizeof(div[0])); 1602 CeedCall(ceed->BasisCreateHdiv(topo, dim, P, Q, interp, div, q_ref, q_weight, *basis)); 1603 return CEED_ERROR_SUCCESS; 1604 } 1605 1606 /** 1607 @brief Create a non tensor-product basis for \f$H(\mathrm{curl})\f$ discretizations 1608 1609 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1610 @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below 1611 @param[in] num_comp Number of components (usually 1 for vectors in \f$H(\mathrm{curl})\f$ bases) 1612 @param[in] num_nodes Total number of nodes (DoFs per element) 1613 @param[in] num_qpts Total number of quadrature points 1614 @param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of basis functions at quadrature points 1615 @param[in] curl Row-major (`curl_comp * num_qpts * num_nodes`, `curl_comp = 1` if `dim < 3` otherwise `curl_comp = dim`) matrix expressing curl of basis functions at quadrature points 1616 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element 1617 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 1618 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1619 1620 @return An error code: 0 - success, otherwise - failure 1621 1622 @ref User 1623 **/ 1624 int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 1625 const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { 1626 CeedInt Q = num_qpts, P = num_nodes, dim = 0, curl_comp = 0; 1627 1628 if (!ceed->BasisCreateHcurl) { 1629 Ceed delegate; 1630 1631 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1632 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHcurl"); 1633 CeedCall(CeedBasisCreateHcurl(delegate, topo, num_comp, num_nodes, num_qpts, interp, curl, q_ref, q_weight, basis)); 1634 CeedCall(CeedDestroy(&delegate)); 1635 return CEED_ERROR_SUCCESS; 1636 } 1637 1638 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1639 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1640 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1641 1642 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 1643 curl_comp = (dim < 3) ? 1 : dim; 1644 1645 CeedCall(CeedCalloc(1, basis)); 1646 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1647 (*basis)->ref_count = 1; 1648 (*basis)->is_tensor_basis = false; 1649 (*basis)->dim = dim; 1650 (*basis)->topo = topo; 1651 (*basis)->num_comp = num_comp; 1652 (*basis)->P = P; 1653 (*basis)->Q = Q; 1654 (*basis)->fe_space = CEED_FE_SPACE_HCURL; 1655 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); 1656 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); 1657 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); 1658 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); 1659 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); 1660 CeedCall(CeedMalloc(curl_comp * Q * P, &(*basis)->curl)); 1661 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); 1662 if (curl) memcpy((*basis)->curl, curl, curl_comp * Q * P * sizeof(curl[0])); 1663 CeedCall(ceed->BasisCreateHcurl(topo, dim, P, Q, interp, curl, q_ref, q_weight, *basis)); 1664 return CEED_ERROR_SUCCESS; 1665 } 1666 1667 /** 1668 @brief Create a `CeedBasis` for projection from the nodes of `basis_from` to the nodes of `basis_to`. 1669 1670 Only @ref CEED_EVAL_INTERP will be valid for the new basis, `basis_project`. 1671 For \f$H^1\f$ spaces, @ref CEED_EVAL_GRAD will also be valid. 1672 The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization. 1673 The gradient (for the \f$H^1\f$ case) is given by `grad_project = interp_to^+ * grad_from`. 1674 1675 Note: `basis_from` and `basis_to` must have compatible quadrature spaces. 1676 1677 Note: `basis_project` will have the same number of components as `basis_from`, regardless of the number of components that `basis_to` has. 1678 If `basis_from` has 3 components and `basis_to` has 5 components, then `basis_project` will have 3 components. 1679 1680 Note: If either `basis_from` or `basis_to` are non-tensor, then `basis_project` will also be non-tensor 1681 1682 @param[in] basis_from `CeedBasis` to prolong from 1683 @param[in] basis_to `CeedBasis` to prolong to 1684 @param[out] basis_project Address of the variable where the newly created `CeedBasis` will be stored 1685 1686 @return An error code: 0 - success, otherwise - failure 1687 1688 @ref User 1689 **/ 1690 int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasis *basis_project) { 1691 Ceed ceed; 1692 bool create_tensor; 1693 CeedInt dim, num_comp; 1694 CeedScalar *interp_project, *grad_project; 1695 1696 CeedCall(CeedBasisGetCeed(basis_to, &ceed)); 1697 1698 // Create projection matrix 1699 CeedCall(CeedBasisCreateProjectionMatrices(basis_from, basis_to, &interp_project, &grad_project)); 1700 1701 // Build basis 1702 { 1703 bool is_tensor_to, is_tensor_from; 1704 1705 CeedCall(CeedBasisIsTensor(basis_to, &is_tensor_to)); 1706 CeedCall(CeedBasisIsTensor(basis_from, &is_tensor_from)); 1707 create_tensor = is_tensor_from && is_tensor_to; 1708 } 1709 CeedCall(CeedBasisGetDimension(basis_to, &dim)); 1710 CeedCall(CeedBasisGetNumComponents(basis_from, &num_comp)); 1711 if (create_tensor) { 1712 CeedInt P_1d_to, P_1d_from; 1713 1714 CeedCall(CeedBasisGetNumNodes1D(basis_from, &P_1d_from)); 1715 CeedCall(CeedBasisGetNumNodes1D(basis_to, &P_1d_to)); 1716 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_from, P_1d_to, interp_project, grad_project, NULL, NULL, basis_project)); 1717 } else { 1718 // Even if basis_to and basis_from are not H1, the resulting basis is H1 for interpolation to work 1719 CeedInt num_nodes_to, num_nodes_from; 1720 CeedElemTopology topo; 1721 1722 CeedCall(CeedBasisGetTopology(basis_from, &topo)); 1723 CeedCall(CeedBasisGetNumNodes(basis_from, &num_nodes_from)); 1724 CeedCall(CeedBasisGetNumNodes(basis_to, &num_nodes_to)); 1725 CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, NULL, NULL, basis_project)); 1726 } 1727 1728 // Cleanup 1729 CeedCall(CeedFree(&interp_project)); 1730 CeedCall(CeedFree(&grad_project)); 1731 CeedCall(CeedDestroy(&ceed)); 1732 return CEED_ERROR_SUCCESS; 1733 } 1734 1735 /** 1736 @brief Copy the pointer to a `CeedBasis`. 1737 1738 Note: If the value of `*basis_copy` passed into this function is non-`NULL`, then it is assumed that `*basis_copy` is a pointer to a `CeedBasis`. 1739 This `CeedBasis` will be destroyed if `*basis_copy` is the only reference to this `CeedBasis`. 1740 1741 @param[in] basis `CeedBasis` to copy reference to 1742 @param[in,out] basis_copy Variable to store copied reference 1743 1744 @return An error code: 0 - success, otherwise - failure 1745 1746 @ref User 1747 **/ 1748 int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) { 1749 if (basis != CEED_BASIS_NONE) CeedCall(CeedBasisReference(basis)); 1750 CeedCall(CeedBasisDestroy(basis_copy)); 1751 *basis_copy = basis; 1752 return CEED_ERROR_SUCCESS; 1753 } 1754 1755 /** 1756 @brief View a `CeedBasis` 1757 1758 @param[in] basis `CeedBasis` to view 1759 @param[in] stream Stream to view to, e.g., `stdout` 1760 1761 @return An error code: 0 - success, otherwise - failure 1762 1763 @ref User 1764 **/ 1765 int CeedBasisView(CeedBasis basis, FILE *stream) { 1766 bool is_tensor_basis; 1767 CeedElemTopology topo; 1768 CeedFESpace fe_space; 1769 1770 // Basis data 1771 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 1772 CeedCall(CeedBasisGetTopology(basis, &topo)); 1773 CeedCall(CeedBasisGetFESpace(basis, &fe_space)); 1774 1775 // Print FE space and element topology of the basis 1776 fprintf(stream, "CeedBasis in a %s on a %s element\n", CeedFESpaces[fe_space], CeedElemTopologies[topo]); 1777 if (is_tensor_basis) { 1778 fprintf(stream, " P: %" CeedInt_FMT "\n Q: %" CeedInt_FMT "\n", basis->P_1d, basis->Q_1d); 1779 } else { 1780 fprintf(stream, " P: %" CeedInt_FMT "\n Q: %" CeedInt_FMT "\n", basis->P, basis->Q); 1781 } 1782 fprintf(stream, " dimension: %" CeedInt_FMT "\n field components: %" CeedInt_FMT "\n", basis->dim, basis->num_comp); 1783 // Print quadrature data, interpolation/gradient/divergence/curl of the basis 1784 if (is_tensor_basis) { // tensor basis 1785 CeedInt P_1d, Q_1d; 1786 const CeedScalar *q_ref_1d, *q_weight_1d, *interp_1d, *grad_1d; 1787 1788 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 1789 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 1790 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); 1791 CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d)); 1792 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 1793 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 1794 1795 CeedCall(CeedScalarView("qref1d", "\t% 12.8f", 1, Q_1d, q_ref_1d, stream)); 1796 CeedCall(CeedScalarView("qweight1d", "\t% 12.8f", 1, Q_1d, q_weight_1d, stream)); 1797 CeedCall(CeedScalarView("interp1d", "\t% 12.8f", Q_1d, P_1d, interp_1d, stream)); 1798 CeedCall(CeedScalarView("grad1d", "\t% 12.8f", Q_1d, P_1d, grad_1d, stream)); 1799 } else { // non-tensor basis 1800 CeedInt P, Q, dim, q_comp; 1801 const CeedScalar *q_ref, *q_weight, *interp, *grad, *div, *curl; 1802 1803 CeedCall(CeedBasisGetNumNodes(basis, &P)); 1804 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &Q)); 1805 CeedCall(CeedBasisGetDimension(basis, &dim)); 1806 CeedCall(CeedBasisGetQRef(basis, &q_ref)); 1807 CeedCall(CeedBasisGetQWeights(basis, &q_weight)); 1808 CeedCall(CeedBasisGetInterp(basis, &interp)); 1809 CeedCall(CeedBasisGetGrad(basis, &grad)); 1810 CeedCall(CeedBasisGetDiv(basis, &div)); 1811 CeedCall(CeedBasisGetCurl(basis, &curl)); 1812 1813 CeedCall(CeedScalarView("qref", "\t% 12.8f", 1, Q * dim, q_ref, stream)); 1814 CeedCall(CeedScalarView("qweight", "\t% 12.8f", 1, Q, q_weight, stream)); 1815 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp)); 1816 CeedCall(CeedScalarView("interp", "\t% 12.8f", q_comp * Q, P, interp, stream)); 1817 if (grad) { 1818 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp)); 1819 CeedCall(CeedScalarView("grad", "\t% 12.8f", q_comp * Q, P, grad, stream)); 1820 } 1821 if (div) { 1822 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp)); 1823 CeedCall(CeedScalarView("div", "\t% 12.8f", q_comp * Q, P, div, stream)); 1824 } 1825 if (curl) { 1826 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp)); 1827 CeedCall(CeedScalarView("curl", "\t% 12.8f", q_comp * Q, P, curl, stream)); 1828 } 1829 } 1830 return CEED_ERROR_SUCCESS; 1831 } 1832 1833 /** 1834 @brief Check input vector dimensions for CeedBasisApply[Add] 1835 1836 @param[in] basis `CeedBasis` to evaluate 1837 @param[in] num_elem The number of elements to apply the basis evaluation to; 1838 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1839 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points; 1840 @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes 1841 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, 1842 @ref CEED_EVAL_INTERP to use interpolated values, 1843 @ref CEED_EVAL_GRAD to use gradients, 1844 @ref CEED_EVAL_DIV to use divergence, 1845 @ref CEED_EVAL_CURL to use curl, 1846 @ref CEED_EVAL_WEIGHT to use quadrature weights 1847 @param[in] u Input `CeedVector` 1848 @param[out] v Output `CeedVector` 1849 1850 @return An error code: 0 - success, otherwise - failure 1851 1852 @ref Developer 1853 **/ 1854 static int CeedBasisApplyCheckDims(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 1855 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; 1856 CeedSize u_length = 0, v_length; 1857 1858 CeedCall(CeedBasisGetDimension(basis, &dim)); 1859 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 1860 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 1861 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 1862 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 1863 CeedCall(CeedVectorGetLength(v, &v_length)); 1864 if (u) CeedCall(CeedVectorGetLength(u, &u_length)); 1865 1866 // Check vector lengths to prevent out of bounds issues 1867 bool has_good_dims = true; 1868 switch (eval_mode) { 1869 case CEED_EVAL_NONE: 1870 case CEED_EVAL_INTERP: 1871 case CEED_EVAL_GRAD: 1872 case CEED_EVAL_DIV: 1873 case CEED_EVAL_CURL: 1874 has_good_dims = ((t_mode == CEED_TRANSPOSE && u_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_qpts * (CeedSize)q_comp && 1875 v_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes) || 1876 (t_mode == CEED_NOTRANSPOSE && v_length >= (CeedSize)num_elem * (CeedSize)num_qpts * (CeedSize)num_comp * (CeedSize)q_comp && 1877 u_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes)); 1878 break; 1879 case CEED_EVAL_WEIGHT: 1880 has_good_dims = v_length >= (CeedSize)num_elem * (CeedSize)num_qpts; 1881 break; 1882 } 1883 CeedCheck(has_good_dims, CeedBasisReturnCeed(basis), CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); 1884 return CEED_ERROR_SUCCESS; 1885 } 1886 1887 /** 1888 @brief Apply basis evaluation from nodes to quadrature points or vice versa 1889 1890 @param[in] basis `CeedBasis` to evaluate 1891 @param[in] num_elem The number of elements to apply the basis evaluation to; 1892 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1893 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points; 1894 @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes 1895 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, 1896 @ref CEED_EVAL_INTERP to use interpolated values, 1897 @ref CEED_EVAL_GRAD to use gradients, 1898 @ref CEED_EVAL_DIV to use divergence, 1899 @ref CEED_EVAL_CURL to use curl, 1900 @ref CEED_EVAL_WEIGHT to use quadrature weights 1901 @param[in] u Input `CeedVector` 1902 @param[out] v Output `CeedVector` 1903 1904 @return An error code: 0 - success, otherwise - failure 1905 1906 @ref User 1907 **/ 1908 int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 1909 CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v)); 1910 CeedCheck(basis->Apply, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedBasisApply"); 1911 CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v)); 1912 return CEED_ERROR_SUCCESS; 1913 } 1914 1915 /** 1916 @brief Apply basis evaluation from quadrature points to nodes and sum into target vector 1917 1918 @param[in] basis `CeedBasis` to evaluate 1919 @param[in] num_elem The number of elements to apply the basis evaluation to; 1920 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1921 @param[in] t_mode @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes; 1922 @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAdd()` 1923 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, 1924 @ref CEED_EVAL_INTERP to use interpolated values, 1925 @ref CEED_EVAL_GRAD to use gradients, 1926 @ref CEED_EVAL_DIV to use divergence, 1927 @ref CEED_EVAL_CURL to use curl, 1928 @ref CEED_EVAL_WEIGHT to use quadrature weights 1929 @param[in] u Input `CeedVector` 1930 @param[out] v Output `CeedVector` to sum into 1931 1932 @return An error code: 0 - success, otherwise - failure 1933 1934 @ref User 1935 **/ 1936 int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 1937 CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAdd only supports CEED_TRANSPOSE"); 1938 CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v)); 1939 CeedCheck(basis->ApplyAdd, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedBasisApplyAdd"); 1940 CeedCall(basis->ApplyAdd(basis, num_elem, t_mode, eval_mode, u, v)); 1941 return CEED_ERROR_SUCCESS; 1942 } 1943 1944 /** 1945 @brief Apply basis evaluation from nodes to arbitrary points 1946 1947 @param[in] basis `CeedBasis` to evaluate 1948 @param[in] num_elem The number of elements to apply the basis evaluation to; 1949 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1950 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 1951 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 1952 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes 1953 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 1954 @ref CEED_EVAL_GRAD to use gradients, 1955 @ref CEED_EVAL_WEIGHT to use quadrature weights 1956 @param[in] x_ref `CeedVector` holding reference coordinates of each point 1957 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 1958 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 1959 1960 @return An error code: 0 - success, otherwise - failure 1961 1962 @ref User 1963 **/ 1964 int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, 1965 CeedVector x_ref, CeedVector u, CeedVector v) { 1966 CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 1967 if (basis->ApplyAtPoints) { 1968 CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 1969 } else { 1970 CeedCall(CeedBasisApplyAtPoints_Core(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 1971 } 1972 return CEED_ERROR_SUCCESS; 1973 } 1974 1975 /** 1976 @brief Apply basis evaluation from nodes to arbitrary points and sum into target vector 1977 1978 @param[in] basis `CeedBasis` to evaluate 1979 @param[in] num_elem The number of elements to apply the basis evaluation to; 1980 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1981 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 1982 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 1983 @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAddAtPoints()` 1984 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 1985 @ref CEED_EVAL_GRAD to use gradients, 1986 @ref CEED_EVAL_WEIGHT to use quadrature weights 1987 @param[in] x_ref `CeedVector` holding reference coordinates of each point 1988 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 1989 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 1990 1991 @return An error code: 0 - success, otherwise - failure 1992 1993 @ref User 1994 **/ 1995 int CeedBasisApplyAddAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, 1996 CeedVector x_ref, CeedVector u, CeedVector v) { 1997 CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAddAtPoints only supports CEED_TRANSPOSE"); 1998 CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 1999 if (basis->ApplyAddAtPoints) { 2000 CeedCall(basis->ApplyAddAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2001 } else { 2002 CeedCall(CeedBasisApplyAtPoints_Core(basis, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2003 } 2004 return CEED_ERROR_SUCCESS; 2005 } 2006 2007 /** 2008 @brief Get the `Ceed` associated with a `CeedBasis` 2009 2010 @param[in] basis `CeedBasis` 2011 @param[out] ceed Variable to store `Ceed` 2012 2013 @return An error code: 0 - success, otherwise - failure 2014 2015 @ref Advanced 2016 **/ 2017 int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) { 2018 *ceed = NULL; 2019 CeedCall(CeedReferenceCopy(CeedBasisReturnCeed(basis), ceed)); 2020 return CEED_ERROR_SUCCESS; 2021 } 2022 2023 /** 2024 @brief Return the `Ceed` associated with a `CeedBasis` 2025 2026 @param[in] basis `CeedBasis` 2027 2028 @return `Ceed` associated with the `basis` 2029 2030 @ref Advanced 2031 **/ 2032 Ceed CeedBasisReturnCeed(CeedBasis basis) { return basis->ceed; } 2033 2034 /** 2035 @brief Get dimension for given `CeedBasis` 2036 2037 @param[in] basis `CeedBasis` 2038 @param[out] dim Variable to store dimension of basis 2039 2040 @return An error code: 0 - success, otherwise - failure 2041 2042 @ref Advanced 2043 **/ 2044 int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) { 2045 *dim = basis->dim; 2046 return CEED_ERROR_SUCCESS; 2047 } 2048 2049 /** 2050 @brief Get topology for given `CeedBasis` 2051 2052 @param[in] basis `CeedBasis` 2053 @param[out] topo Variable to store topology of basis 2054 2055 @return An error code: 0 - success, otherwise - failure 2056 2057 @ref Advanced 2058 **/ 2059 int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) { 2060 *topo = basis->topo; 2061 return CEED_ERROR_SUCCESS; 2062 } 2063 2064 /** 2065 @brief Get number of components for given `CeedBasis` 2066 2067 @param[in] basis `CeedBasis` 2068 @param[out] num_comp Variable to store number of components 2069 2070 @return An error code: 0 - success, otherwise - failure 2071 2072 @ref Advanced 2073 **/ 2074 int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) { 2075 *num_comp = basis->num_comp; 2076 return CEED_ERROR_SUCCESS; 2077 } 2078 2079 /** 2080 @brief Get total number of nodes (in `dim` dimensions) of a `CeedBasis` 2081 2082 @param[in] basis `CeedBasis` 2083 @param[out] P Variable to store number of nodes 2084 2085 @return An error code: 0 - success, otherwise - failure 2086 2087 @ref Utility 2088 **/ 2089 int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) { 2090 *P = basis->P; 2091 return CEED_ERROR_SUCCESS; 2092 } 2093 2094 /** 2095 @brief Get total number of nodes (in 1 dimension) of a `CeedBasis` 2096 2097 @param[in] basis `CeedBasis` 2098 @param[out] P_1d Variable to store number of nodes 2099 2100 @return An error code: 0 - success, otherwise - failure 2101 2102 @ref Advanced 2103 **/ 2104 int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d) { 2105 CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply P_1d for non-tensor CeedBasis"); 2106 *P_1d = basis->P_1d; 2107 return CEED_ERROR_SUCCESS; 2108 } 2109 2110 /** 2111 @brief Get total number of quadrature points (in `dim` dimensions) of a `CeedBasis` 2112 2113 @param[in] basis `CeedBasis` 2114 @param[out] Q Variable to store number of quadrature points 2115 2116 @return An error code: 0 - success, otherwise - failure 2117 2118 @ref Utility 2119 **/ 2120 int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) { 2121 *Q = basis->Q; 2122 return CEED_ERROR_SUCCESS; 2123 } 2124 2125 /** 2126 @brief Get total number of quadrature points (in 1 dimension) of a `CeedBasis` 2127 2128 @param[in] basis `CeedBasis` 2129 @param[out] Q_1d Variable to store number of quadrature points 2130 2131 @return An error code: 0 - success, otherwise - failure 2132 2133 @ref Advanced 2134 **/ 2135 int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q_1d) { 2136 CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply Q_1d for non-tensor CeedBasis"); 2137 *Q_1d = basis->Q_1d; 2138 return CEED_ERROR_SUCCESS; 2139 } 2140 2141 /** 2142 @brief Get reference coordinates of quadrature points (in `dim` dimensions) of a `CeedBasis` 2143 2144 @param[in] basis `CeedBasis` 2145 @param[out] q_ref Variable to store reference coordinates of quadrature points 2146 2147 @return An error code: 0 - success, otherwise - failure 2148 2149 @ref Advanced 2150 **/ 2151 int CeedBasisGetQRef(CeedBasis basis, const CeedScalar **q_ref) { 2152 *q_ref = basis->q_ref_1d; 2153 return CEED_ERROR_SUCCESS; 2154 } 2155 2156 /** 2157 @brief Get quadrature weights of quadrature points (in `dim` dimensions) of a `CeedBasis` 2158 2159 @param[in] basis `CeedBasis` 2160 @param[out] q_weight Variable to store quadrature weights 2161 2162 @return An error code: 0 - success, otherwise - failure 2163 2164 @ref Advanced 2165 **/ 2166 int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) { 2167 *q_weight = basis->q_weight_1d; 2168 return CEED_ERROR_SUCCESS; 2169 } 2170 2171 /** 2172 @brief Get interpolation matrix of a `CeedBasis` 2173 2174 @param[in] basis `CeedBasis` 2175 @param[out] interp Variable to store interpolation matrix 2176 2177 @return An error code: 0 - success, otherwise - failure 2178 2179 @ref Advanced 2180 **/ 2181 int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) { 2182 if (!basis->interp && basis->is_tensor_basis) { 2183 // Allocate 2184 CeedCall(CeedMalloc(basis->Q * basis->P, &basis->interp)); 2185 2186 // Initialize 2187 for (CeedInt i = 0; i < basis->Q * basis->P; i++) basis->interp[i] = 1.0; 2188 2189 // Calculate 2190 for (CeedInt d = 0; d < basis->dim; d++) { 2191 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) { 2192 for (CeedInt node = 0; node < basis->P; node++) { 2193 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 2194 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 2195 2196 basis->interp[qpt * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p]; 2197 } 2198 } 2199 } 2200 } 2201 *interp = basis->interp; 2202 return CEED_ERROR_SUCCESS; 2203 } 2204 2205 /** 2206 @brief Get 1D interpolation matrix of a tensor product `CeedBasis` 2207 2208 @param[in] basis `CeedBasis` 2209 @param[out] interp_1d Variable to store interpolation matrix 2210 2211 @return An error code: 0 - success, otherwise - failure 2212 2213 @ref Backend 2214 **/ 2215 int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_1d) { 2216 bool is_tensor_basis; 2217 2218 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 2219 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "CeedBasis is not a tensor product CeedBasis"); 2220 *interp_1d = basis->interp_1d; 2221 return CEED_ERROR_SUCCESS; 2222 } 2223 2224 /** 2225 @brief Get gradient matrix of a `CeedBasis` 2226 2227 @param[in] basis `CeedBasis` 2228 @param[out] grad Variable to store gradient matrix 2229 2230 @return An error code: 0 - success, otherwise - failure 2231 2232 @ref Advanced 2233 **/ 2234 int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad) { 2235 if (!basis->grad && basis->is_tensor_basis) { 2236 // Allocate 2237 CeedCall(CeedMalloc(basis->dim * basis->Q * basis->P, &basis->grad)); 2238 2239 // Initialize 2240 for (CeedInt i = 0; i < basis->dim * basis->Q * basis->P; i++) basis->grad[i] = 1.0; 2241 2242 // Calculate 2243 for (CeedInt d = 0; d < basis->dim; d++) { 2244 for (CeedInt i = 0; i < basis->dim; i++) { 2245 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) { 2246 for (CeedInt node = 0; node < basis->P; node++) { 2247 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 2248 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 2249 2250 if (i == d) basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->grad_1d[q * basis->P_1d + p]; 2251 else basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p]; 2252 } 2253 } 2254 } 2255 } 2256 } 2257 *grad = basis->grad; 2258 return CEED_ERROR_SUCCESS; 2259 } 2260 2261 /** 2262 @brief Get 1D gradient matrix of a tensor product `CeedBasis` 2263 2264 @param[in] basis `CeedBasis` 2265 @param[out] grad_1d Variable to store gradient matrix 2266 2267 @return An error code: 0 - success, otherwise - failure 2268 2269 @ref Advanced 2270 **/ 2271 int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d) { 2272 bool is_tensor_basis; 2273 2274 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 2275 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "CeedBasis is not a tensor product CeedBasis"); 2276 *grad_1d = basis->grad_1d; 2277 return CEED_ERROR_SUCCESS; 2278 } 2279 2280 /** 2281 @brief Get divergence matrix of a `CeedBasis` 2282 2283 @param[in] basis `CeedBasis` 2284 @param[out] div Variable to store divergence matrix 2285 2286 @return An error code: 0 - success, otherwise - failure 2287 2288 @ref Advanced 2289 **/ 2290 int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div) { 2291 *div = basis->div; 2292 return CEED_ERROR_SUCCESS; 2293 } 2294 2295 /** 2296 @brief Get curl matrix of a `CeedBasis` 2297 2298 @param[in] basis `CeedBasis` 2299 @param[out] curl Variable to store curl matrix 2300 2301 @return An error code: 0 - success, otherwise - failure 2302 2303 @ref Advanced 2304 **/ 2305 int CeedBasisGetCurl(CeedBasis basis, const CeedScalar **curl) { 2306 *curl = basis->curl; 2307 return CEED_ERROR_SUCCESS; 2308 } 2309 2310 /** 2311 @brief Destroy a @ref CeedBasis 2312 2313 @param[in,out] basis `CeedBasis` to destroy 2314 2315 @return An error code: 0 - success, otherwise - failure 2316 2317 @ref User 2318 **/ 2319 int CeedBasisDestroy(CeedBasis *basis) { 2320 if (!*basis || *basis == CEED_BASIS_NONE || --(*basis)->ref_count > 0) { 2321 *basis = NULL; 2322 return CEED_ERROR_SUCCESS; 2323 } 2324 if ((*basis)->Destroy) CeedCall((*basis)->Destroy(*basis)); 2325 CeedCall(CeedTensorContractDestroy(&(*basis)->contract)); 2326 CeedCall(CeedFree(&(*basis)->q_ref_1d)); 2327 CeedCall(CeedFree(&(*basis)->q_weight_1d)); 2328 CeedCall(CeedFree(&(*basis)->interp)); 2329 CeedCall(CeedFree(&(*basis)->interp_1d)); 2330 CeedCall(CeedFree(&(*basis)->grad)); 2331 CeedCall(CeedFree(&(*basis)->grad_1d)); 2332 CeedCall(CeedFree(&(*basis)->div)); 2333 CeedCall(CeedFree(&(*basis)->curl)); 2334 CeedCall(CeedVectorDestroy(&(*basis)->vec_chebyshev)); 2335 CeedCall(CeedBasisDestroy(&(*basis)->basis_chebyshev)); 2336 CeedCall(CeedDestroy(&(*basis)->ceed)); 2337 CeedCall(CeedFree(basis)); 2338 return CEED_ERROR_SUCCESS; 2339 } 2340 2341 /** 2342 @brief Construct a Gauss-Legendre quadrature 2343 2344 @param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-1` exactly) 2345 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]` 2346 @param[out] q_weight_1d Array of length `Q` to hold the weights 2347 2348 @return An error code: 0 - success, otherwise - failure 2349 2350 @ref Utility 2351 **/ 2352 int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) { 2353 CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0 * atan(1.0); 2354 2355 // Build q_ref_1d, q_weight_1d 2356 for (CeedInt i = 0; i <= Q / 2; i++) { 2357 // Guess 2358 xi = cos(PI * (CeedScalar)(2 * i + 1) / ((CeedScalar)(2 * Q))); 2359 // Pn(xi) 2360 P0 = 1.0; 2361 P1 = xi; 2362 P2 = 0.0; 2363 for (CeedInt j = 2; j <= Q; j++) { 2364 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2365 P0 = P1; 2366 P1 = P2; 2367 } 2368 // First Newton Step 2369 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2370 xi = xi - P2 / dP2; 2371 // Newton to convergence 2372 for (CeedInt k = 0; k < 100 && fabs(P2) > 10 * CEED_EPSILON; k++) { 2373 P0 = 1.0; 2374 P1 = xi; 2375 for (CeedInt j = 2; j <= Q; j++) { 2376 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2377 P0 = P1; 2378 P1 = P2; 2379 } 2380 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2381 xi = xi - P2 / dP2; 2382 } 2383 // Save xi, wi 2384 wi = 2.0 / ((1.0 - xi * xi) * dP2 * dP2); 2385 q_weight_1d[i] = wi; 2386 q_weight_1d[Q - 1 - i] = wi; 2387 q_ref_1d[i] = -xi; 2388 q_ref_1d[Q - 1 - i] = xi; 2389 } 2390 return CEED_ERROR_SUCCESS; 2391 } 2392 2393 /** 2394 @brief Construct a Gauss-Legendre-Lobatto quadrature 2395 2396 @param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-3` exactly) 2397 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]` 2398 @param[out] q_weight_1d Array of length `Q` to hold the weights 2399 2400 @return An error code: 0 - success, otherwise - failure 2401 2402 @ref Utility 2403 **/ 2404 int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) { 2405 CeedScalar P0, P1, P2, dP2, d2P2, xi, wi, PI = 4.0 * atan(1.0); 2406 2407 // Build q_ref_1d, q_weight_1d 2408 // Set endpoints 2409 CeedCheck(Q > 1, NULL, CEED_ERROR_DIMENSION, "Cannot create Lobatto quadrature with Q=%" CeedInt_FMT " < 2 points", Q); 2410 wi = 2.0 / ((CeedScalar)(Q * (Q - 1))); 2411 if (q_weight_1d) { 2412 q_weight_1d[0] = wi; 2413 q_weight_1d[Q - 1] = wi; 2414 } 2415 q_ref_1d[0] = -1.0; 2416 q_ref_1d[Q - 1] = 1.0; 2417 // Interior 2418 for (CeedInt i = 1; i <= (Q - 1) / 2; i++) { 2419 // Guess 2420 xi = cos(PI * (CeedScalar)(i) / (CeedScalar)(Q - 1)); 2421 // Pn(xi) 2422 P0 = 1.0; 2423 P1 = xi; 2424 P2 = 0.0; 2425 for (CeedInt j = 2; j < Q; j++) { 2426 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2427 P0 = P1; 2428 P1 = P2; 2429 } 2430 // First Newton step 2431 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2432 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi); 2433 xi = xi - dP2 / d2P2; 2434 // Newton to convergence 2435 for (CeedInt k = 0; k < 100 && fabs(dP2) > 10 * CEED_EPSILON; k++) { 2436 P0 = 1.0; 2437 P1 = xi; 2438 for (CeedInt j = 2; j < Q; j++) { 2439 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2440 P0 = P1; 2441 P1 = P2; 2442 } 2443 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2444 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi); 2445 xi = xi - dP2 / d2P2; 2446 } 2447 // Save xi, wi 2448 wi = 2.0 / (((CeedScalar)(Q * (Q - 1))) * P2 * P2); 2449 if (q_weight_1d) { 2450 q_weight_1d[i] = wi; 2451 q_weight_1d[Q - 1 - i] = wi; 2452 } 2453 q_ref_1d[i] = -xi; 2454 q_ref_1d[Q - 1 - i] = xi; 2455 } 2456 return CEED_ERROR_SUCCESS; 2457 } 2458 2459 /// @} 2460