Problema_Cargas_Presa_Elto3nodos.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ejercicio: Presa de Gravedad con elemento de 3 nodos.\n",
    "\n",
    "Análisis de Estructuras. Grado en Ingeniería Civil.\n",
    "\n",
    "Autores:\n",
    "\n",
    "- Alejandro E. Martínez-Castro, Esther Puertas García, Rafael Gallego.\n",
    "- Departamento de Mecánica de Estructuras e Ingeniería Hidráulica. Universidad de Granada.\n",
    "\n",
    "\n",
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\"><img alt=\"Licencia de Creative Commons\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc-sa/4.0/88x31.png\" /></a><br /><span xmlns:dct=\"http://purl.org/dc/terms/\" property=\"dct:title\">Problema MEF Presa Elementos Triangulares</span> by <span xmlns:cc=\"http://creativecommons.org/ns#\" property=\"cc:attributionName\">Alejandro Martínez Castro, Esther Puertas García, Rafael Gallego Sevilla</span> is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc-sa/4.0/\">Creative Commons Reconocimiento-NoComercial-CompartirIgual 4.0 Internacional License</a>.\n",
    "\n",
    "\n",
    "## 1. Enunciado\n",
    "\n",
    "La figura representa una sección de presa de gravedad, la cual se va a analizar mediante el MEF como un problema de elasticidad bidimensional, de deformación plana. \n",
    "\n",
    "<img src=\"presa_triangulares.PNG\" width=\"50%\" height=\"50%\"/>\n",
    "\n",
    "En la figura, el parámetro h3 = h2/2. Se han tomado 8 nodos y 6 elementos, referenciados mediante las letras a,b,c,d,e,f.  \n",
    "\n",
    "La figura muestra la presión correspondiente a la carga hidráulica. Esta presión está aplicada sobre el paramento vertical, variable en función de $y$. El valor máximo es $Q$.\n",
    "\n",
    "Determinar: \n",
    "- Vector de fuerzas nodales debidas a la presión de agua, en el nodo 5. \n",
    "-  Vector de fuerzas nodales debido a peso propio, en el nodo 4. \n",
    "-  Determinar el término de la matriz de rigidez $\\mathbf{K_{58}}$.\n",
    "\n",
    "\n",
    "**Parámetros**\n",
    "- h1: 10 m \n",
    "- h2: 50 m\n",
    "- h3 = 25 m (puesto que vale h2/2)\n",
    "- b1 = 50 m. \n",
    "- b2 = 5 m. \n",
    "- Peso específico del agua: 10 kN/m$^3$. \n",
    "- Peso específico del hormigón: 25 kN/m$^3$. \n",
    "- Módulo de elasticidad $E = 30$ GPa. \n",
    "- Coeficiente de Poisson: 0.25\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPython console for SymPy 1.0 (Python 2.7.13-32-bit) (ground types: python)\n",
      "\n",
      "These commands were executed:\n",
      ">>> from __future__ import division\n",
      ">>> from sympy import *\n",
      ">>> x, y, z, t = symbols('x y z t')\n",
      ">>> k, m, n = symbols('k m n', integer=True)\n",
      ">>> f, g, h = symbols('f g h', cls=Function)\n",
      ">>> init_printing()\n",
      "\n",
      "Documentation can be found at http://docs.sympy.org/1.0/\n"
     ]
    }
   ],
   "source": [
    "from sympy import init_session\n",
    "init_session()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A continuación se definen los símbolos y las funciones de forma (tomadas del cuaderno del MEF-Elemento Triangular)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x1,y1,x2,y2,x3,y3,A = symbols('x1 y1 x2 y2 x3 y3 A',real=True) # Se definen como reales para evaluación posterior.\n",
    "N1 = (x*(y2 - y3) + x2*y3 - x3*y2 - y*(x2 - x3))/(2*A) \n",
    "N2 = (-x*(y1 - y3) - x1*y3 + x3*y1 + y*(x1 - x3))/(2*A)\n",
    "N3 = (x*(y1 - y2) + x1*y2 - x2*y1 - y*(x1 - x2))/(2*A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Determinación de las fuerzas nodales en el elemento triangular.\n",
    "\n",
    "El vector de fuerzas debido a tensiones en el contorno se calcula mediante una integral de línea, que sólo afecta a los lados cargados (es decir, en los que esté definida la tensión exterior $\\mathbf{t}$). \n",
    "\n",
    "Para un elemento dado, la fuerza nodal equivalente se obtiene evaluando la siguiente integral:\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_t = \\int_{\\Gamma} \\mathbf{N}^T \\cdot \\mathbf{t} \\,d \\Gamma\n",
    "\\end{equation}\n",
    "\n",
    "Como se ha indicado anteriormente, $\\Gamma$ afecta únicamente a los lados en los que se define una tensión exterior aplicada $\\mathbf{t}$. \n",
    "\n",
    "Considerando un sistema de referencia local, en el cual los nodos está numerados como 1-2-3, la matriz de funciones de forma puede escribirse en función de las matrices de los nodos 1, 2 y 3, como sigue: \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{N} = \\left( \\begin{matrix} \\mathbf{N_1} & \\mathbf{N_2} & \\mathbf{N_3} \\end{matrix} \\right) \n",
    "\\end{equation}\n",
    "\n",
    "Y por tanto, \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_t = \\int_{\\Gamma} \\left( \\begin{matrix} \\mathbf{N_1} & \\mathbf{N_2} & \\mathbf{N_3} \\end{matrix} \\right)^T \\cdot \\mathbf{t} \\,d \\Gamma\n",
    "\\end{equation}\n",
    "\n",
    "El vector $\\mathbf{f}_t$ puede a su vez escribirse en función de las componentes del nodo 1, 2 y 3. \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_t = \\begin{pmatrix} \\mathbf{f}_{1} \\\\ \\mathbf{f}_{2} \\\\  \\mathbf{f}_{3}\\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "Finalmente, las componentes vectoriales de $\\mathbf{f}_t$ pueden calcularse como: \n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{pmatrix} \\mathbf{f}_{1} \\\\ \\mathbf{f}_{2} \\\\  \\mathbf{f}_{3}\\end{pmatrix} = \n",
    "\\begin{pmatrix} \\int_{\\Gamma} \\mathbf{N_1}^T \\cdot \\mathbf{t} \\, d\\ \\Gamma \\\\ \n",
    "\\int_{\\Gamma} \\mathbf{N_2}^T \\cdot \\mathbf{t}\\,d\\ \\Gamma \\\\ \\int_{\\Gamma} \\mathbf{N_3}^T  \\cdot \\mathbf{t}\\,d \\Gamma\n",
    "\\end{pmatrix}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Componentes del vector tensión\n",
    "\n",
    "El vector $\\mathbf{t}$ puede escribirse como: \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{t} = \\begin{pmatrix} t_x \\\\ t_y \\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "La figura muestra la carga hidrostática. Claramente, la componente $t_y$ es nula. La componente $t_x$ varía en función de $y$. Considerando el sistema de referencia $x,y$ de la figura, y teniendo en cuenta la variación lineal de la carga, puede definirse la función $t_x(y)$ como sigue:  \n",
    "\n",
    "\\begin{equation}\n",
    "t_x(y) = Q - \\dfrac{Q}{h1 + h2}\\,y\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIIAAAAtBAMAAACQZ671AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMiKJmWaru+9U\ndkQJjQrqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAChElEQVRIDZ2WTWsTURSG38nHNJmmkyC4qC4a\nRqW4E7QuFDXoD+is3Cg0KH6ABGMVwSI0XYgrMbjTVReKX5uo4ErRpYiS2QlCMP8gqUZBqI3nTuZO\nbpjbuTNzNjn3nPM+uXfCOzeANE7fvNSQNqIW9UPAo6jDsjnNKgG7q7JWxFrqFQ3m+xGnZWNXG1TV\nf8ta0WrZX2xualP75uAMnSd+6H+YxhjohRZ+xJeTIj1gsmL/7FQZLxIRihtM1m6Wph38S0RIt5js\nRAk9O+NuJzYlzwgz9O3LmHFhsQnmhtapLNaB5zCasdWuoJPDZfZ7rKDXSEbIzNaO2BlgtrZqJyOQ\n6jHuMO1yKCDUwU8L68iXta0wQriDV4cVpBuGIyF89WoKB596B5idGxIAOCG5gzkhuYM9AnewbJuK\nmkfgDlZMy9oegTsYOw6y2AsMlcFw5yzrqGXtp4w7WPYlihrfQ4vNkYPjh0fwHZyY4Ds4MQHcwckJ\n3MFjwq4341zItJcHhBWl3nNgRc/Bfn/ezyaSVH1iiZ3jJXOwGJviAgXHWxY/TtTFBXOwEKb7kvcL\nPqG3/U/OHCyEPnlen3BemAlPU7fvPREmfML8p++2UA9J5/bhIZApeyOcoP20803oK9dCpF6rXaFb\nMnOLn4UTCgMY67iAB2oCvZffAwWXoFvWnteW1SQVXb90fy5gUX2Uw8jSLTkikJLvIV3FUgXHsdSl\nYmjQm11vSQhzXVxhwnYpVE7NXB/T1bXgHu4Db3OkPqYC0P2CYrcaJNAT3NLpb1VdSTAa0GuSUzwD\nLjrAXSWADwSe5KhhVs0KH1F80r+vUWi2OHnyy+eSuN4+z17/S1sOxofhMFh0K/8BQc26mhjYEIAA\nAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- \\frac{Q y}{h_{1} + h_{2}} + Q$$"
      ],
      "text/plain": [
       "    Q⋅y      \n",
       "- ─────── + Q\n",
       "  h₁ + h₂    "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h1, h2, Q  = symbols('h1 h2 Q',real=True)\n",
    "tx = Q - Q / (h1 + h2) * y\n",
    "tx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determinación del vector de fuerzas en el nodo 5. \n",
    "\n",
    "El vector de fuerzas en el nodo 5 incluye componentes en los elementos $d$ y $f$, pues son  los elementos que tienen el nodo 5 en común y que tienen un lado cargado.\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_5 = \\mathbf{f}_5^{d + f} = \\mathbf{f}_5^d +  \\mathbf{f}_5^f\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1. Determinación del vector de fuerzas en 5 debido al elemento _d_\n",
    "\n",
    "Para este elemento se debe plantear una referencia local, en la que se renumeran los nodos como 1-2-3. En este caso, se ha elegido la numeración 3-6-5 como la 1-2-3, según se muestra la figura. Nótese que esta elección es arbitraria, y cualquiera de los tres nodos (3,6,5) puede elegirse como nodo 1. Lo importante es que, una vez fijado el nodo al que se asigna el 1, el resto de nodos respeten el sentido mostrado en la figura (pues sólo así las funciones de forma pueden expresarse en función del área $A$ en valor absoluto). \n",
    "\n",
    "<img src=\"presa_elementod.PNG\" width=\"50%\" height=\"50%\"/>\n",
    "\n",
    "\n",
    "El nodo 5 del elemento _d_ es el nodo 3 en la referencia local. Por tanto, se determinará la integral siguiente: \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_5^d = \\int_{\\Gamma_d} \\mathbf{N_{3}}^T \\cdot \\mathbf{t}\\,d \\Gamma\n",
    "\\end{equation}\n",
    "\n",
    "En esta integral, $\\Gamma_d$ se refiere a la línea del elemento _d_ en la que está definido $\\mathbf{t}$. \n",
    "\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_5^d = \\int_{\\Gamma_d} \\left( \\begin{matrix} N_3 & 0 \\\\ 0 & N_3 \\end{matrix} \\right) \\cdot \\begin{pmatrix} t_x \\\\ t_y \\end{pmatrix} \\,d \\Gamma\n",
    "\\end{equation}\n",
    "\n",
    "Teniendo en cuenta que $t_y = 0$ y que $t_x$ es una función de $y$ (pues corresponde con la carga hidráulica}, definida anteriormente, es fácil determinar la integral anterior, la cual se extiende al lado 1-3 en coordenadas locales, pues es el lado que está cargado. \n",
    "\n",
    "La carga nodal se obtiene mediante la siguiente integral, en $y$, particularizando la función de forma $N_3(x,y)$ en un punto de la línea vertical del paramento cargado, es decir, en un punto de la forma $(0,y)$. \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_5^d = \\int_{y=h3}^{y = 2\\, h3} \\left( \\begin{matrix} N_3(0,y) & 0 \\\\ 0 & N_3(0,y) \\end{matrix} \\right) \\cdot \\begin{pmatrix} t_x(y) \\\\ 0 \\end{pmatrix} \\,d y = \\begin{pmatrix} \\int_{y=h3}^{y =2 \\, h3} N_3(0,y) \\,t_x(y) \\,dy \\\\ 0 \\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "A continuación se realizan estas operaciones. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACgAAAArBAMAAAD4RNfwAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2rVESJdpki\nZjI6QXVuAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABTklEQVQoFcWSMUvDQBTH/6ltzSWxCRZcI3V1\nUHQSxI4qgp0cVGgFdbW4ODi036A3iGDsJ3CKghYsYgZ346pLFmctKNZB4yu5w9Jk1rfc7/3y3rs7\nLlAmiohidFYALZmK5A0JgFmV/CEBaFiC1c6vPJLIipKAtbNWlGRaj5NCK++W7bGLc8DexC7gXNG4\n4Q5M3sQOUAtQgO7r44Begs1X8OLhEpiDzkde6ZguGtVF3NexBOObZvYqTR+3NKVmKV2wKZL5gGRd\nmSdcBrUN8Qdo65Rl/ZQLMBfZEhXwqF25vqFPT0DKB3MIDXlblasB5RT5wPiKCM8H+7RZL9KeNh0R\nFsJQUM459gT+2xImxN8fxmgXgtiuY2DdmNwCZmKSfs69WH/ZSpDUuiqepH+K+tmfCU5XEmQ7wWlu\ngjyE4Q3qnAstJu9OTrcHC1EOwzcpfwCeg1+py9RxhwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{b_{2} h_{3}}{2}$$"
      ],
      "text/plain": [
       "b₂⋅h₃\n",
       "─────\n",
       "  2  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b2, h3 = symbols('b2 h3',real=True)\n",
    "Area_d = Rational(1,2) * b2 * h3\n",
    "Area_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Se definen las coordenadas (x1,y1,x2,y2,x3,y3) de los vértices del triángulo _d_. Con las coordenadas, y el valor del área, se particulariza la expresión de la función de forma $N_3(x,y)$ al triángulo _d_.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOQAAAAuBAMAAADKCqhsAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAVO8Qq5l2zWa73USJ\nIjJt8O9gAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD2ElEQVRYCe1WO2gUQRj+Lpt77+VOK7WJIKiI\n6BLUYHUHPgtFRdEEBA/xFEFIxMRARFkQtTT4QgsxlYoauMpCRBfRFKKwpWCR62yERFJEiPH8Z3Zm\nZ+buQpqsafJDdr7/MfPt/Ds3X4AltJXOfya39hz635RA3zJlNJ95ubHR9HX5xEbU1yVo7LGtA66x\nm2+Gt1hOyp9/pdxImNu5SUKFZISN1p1R3TVwU+6gkTacpPY6A2FmMEQ6SJZ0z8TJw6bfbbq6d0Nz\nfoVYIRaKlYNEvhCMrZ6NOVvMaVE7pmLZaYkV4hFJ+d6RBc1jY85a21wTRHKjKmOHWCGelZSfVW0T\naspdaSoRgYT2CZJ7K6JzySM3H2szJOXA2SNBdMdx5NZpBQQHV533RfLM9wtl4JOZV17aVbh4yZrD\np/vUveJV7AbGjxFkJiitGafo271HAW8lElNBTjytGb+9BLuyC5aXeprpB24Zec3JFJTT6eCP7bdT\neWcNjxB34yNBUlDGppH3vmIfcoUTSIucmB77yyI96IXt2NMxH1iDeguj8rYq0P2C2TWcgPUn43fM\ngRDWI+6lpqjCHhoa3jA0VALi/Sh6Y5j0Lcrmy/RDHWYTn9XAc21lnMRtWEjyt5mkaEtjlNIGkf2V\nqeZmgVfI/aVowy7pp/e+8BITVWTnMFmT0/iY8TDBIr30x94GmJdSb+xvxOlnQu2zZmGvpVldNXqQ\nicbmXZxyqOsO2+/pICWfxSo+ELa20CN4mzcy1Tjqx2cWXS6QKYMa2uZdR+KyqJaUVWsjRZ4C6RLu\nmSvdBjanHJwrI+b0IUFvNu/xsQ+rqftxnxw63O39yFe9psa2u2mqZjPiI9nfah5D+4BZm8Ze6v0T\nMKk4ycKtLDuqotsrLpCiFWnr9jjFc+LeE7u0jvdQkN2Q1oEKHTLdHgJfyhSYKGw7u6PHJXRNkwFT\nJsYoq9tbvBNuVy0ntiIoeTzrZWsE6Gs2G33JCVeE2bWWLAkHhkywl9as4+KKu8LN+Alxw1i+qvix\nYrVDbaTj0sLGcNoRYbtMR7cgawyZSGvLUUFbvT4l6jrGP5o5Ht9arwPP8UAuZow/LxyUPhOv8Ko3\nZUKTaFm98Pi6UluoiN3R4VXfIBPsdC2+8X9EQhlokInFp5MrchngTiATMh7hyGWAK08nkwmlTNFx\nchlgyiNkIlSm6Ci5DDDlCWRCKFN0fLQylwGmPFImpDJFx5oPZKDTETIRKlN0lIEMkPIImQiVKTrK\nQAZIeaRMRN/YYDPaJS6VKbpd8pWF8hBWyhQxJVcezqGUKWJKrjyco7UyLSr9Px7tD5RmOniAAAAA\nAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{1}{b_{2} h_{3}} \\left(- b_{2} h_{3} + b_{2} y - h_{3} x\\right)$$"
      ],
      "text/plain": [
       "-b₂⋅h₃ + b₂⋅y - h₃⋅x\n",
       "────────────────────\n",
       "       b₂⋅h₃        "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1_d = 0\n",
    "y1_d = h3\n",
    "x2_d = b2\n",
    "y2_d = 2*h3\n",
    "x3_d = 0\n",
    "y3_d = 2*h3\n",
    "N3_d = N3.subs([(x1,x1_d),(x2,x2_d),(x3,x3_d),(y1,y1_d),(y2,y2_d),(y3,y3_d),(A,Area_d)])\n",
    "N3_d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A continuación se determina la integral para definir la componente $x$ de $\\mathbf{f}_5^d$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAAyBAMAAADBxu0RAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMquZiWYiRHa7\nVO+0fdwUAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAIEklEQVRoBd1aXWgcVRT+NrvZySbZZIn0QZQm\nblWoog1qLbW1HRX0sVEREcRuK036EOliMVEiTeyDKD50odBSKCSg6INgQ5WKVHSFUqxUGxBEsCWp\nCoWqNWlTa1vCes79mbkzk8zMZre29ELu3LnnO39z786ce06AqtvQllzVPBEMVv+mCMT/Tk6PpMv1\nVvoGnqm3yFrlpYtN07XK8POvxoaSf+6636fL9TbhEYyO1VtmzfKW2jWLCAiYFL/3zEiAUPNEU4x9\nklrZNzhImvb2DZAFt//GSlvu5b72lhl8cwKJH7qFpDWi38+9qUvMervMwOBWOHxeGt8F2F8IYvwz\nTZVK5SJgPQS8SLS7BT1dFpdauz3I/gs0FFiOJXos56FHFxO9bbuNxrLm85L4LsjeEQT5Z9K04j1I\n5GlLLS0Cc4KekRc/ttr7rydwAWifYL63BXNTgZbJq0vMU5cdV6OOOwhzRvNpMn5Ro3nYLc3qgAMD\nix6PjYbTRGicRvMMXZfaGVqEOrRPxhLngQM5EtVcbLbpkizRKpm6XC2Om/czfqPic+nazXnYE1Mu\nbOHRx8DwCJGty7C66ZostbAl9Wi8abezoMdPnWTrd9OfRxfTZNNuWuIRD5ckn6bCWc352O9xYQuP\nysjQ3gLSc2h45b1jSA29Rg+9Lu2WIv3cf/6yhL8rFRZ4CD5drhbt5ig/aWybkHwuXa2m11RFHnBh\nC46sEViXmNo6i87leG5BYPWEvQ/nkDhfauzRrN/SnvHqSnUpmnZzmJ4McMKWfG6QqNz0scu4dKeS\nEXbppG06y4D2aUzaOB6GrZZmnUV2Fq1lzfeoX1dqi1g8oms319sMrig+N0hUbnpNVXHpj8wS0Q6S\ng/zmwWQPNgN/AHUM3c/l0l1oG4d16GlWcTagKyvctPL5Zb/m87zqlRx12YuQfCpIfDWfX5fP30UU\nr6kqLt1AhKi2lp7wFIPW57AKmVnULXTfR5tvIlnEqI2t6GcV5KZXF6SbRNKruY5xyRlIPjdI1Ks5\nxXRtKr1QynQbw80EBQeNU4RtuYrEFVhTqFvoTgtzYqJzDDuAPuwiFaBN69UVdPMlxn09DslHe4xX\nl5py08dOXz+ivc+A8Jb6h75qM4nD9oYC6GTSVnxLPaJwtjjUP+mxl8i9M01kqlhNegX5dAVW80Du\nts/YJs23RmlSbvrYZVy6M9qaLInE4Sa8Tm/Axi60jxU5QIjmi4F4B9mrfM68QjFI4ndm4De/V1fA\nzZbxEWsX2aD4rALzUVNu+tiRLhOtj/4iWuYvAqSWDK0tpdBK35ah+oXuzb39Y6AFHBwHDtKf3F1e\nXQE38djznx6DrflkkEis2k0vO0RceifLjteO4F0HmC47w9BB1FnCZBabloM9aoauoJv0AbebJM4J\nEolHu+lhl3FpppuFxmtHs2UHGDd0jzhLOPJ4MDpBHYfu1Axd2SkxQ7tae0ZuFk6qSR0k0u2tespk\nl3FpjNDdYd5WsdU4dugecZZwRItf5ugI3x8Tk66uzMsXx8WU2aUrIlzhKRUkmlSKBB1TZVza4SWH\n3j3Bv1HRYofuEWcJLY+vq/FVjq/7uYOrS9z6u9Q6XvmFm59dPruF8fNT4obuUWcJU/pTn8sj/vVK\nkpi2VDmOOEtUKe2GhfvOEjJqJWuNd+INa3sVhq23GazPEipqpZnr42bHg9zoC0rJrRob+ZBYxtIe\nsNlBfqk4ZwkVtbpnCbTXqK0a9mky5Ro171lCRa2k6/qs5jVyEvCdJWTUevO56TtLyKj15nPTd5aQ\nR62bz03fWUJGrTeAm3sHv9dlErLGMxT3qosJY7R7lqBfpoxao15Be/qPeDTLso2pX4wzvcdts46y\nACzARxnSAiZHdJlE0GXFxAeNCRNcxllCR62es4RPNN0mejA8Fm0ElsC6Qhn6giNhXlsdqjFomUN7\nty6TiPk5g+oknaJgJo9xltBRq0meZ9yYQ9tUtBGgROA5yuNNOCK8tjrTwUHzWXR2yTKJJIqKiYtT\nubUomMtAWYeIs4SJleO2bjTORhsBqkButGX9RTD6bA0KNmdo04oyiZwzhjyhU4iUTQuDmfKqHzfM\nsJuGZmNoGrEix26K+otQ4oOFKk6sou3OZRLZjCFPOG6GwxTz4i+UhDU0G0OPESIZK+ovQpMPFqY9\ns30TVJlElC5kxcTJsWs3DdjAUYpbfbAwDbFoB4oxjKA00CVVfxF1lKqM2DEmyySydDHJFRM3x67d\nBDTMKjV20Rb2wmK5EgZaqURGGJHsMesoVRnRdlqVSURCbTNXTFSO3ahiABqWLKUuQ3LULRWPdEGJ\nlBn3BY3opbS4qL+IOooJC3uG9N0soXVOlkmkBlExcXPsajUNWHIsQx8vHyxcSTT1EEGkSPGsfdKd\nLdVSoFKKqL+IOooPFqKGikmtc7JMItxUQyfHrjSYMGRnVGGlbqn4bAH7oo0AfqJlcesofltD3Gwr\nI3lZlUn4Qaqh8+8xyk0ThuR4ABaiIQ7pSeCLaCOQKqDFqb+sqcaIphFMjtM7hcsk7KYaIl2W5ik3\nTRgXPfwwCV5sn1k1uKU72gjsHhx41qijVGPEh73f6TIJu9kqKyayDEFm65+FAeMsuR+2WAclXwOl\nPrqVyDAjVlQqF4w6yiKNYA2yuTl27aam0PUjfBCAGeQahzGNUP9sQ8pcW2NqdkoXxr/HGFUMJSV1\n36lv1DB2Kj6mAQyLZ4RRR6nWCKN0EZZjb6tUppXZYbAqPDOhMY0w6ighRvwHFVyS975fvwUAAAAA\nSUVORK5CYII=\n",
      "text/latex": [
       "$$- \\frac{7 Q h_{3}^{3}}{3 h_{1} h_{3} + 3 h_{2} h_{3}} - Q h_{3} + \\frac{3 h_{3}^{2} \\left(Q h_{1} + Q h_{2} + Q h_{3}\\right)}{2 h_{1} h_{3} + 2 h_{2} h_{3}}$$"
      ],
      "text/plain": [
       "             3                   2                     \n",
       "       7⋅Q⋅h₃                3⋅h₃ ⋅(Q⋅h₁ + Q⋅h₂ + Q⋅h₃)\n",
       "- ───────────────── - Q⋅h₃ + ──────────────────────────\n",
       "  3⋅h₁⋅h₃ + 3⋅h₂⋅h₃              2⋅h₁⋅h₃ + 2⋅h₂⋅h₃     "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f5d_x = integrate(N3_d.subs(x,0) * tx,(y,h3,2*h3))\n",
    "f5d_x.simplify()\n",
    "f5d_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 2.2. Determinación del vector de fuerzas en 5 debido al elemento elemento _f_\n",
    "\n",
    "Para este elemento se plantea nuevamente una referencia local, en la que se renumeran los nodos como 1-2-3. En este caso, se ha elegido la numeración 5-8-7 como la 1-2-3, según se muestra la figura. \n",
    "\n",
    "<img src=\"presa_elementof.PNG\" width=\"50%\" height=\"50%\"/>\n",
    "\n",
    "\n",
    "El nodo 5 del elemento _f_ es el nodo 1 en la referencia local. Por tanto, se determinará la integral siguiente: \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_5^f = \\int_{\\Gamma_f} \\mathbf{N_{1}}^T \\cdot \\mathbf{t}\\,d \\Gamma\n",
    "\\end{equation}\n",
    "\n",
    "Es decir\n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_5^f = \\int_{\\Gamma_f} \\left( \\begin{matrix} N_1 & 0 \\\\ 0 & N_1 \\end{matrix} \\right) \\cdot \\begin{pmatrix} t_x \\\\ t_y \\end{pmatrix} \\,d \\Gamma\n",
    "\\end{equation}\n",
    "\n",
    "Y teniendo nuevamente en cuenta que $t_y = 0$ y que $t_x$ es una función de $y$, la componente no nula a evaluar de $\\mathbf{f}_5$ será\n",
    "\n",
    "\\begin{equation}\n",
    "f_{5x}^d = \\int_{h2}^{h1} N_1(0,y) \\, t_x(y) \\,dy\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAACgAAAArBAMAAAD4RNfwAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2rVESJdpki\nZjI6QXVuAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABQklEQVQoFcWSv0vDQBTHv6mmTWpMgoNroK4O\ngk6CmFFFMJOLQiuoq8XFwaH5D5rBxR+Ds1MUtWARO7gbV12yOGvBH3UxPpt7GHqd9S3vc597797B\nHZQxF2mMTAmgpFaYVxgAq8r8wQDUbcFa+1fuMeouE7B80kgXauNhXGjlzXZa+tkp4KxiEzA92ii0\nYQVH2ABqMUowL12SQx6cYBHPLVwA01T0I9UQ9eoc7nzMw/gS0opwQ5eq2UoH+gRLX5mhhgUMv2Ag\nuE/b81EuBPQQeQ+WH6RSubqmwkcgF0E/EO1kAC3Q4i5k5dPONg3rRoGGpTGbJIKM8/djgf+Wkj7x\n95cxmqVYmjoKvSPJNWBSkvQ5t6T+st1HUusSP0nmGO0zs2AcrDBlcjPDjEX6KFLswmj1SjNEUZK3\n+4frvYUoJ8kry2+7p1niEPilzAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{b_{2} h_{1}}{2}$$"
      ],
      "text/plain": [
       "b₂⋅h₁\n",
       "─────\n",
       "  2  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h1 = symbols('h1',real=True)\n",
    "Area_f = Rational(1,2) * b2 * h1\n",
    "Area_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOwAAAAuBAMAAADZ3eiYAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAVO8Qq5l2zWa73USJ\nIjJt8O9gAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAD80lEQVRYCe1WTWhjVRT+ktc2TfKStBuZ6aYB\noSMqGGREXEiCDGM3QovF2kFsNi24suKUgDLlgVSXFp2CsxjIuJjFzBSyFjFBaAVRiLPquGl3LjtD\nxQrtTPzOve8vLy+vQtJ20wMn79zvfvecew/3vS/AGVvs+bPYwOWZ/bMoi6Hzsiff9/Mmn3yPz2/y\nKfT4jJo8+NbBZNvptttG/2eQsY5nDTcB8W6WWXNn3njRDYOB8U3Vgx4BXag+2jvki3exhG9L17tw\nCCdK3twNhsve0B8lppzR6wzEu9iqD3/ii+3wOfuZG3HnhiV1OzVW7qCZhMS72IaHp0O0ySlbL7i8\nVBMIUN2yHs3IA+Lhlql6uOmLHdQp+6sDAH8yDFDdsj7a56SJh9qQdMy2xNWFERUaC2U81Mdzyl5f\nnHZoCwwS01/edsZ8umWXL3zEXijb5K94qCUtDy5+ajzF5noBZiyPjzVulzUOCsWmOfeugN/Rizdw\nBcguaZJb1jhoDpSgeV9xTjzUUvp8am68gCOzObCEPwaXcEvT7bKxfeQav+NtAe/Qx3dJyM5UZUxz\nTht7huQaNG+M+BhaIcaJeI33/J7YCmZhHKWa2acoxMv4h5O/VCrfVyryWnEfxcYG9qSDL9FngQlW\nqzKEWal8dqlSKTEkjUs1b49j8VCTso4tI/0kVcscAvVm9plG7dPyfayP3MeOsKXsA2RI0GU5dE6b\namBnF5oXWdbf5H8xyFeI/eRhhvLMRrPL5izMF9hbujTZOISZDylbrOFnWSW8H/gUDzX/lTrEqxaQ\nKgM3kSxpulO2ZrxA4K6AvFLDjxFvfNF52q+Bl4dZUniRV8qckkzaJrHOQC79HOqWxuyyA1aSRE0m\ngdcuV2t0luWVOzRt3gdMIB5q6aoHv7Zg8SCyj8tb802N22WN997neFVh7Bx7ZG6FNJmN+K1s81ZI\nXoEnDu3qsaFSeT8/4icZ8Koqs8uqON1I7zKQj6O2jiulYcWTz594oqTBgHroEzhTyH4y+i17yDuj\n7U13Bvhr9GKBQ9UPBcfy6sEr5u5EAMUzy2w23dOQNvVItq1AvNV6jJSV5IIOe6XVUpjzWcxc+7vc\nQSKgeCJ64nXZqVi7evhkXk8LY/FDNw4JKPPHm+xN3BWHgHpsH58iwMhYASBk6P6pcTUkoB4ha/oI\nKQ1R+bR69DF1VCqlIUqSxkU9PNGKWtT7nNIQJUl8JSd8otV76qgMSkOUJAXUI2pR73NKQ0SSgurR\ne+qoDDmtIeOFoHpELep9TmsIJSmoHr2njsqgNYTSFVSPqEX9mvM+/K569Ct1RB5buoRxmmVt6VJl\n8/J7OuZIFxW4i2j1eR//AarUFaVA/NPwAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{1}{b_{2} h_{1}} \\left(- b_{2} y + b_{2} \\left(h_{1} + h_{2}\\right)\\right)$$"
      ],
      "text/plain": [
       "-b₂⋅y + b₂⋅(h₁ + h₂)\n",
       "────────────────────\n",
       "       b₂⋅h₁        "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1_f = 0\n",
    "y1_f = h2\n",
    "x2_f = b2\n",
    "y2_f = h2+h1\n",
    "x3_f = 0\n",
    "y3_f = h2+h1\n",
    "N1_f = N1.subs([(x1,x1_f),(x2,x2_f),(x3,x3_f),(y1,y1_f),(y2,y2_f),(y3,y3_f),(A,Area_f)])\n",
    "N1_f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2kAAAA3BAMAAABkyh+PAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMiKJmWaru+9U\ndkQJjQrqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAMtElEQVR4Ae1bXYhdVxVed+7P3JnMHxEfapCM\nU0OQqgxqKxhjrhX0MRfEtPgyg8GMWkKmEdImCB2Foi/SqihqEEcwVFRsqtIHrelAXqRWZlCxoJQO\nVPTBh0zSVltIM661f9f+PWfP3DSjzoY5Z5/1833r7HXOufeuNQfg1o6TZ6dubQC77OUrMLw0vFzu\ntetxa1dgeL67cWsj+B9n7zz3tZtxhrl7rXPu/ptB+f+E+SCcuxmnu6+XRv00fD2t3NXUWYEz8Egd\ns0Kb8XdlHO6Eo/2MeldVZwUK7rX2UjVgV2Qk94T8EMytV+PsWngr0F448SBA4yezJG/80tOGh8fO\nLCwJ6ScB3vKLUI8SBUa6b9KmfZ22djiUsDplNWJmGDx5/UMWQH2n/yrL4z0YWQYYWqSo/7BWFXvn\nfQDfFUYHcHswbi7BSLcXYF+v/apr5lDCYVcJlsFTFBzaAAqcKkzHnqgweAPVe9+G98aTAJMrgrTq\nCdmYwTtj3zzadinN7k00pnOuwNCgswbN/rh7vi5lh3DYsAxMWDq1AZR6Ju1vO3s1qXvDFXdgFuD7\nAFdoDzC3InbJzRAlYGQDN80+wIR7IiZrCgytGhehdfKzaMqGS/llpqGpZfAUJYc2gBKvvO24e7J5\n45ur7YiH16k+HEce/FSbW8rznSJ95xXcfJwms7ixw2SNwNR4t56YvUMJE/MTPaOiiWVwxGUHLIAy\nx4z1DsranFj3h1bg4G//0oc74a9Tmbjxe8VLpB6mB+Nz+Df0ha9eIIEaJmsCTAoXtNLsXcq7n33G\noWQMxqN8wgIod0547KCsnZqnGB/vNa71Rx6Fj/5ZfOdLhI3izr9IN/oybn6Ef/sPwHcAWtM4paGz\nJsFk1eNLUsW2LuXTm5tM5zA48qIDHkCRY854B2XtSI8C3YSxl2F0mab50aR8weQGbj6If6s9+DG0\nzor7FQ911iSYrHpcRrk78pSMwXXLH7mFMR5A3q9Au4OytkmPp7HXYHga9qxVn8Kk+ERefRQtL+Hf\naYDfoLvIWmdm5vafz8yQSoLJqsfRADRPyRgCz4zALYzxADJOZaodlLVDFHnzKjTnYa4H7hVLKm80\nL5LgyBRuKGt3QRtvPpk1PNT3WlOAyapHmLU8JWNAxNrDLYw1WQC1IaoMd1DWzlOsL6zB/nW4D8C9\nYiOnMUJZG3+dNPiEbNyADgqCrEkwfIBidl8MUPKUjCHwzAjcwhgPIONUptpBWbsy9dY/tf4NVDV+\nsjvlXrGRk5q42ni+d3SRNPhtpLsBe+a/GGZNgoGoeoTfRvKUjCHCnxHRJaIHD0DLtr3fQVkbX1vq\nPDIP1C650QH3io2d5/Nd+Iz4Hgn4lX5kGibX58OsSTAQVY8zAUoFpWUIPLOCw0zLA2Di7U13UNbg\nw9/64wXo0Z1zYg3Pil+xsZNs3XbyA/0WafDRN7oEnZORJ6QCE1WPt4coeUrLQJ73hu6B5GGSOIUx\nHgDq6nQn4qhMOnz+2rfNYe2wBkRuiM1ktNftmwN+xRqhO/khiGWiipYcweeaFIuqR3tWW/F9BaVi\nQI/2Mm8gcAhW1h+hSPzCGIpM2QW7EwCyk5DoUiATb32QvUClSWS0a4cFMHByGc/o4jMmMOeKNVJ3\n8tjYMgm6i1o8dlHNGrR8ZoiqB1aPI6OCUjGg4xAiJur3VtyYZhlibKbscgCFupNwkBnwqduHQE1j\nGjeJUTss9B84uYxpeFP8dBYHkStWBa6/1QM8tNkTsgtK0/78a9HUiKrHXmVEOwtRQakZ5B3k1e9/\nrxCZGJ/SJkOMT5dd6PoynYREl8LtQwgQRMWh6cSB3tAyMX4SazsmFgAZcuOTJNd8kX3r0IqWmmeK\nFti9XfKPPCWl4ta3BvGZTi1pLUQFpWYALIyaboTG18tzZUpLYB+AzpCRsQk9y3UnIdWluIPQbOuD\nvBEVh6YTB3pTOyzZGkmQG+wkuebL72NXrPKwS64h6nzGd3HBzAghUJWhxA+bWbQ4bgDERC8jEw+t\nuybeEXUndCch0aVw+xDSX6JqOo5ZElaGXGctTc5J0/PMFRtd8jRSTBOFyFDi7/lFxPHq93oZmbji\n0xi7E6aTkOhSuH0IGbxE1XT8hErCypDrrKXJOemW5tElL0Mqh8CfFiDr95ZJLSMXd6etOjLD3wGm\nVyG7FMZIh+T2IaRaosayVhJWhlxnLU1u4tzqRJ/fVv3RrxyiuYJeshtxtyZWyyjFsnBqfnpoG3eP\ntTfTSVilLkXYWzrSIxfd+uCosayVhJUh11lLkO99L43I71337NJHrJ6fNspr8hCJCPes6wbCx74h\n4D81M3NoZuYdOJdlfVk4bW0IZeN2Os339MQB22Cd23QSTlOXQvWWWEhuH8KgWjoGB1ASVoocLHaC\nXFNObm5piDWRN8qW/MmJQpAQkSDkoosoGQEe0/I0Rf0eLuOhGOribwqxLJyqrCkDzJFFIWhcONNJ\nSHQpDpGvaX1wVEHnAZaElSHX91qa3JzRViflj7eAqRyCHkWqfu9nTYpl4XRiNuDiAnxI6U5Cqktx\nnuxN64Ojpp6QdcPKkOuspcn5WYh5ZSvN84gs+cJjU55R/jACkXegQqfqRgT3mi7rU+G0+tuI7iR0\nE10Ktw8hy7H5byOK37+YgrDw20iKXGctTR4sT2UrzfMIl7zTH5n2jPKHIUTeXpaEVf3eXx4lFt2g\n4fks0AJqVScBI452Kdw+hOwxSdTYvUa/CeqGlSHXWUuTB6dV2UrzPMIlb/Zb9B929UcIUeVLjz68\nWKkb4WdNicWPqqFeFuhF1KpOQqpL4bU+GGosayVhZch11tLkwWlVt9Jcl3DJm+vtG65NxVEIUeEg\nK1rSyM+achWFU6o/ZIbqTthOQtgRRG/eh2CosayJQpskrAwrQ26yliSXHO0Tn1sxrQ/6RODj2Imf\nOu9buE0NXc/nVmPef+M68AjNuilEpCFCDWnVcFnF8gmNWZ43a0vay8LpSS6iuUtMBVwctpNgs8ZC\nYn0IjirpXEDWFqoMK0MO9lQS5CJuuAfGXjWtj8NSprftRVhdct63OKhVbO9YNfHJxYcDTwrbTeFm\neY3DOtrXjr/TE2cvqpiNaUdGBy7xBaG3nQQwvSXmyPoQEVQXsCSs7ZO/sAIv6R6DeHZT1OrBNX4d\nJmed9y2us3PSU8dqQUn1o4/DC9Xkivbz92mNw9peVo73HPmeD0HHT9MvQbuGxsSFl90J20mI9pZY\nHyKC6gK2lxVTjbC2T/7P9cY13fowDx+15hOXYP80f9/Ca2rIOLmVuvlN3oHDC/MrU+rsgl1S47He\nG3iGgodDkQtfpzsRYLioLmD9f4zAwvVSAF0tcMnFE/I4edlWmr5T8FfKEn/fwmtqGCpj9Xf4mxRa\nBHoAC3ipYFPjXqFJsXr+FYdp4grHlHrggCmiqPxN86r1YftaZs0bd+EHgn3fgk05lLFqvfPZH0iF\nQQADLxWsm8IhcJ7UJFg996rDJHyVY0o/cMAUUUx+7P1TqvUhnt3CRK95+/j9zvsWsqnhF0Cs1Z7N\nzQ3JoRHAwgPV6NMvSag+iynk21i9Fz6somSWhi9BYbYDB2TYNaadS7r1YYzNmsN967Bq37cQ00gB\nRFsZAP19BgUGXtToZTclVoWRGlXIt0A4YwE48qKDNHwRjDUeOKCFrjX79ZR6EUNe56xbgSX2J/j7\nFqepqREpgBir1jQxOghg4C+jKv2ShIqBjLwhWO0PKk9b8zANXxPANxs4oE+QO/4Evqq20hQ9Dnud\nq3ut3YfR6/x9C9HU8AsgzCp8A4rBiwqUZIpVYaTGlqnADK+VYuRFk6bX3SlyjhkPHDBGkpJh/+3x\nFb/HoLKGjcPR66yToabgFkB8K8Gkn7E+vGKK/EOzH4MJmAVgZOWTJHw5lPQYOGBJIL/Ct5j6fo9B\nrfmeZWi+0rWdDDUFtwDiWwl2nTUfXrctvCoM+vgxmJNgARhZ+SQJXw4lPQYOWBLIV2Ds9aDHoNa8\nuwSra6yToab0FgYbvpVQ6az58KqbYaowFkdpws81FoC1Lp4l4YuRlMPAAUsCmXjg3HrQ+tBr/o8H\nfsbft1BNDVMAUTyelZBqBB9edVNMFcaGqjRh1vxWinUpmSXhS0C47cABOXjJ3KyYXnPmzP7pyRRA\nmFpNjVUEgX3RsFWYAMLE4GsMtK8oOk7CF6Ew44EDMux6UxMB61ZoT1sTtwUQrbN7s7QRBLA1eluF\nsa5qZo08lQ3AUxQdJuGLUJjxwAEZdr2pyVpozmritgASmpmshSpgxXBbhfHtmJGrYgG4iqKjJHwR\nCjMeOCDDrjnd/nUzmBuiZri7ZrQC279uBnND7GYjsQL/AXNKYN1n9+dpAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$- \\frac{Q h_{2}^{3}}{3 h_{1}^{2} + 3 h_{1} h_{2}} + \\frac{Q \\left(h_{1} + h_{2}\\right)^{3}}{3 h_{1}^{2} + 3 h_{1} h_{2}} + \\frac{Q h_{2}^{2}}{h_{1}} - \\frac{Q}{h_{1}} \\left(h_{1} + h_{2}\\right)^{2} - \\frac{h_{2}}{h_{1}} \\left(Q h_{1} + Q h_{2}\\right) + \\frac{1}{h_{1}} \\left(h_{1} + h_{2}\\right) \\left(Q h_{1} + Q h_{2}\\right)$$"
      ],
      "text/plain": [
       "           3                     3        2              2                    \n",
       "       Q⋅h₂           Q⋅(h₁ + h₂)     Q⋅h₂    Q⋅(h₁ + h₂)    h₂⋅(Q⋅h₁ + Q⋅h₂) \n",
       "- ─────────────── + ─────────────── + ───── - ──────────── - ──────────────── \n",
       "      2                 2               h₁         h₁               h₁        \n",
       "  3⋅h₁  + 3⋅h₁⋅h₂   3⋅h₁  + 3⋅h₁⋅h₂                                           \n",
       "\n",
       "                         \n",
       "  (h₁ + h₂)⋅(Q⋅h₁ + Q⋅h₂)\n",
       "+ ───────────────────────\n",
       "             h₁          \n",
       "                         "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f5f_x = integrate(N1_f.subs(x,0) * tx,(y,h2,h2+h1))\n",
    "f5f_x.simplify()\n",
    "f5f_x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3. Evaluación numérica para unos parámetros establecidos \n",
    "\n",
    "A continuación se definene los valores numéricos de los parámetros. Para distinguirlos de los símbolos se le añade una \"v\", para indicar que son valores de los parámetros. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKAAAAAPBAMAAACRq9klAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABxUlEQVQ4Ea3Uv2sTYRzH8Xd+XFJj2mbVKWao\nFZdiJ7dbRHDpLRZKh4RiKe2goYhZirn/wDh0UCk4delgBhEnU0RwcDD0H+itghBLKS0NeH6f7+NV\nfHxGb3gd5HP58M33HkKuMR+iBNef1PhzdT6H4FDeegyKL7m30um0WaR8ZllKcs9hctt2Lsf5Pg7B\ne26ieJNRmqYtVuGFZR1WWJw5sYW3KMzhUKgzj+JNrkEJXsKjUBlDF6q2cOJIeh14ZlaieJM2yEYG\nkRQa7hzD1zArLPTkuw7cNsMr3kQm2DdP8CZSXsmEcVY4/XZrFgfGdz/GFm8iE2jfxKncBPm9gyQr\nbH7g0tAhGCfsKPgSaXmqhZWW3ITL+5Pd4UXhCfnt5t8EacSN+4bIl0QEP7SwYTR8/9S92OH0a6bO\nHDiHUax4E6p9U1WtZ8j7jrIJK32mfjqwJ4WJ4k2ozJnCA3I1C7yTfntsSj2Z0AE5qaNY8SY0W9JX\nrFOtKd9qRTl8vwuLsr6eAwuyPos3YaEuhVc7mxuWL9FyYgsDOeNrLA1dSu1gB8WbcGgKB2l6bik3\nZiG/Pn4Iu/IX0HjwD2zOhFic2H54ZSiF//v6BSWpD0zuU7eVAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$2291.66666666667$$"
      ],
      "text/plain": [
       "2291.66666666667"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h1_v = 10.\n",
    "h2_v = 50.\n",
    "h3_v = h2_v / 2.\n",
    "Q_v = 10 * (h2_v + h1_v) # Se ha tomado el peso específico del agua como 10 kN/m3\n",
    "b1_v = 50.\n",
    "b2_v = 5.\n",
    "f5dx= f5d_x.evalf(subs={Q:Q_v,h1:h1_v,h2:h2_v,h3:h3_v})\n",
    "f5dx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKAAAAAOBAMAAABa9wqAAAAALVBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAIom7VJlmdt1E7xDNqxlC\nciAAAAAJcEhZcwAADsQAAA7EAZUrDhsAAAEVSURBVDgRjZK/SgNBGMTnYjD/IMSAvcRGrCzS2Mg9\nwRGQVDZXWfsEEgRL5RobW6tokzKV+CqBEIWIz+Bk0ZvZbpsf38x8O9yyBwxHYxjgIhsd5xBsTWaY\nPJmgKLMa7DNxiOYGQlLSWqM/E1ho4jbHF4SkpDPH+4fAQhPPVe8bQlLChqKMEAleGbt7/8F2ZYbJ\nkt4d14VYLI+oBVuTGSYl2ckVH6IGz5sY3i/46jWSEuC0YosQicacUrA1mWGyBN03HhJi8cpPhGBr\nMsP0n2QD7K0FnjZxAFzkQlLS37JQYKGJnwULhaSke4b2Stjdt3bwAjwOhKRkv0RxI/BHlsA1mp+G\npARP03MYHlx0ppcVBFuTGSYlv4Ls6RU8uwihAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$333.333333333333$$"
      ],
      "text/plain": [
       "333.333333333333"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f5fx= f5f_x.evalf(subs={Q:Q_v,h1:h1_v,h2:h2_v,h3:h3_v})\n",
    "f5fx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resultado final para la componente $x$ de $\\mathbf{f}_5^{d+f}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAADcAAAAOBAMAAAB5gJkQAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABIElEQVQYGTWQMUsDQRCFvyPZRJMgwdLqSGGQ\nIBYBQavrLLOVrSJYqIUhiBYGzD8Igo1KCmsLLe3UgNaHfyDpxMaAhQmmON/ukYWFb+btzLxZgko9\ngnz7BE/Ft+wqOqb6qfQ2+QnmiVpK80nSdGIxNnewD9dkQuopFU9fncYXnMMNHEfdsmJPOS/BOgws\nz1bihkt5mom/0I1d9sFOt/qdlHLf/UhkfiS+CObGZjqk54mMDSaiYAwN56ywaxLLinWkkE3dQJVe\nrMAfjDog0lmTu1nbUgj3Eoc4OoKLjl7I0ECGPgjKB77S0a0qNUBrM7JkQ0rlhmampGK18Z9Qg6Wz\n1iG5pumltMjCJTuPFGJzpdWTRHZay1FKQXXPOtem/R7xDwppWAlpEW9TAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$2625.0$$"
      ],
      "text/plain": [
       "2625.00000000000"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f5x = f5dx + f5fx\n",
    "f5x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "El vector de cargas nodales equivalentes en 5 debido a la acción del agua será: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEYAAAAyBAMAAAATs7BgAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhAimburRN3N\ndmbBWFV7AAAACXBIWXMAAA7EAAAOxAGVKw4bAAABwUlEQVQ4Ee2VPU/CUBSG3wLSFklpos4S4uKg\n4D/ABMMKxsRBE5hkpIPGj4lB49oYVg3/gEYntkZddXDQycEf4OBHHYxaz73tRZRwQ9xMvEk5b+99\ncnrOIW8uJv1HyFbWf8BYYUGGYKkwh3EpwQ7nOVPb2iWducqDq4tSMccOl2enWODMiI1sC6st1Q2U\nKPEMK1XB6CUoJvahNQJ1esiOYJiIuBR5HsWB/px4oVeuYJOipdmIehQ5k/SI0Ux65UowKRuJJ8FQ\nHPVS5zM7JEjh6H6dqYqFxDvFoC+gna5cQmHdtNPYQ71Fqp5GjM1XMB1UPKgTtNOhB5pDP3XrG6Na\nSDkwXgFStFTWwI9v3VAdLowPgJTuIkI0qOboV81RC2u6yfIwpbhQWc9aA0a3dxSB6xjVY3Kl8qR8\nhqpLLK85fpzZdnCAco6ruIW76sgbcIJyXjBJ3/cdRDabCFRteoP3t3jbJKTbO9ODlpjPoHO2/8/I\npvPb+YTG60ndP+fQeDJGGE/GCOPJGGE8GSPMIGOE8aRMaDwZM8y3hPFkeYTxZAzNkBtPxgjjSZnQ\neFKm5zCU/f/7X2CGuZvEBdLfT7BDd9wnyfGNrQ1fC8sAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}2625.0\\\\0\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡2625.0⎤\n",
       "⎢      ⎥\n",
       "⎣  0   ⎦"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f5_agua = Matrix([f5x,0])\n",
    "f5_agua"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Observación:**\n",
    "En el elemento lineal triangular, puede demostrarse (mediante el PTV) que las fuerzas nodales equivalen a las que se obtendrían mediante un reparto coherente con las ecuaciones de la estática. \n",
    "\n",
    "Este resultado es práctico. Observe: \n",
    "\n",
    "- Para el caso de la componente debida al elemento _f_ , la carga triangular equivalente para el nodo 5 puede obtenerse tomando momentos en 7 e igualando el momento producido por la fuerza nodal $f_{5x}^f$ respecto a dicho nodo (actuando con un brazo de 10 m). \n",
    "\n",
    "\\begin{equation}\n",
    "f_{5x}^f \\cdot 10 = \\dfrac{100 \\cdot 10}{2} \\,\\dfrac{2}{3} \\cdot 10 \\Rightarrow f_{5x}^f = 333.3 kN\n",
    "\\end{equation}\n",
    "\n",
    "- Para el caso de la componente debida al elemento _d_ se obtiene el mismo resultado: transfiriendo la carga trapezoidal variable, desde 100 a 350 kN/m$^2$; observe que al nodo 5 le corresponde (descomponiendo la carga en rectangular más triangular):\n",
    "\n",
    "\\begin{equation}\n",
    "f_{5x}^d = \\dfrac{100 \\cdot 25}{2} + \\dfrac{250 \\cdot 25}{2} \\cdot \\dfrac{1}{3} = 2291.667 kN\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## 3. Evaluación de las fuerzas nodales debido a cargas de peso propio\n",
    "\n",
    "En este apartado se va a evaluar el vector de fuerzas por unidad de volumen debido a cargas de peso propio, en el nodo 4. \n",
    "\n",
    "Tal y como se vio en el cuaderno de MEF del Elemento Triangular, la expresión dentro de un elemento de este vector es la siguiente:  \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}^{(e)}_b = \\dfrac{A^{(e)}}{3} \\,\\begin{pmatrix} 0 \\\\ -\\rho \\,g \\\\ 0 \\\\- \\rho \\,g  \\\\ 0 \\\\ -\\rho \\,g \\end{pmatrix}\n",
    "\\end{equation}\n",
    "\n",
    "Donde el superíndice (e) se ha empleado para indicar que los datos son los de un elemento en concreto. \n",
    "\n",
    "Observando la figura de la presa, vemos que el nodo 4 pertenece a los elementos (a,b,c). \n",
    "\n",
    "El vector de fuerzas volumétricas existe para todos los elementos. En este caso, \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{f}_4 = \\mathbf{f}_4^{a+b+c} = \\mathbf{f}_4^{a} + \\mathbf{f}_4^b + \\mathbf{f}_4^c\n",
    "\\end{equation}\n",
    "\n",
    "Observando la expresión de $\\mathbf{f}_b$, la única dependencia del elemento es el área de cada elemento. Por tanto, evaluar el vector de fuerzas por unidad de volumen es sencillo. Basta con evaluar las áreas de cada triángulo, y \"enviar\" 1/3 del valor total al nodo 4. Finalmente, considerar el peso específico, que en este caso es directamente el producto $\\rho\\,g$ y toma el valor 25 kN/m$^3$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Area triángulo a 625.0\n",
      "Area triángulos b y c (áreas iguales) 343.75\n",
      "Fuerza vertical en 4 debida a peso propio -10937.5\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAGAAAAAyBAMAAABR8MP3AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhCZZs3dIrur\nRHbLQ9+lAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACH0lEQVRIDe2WTWsTURSGn0xyO8nkq1BFEaQx\niB9QbOofmFHqulGk4CohYnEhGN24NCBFBBHRTf1YjC7cGnAlKsxCirjRlUvtRhGRWhujIOp4blYy\nmaG52+JZnDvMfZ/7nnPvzDBMhuuMHNPhVyZmj4ys5/jsDFtGl2vloQRAtc57sSslAdmaum0EnIN9\nRsBTcOtxRFJJv+BdzQBQ3wVYMQAcOcy5igkgDkaAcUlI065J0zyBi0bbKgd3Oq7nxGcpV1PxT2XS\nwakDnz0jh1ixvpnk8B/4dwc26S7N790N89WXcL/6iMzdVrWqu86+z5yUIdL0jLwDlzjhq23c850m\nbtsOw1DeJkiFYSUKqOVXdexxrE5hnEKz1KN8fcyHQIRkpz7pIeLwpU4+IN0vd7DWijtpdCywPBFS\n0CkOKAcUe40upW8y7bYlvdXCZKDRpPgnJQ6/Qd3Q0q5OFC588GQYLulhhcy6tUb2B87+UyKxtAvk\n644sEQc0BeAaB3syves1NLR8EI8lDxzUlVsSSx5I04OSsFcP6x5SW+HjQKzTWT/OQZpO68ULfccn\nL1d3tBZewLT4DfeQ72L3RZANyn0NqMGxwXZxqMcBcnBjHXsHc36qS+4nGfku62jCkgzDDtzkqFda\nVJex27grpDXQ6PAAWSUKnJm82ubY6nPYsyD1vll4Bs6EyHIBTmtZKoo6yJ0NIlLSBmqZ3hyA8e+P\n6Q/WX5lhmoFt71VKAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\left[\\begin{matrix}0\\\\-10937.5\\end{matrix}\\right]$$"
      ],
      "text/plain": [
       "⎡   0    ⎤\n",
       "⎢        ⎥\n",
       "⎣-10937.5⎦"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "gamma_h = 25\n",
    "A_a = 1/2. * b1_v * h3_v\n",
    "A_b = 1/2. * (b2_v + b1_v)/2 * h3_v\n",
    "A_c = A_b\n",
    "print \"Area triángulo a\", A_a\n",
    "print \"Area triángulos b y c (áreas iguales)\", A_b\n",
    "f4y = -gamma_h * 1/3.*(A_a + A_b + A_c)\n",
    "print \"Fuerza vertical en 4 debida a peso propio\", f4y\n",
    "f4_peso = Matrix([0,f4y])\n",
    "f4_peso"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Evaluación de un término de la matriz de rigidez completa.\n",
    "\n",
    "A continuación se evalúa un término concreto de la matriz de rigidez. Se evaluará $\\mathbf{K_{58}}$. \n",
    "\n",
    "Los elementos $e$ y $f$ contribuyen a este término. Por tanto: \n",
    "\n",
    "\\begin{equation}\n",
    "\\mathbf{K_{58}} = \\mathbf{K_{58}}^{e + f} = \\mathbf{K_{58}}^{e}+ \\mathbf{K_{58}}^{f}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Problema elástico de deformación plana. \n",
    "\n",
    "Evaluaremos en primer lugar los vales de d11, d12, d33, correspondientes a deformación plana: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "E = 30E6\n",
    "n = 0.25\n",
    "G = E / (2 * (1 + n))\n",
    "d33 = G\n",
    "d11 = E * (1 - n) / ( (1 + n) * (1 - 2 * n))\n",
    "d12 = n / (1 - n) * d11"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### 4.1. Componente de la matriz en el elemento _e_\n",
    "\n",
    "En primer lugar es necesario definir un sentido de numeración. Eligiendo el nodo 5 como nodo 1 local, planteamos:\n",
    "\n",
    "    (5,6,8) en global se referencia como (1,2,3) en local. (Sentido de recorrido antihorario). \n",
    "\n",
    "El término $\\mathbf{K_{58}}^e$ es $\\mathbf{K_{13}}$ en locales. \n",
    "\n",
    "Su expresión puede tomarse del cuaderno del Elemento de 3 nodos\n",
    "\n",
    "\\begin{equation}\n",
    "\\left[\\begin{matrix}\\frac{b_{1} b_{3} d_{11}}{4 A} + \\frac{c_{1} c_{3} d_{33}}{4 A} & \\frac{b_{1} c_{3} d_{12}}{4 A} + \\frac{b_{3} c_{1} d_{33}}{4 A}\\\\\\frac{b_{1} c_{3} d_{33}}{4 A} + \\frac{b_{3} c_{1} d_{12}}{4 A} & \\frac{b_{1} b_{3} d_{33}}{4 A} + \\frac{c_{1} c_{3} d_{11}}{4 A}\\end{matrix}\\right]\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(x1,y1) = 0, h2_v\n",
    "(x2,y2) = b2_v,h2_v\n",
    "(x3,y3) = b2_v,h2_v + h1_v\n",
    "A = 1/2. * b2_v * h1_v\n",
    "\n",
    "(b1, b2, b3, c1, c2, c3) = y2 - y3, y3-y1, y1-y2, x3-x2, x1-x3, x2-x1\n",
    "\n",
    "k58_c = np.array([[b1*b3*d11/(4*A) + c1*c3*d33/(4*A), b1*c3*d12/(4*A) + b3*c1*d33/(4*A)], [b1*c3*d33/(4*A) + b3*c1*d12/(4*A), b1*b3*d33/(4*A) + c1*c3*d11/(4*A)]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[       0. -6000000.]\n",
      " [-6000000.        0.]]\n"
     ]
    }
   ],
   "source": [
    "print k58_c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4.2. Componente de la matriz en el elemento _f_\n",
    "\n",
    "En primer lugar es necesario definir un sentido de numeración. Eligiendo el nodo 5 como nodo 1 local, planteamos\n",
    "\n",
    "    (5,8,7) en global se referencia como (1,2,3) en local. (Sentido de recorrido antihorario).\n",
    "\n",
    "El término $\\mathbf{K_{58}}^f$ será $\\mathbf{K_{12}}$ en planteamiento local.\n",
    "\n",
    "Su expresión es: \n",
    "\n",
    "\\begin{equation}\n",
    "\\left[\\begin{matrix}\\frac{b_{1} b_{2} d_{11}}{4 A} + \\frac{c_{1} c_{2} d_{33}}{4 A} & \\frac{b_{1} c_{2} d_{12}}{4 A} + \\frac{b_{2} c_{1} d_{33}}{4 A}\\\\\\frac{b_{1} c_{2} d_{33}}{4 A} + \\frac{b_{2} c_{1} d_{12}}{4 A} & \\frac{b_{1} b_{2} d_{33}}{4 A} + \\frac{c_{1} c_{2} d_{11}}{4 A}\\end{matrix}\\right]\n",
    "\\end{equation}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "(x1,y1) = 0, h2_v\n",
    "(x2,y2) = b2_v,h2_v + h1_v\n",
    "(x3,y3) = 0,   h2_v + h1_v\n",
    "A = 1/2. * b2_v * h1_v\n",
    "\n",
    "(b1, b2, b3, c1, c2, c3) = y2-y3, y3-y1, y1-y2, x3-x2, x1-x3, x2-x1\n",
    "\n",
    "k58_f = np.array([[b1*b2*d11/(4*A) + c1*c2*d33/(4*A), b1*c2*d12/(4*A) + b2*c1*d33/(4*A)], [b1*c2*d33/(4*A) + b2*c1*d12/(4*A), b1*b2*d33/(4*A) + c1*c2*d11/(4*A)]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[       0. -6000000.]\n",
      " [-6000000.        0.]]\n"
     ]
    }
   ],
   "source": [
    "print k58_f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finalmente, sumando ambas matrices, se obtiene: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[        0. -12000000.]\n",
      " [-12000000.         0.]]\n"
     ]
    }
   ],
   "source": [
    "k58 = k58_c + k58_f\n",
    "print k58"
   ]
  },
  {
   "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}