diff options
Diffstat (limited to 'basis.c')
-rw-r--r-- | basis.c | 87 |
1 files changed, 80 insertions, 7 deletions
@@ -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; +} |