# rolling_cylinder.mbd

#-----------------------------------------------------------------------------
# [Data Block]

begin: data;
   problem: initial value;
end: data;

#-----------------------------------------------------------------------------
# [<Problem> Block]

begin: initial value;
   initial time:   0.;
   final time:     2.;
   time step:      1.e-3;
   max iterations: 10;
   tolerance:      1.e-7;
end: initial value;

#-----------------------------------------------------------------------------
# [Control Data Block]

begin: control data;
   output frequency: 5;
   structural nodes: 3;
   abstract nodes:   1;
   rigid bodies:     1;
   joints:           3;
   forces:           3;
   genels:           1;
   gravity;
end: control data;

#-----------------------------------------------------------------------------
# Design Variables
set: real Cylinder_Mass           = 5.;    #[kg]
set: real Cylinder_Radius         = 0.2;   #[m]
set: real Surface_Angle           = pi/6.; #[rad]
set: real Offset_Cylinder_Surface = 0.4;   #[m] Initial Offset of Cylinder Center
                                           #       from Surface (>=Cylinder_Radius)

set: real Kc = 100000.; #[N/m]  Contact Stiffness
set: real Ec = 1.2;     #[-]    Contact Stiffness Exponent
set: real Cc = 10.;     #[Ns/m] Contact Damping
set: real Dc = 0.0001;  #[m]    Contact Damping Activation Depth

set: real Mu = 0.1;   #[-]   Friction Coefficient
set: real Vt = 0.001; #[m/s] Friction Threshold Velocity

#-----------------------------------------------------------------------------
# Intermediate Variables
set: real Cylinder_Inertia = 1./2.*Cylinder_Mass*Cylinder_Radius^2.; #[kg-m^2]

#-----------------------------------------------------------------------------
# Reference Labels
set: integer Ref_Surface = 1;

# Node Labels
set: integer NoSta_Surface      = 1;
set: integer NoDyn_Cylinder     = 2;
set: integer NoDum_Cylinder     = 3;
set: integer NoAbs_Force_Normal = 4;

# Body Labels
set: integer Body_Cylinder = 1;

# Joint Labels
set: integer JoClamp_Surface           = 1;
set: integer JoInp_Surface_Cylinder    = 2;
set: integer JoRevrot_Surface_Cylinder = 3;

# Force Labels
set: integer FoStrin_Normal_Cylinder_Surface   = 1;
set: integer FoStrin_Friction_Cylinder_Surface = 2;
set: integer CoStrin_Friction_Cylinder_Surface = 3;

# Genel Labels
set: integer GeClamp_NoAbs_Force_Normal = 1;

#-----------------------------------------------------------------------------
# Scalar Functions
scalar function: "cubstep",
   cubicspline, do not extrapolate,
      -0.03, 0.00,
      -0.02, 0.00,
      -0.01, 0.00,
       0.00, 0.00,
       1.00, 1.00,
       1.01, 1.00,
       1.02, 1.00,
       1.03, 1.00;

scalar function: "cubsign",
   cubicspline, do not extrapolate,
      -1.03, -1.00,
      -1.02, -1.00,
      -1.01, -1.00,
      -1.00, -1.00,
       1.00,  1.00,
       1.01,  1.00,
       1.02,  1.00,
       1.03,  1.00;
   
#-----------------------------------------------------------------------------
# References
reference: Ref_Surface,
   reference, global, null,                         # absolute position
   reference, global, euler, Surface_Angle, 0., 0., # absolute orientation
   reference, global, null,                         # absolute velocity
   reference, global, null;                         # absolute angular velocity        

#-----------------------------------------------------------------------------
# [Nodes Block]

