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