saving modif in depot

python/scontact.py

@@ -660,9 +660,9 @@ class ParticleProblem :
 
     def remove_bodies_flag(self,flag) :
       """Removes particles based on given flag."""
-      if flag.shape != (self.n_particles(),) :
+      if flag.shape != (self.n_bodies(),) :
           raise NameError("size of flag array should be the number of particles")
-      self._lib.particle_roblem_remove_bodies(self._p, _np2c(flag,np.int32))
+      self._lib.particle_problem_remove_bodies(self._p, _np2c(flag,np.int32))
 
 
     def load_msh_boundaries(self, filename, tags=None, shift=None,material="default") :

testcases/depot-2d/depot.py

@@ -17,6 +17,25 @@
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with this program (see COPYING and COPYING.LESSER files).  If not,
+# MigFlow - Copyright (C) <2010-2020>
+# <Universite catholique de Louvain (UCL), Belgium
+#  Universite de Montpellier, France>
+#
+# List of the contributors to the development of MigFlow: see AUTHORS file.
+# Description and complete License: see LICENSE file.
+#
+# This program (MigFlow) is free software:
+# you can redistribute it and/or modify it under the terms of the GNU Lesser General
+# Public License as published by the Free Software Foundation, either version
+# 3 of the License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program (see COPYING and COPYING.LESSER files).  If not,
 # see <http://www.gnu.org/licenses/>.
 
 #!/usr/bin/env python
@@ -31,8 +50,9 @@ from migflow import time_integration
 
 import numpy as np
 import os
-import time
 import shutil
+import matplotlib.pyplot as plt
+import time
 import random
 
 def genInitialPosition(filename, r, H, ly, lx, rhop) :
