Full script for defining and running MBDyn model

main.py

#!/usr/bin/env python3

# Extracted from seat-dynamics @ 240d093c23b523e262fb7a617f84673ca0fdaa74
# Written for tutorial chapter 3 "Rigid Pendulum", see:
# https://raw.githubusercontent.com/mmorandi/MBDyn-web/main/userfiles/documents/tutorials.pdf

from argparse import ArgumentDefaultsHelpFormatter, ArgumentParser
from dataclasses import dataclass
from subprocess import run
from os import makedirs
from preprocess.contrib.PythonPreprocessor.MBDynLib import *
from preprocess.contrib.PythonPreprocessor.MBDynModel import *


@dataclass
class Arguments:
    verbose: bool
    name: str


def main():
    parser = ArgumentParser(formatter_class=ArgumentDefaultsHelpFormatter)
    _ = parser.add_argument(
        "-v", "--verbose", action="store_true", help="print more information"
    )
    _ = parser.add_argument("-n", "--name", help="mbdyn file name", default="seat")
    args = Arguments(**vars(parser.parse_args()))  # pyright: ignore[reportAny]

    model = create_model()
    if args.verbose:
        print_lines(str(model))

    makedirs("./out", exist_ok=True)
    # When passing the model over stdin, the output files didn't exist in the container-mounted directory after running
    with open(f"./out/{args.name}.mbd", "w") as out_file:
        _ = out_file.write(str(model))

    _ = run(
        [
            "mbdyn",
            "-f",
            f"./out/{args.name}.mbd",
        ],
        check=True,
    )


