Skip to content
Snippets Groups Projects
Commit 7ac4cdf8 authored by cmaffeo2's avatar cmaffeo2
Browse files

Grid-Grid torque calculation now correctly handles wrapping

parent bd134e1d
No related branches found
No related tags found
No related merge requests found
......@@ -2,7 +2,6 @@
#include "ComputeGridGrid.cuh"
#include "RigidBodyGrid.h"
#include "CudaUtil.cuh"
//RBTODO handle periodic boundaries
//RBTODO: add __restrict__, benchmark (Q: how to restrict member data?)
class GridPositionTransformer {
......@@ -17,6 +16,7 @@ private:
const Vector3 c;
const BaseGrid* s;
};
//class PmfPositionTransformer : public BasePositionTransformer {
class PmfPositionTransformer {
public:
......@@ -52,8 +52,7 @@ inline void common_computeGridGridForce(const RigidBodyGrid* rho, const RigidBod
Vector3 r_pos= rho->getPosition(r_id); /* i,j,k value of voxel */
r_pos = basis_rho.transform( r_pos );
r_pos = transformer(r_pos);
const Vector3 u_ijk_float = basis_u_inv.transform( r_pos );
const Vector3 u_ijk_float = basis_u_inv.transform( transformer( r_pos ) );
// RBTODO: Test for non-unit delta
/* Vector3 tmpf = Vector3(0.0f); */
/* float tmpe = 0.0f; */
......@@ -71,7 +70,8 @@ inline void common_computeGridGridForce(const RigidBodyGrid* rho, const RigidBod
const float r_val = rho->val[r_id]; /* maybe move to beginning of function? */
force[tid].f = basis_u_inv.transpose().transform( r_val*(force[tid].f) ); /* transform to lab frame, with correct scaling factor */
force[tid].e = r_val;
// Calculate torque about origin_u in the lab frame
// Calculate torque about origin_rho in the lab frame
torque[tid].f = r_pos.cross(force[tid].f);
}
......
......@@ -902,10 +902,10 @@ void RigidBodyForcePair::processGPUForces(BaseGrid* sys) {
tmpT = tmpT + torques[i][j];
}
// tmpT is the torque calculated about the origin of grid k2 (e.g. c2)
// tmpT is the torque calculated about the origin of density grid
// so here we transform torque to be about rb1
Vector3 o2 = getOrigin2(i);
tmpT = tmpT - sys->wrapDiff(rb1->getPosition() - o2).cross( tmpF.f );
Vector3 o1 = getOrigin1(i);
tmpT = tmpT - (rb1->getPosition() - o1).cross( tmpF.f );
// clear forces on GPU
gpuErrchk(cudaMemset((void*)(forces_d[i]),0,nb*sizeof(ForceEnergy)));
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment