diff --git a/BaseGrid.cu b/BaseGrid.cu
index 6db492baa1966661e7286c50dfb1c3260e905056..b5cb8d6a2e7f63bdfbcc48fb33f85a2ba538968a 100644
--- a/BaseGrid.cu
+++ b/BaseGrid.cu
@@ -622,83 +622,6 @@ bool BaseGrid::crop(int x0, int y0, int z0, int x1, int y1, int z1, bool keep_or
return true;
}
-// Added by Rogan for times when simpler calculations are required.
-float BaseGrid::interpolatePotentialLinearly(Vector3 pos) const {
- // Find the home node.
- Vector3 l = basisInv.transform(pos - origin);
- int homeX = int(floorf(l.x));
- int homeY = int(floorf(l.y));
- int homeZ = int(floorf(l.z));
-
- // Get the array jumps.
- int jump[3];
- jump[0] = nz*ny;
- jump[1] = nz;
- jump[2] = 1;
-
- // Shift the indices in the home array.
- int home[3];
- home[0] = homeX;
- home[1] = homeY;
- home[2] = homeZ;
-
- // Get the grid dimensions.
- int g[3];
- g[0] = nx;
- g[1] = ny;
- g[2] = nz;
-
- // Get the interpolation coordinates.
- float w[3];
- w[0] = l.x - homeX;
- w[1] = l.y - homeY;
- w[2] = l.z - homeZ;
-
- // Find the values at the neighbors.
- float g1[2][2][2];
- for (int ix = 0; ix < 2; ix++) {
- for (int iy = 0; iy < 2; iy++) {
- for (int iz = 0; iz < 2; iz++) {
- // Wrap around the periodic boundaries.
- int jx = ix + home[0];
- jx = wrap(jx, g[0]);
- int jy = iy + home[1];
- jy = wrap(jy, g[1]);
- int jz = iz + home[2];
- jz = wrap(jz, g[2]);
-
- int ind = jz*jump[2] + jy*jump[1] + jx*jump[0];
- g1[ix][iy][iz] = val[ind];
- }
- }
- }
-
- // Mix along x.
- float g2[2][2];
- for (int iy = 0; iy < 2; iy++) {
- for (int iz = 0; iz < 2; iz++) {
- // p = w[0] * g[0][iy][iz] + (1-w[0]) * g[1][iy][iz]
- g2[iy][iz] = (1.0f-w[0])*g1[0][iy][iz] + w[0]*g1[1][iy][iz];
- }
- }
-
- // Mix along y.
- float g3[2];
- for (int iz = 0; iz < 2; iz++) {
- g3[iz] = (1.0f-w[1])*g2[0][iz] + w[1]*g2[1][iz];
- }
-
- // DEBUG
- //printf("(0,0,0)=%.1f (0,0,1)=%.1f (0,1,0)=%.1f (0,1,1)=%.1f (1,0,0)=%.1f (1,0,1)=%.1f (1,1,0)=%.1f (1,1,1)=%.1f ",
- // g1[0][0][0], g1[0][0][1], g1[0][1][0], g1[0][1][1], g1[1][0][0], g1[1][0][1], g1[1][1][0], g1[1][1][1] );
- //printf ("%.2f\n",(1.0-w[2])*g3[0] + w[2]*g3[1]);
-
- // Mix along z
- return (1.0f-w[2])*g3[0] + w[2]*g3[1];
-}
-
-
-
Vector3 BaseGrid::wrapDiffNearest(Vector3 r) const {
Vector3 l = basisInv.transform(r);
l.x = wrapDiff(l.x, nx);
diff --git a/BaseGrid.h b/BaseGrid.h
index fe4c127cf1d838dad4e1cfc4aec23e630142260f..ce171dbecb64bc641bf7c3061d57aa7098ede86f 100644
--- a/BaseGrid.h
+++ b/BaseGrid.h
@@ -31,6 +31,14 @@ public:
float v[3][3][3];
};
+class ForceEnergy {
+public:
+ DEVICE ForceEnergy(Vector3 &f, float &e) :
+ f(f), e(e) {};
+ Vector3 f;
+ float e;
+};
+
class BaseGrid {
friend class SparseGrid;
@@ -192,7 +200,7 @@ public:
bool crop(int x0, int y0, int z0, int x1, int y1, int z1, bool keep_origin);
// Added by Rogan for times when simpler calculations are required.
- virtual float interpolatePotentialLinearly(Vector3 pos) const;
+ // virtual float interpolatePotentialLinearly(Vector3 pos) const;
HOST DEVICE inline float interpolateDiffX(Vector3 pos, float w[3], float g1[4][4][4]) const {
float a0, a1, a2, a3;
@@ -308,91 +316,51 @@ public:
return -(3.0f*a3*w[2]*w[2] + 2.0f*a2*w[2] + a1);
}
- HOST DEVICE inline float interpolatePotential(Vector3 pos) const {
+ HOST DEVICE inline float interpolatePotential(const Vector3& pos) const {
// Find the home node.
Vector3 l = basisInv.transform(pos - origin);
- int homeX = int(floor(l.x));
- int homeY = int(floor(l.y));
- int homeZ = int(floor(l.z));
-
- // out of grid? return 0
- // RBTODO
-
- // Get the array jumps.
- int jump[3];
- jump[0] = nz*ny;
- jump[1] = nz;
- jump[2] = 1;
-
- // Shift the indices in the home array.
- int home[3];
- home[0] = homeX;
- home[1] = homeY;
- home[2] = homeZ;
-
- // Get the grid dimensions.
- int g[3];
- g[0] = nx;
- g[1] = ny;
- g[2] = nz;
-
- // Get the interpolation coordinates.
- float w[3];
- w[0] = l.x - homeX;
- w[1] = l.y - homeY;
- w[2] = l.z - homeZ;
-
- // Find the values at the neighbors.
- float g1[4][4][4];
- for (int ix = 0; ix < 4; ix++) {
- int jx = ix-1 + home[0];
- jx = wrap(jx, g[0]);
- for (int iy = 0; iy < 4; iy++) {
- int jy = iy-1 + home[1];
- jy = wrap(jy, g[1]);
- for (int iz = 0; iz < 4; iz++) {
- // Wrap around the periodic boundaries.
- int jz = iz-1 + home[2];
- jz = wrap(jz, g[2]);
-
- int ind = jz*jump[2] + jy*jump[1] + jx*jump[0];
- g1[ix][iy][iz] = val[ind];
- }
- }
- }
- float a0, a1, a2, a3;
-
- // Mix along x.
- float g2[4][4];
- for (int iy = 0; iy < 4; iy++) {
- for (int iz = 0; iz < 4; iz++) {
- a3 = 0.5f*(-g1[0][iy][iz] + 3.0f*g1[1][iy][iz] - 3.0f*g1[2][iy][iz] + g1[3][iy][iz]);
- a2 = 0.5f*(2.0f*g1[0][iy][iz] - 5.0f*g1[1][iy][iz] + 4.0f*g1[2][iy][iz] - g1[3][iy][iz]);
- a1 = 0.5f*(-g1[0][iy][iz] + g1[2][iy][iz]);
- a0 = g1[1][iy][iz];
-
- g2[iy][iz] = a3*w[0]*w[0]*w[0] + a2*w[0]*w[0] + a1*w[0] + a0;
- }
- }
- // Mix along y.
- float g3[4];
- for (int iz = 0; iz < 4; iz++) {
- a3 = 0.5f*(-g2[0][iz] + 3.0f*g2[1][iz] - 3.0f*g2[2][iz] + g2[3][iz]);
- a2 = 0.5f*(2.0f*g2[0][iz] - 5.0f*g2[1][iz] + 4.0f*g2[2][iz] - g2[3][iz]);
- a1 = 0.5f*(-g2[0][iz] + g2[2][iz]);
- a0 = g2[1][iz];
-
- g3[iz] = a3*w[1]*w[1]*w[1] + a2*w[1]*w[1] + a1*w[1] + a0;
- }
-
- // Mix along z.
- a3 = 0.5f*(-g3[0] + 3.0f*g3[1] - 3.0f*g3[2] + g3[3]);
- a2 = 0.5f*(2.0f*g3[0] - 5.0f*g3[1] + 4.0f*g3[2] - g3[3]);
- a1 = 0.5f*(-g3[0] + g3[2]);
- a0 = g3[1];
+ const int homeX = int(floor(l.x));
+ const int homeY = int(floor(l.y));
+ const int homeZ = int(floor(l.z));
+ const float wx = l.x - homeX;
+ const float wy = l.y - homeY;
+ const float wz = l.z - homeZ;
+ const float wx2 = wx*wx;
+ const float wy2 = wy*wy;
+ const float wz2 = wz*wz;
+
+ float g3[4];
+ for (int iz = 0; iz < 4; iz++) {
+ float g2[4];
+ const int jz = (iz + homeZ - 1);
+ for (int iy = 0; iy < 4; iy++) {
+ float v[4];
+ const int jy = (iy + homeY - 1);
+ for (int ix = 0; ix < 4; ix++) {
+ const int jx = (ix + homeX - 1);
+ const int ind = jz + jy*nz + jx*nz*ny;
+ v[ix] = jz < 0 || jz >= nz || jy < 0 || jy >= ny || jx < 0 || jx >= nx ?
+ 0 : val[ind];
+ }
+ g2[iy] = 0.5f*(-v[0] + 3.0f*v[1] - 3.0f*v[2] + v[3])*wx2*wx +
+ 0.5f*(2.0f*v[0] - 5.0f*v[1] + 4.0f*v[2] - v[3]) *wx2 +
+ 0.5f*(-v[0] + v[2]) *wx +
+ v[1];
+ }
- return a3*w[2]*w[2]*w[2] + a2*w[2]*w[2] + a1*w[2] + a0;
+ // Mix along y.
+ g3[iz] = 0.5f*(-g2[0] + 3.0f*g2[1] - 3.0f*g2[2] + g2[3])*wy2*wy +
+ 0.5f*(2.0f*g2[0] - 5.0f*g2[1] + 4.0f*g2[2] - g2[3]) *wy2 +
+ 0.5f*(-g2[0] + g2[2]) *wy +
+ g2[1];
+ }
+ // Mix along z.
+ const float e = 0.5f*(-g3[0] + 3.0f*g3[1] - 3.0f*g3[2] + g3[3])*wz2*wz +
+ 0.5f*(2.0f*g3[0] - 5.0f*g3[1] + 4.0f*g3[2] - g3[3]) *wz2 +
+ 0.5f*(-g3[0] + g3[2]) *wz +
+ g3[1];
+ return e;
}
HOST DEVICE inline static int wrap(int i, int n) {
@@ -406,58 +374,168 @@ public:
}
return i;
}
+ HOST DEVICE float interpolatePotentialLinearly(const Vector3& pos) const {
+ Vector3 f;
+ const Vector3 l = basisInv.transform(pos - origin);
+
+ // Find the home node.
+ const int homeX = int(floor(l.x));
+ const int homeY = int(floor(l.y));
+ const int homeZ = int(floor(l.z));
+
+ const float wx = l.x - homeY;
+ const float wy = l.y - homeY;
+ const float wz = l.z - homeZ;
+
+ float v[2][2][2];
+ for (int iz = 0; iz < 2; iz++) {
+ int jz = (iz + homeZ);
+ for (int iy = 0; iy < 2; iy++) {
+ int jy = (iy + homeY);
+ for (int ix = 0; ix < 2; ix++) {
+ int jx = (ix + homeX);
+ int ind = jz + jy*nz + jx*nz*ny;
+ v[ix][iy][iz] = jz < 0 || jz >= nz || jy < 0 || jy >= ny || jx < 0 || jx >= nx ?
+ 0 : val[ind];
+ }
+ }
+ }
+
+ float g3[2];
+ for (int iz = 0; iz < 2; iz++) {
+ float g2[2];
+ for (int iy = 0; iy < 2; iy++) {
+ g2[iy] = wx * (v[1][iy][iz] - v[0][iy][iz]) + v[0][iy][iz];
+ }
+ // Mix along y.
+ g3[iz] = wy * (g2[1] - g2[0]) + g2[0];
+ }
+ // Mix along z.
+ float e = wz * (g3[1] - g3[0]) + g3[0];
+ return e;
+ }
+
/** interpolateForce() to be used on CUDA Device **/
- DEVICE inline Vector3 interpolateForceD(Vector3 pos) const {
+ DEVICE inline ForceEnergy interpolateForceD(const Vector3& pos) const {
Vector3 f;
- Vector3 l = basisInv.transform(pos - origin);
- int homeX = int(floor(l.x));
- int homeY = int(floor(l.y));
- int homeZ = int(floor(l.z));
- // Get the array jumps with shifted indices.
- int jump[3];
- jump[0] = nz*ny;
- jump[1] = nz;
- jump[2] = 1;
- // Shift the indices in the home array.
- int home[3];
- home[0] = homeX;
- home[1] = homeY;
- home[2] = homeZ;
+ const Vector3 l = basisInv.transform(pos - origin);
+
+ const int homeX = int(floor(l.x));
+ const int homeY = int(floor(l.y));
+ const int homeZ = int(floor(l.z));
+ const float wx = l.x - homeX;
+ const float wy = l.y - homeY;
+ const float wz = l.z - homeZ;
+ const float wx2 = wx*wx;
+
+ /* f.x */
+ float g3[3][4];
+ for (int iz = 0; iz < 4; iz++) {
+ float g2[2][4];
+ const int jz = (iz + homeZ - 1);
+ for (int iy = 0; iy < 4; iy++) {
+ float v[4];
+ const int jy = (iy + homeY - 1);
+ for (int ix = 0; ix < 4; ix++) {
+ const int jx = (ix + homeX - 1);
+ const int ind = jz + jy*nz + jx*nz*ny;
+ v[ix] = jz < 0 || jz >= nz || jy < 0 || jy >= ny || jx < 0 || jx >= nx ?
+ 0 : val[ind];
+ }
+ const float a3 = 0.5f*(-v[0] + 3.0f*v[1] - 3.0f*v[2] + v[3])*wx2;
+ const float a2 = 0.5f*(2.0f*v[0] - 5.0f*v[1] + 4.0f*v[2] - v[3])*wx;
+ const float a1 = 0.5f*(-v[0] + v[2]);
+ g2[0][iy] = 3.0f*a3 + 2.0f*a2 + a1; /* f.x (derivative) */
+ g2[1][iy] = a3*wx + a2*wx + a1*wx + v[1]; /* f.y & f.z */
+ }
- // Shift the indices in the grid dimensions.
- int g[3];
- g[0] = nx;
- g[1] = ny;
- g[2] = nz;
+ // Mix along y.
+ {
+ g3[0][iz] = 0.5f*(-g2[0][0] + 3.0f*g2[0][1] - 3.0f*g2[0][2] + g2[0][3])*wy*wy*wy +
+ 0.5f*(2.0f*g2[0][0] - 5.0f*g2[0][1] + 4.0f*g2[0][2] - g2[0][3]) *wy*wy +
+ 0.5f*(-g2[0][0] + g2[0][2]) *wy +
+ g2[0][1];
+ }
- // Get the interpolation coordinates.
- float w[3];
- w[0] = l.x - homeX;
- w[1] = l.y - homeY;
- w[2] = l.z - homeZ;
- // Find the values at the neighbors.
- float g1[4][4][4];
- for (int ix = 0; ix < 4; ix++) {
- for (int iy = 0; iy < 4; iy++) {
- for (int iz = 0; iz < 4; iz++) {
- // Wrap around the periodic boundaries.
- int jx = ix-1 + home[0];
- jx = wrap(jx, g[0]);
- int jy = iy-1 + home[1];
- jy = wrap(jy, g[1]);
- int jz = iz-1 + home[2];
- jz = wrap(jz, g[2]);
- int ind = jz*jump[2] + jy*jump[1] + jx*jump[0];
- g1[ix][iy][iz] = val[ind];
+ {
+ const float a3 = 0.5f*(-g2[1][0] + 3.0f*g2[1][1] - 3.0f*g2[1][2] + g2[1][3])*wy*wy;
+ const float a2 = 0.5f*(2.0f*g2[1][0] - 5.0f*g2[1][1] + 4.0f*g2[1][2] - g2[1][3])*wy;
+ const float a1 = 0.5f*(-g2[1][0] + g2[1][2]);
+ g3[1][iz] = 3.0f*a3 + 2.0f*a2 + a1; /* f.y */
+ g3[2][iz] = a3*wy + a2*wy + a1*wy + g2[1][1]; /* f.z */
+ }
+ }
+
+ // Mix along z.
+ f.x = -0.5f*(-g3[0][0] + 3.0f*g3[0][1] - 3.0f*g3[0][2] + g3[0][3])*wz*wz*wz +
+ -0.5f*(2.0f*g3[0][0] - 5.0f*g3[0][1] + 4.0f*g3[0][2] - g3[0][3])*wz*wz +
+ -0.5f*(-g3[0][0] + g3[0][2]) *wz -
+ g3[0][1];
+ f.y = -0.5f*(-g3[1][0] + 3.0f*g3[1][1] - 3.0f*g3[1][2] + g3[1][3])*wz*wz*wz +
+ -0.5f*(2.0f*g3[1][0] - 5.0f*g3[1][1] + 4.0f*g3[1][2] - g3[1][3])*wz*wz +
+ -0.5f*(-g3[1][0] + g3[1][2]) *wz -
+ g3[1][1];
+ f.z = -1.5f*(-g3[2][0] + 3.0f*g3[2][1] - 3.0f*g3[2][2] + g3[2][3])*wz*wz -
+ (2.0f*g3[2][0] - 5.0f*g3[2][1] + 4.0f*g3[2][2] - g3[2][3]) *wz -
+ 0.5f*(-g3[2][0] + g3[2][2]);
+ float e = 0.5f*(-g3[2][0] + 3.0f*g3[2][1] - 3.0f*g3[2][2] + g3[2][3])*wz*wz*wz +
+ 0.5f*(2.0f*g3[2][0] - 5.0f*g3[2][1] + 4.0f*g3[2][2] - g3[2][3]) *wz*wz +
+ 0.5f*(-g3[2][0] + g3[2][2]) *wz +
+ g3[2][1];
+
+ f = basisInv.transpose().transform(f);
+ return ForceEnergy(f,e);
+
+ }
+
+ DEVICE inline ForceEnergy interpolateForceDLinearly(const Vector3& pos) const {
+ Vector3 f;
+ const Vector3 l = basisInv.transform(pos - origin);
+
+ // Find the home node.
+ const int homeX = int(floor(l.x));
+ const int homeY = int(floor(l.y));
+ const int homeZ = int(floor(l.z));
+
+ const float wx = l.x - homeY;
+ const float wy = l.y - homeY;
+ const float wz = l.z - homeZ;
+
+ float v[2][2][2];
+ for (int iz = 0; iz < 2; iz++) {
+ int jz = (iz + homeZ);
+ for (int iy = 0; iy < 2; iy++) {
+ int jy = (iy + homeY);
+ for (int ix = 0; ix < 2; ix++) {
+ int jx = (ix + homeX);
+ int ind = jz + jy*nz + jx*nz*ny;
+ v[ix][iy][iz] = jz < 0 || jz >= nz || jy < 0 || jy >= ny || jx < 0 || jx >= nx ?
+ 0 : val[ind];
}
}
- }
- f.x = interpolateDiffX(pos, w, g1);
- f.y = interpolateDiffY(pos, w, g1);
- f.z = interpolateDiffZ(pos, w, g1);
- Vector3 f1 = basisInv.transpose().transform(f);
- return f1;
+ }
+
+ float g3[3][2];
+ for (int iz = 0; iz < 2; iz++) {
+ float g2[2][2];
+ for (int iy = 0; iy < 2; iy++) {
+ g2[0][iy] = (v[1][iy][iz] - v[0][iy][iz]); /* f.x */
+ g2[1][iy] = wx * (v[1][iy][iz] - v[0][iy][iz]) + v[0][iy][iz]; /* f.y & f.z */
+ }
+ // Mix along y.
+ g3[0][iz] = wy * (g2[0][1] - g2[0][0]) + g2[0][0];
+ g3[1][iz] = (g2[1][1] - g2[1][0]);
+ g3[2][iz] = wy * (g2[1][1] - g2[1][0]) + g2[1][0];
+ }
+ // Mix along z.
+ f.x = -(wz * (g3[0][1] - g3[0][0]) + g3[0][0]);
+ f.y = -(wz * (g3[1][1] - g3[1][0]) + g3[1][0]);
+ f.z = - (g3[2][1] - g3[2][0]);
+
+ f = basisInv.transpose().transform(f);
+ float e = wz * (g3[2][1] - g3[2][0]) + g3[2][0];
+ return ForceEnergy(f,e);
}
inline virtual Vector3 interpolateForce(Vector3 pos) const {
diff --git a/GrandBrownTown.cuh b/GrandBrownTown.cuh
index 5f8c1a7e4dd057acb324b34a166ad15505b2ada6..fb701140aba0367fadea9650d4e4b6177f3e4dcf 100644
--- a/GrandBrownTown.cuh
+++ b/GrandBrownTown.cuh
@@ -33,65 +33,21 @@ void updateKernel(Vector3 pos[], Vector3 forceInternal[],
// Compute external force
Vector3 forceExternal = Vector3(0.0f, 0.0f, pt.charge * electricField);
- BaseGrid* pmf = pt.pmf;
-
- // Find the home node.
- Vector3 l = pmf->getInverseBasis().transform(p - pmf->getOrigin());
- int homeX = int(floorf(l.x));
- int homeY = int(floorf(l.y));
- int homeZ = int(floorf(l.z));
-
- // Shift the indices in the grid dimensions.
- int gX = pmf->getNx();
- int gY = pmf->getNy();
- int gZ = pmf->getNz();
-
- // Get the interpolation coordinates.
- float w[3];
- w[0] = l.x - homeX;
- w[1] = l.y - homeY;
- w[2] = l.z - homeZ;
-
- // Find the values at the neighbors.
- float g1[4][4][4];
- for (int ix = 0; ix < 4; ix++) {
- int jx = ix-1 + homeX;
- jx = pmf->wrap(jx, gX);
-
- for (int iy = 0; iy < 4; iy++) {
- int jy = iy-1 + homeY;
- jy = pmf->wrap(jy, gY);
-
- for (int iz = 0; iz < 4; iz++) {
- int jz = iz-1 + homeZ;
- jz = pmf->wrap(jz, gZ);
- // Wrap around the periodic boundaries.
- int ind = jz + jy*gZ + jx*gZ*gY;
- g1[ix][iy][iz] = pmf->val[ind];
- }
- }
- }
-
- // Interpolate this particle's change in X, Y, and Z axes
- Vector3 f;
- f.x = pmf->interpolateDiffX(p, w, g1);
- f.y = pmf->interpolateDiffY(p, w, g1);
- f.z = pmf->interpolateDiffZ(p, w, g1);
- Vector3 forceGrid = pmf->getInverseBasis().transpose().transform(f);
-
+ // Compute PMF
+ ForceEnergy fe = pt.pmf->interpolateForceD(p);
#ifndef FORCEGRIDOFF
// Add a force defined via 3D FORCE maps (not 3D potential maps)
- if (pt.forceXGrid != NULL) forceGrid.x += pt.forceXGrid->interpolatePotential(p);
- if (pt.forceYGrid != NULL) forceGrid.y += pt.forceYGrid->interpolatePotential(p);
- if (pt.forceZGrid != NULL) forceGrid.z += pt.forceZGrid->interpolatePotential(p);
+ if (pt.forceXGrid != NULL) fe.f.x += pt.forceXGrid->interpolatePotential(p);
+ if (pt.forceYGrid != NULL) fe.f.y += pt.forceYGrid->interpolatePotential(p);
+ if (pt.forceZGrid != NULL) fe.f.z += pt.forceZGrid->interpolatePotential(p);
#endif
// Compute total force:
// Internal: interaction between particles
// External: electric field (now this is basically a constant vector)
// forceGrid: ADD force due to PMF or other potentials defined in 3D space
- Vector3 force = forceInternal[idx] + forceExternal + forceGrid;
+ Vector3 force = forceInternal[idx] + forceExternal + fe.f;
if (idx == 0)
forceInternal[idx] = force; // write it back out for force0 in run()
@@ -104,7 +60,7 @@ void updateKernel(Vector3 pos[], Vector3 forceInternal[],
p = step(p, kTlocal, force, pt.diffusion, timestep, sys, randoGen, num);
} else {
float diffusion = pt.diffusionGrid->interpolatePotential(p);
- Vector3 diffGrad = pt.diffusionGrid->interpolateForceD(p);
+ Vector3 diffGrad = (pt.diffusionGrid->interpolateForceD(p)).f;
p = step(p, kTlocal, force, diffusion, -diffGrad, timestep, sys, randoGen, num);
}
}
diff --git a/OverlordGrid.h b/OverlordGrid.h
index 22cbf13f2502c399e2bef0465e08b54f397fb94d..060f909643a234a33382db97b392c3562ebe9e16 100644
--- a/OverlordGrid.h
+++ b/OverlordGrid.h
@@ -8,6 +8,14 @@
#include "BaseGrid.h"
#include "useful.h"
+#ifdef __CUDACC__
+ #define HOST __host__
+ #define DEVICE __device__
+#else
+ #define HOST
+ #define DEVICE
+#endif
+
class OverlordGrid : public BaseGrid {
public:
OverlordGrid(const BaseGrid& grid) : BaseGrid(grid) {
@@ -222,7 +230,7 @@ public:
return subgrid[j]->getPotential(r);
}
- virtual float interpolatePotential(Vector3 pos) const {
+ DEVICE virtual float interpolatePotential(Vector3 pos) const {
// Find the nearest node.
int j = nearestIndex(pos);
@@ -235,7 +243,7 @@ public:
return subgrid[j]->interpolatePotential(r);
}
- virtual float interpolatePotentialLinearly(Vector3 pos) const {
+ DEVICE virtual float interpolatePotentialLinearly(Vector3 pos) const {
// Find the nearest node.
int j = nearestIndex(pos);
diff --git a/RigidBodyGrid.h b/RigidBodyGrid.h
index 6206b9dc27e460acac7d688cdceee52e291ba005..bcc44d1294862511b17d5d81b068a8cb73fe9b26 100644
--- a/RigidBodyGrid.h
+++ b/RigidBodyGrid.h
@@ -26,13 +26,13 @@ using namespace std;
#define STRLEN 512
-class ForceEnergy {
-public:
- DEVICE ForceEnergy(Vector3 &f, float &e) :
- f(f), e(e) {};
- Vector3 f;
- float e;
-};
+/* class ForceEnergy { */
+/* public: */
+/* DEVICE ForceEnergy(Vector3 &f, float &e) : */
+/* f(f), e(e) {}; */
+/* Vector3 f; */
+/* float e; */
+/* }; */
class RigidBodyGrid {
friend class SparseGrid;