reload segments and mu

python/scontact.py

@@ -374,7 +374,8 @@ class ParticleProblem :
                 ("Velocity",v),
                 ("RelativePosition",xrel)]
         tags = list (s.encode() for s in self.id2material)
-        VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags))])
+        mu = self.mu()
+        VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags)),(b"FrictionCoefficients",mu.reshape([-1,1]))])
 
         segments = self.segments()
         triangles = self.triangles()
@@ -398,8 +399,10 @@ class ParticleProblem :
         offsets = np.cumsum(np.repeat([2,3],[segments.shape[0],triangles.shape[0]]),dtype=np.int32)
         tags = np.array(np.hstack([segments["tag"],triangles["tag"]]))
         materials = np.array(np.hstack([segments["material"],triangles["material"]]))
+        bodies = np.array(np.hstack([segments["body"],triangles["body"]]))
         tagnames = list(s.encode() for s in self.id2tag)
-        VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,data=[("RelativePosition",xrel)],cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
+        cell_data = [("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1])),("Body",bodies.reshape([-1,1]))]
+        VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,data=[("RelativePosition",xrel)],cell_data=cell_data,field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
         periodicEntity = self.periodic_entity()
         periodicSegments = self.periodic_segments()
         periodicTriangles = self.periodic_triangles()
@@ -447,6 +450,8 @@ class ParticleProblem :
         x,_,data,cdata,fdata = VTK.read(dirname+"/particles_%05d.vtp"%i)
         self.id2material = list(s.decode() for s in VTK.string_array_decode(fdata["MaterialNames"]))
         self.material2id = dict((name,i) for i,name in enumerate(self.id2material))
+        self._lib.particle_problem_allocate_material(self._p, c_int(len(self.material2id)))
+        self.mu()[:] = fdata["FrictionCoefficients"].reshape(self.mu().shape)
         self._lib.particle_problem_allocate_particles(self._p, c_size_t(x.shape[0]))
         particles = self._get_array("particle",self._particletype)
         particles["r"][:] = data["Radius"][:,0]
@@ -462,8 +467,8 @@ class ParticleProblem :
         self._lib.particle_problem_allocate_segments(self._p, c_size_t(segidx.shape[0]))
         segxidx = np.column_stack([el[offsets[segidx]], el[offsets[segidx]+1]])
         segments = self.segments()
-        segments["p"] = x[el[segxidx],:self._dim]
-        segments["v"] = 0
+        segments["x"][:] = d["RelativePosition"][segxidx,:self._dim]
+        segments["body"][:] = cdata["Body"][:,0]
         segments["material"] = cdata["Material"][segidx][:,0]
         segments["tag"] = cdata["Tag"][segidx][:,0]
         if self._dim == 3:
@@ -503,7 +508,7 @@ class ParticleProblem :
         ks = ("particle_particle","particle_segment","particle_triangle")
         v = np.vstack(list([fdata[k] for k in ks]))
         o = np.vstack([fdata[k+"_idx"] for k in ks])
-        t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
+        t = np.repeat([0,2,3],[fdata[k].shape[0] for k in ks])
         self._lib.particle_problem_allocate_contacts(self._p, c_size_t(t.size))
         self.contacts()["reaction"] = v
         self.contacts()["o"] = o