{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Optimal lumped mass matrix for the FEM"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from __future__ import division\n",
"from sympy import *\n",
"from sympy import symbols\n",
"from sympy.matrices import *\n",
"import numpy as np\n",
"\n",
"r, s, t, lamda = symbols('r s t lambda') \n",
"x0, x1, x2, x3, x4, x5, x6, x7 = symbols('x0 x1 x2 x3 x4 x5 x6 x7')\n",
"y0, y1, y2, y3, y4, y5, y6, y7 = symbols('y0 y1 y2 y3 y4 y5 y6 y7')\n",
"init_printing()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def mat_fun(M, fun):\n",
" \"\"\"Compute the matrix-function for M\n",
" \n",
" Parameters\n",
" ----------\n",
" M : (n,n) matrix\n",
" (Invertible) Square matrix.\n",
" fun : Python function\n",
" Function to be applied to the matrix.\n",
" \n",
" >>> M = Matrix([\n",
" ...[4, 0, 2, 0, 1, 0, 2, 0],\n",
" ...[0, 4, 0, 2, 0, 1, 0, 2],\n",
" ...[2, 0, 4, 0, 2, 0, 1, 0],\n",
" ...[0, 2, 0, 4, 0, 2, 0, 1],\n",
" ...[1, 0, 2, 0, 4, 0, 2, 0],\n",
" ...[0, 1, 0, 2, 0, 4, 0, 2],\n",
" ...[2, 0, 1, 0, 2, 0, 4, 0],\n",
" ...[0, 2, 0, 1, 0, 2, 0, 4]])\n",
" >>> print mat_fun(M, exp) - exp(M)\n",
" Matrix([\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0],\n",
" [0, 0, 0, 0, 0, 0, 0, 0]])\n",
" \"\"\"\n",
" sol = M.eigenvects()\n",
" n = M.shape[0]\n",
" vals = zeros(n)\n",
" vecs = zeros(n)\n",
"\n",
" cont = 0\n",
" for i in range(len(sol)):\n",
" for k in range(sol[i][1]):\n",
" vals[cont, cont] = fun(sol[i][0])\n",
" vecs[:, cont] = sol[i][2][k]\n",
" cont = cont + 1\n",
" \n",
" return vecs*vals*vecs.inv()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mat_norm(A):\n",
" r\"\"\"Compute the norm of a matrix\n",
" \n",
" The norm is given by\n",
" .. math::\n",
" \\Vert A\\Vert = [\\text{tr}{A^T A}]^{1/2}\n",
" \n",
" Parameters\n",
" ----------\n",
" A : (n,n) Matrix\n",
" Real matrix.\n",
"\n",
" Returns\n",
" -------\n",
" norm : nonnegative\n",
" Norm of the matrix.\n",
"\n",
" \"\"\"\n",
" \n",
" norm = sqrt((A.T*A).trace())\n",
" return simplify(norm)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def mat_dist(A, B, dist=\"frob\"):\n",
" r\"\"\"Compute the distant function for the tensor `A` and `B`\n",
"\n",
" The distance functions are\n",
" \n",
" .. math::\n",
" \\begin{align}\n",
" &d_F(A,B) = \\Vert A - B\\Vert\\\\\n",
" &d_L(A,B) = \\Vert \\log{A} - \\log{B}\\Vert\\\\\n",
" &d_R(A,B) = \\Vert \\log{A^{-1/2}BA^{1/2}}\\Vert\n",
" \\end{align}\n",
"\n",
" where :math:`\\Vert M\\Vert = [\\text{tr}{M^T M}]^{1/2}`. \n",
"\n",
" References\n",
" ----------\n",
" .. [1] Norris, Andrew. \"The isotropic material closest to a given\n",
" anisotropic material.\" Journal of Mechanics of Materials and\n",
" Structures 1.2 (2006): 223-238.\n",
"\n",
" .. [2] Moakher, Maher, and Andrew N. Norris. \"The closest elastic\n",
" tensor of arbitrary symmetry to an elasticity tensor of lower\n",
" symmetry.\" Journal of Elasticity 85.3 (2006): 215-263.\n",
"\n",
" \"\"\"\n",
"\n",
" if dist==\"frob\":\n",
" C = A - B\n",
"\n",
" if dist==\"riemman\": # This is too slow\n",
" C = B*mat_fun(A, sqrt)\n",
" C = mat_fun(A, lambda x: S(1)/sqrt(x))*C\n",
" C = mat_fun(C, log)\n",
" \n",
" if dist==\"log\":\n",
" C = mat_fun(A, log) - mat_fun(B, log)\n",
" \n",
" return mat_norm(C)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The mass matrix is computed below"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAABmCAMAAADfytJlAAAAV1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAABcPecEAAAAHHRSTlMAMquZdlQQQO0wRLuJZiLvzd3B5am3atHZ83xswq/p\nLQAAAAlwSFlzAAAOxAAADsQBlSsOGwAABA5JREFUaAXtm9uaqjAMRosUdAviuM8H3v85dwsCzd9M\nU/xSr+BihE5Mli0Cy4IxsHSwDZt93dUttMFmJaQwUo5mTFboB2NOH1CUbFb16UEaoo04RzVOy2kO\nbfs0w8MxVKP7k1gagWHN8ZhL+4y1dUszJ7VtmuF+NWYY3Z/EIjGsOW6+8Gn0DNWWr2oEBh8qDJeR\nGGiOGzL0JoPhYTdobi2HYcuBDLchg6GvucJBWwZDkAMYWvcBxX64SQgZYxHmAIamdstY34LPFK1e\ne2OuwR4UBbj9RfheGJIDGHy+If29qD6apumShxDT3BmwoInmiBlu3Xh6fk2Dd22r9+krvW3Ha0P9\nGE/JvZbmiBninKVbDoa5h6Ef5sO3+/tp/xeIAIZPS5f8R8TQ1L37/qcW9QhksO4MbpPHQf0IZPAn\n5fRBSj8CGIbpNJ66OigQAQztzJA4ThaIAAZzd+WrMXXO0o9AhsZdEdvUWBj9CGQwlbXCFat6BDL4\nc/I1aQf6EcjQfa27b+fUIapzZ+VL8hImJ0fVdWsOZPhe9/WPpMFc+75LX+XIOX7++v1nMxRkWO0j\n0RXaFoQMq30kGLQtCBl8acFg1C2IY9jsg+0MdQtiGAL74Bj0LShmCO2DYShgQREDsQ+GoYAFIQO1\nD4bBN6UvMExOjtCCkIHaB8+gbUHIwFct23owzP3r++E8rifKAgbDjCNU+Yu/BTFvKd0U7Q/qBsN+\nAlIFGfQNhmOgVZBB32A4BloFGAoYDMMAVYChgMEwDFAFGA7HWbpM3WCWxOSVVMGx0DcYUvu5Qasg\nQ2dD++Den2Mw+yzIM7TWc80LtY+lNXyVDWavBQ3O3Mj8hfhbs74F4VjIv3frW9B+Bj8uuhb0GoOu\nBb3EoGxBrzBoW9ALDOoWxDDsmoMJjxzB+i4LQoa9czBB3WB1nwUhQ5DobasHw9zVvh/Oly9Lv4N9\nLM3Ba4GIfxc8ZwX13rQa7Q/EPlgI9QhkoPbBMehHIAO1D45BPwIYwD4YhgIRwAD2wTAUiACGw3GW\nbif2sTSSV/UIHAtqH6T2c0M/AhlyDEbbgpBBNhh9C0KGHIOR7rjJyRHOBSFDjsFIDDk5wrkgZPA7\nnmAwog1m5CBzQRyDYDBZDEIOMhfEMAgGIxup6wchB50LihkEg/FjJd0FZoQcMBcUMUgGk8Mg5YC5\nIGTYOwfj9z9ccnKEc0HIIM/j6FsQMuBnesf2wTD3MvRDAYNhRhOqAAPzhvJNEYO6wbCfgVRBBn2D\n4RhoFWTQNxiOgVYBhgIGwzBAFWAoYDAMA1QBhsNxli5TN5glMXklVXAs9A2G1H5u0Cozw3TwnJ5D\n8XeibXeRce9XjVifQ2n9wyDWTo+WuDmYeruLjGNQjZieQ7HW/Af65m5zAJ+iKQAAAABJRU5ErkJg\ngg==\n",
"text/latex": [
"$$\\left[\\begin{matrix}\\frac{4}{9} & \\frac{2}{9} & \\frac{1}{9} & \\frac{2}{9}\\\\\\frac{2}{9} & \\frac{4}{9} & \\frac{2}{9} & \\frac{1}{9}\\\\\\frac{1}{9} & \\frac{2}{9} & \\frac{4}{9} & \\frac{2}{9}\\\\\\frac{2}{9} & \\frac{1}{9} & \\frac{2}{9} & \\frac{4}{9}\\end{matrix}\\right]$$"
],
"text/plain": [
"⎡4/9 2/9 1/9 2/9⎤\n",
"⎢ ⎥\n",
"⎢2/9 4/9 2/9 1/9⎥\n",
"⎢ ⎥\n",
"⎢1/9 2/9 4/9 2/9⎥\n",
"⎢ ⎥\n",
"⎣2/9 1/9 2/9 4/9⎦"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"H = S(1)/4*Matrix([(1 + r)*(1 + s),\n",
" (1 - r)*(1 + s),\n",
" (1 - r)*(1 - s),\n",
" (1 + r)*(1 - s)])\n",
" \n",
"M_int = H*H.T\n",
"\n",
"M = zeros(4,4)\n",
"for i in range(4):\n",
" for j in range(4):\n",
" M[i,j] = integrate(M_int[i,j],(r,-1,1), (s,-1,1))\n",
"\n",
"M"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Ms = symbols('M0:%d'%4)\n",
"\n",
"def f(i,j):\n",
" if i == j:\n",
" return Ms[i]\n",
" else:\n",
" return 0\n",
" \n",
"M_diag = Matrix(4, 4, f)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Unconstrained optimization"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Frobenius distance"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dist_sq = mat_dist(M_diag, M)**2"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"grad = [diff(dist_sq, x) for x in Ms]"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAAAyBAMAAABlgJeHAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJlUzSJmiTKrRN3v\ndrsdCiq5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEEklEQVRoBdWYMYgcVRjH/ztzM7ubvcstUdJk\n5UYbEQMJ2FgksqCNhWSqFAZyU5giELnluoDhprCxuksZEG+7BFJEQQSR4GliYUAziBAL4VbtbDxN\nlDvxcr735s3Me5vZed8U9/byijfvve+33/vt7MybNwuwciTg9VNSzg8zUfdRmDWfgqP/uCstlwaV\nuq2gMjwW3H/63CfplK2dsanHup3R2EBld/9pT/q2/6wUwc+1vC3QJ2IhPN8Xh0mV985oUqhk3AZ9\nKxITL6yVzF8MNWdGRcfYskEvBkLD4H2hlrcNWgpXezeSOt5WaJJ3E3W8rdAk71dreVuhKd6NQR1v\nOzTF2712bf3bxLiMSMAOTfFmQrMjqjbnLNBE7/la3hZomnfzj52In0hasUHr3hf3TjK1uf9+IwlO\nk9a9sc63sx/8S9LGNGnd27vE9ykXt2je06Sl9yo/z4B/vs+q3oi3zWWatLz1l4ZC03E2gLcXE7Mz\nJ6ZJzzBRVr4QNZw5dokkS3HaAzprWavsOE3a/YcbNf9Ovc6623C768CDB2Kg10/Hy+ucbo0EQKSX\nf+1ynEhfuVxOr4csx3tDngm44D3EM/ge/sCtPNM63Xrp43Sgss5yN8N2vxIUwYx2ImdUSrfvAceu\nylCAXSTeLjoxlBfO+T1envx0RsNVvU10J2zJX5dnNNHOwN8qZlbpUwmOPS9DAb52Q38bCyHuFPSk\nVkbr3iZ6dugRng95bjj98ow3E+B4LGJehDeOor2B1S4+K6eV0ZwmeRc0XOW3VPKpTYXuxWogbzvs\nOoGzLfp+iOMJ380x79s5MamR0yTvgkaHnSlDKejn3ipHV2I+/qkIdoDNGKsJVkN8x0e0e77xSEB5\nldOZN5HGskhBpdP7UqPZ5w/tiiSbQ3Zo37+N03jt/ufdhRjiGtTX7y+7gpVVQWfeRNoPRAYijdZD\njus0cDi9WxcHIldWdYbqeoLGiz+ELHSIV2VFW0+M9HXcUJIYcvdiTzxh5CckzdahvhhZCGQgPcwF\n2vrdi7xvWMDXGKXjbigdGOjW66/8ouAGuhP6Xz1JA/p+MCeWl7t5G/gQ+Ih1jypDatN7eSdR+gZ6\ndm9vi063rlxWf2WZe6K3kpg12Q25wr5HpI9O6lmgad6Nx2ydGU7SHB+3Q0+4TnSZu+x8R/pQRc8K\nTfJeAc6oV3CFNAtZoUneM4lb43xboUneePenlZPVJ1mN2qBp3uzVna0n9LL/NNX7Tbo0I/efJnlf\n0rb6pm9ghSZ5/9h9NjLZFnErtNyYSP1icq3VfOF3rV/dsUIvBkJitl/tcuCi8o+qw+ZXp4Olnr7u\nwPvrYGmZbPgmSJRN+nPclNNG/Fy2Kfd3aj1YbLhVzOHezW2PrFVwBy30fsSN/geKBseO8yb+kgAA\nAABJRU5ErkJggg==\n",
"text/latex": [
"$$\\left \\{ M_{0} : \\frac{4}{9}, \\quad M_{1} : \\frac{4}{9}, \\quad M_{2} : \\frac{4}{9}, \\quad M_{3} : \\frac{4}{9}\\right \\}$$"
],
"text/plain": [
"{M₀: 4/9, M₁: 4/9, M₂: 4/9, M₃: 4/9}"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(grad, Ms)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Log-Euclidean distance"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dist_sq = mat_dist(M_diag, M, dist='log')**2"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grad = [diff(dist_sq, x) for x in Ms]"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAMIAAAAyBAMAAAD8XRkbAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhBEmau73WYi\nds1/9lIHAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADa0lEQVRYCbVYv28TMRR+aZpEd0mTokKBqaUL\ngqXpypIWlbkZGIoYekwsSA1Ty5QyI9QiJBiTDbZmRUg0YkFIDPkPmgXB2FK6lB/B9nvvYl9s6ZQa\nD2e/7773fbF9ZzsHc4MfoMpEhLWX6zqpLA6OYWb1DkbXvUiTSL6PjburS3CRsIkd2SjXKTSrr2ZI\nkYOs4A+csRI7LEntUsvq8PmU+XrtICMczBM1dsjNCiT3etfmcG/b5uAgE5ybSToEDYUc2hxgyuYA\nYCcTvJR0qKG0PWkch3AeLeJReoOxP4fsb9OhSLE/B6DfzH0I93z3AQ5w3NlhLfLusNFTkuzwGA0c\nj8c4Mw3lvuFwxb/DFE4t9WHyl3+H0h+9D+TneonGGiW4qjsE09iH7bm3HeqNVhWenj3TQm7ayRDD\nLxSRRimzx2ke68WuFCOHjYZHZZY6iGSLHDZ7DHus15QoOdSqHpVZCgeGHA6XGfZYV3akGDk0uzLw\nXCp7UvB/OmT6moN185T3z1MybZlNfdhFpYEsR9geXisKHsbUspNhCIfTkmk6jKicCzAcmvVzadmT\nw7bEqQ/4gtuJY6PGPBjvQ+6vVdQBb/StbAEazxIuIUzdt4+ZHQ7tZymhZbwPtQjVPy68FI18V1y0\nUlq4URVhvisuWkEyPNIg2SQYwFg1aOUrNaDWASgmcj5BVu6BCZjIIBL0wjCAsfLR6j11ApU2wH09\nQ7QPq3AmqgRM5FzCgWDBx+WUniWcFMjPgpy4RA48iNQsJ2AiJ3rGGsLB2IGCI4GoIkdptKhRGoXt\nZDXSgowLBfVhuNV/G9URSLlhhe1kINg4CUyeoMLkTavUbTqDmjYOMsPmaQb4RAbXlk0ZjCb4v5h5\n004mDfNEBk1OzVzillFfqBshBQ4ywkFbsWgeYDOSYakLAY2XukuXLwCtqg6otoMcw5WeorFDqMa/\ncgoBny81xUEdWh0txqaDHMOJ031RrS7iYBbiedbQEw/FftdAZOAgx/B3zOA+wGUZFztQw77hXbq+\nhyyecg3UQWY4Tz81dqjVZfrDrVeGCgb5rSeRBbaTWYP+3PIO5HynLMIpIfUFQHDjPpR+psxMSVNf\nACQ3doBbKVNT0niQNIfCTsrcdLR3TFsZfv15zpiPWv/6E3/BKkQ+pEljXT2aYo8YHP8DJrYKCVT6\nqMkAAAAASUVORK5CYII=\n",
"text/latex": [
"$$\\left [ \\left ( \\frac{1}{3}, \\quad \\frac{1}{3}, \\quad \\frac{1}{3}, \\quad \\frac{1}{3}\\right )\\right ]$$"
],
"text/plain": [
"[(1/3, 1/3, 1/3, 1/3)]"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(grad, Ms)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Riemmanian distance"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"dist_sq = mat_dist(M_diag, M, dist='riemman')**2"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grad = [diff(dist_sq, x) for x in Ms]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constrained optimization"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f = mat_dist(M_diag, M)**2 + lamda*(4 - sum(Ms))\n",
"\n",
"var = list(Ms)\n",
"var.append(lamda)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"grad = [diff(f, x) for x in var]"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false,
"scrolled": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAAyBAMAAAD2CCIlAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJlUzSJmiTKrRN3v\ndrsdCiq5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFC0lEQVRoBe1aT2hcRRz+3r68/ZttHlXqIcFd\nRayFoBHxZLVbKoIouHjIwUP3ecihEMmSs5JFBBHErAdFMLp7KXoQFPwDoZYGtaCHJosU4kU2tp5E\naLRVU2Mbf/Pnzc68zNvk4EL34RxmfvObb7+Z783Mb2b3LUDpYJnlCUneL6EQ91o1NJNQzv0kVczV\nh1/OiE/L7qX7mBDvDNmUstd5MdSZd4W0nMB4lamoNbmW3G9DLYkN3rnU8pFeQaZCFaS2WI6xCi+G\nO+v6ONCEy6co/yfXUhLTNvS6xprIX2UqvN8TpatURv4GU+SIgJGU+arVkf3rf118rd7yGcWNUj2Z\n63CsDJfHjYTtr8IG0jzOJ0xXcQWpCtswCdOFTzHRGHJdzgwT0EtH1pYD3PHjO9xjztepnSnyFv+5\n3ENz63ikTtV0i+/P7ir7hJbYpTqaBkQ714521KubutCqU9Nrf/fauXXSdivOsfMv22qbWH6pNl1U\nGwwtBYnYZOryZpuEPLVpwu88YtM1/R2h0i0Tyi/VpotqA6It2oYlO5e6Ftk80TBPViib2GC2lqwE\nMz8T4vlVDcZNOhx3pQHRevwmuKs34bjJC7lUU6kV4LlaJ4K16ipP0kSc3kVt0zUgWjwVGadefZpX\nzghXqkhLsDPXCAGFJrfsuroUYdLi6xtDLYsPWXUNhhZyoLJn0X+YHwvIyvwhqtPuFlyf9sz6OndM\nVHhh0+UFxxp4MEfP4Q2Ocb7hBWy6FO38JZ+h9kObWX+TQfvRstugBmF2L+WYpBfbwjHjXcVtWEW6\n7oqZEm6brnS1NuVMjW7g5XcFSOY2XSFtppqraOC+tPfifg1qfVx4fUWHROylNsb5oyF/GTfQoe1Y\naECPNbYBpDDWKYJtxUcMQpsuRVvNyoXBP9KX9jDmdF4bbfayfiDtqLTJP+hewPhdkqKMr9wq7ZlS\nFV9rrPYBjJYf5it8P7oE7Wjb00fSnxb3aCOwztcrPjtq4lL+LDDZ4K1egMcPIbeCRR9faHjbAKYx\ncjzAEqH21KVoIX9UkdT9aZ1vtRHYdOXP4wEdErEn2+pXqXQVkx3QniFdpFYl2wBmkLoA2or70KVo\ngQKtW5X60joXdahN10cNLE4psqiR4ZH6M+4uAF0Cd7BYBQ3aCFylCoeo7BMUt5Blx5eYLy1wOdcU\nihuKFpjnDi0extOKuBFPy3534jepEGJ2+jF/LN02eXNrZ/EoHltb9ksN8I0gzq/U0vbbKOiBBFja\n7rgd99dtijhyHYpThF+qz/l6Hz1apMu8YV+0cjJiaQ/QwF0ehgRE75LsK7xeqxvuQtuIh7zttIHQ\nKlKX5slXtYpufoAP9SqzY2hpd9WCCDaWNsQ5h79XPYs1UyqHbbwslo3zi/sCA6FVjmq2MNO7PMKR\nPfFQ+J5DIQJlmcbnmPVND+JoFWwi8M7LinmfV4j5+SipE6hG03j23HumAzgUdcj66M7OZqQplvbV\n9bci0FhahaPo/L6sxOhSUGXs+bAUEgg0ew/zP6Wltbcg52PfuvYY363Q7NykeN4WI0mSLtDVeyFI\noK4F4JlOAnWNdNxEzhde+GFhKoHzRZJaCYyHbKaeFNOVqN+xMdv73pqoOH/Rvz0I5ytB75eRuVv9\nf0i+Nx+tSJ1JKcRbZhwwv10NvzrxdTn8u8Pw6wkVlDaE1ZX3j9A/5KX3pTzH0telMeSC5PCPqjd2\nB5vJUMRVZJ9gxb8usKv5iKiyUgAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$\\left \\{ M_{0} : 1, \\quad M_{1} : 1, \\quad M_{2} : 1, \\quad M_{3} : 1, \\quad \\lambda : \\frac{10}{9}\\right \\}$$"
],
"text/plain": [
"{M₀: 1, M₁: 1, M₂: 1, M₃: 1, λ: 10/9}"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(grad, var)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Distorted element"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"coords = Matrix([\n",
" [x0,x1,x2,x3], \n",
" [y0,y1,y2,y3]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dNdr = Matrix(4,2, lambda i,j: diff(H[i], [r,s][j]))\n",
"jac = coords*dNdr\n",
"detjac = jac.det()\n",
"\n",
"area = integrate(detjac,(r,-1,1), (s,-1,1))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"M_int = H*H.T\n",
"\n",
"M = zeros(4,4)\n",
"for i in range(4):\n",
" for j in range(4):\n",
" M[i,j] = integrate(M_int[i,j]*detjac,(r,-1,1), (s,-1,1))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"f = mat_dist(M_diag, M)**2 + lamda*(area - sum(Ms))\n",
"\n",
"var = list(Ms)\n",
"var.append(lamda)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grad = [diff(f, x) for x in var]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the solution is given by"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAEIEAAAAyBAMAAAAnJmFzAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJlUzSJmiTKrRN3v\ndrsdCiq5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAgAElEQVR4Ae1df4zlV1U/783Oj92d2Z0QlBAo\nM5BiKWK7SDAxojvFjRAJMqmKiRB2JAFihOy4UZZGSkdCQgwJu1SoJpTsxIiaCOlqClqxsISSYIV2\naCAbEpsdUKMo0i39YUsp47n3fr/3x/mee+85M9+3rMn3+8d793u/n+85n3M+595335s37wLg8Yxl\n8zgcQwaGDAwZGDIwZGDIwJCBagbetNlCph5dbZvD85CBIQNDBoYMDBkYMjBkoJiBmR8uNtdPrBeB\nw8UhA0MGhgwMGRgyMGRgyEDIwBv+xrXnngx9TGufWWjYB+Zip2ty6Lm7rjl1quNvdx1iW2IgwA2G\nin3YHSdyl8IzuTOc6sQI9+VaH7j11CnhalOVjb556sSQRiVWRAzM5RnEFsRAaZBISSVdNoTogs6g\nBq2IKuKTa2o852zIpctaSC5oOEmzIS6ahAl/orClG+WawPdYs/u+OvdiPjrS+9rF926SLv5UChR7\nFgMvq+fpZuWw/2E+Da53+iFcQdiHEqq9NkH0zM7OTnmt05KoP4ttiYHwJpNF+1B3L0HIPWet6cTI\nmgkXTqAEZ8JpoaXKRu88dWJIoxIrIgZmEyi2IAZKg2xquMdC1omhQ8ujyqY6XFDVbLiNtMSKkPv4\nUxUnaTb6pCi3pRvlqsB1RdPJ9PzOzlank+u4f+dprrvbJwWKPYuBl9fz9Rs29sMr3RT4ntG3zi6C\nffBdhcYk0eNVANnLV4Fgc0lsSwx83rW4grAPde8ihNhz1ppOjKyZ6MILAfZHp/mmKhv981S5B2lU\nYkXEwGwCxRbEQGmQrob7LGSdGDq0OKpspsMFnedwH2mJFSH3sac6TtJs9ElRbEs3ynWB69CdTO+7\n6eudPrbjVR/D1yDJIQWKPYuBl9fzHWs2G0vl1+WLi4iyD5LcTQ49CzC7IaEgwIhtiYGwYD6DsA8C\n/wKI3HPemE6MvJ32yjrAu9p2+VmXjb556sSQRiVWRAzM5lBsQQyUBtnUcI+FrBNDh5ZHlU11uKCr\n2XBf2hIrkt6WOVNxkmajT4oKW6pRrgpcVzSdVMveFeFtZzq3ZjqkQLFnMfDyej6+bBPw/2QFgVw/\nmRFsN91iW0KgruRFjIWe87ZUYzZvJrpyYCs6KTR12eidp849SKOSF+CepevflTRIZe4KRdBc0hnU\noeXS9c2zYG/v4nvjk8pGjxTFlaoa5brAdWif3abR+8uzeKkh9iwGXt4VRLN0kK4grvrGb16ovoD4\nMpkAGrYBBGZHD37kkzfb4hid3II3LtKCcecyW4gVAn0RB/fOEfcYML1Q5FyYPpUYIbPP+imY/jxv\n8xBA4M5DbK8qG/3z1LkHaVTSWhAXTSGD26JaRwNSoDRIVe5C0eQLWWUwvJ0U1ZlUusAzX9w6nnuX\nbgKcpBoLiyZIkFdXXoCq2Ugnhg5Npdv/1gc2RNPabe94kWxUSoFiz2Lg5fWsW0GM1mfuPrgCFy5Q\nAZLztkx2hZ7bTmzRk9k1aMwWgc+C39m6E+CD+FePqfNwEW6+iVlEBFsnv8VcDq6lQF/Ezv3shQ8F\nG51WwDiKbDzBMxtBx2i3QyNGJNj6K2Dh0ux1z+8ahHcCOO4LN68zl30XyYYRo3AQnmw22tsbnkY1\ny5MlQtzzsbQWIURVAYoVEQMDB9IKFioDTgz00lXEILkroqPhmB9rxGAlxzp0JF1pwDU8zTiyRcMG\nRT2XDBK1klOhIqIBRzixxR45l2ocKBbHWjNN4Uzg1GVLMdiqzKRklJfRJPDydBqWnWFijbJSax5a\nnP6+aFr7JhzfjGTL242AxQRHntmS9B4iYHn8RJ6LCY4NFgu9CNStIGYXZx+eWp1Znzrj42IabZns\nBj137TnGYuhawjWBJVEG/i7csvhyeO9HAX51vALfHa+Nt4ONtuVtza7uX2k7uWcp0Je8c/8TcB1n\nrekLGEuRj8d75iMomG8vacQIgk0fuRoObf8+fLg1E55HlwAc97fBz4bubotkw4hROFKefDba2x1P\no5rjyRIh7tlYWoP4qYqPqgIUKyIGBhKk5S3UBpwUGIKsiEFyV0RHwzE/1ojBSo516BBVccA5nmYc\nuaJhgyKeiwaJWsmpUBHRgCOc2GIPzkM22AAD0FMsjzU31I1gVl2+FL2t2kyajvIKmgReEUOHDnnw\nrf+WTWtw6FyQzd/MNVpgOcF4Z+u5olgAVsaPp1hJcDBYyW0JqFtBjGDfNsDBDTBfGcwebZnsBg1T\n5RXE13CytyTKwEX4FUPv5wAW57fg6fH6DL5C0MPbOrg69xi9GJ9Lgb6Inftr4ERshbQDxlLk4/Ge\n+QiISe5UI0YQbARfgcNbn4LTqx2bCyuYVJvev4fDi53LoYNkw4hROFKefDba2xueqJrjyRIh7tlY\nWoP43cGVNqoKUKyIGBhIkJa3UBtwUmAIsiIGyV0RHQ3H/FgjBis5VqK9dMUB53iaceSKhg2KeC4a\nJGolp0JFZAOu/X62bNT5bLABBpaeYnmsOadGMKsuX4reVm0mTUd5Ba0TQ4cOefCt61dFCYYxzjqS\nl0EPLCcY/beeK4oFYGX8eM+VBAeD1UJvKXaBuhUE4MsKwNIqfMHnnWm0ZbIbdC3bnzX+DIka0P7X\nrnnRumNx7gmA8Qo2yeFtzW/iR1iFQwr0RQzg/mn46oLRCGMpsksn75mPoGjeXdSJ4eU98DCc3rgX\nP7HruDi4bbowvtFjMH+kczl00GzIVhAydQ1Pq5rh+YcsEeKejSWQ9VFBBShWRAwMJEjLW6gNOCkw\nBFmZq0juKugwHLNjjRis5FiHjqKC4oCzPGG8Ara42aCIZygbJHpFp1JFwnyaH3App9qoi7JRHnCe\nomQmtYIZdflS9LZqM2k6G1XQaeCY3aK6OnSklW1+GOCWTcG0hh9Wz+JrhZ8nqR1/HgHLCQ6e7Tte\nb6DTiIDF8RN5Lic4MljObRGoXUGc3oDF04vw6U58UUcoEz26nG0Y2d+TMmYrQFw1jJwiV8PCeYDn\nbkQEXTOyBVOlD1XEwFDEzv3oSx2nUUfAWIrcCiL2zEQQGcs2dWJ4wcYr8C9o8+hix7CtGMN97mE4\nuNa5HDpoNsoTWsqzoi6ubrACjGpj5MkToe65WAJZHxV2MUF7oFgRMdCbpo1goTLgxMAoyLIYNHdl\ndBiO2bFGDRZzHP6kDW6MlNFRVOUB54oGx9HYFjcXFOVZNkgV8+diRWwhW3kdJy7SlBNf7N4zRNng\nAvTAQLEy1oIERl22FCNb5Zk0fF86jGDPqNNIA8e3LMXpVIemzu7B9+OLttzKCZ5fgTG+Vvh5ktrx\n5xGwnODgubKCiIFcnbSuY89FOWKDxdwWgbbg8GOFM61/9rmd3qcWL8KCWUHczcKazr2gy9mGue/h\n64YlUQa+dOESfNUp8m04dAauem2XcGQLDuJ72uwhBvoidu5HXytZjTCWIreCiDxzEWQZRxc0YkTy\nYs6+jVYYnU8vA1juZqiVAiTZYN/zBaIpz7K6TQUY94YnT4S652IJ7n1UfNAeKFZEDPSmaSNYqAw4\nMTAKsvjqEl7Co3FE6bXn8XDMjjWdGDp0iKo44BqeZhy54uZSQDwXDbYJYJ6liogGXMqJL/bAIWSj\nPOACxfJYayrAjB6jLluKka3ylOBXENEIDtRpKw0cKmLo0NTXMsCdkmltjH/EX4lko3b8eQQsJzh4\nrqwgYmBpMos9F+WIDJZzWwTqVhB3HHkIXm/+OP4VKBzti8Fu0OVswxSuIJzZIvDAI4cuTa85Ra6D\nO7DFfJMysgUnC+HETsvAtohb96VvUsYYS5FbQcQUmQhKpNtrGjEiwcbbBx7HL60ut2bC89FlcNzx\nuyPzkhVEGyk3YwezKc+ium0FGDEMT54IFYOLJbj3UfFBe6BYETHQm6aNYKEy4MTAKMiyGDR3JXQ8\nHLNjjRosi6FDR1GVBlzD08wErri5oKjnkkGqV3QuVUQ04FJOfLEH31E2uAA9MFAsjrV28BrBjLps\nKUa2yhOkX0G0YhTn3TRwZH6dZ880dGhq4EaY+aJkWptehrevRrJRO/48AhYTDMFzZQURA0vjJ/Zc\nlCM2WMxtEShaQVx7/11rYB6e/eBV71yDpQ3APwVljz2hy9mGafwChiNRBI4+9tULb0WCZgw9++br\nVwHmHunwjWzNLHeuRh1S4Pj2p24D89C6h9NHIjNpM8ZYitwKIvLMRZBaZM9UYkTyjv7kAfxy6ZsZ\nm288Ao776Ak4vMkAmq5ONooTGuFZVLetAKOa4ckS6bjnYgnsfVR80B4oVkQM9KZpI1ioDDgxMAqy\nKEYndyV0PBxzY61jsCiGDg1RVKUB1/A048gVNxNUx3PJINUrOpcqIhlwhBNb7JHrKBtMgAEYKBbH\nWjtNGcGMukvc3B/ZKs+kdJQX0SRwJFCYTu28K5580daP44tCdExfc+uiaFq78abnN/MPm4pgMgCL\nCYbgubKCiIHF8RN5LiY4NljMbRGYrCBmztrX2Yv35V/4TIYObpb/FyNkUY8uZzu2zL3ixtdN242h\nq/FbENP4jjp//CX8Vf5ifEUMxJvwT0vH1+KbmXaLubr2zdBaBD1K18g7xm+crR/YYEi3Xf9Q/l+M\nFtY+Fye0FtQ+S9R1YoxXQEKkEkvrF2pBe2BNETGwLp14wImB5U+4PfW2IZCuEaw+1oxNsRiWgBTd\nDqaWNffcqDZeaScHDtT2VQ32I514wFlekmJ3AQhUs0DBWGskuLo692smSAANuipGq5p9rqJnsm+C\npQmWDzZBgnWKSUeENMHVbLW5ZYDJCgL2/y9C585utjfwzwvL5d+DSO/SoWHqfHp79kwCfAVuA7Uy\n+j4cXJ3J/MCiNT/3iy/7ZtZPfEEMtDd9Ct6+GN/NtQ3GUqwEXokA4+xNOivYh2FpE973sl8oBfBA\n+fcgaLgohvwQqGvFsDwlRCqxBGZSYFWR1mQVWJVOPITEQFCJIUEbwSRjzWZFmmOXQjFaMOCsGLZo\nJEFVDfYinXjA2XRIit3lTaqxYKwZCay6lQrTTZA6dFUMF3XzWEV/NoFHJ9IEV1IRWRQk2KGliklH\nhDjB1Wy10XSBzQri9LqFvP6f8GnmbAvPPZ88WXphoXep0NM/+eQWNQDwY90uFjhaTYA3fu5j+FWS\nQ1swd/NN6ZUEBrhv6qW0B8+ILXtdDLTo91/4SMcq7TAYS5GNJ6ArEQD0KJ0R7G74H4D7d3YCg25r\n9hruJyu7ONdjxOgcjLAWw2cjRVsxLE8JETYWTmMWmHq2FKuKtMFWgXXpxENIClSJASy6Dc89W8Ek\nY83C2RxzYujQ7IBLpbNi2KIRBAWswTjwXUqXcgLpgLOeJcVugWyAxLMB8mPNmvAPRjCrruXquzsN\n3QSpQ1fFSNhU0faLV8ktzYk4wdLBxiaYq3ZWMY4iO34YacUJrmarZdEFHt62105s2qe3/Bs+/dZ9\ntn0lPfy8kMxUd/Vx4MEXCW8mMMYWQTSnYiB/O36mu1uKscGepfvkAxux9Um1pcI6/wx6bzzF0jGe\n+0tJz9Ltmtieg9xbIYvFsAGK0UxUeyuaOL+7lG6inGJ+nTbjuYPJdOxaXbFU1rEOneEq6D7efa0Q\n3NUbpP8wpdL27HnfeZuTz7jMLL8El6Qff7q3NPVl6EeRHHGixcC+ssHauVKlY8n6Tqmw7gYd2jvJ\nN8TS9e455nSlSDfRIOOA+bZYDHu7GD3RqHYp3UQ58dlten8UnsVS6YQtxlm/uO9cHTNBhC4pEiJS\naXv2PGW/YTj7mKO4fBF/33cGf8KxOe5qG9xzfoe7PaPpNnDZ5BDgnpIjtiUG2jzo0sSlLvQRz+EC\nwASli91027r4CDorrPOjQ3e5dXpIAvP1IvVMDHYc+o4CUChdwYL3YhtioDRIZ56gU5+7OSM882Jw\nUeXRhGelxATECc/4jl1Kt3dOJMiYU9omwKznQpCpQcEZsZWXytraE7rChtiO0c3L3MEzcWdok7yF\nC6RV8JAiCTCfFKlnIECptH17PruKgb5700Y7vXZ0A166/xI0+7+N7kmTkJ7ld7hLce7MbfZY3AIs\n3Ea3gcsmxwH9DmT55ATb2ZazNbfdAPK2HNBvlpkHWkv53SOzTLIXSLgRzkvX7O/Wp3SRG6apEpbu\npZkV1jlytv0OoRU0Q452iTUmUWU9O4PshoWpb4F0zY57OekI99R8fOaAvj7jS2mblGY2SHcXESM1\ntZszElBlHJGo8mhdVALizrOfZKI7/Khr9o2USlfJdOQh1yT1mYOBeMDlg8zbzl2ZkLDWnbMtGHAR\nmpMORuar5/hjxSv2qfNgE1zeLNPe4/gIBhshXqnf8u6bji2hmC0qQrHiubb3q/Htqq+Zl/ffC/Cc\nDzlKM6vHj4yOzG9DZf83iy7scOesJY8OXd0CrLmHbgOXTY4Fhh3I8slJyPAn1lbYRy1vywLDZpl5\noPFT2D2Sp1HqJeFGUC9ddX83z0mCjDxkmjphaTaywlpvDh02oiujMwSTbqnGNKqsZ2uQ37Awcex2\nNgyVGl300oVAo6uhSbiHC7RlgaE+6WV/vjcxvJldN0hA5XFEN9nNonVRSchbnmXpKhMbiRSy5SSh\nYzC0PrP3ibORDzJrO3uBhJuVyhmw6FCuErRkwFnjhajmvpzljxdc3gRTJGGft2mBgXg2TOe5MhcY\nN5RitqgIxYpnwWs/ofiKLXhO8536MRzeWgD8ikm62eEOd+C2E2b7xuYfy5/5GXO82gTGgXcadLOz\nVxVNt4Frk0Ntu503ww5kbXIOU6A5N+QKnp2t8HsUra1OQHS7zxbIOm02AmwSGrln0chw9BqTyTs3\nDFsaBQ3XYJrDS5f+pz+1YM9V0l1CBxEnSruJryssBVrPaTbQck7YmGfYiK6ErvB07nMaU7JpVAWe\njSJus9ooTbuSLgRqRLWMo4eUe1RJEcY2HbDZzDXilAkyjPU2vRRorTYpaTgW3Bt0r2I490lUmJ7c\nqFOVmIJnmGSMOM3hR126ZaFlHD2k0uG9baYjjG+KOLnJtz7qxNlwFJsgC0VjaZroCxVAw22l6pS0\nM0Z2dG3RfBU2PDsDrohmt17+4NFNE4c76O1N3uqvbo5Pu3NyNSlha9NcmNlxVqaIYeSKKqUoHTlR\nKBnPfrr66y2Al2zYRI5hfvnlcGKjtv+by3p+hzt3PX106PL2auGOeBu43zt27B+PHXtluBi3DLDZ\ngWz22LHX3HPs2HJ8WdW2Owu6X/8o27JAGK+g9TIQAXYjQH73SHNVeUThpne20kl+KMdx2kjXGqk5\nxZlOWINus1EU1jKwtpuN6OpoCWmpxnFURc/G4FJls1pLrC5dccc9tBFxL4dqgTBeKaPS0iwGaQ3F\nYtQsy65HAdXHUbzJbhGtKzEJ07p05b0NU+nqmRZwiuuzCBdnIxtk0Tx/cTLCOl/GtmjAWXg2qrkb\n5s/w5F2vS7DgR7lsrDBecbcVHiPi1frlt0UmxmOKxaKKKFY9t5Mz8ZWeWs/tdDXGv2LA2H158vWw\n74Y1uB3K25O1xsbZHe5aRPzs0PibVqKDbAPXLq+697rdANutNduFXRcn6bG2/O+HFWw5p+1mmQWg\n8TrGNMkSatCVIw03ArfSyVYQY5V0kRum6WxJhTXoKBt5Ya0nZ9tvRFdB21sqD1KNnWcfVd6zMYjb\nDX264te8xCycxzcA3U1gg3Q+UN5Yyp3H2F7nqq3PPHC8FzHyZsVX0oAq46hJYBtVHq2LSkJWIJ2v\nFN5eGql/u8iDJb1jO4IrXo0hAxQNuGyQxorySMPNS+XMOs91YQNaNOAsPBvVJ1bH5tOe7DG2CRas\nIFL2WXtuAgjE80lxnitzgfHjgJ5ieY4K21OXPUe1Ynzwx9gmp6F4y4YB3Wke4C0wxi2z7rMrCPxk\nonLkd7jjbrTo8hZg0W1kG7h8cuxugH4HsnxyItvZprUlWUFYoN8ss+LUBG5kqSc0SyxcSMMN/V46\n2QpCJ13khmnqhCXZyAtrPTmefiO6Cpoh1+mSakyiyns2BnFeuLvjiXYIpKtYSblT+9G5Bfr6jC6Q\n5p7EILZ2c5oGVBlHdmvIsMluHq2LSsK7Ll1tYksj7WEFQeozH4U4G9kg87azV9Jw81I5Axbty1WC\nFg04azwX1YEvVTYasgmWfEybss+mhG5tmg/Tea7MBcYPoVieo3YxcvKxxBQPPG1xFzfN0zlYeALm\nnoZ2/7fcV4vtLYUd7tz19NHth1fcAiy6gWwDl0+O2w2w3eItL0tkO9u0tiQrCOd0vO0sVZzmd4/M\nEslf8OEuraSgVrp2BdGndKkneqYTlmQjL6x142z7jegqaMqMO5dqTKLKezYG2Q0LqXeBdC7QrHQp\nd2o/Oneu2vqMLpDmnsQgtnZzmgZUGUd2a0h847XtPOXRuqgkvAXSuYlNKl2+nCR0DIbUZ/42cTZ8\nkHlb4iuTEda5N7ZFA87CfVRLKwn7O44A2N+1zv0vhkuwf4Of3JycOA9tWSaX0pOUeLl+KzsEO8OE\nYr6oUoplz+1rf8qdnDnPbro65D7KOb6OmNuf2pramvrOUx/y+78Vfw8iu8MdcedOHbq4BVh0H9kG\nLp8cuxug34Esn5zIdrZpbUlWEBbo17AVp9ndI7M8Chd8uAfTT8O9dO0KAnqUrsAHL+mEJdnIC2ud\nOtt+I7oKuszTXZVqTKLKezYGl7gNCykbgXRNoDnpUu7UfnRugb4+owukuScxiK3dnKYBVcYR2WQ3\nj9ZFJeEtkK6Z2ITS5ctJQsdgSH3mbxNnwweZtyW+MhlhnXtjWzTgLNxHRSZM/M1+OLqKDwfPOLP0\n0SVYsIKwHgSDjWxtWq7fyg7BjiyhmC+qlGLZs3/tpwmJz51nN10dXrFXlpZjALvZ4QduPXVqffoF\n390IyPGKf8EKna7FovHPdsfXKNKcM+hkG7goOR0obmvpt3ijybkBh9qDJzc5l6Zv7q5rTp2CG6/7\nzxYQb5EZ22Kchq0ZYyBj0/69Kr/fW4Ui49mF+/GWM31mir5jxHLKfkS3bxHgqgsvpob9OTU3XtEI\nO4730oyERfPUMHYhOmxEl6KhxpMVw3yXl1sldmohiSry3AEmlYqU26MTTVupOelCoI0JxlXg3rox\nz4yrUJ8l4Hj3YsRmsd2rGGivEzv2xRt/RqOuAx3no9o9TzvJZKTr1n+HU1x29K8Y9dx1BDYjo+sV\nc9RNXC4bHZttfRor9KhMU8VwI6nQbAeKfRlhs+jctpi87cqE2Yb6G6f+FZunLmw0HeOVzKsbBcbs\nm3vt0/Sp9xzBDxNuflfbGc8UcVIocBxPeu3N5pl6RmCYxaM5ihrMJpgC0SD/UsV4bqardG9OJGkP\nbnuyE/i/N2eeDbPfb0HFHe54dHdnL2eMouk2cH/R+gRIoRYYdiA7sBaA2HoTvlEfA/xx0hmdzGBI\nT47OwC2bptPaCltkxrYYp2HHxRiI+5IlNtFseffIGkXGswt3zVDmDmZ/t9QI3lSUbvohXEEsw36s\nfv5IzVlbUmFpNiJh0VdqGDus7bARXYqu8mTEyGtMdCNRRZ5TIK3UKGNpNFGlrkWguBkCbXoZV6E+\n4zsZV6E+80AbZBjrUZB4T2oTO4gYsVn85/Ra0ajEQNtp7NhhExiiikYdgRaj2jXP8qjr1H/KiZQd\npJkWcOLF6Hg1mqSeCwMutRnVp7FCjto0lTol4UZSodkUih15YbPo3LaYvO2ydG2k08twdA2euTle\nsT22itj9W1MgYW/vbR5uBPzF59Evm/cseFhgIB4nJQXmx1nq2UkbKEZFlRokFCuew3zQhGGeOM/N\ndMWvILjtyV6ISQD86Sn/QxzFHe54dHdnL8eSogvbwKVQC2R3IEPDz7sWVxA/A/Brzkn3cbwKcGb/\nIsyfN9esLXYfNWCc5nZcTG2i2eLukVWKjGcb7mjNUGYObn+31AjeVJJu9K2zizC1BtMrCGSP1Jy1\nJRW2nI3UMPq2tu/ndwit88TlYyIwGsxrTKD5qFLgLio1K10nUMbVHuszzbANkhvrRvYUennFQG9p\n7NhhM82OOgItRbX7oimPuk79p5zyZYeBSTjxYnS8ojWauHw2UpuW4uRn0g6/krDdaBp0bofQNOt4\ne4gqO+pMzsyx8AgcPgevgkPb9tTmjZtMCTBflnDxCDwF+9bhu9agBfLEU6D13JkLjA0xxdSgmGK+\nVjjPDUV+BWFDJg/rAO+CfwY4u9FcKO5wtzd0YRu41HABiCwXcAXxHvzJTRKJP53Fz5g25s/B/idM\nV8GWwmlqE80W01SlmPM844MQNFIjVU647+3sNozO5Syn5orxpdCaZx0aajxR3URgdJ/XmEDzUaXA\nvEFIowlAuXR7dtVRMOWUDxJvTKE16foVA72lsWNHSCCeJAeBFqPaG88JSGciqReyQgxxNlKb+fQa\nhr3NpCph0TGJplQGZXRVugP3wtLKTLNBVGnaToGFvH1ic/QonFg1+cNDDMzXb+q58MpyOT3LVxCY\ngy143WJYQbjEZB8nhpYbNnV/xzfhnVmOWCiw7+FmBVFA6WKX2jQO6xTl4RYCUBnBCe3ADxdn1/L2\nFOYUUPSnQtd5okWFGGKoFKiKhs92764UnBRQ+ypYKRqdGDq0NE0T4ckrp6g8SSFrxJBmQ26zPk0p\nwtVAdWWgRqfSHV07dD7tyZxJgfhXjO9kTKTdYuCV51mzgjhko/7cYhp87mxiaLlhU/f7d962lqOI\n/dvmGvmerumih9yp3KbxUaeo8UxZ+3OVEZzQ4OzjpXWXwpwCimxV6DpPtLhtciAQWAOV2lRFY3h2\nj23TJaAvdiUG/ojFwLClsaugfRcNOs8c26ZfIJ1d1VQGnE4MqWd5LdSnqcmoZVIojcZgtWh3T/v4\nSjj86pvz3yBvYQBS4L5lePT9D2yGG3MtMfDK86xZQdjXlANuS7NcKkL/xNByw6bu4fpHC0se90b7\njvVAO9OSOwWxTeOqTlHhOcMdu9bT4h8AAAVBSURBVFVGzCQ73sFvX2cPhTkFtH+e+CnomglCILAC\nKrapi90QpUf/rhScFFD3KlguGkWGbRrEsasM14tbZY4KFs4V9CWcFGKIPctt1qcpsVNtehWGMfk6\ndFALW3N/sAxLvwTzR5Je5kQMfP/fwujRNfgCYyPtkgKvRM+KFcTokgn74HYafO5sYmiFYVP3My96\n3b05jgBL9pL9aZE8CK8onIptWodVihrP2RB0RsyE9ub3/WCzD3M6zzp0lScGIBVYA5Xa1EXDZrt3\nVwpOCihy71kMjR4qlfvnyQrXNyeNGL0XjXuj089MqhN2wmgi3bs3lh6G8edJL3MqBc5+abSzCC/B\n+bRyiIFXnmfFCmJhxaThukou2ssTQysMm5fn34aZH+QV/JqhO95uSWefFU5BatM6q1LUeO6Fvn0x\nmDkDJ/IjScFJAUX2OjS+GJR5okWFGGKoFKiLhhWvd1cKTgoocu9ZjAlJNwmerHCqyhPkTiNG70Xj\nVhD9zKQ6YSeMJtLNfx7/G2PG/1wBuRqdioFfXnwK/y96M7o105QCrzzPihWE/fRh4UwmBbR7YmiF\nYfPyjK+Dp49Qcv7cfvrwDX+abSicuh9LFdi0zqoUNZ57oW8n2fnF0q+sKTgpoMheh8aJt8wTLUoF\n1kClNnXRsOL17krBSQFF7j2LodFDpXL/PFnh+uakEaP3onEriH5mUp2wE0bH0k2vwqFH5ldg5vG4\nl2mLgb8O8NCRe3AFscZYibukwCvSs2IFYaE/DRiG5JgYWmEYX57Np3/713KER0/ilallQAnLh8Kp\n2Kb1WKWo8JwPQWcEJ1nzQehtWXsKcwooutOhqzzxr09SgRVQsU1dNFy2+3el4KSAIvd+xUCD4thV\n0P55crrp6Es4KcQQJ05hszpNiZ3q1JosOpHu8MNw6HH8X4zqZxBi4M4iPLT2HcFnEFLgFenZVlEz\ndScJ7Z6cXsZvmyzDgmwFMTG0wrB5g/9F/BpzlvHc9zDM9wJ8vRtt2qNwKrZpPVQpKjynlOMznRF8\nMZjHlP17bCFpK8wpoOhCh67yxHqVCqyAim3qokky3Jz070rBSQFFuv2KgQbFsaug/fPkdNPRl3BS\niCFOnMJmdZoSO9WpNVl0It38Nhx8bAq/B3E+6e6eiIH3AXxu9Sh+D6JrI+2RAq9Iz0vLNphmIZEG\nRs6OIvQNp05+lHRnTieGVhg2df/ni/CODEX8+AFfYKb/7tS157KI5oLCqdimNV2lqPCcD0JnBF8M\nDtwAM1tZewpzCii606GrPBUCK6BigXXRcNnu35WCkwKK3PsVAw2KY1dB++fJ6aajL+GkEEOcOIXN\n6jQldqpTa7LoRLqZNTi6Bf8Bzz2SdHdPxMC3wtTTcGh9VP1fDCnwivR8fNmmaH6lmyna80ZM7et2\ndvCrIZJjYmi54fHtT90Gc7fmd9aCaRR3H25kUV1ByJ3KbZos1ikqPOdlURm59v671uCPrin8X7TC\nnAKK9FXoOk9cHkoFVkDFNlXRsOL170rBSQEF6FkMzIY4dhW0f56scL1zUoghTpzcZn2aEjtVZUZX\nBmp0Kt2fveC/8L9BX/CnaS9zJgUeeMF1mwAnr91gbCRdYuCV6Pn0uo3lEL5bH44hA0MGhgwMGRgy\nMGRgyIAwA7dsWOC0+XvxcAwZGDIwZGDIwJCBIQNDBkQZGP2wgV3cEuEH0JCBIQNDBoYMDBkYMjBk\nAL8Z2f5s0MyTi0M+hgwMGRgyMGRgyMCQgSEDogxM3ePXDc84I7pjAA0ZGDIwZGDIwJCBIQNDBt63\nZnLwf12kZo0pyn4vAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left \\{ M_{0} : \\frac{11 x_{0}}{72} y_{1} - \\frac{11 x_{0}}{72} y_{3} - \\frac{11 x_{1}}{72} y_{0} + \\frac{7 x_{1}}{72} y_{2} + \\frac{x_{1} y_{3}}{18} - \\frac{7 x_{2}}{72} y_{1} + \\frac{7 x_{2}}{72} y_{3} + \\frac{11 x_{3}}{72} y_{0} - \\frac{x_{3} y_{1}}{18} - \\frac{7 x_{3}}{72} y_{2}, \\quad M_{1} : \\frac{11 x_{0}}{72} y_{1} - \\frac{x_{0} y_{2}}{18} - \\frac{7 x_{0}}{72} y_{3} - \\frac{11 x_{1}}{72} y_{0} + \\frac{11 x_{1}}{72} y_{2} + \\frac{x_{2} y_{0}}{18} - \\frac{11 x_{2}}{72} y_{1} + \\frac{7 x_{2}}{72} y_{3} + \\frac{7 x_{3}}{72} y_{0} - \\frac{7 x_{3}}{72} y_{2}, \\quad M_{2} : \\frac{7 x_{0}}{72} y_{1} - \\frac{7 x_{0}}{72} y_{3} - \\frac{7 x_{1}}{72} y_{0} + \\frac{11 x_{1}}{72} y_{2} - \\frac{x_{1} y_{3}}{18} - \\frac{11 x_{2}}{72} y_{1} + \\frac{11 x_{2}}{72} y_{3} + \\frac{7 x_{3}}{72} y_{0} + \\frac{x_{3} y_{1}}{18} - \\frac{11 x_{3}}{72} y_{2}, \\quad M_{3} : \\frac{7 x_{0}}{72} y_{1} + \\frac{x_{0} y_{2}}{18} - \\frac{11 x_{0}}{72} y_{3} - \\frac{7 x_{1}}{72} y_{0} + \\frac{7 x_{1}}{72} y_{2} - \\frac{x_{2} y_{0}}{18} - \\frac{7 x_{2}}{72} y_{1} + \\frac{11 x_{2}}{72} y_{3} + \\frac{11 x_{3}}{72} y_{0} - \\frac{11 x_{3}}{72} y_{2}, \\quad \\lambda : \\frac{5 x_{0}}{36} y_{1} - \\frac{5 x_{0}}{36} y_{3} - \\frac{5 x_{1}}{36} y_{0} + \\frac{5 x_{1}}{36} y_{2} - \\frac{5 x_{2}}{36} y_{1} + \\frac{5 x_{2}}{36} y_{3} + \\frac{5 x_{3}}{36} y_{0} - \\frac{5 x_{3}}{36} y_{2}\\right \\}$$"
],
"text/plain": [
"⎧ 11⋅x₀⋅y₁ 11⋅x₀⋅y₃ 11⋅x₁⋅y₀ 7⋅x₁⋅y₂ x₁⋅y₃ 7⋅x₂⋅y₁ 7⋅x₂⋅y₃ 11\n",
"⎨M₀: ──────── - ──────── - ──────── + ─────── + ───── - ─────── + ─────── + ──\n",
"⎩ 72 72 72 72 18 72 72 \n",
"\n",
"⋅x₃⋅y₀ x₃⋅y₁ 7⋅x₃⋅y₂ 11⋅x₀⋅y₁ x₀⋅y₂ 7⋅x₀⋅y₃ 11⋅x₁⋅y₀ 11⋅x₁⋅y₂\n",
"────── - ───── - ───────, M₁: ──────── - ───── - ─────── - ──────── + ────────\n",
" 72 18 72 72 18 72 72 72 \n",
"\n",
" x₂⋅y₀ 11⋅x₂⋅y₁ 7⋅x₂⋅y₃ 7⋅x₃⋅y₀ 7⋅x₃⋅y₂ 7⋅x₀⋅y₁ 7⋅x₀⋅y₃ 7⋅x\n",
" + ───── - ──────── + ─────── + ─────── - ───────, M₂: ─────── - ─────── - ───\n",
" 18 72 72 72 72 72 72 \n",
"\n",
"₁⋅y₀ 11⋅x₁⋅y₂ x₁⋅y₃ 11⋅x₂⋅y₁ 11⋅x₂⋅y₃ 7⋅x₃⋅y₀ x₃⋅y₁ 11⋅x₃⋅y₂ \n",
"──── + ──────── - ───── - ──────── + ──────── + ─────── + ───── - ────────, M₃\n",
"72 72 18 72 72 72 18 72 \n",
"\n",
" 7⋅x₀⋅y₁ x₀⋅y₂ 11⋅x₀⋅y₃ 7⋅x₁⋅y₀ 7⋅x₁⋅y₂ x₂⋅y₀ 7⋅x₂⋅y₁ 11⋅x₂⋅y₃ \n",
": ─────── + ───── - ──────── - ─────── + ─────── - ───── - ─────── + ──────── \n",
" 72 18 72 72 72 18 72 72 \n",
"\n",
" 11⋅x₃⋅y₀ 11⋅x₃⋅y₂ 5⋅x₀⋅y₁ 5⋅x₀⋅y₃ 5⋅x₁⋅y₀ 5⋅x₁⋅y₂ 5⋅x₂⋅y₁ 5⋅\n",
"+ ──────── - ────────, λ: ─────── - ─────── - ─────── + ─────── - ─────── + ──\n",
" 72 72 36 36 36 36 36 \n",
"\n",
"x₂⋅y₃ 5⋅x₃⋅y₀ 5⋅x₃⋅y₂⎫\n",
"───── + ─────── - ───────⎬\n",
" 36 36 36 ⎭"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol = solve(grad, var)\n",
"sol"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can check that we recover the solution for the undistorted case"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"coords_undist = {x0:-1, x1:1, x2:1, x3:-1, y0:-1, y1:-1, y2:1, y3:1}"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAAyBAMAAAD2CCIlAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJlUzSJmiTKrRN3v\ndrsdCiq5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFC0lEQVRoBe1aT2hcRRz+3r68/ZttHlXqIcFd\nRayFoBHxZLVbKoIouHjIwUP3ecihEMmSs5JFBBHErAdFMLp7KXoQFPwDoZYGtaCHJosU4kU2tp5E\naLRVU2Mbf/Pnzc68zNvk4EL34RxmfvObb7+Z783Mb2b3LUDpYJnlCUneL6EQ91o1NJNQzv0kVczV\nh1/OiE/L7qX7mBDvDNmUstd5MdSZd4W0nMB4lamoNbmW3G9DLYkN3rnU8pFeQaZCFaS2WI6xCi+G\nO+v6ONCEy6co/yfXUhLTNvS6xprIX2UqvN8TpatURv4GU+SIgJGU+arVkf3rf118rd7yGcWNUj2Z\n63CsDJfHjYTtr8IG0jzOJ0xXcQWpCtswCdOFTzHRGHJdzgwT0EtH1pYD3PHjO9xjztepnSnyFv+5\n3ENz63ikTtV0i+/P7ir7hJbYpTqaBkQ714521KubutCqU9Nrf/fauXXSdivOsfMv22qbWH6pNl1U\nGwwtBYnYZOryZpuEPLVpwu88YtM1/R2h0i0Tyi/VpotqA6It2oYlO5e6Ftk80TBPViib2GC2lqwE\nMz8T4vlVDcZNOhx3pQHRevwmuKs34bjJC7lUU6kV4LlaJ4K16ipP0kSc3kVt0zUgWjwVGadefZpX\nzghXqkhLsDPXCAGFJrfsuroUYdLi6xtDLYsPWXUNhhZyoLJn0X+YHwvIyvwhqtPuFlyf9sz6OndM\nVHhh0+UFxxp4MEfP4Q2Ocb7hBWy6FO38JZ+h9kObWX+TQfvRstugBmF2L+WYpBfbwjHjXcVtWEW6\n7oqZEm6brnS1NuVMjW7g5XcFSOY2XSFtppqraOC+tPfifg1qfVx4fUWHROylNsb5oyF/GTfQoe1Y\naECPNbYBpDDWKYJtxUcMQpsuRVvNyoXBP9KX9jDmdF4bbfayfiDtqLTJP+hewPhdkqKMr9wq7ZlS\nFV9rrPYBjJYf5it8P7oE7Wjb00fSnxb3aCOwztcrPjtq4lL+LDDZ4K1egMcPIbeCRR9faHjbAKYx\ncjzAEqH21KVoIX9UkdT9aZ1vtRHYdOXP4wEdErEn2+pXqXQVkx3QniFdpFYl2wBmkLoA2or70KVo\ngQKtW5X60joXdahN10cNLE4psqiR4ZH6M+4uAF0Cd7BYBQ3aCFylCoeo7BMUt5Blx5eYLy1wOdcU\nihuKFpjnDi0extOKuBFPy3534jepEGJ2+jF/LN02eXNrZ/EoHltb9ksN8I0gzq/U0vbbKOiBBFja\n7rgd99dtijhyHYpThF+qz/l6Hz1apMu8YV+0cjJiaQ/QwF0ehgRE75LsK7xeqxvuQtuIh7zttIHQ\nKlKX5slXtYpufoAP9SqzY2hpd9WCCDaWNsQ5h79XPYs1UyqHbbwslo3zi/sCA6FVjmq2MNO7PMKR\nPfFQ+J5DIQJlmcbnmPVND+JoFWwi8M7LinmfV4j5+SipE6hG03j23HumAzgUdcj66M7OZqQplvbV\n9bci0FhahaPo/L6sxOhSUGXs+bAUEgg0ew/zP6Wltbcg52PfuvYY363Q7NykeN4WI0mSLtDVeyFI\noK4F4JlOAnWNdNxEzhde+GFhKoHzRZJaCYyHbKaeFNOVqN+xMdv73pqoOH/Rvz0I5ytB75eRuVv9\nf0i+Nx+tSJ1JKcRbZhwwv10NvzrxdTn8u8Pw6wkVlDaE1ZX3j9A/5KX3pTzH0telMeSC5PCPqjd2\nB5vJUMRVZJ9gxb8usKv5iKiyUgAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$\\left \\{ M_{0} : 1, \\quad M_{1} : 1, \\quad M_{2} : 1, \\quad M_{3} : 1, \\quad \\lambda : \\frac{10}{9}\\right \\}$$"
],
"text/plain": [
"{M₀: 1, M₁: 1, M₂: 1, M₃: 1, λ: 10/9}"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol_undist = {key:sol[key].subs(coords_undist) for key in sol.keys() for val in sol.values()}\n",
"sol_undist"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And also for some distorted cases..."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coords_dist = {x0:-2, x1:1, x2:1, x3:-1, y0:-2, y1:-1, y2:1, y3:1}"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAc0AAAAyBAMAAADM2J1WAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJlUzSJmiTKrRN3v\ndrsdCiq5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAH9klEQVRoBd1ZXYwT1xX+7PHM2mt7d0QjpAqa\nnaKqQIrAbdU+tV2TVor6UK2FCpGKxE4jkQoJtKt9CEJFrIuQEC9dl7a8lNZW1SiVGmndKEmFqoSl\nBYkQAi4i2vAQrSGJVOWPJVD+FnDPvfPjmfFc3zuP7H2Ye+453zn3fuPxmXvnANRWWOy6TNuOhkdM\nu1XxxGXYG49Ml9XEZB96+prPq0itvZTgVuh7D5T6RIyatmz8L/AT82AjahCPn937gdjoWTL/yT7F\n5O1/dzTZe54lrv8yBu5jta2fiTPG67ZAux1vidOmaphu4GLnYZwxXqdbGLXjTQFtodNpsaHu8svd\nCBh7xN8A53Ec+FOPRahYKGFJaOwx5EwU5vDMHxM8MMWbGG72BIoqMvvedVSbqrwfLjvD+OvbQL16\nC5g24+0x2pcbKXJQbYUmcndRU4Uz3OA5jJSlDjkPMWtzaaTvHGMm6r96BCw0PDeFPslzm7mRmCet\nQOG59XmOW3zJ/XkS5KR5mn5Pm4PVLhlLDeei8jdwbM+6RC5Py9G55y9VOcolKOM5eAfTwFhLHtlD\nHH7Fk9T62UlcxXhDDcxQ2V9acvCQqd/nKEWe+TYyLS3R7zlwVr6MAOJNJg81AxqpuL8qhRDgEw5S\n5LmRwHvemy6pBPYw501PUujTbQZK/08B6kMKp3yxj7CpwoxqPIs1HqiuvvKtwPVSn9mjpitArsxe\n06pNr2DophT8O0oqDYZS4/ktUFjgx8xDrXVMXLfVoAylWdhaKCPd9zUeCjd8A0PynQglz038x1Hi\nmbVQrOyuZBM8VRcoRbNbo9gOAu+mgXxZEQ8U2sjL12MBr/KQSjy37536Ay6bT9jKq8Dz0JJs4v6x\nd32TdnK71W+NYWO0JV3PFhjOZlWJ51ins4SBNR9Lw3YBg2s2NrojmZTpdDpNbNn3VRkwYP9z7HpW\nhu+UvvaoyX3CPHd1SqQtPlA4CXgzJnNJhuZzGHWebhYusKVJmyHIZGGeqE9SpCMCbPwkyVySofmM\nuTvUZeuN+OmjWv4ejiqj+VbfXSPIrsVenFCTzCUZ2pl021vUG3XhCsKGBTM8dkfu7znDfkeKtqNM\nl9VtJiu2ZC7J0M4Sdn5I/c8pgSu18VYsbLjN1RMN3qXTc8DPBFCO6Lkkc0mGdiazNtA5+UXV7J1p\n9iyRKTJzXP1PfkW6uAi0JqrOiN5nNU8S9slcomicEAb2DdYCpUbjrj/u7+ICoyt3TokD7ht3m3YX\nmkl/hfl5HnZ12Y8uEnyXbJtD+rv46KlrJoOnaMMiabo9WsU3c4sYmKfvGlKXFMtaQM8y6hXS7m8w\nG7BTv4kv4QKMSU3+S4ZdsuubjqLv1ZtgoEK7WbVmVMZLqRJtf74OdpqQtez5eETuHLCK3yiyW3iI\nlv4Q+SoCO81heou7rTeG5wItyNNxEaPz4U2kF77TWezxSWO4VQTljLWYCBi5Sy8avx5tdFHBZXyv\nhVXePsTCv7QKPeEjFfy7i+4reS5hniIXD11ouOdfEbCrT6NgfQc8Z3ytqxVJ2c2FWrztby1gQ5Xb\ndBs/WoncHGZMvB6Pjmp9FyWeXTS0wAMTDRoab0Nms82+NiKlcHJ/uZKO+ZHJOU3PLdJONjMq2NBi\nJwHi+UZoMuHAd1Hi2UUjT7dXqe1E+h1QzkDqstxl8Cyy8afS6SqbzTm85OmTXhUzLcxUQLFjshai\nXyt9F49nKNGJ0Zhi8ePybfRrZRPFu8jy1yfPQ+EUHUHPloA3WeDQMmg8yAM4XyxzF9/A9/GDiyfM\nkSr4Hjf6FiKHkyYL47Wui8cz7CJEG5YTovdlSB/+gu34UktraZ8usVQ5QywQfuVG0J+RfbRCl/Ay\n6MvTInPF+CTvvEu+Ecy3pN1Mt/79qQZJgyxKXAvlW/CSTB/0S/hrKEjGBL4y/xTTvRgy+AP6d47b\n/gjgJRkh2gX6ZR73Q/yIFYhA2w8r/P7cQfeYjvu/JZARAgYG2lxgAF6SEaOzP/z21SBcv048LeRK\npLSDhoD8GnYTyGtOSUaMdnB+mSdyLvOiTE0FQuLJ9cTzu8BPybzSg0R6/Rv3WgEVPWbshS1CU31n\nMYBOXaub0GzoZXpubMS3w/O/DxhyvCQjRjtQv8wj4BkIyMQi8TxA+wkSbTaWt7dZSUYZTfnPxEAb\nqab4EYjMWeAlGfHj5cD9Mo86z9mreCEyV5/hGJVkqn3sURPxHHxkDthRvXDslGSEZt/glnnUeeY6\nv7B9bxXhZPDRlzkQT9RvJ7iRFDCSa+OnyFhcr84Tm24lWTioJJOgMZ7pDqX1BG12Ug72yjzqPI11\nY2zrpNzybWUoARnP5w49aCTx4RsCmYNb5lHnuQvGgyQ/qMohqrtI4mnUMHGqq5FK6bYUwgDn+aLV\nedISnP2IUnQUa2o4F0U8C6Zobxof6Uq8OqTd6pV5lHmmFqnSY4eC9B24JZm+mICReI7Q8FhAJRE1\nKslIIEDHK/Mo88QZ2iRWpIE9AC/JeAOFnv2eFP0jBagLOUglGSmajjlOmUed519M7JHG9QG8JOOP\n5ALxHNwMoyVHugidlWSkaL/M425sXboCv/TxpWPIHuX7eAEkquYlmahSPF5/8YSNI2v5Pl6MClqc\nkkxQEyf7ZZ5xi5sL5TjUMtLNTHIyQ7SBXdZtusrp6V8sa5ZIPXL5LahngMfxjmz3NiDGPfNxXL/i\nmrXTPrsVNUWfxxF2yGar/j+S/khhkS01WwAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$\\left \\{ M_{0} : \\frac{29}{18}, \\quad M_{1} : \\frac{3}{2}, \\quad M_{2} : \\frac{25}{18}, \\quad M_{3} : \\frac{3}{2}, \\quad \\lambda : \\frac{5}{3}\\right \\}$$"
],
"text/plain": [
"⎧ 29 25 ⎫\n",
"⎨M₀: ──, M₁: 3/2, M₂: ──, M₃: 3/2, λ: 5/3⎬\n",
"⎩ 18 18 ⎭"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sol_dist = {key:sol[key].subs(coords_dist) for key in sol.keys() for val in sol.values()}\n",
"sol_dist"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the sum of the elements is the same the area"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAAoAAAAOBAMAAADkjZCYAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXYyEM1EmbtmIu9U\n3auvYvmWAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAVElEQVQIHWNgYBBUNGBgcE1gD2BgCGTgaGDg\n/MUABBwbQCR/sFAZA4N/JAP3AQb/LwzsWxj4LzBwfmPgDmDg/AlSA2TzAMWBihsZvA8wMDCWP2YA\nAHWLEd5O6O0DAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$6$$"
],
"text/plain": [
"6"
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"area.subs(coords_dist)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The diagonal scaling method would give"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAKMAAAAyBAMAAADCeHZVAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhDdq5lmIkR2\nu82aysa+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC0klEQVRYCe2XO2wTQRCGfz/OxudXIl4SCDAn\nxKNAuKCgswUB0eEOUSWiooG4gaQiJxoiUBS3oQBDgSKq0CEaTIFAabAEokEoKekSEQziIZbZw+ec\nd2YbKxQoTOFkP//73d7cJXuHfeoTNq7iSg1h68jpjTPCGTk1hG1aeGbS80rc/MB7wyEKL5K3Bex6\nR8uEk13lPaVUncVcH5UGo4gp1eIU80h8IxwqbwFpnsqtoTjGceHQIofA6zK+EA+VJeAAj2V3Y3SZ\nY+HgOvR4wflFP0IlkG1pzEo6cYsS/ScO5JktAB8FnJ57VRUwUPAJr6/yiBRKHtQhs/I194fJ9Pjk\nJf3ZUzqreshrf5UzIu9EGtc3ZE+ZWxZDiG0X+WxTxMO1iDLVFjJuE/k1zq8DxxY4fgQslSPK0TrP\noNhBXt+9Ru0FZmsGo6GqYakRUZ7zeQaxNlLfOafoFU6xBxhvRpQVyrHKNFBpMUp/eZldnGIOiZ+E\ne5fnPnWB15Pp9xzCvTEpnDey0zMLlO4phZkDov/KARsnTNvsvSzSbqYUa0xAV01sCYd4s/ey2y5p\nh6WvZBzssGafadzfy1RHiAAWPC79O4oqkz7pnhrK+Zm3mpgYSZ9otkkfkXKnP1RpuL7KYPNpRBL0\nq1P/s8kYmNbdpm8z/WGcR1zvxevK55RyjLnpGmJDHAM6jAuG8iowHFW6dyllHjc2hvRnjhGEYRwf\nO4CpamSV8YI+sFGFTqA0KCCGsVLrV96RlOQSr7YtDH0ThL10WhblsxJbI6zhrH7qD5VxWJSXuRHW\ncHAjhMrjFuUWocOwhTGjj99VOiWL8qWwSGs4V48oExMTUztbfHrCBz08GWUL4wTouSxcJU2i5x9e\nZ4FFTi3hpI9cM6osCkr3ondzTFJKYTz0Dl+jcHh5EF/62mCTC7RTSEoxjBWloi8pzDY46K1ycIU5\n859R/oWX541/xf8N4RG6A7bFd6IAAAAASUVORK5CYII=\n",
"text/latex": [
"$$\\left [ \\frac{7}{4}, \\quad \\frac{3}{2}, \\quad \\frac{5}{4}, \\quad \\frac{3}{2}\\right ]$$"
],
"text/plain": [
"[7/4, 3/2, 5/4, 3/2]"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M2 = M.subs(coords_dist)\n",
"M2 = M2*area.subs(coords_dist)/M2.trace()\n",
"[M2[i,i] for i in range(4)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And the row summing gives"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAKMAAAAyBAMAAADCeHZVAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhDNdplEqyJm\nu91CXIoXAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC0ElEQVRYCe2XO2gUURSG/327O9lH8AEqatxK\nBDWFhV0WjWA5tQha2mUtBDGFi9UqLImdYGEaRVKYxSoQkK1sTLGgjYXuYiE2EjGu4gPHe2fuzeyc\ncwZhWQuJp8r95j//3Dk75P6DQ94njK+SnlfBztmz43NEZvZMBbu0YfZq+pZg/KT6QqBAuiFgp3ps\nWuG0sUx4Xp2rnAZmOhwDpa5AV5H6qbC1zB6dF0S5TZT7AscVyfLGNH4NWRakRuT34uIl4YrzTLK8\ns5j58kdLJRAfPJmVLBF98ML6tZqwHeCyRJ/HWGYbSm1nWXSdb0Jz+rgW0crUZcvT77XSWqo/39BO\nf324xnESsiWS+oUcslxo82YgsZvTk3GWmHRDy4fAiUXW7LRR3GQ0MyVbrgG96dByD7Cg7kCqPEBR\nv73RSi0vr+yrR5leeS56ndCyAbzTOFqJLko/oihYqQu8DgAb7dByFRP7uWiigxlhO0BZslxH6qvy\nsD+Pc/ORyy1xt/lKoEj2vnc4zzdbi4paSy4Ymfy3HHl0rHG7z7Ls6WJj8elHimPEFm/3WZpxiScs\n5IMX/glL56zW0VmWBoJEZQsZb0j/u0JLk2fuEcvV1ktNCDbifDuqdpqva4rYXZo804mKMkvBiUSw\nEU9ExTiPpD64rWWQZzKkt+AiUQEoNuHnArF8AEwOWQZ5ht430UfhM0CxCT/k/jgIrNTCXZo8Q+6b\nHfiWhKqln/oonnOJpZ9nqCru144RQ78EdpYI8gy3vD/FWaw4r1P/lmWQZ3j7B44U8cMPv+KH49DS\nzzNMtaPLkA90+OHV0shargV5hqmuM6JAnDi3pNXW0uQZjYYr1dD9tGLEOAWVy7YsTZ6hzeeAecqA\nGHG6gdyQpckzpN15W73dJ0wtZTGeVo88Vlftg5s8Q9qz6qToE6aWshhznjf8kcL7RiZ2lyMb8MZ/\nxvIvfDyP/xP/N49axR9f/zwJAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left [ \\frac{5}{3}, \\quad \\frac{3}{2}, \\quad \\frac{4}{3}, \\quad \\frac{3}{2}\\right ]$$"
],
"text/plain": [
"[5/3, 3/2, 4/3, 3/2]"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M2 = M.subs(coords_dist)\n",
"[sum(M2[i,:]) for i in range(4)]\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see that some values are the same..."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Eight nodes element"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"Haux = S(1)/2*Matrix([(1-r**2)*(1+s),(1-s**2)*(1-r),(1-r**2)*(1-s),(1-s**2)*(1+r)])\n",
"H = S(1)/4*Matrix(\n",
" [(1 + r)*(1 + s) - S(1)/2*Haux[0] - S(1)/2*Haux[3],\n",
" (1 - r)*(1 + s)- S(1)/2*Haux[0] - S(1)/2*Haux[1],\n",
" (1 - r)*(1 - s)- S(1)/2*Haux[1] - S(1)/2*Haux[2],\n",
" (1 + r)*(1 - s)- S(1)/2*Haux[2] -S(1)/2*Haux[3],\n",
" Haux[0], Haux[1], Haux[2], Haux[3]])\n",
" \n",
"M_int = H*H.T"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAADKCAMAAAC8A+bbAAAAqFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADaTSGwAAAAN3RSTlMAMquZdlQQ\nQO0wRIlmu83vIt290d+v+YHDwXzXoSDp8ZHT4celp8X3hdnzqeXbYMmxz0614+tsH9aEbwAAAAlw\nSFlzAAAOxAAADsQBlSsOGwAAGKpJREFUeAHtXet6I7lxJcXLOKKkkXZnN+Pdib2etZM4ceI4N77/\nm6UK6INGAygAB1+T4uQjf6i6gTqoG9kkGkfozSZ5nfanh6O2bZ+THvN0AOLG6rcg6oNGAsz0ftEB\nlyAXneUTqEKWtdLWbr92T5vNTkqx3T28pYMY5wMQHYmwINqDRgAzXE+a4RJk0l06hSpkSafQVvFr\ne3avhwm1O202T2c9OXQXhId4Y90WRH3ALzUyw7zJ1l+4BNnSl36oQnZASn69+TrIG/W828vrEI2z\ne9ETzgAPIS2IRwN+aSAepkftF4KGbCNCngiIH3Xp15NW4UE+CtvzNrF6/OhaGAMDEDHKWBD1MSOA\nJUEap3AJ0lCLm6EKGffVjqdwFip6bUoLcty/7J1Sv4EBCGlB1AeNBNgicPMEQUOainMHVCHnntpR\n2a9SQWQU/VIn378DENLCmF8aiPdNj9ovZBWyjQifcwLiRy34ZRRke9ZvFMrAAIS0IB4NGQGsI7VR\n0ETwUIXsMxTCWajnBXEXttez/PbtLsgAZHKCCGHQSARbBG6ewCVIU3HugCrk3GMfWX7lBdm8vcrP\n3o86Mzx8tAdc9AxAPL7bgqgPGplhC5fNE7gEaSrOHVCFnHsqR4ZfhYJsdQopRXndvZ0f/Ld7ZVzX\nNQBRHGFBtAeNAOYcbf6BS5BNwBwFAdFRDb8KBenw4a5ysQzcC3Kx1I4NfC/IWN4uhloWxN9MKfy1\n7ReUcUvGBm02GaymnGsDXkXlsKo6xsykjcpU0WBDpAdKQcbay4LEPffjd8lAVpDpnv5p97zTX757\n+cnVcCzc1z+d9jp76YBMI7rbAX3qt2pk2C8r9qQguKd/kp+9D3L79kmqccB9+XJhwn39l8fNq0xc\nOiDTQIezlLxP/VaNDPtlxp4UBLNznbVsz6+bN73rq3mzX1hv2EvdjjJt6YD4wY4nHbhP/VaNjPpl\nx24U5OOjTHjOj0epyWbzMV4pKVZG7+tPWr0QubQdpSD96m5R4xaNjARvx24URNMun6qt+3C86TdD\n7eXuy5wP+7187XRC5AN40IJ0q/tFjRs0MhJ8JfZKQd72SFf9/om/r/963m02j299EK3uST4doSB1\nC6J9q0bG/KrEbhfkJAnWrxG5GDXTJff1j+6G/fmxF/L0OhWkz4K4catGeL9qsZsFeZJ6SMrcl7p8\nodRfukxxVq3zUydEv/79d0ifBRn7Ro3wflVjtwryKL93H7ebN5fl2q+scF//Rb9ppCptiBb3sJPX\nWX41dqnfqpExv6qx5wX5qIS0774/HA7PH5738jlpzEPcff0f//70+bei+rLtgWhFxMru/A/yW2u3\nfW5YEGUsHjw+i7b8pOvwyxnR34q/7E9/+NIMIzbydX/69fsuI96vD8+PMs9jgt/+o/5cKsWeFETv\n6X/6YbP50d1nEaLcaXdyM/YpvILQ+/qffpJL/M+70w9fuiA6yu7p6fn8u8Nm+/2n3zcsiDYWD7af\nJ+22X95VMfJy2Oy+NsOIjOz+uDv9LG+vzuC/+5ME/cQEv929/NP5oRx7UhCNApMdqfkFiXIDVtS5\nfp+ctvyZDaGlJWkE8gTZMjBHkSMKBdHhPIMrV6+Yug6EL8gcTsX5rMtHkzWXG5AnyLJW3ApNyLmv\nXBD3ZTXXcVa3j64D4XyCt5NvOO2QHAJphWwbgCbkjCgVJDC4cvUZuDy6DkRt9vsED4NvaGhKFgGf\nIJsGQhQ5olQQGc8zuHL1iqnrQEIoFU/yLu9b3m63UAjkCdIeFj3QhES73P0uUEmld4CQdh3IWEG8\nb3PU7SMKgbRCtoeHJuSMyAviLp8XJ8oNWPFO5yHMwZSOIkOl7kIbj4BPkIVBkyZoQs7deUHCJEyu\n1xckymGqR1hxXnf7hBhnQ2hpSRoBnyBbBuaYc0ShIJiEEcSv60AkTsInZAW+4bwtSQR8gmwbgCZk\njCgUJO6+H187A/eCXDvjDXv3gjQSdO1uLcjxgPvrgbqVHthupZrh3IQEjfTARLiOVFscb714CInI\n1NFgegaFTALxKgSG9F/a0HeX75CB9JI1QvyiuXUjRnxuLHpZJXM8hEHwoTQQSUEGiF88t27AyJRw\nk15mF4SHMAg+lBYiLci0AwDFevPTTZ1O9XHrsN5AGdGU2/QysyA8hELwobQQSUE0MJr45QvCceto\nI+qZTS/T3uKLh9AIPpQaIi+Iu5VDEdLmGzLd3DreiKS7Qi8rFmMIQhvhQ6ki0oL4pQCO9TYXpJNb\nN2JEsluhl1kF4SEkgg+lgUgLIpHJUgDHegsFIbh1tBFZKaCpdQOQASN8KDVEoSC6FECx3lAQhltH\nG3HMepJaV2WkFT9UPMKvHVH5qlLrlgVxVzddC6FYb1NBerl1Y0ZIap1Ld5WRViwIieBDaSKWBQlr\nIUL82v6hl5AmN/VlttPPrfPrDb/s//V3R2GXbbb/1mbjuexR1Lop3zYjrVgQbSSodXO+url1Pnib\nWpcUJCwF/PPPf/4XvcPVJqTpTX2OW+eM/EmYp5//ctp9kH3r2kY0URy1ThGyZZ3NSHMKhT8MtS4Q\n+Pq5dRp8jVqnBflw/lBwDN8Mha6sCbMdghTCQ3iE+IkoIDPXCw2zoUJnuYmFwB3IedSfyiQHUch1\nZ1TpSGc7JIqH0AhEAVnyvNTmDZV6zDYGAncg50GTS9bcwRbEfVtxBeEhPAIRQ0YB1g4nQzWVtI+C\nwB3IebCVCuJnOzJsbmG2tTziITxi9qffL/EyGFp6XDsjIXAHch56pYLIgDLbkVduYbaVHfEQGgF/\nIDMfjAZvyOgsNzMQuAM5j7heQTy1LLcw28qOeAiNgD+QmQ9GgzdkdJabGQjcgZxHXKUg7vJ5aW7d\ngBEJExFDzpGbR5EhUyfpoCFwB3Ier1KQXpacDDZTy3Lm12xrccRDeIQYhD+QCx+Mk9mQoZA3sxC4\nAzmPaBWkxOGaUekR5pMEiofwiMCrI/yS0GAojbJyzkHgDmQ8sFWQWOd+fMUM3AtyxWT3mLoXpCdL\nV9TRgizvZZkUrrJXmToayurSCoVMmogrQUq+2T6hhwijZABwDGffy4LGXV41A8Yl67JcMTy657Lb\n1jUYaXaameCnUQhIw61yQS7LFeOpdeHpOpv+betajDSzIEzw0yAEpOVWsSAX5oph/qyzqT5qXfjv\nf4JbhyUKAuLSSwXvC8JAWm4VC3Jxrpi/Y8BR63gC3wDnTxJMB89DdOXE2huvVJDLc8XmWzjd1LqR\nHeXcLSaK8ycfWXKjOykhC6m6VSrIxbli4Zaf3AW73LZ1fomC4/xdgY3XcKtQkCtwxcInhKDW8QS+\nAc7fALVuACIrJyYVMS/INbhiKAhDravSy/yXa/ZXlygoDttA8AOQmlt5QS7OFZO0TQXppdb57w+O\nwBeWKDjOH7XRnXsDcPlqupUXRC6jzF5vYS2E2exNlwFO+y8f/122rTv2EOXcegO3bZ1forAZadnH\naWoYoNYx3DqsnDw+n/b/Ic8xSvfG04Ic9/inT+cTSUjDUkD/Zm9uGeCvT5uP7j5OF1FOjZDb1imk\nxkgrF2SEWrdhuHVI1+bz59Pui8zwE5bgq+wOmvzTJ6Yu4cpSdj1vxTdD3lNo4a3wCDELnyALniRN\n0IRMusuns2/l/kKrASldsgRNE9LUIhWCAngrNAI+QarZ+guakHXtqNf7FjW0DwuQckHcVw+dYTYE\n3gqPgE+Q7SRBE7KN8BqTb73qqleClAoSSF+kT5w6b4VHzO+pft+gCdmX4eBbn7orB55wu4CUCiIK\nNCFNB+VCGLJC+wWfINXP+guakHXtqJchyk2wAsQoiCd9kT6R6vprgn3AK42AT5BRAo1DaEIaanmz\n9y1vr7QUIHlB3IWNY71NJpkQeCs8Qt2CT5CV7CRR9COmLwOfs7YBrxGFs4DkBQkzPYmG4MrJqJQ6\nJkj9MB4R+dTvGzQhF+myTmbfLI2s3YAUCoKpS4nGlQ07N5DqvBUeMUCUQxSQc3zVI/hWVVp2GpBC\nQZa4+9l1M3AvyHXz3bR2L0gzRddV0ILMRDmwtjJZd4pTz7TRYBuBRibXhGSDo8E2oj3QCrKu3kLc\niXLN/F1XIb1kNWhcFefW44oVjdDb1ulmBg9uXaH7kbA8Aq4SwU8QC5EUpEXjggO5XJErlg8uOwDI\nU2zkddFHwl4l+Ck4M11pQUZ2lFMba3LFJp+Xwk+cdTZ1MW4dVihYah0VvI/KTldSEFXXm/QWjWuZ\npOiMpZfRRnxBLs6to/3SHLDB1xB5QdxNFpJdti5XLKrzfDjfWrokt+4qwVepdWlB/H19ll1GbvY2\nYmQuyOW4dSN+yVuGJBbWEWlBRFtu0ps0rvkduzjiuXW0kVCQy3LraL9GdrqrpatQEL1JT7HLRjZ7\no42gIBfm1tF+DQRfpdYtCxJu0lPsMnKztzEjU0EuyK0b84sMXi8sVWrdsiBhLcSicS0uVIsTYh+2\naSVg+zeZHHXvKKfLE9vn7Qv3SFhhoj1+3T0/dz1IFSsUDOXP5+BVHj3f/cxZBxFEeSs9LUhElAs3\n6Q0a16IG8QnNFftPmerJ80o7d5TT5Ynv/vrw5rl1Xc93RSinn3anT1KQlJEWe++Pgein/E1jSPCf\nfu15TG2wKYjvnotb6RWIcg6FSVIYonlAI/CtANm0MLAeK2NiLtkeftbodylg1oIkl6wwvh4UaFyL\n/vyEQiAEyHy4tAWakGl/8XyaSxb7rEbKgB9kLUilIO5LznK52M4hEAJkcchFIzQhF521E3fnqKaQ\n9NEG5Muw90HOs6kixCwIz/xiEfAHcvbUOoImpKWXtctcknrRBq5QEAmgQONqhEUhEDVkY2zphiZk\nG+E1dC5JvVgDMvhaEPMTIjYKNK5GWBQCIUA2xo6C7ke4Md1csj16pEEaUORaEKMgYZIUeVk/5BEI\nAbI+vvZCE7KNUI1pLtmn7LU4Aw6zFsQoSJgh9oeBaVU3Akw0yDYQmpBthGj4ueTiX5KaOMqAH20t\niFUQTJKavgcFEgEmGmQYxzyAJqSpuOzwc8llW/2MNKCDrQexClJ3+d57sQzcC3Kx1I4NfC/IWN4u\nhtKCzEQ5NRMIXzho24ZmkFVI0EoPbFSqGc5tiPYENRzU1YcQaxu5E+WaRbqugnHJsmhcFecIyAAh\nbQAy+ap+bR+2r0/yBJmOV2qoAyIq6xkpF8SkcdnuEZABQtoAZHLV+fUo16wH2/eoJzUUdVUOVzRS\nLIhN4zKdYiBYOSEIaQMQ76r36/Gw7ZwZpobMgOOONY0UC8ITv2iu2AAhbQAyMdI2j49x/lrHkaGW\nqu93+VrJSKkg7A5p4hULcXe+ODbeAAR+bR53+/2u8zOy2LquqyA++JWMlArCE79IrtgAIW0Aosn0\noWyOQgk+uF8dzQwvDTXVVzdSKEiNxmU4yEMGCGkDkMkv5/XxLFXpekWGevTXNZIXpErjKjs4AOEJ\naW55ZozAt3mWWhwdthzAsjXybdlRPEPwKxnJC1KlcRVdqjO/MkhYOeln4w1A1OwcihREmVDNV2ao\niVjbSF4QceGdN1UrJMEvtvDbw7k6bJ9PUou+dVxn6Mcvz//1W6mnezJjwZusyRHlSCMfj2AKLobT\ngkREOe17903VFg76E11s4beHEyw4bLJvW2HYvEkN/VkefPrw30Li6/hIuRFGjDyAKZgYKRLlsBwJ\nmftdaMGUqtBlNbEQ+ANpjZu0k+pD1LqwvJzYrp0W/SpesqAJWRt10UcR5TySgcAfyIVt+4RU34xQ\n626zIO4r0U5MqYeCILOQpQELbaS6G8HdnSqMZTYNGClCVvyEsEQ5iY2EIAJIMzvLDlLdgVlq3W1+\nQgaodRQEmYVc5t08I9V1nL6fZLHFASNFyIqfEHGPIsr5cBgIIoCME1I5JtVlJJ5ad3ufkDClqmQm\n6aIhyCxkMp51SqoPUesuXJBpJzmK/EUT5eIntlq5XLbDH8hlr3lGqg9R67jt9LyrRb8KlyyQviDN\nQJcdJFFOwRwE/kAuzZtnpLqMw1Pr7kQ5M/3ffEfhE/LNx/RNB3AvyI2VTwsyE+VAKMtk3WtSXQfj\nIJk2Gtb2C+MGWTfAR9JE3Ily7ZRfVSO9ZKU8sb009DlEEOUwIAEZ9CuFwXJVkqy3bOu66uDotIwk\nBUl5Yk9SjUMXxYwgysEnAjLoVwqD5aokWW/ZY2Grg6PTNJIW5Bo7yk1OXZhbp1aw4AJKHtJRkyzr\nbWUjSUFcFNfYUU4MsXS8iL/mGSTTtne17EpfBGto+u4R1tuKRvKCuDtMHIdN5tzv/MBSO9VxOLbW\n3DPCeqO5dRUjaUGWPLGtY2p0rA3cIrdOk7wMZ067fcRT69Y1khZEPI14YroBqNzbae6DwBPlVn5g\nqZ3gOJyKFrrGWG8rGikUJOKJHc9b8bTJMANXrFNdgx+A8H65LEcwd179A7dI1ltM4quO7zqrRpYF\nCSsU/Rw2NTET0pRl3kNIIyGDfmUwl4/anzmSXmpd/ljY2viur2pkWZBpv4ABQtq1uHVf96dfu7aH\nc6Ev+XXNTHkFkvXmc8Zy69yOcl8L/D0tSESUGySkXY1bt/vj7vSzkArb28O55Mbh9NLeWNabW9Vh\nuXXYUS7n761FlMMyKWTXuxHzti5lp8QjBEb55MzQCP0o+l9A/bGU/UouWX44+APZYQSqkB0Qr8IQ\n5QYRtE/lVNVCWo9b994Fcd+7tVCzPh5xhYKol+72VOZupaHo1/sWhCTKSXA84iqXLE17x/x5WZ0b\nLIg4KNNQ8kUjioFXjfKI1bh17/sJkawwRDmfRBrBp5dHrMate8+ChHlb9d0ad/IIRfPp5RED29aV\n/SoXhCfKgfQFGafRPOa5dTxCC8I9QXYEMbJtXdmvQkHALIM08zl3QBVy7qkecUQ5HYpHDHDYyCjU\nr/W4dYWCqIH7670ycC/Ie2XesHsviJGY92rWgnz4zd9N5gM/LD2w3Us1w7kN0Z6ghoO6+hDimzTy\nv7/RqYAuQ91fN5GB9JKVMss6iHIDkCl0gig3jFhxs7dKwSzamwFJUxapJQVJmWUdRLkByGSfIMoN\nI/wtv293RzmsN4BZ9ubW1JPNBqJ6yuEAxA/AEOVGEdNTUb/tHeUiZlkvIW0AwhPlRhCejbfSZm/L\nN2J0NsKti1IWjST/byq/spZf6u5+EUeUG4CIVZpbxyMmIytt9rZIXHRSob1FWsvDOGWLnrQgS2ZZ\nF1FuAKIu8Nw6HvH/Y0e5AaLcAITn1vEIsPHcW/Cb3VHOrVBwO7cNQEAW6+fW8YjAxiNpbxS1bm0j\ny0tWWG/oJ8oNQPQNO5PFerl1PCIy0kt7y6JxH67qn9mvVYwsCxIerEM8enVJRiOeQDrArSOe74ok\nkrS3ZTQYpCWHjPyy+9tfjtm2dUlBwnpD/6NXYzLaD1+6noqq8Y1w64Re9nLY7L5eeLM3bF3XS60b\n2rbuYScfqJf/ybatSwoS3gqY7oWGxgEWPSEb6toNVcgOSJiF9ugGHcaAB/GIEE2w2jrQz2KBW2cV\nRIdjOGwIAbLljvRDFbID4lQYtxyANTC71uuS6NFGDG5dpSDuC67XJfgD2YGDKmQHRFUot9yYpAHB\n8IgRiNjJ/1PALAjJSEMIkC4T9T9QhaxrT72kWw5FGRhEjBWkwK0zCyKOMYw0BA055a8moApZ0437\nGLccjjVwtU9Iad+6WkEYRhqChowzaBxDFdJQy5oZtxyYNXCtghT3rTMKEiZIWTqMBgQNaajFzVCF\njPuMY9otNw5hYLLLIwYuWWVunVGQMEM0MpM1g4wGmSnkDVCFzDXyllslysmHimXjGdw6qyBhhpgn\npdACahlkQSVtgipk2l8+59xyY3AGFMIjRiAGt84qSDkf99aLZ+BekIunmDNwLwiXr4tr5wUBcy3I\nug9BLT0wYaliODcRriOo4aCurr3QDLIJCZo4aCLWNpIXpO3CXeOCGTAKwlDYUtZXm1vHI5ABxi/F\npJYwTlWSrLd1jZQLwlDYeKIcj0D+GL8Uk1rCOFXpjPRT61Y2UiwIRWHDykk/t45HTAmk/FJMaqla\niIWRfmrdykaKBWH3eos3VOvk1kU8sU6E5Iv2S3McWZpSXhcjrLcVjZQKQhPS3C0milvHIySLtF+a\n+dhSvRK+d4T1tqaRUkFIQhpPlOMRLlmkX4pZWvIpr//1RpiHta5rpFCQAUIaT5TjEWC99e5yh7xH\nltBUkTe4o9wAIW2AKBeR0Tq5ciN+aeYjS5VCTF0wQlLrVjSSf0Jm4lcXhS0sUXRz63iEyxbpl2Iy\nS1PabTEb6WW9rW3EF8TdKIg2sNb9zjb7ndzljxqNMLBE8SjzKX1YaRvGI2Ca8UsxqSWMU5XOSP9z\nVNcz8ubv1+g3n3vpB8K/hJD2cGB2bXuQ95Ps9uFZX7JLWp1iFhY1uhEjfikmszQNVBMIPt/szUCt\nZuTJ12Hzf6CIvwdIzzaeAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left[\\begin{matrix}\\frac{31}{120} & \\frac{31}{360} & \\frac{1}{40} & \\frac{31}{360} & \\frac{3}{40} & \\frac{11}{360} & \\frac{11}{360} & \\frac{3}{40}\\\\\\frac{31}{360} & \\frac{31}{120} & \\frac{31}{360} & \\frac{1}{40} & \\frac{3}{40} & \\frac{3}{40} & \\frac{11}{360} & \\frac{11}{360}\\\\\\frac{1}{40} & \\frac{31}{360} & \\frac{31}{120} & \\frac{31}{360} & \\frac{11}{360} & \\frac{3}{40} & \\frac{3}{40} & \\frac{11}{360}\\\\\\frac{31}{360} & \\frac{1}{40} & \\frac{31}{360} & \\frac{31}{120} & \\frac{11}{360} & \\frac{11}{360} & \\frac{3}{40} & \\frac{3}{40}\\\\\\frac{3}{40} & \\frac{3}{40} & \\frac{11}{360} & \\frac{11}{360} & \\frac{2}{45} & \\frac{1}{36} & \\frac{1}{45} & \\frac{1}{36}\\\\\\frac{11}{360} & \\frac{3}{40} & \\frac{3}{40} & \\frac{11}{360} & \\frac{1}{36} & \\frac{2}{45} & \\frac{1}{36} & \\frac{1}{45}\\\\\\frac{11}{360} & \\frac{11}{360} & \\frac{3}{40} & \\frac{3}{40} & \\frac{1}{45} & \\frac{1}{36} & \\frac{2}{45} & \\frac{1}{36}\\\\\\frac{3}{40} & \\frac{11}{360} & \\frac{11}{360} & \\frac{3}{40} & \\frac{1}{36} & \\frac{1}{45} & \\frac{1}{36} & \\frac{2}{45}\\end{matrix}\\right]$$"
],
"text/plain": [
"⎡ 31 31 31 11 11 ⎤\n",
"⎢─── ─── 1/40 ─── 3/40 ─── ─── 3/40⎥\n",
"⎢120 360 360 360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 31 31 31 11 11 ⎥\n",
"⎢─── ─── ─── 1/40 3/40 3/40 ─── ─── ⎥\n",
"⎢360 120 360 360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 31 31 31 11 11 ⎥\n",
"⎢1/40 ─── ─── ─── ─── 3/40 3/40 ─── ⎥\n",
"⎢ 360 120 360 360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 31 31 31 11 11 ⎥\n",
"⎢─── 1/40 ─── ─── ─── ─── 3/40 3/40⎥\n",
"⎢360 360 120 360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 11 11 ⎥\n",
"⎢3/40 3/40 ─── ─── 2/45 1/36 1/45 1/36⎥\n",
"⎢ 360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 11 11 ⎥\n",
"⎢─── 3/40 3/40 ─── 1/36 2/45 1/36 1/45⎥\n",
"⎢360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 11 11 ⎥\n",
"⎢─── ─── 3/40 3/40 1/45 1/36 2/45 1/36⎥\n",
"⎢360 360 ⎥\n",
"⎢ ⎥\n",
"⎢ 11 11 ⎥\n",
"⎢3/40 ─── ─── 3/40 1/36 1/45 1/36 2/45⎥\n",
"⎣ 360 360 ⎦"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = zeros(8,8)\n",
"for i in range(8):\n",
" for j in range(8):\n",
" M[i,j] = integrate(M_int[i,j],(r,-1,1), (s,-1,1))\n",
"\n",
"M"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"Ms = symbols('M0:%d'%8)\n",
"\n",
"def f(i,j):\n",
" if i == j:\n",
" return Ms[i]\n",
" else:\n",
" return 0\n",
" \n",
"M_diag = Matrix(8, 8, f)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dist_sq = mat_dist(M_diag, M)**2"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grad = [diff(dist_sq, x) for x in Ms]"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAAAyBAMAAABc9AwEAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJlUzSJmiTKrRN3v\ndrsdCiq5AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKMElEQVR4Ae1cbYhcVxl+Zu7c+Z6dacQGSepO\ngxCjYTNRBAs2O1FBENoMASMYyF6LGylEdwxKCHabMVSkv7IRCQZTdqj0TxSyfoQapclGKrRpzE4D\n/gi07NQP8AuTbVqb7TYZ33Pu951775xzwOku7IG9H+e8z3neeZ97zz3n3jcBqGyosu16WcURONC2\nndPeatiH6/tVGoH0vYrl2VRzlbq47pYbgf2/MI+zy25d6NGXj/yF6lO2uKE2YZX6kSdrVL87rC2u\nThUHZUe3/KcV51BImy4Pwd6xv4f0FFvlQnRLpdxSLAB6FeMG9JvSau2F9l/gwIDu+8lVccqOfgiZ\nd/vdiK1RgCRmcKwd22lfoxeyo8Wby/U+K19F6TbKc4k/z0qrtVjDCj68TVotVZyyoyeBq76fPPhE\nAZKroDg/uGevhRdyzuAtozNeg/7j/BWM1oFFabV+3k68BZSk1VLFKTv6CjDb6v/hcTUKkOIccnfi\n+uxv80Imqrx9kFpkRCOhglrgI6G8Wuo4RUf3VKTVUoCklqTV8kIsmQTU+iwpJn9v0cykqnJvqeOg\n6iguSY8cKpCC9DgDByKqVva7FHQVtX7wS8Ip3FuqOGVHkX+HPJUrChCca8pxkLUDEVULeKKlpBYy\nLymppYxTdrTQlQ6jAgQXpVlciLhaxctqauFqReXeoimaIk7V0TH5MCpAkvLXhAsRVEtvYOS2ilpf\nAm7WFNRSxSk7itKAaXGIlgoQ3AjpZ0CVCxFUq7yEEVrlys8yehXcNBTUUsUpO4pPgJSWKwoQrQq6\nEKWKByKoVrGLwtsqal0DLjUU1FLFKTuaraIkqZYCBMeBP0lpBS9EUK20gfGOiloHod1VmWWo4pQd\n3X/k8E8kw6gA0Z8/sm1OjsYLEVQLz275J7Bt4YIhR4X8lrE2kmdWTg0Jp+zonl5vRdJHBUiq1+tJ\nquWF+NRKz9JMgu6gazVJvx/vMUTpPfaeXqasBUIVHxXiIQTxqYUcWyBmZ9syEee2s03aPf2uNG4t\nEKr4qBAPEYil1gkWbWDfy7RJz/JjmY1+aIbMH78lg+G2a4FQwUeFeAhByl0etak2303+lXZfpfmY\nZEkfqJPKm7uSMGAtECr4qBAPIUhqnkf4d2acq9sB/TmaxEmWZHIe+MpERxIGrAVCBR8V4iEEYV92\ngQytpVipLtJcIe1+gLlgVg/cJks0CHamWrZhgQ2MImUtECr4qBAPMchsg6L6RJuHVjfGW9iZu4XS\ndJNVJF7k1YM3+7Q70Cr0vMt2ufHm+mAMs3AIM2MPsnN5wumjFQaUJgS6DCdA6Pj49PcPCkLgxCMz\nbTCMgHs2JLPzvphflLsCbKJv1qykGxO1RI1eB3wdD5k1gttJ/TY+gGvISi79HMLv4EeCVKaZTZg0\n3FeeIh04hNBEPzM5kHO9GREKZmO7hy/q84IYG5Lr9e7FQD7TwSZ+YQNJlDsl0OPnNyhzgU0Yreh8\npb+zKu6io9PzTpvzNJZ9qF7vlqfNPHQIz+NEw20N4Hpui31kEyabaU+vEoTY6FXLJfT0ZnE5Pv70\nezY728dBYLtXmkl0XJDHPbfSPrIhubNn37DrgD7Iz6i/7S1ukESx+ilMtRJvo1hzIYOPqvi91mDP\nO59ag3E2Ia5goj3Y3LVwCJGsu7WDjxxCvOpVKw7oQM7GWfnbbPckomhDNGhtf2fesySNhEiaE4t9\nSO02cAbZJRQMr9GAY93A5+9Hbl5aLZuQ+h+vDCDxNruE2NzyNgw6dgizO0XVciDPvm4NQYNI6HFs\nxaP88HXBn+VAgI0x/R9rscZfc4tJJP8I9vghteiOEy7pBrZ3QM872XvLJiSmF4TZyNAlfOARGRwc\nwkxCVC0H8kHsEuRy3Juo5ufEMA4EOBmNyJuLq8U2M5lD6Q6yd5GlkZCpFTVjYkln3lKgd4stnCCI\nNRJGzYFG616YS0hriCprkSeEOcuQJpw01QohjPYRo41wHyPjMVoDfxvX714kJHZaPXKLRQkTTdqc\nWeloHe3fKycTd1Bus+qo9dalCmu1S27hBTyMXQsXKrZaUestJ3nHhDqEwGNmjTQhsvxFtCxhomPd\nW/2EkT7+DaNt5mU/JJD65MaDosjXsyHuRYYww6QAQiA06ajzttEq39mb3/rmhGb+O0/H3jD9UWaT\nb9iWgX1wlrGb1lQsWzyx9XqLmT4XsLdP8818yz5m+xRdDiKEm1s6D4eDFSTUtmx503BA7IA+Ez5S\nOd6moygfz2PCd5F6IFHxSBnmvcUIeKEPmKlXsx+j4ygIUh3L1tq5EMD/Dt62u+5fb7H8dzMd+3PY\nxIRK24bBvTbvq2H57zxbPFVLnGcthq/ZPXnqk7u8kRAmLDTSl91ezIR7IULogefWRWChxx8LUT5+\nC4e9VGCpSBYkKh75rlb3YShJqtjrMT2iICgaPgRcSJRama0PeiA8/z3HMrjT88jUqeV++gsr+seX\nvVcGz3+nh+ZVHAK+TYCEEYaiuoWed1ElTpidPsquHruIE+KbyzUbxfYlepX9hWdYX5E+Zg43vQgP\nJDIeeOxwxYd5rYvUUf6pPyqEKLR9CLiQKLX89vwLf5FlcI/MmO8AjKBBxDnL+nyFJZj/w5yiR15Q\nQfwiXRtDJdz3KDDDvRD20YVEDhnBn6V/rYucWWkE2yLOPRAJtXg6dnkGef5Uj+g6WM3U2sMSzJcp\n4bQWbI05J7WGS9h01Irxyt+kAMmkHLX8fUWf+SDhz60gmIJHpbBEkxFryh+0CD+3Mqov3fcmqdUJ\ntwmtHTZhukFqnfoGn0KFOtRfqQDBJFProDnh6u8xtMYHkVHrXJMm+lmZbHFTrfw7OoGsz9OhLvVV\nmmoNj/ABkFpvSL3+UoAkOqTWSEXna7C+3xxa4YfIqHWRLcsU1Cp0FdUaHmGTqQWMzIWGLLRSAZIB\nqUXlX6Edhlb6IRJqJbsYbaqMhGNI0L0lPxIOj1BrmWolra+yoYHzVypA8GlLrR0Nf18xZ36IhFo3\naC1dhSY7y+AJ5vTcGq/FOBVs4iPh8AjpReqjyNVl/umxAiTRZGrRh7xj7eDvjToPQMTV0iiDm1Zq\n6cCqMoqH1/PnFssWp4XkVCXW1N/I1Boi4UOnT7/3TLGOpPiPU4Bop0/Pvtyh7/E7hEMRgIirdZwy\nuEvzSNb9cY09Y2rRm5NS4zXgI7GWgUam1lAJ6ctDkma99YAfsacKkGIXVfuLR2zfbqMXQnNyVizR\n+HHIhoJnpmP/Su57ElOLZ4sXaomXQvqNrBo6Ie6x/2XiUCPSo5AGBUi5i71I/yGks8gqL2Siys1o\nFIgpLP/dTMfe+PqPY+yCTTz/fQ9LME9M32gFW2POh06IJ3unsPeo931bjHtmkwIkc3PZ0Lf+UHgg\nJCIfxFoFjYgP2AN/xbrB/y0C5qdj6DRpWy+rPQKJe5aHi53V7uq6f9h/2QpCellmMF2P3PsRAe1F\nR6MNM++HA+ucEhF4ymDG/wPM6GydeLOv/gAAAABJRU5ErkJggg==\n",
"text/latex": [
"$$\\left \\{ M_{0} : \\frac{31}{120}, \\quad M_{1} : \\frac{31}{120}, \\quad M_{2} : \\frac{31}{120}, \\quad M_{3} : \\frac{31}{120}, \\quad M_{4} : \\frac{2}{45}, \\quad M_{5} : \\frac{2}{45}, \\quad M_{6} : \\frac{2}{45}, \\quad M_{7} : \\frac{2}{45}\\right \\}$$"
],
"text/plain": [
"⎧ 31 31 31 31 ⎫\n",
"⎨M₀: ───, M₁: ───, M₂: ───, M₃: ───, M₄: 2/45, M₅: 2/45, M₆: 2/45, M₇: 2/45⎬\n",
"⎩ 120 120 120 120 ⎭"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"solve(grad, Ms)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Constrained optimization"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"f = mat_dist(M_diag, M)**2 + lamda*(4 - sum(Ms))\n",
"\n",
"var = list(Ms)\n",
"var.append(lamda)"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"grad = [diff(f, x) for x in var]"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"sol = solve(grad, var)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAAyBAMAAACg3s9TAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhCZ3bt2ImbN\nq0Q16XkdAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADg0lEQVRoBe2bv2sUQRTHv16yXnJJLgkqiqAe\nh4RILLYQG4scGrHMNmIhkkMkBBHsRGxuG8FGcmkEC/EsVBCL085fJP9BAmKsQmIrQmJAg5i4vpnJ\nTu72bUiTOcjsTiHs8OHe931vZndvvhEnglUkdWSCoA8HRi4mtX84Ixf6cFC1315Grjjo4tJ8sViI\nc8Q4kBt/UILzeGBJCeEaTADt2oCsjxdo+4vpIAjmeHHAOHANmX/YD3xWQrgGE8CWAVM+nrpYw02g\nh9emGePAF6Af54DXSggXYQLQBuQ++vi25KyjAJzktQHzwGFgsTQIXFFCuAgTgDYg0+1TRdoCQL7K\nawPmgVmPDFi+j6FQSFSGCUAb8F4a0F2moh3RwvK6BQBQ8XqCt3WqJ4XE6Nh1IDTAqQoDzt8RNcU3\nwEYLAFp79EZye90LhTARBoDQgAzkCsjQM9FZ4YVpB5gH5IOm68bscaovhMQMehLtMhAacEb1h34P\nnTMxldECAJgAPqHrJy0BISRm7D6waYBTEAa8AiZd8bznowUA0DkHHAJGXSWEqzAAbBrQtrCweKwa\neJisY4xksNECADiL3PMVeg2pKyFMhAkg3ALAPh9H6DZcw2iZVxYzxgF61+6sHaWnUE0JYTJMAFsG\n9Pr4gLY/wHCZVZYTxoGXxYGv+OHh3aYQJsMEoA3ITG7U8+MTS8Azl1UWE+aB2SBYQ/s8/RhSQpgM\nE4A2gFVLyERqgP45nJBvPNpmugLSFRAeiUXXRkKu0y2QbgG5BXrpJLRx0Bt587AWSLdAugXSp0Ds\n4ROdj+8wxmZsALDtPaDi7dBf9pcNgDZAJoIqe7v+qESt5WvN/TEA360AwhUwTQ/BOZm9dbvOPWqt\nq7k9FRk2Aqg3E+wT9gYQGiATQZm9vQFOU2tXm9tTkWEj4EQMYJ+wJwC9BQoiEZTZ20M6FKP9H2lP\nRYaNQHSJsE/YE4A2QCaCMnvbAJbdyNcvLykytBAIt0CYCFYu/yYDqnEGdMhJ24AtA4ZEf/nVHMVz\nowXZa+QfOwFtgEoEs/62BlgKaANUIjgBh1ZA7BawFNAGyERQZG90DxiOuwlaCmgDZCJI4VztLjDt\nRba/uLQU0AaIRFBmb1PArZj+ZWRoIaANEImgzN6yrhP7+9BSQBsgEkGZvTmnnpTiVoClgDYgruck\nzKUGpGeC6Zlg4v9cPun/YeI/tMrYvNXXsfIAAAAASUVORK5CYII=\n",
"text/latex": [
"$$\\left [ \\frac{437}{720}, \\quad \\frac{437}{720}, \\quad \\frac{437}{720}, \\quad \\frac{437}{720}, \\quad \\frac{283}{720}, \\quad \\frac{283}{720}, \\quad \\frac{283}{720}, \\quad \\frac{283}{720}\\right ]$$"
],
"text/plain": [
"⎡437 437 437 437 283 283 283 283⎤\n",
"⎢───, ───, ───, ───, ───, ───, ───, ───⎥\n",
"⎣720 720 720 720 720 720 720 720⎦"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[sol[key] for key in Ms]"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAAyBAMAAACg3s9TAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhCZdiJmq7vN\nRN2DtcTXAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADuklEQVRoBe2bT2gTQRTGv6SdJmnSJtAKFZTW\nIs1FTBA8J9R6bnoRxEKqhyJVIYpUEMHo1cv25FELXsRDo4IHRQxilQpiTuKtHkUPFkv8W41vtqa7\ns9m320sXstk5hO7k92bn+/bNDJnZYrjxFZ1awo1GCgMTRztVP8TEkRQG/8s/lj4NXB99yLrhAYBE\ngW4/eiHPdsIVwDLFhtP72RZUoHvLALEb18rRGeQ0JtQDANFnZMB0rafK9MEdwI06IN7gLNeCBTAM\n6E2htxJZR/IDE+oBIB4tkQHPEa8wfXAFcGKODEiM4AnTghUwDEhWEa737UWxyoR6AABXC4j9YDqg\nV7sCETLgVtmhCRUwDChWEPlDcewQ8ADQDYinHHovHXIGpL77Ti2ogGFAiDLgJwW+5II9AHQDki/S\n57g+bAOQ+jbG52tsEypgGBBeQ+I7ug/OcJEeALq+4gpCWa4TMgOcAdInNjQMsC2ogGEAXuGQzIB9\neS7UA0DXV0fPLq4P7oDU1yjgZIFrQgVMBsQuja9TUIi9uQeA1EfLUExORvbFFSB9+EvzIDsGVMBk\nANBbj5YRly4wZecB0kdzTUymon1xBaS+PWSAZh8PqIBiQOJ2so44f3PsPED6aJJ3zgBnQOqjVcA5\nAwzAMCA2hMlyqIJ+bhn2ANAnwS6aA1Lc45NjxBmQBuRoDmBbUAHDgMiyWEFMQ67EhHoA6AbgAaay\nTB+2AUh98RHhuAqYAMMAjJ2pAe9mn7L39gCYG/6kITx7h+2DK9Cz+PcuMDaX55qwACYDuAh/1wcG\nbP0c9veDZtUFGRBkwNaWGJslvv4iGALBENCHQLKhljVr2vsWCIZAMASCVaB5MqQMe7GhXNpcFKs2\nleaqtgDAzgFfCmYxNn/3089Ox9IWgMkAOlITp+bz8qNMwvrkh1IsAN4r39JFWwJbGSCP1BJZ8RFT\nWnSI1MSs8qwANAvRnkDTAP1IbQE4gEXgHkk7bpHXAgiLAe0JGENA7iQ9lrtpNPtlaPxb5G1uppqB\nlhSxttAWgGrAL2D18G/ad6tZHr+8lPp8CDSHgNQnvpEBJdpTz7Q8/k0D/AiYDYjSu0KTIxngbck+\nA/wItBqQKHVxGbBpgL8AswGCMmC1hJsXM1n7DPAjYDYANAfkpPalgr0BfgQUAz7TkZrU/tpGv74K\n+BBQDDgPerlqodxtezool0EfAooB/VkxiCuFaY3LAB8CTQP0EzORvpynkznb00G/Ak0D7J55R9QF\nBgR7gsGeYMe/Lt/p/zDxD5AzdYNQEvVBAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left [ \\frac{93}{109}, \\quad \\frac{93}{109}, \\quad \\frac{93}{109}, \\quad \\frac{93}{109}, \\quad \\frac{16}{109}, \\quad \\frac{16}{109}, \\quad \\frac{16}{109}, \\quad \\frac{16}{109}\\right ]$$"
],
"text/plain": [
"⎡ 93 93 93 93 16 16 16 16⎤\n",
"⎢───, ───, ───, ───, ───, ───, ───, ───⎥\n",
"⎣109 109 109 109 109 109 109 109⎦"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M2 = M*4/M.trace()\n",
"[M2[i,i] for i in range(8)]"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAVkAAAAyBAMAAAAEmLXFAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhAimburRN3N\ndmbBWFV7AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC5UlEQVRoBe2avWtTURjGn3zdtvloUlQQHapZ\nFBzM5tigFdc76qLBwQou2ZJsBQXFpXVwEmxcdLS4qFNdxLFOglPzF4hY42DV47lJTmzued5yCz3S\nC/dCIPnx8pzfee+5uR8J5tU3xGJLK1XBkcXLsZBFbvFSBUcDV6/1oE6cl9r3CAXHQgaKPsuQ8Ps9\nirMj22tI/7TLcqs4342MeQa8FWor4Nt9ezzdzGGGsX0MzNllMz5SlciYZ+TurzFbAV/tMFtTbGxP\nAmt1yyu1jpnvFoWAeQawxWwlXGK2ptjYbvrMttintgLmGWYka9J8EpFsddYGbUGBz1XANINrCS2P\naJvn37ovT1lNCQDHPMOFbWGZar2hFBzzDBe2bao1xecgYJ7hwLa0Sm0/UwqOhQwHthfgNWyzTBO3\nbAoB8wzhcBJwpKMs20SJ2F4BPhFbjoUMQUvAkWzvVM88sbW8V9XOemTMM1x8324qtWNrFZVSxFbA\nPAOd+bc9O1rAU492nsrF5lxGKg4hSmzd7ZSkt0lvhx1IVoLzlVDWpwGlrGEG9GtE/B8ykpVg7YsD\nA3Hube4XbcONdwwLxQI+iAyEervhM6/93ODu896ZDyjcf49tverZmhbNNyZtb1bvBuD5JMUQh4vd\nZ5jeLiHzRyvNTmp5TSz0NApeu7YRDhXDfYax3aohuBq/vstJvy1to6yvxnM9ikPFcJ4xXgnPuoOj\nI6SVP4Hg6Ag3cYRDxXCeMbbFcCVM9jD4NFgJUfFgNUUtFqL3yjArASg27WE0+UCpgF1njG0vvmZa\n2XN0DgJ2njG2RXrwRN9SPl23UAA4dp3xzxZzPvNKHWMUAnacYWz186OVmuWlHzZNb1s0eAbFsPsM\nY6t8rPQsr3If08FJI7QJ2H2GsT2uz++NkBSQWkbhh0Ul7D7D2H5B5retNdvDwsfI2H2Gsc232l1b\nCy9aDwkVsPsMY8ucDh9LbN3tk9j1Nl6/+Mfq3xR/Ab6NbifeIncBAAAAAElFTkSuQmCC\n",
"text/latex": [
"$$\\left [ \\frac{2}{3}, \\quad \\frac{2}{3}, \\quad \\frac{2}{3}, \\quad \\frac{2}{3}, \\quad \\frac{1}{3}, \\quad \\frac{1}{3}, \\quad \\frac{1}{3}, \\quad \\frac{1}{3}\\right ]$$"
],
"text/plain": [
"[2/3, 2/3, 2/3, 2/3, 1/3, 1/3, 1/3, 1/3]"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[sum(M[i,:]) for i in range(8)]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Distorted element"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"coords = Matrix([\n",
" [x0,x1,x2,x3,x4,x5,x6,x7], \n",
" [y0,y1,y2,y3,y4,y5,y6,y7]\n",
"])"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"dNdr = Matrix(8,2, lambda i,j: diff(H[i], [r,s][j]))\n",
"jac = coords*dNdr\n",
"detjac = jac.det()"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"M_int = (H*H.T)*detjac"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"area = integrate(detjac,(r,-1,1), (s,-1,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"M_int = simplify((H*H.T)*detjac)\n",
"\n",
"M = zeros(8,8)\n",
"for i in range(8):\n",
" for j in range(8):\n",
" \n",
" M[i,j] = integrate(M_int[i,j],(r,-1,1), (s,-1,1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"f = mat_dist(M_diag, M)**2 + lamda*(area - sum(Ms))\n",
"\n",
"var = list(Ms)\n",
"var.append(lamda)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"grad = [diff(f, x) for x in var]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"sol = solve(grad, var)\n",
"sol"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = open('./styles/custom_barba.css', 'r').read()\n",
" return HTML(styles)\n",
"css_styling()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}