Skip to content

Commit bd85a55

Browse files
author
Roberto Di Remigio
committed
Slight twist in the design of the integrator classes: this helps with the fake run-time selection metafunction.
The integrator classes return now the full matrix representation of the single and double layer operators, not only the diagonal. I have implemented only the collocation integrator, so maybe there is still some refactoring down the road... Note that not all parts compile in this commit.
1 parent 658e2e7 commit bd85a55

18 files changed

+291
-236
lines changed

‎src/bi_operators/CollocationIntegrator.hpp‎

Lines changed: 97 additions & 63 deletions
Original file line numberDiff line numberDiff line change
@@ -1,26 +1,28 @@
11
#ifndef COLLOCATIONINTEGRATOR_HPP
22
#define COLLOCATIONINTEGRATOR_HPP
33

4+
#include <cmath>
5+
#include <functional>
46
#include <iosfwd>
5-
#include <tuple>
7+
#include <stdexcept>
68

79
#include "Config.hpp"
810

911
#include <Eigen/Dense>
1012

11-
#include "DerivativeTypes.hpp"
12-
#include "AnisotropicLiquid.hpp"
13+
#include "IntegratorHelperFunctions.hpp"
1314
#include "Element.hpp"
15+
#include "AnisotropicLiquid.hpp"
1416
#include "IonicLiquid.hpp"
1517
#include "SphericalDiffuse.hpp"
1618
#include "UniformDielectric.hpp"
1719
#include "Vacuum.hpp"
1820

1921
/*! \file CollocationIntegrator.hpp
20-
* \class CollocationIntegrator
21-
* \brief Implementation of diagonal elements of S and D using approximate collocation
22+
* \struct CollocationIntegrator
23+
* \brief Implementation of the single and double layer operators matrix representation using one-point collocation
2224
* \author Roberto Di Remigio
23-
* \date 2014
25+
* \date 2015
2426
*
2527
* Calculates the diagonal elements of S as:
2628
* \f[
@@ -32,80 +34,112 @@
3234
* \f]
3335
*/
3436

