/* * Copyright 2014-2015 Anton Khirnov * * 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 . */ /** * @file * definitions of the basis functions for the spectral method */ #include #include #include "basis.h" #include "internal.h" /* 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(s->sf, coord); idx *= 2; // even only return sin((idx + 1) * val); } 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 return - sf * (idx + 1) * cos((idx + 1) * val) / (SQR(sf) + SQR(coord)); } 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 return sf * (idx + 1) * (2 * coord * cos((idx + 1) * val) - sf * (idx + 1) * sin((idx + 1) * val)) / SQR(SQR(sf) + SQR(coord)); } static void sb_eval_multi(const BasisSetContext *s, double coord, double **out, int order_start, int order_end) { const double sf = s->sf; const double val = atan2(sf, coord); for (int i = order_start; i < order_end; i++) { const int idx = 2 * i; const double sval = sin((idx + 1) * val); const double cval = cos((idx + 1) * val); out[0][i] = sval; out[1][i] = - sf * (idx + 1) * cval / (SQR(sf) + SQR(coord)); out[2][i] = sf * (idx + 1) * (2 * coord * cval - sf * (idx + 1) * sval) / SQR(SQR(sf) + SQR(coord)); } } 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 + 3); return s->sf / tan(t); } /* * 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 odd powers of x in infinity. */ static const BasisSet sb_even_basis = { .eval = sb_even_eval, .eval_diff1 = sb_even_eval_diff1, .eval_diff2 = sb_even_eval_diff2, .eval_multi = sb_eval_multi, .colloc_point = sb_even_colloc_point, }; static double cos_even_eval(const BasisSetContext *s, double coord, int idx) { return cos(2 * idx * coord); } static double cos_even_eval_diff1(const BasisSetContext *s, double coord, int idx) { return -2 * idx * sin(2 * idx * coord); } static double cos_even_eval_diff2(const BasisSetContext *s, double coord, int idx) { return -4 * SQR(idx) * cos(2 * idx * coord); } static void cos_even_eval_multi(const BasisSetContext *s, double coord, double **out, int order_start, int order_end) { for (int i = order_start; i < order_end; i++) { const double sval = sin(2 * i * coord); const double cval = cos(2 * i * coord); out[0][i] = cval; out[1][i] = -2 * i * sval; out[2][i] = -4 * SQR(i) * cval; } } static double cos_even_colloc_point(const BasisSetContext *s, int order, int idx) { return M_PI * idx / (2 * order); } static const BasisSet cos_even_basis = { .eval = cos_even_eval, .eval_diff1 = cos_even_eval_diff1, .eval_diff2 = cos_even_eval_diff2, .eval_multi = cos_even_eval_multi, .colloc_point = cos_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_eval_multi(BasisSetContext *s, double coord, double **out, int order_start, int order_end) { if (s->bs->eval_multi) s->bs->eval_multi(s, coord, out, order_start, order_end); else { for (int i = order_start; i < order_end; i++) { out[0][i] = s->bs->eval(s, coord, i); out[1][i] = s->bs->eval_diff1(s, coord, i); out[2][i] = s->bs->eval_diff2(s, coord, i); } } } 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; case BASIS_FAMILY_COS_EVEN: s->bs = &cos_even_basis; break; default: free(s); return NULL; } return s; }