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