# half_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-4;
   max iterations: 10;
   tolerance:      1.e-6;
end: initial value;

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

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

#-----------------------------------------------------------------------------
# Design Variables
set: real m     = 5.;    #[kg]  Mass of the half-cylinder
set: real R     = 0.4;   #[m]   Radius of the half-cylinder
set: real theta = pi/6.; #[rad] Surface angle

set: real Kc = 100000.; #[N/m]  Contact stiffness
set: real Ec = 1.2;     #[-]    Contact stiffness exponent
set: real Cc = 100.;    #[Ns/m] Contact damping
set: real Dc = 0.0001;  #[m]    Contact damping activation depth

set: real mu = 0.3;   #[-]   Friction coefficient
set: real Vt = 0.005; #[m/s] Friction threshold velocity

#-----------------------------------------------------------------------------
# Intermediate Variables
set: real Ic = 1./2.*m*R^2.; #[kg-m^2] Moment of inertia about the cylinder axis
set: real d  = 4./3.*R/pi;   #[m]      Distance from the cylinder axis to the center of mass
set: real I  = Ic-m*d^2;     #[kg-m~2] Moment of inertia about the axis through the center of mass

#-----------------------------------------------------------------------------
# Reference Labels
set: integer Ref_Surface        = 1;
set: integer Ref_CylinderCenter = 2;
set: integer Ref_CylinderEdgeR  = 3;
set: integer Ref_CylinderEdgeL  = 4;

# Structural Node Labels
set: integer NoSta_Surface                     = 1;
set: integer NoDyn_Cylinder                    = 2;
set: integer NoDum_Cylinder_RelTo_Surface      = 3;
set: integer NoDum_CylinderEdgeR               = 4;
set: integer NoDum_CylinderEdgeL               = 5;
set: integer NoDum_CylinderEdgeR_RelTo_Surface = 6;
set: integer NoDum_CylinderEdgeL_RelTo_Surface = 7;

# Abstract Node Labels
set: integer NoAbs_FN_Cylinder       = 1;
set: integer NoAbs_FN_CylinderEdgeR  = 2;
set: integer NoAbs_FN_CylinderEdgeL  = 3;

# 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;
set: integer FoStrin_Normal_CylinderEdgeR_Surface   = 4;
set: integer FoStrin_Normal_CylinderEdgeL_Surface   = 5;
set: integer FoStrin_Friction_CylinderEdgeR_Surface = 6;
set: integer FoStrin_Friction_CylinderEdgeL_Surface = 7;

# Genel Labels
set: integer GeClamp_FN_Cylinder      = 1;
set: integer GeClamp_FN_CylinderEdgeR = 2;
set: integer GeClamp_FN_CylinderEdgeL = 3;

#-----------------------------------------------------------------------------
# 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,
   null,                 # absolute position
   euler, theta, 0., 0., # absolute orientation
   null,                 # absolute velocity
   null;                 # absolute angular velocity

reference: Ref_CylinderCenter,
   reference, Ref_Surface, 0., 0., R, # absolute position
   reference, global, eye,            # absolute orientation
   reference, global, null,           # absolute velocity
   reference, global, null;           # absolute angular velocity
   
reference: Ref_CylinderEdgeR,
   reference, Ref_CylinderCenter, 0., -R, 0., # absolute position
   reference, Ref_CylinderCenter, eye,        # absolute orientation
   reference, Ref_CylinderCenter, null,       # absolute velocity
   reference, Ref_CylinderCenter, null;       # absolute angular velocity
   
reference: Ref_CylinderEdgeL,
   reference, Ref_CylinderCenter, 0., R, 0., # absolute position
   reference, Ref_CylinderCenter, eye,       # absolute orientation
   reference, Ref_CylinderCenter, null,      # absolute velocity
   reference, Ref_CylinderCenter, 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,
      reference, Ref_CylinderCenter, null, # absolute position
      reference, Ref_CylinderCenter, eye,  # absolute orientation
      reference, Ref_CylinderCenter, null, # absolute velocity
      reference, Ref_CylinderCenter, null, # absolute angular velocity
      accelerations, yes;
      
   structural: NoDum_Cylinder_RelTo_Surface, dummy, NoDyn_Cylinder,
      relative frame, NoSta_Surface;

   structural: NoDum_CylinderEdgeR, dummy, NoDyn_Cylinder,
      offset,
         reference, Ref_CylinderEdgeR, null, # relative offset
         reference, Ref_CylinderEdgeR, eye;  # relative orientation
         
   structural: NoDum_CylinderEdgeL, dummy, NoDyn_Cylinder,
      offset,
         reference, Ref_CylinderEdgeL, null, # relative offset
         reference, Ref_CylinderEdgeL, eye;  # relative orientation
         
   structural: NoDum_CylinderEdgeR_RelTo_Surface, dummy, NoDum_CylinderEdgeR,
      relative frame, NoSta_Surface;
      
   structural: NoDum_CylinderEdgeL_RelTo_Surface, dummy, NoDum_CylinderEdgeL,
      relative frame, NoSta_Surface;
      
   abstract: NoAbs_FN_Cylinder, algebraic;
   
   abstract: NoAbs_FN_CylinderEdgeR, algebraic;
   
   abstract: NoAbs_FN_CylinderEdgeL, algebraic;
      
end: nodes;

