Metavoxel GPU Validation Snake.ipynb

    "        \"E_angle2vGPU\" =>E_angle2vGPU,                   \n",
    "        \"E_angle1GPU\" =>E_angle1GPU,                    \n",
    "        \"E_angle2GPU\" =>E_angle2GPU,                    \n",
    "        \"E_currentTransverseStrainSumGPU\" =>E_currentTransverseStrainSumGPU,\n",
    "        \"E_intForce1GPU\" =>E_intForce1GPU,                 \n",
    "        \"E_intMoment1GPU\" =>E_intMoment1GPU,                \n",
    "        \"E_intForce2GPU\" =>E_intForce2GPU,                 \n",
    "        \"E_intMoment2GPU\" =>E_intMoment2GPU,                \n",
    "        \"E_dampGPU\" =>E_dampGPU                      \n",
    "    )\n",
    "\n",
    "    #########################################\n",
    "    \n",
    "\n",
    "    dt=0.0251646\n",
    "    E =2000  # MPa\n",
    "    s=2.38\n",
    "    mass=1  \n",
    "    \n",
    "    \n",
    "    \n",
    "    MaxFreq2=E*s/mass\n",
    "    dt= 1/(6.283185*sqrt(MaxFreq2))\n",
    "#     dt=0.01\n",
    "    println(\"dt: $dt\")\n",
    "    \n",
    "    append!(displacements,[Array(metavoxel[\"N_displacementGPU\"])])\n",
    "    \n",
    "    t=@timed doTimeStep!(metavoxel,dt,0)\n",
    "    append!(displacements,[Array(metavoxel[\"N_displacementGPU\"])])\n",
    "    time=t[2]\n",
    "    println(\"first timestep took $time seconds\")\n",
    "    t=@timed simulateParallel!(metavoxel,numTimeSteps-1,dt,returnEvery)\n",
    "    time=t[2]\n",
    "    \n",
    "    if save\n",
    "        updateDataAndSave!(metavoxel,setup,\"../json/trialJuliaParallelGPU.json\")\n",
    "        updateDataAndSave!(metavoxel,setup,\"../json/trialJuliaParallelGPUDynamic.json\",displacements)\n",
    "    end\n",
    "    println(\"ran latticeSize $latticeSize with $voxCount voxels and $linkCount edges for $numTimeSteps time steps took $time seconds\")\n",
    "    return\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fea (generic function with 1 method)"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fea(setup)\n",
    "    #######################################################\n",
    "    function points(element, properties)\n",
    "        elements = properties[\"elements\"]\n",
    "        nodes = properties[\"nodes\"]\n",
    "        degrees_of_freedom = properties[\"degrees_of_freedom\"]\n",
    "\n",
    "        # find the nodes that the lements connects\n",
    "        fromNode = elements[element][1]\n",
    "        toNode = elements[element][2]\n",
    "\n",
    "        # the coordinates for each node\n",
    "        fromPoint = nodes[fromNode]\n",
    "        toPoint = nodes[toNode]\n",
    "\n",
    "        # find the degrees of freedom for each node\n",
    "        dofs = degrees_of_freedom[fromNode]\n",
    "        dofs=vcat(dofs,degrees_of_freedom[toNode])\n",
    "\n",
    "        return fromPoint, toPoint, dofs\n",
    "    end\n",
    "\n",
    "    function direction_cosine(vec1, vec2)\n",
    "        return dot(vec1,vec2) / (norm(vec1) * norm(vec2))\n",
    "    end\n",
    "\n",
    "    function rotation_matrix(element_vector, x_axis, y_axis,z_axis)\n",
    "        # find the direction cosines\n",
    "        x_proj = direction_cosine(element_vector, x_axis)\n",
    "        y_proj = direction_cosine(element_vector, y_axis)\n",
    "        z_proj = direction_cosine(element_vector, z_axis);\n",
    "        return [[x_proj y_proj z_proj 0 0 0];[0 0 0 x_proj y_proj z_proj]]\n",
    "    end\n",
    "\n",
    "    function rotation_matrix(element_vector, x_axis, y_axis,z_axis)\n",
    "        # find the direction cosines\n",
    "        L=norm(element_vector)\n",
    "        l = (element_vector[1])/L\n",
    "        m = (element_vector[2])/L\n",
    "        n = (element_vector[3])/L\n",
    "        D = ( l^2+ m^2+n^2)^0.5\n",
    "\n",
    "        transMatrix=[[l m  n  0  0  0  0  0  0  0  0  0];[-m/D l/D  0  0  0  0  0  0  0  0  0  0];[ -l*n/D  -m*n/D  D  0  0  0  0  0  0  0  0  0];[ 0  0  0       l       m  n  0  0  0  0  0  0];[ 0  0  0    -m/D     l/D  0  0  0  0  0  0  0];[ 0  0  0  -l*n/D  -m*n/D  D  0  0  0  0  0  0];[ 0  0  0  0  0  0       l       m  n  0  0  0];[ 0  0  0  0  0  0    -m/D     l/D  0  0  0  0];[ 0  0  0  0  0  0  -l*n/D  -m*n/D  D  0  0  0];[ 0  0  0  0  0  0  0  0  0       l       m  n];[ 0  0  0  0  0  0  0  0  0    -m/D     l/D  0];[ 0  0  0  0  0  0  0  0  0  -l*n/D  -m*n/D  D]]\n",
    "\n",
    "        return transMatrix\n",
    "    end\n",
    "    \n",
    "    #######################################################\n",
    "    function get_matrices(setup)\n",
    "\n",
    "        nodes      = setup[\"nodes\"]\n",
    "        edges      = setup[\"edges\"]\n",
    "        ndofs      = length(nodes)*6\n",
    "\n",
    "        x_axis     = [1 0 0]\n",
    "        y_axis     = [0 1 0]\n",
    "        z_axis     = [0 0 1]\n",
    "\n",
    "        M = zeros((ndofs,ndofs))\n",
    "        K = zeros((ndofs,ndofs))\n",
    "        \n",
    "        \n",
    "        for edge in edges\n",
    "            #degrees_of_freedom = properties[\"degrees_of_freedom\"]\n",
    "\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",
    "            # the coordinates for each node\n",
    "            fromPoint = [fromNode[\"position\"][\"x\"]*15.0 fromNode[\"position\"][\"y\"]*15.0 fromNode[\"position\"][\"z\"]*15.0]\n",
    "            toPoint = [toNode[\"position\"][\"x\"]*15.0 toNode[\"position\"][\"y\"]*15.0 toNode[\"position\"][\"z\"]*15.0]\n",
    "            \n",
    "            fromPoint = fromPoint .+ [fromNode[\"displacement\"][\"x\"]*15.0 fromNode[\"displacement\"][\"y\"]*15.0 fromNode[\"displacement\"][\"z\"]*15.0]\n",
    "            toPoint = toPoint .+ [toNode[\"displacement\"][\"x\"]*15.0 toNode[\"displacement\"][\"y\"]*15.0 toNode[\"displacement\"][\"z\"]*15.0]\n",
    "\n",
    "            # find the degrees of freedom for each node\n",
    "            dofs = convert(Array{Int}, fromNode[\"degrees_of_freedom\"])\n",
    "            dofs=vcat(dofs,convert(Array{Int}, toNode[\"degrees_of_freedom\"]))\n",
    "\n",
    "            element_vector=toPoint-fromPoint\n",
    "\n",
    "            # find element mass and stifness matrices\n",
    "            length   = norm(element_vector)\n",
    "            rho      = edge[\"density\"]\n",
    "            area     = edge[\"area\"]\n",
    "            E        = edge[\"stiffness\"]# youngs modulus\n",
    "\n",
    "            A = edge[\"area\"]\n",
    "            G=1.0#todo shear_modulus\n",
    "            ixx = 1.0#todo section ixx\n",
    "            iyy = 1.0#todo section.iyy#\n",
    "            l0=length\n",
    "            j=1.0;#todo check\n",
    "            l02 = l0 * l0\n",
    "            l03 = l0 * l0 * l0\n",
    "            \n",
    "            # find element mass and stifness matrices\n",
    "            length   = norm(element_vector)\n",
    "            rho      = edge[\"density\"]\n",
    "            area     = edge[\"area\"]\n",
    "            E        = edge[\"stiffness\"]# youngs modulus\n",
    "            \n",
    "            \n",
    "    \n",
    "            A = edge[\"area\"]\n",
    "            G=1.0#todo shear_modulus\n",
    "            ixx = 1.0#todo section ixx\n",
    "            iyy = 1.0#todo section.iyy#\n",
    "            j=1.0;#todo check\n",
    "            \n",
    "            \n",
    "            h = 2.38 # mm\n",
    "            b = 2.38  # mm\n",
    "            E =2000  # MPa\n",
    "            rho = 7.85e-9 / 3  # kg/mm^3\n",
    "            G = E * 1 / 3  # MPa\n",
    "            A=h*b\n",
    "            Q = 1 / 3 - 0.2244 / (min(h / b, b / h) + 0.1607)\n",
    "            J = Q * min(h * b^3, b * h^3)\n",
    "            I= b*h^3/12\n",
    "            ixx=I\n",
    "            iyy=I\n",
    "            j=J\n",
    "            \n",
    "            l0=length\n",
    "            l02 = l0 * l0\n",
    "            l03 = l0 * l0 * l0\n",
    "\n",
    "            # Cm = rho * area * length /6.0\n",
    "            # Ck= E * area / length \n",
    "\n",
    "            # m = [[2 1];[1 2]]\n",
    "            # k = [[1 -1];[-1 1]]\n",
    "\n",
    "            k = [[E*A/l0  0  0  0  0  0  -E*A/l0  0  0  0  0  0];[0  12*E*ixx/l03  0  0  0  6*E*ixx/l02  0  -12*E*ixx/l03  0  0  0  6*E*ixx/l02];[0  0  12*E*iyy/l03  0  -6*E*iyy/l02  0  0  0  -12*E*iyy/l03  0  -6*E*iyy/l02  0];[0  0  0  G*j/l0  0  0  0  0  0  -G*j/l0  0  0];[0  0  -6*E*iyy/l02  0  4*E*iyy/l0  0  0  0  6*E*iyy/l02  0  2*E*iyy/l0  0];[0  6*E*ixx/l02  0  0  0  4*E*ixx/l0  0  -6*E*ixx/l02  0  0  0  2*E*ixx/l0];[-E*A/l0  0  0  0  0  0  E*A/l0  0  0  0  0  0];[0  -12*E*ixx/l03  0  0  0  -6*E*ixx/l02  0  12*E*ixx/l03  0  0  0  -6*E*ixx/l02];[0  0  -12*E*iyy/l03  0  6*E*iyy/l02  0  0  0  12*E*iyy/l03  0  6*E*iyy/l02  0];[0  0  0  -G*j/l0  0  0  0  0  0  G*j/l0  0  0];[0  0  -6*E*iyy/l02  0  2*E*iyy/l0  0  0  0  6*E*iyy/l02  0  4*E*iyy/l0  0];[0  6*E*ixx/l02  0  0  0  2*E*ixx/l0  0  -6*E*ixx/l02  0  0  0  4*E*ixx/l0]]\n",
    "\n",
    "            \n",
    "            ################################\n",
    "\n",
    "            # find rotated mass and stifness matrices\n",
    "            tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis)\n",
    "\n",
    "            # m_r=transpose(tau)*m*tau\n",
    "            k_r=transpose(tau)*k*tau\n",
    "\n",
    "            # change from element to global coordinate\n",
    "            index= dofs.+1\n",
    "\n",
    "            B=zeros((12,ndofs))\n",
    "            for i in 1:12\n",
    "                  B[i,index[i]]=1.0\n",
    "            end\n",
    "\n",
    "\n",
    "            # M_rG= transpose(B)*m_r*B\n",
    "            K_rG= transpose(B)*k_r*B\n",
    "\n",
    "            # M += Cm .* M_rG\n",
    "            # K += Ck .* K_rG\n",
    "            K +=  K_rG\n",
    "\n",
    "        end\n",
    "        \n",
    "        \n",
    "        # construct the force vector\n",
    "        F=zeros(ndofs)\n",
    "        remove_indices=[];\n",
    "        for node in nodes\n",
    "            #insert!(F,i, value);\n",
    "            #F=vcat(F,value)\n",
    "            \n",
    "            \n",
    "            i=parse(Int,node[\"id\"][2:end])\n",
    "            f=node[\"force\"]\n",
    "            \n",
    "            # println(f)\n",
    "            F[(i)*6+1]+=f[\"x\"]\n",
    "            F[(i)*6+2]+=f[\"y\"]\n",
    "            F[(i)*6+3]+=f[\"z\"]\n",
    "            F[(i)*6+4]+=0\n",
    "            F[(i)*6+5]+=0\n",
    "            F[(i)*6+6]+=0\n",
    "            Load+=f[\"y\"]\n",
    "            \n",
    "            if f[\"z\"]==-1.0\n",
    "                currPos=Vector3((node[\"position\"][\"x\"]+node[\"displacement\"][\"x\"])*15.0,(node[\"position\"][\"y\"]+node[\"displacement\"][\"y\"])*15.0,(node[\"position\"][\"z\"]+node[\"displacement\"][\"z\"])*15.0)\n",
    "                currPos=Vector3(0.0,(node[\"position\"][\"y\"]+node[\"displacement\"][\"y\"])*15.0,(node[\"position\"][\"z\"]+node[\"displacement\"][\"z\"])*15.0)\n",
    "\n",
    "                ff=normalizeVector3(currPos)\n",
    "#                 F[(i)*6+1]+=f[\"x\"]*-0.1\n",
    "#                 F[(i)*6+2]+=f[\"y\"]*-0.1\n",
    "#                 F[(i)*6+3]+=f[\"z\"]*-0.1\n",
    "                F[(i)*6+3]=-0.1\n",
    "            else\n",
    "                F[(i)*6+3]=0.0\n",
    "            end\n",
    "            \n",
    "            if (F[(i)*6+2]!=0)\n",
    "                append!(topNodesIndices,i+1)\n",
    "            end\n",
    "            \n",
    "            dofs = convert(Array{Int}, node[\"degrees_of_freedom\"]).+1\n",
    "            restrained_dofs=node[\"restrained_degrees_of_freedom\"]\n",
    "            for (index, value) in enumerate(dofs)\n",
    "                if restrained_dofs[index]\n",
    "                    append!( remove_indices, value)\n",
    "                end\n",
    "            end\n",
    "            \n",
    "        end\n",
    "\n",
    "        #println(remove_indices)\n",
    "        #print(K)\n",
    "        #print(F)\n",
    "        \n",
    "\n",
    "        #M = M[setdiff(1:end, remove_indices), :]\n",
    "        K = K[setdiff(1:end, remove_indices), :]\n",
    "\n",
    "        #M = M[:,setdiff(1:end, remove_indices)]\n",
    "        K = K[:,setdiff(1:end, remove_indices)]\n",
    "\n",
    "        F = F[setdiff(1:end, remove_indices)]\n",
    "        \n",
    "        U=zeros(ndofs)\n",
    "        \n",
    "        return M,K,F,U,remove_indices\n",
    "    end\n",
    "\n",
    "    \n",
    "    function updateDisplacement(setup, X)\n",
    "        nodes= setup[\"nodes\"]\n",
    "        i=0\n",
    "        for node in nodes\n",
    "            \n",
    "#             if !node[\"restrained_degrees_of_freedom\"][2]\n",
    "                #i=parse(Int,node[\"id\"][2:end])\n",
    "                node[\"displacement\"][\"x\"]+=X[(i)*6+1]/15\n",
    "                node[\"displacement\"][\"y\"]+=X[(i)*6+2]/15\n",
    "                node[\"displacement\"][\"z\"]+=X[(i)*6+3]/15\n",
    "                node[\"angle\"][\"x\"]+=X[(i)*6+4]\n",
    "                node[\"angle\"][\"y\"]+=X[(i)*6+5]\n",
    "                node[\"angle\"][\"z\"]+=X[(i)*6+6]\n",
    "                append!(displacementFEA,[Vector3(node[\"displacement\"][\"x\"]*15,node[\"displacement\"][\"y\"]*15,node[\"displacement\"][\"z\"]*15)])\n",
    "                i=i+1\n",
    "#             else\n",
    "#                 append!(displacementFEA,[Vector3(0,0,0)])\n",
    "#             end\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    #######################################################\n",
    "\n",
    "    function get_stresses(setup)\n",
    "        nodes      = setup[\"nodes\"]\n",
    "        edges      = setup[\"edges\"]\n",
    "        ndofs      = length(nodes)*6\n",
    "\n",
    "        x_axis     = [1 0 0]\n",
    "        y_axis     = [0 1 0]\n",
    "        z_axis     = [0 0 1]\n",
    "\n",
    "        # find the stresses in each member\n",
    "        stresses=zeros(length(edges))\n",
    "        max11=-10e6\n",
    "        min11=10e6\n",
    "        for edge in edges\n",
    "            #degrees_of_freedom = properties[\"degrees_of_freedom\"]\n",
    "\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",
    "            # the coordinates for each node\n",
    "            fromPoint = [fromNode[\"position\"][\"x\"]*15.0 fromNode[\"position\"][\"y\"]*15.0 fromNode[\"position\"][\"z\"]*15.0]\n",
    "            toPoint = [toNode[\"position\"][\"x\"]*15.0 toNode[\"position\"][\"y\"]*15.0 toNode[\"position\"][\"z\"]*15.0]\n",
    "            \n",
    "            fromPoint = fromPoint .+ [fromNode[\"displacement\"][\"x\"]*15.0 fromNode[\"displacement\"][\"y\"]*15.0 fromNode[\"displacement\"][\"z\"]*15.0]\n",
    "            toPoint = toPoint .+ [toNode[\"displacement\"][\"x\"]*15.0 toNode[\"displacement\"][\"y\"]*15.0 toNode[\"displacement\"][\"z\"]*15.0]\n",
    "\n",
    "\n",
    "            # find the degrees of freedom for each node\n",
    "            dofs = convert(Array{Int}, fromNode[\"degrees_of_freedom\"])\n",
    "            dofs=vcat(dofs,convert(Array{Int}, toNode[\"degrees_of_freedom\"]))\n",
    "\n",
    "            element_vector=toPoint-fromPoint\n",
    "\n",
    "\n",
    "            # find rotated mass and stifness matrices\n",
    "            tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis)\n",
    "\n",
    "            # i1=parse(Int,fromNode[\"id\"][2:end])\n",
    "            # i2=parse(Int,toNode[\"id\"][2:end])\n",
    "\n",
    "            # global_displacements=[X[(i1)*6+1] X[(i1)*6+2] X[(i1)*6+3] X[(i1)*6+4] X[(i1)*6+5] X[(i1)*6+6] X[(i2)*6+1] X[(i2)*6+2] X[(i2)*6+3] X[(i2)*6+4] X[(i2)*6+5] X[(i2)*6+6]] # todo change\n",
    "            global_displacements=[fromNode[\"displacement\"][\"x\"]*15 fromNode[\"displacement\"][\"y\"]*15 fromNode[\"displacement\"][\"z\"]*15 fromNode[\"angle\"][\"x\"] fromNode[\"angle\"][\"y\"] fromNode[\"angle\"][\"z\"] toNode[\"displacement\"][\"x\"]*15 toNode[\"displacement\"][\"y\"]*15 toNode[\"displacement\"][\"z\"]*15 toNode[\"angle\"][\"x\"] toNode[\"angle\"][\"y\"] toNode[\"angle\"][\"z\"]] # todo change\n",
    "\n",
    "            # nodal displacement\n",
    "\n",
    "            q=tau*transpose(global_displacements)\n",
    "            # println(q)\n",
    "            # calculate the strain and stresses\n",
    "            strain =(q[7]-q[1])/norm(element_vector)\n",
    "            E = edge[\"stiffness\"]# youngs modulus\n",
    "            E =2000\n",
    "            stress=E.*strain\n",
    "            edge[\"stress\"]=stress\n",
    "            if stress>max11\n",
    "                max11=stress\n",
    "            end\n",
    "            if stress<min11\n",
    "                min11=stress\n",
    "            end\n",
    "            # println(element)\n",
    "            # println(stress)\n",
    "        end\n",
    "\n",
    "\n",
    "\n",
    "        setup[\"viz\"][\"minStress\"]=min11\n",
    "        setup[\"viz\"][\"maxStress\"]=max11\n",
    "        setup[\"animation\"][\"exaggeration\"]=20.0\n",
    "        return stresses\n",
    "    end\n",
    "    \n",
    "    function initialize(setup)\n",
    "        nodes      = setup[\"nodes\"]\n",
    "        ndofs      = length(nodes)*6\n",
    "        \n",
    "        i=0\n",
    "        for node in nodes\n",
    "            dg=[]\n",
    "            for ii in 0:5\n",
    "                append!(dg,i+ii) \n",
    "            end\n",
    "            i+=6\n",
    "            node[\"degrees_of_freedom\"]=dg\n",
    "        end\n",
    "    end\n",
    "\n",
    "    #######################################################\n",
    "    function solveFea(setup)\n",
    "        // # determine the global matrices\n",
    "        initialize(setup)\n",
    "        \n",
    "        M,K,F,U,ind=get_matrices(setup)\n",
    "        \n",
    "        #println(M)\n",
    "        #println(K)\n",
    "        #println(F)\n",
    "\n",
    "        #evals=eigvals(K,M)\n",
    "        #evecs=eigvecs(K,M)\n",
    "        #frequencies=sqrt.(evals)\n",
    "        X=inv(K)*F\n",
    "        U[setdiff(1:end, ind)]=X\n",
    "\n",
    "        updateDisplacement(setup, U)\n",
    "\n",
    "        # determine the stresses in each element\n",
    "        stresses=get_stresses(setup)\n",
    "    end\n",
    "    #######################################################\n",
    "    displacementFEA=[]\n",
    "    Load=0\n",
    "    topNodesIndices=[]\n",
    "    solveFea(setup)\n",
    "    return displacementFEA,Load,topNodesIndices\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "getYoungsModulus (generic function with 1 method)"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function getYoungsModulus(latticeSize,voxelSize,disp,Load,topNodesIndices)\n",
    "    F=-Load\n",
    "    l0=voxelSize*latticeSize\n",
    "    A=l0*l0\n",
    "\n",
    "    δl1=-mean( x.y for x in disp[topNodesIndices])\n",
    "        \n",
    "    stresses=F/A\n",
    "    strain=δl1/l0\n",
    "    println(\"Load=$Load\")\n",
    "    println(\"stress=$stresses\")\n",
    "\n",
    "    E=stresses/strain \n",
    "\n",
    "    return E\n",
    "end\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "getSetup (generic function with 1 method)"
      ]
     },
     "execution_count": 97,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function getSetup(latticeSize)\n",
    "    setup = Dict()\n",
    "    name=string(\"../json/setupTestUni$latticeSize\",\".json\")\n",
    "    name=string(\"../json/canteliver\",\".json\")\n",
    "    name=string(\"../json/snake\",\".json\")\n",
    "    \n",
    "#     open(\"../json/setupValid2.json\", \"r\") do f\n",
    "#     open(\"../json/setupTest.json\", \"r\") do f\n",
    "    # open(\"../json/trialJulia.json\", \"r\") do f\n",
    "#     open(\"../json/setupTestUni4.json\", \"r\") do f\n",
    "    # open(\"../json/setupChiral.json\", \"r\") do f\n",
    "#     open(\"../json/setupTestCubeUni10.json\", \"r\") do f\n",
    "    open(name, \"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\"]\n",
    "    return setup\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "44\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "104"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latticeSize=1\n",
    "setup=getSetup(latticeSize)\n",
    "displacementFEA,Load,topNodesIndices=fea(setup)\n",
    "topNodesIndices\n",
    "println(length(setup[\"nodes\"]))\n",
    "length(setup[\"edges\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DDisplacements=[[],[],[],[],[]]\n",
    "DDisplacementsFEA=[[],[],[],[],[]]\n",
    "Loads=[0.0,0,0,0,0]\n",
    "EsFEA=[0.0,0,0,0,0]\n",
    "Es=[0.0,0,0,0,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dict{String,Any} with 7 entries:\n",
       "  \"nodes\"      => Any[Dict{String,Any}(\"degrees_of_freedom\"=>Any[0, 1, 2, 3, 4,…\n",
       "  \"voxelSize\"  => 5\n",
       "  \"hierarchical\" => false\n",
       "  \"animation\"  => Dict{String,Any}(\"speed\"=>3,\"exaggeration\"=>20.0,\"showDisplace…\n",
       "  \"viz\"        => Dict{String,Any}(\"colorMap\"=>0,\"colorMaps\"=>Any[Any[Any[0, An…\n",
       "  \"edges\"      => Any[Dict{String,Any}(\"source\"=>0,\"area\"=>5.6644,\"density\"=>0.…\n",
       "  \"ndofs\"      => 264"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "44-element Array{Any,1}:\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.010613496971641401, y:-0.01046280778290931, z:-0.01059129790553449    \n",
       " ⋮                                                                         \n",
       " x:-7.762493357730787e-5, y:2.9952846496351806e-15, z:-0.06302666774019734 \n",
       " x:-0.010957579507126305, y:0.010158333207868608, z:-0.052980102215182984  \n",
       " x:-3.500636791186467e-5, y:6.963266467751947e-5, z:-0.0640476197001775    \n",
       " x:-2.3096342910205877e-5, y:2.9823993790924887e-15, z:-0.06309994346837483\n",
       " x:0.002222557388210866, y:-0.0033974994957394963, z:-0.06642143213658186  \n",
       " x:-0.017667638965804384, y:-0.018868777070104946, z:-0.08187400904630208  \n",
       " x:-0.015410075342964763, y:-0.0002086089466902828, z:-0.08403859895720674 \n",
       " x:0.002222557388210005, y:0.0033974994957447547, z:-0.06642143213658117   \n",
       " x:-0.015582059116631436, y:2.2633004501554508e-15, z:-0.06985538723491842 \n",
       " x:-0.01766763896580353, y:0.018868777070110178, z:-0.08187400904630134    \n",
       " x:-0.015410075342964766, y:0.00020860894669479622, z:-0.08403859895720513 \n",
       " x:-0.015207250272737447, y:2.2496589914265297e-15, z:-0.10102732897327556 "
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "latticeSize=1\n",
    "setup=getSetup(latticeSize)\n",
    "\n",
    "for i in 1:1\n",
    "    displacementFEA,Load,topNodesIndices=fea(setup)\n",
    "end\n",
    "# displacementFEA,Load,topNodesIndices=feaDisplacement(setup,latticeSize)\n",
    "updateDataAndSaveFEA!(setup,\"../json/trialJuliaParallelGPU.json\")\n",
    "displacementFEA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dt: 0.0023068357674191054\n",
      "first timestep took 0.9859363 seconds\n",
      "ran latticeSize 1 with 44 voxels and 104 edges for 50000 time steps took 16.158550899 seconds\n",
      "FEA displacement= -0.10102732897327556,converged displacement= -0.08487774089695395\n"
     ]
    },
    {
     "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=\"clip9800\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip9800)\" 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=\"clip9801\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip9800)\" d=\"\n",
       "M270.627 1425.62 L2352.76 1425.62 L2352.76 121.675 L270.627 121.675  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip9802\">\n",
       "    <rect x=\"270\" y=\"121\" width=\"2083\" height=\"1305\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  325.626,1425.62 325.626,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  718.481,1425.62 718.481,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1111.34,1425.62 1111.34,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1504.19,1425.62 1504.19,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1897.04,1425.62 1897.04,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2289.9,1425.62 2289.9,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  270.627,1248.06 2352.76,1248.06 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  270.627,1028.73 2352.76,1028.73 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  270.627,809.405 2352.76,809.405 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  270.627,590.075 2352.76,590.075 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  270.627,370.745 2352.76,370.745 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  270.627,151.415 2352.76,151.415 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,1425.62 2352.76,1425.62 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,1425.62 270.627,121.675 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  325.626,1425.62 325.626,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  718.481,1425.62 718.481,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1111.34,1425.62 1111.34,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1504.19,1425.62 1504.19,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1897.04,1425.62 1897.04,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2289.9,1425.62 2289.9,1409.97 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,1248.06 295.612,1248.06 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,1028.73 295.612,1028.73 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,809.405 295.612,809.405 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,590.075 295.612,590.075 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,370.745 295.612,370.745 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  270.627,151.415 295.612,151.415 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 325.626, 1479.62)\" x=\"325.626\" y=\"1479.62\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 718.481, 1479.62)\" x=\"718.481\" y=\"1479.62\">100</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 1111.34, 1479.62)\" x=\"1111.34\" y=\"1479.62\">200</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 1504.19, 1479.62)\" x=\"1504.19\" y=\"1479.62\">300</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 1897.04, 1479.62)\" x=\"1897.04\" y=\"1479.62\">400</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 2289.9, 1479.62)\" x=\"2289.9\" y=\"1479.62\">500</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 246.627, 1265.56)\" x=\"246.627\" y=\"1265.56\">-1.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 246.627, 1046.23)\" x=\"246.627\" y=\"1046.23\">-0.75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 246.627, 826.905)\" x=\"246.627\" y=\"826.905\">-0.50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 246.627, 607.575)\" x=\"246.627\" y=\"607.575\">-0.25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 246.627, 388.245)\" x=\"246.627\" y=\"388.245\">0.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 246.627, 168.915)\" x=\"246.627\" y=\"168.915\">0.25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 1311.69, 73.2)\" x=\"1311.69\" y=\"73.2\">1 Voxel Convergence Study</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 1311.69, 1559.48)\" x=\"1311.69\" y=\"1559.48\">timestep</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip9800)\">\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\">displacement</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  329.555,370.745 333.483,370.745 337.412,370.872 341.34,371.472 345.269,372.771 349.197,374.94 353.126,378.087 357.055,382.271 360.983,387.525 364.912,393.861 \n",
       "  368.84,401.279 372.769,409.77 376.697,419.32 380.626,429.908 384.554,441.513 388.483,454.111 392.411,467.674 396.34,482.176 400.269,497.589 404.197,513.883 \n",
       "  408.126,531.027 412.054,548.99 415.983,567.74 419.911,587.241 423.84,607.459 427.768,628.354 431.697,649.887 435.625,672.016 439.554,694.698 443.483,717.884 \n",
       "  447.411,741.527 451.34,765.576 455.268,789.977 459.197,814.674 463.125,839.611 467.054,864.727 470.982,889.961 474.911,915.25 478.839,940.53 482.768,965.735 \n",
       "  486.697,990.798 490.625,1015.65 494.554,1040.23 498.482,1064.46 502.411,1088.29 506.339,1111.63 510.268,1134.42 514.196,1156.61 518.125,1178.12 522.053,1198.88 \n",
       "  525.982,1218.85 529.911,1237.95 533.839,1256.14 537.768,1273.35 541.696,1289.54 545.625,1304.65 549.553,1318.64 553.482,1331.46 557.41,1343.07 561.339,1353.44 \n",
       "  565.267,1362.52 569.196,1370.3 573.125,1376.74 577.053,1381.82 580.982,1385.52 584.91,1387.82 588.839,1388.71 592.767,1388.19 596.696,1386.24 600.624,1382.87 \n",
       "  604.553,1378.07 608.481,1371.85 612.41,1364.23 616.339,1355.21 620.267,1344.8 624.196,1333.04 628.124,1319.93 632.053,1305.5 635.981,1289.78 639.91,1272.81 \n",
       "  643.838,1254.6 647.767,1235.2 651.695,1214.65 655.624,1192.98 659.553,1170.24 663.481,1146.46 667.41,1121.7 671.338,1095.99 675.267,1069.39 679.195,1041.95 \n",
       "  683.124,1013.71 687.052,984.726 690.981,955.049 694.909,924.731 698.838,893.823 702.767,862.379 706.695,830.451 710.624,798.092 714.552,765.355 718.481,732.291 \n",
       "  722.409,698.954 726.338,665.394 730.266,631.773 734.195,598.49 738.123,565.823 742.052,533.996 745.981,503.167 749.909,473.444 753.838,444.906 757.766,417.61 \n",
       "  761.695,391.599 765.623,366.907 769.552,343.558 773.48,321.574 777.409,300.968 781.337,281.753 785.266,263.936 789.195,247.524 793.123,232.52 797.052,218.924 \n",
       "  800.98,206.735 804.909,195.947 808.837,186.554 812.766,178.545 816.694,171.908 820.623,166.624 824.551,162.675 828.48,160.036 832.409,158.681 836.337,158.579 \n",
       "  840.266,159.696 844.194,161.995 848.123,165.433 852.051,169.967 855.98,175.549 859.908,182.128 863.837,189.652 867.766,198.063 871.694,207.303 875.623,217.313 \n",
       "  879.551,228.029 883.48,239.388 887.408,251.325 891.337,263.774 895.265,276.67 899.194,289.944 903.122,303.531 907.051,317.365 910.98,331.38 914.908,345.512 \n",
       "  918.837,359.697 922.765,373.873 926.694,387.981 930.622,401.961 934.551,415.759 938.479,429.321 942.408,442.596 946.336,455.535 950.265,468.094 954.194,480.229 \n",
       "  958.122,491.903 962.051,503.078 965.979,513.722 969.908,523.805 973.836,533.303 977.765,542.191 981.693,550.45 985.622,558.066 989.55,565.024 993.479,571.317 \n",
       "  997.408,576.937 1001.34,581.883 1005.26,586.153 1009.19,589.752 1013.12,592.686 1017.05,594.963 1020.98,596.596 1024.91,597.597 1028.84,597.984 1032.76,597.776 \n",
       "  1036.69,596.992 1040.62,595.656 1044.55,593.793 1048.48,591.429 1052.41,588.59 1056.34,585.306 1060.26,581.607 1064.19,577.525 1068.12,573.089 1072.05,568.334 \n",
       "  1075.98,563.291 1079.91,557.994 1083.84,552.475 1087.76,546.769 1091.69,540.909 1095.62,534.926 1099.55,528.853 1103.48,522.723 1107.41,516.565 1111.34,510.411 \n",
       "  1115.26,504.29 1119.19,498.229 1123.12,492.257 1127.05,486.398 1130.98,480.678 1134.91,475.12 1138.84,469.746 1142.76,464.575 1146.69,459.628 1150.62,454.922 \n",
       "  1154.55,450.471 1158.48,446.29 1162.41,442.391 1166.34,438.786 1170.26,435.483 1174.19,432.49 1178.12,429.812 1182.05,427.454 1185.98,425.419 1189.91,423.708 \n",
       "  1193.83,422.32 1197.76,421.253 1201.69,420.505 1205.62,420.07 1209.55,419.944 1213.48,420.118 1217.41,420.585 1221.33,421.335 1225.26,422.358 1229.19,423.642 \n",
       "  1233.12,425.177 1237.05,426.948 1240.98,428.942 1244.91,431.144 1248.83,433.541 1252.76,436.117 1256.69,438.856 1260.62,441.743 1264.55,444.76 1268.48,447.893 \n",
       "  1272.41,451.125 1276.33,454.439 1280.26,457.82 1284.19,461.251 1288.12,464.717 1292.05,468.201 1295.98,471.69 1299.91,475.168 1303.83,478.62 1307.76,482.034 \n",
       "  1311.69,485.396 1315.62,488.692 1319.55,491.913 1323.48,495.045 1327.41,498.078 1331.33,501.003 1335.26,503.81 1339.19,506.491 1343.12,509.039 1347.05,511.446 \n",
       "  1350.98,513.708 1354.91,515.817 1358.83,517.771 1362.76,519.566 1366.69,521.198 1370.62,522.666 1374.55,523.969 1378.48,525.107 1382.41,526.079 1386.33,526.886 \n",
       "  1390.26,527.531 1394.19,528.015 1398.12,528.342 1402.05,528.515 1405.98,528.538 1409.9,528.416 1413.83,528.154 1417.76,527.757 1421.69,527.232 1425.62,526.585 \n",
       "  1429.55,525.823 1433.48,524.952 1437.4,523.98 1441.33,522.914 1445.26,521.763 1449.19,520.534 1453.12,519.234 1457.05,517.872 1460.98,516.455 1464.9,514.993 \n",
       "  1468.83,513.491 1472.76,511.96 1476.69,510.405 1480.62,508.835 1484.55,507.258 1488.48,505.68 1492.4,504.108 1496.33,502.549 1500.26,501.009 1504.19,499.496 \n",
       "  1508.12,498.014 1512.05,496.568 1515.98,495.165 1519.9,493.809 1523.83,492.505 1527.76,491.255 1531.69,490.065 1535.62,488.937 1539.55,487.875 1543.48,486.88 \n",
       "  1547.4,485.955 1551.33,485.102 1555.26,484.322 1559.19,483.615 1563.12,482.983 1567.05,482.425 1570.98,481.942 1574.9,481.533 1578.83,481.196 1582.76,480.932 \n",
       "  1586.69,480.738 1590.62,480.613 1594.55,480.554 1598.48,480.56 1602.4,480.627 1606.33,480.754 1610.26,480.937 1614.19,481.173 1618.12,481.459 1622.05,481.791 \n",
       "  1625.97,482.167 1629.9,482.581 1633.83,483.031 1637.76,483.512 1641.69,484.022 1645.62,484.555 1649.55,485.109 1653.47,485.679 1657.4,486.261 1661.33,486.852 \n",
       "  1665.26,487.448 1669.19,488.045 1673.12,488.64 1677.05,489.229 1680.97,489.809 1684.9,490.377 1688.83,490.929 1692.76,491.462 1696.69,491.975 1700.62,492.463 \n",
       "  1704.55,492.926 1708.47,493.36 1712.4,493.763 1716.33,494.133 1720.26,494.469 1724.19,494.77 1728.12,495.033 1732.05,495.258 1735.97,495.443 1739.9,495.589 \n",
       "  1743.83,495.694 1747.76,495.758 1751.69,495.78 1755.62,495.762 1759.55,495.702 1763.47,495.602 1767.4,495.461 1771.33,495.281 1775.26,495.062 1779.19,494.805 \n",
       "  1783.12,494.511 1787.05,494.181 1790.97,493.817 1794.9,493.42 1798.83,492.991 1802.76,492.533 1806.69,492.046 1810.62,491.533 1814.55,490.994 1818.47,490.433 \n",
       "  1822.4,489.851 1826.33,489.25 1830.26,488.632 1834.19,487.998 1838.12,487.351 1842.04,486.692 1845.97,486.024 1849.9,485.349 1853.83,484.668 1857.76,483.983 \n",
       "  1861.69,483.296 1865.62,482.608 1869.54,481.922 1873.47,481.239 1877.4,480.561 1881.33,479.888 1885.26,479.223 1889.19,478.566 1893.12,477.92 1897.04,477.285 \n",
       "  1900.97,476.662 1904.9,476.052 1908.83,475.456 1912.76,474.875 1916.69,474.31 1920.62,473.761 1924.54,473.228 1928.47,472.713 1932.4,472.215 1936.33,471.734 \n",
       "  1940.26,471.272 1944.19,470.827 1948.12,470.4 1952.04,469.99 1955.97,469.599 1959.9,469.224 1963.83,468.866 1967.76,468.525 1971.69,468.2 1975.62,467.89 \n",
       "  1979.54,467.595 1983.47,467.315 1987.4,467.048 1991.33,466.794 1995.26,466.553 1999.19,466.322 2003.12,466.102 2007.04,465.892 2010.97,465.691 2014.9,465.497 \n",
       "  2018.83,465.311 2022.76,465.131 2026.69,464.956 2030.62,464.786 2034.54,464.619 2038.47,464.455 2042.4,464.292 2046.33,464.131 2050.26,463.97 2054.19,463.808 \n",
       "  2058.11,463.645 2062.04,463.48 2065.97,463.312 2069.9,463.14 2073.83,462.965 2077.76,462.785 2081.69,462.6 2085.61,462.409 2089.54,462.212 2093.47,462.009 \n",
       "  2097.4,461.799 2101.33,461.581 2105.26,461.357 2109.19,461.124 2113.11,460.884 2117.04,460.636 2120.97,460.379 2124.9,460.115 2128.83,459.843 2132.76,459.562 \n",
       "  2136.69,459.274 2140.61,458.977 2144.54,458.673 2148.47,458.362 2152.4,458.043 2156.33,457.718 2160.26,457.385 2164.19,457.047 2168.11,456.702 2172.04,456.352 \n",
       "  2175.97,455.997 2179.9,455.637 2183.83,455.272 2187.76,454.904 2191.69,454.532 2195.61,454.157 2199.54,453.78 2203.47,453.401 2207.4,453.02 2211.33,452.639 \n",
       "  2215.26,452.257 2219.19,451.875 2223.11,451.493 2227.04,451.113 2230.97,450.734 2234.9,450.357 2238.83,449.982 2242.76,449.61 2246.69,449.242 2250.61,448.877 \n",
       "  2254.54,448.516 2258.47,448.16 2262.4,447.809 2266.33,447.463 2270.26,447.122 2274.18,446.787 2278.11,446.459 2282.04,446.136 2285.97,445.82 2289.9,445.512 \n",
       "  2293.83,445.21 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  329.555,459.378 333.483,459.378 337.412,459.378 341.34,459.378 345.269,459.378 349.197,459.378 353.126,459.378 357.055,459.378 360.983,459.378 364.912,459.378 \n",
       "  368.84,459.378 372.769,459.378 376.697,459.378 380.626,459.378 384.554,459.378 388.483,459.378 392.411,459.378 396.34,459.378 400.269,459.378 404.197,459.378 \n",
       "  408.126,459.378 412.054,459.378 415.983,459.378 419.911,459.378 423.84,459.378 427.768,459.378 431.697,459.378 435.625,459.378 439.554,459.378 443.483,459.378 \n",
       "  447.411,459.378 451.34,459.378 455.268,459.378 459.197,459.378 463.125,459.378 467.054,459.378 470.982,459.378 474.911,459.378 478.839,459.378 482.768,459.378 \n",
       "  486.697,459.378 490.625,459.378 494.554,459.378 498.482,459.378 502.411,459.378 506.339,459.378 510.268,459.378 514.196,459.378 518.125,459.378 522.053,459.378 \n",
       "  525.982,459.378 529.911,459.378 533.839,459.378 537.768,459.378 541.696,459.378 545.625,459.378 549.553,459.378 553.482,459.378 557.41,459.378 561.339,459.378 \n",
       "  565.267,459.378 569.196,459.378 573.125,459.378 577.053,459.378 580.982,459.378 584.91,459.378 588.839,459.378 592.767,459.378 596.696,459.378 600.624,459.378 \n",
       "  604.553,459.378 608.481,459.378 612.41,459.378 616.339,459.378 620.267,459.378 624.196,459.378 628.124,459.378 632.053,459.378 635.981,459.378 639.91,459.378 \n",
       "  643.838,459.378 647.767,459.378 651.695,459.378 655.624,459.378 659.553,459.378 663.481,459.378 667.41,459.378 671.338,459.378 675.267,459.378 679.195,459.378 \n",
       "  683.124,459.378 687.052,459.378 690.981,459.378 694.909,459.378 698.838,459.378 702.767,459.378 706.695,459.378 710.624,459.378 714.552,459.378 718.481,459.378 \n",
       "  722.409,459.378 726.338,459.378 730.266,459.378 734.195,459.378 738.123,459.378 742.052,459.378 745.981,459.378 749.909,459.378 753.838,459.378 757.766,459.378 \n",
       "  761.695,459.378 765.623,459.378 769.552,459.378 773.48,459.378 777.409,459.378 781.337,459.378 785.266,459.378 789.195,459.378 793.123,459.378 797.052,459.378 \n",
       "  800.98,459.378 804.909,459.378 808.837,459.378 812.766,459.378 816.694,459.378 820.623,459.378 824.551,459.378 828.48,459.378 832.409,459.378 836.337,459.378 \n",
       "  840.266,459.378 844.194,459.378 848.123,459.378 852.051,459.378 855.98,459.378 859.908,459.378 863.837,459.378 867.766,459.378 871.694,459.378 875.623,459.378 \n",
       "  879.551,459.378 883.48,459.378 887.408,459.378 891.337,459.378 895.265,459.378 899.194,459.378 903.122,459.378 907.051,459.378 910.98,459.378 914.908,459.378 \n",
       "  918.837,459.378 922.765,459.378 926.694,459.378 930.622,459.378 934.551,459.378 938.479,459.378 942.408,459.378 946.336,459.378 950.265,459.378 954.194,459.378 \n",
       "  958.122,459.378 962.051,459.378 965.979,459.378 969.908,459.378 973.836,459.378 977.765,459.378 981.693,459.378 985.622,459.378 989.55,459.378 993.479,459.378 \n",
       "  997.408,459.378 1001.34,459.378 1005.26,459.378 1009.19,459.378 1013.12,459.378 1017.05,459.378 1020.98,459.378 1024.91,459.378 1028.84,459.378 1032.76,459.378 \n",
       "  1036.69,459.378 1040.62,459.378 1044.55,459.378 1048.48,459.378 1052.41,459.378 1056.34,459.378 1060.26,459.378 1064.19,459.378 1068.12,459.378 1072.05,459.378 \n",
       "  1075.98,459.378 1079.91,459.378 1083.84,459.378 1087.76,459.378 1091.69,459.378 1095.62,459.378 1099.55,459.378 1103.48,459.378 1107.41,459.378 1111.34,459.378 \n",
       "  1115.26,459.378 1119.19,459.378 1123.12,459.378 1127.05,459.378 1130.98,459.378 1134.91,459.378 1138.84,459.378 1142.76,459.378 1146.69,459.378 1150.62,459.378 \n",
       "  1154.55,459.378 1158.48,459.378 1162.41,459.378 1166.34,459.378 1170.26,459.378 1174.19,459.378 1178.12,459.378 1182.05,459.378 1185.98,459.378 1189.91,459.378 \n",
       "  1193.83,459.378 1197.76,459.378 1201.69,459.378 1205.62,459.378 1209.55,459.378 1213.48,459.378 1217.41,459.378 1221.33,459.378 1225.26,459.378 1229.19,459.378 \n",
       "  1233.12,459.378 1237.05,459.378 1240.98,459.378 1244.91,459.378 1248.83,459.378 1252.76,459.378 1256.69,459.378 1260.62,459.378 1264.55,459.378 1268.48,459.378 \n",
       "  1272.41,459.378 1276.33,459.378 1280.26,459.378 1284.19,459.378 1288.12,459.378 1292.05,459.378 1295.98,459.378 1299.91,459.378 1303.83,459.378 1307.76,459.378 \n",
       "  1311.69,459.378 1315.62,459.378 1319.55,459.378 1323.48,459.378 1327.41,459.378 1331.33,459.378 1335.26,459.378 1339.19,459.378 1343.12,459.378 1347.05,459.378 \n",
       "  1350.98,459.378 1354.91,459.378 1358.83,459.378 1362.76,459.378 1366.69,459.378 1370.62,459.378 1374.55,459.378 1378.48,459.378 1382.41,459.378 1386.33,459.378 \n",
       "  1390.26,459.378 1394.19,459.378 1398.12,459.378 1402.05,459.378 1405.98,459.378 1409.9,459.378 1413.83,459.378 1417.76,459.378 1421.69,459.378 1425.62,459.378 \n",
       "  1429.55,459.378 1433.48,459.378 1437.4,459.378 1441.33,459.378 1445.26,459.378 1449.19,459.378 1453.12,459.378 1457.05,459.378 1460.98,459.378 1464.9,459.378 \n",
       "  1468.83,459.378 1472.76,459.378 1476.69,459.378 1480.62,459.378 1484.55,459.378 1488.48,459.378 1492.4,459.378 1496.33,459.378 1500.26,459.378 1504.19,459.378 \n",
       "  1508.12,459.378 1512.05,459.378 1515.98,459.378 1519.9,459.378 1523.83,459.378 1527.76,459.378 1531.69,459.378 1535.62,459.378 1539.55,459.378 1543.48,459.378 \n",
       "  1547.4,459.378 1551.33,459.378 1555.26,459.378 1559.19,459.378 1563.12,459.378 1567.05,459.378 1570.98,459.378 1574.9,459.378 1578.83,459.378 1582.76,459.378 \n",
       "  1586.69,459.378 1590.62,459.378 1594.55,459.378 1598.48,459.378 1602.4,459.378 1606.33,459.378 1610.26,459.378 1614.19,459.378 1618.12,459.378 1622.05,459.378 \n",
       "  1625.97,459.378 1629.9,459.378 1633.83,459.378 1637.76,459.378 1641.69,459.378 1645.62,459.378 1649.55,459.378 1653.47,459.378 1657.4,459.378 1661.33,459.378 \n",
       "  1665.26,459.378 1669.19,459.378 1673.12,459.378 1677.05,459.378 1680.97,459.378 1684.9,459.378 1688.83,459.378 1692.76,459.378 1696.69,459.378 1700.62,459.378 \n",
       "  1704.55,459.378 1708.47,459.378 1712.4,459.378 1716.33,459.378 1720.26,459.378 1724.19,459.378 1728.12,459.378 1732.05,459.378 1735.97,459.378 1739.9,459.378 \n",
       "  1743.83,459.378 1747.76,459.378 1751.69,459.378 1755.62,459.378 1759.55,459.378 1763.47,459.378 1767.4,459.378 1771.33,459.378 1775.26,459.378 1779.19,459.378 \n",
       "  1783.12,459.378 1787.05,459.378 1790.97,459.378 1794.9,459.378 1798.83,459.378 1802.76,459.378 1806.69,459.378 1810.62,459.378 1814.55,459.378 1818.47,459.378 \n",
       "  1822.4,459.378 1826.33,459.378 1830.26,459.378 1834.19,459.378 1838.12,459.378 1842.04,459.378 1845.97,459.378 1849.9,459.378 1853.83,459.378 1857.76,459.378 \n",
       "  1861.69,459.378 1865.62,459.378 1869.54,459.378 1873.47,459.378 1877.4,459.378 1881.33,459.378 1885.26,459.378 1889.19,459.378 1893.12,459.378 1897.04,459.378 \n",
       "  1900.97,459.378 1904.9,459.378 1908.83,459.378 1912.76,459.378 1916.69,459.378 1920.62,459.378 1924.54,459.378 1928.47,459.378 1932.4,459.378 1936.33,459.378 \n",
       "  1940.26,459.378 1944.19,459.378 1948.12,459.378 1952.04,459.378 1955.97,459.378 1959.9,459.378 1963.83,459.378 1967.76,459.378 1971.69,459.378 1975.62,459.378 \n",
       "  1979.54,459.378 1983.47,459.378 1987.4,459.378 1991.33,459.378 1995.26,459.378 1999.19,459.378 2003.12,459.378 2007.04,459.378 2010.97,459.378 2014.9,459.378 \n",
       "  2018.83,459.378 2022.76,459.378 2026.69,459.378 2030.62,459.378 2034.54,459.378 2038.47,459.378 2042.4,459.378 2046.33,459.378 2050.26,459.378 2054.19,459.378 \n",
       "  2058.11,459.378 2062.04,459.378 2065.97,459.378 2069.9,459.378 2073.83,459.378 2077.76,459.378 2081.69,459.378 2085.61,459.378 2089.54,459.378 2093.47,459.378 \n",
       "  2097.4,459.378 2101.33,459.378 2105.26,459.378 2109.19,459.378 2113.11,459.378 2117.04,459.378 2120.97,459.378 2124.9,459.378 2128.83,459.378 2132.76,459.378 \n",
       "  2136.69,459.378 2140.61,459.378 2144.54,459.378 2148.47,459.378 2152.4,459.378 2156.33,459.378 2160.26,459.378 2164.19,459.378 2168.11,459.378 2172.04,459.378 \n",
       "  2175.97,459.378 2179.9,459.378 2183.83,459.378 2187.76,459.378 2191.69,459.378 2195.61,459.378 2199.54,459.378 2203.47,459.378 2207.4,459.378 2211.33,459.378 \n",
       "  2215.26,459.378 2219.19,459.378 2223.11,459.378 2227.04,459.378 2230.97,459.378 2234.9,459.378 2238.83,459.378 2242.76,459.378 2246.69,459.378 2250.61,459.378 \n",
       "  2254.54,459.378 2258.47,459.378 2262.4,459.378 2266.33,459.378 2270.26,459.378 2274.18,459.378 2278.11,459.378 2282.04,459.378 2285.97,459.378 2289.9,459.378 \n",
       "  2293.83,459.378 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip9800)\" d=\"\n",
       "M1853.56 386.635 L2280.76 386.635 L2280.76 205.195 L1853.56 205.195  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1853.56,386.635 2280.76,386.635 2280.76,205.195 1853.56,205.195 1853.56,386.635 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1877.56,265.675 2021.56,265.675 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 2045.56, 283.175)\" x=\"2045.56\" y=\"283.175\">Dynamic</text>\n",
       "</g>\n",
       "<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1877.56,326.155 2021.56,326.155 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip9800)\">\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, 2045.56, 343.655)\" x=\"2045.56\" y=\"343.655\">FEA</text>\n",
       "</g>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "setup=getSetup(latticeSize)\n",
    "numTimeSteps=50000 #30000\n",
    "displacements=[]\n",
    "save=true\n",
    "returnEvery=100\n",
    "runMetavoxelGPU!(setup,numTimeSteps,latticeSize,displacements,returnEvery,true)\n",
    "\n",
    "numTimeStepsRecorded=length(displacements)\n",
    "d=[]\n",
    "dFEA=[]\n",
    "j=length(displacements[end])\n",
    "step=1\n",
    "for i in 1:step:numTimeStepsRecorded\n",
    "    append!(d,displacements[i][j].z)\n",
    "    append!(dFEA,displacementFEA[j].z)\n",
    "end\n",
    "\n",
    "\n",
    "println(\"FEA displacement= $(displacementFEA[j].z),converged displacement= $(displacements[numTimeStepsRecorded][j].z)\")\n",
    "plot(1:step:numTimeStepsRecorded,d,label=\"Dynamic\",xlabel=\"timestep\",ylabel=\"displacement\",title=\"$latticeSize Voxel Convergence Study\")\n",
    "plot!(1:step:numTimeStepsRecorded,dFEA,label=\"FEA\")\n",
    "# savefig(\"4_voxel_convergence\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0\n",
      "0.03141075907812829\n",
      "0.06279051952931337\n",
      "0.09410831331851431\n",
      "0.12533323356430426\n",
      "0.15643446504023087\n",
      "0.1873813145857246\n",
      "0.21814324139654256\n",
      "0.2486898871648548\n",
      "0.2789911060392293\n",
      "0.3090169943749474\n",
      "0.33873792024529137\n",
      "0.3681245526846779\n",
      "0.3971478906347806\n",
      "0.4257792915650727\n",