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