Surface area and overlap

Hi,

you can iterate over interactions and find read penetration depth. Based on this, you can compute overlapping volume manually. Please note it is not a ready-for-use code, but just schematic pseudocode:

###
totalOverlapVolume = 0

for ii in O.interactions: #iterate over all the interactions
    b1 = O.bodies[ii.id1]
    b2 = O.bodies[ii.id2]
    penetrationDepth = ii.geom.penetrationDepth # Should work for Mindlin law. You can inspect interactions manually to find out whether .penetrationDepth is a right property to read.
    overlapVolume = yourFunctionComputingOverlap(b1,b2,penetrationDepth)
    totalOverlapVolume += overlapVolume

###
yourFunctionComputingOverlap() should check the shape of the particles b1 and b2 (whether it is a sphere, wall, wall, etc), extract some shape properties (e.g. radii) and compute volume based on geometrical relationships.

Best wishes,
Karol

Hello,
Alternative approach, based on each particle. See also [2]:
### (not tested)
def computeSurfaceOfOneSphere(body):
    radius = b.shape.radius
    surface = 4*pi*pow(radius,2)
    interactions = O.intrs() # [1], interactions of the particle
    contactPoints = [i.geom.contactPoint for i in interactions]
    center = b.state.pos
    capHeights = [radius - (cp-center).norm() for cp in contactPoints] # [2]
    caps = [2*pi*radius*h for h in capHeights] # [2]
    surfaceActual = surface - sum(caps)
    return surfaceActual

surfaces = [computeSurfaceOfOneSphere(b) for b in O.bodies]
surfaceTotal = sum(surfaces)
###

Cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Body.intrs
[2] https://en.wikipedia.org/wiki/Spherical_cap

Hi,

Thank you for your help, I am getting error message saying:

"unsupported operand type(s) for -: 'Vector3' and 'builtin_function_or_method"

I probably need do some modification on your code. I was wondering how is the contactPoints defined, where on the on sphere is it measured.

Best regards,
Mithushan

The contactPoint is probably not a good idea for different radii.. Better use [3].

A complete working example:
###
O.bodies.append([
    sphere((1,2,3),4,fixed=True),
    sphere((4,5,6),3,fixed=True),
])

def computeCapSurface(body,interaction):
    radius = body.shape.radius
    surface = 4*pi*pow(radius,2)
    center = body.state.pos
    if interaction.id2 == body.id: id2 = interaction.id1
    else: id2 = interaction.id2
    body2 = O.bodies[id2]
    radius2 = body2.shape.radius
    center2 = body2.state.pos
    d = (center2-center).norm()
    d1 = (pow(d,2)-pow(radius2,2)+pow(radius,2)) / (2*d) # [3]
    print(radius,d,d1)
    capHeight = radius - d1 # [2]
    capSurface = 2*pi*radius*capHeight # [2]
    return capSurface

def computeSurfaceOfOneSphere(body):
    interactions = body.intrs() # [1]
    radius = body.shape.radius
    surface = 4*pi*pow(radius,2)
    capSurfaces = [computeCapSurface(body,interaction) for interaction in interactions]
    surfaceActual = surface - sum(capSurfaces)
    return surfaceActual

O.step()

surfaces = [computeSurfaceOfOneSphere(b) for b in O.bodies]
surfaceTotal = sum(surfaces)
###

Double check the math and equations, I give no warranty that they are correct :-)

Cheers
Jan

[3] https://mathworld.wolfram.com/Circle-CircleIntersection.html

> How can I measure the average overlap ratio between the particles?

What is definition "overlap ratio" and "the average overlap ratio"?