def get_dummy_particles():
    x, y = numpy.mgrid[-5 * dx : box_length + 5 * dx + 1e-10 : dx, -5 * dx : box_height + 5 * dx + 1e-10 : dx]

    xd, yd = x.ravel(), y.ravel()

    md = numpy.ones_like(xd) * m
    hd = numpy.ones_like(xd) * h

    rhod = numpy.ones_like(xd) * ro
    cd = numpy.ones_like(xd) * co
    pd = numpy.zeros_like(xd)

    dummy_fluid = base.get_particle_array(name="dummy_fluid", type=Fluid, x=xd, y=yd, h=hd, rho=rhod, c=cd, p=pd)

    # remove indices within the square

    indices = []

    np = dummy_fluid.get_number_of_particles()
    x, y = dummy_fluid.get("x", "y")

    for i in range(np):
        if -dx / 2 <= x[i] <= box_length + dx / 2:
            if -dx / 2 <= y[i] <= box_height + dx / 2:
                indices.append(i)

    to_remove = base.LongArray(len(indices))
    to_remove.set_data(numpy.array(indices))

    dummy_fluid.remove_particles(to_remove)

    return dummy_fluid

def get_boundary_particles(**kwargs):

    # left boundary
    x = numpy.ones(50)
    for i in range(50):
        x[i] = -0.6 - (i + 1) * dxl

    m = numpy.ones_like(x) * dxl
    h = numpy.ones_like(x) * 2 * dxr
    rho = numpy.ones_like(x)
    u = numpy.zeros_like(x)
    e = numpy.ones_like(x) * 2.5

    p = (0.4) * rho * e
    cs = numpy.sqrt(1.4 * p / rho)

    left = base.get_particle_array(name="left", type=Boundary, x=x, m=m, h=h, rho=rho, u=u, e=e, cs=cs, p=p)

    # right boundary
    for i in range(50):
        x[i] = 0.6 + (i + 1) * dxr

    m = numpy.ones_like(x) * dxl
    h = numpy.ones_like(x) * 2 * dxr
    rho = numpy.ones_like(x) * 0.25
    u = numpy.zeros_like(x)
    e = numpy.ones_like(x) * 1.795

    p = (0.4) * rho * e
    # cs = numpy.sqrt(0.4*e)
    cs = numpy.sqrt(1.4 * p / rho)

    right = base.get_particle_array(name="right", type=Boundary, x=x, m=m, h=h, rho=rho, u=u, e=e, cs=cs, p=p)

    return [left, right]

    def setUp(self):
        """ The setup consists of points randomly distributed in a
        cube [-1,1] X [-1.1] X [1-,1]

        The smoothing length for the points is proportional to the
        number of particles.

        """

        self.cl_precision = cl_precision = "single"
        self.np = np = 1<<14

        self.x = x = random.random(np) * 2.0 - 1.0
        self.y = y = random.random(np) * 2.0 - 1.0
        self.z = z = random.random(np) * 2.0 - 1.0

        vol_per_particle = numpy.power(8.0/np, 1.0/3.0)

        self.h = h = numpy.ones_like(x) * vol_per_particle

        self.cy_pa = base.get_particle_array(name='test', x=x, y=y, z=z, h=h)

        self.cl_pa = base.get_particle_array(
            name="test", cl_precision=self.cl_precision,
            x=x, y=y, z=z, h=h)

        # the scale factor for the cell sizes
        self.kernel_scale_factor = kernel_scale_factor = 2.0

        self._setup()

 def t(self):
     ret = {}
     da = DoubleArray()
     pa = ParticleArray()
     kernel = kernels.CubicSplineKernel(3)
     get_time = time.time
     for N in Ns:
         x = numpy.arange(N)
         z = y = numpy.zeros(N)
         mu = m = rho = numpy.ones(N)
         h = 2*m
         da = DoubleArray(N)
         da2 = DoubleArray(N)
         da.set_data(z)
         da2.set_data(z)
         pa = get_particle_array(x=x, y=y, z=z, h=h, mu=mu, rho=rho, m=m, tmp=z,
                                 tx=z, ty=m, tz=z, nx=m, ny=z, nz=z, u=z, v=z, w=z,
                                 ubar=z, vbar=z, wbar=z, q=m)
         pb = get_particle_array(x=x+0.1**0.5, y=y, z=z, h=h, mu=mu, rho=rho, m=m, tmp=z,
                                 tx=m, ty=z, tz=z, nx=z, ny=m, nz=z, u=z, v=z, w=z,
                                 ubar=z, vbar=z, wbar=z, q=m)
         particles = Particles(arrays=[pa, pb])
         
         func = func_getter.get_func(pa, pb)
         calc = SPHCalc(particles, [pa], pb, kernel, [func], ['tmp']*func.num_outputs)
         print cls.__name__
         t = get_time()
         calc.sph('tmp', 'tmp', 'tmp')
         t = get_time() - t
         
         nam = '%s'%(cls.__name__)
         ret[nam +' /%d'%(N)] = t/N
     return ret

    def setUp(self):
        """ Setup for SPHOperationTestCase

        Setup:
        ------
        Create two particle arrays, one Fluid and one Solid.
        Instantiate the class with a default function `SPHRho` with 
        various combinations of the from and on types and check for the 
        filtering of the arrays.

        """

        x = numpy.linspace(0,1,11)
        h = numpy.ones_like(x) * 2 * (x[1] - x[0])

        #create the fluid particle array

        self.fluid = base.get_particle_array(name='fluid', type=Fluid, x=x,h=h)

        #create the solid particle array
        
        self.solid = base.get_particle_array(name="solid", type=Solid, x=x,h=h)

        #create the particles
        
        self.particles = particles = base.Particles(arrays=[self.fluid,
                                                            self.solid])

        #set the kernel
        
        self.kernel = base.CubicSplineKernel()

    def setUp(self):
        """ The setup consists of two fluid particle arrays, each
        having one particle. The fluids are acted upon by an external
        vector force and gravity.

        Comparison is made with the PySPH integration of the system.
        
        """

        x1 = numpy.array([-0.5])
        y1 = numpy.array([1.0])

        x2 = numpy.array([0.5])
        y2 = numpy.array([1.0])

        tmpx1 = numpy.ones_like(x1)
        tmpx2 = numpy.ones_like(x2)

        self.f1 = base.get_particle_array(name="fluid1", x=x1, y=y1, tmpx=tmpx1)
        self.f2 = base.get_particle_array(name="fluid2", x=x2, y=y2, tmpx=tmpx2)

        self.particles = base.Particles(arrays=[self.f1, self.f2])
        self.kernel = kernel = base.CubicSplineKernel(dim=2)

        gravity = solver.SPHIntegration(
            sph.GravityForce.withargs(gy=-10.0), on_types=[Fluid], updates=["u", "v"], id="gravity"
        )

        force = solver.SPHIntegration(
            sph.GravityForce.withargs(gx=-10.0), on_types=[Fluid], updates=["u", "v"], id="force"
        )

        position = solver.SPHIntegration(sph.PositionStepping, on_types=[Fluid], updates=["x", "y"], id="step")

        gravity.calc_type = sph.CLCalc
        force.calc_type = sph.CLCalc
        position.calc_type = sph.CLCalc

        gravity_calcs = gravity.get_calcs(self.particles, kernel)
        force_calcs = force.get_calcs(self.particles, kernel)
        position_calcs = position.get_calcs(self.particles, kernel)

        self.calcs = calcs = []
        calcs.extend(gravity_calcs)
        calcs.extend(force_calcs)
        calcs.extend(position_calcs)

        self.integrator = CLIntegrator(self.particles, calcs)

        self.ctx = ctx = cl.create_some_context()
        self.integrator.setup_integrator(ctx)
        self.queue = calcs[0].queue

        self.dt = 0.1
        self.nsteps = 10

