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