ceed-basis.c (d1ca5d5ed7eef5ba04c025231ef3837c27c70d10) ceed-basis.c (c8c3fa7d27bffffddcff68a8a1d51314e0358a98)
1// Copyright (c) 2017-2022, 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>

--- 244 unchanged lines hidden (view full) ---

253 @param[in] basis CeedBasis
254 @param[out] collo_grad_1d Row-major (Q_1d * Q_1d) matrix expressing derivatives of basis functions at quadrature points
255
256 @return An error code: 0 - success, otherwise - failure
257
258 @ref Backend
259**/
260int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) {
1// Copyright (c) 2017-2022, 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>

--- 244 unchanged lines hidden (view full) ---

253 @param[in] basis CeedBasis
254 @param[out] collo_grad_1d Row-major (Q_1d * Q_1d) matrix expressing derivatives of basis functions at quadrature points
255
256 @return An error code: 0 - success, otherwise - failure
257
258 @ref Backend
259**/
260int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) {
261 int i, j, k;
262 Ceed ceed;
263 CeedInt P_1d = (basis)->P_1d, Q_1d = (basis)->Q_1d;
264 CeedScalar *interp_1d, *grad_1d, *tau;
265
266 CeedCall(CeedMalloc(Q_1d * P_1d, &interp_1d));
267 CeedCall(CeedMalloc(Q_1d * P_1d, &grad_1d));
268 CeedCall(CeedMalloc(Q_1d, &tau));
269 memcpy(interp_1d, (basis)->interp_1d, Q_1d * P_1d * sizeof(basis)->interp_1d[0]);
270 memcpy(grad_1d, (basis)->grad_1d, Q_1d * P_1d * sizeof(basis)->interp_1d[0]);
271
272 // QR Factorization, interp_1d = Q R
273 CeedCall(CeedBasisGetCeed(basis, &ceed));
274 CeedCall(CeedQRFactorization(ceed, interp_1d, tau, Q_1d, P_1d));
275 // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure.
276
261 Ceed ceed;
262 CeedInt P_1d = (basis)->P_1d, Q_1d = (basis)->Q_1d;
263 CeedScalar *interp_1d, *grad_1d, *tau;
264
265 CeedCall(CeedMalloc(Q_1d * P_1d, &interp_1d));
266 CeedCall(CeedMalloc(Q_1d * P_1d, &grad_1d));
267 CeedCall(CeedMalloc(Q_1d, &tau));
268 memcpy(interp_1d, (basis)->interp_1d, Q_1d * P_1d * sizeof(basis)->interp_1d[0]);
269 memcpy(grad_1d, (basis)->grad_1d, Q_1d * P_1d * sizeof(basis)->interp_1d[0]);
270
271 // QR Factorization, interp_1d = Q R
272 CeedCall(CeedBasisGetCeed(basis, &ceed));
273 CeedCall(CeedQRFactorization(ceed, interp_1d, tau, Q_1d, P_1d));
274 // Note: This function is for backend use, so all errors are terminal and we do not need to clean up memory on failure.
275
277 // Apply Rinv, collo_grad_1d = grad_1d Rinv
278 for (i = 0; i < Q_1d; i++) { // Row i
276 // Apply R_inv, collo_grad_1d = grad_1d R_inv
277 for (CeedInt i = 0; i < Q_1d; i++) { // Row i
279 collo_grad_1d[Q_1d * i] = grad_1d[P_1d * i] / interp_1d[0];
278 collo_grad_1d[Q_1d * i] = grad_1d[P_1d * i] / interp_1d[0];
280 for (j = 1; j < P_1d; j++) { // Column j
279 for (CeedInt j = 1; j < P_1d; j++) { // Column j
281 collo_grad_1d[j + Q_1d * i] = grad_1d[j + P_1d * i];
280 collo_grad_1d[j + Q_1d * i] = grad_1d[j + P_1d * i];
282 for (k = 0; k < j; k++) collo_grad_1d[j + Q_1d * i] -= interp_1d[j + P_1d * k] * collo_grad_1d[k + Q_1d * i];
281 for (CeedInt k = 0; k < j; k++) collo_grad_1d[j + Q_1d * i] -= interp_1d[j + P_1d * k] * collo_grad_1d[k + Q_1d * i];
283 collo_grad_1d[j + Q_1d * i] /= interp_1d[j + P_1d * j];
284 }
282 collo_grad_1d[j + Q_1d * i] /= interp_1d[j + P_1d * j];
283 }
285 for (j = P_1d; j < Q_1d; j++) collo_grad_1d[j + Q_1d * i] = 0;
284 for (CeedInt j = P_1d; j < Q_1d; j++) collo_grad_1d[j + Q_1d * i] = 0;
286 }
287
288 // Apply Q^T, collo_grad_1d = collo_grad_1d Q^T
289 CeedCall(CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, P_1d, 1, Q_1d));
290
291 CeedCall(CeedFree(&interp_1d));
292 CeedCall(CeedFree(&grad_1d));
293 CeedCall(CeedFree(&tau));

--- 651 unchanged lines hidden (view full) ---

945 @param[out] basis Address of the variable where the newly created CeedBasis will be stored.
946
947 @return An error code: 0 - success, otherwise - failure
948
949 @ref User
950**/
951int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, CeedBasis *basis) {
952 // Allocate
285 }
286
287 // Apply Q^T, collo_grad_1d = collo_grad_1d Q^T
288 CeedCall(CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, P_1d, 1, Q_1d));
289
290 CeedCall(CeedFree(&interp_1d));
291 CeedCall(CeedFree(&grad_1d));
292 CeedCall(CeedFree(&tau));

