The Sod shock IC script now explicitly removes particle pairs located at exactly the same position.

Former-commit-id: 5bc712f89d8cf024e193288a87b4266157d06c67

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)