{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# MetaVoxel Tutorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Amira Abdel-Rahman\n",
    "# (c) Massachusetts Institute of Technology 2020\n",
    "\n",
    "# tested using julia 1.2.0 and windows Nvidia geforce gtx 1070 Ti"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import Julia Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "# using Plots\n",
    "import JSON\n",
    "using StaticArrays, BenchmarkTools\n",
    "using Base.Threads\n",
    "using CUDAnative\n",
    "using CuArrays,CUDAdrv \n",
    "import Base: +, * , -, ^"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "axialStrain (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "include(\"./julia/include/vector.jl\") #utils for vectors and quaternions\n",
    "include(\"./julia/include/material.jl\") #utils for node and edge material\n",
    "include(\"./julia/include/export.jl\") #export simulation data to json\n",
    "include(\"./julia/include/run.jl\") #turn setup to cuda arrays and run simulation\n",
    "include(\"./julia/include/updateEdges.jl\") #edges properties update\n",
    "include(\"./julia/include/externalForces.jl\") #external forces applied to the system\n",
    "include(\"./julia/include/forces.jl\") #force integration\n",
    "include(\"./julia/include/updateNodes.jl\") #nodes properties update"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Voxel Design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"tutorial\""
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set name for simulation\n",
    "name= \"tutorial\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.a. Import lines from Rhino (.3dm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhino=true\n",
    "\n",
    "setup = Dict()\n",
    "\n",
    "setup[\"rhino\"]=rhino\n",
    "setup[\"rhinoFileName\"]=\"./julia/examples/trial.3dm\";\n",
    "setup[\"layerIndex\"]=\"1\"; #layer index to import, only lines from these layer will get imported\n",
    "setup[\"voxelList\"]=false\n",
    "\n",
    "# make sure to divide curves into smaller lines, it will only add nodes at the start and end of each line/curve\n",
    "\n",
    "voxelSize=75 #in case you want to array the base rhino curve, what is the size of the voxel\n",
    "latticeSizeX=2 # if you don't want to copy/array make this 1\n",
    "latticeSizeY=2 # if you don't want to copy/array make this 1\n",
    "latticeSizeZ=2 # if you don't want to copy/array make this 1\n",
    "\n",
    "setup[\"latticeSizeX\"]=latticeSizeX\n",
    "setup[\"latticeSizeY\"]=latticeSizeY\n",
    "setup[\"latticeSizeZ\"]=latticeSizeZ\n",
    "\n",
    "gridSize=10 #lattice size\n",
    "setup[\"gridSize\"]=gridSize\n",
    "\n",
    "#scaling params\n",
    "setup[\"voxelSize\"]=voxelSize; #voxel size\n",
    "setup[\"scale\"]=1e4; #scale for visualization\n",
    "setup[\"hierarchical\"]=false; #hierachical simualtion\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.b Draw Lattice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhino=false\n",
    "\n",
    "voxelSize=0.001\n",
    "latticeSizeX=7\n",
    "latticeSizeY=2\n",
    "latticeSizeZ=2\n",
    "\n",
    "gridSize=10\n",
    "\n",
    "setup = Dict()\n",
    "setup[\"gridSize\"]=gridSize\n",
    "\n",
    "setup[\"rhino\"]=false\n",
    "setup[\"voxelList\"]=false\n",
    "\n",
    "\n",
    "setup[\"latticeSizeX\"]=latticeSizeX\n",
    "setup[\"latticeSizeY\"]=latticeSizeY\n",
    "setup[\"latticeSizeZ\"]=latticeSizeZ\n",
    "\n",
    "#scaling params\n",
    "setup[\"voxelSize\"]=voxelSize; #voxel size\n",
    "setup[\"scale\"]=1e4; #scale for visualization\n",
    "setup[\"hierarchical\"]=true; #hierachical simualtion \n",
    "# if setup[\"hierarchical\"] is true it will assume each voxel is one node, \n",
    "# else it will assume each voxel is a detailed cuboct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.c. Fill mesh with voxels (wip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# rhino=false\n",
    "# rhinoFileName= \"./trial.stl\"\n",
    "# voxelSize=5.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.c. Draw from voxel list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "setup[\"useVoxelList\"]=false\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Boundary Conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.a. Global Settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#simulation params\n",
    "setup[\"numTimeSteps\"]=5000; #num of saved timesteps for simulation\n",
    "\n",
    "setup[\"poisson\"]=false; # account for poisson ration (only for hierarchical)\n",
    "setup[\"linear\"]=true; # linear vs non-linear\n",
    "setup[\"thermal\"]=true; #if there is change in temperature\n",
    "setup[\"globalDamping\"]=0.15; # (usually from 0.1 to 0.4)\n",
    "\n",
    "\n",
    "#visualization params\n",
    "setup[\"maxNumFiles\"]=200; #num of saved timesteps for visualization, make sure it's bigger than numTimeSteps\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.b. Materials"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#default material\n",
    "material1= Dict()\n",
    "material1[\"area\"]=voxelSize*voxelSize\n",
    "material1[\"density\"]=1e3\n",
    "material1[\"stiffness\"]=1e6\n",
    "material1[\"poissonRatio\"]=0.0\n",
    "material1[\"cTE\"]=0.0 #coefficient of thermal expansion\n",
    "\n",
    "#large bounding box for default material\n",
    "boundingBoxMaterial1=Dict()\n",
    "boundingBoxMaterial1[\"min\"]=Dict()\n",
    "boundingBoxMaterial1[\"max\"]=Dict()\n",
    "\n",
    "boundingBoxMaterial1[\"min\"][\"x\"]=-voxelSize*gridSize;\n",
    "boundingBoxMaterial1[\"min\"][\"y\"]=-voxelSize*gridSize;\n",
    "boundingBoxMaterial1[\"min\"][\"z\"]=-voxelSize*gridSize;\n",
    "\n",
    "boundingBoxMaterial1[\"max\"][\"x\"]= voxelSize*gridSize;\n",
    "boundingBoxMaterial1[\"max\"][\"y\"]= voxelSize*gridSize;\n",
    "boundingBoxMaterial1[\"max\"][\"z\"]= voxelSize*gridSize;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#second material\n",
    "material2= Dict()\n",
    "material2[\"area\"]=voxelSize*voxelSize\n",
    "material2[\"density\"]=1e3\n",
    "material2[\"stiffness\"]=1e6\n",
    "material2[\"poissonRatio\"]=0.0\n",
    "material2[\"cTE\"]=0.1 #coefficient of thermal expansion\n",
    "\n",
    "#bounding box material 2\n",
    "boundingBoxMaterial2=Dict()\n",
    "boundingBoxMaterial2[\"min\"]=Dict()\n",
    "boundingBoxMaterial2[\"max\"]=Dict()\n",
    "\n",
    "\n",
    "boundingBoxMaterial2[\"min\"][\"x\"]=0;\n",
    "boundingBoxMaterial2[\"min\"][\"y\"]=voxelSize;\n",
    "boundingBoxMaterial2[\"min\"][\"z\"]=0;\n",
    "\n",
    "boundingBoxMaterial2[\"max\"][\"x\"]= voxelSize*(latticeSizeX);\n",
    "boundingBoxMaterial2[\"max\"][\"y\"]= voxelSize*(latticeSizeY);\n",
    "boundingBoxMaterial2[\"max\"][\"z\"]= voxelSize*(latticeSizeZ);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "setup[\"materials\"]=[\n",
    "    [boundingBoxMaterial1,material1],\n",
    "    [boundingBoxMaterial2,material2]\n",
    "];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.c. Supports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "#x,y,z,rx,ry,rz (default is pinned joing i.e [false,false,false,true,true,true])\n",
    "dof=[true,true,true,true,true,true]\n",
    "\n",
    "boundingBoxSupport1=Dict()\n",
    "boundingBoxSupport1[\"min\"]=Dict()\n",
    "boundingBoxSupport1[\"max\"]=Dict()\n",
    "\n",
    "\n",
    "boundingBoxSupport1[\"min\"][\"x\"]= 0;\n",
    "boundingBoxSupport1[\"min\"][\"y\"]= 0;\n",
    "boundingBoxSupport1[\"min\"][\"z\"]= 0;\n",
    "\n",
    "boundingBoxSupport1[\"max\"][\"x\"]= voxelSize;\n",
    "boundingBoxSupport1[\"max\"][\"y\"]= voxelSize*(latticeSizeY);\n",
    "boundingBoxSupport1[\"max\"][\"z\"]= voxelSize*(latticeSizeZ);\n",
    "\n",
    "setup[\"supports\"]=[\n",
    "        [boundingBoxSupport1,dof]\n",
    "    ];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.d. Loads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.d.1 Static Loads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "load1=Dict()\n",
    "load1[\"x\"]=0.0\n",
    "load1[\"y\"]=0.0\n",
    "load1[\"z\"]=0.0\n",
    "\n",
    "boundingBoxLoad1=Dict()\n",
    "boundingBoxLoad1[\"min\"]=Dict()\n",
    "boundingBoxLoad1[\"max\"]=Dict()\n",
    "\n",
    "boundingBoxLoad1[\"min\"][\"x\"]=voxelSize*(latticeSizeX-1);\n",
    "boundingBoxLoad1[\"min\"][\"y\"]=0;\n",
    "boundingBoxLoad1[\"min\"][\"z\"]=0;\n",
    "\n",
    "boundingBoxLoad1[\"max\"][\"x\"]=voxelSize*(latticeSizeX);\n",
    "boundingBoxLoad1[\"max\"][\"y\"]=voxelSize*(latticeSizeY);\n",
    "boundingBoxLoad1[\"max\"][\"z\"]=voxelSize*(latticeSizeZ);\n",
    "\n",
    "\n",
    "setup[\"loads\"]=[\n",
    "        [boundingBoxLoad1,load1]\n",
    "    ];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.d.2 Fixed Displacements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "setup[\"fixedDisplacements\"]=[];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2.d.3 Dynamic Loads"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "updateTemperature (generic function with 1 method)"
      ]
     },
     "execution_count": 73,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function floorEnabled()\n",
    "    return false\n",
    "end\n",
    "\n",
    "function gravityEnabled()\n",
    "    return false\n",
    "end\n",
    "\n",
    "function externalDisplacement(currentTimeStep,N_position,N_fixedDisplacement)\n",
    "    return N_fixedDisplacement\n",
    "end\n",
    "\n",
    "function externalForce(currentTimeStep,N_position,N_force)\n",
    "    if currentTimeStep>1000\n",
    "        return Vector3(0,0,0)\n",
    "    else\n",
    "        return N_force\n",
    "    end\n",
    "end\n",
    "\n",
    "# if no temperature:\n",
    "# function updateTemperature(currentRestLength,currentTimeStep,mat)\n",
    "#     return currentRestLength\n",
    "# end\n",
    "\n",
    "function updateTemperature(currentRestLength,currentTimeStep,mat)\n",
    "    if currentTimeStep<1000\n",
    "        temp=-5.0*currentTimeStep/1000\n",
    "        currentRestLength=0.5*mat.L*(2.0+temp*mat.cTE)\n",
    "    elseif currentTimeStep==2500\n",
    "        temp=0\n",
    "        currentRestLength=0.5*mat.L*(2.0+temp*mat.cTE)\n",
    "    end\n",
    "    return currentRestLength\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Export setup to json"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "axialStrain (generic function with 1 method)"
      ]
     },
     "execution_count": 234,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# alternativly you can get a saved setup from an external julia file\n",
    "# include(\"./julia/examples/thermalTest.jl\") #template for multimaterial hierarchical voxels with different thermal coefficient of thermal expansion \n",
    "# include(\"./julia/examples/poissonTest.jl\") #template for hierarchical voxels with global poisson ratio\n",
    "# include(\"./julia/examples/latticeTest.jl\") #template for lattice voxel (small scale)\n",
    "# include(\"./julia/examples/latticeTest2.jl\") #template for lattice voxel (big scale with real params)\n",
    "# include(\"./julia/examples/rhinoTest.jl\") #template for importing geometry from rhino\n",
    "# include(\"./julia/examples/rhinoTestChiral.jl\") #template for importing chiral array\n",
    "include(\"./julia/examples/rover.jl\") #template for importing chiral array\n",
    "\n",
    "\n",
    "## rerun these just for sanity check for dynamic loads\n",
    "include(\"./julia/include/run.jl\") #turn setup to cuda arrays and run simulation\n",
    "include(\"./julia/include/updateEdges.jl\") #edges properties update\n",
    "include(\"./julia/include/forces.jl\") #force integration\n",
    "include(\"./julia/include/updateNodes.jl\") #nodes properties update"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2097"
      ]
     },
     "execution_count": 235,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#export prev. settings to json\n",
    "fileName=\"./json/$(name)Init.json\"\n",
    "setup1=Dict()\n",
    "setup1[\"setup\"]=setup\n",
    "stringdata = JSON.json(setup1)\n",
    "open(fileName, \"w\") do f\n",
    "        write(f, stringdata)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loaded rhino3dm.\n",
      "Success!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Process(`\u001b[4mnode\u001b[24m \u001b[4mapp1.js\u001b[24m \u001b[4mtutorial\u001b[24m`, ProcessExited(0))"
      ]
     },
     "execution_count": 236,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#run node.js to draw the gerometry using rhino3dm\n",
    "mycommand = `node app1.js $(name)`\n",
    "run(mycommand)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Run Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 237,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dt: 5.032921456503956e-6\n",
      "first timestep took 29.999862 seconds\n",
      "ran 17 nodes and 16 edges for 10000 time steps took 8.470833901 seconds\n"
     ]
    }
   ],
   "source": [
    "folderPath=\"./json/tutorial/\" # make sure this folder exists\n",
    "setupSim=getSetup(name);\n",
    "runMetavoxelGPULive!(setupSim,folderPath)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Visualize \n",
    "(only need to run it once to open the server then press stop, the server will keep running and other changes will update automatically.. will change later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Server listening on port 8080\n",
      "Open http://localhost:8080/demos/indexTutorial.html in your browser\n"
     ]
    },
    {
     "ename": "InterruptException",
     "evalue": "InterruptException:",
     "output_type": "error",
     "traceback": [
      "InterruptException:",
      "",
      "Stacktrace:",
      " [1] try_yieldto(::typeof(Base.ensure_rescheduled), ::Base.RefValue{Task}) at .\\task.jl:517",
      " [2] wait() at .\\task.jl:592",
      " [3] wait(::Base.GenericCondition{Base.AlwaysLockedST}) at .\\condition.jl:104",
      " [4] stream_wait(::Base.Process, ::Base.GenericCondition{Base.AlwaysLockedST}) at .\\stream.jl:47",
      " [5] wait at .\\process.jl:956 [inlined]",
      " [6] success at .\\process.jl:771 [inlined]",
      " [7] #run#536(::Bool, ::typeof(run), ::Cmd) at .\\process.jl:728",
      " [8] run(::Cmd) at .\\process.jl:726",
      " [9] top-level scope at In[49]:3"
     ]
    }
   ],
   "source": [
    "#run node.js to serve the indexTutorial.html for visualizarion\n",
    "mycommand = `node serve.js $(name)`\n",
    "run(mycommand)\n",
    "\n",
    "# vis 1 stable\n",
    "# http://localhost:8080/demos/indexTutorial.html\n",
    "\n",
    "# vis 2 faster for larger simulations\n",
    "# http://localhost:8080/demos/indexTutorialGraph.html\n",
    "\n",
    "# vis 3 (GPU Shaders) even faster (max 40 timesteps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.2.0",
   "language": "julia",
   "name": "julia-1.2"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}