cas test clusters

python/scontact.py

@@ -114,8 +114,8 @@ class ParticleProblem :
         self._p = c_void_p(self._lib.particleProblemNew())
         if self._p == None :
             raise NameError("Cannot create particle problem.")
-
         bndtype =[('material',np.int32),('tag',np.int32)]
+        self._particletype = np.dtype([('r',np.float64),('body',np.int64),('x',np.float64,dim)])
         self._disktype = np.dtype(bndtype+[('x',np.float64,dim),('v',np.float64,dim),('r',np.float64)])
         self._segmenttype = np.dtype(bndtype+[('p',np.float64,(2,dim)),('v',np.float64,(dim))])
         self._triangletype = np.dtype(bndtype+[('p',np.float64,(3,dim)),('v',np.float64,(dim))])
@@ -131,21 +131,32 @@ class ParticleProblem :
     def volume(self):
         """Returns the list of particle volumes."""
         if self._dim == 2 :
-            return np.pi * self._get_matrix("Particle", 2)[:, [0]]**2
+            return np.pi * self._get_array("Particle", self._particletype)["r"].reshape(-1,1)**2
         else :
-            return 4./3.*np.pi * self._get_matrix("Particle", 2)[:, [0]]**3
+            return 4./3.*np.pi * self._get_array("Particle", self._particletype)["r"].reshape(-1,1)**3
 
     def r(self) :
         """Returns the list of particle radii."""
-        return self._get_matrix("Particle", 2)[:, 0][:,None]
+        return self._get_array("Particle", self._particletype)["r"].reshape(-1,1)
 
     def mass(self):
-        """Returns the list of particle masses."""
-        return self._get_matrix("Particle", 2)[:, 1][:,None]
+        """Returns the list of body masses."""
+        return self._get_matrix("Body", 2)[:, 1][:,None]
 
     def position(self):
         return self._get_matrix("Position", self._dim)
 
+    def particle_position(self):
+        ppos = np.zeros((self.n_particles(),self._dim))
+        self._lib.getParticlePosition(self._p,_np2c(ppos))
+        return ppos
+
+    def particle_velocity(self):
+        pvel = np.zeros((self.n_particles(),self._dim))
+        self._lib.getParticleVelocity(self._p,_np2c(pvel))
+        return pvel
+
+
     def velocity(self):
         return self._get_matrix("Velocity", self._dim)
 
@@ -154,7 +165,7 @@ class ParticleProblem :
         if self._friction_enabled and self._rotation_enabled :
             return(self._get_matrix("Omega", d))
         else :
-            return np.zeros((self.n_particles(),d))
+            return np.zeros((self.n_bodies(),d))
 
     def save_state(self) :
         self._saved_velocity = np.copy(self.velocity())
@@ -202,6 +213,12 @@ class ParticleProblem :
         self._lib.particleProblemNParticle.restype = c_size_t
         return self._lib.particleProblemNParticle(self._p)
 
+    def n_bodies(self) :
+        """Returns the number of bodies."""
+        self._lib.particleProblemNBodies.restype = c_size_t
+        return self._lib.particleProblemNBodies(self._p)    
+
+
     def mu(self) :
         """Returns the matrix of the friction coefficients between materials."""
         return self._get_matrix("Mu", self._lib.particleProblemNMaterial(self._p))
@@ -345,16 +362,16 @@ class ParticleProblem :
         odir -- Directory in which to write the file
         i -- Number of the fiel to write
         t -- Time at which the simulation is 
-        """
-        v = self.velocity()
-        omega = self.omega()
-        x = self.position()
-        material = self._get_idx("ParticleMaterial").reshape(-1,1)
+        """      
+        v = self.particle_velocity()
+        x = self.particle_position()
+        material = self._get_idx("BodyMaterial").reshape(-1,1)
+        bodyId = self._get_array("Particle", self._particletype)["body"].reshape(-1,1)
         forced = self._get_idx("ForcedFlag").reshape(-1,1)
         if self._dim == 2 :
             v = np.insert(v,self._dim,0,axis=1)
             x = np.insert(x,self._dim,0,axis=1)
