<jupyter_start><jupyter_text>STANDALONE STUDY **SPATIAL FOURBAR**-----------------  Environment Setup--------------------  Colab Session*Run if in a Colab Session Only!*  Setting-up Colab Virtual MachineThis will clone the package repository and install the needed modules into the Virtual Machine dedicated for the running notebook session.<jupyter_code>!git clone https://github.com/khaledghobashy/smbd.git
!cd smbd/smbd/numenv/python/numerics/core/math_funcs/_cython_definitions && python setup.py build_ext -if
!pip install sympy==1.4<jupyter_output><empty_output><jupyter_text>Adding package path to the python system path and creating a new folder, **smbd_project**, where we can store the model files and results.<jupyter_code>import os
import sys

pkg_path = '/content/smbd'
sys.path.append(pkg_path)

os.mkdir('smbd_project')
os.chdir('smbd_project')<jupyter_output><empty_output><jupyter_text>Rendering SymPy in ColabRendering SymPy equations in Colab requires the MathJax javascript library to be available within each cell output. The following defines an execution hook that loads it automatically each time you execute a cell.<jupyter_code>from IPython.display import Math, HTML

def load_mathjax_in_cell_output():
  display(HTML("<script src='https://www.gstatic.com/external_hosted/"
               "mathjax/latest/MathJax.js?config=default'></script>"))
