Use std::tuple to return multiple values from function

Author: Roberto Di Remigio

src/bin/plot_green_spherical.cpp

@@ -11,15 +11,15 @@
 int main()
 {
     int nPoints = 10000;
-    double epsInside = 2.0;
+    double epsInside = 80.0;
     double epsOutside = 2.0;
     Eigen::Vector3d sphereCenter;
     sphereCenter << 0.0, 0.0, 0.0;
     double sphereRadius = 100.0;
     double width = 10.0;
 
     Eigen::Vector3d source;
-    source << 4.0, 0.0, 100.0;
+    source << 0.0, 0.0, 1.0;
     Eigen::Vector3d probe;
     probe << 0.0, 0.0, 0.0;
 

src/green/InterfacesDef.hpp

@@ -29,6 +29,7 @@
 #include <array>
 #include <cmath>
 #include <functional>
+#include <tuple>
 #include <vector>
 
 #include "Config.hpp"
@@ -50,7 +51,7 @@ typedef std::array<StateType, 3> RadialFunction;
 /*! \typedef ProfileEvaluator
  *  \brief sort of a function pointer to the dielectric profile evaluation function
  */
-typedef std::function<void(double &, double &, const double)> ProfileEvaluator;
+typedef std::function< std::tuple<double, double>(const double) > ProfileEvaluator;
 
 /*! \struct IntegratorParameters
  *  \brief holds parameters for the integrator
@@ -97,7 +98,7 @@ class LnTransformedRadial
         {
             // Evaluate the dielectric profile
             double eps = 0.0, epsPrime = 0.0;
-            eval_(eps, epsPrime, r);
+            std::tie(eps, epsPrime) = eval_(r);
             if (numericalZero(eps)) throw std::domain_error("Division by zero!");
             double gamma_epsilon = epsPrime / eps;
             // System of equations is defined here

src/green/SphericalDiffuse.hpp

@@ -31,6 +31,7 @@
 #include <functional>
 #include <iosfwd>
 #include <string>
+#include <tuple>
 #include <vector>
 
 #include "Config.hpp"
@@ -120,8 +121,8 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
                               const Eigen::Vector3d & p1, const Eigen::Vector3d & p2) const
     {
         /*
-        double eps_r2 = 0.0, epsPrime_r2 = 0.0;
-        this->profile_(eps_r2, epsPrime_r2, p2.norm());
+        double eps_r2 = 0.0;
+        std::tie(eps_r2, std::ignore) = this->profile_(p2.norm());
 
         return (eps_r2 * this->derivativeProbe(direction, p1, p2));
         */
@@ -131,8 +132,8 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
         double imageDeriv = threePointStencil(std::bind(&SphericalDiffuse::imagePotential, this, _1, _2),
                 p2, p1, direction, this->delta_);
 
-        double eps_r2 = 0.0, epsPrime_r2 = 0.0;
-        this->profile_(eps_r2, epsPrime_r2, p2.norm());
+        double eps_r2 = 0.0;
+        std::tie(eps_r2, std::ignore) = this->profile_(p2.norm());
 
         return (eps_r2 * CoulombDeriv + imageDeriv);
     }
@@ -208,7 +209,7 @@ class SphericalDiffuse : public GreensFunction<Numerical, ProfilePolicy>
     double imagePotentialDerivative(const Eigen::Vector3d & direction, const Eigen::Vector3d & p1, const Eigen::Vector3d & p2) const;
     /*! Handle to the dielectric profile evaluation
      */
-    void epsilon(double & v, double & d, double point) const { this->profile_(v, d, point); }
+    void epsilon(double & v, double & d, double point) const { std::tie(v, d) = this->profile_(point); }
 private:
     /*! Evaluates the Green's function given a pair of points
      *

src/green/TanhSphericalDiffuse.hpp

@@ -30,6 +30,7 @@
 #include <functional>
 #include <iosfwd>
 #include <string>
+#include <tuple>
 
 #include "Config.hpp"
 
@@ -113,7 +114,7 @@ inline void TanhSphericalDiffuse::initSphericalDiffuse()
     double r_infinity_  = profile_.center() + 200.0; /*! Upper bound of the integration interval */
     double observer_step_ = 1.0e-03; /*! Time step between observer calls */
     IntegratorParameters params_(r_0_, r_infinity_, observer_step_);