-        data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material),("ForcedFlag",forced),("Contacts",self.contact_forces())]
+        data = [("Radius",self.r()),("Body",bodyId),("Material",material[bodyId]),("Velocity",v),("ForcedFlag",forced[bodyId])]
         nmat = self._lib.particleProblemNMaterial(self._p)
         self._lib.particleProblemGetMaterialTagName.restype = c_char_p
         tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
@@ -570,8 +587,10 @@ class ParticleProblem :
         self._lib.particleProblemAddBoundaryPeriodicTriangle.restype = c_size_t
         return self._lib.particleProblemAddBoundaryPeriodicTriangle(self._p, self._coord_type(*x0), self._coord_type(*x1),  self._coord_type(*x2), c_int(etag))
 
-    def add_particle(self, x, r, m, tag="default",forced = 0):
-        """Adds a particle in the particle problem.
+
+
+    def add_body(self, positions,radii,masses,material="default",forced=0):
+        """Adds a body in the particle problem.
 
         Keyword arguments:
         x -- Tuple to set the centre particle position
@@ -580,32 +599,30 @@ class ParticleProblem :
         tag -- Particle material
         forced -- Flag indicating whether the particle is forced or not
         """
-        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_double(m), tag.encode(), c_int(forced))
+        idd = self.n_bodies()
+        if positions.shape[0] != radii.shape[0] or positions.shape[0]!=masses.shape[0] or masses.shape[0]!=radii.shape[0]: 
+          raise NameError("Wrong lengths of particle data vectors.")
+
+        x = np.sum(positions*masses.reshape((-1,1)),axis=0)/np.sum(masses)
+        I = np.sum(masses.reshape((-1,1))*radii.reshape((-1,1))**2/2 + np.sum((positions-x)**2,axis=1)*masses.reshape((-1,1)))
+        self._lib.particleProblemAddBody(self._p, self._coord_type(*x), c_double(np.sum(masses)), c_double(I),material.encode(),c_int(forced))
+        for i in range(positions.shape[0]):
+          self.add_particle(positions[i,:]-x,radii[i],idd)
+
+
+    def add_particle(self, x, r, body=0):
+        """Adds a particle in the particle problem.
 
-    def add_particles(self, x, r, m, v=None, omega=None, tag="default", forced=None, contact_forces=None):
-        """Adds particles in the particle problem.
         Keyword arguments:
-        x -- Array with the particles's centers position
-        r -- Array with the particles's radii
-        m -- Array with the particles's masses
-        v -- Array with the particles's velocities
-        omega -- Array with the particles's angular velocities
+        x -- Tuple to set the centre particle position
+        r -- Particle radius
+        m -- Particle mass
         tag -- Particle material
-        forced -- Array of flags indicating whether the particles are forced or not
-        contact_forces -- Array of contact forces between the particles
+        forced -- Flag indicating whether the particle is forced or not
         """
-        n_particles = len(m)
-        tags = np.ones(n_particles) * self.get_tag_id(tag)
-        if omega is None :
-         omega = np.zeros((1 if self._dim == 2 else 3)*n_particles)
-        if v is None :
-         v = np.zeros(self._dim*n_particles)
-        if forced is None :
-         forced = np.zeros(n_particles)
-        if contact_forces is None :
-         contact_forces = np.zeros(self._dim*n_particles)
-        self._lib.particleProblemAddParticles(self._p, c_int(n_particles), _np2c(x), _np2c(m), _np2c(r),_np2c(v), _np2c(omega), _np2c(tags,dtype=np.int32), _np2c(forced,dtype=np.int32), _np2c(contact_forces))
+        self._lib.particleProblemAddParticle(self._p, self._coord_type(*x), c_double(r), c_int(body))
 
+    
     def set_use_queue(self, use_queue=1):
       """Enables the use of the queue if 1 and disables it if 0."""
       self._lib.particleProblemSetUseQueue(self._p, c_int(use_queue))
