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