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