Draft: add basis transformation procedure

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,

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 <ceed/ceed.h>
 #include <ceed/backend.h>
 #include <ceed-impl.h>
@@ -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<Q; i++) {
+    CeedScalar g, g_prime;
+    ierr = CeedHaleTrefethenStripMap(rho, q_ref_1d[i], &g, &g_prime); CeedChk(ierr);
+    q_ref_1d[i] = g;
+    q_weight_1d[i] *= g_prime;
+  }
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Construct a Gauss-Legendre quadrature
 
@@ -1270,6 +1319,14 @@ int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
     q_ref_1d[i] = -xi;
     q_ref_1d[Q-1-i]= xi;
   }
+
+  return CEED_ERROR_SUCCESS;
+}
+
+int CeedGaussHaleTrefethenQuadrature(CeedScalar rho, CeedInt Q, CeedScalar *q_ref_1d,
+                        CeedScalar *q_weight_1d) {
+  CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d);
+  CeedTransformQuadrature(rho, Q, q_ref_1d, q_weight_1d);
   return CEED_ERROR_SUCCESS;
 }
 
@@ -1342,6 +1399,7 @@ int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
     q_ref_1d[i] = -xi;
     q_ref_1d[Q-1-i]= xi;
   }
+
   return CEED_ERROR_SUCCESS;
 }
 

tests/t333-basis.c

@@ -0,0 +1,18 @@
+/// @file
+/// Test that Hale-Trefethen transform  works
+/// \test Test that Hale-Trefethen transform works
+#include <ceed.h>
+#include <math.h>
+
+int main(int argc, char **argv) {
+    Ceed ceed;
+    CeedInt Q = 16;
+    CeedScalar rho;
+    CeedScalar *g, *g_prime;
+
+    CeedHaleTrefethenStripMap();
+
+
+    CeedDestroy(&ceed);
+    return 0;
+}

tests/t334-basis.c

@@ -0,0 +1,21 @@
+/// @file
+/// Test that quadrature formula is transformed
+/// \test Test that quadrature formula is transformed
+#include <ceed.h>
+#include <math.h>
+
+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<Q; i++)
+        printf("%g (%g) -> %g (%g)\n", q_gauss[i], w_gauss[i], q_ht[i], w_ht[i]);
+
+    CeedScalar w_sum = 0;
+    for (CeedInt i=0; i<Q; i++) w_sum += w_ht[i];
+    if (fabs(w_sum - 2) > 1e-12) printf("Unexpected sum of weights %g\n", w_sum);
+    return 0;
+}