GridForceGrid.inl

Go to the documentation of this file.
00001 
00007 #ifndef GRIDFORCEGRID_INL
00008 #define GRIDFORCEGRID_INL
00009 
00010 #include "GridForceGrid.h"
00011 
00012 inline int GridforceFullBaseGrid::compute_VdV(Position pos, float &V, Vector &dV) const
00013 {
00014     //SimParameters *simParams = Node::Object()->simParameters;
00015     int inds[3];
00016     Vector g, dg;
00017     Vector gapscale = Vector(1, 1, 1);
00018     
00019     int err = get_inds(pos, inds, dg, gapscale);
00020     if (err) {
00021         return -1;
00022     }
00023     
00024     DebugM(1, "gapscale = " << gapscale << "\n");
00025     DebugM(1, "dg = " << dg << "\n");
00026     DebugM(1, "ind + dg = " << inds[0]+dg[0] << " " << inds[1]+dg[1] << " " << inds[2]+dg[2] << "\n");
00027     DebugM(3, "compute_VdV: generation = " << generation << "\n" << endi);
00028     
00029     // Pass to subgrid if one exists here
00030     for (int i = 0; i < numSubgrids; i++) {
00031         if (((inds[0] >= subgrids[i]->pmin[0] && inds[0] <= subgrids[i]->pmax[0]) || subgrids[i]->cont[0]) &&
00032             ((inds[1] >= subgrids[i]->pmin[1] && inds[1] <= subgrids[i]->pmax[1]) || subgrids[i]->cont[1]) &&
00033             ((inds[2] >= subgrids[i]->pmin[2] && inds[2] <= subgrids[i]->pmax[2]) || subgrids[i]->cont[2]))
00034         {
00035             return subgrids[i]->compute_VdV(pos, V, dV);
00036         }
00037     }
00038     
00039     // Compute b
00040     float b[64];        // Matrix of values at 8 box corners
00041     compute_b(b, inds, gapscale);
00042     for (int j = 0; j < 64; j++) DebugM(1, "b[" << j << "] = " << b[j] << "\n" << endi);
00043     
00044     // Compute a
00045     float a[64];
00046     compute_a(a, b);
00047     for (int j = 0; j < 64; j++) DebugM(1, "a[" << j << "] = " << a[j] << "\n" << endi);
00048             
00049     // Calculate powers of x, y, z for later use
00050     // e.g. x[2] = x^2
00051     float x[4], y[4], z[4];
00052     x[0] = 1; y[0] = 1; z[0] = 1;
00053     for (int j = 1; j < 4; j++) {
00054         x[j] = x[j-1] * dg.x;
00055         y[j] = y[j-1] * dg.y;
00056         z[j] = z[j-1] * dg.z;
00057     }
00058     
00059     V = compute_V(a, x, y, z);
00060     dV = Tensor::diagonal(gapscale) * (compute_dV(a, x, y, z) * inv);
00061     
00062     return 0;
00063 }
00064 
00065 
00066 inline int GridforceLiteGrid::compute_VdV(Position pos, float &V, Vector &dV) const
00067 {
00068     int inds[3];
00069     Vector g, dg;
00070     
00071     int err = get_inds(pos, inds, dg);
00072     if (err) {
00073         return -1;
00074     }
00075     
00076     float wts[4][8];
00077     float results[4];
00078     
00079     // compute_wts(wts, dg);
00080     // wts[0][0] = (1-dg.x) * (1-dg.y) * (1-dg.z);
00081     // wts[0][1] = (1-dg.x) * (1-dg.y) *   dg.z;
00082     // wts[0][2] = (1-dg.x) *   dg.y   * (1-dg.z);
00083     // wts[0][3] = (1-dg.x) *   dg.y   *   dg.z;
00084     // wts[0][4] =   dg.x   * (1-dg.y) * (1-dg.z);
00085     // wts[0][5] =   dg.x   * (1-dg.y) *   dg.z;
00086     // wts[0][6] =   dg.x   *   dg.y   * (1-dg.z);
00087     // wts[0][7] =   dg.x   *   dg.y   *   dg.z;
00088 
00089     int i = 1;
00090     wts[i][0] = -(1-dg.y) * (1-dg.z);
00091     wts[i][1] = -(1-dg.y) *   dg.z;
00092     wts[i][2] = -  dg.y   * (1-dg.z);
00093     wts[i][3] = -  dg.y   *   dg.z;
00094     for (int j=0; j<4; j++) wts[i][j+4] = -wts[i][j];
00095 
00096     i = 2;
00097     wts[i][0] = -(1-dg.x) * (1-dg.z);
00098     wts[i][1] = -(1-dg.x) *   dg.z;
00099     wts[i][2] = -wts[i][0];
00100     wts[i][3] = -wts[i][1];
00101     wts[i][4] =   - dg.x  * (1-dg.z);
00102     wts[i][5] =   - dg.x  *   dg.z;
00103     wts[i][6] = -wts[i][4];
00104     wts[i][7] = -wts[i][5];
00105 
00106     i = 3;
00107     wts[i][0] = - (1-dg.x) * (1-dg.y);
00108     wts[i][1] = -wts[i][0];
00109     wts[i][2] = - (1-dg.x) *   dg.y  ;
00110     wts[i][3] = -wts[i][2];
00111     wts[i][4] = - dg.x     * (1-dg.y);
00112     wts[i][5] = -wts[i][4];
00113     wts[i][6] = - dg.x     *   dg.y  ;
00114     wts[i][7] = -wts[i][6];
00115 
00116     i = 0;
00117     for (int j=0; j<4; j++) wts[i][j]   = (1-dg.x) * wts[i+1][j+4];
00118     for (int j=0; j<4; j++) wts[i][j+4] =   dg.x   * wts[i+1][j+4];    
00119 
00120     for (i = 0; i < 4; i++) {
00121         results[i] = linear_interpolate(inds[0], inds[1], inds[2], 0, wts[i]);
00122     }
00123     
00124     V = results[0];
00125     dV = Vector(results[1], results[2], results[3]) * inv;
00126     
00127     return 0;
00128 }
00129 
00130 
00131 inline int GridforceFullBaseGrid::get_inds(Position pos, int *inds, Vector &dg, Vector &gapscale) const
00132 {
00133     Vector p = pos - origin;
00134     Vector g;
00135     
00136     g = inv * p;
00137     
00138     for (int i = 0; i < 3; i++) {
00139         inds[i] = (int)floor(g[i]);
00140         dg[i] = g[i] - inds[i];
00141     }
00142     
00143     for (int i = 0; i < 3; i++) {
00144         if (inds[i] < 0 || inds[i] >= k[i]-1) {
00145             if (cont[i]) inds[i] = k[i]-1;
00146             else return -1;     // Outside potential and grid is not continuous
00147         }
00148         if (cont[i] && inds[i] == k[i]-1) {
00149             // Correct for non-unit spacing between continuous grid images
00150             gapscale[i] *= gapinv[i];
00151             if (g[i] < 0.0) dg[i] = 1.0 + g[i]*gapinv[i]; // = (gap[i] + g[i]) * gapinv[i]
00152             else dg[i] = (g[i] - inds[i]) * gapinv[i];
00153         }
00154     }
00155     
00156     return 0;
00157 }
00158 
00159 
00160 inline float GridforceFullBaseGrid::compute_V(float *a, float *x, float *y, float *z) const
00161 {
00162     float V = 0.0;
00163     long int ind = 0;
00164     for (int l = 0; l < 4; l++) {
00165         for (int k = 0; k < 4; k++) {
00166             for (int j = 0; j < 4; j++) {
00167                 V += a[ind] * x[j] * y[k] * z[l];
00168                 ind++;
00169             }
00170         }
00171     }
00172     return V;
00173 }
00174 
00175 
00176 inline Vector GridforceFullBaseGrid::compute_dV(float *a, float *x, float *y, float *z) const
00177 {
00178     Vector dV = 0;
00179     long int ind = 0;
00180     for (int l = 0; l < 4; l++) {
00181         for (int k = 0; k < 4; k++) {
00182             for (int j = 0; j < 4; j++) {
00183                 if (j > 0) dV.x += a[ind] * j * x[j-1] * y[k]   * z[l];         // dV/dx
00184                 if (k > 0) dV.y += a[ind] * k * x[j]   * y[k-1] * z[l];         // dV/dy
00185                 if (l > 0) dV.z += a[ind] * l * x[j]   * y[k]   * z[l-1];       // dV/dz
00186                 ind++;
00187             }
00188         }
00189     }
00190     return dV;
00191 }
00192 
00193 
00194 inline Vector GridforceFullBaseGrid::compute_d2V(float *a, float *x, float *y, float *z) const
00195 {
00196     Vector d2V = 0;
00197     int ind = 0;
00198     for (int l = 0; l < 4; l++) {
00199         for (int k = 0; k < 4; k++) {
00200             for (int j = 0; j < 4; j++) {
00201                 if (j > 0 && k > 0) d2V.x += a[ind] * j * k * x[j-1] * y[k-1] * z[l];   // d2V/dxdy
00202                 if (j > 0 && l > 0) d2V.y += a[ind] * j * l * x[j-1] * y[k]   * z[l-1]; // d2V/dxdz
00203                 if (k > 0 && l > 0) d2V.z += a[ind] * k * l * x[j]   * y[k-1] * z[l-1]; // d2V/dydz
00204                 ind++;
00205             }
00206         }
00207     }
00208     return d2V;
00209 }
00210 
00211 
00212 inline float GridforceFullBaseGrid::compute_d3V(float *a, float *x, float *y, float *z) const
00213 {
00214     float d3V = 0.0;
00215     long int ind = 0;
00216     for (int l = 0; l < 4; l++) {
00217         for (int k = 0; k < 4; k++) {
00218             for (int j = 0; j < 4; j++) {
00219                 if (j > 0 && k > 0 && l > 0) d3V += a[ind] * j * k * l * x[j-1] * y[k-1] * z[l-1];      // d3V/dxdydz
00220                 ind++;
00221             }
00222         }
00223     }
00224     return d3V;
00225 }
00226 
00227 
00228 inline void GridforceFullBaseGrid::compute_a(float *a, float *b) const
00229 {
00230     // Static sparse 64x64 matrix times vector ... nicer looking way than this?
00231     a[0] = b[0];
00232     a[1] = b[8];
00233     a[2] = -3*b[0] + 3*b[1] - 2*b[8] - b[9];
00234     a[3] = 2*b[0] - 2*b[1] + b[8] + b[9];
00235     a[4] = b[16];
00236     a[5] = b[32];
00237     a[6] = -3*b[16] + 3*b[17] - 2*b[32] - b[33];
00238     a[7] = 2*b[16] - 2*b[17] + b[32] + b[33];
00239     a[8] = -3*b[0] + 3*b[2] - 2*b[16] - b[18];
00240     a[9] = -3*b[8] + 3*b[10] - 2*b[32] - b[34];
00241     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]
00242         + 6*b[16] - 6*b[17] + 3*b[18] - 3*b[19] + 4*b[32] + 2*b[33] + 2*b[34] + b[35];
00243     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]
00244         - 4*b[16] + 4*b[17] - 2*b[18] + 2*b[19] - 2*b[32] - 2*b[33] - b[34] - b[35];
00245     a[12] = 2*b[0] - 2*b[2] + b[16] + b[18];
00246     a[13] = 2*b[8] - 2*b[10] + b[32] + b[34];
00247     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]
00248         - 3*b[16] + 3*b[17] - 3*b[18] + 3*b[19] - 2*b[32] - b[33] - 2*b[34] - b[35];
00249     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]
00250         + 2*b[16] - 2*b[17] + 2*b[18] - 2*b[19] + b[32] + b[33] + b[34] + b[35];
00251     a[16] = b[24];
00252     a[17] = b[40];
00253     a[18] = -3*b[24] + 3*b[25] - 2*b[40] - b[41];
00254     a[19] = 2*b[24] - 2*b[25] + b[40] + b[41];
00255     a[20] = b[48];
00256     a[21] = b[56];
00257     a[22] = -3*b[48] + 3*b[49] - 2*b[56] - b[57];
00258     a[23] = 2*b[48] - 2*b[49] + b[56] + b[57];
00259     a[24] = -3*b[24] + 3*b[26] - 2*b[48] - b[50];
00260     a[25] = -3*b[40] + 3*b[42] - 2*b[56] - b[58];
00261     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]
00262         + 6*b[48] - 6*b[49] + 3*b[50] - 3*b[51] + 4*b[56] + 2*b[57] + 2*b[58] + b[59];
00263     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]
00264         - 4*b[48] + 4*b[49] - 2*b[50] + 2*b[51] - 2*b[56] - 2*b[57] - b[58] - b[59];
00265     a[28] = 2*b[24] - 2*b[26] + b[48] + b[50];
00266     a[29] = 2*b[40] - 2*b[42] + b[56] + b[58];
00267     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]
00268         - 3*b[48] + 3*b[49] - 3*b[50] + 3*b[51] - 2*b[56] - b[57] - 2*b[58] - b[59];
00269     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]
00270         + 2*b[48] - 2*b[49] + 2*b[50] - 2*b[51] + b[56] + b[57] + b[58] + b[59];
00271     a[32] = -3*b[0] + 3*b[4] - 2*b[24] - b[28];
00272     a[33] = -3*b[8] + 3*b[12] - 2*b[40] - b[44];
00273     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]
00274         + 6*b[24] - 6*b[25] + 3*b[28] - 3*b[29] + 4*b[40] + 2*b[41] + 2*b[44] + b[45];
00275     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]
00276         - 4*b[24] + 4*b[25] - 2*b[28] + 2*b[29] - 2*b[40] - 2*b[41] - b[44] - b[45];
00277     a[36] = -3*b[16] + 3*b[20] - 2*b[48] - b[52];
00278     a[37] = -3*b[32] + 3*b[36] - 2*b[56] - b[60];
00279     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]
00280         + 6*b[48] - 6*b[49] + 3*b[52] - 3*b[53] + 4*b[56] + 2*b[57] + 2*b[60] + b[61];
00281     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]
00282         - 4*b[48] + 4*b[49] - 2*b[52] + 2*b[53] - 2*b[56] - 2*b[57] - b[60] - b[61];
00283     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]
00284         + 6*b[24] - 6*b[26] + 3*b[28] - 3*b[30] + 4*b[48] + 2*b[50] + 2*b[52] + b[54];
00285     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]
00286         + 6*b[40] - 6*b[42] + 3*b[44] - 3*b[46] + 4*b[56] + 2*b[58] + 2*b[60] + b[62];
00287     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]
00288         - 18*b[8] - 9*b[9] + 18*b[10] + 9*b[11] + 18*b[12] + 9*b[13] - 18*b[14] - 9*b[15]
00289         - 18*b[16] + 18*b[17] - 9*b[18] + 9*b[19] + 18*b[20] - 18*b[21] + 9*b[22] - 9*b[23]
00290         - 18*b[24] + 18*b[25] + 18*b[26] - 18*b[27] - 9*b[28] + 9*b[29] + 9*b[30] - 9*b[31]
00291         - 12*b[32] - 6*b[33] - 6*b[34] - 3*b[35] + 12*b[36] + 6*b[37] + 6*b[38] + 3*b[39]
00292         - 12*b[40] - 6*b[41] + 12*b[42] + 6*b[43] - 6*b[44] - 3*b[45] + 6*b[46] + 3*b[47]
00293         - 12*b[48] + 12*b[49] - 6*b[50] + 6*b[51] - 6*b[52] + 6*b[53] - 3*b[54] + 3*b[55]
00294         - 8*b[56] - 4*b[57] - 4*b[58] - 2*b[59] - 4*b[60] - 2*b[61] - 2*b[62] - b[63];
00295     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]
00296         + 9*b[8] + 9*b[9] - 9*b[10] - 9*b[11] - 9*b[12] - 9*b[13] + 9*b[14] + 9*b[15]
00297         + 12*b[16] - 12*b[17] + 6*b[18] - 6*b[19] - 12*b[20] + 12*b[21] - 6*b[22] + 6*b[23]
00298         + 12*b[24] - 12*b[25] - 12*b[26] + 12*b[27] + 6*b[28] - 6*b[29] - 6*b[30] + 6*b[31]
00299         + 6*b[32] + 6*b[33] + 3*b[34] + 3*b[35] - 6*b[36] - 6*b[37] - 3*b[38] - 3*b[39]
00300         + 6*b[40] + 6*b[41] - 6*b[42] - 6*b[43] + 3*b[44] + 3*b[45] - 3*b[46] - 3*b[47]
00301         + 8*b[48] - 8*b[49] + 4*b[50] - 4*b[51] + 4*b[52] - 4*b[53] + 2*b[54] - 2*b[55]
00302         + 4*b[56] + 4*b[57] + 2*b[58] + 2*b[59] + 2*b[60] + 2*b[61] + b[62] + b[63];
00303     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]
00304         - 4*b[24] + 4*b[26] - 2*b[28] + 2*b[30] - 2*b[48] - 2*b[50] - b[52] - b[54];
00305     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]
00306         - 4*b[40] + 4*b[42] - 2*b[44] + 2*b[46] - 2*b[56] - 2*b[58] - b[60] - b[62];
00307     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]
00308         + 12*b[8] + 6*b[9] - 12*b[10] - 6*b[11] - 12*b[12] - 6*b[13] + 12*b[14] + 6*b[15]
00309         + 9*b[16] - 9*b[17] + 9*b[18] - 9*b[19] - 9*b[20] + 9*b[21] - 9*b[22] + 9*b[23]
00310         + 12*b[24] - 12*b[25] - 12*b[26] + 12*b[27] + 6*b[28] - 6*b[29] - 6*b[30] + 6*b[31]
00311         + 6*b[32] + 3*b[33] + 6*b[34] + 3*b[35] - 6*b[36] - 3*b[37] - 6*b[38] - 3*b[39]
00312         + 8*b[40] + 4*b[41] - 8*b[42] - 4*b[43] + 4*b[44] + 2*b[45] - 4*b[46] - 2*b[47]
00313         + 6*b[48] - 6*b[49] + 6*b[50] - 6*b[51] + 3*b[52] - 3*b[53] + 3*b[54] - 3*b[55]
00314         + 4*b[56] + 2*b[57] + 4*b[58] + 2*b[59] + 2*b[60] + b[61] + 2*b[62] + b[63];
00315     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]
00316         - 6*b[8] - 6*b[9] + 6*b[10] + 6*b[11] + 6*b[12] + 6*b[13] - 6*b[14] - 6*b[15]
00317         - 6*b[16] + 6*b[17] - 6*b[18] + 6*b[19] + 6*b[20] - 6*b[21] + 6*b[22] - 6*b[23]
00318         - 8*b[24] + 8*b[25] + 8*b[26] - 8*b[27] - 4*b[28] + 4*b[29] + 4*b[30] - 4*b[31]
00319         - 3*b[32] - 3*b[33] - 3*b[34] - 3*b[35] + 3*b[36] + 3*b[37] + 3*b[38] + 3*b[39]
00320         - 4*b[40] - 4*b[41] + 4*b[42] + 4*b[43] - 2*b[44] - 2*b[45] + 2*b[46] + 2*b[47]
00321         - 4*b[48] + 4*b[49] - 4*b[50] + 4*b[51] - 2*b[52] + 2*b[53] - 2*b[54] + 2*b[55]
00322         - 2*b[56] - 2*b[57] - 2*b[58] - 2*b[59] - b[60] - b[61] - b[62] - b[63];
00323     a[48] = 2*b[0] - 2*b[4] + b[24] + b[28];
00324     a[49] = 2*b[8] - 2*b[12] + b[40] + b[44];
00325     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]
00326         - 3*b[24] + 3*b[25] - 3*b[28] + 3*b[29] - 2*b[40] - b[41] - 2*b[44] - b[45];
00327     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]
00328         + 2*b[24] - 2*b[25] + 2*b[28] - 2*b[29] + b[40] + b[41] + b[44] + b[45];
00329     a[52] = 2*b[16] - 2*b[20] + b[48] + b[52];
00330     a[53] = 2*b[32] - 2*b[36] + b[56] + b[60];
00331     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]
00332         - 3*b[48] + 3*b[49] - 3*b[52] + 3*b[53] - 2*b[56] - b[57] - 2*b[60] - b[61];
00333     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]
00334         + 2*b[48] - 2*b[49] + 2*b[52] - 2*b[53] + b[56] + b[57] + b[60] + b[61];
00335     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]
00336         - 3*b[24] + 3*b[26] - 3*b[28] + 3*b[30] - 2*b[48] - b[50] - 2*b[52] - b[54];
00337     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]
00338         - 3*b[40] + 3*b[42] - 3*b[44] + 3*b[46] - 2*b[56] - b[58] - 2*b[60] - b[62];
00339     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]
00340         + 12*b[8] + 6*b[9] - 12*b[10] - 6*b[11] - 12*b[12] - 6*b[13] + 12*b[14] + 6*b[15]
00341         + 12*b[16] - 12*b[17] + 6*b[18] - 6*b[19] - 12*b[20] + 12*b[21] - 6*b[22] + 6*b[23]
00342         + 9*b[24] - 9*b[25] - 9*b[26] + 9*b[27] + 9*b[28] - 9*b[29] - 9*b[30] + 9*b[31]
00343         + 8*b[32] + 4*b[33] + 4*b[34] + 2*b[35] - 8*b[36] - 4*b[37] - 4*b[38] - 2*b[39]
00344         + 6*b[40] + 3*b[41] - 6*b[42] - 3*b[43] + 6*b[44] + 3*b[45] - 6*b[46] - 3*b[47]
00345         + 6*b[48] - 6*b[49] + 3*b[50] - 3*b[51] + 6*b[52] - 6*b[53] + 3*b[54] - 3*b[55]
00346         + 4*b[56] + 2*b[57] + 2*b[58] + b[59] + 4*b[60] + 2*b[61] + 2*b[62] + b[63];
00347     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]
00348         - 6*b[8] - 6*b[9] + 6*b[10] + 6*b[11] + 6*b[12] + 6*b[13] - 6*b[14] - 6*b[15]
00349         - 8*b[16] + 8*b[17] - 4*b[18] + 4*b[19] + 8*b[20] - 8*b[21] + 4*b[22] - 4*b[23]
00350         - 6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 6*b[28] + 6*b[29] + 6*b[30] - 6*b[31]
00351         - 4*b[32] - 4*b[33] - 2*b[34] - 2*b[35] + 4*b[36] + 4*b[37] + 2*b[38] + 2*b[39]
00352         - 3*b[40] - 3*b[41] + 3*b[42] + 3*b[43] - 3*b[44] - 3*b[45] + 3*b[46] + 3*b[47]
00353         - 4*b[48] + 4*b[49] - 2*b[50] + 2*b[51] - 4*b[52] + 4*b[53] - 2*b[54] + 2*b[55]
00354         - 2*b[56] - 2*b[57] - b[58] - b[59] - 2*b[60] - 2*b[61] - b[62] - b[63];
00355     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]
00356         + 2*b[24] - 2*b[26] + 2*b[28] - 2*b[30] + b[48] + b[50] + b[52] + b[54];
00357     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]
00358         + 2*b[40] - 2*b[42] + 2*b[44] - 2*b[46] + b[56] + b[58] + b[60] + b[62];
00359     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]
00360         - 8*b[8] - 4*b[9] + 8*b[10] + 4*b[11] + 8*b[12] + 4*b[13] - 8*b[14] - 4*b[15]
00361         - 6*b[16] + 6*b[17] - 6*b[18] + 6*b[19] + 6*b[20] - 6*b[21] + 6*b[22] - 6*b[23]
00362         - 6*b[24] + 6*b[25] + 6*b[26] - 6*b[27] - 6*b[28] + 6*b[29] + 6*b[30] - 6*b[31]
00363         - 4*b[32] - 2*b[33] - 4*b[34] - 2*b[35] + 4*b[36] + 2*b[37] + 4*b[38] + 2*b[39]
00364         - 4*b[40] - 2*b[41] + 4*b[42] + 2*b[43] - 4*b[44] - 2*b[45] + 4*b[46] + 2*b[47]
00365         - 3*b[48] + 3*b[49] - 3*b[50] + 3*b[51] - 3*b[52] + 3*b[53] - 3*b[54] + 3*b[55]
00366         - 2*b[56] - b[57] - 2*b[58] - b[59] - 2*b[60] - b[61] - 2*b[62] - b[63];
00367     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]
00368         + 4*b[8] + 4*b[9] - 4*b[10] - 4*b[11] - 4*b[12] - 4*b[13] + 4*b[14] + 4*b[15]
00369         + 4*b[16] - 4*b[17] + 4*b[18] - 4*b[19] - 4*b[20] + 4*b[21] - 4*b[22] + 4*b[23]
00370         + 4*b[24] - 4*b[25] - 4*b[26] + 4*b[27] + 4*b[28] - 4*b[29] - 4*b[30] + 4*b[31]
00371         + 2*b[32] + 2*b[33] + 2*b[34] + 2*b[35] - 2*b[36] - 2*b[37] - 2*b[38] - 2*b[39]
00372         + 2*b[40] + 2*b[41] - 2*b[42] - 2*b[43] + 2*b[44] + 2*b[45] - 2*b[46] - 2*b[47]
00373         + 2*b[48] - 2*b[49] + 2*b[50] - 2*b[51] + 2*b[52] - 2*b[53] + 2*b[54] - 2*b[55]
00374         + b[56] + b[57] + b[58] + b[59] + b[60] + b[61] + b[62] + b[63];
00375 }
00376 
00377 
00378 inline int GridforceLiteGrid::get_inds(Position pos, int *inds, Vector &dg) const
00379 {
00380     Vector p = pos - origin;
00381     Vector g;
00382     
00383     g = inv * p;
00384     
00385     for (int i = 0; i < 3; i++) {
00386         inds[i] = (int)floor(g[i]);
00387         dg[i] = g[i] - inds[i];
00388     }
00389     
00390     for (int i = 0; i < 3; i++) {
00391         if (inds[i] < 0 || inds[i] >= k[i]-1) {
00392             return -1;  // Outside potential and grid is not continuous
00393         }
00394     }
00395     
00396     return 0;
00397 }
00398 
00399 
00400 inline void GridforceLiteGrid::compute_wts(float *wts, const Vector &dg) const
00401 {
00402     wts[0] = (1-dg.x) * (1-dg.y) * (1-dg.z);
00403     wts[1] = (1-dg.x) * (1-dg.y) *   dg.z;
00404     wts[2] = (1-dg.x) *   dg.y   * (1-dg.z);
00405     wts[3] = (1-dg.x) *   dg.y   *   dg.z;
00406     wts[4] =   dg.x   * (1-dg.y) * (1-dg.z);
00407     wts[5] =   dg.x   * (1-dg.y) *   dg.z;
00408     wts[6] =   dg.x   *   dg.y   * (1-dg.z);
00409     wts[7] =   dg.x   *   dg.y   *   dg.z;
00410     DebugM(2, "dg = " << dg << "\n" << endi);
00411 }
00412 
00413 
00414 inline float GridforceLiteGrid::linear_interpolate(int i0, int i1, int i2, int i3, const float *wts) const
00415 {
00416 #ifdef DEBUGM
00417     float vals[8];
00418     vals[0] = get_grid(i0,   i1,   i2,   i3);
00419     vals[1] = get_grid(i0,   i1,   i2+1, i3);
00420     vals[2] = get_grid(i0,   i1+1, i2,   i3);
00421     vals[3] = get_grid(i0,   i1+1, i2+1, i3);
00422     vals[4] = get_grid(i0+1, i1,   i2,   i3);
00423     vals[5] = get_grid(i0+1, i1,   i2+1, i3);
00424     vals[6] = get_grid(i0+1, i1+1, i2,   i3);
00425     vals[7] = get_grid(i0+1, i1+1, i2+1, i3);
00426     
00427     switch (i3) {
00428     case 0:
00429         DebugM(2, "V\n" << endi);
00430         break;
00431     case 1:
00432         DebugM(2, "dV/dx\n" << endi);
00433         break;
00434     case 2:
00435         DebugM(2, "dV/dy\n" << endi);
00436         break;
00437     case 3:
00438         DebugM(2, "dV/dz\n" << endi);
00439         break;
00440     }
00441     
00442     for (int i = 0; i < 8; i++) {
00443         DebugM(2, "vals[" << i << "] = " << vals[i] << " wts[" << i << "] = " << wts[i] << "\n" << endi);
00444     }
00445 #endif
00446     
00447     float result =
00448         wts[0] * get_grid(i0,   i1,   i2,   i3) +
00449         wts[1] * get_grid(i0,   i1,   i2+1, i3) +
00450         wts[2] * get_grid(i0,   i1+1, i2,   i3) +
00451         wts[3] * get_grid(i0,   i1+1, i2+1, i3) +
00452         wts[4] * get_grid(i0+1, i1,   i2,   i3) +
00453         wts[5] * get_grid(i0+1, i1,   i2+1, i3) +
00454         wts[6] * get_grid(i0+1, i1+1, i2,   i3) +
00455         wts[7] * get_grid(i0+1, i1+1, i2+1, i3);
00456     
00457     DebugM(2, "result = " << result << "\n" << endi);
00458     
00459     return result;
00460 }
00461 
00462 
00463 inline Position GridforceGrid::wrap_position(const Position &pos, const Lattice &lattice)
00464 {
00465     // Wrap 'pos' about grid center, using periodic cell information in 'lattice'
00466     // Position pos_wrapped = pos;
00467     // Position center = get_center();
00468     // pos_wrapped += lattice.wrap_delta(pos);
00469     // pos_wrapped += lattice.delta(pos_wrapped, center) - (pos_wrapped - center);
00470     
00471     Position pos_wrapped = pos + lattice.wrap_delta(pos - get_center() + lattice.origin());
00472     
00473     return pos_wrapped;
00474 }
00475 
00476 #endif

Generated on Tue Sep 19 01:17:12 2017 for NAMD by  doxygen 1.4.7