@@ -48,45 +68,51 @@ def genInitialPosition(filename, r, H, ly, lx, rhop) :
     """
     # Particles structure builder
     p = scontact.ParticleProblem(2)
-    #p2 = scontact.ParticleProblem(2)
-    # Load mesh.msh file specifying physical boundaries names
-    #p2.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
-    #for i in range(p2.segments().shape[0]):
-    #  p.add_boundary_segment((0,0),(1,1),0)
-    #p.segments()[:] = p2.segments()[:]
-    #print(p2.segments()[:],p.segments()[:])
-    p.add_boundary_segment((-0.2,-0.3),(-0.2,0.3),0)
-    p.add_boundary_segment((-0.2,0.3),(0.2,0.3),0)
-    p.add_boundary_segment((0.2,0.3),(0.2,-0.3),0)
-    p.add_boundary_segment((0.2,-0.3),(-0.2,-0.3),0)
-    #Definition of the points where the particles are located
-    x = np.arange(-lx/2+r-1e-8, lx/2-r+1e-8, 2.05*r)
-    y = np.arange(H/2-r, H/2-ly+r, -2.05*r)
-    x, y = np.meshgrid(x, y)
-    x = x.flat
-    y = y.flat
-    #Add a grain at each centre position
-    for i in range(len(x)) :
-        p.add_particle(np.array([x[i], y[i]]), r, r**2 * np.pi * rhop)
-    #p.add_particle_body(np.array([[0,0],[6*r,0],[12*r,0]]),np.array([3*r,3*r,3*r]),np.pi*rhop*np.array([3*r,3*r,3*r])**2)
+
+    p.add_boundary_segment((-0.07,-0.21),(-0.07,0.23/5),0,"Steel")
+    p.add_boundary_segment((-0.07,0.23/5),(0.07,0.23/5),0,"Steel")
+    p.add_boundary_segment((0.07,0.23/5),(0.07,-0.21),0,"Steel")
+    p.add_boundary_segment((0.07,-0.21),(-0.07,-0.21),0,"Steel")
+    R = np.ones((8))*r
+    angle = 5
+    Dx = 2*r*np.cos(np.deg2rad(angle))
+    Dy = 2*r*np.sin(np.deg2rad(angle))
+    pos = np.array([[s*Dx,s*Dy] for s in range(8)])
+    material = np.array(["Sand" for i in range(8)],dtype=object)
+    pos2 = np.copy(pos)
+    R2 = np.copy(R)
+    mat2 = np.copy(material)
+    for i in range (1):
+      pos = np.concatenate((pos,pos2+np.array([-(i+1)*Dy,(i+1)*Dx])))
+      R = np.concatenate((R,R2))
+      material = np.concatenate((material,mat2))
+    #R = np.concatenate((R,np.array([r,r,r,r,r,r])))
+    #material = np.concatenate((material,np.array(["Sand","Sand","Sand","Sand","Sand","Sand"],dtype=object)))
+    #pos3 = np.array([[s*Dx,s*Dy+0.013] for s in [0,15]])
+    #for i in range (3): 
+    #  pos = np.concatenate((pos,pos3+np.array([-(i+1)*Dy,(i+1)*Dx])))
+    #for i in range(1,4):
+    #  for j in range(1,15):
+    #    material[16*i + j] = "ignore"
+    p.add_particle_body(pos,R,np.pi*rhop*R**2,material)
     p.write_vtk(filename,0,0)
 
 # Define output directory
-outputdir = "output"
+outputdir = "outputLowDragHighI*"
 if not os.path.isdir(outputdir) :
     os.makedirs(outputdir)
 
 # Physical parameters
 g = np.array([0,-9.81])                         # gravity
-r = 1.5e-3                                        # particles radius
-rhop = 1500                                     # particles density
+r = 2.15e-3/4                                   # particles radius
+rhop = 1840                                     # particles density
 rho = 1000                                      # fluid density
 nu = 1e-6                                       # kinematic viscosity
 
 # Numerical parameters
-outf = 10                                        # number of iterations between output files
-dt = 1e-3                                     # time step
-tEnd = 100                                      # final time
+outf = 1                                        # number of iterations between output files
+dt = 1e-3                                       # time step
+tEnd = 2.5                                        # final time
 
 # Geometrical parameters
 ly = 1e-1                                       # particles area height
@@ -101,17 +127,18 @@ genInitialPosition(outputdir, r, H, ly, lx, rhop)
 p = scontact.ParticleProblem(2)
 p.read_vtk(outputdir,0)
 
-
+p.set_friction_coefficient(0.,"Sand","Steel")
 # Initial time and iteration
 t         =  0
 ii        =  0
-
+p.body_invert_mass()[4] = 1/(rhop*6.92e-3*1.25*2.15e-3)
+p.body_invert_inertia()[4] = 12*p.body_invert_mass()[4]/(6.92e-3**2+(1.25*2.15e-3)**2)
 #
 # FLUID PROBLEM
 #
 fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
 # Set the mesh geometry for the fluid computation
-fluid.load_msh("mesh.msh")
+fluid.load_msh("rect.msh")
 fluid.set_wall_boundary("Bottom")
 fluid.set_wall_boundary("Lateral")
 fluid.set_wall_boundary("Top",pressure=0)
@@ -124,8 +151,60 @@ tic = time.time()
 #
 # COMPUTATION LOOP
 #
+
+xdata = []
+ydata1 = []
+ydata2 = []
+ydata3 = []
+ydata4 = []
+fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2,2)
+
+
+ax1.plot(np.linspace(0,6,1000),np.zeros(1000))
+ax1.set_xlim(0, 6)
+ax1.set_ylim(-np.pi, np.pi)
+ax1.set_title('Theta')
+line1, = ax1.plot(xdata, ydata1, 'r-')
+
+ax2.plot(np.linspace(0,6,1000),np.zeros(1000))
+ax2.set_xlim(0, 6)
+ax2.set_ylim(-30, 30)
+ax2.set_title('Omega')
+line2, = ax2.plot(xdata, ydata2, 'c-')
+
+ax3.plot(np.linspace(0,6,1000),np.zeros(1000))
+ax3.set_xlim(0, 6)
+ax3.set_ylim(-0.1, 0.1)
+ax3.set_title('Vx')
+line3, = ax3.plot(xdata, ydata3, 'k-')
+
+ax4.plot(np.linspace(0,6,1000),np.zeros(1000))
+ax4.set_xlim(0, 6)
+ax4.set_ylim(-0.2, 0.1)
+ax4.set_title('Vy')
+line4, = ax4.plot(xdata, ydata4, 'g-')
+
 while t < tEnd :
-    time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=g*p.r()**2*rhop*np.pi)
+    #xdata.append(t)
+    #ydata1.append(p.body_theta()[4][0])
+    #line1.set_xdata(xdata)
+    #line1.set_ydata(ydata1)
+
+    #ydata2.append(p.body_omega()[4][0])
+    #line2.set_xdata(xdata)
+    #line2.set_ydata(ydata2)
+
+    #ydata3.append(p.body_velocity()[4][0])
+    #line3.set_xdata(xdata)
+    #line3.set_ydata(ydata3)
+
+    #ydata4.append(p.body_velocity()[4][1])
+    #line4.set_xdata(xdata)
+    #line4.set_ydata(ydata4)
+
+    #plt.draw()
+    #plt.pause(1e-21)
+    time_integration.iterate(fluid, p, dt, min_nsub=1, external_particles_forces=g*p.r()**2*rhop*np.pi)
     t += dt
 
     # Output files writting
@@ -135,3 +214,6 @@ while t < tEnd :
         fluid.write_vtk(outputdir, ioutput, t)
     ii += 1
     print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
+    print("\n\nCluster velocity : %f\n\n"%(np.linalg.norm(p.body_velocity()[4,:])))
+
+plt.show()

testcases/depot-2d/mesh.geo

-L = 0.2;
+L = 0.2;
 H = 0.3;
 y = 0;
 lc = 0.01;
@@ -11,16 +11,16 @@ Line(1) = {1, 2};
 Line(2) = {2, 3};
 Line(3) = {3, 4};
 Line(4) = {4, 1};
-Point(5) = {0,0.52,0};
+Point(5) = {0,0.52,lc};
 
 Line Loop(1) = {1:4};
 
 Plane Surface(1) = {1};
+
 Physical Line("Bottom") = {2};
 Physical Line("Lateral") = {1,3};
 Physical Line("Top") = {4};
 Physical Surface("Domain") = {1};
 Physical Point("PtFix") = {1};
-
 Mesh.CharacteristicLengthFromPoints = 1;
 Mesh.MshFileVersion = 2;