def get_boundary_particles():
    """ Get the particles corresponding to the dam and fluids """

    # get the tank
    xt1, yt1 = geom.create_2D_tank(x1=0, y1=0,
                                   x2=tank_length, y2=tank_height,
                                   dx=dx)

    xt2,  yt2 = geom.create_2D_tank(x1=-dx/2, y1=-dx/2,
                                    x2=tank_length + dx/2, y2=tank_height+dx/2,
                                    dx=dx)

    x = numpy.concatenate( (xt1, xt2) )
    y = numpy.concatenate( (yt1, yt2) )

    h = numpy.ones_like(x) * h0
    m = numpy.ones_like(x) * ro*dx*dx*0.5
    rho = numpy.ones_like(x) * ro
    cs = numpy.ones_like(x) * co

    tank = base.get_particle_array(cl_precision="single", name="tank",
                                   type=Solid, x=x,y=y,m=m,rho=rho,h=h,cs=cs)
    np = tank.get_number_of_particles()

    # create the gate
    y1 = numpy.arange(dx/2, tank_height+1e-4, dx/2)
    x1 = numpy.ones_like(y1)*(0.38-dx/2)

    y2 = numpy.arange(dx/2+dx/4, tank_height+1e-4, dx/2)
    x2 = numpy.ones_like(y2)*(0.38-dx)

    y3 = numpy.arange(dx/2, tank_height+1e-4, dx/2)
    x3 = numpy.ones_like(y3)*(0.38-1.5*dx)

    x = numpy.concatenate( (x1, x2, x3) )
    y = numpy.concatenate( (y1, y2, y3) )

    h = numpy.ones_like(x) * h0
    m = numpy.ones_like(x) * 0.5 * dx/2 * dx/2 * ro
    rho = numpy.ones_like(x) * ro
    cs = numpy.ones_like(x) * co
    v = numpy.ones_like(x) * 1.5

    gate = base.get_particle_array(cl_precision="single", name="gate",
                                   x=x, y=y, m=m, rho=rho, h=h, cs=cs,
                                   v=v,
                                   type=Solid)

    np += gate.get_number_of_particles()
    print "Number of solid particles = %d"%(np)

    return [tank, gate]

    def setUp(self):

        self.np = 25

        self.x1 = x1 = numpy.array([-0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,]
                                   )
        
        self.y1 = y1 = numpy.array([-0.125, -0.125, -0.125, -0.125, -0.125,
                                    0.125, 0.125, 0.125, 0.125, 0.125,
                                    0.375, 0.375, 0.375, 0.375, 0.375,
                                    0.625, 0.625, 0.625, 0.625, 0.625,
                                    0.875, 0.875, 0.875, 0.875, 0.875]
                                   )

        self.z1 = z1 = numpy.zeros_like(x1)

        self.h1 = h1 = numpy.ones_like(x1) * 0.1

        self.x2 = x2 = numpy.array([-0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,
                                    -0.125, 0.125, 0.375, 0.625, 0.875,]
                                   )
        
        self.y2 = y2 = numpy.array([-0.125, -0.125, -0.125, -0.125, -0.125,
                                    0.125, 0.125, 0.125, 0.125, 0.125,
                                    0.375, 0.375, 0.375, 0.375, 0.375,
                                    0.625, 0.625, 0.625, 0.625, 0.625,
                                    0.875, 0.875, 0.875, 0.875, 0.875]
                                   )
       
        self.z2 = z2 = numpy.zeros_like(x2)
        
        self.h2 = h2 = numpy.ones_like(x2) * 0.15

        self.pa1 = pa1 = base.get_particle_array(name='test1',
                                                 x=x1, y=y1, z=z1, h=h1)

        self.pa2 = pa2 = base.get_particle_array(name='test2',
                                                 x=x2, y=y2, z=z2, h=h2)

        self.periodic_domain = periodic_domain = PeriodicDomain(xmin=-0.2,
                                                                xmax=1.0)

    def setUp(self):

        self.np = np = 100
        
        x = numpy.random.random(np)
        y = numpy.random.random(np)
        z = numpy.random.random(np)
        m = numpy.ones_like(x)

        cy_pa = base.get_particle_array(name="test", x=x, y=y, z=z, m=m)
        cl_pa = base.get_particle_array(name="test", cl_precision="double",
                                        x=x, y=y, z=z, m=m)

        cy_particles = base.Particles(
            [cy_pa,], locator_type=CYLoctor.NSquareNeighborLocator)

        cl_particles = base.CLParticles( [cl_pa,] )

        cy_solver = solver.Solver(
            dim=3, integrator_type=solver.EulerIntegrator)

        cl_solver = solver.Solver(
            dim=3, integrator_type=solver.EulerIntegrator)

        self.cy_solver = cy_solver
        self.cl_solver = cl_solver

        cy_solver.add_operation(solver.SPHIntegration(

            sph.NBodyForce.withargs(), on_types=[0], from_types=[0],
            updates=['u','v','w'], id='nbody_force')

                                )

        cy_solver.add_operation_step([0,])

        cl_solver.add_operation(solver.SPHIntegration(

            sph.NBodyForce.withargs(), on_types=[0], from_types=[0],
            updates=['u','v','w'], id='nbody_force')

                                )

        cl_solver.add_operation_step([0,])
        cl_solver.set_cl(True)
        
        cy_solver.setup(cy_particles)
        cl_solver.setup(cl_particles)

    def setUp(self):

        if not solver.HAS_CL:
            try:
                import nose.plugins.skip as skip
                reason = "PyOpenCL not installed!"
                raise skip.SkipTest(reason)

            except ImportError:
                pass
        
        self.np = np = 101
        self.x = x = numpy.linspace(0, 1, np)
        self.m = m = numpy.ones_like(x) * (x[1] - x[0])
        self.h = h = 2*self.m

        pa = base.get_particle_array(name="test", cl_precision="single",
                                     x=x, m=m, h=h)

        particles = base.CLParticles([pa,])
        kernel = base.CubicSplineKernel(dim=1)

        func = sph.GravityForce.withargs(gx=-1, gy=-1, gz=-1).get_func(pa,pa)

        self.calc = sph.CLCalc(particles, sources=[pa,], dest=pa,
                               updates=['u','v','w'], kernel=kernel,
                               funcs=[func,])

        if solver.HAS_CL:
            self.context = solver.create_some_context()

