Draft: add basis transformation procedure

interface/ceed-basis.c

@@ -1219,6 +1219,69 @@ int CeedBasisDestroy(CeedBasis *basis) {
   return CEED_ERROR_SUCCESS;
 }
 
+/**
+  @brief Hale and Trefethen strip transformation
+
+  @param[in] Q    Number of quadrature points
+  @param[in] rho  sum of semiminor and major axis
+  @param[out] g   conformal map
+  @param[out] g_prime   derivative of conformal map g
+
+  @return An error code: 0 - success, otherwise - failure
+
+  @ref Utility
+**/
+int CeedHaleTrefethenStripMap(CeedInt Q, CeedScalar rho, CeedScalar *g, CeedScalar *g_prime) {
+  CeedScalar i, tau, d, C, PI2, PI = 4.0*atan(1.0);
+  CeedScalar u[Q];
+  tau = PI / log(rho);
+  d = .5 + 1. / (exp(tau * PI) + 1.0)
+  PI2 = PI / 2.0;
+  // u = asin(s)
+  for (int i = 0; i < Q; i++) {
+    u[i] = (asin(PI) * i) / Q;
+    // Unscaled function of u
+    g_u[i] = log(1.0 + exp(-tau * (PI2 + u[i]))) - log(1.0 + exp(-tau * (PI2 - u[i]))) + d * tau * u[i];
+    g_prime_u[i] = 1.0 / (exp(tau * (PI2 + u[i])) + 1) + 1.0 / (exp(tau * (PI2 - u[i])) + 1) - d;
+    // Normalizing factor and scaled function of s
+    C = 1.0 / log(1.0 + exp(-tau * (PI2 + PI2))) - log(1.0 + exp(-tau * (PI2 - PI2))) + d * tau * PI2;
+    g[i] = C * g_u[i];
+    g_prime[i] = -tau * C / sqrt(1.0 - sin(u[i]) * sin(u[i])) * g_prime_u[u[i]];
+  }
+  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(CeeddInt Q, CeedScalar *q_weight_1d) {
+  CeedScalar points, m_points, m_weights;
+  if (!mapping) {
+    if (!q_weight_1d)
+      points = Q;
+      else
+        points = Q;
+        weights = q_weight_1d;
+  }
+
+  m_points = ; // apply map g on quadrature points
+  if (q_weight_1d) {
+    m_weights = ; //apply derivative of map g on quadrature weights evaluated at quadrature points
+    m_points =  ;
+    m_weights = ;
+  }
+  else {
+    m_points = ;
+  }
+  return CEED_ERROR_SUCCESS;
+}
+
 /**
   @brief Construct a Gauss-Legendre quadrature
 
@@ -1270,6 +1333,8 @@ int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
     q_ref_1d[i] = -xi;
     q_ref_1d[Q-1-i]= xi;
   }
+  // Call transformed Gauss-Legendre quadrature
+
   return CEED_ERROR_SUCCESS;
 }
 
@@ -1342,6 +1407,8 @@ int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d,
     q_ref_1d[i] = -xi;
     q_ref_1d[Q-1-i]= xi;
   }
+  // Call transformed Gauss-Legendre-Lobatto quadrature
+
   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,16 @@
+/// @file
+/// Test that quadrature formula is transformed
+/// \test Test that quadrature formula is transformed
+#include <ceed.h>
+
+int main(int argc, char **argv) {
+    Ceed ceed;
+    CeedInt Q = 16;
+    CeedScalar *q_weight_1d;
+    CeedBasis basis;
+
+    CeedTransformQuadrature();
+
+    CeedDestroy(&ceed)
+    return 0;
+}