Skip to content

Commit 90dd095

Browse files
author
Roberto Di Remigio
committed
Numerical integration (from Gamess written by B. Mennucci) code compiles. NOT TESTED!
1 parent ceffe63 commit 90dd095

File tree

2 files changed

+292
-3
lines changed

2 files changed

+292
-3
lines changed

‎src/operators_diagonal/NumericalIntegrator.cpp‎

Lines changed: 176 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -25,17 +25,20 @@
2525

2626
#include "NumericalIntegrator.hpp"
2727

28+
#include <algorithm>
2829
#include <cmath>
2930
#include <stdexcept>
3031

3132
#include "Config.hpp"
3233

3334
#include "EigenPimpl.hpp"
3435

35-
#include <boost/numeric/quadrature/kronrodgauss.hpp>
36-
#include <boost/numeric/quadrature/error_estimator.hpp>
36+
//#include <boost/numeric/quadrature/kronrodgauss.hpp>
37+
//#include <boost/numeric/quadrature/error_estimator.hpp>
3738

3839
#include "Element.hpp"
40+
#include "MathUtils.hpp"
41+
#include "Sphere.hpp"
3942

4043
double NumericalIntegrator::computeS(const Vacuum<double> * gf, const Element & e) const {
4144
return 0.0;
@@ -116,7 +119,177 @@ double NumericalIntegrator::computeD(const IonicLiquid<AD_hessian> * gf, const E
116119
}
117120

118121
double NumericalIntegrator::computeS(const AnisotropicLiquid<double> * gf, const Element & e) const {
119-
return 0.0;
122+
double Sii = 0.0, S64 = 0.0;
123+
// Extract relevant data from Element
124+
int nVertices = e.nVertices();
125+
double area = e.area();
126+
Eigen::Vector3d center = e.center();
127+
Eigen::Vector3d normal = e.normal();
128+
Sphere sph = e.sphere();
129+
Eigen::Matrix3Xd vertices = e.vertices();
130+
Eigen::Matrix3Xd arcs = e.arcs();
131+
132+
double xx = 0.0,
133+
xy = 0.0,
134+
xz = normal(0),
135+
yy = 0.0,
136+
yx = 0.0,
137+
yz = normal(1),
138+
zz = normal(2),
139+
zx = 0.0,
140+
zy = 0.0;
141+
142+
double rmin = 0.99; // Some kind of threshold
143+
if (std::abs(xz) <= rmin) {
144+
rmin = std::abs(xz);
145+
xx = 0.0;
146+
yx = -zz / std::sqrt(1.0 - std::pow(xz, 2));
147+
zx = yz / std::sqrt(1.0 - std::pow(xz, 2));
148+
}
149+
if (std::abs(yz) <= rmin) {
150+
rmin = std::abs(yz);
151+
xx = yz / std::sqrt(1.0 - std::pow(yz, 2));
152+
yx = 0.0;
153+
zx = -xz / std::sqrt(1.0 - std::pow(yz, 2));
154+
}
155+
if (std::abs(zz) <= rmin) {
156+
rmin = std::abs(yz);
157+
xx = yz / std::sqrt(1.0 - std::pow(zz, 2));
158+
yx = -xz / std::sqrt(1.0 - std::pow(zz, 2));
159+
zx = 0.0;
160+
}
161+
Eigen::Vector3d tangent;
162+
tangent << xx, yx, zx; // Not sure this is the tangent vector...
163+
164+
/*xy = yz * zx - yx * zz;
165+
yy = zz * xx - zx * xz;
166+
zy = xz * yx - xx * yz;
167+
Eigen::Vector3d bitangent << xy, yy, zy; // Not sure this is the bitangent vector...*/
168+
Eigen::Vector3d bitangent = normal.cross(tangent);
169+
170+
std::vector<double> theta(nVertices), phi(nVertices), phinumb(nVertices+1);
171+
std::vector<int> numb(nVertices+1);
172+
// Clean-up heap crap
173+
std::fill_n(theta.begin(), nVertices, 0.0);
174+
std::fill_n(phi.begin(), nVertices, 0.0);
175+
std::fill_n(numb.begin(), nVertices+1, 0);
176+
std::fill_n(phinumb.begin(), nVertices+1, 0.0);
177+
// Calculate a number of angles
178+
for (int i = 0; i < nVertices; ++i) {
179+
Eigen::Vector3d vertex_normal = (vertices.col(i) - sph.center()) / sph.radius();
180+
double scal1 = vertex_normal.dot(normal);
181+
if (scal1 >= 1.0) scal1 = 1.0;
182+
if (scal1 <= -1.0) scal1 = -1.0;
183+
theta[i] = std::acos(scal1);
184+
double scal2 = vertex_normal.dot(tangent) / std::sin(theta[i]);
185+
if (scal2 >= 1.0) scal2 = 1.0;
186+
if (scal2 <= -1.0) scal2 = -1.0;
187+
phi[i] = std::acos(scal2);
188+
double sin_phi = vertex_normal.dot(bitangent) / std::sin(theta[i]);
189+
if (sin_phi <= 0.0) phi[i] = 2 * M_PI - phi[i];
190+
if (i != 0) {
191+
phi[i] -= phi[0];
192+
if (phi[i] < 0.0) phi[i] += 2 * M_PI;
193+
}
194+
}
195+
// Recalculate tangent and bitangent vectors
196+
tangent *= std::cos(phi[0]);
197+
tangent += bitangent * std::sin(phi[0]);
198+
bitangent = normal.cross(tangent);
199+
200+
// Populate numb and phinumb arrays
201+
phi[0] = 0.0;
202+
numb[0] = 1; numb[1] = 2;
203+
phinumb[0] = phi[0]; phinumb[1] = phinumb[1];
204+
for (int i = 2; i < nVertices; ++i) { // This loop is 2-based
205+
for (int j = 1; j < (i - 1); ++j) { // This loop is 1-based
206+
if (phi[i] < phinumb[j]) {
207+
for (int k = 0; k < (i - j); ++k) {
208+
numb[i - k + 1] = numb[i - k];
209+
phinumb[i - k + 1] = phinumb[i - k];
210+
}
211+
numb[j] = i;
212+
phinumb[j] = phi[i];
213+
goto jump; // Ugly, to be refactored!!!
214+
}
215+
numb[i] = i;
216+
phinumb[i] = phi[i];
217+
}
218+
jump:
219+
std::cout << "jumped" << std::endl;
220+
}
221+
numb[nVertices] = numb[0];
222+
phinumb[nVertices] = 2 * M_PI;
223+
224+
for (int i = 0; i < nVertices; ++i) { // Loop on edges
225+
double pha = phinumb[i];
226+
double phb = phinumb[i+1];
227+
double aph = (pha - phb) / 2.0;
228+
double bph = (pha + phb) / 2.0;
229+
double tha = theta[numb[i]];
230+
double thb = theta[numb[i+1]];
231+
double thmax = 0.0;
232+
Eigen::Vector3d oc = (arcs.col(i) - sph.center()) / sph.radius();
233+
double oc_norm = oc.norm();
234+
double oc_norm2 = std::pow(oc_norm, 2);
235+
for (int j = 0; j < 32; ++j) { // Loop on Gaussian points (64-points rule)
236+
for (int k = 0; k <= 1; ++k) {
237+
double ph = (2*k - 1) * aph * gauss64Abscissa(j) + bph;
238+
double cos_ph = std::cos(ph);
239+
double sin_ph = std::sin(ph);
240+
if (oc_norm2 < 1.0e-07) {
241+
double cotg_thmax = std::sin(ph-pha) / std::tan(thb) + std::sin(phb-ph) / std::tan(tha);
242+
cotg_thmax /= std::sin(phb-pha);
243+
thmax = std::atan(1.0 / cotg_thmax);
244+
} else {
245+
Eigen::Vector3d scratch;
246+
scratch << tangent.dot(oc), bitangent.dot(oc), normal.dot(oc);
247+
double aa = std::pow(tangent.dot(oc)*cos_ph + bitangent.dot(oc)*sin_ph, 2) + std::pow(normal.dot(oc), 2);
248+
double bb = -normal.dot(oc) * oc_norm2;
249+
double cc = std::pow(oc_norm2, 2) - std::pow(tangent.dot(oc)*cos_ph + bitangent.dot(oc)*sin_ph, 2);
250+
double ds = std::pow(bb, 2) - aa*cc;
251+
if (ds < 0.0) ds = 0.0;
252+
double cs = (-bb + std::sqrt(ds)) / aa;
253+
if (cs > 1.0) cs = 1.0;
254+
if (cs < -1.0) cs = 1.0;
255+
thmax = std::acos(cs);
256+
}
257+
double ath = thmax / 2.0;
258+
259+
double S16 = 0.0;
260+
if (thmax < 1.0e-08) goto jump1; // Ugly, to be refactored!!!
261+
for (int l = 0; l < 8; ++l) { // Loop on Gaussian points (16-points rule)
262+
for (int m = 0; m <= 1; ++m) {
263+
double th = (2*m - 1) * ath * gauss16Abscissa(l) + ath;
264+
double cos_th = std::cos(th);
265+
double sin_th = std::sin(th);
266+
double vx = tangent(0) * sin_th * cos_ph
267+
+ bitangent(0) * sin_th * sin_ph
268+
+ normal(0) * (cos_th - 1.0);
269+
double vy = tangent(1) * sin_th * cos_ph
270+
+ bitangent(1) * sin_th * sin_ph
271+
+ normal(1) * (cos_th - 1.0);
272+
double vz = tangent(2) * sin_th * cos_ph
273+
+ bitangent(2) * sin_th * sin_ph
274+
+ normal(2) * (cos_th - 1.0);
275+
Eigen::Vector3d evaluation;
276+
evaluation << vx, vy, vz;
277+
double rth = std::sqrt(2 * (1.0 - cos_th));
278+
double rtheps = gf->function(evaluation, Eigen::Vector3d::Zero()); // Evaluate Green's function at Gaussian points
279+
S16 += sph.radius() * rtheps * sin_th * ath * gauss16Weight(l);
280+
}
281+
}
282+
S64 += S16 * aph * gauss64Weight(j);
283+
284+
jump1:
285+
std::cout << "jumped!" << std::endl;
286+
}
287+
}
288+
}
289+
290+
Sii = S64 * area;
291+
292+
return Sii;
120293
}
121294
double NumericalIntegrator::computeS(const AnisotropicLiquid<AD_directional> * gf, const Element & e) const {
122295
return 0.0;

‎src/utils/MathUtils.hpp‎

Lines changed: 116 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -194,4 +194,120 @@ inline void eulerRotation(Eigen::Matrix3d & R_, const Eigen::Vector3d & eulerAng
194194
* Eigen::AngleAxis<double>(phi, Eigen::Vector3d::UnitZ());
195195
}
196196

197+
/*! Abscissae for 16-point Gaussian quadrature rule */
198+
inline double gauss16Abscissa(int i)
199+
{
200+
static std::vector<double> x16(8);
201+
202+
x16[0] = 0.9894009349916499325961542;
203+
x16[1] = 0.9445750230732325760779884;
204+
x16[2] = 0.8656312023878317438804679;
205+
x16[3] = 0.7554044083550030338951012;
206+
x16[4] = 0.6178762444026437484466718;
207+
x16[5] = 0.4580167776572273863424194;
208+
x16[6] = 0.2816035507792589132304605;
209+
x16[7] = 0.0950125098376374401853193;
210+
211+
return x16[i];
212+
}
213+
214+
/*! Weights for 16-point Gaussian quadrature rule */
215+
inline double gauss16Weight(int i)
216+
{
217+
static std::vector<double> w16(8);
218+
219+
w16[0] = 0.0271524594117540948517806;
220+
w16[1] = 0.0622535239386478928628438;
221+
w16[2] = 0.0951585116824927848099251;
222+
w16[3] = 0.1246289712555338720524763;
223+
w16[4] = 0.1495959888165767320815017;
224+
w16[5] = 0.1691565193950025381893121;
225+
w16[6] = 0.1826034150449235888667637;
226+
w16[7] = 0.1894506104550684962853967;
227+
228+
return w16[i];
229+
}
230+
231+
/*! Abscissae for 64-point Gaussian quadrature rule */
232+
inline double gauss64Abscissa(int i)
233+
{
234+
static std::vector<double> x64(32);
235+
236+
x64[0] = 0.9993050417357721394569056;
237+
x64[1] = 0.9963401167719552793469245;
238+
x64[2] = 0.9910133714767443207393824;
239+
x64[3] = 0.9833362538846259569312993;
240+
x64[4] = 0.9733268277899109637418535;
241+
x64[5] = 0.9610087996520537189186141;
242+
x64[6] = 0.9464113748584028160624815;
243+
x64[7] = 0.9295691721319395758214902;
244+
x64[8] = 0.9105221370785028057563807;
245+
x64[9] = 0.8893154459951141058534040;
246+
x64[10] = 0.8659993981540928197607834;
247+
x64[11] = 0.8406292962525803627516915;
248+
x64[12] = 0.8132653151227975597419233;
249+
x64[13] = 0.7839723589433414076102205;
250+
x64[14] = 0.7528199072605318966118638;
251+
x64[15] = 0.7198818501716108268489402;
252+
x64[16] = 0.6852363130542332425635584;
253+
x64[17] = 0.6489654712546573398577612;
254+
x64[18] = 0.6111553551723932502488530;
255+
x64[19] = 0.5718956462026340342838781;
256+
x64[20] = 0.5312794640198945456580139;
257+
x64[21] = 0.4894031457070529574785263;
258+
x64[22] = 0.4463660172534640879849477;
259+
x64[23] = 0.4022701579639916036957668;
260+
x64[24] = 0.3572201583376681159504426;
261+
x64[25] = 0.3113228719902109561575127;
262+
x64[26] = 0.2646871622087674163739642;
263+
x64[27] = 0.2174236437400070841496487;
264+
x64[28] = 0.1696444204239928180373136;
265+
x64[29] = 0.1214628192961205544703765;
266+
x64[30] = 0.0729931217877990394495429;
267+
x64[31] = 0.0243502926634244325089558;
268+
269+
return x64[i];
270+
}
271+
272+
/*! Weights for 64-point Gaussian quadrature rule */
273+
inline double gauss64Weight(int i)
274+
{
275+
static std::vector<double> w64(32);
276+
277+
w64[0] = 0.0017832807216964329472961;
278+
w64[1] = 0.0041470332605624676352875;
279+
w64[2] = 0.0065044579689783628561174;
280+
w64[3] = 0.0088467598263639477230309;
281+
w64[4] = 0.0111681394601311288185905;
282+
w64[5] = 0.0134630478967186425980608;
283+
w64[6] = 0.0157260304760247193219660;
284+
w64[7] = 0.0179517157756973430850453;
285+
w64[8] = 0.0201348231535302093723403;
286+
w64[9] = 0.0222701738083832541592983;
287+
w64[10] = 0.0243527025687108733381776;
288+
w64[11] = 0.0263774697150546586716918;
289+
w64[12] = 0.0283396726142594832275113;
290+
w64[13] = 0.0302346570724024788679741;
291+
w64[14] = 0.0320579283548515535854675;
292+
w64[15] = 0.0338051618371416093915655;
293+
w64[16] = 0.0354722132568823838106931;
294+
w64[17] = 0.0370551285402400460404151;
295+
w64[18] = 0.0385501531786156291289625;
296+
w64[19] = 0.0399537411327203413866569;
297+
w64[20] = 0.0412625632426235286101563;
298+
w64[21] = 0.0424735151236535890073398;
299+
w64[22] = 0.0435837245293234533768279;
300+
w64[23] = 0.0445905581637565630601347;
301+
w64[24] = 0.0454916279274181444797710;
302+
w64[25] = 0.0462847965813144172959532;
303+
w64[26] = 0.0469681828162100173253263;
304+
w64[27] = 0.0475401657148303086622822;
305+
w64[28] = 0.0479993885964583077281262;
306+
w64[29] = 0.0483447622348029571697695;
307+
w64[30] = 0.0485754674415034269347991;
308+
w64[31] = 0.0486909570091397203833654;
309+
310+
return w64[i];
311+
}
312+
197313
#endif // MATHUTILS_HPP

0 commit comments

Comments
 (0)