13
13
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=interactionRadius,label='is2aabb'),Bo1_Facet_Aabb()]),
15
15
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=interactionRadius,label='ss2d3dg'),Ig2_Facet_Sphere_ScGeom()],
16
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,alpha=1,tensileStrength=1e6,cohesion=1e6,jointNormalStiffness=1,jointShearStiffness=1,jointTensileStrength=1e6,jointCohesion=1e6,jointFrictionAngle=1,label='interactionPhys')],
16
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1,label='interactionPhys')],
17
17
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=True,label='interactionLaw')]
19
19
NewtonIntegrator(damping=1)
47
47
nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
48
48
normalFacetSphere=i.geom.normal # geom.normal is oriented from id1 to id2 -> normalFacetSphere from facet (i.id1) to sphere (i.id2)
50
if O.bodies[i.id2].mat.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
51
O.bodies[i.id2].mat.onJoint=True
52
O.bodies[i.id2].mat.joint=1
50
if O.bodies[i.id2].state.onJoint==False : ## particles has not yet been identified as belonging to a joint plane
51
O.bodies[i.id2].state.onJoint=True
52
O.bodies[i.id2].state.joint=1
53
53
O.bodies[i.id2].shape.color=jointcolor1
54
54
if nRef.dot(normalFacetSphere)>=0 :
55
O.bodies[i.id2].mat.jointNormal1=nRef
55
O.bodies[i.id2].state.jointNormal1=nRef
56
56
elif nRef.dot(normalFacetSphere)<0 :
57
O.bodies[i.id2].mat.jointNormal1=-nRef
58
elif O.bodies[i.id2].mat.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
59
if O.bodies[i.id2].mat.joint==1 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
60
O.bodies[i.id2].mat.joint=2
57
O.bodies[i.id2].state.jointNormal1=-nRef
58
elif O.bodies[i.id2].state.onJoint==True : ## particles has already been identified as belonging to, at least, 1 facet
59
if O.bodies[i.id2].state.joint==1 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) : ## particles has already been identified as belonging to only 1 facet
60
O.bodies[i.id2].state.joint=2
61
61
O.bodies[i.id2].shape.color=jointcolor2
62
62
if nRef.dot(normalFacetSphere)>=0 :
63
O.bodies[i.id2].mat.jointNormal2=nRef
63
O.bodies[i.id2].state.jointNormal2=nRef
64
64
elif nRef.dot(normalFacetSphere)<0 :
65
O.bodies[i.id2].mat.jointNormal2=-nRef
66
elif O.bodies[i.id2].mat.joint==2 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
67
O.bodies[i.id2].mat.joint=3
65
O.bodies[i.id2].state.jointNormal2=-nRef
66
elif O.bodies[i.id2].state.joint==2 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05): ## particles has already been identified as belonging to more than 1 facet
67
O.bodies[i.id2].state.joint=3
68
68
O.bodies[i.id2].shape.color=jointcolor3
69
69
if nRef.dot(normalFacetSphere)>=0 :
70
O.bodies[i.id2].mat.jointNormal3=nRef
70
O.bodies[i.id2].state.jointNormal3=nRef
71
71
elif nRef.dot(normalFacetSphere)<0 :
72
O.bodies[i.id2].mat.jointNormal3=-nRef
73
elif O.bodies[i.id2].mat.joint==3 and ((O.bodies[i.id2].mat.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].mat.jointNormal3.cross(nRef)).norm()>0.05):
74
O.bodies[i.id2].mat.joint=4
72
O.bodies[i.id2].state.jointNormal3=-nRef
73
elif O.bodies[i.id2].state.joint==3 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal2.cross(nRef)).norm()>0.05) and ((O.bodies[i.id2].state.jointNormal3.cross(nRef)).norm()>0.05):
74
O.bodies[i.id2].state.joint=4
75
75
O.bodies[i.id2].shape.color=jointcolor5
77
77
#### second step -> find spheres interacting with spheres on facet (could be executed in the same timestep as step 1?)
82
82
vertices=O.bodies[j.id1].shape.vertices
83
83
normalRef=vertices[0].cross(vertices[1]) # defines the normal to the facet normalRef
84
84
nRef=normalRef/(normalRef.norm()) ## normalizes normalRef
85
if ((O.bodies[j.id2].mat.jointNormal1.cross(nRef)).norm()<0.05) :
86
jointNormalRef=O.bodies[j.id2].mat.jointNormal1
87
elif ((O.bodies[j.id2].mat.jointNormal2.cross(nRef)).norm()<0.05) :
88
jointNormalRef=O.bodies[j.id2].mat.jointNormal2
89
elif ((O.bodies[j.id2].mat.jointNormal3.cross(nRef)).norm()<0.05) :
90
jointNormalRef=O.bodies[j.id2].mat.jointNormal3
85
if ((O.bodies[j.id2].state.jointNormal1.cross(nRef)).norm()<0.05) :
86
jointNormalRef=O.bodies[j.id2].state.jointNormal1
87
elif ((O.bodies[j.id2].state.jointNormal2.cross(nRef)).norm()<0.05) :
88
jointNormalRef=O.bodies[j.id2].state.jointNormal2
89
elif ((O.bodies[j.id2].state.jointNormal3.cross(nRef)).norm()<0.05) :
90
jointNormalRef=O.bodies[j.id2].state.jointNormal3
92
92
facetCenter=O.bodies[j.id1].state.pos
93
93
#### seek for each sphere interacting with the identified sphere j.id2
102
102
facetSphereDir=(O.bodies[othSph].state.pos-facetCenter)
103
if O.bodies[othSph].mat.onJoint==True :
104
if O.bodies[othSph].mat.joint==3 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal3.cross(jointNormalRef)).norm()>0.05):
105
O.bodies[othSph].mat.joint=4
103
if O.bodies[othSph].state.onJoint==True :
104
if O.bodies[othSph].state.joint==3 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal3.cross(jointNormalRef)).norm()>0.05):
105
O.bodies[othSph].state.joint=4
106
106
O.bodies[othSph].shape.color=jointcolor5
107
elif O.bodies[othSph].mat.joint==2 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].mat.jointNormal2.cross(jointNormalRef)).norm()>0.05):
108
O.bodies[othSph].mat.joint=3
109
if facetSphereDir.dot(jointNormalRef)>=0:
110
O.bodies[othSph].mat.jointNormal3=jointNormalRef
111
elif facetSphereDir.dot(jointNormalRef)<0:
112
O.bodies[othSph].mat.jointNormal3=-jointNormalRef
113
elif O.bodies[othSph].mat.joint==1 and ((O.bodies[othSph].mat.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
114
O.bodies[othSph].mat.joint=2
115
if facetSphereDir.dot(jointNormalRef)>=0:
116
O.bodies[othSph].mat.jointNormal2=jointNormalRef
117
elif facetSphereDir.dot(jointNormalRef)<0:
118
O.bodies[othSph].mat.jointNormal2=-jointNormalRef
119
elif O.bodies[othSph].mat.onJoint==False :
120
O.bodies[othSph].mat.onJoint=True
121
O.bodies[othSph].mat.joint=1
107
elif O.bodies[othSph].state.joint==2 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) and ((O.bodies[othSph].state.jointNormal2.cross(jointNormalRef)).norm()>0.05):
108
O.bodies[othSph].state.joint=3
109
if facetSphereDir.dot(jointNormalRef)>=0:
110
O.bodies[othSph].state.jointNormal3=jointNormalRef
111
elif facetSphereDir.dot(jointNormalRef)<0:
112
O.bodies[othSph].state.jointNormal3=-jointNormalRef
113
elif O.bodies[othSph].state.joint==1 and ((O.bodies[othSph].state.jointNormal1.cross(jointNormalRef)).norm()>0.05) :
114
O.bodies[othSph].state.joint=2
115
if facetSphereDir.dot(jointNormalRef)>=0:
116
O.bodies[othSph].state.jointNormal2=jointNormalRef
117
elif facetSphereDir.dot(jointNormalRef)<0:
118
O.bodies[othSph].state.jointNormal2=-jointNormalRef
119
elif O.bodies[othSph].state.onJoint==False :
120
O.bodies[othSph].state.onJoint=True
121
O.bodies[othSph].state.joint=1
122
122
O.bodies[othSph].shape.color=jointcolor4
123
123
if facetSphereDir.dot(jointNormalRef)>=0:
124
O.bodies[othSph].mat.jointNormal1=jointNormalRef
124
O.bodies[othSph].state.jointNormal1=jointNormalRef
125
125
elif facetSphereDir.dot(jointNormalRef)<0:
126
O.bodies[othSph].mat.jointNormal1=-jointNormalRef
126
O.bodies[othSph].state.jointNormal1=-jointNormalRef
128
128
#### for visualization:
132
132
#for o in O.bodies:
133
#if o.mat.onJoint==True : # or o.shape.name=='Facet':
133
#if o.state.onJoint==True : # or o.shape.name=='Facet':
134
134
##if o.shape.name=='Facet':
135
135
##o.shape.wire=True
136
136
##o.state.pos+=(0,50,0)
138
#if o.mat.jointNormal1.dot(hor)>0 :
138
#if o.state.jointNormal1.dot(hor)>0 :
139
139
##o.state.pos+=(0,50,0)
140
140
#o.shape.color=jointcolor1
141
#elif o.mat.jointNormal1.dot(hor)<0 :
141
#elif o.state.jointNormal1.dot(hor)<0 :
142
142
##o.state.pos+=(0,55,0)
143
143
#o.shape.color=jointcolor2
144
144
#if o.mat.type>2 :
146
146
#o.shape.color=jointcolor5
147
##print o.mat.jointNormal.dot(vert)
147
##print o.state.jointNormal.dot(vert)
150
150
##### to delete facets (should be OK to start the simulation after that!