begin: nodes;

   #-----------------------------------------------------------------------------
   # Nodes
   structural: NoSta_Surface, static,
      reference, Ref_Surface, null, # absolute position
      reference, Ref_Surface, eye,  # absolute orientation
      reference, Ref_Surface, null, # absolute velocity
      reference, Ref_Surface, null; # absolute angular velocity

   structural: NoDyn_Cylinder, dynamic,
      0., 0., Offset_Cylinder_Surface/cos(Surface_Angle), # absolute position
      euler, 0., pi/2., 0.,                               # absolute orientation
      null,                                               # absolute velocity
      null,                                               # absolute angular velocity
      accelerations;
      
   structural: NoDum_Cylinder, dummy, NoDyn_Cylinder,
      relative frame, NoSta_Surface;
      
   abstract: NoAbs_Force_Normal, algebraic;
      
end: nodes;

#-----------------------------------------------------------------------------
# Plugin Variables
set: [node, W, NoDyn_Cylinder, structural, string="omega[3]"];
set: [node, DZ_Dum, NoDum_Cylinder, structural, string="X[3]"];
set: [node, VY_Dum, NoDum_Cylinder, structural, string="XP[2]"];
set: [node, VZ_Dum, NoDum_Cylinder, structural, string="XP[3]"];

#-----------------------------------------------------------------------------
# [Elements Block]

begin: elements;

   #-----------------------------------------------------------------------------
   # Bodies
   body: Body_Cylinder, NoDyn_Cylinder,
      Cylinder_Mass,                  # mass
      null,                           # relative center of mass
      diag, 1., 1., Cylinder_Inertia; # inertia matrix
   
   #-----------------------------------------------------------------------------
   # Joints   
   joint: JoClamp_Surface,
      clamp,
         NoSta_Surface,
            position, node,
            orientation, node;
            
   joint: JoInp_Surface_Cylinder,
      in plane,
         NoSta_Surface,
            null,       # relative plane position
            1., 0., 0., # relative normal direction
         NoDyn_Cylinder;
         
   joint: JoRevrot_Surface_Cylinder,
      revolute rotation,
         NoSta_Surface,
            hinge,
               euler, 0., pi/2., 0., # relative orientation matrix
         NoDyn_Cylinder;
     
   #-----------------------------------------------------------------------------
   # Genels     
   genel: GeClamp_NoAbs_Force_Normal,
      clamp,
         NoAbs_Force_Normal, abstract,
         string, "max(0,Kc*sign(Cylinder_Radius-DZ_Dum)*abs(Cylinder_Radius-DZ_Dum)^Ec\
                       -Cc*VZ_Dum*model::sf::cubstep((Cylinder_Radius-DZ_Dum)/Dc))";
    
   #-----------------------------------------------------------------------------
   # Plugin Variables        
   set: [dof, FN, NoAbs_Force_Normal, abstract, algebraic];
   
   #-----------------------------------------------------------------------------
   # Forces      
   force: FoStrin_Normal_Cylinder_Surface,
      absolute internal,
         NoDyn_Cylinder,
            position, null, # relative arm
         NoSta_Surface,
            position, null, # relative arm
         single, reference, Ref_Surface, 0., 0., 1., string, "FN"; # force value
         
   force: FoStrin_Friction_Cylinder_Surface,
      absolute internal,
         NoDyn_Cylinder,
            position, null, # relative arm
         NoSta_Surface,
            position, null, # relative arm
         single, reference, Ref_Surface, 0., 1., 0.,
            string, "-model::sf::cubsign((VY_Dum+Cylinder_Radius*W)/Vt)*Mu*FN"; # force value
                           
   couple: CoStrin_Friction_Cylinder_Surface,
      absolute internal,
         NoDyn_Cylinder,
         NoSta_Surface,
         single, reference, Ref_Surface, 1., 0., 0.,
            string, "-Cylinder_Radius\
               *model::sf::cubsign((VY_Dum+Cylinder_Radius*W)/Vt)*Mu*FN"; # couple value
                           
   #-----------------------------------------------------------------------------
   # Gravity 
   gravity: 0., 0., -1., const, 9.81;
   
end: elements;