{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra\n",
"using Plots\n",
"import JSON\n",
"using StaticArrays, BenchmarkTools\n",
"using Base.Threads\n",
"using CUDAnative\n",
"using CuArrays,CUDAdrv \n",
"# using Test\n",
"import Base: +, * , -, ^"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"struct Vector3\n",
" x::Float64\n",
" y::Float64\n",
" z::Float64\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::Float64\n",
" y::Float64\n",
" z::Float64\n",
" w::Float64\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::Float64\n",
" te2::Float64\n",
" te3::Float64\n",
" te4::Float64\n",
" te5::Float64\n",
" te6::Float64\n",
" te7::Float64\n",
" te8::Float64\n",
" te9::Float64\n",
" te10::Float64\n",
" te11::Float64\n",
" te12::Float64\n",
" te13::Float64\n",
" te14::Float64\n",
" te15::Float64\n",
" te16::Float64\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",
"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::Float64,fy::Float64,fz::Float64,fw::Float64)\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": [
"struct material\n",
" E::Float64\n",
" mass::Float64\n",
" nu::Float64\n",
" rho::Float64\n",
" b::Float64\n",
" h::Float64\n",
" L::Float64\n",
" area::Float64\n",
" I::Float64\n",
" J::Float64\n",
" G::Float64\n",
" a1::Float64\n",
" a2::Float64\n",
" b1::Float64\n",
" b2::Float64\n",
" b3::Float64\n",
" massInverse::Float64\n",
" momentInertiaInverse::Float64\n",
" inertia::Float64\n",
" zeta::Float64\n",
" zetaCollision::Float64\n",
" muStatic::Float64\n",
" muKinetic::Float64\n",
" nomSize::Float64\n",
" sqA1::Float64\n",
" sqA2xIp::Float64 \n",
" sqB1::Float64\n",
" sqB2xFMp::Float64\n",
" sqB3xIp::Float64\n",
" _2xSqMxExS::Float64\n",
" function material()\n",
" E=2000\n",
" mass=10\n",
" nu=0.35\n",
" rho=7.85e-9\n",
" b=2.38\n",
" h=2.38\n",
" L=75\n",
" area=b*h\n",
" I=b*h^3/12\n",
" J=b*h*(b*b+h*h)/12\n",
" G = E * 1 / 3\n",
" a1=E*b*h/L\n",
" a2=G*J/L\n",
" b1=12*E*I/(L^3)\n",
" b2=6*E*I/(L^2)\n",
" b3=2*E*I/(L)\n",
" massInverse=1.0/mass\n",
" momentInertiaInverse=1.92e-6\n",
" inertia=1/momentInertiaInverse\n",
" zeta=1.0\n",
" zetaCollision=0.0\n",
" muStatic= 2.0\n",
" muKinetic= 0.1\n",
" nomSize=1.0\n",
" sqA1=sqrt(a1) \n",
" sqA2xIp=sqrt(a2*L*L/6.0) \n",
" sqB1=sqrt(b1) \n",
" sqB2xFMp=sqrt(b2*L/2) \n",
" sqB3xIp=sqrt(b3*L*L/6.0)\n",
" _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))\n",
" new(E,mass,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,massInverse,momentInertiaInverse,inertia,zeta,zetaCollision,muStatic,muKinetic,nomSize,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,_2xSqMxExS)\n",
" end\n",
" function material(E,mass,nu,rho,b,h,L,momentInertiaInverse,zeta,zetaCollision,muStatic,muKinetic,nomSize)\n",
" area=b*h\n",
" I=b*h^3/12\n",
" J=b*h*(b*b+h*h)/12\n",
" G = E * 1 / 3\n",
" a1=E*b*h/L\n",
" a2=G*J/L\n",
" b1=12*E*I/(L^3)\n",
" b2=6*E*I/(L^2)\n",
" b3=2*E*I/(L)\n",
" massInverse=1.0/mass\n",
" inertia=1/momentInertiaInverse\n",
" sqA1=sqrt(a1) \n",
" sqA2xIp=sqrt(a2*L*L/6.0) \n",
" sqB1=sqrt(b1) \n",
" sqB2xFMp=sqrt(b2*L/2) \n",
" sqB3xIp=sqrt(b3*L*L/6.0)\n",
" _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))\n",
" new(E,mass,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,massInverse,momentInertiaInverse,inertia,zeta,zetaCollision,muStatic,muKinetic,nomSize,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,_2xSqMxExS)\n",
" end\n",
"end\n",
"\n",
"struct voxelMaterial\n",
" E::Float64\n",
" mass::Float64\n",
" nu::Float64\n",
" rho::Float64\n",
" massInverse::Float64\n",
" momentInertiaInverse::Float64\n",
" inertia::Float64\n",
" zeta::Float64\n",
" zetaCollision::Float64\n",
" muStatic::Float64\n",
" muKinetic::Float64\n",
" nomSize::Float64\n",
" _2xSqMxExS::Float64\n",
" function voxelMaterial()\n",
" E=2000\n",
" mass=10\n",
" nu=0.35\n",
" rho=7.85e-9\n",
" massInverse=1.0/mass\n",
" momentInertiaInverse=1.92e-6\n",
" inertia=1/momentInertiaInverse\n",
" zeta=1.0\n",
" zetaCollision=0.0\n",
" muStatic= 2.0\n",
" muKinetic= 0.1\n",
" nomSize=1.0\n",
" _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))\n",
" new(E,mass,nu,rho,massInverse,momentInertiaInverse,inertia,zeta,zetaCollision,muStatic,muKinetic,nomSize,_2xSqMxExS)\n",
" end\n",
" function voxelMaterial(E,mass,nu,rho,momentInertiaInverse,zeta,zetaCollision,muStatic,muKinetic,nomSize)\n",
" massInverse=1.0/mass\n",
" inertia=1/momentInertiaInverse\n",
" _2xSqMxExS = (2.0*sqrt(mass*E*nomSize))\n",
" new(E,mass,nu,rho,massInverse,momentInertiaInverse,inertia,zeta,zetaCollision,muStatic,muKinetic,nomSize,_2xSqMxExS)\n",
" end\n",
"end\n",
"\n",
"struct edgeMaterial\n",
" E::Float64\n",
" nu::Float64\n",
" rho::Float64\n",
" b::Float64\n",
" h::Float64\n",
" L::Float64\n",
" area::Float64\n",
" I::Float64\n",
" J::Float64\n",
" G::Float64\n",
" a1::Float64\n",
" a2::Float64\n",
" b1::Float64\n",
" b2::Float64\n",
" b3::Float64\n",
" sqA1::Float64\n",
" sqA2xIp::Float64 \n",
" sqB1::Float64\n",
" sqB2xFMp::Float64\n",
" sqB3xIp::Float64\n",
" dampingM::Float64\n",
" function edgeMaterial()\n",
" E=2000\n",
" mass=10.0\n",
" nu=0.35\n",
" rho=7.85e-9\n",
" b=2.38\n",
" h=2.38\n",
" L=75\n",
" area=b*h\n",
" I=b*h^3/12\n",
" J=b*h*(b*b+h*h)/12\n",
" G = E * 1 / 3\n",
" a1=E*b*h/L\n",
" a2=G*J/L\n",
" b1=12*E*I/(L^3)\n",
" b2=6*E*I/(L^2)\n",
" b3=2*E*I/(L)\n",
" sqA1=sqrt(a1) \n",
" sqA2xIp=sqrt(a2*L*L/6.0) \n",
" sqB1=sqrt(b1) \n",
" sqB2xFMp=sqrt(b2*L/2) \n",
" sqB3xIp=sqrt(b3*L*L/6.0)\n",
" zeta=1\n",
" dampingM=2.0*sqrt(mass)*zeta\n",
" new(E,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,dampingM)\n",
" end\n",
" function edgeMaterial(E,mass,nu,rho,b,h,L)\n",
" area=b*h\n",
" I=b*h^3/12\n",
" J=b*h*(b*b+h*h)/12\n",
" G = E * 1 / 3\n",
" a1=E*b*h/L\n",
" a2=G*J/L\n",
" b1=12*E*I/(L^3)\n",
" b2=6*E*I/(L^2)\n",
" b3=2*E*I/(L)\n",
" sqA1=sqrt(a1) \n",
" sqA2xIp=sqrt(a2*L*L/6.0) \n",
" sqB1=sqrt(b1) \n",
" sqB2xFMp=sqrt(b2*L/2) \n",
" sqB3xIp=sqrt(b3*L*L/6.0)\n",
" zeta=1\n",
" dampingM=2.0*sqrt(mass)*zeta\n",
" new(E,nu,rho,b,h,L,area,I,J,G,a1,a2,b1,b2,b3,sqA1,sqA2xIp,sqB1,sqB2xFMp,sqB3xIp,dampingM)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"doTimeStep! (generic function with 1 method)"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function doTimeStep!(metavoxel,dt,currentTimeStep)\n",
" # update forces: go through edges, get currentposition from nodes, calc pos2 and update stresses and interior forces of nodes\n",
" run_updateEdges!(dt,\n",
" metavoxel[\"E_sourceGPU\"], \n",
" metavoxel[\"E_targetGPU\"],\n",
" metavoxel[\"E_areaGPU\"],\n",
" metavoxel[\"E_densityGPU\"],\n",
" metavoxel[\"E_stiffnessGPU\"],\n",
" metavoxel[\"E_stressGPU\"],\n",
" metavoxel[\"E_axisGPU\"],\n",
" metavoxel[\"E_currentRestLengthGPU\"],\n",
" metavoxel[\"E_pos2GPU\"],\n",
" metavoxel[\"E_angle1vGPU\"],\n",
" metavoxel[\"E_angle2vGPU\"],\n",
" metavoxel[\"E_angle1GPU\"],\n",
" metavoxel[\"E_angle2GPU\"],\n",
" metavoxel[\"E_intForce1GPU\"],\n",
" metavoxel[\"E_intMoment1GPU\"],\n",
" metavoxel[\"E_intForce2GPU\"],\n",
" metavoxel[\"E_intMoment2GPU\"],\n",
" metavoxel[\"E_dampGPU\"],\n",
" metavoxel[\"E_materialGPU\"],\n",
" metavoxel[\"N_currPositionGPU\"],\n",
" metavoxel[\"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",
" metavoxel[\"N_positionGPU\"], \n",
" metavoxel[\"N_restrainedGPU\"],\n",
" metavoxel[\"N_displacementGPU\"],\n",
" metavoxel[\"N_angleGPU\"],\n",
" metavoxel[\"N_currPositionGPU\"],\n",
" metavoxel[\"N_linMomGPU\"],\n",
" metavoxel[\"N_angMomGPU\"],\n",
" metavoxel[\"N_intForceGPU\"],\n",
" metavoxel[\"N_intMomentGPU\"],\n",
" metavoxel[\"N_forceGPU\"],\n",
" metavoxel[\"N_momentGPU\"],\n",
" metavoxel[\"N_orientGPU\"],\n",
" metavoxel[\"N_edgeIDGPU\"], \n",
" metavoxel[\"N_edgeFirstGPU\"],\n",
" metavoxel[\"N_materialGPU\"], \n",
" metavoxel[\"E_intForce1GPU\"],\n",
" metavoxel[\"E_intMoment1GPU\"],\n",
" metavoxel[\"E_intForce2GPU\"],\n",
" metavoxel[\"E_intMoment2GPU\"])\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"run_updateEdges! (generic function with 1 method)"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function updateEdges!(dt,E_source,E_target,E_area,E_density,E_stiffness,E_stress,E_axis,\n",
" E_currentRestLength,E_pos2,E_angle1v,E_angle2v,\n",
" E_angle1,E_angle2,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,E_damp,E_material,\n",
" 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=Vector3(E_pos2[i].x,E_pos2[i].y,E_pos2[i].z) #?copy?\n",
" @inbounds oldAngle1v = Vector3(E_angle1v[i].x,E_angle1v[i].y,E_angle1v[i].z)\n",
" @inbounds oldAngle2v = Vector3(E_angle2v[i].x,E_angle2v[i].y,E_angle2v[i].z)# remember the positions/angles from last timestep to calculate velocity\n",
" \n",
" \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",
" \n",
" @inbounds dPos2 = Vector3(0.5,0.5,0.5) * (E_pos2[i]-oldPos2) #deltas for local damping. velocity at center is half the total velocity\n",
" @inbounds dAngle1 = Vector3(0.5,0.5,0.5) *(E_angle1v[i]-oldAngle1v)\n",
" @inbounds dAngle2 = Vector3(0.5,0.5,0.5) *(E_angle2v[i]-oldAngle2v)\n",
" \n",
" \n",
" @inbounds strain=(E_pos2[i].x/E_currentRestLength[i])\n",
" \n",
" positiveEnd=true\n",
" if axialStrain( positiveEnd,strain)>100.0\n",
" diverged=true\n",
" @cuprintln(\"DIVERGED!!!!!!!!!!\")\n",
" return \n",
" end\n",
" \n",
" @inbounds E = E_stiffness[i]\n",
" \n",
" \n",
" \n",
" @inbounds l = E_currentRestLength[i]\n",
" \n",
"\n",
" \n",
" @inbounds nu = convert(Float64,E_material[i].nu)\n",
" \n",
" # Cross Section inputs, must be floats \n",
" @inbounds E = convert(Float64,E_material[i].E) # MPa\n",
" @inbounds h = convert(Float64,E_material[i].h) # mm\n",
" @inbounds b = convert(Float64,E_material[i].b) # mm\n",
" \n",
"\n",
" @inbounds a1=convert(Float64,E_material[i].a1)\n",
" @inbounds a2=convert(Float64,E_material[i].a2)\n",
" @inbounds b1=convert(Float64,E_material[i].b1)\n",
" @inbounds b2=convert(Float64,E_material[i].b2)\n",
" @inbounds b3=convert(Float64,E_material[i].b3)\n",
" \n",
" \n",
"\n",
" \n",
" #inbounds currentTransverseArea=25.0 #?? change!!!!!! E_area[i]\n",
" @inbounds currentTransverseArea= b*h\n",
" @inbounds _stress=updateStrain(strain,E)\n",
" \n",
" #@inbounds currentTransverseArea= E_area[i]\n",
" #@inbounds _stress=updateStrain(strain,E_stiffness[i])\n",
" \n",
" @inbounds E_stress[i]=_stress\n",
" \n",
" #@cuprintln(\"_stress $_stress\")\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",
" \n",
" x=convert(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,z)\n",
" \n",
" # Use Curstress instead of -a1*Pos2.x to account for non-linear deformation \n",
" forceNeg = Vector3(x,y,z)\n",
" \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(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,z)\n",
" momentNeg = Vector3(x,y,z)\n",
" \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(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,z)\n",
" momentPos = Vector3(x,y,z)\n",
" \n",
" \n",
" ### damping\n",
" @inbounds if E_damp[i] #first pass no damping\n",
" @inbounds sqA1 =convert(Float64,E_material[i].sqA1)\n",
" @inbounds sqA2xIp =convert(Float64,E_material[i].sqA2xIp)\n",
" @inbounds sqB1 =convert(Float64,E_material[i].sqB1)\n",
" @inbounds sqB2xFMp =convert(Float64,E_material[i].sqB2xFMp)\n",
" @inbounds sqB3xIp =convert(Float64,E_material[i].sqB3xIp)\n",
" \n",
" dampingMultiplier=Vector3(28099.3,28099.3,28099.3) # 2*mat->_sqrtMass*mat->zetaInternal/previousDt;?? todo link to material\n",
" \n",
" zeta=1\n",
" dampingM= convert(Float64,E_material[i].dampingM)/dt\n",
" dampingMultiplier=Vector3(dampingM,dampingM,dampingM)\n",
" \n",
" posCalc=Vector3(sqA1*dPos2.x, sqB1*dPos2.y - sqB2xFMp*(dAngle1.z+dAngle2.z),sqB1*dPos2.z + sqB2xFMp*(dAngle1.y+dAngle2.y))\n",
" \n",
" \n",
" forceNeg =forceNeg + (dampingMultiplier*posCalc);\n",
" forcePos =forcePos - (dampingMultiplier*posCalc);\n",
"\n",
" momentNeg -= Vector3(0.5,0.5,0.5)*dampingMultiplier*Vector3(-sqA2xIp*(dAngle2.x - dAngle1.x),\n",
" sqB2xFMp*dPos2.z + sqB3xIp*(2*dAngle1.y + dAngle2.y),\n",
" -sqB2xFMp*dPos2.y + sqB3xIp*(2*dAngle1.z + dAngle2.z));\n",
" momentPos -= Vector3(0.5,0.5,0.5)*dampingMultiplier*Vector3(sqA2xIp*(dAngle2.x - dAngle1.x),\n",
" sqB2xFMp*dPos2.z + sqB3xIp*(dAngle1.y + 2*dAngle2.y),\n",
" -sqB2xFMp*dPos2.y + sqB3xIp*(dAngle1.z + 2*dAngle2.z));\n",
"\n",
" else\n",
" @inbounds E_damp[i]=true \n",
" end\n",
"\n",
" smallAngle=true\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",
" @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",
"\n",
" @inbounds E_intForce1[i] =forceNeg\n",
" @inbounds E_intForce2[i] =forcePos\n",
" \n",
"\n",
"\n",
" @inbounds x= momentNeg.x\n",
" @inbounds y= momentNeg.y\n",
" @inbounds z= momentNeg.z \n",
" x=convert(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,z)\n",
" \n",
" @inbounds E_intMoment1[i]=Vector3(x,y,z)\n",
"\n",
" @inbounds x= momentNeg.x\n",
" @inbounds y= momentNeg.y\n",
" @inbounds z= momentNeg.z\n",
" x=convert(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,z)\n",
" \n",
" @inbounds E_intMoment2[i]=Vector3(x,y,z)\n",
" \n",
" #x=E_pos2[i].x*10000000000\n",
" #y=E_pos2[i].y*10000000000\n",
" #z=E_pos2[i].z*10000000000\n",
" #@cuprintln(\"pos2 x $x, y $y, z $z \")\n",
" ##x=E_intMoment2[i].x*10000000000\n",
" #y=E_intMoment2[i].y*10000000000\n",
" #z=E_intMoment2[i].z*10000000000\n",
" #@cuprintln(\"E_intMoment2 x $x, y $y, z $z \")\n",
"\n",
" \n",
" \n",
" end\n",
" return\n",
"end\n",
"\n",
"function run_updateEdges!(dt,E_source,E_target,E_area,E_density,E_stiffness,\n",
" E_stress,E_axis,E_currentRestLength,E_pos2,E_angle1v,E_angle2v,\n",
" E_angle1,E_angle2,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,\n",
" E_damp,E_material,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!(dt,E_source,E_target,E_area,E_density,\n",
" E_stiffness,E_stress,E_axis,E_currentRestLength,E_pos2,E_angle1v,\n",
" E_angle2v,E_angle1,E_angle2,E_intForce1,E_intMoment1,E_intForce2,\n",
" E_intMoment2,E_damp,E_material,N_currPosition,N_orient)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"floorForce! (generic function with 1 method)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function floorPenetration(y)\n",
" floor=0.0\n",
" p=0.0\n",
" if(y<floor)\n",
" p=floor-y\n",
" end\n",
" return p\n",
"end\n",
"#Returns the interference (in meters) between the collision envelope of this voxel and the floor at Z=0. Positive numbers correspond to interference. If the voxel is not touching the floor 0 is returned.\n",
"\n",
"function penetrationStiffness(E,nomSize)\n",
" return (2.0*E*nomSize)\n",
"end \n",
"#!< returns the stiffness with which this voxel will resist penetration. This is calculated according to E*A/L with L = voxelSize/2.\n",
"\n",
"\n",
"function globalDampingTranslateC(_2xSqMxExS,zetaGlobal)\n",
" return zetaGlobal*_2xSqMxExS;\n",
"end #!< Returns the global material damping coefficient (translation)\n",
"\n",
"function collisionDampingTranslateC(_2xSqMxExS,zetaCollision)\n",
"# _2xSqIxExSxSxS = (2.0f*CUDAnative.sqrt(_momentInertia*E*nomSize*nomSize*nomSize));\n",
" return zetaCollision*_2xSqMxExS;\n",
"end #!< Returns the global material damping coefficient (translation)\n",
"\n",
"function floorForce!(dt,pTotalForce,y,linMom,FloorStaticFriction,N_material)\n",
" E=convert(Float64,N_material.E)\n",
" nomSize=convert(Float64,N_material.nomSize)\n",
" mass=convert(Float64,N_material.mass)\n",
" massInverse=convert(Float64,N_material.massInverse)\n",
" muStatic=convert(Float64,N_material.muStatic)\n",
" muKinetic=convert(Float64,N_material.muKinetic)\n",
" _2xSqMxExS=convert(Float64,N_material._2xSqMxExS)\n",
" zetaCollision=convert(Float64,N_material.zetaCollision)\n",
" \n",
" CurPenetration = floorPenetration(convert(Float64,y)); #for now use the average.\n",
"\n",
"\n",
" if (CurPenetration>=0.0)\n",
" vel = linMom*Vector3((massInverse),(massInverse),(massInverse)) #Returns the 3D velocity of this voxel in m/s (GCS)\n",
" horizontalVel= Vector3(convert(Float64,vel.x), 0.0, convert(Float64,vel.z));\n",
" \n",
" normalForce = penetrationStiffness(E,nomSize)*CurPenetration;\n",
" pTotalForce=Vector3( pTotalForce.x, convert(Float64,pTotalForce.y) + normalForce - collisionDampingTranslateC(_2xSqMxExS,zetaCollision)*convert(Float64,vel.y),pTotalForce.z)\n",
" #in the z direction: k*x-C*v - spring and damping\n",
"\n",
" if (FloorStaticFriction) #If this voxel is currently in static friction mode (no lateral motion) \n",
" # assert(horizontalVel.Length2() == 0);\n",
" surfaceForceSq = convert(Float64,(pTotalForce.x*pTotalForce.x + pTotalForce.z*pTotalForce.z)); #use squares to avoid a square root\n",
" frictionForceSq = (muStatic*normalForce)*(muStatic*normalForce);\n",
"\n",
"\n",
" if (surfaceForceSq > frictionForceSq) \n",
" FloorStaticFriction=false; #if we're breaking static friction, leave the forces as they currently have been calculated to initiate motion this time step\n",
" end\n",
"\n",
" else #even if we just transitioned don't process here or else with a complete lack of momentum it'll just go back to static friction\n",
" #add a friction force opposing velocity according to the normal force and the kinetic coefficient of friction\n",
" leng=CUDAnative.sqrt((convert(Float64,vel.x) * convert(Float64,vel.x)) + (0.0 * 0.0) + (convert(Float64,vel.z) * convert(Float64,vel.z)))\n",
" if(leng>0)\n",
" horizontalVel= Vector3(convert(Float64,vel.x)/(leng),0.0,convert(Float64,vel.z)/(leng))\n",
" else\n",
" horizontalVel= Vector3(convert(Float64,vel.x)*(leng),0.0,convert(Float64,vel.z)*(leng))\n",
"\n",
" end \n",
" pTotalForce = pTotalForce- Vector3(muKinetic*normalForce,muKinetic*normalForce,muKinetic*normalForce) * horizontalVel;\n",
" end\n",
"\n",
" else \n",
" FloorStaticFriction=false;\n",
" end\n",
" \n",
" \n",
" \n",
" return pTotalForce,FloorStaticFriction\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"run_updateNodes! (generic function with 1 method)"
]
},
"execution_count": 24,
"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,N_material,\n",
" 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",
" if false\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",
" \n",
" isFloorEnabled=true;#todo make these arguments ??\n",
" static=false;#todo make these arguments ??\n",
" FloorStaticFriction=false; #todo make these arguments ??\n",
" \n",
" #get properties\n",
" \n",
" @inbounds E=N_material[i].E\n",
" @inbounds nomSize=N_material[i].nomSize\n",
" @inbounds mass=N_material[i].mass\n",
" @inbounds massInverse=N_material[i].massInverse\n",
" \n",
" @inbounds curForce = force(N_intForce[i],N_orient[i],N_force[i],static,currentTimeStep,mass)\n",
" \n",
" \n",
" \n",
" \n",
" fricForce = Vector3(curForce.x,curForce.y,curForce.z);\n",
" if (isFloorEnabled)\n",
" @inbounds curForce,FloorStaticFriction=floorForce!(dt,curForce,N_currPosition[i].y,\n",
" Vector3(N_linMom[i].x,N_linMom[i].y ,N_linMom[i].z),FloorStaticFriction,N_material[i])\n",
" \n",
" end\n",
" \n",
" fricForce = curForce-fricForce;\n",
" \n",
" \n",
" #########################################\n",
" \n",
" @inbounds N_intForce[i]=Vector3(0,0,0)\n",
" \n",
" \n",
" @inbounds N_linMom[i]=N_linMom[i]+curForce*Vector3(dt,dt,dt) #todo make sure right\n",
" @inbounds translate=N_linMom[i]*Vector3((dt*massInverse),(dt*massInverse),(dt*massInverse)) # ??massInverse\n",
" \n",
" \n",
" ############################\n",
" \n",
" \n",
" #we need to check for friction conditions here (after calculating the translation) and stop things accordingly\n",
" @inbounds if (isFloorEnabled && floorPenetration( N_currPosition[i].y)>= 0.0 && !static)#we must catch a slowing voxel here since it all boils down to needing access to the dt of this timestep.\n",
" work =convert(Float64,fricForce.x*translate.x + fricForce.z*translate.z); #F dot disp\n",
" @inbounds hKe = 0.5*massInverse*convert(Float64,(N_linMom[i].x*N_linMom[i].x + N_linMom[i].z*N_linMom[i].z)); #horizontal kinetic energy\n",
"\n",
" if((hKe + work) <= 0.0) \n",
" FloorStaticFriction=true; #this checks for a change of direction according to the work-energy principle\n",
" end\n",
"\n",
" if(FloorStaticFriction)\n",
" #if we're in a state of static friction, zero out all horizontal motion\n",
" N_linMom[i]=Vector3(0.0 ,convert(Float64,N_linMom[i].y) ,0.0)\n",
" translate=Vector3(0.0 ,convert(Float64,translate.y) ,0.0)\n",
" end\n",
" else\n",
" FloorStaticFriction=false\n",
" end\n",
" \n",
" \n",
" @inbounds N_currPosition[i]=N_currPosition[i]+translate\n",
" @inbounds N_displacement[i]=N_displacement[i]+translate\n",
" \n",
" \n",
" \n",
" # Rotation\n",
" @inbounds curMoment = moment(N_intMoment[i],N_orient[i],N_moment[i]) \n",
" \n",
" \n",
" \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",
" \n",
" @inbounds momentInertiaInverse=N_material[i].momentInertiaInverse\n",
" \n",
" @inbounds temp=FromRotationVector(N_angMom[i]*Vector3((dt*momentInertiaInverse),(dt*momentInertiaInverse),(dt*momentInertiaInverse)))\n",
"\n",
" \n",
" @inbounds N_orient[i]=multiplyQuaternions(temp,N_orient[i])\n",
" \n",
" \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",
" N_linMom,N_angMom,N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,N_material,\n",
" E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\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",
" N_angle,N_currPosition,N_linMom,N_angMom,N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,N_material,\n",
" E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"orientLink! (generic function with 1 method)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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",
" #x=pos2.x*10000000000\n",
" #y=pos2.y*10000000000\n",
" #z=pos2.z*10000000000\n",
" #@cuprintln(\"pos2 x $x, y $y, z $z \")\n",
" \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",
" \n",
" \n",
" # pos2,angle1v,angle2v,angle1,angle2,\n",
" return pos2,angle1v,angle2v,angle1,angle2,totalRot\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"toAxisOriginalQuat (generic function with 1 method)"
]
},
"execution_count": 26,
"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",
" vector=normalizeVector3(axis)\n",
" q=setFromUnitVectors(vector,xaxis)\n",
" \n",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" return Vector3(v.x,v.y,v.z)\n",
"end\n",
"\n",
"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",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" return Vector3(v.x,v.y,v.z)\n",
"end\n",
"\n",
"function toAxisXQuat(pQ::Quaternion,axis::Vector3)\n",
" xaxis=Vector3(1.0,0.0,0.0)\n",
"\n",
" vector=normalizeVector3(axis)\n",
"\n",
" q=setFromUnitVectors(vector,xaxis)\n",
" \n",
" pV=Vector3(pQ.x,pQ.y,pQ.z)\n",
" \n",
" v=applyQuaternion1( pV ,q )\n",
"\n",
" return Quaternion(v.x,v.y,v.z,1.0)\n",
" \n",
"end\n",
"\n",
"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",
" pV=Vector3(pQ.x,pQ.y,pQ.z)\n",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" return Quaternion(v.x,v.y,v.z,1.0)\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"applyQuaternion1 (generic function with 1 method)"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function setFromUnitVectors(vFrom::Vector3, vTo::Vector3)\n",
" # assumes direction vectors vFrom and vTo are normalized\n",
" EPS = 0.000000001;\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.0\n",
" qw = r\n",
" else \n",
" qx = 0.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) ? 0.0 : qx\n",
" qy= (qy==-0.0) ? 0.0 : qy\n",
" qz= (qz==-0.0) ? 0.0 : qz\n",
" qw= (qw==-0.0) ? 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(Float64,mm)#??????????????????? todo check later\n",
" \n",
" l=CUDAnative.sqrt(mm)\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",
" \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\n",
"\n",
"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.0 - ( yy + zz ) ) * sx\n",
" te2 = ( xy + wz ) * sx\n",
" te3 = ( xz - wy ) * sx\n",
" te4 = 0.0\n",
"\n",
" te5 = ( xy - wz ) * sy\n",
" te6 = ( 1.0 - ( xx + zz ) ) * sy\n",
" te7 = ( yz + wx ) * sy\n",
" te8 = 0.0\n",
"\n",
" te9 = ( xz + wy ) * sz\n",
" te10 = ( yz - wx ) * sz\n",
" te11 = ( 1.0 - ( xx + yy ) ) * sz\n",
" te12 = 0.0\n",
"\n",
" te13 = 0.0 #position.x;\n",
" te14 = 0.0 #position.y;\n",
" te15 = 0.0 #position.z;\n",
" te16 = 1.0\n",
" \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\n",
"\n",
"function setFromRotationMatrix(m::RotationMatrix)\n",
"\n",
" m11 = convert(Float64,m.te1 )\n",
" m12 = convert(Float64,m.te5 )\n",
" m13 = convert(Float64,m.te9 )\n",
" m21 = convert(Float64,m.te2 )\n",
" m22 = convert(Float64,m.te6 )\n",
" m23 = convert(Float64,m.te10)\n",
" m31 = convert(Float64,m.te3 )\n",
" m32 = convert(Float64,m.te7 )\n",
" m33 = convert(Float64,m.te11)\n",
" \n",
"\n",
" y = CUDAnative.asin( clamp( m13, -1.0, 1.0 ) ) ##check if has to be changed to cuda\n",
"\n",
" if ( abs( m13 ) < 0.9999999999 ) \n",
" \n",
" x = CUDAnative.atan2( - m23, m33 )\n",
" z = CUDAnative.atan2( - m12, m11 )#-m12, m11\n",
"\n",
"\n",
" else\n",
"\n",
" x = CUDAnative.atan2( m32, m22 )\n",
" z = 0.0;\n",
"\n",
" end\n",
" \n",
" \n",
" return Vector3(x,y,z)\n",
" \n",
"end\n",
"\n",
"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\n",
"\n",
"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",
" d=15\n",
"\n",
" return Vector3(xx,yy,zz)\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"multiplyQuaternions (generic function with 1 method)"
]
},
"execution_count": 28,
"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(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,z)\n",
" w=convert(Float64,w)\n",
" return Quaternion(x,y,z,w)\n",
"end\n",
"\n",
"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",
"\n",
"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",
"\n",
"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(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,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(Float64,x)\n",
" y=convert(Float64,y)\n",
" z=convert(Float64,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",
"\n",
"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(Float64,qx)\n",
" qy=convert(Float64,qy)\n",
" qz=convert(Float64,qz)\n",
" qw=convert(Float64,qw)\n",
" \n",
" return Quaternion(qx,qy,qz,qw)\n",
"end\n",
"\n",
"function multiplyQuaternions(q::Quaternion,f::Quaternion)\n",
" x=q.x\n",
" y=q.y\n",
" z=q.z\n",
" w=q.w\n",
" x1=w*f.x + x*f.w + y*f.z - z*f.y \n",
" y1=w*f.y - x*f.z + y*f.w + z*f.x\n",
" z1=w*f.z + x*f.y - y*f.x + z*f.w\n",
" w1=w*f.w - x*f.x - y*f.y - z*f.z\n",
"\n",
"\treturn Quaternion(x1,y1,z1,w1 ); #!< overload quaternion multiplication.\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"axialStrain (generic function with 1 method)"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function updateStrain( axialStrain,E) # ?from where strain\n",
" strain = axialStrain # redundant?\n",
" currentTransverseStrainSum=0.0 # ??? todo\n",
" linear=true\n",
" maxStrain=1000000000000000;# ?? 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.0-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\n",
"\n",
"function stress( strain , E ) #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",
" 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 \n",
"\n",
"function axialStrain( positiveEnd,strain)\n",
"\t#strainRatio = pVPos->material()->E/pVNeg->material()->E;\n",
"\tstrainRatio=1.0;\n",
"\treturn positiveEnd ? 2.0 *strain*strainRatio/(1.0+strainRatio) : 2.0*strain/(1.0+strainRatio)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"moment (generic function with 1 method)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function force(N_intForce,N_orient,N_force,static,currentTimeStep,mass) \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\n",
" #massInverse=1.0/mass #\n",
" #vel = linMom*Vector3((massInverse),(massInverse),(massInverse))\n",
" #totalForce =totalForce- vel*globalDampingTranslateC(); #global damping f-cv\n",
" \n",
" gravity=true;\n",
" grav=-mass*9.80665*1.0;\n",
" if(gravity && !static)\n",
" totalForce =totalForce +Vector3(0,grav,0);\n",
" end\n",
"\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\n",
"\n",
"\n",
"function moment(intMoment,orient,moment) \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+intMoment\n",
" \n",
" \n",
"\n",
" totalMoment = RotateVec3D(orient,totalMoment);\n",
" \n",
" \n",
"\n",
" totalMoment=totalMoment+moment\n",
"\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": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"updateDataAndSave! (generic function with 1 method)"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function updateDataAndSave!(metavoxel,setup,fileName,displacements)\n",
" nodes = setup[\"nodes\"]\n",
" edges = setup[\"edges\"]\n",
" \n",
" setup[\"animation\"][\"showDisplacement\"]=true\n",
" voxCount=size(nodes)[1]\n",
" linkCount=size(edges)[1]\n",
" \n",
" N_displacement=Array(metavoxel[\"N_displacementGPU\"])\n",
" N_angle=Array(metavoxel[\"N_angleGPU\"])\n",
" E_stress=Array(metavoxel[\"E_stressGPU\"])\n",
" \n",
" setup[\"viz\"][\"maxStress\"]=maximum(E_stress)\n",
" setup[\"viz\"][\"minStress\"]=minimum(E_stress)\n",
" setup[\"animation\"][\"exaggeration\"]=3.0\n",
"\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[\"posTimeSteps\"]=[]\n",
" node[\"angTimeSteps\"]=[]\n",
" node[\"degrees_of_freedom\"]=\"\"\n",
" \n",
" \n",
" node[\"displacement\"][\"x\"]=N_displacement[i].x/15\n",
" node[\"displacement\"][\"y\"]=N_displacement[i].y/15\n",
" node[\"displacement\"][\"z\"]=N_displacement[i].z/15\n",
" \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",
" \n",
"# for j in 1:length(displacements)\n",
"# i=1 \n",
"# for node in nodes\n",
"# d=displacements[j][i]\n",
"# dis = Dict{String, Float64}(\"x\" => d.x/15, \"y\" => d.y/15,\"z\" => d.z/15) \n",
"# append!(node[\"posTimeSteps\"],[dis])\n",
"# ang = Dict{String, Float64}(\"x\" => N_angle[i].x, \"y\" => N_angle[i].y,\"z\" => N_angle[i].z) \n",
"# append!(node[\"angTimeSteps\"],[ang])\n",
"# i=i+1\n",
"# end\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"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"runMetavoxelGPU! (generic function with 1 method)"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function runMetavoxelGPU!(setup,numTimeSteps,latticeSize,displacements,returnEvery,save)\n",
" function initialize!(setup)\n",
" nodes = setup[\"nodes\"]\n",
" edges = setup[\"edges\"]\n",
"\n",
" i=1\n",
" # pre-calculate current position\n",
" for node in nodes\n",
" node[\"position\"][\"y\"]=node[\"position\"][\"y\"]+10.0\n",
" # element=parse(Int,node[\"id\"][2:end])\n",
" N_position[i]=Vector3(node[\"position\"][\"x\"]*15.0,node[\"position\"][\"y\"]*15.0,node[\"position\"][\"z\"]*15.0)\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\"]*15,node[\"displacement\"][\"y\"]*15,node[\"displacement\"][\"z\"]*15)\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\"]*15.0,node[\"position\"][\"y\"]*15.0,node[\"position\"][\"z\"]*15.0)\n",
"\n",
" # for dynamic simulations\n",
" # append!(N_posTimeSteps,[[]])\n",
" # append!(N_angTimeSteps,[[]])\n",
"\n",
" i=i+1\n",
" end \n",
"\n",
" i=1\n",
" # pre-calculate the axis\n",
" for 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\"]*15.0 fromNode[\"position\"][\"y\"]*15.0 fromNode[\"position\"][\"z\"]*15.0]\n",
" node2 = [toNode[\"position\"][\"x\"]*15.0 toNode[\"position\"][\"y\"]*15.0 toNode[\"position\"][\"z\"]*15.0]\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 #?????? todo change \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",
" end \n",
" end\n",
" function simulateParallel!(metavoxel,numTimeSteps,dt,returnEvery)\n",
" # initialize(setup)\n",
"\n",
" for i in 1:numTimeSteps\n",
" #println(\"Timestep:\",i)\n",
" doTimeStep!(metavoxel,dt,i)\n",
" if(mod(i,returnEvery)==0)\n",
" append!(displacements,[Array(metavoxel[\"N_displacementGPU\"])])\n",
" end\n",
" end\n",
" end\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",
" maxNumEdges=10\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_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_material=fill(voxelMaterial(),voxCount)\n",
" #voxelMaterial(E,mass,nu,rho,momentInertiaInverse,zeta,zetaCollision,muStatic,muKinetic,nomSize)\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",
" E_damp=fill(false,linkCount)\n",
" E_material=fill(edgeMaterial(),linkCount)\n",
" E_currentTransverseStrainSum=fill(0.0f0,linkCount)# TODO remove ot incorporate\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_materialGPU= CuArray(N_material) \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_dampGPU= CuArray(E_damp)\n",
" E_materialGPU= CuArray(E_material)\n",
"\n",
"\n",
"\n",
" #########################################\n",
" metavoxel = Dict(\n",
" \"N_positionGPU\" => N_positionGPU, \n",
" \"N_restrainedGPU\" => N_restrainedGPU, \n",
" \"N_displacementGPU\" => N_displacementGPU,\n",
" \"N_angleGPU\" => N_angleGPU, \n",
" \"N_currPositionGPU\" => N_currPositionGPU,\n",
" \"N_linMomGPU\" => N_linMomGPU, \n",
" \"N_angMomGPU\" => N_angMomGPU, \n",
" \"N_intForceGPU\" => N_intForceGPU, \n",
" \"N_intMomentGPU\" => N_intMomentGPU, \n",
" \"N_momentGPU\" => N_momentGPU, \n",
" \"N_forceGPU\" => N_forceGPU, \n",
" \"N_orientGPU\" => N_orientGPU, \n",
" \"N_edgeIDGPU\" => N_edgeIDGPU, \n",
" \"N_edgeFirstGPU\" => N_edgeFirstGPU,\n",
" \"N_materialGPU\"=> N_materialGPU,\n",
" \"E_sourceGPU\" =>E_sourceGPU, \n",
" \"E_targetGPU\" =>E_targetGPU, \n",
" \"E_areaGPU\" =>E_areaGPU, \n",
" \"E_densityGPU\" =>E_densityGPU, \n",
" \"E_stiffnessGPU\" =>E_stiffnessGPU, \n",
" \"E_stressGPU\" =>E_stressGPU, \n",
" \"E_axisGPU\" =>E_axisGPU, \n",
" \"E_currentRestLengthGPU\" =>E_currentRestLengthGPU, \n",
" \"E_pos2GPU\" =>E_pos2GPU, \n",
" \"E_angle1vGPU\" =>E_angle1vGPU, \n",
" \"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",
" \"E_materialGPU\" =>E_materialGPU\n",
" )\n",
"\n",
" #########################################\n",
" \n",
"\n",
" dt=0.0251646\n",
" E = 2000 # MPa\n",
" s=2.38\n",
" mass=10 \n",
" \n",
" \n",
" \n",
" MaxFreq2=E*s/mass\n",
" dt= 1/(6.283185*sqrt(MaxFreq2))\n",
"# dt=0.0001646\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/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": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"runMetavoxelGPULive! (generic function with 1 method)"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function runMetavoxelGPULive!(setup,maxNumTimeSteps,saveEvery,maxNumFiles)\n",
" function initialize!(setup)\n",
" nodes = setup[\"nodes\"]\n",
" edges = setup[\"edges\"]\n",
"\n",
" i=1\n",
" # pre-calculate current position\n",
" for node in nodes\n",
" # element=parse(Int,node[\"id\"][2:end])\n",
" node[\"position\"][\"y\"]=node[\"position\"][\"y\"]+10.0\n",
" N_position[i]=Vector3(node[\"position\"][\"x\"]*15.0,node[\"position\"][\"y\"]*15.0,node[\"position\"][\"z\"]*15.0)\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\"]*15,node[\"displacement\"][\"y\"]*15,node[\"displacement\"][\"z\"]*15)\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\"]*15.0,node[\"position\"][\"y\"]*15.0,node[\"position\"][\"z\"]*15.0)\n",
" #voxelMaterial(E,mass,nu,rho,momentInertiaInverse,zeta,zetaCollision,muStatic,muKinetic,nomSize)\n",
" i=i+1\n",
" end \n",
"\n",
" i=1\n",
" # pre-calculate the axis\n",
" for 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\"]*15.0 fromNode[\"position\"][\"y\"]*15.0 fromNode[\"position\"][\"z\"]*15.0]\n",
" node2 = [toNode[\"position\"][\"x\"]*15.0 toNode[\"position\"][\"y\"]*15.0 toNode[\"position\"][\"z\"]*15.0]\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 #?????? todo change\n",
" \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",
" #edgeMaterial(E,mass,nu,rho,b,h,L)\n",
"\n",
" i=i+1\n",
" end \n",
" end\n",
" function simulateParallel!(metavoxel,maxNumTimeSteps,dt,returnEvery)\n",
" # initialize(setup)\n",
"\n",
" for i in 1:maxNumTimeSteps\n",
" doTimeStep!(metavoxel,dt,i)\n",
" if(mod(i,saveEvery)==0)\n",
" #append!(displacements,[Array(metavoxel[\"N_displacementGPU\"])])\n",
" updateDataAndSave!(metavoxel,setup,\"../json/live/$(numFile).json\",displacements)\n",
" numFile+=1\n",
" if numFile>maxNumFiles\n",
" numFile=0\n",
" end\n",
" end\n",
" end\n",
" end\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",
" maxNumEdges=10\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_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_material=fill(voxelMaterial(),voxCount)\n",
" #voxelMaterial(E,mass,nu,rho,momentInertiaInverse,zeta,zetaCollision,muStatic,muKinetic,nomSize)\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",
" E_damp=fill(false,linkCount)\n",
" E_material=fill(edgeMaterial(),linkCount)\n",
" E_currentTransverseStrainSum=fill(0.0f0,linkCount)# TODO remove ot incorporate\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_materialGPU= CuArray(N_material) \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_dampGPU= CuArray(E_damp)\n",
" E_materialGPU= CuArray(E_material)\n",
"\n",
"\n",
" #########################################\n",
" metavoxel = Dict(\n",
" \"N_positionGPU\" => N_positionGPU, \n",
" \"N_restrainedGPU\" => N_restrainedGPU, \n",
" \"N_displacementGPU\" => N_displacementGPU,\n",
" \"N_angleGPU\" => N_angleGPU, \n",
" \"N_currPositionGPU\" => N_currPositionGPU,\n",
" \"N_linMomGPU\" => N_linMomGPU, \n",
" \"N_angMomGPU\" => N_angMomGPU, \n",
" \"N_intForceGPU\" => N_intForceGPU, \n",
" \"N_intMomentGPU\" => N_intMomentGPU, \n",
" \"N_momentGPU\" => N_momentGPU, \n",
" \"N_forceGPU\" => N_forceGPU, \n",
" \"N_orientGPU\" => N_orientGPU, \n",
" \"N_edgeIDGPU\" => N_edgeIDGPU, \n",
" \"N_edgeFirstGPU\" => N_edgeFirstGPU,\n",
" \"N_materialGPU\"=> N_materialGPU,\n",
" \"E_sourceGPU\" =>E_sourceGPU, \n",
" \"E_targetGPU\" =>E_targetGPU, \n",
" \"E_areaGPU\" =>E_areaGPU, \n",
" \"E_densityGPU\" =>E_densityGPU, \n",
" \"E_stiffnessGPU\" =>E_stiffnessGPU, \n",
" \"E_stressGPU\" =>E_stressGPU, \n",
" \"E_axisGPU\" =>E_axisGPU, \n",
" \"E_currentRestLengthGPU\" =>E_currentRestLengthGPU, \n",
" \"E_pos2GPU\" =>E_pos2GPU, \n",
" \"E_angle1vGPU\" =>E_angle1vGPU, \n",
" \"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",
" \"E_materialGPU\" =>E_materialGPU\n",
" )\n",
"\n",
" #########################################\n",
" \n",
"\n",
" dt=0.0251646\n",
" E = 2000 # MPa\n",
" s=2.38\n",
" mass=10 \n",
" MaxFreq2=E*s/mass\n",
" dt= 1/(6.283185*sqrt(MaxFreq2))\n",
" println(\"dt: $dt\")\n",
" \n",
" numFile=0\n",
" numSaves=0\n",
" \n",
" \n",
" t=@timed doTimeStep!(metavoxel,dt,0)\n",
" time=t[2]\n",
" println(\"first timestep took $time seconds\")\n",
" t=@timed simulateParallel!(metavoxel,maxNumTimeSteps-1,dt,returnEvery)\n",
" time=t[2]\n",
" \n",
" setup[\"maxNumFiles\"]=numSaves\n",
" \n",
"\n",
" println(\"ran $voxCount voxels and $linkCount edges for $maxNumTimeSteps time steps took $time seconds\")\n",
" return\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"getYoungsModulus (generic function with 1 method)"
]
},
"execution_count": 34,
"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": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"getSetup (generic function with 1 method)"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function getSetup(latticeSize)\n",
" setup = Dict()\n",
" name=string(\"../json/setupTestUni$latticeSize\",\".json\")\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",
" setup[\"viz\"][\"colorMaps\"]=\"\"\n",
" return setup\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 36,
"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": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"DDisplacements=[[],[],[],[],[]]\n",
"Es=[0.0,0,0,0,0]"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"dt: 0.007294855212986816\n",
"first timestep took 0.406329599 seconds\n",
"ran latticeSize 2 with 54 voxels and 144 edges for 10000 time steps took 4.5745323 seconds\n",
"num displacements 10001\n",
"converged displacement= -127.03975744706469\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=\"clip7700\">\n",
" <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip7700)\" 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=\"clip7701\">\n",
" <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip7700)\" 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=\"clip7702\">\n",
" <rect x=\"270\" y=\"121\" width=\"2083\" height=\"1305\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 329.358,1425.62 329.358,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 820.427,1425.62 820.427,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1311.49,1425.62 1311.49,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1802.56,1425.62 1802.56,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 2293.63,1425.62 2293.63,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,1218.47 2352.76,1218.47 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,953.496 2352.76,953.496 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,688.524 2352.76,688.524 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,423.551 2352.76,423.551 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7702)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,158.579 2352.76,158.579 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" 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(#clip7700)\" 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(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 329.358,1425.62 329.358,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 820.427,1425.62 820.427,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1311.49,1425.62 1311.49,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1802.56,1425.62 1802.56,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2293.63,1425.62 2293.63,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,1218.47 295.612,1218.47 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,953.496 295.612,953.496 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,688.524 295.612,688.524 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,423.551 295.612,423.551 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,158.579 295.612,158.579 \n",
" \"/>\n",
"<g clip-path=\"url(#clip7700)\">\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, 329.358, 1479.62)\" x=\"329.358\" y=\"1479.62\">0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 820.427, 1479.62)\" x=\"820.427\" y=\"1479.62\">2500</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 1311.49, 1479.62)\" x=\"1311.49\" y=\"1479.62\">5000</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 1802.56, 1479.62)\" x=\"1802.56\" y=\"1479.62\">7500</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 2293.63, 1479.62)\" x=\"2293.63\" y=\"1479.62\">10000</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 1235.97)\" x=\"246.627\" y=\"1235.97\">-10.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 970.996)\" x=\"246.627\" y=\"970.996\">-7.5</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 706.024)\" x=\"246.627\" y=\"706.024\">-5.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 441.051)\" x=\"246.627\" y=\"441.051\">-2.5</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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, 176.079)\" x=\"246.627\" y=\"176.079\">0.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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\">2 Voxel Convergence Drop</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7700)\">\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(#clip7700)\">\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(#clip7702)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 329.555,158.579 349.197,177.201 368.84,232.696 388.483,325.066 408.126,454.311 427.768,620.429 447.411,823.422 467.054,1062.63 486.697,1269.21 506.339,1388.71 \n",
" 525.982,1379.31 545.625,1256.58 565.267,1086.2 584.91,934.564 604.553,845.784 624.196,833.111 643.838,884.38 663.481,974.262 683.124,1078.4 702.767,1158.18 \n",
" 722.409,1196.34 742.052,1182.59 761.695,1127.42 781.337,1057.51 800.98,998.91 820.623,968.706 840.266,971.808 859.908,1002.62 879.551,1047.14 899.194,1088.58 \n",
" 918.837,1113.1 938.479,1114.69 958.122,1096.14 977.765,1067.37 997.408,1040.34 1017.05,1023.58 1036.69,1020.81 1056.34,1030.73 1075.98,1047.98 1095.62,1065.21 \n",
" 1115.26,1076.48 1134.91,1078.82 1154.55,1072.86 1174.19,1062.05 1193.83,1051.22 1213.48,1044.59 1233.12,1043.9 1252.76,1048.21 1272.41,1054.96 1292.05,1061.32 \n",
" 1311.69,1064.9 1331.33,1064.69 1350.98,1061.33 1370.62,1056.66 1390.26,1052.73 1409.9,1050.95 1429.55,1051.66 1449.19,1054.16 1468.83,1057.19 1488.48,1059.47 \n",
" 1508.12,1060.22 1527.76,1059.38 1547.4,1057.53 1567.05,1055.54 1586.69,1054.19 1606.33,1053.9 1625.97,1054.59 1645.62,1055.83 1665.26,1057.04 1684.9,1057.74 \n",
" 1704.55,1057.74 1724.19,1057.14 1743.83,1056.26 1763.47,1055.48 1783.12,1055.08 1802.76,1055.16 1822.4,1055.59 1842.04,1056.15 1861.69,1056.6 1881.33,1056.77 \n",
" 1900.97,1056.64 1920.62,1056.31 1940.26,1055.93 1959.9,1055.65 1979.54,1055.57 1999.19,1055.67 2018.83,1055.9 2038.47,1056.13 2058.11,1056.27 2077.76,1056.29 \n",
" 2097.4,1056.2 2117.04,1056.05 2136.69,1055.92 2156.33,1055.85 2175.97,1055.87 2195.61,1055.96 2215.26,1056.08 2234.9,1056.19 2254.54,1056.25 2274.18,1056.26 \n",
" 2293.83,1056.23 \n",
" \"/>\n",
"<path clip-path=\"url(#clip7700)\" d=\"\n",
"M1853.56 326.155 L2280.76 326.155 L2280.76 205.195 L1853.56 205.195 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<polyline clip-path=\"url(#clip7700)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1853.56,326.155 2280.76,326.155 2280.76,205.195 1853.56,205.195 1853.56,326.155 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7700)\" 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(#clip7700)\">\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",
"</svg>\n"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"latticeSize =2\n",
"setup=getSetup(latticeSize)\n",
"numTimeSteps=10000\n",
"displacements=[]\n",
"save=true\n",
"returnEvery=1\n",
"runMetavoxelGPU!(setup,numTimeSteps,latticeSize,displacements,returnEvery,true)\n",
"mmm=length(displacements)\n",
"println(\"num displacements $mmm\")\n",
"numTimeStepsRecorded=length(displacements)\n",
"d=[]\n",
"dFEA=[]\n",
"j=length(displacements[end])\n",
"step=100\n",
"for i in 1:step:numTimeStepsRecorded\n",
" append!(d,displacements[i][j].y/15)\n",
"end\n",
"DDisplacements[latticeSize]=d\n",
"\n",
"# E2=getYoungsModulus(latticeSize,75,displacements[end],Load,topNodesIndices)\n",
"\n",
"# Es[latticeSize]=E2\n",
"\n",
"\n",
"\n",
"println(\"converged displacement= $(displacements[numTimeStepsRecorded][j].y)\")\n",
"plot(1:step:numTimeStepsRecorded,d,label=\"Dynamic\",xlabel=\"timestep\",ylabel=\"displacement\",title=\"$latticeSize Voxel Convergence Drop\")\n",
"# savefig(\"$(latticeSize)_voxel_convergence_drop\")\n"
]
},
{
"cell_type": "code",
"execution_count": 21,
"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=\"clip7300\">\n",
" <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip7300)\" 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=\"clip7301\">\n",
" <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip7300)\" d=\"\n",
"M257.245 1425.62 L2352.76 1425.62 L2352.76 121.675 L257.245 121.675 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<defs>\n",
" <clipPath id=\"clip7302\">\n",
" <rect x=\"257\" y=\"121\" width=\"2097\" height=\"1305\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 279.252,1425.62 279.252,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 652.252,1425.62 652.252,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1025.25,1425.62 1025.25,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1398.25,1425.62 1398.25,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1771.25,1425.62 1771.25,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 2144.25,1425.62 2144.25,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 257.245,1219.47 2352.76,1219.47 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 257.245,981.793 2352.76,981.793 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 257.245,744.112 2352.76,744.112 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 257.245,506.432 2352.76,506.432 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7302)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 257.245,268.752 2352.76,268.752 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,1425.62 2352.76,1425.62 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,1425.62 257.245,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 279.252,1425.62 279.252,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 652.252,1425.62 652.252,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1025.25,1425.62 1025.25,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1398.25,1425.62 1398.25,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1771.25,1425.62 1771.25,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2144.25,1425.62 2144.25,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,1219.47 282.391,1219.47 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,981.793 282.391,981.793 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,744.112 282.391,744.112 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,506.432 282.391,506.432 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.245,268.752 282.391,268.752 \n",
" \"/>\n",
"<g clip-path=\"url(#clip7300)\">\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, 279.252, 1479.62)\" x=\"279.252\" y=\"1479.62\">0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 652.252, 1479.62)\" x=\"652.252\" y=\"1479.62\">10</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 1025.25, 1479.62)\" x=\"1025.25\" y=\"1479.62\">20</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 1398.25, 1479.62)\" x=\"1398.25\" y=\"1479.62\">30</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 1771.25, 1479.62)\" x=\"1771.25\" y=\"1479.62\">40</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 2144.25, 1479.62)\" x=\"2144.25\" y=\"1479.62\">50</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 233.245, 1236.97)\" x=\"233.245\" y=\"1236.97\">-126</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 233.245, 999.293)\" x=\"233.245\" y=\"999.293\">-123</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 233.245, 761.612)\" x=\"233.245\" y=\"761.612\">-120</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 233.245, 523.932)\" x=\"233.245\" y=\"523.932\">-117</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 233.245, 286.252)\" x=\"233.245\" y=\"286.252\">-114</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 1305, 73.2)\" x=\"1305\" y=\"73.2\">Node Displacement</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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, 1305, 1559.48)\" x=\"1305\" y=\"1559.48\">Node ID</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip7300)\">\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",
"<circle clip-path=\"url(#clip7302)\" cx=\"316.552\" cy=\"158.664\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"353.852\" cy=\"158.703\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"391.152\" cy=\"158.603\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"428.452\" cy=\"158.724\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"465.752\" cy=\"918.136\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"503.052\" cy=\"546.831\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"540.352\" cy=\"531.095\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"577.652\" cy=\"877.189\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"614.952\" cy=\"905.468\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"652.252\" cy=\"909.57\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"689.552\" cy=\"588.292\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"726.852\" cy=\"538.716\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"764.152\" cy=\"158.689\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"801.451\" cy=\"158.579\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"838.751\" cy=\"158.602\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"876.051\" cy=\"911.192\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"913.351\" cy=\"539.306\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"950.651\" cy=\"893.175\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"987.951\" cy=\"869.282\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1025.25\" cy=\"579.163\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1062.55\" cy=\"1380.66\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1099.85\" cy=\"1154.08\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1137.15\" cy=\"1198.63\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1174.45\" cy=\"1298.83\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1211.75\" cy=\"1307.78\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1249.05\" cy=\"1364.73\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1286.35\" cy=\"1035.18\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1323.65\" cy=\"1129.77\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1360.95\" cy=\"1367.02\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1398.25\" cy=\"1134.53\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1435.55\" cy=\"1281.64\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1472.85\" cy=\"1293.13\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1510.15\" cy=\"1016.58\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1547.45\" cy=\"158.698\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1584.75\" cy=\"158.67\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1622.05\" cy=\"158.7\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1659.35\" cy=\"922.286\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1696.65\" cy=\"603.784\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1733.95\" cy=\"550.292\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1771.25\" cy=\"894.974\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1808.55\" cy=\"919.93\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1845.85\" cy=\"158.661\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1883.15\" cy=\"158.628\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1920.45\" cy=\"890.765\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1957.75\" cy=\"595.307\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"1995.05\" cy=\"907.419\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2032.35\" cy=\"1329.44\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2069.65\" cy=\"1064.94\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2106.95\" cy=\"1163.5\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2144.25\" cy=\"1333.81\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2181.55\" cy=\"1388.71\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2218.85\" cy=\"1328.66\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2256.15\" cy=\"1046.24\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip7302)\" cx=\"2293.45\" cy=\"1301.85\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<path clip-path=\"url(#clip7300)\" d=\"\n",
"M1853.56 326.155 L2280.76 326.155 L2280.76 205.195 L1853.56 205.195 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<polyline clip-path=\"url(#clip7300)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1853.56,326.155 2280.76,326.155 2280.76,205.195 1853.56,205.195 1853.56,326.155 \n",
" \"/>\n",
"<circle clip-path=\"url(#clip7300)\" cx=\"1961.56\" cy=\"265.675\" r=\"21\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<g clip-path=\"url(#clip7300)\">\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\">Dyanmic</text>\n",
"</g>\n",
"</svg>\n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n=[]\n",
"nFEA=[]\n",
"j=length(displacements[end])\n",
"for i in 1:j\n",
" append!(n,displacements[end][i].y)\n",
"end\n",
"scatter(1:j,n,label=\"Dyanmic\",xlabel=\"Node ID\",ylabel=\"displacement\",title=\"Node Displacement\")\n",
"# savefig(\"node_displacement_four_voxel\")"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"dt: 0.007294855212986816\n",
"first timestep took 0.0005791 seconds\n",
"ran 144 voxels and 432 edges for 10000 time steps took 6.5753683 seconds\n"
]
},
{
"data": {
"text/plain": [
"199.0"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"latticeSize=3\n",
"setup=getSetup(latticeSize)\n",
"maxNumTimeSteps=10000\n",
"maxNumFiles=200\n",
"\n",
"saveEvery=round(maxNumTimeSteps/maxNumFiles)\n",
"maxNumFiles=round(maxNumTimeSteps/saveEvery)-1\n",
"setup[\"maxNumFiles\"]=maxNumFiles\n",
"\n",
"runMetavoxelGPULive!(setup,maxNumTimeSteps,saveEvery,maxNumFiles)\n",
"maxNumFiles"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.2.0",
"language": "julia",
"name": "julia-1.2"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}