{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Plots\n",
    "import JSON\n",
    "# using Quaternions\n",
    "using StaticArrays, Rotations\n",
    "using Distributed\n",
    "using StaticArrays, BenchmarkTools\n",
    "using Base.Threads\n",
    "using CUDAnative\n",
    "using CuArrays,CUDAdrv \n",
    "using Test\n",
    "import Base: +, * , -, ^\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "struct Vector3\n",
    "    x::Float32\n",
    "    y::Float32\n",
    "    z::Float32\n",
    "    function Vector3()\n",
    "        x=0.0\n",
    "        y=0.0\n",
    "        z=0.0\n",
    "        new(x,y,z)\n",
    "    end\n",
    "    function Vector3(x,y,z)\n",
    "       new(x,y,z)\n",
    "    end\n",
    "end\n",
    "struct Quaternion\n",
    "    x::Float32\n",
    "    y::Float32\n",
    "    z::Float32\n",
    "    w::Float32\n",
    "    function Quaternion()\n",
    "        x=0.0\n",
    "        y=0.0\n",
    "        z=0.0\n",
    "        w=1.0\n",
    "        new(x,y,z,w)\n",
    "    end\n",
    "    function Quaternion(x,y,z,w)\n",
    "        new(x,y,z,w)\n",
    "    end\n",
    "end\n",
    "struct RotationMatrix\n",
    "    te1::Float32\n",
    "    te2::Float32\n",
    "    te3::Float32\n",
    "    te4::Float32\n",
    "    te5::Float32\n",
    "    te6::Float32\n",
    "    te7::Float32\n",
    "    te8::Float32\n",
    "    te9::Float32\n",
    "    te10::Float32\n",
    "    te11::Float32\n",
    "    te12::Float32\n",
    "    te13::Float32\n",
    "    te14::Float32\n",
    "    te15::Float32\n",
    "    te16::Float32\n",
    "    function RotationMatrix()\n",
    "        te1 =0.0\n",
    "        te2 =0.0\n",
    "        te3 =0.0\n",
    "        te4 =0.0\n",
    "        te5 =0.0\n",
    "        te6 =0.0\n",
    "        te7 =0.0\n",
    "        te8 =0.0\n",
    "        te9 =0.0\n",
    "        te10=0.0\n",
    "        te11=0.0\n",
    "        te12=0.0\n",
    "        te13=0.0\n",
    "        te14=0.0\n",
    "        te15=0.0\n",
    "        te16=0.0\n",
    "        new(te1,te2,te3,te4,te5,te6,te7,te8,te9,te10,te11,te12,te13,te14,te15,te16)\n",
    "    end\n",
    "    function RotationMatrix(te1,te2,te3,te4,te5,te6,te7,te8,te9,te10,te11,te12,te13,te14,te15,te16)\n",
    "        new(te1,te2,te3,te4,te5,te6,te7,te8,te9,te10,te11,te12,te13,te14,te15,te16)\n",
    "    end\n",
    "end\n",
    "\n",
    "+(f::Vector3, g::Vector3)=Vector3(f.x+g.x , f.y+g.y,f.z+g.z )\n",
    "-(f::Vector3, g::Vector3)=Vector3(f.x-g.x , f.y-g.y,f.z-g.z )\n",
    "*(f::Vector3, g::Vector3)=Vector3(f.x*g.x , f.y*g.y,f.z*g.z )\n",
    "\n",
    "+(f::Vector3, g::Number)=Vector3(f.x+g , f.y+g,f.z+g )\n",
    "-(f::Vector3, g::Number)=Vector3(f.x-g , f.y-g,f.z-g )\n",
    "*(f::Vector3, g::Number)=Vector3(f.x*g , f.y*g,f.z*g )\n",
    "\n",
    "+(g::Vector3, f::Number)=Vector3(f.x+g , f.y+g,f.z+g )\n",
    "-(g::Vector3, f::Number)=Vector3(g-f.x , g-f.y,g-f.z )\n",
    "*(g::Vector3, f::Number)=Vector3(f.x*g , f.y*g,f.z*g )\n",
    "\n",
    "addX(f::Vector3, g::Number)=Vector3(f.x+g , f.y,f.z)\n",
    "addY(f::Vector3, g::Number)=Vector3(f.x , f.y+g,f.z )\n",
    "addZ(f::Vector3, g::Number)=Vector3(f.x , f.y,f.z+g )\n",
    "\n",
    "function normalizeVector3(f::Vector3)\n",
    "    leng=sqrt((f.x * f.x) + (f.y * f.y) + (f.z * f.z))\n",
    "    return Vector3(f.x/leng,f.y/leng,f.z/leng)\n",
    "    \n",
    "end\n",
    "\n",
    "\n",
    "function normalizeQuaternion(f::Quaternion)\n",
    "    l = sqrt((f.x * f.x) + (f.y * f.y) + (f.z * f.z)+ (f.w * f.w))\n",
    "    if l === 0 \n",
    "        qx = 0\n",
    "        qy = 0\n",
    "        qz = 0\n",
    "        qw = 1\n",
    "    else \n",
    "        l = 1 / l\n",
    "        qx = f.x * l\n",
    "        qy = f.y * l\n",
    "        qz = f.z * l\n",
    "        qw = f.w * l\n",
    "    end\n",
    "    return Quaternion(qx,qy,qz,qw)\n",
    "end\n",
    "\n",
    "function normalizeQuaternion1!(fx::Float32,fy::Float32,fz::Float32,fw::Float32)\n",
    "    l = sqrt((fx * fx) + (fy * fy) + (fz * fz)+ (fw * fw))\n",
    "    if l === 0 \n",
    "        qx = 0.0\n",
    "        qy = 0.0\n",
    "        qz = 0.0\n",
    "        qw = 1.0\n",
    "    else \n",
    "        l = 1.0 / l\n",
    "        qx = fx * l\n",
    "        qy = fy * l\n",
    "        qz = fz * l\n",
    "        qw = fw * l\n",
    "    end\n",
    "    return qx,qy,qz,qw\n",
    "end\n",
    "\n",
    "\n",
    "function dotVector3(f::Vector3, g::Vector3)\n",
    "    return (f.x * g.x) + (f.y * g.y) + (f.z * g.z)\n",
    "end\n",
    "\n",
    "function Base.show(io::IO, v::Vector3)\n",
    "    print(io, \"x:$(v.x), y:$(v.y), z:$(v.z)\")\n",
    "end\n",
    "\n",
    "function Base.show(io::IO, v::Quaternion)\n",
    "    print(io, \"x:$(v.x), y:$(v.y), z:$(v.z), w:$(v.z)\")\n",
    "end\n",
    "\n",
    "Base.Broadcast.broadcastable(q::Vector3) = Ref(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "mutable struct Node\n",
    "    id::Int32\n",
    "    position::Vector3\n",
    "    restrained::Bool\n",
    "    displacement::Vector3\n",
    "    angle::Vector3\n",
    "    force::Vector3\n",
    "    currPosition::Vector3\n",
    "    orient::Quaternion\n",
    "    linMom::Vector3\n",
    "    angMom::Vector3\n",
    "    intForce::Vector3\n",
    "    intMoment::Vector3\n",
    "    moment::Vector3\n",
    "    function Node()\n",
    "        id=0\n",
    "        position=Vector3()\n",
    "        restrained=false\n",
    "        displacement=Vector3()\n",
    "        angle=Vector3()\n",
    "        force=Vector3()\n",
    "        currPosition=Vector3()\n",
    "        orient=Quaternion()\n",
    "        linMom=Vector3()\n",
    "        angMom=Vector3()\n",
    "        intForce=Vector3()\n",
    "        intMoment=Vector3()\n",
    "        moment=Vector3()\n",
    "        new(id,position,restrained,displacement,angle,force,currPosition,orient,linMom,angMom,intForce,intMoment,moment)\n",
    "    end\n",
    "end\n",
    "\n",
    "\n",
    "mutable struct Edge\n",
    "    id::Int32\n",
    "    source::Int32 #change to Int32\n",
    "    target::Int32\n",
    "    area::Float32\n",
    "    density::Float32\n",
    "    stiffness::Float32\n",
    "    stress::Float32\n",
    "    axis::Vector3\n",
    "    currentRestLength::Float32\n",
    "    pos2::Vector3\n",
    "    angle1v::Vector3\n",
    "    angle2v::Vector3\n",
    "    angle1::Quaternion\n",
    "    angle2::Quaternion\n",
    "    currentTransverseStrainSum::Float32\n",
    "    ## add pos node and negative node\n",
    "    ## add memory cuda??\n",
    "    function Edge()\n",
    "        id=0\n",
    "        source=0\n",
    "        target=0\n",
    "        area=0.0\n",
    "        density=0.0\n",
    "        stiffness=0.0\n",
    "        stress=0.0\n",
    "        axis=Vector3(1.0,0.0,0.0)\n",
    "        currentRestLength=0.0\n",
    "        pos2=Vector3()\n",
    "        angle1v=Vector3()\n",
    "        angle2v=Vector3()\n",
    "        angle1=Quaternion()\n",
    "        angle2=Quaternion()\n",
    "        currentTransverseStrainSum=0.0\n",
    "        \n",
    "        new(id,source,target,area,density,stiffness,stress,axis,currentRestLength,pos2,angle1v,angle2v,angle1,angle2,currentTransverseStrainSum)\n",
    "    end\n",
    "end\n",
    "\n",
    "function Base.show(io::IO, v::Node)\n",
    "    print(io, \"node:$(v.id), position:($(v.position)), restrained:$(v.restrained)\")\n",
    "end\n",
    "\n",
    "function Base.show(io::IO, v::Edge)\n",
    "    print(io, \"edge:$(v.id), source:$(v.source), target:$(v.target), stress:$(v.stress), axis:($(v.axis))\")\n",
    "end\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element CuArray{Vector3,1,Nothing}:\n",
       " x:2.0, y:2.0, z:2.0\n",
       " x:2.0, y:2.0, z:2.0\n",
       " x:2.0, y:2.0, z:2.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# allowscalar(false)\n",
    "m=Vector3(1.0,1.0,1.0)\n",
    "mg=CuArray([m,m,m])\n",
    "mm=mg.+m\n",
    "# broadcast(+, mg, m)\n",
    "# broadcast(addX, mg, 1)\n",
    "\n",
    "# m=node\n",
    "# println(m)\n",
    "# m=Node()\n",
    "# mg=CuArray([m,m,m])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RotationMatrix(0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0, 0.0f0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mat=RotationMatrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{String,Any} with 8 entries:\n",
       "  \"nodes\"        => Any[Dict{String,Any}(\"degrees_of_freedom\"=>Any[0, 1, 2, 3, …\n",
       "  \"voxelSize\"    => 5\n",
       "  \"numTimeSteps\" => 100\n",
       "  \"hierarchical\"   => false\n",
       "  \"animation\"    => Dict{String,Any}(\"speed\"=>3,\"exaggeration\"=>2000,\"showDispla…\n",
       "  \"viz\"          => Dict{String,Any}(\"colorMap\"=>0,\"colorMaps\"=>Any[Any[Any[0, …\n",
       "  \"edges\"        => Any[Dict{String,Any}(\"source\"=>0,\"area\"=>1,\"density\"=>0.028…\n",
       "  \"ndofs\"        => 1800"
      ]
     },
     "execution_count": 195,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "###############################################################################################\n",
    "setup = Dict()\n",
    "open(\"../json/setupTestUni4.json\", \"r\") do f\n",
    "    global setup\n",
    "    dicttxt = String(read(f))  # file information to string\n",
    "    setup=JSON.parse(dicttxt)  # parse and transform data\n",
    "end\n",
    "\n",
    "setup=setup[\"setup\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "initialize (generic function with 1 method)"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function initialize(setup)\n",
    "\tnodes      = setup[\"nodes\"]\n",
    "    edges      = setup[\"edges\"]\n",
    "    \n",
    "    i=1\n",
    "\t# pre-calculate current position\n",
    "\tfor node in nodes\n",
    "        # element=parse(Int,node[\"id\"][2:end])\n",
    "        N_position[i]=Vector3(node[\"position\"][\"x\"],node[\"position\"][\"y\"],node[\"position\"][\"z\"])\n",
    "        N_restrained[i]=node[\"restrained_degrees_of_freedom\"][1] ## todo later consider other degrees of freedom\n",
    "        N_displacement[i]=Vector3(node[\"displacement\"][\"x\"],node[\"displacement\"][\"y\"],node[\"displacement\"][\"z\"])\n",
    "        N_angle[i]=Vector3(node[\"angle\"][\"x\"],node[\"angle\"][\"y\"],node[\"angle\"][\"z\"])\n",
    "        N_force[i]=Vector3(node[\"force\"][\"x\"],node[\"force\"][\"y\"],node[\"force\"][\"z\"])\n",
    "        N_currPosition[i]=Vector3(node[\"position\"][\"x\"],node[\"position\"][\"y\"],node[\"position\"][\"z\"])\n",
    "\n",
    "        # for dynamic simulations\n",
    "        # append!(N_posTimeSteps,[[]])\n",
    "        # append!(N_angTimeSteps,[[]])\n",
    "        \n",
    "        i=i+1\n",
    "\tend \n",
    "    \n",
    "    i=1\n",
    "\t# pre-calculate the axis\n",
    "\tfor edge in edges\n",
    "        # element=parse(Int,edge[\"id\"][2:end])\n",
    "        \n",
    "        # find the nodes that the lements connects\n",
    "        fromNode = nodes[edge[\"source\"]+1]\n",
    "        toNode = nodes[edge[\"target\"]+1]\n",
    "\n",
    "        \n",
    "        node1 = [fromNode[\"position\"][\"x\"] fromNode[\"position\"][\"y\"] fromNode[\"position\"][\"z\"]]\n",
    "        node2 = [toNode[\"position\"][\"x\"] toNode[\"position\"][\"y\"] toNode[\"position\"][\"z\"]]\n",
    "        \n",
    "        length=norm(node2-node1)\n",
    "        axis=normalize(collect(Iterators.flatten(node2-node1)))\n",
    "        \n",
    "        E_source[i]=edge[\"source\"]+1\n",
    "        E_target[i]=edge[\"target\"]+1\n",
    "        E_area[i]=edge[\"area\"]\n",
    "        E_density[i]=edge[\"density\"]\n",
    "        E_stiffness[i]=edge[\"stiffness\"]\n",
    "        E_axis[i]=Vector3(axis[1],axis[2],axis[3])\n",
    "        E_currentRestLength[i]=length\n",
    "        \n",
    "        N_edgeID[E_source[i],N_currEdge[E_source[i]]]=i\n",
    "        N_edgeFirst[E_source[i],N_currEdge[E_source[i]]]=true\n",
    "        N_currEdge[E_source[i]]+=1\n",
    "        \n",
    "        N_edgeID[E_target[i],N_currEdge[E_target[i]]]=i\n",
    "        N_edgeFirst[E_target[i],N_currEdge[E_target[i]]]=false\n",
    "        N_currEdge[E_target[i]]+=1\n",
    "       \n",
    "        \n",
    "        # for dynamic simulations\n",
    "        # append!(E_stressTimeSteps,[[]])\n",
    "        \n",
    "        i=i+1\n",
    "\tend \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "960-element CuArray{Vector3,1,Nothing}:\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " ⋮                  \n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0\n",
       " x:0.0, y:0.0, z:0.0"
      ]
     },
     "execution_count": 217,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "########\n",
    "voxCount=0\n",
    "linkCount=0\n",
    "nodes      = setup[\"nodes\"]\n",
    "edges      = setup[\"edges\"]\n",
    "voxCount=size(nodes)[1]\n",
    "linkCount=size(edges)[1]\n",
    "strain =0 #todooo moveeee\n",
    "maxNumEdges=10\n",
    "\n",
    "\n",
    "# Nodes=fill(Node(),voxCount)\n",
    "# Edges=fill(Edge(),linkCount)\n",
    "\n",
    "########\n",
    "voxCount=0\n",
    "linkCount=0\n",
    "nodes      = setup[\"nodes\"]\n",
    "edges      = setup[\"edges\"]\n",
    "voxCount=size(nodes)[1]\n",
    "linkCount=size(edges)[1]\n",
    "strain =0 #todooo moveeee\n",
    "\n",
    "############# nodes\n",
    "N_position=fill(Vector3(),voxCount)\n",
    "N_restrained=zeros(Bool, voxCount)\n",
    "N_displacement=fill(Vector3(),voxCount)\n",
    "N_angle=fill(Vector3(),voxCount)\n",
    "N_currPosition=fill(Vector3(),voxCount)\n",
    "N_linMom=fill(Vector3(),voxCount)\n",
    "N_angMom=fill(Vector3(),voxCount)\n",
    "N_intForce=fill(Vector3(),voxCount)\n",
    "N_intMoment=fill(Vector3(),voxCount)\n",
    "N_moment=fill(Vector3(),voxCount)\n",
    "# N_posTimeSteps=[]\n",
    "# N_angTimeSteps=[]\n",
    "N_force=fill(Vector3(),voxCount)\n",
    "N_orient=fill(Quaternion(),voxCount)\n",
    "N_edgeID=fill(-1,(voxCount,maxNumEdges))\n",
    "N_edgeFirst=fill(true,(voxCount,maxNumEdges))\n",
    "N_currEdge=fill(1,voxCount)\n",
    "\n",
    "############# edges\n",
    "E_source=fill(0,linkCount)\n",
    "E_target=fill(0,linkCount)\n",
    "E_area=fill(0.0f0,linkCount)\n",
    "E_density=fill(0.0f0,linkCount)\n",
    "E_stiffness=fill(0.0f0,linkCount)\n",
    "E_stress=fill(0.0f0,linkCount)\n",
    "E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)\n",
    "E_currentRestLength=fill(0.0f0,linkCount)\n",
    "E_pos2=fill(Vector3(),linkCount)\n",
    "E_angle1v=fill(Vector3(),linkCount)\n",
    "E_angle2v=fill(Vector3(),linkCount)\n",
    "E_angle1=fill(Quaternion(),linkCount)\n",
    "E_angle2=fill(Quaternion(),linkCount)\n",
    "\n",
    "E_intForce1=fill(Vector3(),linkCount)\n",
    "E_intMoment1=fill(Vector3(),linkCount) \n",
    "\n",
    "E_intForce2=fill(Vector3(),linkCount)\n",
    "E_intMoment2=fill(Vector3(),linkCount) \n",
    "\n",
    "E_currentTransverseStrainSum=fill(0.0f0,linkCount)# TODO remove ot incorporate\n",
    "# E_stressTimeSteps=[]\n",
    "\n",
    "\n",
    "#################################################################\n",
    "initialize(setup)\n",
    "#################################################################\n",
    "\n",
    "########################## turn to cuda arrays\n",
    "############# nodes\n",
    "N_positionGPU=    CuArray(N_position)  \n",
    "N_restrainedGPU=  CuArray(N_restrained)  \n",
    "N_displacementGPU=CuArray(N_displacement)   \n",
    "N_angleGPU=       CuArray(N_angle)  \n",
    "N_currPositionGPU=CuArray(N_currPosition)    \n",
    "N_linMomGPU=      CuArray(N_linMom)    \n",
    "N_angMomGPU=      CuArray(N_angMom)  \n",
    "N_intForceGPU=    CuArray(N_intForce) \n",
    "N_intMomentGPU=   CuArray(N_intMoment)   \n",
    "N_momentGPU=      CuArray(N_moment)\n",
    "N_forceGPU=       CuArray(N_force)  \n",
    "N_orientGPU=      CuArray(N_orient)\n",
    "N_edgeIDGPU=      CuArray(N_edgeID)\n",
    "N_edgeFirstGPU=   CuArray(N_edgeFirst)\n",
    " \n",
    "\n",
    "############# edges\n",
    "E_sourceGPU=                    CuArray(E_source)   \n",
    "E_targetGPU=                    CuArray(E_target)\n",
    "E_areaGPU=                      CuArray(E_area)                             \n",
    "E_densityGPU=                   CuArray(E_density)\n",
    "E_stiffnessGPU=                 CuArray(E_stiffness)\n",
    "E_stressGPU=                    CuArray(E_stress)\n",
    "E_axisGPU=                      CuArray(E_axis)          \n",
    "E_currentRestLengthGPU=         CuArray(E_currentRestLength)\n",
    "E_pos2GPU=                      CuArray(E_pos2)\n",
    "E_angle1vGPU=                   CuArray(E_angle1v)\n",
    "E_angle2vGPU=                   CuArray(E_angle2v)\n",
    "E_angle1GPU=                    CuArray(E_angle1)\n",
    "E_angle2GPU=                    CuArray(E_angle2)\n",
    "E_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)\n",
    "E_intForce1GPU=      CuArray(E_intForce1) \n",
    "E_intMoment1GPU=     CuArray(E_intMoment1)  \n",
    "E_intForce2GPU=      CuArray(E_intForce2) \n",
    "E_intMoment2GPU=     CuArray(E_intMoment2) \n",
    "# E_stressTimeSteps=[]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "voxCount 300 linkCount 960\n"
     ]
    }
   ],
   "source": [
    "N,M=size(N_edgeIDGPU)\n",
    "println(\"voxCount $voxCount linkCount $linkCount\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(nothing, 0.0904174, 3580399, 0.0, Base.GC_Diff(3580399, 0, 0, 61035, 0, 0, 0, 0, 0))"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function simulateParallel(numTimeSteps,dt,static=true)\n",
    "#     initialize(setup)\n",
    "    \n",
    "    for i in 1:numTimeSteps\n",
    "#         println(\"Timestep:\",i)\n",
    "        doTimeStep(dt,static,i)\n",
    "    end\n",
    "end\n",
    "dt=0.0251646\n",
    "numTimeSteps=100\n",
    "@timed simulateParallel(numTimeSteps,dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "244018"
      ]
     },
     "execution_count": 200,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function updateDataAndSave(setup,fileName)\n",
    "    nodes      = setup[\"nodes\"]\n",
    "    edges      = setup[\"edges\"]\n",
    "    \n",
    "    setup[\"animation\"][\"showDisplacement\"]=false\n",
    "    voxCount=size(nodes)[1]\n",
    "    linkCount=size(edges)[1]\n",
    "    \n",
    "    N_displacement=Array(N_displacementGPU)\n",
    "    N_angle=Array(N_angleGPU)\n",
    "    E_stress=Array(E_stressGPU)\n",
    "    \n",
    "    setup[\"viz\"][\"maxStress\"]=maximum(E_stress)\n",
    "    setup[\"viz\"][\"minStress\"]=minimum(E_stress)  \n",
    "\n",
    "    i=1\n",
    "\tfor edge in edges\n",
    "        edge[\"stress\"]=E_stress[i]\n",
    "        i=i+1\n",
    "\n",
    "    end\n",
    "    \n",
    " \n",
    "    i=1          \n",
    "\tfor node in nodes\n",
    "        node[\"displacement\"][\"x\"]=N_displacement[i].x\n",
    "        node[\"displacement\"][\"y\"]=N_displacement[i].y\n",
    "        node[\"displacement\"][\"z\"]=N_displacement[i].z\n",
    "        \n",
    "        node[\"angle\"][\"x\"]=N_angle[i].x\n",
    "        node[\"angle\"][\"y\"]=N_angle[i].y\n",
    "        node[\"angle\"][\"z\"]=N_angle[i].z\n",
    "        i=i+1\n",
    "\n",
    "    end\n",
    "    \n",
    "    # pass data as a json string (how it shall be displayed in a file)\n",
    "    stringdata = JSON.json(setup)\n",
    "\n",
    "    # write the file with the stringdata variable information\n",
    "    open(fileName, \"w\") do f\n",
    "            write(f, stringdata)\n",
    "         end\n",
    "    \n",
    "end\n",
    "updateDataAndSave(setup,\"../json/trialJuliaParallelGPU.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function addScalar(vec,sc)\n",
    "    return vec+Vector3(sc,sc,sc)\n",
    "end\n",
    "\n",
    "function updateEdges1!(E_source, E_target, E_stress,E_intForce1,E_intForce2,N_intForce)\n",
    "    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x\n",
    "    stride = blockDim().x * gridDim().x\n",
    "    ## @cuprintln(\"thread $index, block $stride\")\n",
    "    for i = index:stride:length(E_source)\n",
    "        @inbounds E_stress[i]=10.0\n",
    "#         @inbounds s=E_source[i]\n",
    "#         @inbounds t=E_target[i]\n",
    "        #N_intForce[s]=N_intForce[s]+Vector3(0.1,0.1,0.1)\n",
    "        #N_intForce[t]=N_intForce[t]+Vector3(0.1,0.1,0.1)\n",
    "        @inbounds E_intForce1[i]=Vector3(0.1,0.1,0.1)\n",
    "        @inbounds E_intForce2[i]=Vector3(0.1,0.1,0.1)\n",
    "        #@inbounds  N_positionVGPU[i]=ifelse(N_restrainedGPU[i], N_positionVGPU[i], N_positionVGPU[i] +vecGPU[i] )\n",
    "    end\n",
    "    return\n",
    "end\n",
    "\n",
    "\n",
    "function bench_updateEdges!(E_source, E_target, E_stress,E_intForce1,E_intForce2,N_intForce)\n",
    "    N=length(E_source)\n",
    "    numblocks = ceil(Int, N/256)\n",
    "    CuArrays.@sync begin\n",
    "        @cuda threads=256 blocks=numblocks updateEdges1!(E_source, E_target, E_stress,E_intForce1,E_intForce2,N_intForce)\n",
    "    end\n",
    "end\n",
    "@btime bench_updateEdges!(E_sourceGPU, E_targetGPU, E_stressGPU,E_intForce1GPU,E_intForce2GPU,N_intForceGPU)\n",
    "# @timed bench_updateEdges!(E_sourceGPU, E_targetGPU, E_stressGPU,E_intForce1GPU,E_intForce2GPU,N_intForceGPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function updateNodes1!(N_position, N_restrained,N_intForce, N_edgeID, N_edgeFirst, E_intForce1, E_intForce2)\n",
    "    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x\n",
    "    stride = blockDim().x * gridDim().x\n",
    "    ## @cuprintln(\"thread $index, block $stride\")\n",
    "    N,M=size(N_edgeID)\n",
    "    for i = index:stride:N\n",
    "        for j in 1:M\n",
    "            @inbounds if (N_edgeID[i,j]!=-1)\n",
    "                #@cuprintln(\"i $i, j $j, N_edgeID[i,j] $temp\")\n",
    "                @inbounds N_intForce[i]=ifelse(N_edgeFirst[i,j], N_intForce[i]+E_intForce1[N_edgeID[i,j]], N_intForce[i]+E_intForce2[N_edgeID[i,j]] ) \n",
    "            end\n",
    "        end\n",
    "        @inbounds  N_position[i]=ifelse(N_restrained[i], N_position[i], N_position[i] + N_position[i]*N_intForce[i] )\n",
    "    end\n",
    "    return\n",
    "end\n",
    "\n",
    "\n",
    "function bench_updateNodes!(N_position, N_restrained, N_intForce, N_edgeID, N_edgeFirst, E_intForce1, E_intForce2)\n",
    "    N=length(N_position)\n",
    "    numblocks = ceil(Int, N/256)\n",
    "    CuArrays.@sync begin\n",
    "        @cuda threads=256 blocks=numblocks updateNodes1!(N_position, N_restrained,N_intForce, N_edgeID, N_edgeFirst, E_intForce1, E_intForce2)\n",
    "    end\n",
    "end\n",
    "@btime bench_updateNodes!(N_positionGPU, N_restrainedGPU, N_intForceGPU, N_edgeIDGPU, N_edgeFirstGPU, E_intForce1GPU, E_intForce2GPU)\n",
    "# @timed bench_updateNodes!(N_positionGPU, N_restrainedGPU, N_intForceGPU, N_edgeIDGPU, N_edgeFirstGPU, E_intForce1GPU, E_intForce2GPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# N_positionGPU N_restrainedGPU N_intForceGPU N_edgeIDGPU N_edgeFirstGPU E_intForce1 E_intForce2\n",
    "function addVectors(f,g)\n",
    "    return Vector3(f.x+g.x , f.y+g.y,f.z+g.z )\n",
    "end\n",
    "\n",
    "function updateIntForces1!(N_intForce, N_edgeID, N_edgeFirst, E_intForce1, E_intForce2)\n",
    "    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x\n",
    "    stride = blockDim().x * gridDim().x\n",
    "    ## @cuprintln(\"thread $index, block $stride\")\n",
    "    N,M=size(N_edgeID)\n",
    "    for i = index:stride:N\n",
    "        for j in 1:M\n",
    "            temp=N_edgeID[i,j]\n",
    "            @inbounds if (N_edgeID[i,j]!=-1)\n",
    "                #@cuprintln(\"i $i, j $j, N_edgeID[i,j] $temp\")\n",
    "                @inbounds N_intForce[i]=ifelse(N_edgeFirst[i,j], N_intForce[i]+E_intForce1[N_edgeID[i,j]], N_intForce[i]+E_intForce2[N_edgeID[i,j]] ) \n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    return\n",
    "end\n",
    "\n",
    "\n",
    "function bench_updateIntForces!(N_intForce, N_edgeID, N_edgeFirst, E_intForce1, E_intForce2)\n",
    "    N=length(N_intForce)\n",
    "    numblocks = ceil(Int, N/256)\n",
    "    CuArrays.@sync begin\n",
    "        @cuda threads=256 blocks=numblocks updateIntForces1!(N_intForce, N_edgeID, N_edgeFirst, E_intForce1, E_intForce2)\n",
    "    end\n",
    "end\n",
    "# @btime bench_updateIntForces!(N_intForceGPU, N_edgeIDGPU, N_edgeFirstGPU, E_intForce1GPU, E_intForce2GPU)\n",
    "@timed bench_updateIntForces!(N_intForceGPU, N_edgeIDGPU, N_edgeFirstGPU, E_intForce1GPU, E_intForce2GPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "doTimeStep (generic function with 2 methods)"
      ]
     },
     "execution_count": 180,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function doTimeStep(dt,static,currentTimeStep)\n",
    "    # update forces: go through edges, get currentposition from nodes, calc pos2 and update stresses and interior forces of nodes\n",
    "    run_updateEdges!(\n",
    "        E_sourceGPU, \n",
    "        E_targetGPU,\n",
    "        E_areaGPU,\n",
    "        E_densityGPU,\n",
    "        E_stiffnessGPU,\n",
    "        E_stressGPU,\n",
    "        E_axisGPU,\n",
    "        E_currentRestLengthGPU,\n",
    "        E_pos2GPU,\n",
    "        E_angle1vGPU,\n",
    "        E_angle2vGPU,\n",
    "        E_angle1GPU,\n",
    "        E_angle2GPU,\n",
    "        E_intForce1GPU,\n",
    "        E_intMoment1GPU,\n",
    "        E_intForce2GPU,\n",
    "        E_intMoment2GPU,\n",
    "        N_currPositionGPU,\n",
    "        N_orientGPU)\n",
    "    \n",
    "    # update forces: go through nodes and update interior force (according to int forces from edges), integrate and update currpos\n",
    "    run_updateNodes!(dt,currentTimeStep,\n",
    "        N_positionGPU, \n",
    "        N_restrainedGPU,\n",
    "        N_displacementGPU,\n",
    "        N_angleGPU,\n",
    "        N_currPositionGPU,\n",
    "        N_linMomGPU,\n",
    "        N_angMomGPU,\n",
    "        N_intForceGPU,\n",
    "        N_intMomentGPU,\n",
    "        N_forceGPU,\n",
    "        N_momentGPU,\n",
    "        N_orientGPU,\n",
    "        N_edgeIDGPU, \n",
    "        N_edgeFirstGPU, \n",
    "        E_intForce1GPU,\n",
    "        E_intMoment1GPU,\n",
    "        E_intForce2GPU,\n",
    "        E_intMoment2GPU)\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "run_updateEdges! (generic function with 1 method)"
      ]
     },
     "execution_count": 181,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function updateEdges!(E_source,E_target,E_area,E_density,E_stiffness,E_stress,E_axis,E_currentRestLength,\n",
    "        E_pos2,E_angle1v,E_angle2v,E_angle1,E_angle2,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,N_currPosition,N_orient)\n",
    "    \n",
    "    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x\n",
    "    stride = blockDim().x * gridDim().x\n",
    "    ## @cuprintln(\"thread $index, block $stride\")\n",
    "    N=length(E_source)\n",
    "    for i = index:stride:N\n",
    "        \n",
    "        @inbounds pVNeg=N_currPosition[E_source[i]]\n",
    "        @inbounds pVPos=N_currPosition[E_target[i]]\n",
    "        \n",
    "        @inbounds oVNeg=N_orient[E_source[i]]\n",
    "        @inbounds oVPos=N_orient[E_target[i]]\n",
    "        \n",
    "        @inbounds oldPos2=E_pos2[i] #?copy?\n",
    "        @inbounds oldAngle1v = E_angle1v[i]\n",
    "        @inbounds oldAngle2v =  E_angle2v[i]# remember the positions/angles from last timestep to calculate velocity\n",
    "        \n",
    "#         E_pos2[i],E_angle1v[i],E_angle2v[i],E_angle1[i],E_angle2[i],\n",
    "        @inbounds E_pos2[i],E_angle1v[i],E_angle2v[i],E_angle1[i],E_angle2[i],totalRot= orientLink!(E_currentRestLength[i],pVNeg,pVPos,oVNeg,oVPos,E_axis[i])\n",
    "        @inbounds m=(E_pos2[i].x/E_currentRestLength[i])\n",
    "        @inbounds _stress=updateStrain(m,E_stiffness[i])\n",
    "        @inbounds E_stress[i]=_stress\n",
    "#         @cuprintln(_stress)\n",
    "        \n",
    "        @inbounds l   = E_currentRestLength[i]\n",
    "        @inbounds E = E_stiffness[i]\n",
    "        \n",
    "        nu=0\n",
    "        L = 5.0\n",
    "        a1 = E*L # EA/L : Units of N/m\n",
    "        a2 = E * L*L*L / (12.0*(1+nu)) # GJ/L : Units of N-m\n",
    "        b1 = E*L # 12EI/L^3 : Units of N/m\n",
    "        b2 = E*L*L/2.0 # 6EI/L^2 : Units of N (or N-m/m: torque related to linear distance)\n",
    "        b3 = E*L*L*L/6.0 # 2EI/L : Units of N-m\n",
    "        @inbounds currentTransverseArea=E_area[i]\n",
    "        \n",
    "        x=(_stress*currentTransverseArea)\n",
    "        @inbounds y=(b1*E_pos2[i].y-b2*(E_angle1v[i].z + E_angle2v[i].z))\n",
    "        @inbounds z=(b1*E_pos2[i].z + b2*(E_angle1v[i].y + E_angle2v[i].y))\n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    "        \n",
    "        # Use Curstress instead of -a1*Pos2.x to account for non-linear deformation \n",
    "        forceNeg = Vector3(x,y,z)\n",
    "        # forceNeg = Vector3((_stress*currentTransverseArea),(b1*E_pos2[edge].y-b2*(E_angle1v[i].z + E_angle2v[i].z)),(b1*E_pos2[i].z + b2*(E_angle1v[i].y + E_angle2v[i].y))) # Use Curstress instead of -a1*Pos2.x to account for non-linear deformation \n",
    "        forcePos = Vector3(-x,-y,-z)\n",
    "        \n",
    "        @inbounds x= (a2*(E_angle2v[i].x-E_angle1v[i].x))\n",
    "        @inbounds y= (-b2*E_pos2[i].z-b3*(2.0*E_angle1v[i].y+E_angle2v[i].y))\n",
    "        @inbounds z=(b2*E_pos2[i].y - b3*(2.0*E_angle1v[i].z + E_angle2v[i].z))  \n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    "        momentNeg = Vector3(x,y,z)\n",
    "        # momentNeg = Vector3((a2*(E_angle2v[i].x-E_angle1v[i].x)) , (-b2*E_pos2[i].z-b3*(2.0*E_angle1v[i].y+E_angle2v[i].y)),(b2*E_pos2[i].y - b3*(2.0*E_angle1v[i].z + E_angle2v[i].z)))\n",
    "\n",
    "        @inbounds x= (a2*(E_angle1v[i].x-E_angle2v[i].x))\n",
    "        @inbounds y= (-b2*E_pos2[i].z- b3*(E_angle1v[i].y+2.0*E_angle2v[i].y))\n",
    "        @inbounds z=(b2*E_pos2[i].y - b3*(E_angle1v[i].z + 2.0*E_angle2v[i].z))\n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    "        momentPos = Vector3(x,y,z)\n",
    "#         momentPos = Vector3((a2*(E_angle1v[i].x-E_angle2v[i].x)) , (-b2*E_pos2[i].z- b3*(E_angle1v[i].y+2.0*E_angle2v[i].y)),(b2*E_pos2[i].y - b3*(E_angle1v[i].z + 2.0*E_angle2v[i].z)))\n",
    "        smallAngle=false\n",
    "        if !smallAngle # ?? check\n",
    "            @inbounds forceNeg = RotateVec3DInv(E_angle1[i],forceNeg)\n",
    "            @inbounds momentNeg = RotateVec3DInv(E_angle1[i],momentNeg)\n",
    "        end\n",
    "        \n",
    "        @inbounds forcePos = RotateVec3DInv(E_angle2[i],forcePos)\n",
    "        @inbounds momentPos = RotateVec3DInv(E_angle2[i],momentPos)\n",
    "\n",
    "        # println(momentPos)\n",
    "\n",
    "        @inbounds forceNeg =toAxisOriginalVector3(forceNeg,E_axis[i])\n",
    "        @inbounds forcePos =toAxisOriginalVector3(forcePos,E_axis[i])\n",
    "\n",
    "        @inbounds momentNeg=toAxisOriginalQuat(momentNeg,E_axis[i])# TODOO CHECKKKKKK\n",
    "        @inbounds momentPos=toAxisOriginalQuat(momentPos,E_axis[i])\n",
    "\n",
    "        # println(momentPos[2],\" \",momentPos[3],\" \",momentPos[4],\" \",momentPos[1],\" \")\n",
    "\n",
    "        @inbounds E_intForce1[i] =forceNeg\n",
    "        @inbounds E_intForce2[i] =forcePos\n",
    "        #@inbounds E_intForce1[i] =E_intForce1[i]+forceNeg\n",
    "        #@inbounds E_intForce2[i] =E_intForce2[i]+forcePos\n",
    "        \n",
    "        \n",
    "        #@inbounds x= E_intMoment1[i].x+momentNeg.x\n",
    "        #@inbounds y= E_intMoment1[i].y+momentNeg.y\n",
    "        #@inbounds z= E_intMoment1[i].z+momentNeg.z\n",
    "        @inbounds x= momentNeg.x\n",
    "        @inbounds y= momentNeg.y\n",
    "        @inbounds z= momentNeg.z  \n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    "        \n",
    "        @inbounds E_intMoment1[i]=Vector3(x,y,z)\n",
    "        \n",
    "        #@inbounds x= E_intMoment2[i].x+momentNeg.x\n",
    "        #@inbounds y= E_intMoment2[i].y+momentNeg.y\n",
    "        #@inbounds z= E_intMoment2[i].z+momentNeg.z\n",
    "        @inbounds x= momentNeg.x\n",
    "        @inbounds y= momentNeg.y\n",
    "        @inbounds z= momentNeg.z\n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    "        \n",
    "        @inbounds E_intMoment2[i]=Vector3(x,y,z)\n",
    "        \n",
    "        \n",
    "    end\n",
    "    return\n",
    "end\n",
    "\n",
    "# function trial(f,g)\n",
    "# #     xaxis=Vector3(1,0,0)\n",
    "# #     xaxis=Vector{Float64}(1.0, 3)\n",
    "#     m=normalizeVector3(f)\n",
    "#     return m\n",
    "# end\n",
    "\n",
    "\n",
    "function run_updateEdges!(E_source,E_target,E_area,E_density,E_stiffness,E_stress,E_axis,E_currentRestLength,\n",
    "        E_pos2,E_angle1v,E_angle2v,E_angle1,E_angle2,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,N_currPosition,N_orient)\n",
    "    N=length(E_source)\n",
    "    numblocks = ceil(Int, N/256)\n",
    "    CuArrays.@sync begin\n",
    "        @cuda threads=256 blocks=numblocks updateEdges!(E_source,E_target,E_area,E_density,E_stiffness,E_stress,E_axis,E_currentRestLength,\n",
    "            E_pos2,E_angle1v,E_angle2v,E_angle1,E_angle2,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,N_currPosition,N_orient)\n",
    "    end\n",
    "end\n",
    "\n",
    "\n",
    "# @btime run_updateEdges!(E_sourceGPU, E_targetGPU, E_stressGPU,E_intForce1GPU,E_intForce2GPU,N_intForceGPU)\n",
    "# @timed run_updateEdges!(\n",
    "#         E_sourceGPU, \n",
    "#         E_targetGPU,\n",
    "#         E_areaGPU,\n",
    "#         E_densityGPU,\n",
    "#         E_stiffnessGPU,\n",
    "#         E_stressGPU,\n",
    "#         E_axisGPU,\n",
    "#         E_currentRestLengthGPU,\n",
    "#         E_pos2GPU,\n",
    "#         E_angle1vGPU,\n",
    "#         E_angle2vGPU,\n",
    "#         E_angle1GPU,\n",
    "#         E_angle2GPU,\n",
    "#         E_intForce1GPU,\n",
    "#         E_intMoment1GPU,\n",
    "#         E_intForce2GPU,\n",
    "#         E_intMoment2GPU,\n",
    "#         N_currPositionGPU,\n",
    "#         N_orientGPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "run_updateNodes! (generic function with 3 methods)"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function updateNodes!(dt,currentTimeStep,N_position, N_restrained,N_displacement,N_angle,N_currPosition,N_linMom,N_angMom,\n",
    "        N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\n",
    "    \n",
    "    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x\n",
    "    stride = blockDim().x * gridDim().x\n",
    "    ## @cuprintln(\"thread $index, block $stride\")\n",
    "    N,M=size(N_edgeID)\n",
    "    for i = index:stride:N\n",
    "        @inbounds if N_restrained[i]\n",
    "            return\n",
    "        else\n",
    "            for j in 1:M\n",
    "                temp=N_edgeID[i,j]\n",
    "                @inbounds if (N_edgeID[i,j]!=-1)\n",
    "                    #@cuprintln(\"i $i, j $j, N_edgeID[i,j] $temp\")\n",
    "                    @inbounds N_intForce[i]=ifelse(N_edgeFirst[i,j], N_intForce[i]+E_intForce1[N_edgeID[i,j]], N_intForce[i]+E_intForce2[N_edgeID[i,j]] )\n",
    "                    @inbounds N_intMoment[i]=ifelse(N_edgeFirst[i,j], N_intMoment[i]+E_intMoment1[N_edgeID[i,j]], N_intMoment[i]+E_intMoment2[N_edgeID[i,j]] )\n",
    "                end\n",
    "            end\n",
    "            @inbounds curForce = force(N_intForce[i],N_orient[i],N_force[i],true,currentTimeStep)\n",
    "            \n",
    "            @inbounds N_force[i]=Vector3(0,0,0) ##????\n",
    "            \n",
    "            @inbounds N_intForce[i]=Vector3(0,0,0)\n",
    "            \n",
    "            @inbounds N_linMom[i]=N_linMom[i]+curForce*Vector3(dt,dt,dt) #todo make sure right\n",
    "            massInverse=8e-6 # todo ?? later change\n",
    "            @inbounds translate=N_linMom[i]*Vector3((dt*massInverse),(dt*massInverse),(dt*massInverse)) # ??massInverse\n",
    "            \n",
    "            @inbounds N_currPosition[i]=N_currPosition[i]+translate\n",
    "            @inbounds N_displacement[i]=N_displacement[i]+translate\n",
    "            \n",
    "            # Rotation\n",
    "            @inbounds curMoment = moment(N_intMoment[i],N_orient[i],N_moment[i]) \n",
    "            @inbounds N_intMoment[i]=Vector3(0,0,0) # do i really need it?\n",
    "            \n",
    "            @inbounds N_angMom[i]=N_angMom[i]+curMoment*Vector3(dt,dt,dt)\n",
    "            momentInertiaInverse=1.0 # todo ?? later change\n",
    "            @inbounds temp=FromRotationVector(N_angMom[i]*Vector3((dt*momentInertiaInverse),(dt*momentInertiaInverse),(dt*momentInertiaInverse)))\n",
    "            \n",
    "            @inbounds x= N_orient[i].x*temp.x\n",
    "            @inbounds y= N_orient[i].y*temp.y\n",
    "            @inbounds z= N_orient[i].z*temp.z\n",
    "            @inbounds w= N_orient[i].w*temp.w\n",
    "            x=convert(Float32,x)\n",
    "            y=convert(Float32,y)\n",
    "            z=convert(Float32,z)\n",
    "            w=convert(Float32,w)\n",
    "            \n",
    "            @inbounds N_orient[i]=Quaternion(x,y,z,w)\n",
    "        end\n",
    "    end\n",
    "    return\n",
    "end\n",
    "\n",
    "\n",
    "function run_updateNodes!(dt,currentTimeStep,N_position, N_restrained,N_displacement, N_angle,N_currPosition,N_linMom,N_angMom,\n",
    "        N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\n",
    "    \n",
    "    N=length(N_intForce)\n",
    "    numblocks = ceil(Int, N/256)\n",
    "    CuArrays.@sync begin\n",
    "        @cuda threads=256 blocks=numblocks updateNodes!(dt,currentTimeStep,N_position, N_restrained,N_displacement, N_angle,N_currPosition,N_linMom,N_angMom,\n",
    "                                                N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\n",
    "    end\n",
    "end\n",
    "\n",
    "# dt=0.1\n",
    "# currentTimeStep=0.0\n",
    "# @timed run_updateNodes!(dt,currentTimeStep,\n",
    "#         N_positionGPU, \n",
    "#         N_restrainedGPU,\n",
    "#         N_displacementGPU,\n",
    "#         N_angleGPU,\n",
    "#         N_currPositionGPU,\n",
    "#         N_linMomGPU,\n",
    "#         N_angMomGPU,\n",
    "#         N_intForceGPU,\n",
    "#         N_intMomentGPU,\n",
    "#         N_forceGPU,\n",
    "#         N_momentGPU,\n",
    "#         N_orientGPU,\n",
    "#         N_edgeIDGPU, \n",
    "#         N_edgeFirstGPU, \n",
    "#         E_intForce1GPU,\n",
    "#         E_intMoment1GPU,\n",
    "#         E_intForce2GPU,\n",
    "#         E_intMoment2GPU)\n",
    "\n",
    "# @timed bench_updateIntForces!(N_intForceGPU, N_edgeIDGPU, N_edgeFirstGPU, E_intForce1GPU, E_intForce2GPU)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "orientLink! (generic function with 1 method)"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# function orientLink!(currentRestLength::Vector3,pVNeg::Vector3,pVPos::Vector3,oVNeg::Quaternion,oVPos::Quaternion,axis::Vector3)  # updates pos2, angle1, angle2, and smallAngle //Quat3D<double> /*double restLength*/\n",
    "function orientLink!(currentRestLength,pVNeg,pVPos,oVNeg,oVPos,axis)  # updates pos2, angle1, angle2, and smallAngle //Quat3D<double> /*double restLength*/\n",
    "        \n",
    "    pos2 = toAxisXVector3(pVPos-pVNeg,axis) # digit truncation happens here...\n",
    "    angle1 = toAxisXQuat(oVNeg,axis)\n",
    "    angle2 = toAxisXQuat(oVPos,axis)\n",
    "    \n",
    "    \n",
    "    # println(angle1[2],\" \",angle1[3],\" \",angle1[4],\" \",angle1[1])\n",
    "    \n",
    "    totalRot = conjugate(angle1) #keep track of the total rotation of this bond (after toAxisX()) # Quat3D<double>\n",
    "    pos2 = RotateVec3D(totalRot,pos2)\n",
    "    \n",
    "    \n",
    "    angle2 = Quaternion(angle2.x*totalRot.x,angle2.y*totalRot.y,angle2.z*totalRot.z,angle2.w*totalRot.w)\n",
    "    angle1 = Quaternion(0.0,0.0,0.0,1.0)#new THREE.Quaternion() #zero for now...\n",
    "\n",
    "    smallAngle=true #todo later remove\n",
    "    \n",
    "    \n",
    "    if (smallAngle)\t #Align so Angle1 is all zeros\n",
    "        #pos2[1] =pos2[1]- currentRestLength #only valid for small angles\n",
    "        pos2=Vector3(pos2.x-currentRestLength,pos2.y,pos2.z)\n",
    "    else  #Large angle. Align so that Pos2.y, Pos2.z are zero.\n",
    "        # FromAngleToPosX(angle1,pos2) #get the angle to align Pos2 with the X axis\n",
    "        # totalRot = angle1.clone().multiply(totalRot)  #update our total rotation to reflect this\n",
    "        # angle2 = angle1.clone().multiply(  angle2) #rotate angle2\n",
    "        # pos2 = new THREE.Vector3(pos2.length() - currentRestLength, 0, 0);\n",
    "    end\n",
    "    \n",
    "    angle1v = ToRotationVector(angle1)\n",
    "    angle2v = ToRotationVector(angle2)\n",
    "#     pos2,angle1v,angle2v,angle1,angle2,\n",
    "    return pos2,angle1v,angle2v,angle1,angle2,totalRot\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "toAxisXVector3 (generic function with 1 method)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function toAxisXVector3(pV::Vector3,axis::Vector3) #TODO CHANGE\n",
    "\n",
    "    xaxis=Vector3(1.0,0.0,0.0)\n",
    "\n",
    "    vector=normalizeVector3(axis)\n",
    "    q=setFromUnitVectors(vector,xaxis)\n",
    "    \n",
    "    # d=17\n",
    "    # qw=round(q[1], digits=d)\n",
    "    # qx=round(q[2], digits=d)\n",
    "    # qy=round(q[3], digits=d)\n",
    "    # qz=round(q[4], digits=d)\n",
    "    \n",
    " \n",
    "    rot=setFromRotationMatrix(quatToMatrix( q  ))\n",
    "    \n",
    "    return applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "toAxisOriginalVector3 (generic function with 1 method)"
      ]
     },
     "execution_count": 71,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function toAxisOriginalVector3(pV::Vector3,axis::Vector3)\n",
    "    \n",
    "    xaxis=Vector3(1.0,0.0,0.0)\n",
    "\n",
    "    vector=normalizeVector3(axis)\n",
    "\n",
    "    q=setFromUnitVectors(xaxis,vector)\n",
    "    \n",
    "\n",
    "    rot=setFromRotationMatrix(quatToMatrix( q  ))\n",
    "\n",
    "    return applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "toAxisXQuat (generic function with 1 method)"
      ]
     },
     "execution_count": 201,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function  toAxisXQuat(pQ::Quaternion,axis::Vector3)\n",
    "    \n",
    "    xaxis=Vector3(1.0,0.0,0.0)\n",
    "\n",
    "    vector=normalizeVector3(axis)\n",
    "\n",
    "\n",
    "    q=setFromUnitVectors(vector,xaxis)\n",
    "        \n",
    "    #d=17\n",
    "    #qw=round(q[1], digits=d)\n",
    "    #qx=round(q[2], digits=d)\n",
    "    #qy=round(q[3], digits=d)\n",
    "    #qz=round(q[4], digits=d)\n",
    "    #\n",
    "\n",
    "    rot=setFromRotationMatrix(quatToMatrix( q  ))\n",
    "    \n",
    "    pV=Vector3(pQ.x,pQ.y,pQ.z)\n",
    "    v=applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
    "    \n",
    "    return Quaternion(v.x,v.y,v.z,1.0)\n",
    "    \n",
    "    # return [1.0 v[1] v[2] v[3]]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "toAxisOriginalQuat (generic function with 2 methods)"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function toAxisOriginalQuat(pQ::Vector3,axis::Vector3)\n",
    "    xaxis=Vector3(1.0,0.0,0.0)\n",
    "\n",
    "    vector=normalizeVector3(axis)\n",
    "    \n",
    "    q=setFromUnitVectors(xaxis,vector)\n",
    "    \n",
    "\n",
    "    rot=setFromRotationMatrix(quatToMatrix( q  ))\n",
    "    \n",
    "    pV=Vector3(pQ.x,pQ.y,pQ.z)\n",
    "    v=applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
    "    \n",
    "    return Quaternion(v.x,v.y,v.z,1.0)\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "setFromUnitVectors (generic function with 1 method)"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function setFromUnitVectors(vFrom::Vector3, vTo::Vector3)\n",
    "    # assumes direction vectors vFrom and vTo are normalized\n",
    "    EPS = 0.000001;\n",
    "    r= dotVector3(vFrom,vTo)+1.0\n",
    "    # r =  dot(vFrom,vTo)+1\n",
    "\n",
    "    if r < EPS\n",
    "        r = 0;\n",
    "        if abs( vFrom.x ) > abs( vFrom.z ) \n",
    "            qx = - vFrom.y\n",
    "            qy = vFrom.x\n",
    "            qz = 0\n",
    "            qw = r\n",
    "        else \n",
    "            qx = 0\n",
    "            qy = - vFrom.z\n",
    "            qz = vFrom.y\n",
    "            qw = r\n",
    "        end\n",
    "   else \n",
    "        # crossVectors( vFrom, vTo ); // inlined to avoid cyclic dependency on Vector3\n",
    "        qx = vFrom.y * vTo.z - vFrom.z * vTo.y\n",
    "        qy = vFrom.z * vTo.x - vFrom.x * vTo.z\n",
    "        qz = vFrom.x * vTo.y - vFrom.y * vTo.x\n",
    "        qw = r\n",
    "\n",
    "    end\n",
    "    qx= (qx==-0) ? 0 : qx\n",
    "    qy= (qy==-0) ? 0 : qy\n",
    "    qz= (qz==-0) ? 0 : qz\n",
    "    qw= (qw==-0) ? 0 : qw\n",
    "        \n",
    "    \n",
    "    mx=qx*qx\n",
    "    my=qy*qy\n",
    "    mz=qz*qz\n",
    "    mw=qw*qw\n",
    "    mm=mx+my\n",
    "    mm=mm+mz\n",
    "    mm=mm+mw\n",
    "    mm=convert(Float32,mm)#??????????????????? todo check later\n",
    "    \n",
    "    l=CUDAnative.sqrt(mm)\n",
    "    #@cuprintln(CUDAnative.sqrt(4.0))\n",
    "    \n",
    "    #l = sqrt((qx * qx) + (qy * qy) + (qz * qz)+ (qw * qw))\n",
    "    if l === 0 \n",
    "        qx = 0.0\n",
    "        qy = 0.0\n",
    "        qz = 0.0\n",
    "        qw = 1.0\n",
    "    else \n",
    "        l = 1.0 / l\n",
    "        qx = qx * l\n",
    "        qy = qy * l\n",
    "        qz = qz * l\n",
    "        qw = qw * l\n",
    "    end\n",
    "\n",
    "    # return qx,qy,qz,qw\n",
    "    return Quaternion(qx,qy,qz,qw)\n",
    "    \n",
    "    # return normalizeQ(Quat(qw,qx,qy,qz))\n",
    "    # return Quat(nn[1], nn[2], nn[3], nn[4])\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "quatToMatrix (generic function with 1 method)"
      ]
     },
     "execution_count": 204,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function quatToMatrix( quaternion::Quaternion)\n",
    "\n",
    "    #te = RotationMatrix()\n",
    "    \n",
    "    x = quaternion.x\n",
    "    y = quaternion.y\n",
    "    z = quaternion.z\n",
    "    w = quaternion.w\n",
    "    \n",
    "    x2 = x + x\n",
    "    y2 = y + y\n",
    "    z2 = z + z\n",
    "    xx = x * x2\n",
    "    xy = x * y2\n",
    "    xz = x * z2\n",
    "    yy = y * y2\n",
    "    yz = y * z2\n",
    "    zz = z * z2\n",
    "    wx = w * x2\n",
    "    wy = w * y2\n",
    "    wz = w * z2\n",
    "\n",
    "    sx = 1.0\n",
    "    sy = 1.0\n",
    "    sz = 1.0\n",
    "\n",
    "    te1 = ( 1 - ( yy + zz ) ) * sx\n",
    "    te2 = ( xy + wz ) * sx\n",
    "    te3 = ( xz - wy ) * sx\n",
    "    te4 = 0;\n",
    "\n",
    "    te5 = ( xy - wz ) * sy\n",
    "    te6 = ( 1 - ( xx + zz ) ) * sy\n",
    "    te7 = ( yz + wx ) * sy\n",
    "    te8 = 0;\n",
    "\n",
    "    te9 = ( xz + wy ) * sz\n",
    "    te10 = ( yz - wx ) * sz\n",
    "    te11 = ( 1 - ( xx + yy ) ) * sz\n",
    "    te12 = 0\n",
    "\n",
    "    te13 = 0 #position.x;\n",
    "    te14 = 0 #position.y;\n",
    "    te15 = 0 #position.z;\n",
    "    te16 = 1\n",
    "    \n",
    "    te= RotationMatrix(te1,te2,te3,te4,te5,te6,te7,te8,te9,te10,te11,te12,te13,te14,te15,te16)\n",
    "\n",
    "    return te\n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "setFromRotationMatrix (generic function with 1 method)"
      ]
     },
     "execution_count": 205,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function  setFromRotationMatrix(m::RotationMatrix)\n",
    "    #te = m\n",
    "    #m11 = (te[ 1 ]== -0.0) ? 0.0 : te[ 1 ]\n",
    "    #m12 = (te[ 5 ]== -0.0) ? 0.0 : te[ 5 ]\n",
    "    #m13 = (te[ 9 ]== -0.0) ? 0.0 : te[ 9 ]\n",
    "    #m21 = (te[ 2 ]== -0.0) ? 0.0 : te[ 2 ]\n",
    "    #m22 = (te[ 6 ]== -0.0) ? 0.0 : te[ 6 ]\n",
    "    #m23 = (te[ 10]== -0.0) ? 0.0 : te[ 10]\n",
    "    #m31 = (te[ 3 ]== -0.0) ? 0.0 : te[ 3 ]\n",
    "    #m32 = (te[ 7 ]== -0.0) ? 0.0 : te[ 7 ]\n",
    "    #m33 = (te[ 11]== -0.0) ? 0.0 : te[ 11]\n",
    "\n",
    "    m11 = m.te1 \n",
    "    m12 = m.te5 \n",
    "    m13 = m.te9 \n",
    "    m21 = m.te2 \n",
    "    m22 = m.te6 \n",
    "    m23 = m.te10\n",
    "    m31 = m.te3 \n",
    "    m32 = m.te7 \n",
    "    m33 = m.te11\n",
    "\n",
    "\n",
    "\n",
    "    y = asin( clamp( m13, - 1, 1 ) ) ##check if has to be changed to cuda\n",
    "\n",
    "    if ( abs( m13 ) < 0.9999999 ) \n",
    "        \n",
    "        x = atan( - m23, m33 )\n",
    "        z = atan( - m12, m11 )#-m12, m11\n",
    "        # if(m23==0.0)\n",
    "        #     x = atan( m23, m33 )\n",
    "        # end\n",
    "        # if(m12==0.0)\n",
    "        #     z = atan( m12, m11 )\n",
    "        # end\n",
    "\n",
    "    else\n",
    "\n",
    "        x = atan( m32, m22 )\n",
    "        z = 0;\n",
    "\n",
    "    end\n",
    "    \n",
    "    return Vector3(x,y,z)\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "setQuaternionFromEuler (generic function with 1 method)"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function setQuaternionFromEuler(euler::Vector3)\n",
    "    x=euler.x\n",
    "    y=euler.y\n",
    "    z=euler.z\n",
    "    \n",
    "    \n",
    "    c1 = CUDAnative.cos( x / 2.0 )\n",
    "    c2 = CUDAnative.cos( y / 2.0 )\n",
    "    c3 = CUDAnative.cos( z / 2.0 )\n",
    "\n",
    "    s1 = CUDAnative.sin( x / 2.0 )\n",
    "    s2 = CUDAnative.sin( y / 2.0 )\n",
    "    s3 = CUDAnative.sin( z / 2.0 )\n",
    "    \n",
    "   \n",
    "    x = s1 * c2 * c3 + c1 * s2 * s3\n",
    "    y = c1 * s2 * c3 - s1 * c2 * s3\n",
    "    z = c1 * c2 * s3 + s1 * s2 * c3\n",
    "    w = c1 * c2 * c3 - s1 * s2 * s3\n",
    "        \n",
    "    return Quaternion(x,y,z,w)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "applyQuaternion1 (generic function with 1 method)"
      ]
     },
     "execution_count": 207,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function applyQuaternion1(e::Vector3,q2::Quaternion)\n",
    "    x = e.x\n",
    "    y = e.y\n",
    "    z = e.z\n",
    "\n",
    "    qx = q2.x\n",
    "    qy = q2.y\n",
    "    qz = q2.z\n",
    "    qw = q2.w\n",
    "\n",
    "    # calculate quat * vector\n",
    "\n",
    "    ix = qw * x + qy * z - qz * y\n",
    "    iy = qw * y + qz * x - qx * z\n",
    "    iz = qw * z + qx * y - qy * x\n",
    "    iw = - qx * x - qy * y - qz * z\n",
    "\n",
    "    # calculate result * inverse quat\n",
    "\n",
    "    xx = ix * qw + iw * - qx + iy * - qz - iz * - qy\n",
    "    yy = iy * qw + iw * - qy + iz * - qx - ix * - qz\n",
    "    zz = iz * qw + iw * - qz + ix * - qy - iy * - qx\n",
    "\n",
    "    return Vector3(xx,yy,zz)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "conjugate (generic function with 1 method)"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function conjugate(q::Quaternion)\n",
    "    x= (-q.x==-0) ? 0.0 : -q.x\n",
    "    y= (-q.y==-0) ? 0.0 : -q.y\n",
    "    z= (-q.z==-0) ? 0.0 : -q.z\n",
    "    w=q.w\n",
    "    x=convert(Float32,x)\n",
    "    y=convert(Float32,y)\n",
    "    z=convert(Float32,z)\n",
    "    w=convert(Float32,w)\n",
    "    return Quaternion(x,y,z,w)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RotateVec3D (generic function with 1 method)"
      ]
     },
     "execution_count": 209,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function RotateVec3D(a::Quaternion, f::Vector3)   \n",
    "    fx= (f.x==-0) ? 0 : f.x\n",
    "    fy= (f.y==-0) ? 0 : f.y\n",
    "    fz= (f.z==-0) ? 0 : f.z\n",
    "    # fx= f.x\n",
    "    # fy= f.y\n",
    "    # fz= f.z\n",
    "    tw = fx*a.x + fy*a.y + fz*a.z\n",
    "    tx = fx*a.w - fy*a.z + fz*a.y\n",
    "    ty = fx*a.z + fy*a.w - fz*a.x\n",
    "    tz = -fx*a.y + fy*a.x + fz*a.w\n",
    "\n",
    "    return Vector3((a.w*tx+a.x*tw+a.y*tz-a.z*ty),(a.w*ty-a.x*tz+a.y*tw+a.z*tx),(a.w*tz+a.x*ty-a.y*tx+a.z*tw))\n",
    "end\n",
    "#!< Returns a vector representing the specified vector \"f\" rotated by this quaternion. @param[in] f The vector to transform.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "RotateVec3DInv (generic function with 1 method)"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function RotateVec3DInv(a::Quaternion, f::Vector3)  \n",
    "    fx=f.x\n",
    "    fy=f.y\n",
    "    fz=f.z\n",
    "    tw = a.x*fx + a.y*fy + a.z*fz\n",
    "    tx = a.w*fx - a.y*fz + a.z*fy\n",
    "    ty = a.w*fy + a.x*fz - a.z*fx\n",
    "    tz = a.w*fz - a.x*fy + a.y*fx\n",
    "    return Vector3((tw*a.x + tx*a.w + ty*a.z - tz*a.y),(tw*a.y - tx*a.z + ty*a.w + tz*a.x),(tw*a.z + tx*a.y - ty*a.x + tz*a.w))\n",
    "end\n",
    "#!< Returns a vector representing the specified vector \"f\" rotated by the inverse of this quaternion. This is the opposite of RotateVec3D. @param[in] f The vector to transform.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ToRotationVector (generic function with 1 method)"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ToRotationVector(a::Quaternion)  \n",
    "    if (a.w >= 1.0 || a.w <= -1.0) \n",
    "        return Vector3(0.0,0.0,0.0)\n",
    "    end\n",
    "    squareLength = 1.0-a.w*a.w; # because x*x + y*y + z*z + w*w = 1.0, but more susceptible to w noise (when \n",
    "    SLTHRESH_ACOS2SQRT= 2.4e-3; # SquareLength threshhold for when we can use square root optimization for acos. From SquareLength = 1-w*w. (calculate according to 1.0-W_THRESH_ACOS2SQRT*W_THRESH_ACOS2SQRT\n",
    "\n",
    "    if (squareLength < SLTHRESH_ACOS2SQRT) # ???????\n",
    "        x=a.x*(2.0*CUDAnative.sqrt((2-2*a.w)/squareLength))\n",
    "        y=a.y*(2.0*CUDAnative.sqrt((2-2*a.w)/squareLength))\n",
    "        z=a.z*(2.0*CUDAnative.sqrt((2-2*a.w)/squareLength))\n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    " \n",
    "        return Vector3(x,y,z) ; # acos(w) = sqrt(2*(1-x)) for w close to 1. for w=0.001, error is 1.317e-6\n",
    "    else \n",
    "        x=a.x*(2.0*CUDAnative.acos(a.w)/CUDAnative.sqrt(squareLength))\n",
    "        y=a.y*(2.0*CUDAnative.acos(a.w)/CUDAnative.sqrt(squareLength))\n",
    "        z=a.z*(2.0*CUDAnative.acos(a.w)/CUDAnative.sqrt(squareLength))\n",
    "        x=convert(Float32,x)\n",
    "        y=convert(Float32,y)\n",
    "        z=convert(Float32,z)\n",
    "\n",
    "        return Vector3(x,y,z)\n",
    "    end                                    \n",
    "end \n",
    "# !< Returns a rotation vector representing this quaternion rotation. Adapted from http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToAngle/\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "updateStrain (generic function with 1 method)"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function updateStrain( axialStrain::Float32,E::Float32) # ?from where strain\n",
    "    strain = axialStrain # redundant?\n",
    "    currentTransverseStrainSum=0.0 # ??? todo\n",
    "    linear=true\n",
    "    maxStrain=100000000000000000;# ?? todo later change\n",
    "    if linear\n",
    "        if axialStrain > maxStrain\n",
    "            maxStrain = axialStrain # remember this maximum for easy reference\n",
    "        end\n",
    "        return stress(axialStrain,E)\n",
    "    else \n",
    "        if (axialStrain > maxStrain) # if new territory on the stress/strain curve\n",
    "            maxStrain = axialStrain # remember this maximum for easy reference\n",
    "            returnStress = stress(axialStrain,E) # ??currentTransverseStrainSum\n",
    "            if (nu != 0.0) \n",
    "                strainOffset = maxStrain-stress(axialStrain,E)/(_eHat*(1-nu)) # precalculate strain offset for when we back off\n",
    "            else \n",
    "                strainOffset = maxStrain-returnStress/E # precalculate strain offset for when we back off\n",
    "            end\n",
    "        else  # backed off a non-linear material, therefore in linear region.\n",
    "            relativeStrain = axialStrain-strainOffset #  treat the material as linear with a strain offset according to the maximum plastic deformation\n",
    "            if (nu != 0.0) \n",
    "                returnStress = stress(relativeStrain,E)\n",
    "            else \n",
    "                returnStress = E*relativeStrain\n",
    "            end\n",
    "        end\n",
    "        return returnStress\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "stress (generic function with 1 method)"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function stress( strain::Float32 , E::Float32 ) #end,transverseStrainSum, forceLinear){\n",
    "    #  reference: http://www.colorado.edu/engineering/CAS/courses.d/Structures.d/IAST.Lect05.d/IAST.Lect05.pdf page 10\n",
    "    #  if (isFailed(strain)) return 0.0f; //if a failure point is set and exceeded, we've broken!\n",
    "    #   var E =setup.edges[0].stiffness; //todo change later to material ??\n",
    "    #   var E=1000000;//todo change later to material ??\n",
    "    #   var scaleFactor=1;\n",
    "    #     @cuprintln(E*strain)\n",
    "    return E*strain;\n",
    "\n",
    "    #  #   if (strain <= strainData[1] || linear || forceLinear){ //for compression/first segment and linear materials (forced or otherwise), simple calculation\n",
    "\n",
    "        #   if (nu==0.0) return E*strain;\n",
    "        #   else return _eHat*((1-nu)*strain + nu*transverseStrainSum); \n",
    "        #  else return eHat()*((1-nu)*strain + nu*transverseStrainSum); \n",
    "    #  #  }\n",
    "\n",
    "      #//the non-linear feature with non-zero poissons ratio is currently experimental\n",
    "      #int DataCount = modelDataPoints();\n",
    "      #for (int i=2; i<DataCount; i++){ //go through each segment in the material model (skipping the first segment because it has already been handled.\n",
    "      #  if (strain <= strainData[i] || i==DataCount-1){ //if in the segment ending with this point (or if this is the last point extrapolate out) \n",
    "      #      float Perc = (strain-strainData[i-1])/(strainData[i]-strainData[i-1]);\n",
    "      #      float basicStress = stressData[i-1] + Perc*(stressData[i]-stressData[i-1]);\n",
    "      #      if (nu==0.0f) return basicStress;\n",
    "      #      else { //accounting for volumetric effects\n",
    "      #          float modulus = (stressData[i]-stressData[i-1])/(strainData[i]-strainData[i-1]);\n",
    "      #          float modulusHat = modulus/((1-2*nu)*(1+nu));\n",
    "      #          float effectiveStrain = basicStress/modulus; //this is the strain at which a simple linear stress strain line would hit this point at the definied modulus\n",
    "      #          float effectiveTransverseStrainSum = transverseStrainSum*(effectiveStrain/strain);\n",
    "      #          return modulusHat*((1-nu)*effectiveStrain + nu*effectiveTransverseStrainSum);\n",
    "      #      }\n",
    "      #  }\n",
    "      #}\n",
    "\n",
    "    #  assert(false); //should never reach this point\n",
    "    #  return 0.0f;\n",
    "end "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "force (generic function with 2 methods)"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function force(N_intForce::Vector3,N_orient::Quaternion,N_force::Vector3,static::Bool,currentTimeStep::Float32) \n",
    "    # forces from internal bonds\n",
    "    totalForce=Vector3(0,0,0)\n",
    "    # new THREE.Vector3(node.force.x,node.force.y,node.force.z);\n",
    "    #  todo \n",
    "\n",
    "\n",
    "    totalForce=totalForce+N_intForce\n",
    "\n",
    "    #  for (int i=0; i<6; i++){ \n",
    "    #  \tif (links[i]) totalForce += links[i]->force(isNegative((linkDirection)i)); # total force in LCS\n",
    "    #  }\n",
    "    totalForce = RotateVec3D(N_orient,totalForce); # from local to global coordinates\n",
    "\n",
    "\n",
    "    # assert(!(totalForce.x != totalForce.x) || !(totalForce.y != totalForce.y) || !(totalForce.z != totalForce.z)); //assert non QNAN\n",
    "\n",
    "    # other forces\n",
    "    if(static)\n",
    "        totalForce=totalForce+N_force\n",
    "    #  }else if(currentTimeStep<50){\n",
    "    #  \ttotalForce.add(new THREE.Vector3(node.force.x,node.force.y,node.force.z));\n",
    "    else\n",
    "        #  var ex=0.1;\n",
    "        #  if(node.force.y!=0){\n",
    "        #  \tvar f=400*Math.sin(currentTimeStep*ex);\n",
    "        #  \ttotalForce.add(new THREE.Vector3(0,f,0));\n",
    "\n",
    "        #  }\n",
    "        #x=N_position[node][3]\n",
    "        #t=currentTimeStep\n",
    "        #wave=getForce(x,t)\n",
    "        #totalForce=totalForce+[0 wave 0]\n",
    "    end\n",
    "\n",
    "\n",
    "    #  if (externalExists()) totalForce += external()->force(); //external forces\n",
    "    #  totalForce -= velocity()*mat->globalDampingTranslateC(); //global damping f-cv\n",
    "    #  totalForce.z += mat->gravityForce(); //gravity, according to f=mg\n",
    "\n",
    "    #  if (isCollisionsEnabled()){\n",
    "    #  \tfor (std::vector<CVX_Collision*>::iterator it=colWatch->begin(); it!=colWatch->end(); it++){\n",
    "    #  \t\ttotalForce -= (*it)->contactForce(this);\n",
    "    #  \t}\n",
    "    #  }\n",
    "    # todo make internal forces 0 again\n",
    "    # N_intForce[node]=[0 0 0] # do i really need it?\n",
    "\n",
    "    #  node.force.x=0;\n",
    "    #  node.force.y=0;\n",
    "    #  node.force.z=0;\n",
    "\n",
    "\n",
    "    return totalForce\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "moment (generic function with 1 method)"
      ]
     },
     "execution_count": 215,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function moment(N_intMoment::Vector3,N_orient::Quaternion,N_moment::Vector3) \n",
    "    #moments from internal bonds\n",
    "    totalMoment=Vector3(0,0,0)\n",
    "    # for (int i=0; i<6; i++){ \n",
    "    # \tif (links[i]) totalMoment += links[i]->moment(isNegative((linkDirection)i)); //total force in LCS\n",
    "    # }\n",
    "\n",
    "    totalMoment=totalMoment+N_intMoment\n",
    "\n",
    "    totalMoment = RotateVec3D(N_orient,totalMoment);\n",
    "\n",
    "    totalMoment=totalMoment+N_moment\n",
    "\n",
    "    #other moments\n",
    "    # if (externalExists()) totalMoment += external()->moment(); //external moments\n",
    "    # totalMoment -= angularVelocity()*mat->globalDampingRotateC(); //global damping\n",
    "\n",
    "    return totalMoment\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FromRotationVector (generic function with 1 method)"
      ]
     },
     "execution_count": 216,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function FromRotationVector(VecIn::Vector3)\n",
    "    theta=VecIn*Vector3(0.5,0.5,0.5)\n",
    "    ntheta=CUDAnative.sqrt((theta.x * theta.x) + (theta.y * theta.y) + (theta.z * theta.z))\n",
    "    thetaMag2=ntheta*ntheta\n",
    "    \n",
    "    DBL_EPSILONx24 =5.328e-15\n",
    "    if thetaMag2*thetaMag2 < DBL_EPSILONx24\n",
    "        qw=1.0 - 0.5*thetaMag2\n",
    "\t\ts=1.0 - thetaMag2 / 6.0\n",
    "    else\n",
    "        thetaMag = CUDAnative.sqrt(thetaMag2)\n",
    "\t\tqw=CUDAnative.cos(thetaMag)\n",
    "\t\ts=CUDAnative.sin(thetaMag) / thetaMag\n",
    "    end\n",
    "    qx=theta.x*s\n",
    "    qy=theta.y*s\n",
    "    qz=theta.z*s\n",
    "    \n",
    "    qx=convert(Float32,qx)\n",
    "    qy=convert(Float32,qy)\n",
    "    qz=convert(Float32,qz)\n",
    "    qw=convert(Float32,qw)\n",
    "    \n",
    "    return Quaternion(qx,qy,qz,qw)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip6700\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6700)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6701\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6700)\" d=\"\n",
       "M202.373 1425.62 L2352.76 1425.62 L2352.76 121.675 L202.373 121.675  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6702\">\n",
       "    <rect x=\"202\" y=\"121\" width=\"2151\" height=\"1305\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  242.741,1425.62 242.741,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  755.03,1425.62 755.03,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1267.32,1425.62 1267.32,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1779.61,1425.62 1779.61,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2291.9,1425.62 2291.9,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.373,1226.97 2352.76,1226.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.373,995.432 2352.76,995.432 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.373,763.893 2352.76,763.893 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.373,532.354 2352.76,532.354 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  202.373,300.815 2352.76,300.815 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,1425.62 2352.76,1425.62 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,1425.62 202.373,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  242.741,1425.62 242.741,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  755.03,1425.62 755.03,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1267.32,1425.62 1267.32,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1779.61,1425.62 1779.61,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2291.9,1425.62 2291.9,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,1226.97 228.178,1226.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,995.432 228.178,995.432 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,763.893 228.178,763.893 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,532.354 228.178,532.354 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  202.373,300.815 228.178,300.815 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 242.741, 1479.62)\" x=\"242.741\" y=\"1479.62\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 755.03, 1479.62)\" x=\"755.03\" y=\"1479.62\">250</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1267.32, 1479.62)\" x=\"1267.32\" y=\"1479.62\">500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1779.61, 1479.62)\" x=\"1779.61\" y=\"1479.62\">750</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2291.9, 1479.62)\" x=\"2291.9\" y=\"1479.62\">1000</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 178.373, 1244.47)\" x=\"178.373\" y=\"1244.47\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 178.373, 1012.93)\" x=\"178.373\" y=\"1012.93\">20</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 178.373, 781.393)\" x=\"178.373\" y=\"781.393\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 178.373, 549.854)\" x=\"178.373\" y=\"549.854\">40</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 178.373, 318.315)\" x=\"178.373\" y=\"318.315\">50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1277.56, 73.2)\" x=\"1277.56\" y=\"73.2\">Benchmark 1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1277.56, 1559.48)\" x=\"1277.56\" y=\"1559.48\">Timesteps</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 89.2861, 773.647)\" x=\"89.2861\" y=\"773.647\">time(sec)</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  263.233,1154.27 447.657,1153.23 652.572,1152.17 1062.4,1150.24 2291.9,1143.61 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6702)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  263.233,1388.71 447.657,1280.21 652.572,1159.03 1062.4,906.472 2291.9,158.579 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip6700)\" d=\"\n",
       "M1936.45 386.635 L2280.76 386.635 L2280.76 205.195 L1936.45 205.195  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1936.45,386.635 2280.76,386.635 2280.76,205.195 1936.45,205.195 1936.45,386.635 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1960.45,265.675 2104.45,265.675 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2128.45, 283.175)\" x=\"2128.45\" y=\"283.175\">GPU</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip6700)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1960.45,326.155 2104.45,326.155 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6700)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2128.45, 343.655)\" x=\"2128.45\" y=\"343.655\">CPU</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot1 4*4*4 voxels, 300 nodes 960 edges\n",
    "numStep=[10,100,200,400,1000]\n",
    "firstStepAvg=13.1276\n",
    "performanceGPU=[0.012447,0.05713079,0.102851299,0.1864829,0.4725757]\n",
    "performanceCPU=[3.014437901,7.700779099,12.9343908,23.8421247,56.1430382]\n",
    "plot(numStep,[firstStepAvg.+performanceGPU,performanceCPU],label=[\"GPU\" \"CPU\"],title=\"Benchmark 1\",xlabel=\"Timesteps\",ylabel=\"time(sec)\")\n",
    "# png(\"../../../02_Presentation/GPUBench1.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip6300\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6300)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6301\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip6300)\" d=\"\n",
       "M229.135 1425.62 L2352.76 1425.62 L2352.76 121.675 L229.135 121.675  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip6302\">\n",
       "    <rect x=\"229\" y=\"121\" width=\"2125\" height=\"1305\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  289.238,1425.62 289.238,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  623.14,1425.62 623.14,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  957.043,1425.62 957.043,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1290.95,1425.62 1290.95,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1624.85,1425.62 1624.85,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1958.75,1425.62 1958.75,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2292.65,1425.62 2292.65,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  229.135,1258.01 2352.76,1258.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  229.135,1028.23 2352.76,1028.23 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  229.135,798.46 2352.76,798.46 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  229.135,568.686 2352.76,568.686 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  229.135,338.913 2352.76,338.913 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,1425.62 2352.76,1425.62 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,1425.62 229.135,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  289.238,1425.62 289.238,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  623.14,1425.62 623.14,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  957.043,1425.62 957.043,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1290.95,1425.62 1290.95,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1624.85,1425.62 1624.85,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1958.75,1425.62 1958.75,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2292.65,1425.62 2292.65,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,1258.01 254.619,1258.01 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,1028.23 254.619,1028.23 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,798.46 254.619,798.46 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,568.686 254.619,568.686 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  229.135,338.913 254.619,338.913 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 289.238, 1479.62)\" x=\"289.238\" y=\"1479.62\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 623.14, 1479.62)\" x=\"623.14\" y=\"1479.62\">5</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 957.043, 1479.62)\" x=\"957.043\" y=\"1479.62\">6</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1290.95, 1479.62)\" x=\"1290.95\" y=\"1479.62\">7</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1624.85, 1479.62)\" x=\"1624.85\" y=\"1479.62\">8</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1958.75, 1479.62)\" x=\"1958.75\" y=\"1479.62\">9</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2292.65, 1479.62)\" x=\"2292.65\" y=\"1479.62\">10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 205.135, 1275.51)\" x=\"205.135\" y=\"1275.51\">30</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 205.135, 1045.73)\" x=\"205.135\" y=\"1045.73\">60</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 205.135, 815.96)\" x=\"205.135\" y=\"815.96\">90</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 205.135, 586.186)\" x=\"205.135\" y=\"586.186\">120</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 205.135, 356.413)\" x=\"205.135\" y=\"356.413\">150</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1290.95, 73.2)\" x=\"1290.95\" y=\"73.2\">Benchmark 2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1290.95, 1559.48)\" x=\"1290.95\" y=\"1559.48\">Lattice Size</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 89.2861, 773.647)\" x=\"89.2861\" y=\"773.647\">time(sec)</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  289.238,1386.45 623.14,1386.45 957.043,1386.43 1290.95,1386.44 1624.85,1386 1958.75,1385.95 2292.65,1385.64 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6302)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  289.238,1388.71 623.14,1311.64 957.043,1190.31 1290.95,1021.06 1624.85,816.944 1958.75,501.963 2292.65,158.579 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip6300)\" d=\"\n",
       "M1936.45 386.635 L2280.76 386.635 L2280.76 205.195 L1936.45 205.195  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1936.45,386.635 2280.76,386.635 2280.76,205.195 1936.45,205.195 1936.45,386.635 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1960.45,265.675 2104.45,265.675 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2128.45, 283.175)\" x=\"2128.45\" y=\"283.175\">GPU</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip6300)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1960.45,326.155 2104.45,326.155 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip6300)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2128.45, 343.655)\" x=\"2128.45\" y=\"343.655\">CPU</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot2 200 timesteps\n",
    "latticeSizes=[4,5,6,7,8,9,10]\n",
    "performanceGPU=[0.102851299,0.1021396,0.1044742,0.1043413,0.1611617,0.1674361,0.2076308]\n",
    "performanceCPU=[12.934390801,22.9971574,38.838537,60.9359617,87.5866625,128.7116549,173.5449189]\n",
    "plot(latticeSizes,[firstStepAvg.+performanceGPU,performanceCPU],label=[\"GPU\" \"CPU\"],title=\"Benchmark 2\",xlabel=\"Lattice Size\",ylabel=\"time(sec)\")\n",
    "# png(\"../../../02_Presentation/GPUBench2.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "CuDevice(0): GeForce GTX 1070 Ti"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CuDevice(0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "using WebSockets, WebIO, JSCall, Interact ,JSExpr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JSCall\n",
    "# using Pkg; Pkg.add(\"JSCall\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.webio.node+json": {
       "children": [
        {
         "children": [],
         "instanceArgs": {
          "namespace": "html",
          "tag": "span"
         },
         "nodeType": "DOM",
         "props": {},
         "type": "node"
        }
       ],
       "instanceArgs": {
        "handlers": {},
        "id": "2719943122382095872",
        "imports": {
         "data": [
          {
           "name": null,
           "type": "js",
           "url": "https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js"
          }
         ],
         "type": "async_block"
        },
        "mount_callbacks": [
         "function () {\n    var handler = (    function (mod){\n        WebIO.object_pool = {}\n        WebIO.object_pool[\"14338817269949188462\"] = mod\n        WebIO.object_pool[\"796947069841572863\"] = document\n        WebIO.object_pool[\"13844523268319912121\"] = window\n        WebIO.object_pool[\"18328825248940393001\"] = this\n         // execute all queued statements in onimport\n    }\n);\n    (WebIO.importBlock({\"data\":[{\"name\":null,\"type\":\"js\",\"url\":\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\"}],\"type\":\"async_block\"})).then((imports) => handler.apply(this, imports));\n}\n"
        ],
        "observables": {
         "_jscall_value_comm": {
          "id": "ob_11",
          "sync": false,
          "value": ""
         }
        },
        "systemjs_options": null
       },
       "nodeType": "Scope",
       "props": {},
       "type": "node"
      },
      "text/plain": [
       "JSModule(Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(Task (runnable) @0x000000001156fc70, Task (runnable) @0x000000001156fc70), Base.AlwaysLockedST(1))), JSString[]), JSObject(:THREE, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(Task (runnable) @0x000000001156fc70, Task (runnable) @0x000000001156fc70), Base.AlwaysLockedST(1))), JSString[]), :module, 0xc6fdb656328ec96e), JSObject(:document, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(Task (runnable) @0x000000001156fc70, Task (runnable) @0x000000001156fc70), Base.AlwaysLockedST(1))), JSString[]), :module, 0x0b0f53233ab203ff), JSObject(:window, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(Task (runnable) @0x000000001156fc70, Task (runnable) @0x000000001156fc70), Base.AlwaysLockedST(1))), JSString[]), :module, 0xc021a08fa4d6a0b9), JSObject(:window, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(Task (runnable) @0x000000001156fc70, Task (runnable) @0x000000001156fc70), Base.AlwaysLockedST(1))), JSString[]), :module, 0xfe5d1174214d8229), identity)"
      ]
     },
     "execution_count": 27,
     "metadata": {
      "application/vnd.webio.node+json": {
       "kernelId": "af013df4-dfe2-458b-ae4a-529a432f5d4e"
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "THREE=JSModule(:THREE,\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.webio.node+json": {
       "children": [
        {
         "children": [],
         "instanceArgs": {
          "namespace": "html",
          "tag": "span"
         },
         "nodeType": "DOM",
         "props": {},
         "type": "node"
        }
       ],
       "instanceArgs": {
        "handlers": {},
        "id": "2719943122382095872",
        "imports": {
         "data": [
          {
           "name": null,
           "type": "js",
           "url": "https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js"
          }
         ],
         "type": "async_block"
        },
        "mount_callbacks": [
         "function () {\n    var handler = (    function (mod){\n        WebIO.object_pool = {}\n        WebIO.object_pool[\"14338817269949188462\"] = mod\n        WebIO.object_pool[\"796947069841572863\"] = document\n        WebIO.object_pool[\"13844523268319912121\"] = window\n        WebIO.object_pool[\"18328825248940393001\"] = this\n         // execute all queued statements in onimport\n    }\n);\n    (WebIO.importBlock({\"data\":[{\"name\":null,\"type\":\"js\",\"url\":\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\"}],\"type\":\"async_block\"})).then((imports) => handler.apply(this, imports));\n}\n",
         "function () {\n    var handler = (    function (mod){\n        WebIO.object_pool = {}\n        WebIO.object_pool[\"14338817269949188462\"] = mod\n        WebIO.object_pool[\"796947069841572863\"] = document\n        WebIO.object_pool[\"13844523268319912121\"] = window\n        WebIO.object_pool[\"18328825248940393001\"] = this\n         // execute all queued statements in onimport\n    }\n);\n    (WebIO.importBlock({\"data\":[{\"name\":null,\"type\":\"js\",\"url\":\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\"}],\"type\":\"async_block\"})).then((imports) => handler.apply(this, imports));\n}\n"
        ],
        "observables": {
         "_jscall_value_comm": {
          "id": "ob_11",
          "sync": false,
          "value": ""
         }
        },
        "systemjs_options": null
       },
       "nodeType": "Scope",
       "props": {},
       "type": "node"
      },
      "text/plain": [
       "JSModule(Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[WebIO.IJuliaConnection(IJulia.CommManager.Comm{:webio_comm}(\"86cde06edea04e98a20b02b99b801168\", false, getfield(WebIO, Symbol(\"##93#94\")){WebIO.IJuliaConnection}(WebIO.IJuliaConnection(#= circular reference @-3 =#)), IJulia.CommManager.noop_callback))]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(nothing, nothing), Base.AlwaysLockedST(1))), JSString[JSString(\"function () {\\n    var handler = (    function (mod){\\n        WebIO.object_pool = {}\\n        WebIO.object_pool[\\\"14338817269949188462\\\"] = mod\\n        WebIO.object_pool[\\\"796947069841572863\\\"] = document\\n        WebIO.object_pool[\\\"13844523268319912121\\\"] = window\\n        WebIO.object_pool[\\\"18328825248940393001\\\"] = this\\n         // execute all queued statements in onimport\\n    }\\n);\\n    (WebIO.importBlock({\\\"data\\\":[{\\\"name\\\":null,\\\"type\\\":\\\"js\\\",\\\"url\\\":\\\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\\\"}],\\\"type\\\":\\\"async_block\\\"})).then((imports) => handler.apply(this, imports));\\n}\\n\")]), JSObject(:THREE, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[WebIO.IJuliaConnection(IJulia.CommManager.Comm{:webio_comm}(\"86cde06edea04e98a20b02b99b801168\", false, getfield(WebIO, Symbol(\"##93#94\")){WebIO.IJuliaConnection}(WebIO.IJuliaConnection(#= circular reference @-3 =#)), IJulia.CommManager.noop_callback))]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(nothing, nothing), Base.AlwaysLockedST(1))), JSString[JSString(\"function () {\\n    var handler = (    function (mod){\\n        WebIO.object_pool = {}\\n        WebIO.object_pool[\\\"14338817269949188462\\\"] = mod\\n        WebIO.object_pool[\\\"796947069841572863\\\"] = document\\n        WebIO.object_pool[\\\"13844523268319912121\\\"] = window\\n        WebIO.object_pool[\\\"18328825248940393001\\\"] = this\\n         // execute all queued statements in onimport\\n    }\\n);\\n    (WebIO.importBlock({\\\"data\\\":[{\\\"name\\\":null,\\\"type\\\":\\\"js\\\",\\\"url\\\":\\\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\\\"}],\\\"type\\\":\\\"async_block\\\"})).then((imports) => handler.apply(this, imports));\\n}\\n\")]), :module, 0xc6fdb656328ec96e), JSObject(:document, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[WebIO.IJuliaConnection(IJulia.CommManager.Comm{:webio_comm}(\"86cde06edea04e98a20b02b99b801168\", false, getfield(WebIO, Symbol(\"##93#94\")){WebIO.IJuliaConnection}(WebIO.IJuliaConnection(#= circular reference @-3 =#)), IJulia.CommManager.noop_callback))]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(nothing, nothing), Base.AlwaysLockedST(1))), JSString[JSString(\"function () {\\n    var handler = (    function (mod){\\n        WebIO.object_pool = {}\\n        WebIO.object_pool[\\\"14338817269949188462\\\"] = mod\\n        WebIO.object_pool[\\\"796947069841572863\\\"] = document\\n        WebIO.object_pool[\\\"13844523268319912121\\\"] = window\\n        WebIO.object_pool[\\\"18328825248940393001\\\"] = this\\n         // execute all queued statements in onimport\\n    }\\n);\\n    (WebIO.importBlock({\\\"data\\\":[{\\\"name\\\":null,\\\"type\\\":\\\"js\\\",\\\"url\\\":\\\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\\\"}],\\\"type\\\":\\\"async_block\\\"})).then((imports) => handler.apply(this, imports));\\n}\\n\")]), :module, 0x0b0f53233ab203ff), JSObject(:window, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[WebIO.IJuliaConnection(IJulia.CommManager.Comm{:webio_comm}(\"86cde06edea04e98a20b02b99b801168\", false, getfield(WebIO, Symbol(\"##93#94\")){WebIO.IJuliaConnection}(WebIO.IJuliaConnection(#= circular reference @-3 =#)), IJulia.CommManager.noop_callback))]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(nothing, nothing), Base.AlwaysLockedST(1))), JSString[JSString(\"function () {\\n    var handler = (    function (mod){\\n        WebIO.object_pool = {}\\n        WebIO.object_pool[\\\"14338817269949188462\\\"] = mod\\n        WebIO.object_pool[\\\"796947069841572863\\\"] = document\\n        WebIO.object_pool[\\\"13844523268319912121\\\"] = window\\n        WebIO.object_pool[\\\"18328825248940393001\\\"] = this\\n         // execute all queued statements in onimport\\n    }\\n);\\n    (WebIO.importBlock({\\\"data\\\":[{\\\"name\\\":null,\\\"type\\\":\\\"js\\\",\\\"url\\\":\\\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\\\"}],\\\"type\\\":\\\"async_block\\\"})).then((imports) => handler.apply(this, imports));\\n}\\n\")]), :module, 0xc021a08fa4d6a0b9), JSObject(:window, Scope(Node{WebIO.DOM}(WebIO.DOM(:html, :span), Any[], Dict{Symbol,Any}()), Dict{String,Tuple{Observables.AbstractObservable,Union{Nothing, Bool}}}(\"_jscall_value_comm\" => (Observable{Any} with 0 listeners. Value:\n",
       "JSCall.NoValue(), nothing)), Set(String[]), nothing, Asset[Asset(\"js\", nothing, \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\")], Dict{Any,Any}(), WebIO.ConnectionPool(Channel{Any}(sz_max:32,sz_curr:0), Set(AbstractConnection[WebIO.IJuliaConnection(IJulia.CommManager.Comm{:webio_comm}(\"86cde06edea04e98a20b02b99b801168\", false, getfield(WebIO, Symbol(\"##93#94\")){WebIO.IJuliaConnection}(WebIO.IJuliaConnection(#= circular reference @-3 =#)), IJulia.CommManager.noop_callback))]), Base.GenericCondition{Base.AlwaysLockedST}(Base.InvasiveLinkedList{Task}(nothing, nothing), Base.AlwaysLockedST(1))), JSString[JSString(\"function () {\\n    var handler = (    function (mod){\\n        WebIO.object_pool = {}\\n        WebIO.object_pool[\\\"14338817269949188462\\\"] = mod\\n        WebIO.object_pool[\\\"796947069841572863\\\"] = document\\n        WebIO.object_pool[\\\"13844523268319912121\\\"] = window\\n        WebIO.object_pool[\\\"18328825248940393001\\\"] = this\\n         // execute all queued statements in onimport\\n    }\\n);\\n    (WebIO.importBlock({\\\"data\\\":[{\\\"name\\\":null,\\\"type\\\":\\\"js\\\",\\\"url\\\":\\\"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\\\"}],\\\"type\\\":\\\"async_block\\\"})).then((imports) => handler.apply(this, imports));\\n}\\n\")]), :module, 0xfe5d1174214d8229), identity)"
      ]
     },
     "execution_count": 29,
     "metadata": {
      "application/vnd.webio.node+json": {
       "kernelId": "af013df4-dfe2-458b-ae4a-529a432f5d4e"
      }
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "THREE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "ename": "MethodError",
     "evalue": "MethodError: no method matching iterate(::JSModule)\nClosest candidates are:\n  iterate(!Matched::Core.SimpleVector) at essentials.jl:604\n  iterate(!Matched::Core.SimpleVector, !Matched::Any) at essentials.jl:604\n  iterate(!Matched::ExponentialBackOff) at error.jl:214\n  ...",
     "output_type": "error",
     "traceback": [
      "MethodError: no method matching iterate(::JSModule)\nClosest candidates are:\n  iterate(!Matched::Core.SimpleVector) at essentials.jl:604\n  iterate(!Matched::Core.SimpleVector, !Matched::Any) at essentials.jl:604\n  iterate(!Matched::ExponentialBackOff) at error.jl:214\n  ...",
      "",
      "Stacktrace:",
      " [1] indexed_iterate(::JSModule, ::Int64) at .\\tuple.jl:66",
      " [2] top-level scope at In[16]:1"
     ]
    }
   ],
   "source": [
    "THREE, document, window = JSCall.JSModule(\n",
    "  :THREE,\n",
    "  \"https://cdnjs.cloudflare.com/ajax/libs/three.js/103/three.js\",\n",
    ")\n",
    "scene = THREE.new.Scene()\n",
    "width, height = 650, 200\n",
    "# Create a basic perspective camera\n",
    "camera = THREE.new.PerspectiveCamera(75, width / height, 0.1, 1000)\n",
    "camera.position.z = 4\n",
    "renderer = THREE.new.WebGLRenderer(Dict(:antialias => true))\n",
    "renderer.setSize(width, height)\n",
    "renderer.setClearColor(\"#fcfcfc\")\n",
    "geometry = THREE.new.BoxGeometry(1, 1, 1)\n",
    "material = THREE.new.MeshBasicMaterial(color = \"#433F81\")\n",
    "cube = THREE.new.Mesh(geometry, material);\n",
    "scene.add(cube)\n",
    "container = document.querySelector(\"#container\")\n",
    "container.appendChild(renderer.domElement);\n",
    "slider = @manipulate for r in LinRange(0.0, 2pi, 200)\n",
    "  cube.rotation.x = r\n",
    "\tcube.rotation.y = r\n",
    "\trenderer.render(scene, camera)\n",
    "  nothing\n",
    "end\n",
    "\n",
    "# Output will show WebIO not detected, as long as you don't start the notebook\n",
    "# more information below\n",
    "vbox(\n",
    "  slider, \n",
    "  JSCall.scope(THREE)(dom\"div#container\"(\n",
    "      style = Dict(:height => height, :width => width)\n",
    "  ))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": "e34f50e635f24b1f82de37ae2d08477b",
   "lastKernelId": "91178c54-c963-4beb-9d0e-1540cf441fa4"
  },
  "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
}