diff --git a/doc/docs/Python_Tutorials/Adjoint_Solver.md b/doc/docs/Python_Tutorials/Adjoint_Solver.md index 1e100a2b4..117b8f8c5 100644 --- a/doc/docs/Python_Tutorials/Adjoint_Solver.md +++ b/doc/docs/Python_Tutorials/Adjoint_Solver.md @@ -8,7 +8,7 @@ Meep contains an adjoint-solver module for efficiently computing the gradient of This interface to this functionality is implemented entirely in Python using [autograd](https://github.com/HIPS/autograd) and [JAX](https://github.com/google/jax). The adjoint solver supports inverse design and [topology optimization](https://en.wikipedia.org/wiki/Topology_optimization) by providing the functionality to wrap an optimization library around the gradient computation. -There are six Jupyter notebooks that demonstrate the main features of the adjoint solver. +There are seven Jupyter notebooks that demonstrate the main features of the adjoint solver. - [Introduction](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/adjoint_optimization/01-Introduction.ipynb) @@ -22,4 +22,6 @@ There are six Jupyter notebooks that demonstrate the main features of the adjoin - [Near2Far Optimization with Epigraph Formulation](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/adjoint_optimization/06-Near2Far-Epigraph.ipynb) +- [Connectivity Constraint](https://nbviewer.jupyter.org/github/NanoComp/meep/blob/master/python/examples/adjoint_optimization/07-Connectivity-Constraint.ipynb) + More documentation will be available soon. diff --git a/python/examples/adjoint_optimization/07-Connectivity-Constraint.ipynb b/python/examples/adjoint_optimization/07-Connectivity-Constraint.ipynb new file mode 100644 index 000000000..db2b23446 --- /dev/null +++ b/python/examples/adjoint_optimization/07-Connectivity-Constraint.ipynb @@ -0,0 +1,391 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4b5849aa", + "metadata": {}, + "source": [ + "# Connectivity Constraint " + ] + }, + { + "cell_type": "markdown", + "id": "890e58ab", + "metadata": {}, + "source": [ + "For manufacturability, it is often necessary to impose connectivity constraint during topological optimization, so that all solid pixels are connected to one side (e.g. the supporting substrate). Meep adjoint supports this feature and it is illustrated in this simple tutorial. Please also note that the connectivity constraint is independent from the rest of Meep, and may be used separately.\n", + "\n", + "The basic underlying idea is to solve for an auxilary artificial temperature field. For connected structure, heat can freely flow out, so the temperature in the region should be low overall; on the other hand, for disconnected structure, heat will be trapped in islands, so temperature would be high in some area." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "a6964ee3", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using MPI version 3.1, 1 processes\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/opt/miniconda3/lib/python3.9/site-packages/numpy/core/getlimits.py:500: UserWarning: The value of the smallest subnormal for type is zero.\n", + " setattr(self, word, getattr(machar, word).flat[0])\n", + "/opt/miniconda3/lib/python3.9/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for type is zero.\n", + " return self._float_to_str(self.smallest_subnormal)\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import meep.adjoint as mpa" + ] + }, + { + "cell_type": "markdown", + "id": "4210a705", + "metadata": {}, + "source": [ + "Given (filtered and projected) design variable `v` with size `(nz, ny, nx)`, `mpa.constraint_connectivity(v, nx, ny, nz)` computes and simultaneously outputs the constraint value (negative for connected structure, positive for disconnected structure) and the gradient of the constraint with respect to `v`. Additional chain rule needs to be applied for gradient with repect to the unprojected and unfiltered design variable.\n", + "\n", + "The temperature field is also returned for reference, but it can be ignored when using the constraint.\n", + "#### IMPORTANT: The solver will reshape input `v` via `np.reshape(v, (nz, ny, nx))`. For 2D structure, please specify `ny=1`. The solver assumes the last index of the 0-axis of the array (bottom side of the image above) is where the structure should be connected to. So the user must first rotate the structure properly to match this condition; the gradient must also be post-processed and rotated so that the gradient corresponds correctly.\n", + "\n", + "There are two important hyperparameter, `p` and `cond_s`. The default values should be good in general but it is advised to check and tune them for each individual problem setup so that the constraint behaves as expected." + ] + }, + { + "cell_type": "markdown", + "id": "d86eeb17", + "metadata": {}, + "source": [ + "### disconnected structure" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7299cb1d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD8CAYAAADOpsDvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAARL0lEQVR4nO3df6zVd33H8eeLC5RSrf2BaxigJRG3EZ2pIbSm2axtzWg10GSNKXMduk72h7i6GhW3pS6df9i5Wd3CdMQ6mXEiq40lE0WHbYzLWqG26wSGvcMfwKiUlvqrVuCe1/74frGnt9x7vvfec+49l8/rkXzD+X7P53w+n3N73v38+p7zkW0iSjZjqisQMdUSBFG8BEEUL0EQxUsQRPESBFG8BEFMG5I+KemIpG+P8Lwk/Z2kQUmPSHp1k3wnFASSVkjaVxe6fiJ5RTTwKWDFKM9fAyypj7XAx5pkOu4gkDQAbKgLXgqslrR0vPlFdGL768CToyRZBfyzK/cD50ma3ynfmROo03Jg0PZ+AEmb60rsGekFs3WW53DOBIqMyfYTjh21/eLxvv53XneOn3hyqFHaBx/5xW7gmbZLG21vHENxC4ADbecH62uHR3vRRILgdAVeOjyRpLVUTRNzmMulumoCRcZk+3ff9f2JvP6JJ4f45vaXNEo7MP/RZ2wvm0h54zGRIGikjuSNAOfqgtyoVBgDLVqTVdwhYFHb+cL62qgmMjAeV4FRFmNOeKjR0QVbgT+oZ4kuA35ke9SuEEysJdgJLJG0mOrDfwPwexPIL85Q3WoJJH0WuAKYJ+kg8H5gFoDtjwPbgGuBQeBp4K1N8h13ENg+KWkdsB0YAD5pe/d484szkzFDXbpd3/bqDs8bePtY853QmMD2NqroixhRi/4eCvZ8YBxlMzCUIIjSpSWIohk40edf4U0QRE8ZpzsUhTMM9XcMJAiit6oV4/6WIIgeE0NoqisxqgRB9FQ1ME4QRMGqdYIEQRSulZYgSpaWIIpnxFCf/55DgiB6Lt2hKJoRxz0w1dUYVYIgeqpaLEt3KAqXgXEUzRZDTksQhWulJYiSVQPj/v6Y9XftYtrLwDgCGMo6QZQsK8YRQCuzQ1Gy6ga6BEEUzIgTuW0iSmaTxbIonbJYFmUzaQkiMjCOshnlSzVRtuonV/r7Y9bftYszQH58KwpnsmIc0fctQX+HaEx7tmh5RqOjCUkrJO2TNChp/Wmef4mkeyU9JOkRSdd2yrNjyZIW1ZnukbRb0s319QskfVXSo/W/5zd6F1GUamA80OjoRNIAsAG4BlgKrJa0dFiyvwC22L6EakfVf+iUb5PwOwm8y/ZS4DLg7XXB64EdtpcAO+rziGGq7xg3ORpYDgza3m/7OLAZWDUsjYFz68cvAv6vU6YdxwT1ZsiH68c/kbQXWFAXfkWdbBNwH/DeTvlFWaqBceMxwTxJu9rON9re2Ha+ADjQdn4QuHRYHn8JfEXSO4BzgKs7FTqmgbGki4FLgAeAi9p2C38MuGiE16wF1gLMYe5YioszxBhWjI/aXjbB4lYDn7L9t5JeA3xa0itsj7hXSOMgkPQC4PPAO23/WHo2um1b0mk35akjeSPAubqgzzfuiW7r8orxIWBR2/nC+lq7m4AVALb/U9IcYB5wZKRMG4WopFlUAfAZ23fXl38oaX79/PzRComytZjR6GhgJ7BE0mJJs6kGvluHpfkBcBWApN8A5gCPj5Zpk9khAXcCe21/uO2prcCa+vEa4J4GbyIKY8OJ1oxGR+e8fBJYB2wH9lLNAu2WdJuklXWydwFvk/RfwGeBt9ij7yHbpDt0OXAj8N+SHq6v/RnwQWCLpJuA7wNvapBXFKbqDnVvOcr2NmDbsGu3tj3eQ/WZbazJ7NA3YMQlv6vGUliUqd9XjHPbRPTUGKdIp0SCIHqsu92hXkgQRM/lO8ZRtGp2KD+5EgXL1ysjSHcoCpfZoQjy9coonC1OJgiidOkORdEyJoggQRCFyzrBdDZjAA0MwIwu/gdsGQ8NQWuoe3lOA1knmKaeXrWMgytPooHufSPUQ2LRFwY4+55vdi3PfmfDyQZfmJlKCYIRHH3lALtf//fMnTG7a3n+tPUMywdvYVFh38FLdyiKljFBBNWCWT9LEETPZWAcRbMzJojiiaHMDkXpMiaIouXeoQhX44J+liCInsvsUBTNGRhHpDsUkdmhKJudIIjIFGlExgRRNCNamR2K0vV5Q9B8b82IcakHxk2OJiStkLRP0qCk024gL+lNkvZI2i3pXzrlOZYtXAeAXcAh22+UtJhqR/ELgQeBG+tdxiOeq0tNQf0Z3AC8nmoj752Sttb7lJ1KswR4H3C57WOSfqVTvmNpCW6m2jHwlNuBO2y/DDhGtX9sxPN0sSVYDgza3l//D3czsGpYmrcBG2wfq8p2x62Fm+5jvBB4A/CJ+lzAlcBddZJNwHVN8oqyGGi11OgA5kna1XasHZbdAuBA2/nB+lq7lwMvl/Qfku6XtKJTHZt2hz4CvAd4YX1+IfBUva/sSJUBoH4jawHmMLdhcXHGMNB8neCo7WUTLHEmsAS4gmrH+69LeqXtp0Z6QZPNvN8IHLH94HhqZHuj7WW2l83irPFkEdOc3exo4BCwqO18YX2t3UFgq+0Ttr8LfIcqKEbUpDt0ObBS0veo+mBXAh8FzpN0qiU5XWUiKm54dLYTWCJpsaTZwA3A1mFpvkDVCiBpHlX3aP9omXYMAtvvs73Q9sV1oV+z/WbgXuD6OtkaoLCflIpmmg2KmwyM6+73OmA71STNFtu7Jd0maWWdbDvwhKQ9VJ/Rd9t+YrR8J7JY9l5gs6QPAA8Bd04grziTdXG1zPY2YNuwa7e2PTZwS300MqYgsH0fcF/9eD/VlFXEyAxu5Qa6KF6CIErX5zcPJQii9xIEUbSxLZZNiQRB9Fy+VBOR2aEondISRNGa3xIxZRIE0WPKwDgiLUFEa6orMLoEQfRW1gkiMjsU0fdjgvzuUBQvLUH0XLpDUTaT2yYi+n1MkCCInkt3KCJBEMVLEETJ5HSHIjI7FJGWICJBEEXLmCCCtAQR6vMv1eQu0iheWoLovXSHomgZGEeQliAiQRBFE5kditL52ZvoOh1NSFohaZ+kQUnrR0n3u5IsqeO+yE13tD9P0l2S/kfSXkmvkXSBpK9KerT+9/xmbyOK06UtXCUNABuAa4ClwGpJS0+T7oXAzcADTarXtCX4KPBl278OvIpq+8z1wA7bS4Ad9XnE83VvH+PlwKDt/baPU+2rveo06f4KuB14pkmmTXa0fxHw29RbtNo+bvupuvBNdbJNwHVNCozyjKE7NE/SrrZj7bCsFgAH2s4P1teeLUt6NbDI9heb1q/JwHgx8DjwT5JeBTxI1dRcZPtwneYx4KLTvbh+I2sB5jC3ab3iTNJ8duio7Y59+JFImgF8GHjLWF7XpDs0E3g18DHblwA/Y1jXp95A+bRv1fZG28tsL5vFWWOpW5wJXM0ONTkaOAQsajtfWF875YXAK4D7JH0PuAzY2mlw3CQIDgIHbZ8aZNxFFRQ/lDQfoP73SIO8okTdGxPsBJZIWixpNnADsPWXxdg/sj3P9sW2LwbuB1ba3jVaph2DwPZjwAFJv1ZfugrYUxe+pr62Brin0duI4nRritT2SWAdsJ1qcmaL7d2SbpO0crz1a7pY9g7gM3X07QfeShVAWyTdBHwfeNN4KxFnuC6uGNveBmwbdu3WEdJe0STPRkFg+2HgdP2qq5q8PgqWPcuidCJ3kUYkCCLSHZqmzt/XYtn9f8jAQPdugTx5coALvjPUtfymjQTB9HTu3d/ivC/N6Xq+rZ8/0++fie7KN8umL584ztCJ41NdjTNDgiBK1+9fqkkQRM+lOxRly2JZBAmCKFtWjCMAtfo7ChIE0VsZE0SkOzRttX7rEg699mzcxV9mUgsW3Pc0M77xcPcynQ4SBNPTodeezdf/+EPM1ayu5flTn+B1Q+9m4Te6luW0kJZgmvIMmKtZzJ0xu2t5tlqtrrYs00aCIIrm3DYRhcs6QQSA+zsKEgTRc2kJomxZLIvIwDgiQRCFMxkYR2RgHJEgiJJlsSzCzpdqItIdiuKlOxRlM5DuUBSvv2Og8WbeEePWrT3LACStkLRP0qCk520gL+kWSXskPSJph6SXdsqzURBI+lNJuyV9W9JnJc2pdxB8oK7M5+r9zCKeRy03OjrmIw0AG4BrgKXAaklLhyV7CFhm+zepdlr96075NtnRfgHwJ3XGrwAGqLbOvB24w/bLgGPATR3fRZSn6fatzVqC5cCg7f22jwObgVXPKc6+1/bT9en9VHsdj6ppd2gmcLakmcBc4DBwJVWkAWwCrmuYVxSkWixzowOYJ2lX27F2WHYLgANt5wfrayO5CfhSpzp2HBjbPiTpb4AfAD8HvgI8CDxV7ys7amXqN7IWYA5zOxUXZ6Lmd5EetT3q7vNNSfp9qh1XX9spbZPu0PlUTc5i4FeBc4AVTStje6PtZbaXzeKspi+LM8gYWoJODgGL2s4X1teeW550NfDnVLvZ/6JTpk26Q1cD37X9uO0TwN3A5cB5dfdoxMpEdHlMsBNYUk/KzKYam25tTyDpEuAfqQLgSJNMmwTBD4DLJM2VJKoNvPcA9wLX12nWAPc0ehtRmGYzQ01mh+ru9zpgO7AX2GJ7t6TbJK2sk30IeAHwr5IelrR1hOx+qcmY4AFJdwHfAk5STUFtBL4IbJb0gfranR3fRZSpi1+qsb0N2Dbs2q1tj68ea56NVoxtvx94/7DL+6mmrCJGlh/fiiBfr4zo93uHEgTRc2r1d38oQRC9ZcayWDYlEgTRU6LxQtiUSRBE7yUIongJgihaxgQRmR2K4jndoShcfpA3gowJIrJOEJEgiKLZMNTf/aEEQfReWoIoXoIgipYf5I0wOGOCKJnJwDgiY4KIBEGULTfQRekM5FbqKF5agihbbpuI0hmcdYIoXlaMo3gZE0TR7MwORaQliMIZDw1NdSVGlSCI3sqt1BH0/a3UTTfzjhgXA2650dGEpBWS9kkalLT+NM+fJelz9fMPSLq4U54Jgugt11+qaXJ0IGkA2ABcAywFVktaOizZTcAx2y8D7gBu75RvgiB6zkNDjY4GlgODtvfbPg5sptpovt0qYFP9+C7gqnrr4RHJkzh9Jelx4GfA0UkrdGLmMX3qCr2p70ttv3i8L5b0Zap6NTEHeKbtfKPtjW15XQ+ssP1H9fmNwKW217Wl+Xad5mB9/r91mhH/LpM6MLb9Ykm7bC+bzHLHazrVFfqzvrZXTHUdOkl3KKaTQ8CitvOF9bXTppE0E3gR8MRomSYIYjrZCSyRtFjSbOAGYOuwNFuBNfXj64GvuUOffyrWCTZ2TtI3plNdYfrVd0xsn5S0DtgODACftL1b0m3ALttbgTuBT0saBJ6kCpRRTerAOKIfpTsUxUsQRPEmLQg6LXdPNUmLJN0raY+k3ZJurq9fIOmrkh6t/z1/qut6iqQBSQ9J+rf6fHF9q8BgfevA7Kmu43QwKUHQcLl7qp0E3mV7KXAZ8Pa6juuBHbaXADvq835xM7C37fx24I76loFjVLcQRAeT1RI0We6eUrYP2/5W/fgnVB+uBTx3GX4TcN2UVHAYSQuBNwCfqM8FXEl1qwD0UV373WQFwQLgQNv5wfpaX6rvPLwEeAC4yPbh+qnHgIumql7DfAR4D89ui3ch8JTtk/V5X/+N+0kGxsNIegHweeCdtn/c/ly96DLlc8qS3ggcsf3gVNflTDBZi2VNlrunnKRZVAHwGdt315d/KGm+7cOS5gNHpq6Gv3Q5sFLStVQ3nZ0LfBQ4T9LMujXoy79xP5qslqDJcveUqvvUdwJ7bX+47an2Zfg1wD2TXbfhbL/P9kLbF1P9Lb9m+83AvVS3CkCf1HU6mJQgqP/PdGq5ey+wxfbuySh7DC4HbgSulPRwfVwLfBB4vaRHgavr8371XuCW+paBC6mCOjrIbRNRvAyMo3gJgihegiCKlyCI4iUIongJgihegiCK9/8a9sarxBDseQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "nz, ny, nx = 100, 1, 50\n", + "v = np.zeros((nz,ny,nx))\n", + "v[30:65,:,20:30]=1\n", + "v[70:,:,20:30]=1\n", + "plt.imshow(v[:,0,:])\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "markdown", + "id": "2ad6e653", + "metadata": {}, + "source": [ + "The structure above has one part that is connected to the bottom, but one part is floating above and is disconnected. We can see it has a positive constraint value of around 104" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "8e8f0e8d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "104.04171048992157" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "T, f, df_dv = mpa.constraint_connectivity(v, nx, ny, nz, p=6, cond_s=1e6)\n", + "f" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "df183225", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAL0AAAD7CAYAAAA/1pVBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7c0lEQVR4nO2da6w1WVnnf8+qvd/37W4uTQMahsahE4iGmIBOx8EwMcplBtGIH4gBDXEcJj0fFEFMRP3CfJgPmBiRD2ZmOqLDJIzAcMkYR0GDGOOXHq5RARlaFOkerg0NfXnfc/au9cyHdalVa6/au/Y5VbXr9Nn/k5196rLrsupf//qvZz1rlagqRxxxmWAOfQBHHDE1jqQ/4tLhSPojLh2OpD/i0uFI+iMuHY6kP+LS4VykF5GXishnReReEfnVoQ7qiCPGhJw1Ti8iFfB/gZcA9wEfAV6lqp8e7vCOOGJ4LM7x2x8A7lXVzwOIyDuBlwOdpH/ybUa/6xnn2eURU+OTf736uqo+9ay//zc/cos+8I2617of++uTD6rqS8+6r744DwOfDnwxmb4P+Jf5SiJyF3AXwO1Pr/jzP/mOc+xyOypktG3PFTXjtqjf9vT7v3Ce3z/wjZr/88Hv6rVu9bTPPeU8++qL0Suyqnq3qt6pqnc+5cnHevNlgwK2599UOI/S3w88I5m+3c874ogIRVlpP3szFc5D+o8AzxaRO3BkfyXw04Mc1RGPKUyp4n1wZtKr6lpEfgH4IFABv6eqnxrsyI54TEBR6pll8p4rlKKqfwz88UDHcsRjFHbkyva+OMYPjxgVyvgRpn1xJP0Ro+Oo9EdcKiiweix5+iOO2AVFj/bmiEsGhXpenD+S/ohx4Vpk54Uj6Y8YGUI9s5yoI+mPGBWuInsk/Wio0UuVaTm3CmIJLk4/r2vymCL9EfOEPSr9EZcJR6U/4tJBEeqZjT9wJP0Ro2Nu9mZet+AAuAiVuyFwUc5TEU616vXpAxH5JRH5lIj8rYj8gYhcE5E7ROQePyrHu0TkyrZtPOZIf8S84BqnTK/PLojI04FfBO5U1e/F9eN4JfAbwFtU9VnAN4HXbNvOkfRHjI7aN1Dt+vTEArhJRBbAzcCXgBcC7/HL3w785K4NHHHEaFAVau2trU8RkY8m03er6t3NtvR+EflN4J+A68CfAh8DHlTVtV/tPtxIHZ04kv6I0WH7q/jXVfXOroUi8iTc2Ep3AA8C/xPYe5ycI+mPGBWuIjsYzV4M/IOqfg1ARN4HvAC4VUQWXu13jspx9PRHjIohK7I4W/N8EblZRAR4EW5EvQ8Dr/Dr/Czwv7Zt5Ej6I0ZHrdLrswuqeg+uwvpx4G9w/L0beCPwBhG5F3gy8LZt2znamyNGxdAtsqr6JuBN2ezP48ZW7YUj6Y8YHbZ/9GYSHEl/xKhwCWdH0h9xiaAIq54pBlPhSPojRoUq+zROTYIj6Y8YGbJP49QkOJL+DDA9POrcRuo9FJSj0k+CofvK9iH5rt8MeRNclLTigGNF9ohLBUVm14lkUtIrYAvjGhqZV6EEnEXhd21rrrandF2GgBsCZF7aOoujCQU+J/IPSfh8u3Mi/lhkb3Ac7Gkr5kD+sche2schyT8+2R3c0/3o6XfCqs5K9R9rmIrwAUel74lDqP4UKp/vb0q1n5rs4HpOzU3pdx6NiDxDRD4sIp/2vdBf5+ffJiJ/JiKf899PGuMAp7pQUxN+6v0egvAQKrJVr89U6FPia+CXVfU5wPOBnxeR5wC/CnxIVZ8NfMhPj4KzXLCLFsvui7Oc16EI7+D6yPb5TIWde1LVL6nqx/3/DwGfwXW8fTmu5zn06IF+Xox14Yz/OyTG3P9hCR8qstLrMxX2Km0ReSbwfcA9wHeq6pf8oi8D39nxm7tE5KMi8tEHHjiffz30BRwTYxB/LuVVY3p9dkFEvltEPpl8vi0ir9/XavcuaRF5HPBe4PWq+u10maoqlJ+7qnq3qt6pqnc++cnnv7BzuJCVyM7PoTGHcoKmRXYIpVfVz6rq81T1ecC/AB4F3s+eVrsXC0VkiSP8O1T1fX72V0TkaX7504Cv9tnWRcY+hJ4L+eeAATuGp3gR8Peq+gX2tNp9ojeC62j7GVX9rWTRH+J6nkOPHuhDYggV28fLn4fA+/x2qPrFXFQeXD79yppeH/xgT8nnri2bfiXwB/7/XlY7oE+c/gXAq4G/EZFP+nm/DrwZeLeIvAb4AvBTPbbVij5chLeGDKXWlQj1jMjYhaGjXs7e9L6Rtw72FOAHaP0J4Nc29qeqIrL1JHaSXlX/CjrZ+aJdv9+GUgH3vRHGbrUdw5qEbY5J/n1Ufqqw7ggtsj8KfFxVv+KnvyIiT1PVL/Wx2vNqKpsJxvbil8nrjxSyfBWNtYE9rfbs0hD2sT9jqP1UhBzD7vRR+ekb7YZNQxCRW4CXAP8hmb2X1Z4d6VP06QF1FuIPUVncto3z5tOcJSdnnoR3GLKPrKo+ghvFLJ33AHtY7VmTHqZ7TWYfhe97s6TrbSPvlJXbQxHeRW+OQ4DsjSGJP0R/17P8tov8XcQfMgPzkHlIx+6C0LqM+1BpG/HHiuRclO6C26zNvoQfI9H5OARIgryAd1FsX+L3fUKUrM2U3QX72pwSgc9L+LGz+UP0Zk6Ylb2xnI/4Z8FZCW+2HIPdQbbzEL8v5kD4uJ+ZdSKZFemhuRDbimnMym2J8NsIXt7G5vr5jTBmr6lthJ+6V66qsD6Svh/6kH8K7Ev4OeNQ3dAvt71RjY/wvo1AXeQvqf15K7Spyg9N9nR7QfXPq/YlP1/0/Xtud0ibdfT0CdKCPWsr6HltTrrffSqulexet9Z+VEuJv4+v70v4vhg3H+hI+g30uQH6Kn6u9me5MYIq9yF3F7p+W6vFIDsruxu/27H+WRR+ioaxSx+n74N97c/QmMrDn4X4Q2HqFOdjnL4nushfCmueRc13WZvzqPz2/Zqt1ucsocu+Kn+IfH5VWNtDhyPamLxFti7M35aZUav2In5r+RkrtEPYmj4IxA9qv2+Fdldy2VkIX7ouQ+FobwpIC7x0A/Qh/lSJaXNDrvL7EH5MogccPX0PhAuRk79kd7Yp/mN1PMxtKp8T/pBkT6EzI/3kZqvvQ/wsF6YrwnHRRzvrc159yrWmf7nantvsty3p9ekDEblVRN4jIn8nIp8RkR8cbdybIWE7PjlKF6lOGrjCttq/udgE74tt55mXkVu/jL7X4qxQHby74FuBD6jq9wDPxY24N/y4N1Ohq7Cnfhw/1rCN8ONDqK3p9dm5JZEnAj+EG5IGVT1V1QfZc9yb6aM3iQBVhZu7K98+9/pp5bar4Wqbry9VjueAbVGWkp/f5uNLZN/dYLVjhTNgD0//FBH5aDJ9t6renUzfAXwN+H0ReS7wMeB1jDDuzWgIBVwiP+wOTZ6FuBZ78AFb++AsOTm7wpLbtjgG2WHv3Jtd494sgO8HXquq94jIW8msTJ9xb2Zx9Wttf1LkHjP3+SV/X6PR8wZ1vKhev3QecV66XqbwaRmVfPq2Mh8U6nx9n08P3Afcp6r3+On34G6CvYaYPEj05iyP2Pm8muxioRy33/2bOUZvVPXLwBdF5Lv9rBcBn+aijHvTVaDhLixZn9Tu1LT9vVtXtlqi0IAVLM7cfH04j2Bt+nQG2RalKal7aRtjQn1FdkC8FniHH9rv88DP4S75xR33ZqOlNSN/F/G3bnNHQ9Whff4u/953qL4uwh/6yTlkyo+qfhIo+f55jnujW94pWnU0tKTK35f4YfljLTWh5OWbZRSXpYTvDgmPW0Zza5GdjdKnBZ/fAH2Jn4cxH0vEzwm/KzS5y8pM9ZpLV0mdV/kfIE6/WQCV5L5UthLf/WY38VNc1FycXdamZGnOQvjSdRkKx4SzAvICr0Q3lD9vgAqq32fYkNa+to2dg07aEXxbJ5KzDtLUZWdyoo9J8hxzG5Z/ctKXQlMmu8C1Skv9U+UvkbyxMm21zy1OqvaHrrx2Ia3U5m0MubXp8vH5vD6EH6t3kyLYYyeStnWBzQI3aOvCpMofVN9QtjpdxE8xd5+/rSdUF+FLCp+SPSf6LpIP6flnJvSHsTd5gZZuglT9c+V368ykOXkgDNWVr9wYtZvwo1VsZ1iR3eeVmpWIfEJE/shP3yEi94jIvSLyLt9YsBOlVNLahzLrJKSZt9TVKu7TWsd9QlN6nA7bbUU4SuppZ/UeqFK8vm7d/I3Kd527W+7LyZeZ23ZWnoUyh8004EGgPT8TYR+xfB0udzngN4C3qOqzgG8Cr9l3512Fm6vOWZSpq6FmTm/e24b0OPPco20olUtefhvlO/Jbu1Wl12cq9H2P7O3AjwG/66cFeCEu4Qd65DBD0zhVVPZMYbqUP1WvVPGh3PK4S+3zdabGtqdNSeVbyzVsI6zffip2KXtJzUvXZQjLo4C10uszFfoq/W8Dv0JTvk8GHlTVtZ++D3h66Ycicld4L+iD3yjrVKmAc/sT5ycXNl2WEn9XE7tV3bgBxhpMdV/UaK/RDrYR3s3vKL/CU3XUhioFVPp9JkKflyf/OPBVVf3YWXagqner6p2qeuett1VbFWWb+ueqD5vET5F7+3lQen84gre9fI4uwpeUPZ1fKvehlR4GTS0eBH1fnvwTIvIy4BrwBFw/xVtFZOHV/nbg/l0bUtgYq9xIyChsN0alIUpwN4DxoUsXtpQY2mxCmhpDmSFdIc3Nsf4oqnixdVYRoJadQTu9fKry2wgP3U/MolCMNaT2zKpRO89SVX9NVW9X1WfiXk3+56r6M8CHgVf41XbmMAdsRgtM65OuU1L+sBw2Fd/N8/vpYXPc+vN4BvQ5jpKtgd2E7y7PzbLP1z8/+lViZ1eR7cAbgTeIyL04j/+2XT/QTGl2XYjSeqndgTbx4zrJPoPN6VOhnQvyCmxua4LCh7Ckm9cQvBQM6Fu+7afBQMo/s5DlXo1TqvoXwF/4/z8P/MC+OywVZEXZ4lg1LfsTbI9bJiC0rM4+2N5p3I46tN+2sSzPElLNCR/+T7/dts3GvHhMYxk9BR0wMiMi/wg8hNOytareKSK3Ae8Cngn8I/BTqvrNrm1Ma2l185Fq1bDSBTUm+wRlNxsXK318p4qfNmDlViCofTlPxba+p8K2/YYKbKrywbKlDU9dhC+pezO/XdYrXRSvy3AeX3p+euNHVPV5SSfy+Y57E3JvyiHK3Fe2H7/phQvz3e+2e89t/UEP2VBVyrDc+qbALYdaUni3vTbZ3XKzsTzf1tDRmwnszcuZ67g3INQtYjtUkihdstzS2Juod2pAMmVU2kIRp9vRmRDJsST9a300p0lQmza9uAlHbrbA5q0aqcrnCr8pBmaD6CWCN8cxov71J/SucW/C1v7UD/PxX/3y+Y5703XuocArsRuhyhJSr9/Mazx+a37YsYTMTEfuuWZaBvI3ncS7I1Gpwodp2B167Et2O4QRCI1T/bBr3BuAf6Wq94vIdwB/JiJ/19pdj3FvJs+yXGU9Wls9pNSN1A5ENa+18pVVR3Sn1hbUrZmS3yTx/UD00s1Tz7QXVYnw7eWNyufKHixgXC9T9xLRc1KP1TI7cMfw+/33V0Xk/bhgyldE5Gmq+qVZjnuTYyNuj8FiqNXEi1NSsHLjSrqdhAyU1XLuoUtoV16hnV7QOt8CuTcCAL5MQxnHfQzp30uw0u+zAyJyi4g8PvwP/Gvgb5nzuDeKsNJmlyah4aa/N+FHDtJ4+sbrW2qtQBL3myh8KQ8/9fVN/9rDtsz28fPQTrtI4/C5fy+R3W27PD/FIJYmw3azsRe+E3i/y3dkAfwPVf2AiHyEuY57o74iG4idF7DBNhdiC/lDDN+R28aLXEntLmhasU0qtSgYyVITZmJ1bKvxzH9nIUobSS4dT7+s0akH2btIPljFdsCGJ9829NzC/AeY77g34WI4yuV+u6a5ITYUPiW/pJEcCPoY3xnlG68qcQ1XXSMrpFEc6+sJZkL+u3aDpE9ssiy3YjaSvNz4FAifEruL7HH+lmjOcJBJMyj7YPqO4S2VYSMKE0OWqVJDm/zpvILi16Kt1tvaR3ZcspnD3Ib0C0gb0FKVh3Kl1U0n8fiE2LvIXmorGQUzqzodYISzTF3UxDQEoK3i0EH+Jsrj7Eyj+GmGZmpzgtrX0fM3sAe+AUr5/e2O3m1bk5I/V3i3/uYN4La5LTAwYq1mHjl9EdOSXmGlVasCC1BTxQpnuAnCUyCELNse3VJTNSHORPEh/N/YpEoatQ85OrU/ICPMJl4f0g4C4VsJZYWKa67wubqnZM+JnpK8PCTIkJ5+HuUbcJDREIoF6pOwXG58499NaLBKlV9pVYYr3DpuCBBDRY1VQ+WjOm5ZEvv2Fdrysbktjo2dg7Ym/5fDkJlqt0jchH1b66gpkn2MiE2KAaM3g2Bye3Oqm7t0pA2P4lApdWpeJfONV/Qa97LhGuvIr4v49KhQTr3HN4nNsSJR9dORAyq/grspLNUEhA9w5+CsTUiKa41qEPq5bqQJlxW+RPaU6DnJ66nszWUmPbQbVADfG8qHHFshS7+C2kb9vbVxlVPfQusVPVd8ty/TrJ/44G0pDmNjW1pxaz02W15LjXMbfn0L4XeRfW5jTo6FgzZOQbuBqt2ZubkRrG9DDMrfqP6CSiynwe54xQ9pC8jaPRlEm/yb4O01qScUMHZOfQlxHBuSFGlt+/hQfl0Kn5I9VfY6uQHa+TrtrMsxcMntTdufxkYmgndPiO77xIJXpkT5kUb1YwVWiYpfe+IHbFP3WhUOHLpMO3+XkIcW3VOgUe0uwpeUvYnzl8k+eKxe6ZViMCUOoPSNZ27llGsVY/bhYjlNby5ShQUBq4HIhuZVNU2mplHTtGqK5YoS1T719igs/fVIw5Zjpxen591qiQ0tsJmXt2o49dGqU19+QeGD8q+0KpLd+f/uFIVt6caD4TIrfeg5FdHZMJW3uNJS+EpCjKWm9jH7aGk8odt5PRJJvNxyeDW6dfnQ2JXw1tnhhiZxrFlvMyYfCJ9mXG6GLsdvoLrk9qat9CYpYEdkP58wzEeTWmA9ua2/wKFBy5F7UVB8r+g4D1+p+vh+2rqpMV6/lBDf1klSEVod1ZNugWkFNrUrK62cslNtVfic7F1Eb1maViv5UekHhZJHDhwqsVh/MxgfmbFis8apxr87G+AvjhCfAk3YslH8kJ5Qi9PA0NkkV/R64h5T6X5hMy6f5tgEgqYKn+bc1CobNwE0yp4Tv2mpTSq0Y6UgwGUnfVvpYwNTSniam8D4jMyg+MFrxxZadRf9iqyxvmArH/tGoMZFdJA1BkOlGr197Z8AwV9PaWtyxPqHj9q0IzZN1CZUYGs1rSeA++2muqffECq9eT7OuMQXPdqbdkU2CZ25eL2f71tQa5qLYoL6J+QPlgdCY1Si+GowuDi9UeNi9dI0WFnEk7+xM3ZHGHNo5JXYNH24RmLlNURm6uT/WKFNyF8iexfRSxEdGKmR6jJHb1xFNk2M8govGvNNnB2pMKI+9p5UbvP/JXQbdGSvMFj13QnFtmxO7W+RSjWmJzTH0U4znhJpCyw0ldeW/85szUoXSSXWKfw2wudkL4ctxyPm0EovIhXwUeB+Vf1xEbkDeCdu0LGPAa9W1dOu309KeotwYpOeU6E0AlHxIUWc1Vlp5RTc58UDLKk3VD9sI1xDF80JMf8m/AmhEcvVE5ah8uwbrlw6wDT0DykI0IQqQ5iyaYhyldZTraK9WemCWpt0jhPrjJnz8m2yr5KnQTuKk+fmj0z+4e3N63DvSniCnw7vSniniPwX3LsS/nPXj6dtcoQ41orVoFASL0ijWu35VpNHe6Je7QGMkt+FaEYrKlHILmz52cMYz/Zwg22bkYYlIVRC3fmVhu7L1T2WV1KOu8t84MGetPH1uz59MMS7EqZVenVKb9I4vEdILQ7ZkGEd4yMwRtQpvm+pdY1Q7ttKMyRIhZPNZcywbF+4kJC2pHaeXnyNOD1Oxs20zDMsY6hS217e0qj8DV3GKM2qpfxu3fBUSC1M3bo5GlvTRIOC1y+kFg9ZqR1WT34b966Ex/vp3u9KCDhIwlnw7NkS9xUbo0wzX0xMPW7+bzx9kzbsQptLQgtuaMByrbc1aQpyc8ErX4me2tOHzMpwHJv58tIKT1rawx7GEcw2noBlwneRPSf40BYnb4Pcgq2DPaXvShCRHz7r8UwevYlKn3E+JHc1St8ovmthVYzYlvIvpY5+36g6s6YVGN8YZXzEIrmmIQfHVXZrmkzG0LlkmlbZND6fNkgFb+4UfdEovV1SYzixS2okevkbdtkiu1P+NtHDOJ+R9PlIx2MrfX/sGuxpkHclTHpmiivwta1Y26rlHfP3JKUes8lFaauYe8xXrYpb6vWjb9XG77a60cXpMDbOtL4+DumXxOXddJI9GcomxuGbc1lpFQkfymIb4Zs6j8T5zTL3aa5Ne/S0c0F7fnZtZqB3JUzcXVA4tUkaQmJxzFalN17ptaX8C1NTqcWKz8g0/vfWR2mM9+5cAVwEq9IFCK4FU9bu4opQJ8ln4KzSGHYnRJvi28C1idpEhaeK3v1UK27olUjqE7tsEf6GXbZszNpWm0TfofTtQaOGz7KcoHHqjcA7ReQ/AZ9gx7sSJk9DKBVwqKCG/8GR1/g4pFVPfhSjgvXkx7oQ59I745Vtx/etGlb42L8Y6tD5BEMt4r5pktFqpm2garXEEqIpaSOUiZGa2PiUPNXSqEtToTVbyR5usDAdMGoHkhFIf553JUzfXdBmnUhCBw/akZtG6d33QoQ1OHWP5HdRmxi9McRtVf6mWJp1jOhgmvaAlfrojgiVihtZzh/T2Eln6ViVFhKVr6KXX2nFDbtM5lWs7IITXVCr4cQudpJ97Z+q61zhC5XadPnguNRpCApru1mNyAnuBl1qKq9AtDDWSrQ6S1O7nlDGZ1Natw2jPrRpvNUxfvvapBzXiGvkUtezqkZYsjkUx1iIDVNZjD2Q/9QnjqWEd99hvlP6tVaR4BZhFfx4q17URGlSpXfzNlV/SAh7RW8mweQhy1JSkyavjJTsBjBiIpHDEyDcEI7kweZ4T68+s1J8RAeiyhu1VIRsROftT6lcvF7XMVFtKoR8m5VWnFJFL3+aqntL8ZvKaiD5SR3Ox1dEs0a/lOg5wdOXm42SUgxTefq90Iv0InIrrgXse3EPq38HfJY93vMDrpBX9Wb1UFoV2oz0KBKILu0ngDWCEZeSAMSGK/f7JK3BEL1LZdw/IWFrpZXP5JTGboycimD9HiywwrSsTUr4E7vkhl1yogtWNlRk3atygr0J34HYa2s2iG7Vvb0vr7im6j7q2/0uIulxsdAPqOorROQKcDPw67j3/LxZRH4V956fN27biLJpb5ynbyqxtQqStLy6ELqgqr4CG5S+sTMharNS4/2+oaJR+tByG4gFtDpY19J+RQ2MR/zY2UU1C6uaTOEXxZBktDM+lBn8uiO7szvbyB7mxeOZojJ70UgvIk8Efgj4twA+e+1URF4O/LBf7e24mvR20qtwWpdaYzdtTSj+MB0UOrU3a+OUf2Fdh5OFrwvEVlYV1tKQYGnqOHjU0tbUIqxMBRZWcuqHFqknsTg1zpLUPux4Q5fexriUg5VWPGqvsLJO3U/skpUaTq1POPMKH5Q92JrabtqZ2rYrsqH0SzaHbNkQuIj25g7ga8Dvi8hzcambr6Pne35E5C7gLoArT30Cqs1QdW55qLxKJGyYn4Yysca3voYfunnRy6cVWSr3VCD4/MrVCVRZ2YVrqfVpCitdxDeOu8/4rbIbldgQmqTJjV/ZRUvhV77hKFRUc2VfWxNVvaXy4X+/721KP5rFuYCkXwDfD7xWVe8RkbeSvbJw23t+fO7E3QA3Peuf6em63Dglhf+NtJ8A4UYQGsUXUWq1LeUP1siqsEhaG0Oq7VIrllKzkoprsgIDN3RJhbJSF6sfq1Jbq2LV7Welxim8V/lH7FWn8LX7fri+ytpWXK+XkfSnvsU0V/ZI+oTkKfHdtRBCDzPtIPsojVMXMHpzH3Cfqt7jp9+DI/1e7/kJSO+MqOrQ8vXqfX2tzZA0Km6eeqKHi9qq+PoojhGv7BJaZ53Px+LsDMQb4FQrjIZ+tNP28InDe/gQ5Sr7rK1TeEdqR/QQtQkkb5G/QPYS0cN3S+XHPNGLpvSq+mUR+aKIfLeqfhb3xodP+8/PAm+mR74D+Irsul05zB8QgXdhfvymUXun+I0FEm99Kj9d+8zLtRoWYlmr4UQWXK3WWIQrZk0llqXUXJMVFsMj5gSjlhWnSWeS4RE6j9zwKv+IXuFRe5VH7FUerq9FhT+1C67XS07qRUvdV7WPx9c+zyazMrVtiK3atjY56fNUo7HszUX09ACvBd7hIzefB34OFwjs/Z4fAFSwNnh2Pyvz9+FCNORvPL4C4tXdKX77CYCxUb0qY6EGG5pWqzXGVizEWaEbdhltjVFNoiZmNMIHhF5SIUqTxuFv2CWn1ocoE3U/tRW1dRXW2hqn9NAovFd0m1RalYb8rqzx39t9/OB5dxeR9Kr6SaCU8tn7PT8Bti4kdkYpaN8Q7n+N84LK59OxA4qxGIG1sVTGJaZVxlIbw6mtuGJqp/rGN+yYhfP2WnFNXELareY6p1pzdSR5qlFOVbmhFQ/Zm/h2fY2H7DUera/yrfVNnNoFD6+vcGoX3FgvI9nDd+2tjCN7Q/Lak9t93L5K0wEbxB6zEnsRST8UVDsqSqm3x3dZ9aTTqPhuDfEKD47kIGhMMwgVWKeAi6qOZFhSN40ylbAUS63Co+YKtRpuVFe4YmuvuE0W5NAIldgbXuFv6BUera9y4u3M2lbcqJec1gtu1D6CU1esQrSmriLZW62sNny7m0C1YGMKIcqx+ShcXHszEAS7Nhs+PkDdKnHdqPhZ7F7EtW1bn1BlTKP0Io78AtTWV4gry8oalsbG9GarwhVTU4lyxay5ZlbYSniwupmbzYprsffZsLihlod0wUP2Jh6oH8fD9TW+ub6ZU7vgodU1Tm3Fo+sr1NZwUi9YeVVf1a7iWntbY61xIhLsjE0InlqaWLgFH7+FjEP6+8tNegWtBU0JnSMpoeYmkOjpm+lmXfXhieDfjXE3Vu2/rbUYY1kZZW0NC+/9F77B66ZqxVWzplbDbYuHeby9zuNlLNLDQ/YKD9Y38631zXyrvomHVte4Xi/59uoaa2u4sV6ytobTuqK2Qd29b7fNtzv3LBKT2Jk4UERclh1MB7GPnn5IKISBf1rlIIWJXB6kPV+lmVWLvwnUxfDVuieBMc4KWSMYH8OvrVAZ/8pN4xLYTm3FwtSc2IVTfJ/LXsnDGJpkhLO8jC32jsJVYL+wfgJfXj+R/7d6El9dPZ5H1lf51uoap3bBIyun8DfWixbZrZVE2SUqufqC0MQ3q0pTuB1RmiLZxyTmpSY9wDqtpW5ZLyOYppXdRtT8LG0pP76Sa41fRoUYizHKylTue11hjOW0di95u1EvuVatuF4v+cryCXz5yhP52vIbLGXNLeYEaPrXbrwGtIB8VGDX8LTgC6dP4SurJ/Dg6iYeOLmFG/WSb59co1bh+ukyJuVZK9S1I7p6S4MKajeJHSyNm6b93SqowrI+888DHc7eiMg14C+BqzjuvkdV3zTrwZ5IL06Y7kKX0kNi7pv4pqbriGvVcooviFFEK9Qq1vt+9y6GClWn/EaUk8oVx6ldsLauBTfE8kOH9H0R+quGfJovnzyRb6xu4aHVVb51ehMn6wXXV0tqK5yuF67PwbpqyO5JrVZi+W0lecnKTBWp6cJwN9MJ8EJVfVhElsBficifAG9gj8GeJq/IylryWZ3rBmhJKiRZJ99Gy//7axwyNr3lWZkFIsqNaokR5eHFVYyxfH1xC4vKcmWx5ublCoOyrFwvqzTdeRfyNN7QqPToasnpesHpumJdG6w11Gsfa699BdR6YidE7yQ32fz8/1gkPYg+kg0ZKg1BXceLh/3k0n8UN9jTT/v5bwf+I7MhfSlmu6ugpeOCpQ+MTv+frS84QoUngYBY67oKWlfpXa8rqspyo1pw/XTp2gFMk+ezL9KW0mBd6tpQ1wbrIzK2Fl/JN00ZBaJ3PR2TQ9konyntSw/sUWxbx70BwjiWHwOeBfwO8PfMfbCnDaUfYpuJzy+hKHLRGrlqahhYtxY4TeoGgyadJTe9FNRbulScHcSZWUWxhf0ap3aNe4Oq1sDzfMem9wPfs+8hTa70h8i4K/N2mzrmNeUhD6bLqj2GMcJNqaoPisiHgR9kz8Gepif9vkrfo8A2tljwt5KqaogoqMSbUGriTSlu8DNMna2/xzG17RfxqWErQStQ4z4IhCH71Ueb0vXj7/MT7elouo5pKgzZIisiTwVWnvA3AS/BjVj8YdxgT+9kboM99SqAPgWUW/guQqZkVfFEp3niBJIryNr9b9ae8GvFrNrrnEmxxNtzT3JdgF044ttFMy+sg0i8GcJvWykZXcTfqMwXjiU//oluArEDsR6eBrzd+3oDvFtV/0hEPs1cB3vCk6vvuiW0CS7teSWPnJKchORexc3KTZtTRWpYnDiyL04s1XWL1Iqp1d9AZ6jIiiQKL9Q3GdZXDXYJ66ue/Fcc0e2CqPzxJgDUSOsmcCfS/k7bLLLi2cSUin9WsShtSvWvge8rzJ/vYE8AUvcr8Y0nQsd0S4EL8yRX9UD2tVMgs3LzqxPFrGH5qKU6sSweWVN96wZiLaz8nepJ30e5NKQ0i286XlSoMdRPvMb6lgX1VYO52WAXULsXm1Nf8TdISv5ggYxXfdgkv5+vYWEevMqt0MQ259Ln3mxUZPsoer5eyacn03EfNiF94tVFE7KfKmJhcV0xK2X58JrFIyvMIyfItx6GukbrMJ6+3XbIHSfi5Foqg1QVYi3m9CrrW5agC+wVZ7u0GXncfUvj/TEN6Ynq778T8ks23RxDfkxdx7rPie2BS016ttuboiJ0+feShcm/g8LXTp0j6S1Uq8bOSK0sHqmdwn/rOubbj6IPPcz6gW+c82zzExSqG0+ievzjkJObkfom6qsGqSvsQpDaqXtdu3TpWOGtpKX4UemT7+j/YdPyxP23vsrrjIBLrfRSUvod60d02JmWlUl8exOJ0VZEJpDeKb1Td6mVamWRtUVWNaxrqNsvYxsEqm676xpZ1cjaUhmwKxNVXlWin7f+fFXVV4R9bD9R/PidRXw2/vfllZN8EkJeZtKDDwN2YZuF8fPKoce80qrtymodSK/Ry1eh4nq9dqR/ZIWcrpHrJ+j16+jparBzbp3iag3XryOLimpZoeslKoJdOBOvC0VqQY246I5XeVsBRqPN2ajcZqoP7WiPpiq/zfoMfsL7Cd0UOGwaQqmdpod330Z2oInQ2ORTa1R+ExS/VneDrC3UTvGpLQwXYtuE6z4GtYsM6do9YYwYbO1k3Bg3jr1YIUi7gSZFWEIh4PsI01L6QHYJTwWyb78dTf6PGPgmGDJOPxSmV/rU02/x8F2V1EbhdZP8acVVwdTaIr2pcXbGKubEYmrF3Kgxa4s5WbkozWqFrtfoGPYG3HbXa2S1gtOVG3pw4Qhf+Ri9qMEYwS6dwqeNWbaS6O1dxVY3lD6ESVOVL3n9YveAMZT/QG9u7MIBOpFszu6K1JQVvk12yDx8HqLMVF6sRl/v/rdgrVdeG0dQHhuqbt+6ULDWtwC79xq6FyP6Yw4E9pCEzASL72+CGAVSbZE/Kr739Ol0PJ7E9w+NS630opnSp9B0PW3Ny+PuG34+kl3b057splbf4qqYtSJrxZw6L29Oa1hbV8Fcrd13XY9ncaw2+6hrWAlyUiELRSuDqdx+deEYLQZ0oa5xK4YzkyiOkfZ0y+aEQttUd836J4xm7XNLOwNM33NKE1IXlgWkKt8VsWmTXrdHbazz7LLW5n+v8mKdyvshBMY68+xcva83br8ajgHjbt41iPFyXKed5Ml8utJKpU7Dl/lvWmrf3BA5KbXoe86OS1+RNetNUhVzZwoWZ4Psweq0bE1DdrP2y30qgawVs7Je8S34SiTr2lucoPKW+N7awcvAtvdVGRfCVIOsrDs3I7FFVxduECej3tuH0GaSp6OiidJL5u+JSr7h5RNua0v4h73xLzfpKRRAVr5djU9t8jdkD9tM7U3MlbHqRjJOFT74+LVTeILCh2+ro/t69R3Y436t9edhUfdiLSfMXvHDOyVEXYqCqO+ZFcx5KsxGvapLm+yFuH2b6DCKx1Eud0U2xs4L8yM6KrExHJncBFHZ/TxT+3XW3u4EW7NWTB3CktZ7fevfZ9lYm5huMAWs259UvkIL7viwmLVFK2dppHJDebuQpXvfloo0b0mvfEU2GfkgZGZGUvunhmRklw6lHxqXuiILjSpvIJK9UIn10217o62UX1GNZI+hy7X1FVmLrJ2Kx5j82jaED3ZDffTGahxLZ2io9aMuqyJqo50S8LZKUGMRFRe5UXCjaxqn4pgk/EK8CUjKTUJXyODt/bLOymyhIjvoTXAkfUdFtkvho9KH+Y2yxwpsHdZpbE0apnSK7m1NUPnQSOQrrxoqsXYitffH4MivjdWpDKKK1iDGou4NEj6nXtwxpzalFgSfp9MaDQIQN9KbpqkKrW9tT6cYqDJ7bJzaoyLbkL8hOSQ3TbQ3GdkV599VI+El2hrb2Jp13ZA8tMKGm2AqhH0GWxVuuLVFKmneruxdvVR+kCcjMec+5uKoHwUuqdCGMT5DmDOMAN3H3gxWmQ3XYQCIyDOA/457643iOo6/VURuY4+X/hWGEB4P4a6PH9sQdzPeHiqetDx8K2LjK6ZxeVB0vw4hzSAUfKwDNArfqP3EIYZoR7KnjrXN8YVzTG7gpoKusRyjvWuVV2b/bFauWih7m12foe5/7fnZjTXwy6r6HOD5wM+LyHNwLwn5kKo+G/gQ2ZtyckyfT1/rDk+fTyeqHraR+HZU4zbjDVBrM7+2TX6N9SHKOonLRz+fEH/sGyDdj0rM6JTaggjqQ5huHf9txIU3cesYdYUVG6eqYOI15uPHboaJ+jdRGt1IT2gd4oCefqibx7/j7Ev+/4dE5DO44T5eDvywX+3t7Hjp3+FDlnTZG43TqXq1W2Pbateskyi8Xy8oJ6nK51bG6rSe3mr7rZ1u2LX28Yo/Hw3zBcWdl/o8hNBGpbgktdivNhDdx/aDZfEB04brWrA3Q6p8f3uzc9ybABF5Jq7r4D30fOlfwMRKrzGs2Mxr/k3VHBKS+3kp0cPjvfHy7vHfLPPqHr5jmNLGSiS1U/kYsUlvgjEbpzzT1d+cLg+nBuNVPb6tRXBvzLXoAiQ1417ppRIfvQnRGkW8+ruYpyO/1CRqr5Hx4UEyWpwe+loX6DHuDYCIPA54L/B6Vf22JJWTbS/9CzhIGkJAS/XT+RnZ3Tw2CN/E7jO7461Nqu6pjw/eNsLqeCTfBQ2KH1iozbc/blXceRjv5Y1/KOAGdBXTPAExbgxPMf6JkFVJYzJaUHt/M8TDGaGWN2T0xo9h+V7gHar6Pj97r5f+9SK9iPwS8O9xxfo3uHdOPY09RooFd/Kt6E1eGBu+PvH/qbLHTEuNlb1IZtuQWmLjU/DuiZdPG6QC8ezEWZbWenm2LkRZOY8vxnv78L4tDFB7Rlr3JDD+48e+xPeljYofNu2VH2jUPy3k3NaMcO8PGL0R3PAen1HV30oW/SF7vPRv530tIk8HfhG4U1W/F/dsfiVukJ23qOqzgG/iRordjqBY4RPtSUfkJqh6rRu+PtqZGH0huTES/x6/C14++PeM6GM1THVuv3Q8eXQpO59QNqGZP4Zn0/IJZRujWuUIThrFaV2fIW4A3eOzGy8AXg28UEQ+6T8vw5H9JSLyOeDFfroTfe3NArhJRFbAzbga9F4jxQakd/3GYy9NNYDm4qX/t1IPUsXPFD7eFJmPjw1B7YjNVAqfo8nD8bXNkISGP7+QlwPR37vKrPVxd+ssjQ/eBK8fhiiPOfV+F+3kM21SFMLxDOzrXZh6mLJV1b+iu+bR+6V/fd4je7+I/CbwT8B14E9xdubBPiPFishdwF0AV6/d6giZLk8nNZ2v7QqtV51406R2JtwAQems9banTXRJksr8ubWmW/+PjbReYdUN86HqiOynpbZoZZK6iG+hNYpgmpCj4p6/io/SuEhPdPQGV4ENfl6kZSED4vYGP9cxNnp27CS9iDwJFwe9A3gQ+J/AS/vuwIec7gZ4/BNu13LIMiGaZvNSZXcbbKu70qhhEuqTlNBp5TW1OXGfdrpQZY6YxhwGtEzsjPjojmhzI6i6FzSIe0JosDfJTRDTiRO3JEL0/OX+srL55B0IQyn9UOhjb14M/IOqfg1ARN6H81Z7jRQLNOmyYXojTu6/Nbc5ueq34+/x6ZErvFd2ye1NEqZsKfvUxLcWKk92fyxqcOHLqor2JnIzDGPoQ5lB8YHG7gR7E8KW0HQxjOkIyTHEkKV/+g3cgWQPvz4Z+gSo/gl4vojc7GvPLwI+DXwYN1Is9KgxA6ThRFe5ovmEVtSswhUqa6JJWnBiZ5oKni0qfLvimv+fRG4OidKx5LYrf3KlldvwfyhDpWmNTir5sXzTssmuQ7g2EvZxbjT73PWZCn08/T0i8h7g47jch0/g7Mr/Zo+RYgOkkHDWthp+vZSc0FL2FrkhiVqEC9moeh6ibBTfZvsNBJtI7dXijXx2DL5CG8KZIfMy+HuffgyJ4oOr4IKr5IoLYbaUH2KldTP9QDP5Z1gffgHtDar6JuBN2ezPs8dIsW5DlAsgmVUkO2wSPswjm07tk22vF7bRitTkLbEHgIvZk9j6ENHxZPTfYtXl2EDL44dVAbde8P0hy9KzO4RKJQSHWtXWkcpAx4n9nwcHyL1pkzrOTyczUm8QPSd/8OKpwoc8mujdbdP4lFqngLRhapKEs4rYQJX6ekPj7WMdo0k0E/CjmwlpfDF6e2sTlXe1gaj84H5n/U0UfXx2fGZoX38BlX5QbIvTl+xGB+E312vszkaaQX4TJcehwfIcEtaixvjeU7pZ01IlehH/VNBKmhvXNE8Dtz2iyreVn/bTIy5PdiUMX8eZF+enTziTrn6oyezcvqREb3n5QoW1XUnriNj0SSEeS51Scm4ss857qAVrXERG1T0JfNpxaLiS2jrFDx5fXFpCTL5SSRqgvPJbYgtWE6Vp32RjxOnl0KKS4QBKX54t25T5LIRvqXqy01CBTfJt3OZ0eIXbhTTRzLoKa8yxJzH5qQVKziukDDfbwxE45NWDI3eq/OG3yf9qx6B62DgXr3FqSLhQZEcJtEia/qZA9PB/sDN5d7t0OumLWvTy2b7HzrtJ9yOVr7QG9Q9WJff24NITTJJV5ltq/dCuScttEojfUP6kNSq5CaTryTMAhKFCn8PhAKnFhQIoWZt03YJvJ78Zun5DouKpnZlB1CZHHsUJNic+DdxKrYhOsQwCgvJ7xJsLXNJmh8UZHDMr54PYm+5h/UoV2YK/zyusYVl8TU5Q9CQW3+qLmtqdeVRkYwcSaJNdLarSDBESO5pAVHxouhqmlVVwjU9B4bMKb34TpBi0ZfZSkz6E6LJ5LdgC8bNv2aX4CeE3vHqaVWnbyj995/BMyWMURxpv721PtDkbN0jm1T3JNY/UiLRvAGhugvjjJAQ6FFEvu6cHMu9eKNgtlkUKfn7D50PBx9v2eofMt+nCRrw+8fb4aE6J+OFcMnsSk81K4c80tNleOHyMnvlFbyYdAgTF90u17XBiqHDGxqR2GkEcfzJbFserKWyjSHgfsdF0O3Bwb6+tGzhUwpvjTaNN2llmulmmpfJLzz3fRijT9HP+s2u2veuzAyLyeyLyVRH522TebSLyZyLyOf/9pF3bmZb0fQsANhOj8k9J9f30ZC2rI2HXTaiFc96o7Bc+rcSunGxnJOLuk9my7f339d/YTGvfa8wbOMjLk3uELPPp9IKm0yU7E9ZLE8r6qPwhO4aHlAR3QPG4xALG+34DaI3z+SZ69lgZNaG3uHRXXoOgbIz3MWKcHgbz9Kr6l37ojxQvZ48xb+AQnh6239UlskOZ8B19XPPQZCehS15z0izLqj0vePW8sgqNzy9uS5vfpARPY/+hjLIUhPj7HAPeCCPH6fca8wYO8s6pLFxYQimCk/6fkX2nwvv/WxGbPHIDTcPU2B7fEzIM1dEidKjQhkhOUHvoVnySBqzYMJUpP2yP3ARsuxHOc7790Huwp/Judo95A5Mr/RbvtqWVtK34BcJHsm4SOfyvXdubIwJZSW1OQXlDyBNcZKel8Inyw6bCp+WTbnvosgmV437oNdhThr3GvIFDKP2uAsgLvWBhWsoODdlL81PCp41VZDfMoRFi8anSeuJu+HtIQhDNU0ATDy8izVvPU/UP207n56PODY1xRWavMW/gUJ4etih+d+NV3vnDzSwQPrE0xe0WBnWaKudmF2KOfYzDZzH5ltVJK6nWZWh6qxTOT1KFz8s2rzcEDN5PdpiyFZE/wFVanyIi9+E6Nr0ZeLeIvAb4AvBTu7Yzvb3Z9lLiQuFsqDq0bUyX6ucKH9Y9ZOW1AOfrbdOdKSCQ3VsV9QSNVke3qH5N3F7aOLVB5bQOMBaUwZ6mqvqqjkW9x7yBQ7fIbiwqLOtL9nR5SeG7tj2XWH4gce7dc48OSUQmUXdoyB+UH5obA4rjVI5MeVw9biZl7DG5p+9seOlSgxLR8/UzshcV3k8XnxxzQurt8zBkrvjQ2J0S+d0Gk383Wa8jpR4kOxioZXc4HLS7ILBdBTYqtR2q30X2dBt98j8OkXCW25qwKFRGC8QHuskPbXKnm9cOOxMcZ8exnBszi5YdIGTZs0U2YNtNkiwrPkEywndFfQ5diY3x+kJltEj89LdxnTT+XlD/uGxbbH6km/5ykx5fmexRCKULkP2uHc3Zru47CZ9HgaaCJ3ov4kM7qhM2kURh2sNt1G3r0qnohXMezPJsaZs5EA7QItvDu6corL+h6tsqq6X154CCanevqptd+loNTU0ltxWqhPYTIP42K+uSrRnq6af0s5YTYl4hy3ztYsSlO46fL98a1+fwtiZHS+0hCUkmcff0/NMW1yz8qGF5QJ3dDC30vyZnwsxE53CNUxl2qnExvr4n4dOfzozwAZH4KZLcnJbql3pRpaTOLRHlch6zY7gTusus9Aq6h9K73/RsuaVkezqiPa0fzeCCpBXPrvmJTUnPs9grKW9pzdMR8t3ve7z7QEHnUMYJZqP0wO7HYIc33Ep2mLWlKaEZczKJqoQbIvPoGx4e2uW0pT+t/+Egx7wVMyvzA0Vv9r/zt9qfYl/bwpNgVxvB5JGbzI5kcfuW1UmPNQ9FUm5tBaCuz2ZfuvJyzoJL7el3EL53lGWbcvQhe2m9Q12YHsQH2j6/ZIfyc0zW31aunTfEUBGXM4rcmJiU9MoZ0hA6N7a9IDttzMz8ZREFUqfnU1T/FIUnQeeuxk5DgEuu9LC3SvfBTp9+iIFa90Ge3x7nZ5Ymzi5EYLbF4bdhjxvkbND9gxcjYz5pCF2/OMsF6bOPOZA9Rxf5YWeDUp9y2giFlrY7NLY1SB4I04csz5qCcK79zqvQd6JPlOUMZaSp4I6VXFbc8bws5QGiNyMWwEUj9z7oOrezhhwnIqIyvzDxrDqRHHEGzL08dX9LOzbm1Th1xGMSc6vIypQZiCLyNeAR4OuT7fR8eAoX51hhnOP956r61LP+WEQ+gDuuPvi6qvZ+G/1ZMSnpAUTko2cY2+QguEjHChfveA+FiQdwPeKIw+NI+iMuHQ5B+t5jE84AF+lY4eId70Ewuac/4ohD42hvjrh0OJL+iEuHyUgvIi8Vkc+KyL0isvMVKVNDRJ4hIh8WkU+LyKdE5HV+/t7vNJoKIlKJyCdE5I/89B0ico8v43eJyJVDH+McMQnpRaQCfgf4UeA5wKtE5DlT7HsPrIFfVtXnAM8Hft4f497vNJoQrwM+k0z/BvAWVX0W8E3gNQc5qpljKqX/AeBeVf28qp4C78S9K2g2UNUvqerH/f8P4cj0dNxxvt2v9nbgJw9ygBlE5Hbgx4Df9dMCvBB4j19lNsc6N0xF+qcDX0ym7/PzZgn/Mq/vA+7hDO80mgi/DfwKzQitTwYeVNW1n551GR8Sx4psBhF5HPBe4PWq+u10mbr47sFjvCLy48BXVfVjhz6Wi4ipsizvB56RTN/u580KIrLEEf4dqvo+P3vvdxpNgBcAPyEiLwOuAU8A3grcKiILr/azLOM5YCql/wjwbB9duAK8EveuoNnAe+K3AZ9R1d9KFoV3GkHPdxqNDVX9NVW9XVWfiSvLP1fVnwE+DLzCrzaLY50jJiG9V55fAD6IqyC+W1U/NcW+98ALgFcDLxSRT/rPy3DvNHqJiHwOeLGfniveCLxBRO7Fefy3Hfh4ZoljGsIRlw7HiuwRlw5H0h9x6XAk/RGXDkfSH3HpcCT9EZcOR9IfcelwJP0Rlw7/HxAHoDMGQn+fAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(T.reshape(nz, nx))\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "markdown", + "id": "e7367fc3", + "metadata": {}, + "source": [ + "The temperature in the region is plotted above. As we can see only the floating part is heated up, while the bottom part is cooled down. We also plotted the gradient below, which has large negative values between the two components, indicating tendency to connect the structure in order to lower the temperature." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "fb8073ec", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMkAAAEDCAYAAACfwiBjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAfgklEQVR4nO2de5QkZ3mfn19V98zsRdJKu0JIu2ukg4Sx4hgwawHmcImAEwX7GOMDGBsTkShRzrFxIMYxOJwT+yT5A+KES46Jk41wLAyxUIRwFFtBEBDBTmDRBSGQ1iBZ3FbWBd1Xe5ueqjd/VNVMTU93V/V0V3VVz/vMqTNd1VVffd39/ep93+/76i2ZGY7jDCeYdQUcp+m4SBynABeJ4xTgInGcAlwkjlOAi8RxCnCROBMj6Q8lPSzpmyX3f5OkuyXdJem/VV2/SZGPkziTIunlwNPAx8zsxwv2vQi4FrjUzB6X9Awze7iOem4WtyTOxJjZl4DH8tskPVvSZyTdJukvJD03fesfAx8xs8fTYxstEHCRONVxEPh1M3sh8JvAf0y3Pwd4jqT/K+krki6bWQ1L0pl1BZz5Q9JO4KeB/y4p27yY/u8AFwGvBPYBX5L0t83siZqrWRoXiVMFAfCEmT1/wHtHgENm1gO+I+nbJKK5pcb6jYW7W87UMbOnSATwRgAlPC99+09JrAiS9pC4X/fNoJqlcZE4EyPpT4AvAz8q6YikK4C3AFdI+jpwF/C6dPebgEcl3Q3cDPxzM3t0FvUuy0RdwGnQ9WEgBK4ys/dNq2KO0xQ2LRJJIfBt4DUkfuYtwC+Z2d3Tq57jzJ5JAvdLgHvN7D4ASdeQmNShItlzVmjP2u99BW3i9juXHzGzszd7/N/9Ozvs0ceiUvveduepm8xsZJdwkfciaRH4GPBC4FHgF83su5uo+iqTtNi9wA9y60eAF/XvJOlK4EqA/XtD/t9n9k5wSqduls77zvcmOf7RxyK+etOPlNo3PPeePaPeT72Xj5DzXiTd0Oe9XAE8bmYXSnoz8H7gFzdV+ZTKA3czO2hmB8zswNm7w6pP5zQMA+KSfyVY9V7MbBnIvJc8rwOuTl9fB7xKucGazTCJJbkf2J9b35duc5xVDKNn5dwtYI+kW3PrB83sYG69jPeyuo+ZrUh6EtgNPDJWxXNMIpJbgIskXUAijjcDvzxBec6cUtJKADxiZgeqrMtm2LRIUpW+naTfOwT+0MzumlrNnLnAMKLpzTQv471k+xyR1AHOIAngN81EXU1mdiNw4yRlOPNPzNREUsZ7uQG4nGRw8w3AF2zC+0G8P9apFAOiKYlkmPci6V8Bt5rZDcBHgT+WdC/J9P03T3peF4lTOVO0JAO9FzP7l7nXJ4E3Tu2EuEicijGg1/K7X10kTqUYNjV3a1a4SJxqMYjarREXiVMtyYh7u3GROBUjIiaaFTJzXCROpSSBu4vEcYaSjJO4SBxnJLFbEscZjlsSxynAEFHL8424SJzKcXfLcUZgiGVr9x2pLhKnUpLBRHe3HGckHrg7zgjMRGRuSRxnJLFbEscZThK4t7uZtbv2TuPxwN1xShC1fJyk3RJ3Gk824l5mmQRJZ0n6nKR70v9nDtjn+ZK+nD71905JpdKfukicyoktKLVMyHuAz5vZRcDn0/V+jgN/38z+FnAZ8CFJu4oKdpE4lZJMcKzekrA+B/DVwM9vqIvZt83snvT13wAPA4UZ8z0mcSrFEL3y01KKcgGP4hwzeyB9/SBwzqidJV0CLAB/XVSwi8SpFDPGGUwcmQtY0v8GnjngrfeuP6eZpKHpJySdC/wxcLmZFd6C7yJxKkZTG0w0s1cPPYv0kKRzzeyBVAQPD9nvdODPgfea2VfKnNdjEqdSjMSSlFkmJMsBTPr/f/TvIGkB+DTwMTO7rmzBLhKncmoK3N8HvEbSPcCr03UkHZB0VbrPm4CXA2+TdEe6PL+oYHe3nEoxVMtNV+ljrl81YPutwD9KX38c+Pi4ZbtInEpJUgq1u5m1u/ZOC/DkdI4zEoNpjKbPlFpFkj2JFSDwPoPGMsYzDkvhlmST9P8QLprZMW1R5DFT6y1JYe0l7Zd0s6S709mT70i3F866HIcxnuXtTIk6vvMkcA9LLU2ljMRXgHeZ2cXAi4Ffk3Qx5WZdjo2LpXrq/Y5V12BiZRTWzMweMLPb09dHgcMkD5QvnHU5CS6Uaqj7e00Cd5VamspYMYmk84EXAIcoOetS0pXAlQD7945nUmNij1WmyKwuPG1Pc1q69pJ2Ap8C3mlmT+XfS5+TPXDWpZkdNLMDZnZgz+7x/U63KNNhVt9jNuI+95ZEUpdEIJ8ws+vTzaVmXTpO2xNBlOndEskD5A+b2QdybxXOupwWbk0mY5bfnxn04qDU0lTKWJKXAm8FviHpjnTbvyCZZXmtpCuA75HMsHScdSTuVnMFUIZCkZjZX8LQIdMNsy6rwoP4zdEEK+wj7o4zgqwLuM20SiRuTcajCVaEreBuOc6keMJsxxlB0rvV3HlZZWidSNzlKkczXK36bt+tEm9tTuXEaVqhomUSxpmVLul0SUck/X6Zsl0kTqXUOMFxnFnp/xr4UtmCXSRO5dSUMLvUrHRJLySZjPvZsgW3LiZx2oWZWCkvgEpzAUsKgH8P/ApJbq5SuEicyhnDlao6F/CvAjea2ZFkSmI5WimSqnu4giFBZDz4boDayipLU3q2YLoj7lPIBfwS4GWSfhXYCSxIetrMRt5V20qRVMmwRp29N07jnmZZbaamLuBsVvr7GDIr3czekr2W9DbgQJFAwAP3dfQ36lABoYKR+9RRVpup8aarMrmAN4VbkpRBjboKQgVEuUdibAWLUse0lDK5gPu2/xHwR2XKdpEMoL8hb6YRbzgmV15/+fOMGaw0+IaqMrhIWG9F8hZkmlf4rKzsXHmhzLs1afu0FBdJH5HFlbla/efJmGeBzMPcLRcJSSPNW5PI4kob7lZxtTLMReKMQ78gtwJ+P4njjMDMYxLHKUBE3rs1H8zKDZp27NOkKSkZHpPMCL9DsR14tpQ5o25rMs9dv6tYEpe0GReJUzneuzVn1GVNtoQVIRlM9MB9Dqm6AW8VgWS4u+U4BXjvluOMwMxF4jiFeBew4xTgMckM8QHF9TRytB0Re++W44ym5YbEL8NOxaSBe5llEsrmApb0I5I+K+mwpLvTx66PZJxHVIeSvibpz9L1CyQdknSvpE9KWij9iaZIE12MWdDo78FKLpNRNhfwx4DfM7MfAy6hxFOjx7Ek7wAO59bfD3zQzC4EHgeuGKMsZwtRhyWhRC5gSRcDHTP7XFIve9rMjhcVXEokkvYBPwNcla4LuBS4blSlHMeAOFaphTQXcG65coxTFeYCBp4DPCHp+tQr+j1JhU8YKhu4fwj4LeC0dH038ISZraTrR4C9gw5MP+iVAPv3tvuJR84mMKA5uYA7wMuAFwDfBz4JvA346KhKFYpE0s8CD5vZbZJeWbR/P2lW8IMAP/m8xbZ3dDibYFrjJFPIBXwEuMPM7kuP+VPgxRSIpIy79VLg5yR9F7iGxM36MLBLUiayfcD9JcpytiL1BO5ZLmAYkgsYuIWk3Z6drl8K3F1UcKFIzOy3zWyfmZ0PvBn4Qpp4+GbgDQWVqoVG9+zUQLM/f7mgfQqBe2EuYDOLgN8EPi/pG4CA/1JU8CSDie8GrpH0b4CvUWCyqmarjr43WyApNTjZZXMBpz1bPzFO2WOJxMy+CHwxfX0fST+z4wzHwGKf4Og4BbhIHGc0Le/TnCuRbLW4pBXxCLhIHGck4w0mNpK5E8lWsSatsSL4TVeOU4z3bjWPebcmbbIiAANnUbWIuRSJ0yCmM+VkpsytSObVmrTNioA8cHecQtySNJd5sybtsyIpLa12Ru0iifr6A0NVa4rnRSh1CKT/t5kKPk4yOdkPU7VYnOFUIo4c3rs1Jaq0MG23JtO2IlWLYgMtF0l7W86YtNWfb2u954nGWJJ+qrAsbbMo0xJI7ZajD3e3WkZbhDI3FsTwaSl1kb8aTmpVmi6UaQhk1tZjHQ2qymZojUjyRGZzK5S5Ewjtd7ea10pKEplN3Bia5tJMWp9pfCeVUENKoTESZv9bSXelCbP/Q5qNdCStFUnGpA0jTv9myaR1aKw4MhqSMFvST5PkkfsJ4MeBnwJeUVRw60WSMQ2xzIK5FgeJq1V2mZDChNkkUlwCFoBFoAs8VFRwK2OSUUwygp812DpilUnF0SrK927tkXRrbv1gmia3DIUJs83sy5JuBh4gSeHy+2Z2uH+/fuZOJBlNFcuWEkfKGFai0oTZki4EfowkLS/A5yS9zMz+YlSl5lYkGdPsOp4FbRXGOpqTMPv1wFfM7On0mP8FvAQYKZK5iUnK0LYG17b6DqS+mKRMwuzvA6+Q1JHUJQnaC92tLSUSWAt2mxj0NrluE1FP71ZhwmySh079NfAN4OvA183sfxYVPPfuVhFNmKo/V4IYgGroOCyTMDvNKv9Pxi27VpEYVklX6zQC7DpvBqtCFLMe65ln5sKSDGsgTZx2MimtFEPLDeVciGQYZRtUE8TUysZfhukE5TNlrkVSlnwDzQQTmRFU4HHlXa25FUY/LpL5oo6BxC0jjgwXSTOJNvnLhDU+cKYNdZwUUU/vVpXMhUg229hGlxUTo6lblJh44voOOr6xwtkqMYmkXcBVJNOLDfiHwLdIHhZ/PvBd4E1m9ngVlexnM40sHrPbtcou4HHrEpSoS/930ijRtFwkZS+THwY+Y2bPBZ5HMpRfOH9/mkTY6lKW2Gx1KVPuZs6xGcY9Z5nPMeocM6eeEffKKLQkks4AXg68DcDMloFlSa8DXpnudjXJU3nfPc3KTdNiDCprlKscYXTHPnsxETb0vDE28KqVtwr5z1fGwmTnHFRWXbTd3SpjSS4Afgj8V0lfk3SVpB2UmL8PIOlKSbdKuvXRR8tFcJu9ApYVSMyaQPrnS2XLuC7ROHUcds7+ug2rf76scZmJZWm5JSkjkg7wk8AfmNkLgGP0uVZmNvRjmtlBMztgZgd27x58uklcnWGuSH95WePr5RrlshkR0EuXk7a2VNWYeti682TnjiCpT7r0zNaLecj3M6krVrlLZknvVpmlqZQRyRHgiJkdStevIxHNQ+m8fUbM3x/JpD/QOJZj4z5r70WWLFmjHHQ1nyb5c+TPna9X//55pmlZ8mVWJpaWW5LCmMTMHpT0A0k/ambfIplpeXe6XE4yJXnY/P0NTOuHGNQghokjc2X6hREDESIyESMiRJxmQO+qxxlTqel6TppxNE6inUBGiBGQ/A/TmCRWcvWKzAizA9P4I7uqZZ+1P8aIzUrHKoOoQihtj0nKjpP8OvAJSQvAfcA/IPm9rpV0BfA94E3VVHEjZa6YG6++a9sHCaRnQbKOiC2gZ70p1zqhZ3DSOgSKCc2IMbqKkyupAGz1dZDWO8wdH1Ns/icVytTZCiIxszuAQfceb5i/XzWbtSDZtl6667IFxIiTFhIhjseLLFvIcVvkZNylq4dWb4SeJo/Gi3y3t4eloMd2nWJBEduDU4QYS0REGAuKiQ26+Xaeu+8lL5SqLMrUaLgrVYZWjbhPYkEgsSDJNtEjIDbRs5BlQo7ZAifjBY7GSxyNtvHMzpPAqWlVfZXj8SI/XDmd08ITREHAUrBMYDELRIm7pTXXKzIItfY52mhRxNZxt2ZOmSA9L5Ao7bnKtmcuVmZBMvfqmC1w0ro82NvF0XiJrz5xAd87eiZPnreNSxbvnOpM4MiM//P0c/nsA8/lgtMf48AZ32VXeJywY/S0suqChWkQvZBzwzLXi1yK136LMmgMpBFCcZHMjmFB5vrp6OtjkCxAX7Ug8SInrcsjK6fxSG8ndz50Hie+fxpf3XYczrpzqvWNiTn02Pk8+FfP4KlnLXHu0pOc6nbZHpxiST26ilggSmKUtL79MUr2+QZNmxkmlJnjIqmezbhZg4hsrQdrLUDv0LMOJ+MuPQsJg5h4Keab3zuPFz9+OUvdFfZsP0YniOkoIlDiEhXXOTnPioWsxAGPHN/ByV6Hp5/ahi3FdIKYnoXp0iHEiC0gUpz0sAkCE+EULsMztyY1iETSG4HfJcmrdUl6b/ug/S4jmWYVAleZ2fuKym68SMaZZgIbg/V+KxKbWCYRx7KFnLQuJ+Mup+IOy3GH7YvLHDutx46vbmPf9Q8Tn7mToxfuJ+qKaFGYwALILtj5Z2autufcAFq4bITLxu57niQ4+jRHXn8mx37qBNsXl1mOE3GejLuEQcyyJVHHsiIWLCaWNlgTSH7dzJr0xyaNc7usNnfrm8AvAP952A6SQuAjwGtIxv9ukXSDmd09quDGi6QMo4L10celdyESJIuJQIZkxAtgp+8gXuqi2AgiYDkVhUY/UFZpj44MgsiQGfH2LrCDeAGy5IKRZVYtOX+8FnmMZDNB/EypQSRZutKCJPGXAPea2X3pvteQ5BCef5FkDMtCko2H5MdBltOu38zlSdyjAAFhJ+L4uTEPvGI3wYoR9NZPnVCB+2cCS3+slUVhARx7xk7ijjh+bky3ExHKUncvXF0ixLKFLFlAJCNAI12uaTynpQ7GmHIySS7gMuwFfpBbPwK8qOiguRLJuEQWEFlyDc5iiAzrGis7IFgW3dhWXR0ZGBopFJNWrY2FybKyTcSLSbkZcZ85iiygq42WJEIELY5+68gFbGalZnxshrkWSTTgx1m1JjkHJbb1zoqUimQbhCEEPRFERmBAnApFg4WSWZAsdom7EIci2gbRIljHyF/8Mwu2KlaCVBTaMFU/MipJTlEpUxxMHJULuCT3A/tz6/vSbSNptCtbB1EuLon7A95gLQZZpX99EAOOWQv4893TbWvxm6Q5ExxvAS6SdEE6xerNJDmER7LlRRL2hf0GxLFQT3SOQ+c4hKeSuCSIQBEoToLxQT+0zJL3o2T/oJf0cHVOQPcYaDkgjrXB1ZpXshH3qhNmS3q9pCMkWeL/XNJN6fbzJN0IYGYrwNuBm0jurr3WzO4qKnuu3a3NEJswA62I4FTSwBNhsCoCYPiVL+2qlVkSu0TJtmA5PT5Kyh92eDCH6YYUV28mzOzTwKcHbP8b4LW59RuBG8cpe65FEirx45NpHsmVO8TS22RjQolQcTIdJG2cURwQrYRseyjgnK8eI14MOXVWN4kvOtrofg0itSqBGcFKIrLF7/QITkWsbN/ByT0hUbxmxAPFhOmS1TFbsvXs87SOOZjg2Hh3a5wBsKLu0DD3a4WK0/s31gQSYESxsEgsPWLo0DdZ/Kv7WXhyhc6JOOkOzlmVgaQDiYmrZXSOxyw8tcLC4SMEt9zF0qOGRSKKA4LVCY1xWg9bFe00Pm+eWY641/R8kspohSUJpA0j7yFaHXUPWD+gGLI2JBcKYiMdl1BiVdKJhBExoYyuIjpBTCDj5HIXnu7y5HPg+HtfxMp2o3fWCoSGur3ElQrW6pJve/kqWqxkoLwXQBTQffWzCU9cSO80g6e7nDy9Q5Ceu6uIBUWrQkk+Q3ozljZakfxAYvb589/LoO9vpjRYAGVohUjKkAkllNbd0bc2Gp3MrA1khGZ0tQLAknpEUtJYg4hTp7p0ngx45gsf5LqLP86i1hvbcSYQ9k+dOWUxv3DXW3n49nM4dVaXbpCcs6uVVaF0tcIC8er8sOxuxewzZp8r7LtTsck02UqUodUiyVuTUvsrTZqQ3a+hmK5WWArEGZ3jxCbOPvMoD54KefYZj7BTXbrqv25PUt8e55/+GPefeybnnXmUXZ3j7AxPshT06GZT5XOWZNzJjY20IuCWZFzGbdgZ2Y89yu1a3ZZaE0ivtOkcwcyadBUTpSMjpwcnWbIecfcxzu4c5Yx9J3jsmTt4yc57pz7lIyDg9XtuZ++2Jzirc4xzuk+ypB67guMsKGJJPcK0ftl97wGJq7XOpRpQr2kLZGpT7q3ZmVDKMBNLkv0A0046MCo2Wd2WzCtZvWIHilkAloIegcWc2TnGYtDjtODEVOuWnFvsCo9x7sITbA+WV+8jydystWQQg63IqFhkanWc8gCn35k4If0/SBnRlA3i+2MTlDYyW2uAMSLkVBqTrNCzDruC4/QIOTs8xsZmOTlnh8e4aPFBuqzFIEupQDILspBeerN73ItikUmsSC03abX8mZCNiknKumJFQlndNiSIT0SUWJQg/d9Vdo95zIJFLKl/Qvp06BKzQ8upBYtW45BA6y3IqGC9/3P30yiB4JZk6pR1xUYJJWtgoyxKkFqUbKCxS0wksWQrIFiq6JfdEcTsshPpZ1jvXuVjEKjOgtR6i+8cDCY2TiQZZazKpl0vWBfMh0qEEqS3zSbPJamGAFbdKmBVHNl7+SC99QJJ8cC9QspYlXGEAqzzjwOtBfqBGV2tZXxfqKjrdEFKXbmcGHKWI6l/tn2tDpMKZJYJIlwkNVBkVQZ1Dw8bkR9oVWDd/ePZMVWQWJKN2/LnhvICaar1WMXwwL0uNmNV8sf0xymrZJkRWd9dXFXDCvvuMhwmDGi39cjjgXvNFIllmFXJjsk3vH7BZFlIqmZUb1W/BWuzOFZxkcyGMi5YUTdxv2AgFyRXFJME0kaXaci+8yAQH0ycMfkGMUgw+QYVr7pVGxtRv4WpmqJHvm3Yv2XCWEd6p2ababVI8mwmuM8fW2bbNAjTJHllaHxQXpZ2a6RekYji56LHE9y+Wmaay7CGV9UzEssw1o1lUxTFtJ9RP4y2u1uNux0hSP+mQchG/3/oeaV1S5Vs5lzjfJbC80/xOy7ESO56K7NMgKQ3SrpLUixpYO4uSfsl3Szp7nTfd5Qpu7Hu1rAfcTOWZlgcspnzT0pAQFSQynSaYmgE9ViSwlzAwArwLjO7XdJpwG2SPjd3uYCzH34StwxKumY1NbJpxxWNEUdKHe5WmVzA6SPVH0hfH5V0mCT16UiRlPo2Jf2z1Dx9U9KfSFpKE3wdknSvpE+myb5qIxjwNwmZOzNNt6aOc037e6gCxVZqIc0FnFuurKxO0vnAC4BDBbsWWxJJe4F/ClxsZickXUuS+e61wAfN7BpJ/wm4AviDSSo+KdOyMm2giWIYyHizgGvJBSxpJ/Ap4J1m9lTR/mXdrQ6wTVIP2E5isi4Ffjl9/2rgd5mxSDJGNaBZCqhMPNK/f9tJBhOn429NIRcwkrokAvmEmV1f5pgyz3G/X9K/A74PnAA+C9wGPJGmjYQkhf3eIZW6ErgSYP/e6d/ENC5NaHhNqEOtNMSwKwlYPgocNrMPlD2u8NeSdCbJg04uAM4DdgCXlT2BmR00swNmdmDP7tmLpCxteO5HW5BZqWWic5TIBQy8FHgrcKmkO9LltUOKXKWMu/Vq4Dtm9sP0pNenJ9slqZNak1Ip7J31mVy2BDXdmVgmF7CZ/SXFSWo3UMbufx94saTtqbl6FUmX2c3AG9J9LgdKB07OVqJcz1aT53cVisTMDgHXAbcD30iPOQi8G/gNSfcCu0l8PcfZiFm5paGU6t0ys98Bfqdv830kD2p0nOF4cjrHKUGDrUQZ6k9z6r1GW+87aLdG3JI41aO43f6Wi8SpFqMxg4mbxUXiVIqYfKBw1rhInOpxkThOAS4SxxmBxySOU4z3bjnOSJo95aQMLhKnWjxhtuOUoN3elovEqR4fJ3GcIlwkjjMCM4ja7W+5SJzqabkl2WJpO5yZUMOdiWVyAef2DSV9TdKflSnbReJUS00Js1nLBfylEvu+AzhctmAXiVMxBhaXWyY5i9lhM/tW0X6S9gE/A1xVtmyPSZxqMcYJ3PdIujW3ftDMDk65Rh8Cfgs4rewBLhKnesrHG5XmApb0s8DDZnabpFeWrZSLxKme5uQCfinwc2nWxiXgdEkfN7NfGXWQxyRDqPIhPluLkj1bNXQTm9lvm9k+Mzuf5MkIXygSCLhInKoxII7LLRNQMhfwpnB3y6meeqxEYS7gvu1fBL5YpmwXiVMxPi1lLqg7Thh0vrl9OpeBTTgGMmu2vEiaEkgHBPMrlAZnjC/DlhVJU8SRZ26f+djyCY5bViROTZhN3HM1a7akSJpoRfLMnevllqRdNF0gGfMjFMOi8k8cbiJbTiROzWRT5VuMi8SpHu8Cbg9tcbUy5sHlMsDckjSXtoliEP2foXWiMXNL0lTmQSCDaKN1aXvgLquxe07SD4FjwCO1nXQy9tCeukI19X2WmZ292YMlfYakXmV4xMwu2+y5qqJWkQBIunXU3WdNok11hfbVty3Mp0/iOFPEReI4BcxCJNPOflElbaortK++raD2mMRx2oa7W45TgIvEcQqoTSSSLpP0LUn3SnpPXecti6T9km6WdHeaePkd6fazJH1O0j3p/zNnXdeM/sTPki6QdCj9jj8paWHWdZwHahGJpBD4CPD3gIuBX5J0cR3nHoMV4F1mdjHwYuDX0jq+B/i8mV0EfD5dbwr9iZ/fD3zQzC4EHgeumEmt5oy6LMklwL1mdp+ZLQPXAK+r6dylMLMHzOz29PVRksa3l6SeV6e7XQ38/Ewq2Ed/4mdJAi4Frkt3aUxd205dItkL/CC3fiTd1kgknQ+8ADgEnGNmD6RvPQicM6t69fEhksTP2USu3cATZraSrjf6O24THrj3IWkn8CngnWb2VP49S/rLZ95nnk/8POu6bAXqmgV8P7A/t74v3dYoJHVJBPIJM7s+3fyQpHPN7AFJ5wIPz66Gq2xI/Ax8GNglqZNak0Z+x22kLktyC3BR2vuyQJKs+Iaazl2K1Kf/KHDYzD6Qe+sG4PL09eVAYYr/qhmS+PktwM3AG9LdGlHXeaAWkaRXtrcDN5EExNea2V11nHsMXgq8FbhU0h3p8lrgfcBrJN0DvDpdbyrvBn5D0r0kMcpHZ1yfucCnpThOAR64O04BLhLHKcBF4jgFuEgcpwAXieMU4CJxnAJcJI5TwP8H0irjYbujh2IAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(df_dv.reshape(nz, nx))\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "markdown", + "id": "2b886600", + "metadata": {}, + "source": [ + "### connected structure" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "ee852785", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMEAAAD8CAYAAADOpsDvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAR70lEQVR4nO3df5BdZX3H8fcnm4QQEJIQSzGJkhnT1ozWwckAFlsRcAzoBDtlHKLFaKlpO8ZiYdTYdtChdir9IdpOqt0pVupQY4qMZDQaFWEcrWCCUGoSkW0skgiG8MNfiMnufvrHOZHLmt17Nnvv7t19Pq+ZM3vPued+z3M3+83z69z7yDYRJZs11QWImGpJgihekiCKlySI4iUJonhJgihekiCmDUkflXRA0rdGeV6S/lHSgKR7Jb2kSdwJJYGk1ZLuqy+6cSKxIhr4GLB6jOcvBFbU23rgw02CHnMSSOoDNtUXXgmslbTyWONFtGP7K8BjY5xyMfDvrtwBLJB0Wru4sydQpjOBAdt7ASRtrguxe7QXzNVxnscJE7hkTLYf8/hB288+1te/6hUn+NHHhhqde9e9P98FPNVyqN92/zgutwR4sGV/X33sobFeNJEkONoFzxp5kqT1VFUT85jPWTp/ApeMyfYl3/TARF7/6GNDfGP7cxud23fa/U/ZXjWR6x2LiSRBI3Um9wOcpEW5UakwBoYZnqzL7QeWtewvrY+NaSId42O6YJTFmMMearR1wFbgjfUo0dnAD22P2RSCidUEO4AVkpZT/fFfCrx+AvFihupUTSDpE8C5wGJJ+4D3AHMAbH8E2AZcBAwATwJvbhL3mJPA9qCkDcB2oA/4qO1dxxovZiZjhjp0u77ttW2eN/DW8cadUJ/A9jaq7IsY1TC93RXsesc4ymZgKEkQpUtNEEUzcLjHP8KbJIiuMk5zKApnGOrtHEgSRHdVM8a9LUkQXSaG0FQXYkxJguiqqmOcJIiCVfMESYIo3HBqgihZaoIonhFDPf59DkmC6Lo0h6JoRhxy31QXY0xJguiqarIszaEoXDrGUTRbDDk1QRRuODVBlKzqGPf2n1lvly6mvXSMI4ChzBNEyTJjHAEMZ3QoSlbdQJckiIIZcTi3TUTJbDJZFqVTJsuibCY1QUQ6xlE2o3yoJspWfeVKb/+Z9XbpYgbIl29F4UxmjCN6vibo7RSNac8Ww57VaGtC0mpJ90kakLTxKM8/V9Jtku6WdK+ki9rFbHtlScvqoLsl7ZJ0RX18kaQvSrq//rmw0buIolQd475GWzuS+oBNwIXASmCtpJUjTvtLYIvtM6hWVP3ndnGbpN8gcJXtlcDZwFvrC28EbrW9Ari13o8YofqMcZOtgTOBAdt7bR8CNgMXjzjHwEn145OB77cL2rZPUC+G/FD9+MeS9gBL6oufW592A3A78K528aIsVce4cZ9gsaSdLfv9tvtb9pcAD7bs7wPOGhHjvcAXJL0NOAG4oN1Fx9UxlnQ6cAZwJ3Bqy2rhDwOnjvKa9cB6gHnMH8/lYoYYx4zxQdurJni5tcDHbP+DpJcCH5f0QtujrhXSOAkknQh8Cni77R9JT2e3bUs66qI8dSb3A5ykRT2+cE90WodnjPcDy1r2l9bHWl0OrAaw/XVJ84DFwIHRgjZKUUlzqBLgRts314d/IOm0+vnTxrpIlG2YWY22BnYAKyQtlzSXquO7dcQ53wPOB5D0AmAe8MhYQZuMDgm4Hthj+wMtT20F1tWP1wG3NHgTURgbDg/ParS1j+VBYAOwHdhDNQq0S9I1ktbUp10FvEXSfwOfAN5kj72GbJPm0DnAZcD/SLqnPvbnwPuBLZIuBx4AXtcgVhSmag51bjrK9jZg24hjV7c83k31N9tYk9Ghr8KoU37nj+diUaZenzHObRPRVeMcIp0SSYLoss42h7ohSRBdl88YR9Gq0aF85UoULB+vjCDNoShcRociyMcro3C2GEwSROnSHIqipU8QQZIgCpd5gulsVh/q64NZHfwHHDYeGoLhoc7FrGn2bNCsrsWfiMwTTFNPXryKfWsGUV/nPhHqIbHs030cf8s3OhYToG/hQh744xfw5IpDPPfmWcz7TGfjT4QNgw0+MDOVkgSjOPiiPna98p+YP2tux2L+ZPgpzhy4kmUd/gyeTpjP8ld9lxuffzO/df9VLP1MZ+NPVJpDUbT0CSKoJsx6WZIgui4d4yianT5BFE8MZXQoSpc+QRQt9w5FuOoX9LIkQXRdRoeiaE7HOCLNoYiMDkXZ7CRBRIZII9IniKIZMZzRoShdj1cEzdfWjDgmdce4ydaEpNWS7pM0IOmoC8hLep2k3ZJ2SfqPdjHHs4RrH7AT2G/7NZKWU60ofgpwF3BZvcp4xDN1qCqo/wY3Aa+kWsh7h6St9TplR85ZAbwbOMf245J+pV3c8TSHrqBaMfCkev9a4DrbmyV9hGr92A+PI16RBk80s5cu6WjMoV9dyII53+9ozE7q4BDpmcCA7b0AkjYDFwO7W855C7DJ9uPVtd12aeFGSSBpKfBq4K+BK+tlXc8DXl+fcgPwXpIEYzpec3nbmm1sf+nKjsZdMOf7vOM5n+9ozE4xMDzcOAkWS9rZst9fLwZ/xBLgwZb9fcBZI2L8GoCkrwF9wHttj/nLaVoTfBB4J/Csev8U4Il6XdkjhTnqf2+S1gPrAeYxv+HlZqY+zWL9ggHeePK3OxsXcbzm8hP31vcNAVUWNK8JDtpeNcErzgZWAOdSrXj/FUkvsv3EWC8Yk6TXAAds3yXp3PGWqM7kfoCTtKjXBwq66uc+zMvufgNPfm1xR+MOnmD+5Hc/x5tP3tPRuJ3SwXmC/cCylv2l9bFW+4A7bR8GvivpO1RJsWO0oE0X814j6SJgHlWf4EPAAkmz69rgaIWJEQ57iJ9+fTHL/ua/Ohp39tIlfOllL+jZJOjgGOkOYEU9KLMfuJSnm+RHfBpYC/ybpMVUzaO9YwVtO0Rq+922l9o+vb7ol22/AbgNuKQ+bR3Q4a+Uipmh2fBok85z/R/uBmA71SDNFtu7JF0jaU192nbgUUm7qf5G32H70bHiTmSy7F3AZknvA+4Grp9ArJjJOtgItr0N2Dbi2NUtjw1cWW+NjCsJbN8O3F4/3ks1ZBUxOoObjw5Nidw2EZMgSRCl6/ExwSRBdF+SIIo2vsmyKZEkiK7Lh2oiMjoUpVNqgiiaScc4Sqd0jCNSE0QMT3UBxpYkiO7KPEFERocier5PkO8diuKlJoiuS3MoymZy20REr/cJkgTRdWkORSQJonhJgiiZnOZQREaHIlITRCQJomjpE0SQmiBCPf6hmtxFGsVLTRDdl+ZQFC0d4whSE0QkCaJoIqNDUTo/fRNdu60JSasl3SdpQNLGMc77PUmW1HZd5EZJIGmBpJskfVvSHkkvlbRI0hcl3V//XNjsbURx3HBrQ1IfsAm4EFgJrJW08ijnPQu4ArizSfGa1gQfAj5v+zeAF1Mtn7kRuNX2CuDWej/il3UoCagWihywvdf2IWAzcPFRzvsr4FrgqSZB2yaBpJOB36FeotX2IdtP1Be/oT7tBuC1TS4Y5RlHc2ixpJ0t2/oRoZYAD7bs76uPPX0t6SXAMtufbVq+Jh3j5cAjVCuEvxi4i6qqOdX2Q/U5DwOnHu3F9RtZDzCP+U3LFTNJ89Ghg7bbtuFHI2kW8AHgTeN5XZPm0GzgJcCHbZ8B/JQRTZ96AeWjvlXb/bZX2V41h+PGU7aYCVyNDjXZGtgPLGvZX1ofO+JZwAuB2yX9H3A2sLVd57hJEuwD9tk+0sm4iSopfiDpNID654EGsaJEnesT7ABWSFouaS5wKbD1F5exf2h7se3TbZ8O3AGssb1zrKBtk8D2w8CDkn69PnQ+sLu++Lr62DrglkZvI4rTqSFS24PABmA71eDMFtu7JF0jac2xlq/pZNnbgBvr7NsLvJkqgbZIuhx4AHjdsRYiZrgOzhjb3gZsG3Hs6lHOPbdJzEZJYPse4GjtqvObvD4KljXLonQid5FGJAki0hyaphbeN8yqO/6Avr7O3QI5ONjHou8MdSzetJEkmJ5OuvmbLPjcvI7HHf7ZU73+N9FZ+WTZ9OXDhxg6fGiqizEzJAmidL3+oZokQXRdmkNRtkyWRZAkiLJlxjgC0HBvZ0GSILorfYKINIemreHfPoP9Lz8ed/CbmTQMS25/kllfvadzQaeDJMH0tP/lx/OVP/o75mtOx2L+xId5xdA7WPrVjoWcFlITTFOeBfM1h/mz5nYs5vDwcEdrlmkjSRBFc26biMJlniACwL2dBUmC6LrUBFG2TJZFpGMckSSIwpl0jCPSMY5IEkTJMlkWYedDNRFpDkXx0hyKshlIcyiK19s50Hgx74hj1qk1ywAkrZZ0n6QBSb+0gLykKyXtlnSvpFslPa9dzEZJIOnPJO2S9C1Jn5A0r15B8M66MJ+s1zOL+CUadqOtbRypD9gEXAisBNZKWjnitLuBVbZ/k2ql1b9tF7fJivZLgD+tA78Q6KNaOvNa4DrbzwceBy5v+y6iPE2Xb21WE5wJDNjea/sQsBm4+BmXs2+z/WS9ewfVWsdjatocmg0cL2k2MB94CDiPKtMAbgBe2zBWFKSaLHOjDVgsaWfLtn5EuCXAgy37++pjo7kc+Fy7MrbtGNveL+nvge8BPwO+ANwFPFGvKztmYeo3sh5gHvPbXS5mouZ3kR60Pebq801J+n2qFVdf3u7cJs2hhVRVznLgOcAJwOqmhbHdb3uV7VVzOK7py2IGGUdN0M5+YFnL/tL62DOvJ10A/AXVavY/bxe0SXPoAuC7th+xfRi4GTgHWFA3j0YtTESH+wQ7gBX1oMxcqr7p1tYTJJ0B/AtVAhxoErRJEnwPOFvSfEmiWsB7N3AbcEl9zjrglkZvIwrTbGSoyehQ3fzeAGwH9gBbbO+SdI2kNfVpfwecCPynpHskbR0l3C806RPcKekm4JvAINUQVD/wWWCzpPfVx65v+y6iTB38UI3tbcC2Eceubnl8wXhjNpoxtv0e4D0jDu+lGrKKGF2+fCuCfLwyotfvHUoSRNdpuLfbQ0mC6C4znsmyKZEkiK4SjSfCpkySILovSRDFSxJE0dIniMjoUBTPaQ5F4fKFvBGkTxCReYKIJEEUzYah3m4PJQmi+1ITRPGSBFG0fCFvhMHpE0TJTDrGEekTRCQJomy5gS5KZyC3UkfxUhNE2XLbRJTO4MwTRPEyYxzFS58gimZndCgiNUEUznhoaKoLMaYkQXRXbqWOoOdvpW66mHfEMTHgYTfampC0WtJ9kgYkbTzK88dJ+mT9/J2STm8XM0kQ3eX6QzVNtjYk9QGbgAuBlcBaSStHnHY58Ljt5wPXAde2i5skiK7z0FCjrYEzgQHbe20fAjZTLTTf6mLghvrxTcD59dLDo5IncfhK0iPAT4GDk3bRiVnM9CkrdKe8z7P97GN9saTPU5WriXnAUy37/bb7W2JdAqy2/Yf1/mXAWbY3tJzzrfqcffX+/9bnjPp7mdSOse1nS9ppe9VkXvdYTaeyQm+W1/bqqS5DO2kOxXSyH1jWsr+0PnbUcyTNBk4GHh0raJIgppMdwApJyyXNBS4Fto44Zyuwrn58CfBlt2nzT8U8QX/7U3rGdCorTL/yjovtQUkbgO1AH/BR27skXQPstL0VuB74uKQB4DGqRBnTpHaMI3pRmkNRvCRBFG/SkqDddPdUk7RM0m2SdkvaJemK+vgiSV+UdH/9c+FUl/UISX2S7pb0mXp/eX2rwEB968DcqS7jdDApSdBwunuqDQJX2V4JnA28tS7jRuBW2yuAW+v9XnEFsKdl/1rguvqWgcepbiGINiarJmgy3T2lbD9k+5v14x9T/XEt4ZnT8DcAr52SAo4gaSnwauBf630B51HdKgA9VNZeN1lJsAR4sGV/X32sJ9V3Hp4B3Amcavuh+qmHgVOnqlwjfBB4J08vi3cK8ITtwXq/p3/HvSQd4xEknQh8Cni77R+1PldPukz5mLKk1wAHbN811WWZCSZrsqzJdPeUkzSHKgFutH1zffgHkk6z/ZCk04ADU1fCXzgHWCPpIqqbzk4CPgQskDS7rg168nfciyarJmgy3T2l6jb19cAe2x9oeap1Gn4dcMtkl20k2++2vdT26VS/yy/bfgNwG9WtAtAjZZ0OJiUJ6v+Zjkx37wG22N41Gdceh3OAy4DzJN1TbxcB7wdeKel+4IJ6v1e9C7iyvmXgFKqkjjZy20QULx3jKF6SIIqXJIjiJQmieEmCKF6SIIqXJIji/T99ILvGL5Ov6gAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "v[30:,:,40]=1\n", + "v[40,:,30:40]=1\n", + "plt.imshow(v[:,0,:])\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "6fd7737b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "-0.967774527264677" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "T, f, df_dv = mpa.constraint_connectivity(v, nx, ny, nz, p=6, cond_s=1e6)\n", + "f" + ] + }, + { + "cell_type": "markdown", + "id": "c02f16f7", + "metadata": {}, + "source": [ + "If we connect the structure so that all solid pixel is connected to the bottom, we can see that the constraint value become negative, and the tempeature is low overall as heat is leaked through the bottom side. The gradient values are relatively small in generall as well, indicating no more change is needed to connect the structure." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "50dd1f7f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM0AAAD7CAYAAADJnxDZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA85UlEQVR4nO2de6ws2VXef2tXn3PvzODx+AHWeOzEluwkMiSIeASWHEVgg3HAYSyFEINlnMRh/uARE5CwicRDkZFAigATEMqAIQaBjGVDPOJlwNiKiMLAGFsxtnE8GAwexu8Zex73ntNde+WP/ahdu6uqq7urq7vPqe+qbnfVqa7aXV1ffWutvfbaoqpMmDChP8y+GzBhwrFhIs2ECWtiIs2ECWtiIs2ECWtiIs2ECWtiIs2ECWtiK9KIyItF5EMicp+IvHaoRk2YcMiQTftpRKQA/h/wNcDHgD8FvllVPzBc8yZMODzMtvjslwP3qepHAETkTcAdQCtpnvhEo7c9rdjilPuFImvsO8T5dos+3+cv3nf+aVX9wk3P8bVfdZN+5rNlr33f/X/P3q6qL970XGNhG9LcBvxtsv4x4CvynUTkTuBOgKfeZvifv/XkLU65PewWny3XIE2pq/e1K47X93y2x7k2Pf6XP+NvPrrRwT0+89mSP3n73+u1b3Hrh/d7c/TEzgMBqnqXqt6uqrc/8YlT3OGyQQHb89+xYBuluR94erL+NL9twoQIRZlrP/PsWLANaf4UeLaIPBNHlpcB3zJIqyb0xqam2Zg4JhXpg41Jo6oLEflO4O1AAfyCqr5/sJZN2CnW8c+2gaKUFyyTfhulQVV/G/jtgdoyYSSMRZgAu/M44LjYijTHCMPmEbQCHf2GO3YoUE6kmTBhPUxKcwGwjdocM/ahkgrMJ5/mcmMy0daDopN5dlFwEdRmnXDz3oiuUF4szlxe0kwYBy4j4GLhUpNml2pTiPbKP9sUR6EyAMiFM2cvNWk2xTH5NftupwsEHMe16otLT5qL4Nu0Yd+EgdBPs/92DIlLTxrYjDh91GYbE+0i3WjHkB+3DibSXFAcCukmpbnACCN9jsFUW/XkHuImHepGV4TygtVvmUizBY4pILAKu/wek3l2wbGufzM2cYZUmTHarQjnerx1IZowkaYBFzmiFjDeeBqwk3l2OTAUcTaJoG16Q/f53D7MyYtiwgZcrEfAnlAcQELioRJGVSjV9Fr6YFWBShG5IiK/5v9+j4g8w2//GhF5t4i8z7++IPnMu/wx3+uXL+pqw6Q0HRgqotamNgZdWcYpxS5LNe0S63zHLvgClT9DUqBSRO7OClS+CnhQVZ8lIi8Dfgz4N8CngX+pqn8nIl+CG6Z/W/K5l6vqvX3aMZFmIFykSNqQcIGAwW6zPgUq7wB+2L9/C/DTIiKq+p5kn/cDN4jIFVU9W7cRk3nWA4Z+F2ofZlrZkRDZ9bdVsGqwPU2mLoRAQJ+lB5oKVN7Wto+qLoDPAU/K9vlXwJ9lhPlFb5r9gIh0XrSjU5qi4etctPEaTdh1X8cQBGnDGoGQJ4tIaiLdpap3DdkWEflinMn2omTzy1X1fhF5HPBW4BXAL7Ud4+BI00SKIT4zBLG28XF2PVQgRx+F2SVRAtbMCPi0qt7e8fc+BSrDPh8TkRnweOAzACLyNOA3gG9V1b+MbVS9378+LCK/ijMDW0lzMOZZIZsRZt3j7/IcMK6J1mWWdWEo06svwvlWLT0QC1SKyCmuQOXd2T53A6/0778R+ENVVRG5Bfgt4LWq+r/DziIyE5En+/cnwEuAP+9qxF6UZtc37jrn30SBVinOGEGBTQgzJlECXMLmMOdtK1ApIv8FuFdV7wbeAPyyiNwHfBZHLIDvBJ4F/KCI/KDf9iLgUeDtnjAF8AfAz3W1Y1TSCPsnTI5tCbTeuZZNtD5h5yH8mX0QBpx5Nh8wjaapQKWq/mDy/jrwrxs+9zrgdS2Hfe46bTg4n2afCATqS55DS7dpU5l1CTOkSqrSu+PyWDCRpgHrkKeNOMG32YWZ1nTMpm19yLL7viUZrHPzUDCRpgN9yXNoitMX42Q5T0pzKdGHPH2Js23oOb/R11WY/SRsTqS5tChkfeIcUnrNfsrSyjQI7bJj3WABLBNnHbVZ94ZrUplVZNmlErgSThfrNrtY32ZEtKnOmP5NToZ1CTOO2TQVC7wQ6LpV1hrq3KI6uyJOevNtQ5gxfQxlf31Eu8KlIU3fny3db8gbf2zfpjksvZ+bd1KaI8K2t0j++da0mQbFaUq1SYmT+jVtWQGpP7NOesxyhK37SuxSCVTlwinNym8jIk8XkXeKyAdE5P0i8mq//Yki8vsi8mH/+oTdN7cf+o5/Gfq4Y6UIpaTYlDBrJkpuDBcIKHotx4I+V2wBfK+qPgd4HvAdIvIc4LXAO1T12cA7/PqlQBd5cuKM+YztQ5ixM5xh2BoBh4CV5pmqPgA84N8/LCIfxI2OuwP4Sr/bG4F3Aa/pPBbOhNnVE3nsy96W7ZxH1tLAQJuJ1hdtKlMPEjSTpc8xh4YLBFxin8ZX9vgy4B7gKZ5QAB8HntLymTuBOwGeepv74cINdXAZz9l62fNzQ0TLurKdN80rW7XvWA76pc0IEJEvwA0F/W5V/Xw6jNoP8mns7vPDVe8C+Mf/5KS2z1Cqs8lP0seCbtuniUy56uTBgTa12RZdCrPK3+k87kDm0qXNCPADdN4K/Iqq/rrf/AkRuVVVHxCRW4FPbtKAXZprTRjC3QzH6KtEAY1pNmuYaOvc9GsPB9ihT3HRKmz2iZ4JbjTcB1X1x5M/pcNKXwm8bdNGjFEYo2AYwjQdMz1uHiRoG2K9zbDoQIgmlckd/VCRprbviE64Ksyt6bUcC/oozfNx1TneJyLv9dv+M/CjwJtF5FXAR4Fv2kkLV+BQLnVBXXl2lRXQ3PPf32cZO0rlzLND+ZWGQZ/o2R9Bq13wwqEasqvgQC/fpbvMVStKratFbrbV/BgfUduUTG1DAlKFad235aYdy2yaMgIuGDYlTPrZJvIMRZx1b7g+hBnTx7j0IecxsE5gYNVP36UyfcjSdvxlZz71F7R27pLuIdEl0jsYkPsyucrkfkvV3n0OSruE5tk+sOuIWhdh+vy8XUmdhciS8oTPWJbVZtPw8zaEGdtcmmoEHAnaVKaNMJs+C5sIlJptTYqzagRojq4hAem2LrKsHIi2IzVw0bPjySvrg4MlzTZqsw5hOhMwe835Ut39yx2czaqT7l/tq6DtT+W6o2+WTTW/npKldTzNiObSpe3c3BeGNNPWIUwfsjTtGwiUkieeNyHPJmrThT6E6UuUXQQJJvPsgiC/NbqIYjp8IJuQISdQXUkq1Um3u3Y0+zVNYeYyI0ZOmDbfJt1nTEzRsyPFJmHlLqK07ZcTKBAHEtVRXTv9JodVU892ptlUy/++LoYy46bo2cgYOpK2SmGayNLXt1n6rNZXUuIAkJhoIYpmMtUpfcg2bIu+TKIwXX5Nun0fUBUWE2m2g2X9SNVYSZ35Tb+Jb1MLDIjU1Kc5YRNsT9+mbDLBugIBB3KzTubZAMhvnE1/2j4E7DLNUlKkhKltX7t12bcLh4rEUEjUJnyHdTMEokmWKk5HKs2+cBF9moN4FPW5YYacyQx2RRj3mfRzBUKBYETieQzNZC6ahyQtoczCzYEwTRnNhwCr0mvpgx1Nif5cv/0+Efmpo5lzsx5NasYuzLQmwoSbfqu8NN9bVI+YWRCDVcWy3PHZFiAIUTObpPGHyFmoXjnXIppvxQGVYx+yn2aHU6L/LPBtuBHJvw28GPidtnYchNKkWPfnHur26CKM6fmv8bgi8VjpPvmYmxxdN1qTKVZiYv5aielcmrDOvuvC+uk2Vi09EKdEV9VzIEyJnuIOXM0KcFOivzBMia6qf+e3xynR/QDKm1X1j1VVcXNtvrSrEQejNCk2CRasQptp1oSmm7zfOTr2F09vLbBiQS1WANWNMp6DTxMUJpRACkEFkxw1N/tGrbCpsOg/wGzV7M5NU6J/RXaM2pToIvI53JTon072iVOii8ht/jjpMfNp1ms4SNJAu7m2jom2yrwKplmqMpsSZhXi8cRi1Z+7I8UmIISYbUuIOShMeB9NM6mIM+as0k1YwzxbNbvz1miZEn0tjEoabXFSx5wReZ0w8q5gMJTheS/CvIU8ncXLsw7OJgQ1Mnv0cQbOPdvFlOj3++N0HbOGg/BpuiI+lmW/ZZBIWoPKxL+NfFlcNM29b4qgtaXOgFORxvT/LON502UIqEqvpQcGnxLdlyH7vIg8z0fNvpUV9S4OgjQBm4ZK256jXRnGx4g+nZU1gh1I5cqhAgGqusBNbf524IPAm8OU6CLyDX63NwBP8lOifw9V5dd0SvT3+uWL/N++Hfh54D7gL+mInMEB+jSx+mRmsuXBgdy32UXw4FhR823YdxrNsJ2bu5gSXVXvBb6kbxtG9mnqTmlXZ16JrCTOZUJTxvNxJEIK5RGVZ+qDvSpNHtVZDo12E2fsQoPHhEMqBdvTXzkaHJR5Vqr0Ik79M8vEKRm+MOCEzXARc8/2kOUcolUtYVZ/gVPy5H5Om5nWtL1U3SodZmz0qkqzpYqMehNrr+6oo8LelCaNljQRqEl16p9fkadGXW0usz8UsK8n/jTceQdoI1CuOm2mWmqirSaTHkQH56ZYJxLWRZKxfB6dAgG7h/WjF1OkqpMS5zIEBZryzbY73upjDB2Vm8yzLaA0d2AuR8iW/Z4+xDlWWIarTBOPqXmIevkqjRWynqJnO0BKpJRAueq0ESdHIFLwa5qCAVYbxvQfKLpu7j7qkxJm7FnRVCfSbI30CWga86zySFlddZoCBJEkK0y0Y1OldZSgTgxp2Z4V3BjJt5tCzgMiv5imM8wsS8RpMtMCcboIss9gQKnqRnAm6zmGusn6TMMxhok2+TRbQFelsqs0qk/8e1OQYFXnJ90dnWMSKCWLW2/HLlSgD2GGPq8i2Cl6th26fpQCXTLf2hSnyb/pVJc9+TU5UdJz53Bp/tu1J1y/VGX6TKPe9PehcMGE5jACAQG5aqTKUzfFlolTO04PE+2QUOrmtQ4MtjUYkBOgaVqOpv0GxQUMBPS+WiJSiMh7ROQ3/fozfYmc+3zJnNM+x0nLDOVL/neoyv+kf4MqQFAmf3Pb/XkSLqVVXvZdp8ViKVG3DFCidhVCNnRQnLRSZ9zesQwC7bkcCda5Kq/GDfwJ+DHgJ1T1WcCDuNI5nVh1XdrmlYTNneOUJGMPSmszzXaNVR2YaWnb5brPa1eIWYkBR24eBHqRxo+t/nrc6LYwTfoLcCVywJXMeWmfY6VPuHyBZaXJVSfdFn5YV1hCltSmL0p05zd4qUFZtNGf2QRGqjbntc6CsqQKkytLTpC5FoMrjQLWSq/lWND3qvwk8H1U9+OTgIf88FPoKHsjIneKyL0icu9Dn7GdVSC7+hFy4rShPtVEuj05D44o+4TNXrdBWjgjrfecXs8mssTP7HJYtAIq/ZYjwcorJSIvAT6pqu/e5ASqepeq3q6qt9/ypGVS5Eub8oT3UFccWPZv3DYa368Dg8Rls88n5WmXCnlQe4Xke9SuiYnFNNLSTdZX17S417kWzLXw24uoMHMtmNtZjSypksy1aByvH443BFT7LceCPtGz5wPfICJfB1wFbgZeD9wiIjOvNivL3kAo4bTM09p49iTEbNVEEyREz2J42UfWqnBzFVFDqhD1UImcBsFuqU6FSKx51uZftVflMbXXUO/sXGdYFeY6qzn5S0pDvYpNrjQ7xRERog9WXi1V/X5VfZqqPgNXMucPVfXlwDtxJXLAlczpLHvjDsaSkoC/AbIyqLnyhG3pa644KZzPk2+rBwPWNdE2UZ28HFRTeShnLi6HgdP6zUFtzr0ChNJNVis1CioTfRm/HnyVVElStdm2KHk3+gUBjikQsE0/zWuAN4nI64D34ErndCJkOcf+luwJZ8TWiFNgq5sp7FtzfuuKA20D2vAVJ5exSQdnIE5f5TFRJaxvtzeB2kav5opBIMAskskpS2WapcoSCNL03h23OTdtZ7hgSrMWaVT1XcC7/PuP4ApSr438qdpEokCgpqiQ+1tqxlWdoBYBrQautQ2B3keGc6is2QZ349dnPitT4gQfR8MsArKkLEsEScyy9kTOHV4LBT2iyFgfjJwRIM32s9SJEfwZ/xEPW/tbE3Hqvkx9vTqCy0UbIltgXcXJ0TVlejS5SM0oE02t1OkPZElJA0TTKxwvIGxrIstulGcizcZQnCOb1xYORCoypx9SZTHRXGsLEoT3UDfTrD/5rkZ2rhMkCMEAaO5srUULtXL+K5VJfJ2EKOlrwHIeWjNZdm6iXWbzbAi4H9+FMpdGbGpFqJKKRLSoDmpAbIyqoSyZabna1M6HV50BKtasozrOx8kfHFK7ya26WsrnQUm0Msuic+9DwxZhbov42fj9MpJ0mWS7zj+7SNhvscBEtqtkzIQ4aihkWV1SX8cmxHHry2ZZPJ+CkfbhAhYbZzDbFJuYbE2duPnfS0yMlqUh5VRd2hSljSzdI0IHkuXQuXmBMPp4mrnOlswzqCtLSVHd9OHea1SduuIE1KYV17Cvxue7G6y223poTaHpJiLVcuOyqFnwaYLvcq6z2B8zt7OoMCUmTpzUrDTL+WXxnCPc0MfUcdkHe5rdOe+7qCsLpGNknPLkquP2qStOIaV/3z2YLcWYg9BMQ99RgKV+g9dywLLrlSpMiKIBtffpsdLju2OvTmEaFBcsejbqcBPNOtWqVA6XDhJSQsL282C3U9SiSPV09/qsxlWnaJXUWcb3DtXwgc07OneBpc7bmEITJqUtap2Vcy1YWPd6Zk9Y2IKFVstcTVxKdao0ty7qtrBFbQmdmQs1tWUIiPZbeh1r89mdnyQi7xSRR0Tkp7PPvMsfM5+CoxEHU8IpjaClM3ilCkOiOsGPcT6N6whNI2ptOMT6aOmTPyUMBFVpjpiFbACr1cMhxfIgtDQ3b0Xi61Bm24BjZbac3fk68AO4KTWaptV4uZ9yYyX2UI2mOY0kho818WfC71bza0y1zX/GlWvyHkvq32h2DEmq2lAFA0p046TMoZD36AdyhPWonJ4wwZdx6uDTarRunsEyOUY3zRg0gznO7gwgImF255Q0dwA/7N+/BfhpP7vzo8Aficiztm3EHubcbPmB0qdRW6g5IQoY3/di448eFMmRqErsTIdEhz4bcyBqkxcJDFEy9z7klWV5aN4sXVgTCbOwRSdBGh9WY0W1+ivNWLM7N+EXRaQE3gq8zk+P3ohxlUY7ZF9MbT+3rSJPNNuUmrkWiOT6ajx5VGPGQOzD8cetz0bg1Ca9nVyKzfZftQ9CQo3NfBnXjqpvBtI8NKmZZlFhEvOsyyxLz9eFQQnVf2zGzmd3bsHLVfV+EXkcjjSvAH6pbee9hJwDisT/sJpE0ShcXTN10bGQBVygIEknaDDXggKl6TWZmZarzaEV3bBL5DBRdfKxMnNb+CBAERVm7p35eLyMGE0kGEVphu2n2Wp259Ymqt7vXx8WkV/FmYGtpBk5elYNA4B6WBUgnVG4KpiRDRdIU0uyLOBqeO/ydOExknbAeVDpdYnbsvdpmDkd7p2uByUaYhkCA0bPNp7dubVtIjMRebJ/fwK8BPjzrkbsrSytTXrejSTRM3yPv9ioPq6zs1l1Kj+nrjiIM9OqYEE9kdOo+3ypCiJYv74P9WkanRpCysGXiSMsY2i4iGHkECpuu9FH6fnvwkDRM++jhNmdC+AXwuzOwL2qejduiMov+9mdP4sjFgAi8te4QZSnIvJS4EXAR4G3e8IUwB8AP9fVjr0HAlyoOPnhgm8T/ZrwYQtSdWZWJhiRLK6nv+rkDPsUUsb0mrYZ2HKMFU2rDcvOEjQD6sMEqghZiJalZEqPs3yu9u90LPWWN53d2f/tGS2Hfe46bRi9nyYfd55mBxixMZpUOceeKOL8HqcyYZvNiOP7fLCUWsRUm0gABSsukfNElFLhpONeKWQ73Sm13QNOO1Nr41zUJEGA+hCAtDNzbovYAWlVOLez1ht/FSF2TZi+HZfHgtHH0ywPNHM/WDqGJk/KLIPp5Q4RVWfJHPNDB0ySjgP1FJu4jXRCKOVkB3loOelyEuV5Z8t+WH2Yd/xc6Ldp8D26CNCHHIMPO1YuXBrNuEqjdaVJTaA0mbHKFK4yni0+n0xDX4b1PokBWThlwflAhd8vHANZxGPWz1mFncfo4CzEsNA6eZcIkYydCSMzU6WJQwGSVBiLcF4WnTd8X99lJ6ozKc3mCLlnjZCWkZlpB2fu1/iOy7i/V5zS60g6VDpXmzTxc99IhytX75tMtupvoW8m+DOq3RUx1yHD0MSZzLMtkPo0LjqWqE5wgmO0K6iNxvWgPNY/iV00zUY/B1m4vhp8H43O/LZZVJsiRs/cazj+CcQI2r5R9c/MKnIkYeVF8GesW9LoWRf2VvFlIs12qFLgi1okK4ynCdvjMIE0pabm+Kd+jVsvCImcIXOgjNG0EkOhShlMt0MgR3Izpf1NtX2ohgDEDICsH0aT1xx9VWPXxTUuEsbPCPDDcl3ZpQo1wqTrCXlKLbDe90hVB/BmWdrf415Pg7mmhlJcJC0UF3T+kc8yIBTcCO/GQVqDuiKO1Ey1Wv+MLZIAgEG9yqhKHIjWhn3cu+uk/R8LRg4EOMfW+HBvgOurqcgElSIZf4MH8w2zaFCdKtJkqCvOeYy6LTCJ2px0tNNluY1HHMhGXDYkaIYOzWieJSHn0ptp4av2VY3RzLUperY5FGoh0oogRTV5k4YytSEA4Ah14p34uZ1hvWIgZdUfo77QoFecoCImCz0jZS28m+pKSUUmi+6MNpb6HAVNozXTQXQhOdPt05wqo0DZc5q+sTsyJ6XZAi5hs0g3AMRQchy67DsgjVTTCZbizLSa8iT3SFrj2L023Biy4JyCUyWGrMER1RFs3JvJJY4mvftJuDlVmpCcOU9GWS6sU5jS+lGZK0LO8Zz7uIEn0myH5hSP0BdTmW4xXBzTarz6RLOsGkdjUZ9GY2qKA2WMrJVYjBpOEqUJBTj6ptbsCk0ETytm1lJlskhaPX2m3/lGjaJNPs12sCqc2eqUJo2MQeyUNGhNbQpskv5SOgXyUaViKURdKU70b/DDDgQK9Z2iPvxcEqZXrzo5u/ydoZEGAWpVNL0JFkZolho6NI0PBrj3pYrzaUqzkgx7yy+bSLMd6ikfRVY1JqgLVacmrhrZiY+1lWKWVCcoRa441u97InjXPq2VlheyGP+Xrc5fH6kZ2pjWYG6qLJNX3e8dABjsG/TDirINR4fRfZolpfG/oFMZR4AFVZDAiK33+lt38xvRSoWMYFSj34PximOdMqVxhZDQeaIFSBlDz1bUDxMY51q4SjhJwfNkpGZtQqZYXcYpzLmduaTNskj8GpeH1qY0xzSNxTFgL9GzgFRp3GRHNppL1ePQYBVOTInVghml79SsfB7jx81Yn3YTqm663AClCE41SimhhkAIP+/fp4Gk+gz1CZliICRJk0kjZlbritOFvRXt2//lHRSj99Oc2/rgs3oErUNpSj8PjThSGR9ONuL8HyNaRdMsUXHiubxpVugMBE615Fzgqi6wIoSStWOm0oSOzbCEus1pXbM4DMAWMW2m1MqXsd6vsfZAlWYKBGwHBRa2iAGAtL8m77sxqFcAE1WoQDEqWHGEcQ57lRUwt0UkUVCcOS4SN2cW6w4YTKI4y7ObjVlcw52v8l+q4ueV6ZZGykpbn7ms1PbZxA6mHOyhtGMgjO7TLNQkTn6mNlTqUimNe52JsABmpkzI4/yc6rUqR1uIdT6NWVTZzBYK4wg71xBYcFE0K+sUTdkOpU8VCjUL0uqiVuvVNKthACZmBJSJL6NeZbqVZqQv1oZ9n39gjD6epik/KieIqypTmWBAJIS14rcrJ6Z02comnTVZMepnRzPqiVId/1zDlBSug9SqofTTD5qQ0DkS8nBzmA4wzjejstShGZYYMbOCtQZrZeOqL7s034QperY1msbThGIhJSC5wng7yWjlwwRCYX0EzrroF8btF/p1jFYK4/5mKfxXnnvf5pzCpdzogpMRCRPrStfG0chS1CwOB/DvS59jVtrKp1EFWgIBe3/IX1afRkRuAX4eVwNXgX8PfAj4NeAZwF8D36SqD3YdR326x/Lxk1GbCWkMdV/HSF2BrBGMGJeLJq5ARxpAAJ+SY3ynZmKeBcWZa+ELEQqhHpqNftJukzar/LLQqTlLHH/D3Bac2RnnZbUeAgDRt7GCqvcPJ/NsFPStHPF64HdV9R8BXwp8EHgt8A5VfTbwDr/eCRcIqA+esslNoImj62z4pAxr8jmX3VvU162rkr9Q41+L2hDhapnF1xCtSmce2DXCTAUu360+JWC9blvq8GfOv1cbZ6KBWuOVhvpifWRet1iGwJjnGgErlUZEHg/8c+DfAqjqOXAuIncAX+l3eyNu1ufXdB1LaR7z0aQ0kqwb0bjPzLikTRFl4ZVnZp3SzLwvFIZBlyospJpG4opZxL6gE1uCgTkFWJjLuZ+9oFw5dGBb2GiOSo3QZ/bEv844tzP/WnDukzXnZcGiLJxpVnrTLPg3loOdRfkymmfPBD6FKxD9pcC7gVcDT1HVB/w+Hwee0vRhEbkTuBPg9AtvrqWvR7IkZlWpgniiQNWRF8izsJWJBr5mmnG5Z8GHWdgCTOme0KJuDI+6AMHcziiMUs1d6UPRsb9kHFhSf8Zw7uflmQffJVGW9s5M/5DWzJk/tAyAS0iaGfBPge9S1XtE5PVkppiqqkjz88RXfb8L4MZnP1UXZbfSVIGAfF0jmYS6ApXq1GdhHKFOzQJjCxamZGZKrgQC+ZvJDYSzzKXgqszBwHU98SM6nRKEEZ27QOmzq88pmFNw3Z5wZk840xnX7QnXylPObMH1csb18oS5LTgvK6Wx6pRGraClJ1C5efRsp9DLGT37GPAxVb3Hr78FR5pPiMitqvqAiNwKfHLVgVTrswbEd4mzr15pSvWEUXE99X4/9UQJNDMQ1UdCqFm8sogPOYti1GC0YKaWQqwLAOBC0EatDz27zIBdujZh+FnMMUt9Kk19PROJHvw8DWpjJaqLqni5kcN9oh9quzbEStKo6sdF5G9F5B+q6oeAF+Im0fkArtD0j/rXt608FlA2KA2ASDX1RFCX+Eri60TFISpPYVysqzTGE86n46hhJtYFB0zBwrgaYmfGh51lwYkssBgeNWcYtcw5p0SzsZXDYq7qFEZPuK4nPGZPeay8wiPlFc7tjGvlCdfLGdcWJ5zbgrOy4HzhVGaxKGLfDCpgxfkylmGUZgc3+GX0aQC+C/gVX6n9I8C/wz3k3ywir8IVkf6mlUdR/2NnENEYFhWvNu69RFUJ1128f+MUB1Qkfg6vOODNOev7bsJ2O2NmXNbz3BYuEKAz16kZc76EdaYzXxehHG3VJzOr+mVsWm62iOP/1V+3qDY+vKzWh5OHUppdfe3LSBpVfS/QNNnOC9c9YbkIY//TKxlu/GRLVJZ83ZElKI5TIacwhVH/6jpBS2MojPNdzm3JqVlgVRxxgCtmzomUXNEZN5kzAB4z13icnnNlR49Hq8ocuK4zHi5v4OHyKo+UV3msPOXR8pSzcsZji1PmZcHZwg0DOJvPmJdOYcqFN99KpzIEX2YopRkaRxZO7oNx56fRKhpmbRIy9YsNzq0VbOlSQ8pS/KuJZklZJkvaf1OGxTnNc2s4W8y8ieOc6uvlLJpA1+wpj5WnnNkTHrVXuK4n/qlfhYWHRokyV5eJEMwzFwiYcVZWbT0rw7iZKou5LL3aWG+WBYVpI8w2/TMD3ewCQ85Ps6vZnZ8rIu/zn/kpke5iEaOn0WiZZeIuNU9qihOuZtgk8QpLojju1YYcs+DjWBdECP02i8KR7Ny6qTeulScYlDM744qZY1W4xTzGjWbOVV3s4NvDdVUetic8VN7IQ+WNfG5xIw/Nb+DcznhkcYXzsuDa/IRShbO5I05ZGhYL/3ApJT5YYgpDDAb0a4OMrEhDifYOZ3f+WeDbgHtw03i8GPidtnaMPOemoE2BgOyqKkSWBAKFbfEh4E0zt48jUiCPMc4Xsj4StygNs8K6J7cpOPGkmXkz7oo54YpZMC9mPHH2CDfb68zNbkgzV3hUT3jY3sDnFjfy+cVVHp5f5Xo549F5YpaVbnSmU1TBloXzYbwaO7MMZ6KpJ8KhmkHDtWvw2Z195PdmVf1jv/5LwEs5GNLQ1mtdEaG2VVwWQfo3TchUvZdIIhHFWkcg430ca0KWgKH05FEV5/ugXDcuQHBmZ9xYPAGrhnMeBB7B4DLQii3KO5UahgLARxc3c//iCfzd/BY+df44Hi1P+fz8KvOy4Ppi5klTYK1TF+t7/oMPo5EsRIWRYKKtwr5I1f+8+5jd+TZ/nPSYt3U1cvTZnVm03HwS/6t2T0lUixIE5UmCBaKJOrn3YioiiVGMUYyxGKOcLwoKo5wtZsyM5dzOuFrMuWZP+eTsZj5++ng+e/ppTmTBVZm7MTxr9tLFIoCI911O+ej5k/nE/GY+e34Tn7r+BVwvZzx8fsWRZj7DWsN8XqC26sCskSU4/t4ck7C9w+Taa8h3DX+F/c3uvBbGJ02uNDGVpmH/jChLn0EcgQLhwvZYKk08YQRUUasuM9pUnahWiVG30H9zVs7ivDAnUnLFzIGqFlvRgzz1kZgSgwwfP3s8n53fxOfOr/Lw/ApnixnXzk8orWE+n0WHP/b2e9KQ9sUoiE38GJX2G/MQTLbh2rCL2Z3v98fpOmYNo8+EJmX+RHTruupxJNlKo5VXN+EcobwCCZTe78ErDwJnxQmI8sjsCsYon5rdxOms5LQouWE2R0Q5LVyI2qyhNlUhQ9eY89JlVF+bn7iI3qJg7jsqF4vC9WEtfP9Lm6qkZpj1REkCAIfaiThgGk2c3Rl3Y78M+JZsnzC78/+hx+zOPqPl8yLyPFwg4FuB/9bViNF9mranTmtEJ2zOPlc33eJGlkqYSfJ3EU8iQY0j0MKrkfUdoGVpOF9YCmO5NjuJ/T6QjCylmbNNXy90SLpUGKJzv1gUNUWJZMkUpebk58oS/95y8gPBUGTexezOPvL27cD/AG7ABQBagwCwB/NM2nyaNSFRobp2cqHp7IP1V6r7sxRlLssfGeRXr2Uh00iG8OBYIkGuJPG1OuahqszQnZu7mN1ZVe9lOQzdij34NMMespUzDYGFVoSbsmH/vdyM+TmzJ8NSmw6VMAGH3r41Mbp5ZjZVmj4XPttH8u3pEzs83X1gQkofgCvD9mSbzY6zDkI/UwhOFGALtz0W5vG5l4j61+yzDeq4/CW3xI5u7BDYvEg4WKVZeS80/RCtZky1PZAEWyeJCa8LkIXbXsw1Eqg1dJpvyxoeSKDGvZYnghaePDNcCd1IIqntG4i0ZFKmZNLsb13Y080re5nfY3cYlTSCvznXwYobtXYjpz5BGlVKb3rNyLJwfzdzdUQ5c2QpzpTZNesUZ2Hr5+p7DyQdsDozqIHFDYbyilCeCOUVT54TQY0nkbhtGsmThNWFZTJ1qVAbxryHB/ZpDgF7VZpest1lz3eYXVCZVdHECmaX9WqidbKYUpldU4rrlpPHFsw+dwbWRtL4jNP+39cN+nEfnRkwhsXjrzC/cUZ51bC4QbAzobyiqIHy1BHEzhKzzZtukiiQT+12gUKThdkPEJN5tiXMUj+NR4/OuTbC5AqQKkpKFGz13izUvZ57slxXzEI5ebikuLageOQM87lHHUmsJ43dIIphnDMjxoARZqrI/JTFfIZoQXni2KAFoI48qFMeLXAzigTfJwQDvVnmIuz+enY5D/tm1ESaLaDuCd9nv4BGouR/07rKpGSpOfVWnXlooZi79eLMbZs9ZjHnlpMHr2MeuYY8/Cjlgw8NO9+eEYqzc+RxN2FuvgGxVzFXDKJgC0GsI0tpHXl0VpEnBA2CAkUCBfII0JYft+ebdlKaLVCLRLWh5QLn5Kmpiy6rS80Us1o3y1KlWbj35txizktkXiLzBVra4SeotIqWFlm485i5RQ2YueuvcZE1cRNYF75qaItKqHHfdak/9xBv0ENs0xYY3TxbFb5t7YPI/RZa1MX7TaJaM8tMWSdLce6GQxfXLGahFNcWmPMFcv0Mrp/B/Hyw71zDYgHXz5CTGcXpDLEzMM63AYMpXCFB9dVEtNBY06ymOkUWih51OOEa6POgPDLsjzQp1ggfp6RpIos7hy6Fk52JptFMi+9L994snMMvpaUjVWl7qEWtdf5RqUhpXRvE+VYWfH6euvB4EjFzjj/VXL3hUhxwhGrqp9kWLT5NZw/3SoXRKoyc+DC1vpeELNEcmyuysBTXS3fTni+Q8wUsSlh482wXsAqlMwHN+cLVZp8ZxBq0EDeZtbgyTeBmSXClZZ2fYxDCbPDx9YAjZwD7LyY9LMZXmq6nYlMAoElxMrLUzbKUPJr4NpqYZ5r4Oe5p74or253/wKquPluIyom1rp2+PUFhRMJrpS4I+Mmpq2ia+i6advdn75iUZgtIpjT9FEbjemu0rIEswTyLZCkVWVSkMefOHDPnpTOT5gvEqwxluVl4uQ+srcLYixIxBjN3VdZMYV3EbO4JJOJJLLGLyODKXYmfIj413do6OPc9CO1QTcdNMb55tsKnSUni1llez53/SBpdWo+EiQTy/szCP+FLRYKPobpbfyZ+3+p84TW0Dd9OELdN3Ksacb5M7sNZXCeotnDmAG7YKRCwJaR0v+KqPK7GkDKwZJKlytKiNGaR+DRzG4ljFhYpvU9T2kQFdPd2uLpzqqkyDtyr8QEMRUo/87SvgaC+01M0TGIl0Vw7ZEyk2QaJCjT9LX+/5PyHY2RkqQcCNEm+9E/v4LsEf2bhCIN1kav41C+ti27tNHrmSBl9m6A2pR9YVmpNaRCJoXMV/0q4ERXRijgH6dMoUyBgW0iSsNnm03SFmRujZf6HiWZZqVWIuYEw0QkPhCkTc2msHziomXpSGx8MANfGmet4tYAYoupomE4hhtchFBrZomDOTjEFAraAaEOWsy7vk26vRcn89nal0agwBPKUwa/J+kUWibqkznl83Z1NoeoLU9nKLFSLM9cw0behEGdWGnyFHYmmTjTLrCdLIM4h3qCH2KYtsP/OzaYgQLK95r9Ao7KEEHSMzOWKEx1+rZljUmp0wkdVGaCqXk6lNlbBhJB08j0VZ4aFBwPOlAvrIUhwiM7N1Lm5LbQKBEA/86x6X1eVWl9NSe1GC/0eJvgHQXHKyo8J5IkBgLJ0voxmrM7Xt4EkuS6pP2MVEedLSenbSxIYME5NKvPMZwtoOKSieTnfQ0F4KF0g7F9pyMizRBytfa6mLGXYpnUFKrXqHwg+TXhil6kvkSlM+HF3qTjWQpHNGp1O3S6uneofAI7Igqiv1Ra+lyZPcf+9dUOpkV0r7MXizCH4NLq0j9uerGv1w8Yx/FoRIzzJzCLs402cYJJZ78NEdVGXLpOQR631/Scj/MLp+Bxrom8l1vk0MSBgfajceH9GAul9wUOf1BmvWfL0OSST6JDaMgT2UPdsXfNMq/eWZcIEhQkdm5r4BGFb8F389pov0eLLNJpqQyLtC8peJahG8v00fjepJWmKvxQxnab2JXbX/N4Iv80Fwug+TVSa7DouZQLY+n45UaKPE3LIgu0cyVJ3/ONrEmJ2OWehZ959ZtcZzks5/L6dKtZHwiyUxs30VjomhMwGSpfQWfNtSqkOeYg+DRwGeQfEyErT4hQmm7rI4tapIkoxu0Cj/1KpB9XTPGxP/l5rRx5e3mG4OUXs4PRwfkvyPrbHmWTxOsQHR+rTHO69eSnNMxH5T8B/wP0u78PNuXkr8CbcNAbvBl6hqp0jt/KEzfrf6krTFBCIypL6Oan55YkgC+v39U/vMHbFpsqSdGYGU2xMMyI1/dJ2WHEqiEEKr5wmZAeEh0rVXxNC0G5ls6bs+qa+aNGzleP9ROQ24D8Ct6vql+Bq6IYZpn5CVZ8FPIibgWolgs8Rl5CoGJWCypZfMrO05tfE7Wk4OpJK4whOd4yMJFAvmhG3jaAy+fmTvqLod4X9alkPRFM09FPl6tO52OYlPcbSsvV3XWM5EvQ1z2bADSIyB24EHgBeQFWx/Y242ad+tvMonhjNuWfB1PLrNltPlSUxtSJBFNfbnytM0okZVKbmy6R+zC4d/y6kIXBxxHHRPFPPfhY8oQQpcEXcbRi/GcYHNB1/rC+yDBcWPyJG9MBK0qjq/SLyX4G/Aa4Bv4czxx5SjRNTts4eJSJ3AncCXLl6S2VC5dcxkfClXDNYIgtUsl8N4gpP8ORpHDIAEtVKo1UabtaRTIg0hSa0oWEn3FQaGqNmoY9GE8WBEDnEbxzjG2yAPT2LdoU+5tkTcPMYPhN4KnATbiLPXlDVu1T1dlW9/eTkpqQ3XuP7Kt0lLHXTrTLPqIIAceyM/3xQnGR8TFSYxAxqDDMHhakFB0YgUXrerG1xdGdU52q9Mm+Ta9Zhfm2zDIElk7xl6XWsDWd39n/7fr/9QyLytcn2v/azO783m76wEX3Ms68G/kpVP+VP8OvA84FbRGTm1Wbl7FEBpmy5OLmypHY9VdSstl/wX7wfUplhyVM8McMkTcrU5D1UN/Cuw84ZovKkJlrIDPC5aGFbHJCG+GwI79eEyjVygJGqAf2VbWZ3FpHn4HzxL8Y9/P9ARP6BqoZOkK9S1bZ5OWvoU/jnb4DniciNfn71F+Jm030nbqYpcDNPvW3VgcJY9nyp0lz8DZOGhrVOGAk+SQwgNH1Go3lWI0xKhqg8Nm4fkyz+hPXX/D3UOkGldj2oX78+Tv0my/Zfsm4xdCw9EGd39pHaMLtzijtwPja42Z1f6O/bO4A3qeqZqv4VcJ8/3tro49PcIyJvAf4MWADvAe4Cfgt4k4i8zm97w8qzBR/Dv09Rzz8L+4S/aT0wEJ/K/u8h4hXSZ+J64viH7U0BABg33Gxt3beJppfG76MYl/EsLrs5JnGGvqkQEAj1nTfwaUZTpeEeRtvM7nwb8MfZZ4MfrsDviYgC/z2bUXoJvaJnqvpDwA9lmz/CJkzNo2TZdvcev08w1bTWX1MLJ9eOmahJnozZ9sMlhNw3xCpaZHd+VM0kfcZ3dgbliSk3655wjK+ta/lGq6ZE3xX+mQ94fRHw+yLyF6r6v9p23uvQgABpeOLnWc4VGaibXiTHzBSm5sek/o1/0qdRM40+zoihHuvvKGuJlWdiuy1YiYqDqA87K5D4NgugcNuW5hs9FPRXmlVTom8zu3PrZ1U1vH5SRH4DJwatpBm9mGlj1MRWS92HSRbvt0jDtroTnZo5mfJARSCrifN/ODFRSYgct2n+mlwnQvRMa4TpG7EaIqK1Em3+Ur6sRpzdWUROcY793dk+YXZnqM/ufDfwMh9deybwbOBPROQmEXkcgIjcBLwI+POuRuxXafIfpaY47iU3wWqmWW1bnWQ1hcn+FjOYw2uOXZprVqGQuqJlyhnfl4oQhkKDYtz1M4lvUxLLbMqQTsqAh5KB1Hub2Z39fm/GBbEWwHeoaikiTwF+w8UKmAG/qqq/29WO0RM2OyNF0E0W/3fJj9FEmCZCtv14KcHGhPo4cVxX4vBLq+628KiGAHhTTN0+YvwwaDYfhOaOv6Pvrgzaubnl7M4/AvxItu0jwJeu04bxx9PUhjs3/FANJKnt2+HsS6i/nPo2YT0tBphGrJZOH86zO5MtRs7S9TAgTTyBjYlRNBH1QwUsiHFJM963CR2eanoozR78HWFAM+9AMH5Z2vSmz9HbPMsI1BQha1WqzCSz43ZmtqKtveG9KrG4hhLNVLUgxs9v43c/vA7OQ2vQdhhfaZbGrtRXl0yv9LWFLEvVZHKTzZ+zNhozJ0+X+bYr+ATN6OcEf82GaQSdAqkxVZaAuvlqKEEKXMJm+N6JfB0UcSbSbAHVxpBz+Fvjep7YmLxKm+IEsywQITXL4v4JkdqaO3BAQEzmv6TnCtEqH3oW689vIIzoDKYaBaCClupMNyNxlGfjefd5zw7s0xwC9lojoHVberO2mC15tvKSj9JEmBAxi+Hn8DrOr6pWkWJpoydwUbU7IU+ovqnetHVlnnDDAqj2USOHpS4JhoqeHQr2EAhouYBtSpNsX1KW8D5Xo7KM670JY1vMth3CBQSSaJkPANTMyrBufTZm6ZTG1RWgGmezAcZx0Dt82CPF6P00nRewJRwtHSZabf9svdX0ajHN9hYQ0AbTKlPR+tgafLgZZ6ahVX/NKoz90F/1mx8hxu+n6ZDqpUzXdYICUB17yYdZzgLQFoLuHBpCy9YVDQz9MWpBDbEWmkkCAl6dhVClhthxo6JVLefkaxyUqXaxrLM9mGcJmivTrDbTGkc9Jp2TdUJkZln6+abzjQmbEkeq/hvvp2B8QCAoSEN4PdZIIzH1oLbPvjH102wDTToga9s7yNN2g6ch4pwsTQoT1kO/TKp4dsdz0nSgThT1kzblqurbWrg0GhV8hycI1g0fKNbI1hw7q3sizZZo9TM6Imjp3/MbKf45IYzWyeTer2cjDB1uXnEyajkzcXv+kFDXV0NFNk0iaH0xqumm2h78OVKM79OsuoBt5hksKYt720GWTGHi/ulxaoozcuQsNznF+zOlD4m5qQLcvDWFqdpqBEpx5tsCVKTq41k+0SjfpxOH0IYBsT+fpm8ULQ8cNPktTaqQj/3PP9NEvH0gjKUxlqaStWm9gJAlAL5SjQ9DC7j3Jr0m6zVjp37HRJotoHR3JDb9re3mbiDFSoXJz5H6Mgnxdm6adZpjFlWp/Jzw6qtuqlRqLVINUnNTpEv/G3SsG9mHxy8SDifk3PIjtirKKrIk6zXCtPhEe4HvvHRZzjhnw1AFBELnZkinqYXX3TY3pQhoQfM13PvXHM/kHQsHEQhoNY+WggENREn3a6hf1qgwbecZ4cdVPwYmJmk2tSX4NmGYgIi/+a1TIAn1AkIN51DSqevEe3raK1MgYCtocMR7/IBL0/h1RNf6kCVbbwxPk5lmuyCRZH6LJ4Uag7g8/2UzrSx9fw4x9CWlrWZ6FnHBgTbs2zyafJot0aUeORrNjZabug9hDgyq3hcJsOoDAkW930YSXyWYaMa4js8QPy4bomervvdofs1hXv9NsQefZpUdQfOTMSdXE0ni3xrIkitM23nazjcgnImWZTb7cTOSqEkw0xRbDymHSuhS5QCokSozugX7mfJCJ9IMhk1MtIbPaQMxlj+z4jhNptlIqGUEpO1JzDSsqfyfqKJUwQGI6TadxNjHzbsqYnqEGD/kXK5Qmo4fvVVRqh1a/74qCrdzX6YJVn0d5hAls37SjAIogcJ33Vi0BEzSiWmMCwSEY0g10VP1PXoGWHaNSWm2Qc8gQPqJRr+mP1mWjrEPcvRBMlV69HVCRrQnDtY4NQn7Gx9UaOgXddtbrvWoN7FO0bOtoC0kSLFKynuQaFmR2oMPSybZronkTa/o13giQNIrn/o3lE5RxFfaVEHFLKuOpd2n6bqmuyaQgh7Sw2kA7L+wxip0/airyALrRev2jZp6eOIE53+V6kA08WrouH6jpQ/tO+Q9MA6ic3MJPYjVu0MUGokyusLk5wr9NYEM4AeclY4oiY+DBvK0qA5UvlFAz+swCiafZgvkWcW9PtKHZG22ew+y7BExOyAlUUCD6gB+EqdANK1S2Gzy+VUdw/EcI1yLDX7zQ8eopFHWNAn6/qgrnqCtRFm3Y3UIpOZW9G1kWXFIqriYarhApTxaO049w7me5VDDPh4ak9JsiXV/tKEHj6063s4dY20kDuA7PKnIgw8nB5MNIoHUm3C1YEFAnqiaYvSnvqKruhmODHsogL75j7aRadX3fGM+DRuI495m5IEl9an5PIE8tR3C51akEI1FHmUKBGwF3cKnGMqJPRRTIbQj9Mekf0r9E0pnviUPaxXj+j58kqbkoea85sG+b9pDjlhugD1Ez0a4gIdCjD5IyRO3dZEIakQSU5tkY6l01Z6hHFbwZQgcZsj5MqJz+PeyT7DEB505ApVlv6TYsaDbmeSHiL3WPZswHGKGwQHeoBctECBjFpUQkU8BjwKfHu2k2+HJHE9bYTft/fuq+oWbflhEfhfXrj74tKq+eNNzjYVRSQMgIveumMH3YHBMbYXja++xYvTZnSdMOHZMpJkwYU3sgzR37eGcm+KY2grH196jxOg+zYQJx47JPJswYU1MpJkwYU2MRhoRebGIfEhE7hOR14513r4QkaeLyDtF5AMi8n4RebXf/kQR+X0R+bB/fcK+2xogIoWIvEdEftOvP1NE7vHX+NdE5HTfbbyIGIU0IlIAPwP8C+A5wDeLyHPGOPcaWADfq6rPAZ4HfIdv42uBd6jqs4F3+PVDwauBDybrPwb8hKo+C3gQeNVeWnXBMZbSfDlwn6p+RFXPgTcBd4x07l5Q1QdU9c/8+4dxN+NtuHa+0e/2RuCle2lgBhF5GvD1wM/7dQFeALzF73Iwbb1oGIs0twF/m6x/zG87SIjIM4AvA+4BnqKqD/g/fRx4yr7aleEnge+jmhfgScBDqrrw6wd9jY8ZUyAgg4h8AfBW4LtV9fPp39TF5/ceoxeRlwCfVNV377stlxFjZTnfDzw9WX+a33ZQEJETHGF+RVV/3W/+hIjcqqoPiMitwCf318KI5wPfICJfB1wFbgZeD9wiIjOvNgd5jS8CxlKaPwWe7aM7p8DLgLtHOncveJ/gDcAHVfXHkz/dDbzSv38l8Lax25ZDVb9fVZ+mqs/AXcs/VNWXA+8EvtHvdhBtvYgYhTT+yfedwNtxDvabVfX9Y5x7DTwfeAXwAhF5r1++DvhR4GtE5MPAV/v1Q8VrgO8RkftwPs4b9tyeC4kpjWbChDUxBQImTFgTE2kmTFgTE2kmTFgTE2kmTFgTE2kmTFgTE2kmTFgTE2kmTFgT/x891D9wYiXyNQAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(T.reshape(nz, nx))\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "df47301c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMYAAAD7CAYAAAAxf+suAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlRklEQVR4nO2de6wkV33nP9/qvnfueMb2+IWZeEzGgGHXWW1g12KJnES8sjgkimHFIlCEvLuWvCtBloSsgtmVlkhRJFglIWgVsZlgEq/EBhCPxSIIBxxYlD/WwQQLbIyXiRPA48fgFx7PeO693fXbP+pUd3XdenV3VXVVz/lclW7X65zT1edbv/P7nVOnZGZ4PJ5ZglUXwOPpIl4YHk8GXhgeTwZeGB5PBl4YHk8GXhgeTwZLCUPS9ZIekHRc0i11FcrjWTVatB9D0gD4f8AvAA8BXwfeZmbfqa94Hs9qGC5x7iuA42b2IICkjwM3ALnCGJx/wIaXXLRElp622fn+icfN7LJFz3/9qw/YE0+OKx37jW9t32Fm1y+aV50sI4wrgB8m1h8C/kX6IEk3AzcDDC45xPP/668tkaWnbX5w0y3fX+b8J54c8zd3vKDSsYPD37t0mbzqZBlhVMLMjgHHAPYdPeLHn5xjGBASrroYc7OMME4AVybWj7htHs8Ew9i1ak2pLrFMVOrrwNWSrpK0CbwVuL2eYnnWibDiXxmSPirppKR7E9t+W9IJSfe45Q11lHlhYZjZCHgncAdwP/BJM7uvjkJ51gfDGFu1pQJ/BmQ55x80s5e55Qt1lHspH8MVopaCeNaXkHpcSzP7mqSjtSRWgu/57gqm4qWnGDDGKi3ApZLuTiw3V8zmnZK+5ZpatfQHNB6V8hQwT4Wvcqy6GfSbw2I8bmbXzpn8h4HfIdLg7wC/D/y7OdPYgxfGqmjCCnTQshiw2+BTomb2WPxZ0p8An68jXS+MVdDBCtwUNm0mNYKkw2b2iFt9E3Bv0fFV8cJom3NIFAAYjGvShaQ/B15F5Is8BLwPeJWkl0U58Q/Av68jLy8MT6NEPd81pWX2tozNt9aU/AxeGJ6GEWP6ZyW9MDyNEjnfXhgezwxRP4YXhsezh9BbDI9nFm8xPJ4MDDHu4cgjLwxP4/imlMeTwhA7Nlh1MebGC8PTKFEHn29KeTx78M63p1+0MErdTIzNWwxPV+jQoxmhtxieUmTNjrDtkCAgdr77V836V2LPLB0TQhrvfHvap+OiiBn7fgxPJZZtTvVEEOB7vj1t0SNRxIQ+KuVplCZE0bDQokGEXhj9pUrLps5KNG9zapm8V2hlDLHrh4T0iEWa+Olz2qpwi+TTkSaXGb6Dr/PUHRxpQyjzpNkRMcyiXnbw9U/Ki9LGb7PK37+TonA+hgWVljJyZju/WNKXJH3P/a9lis71F4Zot8KK6nmWTalZpbJbxeNK02luztwxQaWlAn/G3tnObwHuNLOrgTvd+tKstzBWbcGbzH8ZQbQ4abQhQqu2lKZl9jXgydTmG4Db3OfbgDfWUe719TFWLYqYuBx5lXiRzr5KlqQbFyCaPqfRanZ5YorOR4HL60h0fYWxrvRIFBFzTbh2qaS7E+vH3DscK2FmJtUz5ft6CmOZerHMuUU/SZnlmCetzOO7JIYpxlw934u8BuCxeGJnSYeBk3Oen8l6+xjzUIeT3rSjnyeWjooiZuysRtmyILcDN7rPNwKfq6PM62Ux6ui0q4NkmunKnGU5qvgZWaLouCAgeoKvrrFSObOdvx/4pKSbgO8Db6kjr1JhSLoS+J9ETo0Rtfs+JOli4BPAUaLp199iZk/VUahWaKtOVW1CzT1EpPuigNj5rmdISM5s5wCvrSWDBFWkPAJ+08yuAV4JvEPSNSwaP+7Cu+WqZiubf6ma5zxfPZ1sT0QRodo6+Nqk1GK4UNgj7vMpSfcDVxDFj1/lDrsN+CrwnrlLEP/Ibb0/rqxOLVuO9PnJSixSTSiadbI70BseOd99EnLEXD6Ge5Xsy4G7qBg/dm/evBlgcPGh/MSXFUilnuaifdXzVSqdwlfMxelOvh/Li6OMDggiyVoPO5d0EPg08Otm9owStaMofuzi0McA9h09Uv6TmZqxHnmiKMgrLYB5j5sRTFIg84Ruk8cU3Xk7JoaYuOe7b1QShqQNIlF8zMw+4zY3Ej8G5hdHafMob3t2HlUFUUacTuWXltbdC14l3xbq7FpOhqDINNwK3G9mf5DYFceP30+N8eMJdVmOiqLIFEMd+Ztm0jZSlsOo1pxKC2aPQ75o+RY8r2ryBrvhGgoDuA54O/BtSfe4bf+ZhuLHrVClwpcdkyWkrFPSPkZye5Y4mho7tSKiptQaCsPM/pr8+27t8ePZzJe0GhWswIylSOdVUD8zXSpFHVrZmDsk2m/ziCDPWiSKoCXa8dZwRNDPXdt1lhRFlfFp6WMsHY1yYpAS4kDzj6GyuNgZBU1HvcrK7NJoQiDnRLh2JSxqNfZ0qOWIIrk92eUwh/UoLIZsWlEnFmVqPSyxD5g5NvcF2ZYQRKlvklWonLI2UoHXtCnVCWqOUmU72sn9tmdb4bkFWOw/QFShZVjCQkjOIa/UpKJgIGHBeemk02JsmD4+890PYSxLQlRFliItiD0iWMByTdJwjrYZCCeOmQhVyhkvStM026QqE0Ve06oFgURRKT99TvfIEsXM/njfVBSZ4qFa/cn3McwFn+LKb1ioSbkmliXtlCebTJYQRZxNOC1Vti+fsn5JC0bBtppY6w6+WmnZjGcyI4IMKyFLbnLbS2/j2ZsT4dq43k8FYvGuqb9RZjEsuWia7cQZjzN2IplYI8rFQcb2GvBNqXmYM3KS62fsae4kPxc0oWaOm90vXIWu0gmYzj5XQBZZDzNkir5+kBAD5kSRE2VKWQvFliJMfK04wOXEFYkc12yLEyJhmTKLWas4fFRqGRo05UVElT9emVqJWBSxEOZ5jDgpnr3NqulCOGs5rMoFSFiKqUDifZq6Ke57WZEAyvKp8ffwUallqGLKFwzdzliLZGVPWIqkIIIgEsW8z9XHxycFlU4jDAPCUIQShMLCAJOQu8VbfKtPR4snYhAKnSDC6DOxcTNQACZhgXPyE5fMSF2/Fm5IZmLkhbFiin7kwoeIbMZKxKIIgr13/Oj/bEbJyh+fG28bTMQS/R8B8fNh5mqt5JpWSe9ZiVqbcLYnWYWRhZgIYxztC4dR8oot06JWo0Z8U6oO6nICU3fvmV0JazG5ywdRZR4MQoLACGRsDseE5u7yzi+IO+iS4lCi8sdbB0GIZGwM3P8gJJCxPR6wMxowGkchzFByack1g3J8byOKQDlrEQsh2BEaw2An2jbegvE+w4aR5VBozP1kqfcxOiiMJkkKIeULKGEtAhnDwZjN4YhxGLBL1JQZh0F59Ino/EEstCBkIGP/cJdBEAKbmEViGAcBZnKdfnEB997ekz3SSkSjYoEEIxg8B8FuZDG0ARbGzSp3YlWr0URUygujRpYx/+nKnze8Q4aCMBJDEDWdNocjNodjLtx3lsv2P8uZ0QYnz5zPKAzYHUeWYhwGico8TS5QLAYniCDkvI1d9g1GvPT8x7hk4zTfPX05j525gGeDzYkowlAEYTRWyfI6+eJmUxhZiGAXNBbDMzA8C5d++zk2Tj7LE9deyqmjitIakAhRZV2nkvUa6Gs/Rv+8okXZUwlm40BKVOphELJ/uMvFm6c5uLHN0FXysiTjdCbpARvBmM1gxPM2T3Fk8wku3jzD/uEum4MxgWtuzfSl5IVQ3f84fqBQaBQ1oQZnjY2Hf0z44A/YfDaMmllhRkKFflb+rmUJUaWlS6yugw/Kf4yancZ0JCraFlmLwSBkGIRcsHWWS7ZO88TZA/z1wy9kdzxge3sI8cTD8UjUrC6VjL6QxwcHGAxCHjp1iPM2djly8GleesFjfP/MxZwdDTETu6I82pboAdeYSBTbsP9HIYMd44lXXk74c5czOk8MzoINYDwpWKJJlRZf41EpGNX4oJKkfwBOEX290QIzF1ZitU2phWLsi4Vs08QVd3KHBwaBcf7GNs/b9ywnnr2Qpx67oHTcUlXOcB4ExqGrn+NFF53k9GgfDwcXMgpi/ybnQiTyn3TiOXEMdmDrqTEKjcf/6YDt5++y9fAGW0/AeF/ypMSXTtLSTbqBptSrzezxuhNNsnofo8GhCFG61Wq2AeNQ/Hhni0AXcGZ7s/6ymHjy7Hnce/oKHn7uAnbGA0bjwDnjOedo5vQoahU4i7AFZ543BINwnwsH7ze2LxTj/YYNLTo2sFlL0WKrpa8+xuqFEdPkEIWC86cOcMAIeObsPp7b3eC57Y0lM83mqTP7uW9wmLOjYRS2dR1+k97sIuLKHUC4EQnk9OHIWR/vC6P/B0LCTUEsCJh6kjMiq8kUViD/qcY9VJnt3IC/dLPS/PE8s6HPQ3eEASvriJoO1RChC89GlbX+vMbjgO3RkJ3RkFE4tRaW8CH2MPFBor4OBDaIrEA4joQx6Vx21mTqpTPtJ2xRDEnmcKyrzHb+s2Z2QtLzgC9J+q57oUytdEsYeSwgmHSTPWt4RzyoDyJBmEXhUoUB1tDMFqNRwOntTUZhwHgcWQsLg0TnIbnisCDqiDRF7XZZopk0iE6KBJOIclWhQb2Y1etjmNkJ9/+kpM8CrwBqF0b3wrWrualNs3d376aKYWHA2HXwxUu0I/U/eU5C1JawAhb7G+lfMXBLymLkUuWYhYn6faospSlJBySdH38G/iVwb/FZi9GuxXC9tpVeyjjvD5XXmVeWhUUD98JQrudb7jntRnq7okF1o0E0kDC2FiYsjPfHx6ZNHhBE5TIjshBxZ6AS/sTMIEHNnp9TptJjlmQOH6OMy4HPugjeEPhfZvbFuhJPspqm1MxAuQZZMP3CKNGSmEX+i2X1iZT5GQktxH7HnhGzyV7zxENSuWT1stdInWOlzOxB4KdrSayEFfdjtCSQvOyZtRqN1xKXadKviKNREzFmVaK4WLEYYmUY2f5EQOpr5FiRdPpNkLSCPaK7zndRc6qmTj6YhmujzwBxU6qW5DPzI/ZhkhGp0o5vY/q8eMJyREWeku6rSOs9r5nYoDi6NtyjCt0QRl5FT4qj6VCuK0MsjubyYTK0xMKEKOKmVV4zCjBcD7mA0KbPdKeLnGxCJa1NlX6SmsVhzvnuG90QBtRqBebOOmE1ms9MUwc6sS29nkuyokcHFxxr2eIoCi4kb0Q14ZtSy1LXOKhFb/gzVqM5zI13iqfPmTrfTCtkXl8GrkmVR9yxl246zSMOqNV61BiVao1uCaNNSkTY1I9prv8tGZUqrqCpik3C3yiijopdwyWY6dXvEd0TRlGFbdrPKMu/riwSvkSmtcgj8d2NaBqe7IdCmM86NIwfRFgXycrZgBjM9ja3WvMz4mEoyfDsTBOqJEKVFEcyepb1lOKeiBQrEYz3MXpCXa8SWwgnhlxRZFHFH8j6PE8aDWFuHFrf6IcwGmxCZSbdcAXaIwqbaSNlk67Y6f6KrGNb6K+sQgeKMDf9EEaDtH4PzbIUyX1FzHPXX7FfMaGnzndlGydpIOmbkj7v1q+SdJek45I+IaneR95WeTGb9DVM0XyzoaZLnuOtxJIsW+YyPaewqbiKy2oVlw4xT+PvXcD9ifUPAB80sxcDTwE31VmwtWYR0StjyT02UcvKsmoh4JAcYl+0dIlKwpB0BPgl4CNuXcBrgE+5Q24D3lgpx3nuEEUXK3O8T0aieY5p6tj2er6p7w6Za1EKrEbLd2aDyRD7sqVLVPUx/hD4LeB8t34J8LSZjdz6Q8AVWSdKuhm4GWBw8aFFyzkl7S230bfhWZx0cKEnlFoMSb8MnDSzbyySgZkdM7NrzezawcEDqZ20fgcrDdUmJkxrrhAl++ZZ0ljkt8x0HMbbV0Ty9QdFS5eoYjGuA35F0huALeAC4EPAIUlDZzWOACcWLkWVu/66WIZkmDVdGYr6INLkVfS4kmV1FBaNw2qSjlX6KpRaDDN7r5kdMbOjwFuBvzKzXwW+ArzZHXYj8LnGSrlupKfkLHOmJ0PSEwvk+ytJUay8KVPN8e6l853De4B3SzpO5HPculRJ8u4qNVywsmbRdO7YuBm1dJYFmbnQahxmnZesEGe6g3BVliGPHoZr5+rgM7OvAl91nx8kmrpkPoqaRC03l1Z6j5pnvFJWpUmLIisokdy/KozJ8Po+saLJENz/qtcrHlS4jHDyhiG1aS2SZTHKxVEkCMg5J+G8VBlq0gr9E8ZqR3et6sdKvQKg3bxT//NYRBTzpNUmNTalJF0v6QE34uKWRsrLqoUB+ZGTec6pel4Ji7yldbGMMnqmi/JMV5wqoiiIWrVOTcKQNAD+CPhF4BrgbZKuaaLIqxcGLCaOOinrMW+bou+/ap9hXiah4wpLOa8AjpvZg2a2A3wcuKGJYrcsDM28T26GOsRQcnEVR4Py9tHO0JBJOarUhVKfIuucrvgWrgjVO/gulXR3Yrk5ldQVwA8T67kjLpZlJc63TO3MvF1S8bRoyLSvrOqrVo9KVZntvBVW9jxGbDlmBFIYys2ITBU9n71Ahc96o2vtpNOefJeMY6v6FXnXrSOar/HecwK4MrG+3IiLAlbuYzTatCrKt6RJ1WzmVk8+Han4hVR1vKt9l68DV7tngTaJRmLcXn+hu/gE34xFoMCCFOwrYc9NW7NCUV0Vt7QgFTr3oLpvEW+ft+yNPu1X2bEuxcxGkt4J3AEMgI+a2X21JJ6iE8LY43PkNg0Kmk6L5q3sz62TNahwGbpkTWosi5l9AfhCfSlms5rXGWdUwNYc8tlM925apTO+yJ27rmE0TVqNva9I7zydGhJSWRxVHPA5rEu66dTssxhxmNYSAYUFxLDM/qKy1S2OuB+jZ3T/Pd81HTPpp5hs6Ji1WGP6eFlX72OkKvXEatTpZxSIJi2GifVoqaOv9Mm19N22oeEwE/JeVbYMPRTGysO1QHGEpWh/5jlFKshPqJfWouki9/Ga1MQKXk5JfmdU3Nyp4mvk+ClVSHfkJYeDpEO3TSBZ555Ya5I+6qtTzvfSaZ47da0/GPMMCekM3XK+syp35rbl+zOKrIXaeF4j2W+xSB9Gn+7CfSqro3POd0zlTr8a6KV/0SP6eHm753wnPk/GURVd2GSc3B23Z06lHGaHgezdtlI6UoxaqG+sVGus3mJkUcFBX5Ss8Gz6szKO8yxBDy/laoaEQDVfojAd1zBfRCgl89YGXhS10ddHXlZnMSo62qUdfpNzM57XKCC7Y695pzvqPKTaXXRdQro+KjUnRWHbqhakAae8s0/2dbBIVejipSyj0853/HnGCS8aEpEMYuWlmWAanmViLVbajFoXC5HGO981UdlaLD+yNk3cpFpL6n7mowodNb5ldEcY8zSJivyTOZtWyVG3sbUINKcv4Cmmh9dwBVGpgrt5QYXPHT9VmF41y5EVsm0ayVxgbYnnc3tS4dTDB5VW42NUmWAr2e5M+hpZnXeTbYk0SzoH04MFA9eECtTh5tSeh9Vzts+ThieTFUelMu7odUSZMtIwmzabyvD9GDXTw8u5eh9jEXEU+hh70zNTwXQ5lrPM9zXmYc8DSumi9bAi5dJT57sj4dqCWpjVnMo7P1XZJq/cqkigyFr0ymL0oWnUw3BtN4RRRp5fMbOuvcfvSWdvLYr7L6YPL82ut04b2bYtphaEIem3JZ2QdI9b3rBMeqtvSsWkm0AFzamZCFX6uOR6Ms2S5lnsdPfOYixKS1Et0WpU6oNm9nt1JNQdYcBc4sjcn/Qz3ABDi9u4cWsrw9+Y9l9Me78bsxhFoepK589xbBdYZx9D0iFJn5L0XUn3S/oZSRdL+pKk77n/FzVa0iq+RvK49GfYMyNHctBg7HQHMgZBuNpZCcvIC9t2lepNqbLXAJTxTknfkvTRZetjVR/jQ8AXzewfAT8N3A/cAtxpZlcDd7r1Zsm688z0dyTHVGm6Lf263AxRBYlmVKMWA+q5hcY988um0QbVhfG4mV2bWI7NFFf6sqR7M5YbgA8DLwJeBjwC/P4yRS5tSkm6EPh54N8AuDfZ7LjCvModdhvR21zfs0xhgOzmFGQ2mTJ9jeTwc7dt0pyy7NZZcsj5IAgZKFyN850XTCiiDnE0/FXrupRm9rpK+Ul/Anx+mbyqWIyrgB8Bfyrpm5I+IukAcLmZPeKOeRS4PKeQN8emcfzs6WqlqjQDeMV9Bcelp7CJrcUwCNuZDOFcoZ2o1OHE6puAe5dJr4owhsA/Az5sZi8HTpNqNplZ7lczs2OxaRwcPLBMWV2Cs6uTYSLpJpWl/7vmlPsf+xtxs2ngfIuNYMxmMGIQNBdK6Zwm6miW5WFRVKrKsiT/TdK3JX0LeDXwG8skViUq9RDwkJnd5dY/RSSMxyQdNrNHnFpPLlOQUmbCsCw0kraMpMVohZyOyZXRVLOqhe9kZm+vM71Si2FmjwI/lPRSt+m1wHeI3mRzo9t2I/C5OgtW+rKURJRKplnLEZ8fTv9bKAiFhTiLES2hy2cYhAwVsjXY5cBwh+GgIXHkCbkkONAaSi11JGnVli5RtR/j14CPudc7PQj8WyJRfVLSTcD3gbc0U8Q5KerwS27OCt06azFU2GxTSjWbuq7TsUpfhUrCMLN7gKy3ab621tJUoaBJNTNxAkx3htF/Q26IrUEgwnB6aNyEOn+4zWWbp9gaXtpI8ePnyRWPoE8Owe9hBSqlBsd6FXSr57sqZeJwvd7R/jh86+K2zvG2eJ3ocyBjqJD9gx0uHD7HZjBurPiZnYfLVp6OGiDRvWZSFfoxiLCMVK/45HO8Hrrbc6ipr2HCwoAwDBiHAaGJQCFHt57gZw88wE8c+HEjv6gCGAxCN+M5e++o8/oXXWygp+ijj9FdYZRdqfTuMmd8Io7IEQ/HkTjG42DigAcyXrzvUa7bCnjB/ifr+y4xbqDicDje24EYh5Urp9XB2pRHC/0YddO9ptQ8P3aeo51uVuG2h3L/I78jBMZjAQNO72wC8M0zRzkQPADAlVc+wemdDU6d3oqsixNTPIHYPM97KIgqsgHb2xuMx0GUTizYqvRFDEl6WOT2hVH3Dxsnp9T6ZHNU8UzO7xhH4jCAQIw1IBzDM9riuZ0NvsaLue/UYV5y8CT/6YV38MD2Yb782D/m9O4mz5zdRxgG7O4OIj8lDCYO9MxE0gmtxK8tGwxCgsDY3Rmyc2YjEsRY81mKnoqij8XunsVYlBLrMbvN9XFg2DjAzBiNAszEM2e3CE0cGO7w4MbzeHT7Qp4bbbAzHkTNrjCImmKTXvQo3dlBiolyuNBsSICFRjjWuWMpYnpY9PURBuT/AAaKQ0Gx9XBNIhsbBLA7ChgNQnZ3Bzwz2OJHpw5y7+ZhRuOAM2c3sTDyR6bOO4khJ8m8Zu/+8SDHyfStoabvva7DUnQ0GpWkj9PnrJcw0iSbWYm6JRMWuP6OQOAis2ZB1FkeinAcMBoFhKEIx67pNAoiQYRTkZX2QWS9dXVeJ7vn9PErrLkwMipt/Gx3XLnjOjpyvocCLIhu6qMgceIyIdVF6bmlADoZcarCCpzvjG11X7jCWUcy7uAJJlNzxiuaNoc68+RcX0QR44WxIGU/dNULmyOImRtvqjMwfoApPleJZlLsp9swmk7TBu6EYFrmXNHkfYfkcUXn5FmLnomirz3f3RBGGTmh2NLTKgpCruNPuP+TRYQDIzShAEI3xmoy5EQuHAzT/pIy0tGz0oInt1fLomso7J8y+iGMmAUFMiFLFKHrRB4DIQS7QmMYbMPwDIy3xM4hwwZgkykEnfOeKFfepNN73u0RkHlcKT0VhfcxqpA1xn+Ri5YnkDgMm8ZmD0laiom1GIPGItiFYBc2Txn7njZ2DorRfmEbEA7cCN3AUFhdHEAiREux1ehju6OEPn6l1VuMvApS5WKmwrDRthxxZByWvJspjKzH8CwMTxvjfeLUC8R4E3bPdxZjwyIxOB9j0pte8D0Mq9Z8mhQqa3uFc7uMF0aNVLUsWeKIT9Gk47m48881pTSCjWeNfU+HPPOTA5590cjtXOQL5DBWZG2q3kb7Lgq8xWiWIv8iva+C1TA3njAeN2UBMIDRltDBgHCT+kUR55euKetqKWK8MCqwTNMpff6eZhR7xDFjNUi2auIn+qLuinDD0AB2Dondg2L3goZ+TUWRrZn1dcb6OSSkO89jKLVUPafCtj19eqm+hNhiWBD1VdjQWZCm79jr0LNdQtyP0fSDSpL+taT7JIWSrk3te6+k45IekPT6Kul1tylVtYc8y4IotSG2HMlDA5A5q+HOiYQROcvhsME7+TkgiBnSM080w73AvwL+OLlR0jXAW4GfAn4C+LKkl5hZ4bPL3bEYVZjHmsycF/0wljrfZm5nuPY/0VVp6soU9nY3lOeKacNimNn9ZvZAxq4bgI+b2baZ/T1wHHhFWXotWwwrjMbkzmC+58CZJPduT1sOYyqO2K9woVqLBZCYXcQs1UfRNGsqCGDeDr5LJd2dWD+Wnth5Aa4A/m9i/SG3rZBONaWSoplLJFlOOMwKpKQZZm5dsYga9y8aTr9DzOF8P25mWdM0RelIXwaen7Hrv5jZ5xYoWi6dEkaSPMuSKZi8Po8ZizG7c/KO7XTUao+f0gDnkCigvqhU1dnOU5wArkysH3HbCumXj0EkmHjJJWt4uLI+R43beM4zUyIadY5V3sYwcBN5lS/NcDvwVkn7JF0FXA38TdlJnbUYVShsehVFq/YM0Zj9USyrE86zMG1cSklvAv47cBnwF5LuMbPXm9l9kj5JNN/yCHhHWUQKutTBl0fFizqXSLw1aJcWhGFmnwU+m7Pvd4HfnSe97luMBUbjVhJJVjSrKF/PQsQR8b7RfWGkmXNISbpJNBGKr/jtYOYfVFopRWOokrviyRDKwsFeOPXRP110zMeo6wJWsCreuW6PPl7qblmMZUfezpN+D3+sXmJEcwX3jG4JI4+qzZp5rr9vKrVH/3TRE2FUZYEIlqd5+tiUqtTzLek33Fj3eyX9uaQtSVdJusuNc/+Eez9ft0g/4zHP4qkNhVZp6RKlwpB0BfAfgWvN7J8AA6Lx7R8APmhmLwaeAm5qsqCenmJzLB2i6lipIbBf0hA4D3gEeA3RO78BbgPeWHvpPL0n6uCzSkuXqPKe7xPA7wE/IBLEj4FvAE+b2cgdljvGXdLNku6WdPf42dP1lLoN+tgw7iphxaVDVGlKXUT0FNRVRI8GHgCur5qBmR0zs2vN7NrBwQMLF9TTX/poMapEpV4H/L2Z/QhA0meA64BDkobOalQa4+45B+mg/1CFKj7GD4BXSjpP0WuJXks0hPcrwJvdMTcCtT5B5VkXqkWkeheVMrO7iJzsvwW+7c45BrwHeLek48AlwK0NltPTZ1b7oNJCVOrgM7P3Ae9LbX6QCrMteM5xrL5HW9tkvXq+y/CRptXQMWtQhfUShq/43aSHP0t3hOEr9dqisH9tqRW8OMYL4JzC6FznXRW6YzE8a4noXuddFXo3r1Rr+BG29dFCuDZvtnNJRyU9J+ket/yPKul5i+FpnhXOdu74OzN72TyJeWF4mqUlH8PM7geIBmcsj29KeRpHYVhpaZCrJH1T0v+R9HNVTvAWw9Mwc/kPha8BWHC280eAF5jZE5L+OfC/Jf2UmT1TVBAvDE+zxJM6V6PwNQCLzHZuZtvAtvv8DUl/B7wEuLvoPC8MH31qnhX2Y0i6DHjSzMaSXkg02/mDZeedmz6Gn/SgVdp4UEnSmyQ9BPwM0Wznd7hdPw98S9I9RKPE/4OZPVmW3rlhMbwAVksL4dq82c7N7NPAp+dNb72F4QWxesxg3L8xIesnDC+G7tHDISHrJQwvim7ihbEivCC6i5/UeQV4QfQAA/M+Rjt4QfQHwzvfreBF0T+8j9EwXhT9xAujIbwgekz35oyqQveF4UXRbwzwkyHUiBfE+uAthseTxg8JqQdvKdYLA/P9GB5PBr7newm8pVhfvI/h8aQw81GphfHWYr3xFmMBvCjWHMPG41UXYm5WKwwvivXHDzv3eHLw4doKeCtxTmGAeYvh8aQw/6CSx5NJH51vWYuhNEk/Ak4Dj7eW6XJcSn/KCs2U9yfN7LJFT5b0RaJyVeFxM7t+0bzqpFVhAEi6u2h+0i7Rp7JC/8rbZc7NKTo9nhK8MDyeDFYhjGPlh3SGPpUV+lfeztK6j+Hx9AHflPJ4MvDC8HgyaE0Ykq6X9ICk45JuaSvfqki6UtJXJH3HvS/6XW77xZK+JOl77v9Fqy5rjKSBe+ni5936VZLuctf4E5I2V13GvtKKMCQNgD8CfhG4BnibpGvayHsORsBvmtk1wCuBd7gy3gLcaWZXA3e69a7wLuD+xPoHgA+a2YuBp4CbVlKqNaAti/EK4LiZPWhmO8DHgRtayrsSZvaImf2t+3yKqMJdQVTO29xhtwFvXEkBU0g6AvwS8BG3LuA1RK/Tgg6VtY+0JYwrgB8m1h9y2zqJpKPAy4G7gMvN7BG361Hg8lWVK8UfAr/F9NWPlwBPm9nIrXf6Gncd73ynkHSQ6J1tv55+F7RFse2Vx7cl/TJw0sy+seqyrCttja49AVyZWD/itnUKSRtEoviYmX3GbX5M0mEze0TSYeDk6ko44TrgVyS9AdgCLgA+BBySNHRWo5PXuC+0ZTG+DlztoiabwFuB21vKuxKujX4rcL+Z/UFi1+3Aje7zjcDn2i5bGjN7r5kdMbOjRNfyr8zsV4GvAG92h3WirH2lFWG4O9g7gTuInNpPmtl9beQ9B9cBbwdeI+ket7wBeD/wC5K+B7zOrXeV9wDvlnScyOe4dcXl6S1+SIjHk4F3vj2eDLwwPJ4MvDA8ngy8MDyeDLwwPJ4MvDA8ngy8MDyeDP4/Hzz5fr+9lMwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.imshow(df_dv.reshape(nz, nx))\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b3eb0aa7", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}