generateE.cu

Go to the documentation of this file.

#include "generateE.h"


__device__
real dhRomaDeviceE(real x, real h)
{
    real r = fabs(x)/h;
    
    if(r>1.5)
        return 0.0;
    else if(r>0.5 && r<=1.5)
        return 1.0/(6*h)*( 5.0 - 3.0*r - sqrt(-3.0*(1-r)*(1-r) + 1.0) );
    else
        return 1.0/(3*h)*( 1.0 + sqrt(-3.0*r*r + 1.0) );
} // dhRomaDeviceE


__device__
real deltaDeviceE(real x, real y, real h)
{
    return dhRomaDeviceE(x, h) * dhRomaDeviceE(y, h);
} // deltaDeviceE


namespace kernels
{

__global__
void generateE(int *ERows,  int *ECols,  real *EVals,
               int nx, int ny, real *x, real *y, real *dx,
               int totalPoints, real *xB, real *yB, int *I, int *J)
{
    int bodyIdx = threadIdx.x + blockIdx.x*blockDim.x;
    
    if(bodyIdx < totalPoints)
    {
        int Ib=I[bodyIdx],
            Jb=J[bodyIdx],
            EIdx  = bodyIdx*12,
            i, j;

        real Dx = dx[Ib];
    
        // uB = integral u * delta * dxdy = Ehat * u
        // E = Ehat * R^-1 => divide by Dx
        // E = Dx * delta
    
        // populate x-components
        for(j=Jb-1; j<=Jb+1; j++)
        {
            for(i=Ib-2; i<=Ib+1; i++)
            {
                ERows[EIdx] = bodyIdx;
                ECols[EIdx] = j*(nx-1) + i;
                EVals[EIdx] = Dx*deltaDeviceE(x[i+1]-xB[bodyIdx], 0.5*(y[j]+y[j+1])-yB[bodyIdx], Dx);
                EIdx++;
            }
        }

        // populate y-components
        for(j=Jb-2; j<=Jb+1; j++)
        {
            for(i=Ib-1; i<=Ib+1; i++)
            {
                ERows[EIdx+12*totalPoints-12] = bodyIdx + totalPoints;
                ECols[EIdx+12*totalPoints-12] = j*nx + i + (nx-1)*ny;
                EVals[EIdx+12*totalPoints-12] = Dx*deltaDeviceE(0.5*(x[i]+x[i+1])-xB[bodyIdx], y[j+1]-yB[bodyIdx], Dx);
                EIdx++;
            }
        }
    }
} // generateE

} // End of namespace kernels