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