def get_fluid_particles():
    
    xf1, yf1 = geom.create_2D_filled_region(x1=dx, y1=dx,
                                            x2=fluid_column_width,
                                            y2=fluid_column_height,
                                            dx=dx)

    xf2, yf2 = geom.create_2D_filled_region(x1=dx/2, y1=dx/2,
                                            x2=fluid_column_width,
                                            y2=fluid_column_height,
                                            dx=dx)
    

    x = numpy.concatenate((xf1, xf2))
    y = numpy.concatenate((yf1, yf2))

    print 'Number of fluid particles: ', len(x)

    hf = numpy.ones_like(x) * h
    mf = numpy.ones_like(x) * dx * dy * ro * 0.5
    rhof = numpy.ones_like(x) * ro
    csf = numpy.ones_like(x) * co
    
    fluid = base.get_particle_array(cl_precision="single",
                                    name="fluid", type=Fluid,
                                    x=x, y=y, h=hf, m=mf, rho=rhof,
                                    cs=csf)

    return fluid

    def test_load_output(self):
        
        self.setup_solver()
        s = self.solver
        s.particles = self.particles
        d = tempfile.mkdtemp()
        s.output_directory = d
        s.detailed_output = True
        s.fname = 'temp_solver'
        s.dt = 0

        x = numpy.arange(10)
        pa = base.get_particle_array(name='pa', x=x)
        #pa = base.ParticleArray(name='pa')
        pa.add_property(dict(name='q', data=-x))
        s.particles.arrays[0] = pa
        
        old_props = {}
        for name,prop in s.particles.arrays[0].properties.iteritems():
            old_props[name] = prop.get_npy_array().copy()

        s.dump_output(s.dt)
        pa.q = pa.x = pa.z
        ret = s.load_output('?')
        self.assertEqual(ret, ["0"])
        s.load_output('0')

        try:
            for name,prop in s.particles.arrays[0].properties.iteritems():
                self.assertTrue(numpy.allclose(prop,old_props[name]),
                                msg='prop:%s\nold:%s, new:%s'%(name,old_props[name], pa.get(name)))
        finally:
            shutil.rmtree(d, True)

