1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75

/*!
 * 通過單元內粒子見得距離計算密度
 * @param[in] gridPos 網格位置
 * @param[in] index 粒子索引
 * @param[in] pos 計算座標
 * @param[in] cell 粒子網格資料
 * @return 密度值
 */
__device__
float calDensityCell(int3 gridPos, uint i, float3 pos0, rxParticleCell cell)
{
    uint gridHash = calcGridHash(gridPos);
 
    // cell內粒子的開始索引
    uint startIndex = cell.dCellStart[gridHash];
 
    float h = params.EffectiveRadius;
    float dens = 0.0f;
    if(startIndex != 0xffffffff){    // cell非空
        // cell內粒子迴圈
        uint endIndex = cell.dCellEnd[gridHash];
        for(uint j = startIndex; j < endIndex; ++j){
            //if(j == i) continue;
 
            float3 pos1 = make_float3(cell.dSortedPos[j]);
 
            float3 rij = pos0-pos1;
            float r = length(rij);
 
            if(r <= h){
                float q = h*h-r*r;
                dens += params.Mass*params.Wpoly6*q*q*q;
            }
        }
    }
 
    return dens;
}
 
/*!
 * 粒子密度計算（kernel函式）
 * @param[out] newDens 粒子密度
 * @param[out] newPres 粒子壓力
 * @param[in]  cell 粒子網格資料
 */
__global__
void sphCalDensity(float* newDens, float* newPres, rxParticleCell cell)
{
    // 粒子檢索
    uint index = __mul24(blockIdx.x,blockDim.x)+threadIdx.x;
    if(index >= cell.uNumParticles) return;    
    
    float3 pos = make_float3(cell.dSortedPos[index]);    // 粒子位置
    float h = params.EffectiveRadius;
 
    // 粒子周圍的網格
    int3 grid_pos0, grid_pos1;
    grid_pos0 = calcGridPos(pos-make_float3(h));
    grid_pos1 = calcGridPos(pos+make_float3(h));
 
    // 包含周圍網格在內的鄰域搜尋，密度計算
    float dens = 0.0f;
    for(int z = grid_pos0.z; z <= grid_pos1.z; ++z){
        for(int y = grid_pos0.y; y <= grid_pos1.y; ++y){
            for(int x = grid_pos0.x; x <= grid_pos1.x; ++x){
                int3 n_grid_pos = make_int3(x, y, z);
                dens += calDensityCell(n_grid_pos, index, pos, cell);
            }
        }
    }
 
    // 密度值儲存
    uint oIdx = cell.dSortedIndex[index];
    newDens[oIdx] = dens;
}