35-
template <typename DerivativeTraits,
36-
typename ProfilePolicy>
37+
using std::placeholders::_1;
38+
using std::placeholders::_2;
39+
using std::placeholders::_3;
40+
using integrator::KernelS;
41+
using integrator::KernelD;
42+
3743
struct CollocationIntegrator
3844
{
3945
CollocationIntegrator() : factor_(1.07) {}
4046
~CollocationIntegrator() {}
4147

42-
double computeS(const Vacuum<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & e) const {
43-
double area = e.area();
44-
return (factor_ * std::sqrt(4 * M_PI / area));
48+
/**@{ Single and double layer potentials for a Vacuum Green's function by collocation */
49+
template <typename DerivativeTraits>
50+
Eigen::MatrixXd singleLayer(const Vacuum<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
51+
auto f = this->factor_;
52+
auto diagonal = [f] (double area) -> double { return (f * std::sqrt(4 * M_PI / area)); };
53+
KernelS kernelS = std::bind(&Vacuum<DerivativeTraits, CollocationIntegrator>::kernelS, gf, _1, _2);
54+
return integrator::singleLayer(e, diagonal, kernelS);
4555
}
46-
double computeD(const Vacuum<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & e) const {
47-
double area = e.area();
48-
double radius = e.sphere().radius();
49-
return (-factor_ * std::sqrt(M_PI/ area) * (1.0 / radius));
56+
template <typename DerivativeTraits>
57+
Eigen::MatrixXd doubleLayer(const Vacuum<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
58+
auto f = this->factor_;
59+
auto diagonal = [f] (double area, double radius) -> double { return (-f * std::sqrt(M_PI/ area) * (1.0 / radius)); };
60+
KernelD kernelD = std::bind(&Vacuum<DerivativeTraits, CollocationIntegrator>::kernelD, gf, _1, _2, _3);
61+
return integrator::doubleLayer(e, diagonal, kernelD);
5062
}
51-
52-
double computeS(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & gf, const Element & e) const {
53-
double epsInv = 1.0 / gf.epsilon();
54-
double area = e.area();
55-
return (factor_ * std::sqrt(4 * M_PI / area) * epsInv);
63+
/**@}*/
64+
65+
/**@{ Single and double layer potentials for a UniformDielectric Green's function by collocation */
66+
template <typename DerivativeTraits>
67+
Eigen::MatrixXd singleLayer(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
68+
auto f = this->factor_;
69+
auto epsInv = 1.0 / gf.epsilon();
70+
auto diagonal = [f, epsInv] (double area) -> double { return (f * std::sqrt(4 * M_PI / area) * epsInv); };
71+
KernelS kernelS = std::bind(&UniformDielectric<DerivativeTraits, CollocationIntegrator>::kernelS, gf, _1, _2);
72+
return integrator::singleLayer(e, diagonal, kernelS);
5673
}
57-
double computeD(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & e) const {
58-
double area = e.area();
59-
double radius = e.sphere().radius();
60-
return (-factor_ * std::sqrt(M_PI/ area) * (1.0 / radius));
74+
template <typename DerivativeTraits>
75+
Eigen::MatrixXd doubleLayer(const UniformDielectric<DerivativeTraits, CollocationIntegrator> & gf, const std::vector<Element> & e) const {
76+
auto f = this->factor_;
77+
auto diagonal = [f] (double area, double radius) -> double { return (-f * std::sqrt(M_PI/ area) * (1.0 / radius)); };
78+
KernelD kernelD = std::bind(&UniformDielectric<DerivativeTraits, CollocationIntegrator>::kernelD, gf, _1, _2, _3);
79+
return integrator::doubleLayer(e, diagonal, kernelD);
6180
}
81+
/**@}*/
6282

63-
double computeS(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & /* e */) const {
64-
return 1.0;
83+
/**@{ Single and double layer potentials for a IonicLiquid Green's function by collocation */
84+
template <typename DerivativeTraits>
85+
Eigen::MatrixXd singleLayer(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
86+
throw std::runtime_error("Not implemented yet!");
6587
}
66-
double computeD(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & /* e */) const {
67-
return 1.0;
88+
template <typename DerivativeTraits>
89+
Eigen::MatrixXd doubleLayer(const IonicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
90+
throw std::runtime_error("Not implemented yet!");
6891
}
92+
/**@}*/
6993

70-
double computeS(const AnisotropicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & /* e */) const {
71-
return 1.0;
94+
/**@{ Single and double layer potentials for an AnisotropicLiquid Green's function by collocation */
95+
template <typename DerivativeTraits>
96+
Eigen::MatrixXd singleLayer(const AnisotropicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
97+
throw std::runtime_error("Not implemented yet!");
7298
}
73-
double computeD(const AnisotropicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const Element & /* e */) const {
74-
return 1.0;
99+
template <typename DerivativeTraits>
100+
Eigen::MatrixXd doubleLayer(const AnisotropicLiquid<DerivativeTraits, CollocationIntegrator> & /* gf */, const std::vector<Element> & /* e */) const {
101+
throw std::runtime_error("Not implemented yet!");
75102
}
76-
77-
double computeS(const SphericalDiffuse<CollocationIntegrator, ProfilePolicy> & gf, const Element & e) const {
78-
// The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
79-
double area = e.area();
80-
// Diagonal of S inside the cavity
81-
double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
82-
// "Diagonal" of Coulomb singularity separation coefficient
83-
double coulomb_coeff = gf.coefficientCoulomb(e.center(), e.center());
84-
// "Diagonal" of the image Green's function
85-
double image = gf.imagePotential(e.center(), e.center());
86-
87-
return (Sii_I / coulomb_coeff + image);
103+
/**@}*/
104+
105+
/**@{ Single and double layer potentials for a SphericalDiffuse Green's function by collocation */
106+
template <typename ProfilePolicy>
107+
Eigen::MatrixXd singleLayer(const SphericalDiffuse<CollocationIntegrator, ProfilePolicy> & /* gf */, const std::vector<Element> & e) const {
108+
// // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
109+
// double area = e.area();
110+
// // Diagonal of S inside the cavity
111+
// double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
112+
// // "Diagonal" of Coulomb singularity separation coefficient
113+
// double coulomb_coeff = gf.coefficientCoulomb(e.center(), e.center());
114+
// // "Diagonal" of the image Green's function
115+
// double image = gf.imagePotential(e.center(), e.center());
116+
117+
// return (Sii_I / coulomb_coeff + image);
118+
return Eigen::MatrixXd::Zero(e.size(), e.size());
88119
}
89-
double computeD(const SphericalDiffuse<CollocationIntegrator, ProfilePolicy> & gf, const Element & e) const {
90-
// The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
91-
double area = e.area();
92-
double radius = e.sphere().radius();
93-
// Diagonal of S inside the cavity
94-
double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
95-
// Diagonal of D inside the cavity
96-
double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
97-
// "Diagonal" of Coulomb singularity separation coefficient
98-
double coulomb_coeff = gf.coefficientCoulomb(e.center(), e.center());
99-
// "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
100-
double coeff_grad = gf.coefficientCoulombDerivative(e.normal(), e.center(), e.center()) / std::pow(coulomb_coeff, 2);
101-
// "Diagonal" of the directional derivative of the image Green's function
102-
double image_grad = gf.imagePotentialDerivative(e.normal(), e.center(), e.center());
103-
104-
double eps_r2 = 0.0;
105-
std::tie(eps_r2, std::ignore) = gf.epsilon(e.center());
106-
107-
return eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
120+
template <typename ProfilePolicy>
121+
Eigen::MatrixXd doubleLayer(const SphericalDiffuse<CollocationIntegrator, ProfilePolicy> & /* gf */, const std::vector<Element> & e) const {
122+
// // The singular part is "integrated" as usual, while the nonsingular part is evaluated in full
123+
// double area = e.area();
124+
// double radius = e.sphere().radius();
125+
// // Diagonal of S inside the cavity
126+
// double Sii_I = factor_ * std::sqrt(4 * M_PI / area);
127+
// // Diagonal of D inside the cavity
128+
// double Dii_I = -factor_ * std::sqrt(M_PI/ area) * (1.0 / radius);
129+
// // "Diagonal" of Coulomb singularity separation coefficient
130+
// double coulomb_coeff = gf.coefficientCoulomb(e.center(), e.center());
131+
// // "Diagonal" of the directional derivative of the Coulomb singularity separation coefficient
132+
// double coeff_grad = gf.coefficientCoulombDerivative(e.normal(), e.center(), e.center()) / std::pow(coulomb_coeff, 2);
133+
// // "Diagonal" of the directional derivative of the image Green's function
134+
// double image_grad = gf.imagePotentialDerivative(e.normal(), e.center(), e.center());
135+
136+
// double eps_r2 = 0.0;
137+
// std::tie(eps_r2, std::ignore) = gf.epsilon(e.center());
138+
139+
// return eps_r2 * (Dii_I / coulomb_coeff - Sii_I * coeff_grad + image_grad);
140+
return Eigen::MatrixXd::Zero(e.size(), e.size());
108141
}
142+
/**@}*/
109143

110144
/// Scaling factor for the collocation formulas
111145
double factor_;

‎src/bi_operators/IntegratorForward.hpp‎

Lines changed: 31 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,31 @@
1+
/* pcmsolver_copyright_start */
2+
/*
3+
* PCMSolver, an API for the Polarizable Continuum Model
4+
* Copyright (C) 2013 Roberto Di Remigio, Luca Frediani and contributors
5+
*
6+
* This file is part of PCMSolver.
7+
*
8+
* PCMSolver is free software: you can redistribute it and/or modify
9+
* it under the terms of the GNU Lesser General Public License as published by
10+
* the Free Software Foundation, either version 3 of the License, or
11+
* (at your option) any later version.
12+
*
13+
* PCMSolver is distributed in the hope that it will be useful,
14+
* but WITHOUT ANY WARRANTY; without even the implied warranty of
15+
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16+
* GNU Lesser General Public License for more details.
17+
*
18+
* You should have received a copy of the GNU Lesser General Public License
19+
* along with PCMSolver. If not, see <http://www.gnu.org/licenses/>.
20+
*
21+
* For information on the complete list of contributors to the
22+
* PCMSolver API, see: <http://pcmsolver.github.io/pcmsolver-doc>
23+
*/
24+
/* pcmsolver_copyright_end */
25+
26+
#ifndef INTEGRATORFORWARD_HPP
27+
#define INTEGRATORFORWARD_HPP
28+
29+
struct CollocationIntegrator;
30+
31+
#endif // INTEGRATORFORWARD_HPP

‎src/solver/SolverHelperFunctions.hpp‎ renamed to ‎src/bi_operators/IntegratorHelperFunctions.hpp‎

Lines changed: 17 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -23,8 +23,8 @@
2323
*/
2424
/* pcmsolver_copyright_end */
2525

26-
#ifndef SOLVERHELPERFUNCTIONS_HPP
27-
#define SOLVERHELPERFUNCTIONS_HPP
26+
#ifndef INTEGRATORHELPERFUNCTIONS_HPP
27+
#define INTEGRATORHELPERFUNCTIONS_HPP
2828

2929
#include <cmath>
3030
#include <functional>
@@ -36,16 +36,22 @@
3636

3737
#include "Element.hpp"
3838

39-
/*! \typedef Diagonal
40-
* \brief functor handle to the diagonalS and diagonalD methods in IGreensFunction
39+
namespace integrator {
40+
/*! \typedef DiagonalS
41+
* \brief functor handle to the calculation of the diagonal of S
4142
*/
42-
typedef std::function<double(const Element &)> Diagonal;
43+
typedef std::function<double(double)> DiagonalS;
4344

4445
/*! \typedef KernelS
4546
* \brief functor handle to the kernelS method in IGreensFunction
4647
*/
4748
typedef std::function<double(const Eigen::Vector3d &, const Eigen::Vector3d &)> KernelS;
4849

50+
/*! \typedef DiagonalD
51+
* \brief functor handle to the calculation of the diagonal of D
52+
*/
53+
typedef std::function<double(double, double)> DiagonalD;
54+
4955
/*! \typedef KernelD
5056
* \brief functor handle to the kernelD method in IGreensFunction
5157
*/
@@ -57,13 +63,13 @@ typedef std::function<double(const Eigen::Vector3d &, const Eigen::Vector3d &, c
5763
* \param[in] kern function for the evaluation of the off-diagonal of S
5864
*/
5965
inline Eigen::MatrixXd singleLayer(const std::vector<Element> & elements,
60-
const Diagonal & diag, const KernelS & kern)
66+
const DiagonalS & diag, const KernelS & kern)
6167
{
6268
size_t mat_size = elements.size();
6369
Eigen::MatrixXd S = Eigen::MatrixXd::Zero(mat_size, mat_size);
6470
for (size_t i = 0; i < mat_size; ++i) {
6571
// Fill diagonal
66-
S(i, i) = diag(elements[i]);
72+
S(i, i) = diag(elements[i].area());
6773
Eigen::Vector3d source = elements[i].center();
6874
for (size_t j = 0; j < mat_size; ++j) {
6975
// Fill off-diagonal
@@ -80,13 +86,13 @@ inline Eigen::MatrixXd singleLayer(const std::vector<Element> & elements,
8086
* \param[in] kern function for the evaluation of the off-diagonal of D
8187
*/
8288
inline Eigen::MatrixXd doubleLayer(const std::vector<Element> & elements,
83-
const Diagonal & diag, const KernelD & kern)
89+
const DiagonalD & diag, const KernelD & kern)
8490
{
8591
size_t mat_size = elements.size();
8692
Eigen::MatrixXd D = Eigen::MatrixXd::Zero(mat_size, mat_size);
8793
for (size_t i = 0; i < mat_size; ++i) {
8894
// Fill diagonal
89-
D(i, i) = diag(elements[i]);
95+
D(i, i) = diag(elements[i].area(), elements[i].sphere().radius());
9096
Eigen::Vector3d source = elements[i].center();
9197
for (size_t j = 0; j < mat_size; ++j) {
9298
// Fill off-diagonal
@@ -98,5 +104,6 @@ inline Eigen::MatrixXd doubleLayer(const std::vector<Element> & elements,
98104
}
99105
return D;
100106
}
107+
} // namespace integrator
101108

102-
#endif // SOLVERHELPERFUNCTIONS_HPP
109+
#endif // INTEGRATORHELPERFUNCTIONS_HPP

‎src/bi_operators/IntegratorTypes.hpp‎

Lines changed: 1 addition & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -30,52 +30,8 @@
3030

3131
#include <boost/mpl/vector.hpp>
3232

33-
#include "DerivativeTypes.hpp"
3433
#include "CollocationIntegrator.hpp"
35-
#include "NumericalIntegrator.hpp"
3634

37-
typedef boost::mpl::vector<
38-
CollocationIntegrator<Numerical, Uniform>,
39-
CollocationIntegrator<AD_directional, Uniform>,
40-
CollocationIntegrator<AD_gradient, Uniform>,
41-
CollocationIntegrator<AD_hessian, Uniform>,
42-
NumericalIntegrator<Numerical, Uniform>,
43-
NumericalIntegrator<AD_directional, Uniform>,
44-
NumericalIntegrator<AD_gradient, Uniform>,
45-
NumericalIntegrator<AD_hessian, Uniform>
46-
> integrator_types_uniform;
47-
48-
typedef boost::mpl::vector<
49-
CollocationIntegrator<Numerical, Yukawa>,
50-
CollocationIntegrator<AD_directional, Yukawa>,
51-
CollocationIntegrator<AD_gradient, Yukawa>,
52-
CollocationIntegrator<AD_hessian, Yukawa>,
53-
NumericalIntegrator<Numerical, Yukawa>,
54-
NumericalIntegrator<AD_directional, Yukawa>,
55-
NumericalIntegrator<AD_gradient, Yukawa>,
56-
NumericalIntegrator<AD_hessian, Yukawa>
57-
> integrator_types_yukawa;
58-
59-
typedef boost::mpl::vector<
60-
CollocationIntegrator<Numerical, Anisotropic>,
61-
CollocationIntegrator<AD_directional, Anisotropic>,
62-
CollocationIntegrator<AD_gradient, Anisotropic>,
63-
CollocationIntegrator<AD_hessian, Anisotropic>,
64-
NumericalIntegrator<Numerical, Anisotropic>,
65-
NumericalIntegrator<AD_directional, Anisotropic>,
66-
NumericalIntegrator<AD_gradient, Anisotropic>,
67-
NumericalIntegrator<AD_hessian, Anisotropic>
68-
> integrator_types_anisotropic;
69-
70-
typedef boost::mpl::vector<
71-
CollocationIntegrator<Numerical, TanhDiffuse>,
72-
CollocationIntegrator<AD_directional, TanhDiffuse>,
73-
CollocationIntegrator<AD_gradient, TanhDiffuse>,
74-
CollocationIntegrator<AD_hessian, TanhDiffuse>,
75-
NumericalIntegrator<Numerical, TanhDiffuse>,
76-
NumericalIntegrator<AD_directional, TanhDiffuse>,
77-
NumericalIntegrator<AD_gradient, TanhDiffuse>,
78-
NumericalIntegrator<AD_hessian, TanhDiffuse>
79-
> integrator_types_tanhdiffuse;
35+
typedef boost::mpl::vector<CollocationIntegrator> integrator_types;
8036

8137
#endif // INTEGRATORTYPES_HPP

0 commit comments

Comments
 (0)