def get_fluid():
    """ Get the fluid particle array """

    x, y = numpy.mgrid[dx : box_length - 1e-10 : dx, dx : box_height - 1e-10 : dx]

    xf, yf = x.ravel(), y.ravel()

    mf = numpy.ones_like(xf) * m
    hf = numpy.ones_like(xf) * h

    rhof = numpy.ones_like(xf) * ro
    cf = numpy.ones_like(xf) * co
    pf = numpy.zeros_like(xf)

    fluid = base.get_particle_array(name="fluid", type=Fluid, x=xf, y=yf, h=hf, rho=rhof, c=cf, p=pf)

    # remove indices within the square

    indices = []

    np = fluid.get_number_of_particles()
    x, y = fluid.get("x", "y")

    for i in range(np):
        if 1.0 - dx / 2 <= x[i] <= 2.0 + dx / 2:
            if 2.0 - dx / 2 <= y[i] <= 3.0 + dx / 2:
                indices.append(i)

    to_remove = base.LongArray(len(indices))
    to_remove.set_data(numpy.array(indices))

    fluid.remove_particles(to_remove)

    return fluid

    def setUp(self):
        """ A Dummy fluid solver is created with the following operations

        (i)  -- Equation of State
        (ii) -- Density Rate
        (iii)-- Momentum Equation Pressure Gradient
        (iv) -- Momentum Equation Viscosity 
        
        """
        self.kernel = kernel = base.CubicSplineKernel(dim = 2)

        self.solver = s = solver.Solver(dim=2,
                                        integrator_type=solver.EulerIntegrator)

        s.default_kernel = kernel

        self.particles = base.Particles(arrays=[base.get_particle_array()])

        # Create some replacement operations
        
        self.visc = solver.SPHIntegration(

            sph.MorrisViscosity, on_types=[Fluids],
            from_types = [Fluids, Solids], updates=['u','v'], id='visc'
            
            )

        self.summation_density = solver.SPHOperation(

            sph.SPHRho, on_types=[Fluids, Solids],
            updates=['rho'], id='sd'

            )

