{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example: More bars"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example is from Calfem manual (exs4.py).  \n",
    "\n",
    "**Purpose:**  \n",
    "\n",
    "Analysis of a plane truss.\n",
    "\n",
    "**Description:**  \n",
    "\n",
    "Consider a plane truss, loaded by a single force P = 0.5 MN at 30 $^\\circ$ from normal.\n",
    "\n",
    "![more bars](../images/bars3.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The corresponding finite element model consists of ten elements and twelve degrees\n",
    "of freedom.\n",
    "\n",
    "![more bars2](../images/bars2.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, import necessart modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import calfem.core as cfc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The topology is defined by the matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "Edof = np.array([\n",
    "    [1, 2, 5, 6],\n",
    "    [3, 4, 7, 8],\n",
    "    [5, 6, 9, 10],\n",
    "    [7, 8, 11, 12],\n",
    "    [7, 8, 5, 6],\n",
    "    [11, 12, 9, 10],\n",
    "    [3, 4, 5, 6],\n",
    "    [7, 8, 9, 10],\n",
    "    [1, 2, 7, 8],\n",
    "    [5, 6, 11, 12]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A global stiffness matrix K and a load vector f are defined. The load P is divided\n",
    "into x and y components and inserted in the load vector f."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = np.zeros([12,12])\n",
    "f = np.zeros([12,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The element matrices `Ke` are computed by the function `bar2e`. These matrices are\n",
    "then assembled in the global stiffness matrix using the function assem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 25.0e-4\n",
    "E = 2.1e11\n",
    "ep = [E,A]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Element coordinates are defined as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex = np.array([\n",
    "    [0., 2.],\n",
    "    [0., 2.],\n",
    "    [2., 4.],\n",
    "    [2., 4.],\n",
    "    [2., 2.],\n",
    "    [4., 4.],\n",
    "    [0., 2.],\n",
    "    [2., 4.],\n",
    "    [0., 2.],\n",
    "    [2., 4.]\n",
    "])\n",
    "\n",
    "ey = np.array([\n",
    "    [2., 2.],\n",
    "    [0., 0.],\n",
    "    [2., 2.],\n",
    "    [0., 0.],\n",
    "    [0., 2.],\n",
    "    [0., 2.],\n",
    "    [0., 2.],\n",
    "    [0., 2.],\n",
    "    [2., 0.],\n",
    "    [2., 0.]\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All the element matrices are computed and assembled in the loop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "for elx, ely, eltopo in zip(ex, ey, Edof):\n",
    "    Ke = cfc.bar2e(elx, ely, ep)\n",
    "    cfc.assem(eltopo, K, Ke)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The displacements in a and the support forces in r are computed by solving the\n",
    "system of equations considering the boundary conditions in bc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "bc = np.array([1,2,3,4])\n",
    "\n",
    "f[10]=0.5e6*np.sin(np.pi/6)\n",
    "f[11]=-0.5e6*np.cos(np.pi/6)\n",
    "\n",
    "a, r = cfc.solveq(K,f,bc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Displacement: \n",
      "[[ 0.        ]\n",
      " [ 0.        ]\n",
      " [ 0.        ]\n",
      " [ 0.        ]\n",
      " [ 0.00238453]\n",
      " [-0.0044633 ]\n",
      " [-0.00161181]\n",
      " [-0.00419874]\n",
      " [ 0.00303458]\n",
      " [-0.01068377]\n",
      " [-0.00165894]\n",
      " [-0.01133382]]\n"
     ]
    }
   ],
   "source": [
    "print(\"Displacement: \") \n",
    "print(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Support force: \n",
      "[[-8.66025404e+05]\n",
      " [ 2.40086918e+05]\n",
      " [ 6.16025404e+05]\n",
      " [ 1.92925784e+05]\n",
      " [ 4.65661287e-10]\n",
      " [-2.91038305e-11]\n",
      " [ 1.16415322e-10]\n",
      " [ 1.16415322e-10]\n",
      " [-4.65661287e-10]\n",
      " [ 3.49245965e-10]\n",
      " [-2.03726813e-10]\n",
      " [ 4.65661287e-10]]\n"
     ]
    }
   ],
   "source": [
    "print(\"Support force: \")\n",
    "print(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The displacement at the point of loading is $−1.7 · 10^{−3} m$ in the x-direction and\n",
    "$−11.3 · 10^{−3} ~m$ in the y-direction. At the upper support the horizontal force is\n",
    "$−0.866~MN$ and the vertical $0.240 ~MN$. At the lower support the forces are \n",
    "$0.616~MN$ and $0.193~MN$, respectively.\n",
    "Normal forces are evaluated from element displacements. These are obtained from\n",
    "the global displacements a using the function `extract`. The normal forces are evaluated\n",
    "using the function `bar2s`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Element forces:\n",
      "N1 = 625938\n",
      "N2 = -423100\n",
      "N3 = 170640\n",
      "N4 = -12372.8\n",
      "N5 = -69447\n",
      "N6 = 170640\n",
      "N7 = -272838\n",
      "N8 = -241321\n",
      "N9 = 339534\n",
      "N10 = 371051\n"
     ]
    }
   ],
   "source": [
    "ed = cfc.extractEldisp(Edof,a) \n",
    "\n",
    "N = np.zeros([Edof.shape[0]]) \n",
    "\n",
    "print(\"Element forces:\")\n",
    "\n",
    "i = 0\n",
    "for elx, ely, eld in zip(ex, ey, ed):\n",
    "    N[i] = cfc.bar2s(elx, ely, ep, eld)\n",
    "    print(\"N%d = %g\" % (i+1,N[i]))\n",
    "    i+=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The largest normal force N = 0.626 MN is obtained in element 1 and is equivalent\n",
    "to a normal stress σ = 250 MPa."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To reduce the quantity of input data, the element coordinates matrices `Ex` and\n",
    "`Ey` can alternatively be created from a global coordinate matrix `Coord` and a global\n",
    "topology matrix Dof using the function `coordxtr`, i.e."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Coord = np.array([\n",
    "[0, 2],\n",
    "[0, 0],\n",
    "[2, 2],\n",
    "[2, 0],\n",
    "[4, 2],\n",
    "[4, 0]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "Dof = np.array([\n",
    "[1, 2],\n",
    "[3, 4],\n",
    "[5, 6],\n",
    "[7, 8],\n",
    "[9, 10],\n",
    "[11, 12]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "ex, ey = cfc.coordxtr(Edof, Coord, Dof)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 2.],\n",
       "       [0., 2.],\n",
       "       [2., 4.],\n",
       "       [2., 4.],\n",
       "       [2., 2.],\n",
       "       [4., 4.],\n",
       "       [0., 2.],\n",
       "       [2., 4.],\n",
       "       [0., 2.],\n",
       "       [2., 4.]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(ex)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2., 2.],\n",
       "       [0., 0.],\n",
       "       [2., 2.],\n",
       "       [0., 0.],\n",
       "       [0., 2.],\n",
       "       [0., 2.],\n",
       "       [0., 2.],\n",
       "       [0., 2.],\n",
       "       [2., 0.],\n",
       "       [2., 0.]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(ey)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen above, `ex` and `ey` is same as previous."
   ]
  }
 ],
 "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
}