@@ -628,11 +645,11 @@ class ParticleProblem :
       """Restores the saved contacts from the previous time step."""
       self._lib.particleProblemRestoreContacts(self._p)
 
-    def remove_particles_flag(self,flag) :
+    def remove_bodies_flag(self,flag) :
       """Removes particles based on given flag."""
       if flag.shape != (self.n_particles(),) :
           raise NameError("size of flag array should be the number of particles")
-      self._lib.particleProblemRemoveParticles(self._p, _np2c(flag,np.int32))
+      self._lib.particleProblemRemoveBodies(self._p, _np2c(flag,np.int32))
 
 
     def load_msh_boundaries(self, filename, tags=None, shift=None,material="default") :

scontact/scontact.h

@@ -33,7 +33,8 @@ void particleProblemLoad(ParticleProblem *p, const char *filename);
 void particleProblemWrite(const ParticleProblem *p, const char *filename);
 int particleProblemIterate(ParticleProblem *p, double alert, double dt, double tol, int maxit,int force_motion, double *externalForces);
 int particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
-void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m, const char *material, int forced);
+void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, int body);
+void particleProblemAddBody(ParticleProblem *p, const double x[DIMENSION], double m, double I, const char *material, int forced);
 size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag, const char *material);
 size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag, const char *material);
 size_t particleProblemAddBoundaryPeriodicSegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const int etag);
@@ -51,6 +52,7 @@ double *particleProblemTriangle(ParticleProblem *p);
 double *particleProblemSegment(ParticleProblem *p);
 double *particleProblemPeriodicSegment(ParticleProblem *p);
 double *particleProblemParticle(ParticleProblem *p);
+double *particleProblemBody(ParticleProblem *p);
 double *particleProblemVelocity(ParticleProblem *p);
 double *particleProblemPosition(ParticleProblem *p);
 double *particleProblemContactForces(ParticleProblem *p);
@@ -69,16 +71,14 @@ int *particleProblemTriangleTag(ParticleProblem *p);
 int *particleProblemTriangleMaterial(ParticleProblem *p);
 double *particleProblemTriangle(ParticleProblem *p);
 double *particleProblemPeriodicTriangle(ParticleProblem *p);
-
 #endif
-#if FRICTION_ENABLED
-void particleProblemSetFrictionCoefficient(ParticleProblem *p, double mu, const char *mat0, const char *mat1);
-double *particleProblemMu(ParticleProblem *p);
 #if ROTATION_ENABLED
 double *particleProblemOmega(ParticleProblem *p);
 #endif
+#if FRICTION_ENABLED
+void particleProblemSetFrictionCoefficient(ParticleProblem *p, double mu, const char *mat0, const char *mat1);
+double *particleProblemMu(ParticleProblem *p);
 #endif
 int particleProblemNTag(const ParticleProblem *p);
-void particleProblemRemoveParticles(ParticleProblem *p, const int *remove_flag);
 int particleProblemNMaterial(const ParticleProblem *p);
 #endif

testcases/bodies/test.py

