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