{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: Bars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is from the CALFEM manual (exs2.py)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Purpose:**  \n",
    "\n",
    "Analysis the displacement at point of loading on a plane truss\n",
    "\n",
    "**Description:**  \n",
    "\n",
    "Consider a plane truss consisting of three bars with the properties $E = 200 GPa$,\n",
    "$A_1 = 6.0 · 10^−4 m^2$ , $A_2 = 3.0 · 10^−4 m^2$ and $A_3 = 10.0 · 10^−4 m^2$ , and loaded by a\n",
    "single force $P = 80 kN$. The corresponding finite element model consists of three\n",
    "elements and eight degrees of freedom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![bars1](../images/bars1.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import calfem.core as cfc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Topology matrix Edof**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define element topology. Element number are not used in the **Edof** matrix in CALFEM for Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Edof = np.array([\n",
    "    [1,2,5,6],\n",
    "    [5,6,7,8],\n",
    "    [3,4,5,6]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Stiffness matrix K and load vector f**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = np.matrix(np.zeros((8,8)))\n",
    "f = np.matrix(np.zeros((8,1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element properties**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "E = 2.0e11\n",
    "A1 = 6.0e-4\n",
    "A2 = 3.0e-4\n",
    "A3 = 10.0e-4\n",
    "ep1 = [E,A1]\n",
    "ep2 = [E,A2]\n",
    "ep3 = [E,A3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element coordinates**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex1 = np.array([0., 1.6])\n",
    "ex2 = np.array([1.6, 1.6])\n",
    "ex3 = np.array([0., 1.6])\n",
    "\n",
    "ey1 = np.array([0., 0.])\n",
    "ey2 = np.array([0., 1.2])\n",
    "ey3 = np.array([1.2, 0.])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element stiffness matrices**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ke1 = cfc.bar2e(ex1,ey1,ep1)\t \n",
    "Ke2 = cfc.bar2e(ex2,ey2,ep2)\n",
    "Ke3 = cfc.bar2e(ex3,ey3,ep3)\t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assemble Ke into K**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the topology information, the global stiffness matrix can be generated by\n",
    "assembling the element stiffness matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 7.50e+07,  0.00e+00,  0.00e+00,  0.00e+00, -7.50e+07,\n",
       "          0.00e+00,  0.00e+00,  0.00e+00],\n",
       "        [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,\n",
       "          0.00e+00,  0.00e+00,  0.00e+00],\n",
       "        [ 0.00e+00,  0.00e+00,  6.40e+07, -4.80e+07, -6.40e+07,\n",
       "          4.80e+07,  0.00e+00,  0.00e+00],\n",
       "        [ 0.00e+00,  0.00e+00, -4.80e+07,  3.60e+07,  4.80e+07,\n",
       "         -3.60e+07,  0.00e+00,  0.00e+00],\n",
       "        [-7.50e+07,  0.00e+00, -6.40e+07,  4.80e+07,  1.39e+08,\n",
       "         -4.80e+07,  0.00e+00,  0.00e+00],\n",
       "        [ 0.00e+00,  0.00e+00,  4.80e+07, -3.60e+07, -4.80e+07,\n",
       "          8.60e+07,  0.00e+00, -5.00e+07],\n",
       "        [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,\n",
       "          0.00e+00,  0.00e+00,  0.00e+00],\n",
       "        [ 0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,  0.00e+00,\n",
       "         -5.00e+07,  0.00e+00,  5.00e+07]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cfc.assem(Edof[0,:],K,Ke1)\n",
    "cfc.assem(Edof[1,:],K,Ke2)\n",
    "cfc.assem(Edof[2,:],K,Ke3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Stiffness matrix K:\n",
      "[[ 7.50e+07  0.00e+00  0.00e+00  0.00e+00 -7.50e+07  0.00e+00  0.00e+00\n",
      "   0.00e+00]\n",
      " [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00\n",
      "   0.00e+00]\n",
      " [ 0.00e+00  0.00e+00  6.40e+07 -4.80e+07 -6.40e+07  4.80e+07  0.00e+00\n",
      "   0.00e+00]\n",
      " [ 0.00e+00  0.00e+00 -4.80e+07  3.60e+07  4.80e+07 -3.60e+07  0.00e+00\n",
      "   0.00e+00]\n",
      " [-7.50e+07  0.00e+00 -6.40e+07  4.80e+07  1.39e+08 -4.80e+07  0.00e+00\n",
      "   0.00e+00]\n",
      " [ 0.00e+00  0.00e+00  4.80e+07 -3.60e+07 -4.80e+07  8.60e+07  0.00e+00\n",
      "  -5.00e+07]\n",
      " [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00\n",
      "   0.00e+00]\n",
      " [ 0.00e+00  0.00e+00  0.00e+00  0.00e+00  0.00e+00 -5.00e+07  0.00e+00\n",
      "   5.00e+07]]\n"
     ]
    }
   ],
   "source": [
    "print(\"Stiffness matrix K:\")\n",
    "print(K)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Solve the system of equations**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considering the prescribed displacements in bc, the system of equations is solved\n",
    "using the function `solveq`, yielding displacements `a` and support forces `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Displacements a:\n",
      "[[ 0.        ]\n",
      " [ 0.        ]\n",
      " [ 0.        ]\n",
      " [ 0.        ]\n",
      " [-0.00039793]\n",
      " [-0.00115233]\n",
      " [ 0.        ]\n",
      " [ 0.        ]]\n",
      "Reaction forces r:\n",
      "[[ 29844.55958549]\n",
      " [     0.        ]\n",
      " [-29844.55958549]\n",
      " [ 22383.41968912]\n",
      " [     0.        ]\n",
      " [     0.        ]\n",
      " [     0.        ]\n",
      " [ 57616.58031088]]\n"
     ]
    }
   ],
   "source": [
    "bc = np.array([1,2,3,4,7,8])\n",
    "f[5] = -80e3\n",
    "a, r = cfc.solveq(K,f,bc)\n",
    "\n",
    "print(\"Displacements a:\")\n",
    "print(a)\n",
    "\n",
    "print(\"Reaction forces r:\")\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Element forces**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The vertical displacement at the point of loading is 1.15 mm. The negative values shows the direction is to the bottom, against the direction of node 6. The section forces es1,\n",
    "`es2` and `es3` are calculated using `bar2s` from element displacements `ed1`, `ed2` and `ed3`\n",
    "obtained using extract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N1= -29844.559585492225\n",
      "N2= 57616.580310880825\n",
      "N3= 37305.69948186528\n"
     ]
    }
   ],
   "source": [
    "ed1 = cfc.extract_eldisp(Edof[0,:],a);\n",
    "N1 = cfc.bar2s(ex1,ey1,ep1,ed1)\n",
    "ed2 = cfc.extract_eldisp(Edof[1,:],a);\n",
    "N2 = cfc.bar2s(ex2,ey2,ep2,ed2)\n",
    "ed3 = cfc.extract_eldisp(Edof[2,:],a);\n",
    "N3 = cfc.bar2s(ex3,ey3,ep3,ed3)\n",
    "print(\"N1=\",N1)\n",
    "print(\"N2=\",N2)\n",
    "print(\"N3=\",N3)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.6.10 64-bit",
   "language": "python",
   "name": "python361064bit3b840f9918f246278fc4b65bf6247be2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}