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