examples/Binder/1_Timoshenko_Beam.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Timoshenko Beam Example\n",
    "\n",
    "This Elastica tutorial explains the basics of setting up and running a simple simulation of rods in Elastica. Elastica simulates Cosserat Rods, which are thin, 1-dimensional rods that undergo all possible modes of deformation. This example considers a Timoshenko beam, which is the deformation of a beam under a constant applied force while accounting for shear deforation and rotational bending. This is a good example of the capabilities of Elastica and Cosserat Rods as it requires accounting for the effects of shear deformation, something that the classical Euler-Bernoulli beam solution does not.\n",
    "\n",
    "![timoshenko_beam_figure.png](../../assets/timoshenko_beam_figure.png)\n",
    "\n",
    "## Getting Started\n",
    "To set up the simulation, the first thing you need to do is import the necessary classes. Here we will only import the classes that we need. The `elastica.wrappers` classes make it easy to construct different simulation systems. Along with these wrappers, we need to import a rod class, classes for the boundary conditions, and time-stepping functions. As a note, this method of explicitly importing all classes can be a bit cumbersome. Future releases will simplify this step. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Import Wrappers\n",
    "from elastica.wrappers import BaseSystemCollection, Constraints, Forcing\n",
    "\n",
    "# Import Cosserat Rod Class\n",
    "from elastica.rod.cosserat_rod import CosseratRod\n",
    "\n",
    "# Import Boundary Condition Classes\n",
    "from elastica.boundary_conditions import OneEndFixedRod, FreeRod\n",
    "from elastica.external_forces import EndpointForces\n",
    "\n",
    "# Import Timestepping Functions\n",
    "from elastica.timestepper.symplectic_steppers import PositionVerlet\n",
    "from elastica.timestepper import integrate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have imported all the necessary classes, we want to create our beam system. We do this by combining all the wrappers we need to represent the physics that we to include in the simulation. In this case, that is the `BaseSystemCollection`, `Constraint` and `Forcings` because the simulation will consider a rod that is fixed in place on one end, and subject to an applied force on the other end. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing):\n",
    "    pass\n",
    "\n",
    "\n",
    "timoshenko_sim = TimoshenkoBeamSimulator()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating Rods\n",
    "With our simulator set up, we can now define the numerical, material, and geometric properties. \n",
    "\n",
    "First we define the number of elements in the rod. Next, the material properties are defined for every rod. These are the Young's modulus, the Poisson ratio, the density and the viscous damping coefficient. Finally, the geometry of the rod also needs to be defined by specifying the location of the rod and its orientation, length and radius. \n",
    "\n",
    "All of the values defined here are done in SI units, though this is not strictly necessary. You can rescale properties however you want, as long as you use consistent units throughout the simulation. See [here](https://info.simuleon.com/blog/units-in-abaqus) for an example of consistent units.\n",
    "\n",
    "In order to make the difference between a shearable and unshearable rod more clear, we are using a Poisson ratio of 99. This is an unphysical value, as Poisson ratios can not exceed 0.5, however, it is used here for demonstration purposes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting up test params\n",
    "n_elem = 100\n",
    "\n",
    "density = 1000\n",
    "nu = 0.1\n",
    "E = 1e6\n",
    "# For shear modulus of 1e4, nu is 99!\n",
    "poisson_ratio = 99\n",
    "shear_modulus = E / (poisson_ratio + 1.0)\n",
    "\n",
    "start = np.zeros((3,))\n",
    "direction = np.array([0.0, 0.0, 1.0])\n",
    "normal = np.array([0.0, 1.0, 0.0])\n",
    "base_length = 3.0\n",
    "base_radius = 0.25\n",
    "base_area = np.pi * base_radius ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With all of the rod's parameters set, we can now create a rod with the specificed properties and add the rod to the simulator system. **Important:** Make sure that any rods you create get added to the simulator system (`timoshenko_sim`), otherwise they will not be included in your simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shearable_rod = CosseratRod.straight_rod(\n",
    "    n_elem,\n",
    "    start,\n",
    "    direction,\n",
    "    normal,\n",
    "    base_length,\n",
    "    base_radius,\n",
    "    density,\n",
    "    nu,\n",
    "    E,\n",
    "    shear_modulus=shear_modulus,\n",
    ")\n",
    "\n",
    "timoshenko_sim.append(shearable_rod)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding Boundary Conditions\n",
    "With the rod added to the system, we need to apply boundary conditions. The first condition we will apply is fixing the location of one end of the rod. We do this using the `.constrain()` option and the `OneEndFixedRod` boundary condition. We are modifying the `timoshenko_sim` simulator to `constrain` the `shearable_rod` object using the `OneEndFixedRod` type of constraint. \n",
    "\n",
    "We also need to define which node of the rod is being constrained. We do this by passing the index of the nodes that we want to constain to `constrained_position_idx`. Here we are fixing the first node in the rod. In order to keep the rod from rotating around the fixed node, we also need to constrain an element between two nodes. This fixes the orientation of the rod. We do this by passing the index of the element that we want to fix to `constrained_director_idx`. Like with the position, we are fixing the first element of the rod. Together, this contrains the position and orientation of the rod at the origin. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timoshenko_sim.constrain(shearable_rod).using(\n",
    "    OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)\n",
    ")\n",
    "print(\"One end of the rod is now fixed in place\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next boundary condition that we want to apply is the endpoint force. Similarly to how we constrained one of the points, we want the `timoshenko_sim` simulator to `add_forcing_to` the `shearable_rod` object using the `EndpointForces` type of forcing. This `EndpointForces` applies forces to both ends of the rod. We want to apply a negative force in the $d_1$ direction, but only at the end of the rod. We do this by specifying the force vector to be applied at each end as `origin_force` and `end_force`. We also want to ramp up the force over time, so we make the force take some `ramp_up_time` to reach its steady-state value. This helps avoid numerical errors due to discontinuities in the applied force. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "origin_force = np.array([0.0, 0.0, 0.0])\n",
    "end_force = np.array([-10.0, 0.0, 0.0])\n",
    "ramp_up_time = 5.0\n",
    "\n",
    "timoshenko_sim.add_forcing_to(shearable_rod).using(\n",
    "    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time\n",
    ")\n",
    "print(\"Forces added to the rod\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add Unshearable Rod\n",
    "\n",
    "Along with the shearable rod, we also want to add an unshearable rod to be able to compare the difference between the two. We do this the same way we did for the first rod, however, because this rod is unsherable, we need to change the Poisson ratio to make the rod unsherable. For a truely unsheraable rod, you would need a Poisson ratio of -1.0, however, this causes the system to be numerically unstable, so instead we make the system nearly unshearable by using a Poisson ratio of -0.85. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start into the plane\n",
    "unshearable_start = np.array([0.0, -1.0, 0.0])\n",
    "unshearable_rod = CosseratRod.straight_rod(\n",
    "    n_elem,\n",
    "    unshearable_start,\n",
    "    direction,\n",
    "    normal,\n",
    "    base_length,\n",
    "    base_radius,\n",
    "    density,\n",
    "    nu,\n",
    "    E,\n",
    "    # Unshearable rod needs G -> inf, which is achievable with a poisson ratio of -1.0\n",
    "    shear_modulus=E / (-0.85 + 1.0),\n",
    ")\n",
    "\n",
    "timoshenko_sim.append(unshearable_rod)\n",
    "timoshenko_sim.constrain(unshearable_rod).using(\n",
    "    OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)\n",
    ")\n",
    "\n",
    "timoshenko_sim.add_forcing_to(unshearable_rod).using(\n",
    "    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time\n",
    ")\n",
    "print(\"Unshearable rod set up\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System Finalization\n",
    "\n",
    "We have now added all the necessary rods and boundary conditions to our system. The last thing we need to do is finalize the system. This goes through the system, rearranges things, and precomputes useful quantities to prepare the system for simulation. \n",
    "\n",
    "As a note, if you make any changes to the rod after calling finalize, you will need to re-setup the system. This requires rerunning all cells above this point. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "timoshenko_sim.finalize()\n",
    "print(\"System finalized\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Simulation Time\n",
    "\n",
    "The last thing we need to do deceide how long we want the simulation to run for and what timestepping method to use. Currently, the PositionVerlet algorithim is suggested default method. \n",
    "\n",
    "In this example, we are trying to match a steady-state solution by temporally evolving our system to reach equillibrium. As such, there is a tradeoff between letting the simulation run long enough to each the equillibrium and waiting around for the simulation to be done. Here we are running the simulation for 10 seconds, this produces reasonable agreement with the analytical solution without taking to long to finish. If you run the simulation for longer, you will get better agreement with the analytical solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "final_time = 10.0\n",
    "dl = base_length / n_elem\n",
    "dt = 0.01 * dl\n",
    "total_steps = int(final_time / dt)\n",
    "print(\"Total steps to take\", total_steps)\n",
    "\n",
    "timestepper = PositionVerlet()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Simulation\n",
    "We are now ready to perform the simulation. To run the simulation, we `integrate` the `timoshenko_sim` system using the `timestepper` method until `final_time` by taking `total_steps`. As currently setup, the beam simulation takes about 1 minute to run. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "integrate(timestepper, timoshenko_sim, final_time, total_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Post Processing Results\n",
    "Now that we have finished the simulation, we want to post-process the results. We will do this by comparing the solutions for the shearable and unshearable beams with the analytical Timoshenko and Euler-Bernoulli beam results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute beam position for sherable and unsherable beams.\n",
    "def analytical_result(arg_rod, arg_end_force, shearing=True, n_elem=500):\n",
    "    base_length = np.sum(arg_rod.rest_lengths)\n",
    "    arg_s = np.linspace(0.0, base_length, n_elem)\n",
    "    if type(arg_end_force) is np.ndarray:\n",
    "        acting_force = arg_end_force[np.nonzero(arg_end_force)]\n",
    "    else:\n",
    "        acting_force = arg_end_force\n",
    "    acting_force = np.abs(acting_force)\n",
    "    linear_prefactor = -acting_force / arg_rod.shear_matrix[0, 0, 0]\n",
    "    quadratic_prefactor = (\n",
    "        -acting_force\n",
    "        / 2.0\n",
    "        * np.sum(arg_rod.rest_lengths / arg_rod.bend_matrix[0, 0, 0])\n",
    "    )\n",
    "    cubic_prefactor = (acting_force / 6.0) / arg_rod.bend_matrix[0, 0, 0]\n",
    "    if shearing:\n",
    "        return (\n",
    "            arg_s,\n",
    "            arg_s * linear_prefactor\n",
    "            + arg_s ** 2 * quadratic_prefactor\n",
    "            + arg_s ** 3 * cubic_prefactor,\n",
    "        )\n",
    "    else:\n",
    "        return arg_s, arg_s ** 2 * quadratic_prefactor + arg_s ** 3 * cubic_prefactor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to plot the results. The one thing to point out in this function is how to access the position of the rods. They are located in `rod.position_collection[dim, n_elem]`. In this case, we are plotting the x- and z-dimensions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_timoshenko(shearable_rod, unshearable_rod, end_force):\n",
    "    import matplotlib.pyplot as plt\n",
    "\n",
    "    analytical_shearable_positon = analytical_result(\n",
    "        shearable_rod, end_force, shearing=True\n",
    "    )\n",
    "    analytical_unshearable_positon = analytical_result(\n",
    "        unshearable_rod, end_force, shearing=False\n",
    "    )\n",
    "\n",
    "    fig = plt.figure(figsize=(5, 4), frameon=True, dpi=150)\n",
    "    ax = fig.add_subplot(111)\n",
    "    ax.grid(b=True, which=\"major\", color=\"grey\", linestyle=\"-\", linewidth=0.25)\n",
    "\n",
    "    ax.plot(\n",
    "        analytical_shearable_positon[0],\n",
    "        analytical_shearable_positon[1],\n",
    "        \"k--\",\n",
    "        label=\"Timoshenko\",\n",
    "    )\n",
    "    ax.plot(\n",
    "        analytical_unshearable_positon[0],\n",
    "        analytical_unshearable_positon[1],\n",
    "        \"k-.\",\n",
    "        label=\"Euler-Bernoulli\",\n",
    "    )\n",
    "\n",
    "    ax.plot(\n",
    "        shearable_rod.position_collection[2, :],\n",
    "        shearable_rod.position_collection[0, :],\n",
    "        \"b-\",\n",
    "        label=\"n=\" + str(shearable_rod.n_elems),\n",
    "    )\n",
    "    ax.plot(\n",
    "        unshearable_rod.position_collection[2, :],\n",
    "        unshearable_rod.position_collection[0, :],\n",
    "        \"r-\",\n",
    "        label=\"n=\" + str(unshearable_rod.n_elems),\n",
    "    )\n",
    "\n",
    "    ax.legend(prop={\"size\": 12})\n",
    "    ax.set_ylabel(\"Y Position (m)\", fontsize=12)\n",
    "    ax.set_xlabel(\"X Position (m)\", fontsize=12)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "plot_timoshenko(shearable_rod, unshearable_rod, end_force)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the sake of time, we are stopping this simulation early. This leads to some disagreement between the analytical solution and the Elastica solution as there are still some transient effects in the Elastica solution. Allowing the simulation to run longer will lead to a closer result between the analytical and Elastica solutions. \n",
    "\n",
    "## Dynamic Plotting \n",
    "To illustrate how the system evolves over time, we can also plot the system in time. To do this, we need to recreate the system, which we now call `BeamSimulator`. It is the same as the previous system so we will just write everything very compactly to save space. We also slightly modify our plotting and integrating functions to allow the output to be plotting during the simulation. \n",
    "\n",
    "Since we will be plotting the system over time, we also need to initalize the time at 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "time = 0.0\n",
    "\n",
    "\n",
    "class BeamSimulator(BaseSystemCollection, Constraints, Forcing):\n",
    "    pass\n",
    "\n",
    "\n",
    "dynamic_update_sim = BeamSimulator()\n",
    "\n",
    "shearable_rod_new = CosseratRod.straight_rod(\n",
    "    n_elem,\n",
    "    start,\n",
    "    direction,\n",
    "    normal,\n",
    "    base_length,\n",
    "    base_radius,\n",
    "    density,\n",
    "    nu,\n",
    "    E,\n",
    "    poisson_ratio,\n",
    ")\n",
    "dynamic_update_sim.append(shearable_rod_new)\n",
    "dynamic_update_sim.constrain(shearable_rod_new).using(\n",
    "    OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)\n",
    ")\n",
    "dynamic_update_sim.add_forcing_to(shearable_rod_new).using(\n",
    "    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time\n",
    ")\n",
    "\n",
    "unshearable_rod_new = CosseratRod.straight_rod(\n",
    "    n_elem,\n",
    "    unshearable_start,\n",
    "    direction,\n",
    "    normal,\n",
    "    base_length,\n",
    "    base_radius,\n",
    "    density,\n",
    "    nu,\n",
    "    E,\n",
    "    poisson_ratio=-0.85,\n",
    ")\n",
    "dynamic_update_sim.append(unshearable_rod_new)\n",
    "dynamic_update_sim.constrain(unshearable_rod_new).using(\n",
    "    OneEndFixedRod, constrained_position_idx=(0,), constrained_director_idx=(0,)\n",
    ")\n",
    "dynamic_update_sim.add_forcing_to(unshearable_rod_new).using(\n",
    "    EndpointForces, origin_force, end_force, ramp_up_time=ramp_up_time\n",
    ")\n",
    "\n",
    "dynamic_update_sim.finalize()\n",
    "\n",
    "\n",
    "def run_and_update_plot(simulator, dt, start_time, stop_time, ax):\n",
    "    from elastica.timestepper import extend_stepper_interface\n",
    "    from elastica.timestepper.symplectic_steppers import PositionVerlet\n",
    "\n",
    "    timestepper = PositionVerlet()\n",
    "    do_step, stages_and_updates = extend_stepper_interface(timestepper, simulator)\n",
    "\n",
    "    n_steps = int((stop_time - start_time) / dt)\n",
    "    time = start_time\n",
    "    for i in range(n_steps):\n",
    "        time = do_step(timestepper, stages_and_updates, simulator, time, dt)\n",
    "    plot_timoshenko_dynamic(shearable_rod_new, unshearable_rod_new, end_force, time, ax)\n",
    "    return time\n",
    "\n",
    "\n",
    "def plot_timoshenko_dynamic(shearable_rod, unshearable_rod, end_force, time, ax):\n",
    "    import matplotlib.pyplot as plt\n",
    "    from IPython import display\n",
    "\n",
    "    analytical_shearable_positon = analytical_result(\n",
    "        shearable_rod, end_force, shearing=True\n",
    "    )\n",
    "    analytical_unshearable_positon = analytical_result(\n",
    "        unshearable_rod, end_force, shearing=False\n",
    "    )\n",
    "\n",
    "    ax.clear()\n",
    "    ax.grid(b=True, which=\"major\", color=\"grey\", linestyle=\"-\", linewidth=0.25)\n",
    "    ax.plot(\n",
    "        analytical_shearable_positon[0],\n",
    "        analytical_shearable_positon[1],\n",
    "        \"k--\",\n",
    "        label=\"Timoshenko\",\n",
    "    )\n",
    "    ax.plot(\n",
    "        analytical_unshearable_positon[0],\n",
    "        analytical_unshearable_positon[1],\n",
    "        \"k-.\",\n",
    "        label=\"Euler-Bernoulli\",\n",
    "    )\n",
    "\n",
    "    ax.plot(\n",
    "        shearable_rod.position_collection[2, :],\n",
    "        shearable_rod.position_collection[0, :],\n",
    "        \"b-\",\n",
    "        label=\"shearable rod\",\n",
    "    )\n",
    "    ax.plot(\n",
    "        unshearable_rod.position_collection[2, :],\n",
    "        unshearable_rod.position_collection[0, :],\n",
    "        \"r-\",\n",
    "        label=\"unshearable rod\",\n",
    "    )\n",
    "\n",
    "    ax.legend(prop={\"size\": 12}, loc=\"lower left\")\n",
    "    ax.set_ylabel(\"Y Position (m)\", fontsize=12)\n",
    "    ax.set_xlabel(\"X Position (m)\", fontsize=12)\n",
    "    ax.set_title(\"Simulation Time: %0.2f seconds\" % time)\n",
    "    ax.set_xlim([-0.1, 3.1])\n",
    "    ax.set_ylim([-0.045, 0.002])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can run the simulation for a time interval `evolve_for_time` and have the system be plotted every `update_interval`. If you run the cell multiple times in a row, you will see that the that the system continues to evolve in time, this is because you are continually updating `dynamic_update_sim`. If you want to reset the system back to its original configuration, run the cell above this one. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython import display\n",
    "\n",
    "evolve_for_time = 10.0\n",
    "update_interval = 1.0e-1\n",
    "\n",
    "# update the plot every 1 second\n",
    "fig = plt.figure(figsize=(5, 4), frameon=True, dpi=150)\n",
    "ax = fig.add_subplot(111)\n",
    "first_interval_time = update_interval + time\n",
    "last_interval_time = time + evolve_for_time\n",
    "for stop_time in np.arange(\n",
    "    first_interval_time, last_interval_time + dt, update_interval\n",
    "):\n",
    "    time = run_and_update_plot(dynamic_update_sim, dt, time, stop_time, ax)\n",
    "    display.clear_output(wait=True)\n",
    "    display.display(plt.gcf())\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Important note on saving data:\n",
    "This current method of plotting data during the simulation is helpful for visualizing how the system evolves, but it is computationally inefficient as we are constantly pausing the simulation to plot. It also does not save data for additional post-processing later. A better method for saving data from a simulation is to use call-back functions. There is information on how to use these functions in the [snake tutorial](./2_Slithering_Snake.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}