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