Merge branch 'Rob-green_spherical' into Rob-green_metalnp

Conflicts:
	src/bi_operators/CollocationIntegrator.hpp
	src/bi_operators/NumericalIntegrator.hpp
	src/green/RegisterToFactory.hpp
	src/green/SphericalDiffuse.hpp
	src/green/dielectric_profile/Sharp.hpp
	tests/green/CMakeLists.txt

Author: Roberto Di Remigio

src/bi_operators/CollocationIntegrator.hpp

@@ -35,11 +35,12 @@
 
 #include "Config.hpp"
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "IntegratorHelperFunctions.hpp"
 #include "Element.hpp"
 #include "Exception.hpp"
+#include "AlternateSphericalDiffuse.hpp"
 #include "AnisotropicLiquid.hpp"
 #include "IonicLiquid.hpp"
 #include "MetalNP.hpp"
@@ -214,6 +215,69 @@ struct CollocationIntegrator
     }
     /**@}*/
 
+    /**@{ Single and double layer potentials for a SphericalDiffuse Green's function with alternative handling of Coulomb singularty by collocation */
+    /*! \tparam ProfilePolicy the permittivity profile for the diffuse interface
+     *  \param[in] gf Green's function
+     *  \param[in] e  list of finite elements
+     */
+    template <typename ProfilePolicy>
+    Eigen::MatrixXd singleLayer(const AlternateSphericalDiffuse<CollocationIntegrator, ProfilePolicy> & gf, const std::vector<Element> & e) const {
+        // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
+        size_t mat_size = e.size();
+        Eigen::MatrixXd S = Eigen::MatrixXd::Zero(mat_size, mat_size);
+        for (size_t i = 0; i < mat_size; ++i) {
+            // Fill diagonal
+            // Diagonal of S inside the cavity
+            double Sii_I = factor_ * std::sqrt(4 * M_PI / e[i].area());
+            double eps_r2 = 0.0;
+            std::tie(eps_r2, std::ignore) = gf.epsilon(e[i].center());
+            S(i, i) = Sii_I / eps_r2;
+            Eigen::Vector3d source = e[i].center();
+            for (size_t j = 0; j < mat_size; ++j) {
+                // Fill off-diagonal
+                Eigen::Vector3d probe = e[j].center();
+                if (i != j) S(i, j) = gf.kernelS(source, probe);
+            }
+        }
+        return S;
+    }
+    /*! \tparam ProfilePolicy the permittivity profile for the diffuse interface
+     *  \param[in] gf Green's function
+     *  \param[in] e  list of finite elements
+     */
+    template <typename ProfilePolicy>
+    Eigen::MatrixXd doubleLayer(const AlternateSphericalDiffuse<CollocationIntegrator, ProfilePolicy> & gf, const std::vector<Element> & e) const {
+        // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
+        size_t mat_size = e.size();
+        Eigen::MatrixXd D = Eigen::MatrixXd::Zero(mat_size, mat_size);
+        for (size_t i = 0; i < mat_size; ++i) {
+            // Fill diagonal
+            double area = e[i].area();
+            double radius = e[i].sphere().radius();
+            // Diagonal of S inside the cavity
+            double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
+            // Diagonal of D inside the cavity
+            double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
+            // "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
+            double invE_grad = gf.inverseEDerivative(e[i].normal(), e[i].center(), e[i].center());
+
+            double eps_r2 = 0.0;
+            std::tie(eps_r2, std::ignore) = gf.epsilon(e[i].center());
+
+            D(i, i) = Dii_I + eps_r2 * Sii_I * invE_grad;
+            Eigen::Vector3d source = e[i].center();
+            for (size_t j = 0; j < mat_size; ++j) {
+                // Fill off-diagonal
+                Eigen::Vector3d probe = e[j].center();
+                Eigen::Vector3d probeNormal = e[j].normal();
+                probeNormal.normalize();
+                if (i != j) D(i, j) = gf.kernelD(probeNormal, source, probe);
+            }
+        }
+        return D;
+    }
+    /**@}*/
+
     /**@{ Single and double layer potentials for a MetalNP Green's function by collocation */
     template <typename DerivativeTraits>
     Eigen::MatrixXd singleLayer(const MetalNP<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {

src/bi_operators/IntegratorHelperFunctions.hpp

@@ -32,7 +32,11 @@
 
 #include "Config.hpp"
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
+
+#include <boost/mpl/at.hpp>
+#include <boost/mpl/int.hpp>
+#include <boost/mpl/map.hpp>
 
 #include "Element.hpp"
 #include "QuadratureRules.hpp"

src/bi_operators/NumericalIntegrator.hpp

@@ -34,10 +34,11 @@
 
 #include "Config.hpp"
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "IntegratorHelperFunctions.hpp"
 #include "Element.hpp"
+#include "AlternateSphericalDiffuse.hpp"
 #include "AnisotropicLiquid.hpp"
 #include "SphericalDiffuse.hpp"
 #include "IonicLiquid.hpp"
@@ -186,6 +187,45 @@ struct NumericalIntegrator
 //      double eps_r2 = 0.0;
 //      std::tie(eps_r2, std::ignore) = gf.epsilon(e.center());
 
+//      return eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
+        return Eigen::MatrixXd::Zero(e.size(), e.size());
+    }
+    /**@}*/
+
+    /**@{ Single and double layer potentials for a SphericalDiffuse Green's function with alternative handling of Coulomb singularty by collocation: numerical integration of diagonal */
+    template <typename ProfilePolicy>
+    Eigen::MatrixXd singleLayer(const AlternateSphericalDiffuse<NumericalIntegrator, ProfilePolicy> & /* gf */, const std::vector<Element> & e) const {
+//      // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
+//      double area = e.area();
+//      // Diagonal of S inside the cavity
+//      double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
+//      // "Diagonal" of Coulomb singularity separation coefficient
+//      double coulomb_coeff = gf.coefficientCoulomb(e.center(), e.center());
+//      // "Diagonal" of the image Green's function
+//      double image = gf.imagePotential(e.center(), e.center());
+
+//      return (Sii_I / coulomb_coeff + image);
+        return Eigen::MatrixXd::Zero(e.size(), e.size());
+    }
+    template <typename ProfilePolicy>
+    Eigen::MatrixXd doubleLayer(const AlternateSphericalDiffuse<NumericalIntegrator, ProfilePolicy> & /* gf */, const std::vector<Element> & e) const {
+//      // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
+//      double area = e.area();
+//      double radius = e.sphere().radius();
+//      // Diagonal of S inside the cavity
+//      double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
+//      // Diagonal of D inside the cavity
+//      double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
+//      // "Diagonal" of Coulomb singularity separation coefficient
+//      double coulomb_coeff = gf.coefficientCoulomb(e.center(), e.center());
+//      // "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
+//      double coeff_grad = gf.coefficientCoulombDerivative(e.normal(), e.center(), e.center()) / std::pow(coulomb_coeff, 2);
+//      // "Diagonal" of the directional derivative of the image Green's function
+//      double image_grad = gf.imagePotentialDerivative(e.normal(), e.center(), e.center());
+
+//      double eps_r2 = 0.0;
+//      std::tie(eps_r2, std::ignore) = gf.epsilon(e.center());
+
 //      return eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
         return Eigen::MatrixXd::Zero(e.size(), e.size());
     }

src/bi_operators/PurisimaIntegrator.hpp

@@ -35,11 +35,12 @@
 
 #include "Config.hpp"
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "IntegratorHelperFunctions.hpp"
 #include "Element.hpp"
 #include "Exception.hpp"
+#include "AlternateSphericalDiffuse.hpp"
 #include "AnisotropicLiquid.hpp"
 #include "IonicLiquid.hpp"
 #include "MetalNP.hpp"
@@ -223,6 +224,69 @@ struct PurisimaIntegrator
         return D;
     }
     /**@}*/