def test_sph_calc():

    x = numpy.array([0,])
    y = numpy.array([0,])
    z = numpy.array([0,])
    h = numpy.ones_like(x)

    pa = base.get_particle_array(name="test", x=x, y=y, z=z,h=h)
    particles = base.Particles(arrays=[pa,])
    kernel = base.CubicSplineKernel(dim=1)

    vector_force1 = sph.VectorForce.withargs(force=base.Point(1,1,1))
    vector_force2 = sph.VectorForce.withargs(force=base.Point(1,1,1))

    func1 = vector_force1.get_func(pa,pa)
    func2 = vector_force2.get_func(pa,pa)

    calc = sph.SPHCalc(particles=particles, sources=[pa,pa], dest=pa,
                       kernel=kernel, funcs=[func1, func2],
                       updates=['u','v','w'], integrates=True)

    # evaluate the calc. Accelerations are stored in _tmpx, _tmpy and _tmpz

    calc.sph('_tmpx', '_tmpy', '_tmpz')

    tmpx, tmpy, tmpz = pa.get('_tmpx', '_tmpy', '_tmpz')

    # the acceleration should be 2 in each direction

    assert ( abs(tmpx[0] - 2.0) < 1e-16 )
    assert ( abs(tmpy[0] - 2.0) < 1e-16 )
    assert ( abs(tmpz[0] - 2.0) < 1e-16 )

    def setUp(self):

        # create a dummy particle array
        self.pa = pa = base.get_particle_array(name="test",
                                               cl_precision="single",
                                               type=1)

        # create a simple solver with one operation
        self.s = s = solver.Solver(3, solver.EulerIntegrator)
        
        # add the velocity gradient operation 
        s.add_operation(solver.SPHOperation(

            sph.VelocityGradient3D.withargs(), on_types=[1,],
            from_types=[1,], id="vgrad")

                             )
        # add the stress rate function
        s.add_operation(solver.SPHIntegration(

            sph.HookesDeviatoricStressRate3D.withargs(),
            on_types=[1,], id="stress_rate")

                             )
        particles = base.Particles(arrays=[pa,])
        s.setup(particles)

        self.integrator = s.integrator

