merge and clean friction

python/scontact.py

@@ -42,6 +42,7 @@ signal.signal(signal.SIGINT, signal.SIG_DFL)
 
 dir_path = os.path.dirname(os.path.realpath(__file__))
 lib2 = np.ctypeslib.load_library("libscontact2",dir_path)
+lib2f = np.ctypeslib.load_library("libscontact2f",dir_path)
 lib3 = np.ctypeslib.load_library("libscontact3",dir_path)
 
 is_64bits = sys.maxsize > 2**32
@@ -71,7 +72,7 @@ class ParticleProblem :
         buf = (size*c_int).from_address(ptr)
         return np.ctypeslib.as_array(buf)
 
-    def __init__(self, dim) :
+    def __init__(self, dim, friction_enabled=False) :
         """Build the particles contanier structure.
 
         Keyword arguments:
@@ -82,8 +83,9 @@ class ParticleProblem :
         NameError -- if C builder cannot be found
         """
         self._dim = dim
+        self._friction_enabled = friction_enabled
         if dim == 2 :
-            self._lib = lib2
+            self._lib = lib2f if friction_enabled else lib2
             self._coord_type = c_double*2
         elif dim == 3 :
             self._lib = lib3
@@ -125,13 +127,14 @@ class ParticleProblem :
 
     def omega(self):
         """Return the list of particle angular velocity."""
+        assert self._friction_enabled
         return self._get_matrix("Omega", 1)
 
     def disk_tag(self):
         return self._get_idx("DiskTag")
 
     def segments(self):
-        return self._get_matrix("Segment",3*self._dim+1)
+        return self._get_matrix("Segment",3*self._dim+(1 if self._friction_enabled else 0))
 
     def segment_tag(self):
         return self._get_idx("SegmentTag")
@@ -163,12 +166,15 @@ class ParticleProblem :
             triangles[triangleTag == tag, 9:12] = v
 
     def set_tangent_boundary_velocity(self, tag, v):
+        assert self._friction_enabled
         self._lib.particleProblemSetBoundaryTangentVelocity(self._p,tag.encode(),c_double(v))
     
     def set_friction_relaxation(self, relaxation):
+        assert self._friction_enabled
         self._lib.particleProblemSetFrictionRelaxation(self._p,c_double(relaxation))
 
     def add_friction(self, static) :
+        assert self._friction_enabled
         arg_type = c_double*len(static)
         self._lib.particleProblemAddFriction(self._p, arg_type(*static),len(static))
 
@@ -187,7 +193,10 @@ class ParticleProblem :
         """Write output files for post-visualization."""
 
         v = self.velocity()
-        omega = self.omega()
+        if self._friction_enabled :
+            omega = self.omega()
+        else :
+            omega = np.zeros((v.shape[0],1))
         x = self.position()
         tag = self._get_idx("ParticleTag").reshape(-1,1)
         if self._dim == 2 :
@@ -196,7 +205,7 @@ class ParticleProblem :
         data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Tag",tag)]
         VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp")
 
-        xdisk = self._get_matrix("Disk",2*self._dim+2)[:,:self._dim]
+        xdisk = self._get_matrix("Disk",2*self._dim+(2 if self._friction_enabled else 1))[:,:self._dim]
         nd = xdisk.shape[0]
         xsegment = self.segments()[:,:self._dim*2].reshape([-1,self._dim])
         ns = xsegment.shape[0]/2
@@ -229,8 +238,8 @@ class ParticleProblem :
     def read_vtk(self,dirname,i) :
         """Read an existing output file to set the particle data."""
         x,_,d,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
-        print(d["Omega"])
         self._lib.particleProblemClearParticles(self._p)
+
         self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
                 _np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(d["Tag"]))
 

scontact/CMakeLists.txt

@@ -28,6 +28,11 @@ include_directories(. ../tools)
 add_library(scontact2 SHARED ${SRC})
 target_compile_definitions(scontact2 PUBLIC "-DDIMENSION=2")
 target_link_libraries(scontact2 mbtools2)
+
+add_library(scontact2f SHARED ${SRC})
+target_compile_definitions(scontact2f PUBLIC -DDIMENSION=2 -DFRICTION_ENABLED=1)
+target_link_libraries(scontact2f mbtools2)
+
 add_library(scontact3 SHARED ${SRC})
 target_link_libraries(scontact3 mbtools3)
 target_compile_definitions(scontact3 PUBLIC "-DDIMENSION=3")

scontact/scontact.h

@@ -35,7 +35,6 @@ void particleProblemIterate(ParticleProblem *p, double alert, double dt, double
 void particleProblemSolve(ParticleProblem *p, double alert, double dt, double tol, int maxit);
 double particleProblemMaxDt(const ParticleProblem *p, double alert);
 void particleProblemAddParticle(ParticleProblem *p, const double x[DIMENSION], double r, double m, int tag);
-void particleProblemAddFriction(ParticleProblem *p, double *stat, int len);
 size_t particleProblemAddBoundaryDisk(ParticleProblem *p, const double x0[DIMENSION], double r, const char *tag);
 size_t particleProblemAddBoundarySegment(ParticleProblem *p, const double x0[DIMENSION], const double x1[DIMENSION], const char *tag);
 size_t particleProblemGetTagId(ParticleProblem *p, const char *tag);
@@ -46,12 +45,10 @@ void particleProblemAddBoundaryTriangle(ParticleProblem *p, const double x0[DIME
 double *particleProblemDisk(ParticleProblem *p);
 double *particleProblemSegment(ParticleProblem *p);
 double *particleProblemParticle(ParticleProblem *p);
-double *particleProblemOmega(ParticleProblem *p);
 double *particleProblemVelocity(ParticleProblem *p);
 double *particleProblemPosition(ParticleProblem *p);
 double *particleProblemExternalForces(ParticleProblem *p);
 void particleProblemSetUseQueue(ParticleProblem *p, int use_queue);
-void particleProblemSetFrictionRelaxation(ParticleProblem *p, double friction_relaxation);
 unsigned long int particleProblemNParticle(ParticleProblem *p);
 double *particleProblemRadius(ParticleProblem *p);
 int *particleProblemDiskTag(ParticleProblem *p);
@@ -62,3 +59,9 @@ double *particleProblemTriangle(ParticleProblem *p);
 #endif
 void ParticleToMesh(size_t nParticles, double *particles, int nElements, double *elements, int *elid, double *xi);
 #endif
+
+#ifdef FRICTION_ENABLED
+void particleProblemAddFriction(ParticleProblem *p, double *stat, int len);
+double *particleProblemOmega(ParticleProblem *p);
+void particleProblemSetFrictionRelaxation(ParticleProblem *p, double friction_relaxation);
+#endif

testcases/mill-2d/mesh.geo

-rout = 0.1;
+rout = 0.11;
 lcout = rout/20;
 
 Point(1) = {0, 0, 0};

testcases/mill-2d/mill.py

@@ -83,7 +83,7 @@ if use_lmgc90 :
     lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
     p = lmgc90Interface.ParticleProblem(2)
 else :
-    p = scontact.ParticleProblem(2)
+    p = scontact.ParticleProblem(2,True)
     p.read_vtk(outputdir,0)
 
 outf = 5                                # number of iterations between output files