+from migflow import scontact
+import os
+import numpy as np
+
+os.makedirs("output",exist_ok=True)
+
+p = scontact.ParticleProblem(2,friction_enabled=True,rotation_enabled=True)
+a = 1
+g = np.array([0,-9.81])
+r = 0.05
+rho = 1000
+p.add_boundary_segment([-a,-a],[-a,a],"bnd",material="Steel")
+p.add_boundary_segment([-a,a],[a,a],"bnd",material="Steel")
+p.add_boundary_segment([a,a],[a,-a],"bnd",material="Steel")
+p.add_boundary_segment([a,-a],[-a,-a],"bnd",material="Steel")
+R = np.array([r,r,r])
+x = np.arange(-15*r,a-2*r,3*r)
+y = np.arange(-10*r,a-2*r,3*r)
+#for i in range(len(x)):
+ #for j in range(len(y)):
+    #p.add_body(np.array([[x[i],y[j]]]),R,np.pi*R**2*rho,material="Sand")
+
+
+p.add_body(np.array([[0,5*r],[1.5*r,6.5*r],[3*r,8*r]]),R,np.pi*R**2*rho,material="Sand",forced=0)
+#p.add_body(np.array([[1.5*r,9*r],[3*r,9*r]]),R,np.pi*R**2*rho,material="Sand",forced=0)
+#p.add_body(np.array([[0,0*r],[1.5*r,0*r]]),R,np.pi*R**2*rho,material="Sand",forced=0)
+#p.add_body(np.array([[0,-2*r],[1.5*r,-2*r]]),R,np.pi*R**2*rho,material="Sand",forced=0)
+#p.add_body(np.array([[0,-4*r],[1.5*r,-4*r]]),R,np.pi*R**2*rho,material="Sand",forced=1)
+#R = np.array([r])
+#p.add_body(np.array([[-14*r,-15*r]]),R,np.pi*R**2*rho,material="Sand",forced=0)
+#p.add_body(np.array([[1.5*r,9*r]]),R,np.pi*R**2*rho,material="Sand",forced=0)
+#p.add_body(np.array([[0,4*r]]),R,np.pi*R**2*rho,material="Sand",forced=1)
+#p.add_body(np.array([[0,-4*r]]),R,np.pi*R**2*rho,material="Sand",forced=1)
+#p.add_body(np.array([[0,9*r]]),R,np.pi*R**Sand*2,rho="material",forced=0)
+
+tend=0.592
+p.set_friction_coefficient(.2,"Sand","Sand")
+p.set_friction_coefficient(.5,"Sand","Steel")
+i = 0
+f = g*(np.pi*p.r()**2*rho)
+p.write_vtk("output",0,0)
+dt = 0.001
+outf=1
+t = 0
+while t<tend :
+    print(i)
+    print(p.omega())
+    t += dt
+    p.iterate(dt,f,tol=1e-8,force_motion=1)
+    if i %outf == 0 :
+         ioutput = int(i/outf) + 1
+         p.write_vtk("output", ioutput, t)
+    i += 1

validation/inclined2d/inclined2d.py

@@ -49,25 +49,34 @@ class Inclined2d(unittest.TestCase) :
         p = scontact.ParticleProblem(2,True,True)
         p.add_boundary_segment((-10,-10),(10,-10),"bottom")
         #p.add_boundary_disk([0,0],0,"fake")
-        p.add_particle((-10+r,-10+r), r, r**2 * np.pi * rhop)
+        R = np.array([r])
+        p.add_body(np.array([[-10+r,-10+r]]),R,np.pi*R**2*rhop)
         p.set_friction_coefficient(0.5)
         # Initial time and iteration
         t         =  0
         ii        =  0
+        i=0
+        outf = 10
         #
         # COMPUTATION LOOP
         #
         vit = []
         ref = np.array([1,2,3,4,5,6,7,8,9,10])*0.55278468635
-        G = np.zeros_like(p.velocity())
-        G[:,1] = g[1]*p.mass()[:,0]*cos(radians(25))
-        G[:,0] = -g[1]*p.mass()[:,0]*sin(radians(25))
+        G = np.zeros((p.n_particles(),2))
+        
+        G[:,1] = g[1]*(np.pi*p.r()**2*rhop)*cos(radians(25))
+        G[:,0] = -g[1]*(np.pi*p.r()**2*rhop)*sin(radians(25))
         while t < tEnd :
-            p.iterate(dt, G, tol=1e-8)
+            p.iterate(dt, G, tol=1e-8,force_motion=1)
             t += dt
             ii += 1
             if ii%200==0:
               vit.append(np.sqrt(np.sum(p.velocity()**2)))
+            if i %outf == 0 :
+              ioutput = int(i/outf) + 1
+              p.write_vtk("output", ioutput, t)
+            i+=1
             print("%i : %.5g/%.5g" % (ii, t, tEnd))
+        print(vit-ref)
         error =  max(abs(ref - np.array(vit))/ref)
-        self.assertLess(error,.5/100., "error is too large in inclined2d")
+        self.assertLess(error,.01/100., "error is too large in inclined2d")