diff --git a/include/ceed/ceed.h b/include/ceed/ceed.h index 3ab9418823..0f647832c6 100644 --- a/include/ceed/ceed.h +++ b/include/ceed/ceed.h @@ -137,7 +137,7 @@ typedef enum { /// Double precision CEED_SCALAR_FP64 } CeedScalarType; -/// Base scalar type for the library to use: change which header is +/// Base scalar type for the library to use: change which header is /// included to change the precision. #include "ceed-f64.h" @@ -570,6 +570,10 @@ CEED_EXTERN int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d); CEED_EXTERN int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d); +CEED_EXTERN int CeedHaleTrefethenStripMap(CeedScalar rho, CeedScalar s, + CeedScalar *g, CeedScalar *g_prime); +CEED_EXTERN int CeedGaussHaleTrefethenQuadrature(CeedScalar rho, CeedInt Q, + CeedScalar *q_ref_1d, CeedScalar *q_weight_1d); CEED_EXTERN int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n); CEED_EXTERN int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c index bfa7f62e79..8b6cefa4ad 100644 --- a/interface/ceed-basis.c +++ b/interface/ceed-basis.c @@ -14,6 +14,8 @@ // software, applications, hardware, advanced system engineering and early // testbed platforms, in support of the nation's exascale computing imperative. +#define _XOPEN_SOURCE + #include #include #include @@ -1219,6 +1221,53 @@ int CeedBasisDestroy(CeedBasis *basis) { return CEED_ERROR_SUCCESS; } +/** + @brief Hale and Trefethen strip transformation + + @param[in] rho sum of semiminor and major axis + @param[in] s input coordinate in [-1, 1] + @param[out] g conformal map g(s); will be in [-1, 1] + @param[out] g_prime derivative of conformal map g + + @return An error code: 0 - success, otherwise - failure + + @ref Utility +**/ +int CeedHaleTrefethenStripMap(CeedScalar rho, CeedScalar s, CeedScalar *g, CeedScalar *g_prime) { + const CeedScalar u = asin(s); + const CeedScalar tau = M_PI / log(rho); + const CeedScalar + e_plus = exp(-tau*(M_PI_2 + u)), + e_minus = exp(-tau*(M_PI_2 - u)), + c2 = .5 + 1/(1 + exp(tau*M_PI)), + //C = log1p(exp(-tau*M_PI) - log1p(1) + c2*tau*M_PI_2); + C = 1 / log1p(exp(-tau*M_PI) - log1p(1) + c2*tau*M_PI_2); + //C = 1 / log1p(exp(-tau*M_PI) - log(2) + c2*tau*M_PI_2); + *g = C * (log1p(e_plus) - log1p(e_minus) + c2*tau*u); // g(s) + *g_prime = -tau*C/sqrt(1 - s*s) * (1 / (1 + e_plus) + 1 / (1 + e_minus) - c2); // g'(s) + return CEED_ERROR_SUCCESS; +} + +/** + @brief Transform quadrature by applying a smooth mapping = (g, g_prime) + + @param Q Number of quadrature points + + @return An error code: 0 - success, otherwise - failure + + @ref Utility +**/ +int CeedTransformQuadrature(CeedScalar rho, CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d) { + int ierr; + for (CeedInt i=0; i +#include + +int main(int argc, char **argv) { + Ceed ceed; + CeedInt Q = 16; + CeedScalar rho; + CeedScalar *g, *g_prime; + + CeedHaleTrefethenStripMap(); + + + CeedDestroy(&ceed); + return 0; +} \ No newline at end of file diff --git a/tests/t334-basis.c b/tests/t334-basis.c new file mode 100644 index 0000000000..afe5800783 --- /dev/null +++ b/tests/t334-basis.c @@ -0,0 +1,21 @@ +/// @file +/// Test that quadrature formula is transformed +/// \test Test that quadrature formula is transformed +#include +#include + +int main(int argc, char **argv) { + CeedInt Q = 5; + CeedScalar w_gauss[Q], w_ht[Q]; + CeedScalar q_gauss[Q], q_ht[Q]; + + CeedGaussQuadrature(Q, q_gauss, w_gauss); + CeedGaussHaleTrefethenQuadrature(1.4, Q, q_ht, w_ht); + for (CeedInt i=0; i %g (%g)\n", q_gauss[i], w_gauss[i], q_ht[i], w_ht[i]); + + CeedScalar w_sum = 0; + for (CeedInt i=0; i 1e-12) printf("Unexpected sum of weights %g\n", w_sum); + return 0; +} \ No newline at end of file