+    
+    /**@{ Single and double layer potentials for a SphericalDiffuse Green's function with alternative handling of Coulomb singularty by collocation */
+    /*! \tparam ProfilePolicy the permittivity profile for the diffuse interface
+     *  \param[in] gf Green's function
+     *  \param[in] e  list of finite elements
+     */
+    template <typename ProfilePolicy>
+    Eigen::MatrixXd singleLayer(const AlternateSphericalDiffuse<PurisimaIntegrator, ProfilePolicy> & gf, const std::vector<Element> & e) const {
+        // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
+        size_t mat_size = e.size();
+        Eigen::MatrixXd S = Eigen::MatrixXd::Zero(mat_size, mat_size);
+        for (size_t i = 0; i < mat_size; ++i) {
+            // Fill diagonal
+            // Diagonal of S inside the cavity
+            double Sii_I = factor_ * std::sqrt(4 * M_PI / e[i].area());
+            double eps_r2 = 0.0;
+            std::tie(eps_r2, std::ignore) = gf.epsilon(e[i].center());
+            S(i, i) = Sii_I / eps_r2;
+            Eigen::Vector3d source = e[i].center();
+            for (size_t j = 0; j < mat_size; ++j) {
+                // Fill off-diagonal
+                Eigen::Vector3d probe = e[j].center();
+                if (i != j) S(i, j) = gf.kernelS(source, probe);
+            }
+        }
+        return S;
+    }
+    /*! \tparam ProfilePolicy the permittivity profile for the diffuse interface
+     *  \param[in] gf Green's function
+     *  \param[in] e  list of finite elements
+     */
+    template <typename ProfilePolicy>
+    Eigen::MatrixXd doubleLayer(const AlternateSphericalDiffuse<PurisimaIntegrator, ProfilePolicy> & gf, const std::vector<Element> & e) const {
+        // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
+        size_t mat_size = e.size();
+        Eigen::MatrixXd D = Eigen::MatrixXd::Zero(mat_size, mat_size);
+        for (size_t i = 0; i < mat_size; ++i) {
+            // Fill diagonal
+            double area = e[i].area();
+            double radius = e[i].sphere().radius();
+            // Diagonal of S inside the cavity
+            double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
+            // Diagonal of D inside the cavity
+            double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
+            // "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
+            double invE_grad = gf.inverseEDerivative(e[i].normal(), e[i].center(), e[i].center());
+
+            double eps_r2 = 0.0;
+            std::tie(eps_r2, std::ignore) = gf.epsilon(e[i].center());
+
+            D(i, i) = Dii_I + eps_r2 * Sii_I * invE_grad;
+            Eigen::Vector3d source = e[i].center();
+            for (size_t j = 0; j < mat_size; ++j) {
+                // Fill off-diagonal
+                Eigen::Vector3d probe = e[j].center();
+                Eigen::Vector3d probeNormal = e[j].normal();
+                probeNormal.normalize();
+                if (i != j) D(i, j) = gf.kernelD(probeNormal, source, probe);
+            }
+        }
+        return D;
+    }
+    /**@}*/
 
     /**@{ Single and double layer potentials for a MetalNP Green's function by collocation */
     template <typename DerivativeTraits>

