1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <ceed-impl.h> 20 #include <math.h> 21 #include <stdbool.h> 22 #include <stdio.h> 23 #include <string.h> 24 25 /// @file 26 /// Implementation of CeedBasis interfaces 27 28 /// @cond DOXYGEN_SKIP 29 static struct CeedBasis_private ceed_basis_collocated; 30 /// @endcond 31 32 /// @addtogroup CeedBasisUser 33 /// @{ 34 35 /// Indicate that the quadrature points are collocated with the nodes 36 const CeedBasis CEED_BASIS_COLLOCATED = &ceed_basis_collocated; 37 38 /// @} 39 40 /// ---------------------------------------------------------------------------- 41 /// CeedBasis Library Internal Functions 42 /// ---------------------------------------------------------------------------- 43 /// @addtogroup CeedBasisDeveloper 44 /// @{ 45 46 /** 47 @brief Compute Householder reflection 48 49 Computes A = (I - b v v^T) A 50 where A is an mxn matrix indexed as A[i*row + j*col] 51 52 @param[in,out] A Matrix to apply Householder reflection to, in place 53 @param v Householder vector 54 @param b Scaling factor 55 @param m Number of rows in A 56 @param n Number of columns in A 57 @param row Row stride 58 @param col Col stride 59 60 @return An error code: 0 - success, otherwise - failure 61 62 @ref Developer 63 **/ 64 static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v, 65 CeedScalar b, CeedInt m, CeedInt n, 66 CeedInt row, CeedInt col) { 67 for (CeedInt j=0; j<n; j++) { 68 CeedScalar w = A[0*row + j*col]; 69 for (CeedInt i=1; i<m; i++) 70 w += v[i] * A[i*row + j*col]; 71 A[0*row + j*col] -= b * w; 72 for (CeedInt i=1; i<m; i++) 73 A[i*row + j*col] -= b * w * v[i]; 74 } 75 return CEED_ERROR_SUCCESS; 76 } 77 78 /** 79 @brief Apply Householder Q matrix 80 81 Compute A = Q A where Q is mxm and A is mxn. 82 83 @param[in,out] A Matrix to apply Householder Q to, in place 84 @param Q Householder Q matrix 85 @param tau Householder scaling factors 86 @param t_mode Transpose mode for application 87 @param m Number of rows in A 88 @param n Number of columns in A 89 @param k Number of elementary reflectors in Q, k<m 90 @param row Row stride in A 91 @param col Col stride in A 92 93 @return An error code: 0 - success, otherwise - failure 94 95 @ref Developer 96 **/ 97 int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q, 98 const CeedScalar *tau, CeedTransposeMode t_mode, 99 CeedInt m, CeedInt n, CeedInt k, 100 CeedInt row, CeedInt col) { 101 int ierr; 102 CeedScalar *v; 103 ierr = CeedMalloc(m, &v); CeedChk(ierr); 104 for (CeedInt ii=0; ii<k; ii++) { 105 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k-1-ii; 106 for (CeedInt j=i+1; j<m; j++) 107 v[j] = Q[j*k+i]; 108 // Apply Householder reflector (I - tau v v^T) collo_grad_1d^T 109 ierr = CeedHouseholderReflect(&A[i*row], &v[i], tau[i], m-i, n, row, col); 110 CeedChk(ierr); 111 } 112 ierr = CeedFree(&v); CeedChk(ierr); 113 return CEED_ERROR_SUCCESS; 114 } 115 116 /** 117 @brief Compute Givens rotation 118 119 Computes A = G A (or G^T A in transpose mode) 120 where A is an mxn matrix indexed as A[i*n + j*m] 121 122 @param[in,out] A Row major matrix to apply Givens rotation to, in place 123 @param c Cosine factor 124 @param s Sine factor 125 @param t_mode @ref CEED_NOTRANSPOSE to rotate the basis counter-clockwise, 126 which has the effect of rotating columns of A clockwise; 127 @ref CEED_TRANSPOSE for the opposite rotation 128 @param i First row/column to apply rotation 129 @param k Second row/column to apply rotation 130 @param m Number of rows in A 131 @param n Number of columns in A 132 133 @return An error code: 0 - success, otherwise - failure 134 135 @ref Developer 136 **/ 137 static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s, 138 CeedTransposeMode t_mode, CeedInt i, CeedInt k, 139 CeedInt m, CeedInt n) { 140 CeedInt stride_j = 1, stride_ik = m, num_its = n; 141 if (t_mode == CEED_NOTRANSPOSE) { 142 stride_j = n; stride_ik = 1; num_its = m; 143 } 144 145 // Apply rotation 146 for (CeedInt j=0; j<num_its; j++) { 147 CeedScalar tau1 = A[i*stride_ik+j*stride_j], tau2 = A[k*stride_ik+j*stride_j]; 148 A[i*stride_ik+j*stride_j] = c*tau1 - s*tau2; 149 A[k*stride_ik+j*stride_j] = s*tau1 + c*tau2; 150 } 151 return CEED_ERROR_SUCCESS; 152 } 153 154 /** 155 @brief View an array stored in a CeedBasis 156 157 @param[in] name Name of array 158 @param[in] fp_fmt Printing format 159 @param[in] m Number of rows in array 160 @param[in] n Number of columns in array 161 @param[in] a Array to be viewed 162 @param[in] stream Stream to view to, e.g., stdout 163 164 @return An error code: 0 - success, otherwise - failure 165 166 @ref Developer 167 **/ 168 static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m, 169 CeedInt n, const CeedScalar *a, FILE *stream) { 170 for (int i=0; i<m; i++) { 171 if (m > 1) 172 fprintf(stream, "%12s[%d]:", name, i); 173 else 174 fprintf(stream, "%12s:", name); 175 for (int j=0; j<n; j++) 176 fprintf(stream, fp_fmt, fabs(a[i*n+j]) > 1E-14 ? a[i*n+j] : 0); 177 fputs("\n", stream); 178 } 179 return CEED_ERROR_SUCCESS; 180 } 181 182 /// @} 183 184 /// ---------------------------------------------------------------------------- 185 /// Ceed Backend API 186 /// ---------------------------------------------------------------------------- 187 /// @addtogroup CeedBasisBackend 188 /// @{ 189 190 /** 191 @brief Return collocated grad matrix 192 193 @param basis CeedBasis 194 @param[out] collo_grad_1d Row-major (Q_1d * Q_1d) matrix expressing derivatives of 195 basis functions at quadrature points 196 197 @return An error code: 0 - success, otherwise - failure 198 199 @ref Backend 200 **/ 201 int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { 202 int i, j, k; 203 Ceed ceed; 204 CeedInt ierr, P_1d=(basis)->P_1d, Q_1d=(basis)->Q_1d; 205 CeedScalar *interp_1d, *grad_1d, *tau; 206 207 ierr = CeedMalloc(Q_1d*P_1d, &interp_1d); CeedChk(ierr); 208 ierr = CeedMalloc(Q_1d*P_1d, &grad_1d); CeedChk(ierr); 209 ierr = CeedMalloc(Q_1d, &tau); CeedChk(ierr); 210 memcpy(interp_1d, (basis)->interp_1d, Q_1d*P_1d*sizeof(basis)->interp_1d[0]); 211 memcpy(grad_1d, (basis)->grad_1d, Q_1d*P_1d*sizeof(basis)->interp_1d[0]); 212 213 // QR Factorization, interp_1d = Q R 214 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 215 ierr = CeedQRFactorization(ceed, interp_1d, tau, Q_1d, P_1d); CeedChk(ierr); 216 // Note: This function is for backend use, so all errors are terminal 217 // and we do not need to clean up memory on failure. 218 219 // Apply Rinv, collo_grad_1d = grad_1d Rinv 220 for (i=0; i<Q_1d; i++) { // Row i 221 collo_grad_1d[Q_1d*i] = grad_1d[P_1d*i]/interp_1d[0]; 222 for (j=1; j<P_1d; j++) { // Column j 223 collo_grad_1d[j+Q_1d*i] = grad_1d[j+P_1d*i]; 224 for (k=0; k<j; k++) 225 collo_grad_1d[j+Q_1d*i] -= interp_1d[j+P_1d*k]*collo_grad_1d[k+Q_1d*i]; 226 collo_grad_1d[j+Q_1d*i] /= interp_1d[j+P_1d*j]; 227 } 228 for (j=P_1d; j<Q_1d; j++) 229 collo_grad_1d[j+Q_1d*i] = 0; 230 } 231 232 // Apply Qtranspose, collo_grad = collo_grad Q_transpose 233 ierr = CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, 234 Q_1d, Q_1d, P_1d, 1, Q_1d); CeedChk(ierr); 235 236 ierr = CeedFree(&interp_1d); CeedChk(ierr); 237 ierr = CeedFree(&grad_1d); CeedChk(ierr); 238 ierr = CeedFree(&tau); CeedChk(ierr); 239 return CEED_ERROR_SUCCESS; 240 } 241 242 /** 243 @brief Get tensor status for given CeedBasis 244 245 @param basis CeedBasis 246 @param[out] is_tensor Variable to store tensor status 247 248 @return An error code: 0 - success, otherwise - failure 249 250 @ref Backend 251 **/ 252 int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor) { 253 *is_tensor = basis->tensor_basis; 254 return CEED_ERROR_SUCCESS; 255 } 256 257 /** 258 @brief Get backend data of a CeedBasis 259 260 @param basis CeedBasis 261 @param[out] data Variable to store data 262 263 @return An error code: 0 - success, otherwise - failure 264 265 @ref Backend 266 **/ 267 int CeedBasisGetData(CeedBasis basis, void *data) { 268 *(void **)data = basis->data; 269 return CEED_ERROR_SUCCESS; 270 } 271 272 /** 273 @brief Set backend data of a CeedBasis 274 275 @param[out] basis CeedBasis 276 @param data Data to set 277 278 @return An error code: 0 - success, otherwise - failure 279 280 @ref Backend 281 **/ 282 int CeedBasisSetData(CeedBasis basis, void *data) { 283 basis->data = data; 284 return CEED_ERROR_SUCCESS; 285 } 286 287 /** 288 @brief Increment the reference counter for a CeedBasis 289 290 @param basis Basis to increment the reference counter 291 292 @return An error code: 0 - success, otherwise - failure 293 294 @ref Backend 295 **/ 296 int CeedBasisReference(CeedBasis basis) { 297 basis->ref_count++; 298 return CEED_ERROR_SUCCESS; 299 } 300 301 /** 302 @brief Get dimension for given CeedElemTopology 303 304 @param topo CeedElemTopology 305 @param[out] dim Variable to store dimension of topology 306 307 @return An error code: 0 - success, otherwise - failure 308 309 @ref Backend 310 **/ 311 int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) { 312 *dim = (CeedInt) topo >> 16; 313 return CEED_ERROR_SUCCESS; 314 } 315 316 /** 317 @brief Get CeedTensorContract of a CeedBasis 318 319 @param basis CeedBasis 320 @param[out] contract Variable to store CeedTensorContract 321 322 @return An error code: 0 - success, otherwise - failure 323 324 @ref Backend 325 **/ 326 int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) { 327 *contract = basis->contract; 328 return CEED_ERROR_SUCCESS; 329 } 330 331 /** 332 @brief Set CeedTensorContract of a CeedBasis 333 334 @param[out] basis CeedBasis 335 @param contract CeedTensorContract to set 336 337 @return An error code: 0 - success, otherwise - failure 338 339 @ref Backend 340 **/ 341 int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract) { 342 int ierr; 343 basis->contract = contract; 344 ierr = CeedTensorContractReference(contract); CeedChk(ierr); 345 return CEED_ERROR_SUCCESS; 346 } 347 348 /** 349 @brief Return a reference implementation of matrix multiplication C = A B. 350 Note, this is a reference implementation for CPU CeedScalar pointers 351 that is not intended for high performance. 352 353 @param ceed A Ceed context for error handling 354 @param[in] mat_A Row-major matrix A 355 @param[in] mat_B Row-major matrix B 356 @param[out] mat_C Row-major output matrix C 357 @param m Number of rows of C 358 @param n Number of columns of C 359 @param kk Number of columns of A/rows of B 360 361 @return An error code: 0 - success, otherwise - failure 362 363 @ref Utility 364 **/ 365 int CeedMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, 366 const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, 367 CeedInt n, CeedInt kk) { 368 for (CeedInt i=0; i<m; i++) 369 for (CeedInt j=0; j<n; j++) { 370 CeedScalar sum = 0; 371 for (CeedInt k=0; k<kk; k++) 372 sum += mat_A[k+i*kk]*mat_B[j+k*n]; 373 mat_C[j+i*n] = sum; 374 } 375 return CEED_ERROR_SUCCESS; 376 } 377 378 /// @} 379 380 /// ---------------------------------------------------------------------------- 381 /// CeedBasis Public API 382 /// ---------------------------------------------------------------------------- 383 /// @addtogroup CeedBasisUser 384 /// @{ 385 386 /** 387 @brief Create a tensor-product basis for H^1 discretizations 388 389 @param ceed A Ceed object where the CeedBasis will be created 390 @param dim Topological dimension 391 @param num_comp Number of field components (1 for scalar fields) 392 @param P_1d Number of nodes in one dimension 393 @param Q_1d Number of quadrature points in one dimension 394 @param interp_1d Row-major (Q_1d * P_1d) matrix expressing the values of nodal 395 basis functions at quadrature points 396 @param grad_1d Row-major (Q_1d * P_1d) matrix expressing derivatives of nodal 397 basis functions at quadrature points 398 @param q_ref_1d Array of length Q_1d holding the locations of quadrature points 399 on the 1D reference element [-1, 1] 400 @param q_weight_1d Array of length Q_1d holding the quadrature weights on the 401 reference element 402 @param[out] basis Address of the variable where the newly created 403 CeedBasis will be stored. 404 405 @return An error code: 0 - success, otherwise - failure 406 407 @ref User 408 **/ 409 int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, 410 CeedInt P_1d, CeedInt Q_1d, 411 const CeedScalar *interp_1d, 412 const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, 413 const CeedScalar *q_weight_1d, CeedBasis *basis) { 414 int ierr; 415 416 if (!ceed->BasisCreateTensorH1) { 417 Ceed delegate; 418 ierr = CeedGetObjectDelegate(ceed, &delegate, "Basis"); CeedChk(ierr); 419 420 if (!delegate) 421 // LCOV_EXCL_START 422 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 423 "Backend does not support BasisCreateTensorH1"); 424 // LCOV_EXCL_STOP 425 426 ierr = CeedBasisCreateTensorH1(delegate, dim, num_comp, P_1d, 427 Q_1d, interp_1d, grad_1d, q_ref_1d, 428 q_weight_1d, basis); CeedChk(ierr); 429 return CEED_ERROR_SUCCESS; 430 } 431 432 if (dim<1) 433 // LCOV_EXCL_START 434 return CeedError(ceed, CEED_ERROR_DIMENSION, 435 "Basis dimension must be a positive value"); 436 // LCOV_EXCL_STOP 437 CeedElemTopology topo = dim == 1 ? CEED_TOPOLOGY_LINE 438 : dim == 2 ? CEED_TOPOLOGY_QUAD 439 : CEED_TOPOLOGY_HEX; 440 441 ierr = CeedCalloc(1, basis); CeedChk(ierr); 442 (*basis)->ceed = ceed; 443 ierr = CeedReference(ceed); CeedChk(ierr); 444 (*basis)->ref_count = 1; 445 (*basis)->tensor_basis = 1; 446 (*basis)->dim = dim; 447 (*basis)->topo = topo; 448 (*basis)->num_comp = num_comp; 449 (*basis)->P_1d = P_1d; 450 (*basis)->Q_1d = Q_1d; 451 (*basis)->P = CeedIntPow(P_1d, dim); 452 (*basis)->Q = CeedIntPow(Q_1d, dim); 453 (*basis)->Q_comp = 1; 454 (*basis)->basis_space = 1; // 1 for H^1 space 455 ierr = CeedCalloc(Q_1d, &(*basis)->q_ref_1d); CeedChk(ierr); 456 ierr = CeedCalloc(Q_1d, &(*basis)->q_weight_1d); CeedChk(ierr); 457 if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d*sizeof(q_ref_1d[0])); 458 if (q_weight_1d) memcpy((*basis)->q_weight_1d, q_weight_1d, 459 Q_1d*sizeof(q_weight_1d[0])); 460 ierr = CeedCalloc(Q_1d*P_1d, &(*basis)->interp_1d); CeedChk(ierr); 461 ierr = CeedCalloc(Q_1d*P_1d, &(*basis)->grad_1d); CeedChk(ierr); 462 if (interp_1d) memcpy((*basis)->interp_1d, interp_1d, 463 Q_1d*P_1d*sizeof(interp_1d[0])); 464 if (grad_1d) memcpy((*basis)->grad_1d, grad_1d, Q_1d*P_1d*sizeof(grad_1d[0])); 465 ierr = ceed->BasisCreateTensorH1(dim, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, 466 q_weight_1d, *basis); CeedChk(ierr); 467 return CEED_ERROR_SUCCESS; 468 } 469 470 /** 471 @brief Create a tensor-product Lagrange basis 472 473 @param ceed A Ceed object where the CeedBasis will be created 474 @param dim Topological dimension of element 475 @param num_comp Number of field components (1 for scalar fields) 476 @param P Number of Gauss-Lobatto nodes in one dimension. The 477 polynomial degree of the resulting Q_k element is k=P-1. 478 @param Q Number of quadrature points in one dimension. 479 @param quad_mode Distribution of the Q quadrature points (affects order of 480 accuracy for the quadrature) 481 @param[out] basis Address of the variable where the newly created 482 CeedBasis will be stored. 483 484 @return An error code: 0 - success, otherwise - failure 485 486 @ref User 487 **/ 488 int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, 489 CeedInt P, CeedInt Q, CeedQuadMode quad_mode, 490 CeedBasis *basis) { 491 // Allocate 492 int ierr, ierr2, i, j, k; 493 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, 494 *q_weight_1d; 495 496 if (dim<1) 497 // LCOV_EXCL_START 498 return CeedError(ceed, CEED_ERROR_DIMENSION, 499 "Basis dimension must be a positive value"); 500 // LCOV_EXCL_STOP 501 502 // Get Nodes and Weights 503 ierr = CeedCalloc(P*Q, &interp_1d); CeedChk(ierr); 504 ierr = CeedCalloc(P*Q, &grad_1d); CeedChk(ierr); 505 ierr = CeedCalloc(P, &nodes); CeedChk(ierr); 506 ierr = CeedCalloc(Q, &q_ref_1d); CeedChk(ierr); 507 ierr = CeedCalloc(Q, &q_weight_1d); CeedChk(ierr); 508 ierr = CeedLobattoQuadrature(P, nodes, NULL); 509 if (ierr) { goto cleanup; } CeedChk(ierr); 510 switch (quad_mode) { 511 case CEED_GAUSS: 512 ierr = CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); 513 break; 514 case CEED_GAUSS_LOBATTO: 515 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); 516 break; 517 } 518 if (ierr) { goto cleanup; } CeedChk(ierr); 519 520 // Build B, D matrix 521 // Fornberg, 1998 522 for (i = 0; i < Q; i++) { 523 c1 = 1.0; 524 c3 = nodes[0] - q_ref_1d[i]; 525 interp_1d[i*P+0] = 1.0; 526 for (j = 1; j < P; j++) { 527 c2 = 1.0; 528 c4 = c3; 529 c3 = nodes[j] - q_ref_1d[i]; 530 for (k = 0; k < j; k++) { 531 dx = nodes[j] - nodes[k]; 532 c2 *= dx; 533 if (k == j - 1) { 534 grad_1d[i*P + j] = c1*(interp_1d[i*P + k] - c4*grad_1d[i*P + k]) / c2; 535 interp_1d[i*P + j] = - c1*c4*interp_1d[i*P + k] / c2; 536 } 537 grad_1d[i*P + k] = (c3*grad_1d[i*P + k] - interp_1d[i*P + k]) / dx; 538 interp_1d[i*P + k] = c3*interp_1d[i*P + k] / dx; 539 } 540 c1 = c2; 541 } 542 } 543 // Pass to CeedBasisCreateTensorH1 544 ierr = CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, interp_1d, grad_1d, 545 q_ref_1d, q_weight_1d, basis); CeedChk(ierr); 546 cleanup: 547 ierr2 = CeedFree(&interp_1d); CeedChk(ierr2); 548 ierr2 = CeedFree(&grad_1d); CeedChk(ierr2); 549 ierr2 = CeedFree(&nodes); CeedChk(ierr2); 550 ierr2 = CeedFree(&q_ref_1d); CeedChk(ierr2); 551 ierr2 = CeedFree(&q_weight_1d); CeedChk(ierr2); 552 CeedChk(ierr); 553 return CEED_ERROR_SUCCESS; 554 } 555 556 /** 557 @brief Create a non tensor-product basis for H^1 discretizations 558 559 @param ceed A Ceed object where the CeedBasis will be created 560 @param topo Topology of element, e.g. hypercube, simplex, ect 561 @param num_comp Number of field components (1 for scalar fields) 562 @param num_nodes Total number of nodes 563 @param num_qpts Total number of quadrature points 564 @param interp Row-major (num_qpts * num_nodes) matrix expressing the values of 565 nodal basis functions at quadrature points 566 @param grad Row-major (num_qpts * dim * num_nodes) matrix expressing 567 derivatives of nodal basis functions at quadrature points 568 @param q_ref Array of length num_qpts holding the locations of quadrature 569 points on the reference element 570 @param q_weight Array of length num_qpts holding the quadrature weights on the 571 reference element 572 @param[out] basis Address of the variable where the newly created 573 CeedBasis will be stored. 574 575 @return An error code: 0 - success, otherwise - failure 576 577 @ref User 578 **/ 579 int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, 580 CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 581 const CeedScalar *grad, const CeedScalar *q_ref, 582 const CeedScalar *q_weight, CeedBasis *basis) { 583 int ierr; 584 CeedInt P = num_nodes, Q = num_qpts, dim = 0; 585 586 if (!ceed->BasisCreateH1) { 587 Ceed delegate; 588 ierr = CeedGetObjectDelegate(ceed, &delegate, "Basis"); CeedChk(ierr); 589 590 if (!delegate) 591 // LCOV_EXCL_START 592 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 593 "Backend does not support BasisCreateH1"); 594 // LCOV_EXCL_STOP 595 596 ierr = CeedBasisCreateH1(delegate, topo, num_comp, num_nodes, 597 num_qpts, interp, grad, q_ref, 598 q_weight, basis); CeedChk(ierr); 599 return CEED_ERROR_SUCCESS; 600 } 601 602 ierr = CeedCalloc(1, basis); CeedChk(ierr); 603 604 ierr = CeedBasisGetTopologyDimension(topo, &dim); CeedChk(ierr); 605 606 (*basis)->ceed = ceed; 607 ierr = CeedReference(ceed); CeedChk(ierr); 608 (*basis)->ref_count = 1; 609 (*basis)->tensor_basis = 0; 610 (*basis)->dim = dim; 611 (*basis)->topo = topo; 612 (*basis)->num_comp = num_comp; 613 (*basis)->P = P; 614 (*basis)->Q = Q; 615 (*basis)->Q_comp = 1; 616 (*basis)->basis_space = 1; // 1 for H^1 space 617 ierr = CeedCalloc(Q*dim, &(*basis)->q_ref_1d); CeedChk(ierr); 618 ierr = CeedCalloc(Q, &(*basis)->q_weight_1d); CeedChk(ierr); 619 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q*dim*sizeof(q_ref[0])); 620 if(q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q*sizeof(q_weight[0])); 621 ierr = CeedCalloc(Q*P, &(*basis)->interp); CeedChk(ierr); 622 ierr = CeedCalloc(dim*Q*P, &(*basis)->grad); CeedChk(ierr); 623 if(interp) memcpy((*basis)->interp, interp, Q*P*sizeof(interp[0])); 624 if(grad) memcpy((*basis)->grad, grad, dim*Q*P*sizeof(grad[0])); 625 ierr = ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, 626 q_weight, *basis); CeedChk(ierr); 627 return CEED_ERROR_SUCCESS; 628 } 629 630 /** 631 @brief Create a non tensor-product basis for H(div) discretizations 632 633 @param ceed A Ceed object where the CeedBasis will be created 634 @param topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), 635 dimension of which is used in some array sizes below 636 @param num_comp Number of components (usually 1 for vectors in H(div) bases) 637 @param num_nodes Total number of nodes (dofs per element) 638 @param num_qpts Total number of quadrature points 639 @param interp Row-major (dim*num_qpts * num_nodes) matrix expressing the values of 640 nodal basis functions at quadrature points 641 @param div Row-major (num_qpts * num_nodes) matrix expressing 642 divergence of nodal basis functions at quadrature points 643 @param q_ref Array of length num_qpts holding the locations of quadrature 644 points on the reference element 645 @param q_weight Array of length num_qpts holding the quadrature weights on the 646 reference element 647 @param[out] basis Address of the variable where the newly created 648 CeedBasis will be stored. 649 650 @return An error code: 0 - success, otherwise - failure 651 652 @ref User 653 **/ 654 int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, 655 CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 656 const CeedScalar *div, const CeedScalar *q_ref, 657 const CeedScalar *q_weight, CeedBasis *basis) { 658 int ierr; 659 CeedInt Q = num_qpts, P = num_nodes, dim = 0; 660 ierr = CeedBasisGetTopologyDimension(topo, &dim); CeedChk(ierr); 661 if (!ceed->BasisCreateHdiv) { 662 Ceed delegate; 663 ierr = CeedGetObjectDelegate(ceed, &delegate, "Basis"); CeedChk(ierr); 664 665 if (!delegate) 666 // LCOV_EXCL_START 667 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 668 "Backend does not implement BasisCreateHdiv"); 669 // LCOV_EXCL_STOP 670 671 ierr = CeedBasisCreateHdiv(delegate, topo, num_comp, num_nodes, 672 num_qpts, interp, div, q_ref, 673 q_weight, basis); CeedChk(ierr); 674 return CEED_ERROR_SUCCESS; 675 } 676 677 ierr = CeedCalloc(1, basis); CeedChk(ierr); 678 679 (*basis)->ceed = ceed; 680 ierr = CeedReference(ceed); CeedChk(ierr); 681 (*basis)->ref_count = 1; 682 (*basis)->tensor_basis = 0; 683 (*basis)->dim = dim; 684 (*basis)->topo = topo; 685 (*basis)->num_comp = num_comp; 686 (*basis)->P = P; 687 (*basis)->Q = Q; 688 (*basis)->Q_comp = dim; 689 (*basis)->basis_space = 2; // 2 for H(div) space 690 ierr = CeedMalloc(Q*dim, &(*basis)->q_ref_1d); CeedChk(ierr); 691 ierr = CeedMalloc(Q, &(*basis)->q_weight_1d); CeedChk(ierr); 692 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q*dim*sizeof(q_ref[0])); 693 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q*sizeof(q_weight[0])); 694 ierr = CeedMalloc(dim*Q*P, &(*basis)->interp); CeedChk(ierr); 695 ierr = CeedMalloc(Q*P, &(*basis)->div); CeedChk(ierr); 696 if (interp) memcpy((*basis)->interp, interp, dim*Q*P*sizeof(interp[0])); 697 if (div) memcpy((*basis)->div, div, Q*P*sizeof(div[0])); 698 ierr = ceed->BasisCreateHdiv(topo, dim, P, Q, interp, div, q_ref, 699 q_weight, *basis); CeedChk(ierr); 700 return CEED_ERROR_SUCCESS; 701 } 702 703 /** 704 @brief Copy the pointer to a CeedBasis. Both pointers should 705 be destroyed with `CeedBasisDestroy()`; 706 Note: If `*basis_copy` is non-NULL, then it is assumed that 707 `*basis_copy` is a pointer to a CeedBasis. This CeedBasis 708 will be destroyed if `*basis_copy` is the only 709 reference to this CeedBasis. 710 711 @param basis CeedBasis to copy reference to 712 @param[out] basis_copy Variable to store copied reference 713 714 @return An error code: 0 - success, otherwise - failure 715 716 @ref User 717 **/ 718 int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) { 719 int ierr; 720 721 ierr = CeedBasisReference(basis); CeedChk(ierr); 722 ierr = CeedBasisDestroy(basis_copy); CeedChk(ierr); 723 *basis_copy = basis; 724 return CEED_ERROR_SUCCESS; 725 } 726 727 /** 728 @brief View a CeedBasis 729 730 @param basis CeedBasis to view 731 @param stream Stream to view to, e.g., stdout 732 733 @return An error code: 0 - success, otherwise - failure 734 735 @ref User 736 **/ 737 int CeedBasisView(CeedBasis basis, FILE *stream) { 738 int ierr; 739 CeedFESpace FE_space = basis->basis_space; 740 CeedElemTopology topo = basis->topo; 741 // Print FE space and element topology of the basis 742 if (basis->tensor_basis) { 743 fprintf(stream, "CeedBasis (%s on a %s element): dim=%d P=%d Q=%d\n", 744 CeedFESpaces[FE_space], CeedElemTopologies[topo], 745 basis->dim, basis->P_1d, basis->Q_1d); 746 } else { 747 fprintf(stream, "CeedBasis (%s on a %s element): dim=%d P=%d Q=%d\n", 748 CeedFESpaces[FE_space], CeedElemTopologies[topo], 749 basis->dim, basis->P, basis->Q); 750 } 751 // Print quadrature data, interpolation/gradient/divergene/curl of the basis 752 if (basis->tensor_basis) { // tensor basis 753 ierr = CeedScalarView("qref1d", "\t% 12.8f", 1, basis->Q_1d, basis->q_ref_1d, 754 stream); CeedChk(ierr); 755 ierr = CeedScalarView("qweight1d", "\t% 12.8f", 1, basis->Q_1d, 756 basis->q_weight_1d, stream); CeedChk(ierr); 757 ierr = CeedScalarView("interp1d", "\t% 12.8f", basis->Q_1d, basis->P_1d, 758 basis->interp_1d, stream); CeedChk(ierr); 759 ierr = CeedScalarView("grad1d", "\t% 12.8f", basis->Q_1d, basis->P_1d, 760 basis->grad_1d, stream); CeedChk(ierr); 761 } else { // non-tensor basis 762 ierr = CeedScalarView("qref", "\t% 12.8f", 1, basis->Q*basis->dim, 763 basis->q_ref_1d, 764 stream); CeedChk(ierr); 765 ierr = CeedScalarView("qweight", "\t% 12.8f", 1, basis->Q, basis->q_weight_1d, 766 stream); CeedChk(ierr); 767 ierr = CeedScalarView("interp", "\t% 12.8f", basis->Q_comp*basis->Q, basis->P, 768 basis->interp, stream); CeedChk(ierr); 769 if (basis->grad) { 770 ierr = CeedScalarView("grad", "\t% 12.8f", basis->dim*basis->Q, basis->P, 771 basis->grad, stream); CeedChk(ierr); 772 } 773 if (basis->div) { 774 ierr = CeedScalarView("div", "\t% 12.8f", basis->Q, basis->P, 775 basis->div, stream); CeedChk(ierr); 776 } 777 } 778 return CEED_ERROR_SUCCESS; 779 } 780 781 /** 782 @brief Apply basis evaluation from nodes to quadrature points or vice versa 783 784 @param basis CeedBasis to evaluate 785 @param num_elem The number of elements to apply the basis evaluation to; 786 the backend will specify the ordering in 787 CeedElemRestrictionCreateBlocked() 788 @param t_mode \ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature 789 points, \ref CEED_TRANSPOSE to apply the transpose, mapping 790 from quadrature points to nodes 791 @param eval_mode \ref CEED_EVAL_NONE to use values directly, 792 \ref CEED_EVAL_INTERP to use interpolated values, 793 \ref CEED_EVAL_GRAD to use gradients, 794 \ref CEED_EVAL_WEIGHT to use quadrature weights. 795 @param[in] u Input CeedVector 796 @param[out] v Output CeedVector 797 798 @return An error code: 0 - success, otherwise - failure 799 800 @ref User 801 **/ 802 int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, 803 CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 804 int ierr; 805 CeedInt u_length = 0, v_length, dim, num_comp, num_nodes, num_qpts; 806 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 807 ierr = CeedBasisGetNumComponents(basis, &num_comp); CeedChk(ierr); 808 ierr = CeedBasisGetNumNodes(basis, &num_nodes); CeedChk(ierr); 809 ierr = CeedBasisGetNumQuadraturePoints(basis, &num_qpts); CeedChk(ierr); 810 ierr = CeedVectorGetLength(v, &v_length); CeedChk(ierr); 811 if (u) { 812 ierr = CeedVectorGetLength(u, &u_length); CeedChk(ierr); 813 } 814 815 if (!basis->Apply) 816 // LCOV_EXCL_START 817 return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, 818 "Backend does not support BasisApply"); 819 // LCOV_EXCL_STOP 820 821 // Check compatibility of topological and geometrical dimensions 822 if ((t_mode == CEED_TRANSPOSE && (v_length%num_nodes != 0 || 823 u_length%num_qpts != 0)) || 824 (t_mode == CEED_NOTRANSPOSE && (u_length%num_nodes != 0 || 825 v_length%num_qpts != 0))) 826 // LCOV_EXCL_START 827 return CeedError(basis->ceed, CEED_ERROR_DIMENSION, 828 "Length of input/output vectors " 829 "incompatible with basis dimensions"); 830 // LCOV_EXCL_STOP 831 832 // Check vector lengths to prevent out of bounds issues 833 bool bad_dims = false; 834 switch (eval_mode) { 835 case CEED_EVAL_NONE: 836 case CEED_EVAL_INTERP: bad_dims = 837 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 838 v_length < num_elem*num_comp*num_nodes)) || 839 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 840 u_length < num_elem*num_comp*num_nodes))); 841 break; 842 case CEED_EVAL_GRAD: bad_dims = 843 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts*dim || 844 v_length < num_elem*num_comp*num_nodes)) || 845 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp*dim || 846 u_length < num_elem*num_comp*num_nodes))); 847 break; 848 case CEED_EVAL_WEIGHT: 849 bad_dims = v_length < num_elem*num_qpts; 850 break; 851 // LCOV_EXCL_START 852 case CEED_EVAL_DIV: bad_dims = 853 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 854 v_length < num_elem*num_comp*num_nodes)) || 855 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 856 u_length < num_elem*num_comp*num_nodes))); 857 break; 858 case CEED_EVAL_CURL: bad_dims = 859 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 860 v_length < num_elem*num_comp*num_nodes)) || 861 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 862 u_length < num_elem*num_comp*num_nodes))); 863 break; 864 // LCOV_EXCL_STOP 865 } 866 if (bad_dims) 867 // LCOV_EXCL_START 868 return CeedError(basis->ceed, CEED_ERROR_DIMENSION, 869 "Input/output vectors too short for basis and evaluation mode"); 870 // LCOV_EXCL_STOP 871 872 ierr = basis->Apply(basis, num_elem, t_mode, eval_mode, u, v); CeedChk(ierr); 873 return CEED_ERROR_SUCCESS; 874 } 875 876 /** 877 @brief Get Ceed associated with a CeedBasis 878 879 @param basis CeedBasis 880 @param[out] ceed Variable to store Ceed 881 882 @return An error code: 0 - success, otherwise - failure 883 884 @ref Advanced 885 **/ 886 int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) { 887 *ceed = basis->ceed; 888 return CEED_ERROR_SUCCESS; 889 } 890 891 /** 892 @brief Get dimension for given CeedBasis 893 894 @param basis CeedBasis 895 @param[out] dim Variable to store dimension of basis 896 897 @return An error code: 0 - success, otherwise - failure 898 899 @ref Advanced 900 **/ 901 int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) { 902 *dim = basis->dim; 903 return CEED_ERROR_SUCCESS; 904 } 905 906 /** 907 @brief Get topology for given CeedBasis 908 909 @param basis CeedBasis 910 @param[out] topo Variable to store topology of basis 911 912 @return An error code: 0 - success, otherwise - failure 913 914 @ref Advanced 915 **/ 916 int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) { 917 *topo = basis->topo; 918 return CEED_ERROR_SUCCESS; 919 } 920 921 /** 922 @brief Get number of Q-vector components for given CeedBasis 923 924 @param basis CeedBasis 925 @param[out] Q_comp Variable to store number of Q-vector components of basis 926 927 @return An error code: 0 - success, otherwise - failure 928 929 @ref Advanced 930 **/ 931 int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedInt *Q_comp) { 932 *Q_comp = basis->Q_comp; 933 return CEED_ERROR_SUCCESS; 934 } 935 936 /** 937 @brief Get number of components for given CeedBasis 938 939 @param basis CeedBasis 940 @param[out] num_comp Variable to store number of components of basis 941 942 @return An error code: 0 - success, otherwise - failure 943 944 @ref Advanced 945 **/ 946 int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) { 947 *num_comp = basis->num_comp; 948 return CEED_ERROR_SUCCESS; 949 } 950 951 /** 952 @brief Get total number of nodes (in dim dimensions) of a CeedBasis 953 954 @param basis CeedBasis 955 @param[out] P Variable to store number of nodes 956 957 @return An error code: 0 - success, otherwise - failure 958 959 @ref Utility 960 **/ 961 int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) { 962 *P = basis->P; 963 return CEED_ERROR_SUCCESS; 964 } 965 966 /** 967 @brief Get total number of nodes (in 1 dimension) of a CeedBasis 968 969 @param basis CeedBasis 970 @param[out] P_1d Variable to store number of nodes 971 972 @return An error code: 0 - success, otherwise - failure 973 974 @ref Advanced 975 **/ 976 int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d) { 977 if (!basis->tensor_basis) 978 // LCOV_EXCL_START 979 return CeedError(basis->ceed, CEED_ERROR_MINOR, 980 "Cannot supply P_1d for non-tensor basis"); 981 // LCOV_EXCL_STOP 982 983 *P_1d = basis->P_1d; 984 return CEED_ERROR_SUCCESS; 985 } 986 987 /** 988 @brief Get total number of quadrature points (in dim dimensions) of a CeedBasis 989 990 @param basis CeedBasis 991 @param[out] Q Variable to store number of quadrature points 992 993 @return An error code: 0 - success, otherwise - failure 994 995 @ref Utility 996 **/ 997 int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) { 998 *Q = basis->Q; 999 return CEED_ERROR_SUCCESS; 1000 } 1001 1002 /** 1003 @brief Get total number of quadrature points (in 1 dimension) of a CeedBasis 1004 1005 @param basis CeedBasis 1006 @param[out] Q_1d Variable to store number of quadrature points 1007 1008 @return An error code: 0 - success, otherwise - failure 1009 1010 @ref Advanced 1011 **/ 1012 int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q_1d) { 1013 if (!basis->tensor_basis) 1014 // LCOV_EXCL_START 1015 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1016 "Cannot supply Q_1d for non-tensor basis"); 1017 // LCOV_EXCL_STOP 1018 1019 *Q_1d = basis->Q_1d; 1020 return CEED_ERROR_SUCCESS; 1021 } 1022 1023 /** 1024 @brief Get reference coordinates of quadrature points (in dim dimensions) 1025 of a CeedBasis 1026 1027 @param basis CeedBasis 1028 @param[out] q_ref Variable to store reference coordinates of quadrature points 1029 1030 @return An error code: 0 - success, otherwise - failure 1031 1032 @ref Advanced 1033 **/ 1034 int CeedBasisGetQRef(CeedBasis basis, const CeedScalar **q_ref) { 1035 *q_ref = basis->q_ref_1d; 1036 return CEED_ERROR_SUCCESS; 1037 } 1038 1039 /** 1040 @brief Get quadrature weights of quadrature points (in dim dimensions) 1041 of a CeedBasis 1042 1043 @param basis CeedBasis 1044 @param[out] q_weight Variable to store quadrature weights 1045 1046 @return An error code: 0 - success, otherwise - failure 1047 1048 @ref Advanced 1049 **/ 1050 int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) { 1051 *q_weight = basis->q_weight_1d; 1052 return CEED_ERROR_SUCCESS; 1053 } 1054 1055 /** 1056 @brief Get interpolation matrix of a CeedBasis 1057 1058 @param basis CeedBasis 1059 @param[out] interp Variable to store interpolation matrix 1060 1061 @return An error code: 0 - success, otherwise - failure 1062 1063 @ref Advanced 1064 **/ 1065 int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) { 1066 if (!basis->interp && basis->tensor_basis) { 1067 // Allocate 1068 int ierr; 1069 ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr); 1070 1071 // Initialize 1072 for (CeedInt i=0; i<basis->Q*basis->P; i++) 1073 basis->interp[i] = 1.0; 1074 1075 // Calculate 1076 for (CeedInt d=0; d<basis->dim; d++) 1077 for (CeedInt qpt=0; qpt<basis->Q; qpt++) 1078 for (CeedInt node=0; node<basis->P; node++) { 1079 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 1080 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 1081 basis->interp[qpt*(basis->P)+node] *= basis->interp_1d[q*basis->P_1d+p]; 1082 } 1083 } 1084 *interp = basis->interp; 1085 return CEED_ERROR_SUCCESS; 1086 } 1087 1088 /** 1089 @brief Get 1D interpolation matrix of a tensor product CeedBasis 1090 1091 @param basis CeedBasis 1092 @param[out] interp_1d Variable to store interpolation matrix 1093 1094 @return An error code: 0 - success, otherwise - failure 1095 1096 @ref Backend 1097 **/ 1098 int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_1d) { 1099 if (!basis->tensor_basis) 1100 // LCOV_EXCL_START 1101 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1102 "CeedBasis is not a tensor product basis."); 1103 // LCOV_EXCL_STOP 1104 1105 *interp_1d = basis->interp_1d; 1106 return CEED_ERROR_SUCCESS; 1107 } 1108 1109 /** 1110 @brief Get gradient matrix of a CeedBasis 1111 1112 @param basis CeedBasis 1113 @param[out] grad Variable to store gradient matrix 1114 1115 @return An error code: 0 - success, otherwise - failure 1116 1117 @ref Advanced 1118 **/ 1119 int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad) { 1120 if (!basis->grad && basis->tensor_basis) { 1121 // Allocate 1122 int ierr; 1123 ierr = CeedMalloc(basis->dim*basis->Q*basis->P, &basis->grad); 1124 CeedChk(ierr); 1125 1126 // Initialize 1127 for (CeedInt i=0; i<basis->dim*basis->Q*basis->P; i++) 1128 basis->grad[i] = 1.0; 1129 1130 // Calculate 1131 for (CeedInt d=0; d<basis->dim; d++) 1132 for (CeedInt i=0; i<basis->dim; i++) 1133 for (CeedInt qpt=0; qpt<basis->Q; qpt++) 1134 for (CeedInt node=0; node<basis->P; node++) { 1135 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 1136 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 1137 if (i == d) 1138 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *= 1139 basis->grad_1d[q*basis->P_1d+p]; 1140 else 1141 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *= 1142 basis->interp_1d[q*basis->P_1d+p]; 1143 } 1144 } 1145 *grad = basis->grad; 1146 return CEED_ERROR_SUCCESS; 1147 } 1148 1149 /** 1150 @brief Get 1D gradient matrix of a tensor product CeedBasis 1151 1152 @param basis CeedBasis 1153 @param[out] grad_1d Variable to store gradient matrix 1154 1155 @return An error code: 0 - success, otherwise - failure 1156 1157 @ref Advanced 1158 **/ 1159 int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d) { 1160 if (!basis->tensor_basis) 1161 // LCOV_EXCL_START 1162 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1163 "CeedBasis is not a tensor product basis."); 1164 // LCOV_EXCL_STOP 1165 1166 *grad_1d = basis->grad_1d; 1167 return CEED_ERROR_SUCCESS; 1168 } 1169 1170 /** 1171 @brief Get divergence matrix of a CeedBasis 1172 1173 @param basis CeedBasis 1174 @param[out] div Variable to store divergence matrix 1175 1176 @return An error code: 0 - success, otherwise - failure 1177 1178 @ref Advanced 1179 **/ 1180 int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div) { 1181 if (!basis->div) 1182 // LCOV_EXCL_START 1183 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1184 "CeedBasis does not have divergence matrix."); 1185 // LCOV_EXCL_STOP 1186 1187 *div = basis->div; 1188 return CEED_ERROR_SUCCESS; 1189 } 1190 1191 /** 1192 @brief Destroy a CeedBasis 1193 1194 @param basis CeedBasis to destroy 1195 1196 @return An error code: 0 - success, otherwise - failure 1197 1198 @ref User 1199 **/ 1200 int CeedBasisDestroy(CeedBasis *basis) { 1201 int ierr; 1202 1203 if (!*basis || --(*basis)->ref_count > 0) return CEED_ERROR_SUCCESS; 1204 if ((*basis)->Destroy) { 1205 ierr = (*basis)->Destroy(*basis); CeedChk(ierr); 1206 } 1207 if ((*basis)->contract) { 1208 ierr = CeedTensorContractDestroy(&(*basis)->contract); CeedChk(ierr); 1209 } 1210 ierr = CeedFree(&(*basis)->interp); CeedChk(ierr); 1211 ierr = CeedFree(&(*basis)->interp_1d); CeedChk(ierr); 1212 ierr = CeedFree(&(*basis)->grad); CeedChk(ierr); 1213 ierr = CeedFree(&(*basis)->div); CeedChk(ierr); 1214 ierr = CeedFree(&(*basis)->grad_1d); CeedChk(ierr); 1215 ierr = CeedFree(&(*basis)->q_ref_1d); CeedChk(ierr); 1216 ierr = CeedFree(&(*basis)->q_weight_1d); CeedChk(ierr); 1217 ierr = CeedDestroy(&(*basis)->ceed); CeedChk(ierr); 1218 ierr = CeedFree(basis); CeedChk(ierr); 1219 return CEED_ERROR_SUCCESS; 1220 } 1221 1222 /** 1223 @brief Construct a Gauss-Legendre quadrature 1224 1225 @param Q Number of quadrature points (integrates polynomials of 1226 degree 2*Q-1 exactly) 1227 @param[out] q_ref_1d Array of length Q to hold the abscissa on [-1, 1] 1228 @param[out] q_weight_1d Array of length Q to hold the weights 1229 1230 @return An error code: 0 - success, otherwise - failure 1231 1232 @ref Utility 1233 **/ 1234 int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, 1235 CeedScalar *q_weight_1d) { 1236 // Allocate 1237 CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0*atan(1.0); 1238 // Build q_ref_1d, q_weight_1d 1239 for (int i = 0; i <= Q/2; i++) { 1240 // Guess 1241 xi = cos(PI*(CeedScalar)(2*i+1)/((CeedScalar)(2*Q))); 1242 // Pn(xi) 1243 P0 = 1.0; 1244 P1 = xi; 1245 P2 = 0.0; 1246 for (int j = 2; j <= Q; j++) { 1247 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1248 P0 = P1; 1249 P1 = P2; 1250 } 1251 // First Newton Step 1252 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1253 xi = xi-P2/dP2; 1254 // Newton to convergence 1255 for (int k=0; k<100 && fabs(P2)>10*CEED_EPSILON; k++) { 1256 P0 = 1.0; 1257 P1 = xi; 1258 for (int j = 2; j <= Q; j++) { 1259 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1260 P0 = P1; 1261 P1 = P2; 1262 } 1263 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1264 xi = xi-P2/dP2; 1265 } 1266 // Save xi, wi 1267 wi = 2.0/((1.0-xi*xi)*dP2*dP2); 1268 q_weight_1d[i] = wi; 1269 q_weight_1d[Q-1-i] = wi; 1270 q_ref_1d[i] = -xi; 1271 q_ref_1d[Q-1-i]= xi; 1272 } 1273 return CEED_ERROR_SUCCESS; 1274 } 1275 1276 /** 1277 @brief Construct a Gauss-Legendre-Lobatto quadrature 1278 1279 @param Q Number of quadrature points (integrates polynomials of 1280 degree 2*Q-3 exactly) 1281 @param[out] q_ref_1d Array of length Q to hold the abscissa on [-1, 1] 1282 @param[out] q_weight_1d Array of length Q to hold the weights 1283 1284 @return An error code: 0 - success, otherwise - failure 1285 1286 @ref Utility 1287 **/ 1288 int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, 1289 CeedScalar *q_weight_1d) { 1290 // Allocate 1291 CeedScalar P0, P1, P2, dP2, d2P2, xi, wi, PI = 4.0*atan(1.0); 1292 // Build q_ref_1d, q_weight_1d 1293 // Set endpoints 1294 if (Q < 2) 1295 // LCOV_EXCL_START 1296 return CeedError(NULL, CEED_ERROR_DIMENSION, 1297 "Cannot create Lobatto quadrature with Q=%d < 2 points", Q); 1298 // LCOV_EXCL_STOP 1299 wi = 2.0/((CeedScalar)(Q*(Q-1))); 1300 if (q_weight_1d) { 1301 q_weight_1d[0] = wi; 1302 q_weight_1d[Q-1] = wi; 1303 } 1304 q_ref_1d[0] = -1.0; 1305 q_ref_1d[Q-1] = 1.0; 1306 // Interior 1307 for (int i = 1; i <= (Q-1)/2; i++) { 1308 // Guess 1309 xi = cos(PI*(CeedScalar)(i)/(CeedScalar)(Q-1)); 1310 // Pn(xi) 1311 P0 = 1.0; 1312 P1 = xi; 1313 P2 = 0.0; 1314 for (int j = 2; j < Q; j++) { 1315 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1316 P0 = P1; 1317 P1 = P2; 1318 } 1319 // First Newton step 1320 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1321 d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi); 1322 xi = xi-dP2/d2P2; 1323 // Newton to convergence 1324 for (int k=0; k<100 && fabs(dP2)>10*CEED_EPSILON; k++) { 1325 P0 = 1.0; 1326 P1 = xi; 1327 for (int j = 2; j < Q; j++) { 1328 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1329 P0 = P1; 1330 P1 = P2; 1331 } 1332 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1333 d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi); 1334 xi = xi-dP2/d2P2; 1335 } 1336 // Save xi, wi 1337 wi = 2.0/(((CeedScalar)(Q*(Q-1)))*P2*P2); 1338 if (q_weight_1d) { 1339 q_weight_1d[i] = wi; 1340 q_weight_1d[Q-1-i] = wi; 1341 } 1342 q_ref_1d[i] = -xi; 1343 q_ref_1d[Q-1-i]= xi; 1344 } 1345 return CEED_ERROR_SUCCESS; 1346 } 1347 1348 /** 1349 @brief Return QR Factorization of a matrix 1350 1351 @param ceed A Ceed context for error handling 1352 @param[in,out] mat Row-major matrix to be factorized in place 1353 @param[in,out] tau Vector of length m of scaling factors 1354 @param m Number of rows 1355 @param n Number of columns 1356 1357 @return An error code: 0 - success, otherwise - failure 1358 1359 @ref Utility 1360 **/ 1361 int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, 1362 CeedInt m, CeedInt n) { 1363 CeedScalar v[m]; 1364 1365 // Check m >= n 1366 if (n > m) 1367 // LCOV_EXCL_START 1368 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1369 "Cannot compute QR factorization with n > m"); 1370 // LCOV_EXCL_STOP 1371 1372 for (CeedInt i=0; i<n; i++) { 1373 if (i >= m-1) { // last row of matrix, no reflection needed 1374 tau[i] = 0.; 1375 break; 1376 } 1377 // Calculate Householder vector, magnitude 1378 CeedScalar sigma = 0.0; 1379 v[i] = mat[i+n*i]; 1380 for (CeedInt j=i+1; j<m; j++) { 1381 v[j] = mat[i+n*j]; 1382 sigma += v[j] * v[j]; 1383 } 1384 CeedScalar norm = sqrt(v[i]*v[i] + sigma); // norm of v[i:m] 1385 CeedScalar R_ii = -copysign(norm, v[i]); 1386 v[i] -= R_ii; 1387 // norm of v[i:m] after modification above and scaling below 1388 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1389 // tau = 2 / (norm*norm) 1390 tau[i] = 2 * v[i]*v[i] / (v[i]*v[i] + sigma); 1391 for (CeedInt j=i+1; j<m; j++) 1392 v[j] /= v[i]; 1393 1394 // Apply Householder reflector to lower right panel 1395 CeedHouseholderReflect(&mat[i*n+i+1], &v[i], tau[i], m-i, n-i-1, n, 1); 1396 // Save v 1397 mat[i+n*i] = R_ii; 1398 for (CeedInt j=i+1; j<m; j++) 1399 mat[i+n*j] = v[j]; 1400 } 1401 return CEED_ERROR_SUCCESS; 1402 } 1403 1404 /** 1405 @brief Return symmetric Schur decomposition of the symmetric matrix mat via 1406 symmetric QR factorization 1407 1408 @param ceed A Ceed context for error handling 1409 @param[in,out] mat Row-major matrix to be factorized in place 1410 @param[out] lambda Vector of length n of eigenvalues 1411 @param n Number of rows/columns 1412 1413 @return An error code: 0 - success, otherwise - failure 1414 1415 @ref Utility 1416 **/ 1417 CeedPragmaOptimizeOff 1418 int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, 1419 CeedScalar *lambda, CeedInt n) { 1420 // Check bounds for clang-tidy 1421 if (n<2) 1422 // LCOV_EXCL_START 1423 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1424 "Cannot compute symmetric Schur decomposition of scalars"); 1425 // LCOV_EXCL_STOP 1426 1427 CeedScalar v[n-1], tau[n-1], mat_T[n*n]; 1428 1429 // Copy mat to mat_T and set mat to I 1430 memcpy(mat_T, mat, n*n*sizeof(mat[0])); 1431 for (CeedInt i=0; i<n; i++) 1432 for (CeedInt j=0; j<n; j++) 1433 mat[j+n*i] = (i==j) ? 1 : 0; 1434 1435 // Reduce to tridiagonal 1436 for (CeedInt i=0; i<n-1; i++) { 1437 // Calculate Householder vector, magnitude 1438 CeedScalar sigma = 0.0; 1439 v[i] = mat_T[i+n*(i+1)]; 1440 for (CeedInt j=i+1; j<n-1; j++) { 1441 v[j] = mat_T[i+n*(j+1)]; 1442 sigma += v[j] * v[j]; 1443 } 1444 CeedScalar norm = sqrt(v[i]*v[i] + sigma); // norm of v[i:n-1] 1445 CeedScalar R_ii = -copysign(norm, v[i]); 1446 v[i] -= R_ii; 1447 // norm of v[i:m] after modification above and scaling below 1448 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1449 // tau = 2 / (norm*norm) 1450 tau[i] = i == n - 2 ? 2 : 2 * v[i]*v[i] / (v[i]*v[i] + sigma); 1451 for (CeedInt j=i+1; j<n-1; j++) 1452 v[j] /= v[i]; 1453 1454 // Update sub and super diagonal 1455 for (CeedInt j=i+2; j<n; j++) { 1456 mat_T[i+n*j] = 0; mat_T[j+n*i] = 0; 1457 } 1458 // Apply symmetric Householder reflector to lower right panel 1459 CeedHouseholderReflect(&mat_T[(i+1)+n*(i+1)], &v[i], tau[i], 1460 n-(i+1), n-(i+1), n, 1); 1461 CeedHouseholderReflect(&mat_T[(i+1)+n*(i+1)], &v[i], tau[i], 1462 n-(i+1), n-(i+1), 1, n); 1463 1464 // Save v 1465 mat_T[i+n*(i+1)] = R_ii; 1466 mat_T[(i+1)+n*i] = R_ii; 1467 for (CeedInt j=i+1; j<n-1; j++) { 1468 mat_T[i+n*(j+1)] = v[j]; 1469 } 1470 } 1471 // Backwards accumulation of Q 1472 for (CeedInt i=n-2; i>=0; i--) { 1473 if (tau[i] > 0.0) { 1474 v[i] = 1; 1475 for (CeedInt j=i+1; j<n-1; j++) { 1476 v[j] = mat_T[i+n*(j+1)]; 1477 mat_T[i+n*(j+1)] = 0; 1478 } 1479 CeedHouseholderReflect(&mat[(i+1)+n*(i+1)], &v[i], tau[i], 1480 n-(i+1), n-(i+1), n, 1); 1481 } 1482 } 1483 1484 // Reduce sub and super diagonal 1485 CeedInt p = 0, q = 0, itr = 0, max_itr = n*n*n*n; 1486 CeedScalar tol = CEED_EPSILON; 1487 1488 while (itr < max_itr) { 1489 // Update p, q, size of reduced portions of diagonal 1490 p = 0; q = 0; 1491 for (CeedInt i=n-2; i>=0; i--) { 1492 if (fabs(mat_T[i+n*(i+1)]) < tol) 1493 q += 1; 1494 else 1495 break; 1496 } 1497 for (CeedInt i=0; i<n-q-1; i++) { 1498 if (fabs(mat_T[i+n*(i+1)]) < tol) 1499 p += 1; 1500 else 1501 break; 1502 } 1503 if (q == n-1) break; // Finished reducing 1504 1505 // Reduce tridiagonal portion 1506 CeedScalar t_nn = mat_T[(n-1-q)+n*(n-1-q)], 1507 t_nnm1 = mat_T[(n-2-q)+n*(n-1-q)]; 1508 CeedScalar d = (mat_T[(n-2-q)+n*(n-2-q)] - t_nn)/2; 1509 CeedScalar mu = t_nn - t_nnm1*t_nnm1 / 1510 (d + copysign(sqrt(d*d + t_nnm1*t_nnm1), d)); 1511 CeedScalar x = mat_T[p+n*p] - mu; 1512 CeedScalar z = mat_T[p+n*(p+1)]; 1513 for (CeedInt k=p; k<n-q-1; k++) { 1514 // Compute Givens rotation 1515 CeedScalar c = 1, s = 0; 1516 if (fabs(z) > tol) { 1517 if (fabs(z) > fabs(x)) { 1518 CeedScalar tau = -x/z; 1519 s = 1/sqrt(1+tau*tau), c = s*tau; 1520 } else { 1521 CeedScalar tau = -z/x; 1522 c = 1/sqrt(1+tau*tau), s = c*tau; 1523 } 1524 } 1525 1526 // Apply Givens rotation to T 1527 CeedGivensRotation(mat_T, c, s, CEED_NOTRANSPOSE, k, k+1, n, n); 1528 CeedGivensRotation(mat_T, c, s, CEED_TRANSPOSE, k, k+1, n, n); 1529 1530 // Apply Givens rotation to Q 1531 CeedGivensRotation(mat, c, s, CEED_NOTRANSPOSE, k, k+1, n, n); 1532 1533 // Update x, z 1534 if (k < n-q-2) { 1535 x = mat_T[k+n*(k+1)]; 1536 z = mat_T[k+n*(k+2)]; 1537 } 1538 } 1539 itr++; 1540 } 1541 1542 // Save eigenvalues 1543 for (CeedInt i=0; i<n; i++) 1544 lambda[i] = mat_T[i+n*i]; 1545 1546 // Check convergence 1547 if (itr == max_itr && q < n-1) 1548 // LCOV_EXCL_START 1549 return CeedError(ceed, CEED_ERROR_MINOR, 1550 "Symmetric QR failed to converge"); 1551 // LCOV_EXCL_STOP 1552 return CEED_ERROR_SUCCESS; 1553 } 1554 CeedPragmaOptimizeOn 1555 1556 /** 1557 @brief Return Simultaneous Diagonalization of two matrices. This solves the 1558 generalized eigenvalue problem A x = lambda B x, where A and B 1559 are symmetric and B is positive definite. We generate the matrix X 1560 and vector Lambda such that X^T A X = Lambda and X^T B X = I. This 1561 is equivalent to the LAPACK routine 'sygv' with TYPE = 1. 1562 1563 @param ceed A Ceed context for error handling 1564 @param[in] mat_A Row-major matrix to be factorized with eigenvalues 1565 @param[in] mat_B Row-major matrix to be factorized to identity 1566 @param[out] mat_X Row-major orthogonal matrix 1567 @param[out] lambda Vector of length n of generalized eigenvalues 1568 @param n Number of rows/columns 1569 1570 @return An error code: 0 - success, otherwise - failure 1571 1572 @ref Utility 1573 **/ 1574 CeedPragmaOptimizeOff 1575 int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, 1576 CeedScalar *mat_B, CeedScalar *mat_X, 1577 CeedScalar *lambda, CeedInt n) { 1578 int ierr; 1579 CeedScalar *mat_C, *mat_G, *vec_D; 1580 ierr = CeedCalloc(n*n, &mat_C); CeedChk(ierr); 1581 ierr = CeedCalloc(n*n, &mat_G); CeedChk(ierr); 1582 ierr = CeedCalloc(n, &vec_D); CeedChk(ierr); 1583 1584 // Compute B = G D G^T 1585 memcpy(mat_G, mat_B, n*n*sizeof(mat_B[0])); 1586 ierr = CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n); CeedChk(ierr); 1587 1588 // Sort eigenvalues 1589 for (CeedInt i=n-1; i>=0; i--) 1590 for (CeedInt j=0; j<i; j++) { 1591 if (fabs(vec_D[j]) > fabs(vec_D[j+1])) { 1592 CeedScalar temp; 1593 temp = vec_D[j]; vec_D[j] = vec_D[j+1]; vec_D[j+1] = temp; 1594 for (CeedInt k=0; k<n; k++) { 1595 temp = mat_G[k*n+j]; mat_G[k*n+j] = mat_G[k*n+j+1]; mat_G[k*n+j+1] = temp; 1596 } 1597 } 1598 } 1599 1600 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T 1601 // = D^-1/2 G^T A G D^-1/2 1602 // -- D = D^-1/2 1603 for (CeedInt i=0; i<n; i++) 1604 vec_D[i] = 1./sqrt(vec_D[i]); 1605 // -- G = G D^-1/2 1606 // -- C = D^-1/2 G^T 1607 for (CeedInt i=0; i<n; i++) 1608 for (CeedInt j=0; j<n; j++) { 1609 mat_G[i*n+j] *= vec_D[j]; 1610 mat_C[j*n+i] = mat_G[i*n+j]; 1611 } 1612 // -- X = (D^-1/2 G^T) A 1613 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_C, 1614 (const CeedScalar *)mat_A, mat_X, n, n, n); 1615 CeedChk(ierr); 1616 // -- C = (D^-1/2 G^T A) (G D^-1/2) 1617 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_X, 1618 (const CeedScalar *)mat_G, mat_C, n, n, n); 1619 CeedChk(ierr); 1620 1621 // Compute Q^T C Q = lambda 1622 ierr = CeedSymmetricSchurDecomposition(ceed, mat_C, lambda, n); CeedChk(ierr); 1623 1624 // Sort eigenvalues 1625 for (CeedInt i=n-1; i>=0; i--) 1626 for (CeedInt j=0; j<i; j++) { 1627 if (fabs(lambda[j]) > fabs(lambda[j+1])) { 1628 CeedScalar temp; 1629 temp = lambda[j]; lambda[j] = lambda[j+1]; lambda[j+1] = temp; 1630 for (CeedInt k=0; k<n; k++) { 1631 temp = mat_C[k*n+j]; mat_C[k*n+j] = mat_C[k*n+j+1]; mat_C[k*n+j+1] = temp; 1632 } 1633 } 1634 } 1635 1636 // Set X = (G D^1/2)^-T Q 1637 // = G D^-1/2 Q 1638 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_G, 1639 (const CeedScalar *)mat_C, mat_X, n, n, n); 1640 CeedChk(ierr); 1641 1642 // Cleanup 1643 ierr = CeedFree(&mat_C); CeedChk(ierr); 1644 ierr = CeedFree(&mat_G); CeedChk(ierr); 1645 ierr = CeedFree(&vec_D); CeedChk(ierr); 1646 return CEED_ERROR_SUCCESS; 1647 } 1648 CeedPragmaOptimizeOn 1649 1650 /// @} 1651