-    ProfileEvaluator eval_ = std::bind(&TanhDiffuse::operator(), this->profile_, _1, _2, _3);
+    ProfileEvaluator eval_ = std::bind(&TanhDiffuse::operator(), this->profile_, _1);
 
     LOG("Computing coefficient for the separation of the Coulomb singularity");
     LOG("Computing first radial solution L = " + std::to_string(maxLC_));
@@ -175,8 +176,8 @@ inline double TanhSphericalDiffuse::coefficient(double r1, double r2) const
     double d_omega2 = linearInterpolation(r2, omegaC_[0], omegaC_[2]);
 
     double tmp = 0.0, coeff = 0.0;
-    double eps_r2 = 0.0, epsPrime_r2 = 0.0;
-    this->profile_(eps_r2, epsPrime_r2, r2);
+    double eps_r2 = 0.0;
+    std::tie(eps_r2, std::ignore) = this->profile_(r2);
 
     double denominator = (d_zeta2 - d_omega2) * std::pow(r2, 2) * eps_r2;
 
@@ -217,8 +218,8 @@ inline double TanhSphericalDiffuse::functionSummation(int L, double r1, double r
     /* Value of derivative of omega_[L] at point with index 2 */
     double d_omega2 = linearInterpolation(r2, omega_[L][0], omega_[L][2]);
 
-    double eps_r2 = 0.0, epsPrime_r2 = 0.0;
-    this->profile_(eps_r2, epsPrime_r2, r2);
+    double eps_r2 = 0.0;
+    std::tie(eps_r2, std::ignore) = this->profile_(r2);
 
     double denominator = (d_zeta2 - d_omega2) * std::pow(r2, 2) * eps_r2;
 

src/green/dielectric_profile/TanhDiffuse.hpp

@@ -27,6 +27,7 @@
 #define TANHDIFFUSE_HPP
 
 #include <iosfwd>
+#include <tuple>
 
 #include "Config.hpp"
 
@@ -40,16 +41,26 @@
 class TanhDiffuse
 {
 private:
+    /// Dielectric constant on the left of the interface
     double epsilon1_;
+    /// Dielectric constant one the right of the interface
     double epsilon2_;
+    /// Width of the transition layer
     double width_;
+    /// Center of the transition layer
     double center_;
+    /*! Returns value of dielectric profile at given point
+     *  \param[in] point where to evaluate the profile
+     */
     double value(double point) const {
         double epsPlus = (epsilon1_ + epsilon2_) / 2.0;
         double epsMinus = (epsilon2_ - epsilon1_) / 2.0;
         double tanh_r = std::tanh((point - center_) / width_);
         return (epsPlus + epsMinus * tanh_r); //epsilon(r)
     }
+    /*! Returns value of derivative of dielectric profile at given point
+     *  \param[in] point where to evaluate the derivative
+     */
     double derivative(double point) const {
         double factor = (epsilon1_ - epsilon2_) / (2.0 * width_);
         double tanh_r = std::tanh((point - center_) / width_);
@@ -59,24 +70,20 @@ class TanhDiffuse
     {
         os << "Permittivity inside  = " << epsilon1_ << std::endl;
         os << "Permittivity outside = " << epsilon2_ << std::endl;
-        os << "Profile width        = " << width_    << std::endl;
-        os << "Profile center       = " << center_;
+        os << "Profile width        = " << width_    << " AU" << std::endl;
+        os << "Profile center       = " << center_   << " AU";
         return os;
     }
 public:
     TanhDiffuse() {}
     TanhDiffuse(double e1, double e2, double w, double c) :
         epsilon1_(e1), epsilon2_(e2), width_(w), center_(c) {}
-    /*! The permittivity profile of the transition layer
-     *  \param[out]  e the value of the dielectric constant at point r
-     *  \param[out] de the value of the derivative of the dielectric constant
-     *                 at point r
+    /*! Returns a tuple holding the permittivity and its derivative
      *  \param[in]   r evaluation point
      */
-    void operator()(double & e, double & de, const double r) const
+    std::tuple<double, double> operator()(const double r) const
     {
-        e = value(r);
-        de = derivative(r);
+        return std::make_tuple(value(r), derivative(r));
     }
     double epsilon1() const { return epsilon1_; }
     double epsilon2() const { return epsilon2_; }