src/bin/check_Coulomb_coefficient.cpp

@@ -2,7 +2,7 @@
 #include <ostream>
 #include <fstream>
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "DerivativeTypes.hpp"
 #include "UniformDielectric.hpp"

src/bin/debug_wavcav.cpp

@@ -27,7 +27,7 @@
 #include <fstream>
 #include <string>
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "PWCSolver.hpp"
 #include "Vacuum.hpp"

src/bin/plot_green_spherical-CASE1.cpp

@@ -2,7 +2,7 @@
 #include <ostream>
 #include <fstream>
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "CollocationIntegrator.hpp"
 #include "DerivativeTypes.hpp"
@@ -19,6 +19,7 @@ int main()
     sphereCenter << 0.0, 0.0, 0.0;
     double sphereRadius = 100.0;
     double width = 10.0;
+    int maxL = 30;
 
     Eigen::Vector3d source;
     source << 4.0, 0.0, 85.0;
@@ -29,7 +30,7 @@ int main()
     double zMax = 300.0;
     double step = (zMax - zMin) / nPoints;
 
-    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter);
+    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
     LOG(gf);
     std::ofstream out;
     out.open("gf_spherical_CASE1.dat");

src/bin/plot_green_spherical-CASE2.cpp

@@ -2,7 +2,7 @@
 #include <ostream>
 #include <fstream>
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "CollocationIntegrator.hpp"
 #include "DerivativeTypes.hpp"
@@ -19,6 +19,7 @@ int main()
     sphereCenter << 0.0, 0.0, 0.0;
     double sphereRadius = 100.0;
     double width = 10.0;
+    int maxL = 30;
 
     Eigen::Vector3d source;
     source << 4.0, 0.0, 150.0;
@@ -29,7 +30,7 @@ int main()
     double zMax = 300.0;
     double step = (zMax - zMin) / nPoints;
 
-    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter);
+    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
     LOG(gf);
     std::ofstream out;
     out.open("gf_spherical_CASE1.dat");

src/bin/plot_green_spherical-CASE3.cpp

@@ -2,7 +2,7 @@
 #include <ostream>
 #include <fstream>
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "CollocationIntegrator.hpp"
 #include "DerivativeTypes.hpp"
@@ -19,6 +19,7 @@ int main()
     sphereCenter << 0.0, 0.0, 0.0;
     double sphereRadius = 100.0;
     double width = 10.0;
+    int maxL = 30;
 
     Eigen::Vector3d source;
     source << 4.0, 0.0, 100.0;
@@ -29,7 +30,7 @@ int main()
     double zMax = 300.0;
     double step = (zMax - zMin) / nPoints;
 
-    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter);
+    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
     LOG(gf);
     std::ofstream out;
     out.open("gf_spherical_CASE1.dat");

src/bin/plot_green_spherical-CASE4.cpp

@@ -3,7 +3,7 @@
 #include <ostream>
 #include <fstream>
 
-#include <Eigen/Dense>
+#include <Eigen/Core>
 
 #include "CollocationIntegrator.hpp"
 #include "DerivativeTypes.hpp"
@@ -20,6 +20,7 @@ int main()
     sphereCenter << 0.0, 0.0, 0.0;
     double sphereRadius = 100.0;
     double width = 10.0;
+    int maxL = 30;
 
     Eigen::Vector3d source;
     source << 0.0, 0.0, 1.0;
@@ -30,7 +31,7 @@ int main()
     double zMax = 300.0;
     double step = (zMax - zMin) / nPoints;
 
-    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter);
+    SphericalDiffuse<CollocationIntegrator, OneLayerTanh> gf(epsInside, epsOutside, width, sphereRadius, sphereCenter, maxL);
     LOG(gf);
     std::ofstream out;
     out.open("gf_spherical_CASE1.dat");