From 7490aa8e933f2403fa23d1f35ac6f7d1c05e95d9 Mon Sep 17 00:00:00 2001 From: Adrian Kummerlaender Date: Mon, 10 Jun 2019 14:06:02 +0200 Subject: Add fixed velocity boundaries to generated LBM kernel Interestingly this increased performance to ~750 MLUPS compared to ~665 MLUPS. --- codegen_lbm.py | 84 ++- lbm_codegen.ipynb | 1872 ++++++++++++----------------------------------------- 2 files changed, 467 insertions(+), 1489 deletions(-) diff --git a/codegen_lbm.py b/codegen_lbm.py index cd93649..8e7b0a7 100644 --- a/codegen_lbm.py +++ b/codegen_lbm.py @@ -51,50 +51,48 @@ __kernel void collide_and_stream(__global __write_only float* f_a, const float f_curr_7 = f_i(f_b, cell.x , cell.y+1, 0,-1); const float f_curr_8 = f_i(f_b, cell.x-1, cell.y+1, 1,-1); + const float ux0 = f_curr_3 + f_curr_6; + const float ux1 = f_curr_1 + f_curr_2; + const float ux2 = 1.0/(f_curr_0 + f_curr_4 + f_curr_5 + f_curr_7 + f_curr_8 + ux0 + ux1); + const float ux3 = f_curr_0 - f_curr_8; + + float u_x = -ux2*(-f_curr_2 - f_curr_5 + ux0 + ux3); + float u_y = ux2*(-f_curr_6 - f_curr_7 + ux1 + ux3); + + if ( m == 2 ) { + u_x = 0.0; + u_y = 0.0; + } + const float x0 = f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f_curr_7 + f_curr_8; - const float x1 = 2*f_curr_0; - const float x2 = 2*f_curr_8; - const float x3 = -f_curr_3 + f_curr_5; - const float x4 = pow(x0, -2); - const float x5 = 9*x4; - const float x6 = f_curr_0 - f_curr_8; - const float x7 = f_curr_1 - f_curr_7; - const float x8 = f_curr_2 - f_curr_6; - const float x9 = x6 + x7 + x8; - const float x10 = 6/x0; - const float x11 = x10*x9; - const float x12 = f_curr_3 - f_curr_5; - const float x13 = -f_curr_2 + f_curr_6 + x12 + x6; - const float x14 = pow(x13, 2); - const float x15 = 3*x4; - const float x16 = -x14*x15 + 2; - const float x17 = x11 + x16; - const float x18 = pow(x9, 2); - const float x19 = x15*x18; - const float x20 = -x19; - const float x21 = x10*x13; - const float x22 = x20 + x21; - const float x23 = 1.0/$tau; - const float x24 = (1.0/72.0)*x23; - const float x25 = 6*x4; - const float x26 = x18*x25; - const float x27 = (1.0/18.0)*x23; - const float x28 = x5*pow(2*f_curr_2 - 2*f_curr_6 + x3 + x7, 2); - const float x29 = x20 - x21; - const float x30 = x14*x25 + 2; - const float x31 = -f_curr_0 + f_curr_8 + x3 + x8; - const float x32 = x15*pow(x31, 2) - 2; - const float x33 = x19 + x32; - - f_a[0*$nCells + gid] = f_curr_0 - x24*(72*f_curr_0 - x0*(x17 + x22 + x5*pow(-f_curr_1 + f_curr_7 - x1 + x2 + x3, 2))); - f_a[1*$nCells + gid] = f_curr_1 - x27*(18*f_curr_1 - x0*(x17 + x26)); - f_a[2*$nCells + gid] = f_curr_2 - x24*(72*f_curr_2 - x0*(x17 + x28 + x29)); - f_a[3*$nCells + gid] = f_curr_3 - x27*(18*f_curr_3 - x0*(x22 + x30)); - f_a[4*$nCells + gid] = f_curr_4 - 1.0/9.0*x23*(9*f_curr_4 + 2*x0*x33); - f_a[5*$nCells + gid] = f_curr_5 - x27*(18*f_curr_5 - x0*(x29 + x30)); - f_a[6*$nCells + gid] = f_curr_6 - x24*(72*f_curr_6 + x0*(x10*x31 + x11 - x28 + x33)); - f_a[7*$nCells + gid] = f_curr_7 - x27*(18*f_curr_7 + x0*(x11 - x26 + x32)); - f_a[8*$nCells + gid] = f_curr_8 - x24*(72*f_curr_8 - x0*(-x11 + x16 + x29 + x5*pow(x1 + x12 - x2 + x7, 2))); + const float x1 = 6*u_y; + const float x2 = 6*u_x; + const float x3 = pow(u_y, 2); + const float x4 = 3*x3; + const float x5 = pow(u_x, 2); + const float x6 = 3*x5; + const float x7 = x6 - 2; + const float x8 = x4 + x7; + const float x9 = x2 + x8; + const float x10 = 1.0/$tau; + const float x11 = (1.0/72.0)*x10; + const float x12 = 6*x3; + const float x13 = x1 - x6 + 2; + const float x14 = (1.0/18.0)*x10; + const float x15 = -x4; + const float x16 = 9*pow(u_x + u_y, 2); + const float x17 = -x2; + const float x18 = x15 + 6*x5 + 2; + + f_a[0*$nCells + gid] = f_curr_0 - x11*(72*f_curr_0 + x0*(-x1 + x9 - 9*pow(-u_x + u_y, 2))); + f_a[1*$nCells + gid] = f_curr_1 - x14*(18*f_curr_1 - x0*(x12 + x13)); + f_a[2*$nCells + gid] = f_curr_2 - x11*(72*f_curr_2 - x0*(x13 + x15 + x16 + x2)); + f_a[3*$nCells + gid] = f_curr_3 - x14*(18*f_curr_3 - x0*(x17 + x18)); + f_a[4*$nCells + gid] = f_curr_4 - 1.0/9.0*x10*(9*f_curr_4 + 2*x0*x8); + f_a[5*$nCells + gid] = f_curr_5 - x14*(18*f_curr_5 - x0*(x18 + x2)); + f_a[6*$nCells + gid] = f_curr_6 - x11*(72*f_curr_6 + x0*(x1 - x16 + x9)); + f_a[7*$nCells + gid] = f_curr_7 - x14*(18*f_curr_7 + x0*(x1 - x12 + x7)); + f_a[8*$nCells + gid] = f_curr_8 - x11*(72*f_curr_8 + x0*(x1 + x17 + x8 - 9*pow(u_x - u_y, 2))); moments[gid] = x0; }""" diff --git a/lbm_codegen.ipynb b/lbm_codegen.ipynb index 7d13f70..34e9e6b 100644 --- a/lbm_codegen.ipynb +++ b/lbm_codegen.ipynb @@ -14,28 +14,30 @@ }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "c = [Matrix((x,)) for x in [(-1, 1), ( 0, 1), ( 1, 1), (-1, 0), ( 0, 0), ( 1, 0), (-1,-1), ( 0, -1), ( 1, -1)]]" + "c = [Matrix(x) for x in [(-1, 1), ( 0, 1), ( 1, 1), (-1, 0), ( 0, 0), ( 1, 0), (-1,-1), ( 0, -1), ( 1, -1)]]" ] }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/latex": [ - "$$\\left [ \\left[\\begin{matrix}-1 & 1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0 & 1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1 & 1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1 & 0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0 & 0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1 & 0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1 & -1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0 & -1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1 & -1\\end{matrix}\\right]\\right ]$$" + "$$\\left [ \\left[\\begin{matrix}-1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right], \\quad \\left[\\begin{matrix}-1\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}0\\\\-1\\end{matrix}\\right], \\quad \\left[\\begin{matrix}1\\\\-1\\end{matrix}\\right]\\right ]$$" ], "text/plain": [ - "[[-1 1], [0 1], [1 1], [-1 0], [0 0], [1 0], [-1 -1], [0 -1], [1 -1]]" + "⎡⎡-1⎤ ⎡0⎤ ⎡1⎤ ⎡-1⎤ ⎡0⎤ ⎡1⎤ ⎡-1⎤ ⎡0 ⎤ ⎡1 ⎤⎤\n", + "⎢⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥, ⎢ ⎥⎥\n", + "⎣⎣1 ⎦ ⎣1⎦ ⎣1⎦ ⎣0 ⎦ ⎣0⎦ ⎣0⎦ ⎣-1⎦ ⎣-1⎦ ⎣-1⎦⎦" ] }, - "execution_count": 31, + "execution_count": 3, "metadata": {}, "output_type": "execute_result" } @@ -46,7 +48,7 @@ }, { "cell_type": "code", - "execution_count": 32, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -55,7 +57,7 @@ }, { "cell_type": "code", - "execution_count": 33, + "execution_count": 5, "metadata": {}, "outputs": [ { @@ -68,7 +70,7 @@ "[1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36]" ] }, - "execution_count": 33, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } @@ -79,7 +81,7 @@ }, { "cell_type": "code", - "execution_count": 34, + "execution_count": 6, "metadata": {}, "outputs": [ { @@ -92,7 +94,7 @@ "1" ] }, - "execution_count": 34, + "execution_count": 6, "metadata": {}, "output_type": "execute_result" } @@ -103,7 +105,7 @@ }, { "cell_type": "code", - "execution_count": 35, + "execution_count": 7, "metadata": {}, "outputs": [ { @@ -118,7 +120,7 @@ "3 " ] }, - "execution_count": 35, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -130,7 +132,7 @@ }, { "cell_type": "code", - "execution_count": 36, + "execution_count": 22, "metadata": {}, "outputs": [], "source": [ @@ -139,7 +141,7 @@ }, { "cell_type": "code", - "execution_count": 37, + "execution_count": 23, "metadata": {}, "outputs": [ { @@ -149,7 +151,7 @@ " f_next_6, f_next_7, f_next_8], dtype=object)" ] }, - "execution_count": 37, + "execution_count": 23, "metadata": {}, "output_type": "execute_result" } @@ -161,7 +163,7 @@ }, { "cell_type": "code", - "execution_count": 38, + "execution_count": 24, "metadata": {}, "outputs": [ { @@ -171,7 +173,7 @@ " f_curr_6, f_curr_7, f_curr_8], dtype=object)" ] }, - "execution_count": 38, + "execution_count": 24, "metadata": {}, "output_type": "execute_result" } @@ -183,7 +185,7 @@ }, { "cell_type": "code", - "execution_count": 39, + "execution_count": 25, "metadata": {}, "outputs": [ { @@ -197,7 +199,7 @@ "_curr_7 + f_curr_8" ] }, - "execution_count": 39, + "execution_count": 25, "metadata": {}, "output_type": "execute_result" } @@ -209,1617 +211,595 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 55, "metadata": {}, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAAlCAYAAADiFfbCAAAABHNCSVQICAgIfAhkiAAACCJJREFUeJztnVuM1FQYx3/LIoKg3BRCZA1e42VXN4KioogEifeAGHnwZY0YEjEYHwzGROXFqC9qDPEaSdUoRAWioJiIKBqVi2YBkUQQXQRUBEHAIAq4PnynbLd0mE57Oj2d/X5Js3O6nfb/9fs6p+2Z6R8URVEURak5pgA/AoeA53PWkpSix6D686foMRRdv6IAcC5SxJOAIUDvfOUkougxqP78KXoMRdevKEd4APgmbxEhPGBmBcu7GEMluKjfI34OXNRfKS7G4NG1cqAobADaA9O8fOUcwSP+wehqDHFxVb9HvBy4qr8SXI3Bo+vkoCbonreAGuBK4HPgNeAVYH++chJR9BhUf/4UPYai668ZuuUtoAbYC5wBfAH8Ztq2uQH4HtgI3JPB+rOMoQH4FFgPrAFutbhun2rk4EtE/zrgEcvrrob+NmAtsBr4JIP1Zx3DMGApUkffASdbXn81cjANqZ91SMdXn8E2lBplJp0v5aOmMWbZEabdz9K2w4V6PNIZNQAnIAflkIj3PQT8FZgOAv+G5l1VYps2YwjrPxW4yLweBGxB4ohiJvH3e5Csc1APnBR4/RXQHPG+pDmohv42oE+Z980k2f6HbGuoHlgGjDbtvshxEYWrORgMbAJ6AnXAQuAmS9uqKfT2XTSzgLlllvnZ/G0GNgN/Bv7XYNZxGtADKb6fgFZgPLADGI4cQJOQAl2PHPCzzfJ++wPzeotZ9wLz/5dDel4A3gq0nwS2Ac8G5m0rEUvcGOYn0D8beNGs83dgN3KW6++/IJXsd5v6qSCGHmaKImkOsq6h2SX0hkm6/yuJIUkNvYl0Lp+Z9e45hj5XczAf+bztCRwGegHbjxFHl0U7pWh2mikOzcgtEZ8eSEcyHblN0hc4gOzr/kghAzQhl/EAjcgtrpGm/WCgfRtwSmD9W5GrjzC7zOSzz7R/sBhDEv1BRgDH0dHBhqlkv9vUD/FiWAGcDzwX2p5P0hxkXUN+exnwH/AM8EaEjqT7v5IYkuifgOzLd5EOYiGlb6G6nIP9SKd+EJgDrCqjqUuiY0rpCRfzRGAlHfft9wD/AOfQ+cC4ECnmE5ED4CkzP9yui9hmuw3hAeLEcCbJ9PsMRAaR78It/RA/hpHICUEz8gFki6xrCGAUclZ/C3Jm32RRP2RbQ92Rq43pwKVIHBNz0J8mB/2BG5GxsaFI/YyxG0JtoJ1SOuqQogwWcxNSzGEa6XxmPsK0G5FB9PbAcsH2NuQ2gs9Q4Je0wgPEjSGpfpD7/wuAx83/bJJWv/+/cjH47EUG3K9LpbqDatQQdNTMr8gVwPC0wgNkXUNbkd8PbUauMhYRPaaXlGrkYBwyprQL+Bt4H+lglRB6+y4d7XQMgPtsp+MsuhtyhvQHMAApRoDLkA+FjcCdyLeifBpD7ZXABUjHtBM5gxsXQ1uL5RiS6q9DfiuyFHg9pqZKSKsfjtYcbvdFznp3IGMC44GnY2hrsag/TQ31NuvZh3zZYSydx13SknUNrUKutAciH+pXIyc5cWixqD9NDrYAlyP1cxC5SnopXghdC71Sso+H3KZYh5zdnWfmL0YGPucgX/HeiDzSpInOxRtuHwLuBz4GvkWex2XzSikKj6NjSKp/FDAZGRdYbSbbt47CeMTXH6U53O4PfGjmfY2MzSzKTL39GhqM/AZnDbAcuY2a9XiGh70aOgzMQG6lrUVOzmx2qlF42M3BcvPeVjN/E/BeZuoVRVEURVEURVEURVEURVEURVEUpYbxfwNj+3cjiqIoiqIoiqIotUXR7YBVf/4UPQbVnz9Fj6Ho+sGRGIpuB6z686foMaj+/Cl6DEXXDw7F4JodsEfxLb09im3H7KE5yBuPrqMfih9D0fVDjjEEHzO0ATjbvG5HHrU+6ah3uIvqz5+ix6D686foMRRdP+QcQ7BTKrodsOrPn6LHoPrzp+gxFF0/5BxD8Nl3WdsBF9nSG+Q5VbuBdyyv1ydr/bVgS56lJTlUxxIb5Lhbhf1aUlv18gyj2LbqTXQ8Q3I18nDYCZa3kXUMsW3hs7QzLrqldz1wDXAz5T9IXLVjjmtL7nIO4liSg7s58NvTEIfXUrXksv42ytuqu1xDRbdVD7b7IA+njfoSgqs5qMgWfgpScD4NiNNjK3JGcbqZ30qHE+pwYJ55vRCxHl4BTA21HzXr8nkMuDtCwwDgrMA0D7ExDs7rFVN/qRiS6J9qlhlD+U4paQxx9ZMyBp+1iIunLf2VxJBWfy9kELZUp5R1DkrpjxPDIOSJ72MpXUsu11Ab5TslV2voPmBJGe1pY6jmcXwHpe3rXc3BDMQXqx/iQr0EuCS4oeCYUtB5US29S1t6l8NlO2afY9mSu56DcpbkaWKohq36q8DDZt2lcLmG4tiqu1pDtWKr7nM78q06m/oriSETW/jgmFJQiFp6R9thZ0k1LLHBXVvyvC3JIXtb9dHIPrftvuujtuq1b6vu0xcxDVxsUbtPljkoawvvd0phO2C19C5th50F1bLEdtmWPE9LcqiOrfoVwLXIrZG5wPXIQK8N1Fa9a9iq+0xAjCcPpFbdmaxzUNYW3r99F7YDVkvvzu00tMRYphp2zEltyVtiLpd1DpJakoPdHKSxVX/CTCBnh/ciV6zlaImxjMu26nH0VxJDV7ZV95kMzIqhyacl5nJZ56CsLXwpO3QPtfQOF8IS4G3zvq3Ijs0SD7s5cN2W3DVLcrBvq15tPNRWvdZs1UE6hIuBj7ISHcLDXg7UFl5RFEVRFEVRFEVRFEVRFEVRFEVJzf8pYN+WylHVPgAAAABJRU5ErkJggg==\n", - "text/latex": [ - "$$\\frac{- f_{curr 0} + f_{curr 2} - f_{curr 3} + f_{curr 5} - f_{curr 6} + f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}}$$" - ], - "text/plain": [ - " -f_curr_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_curr\n", - "──────────────────────────────────────────────────────────────────────────────\n", - "f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f\n", - "\n", - "_8 \n", - "──────────────────\n", - "_curr_7 + f_curr_8" - ] - }, - "execution_count": 40, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], + "source": [ + "#u_x = sum([ (c_i*f_curr[i])[0] for i, c_i in enumerate(c) ]) / rho\n", + "#u_x" + ] + }, + { + "cell_type": "code", + "execution_count": 56, + "metadata": {}, + "outputs": [], "source": [ - "u_x = sum([ (c_i*f_curr[i])[0] for i, c_i in enumerate(c) ]) / rho\n", - "u_x" + "#u_y = sum([ (c_i*f_curr[i])[1] for i, c_i in enumerate(c) ]) / rho\n", + "#u_y" ] }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 26, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAAlCAYAAADiFfbCAAAABHNCSVQICAgIfAhkiAAAB3pJREFUeJztnWuIFWUYx3+blpqlrZomaXSny25JWlaWbWLRnbWiPvRlI0PIMPoQSVHtl6i+VEQUBMLJCKWyMC2DSovIUgsvWZBmrK3dcDPSyC7W9uF5hzM7zpwzl/edi/v8YOC8s+fM+T/7PDPvO/POmT8oiqIoilJp5gLfAgeA5wvWkgbVXzxVj0H1K0pJOAMp5JuAicDIYuUkRvUXT9VjUP2KUiLuAz4vWoSPGtCd4P1l0w/JYiij/qRUPQbVryglYRvQ71uWFSsHSHZAL6N+iB9DWfUnoeoxqH5FKRHjga+BB4HjgFHFygGSdUpl1A/xYyir/iRUPQbVrzjlsKIFVIy9wMnAx8BPpm2Ta5AdZjtwl+Vtg3v9bwK/Aq9Z3q6Ha/2TgQ+Ar4DNwI2Wtw/uYzgRWI3E8CUwzvL2XetvBzb5lv1Ap8Xtu9YPMB/YapZFwBAH36EoAExDTvmPsbCtYKEOQzqjycCRyEFlYsjnHgB+9y3/AH8H1l0a8Z0u9Q8BLgeup3mnlDYGm/rh4BiOB841r8cDvUgugnQz8BJQ2NIR8Z2uc/AhMNO0RyN1FUY36WJwnQN/+yigj/AbEbopp/4JwA5gONACrACus/RdinIQc4EeX3sysBzYiIxKTzLrNwLHmtdTqV+3XgE8AawD5gXaj5hteTwK3BmiYQxwqm9ZBjwTWDcipv6oGNLon2fe00HzTiltDHH1kzEGjy3ACSE6xiF3cDVawjqzJDGk0X8P8F7E99qKIav+ZjH4c3AbsLQg/TSIoZH++4GdSKd3OJKP8yNiUEIYWrSAijEFuaQAcATwNrAAWIOMSv9E/qetwG7zvnbkNB6gDbk8NN20F/raN1PfAQB2ISP3IHvM4rHPtL9JqL9RDGn0JyFtDHH1p82Bn2nIQaU3REefWdLgMgedyP9yOdKZrgAejtCRNoas+pvF4OcWZL4xDNf6s9TQH8B3yBWAJcCGFDoHLTqnlAx/Qc8B1iOFDPAb8BdwOgMPrucgxXw0sgM8adYH2y0h39dvS7ghuEOGxXAK6fTnQRz9WXLgMRZYDNxBtXIwFDlTXQBcgIzu55RIf5jmqByMBi4CVtkSbnBdQ63Atcjc3iSkA+uwG8KhjXZK8WlBCtMr6HakmIO0MXBUOM2024C11A9ywfb3yGUEj0nADzaEG4L6ITyGtPpdE1c/ZIthGPAG8Jj5m01c52AX8vubncgofSVyELZFVv1hmqPqqBN4BzlrsUUeNTQbmVPag9yk8RYyQFBiopfv4tPPwNtHf0YKEqRzbwV+QeZL9pv1FyKj1e3A7cgchUdboL0eOBvpmPqQEdzsGLq6UuqPiiGt/ix0xXhPXP1ZctCCXC5aDbwUV3wCXOdgA3KWNxY5KF6GdLC2yKo/THNUHd0KPJtd8gDyqKFe5AxvODIw6ABesBXAYEDPlNJTQy5TbEVGp2ea9auQu22WILd4b0ceadLOwOINtg8A9wLvA18gz+OyeaYURo2DY0irH2RS91XzuV3IzumSGnZzMAM5GHZSvyW53Zl6oYa9HPyLTLSvMev7gFecqk+mP0xzWB2NAc4D3nUl2kcNuzX0qfnsRrN+B/JTCUVRFEVRFEVRFEVRFEVRFEVRFEVRLOH9Niav23oVRVEURVEURVGqSdVtglV/8VQ9BtVfPFWPoer6oSQxVN0mWPUXT9VjUP3FU/UYqq4fShRD2WyCa6jVd9HU0BwUTY3Box+qH0PV9UOBMfgfM7QNOM287gdeR3rJqqD6i6fqMaj+4ql6DFXXDwXH4O+ULgE+Qp6OvAh5/HqVUP3FU/UYVH/xVD2GquuHgmPwP/tOrb4bo1bfzXEdw1pE+1aifYKykIdVNsh+twH7tZSH/h7kmW6bqNs92ETt4hvj2i4e3McQ2y5erb7rDGar7zLnYJTv9SdE2zKUNQdeez7iqBpVS2XW34PYlDeizDUU1y6+zDnwaGQXX9YcJLKLV6tvtfrOoj9JDFn1j0AmYaM6Jdc5yGL1PR55EvwsomupzDXUQ/NOqaw1lMQuvsw58GhkF1/WHDS1i/fPKanVt1p9Z9GfJIYs+tcBZwHPBb7LRgx5WH2/CDxkth1FmWtoIXK28R/wNPCyRf1JYnBtF1/mHHg0sosvaw6giV28f05Jrb7V6jsreeRgOjKgmULdnM0Wrq2+ZyL/c9uOth551NAMZFR8A3KJyLbf1GCwi7exH7uyiwe3OWhqF+91Smr1rVbfWckzB3uRyeqrMquuk4fV98XAFcilkaXA1chErw3yqiFvv/0RGT1PzSrcx2Cxi/e0ZdkPXNjFg/scNLWL9y7fqdW3Wn1H0RXzfa5zMBoZce1GJkmvBJ6Kqa0rB/1hmoPtx80CMjq8GzljbUZXjPfkUUMjzXb2IfNKs4jnbBtHf5IYirCL77KoP0sOPJLaxXfFfJ/rHDS1i4+yQ6+hVt9q9Z2dGvZy0IqMDLcAnyHzGiudqndj9Z0nNezW0ATk9yubEdvvxQTmAxxQ49Czi4dsx6I87eLBbg7ULl5RFEVRFEVRFEVRFEVRFEVRFEXJzP/mgKdK6IgqEQAAAABJRU5ErkJggg==\n", "text/latex": [ - "$$\\frac{f_{curr 0} + f_{curr 1} + f_{curr 2} - f_{curr 6} - f_{curr 7} - f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}}$$" + "$$\\left[\\begin{matrix}u_{x}\\\\u_{y}\\end{matrix}\\right]$$" ], "text/plain": [ - " f_curr_0 + f_curr_1 + f_curr_2 - f_curr_6 - f_curr_7 - f_curr_\n", - "──────────────────────────────────────────────────────────────────────────────\n", - "f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + f\n", - "\n", - "8 \n", - "──────────────────\n", - "_curr_7 + f_curr_8" + "⎡uₓ ⎤\n", + "⎢ ⎥\n", + "⎣u_y⎦" ] }, - "execution_count": 41, + "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "u_y = sum([ (c_i*f_curr[i])[1] for i, c_i in enumerate(c) ]) / rho\n", - "u_y" + "u = Matrix([u_x, u_y])\n", + "u" ] }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 27, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAKwYAAAAfCAYAAABhqLSIAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJzt3Xn4LFld5/k3RQHFAAKFyDJtg7aySLGIVxQc6AQRXKAFlXZ0BknbYnSARppHuwV1uOOI0NIOOzNqd3uHwRUQbFrZdCwElUUBGwQaWX6IUCCICAoUFNz+IyLn5s2bS2REnDjfE+f9ep56bt1cfr/Ic058M+KT34wLkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ0kzdDXgH8DLgxcByz2MvmmKDEnkc8AbgE8BHgJcAl2XdIkXkOpHych88p+T3XEmSJElSvUo9nzWTGMbxi6/UfVOSJEmS+ij5HMhzbHXhOpGGcR9qlPx+KUmSJEmqV47zWbOE+XJuzzErkiRJklSTUs+BPI8dxvFrlLr+JUmSJEl1s19EY3JuzzErkiRJklSTUs+BPI9VV64VKS/3wUap77eSJEmSpLqVej5rHjGM41eGUvdPSZIkSeqj1HMgz7HVlWtF6s/9p1Hqe6UkSZIkaTo/BLwYeAXwDuAO+x68AM4Cpw780NsA3zfCxuXycuD7acKEOwIvAj4EXJpzoxSO60TKy33wnK8G7p17IyRJkiRJOkLJGbKZxDCOX3xmTZIkSZJqUXI+AZ5jqxvXiTSM+1DDzFCSJEmSVJpc+a9Zwnw5t+eYFUmSJEmqRcn9ZZ7HDuP4NcwAJEmSJEmlsV9EY3NuzzErkiRJklQL+0VUA9eKlJf7YMPMUZIkSZJUGvPjejl+ZTBvkiRJklQLMwrVwLUi9ef+0zAvlCRJkiR1tQDOAqeGPugGwH8GLhlpwyK4PvB54IG5N0ShuU7q8lzgr4Hr5d6QGfsamvecH+j4+Nr3wWfQfIDaxaH1+0XtzzsBPkszD/9m4PZpGOckHuckJuclHuckHuckJuclHuckHuekm2PPpXOZW4ZceyYxlOMX0zFZE+zPm6zh6TnGaTm+6TnG6TnGaTm+6TnGaTm+6dU4xiVkFHPLJ8BzbHXjOqmLPW6H2aPWnf1p8+acxOS8xOOcxOOcxOS8xOOcxOOcxOOcdFNC9gux8t+as4S5q31uzYrK4fim5xin5fim5xin5xin5fim5xinVeP4mi/kUft57FA1j9+YGUAfNdbJ6JyTeJyTmJyXeJyTeJyTmJyXeJyTeGqcE/Oc49WcJcxd7XNrv0g5HN/0HOO0HN/0HOP0HOO0HN/0HOO0ahxf84U8aj+PVXeulbp4PaL0vJ5Rd2aO8+acxOS8xOOcxOOcxOS8xOOcxOOcxOS8HGZ+nEfNecQYHL+4zJvK4fim5xin5fim5xin5fim5xin5xinVeP4mlHk4Tm2unKt1MUet/3sT+vOvHD+nJd4nJN4nJOYnJd4nJN4nJN4nJOYnJduDp1LL9r7T+37IV0e9AvAfY/evO7+LfB7wPuBTwMfA94EPAG4SaLfeQua1/0NO+6/GU3w8Iz27zcBLgdeBLyr3c6/A15DMwEXJdrOQ3KMXS4R1wm4ViK6BvAvgNcCnwQ+RfNaHw1cc8dzTgFfAB47xQZW7kXAlTQB7yFz2ge/C3gm8GrgEzSv63kHnnNT4FUc3uYu6/cl7e/8beCngdPA7Q9tdEab8wrD6tg9gRfSrL2r2j9fAXzrqFt9nNLmBC6cF+ckvzH3lYhzAuXNi/UrJutXPNaveKxfMVm/4slVv445l84lZYYcMRcsJY/IlQl2yXRSqSkHPTZv6po1weG8aS41vE+WPJXSxtjxTc8xTm9zjCOPL5Q3xq7h9FzDabmG03OM0xprfKNnFHPLJ8CMokQRsyw4f61EXSdQ11qxxy2tMXvUSqm19qcdz8/i4rGXIB57CeJxTmKyfsXjvhKPcxKTx8Tx5NpXome/sDv/jZrHpVJTfhZ1bs2K5pkV+Vlleo5xWo5veo5xevY+peUaTs81nJZrOL0xxrjkfGEMJZzHRj2HhZjjZwZwXK/SQ9vffZZmnPYprU6W9j60Muc5gfKOv1a6zssc5gTKmJc57yvOSUzWr3jcV+JxTmIqsX7VNidQV54TMUtIyX6R/HNrVmS/SIT3uNLGFxzj1Bzf9Bzj9Eo71yhtjF3D6bmG03INp1d7vjCGEs5jo57DQj2Zh9etGqaWdQL93+e8HtF0vJ6RmWMXfnc7Hq8HEo/XnojJ+hWP9Sse61dM1q94rF/xeD2j3cyPY+YREHf8UqkpBzVvOo6fmabnGKfl+KbnGKfl+KbnGKdnj1taruH0xljDZhQxz7HNKGKxx22YWtYJ2OOWmv1p5oVd+FlcPPYSxOOcxGT9isd9JR7nJCb70+JxX4kn55zsO5de0IzjqX2/7NCDvgZ464ENHuqzNGHLfwSeTHOi/IZ2uz4AfGmC3/nrNBOyK9D5X9rff+/27z/U/v2DwC8DT2q39+Pt7S+gCY6mlmPscom4TsC1EtFzaV7Xh4F/Dzwd+HP2j/8raObouhNtY83uRjMXj+/w2Dntg29uf/8ngbfTLQgG+HngfzrwmEPr93bt73tZpy2NYXNeoX8d+4n2MR8Bfgn4GZoPp98A/GyCbe+ixDmBC+fFOclvrH0l4pxAmfNi/YrJ+hWP9Sse61dM1q94ctWvY86lc0idIUfMBUvJI3Jlgl0ynVRqykH75E1dsibYnzfNqYb3yZKnUOIYO77pOcbpbY5x1PGFMsfYNZyeazgt13B6jnFaY41v5IxijvkEmFGUKGKWBeevlajrBOpaK/a4pTVmj1optdb+tOP5WVw89hLEYy9BPM5JTNaveNxX4nFOYvKYOJ5c+0rk7Bf2579R87hUasrPos6tWdE8syI/q0zPMU7L8U3PMU7P3qe0XMPpuYbTcg2nN8YYl5wvjKGE89io57AQc/zMALr3Kn1p+7hPtr//8j2PLbFOlvQ+tDL3OYGyjr9Wus7LXOYE4s/L3PcV5yQm61c87ivxOCcxlVa/apwTqCvPiZglpGS/SP65NSuyXyT3e1yJ4wuOcWqOb3qOcXolnWuUOMau4fRcw2m5htOrPV8YQwnnsVHPYaGezMPrVg1TyzqB/u9zXo9oOl7PyMyxC7+7HY/XA4nHa0/EZP2Kx/oVj/UrJutXPNaveLye0Xbmx3HzCIg7fqnUlIOaNx3Hz0zTc4zTcnzTc4zTcnzTc4zTs8ctLddwemOsYTOKmOfYZhSx2OM2TC3rBOxxS83+NPPCLvwsLh57CeJxTmKyfsXjvhKPcxKT/WnxuK/Ek3NO9p1LL9r7Tu3b+EMPehHwr/f9gBFcsuP2J9Js23NG/n1PAa4EvmLPY14KfJRzwcN9gAcCF2087ubAX9Js53eOu5mdTD12OUVcJzDdWlm2z130eC6UtVaW9H+tD2qf+x7gi9duvxZNPTvb/vx1twG+QFNs52bJsHWTytuB97H/g7S57YP3Br6SJoRe0D0Ivgfwbpo1vE2X9fuI9vc9vOO2RrA5r9Cvjj2kve+VwA223L9rXPdZMny/KnFO4MJ5iTInMHxe5jIncPy8RJ0TKHNerF8xWb/isX7FY/2KyfoVT8761eVcOpfUGXLEXHDK/HhJ//0tRybYNdPZZYk5aFd98qZDWRMczpvmUsP7ZMldLKnz2GWq8QWPQ1zD6ayPceQ1DGWOsWs4PddwWr7XpWedSGvM8Y2aUcwtnwAzilJFzLLg/LUSdZ1AWWtliT1uELc/DcbrUbM/bZ79aeBncRHZSxCPvQTxRJ4TsH5Zv+KIvK8scU5WnJP8PCaOJ+e+EjX7hf35b9Q8LpWS8rOhos6tWdE8syI/D07Pz4PTcg2n5xpOz96ntFzD6bmG0/K9Lr2xxrjUfGEMJZzHRr1+GsQcv1J683JmALS/93fbn/WU9vdfvufxJdbJkt6HoI45gbKOv+C4eZnDnEDs8w6oY1+xfsVk/YrH+hWP9SumkupXjXMC9eU5EbOEQ5b0X9/2i+SfW/tF7BfJ/R5X4vhCWcfbJY6xazg913B6JZ1rlDjGruH0XMNp+V6XXu35whhKOI+1XyQ/r1s1TEnrZEn/95O+NTjq9YjmrLTrGUG5/WlzOcbyu9v5eT2QeLz2REzWr3isX/FYv2KyfsVj/YrH6xltZ34cN4+AuOOXSi05KJg3HcvPTNOztyIt13B6ruG0XMPpuYbTs8ctLddwemOtYTOKC+U+xzajiMUet2FKWidL7HGLzv4088JD/CwuHnsJ4ok8J2AOZP2KI/K+ssQ5WXFO8rM/LR73lXhyz8muc+lF+/NO7XjewQf9I+DzwJft+wEJ3ZlzA7Lu8e3tD97ynFu19/3mjp/5c8CHga/a83tvCFwF/FLH7VxtzzM7Pn4Ku8YOho1fRClea5d1AtOulSXDC+w2EdfKkv6v9bntcx+55b7L2vv+dOP2J7e3f2OP39dViWMJ6bb7Ce1z77/j/rnvgwu6B8EAfwV894779q3f72zv2/bf7bpv7uSOnddddewimg90/gG46WhbN2wtlDoncNy8TD0n0H9eapkT2D4vEecEyp0X61dM1q94rF/xWL9isn7Fk7t+HTqXziVnhlxLfrxk/FwwVSbYNdPZZ8m0OWjOvHxJnrxpX9YEu/OmudXwPllyF0vqO3aZcnzB45AV1/C4Nsc44hqGcsfYNZyeazgt3+vSs06kNfb4RswoIuYTYEbRRcS+pVRK6XGLuE5gftlOih43+9MudOg9K1qtXZInL4T6+tPAz+Iiyv1Z3CFLPJc4xF6C9KLPCVi/urB+pRd9X1ninBzinEzDY+J4cu8rEbNf6J//5szjcjBrbUTKWte3x6wodg328+D0/Dw4Lddweq7h9Ox9Sss1nJ5rOC3f69Ibc4znli+MoZTz2IjXTwO/P3isBdNkAOt+mOYC+PcCTrfPuXzL40qtk6W8D62b+5xAOcdf67rMy5zmBGKed6yb+75i/YrJ+hWP9Sse61dMpdWvGucEzHNWIveLLBnnvWCd/SKNSDnb+vbYLxK7npZy3FHq+EI5x9uljrFrOD3XcHqlnGuUOsau4fRcw2n5Xpee+UJapZzH1tQvElFN161KIWI2tqT/e3bfGpz639wrdb/zekbnWzJO7rjAzPEQv7sdj9cDiSf39+kPWVLfnID1KyLrVzzWr5isX/FYv+LJXb/Mjy9US368ZJw8YlOu7xvmMrccdNMC86Z9/Mw0PXsr0nINp+caTss1nJ5rOD173NJyDac35ho2o7hQznNsM4py2OM2zNyynRQ9bqXuM/annW/JOLV2gXnhIX4WF0/uz+IOWeKcHGIvwTSsX/FE31eWOCeHOCfTsD8tHveVeCLMya5z6UV7+6nVDRcf+YO/G/g74L1HPm8sD2z//C8bt9+1/XNbELN6sW/act8zaF7TvYG37fm93wZcm+7hxufaP6/u+Pgp7Bo76D9+UY39WruuE3CtQLy1cvP2z/dsuW91212BGwEfb/9+X5oPq16bcLtKHEtIt91/2P75TcDLN+6rbR/s4o3A9wC/vuW+fev3SuB/Bx4B3Bh4Ynv7WeAvxt/M0Rw7r7vq2D1oPoR+AfC37c+9DPgM8Hrgjwdv6fFKnRM4bl6ck2mMsa9EnBMod16sXzFZv+KxfsVj/YrJ+hVP7vq171w6p5wZsvlxfykywWMynRzGXi8l25c1we68aW41vE+WnFqpY+z4pucYp7c5xhHHF8odY9dweq7htFzD6TnGaY09vhEzioj5BJhRdDG3vqV9Sulxi7hOYH7ZTooet1LHIuV2j9GjNoda20Vt/WngZ3ER5f4sLqVS58Vegnick5isX/G4r8TjnMTkMXE8ufeViNkv9M9/c+ZxOZi1NiJlrWBWVEoN9rPK9BzjtBzf9Bzj9Ox9Sss1nJ5rOC3XcHpjjvHc8oUxlHIeG/Uc1u8PptU3A1i5Pc0F8J8O/AFwnz2PLbVOlvI+tFLDnEA5x18rXedlTnMCseelhn3FOYnJ+hWP+0o8zklMJdWvWucEzHNW7Bc5x36Rhv0i3dkvEvM9rtTxBcc4Ncc3Pcc4vVLONUodY9dweq7htFzD6ZkvpFXKeWzUc9i5XWNmF69bNczcsrG+NTj1v7lX4liC1zOKorbMEfzudkReDySe3N+nT6nUOQHrV0TWr3isXzFZv+KxfsWTu36ZH1/I/HiYXN83zGVuOehQteVNfmaanmOcluObnmOcluObnmOcnj1uabmG0xtzDZtRXCjnObYZRTnscRtmbtlOih63EscB7E+Lora8EPwsLqLcn8Wl5Jw4J2OyfsXjvhKPcxKT/WnxuK/EE2FOBuW/C5pBPrXlvpcCv3fsDxzgR4DTwFOBV9Ns158BN9143HuBj+z4GU9qn/eAjdufA3yC5qJ8N1/77/pbfsbzgb8HLumwzRcDb2l/5/07PD6VrmMH/cYvkpSv9Zh1AtOulWX73EWP564rYa0s6f9af6V97iO23HdZe99Z4Ovb265HEwq+pcfvOkaJYwnptvuG7XNfv3F7DfsgnHvvfV7Hx/8UzQeqF23c3mX9XhP4FNs/qInq0Lx2rWP/qr3vWTSv/+zGf6/a8pwulgxbCyXOCeyfl9xzAsPmZY5zAt3mJeqcQJnzYv2KyfoVj/UrHutXTNaveHLXr13n0rlNmSHXmh8vGZ5JpM4Ej8109lkybQ6aMy9fkidv2pU1weG8aU41/Ngsuasl9R27TDm+4HHIimt4XJtjHHUNQ5lj7BpOzzWclu916Vkn0hp7fCNmFBHzCTCj2KaEvqWxlNjjFmWdwPyznRQ9brnGYsmwOU+53WP0qNmfNs/+NPCzuIhyfxZ3yBLPJTbl3leWOCebcs8JWL+2sX5NL/q+ssQ52eSc5OExcTy595WI2S90z38j5XFTMGuNnbWCWVFJNdjPg9Pz8+C0XMPpuYbTs/cpLddweq7htHyvS2/MMS49XxhDieexUa6fBn5/cNHz+SsLpssALgb+BPivwHXb2063v//yHc8psU6W9D5Uy5xAWcdfx87LXOYE4p531LKvWL9isn7FY/2Kx/oVUyn1q+Y5gXrznEhZzCFLhp//2y8SO2cD+0VKqqclHXeUOL5Q1vF2iWPsGk7PNZxeKecaUOYYu4bTcw2n5XtdeuYL4yrxPLaWfpFIar5u1RhKyMaW9H8/6VODp/g390rd71Jud2nXM4Iy+9PmdIy14ne38/B6IPHk/j79IUvqmxOwfkVk/YrH+hWT9Sse61c8ueuX+XG9+fGScfKIKN83nMrcc9BNC8yb9vEz0/TsrUjLNZyeazgt13B6ruH07HFLyzWc3phr2Iwi1jm2GUVc9rgNM/dsJ0WPW6n7TMrttj/NvHAfP4uLJ/dncYcscU425Z4TMAfaxvo1vej7yhLnZJNzkof9afG4r8QTYU52nUsv2ttP7XsB+x70MeDp+548sg9x/ot+KXCzjcdc2t73sh0/45Xt/bfcuH1zQFf/nd543CXAJ2kmtot/1/6c3+74+FS6jB30H79IUr7WrusEpl8rS8YJoEpYK0v6v9bvbZ/7LprXsHIx8ELOve5vaW+/Tfv3V/Tb1E5KHcvU2/1pmvW4roZ9EI4Pgv95+/jLNm7vsn7v2D7mzFFbmE+Xee1ax1YfVlwN/AXwjTQfKNyBZl2fBa7osY1Lhq2F0uYEDs9L7jmBYfMyxzmBbvMSdU6gvHmxfsVk/YrH+hWP9Ssm61c8UerXtnPp3KbMkGvNj5cMzyRSZ4LHZDqHLJkuB82dly/JkzftyprgcN40pxp+bJbc1ZK6jl2mHl/wOGTFNTyebWMcdQ1DeWPsGk7PNZyW73XpWSfSSjW+0TKKaPkEmFHsUkLf0lhK7HGLsk5g/tnO2D1u9qftNqRHzf60efangZ/FRRTls7h9lngusSn3vrLEOdmUe07A+rWN9WtaJewrS5yTTc7J9DwmjifKvhIt+4Xu+W+UPG4qZq0XipS1glnRmaO2MB8/D07Pz4PTcg2n5xpOz96ntFzD6bmG0/K9Lr0UY1xyvjCGEs9jo1w/Dfz+4GLAz4BpM4CfAj4P3H3tttPt8y7f8ZzS6mRp70M1zAmUd/x17LzMZU4g7nlHDfuK9Ssm61c81q94rF8xlVS/ap4TqDfPiZLFdLFk+Pm//SIXipSzgf0iZ47awnxKO+4obXyhvOPt0sbYNZyeazi9ks41oLwxdg2n5xpOy/e69MwXxlfieWwt/SKR1HzdqjGUkI0t6f9+0qcGp/4390rd77ye0YWWjJM7LjBz3Mfvbsfj9UDiifJ9+n2W1DUnYP2KyPoVj/UrJutXPNaveKLUL/PjOvPjJePkERG+bzilueegmxaYN+3iZ6bp2VuRlms4PddwWq7h9FzD6dnjlpZrOL0Ua9iMIsY5thlFvF6bdfa4DTP3bGfsHrdS9xn70y60ZJxau8C8cB8/i4snymdx+yxxTjblnhMwB9rG+jWtEvaVJc7JJudkevanxeO+Ek+kOdl2Lr1on3dq34vY9aAvZnfIue6E3Sft2/7rctJ7M+DBwH8FPgjcde2++7Y/56d3PPdvgA93+B27/LP2539vh8c+un3s2zk/LNrnhPHHa92+sYP047dyQtrXCflfa8q1csJx43fmqC1v5B6/lRPGfa0XAb/TPvZDwC8ATwPeSlMo39ned7/28Xdv//7rI27j5noudSxTb/cHaN70+ip5H1xwXO37pvbxD9q4vcv6fVj7mEcft4nZHDOvh+rYz7Y/6/PAnTfuuy7w/vb+u7PbCeOvhdLmBLrPyxRzAuPPy5znBPbPS9Q5gfLmxfoVk/UrHutXPNavmKxf8USpX0PPpcfWJUM+4bj1Y348fh1cFyUTXDkh7euFvOtl3QlpX+uC7vvQrqwJDudNc6rhx2bJ25zgsUvK8QWPQ8A1PIVtYxx1DUN5Y+waTs81nJbvdelZJ9JKNb6RMopcPW65z7HNKOxx66rrWom4TmC+2c7YPW5R95kzB37eFNs95D3L/rR59qeBn8VFFOWzuJUTxt9HS5sXewniiTYnYP0C61dE0faVE5wT5yQmj4njibKvRMp+oXv+uy53HrfphOPWd6nfJ4b5v1azonlmRX4enJ6fB6flGk7PNZyevU9puYbTcw2n5XtdeinGuMR84YTj1sLQ779BnPPYiNdPA78/eGx9X1lw3Jj3zQDuRrOf/+zG7afb512+43ml1cmS3odqmRMo6/irz7zMZU4g5nlHLfuK9Ssm61c81q94rF8xlVK/ap8TqC/P2ZQ7i9l0wrjvBZuivN4Tjnudc/5uHtgvUlI9Lem4A8obXyjreBvKG2PXcHqu4fRKOddYKW2MXcPpuYbT8r0uPfOFxgnHrQX7RcYfr02lXmPG78ecc0LedQJxs7EzB35enxp8KLc5dhtz/Xt7Y2zrupqvZwTz6E+DeR1jbfK729PxeiDxRPk+/coJzglYvyKyfsVj/YrJ+hWP9SueKPXL/LiRO9eKlB+fOWrLG1HyYzAHHXNuFxw3RjXlTX5mmp69FWm5htNzDaflGk7PNZyePW5puYbTS7GGzSgauc+xzSjGccL4a2Nd7mzCHreY2c7YPW6l7jP2p6WrtQuO2+9qygvBz+IiivJZ3MoJzkm0OQFzILB+RRRtXznBOXFOYrI/LR73lXgizcm2c+lF+5xTqxsu7rChK/99++cnDjzu3cBnjvi5H+zwmA8DLwLeSBO4PBe4rL3va9o//3TL876c5gT/ZUdsz6YHA58FfvvA4x4JPB14G/CNwMc6/vwU47Vu39hB+vFbSf06If9rTblWngbcaOO2uwDfDvw/NAV43Zs7/MxNucdvZezX+gWa4vzDwEPb/z4H/BHNm8yzgK8E/rp9/KfbPy/Z8zOHrudSxzL1dl+Xc+PfR+n74DFW78W33Li9y/pdHXC8adQt6ucRwI8CtwD+HHgM8OqNx3SdVzhcx/62/fM9wJ9tPPfTwMuBH6C5kO8f7/gdKdZClDm5F818fA3NnDwEeMGOx3adlynmBMaflyhzAvA44DuA2wJXAa9tb3vrxuPG2leizgnEmZdHAj8I3Lr9+5/TfEi6OfbWrzweDzwReDbwqC33W7+mcRp4wsZtHwZuvuWx1q9p3QJ4MvCtwPWBd9HsK69ae4z1azonwK223P4cmvebddavaVyTpob9zzT7y5XAL7e3bQatUerX0HPpsXXJkM2Pj8uPU2cSUTLBldw56JSvN1LetCtrgsN5U5QaDt2yJti93x6bJW8z52OXrtlRyvGF+R6HwPAcyDW8X9dMB7aPcdQ1DHHGeN2+jMY13M9phuUtruFuhmQnUcc4yvieMDwHsU7s1zXXSDW+kTKKXD1uuc+xzSjsceuqy1qJuk5gvtnO2D1u9qftNuQ9y/60efangZ/FTSVFfxrYSzDU1P1pYC/BIaX2p8F869e6sfrTwPo11Gm6ZabR9pU5zwnYn7YSZU5OsD9tJcqcwPAsd5uU8xIp+4Xu+e+63Hncplq+Twzzf61mRfPMivw8uJ8xch8/D94vQm8ZzHcNw/BMxjXcXZ98JeoahjhjfJphWYlr+LAuuQe4hvs6oVuO4XtdP2N8Z27IGJeYL0z9/TeIcR4b9fpp4PcH1+X4jtm+DOBi4P+lmZefPPL3RamTMK/cYA5zMjRriDYn0H9eoswJxMgnIswJxJmXCHkGWL92mUsGUnr9Os3wayBZv9KY0/f45jAnJwzrMXFOxmeec7459Iv43bzt5vrdPLBfJEo9hXnlPitRxtd+kfQi5DEw3zUMMfKVOa/hdXPJSlaijPFp7BdJzX6RtE6wXyQl84Xz2S+ym/0ih/n9mEaftZJ7nUC51+DpU4MP5Tal/Ht7MO7aqfl6RhArd6wpcwS/uz2loddN38brgQzj9YzOiTInQ6/Nvo31azxezyjOvJzG6xmtRJkTGJbTbmP9GuaE8a9nBNavobye0TDmx7tFziOi5MdgDmqP23D2uKVjj1t69rilNUam4xruLte/twc9PiGmAAAgAElEQVTzXcNgj9sU7HFL64T8PW5zHt+hPW5z+vf2IGZGMcX5phnFOOxxa9jjtl/0HrdS9xn70+JkhjXlheBncVOauj8N7CU4pNT+NJhvDpSiPw2sX2Pqk39uY/0a5jT2p61EmROwP20l0pycMCyX3cZ9ZZjS+tNgwLn0AjgLnNq4/R7t7T/Q54eO6E3tdnxx+/dfa/++bad5VHvfz/T8XdcEPgq89MDjHtP+nrcAX9Lzd01hc+wg7fjlNPVrzbFWlu3PWgz8OdtEWytL0rzWVaH8FHCt9rZbtr/rNSP/rnWljmXK7b6IJnB/d8/nl74PLtqf9byOj79d+/gf3bi9y/p9Nc1Y3+C4TezlWntu/26aD1QeDtweeAbw98A/Xntc13ndZlsd+472tjfseM5T2vt/7MjftWTYWphyTmD7vFwL+BaaMGs1Tt+14/l952XKOYFh8xJhTla3vxz4fpoThTvSnEB8iObDt5Ux95WocwJx6te30wQoXwnchiZw/Bxwp7XHWb/S2FW/Vr4eeC/NCfSztjzW+jW+XfvKE4F30IS+q/9uuuVx1q/x7atfN6IJmZ5LEyZ9GU1jxO3XHmf9SmNX/bop5+8n92X767J+jW/XvvIEmmahB9J86PvP2r9v/iMfUerX0HPpFCJkyDXkx0vSZGXRMsGVJdPloLlf75I8edOurAkO501RaniXrAn61/BtWXJXS+Zx7NIlO8oxvjCP45DUOZBruFumA/3GOOcahjjnnSv7MhrX8GE58hbXcPrsxPe69DlILXVi3xp+PIdzjVTjGy2jiJhPgBnFMaJmFClEXCulrBOYd7azrk+Pm/1p2w15z7I/rRGtPw12ZwGpM0M/i9tt13HrA0nbnwb2EuwSqT8N7CWAevrTYB71a3V7qv40sH7tsm9OTnM4My1lX1kyjzmxP+2cKPXL/rRzosxJ6iwXxp2XaNkvDM9/a8oeoa7XGzFrBbMiiHOskvK7jH4enD738fPguL1lMI81nDqTqWUNQ57v/9n7lD4rcQ13yz3ANXzIvjXcJcfwve6wXN+Z2zfGc8wXxhDxPLaU66eB3x88xoL0GcCN2vu6/Pe0jedGqJNzzA1KnxNImzXkOjboOy9R5iRyPrFk2jmBOMfRUfMMqLN+zTUDKb1+pb4GkvVrt5q+x1fKnECeayA5J/uZ53Qz136RJeOc/2+K+npTiJizgf0iEOe4Y265z0qU9yv7RcZRYh4D81jDkfOVJeWsYagrK1mJUidOY7/IGOwXSct+kfTMF7qxX2Q7+0XiibhOoJy1Ei0bW5ImA9xVg1P/m3ul7ndez+hCS8ZZmwvmmTmC1zPqIsKcrG5Ped108Hogu3g9o24i7Ctds1CwfqXg9Yy6iVK/TuP1jFYi7Cupc1qwfu2T63pGYP3axesZpWN+vF3UPGKbWvLjlTnnoAvmmTfZ45aePW5p2eOW3q7XlzrTqWUNw7AeN9fwYfa4pWePW1ol9rgtt2zHMSKsYUjb41bav7cHMTOK1OebZhRlZhQRsyywx62vJXF63ErdZ+xPu9AS+9MO2XV8lLo/DfwsbpdI/WlgLwHU058G88iBUmeZYP3a5VB/GqTJP8H6tYv9ad1EqV83xf60lShzkjqXBfeVXebSnwa7z6UX7XNOrT+wq2u0f5494jkp3LL98/Ptn7ejOfB638bjrgP8YPv/b+z5u+4F3ITmgHuXfwM8FXgzcG/gr3v+rilsjh30G7+HAFdxfhj0dOAv2H6h0xzGeq1duVYa6+P3V8BjNx57F5q181XjbOZgDwUuAX6D5rUBXAl8BLhtwt9b6n7XZ7u7roPb0rzPvLnnts1tHzzkC+2f19i4/dD6vQZwZ5p188kt9/8j4AxNyPQZ4K3AN7X/fZbzD5a+lOaY4J+sPfcs8L3AH9DM8ffsuf2x7e/6ReDtwKPb7f9f135Hl3ndZVsd+wPgaprw5tpbnnNZ++dJj9/XV985gWHzsu22lwI/AfzmgW3uOy/OybnnH7Ov3B/4pfb3vIXmveumwDes/Y4x95WIcwKx6tdvAb/Tbss7gR9vt+nua7/D+jVt/QK4IfDLNA1Hf7tjm61f09avq9vftfrvI1u2yfo1bf361zTHW98HvJ7mQ5PfozkWW7F+TVu/PsL5+8kDaMK8V21sk/Vruvr1BOA/AS+hGYPV/3/dxjZFqV9Dz6VTiJAhmx/3N3WmGsEY66WEHPSQXVkT7M+bItXwJ3M4a4L+NXxbljyFSMcuXbIjxzduDlTbGPfNdKDfGOcaX4h13gmHMxrXcMy8xTWcPjtxDTdfRkiZg8xxjI9dw/fgcK6RanyjZRQR8wkwoziGPW751kpJ6wTqyXb69LiVus+knsMh71lzq7WHTN2fBuNnAV3608DP4qY8L7gxafvTwF6C1XOj5JLgvtLn/M7+tEak+vU9pO1PA+vX6rnHZtKHMlP3lWnrl/1pjUj1K3UuO5c5gWnrV8osF8adl2jZLwzPf2vrjzJrNZcfi99l3G+On1Wunh8l95njGEfKcOY4vqvnR8lkahvjqb//Z+9T+qzENdwt9wDX8JA13OV7dq7huN+Z2zfGc8wXxhDtPLa0c1i/P5hOnwzgKuA/7PjvTe1jXtP+/Y/XnhfpvX5uuUFpczJ11pDr2KDPvET7jGNu+UTEfaX2PCPinKyeX3MGkqJ+zeU7eVDHvuL3+A6bQ/1KeQ0k58Q8Zwz2ixzHfhGzorHYL7Kfnz/YL2K/SD+R1vAc8xWIda4B88tKINZ7HdgvYr/I8SKtYftFzBemZL/IhUo7hx2jX8S840JzyzvGGr/o1yPaVYNT/5t7pe53Xs8oDq9ndNhcvrs9dcaT8rrp4PVAVs+Nkm2C+0q0a7ODc+L1jIaJVL/A6xlBrPqVMqcF5yTi9YzA+rV6rtczmo758YVKyyOm/L5hBGOsl+gZaFf2uO3nZ6b2uNnj1k+kNXxD0mY6tY1x3x63uY0vzCdviTrGc8lO7HGzx630HrdD42tGsd3U59hmFDF7bQ6JlmVBWWullmynT49bqfuM/WlxlNaftu32J5O2Pw38LG713Mj9aRBzXmrJJcE5iZhlgvVr9dxjs7lU+Se4r/SpX2B/GsSqX6/D/jSI99lCylwW5jMvU9avkvrTYGD+u6AZhFMbt9+lvf0xfX7oEW4H3HzL7RcBT2y34Q/Xbn99e9tt1m67Hs0BwFnO3xmP9UyaSbnZjvt/sv35fwJc2vN3jOnYsYN+43cNmtf8i+3ffwT48JbHpTTVa+0qx1pZtj9z0eO5U4zf84Ff3fgZvw88q8f2Lun/WgG+aMttXwt8jOYN5ss37ntB+/u+oufvOyTnfrek/1j22e6u6+D72+c/qsd2QXn74KZF+7Oe1/HxX9s+/pFb7tu3fm/b3vdrW+77Uprg/j8B/wPNPP8L4OuBH+XCg4oHAp/gXBj9gPZnvxH4Zpq1cOmO27+E5iDnIRs/89mcf8C7b1771DFoxvgs8NMbt38TTcD+ceBGW563z5L+a6HvnED/efm+Lbdt7hdnge/asc275iXSnED/eckxJ7v2lW316hbt49aD4LH3lWhzAnHq1+acXJPmhPGzwB3Xbrd+TV+/fh34t+3/X8H2Y07r13T16xnAp4AP0ISNvwrcest2Wb+mrV9vA34O+BWapog30xx/rzcXWL/yHX9dG/go8Pgt22X9mq5+PR74S5qxhaaB6f3AIzZ+R5T6te9c+kx733LLfSlNkSGbH/ff36Jlql0tmS4HTZl/dbEkT960L2uC3XlTlBp+W5p6eShrgsP77bFZchdL5nfssis7yjG+MI/jkDFyIHANd13DuzId2D/G0dYwxDrvhMMZjWs4b97iGk6bnUC8MY60htdrcZ8cZKWmOnHs+HbJNVKNb7SMImI+AWYU6+xxi7lWoq0TqC/bGbPHzf607XM4pEfN/rRzUvSnwbhZQNf+NPCzuFznBbsyRojzWdwhS+ZxPrxpiv40qKuXoOb+NJhP/RrSnwbWrxRzcprDmWkp+8qSedQv+9Ma0erXiv1p2+cEpq9fQ7PcKeclWvYLh/PfaHlcamatMbNWMCtaiXL8OMZ3GcHPg1PnPuDnwSX2lsH81jD0y2SgrjUM037/D+KtYYjzXjdGVgKu4X1ruEvuAa7hserwrhzD97q835mDfmNcYr4whpLOY6NdPw38/iCMlwMsmCYD2OV0+5zLt9wX5b1+rrnBLqeJNScps4ZS5gR2z0u0zzg2Rcknlkw3JxA3M4qUZ8D861dtGcgupzm+fpX4nTyou375Pb5hTlNO/Vrp+90/58Q8J6W59oss6be+7ReJmbOB/SIrUY475pr7RHq/sl/EfpE5rWGIla8sibmGwawE4rzX2S9iv8gc1vA6+0XMF1KzX+R8tfaLmHdcyH+bcfv4jXU9oiXDsp0+NbhPj09Xpe53Xs/oQkvGyR0XxM4cwesZHbKk/OPezX1jV1bp9UDK6jWLtK8smUfOsW5XFgrWL69nZP06jdczglj1a2hOG2lOYD71a6Vvr5n1y+sZmR+fEzEXjJhHRPq+4RSmWC+Rrsm+bkHsvMket26iHrvsyo3A3gp73M4XdQ33yXRWalrDMH6P25zWMJSXt4Br2B6389nj1lhS1hpO2ePWd3zNKM6X6xzbjCJmr01JWRbU2eMWKdsZs8et1H3G/rQLLRknM1wQOy+EcY+Pbkvz+cqQ/jTws7i59KdBXb0ENfengVnmivUrzZwMyT8j7StL5lG/TmN/GsSqX/anNaJ9trCuTy4baV6WxNxX5tyfBrvPpRft7adWN1y85cm7fLj9c1sYMqZvBp4C/AHwbuBvaAb1n9IELR8CHr72+JfTnBC/CngRcH3gG4H/AlwJ/HfAezZ+x0NoBvc2wPva254OfCtwD+Aj7W0PAv6Ic6993cOAn6KZ9FcDj97ymBOawH0qx44d9Bu/szQ7xW+3v+fxwH3a/5/KFK/1r4D/s/1v5S7A64CvpnkTXXGtXDh+f8z5RfI7gTuzu3kppVcCnwbeShPa3oFmf78K+A4uXOMvpNne+wPvSrA9pe53fba76zq4H80+8ls9t620fRCabX5Q+/+rg4S7r23HR2kC/21W78VXbrlv3/q9a/vnG7c87xdo3hMfRHOQAfDO9s9/yYUHVl9NM/dn27/fmWY/+y7OXwfbbr8lTXiyOV8fBu679vd989qnjgE8Fvg64MeBe9F8wHEr4ME0a+ThNAdYU+k7J9C8J/WZl/9xy23H2DUvzkn/Odl3+6antb/jtWu3jb2vRJsTiFO/Vu5I8x53Cc2xxYOAt6zdb/2atn49nOYD0Ice2Gbr13T166U0x4zvpBnLx9HsM3egaaRYsX5NW7++nKaR4Kk0x4p3oQmzztJ8IA/Wr5zHXw+iCfnO7LjP+jVN/XoycAOazOfzNNnpE4HnbPyOKPVr37n0Re2fV2+5L6UpMmTz4/6myo+PyVVTm2K9RMpB++ZN+7Im2J03Ranht6T58O9Q1gT791s4PktOLeqxyy6Ob94cCBxj2L+GD2U6sH+Mo40vxDrv7JLRuIbz5i3RxhdireGh2QnEG+NIa3hdnxxkpaYxPnZ8u+QaqcY3WkYRMZ8AM4p19rjF63GLuE6gvmxnzB63UveZ1HM4pEettFoLZfWnwbhZQNf+NPCzuNwZ2DZRPotLKdL58KYp+tMg3rxEOr8D+9MgVv0a2p8G1q8U9et1NHOyLzN1X5m2ftmf1ohUv9bZn7Z9TmD6+jU0y51yXqJlv3A4/50ij+ua207BrDVf75tZ0Tm1fJcR6vqsEqbPfaCuMba3LOYaHpLJgGMM6b7/B/HGF+K818HwrATijXGkNdwl9wDX8FjHErtyDNdw3u/MQb8xLjFfGEO0jMKeITOA1BlAH1He6+eaG/QR6b2rq9KOv44V7TOOTXPIJ/qIlhnNLc/oI8p7Csw3AzlWpPeUMb6TB3XPi9/jSydS/VrXp8cEnBPznLTsFzmf/SL2ixzDfpH55T6R3q+6coztF1kXcQ3PLV+JdK4x16wkynsd2C8C9ov0EWkNr7NfxHwhNftFzol6DjtFv4h5R/l5x1TjF+V6RH1qcMp/c6/U/c7rGY2rpMwRvJ5RKpGOe7vat795PZB4vWbuK/Guze6ceD2jlCLVL69n1IhUv4bmtM5JvOsZgfXL6xmZH5sfHy/S9w2nMMV6iZKBQll5kz1u3UQ9dtnHMbbHbV20NTz039sDxxiG9bjNaXyhvLwFyhrj0rITiDe+EOd4bZM9bjF73Oby7+1BzIwi5b+3B2YUUXtt7HEbprZsZ8wet1L3GfvTxlVSXgjjHh/dErgGw/rTwM/ipuiF2ibKZ3EpzSmXdE5iZpnWr/HnZGj+6b4yfv2yP60RqX7Zn9aI9tnCuj65bA3zMnX9Kqk/DYadS7OgGaRTW+77e84Ps1K4jKYAvZnmJPhq4O+ANwCngUs3Hn8JTQD7QeBTwJ8AP0iz43wBuGLL77hG+7hfbP/+IzQT90/WHvO1NOPw2B3bebq9f99/2353SseOHfQbv5U/an/HN4+x8Uea4rU+H/jVjdt+H3jWxm251sqyfe6ix3OnGL97tNt3KXAdmuL3mB7bCsNeK8CPAn9KUzCvAt4L/N/ArXc8/to0hfl1PX/fITn3uyX9x7LPdndZBzekefN9cY9tgjL3wS7bdbLnud/RPuZOW+7bt36f0j5vM2y9VXv73Xb8vrcB/2rjtt/k/BOG5wO/seW5226/Zfv77rlx+xOAd7T/f2he+9SxlUtpjmfeC3yW5sDst4Cv3/OcfZb0Xwt95wT6z8uuuVp3lu0f2Oybl0hzAv3nJcec7Lt9c9uupAm7VlLtK5HmBOLUr5Vr08zDKeBJNGN7WXuf9Wva+nVbmiac263ddgXHHbNHmhOYZ/26Hs059/r4W7+mr1+fpQnj1z2Jc81J1q+8x18vB16y5Xbr17T1658Df0nTtH5Hmg+0Pgb8wNpjotSvQ+fSbwI+Adx4x/0ppc6QzY/7729T5cddc9WulkyXg6bKv7pakidv2pc1we68KUoN75I1weH9Fo7PkrtYMr9jl23ZUa7xhXkch2zbtmNzIHAN77ptZV+mA4fHONoahjjnnV0yGtdw/rzFNbz/9qHZCcQb46hruE8OslJTnTh2fA/lGqnGN2pGES2fADOKdfa4xetxO028dQL1ZTtj9rjZn3ahIT1q9qedL0V/GoybBYyRGfpZXNrzgm0ZI8T5LK6LJfM4H97crqn606CeXoKa+9Og/Po1Rn8aWL9S1y+4MDMtaV9ZMo/6ZX9aI0r92mR/2oVzAnnq19Asd6p5iZr9wv78d4o8rktuOxWz1jy9b2ZF56vlu4zg58GHbl/pm/uAnwfvu31zu6L0lsG81vCQTAbqWsMw7ff/IN4ahjjvddscm5WAa3jXbXA49wDX8JhreFuO4Xtd/u/MwfFjXGq+MIZIGUWunqEhpugZipwBwLD6fmi7TvY8t28GcGhbLt9yX5T3+rnmBrucJtac7Lt9pW/WUMqcwO55ifYZx+a2Rcknlkw3JxAvM4qYZ8D861dtGcgupzmufuV+T4F+38mDuuuX3+Mb5jTx69emvt/9c07Mc1KbY7/Ikn7r234R+0WOcWi7TvY8136RRrT3uEjvV+v6ZjgQ73g72rn/+nZFyWNgXms4Yr6yJOYaBrMSiPNet439IvaLdBF1DdsvYr4wBftFGqdJ1+8wxBT9IivmHY2cvTF9TTV+Y12PaMmwbKfP+1zKf3Ov1P3O6xldaEn/tXlou072PNfrGc3/u9u5c0rYnlV6PZC853t9ss1I+8qSeeQcMOza7JHmBMqvX17P6HxR6xd4PaMI9WtoThtpTqD8+rWpb6+Z9cvrGZkf15kfL+m/v0X5vuFUplgvka7Jfppy8iZ73LqJeuyyLTcCeyvAHrdN0dbw0H9vD+pawzBuj9uc1jCUmbeAa3jXbWCP26aaetyWlLWGd90+Ro9bn/E1o0h7jm1GcU6pvTb2uA1TW7YzZo9bqfuM/WkXWtJ/XR3arpM9zy2pP23b7WP0p4Gfxc2pPw3q6SWouT8NzDJXrF/jzskY+WekfWXJ/OoX2J8WoX7Zn9aI/NlCn1w20rwsibmvzLU/DfafSy/a7Ty15/ftfdDvAy/d9+SC3A/4HPBjNIH45uv9GZpx+LKJt6sU9wH+gSb8uWvmbUnlscC71v7+nTQ7/k02Huda2e46NIHp/Wj2s3cA18q6Rcd5HM28fnXuDVlT4n7XZR38S7aHgl3VuA/+BM2HqhfvuP/Y9fvtNAcnF22577rtfffeuP39wMPX/v5O4BFbnr/t9mu3P/MhG7c/G3hV+/81zuu6fXMCw+Zl11yt2xYCQ93zknJO9t2+8nM0AcpXbdxe85zA9PVrl98F/n37/87JtPVrSTPeV6/9d5bmeOlqmmMRqHtectevld8H/q+1v9c8J5Cnfr2Pc7Vq5aE05xjgnOQ8/roV8Pl2GzbVPC856tf7gR/euO0naJqaVqLMyb5z6RvRrKmfnXSLzplLhmx+3F/XXHUuSs9B4XDWBMflTVPX8C5ZE8xrv42YHTm+5kDHyL2GV9YzHahnjFOs4SWHM5paxhfMW8ZgdpJWzjVcSw4yxhgfO76Hco1U4xs1o5hLPgFmFEOV2GtzLHvchik92xm7x63EfabrHA7pUatx/5myPw3GzwJqzAyPlfvc1v607XL2QtWQS/aR4/xuG/vTzpm6fi2xP62LKL2c65mpczJ9/aopY+3DXDamHPUrV5Z7rKjZL8TIfw/ltnNUYm7UR5e8Ncp+OiW/y1gWc5/0cuYFNWQ4udfwurlmMn7/Lz2zkrRy9KAeyj2gnjFOvYZ35Ri1jC/M6ztz5gvpdcko5rT/jMkMYLuxM4A+IvaY1LgW1uU+T6sha+jDHpN47DGJxwwkntzvKStz/k5eH36PLx6vgRSPeU7DPMd+EftF6quD9ouUJffxdg0Zjv0iaeVew+vmmq+YlaRnv0ha9oukZ79IWuYLDfOFcdgvMox5xzmuk928HtG4StzvvJ5RGl7PaP4iZDxzv256H/aaxeO12ePxekYxmbnGY04bj9czisnrGe1mfpye+XF//nt7ZWWgK/a4lSV3dmSPmz1uQ+VewytzznQi9rjNaXzBvCU1s5P07HFLyx63hhnFOMwohimx16YPe9z6m0O2M2aPW4n7jP1paZTUn7btdvvTusl9bmt/2oXMJeMxy4xp6vq1pL7881hRrhs417y0jxz1q7aMtQ9z2Xhy1K9S+tNg/7n0or1v73XD9j3ofwM+NGjzYvkjmsXyzVvuezvw5mk3pxh3Bj4OPAx4MfMI9be5B82+cCnNweO7gcdseZxrZbfXAs+k+aDkAZm35ViX0BwovST3hrRK3u/2rYPrAh8EXjDg59e4D74Q+J099x+7fr+V5mT5i7bcdxlNLbzF2m33am+7W/v369EctN5947m7bgd4HfALG7e9E3hS+/81zuu6fXMC/edl35ys29VsWvO8pJqTfbevPIPtITDUPSeQp35t8/8Bz2v/3zmZtn7dqP2Z6/+9AfiV9v+v0T6u5nnJWb9WLgGupMkUVmqeE8hTv34FePXGbf8H8Lb2/52TfMdfp2n2kW2NHjXPS4769TfAozZuexznN/RGmJND59IPBD4D3HyyLTrfnDJk8+N+uuaqc1JyDgqHsyY4Lm/KUcMPZU0wr/02Ynbk+JoDHSP3Gl5Zz3SgnjFOsYa7ZDS1jC+Yt4zB7CStnGv4NHXkIEPHuM/4Hso1Uoxv5IxiTvkEmFH0VXKvzTHscRuu5GxnzB63kveZQ3M4tEetxv1nyv40SJMF1JYZHiv3ua39advl6oWqJZfsI8f53Tb2p50zdf2yP62bCL2cm5mpczJ9/aopY+3DXDamHPUrR5Z7rMjZL8TJf/fltnNTcm50rC55a4T9dGp+l7Es5j7p5coLaslwcq/hdXPNZPz+X3pmJWnl+NzpUO4B9Yxx6jV8mu05Ri3jC/P5zpz5wnQOZRRz2n/GZAaw3dgZQB8Re0xqXAvrcp+n1ZA19GGPSTz2mMRjBhJP7vcUmP938vrwe3zxeA2keMxzGuY5DftF5smsaDv7RcqS+3i7hgzHfpG0cq/hdXPNV8xK0rNfJC37RdKzXyQt84WG+cJ47Bfpx7zjfK6T/bwe0ThK3u+8ntH4vJ7R/EXIeOZ+3fQ+7DWLx2uzx+P1jGIyc43HnDYer2cUk9cz2s78eDrmx/347+2VyR63suTOjuxxs8dtqNxreGXOmU7EHrc5jS+Yt6RmdpKePW5p2ePWMKMYjxlFPyX32hzLHrdhSs92xupxK3mfsT9tfKX0p+273f60w3Kf29qfdiFzyXjMMmOaun7VmH8eK8J1A+ecl/aRo37VlrH2YS4bT476VUJ/Ghw+l17QjMOpfT9k34Nu0953y75bGMh9gH+gWUx3zbwtJbkV8AHgx9u/34lmDO+ZbYvSuQ5wFXA/4MeAdwDXyrpF5Xkqzfp4ee4N6elewBNo3iByKn2/27cObk9zwHPrCbdnDt5D86HAPses35sAH6M5ObgDcDvg4TQHVTenmb8Ht4/9Opp6+HmaAw9oDp4+v+V37bod4LuBzwKX06yDpwF/T7PetX9OoP+87JuT6wN3af87C/xI+///eIwXNAOp5mTf7QDPofkg7T7t71j9d/2hL2gmctSvJ9O8B98auCPNB1hfAL5lhNczBznq16YrgGf1fQEzlKN+/TvgnwJf1v68/0xTy3yfPydH/fpa4HM05xVfATwE+DvgkWO8oBnIVb8uomniePIYL2JmctSvM8BfAd9G817/YOAjwM8Nfzmjin4uPZcM2fy4vxpz1dJz0C5ZE3TPm3LU8NqyJrOjtMyB0suxhmvLdHKcd266gvlmNOYt6ZmdpJXrva6mHGToGPcZ3zNMn2tEzvw3Gz8AAAmxSURBVCjmkk+AGUVfpffaHKPGLGZspWc7Y/S4lb7PHJrDyO9ZUU3ZnwZpsoDaMsNj5TgvqClj7CvH+XBNuWQfOc7vassyj5Ur11h3BfPNPvvKUb9qy0yPlaN+1ZSx9mEuG1OO+nWG+D1q0c+jI+S/NeW2pedGxzJv3c7vMpbF3Ce9HMfbNWU4uY6da8pk/P5femYlaeVYw7XlHrl6UGvJMWr5zpz5wjRqyijGZgaw3dgZQB/2mMRj1hCTPSbx2GMSjxlIPH4nLya/xxeP10CKxzwnhgh5Tk1ZjP0iZkVgv0hpzHDSs18kLftF0jMrSc9+kbTsF0nPfpG0zBdiiJAvjKGmjGJM5h3mHcfyekTDlb7feT2j8Xk9o/nLlfHUllUey16zeLw2ezxezygmM9d4zGnj8XpGMXk9o+2in0ebH6vGTLX0DBTscSuNPW7p2eOWVo41XFumY49beuYtaZmdpGePW1r2uMVgRlG30nttjlVjHjOmOWQ7Q3vcSt9n7E8bXyn9afturykv7CvHeUFtOeOxzCXjMcuMKVeuse4K5p1/HitH/aopL+0jR/2qLWPtw1w2nhz16wzx+9Pg8Ln0guaY/tS+H3LoQb8HPKrHxkVyZ+DjNCfRLwZemndzinEp8Hbg5zdu/w3gNdNvziReCzyT5oDlAZm3pUQPA66mKdbqZw77netgXHcGPsr4X57/BuAPacLYvwVeAXxJe9+P0Rx8fQB4Lk3A/Pa15/4QzcHWpl23rzwCOKH5sOxPacJrnbNvTqDfvOybkwXNMeDmf2f6v4TZSTEn+26H7XNyluagX42p69cZmpP1q4C/Bn4XuP+QFzBDU9evTVdgCLxp6vr1a8AHaT70/QDwQuCrBr2Cecpx/PVtwJ8BnwHeCTwauEbvVzA/OerX/Wje228zZMNnbOr6dQOaJpX3AZ+macL5GeCSIS+iUqVnyObHw9WWq5acf+XImiDN8U5tWZPZUVrmQOlNvYbPUF+mk+O8c90VzDujMW9Jz+wkrRzvdbXlIEPGuM/4mmtcqPR8Aswo+ppDr82xastixlZytjOGOewztc/h2ErKDO1PG2bq84IFdWWMfU19PlxbLtnH1Od3Z6gvyzxWjlxj3RXMO/vsa+r6VWNmeqyp6xfUlbH2YS4b09T1yyx3HDnz35py2znkRn2Yt57P7zKWydwnvamPt2vLcHIcO5+hrkzG7/+lZ1aSVo41XFvukeN4raYcw+/MxVB6f1lNGUUqZgDnS5UB9GGPSTxmDTHZYxKPPSbxmIHE43fyYvJ7fPF4DaR4zHNisF9kGvaLmBWB/SKlMsNJz36RtOwXSc+sJD37RdKyXyQ9+0XSMl+IwX6ROpl3mHf04bVshpnDfucaGJfXM6pHjoxnQX1Z5bHsNYvHa7PH4/WMYjJzjcecNh6vZxST1zMqk/mxastUS8+/7HErkz1u6dnjltbUa/gM9WU69rilZ96SltlJeva4pWWPWwxmFHWaQ69NH7XlMWMqPdsZag77TO1zOLaS+tP23Q515YV9TX1esKC+nPFY5pLxmGXGlCPXWHcF888/jzV1/aotL+0jx/FXbRlrH+ay8Uxdv+aS5S5o1tWpIQ+6N/CGMbdqYreiWRg/3v79TsAXgHtm2yJF9lSa9fHy3BtSqFcCz869EcrOdTCupwI/mXsjJEmSJEnao+QM2fx4HLXlqiXnX2ZNkiRJkuaq5HwCzCh0nNqymLGVnO2o4RyOy8xQkiRJkhRdrvzX3LYO5q3nMyuSJEmSNFcl95eZUYzDDOB8ZgCSJEmSpOjsF1FKZkXnMyuSJEmSNFf2i6gm5h3DeC0buQbGZeYoSZIkSYrO/Fi1Zaql51/mTZIkSZLmyoxCNaktjxlT6dmOnMOxmRdKkiRJkg5ZAGeBU0Mf9HzgbmNt1YQuBd4O/PzG7b8BvGb6zVEBHgZcDdwh94YU5CLgZsDjgSuBG+fdHGXiOkjjBsDrgevm3hBJkiRJkg4oMUM2Px5PDbnqHPIvsyZJkiRJc1diPgFmFDpeDVnM2OaQ7dTOOUzDzFCSJEmSVIqp819z23qYt55jViRJkiRp7krsLzOjGI8ZwDlmAJIkSZKkUtgvolTMis4xK5IkSZI0d/aLqBbmHcfzWjZyDaRh5ihJkiRJKoX5cd1qyFTnkn+ZN0mSJEmaOzMK1aKGPGZMc8l2auYcpmFeKEmSJEnqYgGcBU4NfdBNgZcAF4+0YVJUrwSenXsjCrMAvkDzgck98m6KMlrgOkjhacA9c2+EJEmSJEkdmCHXrYZcdUH5+ZdZkyRJkqS5M59QLWrIYsa2oPxsp3YLnMMUzAwlSZIkSaUw/1Uq5q3nmBVJkiRJmjvzhbqZAZxjBiBJkiRJKoV5jlIxKzrHrEiSJEnS3JkvqBbmHcdb4LVsarfANZCCmaMkSZIkqRTmx3WrIVNdMI/8y7xJkiRJ0tyZUagWNeQxY1owj2ynZgucwxTMCyVJkiRJXSyAs8CpwQ9q7/++MbZKCuYi4GbA44ErgRvn3RxJAuBOwOW5N0KSJEmSpCOYIdfFXLUsZk2SJEmSamE+obkyi5E0NjNDSZIkSVJpzH81FvPWC5kVSZIkSaqF+UJdzAAuZAYgSZIkSSqNeY7GYlZ0IbMiSZIkSbUwX9BcmXdIisbMUZIkSZJUGvPjupiplse8SZIkSVItzCg0V+YxksZkXihJkiRJ6moBnKXJXne6G/AO4GXAi4HlnsdeZ6QNkyJZAF8A3g7cI++mSNL/z/dcSZIkSVKJPJ+txwJz1ZK4b0qSJEmqiedAmqMFZjGSxuX7pSRJkiSpRJ7PagwLzFs3uW9JkiRJqonnQPVYYAawyfUvSZIkSSqR57MawwKzok3uW5IkSZJq4jmQ5miBeYekWHy/lSRJkiSVyPPZeiwwUy2N+6ckSZKkmngOpDlaYB4jaTy+V0qSJEmSDvkh4MXAK4B3AHfIuzmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJElSXP8Nqm1iRteF1woAAAAASUVORK5CYII=\n", "text/latex": [ - "$$\\left[\\begin{matrix}\\frac{- f_{curr 0} + f_{curr 2} - f_{curr 3} + f_{curr 5} - f_{curr 6} + f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}}\\\\\\frac{f_{curr 0} + f_{curr 1} + f_{curr 2} - f_{curr 6} - f_{curr 7} - f_{curr 8}}{f_{curr 0} + f_{curr 1} + f_{curr 2} + f_{curr 3} + f_{curr 4} + f_{curr 5} + f_{curr 6} + f_{curr 7} + f_{curr 8}}\\end{matrix}\\right]$$" + "$$\\left [ \\left(- \\frac{3 u_{x}^{2}}{2} - 3 u_{x} - \\frac{3 u_{y}^{2}}{2} + 3 u_{y} + \\frac{9 \\left(- u_{x} + u_{y}\\right)^{2}}{2} + 1\\right) \\left(\\frac{f_{curr 0}}{36} + \\frac{f_{curr 1}}{36} + \\frac{f_{curr 2}}{36} + \\frac{f_{curr 3}}{36} + \\frac{f_{curr 4}}{36} + \\frac{f_{curr 5}}{36} + \\frac{f_{curr 6}}{36} + \\frac{f_{curr 7}}{36} + \\frac{f_{curr 8}}{36}\\right), \\quad \\left(- \\frac{3 u_{x}^{2}}{2} + 3 u_{y}^{2} + 3 u_{y} + 1\\right) \\left(\\frac{f_{curr 0}}{9} + \\frac{f_{curr 1}}{9} + \\frac{f_{curr 2}}{9} + \\frac{f_{curr 3}}{9} + \\frac{f_{curr 4}}{9} + \\frac{f_{curr 5}}{9} + \\frac{f_{curr 6}}{9} + \\frac{f_{curr 7}}{9} + \\frac{f_{curr 8}}{9}\\right), \\quad \\left(- \\frac{3 u_{x}^{2}}{2} + 3 u_{x} - \\frac{3 u_{y}^{2}}{2} + 3 u_{y} + \\frac{9 \\left(u_{x} + u_{y}\\right)^{2}}{2} + 1\\right) \\left(\\frac{f_{curr 0}}{36} + \\frac{f_{curr 1}}{36} + \\frac{f_{curr 2}}{36} + \\frac{f_{curr 3}}{36} + \\frac{f_{curr 4}}{36} + \\frac{f_{curr 5}}{36} + \\frac{f_{curr 6}}{36} + \\frac{f_{curr 7}}{36} + \\frac{f_{curr 8}}{36}\\right), \\quad \\left(3 u_{x}^{2} - 3 u_{x} - \\frac{3 u_{y}^{2}}{2} + 1\\right) \\left(\\frac{f_{curr 0}}{9} + \\frac{f_{curr 1}}{9} + \\frac{f_{curr 2}}{9} + \\frac{f_{curr 3}}{9} + \\frac{f_{curr 4}}{9} + \\frac{f_{curr 5}}{9} + \\frac{f_{curr 6}}{9} + \\frac{f_{curr 7}}{9} + \\frac{f_{curr 8}}{9}\\right), \\quad \\left(- \\frac{3 u_{x}^{2}}{2} - \\frac{3 u_{y}^{2}}{2} + 1\\right) \\left(\\frac{4 f_{curr 0}}{9} + \\frac{4 f_{curr 1}}{9} + \\frac{4 f_{curr 2}}{9} + \\frac{4 f_{curr 3}}{9} + \\frac{4 f_{curr 4}}{9} + \\frac{4 f_{curr 5}}{9} + \\frac{4 f_{curr 6}}{9} + \\frac{4 f_{curr 7}}{9} + \\frac{4 f_{curr 8}}{9}\\right), \\quad \\left(3 u_{x}^{2} + 3 u_{x} - \\frac{3 u_{y}^{2}}{2} + 1\\right) \\left(\\frac{f_{curr 0}}{9} + \\frac{f_{curr 1}}{9} + \\frac{f_{curr 2}}{9} + \\frac{f_{curr 3}}{9} + \\frac{f_{curr 4}}{9} + \\frac{f_{curr 5}}{9} + \\frac{f_{curr 6}}{9} + \\frac{f_{curr 7}}{9} + \\frac{f_{curr 8}}{9}\\right), \\quad \\left(- \\frac{3 u_{x}^{2}}{2} - 3 u_{x} - \\frac{3 u_{y}^{2}}{2} - 3 u_{y} + \\frac{9 \\left(- u_{x} - u_{y}\\right)^{2}}{2} + 1\\right) \\left(\\frac{f_{curr 0}}{36} + \\frac{f_{curr 1}}{36} + \\frac{f_{curr 2}}{36} + \\frac{f_{curr 3}}{36} + \\frac{f_{curr 4}}{36} + \\frac{f_{curr 5}}{36} + \\frac{f_{curr 6}}{36} + \\frac{f_{curr 7}}{36} + \\frac{f_{curr 8}}{36}\\right), \\quad \\left(- \\frac{3 u_{x}^{2}}{2} + 3 u_{y}^{2} - 3 u_{y} + 1\\right) \\left(\\frac{f_{curr 0}}{9} + \\frac{f_{curr 1}}{9} + \\frac{f_{curr 2}}{9} + \\frac{f_{curr 3}}{9} + \\frac{f_{curr 4}}{9} + \\frac{f_{curr 5}}{9} + \\frac{f_{curr 6}}{9} + \\frac{f_{curr 7}}{9} + \\frac{f_{curr 8}}{9}\\right), \\quad \\left(- \\frac{3 u_{x}^{2}}{2} + 3 u_{x} - \\frac{3 u_{y}^{2}}{2} - 3 u_{y} + \\frac{9 \\left(u_{x} - u_{y}\\right)^{2}}{2} + 1\\right) \\left(\\frac{f_{curr 0}}{36} + \\frac{f_{curr 1}}{36} + \\frac{f_{curr 2}}{36} + \\frac{f_{curr 3}}{36} + \\frac{f_{curr 4}}{36} + \\frac{f_{curr 5}}{36} + \\frac{f_{curr 6}}{36} + \\frac{f_{curr 7}}{36} + \\frac{f_{curr 8}}{36}\\right)\\right ]$$" ], "text/plain": [ - "⎡ -f_curr_0 + f_curr_2 - f_curr_3 + f_curr_5 - f_curr_6 + f_cur\n", - "⎢─────────────────────────────────────────────────────────────────────────────\n", - "⎢f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + \n", - "⎢ \n", - "⎢ f_curr_0 + f_curr_1 + f_curr_2 - f_curr_6 - f_curr_7 - f_curr\n", - "⎢─────────────────────────────────────────────────────────────────────────────\n", - "⎣f_curr_0 + f_curr_1 + f_curr_2 + f_curr_3 + f_curr_4 + f_curr_5 + f_curr_6 + \n", - "\n", - "r_8 ⎤\n", - "───────────────────⎥\n", - "f_curr_7 + f_curr_8⎥\n", - " ⎥\n", - "_8 ⎥\n", - "───────────────────⎥\n", - "f_curr_7 + f_curr_8⎦" + "⎡⎛ 2 2 2 ⎞ \n", + "⎢⎜ 3⋅uₓ 3⋅u_y 9⋅(-uₓ + u_y) ⎟ ⎛f_curr_0 f_curr_1 \n", + "⎢⎜- ───── - 3⋅uₓ - ────── + 3⋅u_y + ────────────── + 1⎟⋅⎜──────── + ──────── +\n", + "⎣⎝ 2 2 2 ⎠ ⎝ 36 36 \n", + "\n", + " \n", + " f_curr_2 f_curr_3 f_curr_4 f_curr_5 f_curr_6 f_curr_7 f_curr_8⎞ \n", + " ──────── + ──────── + ──────── + ──────── + ──────── + ──────── + ────────⎟, \n", + " 36 36 36 36 36 36 36 ⎠ \n", + "\n", + "⎛ 2 ⎞ \n", + "⎜ 3⋅uₓ 2 ⎟ ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3 f_\n", + "⎜- ───── + 3⋅u_y + 3⋅u_y + 1⎟⋅⎜──────── + ──────── + ──────── + ──────── + ──\n", + "⎝ 2 ⎠ ⎝ 9 9 9 9 \n", + "\n", + " ⎛ 2 2 \n", + "curr_4 f_curr_5 f_curr_6 f_curr_7 f_curr_8⎞ ⎜ 3⋅uₓ 3⋅u_y \n", + "────── + ──────── + ──────── + ──────── + ────────⎟, ⎜- ───── + 3⋅uₓ - ────── \n", + " 9 9 9 9 9 ⎠ ⎝ 2 2 \n", + "\n", + " 2 ⎞ \n", + " 9⋅(uₓ + u_y) ⎟ ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3 f_cu\n", + "+ 3⋅u_y + ───────────── + 1⎟⋅⎜──────── + ──────── + ──────── + ──────── + ────\n", + " 2 ⎠ ⎝ 36 36 36 36 3\n", + "\n", + " ⎛ 2 ⎞\n", + "rr_4 f_curr_5 f_curr_6 f_curr_7 f_curr_8⎞ ⎜ 2 3⋅u_y ⎟\n", + "──── + ──────── + ──────── + ──────── + ────────⎟, ⎜3⋅uₓ - 3⋅uₓ - ────── + 1⎟\n", + "6 36 36 36 36 ⎠ ⎝ 2 ⎠\n", + "\n", + " \n", + " ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3 f_curr_4 f_curr_5 f_curr_6 \n", + "⋅⎜──────── + ──────── + ──────── + ──────── + ──────── + ──────── + ──────── +\n", + " ⎝ 9 9 9 9 9 9 9 \n", + "\n", + " ⎛ 2 2 ⎞ \n", + " f_curr_7 f_curr_8⎞ ⎜ 3⋅uₓ 3⋅u_y ⎟ ⎛4⋅f_curr_0 4⋅f_curr_1 4⋅f_c\n", + " ──────── + ────────⎟, ⎜- ───── - ────── + 1⎟⋅⎜────────── + ────────── + ─────\n", + " 9 9 ⎠ ⎝ 2 2 ⎠ ⎝ 9 9 9\n", + "\n", + " \n", + "urr_2 4⋅f_curr_3 4⋅f_curr_4 4⋅f_curr_5 4⋅f_curr_6 4⋅f_curr_7 4⋅f_c\n", + "───── + ────────── + ────────── + ────────── + ────────── + ────────── + ─────\n", + " 9 9 9 9 9 9\n", + "\n", + " ⎛ 2 ⎞ \n", + "urr_8⎞ ⎜ 2 3⋅u_y ⎟ ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3\n", + "─────⎟, ⎜3⋅uₓ + 3⋅uₓ - ────── + 1⎟⋅⎜──────── + ──────── + ──────── + ────────\n", + " ⎠ ⎝ 2 ⎠ ⎝ 9 9 9 9 \n", + "\n", + " ⎛ 2 \n", + " f_curr_4 f_curr_5 f_curr_6 f_curr_7 f_curr_8⎞ ⎜ 3⋅uₓ 3⋅\n", + " + ──────── + ──────── + ──────── + ──────── + ────────⎟, ⎜- ───── - 3⋅uₓ - ──\n", + " 9 9 9 9 9 ⎠ ⎝ 2 \n", + "\n", + " 2 2 ⎞ \n", + "u_y 9⋅(-uₓ - u_y) ⎟ ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3 \n", + "──── - 3⋅u_y + ────────────── + 1⎟⋅⎜──────── + ──────── + ──────── + ──────── \n", + "2 2 ⎠ ⎝ 36 36 36 36 \n", + "\n", + " ⎛ 2 \n", + " f_curr_4 f_curr_5 f_curr_6 f_curr_7 f_curr_8⎞ ⎜ 3⋅uₓ 2 \n", + "+ ──────── + ──────── + ──────── + ──────── + ────────⎟, ⎜- ───── + 3⋅u_y - 3\n", + " 36 36 36 36 36 ⎠ ⎝ 2 \n", + "\n", + " ⎞ \n", + " ⎟ ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3 f_curr_4 f_curr_5 f\n", + "⋅u_y + 1⎟⋅⎜──────── + ──────── + ──────── + ──────── + ──────── + ──────── + ─\n", + " ⎠ ⎝ 9 9 9 9 9 9 \n", + "\n", + " ⎛ 2 2 \n", + "_curr_6 f_curr_7 f_curr_8⎞ ⎜ 3⋅uₓ 3⋅u_y 9⋅(uₓ - u_y\n", + "─────── + ──────── + ────────⎟, ⎜- ───── + 3⋅uₓ - ────── - 3⋅u_y + ───────────\n", + " 9 9 9 ⎠ ⎝ 2 2 2 \n", + "\n", + " 2 ⎞ \n", + ") ⎟ ⎛f_curr_0 f_curr_1 f_curr_2 f_curr_3 f_curr_4 f_curr_5 f_c\n", + "── + 1⎟⋅⎜──────── + ──────── + ──────── + ──────── + ──────── + ──────── + ───\n", + " ⎠ ⎝ 36 36 36 36 36 36 \n", + "\n", + " ⎤\n", + "urr_6 f_curr_7 f_curr_8⎞⎥\n", + "───── + ──────── + ────────⎟⎥\n", + "36 36 36 ⎠⎦" ] }, - "execution_count": 42, + "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "u = Matrix([u_x, u_y])\n", - "u" + "f_eq = []\n", + "\n", + "for i, c_i in enumerate(c):\n", + " f_eq_i = w[i] * rho * ( 1\n", + " + c_i.dot(u) / c_s**2\n", + " + c_i.dot(u)**2 / (2*c_s**4)\n", + " - u.dot(u) / (2*c_s**2) )\n", + " f_eq.append(f_eq_i)\n", + "\n", + "f_eq" ] }, { "cell_type": "code", - "execution_count": 43, + "execution_count": 28, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAZU8AAAAqCAYAAADTjjwVAAAABHNCSVQICAgIfAhkiAAAIABJREFUeJzs3Xu0LVtd2PnvhXsVBEVAedhCrrb4aC/qULAFL2YrMUNJ2/Ic2sbHSsTWiJp02rSRNs0xIwTQmA4+SNqYpEQdRtpHG1REwLeCMhRBhKEgbuUtCgZUlNfpP2qtnHXWWY+qWjVr/n5V388YZ5xz9lp779+sWat+q+b8zblAkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJ6uEhwOuA3wRetP7z4A7fd4eSQUmSJEmSJElSRfcDfh54BfBS4NFVo5EkSZIkReS9oyRJmoubagcgSZIkSdKZXO8qSdrlfK4kSd2YMyVJUlbWPEmSJElaKutkJEmSJEmSJEVkTaokSZIk6RTvHSVJ0ly4vlGSJEmSlJ1rVSUpL+ddJUmSJEn7eL84b9YrSZIkSVoy61wkSdscC5UkaTjz6Pw5ryhJkiRJys75YUnKyzFoSZIkSdIu7xXnz3olSZIkSUtljYskSZIkSZKkObLeR5IkSZLUhfePkiRpLlwjKUmSJEnKzvWukqRtzuVKktSdeVOSJGVlzZMkSZKkpbJORpIkSZIkSVI01qNKkiRJkrrw/lGSJM2F6xslSZIkSV18GvCi9Z+XAK8DHnzoyRfAVeD2Hr/gkcBnDI9PkiRJkiRJkkK7L/CJ63/fC3gt8AH1wpEkSZIkBeS9oyRJmosnATfXDkKSJEmSpIFc7ypJ2sf5XEmSujFnSpKkrKx5kiRJkrRE1slIkiRJkiRJisqaVEmSJEnSKd47SpKkuXB9oyRJkiQpM9eqSlJuzrtKkiRJkvbxfnHerFeSJEmStFTWuUiSdjkWKknScObR+XNeUZIkSZKUmfPDkpSbY9CSJEmSpF3eK86f9UqSJEmSlsgaF0mSJEmSJElzZb2PJEmSJKkL7x8lSdJcuEZSkiRJkpSZ610lSbucy5UkqTvzpiRJysqaJ0mSJElLZJ2MJEmSJEmSpIisR5UkSZIkdeH9oyRJmgvXN0qSJEmS+roArgK3D37Cjs8FnnJuVCfcqfDPVz+PB14DvAf4t5VjGSp7G7LHD/nbYPzS/MzhdZG9DcZfX/Y2ZI9/7m4BbqodhFK5O/Bm4L8/8+dkvzZkjx/yt8H466vRhh8G/vGJ57wMuP8EsUhdmDdb2eOH/G0w/vqi5k0pimM5M/s1wPjry96G7PFDzDzovaMyM2/GlT1+yN+G7PFDzLypWB4EfE/tICRJkiRJGmCK9a7qxpqmVvb4IX8bsscP+dsQdUza+VxldSjHZr9WQP42GH992dtgzpTOZ56MK3v8kL8Nxj+MNU+5WPMkSZIkaWncFz4fx2jqyh4/5G9D9vhhHm2QxjSH10T2Nhh/fdnbkD3+uXNPd/XlOo5W9vghfxuyxw/522BNqnSaebOVPX7I34bs8UP+NkTNm1IU7k0bV/b4IX8bsscPMfOg947KzLwZV/b4IX8bsscPMfOmYnF9oyRJkiQpK/d078/ao2uytyF7/JC/DcY/jPOuisw8eU32NmSPH/K3wfiHcZ5VmbifelzGX1/2NkStNfJ+cV6sV5IkSZK0RNa5xOCc4DXZ25A9fsjfBuMfxtoZLYl7GMSVPX7I34bs8UPMeUXz6Pw4ryhJkiRJysr54f6cy21ljx/yt8H464s4/gyOQasu82Qre/yQvw3GX1/UPClF4R4IcWWPH/K3wfiHsV5peaxXkiRJkrQ0U9S43Knwz1c/jhPVl70N2eOH/G3IHr9UQvbXhfHXl70N2eOHebRhzm4BbqodhNKw5vWa7G3IHj/kb4PxD2O9jzIxb16TvQ3Z44f8bTD+YVwvokxcLxKX8deXvQ1R10x6/6jM5po3s8cP+duQPX7I34aoeVOxuEZSkiRJkpSVe7rHMFZNEzgeV1v2+CF/G7LHDzHHpZ3LVVZ+vnds2dtg/PVFzJlg3lQe5snYsrfB+OuLmicVizVPkiRJkpbGfeFzcYymPuOvL3sbsscvlZD9dZE9fsjfBuOvbw5tmDP3dFcfruO4xvjry96G7PFDzFob61EViXnzGuOvL3sbsscPMfOmFMlc96aF/G0w/vqytyFqDvT+UZmZN+My/vqytyFq3lQsrm+UJEmSJPV1AVwFbh/8hC0fCfwK8H4jBHbIo4HPLvjz1c/H0g5YPQa4L3CXuuEMkr0N2eOH/G0wfml+5vC6yN4G468vexuyx78EHwg8qXYQCuEJtAV7b1//eSHwd/Y871uB7z3zd2W/NmSPH/K3wfjrq9WGTwDeCtztwOMPAl6JmzKoPPNmd9njh/xtMP76ouZNaQrn5szs1wDjry97G7LHDzHzoPeOisq8afy1ZW9D9vghZt5UTFeAf1A7CEmSJEmSephivausaeoje/yQvw3Z44f8bYg6Ju18riI6J8dmv1ZA/jYYf33Z22DOlI4zT+ZuQ/b4IX8bjH84a57yuYI1T5IkSZKWwX3h83GMpq7s8UP+NmSPH+bRBmlMc3hNZG+D8deXvQ3Z418C93TXhus4usseP+RvQ/b4IX8brEnV0pk3u8seP+RvQ/b4IX8bouZNaQruTWv8tWVvQ/b4IWYe9N5RUZk3jb+27G3IHj/EzJuK6Qqub5QkSZIk5eKe7tez9qif7G3IHj/kb4PxD+e8q2owT/aTvQ3Z44f8bTD+4ZxnVQTup567DcZfX/Y2RK018n5xnq5gvZIkSZKk5bDOpTznBPvJ3obs8UP+Nhj/cNbOaA7cw8D4a8vehuzxQ8x5RfPofF3BeUVJkiRJUi7OD1/PudzusscP+dtg/PVFHH8Gx6BVjnmyu+zxQ/42GH99UfOkNAX3QMjdhuzxQ/42GP9w1ist0xWsV5IkSZK0DFPUuDwa+OyCP1/9OE5UX/Y2ZI8f8rche/xSCdlfF8ZfX/Y2ZI8f5tGGuftA4Em1g1B11rz2k70N2eOH/G0w/uGs91EE5s1+srche/yQvw3GP5zrRRSB60Vyt8H468vehqhrJr1/VFRLzpvZ44f8bcgeP+RvQ9S8qZiu4BpJSZIkSVIu7uk+jS7jzGPUNIHjcbVljx/ytyF7/BBzXNq5XEXk53vnjh/yt8H464uYM8G8qTjOuRfNfo3IHj/kb4Px1xc1TyqmK1jzJEmSJGkZ3Bc+F8do6jP++rK3IXv8UgnZXxfZ44f8bTD++ubQhrlzT3fBtHvTQv5rg/HXl70N2eOHmLU21qNqKubNfoy/vuxtyB4/xMyb0lSWvDct5G+D8deXvQ1Rc6D3j4rKvJm7DcZfX/Y2RM2biukKrm+UJEmSJHV3AVwFbh/8hC0vAD5njKgO+JvANxf8+ervnwC/UTuIHQ3tAElXEdvQR8T4G+yD2hq690HE+KXaIr4uGnJfWxtyx99X9vghfxuyx78UjwC+oXYQqu7zac+FBwAfDTwZeDdtwcLGBwB/BjzszN8V8drQsKz37hHb0LCcPogYf0Oe94m/QVvIueuewCuAh04bjhbKvLmcazbEbEPDcvogYvwN+fOmNJVzc2b2a4Dx15e9Ddnjh3h50HtHRTa3vNmQ533zIQ3mzUyyxw/x8qbiujPwu7Q5Q5IkSZKkDEqvd1XLmibHdGtrsA9qa8jRB9YCK5tzcmz2awXkb0PE+PvIHj/kb4M5UzrOPJm7DRHj7yt7G4z/PNY85WLNkyRJkqSlcF/4fGqPcezTsJxxvojxN+QfK+4jYvwNy+oDaWwRXxMN+V/XDebnTLK3IXv8S+Ge7gLXcTQsJz9CzDY02Ae1NeToA2tSFYF5M8f1YiwR29BgH9TWkKMPrE9VbXPbmxbyvP4Pacgdf1/Z25A9foiXB713VGTmTeOvLXsbsscP8fKm4nJ9oyRJkiQpG/d0v561R/nX5fQRMf4G+6C2hjxzEM67amrmSa/RtTXYB7U15M6T0pTcTz1WGxpyx99X9vghfxsi1hp5vzhf1itJkiRJWhLrXMpzTtCxxNoa7IPaGnLPCToWqkzcw8D4a8vehuzxQ7x5RfPovDmvKEmSJEnKxvnh6zmX6/htTQ3O49bWkKcPXNeiGsyT5snaGpbTBxHjb8ifJ6WpuAdC7jZEjL+v7G0w/vNYr7Q81itJkiRJWorSNS5/E/jmgj9f/dUeZ9mnwbHG2hrsg9oaltUH0tgivi4acr+uG3LH31f2NmSPH+bRhiV4BPANtYNQVda85q/36SNi/A32QW0Ned4nWu+j2sybXrNra7APamvInTelKbleJFYbGnLH31f2+CF/GyKumfT+UZEtOW9mjx9itqGP7PFD/jZEzJuKyzWSkiRJkqRs3NN9GqfGmceqaYKY43ENucel+4gYf4Pj6rU15OkDa4GViZ/vnTt+iNmGPoy/vmg5E8ybiuWce9GI14gG82Qmxl9fxDypuKx5kiRJkrQU7gufS8QxmgbH+WpryD1W3FfENjQsqw+ksUV8XTTkfl03mJ8zyR4/zKMNS+Ce7ppyb1qIeW1oyJ3j+4gYf4PvUWpryNMHruNQbUvPmw15rhdjiBh/g31QW0OePrBGVbUteW9ayN+GiPH3kT1+yN+GiDnQ+0dFZt7M3YaI8feRPX7I34aIeVNxub5RkiRJktTHBXAVuH3wE9a+EPj1saLa40OBXwLev+DvUD+/R3tubP78SN1w/puG7oO3UdvQVdT4G+yD2hq69UHU+KWaor4uGnJfWxtyx99H9vghfxuyx78038/p+00tz1uBr9z6/+PWX7vpjJ8Z9drQsJz37lHb0LCMPogaf0Oe94lPAn5l52vvD/wi8CUTxyJtM2/eKGr8fURtQ8My+iBq/A2586ZUW9ecWfv1c0hD7mtwQ+74+8jehuzxQ/027OZB7x2VUea82ZDnffMhDebNLLLHD/Xb4P1jPl8FPKd2EJIkSZIkdVB6vauOs6bpRlHj7yNqGxrsg9oacvSBtcCagy45tvZr7ZAG53KzyB4/5G9D7fjNmcrKPFlXwzLyJORvg/Gfz5qnfKx5kiRJkjR37gufT4Qxjn0aljHOFzX+hvxjxV1Fjb9hOX0gjS3qa6Ih/+u6wfycRfY2ZI9/adzTXfu4juNGUePvI2obGuyD2hpy9IE1qYrKvHmjqPH3EbUNDfZBbQ05+sD6VEWUeW9ayPP6P6Qhd/x9ZG9D9vihfht286D3jsrIvFlXQ+74+8jehuzxQ/02eP+Yj+sbJUmSJElZuKd7N9Ye7Re1DV1Fjb/BPqitIc8chPOuisA8uV/UNnQVNf4G+6C2hrx5UorA/dTracgdfx/Z44f8bagdv2s8l8l6JUmSJElLYJ1LPc4J7he1DV1Fjb/BPqitIe+coGOhmgP3MKirIXf8fWRvQ/b4oX4bzKPL5LyiJEmSJCkL54e7cS73RlHj7yNiGxqcx62tIU8fuK5FUZgnbxQ1/j6itqFhGX0QNf6G3HlSqs09EOpqWEYOgfxtMP7zWa+0TNYrSZIkSZq70jUuHwr8Eu19tGKIMM6yT4NjjbU12Ae1NSynD6SxRX1dNOR+XTfkjr+P7G3IHj/Mow1L8v3A7bWDUCjWvO4XtQ1dRY2/wT6orSHP+0TrfRSReXO/qG3oKmr8DfZBbQ1586YUgetF6mnIHX8f2eOH/G2oHf++HOj9ozJaQt7MHj/EbUNX2eOH/G2oHb/3jjm5RlKSJEmSlIV7ute1Pc48Rk0T1B/POqQh97h0V1Hjb3BcvbaGPH1gLbCy8/O962rIc707l/HXV7sN1kApqy73orVfX4c0mCezMP76arfBmqecrHmSJEmSNHfuC59L7fGNQxoc56utIfdYcR9R29CwnD6Qxhb1ddGQ+3XdYH7OInv8MI82LIl7umtXib1pIe61oSF3ju8qavwNvkeprSFPH7iOQxEtKW825LlenCtq/A32QW0NefrAGlVFtIS9aSF/G6LG31X2+CF/G2rH73pGzYV5s64G82YW2dtQO37vHXNyfaMkSZIkqasL2jGHg3XCJ59AOyj5SuCLRwxs1/fSLmxVHPcCfhf4P4H7AB9UN5z/pqH7AHTUNnQVNf4G+6C2hm59EDV+qaaor4uG3NfWhtzx95E9fsjfhuzxL839gZcCtxx4/BuBFwNvB94CPBu4bZrQVMEdgf8FeBfwwK2vPx143pk/O+q1oWE5792jtqFhGX0QNf6GPO8TP4f2+nTn9f9vAn6Q7vFLYzNvHhY1/j6itqFhGX0QNf6GvHlTqqlvzqz9+jmkIfc1uCF3/H1kb0P2+KF+G7bzoPeOymYOebMhz/vmQxrMm1lkjx/qt8H7x3zuBPwx8PDagUiSJEmSdMQU6121nzVNh0WNv4+obWiwD2pryNEH1gIrsz45tvZr7ZAG53KzyB4/5G9D7fjNmcrGPBlDwzLyJORvg/Gfz5qnfKx5kiRJkjRn7gufU4Qxjn0aljHOFzX+hvxjxV1Fjb9hOX0gjS3qa6Ih/+u6wfycRfY2ZI9/adzTXdtcx3FY1Pj7iNqGBvugtoYcfWBNqqIxbx4WNf4+orahwT6orSFHH1ifqkjmsDct5Hn9H9KQO/4+srche/xQvw3u6a7MzJsxNOSOv4/sbcgeP9Rvg/eP+bi+UZIkSZKUgXu6n2bt0XFR29BV1Pgb7IPaGvLMQTjvqprMk8dFbUNXUeNvsA9qa8iZJ6Xa3E+9vobc8feRPX7I34ba8bvGc5msV5IkSZI0d9a51OGc4HFR29BV1Pgb7IPaGnLOCToWquzcwyCGhtzx95G9Ddnjh/ptMI8uk/OKkiRJkqQMnB8+zbncw6LG30fENjQ4j1tbQ54+cF2LajNPHhY1/j6itqFhGX0QNf6GvHlSqsk9EGJoWEYOgfxtMP7zWa+0TNYrSZIkSZqzKWpcvhf4woI/X/1FGGfZp8Gxxtoa7IPaGpbTB9LYor4uGnK/rhtyx99H9jZkjx/m0YYluT/wUuCWA49/I/Bi4O3AW4BnA7dNE5omZs3rcVHb0FXU+Bvsg9oa8rxPtN5HkZg3j4vahq6ixt9gH9TWkDNvSrW5XqS+htzx95E9fsjfhtrxu7eAsltS3sweP8RtQ1fZ44f8bagdv/eOOblGUpIkSZKUgXu617NvnHmMmiaoP551SEPucemuosbf4Lh6bQ15+sBaYGXl53vH0JDnencu46+vdhusgVI2fe5Fa7++DmkwT2Zh/PXVboM1TzlZ8yRJkiRpztwXPp/a4xuHNDjOV1tD7rHiPqK2oWE5fSCNLerroiH367rB/JxF9vhhHm1YEvd010bJvWkh7rWhIXeO7ypq/A2+R6mtIU8fuI5DkSwxbzbkuV6cK2r8DfZBbQ15+sAaVUWypL1pIX8bosbfVfb4IX8basfvekZlZ96MocG8mUX2NtSO33vHnFzfKEmSJEnq6gK4CtwOcIeBP+TzgQ8DfnScmG7wSbQBPqvQz9cwbwc+EvgV4E3r/4/tEbSDY68CvrrAzy/ZhvsBPw+8grYQ/9Ej/uyNKfr