Self Avoiding Random Cammina in un reticolo 3d usando l'algoritmo pivot in python

Question

Stavo lavorando su un problema degli ultimi giorni, È legato alla creazione di un auto che evita le passeggiate casuali usando l'algoritmo pivot e quindi per implementare un altro codice che pone inclusioni sferiche in un reticolo 3d. Ho scritto due codici: un codice genera le coordinate delle sfere usando l'algoritmo Self avoiding random walk e quindi un altro codice usa le coordinate generate dal primo programma e quindi crea un reticolo 3d in abaqus con inclusioni sferiche in esso.

Il risultato che viene generato dopo aver eseguito i codici:

Ora la sezione ingrandita (segnato in riquadro rosso)

Problema: Le sfere generati si fondono con l'altro, Come evitare questo, ho finito le idee. Ho solo bisogno di qualche direzione o algoritmo per lavorare su.

Il codice che ho scritto è il seguente: codice 1: genera le coordinate delle sfere

# Functions of the code: 1) This code generates the coordinates for the spherical fillers using self avoiding random walk which was implemented by pivot algorithm 
# Algorithm of the code: 1)Prepare an initial configuration of a N steps walk on lattice(equivalent N monomer chain) 
#      2)Randomly pick a site along the chain as pivot site 
#      3)Randomly pick a side(right to the pivot site or left to it), the chain on this side is used for the next step. 
#      4)Randomly apply a rotate operation on the part of the chain we choose at the above step. 
#      5)After the rotation, check the overlap between the rotated part of the chain and the rest part of the chain. 
#      6)Accept the new configuration if there is no overlap and restart from 2th step. 
#      7)Reject the configuration and repeat from 2th step if there are overlaps. 
################################################################################################################################################################ 
# Modules to be imported are below 
import numpy as np 
import timeit 
from scipy.spatial.distance import cdist 
import math 
import random 
import sys 
import ast 

# define a dot product function used for the rotate operation 
def v_dot(a):return lambda b: np.dot(a,b) 

def distance(x, y): #The Euclidean Distance to Spread the spheres initiallly 
    if len(x) != len(y): 
     raise ValueError, "vectors must be same length" 
    sum = 0 
    for i in range(len(x)): 
     sum += (x[i]-y[i])**2 
    return math.sqrt(sum) 

x=1/math.sqrt(2) 
class lattice_SAW: # This is the class which creates the self avoiding random walk coordinates 
    def __init__(self,N,l0): 
     self.N = N #No of spheres 
     self.l0 = l0 #distance between spheres 
     # initial configuration. Usually we just use a straight chain as inital configuration 
     self.init_state = np.dstack((np.arange(N),np.zeros(N),np.zeros(N)))[0] #initially set all the coordinates to zeros 
     self.state = self.init_state.copy() 

     # define a rotation matrix 
     # 9 possible rotations: 3 axes * 3 possible rotate angles(45,90,135,180,225,270,315) 
     self.rotate_matrix = np.array([[[1,0,0],[0,0,-1],[0,1,0]],[[1,0,0],[0,-1,0],[0,0,-1]] 
            ,[[1,0,0],[0,0,1],[0,-1,0]],[[0,0,1],[0,1,0],[-1,0,0]] 
            ,[[-1,0,0],[0,1,0],[0,0,-1]],[[0,0,-1],[0,1,0],[-1,0,0]] 
            ,[[0,-1,0],[1,0,0],[0,0,1]],[[-1,0,0],[0,-1,0],[0,0,1]] 
            ,[[0,1,0],[-1,0,0],[0,0,1]],[[x,-x,0],[x,x,0],[0,0,1]] 
            ,[[1,0,0],[0,x,-x],[0,x,x]] 
            ,[[x,0,x],[0,1,0],[-x,0,x]],[[-x,-x,0],[x,-x,0],[0,0,1]] 
            ,[[1,0,0],[0,-x,-x],[0,x,-x]],[[-x,0,x],[0,1,0],[-x,0,-x]] 
            ,[[-x,x,0],[-x,-x,0],[0,0,1]],[[1,0,0],[0,-x,x],[0,-x,-x]] 
            ,[[-x,0,-x],[0,1,0],[x,0,-x]],[[x,x,0],[-x,x,0],[0,0,1]] 
            ,[[1,0,0],[0,x,x],[0,-x,x]],[[x,0,-x],[0,1,0],[x,0,x]]]) 

    # define pivot algorithm process where t is the number of successful steps 
    def walk(self,t): #this definitions helps to start walking in 3d 
     acpt = 0 
     # while loop until the number of successful step up to t 
     while acpt <= t: 
      pick_pivot = np.random.randint(1,self.N-1)  # pick a pivot site 
      pick_side = np.random.choice([-1,1])   # pick a side 
      if pick_side == 1: 
       old_chain = self.state[0:pick_pivot+1] 
       temp_chain = self.state[pick_pivot+1:] 
      else: 
       old_chain = self.state[pick_pivot:] # variable to store the coordinates of the old chain 
       temp_chain = self.state[0:pick_pivot]# for the temp chain 
      # pick a symmetry operator 
      symtry_oprtr = self.rotate_matrix[np.random.randint(len(self.rotate_matrix))]  
      # new chain after symmetry operator 
      new_chain = np.apply_along_axis(v_dot(symtry_oprtr),1,temp_chain - self.state[pick_pivot]) + self.state[pick_pivot] 
      # use cdist function of scipy package to calculate the pair-pair distance between old_chain and new_chain 
      overlap = cdist(new_chain,old_chain) #compare the old chain and the new chain to check if the new chain overlaps on any previous postion 

      overlap = overlap.flatten() # just to combine the coordinates in a list which will help to check the overlap 

      # determinte whether the new state is accepted or rejected 
      if len(np.nonzero(overlap)[0]) != len(overlap): 
       continue 
      else: 
       if pick_side == 1: 
        self.state = np.concatenate((old_chain,new_chain),axis=0) 
       elif pick_side == -1: 
        self.state = np.concatenate((new_chain,old_chain),axis=0) 
       acpt += 1 

     # place the center of mass of the chain on the origin, so then we can translate these coordinates to the initial spread spheres 
     self.state = self.l0*(self.state - np.int_(np.mean(self.state,axis=0))) 
