{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Linear Elastostatics of a Console in 1d\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "A short comparison of different modeling techniques of a simple console. The beam has a rectangular prismatic cross-section and a linear elastic material model, governed by the following parameters:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "L = 100.  # length of the console\nw, h = 10., 10.  # width and height of the rectangular cross section\nF = -100.  # value of the vertical load at the free end\nE = 210000.0  # Young's modulus\nnu = 0.3  # Poisson's ratio"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The resistance of a rectangular section can be characterized by the following constants (the cross section lies in the y-z plane):\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "cross section\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A = w * h  # area\nIy = w * h**3 / 12  # second moment of inertia around the y axis\nIz = w * h**3 / 12  # second moment of inertia around the z axis\nIx = Iy + Iz  # torsional inertia"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 1d Solution - Exact Analytical\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Analytic solution for the tip displacement of an 1d Euler-Bernoulli beam:\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Bernoulli solution\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nEI = E * Iy\nsol_exact = F * L**3 / (3 * EI)\ntol = np.abs(sol_exact / 1000)\nsol_exact"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# 1d Solution - Approximate Numerical\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from sigmaepsilon.solid import Structure, LineMesh, PointData\nfrom neumann.linalg import linspace, Vector\nfrom polymesh.space import StandardFrame, PointCloud, frames_of_lines\nfrom sigmaepsilon.solid.fem.cells import B2 as Beam\n\nimport numpy as np\n\n# model stiffness matrix\nG = E / (2 * (1 + nu))\nHooke = np.array([\n    [E*A, 0, 0, 0],\n    [0, G*Ix, 0, 0],\n    [0, 0, E*Iy, 0],\n    [0, 0, 0, E*Iz]\n])\n\n# space\nGlobalFrame = StandardFrame(dim=3)\n\n# mesh\nnElem = 20  # number of finite elements to use\np0 = np.array([0., 0., 0.])\np1 = np.array([L, 0., 0.])\ncoords = linspace(p0, p1, nElem+1)\ncoords = PointCloud(coords, frame=GlobalFrame).show()\ntopo = np.zeros((nElem, 2), dtype=int)\ntopo[:, 0] = np.arange(nElem)\ntopo[:, 1] = np.arange(nElem) + 1\n\n# support at the leftmost, load at the rightmost node\nloads = np.zeros((coords.shape[0], 6))\nfixity = np.zeros((coords.shape[0], 6)).astype(bool)\nglobal_load_vector = Vector([0., 0, F], frame=GlobalFrame).show()\nloads[-1, :3] = global_load_vector\nfixity[0, :] = True\n\n# pointdata\npd = PointData(coords=coords, frame=GlobalFrame,\n               loads=loads, fixity=fixity)\n\n# celldata\nframes = frames_of_lines(coords, topo)\ncd = Beam(topo=topo, frames=frames)\n\n# set up mesh and structure\nmesh = LineMesh(pd, cd, model=Hooke, frame=GlobalFrame)\nstructure = Structure(mesh=mesh)\nstructure.linsolve()\n\n# postproc\n# 1) displace the mesh\nstructure.nodal_dof_solution(store='dofsol')\ndofsol = structure.mesh.pd['dofsol'].to_numpy()[:, :3]\nlocal_dof_solution = dofsol[-1, :3]\nsol_fem_1d_B2 = local_dof_solution[2]\nsol_fem_1d_B2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh.config['pyvista', 'plot', 'scalars'] = dofsol[:, 2]\nmesh.config['pyvista', 'plot', 'line_width'] = 4\nmesh.pvplot(notebook=True, jupyter_backend='static', window_size=(600, 400),\n            config_key=('pyvista', 'plot'), cmap='plasma')"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "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.8.10"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}