{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n\n# Mesh Refinement\nPerform a mesh refinement study.\n\nIn this example the convergence of the torsion constant is investigated through an analysis of an\nI Section. The mesh is refined both by modifying the mesh size and by specifying the number of\npoints making up the root radius. The figure below the example code shows that mesh refinement\nadjacent to the root radius is a far more efficient method in obtaining fast convergence when\ncompared to reducing the mesh area size for the entire section.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 1\n\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport sectionproperties.pre.library.steel_sections as steel_sections\nfrom sectionproperties.analysis.section import Section"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Define mesh sizes\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh_size_list = [200, 100, 50, 20, 10, 5]\nnr_list = [4, 8, 12, 16, 20, 24, 32]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialise result lists\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh_results = []\nmesh_elements = []\nnr_results = []\nnr_elements = []"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calculate reference solution\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=32)\ngeometry.create_mesh(mesh_sizes=[5])  # create mesh\nsection = Section(geometry)  # create a Section object\nsection.calculate_geometric_properties()\nsection.calculate_warping_properties()\nj_reference = section.get_j()  # get the torsion constant"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run through mesh_sizes with n_r = 8\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for mesh_size in mesh_size_list:\n    geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=8)\n    geometry.create_mesh(mesh_sizes=[mesh_size])  # create mesh\n    section = Section(geometry)  # create a Section object\n    section.calculate_geometric_properties()\n    section.calculate_warping_properties()\n\n    mesh_elements.append(len(section.elements))\n    mesh_results.append(section.get_j())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Run through n_r with mesh_size = 10\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for n_r in nr_list:\n    geometry = steel_sections.i_section(d=203, b=133, t_f=7.8, t_w=5.8, r=8.9, n_r=n_r)\n    geometry.create_mesh(mesh_sizes=[10])  # create mesh\n    section = Section(geometry)  # create a Section object\n    section.calculate_geometric_properties()\n    section.calculate_warping_properties()\n\n    nr_elements.append(len(section.elements))\n    nr_results.append(section.get_j())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Convert results to a numpy array and compute the error\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh_results = np.array(mesh_results)\nnr_results = np.array(nr_results)\nmesh_error_vals = (mesh_results - j_reference) / mesh_results * 100\nnr_error_vals = (nr_results - j_reference) / nr_results * 100"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plot the results\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "(fig, ax) = plt.subplots()\nax.loglog(mesh_elements, mesh_error_vals, \"kx-\", label=\"Mesh Size Refinement\")\nax.loglog(nr_elements, nr_error_vals, \"rx-\", label=\"Root Radius Refinement\")\nplt.xlabel(\"Number of Elements\")\nplt.ylabel(\"Torsion Constant Error [%]\")\nplt.legend(loc=\"center left\", bbox_to_anchor=(1, 0.5))\nplt.tight_layout()\nplt.show()"
      ]
    }
  ],
  "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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}