get_ipython().events.register('pre_run_cell', load_mathjax_in_cell_output)<jupyter_output><empty_output><jupyter_text>---------------------------------------------------------------  Local Session*Run if on your Local Machine!*The ```pkg_path``` variable is where you have the repository on your machine and should be changed to where you cloned the repository.*Uncomment the code below and change the ```pkg_path``` variable!*<jupyter_code>#import sys
#pkg_path = '/path/to/smbd'
#sys.path.append(pkg_path)<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------------------------------ Creating Project Directory<jupyter_code>from smbd.systems import standalone_project
project = standalone_project()
project.create()<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------------------------------  **SYSTEM DESCRIPTION**-------------------------  **Important Note**The tool uses [euler-parameters](https://en.wikibooks.org/wiki/Multibody_Mechanics/Euler_Parameters) -which is a 4D unit quaternion- to represents bodies orientation in space. This makes the generalized coordinates used to fully define a body in space to be **7,** instead of **6**, it also adds an algebraic equation to the constraints that ensures the unity/normalization of the body quaternion. This is an important remark as the calculations of the degrees-of-freedom depends on it.  **Summary**A four-bar linkage is the simplest movable closed chain linkage. It consists of four bodies, called bars or links, connected in a loop by four joints that gives the mechanism only 1 degree-of-freedom. More general information can be found on [wikipedia](https://en.wikipedia.org/wiki/Four-bar_linkage).  **Topology Layout**The mechanism consists of 3 Bodies + 1 Ground. Therefore, total system coordinates -including the ground- is $$n=n_b\times7 = 4\times7 = 28$$ where $n_b$ is the total number of bodies. The list of bodies is given below:- Crank $l_1$.- Connecting-Rod $l_2$.- Rocker $l_3$.The system connectivity is as follows:- Crank $l_1$ is connected to the ground by a revolute joint, resulting in constraint equations $n_{c,rev} = 5$- Connecting-Rod $l_2$ is connected to the Crank $l_1$ by a spherical joint, resulting in constraint equations $n_{c,sph} = 3$- Rocker $l_3$ is connected to Connecting-Rod $l_2$ by a universal joint, resulting in constraint equations $n_{c,uni} = 4$- Rocker $l_3$ is connected to the ground by a revolute joint, resulting in constraint equations $n_{c,tri} = 5$- Total free DOF is $$n-( n_{c,rev}+n_{c,sph}+n_{c,uni}+n_{c,P}+n_{c,g}) = 28-(10+3+4+(3\times1)+7) = 28-27 = 1$$, where the $n_{c,P}$ and $n_{c,g}$ represents the constraints due to euler-parameters normalization equations and the ground-constraints respectively.- The Crank revolute joint is actuated by a rotational actuator, resulting in constraint equations $n_{c,act} = 1$| Joint Name  | Body i         | Body j         | Joint Type | $n_c$ ||:-----------:|:-------------- |:-------------- | ---------- | ----- || a           | Ground         | Crank          | Revolute   | 5     || b           | Crank          | Connecting-Rod | Spherical  | 3     || c           | Connecting-Rod | Rocker         | Universal  | 4     || d           | Rocker         | Ground         | Revolute   | 5     | ------------------------------------------------------------------------------------------------------------------------------  **SYMBOLIC TOPOLOGY**------------------------ In this section, we create the symbolic topology that captures the topological layout that we just discussed.Defining the topology is very simple. We start by importing the ```standalone_topology``` class and create a new instance that represents our symbolic model. Then we start adding the components we discussed earlier, starting by the bodies, then the joints, actuators and forces, and thats it.These components will be represented symbolically, and therefore there is no need for any numerical inputs at this step.The system is stored in a form of a network graph that stores all the data needed for the assemblage of the system equations later. But even before the assemblage process, we can gain helpful insights about our system as well be shown.<jupyter_code>from smbd.systems import standalone_topology, configuration
project_name = 'spatial_fourbar'
sym_model = standalone_topology(project_name)<jupyter_output><empty_output><jupyter_text>Adding Bodies<jupyter_code>sym_model.add_body('l1')
sym_model.add_body('l2')
sym_model.add_body('l3')<jupyter_output><empty_output><jupyter_text>Adding Joints<jupyter_code>sym_model.add_joint.revolute('a','ground','rbs_l1')
sym_model.add_joint.spherical('b','rbs_l1','rbs_l2')
sym_model.add_joint.universal('c','rbs_l2','rbs_l3')
sym_model.add_joint.revolute('d','rbs_l3','ground')<jupyter_output><empty_output><jupyter_text>Adding Actuators<jupyter_code>sym_model.add_actuator.rotational_actuator('act', 'jcs_a')<jupyter_output><empty_output><jupyter_text>Adding Forces<jupyter_code>#model.add_force.internal_force('force_name', 'body_1', 'body_2', mirror=False)<jupyter_output><empty_output><jupyter_text>**SYMBOLIC CHARACTERISTICS**  Topology GraphVisualizing the connectivity of the system as a network graph, where the nodes represent the bodies, and the edges represent the joints, forces and/or actuators between the bodies. Checking the system's number of generalized coordinates $n$ and number of constraints $n_c$.  **ASSEMBLING** This is the last step of the symbolic building process, where we make the system starts the assemblage process of the governing equations, which will be used then in the code generation for the numerical simulation, and also can be used for further symbolic manipulations.*Note: The equations' notations will be discussed in the tool documentation files.*<jupyter_code>sym_model.assemble()<jupyter_output><empty_output><jupyter_text>Checking the System Equations ------------------------------------------------------------------------------------------------------------------------------  **SYSTEM CONFIGURATION**--------------------------- In this step we define a symbolic configuration of our symbolic topology. As you may have noticed in the symbolic topology building step, we only cared about the **_topology_**, thats is the system bodies and their connectivity, and we did not care explicitly with how these components are configured in space.In order to create a valid numerical simulation session, we have to provide the system with its numerical configuration needed, for example, the joints' locations and orientations. The symbolic topology in its raw form will require you to manually enter all these numerical arguments, which can be cumbersome even for smaller systems. This can be checked by checking the configuration inputs of the symbolic configuration as ```sym_config.config.input_nodes```Here we start by stating the symbolic inputs we wish to use instead of the default inputs set, and then we define the relation between these newly defined arguments and the original ones. **_The details of this process will be provided in the documentation._**<jupyter_code># Symbolic configuration name.
config_name = '%s_cfg'%project_name

# Symbolic configuration instance.
sym_config = configuration(config_name, sym_model)<jupyter_output><empty_output><jupyter_text>CONFIGURATION INPUTS<jupyter_code># Adding the desired set of UserInputs
# ====================================

sym_config.add_point.UserInput('a')
sym_config.add_point.UserInput('b')
sym_config.add_point.UserInput('c')
sym_config.add_point.UserInput('d')

sym_config.add_vector.UserInput('x')
sym_config.add_vector.UserInput('y')
sym_config.add_vector.UserInput('z')<jupyter_output><empty_output><jupyter_text>CONFIGURATION RELATIONS<jupyter_code># Defining Relations between original topology inputs
# and our desired UserInputs.
# ===================================================

# Revolute Joint (a) location and orientation
sym_config.add_relation.Equal_to('pt1_jcs_a', ('hps_a',))
sym_config.add_relation.Equal_to('ax1_jcs_a', ('vcs_x',))

# Spherical Joint (b) location and orientation
sym_config.add_relation.Equal_to('pt1_jcs_b', ('hps_b',))
sym_config.add_relation.Equal_to('ax1_jcs_b', ('vcs_z',))

# Universal Joint (c) location and orientation
sym_config.add_relation.Equal_to('pt1_jcs_c', ('hps_c',))
sym_config.add_relation.Oriented('ax1_jcs_c', ('hps_b', 'hps_c'))
sym_config.add_relation.Oriented('ax2_jcs_c', ('hps_c', 'hps_b'))

# Revolute Joint (d) location and orientation
sym_config.add_relation.Equal_to('pt1_jcs_d', ('hps_d',))
sym_config.add_relation.Equal_to('ax1_jcs_d', ('vcs_y',))<jupyter_output><empty_output><jupyter_text>GEOMETRIES Here we start defining basic geometric shapes that can represents the shapes of the bodies in our system. This serves two points:- Visualization and Animation in blender.- Evaluating the bodies inertia properties from these basic geometries instead of explicit definition.<jupyter_code># links radius
sym_config.add_scalar.UserInput('links_ro')

# Link 1 geometry
sym_config.add_geometry.Cylinder_Geometry('l1', ('hps_a','hps_b','s_links_ro'))
sym_config.assign_geometry_to_body('rbs_l1', 'gms_l1')

# Link 2 geometry
sym_config.add_geometry.Cylinder_Geometry('l2', ('hps_b','hps_c','s_links_ro'))
sym_config.assign_geometry_to_body('rbs_l2', 'gms_l2')

# Link 3 geometry
sym_config.add_geometry.Cylinder_Geometry('l3', ('hps_c','hps_d','s_links_ro'))
sym_config.assign_geometry_to_body('rbs_l3', 'gms_l3')<jupyter_output><empty_output><jupyter_text>Now lets check the current ```input_nodes```.  ASSEMBLING<jupyter_code>sym_config.export_JSON_file()<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------------------------------ ------------------------------------------------------------------------------------------------------------------------------  **CODE GENERATION**---------------------- This step aims to create a valid code that can be used for numerical simulation. We will use both the ***python code generator*** and the **C++ code generator** to create a valid numerical simulation environment in python and C++.Theoretically, the symbolic environment is uncoupled from the simulation environment, which opens the door to create various simulation environments that can be in any language.  Generating Python Numerical Environment<jupyter_code>from smbd.numenv.python.codegen import projects as py_numenv
# Creating the needed directories structure for the source code.
py_project = py_numenv.standalone_project(parent_dir='')
py_project.create_dirs()
# Generating the source code.
py_project.write_topology_code(sym_model.topology)
py_project.write_configuration_code(sym_config.config)<jupyter_output>File full path : numenv\python\src\spatial_fourbar.py
File full path : numenv\python\src\spatial_fourbar_cfg.py<jupyter_text>Generating C++ Numerical Environment<jupyter_code>from smbd.numenv.cpp_eigen.codegen import projects as cpp_numenv
cpp_project = cpp_numenv.standalone_project(parent_dir='')
cpp_project.create_dirs()
# Generating the source code.
cpp_project.write_topology_code(sym_model.topology)
cpp_project.write_configuration_code(sym_config.config)
# Generating the makefile and main.cpp files.
cpp_project.write_mainfile()
cpp_project.write_makefile()<jupyter_output>File full path : numenv\cpp_eigen\main.cpp
File full path : numenv\cpp_eigen\Makefile<jupyter_text>Generating Blender Script<jupyter_code>from smbd.utilities.blender import codegen
blender_code = codegen.script_generator(sym_config.config)
blender_code.write_code_file('numenv/')<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------------------------------  **PYTHON NUMERICAL SIMULATION**--------------------------- The **numerical simulation** step is also a straight forward one. We start by importing helper libraries like **numpy** and **pandas**, then we import the ```multibody_system``` class and the ```simulation``` class that will be used to create a numerical system and a simulation instance respectively. And finally we import the generated code files from the previous step.The steps then can be summarized as follows:1. Create a numerical model instance by passing the generated topology code to the ```multibody_system``` class.2. Make the configuration of that instance refers to the generated configuration code.3. Setting the numerical data of the configuration.4. Setting the **user-defined functions** that will be used for the actuators/forces.5. Creating a **simulation instance**, set the desired time-sampling and run the simulation.6. Plotting the results.   Imports<jupyter_code>import numpy as np
from smbd.numenv.python.numerics.core.systems import multibody_system, simulation, configuration
# Importing the generated modules/souce-code of the model.
from numenv.python.src import spatial_fourbar<jupyter_output><empty_output><jupyter_text>Creating Numerical Model<jupyter_code>num_model  = multibody_system(spatial_fourbar)

num_config = configuration('base')
num_config.construct_from_json('spatial_fourbar_cfg.json')

num_model.topology.config = num_config<jupyter_output><empty_output><jupyter_text>Setting Configuration Numerical Data<jupyter_code>num_config.hps_a.flat[:] = 0, 0, 0
num_config.hps_b.flat[:] = 0, 0, 200
num_config.hps_c.flat[:] = -750, -850, 650
num_config.hps_d.flat[:] = -400, -850, 0

num_config.vcs_x.flat[:] = 1, 0, 0
num_config.vcs_y.flat[:] = 0, 1, 0
num_config.vcs_z.flat[:] = 0, 0, 1

num_config.s_links_ro = 20
num_config.assemble()<jupyter_output><empty_output><jupyter_text>Setting User Functions.  Actuators Here we set the actuation functions of the system actuators.<jupyter_code>num_config.UF_mcs_act = lambda t : -np.deg2rad(360)*t<jupyter_output><empty_output><jupyter_text>Assembling the numerical configuration.Here, we pass the numerical inputs instance ```num_inputs``` to the ```assemble()``` method of the ```num_config``` instance of our model.  Creating Simulation Instance A simulation instance takes in three arguments, the simulation ```name```, the numerical ```model``` and the simulation ```typ```. The simulation ```typ``` argument should be either *kds* for kinematically driven systems or *dds* for dynamically driven systems. For our case, the system is fully defined kinematically, and therefore we will set the ```typ``` as *kds*.<jupyter_code>sim = simulation('sim', num_model, 'kds')
sim.set_time_array(5, 2e-2)
sim.solve()
sim.soln.pos_dataframe.plot(x='time', y='rbs_l3.x', grid=True, figsize=(10,4))
sim.soln.vel_dataframe.plot(x='time', y='rbs_l3.x', grid=True, figsize=(10,4))
sim.soln.acc_dataframe.plot(x='time', y='rbs_l3.x', grid=True, figsize=(10,4))<jupyter_output><empty_output><jupyter_text>Evaluating System Reactions Here we take the simulation a step further by making it an inverse dynamics problem and evaluating the reactions resulted due the imposed constraints, either by joints or actuators. This step is only valid if you have defined the numerical values of the inertia properties of the system.<jupyter_code>sim.eval_reactions()
sim.soln.reactions_dataframe.plot(x='time', y=['T_ground_mcs_act.x', 'T_ground_jcs_a.x'],
                                   grid=True, figsize=(10,4))<jupyter_output><empty_output><jupyter_text>------------------------------------------------------------------------------------------------------------------------------  **C++ NUMERICAL SIMULATION**_Only for Linux Machines_------------------------------- ***TO DO BRIEF*** The C++ numerical environment provides samples of **Makefile** and **main.cpp** files that can be used to create, compile and run an executable of the created model.In order to create a valid environment, we have to edit the **Makefile** and the **main.cpp** files to match the model variables as shown below. These edits can be carried out on the text files using any text editor normally, but for completeness and compatibility with the **Jupyter Notebook** environment the files can be edited here locally using the *Jupyter magic functions*.First, we check the current working directory, it should be the parent directory that contains [numenv, config_inputs, results, ..], then we use the command```%cd numenv/cpp_eigen/``` to change the working directory.<jupyter_code># Changing the current working directory to 'numenv/cpp_eigen/'
%cd numenv/cpp_eigen/<jupyter_output>/home/khaledghobashy/Documents/smbd/examples/standalone_models/notebooks/spatial_fourbar/numenv/cpp_eigen<jupyter_text>---------------------------------------------------  **Editing the main.cpp file** Using the ```%load ``` magic function.-------------------------------------------Then we use the command ```%load main.cpp```, this will load the generated **main.cpp** content in the current cell where we can edit the content. The edits are as follows:- Remove the first line ``` %load main.cpp``` and insert the line ```%%writefile new_main.cpp```.- Change the variables topology and configuration in the include directive to match the names of the model generated source files.- Start specifying the configuration numerical data via the ```Config``` class.    After these edits, just run the cell like any normal cell via Shift+Enter.  Using the ```%cat ``` magic function.-------------------------------------------Then we use the command ```%cat main.cpp```, this will load the generated **main.cpp** content as the output of the current cell, we can then copy the content and past and edit it in another cell. The edits are as follows:- Insert the line ```%%writefile new_main.cpp``` at the top of the cell.- Change the variables topology and configuration in the include directive to match the names of the model generated source files.- Start specifying the configuration numerical data via the ```Config``` class.    After these edits, just run the cell like any normal cell via Shift+Enter.<jupyter_code>%cat main.cpp
%%writefile new_main.cpp

#include <iostream>

#include "smbd/solvers.hpp"

#include "src/spatial_fourbar.hpp"
#include "src/spatial_fourbar_cfg.hpp"


int main()
{
    Topology model("");
    auto Config = ConfigurationInputs<Configuration>(model.config);

    // assign the configuration inputs needed ...
    
    Config.R_ground << 0, 0, 0 ;
	Config.P_ground << 1, 0, 0, 0 ;

	Config.hps_a << 0, 0, 0 ;
	Config.hps_b << 0, 0, 200 ;
	Config.hps_c << -750, -850, 650 ;
	Config.hps_d << -400, -850, 0 ;

	Config.vcs_x << 1, 0, 0 ;
	Config.vcs_y << 0, 1, 0 ;
	Config.vcs_z << 0, 0, 1 ;

	Config.s_links_ro = 20 ;
	Config.UF_mcs_act = [](double t)->double{return 2*(22/7)*t;};

    Config.assemble();

    Solver<Topology> Soln(model);
    Soln.set_time_array(5, 250);
    Soln.Solve();
    Soln.ExportResultsCSV("../../results/", 0);

};<jupyter_output>Overwriting new_main.cpp<jupyter_text>---------------------------------------------------  **Editing the Makefile file**Following the same steps done for the main.cpp file, then apply the following edits:- Change the variables MODEL, CONFIG and MAIN to match the generated source files of the topology, configuration and the edited new_main.cpp file.    After these edits, just run the cell like any normal cell via Shift+Enter.<jupyter_code>%cat Makefile
%%writefile Makefile

# Change MODEL, CONFG and MAIN to match the source files you want to build
# ========================================================================
MODEL := spatial_fourbar
CONFG := spatial_fourbar_cfg
MAIN := new_main.cpp
# ========================================================================

M_BUILD := build/
M_SRC := src/
M_BIN := bin/

NUM_DIR := /home/khaledghobashy/Documents/smbd/smbd/numenv/cpp_eigen/numerics

SMBD_SRC := $(NUM_DIR)/src
SMBD_BUILD := $(NUM_DIR)/build

SMBD_OBJS = $(SMBD_BUILD)/*.o

DEPS := $(M_BUILD)$(MODEL).o $(MAIN) $(M_SRC)$(CONFG).hpp $(SMBD_SRC)/smbd/solvers.hpp

INC := -I $(SMBD_SRC)
CC := g++


$(M_BIN)$(MODEL): $(DEPS) $(SMBD_OBJS)
	$(CC) $(INC) $(M_BUILD)$(MODEL).o $(MAIN) $(SMBD_OBJS) -o $@

$(M_BUILD)$(MODEL).o: $(M_SRC)$(MODEL).cpp $(M_SRC)$(MODEL).hpp
	$(CC) $(INC) -c -o $@ $<


$(SMBD_BUILD)/%.o: $(SMBD_SRC)/smbd/%.cpp $(SMBD_SRC)/smbd/%.hpp
	cd $(SMBD_SRC)/../ && make


clear:
	rm $(M_BUILD)*.o $(M_BIN)$(MODEL)<jupyter_output>Overwriting Makefile<jupyter_text>Running ```make``` command<jupyter_code>!make<jupyter_output>g++ -I /home/khaledghobashy/Documents/smbd/smbd/numenv/cpp_eigen/numerics/src -c -o build/spatial_fourbar.o src/spatial_fourbar.cpp
g++ -I /home/khaledghobashy/Documents/smbd/smbd/numenv/cpp_eigen/numerics/src build/spatial_fourbar.o new_main.cpp /home/khaledghobashy/Documents/smbd/smbd/numenv/cpp_eigen/numerics/build/*.o -o bin/spatial_fourbar<jupyter_text>Running the executable.<jupyter_code>!bin/spatial_fourbar<jupyter_output>Starting Solver ...
[=================================================>] 249 

Positions results saved as : ../../results/Positions.csv<jupyter_text>Loading and Plotting Results.<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
data = pd.read_csv('../../results/Positions.csv')
data.plot(x='time', y='rbs_l3.x', grid=True, figsize=(10,4))
data.plot(x='time', y='rbs_l1.y', grid=True, figsize=(10,4))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17592592592592593, 'fraction_numerical': 0.008822501747030049, 'mean_word_length': 3.9786910197869103, 'pattern_counts': {'":': 0, '<': 142, '<?xml version=': 0, '>': 125, 'https://': 4, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '31', 'dir': 'jupyter-structured-clean-dedup', 'id': '494435', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 6913, 'n_tokens_neox': 6581, 'n_words': 2092, 'path': 'uraeus/smbd/utilities/serialization/structural/test/spatial_fourbar.ipynb', 'repo_name': 'khaledghobashy/asurt_mbdt', 'size': '22896'}