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