{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Visualize Hertzian Contact Stress {#hertzian_contact_example}\n\nThe following example demonstrates how to use PyVista to visualize\nHertzian contact stress between a cylinder and a flat plate.\n\nThis example loads a dataset, constructs a line to represent the point\nof contact between the cylinder and the block, and samples the stress\nalong that line. Finally, it plots the dataset and the stress\ndistribution.\n\n**Background** Hertzian contact stress refers to the stress that occurs\nbetween two curved surfaces that are in contact with each other. It is\nnamed after Heinrich Rudolf Hertz, a German physicist who first\ndescribed the phenomenon in the late 1800s. Hertzian contact stress is\nan important concept in materials science, engineering, and other fields\nwhere the behavior of materials under stress is a critical\nconsideration.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pyvista as pv\nfrom pyvista import examples"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Load the dataset\n\nStart by loading the dataset using `pyvista.examples`{.interpreted-text\nrole=\"mod\"} module. This module provides access to a range of datasets,\nincluding FEA (finite element analysis) datasets that are useful for\nstress analysis.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh = examples.download_fea_hertzian_contact_cylinder()\nmesh"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Plot the Dataset\n\nPlot the dataset by part ID.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "mesh.plot(scalars='PartID', cmap=['green', 'blue'], show_scalar_bar=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Creating a Line to Denote the Point of Contact\n\nConstruct a line to represent the point of contact between the cylinder\nand the plate.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ypos = 0.024\na = [0.1, ypos, 0.2 - 1e-4]\nb = [0.095, ypos, 0.2 - 1e-4]\nline = pv.Line(a, b, resolution=100)\nline.clear_data()\nline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Sampling the Stress along the Line\n\nWe can sample the Z component stress along the contact edge and compare\nit with expected pressure.\n\nThe expected values array is the Hertzian contact pressure and is the\nanalytical solution to the non-adhesive contact problem. Computation of\nthese values is an exercise left up to the reader (the radius of the\ncylinder is 0.05). See [Contact\nMechanics](https://en.wikipedia.org/wiki/Contact_mechanics)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Sample the stress\nsampled = line.sample(mesh, tolerance=1e-3)\nx_coord = 0.1 - sampled.points[:, 0]\nsamp_z_stress = -sampled['Stress'][:, 2]\n\n# Expected Hertzian contact pressure\nh_pressure = np.array(\n    [\n        [0.0000, 1718094092],\n        [0.0002, 1715185734],\n        [0.0004, 1703502649],\n        [0.0006, 1683850714],\n        [0.0008, 1655946243],\n        [0.001, 1619362676],\n        [0.0012, 1573494764],\n        [0.0014, 1517500856],\n        [0.0016, 1450208504],\n        [0.0018, 1369953775],\n        [0.002, 1274289906],\n        [0.0022, 1159408887],\n        [0.0024, 1018830677],\n        [0.0026, 839747409.8],\n        [0.0028, 587969605.2],\n        [0.003, 0],\n        [0.005, 0],\n    ],\n)\n\nplt.plot(x_coord, samp_z_stress, '.', label='Z Component Stress')\nplt.plot(h_pressure[:, 0], h_pressure[:, 1], label='Hertzian contact pressure')\nplt.legend()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Visualizing the Z Stress Distribution\n\nYou can now visualize the Z stress distribution. Use\n`pyvista.Plotter`{.interpreted-text role=\"class\"} to create a plot\nwindow and add the dataset to it.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pl = pv.Plotter()\nz_stress = np.abs(mesh['Stress'][:, 2])\npl.add_mesh(\n    mesh,\n    scalars=z_stress,\n    clim=[0, 1.2e9],\n    cmap='gouldian',\n    scalar_bar_args={'title': 'Z Component Stress (Pa)', 'color': 'w'},\n    lighting=True,\n    show_edges=False,\n    ambient=0.2,\n)\npl.camera_position = 'xz'\npl.set_focus(a)\npl.camera.zoom(2.5)\npl.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.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}