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