contacts-extract.py 4.26 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
# Marblesbag - Copyright (C) <2010-2018>
2
3
# <Universite catholique de Louvain (UCL), Belgium
#  Universite de Montpellier, France>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
# 	
5
# List of the contributors to the development of Marblesbag: see AUTHORS file.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
6
7
8
# Description and complete License: see LICENSE file.
# 	
# This program (Marblesbag) is free software: 
9
# you can redistribute it and/or modify it under the terms of the GNU Lesser General 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
10
11
12
13
14
15
# 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
16
# GNU Lesser General Public License for more details.
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
17
# 
18
19
# 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, 
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
20
21
# see <http://www.gnu.org/licenses/>.

22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
Name = 'MBExtractContact'
Label = 'Marblesbag extract contacts'
Help = 'Generate tubes to represent contacts from marblesbag scontact output'

NumberOfInput = 1

InputDataType = 'vtkMultiBlockDataSet'
OutputDataType = 'vtkUnstructuredGrid'
ExtraXml = ''

Properties = dict(nf = 5, rf = 2.5e-4)


def RequestData():
    import numpy as np

    #input
    #nf = 5
    #rf = 2.5e-4

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

    #split bnd cells
    n_disks, n_segments, n_triangles = in_bnd.FieldData["sizes"]
    if n_disks :
        disks = in_bnd.Cells[1:2*n_disks:2]
    if n_segments :
        segments = in_bnd.Cells[2*n_disks:2*n_disks+3*n_segments].reshape([-1,3])[:,1:]
    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 = []
    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
    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
    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],:]
    points_s[:,1,:] = np.mean(in_bnd.Points[segments[contacts[:,0],:]],axis=1)
    points.append(points_s)

    #particle-triangle contacts
    if "particle_triangle" in in_contacts.FieldData :
        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)

    #merge everything
    points = np.vstack(points)
    vals = np.hstack(vals)

    #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(t[:,1,None]>t[:,2,None],np.array([[0,0,1]]),np.array([[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")
    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]))