{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b0ab0f39-0aff-4bf0-ab68-b632e5fd77fb",
   "metadata": {},
   "source": [
    "Please read **carefully** the explanations written between the code cells. Edit or add your own code when you see `#####` symbols."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17e4b0c9-8fc8-4f4f-b41f-0b05c5ee22e5",
   "metadata": {},
   "source": [
    "### Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4a5e7e4-ebc4-4555-a3ff-e095c8ac03bf",
   "metadata": {},
   "source": [
    "`tamaas` should already be installed from last week session. But if that is not the case, run the cell bellow (clic on it and press `shift+enter`), and when the installation is finished ('`Successfully installed ...`'), restart the Python kernel from the menu: `Kernel > Restart`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de9aabcf-4e75-47d9-b0ec-058be0b7f5ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install --user --no-warn-script-location tamaas==2.3.0.post1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62d96db4-132d-40dc-b0e3-2c8eeceefec4",
   "metadata": {},
   "source": [
    "We import the usual libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12ef8dc7-8413-4842-8938-f76a977732d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D # to make 3d plots\n",
    "import tamaas as tm\n",
    "tm.initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fdb80fa-43ad-43e9-8f6a-947977dba487",
   "metadata": {},
   "source": [
    "# Rough contact mechanics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58344f5f-dad6-4103-be18-7376cf8c0496",
   "metadata": {},
   "source": [
    "Tamaas is made for solving the contact between a rigid surface having a given shape (near to flat) and a flat elastic surface, knowing the total applied normal load or **normal pressure** $p_\\text{N}$. You will generate some rough surfaces, solve for the contact, and look at the **real contact area** $A$ between the two contacting surfaces, which is **a fraction** of the total surface area $A_0 = L^2$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cf14371-8c32-4fcb-8197-aabcea62dbc0",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Definition of some functions (run but don't edit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "666dd09f-91a2-46a8-b088-d751f162b568",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_surf(L, n, l_smallest, l_largest, dh_RMS, Hurst, seed=1):\n",
    "    # Create PSD object\n",
    "    spectrum = tm.Isopowerlaw2D()\n",
    "\n",
    "    # Set spectrum parameters\n",
    "    spectrum.q0 = round(L / l_largest) # round: must be integers\n",
    "    spectrum.q1 = round(L / l_largest) # q1 is the \"roll-off\" frequency, usually equal to q0\n",
    "    spectrum.q2 = round(L / l_smallest)\n",
    "    spectrum.hurst = Hurst\n",
    "\n",
    "    # Create surface generator\n",
    "    generator = tm.SurfaceGeneratorFilter2D([n, n])\n",
    "    generator.spectrum = spectrum\n",
    "    generator.random_seed = seed\n",
    "\n",
    "    # Generate surface\n",
    "    h = generator.buildSurface()\n",
    "    h *= dh_RMS / tm.Statistics2D.computeSpectralRMSSlope(h) * L\n",
    "    \n",
    "    x = np.linspace(-L/2, L/2, n)\n",
    "    y = x.copy()\n",
    "    x, y = np.meshgrid(x, y, indexing=\"ij\")\n",
    "    \n",
    "    return x, y, h\n",
    "\n",
    "\n",
    "def plot_surf(x, y, h):\n",
    "    fig = plt.figure(figsize=(5, 5), dpi=120)\n",
    "    ax = fig.add_subplot(111, projection=\"3d\")\n",
    "    surf = ax.plot_surface(x, y, h, rstride=1, cstride=1, cmap=plt.cm.inferno)\n",
    "\n",
    "    z_ratio = 0.2\n",
    "    ax.set_zlim([-L/2 * z_ratio, L/2 * z_ratio])\n",
    "    ax.set_box_aspect([1, 1, z_ratio])\n",
    "    plt.locator_params(axis=\"z\", nbins=2) # reduce number of ticks on z axis\n",
    "\n",
    "    ax.set_xlabel(\"$x$\")\n",
    "    ax.set_ylabel(\"$y$\")\n",
    "    ax.set_zlabel(\"$h$\")\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "import sys, os\n",
    "sys.stderr = open(os.devnull, \"w\") # To get rid of Tamaas console output\n",
    "\n",
    "\n",
    "def solve(h, pressure):\n",
    "    model = tm.ModelFactory.createModel(tm.model_type.basic_2d, [L, L], [n, n])\n",
    "    model.E = E\n",
    "    model.nu = nu\n",
    "    solver = tm.PolonskyKeerRey(model, h, 1e-12)\n",
    "    solver.solve(pressure)\n",
    "    u = model.displacement.copy()\n",
    "    p = model.traction.copy()\n",
    "    return u, p\n",
    "\n",
    "\n",
    "def plot_field(x, y, f, title=\"\"):\n",
    "    plt.figure(figsize=(3.65, 3), dpi=120)\n",
    "    plt.pcolormesh(x, y, f, cmap=plt.cm.inferno, shading=\"gouraud\")\n",
    "    plt.colorbar()\n",
    "    plt.axis(\"equal\")\n",
    "    \n",
    "    plt.xlabel(\"$x$\")\n",
    "    plt.ylabel(\"$y$\")\n",
    "    \n",
    "    if title != \"\":\n",
    "        plt.title(title)\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "    \n",
    "def comp_area_ratio(p):\n",
    "    return np.sum(p > 0) / p.size"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11753f0a-b056-42f9-8398-c8f6f81cadec",
   "metadata": {},
   "source": [
    "## Playground"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a14bd0e-41e4-4f12-b425-7f7af600ec76",
   "metadata": {},
   "source": [
    "We start by defining the geometric and physical parameters of our model. The values below can be kept as is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f9f39b6-9691-4ef5-9d7e-c53efcc54ee9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Surface geometry\n",
    "L = 1                  # side length\n",
    "n = 256                # number of points (discretization)\n",
    "l_smallest = L / (n/2) # smallest wavelength (not smaller than L/n)\n",
    "l_largest  = L / 2     # largest wavelength\n",
    "Hurst = 0.8\n",
    "\n",
    "# Physical parameters\n",
    "E  = 1.0               # Young's modulus\n",
    "nu = 0.3               # Poisson's ratio\n",
    "Es = E / (1 - nu**2)   # Normalized Young's modulus"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c3bd3b0-83b9-4579-ba0d-d05fa4f21fc7",
   "metadata": {},
   "source": [
    "Next, we can generate some rough surfaces. One important parameter is the **RMS** (root mean square) **of the slopes** (derivative of the height distribution) $h'_\\text{RMS} = \\sqrt{\\langle |\\nabla h|^2\\rangle}$. Change this parameter and see how it affects the resulting rough surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d5c95eee-d967-4477-a066-2dc0f9567624",
   "metadata": {},
   "outputs": [],
   "source": [
    "dh_RMS = 0.5 # CHANGE THIS VALUE ###############################################\n",
    "\n",
    "x, y, h = gen_surf(L, n, l_smallest, l_largest, dh_RMS, Hurst)\n",
    "plot_surf(x, y, h)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e82a578e-a199-4d5a-89dd-7deb927ba330",
   "metadata": {},
   "source": [
    "With the defined rough surface, we can solve for the contact for a given normal pressure $p_\\text{N}$. As a result, we can look at the deformation of the elastic surface and at the local pressures. We can also see which parts of the two surfaces are in contact by searching where the local pressure is non-zero."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4d9b0b8-edb7-4694-bfa7-a58397a486c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "pressure = 2e-2 # CHANGE THIS VALUE ############################################\n",
    "\n",
    "u, p = solve(h, pressure)\n",
    "plot_field(x, y, u, \"Displacement\")\n",
    "plot_field(x, y, p, \"Pressure\")\n",
    "plot_field(x, y, p > 0, \"Contact zones\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d74c040-0884-48a6-bb26-5d3b5b40e447",
   "metadata": {},
   "source": [
    "Finally, from the map of the contact zones, we can deduce the fraction of real contact area $A / A_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e01c76fd-8d10-4f9a-8a85-71f118e64910",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Fraction of area in contact: {:.2f}%\".format(comp_area_ratio(p) * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2412022-65d0-4ef4-a3b4-6452a62e64e8",
   "metadata": {},
   "source": [
    "## Fraction of area in contact"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94be5dee-5af6-4c9f-87e8-0accdb4ef79d",
   "metadata": {},
   "source": [
    "The fraction of area in contact depends on several parameters. Let's see which ones affect this quantity and how."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a092272f-3ad0-4e4f-9698-d8b149b2bb8a",
   "metadata": {},
   "source": [
    "### Effect of normal pressure"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab6b932-cd5e-4194-b37e-3979c31570cc",
   "metadata": {},
   "source": [
    "We can run a loop with several values for the normal pressure and record the resulting area ratio. Choose a list of normal pressures to be tested and try to measure a wide range of area ratio, going from 0.1% to 100%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6d1f6388-00e7-4f68-8291-5e2ad4294de0",
   "metadata": {},
   "outputs": [],
   "source": [
    "dh_RMS = 0.5 # CHOOSE A FIXED VALUE FOR h'RMS ##################################\n",
    "_, _, h = gen_surf(L, n, l_smallest, l_largest, dh_RMS, Hurst)\n",
    "\n",
    "pressure_list = np.array([1]) # CREATE A LIST OF VALUES FOR THE PRESSURE #######\n",
    "area_list_1 = []\n",
    "\n",
    "for pressure_i in pressure_list:\n",
    "    _, p = solve(h, pressure_i)\n",
    "    area = comp_area_ratio(p)\n",
    "    area_list_1.append(area)\n",
    "    print(\"pN = {:.2e}, A/A0 = {:.2f}%\".format(pressure_i, area * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "587ff5b4-d8db-4b84-bdae-b5e33dec0242",
   "metadata": {},
   "source": [
    "Now look at a logarithmic plot of $A/A_0$ versus $p_\\text{N}/E^*$. Find the slope of the curve for smaller values of $A/A_0$ (under 0.1, or 10%). Can you state a law for the evolution of $A/A_0$ as a function of $p_\\text{N}/E^*$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "723888ea-826e-4894-8b56-841460417b52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIND THE RIGHT FACTOR AND SLOPE ##############################################\n",
    "factor = 1.0\n",
    "slope  = 0\n",
    "################################################################################\n",
    "\n",
    "plt.figure(figsize=(4, 3), dpi=120)\n",
    "plt.loglog(pressure_list, area_list_1 ,\".-\", label=\"simulations\")\n",
    "plt.loglog(pressure_list, factor * (pressure_list / Es)**slope / dh_RMS, \"k--\", label=\"$A/A_0 \\propto (p_N/E^*)^{\" + str(slope) + \"}$\")\n",
    "\n",
    "plt.xlabel(\"$p_N/E^*$\")\n",
    "plt.ylabel(\"$A/A_0$\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "feefe8f4-742c-4908-9b4f-c888085b211f",
   "metadata": {},
   "source": [
    "### Effect of roughness"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2242fe7c-659e-4a18-b459-de13dfd02761",
   "metadata": {},
   "source": [
    "Now, do the same process but by varying the roughness of the rigid surface. Choose a list of $h'_\\text{RMS}$ to be tested and try to measure a wide range of area ratio, going from 0.1% to 100%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a84d81c8-d22d-4958-a622-006e016df3e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pressure = 2e-3 # CHOOSE A FIXED VALUE FOR THE PRESSURE ###########################\n",
    "\n",
    "dh_RMS_list = np.array([1]) # CREATE A LIST OF VALUES FOR h'RMS ###################\n",
    "area_list_2 = []\n",
    "\n",
    "for dh_RMS_i in dh_RMS_list:\n",
    "    _, _, h = gen_surf(L, n, l_smallest, l_largest, dh_RMS_i, Hurst)\n",
    "    _, p = solve(h, pressure)\n",
    "    area = comp_area_ratio(p)\n",
    "    area_list_2.append(area)\n",
    "    print(\"h'RMS = {:.2e}, A/A0 = {:.2f}%\".format(dh_RMS_i, area * 100))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "154633c9-8b18-45a7-a710-8a3e485c9e8e",
   "metadata": {},
   "source": [
    "Again, look at a logarithmic plot of $A/A_0$ versus $h'_\\text{RMS}$ and find a relation between the two."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccbc62d1-cf55-4a29-8ac2-005d04dc71d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# FIND THE RIGHT FACTOR AND SLOPE ##############################################\n",
    "factor = 1.0\n",
    "slope  = 0\n",
    "################################################################################\n",
    "\n",
    "plt.figure(figsize=(4, 3), dpi=120)\n",
    "plt.loglog(dh_RMS_list, area_list_2 ,\".-\", label=\"simulations\")\n",
    "plt.loglog(dh_RMS_list, factor * pressure * dh_RMS_list**slope / Es, \"k--\", label=\"$A/A_0 \\propto h'_{RMS}^{\" + str(slope) + \"}$\")\n",
    "\n",
    "plt.xlabel(\"$h'_{RMS}$\")\n",
    "plt.ylabel(\"$A/A_0$\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1ef077b-2fbc-4bbd-bd9b-5d1aacfae6cd",
   "metadata": {},
   "source": [
    "### Theory"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b01f5263-39f5-495d-b875-4868a0e01243",
   "metadata": {},
   "source": [
    "From the two power laws you found above, come up with a theoretical estimate for the ratio $A/A_0$ as a function of $p_\\text{N}/E^*$ and $h'_\\text{RMS}$ and see if it maches the simultions for small values of $A/A_0$ (under 0.1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31f9e3ee-1a69-497e-a837-aad5247bd429",
   "metadata": {},
   "outputs": [],
   "source": [
    "def area_ratio_theory(pressure, dh_RMS):\n",
    "    return 0 # DEFINE THE THEORETICAL CONTACT AREA RATIO #######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aff04b62-fee9-4ab9-9a49-1e83e5b95803",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Varying pressure\n",
    "plt.figure(figsize=(4, 3), dpi=120)\n",
    "plt.loglog(pressure_list, area_list_1 ,\".-\", label=\"simulations\")\n",
    "plt.loglog(pressure_list, area_ratio_theory(pressure_list, dh_RMS), \"k--\", label=\"theory\")\n",
    "\n",
    "plt.xlabel(\"$p_N/E^*$\")\n",
    "plt.ylabel(\"$A/A_0$\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# Varying h'RMS\n",
    "plt.figure(figsize=(4, 3), dpi=120)\n",
    "plt.loglog(dh_RMS_list, area_list_2 ,\".-\", label=\"simulations\")\n",
    "plt.loglog(dh_RMS_list, area_ratio_theory(pressure, dh_RMS_list), \"k--\", label=\"theory\")\n",
    "\n",
    "plt.xlabel(\"$h'_{RMS}$\")\n",
    "plt.ylabel(\"$A/A_0$\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python",
   "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": 5
}