remove contact basis

paraview_utils/MigFlow.xml

@@ -63,7 +63,7 @@
        name="Script"
        command="SetScript"
        number_of_elements="1"
-       default_values="import numpy as np

in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)

vals = []
vals_t = []
vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
vals_t.append(in_contacts.FieldData["particle_particle_t"])
if vals_s is not None:
    vals_s.append(in_contacts.FieldData["particle_particle_s"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
    disks = in_bnd.Cells[1:2*n_disks:2]
    vals.append(in_contacts.FieldData["particle_disk"])
    vals_t.append(in_contacts.FieldData["particle_disk_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_disk_s"])
    contacts = in_contacts.FieldData["particle_disk_idx"]
    points_d = np.ndarray((contacts.shape[0],2,3))
    points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
    points.append(points_d)

#particle-segments contacts

if n_segments :
    segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    vals.append(in_contacts.FieldData["particle_segment"])
    vals_t.append(in_contacts.FieldData["particle_segment_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_segment_s"])
    contacts = in_contacts.FieldData["particle_segment_idx"]
    points_s = np.ndarray((contacts.shape[0],2,3))
    points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
    s = in_bnd.Points[segments[contacts[:,0],:]]
    t = s[:,1,:]-s[:,0,:]
    t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
    d = points_s[:,0,:]-s[:,0,:]
    l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
    points_s[:,1,:] = s[:,0,:]+l[:,None]*t
    points.append(points_s)

#particle-triangle contacts
if n_triangles :
    triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
    vals.append(in_contacts.FieldData["particle_triangle"])
    vals_t.append(in_contacts.FieldData["particle_triangle_t"])
    if vals_s is not None:
        vals_s.append(in_contacts.FieldData["particle_triangle_s"])
    contacts = in_contacts.FieldData["particle_triangle_idx"]
    points_t = np.ndarray((contacts.shape[0],2,3))
    points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
    points.append(points_t)

#merge everything

points = np.vstack(points)
vals = np.hstack(vals)
vals_t = np.hstack(vals_t)
if vals_s is not None :
    vals_s = np.hstack(vals_s)
#generate tubes
t = points[:,0,:]-points[:,1,:]
t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
n1 = np.cross(t,ez)
n2 = np.cross(t,n1)
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals**(1./2)*1

opoints = points[:,None,:,:] \
    +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
    +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
if vals_s is not None :
    output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int)
for i in range(nf) :
    j = (i+1)%nf
    pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
+       default_values="import numpy as np

in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])

n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)

vals = []
points = []

#particle-particle contacts
vals.append(in_contacts.FieldData["particle_particle"])
contacts = in_contacts.FieldData["particle_particle_idx"]
points.append(in_particles.Points[contacts,:])
#particle-disk contacts
if n_disks :
    disks = in_bnd.Cells[1:2*n_disks:2]
    vals.append(in_contacts.FieldData["particle_disk"])
    contacts = in_contacts.FieldData["particle_disk_idx"]
    points_d = np.ndarray((contacts.shape[0],2,3))
    points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_d[:,1,:] = in_bnd.Points[disks[contacts[:,0]],:]
    points.append(points_d)

#particle-segments contacts

if n_segments :
    segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    vals.append(in_contacts.FieldData["particle_segment"])
    contacts = in_contacts.FieldData["particle_segment_idx"]
    points_s = np.ndarray((contacts.shape[0],2,3))
    points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
    s = in_bnd.Points[segments[contacts[:,0],:]]
    t = s[:,1,:]-s[:,0,:]
    t /= ((t[:,0]**2+t[:,1]**2+t[:,2]**2)**0.5)[:,None]
    d = points_s[:,0,:]-s[:,0,:]
    l = d[:,0]*t[:,0]+d[:,1]*t[:,1]+d[:,2]*t[:,2]
    points_s[:,1,:] = s[:,0,:]+l[:,None]*t
    points.append(points_s)

#particle-triangle contacts
if n_triangles :
    triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
    vals.append(in_contacts.FieldData["particle_triangle"])
    contacts = in_contacts.FieldData["particle_triangle_idx"]
    points_t = np.ndarray((contacts.shape[0],2,3))
    points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
    points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
    points.append(points_t)

