en averaging for Law2_ScGeom_ViscElPhys_Basic, bug or incorrect documentation?

The documentation for Law2_ScGeom_ViscElPhys_Basic (https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.Law2_ScGeom_ViscElPhys_Basic) states the following:

  The damping constant is computed at each contact in order to fulfill the normal restitution coefficient en=(en1 en2)/(en1+en2)

However, this does not appear to be what is actually implemented. It looks to me like the code actually does en = (en1+en2)/2

For example, if I test two materials, one with en=0.9 and the other with en=0.1, the formula stated in the documentation would result in en = 0.9 * 0.1 / (0.9 + 0.1) = 0.09. Example code with dropping a ball onto a wall is below. The expected rebound height of the ball is square root of coefficient of restitution times initial height. i.e. expect a rebound height of 1% of the initial height. The actual rebound height that I get is 0.25, which implies the coefficient of restitution is 0.5. That is consistent with an arithmetic average (0.9+0.1)/2.

So, my question: should it be considered as the documentation is the intended behavior and this is an issue with the code? or the code is the intended behavior and this is an issue with the documentation? or neither one is right, and it should be something else?

The behavior in the documentation does not make physical sense to me. imagine two materials both with en=1. You would end up with (1*1)/(1+1) = 0.5. That doesn't make sense. en=1 implies elastic behavior, but that's not what would happen. So to me, I would suggest to keep current arithmetic averaging and change the documentation to match.

In reality, of course, coefficient of restitution is not really defined for an individual material, but only really for pair of materials (and even then only empirically determined). I looked around to see if I could find any articles of someone suggesting one averaging method or the other, but I could not find any.

###################################################################

import matplotlib.pyplot as pyplot
from yade import qt, plot
qt.View() #open the controlling and visualization interfaces

box_x = 0.05
box_y = 0.05
box_z = 0.05

particle_dia = 0.005

young2 = 200e9
rho = 8230

mn = Vector3(0, box_y, 0)
mx = Vector3(box_x,2*box_y, box_z) #corners of the initial packing
thick = 2*particle_dia # the thickness of the walls

bigmx = (mx[0], 3 * mx[1], mx[2])

restitution_ball = 0.9
restitution_wall = 0.1

O.materials.append(ViscElMat(young=young2,poisson=0.7,frictionAngle=0, en = restitution_ball , et=1, density=rho , label='ve_ball'))
O.materials.append(ViscElMat(young=young2,poisson=0.7,frictionAngle=0, en = restitution_wall , et=1, density=rho , label='ve_wall'))

ball = sphere( (mx[0]/2, 2*mx[1] , mx[2]/2 ) , particle_dia/2, material='ve_ball' )
ballIds = O.bodies.append(ball)

walls=utils.aabbWalls([mn,bigmx],thickness=thick,oversizeFactor=1.5,material='ve_wall')
wallIds=O.bodies.append(walls)

#turn on gravity and let it settle
O.engines=[
 ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()] ),
        InteractionLoop(
                [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
                [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
                [Law2_ScGeom_ViscElPhys_Basic()] ),
 NewtonIntegrator(gravity=(0,-9.81,0),damping=0.0),
        PyRunner(command='addPlotData()', iterPeriod=1000),
]
O.dt=.05*PWaveTimeStep()

def addPlotData():
        plot.addData(time=O.time , pos = (O.bodies[ballIds].state.pos[1]-box_y)/0.15 )

plot.plots={'time':('pos')}
plot.plot()

#O.run(-1, True)