tutorial5_1.ipynb

{"metadata":{"interpreter":{"hash":"8661b582be0f012a274b7385348a70dad23ce1ee0d922b32977c503b98c884f2"},"kernelspec":{"name":"python3","display_name":"Python 3 (ipykernel)","language":"python"},"language_info":{"name":"python","version":"3.7.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Tutorial 7: Modal superposition of a jacket wind turbine\n\nIn this tutorial we will use the model defined in tutorial 6 to apply the modal superposition technique in a simplified jacket wind turbine. We assume that the structure can be modelled as a space frame with structural elements that are subject to axial displacement and bending. Therefore, the equations of motion at a given element are given by: \n\n$$ \\rho A \\frac{\\partial^2 u (x,t)}{\\partial t^2} - E A \\frac{\\partial^2 u(x,t)}{\\partial x^2} = q_u(x) $$\n\n$$ \\rho A \\frac{\\partial^2 v (x,t)}{\\partial t^2} + E I \\frac{\\partial^4 v(x,t)}{\\partial x^4} = q_v(x) $$\n\nIn the following figure we have a sketch of the simplified geometry that will be used in this example.\n\n![figure](tutorial8_figure.png)\n","metadata":{}},{"cell_type":"markdown","source":"## Step 0: input parameters\n\nLet's start by setting some parameters:","metadata":{}},{"cell_type":"code","source":"pip install scipy","metadata":{"trusted":true},"execution_count":1,"outputs":[{"name":"stdout","text":"Collecting scipy\n  Downloading scipy-1.7.3-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (38.1 MB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m38.1/38.1 MB\u001b[0m \u001b[31m30.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m00:01\u001b[0m00:01\u001b[0m\n\u001b[?25hCollecting numpy<1.23.0,>=1.16.5\n  Downloading numpy-1.21.6-cp37-cp37m-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (15.7 MB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m15.7/15.7 MB\u001b[0m \u001b[31m61.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m00:01\u001b[0m00:01\u001b[0m\n\u001b[?25hInstalling collected packages: numpy, scipy\nSuccessfully installed numpy-1.21.6 scipy-1.7.3\nNote: you may need to restart the kernel to use updated packages.\n","output_type":"stream"}]},{"cell_type":"code","source":"pip install ipympl","metadata":{"trusted":true},"execution_count":2,"outputs":[{"name":"stdout","text":"Collecting ipympl\n  Downloading ipympl-0.9.2-py2.py3-none-any.whl (510 kB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m510.3/510.3 kB\u001b[0m \u001b[31m10.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0ma \u001b[36m0:00:01\u001b[0m\n\u001b[?25hRequirement already satisfied: numpy in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipympl) (1.21.6)\nCollecting pillow\n  Downloading Pillow-9.2.0-cp37-cp37m-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (3.1 MB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m3.1/3.1 MB\u001b[0m \u001b[31m50.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0ma \u001b[36m0:00:01\u001b[0m\n\u001b[?25hRequirement already satisfied: ipython-genutils in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipympl) (0.2.0)\nRequirement already satisfied: traitlets<6 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipympl) (5.3.0)\nCollecting matplotlib<4,>=3.4.0\n  Downloading matplotlib-3.5.3-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.whl (11.2 MB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m11.2/11.2 MB\u001b[0m \u001b[31m60.0 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m00:01\u001b[0m00:01\u001b[0m\n\u001b[?25hRequirement already satisfied: ipython<9 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipympl) (7.33.0)\nRequirement already satisfied: ipywidgets<9,>=7.6.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipympl) (7.7.0)\nRequirement already satisfied: setuptools>=18.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (63.3.0)\nRequirement already satisfied: matplotlib-inline in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (0.1.3)\nRequirement already satisfied: decorator in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (5.1.1)\nRequirement already satisfied: pickleshare in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (0.7.5)\nRequirement already satisfied: pexpect>4.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (4.8.0)\nRequirement already satisfied: prompt-toolkit!=3.0.0,!=3.0.1,<3.1.0,>=2.0.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (3.0.30)\nRequirement already satisfied: jedi>=0.16 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (0.18.1)\nRequirement already satisfied: backcall in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (0.2.0)\nRequirement already satisfied: pygments in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipython<9->ipympl) (2.12.0)\nRequirement already satisfied: nbformat>=4.2.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipywidgets<9,>=7.6.0->ipympl) (5.4.0)\nRequirement already satisfied: widgetsnbextension~=3.6.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipywidgets<9,>=7.6.0->ipympl) (3.6.1)\nRequirement already satisfied: jupyterlab-widgets>=1.0.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipywidgets<9,>=7.6.0->ipympl) (1.1.1)\nRequirement already satisfied: ipykernel>=4.5.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipywidgets<9,>=7.6.0->ipympl) (6.15.1)\nCollecting cycler>=0.10\n  Downloading cycler-0.11.0-py3-none-any.whl (6.4 kB)\nRequirement already satisfied: packaging>=20.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib<4,>=3.4.0->ipympl) (21.3)\nCollecting fonttools>=4.22.0\n  Downloading fonttools-4.37.4-py3-none-any.whl (960 kB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m960.8/960.8 kB\u001b[0m \u001b[31m57.3 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n\u001b[?25hCollecting kiwisolver>=1.0.1\n  Downloading kiwisolver-1.4.4-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.whl (1.1 MB)\n\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m1.1/1.1 MB\u001b[0m \u001b[31m65.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n\u001b[?25hRequirement already satisfied: python-dateutil>=2.7 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib<4,>=3.4.0->ipympl) (2.8.2)\nRequirement already satisfied: pyparsing>=2.2.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from matplotlib<4,>=3.4.0->ipympl) (3.0.9)\nRequirement already satisfied: jupyter-client>=6.1.12 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (7.3.4)\nRequirement already satisfied: nest-asyncio in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (1.5.5)\nRequirement already satisfied: psutil in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (5.9.1)\nRequirement already satisfied: debugpy>=1.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (1.6.0)\nRequirement already satisfied: tornado>=6.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (6.2)\nRequirement already satisfied: pyzmq>=17 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (23.2.0)\nRequirement already satisfied: parso<0.9.0,>=0.8.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jedi>=0.16->ipython<9->ipympl) (0.8.3)\nRequirement already satisfied: typing-extensions in /srv/conda/envs/notebook/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib<4,>=3.4.0->ipympl) (4.3.0)\nRequirement already satisfied: jupyter-core in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (4.11.1)\nRequirement already satisfied: jsonschema>=2.6 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (4.9.0)\nRequirement already satisfied: fastjsonschema in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (2.16.1)\nRequirement already satisfied: ptyprocess>=0.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from pexpect>4.3->ipython<9->ipympl) (0.7.0)\nRequirement already satisfied: wcwidth in /srv/conda/envs/notebook/lib/python3.7/site-packages (from prompt-toolkit!=3.0.0,!=3.0.1,<3.1.0,>=2.0.0->ipython<9->ipympl) (0.2.5)\nRequirement already satisfied: six>=1.5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from python-dateutil>=2.7->matplotlib<4,>=3.4.0->ipympl) (1.16.0)\nRequirement already satisfied: notebook>=4.4.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (6.4.10)\nRequirement already satisfied: importlib-resources>=1.4.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jsonschema>=2.6->nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (5.9.0)\nRequirement already satisfied: attrs>=17.4.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jsonschema>=2.6->nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (22.1.0)\nRequirement already satisfied: pyrsistent!=0.17.0,!=0.17.1,!=0.17.2,>=0.14.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jsonschema>=2.6->nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (0.18.1)\nRequirement already satisfied: importlib-metadata in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jsonschema>=2.6->nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (4.11.4)\nRequirement already satisfied: pkgutil-resolve-name>=1.3.10 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jsonschema>=2.6->nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (1.3.10)\nRequirement already satisfied: entrypoints in /srv/conda/envs/notebook/lib/python3.7/site-packages (from jupyter-client>=6.1.12->ipykernel>=4.5.1->ipywidgets<9,>=7.6.0->ipympl) (0.4)\nRequirement already satisfied: terminado>=0.8.3 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.15.0)\nRequirement already satisfied: prometheus-client in /srv/conda/envs/notebook/lib/python3.7/site-packages (from notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.14.1)\nRequirement already satisfied: nbconvert>=5 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (6.5.0)\nRequirement already satisfied: argon2-cffi in /srv/conda/envs/notebook/lib/python3.7/site-packages (from notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (21.3.0)\nRequirement already satisfied: jinja2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (3.1.2)\nRequirement already satisfied: Send2Trash>=1.8.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (1.8.0)\nRequirement already satisfied: zipp>=3.1.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from importlib-resources>=1.4.0->jsonschema>=2.6->nbformat>=4.2.0->ipywidgets<9,>=7.6.0->ipympl) (3.8.0)\nRequirement already satisfied: beautifulsoup4 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (4.11.1)\nRequirement already satisfied: nbclient>=0.5.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.6.6)\nRequirement already satisfied: bleach in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (5.0.1)\nRequirement already satisfied: tinycss2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (1.1.1)\nRequirement already satisfied: jupyterlab-pygments in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.2.2)\nRequirement already satisfied: defusedxml in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.7.1)\nRequirement already satisfied: mistune<2,>=0.8.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.8.4)\nRequirement already satisfied: MarkupSafe>=2.0 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (2.1.1)\nRequirement already satisfied: pandocfilters>=1.4.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (1.5.0)\nRequirement already satisfied: argon2-cffi-bindings in /srv/conda/envs/notebook/lib/python3.7/site-packages (from argon2-cffi->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (21.2.0)\nRequirement already satisfied: cffi>=1.0.1 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from argon2-cffi-bindings->argon2-cffi->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (1.15.1)\nRequirement already satisfied: soupsieve>1.2 in /srv/conda/envs/notebook/lib/python3.7/site-packages (from beautifulsoup4->nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (2.3.2.post1)\nRequirement already satisfied: webencodings in /srv/conda/envs/notebook/lib/python3.7/site-packages (from bleach->nbconvert>=5->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (0.5.1)\nRequirement already satisfied: pycparser in /srv/conda/envs/notebook/lib/python3.7/site-packages (from cffi>=1.0.1->argon2-cffi-bindings->argon2-cffi->notebook>=4.4.1->widgetsnbextension~=3.6.0->ipywidgets<9,>=7.6.0->ipympl) (2.21)\nInstalling collected packages: pillow, kiwisolver, fonttools, cycler, matplotlib, ipympl\nSuccessfully installed cycler-0.11.0 fonttools-4.37.4 ipympl-0.9.2 kiwisolver-1.4.4 matplotlib-3.5.3 pillow-9.2.0\nNote: you may need to restart the kernel to use updated packages.\n","output_type":"stream"}]},{"cell_type":"code","source":"import numpy as np\nimport matplotlib.pyplot as plt\nimport scipy.linalg as scp\n\n# Input parameters\n# Nacelle\nMass_Nacelle = 5000             # [kg]\n\n# Pile\nm_Pile = 1000                   # [kg/m]\nEI_Pile = 1e9                   # [N.m2]\nEA_Pile = 1e8                   # [N]\nH_Pile = 100                    # [m]\n\n# Jacket\nm_Jacket = 100                  # [kg/m]\nEI_Jacket = 1e7                 # [N.m2]\nEA_Jacket = 1e6                 # [N]\nH_Jacket = 70                   # [m]\nBase_Jacket = 30.23             # [m]\nTop_Jacket = 16.13              # [m]\n\n# Define load parameters\nf0 = 2                          # [Hz]\nA0 = 0.1                        # [m]\nT0 = 20                         # [s]\n\n# Define output time vector\ndt = 0.01                       # [s]\nT = np.arange(0, 5*T0, dt)\nnT = len(T)","metadata":{"trusted":true},"execution_count":3,"outputs":[]},{"cell_type":"markdown","source":"## Step 1: discretize the domain \n\nAs done in tutorial 6, we start by discretizing the domain. We will discretize the structure using 21 nodes, as depicted in the following figure.\n\n![figure](tutorial8_figure2.png)\n\nHere we first define a matrix with the nodal coordinates.\n","metadata":{}},{"cell_type":"code","source":"# Define node coordinates\nNodeCoord = ([-(Base_Jacket*4+Top_Jacket*0)/4/2,      -4*H_Jacket/4],  # Point 1\n              [(Base_Jacket*4+Top_Jacket*0)/4/2,      -4*H_Jacket/4],  # Point 2\n             [-(Base_Jacket*3+Top_Jacket*1)/4/2,      -3*H_Jacket/4],  # Point 3\n              [(Base_Jacket*3+Top_Jacket*1)/4/2,      -3*H_Jacket/4],  # Point 4\n             [-(Base_Jacket*2+Top_Jacket*2)/4/2,      -2*H_Jacket/4],  # Point 5\n              [(Base_Jacket*2+Top_Jacket*2)/4/2,      -2*H_Jacket/4],  # Point 6\n             [-(Base_Jacket*1+Top_Jacket*3)/4/2,      -1*H_Jacket/4],  # Point 7\n              [(Base_Jacket*1+Top_Jacket*3)/4/2,      -1*H_Jacket/4],  # Point 8\n             [-(Base_Jacket*0+Top_Jacket*4)/4/2,       0*H_Jacket/4],  # Point 9\n                                             [0,       0*H_Jacket/4],  # Point 10\n              [(Base_Jacket*0+Top_Jacket*4)/4/2,       0*H_Jacket/4],  # Point 11\n                                             [0,                 10],  # Point 12\n                                             [0,                 20],  # Point 13\n                                             [0,                 30],  # Point 14\n                                             [0,                 40],  # Point 15\n                                             [0,                 50],  # Point 16\n                                             [0,                 60],  # Point 17\n                                             [0,                 70],  # Point 18\n                                             [0,                 80],  # Point 19\n                                             [0,                 90],  # Point 20\n                                             [0,                100])  # Point 21\nnNode = len(NodeCoord)","metadata":{"trusted":true},"execution_count":4,"outputs":[]},{"cell_type":"markdown","source":"Once we have the node coordinates we proceed to define the elemental connectivities. Here, we will use the same array to assign the material properties to each element. Note that they are different depending on which part of the structure they belong to.","metadata":{}},{"cell_type":"code","source":"# Define elements (and their properties\n#             NodeLeft    NodeRight          m         EA        EI\nElements = (         [1,            3,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [1,            4,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [2,            4,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [3,            4,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [3,            5,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [5,            4,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [4,            6,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [5,            6,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [5,            7,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [5,            8,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [6,            8,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [7,            8,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [7,            9,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [9,            8,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [8,           11,   m_Jacket,  EA_Jacket, EI_Jacket],\n                     [9,           10,   m_Jacket,  EA_Jacket, EI_Jacket],\n                    [10,           11,   m_Jacket,  EA_Jacket, EI_Jacket],\n                    [10,           12,   m_Pile,    EA_Pile,   EI_Pile],\n                    [12,           13,   m_Pile,    EA_Pile,   EI_Pile],\n                    [13,           14,   m_Pile,    EA_Pile,   EI_Pile],\n                    [14,           15,   m_Pile,    EA_Pile,   EI_Pile],\n                    [15,           16,   m_Pile,    EA_Pile,   EI_Pile],\n                    [16,           17,   m_Pile,    EA_Pile,   EI_Pile],\n                    [17,           18,   m_Pile,    EA_Pile,   EI_Pile],\n                    [18,           19,   m_Pile,    EA_Pile,   EI_Pile],\n                    [19,           20,   m_Pile,    EA_Pile,   EI_Pile],\n                    [20,           21,   m_Pile,    EA_Pile,   EI_Pile])\nnElem = len(Elements)","metadata":{"trusted":true},"execution_count":5,"outputs":[]},{"cell_type":"markdown","source":"Let's plot the structure to make sure that it looks like the model in the figure.","metadata":{}},{"cell_type":"code","source":"plt.figure()\nfor iElem in np.arange(0, nElem):\n    NodeLeft = Elements[iElem][0]-1\n    NodeRight = Elements[iElem][1]-1\n    m = Elements[iElem][2]\n    EA = Elements[iElem][3]\n    EI = Elements[iElem][4]\n    if m == m_Jacket and EA == EA_Jacket and EI == EI_Jacket:\n        c = 'b'\n    elif m == m_Pile and EA == EA_Pile and EI == EI_Pile:\n        c = 'r'\n    else:\n        print(\"ERROR: unknown material. Check your inputs.\")\n        break\n    plt.plot([NodeCoord[NodeLeft][0], NodeCoord[NodeRight][0]], [NodeCoord[NodeLeft][1], NodeCoord[NodeRight][1]], c=c)\nplt.axis('equal')\nplt.grid()","metadata":{"trusted":true},"execution_count":6,"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbbElEQVR4nO3deZSU1Z3/8feXNQrkgIJlgwt6JO4TIj38TGJiE01UxhPUEZUTl1EnrVF/JzHORIhLHI1Z3J1xz4SoGRUdxagJiWMIhRmPqI3iLpkWYcIiRKUJjcjSfH9/3KofZdtLLU/1033r8zqnT1U9z1PPvfdU8vFyn+fex9wdERGJU7+0KyAiItWjkBcRiZhCXkQkYgp5EZGIKeRFRCI2IO0KFBo5cqSPHTs27WokYsOGDQwZMiTtavSIWmmr2hmXmNq5cOHC99x9VEf7elXIjx07lqamprSrkYhsNktDQ0Pa1egRtdJWtTMuMbXTzJZ1tk/DNSIiEVPIi4hETCEvIhIxhbyISMQU8iIiESsp5M1sppmtMbPXCrbtZGZPmdn/5F5H5Labmf2rmTWb2StmdkjSlRcRka6V2pO/Gzi63bbpwFx3HwfMzX0GOAYYl/trBG4vv5oiIlKOkkLe3Z8GPmi3eQpwT+79PcBxBdvv9WABMNzM6iqoq0h6vvMd9rnllrRrIVKyJCZDZdx9Ve79u0Am934M8OeC45bntq0q2IaZNRJ6+mQyGbLZbAJVSl9ra2s0belOLbR1fDbLDm1t0bcTauP3hNppZ6IzXt3dzaykp5C4+13AXQD19fUeywy0mGbTdacm2jp8OC0tLfG3kxr5PamddiZxd83q/DBM7nVNbvsKYPeC43bLbRMRkR6SRMg/DpyRe38G8FjB9tNzd9kcCqwrGNYREZEeUNJwjZk9ADQAI81sOfAD4CfAQ2Z2NrAMOCl3+BxgMtAMfAicmVCdRUSkSCWFvLtP62TXER0c68D55VRKRESSoRmvIiIRU8iLiERMIS8iEjGFvIhIxBTyIiIR61XPeBXptZqb2WHTprRrIVIy9eRFRCKmnrxIMfbZh40tLQxOux4iJVJPXkQkYgp5EZGIKeRFRCKmkBcRiZhCXkQkYgp5EZGIKeRFRCKmkBcRiZgmQ4kUQ8saSB+lnryISMTUkxcphpY1kD5KPXkRkYhV3JM3s32BBws27Q1cDgwHvgn8Jbf9++4+p9LyRESkeBWHvLsvBsYDmFl/YAXwKHAmcKO7X1dpGSIiUp6kh2uOAN5292UJn1dERMqQ9IXXU4AHCj5fYGanA03ARe6+tv0XzKwRaATIZDJks9mEq5SO1tbWaNrSnVpo6/iWFtra2qJvJ9TG7wm1005z92ROZDYIWAkc6O6rzSwDvAc4cBVQ5+5ndXWO+vp6b2pqSqQ+actmszQ0NKRdjR5RE21taKClpYXhixalXZOqq4nfk7jaaWYL3b2+o31JDtccA7zo7qsB3H21u7e5+zbgZ8DEBMsSEZEiJBny0ygYqjGzuoJ9xwOvJViWiIgUIZExeTMbAnwVOKdg8zVmNp4wXLO03T4REekBiYS8u28Adm637bQkzi3SK2jtGumjNONVRCRiWrtGpBhau0b6KPXkRUQippAXEYmYQl5EJGIKeRGRiCnkRUQippAXEYmYQl5EJGIKeRGRiCnkRUQiphmvIsXQ2jXSR6knLyISMfXkRYqhtWukj1JPXkQkYgp5EZGIKeRFRCKmkBcRiZhCXkQkYgp5EZGIJXYLpZktBdYDbcBWd683s52AB4GxwFLgJHdfm1SZIiLStaR78pPcfby71+c+Twfmuvs4YG7us4iI9JBqD9dMAe7Jvb8HOK7K5YlUR3MzO6xYkXYtREpm7p7MiczeAdYCDtzp7neZWYu7D8/tN2Bt/nPB9xqBRoBMJjNh1qxZidQnba2trQwdOjTtavSIWmjr56dOZZs7zz38cNpVqbpa+D0hrnZOmjRpYcEIysckGfJj3H2Fme0CPAX8X+DxwlA3s7XuPqKzc9TX13tTU1Mi9UlbNpuloaEh7Wr0iJpoa0MDLS0tDF+0KO2aVF1N/J7E1U4z6zTkExuucfcVudc1wKPARGC1mdXlKlEHrEmqPBER6V4iIW9mQ8xsWP498DXgNeBx4IzcYWcAjyVRnoiIFCepWygzwKNh2J0BwP3u/jszewF4yMzOBpYBJyVUnoiIFCGRkHf3JcBnO9j+PnBEEmWIiEjpNONVRCRiCnkRkYgp5EVEIqaQFxGJmEJeRCRiepC3SDGam9lh06a0ayFSMvXkRUQipp68SDH22YeNLS0MTrseIiVST15EJGIKeRGRiCnkRUQippAXEYmYQl5EJGIKeRGRiCnkRUQippAXEYmYQl5EJGKa8SpSDK1dI32UevIiIhFTT16kGFq7RvqoinvyZra7mc0zszfM7HUz+3Zu+xVmtsLMFuX+JldeXRERKUUSPfmtwEXu/qKZDQMWmtlTuX03uvt1CZQhIiJlqDjk3X0VsCr3fr2ZvQmMqfS8IiJSuUTH5M1sLPA54Dngi8AFZnY60ETo7a/t4DuNQCNAJpMhm80mWaXUtLa2RtOW7tRCW8e3tNDW1hZ9O6E2fk+onXaauydzIrOhwHzganefbWYZ4D3AgauAOnc/q6tz1NfXe1NTUyL1SVs2m6WhoSHtavSImmhrQwMtLS0MX7Qo7ZpUXU38nsTVTjNb6O71He1L5BZKMxsIPALc5+6zAdx9tbu3ufs24GfAxCTKEhGR4iVxd40BPwfedPcbCrbXFRx2PPBapWWJiEhpkhiT/yJwGvCqmS3Kbfs+MM3MxhOGa5YC5yRQloiIlCCJu2v+G7AOds2p9NwivYaWNZA+SssaiIhETMsaiBRDyxpIH6WevIhIxBTyIiIRU8iLiERMIS8iEjGFvIhIxBTyIiIRU8iLiERMIS8iEjGFvIhIxDTjVaQYWrtG+ij15EVEIqaevEgxtHaN9FHqyYuIREwhLyISMYW8iEjEFPIiIhFTyIuIREwhLyISsaqHvJkdbWaLzazZzKZXuzwREdmuqiFvZv2BW4FjgAOAaWZ2QDXLFBGR7ardk58INLv7EnffDMwCplS5TJGi/e//ws03QxIrFixZAi+8UPl5RJJU7RmvY4A/F3xeDvyfKpcpUrRJk0I433wz7LFH58f9+FnYtu3TXNLQ+TFvvAFbtsDatYlXU6RsqS9rYGaNQCNAJpMhm82mW6GEtLa2RtOW7vTltg4ceAgwjHfeMTZv/pCRIzd3eNwr/Q7G+zktLS2dnsv9U6xfP4h5857GrDr17Ql9+fcsRa20E3ev2h/weeDJgs8zgBmdHT9hwgSPxbx589KuQo/py209+2z3XXd1P+44d3C//vrOj+2unddeG86xbl2ydexpffn3LEVM7QSavJNcrfaY/AvAODPby8wGAacAj1e5TJGirVwJo0fDQw/B1Klw0UXwk5+Ud666uu3nFOktqjpc4+5bzewC4EmgPzDT3V+vZpkipVi1CnbbDQYOhPvvhwEDYMaMMLZ+2WWlnWv06O3n3G+/5OsqUo6qj8m7+xxgTrXLESnHypUwcWJ4P2AA/PKX4fXyy2HzZrjySooeX8+HvHry0pukfuFVJC1btsCaNdvDGaB/f/jFL2DQIPjhD8MxP/5xcUGvkJfeSCEvNevdd8NrYchDCPq77gpDOD/9aQj6667r/nzDhsHQoQp56V0U8lKzVq0Kr/kLpoX69YPbbgs9+htuCEM3J5zQ/Tnr6rafV6Q3UMhLzcr3uNv35PPM4KabQo/++uth2bLPcPjh4T8AnRk9Wj156V20CqXUrO5CHkLQX3stTJ8OTzwxmm9+E9raOj9eIS+9jXryUrNWrgzj76NGdX2cGfzoR7By5VJmzhzLli3h4mz//p88Nh/y7sXflSNSTQp5qVkrV0Im03FYt2cGZ565lHHjxnLZZeFibP52y0J1dbBxI6xbB8OHV6XaIiVRyEvNWrWq66Gajlx6abgYe/HFIegfeCCM2ecVTohSyEtvoDF5qVn5JQ1K9b3vhTtuHnkkLIVQuEyx7pWX3kYhLzWr3JAHuPBCuOUWeOwx+Pu/h48+CtsV8tLbKOSlJm3eDO+9V37IA5x/Ptx5J/zmNzBlShiL1yJl0tso5KUm5We7djQRqhSNjTBzJjz1FBx7bLhAO2yYJkRJ76ELr1KTirlHvlhnnhkuvp5xBkyeDLvuqp689B4KealJSYY8wKmnhqD/xjdgyBBYtiyZ84pUSiEvNWnGjPDa2Ag77ljcd1paxnd7W+RnPgNvvgnPPw/vvw8771xRNUUqpjF5qUnhaZSwww7JnXPLlo+PxQ8blty5RcqlnrzUpC98IdwNM39+8d/JZhfR0NDQ4b7XXgt32Hz4IUybFiZJtbbCTjslU1+RcqknLzWpknvk23vsMfj850PAz58fwj5fhkjaFPJSk1atqvz2Sffw9KjjjoP994emJjj00I8vbSCSNoW81KRKe/IbNsDJJ4eHfZ96aujBjxkT9mnWq/QmFYW8mV1rZm+Z2Stm9qiZDc9tH2tmG81sUe7vjkRqK5KAjz6CDz4oP+SXLYPDDgtr11x7Ldx778cv4GrWq/QmlV54fQqY4e5bzeynwAzg4ty+t919fIXnF0lcfhilnJD/4x/DWjWbN8Ovfw3HHPPJY3bcMaxAqZCX3qCinry7/5e7b819XADsVnmVRKqr3JB/4ok6vvIVGDECnnuu44DP07NepbdIckz+LOC3BZ/3MrOXzGy+mX0pwXJEKlLqbNctW8JiZDfcsC9HHhkCft99u/6OHgMovUW3wzVm9ntg1w52XeLuj+WOuQTYCtyX27cK2MPd3zezCcCvzOxAd/9rB+dvBBoBMpkM2Wy2rIb0Nq2trdG0pTt9ra3z548BxrFkyTN88MGWLo9dt24gV1xxAIsWjeCEE97mvPP+zKJF3ZfRr99+LFkynGx2QSJ17kl97fcsV620E3ev6A/4B+BZYMcujskC9d2da8KECR6LefPmpV2FHtPX2nrxxe4DB7q3tXV93Msvu48d6z54sPsvf1laO/NlbNtWWV3T0Nd+z3LF1E6gyTvJ1Urvrjka+B7wdXf/sGD7KDPrn3u/NzAOWFJJWSJJWbkyjJn36+J//bNnh1mxmzfD00+H2yRLUVcXhnnef7+yuopUqtIx+VuAYcBT7W6V/DLwipktAh4GznX3DyosSyQRXT3bdds2+Jd/CXfQHHRQmOA0cWLpZWhClPQWFd1C6e77dLL9EeCRSs4tUi0rV8J++31ye2trWBN+9uzwescd8KlPlVdG4YSogw8uv64ildKMV6k5Hc12XboUvvhF+NWvwkO6f/GL8gMeNOtVeg+tQik1ZeNGaGn5+Lo12SyceCK0tcFvfwtf+1rl5WjWq/QW6slLTWk/Eer22+GrX4VRo8KDPpIIeAj/ChgxQmPykj6FvNSUfM961Cg491w47zw46ihYsADGjUu2LE2Ikt5AwzVSU/Khe/nl8OKLMH16WC64f//ky1LIS2+gkJea8k//FF5feimsAf/ss3DEEcV9t5hnvBZ6+eXwIBGRNGm4RmrK6tXhdfBgGDq0umW1tYXJVG1t1S1HpCvqyUtNmTo13EFjFp7Leu21YVzerPvvdvWM147cfDN85zth7fpRo8quskhF1JOXmrJqVZgI9eqrcPjhcMEF8Hd/B+++m3xZu+UW3ta4vKRJIS81Jf9s17o6mDMHbrkF5s0LSxg8+miyZeleeekNFPJSUwpnu5qFdeJffBH23BNOOAH+8R9h/fpkytL6NdIbKOSlZmzYAOvWfXJJg/xdNjNmwMyZMH58+Fwp9eSlN1DIS83o6rF/gwbBj34E8+eHu2EOOyzcS7+l62eKdGnwYNh5Z4W8pEshLzUjH/KF69a096UvwSuvwGmnwVVXhUXL/vSn8susq1PIS7oU8lIzin2266c/DXffDQ89BM3N8LnPwZ13QnjIWWk061XSppCXmlHqA7ynTg23Wn7hC2Gdm0suOej/T6Yq1ujRuvAq6VLIS81YuTKsDlnK0gRjxsCTT8JNN0FT004cfDA88UTx38+H/LZtpdZWJBkKeakZ+We7FjO7tVC/fvDtb8Oddy5k9Gj4+tdDz37Dhu6/W1cXLuT+5S/l1VmkUgp5qRldPdu1GHvttYHnnoN//me4664wVv/8811/R0+IkrQp5KVmdPTYv1INHgzXXAN/+AN89FEYr7/ySti6tePjFfKSNoW81IwkQj6voSHcannyyfCDH4RbL99++5PHKeQlbRWFvJldYWYrzGxR7m9ywb4ZZtZsZovN7KjKqypSvtbWsFxBUiEP4QLufffB/ffDm2/CZz8LP//5x2+13HXX8Ko7bCQtSfTkb3T38bm/OQBmdgBwCnAgcDRwm5lV4dk7IsUpZiJUuaZNC736v/3bsPbNCSdsv9A6aBCMHKmevKSnWsM1U4BZ7r7J3d8BmoGJVSpLpFul3iNfqj32gLlz4brrwuqWBx8c1q3Pl6mQl7Qk8dCQC8zsdKAJuMjd1wJjgAUFxyzPbfsEM2sEGgEymQzZbDaBKqWvtbU1mrZ0py+0de7cXYADWL78ebLZ8p7JV0w7J0yAW28dwtVXH8DkyUOYMmUFgwbtwOLFA8lmF5ZVbk/rC79nEmqlnbh7l3/A74HXOvibAmSA/oR/EVwNzMx95xbg1IJz/Bw4sbuyJkyY4LGYN29e2lXoMX2hrZde6g7ua9eWf45S2rlxo/uFF4Yyhw1zHzq0/HJ7Wl/4PZMQUzuBJu8kV7vtybv7kcX8x8LMfgb8OvdxBbB7we7dcttEUnHnneF1ypTSJ0Pllfogb4C994YlS8L7piaory+vbJFyVXp3TeFlrOMJPXyAx4FTzGywme0FjAO6mTYiUj077ggDB5Yf8KXasgXeeisEfL/c/8symZ4pW6RQpWPy15jZeMCBpcA5AO7+upk9BLwBbAXOd3c9s15Ss+ee4eJoJUOwxTzIe9s2+Pd/h+nTw22b3/8+/M3fwCmnwJo1sPvuXX5dJHEVhby7n9bFvqsJ4/QiqVu1Cg45pLplvPwyfOtb4alShx8Ot98enjqVX/pA98pLGjTjVaLnnuxs1/bWr4fvfjfcWdPcDPfeGx4Ovv/+Yb9mvUqakriFUqRXW78+rBiZdMi7w+zZYYXKFSvgnHPCIwR32unjx2Uy4VqAQl7SoJ68RK8aE6GWLIFjj4UTTwwzWp99Fu6445MBD+GC7y67KOQlHQp5iV6SSxps2gRXXw0HHghPPw033hhujTz00K6/V1enMXlJh4ZrJHpJ9eRfemk43/pWuDXyxBNDwO+2W3Hf1dIGkhb15CV6lYb86tVw2mnw3e+OZ9OmsDbNf/5n8QGfL1shL2lQyEv0Vq6EoUNh2LDSvrdtW5gpu99+8OCDcNppS3n9dTjmmNLrMHp0+I9FZw8XEakWhbxEL/9s11K89FJ46tO554bH/L3yCpx11lJ22KG8OtTVhbtxVq8u7/si5VLIS/RKebbr+vVw4YVhjZl33oH/+I+whPB++1VWh3z5uvgqPU0hL9ErZiKUOzz8cAjzm28O97y/9RZ84xvJrHejCVGSFoW8RK2Y2a5vvw2TJ8PUqWHi0oIFcNttMGJEcvVQyEtaFPIStXXrYOPGjkN+0yb44Q/hoIPgmWfgppvCOjMTq/AMs112CatRKuSlp+k+eYlaPlTbX3j9wx/gvPNg8WI46SS44QYY0+Gzy5IxYIBmvUo61JOXqOUvdOZ78qtXw6mnwhFHhNsZf/e7cHtkNQM+b/RoXXiVnqeQl6jle86ZTFj6d999w0Smyy6DV1+Fo47qubpoQpSkQcM1ErV8qJ56KixcGHrwt94awr6njR69fW15kZ6ikJeozZ4dXpcvh/vug2nTeu4RgO3V1YWnQ23ZElamFOkJCnmJ2v77h8fwPfMMJT+EO2n19eE/Mhs3KuSl5yjkJWp33512DbY79tjwJ9KTdOFVRCRiFfXkzexBIH8JazjQ4u7jzWws8CawOLdvgbufW0lZIiJSuopC3t1Pzr83s+uBdQW733b38ZWcX0REKpPImLyZGXAS8JUkziciIskwd6/8JGZfBm5w9/rc57HA68CfgL8Cl7r7Hzv5biPQCJDJZCbMmjWr4vr0Bq2trQwdOjTtavSIWmmr2hmXmNo5adKkhfn8ba/bkDez3wO7drDrEnd/LHfM7UCzu1+f+zwYGOru75vZBOBXwIHu/teuyqqvr/empqbu2tMnZLNZGhoa0q5Gj6iVtqqdcYmpnWbWach3O1zj7kd2c/IBwAnAhILvbAI25d4vNLO3gc8AcSS4iEgfkcQtlEcCb7n78vwGMxtlZv1z7/cGxgFLEihLRERKkMSF11OAB9pt+zJwpZltAbYB57r7B92daOHChe+Z2bIE6tQbjATeS7sSPaRW2qp2xiWmdu7Z2Y5ELrzKJ5lZU2djZLGplbaqnXGplXZqxquISMQU8iIiEVPIV89daVegB9VKW9XOuNREOzUmLyISMfXkRUQippAXEYmYQj4BZjbVzF43s21mVt9u3wwzazazxWZ2VMH2o3Pbms1ses/XujJmdoWZrTCzRbm/yQX7OmxzX9XXf6uumNlSM3s19xs25bbtZGZPmdn/5F5HpF3PcpjZTDNbY2avFWzrsG0W/GvuN37FzA5Jr+bJUsgn4zXC0g5PF240swMIk8UOBI4GbjOz/rnZwLcCxwAHANNyx/Y1N7r7+NzfHOi8zWlWshIR/VZdmZT7DfMdlOnAXHcfB8zNfe6L7ib8b7BQZ207hjAzfxxhwcTbe6iOVaeQT4C7v+nuizvYNQWY5e6b3P0doBmYmPtrdvcl7r4ZmJU7Ngadtbmvivm36swU4J7c+3uA49KrSvnc/Wmg/Uz7zto2BbjXgwXAcDOr65GKVplCvrrGAH8u+Lw8t62z7X3NBbl/2s4s+Cd9LG3Li6097TnwX2a2MLfsN0DG3Vfl3r8LZNKpWlV01rZof2c9yLtIxSy5HJuu2kz45+xVhJC4CrgeOKvnaicJOczdV5jZLsBTZvZW4U53dzOL8j7rmNtWSCFfpO6WXO7ECmD3gs+75bbRxfZeo9g2m9nPgF/nPnbV5r4otvZ8jLuvyL2uMbNHCcNTq82szt1X5YYs1qRayWR11rZof2cN11TX48ApZjbYzPYiXNR5HngBGGdme5nZIMKFysdTrGfJ2o1XHk+4+Aydt7mv6vO/VWfMbIiZDcu/B75G+B0fB87IHXYGENO/VDtr2+PA6bm7bA4F1hUM6/Rp6sknwMyOB/4NGAX8xswWuftR7v66mT0EvAFsBc5397bcdy4AngT6AzPd/fWUql+ua8xsPGG4ZilwDkBXbe6L3H1rBL9VZzLAo+ERzQwA7nf335nZC8BDZnY2sIzw/OY+x8weABqAkWa2HPgB8BM6btscYDLhRoEPgTN7vMJVomUNREQipuEaEZGIKeRFRCKmkBcRiZhCXkQkYgp5EZGIKeRFRCKmkBcRidj/A+y/aSyO3QYLAAAAAElFTkSuQmCC\n"},"metadata":{"needs_background":"light"}}]},{"cell_type":"markdown","source":"## Step 2: define the shape functions\n\nHere we will use the exact same functions as in tutorial 4 and 5. Linear shape functions for the axial displacement and cubic shape functions for the deflection and rotations. Since we already know its expression and we already have the value of the elemental matrices, we skip this step in this tutorial. ","metadata":{}},{"cell_type":"markdown","source":"## Step 3: computation of the elemental matrices\n\nIn the theory we have seen that the mass and stiffness elemental matrices for the space frame using linear and cubic shape functions are given by:\n\n$$ \\bold{M} = \\frac{mL}{420} \\begin{bmatrix} 140 & 0 & 0 & 70 & 0 & 0 \\\\ 0 & 156 & 22L & 0 & 54 & -13L \\\\ 0 & 22L & 4L^2 & 0 & 13L & -3L^2 \\\\ 70 & 0 & 0 & 140 & 0 & 0 \\\\ 0 & 54 & 13L & 0 & 156 & -22L \\\\ 0 & -13L & -3L^2 & 0 & -22L & 4L^2 \\end{bmatrix} $$\n$$ \\quad \\bold{K} = \\begin{bmatrix} \\frac{EA}{L} & 0 & 0 & \\frac{-EA}{L} &  0 & 0 \\\\ 0 & \\frac{12EI}{L^3} & \\frac{6EI}{L^2} & 0 & \\frac{-12EI}{L^3} & \\frac{6EI}{L^2} \\\\ 0 & \\frac{6EI}{L^2} & \\frac{4EI}{L} & 0 & \\frac{-6EI}{L^2} & \\frac{2EI}{L} \\\\ \\frac{-EA}{L} & 0 & 0 & \\frac{EA}{L} & 0 & 0 \\\\ 0 & \\frac{-12EI}{L^3} & \\frac{-6EI}{L^2} & 0 & \\frac{12EI}{L^3} & \\frac{-6EI}{L^2} \\\\ 0 & \\frac{6EI}{L^2} & \\frac{2EI}{L} & 0 & \\frac{-6EI}{L^2} & \\frac{4EI}{L} \\end{bmatrix}$$\n\nThese matrices are used directly when calling the `BeamMatrices` function within the assembly process.","metadata":{}},{"cell_type":"markdown","source":"## Step 4: global assembly\n\nThe last step is to compute the global matrices and the global forcing vector. We start by initializing the global matrices as 1-dimensional arrays.","metadata":{}},{"cell_type":"code","source":"nDof = 3*nNode              # 3 Degrees of freedom per node\nK = np.zeros((nDof*nDof))\nM = np.zeros((nDof*nDof))","metadata":{"trusted":true},"execution_count":7,"outputs":[]},{"cell_type":"markdown","source":"Then we loop over elements and perform all the elemental operations.","metadata":{}},{"cell_type":"code","source":"from BeamMatrices import BeamMatricesJacket\nfor iElem in np.arange(0, nElem):\n    # Get the nodes of the elements\n    NodeLeft = Elements[iElem][0]-1\n    NodeRight = Elements[iElem][1]-1\n    \n    # Get the degrees of freedom that correspond to each node\n    Dofs_Left = 3*(NodeLeft) + np.arange(0, 3)\n    Dofs_Right = 3*(NodeRight) + np.arange(0, 3)\n\n    # Get the properties of the element\n    m = Elements[iElem][2]\n    EA = Elements[iElem][3]\n    EI = Elements[iElem][4]\n\n    # Calculate the matrices of the element\n    Me, Ke = BeamMatricesJacket(m, EA, EI, ([NodeCoord[NodeLeft][0], NodeCoord[NodeLeft][1]], [NodeCoord[NodeRight][0], NodeCoord[NodeRight][1]]))\n\n    # Assemble the matrices at the correct place\n    nodes = np.append(Dofs_Left, Dofs_Right)\n    for i in np.arange(0, 6):\n        for j in np.arange(0, 6):\n            ij = nodes[j] + nodes[i]*nDof\n            #print(ij)\n            M[ij] = M[ij] + Me[i, j]\n            K[ij] = K[ij] + Ke[i, j]\n            \n# Reshape the global matrix from a 1-dimensional array to a 2-dimensional array\nM = M.reshape((nDof, nDof))\nK = K.reshape((nDof, nDof))","metadata":{"trusted":true},"execution_count":8,"outputs":[]},{"cell_type":"markdown","source":"Now we have the global mass and stiffness matrices. However, in this example we have an additional point mass at the top corresponding to the nacelle. Then, we need to account for this mass adding its value at the corresponding DOFs, in this case the corresponding horizontal and vertical displacements associated to the top node.","metadata":{}},{"cell_type":"code","source":"nacelle_node = nNode\nnacelle_dof_h = 3*(nacelle_node-1)\nnacelle_dof_v = 3*(nacelle_node-1)+1\nM[nacelle_dof_h, nacelle_dof_h] += Mass_Nacelle\nM[nacelle_dof_v, nacelle_dof_v] += Mass_Nacelle","metadata":{"trusted":true},"execution_count":9,"outputs":[]},{"cell_type":"markdown","source":"That completes the filling of the matrices. Let's have a look at the matrices' structure.","metadata":{}},{"cell_type":"code","source":"# Look at the matrix structure\nplt.figure()\nplt.spy(M)\nplt.title(\"Mass matrix\")\nplt.figure()\nplt.spy(K)\nplt.title(\"Stiffness matrix\")\n","metadata":{"trusted":true},"execution_count":10,"outputs":[{"execution_count":10,"output_type":"execute_result","data":{"text/plain":"Text(0.5, 1.0, 'Stiffness matrix')"},"metadata":{}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAPsAAAEMCAYAAAACksCzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAASe0lEQVR4nO3de4xc5X3G8e8TzK3hYi7brcFQU+FCUVIMrLgUmiZ2odwa+w9EQmm7UFduE9oSJVIKVIogihSo2hCQKloXSNxAuMQJsUtSwHUgbdXEsA53DDUYI2xs7zrg2iBKMfz6x3kXhu1ezs6cmTm77/ORRjvnnJk5P3v22fdc3nNeRQRmNv19qNsFmFlnOOxmmXDYzTLhsJtlwmE3y4TDbpYJh91qQ9JVkm7udh3TlXyevV4kbQQOAw6LiO0N8x8F5gFHRcTGrhTXJEkfB26LiNldLiVrbtnr6UXgouEJSR8FfqF75bSfpBndrmG6c9jr6VvAHzZM9wP/1PgCSedJelTSTkkvS7q6Ydk+km6T9HNJOyQ9Iqk3LbtE0gZJuyS9KOni0QqQdLWk76TP2SXpSUm/KulKSYNpnWc1vP5SSevSazdI+pM0/8PAvwCHSXo9PQ5Ln788ff5O4JI077b0vk+l+g5I0+dI2iqpp4L/3yw57PX0U+AASb8maQ/g08BtI17zBsUfhJnAecBnJC1Ky/qBA4EjgEOAPwXeTMG7ETgnIvYHfgN4bJw6fpfiD89BwKPA/RS/M4cDXwb+oeG1g8D5wAHApcD1kk6MiDeAc4BXImK/9HglvWchsDz9G25vXHFE3AX8J3CjpEOAW4A/joihceq1cTjs9TXcup8JrAM2Ny6MiIci4smIeDcingDuAH4rLX6bIuRHR8Q7EbE2InamZe8CH5G0b0RsiYinx6nh3yPi/ojYDXwH6AGujYi3gTuBOZJmpnp+EBEvROHHwAPAb07wb/xJRHw//RveHGX5ZcB84CHgnyPi3gk+z8bhsNfXt4DfAy5hxCY8gKRTJD0oaUjSf1O03oc2vPd+4E5Jr0j6a0l7plb2U+m1WyT9QNKx49SwreH5m8D2iHinYRpgv1TPOZJ+KulVSTuAcxvqGcvL4y2MiB0Uf2Q+AvztBJ9lE3DYayoiXqI4UHcu8L1RXvJtYCVwREQcCPw9oPTetyPimog4jmJT/XzSMYDUUp8JzAKeBf6x1Vol7Q18F/gboDciZgI/HK4HGOuUz7ingiTNA/6IYqvlxlbrzJ3DXm+LgfmpRR5pf+DViPgfSSdTbAUAIOkTkj6a9vd3UmzWvyupV9LCtO/+FvA6xWZ9q/YC9gaGgN2SzgHOali+DThE0oFlP1DSPhTHKa6iOAZwuKTPVlBrthz2Gkv7wANjLP4s8GVJu4AvAXc3LPsligNfOyn2939MsWn/IeDzwCvAqxT7+J+poM5dwF+kGl6j+MOzsmH5sxSt84Z0duCwEh/7VeDliLgpIt4Cfh/4iqS5rdabK3eqMcuEW3azTDjsZplw2M0y4bCbZcJhN8tER8Mu6WxJz0l6XtIVHV73rekCjqca5h0saZWk9ennQR2o44jU8+0ZSU9LuryLtewj6WFJj6darknzj5K0Jn1Pd0naq921pPXukS7uubfLdWxMF/48JmkgzevG9zMzXSz0bLrI6LRW6uhY2FMHj7+juCjiOOAiScd1av3AN4GzR8y7AlgdEXOB1Wm63XYDX0i9204FLkv/D92o5S2KTjvHU1wrf7akU4HrgOsj4miK8+aLO1ALwOUU/QKGdasOgE9ExLyI6EvT3fh+bgDui4hjgeMp/m+aryMiOvIATgPub5i+EriyU+tP65wDPNUw/RwwKz2fBTzXyXrSeldQXOzS1Voorpf/GXAKsB2YMdr31sb1z06/vPOBeym62na8jrSujcChI+Z19PuhuGrxRVJfmCrq6ORm/OF88MKHTWleN/VGxJb0fCvQ28mVS5oDnACs6VYtadP5MYpLVFcBLwA7orjSDTr3PX0d+CLvd989pEt1QNFn/wFJayUtSfM6/f0cRdH9+Btp1+bm1M256Tp8gC6J4k9lx7oTStqP4uKRz8X7l592vJYoLoGdR9GyngyMdxVcW0g6HxiMiLWdXvcYzoiIEyl2OS+T9LHGhR36fmYAJwI3RcQJFPcv+MAm+2Tr6GTYN1PcTGHYbEZco90F2yTNAkg/BzuxUkl7UgT99ogYvqKtK7UMi+Jy0gcpNpdn6v3bRHXiezod+KSK++/dSbEpf0MX6gAgIjann4PAPRR/BDv9/WwCNkXEmjS9nCL8TdfRybA/AsxNR1j3orj7ysoJ3tNuKynu6kL6uaLdK5QkiruurIuIr3W5lp7hm09I2pf3b5TxIHBBp2qJiCsjYnZEzKH4vfhRRFzc6TqguI2WpP2Hn1NcvfcUHf5+ImIr8LKkY9KsBcAzLdXRiQMeDQcXzgX+i2K/8K86vO47gC0Ul3tuojiyewjFQaH1wL8CB3egjjMoNr2eoLgl1GPp/6Ubtfw6xe2mnqD4hf5Smv8rwMPA8xQ3j9i7g9/Tx4F7u1VHWufj6fH08O9pl76fecBA+n6+T3F7sKbr8FVvZpnwATqzTDjsZplw2M0y4bCbZcJhN8tEV8Le0AWxq+pSB9SnlrrUAfWppS51QGu1tBT2Fi5Zrct/Xl3qgPrUUpc6oD611KUOaKGWpsNeg0tWzWwSmu5UI+k04OqI+J00fSVARHx1nPeUWtlJJ51Uuo61a6u/dmIy66/K0NAQPT3dH6C0LnVAfWqpSx0wcS0bN25k+/btGm1ZK2Nij3bJ6iktfN57BgbGGhfh/yu6mldrMus3q5O+vr4xl7US9lLSAYU67fOYZamVsJe6ZDUilgJLofxmvJlVr5Wj8XW8ZNXMxtB0yx4RuyX9GcU44HsAt0bE05VVZmaVammfPSJ+SDEOt5nVnLvLmmXCYTfLREfDftJJJ5W7fY5U+jGJW/yYZc0tu1kmHHazTDjsZplw2M0y4bCbZcJhN8uEw26WCYfdLBMOu1kmHHazTLT9TjWN1q5dW+o2UpPp3jre5433OeMtK3urK3fDtanELbtZJhx2s0x0dDO+0xo3x0ducrfjrrRmdeaW3SwTDrtZJhx2s0xM+X32xn1x74ebjc0tu1kmHHazTNRyM77ZzfEqesVNhnva2VTilt0sExOGXdKtkgYlPdUw72BJqyStTz8Pam+ZZtaqMi37N4GzR8y7AlgdEXOB1WnazGpswrBHxL8Br46YvRBYlp4vAxZVW1Z5ZQeXMMtds/vsvRGxJT3fCvRWVI+ZtUnLR+MjIiSNebhZ0hJgSavrMbPWNNuyb5M0CyD9HBzrhRGxNCL6IqKvyXWZWQWaDftKoD897wdWVFOOmbVLmVNvdwA/AY6RtEnSYuBa4ExJ64HfTtNmVmMT7rNHxEVjLFpQcS1m1ka1HJ99Mo+6nF6rarz5sg+zyXJ3WbNMOOxmmajlVW9ThW+cYVOJW3azTDjsZplw2M0y4bCbZcJhN8uEw26WiSl/6q3szRw7fRPLZm8y6VN41i5u2c0y4bCbZWLKb8a3QxU948YbLtqsG9yym2XCYTfLhMNulgnvs7fAV73ZVOKW3SwTDrtZJrLZjJ/M6a9mNsnb0dOuiiGofdrPhrllN8uEw26WCYfdLBPZ7LPX1ch9b5/Os3Zxy26WiTJjvR0h6UFJz0h6WtLlaf7BklZJWp9+HtT+cs2sWWVa9t3AFyLiOOBU4DJJxwFXAKsjYi6wOk3bKEYOB2XWDROGPSK2RMTP0vNdwDrgcGAhsCy9bBmwqE01mlkFJrXPLmkOcAKwBuiNiC1p0Vagt9rSzKxKpcMuaT/gu8DnImJn47Iotk1H3T6VtETSgKSBoaGhloo1s+aVCrukPSmCfntEfC/N3iZpVlo+Cxgc7b0RsTQi+iKir6enp4qaa6PZ4ZbrOgS1h4ue3socjRdwC7AuIr7WsGgl0J+e9wMrqi/PzKpSplPN6cAfAE9KeizNuwq4Frhb0mLgJeDCtlRoZpWYMOwR8R/AWNtuC6otZ3oaebrNm8LWDe5BZ5YJh90sEw67WSYcdrNMOOxmmXDYzTLhm1eMouor07p5qm28m2OM5Cvypje37GaZcNjNMuGwm2XC++zT3GS66nqffXpzy26WCYfdLBMOu1kmHHazTDjsZpnw0fgOaMdR7saj6lUM7dyuz/QR/vpwy26WCYfdLBMOu1kmvM8+RbVjaGcPFz29uWU3y4TDbpYJh90sEw67WSbKjPW2j6SHJT0u6WlJ16T5R0laI+l5SXdJ2qv95ZpZs8q07G8B8yPieGAecLakU4HrgOsj4mjgNWBx26o0s5ZNGPYovJ4m90yPAOYDy9P8ZcCidhSYm2aGTR5vCOiRqhguerzPmGj91j1lx2ffI43gOgisAl4AdkTE7vSSTcDhbanQzCpRKuwR8U5EzANmAycDx5ZdgaQlkgYkDQwNDTVXpZm1bFJH4yNiB/AgcBowU9JwD7zZwOYx3rM0Ivoioq+np6eVWq0J422qj7drYNNPmaPxPZJmpuf7AmcC6yhCf0F6WT+wok01mlkFyvSNnwUsk7QHxR+HuyPiXknPAHdK+grwKHBLG+s0sxZNGPaIeAI4YZT5Gyj2381sCnAPOrNMOOxmmXDYzTLhm1fUTDeHiy677skMIVV2/e5t135u2c0y4bCbZcJhN8uEw26WCYfdLBMOu1kmfOrNJq2qceCss9yym2XCYTfLhMNulgnvs09z7R4bvlnuVtt5btnNMuGwm2XCYTfLhMNulgmH3SwTPhpvbTXW0XP3tOs8t+xmmXDYzTLhsJtlwvvsNmmT6cVWdW+7dtxAMxdu2c0yUTrsaYz2RyXdm6aPkrRG0vOS7pK0V/vKNLNWTaZlv5xi9NZh1wHXR8TRwGvA4ioLM7NqlQq7pNnAecDNaVrAfGB5eskyYFEb6jOzipRt2b8OfBF4N00fAuyIiN1pehNweLWlmVmVJgy7pPOBwYhY28wKJC2RNCBpYGhoqJmPMLMKlDn1djrwSUnnAvsABwA3ADMlzUit+2xg82hvjoilwFKAvr4+nwuxUfkmlu03YcseEVdGxOyImAN8GvhRRFwMPAhckF7WD6xoW5Vm1rJWzrP/JfB5Sc9T7MPfUk1JZtYOk+pBFxEPAQ+l5xuAk6svyaaTKoaBbnZdvq/dB7kHnVkmHHazTDjsZplw2M0y4bCbZcJhN8uEb15hU4p72jXPLbtZJhx2s0w47GaZ8D671UInu9WOXF8u3WrdsptlwmE3y4TDbpYJh90sEw67WSYcdrNM+NSbTVvuWvtBbtnNMuGwm2XCm/E2pXi46Oa5ZTfLhMNulgmH3SwTDrtZJkodoJO0EdgFvAPsjog+SQcDdwFzgI3AhRHxWnvKNLNWTaZl/0REzIuIvjR9BbA6IuYCq9O0mdVUK5vxC4Fl6fkyYFHL1Zh1QUSM+ZD03mOqKxv2AB6QtFbSkjSvNyK2pOdbgd7KqzOzypTtVHNGRGyW9IvAKknPNi6MiJA0ai+C9MdhCcCRRx7ZUrFm1rxSLXtEbE4/B4F7KIZq3iZpFkD6OTjGe5dGRF9E9PX09FRTtVkJ422ej7WpPvLR7LrG+8xu7RpMGHZJH5a0//Bz4CzgKWAl0J9e1g+saFeRZta6MpvxvcA96a/QDODbEXGfpEeAuyUtBl4CLmxfmWbWqgnDHhEbgONHmf9zYEE7ijKz6rkHnVkmHHazTDjsZpnwzSvMSprq97Rzy26WCYfdLBMOu1kmvM9u2ZsKw0WPfF8z3LKbZcJhN8uEw26WCYfdLBMOu1kmHHazTPjUm1kb1LFrrVt2s0w47GaZ8Ga8WUndHC66is90y26WCYfdLBMOu1kmHHazTDjsZplw2M0y4VNvZl3UyZ52btnNMlEq7JJmSlou6VlJ6ySdJulgSaskrU8/D2p3sWbWvLIt+w3AfRFxLMW4b+uAK4DVETEXWJ2mzaymygzZfCDwMeAWgIj434jYASwElqWXLQMWtadEs6mnk2PDl1WmZT8KGAK+IelRSTencdp7I2JLes1WiqGdzaymyoR9BnAicFNEnAC8wYhN9igOKY56WFHSEkkDkgaGhoZardfMmlQm7JuATRGxJk0vpwj/NkmzANLPwdHeHBFLI6IvIvp6enqqqNnMmjBh2CNiK/CypGPSrAXAM8BKoD/N6wdWtKVCM6tE2U41fw7cLmkvYANwKcUfirslLQZeAi5sT4lmVoVSYY+Ix4C+URYtqLQaM2sb96Azy4TDbpYJh90sE2p1GNhJrUwaojiYdyiwvWMrHltd6oD61FKXOqA+tdSlDpi4ll+OiFHPcXc07O+tVBqIiNEO+GVZB9SnlrrUAfWppS51QGu1eDPeLBMOu1kmuhX2pV1a70h1qQPqU0td6oD61FKXOqCFWrqyz25mnefNeLNMOOxmmXDYzTLhsJtlwmE3y8T/AXfNs8+l9LkCAAAAAElFTkSuQmCC\n"},"metadata":{"needs_background":"light"}},{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 1 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAPsAAAEMCAYAAAACksCzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAATGklEQVR4nO3dfbBcdX3H8fdHIkJBDMFrJhIQLAhFW4K5KgzUBxClVCV/KNWxNWjatNZ2cKRa0LZTR6diO4Mw01YnihCV8mDUJoOMEgPU2ofgjTwYCArSWJImuTdK5KFKG/z2j9/v4nJ7H87dPWf33Pv7vGZ27p5zdvd8d/d+9ncefuccRQRmNv89Y9AFmFl/OOxmhXDYzQrhsJsVwmE3K4TDblYIh72FJH1K0p93DL9b0h5Jj0k6QtLpku7PwysGWOpASDo6v/cDBl3LXCLvZ2+epDOAvwZeDDwJbAPeGxHflnQB8LsRccYUz30m8AhwakTclcdtAjZExBX9qL+fJG0nfR7fGHQt882CQRcw30k6DLgReDdwA3Ag8OvAExVfYjFwEHBPx7gXTBguhqQFEbF/0HXMSRHhW4M3YBjYN8W0XwF+RmrtHxt/HHA18FHgRcDjQOTptwA/AH4O/DSPexZwG/AR4F+AR4Gbged2zOdU4F+BfcBdwKs7pl0APJif9x/A2/P444B/An4C7AWun+I9HJPreyfwEPAw8AfAy4C78zz/tuPxv5zfx4/y614DLMzTPj/hvX2g4/VXAf8JfLNj3AJgEbADeGN+jUOBB4B3DPq7b9tt4AXM9xtwWP7HXgv8BnD4hOkXAN+aMO5q4KP5/lP/2B3TtwOv7Ri+Lf8IvAg4OA9fmqcdmed/Lmkbzdl5eAg4hLSKcEJ+7BLgxfn+tcCH8nMOAs6Y4v2N1/ep/LjXkX7A/hF4Xp7/KPCq/Pjjcg3PyjV8E7h8mvc2/vqfy/UePPEzyfPcnef3aWDdoL/3Nt68ga5hEfEIcAbpn/PTwJikDZIW1zyrqyLi+xHxU9LqwrI8/reBmyLipoj4eURsBEZI4YfUkr5E0sERsSsixlcP/pe0uvD8iPhZRHxrhvl/JD/uZtLSyLURMRoRO4F/Bk4BiIgHImJjRDwREWPAZcCrKry/v4yIx/P7e5o8zy8Cm/L7+v0Kr1cch70PImJbRFwQEUuBlwDPBy6veTa7O+7/N2lxFlJg3yJp3/iN9OOzJCIeB36LtNi9S9JXJZ2Yn/cBQMDtku6R9K4Z5r+n4/5PJxk+FEDSYknXSdop6RHgC8BzK7y/h2aYvob02V4dET+q8HrFcdj7LCLuIy2mv2R8VMOzfAj4fEQs7LgdEhGX5nq+HhFnkxbh7yMtfRARuyPi9yLi+aSW8u8lHVdDPX9Fes+/GhGHkZY81DF9qs9jys8p74JbQ1rU/8Oa6px3HPaGSTpR0kWSlubho4C3Af+eH7IHWCrpwIZK+ALwRkmvl3SApIMkvVrS0tzKnifpENLegcdIi/VIest4zaSNbjE+rUfPzvP5iaQjgfdPmL4HeOEsX/ODpPreBfwN8Dnvg///HPbmPQq8Atgs6XFSyLcCF+Xpt5B2o+2WtLfumUfEQ8B5pECMkVr695O++2cA7wP+C/gxad353fmpL8s1PwZsAC6MiAdrKOnDwEtJW/m/Cnx5wvSPAX+WVzn+ZKYXk7Q8v4d3RMSTwMdJwb+4hlrnFXeqMSuEW3azQjjsZoVw2M0K4bCbFcJhNytEX8Mu6RxJ35P0gKS+7hqR9FlJo5K2doxbJGljPjZ8o6TD+1DHUZJulXRv7pl24QBrOUjS7ZLuyrV8OI8/VtLm/D1d32AfgIn1HCDpDkk3DriO7ZK+K+lOSSN53CC+n4WS1km6T9I2Saf1Ukffwp47Ofwd6WCQk4C3STqpX/Mn9Vo7Z8K4i4FNEXE8qV91P36A9gMXRcRJpKPR3pM/h0HU8gRwZkScTOpLf46kU0n7qj8REceROtSs6kMtABeSjvUfN6g6AF4TEcsiYjgPD+L7uQL4WkScCJxM+my6r6NfR9wApwFf7xi+BLikn0f9kI6W2tox/D1SH3FI3UW/18968nzXk44CG2gtwC8B3yF1ANrLL44oe9r31uD8l+Z/3jNJx/9rEHXkeW2n4xDhQfyvAM8hHXKsuuro52L8kTz9YIYdedwgLY6IXfn+btKJIvpG0jGko8E2D6qWvOh8J+kw1I2kQ2X3xS9OENGv7+ly0sE3411yjxhQHZB64N0saYuk1Xlcv7+fY0k9Hq/Kqzafyd2au67DG+iySD+VfetOKOlQ4Euk01M9MqhaIuLJiFhGallfDpw4/TPqJ+kNwGhEbOn3vKdwRkS8lLTK+R5Jr+yc2KfvZwGpW/EnI+IU0mHDT1tkn20d/Qz7TuCojuGledwg7ZG0BCD/He3HTPN55b4EXBMR433DB1LLuIjYB9xKWlxeKGn8lGX9+J5OB96Uzz93HWlR/ooB1AFApGPwiYhR4CukH8F+fz87gB0RsTkPryOFv+s6+hn2bwPH5y2sBwJvJR1gMUgbgJX5/krS+nOjJAm4EtgWEZcNuJYhSQvz/YNJ2w62kUL/5n7VEhGXRMTSiDiG9H9xS0S8vd91AEg6RNKzx++TzoKzlT5/PxGxG3hI0gl51FnAvT3V0Y8NHh0bF84Fvk9aL/xQn+d9LbCLdAaWHaQtu0eQNgrdD3wDWNSHOsbPWnM3cGe+nTugWn4NuCPXshX4izz+hcDtpHO5fRF4Vh+/p1cDNw6qjjzPu/LtnvH/0wF9P8tIZxW6m3Sar8N7qcNHvZkVwhvozArhsJsVwmE3K4TDblYIh92sEAMJe0cXxIFqSx3QnlraUge0p5a21AG91dJT2Hs4ZLUtH15b6oD21NKWOqA9tbSlDuihlq7D3oJDVs1sFrruVCPpNNL1t16fhy8BiIiPTfOcSjNbvnx55Tq2bKn/2InZzL8uY2NjDA0N9X2+ba0D2lNLW+qAmWvZvn07e/fu1WTTerk++2SHrL6ih9d7ysjISOXHpq7m9ZrN/M3aZHh4eMppvYS9krxBoU3rPGZF6iXslQ5ZjYg1pIvuVV6MN7P69bI1vo2HrJrZFLpu2SNiv6Q/Ar4OHAB8NiLuqa0yM6tVT+vsEXETcFNNtZhZg9xd1qwQDrtZIfoa9uXLl1c7fY5U+TaLU/yYFc0tu1khHHazQjjsZoVw2M0K4bCbFcJhNyuEw25WCIfdrBAOu1khHHazQjR+pppOW7ZsqXQaqdl0b53u9aZ7nemmVT3Vlbvh2lzilt2sEA67WSH6uhjfb52L4xMXuZs4K61Zm7llNyuEw25WCIfdrBBzfp29c13c6+FmU3PLblYIh92sEK1cjO92cbyOXnGz4Z52Npe4ZTcrxIxhl/RZSaOStnaMWyRpo6T789/Dmy3TzHpVpWW/GjhnwriLgU0RcTywKQ+bWYvNGPaI+Cbw4wmjzwPW5vtrgRX1llVd1YtLmJWu23X2xRGxK9/fDSyuqR4za0jPW+MjIiRNublZ0mpgda/zMbPedNuy75G0BCD/HZ3qgRGxJiKGI2K4y3mZWQ26DfsGYGW+vxJYX085ZtaUKrvergX+DThB0g5Jq4BLgbMl3Q+8Ng+bWYvNuM4eEW+bYtJZNddiZg1q5fXZZ3Nry+61uq43X/VmNlvuLmtWCIfdrBCtPOptrvCJM2wucctuVgiH3awQDrtZIRx2s0I47GaFcNjNCjHnd71VPZljv09i2e1JJr0Lz5rilt2sEA67WSHm/GJ8E+roGTfd5aLNBsEtu1khHHazQjjsZoXwOnsPfNSbzSVu2c0K4bCbFaKYxfjZ7P7qZpG8iZ52dVyC2rv9bJxbdrNCOOxmhXDYzQpRzDp7W01c9/buPGuKW3azQlS51ttRkm6VdK+keyRdmMcvkrRR0v357+HNl2tm3arSsu8HLoqIk4BTgfdIOgm4GNgUEccDm/KwTWLi5aDMBmHGsEfEroj4Tr7/KLANOBI4D1ibH7YWWNFQjWZWg1mts0s6BjgF2AwsjohdedJuYHG9pZlZnSqHXdKhwJeA90bEI53TIi2bTrp8Kmm1pBFJI2NjYz0Va2bdqxR2Sc8kBf2aiPhyHr1H0pI8fQkwOtlzI2JNRAxHxPDQ0FAdNbdGt5dbbuslqH256PmtytZ4AVcC2yLiso5JG4CV+f5KYH395ZlZXap0qjkd+B3gu5LuzOM+CFwK3CBpFfBD4PxGKjSzWswY9oj4FjDVsttZ9ZYzP03c3eZFYRsE96AzK4TDblYIh92sEA67WSEcdrNCOOxmhfDJKyZR95Fpg9zVNt3JMSbyEXnzm1t2s0I47GaFcNjNCuF19nluNl11vc4+v7llNyuEw25WCIfdrBAOu1khHHazQnhrfB80sZW7c6t6HZd2buo1vYW/PdyymxXCYTcrhMNuVgivs89RTVza2ZeLnt/cspsVwmE3K4TDblYIh92sEFWu9XaQpNsl3SXpHkkfzuOPlbRZ0gOSrpd0YPPlmlm3qrTsTwBnRsTJwDLgHEmnAh8HPhERxwEPA6saq9LMejZj2CN5LA8+M98COBNYl8evBVY0UWBpurls8nSXgJ6ojstFT/caM83fBqfq9dkPyFdwHQU2Aj8A9kXE/vyQHcCRjVRoZrWoFPaIeDIilgFLgZcDJ1adgaTVkkYkjYyNjXVXpZn1bFZb4yNiH3ArcBqwUNJ4D7ylwM4pnrMmIoYjYnhoaKiXWq0L0y2qT7dqYPNPla3xQ5IW5vsHA2cD20ihf3N+2EpgfUM1mlkNqvSNXwKslXQA6cfhhoi4UdK9wHWSPgrcAVzZYJ1m1qMZwx4RdwOnTDL+QdL6u5nNAe5BZ1YIh92sEA67WSF88oqWGeTloqvOezaXkKo6f/e2a55bdrNCOOxmhXDYzQrhsJsVwmE3K4TDblYI73qzWavrOnDWX27ZzQrhsJsVwmE3K4TX2ee5pq8N3y13q+0/t+xmhXDYzQrhsJsVwmE3K4TDblYIb423Rk219dw97frPLbtZIRx2s0I47GaF8Dq7zdpserHV3duuiRNolsItu1khKoc9X6P9Dkk35uFjJW2W9ICk6yUd2FyZZtar2bTsF5Ku3jru48AnIuI44GFgVZ2FmVm9KoVd0lLgN4HP5GEBZwLr8kPWAisaqM/MalK1Zb8c+ADw8zx8BLAvIvbn4R3AkfWWZmZ1mjHskt4AjEbElm5mIGm1pBFJI2NjY928hJnVoMqut9OBN0k6FzgIOAy4AlgoaUFu3ZcCOyd7ckSsAdYADA8Pe1+ITconsWzejC17RFwSEUsj4hjgrcAtEfF24FbgzflhK4H1jVVpZj3rZT/7nwLvk/QAaR3+ynpKMrMmzKoHXUTcBtyW7z8IvLz+kmw+qeMy0N3Oy+e1ezr3oDMrhMNuVgiH3awQDrtZIRx2s0I47GaF8MkrbE5xT7vuuWU3K4TDblYIh92sEF5nt1boZ7faifMrpVutW3azQjjsZoVw2M0K4bCbFcJhNyuEw25WCO96s3nLXWufzi27WSEcdrNCeDHe5hRfLrp7btnNCuGwmxXCYTcrhMNuVohKG+gkbQceBZ4E9kfEsKRFwPXAMcB24PyIeLiZMs2sV7Np2V8TEcsiYjgPXwxsiojjgU152MxaqpfF+POAtfn+WmBFz9WYDUBETHmT9NRtrqsa9gBulrRF0uo8bnFE7Mr3dwOLa6/OzGpTtVPNGRGxU9LzgI2S7uucGBEhadJeBPnHYTXA0Ucf3VOxZta9Si17ROzMf0eBr5Au1bxH0hKA/Hd0iueuiYjhiBgeGhqqp2qzCqZbPJ9qUX3irdt5Tfeag1o1mDHskg6R9Ozx+8DrgK3ABmBlfthKYH1TRZpZ76osxi8GvpJ/hRYA/xARX5P0beAGSauAHwLnN1emmfVqxrBHxIPAyZOM/xFwVhNFmVn93IPOrBAOu1khHHazQvjkFWYVzfVz2rllNyuEw25WCIfdrBBeZ7fizYXLRU98XjfcspsVwmE3K4TDblYIh92sEA67WSEcdrNCeNebWQPa2LXWLbtZIRx2s0J4Md6sokFeLrqO13TLblYIh92sEA67WSEcdrNCOOxmhXDYzQrhXW9mA9TPnnZu2c0KUSnskhZKWifpPknbJJ0maZGkjZLuz38Pb7pYM+te1Zb9CuBrEXEi6bpv24CLgU0RcTywKQ+bWUtVuWTzc4BXAlcCRMT/RMQ+4DxgbX7YWmBFMyWazT39vDZ8VVVa9mOBMeAqSXdI+ky+TvviiNiVH7ObdGlnM2upKmFfALwU+GREnAI8zoRF9kibFCfdrChptaQRSSNjY2O91mtmXaoS9h3AjojYnIfXkcK/R9ISgPx3dLInR8SaiBiOiOGhoaE6ajazLswY9ojYDTwk6YQ86izgXmADsDKPWwmsb6RCM6tF1U41fwxcI+lA4EHgnaQfihskrQJ+CJzfTIlmVodKYY+IO4HhSSadVWs1ZtYY96AzK4TDblYIh92sEOr1MrCzmpk0RtqY91xgb99mPLW21AHtqaUtdUB7amlLHTBzLS+IiEn3cfc17E/NVBqJiMk2+BVZB7SnlrbUAe2ppS11QG+1eDHerBAOu1khBhX2NQOa70RtqQPaU0tb6oD21NKWOqCHWgayzm5m/efFeLNCOOxmhXDYzQrhsJsVwmE3K8T/AZwbDRCBWBMiAAAAAElFTkSuQmCC\n"},"metadata":{"needs_background":"light"}}]},{"cell_type":"markdown","source":"To apply the boundary conditions, we will remove the rows associated to the fixed DOFs and add the contribution to the right-hand-side. First, we obtain the free and fixed DOFs.","metadata":{}},{"cell_type":"code","source":"DofsP = np.arange(0, 6)          # prescribed DOFs\nDofsF = np.arange(0, nDof)       # free DOFs\nDofsF = np.delete(DofsF, DofsP)  # remove the fixed DOFs from the free DOFs array\n\n# free & fixed array indices\nfx = DofsF[:, np.newaxis]\nfy = DofsF[np.newaxis, :]\nbx = DofsP[:, np.newaxis]\nby = DofsP[np.newaxis, :]","metadata":{"trusted":true},"execution_count":11,"outputs":[]},{"cell_type":"markdown","source":"We can re-order the matrices and vectors in blocks, such that it's easy to operate with the blocks corresponding with the fixed DOFs.","metadata":{}},{"cell_type":"code","source":"# Mass\nM_FF = M[fx, fy]\nM_FP = M[fx, by]\nM_PF = M[bx, fy]\nM_PP = M[bx, by]\n\n# Stiffness\nK_FF = K[fx, fy]\nK_FP = K[fx, by]\nK_PF = K[bx, fy]\nK_PP = K[bx, by]","metadata":{"trusted":true},"execution_count":12,"outputs":[]},{"cell_type":"markdown","source":"## Step 5: modal analysis\n\nUsing the matrices associated to the free DOFs, we can perform a modal analysis to get more information on how the structure will deform and determine the natural frequencies.\n\n$$ ( \\bold{K}_{FF} - \\omega^2 \\bold{M}_{FF} ) \\bold{\\phi} = \\bold{0} $$\n\nTo compute the natural frequencies and mode shapes we use the `eig` command, which is part of the NumPy package. For more information see: https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html","metadata":{}},{"cell_type":"code","source":"mat = np.dot(np.linalg.inv(M_FF), K_FF)\nw2, vr = np.linalg.eig(mat)\nw = np.sqrt(w2.real)\nf = w/2/np.pi","metadata":{"trusted":true},"execution_count":18,"outputs":[]},{"cell_type":"markdown","source":"We sort the frequencies and mode shapes in descending order:","metadata":{}},{"cell_type":"code","source":"idx = f.argsort()[::-1]  \nf = f[idx]\nModalShape = np.zeros((nDof, len(f)))\nModalShape[6:, :] = vr[:, idx]","metadata":{"trusted":true},"execution_count":14,"outputs":[]},{"cell_type":"markdown","source":"Let's see what these modes look like. Here, we plot the first 9 modes. Note that the system will have 19 x 3 = 57 modes (as many as the discrete system DOFs).","metadata":{}},{"cell_type":"code","source":"nMode = 9\nplt.figure()\nnCol = int(np.round(np.sqrt(nMode)))\nnRow = int(np.ceil(nMode/nCol))\nfor iMode in np.arange(1, nMode + 1):\n    plt.subplot(nRow, nCol, iMode)\n    Shape = ModalShape[:, -iMode]\n\n    # Scale the mode such that maximum deformation is 10\n    MaxTranslationx = np.max(np.abs(Shape[0::3]))\n    MaxTranslationy = np.max(np.abs(Shape[1::3]))\n    Shape[0::3] = Shape[0::3]/MaxTranslationx*10\n    Shape[1::3] = Shape[1::3]/MaxTranslationy*10\n\n    # Get the deformed shape\n    DisplacedNode = ([i[0] for i in NodeCoord] + Shape[0::3], [i[1] for i in NodeCoord] + Shape[1::3])\n    for iElem in np.arange(0, nElem):\n        NodeLeft = Elements[iElem][0]-1\n        NodeRight = Elements[iElem][1]-1\n        m = Elements[iElem][2]\n        EA = Elements[iElem][3]\n        EI = Elements[iElem][4]\n        if m == m_Jacket and EA == EA_Jacket and EI == EI_Jacket:\n            c = 'b'\n        elif m == m_Pile and EA == EA_Pile and EI == EI_Pile:\n            c = 'r'\n        else:\n            print(\"ERROR: unknown material. Check your inputs.\")\n            break\n        plt.plot([DisplacedNode[0][NodeLeft], DisplacedNode[0][NodeRight]], \n                    [DisplacedNode[1][NodeLeft], DisplacedNode[1][NodeRight]], c=c)\n    plt.title(\"Mode \"+str(iMode)+\": f = \"+str(round(f[-iMode]))+\" Hz\")\n    plt.axis('equal')\n\n# automatically fix subplot spacing\nplt.tight_layout()","metadata":{"trusted":true},"execution_count":19,"outputs":[{"output_type":"display_data","data":{"text/plain":"<Figure size 432x288 with 9 Axes>","image/png":"iVBORw0KGgoAAAANSUhEUgAAAa8AAAEYCAYAAADrpHnMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABCJ0lEQVR4nO2deZgU1bXAfweGTRZZZkQElR1EXMCJGqOIGwLuiRqMimuMSYya914UNSYmYlyiSVzjEg3ghgbXGBERJa6o4IIgIAOCgKwKyiaLnPfHqXZ6mtl6prurqvv8vq++W31vddWprtN17nLuuaKqOI7jOE6caBC2AI7jOI6TLm68HMdxnNjhxstxHMeJHW68HMdxnNjhxstxHMeJHW68HMdxnNhREMZLRDqLiIpIUQjX/oGIzBWRdSJyYq6v76SP64uTDq4v4RA54yUiC0Rks4gUp+S/HyhI5xzL01dEJojIKhGpy6S4PwJ3qGoLVX06g3J1FpFXRGSDiMwWkSOrOXaUiIys5Puh/OEySQT15SwRmSYiX4vIYhG5Kc3f2PUli0RQX4aJyBwR+UpEVojIaBFplcYpClZfIme8Aj4FTkt8EJG9gB1CkmUL8DhwXh2/vzswM3PifMejwPtAO+AqYJyIlGThOnEgSvqyA3ApUAwcABwB/F8a33d9yT5R0pc3gB+o6o5AV6AIGFn9VypQsPoSVeP1IDA86fNZwJjkA0RkRxEZIyIrRWShiPxWRBoEZQ1F5OagtTQfOKaS794vIktFZImIjBSRhpUJoqpzVPV+6qAgIjIPU8h/B836Jumeo4rz9gT6A79X1Y2q+gTwEfCjOp5vl0C+xLahjq3MsIiSvvxdVV9T1c2qugR4GPhBbW7C9SVnRElfFqnqqqSsb4HutbmJQteXqBqvKUArEdkjeOjDgIdSjrkdSNRWDsWU8Zyg7KfAsUA/oBQ4OeW7o4CtmJL0AwYB59dFUBG5S0TuqqxMVbsBnwHHBc36TZV8/zkRWVPF9lwVl90TmK+qa5PyPgzy00ZVPw/ka6GqLYCngLF1OVdIRFlfBpBU8XF9iQSR0hcROVhEvgLWYgbib0llri/VnDhSG7AAOBL4LXA9MBiYiDWnFegMNAQ2A32SvvczYHKw/zJwYVLZoOC7RUB7YBPQLKn8NOCVGuTqbj9X3e4nw7/RmcCUlLzrgFFVHD8K+AZYk7R9nfhNUo69HJiW/PtEeYuqvgTHnQssBopdX6KxRVxfOgLXAD1dX2reojz4+iDwKtCFlCY9Np7QCFiYlLcQe/gAuwCLUsoS7B58d6mIJPIapBwfddYBqYO6rbCaW1XcrKq/TXwQG5j+NPkAERkCXAIcoKobMyNqzoiUvoh5fl2PvVhWVXdsDnB92Z5I6QuAqi4RkRewVkn/Wt1FdoiFvkS12xBVXYjd/FDgyZTiVZgjxe5JebsBS4L9pcCuKWUJFmE1o2JVbR1srVS1Tk3i+iIi41P6g5O38VV8bSbQVURaJuXtQz0GbkWkFzAaOFVV42TIgWjpi4gMBu7DunM+qsv9VHNu15cMECV9SaEI6FbLY2skr/UljKZ7bZvB2EMsDfa/a9YHnx/C+k5bYko2Gzg/KPs58DHQCWgDTCKpCQs8A9yK1SYaBNc5tAp5BGgK9AnO0RRoUpf7yfDvNAW4OZDnJKypXlJNs35kSl5nyrs6WgW/30/Dfv55oC+HA18AA+p7P64vBaEvpwO7Bfu7A/8FnnR9qXmLbMsLQFXnqerUKop/BawH5gOvA48ADwRl9wETsEHG99i+ZjUcaIwp4GpgHNChiuvsDmykvNaxEZiTKBSRu0Xk7trfVcYYhg0WrwZuAE5W1ZV1PFd/oBfw1+SaWYbkzBkR0ZersYH+5yur4bq+RIeI6Esf4E0RWY+5zc/BHEIA15fqkMBKOo7jOE5siHTLy3Ecx3Eqw42X4ziOEzvceDmO4zixo0bjJSIPiAWMnJGU11ZEJopFM54oIm2CfBGR20SkTESmi0iYcxUcx3GcPKVGhw0RGYBNWhujqn2DvJuAL1X1BhEZAbRR1ctFZCjmpTMUC0p6q6oeUJMQxcXF2rlz5/rdiZM1pk2btkpVax2UU0QewMLnrEjSmbbAY5gL7QJsvsdqsZmct2I6swE4W1Xfq+78ri/RJl19yQWuM9GmLjpTY4QNVX1Vtl8m4ARgYLA/GpiMhf04ATNyCkwRkdYi0kFVl1Z3jc6dOzN1alUeqzFg3TpYtAh69YIG+dcTKyILaz6qAqOAO6gYuWAEMCmpwjMC05khQI9gOwD4e5BWSez1BWDjRvj0U+jTJ2xJMk4d9CXr5IXOJFNWBm3b2pYH1EVn6vqmbZ9kkJZh8bzAwqckz55eTHlIlfxl1Ch7Ca0KOwpQNFDVV4EvU7JPwCo6BOmJSflj1JgCtBaRqubE5A8jRsD++8O334YtSej40EQdOOYYuPDCsKUIlXo3E4JWVtqTxUTkAhGZKiJTV66s69y3iFAStHbjfh/ZpV4VnrzSF4B+/WD9epg7N2xJosAoLEBuMomWeg8sgsWIID+5pX4B1lIvLFThs89gt91qPjaPqavxWp6oHQfpiiB/CRVjfnWiPB5YBVT1XlUtVdXSkpJIdY+nz7ZtluZTt0QWqUuFJ6/0BaBnT0vHpMaELTy8pZ4mZWXwzTewabsVUAqKuhqvZ7EF3AjSZ5LyhwdN+wOBr2oa78oLmgRrwH38cbhyRJt6V3jyin33tXTKlFDFiDD1HprIu9Z6glmzLG3dOlQxwqY2rvKPAm8BvURksYich8W6OkpE5mJr49wQHP48FgusDIv/9YusSB01DjrI0l13rf64wsYrPMnssAOUlvqYVy2o69BE3rXWEyTG1o85pvrj8pzaeBueVkXREZUcq8Av6ytU7GjTxtJx4+Cii8KVJQIEFZ6BQLGILAZ+j1VwHg8qPwuBU4PDn8fc5MswV/lztjthvrJ5M0ybZt3OeeilWk+WJzyVvaWewvPPW9qyZfXH5TlRXowyPjRpAk2bQosWYUsSCbzCU0sOOwymT7cxjMQYmJMg0VK/ge1b6heJyFhsSkVhtNSTKSmBoiLo3TtsSULFq3uZoqQE1la30KjjpHD88ZY++2y4coSMD02kyeTJVllu2DBsSULFW16Z4quv4MMPw5bCiRN7723p88/D//1fuLKEiLfU02T9emjVKmwpQsdbXpni+9+H9u1rPs5xEhQXQ48eBT924aTJpk3QvzDnZifjxitTrFtn4X4cJx0aN4bx42s+znHAvFNXrLCtwHHjlSlatYItW2Dr1rAlceLE3nub3nyZOkfXcSph9WpLhw4NV44I4MYrUxx2mKXLloUrhxMvfvQjSydNClcOJx689pqliag+BYwbr0zx2WeWvv9+uHI48SIRn+6JJ8KVw4kHa9ZY2rVrqGJEATdemSIRZcMnmzrpsN9+0Ly5e485tWNpMKWtX79w5YgA/qbNFDvvbOn06eHK4cSLBg0sQsszz9R8rOO8+aalxcXhyhEB3Hhlih13tNSXuHDSpX17myfoODXRuzc0alReWS5g3Hhlir59LfWlxp10Oflkm7uzMHILEDtR49//ttSHJ9x4ZYzGjUHEB96d9NmyxdLnngtXDif6tGgBbduGLUUkcOOVSdq0Ke8+dJzakljawt2fnZpYuxbatQtbikjgxiuTtG8PX38dthRO3Nh7b+sGSix14ThVsWhR+UTlAscD82aSL7/0gXcnfYqKbFkdDy/m1ERJCRxySNhSRAJveWWSgQO/W5pb01731SloBg40LzLHqY7Fi2HmzLCliARZMV4iMlhE5ohImYiMyMY1IsmaNbB8Oa9O3kabNvDww2EL5MSGJUtgxgwf93KqRwT22CNsKSJBxo2XiDQE7gSGAH2A00SkT6avE0natgVVrvzlV3z1lfcCpUPBVngSDBhgqQforRWFpC/Tp8OQITDyd5th2zaemtuXq6+Gxx4LW7JwyUbLa3+gTFXnq+pmYCxwQhauEzm2HXEUAIs/tnEvDzBfOwq5wjNihM2yuO6NgQD889rFvPSSL8pdHYWmL9ddBy+8AH++diNf0JZJ03Zk5Ej4ZWEvyZkV49URWJT0eXGQl/fc/E9zYT33uFXsthvMnx+yQPGhYCs8U6faNK//zOzMFop47LZlHHUUXH552JJFmoLSl+99z9JvGu/Ikg+/4IrFFwEwbFiIQkWA0Bw2ROQCEZkqIlNXrlwZlhgZ49ln4ZY3DuA3PZ/ht//sRteubrzSoMYKT77pS4LFi22GxfQG+7JozkbOeHAwAKWlIQsWbQqqgrxihfnytG4NZ5wB775r+YnVdAqVbBivJcCuSZ87BXkVUNV7VbVUVUtLSkqyIEbuuPxyOPFE2LRje+b0Op4Tz2nDW2/BG2/AqlVhS5cf5JO+JHjuOZgzBy68ELZsa8jIG4rYfXcrS6yU4tSdfKnwLFwIu+8O//wnfPQRXH+9GbMDDghbsnDJhvF6F+ghIl1EpDEwDHg2C9eJDNOmmWt827awYAFMmWKh6gDWrQtVtLhQqwpPvnHxxZbOmAG/+AWMGmX7AM2ahSZWHCioCvLEiRb74I03YN994Z13TD922CFsycIl48ZLVbcCFwETgFnA46qa1xMTjjjC0nfesTXiVq40j1YoD1vnVEvBVXjAlvFq187CYb74IrRsCXffbWVuvKqloPRl9WrrOrzppvIVl7xSnKUxL1V9XlV7qmo3Vb0uG9eIEl99ZR5jRxxRHvR5+HBLEwssO1VTiBWe+fOtlXXFFTB+vFV4Nm4sfzm58aqaQtOXAQPg+9+3yFDt2pXH5S30KFEeYSMDvPyytbAWLLBB1b32spcS+CoXtaXQKjxHHmnpggXw6qu2Kkr79uXl8+aFIlZsKCR96dQJli2Dc86xKRR33WVz2Z/N27Zm7XDjlQFmzrQxr0GDLMjG/feXr+p+zz2hiuZElMQE9jvugJtvtvGuNWvKu5vLysKSzIkaHTqYV+oLL8Att8Cpp5pDz7hxYUsWLm68MsD3vgc9epgyXXyxBQfv1cvKvOXlpLJtm413JQzVXntZBWjtWvj97y3voovCk8+JFpMnW89O27YwYQIcfzxs3mzeqq+8ErZ04eHGKwN88YXVjIqL4cEH4ZprrBXWty/ss0/Y0jlRo6wM1q+3biCAWbOgf397Ga1ebY4bRb7egxOQmG7z9dcwd66FwWzY0PIS3qmFiBuverJ1K3z8sQ22r1plL6F337VWWLdu8PnnYUvoRI3EJNNLLrFKjgh07AjHHQeTJvlCuU5FTj/ddKRlS6so33tveZdht27hyhYmbrzqyaJF5YHAmze35XYWLDBX1o4drZbkOMlMnmy6suee5WOinTrB+edbTXrxYg8u75SzbBnsvDO8+aZ5oR56aPmqKDGee11v3HjVkzZtLD31VNsmTIBTTrEuxKeesm4g9xxzEqjCAw/YJPbx4y1ywnXXmd4cdphVfnbc0RZWdhyA5cvNePXubQEQeveGCy6wsmXLwpUtTPwvUk+WLrX0uOPMWeOKK+Cgg6zmnCi7447w5HOihaqNVxQVmc7062feZPvvb92IxcXlq6M4DtgE9hkzrKV++OHm2NO4sZVdc02oooWKDwvXkzFjLL3wQhuEB5tIuHUrtGhh3Yc//GF48jnRYs0a8xwbORJ22gn+9CeLDt6li7XSN2wwQ+Y4CYqKbNtjD2uRr1tXPhyxbZtViBKeq4WEt7zqyYIFlh58MNx5p63m/sUXcNJJ5dE2PDivkyDRhdyrF5x9tnkaPvKIxan79lszXi+8EKqITsT43vfMEWzcOPj5z23sq2VL6+XZvBneey9sCcPBjVc96dvX0ksusdr0W2/ZDPhx42zuF1ifteNA+eTj7t2txvzBBxbYOTnUj8etc5Jp2dK6CseMgaOPhl12sbGv3/zGug8feihsCcPBjVc9WbjQPICOOcaiarz9ttWORKxbCNx4OeUkjNcjj9jAe2kp3HabpWPHmhfi+eeHK6MTLVq0sB6es86CQw6xllfnzuYsdswx8OijhblquxuveqBqHoUbN1og3qlTK05KbtTIxr/ceDkJRo2y9PrrbSrFffeZx9gzz8CPfwxNmhTmi8ipmpkz4csvzXiNH2/xUxOceaa9X156KTTxQsMdNuqBCBx7rA2aJl5KqVx6qdWwHQdg8GBbaeCee6z7J5WLL4a99869XE50GTbMJq//85/bO2YMHWrOYpXpUr4jqhq2DJSWlurUqVPDFsOpAhGZpqqRWZje9SXaRE1fwHUm6tRFZyJhvERkJZCLELbFQC58//LtOruramSWog30ZT25ufdckKvnmAuKgeZR0hfY7h2TD793vt1D2u+YSBivXCEiU3NRI8y360SRfLp3v5fcEgcZa8LvwR02HMdxnBjixstxHMeJHYVmvO716+QN+XTvfi+5JQ4y1kTB30NBjXk5juM4+UGhtbwcx3GcPMCNl+M4jhM78tJ4icgpIjJTRLaJSGlK2RUiUiYic0Tk6KT8wUFemYiMqON1632OlPM9ICIrRGRGUl5bEZkoInODtE2QLyJyW3Dt6SLSv77XjwIi8mcRmR3c01Mi0jqpLGvPMhfERc5kRGRXEXlFRD4O/mOXBPmR08t81J2oy5cgJ3qiqnm3AXsAvYDJQGlSfh/gQ6AJ0AWYBzQMtnlAV6BxcEyfNK9Z73NUcs4BQH9gRlLeTcCIYH8EcGOwPxQYDwhwIPB22M8hQ89yEFAU7N+YdL9Ze5Y5uq9YyFmJ3B2A/sF+S+CT4FlETi/zTXeiLl+u9SQvW16qOktV51RSdAIwVlU3qeqnQBmwf7CVqep8Vd0MjA2OTYdMnKMCqvoq8GUl9zA62B8NnJiUP0aNKUBrEelQn+tHAVV9UVUToWqnAJ2C/Ww+y1wQFzkroKpLVfW9YH8tMAvoSAT1Mg91J+ryfUcu9CQvjVc1dAQWJX1eHORVlZ+Jc2ea9qq6NNhfBrTP8fXD5FysdgbZfZa5IC5yVomIdAb6AW8Tfb3MB92JunyVki09iW1UeRF5Cdi5kqKrVPWZXMsTBqqqIhL7uQ61eZYichWwFXg4l7I5lSMiLYAngEtV9WtJCneeS7103Yk22dST2BovVT2yDl9bAuya9LlTkEc1+Zk4dyZZLiIdVHVp0KxekePrZ5yanqWInA0cCxyhQQc52X2WuSC2z0tEGmEvpIdV9ckgOxS9LDDdiZXOZF1Pwh7Yy+bG9g4be1JxoHY+NghaFOx3oXwgdM80r1Xvc1Rx3s5UdNj4MxUHPG8K9o+h4oDnO2H//hl6hoOBj4GSlPysPcsc3Vcs5KxEbgHGAH9LyY+cXuab7kRdvlzrSeg3maUf7iSsz3QTsByYkFR2FeaxMwcYkpQ/FPOImYd1OdTluvU+R8r5HgWWAluC+zkPaAdMAuYCLwFtk5TlzuDaH5FktOO8YYPpi4APgu3uXDzLHN1bLORMkflgQIHpSc9kaBT1Mh91J+ry5VJPPDyU4ziOEzsKzdvQcRzHyQPceDmO4zixw42X4ziOEzvceDmO4zixw42X4ziOEzvceDmO4zixw42X4ziOEzvceDmO4zixw42X4ziOEzvceDmO4zixw42X4ziOEzvceDmO4zixoyCMl4h0FhEVkZyvXyYiJ4nIIhFZJyL9cn19J31cX5x0CFlffiAicwN9OTHX1w+TyBkvEVkgIptFpDgl//1AQTqHJBoiMqkOSnozcJGqtlDV9zMoy74iMk1ENgTpvtUcO1lEzk/JGygiizMlT1hETV9E5GwR+TZ4mSS2gWmcwvUli0RNX4JrdxWR50RkrYisEpGb0vj6H4E7An15OoMydRaRVwJ9mS0iVS76KSKjRGRkJd/PqkGPnPEK+BQ4LfFBRPYCdghPHBCR04FGdfjq7sDMDMvSGHgGeAhoA4wGngnyC5Go6ctbwcsksU1O47uuL9knMvoSPIOJwMvAztgKwg+lcYqM60vAo8D72PpbVwHjRKQkC9epM1E1Xg8Cw5M+n4WtyvkdIrKjiIwRkZUislBEfisiDYKyhiJyc1CLmY+t0pn63ftFZKmILBGRkSLSsCphRGRH4PfAZbW9ARFpIiLrsBVaPxSRebX9bi0YiK2q+jdV3aSqt2GLuR1el5OJyPdTWgrfiMiCzImbdSKlL3XB9SWnRElfzgY+V9W/qOp6Vf1GVafX5iYCHekK/Dt4Dk1qdfc1n7cn0B/4vapuVNUnsAUif1TH8+2Soi8bRKTeC0lG1XhNAVqJyB7BQx/G9rWR24EdsYd3KKaM5wRlPwWOBfoBpcDJKd8dBWwFugfHDALOp2r+BPwdWJZaEDT3R6TmBy+JFsHHfVS1W2UnFpHpIrKmiu2uKuTZE5iuFVcSnR7kp42qftdSwGrmb2M1r7gQNX3pF7zYPhGRq5O7TlxfIkGU9OVAYIGIjA90ZnLQEgRARO6q6rkGOvIZcFzwPDalHhPoW1X68lwVMu0JzFfVtUl5H1J3ffk8uScCeAoYW5dzpZ44UhuwADgS+C1wPTAYa1YXYctKd8Zqp5uBPknf+xkwOdh/GbgwqWxQ8N0ioD2wCWiWVH4a8EoV8pRiS1gXBddWoCiN+1Gge4Z/o6uBsSl5DwPXVHH8ZGADsCZpWwcsruTYvwPPAQ3C1oWY6ktXoAtWMdwL+Bi4wvUlGlsE9eVFYAswBGgM/AaYDzRO534y/BudCUxJybsOGFXF8aOAb1L05WsqeVcClwPTkn+fum45945JgweBV7EXwZiUsmJs/GlhUt5CoGOwvwuwKKUswe7Bd5eKSCKvQcrxAATdBHcBl6jq1qTjw2Yd0ColrxWwtpJjE1ysqv9IfBBzIqhQ2xSRn2FdTAeo6rZMCJpDQtcXAFWdn/TxIxH5I/ZCur62N5IFXF+2JxL6AmwEXlfV8QAicjNmWPfAWjthUBd9uVlVf5v4IOb48mnyASIyBLgE05eN9RUyqt2GqOpC7OaHAk+mFK/Caiu7J+XtBiwJ9pcCu6aUJViE1YyKVbV1sLVS1cqaxK2wltdjIrIMeDfIXywih9ThtrZDRGam9Acnb3dX8bWZwN5S0ZruTT0GboP7uRY4QVW/rut5wiIi+lKpaNj4UkZwfckMEdKX6ZiOZIWgO7IqfRlfxddmAl1FpGVS3j7UT196YY5Cp6pqVYY8PcJqvtemGQx0A0qD/e+a9cHnh7C+05aYks0Gzg/Kfo5113TC+uQnkdSExTyvbsWMU4PgOodWIotgHkCJ7XvBeTpS+2Z9NrqBGmO1vUuAJsBFwedKZcK6gc5PyRtI0A2E/RGXAoPDfv5x1pfg2CFA+2C/NzADG/iu7f24vhSWvvTCumiPxLorfw3MS+P98t39ZPh3moJN22gKnIR1BZZUcewoYGRKXmfKu1JbBb/fTzMpY2RbXgCqOk9Vp1ZR/CtgPdY//DrwCPBAUHYfMAFrdr/H9jWr4dgf+mNgNTAO6FDJ9VVVlyU2YGVQtFxVN8N3NZsr63iLdSK49onYfawBzgVOTMhUB47A+urHJdXKsuF+m1XC1peAI4DpIrIeeD44158Sha4v0SEK+qKqc4AzgLuDY08Ajk96v9xdTYs6mwzDep1WAzcAJ6vqyuq/UiX9MSP91+SWX30FlMBKOo7jOE5siHTLy3Ecx3Eqw42X4ziOEzvceDmO4zixo0bjJSIPiMgKEZmRlNdWRCaKRTOeKCJtgnwRkdtEpCyIBNA/m8I70cR1xkkH1xenLtTosCEiA7BJa2NUtW+QdxPwpareIBbqpo2qXi4iQzEvnaHAAcCtqnpATUIUFxdr586d63cnTtaYNm3aKlWtdVDObOuM60u0iZq+gOtM1ElXZ4DazfPCfPZnJH2eA3QI9jsAc4L9e4DTKjuuum2//fbTvOF3v1M999ywpcgowFRNf55I1nQmr/QlwZYtqoccovrII2FLUm+ipi+aTzrzySeqhx+u+tprYUuSUeqiM3Ud82qvqkuD/WXYnA+wybvJs6cXUx5SpQIicoGITBWRqStX1nX6QASZOBHmZTIgeN5QL53JW31J8NRT8Npr0CQjgcHzAX/HVMXLL8Onn9Z8XJ5Tb4eNwGqmPVlMVe9V1VJVLS0pidQyMXVn2zaYMQN69QpbkkhTF53JS31JsG0bXHMNdO8OJ5wQtjSRw98xSbRta+nChdUfVwDU1XgtF5EOAEG6IshfQsWYX50ojweW/3zwAaxdC6tXhy1JFHGdqYo774SPP4aBA6FhRpcJizOuL5Wx446WTpgQrhwRoK7G61lsATeC9Jmk/OGBR9CBwFdJTf/8Z0nwHzr++HDliCauM1UxeTI0awbXhxl4PnK4vlRGURE0bw677BK2JKFT45IoIvIoFpSzWEQWYysK3wA8LiLnYQE+Tw0Ofx7zAirDgk2es90J85kXXrD04IPDlSNkXGfS4J134Omn4bzzoLg4bGlCwfUlTZo1g0WZCcweZ2o0Xqp6WhVFR1RyrAK/rK9QsWXqVBCB3Xar+dg8xnUmDe6808a8hg0LW5LQcH1Jkw0bYP78mo/Lc6K8GGX8aNcO9tkHGnjgEqcWqNpY1777wuGHhy2NExdKS62SXOD4WzaTTJoEX8duXT4nLCZPttb63nuHLYkTJ+bPh2nTwpYidLzllSlUrfunT5+wJXHiwsxgCaxTT63+OMdJpqQENtd1Kbb8wVtemWLhQti6FfbbL2xJnLjw2mvQqBEcdVTYkjhxYqedTG8KHDdemeKddyxdvz5cOZz48NJL0LixbY5TW+bOheXLw5YidLzbMFMkZr4femi4cjjxoXdva607Tjr07evdhrjxyhyvvWapR0hwassnn5RHTHCc2rJwoUXyKXC82zBTbNxoaevWoYrhxIiiIq9BO+mzcaM5hxU4brwyRccgsHWHDuHK4cSHo46CTZvClsKJG506QY8eYUsROt5tmCnefNNSd9hwasuCBbBiBaxaVbChoZw64HO8AG95ZY7ERNNmzcKVw4kPZ59t6auvhiqGEzPatIGePcOWInTceGWKxNjFV1+FK4cTH04/3ebr3H572JI4ceKrr6Bly7ClCB03XpkiERZq2bJw5XDiQ5Mm5uAzZQp8+23Y0jhxYc0a62oucNx4ZYqhQy1t3jxcOZx4ceut8M038O9/hy2JEwe2brVQdL5moBuvjLEiWOj1pZfClcOJFyedZC7zl14atiROHJg3z9KFC8OVIwK48coUiYC8vlSBkw5Nm8IPfgCLF3vIH6dmvvzS0t13D1eOCODGK1PsuaelS5YAvjKKkwb33mtjXv/3f2FL4kSdL76wdN99QxUjCmTFeInIYBGZIyJlIjIiG9eIHI0bW6vrjTe44AILdThjRthCxYeC1JkEPXvaQqZjx8KWLWFLEwsKVl/eestSn0+aeeMlIg2BO4EhQB/gNBEpjEWu9tyTRd/uwn33WUX6nnvCFigeFLTOYPZq27UjbTD+P/8JW5zIU4j68s03toDySx/vAsDtU0oZMsR8NwqVbLS89gfKVHW+qm4GxgInZOE6kePNL3qxfu4Shg6Fc86Bu++GOXPClioWFKzOqForvcUl57Ok4a7898d3MWwYnHJKuQ+Qsx0Fpy+JxZMnPL0BgCv/visvvAAXXFC4CxNkw3h1BBYlfV4c5FVARC4QkakiMnXlypVZECO3PPYYPL70YOY068cTT8ANN8AOO8Cvfx22ZLGgRp3JN31JcP/9sG4dfLO1iAd3/BUrKOGZJ7YybpxNA3MqpeDeMX362DDXxqJWfMA+rKMFJSXwj3/AMccU5hh7aA4bqnqvqpaqamlJSUlYYmSEH/8Yhg2DB9teyu0HPcqxx9rK7jvsAOPHeySXTJBP+pLggw/goots8W1V6Pr333DKpofp3beIgw7y1VLqS77pzAEHwL3yM7556wOuvVZo1cryX3wROncuvEAt2TBeS4Bdkz53CvLylqVLLW3XLqhFf2M1oYTn8zffhCdbTCg4nZkxA77/fesyfO452GUXeOQR84T+8EM4+uiwJYw0BacvYFO8tmwxvbn66vIpXwCrV8PFF8OGDeHJl2uyYbzeBXqISBcRaQwMA57NwnUiwzPP2BqUP/yhRfo5/3z46COrTbdu7V6ttaDgdGbYMKvUNGxoFZ1hw+D5521ThcMOC1vCSFNw+rJlC8yaBSUl1tV8ww1wyCEV177t2rWw4oJnfEkUVd0qIhcBE4CGwAOqOjPT14kSbdrAEUfAE0/Y4Oktt5jnfO/e1lc9fXrYEkabQtOZSZNg5kwYPBjefde8yK66yl5Qo0bZC+iAA8KWMroUmr4APPigTSG9+mp46CF45RWLRPeLX0BZGbzxhjl0FFSMBFUNfdtvv/007owcqWp1ZtXu3VWLilTff1/1d79TbdBAdePGsCWsO8BUjYCeJLY468vatapduqj26KG6YYPqZ5+pHnig6U3r1qrNm6sedVTYUtaPqOmLxlxnpk6190mjRqYnu+yieuONql9+qTp7tqqI6ogRYUtZP+qiMx5hI0MkFjbt3dtqQn/4g3UX7rGHrdj9ySehiudEgK1brWvn008tpOH999u85COPhO99z4KFr19fsSvIcU47zXQnMX/988/hyivt3ZJYRrAQQ2P6SsoZorTU0kWL4MADLdLPvHnlK6Tcfjvcd1948jnh88QTkPDYvummimWNGkGDBlbR2bQp97I50aVjR6vwbN1q3cs77WRzABcutC7EkhJo3z5sKXOPG68MkYiXuXGjefy0bVsxgksiqotTuAwYUPFzr162IsqRR5redO5srfUXXwxDOieKLF4M//0vXH65OYbddx+8/755p44caccU6kLc3m2YIe6/39Jt28zD8PzzbQLh6NGWf+21oYnmRITEIttXX21OGWVl5rRx5plw440Wc3XkSGuBOQ6Yo4aqvU/GjbOpOMOGWQX5zjttSkXv3mFLGQ7e8soQa9eah+G338LmzWasWraEZwMH3k6dwpXPCZ+E8TrwQHjtNRg0yHTlX/+yLqGOHa3MccCM1qhR5hLfrZvl3XsvnHGGBUFYtgweeCBUEUPFjVcGWLMGnnwSzj3XakInnwzHHWdzdhYvtmPceDmrV1vaurVF1Xj1VTjqKHNvXrfOug0dJ8GFF5qj15YtMHCgrVlaVGSVnP/8x8KHFfJkdu+gyAA//amNWZx6Kpx4ojX1X33VJi0vWGAKt9NOYUvphM0dd1h63nkWUmzMGPj5z60LsUEDmDs3XPmcaNGokaU77GDpN99YRXnzZvv87behiBUZ3HhlgCeftPTww63/efx465eeMMFC/nTo4O7PDuyzj6UNGtig+9/+Br/7nXmOffutpZdfHqqIToS45Zby1tXkyfD66zBihHmsHnqodTWPGxe2lOHhxquevPeeOWkMG2aD7b16wcSJ8OijVr5kibnPawGvu+MYv/+9tbIGDbLuoI0bTX8Sc3XA5vA4DpjhOvBA8zYEmD0bhg+3OYH/+Y/N8/rjH+39U4i48aon998PTZvCXXfBFVfYZMEjj7S8BM2bF1jYFqdSGje2sE+vv26f33rLlrMoK4OHH7buxHHjygM6O86hh1orfdEiG5Jo2tTmCzZvbq32mTPtcyHixqsebNhgL53Bg+G226B7d+s6fPZZqyH95S92XCE37Z2KHHywtbZuvtkG4Zs3t2DOP/mJdRlu3gx//WvYUjpRYcAAa1mdeqpVcv71L9g1iKd/yinW+vrDHwqz9eXGqx5cdZW5Pz/9tClQt242433pUrjnHnPWaNrUak+OAxaJZds2+M1vYMgQC8y7115W1qOHvaTuuqvcM9EpbL7/fRsjnTLFKjzJ75KGDW3O4MyZ5ePuhYQbr3owZYp5BP3hDxa+ZeJEOP30cu+ghQutJVZIyxQ41ZOYr3PCCVbpad26YvkVV9j46NSpuZbMiSKNGlnrvHNnuOSS7ctPPdVaX4UYwcfnedWDl182j5+WLSsvf/ppj1PnVKRvXxvTqmrqxN57m9NGVTrlFBaNGpnTV1FR5ePmDRvC228Xpr6IRsANTkRWAgvDliNDFAOrwhYigxQDzVU1MuuoB/qynvz7nfPlfnaPkr5ATt8xuXqO+XadtHUmEsYrnxCRqapaGrYcmSKq9xNVuepKvt1PoZKr55hv16kLPublOI7jxA43Xo7jOE7scOOVee4NW4AME9X7iapcdSXf7qdQydVzzLfrpI2PeTmO4zixw1tejuM4Tuxw4+U4juPEDjdeGUREBovIHBEpE5ERYctTEyKyq4i8IiIfi8hMEbkkyG8rIhNFZG6QtgnyRURuC+5vuoj0z4GMfxaR2cH1nhKR1kllVwSyzBGRo5PyY/Mc4iSrU46InBL8Z7aJSGlKWdb0MtP6IiIPiMgKEZmRlBeZ/3+1qKpvGdiAhsA8oCvQGPgQ6BO2XDXI3AHoH+y3BD4B+gA3ASOC/BHAjcH+UGA8IMCBwNs5kHEQUBTs35gkS5/gN24CdAl++4Zxeg5xktW37Z7dHkAvYDJQmpSfNb3Mhr4AA4D+wIykvMj8/6vbvOWVOfYHylR1vqpuBsYCJ4QsU7Wo6lJVfS/YXwvMAjpico8ODhsNnBjsnwCMUWMK0FpEOmRZxhdVdWvwcQrQKUmWsaq6SVU/BcqwZxCn5xAnWZ0kVHWWqs6ppCibeplxfVHVV4EvK7mHSPz/q8ONV+boCCxK+rw4yIsFItIZ6Ae8DbRX1aVB0TKgfbAf9j2ei9X8qpMlbBnTIU6yOrUjm3qZK32J6v+/Ah6Y10FEWgBPAJeq6teSFAFUVVVEsjqfQkReAnaupOgqVX0mOOYqYCvwcDZlcZwEtdHLfCcX//+64sYrcywBdk363CnIizQi0ggzXA+ramJVoOUi0kFVlwbdAiuC/Kzco6oeWYOMZwPHAkdo0PlegyxxeQ6x1JlCoSa9rIJs6mWu9CWn//+64t2GmeNdoIeIdBGRxsAw4NmQZaoWsSbW/cAsVf1LUtGzwFnB/lnAM0n5wwOvowOBr5K6F7Il42DgMuB4Vd2QIuMwEWkiIl2AHsA7xOs5xElWp3ZkUy9zpS+R+f9XS5jeIvm2Yd44n2AeQVeFLU8t5D0YUGA68EGwDQXaAZOAucBLQNvgeAHuDO7vI5K8rLIoYxnWz56Q7+6ksqsCWeYAQ+L4HOIkq28VnttJ2JjPJmA5MCGpLGt6mWl9AR4FlgJbgvs5L0r//+o2Dw/lOI7jxA7vNnQcx3Fihxsvx3EcJ3a48XIcx3Fihxsvx3EcJ3a48XIcx3Fihxsvx3EcJ3a48XIcx3Fihxsvx3EcJ3a48XIcx3Fihxsvx3EcJ3a48XIcx3Fihxsvx3EcJ3YUhPESkc4ioiKS8/XLROQHIjJXRNaJyIm5vr6TPiHrSy8R+UBE1orIxbm+vpM+/n4Jh8gZLxFZICKbRaQ4Jf/9QEE651ieuwPFSGybRGRtGqf4I3CHqrZQ1aczKNe1IvKRiGwVkWtqOPYaEXmoknwVke6ZkikMIqgvIiIjRWSJiHwlIpNFZM80TnEZ8IqqtlTV2zIo10UiMjXQ31E1HHu2iLxeSf4CEanLAo2RIYL60kRE/ioin4vIahG5K1ggtrYU7PslcsYr4FPgtMQHEdkL2CEMQVT1wkAxWqhqC2z9m3+lcYrdgZlZEK0Me9H9JwvnjhuR0RfgFOBc4BCgLfAW8GAa38+WvnwOjAQeyMK540aU9GUEUAr0BXoC/YHfpvH9gn2/RNV4PQgMT/p8FjAm+QAR2VFExojIShFZKCK/FZEGQVlDEblZRFaJyHzgmEq+e7+ILA1qyCNFpGFNQolIc+BHwOja3ISIzAO6Av8OWm1NavO92qCqo1V1PJBOK7BKRGRNUutyfRi10HoQJX3pAryuqvNV9VvgIaBPbW5CRF4GDgPuCJ5Dz9p8rzao6pNBzfyLTJxPRD5M6ZFQERmYiXPngCjpy3HAbar6paquBG7DKj81Uujvl6garylAKxHZI3jow7CXQDK3AztiD+9QTBnPCcp+ChwL9MNqNSenfHcUsBXoHhwzCDi/FnL9CFgJvJrICJr5d1V2sKp2Az4DjgtabptSjxGR54IHW9n2XC1kygiq2jqpdXkr8BqwJFfXrydR0pexQDcR6Rl0/5wFvJAoFJERVT1XVT0c+90vCp7FJ6nHBPpWlb5Mr0KmjKOq+yTpy/9gqwa/l6vr15Mo6QvYCsXJ+51EZEfw90tNJ43UBiwAjsSaztcDg4GJQBG2ZH1noCGwGeiT9L2fAZOD/ZeBC5PKBgXfLQLaY0t3N0sqPw0bZ6hJtknANXW5nyz+Xg/VJBNwTfB7rUnZFOiecuyPA5lLwtaFOOoL0Bj7cyr2AvsU6JLG/UwGzs/i7zUSGFXDMWcHsqfqy7ZUXQYOBlYAPcPWhZjqy0jgDaAE2Bl4OzhXh3TuJ4u/V2TfLzn3jkmDB7EWThdSmvRAMdAIWJiUtxDoGOzvAixKKUuwe/DdpSLfVXgapBy/HSKyGzAQq3XFkcdV9YzkDBHRlM/9gDuAQWpdGHEiKvryO+B7wK7AMuAM4GUR2VNVN6RxP2EzRVUPTs4QkQUpn3cFHgfO0kpaiREnKvpyHdAa+AAzevdhrbXltb2RiJDz90tUuw1R1YVYrXUo8GRK8SpgC6YoCXajvBm6FHt5JJclWIQpSbFaU7a1qrZS1Zo8ws4E3lDV+endSfWIyPiUsYPkbXwmr1WDHDsBTwO/VNX3c3XdTBEhfdkXeExVF6vqVlUdBbShluNeNSHbe78mb9kYuK9KjmaYvvxNbWwkVkRFX1R1o6pepKodVbUrNiY5TVW31fXeksnn90tkjVfAecDhqro+OVNtIPxx4DoRaSkiu2P97ol+68eBi0Wkk4i0wTx6Et9dCrwI3CIirUSkgYh0E5FDa5BlONaXnVFUdYgmeTOmbEOq+p6INBKRptgzLBKRptUMCleL2PyUccBDqvp43e4kEkRBX94FThGR9sGxZ2I18bJM3KCmeL+mbFVWwESkKNCXhkDDQF/q0/PyADBbVW+qxznCJnR9EZGOIrKLGAcCVwO/z9QN5vP7JdLGS1XnqerUKop/BawH5gOvA49Q7gZ8HzAB+BAbRE6tWQ3HxiY+BlZjP2yHquQQke8DnajERT6oCd9dy1vKJPcBG7H+9KuC/TPreK5OmGv3pSk1s91q+mKUiIi+3Bic5wOs3//XwI9UdQ2AiFyZyxpvEr/FdGQE1pW5kfRcslMZBpyUoi+HZEDOnBERfekGvBlcazQwQlVfTBT6+6VqJBhEcxzHcZzYEOmWl+M4juNUhhsvx3EcJ3a48XIcx3FiR43GS0QeEJEVIjIjKa+tiEwUi2Y8MfC4SQQlvU1EykRkuoj0z6bwTjRxnXHSwfXFqQs1OmyIyABgHTBGVfsGeTcBX6rqDSIyAmijqpeLyFDMS2cocABwq6oeUJMQxcXF2rlz5/rdiZM1pk2btkpVS2p7fLZ1xvUl2kRNX8B1JuqkqzNA7cJDYSFTZiR9nkMQvgRzAZ0T7N8DnFbZcdVt++23n+YNH3yg+tRTqlu2hC1JxgCmavphZbKmM3mlL19/rfr446rz54ctScaImr5ovulMgn/9S3XOnLClyAh10Zm6jnm1V5uMBxYCp32w35GKYVAWUx5SpQIicoHY+kJTV66MWySiahgzBn7yE2jgw4kp1Etn8lZfVq+GU0+Fl14KW5Ko4e+Y6li3Ds4+G26+OWxJQqPeb9jAaqY9WUxV71XVUlUtLSlJr7UYaebMgR493HhVQ110Jm/1pUMHEIElcQngn3v8HVMJTzwB69fDWWeFLUlo1PUNu1xEOgAE6YogfwkVY351Ij7LamSGyZPdcFWO60xlNGoETZrAczlbnSIuuL5Ux2WXQbt2cNBBYUsSGnV9yz6LrVNEkD6TlD88KU7XV0lN//xnyxarDRUX13xs4eE6UxVNm1o3kJOM60tVrFhh2z77WKu9QKkxMKeIPIotBVIsIouxoJE3AI+LyHnYcgCnBoc/j3kBlQEbKF+8rTCYO9fSk1PXpissXGfSZO+94euvw5YiNFxf0uTJIJTitdeGK0fI1Gi8VPW0KoqOqORYBX5ZX6Fiy/gg3uqGOC3blHlcZ9Jk4cKCHvNyfUmTv/zFWlz77x+2JKES5cUo40ezZpYeUOO0E8cpZ6+9Crrl5aRJu3awbRsUFfbr2z0LMsmsWZZ26RKuHE782LgxbAmcuDBvnhmwAseNVyaZPdtSVywnHVasgE2bwpbCiQurV3tLHTdemSXhttqoUbhyOPHigAOgdeuwpXDiwn77wS67hC1F6BR2p2mmmTzZ0gJ2X3XqwCuveE3aqT1Ll0K+RQypA97yyiQdgyg1W7eGK4cTL1q3hlatwpbCiQv9+8M33xR8hceNVybp3t3SL74IVw4nXmze7BPbndrz05+CKjz6aNiShIobr0yyKIgXWlYWrhxOvJg1C1atClsKJy7st5+lo0eHK0fIuPHKJEOGWOrdhk46tG8PBx4YthROXGjfHvr1sxZ7AePGK5MkPMYmTQpVDCdmLFpU8OMXTpr06QPTppWHpCtA3HhlksTk5M8/D1cOJ15s2WKR5R2ntiS6Dt99N1w5QsSNVybp3h122CFsKZw4sXWrhfpJvIwcpzaceKKly5eHKkaYuPHKJCLw7bfw/PNhS+LEhYRn6syZ4crhxIvddrO0gD0OfZJyptljD5uD4Ti1IRHMuUePcOVw4kViwdvS0nDlCBFveWWa4uKCXt7CSZPGjS1dsyZUMZyYkVh2acuWcOUIETdemWbDBli71roPHacmEsZr/vxw5XDiRaKbuYADOrvxyjTHH2+pL3Hh1IYGDSysWLduYUvixIn337f0ssvClSNEsmK8RGSwiMwRkTIRGZGNa0SNUaOCoa7EvIuPPgpTnNhRiDrzHV9/Xb4Kt1MrClpfAO680xzEevYMW5LQyLjxEpGGwJ3AEKAPcJqI9Mn0daLEXXfBOedY5Xl1686WqRqqTHGiEHWmAp07u76kQcHriyrblq9gSavefPBx47ClCY1stLz2B8pUdb6qbgbGAidk4TqRQBXGjrVK0NKl8Ivb9wBg6uwWHuKw9hSUzqSy4sDj+HbVl4z8/RYOPxwWLw5boshTUPoyYwY0bGhLeA0bBrecP4sGK5Zz7VcXc8klYUsXHtkwXh2BRUmfFwd5FRCRC0RkqohMXRnjtWkeeQReew3OO88MWYu2Finh5+dtokcPWLYsZAHjQY06ky/6ksz998OPfwxXjOpFw21befCP83nlFXjyybAlizwF9Y5Zt87msS9dCo89BvMfeAWA2V2PYfr0csfDQiM0hw1VvVdVS1W1tKSkJCwx6sWaNfC//wv77w933GHRoeZ2HMhtv5zDR+zFIYfAzjuHLWV+kA/6kszmzXD++fD44zC727H8+/oZvDCnK82awSefhC1dfpAvOnPggfZuOeEEuPFGmNTzF3RnLm8u2pU1awp3qlc2jNcSYNekz52CvLzjiitsQdO77rLQdP/7v/DfaS248ame7NKlqQfaqD0FozMJGjeGY46x/Tdnt+VXd+/JlGmNOPJIeO45HwKrgYLTl549LX7zZZfBrNnC1aO6066dlc2aFa5sYZEN4/Uu0ENEuohIY2AY8GwWrhMqv/kN3H239UUPHQq77gp//rOVff65lbVoEa6MMaIgdCaVm26y9Be/gDZt4Cc/MSfVhQsLOt5qbSg4fenVy1rkCxfCmWfC2WdX7C4sxEUJMm68VHUrcBEwAZgFPK6qeRe4beBAS3v3hpNOgqOOqth8v+eegp78nhaFojOp9OljyzK9846tbjF6tLXkAQYMCFe2KFOI+tKunY199eoFTzxhvT4/+AG0bWvlr7wSrnxhkJUxL1V9XlV7qmo3Vb0uG9cIm2OOsW3NGus2fOAB83hu0ACuvNIG3c84o/xl5FRPIehMZZxxBkydak4/118P69db/sEHhytX1Ck0fUkEYOnRA+bMgV//GiZOhLPOsh6eCRPClS8MPMJGPRg2zPqh33rLjNTf/25dP9ddB//zPzYY362bj184VfPjH9s0i4EDYfZsGztt374wX0ZO1dxxB3ToYOvd7rabBZPfuhXOPRcOO6ww9cWNVz04/nho2tTmef31rxYR6pJLzHHj1lttPGzdOp+341TOxo1W2UlUbnbayULVtWpluuM4CVq0sNbW669bWMPRo6F/f+jb1zya58+He+8NW8rc4sarHrRqZV2Hjz1mNaP99zd31r/8xWpE06bZSyjhyOE4CT74wCadvvpq+eoWLVpAy5Y2ruE4qZxzjnmpXncdvPeejbUPHw733WflhTan1I1XPRk2zLoM166Ft9+2WtCUKVYL2mcf8wy67z5YsSJsSZ2ocM01VmtOrIKy995WAfrkExg0CD7+OEzpnKhSXAynnGIOGw0a2JyvsWOheXM48kj43e/CljC3uPGqJ32SIqo1aWLjF3fcAX/6kzlt/OhH1hWUcIt2Cpcvv7Sa8x/+YF2FJSUWmeWcc+Crryxay7ffWhfQ66+HLa0TRUpLbYL7tm3mbXjTTebkc+mlYUuWe3wl5Xqyxx7QqZOtPzl8OHz6KUyeDA89VPG4W24xT6Gf/SwUMZ0IsHlz+fQJEWux339/5cf+6lflq144Dti75S9/sf0mTeDZZ63F1bUrDB4crmxh4C2veiJiq1momnGaOBE++8y8fw46qOKxiRnxTmGy885W0TnjDKs5b9oETz8NRx9tepTMxReHIqITUWbPhkMOMQewW24x3bn+epti8fOfF6aDjxuvDNC3Lxx6qHUXPvigOW4cfbS50BcVmRfZoEFw8slhS+qETdu21n3417/CnnvCiSdaC6trVyu//HJLmzQJTUQnYrz/vhmurVvhv/81r8PSUnMEa9rUnMMKETdeGWLAAGtxDR9eHli1f3/rQlyxAo44IlTxnIgwb5611P/nf6xS849/WNfz/PkWlSXR4lq7Nlw5nfBRNZ0YMMD2zzzTws4deyyUldlYV1FReZSNQsPHvDLApk0WTFXE3J3XrbNB+SuvhH/9y445/PBwZXSiQZMmtuL2O+9YpPChQy2O4YMPwumnl8eoW7cuXDmd8Jk8GS680PbXrYObb7YYmM2alXuqrl9vrbFDDw1LyvDwllcGuOgia9qr2rjGW2+Z22pREbz8ss2K79cvbCmdKHDbbeZR+MUXVqGZNs0qOKefbuWJYM7e8nIOO8zG0QHuvNOcwo46ygJ/n3yyBT/o2dPc5z/7LFxZw8CNVz259lrr+gGLDv7++zbmlWDSJAv9U4gDqs72HHWUpaefboPwzz5rk00TNGhg83a85eUAfPgh7LUXXH21DT2MG2fzux5/HDp2hGeesZ6fH/7QIrYUEm686sH771sLq3Fj+Pe/rXbUvHl5+aef2uZdhk6CkhLYYQdz2njhhcpdnFu08JaXYzRrBqNGmb7Mn286c9ll5d6pvXrBww9bxI2rrw5V1JzjY171oF8/6zK88ELzHEulpMS6hA44IPeyOdFExFycd9ih6nGKe+6xWrXjgDl+XXmlVXQOOWT78mOPtXmlgwblXrYwceNVT26/veqyFi3cPd7ZnprmcJ1wQm7kcOLDdTUs+vKTn+RGjighGoH1OkRkJbAwB5cqBlb5ddJmd1UtycF1akWgL+vJzb3nilw9y1wQKX0B15kYkLbORMJ45QoRmaqqpTUf6deJOvl27/l2P1Ek337jfLufdHGHDcdxHCd2uPFyHMdxYkehGa9crTWab9eJIvl27/l2P1Ek337jfLuftCioMS/HcRwnPyi0lpfjOI6TB7jxchzHcWJHXhovETlFRGaKyDYRKU0pu0JEykRkjogcnZQ/OMgrE5ERdbxuvc+RdK4HRGSFiMxIymsrIhNFZG6QtgnyRURuC647XUT61+faUUNE/iwis4N7e0pEWieVZe155oK4yBkn8llfIF6yZhVVzbsN2APoBUwGSpPy+wAfAk2ALsA8oGGwzQO6Ao2DY/qkec16nyPlfAOA/sCMpLybgBHB/gjgxmB/KDAeEOBA4O2wn0GGn+cgoCjYvzHpvrP2PHN0X7GQM25bvuqL60zFLS9bXqo6S1XnVFJ0AjBWVTep6qdAGbB/sJWp6nxV3QyMDY5Nh0ycI/keXgW+rET+0cH+aODEpPwxakwBWotIh7peO2qo6ouqujX4OAXoFOxn83nmgrjIGSvyWF8gXrJmlbw0XtXQEViU9HlxkFdVfibOnUnaq+rSYH8Z0D6H144K52KtTMju88wFcZEzzuSTvkC8ZM0qsQ3MKyIvATtXUnSVqj6Ta3lyjaqqiOTNPIfaPE8RuQrYCjycS9mc6OH64sTWeKnqkXX42hJg16TPnYI8qsnPxLkzxXIR6aCqS4NuwRU5vHZWqel5isjZwLHAERp0/pPd55kLYv/cwqJA9QVcZ8oJe9AtmxvbO2zsScUB2/nYAGhRsN+F8kHQPdO8Vr3PUck5O1PRYePPVHTYuCnYP4aKDhvvhP3bZ/g5DgY+BkpS8rP2PHN0X7GQM25bvuqL60zKbxG2AFl6wCdhfcGbgOXAhKSyqzBvnTnAkKT8ocAnQdlVdbxuvc+RdK5HgaXAluBezgPaAZOAucBLQNvgWAHuDK77UbLBzocNG1hfBHwQbHfn4nnm6N5iIWectnzWl7jJms3Nw0M5juM4saPQvA0dx3GcPMCNl+M4jhM73Hg5juM4scONl+M4jhM73Hg5juM4scONl+M4jhM73Hg5juM4seP/Aa3D2P+aw7y2AAAAAElFTkSuQmCC\n"},"metadata":{"needs_background":"light"}}]},{"cell_type":"markdown","source":"## Step 6: calculate modal quantities\n\nNow we want to compute the mass and stiffness coefficients that go on the diagonal of the respective mode, $m_{ii}$ and $k_{ii}$. Here, we also compute the contribution from the external force, inluding the effect of the boundary conditions, $f_i^{eq}$. \n\n$$ \\begin{bmatrix} m_{11} \\ddot{\\Gamma}_1 \\\\ m_{22} \\ddot{\\Gamma}_2 \\\\ \\vdots \\\\ m_{NN} \\ddot{\\Gamma}_N \\end{bmatrix} + \\begin{bmatrix} k_{11} \\Gamma_1 \\\\ k_{22} \\Gamma_2 \\\\ \\vdots \\\\ k_{NN} \\Gamma_N \\end{bmatrix} = \\begin{bmatrix} f_1^{eq} \\\\ f_2^{eq} \\\\ \\vdots \\\\ f_N^{eq} \\end{bmatrix} \\rightarrow \\begin{cases} m_{11}\\ddot{\\Gamma}_1 + k_{11}\\Gamma_1 = f_1^{eq} \\\\ m_{22}\\ddot{\\Gamma}_2 + k_{22}\\Gamma_2 = f_2^{eq} \\\\ \\qquad \\qquad \\vdots \\\\ m_{NN}\\ddot{\\Gamma}_N + k_{NN}\\Gamma_N = f_N^{eq} \\end{cases} $$ \n\nAdditionally, we also add damping into the system. In that case we use the modal damping ratio approach, with a damping ratio of 1% ($\\xi_j=0.01$).\n\n$$ c_{jj} = \\xi_j c_{cr,j} = 2 \\xi_j \\sqrt{m_{jj}k_{jj}} $$\n\nThis will lead to a final system:\n\n$$ \\begin{cases} m_{11}\\ddot{\\Gamma}_1 + c_{11}\\dot{\\Gamma}_1 + k_{11}\\Gamma_1 = f_1^{eq} \\\\ m_{22}\\ddot{\\Gamma}_2 + c_{22}\\dot{\\Gamma}_2 + k_{22}\\Gamma_2 = f_2^{eq} \\\\ \\qquad \\qquad \\vdots \\\\ m_{NN}\\ddot{\\Gamma}_N + c_{NN}\\dot{\\Gamma}_N + k_{NN}\\Gamma_N = f_N^{eq} \\end{cases} $$ \n\n**Attention**: we are calculating all the modes, but in general this shouldn't be the case. You should only calculate the response for these modes that you want to consider. Here, we are calculating all modes, because we will compare the \"reduced\" and the \"full\" response. To consider only a few, you should replace `nMode` by the number of modes you wish to consider.","metadata":{}},{"cell_type":"code","source":"nMode = len(f)\nPHI = vr[:,0:nMode+1]\nMm = np.zeros(nMode)\nKm = np.zeros(nMode)\nCm = np.zeros(nMode)\nModalDampRatio = 0.01\n# Compute your \"nMode\" entries of the modal mass, stiffness and damping","metadata":{"trusted":true},"execution_count":24,"outputs":[]},{"cell_type":"code","source":"# Compute the equivalent force (accounting for boundary conditions)","metadata":{"trusted":true},"execution_count":26,"outputs":[]},{"cell_type":"code","source":"# Solve the resulting ODE","metadata":{"trusted":true},"execution_count":27,"outputs":[]},{"cell_type":"code","source":"# Show the result on the original coordinates","metadata":{"trusted":true},"execution_count":28,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}