#generate tubes
points = np.vstack(points)
b0 = points[:,0,:]-points[:,1,:]
b0 /= np.sqrt(b0[:,0,None]**2+b0[:,1,None]**2+b0[:,2,None]**2)
ez = np.where(np.abs(b0[:,1,None])>np.abs(b0[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
b1 = np.cross(b0,ez)
b2 = np.cross(b0,b1)
#merge everything
vals = np.vstack(vals)
if vals.shape[1] == 3:
    vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1]+vals[:,2]*b0[:,2])
    vals_t = np.hypot(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1]+vals[:,2]*b1[:,2],
                      vals[:,0]*b2[:,0]+vals[:,1]*b2[:,1]+vals[:,2]*b2[:,2])
else:
    vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1])
    vals_t = np.abs(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1])
alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
r = rf*vals_n**(1./2)*1

opoints = points[:,None,:,:] \
    +b1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
    +b2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
output.Points = opoints.reshape(-1,3)
output.PointData.append(np.repeat(vals_n,2*nf),"Reaction")
output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
n = points.shape[0]
pattern = np.ndarray([nf,4], np.int32)
for i in range(nf) :
    j = (i+1)%nf
    pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
types = np.full([n*nf],9,np.uint8)
locations = np.arange(0,5*n*nf,5,np.uint32)
cells = np.ndarray([n,nf,5],np.uint32)
cells[:,:,0] = 4
cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
output.SetCells(types, locations, cells.reshape([-1]))
"
        panel_visibility="advanced">
        <Hints>
         <Widget type="multi_line"/>

paraview_utils/contacts-extract.py

@@ -34,33 +34,24 @@ Properties = dict(nf = 5, rf = 2.5e-4)
 
 def RequestData():
     import numpy as np
+
     in_particles,in_bnd,in_periodic,in_contacts = list(inputs[0])
 
-    n_disks = np.count_nonzero(in_bnd.CellTypes == 1) + np.count_nonzero(in_periodic.CellTypes == 1)
-    n_segments  = (np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3) 
-                + np.count_nonzero(in_periodic.CellTypes == 7) + np.count_nonzero(in_periodic.CellTypes == 3))
-    n_triangles = (np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)
-                + np.count_nonzero(in_periodic.CellTypes == 13)+np.count_nonzero(in_periodic.CellTypes==5))
+    n_disks = np.count_nonzero(in_bnd.CellTypes == 1)
+    n_segments = np.count_nonzero(in_bnd.CellTypes == 7)+np.count_nonzero(in_bnd.CellTypes==3)
+    n_triangles = np.count_nonzero(in_bnd.CellTypes == 13)+np.count_nonzero(in_bnd.CellTypes==5)
 
     vals = []
-    vals_t = []
-    vals_s = [] if ("particle_particle_s" in in_contacts.FieldData.keys()) else None
     points = []
 
     #particle-particle contacts
     vals.append(in_contacts.FieldData["particle_particle"])
-    vals_t.append(in_contacts.FieldData["particle_particle_t"])
-    if vals_s is not None:
-        vals_s.append(in_contacts.FieldData["particle_particle_s"])
     contacts = in_contacts.FieldData["particle_particle_idx"]
     points.append(in_particles.Points[contacts,:])
     #particle-disk contacts
     if n_disks :
         disks = in_bnd.Cells[1:2*n_disks:2]
         vals.append(in_contacts.FieldData["particle_disk"])
-        vals_t.append(in_contacts.FieldData["particle_disk_t"])
-        if vals_s is not None:
-            vals_s.append(in_contacts.FieldData["particle_disk_s"])
         contacts = in_contacts.FieldData["particle_disk_idx"]
         points_d = np.ndarray((contacts.shape[0],2,3))
         points_d[:,0,:] = in_particles.Points[contacts[:,1],:]
@@ -72,9 +63,6 @@ def RequestData():
     if n_segments :
         segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
         vals.append(in_contacts.FieldData["particle_segment"])
