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