def create_particles_2d(**kwargs):

    x, y = numpy.mgrid[0.25:0.75+1e-10:dx, 0.25:0.75+1e-10:dx]
    x = x.ravel()
    y = y.ravel()

    np = len(x)
    
    u = numpy.zeros_like(x)
    v = numpy.zeros_like(x)

    m = numpy.ones_like(x) * dx**2

    vol_per_particle = numpy.power(0.5**2/np ,1.0/2.0)
    radius = 2 * vol_per_particle

    print "Using smoothing length: ", radius

    h = numpy.ones_like(x) * radius

    fluid = base.get_particle_array(name="fluid", type=base.Fluid,
                                    x=x, y=y,
                                    u=u, v=v,
                                    m=m,
                                    h=h)

    print "Number of particles: ", fluid.get_number_of_particles()

    return [fluid,]

def create_particles_3d(**kwargs):

    x, y, z = numpy.mgrid[0.25:0.75+1e-10:dx,
                          0.25:0.75+1e-10:dx,
                          0.25:0.75+1e-10:dx]


    x = x.ravel()
    y = y.ravel()
    z = z.ravel()

    np = len(x)

    u = random.random(np) * 0
    v = random.random(np) * 0
    w = random.random(np) * 0

    m = numpy.ones_like(x) * dx**3

    vol_per_particle = numpy.power(0.5**3/np ,1.0/3.0)
    radius = 2 * vol_per_particle

    print "Using smoothing length: ", radius

    h = numpy.ones_like(x) * radius

    fluid = base.get_particle_array(name="fluid", type=base.Fluid,
                                    x=x, y=y, z=z,
                                    u=u, v=v, w=w,
                                    m=m,h=h)

    print "Number of particles: ", fluid.get_number_of_particles()

    return [fluid,]

    def setUp(self):
        """ The setup consists of a 2D square ([0,1] X [0,1]) with 25
        particles.


        """

        xc = numpy.arange(0,1.0, 0.2)
        x, y = numpy.meshgrid(xc,xc)

        self.x = x = x.ravel()
        self.y = y = y.ravel()
        self.h = h = numpy.ones_like(x) * 0.25

        dx = dy = 0.2
        self.dx = dx

        self.block_size = 0.5
        self.cell_size = 0.5

        self.pa = pa = base.get_particle_array(name="test", x=x, y=y, h=h)

        self.cm = cm = parallel.ParallelCellManager(arrays_to_bin=[pa,],
                                                    max_radius_scale=2.0,
                                                    dimension=2.0,
                                                    load_balancing=False,
                                                    initialize=False)

        self.block_000_indices = 0,1,2,5,6,7,10,11,12
        self.block_100_indices = 3,4,8,9,13,14
        self.block_010_indices = 15,16,17,20,21,22
        self.block_110_indices = 18,19,23,24

def standard_shock_tube_data(name="", type=0, cl_precision="double",
                             nl=320, nr=80, smoothing_length=None, **kwargs):
    """ Standard 400 particles shock tube problem """
    
    dxl = 0.6/nl
    dxr = dxl*4
    
    x = numpy.ones(nl+nr)
    x[:nl] = numpy.arange(-0.6, -dxl+1e-10, dxl)
    x[nl:] = numpy.arange(dxr, 0.6+1e-10, dxr)

    m = numpy.ones_like(x)*dxl
    h = numpy.ones_like(x)*2*dxr

    if smoothing_length:
        h = numpy.ones_like(x) * smoothing_length
    
    rho = numpy.ones_like(x)
    rho[nl:] = 0.25
    
    u = numpy.zeros_like(x)
    
    e = numpy.ones_like(x)
    e[:nl] = 2.5
    e[nl:] = 1.795

    p = 0.4*rho*e

    cs = numpy.sqrt(1.4*p/rho)

    idx = numpy.arange(nl+nr)
    
    return base.get_particle_array(name=name,x=x,m=m,h=h,rho=rho,p=p,e=e,
                                   cs=cs,type=type, idx=idx,
                                   cl_precision=cl_precision)