from yade import plot,qt
from yade.gridpfacet import *
import math

O.dt = 0.000001

Rad = 0.1
Factor=1.5
m_viscosity=3
m_roughness=0.02

O.materials.append(
        FrictMat(
            density=2500,
            frictionAngle=math.radians(35.),
            poisson=0.2,
            young=1e7,
            label='Frict'
        ))

S2 = O.bodies.append(sphere(center=[1.5*Rad,0,Factor*2.1*Rad],radius=Rad,material='Frict',fixed=False))
Stest = O.bodies.append(sphere(center=[1.5*Rad,0,3*Factor*2.1*Rad],radius=Rad,material='Frict',fixed=False))
W1 = O.bodies.append(aabbWalls([(4*Rad,-2*Rad,15*Rad),(0,3*Rad,0)],thickness=0.,material='Frict',oversizeFactor=1.0))
#F1=O.bodies.append(facet([[-10,-2*Rad,0],[10,-2*Rad,0],[0,-2*Rad,10]],fixed=True,material='Frict'))


theta_method= 0.55  #1 : Euler Backward, 0.5 : trapezoidal rule, 0 : not used, 0.55 : suggested optimum
resolution=1        #Change normal component resolution method, 0: Iterative exact resolution (theta method, linear contact), 1: Newton-Rafson dimentionless resolution (theta method, linear contact), 2: Newton-Rafson with nonlinear surface deflection (Hertzian-like contact)
NewtonRafsonTol = 1e-6  #Precision of the resolution of Newton-Rafson's method
O.engines=[
    ForceResetter(),
    InsertionSortCollider(
        [
         Bo1_Sphere_Aabb(aabbEnlargeFactor=Factor, label="aabb"),
         Bo1_Box_Aabb(),
#        Bo1_Wall_Aabb(),
#        Bo1_Facet_Aabb(),
        ],
         verletDist=-0.1,
         allowBiggerThanPeriod=False),

    InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom6D(interactionDetectionFactor=Factor,label="Ig2"),
#         Ig2_Wall_Sphere_ScGeom(),
         Ig2_Box_Sphere_ScGeom6D(),
#        Ig2_Facet_Sphere_ScGeom()

        ],

        [
         Ip2_FrictMat_FrictMat_LubricationPhys(eta=m_viscosity,eps=m_roughness),
        ],

        [
         Law2_ScGeom_ImplicitLubricationPhys(activateNormalLubrication=True,
                                             activateTangencialLubrication=True,
                                             activateTwistLubrication=False,
                                             activateRollLubrication=False,
                                             theta=theta_method,
                                             resolution=resolution,
                                             NewtonRafsonTol=NewtonRafsonTol,
                                             warnedOnce=True,
                                             maxSubSteps=20,
                                             debug=False,),

        ]
    ),
    NewtonIntegrator(damping=0.,gravity=[0,-5,-5],label='newton'),
]


qt.View()