Added open boundary and automatic particle removal (#38)

dem/app/app.py

@@ -8,6 +8,7 @@
 from dem.app.silo         import *
 from dem.app.mill         import *
 from dem.app.circular     import *
+from dem.app.silo_open    import *
 
 # Declare factory
 app_factory = factory()
@@ -21,3 +22,4 @@
 app_factory.register("silo",        silo)
 app_factory.register("mill",        mill)
 app_factory.register("circular",    circular)
+app_factory.register("silo_open",   silo_open)

dem/app/base_app.py

@@ -51,6 +51,18 @@ def update(self):
         self.it += 1
         self.t  += self.dt
 
+    ### ************************************************
+    ### Check if particles are still in the domain
+    def check_particles(self, d):
+
+        lst = []
+
+        for i in range(self.p.np):
+            if (not d.is_in(self.p.x[i,:])):
+                lst.append(i)
+
+        self.p.delete(lst)
+
     ### ************************************************
     ### Plot
     def plot(self):

dem/app/silo_open.py

@@ -0,0 +1,127 @@
+# Generic imports
+import os
+import math
+import random
+
+# Custom imports
+from dem.app.base_app import *
+
+### ************************************************
+### Silo case with open boundary
+class silo_open(base_app):
+    ### ************************************************
+    ### Constructor
+    def __init__(self,
+                 name             = 'silo_open',
+                 t_max            = 3.5,
+                 dt               = 2.5e-5,
+                 angle            = 0.0,
+                 plot_freq        = 500,
+                 plot_show        = False,
+                 plot_trajectory  = False,
+                 plot_png         = True):
+        super().__init__()
+
+        self.name             = name
+        self.t_max            = t_max
+        self.dt               = dt
+        self.plot_freq        = plot_freq
+        self.plot_show        = plot_show
+        self.plot_trajectory  = plot_trajectory
+        self.plot_png         = plot_png
+
+        self.nt         = int(self.t_max/self.dt)
+        self.plot_it    = 0
+        self.check_freq = 10 # check particles every 10 iterations
+
+        self.n_row  = 25 # nb of particles on a row at start
+        self.n_col  = 20 # nb of particles on a col at start
+        self.radius = 0.025
+
+        self.p = particles(np          = self.n_row*self.n_col,
+                           nt          = self.nt,
+                           material    = "steel",
+                           radius      = self.radius,
+                           color       = "b",
+                           store       = False,
+                           search      = "nearest",
+                           rad_coeff   = 2.0)
+
+        self.p.e_wall[:] = 0.5
+        self.p.e_part[:] = 0.5
+
+        colors = np.array(['r', 'g', 'b', 'c', 'm', 'y', 'k'])
+        self.p.c = colors[np.random.randint(0,len(colors),size=self.p.np)]
+
+        self.d = domain_factory.create("rectangle",
+                                       x_min         = 0.0,
+                                       x_max         = 3.0,
+                                       y_min         = 0.0,
+                                       y_max         = 5.0,
+                                       angle         = 0.0,
+                                       material      = "steel",
+                                       open_boundary = [True, False, False, False])
+
+        self.o0 = domain_factory.create("rectangle",
+                                        x_min      =-1.0,
+                                        x_max      = 1.4,
+                                        y_min      = 2.5,
+                                        y_max      = 2.6,
+                                        angle      =-30.0,
+                                        plot_fill  = True,
+                                        material   = "steel")
+        self.o1 = domain_factory.create("rectangle",
+                                        x_min      = 1.6,
+                                        x_max      = 4.0,
+                                        y_min      = 2.5,
+                                        y_max      = 2.6,
+                                        angle      = 30.0,
+                                        plot_fill  = True,
+                                        material   = "steel")
+
+        self.d_lst = [self.d, self.o0, self.o1]
+
+        self.path = self.base_path+'/'+self.name
+        os.makedirs(self.path, exist_ok=True)
+
+        self.reset()
+
+    ### ************************************************
+    ### Reset app
+    def reset(self):
+
+        self.it = 0
+        self.t  = 0.0
+        sep  = 4.0*self.radius
+        mid  = 0.5*(self.d.x_min + self.d.x_max)
+        half = 0.5*self.n_row*(self.radius + 0.75*sep)
+
+        for i in range(self.n_row):
+            for j in range(self.n_col):
+                self.p.x[self.n_col*i+j,0] = mid - half + sep*i + self.radius*random.random()
+                self.p.x[self.n_col*i+j,1] = self.d.y_max - 2*sep - sep*j
+
+    ### ************************************************
+    ### Compute forces
+    def forces(self):
+
+        # Removing lost particles requires to recompute the
+        # nearest neighbor lists. Everytime the check is done,
+        # we force the recomputation
+        force_nearest = False
+        if (self.it%self.check_freq == 0):
+            self.check_particles(self.d)
+            force_nearest = True
+
+        self.p.reset_forces()
+        self.p.collisions(self.dt, force_nearest=force_nearest)
+        self.d.collisions(self.p, self.dt)
+        self.o0.collisions(self.p, self.dt)
+        self.o1.collisions(self.p, self.dt)
+        self.p.gravity(self.g)
+
+    ### ************************************************
+    ### Iteration printings (overrides base_app::printings)
+    def printings(self):
+
+        print("# it = "+str(self.it)+", t = {0:.3f}".format(self.t)+" / {0:.3f}".format(self.t_max)+", n_particles = "+str(self.p.np), end='\r')

dem/src/core/particles.py

@@ -86,17 +86,43 @@ def set_material(self, i, m):
     ### Compute maximal velocity
     def max_velocity(self):
 
-        return np.max(np.linalg.norm(self.v, axis=1))
+        if (len(self.v) == 0):
+            return 0.0
+        else:
+            return np.max(np.linalg.norm(self.v, axis=1))
 
     ### ************************************************
     ### Compute maximal radius
     def max_radius(self):
 
-        return np.max(self.r)
+        if (len(self.r) == 0):
+            return 0.0
+        else:
+            return np.max(self.r)
+
+    ### ************************************************
+    ### Remove a list of particles
+    def delete(self, lst):
+
+        self.np     -= len(lst)
+        self.m       = np.delete(self.m,       lst, axis=0)
+        self.r       = np.delete(self.r,       lst, axis=0)
+        self.x       = np.delete(self.x,       lst, axis=0)
+        self.d       = np.delete(self.d,       lst, axis=0)
+        self.v       = np.delete(self.v,       lst, axis=0)
+        self.a       = np.delete(self.a,       lst, axis=0)
+        self.f       = np.delete(self.f,       lst, axis=0)
+        self.e_wall  = np.delete(self.e_wall,  lst, axis=0)
+        self.mu_wall = np.delete(self.mu_wall, lst, axis=0)
+        self.e_part  = np.delete(self.e_part,  lst, axis=0)
+        self.mu_part = np.delete(self.mu_part, lst, axis=0)
+        self.Y       = np.delete(self.Y,       lst, axis=0)
+        self.G       = np.delete(self.G,       lst, axis=0)
+        self.c       = np.delete(np.array(self.c), lst, axis=0).tolist()
 
     ### ************************************************
     ### Compute collisions between particles
-    def collisions(self, dt):
+    def collisions(self, dt, force_nearest=False):
 
         if (self.search == "linear"):
 
@@ -117,7 +143,7 @@ def collisions(self, dt):
             # If list is not valid anymore
             v_max       = self.max_velocity()
             self.d_ngb += 2.0*v_max*dt
-            if (self.d_ngb > 0.99*self.m_rad):
+            if ((self.d_ngb > 0.99*self.m_rad) or force_nearest):
                 self.d_ngb = 0.0
                 self.ngbi, self.ngbj = list_ngbs(self.np, self.x,
                                                  self.r,  self.m_rad)

dem/src/domain/base.py

@@ -9,19 +9,7 @@ def __init__(self):
         pass
 
     ### ************************************************
-    ### Compute distance to particle
-    def distance(self, p):
-
-        raise NotImplementedError
-
-    ### ************************************************
-    ### Detect collision with particle
-    def collide(self, p):
-
-        raise NotImplementedError
-
-    ### ************************************************
-    ### Plot domain
-    def plot(self, ax):
+    ### Check if point is in domain
+    def is_in(self, pt):
 
         raise NotImplementedError

dem/src/domain/rectangle.py

@@ -17,13 +17,14 @@ class rectangle(base_domain):
     ### ************************************************
     ### Constructor
     def __init__(self,
-                 x_min     = 0.0,
-                 x_max     = 1.0,
-                 y_min     = 0.0,
-                 y_max     = 1.0,
-                 angle     = 0.0, # in degrees
-                 plot_fill = False,
-                 material  = "steel"):
+                 x_min         = 0.0,
+                 x_max         = 1.0,
+                 y_min         = 0.0,
+                 y_max         = 1.0,
+                 angle         = 0.0, # in degrees
+                 plot_fill     = False,
+                 material      = "steel",
+                 open_boundary = [False, False, False, False]):
 
         # Type
         self.type = "rectangle"
