diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py index 105acce134ae21069e9ecd1278e541516943827b..8278862d4927e49f0e972bedc974b7fd76237d87 100644 --- a/examples/SodShock/makeIC.py +++ b/examples/SodShock/makeIC.py @@ -41,11 +41,30 @@ glass1 = h5py.File("../Glass/glass_200000.hdf5") pos1 = glass1["/PartType1/Coordinates"][:,:] pos1 = pos1 / 20. # Particles are in [0:0.5, 0:0.5, 0:0.5] +toRemove = array(0) +for i in range(size(pos1)/3-1): + if pos1[i,0] == pos1[i+1,0]: + if pos1[i,1] == pos1[i+1,1]: + if pos1[i,2] == pos1[i+1,2]: + toRemove = append(toRemove, i) + +pos1 = delete(pos1, toRemove, 0) + + #Read in high density glass glass2 = h5py.File("../Glass/glass_50000.hdf5") pos2 = glass2["/PartType1/Coordinates"][:,:] pos2 = pos2 / 20. # Particles are in [0:0.5, 0:0.5, 0:0.5] +toRemove = array(0) +for i in range(size(pos2)/3-1): + if pos2[i,0] == pos2[i+1,0]: + if pos2[i,1] == pos2[i+1,1]: + if pos2[i,2] == pos2[i+1,2]: + toRemove = append(toRemove, i) +pos2 = delete(pos2, toRemove, 0) + + #Generate high density region rho1 = size(pos1) / size(pos2) coord1 = append(pos1, pos1 + [0, 0.5, 0], 0)