Welcome to ask.woodem.org. You may post when you login through your GitHub account.

How to detect the Contact Geometry and Contact Force?

0 votes
asked Oct 13, 2016 by MarcoTar (770 points)

Hi Vaclav.

As your request, I poste the question here.

So I repeat what I asked in the previous questions in order to let everyone have a better view.

Here I provide this Rough Picture.

My question is this:

1. How can I determine the particle whether it contacts the membrane?

2. After determining the number of particles contacting the membrane, how can I make a sum of the conctact force which comes from the contact between particle and membrane?



Here I want to know if the total mass of the spheres = volume*num*density ?

Is there any method I can set the total mass of the particles?


Best Regards

Feixiang Xuan

commented Oct 15, 2016 by MarcoTar (770 points)

Add 1 more question:

I can access the node pos. However, I think the Number of the node is random instead of fixed.

For example, the four corner has fixed the translational D.o.F, which means the Pos[2] should be zero. If I type S.dem.nodes[0].pos[2], the displacement is not zero meaning this node is not what we consider at the corner. Besides, I checked the 'Inspector' for the dem.nodes and I found every time I run the simulation, the nodes position is re-arranged again.

So my question is that is there any command which can make the coordinate of the node fixed? Or could you tell me how does the code arrange the coordinate of the node?

enter image description here

Best Regards

Feixiang Xuan

commented Oct 17, 2016 by MarcoTar (770 points)

add 1 more:

could you tell me the right way to access the node force? I want the reaction force data.

I tried S.dem.node.force[2] it doesn't work.


Feixiang Xuan

1 Answer

0 votes
answered Oct 17, 2016 by eudoxos (49,070 points)

Hi, please post unrelated questions separately, otherwise it gets quite messy. Using your question numbers:

1. to access all contacts with the memberane, you can either loop through all membrane elements and then over all its contacts, or over all contacts and filter out only those with membrane:

# option one
for c in S.dem.con:
   # skip contacts where neither particles is a Membrane
   if not isinstance(c.pA.shape,woo.fem.Membrane) and not isinstance(c.pB.shape,woo.dem.Membrane): continue
   # do something, such as sum forces etc.

# option two
for p in S.dem.par:
   if not isinstance(p.shape,woo.fem.Membrane): continue
   for c in p.con:
      # do something with c, which is a contact between Membrane and something

2. I'd guess every time you run the script, the numbering will be the same -- deterministic but arbitrary (maybe not if python sequences come into play). But node numbering is an implementation detail and is not guaranteed to be anything. You should loop over all nodes and find the corner nodes by comparing coordinates -- something like this (make sure you do it just once and store the result e.g. in S.lab, as it traverses all nodes every time called):

from minieigen import *
def findNodesAround(S.dem,coord,eps=1e-3):
   # construct box around the coordinate
   # return all nodes within the box, the caller should get the item (s)he wants
   return [n for n in S.dem.nodes if n.pos in box]
S.lab.cornerNodes=[findNodesAround(S.dem,coord)[0] for coord in (0,0,0),(0,1,0),(1,0,0),(0,0,1)] # in this example, corners of unit square in xy-plane

3. To set the total mass, find volume and compute density from required mass and volume (to compute the volume, call SpherePack.solidVolume, obviously after calling makeCloud and before calling toSimulation). If you use more sophisticated generators, you can set the mass you require directly (see tutorial on inlets).

4. Nodal force: Node.dem.force (DemData.force, read the docs, really) so it will be S.dem.nodes[i].dem.force if you start with Scene, i being the node number.

Cheers, v.

commented Oct 19, 2016 by MarcoTar (770 points)

Hi Vaclav
Thanks for your kind answer.
But here for the node numbering (2nd question), I tried your method here it gives me the exact node <node@address coord> instead of the node id (number). For instance S.dem.nodes[i].pos[2], I need the i.
Maybe I didn't either understand your method or you misunderstand my meaning. Anyhow, I solve this question by myself through the following code:

    # Find the interesting nodes id through coordinates
for i in range(len(S.dem.nodes)):
    if (S.dem.nodes[i].pos[0]==0) and (S.dem.nodes[i].pos[1]==0): 
    elif (S.dem.nodes[i].pos[0]==0) and (S.dem.nodes[i].pos[1]==2.8): 
    elif (S.dem.nodes[i].pos[0]==2.8) and (S.dem.nodes[i].pos[1]==0): 
    elif (S.dem.nodes[i].pos[0]==2.8) and (S.dem.nodes[i].pos[1]==2.8): 
    else: i=i+1`

Now I'm going to do the contact loop, before doing it, could you give me some details explanation about the 'PA 'and 'PB'?

Thanks very much
Best Regards

Feixiang Xuan


commented Oct 20, 2016 by eudoxos (49,070 points)

Hi, that's correct, you are getting the Node instance, which is what you want anyway. You can store it and re-use it later (if you want to keep it across saves, put it in S.lab) as a variable. Nodes don't have IDs, they are numbered, but the number can change anytime (if nodes are deleted, for example), it is just an implementation detail you should not be relying upon; why do you need it anyway? (Particles, on the other hand, do have guaranteed-stable IDs).

For pA and pB see documentation, it is the way to access particles in the contact (they are always 2).

commented Oct 20, 2016 by MarcoTar (770 points)

I need these nodes numbering because I want to get the results of the nodal displacement and nodal reaction forces varying with time instead of a final answer.

Frankly speaking, I need to plot these results in a diagram. So there are several node position that I'm interested in.


commented Oct 20, 2016 by eudoxos (49,070 points)

You don't need to know number of nodes, they have no meaning, it is just and internal detail. You need to get those corner nodes, which I what showed you before. So when you get the nodes, just ask them for displacement. Provided you have S.lab.cornerNodes as I showed you above, you can do something like:

def addPlotData(S):
   # do whatever you need
   cn=S.lab.cornerNodes # we are lazy
  # plot z-positions and z-forces of corner nodes