@@ -44,6 +45,9 @@ def __init__(self,
         # Material
         self.mat = material_factory.create(material)
 
+        # Open boundaries
+        self.open_boundary = np.array(open_boundary)
+
         # Ficticious parameters for domain
         self.r     = 1.0e8
         self.m     = 1.0e8
@@ -81,7 +85,7 @@ def collisions(self, p, dt):
 
         # Search for collisions linearly and compute forces
         linear_search(self.seg_pts, self.r, self.m,
-                      self.v, self.mat.Y, self.mat.G,
+                      self.v, self.mat.Y, self.mat.G, self.open_boundary,
                       p.x, p.r, p.m, p.v, p.e_wall, p.mu_wall,
                       p.Y, p.G, p.f, p.np, dt)
 
@@ -96,13 +100,34 @@ def rotate(self, p, pc, angle):
         dp[1] = (p[0]-pc[0])*sint + (p[1]-pc[1])*cost + pc[1]
         p[:]  = dp[:]
 
+    ### ************************************************
+    ### Check if point is in domain
+    ### pt is assumed to be an np array of size 2
+    def is_in(self, pm):
+
+        p1p2 = self.p2 - self.p1
+        p1pm =      pm - self.p1
+        p1p4 = self.p4 - self.p1
+
+        p1p2p1pm = np.dot(p1p2, p1pm)
+        p1p2p1p2 = np.dot(p1p2, p1p2)
+        p1p4p1pm = np.dot(p1p4, p1pm)
+        p1p4p1p4 = np.dot(p1p4, p1p4)
+
+        if ((p1p2p1pm > 0.0)      and
+            (p1p2p1p2 > p1p2p1pm) and
+            (p1p4p1pm > 0.0)      and
+            (p1p4p1p4 > p1p4p1pm)): return True
+
+        return False
+
 ### ************************************************
 ### Distance from rectangle domain to given coordinates
 ### Prefix d_ corresponds to domain
 ### Prefix p_ corresponds to particle
 @nb.njit(cache=True)
 def linear_search(d_pts, d_r, d_m,
-                  d_v, d_mat_Y, d_mat_G,
+                  d_v, d_mat_Y, d_mat_G, open_boundary,
                   p_x, p_r, p_m, p_v, p_e_wall, p_mu_wall,
                   p_Y, p_G, p_f, n, dt):
 
@@ -111,6 +136,10 @@ def linear_search(d_pts, d_r, d_m,
 
         # Loop on rectangle sides
         for j in range(4):
+
+            # If this boundary is open, loop
+            if (open_boundary[j]): continue
+
             dx, nrm = p_to_segment(p_x[i], d_pts[j,0], d_pts[j,1])
             dx      = dx - p_r[i]