-        vals_t.append(in_contacts.FieldData["particle_segment_t"])
-        if vals_s is not None:
-            vals_s.append(in_contacts.FieldData["particle_segment_s"])
         contacts = in_contacts.FieldData["particle_segment_idx"]
         points_s = np.ndarray((contacts.shape[0],2,3))
         points_s[:,0,:] = in_particles.Points[contacts[:,1],:]
@@ -90,47 +78,39 @@ def RequestData():
     if n_triangles :
         triangles = in_bnd.Cells[2*n_disks+3*n_segments:2*n_disks+3*n_segments+4*n_triangles].reshape([-1,4])[:,1:]
         vals.append(in_contacts.FieldData["particle_triangle"])
-        vals_t.append(in_contacts.FieldData["particle_triangle_t"])
-        if vals_s is not None:
-            vals_s.append(in_contacts.FieldData["particle_triangle_s"])
         contacts = in_contacts.FieldData["particle_triangle_idx"]
         points_t = np.ndarray((contacts.shape[0],2,3))
         points_t[:,0,:] = in_particles.Points[contacts[:,1],:]
-        d0 = in_bnd.Points[triangles[contacts[:,0],:]][:,1,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
-        d1 = in_bnd.Points[triangles[contacts[:,0],:]][:,2,:] - in_bnd.Points[triangles[contacts[:,0],:]][:,0,:]
-        N = np.cross(d0,d1)
-        N /= np.linalg.norm(N,axis=1)[:,newaxis]
-        dd = in_bnd.Points[triangles[contacts[:,0],:]][:,0,:] - in_particles.Points[contacts[:,1],:]
-        dist = einsum('ij,ij->i', N, dd)
-        points_t[:,1,:] = in_particles.Points[contacts[:,1],:] + N*dist[:,newaxis]
+        points_t[:,1,:] = np.mean(in_bnd.Points[triangles[contacts[:,0],:]],axis=1)
         points.append(points_t)
 
-    #merge everything
-
-    points = np.vstack(points)
-    vals = np.hstack(vals)
-    vals_t = np.hstack(vals_t)
-    if vals_s is not None :
-        vals_s = np.hstack(vals_s)
     #generate tubes
-    t = points[:,0,:]-points[:,1,:]
-    t /= np.sqrt(t[:,0,None]**2+t[:,1,None]**2+t[:,2,None]**2)
-    ez = np.where(np.abs(t[:,1,None])>np.abs(t[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
-    n1 = np.cross(t,ez)
-    n2 = np.cross(t,n1)
+    points = np.vstack(points)
+    b0 = points[:,0,:]-points[:,1,:]
+    b0 /= np.sqrt(b0[:,0,None]**2+b0[:,1,None]**2+b0[:,2,None]**2)
+    ez = np.where(np.abs(b0[:,1,None])>np.abs(b0[:,2,None]),np.array([[0,0,1]]),np.array([[0,1,0]]))
+    b1 = np.cross(b0,ez)
+    b2 = np.cross(b0,b1)
+    #merge everything
+    vals = np.vstack(vals)
+    if vals.shape[1] == 3:
+        vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1]+vals[:,2]*b0[:,2])
+        vals_t = np.hypot(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1]+vals[:,2]*b1[:,2],
+                          vals[:,0]*b2[:,0]+vals[:,1]*b2[:,1]+vals[:,2]*b2[:,2])
+    else:
+        vals_n = np.abs(vals[:,0]*b0[:,0]+vals[:,1]*b0[:,1])
+        vals_t = np.abs(vals[:,0]*b1[:,0]+vals[:,1]*b1[:,1])
     alphas = np.arange(0,2*np.pi, 2*np.pi/nf)
-    r = rf*vals**(1./2)*1
+    r = rf*vals_n**(1./2)*1
 
     opoints = points[:,None,:,:] \
-        +n1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
-        +n2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
+        +b1[:,None,None,:]*r[:,None,None,None]*np.sin(-alphas)[None,:,None,None] \
+        +b2[:,None,None,:]*r[:,None,None,None]*np.cos(-alphas)[None,:,None,None]
     output.Points = opoints.reshape(-1,3)
-    output.PointData.append(np.repeat(vals,2*nf),"Reaction")
+    output.PointData.append(np.repeat(vals_n,2*nf),"Reaction")
     output.PointData.append(np.repeat(vals_t,2*nf),"Reaction_t")
-    if vals_s is not None :
-        output.PointData.append(np.repeat(vals_s,2*nf),"Reaction_s")
     n = points.shape[0]
-    pattern = np.ndarray([nf,4], np.int)
+    pattern = np.ndarray([nf,4], np.int32)
     for i in range(nf) :
         j = (i+1)%nf
         pattern[i,:] = (i*2,i*2+1,j*2+1,j*2)
@@ -139,4 +119,5 @@ def RequestData():
     cells = np.ndarray([n,nf,5],np.uint32)
     cells[:,:,0] = 4
     cells[:,:,1:] = np.arange(0,n*nf*2,nf*2,np.uint32)[:,None,None]+pattern[None,:,:]
-    output.SetCells(types, locations, cells.reshape([-1]))
\ No newline at end of file
+    output.SetCells(types, locations, cells.reshape([-1]))
+

python/scontact.py

@@ -345,9 +345,7 @@ class ParticleProblem :
         v = np.ndarray((n,self._dim),dtype=np.float64,order="C")
         o = np.ndarray((n,2),dtype=np.uint64,order="C")
         self._lib.particleProblemContact(self._p,c_int(ctype),c_void_p(o.ctypes.data),c_void_p(v.ctypes.data))
-        self._lib.particleProblemContactBasis(self._p,c_int(ctype),c_void_p(basis.ctypes.data))
-        basis = basis.reshape(n,self._dim, self._dim)
-        return o,v,basis
+        return o,v
 
     def computeStressTensor(self, nodes, radius):
         """Computes the stress tensor of the granular material
@@ -416,13 +414,9 @@ class ParticleProblem :
         self._lib.particleProblemContactN.restype = c_size_t
         fdata = []
         for name in ("particle_particle","particle_disk","particle_segment","particle_triangle") :
-            o,v,basis = self.get_contacts(name)
+            o,v = self.get_contacts(name)
             nameb = name.encode()
-            fdata.append((nameb,v[:,[0]]))
-            fdata.append((nameb+b"_t",v[:,[1]]))
-            fdata.append((nameb+b"_dir_n",basis[:,0,:]))
-            if self._dim == 3:
-                fdata.append((nameb+b"_s",v[:,[2]]))
+            fdata.append((nameb, v))
             fdata.append((nameb+b"_idx",o))
         x = np.array([[0.,0.,0.]])
         VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
@@ -496,25 +490,10 @@ class ParticleProblem :
 
         _,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
         ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
-        v = []
-        basis = []
-        v.append(np.vstack([fdata[k] for k in ks]))
-        v.append(np.vstack([fdata[k+"_t"] for k in ks]))
-        basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks]))
-        if self._dim == 2:
-            basis.append(np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[1,0]]*np.array([-1, 1]))
-        else:
-            dirt = np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[2,0,1]]*np.array([-1, 1, -1]) - np.vstack([fdata[k+"_dir_n"] for k in ks])*np.einsum("ij,ij->i",np.vstack([fdata[k+"_dir_n"] for k in ks]),np.vstack([fdata[k+"_dir_n"] for k in ks])[:,[2,0,1]]*np.array([-1, 1, -1]))[:,np.newaxis]
-            nnorm = np.linalg.norm(dirt,axis=1)
-            nnorm[nnorm==0] = 1
-            basis.append(dirt/nnorm[:,np.newaxis])
-            basis.append(np.cross(np.vstack([fdata[k+"_dir_n"] for k in ks]),dirt/nnorm[:,np.newaxis]))
-            v.append(np.vstack([fdata[k+"_s"] for k in ks]))
-        basis = np.column_stack(basis)
-        v = np.column_stack(v)*contact_factor
+        v = np.vstack(list([fdata[k] for k in ks]))*contact_factor
         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])
-        self._lib.particleProblemSetContacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(basis),_np2c(t,np.int32))
+        self._lib.particle_problem_set_contacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(t,np.int32))
 
     def add_boundary_disk(self, x0, r, tag, material="default") :
         """Adds a boundary disk.