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