Parallel FEA 2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra\n",
    "using Plots\n",
    "import JSON\n",
    "using StaticArrays, Rotations\n",
    "# BASED ON https://github.com/jonhiller/Voxelyze"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "simulateParallel (generic function with 3 methods)"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function simulateParallel(setup,numTimeSteps,dt,static=true,saveInterval=10)\n",
    "\tinitialize(setup)\n",
    "    for i in 1:numTimeSteps\n",
    "        # println(\"Timestep:\",i)\n",
    "        doTimeStep(setup,dt,static,i,saveInterval);\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "initialize (generic function with 1 method)"
      ]
     },
     "execution_count": 150,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function initialize(setup)\n",
    "\tnodes      = setup[\"nodes\"]\n",
    "    edges      = setup[\"edges\"]\n",
    "\t# pre-calculate current position\n",
    "\tfor node in nodes\n",
    "        # element=parse(Int,node[\"id\"][2:end])\n",
    "\n",
    "        append!(N_position,[[node[\"position\"][\"x\"] node[\"position\"][\"y\"] node[\"position\"][\"z\"]]])\n",
    "        append!(N_degrees_of_freedom,[node[\"degrees_of_freedom\"]])\n",
    "        append!(N_restrained_degrees_of_freedom, [node[\"restrained_degrees_of_freedom\"]])\n",
    "        append!(N_displacement,[[node[\"displacement\"][\"x\"] node[\"displacement\"][\"y\"] node[\"displacement\"][\"z\"]]])\n",
    "        append!(N_angle,[[node[\"angle\"][\"x\"] node[\"angle\"][\"y\"] node[\"angle\"][\"z\"]]])\n",
    "        append!(N_force,[[node[\"force\"][\"x\"] node[\"force\"][\"y\"] node[\"force\"][\"z\"]]])\n",
    "        append!(N_currPosition,[[node[\"position\"][\"x\"] node[\"position\"][\"y\"] node[\"position\"][\"z\"]]])\n",
    "        append!(N_orient,[[1.0 0.0 0.0 0.0]])#quat\n",
    "        append!(N_linMom,[[0 0 0]])\n",
    "        append!(N_angMom,[[0 0 0]])\n",
    "        append!(N_intForce,[[0 0 0]])\n",
    "        append!(N_intMoment,[[0 0 0]])\n",
    "        append!(N_moment,[[0 0 0]])\n",
    "        \n",
    "        # for dynamic simulations\n",
    "        append!(N_posTimeSteps,[[]])\n",
    "        append!(N_angTimeSteps,[[]])\n",
    "        \n",
    "\tend \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",
    "        append!(E_source,[edge[\"source\"]+1])\n",
    "        append!(E_target,[edge[\"target\"]+1])\n",
    "        append!(E_area,[edge[\"area\"]])\n",
    "#         append!(E_density,[edge[\"density\"]])\n",
    "#         append!(E_stiffness,[edge[\"stiffness\"]])\n",
    "        append!(E_density,[1000])\n",
    "        append!(E_stiffness,[1000000])\n",
    "        append!(E_stress,[0])\n",
    "        append!(E_axis,[axis])\n",
    "        append!(E_currentRestLength,[length])\n",
    "        append!(E_pos2,[[0 0 0]])\n",
    "        append!(E_angle1v,[[0 0 0]])\n",
    "        append!(E_angle2v,[[0 0 0]])\n",
    "        append!(E_angle1,[[1.0 0.0 0.0 0.0]]) #quat\n",
    "        append!(E_angle2,[[1.0 0.0 0.0 0.0]]) #quat\n",
    "        append!(E_currentTransverseStrainSum,[0])\n",
    "        \n",
    "        # for dynamic simulations\n",
    "        append!(E_stressTimeSteps,[[]])\n",
    "        \n",
    "\tend \n",
    "\t\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "doTimeStep (generic function with 1 method)"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function doTimeStep(setup,dt,static,currentTimeStep,saveInterval)\n",
    "    \n",
    "    nodes      = setup[\"nodes\"]\n",
    "    edges      = setup[\"edges\"]\n",
    "    voxCount=size(nodes)[1]\n",
    "    linkCount=size(edges)[1]\n",
    "\n",
    "\tif dt==0 \n",
    "\t\treturn true\n",
    "\telseif (dt<0) \n",
    "\t\tdt = recommendedTimeStep()\n",
    "    end\n",
    "\n",
    "\t# if (collisions) updateCollisions();\n",
    "\tcollisions=false\n",
    "\n",
    "\t# Euler integration:\n",
    "\tDiverged = false\n",
    "    #  for edge in edges\n",
    "    \n",
    "\tfor i in 1:linkCount\n",
    "        # fromNode = nodes[edge[\"source\"]+1]\n",
    "        # toNode = nodes[edge[\"target\"]+1]\n",
    "        # node1 = [fromNode[\"position\"][\"x\"] fromNode[\"position\"][\"y\"] fromNode[\"position\"][\"z\"]]\n",
    "        # node2 = [toNode[\"position\"][\"x\"] toNode[\"position\"][\"y\"] toNode[\"position\"][\"z\"]]\n",
    "        # updateForces(setup,edge,node1,node2,static)# element numbers??\n",
    "        updateForces(setup,i,static)# element numbers??\n",
    "        #  todo: update forces and whatever\n",
    "\t\tif axialStrain(true) > 100\n",
    "\t\t\tDiverged = true; # catch divergent condition! (if any thread sets true we will fail, so don't need mutex...\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    \n",
    "    if Diverged\n",
    "\t\tprintln(\"Divergedd!!!!!\")\n",
    "\t\treturn false\n",
    "    end\n",
    "                \n",
    "\tfor i in 1:voxCount\n",
    "\t\ttimeStep(dt,i,static,currentTimeStep)\n",
    "        # timeStep(dt,node,static,currentTimeStep)\n",
    "\t\tif(!static&& currentTimeStep%saveInterval==0)\n",
    "            append!(N_posTimeSteps[i],[N_displacement[i]])\n",
    "            append!(N_angTimeSteps[i],[N_angle[i]])\n",
    "\n",
    "        end\n",
    "\t\t#  todo: update linMom,angMom, orient and whatever\n",
    "    end\n",
    "\n",
    "\t\n",
    "\tcurrentTimeStep = currentTimeStep+dt\n",
    "\treturn true\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,