summaryrefslogtreecommitdiff
path: root/src/basis.c
diff options
context:
space:
mode:
authorAnton Khirnov <anton@khirnov.net>2018-04-07 18:13:49 +0200
committerAnton Khirnov <anton@khirnov.net>2018-04-07 18:13:49 +0200
commitd53f73b7f7c728e96ffc07b55fa30b9e6fc5121c (patch)
treec79d588e646ea3e65b2db1532d65bb691ee88b0e /src/basis.c
Initial commit.
Diffstat (limited to 'src/basis.c')
-rw-r--r--src/basis.c281
1 files changed, 281 insertions, 0 deletions
diff --git a/src/basis.c b/src/basis.c
new file mode 100644
index 0000000..8e5bdcc
--- /dev/null
+++ b/src/basis.c
@@ -0,0 +1,281 @@
+/*
+ * Basis sets for pseudospectral methods
+ * Copyright (C) 2016 Anton Khirnov <anton@khirnov.net>
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <errno.h>
+#include <math.h>
+
+#include "basis.h"
+#include "common.h"
+
+typedef struct BasisSet {
+ /* evaluate the idx-th basis function at the specified point*/
+ double (*eval) (const MDBasisSetContext *s, double coord, unsigned int idx);
+ /* evaluate the first derivative of the idx-th basis function at the specified point*/
+ double (*eval_diff1)(const MDBasisSetContext *s, double coord, unsigned int idx);
+ /* evaluate the second derivative of the idx-th basis function at the specified point*/
+ double (*eval_diff2)(const MDBasisSetContext *s, double coord, unsigned int idx);
+ /**
+ * Get the idx-th collocation point for the specified order.
+ * idx runs from 0 to order - 1 (inclusive)
+ */
+ double (*colloc_point)(const MDBasisSetContext *s, unsigned int order, unsigned int idx);
+} BasisSet;
+
+struct MDBasisSetContext {
+ const BasisSet *bs;
+ double sf;
+};
+
+/*
+ * The basis of even (n = 2 * idx) SB functions (Boyd 2000, Ch 17.9)
+ * SB(x, n) = sin((n + 1) arccot(|x| / L))
+ * They are symmetric wrt origin and decay as 1/x in infinity.
+ */
+static double sb_even_eval(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ double val = atan2(s->sf, coord);
+
+ idx *= 2; // even only
+
+ return sin((idx + 1) * val);
+}
+
+static double sb_even_eval_diff1(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ double val = atan2(s->sf, coord);
+
+ idx *= 2; // even only
+
+ return -s->sf * (idx + 1) * cos((idx + 1) * val) / (SQR(s->sf) + SQR(coord));
+}
+
+static double sb_even_eval_diff2(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ const double sf = s->sf;
+ double val = atan2(sf, coord);
+
+ idx *= 2; // even only
+
+ 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(const MDBasisSetContext *s, unsigned int order, unsigned int idx)
+{
+ double t;
+
+ idx = order - idx - 1;
+ //order *= 2;
+
+ //t = (idx + 2) * M_PI / (order + 4);
+#if MD_POLAR
+ t = (idx + 2) * M_PI / (2 * order + 3);
+#else
+ t = (idx + 2) * M_PI / (2 * order + 2);
+#endif
+ return s->sf / tan(t);
+}
+
+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,
+};
+
+static double sb_odd_eval(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ double val = atan2(s->sf, coord);
+
+ idx = 2 * idx + 2; // odd only
+
+ return sin((idx) * val);
+}
+
+static double sb_odd_eval_diff1(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ double val = atan2(s->sf, coord);
+
+ idx = 2 * idx + 2; // odd only
+
+ return -s->sf * (idx) * cos((idx) * val) / (SQR(s->sf) + SQR(coord));
+}
+
+static double sb_odd_eval_diff2(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ const double sf = s->sf;
+ double val = atan2(sf, coord);
+
+ idx = 2 * idx + 2; // odd only
+
+ return sf * (idx) * (2 * coord * cos((idx) * val) - sf * (idx) * sin((idx) * val)) / SQR(SQR(sf) + SQR(coord));
+}
+
+static double sb_odd_colloc_point(const MDBasisSetContext *s, unsigned int order, unsigned int idx)
+{
+ double t;
+
+ idx = order - idx - 1;
+ //order *= 2;
+
+ //t = (idx + 2) * M_PI / (order + 4);
+#if MD_POLAR
+ t = (idx + 2) * M_PI / (2 * order + 3);
+#else
+ t = (idx + 2) * M_PI / (2 * order + 3);
+#endif
+ return s->sf / tan(t);
+}
+
+static const BasisSet sb_odd_basis = {
+ .eval = sb_odd_eval,
+ .eval_diff1 = sb_odd_eval_diff1,
+ .eval_diff2 = sb_odd_eval_diff2,
+ .colloc_point = sb_odd_colloc_point,
+};
+
+static double tb_even_eval(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(s->sf / fabs(coord));
+
+ idx++;
+ idx *= 2; // even only
+
+ return cos(idx * val) - 1.0;
+}
+
+static double tb_even_eval_diff1(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ double val = (coord == 0.0) ? M_PI_2 : atan(s->sf / fabs(coord));
+
+ idx++;
+ idx *= 2; // even only
+
+ return s->sf * idx * SGN(coord) * sin(idx * val) / (SQR(s->sf) + SQR(coord));
+}
+
+static double tb_even_eval_diff2(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ const double sf = s->sf;
+ double val = (coord == 0.0) ? M_PI_2 : atan(sf / fabs(coord));
+
+ idx++;
+ idx *= 2; // even only
+
+ return -sf * idx * SGN(coord) * (2 * fabs(coord) * sin(idx * val) + sf * idx * cos(idx * val)) / SQR(SQR(sf) + SQR(coord));
+}
+
+static double tb_even_colloc_point(const MDBasisSetContext *s, unsigned int order, unsigned int idx)
+{
+ double t;
+
+ idx = order - idx - 1;
+ //order *= 2;
+
+ //t = (idx + 2) * M_PI / (order + 4);
+ t = (idx + 2) * M_PI / (2 * order + 4);
+ return s->sf / tan(t);
+}
+
+static const BasisSet tb_even_basis = {
+ .eval = tb_even_eval,
+ .eval_diff1 = tb_even_eval_diff1,
+ .eval_diff2 = tb_even_eval_diff2,
+ .colloc_point = tb_even_colloc_point,
+};
+
+static double cos_even_eval(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ return cos(2 * idx * coord);
+}
+
+static double cos_even_eval_diff1(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ return -2 * idx * sin(2 * idx * coord);
+}
+
+static double cos_even_eval_diff2(const MDBasisSetContext *s, double coord, unsigned int idx)
+{
+ return -4 * SQR(idx) * cos(2 * idx * coord);
+}
+
+static double cos_even_colloc_point(const MDBasisSetContext *s, unsigned int order, unsigned int idx)
+{
+ return M_PI * idx / (2 * order - 0);
+}
+
+static const BasisSet cos_even_basis = {
+ .eval = cos_even_eval,
+ .eval_diff1 = cos_even_eval_diff1,
+ .eval_diff2 = cos_even_eval_diff2,
+ .colloc_point = cos_even_colloc_point,
+};
+
+double md_basis_eval(const MDBasisSetContext *s, enum MDBasisEvalType type,
+ double coord, unsigned int order)
+{
+ double (*eval)(const MDBasisSetContext *, double, unsigned int) = NULL;
+
+ switch (type) {
+ case MD_BASIS_EVAL_TYPE_VALUE: eval = s->bs->eval; break;
+ case MD_BASIS_EVAL_TYPE_DIFF1: eval = s->bs->eval_diff1; break;
+ case MD_BASIS_EVAL_TYPE_DIFF2: eval = s->bs->eval_diff2; break;
+ }
+
+ return eval(s, coord, order);
+}
+
+double md_basis_colloc_point(const MDBasisSetContext *s, unsigned int order,
+ unsigned int idx)
+{
+ return s->bs->colloc_point(s, order, idx);
+}
+
+void md_basis_free(MDBasisSetContext **pctx)
+{
+ MDBasisSetContext *ctx = *pctx;
+
+ if (!ctx)
+ return;
+
+ free(ctx);
+ *pctx = NULL;
+}
+
+int md_basis_init(MDBasisSetContext **pctx, enum MDBasisFamily family, double sf)
+{
+ MDBasisSetContext *ctx;
+
+ ctx = calloc(1, sizeof(*ctx));
+ if (!ctx)
+ return -ENOMEM;
+
+ switch (family) {
+ case MD_BASIS_FAMILY_TB_EVEN: ctx->bs = &tb_even_basis; break;
+ case MD_BASIS_FAMILY_SB_EVEN: ctx->bs = &sb_even_basis; break;
+ case MD_BASIS_FAMILY_SB_ODD: ctx->bs = &sb_odd_basis; break;
+ case MD_BASIS_FAMILY_COS_EVEN: ctx->bs = &cos_even_basis; break;
+ default:
+ free(ctx);
+ return -EINVAL;
+ }
+
+ ctx->sf = sf;
+
+ *pctx = ctx;
+ return 0;
+}