#now write the coordinates of the spheres in a file 
sys.stdout = open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions.txt", "w") 
n = 30 
#text_file.write("Number of Spheres: " + str(n) + "\n") 

#Radius of one sphere: is calculated based on the volume fraction here it 2 percent 
r = (3*0.40)/(4*math.pi*n) 
#print r 
r = r**(1.0/3.0) 
#print "Sphere Radius is " + str(r) 
sphereList = [] # list to maintain track of the Spheres using their positions 
sphereInstancesList = [] # to maintain track of the instances which will be used to cut and or translate the spheres 
sphere_pos=[] 

#Create n instances of the sphere 
#After creating n instances of the sphere we trie to form cluster around each instance of the sphere 
for i in range(1, n+1): 
    InstanceName = 'Sphere_' + str(i) # creating a name for the instance of the sphere 
    #print InstanceName 
    #text_file.write(InstanceName) 

    #Maximum tries to distribute sphere 
    maxTries = 10000000 

    while len(sphereList) < i: 
     maxTries -= 1 
     if maxTries < 1: 
      print "Maximum Distribution tries exceded. Error! Restart the Script!" 
      break; 

     # Make sure Spheres dont cut cube sides: this will place the spheres well inside the cube so that they does'nt touch the sides of the cube 
     # This is done to avoid the periodic boundary condition: later in the next versions it will be used 
     vecPosition = [(2*r)+(random.random()*(10.0-r-r-r)),(2*r)+(random.random()*(10.0-r-r-r)),(2*r)+(random.random()*(10.0-r-r-r))] 

     sphere_pos.append(vecPosition) 
     for pos in sphereList: 
      if distance(pos, vecPosition) < 2*r: # checking whether the spheres collide or not 
       break 
     else: 
      sphereList.append(vecPosition) 
      print vecPosition 
      #text_file.write(str(vecPosition) + "\n") 

cluster_Size=[10,12,14,16,18,20,22,24,26,28,30] #list to give the random number of spheres which forms a cluster 
for i in range(1,n+1): 
    Number_of_Spheres_Clustered=random.choice(cluster_Size) #selecting randomly from the list cluster_Size 
    radius_sphr=2*r #Distance between centers of the spheres 
    pivot_steps=1000 # for walking the max number of steps 
    chain = lattice_SAW(Number_of_Spheres_Clustered,radius_sphr) #initializing the object 
    chain.walk(pivot_steps) # calling the walk function to walk in the 3d space 
    co_ordinates=chain.state # copying the coordinates into a new variable for processing and converting the coordinates into lists 
    for c in range(len(co_ordinates)): 
     temp_cds=co_ordinates[c] 
     temp_cds=list(temp_cds)   
     for k in range(len(temp_cds)): 
      temp_cds[k]=temp_cds[k]+sphere_pos[i-1][k] 
     #text_file.write(str(temp_cds) + "\n") 
     print temp_cds #sys.stdout redirected into the file which stores the coordinates as lists 
sys.stdout.flush() 

f2=open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions.txt", "r") 
remove_check=[] 
for line in f2: 
    temp_check=ast.literal_eval(line) 
    if (temp_check[0]>10 or temp_check[0]<-r or temp_check[1]>10 or temp_check[1]<-r or temp_check[2]>10 or temp_check[2]<-r): 
     remove_check.append(str(temp_check)) 
f2.close() 
flag=0 
f2=open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions.txt", "r") 
f3=open("C:\Users\haris\Desktop\CAME_MINE\HIWI\Codes\Temp\Sphere_Positions_corrected.txt", "w") 
for line in f2: 
    line=line.strip() 
    if any(line in s for s in remove_check): 
     flag=flag+1 
    else: 
     f3.write(line+'\n') 
f3.close() 
f2.close() 

L'altro codice non sarebbe necessaria perché non c'è calcolo geometria nel secondo codice. Qualsiasi aiuto o una direzione di come evitare la collisione delle sfere è molto utile, Grazie a tutti

Answer 1

Per ospitare non si intersecano sfere con giri di 45.135.225,315 mila (in realtà, solo 45 e 315 sono problemi), non vi resta che riduci le dimensioni delle sfere. Prendi 3 centri sfera consecutivi, con una svolta di 45 gradi al centro. Nel piano contenente i 3 punti, che fa un triangolo isoscele, con angolo al centro di 45 gradi:

noti che i cerchi di fondo (sfere) si sovrappongono. Per evitare questo, si compatta il raggio moltiplicando per un fattore di 0,76: