|
| 1 | +{ |
| 2 | + "cells": [ |
| 3 | + { |
| 4 | + "cell_type": "code", |
| 5 | + "execution_count": null, |
| 6 | + "metadata": { |
| 7 | + "collapsed": false |
| 8 | + }, |
| 9 | + "outputs": [], |
| 10 | + "source": [ |
| 11 | + "#\n", |
| 12 | + "# Program 6.8: A falling tablecloth (tablecloth.ipynb)\n", |
| 13 | + "# J Wang, Computational modeling and visualization with Python\n", |
| 14 | + "#\n", |
| 15 | + "\n", |
| 16 | + "import ode, vpmnb as vpm, numpy as np, vpython as vp\n", |
| 17 | + "vec = vp.vector\n", |
| 18 | + "\n", |
| 19 | + "def force(r): # force of particle pair, with relative pos r\n", |
| 20 | + " s = np.sqrt(np.sum(r*r, axis=-1)) # distance \n", |
| 21 | + " s3 = np.dstack((s, s, s)) # make (m,n,3) array \n", |
| 22 | + " return -spring_k*(1.0 - spring_l/s3)*r # Hooke's law \n", |
| 23 | + " \n", |
| 24 | + "def cloth(Y, t): # tablecloth\n", |
| 25 | + " r, v, f = Y[0], Y[1], np.zeros((N,M,3))\n", |
| 26 | + " \n", |
| 27 | + " rtop = r[0:-1, :] - r[1:, :] # rel pos to top neighbor \n", |
| 28 | + " rright = r[:, 0:-1] - r[:, 1:] # rel pos to right neighbor\n", |
| 29 | + " ftop, fright = force(rtop), force(rright) # forces from top, right\n", |
| 30 | + " f[0:-1, :] = ftop # force from top \n", |
| 31 | + " f[:, 0:-1] += fright # force from right \n", |
| 32 | + " f[1:, :] -= ftop # below, left: use 3rd law \n", |
| 33 | + " f[:, 1:] -= fright\n", |
| 34 | + " a = (f - damp*v)/mass + gvec\n", |
| 35 | + " v[0,0], v[0,-1], v[-1,0], v[-1,-1]=0, 0, 0, 0 # fixed coners \n", |
| 36 | + " return np.array([v,a])\n", |
| 37 | + " \n", |
| 38 | + "L, M, N = 2.0, 15, 15 # size, (M,N) particle array\n", |
| 39 | + "h, mass, damp = 0.01, 0.004, 0.01 # keep damp between [.01,.1]\n", |
| 40 | + "x, y = np.linspace(0,L,M), np.linspace(0,L,N) # particle grid\n", |
| 41 | + "r, v = np.zeros((N,M,3)), np.zeros((N,M,3))\n", |
| 42 | + "spring_k, spring_l = 50.0, x[1]-x[0] # spring const., relaxed length\n", |
| 43 | + "r[:,:, 0], r[:,:, 1] = np.meshgrid(x,y) # initialize pos\n", |
| 44 | + "Y, gvec = np.array([r, v]), np.array([0,0,-9.8]) # [v,a], g vector\n", |
| 45 | + "\n", |
| 46 | + "scene = vp.canvas(title='Tablecloth', background=vec(.2,.5,1), \n", |
| 47 | + " up=vec(0,0,1), center=vec(L/2,L/2,-L/4), forward=vec(1,2,-1))\n", |
| 48 | + "vp.points(pos=[(0,0,0),(0,L,0),(L,L,0),(L,0,0)]) # corners\n", |
| 49 | + "x, y, z = r[:,:,0], r[:,:,1], r[:,:,2] # mesh points\n", |
| 50 | + "t=0\n", |
| 51 | + "mesh = vpm.mesh(x, y, z, vp.color.red, vp.color.yellow)\n", |
| 52 | + "\n", |
| 53 | + "while (t<6):\n", |
| 54 | + " vp.rate(100), vpm.wait(scene) # pause if key pressed\n", |
| 55 | + " Y,t = ode.RK4(cloth, Y, 0, h), t+h\n", |
| 56 | + " x, y, z = Y[0,:,:,0], Y[0,:,:,1], Y[0,:,:,2]\n", |
| 57 | + " mesh.move(x, y, z)\n", |
| 58 | + "net = vpm.net(x, y, z, vp.color.yellow, 0.005) # mesh net\n", |
| 59 | + "net.move(x, y, z) # can be slow" |
| 60 | + ] |
| 61 | + }, |
| 62 | + { |
| 63 | + "cell_type": "code", |
| 64 | + "execution_count": null, |
| 65 | + "metadata": { |
| 66 | + "collapsed": true |
| 67 | + }, |
| 68 | + "outputs": [], |
| 69 | + "source": [] |
| 70 | + } |
| 71 | + ], |
| 72 | + "metadata": { |
| 73 | + "kernelspec": { |
| 74 | + "display_name": "VPython", |
| 75 | + "language": "python", |
| 76 | + "name": "vpython" |
| 77 | + }, |
| 78 | + "language_info": { |
| 79 | + "codemirror_mode": { |
| 80 | + "name": "ipython", |
| 81 | + "version": 2 |
| 82 | + }, |
| 83 | + "file_extension": ".py", |
| 84 | + "mimetype": "text/x-python", |
| 85 | + "name": "python", |
| 86 | + "nbconvert_exporter": "python", |
| 87 | + "pygments_lexer": "ipython2", |
| 88 | + "version": "2.7.8" |
| 89 | + } |
| 90 | + }, |
| 91 | + "nbformat": 4, |
| 92 | + "nbformat_minor": 1 |
| 93 | +} |
0 commit comments