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