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 CeedSize u_length = 0, v_length; 806 CeedInt dim, num_comp, num_nodes, num_qpts; 807 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 808 ierr = CeedBasisGetNumComponents(basis, &num_comp); CeedChk(ierr); 809 ierr = CeedBasisGetNumNodes(basis, &num_nodes); CeedChk(ierr); 810 ierr = CeedBasisGetNumQuadraturePoints(basis, &num_qpts); CeedChk(ierr); 811 ierr = CeedVectorGetLength(v, &v_length); CeedChk(ierr); 812 if (u) { 813 ierr = CeedVectorGetLength(u, &u_length); CeedChk(ierr); 814 } 815 816 if (!basis->Apply) 817 // LCOV_EXCL_START 818 return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, 819 "Backend does not support BasisApply"); 820 // LCOV_EXCL_STOP 821 822 // Check compatibility of topological and geometrical dimensions 823 if ((t_mode == CEED_TRANSPOSE && (v_length%num_nodes != 0 || 824 u_length%num_qpts != 0)) || 825 (t_mode == CEED_NOTRANSPOSE && (u_length%num_nodes != 0 || 826 v_length%num_qpts != 0))) 827 // LCOV_EXCL_START 828 return CeedError(basis->ceed, CEED_ERROR_DIMENSION, 829 "Length of input/output vectors " 830 "incompatible with basis dimensions"); 831 // LCOV_EXCL_STOP 832 833 // Check vector lengths to prevent out of bounds issues 834 bool bad_dims = false; 835 switch (eval_mode) { 836 case CEED_EVAL_NONE: 837 case CEED_EVAL_INTERP: bad_dims = 838 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 839 v_length < num_elem*num_comp*num_nodes)) || 840 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 841 u_length < num_elem*num_comp*num_nodes))); 842 break; 843 case CEED_EVAL_GRAD: bad_dims = 844 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts*dim || 845 v_length < num_elem*num_comp*num_nodes)) || 846 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp*dim || 847 u_length < num_elem*num_comp*num_nodes))); 848 break; 849 case CEED_EVAL_WEIGHT: 850 bad_dims = v_length < num_elem*num_qpts; 851 break; 852 // LCOV_EXCL_START 853 case CEED_EVAL_DIV: bad_dims = 854 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 855 v_length < num_elem*num_comp*num_nodes)) || 856 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 857 u_length < num_elem*num_comp*num_nodes))); 858 break; 859 case CEED_EVAL_CURL: bad_dims = 860 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 861 v_length < num_elem*num_comp*num_nodes)) || 862 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 863 u_length < num_elem*num_comp*num_nodes))); 864 break; 865 // LCOV_EXCL_STOP 866 } 867 if (bad_dims) 868 // LCOV_EXCL_START 869 return CeedError(basis->ceed, CEED_ERROR_DIMENSION, 870 "Input/output vectors too short for basis and evaluation mode"); 871 // LCOV_EXCL_STOP 872 873 ierr = basis->Apply(basis, num_elem, t_mode, eval_mode, u, v); CeedChk(ierr); 874 return CEED_ERROR_SUCCESS; 875 } 876 877 /** 878 @brief Get Ceed associated with a CeedBasis 879 880 @param basis CeedBasis 881 @param[out] ceed Variable to store Ceed 882 883 @return An error code: 0 - success, otherwise - failure 884 885 @ref Advanced 886 **/ 887 int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) { 888 *ceed = basis->ceed; 889 return CEED_ERROR_SUCCESS; 890 } 891 892 /** 893 @brief Get dimension for given CeedBasis 894 895 @param basis CeedBasis 896 @param[out] dim Variable to store dimension of basis 897 898 @return An error code: 0 - success, otherwise - failure 899 900 @ref Advanced 901 **/ 902 int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) { 903 *dim = basis->dim; 904 return CEED_ERROR_SUCCESS; 905 } 906 907 /** 908 @brief Get topology for given CeedBasis 909 910 @param basis CeedBasis 911 @param[out] topo Variable to store topology of basis 912 913 @return An error code: 0 - success, otherwise - failure 914 915 @ref Advanced 916 **/ 917 int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) { 918 *topo = basis->topo; 919 return CEED_ERROR_SUCCESS; 920 } 921 922 /** 923 @brief Get number of Q-vector components for given CeedBasis 924 925 @param basis CeedBasis 926 @param[out] Q_comp Variable to store number of Q-vector components of basis 927 928 @return An error code: 0 - success, otherwise - failure 929 930 @ref Advanced 931 **/ 932 int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedInt *Q_comp) { 933 *Q_comp = basis->Q_comp; 934 return CEED_ERROR_SUCCESS; 935 } 936 937 /** 938 @brief Get number of components for given CeedBasis 939 940 @param basis CeedBasis 941 @param[out] num_comp Variable to store number of components of basis 942 943 @return An error code: 0 - success, otherwise - failure 944 945 @ref Advanced 946 **/ 947 int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) { 948 *num_comp = basis->num_comp; 949 return CEED_ERROR_SUCCESS; 950 } 951 952 /** 953 @brief Get total number of nodes (in dim dimensions) of a CeedBasis 954 955 @param basis CeedBasis 956 @param[out] P Variable to store number of nodes 957 958 @return An error code: 0 - success, otherwise - failure 959 960 @ref Utility 961 **/ 962 int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) { 963 *P = basis->P; 964 return CEED_ERROR_SUCCESS; 965 } 966 967 /** 968 @brief Get total number of nodes (in 1 dimension) of a CeedBasis 969 970 @param basis CeedBasis 971 @param[out] P_1d Variable to store number of nodes 972 973 @return An error code: 0 - success, otherwise - failure 974 975 @ref Advanced 976 **/ 977 int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d) { 978 if (!basis->tensor_basis) 979 // LCOV_EXCL_START 980 return CeedError(basis->ceed, CEED_ERROR_MINOR, 981 "Cannot supply P_1d for non-tensor basis"); 982 // LCOV_EXCL_STOP 983 984 *P_1d = basis->P_1d; 985 return CEED_ERROR_SUCCESS; 986 } 987 988 /** 989 @brief Get total number of quadrature points (in dim dimensions) of a CeedBasis 990 991 @param basis CeedBasis 992 @param[out] Q Variable to store number of quadrature points 993 994 @return An error code: 0 - success, otherwise - failure 995 996 @ref Utility 997 **/ 998 int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) { 999 *Q = basis->Q; 1000 return CEED_ERROR_SUCCESS; 1001 } 1002 1003 /** 1004 @brief Get total number of quadrature points (in 1 dimension) of a CeedBasis 1005 1006 @param basis CeedBasis 1007 @param[out] Q_1d Variable to store number of quadrature points 1008 1009 @return An error code: 0 - success, otherwise - failure 1010 1011 @ref Advanced 1012 **/ 1013 int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q_1d) { 1014 if (!basis->tensor_basis) 1015 // LCOV_EXCL_START 1016 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1017 "Cannot supply Q_1d for non-tensor basis"); 1018 // LCOV_EXCL_STOP 1019 1020 *Q_1d = basis->Q_1d; 1021 return CEED_ERROR_SUCCESS; 1022 } 1023 1024 /** 1025 @brief Get reference coordinates of quadrature points (in dim dimensions) 1026 of a CeedBasis 1027 1028 @param basis CeedBasis 1029 @param[out] q_ref Variable to store reference coordinates of quadrature points 1030 1031 @return An error code: 0 - success, otherwise - failure 1032 1033 @ref Advanced 1034 **/ 1035 int CeedBasisGetQRef(CeedBasis basis, const CeedScalar **q_ref) { 1036 *q_ref = basis->q_ref_1d; 1037 return CEED_ERROR_SUCCESS; 1038 } 1039 1040 /** 1041 @brief Get quadrature weights of quadrature points (in dim dimensions) 1042 of a CeedBasis 1043 1044 @param basis CeedBasis 1045 @param[out] q_weight Variable to store quadrature weights 1046 1047 @return An error code: 0 - success, otherwise - failure 1048 1049 @ref Advanced 1050 **/ 1051 int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) { 1052 *q_weight = basis->q_weight_1d; 1053 return CEED_ERROR_SUCCESS; 1054 } 1055 1056 /** 1057 @brief Get interpolation matrix of a CeedBasis 1058 1059 @param basis CeedBasis 1060 @param[out] interp Variable to store interpolation matrix 1061 1062 @return An error code: 0 - success, otherwise - failure 1063 1064 @ref Advanced 1065 **/ 1066 int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) { 1067 if (!basis->interp && basis->tensor_basis) { 1068 // Allocate 1069 int ierr; 1070 ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr); 1071 1072 // Initialize 1073 for (CeedInt i=0; i<basis->Q*basis->P; i++) 1074 basis->interp[i] = 1.0; 1075 1076 // Calculate 1077 for (CeedInt d=0; d<basis->dim; d++) 1078 for (CeedInt qpt=0; qpt<basis->Q; qpt++) 1079 for (CeedInt node=0; node<basis->P; node++) { 1080 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 1081 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 1082 basis->interp[qpt*(basis->P)+node] *= basis->interp_1d[q*basis->P_1d+p]; 1083 } 1084 } 1085 *interp = basis->interp; 1086 return CEED_ERROR_SUCCESS; 1087 } 1088 1089 /** 1090 @brief Get 1D interpolation matrix of a tensor product CeedBasis 1091 1092 @param basis CeedBasis 1093 @param[out] interp_1d Variable to store interpolation matrix 1094 1095 @return An error code: 0 - success, otherwise - failure 1096 1097 @ref Backend 1098 **/ 1099 int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_1d) { 1100 if (!basis->tensor_basis) 1101 // LCOV_EXCL_START 1102 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1103 "CeedBasis is not a tensor product basis."); 1104 // LCOV_EXCL_STOP 1105 1106 *interp_1d = basis->interp_1d; 1107 return CEED_ERROR_SUCCESS; 1108 } 1109 1110 /** 1111 @brief Get gradient matrix of a CeedBasis 1112 1113 @param basis CeedBasis 1114 @param[out] grad Variable to store gradient matrix 1115 1116 @return An error code: 0 - success, otherwise - failure 1117 1118 @ref Advanced 1119 **/ 1120 int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad) { 1121 if (!basis->grad && basis->tensor_basis) { 1122 // Allocate 1123 int ierr; 1124 ierr = CeedMalloc(basis->dim*basis->Q*basis->P, &basis->grad); 1125 CeedChk(ierr); 1126 1127 // Initialize 1128 for (CeedInt i=0; i<basis->dim*basis->Q*basis->P; i++) 1129 basis->grad[i] = 1.0; 1130 1131 // Calculate 1132 for (CeedInt d=0; d<basis->dim; d++) 1133 for (CeedInt i=0; i<basis->dim; i++) 1134 for (CeedInt qpt=0; qpt<basis->Q; qpt++) 1135 for (CeedInt node=0; node<basis->P; node++) { 1136 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 1137 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 1138 if (i == d) 1139 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *= 1140 basis->grad_1d[q*basis->P_1d+p]; 1141 else 1142 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *= 1143 basis->interp_1d[q*basis->P_1d+p]; 1144 } 1145 } 1146 *grad = basis->grad; 1147 return CEED_ERROR_SUCCESS; 1148 } 1149 1150 /** 1151 @brief Get 1D gradient matrix of a tensor product CeedBasis 1152 1153 @param basis CeedBasis 1154 @param[out] grad_1d Variable to store gradient matrix 1155 1156 @return An error code: 0 - success, otherwise - failure 1157 1158 @ref Advanced 1159 **/ 1160 int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d) { 1161 if (!basis->tensor_basis) 1162 // LCOV_EXCL_START 1163 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1164 "CeedBasis is not a tensor product basis."); 1165 // LCOV_EXCL_STOP 1166 1167 *grad_1d = basis->grad_1d; 1168 return CEED_ERROR_SUCCESS; 1169 } 1170 1171 /** 1172 @brief Get divergence matrix of a CeedBasis 1173 1174 @param basis CeedBasis 1175 @param[out] div Variable to store divergence matrix 1176 1177 @return An error code: 0 - success, otherwise - failure 1178 1179 @ref Advanced 1180 **/ 1181 int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div) { 1182 if (!basis->div) 1183 // LCOV_EXCL_START 1184 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1185 "CeedBasis does not have divergence matrix."); 1186 // LCOV_EXCL_STOP 1187 1188 *div = basis->div; 1189 return CEED_ERROR_SUCCESS; 1190 } 1191 1192 /** 1193 @brief Destroy a CeedBasis 1194 1195 @param basis CeedBasis to destroy 1196 1197 @return An error code: 0 - success, otherwise - failure 1198 1199 @ref User 1200 **/ 1201 int CeedBasisDestroy(CeedBasis *basis) { 1202 int ierr; 1203 1204 if (!*basis || --(*basis)->ref_count > 0) return CEED_ERROR_SUCCESS; 1205 if ((*basis)->Destroy) { 1206 ierr = (*basis)->Destroy(*basis); CeedChk(ierr); 1207 } 1208 if ((*basis)->contract) { 1209 ierr = CeedTensorContractDestroy(&(*basis)->contract); CeedChk(ierr); 1210 } 1211 ierr = CeedFree(&(*basis)->interp); CeedChk(ierr); 1212 ierr = CeedFree(&(*basis)->interp_1d); CeedChk(ierr); 1213 ierr = CeedFree(&(*basis)->grad); CeedChk(ierr); 1214 ierr = CeedFree(&(*basis)->div); CeedChk(ierr); 1215 ierr = CeedFree(&(*basis)->grad_1d); CeedChk(ierr); 1216 ierr = CeedFree(&(*basis)->q_ref_1d); CeedChk(ierr); 1217 ierr = CeedFree(&(*basis)->q_weight_1d); CeedChk(ierr); 1218 ierr = CeedDestroy(&(*basis)->ceed); CeedChk(ierr); 1219 ierr = CeedFree(basis); CeedChk(ierr); 1220 return CEED_ERROR_SUCCESS; 1221 } 1222 1223 /** 1224 @brief Construct a Gauss-Legendre quadrature 1225 1226 @param Q Number of quadrature points (integrates polynomials of 1227 degree 2*Q-1 exactly) 1228 @param[out] q_ref_1d Array of length Q to hold the abscissa on [-1, 1] 1229 @param[out] q_weight_1d Array of length Q to hold the weights 1230 1231 @return An error code: 0 - success, otherwise - failure 1232 1233 @ref Utility 1234 **/ 1235 int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, 1236 CeedScalar *q_weight_1d) { 1237 // Allocate 1238 CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0*atan(1.0); 1239 // Build q_ref_1d, q_weight_1d 1240 for (int i = 0; i <= Q/2; i++) { 1241 // Guess 1242 xi = cos(PI*(CeedScalar)(2*i+1)/((CeedScalar)(2*Q))); 1243 // Pn(xi) 1244 P0 = 1.0; 1245 P1 = xi; 1246 P2 = 0.0; 1247 for (int j = 2; j <= Q; j++) { 1248 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1249 P0 = P1; 1250 P1 = P2; 1251 } 1252 // First Newton Step 1253 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1254 xi = xi-P2/dP2; 1255 // Newton to convergence 1256 for (int k=0; k<100 && fabs(P2)>10*CEED_EPSILON; k++) { 1257 P0 = 1.0; 1258 P1 = xi; 1259 for (int j = 2; j <= Q; j++) { 1260 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1261 P0 = P1; 1262 P1 = P2; 1263 } 1264 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1265 xi = xi-P2/dP2; 1266 } 1267 // Save xi, wi 1268 wi = 2.0/((1.0-xi*xi)*dP2*dP2); 1269 q_weight_1d[i] = wi; 1270 q_weight_1d[Q-1-i] = wi; 1271 q_ref_1d[i] = -xi; 1272 q_ref_1d[Q-1-i]= xi; 1273 } 1274 return CEED_ERROR_SUCCESS; 1275 } 1276 1277 /** 1278 @brief Construct a Gauss-Legendre-Lobatto quadrature 1279 1280 @param Q Number of quadrature points (integrates polynomials of 1281 degree 2*Q-3 exactly) 1282 @param[out] q_ref_1d Array of length Q to hold the abscissa on [-1, 1] 1283 @param[out] q_weight_1d Array of length Q to hold the weights 1284 1285 @return An error code: 0 - success, otherwise - failure 1286 1287 @ref Utility 1288 **/ 1289 int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, 1290 CeedScalar *q_weight_1d) { 1291 // Allocate 1292 CeedScalar P0, P1, P2, dP2, d2P2, xi, wi, PI = 4.0*atan(1.0); 1293 // Build q_ref_1d, q_weight_1d 1294 // Set endpoints 1295 if (Q < 2) 1296 // LCOV_EXCL_START 1297 return CeedError(NULL, CEED_ERROR_DIMENSION, 1298 "Cannot create Lobatto quadrature with Q=%d < 2 points", Q); 1299 // LCOV_EXCL_STOP 1300 wi = 2.0/((CeedScalar)(Q*(Q-1))); 1301 if (q_weight_1d) { 1302 q_weight_1d[0] = wi; 1303 q_weight_1d[Q-1] = wi; 1304 } 1305 q_ref_1d[0] = -1.0; 1306 q_ref_1d[Q-1] = 1.0; 1307 // Interior 1308 for (int i = 1; i <= (Q-1)/2; i++) { 1309 // Guess 1310 xi = cos(PI*(CeedScalar)(i)/(CeedScalar)(Q-1)); 1311 // Pn(xi) 1312 P0 = 1.0; 1313 P1 = xi; 1314 P2 = 0.0; 1315 for (int j = 2; j < Q; j++) { 1316 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1317 P0 = P1; 1318 P1 = P2; 1319 } 1320 // First Newton step 1321 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1322 d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi); 1323 xi = xi-dP2/d2P2; 1324 // Newton to convergence 1325 for (int k=0; k<100 && fabs(dP2)>10*CEED_EPSILON; k++) { 1326 P0 = 1.0; 1327 P1 = xi; 1328 for (int j = 2; j < Q; j++) { 1329 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1330 P0 = P1; 1331 P1 = P2; 1332 } 1333 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1334 d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi); 1335 xi = xi-dP2/d2P2; 1336 } 1337 // Save xi, wi 1338 wi = 2.0/(((CeedScalar)(Q*(Q-1)))*P2*P2); 1339 if (q_weight_1d) { 1340 q_weight_1d[i] = wi; 1341 q_weight_1d[Q-1-i] = wi; 1342 } 1343 q_ref_1d[i] = -xi; 1344 q_ref_1d[Q-1-i]= xi; 1345 } 1346 return CEED_ERROR_SUCCESS; 1347 } 1348 1349 /** 1350 @brief Return QR Factorization of a matrix 1351 1352 @param ceed A Ceed context for error handling 1353 @param[in,out] mat Row-major matrix to be factorized in place 1354 @param[in,out] tau Vector of length m of scaling factors 1355 @param m Number of rows 1356 @param n Number of columns 1357 1358 @return An error code: 0 - success, otherwise - failure 1359 1360 @ref Utility 1361 **/ 1362 int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, 1363 CeedInt m, CeedInt n) { 1364 CeedScalar v[m]; 1365 1366 // Check m >= n 1367 if (n > m) 1368 // LCOV_EXCL_START 1369 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1370 "Cannot compute QR factorization with n > m"); 1371 // LCOV_EXCL_STOP 1372 1373 for (CeedInt i=0; i<n; i++) { 1374 if (i >= m-1) { // last row of matrix, no reflection needed 1375 tau[i] = 0.; 1376 break; 1377 } 1378 // Calculate Householder vector, magnitude 1379 CeedScalar sigma = 0.0; 1380 v[i] = mat[i+n*i]; 1381 for (CeedInt j=i+1; j<m; j++) { 1382 v[j] = mat[i+n*j]; 1383 sigma += v[j] * v[j]; 1384 } 1385 CeedScalar norm = sqrt(v[i]*v[i] + sigma); // norm of v[i:m] 1386 CeedScalar R_ii = -copysign(norm, v[i]); 1387 v[i] -= R_ii; 1388 // norm of v[i:m] after modification above and scaling below 1389 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1390 // tau = 2 / (norm*norm) 1391 tau[i] = 2 * v[i]*v[i] / (v[i]*v[i] + sigma); 1392 for (CeedInt j=i+1; j<m; j++) 1393 v[j] /= v[i]; 1394 1395 // Apply Householder reflector to lower right panel 1396 CeedHouseholderReflect(&mat[i*n+i+1], &v[i], tau[i], m-i, n-i-1, n, 1); 1397 // Save v 1398 mat[i+n*i] = R_ii; 1399 for (CeedInt j=i+1; j<m; j++) 1400 mat[i+n*j] = v[j]; 1401 } 1402 return CEED_ERROR_SUCCESS; 1403 } 1404 1405 /** 1406 @brief Return symmetric Schur decomposition of the symmetric matrix mat via 1407 symmetric QR factorization 1408 1409 @param ceed A Ceed context for error handling 1410 @param[in,out] mat Row-major matrix to be factorized in place 1411 @param[out] lambda Vector of length n of eigenvalues 1412 @param n Number of rows/columns 1413 1414 @return An error code: 0 - success, otherwise - failure 1415 1416 @ref Utility 1417 **/ 1418 CeedPragmaOptimizeOff 1419 int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, 1420 CeedScalar *lambda, CeedInt n) { 1421 // Check bounds for clang-tidy 1422 if (n<2) 1423 // LCOV_EXCL_START 1424 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1425 "Cannot compute symmetric Schur decomposition of scalars"); 1426 // LCOV_EXCL_STOP 1427 1428 CeedScalar v[n-1], tau[n-1], mat_T[n*n]; 1429 1430 // Copy mat to mat_T and set mat to I 1431 memcpy(mat_T, mat, n*n*sizeof(mat[0])); 1432 for (CeedInt i=0; i<n; i++) 1433 for (CeedInt j=0; j<n; j++) 1434 mat[j+n*i] = (i==j) ? 1 : 0; 1435 1436 // Reduce to tridiagonal 1437 for (CeedInt i=0; i<n-1; i++) { 1438 // Calculate Householder vector, magnitude 1439 CeedScalar sigma = 0.0; 1440 v[i] = mat_T[i+n*(i+1)]; 1441 for (CeedInt j=i+1; j<n-1; j++) { 1442 v[j] = mat_T[i+n*(j+1)]; 1443 sigma += v[j] * v[j]; 1444 } 1445 CeedScalar norm = sqrt(v[i]*v[i] + sigma); // norm of v[i:n-1] 1446 CeedScalar R_ii = -copysign(norm, v[i]); 1447 v[i] -= R_ii; 1448 // norm of v[i:m] after modification above and scaling below 1449 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1450 // tau = 2 / (norm*norm) 1451 tau[i] = i == n - 2 ? 2 : 2 * v[i]*v[i] / (v[i]*v[i] + sigma); 1452 for (CeedInt j=i+1; j<n-1; j++) 1453 v[j] /= v[i]; 1454 1455 // Update sub and super diagonal 1456 for (CeedInt j=i+2; j<n; j++) { 1457 mat_T[i+n*j] = 0; mat_T[j+n*i] = 0; 1458 } 1459 // Apply symmetric Householder reflector to lower right panel 1460 CeedHouseholderReflect(&mat_T[(i+1)+n*(i+1)], &v[i], tau[i], 1461 n-(i+1), n-(i+1), n, 1); 1462 CeedHouseholderReflect(&mat_T[(i+1)+n*(i+1)], &v[i], tau[i], 1463 n-(i+1), n-(i+1), 1, n); 1464 1465 // Save v 1466 mat_T[i+n*(i+1)] = R_ii; 1467 mat_T[(i+1)+n*i] = R_ii; 1468 for (CeedInt j=i+1; j<n-1; j++) { 1469 mat_T[i+n*(j+1)] = v[j]; 1470 } 1471 } 1472 // Backwards accumulation of Q 1473 for (CeedInt i=n-2; i>=0; i--) { 1474 if (tau[i] > 0.0) { 1475 v[i] = 1; 1476 for (CeedInt j=i+1; j<n-1; j++) { 1477 v[j] = mat_T[i+n*(j+1)]; 1478 mat_T[i+n*(j+1)] = 0; 1479 } 1480 CeedHouseholderReflect(&mat[(i+1)+n*(i+1)], &v[i], tau[i], 1481 n-(i+1), n-(i+1), n, 1); 1482 } 1483 } 1484 1485 // Reduce sub and super diagonal 1486 CeedInt p = 0, q = 0, itr = 0, max_itr = n*n*n*n; 1487 CeedScalar tol = CEED_EPSILON; 1488 1489 while (itr < max_itr) { 1490 // Update p, q, size of reduced portions of diagonal 1491 p = 0; q = 0; 1492 for (CeedInt i=n-2; i>=0; i--) { 1493 if (fabs(mat_T[i+n*(i+1)]) < tol) 1494 q += 1; 1495 else 1496 break; 1497 } 1498 for (CeedInt i=0; i<n-q-1; i++) { 1499 if (fabs(mat_T[i+n*(i+1)]) < tol) 1500 p += 1; 1501 else 1502 break; 1503 } 1504 if (q == n-1) break; // Finished reducing 1505 1506 // Reduce tridiagonal portion 1507 CeedScalar t_nn = mat_T[(n-1-q)+n*(n-1-q)], 1508 t_nnm1 = mat_T[(n-2-q)+n*(n-1-q)]; 1509 CeedScalar d = (mat_T[(n-2-q)+n*(n-2-q)] - t_nn)/2; 1510 CeedScalar mu = t_nn - t_nnm1*t_nnm1 / 1511 (d + copysign(sqrt(d*d + t_nnm1*t_nnm1), d)); 1512 CeedScalar x = mat_T[p+n*p] - mu; 1513 CeedScalar z = mat_T[p+n*(p+1)]; 1514 for (CeedInt k=p; k<n-q-1; k++) { 1515 // Compute Givens rotation 1516 CeedScalar c = 1, s = 0; 1517 if (fabs(z) > tol) { 1518 if (fabs(z) > fabs(x)) { 1519 CeedScalar tau = -x/z; 1520 s = 1/sqrt(1+tau*tau), c = s*tau; 1521 } else { 1522 CeedScalar tau = -z/x; 1523 c = 1/sqrt(1+tau*tau), s = c*tau; 1524 } 1525 } 1526 1527 // Apply Givens rotation to T 1528 CeedGivensRotation(mat_T, c, s, CEED_NOTRANSPOSE, k, k+1, n, n); 1529 CeedGivensRotation(mat_T, c, s, CEED_TRANSPOSE, k, k+1, n, n); 1530 1531 // Apply Givens rotation to Q 1532 CeedGivensRotation(mat, c, s, CEED_NOTRANSPOSE, k, k+1, n, n); 1533 1534 // Update x, z 1535 if (k < n-q-2) { 1536 x = mat_T[k+n*(k+1)]; 1537 z = mat_T[k+n*(k+2)]; 1538 } 1539 } 1540 itr++; 1541 } 1542 1543 // Save eigenvalues 1544 for (CeedInt i=0; i<n; i++) 1545 lambda[i] = mat_T[i+n*i]; 1546 1547 // Check convergence 1548 if (itr == max_itr && q < n-1) 1549 // LCOV_EXCL_START 1550 return CeedError(ceed, CEED_ERROR_MINOR, 1551 "Symmetric QR failed to converge"); 1552 // LCOV_EXCL_STOP 1553 return CEED_ERROR_SUCCESS; 1554 } 1555 CeedPragmaOptimizeOn 1556 1557 /** 1558 @brief Return Simultaneous Diagonalization of two matrices. This solves the 1559 generalized eigenvalue problem A x = lambda B x, where A and B 1560 are symmetric and B is positive definite. We generate the matrix X 1561 and vector Lambda such that X^T A X = Lambda and X^T B X = I. This 1562 is equivalent to the LAPACK routine 'sygv' with TYPE = 1. 1563 1564 @param ceed A Ceed context for error handling 1565 @param[in] mat_A Row-major matrix to be factorized with eigenvalues 1566 @param[in] mat_B Row-major matrix to be factorized to identity 1567 @param[out] mat_X Row-major orthogonal matrix 1568 @param[out] lambda Vector of length n of generalized eigenvalues 1569 @param n Number of rows/columns 1570 1571 @return An error code: 0 - success, otherwise - failure 1572 1573 @ref Utility 1574 **/ 1575 CeedPragmaOptimizeOff 1576 int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, 1577 CeedScalar *mat_B, CeedScalar *mat_X, 1578 CeedScalar *lambda, CeedInt n) { 1579 int ierr; 1580 CeedScalar *mat_C, *mat_G, *vec_D; 1581 ierr = CeedCalloc(n*n, &mat_C); CeedChk(ierr); 1582 ierr = CeedCalloc(n*n, &mat_G); CeedChk(ierr); 1583 ierr = CeedCalloc(n, &vec_D); CeedChk(ierr); 1584 1585 // Compute B = G D G^T 1586 memcpy(mat_G, mat_B, n*n*sizeof(mat_B[0])); 1587 ierr = CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n); CeedChk(ierr); 1588 1589 // Sort eigenvalues 1590 for (CeedInt i=n-1; i>=0; i--) 1591 for (CeedInt j=0; j<i; j++) { 1592 if (fabs(vec_D[j]) > fabs(vec_D[j+1])) { 1593 CeedScalar temp; 1594 temp = vec_D[j]; vec_D[j] = vec_D[j+1]; vec_D[j+1] = temp; 1595 for (CeedInt k=0; k<n; k++) { 1596 temp = mat_G[k*n+j]; mat_G[k*n+j] = mat_G[k*n+j+1]; mat_G[k*n+j+1] = temp; 1597 } 1598 } 1599 } 1600 1601 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T 1602 // = D^-1/2 G^T A G D^-1/2 1603 // -- D = D^-1/2 1604 for (CeedInt i=0; i<n; i++) 1605 vec_D[i] = 1./sqrt(vec_D[i]); 1606 // -- G = G D^-1/2 1607 // -- C = D^-1/2 G^T 1608 for (CeedInt i=0; i<n; i++) 1609 for (CeedInt j=0; j<n; j++) { 1610 mat_G[i*n+j] *= vec_D[j]; 1611 mat_C[j*n+i] = mat_G[i*n+j]; 1612 } 1613 // -- X = (D^-1/2 G^T) A 1614 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_C, 1615 (const CeedScalar *)mat_A, mat_X, n, n, n); 1616 CeedChk(ierr); 1617 // -- C = (D^-1/2 G^T A) (G D^-1/2) 1618 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_X, 1619 (const CeedScalar *)mat_G, mat_C, n, n, n); 1620 CeedChk(ierr); 1621 1622 // Compute Q^T C Q = lambda 1623 ierr = CeedSymmetricSchurDecomposition(ceed, mat_C, lambda, n); CeedChk(ierr); 1624 1625 // Sort eigenvalues 1626 for (CeedInt i=n-1; i>=0; i--) 1627 for (CeedInt j=0; j<i; j++) { 1628 if (fabs(lambda[j]) > fabs(lambda[j+1])) { 1629 CeedScalar temp; 1630 temp = lambda[j]; lambda[j] = lambda[j+1]; lambda[j+1] = temp; 1631 for (CeedInt k=0; k<n; k++) { 1632 temp = mat_C[k*n+j]; mat_C[k*n+j] = mat_C[k*n+j+1]; mat_C[k*n+j+1] = temp; 1633 } 1634 } 1635 } 1636 1637 // Set X = (G D^1/2)^-T Q 1638 // = G D^-1/2 Q 1639 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_G, 1640 (const CeedScalar *)mat_C, mat_X, n, n, n); 1641 CeedChk(ierr); 1642 1643 // Cleanup 1644 ierr = CeedFree(&mat_C); CeedChk(ierr); 1645 ierr = CeedFree(&mat_G); CeedChk(ierr); 1646 ierr = CeedFree(&vec_D); CeedChk(ierr); 1647 return CEED_ERROR_SUCCESS; 1648 } 1649 CeedPragmaOptimizeOn 1650 1651 /// @} 1652