1 // Copyright (c) 2017-2025, 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_none; 21 /// @endcond 22 23 /// @addtogroup CeedBasisUser 24 /// @{ 25 26 /// Argument for @ref CeedOperatorSetField() indicating that the field does not require a `CeedBasis` 27 const CeedBasis CEED_BASIS_NONE = &ceed_basis_none; 28 29 /// @} 30 31 /// ---------------------------------------------------------------------------- 32 /// CeedBasis Library Internal Functions 33 /// ---------------------------------------------------------------------------- 34 /// @addtogroup CeedBasisDeveloper 35 /// @{ 36 37 /** 38 @brief Compute Chebyshev polynomial values at a point 39 40 @param[in] x Coordinate to evaluate Chebyshev polynomials at 41 @param[in] n Number of Chebyshev polynomials to evaluate, `n >= 2` 42 @param[out] chebyshev_x Array of Chebyshev polynomial values 43 44 @return An error code: 0 - success, otherwise - failure 45 46 @ref Developer 47 **/ 48 static int CeedChebyshevPolynomialsAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_x) { 49 chebyshev_x[0] = 1.0; 50 chebyshev_x[1] = 2 * x; 51 for (CeedInt i = 2; i < n; i++) chebyshev_x[i] = 2 * x * chebyshev_x[i - 1] - chebyshev_x[i - 2]; 52 return CEED_ERROR_SUCCESS; 53 } 54 55 /** 56 @brief Compute values of the derivative of Chebyshev polynomials at a point 57 58 @param[in] x Coordinate to evaluate derivative of Chebyshev polynomials at 59 @param[in] n Number of Chebyshev polynomials to evaluate, `n >= 2` 60 @param[out] chebyshev_dx Array of Chebyshev polynomial derivative values 61 62 @return An error code: 0 - success, otherwise - failure 63 64 @ref Developer 65 **/ 66 static int CeedChebyshevDerivativeAtPoint(CeedScalar x, CeedInt n, CeedScalar *chebyshev_dx) { 67 CeedScalar chebyshev_x[3]; 68 69 chebyshev_x[1] = 1.0; 70 chebyshev_x[2] = 2 * x; 71 chebyshev_dx[0] = 0.0; 72 chebyshev_dx[1] = 2.0; 73 for (CeedInt i = 2; i < n; i++) { 74 chebyshev_x[0] = chebyshev_x[1]; 75 chebyshev_x[1] = chebyshev_x[2]; 76 chebyshev_x[2] = 2 * x * chebyshev_x[1] - chebyshev_x[0]; 77 chebyshev_dx[i] = 2 * x * chebyshev_dx[i - 1] + 2 * chebyshev_x[1] - chebyshev_dx[i - 2]; 78 } 79 return CEED_ERROR_SUCCESS; 80 } 81 82 /** 83 @brief Compute Householder reflection. 84 85 Computes \f$A = (I - b v v^T) A\f$, where \f$A\f$ is an \f$m \times n\f$ matrix indexed as `A[i*row + j*col]`. 86 87 @param[in,out] A Matrix to apply Householder reflection to, in place 88 @param[in] v Householder vector 89 @param[in] b Scaling factor 90 @param[in] m Number of rows in `A` 91 @param[in] n Number of columns in `A` 92 @param[in] row Row stride 93 @param[in] col Col stride 94 95 @return An error code: 0 - success, otherwise - failure 96 97 @ref Developer 98 **/ 99 static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v, CeedScalar b, CeedInt m, CeedInt n, CeedInt row, CeedInt col) { 100 for (CeedInt j = 0; j < n; j++) { 101 CeedScalar w = A[0 * row + j * col]; 102 103 for (CeedInt i = 1; i < m; i++) w += v[i] * A[i * row + j * col]; 104 A[0 * row + j * col] -= b * w; 105 for (CeedInt i = 1; i < m; i++) A[i * row + j * col] -= b * w * v[i]; 106 } 107 return CEED_ERROR_SUCCESS; 108 } 109 110 /** 111 @brief Compute Givens rotation 112 113 Computes \f$A = G A\f$ (or \f$G^T A\f$ in transpose mode), where \f$A\f$ is an \f$m \times n\f$ matrix indexed as `A[i*n + j*m]`. 114 115 @param[in,out] A Row major matrix to apply Givens rotation to, in place 116 @param[in] c Cosine factor 117 @param[in] s Sine factor 118 @param[in] t_mode @ref CEED_NOTRANSPOSE to rotate the basis counter-clockwise, which has the effect of rotating columns of `A` clockwise; 119 @ref CEED_TRANSPOSE for the opposite rotation 120 @param[in] i First row/column to apply rotation 121 @param[in] k Second row/column to apply rotation 122 @param[in] m Number of rows in `A` 123 @param[in] n Number of columns in `A` 124 125 @return An error code: 0 - success, otherwise - failure 126 127 @ref Developer 128 **/ 129 static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s, CeedTransposeMode t_mode, CeedInt i, CeedInt k, CeedInt m, CeedInt n) { 130 CeedInt stride_j = 1, stride_ik = m, num_its = n; 131 132 if (t_mode == CEED_NOTRANSPOSE) { 133 stride_j = n; 134 stride_ik = 1; 135 num_its = m; 136 } 137 138 // Apply rotation 139 for (CeedInt j = 0; j < num_its; j++) { 140 CeedScalar tau1 = A[i * stride_ik + j * stride_j], tau2 = A[k * stride_ik + j * stride_j]; 141 142 A[i * stride_ik + j * stride_j] = c * tau1 - s * tau2; 143 A[k * stride_ik + j * stride_j] = s * tau1 + c * tau2; 144 } 145 return CEED_ERROR_SUCCESS; 146 } 147 148 /** 149 @brief View an array stored in a `CeedBasis` 150 151 @param[in] name Name of array 152 @param[in] fp_fmt Printing format 153 @param[in] m Number of rows in array 154 @param[in] n Number of columns in array 155 @param[in] a Array to be viewed 156 @param[in] stream Stream to view to, e.g., `stdout` 157 158 @return An error code: 0 - success, otherwise - failure 159 160 @ref Developer 161 **/ 162 static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m, CeedInt n, const CeedScalar *a, FILE *stream) { 163 if (m > 1) { 164 fprintf(stream, " %s:\n", name); 165 } else { 166 char padded_name[12]; 167 168 snprintf(padded_name, 11, "%s:", name); 169 fprintf(stream, " %-10s", padded_name); 170 } 171 for (CeedInt i = 0; i < m; i++) { 172 if (m > 1) fprintf(stream, " [%" CeedInt_FMT "]", i); 173 for (CeedInt j = 0; j < n; j++) fprintf(stream, fp_fmt, fabs(a[i * n + j]) > 1E-14 ? a[i * n + j] : 0); 174 fputs("\n", stream); 175 } 176 return CEED_ERROR_SUCCESS; 177 } 178 179 /** 180 @brief Create the interpolation and gradient matrices for projection from the nodes of `basis_from` to the nodes of `basis_to`. 181 182 The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization. 183 The gradient is given by `grad_project = interp_to^+ * grad_from`, and is only computed for \f$H^1\f$ spaces otherwise it should not be used. 184 185 Note: `basis_from` and `basis_to` must have compatible quadrature spaces. 186 187 @param[in] basis_from `CeedBasis` to project from 188 @param[in] basis_to `CeedBasis` to project to 189 @param[out] interp_project Address of the variable where the newly created interpolation matrix will be stored 190 @param[out] grad_project Address of the variable where the newly created gradient matrix will be stored 191 192 @return An error code: 0 - success, otherwise - failure 193 194 @ref Developer 195 **/ 196 static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis basis_to, CeedScalar **interp_project, CeedScalar **grad_project) { 197 bool are_both_tensor; 198 CeedInt Q, Q_to, Q_from, P_to, P_from; 199 200 // Check for compatible quadrature spaces 201 CeedCall(CeedBasisGetNumQuadraturePoints(basis_to, &Q_to)); 202 CeedCall(CeedBasisGetNumQuadraturePoints(basis_from, &Q_from)); 203 CeedCheck(Q_to == Q_from, CeedBasisReturnCeed(basis_to), CEED_ERROR_DIMENSION, 204 "Bases must have compatible quadrature spaces." 205 " 'basis_from' has %" CeedInt_FMT " points and 'basis_to' has %" CeedInt_FMT, 206 Q_from, Q_to); 207 Q = Q_to; 208 209 // Check for matching tensor or non-tensor 210 { 211 bool is_tensor_to, is_tensor_from; 212 213 CeedCall(CeedBasisIsTensor(basis_to, &is_tensor_to)); 214 CeedCall(CeedBasisIsTensor(basis_from, &is_tensor_from)); 215 are_both_tensor = is_tensor_to && is_tensor_from; 216 } 217 if (are_both_tensor) { 218 CeedCall(CeedBasisGetNumNodes1D(basis_to, &P_to)); 219 CeedCall(CeedBasisGetNumNodes1D(basis_from, &P_from)); 220 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis_from, &Q)); 221 } else { 222 CeedCall(CeedBasisGetNumNodes(basis_to, &P_to)); 223 CeedCall(CeedBasisGetNumNodes(basis_from, &P_from)); 224 } 225 226 // Check for matching FE space 227 CeedFESpace fe_space_to, fe_space_from; 228 229 CeedCall(CeedBasisGetFESpace(basis_to, &fe_space_to)); 230 CeedCall(CeedBasisGetFESpace(basis_from, &fe_space_from)); 231 CeedCheck(fe_space_to == fe_space_from, CeedBasisReturnCeed(basis_to), CEED_ERROR_MINOR, 232 "Bases must both be the same FE space type." 233 " 'basis_from' is a %s and 'basis_to' is a %s", 234 CeedFESpaces[fe_space_from], CeedFESpaces[fe_space_to]); 235 236 // Get source matrices 237 CeedInt dim, q_comp = 1; 238 CeedScalar *interp_to_inv, *interp_from; 239 const CeedScalar *interp_to_source = NULL, *interp_from_source = NULL, *grad_from_source = NULL; 240 241 CeedCall(CeedBasisGetDimension(basis_from, &dim)); 242 if (are_both_tensor) { 243 CeedCall(CeedBasisGetInterp1D(basis_to, &interp_to_source)); 244 CeedCall(CeedBasisGetInterp1D(basis_from, &interp_from_source)); 245 } else { 246 CeedCall(CeedBasisGetNumQuadratureComponents(basis_from, CEED_EVAL_INTERP, &q_comp)); 247 CeedCall(CeedBasisGetInterp(basis_to, &interp_to_source)); 248 CeedCall(CeedBasisGetInterp(basis_from, &interp_from_source)); 249 } 250 CeedCall(CeedMalloc(Q * P_from * q_comp, &interp_from)); 251 CeedCall(CeedCalloc(P_to * P_from, interp_project)); 252 253 // `grad_project = interp_to^+ * grad_from` is computed for the H^1 space case so the 254 // projection basis will have a gradient operation (allocated even if not H^1 for the 255 // basis construction later on) 256 if (fe_space_to == CEED_FE_SPACE_H1) { 257 if (are_both_tensor) { 258 CeedCall(CeedBasisGetGrad1D(basis_from, &grad_from_source)); 259 } else { 260 CeedCall(CeedBasisGetGrad(basis_from, &grad_from_source)); 261 } 262 } 263 CeedCall(CeedCalloc(P_to * P_from * (are_both_tensor ? 1 : dim), grad_project)); 264 265 // Compute interp_to^+, pseudoinverse of interp_to 266 CeedCall(CeedCalloc(Q * q_comp * P_to, &interp_to_inv)); 267 CeedCall(CeedMatrixPseudoinverse(CeedBasisReturnCeed(basis_to), interp_to_source, Q * q_comp, P_to, interp_to_inv)); 268 // Build matrices 269 CeedInt num_matrices = 1 + (fe_space_to == CEED_FE_SPACE_H1) * (are_both_tensor ? 1 : dim); 270 CeedScalar *input_from[num_matrices], *output_project[num_matrices]; 271 272 input_from[0] = (CeedScalar *)interp_from_source; 273 output_project[0] = *interp_project; 274 for (CeedInt m = 1; m < num_matrices; m++) { 275 input_from[m] = (CeedScalar *)&grad_from_source[(m - 1) * Q * P_from]; 276 output_project[m] = &((*grad_project)[(m - 1) * P_to * P_from]); 277 } 278 for (CeedInt m = 0; m < num_matrices; m++) { 279 // output_project = interp_to^+ * interp_from 280 memcpy(interp_from, input_from[m], Q * P_from * q_comp * sizeof(input_from[m][0])); 281 CeedCall(CeedMatrixMatrixMultiply(CeedBasisReturnCeed(basis_to), interp_to_inv, input_from[m], output_project[m], P_to, P_from, Q * q_comp)); 282 // Round zero to machine precision 283 for (CeedInt i = 0; i < P_to * P_from; i++) { 284 if (fabs(output_project[m][i]) < 10 * CEED_EPSILON) output_project[m][i] = 0.0; 285 } 286 } 287 288 // Cleanup 289 CeedCall(CeedFree(&interp_to_inv)); 290 CeedCall(CeedFree(&interp_from)); 291 return CEED_ERROR_SUCCESS; 292 } 293 294 /** 295 @brief Check input vector dimensions for CeedBasisApply[Add] 296 297 @param[in] basis `CeedBasis` to evaluate 298 @param[in] num_elem The number of elements to apply the basis evaluation to; 299 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 300 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points; 301 @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes 302 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, 303 @ref CEED_EVAL_INTERP to use interpolated values, 304 @ref CEED_EVAL_GRAD to use gradients, 305 @ref CEED_EVAL_DIV to use divergence, 306 @ref CEED_EVAL_CURL to use curl, 307 @ref CEED_EVAL_WEIGHT to use quadrature weights 308 @param[in] u Input `CeedVector` 309 @param[out] v Output `CeedVector` 310 311 @return An error code: 0 - success, otherwise - failure 312 313 @ref Developer 314 **/ 315 static int CeedBasisApplyCheckDims(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 316 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; 317 CeedSize u_length = 0, v_length; 318 319 CeedCall(CeedBasisGetDimension(basis, &dim)); 320 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 321 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 322 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 323 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 324 CeedCall(CeedVectorGetLength(v, &v_length)); 325 if (u) CeedCall(CeedVectorGetLength(u, &u_length)); 326 327 // Check vector lengths to prevent out of bounds issues 328 bool has_good_dims = true; 329 switch (eval_mode) { 330 case CEED_EVAL_NONE: 331 case CEED_EVAL_INTERP: 332 case CEED_EVAL_GRAD: 333 case CEED_EVAL_DIV: 334 case CEED_EVAL_CURL: 335 has_good_dims = ((t_mode == CEED_TRANSPOSE && u_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_qpts * (CeedSize)q_comp && 336 v_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes) || 337 (t_mode == CEED_NOTRANSPOSE && v_length >= (CeedSize)num_elem * (CeedSize)num_qpts * (CeedSize)num_comp * (CeedSize)q_comp && 338 u_length >= (CeedSize)num_elem * (CeedSize)num_comp * (CeedSize)num_nodes)); 339 break; 340 case CEED_EVAL_WEIGHT: 341 has_good_dims = v_length >= (CeedSize)num_elem * (CeedSize)num_qpts; 342 break; 343 } 344 CeedCheck(has_good_dims, CeedBasisReturnCeed(basis), CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); 345 return CEED_ERROR_SUCCESS; 346 } 347 348 /** 349 @brief Check input vector dimensions for CeedBasisApply[Add]AtPoints 350 351 @param[in] basis `CeedBasis` to evaluate 352 @param[in] num_elem The number of elements to apply the basis evaluation to; 353 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 354 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 355 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 356 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes 357 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 358 @ref CEED_EVAL_GRAD to use gradients, 359 @ref CEED_EVAL_WEIGHT to use quadrature weights 360 @param[in] x_ref `CeedVector` holding reference coordinates of each point 361 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 362 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 363 364 @return An error code: 0 - success, otherwise - failure 365 366 @ref Developer 367 **/ 368 static int CeedBasisApplyAtPointsCheckDims(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 369 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 370 CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1, total_num_points = 0; 371 CeedSize x_length = 0, u_length = 0, v_length; 372 373 CeedCall(CeedBasisGetDimension(basis, &dim)); 374 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 375 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 376 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 377 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &num_q_comp)); 378 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 379 CeedCall(CeedVectorGetLength(v, &v_length)); 380 if (x_ref != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(x_ref, &x_length)); 381 if (u != CEED_VECTOR_NONE) CeedCall(CeedVectorGetLength(u, &u_length)); 382 383 // Check compatibility coordinates vector 384 for (CeedInt i = 0; i < num_elem; i++) total_num_points += num_points[i]; 385 CeedCheck((x_length >= (CeedSize)total_num_points * (CeedSize)dim) || (eval_mode == CEED_EVAL_WEIGHT), CeedBasisReturnCeed(basis), 386 CEED_ERROR_DIMENSION, 387 "Length of reference coordinate vector incompatible with basis dimension and number of points." 388 " Found reference coordinate vector of length %" CeedSize_FMT ", not of length %" CeedSize_FMT ".", 389 x_length, (CeedSize)total_num_points * (CeedSize)dim); 390 391 // Check CEED_EVAL_WEIGHT only on CEED_NOTRANSPOSE 392 CeedCheck(eval_mode != CEED_EVAL_WEIGHT || t_mode == CEED_NOTRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 393 "CEED_EVAL_WEIGHT only supported with CEED_NOTRANSPOSE"); 394 395 // Check vector lengths to prevent out of bounds issues 396 bool has_good_dims = true; 397 switch (eval_mode) { 398 case CEED_EVAL_INTERP: 399 has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp || 400 v_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)) || 401 (t_mode == CEED_NOTRANSPOSE && (v_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp || 402 u_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp))); 403 break; 404 case CEED_EVAL_GRAD: 405 has_good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp * (CeedSize)dim || 406 v_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp)) || 407 (t_mode == CEED_NOTRANSPOSE && (v_length >= (CeedSize)total_num_points * (CeedSize)num_q_comp * (CeedSize)dim || 408 u_length >= (CeedSize)num_elem * (CeedSize)num_nodes * (CeedSize)num_comp))); 409 break; 410 case CEED_EVAL_WEIGHT: 411 has_good_dims = t_mode == CEED_NOTRANSPOSE && (v_length >= total_num_points); 412 break; 413 // LCOV_EXCL_START 414 case CEED_EVAL_NONE: 415 case CEED_EVAL_DIV: 416 case CEED_EVAL_CURL: 417 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points not supported for %s", 418 CeedEvalModes[eval_mode]); 419 // LCOV_EXCL_STOP 420 } 421 CeedCheck(has_good_dims, CeedBasisReturnCeed(basis), CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode"); 422 return CEED_ERROR_SUCCESS; 423 } 424 425 /** 426 @brief Default implimentation to apply basis evaluation from nodes to arbitrary points 427 428 @param[in] basis `CeedBasis` to evaluate 429 @param[in] apply_add Sum result into target vector or overwrite 430 @param[in] num_elem The number of elements to apply the basis evaluation to; 431 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 432 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 433 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 434 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes 435 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 436 @ref CEED_EVAL_GRAD to use gradients, 437 @ref CEED_EVAL_WEIGHT to use quadrature weights 438 @param[in] x_ref `CeedVector` holding reference coordinates of each point 439 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 440 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 441 442 @return An error code: 0 - success, otherwise - failure 443 444 @ref Developer 445 **/ 446 static int CeedBasisApplyAtPoints_Core(CeedBasis basis, bool apply_add, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, 447 CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u, CeedVector v) { 448 CeedInt dim, num_comp, P_1d = 1, Q_1d = 1, total_num_points = num_points[0]; 449 450 CeedCall(CeedBasisGetDimension(basis, &dim)); 451 // Inserting check because clang-tidy doesn't understand this cannot occur 452 CeedCheck(dim > 0, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Malformed CeedBasis, dim > 0 is required"); 453 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 454 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 455 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 456 457 // Default implementation 458 { 459 bool is_tensor_basis; 460 461 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 462 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 463 "Evaluation at arbitrary points only supported for tensor product bases"); 464 } 465 CeedCheck(num_elem == 1, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 466 "Evaluation at arbitrary points only supported for a single element at a time"); 467 if (eval_mode == CEED_EVAL_WEIGHT) { 468 CeedCall(CeedVectorSetValue(v, 1.0)); 469 return CEED_ERROR_SUCCESS; 470 } 471 if (!basis->basis_chebyshev) { 472 // Build basis mapping from nodes to Chebyshev coefficients 473 CeedScalar *chebyshev_interp_1d, *chebyshev_grad_1d, *chebyshev_q_weight_1d; 474 const CeedScalar *q_ref_1d; 475 Ceed ceed; 476 477 CeedCall(CeedCalloc(P_1d * Q_1d, &chebyshev_interp_1d)); 478 CeedCall(CeedCalloc(P_1d * Q_1d, &chebyshev_grad_1d)); 479 CeedCall(CeedCalloc(Q_1d, &chebyshev_q_weight_1d)); 480 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); 481 CeedCall(CeedBasisGetChebyshevInterp1D(basis, chebyshev_interp_1d)); 482 483 CeedCall(CeedBasisGetCeed(basis, &ceed)); 484 CeedCall(CeedVectorCreate(ceed, num_comp * CeedIntPow(Q_1d, dim), &basis->vec_chebyshev)); 485 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d, Q_1d, chebyshev_interp_1d, chebyshev_grad_1d, q_ref_1d, chebyshev_q_weight_1d, 486 &basis->basis_chebyshev)); 487 488 // Cleanup 489 CeedCall(CeedFree(&chebyshev_interp_1d)); 490 CeedCall(CeedFree(&chebyshev_grad_1d)); 491 CeedCall(CeedFree(&chebyshev_q_weight_1d)); 492 CeedCall(CeedDestroy(&ceed)); 493 } 494 495 // Create TensorContract object if needed, such as a basis from the GPU backends 496 if (!basis->contract) { 497 Ceed ceed_ref; 498 CeedBasis basis_ref = NULL; 499 500 CeedCall(CeedInit("/cpu/self", &ceed_ref)); 501 // Only need matching tensor contraction dimensions, any type of basis will work 502 CeedCall(CeedBasisCreateTensorH1Lagrange(ceed_ref, dim, num_comp, P_1d, Q_1d, CEED_GAUSS, &basis_ref)); 503 // Note - clang-tidy doesn't know basis_ref->contract must be valid here 504 CeedCheck(basis_ref && basis_ref->contract, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, 505 "Reference CPU ceed failed to create a tensor contraction object"); 506 CeedCall(CeedTensorContractReferenceCopy(basis_ref->contract, &basis->contract)); 507 CeedCall(CeedBasisDestroy(&basis_ref)); 508 CeedCall(CeedDestroy(&ceed_ref)); 509 } 510 511 // Basis evaluation 512 switch (t_mode) { 513 case CEED_NOTRANSPOSE: { 514 // Nodes to arbitrary points 515 CeedScalar *v_array; 516 const CeedScalar *chebyshev_coeffs, *x_array_read; 517 518 // -- Interpolate to Chebyshev coefficients 519 CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, basis->vec_chebyshev)); 520 521 // -- Evaluate Chebyshev polynomials at arbitrary points 522 CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); 523 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read)); 524 CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array)); 525 switch (eval_mode) { 526 case CEED_EVAL_INTERP: { 527 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 528 529 // ---- Values at point 530 for (CeedInt p = 0; p < total_num_points; p++) { 531 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; 532 533 for (CeedInt d = 0; d < dim; d++) { 534 // ------ Tensor contract with current Chebyshev polynomial values 535 CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 536 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, 537 d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2])); 538 pre /= Q_1d; 539 post *= 1; 540 } 541 for (CeedInt c = 0; c < num_comp; c++) v_array[c * total_num_points + p] = tmp[dim % 2][c]; 542 } 543 break; 544 } 545 case CEED_EVAL_GRAD: { 546 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 547 548 // ---- Values at point 549 for (CeedInt p = 0; p < total_num_points; p++) { 550 // Dim**2 contractions, apply grad when pass == dim 551 for (CeedInt pass = 0; pass < dim; pass++) { 552 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1; 553 554 for (CeedInt d = 0; d < dim; d++) { 555 // ------ Tensor contract with current Chebyshev polynomial values 556 if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 557 else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 558 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false, 559 d == 0 ? chebyshev_coeffs : tmp[d % 2], tmp[(d + 1) % 2])); 560 pre /= Q_1d; 561 post *= 1; 562 } 563 for (CeedInt c = 0; c < num_comp; c++) v_array[(pass * num_comp + c) * total_num_points + p] = tmp[dim % 2][c]; 564 } 565 } 566 break; 567 } 568 default: 569 // Nothing to do, excluded above 570 break; 571 } 572 CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs)); 573 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read)); 574 CeedCall(CeedVectorRestoreArray(v, &v_array)); 575 break; 576 } 577 case CEED_TRANSPOSE: { 578 // Note: No switch on e_mode here because only CEED_EVAL_INTERP is supported at this time 579 // Arbitrary points to nodes 580 CeedScalar *chebyshev_coeffs; 581 const CeedScalar *u_array, *x_array_read; 582 583 // -- Transpose of evaluation of Chebyshev polynomials at arbitrary points 584 CeedCall(CeedVectorGetArrayWrite(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs)); 585 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read)); 586 CeedCall(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array)); 587 588 switch (eval_mode) { 589 case CEED_EVAL_INTERP: { 590 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 591 592 // ---- Values at point 593 for (CeedInt p = 0; p < total_num_points; p++) { 594 CeedInt pre = num_comp * 1, post = 1; 595 596 for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[c * total_num_points + p]; 597 for (CeedInt d = 0; d < dim; d++) { 598 // ------ Tensor contract with current Chebyshev polynomial values 599 CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 600 CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == (dim - 1), tmp[d % 2], 601 d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); 602 pre /= 1; 603 post *= Q_1d; 604 } 605 } 606 break; 607 } 608 case CEED_EVAL_GRAD: { 609 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d]; 610 611 // ---- Values at point 612 for (CeedInt p = 0; p < total_num_points; p++) { 613 // Dim**2 contractions, apply grad when pass == dim 614 for (CeedInt pass = 0; pass < dim; pass++) { 615 CeedInt pre = num_comp * 1, post = 1; 616 617 for (CeedInt c = 0; c < num_comp; c++) tmp[0][c] = u_array[(pass * num_comp + c) * total_num_points + p]; 618 for (CeedInt d = 0; d < dim; d++) { 619 // ------ Tensor contract with current Chebyshev polynomial values 620 if (pass == d) CeedCall(CeedChebyshevDerivativeAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 621 else CeedCall(CeedChebyshevPolynomialsAtPoint(x_array_read[d * total_num_points + p], Q_1d, chebyshev_x)); 622 CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, 623 (p > 0 || (p == 0 && pass > 0)) && d == (dim - 1), tmp[d % 2], 624 d == (dim - 1) ? chebyshev_coeffs : tmp[(d + 1) % 2])); 625 pre /= 1; 626 post *= Q_1d; 627 } 628 } 629 } 630 break; 631 } 632 default: 633 // Nothing to do, excluded above 634 break; 635 } 636 CeedCall(CeedVectorRestoreArray(basis->vec_chebyshev, &chebyshev_coeffs)); 637 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read)); 638 CeedCall(CeedVectorRestoreArrayRead(u, &u_array)); 639 640 // -- Interpolate transpose from Chebyshev coefficients 641 if (apply_add) CeedCall(CeedBasisApplyAdd(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); 642 else CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v)); 643 break; 644 } 645 } 646 return CEED_ERROR_SUCCESS; 647 } 648 649 /// @} 650 651 /// ---------------------------------------------------------------------------- 652 /// Ceed Backend API 653 /// ---------------------------------------------------------------------------- 654 /// @addtogroup CeedBasisBackend 655 /// @{ 656 657 /** 658 @brief Fallback to a reference implementation for a non tensor-product basis for \f$H^1\f$ discretizations. 659 This function may only be called inside of a backend `BasisCreateH1` function. 660 This is used by a backend when the specific parameters for a `CeedBasis` exceed the backend's support, such as 661 when a `interp` and `grad` matrices require too many bytes to fit into shared memory on a GPU. 662 663 @param[in] ceed `Ceed` object used to create the `CeedBasis` 664 @param[in] topo Topology of element, e.g. hypercube, simplex, etc 665 @param[in] num_comp Number of field components (1 for scalar fields) 666 @param[in] num_nodes Total number of nodes 667 @param[in] num_qpts Total number of quadrature points 668 @param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal basis functions at quadrature points 669 @param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of nodal basis functions at quadrature points 670 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element 671 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 672 @param[out] basis Newly created `CeedBasis` 673 674 @return An error code: 0 - success, otherwise - failure 675 676 @ref User 677 **/ 678 int CeedBasisCreateH1Fallback(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 679 const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) { 680 CeedInt P = num_nodes, Q = num_qpts, dim = 0; 681 Ceed delegate; 682 683 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 684 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateH1"); 685 686 CeedCall(CeedReferenceCopy(delegate, &(basis)->ceed)); 687 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 688 CeedCall(delegate->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, basis)); 689 CeedCall(CeedDestroy(&delegate)); 690 return CEED_ERROR_SUCCESS; 691 } 692 693 /** 694 @brief Return collocated gradient matrix 695 696 @param[in] basis `CeedBasis` 697 @param[out] collo_grad_1d Row-major (`Q_1d * Q_1d`) matrix expressing derivatives of basis functions at quadrature points 698 699 @return An error code: 0 - success, otherwise - failure 700 701 @ref Backend 702 **/ 703 int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { 704 Ceed ceed; 705 CeedInt P_1d, Q_1d; 706 CeedScalar *interp_1d_pinv; 707 const CeedScalar *grad_1d, *interp_1d; 708 709 // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure. 710 CeedCall(CeedBasisGetCeed(basis, &ceed)); 711 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 712 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 713 714 // Compute interp_1d^+, pseudoinverse of interp_1d 715 CeedCall(CeedCalloc(P_1d * Q_1d, &interp_1d_pinv)); 716 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 717 CeedCall(CeedMatrixPseudoinverse(ceed, interp_1d, Q_1d, P_1d, interp_1d_pinv)); 718 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 719 CeedCall(CeedMatrixMatrixMultiply(ceed, grad_1d, (const CeedScalar *)interp_1d_pinv, collo_grad_1d, Q_1d, Q_1d, P_1d)); 720 721 CeedCall(CeedFree(&interp_1d_pinv)); 722 CeedCall(CeedDestroy(&ceed)); 723 return CEED_ERROR_SUCCESS; 724 } 725 726 /** 727 @brief Return 1D interpolation matrix to Chebyshev polynomial coefficients on quadrature space 728 729 @param[in] basis `CeedBasis` 730 @param[out] chebyshev_interp_1d Row-major (`P_1d * Q_1d`) matrix interpolating from basis nodes to Chebyshev polynomial coefficients 731 732 @return An error code: 0 - success, otherwise - failure 733 734 @ref Backend 735 **/ 736 int CeedBasisGetChebyshevInterp1D(CeedBasis basis, CeedScalar *chebyshev_interp_1d) { 737 CeedInt P_1d, Q_1d; 738 CeedScalar *C, *chebyshev_coeffs_1d_inv; 739 const CeedScalar *interp_1d, *q_ref_1d; 740 Ceed ceed; 741 742 CeedCall(CeedBasisGetCeed(basis, &ceed)); 743 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 744 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 745 746 // Build coefficient matrix 747 // -- Note: Clang-tidy needs this check 748 CeedCheck((P_1d > 0) && (Q_1d > 0), ceed, CEED_ERROR_INCOMPATIBLE, "CeedBasis dimensions are malformed"); 749 CeedCall(CeedCalloc(Q_1d * Q_1d, &C)); 750 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); 751 for (CeedInt i = 0; i < Q_1d; i++) CeedCall(CeedChebyshevPolynomialsAtPoint(q_ref_1d[i], Q_1d, &C[i * Q_1d])); 752 753 // Compute C^+, pseudoinverse of coefficient matrix 754 CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_coeffs_1d_inv)); 755 CeedCall(CeedMatrixPseudoinverse(ceed, C, Q_1d, Q_1d, chebyshev_coeffs_1d_inv)); 756 757 // Build mapping from nodes to Chebyshev coefficients 758 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 759 CeedCall(CeedMatrixMatrixMultiply(ceed, chebyshev_coeffs_1d_inv, interp_1d, chebyshev_interp_1d, Q_1d, P_1d, Q_1d)); 760 761 // Cleanup 762 CeedCall(CeedFree(&C)); 763 CeedCall(CeedFree(&chebyshev_coeffs_1d_inv)); 764 CeedCall(CeedDestroy(&ceed)); 765 return CEED_ERROR_SUCCESS; 766 } 767 768 /** 769 @brief Get tensor status for given `CeedBasis` 770 771 @param[in] basis `CeedBasis` 772 @param[out] is_tensor Variable to store tensor status 773 774 @return An error code: 0 - success, otherwise - failure 775 776 @ref Backend 777 **/ 778 int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor) { 779 *is_tensor = basis->is_tensor_basis; 780 return CEED_ERROR_SUCCESS; 781 } 782 783 /** 784 @brief Get backend data of a `CeedBasis` 785 786 @param[in] basis `CeedBasis` 787 @param[out] data Variable to store data 788 789 @return An error code: 0 - success, otherwise - failure 790 791 @ref Backend 792 **/ 793 int CeedBasisGetData(CeedBasis basis, void *data) { 794 *(void **)data = basis->data; 795 return CEED_ERROR_SUCCESS; 796 } 797 798 /** 799 @brief Set backend data of a `CeedBasis` 800 801 @param[in,out] basis `CeedBasis` 802 @param[in] data Data to set 803 804 @return An error code: 0 - success, otherwise - failure 805 806 @ref Backend 807 **/ 808 int CeedBasisSetData(CeedBasis basis, void *data) { 809 basis->data = data; 810 return CEED_ERROR_SUCCESS; 811 } 812 813 /** 814 @brief Increment the reference counter for a `CeedBasis` 815 816 @param[in,out] basis `CeedBasis` to increment the reference counter 817 818 @return An error code: 0 - success, otherwise - failure 819 820 @ref Backend 821 **/ 822 int CeedBasisReference(CeedBasis basis) { 823 basis->ref_count++; 824 return CEED_ERROR_SUCCESS; 825 } 826 827 /** 828 @brief Get number of Q-vector components for given `CeedBasis` 829 830 @param[in] basis `CeedBasis` 831 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 832 @ref CEED_EVAL_GRAD to use gradients, 833 @ref CEED_EVAL_DIV to use divergence, 834 @ref CEED_EVAL_CURL to use curl 835 @param[out] q_comp Variable to store number of Q-vector components of basis 836 837 @return An error code: 0 - success, otherwise - failure 838 839 @ref Backend 840 **/ 841 int CeedBasisGetNumQuadratureComponents(CeedBasis basis, CeedEvalMode eval_mode, CeedInt *q_comp) { 842 CeedInt dim; 843 844 CeedCall(CeedBasisGetDimension(basis, &dim)); 845 switch (eval_mode) { 846 case CEED_EVAL_INTERP: { 847 CeedFESpace fe_space; 848 849 CeedCall(CeedBasisGetFESpace(basis, &fe_space)); 850 *q_comp = (fe_space == CEED_FE_SPACE_H1) ? 1 : dim; 851 } break; 852 case CEED_EVAL_GRAD: 853 *q_comp = dim; 854 break; 855 case CEED_EVAL_DIV: 856 *q_comp = 1; 857 break; 858 case CEED_EVAL_CURL: 859 *q_comp = (dim < 3) ? 1 : dim; 860 break; 861 case CEED_EVAL_NONE: 862 case CEED_EVAL_WEIGHT: 863 *q_comp = 1; 864 break; 865 } 866 return CEED_ERROR_SUCCESS; 867 } 868 869 /** 870 @brief Estimate number of FLOPs required to apply `CeedBasis` in `t_mode` and `eval_mode` 871 872 @param[in] basis `CeedBasis` to estimate FLOPs for 873 @param[in] t_mode Apply basis or transpose 874 @param[in] eval_mode @ref CeedEvalMode 875 @param[in] is_at_points Evaluate the basis at points or quadrature points 876 @param[in] num_points Number of points basis is evaluated at 877 @param[out] flops Address of variable to hold FLOPs estimate 878 879 @ref Backend 880 **/ 881 int CeedBasisGetFlopsEstimate(CeedBasis basis, CeedTransposeMode t_mode, CeedEvalMode eval_mode, bool is_at_points, CeedInt num_points, 882 CeedSize *flops) { 883 bool is_tensor; 884 885 CeedCall(CeedBasisIsTensor(basis, &is_tensor)); 886 CeedCheck(!is_at_points || is_tensor, CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Can only evaluate tensor-product bases at points"); 887 if (is_tensor) { 888 CeedInt dim, num_comp, P_1d, Q_1d; 889 890 CeedCall(CeedBasisGetDimension(basis, &dim)); 891 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 892 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 893 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 894 if (t_mode == CEED_TRANSPOSE) { 895 P_1d = Q_1d; 896 Q_1d = P_1d; 897 } 898 CeedInt tensor_flops = 0, pre = num_comp * CeedIntPow(P_1d, dim - 1), post = 1; 899 900 for (CeedInt d = 0; d < dim; d++) { 901 tensor_flops += 2 * pre * P_1d * post * Q_1d; 902 pre /= P_1d; 903 post *= Q_1d; 904 } 905 if (is_at_points) { 906 bool is_gpu = false; 907 908 { 909 CeedMemType mem_type; 910 911 CeedCall(CeedGetPreferredMemType(CeedBasisReturnCeed(basis), &mem_type)); 912 is_gpu = mem_type == CEED_MEM_DEVICE; 913 } 914 915 CeedInt chebyshev_flops = (Q_1d - 2) * 3 + 1, d_chebyshev_flops = (Q_1d - 2) * 8 + 1; 916 CeedInt point_tensor_flops = 0, pre = CeedIntPow(Q_1d, dim - 1), post = 1; 917 918 for (CeedInt d = 0; d < dim; d++) { 919 point_tensor_flops += 2 * pre * Q_1d * post * 1; 920 pre /= P_1d; 921 post *= Q_1d; 922 } 923 924 switch (eval_mode) { 925 case CEED_EVAL_NONE: 926 *flops = 0; 927 break; 928 case CEED_EVAL_INTERP: { 929 *flops = tensor_flops + num_points * num_comp * (point_tensor_flops + (t_mode == CEED_TRANSPOSE ? CeedIntPow(Q_1d, dim) : 0)); 930 if (dim == 3 && is_gpu) { 931 *flops += num_points * Q_1d * 932 (chebyshev_flops + num_comp * (2 * chebyshev_flops + 2 * Q_1d * Q_1d + (t_mode == CEED_TRANSPOSE ? 2 * Q_1d + 1 : 3 * Q_1d))); 933 } else { 934 *flops += num_points * (is_gpu ? num_comp : 1) * dim * chebyshev_flops; 935 } 936 break; 937 } 938 case CEED_EVAL_GRAD: { 939 *flops = tensor_flops + num_points * num_comp * (point_tensor_flops + (t_mode == CEED_TRANSPOSE ? CeedIntPow(Q_1d, dim) : 0)); 940 if (dim == 3 && is_gpu) { 941 CeedInt inner_flops = 942 dim * (2 * Q_1d * Q_1d + (t_mode == CEED_TRANSPOSE ? 2 : 3) * Q_1d) + (dim - 1) * (2 * chebyshev_flops + d_chebyshev_flops); 943 944 *flops += num_points * Q_1d * (chebyshev_flops + d_chebyshev_flops + num_comp * (inner_flops + (t_mode == CEED_TRANSPOSE ? 1 : 0))); 945 } else { 946 *flops += num_points * (is_gpu ? num_comp : 1) * dim * (d_chebyshev_flops + (dim - 1) * chebyshev_flops); 947 } 948 break; 949 } 950 case CEED_EVAL_DIV: 951 case CEED_EVAL_CURL: { 952 // LCOV_EXCL_START 953 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Tensor basis evaluation for %s not supported at points", 954 CeedEvalModes[eval_mode]); 955 break; 956 // LCOV_EXCL_STOP 957 } 958 case CEED_EVAL_WEIGHT: 959 *flops = num_points; 960 break; 961 } 962 } else { 963 switch (eval_mode) { 964 case CEED_EVAL_NONE: 965 *flops = 0; 966 break; 967 case CEED_EVAL_INTERP: 968 *flops = tensor_flops; 969 break; 970 case CEED_EVAL_GRAD: 971 *flops = tensor_flops * 2; 972 break; 973 case CEED_EVAL_DIV: 974 case CEED_EVAL_CURL: { 975 // LCOV_EXCL_START 976 return CeedError(CeedBasisReturnCeed(basis), CEED_ERROR_INCOMPATIBLE, "Tensor basis evaluation for %s not supported", 977 CeedEvalModes[eval_mode]); 978 break; 979 // LCOV_EXCL_STOP 980 } 981 case CEED_EVAL_WEIGHT: 982 *flops = dim * CeedIntPow(Q_1d, dim); 983 break; 984 } 985 } 986 } else { 987 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts; 988 989 CeedCall(CeedBasisGetDimension(basis, &dim)); 990 CeedCall(CeedBasisGetNumComponents(basis, &num_comp)); 991 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp)); 992 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes)); 993 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts)); 994 switch (eval_mode) { 995 case CEED_EVAL_NONE: 996 *flops = 0; 997 break; 998 case CEED_EVAL_INTERP: 999 case CEED_EVAL_GRAD: 1000 case CEED_EVAL_DIV: 1001 case CEED_EVAL_CURL: 1002 *flops = num_nodes * num_qpts * num_comp * q_comp; 1003 break; 1004 case CEED_EVAL_WEIGHT: 1005 *flops = 0; 1006 break; 1007 } 1008 } 1009 return CEED_ERROR_SUCCESS; 1010 } 1011 1012 /** 1013 @brief Get `CeedFESpace` for a `CeedBasis` 1014 1015 @param[in] basis `CeedBasis` 1016 @param[out] fe_space Variable to store `CeedFESpace` 1017 1018 @return An error code: 0 - success, otherwise - failure 1019 1020 @ref Backend 1021 **/ 1022 int CeedBasisGetFESpace(CeedBasis basis, CeedFESpace *fe_space) { 1023 *fe_space = basis->fe_space; 1024 return CEED_ERROR_SUCCESS; 1025 } 1026 1027 /** 1028 @brief Get dimension for given `CeedElemTopology` 1029 1030 @param[in] topo `CeedElemTopology` 1031 @param[out] dim Variable to store dimension of topology 1032 1033 @return An error code: 0 - success, otherwise - failure 1034 1035 @ref Backend 1036 **/ 1037 int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) { 1038 *dim = (CeedInt)topo >> 16; 1039 return CEED_ERROR_SUCCESS; 1040 } 1041 1042 /** 1043 @brief Get `CeedTensorContract` of a `CeedBasis` 1044 1045 @param[in] basis `CeedBasis` 1046 @param[out] contract Variable to store `CeedTensorContract` 1047 1048 @return An error code: 0 - success, otherwise - failure 1049 1050 @ref Backend 1051 **/ 1052 int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) { 1053 *contract = basis->contract; 1054 return CEED_ERROR_SUCCESS; 1055 } 1056 1057 /** 1058 @brief Set `CeedTensorContract` of a `CeedBasis` 1059 1060 @param[in,out] basis `CeedBasis` 1061 @param[in] contract `CeedTensorContract` to set 1062 1063 @return An error code: 0 - success, otherwise - failure 1064 1065 @ref Backend 1066 **/ 1067 int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract) { 1068 basis->contract = contract; 1069 CeedCall(CeedTensorContractReference(contract)); 1070 return CEED_ERROR_SUCCESS; 1071 } 1072 1073 /** 1074 @brief Return a reference implementation of matrix multiplication \f$C = A B\f$. 1075 1076 Note: This is a reference implementation for CPU `CeedScalar` pointers that is not intended for high performance. 1077 1078 @param[in] ceed `Ceed` context for error handling 1079 @param[in] mat_A Row-major matrix `A` 1080 @param[in] mat_B Row-major matrix `B` 1081 @param[out] mat_C Row-major output matrix `C` 1082 @param[in] m Number of rows of `C` 1083 @param[in] n Number of columns of `C` 1084 @param[in] kk Number of columns of `A`/rows of `B` 1085 1086 @return An error code: 0 - success, otherwise - failure 1087 1088 @ref Utility 1089 **/ 1090 int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n, CeedInt kk) { 1091 for (CeedInt i = 0; i < m; i++) { 1092 for (CeedInt j = 0; j < n; j++) { 1093 CeedScalar sum = 0; 1094 1095 for (CeedInt k = 0; k < kk; k++) sum += mat_A[k + i * kk] * mat_B[j + k * n]; 1096 mat_C[j + i * n] = sum; 1097 } 1098 } 1099 return CEED_ERROR_SUCCESS; 1100 } 1101 1102 /** 1103 @brief Return QR Factorization of a matrix 1104 1105 @param[in] ceed `Ceed` context for error handling 1106 @param[in,out] mat Row-major matrix to be factorized in place 1107 @param[in,out] tau Vector of length `m` of scaling factors 1108 @param[in] m Number of rows 1109 @param[in] n Number of columns 1110 1111 @return An error code: 0 - success, otherwise - failure 1112 1113 @ref Utility 1114 **/ 1115 int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n) { 1116 CeedScalar v[m]; 1117 1118 // Check matrix shape 1119 CeedCheck(n <= m, ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute QR factorization with n > m"); 1120 1121 for (CeedInt i = 0; i < n; i++) { 1122 CeedScalar sigma = 0.0; 1123 1124 if (i >= m - 1) { // last row of matrix, no reflection needed 1125 tau[i] = 0.; 1126 break; 1127 } 1128 // Calculate Householder vector, magnitude 1129 v[i] = mat[i + n * i]; 1130 for (CeedInt j = i + 1; j < m; j++) { 1131 v[j] = mat[i + n * j]; 1132 sigma += v[j] * v[j]; 1133 } 1134 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:m] 1135 const CeedScalar R_ii = -copysign(norm, v[i]); 1136 1137 v[i] -= R_ii; 1138 // norm of v[i:m] after modification above and scaling below 1139 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1140 // tau = 2 / (norm*norm) 1141 tau[i] = 2 * v[i] * v[i] / (v[i] * v[i] + sigma); 1142 for (CeedInt j = i + 1; j < m; j++) v[j] /= v[i]; 1143 1144 // Apply Householder reflector to lower right panel 1145 CeedHouseholderReflect(&mat[i * n + i + 1], &v[i], tau[i], m - i, n - i - 1, n, 1); 1146 // Save v 1147 mat[i + n * i] = R_ii; 1148 for (CeedInt j = i + 1; j < m; j++) mat[i + n * j] = v[j]; 1149 } 1150 return CEED_ERROR_SUCCESS; 1151 } 1152 1153 /** 1154 @brief Apply Householder Q matrix 1155 1156 Compute `mat_A = mat_Q mat_A`, where `mat_Q` is \f$m \times m\f$ and `mat_A` is \f$m \times n\f$. 1157 1158 @param[in,out] mat_A Matrix to apply Householder Q to, in place 1159 @param[in] mat_Q Householder Q matrix 1160 @param[in] tau Householder scaling factors 1161 @param[in] t_mode Transpose mode for application 1162 @param[in] m Number of rows in `A` 1163 @param[in] n Number of columns in `A` 1164 @param[in] k Number of elementary reflectors in Q, `k < m` 1165 @param[in] row Row stride in `A` 1166 @param[in] col Col stride in `A` 1167 1168 @return An error code: 0 - success, otherwise - failure 1169 1170 @ref Utility 1171 **/ 1172 int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n, 1173 CeedInt k, CeedInt row, CeedInt col) { 1174 CeedScalar *v; 1175 1176 CeedCall(CeedMalloc(m, &v)); 1177 for (CeedInt ii = 0; ii < k; ii++) { 1178 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k - 1 - ii; 1179 for (CeedInt j = i + 1; j < m; j++) v[j] = mat_Q[j * k + i]; 1180 // Apply Householder reflector (I - tau v v^T) collo_grad_1d^T 1181 CeedCall(CeedHouseholderReflect(&mat_A[i * row], &v[i], tau[i], m - i, n, row, col)); 1182 } 1183 CeedCall(CeedFree(&v)); 1184 return CEED_ERROR_SUCCESS; 1185 } 1186 1187 /** 1188 @brief Return pseudoinverse of a matrix 1189 1190 @param[in] ceed Ceed context for error handling 1191 @param[in] mat Row-major matrix to compute pseudoinverse of 1192 @param[in] m Number of rows 1193 @param[in] n Number of columns 1194 @param[out] mat_pinv Row-major pseudoinverse matrix 1195 1196 @return An error code: 0 - success, otherwise - failure 1197 1198 @ref Utility 1199 **/ 1200 int CeedMatrixPseudoinverse(Ceed ceed, const CeedScalar *mat, CeedInt m, CeedInt n, CeedScalar *mat_pinv) { 1201 CeedScalar *tau, *I, *mat_copy; 1202 1203 CeedCall(CeedCalloc(m, &tau)); 1204 CeedCall(CeedCalloc(m * m, &I)); 1205 CeedCall(CeedCalloc(m * n, &mat_copy)); 1206 memcpy(mat_copy, mat, m * n * sizeof mat[0]); 1207 1208 // QR Factorization, mat = Q R 1209 CeedCall(CeedQRFactorization(ceed, mat_copy, tau, m, n)); 1210 1211 // -- Apply Q^T, I = Q^T * I 1212 for (CeedInt i = 0; i < m; i++) I[i * m + i] = 1.0; 1213 CeedCall(CeedHouseholderApplyQ(I, mat_copy, tau, CEED_TRANSPOSE, m, m, n, m, 1)); 1214 // -- Apply R_inv, mat_pinv = R_inv * Q^T 1215 for (CeedInt j = 0; j < m; j++) { // Column j 1216 mat_pinv[j + m * (n - 1)] = I[j + m * (n - 1)] / mat_copy[n * n - 1]; 1217 for (CeedInt i = n - 2; i >= 0; i--) { // Row i 1218 mat_pinv[j + m * i] = I[j + m * i]; 1219 for (CeedInt k = i + 1; k < n; k++) mat_pinv[j + m * i] -= mat_copy[k + n * i] * mat_pinv[j + m * k]; 1220 mat_pinv[j + m * i] /= mat_copy[i + n * i]; 1221 } 1222 } 1223 1224 // Cleanup 1225 CeedCall(CeedFree(&I)); 1226 CeedCall(CeedFree(&tau)); 1227 CeedCall(CeedFree(&mat_copy)); 1228 return CEED_ERROR_SUCCESS; 1229 } 1230 1231 /** 1232 @brief Return symmetric Schur decomposition of the symmetric matrix mat via symmetric QR factorization 1233 1234 @param[in] ceed `Ceed` context for error handling 1235 @param[in,out] mat Row-major matrix to be factorized in place 1236 @param[out] lambda Vector of length n of eigenvalues 1237 @param[in] n Number of rows/columns 1238 1239 @return An error code: 0 - success, otherwise - failure 1240 1241 @ref Utility 1242 **/ 1243 CeedPragmaOptimizeOff 1244 int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n) { 1245 // Check bounds for clang-tidy 1246 CeedCheck(n > 1, ceed, CEED_ERROR_UNSUPPORTED, "Cannot compute symmetric Schur decomposition of scalars"); 1247 1248 CeedScalar v[n - 1], tau[n - 1], mat_T[n * n]; 1249 1250 // Copy mat to mat_T and set mat to I 1251 memcpy(mat_T, mat, n * n * sizeof(mat[0])); 1252 for (CeedInt i = 0; i < n; i++) { 1253 for (CeedInt j = 0; j < n; j++) mat[j + n * i] = (i == j) ? 1 : 0; 1254 } 1255 1256 // Reduce to tridiagonal 1257 for (CeedInt i = 0; i < n - 1; i++) { 1258 // Calculate Householder vector, magnitude 1259 CeedScalar sigma = 0.0; 1260 1261 v[i] = mat_T[i + n * (i + 1)]; 1262 for (CeedInt j = i + 1; j < n - 1; j++) { 1263 v[j] = mat_T[i + n * (j + 1)]; 1264 sigma += v[j] * v[j]; 1265 } 1266 const CeedScalar norm = sqrt(v[i] * v[i] + sigma); // norm of v[i:n-1] 1267 const CeedScalar R_ii = -copysign(norm, v[i]); 1268 1269 v[i] -= R_ii; 1270 // norm of v[i:m] after modification above and scaling below 1271 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1272 // tau = 2 / (norm*norm) 1273 tau[i] = i == n - 2 ? 2 : 2 * v[i] * v[i] / (v[i] * v[i] + sigma); 1274 for (CeedInt j = i + 1; j < n - 1; j++) v[j] /= v[i]; 1275 1276 // Update sub and super diagonal 1277 for (CeedInt j = i + 2; j < n; j++) { 1278 mat_T[i + n * j] = 0; 1279 mat_T[j + n * i] = 0; 1280 } 1281 // Apply symmetric Householder reflector to lower right panel 1282 CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1); 1283 CeedHouseholderReflect(&mat_T[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), 1, n); 1284 1285 // Save v 1286 mat_T[i + n * (i + 1)] = R_ii; 1287 mat_T[(i + 1) + n * i] = R_ii; 1288 for (CeedInt j = i + 1; j < n - 1; j++) { 1289 mat_T[i + n * (j + 1)] = v[j]; 1290 } 1291 } 1292 // Backwards accumulation of Q 1293 for (CeedInt i = n - 2; i >= 0; i--) { 1294 if (tau[i] > 0.0) { 1295 v[i] = 1; 1296 for (CeedInt j = i + 1; j < n - 1; j++) { 1297 v[j] = mat_T[i + n * (j + 1)]; 1298 mat_T[i + n * (j + 1)] = 0; 1299 } 1300 CeedHouseholderReflect(&mat[(i + 1) + n * (i + 1)], &v[i], tau[i], n - (i + 1), n - (i + 1), n, 1); 1301 } 1302 } 1303 1304 // Reduce sub and super diagonal 1305 CeedInt p = 0, q = 0, itr = 0, max_itr = n * n * n * n; 1306 CeedScalar tol = CEED_EPSILON; 1307 1308 while (itr < max_itr) { 1309 // Update p, q, size of reduced portions of diagonal 1310 p = 0; 1311 q = 0; 1312 for (CeedInt i = n - 2; i >= 0; i--) { 1313 if (fabs(mat_T[i + n * (i + 1)]) < tol) q += 1; 1314 else break; 1315 } 1316 for (CeedInt i = 0; i < n - q - 1; i++) { 1317 if (fabs(mat_T[i + n * (i + 1)]) < tol) p += 1; 1318 else break; 1319 } 1320 if (q == n - 1) break; // Finished reducing 1321 1322 // Reduce tridiagonal portion 1323 CeedScalar t_nn = mat_T[(n - 1 - q) + n * (n - 1 - q)], t_nnm1 = mat_T[(n - 2 - q) + n * (n - 1 - q)]; 1324 CeedScalar d = (mat_T[(n - 2 - q) + n * (n - 2 - q)] - t_nn) / 2; 1325 CeedScalar mu = t_nn - t_nnm1 * t_nnm1 / (d + copysign(sqrt(d * d + t_nnm1 * t_nnm1), d)); 1326 CeedScalar x = mat_T[p + n * p] - mu; 1327 CeedScalar z = mat_T[p + n * (p + 1)]; 1328 1329 for (CeedInt k = p; k < n - q - 1; k++) { 1330 // Compute Givens rotation 1331 CeedScalar c = 1, s = 0; 1332 1333 if (fabs(z) > tol) { 1334 if (fabs(z) > fabs(x)) { 1335 const CeedScalar tau = -x / z; 1336 1337 s = 1 / sqrt(1 + tau * tau); 1338 c = s * tau; 1339 } else { 1340 const CeedScalar tau = -z / x; 1341 1342 c = 1 / sqrt(1 + tau * tau); 1343 s = c * tau; 1344 } 1345 } 1346 1347 // Apply Givens rotation to T 1348 CeedGivensRotation(mat_T, c, s, CEED_NOTRANSPOSE, k, k + 1, n, n); 1349 CeedGivensRotation(mat_T, c, s, CEED_TRANSPOSE, k, k + 1, n, n); 1350 1351 // Apply Givens rotation to Q 1352 CeedGivensRotation(mat, c, s, CEED_NOTRANSPOSE, k, k + 1, n, n); 1353 1354 // Update x, z 1355 if (k < n - q - 2) { 1356 x = mat_T[k + n * (k + 1)]; 1357 z = mat_T[k + n * (k + 2)]; 1358 } 1359 } 1360 itr++; 1361 } 1362 1363 // Save eigenvalues 1364 for (CeedInt i = 0; i < n; i++) lambda[i] = mat_T[i + n * i]; 1365 1366 // Check convergence 1367 CeedCheck(itr < max_itr || q > n, ceed, CEED_ERROR_MINOR, "Symmetric QR failed to converge"); 1368 return CEED_ERROR_SUCCESS; 1369 } 1370 CeedPragmaOptimizeOn 1371 1372 /** 1373 @brief Return Simultaneous Diagonalization of two matrices. 1374 1375 This solves the generalized eigenvalue problem `A x = lambda B x`, where `A` and `B` are symmetric and `B` is positive definite. 1376 We generate the matrix `X` and vector `Lambda` such that `X^T A X = Lambda` and `X^T B X = I`. 1377 This is equivalent to the LAPACK routine 'sygv' with `TYPE = 1`. 1378 1379 @param[in] ceed `Ceed` context for error handling 1380 @param[in] mat_A Row-major matrix to be factorized with eigenvalues 1381 @param[in] mat_B Row-major matrix to be factorized to identity 1382 @param[out] mat_X Row-major orthogonal matrix 1383 @param[out] lambda Vector of length `n` of generalized eigenvalues 1384 @param[in] n Number of rows/columns 1385 1386 @return An error code: 0 - success, otherwise - failure 1387 1388 @ref Utility 1389 **/ 1390 CeedPragmaOptimizeOff 1391 int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *mat_X, CeedScalar *lambda, CeedInt n) { 1392 CeedScalar *mat_C, *mat_G, *vec_D; 1393 1394 CeedCall(CeedCalloc(n * n, &mat_C)); 1395 CeedCall(CeedCalloc(n * n, &mat_G)); 1396 CeedCall(CeedCalloc(n, &vec_D)); 1397 1398 // Compute B = G D G^T 1399 memcpy(mat_G, mat_B, n * n * sizeof(mat_B[0])); 1400 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n)); 1401 1402 // Sort eigenvalues 1403 for (CeedInt i = n - 1; i >= 0; i--) { 1404 for (CeedInt j = 0; j < i; j++) { 1405 if (fabs(vec_D[j]) > fabs(vec_D[j + 1])) { 1406 CeedScalarSwap(vec_D[j], vec_D[j + 1]); 1407 for (CeedInt k = 0; k < n; k++) CeedScalarSwap(mat_G[k * n + j], mat_G[k * n + j + 1]); 1408 } 1409 } 1410 } 1411 1412 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T 1413 // = D^-1/2 G^T A G D^-1/2 1414 // -- D = D^-1/2 1415 for (CeedInt i = 0; i < n; i++) vec_D[i] = 1. / sqrt(vec_D[i]); 1416 // -- G = G D^-1/2 1417 // -- C = D^-1/2 G^T 1418 for (CeedInt i = 0; i < n; i++) { 1419 for (CeedInt j = 0; j < n; j++) { 1420 mat_G[i * n + j] *= vec_D[j]; 1421 mat_C[j * n + i] = mat_G[i * n + j]; 1422 } 1423 } 1424 // -- X = (D^-1/2 G^T) A 1425 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_C, (const CeedScalar *)mat_A, mat_X, n, n, n)); 1426 // -- C = (D^-1/2 G^T A) (G D^-1/2) 1427 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_X, (const CeedScalar *)mat_G, mat_C, n, n, n)); 1428 1429 // Compute Q^T C Q = lambda 1430 CeedCall(CeedSymmetricSchurDecomposition(ceed, mat_C, lambda, n)); 1431 1432 // Sort eigenvalues 1433 for (CeedInt i = n - 1; i >= 0; i--) { 1434 for (CeedInt j = 0; j < i; j++) { 1435 if (fabs(lambda[j]) > fabs(lambda[j + 1])) { 1436 CeedScalarSwap(lambda[j], lambda[j + 1]); 1437 for (CeedInt k = 0; k < n; k++) CeedScalarSwap(mat_C[k * n + j], mat_C[k * n + j + 1]); 1438 } 1439 } 1440 } 1441 1442 // Set X = (G D^1/2)^-T Q 1443 // = G D^-1/2 Q 1444 CeedCall(CeedMatrixMatrixMultiply(ceed, (const CeedScalar *)mat_G, (const CeedScalar *)mat_C, mat_X, n, n, n)); 1445 1446 // Cleanup 1447 CeedCall(CeedFree(&mat_C)); 1448 CeedCall(CeedFree(&mat_G)); 1449 CeedCall(CeedFree(&vec_D)); 1450 return CEED_ERROR_SUCCESS; 1451 } 1452 CeedPragmaOptimizeOn 1453 1454 /// @} 1455 1456 /// ---------------------------------------------------------------------------- 1457 /// CeedBasis Public API 1458 /// ---------------------------------------------------------------------------- 1459 /// @addtogroup CeedBasisUser 1460 /// @{ 1461 1462 /** 1463 @brief Create a tensor-product basis for \f$H^1\f$ discretizations 1464 1465 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1466 @param[in] dim Topological dimension 1467 @param[in] num_comp Number of field components (1 for scalar fields) 1468 @param[in] P_1d Number of nodes in one dimension 1469 @param[in] Q_1d Number of quadrature points in one dimension 1470 @param[in] interp_1d Row-major (`Q_1d * P_1d`) matrix expressing the values of nodal basis functions at quadrature points 1471 @param[in] grad_1d Row-major (`Q_1d * P_1d`) matrix expressing derivatives of nodal basis functions at quadrature points 1472 @param[in] q_ref_1d Array of length `Q_1d` holding the locations of quadrature points on the 1D reference element `[-1, 1]` 1473 @param[in] q_weight_1d Array of length `Q_1d` holding the quadrature weights on the reference element 1474 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1475 1476 @return An error code: 0 - success, otherwise - failure 1477 1478 @ref User 1479 **/ 1480 int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, 1481 const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis *basis) { 1482 if (!ceed->BasisCreateTensorH1) { 1483 Ceed delegate; 1484 1485 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1486 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateTensorH1"); 1487 CeedCall(CeedBasisCreateTensorH1(delegate, dim, num_comp, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, basis)); 1488 CeedCall(CeedDestroy(&delegate)); 1489 return CEED_ERROR_SUCCESS; 1490 } 1491 1492 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis dimension must be a positive value"); 1493 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1494 CeedCheck(P_1d > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1495 CeedCheck(Q_1d > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1496 1497 CeedElemTopology topo = dim == 1 ? CEED_TOPOLOGY_LINE : dim == 2 ? CEED_TOPOLOGY_QUAD : CEED_TOPOLOGY_HEX; 1498 1499 CeedCall(CeedCalloc(1, basis)); 1500 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1501 (*basis)->ref_count = 1; 1502 (*basis)->is_tensor_basis = true; 1503 (*basis)->dim = dim; 1504 (*basis)->topo = topo; 1505 (*basis)->num_comp = num_comp; 1506 (*basis)->P_1d = P_1d; 1507 (*basis)->Q_1d = Q_1d; 1508 (*basis)->P = CeedIntPow(P_1d, dim); 1509 (*basis)->Q = CeedIntPow(Q_1d, dim); 1510 (*basis)->fe_space = CEED_FE_SPACE_H1; 1511 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_ref_1d)); 1512 CeedCall(CeedCalloc(Q_1d, &(*basis)->q_weight_1d)); 1513 if (q_ref_1d) memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d * sizeof(q_ref_1d[0])); 1514 if (q_weight_1d) memcpy((*basis)->q_weight_1d, q_weight_1d, Q_1d * sizeof(q_weight_1d[0])); 1515 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->interp_1d)); 1516 CeedCall(CeedCalloc(Q_1d * P_1d, &(*basis)->grad_1d)); 1517 if (interp_1d) memcpy((*basis)->interp_1d, interp_1d, Q_1d * P_1d * sizeof(interp_1d[0])); 1518 if (grad_1d) memcpy((*basis)->grad_1d, grad_1d, Q_1d * P_1d * sizeof(grad_1d[0])); 1519 CeedCall(ceed->BasisCreateTensorH1(dim, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, q_weight_1d, *basis)); 1520 return CEED_ERROR_SUCCESS; 1521 } 1522 1523 /** 1524 @brief Create a tensor-product \f$H^1\f$ Lagrange basis 1525 1526 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1527 @param[in] dim Topological dimension of element 1528 @param[in] num_comp Number of field components (1 for scalar fields) 1529 @param[in] P Number of Gauss-Lobatto nodes in one dimension. 1530 The polynomial degree of the resulting `Q_k` element is `k = P - 1`. 1531 @param[in] Q Number of quadrature points in one dimension. 1532 @param[in] quad_mode Distribution of the `Q` quadrature points (affects order of accuracy for the quadrature) 1533 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1534 1535 @return An error code: 0 - success, otherwise - failure 1536 1537 @ref User 1538 **/ 1539 int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, CeedBasis *basis) { 1540 // Allocate 1541 int ierr = CEED_ERROR_SUCCESS; 1542 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d; 1543 1544 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis dimension must be a positive value"); 1545 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1546 CeedCheck(P > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1547 CeedCheck(Q > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1548 1549 // Get Nodes and Weights 1550 CeedCall(CeedCalloc(P * Q, &interp_1d)); 1551 CeedCall(CeedCalloc(P * Q, &grad_1d)); 1552 CeedCall(CeedCalloc(P, &nodes)); 1553 CeedCall(CeedCalloc(Q, &q_ref_1d)); 1554 CeedCall(CeedCalloc(Q, &q_weight_1d)); 1555 if (CeedLobattoQuadrature(P, nodes, NULL) != CEED_ERROR_SUCCESS) goto cleanup; 1556 switch (quad_mode) { 1557 case CEED_GAUSS: 1558 ierr = CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); 1559 break; 1560 case CEED_GAUSS_LOBATTO: 1561 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); 1562 break; 1563 } 1564 if (ierr != CEED_ERROR_SUCCESS) goto cleanup; 1565 1566 // Build B, D matrix 1567 // Fornberg, 1998 1568 for (CeedInt i = 0; i < Q; i++) { 1569 c1 = 1.0; 1570 c3 = nodes[0] - q_ref_1d[i]; 1571 interp_1d[i * P + 0] = 1.0; 1572 for (CeedInt j = 1; j < P; j++) { 1573 c2 = 1.0; 1574 c4 = c3; 1575 c3 = nodes[j] - q_ref_1d[i]; 1576 for (CeedInt k = 0; k < j; k++) { 1577 dx = nodes[j] - nodes[k]; 1578 c2 *= dx; 1579 if (k == j - 1) { 1580 grad_1d[i * P + j] = c1 * (interp_1d[i * P + k] - c4 * grad_1d[i * P + k]) / c2; 1581 interp_1d[i * P + j] = -c1 * c4 * interp_1d[i * P + k] / c2; 1582 } 1583 grad_1d[i * P + k] = (c3 * grad_1d[i * P + k] - interp_1d[i * P + k]) / dx; 1584 interp_1d[i * P + k] = c3 * interp_1d[i * P + k] / dx; 1585 } 1586 c1 = c2; 1587 } 1588 } 1589 // Pass to CeedBasisCreateTensorH1 1590 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, interp_1d, grad_1d, q_ref_1d, q_weight_1d, basis)); 1591 cleanup: 1592 CeedCall(CeedFree(&interp_1d)); 1593 CeedCall(CeedFree(&grad_1d)); 1594 CeedCall(CeedFree(&nodes)); 1595 CeedCall(CeedFree(&q_ref_1d)); 1596 CeedCall(CeedFree(&q_weight_1d)); 1597 return CEED_ERROR_SUCCESS; 1598 } 1599 1600 /** 1601 @brief Create a non tensor-product basis for \f$H^1\f$ discretizations 1602 1603 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1604 @param[in] topo Topology of element, e.g. hypercube, simplex, etc 1605 @param[in] num_comp Number of field components (1 for scalar fields) 1606 @param[in] num_nodes Total number of nodes 1607 @param[in] num_qpts Total number of quadrature points 1608 @param[in] interp Row-major (`num_qpts * num_nodes`) matrix expressing the values of nodal basis functions at quadrature points 1609 @param[in] grad Row-major (`dim * num_qpts * num_nodes`) matrix expressing derivatives of nodal basis functions at quadrature points 1610 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element 1611 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 1612 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1613 1614 @return An error code: 0 - success, otherwise - failure 1615 1616 @ref User 1617 **/ 1618 int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 1619 const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { 1620 CeedInt P = num_nodes, Q = num_qpts, dim = 0; 1621 1622 if (!ceed->BasisCreateH1) { 1623 Ceed delegate; 1624 1625 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1626 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateH1"); 1627 CeedCall(CeedBasisCreateH1(delegate, topo, num_comp, num_nodes, num_qpts, interp, grad, q_ref, q_weight, basis)); 1628 CeedCall(CeedDestroy(&delegate)); 1629 return CEED_ERROR_SUCCESS; 1630 } 1631 1632 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1633 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1634 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1635 1636 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 1637 1638 CeedCall(CeedCalloc(1, basis)); 1639 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1640 (*basis)->ref_count = 1; 1641 (*basis)->is_tensor_basis = false; 1642 (*basis)->dim = dim; 1643 (*basis)->topo = topo; 1644 (*basis)->num_comp = num_comp; 1645 (*basis)->P = P; 1646 (*basis)->Q = Q; 1647 (*basis)->fe_space = CEED_FE_SPACE_H1; 1648 CeedCall(CeedCalloc(Q * dim, &(*basis)->q_ref_1d)); 1649 CeedCall(CeedCalloc(Q, &(*basis)->q_weight_1d)); 1650 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); 1651 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); 1652 CeedCall(CeedCalloc(Q * P, &(*basis)->interp)); 1653 CeedCall(CeedCalloc(dim * Q * P, &(*basis)->grad)); 1654 if (interp) memcpy((*basis)->interp, interp, Q * P * sizeof(interp[0])); 1655 if (grad) memcpy((*basis)->grad, grad, dim * Q * P * sizeof(grad[0])); 1656 CeedCall(ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, q_weight, *basis)); 1657 return CEED_ERROR_SUCCESS; 1658 } 1659 1660 /** 1661 @brief Create a non tensor-product basis for \f$H(\mathrm{div})\f$ discretizations 1662 1663 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1664 @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below 1665 @param[in] num_comp Number of components (usually 1 for vectors in H(div) bases) 1666 @param[in] num_nodes Total number of nodes (DoFs per element) 1667 @param[in] num_qpts Total number of quadrature points 1668 @param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of basis functions at quadrature points 1669 @param[in] div Row-major (`num_qpts * num_nodes`) matrix expressing divergence of basis functions at quadrature points 1670 @param[in] q_ref Array of length `num_qpts` * dim holding the locations of quadrature points on the reference element 1671 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 1672 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1673 1674 @return An error code: 0 - success, otherwise - failure 1675 1676 @ref User 1677 **/ 1678 int CeedBasisCreateHdiv(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 1679 const CeedScalar *div, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { 1680 CeedInt Q = num_qpts, P = num_nodes, dim = 0; 1681 1682 if (!ceed->BasisCreateHdiv) { 1683 Ceed delegate; 1684 1685 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1686 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHdiv"); 1687 CeedCall(CeedBasisCreateHdiv(delegate, topo, num_comp, num_nodes, num_qpts, interp, div, q_ref, q_weight, basis)); 1688 CeedCall(CeedDestroy(&delegate)); 1689 return CEED_ERROR_SUCCESS; 1690 } 1691 1692 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1693 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1694 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1695 1696 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 1697 1698 CeedCall(CeedCalloc(1, basis)); 1699 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1700 (*basis)->ref_count = 1; 1701 (*basis)->is_tensor_basis = false; 1702 (*basis)->dim = dim; 1703 (*basis)->topo = topo; 1704 (*basis)->num_comp = num_comp; 1705 (*basis)->P = P; 1706 (*basis)->Q = Q; 1707 (*basis)->fe_space = CEED_FE_SPACE_HDIV; 1708 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); 1709 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); 1710 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); 1711 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); 1712 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); 1713 CeedCall(CeedMalloc(Q * P, &(*basis)->div)); 1714 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); 1715 if (div) memcpy((*basis)->div, div, Q * P * sizeof(div[0])); 1716 CeedCall(ceed->BasisCreateHdiv(topo, dim, P, Q, interp, div, q_ref, q_weight, *basis)); 1717 return CEED_ERROR_SUCCESS; 1718 } 1719 1720 /** 1721 @brief Create a non tensor-product basis for \f$H(\mathrm{curl})\f$ discretizations 1722 1723 @param[in] ceed `Ceed` object used to create the `CeedBasis` 1724 @param[in] topo Topology of element (`CEED_TOPOLOGY_QUAD`, `CEED_TOPOLOGY_PRISM`, etc.), dimension of which is used in some array sizes below 1725 @param[in] num_comp Number of components (usually 1 for vectors in \f$H(\mathrm{curl})\f$ bases) 1726 @param[in] num_nodes Total number of nodes (DoFs per element) 1727 @param[in] num_qpts Total number of quadrature points 1728 @param[in] interp Row-major (`dim * num_qpts * num_nodes`) matrix expressing the values of basis functions at quadrature points 1729 @param[in] curl Row-major (`curl_comp * num_qpts * num_nodes`, `curl_comp = 1` if `dim < 3` otherwise `curl_comp = dim`) matrix expressing curl of basis functions at quadrature points 1730 @param[in] q_ref Array of length `num_qpts * dim` holding the locations of quadrature points on the reference element 1731 @param[in] q_weight Array of length `num_qpts` holding the quadrature weights on the reference element 1732 @param[out] basis Address of the variable where the newly created `CeedBasis` will be stored 1733 1734 @return An error code: 0 - success, otherwise - failure 1735 1736 @ref User 1737 **/ 1738 int CeedBasisCreateHcurl(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 1739 const CeedScalar *curl, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis *basis) { 1740 CeedInt Q = num_qpts, P = num_nodes, dim = 0, curl_comp = 0; 1741 1742 if (!ceed->BasisCreateHcurl) { 1743 Ceed delegate; 1744 1745 CeedCall(CeedGetObjectDelegate(ceed, &delegate, "Basis")); 1746 CeedCheck(delegate, ceed, CEED_ERROR_UNSUPPORTED, "Backend does not implement BasisCreateHcurl"); 1747 CeedCall(CeedBasisCreateHcurl(delegate, topo, num_comp, num_nodes, num_qpts, interp, curl, q_ref, q_weight, basis)); 1748 CeedCall(CeedDestroy(&delegate)); 1749 return CEED_ERROR_SUCCESS; 1750 } 1751 1752 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 component"); 1753 CeedCheck(num_nodes > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 node"); 1754 CeedCheck(num_qpts > 0, ceed, CEED_ERROR_DIMENSION, "CeedBasis must have at least 1 quadrature point"); 1755 1756 CeedCall(CeedBasisGetTopologyDimension(topo, &dim)); 1757 curl_comp = (dim < 3) ? 1 : dim; 1758 1759 CeedCall(CeedCalloc(1, basis)); 1760 CeedCall(CeedReferenceCopy(ceed, &(*basis)->ceed)); 1761 (*basis)->ref_count = 1; 1762 (*basis)->is_tensor_basis = false; 1763 (*basis)->dim = dim; 1764 (*basis)->topo = topo; 1765 (*basis)->num_comp = num_comp; 1766 (*basis)->P = P; 1767 (*basis)->Q = Q; 1768 (*basis)->fe_space = CEED_FE_SPACE_HCURL; 1769 CeedCall(CeedMalloc(Q * dim, &(*basis)->q_ref_1d)); 1770 CeedCall(CeedMalloc(Q, &(*basis)->q_weight_1d)); 1771 if (q_ref) memcpy((*basis)->q_ref_1d, q_ref, Q * dim * sizeof(q_ref[0])); 1772 if (q_weight) memcpy((*basis)->q_weight_1d, q_weight, Q * sizeof(q_weight[0])); 1773 CeedCall(CeedMalloc(dim * Q * P, &(*basis)->interp)); 1774 CeedCall(CeedMalloc(curl_comp * Q * P, &(*basis)->curl)); 1775 if (interp) memcpy((*basis)->interp, interp, dim * Q * P * sizeof(interp[0])); 1776 if (curl) memcpy((*basis)->curl, curl, curl_comp * Q * P * sizeof(curl[0])); 1777 CeedCall(ceed->BasisCreateHcurl(topo, dim, P, Q, interp, curl, q_ref, q_weight, *basis)); 1778 return CEED_ERROR_SUCCESS; 1779 } 1780 1781 /** 1782 @brief Create a `CeedBasis` for projection from the nodes of `basis_from` to the nodes of `basis_to`. 1783 1784 Only @ref CEED_EVAL_INTERP will be valid for the new basis, `basis_project`. 1785 For \f$H^1\f$ spaces, @ref CEED_EVAL_GRAD will also be valid. 1786 The interpolation is given by `interp_project = interp_to^+ * interp_from`, where the pseudoinverse `interp_to^+` is given by QR factorization. 1787 The gradient (for the \f$H^1\f$ case) is given by `grad_project = interp_to^+ * grad_from`. 1788 1789 Note: `basis_from` and `basis_to` must have compatible quadrature spaces. 1790 1791 Note: `basis_project` will have the same number of components as `basis_from`, regardless of the number of components that `basis_to` has. 1792 If `basis_from` has 3 components and `basis_to` has 5 components, then `basis_project` will have 3 components. 1793 1794 Note: If either `basis_from` or `basis_to` are non-tensor, then `basis_project` will also be non-tensor 1795 1796 @param[in] basis_from `CeedBasis` to prolong from 1797 @param[in] basis_to `CeedBasis` to prolong to 1798 @param[out] basis_project Address of the variable where the newly created `CeedBasis` will be stored 1799 1800 @return An error code: 0 - success, otherwise - failure 1801 1802 @ref User 1803 **/ 1804 int CeedBasisCreateProjection(CeedBasis basis_from, CeedBasis basis_to, CeedBasis *basis_project) { 1805 Ceed ceed; 1806 bool create_tensor; 1807 CeedInt dim, num_comp; 1808 CeedScalar *interp_project, *grad_project; 1809 1810 CeedCall(CeedBasisGetCeed(basis_to, &ceed)); 1811 1812 // Create projection matrix 1813 CeedCall(CeedBasisCreateProjectionMatrices(basis_from, basis_to, &interp_project, &grad_project)); 1814 1815 // Build basis 1816 { 1817 bool is_tensor_to, is_tensor_from; 1818 1819 CeedCall(CeedBasisIsTensor(basis_to, &is_tensor_to)); 1820 CeedCall(CeedBasisIsTensor(basis_from, &is_tensor_from)); 1821 create_tensor = is_tensor_from && is_tensor_to; 1822 } 1823 CeedCall(CeedBasisGetDimension(basis_to, &dim)); 1824 CeedCall(CeedBasisGetNumComponents(basis_from, &num_comp)); 1825 if (create_tensor) { 1826 CeedInt P_1d_to, P_1d_from; 1827 1828 CeedCall(CeedBasisGetNumNodes1D(basis_from, &P_1d_from)); 1829 CeedCall(CeedBasisGetNumNodes1D(basis_to, &P_1d_to)); 1830 CeedCall(CeedBasisCreateTensorH1(ceed, dim, num_comp, P_1d_from, P_1d_to, interp_project, grad_project, NULL, NULL, basis_project)); 1831 } else { 1832 // Even if basis_to and basis_from are not H1, the resulting basis is H1 for interpolation to work 1833 CeedInt num_nodes_to, num_nodes_from; 1834 CeedElemTopology topo; 1835 1836 CeedCall(CeedBasisGetTopology(basis_from, &topo)); 1837 CeedCall(CeedBasisGetNumNodes(basis_from, &num_nodes_from)); 1838 CeedCall(CeedBasisGetNumNodes(basis_to, &num_nodes_to)); 1839 CeedCall(CeedBasisCreateH1(ceed, topo, num_comp, num_nodes_from, num_nodes_to, interp_project, grad_project, NULL, NULL, basis_project)); 1840 } 1841 1842 // Cleanup 1843 CeedCall(CeedFree(&interp_project)); 1844 CeedCall(CeedFree(&grad_project)); 1845 CeedCall(CeedDestroy(&ceed)); 1846 return CEED_ERROR_SUCCESS; 1847 } 1848 1849 /** 1850 @brief Copy the pointer to a `CeedBasis`. 1851 1852 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`. 1853 This `CeedBasis` will be destroyed if `*basis_copy` is the only reference to this `CeedBasis`. 1854 1855 @param[in] basis `CeedBasis` to copy reference to 1856 @param[in,out] basis_copy Variable to store copied reference 1857 1858 @return An error code: 0 - success, otherwise - failure 1859 1860 @ref User 1861 **/ 1862 int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) { 1863 if (basis != CEED_BASIS_NONE) CeedCall(CeedBasisReference(basis)); 1864 CeedCall(CeedBasisDestroy(basis_copy)); 1865 *basis_copy = basis; 1866 return CEED_ERROR_SUCCESS; 1867 } 1868 1869 /** 1870 @brief View a `CeedBasis` 1871 1872 @param[in] basis `CeedBasis` to view 1873 @param[in] stream Stream to view to, e.g., `stdout` 1874 1875 @return An error code: 0 - success, otherwise - failure 1876 1877 @ref User 1878 **/ 1879 int CeedBasisView(CeedBasis basis, FILE *stream) { 1880 bool is_tensor_basis; 1881 CeedElemTopology topo; 1882 CeedFESpace fe_space; 1883 1884 // Basis data 1885 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 1886 CeedCall(CeedBasisGetTopology(basis, &topo)); 1887 CeedCall(CeedBasisGetFESpace(basis, &fe_space)); 1888 1889 // Print FE space and element topology of the basis 1890 fprintf(stream, "CeedBasis in a %s on a %s element\n", CeedFESpaces[fe_space], CeedElemTopologies[topo]); 1891 if (is_tensor_basis) { 1892 fprintf(stream, " P: %" CeedInt_FMT "\n Q: %" CeedInt_FMT "\n", basis->P_1d, basis->Q_1d); 1893 } else { 1894 fprintf(stream, " P: %" CeedInt_FMT "\n Q: %" CeedInt_FMT "\n", basis->P, basis->Q); 1895 } 1896 fprintf(stream, " dimension: %" CeedInt_FMT "\n field components: %" CeedInt_FMT "\n", basis->dim, basis->num_comp); 1897 // Print quadrature data, interpolation/gradient/divergence/curl of the basis 1898 if (is_tensor_basis) { // tensor basis 1899 CeedInt P_1d, Q_1d; 1900 const CeedScalar *q_ref_1d, *q_weight_1d, *interp_1d, *grad_1d; 1901 1902 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d)); 1903 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d)); 1904 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d)); 1905 CeedCall(CeedBasisGetQWeights(basis, &q_weight_1d)); 1906 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d)); 1907 CeedCall(CeedBasisGetGrad1D(basis, &grad_1d)); 1908 1909 CeedCall(CeedScalarView("qref1d", "\t% 12.8f", 1, Q_1d, q_ref_1d, stream)); 1910 CeedCall(CeedScalarView("qweight1d", "\t% 12.8f", 1, Q_1d, q_weight_1d, stream)); 1911 CeedCall(CeedScalarView("interp1d", "\t% 12.8f", Q_1d, P_1d, interp_1d, stream)); 1912 CeedCall(CeedScalarView("grad1d", "\t% 12.8f", Q_1d, P_1d, grad_1d, stream)); 1913 } else { // non-tensor basis 1914 CeedInt P, Q, dim, q_comp; 1915 const CeedScalar *q_ref, *q_weight, *interp, *grad, *div, *curl; 1916 1917 CeedCall(CeedBasisGetNumNodes(basis, &P)); 1918 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &Q)); 1919 CeedCall(CeedBasisGetDimension(basis, &dim)); 1920 CeedCall(CeedBasisGetQRef(basis, &q_ref)); 1921 CeedCall(CeedBasisGetQWeights(basis, &q_weight)); 1922 CeedCall(CeedBasisGetInterp(basis, &interp)); 1923 CeedCall(CeedBasisGetGrad(basis, &grad)); 1924 CeedCall(CeedBasisGetDiv(basis, &div)); 1925 CeedCall(CeedBasisGetCurl(basis, &curl)); 1926 1927 CeedCall(CeedScalarView("qref", "\t% 12.8f", 1, Q * dim, q_ref, stream)); 1928 CeedCall(CeedScalarView("qweight", "\t% 12.8f", 1, Q, q_weight, stream)); 1929 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp)); 1930 CeedCall(CeedScalarView("interp", "\t% 12.8f", q_comp * Q, P, interp, stream)); 1931 if (grad) { 1932 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp)); 1933 CeedCall(CeedScalarView("grad", "\t% 12.8f", q_comp * Q, P, grad, stream)); 1934 } 1935 if (div) { 1936 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_DIV, &q_comp)); 1937 CeedCall(CeedScalarView("div", "\t% 12.8f", q_comp * Q, P, div, stream)); 1938 } 1939 if (curl) { 1940 CeedCall(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_CURL, &q_comp)); 1941 CeedCall(CeedScalarView("curl", "\t% 12.8f", q_comp * Q, P, curl, stream)); 1942 } 1943 } 1944 return CEED_ERROR_SUCCESS; 1945 } 1946 1947 /** 1948 @brief Apply basis evaluation from nodes to quadrature points or vice versa 1949 1950 @param[in] basis `CeedBasis` to evaluate 1951 @param[in] num_elem The number of elements to apply the basis evaluation to; 1952 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1953 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature points; 1954 @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes 1955 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, 1956 @ref CEED_EVAL_INTERP to use interpolated values, 1957 @ref CEED_EVAL_GRAD to use gradients, 1958 @ref CEED_EVAL_DIV to use divergence, 1959 @ref CEED_EVAL_CURL to use curl, 1960 @ref CEED_EVAL_WEIGHT to use quadrature weights 1961 @param[in] u Input `CeedVector` 1962 @param[out] v Output `CeedVector` 1963 1964 @return An error code: 0 - success, otherwise - failure 1965 1966 @ref User 1967 **/ 1968 int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 1969 CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v)); 1970 CeedCheck(basis->Apply, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not support CeedBasisApply"); 1971 CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v)); 1972 return CEED_ERROR_SUCCESS; 1973 } 1974 1975 /** 1976 @brief Apply basis evaluation from quadrature points to nodes and sum into target vector 1977 1978 @param[in] basis `CeedBasis` to evaluate 1979 @param[in] num_elem The number of elements to apply the basis evaluation to; 1980 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 1981 @param[in] t_mode @ref CEED_TRANSPOSE to apply the transpose, mapping from quadrature points to nodes; 1982 @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAdd()` 1983 @param[in] eval_mode @ref CEED_EVAL_NONE to use values directly, 1984 @ref CEED_EVAL_INTERP to use interpolated values, 1985 @ref CEED_EVAL_GRAD to use gradients, 1986 @ref CEED_EVAL_DIV to use divergence, 1987 @ref CEED_EVAL_CURL to use curl, 1988 @ref CEED_EVAL_WEIGHT to use quadrature weights 1989 @param[in] u Input `CeedVector` 1990 @param[out] v Output `CeedVector` to sum into 1991 1992 @return An error code: 0 - success, otherwise - failure 1993 1994 @ref User 1995 **/ 1996 int CeedBasisApplyAdd(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 1997 CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAdd only supports CEED_TRANSPOSE"); 1998 CeedCall(CeedBasisApplyCheckDims(basis, num_elem, t_mode, eval_mode, u, v)); 1999 CeedCheck(basis->ApplyAdd, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "Backend does not implement CeedBasisApplyAdd"); 2000 CeedCall(basis->ApplyAdd(basis, num_elem, t_mode, eval_mode, u, v)); 2001 return CEED_ERROR_SUCCESS; 2002 } 2003 2004 /** 2005 @brief Apply basis evaluation from nodes to arbitrary points 2006 2007 @param[in] basis `CeedBasis` to evaluate 2008 @param[in] num_elem The number of elements to apply the basis evaluation to; 2009 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 2010 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 2011 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 2012 @ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes 2013 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 2014 @ref CEED_EVAL_GRAD to use gradients, 2015 @ref CEED_EVAL_WEIGHT to use quadrature weights 2016 @param[in] x_ref `CeedVector` holding reference coordinates of each point 2017 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 2018 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 2019 2020 @return An error code: 0 - success, otherwise - failure 2021 2022 @ref User 2023 **/ 2024 int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, 2025 CeedVector x_ref, CeedVector u, CeedVector v) { 2026 CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2027 if (basis->ApplyAtPoints) { 2028 CeedCall(basis->ApplyAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2029 } else { 2030 CeedCall(CeedBasisApplyAtPoints_Core(basis, false, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2031 } 2032 return CEED_ERROR_SUCCESS; 2033 } 2034 2035 /** 2036 @brief Apply basis evaluation from nodes to arbitrary points and sum into target vector 2037 2038 @param[in] basis `CeedBasis` to evaluate 2039 @param[in] num_elem The number of elements to apply the basis evaluation to; 2040 the backend will specify the ordering in @ref CeedElemRestrictionCreate() 2041 @param[in] num_points Array of the number of points to apply the basis evaluation to in each element, size `num_elem` 2042 @param[in] t_mode @ref CEED_NOTRANSPOSE to evaluate from nodes to points; 2043 @ref CEED_NOTRANSPOSE is not valid for `CeedBasisApplyAddAtPoints()` 2044 @param[in] eval_mode @ref CEED_EVAL_INTERP to use interpolated values, 2045 @ref CEED_EVAL_GRAD to use gradients, 2046 @ref CEED_EVAL_WEIGHT to use quadrature weights 2047 @param[in] x_ref `CeedVector` holding reference coordinates of each point 2048 @param[in] u Input `CeedVector`, of length `num_nodes * num_comp` for @ref CEED_NOTRANSPOSE 2049 @param[out] v Output `CeedVector`, of length `num_points * num_q_comp` for @ref CEED_NOTRANSPOSE with @ref CEED_EVAL_INTERP 2050 2051 @return An error code: 0 - success, otherwise - failure 2052 2053 @ref User 2054 **/ 2055 int CeedBasisApplyAddAtPoints(CeedBasis basis, CeedInt num_elem, const CeedInt *num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, 2056 CeedVector x_ref, CeedVector u, CeedVector v) { 2057 CeedCheck(t_mode == CEED_TRANSPOSE, CeedBasisReturnCeed(basis), CEED_ERROR_UNSUPPORTED, "CeedBasisApplyAddAtPoints only supports CEED_TRANSPOSE"); 2058 CeedCall(CeedBasisApplyAtPointsCheckDims(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2059 if (basis->ApplyAddAtPoints) { 2060 CeedCall(basis->ApplyAddAtPoints(basis, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2061 } else { 2062 CeedCall(CeedBasisApplyAtPoints_Core(basis, true, num_elem, num_points, t_mode, eval_mode, x_ref, u, v)); 2063 } 2064 return CEED_ERROR_SUCCESS; 2065 } 2066 2067 /** 2068 @brief Get the `Ceed` associated with a `CeedBasis` 2069 2070 @param[in] basis `CeedBasis` 2071 @param[out] ceed Variable to store `Ceed` 2072 2073 @return An error code: 0 - success, otherwise - failure 2074 2075 @ref Advanced 2076 **/ 2077 int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) { 2078 *ceed = NULL; 2079 CeedCall(CeedReferenceCopy(CeedBasisReturnCeed(basis), ceed)); 2080 return CEED_ERROR_SUCCESS; 2081 } 2082 2083 /** 2084 @brief Return the `Ceed` associated with a `CeedBasis` 2085 2086 @param[in] basis `CeedBasis` 2087 2088 @return `Ceed` associated with the `basis` 2089 2090 @ref Advanced 2091 **/ 2092 Ceed CeedBasisReturnCeed(CeedBasis basis) { return basis->ceed; } 2093 2094 /** 2095 @brief Get dimension for given `CeedBasis` 2096 2097 @param[in] basis `CeedBasis` 2098 @param[out] dim Variable to store dimension of basis 2099 2100 @return An error code: 0 - success, otherwise - failure 2101 2102 @ref Advanced 2103 **/ 2104 int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) { 2105 *dim = basis->dim; 2106 return CEED_ERROR_SUCCESS; 2107 } 2108 2109 /** 2110 @brief Get topology for given `CeedBasis` 2111 2112 @param[in] basis `CeedBasis` 2113 @param[out] topo Variable to store topology of basis 2114 2115 @return An error code: 0 - success, otherwise - failure 2116 2117 @ref Advanced 2118 **/ 2119 int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) { 2120 *topo = basis->topo; 2121 return CEED_ERROR_SUCCESS; 2122 } 2123 2124 /** 2125 @brief Get number of components for given `CeedBasis` 2126 2127 @param[in] basis `CeedBasis` 2128 @param[out] num_comp Variable to store number of components 2129 2130 @return An error code: 0 - success, otherwise - failure 2131 2132 @ref Advanced 2133 **/ 2134 int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) { 2135 *num_comp = basis->num_comp; 2136 return CEED_ERROR_SUCCESS; 2137 } 2138 2139 /** 2140 @brief Get total number of nodes (in `dim` dimensions) of a `CeedBasis` 2141 2142 @param[in] basis `CeedBasis` 2143 @param[out] P Variable to store number of nodes 2144 2145 @return An error code: 0 - success, otherwise - failure 2146 2147 @ref Utility 2148 **/ 2149 int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) { 2150 *P = basis->P; 2151 return CEED_ERROR_SUCCESS; 2152 } 2153 2154 /** 2155 @brief Get total number of nodes (in 1 dimension) of a `CeedBasis` 2156 2157 @param[in] basis `CeedBasis` 2158 @param[out] P_1d Variable to store number of nodes 2159 2160 @return An error code: 0 - success, otherwise - failure 2161 2162 @ref Advanced 2163 **/ 2164 int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d) { 2165 CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply P_1d for non-tensor CeedBasis"); 2166 *P_1d = basis->P_1d; 2167 return CEED_ERROR_SUCCESS; 2168 } 2169 2170 /** 2171 @brief Get total number of quadrature points (in `dim` dimensions) of a `CeedBasis` 2172 2173 @param[in] basis `CeedBasis` 2174 @param[out] Q Variable to store number of quadrature points 2175 2176 @return An error code: 0 - success, otherwise - failure 2177 2178 @ref Utility 2179 **/ 2180 int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) { 2181 *Q = basis->Q; 2182 return CEED_ERROR_SUCCESS; 2183 } 2184 2185 /** 2186 @brief Get total number of quadrature points (in 1 dimension) of a `CeedBasis` 2187 2188 @param[in] basis `CeedBasis` 2189 @param[out] Q_1d Variable to store number of quadrature points 2190 2191 @return An error code: 0 - success, otherwise - failure 2192 2193 @ref Advanced 2194 **/ 2195 int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q_1d) { 2196 CeedCheck(basis->is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "Cannot supply Q_1d for non-tensor CeedBasis"); 2197 *Q_1d = basis->Q_1d; 2198 return CEED_ERROR_SUCCESS; 2199 } 2200 2201 /** 2202 @brief Get reference coordinates of quadrature points (in `dim` dimensions) of a `CeedBasis` 2203 2204 @param[in] basis `CeedBasis` 2205 @param[out] q_ref Variable to store reference coordinates of quadrature points 2206 2207 @return An error code: 0 - success, otherwise - failure 2208 2209 @ref Advanced 2210 **/ 2211 int CeedBasisGetQRef(CeedBasis basis, const CeedScalar **q_ref) { 2212 *q_ref = basis->q_ref_1d; 2213 return CEED_ERROR_SUCCESS; 2214 } 2215 2216 /** 2217 @brief Get quadrature weights of quadrature points (in `dim` dimensions) of a `CeedBasis` 2218 2219 @param[in] basis `CeedBasis` 2220 @param[out] q_weight Variable to store quadrature weights 2221 2222 @return An error code: 0 - success, otherwise - failure 2223 2224 @ref Advanced 2225 **/ 2226 int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) { 2227 *q_weight = basis->q_weight_1d; 2228 return CEED_ERROR_SUCCESS; 2229 } 2230 2231 /** 2232 @brief Get interpolation matrix of a `CeedBasis` 2233 2234 @param[in] basis `CeedBasis` 2235 @param[out] interp Variable to store interpolation matrix 2236 2237 @return An error code: 0 - success, otherwise - failure 2238 2239 @ref Advanced 2240 **/ 2241 int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) { 2242 if (!basis->interp && basis->is_tensor_basis) { 2243 // Allocate 2244 CeedCall(CeedMalloc(basis->Q * basis->P, &basis->interp)); 2245 2246 // Initialize 2247 for (CeedInt i = 0; i < basis->Q * basis->P; i++) basis->interp[i] = 1.0; 2248 2249 // Calculate 2250 for (CeedInt d = 0; d < basis->dim; d++) { 2251 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) { 2252 for (CeedInt node = 0; node < basis->P; node++) { 2253 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 2254 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 2255 2256 basis->interp[qpt * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p]; 2257 } 2258 } 2259 } 2260 } 2261 *interp = basis->interp; 2262 return CEED_ERROR_SUCCESS; 2263 } 2264 2265 /** 2266 @brief Get 1D interpolation matrix of a tensor product `CeedBasis` 2267 2268 @param[in] basis `CeedBasis` 2269 @param[out] interp_1d Variable to store interpolation matrix 2270 2271 @return An error code: 0 - success, otherwise - failure 2272 2273 @ref Backend 2274 **/ 2275 int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_1d) { 2276 bool is_tensor_basis; 2277 2278 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 2279 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "CeedBasis is not a tensor product CeedBasis"); 2280 *interp_1d = basis->interp_1d; 2281 return CEED_ERROR_SUCCESS; 2282 } 2283 2284 /** 2285 @brief Get gradient matrix of a `CeedBasis` 2286 2287 @param[in] basis `CeedBasis` 2288 @param[out] grad Variable to store gradient matrix 2289 2290 @return An error code: 0 - success, otherwise - failure 2291 2292 @ref Advanced 2293 **/ 2294 int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad) { 2295 if (!basis->grad && basis->is_tensor_basis) { 2296 // Allocate 2297 CeedCall(CeedMalloc(basis->dim * basis->Q * basis->P, &basis->grad)); 2298 2299 // Initialize 2300 for (CeedInt i = 0; i < basis->dim * basis->Q * basis->P; i++) basis->grad[i] = 1.0; 2301 2302 // Calculate 2303 for (CeedInt d = 0; d < basis->dim; d++) { 2304 for (CeedInt i = 0; i < basis->dim; i++) { 2305 for (CeedInt qpt = 0; qpt < basis->Q; qpt++) { 2306 for (CeedInt node = 0; node < basis->P; node++) { 2307 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 2308 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 2309 2310 if (i == d) basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->grad_1d[q * basis->P_1d + p]; 2311 else basis->grad[(i * basis->Q + qpt) * (basis->P) + node] *= basis->interp_1d[q * basis->P_1d + p]; 2312 } 2313 } 2314 } 2315 } 2316 } 2317 *grad = basis->grad; 2318 return CEED_ERROR_SUCCESS; 2319 } 2320 2321 /** 2322 @brief Get 1D gradient matrix of a tensor product `CeedBasis` 2323 2324 @param[in] basis `CeedBasis` 2325 @param[out] grad_1d Variable to store gradient matrix 2326 2327 @return An error code: 0 - success, otherwise - failure 2328 2329 @ref Advanced 2330 **/ 2331 int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d) { 2332 bool is_tensor_basis; 2333 2334 CeedCall(CeedBasisIsTensor(basis, &is_tensor_basis)); 2335 CeedCheck(is_tensor_basis, CeedBasisReturnCeed(basis), CEED_ERROR_MINOR, "CeedBasis is not a tensor product CeedBasis"); 2336 *grad_1d = basis->grad_1d; 2337 return CEED_ERROR_SUCCESS; 2338 } 2339 2340 /** 2341 @brief Get divergence matrix of a `CeedBasis` 2342 2343 @param[in] basis `CeedBasis` 2344 @param[out] div Variable to store divergence matrix 2345 2346 @return An error code: 0 - success, otherwise - failure 2347 2348 @ref Advanced 2349 **/ 2350 int CeedBasisGetDiv(CeedBasis basis, const CeedScalar **div) { 2351 *div = basis->div; 2352 return CEED_ERROR_SUCCESS; 2353 } 2354 2355 /** 2356 @brief Get curl matrix of a `CeedBasis` 2357 2358 @param[in] basis `CeedBasis` 2359 @param[out] curl Variable to store curl matrix 2360 2361 @return An error code: 0 - success, otherwise - failure 2362 2363 @ref Advanced 2364 **/ 2365 int CeedBasisGetCurl(CeedBasis basis, const CeedScalar **curl) { 2366 *curl = basis->curl; 2367 return CEED_ERROR_SUCCESS; 2368 } 2369 2370 /** 2371 @brief Destroy a @ref CeedBasis 2372 2373 @param[in,out] basis `CeedBasis` to destroy 2374 2375 @return An error code: 0 - success, otherwise - failure 2376 2377 @ref User 2378 **/ 2379 int CeedBasisDestroy(CeedBasis *basis) { 2380 if (!*basis || *basis == CEED_BASIS_NONE || --(*basis)->ref_count > 0) { 2381 *basis = NULL; 2382 return CEED_ERROR_SUCCESS; 2383 } 2384 if ((*basis)->Destroy) CeedCall((*basis)->Destroy(*basis)); 2385 CeedCall(CeedTensorContractDestroy(&(*basis)->contract)); 2386 CeedCall(CeedFree(&(*basis)->q_ref_1d)); 2387 CeedCall(CeedFree(&(*basis)->q_weight_1d)); 2388 CeedCall(CeedFree(&(*basis)->interp)); 2389 CeedCall(CeedFree(&(*basis)->interp_1d)); 2390 CeedCall(CeedFree(&(*basis)->grad)); 2391 CeedCall(CeedFree(&(*basis)->grad_1d)); 2392 CeedCall(CeedFree(&(*basis)->div)); 2393 CeedCall(CeedFree(&(*basis)->curl)); 2394 CeedCall(CeedVectorDestroy(&(*basis)->vec_chebyshev)); 2395 CeedCall(CeedBasisDestroy(&(*basis)->basis_chebyshev)); 2396 CeedCall(CeedDestroy(&(*basis)->ceed)); 2397 CeedCall(CeedFree(basis)); 2398 return CEED_ERROR_SUCCESS; 2399 } 2400 2401 /** 2402 @brief Construct a Gauss-Legendre quadrature 2403 2404 @param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-1` exactly) 2405 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]` 2406 @param[out] q_weight_1d Array of length `Q` to hold the weights 2407 2408 @return An error code: 0 - success, otherwise - failure 2409 2410 @ref Utility 2411 **/ 2412 int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) { 2413 CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0 * atan(1.0); 2414 2415 // Build q_ref_1d, q_weight_1d 2416 for (CeedInt i = 0; i <= Q / 2; i++) { 2417 // Guess 2418 xi = cos(PI * (CeedScalar)(2 * i + 1) / ((CeedScalar)(2 * Q))); 2419 // Pn(xi) 2420 P0 = 1.0; 2421 P1 = xi; 2422 P2 = 0.0; 2423 for (CeedInt j = 2; j <= Q; j++) { 2424 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2425 P0 = P1; 2426 P1 = P2; 2427 } 2428 // First Newton Step 2429 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2430 xi = xi - P2 / dP2; 2431 // Newton to convergence 2432 for (CeedInt k = 0; k < 100 && fabs(P2) > 10 * CEED_EPSILON; k++) { 2433 P0 = 1.0; 2434 P1 = xi; 2435 for (CeedInt j = 2; j <= Q; j++) { 2436 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2437 P0 = P1; 2438 P1 = P2; 2439 } 2440 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2441 xi = xi - P2 / dP2; 2442 } 2443 // Save xi, wi 2444 wi = 2.0 / ((1.0 - xi * xi) * dP2 * dP2); 2445 q_weight_1d[i] = wi; 2446 q_weight_1d[Q - 1 - i] = wi; 2447 q_ref_1d[i] = -xi; 2448 q_ref_1d[Q - 1 - i] = xi; 2449 } 2450 return CEED_ERROR_SUCCESS; 2451 } 2452 2453 /** 2454 @brief Construct a Gauss-Legendre-Lobatto quadrature 2455 2456 @param[in] Q Number of quadrature points (integrates polynomials of degree `2*Q-3` exactly) 2457 @param[out] q_ref_1d Array of length `Q` to hold the abscissa on `[-1, 1]` 2458 @param[out] q_weight_1d Array of length `Q` to hold the weights 2459 2460 @return An error code: 0 - success, otherwise - failure 2461 2462 @ref Utility 2463 **/ 2464 int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) { 2465 CeedScalar P0, P1, P2, dP2, d2P2, xi, wi, PI = 4.0 * atan(1.0); 2466 2467 // Build q_ref_1d, q_weight_1d 2468 // Set endpoints 2469 CeedCheck(Q > 1, NULL, CEED_ERROR_DIMENSION, "Cannot create Lobatto quadrature with Q=%" CeedInt_FMT " < 2 points", Q); 2470 wi = 2.0 / ((CeedScalar)(Q * (Q - 1))); 2471 if (q_weight_1d) { 2472 q_weight_1d[0] = wi; 2473 q_weight_1d[Q - 1] = wi; 2474 } 2475 q_ref_1d[0] = -1.0; 2476 q_ref_1d[Q - 1] = 1.0; 2477 // Interior 2478 for (CeedInt i = 1; i <= (Q - 1) / 2; i++) { 2479 // Guess 2480 xi = cos(PI * (CeedScalar)(i) / (CeedScalar)(Q - 1)); 2481 // Pn(xi) 2482 P0 = 1.0; 2483 P1 = xi; 2484 P2 = 0.0; 2485 for (CeedInt j = 2; j < Q; j++) { 2486 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2487 P0 = P1; 2488 P1 = P2; 2489 } 2490 // First Newton step 2491 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2492 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi); 2493 xi = xi - dP2 / d2P2; 2494 // Newton to convergence 2495 for (CeedInt k = 0; k < 100 && fabs(dP2) > 10 * CEED_EPSILON; k++) { 2496 P0 = 1.0; 2497 P1 = xi; 2498 for (CeedInt j = 2; j < Q; j++) { 2499 P2 = (((CeedScalar)(2 * j - 1)) * xi * P1 - ((CeedScalar)(j - 1)) * P0) / ((CeedScalar)(j)); 2500 P0 = P1; 2501 P1 = P2; 2502 } 2503 dP2 = (xi * P2 - P0) * (CeedScalar)Q / (xi * xi - 1.0); 2504 d2P2 = (2 * xi * dP2 - (CeedScalar)(Q * (Q - 1)) * P2) / (1.0 - xi * xi); 2505 xi = xi - dP2 / d2P2; 2506 } 2507 // Save xi, wi 2508 wi = 2.0 / (((CeedScalar)(Q * (Q - 1))) * P2 * P2); 2509 if (q_weight_1d) { 2510 q_weight_1d[i] = wi; 2511 q_weight_1d[Q - 1 - i] = wi; 2512 } 2513 q_ref_1d[i] = -xi; 2514 q_ref_1d[Q - 1 - i] = xi; 2515 } 2516 return CEED_ERROR_SUCCESS; 2517 } 2518 2519 /// @} 2520