def create_model() -> MBDynModel:
    model = MBDynModel(
        data=Data(problem="initial value"),
        problem=InitialValue(
            initial_time=0.0,
            final_time=1.0,
            time_step=1.0e-3,
            max_iterations=MaxIterations(max_iterations=10),
            tolerance=Tolerance(residual_tolerance=1.0e-6),
            derivatives_coefficient=DerivativesCoefficient(coefficient=1.0e-6),
        ),
        control_data=ControlData(),  # let it compute element counts
    )
    # gravity attempted later, tutorial had just this line:
    # gravity;

    # FIXME: MBDynModel should define all MBVars at the start of the input file,
    # unitl then, they are only usable in "begin preprocess" blocks
    pendulum_id = 1
    mass_id = 2
    mass = 1.0
    length = 0.5
    omega0 = 0.2

    pendulum_ref = Reference(
        idx=pendulum_id,
        position=Position(relative_position=NULL, reference="global"),
        orientation=Position(relative_position=EYE, reference="global"),
        velocity=Position(relative_position=NULL, reference="global"),
        angular_velocity=Position(relative_position=[0, omega0, 0], reference="global"),
    )
    mass_ref = Reference(
        idx=mass_id,
        position=Position(
            relative_position=[0, 0, -length], reference=pendulum_ref
        ),  # No MBVar mathematical operators
        orientation=Position(relative_position=EYE, reference=pendulum_ref),
        velocity=Position(relative_position=NULL, reference=pendulum_ref),
        angular_velocity=Position(relative_position=NULL, reference=pendulum_ref),
    )

    pendulum_node = Node(
        idx=1000 + pendulum_id,
        node_type="dynamic",
        position=Position(relative_position=NULL, reference=pendulum_ref),
        orientation=Position(relative_position=EYE, reference=pendulum_ref),
        velocity=Position(relative_position=NULL, reference=pendulum_ref),
        angular_velocity=Position(relative_position=NULL, reference=pendulum_ref),
    )
    model.add_node(pendulum_node)

    mass_node = Node(
        idx=2000 + mass_id,
        node_type="dynamic",
        position=Position(relative_position=NULL, reference=mass_ref),
        orientation=Position(relative_position=EYE, reference=mass_ref),
        velocity=Position(relative_position=NULL, reference=mass_ref),
        angular_velocity=Position(relative_position=NULL, reference=mass_ref),
    )
    model.add_node(mass_node)

    # no dynamic dofs (it will be fully grounded)
    static_pendulum_node = Node(
        idx=3000 + pendulum_id,
        node_type="static",
        position=Position(relative_position=NULL, reference=pendulum_ref),
        orientation=Position(relative_position=EYE, reference=pendulum_ref),
        velocity=Position(relative_position=NULL, reference="global"),
        angular_velocity=Position(
            relative_position=NULL, reference=pendulum_ref
        ),  # "global" means no angular velocity
    )
    model.add_node(static_pendulum_node)

    static_mass_node = Node(
        idx=3000 + mass_id,
        node_type="dynamic",
        position=Position(relative_position=NULL, reference=mass_ref),
        orientation=Position(relative_position=EYE, reference=mass_ref),
        velocity=Position(relative_position=NULL, reference="global"),
        angular_velocity=Position(
            relative_position=NULL, reference=mass_ref
        ),  # "global" means no angular velocity
    )
    model.add_node(static_mass_node)

    body = Body(
        idx=1000 + mass_id,
        node=pendulum_node,
        mass=mass,
        position=Position(relative_position=NULL, reference=mass_ref),
        inertial_matrix=NULL,
    )
    model.add_element(body)

    joint = RevolutePin(
        idx=1000 + mass_id,
        node_label=1000 + pendulum_id,
        relative_offset=Position(relative_position=NULL, reference=pendulum_ref),
        # This is probably the required argument absolute_pin_position
        # hinge, reference, Pendulum, 1, 1.,0.,0., 3, 0.,1.,0.;
        absolute_pin_position=Position(
            relative_position=[1, 0, 0], reference=pendulum_ref
        ),
        absolute_pin_orientation_mat=Position(
            relative_position=[0, 1, 0], reference=pendulum_ref
        ),  # is this the right reference? what is 3 above?
    )
    model.add_element(joint)

    body = Body(
        idx=2000 + mass_id,
        node=mass_node,
        mass=mass,
        position=Position(relative_position=NULL, reference=mass_ref),
        inertial_matrix=NULL,
    )
    model.add_element(body)

    joint = RevolutePin(
        idx=2000 + mass_id,
        node_label=2000 + mass_id,
        relative_offset=Position(relative_position=NULL, reference=pendulum_ref),
        # This is probably the required argument absolute_pin_position
        # hinge, reference, Pendulum, 1, 1.,0.,0., 3, 0.,1.,0.;
        absolute_pin_position=Position(
            relative_position=[1, 0, 0], reference=pendulum_ref
        ),
        absolute_pin_orientation_mat=Position(
            relative_position=[0, 1, 0], reference=pendulum_ref
        ),  # is this the right reference? what is 3 above?
    )
    model.add_element(joint)

    body = Body(
        idx=3000 + mass_id,
        node=static_mass_node,
        mass=mass,
        position=Position(relative_position=NULL, reference=mass_ref),
        inertial_matrix=EYE,  # Otherwise the problem will be singular... But why?
    )
    model.add_element(body)

    joint = Clamp(
        idx=3000 + pendulum_id,
        node=static_pendulum_node,
        position="node",
        orientation_mat="node",
    )
    model.add_element(joint)

    joint = Distance(
        idx=3000 + mass_id,
        node_1_label=3000 + pendulum_id,
        node_2_label=3000 + mass_id,
        distance=ConstDriveCaller(const_value=length),
    )
    model.add_element(joint)

    # Element not implemented, I assume it can work like this
    gravity = StructuralForce(
        idx=5,
        ftype="absolute",
        node=pendulum_node,  # I guess???
        position=NULL,
        force_orientation=Position(relative_position=[0.0, 0.0, -1.0]),
        force_drive=[
            ConstDriveCaller(
                const_value=9.81,
            ),
        ],
    )
    model.add_element(gravity)

    return model


def print_lines(text: str) -> None:
    for i, line in enumerate(text.splitlines()):
        print(f"{i + 1:>3}: {line}")


if __name__ == "__main__":
    main()