Finding the Minimum Distance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "\n",
    "# Defining symbols\n",
    "xi, eta = sp.symbols(r'\\xi \\eta')\n",
    "psx, psy, psz = sp.symbols('p_{sx} p_{sy} p_{sz}')\n",
    "k, n = sp.symbols('k n')\n",
    "px, py, pz = sp.Function('p_x')(k), sp.Function('p_y')(k), sp.Function('p_z')(k)\n",
    "phi = sp.Function('\\phi')(xi, eta, k)"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.366232Z",
     "start_time": "2024-03-28T22:04:24.999500Z"
    }
   },
   "id": "ac54a455a21b2259",
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "Consider a node that is in proximity to a surface.\n",
    "\n",
    "**insert image**\n",
    "\n",
    "The position on the surface is defined as \n",
    "\n",
    "$$\\vec{p}=\\sum_{k=0}^{n - 1}\\phi_k(\\xi, \\eta)\\vec{p}_k$$\n",
    "\n",
    "The minimum distance can be found by minimizing the function of $\\xi$ and $\\eta$ below.\n",
    "\n",
    "$$d=|\\vec{p}_s - \\vec{p}|$$"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "add22e033124e072"
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "sqrt((p_{sx} - Sum(\\phi(\\xi, \\eta, k)*p_x(k), (k, 0, n - 1)))**2 + (p_{sy} - Sum(\\phi(\\xi, \\eta, k)*p_y(k), (k, 0, n - 1)))**2 + (p_{sz} - Sum(\\phi(\\xi, \\eta, k)*p_z(k), (k, 0, n - 1)))**2)",
      "text/latex": "$\\displaystyle \\sqrt{\\left(p_{sx} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{x}{\\left(k \\right)}\\right)^{2} + \\left(p_{sy} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{y}{\\left(k \\right)}\\right)^{2} + \\left(p_{sz} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{z}{\\left(k \\right)}\\right)^{2}}$"
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = sp.sqrt(\n",
    "    (psx - sp.Sum(phi*px, (k, 0, n - 1)))**2 + (psy - sp.Sum(phi*py, (k, 0, n - 1)))**2 + (psz - sp.Sum(phi*pz, (k, 0, n - 1)))**2)\n",
    "d"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.423035Z",
     "start_time": "2024-03-28T22:04:25.368323Z"
    }
   },
   "id": "5a3e78df592ceaed",
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to find the critical points and set the equations equal to zero. The solution for $\\xi$ and $\\eta$ is the closest point on the surface to the point $\\vec{p}_s$."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "880b1132723b57f0"
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "(-(p_{sx} - Sum(\\phi(\\xi, \\eta, k)*p_x(k), (k, 0, n - 1)))*Sum(p_x(k)*Derivative(\\phi(\\xi, \\eta, k), \\xi), (k, 0, n - 1)) - (p_{sy} - Sum(\\phi(\\xi, \\eta, k)*p_y(k), (k, 0, n - 1)))*Sum(p_y(k)*Derivative(\\phi(\\xi, \\eta, k), \\xi), (k, 0, n - 1)) - (p_{sz} - Sum(\\phi(\\xi, \\eta, k)*p_z(k), (k, 0, n - 1)))*Sum(p_z(k)*Derivative(\\phi(\\xi, \\eta, k), \\xi), (k, 0, n - 1)))/sqrt((p_{sx} - Sum(\\phi(\\xi, \\eta, k)*p_x(k), (k, 0, n - 1)))**2 + (p_{sy} - Sum(\\phi(\\xi, \\eta, k)*p_y(k), (k, 0, n - 1)))**2 + (p_{sz} - Sum(\\phi(\\xi, \\eta, k)*p_z(k), (k, 0, n - 1)))**2)",
      "text/latex": "$\\displaystyle \\frac{- \\left(p_{sx} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{x}{\\left(k \\right)}\\right) \\sum_{k=0}^{n - 1} p_{x}{\\left(k \\right)} \\frac{\\partial}{\\partial \\xi} \\phi{\\left(\\xi,\\eta,k \\right)} - \\left(p_{sy} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{y}{\\left(k \\right)}\\right) \\sum_{k=0}^{n - 1} p_{y}{\\left(k \\right)} \\frac{\\partial}{\\partial \\xi} \\phi{\\left(\\xi,\\eta,k \\right)} - \\left(p_{sz} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{z}{\\left(k \\right)}\\right) \\sum_{k=0}^{n - 1} p_{z}{\\left(k \\right)} \\frac{\\partial}{\\partial \\xi} \\phi{\\left(\\xi,\\eta,k \\right)}}{\\sqrt{\\left(p_{sx} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{x}{\\left(k \\right)}\\right)^{2} + \\left(p_{sy} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{y}{\\left(k \\right)}\\right)^{2} + \\left(p_{sz} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{z}{\\left(k \\right)}\\right)^{2}}}$"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_xi = d.diff(xi)\n",
    "d_xi"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.470279Z",
     "start_time": "2024-03-28T22:04:25.424125Z"
    }
   },
   "id": "b008a2ac45deb2fd",
   "execution_count": 3
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "(-(p_{sx} - Sum(\\phi(\\xi, \\eta, k)*p_x(k), (k, 0, n - 1)))*Sum(p_x(k)*Derivative(\\phi(\\xi, \\eta, k), \\eta), (k, 0, n - 1)) - (p_{sy} - Sum(\\phi(\\xi, \\eta, k)*p_y(k), (k, 0, n - 1)))*Sum(p_y(k)*Derivative(\\phi(\\xi, \\eta, k), \\eta), (k, 0, n - 1)) - (p_{sz} - Sum(\\phi(\\xi, \\eta, k)*p_z(k), (k, 0, n - 1)))*Sum(p_z(k)*Derivative(\\phi(\\xi, \\eta, k), \\eta), (k, 0, n - 1)))/sqrt((p_{sx} - Sum(\\phi(\\xi, \\eta, k)*p_x(k), (k, 0, n - 1)))**2 + (p_{sy} - Sum(\\phi(\\xi, \\eta, k)*p_y(k), (k, 0, n - 1)))**2 + (p_{sz} - Sum(\\phi(\\xi, \\eta, k)*p_z(k), (k, 0, n - 1)))**2)",
      "text/latex": "$\\displaystyle \\frac{- \\left(p_{sx} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{x}{\\left(k \\right)}\\right) \\sum_{k=0}^{n - 1} p_{x}{\\left(k \\right)} \\frac{\\partial}{\\partial \\eta} \\phi{\\left(\\xi,\\eta,k \\right)} - \\left(p_{sy} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{y}{\\left(k \\right)}\\right) \\sum_{k=0}^{n - 1} p_{y}{\\left(k \\right)} \\frac{\\partial}{\\partial \\eta} \\phi{\\left(\\xi,\\eta,k \\right)} - \\left(p_{sz} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{z}{\\left(k \\right)}\\right) \\sum_{k=0}^{n - 1} p_{z}{\\left(k \\right)} \\frac{\\partial}{\\partial \\eta} \\phi{\\left(\\xi,\\eta,k \\right)}}{\\sqrt{\\left(p_{sx} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{x}{\\left(k \\right)}\\right)^{2} + \\left(p_{sy} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{y}{\\left(k \\right)}\\right)^{2} + \\left(p_{sz} - \\sum_{k=0}^{n - 1} \\phi{\\left(\\xi,\\eta,k \\right)} p_{z}{\\left(k \\right)}\\right)^{2}}}$"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_eta = d.diff(eta)\n",
    "d_eta"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.501533Z",
     "start_time": "2024-03-28T22:04:25.471286Z"
    }
   },
   "id": "135d1ddfe063a290",
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "For a standard eight node hex element, the shape functions are defined as\n",
    "\n",
    "$$\\phi_k = \\frac{1}{4}(1-\\xi\\xi_k)(1-\\eta\\eta_k)$$"
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "5b04fa4c38e95d0d"
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "(-(p_{sx} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_x(k)/4, (k, 0, 3)))*Sum(p_x(k)*Derivative((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)/4, \\xi), (k, 0, 3)) - (p_{sy} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_y(k)/4, (k, 0, 3)))*Sum(p_y(k)*Derivative((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)/4, \\xi), (k, 0, 3)) - (p_{sz} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_z(k)/4, (k, 0, 3)))*Sum(p_z(k)*Derivative((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)/4, \\xi), (k, 0, 3)))/sqrt((p_{sx} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_x(k)/4, (k, 0, 3)))**2 + (p_{sy} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_y(k)/4, (k, 0, 3)))**2 + (p_{sz} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_z(k)/4, (k, 0, 3)))**2)",
      "text/latex": "$\\displaystyle \\frac{- \\left(p_{sx} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{x}{\\left(k \\right)}}{4}\\right) \\sum_{k=0}^{3} p_{x}{\\left(k \\right)} \\frac{\\partial}{\\partial \\xi} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right)}{4} - \\left(p_{sy} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{y}{\\left(k \\right)}}{4}\\right) \\sum_{k=0}^{3} p_{y}{\\left(k \\right)} \\frac{\\partial}{\\partial \\xi} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right)}{4} - \\left(p_{sz} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{z}{\\left(k \\right)}}{4}\\right) \\sum_{k=0}^{3} p_{z}{\\left(k \\right)} \\frac{\\partial}{\\partial \\xi} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right)}{4}}{\\sqrt{\\left(p_{sx} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{x}{\\left(k \\right)}}{4}\\right)^{2} + \\left(p_{sy} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{y}{\\left(k \\right)}}{4}\\right)^{2} + \\left(p_{sz} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{z}{\\left(k \\right)}}{4}\\right)^{2}}}$"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Making the substitution\n",
    "xi_k, eta_k = sp.Function(r'\\xi_k')(k), sp.Function(r'\\eta_k')(k)\n",
    "phi_sub = sp.Rational(1, 4)*(1 + xi*xi_k)*(1 + eta*eta_k)\n",
    "d_xi = d_xi.subs([\n",
    "    (phi, phi_sub.simplify()),\n",
    "    (n, 4)\n",
    "])\n",
    "d_xi"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.685690Z",
     "start_time": "2024-03-28T22:04:25.502853Z"
    }
   },
   "id": "5e2d2d19e8dbaf85",
   "execution_count": 5
  },
  {
   "cell_type": "code",
   "outputs": [
    {
     "data": {
      "text/plain": "(-(p_{sx} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_x(k)/4, (k, 0, 3)))*Sum(p_x(k)*Derivative((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)/4, \\eta), (k, 0, 3)) - (p_{sy} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_y(k)/4, (k, 0, 3)))*Sum(p_y(k)*Derivative((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)/4, \\eta), (k, 0, 3)) - (p_{sz} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_z(k)/4, (k, 0, 3)))*Sum(p_z(k)*Derivative((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)/4, \\eta), (k, 0, 3)))/sqrt((p_{sx} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_x(k)/4, (k, 0, 3)))**2 + (p_{sy} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_y(k)/4, (k, 0, 3)))**2 + (p_{sz} - Sum((\\eta*\\eta_k(k) + 1)*(\\xi*\\xi_k(k) + 1)*p_z(k)/4, (k, 0, 3)))**2)",
      "text/latex": "$\\displaystyle \\frac{- \\left(p_{sx} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{x}{\\left(k \\right)}}{4}\\right) \\sum_{k=0}^{3} p_{x}{\\left(k \\right)} \\frac{\\partial}{\\partial \\eta} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right)}{4} - \\left(p_{sy} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{y}{\\left(k \\right)}}{4}\\right) \\sum_{k=0}^{3} p_{y}{\\left(k \\right)} \\frac{\\partial}{\\partial \\eta} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right)}{4} - \\left(p_{sz} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{z}{\\left(k \\right)}}{4}\\right) \\sum_{k=0}^{3} p_{z}{\\left(k \\right)} \\frac{\\partial}{\\partial \\eta} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right)}{4}}{\\sqrt{\\left(p_{sx} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{x}{\\left(k \\right)}}{4}\\right)^{2} + \\left(p_{sy} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{y}{\\left(k \\right)}}{4}\\right)^{2} + \\left(p_{sz} - \\sum_{k=0}^{3} \\frac{\\left(\\eta \\eta_{k}{\\left(k \\right)} + 1\\right) \\left(\\xi \\xi_{k}{\\left(k \\right)} + 1\\right) p_{z}{\\left(k \\right)}}{4}\\right)^{2}}}$"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d_eta = d_eta.subs([\n",
    "    (phi, phi_sub.simplify()),\n",
    "    (n, 4)\n",
    "])\n",
    "d_eta"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.717158Z",
     "start_time": "2024-03-28T22:04:25.686766Z"
    }
   },
   "id": "ac1f855327ce410b",
   "execution_count": 6
  },
  {
   "cell_type": "code",
   "outputs": [],
   "source": [
    "# sol = sp.solve([d_xi.doit(), d_eta.doit()], (xi, eta))\n",
    "# sol"
   ],
   "metadata": {
    "collapsed": false,
    "ExecuteTime": {
     "end_time": "2024-03-28T22:04:25.722301Z",
     "start_time": "2024-03-28T22:04:25.718274Z"
    }
   },
   "id": "c7f859abf43e2cc",
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "The solution is not possible to solve symbolically."
   ],
   "metadata": {
    "collapsed": false
   },
   "id": "5e7425e9138c4bc6"
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}