tmirror particle grid for porosity estimation at frictionless walls - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
 (HTM) git clone git://src.adamsgaard.dk/sphere
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) LICENSE
       ---
 (DIR) commit b58442eb622184dad34d0120f2acc12f5d2e928f
 (DIR) parent f643ec5f6e9db8a8ea620826a84ae46b7f781f7e
 (HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
       Date:   Wed, 15 Jun 2016 14:45:02 -0700
       
       mirror particle grid for porosity estimation at frictionless walls
       
       Diffstat:
         M src/contactsearch.cuh               |      79 +++++++++++++++++++++++++++++++
         M src/darcy.cuh                       |       4 ++--
       
       2 files changed, 81 insertions(+), 2 deletions(-)
       ---
 (DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
       t@@ -80,6 +80,85 @@ __device__ int findDistMod(int3* targetCell, Float3* distmod)
                return 0;
        }
        
       +/* Copy of findDistMod, but modified for mirror BC at frictionless boundaries 
       +   during porosity esimation */
       +__device__ int findDistModPorosity(int3* targetCell, Float3* distmod)
       +{
       +    // Check whether x- and y boundaries are to be treated as periodic
       +    // 1: x- and y boundaries periodic
       +    // 2: x boundaries periodic
       +    if (devC_grid.periodic == 1) {
       +
       +        // Periodic x-boundary
       +        if (targetCell->x < 0) {
       +            //targetCell->x = devC_grid.num[0] - 1;
       +            targetCell->x += devC_grid.num[0];
       +            *distmod += MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
       +        }
       +        if (targetCell->x >= devC_grid.num[0]) {
       +            //targetCell->x = 0;
       +            targetCell->x -= devC_grid.num[0];
       +            *distmod -= MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
       +        }
       +
       +        // Periodic y-boundary
       +        if (targetCell->y < 0) {
       +            //targetCell->y = devC_grid.num[1] - 1;
       +            targetCell->y += devC_grid.num[1];
       +            *distmod += MAKE_FLOAT3(0.0, devC_grid.L[1], 0.0);
       +        }
       +        if (targetCell->y >= devC_grid.num[1]) {
       +            //targetCell->y = 0;
       +            targetCell->y -= devC_grid.num[1];
       +            *distmod -= MAKE_FLOAT3(0.0, devC_grid.L[1], 0.0);
       +        }
       +
       +
       +        // Only x-boundaries are periodic
       +    } else if (devC_grid.periodic == 2) {
       +
       +        // Periodic x-boundary
       +        if (targetCell->x < 0) {
       +            //targetCell->x = devC_grid.num[0] - 1;
       +            targetCell->x += devC_grid.num[0];
       +            *distmod += MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
       +        }
       +        if (targetCell->x >= devC_grid.num[0]) {
       +            //targetCell->x = 0;
       +            targetCell->x -= devC_grid.num[0];
       +            *distmod -= MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
       +        }
       +
       +        // Hande out-of grid cases on y-axis by mirroring the grid
       +        if (targetCell->y < 0) {
       +            targetCell->y += 2;
       +            //*distmod += MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
       +        }
       +        if (targetCell->y >= devC_grid.num[1]) {
       +            targetCell->y -= 2;
       +            //*distmod += MAKE_FLOAT3(devC_grid.L[0], 0.0, 0.0);
       +        }
       +
       +
       +        // No periodic boundaries
       +    } else {
       +
       +        // Don't modify targetCell or distmod.
       +
       +        // Hande out-of grid cases on x- and y-axes
       +        if (targetCell->x < 0 || targetCell->x >= devC_grid.num[0])
       +            return -1;
       +        if (targetCell->y < 0 || targetCell->y >= devC_grid.num[1])
       +            return -1;
       +    }
       +
       +    // Handle out-of-grid cases on z-axis
       +    if (targetCell->z < 0 || targetCell->z >= devC_grid.num[2])
       +        return -1;
       +    else
       +        return 0;
       +}
       +
        
        
        // Find overlaps between particle no. 'idx' and particles in cell 'gridpos'.
 (DIR) diff --git a/src/darcy.cuh b/src/darcy.cuh
       t@@ -398,7 +398,7 @@ __global__ void findDarcyPorositiesLinear(
        
                if (np > 0) {
        
       -            // Cell sphere center position
       +            // Cell center position
                    const Float3 X = MAKE_FLOAT3(
                            x*dx + 0.5*dx,
                            y*dy + 0.5*dy,
       t@@ -767,7 +767,7 @@ __global__ void findDarcyPorosities(
                                // Get distance modifier for interparticle
                                // vector, if it crosses a periodic boundary
                                distmod = MAKE_FLOAT3(0.0, 0.0, 0.0);
       -                        if (findDistMod(&targetCell, &distmod) != -1) {
       +                        if (findDistModPorosity(&targetCell, &distmod) != -1) {
        
                                    // Calculate linear cell ID
                                    cellID = targetCell.x