--- 651 unchanged lines hidden (view full) ---

944 @param[out] basis Address of the variable where the newly created CeedBasis will be stored.
945
946 @return An error code: 0 - success, otherwise - failure
947
948 @ref User
949**/
950int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, CeedInt P, CeedInt Q, CeedQuadMode quad_mode, CeedBasis *basis) {
951 // Allocate
953 int ierr = CEED_ERROR_SUCCESS, i, j, k;
952 int ierr = CEED_ERROR_SUCCESS;
954 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d;
955
956 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "Basis dimension must be a positive value");
957 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 component");
958 CeedCheck(P > 0, ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 node");
959 CeedCheck(Q > 0, ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 quadrature point");
960
961 // Get Nodes and Weights

--- 10 unchanged lines hidden (view full) ---

972 case CEED_GAUSS_LOBATTO:
973 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d);
974 break;
975 }
976 if (ierr != CEED_ERROR_SUCCESS) goto cleanup;
977
978 // Build B, D matrix
979 // Fornberg, 1998
953 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, *q_weight_1d;
954
955 CeedCheck(dim > 0, ceed, CEED_ERROR_DIMENSION, "Basis dimension must be a positive value");
956 CeedCheck(num_comp > 0, ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 component");
957 CeedCheck(P > 0, ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 node");
958 CeedCheck(Q > 0, ceed, CEED_ERROR_DIMENSION, "Basis must have at least 1 quadrature point");
959
960 // Get Nodes and Weights

--- 10 unchanged lines hidden (view full) ---

971 case CEED_GAUSS_LOBATTO:
972 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d);
973 break;
974 }
975 if (ierr != CEED_ERROR_SUCCESS) goto cleanup;
976
977 // Build B, D matrix
978 // Fornberg, 1998
980 for (i = 0; i < Q; i++) {
979 for (CeedInt i = 0; i < Q; i++) {
981 c1 = 1.0;
982 c3 = nodes[0] - q_ref_1d[i];
983 interp_1d[i * P + 0] = 1.0;
980 c1 = 1.0;
981 c3 = nodes[0] - q_ref_1d[i];
982 interp_1d[i * P + 0] = 1.0;
984 for (j = 1; j < P; j++) {
983 for (CeedInt j = 1; j < P; j++) {
985 c2 = 1.0;
986 c4 = c3;
987 c3 = nodes[j] - q_ref_1d[i];
984 c2 = 1.0;
985 c4 = c3;
986 c3 = nodes[j] - q_ref_1d[i];
988 for (k = 0; k < j; k++) {
987 for (CeedInt k = 0; k < j; k++) {
989 dx = nodes[j] - nodes[k];
990 c2 *= dx;
991 if (k == j - 1) {
992 grad_1d[i * P + j] = c1 * (interp_1d[i * P + k] - c4 * grad_1d[i * P + k]) / c2;
993 interp_1d[i * P + j] = -c1 * c4 * interp_1d[i * P + k] / c2;
994 }
995 grad_1d[i * P + k] = (c3 * grad_1d[i * P + k] - interp_1d[i * P + k]) / dx;
996 interp_1d[i * P + k] = c3 * interp_1d[i * P + k] / dx;

--- 351 unchanged lines hidden (view full) ---

1348 CeedSize u_length = 0, v_length;
1349 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts;
1350 CeedCall(CeedBasisGetDimension(basis, &dim));
1351 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
1352 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
1353 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
1354 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
1355 CeedCall(CeedVectorGetLength(v, &v_length));
988 dx = nodes[j] - nodes[k];
989 c2 *= dx;
990 if (k == j - 1) {
991 grad_1d[i * P + j] = c1 * (interp_1d[i * P + k] - c4 * grad_1d[i * P + k]) / c2;
992 interp_1d[i * P + j] = -c1 * c4 * interp_1d[i * P + k] / c2;
993 }
994 grad_1d[i * P + k] = (c3 * grad_1d[i * P + k] - interp_1d[i * P + k]) / dx;
995 interp_1d[i * P + k] = c3 * interp_1d[i * P + k] / dx;

--- 351 unchanged lines hidden (view full) ---

1347 CeedSize u_length = 0, v_length;
1348 CeedInt dim, num_comp, q_comp, num_nodes, num_qpts;
1349 CeedCall(CeedBasisGetDimension(basis, &dim));
1350 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
1351 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &q_comp));
1352 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
1353 CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
1354 CeedCall(CeedVectorGetLength(v, &v_length));
1356 if (u) {
1357 CeedCall(CeedVectorGetLength(u, &u_length));
1358 }
1355 if (u) CeedCall(CeedVectorGetLength(u, &u_length));
1359
1360 CeedCheck(basis->Apply, basis->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support BasisApply");
1361
1362 // Check compatibility of topological and geometrical dimensions
1363 CeedCheck((t_mode == CEED_TRANSPOSE && v_length % num_nodes == 0 && u_length % num_qpts == 0) ||
1364 (t_mode == CEED_NOTRANSPOSE && u_length % num_nodes == 0 && v_length % num_qpts == 0),
1365 basis->ceed, CEED_ERROR_DIMENSION, "Length of input/output vectors incompatible with basis dimensions");
1366

--- 15 unchanged lines hidden (view full) ---

1382 }
1383 CeedCheck(good_dims, basis->ceed, CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode");
1384
1385 CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v));
1386 return CEED_ERROR_SUCCESS;
1387}
1388
1389/**
1356
1357 CeedCheck(basis->Apply, basis->ceed, CEED_ERROR_UNSUPPORTED, "Backend does not support BasisApply");
1358
1359 // Check compatibility of topological and geometrical dimensions
1360 CeedCheck((t_mode == CEED_TRANSPOSE && v_length % num_nodes == 0 && u_length % num_qpts == 0) ||
1361 (t_mode == CEED_NOTRANSPOSE && u_length % num_nodes == 0 && v_length % num_qpts == 0),
1362 basis->ceed, CEED_ERROR_DIMENSION, "Length of input/output vectors incompatible with basis dimensions");
1363

--- 15 unchanged lines hidden (view full) ---

1379 }
1380 CeedCheck(good_dims, basis->ceed, CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode");
1381
1382 CeedCall(basis->Apply(basis, num_elem, t_mode, eval_mode, u, v));
1383 return CEED_ERROR_SUCCESS;
1384}
1385
1386/**
1387 @brief Apply basis evaluation from nodes to arbitrary points
1388
1389 @param[in] basis CeedBasis to evaluate
1390 @param[in] num_points The number of points to apply the basis evaluation to
1391 @param[in] t_mode \ref CEED_NOTRANSPOSE to evaluate from nodes to points;
1392 \ref CEED_TRANSPOSE to apply the transpose, mapping from points to nodes
1393 @param[in] eval_mode \ref CEED_EVAL_INTERP to use interpolated values,
1394 \ref CEED_EVAL_GRAD to use gradients
1395 @param[in] x_ref CeedVector holding reference coordinates of each point
1396 @param[in] u Input CeedVector, of length `num_nodes * num_comp` for `CEED_NOTRANSPOSE`
1397 @param[out] v Output CeedVector, of length `num_points * num_q_comp` for `CEED_NOTRANSPOSE` with `CEED_EVAL_INTERP`
1398
1399 @return An error code: 0 - success, otherwise - failure
1400
1401 @ref User
1402**/
1403int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector x_ref, CeedVector u,
1404 CeedVector v) {
1405 CeedSize x_length = 0, u_length = 0, v_length;
1406 CeedInt dim, num_comp, num_q_comp, num_nodes, P_1d = 1, Q_1d = 1;
1407
1408 CeedCall(CeedBasisGetDimension(basis, &dim));
1409 CeedCall(CeedBasisGetNumNodes1D(basis, &P_1d));
1410 CeedCall(CeedBasisGetNumQuadraturePoints1D(basis, &Q_1d));
1411 CeedCall(CeedBasisGetNumComponents(basis, &num_comp));
1412 CeedCall(CeedBasisGetNumQuadratureComponents(basis, eval_mode, &num_q_comp));
1413 CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
1414 CeedCall(CeedVectorGetLength(x_ref, &x_length));
1415 CeedCall(CeedVectorGetLength(v, &v_length));
1416 CeedCall(CeedVectorGetLength(u, &u_length));
1417
1418 // Check compatibility of topological and geometrical dimensions
1419 CeedCheck((t_mode == CEED_TRANSPOSE && v_length % num_nodes == 0) || (t_mode == CEED_NOTRANSPOSE && u_length % num_nodes == 0), basis->ceed,
1420 CEED_ERROR_DIMENSION, "Length of input/output vectors incompatible with basis dimensions and number of points");
1421
1422 // Check compatibility coordinates vector
1423 CeedCheck(x_length >= num_points * dim, basis->ceed, CEED_ERROR_DIMENSION,
1424 "Length of reference coordinate vector incompatible with basis dimension and number of points");
1425
1426 // Check vector lengths to prevent out of bounds issues
1427 bool good_dims = false;
1428 switch (eval_mode) {
1429 case CEED_EVAL_INTERP:
1430 good_dims = ((t_mode == CEED_TRANSPOSE && (u_length >= num_points * num_q_comp || v_length >= num_nodes * num_comp)) ||
1431 (t_mode == CEED_NOTRANSPOSE && (v_length >= num_points * num_q_comp || u_length >= num_nodes * num_comp)));
1432 break;
1433 case CEED_EVAL_GRAD:
1434 case CEED_EVAL_NONE:
1435 case CEED_EVAL_WEIGHT:
1436 case CEED_EVAL_DIV:
1437 case CEED_EVAL_CURL:
1438 // LCOV_EXCL_START
1439 return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points not supported for %s", CeedEvalModes[eval_mode]);
1440 // LCOV_EXCL_STOP
1441 }
1442 CeedCheck(good_dims, basis->ceed, CEED_ERROR_DIMENSION, "Input/output vectors too short for basis and evaluation mode");
1443
1444 // Backend method
1445 if (basis->ApplyAtPoints) {
1446 CeedCall(basis->ApplyAtPoints(basis, num_points, t_mode, eval_mode, x_ref, u, v));
1447 return CEED_ERROR_SUCCESS;
1448 }
1449
1450 // Default implementation
1451 CeedCheck(basis->is_tensor_basis, basis->ceed, CEED_ERROR_UNSUPPORTED, "Evaluation at arbitrary points only supported for tensor product bases");
1452 if (!basis->basis_chebyshev) {
1453 // Build matrix mapping from quadrature point values to Chebyshev coefficients
1454 CeedScalar *tau, *C, *I, *chebyshev_coeffs_1d;
1455 const CeedScalar *q_ref_1d;
1456
1457 // Build coefficient matrix
1458 // -- Note: Clang-tidy needs this check because it does not understand the is_tensor_basis check above
1459 CeedCheck(P_1d > 0 && Q_1d > 0, basis->ceed, CEED_ERROR_INCOMPATIBLE, "Basis dimensions are malformed");
1460 CeedCall(CeedCalloc(Q_1d * Q_1d, &C));
1461 CeedCall(CeedBasisGetQRef(basis, &q_ref_1d));
1462 for (CeedInt i = 0; i < Q_1d; i++) {
1463 const CeedScalar x = q_ref_1d[i];
1464
1465 C[i * Q_1d + 0] = 1.0;
1466 C[i * Q_1d + 1] = 2 * x;
1467 for (CeedInt j = 2; j < Q_1d; j++) C[i * Q_1d + j] = 2 * x * C[i * Q_1d + j - 1] - C[i * Q_1d + j - 2];
1468 }
1469
1470 // Inverse of coefficient matrix
1471 CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_coeffs_1d));
1472 CeedCall(CeedCalloc(Q_1d * Q_1d, &I));
1473 CeedCall(CeedCalloc(Q_1d, &tau));
1474 // -- QR Factorization, C = Q R
1475 CeedCall(CeedQRFactorization(basis->ceed, C, tau, Q_1d, Q_1d));
1476 // -- chebyshev_coeffs_1d = R_inv Q^T
1477 for (CeedInt i = 0; i < Q_1d; i++) I[i * Q_1d + i] = 1.0;
1478 // ---- Apply R_inv, chebyshev_coeffs_1d = I R_inv
1479 for (CeedInt i = 0; i < Q_1d; i++) { // Row i
1480 chebyshev_coeffs_1d[Q_1d * i] = I[Q_1d * i] / C[0];
1481 for (CeedInt j = 1; j < Q_1d; j++) { // Column j
1482 chebyshev_coeffs_1d[j + Q_1d * i] = I[j + Q_1d * i];
1483 for (CeedInt k = 0; k < j; k++) chebyshev_coeffs_1d[j + Q_1d * i] -= C[j + Q_1d * k] * chebyshev_coeffs_1d[k + Q_1d * i];
1484 chebyshev_coeffs_1d[j + Q_1d * i] /= C[j + Q_1d * j];
1485 }
1486 }
1487 // ---- Apply Q^T, chebyshev_coeffs_1d = R_inv Q^T
1488 CeedCall(CeedHouseholderApplyQ(chebyshev_coeffs_1d, C, tau, CEED_NOTRANSPOSE, Q_1d, Q_1d, Q_1d, 1, Q_1d));
1489
1490 // Build basis mapping from nodes to Chebyshev coefficients
1491 CeedScalar *chebyshev_interp_1d, *chebyshev_grad_1d, *chebyshev_q_weight_1d;
1492 const CeedScalar *interp_1d;
1493
1494 CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_interp_1d));
1495 CeedCall(CeedCalloc(Q_1d * Q_1d, &chebyshev_grad_1d));
1496 CeedCall(CeedCalloc(Q_1d, &chebyshev_q_weight_1d));
1497 CeedCall(CeedBasisGetInterp1D(basis, &interp_1d));
1498 CeedCall(CeedMatrixMatrixMultiply(basis->ceed, chebyshev_coeffs_1d, interp_1d, chebyshev_interp_1d, Q_1d, P_1d, Q_1d));
1499
1500 CeedCall(CeedVectorCreate(basis->ceed, num_comp * CeedIntPow(Q_1d, dim), &basis->vec_chebyshev));
1501 CeedCall(CeedBasisCreateTensorH1(basis->ceed, dim, num_comp, Q_1d, Q_1d, chebyshev_interp_1d, chebyshev_grad_1d, q_ref_1d, chebyshev_q_weight_1d,
1502 &basis->basis_chebyshev));
1503
1504 // Cleanup
1505 CeedCall(CeedFree(&C));
1506 CeedCall(CeedFree(&chebyshev_coeffs_1d));
1507 CeedCall(CeedFree(&I));
1508 CeedCall(CeedFree(&tau));
1509 CeedCall(CeedFree(&chebyshev_interp_1d));
1510 CeedCall(CeedFree(&chebyshev_grad_1d));
1511 CeedCall(CeedFree(&chebyshev_q_weight_1d));
1512 }
1513
1514 // Create TensorContract object if needed, such as a basis from the GPU backends
1515 if (!basis->contract) {
1516 Ceed ceed_ref;
1517 CeedBasis basis_ref;
1518
1519 CeedCall(CeedInit("/cpu/self", &ceed_ref));
1520 // Only need matching tensor contraction dimensions, any type of basis will work
1521 CeedCall(CeedBasisCreateTensorH1Lagrange(ceed_ref, dim, num_comp, Q_1d, Q_1d, CEED_GAUSS, &basis_ref));
1522 CeedCall(CeedTensorContractReference(basis_ref->contract));
1523 basis->contract = basis_ref->contract;
1524 CeedCall(CeedBasisDestroy(&basis_ref));
1525 CeedCall(CeedDestroy(&ceed_ref));
1526 }
1527
1528 // Basis evaluation
1529 switch (t_mode) {
1530 case CEED_NOTRANSPOSE: {
1531 // Nodes to arbitrary points
1532 CeedScalar *v_array;
1533 const CeedScalar *chebyshev_coeffs, *x_array_read;
1534
1535 // -- Interpolate to Chebyshev coefficients
1536 CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, u, basis->vec_chebyshev));
1537
1538 // -- Evaluate Chebyshev polynomials at arbitrary points
1539 CeedCall(CeedVectorGetArrayRead(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs));
1540 CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read));
1541 CeedCall(CeedVectorGetArrayWrite(v, CEED_MEM_HOST, &v_array));
1542 {
1543 CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
1544
1545 // ---- Values at point
1546 for (CeedInt p = 0; p < num_points; p++) {
1547 CeedInt pre = num_comp * CeedIntPow(Q_1d, dim - 1), post = 1;
1548
1549 for (CeedInt d = dim - 1; d >= 0; d--) {
1550 // ------ Compute Chebyshev polynomial values
1551 {
1552 const CeedScalar x = x_array_read[p * dim + d];
1553
1554 chebyshev_x[0] = 1.0;
1555 chebyshev_x[1] = 2 * x;
1556 for (CeedInt j = 2; j < Q_1d; j++) chebyshev_x[j] = 2 * x * chebyshev_x[j - 1] - chebyshev_x[j - 2];
1557 }
1558 // ------ Tensor contract
1559 CeedCall(CeedTensorContractApply(basis->contract, pre, Q_1d, post, 1, chebyshev_x, t_mode, false,
1560 d == (dim - 1) ? chebyshev_coeffs : tmp[d % 2], d == 0 ? &v_array[p * num_comp] : tmp[(d + 1) % 2]));
1561 pre /= Q_1d;
1562 post *= 1;
1563 }
1564 }
1565 }
1566 CeedCall(CeedVectorRestoreArrayRead(basis->vec_chebyshev, &chebyshev_coeffs));
1567 CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read));
1568 CeedCall(CeedVectorRestoreArray(v, &v_array));
1569 break;
1570 }
1571 case CEED_TRANSPOSE:
1572 return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, "CEED_TRANSPOSE unsupported for arbitrary basis point evaluation");
1573 }
1574
1575 return CEED_ERROR_SUCCESS;
1576}
1577
1578/**
1390 @brief Get Ceed associated with a CeedBasis
1391
1392 @param[in] basis CeedBasis
1393 @param[out] ceed Variable to store Ceed
1394
1395 @return An error code: 0 - success, otherwise - failure
1396
1397 @ref Advanced

--- 292 unchanged lines hidden (view full) ---

1690 CeedCall(CeedFree(&(*basis)->q_ref_1d));
1691 CeedCall(CeedFree(&(*basis)->q_weight_1d));
1692 CeedCall(CeedFree(&(*basis)->interp));
1693 CeedCall(CeedFree(&(*basis)->interp_1d));
1694 CeedCall(CeedFree(&(*basis)->grad));
1695 CeedCall(CeedFree(&(*basis)->grad_1d));
1696 CeedCall(CeedFree(&(*basis)->div));
1697 CeedCall(CeedFree(&(*basis)->curl));
1579 @brief Get Ceed associated with a CeedBasis
1580
1581 @param[in] basis CeedBasis
1582 @param[out] ceed Variable to store Ceed
1583
1584 @return An error code: 0 - success, otherwise - failure
1585
1586 @ref Advanced

--- 292 unchanged lines hidden (view full) ---

1879 CeedCall(CeedFree(&(*basis)->q_ref_1d));
1880 CeedCall(CeedFree(&(*basis)->q_weight_1d));
1881 CeedCall(CeedFree(&(*basis)->interp));
1882 CeedCall(CeedFree(&(*basis)->interp_1d));
1883 CeedCall(CeedFree(&(*basis)->grad));
1884 CeedCall(CeedFree(&(*basis)->grad_1d));
1885 CeedCall(CeedFree(&(*basis)->div));
1886 CeedCall(CeedFree(&(*basis)->curl));
1887 CeedCall(CeedVectorDestroy(&(*basis)->vec_chebyshev));
1888 CeedCall(CeedBasisDestroy(&(*basis)->basis_chebyshev));
1698 CeedCall(CeedDestroy(&(*basis)->ceed));
1699 CeedCall(CeedFree(basis));
1700 return CEED_ERROR_SUCCESS;
1701}
1702
1703/**
1704 @brief Construct a Gauss-Legendre quadrature
1705

--- 116 unchanged lines hidden --- 		1889 CeedCall(CeedDestroy(&(*basis)->ceed));
1890 CeedCall(CeedFree(basis));
1891 return CEED_ERROR_SUCCESS;
1892}
1893
1894/**
1895 @brief Construct a Gauss-Legendre quadrature
1896

--- 116 unchanged lines hidden ---