basis: disentangle the API from the global solver API

1 files changed, 80 insertions, 7 deletions

diff --git a/basis.c b/basis.c
index 0735752..face638 100644
--- a/basis.c
+++ b/basis.c
@@ -21,20 +21,47 @@
  */
 
 #include <math.h>
+#include <stdlib.h>
 
+#include "basis.h"
 #include "internal.h"
 
-static double sb_even_eval(double coord, int idx, double sf)
+
+/* a set of basis functions */
+typedef struct BasisSet {
+    /* evaluate the idx-th basis function at the specified point*/
+    double (*eval)      (const BasisSetContext *s, double coord, int idx);
+    /* evaluate the first derivative of the idx-th basis function at the specified point*/
+    double (*eval_diff1)(const BasisSetContext *s, double coord, int idx);
+    /* evaluate the second derivative of the idx-th basis function at the specified point*/
+    double (*eval_diff2)(const BasisSetContext *s, double coord, int idx);
+
+    void (*eval_multi)(const BasisSetContext *s, double coord,
+                       double **out, int order_start, int order_end);
+    /**
+     * Get the idx-th collocation point for the specified order.
+     * idx runs from 0 to order - 1 (inclusive)
+     */
+    double (*colloc_point)(const BasisSetContext *s, int order, int idx);
+} BasisSet;
+
+struct BasisSetContext {
+    const struct BasisSet *bs;
+    double sf;
+};
+
+static double sb_even_eval(const BasisSetContext *s, double coord, int idx)
 {
-    double val = atan2(sf, coord);
+    double val = atan2(s->sf, coord);
 
     idx *= 2;   // even only
 
     return sin((idx + 1) * val);
 }
 
-static double sb_even_eval_diff1(double coord, int idx, double sf)
+static double sb_even_eval_diff1(const BasisSetContext *s, double coord, int idx)
 {
+    const double sf = s->sf;
     double val = atan2(sf, coord);
 
     idx *= 2;   // even only
@@ -42,8 +69,9 @@ static double sb_even_eval_diff1(double coord, int idx, double sf)
     return - sf * (idx + 1) * cos((idx + 1) * val) / (SQR(sf) + SQR(coord));
 }
 
-static double sb_even_eval_diff2(double coord, int idx, double sf)
+static double sb_even_eval_diff2(const BasisSetContext *s, double coord, int idx)
 {
+    const double sf = s->sf;
     double val = atan2(sf, coord);
 
     idx *= 2;   // even only
@@ -51,14 +79,14 @@ static double sb_even_eval_diff2(double coord, int idx, double sf)
     return sf * (idx + 1) * (2 * coord * cos((idx + 1) * val) - sf * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(sf) + SQR(coord));
 }
 
-static double sb_even_colloc_point(int order, int idx, double sf)
+static double sb_even_colloc_point(const BasisSetContext *s, int order, int idx)
 {
     double t;
 
     idx = order - idx - 1;
 
     t = (idx + 2) * M_PI / (2 * order + 2);
-    return sf / tan(t);
+    return s->sf / tan(t);
 }
 
 /*
@@ -66,9 +94,54 @@ static double sb_even_colloc_point(int order, int idx, double sf)
  * SB(x, n) = sin((n + 1) arccot(|x| / L))
  * They are symmetric wrt origin and decay as odd powers of x in infinity.
  */
-const BasisSet bdi_sb_even_basis = {
+static const BasisSet sb_even_basis = {
     .eval         = sb_even_eval,
     .eval_diff1   = sb_even_eval_diff1,
     .eval_diff2   = sb_even_eval_diff2,
     .colloc_point = sb_even_colloc_point,
 };
+
+double bdi_basis_eval(BasisSetContext *s, enum BSEvalType type, double coord, int order)
+{
+    double (*eval)(const BasisSetContext *s, double coord, int order);
+
+    switch (type) {
+    case BS_EVAL_TYPE_VALUE: eval = s->bs->eval;       break;
+    case BS_EVAL_TYPE_DIFF1: eval = s->bs->eval_diff1; break;
+    case BS_EVAL_TYPE_DIFF2: eval = s->bs->eval_diff2; break;
+    }
+
+    return eval(s, coord, order);
+}
+
+double bdi_basis_colloc_point(BasisSetContext *s, int idx, int order)
+{
+    return s->bs->colloc_point(s, order, idx);
+}
+
+void bdi_basis_free(BasisSetContext **ps)
+{
+    BasisSetContext *s = *ps;
+    free(s);
+    *ps = NULL;
+}
+
+BasisSetContext *bdi_basis_init(enum BasisFamily family, double sf)
+{
+    BasisSetContext *s = calloc(1, sizeof(*s));
+
+    if (!s)
+        return NULL;
+
+    switch (family) {
+    case BASIS_FAMILY_SB_EVEN:
+        s->bs = &sb_even_basis;
+        s->sf = sf;
+        break;
+    default:
+        free(s);
+        return NULL;
+    }
+
+    return s;
+}