7 #ifndef GRIDFORCEGRID_INL
8 #define GRIDFORCEGRID_INL
19 int err =
get_inds(pos, inds, dg, gapscale);
24 DebugM(1,
"gapscale = " << gapscale <<
"\n");
25 DebugM(1,
"dg = " << dg <<
"\n");
26 DebugM(1,
"ind + dg = " << inds[0]+dg[0] <<
" " << inds[1]+dg[1] <<
" " << inds[2]+dg[2] <<
"\n");
42 for (
int j = 0; j < 64; j++)
DebugM(1,
"b[" << j <<
"] = " << b[j] <<
"\n" <<
endi);
47 for (
int j = 0; j < 64; j++)
DebugM(1,
"a[" << j <<
"] = " << a[j] <<
"\n" <<
endi);
51 float x[4],
y[4],
z[4];
52 x[0] = 1; y[0] = 1; z[0] = 1;
53 for (
int j = 1; j < 4; j++) {
90 wts[i][0] = -(1-dg.
y) * (1-dg.
z);
91 wts[i][1] = -(1-dg.
y) * dg.
z;
92 wts[i][2] = - dg.
y * (1-dg.
z);
93 wts[i][3] = - dg.
y * dg.
z;
94 for (
int j=0; j<4; j++) wts[i][j+4] = -wts[i][j];
97 wts[i][0] = -(1-dg.
x) * (1-dg.
z);
98 wts[i][1] = -(1-dg.
x) * dg.
z;
99 wts[i][2] = -wts[i][0];
100 wts[i][3] = -wts[i][1];
101 wts[i][4] = - dg.
x * (1-dg.
z);
102 wts[i][5] = - dg.
x * dg.
z;
103 wts[i][6] = -wts[i][4];
104 wts[i][7] = -wts[i][5];
107 wts[i][0] = - (1-dg.
x) * (1-dg.
y);
108 wts[i][1] = -wts[i][0];
109 wts[i][2] = - (1-dg.
x) * dg.
y ;
110 wts[i][3] = -wts[i][2];
111 wts[i][4] = - dg.
x * (1-dg.
y);
112 wts[i][5] = -wts[i][4];
113 wts[i][6] = - dg.
x * dg.
y ;
114 wts[i][7] = -wts[i][6];
117 for (
int j=0; j<4; j++) wts[i][j] = (1-dg.
x) * wts[i+1][j+4];
118 for (
int j=0; j<4; j++) wts[i][j+4] = dg.
x * wts[i+1][j+4];
120 for (i = 0; i < 4; i++) {
125 dV =
Vector(results[1], results[2], results[3]) *
inv;
138 for (
int i = 0; i < 3; i++) {
139 inds[i] = (int)floor(g[i]);
140 dg[i] = g[i] - inds[i];
143 for (
int i = 0; i < 3; i++) {
144 if (inds[i] < 0 || inds[i] >=
k[i]-1) {
145 if (
cont[i]) inds[i] =
k[i]-1;
148 if (
cont[i] && inds[i] ==
k[i]-1) {
151 if (g[i] < 0.0) dg[i] = 1.0 + g[i]*
gapinv[i];
152 else dg[i] = (g[i] - inds[i]) * gapinv[i];
164 for (
int l = 0; l < 4; l++) {
165 for (
int k = 0;
k < 4;
k++) {
166 for (
int j = 0; j < 4; j++) {
167 V += a[ind] * x[j] * y[
k] * z[l];
180 for (
int l = 0; l < 4; l++) {
181 for (
int k = 0;
k < 4;
k++) {
182 for (
int j = 0; j < 4; j++) {
183 if (j > 0) dV.
x += a[ind] * j * x[j-1] * y[
k] * z[l];
184 if (
k > 0) dV.
y += a[ind] *
k * x[j] * y[
k-1] * z[l];
185 if (l > 0) dV.
z += a[ind] * l * x[j] * y[
k] * z[l-1];
198 for (
int l = 0; l < 4; l++) {
199 for (
int k = 0;
k < 4;
k++) {
200 for (
int j = 0; j < 4; j++) {
201 if (j > 0 &&
k > 0) d2V.
x += a[ind] * j *
k * x[j-1] * y[
k-1] * z[l];
202 if (j > 0 && l > 0) d2V.
y += a[ind] * j * l * x[j-1] * y[
k] * z[l-1];
203 if (
k > 0 && l > 0) d2V.
z += a[ind] *
k * l * x[j] * y[
k-1] * z[l-1];
216 for (
int l = 0; l < 4; l++) {
217 for (
int k = 0;
k < 4;
k++) {
218 for (
int j = 0; j < 4; j++) {
219 if (j > 0 &&
k > 0 && l > 0) d3V += a[ind] * j *
k * l * x[j-1] * y[
k-1] * z[l-1];
233 a[2] = -3*b[0] + 3*b[1] - 2*b[8] - b[9];
234 a[3] = 2*b[0] - 2*b[1] + b[8] + b[9];
237 a[6] = -3*b[16] + 3*b[17] - 2*b[32] - b[33];
238 a[7] = 2*b[16] - 2*b[17] + b[32] + b[33];
239 a[8] = -3*b[0] + 3*b[2] - 2*b[16] - b[18];
240 a[9] = -3*b[8] + 3*b[10] - 2*b[32] - b[34];
241 a[10] = 9*b[0] - 9*b[1] - 9*b[2] + 9*b[3] + 6*b[8] + 3*b[9] - 6*b[10] - 3*b[11]
242 + 6*b[16] - 6*b[17] + 3*b[18] - 3*b[19] + 4*b[32] + 2*b[33] + 2*b[34] + b[35];
243 a[11] = -6*b[0] + 6*b[1] + 6*b[2] - 6*b[3] - 3*b[8] - 3*b[9] + 3*b[10] + 3*b[11]
244 - 4*b[16] + 4*b[17] - 2*b[18] + 2*b[19] - 2*b[32] - 2*b[33] - b[34] - b[35];
245 a[12] = 2*b[0] - 2*b[2] + b[16] + b[18];
246 a[13] = 2*b[8] - 2*b[10] + b[32] + b[34];
247 a[14] = -6*b[0] + 6*b[1] + 6*b[2] - 6*b[3] - 4*b[8] - 2*b[9] + 4*b[10] + 2*b[11]
248 - 3*b[16] + 3*b[17] - 3*b[18] + 3*b[19] - 2*b[32] - b[33] - 2*b[34] - b[35];
249 a[15] = 4*b[0] - 4*b[1] - 4*b[2] + 4*b[3] + 2*b[8] + 2*b[9] - 2*b[10] - 2*b[11]
250 + 2*b[16] - 2*b[17] + 2*b[18] - 2*b[19] + b[32] + b[33] + b[34] + b[35];
253 a[18] = -3*b[24] + 3*b[25] - 2*b[40] - b[41];
254 a[19] = 2*b[24] - 2*b[25] + b[40] + b[41];
257 a[22] = -3*b[48] + 3*b[49] - 2*b[56] - b[57];
258 a[23] = 2*b[48] - 2*b[49] + b[56] + b[57];
259 a[24] = -3*b[24] + 3*b[26] - 2*b[48] - b[50];
260 a[25] = -3*b[40] + 3*b[42] - 2*b[56] - b[58];
261 a[26] = 9*b[24] - 9*b[25] - 9*b[26] + 9*b[27] + 6*b[40] + 3*b[41] - 6*b[42] - 3*b[43]
262 + 6*b[48] - 6*b[49] + 3*b[50] - 3*b[51] + 4*b[56] + 2*b[57] + 2*b[58] + b[59];
263 a[27] = -6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 3*b[40] - 3*b[41] + 3*b[42] + 3*b[43]
264 - 4*b[48] + 4*b[49] - 2*b[50] + 2*b[51] - 2*b[56] - 2*b[57] - b[58] - b[59];
265 a[28] = 2*b[24] - 2*b[26] + b[48] + b[50];
266 a[29] = 2*b[40] - 2*b[42] + b[56] + b[58];
267 a[30] = -6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 4*b[40] - 2*b[41] + 4*b[42] + 2*b[43]
268 - 3*b[48] + 3*b[49] - 3*b[50] + 3*b[51] - 2*b[56] - b[57] - 2*b[58] - b[59];
269 a[31] = 4*b[24] - 4*b[25] - 4*b[26] + 4*b[27] + 2*b[40] + 2*b[41] - 2*b[42] - 2*b[43]
270 + 2*b[48] - 2*b[49] + 2*b[50] - 2*b[51] + b[56] + b[57] + b[58] + b[59];
271 a[32] = -3*b[0] + 3*b[4] - 2*b[24] - b[28];
272 a[33] = -3*b[8] + 3*b[12] - 2*b[40] - b[44];
273 a[34] = 9*b[0] - 9*b[1] - 9*b[4] + 9*b[5] + 6*b[8] + 3*b[9] - 6*b[12] - 3*b[13]
274 + 6*b[24] - 6*b[25] + 3*b[28] - 3*b[29] + 4*b[40] + 2*b[41] + 2*b[44] + b[45];
275 a[35] = -6*b[0] + 6*b[1] + 6*b[4] - 6*b[5] - 3*b[8] - 3*b[9] + 3*b[12] + 3*b[13]
276 - 4*b[24] + 4*b[25] - 2*b[28] + 2*b[29] - 2*b[40] - 2*b[41] - b[44] - b[45];
277 a[36] = -3*b[16] + 3*b[20] - 2*b[48] - b[52];
278 a[37] = -3*b[32] + 3*b[36] - 2*b[56] - b[60];
279 a[38] = 9*b[16] - 9*b[17] - 9*b[20] + 9*b[21] + 6*b[32] + 3*b[33] - 6*b[36] - 3*b[37]
280 + 6*b[48] - 6*b[49] + 3*b[52] - 3*b[53] + 4*b[56] + 2*b[57] + 2*b[60] + b[61];
281 a[39] = -6*b[16] + 6*b[17] + 6*b[20] - 6*b[21] - 3*b[32] - 3*b[33] + 3*b[36] + 3*b[37]
282 - 4*b[48] + 4*b[49] - 2*b[52] + 2*b[53] - 2*b[56] - 2*b[57] - b[60] - b[61];
283 a[40] = 9*b[0] - 9*b[2] - 9*b[4] + 9*b[6] + 6*b[16] + 3*b[18] - 6*b[20] - 3*b[22]
284 + 6*b[24] - 6*b[26] + 3*b[28] - 3*b[30] + 4*b[48] + 2*b[50] + 2*b[52] + b[54];
285 a[41] = 9*b[8] - 9*b[10] - 9*b[12] + 9*b[14] + 6*b[32] + 3*b[34] - 6*b[36] - 3*b[38]
286 + 6*b[40] - 6*b[42] + 3*b[44] - 3*b[46] + 4*b[56] + 2*b[58] + 2*b[60] + b[62];
287 a[42] = -27*b[0] + 27*b[1] + 27*b[2] - 27*b[3] + 27*b[4] - 27*b[5] - 27*b[6] + 27*b[7]
288 - 18*b[8] - 9*b[9] + 18*b[10] + 9*b[11] + 18*b[12] + 9*b[13] - 18*b[14] - 9*b[15]
289 - 18*b[16] + 18*b[17] - 9*b[18] + 9*b[19] + 18*b[20] - 18*b[21] + 9*b[22] - 9*b[23]
290 - 18*b[24] + 18*b[25] + 18*b[26] - 18*b[27] - 9*b[28] + 9*b[29] + 9*b[30] - 9*b[31]
291 - 12*b[32] - 6*b[33] - 6*b[34] - 3*b[35] + 12*b[36] + 6*b[37] + 6*b[38] + 3*b[39]
292 - 12*b[40] - 6*b[41] + 12*b[42] + 6*b[43] - 6*b[44] - 3*b[45] + 6*b[46] + 3*b[47]
293 - 12*b[48] + 12*b[49] - 6*b[50] + 6*b[51] - 6*b[52] + 6*b[53] - 3*b[54] + 3*b[55]
294 - 8*b[56] - 4*b[57] - 4*b[58] - 2*b[59] - 4*b[60] - 2*b[61] - 2*b[62] - b[63];
295 a[43] = 18*b[0] - 18*b[1] - 18*b[2] + 18*b[3] - 18*b[4] + 18*b[5] + 18*b[6] - 18*b[7]
296 + 9*b[8] + 9*b[9] - 9*b[10] - 9*b[11] - 9*b[12] - 9*b[13] + 9*b[14] + 9*b[15]
297 + 12*b[16] - 12*b[17] + 6*b[18] - 6*b[19] - 12*b[20] + 12*b[21] - 6*b[22] + 6*b[23]
298 + 12*b[24] - 12*b[25] - 12*b[26] + 12*b[27] + 6*b[28] - 6*b[29] - 6*b[30] + 6*b[31]
299 + 6*b[32] + 6*b[33] + 3*b[34] + 3*b[35] - 6*b[36] - 6*b[37] - 3*b[38] - 3*b[39]
300 + 6*b[40] + 6*b[41] - 6*b[42] - 6*b[43] + 3*b[44] + 3*b[45] - 3*b[46] - 3*b[47]
301 + 8*b[48] - 8*b[49] + 4*b[50] - 4*b[51] + 4*b[52] - 4*b[53] + 2*b[54] - 2*b[55]
302 + 4*b[56] + 4*b[57] + 2*b[58] + 2*b[59] + 2*b[60] + 2*b[61] + b[62] + b[63];
303 a[44] = -6*b[0] + 6*b[2] + 6*b[4] - 6*b[6] - 3*b[16] - 3*b[18] + 3*b[20] + 3*b[22]
304 - 4*b[24] + 4*b[26] - 2*b[28] + 2*b[30] - 2*b[48] - 2*b[50] - b[52] - b[54];
305 a[45] = -6*b[8] + 6*b[10] + 6*b[12] - 6*b[14] - 3*b[32] - 3*b[34] + 3*b[36] + 3*b[38]
306 - 4*b[40] + 4*b[42] - 2*b[44] + 2*b[46] - 2*b[56] - 2*b[58] - b[60] - b[62];
307 a[46] = 18*b[0] - 18*b[1] - 18*b[2] + 18*b[3] - 18*b[4] + 18*b[5] + 18*b[6] - 18*b[7]
308 + 12*b[8] + 6*b[9] - 12*b[10] - 6*b[11] - 12*b[12] - 6*b[13] + 12*b[14] + 6*b[15]
309 + 9*b[16] - 9*b[17] + 9*b[18] - 9*b[19] - 9*b[20] + 9*b[21] - 9*b[22] + 9*b[23]
310 + 12*b[24] - 12*b[25] - 12*b[26] + 12*b[27] + 6*b[28] - 6*b[29] - 6*b[30] + 6*b[31]
311 + 6*b[32] + 3*b[33] + 6*b[34] + 3*b[35] - 6*b[36] - 3*b[37] - 6*b[38] - 3*b[39]
312 + 8*b[40] + 4*b[41] - 8*b[42] - 4*b[43] + 4*b[44] + 2*b[45] - 4*b[46] - 2*b[47]
313 + 6*b[48] - 6*b[49] + 6*b[50] - 6*b[51] + 3*b[52] - 3*b[53] + 3*b[54] - 3*b[55]
314 + 4*b[56] + 2*b[57] + 4*b[58] + 2*b[59] + 2*b[60] + b[61] + 2*b[62] + b[63];
315 a[47] = -12*b[0] + 12*b[1] + 12*b[2] - 12*b[3] + 12*b[4] - 12*b[5] - 12*b[6] + 12*b[7]
316 - 6*b[8] - 6*b[9] + 6*b[10] + 6*b[11] + 6*b[12] + 6*b[13] - 6*b[14] - 6*b[15]
317 - 6*b[16] + 6*b[17] - 6*b[18] + 6*b[19] + 6*b[20] - 6*b[21] + 6*b[22] - 6*b[23]
318 - 8*b[24] + 8*b[25] + 8*b[26] - 8*b[27] - 4*b[28] + 4*b[29] + 4*b[30] - 4*b[31]
319 - 3*b[32] - 3*b[33] - 3*b[34] - 3*b[35] + 3*b[36] + 3*b[37] + 3*b[38] + 3*b[39]
320 - 4*b[40] - 4*b[41] + 4*b[42] + 4*b[43] - 2*b[44] - 2*b[45] + 2*b[46] + 2*b[47]
321 - 4*b[48] + 4*b[49] - 4*b[50] + 4*b[51] - 2*b[52] + 2*b[53] - 2*b[54] + 2*b[55]
322 - 2*b[56] - 2*b[57] - 2*b[58] - 2*b[59] - b[60] - b[61] - b[62] - b[63];
323 a[48] = 2*b[0] - 2*b[4] + b[24] + b[28];
324 a[49] = 2*b[8] - 2*b[12] + b[40] + b[44];
325 a[50] = -6*b[0] + 6*b[1] + 6*b[4] - 6*b[5] - 4*b[8] - 2*b[9] + 4*b[12] + 2*b[13]
326 - 3*b[24] + 3*b[25] - 3*b[28] + 3*b[29] - 2*b[40] - b[41] - 2*b[44] - b[45];
327 a[51] = 4*b[0] - 4*b[1] - 4*b[4] + 4*b[5] + 2*b[8] + 2*b[9] - 2*b[12] - 2*b[13]
328 + 2*b[24] - 2*b[25] + 2*b[28] - 2*b[29] + b[40] + b[41] + b[44] + b[45];
329 a[52] = 2*b[16] - 2*b[20] + b[48] + b[52];
330 a[53] = 2*b[32] - 2*b[36] + b[56] + b[60];
331 a[54] = -6*b[16] + 6*b[17] + 6*b[20] - 6*b[21] - 4*b[32] - 2*b[33] + 4*b[36] + 2*b[37]
332 - 3*b[48] + 3*b[49] - 3*b[52] + 3*b[53] - 2*b[56] - b[57] - 2*b[60] - b[61];
333 a[55] = 4*b[16] - 4*b[17] - 4*b[20] + 4*b[21] + 2*b[32] + 2*b[33] - 2*b[36] - 2*b[37]
334 + 2*b[48] - 2*b[49] + 2*b[52] - 2*b[53] + b[56] + b[57] + b[60] + b[61];
335 a[56] = -6*b[0] + 6*b[2] + 6*b[4] - 6*b[6] - 4*b[16] - 2*b[18] + 4*b[20] + 2*b[22]
336 - 3*b[24] + 3*b[26] - 3*b[28] + 3*b[30] - 2*b[48] - b[50] - 2*b[52] - b[54];
337 a[57] = -6*b[8] + 6*b[10] + 6*b[12] - 6*b[14] - 4*b[32] - 2*b[34] + 4*b[36] + 2*b[38]
338 - 3*b[40] + 3*b[42] - 3*b[44] + 3*b[46] - 2*b[56] - b[58] - 2*b[60] - b[62];
339 a[58] = 18*b[0] - 18*b[1] - 18*b[2] + 18*b[3] - 18*b[4] + 18*b[5] + 18*b[6] - 18*b[7]
340 + 12*b[8] + 6*b[9] - 12*b[10] - 6*b[11] - 12*b[12] - 6*b[13] + 12*b[14] + 6*b[15]
341 + 12*b[16] - 12*b[17] + 6*b[18] - 6*b[19] - 12*b[20] + 12*b[21] - 6*b[22] + 6*b[23]
342 + 9*b[24] - 9*b[25] - 9*b[26] + 9*b[27] + 9*b[28] - 9*b[29] - 9*b[30] + 9*b[31]
343 + 8*b[32] + 4*b[33] + 4*b[34] + 2*b[35] - 8*b[36] - 4*b[37] - 4*b[38] - 2*b[39]
344 + 6*b[40] + 3*b[41] - 6*b[42] - 3*b[43] + 6*b[44] + 3*b[45] - 6*b[46] - 3*b[47]
345 + 6*b[48] - 6*b[49] + 3*b[50] - 3*b[51] + 6*b[52] - 6*b[53] + 3*b[54] - 3*b[55]
346 + 4*b[56] + 2*b[57] + 2*b[58] + b[59] + 4*b[60] + 2*b[61] + 2*b[62] + b[63];
347 a[59] = -12*b[0] + 12*b[1] + 12*b[2] - 12*b[3] + 12*b[4] - 12*b[5] - 12*b[6] + 12*b[7]
348 - 6*b[8] - 6*b[9] + 6*b[10] + 6*b[11] + 6*b[12] + 6*b[13] - 6*b[14] - 6*b[15]
349 - 8*b[16] + 8*b[17] - 4*b[18] + 4*b[19] + 8*b[20] - 8*b[21] + 4*b[22] - 4*b[23]
350 - 6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 6*b[28] + 6*b[29] + 6*b[30] - 6*b[31]
351 - 4*b[32] - 4*b[33] - 2*b[34] - 2*b[35] + 4*b[36] + 4*b[37] + 2*b[38] + 2*b[39]
352 - 3*b[40] - 3*b[41] + 3*b[42] + 3*b[43] - 3*b[44] - 3*b[45] + 3*b[46] + 3*b[47]
353 - 4*b[48] + 4*b[49] - 2*b[50] + 2*b[51] - 4*b[52] + 4*b[53] - 2*b[54] + 2*b[55]
354 - 2*b[56] - 2*b[57] - b[58] - b[59] - 2*b[60] - 2*b[61] - b[62] - b[63];
355 a[60] = 4*b[0] - 4*b[2] - 4*b[4] + 4*b[6] + 2*b[16] + 2*b[18] - 2*b[20] - 2*b[22]
356 + 2*b[24] - 2*b[26] + 2*b[28] - 2*b[30] + b[48] + b[50] + b[52] + b[54];
357 a[61] = 4*b[8] - 4*b[10] - 4*b[12] + 4*b[14] + 2*b[32] + 2*b[34] - 2*b[36] - 2*b[38]
358 + 2*b[40] - 2*b[42] + 2*b[44] - 2*b[46] + b[56] + b[58] + b[60] + b[62];
359 a[62] = -12*b[0] + 12*b[1] + 12*b[2] - 12*b[3] + 12*b[4] - 12*b[5] - 12*b[6] + 12*b[7]
360 - 8*b[8] - 4*b[9] + 8*b[10] + 4*b[11] + 8*b[12] + 4*b[13] - 8*b[14] - 4*b[15]
361 - 6*b[16] + 6*b[17] - 6*b[18] + 6*b[19] + 6*b[20] - 6*b[21] + 6*b[22] - 6*b[23]
362 - 6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 6*b[28] + 6*b[29] + 6*b[30] - 6*b[31]
363 - 4*b[32] - 2*b[33] - 4*b[34] - 2*b[35] + 4*b[36] + 2*b[37] + 4*b[38] + 2*b[39]
364 - 4*b[40] - 2*b[41] + 4*b[42] + 2*b[43] - 4*b[44] - 2*b[45] + 4*b[46] + 2*b[47]
365 - 3*b[48] + 3*b[49] - 3*b[50] + 3*b[51] - 3*b[52] + 3*b[53] - 3*b[54] + 3*b[55]
366 - 2*b[56] - b[57] - 2*b[58] - b[59] - 2*b[60] - b[61] - 2*b[62] - b[63];
367 a[63] = 8*b[0] - 8*b[1] - 8*b[2] + 8*b[3] - 8*b[4] + 8*b[5] + 8*b[6] - 8*b[7]
368 + 4*b[8] + 4*b[9] - 4*b[10] - 4*b[11] - 4*b[12] - 4*b[13] + 4*b[14] + 4*b[15]
369 + 4*b[16] - 4*b[17] + 4*b[18] - 4*b[19] - 4*b[20] + 4*b[21] - 4*b[22] + 4*b[23]
370 + 4*b[24] - 4*b[25] - 4*b[26] + 4*b[27] + 4*b[28] - 4*b[29] - 4*b[30] + 4*b[31]
371 + 2*b[32] + 2*b[33] + 2*b[34] + 2*b[35] - 2*b[36] - 2*b[37] - 2*b[38] - 2*b[39]
372 + 2*b[40] + 2*b[41] - 2*b[42] - 2*b[43] + 2*b[44] + 2*b[45] - 2*b[46] - 2*b[47]
373 + 2*b[48] - 2*b[49] + 2*b[50] - 2*b[51] + 2*b[52] - 2*b[53] + 2*b[54] - 2*b[55]
374 + b[56] + b[57] + b[58] + b[59] + b[60] + b[61] + b[62] + b[63];
385 for (
int i = 0; i < 3; i++) {
386 inds[i] = (int)floor(g[i]);
387 dg[i] = g[i] - inds[i];
390 for (
int i = 0; i < 3; i++) {
391 if (inds[i] < 0 || inds[i] >=
k[i]-1) {
402 wts[0] = (1-dg.
x) * (1-dg.
y) * (1-dg.
z);
403 wts[1] = (1-dg.
x) * (1-dg.
y) * dg.
z;
404 wts[2] = (1-dg.
x) * dg.
y * (1-dg.
z);
405 wts[3] = (1-dg.
x) * dg.
y * dg.
z;
406 wts[4] = dg.
x * (1-dg.
y) * (1-dg.
z);
407 wts[5] = dg.
x * (1-dg.
y) * dg.
z;
408 wts[6] = dg.
x * dg.
y * (1-dg.
z);
409 wts[7] = dg.
x * dg.
y * dg.
z;
419 vals[1] =
get_grid(i0, i1, i2+1, i3);
420 vals[2] =
get_grid(i0, i1+1, i2, i3);
421 vals[3] =
get_grid(i0, i1+1, i2+1, i3);
422 vals[4] =
get_grid(i0+1, i1, i2, i3);
423 vals[5] =
get_grid(i0+1, i1, i2+1, i3);
424 vals[6] =
get_grid(i0+1, i1+1, i2, i3);
425 vals[7] =
get_grid(i0+1, i1+1, i2+1, i3);
442 for (
int i = 0; i < 8; i++) {
443 DebugM(2,
"vals[" << i <<
"] = " << vals[i] <<
" wts[" << i <<
"] = " << wts[i] <<
"\n" <<
endi);
449 wts[1] *
get_grid(i0, i1, i2+1, i3) +
450 wts[2] *
get_grid(i0, i1+1, i2, i3) +
451 wts[3] *
get_grid(i0, i1+1, i2+1, i3) +
452 wts[4] *
get_grid(i0+1, i1, i2, i3) +
453 wts[5] *
get_grid(i0+1, i1, i2+1, i3) +
454 wts[6] *
get_grid(i0+1, i1+1, i2, i3) +
455 wts[7] *
get_grid(i0+1, i1+1, i2+1, i3);
457 DebugM(2,
"result = " << result <<
"\n" <<
endi);
int compute_VdV(Position pos, float &V, Vector &dV) const
float compute_V(float *a, float *x, float *y, float *z) const
int get_inds(Position pos, int *inds, Vector &dg) const
static Tensor diagonal(const Vector &v1)
float linear_interpolate(int i0, int i1, int i2, int i3, const float *wts) const
virtual void compute_b(float *b, int *inds, Vector gapscale) const =0
Vector wrap_delta(const Position &pos1) const
Vector compute_dV(float *a, float *x, float *y, float *z) const
int compute_VdV(Position pos, float &V, Vector &dV) const
std::ostream & endi(std::ostream &s)
void compute_a(float *a, float *b) const
void compute_wts(float *wts, const Vector &dg) const
float get_grid(int i0, int i1, int i2, int i3) const
float compute_d3V(float *a, float *x, float *y, float *z) const
GridforceFullSubGrid ** subgrids
virtual Position get_center(void) const =0
int get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const
Position wrap_position(const Position &pos, const Lattice &lattice)
Vector compute_d2V(float *a, float *x, float *y, float *z) const