# Confused by identificationSpheresOnJoint.py

Dear all,
Because my research is about fractured rock mass, I decided to learn relevant examples[1].When I read identificationspheresonjoint.py, the following code I can't fully understand:

1.###if isinstance(O.bodies[i.id1].shape,Facet) and isinstance(O.bodies[i.id2].shape,Sphere)###
In the above code,if we judge that Id1 corresponds to sphere and Id2 corresponds to facet,is it the same result?

2.There are relevant judgment statements about the number of joints and jointNormal1,2,3 in the code.My understanding is that the number of joints to which a particle belongs is not limited. In that way,why is jointnormal only 1,2,3.
(I noticed the words 'plane' and 'surface'. Is that the difference here?)

3.###if O.bodies[i.id2].state.joint==1 and ((O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05)###
I don't understand the meaning of #(O.bodies[i.id2].state.jointNormal1.cross(nRef)).norm()>0.05#.

4.###if ((O.bodies[j.id2].state.jointNormal1.cross(nRef)).norm()<0.05) :
jointNormalRef=O.bodies[j.id2].state.jointNormal1
elif ((O.bodies[j.id2].state.jointNormal2.cross(nRef)).norm()<0.05) :
jointNormalRef=O.bodies[j.id2].state.jointNormal2
elif ((O.bodies[j.id2].state.jointNormal3.cross(nRef)).norm()<0.05) :
jointNormalRef=O.bodies[j.id2].state.jointNormal3
else : continue ###
In the above code, it seems that there must be a judgment statement < 0.05, which is my confusion.

Best regards.
ps:I'm not sure whether my question follows regulations . If there is anything inappropriate, please point it out and forgive me.

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Ziyu Wang
Solved:
Last query:
 Revision history for this message Luc Scholtès (luc) said on 2021-08-30: #1

Hi,

I'll try to answer even though I wrote this code almost 10 years ago:

1. I think there was a rule (in YADE collider) which defined facet-sphere interactions such as body1 is always the facet and body2 always the sphere. See for instance: https://answers.launchpad.net/yade/+question/690010. You may want to doublecheck if this rule is still effective.

2a. I arbitrarily decided that each particle would, at the maximum, belong to 3 discontinuities at the same time. I believe this is quite reasonable when you consider a fractured rock mass but you can, if you want, change that number (it might complexify the numerical treatment though).

2b. I could not find the word "surface" in the script so I cannot answer your concern about it.

3. The cross product of parallel vectors is a zero vector whose norm is (strictly) equal to 0. Testing the norm of the cross product thus enables to discriminate between different/non parallel joint surfaces: if the norm of the cross product of their respective normals is not equal to 0, then the joint surfaces are different/non parallel. The 0.05 I introduced was a way to take into account some sort of error/margin in the joint surfaces orientation and to avoid the definition of 2 different joint planes for sub-parallel surfaces. This was, again, an arbitrary choice from myself si you can, if you want, modify that part of the code.

I hope it helps.

Luc

 Revision history for this message Ziyu Wang (ziyuwang1) said on 2021-08-31: #3

Thanks Luc,I'm glad you can answer my question.But I'm still a little confused.

>2b. I could not find the word "surface" in the script so I cannot answer your concern about it.

I found the 'plane' and surface in the help document[1] as follows:
joint:Indicates the number of joint surfaces to which the particle belongs;
jointNormal1:Specifies the normal direction to the joint plane 1.
Is there any difference between the two words here?

And I found a new problem about it: ###O.bodies[i.id2].state.jointNormal1=nRef###
nRef represents the normal direction of the facet.Then the meaning of the above code is to set jointnormal1 as the normal direction of the facet.Can I understand joint plane 1 as facet face?If so, which plane is joint plane 2 and 3?I'm very confused.

>The cross product of parallel vectors is a zero vector whose norm is (strictly) equal to 0...
After your answer, I suddenly enlightened.But must one of jointNormal1,2,3 be parallel to nRef?

###if O.bodies.state.joint==1....
O.bodies.state.joint==2###
By the way,I suddenly found that I didn't quite understand why you put joint + 1 in the judgment statement,what is the purpose?

It seems that I have too many questions...Thank you again for your help!

 Revision history for this message Luc Scholtès (luc) said on 2021-08-31: #4

I'll try to answer but, unfortunately, you'll have to go through the code by yourself to understand everything in details:

- in the statements you mentioned, surface and plane have the same meaning: they both represent the discontinuity plane/surface.

- you may have several fracture planes/surfaces in your model. The code is written so that you can have maximum 3 planes associated to each particle. For example, in a 2D model, if you have a fracture plane dipping at 20° and another one dipping at 60°, they might intersect at some location in your model. The particles located around the intersection point must then have in memory the fact that they belong to 2 planes with different inclinations, hence the jointNormal1 and the jointNormal2. Then, the code will be able to know if particles located on each side of the planes belong to plane1 or to plane2 and adjust the geometry (and behavior) according to their respective orientations (normals). For instance, if there is an interaction between a particle belonging to plane1 and a particle belonging to plane1 and plane 2, the interaction geometry will be defined according to the geometry of plane1. If they both belong to plan1e and to plane2, I think the interaction geometry is not affected so that the overall behavior is only affected by the surrounding well defined interactions (according to plane1 or to plane2).

- nRef corresponds to the normal to the facet that is tested within the interaction loop and which represents the actual discontinuity plane (temporary attribute). JointNormali is the particle attribute (permanent attribute): the code defines jointNormali according to nRef when and then move on the next interaction.

Be careful: the code is made so that there is 2 loops:
- the first one which tests interactions between facet and spheres to identify which particles belong to the fracture plane.
- the second one which tests interactions between the particles belonging to the fracture plane (identified in the first loop) and particles in contact with these particles.

It is good that you want to understand but may I ask what is your goal exactly? The code has been verified already and seems to work properly. Do you plan to make some modifications?

Luc

 Revision history for this message Ziyu Wang (ziyuwang1) said on 2021-08-31: #5

Thanks a lot,Luc
After your patient explanation, I carefully read the example again and felt that I understood a lot.
My current ability is not enough to modify the code independently.My goal is just to understand the code thoroughly