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