{
"cells": [
{
"cell_type": "raw",
"metadata": {
"pycharm": {
"name": "#%% raw\n"
}
},
"source": [
"Code provided under BSD 3-Clause license, all other content under a Creative Commons Attribution license, CC-BY 4.0. (c) 2022 Francesco Mario Antonio Mitrotta."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"# Buckling Analysis of Euler's Column\n",
"\n",
"***\n",
"\n",
"This notebook shows the calculation of the buckling load of Euler's column, that is to say a pin-ended column loaded in compression. Four methods are considered: Euler's analytical formula, MSC Nastran SOL 105's linear buckling analysis, MSC Nastran SOL 106's nonlinear buckling method and monitoring the definiteness of the tangent stiffness matrix during a nonlinear static analysis (using MSC Nastran SOL 106). The notebook is freely inspired by [this article](https://simcompanion.hexagon.com/customers/s/article/buckling-analysis-of-column-kb8021539) on the buckling analysis of a column with MSC Nastran.\n",
"* [Euler's buckling load](#euler)\n",
"* [Setup of the numerical model](#numerical-model)\n",
"* [SOL 105 - linear buckling analysis](#linear-buckling)\n",
"* [SOL 106 - nonlinear buckling method](#nonlinear-buckling)\n",
"* [SOL 106 - tangent stiffness matrix](#tangent-stiffness-matrix)\n",
"* [Conclusion](#conclusion)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Euler's buckling load \n",
"\n",
"***\n",
"\n",
"Let's consider a column with diameter $d=20$ mm, length $l=420$ mm and Young's modulus $E=207$ GPa. The column is pinned at one end and has a roller support at the other end, where the load is applied in compression, as shown below."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.display import SVG # class to embed an SVG into the display\n",
"import os # module with miscellaneous operating system interfaces\n",
"\n",
"SVG(filename=os.path.join('resources', '01_EulerColumn.svg')) # sketch of Euler's column with boundary conditions and applied load"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Euler's buckling load $P_{c}$ is given by ([Megson, 2017, Chapter 8](http://www.aerostudents.com/courses/structural-analysis-and-design/Aircraft_Structures_for_Engineering_Students_6th_Edition.pdf)):\n",
"$$P_{c}=\\frac{\\pi^2EI}{l^2},$$\n",
"where the area moment of inertia $I$ for a circular cross-sectional shape is calculated as:\n",
"$$I = \\frac{\\pi d^4}{64}.$$\n",
"\n",
"Let's calculate $P_{c}$ for our column."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Euler's buckling load: 90962 N\n"
]
}
],
"source": [
"import numpy as np # package for working with arrays\n",
"\n",
"d = 20. # [mm]\n",
"l = 420. # [mm]\n",
"E = 207e3 # [MPa]\n",
"I = np.pi*d**4/64 # [mm^4]\n",
"P_c = np.pi**2*E*I/l**2 # [N]\n",
"print(f\"Euler's buckling load: {P_c:.0f} N\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Setup of the numerical model \n",
"\n",
"***\n",
"\n",
"We want to run some linear and nonlinear FE analyses with MSC Nastran, and to do this we need to create an appropriate numerical model together with its corresponding input file. To create, manipulate, run and analyze the results of Nastran models in our jupyter notebooks, we are going to use the [`pyNastran`](https://github.com/SteveDoyle2/pyNastran) package.\n",
"\n",
"We are going to define a base input object for our Nastran analyses. We'll complete this input object each time depending on the analysis we want to run. We'll use the function `create_base_bdf` from the module `column_utils`, where we can find some useful functions to work with a Nastran model of Euler's column. `create_base_bdf` returns a `BDF` object representing a Nastran input file (with file extension `.bdf`) including nodes, beam elements, material properties, boundary conditions, unit compression force applied to the column and some parameters for the output format of the results. The function takes as input:\n",
"- material properties including Young's modulus, Poisson's ratio and density;\n",
"- geometric properties of the column, namely length and diameter;\n",
"- number of beam elements used to discretize the column.\n",
"\n",
"Let's define the missing material properties and let's create a base bdf input with 420 beam elements, so that each element is 1 mm long. Note that we'll be working with the _millimeter - newton - megapascal - ton_ system of units.\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"subcase=0 already exists...skipping\n"
]
}
],
"source": [
"from resources import column_utils # module with functions to work with a Nastran model of Euler's column\n",
"\n",
"# Define missing parameters for function create_base_bdf\n",
"nu = 0.3 # Poisson's ratio\n",
"rho = 7.8e-4 # material density [tons/mm^3]\n",
"no_elements = 420 # number of beam elements used to discretize the column\n",
"\n",
"# Create base bdf input\n",
"base_bdf_input = column_utils.create_base_bdf(young_modulus=E, poisson_ratio=nu, density=rho, diameter=d, length=l, no_elements=no_elements)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"Now let's get a summary of our `BDF` object with the `get_bdf_stats` method. The output of this method is a text with a list of the bulk data cards of our Nastran input, that we can print to the screen."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"---BDF Statistics---\n",
"SOL None\n",
"\n",
"bdf.loads[4]: 1\n",
" FORCE: 1\n",
"\n",
"bdf.spcadds[3]: 1\n",
" SPCADD: 1\n",
"\n",
"bdf.spcs[1]: 1\n",
" SPC1: 1\n",
"\n",
"bdf.spcs[2]: 1\n",
" SPC1: 1\n",
"\n",
"bdf.params: 0\n",
" PARAM : 1\n",
"\n",
"bdf.nodes: 0\n",
" GRID : 421\n",
"\n",
"bdf.elements: 0\n",
" CBEAM : 420\n",
"\n",
"bdf.properties: 0\n",
" PBEAML : 1\n",
"\n",
"bdf.materials: 0\n",
" MAT1 : 1\n",
"\n",
"\n"
]
}
],
"source": [
"print(base_bdf_input.get_bdf_stats())"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Note that no solution sequence is defined in the function `create_base_bdf`. We'll need to define the solution sequence for each analysis that we want to run."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## SOL 105 - linear buckling analysis \n",
"\n",
"***\n",
"\n",
"The first numerical analysis that we are going to run is a linear buckling analysis, using solution sequence SOL 105. To complete our bdf input for such analysis, we execute the following steps:\n",
"- create a deep copy of the base `BDF` object (in order to preserve the original `BDF` object and reuse it later);\n",
"- define SOL 105 as solution sequence;\n",
"- create first subcase where we apply the compression load. To accomplish this we use the function `create_static_load_subcase` from the `pynastran_utils` module, which contains many useful functions to work with `pyNastran` objects;\n",
"- create second subcase where we solve the eigenvalue problem and consequently calculate the buckling load."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"from resources import pynastran_utils # module with useful functions to work with pyNastran objects\n",
"\n",
"# Create a deep copy of the original BDF object and define solution sequence\n",
"linear_buckling_bdf = base_bdf_input.__deepcopy__({})\n",
"linear_buckling_bdf.sol = 105\n",
"\n",
"# Create first subcase for the application of the compression load\n",
"load_application_subcase_id = 1 # define subcase id\n",
"force_set_id = list(linear_buckling_bdf.loads.keys())[0] # retrieve id of the FORCE card representing the compression force\n",
"pynastran_utils.create_static_load_subcase(bdf=linear_buckling_bdf, subcase_id=load_application_subcase_id,\n",
" load_set_id=force_set_id) # create subcase\n",
"\n",
"# Create second subcase for the execution of the eigenvalue analysis\n",
"eigrl_set_id = force_set_id + 1 # define id of the EIGRL card\n",
"linear_buckling_bdf.add_eigrl(sid=eigrl_set_id, v1=0., nd=1) # add EIGRL card to define a real eigenvalue analysis with the Lanczos method and calculate only the first positive eigenvalue\n",
"eigenvalue_calculation_subcase_id = 2 # define subcase id\n",
"linear_buckling_bdf.create_subcases(eigenvalue_calculation_subcase_id) # create subcase\n",
"linear_buckling_bdf.case_control_deck.subcases[eigenvalue_calculation_subcase_id].add_integer_type('METHOD', eigrl_set_id) # add EIGRL card id to case control deck"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"The model for the linear buckling analysis is ready to run now! Let's define the name of the directory where we will run all the analyses of this notebook and the name of our input bdf file. Then we can call the function `run_analysis` from the `pynastran_utils` module, which we imported earlier. `run_analysis` creates the analysis directory if it doesn't exist alreardy, writes the input `BDF` object to a bdf file and runs the analysis with Nastran."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"# Define name of analysis directory\n",
"analysis_directory_name = '01_Buckling_Analysis_of_Euler_Column'\n",
"analysis_directory_path = os.path.join(os.getcwd(), 'analyses', analysis_directory_name)\n",
"\n",
"# Define input name\n",
"base_input_name = 'euler_column'\n",
"linear_buckling_input_name = base_input_name + '_linear_buckling'\n",
"\n",
"# Run analysis\n",
"pynastran_utils.run_analysis(\n",
" directory_path=analysis_directory_path, bdf=linear_buckling_bdf, filename=linear_buckling_input_name, run_flag=False)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"Now we can read Nastran's op2 output file with the `read_op2` function and we can look at the eigenvalue calculated in the second subcase. Since the eigenvalue represents the ratio between the buckling and the applied load and since our applied load is exactly 1 N, the eigenvalue corresponds to the buckling load. Let's print its value and the percentage difference with respect to the analytical result."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SOL 105 buckling load: 90578 N\n",
"Difference with respect to analytical buckling load: -0.42 %\n"
]
}
],
"source": [
"from pyNastran.op2.op2 import read_op2 # function to read Nastran's op2 file\n",
"\n",
"# Read op2 file\n",
"op2_filename = os.path.join(analysis_directory_path, linear_buckling_input_name + '.op2')\n",
"linear_buckling_op2 = read_op2(op2_filename=op2_filename, load_geometry=True, debug=None) # load results and geometry and do not print any debug or info message\n",
"\n",
"# Find eigenvalue and print results to screen\n",
"sol_105_buckling_load = linear_buckling_op2.eigenvectors[eigenvalue_calculation_subcase_id].eigr # retrieve eigenvalue from the OP2 object\n",
"print(f'SOL 105 buckling load: {sol_105_buckling_load:.0f} N')\n",
"print(f'Difference with respect to analytical buckling load: {(sol_105_buckling_load/P_c-1)*100:.2f} %')"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We observe only a negligible difference with respect to Euler's buckling load, which can be most probably ascribable to the discretization error.\n",
"\n",
"Let's visualize the buckling mode using the function `plot_buckling_mode` from the `column_utils` module."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAq4AAAIDCAYAAADFfZw2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAABJ0AAASdAHeZh94AABX+ElEQVR4nO3de3hU5bn38d+EBAIJIYIEHBAQ2JBETloJ9RjYgILiK6RlA7WvysHtG1oU3QVULIcWa6TWigeoVZFQBCsStHJIRQRkU0oKKoIhnBSCDpAYwBwkQJL1/kEnEjMzOcxxrfl+risXstaaWU/Mw8ydZ+7nvm2GYRgCAAAAQlxEsAcAAAAA1AeBKwAAAEyBwBUAAACmQOAKAAAAUyBwBQAAgCkQuAIAAMAUCFwBAABgCgSuAAAAMAUCVwAAAJhCZLAHYFZnzpzRli1bdOWVV6pZs2bBHg4AAIDpnDt3TseOHVNqaqri4+PrvJ7AtZG2bNmikSNHBnsYAAAApvfOO+/orrvuqvM6AtdGuvLKKyVd/B/dvXt3v96rtLRUOTk5SklJUWxsrF/vhdDEHABzAMwBWHEOHDp0SCNHjqyOq+pC4NpIzvSA7t276+qrr/brvYqLi3XixAklJSUpLi7Or/dCaGIOgDkA5gCsPAfqm3bJ5iwAAACYAoErAAAATIHAFQAAAKZA4AoAAABTIHAFAACAKRC4AgAAwBQIXAEAAGAKBK4AAAAwBQJXAAAAmAKBKwAAAEyBwBUAAACmQOAKAAAAUzBV4FpSUqLp06fr1ltvVdu2bWWz2TRnzpx6P76goED33XefLr/8crVo0ULXX3+9Nm7c6L8BAwAAwGdMFbgWFRXpz3/+s86dO6eRI0c26LHnzp3T4MGDtXHjRi1YsEDvvvuu2rVrp2HDhmnLli3+GTAAAAB8JjLYA2iIzp076/Tp07LZbPrmm2/06quv1vuxr732mvbu3at//OMfuv766yVJgwYNUt++fTV9+nTt2LHDX8MGAL/KdRQrIztPO4+c0vmKKkVG2NS8aROdPV+piiqj1t+bRUboui6tNWNYopLtccEePgDUm6kCV5vN1ujHrl69Wj179qwOWiUpMjJSP//5z/X444/r66+/VocOHXwxTADwqVxHseat2auPjzaRbec2RUd9H4RG2KTzlUaN6yuqDJVXVLn/+/lKbTlQqC0HCtUs0qbKKlUHt+crqhQRYdO1nS4jsAUQckwVuHpj7969uvnmm2sd79OnjyTp888/dxu4FhQUqLCwsMaxQ4cOSZJKS0tVXFzs49HWVFZWVuNPhB/mQHjZf7JUr23/SrnHS/Tt2Qp9W17x7zM2qbJKZy9UeXx8Q5yruBj0/jC4dQa28c0jFRcdqeQrWmri9R3Vs12sz+6NhuF1AFacA6WlpQ26PmwC16KiIrVu3brWceexoqIit49duHCh5s6d6/JcTk6OTpw44ZtB1iEnJycg90HoYg5Y09dl0gdfRyi/zKbvLkjfVUqSTZLx7z+D58zZCp05e0H5p8uVnVugmEipRaR0ZYyhIR2q1CEmqMMLS7wOwEpzID8/v0HXh03gKnlONfB0bvLkyRo9enSNY4cOHdLIkSOVkpKipKQkn43RlbKyMuXk5CglJUUxMbxLhCPmgPXsP1mqP246ol35Z6pXPWsLbtD6PVv1n2UVUlmFVFhu08dFEazIBhCvA7DiHNi3b1+Drg+bwLVNmzYuV1VPnTolSS5XY50SEhKUkJDg8lxsbKzi4gKTAxYTExOweyE0MQfM7dJNVN+drwz2cHzi4opsxb9XZAv1o07x+u3I3uTG+hGvA7DSHIiNbdgvu2ETuPbu3Vt79uypddx5rFevXoEeEoAwkOso1qIth7XryCk5vi0PyD2bRdrUomlkraoCZecqam3k8rVd+Wd0+/NbZW8VrR91aa301G4EsQB8JmwC11GjRmny5MnasWOHBgwYIEmqqKjQsmXLNGDAANnt9iCPEICV5DqK9cS7e/Xx0dM+f+6W0ZG6UFFVHZS2bB6l2GaR6tWhVZ2BojOQ3vPVGZWdr1TJ2QvVzxMVGaGS6o1g3nF8Wy7Hbofe2+1gFRaAz5gucF2/fr3KyspUUlIiScrNzdXbb78tSbr99tvVokULTZw4UZmZmTp8+LA6d+4sSZowYYJeeukljR49WhkZGUpISNDChQu1f/9+ffDBB0H7fgBYz5rdDk1Z8Yl8tbbZPCpCkRE2dYi+oN/+9Efq/x+N/0U72R6nF8Zd4/a8M53hk/zTulBRVaPSQGPtyj+jO57fqlt6tKXEFgCvmC5wTU9P19GjR6v/vnLlSq1cuVKS9OWXX6pLly6qrKxUZWWlDOP7t41mzZpp48aNmj59uqZMmaLvvvtO/fr10/r165Wamhrw7wOAtfgyf7V5VISimkTomktqqRYXF2vTpk1+3/yUbI/T0gkp1X/3VSBr6PsSW6zAAmgs0wWuR44cqfOaJUuWaMmSJbWOt2vXTpmZmb4fFICw5Kv81SYRNkWHaDcrfwSy5MECaCzTBa4AEArW7HbowTc/UZUX+QBmXHl0Fcg2NpfXmQe79jOHnh97jUb0Za8BAM8IXAGggdbsduiXKz5p1GM7xEfr2s7WWWVMtscpK/0Gr1afqwxV//8keAXgCYErANSDt2kBZlxdbYhLN301Nt/3lys+0YxVn4VkygSA0EDgCgB1aGxagE0Ky530l6YTNLTCQtn5Sm05UKitBwtJHwBQC4ErAHjQ2LQAq6+w1teIvnZ1bRvb4DxY0gcAuELgCgA/0NiPuq2Wv+or3uTBkj4A4FIErgBwicY2D3hxHB9r1+XSPNiGpF+QPgDAKSLYAwCAUOFMC2hI0BphI2htjBF97Voz5Wbd0qNtvR/jTB+4d3GOch3FfhwdgFDFiiuAsOdMDfjoQGG9H0NagPecm7gaWguW1VcgfBG4AghrDa0YEGETAZOPOXNgG/KzYPMWEJ4IXAGErYZUDAjX0laB5KxA0JDVb4JXILwQuAIIOw1NDaC0VeA0Jn3glys+0cpdX/FLBRAGCFwBhJWGfBxNWkDwXFpCqz6/ZJD3CoQHqgoACBvO1ID6BK2pPdpqzZSbCYKCzLkC++K/y2h54sx7XbPbEYCRAQgGVlwBWF5DUwMobxV6nD+P+uQkkzoAWBcrrgAsbc1uh0a8sLVeQSs1WUPbiL52vTjuGkXY6r52y4FCjXhhK6uvgMWw4grAshpSNSCVigGm0JDKA5TMAqyHwBWAJTUkaGWV1Vycea/1/RkTvALWQaoAAEvJdRTrnsU59QpoSA0wt4akDtAqFrAGVlwBWEZDSl2RGmANDUkdoGQWYH6suAKwhFxHcb2D1hfHXaPMCSkErRbR0JJZD775CSuvgEkRuAKwhKez8+odtLLaZk3O1IG6VBlSRnZeAEYEwNcIXAGYmjOndUsdHxOTzxoe6pv3+tGBQnJeARMixxWAadU3p5V81vBS37xXcl4B8yFwBWBK9S2FxCpreKpvySxqvQLmQqoAANOpb9B6S4+2BCNhbkRfu1J7tK3zul+u+IQuW4AJELgCMBVn9YC6REdF6NFhiQEYEULdjGGJio6q++2OagNA6CNwBWAauY5iTVq6s145rVnpN5LTCkkX0way0m/ULXWsvFYZ0qSlOwlegRBG4ArAFNbsdmjEC1vlOHPW43XUaIUr9a316jhzViNe2EraABCi2JwFIOSxEQu+4pwfbNgCzIkVVwAhrT45rdRoRUPUt9YrOa9A6CFwBRCy6pPTao9vrjVTbiZoRYOM6GvXmik3yx7f3O015LwCoYfAFUBIqk9Oa4RNevWe68hnRaMk2+P06j3XeVx5JecVCC0ErgBCjjOnta7qAc+PvYagFV5Jtsfp+bGeN2w5c14JXoHgI3AFEFLIaUWgkfMKmAeBK4CQ8nR2HjmtCLj65rxmZOcFcFQAfojAFUBIyHUU657FOdpyoNDtNeS0wp/qk/P60YFC3bs4h5VXIEgIXAEEnXMj1kceglaJnFb4X31yXrccKGTDFhAkBK4AgsqZ01rXRixyWhEozpxXT6oMcl6BYCBwBRBUdeW0StItPdoStCKgRvS1K7VHW4/XkPMKBB6BK4CgqE9OqyRFR0Xo0WGJARoV8L0ZwxIVHeX5bZKcVyCwCFwBBFyuo1hpi7bVmdOa2qOtstJvJK8VQZFsj1NW+o26pY6V1y0HCpW2aBvBKxAABK4AAu7p7DyVX6jyeM2L465R5oQUglYEVbI9TksnpNSZ81p+oYq0ASAACFwBBNSa3Y460wPIaUWoqU/O60cHCqk0APgZgSuAgKlPVyxyWhGq6pPzSqUBwL8IXAEETF0VBMhpRSirT85rlSFNWrqT4BXwEwJXAH5XnwoCt/RoS04rQp4z59VT2oDjzFkaFAB+QuAKwK/qU0EgwibSA2AqM4YlemwNS4MCwD8IXAH4VX0qCNDKFWbjbA1bV/BKpQHAtwhcAfhNrqO4zlqtVBCAWY3oa9eaKTfLHt/c7TVbDxSy6gr4EIErAL95OjtPnrq5UkEAZpdsj9Or91znduXVEKuugC8RuALwufpsxqKCAKzCmTbgDm1hAd8hcAXgU/XZjEUFAVhNXQ0KaAsL+AaBKwCfqmszlk1UEIA1zRiWKA97tWgLC/gAgSsAn6lPO9ebe7RlpRWWlGyP89icQKItLOAtAlcAPkE7V4C2sIC/EbgC8AnauQK0hQX8jcAVgNfqqtfKZiyEk/q2hWWzFtBwBK4AvJLrKNakzH+5rddKO1eEq7rawrJZC2g4AlcAjbZmt0MjXtgqx7flbq+hnSvCVX3awrJZC2gYAlcAjeLcjOUpr5V2rgh39WkLy2YtoP4IXAE0Sl2bsaggAFzkbAvrbuG1yqAtLFBfBK4AGqyuzVj2+OZUEAAuUVeN160HCll1BeqBwBVAg+Q6ijVp6U6Pm7Fevec6glbgBzxt1jJEiSygPghcAdRbrqNYaYu2yXHmrNtr2IwFuObcrOUOJbKAuhG4Aqi3p7PzVH6hyu15NmMBno3oa/dY35USWYBnBK4A6mXNboe2eMhrZTMWUD91tYWlRBbgHoErgDo5S1+5w2YsoP6cbWEpkQU0nKkC19LSUk2dOlV2u13R0dHq16+f3nzzzXo9dtOmTRo6dKgSEhIUGxurPn366Pnnn1dlZaWfRw2Yn6fSVzaxGQtoqPqUyGKzFlCbqQLXtLQ0ZWZmavbs2Vq/fr369++vcePGafny5R4f98EHH2jIkCGqqKjQK6+8onfeeUcDBw7UQw89pEceeSRAowfMqa7SVzf3aEvQCjRCXSWy2KwF1BYZ7AHU17p167RhwwYtX75c48aNkyQNGjRIR48e1bRp0zRmzBg1adLE5WOXLFmiqKgorVmzRjExMZKkIUOGaP/+/VqyZIkWLFgQsO8DMJP6lL4irxVovBnDErX1YKHbTzScm7WWTkgJ7MCAEGWaFdfVq1crNjZWo0ePrnF8/Pjxcjgc2rFjh9vHRkVFqWnTpmrevGY+UXx8vKKjo/0yXsDsKH0F+J+zRJa7+q4SzQmAS5lmxXXv3r1KSkpSZGTNIffp06f6/A033ODysf/v//0/rVixQg8++KAef/xxtWjRQu+9955Wr16tp556qs57FxQUqLCw5kelhw4dknQx77a42L8vKGVlZTX+RPgJxhx4cs0ej6Wvbuh6mW65Ktbv8x8X8TpgXbdcFau/TrhGU1bm6kTxuVrnDUnz1uzVH+68ShJzIJxZ8XWgtLS0QdebJnAtKipS165dax1v3bp19Xl3BgwYoA8//FCjR4/WSy+9JElq0qSJnnrqKf3P//xPnfdeuHCh5s6d6/JcTk6OTpw4UZ9vwWs5OTkBuQ9CV6DmwNdl0rYvmki1to4YkmyKijB0U2yhNm3aFJDx4Hu8DljXPV2k+Z+5+ncn/eOLU1q9qVAdYpgDsNYcyM/Pb9D1pglcJclmc/9Ziqdzu3bt0qhRozRgwAC9/PLLiomJ0YcffqgnnnhC5eXl+vWvf+3xvpMnT66VonDo0CGNHDlSKSkpSkpKatg30kBlZWXKyclRSkpKdY4uwkug50D6m3sknXFxxqb2cc30wuhk9WwX6/dx4Hu8DoSHbaV7tO2LMy7O2PTRt5dpXMxp5kAYs+LrwL59+xp0vWkC1zZt2rhcVT116pSk71deXfnFL36hdu3aafXq1dUbuAYNGqSIiAjNmTNHd999t8vVXKeEhAQlJCS4PBcbG6u4uMDk+MXExATsXghN/p4DuY5iZWTnuXnjvLgOtPi+FPJag4jXAWubOaK37nh+q8sNkf88WqIzZyJkTzbU/wrmQDiz0utAbGzDFkFMszmrd+/e2rdvnyoqKmoc37NnjySpV69ebh/76aef6kc/+lGtqgP9+/dXVVVVg6N9wIqcm7EofQUET10lsvK+jdDPM3ezWQthyzSB66hRo1RaWqpVq1bVOJ6ZmSm73a4BAwa4fazdbtfOnTtrNRvYvn27JKljx46+HzBgMk9n53ncjGUTpa+AQJgxLNFtYwJJOldxsUQWEI5ME7gOHz5cQ4cOVXp6ul555RVt2rRJ//3f/63s7GzNnz+/ejV14sSJioyM1NGjR6sf+/DDD2vv3r2688479e6772rDhg169NFHNX/+fA0ZMkR9+/YN1rcFhIS6mgxIrLYCgVLXqqtEiSyEL9MErpKUlZWl//t//69mzZqlYcOGaceOHVqxYoXuvvvu6msqKytVWVkpw/g+Q2jKlClatWqVSkpKNGnSJI0aNUpr1qzR7Nmz9c477wThOwFCy9PZeW6bDEhSdFQEq61AAM0YlqjoKPdv0YbEqivCkmk2Z0kXE3gXLFjgsdPVkiVLtGTJklrH09LSlJaW5sfRAeZU12prao+2mjEskdVWIICS7XHKSr9RGdl5bv99Oldd+beJcGKqFVcAvudptfWWHm2VOYEqAkAwJNvjtHRCilLdpA0YkiYt3UnKAMIKgSsQxjyttrIZCwgNnjZrOc6cVdqibQSvCBsErkCYynUUa9LSnW5XW9mMBYSGujZrlV+gygDCB4ErEIacNVsdZ866PM9qKxBaZgxLVLNI92/ZVBlAuCBwBcJQXTVbWW0FQkuyPU7L7u2ry5q6/oyEKgMIFwSuQJipq4oApa+A0NSzXazuT6x0m+/KqivCAYErEGY8VRGwxzdXVvqNrLYCIapDjHRD13iX51h1RTggcAXCSF1VBF695zqCViDETR10FauuCFsErkCYoIoAYA0928W6rTJAbVdYHYErEAaoIgBYC7VdEa4IXIEwQBUBwFqo7YpwReAKWBxVBABrmjEsUdFR1HZFeCFwBSyOKgKANSXb45SVfqPs8c1dnqfKAKyIwBWwMKoIANaWbI/Tq/dcR5UBhA0CV8CiqCIAhAdP+a5UGYDVELgCFkQVASC8UGUA4YLAFbAgqggA4YUqAwgXBK6AxVBFAAhPVBlAOCBwBSyGKgJAeKLKAMIBgStgIVQRAMIbVQZgdQSugIV4Wm0lrxUID3VVGWDVFWZG4ApYRF2rreS1AuHDU5UBVl1hZgSugAVQsxXApajtCqsicAVMjpqtAFyhtiusiMAVMDlqtgJwhdqusCICV8DEqNkKwBNqu8JqCFwBE6NmKwBPqO0KqyFwBUyKmq0A6oParrASAlfApBZtOUwVAQD1UleVgUVbDgd2QEAjEbgCJpTrKNbm/QUuz1FFAIArnqoMbN5fwKorTIHAFTAZZ/mrkvIKl+dZbQXgiqdV15LyCspjwRQIXAGT8VT+iioCADzxVGWA8lgwAwJXwET2nyx1uyGrZXQkVQQAeOSsMtAyOtLleTZqIdQRuAIm8tymL91uyBrYM4GgFUCdku1xGtgzweU5ymMh1BG4AibxdZn0jy/OuDxnk5Se2i2g4wFgXump3SiPBVMicAVM4m9HIyh/BcAn6iqPxaorQpXPA9dbbrlF//u//+vrpwXC2v6Tpcr71vX6COWvADSGp/JYrLoiVPk8cH3wwQf1wAMPaPjw4frkk098/fRA2Ml1FGvKylzJzVsMq60AGqOuVddJS3cSvCLk+Dxw/elPf6o9e/ZozJgx+ulPf6rRo0crL4+PHIDGcNZsPVF8TnKRKMBqKwBveFp1dZw5S21XhBy/5LhGRETovvvuU15engYOHKhbb71V48ePV35+vj9uB1hWzZqttd9eWG0F4A1Pq64StV0RelwXcvNCWVmZdu/erU8++USffPKJPv74YxUWFmrbtm368Y9/rDFjxujJJ59UixYtfH1rwFJyHcVua7ZKNBsA4BszhiVqx5dFbhubOPNd+SUZocDnK66tWrXS+PHj9b//+79KTEzUs88+q5MnT+rAgQM6cOCAqqqqlJaW5uvbApazaMtht1UE7PHNaTYAwCecTQns8c1dnjd08fUICAU+X3EtKipSq1atXJ6LjY3Vc889p7g43myBuuz56ozL4zZJr95zHUErAJ9Jtsfp1Xuu0x3Pb3X5C7O71yMg0Pyy4uqJzWbT1q1bfX1bwFJyHcXKP/Wdy3PktQLwB0/5rvmnvmOTFkKC14HriBEjNHv2bL3zzjv13nzVr18/b28LWFauo1iTlu5UlYtljwgbVQQA+M+MYYmKcFFmoMqgPBZCg9epAgUFBZo/f77OnTsnm82m1q1b69prr63x1a0brSiB+nCWv3K9ScLQlfHNWW0F4DfJ9jh1at1CR4pqf+LjLI9Ffj2CyevANScnR5WVldq7d6927dqlXbt2acOGDdqwYYNstou/tlVWVno9UCAc1Cx/9UM2JV3RMqDjARB+eneMdxm4St+Xx1o6ISXAowIu8kmOa5MmTdS3b19NmDBBL730kg4cOKCVK1eqU6dO+tOf/uSLWwCWV1f5q6gIQxOv7xjAEQEIR+mp3RQd5T48oB0sgskvDQgk6Sc/+YmmT5+uw4cpoQHUh6fyV+3jmunhXpXq2S42oGMCEH4oj4VQ5rfAVZJGjhyplStX+vMWgGV4Kn/1wuhkdYgJ6HAAhDFneSx37WApj4Vg8Tpw/eUvf6nXX39du3fvrpXLevbsWR0/ftzbWwCWV1f5K1ZaAQQa5bEQirzenPXWW29p4cKFstlsatq0qXr37q1rr71WHTp00LvvvquePXv6YpyApT2dnUf5KwAhZ8awRG09WFjr9anKEJu0EBQ+KYeVn5+vnTt3aufOnfrXv/6lt99+W6dOnVLXrl21dOlSX4wTsCxPm7I6tW6hZHuciotZ2QAQeJ7KYzk3aVEaC4Hkk5avnTp1UqdOnZSWllZ97MKFC4qKivLF0wOW5Ww24G5TVu+O8YEcDgDU4q48lqGLTQloQY1A8tvmLIJWwDNnswHHmbMuz9t0sSwNAARTemo3t5u0nE0JyHdFoPi1qgAA9zw3G7i4KYtVDADB5mmTlvR9UwIgEAhcgSCoq9lAdFQEm7IAhIwZwxJpSoCQQOAKBIGnZgP2+Ob0AgcQUmhKgFBB4AoEgadmA2x0ABCKaEqAUEDgCgRYXc0GCFoBhCqaEiDYGlUOKz8/v0HXd+rUqTG3ASzHWf6KZgMAzMpTUwLKY8HfGhW4dunSRTabuw8LavthK1ggHDnLX7mrJOBsNgAAocxTUwJneSzy9OEvjQpcFy9e3KDAFUDd5a9oNgDALNw1JZC+L49FO1j4Q6MC1/vuu8/HwwCsrT7lr2g2AMAs0lO7aUPuCbe/jNMOFv7C5iwgACh/BcBKKI+FYGnUiusPHTx4UC+//LL27duns2drtq+02WzauHGjL24DmBblrwBYjbM81h3Pb3X5iznlseAPXgeue/fu1Y9//GN16NBBhw4dUp8+ffTNN9/o66+/1pVXXqlu3fj4E+GN8lcArMpZHmuLi1QoZ3ksXuPgS16nCjz++OO67bbb9Pnnn8swDL322ms6duyY3nvvPZWXl2vevHm+GKckqbS0VFOnTpXdbld0dLT69eunN998s96Pf/fdd5Wamqq4uDjFxMTo6quv1p///GefjQ9w5ensPMpfAbCsGcMSFeFiv3aVIWVk5wV+QLA0rwPXjz/+WPfee68iIi4+VVXVxUTtO+64Q7/61a/02GOPeXuLamlpacrMzNTs2bO1fv169e/fX+PGjdPy5cvrfGxGRobS0tLUq1cvvfXWW/rb3/6myZMn6/z58z4bH/BDnjZlUf4KgBU4y2O54tykBfiK16kCp0+fVuvWrRUREaGoqCidPn26+tx1112n3/zmN97eQpK0bt06bdiwQcuXL9e4ceMkSYMGDdLRo0c1bdo0jRkzRk2aNHH52F27dmnmzJl66qmnNH369OrjgwcP9snYAHc8bcqi/BUAq3BXHsu5SeuFcdcEflCwJK9XXDt06KBvvvlGktS9e3d99NFH1ec+++wzxcbGensLSdLq1asVGxur0aNH1zg+fvx4ORwO7dixw+1jX3zxRTVr1kxTpkzxyViA+sh1FGtzXoHLczaJ8lcALCM9tZvcVXffvL+AVVf4jNeB60033aR//OMfkqS7775bGRkZmjRpkiZPnqzHHntMd955p9eDlC5uAktKSlJkZM1F4j59+lSfd+ejjz5SUlKSVq1apZ49e6pJkybq2LGjHn30UVIF4BfOLlkl5ypcnmdTFgArcW7ScqWkvEJpi7YRvMInvE4VmDlzphwOhyRpxowZOnHihN544w3ZbDb913/9l5555hmvBylJRUVF6tq1a63jrVu3rj7vztdff63CwkI9+OCD+u1vf6vk5GRt3LhRGRkZOnbsmN544w2P9y4oKFBhYc08xUOHDkm6uGGsuNi//xjLyspq/InQ9+SaPW4LczeLjNCUmzs2aN4wB8AcQKjPgV/e3FH//KJI5ypqv/aVX6jSvDV79aexvYIwMusI9TnQGKWlpQ263mYYhrsUvJDSo0cPdevWTevXr69x/Pjx47Lb7Xrqqaf06KOPunxs06ZNdeHCBa1YsUJjx46tPv7www/rueee08GDB9W9e3e3954zZ47mzp3r8tzzzz+vTp06NeI7glV9XSbN/6yJVOuDM0PNm0hTrq5Uh5hgjAwA/OvrMun5z5uovFJy9Ro4vQ+vf6gpPz9fDz74oPbu3aurr766zut90oAgENq0aeNyVfXUqVOSvl95dffYEydO6LbbbqtxfPjw4Xruuef08ccfewxcJ0+eXCu39tChQxo5cqRSUlKUlJTUkG+lwcrKypSTk6OUlBTFxPAvPtRNfydPkqtKAjal9myrn49oeAks5gCYAzDLHPisIk/Zua5fAz+ruEI/H0QZwMYyyxxoiH379jXoep8GroWFhbU6Z0nyyYpk7969tWLFClVUVNTIc92zZ48kqVcv9x8/9OnTRydOnKh13LnY7Czl5U5CQoISEhJcnouNjVVcXGByFWNiYgJ2LzRe3knXH+HYJD04JNGrnyFzAMwBhPoceHBIov6eW+iyokreybKQHrtZhPocaIiGbuL3enNWSUmJJk2apJiYGLVv315XXXVVrS9fGDVqlEpLS7Vq1aoaxzMzM2W32zVgwAC3j/3JT34iSbXSDNatW6eIiAj179/fJ2ME6JIFINx52qjl7KYFNJbXK65Tp07V8uXLNXHiRPXp00fNmjXzxbhqGT58uIYOHar09HQVFxere/fuWrFihbKzs7Vs2bLqGq4TJ05UZmamDh8+rM6dO0u6WDLr5Zdf1uTJk/XNN98oOTlZH3zwgV566SVNnjy5+jrAG7mOYk1aupMuWQDC3oxhidp6sLDW62GVIU1aulOv3nMdv8ijUbwOXNeuXauMjAw99NBDvhiPR1lZWZo5c6ZmzZqlU6dOKTExsdaGq8rKSlVWVurSPWdRUVHasGGDHn/8cf3ud7/TqVOndNVVVykjI0OPPPKI38cN63OWv3JXSYAuWQDCibOblqumBI4zZ5W2aJuy0m/kdREN5nWqQHl5uXr37u2LsdQpNjZWCxYs0PHjx3Xu3Dnt3r27RtAqSUuWLJFhGOrSpUuN461bt9af/vQnnThxQufPn9f+/fv1q1/9qs78VqA+ns7Ocxu0SnTJAhB+PL3ulV+oUkZ2XuAGA8vwOmq7/fbbtXXrVl+MBTClXEexPjrgagftRdFREXTJAhB20lO7KTrKfZix9UAh+a5oMK8D1yeeeEJ//etf9eyzz+rgwYM6depUrS/AyhZtOexy96wk2eOb83EYgLCUbI9TVvqNssc3d3ne0MXXT6AhvA5ce/Xqpby8PE2bNk2JiYlq27ZtrS/AyvZ8dcblcZvEBgQAYS3ZHqdX77muVisCJ3evn4A7Xm/OmjVrlmw2d1MSsDbKXwGAZ87yWFtcpFQ5y2PxWon68jpwnTNnjg+GAZjT09l5lL8CgDpQHgu+4tMt9QcOHND27dt18OBBXz4tEJI8bcqi/BUAfM9ZHssVZ3ksNmqhPnwSuK5cuVKdO3dWUlKSbrrpJiUmJqpz5856++23ffH0QEjytCmL8lcAUFNd5bHYqIX68DpwXbduncaOHatWrVopIyNDS5cu1VNPPaVWrVpp7NixtdqsAlaQ6yjW5v0FLs/ZJMpfAcAP1FUei41aqA+vc1yffPJJ3XrrrVq7dm2NYv7Tpk3T8OHDNW/ePA0fPtzb2wAho64uWWzKAoDanOWxJi3dKceZs7XOs1EL9eH1iuunn36qyZMn1+pAZbPZNHnyZO3evdvbWwAhxVOXrOioCDZlAYAbzvJYES6KEVUZopsW6uR14NqkSROdP3/e5bkLFy7QUhWW4mlDVsvoSJoNAEAdPG3UopsW6uJ1VNm/f3/Nnz9fZ8/WXPY/d+6cnnnmGQ0YMMDbWwAhw9OGrIE9EwhaAaAe3G3UopsW6uJ1juvcuXM1ePBgde3aVaNHj1b79u11/PhxZWVlqaioSB9++KEvxgmEBE9dstiQBQD1k57aTWt2O1wuBLBJC554HbjedNNNev/99/Xoo4/qpZdekmEYioiI0IABA7RixQrdcMMNvhgnEHR0yQIA36CbFhrL68BVklJTU7V9+3Z99913On36tC677DK1aOE6fwUwo1xHsSYt3UmXLADwEbppoTEaFbjm5+friiuuUFRUlPLz82ud/+abb2r8vVOnTo0bHRAC6ip/RZcsAGg45yatI0W1P8lydtNiwyt+qFGB61VXXaXt27crJSVFXbp0kc3moq7FJSorKxs1OCAULNpy2G3QKtElCwAaq3fHeJeBq/R9N60Xxl0T4FEhlDUqcF28eLG6detW/d91Ba6AmXnaKBAdFcGmLABopPTUbtqQe8Lt4gAbtfBDjQpc77333ur/vu+++3w1FiDkeNqQZY9vTg4WAHiBblpoKLoDAB4s2nLY7YYsglYA8F5d3bSo64pLeR24fvjhh1q5cmX130+ePKnbb79d7du31z333KPy8nJvbwEEjbuPqdiQBQC+46mbFukCuJTXgeusWbOUm5tb/ffp06dr69atuuGGG/T222/r97//vbe3AILCU5oAG7IAwLfcva460wUAyQeB64EDB3TttddKkioqKrR69Wo9/fTTysrK0m9+8xutWLHC60ECwfB0dp7bNAE2ZAGAb6WndnObLjBp6U6CV0jyQeBaXFys+Ph4SdKuXbtUVlam//N//o8kKSUlxWWdVyDU5TqK9ZGLji4SaQIA4A+e0gWcdV0JXuF14JqQkKCDBw9Kkj744AN17txZHTt2lCSVlJQoKirK21sAAbdoy2GXPbQl0gQAwF88vb4667oivHnd8nXYsGF6/PHH9fnnn2vJkiU1SmXl5eWpS5cu3t4CCKhcR7E25RW4PGcTaQIA4C/UdUVdvF5x/d3vfqd+/frplVde0TXXXKMnnnii+tzy5ct1ww03eHsLIGCc7V1Lz1W4PH9zj7akCQCAnzjrutrjm7s8z0YteL3ievnllys7O9vluU2bNik6OtrbWwAB83R2ntvf9KOjIvTosMQAjwgAwouzruuIF7bW2iBbZUgZ2XlaOiElOIND0Pm1AUFcXJyaNm3qz1sAPuNpQ1bL6Ehlpd/IaisABICnjVpbDxSy6hrGvF5xvVRhYaHOnq3dsq1Tp06+vA3gF542ZA3smUDQCgAB1LtjvI4U1a6lbeji6/UL464J/KAQdF4HriUlJXr44Ye1YsUKt12yKisrvb0N4Hfukv7ZkAUAgZee2k1rdjtcLiiwSSt8eR24Tp06VcuXL9fEiRPVp08fNWvWzBfjAgLKU5csNmQBQOAl2+N0S4+22uIihcu5SYvX5vDjdeC6du1aZWRk6KGHHvLFeICgWLTlsNsuWWzIAoDgmDEsUVsPFrrcpEW6QHjyenNWeXm5evfu7YuxAEGR6yjWZjd1W+mSBQDB42mT1ub9BWzSCkNeB6633367tm7d6ouxAAHnrNta4qZuK12yACC43L0Ol5RX0AY2DHkduD7xxBP661//qmeffVYHDx7UqVOnan0Boaquuq1sygKA4EpP7aboKNfhSvmFKmVk5wV4RAgmrwPXXr16KS8vT9OmTVNiYqLatm1b6wsIRdRtBYDQ5+ym1TLa9bYc6rqGF683Z82aNUs2m80XYwECirqtAGAOyfY4DeyZoPd2O2qdo65rePE6cJ0zZ44PhgEEHnVbAcA8qOsKyYctX7/99lv9/e9/1xtvvKHTp0/76mkBv6BuKwCYi7OuqyvOuq6wPp8Err/97W9lt9s1fPhw3XPPPfryyy8lSYMHD1ZGRoYvbgH41NPZedRtBQCTmTEsUREushOrDGnS0p0Er2HA68B14cKFmjt3riZOnKi1a9fKML6PBkaMGKG1a9d6ewvApzxtyqJuKwCELk91XR1nzlIeKwx4Hbi++OKLeuSRR/T888/r1ltvrXHuP/7jP3Tw4EFvbwH4lKdNWdRtBYDQ5ul1uvxClRZtORy4wSDgvA5cv/jiC912220uz7Vs2VJnzpzx9haAz+Q6irV5v+suWWzKAoDQ56muq8RGLavzOnBt1aqVTp486fLckSNHlJCQ4O0tAJ+o7pJV7rpLFpuyACD0Oeu62uObuzzPRi1r8zpwHTx4sObPn6+ysrLqYzabTRUVFVq0aJHb1Vgg0BZtOeyxSxabsgDAHJLtcXr1nuvcbtQiXcC6vK7j+pvf/Eb9+/dXcnKyRo0aJZvNphdffFGffPKJ8vPz9dZbb/linIDX3H181DI6Un/97+tZbQUAE3Fu1DpSVLu0IekC1uX1imv37t21bds2JSUlaeHChTIMQ0uXLtXll1+urVu3qlOnTr4YJ+AVT3Vb6ZIFAObkbqMW6QLW5ZM6rsnJycrOzlZJSYm++uorFRcX6/3331dSUpIvnh7wmqe6rWzIAgBzSk/t5jZdICM7L/ADgt/5rHOWJDVr1kx2u13Nm7tOmAaCgbqtAGBNnuq6bj1QyKqrBXmd4ypJ77zzjt544w0dPXpU5eXlNc7ZbDbt3r3bF7cBGoW6rQBgXb07xrvMczV08fX/hXHXBH5Q8BuvV1x///vfKy0tTR999JGioqLUpk2bGl+tW7f2xTiBRqFuKwBYW3pqN7nIFpAkbd5fwKqrxXi94rpw4UJNmDBBL7/8spo0aeKLMQE+4azb6q4EFnVbAcD8ku1xuqVHW21xkRJWUl6htEXblJV+I6/3FuH1imtRUZF+9rOfEbQi5FC3FQDCw4xhiW67adEG1lq8DlxvvPFG7du3zxdjAXzKU91WfvsGAOtwdtNqGe36g2TqulqH14Hrc889p5deekl/+9vfdP78eV+MCfAadVsBILwk2+M0sKfrNvPUdbUOnzQgGDJkiEaNGqUWLVooLi6uxlerVq18MU6gQRZtOUzdVgAIM57qupIuYA1eb86aPn26XnzxRfXr109JSUlq2rSpL8YFeMXdx0LUbQUA66INrPV5HbguWbJEM2bM0FNPPeWL8QBe85QmQN1WALA2d3VdnekCLF6Ym9epApWVlRo6dKgvxgJ4LddRrElLd5ImAABhylO6wKSlO8l1NTmvA9dbb71V//znP30xFsArzrqtjjNnXZ4nTQAArM9TG1jHmbNKW7SN4NXEvE4V+PWvf60xY8YoJiZGd9xxh8tOWXTPQiB4qtsqkSYAAOHCXbqA9H1dV1rBmpPXgWvfvn0lSY888ogeeeQRl9dUVlZ6exugTp4S76OjIkgTAIAwkZ7aTRtyT7hdzGCjlnl5HbjOmjVLNpu7LsFAYHjakGWPb65X77mONAEACBPOhgSTlu50mT7GRi3z8jpwnTNnjg+GAXjHU91WglYACD/J9ji9es91GvHC1lrvD866rqQLmI/Xm7OAUEDdVgDAD3naqEW6gDk1asU1Pz9fV1xxhaKiopSfn1/n9Z06dWrMbYB6oW4rAMAd6rpaS6MC16uuukrbt29XSkqKunTpUmeOK5uz4E+0dwUAuJOe2k1rP3OQLmARjQpcFy9erG7dulX/d6A2Z5WWluqJJ57QW2+9pVOnTikxMVGPPvqoxo4d26DneeKJJ/Tkk0/q6quv1t69e/00WgQKaQIAAHdoA2stjQpc77333ur/vu+++3w1ljqlpaXpX//6lzIyMtSjRw8tX75c48aNU1VVlX72s5/V6zk+/fRTPfPMM2rXrp2fR4tAIE0AAFAX0gWswzSbs9atW6cNGzZo4cKFeuCBBzRo0CC98sorGjp0qKZNm1avdISKigqNHz9eDzzwgBITEwMwavgT7V0BAPVBG1jraFTgGhERoSZNmtT7yxdWr16t2NhYjR49usbx8ePHy+FwaMeOHXU+R0ZGhk6dOqUnn3zSJ2NC8NDeFQBQX7SBtY5GpQr8sOnA66+/rtLSUt15551q3769jh8/rjVr1igmJkYTJkzwyUD37t2rpKQkRUbWHHKfPn2qz99www1uH5+bm6t58+YpKytLsbGxDbp3QUGBCgsLaxw7dOiQpIt5t8XF/p3sZWVlNf6E9PwHeR7buya2i/H7zyWQmANgDoA54J3EdjEe28A+/0Ge5o8M7U9jrTgHSktLG3R9owLXS5sO/OEPf1D79u31wQcf1AgIS0pKNGTIELVo4fo3nIYqKipS165dax1v3bp19Xl3qqqqNGHCBKWlpen2229v8L0XLlyouXPnujyXk5OjEydONPg5GyMnJycg9zGDXV84V/Iv/ezHkGRTVIShPpHHtWnT8SCMzL+YA2AOgDnQOH0ipY0RTXShyibn+8X3DO36osA07xtWmgP1Kat6Ka87Zy1cuFC///3va61itmzZUtOnT9evfvUrTZs2zdvbSJLH6gWezj377LM6ePCg/va3vzXqvpMnT66VonDo0CGNHDlSKSkpSkpKatTz1ldZWZlycnKUkpKimJgYv97LDPafLFXRPz9xccam9nHN9MLoZPVs17BV9VDHHABzAMwB7/XvX6opK3N1ovjcD87YVHROsif3D+n3DyvOgX379jXoeq8D16+//rrWx/fVTx4Z6bPVyDZt2rhcVT116pSk71defyg/P1+zZs1SRkaGmjZtqjNnzki6uFGrqqpKZ86cUbNmzdS8eXO3905ISFBCQoLLc7GxsYqLC0wuZUxMTMDuFcqWrnVft3XxfSmWzm1lDoA5AOZA4/WPi9Pi+2LdtoFduvOkXhhnD87gGsBKc6Ch6ZteVxVISkrSs88+qwsXLtQ4fv78ef3hD3/w2e793r17a9++faqoqKhxfM+ePZKkXr16uXzcF198obNnz+qhhx7SZZddVv21bds27du3T5dddpkee+wxn4wRgUHdVgBAY9EG1ty8XnGdN2+eRo4cqa5duyotLU3t27fXiRMnlJWVpRMnTuidd97xwTClUaNG6ZVXXtGqVas0ZsyY6uOZmZmy2+0aMGCAy8f169dPmzZtqnV86tSp+vbbb/X666+rY8eOPhkj/I+6rQAAb1HX1by8DlzvuOMOZWdna+bMmXrppZdUVVUlm82mlJQUvf766xoyZIgvxqnhw4dr6NChSk9PV3Fxsbp3764VK1YoOztby5Ytqy67NXHiRGVmZurw4cPq3Lmz4uPjNXDgwFrPFx8fr4qKCpfnELpo7woA8BZtYM3L68BVkgYPHqzBgwfru+++0+nTp3XZZZf5rJrApbKysjRz5kzNmjWruuXrihUrarR8raysVGVlpQzDRXQDU8t1FGvz/gKX50gTAADUl6c2sJv3F7DqGsJ82jmrRYsW6tChg1+CVuliAu+CBQt0/PhxnTt3Trt3764RtErSkiVLZBiGunTp4vG5Nm/erL179/plnPA9Z8OBkvIKl+dJEwAANIS7942S8goaEoQwnwauhYWFys/Pr/UFeOvpbPcNB6KjIkgTAAA0SHpqN0VHuQ6Dyi9UKSM7L8AjQn14HbiWlJRo0qRJiomJUfv27XXVVVfV+gK8keso1kcHCl2eaxkdqaz0G/lIBwDQIMn2OGWl36iW0a6zJrceKGTVNQR5neM6depULV++XBMnTlSfPn3UrFkzX4wLqLZoy2G5y1ge2DOBoBUA0CjJ9jgN7Jmg93Y7ap0zxEatUOR14Lp27VplZGTooYce8sV4gFrc1dWziUoCAADvpKd205rdDpcLJNR1DT1epwqUl5erd+/evhgLUIunuq0392jLaisAwCvJ9jjd0qOty3POuq4IHV4Hrrfffru2bt3qi7EAtXiq2/roMN90ZQMAhLcZwxIVYat93FnXFaHD61SBJ554Qj/96U/VsmVL3XnnnWrTpk2ta1q3bu3tbRCmaO8KAPA3T3VdSRcILV4Hrr169ZIkTZs2TdOmTXN5TWVlpbe3QRiivSsAIFBoA2sOXgeus2bNks3mYn0d8BLtXQEAgUIbWHPwOnCdM2eOD4YB1ER7VwBAINEG1hx82jnrwIED2r59uw4ePOjLp0WYob0rACAYaAMb+nwSuK5cuVKdO3dWUlKSbrrpJiUmJqpz5856++23ffH0CDOLthymvSsAIODqagNLhYHg8zpwXbduncaOHatWrVopIyNDS5cu1VNPPaVWrVpp7NixWr9+vS/GiTDibgcn7V0BAP5UVxtYKgwEn9c5rk8++aRuvfVWrV27VhER38fB06ZN0/DhwzVv3jwNHz7c29sgTHiqJEB7VwCAv3lqA0uFgeDzesX1008/1eTJk2sErZJks9k0efJk7d6929tbIIw8nZ1HJQEAQFClp3Zz25AgIzsv8ANCNa8D1yZNmuj8+fMuz124cKFWQAu4k+so1kcHCl2eo5IAACBQnBUGXNl6oJBNWkHkdVTZv39/zZ8/X2fPnq1x/Ny5c3rmmWc0YMAAb2+BMLFoy2G5WGyVRCUBAEBguXvfMUQb2GDyOsd17ty5Gjx4sLp27arRo0erffv2On78uLKyslRUVKQPP/zQF+OExXmq22oTaQIAgMBKT+2mNbsdLhdUqOsaPF6vuN500016//331aVLF7300kt64okntGjRInXp0kXvv/++brjhBl+MExZWV93Wm3u05cUBABBQyfY43dKjrctz1HUNHq9XXCUpNTVV27dv13fffafTp0/rsssuU4sWrnNDgB+qq27ro8MSAzwiAACkGcMStePLIpfvUc66rrSCDSyf7pxq0aKFOnToQNCKBqFuKwAgFFHXNfQ0asW1T58+9b7WZrNREgtuUbcVABDKqOsaWhoVuLZu3Vo2m4sCZ5coLS3Vrl276rwO4W3RlsPUbQUAhLT01G5a+5mj1vtVlSHSBQKsUYHr5s2b3Z6rqKjQn//8Z/3mN7+RzWbTz372s8aODWHA3ccs1G0FAIQKZ13XI0W1PyEkXSCwfJrjunLlSiUnJ2vKlCnq27evdu3apb/85S++vAUsxFOaAHVbAQChxN37kjNdAIHhk8B18+bNGjBggMaMGaO4uDi9//77+vvf/65+/fr54ulhUaQJAADMwlMbWBoSBI5XgeuePXt0++23a/DgwSoqKtLy5cu1c+dODR482Ffjg4WRJgAAMAtPbWBJFwicRgWux44d07333qtrr71Wu3bt0nPPPad9+/Zp7Nixvh4fLIo0AQCA2ZAuEHyN2pzVo0cPnT9/XsOGDdP06dPVsmVL7dmzx+311157baMHCGsiTQAAYDZUFwi+RgWu586dkyStX79e2dnZbq8zDEM2m02VlZWNGx0sKddRrE15BS7PkSYAAAhVnqoLbN5fQE3XAGhU4Pr666/7ehwIE7mOYqUt2ua2xStpAgCAUNa7Y7zLwLWkvEJpi7bR8dHPGhW43nvvvb4eB8LEoi2H3Qat0VERpAkAAEJaemo3bcg94fK9rPxCFSkDfubTOq5AXdztvGwZHclvqQCAkJdsj1NW+o1qGe167Y8KA/5F4IqA8VRJYGDPBIJWAIApJNvjNLBngstzVBjwLwJXBAyVBAAAVkFDguAgcEXA0HAAAGAVNCQIDgJXBAQNBwAAVkNDgsAjcIXf5TqKNWnpTtIEAACW4ildYNLSnQSvfkDgCr9y1m11nDnr8jxpAgAAs/KULuA4c1Zpi7YRvPoYgSv8ylPdVok0AQCAuXl6H3PWdYXvELjCrzwlqNNwAABgdump3RQd5T6cYqOWbxG4wm88bciyxzen4QAAwPScDQns8c1dnmejlm8RuMJvPNVtffWe6whaAQCWkGyP06v3XEdd1wAgcIXfULcVABAuqOsaGASu8AvqtgIAwg11Xf2PwBV+QXtXAEC4oQ2s/xG4wi9IEwAAhBvSBfyPwBU+R5oAACBckS7gXwSu8DnSBAAA4Yp0Af8icIVP5TqKtXl/gctzpAkAAKzOU7rA5v0FrLp6icAVPpPrKFbaom0qKa9weZ40AQBAOHD3fldSXqG0RdsIXr1A4AqfWbTlsMovVLk8R3tXAEC48NQGtvxCFSkDXiBwhc+42zHZMjqS9q4AgLDhbAPbMjrS5XkqDDQegSt8wlMlgYE9EwhaAQBhJdkep4E9E1yeo8JA4xG4wieoJAAAQE1UGPA9Alf4BA0HAACoiYYEvkfgCq/RcAAAANdoSOBbBK7wGmkCAAC4RrqAbxG4wis0HAAAwD0aEvgWgSsajYYDAADUjYYEvkPgikaj4QAAAHWjIYHvELii0Wg4AABA3WhI4DsErmgUGg4AAFB/NCTwDQJXNAqVBAAAaBgqDHiPwBWNQsMBAAAahoYE3iNwRYPRcAAAgMahIYF3CFzRYKQJAADQOKQLeIfAFQ1GmgAAAI1DQwLvmCpwLS0t1dSpU2W32xUdHa1+/frpzTffrPNxWVlZGjdunLp3767mzZurS5cuuvvuu3Xw4MEAjNpaSBMAAMA7NCRoPFMFrmlpacrMzNTs2bO1fv169e/fX+PGjdPy5cs9Pu7pp5/Wd999p5kzZyo7O1vz5s3TJ598omuvvVaff/55gEZvDaQJAADgHRoSNJ7rSrghaN26ddqwYYOWL1+ucePGSZIGDRqko0ePatq0aRozZoyaNGni8rHvvfeeEhJq1k77z//8T3Xp0kV//OMf9eqrr/p9/FaQ6yjWh3knXZ4jTQAAgPpxNiQY8+ftLtumU2HAPdOsuK5evVqxsbEaPXp0jePjx4+Xw+HQjh073D72h0GrJNntdnXs2FHHjh3z+VitKNdRrLRF21R2rtLledIEAACoPxoSNI5pVlz37t2rpKQkRUbWHHKfPn2qz99www31fr4vvvhCR48e1ciRI+u8tqCgQIWFhTWOHTp0SNLFvNviYv9OrrKyshp/BsPzH+Sp/EKVy3PNIiN0z3Xt/P7/IZyFwhxAcDEHwBywnnuua6e1nzlqpeBVGRffd+ePTKxx3IpzoLS0tEHXmyZwLSoqUteuXWsdb926dfX5+qqoqNDEiRMVGxurhx9+uM7rFy5cqLlz57o8l5OToxMnTtT73t7IyckJyH1c2fWFMw3j0hoehpo3kaYkn5cj919y5AZjZOElmHMAoYE5AOaAtbRp1kSF5dIP3193fVGgTZuOu3yMleZAfn5+g643TeAqSTabi8Jn9Th3KcMwNHHiRG3dulWrVq3SlVdeWedjJk+eXCtF4dChQxo5cqRSUlKUlJRUr3s3VllZmXJycpSSkqKYmBi/3suV/SdLVfTPT1ycsSm1Z1v9fESii3PwpWDPAQQfcwDMAWta/22esnMLf3DUpqJzkj25v3q2i60+asU5sG/fvgZdb5rAtU2bNi5XVU+dOiXp+5VXTwzD0KRJk7Rs2TJlZmbqrrvuqte9ExISXObJSlJsbKzi4gKzKSkmJiZg97rU0rXuKwk8OCQxKGMKV8GaAwgdzAEwB6zlwSGJen9foct0gaU7T+qFcfZaj7HSHIiNja37okuYZnNW7969tW/fPlVU1Nx9t2fPHklSr169PD7eGbS+/vrrevXVV/Xzn//cb2O1GhoOAADgH54aElBdoDbTBK6jRo1SaWmpVq1aVeN4Zmam7Ha7BgwY4PaxhmHo/vvv1+uvv66XX35Z48eP9/dwLYOGAwAA+Je791OqC9RmmlSB4cOHa+jQoUpPT1dxcbG6d++uFStWKDs7W8uWLauu4Tpx4kRlZmbq8OHD6ty5syTpwQcf1GuvvaYJEyaod+/e+uc//1n9vM2aNdM111wTlO/JDGg4AACAf6WndnNbXWDRlsN6YRxxipNpAlfpYuvWmTNnatasWTp16pQSExO1YsUKjR07tvqayspKVVZWyjC+/+m/9957kqTFixdr8eLFNZ6zc+fOOnLkSEDGbza5jmJtyitweY40AQAAfMOZLnCkqPYnnJv3FyjXUcx77r+ZJlVAupjAu2DBAh0/flznzp3T7t27awStkrRkyRIZhqEuXbpUHzty5IgMw3D5RdDqmrPhQOm52h09JNIEAADwJXfvqyXlFUpbtI2UgX8zVeCKwHk6233DgeioCNIEAADwofTUboqOch2WlV+oUkZ2XoBHFJoIXFFLrqNYHx34YU25i1pGRyor/UY+sgAAwIeS7XHKSr9RLaNdZ3FuPVCo/Scb1mXKighcUcuiLYflYj+WJGlgzwSCVgAA/CDZHqeBPV3XjTckvbb9q8AOKAQRuKIWd3XjbKKSAAAA/pSe2k3ueoHuO14S0LGEIgJX1OCpbuvNPdqy2goAgB8l2+N0S4+2Ls8dO1Our8sCPKAQQ+CKGjzVbX10WGLgBwQAQJiZMSxRES6WXasM6YOvwzt0C+/vHrXQ3hUAgODy1AZ23xlbWG/SInBFNdq7AgAQGty9756ttOnnmbvDtq4rgSuq0d4VAIDQ4Kmu67mKKi3acjjAIwoNBK6QdHG1dfN+2rsCABAK6qrr6i61z+oIXFHd3rWknPauAACECk91XfNPfReW6QIErtCiLYdp7woAQAhKT+3mtsJAOKYLELjC7ccNtHcFACC4PFUYCMd0AQLXMOepkgDtXQEACD53KXvhmC5A4BrmqCQAAEBoI13gewSuYY6GAwAAhDbSBb5H4BrGaDgAAIA5kC5wEYFrGCNNAAAAcyBd4CIC1zBFwwEAAMwj2R6njvHRkmqvOG3eXxA2q64ErmGIhgMAAJhP8hUtJdVedi0pr1Daom1hEbwSuIYhGg4AAGA+E6/vqKgIFzl+ksovVIVFygCBaxii4QAAAObTs12sHu5VqdhmTVyeD4cKAwSuYYaGAwAAmFeHGOmmbq1dnguHCgMErmGGSgIAAJjbxOs7hm2FAQLXMEPDAQAAzK1nu9iwbUhA4BpGaDgAAIA1hGtDAgLXMEKaAAAA1hCuDQkIXMMIaQIAAFhDsj3ObbqAlRsSELiGCdIEAACwFnfv31ZuSEDgGiZIEwAAwFrSU7spOsp1KGfVhgQErmEg11GszfsLXJ4jTQAAAHNKtscpK/1GtYyOdHneihUGCFwtLtdRrLRF21RSXuHyPGkCAACYV7I9TgN7Jrg8Z8UKAwSuFrdoy2GVX6hyeS46KoI0AQAATC6cKgwQuFqcu48JWkZHKiv9RtIEAAAwOU8VBqyWLkDgamGeKgkM7JlA0AoAgEWES0MCAlcLo5IAAADhIVzSBQhcLYyGAwAAhIdwSRcgcLUoGg4AABBewiFdgMDVokgTAAAgvIRDugCBqwXRcAAAgPDjKV1g8/4CS6y6ErhaDA0HAAAIX+7e50vKK5S2aJvpg1cCV4uh4QAAAOErPbWboqNch3flF6pMnzJA4GoxNBwAACB8JdvjlJV+o1pGR7o8b/YKAwSuFkLDAQAAkGyP08CeCS7Pmb3CAIGrhVBJAAAASNatMEDgaiE0HAAAAJJ1GxIQuFoEDQcAAMClrNiQgMDVIkgTAAAAl7JiugCBq0WQJgAAAC5lxYYEBK4WQJoAAABwxWoNCQhcLYA0AQAA4IrVGhIQuJpcrqNYm/cXuDxHmgAAAOHNag0JCFxNLNdRrLRF21RSXuHyPGkCAADASg0JCFxNbNGWwyq/UOXyXHRUBGkCAABAknUqDBC4mpi75f2W0ZHKSr+RNAEAACDJOg0JCFxNylMlgYE9EwhaAQBADVZoSEDgalJUEgAAAA1hhXQBAleTouEAAABoCCukCxC4mhANBwAAQGOYPV2AwNWESBMAAACNYfZ0AQJXk6HhAAAAaCxP6QKb9xeE/KorgauJ7D9ZSsMBAADgFXfxQkl5hdIWbQvp4JXA1URe2/4VDQcAAIBX0lO7KTrKdQhYfqEqpFMGCFxNJPd4icvjNBwAAAD1lWyPU1b6jWoZHenyfChXGCBwNYmvy6SvzpS7PEfDAQAA0BDJ9jgN7Jng8lwoVxggcDWJD76OoJIAAADwGTNWGCBwNYljZS5mlqgkAAAAGseMDQlMFbiWlpZq6tSpstvtio6OVr9+/fTmm2/W67EFBQW67777dPnll6tFixa6/vrrtXHjRj+P2Df2nyzVN66zBKgkAAAAGs1sDQlMFbimpaUpMzNTs2fP1vr169W/f3+NGzdOy5cv9/i4c+fOafDgwdq4caMWLFigd999V+3atdOwYcO0ZcuWAI2+8V7b/pUM1V5xJU0AAAB4w2zpAq63k4WgdevWacOGDVq+fLnGjRsnSRo0aJCOHj2qadOmacyYMWrSpInLx7722mvau3ev/vGPf+j666+vfmzfvn01ffp07dixI2DfR2McLiyTZEg/CF5JEwAAAN5wpgscKardSv7gSdfVjILJNCuuq1evVmxsrEaPHl3j+Pjx4+VwODwGn6tXr1bPnj2rg1ZJioyM1M9//nPl5OTo66+/9tu4vZXrKFZxeYV+GLRKpAkAAADvuYsnissrQi5dwDSB6969e5WUlKTIyJqLxH369Kk+7+mxzutcPfbzzz/34Uh9J9dRrLRF23Sy5HytczQcAAAAvuCuIYHjzNmQ66RlmlSBoqIide3atdbx1q1bV5/39FjndQ19rHRxY1dhYWGNY4cOHZJ0ccNYcbF/fqDPf5DnslNW+7hmemF0sjrGym/3RmgpKyur8SfCD3MAzAH4aw50jJX+ck9fTVmZqxPF52qcK79Qpec/yNP8kYk+vadTaWlpg643TeAqSTab65JQdZ3z9rELFy7U3LlzXZ7LycnRiRMnPD6+sXYfcebsXjo+QxEV5XLk/kuOXL/cFiEsJycn2ENAkDEHwByAv+aArcJ13PHZkQJt2nTcL/fMz89v0PWmCVzbtGnjcmX01KlTkuRyRdUXj5WkyZMn18qtPXTokEaOHKmUlBQlJSXVOf7GWP9tno7nFv7gqE19urTVoEH++c0HoamsrEw5OTlKSUlRTExMsIeDIGAOgDkAf8+BYMQd+/bta9D1pglce/furRUrVqiioqJGnuuePXskSb169fL4WOd1l6rPYyUpISFBCQmu26LFxsYqLs4/O/sfHJKozQeL/p0ucLGqQHRUhB4ckui3eyK0xcTE8LMPc8wBMAfgrzlQM+64yN9xR2xsbIOuN83mrFGjRqm0tFSrVq2qcTwzM1N2u10DBgzw+Ni8vLwalQcqKiq0bNkyDRgwQHa73W/j9kayPU5Z6TdqWHJb2VtIw5LbKiv9RkpgAQAAn3PGHXf2tSuxfUvd2dcecnGHaVZchw8frqFDhyo9PV3FxcXq3r27VqxYoezsbC1btqy6huvEiROVmZmpw4cPq3PnzpKkCRMm6KWXXtLo0aOVkZGhhIQELVy4UPv379cHH3wQzG+rTsn2OM0fmahNm45r0CBWWgEAgP8k2+P0wrhrgj0Mt0wTuEpSVlaWZs6cqVmzZunUqVNKTEzUihUrNHbs2OprKisrVVlZKcMwqo81a9ZMGzdu1PTp0zVlyhR999136tevn9avX6/U1NRgfCsAAABoIFMFrrGxsVqwYIEWLFjg9polS5ZoyZIltY63a9dOmZmZfhwdAAAA/Mk0Oa4AAAAIbwSuAAAAMAUCVwAAAJgCgSsAAABMgcAVAAAApkDgCgAAAFMgcAUAAIApELgCAADAFAhcAQAAYAoErgAAADAFAlcAAACYAoErAAAATCEy2AMwq3PnzkmSDh065Pd7lZaWKj8/X/v27VNsbKzf74fQwxwAcwDMAVhxDjjjKGdcVRcC10Y6duyYJGnkyJHBHQgAAIDJHTt2TNdee22d19kMwzACMB7LOXPmjLZs2aIrr7xSzZo18+u9Dh06pJEjR+qdd95R9+7d/XovhCbmAJgDYA7AinPg3LlzOnbsmFJTUxUfH1/n9ay4NlJ8fLzuuuuugN6ze/fuuvrqqwN6T4QW5gCYA2AOwGpzoD4rrU5szgIAAIApELgCAADAFAhcAQAAYAoEribQtm1bzZ49W23btg32UBAkzAEwB8AcAHOAqgIAAAAwCVZcAQAAYAoErgAAADAFAlcAAACYAoErAAAATIHANYSVlpZq6tSpstvtio6OVr9+/fTmm28Ge1jwUklJiaZPn65bb71Vbdu2lc1m05w5c1xe+/HHH2vIkCGKjY1VfHy80tLS9MUXX7i89oUXXlBiYqKaNWumq666SnPnztWFCxf8+J2gMT788ENNmDBBiYmJiomJUYcOHXTXXXdp165dta7l529Nn376qe644w516tRJzZs3V+vWrXX99ddr2bJlta5lDoSPV199VTabTbGxsbXOMQ8uYSBkDR061IiPjzf+9Kc/GR9++KExadIkQ5LxxhtvBHto8MKXX35ptGrVyrjllluqf6azZ8+udd2+ffuMli1bGjfffLOxdu1aY9WqVcbVV19t2O12o6CgoMa18+bNM2w2m/HYY48ZmzZtMubPn280bdrUuP/++wP0XaG+fvrTnxqDBg0yFi5caGzevNlYuXKl8eMf/9iIjIw0Nm7cWH0dP3/r2rRpk/HAAw8Yf/nLX4wPP/zQeO+994yxY8cakozf/va31dcxB8LHV199ZbRq1cqw2+1GTExMjXPMg5oIXEPU2rVrDUnG8uXLaxwfOnSoYbfbjYqKiiCNDN6qqqoyqqqqDMMwjMLCQreB6+jRo43LL7/c+Pbbb6uPHTlyxIiKijKmT59efeybb74xoqOjjf/+7/+u8fgnn3zSsNlsxueff+6fbwSNcvLkyVrHSkpKjHbt2hmDBw+uPsbPP/wMGDDAuPLKK6v/zhwIHyNGjDDuvPNO4957760VuDIPaiJVIEStXr1asbGxGj16dI3j48ePl8Ph0I4dO4I0MnjLZrPJZrN5vKaiokJr1qzRT37yE8XFxVUf79y5swYNGqTVq1dXH8vOzlZ5ebnGjx9f4znGjx8vwzD0zjvv+HT88E5CQkKtY7GxsUpOTtaxY8ck8fMPV5dffrkiIyMlMQfCybJly7RlyxYtXLiw1jnmQW0EriFq7969SkpKqn4Rc+rTp0/1eVjX4cOHdfbs2eqf96X69OmjQ4cOqby8XNL3c6F37941rrviiit0+eWXM1dM4Ntvv9XHH3+sq6++WhI//3BRVVWliooKFRYWauHChfr73/+uGTNmSGIOhIuCggJNnTpVGRkZ6tixY63zzIPaCFxDVFFRkVq3bl3ruPNYUVFRoIeEAHL+fN3NAcMwdPr06eprmzVrppiYGJfXMldC3y9+8QuVlZVp5syZkvj5h4vJkycrKipKCQkJevjhh/X888/rgQcekMQcCBeTJ09Wz549lZ6e7vI886C2yLovQbB4+ji5ro+aYQ31nQPMFfP69a9/rTfeeEMvvPCCfvSjH9U4x8/f2h5//HFNmjRJBQUFeu+99/TLX/5SZWVl+tWvflV9DXPAulatWqX33ntPn3zySZ0/I+bB9whcQ1SbNm1c/nZ06tQpSa5/+4J1tGnTRpLrlfVTp07JZrMpPj6++try8nJ99913atGiRa1rfxgMIXTMnTtX8+bN05NPPqlf/vKX1cf5+YeHTp06qVOnTpKk22+/XZL02GOP6d5772UOWFxpaal+8YtfaMqUKbLb7Tpz5owk6fz585KkM2fOKCoqinngAqkCIap3797at2+fKioqahzfs2ePJKlXr17BGBYCpFu3bmrevHn1z/tSe/bsUffu3RUdHS3p+3ymH1574sQJffPNN8yVEDV37lzNmTNHc+bM0eOPP17jHD//8JSSkqKKigp98cUXzAGL++abb3Ty5En94Q9/0GWXXVb9tWLFCpWVlemyyy7T3XffzTxwgcA1RI0aNUqlpaVatWpVjeOZmZmy2+0aMGBAkEaGQIiMjNSdd96prKwslZSUVB/Pz8/Xpk2blJaWVn1s2LBhio6O1pIlS2o8x5IlS2Sz2TRy5MgAjRr19dvf/lZz5szRE088odmzZ9c6z88/PG3atEkRERHq2rUrc8Di2rdvr02bNtX6uu222xQdHa1NmzZp3rx5zANXgleJC3UZOnSocdlllxl//vOfjQ8//NC4//77DUnGsmXLgj00eGndunXGypUrjcWLFxuSjNGjRxsrV640Vq5caZSVlRmGcbHodGxsrHHLLbcY69atM7KysoxevXp5LDr9+OOPG5s3bzZ+//vfG82aNbNU0WmreOaZZwxJxrBhw4zt27fX+nLi529d999/v/E///M/xl//+ldj8+bNxttvv22MGTPGkGRMmzat+jrmQPhxVceVeVATgWsIKykpMR588EGjffv2RtOmTY0+ffoYK1asCPaw4AOdO3c2JLn8+vLLL6uv27lzpzF48GCjRYsWRlxcnDFy5Ejj0KFDLp9zwYIFRo8ePYymTZsanTp1MmbPnm2cP38+QN8R6is1NdXtz/6Hawn8/K1p8eLFxs0332xcfvnlRmRkpBEfH2+kpqYaf/nLX2pdyxwIL64CV8NgHlzKZhiGEfBlXgAAAKCByHEFAACAKRC4AgAAwBQIXAEAAGAKBK4AAAAwBQJXAAAAmAKBKwAAAEyBwBUAAACmQOAKAAAAUyBwBQAAgCkQuAIAAMAUCFwBAABgCgSuABAkx48fV2xsrMaOHVvj+Jo1axQVFaWZM2c26PnmzJkjm82mzz77TKNHj1arVq3UunVrPfLII6qoqND+/fs1bNgwtWzZUl26dNH8+fP98hwA4C8ErgAQJFdccYWmT5+ut956S7t27ZIkbd68WaNHj1Z6erqefPLJRj3vf/3Xf6lv375atWqV7r//fv3xj3/Uww8/rJEjR+qOO+7Q6tWr9Z//+Z+aMWOGsrKy/PYcAOBzBgAgaMrKygy73W4MHjzYyMnJMVq2bGmMHz/eqKqqavBzzZ4925Bk/OEPf6hxvF+/foYkIysrq/rYhQsXjLZt2xppaWk+fw4A8BdWXAEgiFq0aKF58+Zp48aNGjRokIYPH65XXnlFNput0c85YsSIGn9PSkqSzWbT8OHDq49FRkaqe/fuOnr0qN+eAwB8jcAVAIKsR48ekiSbzaYlS5aoSZMmXj1f69ata/y9adOmatGihaKjo2sdLy8v99tzAICvEbgCQBB9+umnGjFihG688UaVlpZq8eLFwR4SAIQsAlcACJL9+/frtttu0/XXX69Nmzbprrvu0pw5c/Ttt98Ge2gAEJIIXAEgCI4cOaIhQ4aoZ8+eWrVqlaKiopSRkaHTp0/rd7/7Xa3rbTabBg4cGPiBAkAIIXAFgAA7fvy4hgwZooSEBK1Zs0bNmzeXJCUmJmrChAlasGCBjhw5Un19aWmppIvlswAgnNkMwzCCPQgAgHvr1q3TiBEjtHv3bvXu3TvYwwGAoGHFFQBC3KZNmzR27FiCVgBhjxVXAAAAmAIrrgAAADAFAlcAAACYAoErAAAATIHAFQAAAKZA4AoAAABTIHAFAACAKRC4AgAAwBQIXAEAAGAKBK4AAAAwBQJXAAAAmAKBKwAAAEyBwBUAAACm8P8BQgb5sqya8r0AAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt # plotting package\n",
"\n",
"plt.rcParams['figure.dpi'] = 120 # set default dpi of figures\n",
"column_utils.plot_buckling_mode(op2=linear_buckling_op2)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## SOL 106 - nonlinear buckling method \n",
"\n",
"***\n",
"\n",
"Our next numerical analysis is a nonlinear analysis with solution sequence SOL 106, where we want to use the so-called nonlinear buckling method. This method was described by [Lee & Herting (1985)](https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620211016) and is based on the hypothesis that the components of the tangent stiffness matrix vary linearly in a small range close to the last two converged solutions. If a critical point exists sufficiently close to the last converged solution, then the tangent stiffness matrix at such point can be expressed as:\n",
"\n",
"$$\\boldsymbol{K_{c}}=\\boldsymbol{K_n}+\\mu\\boldsymbol{\\Delta K},$$\n",
"\n",
"where $\\boldsymbol{\\Delta K}=\\boldsymbol{K_n}-\\boldsymbol{K_{n-1}}$ is the differential stiffness matrix based on the last two converged solutions $n$ and $n-1$ and $\\mu$ is a proportionality factor.\n",
"\n",
"Since the tangent stiffness matrix is singular at the critical point, we can formulate the following eigenvalue problem:\n",
"\n",
"$$\\left(\\boldsymbol{K_n}+\\mu\\boldsymbol{\\Delta K}\\right)\\boldsymbol{\\phi}=\\boldsymbol{0},$$\n",
"\n",
"where $\\mu$ and $\\boldsymbol{\\phi}$ respectively represent the eigenvalues and the eigenvectors related to the critical displacements $\\boldsymbol{U_{c}}$:\n",
"\n",
"$$\\boldsymbol{U_{c}}=\\boldsymbol{U_n}+\\mu\\boldsymbol{\\Delta U},$$\n",
"\n",
"with $\\boldsymbol{\\Delta U}=\\boldsymbol{U_n}-\\boldsymbol{U_{n-1}}$. Nastran calculates the vector of critical buckling loads $\\boldsymbol{P_{c}}$ as:\n",
"\n",
"$$\\boldsymbol{P_{c}}=\\boldsymbol{P_n}+\\alpha\\boldsymbol{\\Delta P},$$\n",
"\n",
"where $\\boldsymbol{P_n}$ is the vector of the applied loads, $\\alpha$ is the critical buckling factor and $\\boldsymbol{\\Delta P}$ is the last load increment vector applied in the nonlinear analysis. $\\alpha$ and $\\mu$ are related by the following expression:\n",
"\n",
"$$\\alpha=\\frac{\\mu\\boldsymbol{\\Delta U}^\\intercal\\left(\\boldsymbol{K_n}+\\frac{1}{2}\\mu\\boldsymbol{\\Delta K}\\right)\\boldsymbol{\\Delta U}}{\\boldsymbol{\\Delta U}^\\intercal\\boldsymbol{\\Delta P}}.$$\n",
"\n",
"As suggested by [Lee & Herting (1985)](https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620211016), we need to monitor the value of $\\alpha$ because when it is greater than unity the predicted buckling point is not close to the last solution point and as a consequence the prediction should not be considered reliable.\n",
"\n",
"To complete our bdf input for this analysis, we execute the following steps:\n",
"- create a deep copy of the original `BDF` object;\n",
"- assign SOL 106 as solution sequence;\n",
"- create first subcase where we apply a compression force equal to the buckling load predicted by SOL 105;\n",
"- define the parameters for the eigenvalue computation, this time carried out in the same subcase."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"# Create deep copy of the original BDF object and define solution sequence\n",
"nonlinear_buckling_bdf = base_bdf_input.__deepcopy__({}) # deep copy of BDF object\n",
"nonlinear_buckling_bdf.sol = 106 # solution sequence\n",
"\n",
"# Creat first subcase\n",
"nonlinear_buckling_bdf.loads[force_set_id][0].mag = sol_105_buckling_load # set force magnitude equal to SOL 105 buckling load\n",
"first_subcase_id = 1 # define subcase id\n",
"pynastran_utils.create_static_load_subcase(bdf=nonlinear_buckling_bdf, subcase_id=first_subcase_id,\n",
" load_set_id=force_set_id) # create subcase\n",
"\n",
"# Add EIGB card to define an eigenvalue analysis with inverse power method\n",
"eigb_set_id = force_set_id + 1\n",
"nonlinear_buckling_bdf.add_eigb(sid=eigb_set_id, method='INV', L1=-1e4, L2=1e4, nep=1, ndp='', ndn='', norm='', G='', C='') # define buckling analysis with inverse power method to find one positive and one negative root\n",
"nonlinear_buckling_bdf.case_control_deck.subcases[first_subcase_id].add_integer_type('METHOD', eigb_set_id) # add EIGB card id to case control deck"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"In SOL 106 we need to specify some additional parameters and cards for the nonlinear analysis.\n",
"\n",
"First we define two `PARAM` cards to enable large displacement effects and to flag the use of the nonlinear buckling method."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"PARAM BUCKLE 2"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nonlinear_buckling_bdf.add_param('LGDISP', [1]) # assume all nonlinear structural element types that have a large displacement capability to have large displacement effects\n",
"nonlinear_buckling_bdf.add_param('BUCKLE', [2]) # request nonlinear buckling method in a SOL 106 cold start run"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Successively we define the arc-length method as our iteration method. We first add the `NLPARM` card to set the general parameters of the nonlinear iteration strategy. We use the following settings:\n",
"* `ninc=100` to set the initial load increment $\\Delta\\mu^1 = 1/NINC$ to 1% of the total applied load;\n",
"* `kmethod='ITER'` and `kstep=-1` to update the stiffness matrix at every iteration. The value of -1 means that the stiffness matrix is forced to be updated between the convergence of a load increment and the start of the next load increment. We don't set `KSTEP` to 1 (corresponding to no update between load increments) because in this way the solver would ignore the displacement error in the convergence criteria (see the explanation of the NLPARM card on the *MSC Nastran Quick Reference Guide* for more details);\n",
"* `max_iter=25` to set the maximum number of iterations for each load increment;\n",
"* `conv='PUV'` to select convergence based on the load equilibrium and the displacement errors with vector component method (convergence checking is performed on the maximum vector component of all components in the model);\n",
"* `int_out='YES'` to process the output for every converged load increment;\n",
"* `eps_p=1e-3` and `eps_u=1e-3` to set the error tolerance for the load and the displacement criterion, respectively;\n",
"* `max_bisect=10` to set the maximum number of bisections allowed for each load increment."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"NLPARM 1 100 ITER -1 PUV YES\n",
" .001 .001\n",
" 10"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nlparm_id = 1 # id of the NLPARM and NLPCI cards\n",
"nonlinear_buckling_bdf.add_nlparm(nlparm_id, ninc=100, kmethod='ITER', kstep=-1, max_iter=25, conv='PUV', int_out='YES',\n",
" eps_p=1e-3, eps_u=1e-3, max_bisect=10)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Then we add a `NLPCI` card to define the parameters of the arc-length method:\n",
"* `Type='CRIS'` to set Crisfield constraint type;\n",
"* `minalr=.01` and `maxalr=1.0001` to set the minimum and maximum allowable arc-length adjustment ratio between increments, where $MINALR \\leq\\frac{\\Delta l_{new}}{\\Delta l_{old}}\\leq MAXALR$ and $\\frac{\\Delta l_{new}}{\\Delta l_{old}}$ is the adjument ratio of the arc-length;\n",
"* `desiter=5` to set the desired number of iterations for convergence, which is used to calculate the arc-length for the next increment $\\Delta l_{new}=\\Delta l_{old}\\sqrt{DESITER/I_{max}}$, where $I_{max}$ represents the number of iterations required for convergence in the previous load increment;\n",
"* `mxinc=1000` to set the maximum number of controlled increment steps allowed within a subcase."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"text/plain": [
"NLPCI 1 CRIS .01 1.0001 0. 5 1000"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"nonlinear_buckling_bdf.add_nlpci(nlparm_id, Type='CRIS', minalr=.01, maxalr=1.0001, desiter=5, mxinc=1000) # same identification number of NLPARM card is used becasue NLPCI card must be associated to a NLPARM card"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Finally, we need to select the NLPARM card that we have just defined in the case control deck."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"nonlinear_buckling_bdf.case_control_deck.subcases[0].add_integer_type('NLPARM', nlparm_id)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"The parameters employed here are inspired by the shallow cylindrical shell snap-through example from the *MSC Nastran Demontration Problems Guide - Implicit Nonlinear*.\n",
"\n",
"Now the Nastran input is ready to be run, so let's define a name for our input file and call the `run_analysis` function."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Nastran job euler_column_nonlinear_buckling_method.bdf completed\n",
"Wall time: 6.0 s\n"
]
}
],
"source": [
"nonlinear_buckling_input_name = base_input_name + '_nonlinear_buckling_method' # input file name\n",
"pynastran_utils.run_analysis(directory_path=analysis_directory_path, bdf=nonlinear_buckling_bdf,\n",
" filename=nonlinear_buckling_input_name, run_flag=False) # run analysis"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"We read the nonlinear buckling load $\\boldsymbol{P_c}$ and the critical buckling factor $\\alpha$ using the function `read_nonlinear_buckling_load_from_f06` from the `pynastran_utils` module. The nonlinear buckling load is returned by the function as an array with dimensions _(number of subcases, number of nodes, 6)_, where 6 indicates the forces and moments along the 3 axes. In this simple case the nonlinear buckling load corresponds to the $x$ component of the force at the last node of the first subcase.\n",
"\n",
"We read the op2 and f06 file, find the nonlinear buckling load and critical buckling factor, print the results and plot the buckling mode."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"SOL 106 buckling load: 90578 N\n",
"Critical buckling factor ALPHA = -8.3e-06\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAq4AAAIDCAYAAADFfZw2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAABJ0AAASdAHeZh94AABX+ElEQVR4nO3de3hU5bn38d+EBAIJIYIEHBAQ2JBETloJ9RjYgILiK6RlA7WvysHtG1oU3QVULIcWa6TWigeoVZFQBCsStHJIRQRkU0oKKoIhnBSCDpAYwBwkQJL1/kEnEjMzOcxxrfl+risXstaaWU/Mw8ydZ+7nvm2GYRgCAAAAQlxEsAcAAAAA1AeBKwAAAEyBwBUAAACmQOAKAAAAUyBwBQAAgCkQuAIAAMAUCFwBAABgCgSuAAAAMAUCVwAAAJhCZLAHYFZnzpzRli1bdOWVV6pZs2bBHg4AAIDpnDt3TseOHVNqaqri4+PrvJ7AtZG2bNmikSNHBnsYAAAApvfOO+/orrvuqvM6AtdGuvLKKyVd/B/dvXt3v96rtLRUOTk5SklJUWxsrF/vhdDEHABzAMwBWHEOHDp0SCNHjqyOq+pC4NpIzvSA7t276+qrr/brvYqLi3XixAklJSUpLi7Or/dCaGIOgDkA5gCsPAfqm3bJ5iwAAACYAoErAAAATIHAFQAAAKZA4AoAAABTIHAFAACAKRC4AgAAwBQIXAEAAGAKBK4AAAAwBQJXAAAAmAKBKwAAAEyBwBUAAACmQOAKAAAAUzBV4FpSUqLp06fr1ltvVdu2bWWz2TRnzpx6P76goED33XefLr/8crVo0ULXX3+9Nm7c6L8BAwAAwGdMFbgWFRXpz3/+s86dO6eRI0c26LHnzp3T4MGDtXHjRi1YsEDvvvuu2rVrp2HDhmnLli3+GTAAAAB8JjLYA2iIzp076/Tp07LZbPrmm2/06quv1vuxr732mvbu3at//OMfuv766yVJgwYNUt++fTV9+nTt2LHDX8MGAL/KdRQrIztPO4+c0vmKKkVG2NS8aROdPV+piiqj1t+bRUboui6tNWNYopLtccEePgDUm6kCV5vN1ujHrl69Wj179qwOWiUpMjJSP//5z/X444/r66+/VocOHXwxTADwqVxHseat2auPjzaRbec2RUd9H4RG2KTzlUaN6yuqDJVXVLn/+/lKbTlQqC0HCtUs0qbKKlUHt+crqhQRYdO1nS4jsAUQckwVuHpj7969uvnmm2sd79OnjyTp888/dxu4FhQUqLCwsMaxQ4cOSZJKS0tVXFzs49HWVFZWVuNPhB/mQHjZf7JUr23/SrnHS/Tt2Qp9W17x7zM2qbJKZy9UeXx8Q5yruBj0/jC4dQa28c0jFRcdqeQrWmri9R3Vs12sz+6NhuF1AFacA6WlpQ26PmwC16KiIrVu3brWceexoqIit49duHCh5s6d6/JcTk6OTpw44ZtB1iEnJycg90HoYg5Y09dl0gdfRyi/zKbvLkjfVUqSTZLx7z+D58zZCp05e0H5p8uVnVugmEipRaR0ZYyhIR2q1CEmqMMLS7wOwEpzID8/v0HXh03gKnlONfB0bvLkyRo9enSNY4cOHdLIkSOVkpKipKQkn43RlbKyMuXk5CglJUUxMbxLhCPmgPXsP1mqP246ol35Z6pXPWsLbtD6PVv1n2UVUlmFVFhu08dFEazIBhCvA7DiHNi3b1+Drg+bwLVNmzYuV1VPnTolSS5XY50SEhKUkJDg8lxsbKzi4gKTAxYTExOweyE0MQfM7dJNVN+drwz2cHzi4opsxb9XZAv1o07x+u3I3uTG+hGvA7DSHIiNbdgvu2ETuPbu3Vt79uypddx5rFevXoEeEoAwkOso1qIth7XryCk5vi0PyD2bRdrUomlkraoCZecqam3k8rVd+Wd0+/NbZW8VrR91aa301G4EsQB8JmwC11GjRmny5MnasWOHBgwYIEmqqKjQsmXLNGDAANnt9iCPEICV5DqK9cS7e/Xx0dM+f+6W0ZG6UFFVHZS2bB6l2GaR6tWhVZ2BojOQ3vPVGZWdr1TJ2QvVzxMVGaGS6o1g3nF8Wy7Hbofe2+1gFRaAz5gucF2/fr3KyspUUlIiScrNzdXbb78tSbr99tvVokULTZw4UZmZmTp8+LA6d+4sSZowYYJeeukljR49WhkZGUpISNDChQu1f/9+ffDBB0H7fgBYz5rdDk1Z8Yl8tbbZPCpCkRE2dYi+oN/+9Efq/x+N/0U72R6nF8Zd4/a8M53hk/zTulBRVaPSQGPtyj+jO57fqlt6tKXEFgCvmC5wTU9P19GjR6v/vnLlSq1cuVKS9OWXX6pLly6qrKxUZWWlDOP7t41mzZpp48aNmj59uqZMmaLvvvtO/fr10/r165Wamhrw7wOAtfgyf7V5VISimkTomktqqRYXF2vTpk1+3/yUbI/T0gkp1X/3VSBr6PsSW6zAAmgs0wWuR44cqfOaJUuWaMmSJbWOt2vXTpmZmb4fFICw5Kv81SYRNkWHaDcrfwSy5MECaCzTBa4AEArW7HbowTc/UZUX+QBmXHl0Fcg2NpfXmQe79jOHnh97jUb0Za8BAM8IXAGggdbsduiXKz5p1GM7xEfr2s7WWWVMtscpK/0Gr1afqwxV//8keAXgCYErANSDt2kBZlxdbYhLN301Nt/3lys+0YxVn4VkygSA0EDgCgB1aGxagE0Ky530l6YTNLTCQtn5Sm05UKitBwtJHwBQC4ErAHjQ2LQAq6+w1teIvnZ1bRvb4DxY0gcAuELgCgA/0NiPuq2Wv+or3uTBkj4A4FIErgBwicY2D3hxHB9r1+XSPNiGpF+QPgDAKSLYAwCAUOFMC2hI0BphI2htjBF97Voz5Wbd0qNtvR/jTB+4d3GOch3FfhwdgFDFiiuAsOdMDfjoQGG9H0NagPecm7gaWguW1VcgfBG4AghrDa0YEGETAZOPOXNgG/KzYPMWEJ4IXAGErYZUDAjX0laB5KxA0JDVb4JXILwQuAIIOw1NDaC0VeA0Jn3glys+0cpdX/FLBRAGCFwBhJWGfBxNWkDwXFpCqz6/ZJD3CoQHqgoACBvO1ID6BK2pPdpqzZSbCYKCzLkC++K/y2h54sx7XbPbEYCRAQgGVlwBWF5DUwMobxV6nD+P+uQkkzoAWBcrrgAsbc1uh0a8sLVeQSs1WUPbiL52vTjuGkXY6r52y4FCjXhhK6uvgMWw4grAshpSNSCVigGm0JDKA5TMAqyHwBWAJTUkaGWV1Vycea/1/RkTvALWQaoAAEvJdRTrnsU59QpoSA0wt4akDtAqFrAGVlwBWEZDSl2RGmANDUkdoGQWYH6suAKwhFxHcb2D1hfHXaPMCSkErRbR0JJZD775CSuvgEkRuAKwhKez8+odtLLaZk3O1IG6VBlSRnZeAEYEwNcIXAGYmjOndUsdHxOTzxoe6pv3+tGBQnJeARMixxWAadU3p5V81vBS37xXcl4B8yFwBWBK9S2FxCpreKpvySxqvQLmQqoAANOpb9B6S4+2BCNhbkRfu1J7tK3zul+u+IQuW4AJELgCMBVn9YC6REdF6NFhiQEYEULdjGGJio6q++2OagNA6CNwBWAauY5iTVq6s145rVnpN5LTCkkX0way0m/ULXWsvFYZ0qSlOwlegRBG4ArAFNbsdmjEC1vlOHPW43XUaIUr9a316jhzViNe2EraABCi2JwFIOSxEQu+4pwfbNgCzIkVVwAhrT45rdRoRUPUt9YrOa9A6CFwBRCy6pPTao9vrjVTbiZoRYOM6GvXmik3yx7f3O015LwCoYfAFUBIqk9Oa4RNevWe68hnRaMk2+P06j3XeVx5JecVCC0ErgBCjjOnta7qAc+PvYagFV5Jtsfp+bGeN2w5c14JXoHgI3AFEFLIaUWgkfMKmAeBK4CQ8nR2HjmtCLj65rxmZOcFcFQAfojAFUBIyHUU657FOdpyoNDtNeS0wp/qk/P60YFC3bs4h5VXIEgIXAEEnXMj1kceglaJnFb4X31yXrccKGTDFhAkBK4AgsqZ01rXRixyWhEozpxXT6oMcl6BYCBwBRBUdeW0StItPdoStCKgRvS1K7VHW4/XkPMKBB6BK4CgqE9OqyRFR0Xo0WGJARoV8L0ZwxIVHeX5bZKcVyCwCFwBBFyuo1hpi7bVmdOa2qOtstJvJK8VQZFsj1NW+o26pY6V1y0HCpW2aBvBKxAABK4AAu7p7DyVX6jyeM2L465R5oQUglYEVbI9TksnpNSZ81p+oYq0ASAACFwBBNSa3Y460wPIaUWoqU/O60cHCqk0APgZgSuAgKlPVyxyWhGq6pPzSqUBwL8IXAEETF0VBMhpRSirT85rlSFNWrqT4BXwEwJXAH5XnwoCt/RoS04rQp4z59VT2oDjzFkaFAB+QuAKwK/qU0EgwibSA2AqM4YlemwNS4MCwD8IXAH4VX0qCNDKFWbjbA1bV/BKpQHAtwhcAfhNrqO4zlqtVBCAWY3oa9eaKTfLHt/c7TVbDxSy6gr4EIErAL95OjtPnrq5UkEAZpdsj9Or91znduXVEKuugC8RuALwufpsxqKCAKzCmTbgDm1hAd8hcAXgU/XZjEUFAVhNXQ0KaAsL+AaBKwCfqmszlk1UEIA1zRiWKA97tWgLC/gAgSsAn6lPO9ebe7RlpRWWlGyP89icQKItLOAtAlcAPkE7V4C2sIC/EbgC8AnauQK0hQX8jcAVgNfqqtfKZiyEk/q2hWWzFtBwBK4AvJLrKNakzH+5rddKO1eEq7rawrJZC2g4AlcAjbZmt0MjXtgqx7flbq+hnSvCVX3awrJZC2gYAlcAjeLcjOUpr5V2rgh39WkLy2YtoP4IXAE0Sl2bsaggAFzkbAvrbuG1yqAtLFBfBK4AGqyuzVj2+OZUEAAuUVeN160HCll1BeqBwBVAg+Q6ijVp6U6Pm7Fevec6glbgBzxt1jJEiSygPghcAdRbrqNYaYu2yXHmrNtr2IwFuObcrOUOJbKAuhG4Aqi3p7PzVH6hyu15NmMBno3oa/dY35USWYBnBK4A6mXNboe2eMhrZTMWUD91tYWlRBbgHoErgDo5S1+5w2YsoP6cbWEpkQU0nKkC19LSUk2dOlV2u13R0dHq16+f3nzzzXo9dtOmTRo6dKgSEhIUGxurPn366Pnnn1dlZaWfRw2Yn6fSVzaxGQtoqPqUyGKzFlCbqQLXtLQ0ZWZmavbs2Vq/fr369++vcePGafny5R4f98EHH2jIkCGqqKjQK6+8onfeeUcDBw7UQw89pEceeSRAowfMqa7SVzf3aEvQCjRCXSWy2KwF1BYZ7AHU17p167RhwwYtX75c48aNkyQNGjRIR48e1bRp0zRmzBg1adLE5WOXLFmiqKgorVmzRjExMZKkIUOGaP/+/VqyZIkWLFgQsO8DMJP6lL4irxVovBnDErX1YKHbTzScm7WWTkgJ7MCAEGWaFdfVq1crNjZWo0ePrnF8/Pjxcjgc2rFjh9vHRkVFqWnTpmrevGY+UXx8vKKjo/0yXsDsKH0F+J+zRJa7+q4SzQmAS5lmxXXv3r1KSkpSZGTNIffp06f6/A033ODysf/v//0/rVixQg8++KAef/xxtWjRQu+9955Wr16tp556qs57FxQUqLCw5kelhw4dknQx77a42L8vKGVlZTX+RPgJxhx4cs0ej6Wvbuh6mW65Ktbv8x8X8TpgXbdcFau/TrhGU1bm6kTxuVrnDUnz1uzVH+68ShJzIJxZ8XWgtLS0QdebJnAtKipS165dax1v3bp19Xl3BgwYoA8//FCjR4/WSy+9JElq0qSJnnrqKf3P//xPnfdeuHCh5s6d6/JcTk6OTpw4UZ9vwWs5OTkBuQ9CV6DmwNdl0rYvmki1to4YkmyKijB0U2yhNm3aFJDx4Hu8DljXPV2k+Z+5+ncn/eOLU1q9qVAdYpgDsNYcyM/Pb9D1pglcJclmc/9Ziqdzu3bt0qhRozRgwAC9/PLLiomJ0YcffqgnnnhC5eXl+vWvf+3xvpMnT66VonDo0CGNHDlSKSkpSkpKatg30kBlZWXKyclRSkpKdY4uwkug50D6m3sknXFxxqb2cc30wuhk9WwX6/dx4Hu8DoSHbaV7tO2LMy7O2PTRt5dpXMxp5kAYs+LrwL59+xp0vWkC1zZt2rhcVT116pSk71deXfnFL36hdu3aafXq1dUbuAYNGqSIiAjNmTNHd999t8vVXKeEhAQlJCS4PBcbG6u4uMDk+MXExATsXghN/p4DuY5iZWTnuXnjvLgOtPi+FPJag4jXAWubOaK37nh+q8sNkf88WqIzZyJkTzbU/wrmQDiz0utAbGzDFkFMszmrd+/e2rdvnyoqKmoc37NnjySpV69ebh/76aef6kc/+lGtqgP9+/dXVVVVg6N9wIqcm7EofQUET10lsvK+jdDPM3ezWQthyzSB66hRo1RaWqpVq1bVOJ6ZmSm73a4BAwa4fazdbtfOnTtrNRvYvn27JKljx46+HzBgMk9n53ncjGUTpa+AQJgxLNFtYwJJOldxsUQWEI5ME7gOHz5cQ4cOVXp6ul555RVt2rRJ//3f/63s7GzNnz+/ejV14sSJioyM1NGjR6sf+/DDD2vv3r2688479e6772rDhg169NFHNX/+fA0ZMkR9+/YN1rcFhIS6mgxIrLYCgVLXqqtEiSyEL9MErpKUlZWl//t//69mzZqlYcOGaceOHVqxYoXuvvvu6msqKytVWVkpw/g+Q2jKlClatWqVSkpKNGnSJI0aNUpr1qzR7Nmz9c477wThOwFCy9PZeW6bDEhSdFQEq61AAM0YlqjoKPdv0YbEqivCkmk2Z0kXE3gXLFjgsdPVkiVLtGTJklrH09LSlJaW5sfRAeZU12prao+2mjEskdVWIICS7XHKSr9RGdl5bv99Oldd+beJcGKqFVcAvudptfWWHm2VOYEqAkAwJNvjtHRCilLdpA0YkiYt3UnKAMIKgSsQxjyttrIZCwgNnjZrOc6cVdqibQSvCBsErkCYynUUa9LSnW5XW9mMBYSGujZrlV+gygDCB4ErEIacNVsdZ866PM9qKxBaZgxLVLNI92/ZVBlAuCBwBcJQXTVbWW0FQkuyPU7L7u2ry5q6/oyEKgMIFwSuQJipq4oApa+A0NSzXazuT6x0m+/KqivCAYErEGY8VRGwxzdXVvqNrLYCIapDjHRD13iX51h1RTggcAXCSF1VBF695zqCViDETR10FauuCFsErkCYoIoAYA0928W6rTJAbVdYHYErEAaoIgBYC7VdEa4IXIEwQBUBwFqo7YpwReAKWBxVBABrmjEsUdFR1HZFeCFwBSyOKgKANSXb45SVfqPs8c1dnqfKAKyIwBWwMKoIANaWbI/Tq/dcR5UBhA0CV8CiqCIAhAdP+a5UGYDVELgCFkQVASC8UGUA4YLAFbAgqggA4YUqAwgXBK6AxVBFAAhPVBlAOCBwBSyGKgJAeKLKAMIBgStgIVQRAMIbVQZgdQSugIV4Wm0lrxUID3VVGWDVFWZG4ApYRF2rreS1AuHDU5UBVl1hZgSugAVQsxXApajtCqsicAVMjpqtAFyhtiusiMAVMDlqtgJwhdqusCICV8DEqNkKwBNqu8JqCFwBE6NmKwBPqO0KqyFwBUyKmq0A6oParrASAlfApBZtOUwVAQD1UleVgUVbDgd2QEAjEbgCJpTrKNbm/QUuz1FFAIArnqoMbN5fwKorTIHAFTAZZ/mrkvIKl+dZbQXgiqdV15LyCspjwRQIXAGT8VT+iioCADzxVGWA8lgwAwJXwET2nyx1uyGrZXQkVQQAeOSsMtAyOtLleTZqIdQRuAIm8tymL91uyBrYM4GgFUCdku1xGtgzweU5ymMh1BG4AibxdZn0jy/OuDxnk5Se2i2g4wFgXump3SiPBVMicAVM4m9HIyh/BcAn6iqPxaorQpXPA9dbbrlF//u//+vrpwXC2v6Tpcr71vX6COWvADSGp/JYrLoiVPk8cH3wwQf1wAMPaPjw4frkk098/fRA2Ml1FGvKylzJzVsMq60AGqOuVddJS3cSvCLk+Dxw/elPf6o9e/ZozJgx+ulPf6rRo0crL4+PHIDGcNZsPVF8TnKRKMBqKwBveFp1dZw5S21XhBy/5LhGRETovvvuU15engYOHKhbb71V48ePV35+vj9uB1hWzZqttd9eWG0F4A1Pq64StV0RelwXcvNCWVmZdu/erU8++USffPKJPv74YxUWFmrbtm368Y9/rDFjxujJJ59UixYtfH1rwFJyHcVua7ZKNBsA4BszhiVqx5dFbhubOPNd+SUZocDnK66tWrXS+PHj9b//+79KTEzUs88+q5MnT+rAgQM6cOCAqqqqlJaW5uvbApazaMtht1UE7PHNaTYAwCecTQns8c1dnjd08fUICAU+X3EtKipSq1atXJ6LjY3Vc889p7g43myBuuz56ozL4zZJr95zHUErAJ9Jtsfp1Xuu0x3Pb3X5C7O71yMg0Pyy4uqJzWbT1q1bfX1bwFJyHcXKP/Wdy3PktQLwB0/5rvmnvmOTFkKC14HriBEjNHv2bL3zzjv13nzVr18/b28LWFauo1iTlu5UlYtljwgbVQQA+M+MYYmKcFFmoMqgPBZCg9epAgUFBZo/f77OnTsnm82m1q1b69prr63x1a0brSiB+nCWv3K9ScLQlfHNWW0F4DfJ9jh1at1CR4pqf+LjLI9Ffj2CyevANScnR5WVldq7d6927dqlXbt2acOGDdqwYYNstou/tlVWVno9UCAc1Cx/9UM2JV3RMqDjARB+eneMdxm4St+Xx1o6ISXAowIu8kmOa5MmTdS3b19NmDBBL730kg4cOKCVK1eqU6dO+tOf/uSLWwCWV1f5q6gIQxOv7xjAEQEIR+mp3RQd5T48oB0sgskvDQgk6Sc/+YmmT5+uw4cpoQHUh6fyV+3jmunhXpXq2S42oGMCEH4oj4VQ5rfAVZJGjhyplStX+vMWgGV4Kn/1wuhkdYgJ6HAAhDFneSx37WApj4Vg8Tpw/eUvf6nXX39du3fvrpXLevbsWR0/ftzbWwCWV1f5K1ZaAQQa5bEQirzenPXWW29p4cKFstlsatq0qXr37q1rr71WHTp00LvvvquePXv6YpyApT2dnUf5KwAhZ8awRG09WFjr9anKEJu0EBQ+KYeVn5+vnTt3aufOnfrXv/6lt99+W6dOnVLXrl21dOlSX4wTsCxPm7I6tW6hZHuciotZ2QAQeJ7KYzk3aVEaC4Hkk5avnTp1UqdOnZSWllZ97MKFC4qKivLF0wOW5Ww24G5TVu+O8YEcDgDU4q48lqGLTQloQY1A8tvmLIJWwDNnswHHmbMuz9t0sSwNAARTemo3t5u0nE0JyHdFoPi1qgAA9zw3G7i4KYtVDADB5mmTlvR9UwIgEAhcgSCoq9lAdFQEm7IAhIwZwxJpSoCQQOAKBIGnZgP2+Ob0AgcQUmhKgFBB4AoEgadmA2x0ABCKaEqAUEDgCgRYXc0GCFoBhCqaEiDYGlUOKz8/v0HXd+rUqTG3ASzHWf6KZgMAzMpTUwLKY8HfGhW4dunSRTabuw8LavthK1ggHDnLX7mrJOBsNgAAocxTUwJneSzy9OEvjQpcFy9e3KDAFUDd5a9oNgDALNw1JZC+L49FO1j4Q6MC1/vuu8/HwwCsrT7lr2g2AMAs0lO7aUPuCbe/jNMOFv7C5iwgACh/BcBKKI+FYGnUiusPHTx4UC+//LL27duns2drtq+02WzauHGjL24DmBblrwBYjbM81h3Pb3X5iznlseAPXgeue/fu1Y9//GN16NBBhw4dUp8+ffTNN9/o66+/1pVXXqlu3fj4E+GN8lcArMpZHmuLi1QoZ3ksXuPgS16nCjz++OO67bbb9Pnnn8swDL322ms6duyY3nvvPZWXl2vevHm+GKckqbS0VFOnTpXdbld0dLT69eunN998s96Pf/fdd5Wamqq4uDjFxMTo6quv1p///GefjQ9w5ensPMpfAbCsGcMSFeFiv3aVIWVk5wV+QLA0rwPXjz/+WPfee68iIi4+VVXVxUTtO+64Q7/61a/02GOPeXuLamlpacrMzNTs2bO1fv169e/fX+PGjdPy5cvrfGxGRobS0tLUq1cvvfXWW/rb3/6myZMn6/z58z4bH/BDnjZlUf4KgBU4y2O54tykBfiK16kCp0+fVuvWrRUREaGoqCidPn26+tx1112n3/zmN97eQpK0bt06bdiwQcuXL9e4ceMkSYMGDdLRo0c1bdo0jRkzRk2aNHH52F27dmnmzJl66qmnNH369OrjgwcP9snYAHc8bcqi/BUAq3BXHsu5SeuFcdcEflCwJK9XXDt06KBvvvlGktS9e3d99NFH1ec+++wzxcbGensLSdLq1asVGxur0aNH1zg+fvx4ORwO7dixw+1jX3zxRTVr1kxTpkzxyViA+sh1FGtzXoHLczaJ8lcALCM9tZvcVXffvL+AVVf4jNeB60033aR//OMfkqS7775bGRkZmjRpkiZPnqzHHntMd955p9eDlC5uAktKSlJkZM1F4j59+lSfd+ejjz5SUlKSVq1apZ49e6pJkybq2LGjHn30UVIF4BfOLlkl5ypcnmdTFgArcW7ScqWkvEJpi7YRvMInvE4VmDlzphwOhyRpxowZOnHihN544w3ZbDb913/9l5555hmvBylJRUVF6tq1a63jrVu3rj7vztdff63CwkI9+OCD+u1vf6vk5GRt3LhRGRkZOnbsmN544w2P9y4oKFBhYc08xUOHDkm6uGGsuNi//xjLyspq/InQ9+SaPW4LczeLjNCUmzs2aN4wB8AcQKjPgV/e3FH//KJI5ypqv/aVX6jSvDV79aexvYIwMusI9TnQGKWlpQ263mYYhrsUvJDSo0cPdevWTevXr69x/Pjx47Lb7Xrqqaf06KOPunxs06ZNdeHCBa1YsUJjx46tPv7www/rueee08GDB9W9e3e3954zZ47mzp3r8tzzzz+vTp06NeI7glV9XSbN/6yJVOuDM0PNm0hTrq5Uh5hgjAwA/OvrMun5z5uovFJy9Ro4vQ+vf6gpPz9fDz74oPbu3aurr766zut90oAgENq0aeNyVfXUqVOSvl95dffYEydO6LbbbqtxfPjw4Xruuef08ccfewxcJ0+eXCu39tChQxo5cqRSUlKUlJTUkG+lwcrKypSTk6OUlBTFxPAvPtRNfydPkqtKAjal9myrn49oeAks5gCYAzDLHPisIk/Zua5fAz+ruEI/H0QZwMYyyxxoiH379jXoep8GroWFhbU6Z0nyyYpk7969tWLFClVUVNTIc92zZ48kqVcv9x8/9OnTRydOnKh13LnY7Czl5U5CQoISEhJcnouNjVVcXGByFWNiYgJ2LzRe3knXH+HYJD04JNGrnyFzAMwBhPoceHBIov6eW+iyokreybKQHrtZhPocaIiGbuL3enNWSUmJJk2apJiYGLVv315XXXVVrS9fGDVqlEpLS7Vq1aoaxzMzM2W32zVgwAC3j/3JT34iSbXSDNatW6eIiAj179/fJ2ME6JIFINx52qjl7KYFNJbXK65Tp07V8uXLNXHiRPXp00fNmjXzxbhqGT58uIYOHar09HQVFxere/fuWrFihbKzs7Vs2bLqGq4TJ05UZmamDh8+rM6dO0u6WDLr5Zdf1uTJk/XNN98oOTlZH3zwgV566SVNnjy5+jrAG7mOYk1aupMuWQDC3oxhidp6sLDW62GVIU1aulOv3nMdv8ijUbwOXNeuXauMjAw99NBDvhiPR1lZWZo5c6ZmzZqlU6dOKTExsdaGq8rKSlVWVurSPWdRUVHasGGDHn/8cf3ud7/TqVOndNVVVykjI0OPPPKI38cN63OWv3JXSYAuWQDCibOblqumBI4zZ5W2aJuy0m/kdREN5nWqQHl5uXr37u2LsdQpNjZWCxYs0PHjx3Xu3Dnt3r27RtAqSUuWLJFhGOrSpUuN461bt9af/vQnnThxQufPn9f+/fv1q1/9qs78VqA+ns7Ocxu0SnTJAhB+PL3ulV+oUkZ2XuAGA8vwOmq7/fbbtXXrVl+MBTClXEexPjrgagftRdFREXTJAhB20lO7KTrKfZix9UAh+a5oMK8D1yeeeEJ//etf9eyzz+rgwYM6depUrS/AyhZtOexy96wk2eOb83EYgLCUbI9TVvqNssc3d3ne0MXXT6AhvA5ce/Xqpby8PE2bNk2JiYlq27ZtrS/AyvZ8dcblcZvEBgQAYS3ZHqdX77muVisCJ3evn4A7Xm/OmjVrlmw2d1MSsDbKXwGAZ87yWFtcpFQ5y2PxWon68jpwnTNnjg+GAZjT09l5lL8CgDpQHgu+4tMt9QcOHND27dt18OBBXz4tEJI8bcqi/BUAfM9ZHssVZ3ksNmqhPnwSuK5cuVKdO3dWUlKSbrrpJiUmJqpz5856++23ffH0QEjytCmL8lcAUFNd5bHYqIX68DpwXbduncaOHatWrVopIyNDS5cu1VNPPaVWrVpp7NixtdqsAlaQ6yjW5v0FLs/ZJMpfAcAP1FUei41aqA+vc1yffPJJ3XrrrVq7dm2NYv7Tpk3T8OHDNW/ePA0fPtzb2wAho64uWWzKAoDanOWxJi3dKceZs7XOs1EL9eH1iuunn36qyZMn1+pAZbPZNHnyZO3evdvbWwAhxVOXrOioCDZlAYAbzvJYES6KEVUZopsW6uR14NqkSROdP3/e5bkLFy7QUhWW4mlDVsvoSJoNAEAdPG3UopsW6uJ1VNm/f3/Nnz9fZ8/WXPY/d+6cnnnmGQ0YMMDbWwAhw9OGrIE9EwhaAaAe3G3UopsW6uJ1juvcuXM1ePBgde3aVaNHj1b79u11/PhxZWVlqaioSB9++KEvxgmEBE9dstiQBQD1k57aTWt2O1wuBLBJC554HbjedNNNev/99/Xoo4/qpZdekmEYioiI0IABA7RixQrdcMMNvhgnEHR0yQIA36CbFhrL68BVklJTU7V9+3Z99913On36tC677DK1aOE6fwUwo1xHsSYt3UmXLADwEbppoTEaFbjm5+friiuuUFRUlPLz82ud/+abb2r8vVOnTo0bHRAC6ip/RZcsAGg45yatI0W1P8lydtNiwyt+qFGB61VXXaXt27crJSVFXbp0kc3moq7FJSorKxs1OCAULNpy2G3QKtElCwAaq3fHeJeBq/R9N60Xxl0T4FEhlDUqcF28eLG6detW/d91Ba6AmXnaKBAdFcGmLABopPTUbtqQe8Lt4gAbtfBDjQpc77333ur/vu+++3w1FiDkeNqQZY9vTg4WAHiBblpoKLoDAB4s2nLY7YYsglYA8F5d3bSo64pLeR24fvjhh1q5cmX130+ePKnbb79d7du31z333KPy8nJvbwEEjbuPqdiQBQC+46mbFukCuJTXgeusWbOUm5tb/ffp06dr69atuuGGG/T222/r97//vbe3AILCU5oAG7IAwLfcva460wUAyQeB64EDB3TttddKkioqKrR69Wo9/fTTysrK0m9+8xutWLHC60ECwfB0dp7bNAE2ZAGAb6WndnObLjBp6U6CV0jyQeBaXFys+Ph4SdKuXbtUVlam//N//o8kKSUlxWWdVyDU5TqK9ZGLji4SaQIA4A+e0gWcdV0JXuF14JqQkKCDBw9Kkj744AN17txZHTt2lCSVlJQoKirK21sAAbdoy2GXPbQl0gQAwF88vb4667oivHnd8nXYsGF6/PHH9fnnn2vJkiU1SmXl5eWpS5cu3t4CCKhcR7E25RW4PGcTaQIA4C/UdUVdvF5x/d3vfqd+/frplVde0TXXXKMnnnii+tzy5ct1ww03eHsLIGCc7V1Lz1W4PH9zj7akCQCAnzjrutrjm7s8z0YteL3ievnllys7O9vluU2bNik6OtrbWwAB83R2ntvf9KOjIvTosMQAjwgAwouzruuIF7bW2iBbZUgZ2XlaOiElOIND0Pm1AUFcXJyaNm3qz1sAPuNpQ1bL6Ehlpd/IaisABICnjVpbDxSy6hrGvF5xvVRhYaHOnq3dsq1Tp06+vA3gF542ZA3smUDQCgAB1LtjvI4U1a6lbeji6/UL464J/KAQdF4HriUlJXr44Ye1YsUKt12yKisrvb0N4Hfukv7ZkAUAgZee2k1rdjtcLiiwSSt8eR24Tp06VcuXL9fEiRPVp08fNWvWzBfjAgLKU5csNmQBQOAl2+N0S4+22uIihcu5SYvX5vDjdeC6du1aZWRk6KGHHvLFeICgWLTlsNsuWWzIAoDgmDEsUVsPFrrcpEW6QHjyenNWeXm5evfu7YuxAEGR6yjWZjd1W+mSBQDB42mT1ub9BWzSCkNeB6633367tm7d6ouxAAHnrNta4qZuK12yACC43L0Ol5RX0AY2DHkduD7xxBP661//qmeffVYHDx7UqVOnan0Boaquuq1sygKA4EpP7aboKNfhSvmFKmVk5wV4RAgmrwPXXr16KS8vT9OmTVNiYqLatm1b6wsIRdRtBYDQ5+ym1TLa9bYc6rqGF683Z82aNUs2m80XYwECirqtAGAOyfY4DeyZoPd2O2qdo65rePE6cJ0zZ44PhgEEHnVbAcA8qOsKyYctX7/99lv9/e9/1xtvvKHTp0/76mkBv6BuKwCYi7OuqyvOuq6wPp8Err/97W9lt9s1fPhw3XPPPfryyy8lSYMHD1ZGRoYvbgH41NPZedRtBQCTmTEsUREushOrDGnS0p0Er2HA68B14cKFmjt3riZOnKi1a9fKML6PBkaMGKG1a9d6ewvApzxtyqJuKwCELk91XR1nzlIeKwx4Hbi++OKLeuSRR/T888/r1ltvrXHuP/7jP3Tw4EFvbwH4lKdNWdRtBYDQ5ul1uvxClRZtORy4wSDgvA5cv/jiC912220uz7Vs2VJnzpzx9haAz+Q6irV5v+suWWzKAoDQ56muq8RGLavzOnBt1aqVTp486fLckSNHlJCQ4O0tAJ+o7pJV7rpLFpuyACD0Oeu62uObuzzPRi1r8zpwHTx4sObPn6+ysrLqYzabTRUVFVq0aJHb1Vgg0BZtOeyxSxabsgDAHJLtcXr1nuvcbtQiXcC6vK7j+pvf/Eb9+/dXcnKyRo0aJZvNphdffFGffPKJ8vPz9dZbb/linIDX3H181DI6Un/97+tZbQUAE3Fu1DpSVLu0IekC1uX1imv37t21bds2JSUlaeHChTIMQ0uXLtXll1+urVu3qlOnTr4YJ+AVT3Vb6ZIFAObkbqMW6QLW5ZM6rsnJycrOzlZJSYm++uorFRcX6/3331dSUpIvnh7wmqe6rWzIAgBzSk/t5jZdICM7L/ADgt/5rHOWJDVr1kx2u13Nm7tOmAaCgbqtAGBNnuq6bj1QyKqrBXmd4ypJ77zzjt544w0dPXpU5eXlNc7ZbDbt3r3bF7cBGoW6rQBgXb07xrvMczV08fX/hXHXBH5Q8BuvV1x///vfKy0tTR999JGioqLUpk2bGl+tW7f2xTiBRqFuKwBYW3pqN7nIFpAkbd5fwKqrxXi94rpw4UJNmDBBL7/8spo0aeKLMQE+4azb6q4EFnVbAcD8ku1xuqVHW21xkRJWUl6htEXblJV+I6/3FuH1imtRUZF+9rOfEbQi5FC3FQDCw4xhiW67adEG1lq8DlxvvPFG7du3zxdjAXzKU91WfvsGAOtwdtNqGe36g2TqulqH14Hrc889p5deekl/+9vfdP78eV+MCfAadVsBILwk2+M0sKfrNvPUdbUOnzQgGDJkiEaNGqUWLVooLi6uxlerVq18MU6gQRZtOUzdVgAIM57qupIuYA1eb86aPn26XnzxRfXr109JSUlq2rSpL8YFeMXdx0LUbQUA66INrPV5HbguWbJEM2bM0FNPPeWL8QBe85QmQN1WALA2d3VdnekCLF6Ym9epApWVlRo6dKgvxgJ4LddRrElLd5ImAABhylO6wKSlO8l1NTmvA9dbb71V//znP30xFsArzrqtjjNnXZ4nTQAArM9TG1jHmbNKW7SN4NXEvE4V+PWvf60xY8YoJiZGd9xxh8tOWXTPQiB4qtsqkSYAAOHCXbqA9H1dV1rBmpPXgWvfvn0lSY888ogeeeQRl9dUVlZ6exugTp4S76OjIkgTAIAwkZ7aTRtyT7hdzGCjlnl5HbjOmjVLNpu7LsFAYHjakGWPb65X77mONAEACBPOhgSTlu50mT7GRi3z8jpwnTNnjg+GAXjHU91WglYACD/J9ji9es91GvHC1lrvD866rqQLmI/Xm7OAUEDdVgDAD3naqEW6gDk1asU1Pz9fV1xxhaKiopSfn1/n9Z06dWrMbYB6oW4rAMAd6rpaS6MC16uuukrbt29XSkqKunTpUmeOK5uz4E+0dwUAuJOe2k1rP3OQLmARjQpcFy9erG7dulX/d6A2Z5WWluqJJ57QW2+9pVOnTikxMVGPPvqoxo4d26DneeKJJ/Tkk0/q6quv1t69e/00WgQKaQIAAHdoA2stjQpc77333ur/vu+++3w1ljqlpaXpX//6lzIyMtSjRw8tX75c48aNU1VVlX72s5/V6zk+/fRTPfPMM2rXrp2fR4tAIE0AAFAX0gWswzSbs9atW6cNGzZo4cKFeuCBBzRo0CC98sorGjp0qKZNm1avdISKigqNHz9eDzzwgBITEwMwavgT7V0BAPVBG1jraFTgGhERoSZNmtT7yxdWr16t2NhYjR49usbx8ePHy+FwaMeOHXU+R0ZGhk6dOqUnn3zSJ2NC8NDeFQBQX7SBtY5GpQr8sOnA66+/rtLSUt15551q3769jh8/rjVr1igmJkYTJkzwyUD37t2rpKQkRUbWHHKfPn2qz99www1uH5+bm6t58+YpKytLsbGxDbp3QUGBCgsLaxw7dOiQpIt5t8XF/p3sZWVlNf6E9PwHeR7buya2i/H7zyWQmANgDoA54J3EdjEe28A+/0Ge5o8M7U9jrTgHSktLG3R9owLXS5sO/OEPf1D79u31wQcf1AgIS0pKNGTIELVo4fo3nIYqKipS165dax1v3bp19Xl3qqqqNGHCBKWlpen2229v8L0XLlyouXPnujyXk5OjEydONPg5GyMnJycg9zGDXV84V/Iv/ezHkGRTVIShPpHHtWnT8SCMzL+YA2AOgDnQOH0ipY0RTXShyibn+8X3DO36osA07xtWmgP1Kat6Ka87Zy1cuFC///3va61itmzZUtOnT9evfvUrTZs2zdvbSJLH6gWezj377LM6ePCg/va3vzXqvpMnT66VonDo0CGNHDlSKSkpSkpKatTz1ldZWZlycnKUkpKimJgYv97LDPafLFXRPz9xccam9nHN9MLoZPVs17BV9VDHHABzAMwB7/XvX6opK3N1ovjcD87YVHROsif3D+n3DyvOgX379jXoeq8D16+//rrWx/fVTx4Z6bPVyDZt2rhcVT116pSk71defyg/P1+zZs1SRkaGmjZtqjNnzki6uFGrqqpKZ86cUbNmzdS8eXO3905ISFBCQoLLc7GxsYqLC0wuZUxMTMDuFcqWrnVft3XxfSmWzm1lDoA5AOZA4/WPi9Pi+2LdtoFduvOkXhhnD87gGsBKc6Ch6ZteVxVISkrSs88+qwsXLtQ4fv78ef3hD3/w2e793r17a9++faqoqKhxfM+ePZKkXr16uXzcF198obNnz+qhhx7SZZddVv21bds27du3T5dddpkee+wxn4wRgUHdVgBAY9EG1ty8XnGdN2+eRo4cqa5duyotLU3t27fXiRMnlJWVpRMnTuidd97xwTClUaNG6ZVXXtGqVas0ZsyY6uOZmZmy2+0aMGCAy8f169dPmzZtqnV86tSp+vbbb/X666+rY8eOPhkj/I+6rQAAb1HX1by8DlzvuOMOZWdna+bMmXrppZdUVVUlm82mlJQUvf766xoyZIgvxqnhw4dr6NChSk9PV3Fxsbp3764VK1YoOztby5Ytqy67NXHiRGVmZurw4cPq3Lmz4uPjNXDgwFrPFx8fr4qKCpfnELpo7woA8BZtYM3L68BVkgYPHqzBgwfru+++0+nTp3XZZZf5rJrApbKysjRz5kzNmjWruuXrihUrarR8raysVGVlpQzDRXQDU8t1FGvz/gKX50gTAADUl6c2sJv3F7DqGsJ82jmrRYsW6tChg1+CVuliAu+CBQt0/PhxnTt3Trt3764RtErSkiVLZBiGunTp4vG5Nm/erL179/plnPA9Z8OBkvIKl+dJEwAANIS7942S8goaEoQwnwauhYWFys/Pr/UFeOvpbPcNB6KjIkgTAAA0SHpqN0VHuQ6Dyi9UKSM7L8AjQn14HbiWlJRo0qRJiomJUfv27XXVVVfV+gK8keso1kcHCl2eaxkdqaz0G/lIBwDQIMn2OGWl36iW0a6zJrceKGTVNQR5neM6depULV++XBMnTlSfPn3UrFkzX4wLqLZoy2G5y1ge2DOBoBUA0CjJ9jgN7Jmg93Y7ap0zxEatUOR14Lp27VplZGTooYce8sV4gFrc1dWziUoCAADvpKd205rdDpcLJNR1DT1epwqUl5erd+/evhgLUIunuq0392jLaisAwCvJ9jjd0qOty3POuq4IHV4Hrrfffru2bt3qi7EAtXiq2/roMN90ZQMAhLcZwxIVYat93FnXFaHD61SBJ554Qj/96U/VsmVL3XnnnWrTpk2ta1q3bu3tbRCmaO8KAPA3T3VdSRcILV4Hrr169ZIkTZs2TdOmTXN5TWVlpbe3QRiivSsAIFBoA2sOXgeus2bNks3mYn0d8BLtXQEAgUIbWHPwOnCdM2eOD4YB1ER7VwBAINEG1hx82jnrwIED2r59uw4ePOjLp0WYob0rACAYaAMb+nwSuK5cuVKdO3dWUlKSbrrpJiUmJqpz5856++23ffH0CDOLthymvSsAIODqagNLhYHg8zpwXbduncaOHatWrVopIyNDS5cu1VNPPaVWrVpp7NixWr9+vS/GiTDibgcn7V0BAP5UVxtYKgwEn9c5rk8++aRuvfVWrV27VhER38fB06ZN0/DhwzVv3jwNHz7c29sgTHiqJEB7VwCAv3lqA0uFgeDzesX1008/1eTJk2sErZJks9k0efJk7d6929tbIIw8nZ1HJQEAQFClp3Zz25AgIzsv8ANCNa8D1yZNmuj8+fMuz124cKFWQAu4k+so1kcHCl2eo5IAACBQnBUGXNl6oJBNWkHkdVTZv39/zZ8/X2fPnq1x/Ny5c3rmmWc0YMAAb2+BMLFoy2G5WGyVRCUBAEBguXvfMUQb2GDyOsd17ty5Gjx4sLp27arRo0erffv2On78uLKyslRUVKQPP/zQF+OExXmq22oTaQIAgMBKT+2mNbsdLhdUqOsaPF6vuN500016//331aVLF7300kt64okntGjRInXp0kXvv/++brjhBl+MExZWV93Wm3u05cUBABBQyfY43dKjrctz1HUNHq9XXCUpNTVV27dv13fffafTp0/rsssuU4sWrnNDgB+qq27ro8MSAzwiAACkGcMStePLIpfvUc66rrSCDSyf7pxq0aKFOnToQNCKBqFuKwAgFFHXNfQ0asW1T58+9b7WZrNREgtuUbcVABDKqOsaWhoVuLZu3Vo2m4sCZ5coLS3Vrl276rwO4W3RlsPUbQUAhLT01G5a+5mj1vtVlSHSBQKsUYHr5s2b3Z6rqKjQn//8Z/3mN7+RzWbTz372s8aODWHA3ccs1G0FAIQKZ13XI0W1PyEkXSCwfJrjunLlSiUnJ2vKlCnq27evdu3apb/85S++vAUsxFOaAHVbAQChxN37kjNdAIHhk8B18+bNGjBggMaMGaO4uDi9//77+vvf/65+/fr54ulhUaQJAADMwlMbWBoSBI5XgeuePXt0++23a/DgwSoqKtLy5cu1c+dODR482Ffjg4WRJgAAMAtPbWBJFwicRgWux44d07333qtrr71Wu3bt0nPPPad9+/Zp7Nixvh4fLIo0AQCA2ZAuEHyN2pzVo0cPnT9/XsOGDdP06dPVsmVL7dmzx+311157baMHCGsiTQAAYDZUFwi+RgWu586dkyStX79e2dnZbq8zDEM2m02VlZWNGx0sKddRrE15BS7PkSYAAAhVnqoLbN5fQE3XAGhU4Pr666/7ehwIE7mOYqUt2ua2xStpAgCAUNa7Y7zLwLWkvEJpi7bR8dHPGhW43nvvvb4eB8LEoi2H3Qat0VERpAkAAEJaemo3bcg94fK9rPxCFSkDfubTOq5AXdztvGwZHclvqQCAkJdsj1NW+o1qGe167Y8KA/5F4IqA8VRJYGDPBIJWAIApJNvjNLBngstzVBjwLwJXBAyVBAAAVkFDguAgcEXA0HAAAGAVNCQIDgJXBAQNBwAAVkNDgsAjcIXf5TqKNWnpTtIEAACW4ildYNLSnQSvfkDgCr9y1m11nDnr8jxpAgAAs/KULuA4c1Zpi7YRvPoYgSv8ylPdVok0AQCAuXl6H3PWdYXvELjCrzwlqNNwAABgdump3RQd5T6cYqOWbxG4wm88bciyxzen4QAAwPScDQns8c1dnmejlm8RuMJvPNVtffWe6whaAQCWkGyP06v3XEdd1wAgcIXfULcVABAuqOsaGASu8AvqtgIAwg11Xf2PwBV+QXtXAEC4oQ2s/xG4wi9IEwAAhBvSBfyPwBU+R5oAACBckS7gXwSu8DnSBAAA4Yp0Af8icIVP5TqKtXl/gctzpAkAAKzOU7rA5v0FrLp6icAVPpPrKFbaom0qKa9weZ40AQBAOHD3fldSXqG0RdsIXr1A4AqfWbTlsMovVLk8R3tXAEC48NQGtvxCFSkDXiBwhc+42zHZMjqS9q4AgLDhbAPbMjrS5XkqDDQegSt8wlMlgYE9EwhaAQBhJdkep4E9E1yeo8JA4xG4wieoJAAAQE1UGPA9Alf4BA0HAACoiYYEvkfgCq/RcAAAANdoSOBbBK7wGmkCAAC4RrqAbxG4wis0HAAAwD0aEvgWgSsajYYDAADUjYYEvkPgikaj4QAAAHWjIYHvELii0Wg4AABA3WhI4DsErmgUGg4AAFB/NCTwDQJXNAqVBAAAaBgqDHiPwBWNQsMBAAAahoYE3iNwRYPRcAAAgMahIYF3CFzRYKQJAADQOKQLeIfAFQ1GmgAAAI1DQwLvmCpwLS0t1dSpU2W32xUdHa1+/frpzTffrPNxWVlZGjdunLp3767mzZurS5cuuvvuu3Xw4MEAjNpaSBMAAMA7NCRoPFMFrmlpacrMzNTs2bO1fv169e/fX+PGjdPy5cs9Pu7pp5/Wd999p5kzZyo7O1vz5s3TJ598omuvvVaff/55gEZvDaQJAADgHRoSNJ7rSrghaN26ddqwYYOWL1+ucePGSZIGDRqko0ePatq0aRozZoyaNGni8rHvvfeeEhJq1k77z//8T3Xp0kV//OMf9eqrr/p9/FaQ6yjWh3knXZ4jTQAAgPpxNiQY8+ftLtumU2HAPdOsuK5evVqxsbEaPXp0jePjx4+Xw+HQjh073D72h0GrJNntdnXs2FHHjh3z+VitKNdRrLRF21R2rtLledIEAACoPxoSNI5pVlz37t2rpKQkRUbWHHKfPn2qz99www31fr4vvvhCR48e1ciRI+u8tqCgQIWFhTWOHTp0SNLFvNviYv9OrrKyshp/BsPzH+Sp/EKVy3PNIiN0z3Xt/P7/IZyFwhxAcDEHwBywnnuua6e1nzlqpeBVGRffd+ePTKxx3IpzoLS0tEHXmyZwLSoqUteuXWsdb926dfX5+qqoqNDEiRMVGxurhx9+uM7rFy5cqLlz57o8l5OToxMnTtT73t7IyckJyH1c2fWFMw3j0hoehpo3kaYkn5cj919y5AZjZOElmHMAoYE5AOaAtbRp1kSF5dIP3193fVGgTZuOu3yMleZAfn5+g643TeAqSTabi8Jn9Th3KcMwNHHiRG3dulWrVq3SlVdeWedjJk+eXCtF4dChQxo5cqRSUlKUlJRUr3s3VllZmXJycpSSkqKYmBi/3suV/SdLVfTPT1ycsSm1Z1v9fESii3PwpWDPAQQfcwDMAWta/22esnMLf3DUpqJzkj25v3q2i60+asU5sG/fvgZdb5rAtU2bNi5XVU+dOiXp+5VXTwzD0KRJk7Rs2TJlZmbqrrvuqte9ExISXObJSlJsbKzi4gKzKSkmJiZg97rU0rXuKwk8OCQxKGMKV8GaAwgdzAEwB6zlwSGJen9foct0gaU7T+qFcfZaj7HSHIiNja37okuYZnNW7969tW/fPlVU1Nx9t2fPHklSr169PD7eGbS+/vrrevXVV/Xzn//cb2O1GhoOAADgH54aElBdoDbTBK6jRo1SaWmpVq1aVeN4Zmam7Ha7BgwY4PaxhmHo/vvv1+uvv66XX35Z48eP9/dwLYOGAwAA+Je791OqC9RmmlSB4cOHa+jQoUpPT1dxcbG6d++uFStWKDs7W8uWLauu4Tpx4kRlZmbq8OHD6ty5syTpwQcf1GuvvaYJEyaod+/e+uc//1n9vM2aNdM111wTlO/JDGg4AACAf6WndnNbXWDRlsN6YRxxipNpAlfpYuvWmTNnatasWTp16pQSExO1YsUKjR07tvqayspKVVZWyjC+/+m/9957kqTFixdr8eLFNZ6zc+fOOnLkSEDGbza5jmJtyitweY40AQAAfMOZLnCkqPYnnJv3FyjXUcx77r+ZJlVAupjAu2DBAh0/flznzp3T7t27awStkrRkyRIZhqEuXbpUHzty5IgMw3D5RdDqmrPhQOm52h09JNIEAADwJXfvqyXlFUpbtI2UgX8zVeCKwHk6233DgeioCNIEAADwofTUboqOch2WlV+oUkZ2XoBHFJoIXFFLrqNYHx34YU25i1pGRyor/UY+sgAAwIeS7XHKSr9RLaNdZ3FuPVCo/Scb1mXKighcUcuiLYflYj+WJGlgzwSCVgAA/CDZHqeBPV3XjTckvbb9q8AOKAQRuKIWd3XjbKKSAAAA/pSe2k3ueoHuO14S0LGEIgJX1OCpbuvNPdqy2goAgB8l2+N0S4+2Ls8dO1Our8sCPKAQQ+CKGjzVbX10WGLgBwQAQJiZMSxRES6WXasM6YOvwzt0C+/vHrXQ3hUAgODy1AZ23xlbWG/SInBFNdq7AgAQGty9756ttOnnmbvDtq4rgSuq0d4VAIDQ4Kmu67mKKi3acjjAIwoNBK6QdHG1dfN+2rsCABAK6qrr6i61z+oIXFHd3rWknPauAACECk91XfNPfReW6QIErtCiLYdp7woAQAhKT+3mtsJAOKYLELjC7ccNtHcFACC4PFUYCMd0AQLXMOepkgDtXQEACD53KXvhmC5A4BrmqCQAAEBoI13gewSuYY6GAwAAhDbSBb5H4BrGaDgAAIA5kC5wEYFrGCNNAAAAcyBd4CIC1zBFwwEAAMwj2R6njvHRkmqvOG3eXxA2q64ErmGIhgMAAJhP8hUtJdVedi0pr1Daom1hEbwSuIYhGg4AAGA+E6/vqKgIFzl+ksovVIVFygCBaxii4QAAAObTs12sHu5VqdhmTVyeD4cKAwSuYYaGAwAAmFeHGOmmbq1dnguHCgMErmGGSgIAAJjbxOs7hm2FAQLXMEPDAQAAzK1nu9iwbUhA4BpGaDgAAIA1hGtDAgLXMEKaAAAA1hCuDQkIXMMIaQIAAFhDsj3ObbqAlRsSELiGCdIEAACwFnfv31ZuSEDgGiZIEwAAwFrSU7spOsp1KGfVhgQErmEg11GszfsLXJ4jTQAAAHNKtscpK/1GtYyOdHneihUGCFwtLtdRrLRF21RSXuHyPGkCAACYV7I9TgN7Jrg8Z8UKAwSuFrdoy2GVX6hyeS46KoI0AQAATC6cKgwQuFqcu48JWkZHKiv9RtIEAAAwOU8VBqyWLkDgamGeKgkM7JlA0AoAgEWES0MCAlcLo5IAAADhIVzSBQhcLYyGAwAAhIdwSRcgcLUoGg4AABBewiFdgMDVokgTAAAgvIRDugCBqwXRcAAAgPDjKV1g8/4CS6y6ErhaDA0HAAAIX+7e50vKK5S2aJvpg1cCV4uh4QAAAOErPbWboqNch3flF6pMnzJA4GoxNBwAACB8JdvjlJV+o1pGR7o8b/YKAwSuFkLDAQAAkGyP08CeCS7Pmb3CAIGrhVBJAAAASNatMEDgaiE0HAAAAJJ1GxIQuFoEDQcAAMClrNiQgMDVIkgTAAAAl7JiugCBq0WQJgAAAC5lxYYEBK4WQJoAAABwxWoNCQhcLYA0AQAA4IrVGhIQuJpcrqNYm/cXuDxHmgAAAOHNag0JCFxNLNdRrLRF21RSXuHyPGkCAADASg0JCFxNbNGWwyq/UOXyXHRUBGkCAABAknUqDBC4mpi75f2W0ZHKSr+RNAEAACDJOg0JCFxNylMlgYE9EwhaAQBADVZoSEDgalJUEgAAAA1hhXQBAleTouEAAABoCCukCxC4mhANBwAAQGOYPV2AwNWESBMAAACNYfZ0AQJXk6HhAAAAaCxP6QKb9xeE/KorgauJ7D9ZSsMBAADgFXfxQkl5hdIWbQvp4JXA1URe2/4VDQcAAIBX0lO7KTrKdQhYfqEqpFMGCFxNJPd4icvjNBwAAAD1lWyPU1b6jWoZHenyfChXGCBwNYmvy6SvzpS7PEfDAQAA0BDJ9jgN7Jng8lwoVxggcDWJD76OoJIAAADwGTNWGCBwNYljZS5mlqgkAAAAGseMDQlMFbiWlpZq6tSpstvtio6OVr9+/fTmm2/W67EFBQW67777dPnll6tFixa6/vrrtXHjRj+P2Df2nyzVN66zBKgkAAAAGs1sDQlMFbimpaUpMzNTs2fP1vr169W/f3+NGzdOy5cv9/i4c+fOafDgwdq4caMWLFigd999V+3atdOwYcO0ZcuWAI2+8V7b/pUM1V5xJU0AAAB4w2zpAq63k4WgdevWacOGDVq+fLnGjRsnSRo0aJCOHj2qadOmacyYMWrSpInLx7722mvau3ev/vGPf+j666+vfmzfvn01ffp07dixI2DfR2McLiyTZEg/CF5JEwAAAN5wpgscKardSv7gSdfVjILJNCuuq1evVmxsrEaPHl3j+Pjx4+VwODwGn6tXr1bPnj2rg1ZJioyM1M9//nPl5OTo66+/9tu4vZXrKFZxeYV+GLRKpAkAAADvuYsnissrQi5dwDSB6969e5WUlKTIyJqLxH369Kk+7+mxzutcPfbzzz/34Uh9J9dRrLRF23Sy5HytczQcAAAAvuCuIYHjzNmQ66RlmlSBoqIide3atdbx1q1bV5/39FjndQ19rHRxY1dhYWGNY4cOHZJ0ccNYcbF/fqDPf5DnslNW+7hmemF0sjrGym/3RmgpKyur8SfCD3MAzAH4aw50jJX+ck9fTVmZqxPF52qcK79Qpec/yNP8kYk+vadTaWlpg643TeAqSTab65JQdZ3z9rELFy7U3LlzXZ7LycnRiRMnPD6+sXYfcebsXjo+QxEV5XLk/kuOXL/cFiEsJycn2ENAkDEHwByAv+aArcJ13PHZkQJt2nTcL/fMz89v0PWmCVzbtGnjcmX01KlTkuRyRdUXj5WkyZMn18qtPXTokEaOHKmUlBQlJSXVOf7GWP9tno7nFv7gqE19urTVoEH++c0HoamsrEw5OTlKSUlRTExMsIeDIGAOgDkAf8+BYMQd+/bta9D1pglce/furRUrVqiioqJGnuuePXskSb169fL4WOd1l6rPYyUpISFBCQmu26LFxsYqLs4/O/sfHJKozQeL/p0ucLGqQHRUhB4ckui3eyK0xcTE8LMPc8wBMAfgrzlQM+64yN9xR2xsbIOuN83mrFGjRqm0tFSrVq2qcTwzM1N2u10DBgzw+Ni8vLwalQcqKiq0bNkyDRgwQHa73W/j9kayPU5Z6TdqWHJb2VtIw5LbKiv9RkpgAQAAn3PGHXf2tSuxfUvd2dcecnGHaVZchw8frqFDhyo9PV3FxcXq3r27VqxYoezsbC1btqy6huvEiROVmZmpw4cPq3PnzpKkCRMm6KWXXtLo0aOVkZGhhIQELVy4UPv379cHH3wQzG+rTsn2OM0fmahNm45r0CBWWgEAgP8k2+P0wrhrgj0Mt0wTuEpSVlaWZs6cqVmzZunUqVNKTEzUihUrNHbs2OprKisrVVlZKcMwqo81a9ZMGzdu1PTp0zVlyhR999136tevn9avX6/U1NRgfCsAAABoIFMFrrGxsVqwYIEWLFjg9polS5ZoyZIltY63a9dOmZmZfhwdAAAA/Mk0Oa4AAAAIbwSuAAAAMAUCVwAAAJgCgSsAAABMgcAVAAAApkDgCgAAAFMgcAUAAIApELgCAADAFAhcAQAAYAoErgAAADAFAlcAAACYAoErAAAATCEy2AMwq3PnzkmSDh065Pd7lZaWKj8/X/v27VNsbKzf74fQwxwAcwDMAVhxDjjjKGdcVRcC10Y6duyYJGnkyJHBHQgAAIDJHTt2TNdee22d19kMwzACMB7LOXPmjLZs2aIrr7xSzZo18+u9Dh06pJEjR+qdd95R9+7d/XovhCbmAJgDYA7AinPg3LlzOnbsmFJTUxUfH1/n9ay4NlJ8fLzuuuuugN6ze/fuuvrqqwN6T4QW5gCYA2AOwGpzoD4rrU5szgIAAIApELgCAADAFAhcAQAAYAoEribQtm1bzZ49W23btg32UBAkzAEwB8AcAHOAqgIAAAAwCVZcAQAAYAoErgAAADAFAlcAAACYAoErAAAATIHANYSVlpZq6tSpstvtio6OVr9+/fTmm28Ge1jwUklJiaZPn65bb71Vbdu2lc1m05w5c1xe+/HHH2vIkCGKjY1VfHy80tLS9MUXX7i89oUXXlBiYqKaNWumq666SnPnztWFCxf8+J2gMT788ENNmDBBiYmJiomJUYcOHXTXXXdp165dta7l529Nn376qe644w516tRJzZs3V+vWrXX99ddr2bJlta5lDoSPV199VTabTbGxsbXOMQ8uYSBkDR061IiPjzf+9Kc/GR9++KExadIkQ5LxxhtvBHto8MKXX35ptGrVyrjllluqf6azZ8+udd2+ffuMli1bGjfffLOxdu1aY9WqVcbVV19t2O12o6CgoMa18+bNM2w2m/HYY48ZmzZtMubPn280bdrUuP/++wP0XaG+fvrTnxqDBg0yFi5caGzevNlYuXKl8eMf/9iIjIw0Nm7cWH0dP3/r2rRpk/HAAw8Yf/nLX4wPP/zQeO+994yxY8cakozf/va31dcxB8LHV199ZbRq1cqw2+1GTExMjXPMg5oIXEPU2rVrDUnG8uXLaxwfOnSoYbfbjYqKiiCNDN6qqqoyqqqqDMMwjMLCQreB6+jRo43LL7/c+Pbbb6uPHTlyxIiKijKmT59efeybb74xoqOjjf/+7/+u8fgnn3zSsNlsxueff+6fbwSNcvLkyVrHSkpKjHbt2hmDBw+uPsbPP/wMGDDAuPLKK6v/zhwIHyNGjDDuvPNO4957760VuDIPaiJVIEStXr1asbGxGj16dI3j48ePl8Ph0I4dO4I0MnjLZrPJZrN5vKaiokJr1qzRT37yE8XFxVUf79y5swYNGqTVq1dXH8vOzlZ5ebnGjx9f4znGjx8vwzD0zjvv+HT88E5CQkKtY7GxsUpOTtaxY8ck8fMPV5dffrkiIyMlMQfCybJly7RlyxYtXLiw1jnmQW0EriFq7969SkpKqn4Rc+rTp0/1eVjX4cOHdfbs2eqf96X69OmjQ4cOqby8XNL3c6F37941rrviiit0+eWXM1dM4Ntvv9XHH3+sq6++WhI//3BRVVWliooKFRYWauHChfr73/+uGTNmSGIOhIuCggJNnTpVGRkZ6tixY63zzIPaCFxDVFFRkVq3bl3ruPNYUVFRoIeEAHL+fN3NAcMwdPr06eprmzVrppiYGJfXMldC3y9+8QuVlZVp5syZkvj5h4vJkycrKipKCQkJevjhh/X888/rgQcekMQcCBeTJ09Wz549lZ6e7vI886C2yLovQbB4+ji5ro+aYQ31nQPMFfP69a9/rTfeeEMvvPCCfvSjH9U4x8/f2h5//HFNmjRJBQUFeu+99/TLX/5SZWVl+tWvflV9DXPAulatWqX33ntPn3zySZ0/I+bB9whcQ1SbNm1c/nZ06tQpSa5/+4J1tGnTRpLrlfVTp07JZrMpPj6++try8nJ99913atGiRa1rfxgMIXTMnTtX8+bN05NPPqlf/vKX1cf5+YeHTp06qVOnTpKk22+/XZL02GOP6d5772UOWFxpaal+8YtfaMqUKbLb7Tpz5owk6fz585KkM2fOKCoqinngAqkCIap3797at2+fKioqahzfs2ePJKlXr17BGBYCpFu3bmrevHn1z/tSe/bsUffu3RUdHS3p+3ymH1574sQJffPNN8yVEDV37lzNmTNHc+bM0eOPP17jHD//8JSSkqKKigp98cUXzAGL++abb3Ty5En94Q9/0GWXXVb9tWLFCpWVlemyyy7T3XffzTxwgcA1RI0aNUqlpaVatWpVjeOZmZmy2+0aMGBAkEaGQIiMjNSdd96prKwslZSUVB/Pz8/Xpk2blJaWVn1s2LBhio6O1pIlS2o8x5IlS2Sz2TRy5MgAjRr19dvf/lZz5szRE088odmzZ9c6z88/PG3atEkRERHq2rUrc8Di2rdvr02bNtX6uu222xQdHa1NmzZp3rx5zANXgleJC3UZOnSocdlllxl//vOfjQ8//NC4//77DUnGsmXLgj00eGndunXGypUrjcWLFxuSjNGjRxsrV640Vq5caZSVlRmGcbHodGxsrHHLLbcY69atM7KysoxevXp5LDr9+OOPG5s3bzZ+//vfG82aNbNU0WmreOaZZwxJxrBhw4zt27fX+nLi529d999/v/E///M/xl//+ldj8+bNxttvv22MGTPGkGRMmzat+jrmQPhxVceVeVATgWsIKykpMR588EGjffv2RtOmTY0+ffoYK1asCPaw4AOdO3c2JLn8+vLLL6uv27lzpzF48GCjRYsWRlxcnDFy5Ejj0KFDLp9zwYIFRo8ePYymTZsanTp1MmbPnm2cP38+QN8R6is1NdXtz/6Hawn8/K1p8eLFxs0332xcfvnlRmRkpBEfH2+kpqYaf/nLX2pdyxwIL64CV8NgHlzKZhiGEfBlXgAAAKCByHEFAACAKRC4AgAAwBQIXAEAAGAKBK4AAAAwBQJXAAAAmAKBKwAAAEyBwBUAAACmQOAKAAAAUyBwBQAAgCkQuAIAAMAUCFwBAABgCgSuABAkx48fV2xsrMaOHVvj+Jo1axQVFaWZM2c26PnmzJkjm82mzz77TKNHj1arVq3UunVrPfLII6qoqND+/fs1bNgwtWzZUl26dNH8+fP98hwA4C8ErgAQJFdccYWmT5+ut956S7t27ZIkbd68WaNHj1Z6erqefPLJRj3vf/3Xf6lv375atWqV7r//fv3xj3/Uww8/rJEjR+qOO+7Q6tWr9Z//+Z+aMWOGsrKy/PYcAOBzBgAgaMrKygy73W4MHjzYyMnJMVq2bGmMHz/eqKqqavBzzZ4925Bk/OEPf6hxvF+/foYkIysrq/rYhQsXjLZt2xppaWk+fw4A8BdWXAEgiFq0aKF58+Zp48aNGjRokIYPH65XXnlFNput0c85YsSIGn9PSkqSzWbT8OHDq49FRkaqe/fuOnr0qN+eAwB8jcAVAIKsR48ekiSbzaYlS5aoSZMmXj1f69ata/y9adOmatGihaKjo2sdLy8v99tzAICvEbgCQBB9+umnGjFihG688UaVlpZq8eLFwR4SAIQsAlcACJL9+/frtttu0/XXX69Nmzbprrvu0pw5c/Ttt98Ge2gAEJIIXAEgCI4cOaIhQ4aoZ8+eWrVqlaKiopSRkaHTp0/rd7/7Xa3rbTabBg4cGPiBAkAIIXAFgAA7fvy4hgwZooSEBK1Zs0bNmzeXJCUmJmrChAlasGCBjhw5Un19aWmppIvlswAgnNkMwzCCPQgAgHvr1q3TiBEjtHv3bvXu3TvYwwGAoGHFFQBC3KZNmzR27FiCVgBhjxVXAAAAmAIrrgAAADAFAlcAAACYAoErAAAATIHAFQAAAKZA4AoAAABTIHAFAACAKRC4AgAAwBQIXAEAAGAKBK4AAAAwBQJXAAAAmAKBKwAAAEyBwBUAAACm8P8BQgb5sqya8r0AAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f06_filepath = os.path.join(analysis_directory_path, nonlinear_buckling_input_name + '.f06') # path to .f06 file\n",
"op2_filepath = os.path.join(analysis_directory_path, nonlinear_buckling_input_name + '.op2') # path to .op2 file\n",
"nonlinear_buckling_op2 = read_op2(op2_filename=op2_filepath, load_geometry=True, debug=None) # load results and geometry and do not print any debug or info message\n",
"buckling_load_vectors, critical_buckling_factors = pynastran_utils.read_nonlinear_buckling_load_from_f06(\n",
" f06_filepath=f06_filepath, op2=nonlinear_buckling_op2) # return list of nonlinear buckling loads and list of critical buckling factors\n",
"sol_106_buckling_load = np.abs(buckling_load_vectors[0, -1, 0]) # retrieve buckling load from first subcase, last node, x-component\n",
"print(f\"\"\"\n",
"SOL 106 buckling load: {sol_106_buckling_load:.0f} N\n",
"Critical buckling factor ALPHA = {critical_buckling_factors[0]:.1e}\"\"\") # print results\n",
"column_utils.plot_buckling_mode(op2=nonlinear_buckling_op2) # plot buckling mode"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"We observe that there is no noticeable difference between the linear and the nonlinear buckling load and shape. In fact, given the symmetry of the problem and the absence of material nonlinearities, there is no reason why the buckling load and shape should differ between the two analyses. Furthermore, the absolute value of $\\alpha$ is well below unity, which suggests that the points used to calculate the differential stiffness matrix are close to the critical point."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## SOL 106 - tangent stiffness matrix \n",
"\n",
"***\n",
"\n",
"The stability of an equilibrium point of a structure is assessed by looking at the Hessian of the potential energy $\\Pi$, that corresponds to the tangent stiffness matrix $\\boldsymbol{K}_T$ at that point:\n",
"\n",
"$$\\boldsymbol{H}_\\Pi=\\frac{\\partial^2\\Pi}{\\partial\\boldsymbol{q}^2}=\\boldsymbol{K}_T,$$\n",
"\n",
"where $\\boldsymbol{q}$ is the vector of the generalized coordinates. For stable points the tangent stiffness matrix is positive definite, while at critical points the matrix becomes singular. This means that at least one eigenvalue must be zero. As a consequence, we can assess the buckling load of Euler's column by looking at the load where the lowest eigenvalue $\\lambda$ of the tangent stiffness matrix becomes null.\n",
"\n",
"For this analysis we use the same input used for the nonlinear buckling method. We only need to add a command to calculate the smallest magnitude eigenvalue of the tangent stiffness matrix at every load increment and print it in the f06 file. For this purpose, we include an appropriate DMAP sequence in the executive control statements."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"nonlinear_buckling_bdf.executive_control_lines[1:1] = ['include \\'' + os.path.join(\n",
" os.pardir, os.pardir, 'resources', 'kllrh_eigenvalues_nobuckle.dmap') + '\\''] # include DMAP sequence"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"It should be noted that theoretically in this case we don't need the bulk data entry `PARAM,BUCKLE,2`, as we don't want to use SOL 106's nonlinear buckling method. However, without such entry Nastran does not update the tangent stiffness matrix after the last load increment. For this reason we keep the parameter `PARAM,BUCKLE,2` making sure that the DMAP sequence skips the buckling calculation.\n",
"\n",
"Let's apply a compression load equal to the buckling load found in the previous analyis rounded up to the next thousand, as we want to observe the lowest eigenvalue of the tangent stiffness matrix becoming negative."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"nonlinear_buckling_bdf.loads[force_set_id][0].mag = np.ceil(sol_106_buckling_load/1e3)*1e3 # change magnitude of force card to buckling load found by SOL 106 nonlinear buckling method rounded up to next thousand"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Let's define the input file name and run the analysis."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Nastran job euler_column_tangent_stiffness_matrix.bdf completed\n",
"Wall time: 11.0 s\n"
]
}
],
"source": [
"tangent_stiffness_matrix_input_name = base_input_name + '_tangent_stiffness_matrix' # name of input file\n",
"pynastran_utils.run_analysis(directory_path=analysis_directory_path, bdf=nonlinear_buckling_bdf,\n",
" filename=tangent_stiffness_matrix_input_name, run_flag=False)"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"Now we want to plot the lowest eigenvalue $\\lambda$ of the tangent stiffness matrix for every load increment versus the load history of the nonlinear analysis. We start by finding the load history of the nonlinear analysis using the function `read_load_displacement_history_from_op2` from the `pynastran_utils` module."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"op2_filepath = os.path.join(analysis_directory_path, tangent_stiffness_matrix_input_name + '.op2')\n",
"tangent_stiffness_matrix_op2 = read_op2(op2_filename=op2_filepath, debug=None) # read results from op2 file and do not print any debug or info message\n",
"_, load_history, _ = pynastran_utils.read_load_displacement_history_from_op2(op2=tangent_stiffness_matrix_op2)\n",
"applied_load_history = -[*load_history.values()][0][:,0] # retrieve applied load history"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Successively, we find the lowest eigenvalues in the f06 file using the function `read_kllrh_lowest_eigenvalues_from_f06` from the same module."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [],
"source": [
"f06_filepath = os.path.join(analysis_directory_path, tangent_stiffness_matrix_input_name + '.f06') # path to .f06 file\n",
"lowest_eigenvalues = pynastran_utils.read_kllrh_lowest_eigenvalues_from_f06(f06_filepath) # read the lowest eigenvalue of KLLRH matrices from f06 file"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"Finally, we plot the results nondimensionalizing the applied load by the critical buckling load found with SOL 106's nonlinear buckling method."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAApoAAAIHCAYAAAA/9PxiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAABJ0AAASdAHeZh94AAA85klEQVR4nO3deXxU9b3/8fckhIAzhrAktJFNCGjYZJEAFTSg4NZb0V6qRdEKooCI2gVarwtUq/Zy3ahEUW5FVr0UoS6gBQHBazBlUynhJwlblGJiMIaZS0Ig5/cHHUwgYObMOTPnzLyejwcPHs7MOfnAl0nefudzPsdjGIYhAAAAwGIJ0S4AAAAAsYmgCQAAAFsQNAEAAGALgiYAAABsQdAEAACALQiaAAAAsAVBEwAAALYgaAIAAMAWBE0AAADYolG0C4iU8vJyffDBB2rbtq2Sk5OjXQ4AAICrVFVVqbi4WJdddplSU1MbdEzcBM0PPvhAI0aMiHYZAAAArrZ8+XJdd911DXpt3ATNtm3bSjrxl5OZmWnZef1+v/Lz85WdnS2fz2fZeWEP1st9WDP3Yc3chzVzn2isWWFhoUaMGHEyUzVE3ATN4MflmZmZ6tatm2Xnraio0MGDB5WVlaWUlBTLzgt7sF7uw5q5D2vmPqyZ+0RzzUJpQeRiIAAAANiCoAkAAABbEDQBAABgC4ImAAAAbEHQBAAAgC0ImgAAALAFQRMAAAC2IGgCAADAFgRNAAAA2IKgCQAAAFsQNAEAAGALgiYAAABs0SjaBcSa6uM12rT3G3175KiaNW2sizs0V1IieR4AAMQfxwbNdevWaciQIfU+l5eXpwEDBkS4orOrPl6jme/v0ry8vfraf/Tk42m+ZI0e2F4TcjoROAEAQFxxbNAMevzxx08LnN27d49SNfU7XiPd95cd2lD0jTynPPe1v0pPr/pc24rLNXt0X8ImAACIG44Pmp07d3bc7uWpVh/waEPxN5Ik45Tngv+9ZmeJXlxXpHsu7xzR2gAAAKKF7bUwVR+v0YaDCaftZNbn5Q279eGuUlUfr7G9LgAAgGhzfNC8++671ahRI6WkpOjKK6/Uhx9+GO2S6tj2RYUOV3tO28msT0XlMd3y3/n60RNrNPP9XQROAAAQ0xz70XmzZs107733KicnRy1btlRhYaFmzJihnJwcvfPOO7ryyivPeGxJSYlKS0vrPFZYWChJ8vv9qqiosKzOknJ/yMcE+zY37flaz/w0i77NCAoEAnV+h/OxZu7DmrkPa+Y+0Vgzvz/0zOMxDKMhm3GOUF5erh49eqhFixb65JNPzvi6adOmafr06fU+N3PmTLVr186ymnZ969HzOxJNH39N2+O6so1rlgAAAMSp/fv3a/Lkydq+fbu6devWoGMcu6NZn9TUVP34xz/Wiy++qCNHjqhp06b1vm7ixIkaOXJknccKCws1YsQIZWdnKysry7Ka+lT49equLfI38OPzU31Y2lhD+nZWC2+SerVJYXfTZoFAQPn5+crOzpbX6412OWgA1sx9WDP3Yc3cJxprVlBQEPIxrgqakhTcgPV4znz5TXp6utLT0+t9zufzKSUlxdKaBv+gRiuKze1qVlQe12+W75TEzM1I8nq9lv87gL1YM/dhzdyHNXOfSK6Zz+cL+RhXpZlvvvlGb7/9tnr16qUmTZpEu5yTrsgwdGlmC0lq0NXnZxLs3bxr/mYuFAIAAK7n2B3NUaNGqV27drr44ovVqlUr7dq1S0899ZS++uorzZ07N9rl1ZGYID3z0ywt2lKqeXn7VOqvMnUeZm4CAIBY4tig2bNnT73++ut68cUX5ff71aJFCw0aNEjz589Xv379ol3eaZISE3TP5Z01PqeTNu4u08SFW3S48pjp8728Ybd6t0tV/44t+RgdAAC4kmOD5m9/+1v99re/jXYZIUtKTNDgzmkaN7ijnl71uenzBGdu0rcJAADciuRikwk5nTT0whMXJNG3CQAA4hFB0yZJiQmaPbqvfjWsi1r5kk2f59S+TQAAALdw7EfnsaB23+amvd/oUKBKDyzbrooj1aZmbs7ZsEfnp3nV0pusizs056N0AADgaATNCEhKTNDATi0lSUWlAdO9m99WVmvSoq2SmLkJAACcj4QSYfRuAgCAeEHQjDB6NwEAQLwgaEZBsHfzo98N1fyx2Tq3SXgdDC9v2K0Pd5WyswkAAByFoBlFtWduhiM4c/NHT6zRzPd3ETgBAIAjEDQdgL5NAAAQiwiaDkDfJgAAiEWMN3IIZm4CAIBYQ9B0GGZuAgCAWEHqcDB6NwEAgJsRNB2M3k0AAOBmBE2HY+YmAABwK4KmSzBzEwAAuA1B02Xo2wQAAG5B0HQZ+jYBAIBbMN7IhZi5CQAA3ICg6WLM3AQAAE5GkogR9G4CAACnIWjGCHo3AQCA0xA0YwgzNwEAgJMQNGMQMzcBAIATEDRjGH2bAAAgmgiaMYy+TQAAEE2MN4pxzNwEAADRQtCME8zcBAAAkUY6iEP0bgIAgEggaMYhejcBAEAkEDTjFDM3AQCA3QiacY6ZmwAAwC4ETUiibxMAAFiPoAlJ9G0CAADrETRxUu2+zcXjBmjWqN5q1jTJ1A6nR9K8vH3sagIAEMcImjhNcObmtT0zNHbQ+aYGuxuSSv1Vyl1bpLyiMgInAABxiIHtOKsJOZ20rbhca3aWyCOFHDqfWX1iMDzD3QEAiD/8xMdZWdW7yUVCAADEH4ImvpcVMze5SAgAgPhD0ESDWTVzk+HuAADEB4ImQhbuzE2GuwMAEB8ImggZfZsAAKAhCJow5dSZm/df0Tnkc9C3CQBAbCNoIizBmZsTh2Sqla+x6dtXztmwR29/eoCZmwAAxBDmaMISSYkJunVgBz296nNTx39bWa1Ji7ZKYuYmAACxgp/isEy4FwkF0bsJAEBsIGjCMlZdJETvJgAAsYGgCUtZMdy9NmZuAgDgXgRN2MKq4e7M3AQAwL0ImrAVfZsAAMQvgiZsRd8mAADxi/FGsF2wb3N8Tidt2vuNDgWq9MCy7ao4Un0yQIZizoY9Oj/Nq5beZF3coTkjkAAAcCiCJiImONxdkopKA8zcBAAgxvGTGVFB7yYAALGPoImooHcTAIDYR9BE1DBzEwCA2EbQRNQxcxMAgNhE0IRj0LcJAEBsIWjCMejbBAAgtjDeCI5i9czNlzfsVu92qerfsSXjjwAAiDCCJhzJqpmbwb7N4LzNm/ukWVkmAAA4C7Z44HhW9G4G+zbvX1og2jYBAIgMgiYcz4rezeDH7usLD2n1gXAuNQIAAA1F0IQrWDlzc+2BBL23o1R5RWVclQ4AgI3o0YSr1J65abZv88hxj36zfKck7pUOAICd+MkKV2LmJgAAzkfQhCsxcxMAAOdzVdCcM2eOPB6PfD5ftEuBA3CvdAAAnM01QfPLL7/Ur3/9a2VkZES7FDgM90oHAMCZXBM0x48fr0svvVTDhg2LdilwKPo2AQBwFlcEzQULFuiDDz5Qbm5utEuBg9G3CQCAszh+vFFJSYnuu+8+Pfnkk2rTpk20y4HDWX2v9Dkb9uj8NK9aepN1cYfmjEACACAEjg+aEydO1AUXXKAJEyY0+JiSkhKVlpbWeaywsFCS5Pf7VVFRYVl9gUCgzu9wjm5pSVJakm7pl6FZ6/eZOse3ldWatGirJKmlN0k39c3QmIFtCJwRxHvMfVgz92HN3Ccaa+b3+0M+xtFBc+nSpXrrrbe0detWeTwN77rLzc3V9OnT630uPz9fBw8etKrEOueFM3WqkbqmJmhHeYJOfDBuroOzLHBUs9bv09pP9+iOC2pE1ows3mPuw5q5D2vmPpFcs/3794d8jGODpt/v191336177rlHGRkZKi8vlyQdPXpUklReXq6kpCR5vd7Tjp04caJGjhxZ57HCwkKNGDFC2dnZysrKsqzOQCCg/Px8ZWdn11sLnOGynBq9kveFFm/6UmX/d8zkWU4E1B3lCdqdfL7uHNTOugJxRrzH3Ic1cx/WzH2isWYFBQUhH+PYoPn111/rq6++0lNPPaWnnnrqtOebN2+u6667TsuXLz/tufT0dKWnp9d7Xp/Pp5SUFKvLldfrteW8sM6vr0nV7QPbaM6b6zV/d2P5q46bPte8/C81oHNr9e/Yko/RI4T3mPuwZu7DmrlPJNfMzBxzxwbNH/zgB1q7du1pjz/55JP64IMPtHLlSrVq1SoKlcHNkhITdGGqodv6tzHdtyl9N3OTe6UDAHBmjg2aTZo0UU5OzmmPz507V4mJifU+BzTUmIFtVFByRGt2lsgjmboiXfpu5ua24nLNHt2XsAkAQC38VERcYuYmAAD2c13QnDt3rqnL64FT1b5X+uJxAzRrVG81a5pk+q5Cczbs0dufHlBeURl3FAIAQA7+6ByIlKTEBA3s1FKSVFQa0NOrPjd1ntozN+ndBADAhTuagJ24XzoAANYhaAK10LsJAIB1CJrAKWr3bs4fm61zm4TXYfLyht36cFcpO5sAgLhD0ATOICkxQYM7p2nc4I5hnSc4c/NHT6zRzPd3ETgBAHGDoAl8D/o2AQAwh6AJfA/6NgEAMIegCTSAlTM3PZLm5e1jVxMAEPMImkAIgjM3r+2ZobGDzjd160pDUqm/SrlrixjuDgCIaQxsB0yakNNJ24rLTd8v/ZnVJwbDM9wdABCr+KkGmGRV7yYXCQEAYhVBEwiDFTM3uUgIABCrCJqABayauclwdwBALCFoAhYKd+Ymw90BALGEoAlYiL5NAAC+Q9AELHbqzM37r+gc8jno2wQAxAKCJmCT4MzNiUMy1crX2PTtK+ds2KO3Pz3AzE0AgOswRxOwWVJigm4d2EFPr/rc1PHfVlZr0qKtkpi5CQBwF35SAREQ7kVCQfRuAgDchKAJRIBVFwnRuwkAcBOCJhAhVgx3r42ZmwAApyNoAhFm1XB3Zm4CAJyOoAlECX2bAIBYR9AEooS+TQBArGO8ERBFwb7N8TmdtGnvNzoUqNIDy7ar4kj1yQAZijkb9uj8NK9aepN1cYfmjEACAEQVQRNwgOBwd0kqKg0wcxMAEBP46QM4DL2bAIBYQdAEHIbeTQBArCBoAg7EzE0AQCwgaAIOxsxNAICbETQBF6BvEwDgRgRNwAXo2wQAuBHjjQCXYOYmAMBtCJqAyzBzEwDgFvxEAVyM3k0AgJMRNAEXo3cTAOBkBE3A5Zi5CQBwKoImECOYuQkAcBqCJhBj6NsEADgFQROIMfRtAgCcgvFGQAxi5iYAwAkImkAMY+YmACCa+CkBxAl6NwEAkUbQBOIEvZsAgEgjaAJxhJmbAIBIImgCcYiZmwCASCBoAnGMvk0AgJ0ImkAco28TAGAnxhsBcY6ZmwAAuxA0AUhi5iYAwHp85wdwGno3AQBWIGgCOA29mwAAKxA0AdSLmZsAgHARNAGcFTM3AQBmETQBNAh9mwCAUBE0ATQIfZsAgFARNAE0WO2+zcXjBmjWqN5q1jTJ1A6nR9K8vH3sagJADCNoAghZcObmtT0zNHbQ+aYGuxuSSv1Vyl1bpLyiMgInAMQgBrYDCMuEnE7aVlyuNTtL5JFCDp3PrD4xGJ7h7gAQe/huDiAsVvVucpEQAMQegiaAsFkxc5OLhAAg9hA0AVjGqpmbDHcHgNhAjyYAy4Xbtxkc7p7mS9bP+vxAncibAOBK7GgCsJyVfZuz1u/TnP+XwO4mALgQQROALU6duXn/FZ1DPkdwJ3RHeYJeyfvC2gIBALYjaAKwVXDm5sQhmWrla2zy9pWG5uV/obc/PcDMTQBwEXo0AUREUmKCbh3YQU+v+tzE0R5VVB7XpEVbJTFzEwDcwrHfobdt26Zrr71W7dq1U9OmTdWiRQsNHDhQCxYsiHZpAEyakNNJQy9MlySTO5snMHMTANzBsUGzvLxcbdu21eOPP64VK1Zo3rx56tChg0aPHq3HHnss2uUBMMGqi4SYuQkA7uDYj85zcnKUk5NT57Ef//jH2rNnj1566SU9+OCD0SkMQFiCFwmNz+mkjbvLNHHhFh2uPGb6fC9v2K3e7VLVv2NLPkYHAIdx3XflVq1aqVEjx+ZjAA1k1XD34MzNHz2xRjPf38VH6QDgII4PmjU1NTp27JhKS0uVm5ur9957T1OnTo12WQAsQt8mAMQux28NTpw4UbNnz5YkNW7cWDNnztRdd9111mNKSkpUWlpa57HCwkJJkt/vV0VFhWX1BQKBOr/D2VgvZ5pxXWe9kt5Ur20+oK8D1abOUbtv87n3dujOQe2sKxAh4X3mPqyZ+0Rjzfx+f8jHeAzDCPXucBG1f/9+lZSUqKSkRG+99ZZeeukl/fGPf9Svf/3rMx4zbdo0TZ8+vd7nZs6cqXbt+AEEONHxGmn3YY/81dLruxN05LgU+j6noaaJ0o0da+RLkjqea4jWTQAI3/79+zV58mRt375d3bp1a9Axjg+ap5owYYLmzJmjAwcOKC0trd7XnGlHc8SIEdq4caOysrIsqycQCCg/P1/Z2dnyer2WnRf2YL3cY/aH+zVr/b6wz9PSm6Sb+mZozMA2XCwUIbzP3Ic1c59orFlBQYEGDBgQUtB0/Efnp8rOztaLL76o3bt3nzFopqenKz09vd7nfD6fUlJSLK/L6/Xacl7Yg/Vyvvuu7KqCkiNas7NEJz4YN9fBeShQrVnr96mg5Ihmj+5L2Iwg3mfuw5q5TyTXzOfzhXyM677jrl27VgkJCerYMbwrVQE4W3Dm5qRL2+vcJPPnYeYmAESPY3c077zzTqWkpCg7O1utW7fW119/rSVLluj111/Xb37zmzPuZgKIHUmJCbpzUDt1rCqSt2Mv/XrZTmZuAoCLODZoDhw4UK+88opeffVVlZeXy+fz6aKLLtL8+fN1yy23RLs8ABGUmCANPL+5xg3uaPJe6ScEZ25yr3QAiAzHBs3bb79dt99+e7TLAOAgE3I6aVtxudbsLJFH330sHqrgzM1txeX0bQKAjSwJmsuXL9fChQu1b98+VVZW1nnO4/Hok08+seLLAIhzwb7NF9cVaV7ePpX6q0yd59S+zXsu72xdkQCAk8IOmjNmzNDUqVOVlpamzMxMxiIAsFXte6Vv2vuNDgWq9MCy7ao4Um1qh3POhj06P82rlt5kXdyhObubAGChsINmbm6uxowZo9mzZysxMdGKmgDgeyUlJmhgp5aSpKLSgOnezW8rqzVp0VZJoncTACwW9nfSsrIyjRo1ipAJIGq4XzoAOFPYQfOSSy5RQUGBFbUAgCnB3s1fDeuiVr5k0+dh5iYAWCvsoPnss89q1qxZevPNN3X06FEragKAkAV7Nz/63VDNH5utc5uE1xn08obd+nBXKTubABCGsINmZmamrrjiCl1//fU655xzlJKSUudXs2bNrKgTABokKTFBgzunadzg8O4eFpy5+aMn1mjm+7sInABgQtgXA02ZMkXPP/+8evXqpaysLDVu3NiKugAgLMzcBIDoCztozp07V1OnTtUTTzxhRT0AYAlmbgJA9IUdNI8fP65hw4ZZUQsAWIqZmwAQXWEHzeHDh2vjxo0aOnSoFfUAgOWYuQkA0RF20HzooYd04403yuv16tprr1WLFi1Oe019jwFANNC7CQCRE/Z3xosuukg7d+7UL3/5S11wwQVKS0s77RcAOAUzNwEgcsLe0Xz44Yfl8YRzLw4AiKzavZsbd5dp4sItOlx5zPT5Xt6wW73bpap/x5bsbAJALWEHzWnTpllQBgBEXu2Zm2b7NqXvZm7StwkAdfGdEEDc417pAGCPsHc0JWn58uVauHCh9u3bp8rKyjrPeTweffLJJ1Z8GQCwBTM3AcAeYQfNGTNmaOrUqUpLS1NmZqa8Xq8VdQFARDFzEwCsF3bQzM3N1ZgxYzR79mwlJiZaURMARA0zNwHAOmF/xysrK9OoUaMImQBiDr2bABCesIPmJZdcooKCAitqAQBHYeYmAIQn7KD57LPPatasWXrzzTd19OhRK2oCAMcI9m5+9Luhmj82W+c2Md9x5JE0L28fu5oA4kbYQTMzM1NXXHGFrr/+ep1zzjlKSUmp86tZs2ZW1AkAUVV75qZZhqRSf5XmfbSXsAkgLoR9MdCUKVP0/PPPq1evXsrKylLjxo2tqAsAHMmKe6U/+k6BXvxgNxcIAYh5YQfNuXPnaurUqXriiSesqAcAHM2qmZvBC4S2FZdr9ui+hE0AMSns72zHjx/XsGHDrKgFAFyhdt/m4nEDNGtUbzVrmhTSlelcIAQgHoQdNIcPH66NGzdaUQsAuEpw5ua1PTM0dtD5pj5Gl6SXN+zWh7tK6dsEEHPC/uj8oYce0o033iiv16trr71WLVq0OO019T0GALEknN7NispjuuW/8xnsDiDmhP2d7KKLLtLOnTv1y1/+UhdccIHS0tJO+wUAsc6KmZsMdgcQa8Le0Xz44Yfl8YRzzwwAiA2175c+76O9evSd0G5mcWrf5j2Xd7a+SACIoLCD5rRp0ywoAwBiR1Jigm79UQe98EGRyvxHTfVuztmwR+enedXSm6yLOzTno3QArmTqO9dzzz2nL774wupaACBmJCUm6NaBHUxfIPRtZbUmLdqqn7+8UT96Yo1mvr+Lj9MBuI6poPn444+rffv26t+/v2bMmKGiIkZzAMCpJuR00tAL0yUppNFHp6J3E4BbmQqa//znP7V69Wr169dPzz77rLp06aJevXrpscce044dO6yuEQBcyYoLhCRmbgJwL1NBMyEhQUOGDNHzzz+vL7/8UuvXr9eQIUM0Z84c9ejRQ1lZWXrwwQe1detWq+sFAFepPdx9/thsndskvNZ4Zm4CcBNLussvueQSPfPMM9q7d6/y8vL0k5/8RK+//rouvvhidezYUVOmTLHiywCAayUlJmhw5zSNG9wxrPMEZ27StwnADSy/jDE7O1t//OMftWvXLm3evFk333yz3nnnHau/DAC4En2bAOKJrfMyevXqpUcffVT/+Mc/7PwyAOAa9G0CiCemmoV69uzZ4Nd6PB598sknZr4MAMSk2oPdN+39RocCVXpg2XZVHKlm5iaAmGIqaLZo0eJ77wbk9/u1efNm7hoEAGeQlJiggZ1aSpKKSgN6etXnps4TnLkpifulA3AUU0Fz3bp1Z3zu2LFjeumll/T73/9eHo9Ho0aNMlsbAMSNCTmdtK24XGt2lsgjmR70Huzd3FZcrtmj+xI2AUSVpd+BlixZoq5du+qee+7RRRddpM2bN2v+/PlWfgkAiEn0bgKIRZYEzXXr1ql///668cYblZKSor/97W9677331KtXLytODwBxgZmbAGJNWEHzs88+0zXXXKPLL79cZWVlWrRokTZt2qTLL7/cqvoAIO4wcxNArDAVNIuLi3XbbbepT58+2rx5s5599lkVFBTopptusro+AIhbzNwE4HamPpfp0qWLjh49qquuukpTpkzRueeeq88+++yMr+/Tp4/pAgEgXgX7Nl9cV6R5eftU6q8ydZ5T+zbvubyzdUUCwFmYCppVVSe+2a1cuVLvvvvuGV9nGIY8Ho+OHz9urjoAiHPM3ATgZqaC5iuvvGJ1HQCAs2DmJgA3MhU0b7vtNqvrAAA0EDM3AbgF31UAwGWYuQnALQiaAOBCzNwE4AYETQBwMWZuAnAygiYAxABmbgJwIoImAMQA+jYBOFF4TT0AAMdg5iYApyFoAkCMYeYmAKew9bvF0KFDdcstt2jHjh12fhkAwBnQuwkgmmwNmuvWrdOiRYvUs2dPjR492s4vBQCoB72bAKLJ1qBZU1Ojw4cP680339QPf/hDO78UAOAMmLkJIFpsb7Txer265ppr9J//+Z92fykAwFkwcxNApNHRDQBxhr5NAJFie9AMBAJ2fwkAQAjo2wQQKZaONzp48KC2bdumbdu2aevWrdq2bZuKiop07NgxK78MACBMVs7c9Eial7dP4xl9BOAUpoNmYWGhNm/eXCdYlpaWSpIMw1Bqaqp69+6tf/u3f7OsWACAtayYuWlIKvVXKXdtkbLPb6EuLRItrhKAW5kKmpMnT9asWbMknQiVktSiRQs98MAD6t27t/r06aMOHTpYViQAwH4TcjppW3G51uwskUcKeWfzmdUnQmpLb5L6t/BoEH2bQNwz9RnH3LlzNWTIEK1cuVJ79+7V+PHjdejQIe3YsUOXXXYZIRMAXMiq3s1DgWqtKE7U/UsLuEgIiHOmguaRI0c0bdo0DR8+XO3atVNubq7effdd5efnq3v37lq5cqXVdQIAIsCKmZvBndD1hYe4SAiIc6aCpt/vV3Z2dp3Hhg8fru3bt2vo0KH68Y9/rIkTJ+rIkSOmC1uzZo3GjBmjCy+8UF6vV+edd56uu+46bd682fQ5AQANY9XMTYa7A/HNVNBMTk5W48aNT3u8WbNmWrhwoV5//XUtWbJEvXr10t///ndThb3wwgvau3ev7r33Xq1YsULPPfecSkpKNGDAAK1Zs8bUOQEAoQl35ibD3YH4Zul4o6B///d/1+DBg3XHHXfokksu0dGjR0M+x6xZs5Senl7nsauuukqZmZl6/PHHNXToUKvKBQCcQbBv88V1RZqXt0+l/ipT5wkOd99WXK7Zo/syBgmIE7a901u3bq233npLL7zwgqnjTw2ZkuTz+dS1a1cVFxeHWx4AoIFq920uHjdA91/ROeRzMNwdiE+2/y/l2LFjLTvXt99+qy1btqhbt26WnRMA0DDBmZsTh2Sqla+x6dtXztmwR29/ekB5RWV8lA7EOFs+OrfL3XffrUAgoP/4j/846+tKSkpODo8PKiwslHTiQqaKigrLagreYpNbbboD6+U+rJkz3djnh5q1fp+pY7+trNakRVslnZi5eVPfDI0Z2IaP06OI95n7RGPN/H5/yMe4Jmg+9NBDWrhwof70pz+pb9++Z31tbm6upk+fXu9z+fn5OnjwoOX15efnW35O2If1ch/WzFk61UhdUxO0ozxBJz4YN7e/WRY4qlnr92ntp3t0xwU1ImtGF+8z94nkmu3fvz/kYzxG8NY+DjZ9+nRNmzZNf/jDH/TAAw987+vPtKM5YsQIbdy4UVlZWZbVFggElJ+fr+zsbHm9XsvOC3uwXu7DmjlX9fEavZL3hV7bfEBfB6rDPt+kS9vrzkHtLKgMoeJ95j7RWLOCggINGDBA27dvb3Abo+N3NIMhc9q0aQ0KmdKJC4nqu5hIOnFBUUpKipUlSpK8Xq8t54U9WC/3Yc2c6dfXpOreK7tq4+4yTVy4RYcrj5k+17z8LzWgc2v179iSj9GjhPeZ+0RyzXw+X8jHOPqd/Oijj2ratGl68MEH9cgjj0S7HABAPawa7s7MTSD2OHZH86mnntLDDz+sq666Stdee602btxY5/kBAwZEqTIAQH0m5HTStuJyrdlZonD6Npm5CcQOxwbNt956S5L07rvv6t133z3teRe0lgJAXAkOd3/uvR16NW+vDpts2zx15uY9l4c+txOAMzg2aK5bty7aJQAAQpSUmKA7B7VTx6oiNcvsrUqjkR5Ytl0VR6plZntgzoY9Oj/Nq5beZF3coTm7m4DLODZoAgDcKzFB6tc+VSkpKSoqDejpVZ+bOk/tmZtpvmSNHtheE3I6ETgBl+CdCgCw1YScThp64YlJIGbvJiR917t51/zNXCgEuARBEwBgq2Dv5q+GdVErX7Lp83C/dMB9CJoAANslJSbonss766PfDdX8sdk6t0l4nVsvb9itD3eVsrMJOBxBEwAQMczcBOILQRMAEHH0bQLxgaAJAIg4+jaB+MB4IwBAVAT7NsfndNKmvd/oUKCKmZtAjCFoAgCiKikxQQM7tZQkZm4CMYZ3HwDAMejdBGILQRMA4Bj0bgKxhaAJAHAUZm4CsYOgCQBwJGZuAu5H0AQAOBp9m4B7ETQBAI5G3ybgXow3AgA4HjM3AXciaAIAXIOZm4C78I4CALgSvZuA8xE0AQCuRO8m4HwETQCAazFzE3A2giYAwPWYuQk4E0ETABAz6NsEnIWgCQCIGfRtAs5C0AQAxJTafZuLxw3QrFG91axpkqkdTo+keXn72NUETCJoAgBiUnDm5rU9MzR20PmmBrsbkkr9VcpdW6S8ojICJxAiBrYDAGLehJxO2lZcrjU7S+SRQg6dz6w+MRie4e5AaHiXAABinlW9m1wkBISGoAkAiAtWzNzkIiEgNARNAEBcsWrmJsPdge9H0AQAxKVwZ24y3B34fgRNAEBcom8TsB9BEwAQt06duXn/FZ1DPgd9m8CZETQBAHEvOHNz4pBMtfI1Nn37yjkb9ujtTw8wcxP4F+ZoAgDwL0mJCbp1YAc9vepzU8d/W1mtSYu2SmLmJiCxowkAQB3hXiQURO8mQNAEAKAOqy4SoncTIGgCAHAaK4a718bMTcQrgiYAAGdg1XB3Zm4iXhE0AQD4HvRtAuYQNAEA+B70bQLmMN4IAIAGCPZtjs/ppE17v9GhQJUeWLZdFUeqTwbIUMzZsEfnp3nV0pusizs0ZwQSYhJBEwCAEASHu0tSUWmAmZvAWfCvGQAAk+jdBM6OoAkAgEn0bgJnR9AEACAMzNwEzoygCQCABZi5CZyOoAkAgIXo2wS+Q9AEAMBC9G0C32G8EQAAFmPmJnACQRMAAJswcxPxjn+hAABEAL2biEcETQAAIoDeTcQjgiYAABHCzE3EG4ImAAARxsxNxAuCJgAAUULfJmIdQRMAgCihbxOxjvFGAABEETM3EcsImgAAOAAzNxGL+FcHAIDD0LuJWEHQBADAYejdRKwgaAIA4EBWz9ycs2GP3v70gPKKytjdRMTQowkAgIPVnrlptm9ToncT0cG/LgAAXMCqvk2J3k1EDkETAAAXsKpvU6J3E5FD0AQAwCVq920uHjdAs0b1VrOmSaZ3OD2S5uXtY1cTtqFHEwAAl7Fq5qYhqdRfpdc2H9APyZqwATuaAAC4mBW9mzNW79EjWxI1+8P97G7CUo4OmocPH9aUKVM0fPhwpaWlyePxaNq0adEuCwAAx7Cqd/NwtTRr/T4uEIKlHB00y8rK9NJLL6mqqkojRoyIdjkAADiSNb2bJ17NBUKwkqN7NNu3b69vvvlGHo9HX3/9tebMmRPtkgAAcCyrejfnbNij89O8aulN1sUdmjNrE6Y5Omh6POFOCgMAID5NyOmkbcXlWrOzRB59N9KoIRjuDqvwLwYAgBhkVe8mw90RDkfvaJpVUlKi0tLSOo8VFhZKkvx+vyoqKiz7WoFAoM7vcDbWy31YM/dhzZzltn6tNapPml7bfEAzVu8J+fjaw92fe2+H7hzUztoCYUo03md+vz/kY2IyaObm5mr69On1Ppefn6+DBw9a/jXz8/MtPyfsw3q5D2vmPqyZs/ywRjo3KVGHqyVzg5AMzfnfvaop3a3MFEN8iu4MkXyf7d+/P+RjYjJoTpw4USNHjqzzWGFhoUaMGKHs7GxlZWVZ9rUCgYDy8/OVnZ0tr9dr2XlhD9bLfVgz92HNnKsoeb9mrd9n8miPKo9LuQWJaulN0k19MzRmYBv6NqMkGu+zgoKCkI+JyaCZnp6u9PT0ep/z+XxKSUmx/Gt6vV5bzgt7sF7uw5q5D2vmPPdd2VUFJUdMXSBU26FAtWat36eCkiOaPbovYTOKIvk+8/l8IR/DvwwAAOKEVRcI1e7bZOYmzsbxO5orV65UIBDQ4cOHJUk7duzQX/7yF0nSNddco3POOSea5QEA4CrB4e7jczpp095vdChQpd+98ZkqKqtlpneTmZs4G8cHzQkTJmjfvu/6SZYsWaIlS5ZIkvbs2aMOHTpEqTIAANyr9nD3HV8cMt27ycxNnI3j/xXs3btXhmHU+4uQCQBA+MYMbKOuqSdmZIZzqxRmbuJUjg+aAADAXkmJCbrjghpNurQ9vZuwFEETAAAoMUG6c1A7ffS7oZo/NlvnNgmvu+7lDbv14a5SdjbjHEETAACclJSYoMGd0zRucMewzlNReUy3/He+fvTEGs18fxeBM04RNAEAwGkm5HTS0AtPzKSmbxNmETQBAMBpmLkJKzh+vBEAAIiO+mZuPrBsuyqOVJu6qxAzN+MPQRMAAJxV7ZmbRaUBPb3qc1PnYeZm/GFlAQBAg9G7iVAQNAEAQIPRu4lQEDQBAEBIgr2bzNzE9yFoAgAAU5i5ie9D0AQAAGGhbxNnQtAEAABhoW8TZ8J4IwAAEDZmbqI+BE0AAGAZZm6iNlYLAADYgt5NEDQBAIAt6N0EQRMAANiGmZvxjaAJAABsx8zN+ETQBAAAEUPfZnwhaAIAgIihbzO+MN4IAABEFDM34wdBEwAARAUzN2MfKwAAAKKO3s3YRNAEAABRR+9mbCJoAgAAR7By5qZH0ry8fexqRhlBEwAAOIoVMzcNSaX+Ks37aC9hM4oImgAAwJGs6Nt89J0ChrtHEUETAAA4klV9m1wgFD0ETQAA4Fi1+zYXjxugWaN6q1nTpJB2OLlAKHqYowkAABzPqpmbDHePLIImAABwlQk5nbStuFxrdpbII4V0NyGGu0cWf6sAAMBV6N10D4ImAABwndq9mw9dm2XqHPRu2o+gCQAAXCspMUG3/qiDWvkah3Xrypc37NaHu0rZ2bQYQRMAALhaUmKCbh3YIaRezVNVVB7TLf+dz8xNixE0AQCA61kx3F2ib9NqBE0AAOB6Vl0gRN+mtRhvBAAAYkLwAqHxOZ20ae83OhSo0gPLtqviSLWpj9WZuRk+giYAAIgpVg13Z+Zm+PibAgAAMYvezegiaAIAgJhF72Z0ETQBAEBMqz3cff7YbJ3bJLzOQWZuNhxBEwAAxIWkxAQN7pymcYM7hnUeZm42HEETAADEFfo2I4egCQAA4gp9m5HDeCMAABB3mLkZGQRNAAAQt5i5aa/4/tMDAAD8C72b1iNoAgAAiN5NOxA0AQAA/oWZm9YiaAIAAJyCmZvWIGgCAACcAX2b4SFoAgAAnAF9m+FhvBEAAMBZMHPTPIImAABAAzBzM3Sx9ycCAACwGb2bDUPQBAAACBG9mw1D0AQAADCBmZvfj6AJAAAQBmZunhlBEwAAwAL0bZ6OoAkAAGAB+jZPR9AEAACwSO2+zcXjBmjWqN5q1jTJ1A6nR9K8vH2u3tUkaAIAAFgsOHPz2p4ZGjvofFOD3Q1Jpf4q5a4tUl5RmSsDJwPbAQAAbDQhp5O2FZdrzc4SeaSQQ+czq08Mhq893N0t2NEEAACwkVW9m268SIigCQAAYDMrZm7WvkjolbwvrC3QJo4Omn6/X/fdd58yMjLUpEkT9erVS6+99lq0ywIAADDFqpmbr378hf5fucfxO5uODpo33HCDXn31VT3yyCNauXKl+vXrp5///OdatGhRtEsDAAAwLdyZm4erjiu3IFHDn8939HB3xwbNFStWaNWqVcrNzdVdd92lIUOG6OWXX9awYcP0m9/8RsePH492iQAAAKZY1bdZFqjW06s+189ezNP/HT1mYYXWcGzQXLZsmXw+n0aOHFnn8dtvv10HDhzQxx9/HKXKAAAAwnfqzM37r+hs+lxbi8t18WOrHbe76diguX37dmVlZalRo7rNsj179jz5PAAAgNsFZ25OHJKpVr7Gpm9f+X9HjzvuqnTHztEsKytTx46nN8q2aNHi5PNnUlJSotLS0jqPFRYWSjpxgVFFRYVldQYCgTq/w9lYL/dhzdyHNXMf1sw5buzzQ81avy+sc6zZWaLn3tuhOwe1s6iqE/x+f8jHODZoSpLHc+ZMf7bncnNzNX369Hqfy8/P18GDB8Ourb7zwj1YL/dhzdyHNXMf1iz6OtVIXVMTtKM8QScGGpnZ3zT0at5edawqUqKFn13v378/5GMcGzRbtmxZ767loUOHJH23s1mfiRMnntbbWVhYqBEjRig7O1tZWVmW1RkIBJSfn6/s7Gx5vV7Lzgt7sF7uw5q5D2vmPqyZs1yWU6NX8r7Qa5sP6OtAtYkzeHS4WmqW2Vv92qdaVldBQUHIxzg2aPbo0UOLFy/WsWPH6vRpfvbZZ5Kk7t27n/HY9PR0paen1/ucz+dTSkqKtcVK8nq9tpwX9mC93Ic1cx/WzH1YM+f49TWpuvfKrtq4u0wTF27R4crQryiv9iRZup4+ny/kYxx7MdD1118vv9+vpUuX1nn81VdfVUZGhvr37x+lygAAAOwX7nD3Zk0bW1xR6By7o3n11Vdr2LBhmjBhgioqKpSZmanFixfr3Xff1YIFC5SYmBjtEgEAAGw3IaeTthWXa83Okga93iOplS9ZF3dobm9hDeDYHU1JeuONNzR69Gg9/PDDuuqqq/Txxx9r8eLFuvnmm6NdGgAAQETUHu5+TuPv32gzJN06sL2SrLwSyKToV3AWPp9Pzz33nP75z3+qqqpKn3zyiW666aZolwUAABBRweHumx68Qr3bpv7rUaPOa4LXpw+9MF3jczpFsrwzcnTQBAAAwHfOadxI/zN+oCZd2l7nJtV9rpUvWb8a1kWzR/d1xG6m5OAeTQAAAJwuKTFBdw5qp45VRWqW2VvVniQ1a9pYF3do7piAGUTQBAAAcKHEBKlf+1RHj6RyVuwFAABAzCBoAgAAwBYETQAAANiCoAkAAABbEDQBAABgC4ImAAAAbEHQBAAAgC0ImgAAALAFQRMAAAC2IGgCAADAFgRNAAAA2CJu7nVeVVUlSSosLLT0vH6/X/v371dBQYF8Pp+l54b1WC/3Yc3chzVzH9bMfaKxZsEMFcxUDRE3QbO4uFiSNGLEiOgWAgAA4GLFxcXq06dPg17rMQzDsLkeRygvL9cHH3ygtm3bKjk52bLzFhYWasSIEVq+fLkyMzMtOy/swXq5D2vmPqyZ+7Bm7hONNauqqlJxcbEuu+wypaamNuiYuNnRTE1N1XXXXWfb+TMzM9WtWzfbzg9rsV7uw5q5D2vmPqyZ+0R6zRq6kxnExUAAAACwBUETAAAAtiBoAgAAwBYEzTClpaXpkUceUVpaWrRLQQOwXu7DmrkPa+Y+rJn7uGXN4uaqcwAAAEQWO5oAAACwBUETAAAAtiBoAgAAwBYETQAAANiCoFkPv9+v++67TxkZGWrSpIl69eql1157rUHHlpSU6Be/+IVatWqlc845RwMHDtT7779vc8Uwu2ZvvPGGfv7znyszM1NNmzZVhw4ddPPNN2vXrl0RqDq+hfM+q+3BBx+Ux+NR9+7dbagStYW7Zn/961912WWXKSUlRV6vV926ddNLL71kY8UIZ83Wrl2rYcOGKT09XT6fTz179tTMmTN1/Phxm6uOX4cPH9aUKVM0fPhwpaWlyePxaNq0aQ0+3pEZxMBphg0bZqSmphovvviisWbNGuOOO+4wJBkLFy4863GVlZVG9+7djTZt2hgLFiww/va3vxnXXXed0ahRI2PdunURqj4+mV2z7Oxs4yc/+Ynx5z//2Vi3bp0xf/58Iysry/D5fMb27dsjVH18MrtmtW3dutVITk42WrdubXTr1s3GamEY4a3ZE088YSQkJBgTJ040Vq5caaxevdp4/vnnjT/96U8RqDx+mV2zVatWGQkJCUZOTo6xfPlyY9WqVcY999xjSDImT54coerjz549e4xmzZoZl1566cm1euSRRxp0rFMzCEHzFO+8844hyVi0aFGdx4cNG2ZkZGQYx44dO+Oxs2bNMiQZH3300cnHqqurja5duxrZ2dm21Rzvwlmzr7766rTHvvzySyMpKckYO3as5bXihHDWLKi6utro1auXMXnyZOOyyy4jaNosnDXbtGmTkZCQYPzxj3+0u0zUEs6a3XzzzUZycrLh9/vrPD58+HAjJSXFlnphGDU1NUZNTY1hGIZRWloaUtB0agbho/NTLFu2TD6fTyNHjqzz+O23364DBw7o448/PuuxF1xwgQYOHHjysUaNGumWW25Rfn6+vvzyS9vqjmfhrFl6evppj2VkZKhNmzYqLi62vFacEM6aBT355JM6dOiQ/vCHP9hVJmoJZ82ef/55JScn65577rG7TNQSzpolJSWpcePGatq0aZ3HU1NT1aRJE1vqheTxeOTxeEwd69QMQtA8xfbt25WVlaVGjRrVebxnz54nnz/bscHX1XfsP/7xDwsrRVA4a1af3bt3a9++ferWrZtlNaKucNdsx44deuyxx/TCCy/I5/PZVie+E86arV+/XllZWVq6dKkuuOACJSYmqk2bNvrtb3+ro0eP2lp3PAtnzcaPH6+jR49q8uTJOnDggMrLyzV//nwtW7ZMU6ZMsbVumOPUDELQPEVZWZlatGhx2uPBx8rKymw5FuZZ+fd+7NgxjR07Vj6fT/fff79lNaKucNaspqZGY8aM0Q033KBrrrnGthpRVzhr9uWXX2rXrl2aPHmyJk+erNWrV+sXv/iF/uu//ku33367bTXHu3DWrH///lqzZo2WLVum8847T82bN9ftt9+uP/zhD/rVr35lW80wz6kZpNH3vyT+nG3b+vu2tMM5FuZZ8fduGIbGjh2rDRs2aOnSpWrbtq1V5aEeZtfs6aef1q5du/Tmm2/aURbOwuya1dTU6PDhw1q8eLFuuukmSdKQIUMUCAT07LPPavr06crMzLS8Xphfs82bN+v6669X//79NXv2bHm9Xq1Zs0YPPvigKisr9dBDD9lRLsLkxAxC0DxFy5Yt6039hw4dkqR6/2/BimNhnhV/74Zh6I477tCCBQv06quv6rrrrrO8TnzH7Jrt379fDz/8sJ588kk1btxY5eXlkk7sRNfU1Ki8vFzJycmn9ZUhfOF+bzx48KCuvPLKOo9fffXVevbZZ7VlyxaCpg3CWbO7775brVu31rJly5SYmCjpxP8cJCQkaNq0abr55pvVsWNHewqHKU7NIHx0fooePXqooKBAx44dq/P4Z599JklnndXXo0ePk68L9ViYF86aSd+FzFdeeUVz5szRLbfcYlutOMHsmu3evVtHjhzRvffeq+bNm5/89b//+78qKChQ8+bN9bvf/c72+uNROO+z+vrGpBPvPUlKSOBHkR3CWbNt27apb9++J0NmUL9+/VRTU6OCggLrC0ZYnJpBeHef4vrrr5ff79fSpUvrPP7qq68qIyND/fv3P+uxO3furHMl37Fjx7RgwQL1799fGRkZttUdz8JZM8MwNG7cOL3yyiuaPXs2/WIRYnbNevXqpbVr157266KLLlKHDh20du1aTZo0KRJ/hLgTzvvspz/9qSRp5cqVdR5fsWKFEhIS1K9fP+sLRlhrlpGRoU2bNp02nD0vL0+S1KZNG+sLRlgcm0GiNljJwYYNG2Y0b97ceOmll4w1a9YY48aNMyQZCxYsOPmaMWPGGImJicbevXtPPlZZWWl069bNaNu2rbFw4UJj1apVxvXXXx/1YanxwOyaTZo0yZBkjBkzxsjLy6vza8uWLdH4o8QNs2tWH+ZoRobZNTt69KjRp08fo1mzZsZzzz1nrFq1ypg6daqRmJhoTJo0KRp/lLhhds1mzpxpSDKuvvpqY/ny5cbf/vY3Y+rUqUajRo2MK664Ihp/lLixYsUKY8mSJcaf//xnQ5IxcuRIY8mSJcaSJUuMQCBgGIa7MghBsx6HDx82Jk+ebPzgBz8wGjdubPTs2dNYvHhxndfcdttthiRjz549dR4/ePCgceuttxotWrQwmjRpYgwYMMBYtWpVBKuPT2bXrH379oaken+1b98+sn+IOBPO++xUBM3ICGfNysrKjLvuusto3bq1kZSUZHTp0sWYMWOGcfz48Qj+CeJPOGu2dOlSY9CgQUarVq0Mr9drdOvWzXj00UdPG+IOa53t51JwjdyUQTyG8a8mGQAAAMBC9GgCAADAFgRNAAAA2IKgCQAAAFsQNAEAAGALgiYAAABsQdAEAACALQiaAAAAsAVBEwAAALYgaAIAAMAWBE0AAADYgqAJAAAAWxA0AcAGK1askMfjOfkrMTFR7du316RJk1RRUVHntcePH1d6erqeeeaZkI8FACdrFO0CACAWbdmyRZK0dOlSZWRkqLKyUkuWLNGsWbPk9/s1d+7ck69dv369SktLdcMNN4R8LAA4mccwDCPaRQBArLnhhhv03nvv6fDhw0pI+O7Do65du+qrr75SWVnZycfuvvtu5efn6+9//3vIxwKAk/HROQDYYPPmzerRo0edoChJKSkpCgQCJ//bMAwtW7ZMP/3pT0M+FgCcjqAJABYrKyvT/v37ddFFF9V5vLS0VNu3b1e/fv1OPvbRRx/pn//858mgGcqxAOB0BE0AsFiwx7J79+46duyYAoGAPv74Y40YMUJVVVX6/e9/f/K1f/nLX9SjRw917tw55GMBwOkImgBgsc2bN0uSJk+erKSkJPl8Pg0YMEBHjx7Ve++9pyFDhpx87RtvvHHax+YNPRYAnI6rzgHAYlu2bFGTJk20fv16eTweNW7cWOedd55atmxZ53X5+fnav39/naDZ0GMlac6cOXruuee0d+9edejQQW+//bbat29v+58PABqKoAkAFtuyZYt69uz5vf2US5cuVZcuXdS9e/eQj33sscf017/+Vf/zP/+jLl266KOPPqo3jAJANBE0AcBC3377rXbv3q3hw4d/72uXLl2qn/3sZyEfe/DgQc2YMUObNm062ds5ePDg8AoHABvQowkAFtqyZYsMw1Dfvn3P+rpt27apqKjotI/NG3LsqlWrNGDAgJMhEwCciqAJABYKXjXep0+fs75u6dKlat++fZ1Q2dBjDx06pNTU1PAKBYAI4M5AABAFXbt21dVXX62nnnoq5GPXrVunkSNHasOGDerSpYs+/fRTpaWl6bzzzrOhUgAwj6AJAC702GOPKTc3V4cPH9aFF16ot99+W61bt452WQBQB0ETAAAAtqBHEwAAALYgaAIAAMAWBE0AAADYgqAJAAAAWxA0AQAAYAuCJgAAAGxB0AQAAIAtCJoAAACwBUETAAAAtiBoAgAAwBYETQAAANiCoAkAAABbEDQBAABgi/8PdHFmftQTp10AAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, ax = plt.subplots()\n",
"ax.plot(applied_load_history/sol_106_buckling_load, lowest_eigenvalues[0, :], 'o')\n",
"plt.xlabel(\"$P/P_c$\")\n",
"plt.ylabel(\"$\\lambda$, N/mm\")\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"We can see that the lowest eigenvalue of the tangent stiffness matrix decreases linearly with the applied load and that after the last converged iteration is slightly negative. Let's print its value to the screen."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Lowest eigenvalue at the last converged iteration: -0.02 N/mm\n"
]
}
],
"source": [
"print(f\"Lowest eigenvalue at the last converged iteration: {lowest_eigenvalues[0, -1]:.2f} N/mm\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"The eigenvalue is indeed negative and this means that the structural equilibrium at the last converged iteration is unstable. Considering the observed linear trend, we can interpolate the curve to find the load where the lowest eigenvalue is zero, which corresponds to the critical buckling load."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"pycharm": {
"name": "#%%\n"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Load where tangent stiffness matrix becomes singular: 90578 N\n"
]
}
],
"source": [
"# Fit line to eigenvalue curve\n",
"fit = np.polyfit(lowest_eigenvalues[0, :], applied_load_history, 1)\n",
"line = np.poly1d(fit)\n",
"\n",
"# Find load where eigenvalue is zero and print its value\n",
"tangent_stiffness_matrix_buckling_load = line(0.)\n",
"print(f\"Load where tangent stiffness matrix becomes singular: {tangent_stiffness_matrix_buckling_load:.0f} N\")"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false
},
"source": [
"We observe that there is no evident difference w.r.t. the other numerical results."
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {
"collapsed": false,
"pycharm": {
"name": "#%% md\n"
}
},
"source": [
"## Conclusion \n",
"\n",
"***\n",
"\n",
"In this notebook we analyzed the buckling load of Euler's column with different methods. First we used Euler's analytical formula, and for the considered geometry and material we obtained a buckling load of 90962 N. Successively we used MSC Nastran to compute the buckling load of the column with different types of numerical analysis.\n",
"\n",
"We discretized the column with beam elements and we carried out a linear buckling analysis with SOL 105. This resulted in a buckling load of 90578 N, that is 0.42% lower than the analytical result. Then we switched to Nastran's nonlinear solution sequence, SOL 106. We first used the nonlinear buckling method available with the solution sequence, resulting in the same buckling load obtained with SOL 105. Then we used SOL 106 to find the load where the lowest eigenvalue of the tangent stiffness matrix becomes null, that is to say where the matrix becomes singular. The same value of buckling load was obtained also in this case. As a consequence, the different buckling calculations carried out with Nastran appear to be consistent. This is in line with our expectations because no nonlinear behavior is supposed to be observed up to the buckling point considering the symmetry of the problem.\n",
"\n",
"We will continue our investigation on the Euler's column studying its supercritical pitchfork bifurcation in our [next notebook](02_Supercritical_Pitchfork_Bifurcation_of_Euler_Column.ipynb)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.9.17"
}
},
"nbformat": 4,
"nbformat_minor": 4
}