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