#-----------------------------------------------------------------------------
# Plugin Variables
set: [node, DZ, NoDum_Cylinder_RelTo_Surface, structural, string="X[3]"];
set: [node, VY, NoDum_Cylinder_RelTo_Surface, structural, string="XP[2]"];
set: [node, VZ, NoDum_Cylinder_RelTo_Surface, structural, string="XP[3]"];
set: [node, WX, NoDum_Cylinder_RelTo_Surface, structural, string="omega[1]"];
set: [node, RX, NoDum_Cylinder_RelTo_Surface, structural, string="E[1]"];

set: [node, DZ_R, NoDum_CylinderEdgeR_RelTo_Surface, structural, string="X[3]"];
set: [node, VY_R, NoDum_CylinderEdgeR_RelTo_Surface, structural, string="XP[2]"];
set: [node, VZ_R, NoDum_CylinderEdgeR_RelTo_Surface, structural, string="XP[3]"];

set: [node, DZ_L, NoDum_CylinderEdgeL_RelTo_Surface, structural, string="X[3]"];
set: [node, VY_L, NoDum_CylinderEdgeL_RelTo_Surface, structural, string="XP[2]"];
set: [node, VZ_L, NoDum_CylinderEdgeL_RelTo_Surface, structural, string="XP[3]"];

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

begin: elements;

   #-----------------------------------------------------------------------------
   # Bodies
   body: Body_Cylinder, NoDyn_Cylinder,
      m,               # mass
      0., 0., -d,      # relative center of mass
      diag, I, 1., 1.; # 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
         NoDyn_Cylinder,
            hinge,
               euler, 0., pi/2., 0.; # relative orientation
     
   #-----------------------------------------------------------------------------
   # Genels
   genel: GeClamp_FN_Cylinder,
      clamp,
         NoAbs_FN_Cylinder, abstract,
         string, "max(0,Kc*sign(R-DZ)*abs(R-DZ)^Ec-Cc*VZ*model::sf::cubstep((R-DZ)/Dc))\
               *step(RX+pi/2.)*step(-RX+pi/2.)";
               
   genel: GeClamp_FN_CylinderEdgeR,
      clamp,
         NoAbs_FN_CylinderEdgeR, abstract,
         string, "max(0,Kc*sign(0-DZ_R)*abs(0-DZ_R)^Ec-Cc*VZ_R*model::sf::cubstep((0-DZ_R)/Dc))\
               *(1-step(RX+pi/2.)*step(-RX+pi/2.))";
               
   genel: GeClamp_FN_CylinderEdgeL,
      clamp,
         NoAbs_FN_CylinderEdgeL, abstract,
         string, "max(0,Kc*sign(0-DZ_L)*abs(0-DZ_L)^Ec-Cc*VZ_L*model::sf::cubstep((0-DZ_L)/Dc))\
               *(1-step(RX+pi/2.)*step(-RX+pi/2.))";
    
   #-----------------------------------------------------------------------------
   # Plugin Variables        
   set: [dof, FN, NoAbs_FN_Cylinder, abstract, algebraic];
   set: [dof, FN_R, NoAbs_FN_CylinderEdgeR, abstract, algebraic];
   set: [dof, FN_L, NoAbs_FN_CylinderEdgeL, 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+R*WX)/Vt)*mu*FN"; # force value
                           
   couple: CoStrin_Friction_Cylinder_Surface,
      absolute internal,
         NoDyn_Cylinder,
         NoSta_Surface,
         single, reference, Ref_Surface, 1., 0., 0.,
            string, "-R*model::sf::cubsign((VY+R*WX)/Vt)*mu*FN"; # couple value
            
   force: FoStrin_Normal_CylinderEdgeR_Surface,
      absolute internal,
         NoDyn_Cylinder,
            position, reference, Ref_CylinderEdgeR, null,            # relative arm
         NoSta_Surface,
            position, null,                                          # relative arm
         single, reference, Ref_Surface, 0., 0., 1., string, "FN_R"; # force value
         
   force: FoStrin_Normal_CylinderEdgeL_Surface,
      absolute internal,
         NoDyn_Cylinder,
            position, reference, Ref_CylinderEdgeL, null,            # relative arm
         NoSta_Surface,
            position, null,                                          # relative arm
         single, reference, Ref_Surface, 0., 0., 1., string, "FN_L"; # force value
         
   force: FoStrin_Friction_CylinderEdgeR_Surface,
      absolute internal,
         NoDyn_Cylinder,
            position, reference, Ref_CylinderEdgeR, null,   # relative arm
         NoSta_Surface,
            position, null,                                 # relative arm
         single, reference, Ref_Surface, 0., 1., 0.,
            string, "-model::sf::cubsign(VY_R/Vt)*mu*FN_R"; # force value
            
   force: FoStrin_Friction_CylinderEdgeL_Surface,
      absolute internal,
         NoDyn_Cylinder,
            position, reference, Ref_CylinderEdgeL, null,   # relative arm
         NoSta_Surface,
            position, null,                                 # relative arm
         single, reference, Ref_Surface, 0., 1., 0.,
            string, "-model::sf::cubsign(VY_L/Vt)*mu*FN_L"; # force value
                           
   #-----------------------------------------------------------------------------
   # Gravity 
   gravity: 0., 0., -1., const, 9.81;
   
end: elements;