diff --git a/MathLib/ODE/OdeSolver.h b/MathLib/ODE/OdeSolver.h index 7e1b608b678279a5fed7d447ba087d409d38ab9e..03ef827a5a3cdccb0333613ed4c9996a06b9ac56 100644 --- a/MathLib/ODE/OdeSolver.h +++ b/MathLib/ODE/OdeSolver.h @@ -28,8 +28,7 @@ class OdeSolver { public: using Arr = std::array<double, NumEquations>; - using ConstArrRef = - Eigen::Map<const Eigen::Matrix<double, NumEquations, 1>>; + using ConstArrRef = MappedConstVector<NumEquations>; using Function = MathLib::Function<NumEquations, FunctionArguments...>; using JacobianFunction = MathLib::JacobianFunction<NumEquations, FunctionArguments...>; diff --git a/MathLib/ODE/OdeSolverFactory.h b/MathLib/ODE/OdeSolverFactory.h index de5fba09734a4b4e1f72f4203e95e6f9956a18bc..21b2b5e1d9f206cfe92e46b1021c287ea0399aae 100644 --- a/MathLib/ODE/OdeSolverFactory.h +++ b/MathLib/ODE/OdeSolverFactory.h @@ -44,10 +44,7 @@ struct Handles<N, FunctionArgument> : public MathLib::FunctionHandles // consider omission of data pointer and switch to std::function or // alike if (f) - return f(t, - Eigen::Map<const Eigen::Matrix<double, N, 1>>{y}, - Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot}, - _data); + return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot}, _data); return true; } @@ -56,9 +53,9 @@ struct Handles<N, FunctionArgument> : public MathLib::FunctionHandles { if (df) return df(t, - Eigen::Map<const Eigen::Matrix<double, N, 1>>{y}, - Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot}, - Eigen::Map<Eigen::Matrix<double, N, N>>{jac /*, order*/}, + MappedConstVector<N>{y}, + MappedVector<N>{ydot}, + MappedMatrix<N, N>{jac /*, order*/}, _data); return true; } @@ -85,9 +82,7 @@ struct Handles<N> : public MathLib::FunctionHandles { // auto ydot_ = Eigen::Map<Eigen::Matrix<double, N, 1>>{y}; // auto ydot_ = Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot}; - return f(t, - Eigen::Map<const Eigen::Matrix<double, N, 1>>{y}, - Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot}); + return f(t, MappedConstVector<N>{y}, MappedVector<N>{ydot}); } return true; } @@ -97,9 +92,9 @@ struct Handles<N> : public MathLib::FunctionHandles { if (df) return df(t, - Eigen::Map<const Eigen::Matrix<double, N, 1>>{y}, - Eigen::Map<Eigen::Matrix<double, N, 1>>{ydot}, - Eigen::Map<Eigen::Matrix<double, N, N>>{jac /*, order*/}); + MappedConstVector<N>{y}, + MappedVector<N>{ydot}, + MappedMatrix<N, N>{jac /*, order*/}); return true; } diff --git a/MathLib/ODE/declarations.h b/MathLib/ODE/declarations.h index 753e56ef393dff320319b18c75e1b567e5c9933c..3c4d7d703ce97e2c379996b0c30024240bc89d6f 100644 --- a/MathLib/ODE/declarations.h +++ b/MathLib/ODE/declarations.h @@ -14,19 +14,30 @@ namespace MathLib { +template <int M, int N> +using MappedMatrix = Eigen::Map<Eigen::Matrix<double, M, N>>; + +template <int M, int N> +using MappedConstMatrix = Eigen::Map<const Eigen::Matrix<double, M, N>>; + +template <int N> +using MappedVector = MappedMatrix<N, 1>; + +template <int N> +using MappedConstVector = MappedConstMatrix<N, 1>; + template <unsigned N, typename... FunctionArguments> using Function = bool (*)(const double t, - Eigen::Map<const Eigen::Matrix<double, N, 1>> const y, - Eigen::Map<Eigen::Matrix<double, N, 1>> ydot, + MappedConstVector<N> const y, + MappedVector<N> ydot, FunctionArguments&... arg); template <unsigned N, typename... FunctionArguments> -using JacobianFunction = - bool (*)(const double t, - Eigen::Map<const Eigen::Matrix<double, N, 1>> const y, - Eigen::Map<Eigen::Matrix<double, N, 1>> ydot, - Eigen::Map<Eigen::Matrix<double, N, N>> jac, - FunctionArguments&... arg); +using JacobianFunction = bool (*)(const double t, + MappedConstVector<N> const y, + MappedVector<N> ydot, + MappedMatrix<N, N> jac, + FunctionArguments&... arg); // This is an internal detail class FunctionHandles diff --git a/Tests/MathLib/TestCVode.cpp b/Tests/MathLib/TestCVode.cpp index 73181c776d86dcc5d76c1a493203fd2bc380619c..9175181f00ad6ad97a51775181a0b244402b1e70 100644 --- a/Tests/MathLib/TestCVode.cpp +++ b/Tests/MathLib/TestCVode.cpp @@ -16,8 +16,8 @@ #include <cstdio> bool f(const double, - Eigen::Map<const Eigen::Matrix<double, 1, 1>> const y, - Eigen::Map<Eigen::Matrix<double, 1, 1>> ydot) + MathLib::MappedConstVector<1> const y, + MathLib::MappedVector<1> ydot) { if (y[0] <= 0.0) return false; @@ -26,9 +26,9 @@ bool f(const double, } bool df(const double /*t*/, - Eigen::Map<const Eigen::Matrix<double, 1, 1>> const y, - Eigen::Map<Eigen::Matrix<double, 1, 1>> /*ydot*/, - Eigen::Map<Eigen::Matrix<double, 1, 1>> jac) + MathLib::MappedConstVector<1> const y, + MathLib::MappedVector<1> /*ydot*/, + MathLib::MappedMatrix<1, 1> jac) { if (y[0] <= 0.0) return false; @@ -42,8 +42,8 @@ struct ExtraData }; bool f_extra(const double, - Eigen::Map<const Eigen::Matrix<double, 1, 1>> const y, - Eigen::Map<Eigen::Matrix<double, 1, 1>> ydot, + MathLib::MappedConstVector<1> const y, + MathLib::MappedVector<1> ydot, ExtraData& data) { if (y[0] <= 0.0) return false;