From 63b25eca76d8481a0fdad339464d1eb617826583 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Sun, 28 Aug 2016 10:07:20 +0200 Subject: qfunc: disentangle the API from the global solver API --- eval.c | 36 +++++++++++------------------------- 1 file changed, 11 insertions(+), 25 deletions(-) (limited to 'eval.c') diff --git a/eval.c b/eval.c index 6c70790..98fd0a7 100644 --- a/eval.c +++ b/eval.c @@ -1,5 +1,5 @@ /* - * Copyright 2014-2015 Anton Khirnov + * Copyright 2014-2016 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 @@ -27,6 +27,7 @@ #include "brill_data.h" #include "internal.h" +#include "qfunc.h" int bd_eval_psi(const BDContext *bd, const double *rho, int nb_coords_rho, @@ -172,14 +173,14 @@ do { \ * derivatives as needed */ #define DO_EVAL_Q_WRAPPER(DO_EVAL, eval_drho, eval_dz, eval_d2rho, eval_d2z, eval_d2rhoz) \ do { \ - const double base = psi4 * exp(2 * s->q_func->q(bd, rrho, zz)); \ + const double base = psi4 * exp(2 * bdi_qfunc_eval(s->qfunc, rrho, zz, 0)); \ double drho, dz, d2rho, d2z, d2rho_z; \ \ - if (eval_drho) drho = s->q_func->dq_rho(bd, rrho, zz) + 2 * GRID(dpsi_rho) / ppsi; \ - if (eval_dz) dz = s->q_func->dq_z(bd, rrho, zz) + 2 * GRID(dpsi_z) / ppsi; \ - if (eval_d2rho) d2rho = s->q_func->d2q_rho(bd, rrho, zz) + 2 * GRID(d2psi_rho) / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi); \ - if (eval_d2z) d2z = s->q_func->d2q_z(bd, rrho, zz) + 2 * GRID(d2psi_z) / ppsi - 2 * SQR(GRID(dpsi_z) / ppsi); \ - if (eval_d2rhoz) d2rho_z = s->q_func->d2q_rho_z(bd, rrho, zz) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \ + if (eval_drho) drho = bdi_qfunc_eval(s->qfunc, rrho, zz, 1) + 2 * GRID(dpsi_rho) / ppsi; \ + if (eval_dz) dz = bdi_qfunc_eval(s->qfunc, rrho, zz, 3) + 2 * GRID(dpsi_z) / ppsi; \ + if (eval_d2rho) d2rho = bdi_qfunc_eval(s->qfunc, rrho, zz, 2) + 2 * GRID(d2psi_rho) / ppsi - 2 * SQR(GRID(dpsi_rho) / ppsi); \ + if (eval_d2z) d2z = bdi_qfunc_eval(s->qfunc, rrho, zz, 6) + 2 * GRID(d2psi_z) / ppsi - 2 * SQR(GRID(dpsi_z) / ppsi); \ + if (eval_d2rhoz) d2rho_z = bdi_qfunc_eval(s->qfunc, rrho, zz, 4) + 2 * GRID(d2psi_rho_z) / ppsi - 2 * GRID(dpsi_rho) * GRID(dpsi_z) / psi2; \ \ DO_EVAL; \ } while (0) @@ -286,34 +287,19 @@ int bd_eval_q(const BDContext *bd, const double *rho, int nb_coords_rho, double *out, unsigned int out_stride) { BDPriv *s = bd->priv; - double (*eval)(const struct BDContext *bd, double rho, double z); + int diff_idx; if (diff_order[0] > 2 || diff_order[1] > 2) { bdi_log(bd, 0, "Derivatives of higher order than 2 are not supported.\n"); return -ENOSYS; } - switch (diff_order[0]) { - case 0: - switch (diff_order[1]) { - case 0: eval = s->q_func->q; break; - case 1: eval = s->q_func->dq_z; break; - case 2: eval = s->q_func->d2q_z; break; - } - break; - case 1: - switch (diff_order[1]) { - case 0: eval = s->q_func->dq_rho; break; - case 1: eval = s->q_func->d2q_rho_z; break; - } - break; - case 2: eval = s->q_func->d2q_rho; break; - } + diff_idx = diff_order[1] * 3 + diff_order[0]; for (int i = 0; i < nb_coords_z; i++) { double *dst = out + i * out_stride; for (int j = 0; j < nb_coords_rho; j++) - *dst++ = eval(bd, rho[j], z[i]); + *dst++ = bdi_qfunc_eval(s->qfunc, rho[j], z[i], diff_idx); } return 0; -- cgit v1.2.3