{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra\n",
"using Plots\n",
"import JSON\n",
"# using Quaternions\n",
"using StaticArrays, Rotations\n",
"using Distributed\n",
"using StaticArrays, BenchmarkTools\n",
"using Base.Threads\n",
"using CUDAnative\n",
"using CuArrays,CUDAdrv \n",
"using Test\n",
"import Base: +, * , -, ^"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Todo\n",
"- create struct for material and get for each edge its properties"
]
},
{
"cell_type": "code",
"execution_count": 79,
"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": 80,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"simulateParallel (generic function with 2 methods)"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function simulateParallel(numTimeSteps,dt)\n",
" # initialize(setup)\n",
" \n",
" for i in 1:numTimeSteps\n",
" #println(\"Timestep:\",i)\n",
" doTimeStep(dt,i)\n",
" end\n",
"end\n",
"\n",
"function simulateParallel(metavoxel,numTimeSteps,dt)\n",
" # initialize(setup)\n",
" \n",
" for i in 1:numTimeSteps\n",
" #println(\"Timestep:\",i)\n",
" doTimeStep(metavoxel,dt,i)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"initialize (generic function with 1 method)"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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",
" 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",
"# E_currentRestLength[i]=75/sqrt(2)\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"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"doTimeStep! (generic function with 1 method)"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function doTimeStep(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!(\n",
" E_sourceGPU, \n",
" E_targetGPU,\n",
" E_areaGPU,\n",
" E_densityGPU,\n",
" E_stiffnessGPU,\n",
" E_stressGPU,\n",
" E_axisGPU,\n",
" E_currentRestLengthGPU,\n",
" E_pos2GPU,\n",
" E_angle1vGPU,\n",
" E_angle2vGPU,\n",
" E_angle1GPU,\n",
" E_angle2GPU,\n",
" E_intForce1GPU,\n",
" E_intMoment1GPU,\n",
" E_intForce2GPU,\n",
" E_intMoment2GPU,\n",
" E_dampGPU,\n",
" N_currPositionGPU,\n",
" N_orientGPU)\n",
" \n",
" # update forces: go through nodes and update interior force (according to int forces from edges), integrate and update currpos\n",
" run_updateNodes!(dt,currentTimeStep,\n",
" N_positionGPU, \n",
" N_restrainedGPU,\n",
" N_displacementGPU,\n",
" N_angleGPU,\n",
" N_currPositionGPU,\n",
" N_linMomGPU,\n",
" N_angMomGPU,\n",
" N_intForceGPU,\n",
" N_intMomentGPU,\n",
" N_forceGPU,\n",
" N_momentGPU,\n",
" N_orientGPU,\n",
" N_edgeIDGPU, \n",
" N_edgeFirstGPU, \n",
" E_intForce1GPU,\n",
" E_intMoment1GPU,\n",
" E_intForce2GPU,\n",
" E_intMoment2GPU)\n",
" \n",
"end\n",
"\n",
"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!(\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[\"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[\"E_intForce1GPU\"],\n",
" metavoxel[\"E_intMoment1GPU\"],\n",
" metavoxel[\"E_intForce2GPU\"],\n",
" metavoxel[\"E_intMoment2GPU\"])\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"run_updateEdges! (generic function with 1 method)"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function updateEdges!(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,\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",
" nu=0\n",
"# L = 5.0 #?? change!!!!!!\n",
" L=l\n",
" a1 = E*L # EA/L : Units of N/m\n",
" a2 = E * L*L*L / (12.0*(1+nu)) # GJ/L : Units of N-m\n",
" b1 = E*L # 12EI/L^3 : Units of N/m\n",
" b2 = E*L*L/2.0 # 6EI/L^2 : Units of N (or N-m/m: torque related to linear distance)\n",
" b3 = E*L*L*L/6.0 # 2EI/L : Units of N-m\n",
" \n",
" nu=0.35\n",
" W = 75\n",
"# L = W/sqrt(2)\n",
" l=L\n",
" n_min = 1\n",
" n_max = 7\n",
" # Cross Section inputs, must be floats\n",
" mass=1\n",
" E =2000 # MPa\n",
" G = E * 1 / 3 # MPa\n",
" h = 2.38 # mm\n",
" b = 2.38 # mm\n",
" rho = 7.85e-9 / 3 # kg/mm^3\n",
" S = h * b\n",
" Sy = (S * (6 + 12 * nu + 6 * nu^2)/ (7 + 12 * nu + 4 * nu^2))\n",
" # For solid rectangular cross section (width=b, depth=d & ( b < d )):\n",
" Q = 1 / 3 - 0.2244 / (min(h / b, b / h) + 0.1607)\n",
" Jxx = Q * min(h * b^3, b * h^3)\n",
" s=b\n",
" \n",
" MaxFreq2=E*s/mass\n",
" dt= 1/(6.283185*sqrt(MaxFreq2))\n",
"\n",
"\n",
" ##if voxels\n",
" #nu=0\n",
" #L=l\n",
" #a1 = E*L # EA/L : Units of N/m\n",
" #a2 = E * L*L*L / (12.0*(1+nu)) # GJ/L : Units of N-m\n",
" #b1 = E*L # 12EI/L^3 : Units of N/m\n",
" #b2 = E*L*L/2.0 # 6EI/L^2 : Units of N (or N-m/m: torque related to linear distance)\n",
" #b3 = E*L*L*L/6.0 # 2EI/L : Units of N-m\n",
"\n",
" I= b*h^3/12\n",
" J=b*h*(b*b+h*h)/12\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",
" \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",
" sqA1=CUDAnative.sqrt(a1) \n",
" sqA2xIp=CUDAnative.sqrt(a2*L*L/6.0) \n",
" sqB1=CUDAnative.sqrt(b1) \n",
" sqB2xFMp=CUDAnative.sqrt(b2*L/2) \n",
" sqB3xIp=CUDAnative.sqrt(b3*L*L/6.0)\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= 2*sqrt(mass)*zeta/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!(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,N_currPosition,N_orient)\n",
" N=length(E_source)\n",
" numblocks = ceil(Int, N/256)\n",
" CuArrays.@sync begin\n",
" @cuda threads=256 blocks=numblocks updateEdges!(E_source,E_target,E_area,E_density,\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,N_currPosition,N_orient)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"run_updateNodes! (generic function with 1 method)"
]
},
"execution_count": 84,
"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_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\n",
"\n",
" index = (blockIdx().x - 1) * blockDim().x + threadIdx().x\n",
" stride = blockDim().x * gridDim().x\n",
" ## @cuprintln(\"thread $index, block $stride\")\n",
" N,M=size(N_edgeID)\n",
" for i = index:stride:N\n",
" @inbounds if N_restrained[i]\n",
" return\n",
" else\n",
" for j in 1:M\n",
" temp=N_edgeID[i,j]\n",
" @inbounds if (N_edgeID[i,j]!=-1)\n",
" #@cuprintln(\"i $i, j $j, N_edgeID[i,j] $temp\")\n",
" @inbounds N_intForce[i]=ifelse(N_edgeFirst[i,j], N_intForce[i]+E_intForce1[N_edgeID[i,j]], N_intForce[i]+E_intForce2[N_edgeID[i,j]] )\n",
" @inbounds N_intMoment[i]=ifelse(N_edgeFirst[i,j], N_intMoment[i]+E_intMoment1[N_edgeID[i,j]], N_intMoment[i]+E_intMoment2[N_edgeID[i,j]] )\n",
" end\n",
" end\n",
" @inbounds curForce = force(N_intForce[i],N_currPosition[i],N_orient[i],N_force[i],true,currentTimeStep)\n",
" \n",
" #@inbounds N_force[i]=Vector3(0,0,0) ##????\n",
" \n",
" @inbounds N_intForce[i]=Vector3(0,0,0)\n",
" \n",
"# x=curForce.x\n",
"# y=curForce.y\n",
"# z=curForce.z\n",
"# @cuprintln(\"curForce x $x, y $y, z $z \")\n",
" \n",
" #x=N_linMom[i].x*10000000000\n",
" #y=N_linMom[i].y*10000000000\n",
" #z=N_linMom[i].z*10000000000\n",
" #@cuprintln(\"N_linMom[i] x $x, y $y, z $z \")\n",
" \n",
" \n",
" @inbounds N_linMom[i]=N_linMom[i]+curForce*Vector3(dt,dt,dt) #todo make sure right\n",
" massInverse=8e-6 # todo ?? later change //8e-6\n",
" mass=1\n",
" massInverse=1.0/mass #\n",
" @inbounds translate=N_linMom[i]*Vector3((dt*massInverse),(dt*massInverse),(dt*massInverse)) # ??massInverse\n",
" \n",
" #x=translate.x*10000000000\n",
" #y=translate.y*10000000000\n",
" #z=translate.z*10000000000\n",
" #@cuprintln(\"translate x $x, y $y, z $z \")\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",
" \n",
" \n",
" \n",
" momentInertiaInverse=1.92e-6 # todo ?? later change 1/Inertia (1/(kg*m^2))\n",
" \n",
" \n",
" @inbounds temp=FromRotationVector(N_angMom[i]*Vector3((dt*momentInertiaInverse),(dt*momentInertiaInverse),(dt*momentInertiaInverse)))\n",
" \n",
" \n",
" #x=temp.x*10000000000\n",
" #y=temp.y*10000000000\n",
" #z=temp.z*10000000000\n",
" #@cuprintln(\"temp x $x, y $y, z $z \")\n",
" \n",
" @inbounds N_orient[i]=multiplyQuaternions(temp,N_orient[i])\n",
" \n",
" #@inbounds x= N_orient[i].x*temp.x\n",
" #@inbounds y= N_orient[i].y*temp.y\n",
" #@inbounds z= N_orient[i].z*temp.z\n",
" #@inbounds w= N_orient[i].w*temp.w\n",
" #x=convert(Float64,x)\n",
" #y=convert(Float64,y)\n",
" #z=convert(Float64,z)\n",
" #w=convert(Float64,w)\n",
" \n",
" #@inbounds N_orient[i]=Quaternion(x,y,z,w)\n",
" \n",
" #x=N_orient[i].x*10000000000\n",
" #y=N_orient[i].y*10000000000\n",
" #z=N_orient[i].z*10000000000\n",
" #w=N_orient[i].w\n",
" #@cuprintln(\"N_orient x $x, y $y, z $z, w $w \")\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_linMom,N_angMom,N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,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_angle,N_currPosition,N_linMom,N_angMom,N_intForce,N_intMoment,N_force,N_moment,N_orient,N_edgeID,N_edgeFirst,E_intForce1,E_intMoment1,E_intForce2,E_intMoment2)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"orientLink! (generic function with 1 method)"
]
},
"execution_count": 85,
"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",
" #x=angle1.x*10000000000\n",
" #y=angle1.y*10000000000\n",
" #z=angle1.z*10000000000\n",
" #@cuprintln(\"angle1 x $x, y $y, z $z \")\n",
" #x=oVNeg.x*10000000000\n",
" #y=oVNeg.y*10000000000\n",
" #z=oVNeg.z*10000000000\n",
" #@cuprintln(\"oVNeg x $x, y $y, z $z \")\n",
" \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",
" #x=pos2.x*10000000000\n",
" #y=pos2.y*10000000000\n",
" #z=pos2.z*10000000000\n",
" #@cuprintln(\"pos2 2 x $x, y $y, z $z \")\n",
" \n",
" \n",
" #x=totalRot.x*10000000000\n",
" #y=totalRot.y*10000000000\n",
" #z=totalRot.z*10000000000\n",
" #@cuprintln(\"totalRot x $x, y $y, z $z \")\n",
" \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",
" 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",
" \n",
" \n",
"# pos2,angle1v,angle2v,angle1,angle2,\n",
" return pos2,angle1v,angle2v,angle1,angle2,totalRot\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"toAxisOriginalQuat (generic function with 1 method)"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function toAxisXVector3(pV::Vector3,axis::Vector3) #TODO CHANGE\n",
"\n",
" xaxis=Vector3(1.0,0.0,0.0)\n",
"\n",
" vector=normalizeVector3(axis)\n",
" q=setFromUnitVectors(vector,xaxis)\n",
" \n",
" \n",
" \n",
"# rot=setFromRotationMatrix(quatToMatrix( q ))\n",
" \n",
" \n",
"# v= applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" #d=15\n",
" #vx=round(v.x, digits=d)\n",
" #vy=round(v.y, digits=d)\n",
" #vz=round(v.z, digits=d)\n",
" \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",
"\n",
"# rot=setFromRotationMatrix(quatToMatrix( q ))\n",
"\n",
"# v= applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" #d=15\n",
" #vx=round(v.x, digits=d)\n",
" #vy=round(v.y, digits=d)\n",
" #vz=round(v.z, digits=d)\n",
" \n",
" \n",
" return Vector3(v.x,v.y,v.z)\n",
"end\n",
"\n",
"function toAxisXQuat(pQ::Quaternion,axis::Vector3)\n",
" \n",
" xaxis=Vector3(1.0,0.0,0.0)\n",
"\n",
" vector=normalizeVector3(axis)\n",
"\n",
"\n",
" q=setFromUnitVectors(vector,xaxis)\n",
" \n",
" #d=17\n",
" #qw=round(q[1], digits=d)\n",
" #qx=round(q[2], digits=d)\n",
" #qy=round(q[3], digits=d)\n",
" #qz=round(q[4], digits=d)\n",
" #\n",
"\n",
"# rot=setFromRotationMatrix(quatToMatrix( q ))\n",
" \n",
" pV=Vector3(pQ.x,pQ.y,pQ.z)\n",
" \n",
"# v=applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" #d=15\n",
" #vx=round(v.x, digits=d)\n",
" #vy=round(v.y, digits=d)\n",
" #vz=round(v.z, digits=d)\n",
" \n",
" return Quaternion(v.x,v.y,v.z,1.0)\n",
" \n",
" # return [1.0 v[1] v[2] v[3]]\n",
"end\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",
"\n",
"# rot=setFromRotationMatrix(quatToMatrix( q ))\n",
" \n",
" pV=Vector3(pQ.x,pQ.y,pQ.z)\n",
"# v=applyQuaternion1( pV ,setQuaternionFromEuler(rot) )\n",
" v=applyQuaternion1( pV ,q )\n",
" \n",
" #d=15\n",
" #vx=round(v.x, digits=d)\n",
" #vy=round(v.y, digits=d)\n",
" #vz=round(v.z, digits=d)\n",
" \n",
" return Quaternion(v.x,v.y,v.z,1.0)\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"applyQuaternion1 (generic function with 1 method)"
]
},
"execution_count": 87,
"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",
" #te = m\n",
" #m11 = (te[ 1 ]== -0.0) ? 0.0 : te[ 1 ]\n",
" #m12 = (te[ 5 ]== -0.0) ? 0.0 : te[ 5 ]\n",
" #m13 = (te[ 9 ]== -0.0) ? 0.0 : te[ 9 ]\n",
" #m21 = (te[ 2 ]== -0.0) ? 0.0 : te[ 2 ]\n",
" #m22 = (te[ 6 ]== -0.0) ? 0.0 : te[ 6 ]\n",
" #m23 = (te[ 10]== -0.0) ? 0.0 : te[ 10]\n",
" #m31 = (te[ 3 ]== -0.0) ? 0.0 : te[ 3 ]\n",
" #m32 = (te[ 7 ]== -0.0) ? 0.0 : te[ 7 ]\n",
" #m33 = (te[ 11]== -0.0) ? 0.0 : te[ 11]\n",
"\n",
" m11 = 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",
" # if(m23==0.0)\n",
" # x = atan( m23, m33 )\n",
" # end\n",
" # if(m12==0.0)\n",
" # z = atan( m12, m11 )\n",
" # end\n",
"\n",
" else\n",
"\n",
" x = 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": 88,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"multiplyQuaternions (generic function with 1 method)"
]
},
"execution_count": 88,
"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",
"# x1=convert(Float64,x1)\n",
"# y1=convert(Float64,y1)\n",
"# z1=convert(Float64,z1)\n",
"# w1=convert(Float64,w1)\n",
"\treturn Quaternion(x1,y1,z1,w1 ); #!< overload quaternion multiplication.\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": 89,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"axialStrain (generic function with 1 method)"
]
},
"execution_count": 89,
"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": 90,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"moment (generic function with 1 method)"
]
},
"execution_count": 90,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function force(N_intForce,N_currPosition,N_orient,N_force,static,currentTimeStep) \n",
" # forces from internal bonds\n",
" totalForce=Vector3(0,0,0)\n",
"\n",
" \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",
" static=false\n",
" # other forces\n",
" if(static)\n",
" totalForce=totalForce+N_force\n",
"\n",
" else\n",
" if convert(Float64,N_force.z)==-1.0\n",
" \n",
" if(currentTimeStep<10000.0)# ||(currentTimeStep>20000&¤tTimeStep<30000))\n",
"# ax=normalizeVector3(N_currPosition)\n",
"# totalForce+=Vector3(-0.1*ax.x,-0.1*ax.y,-0.1*ax.z)\n",
" totalForce+=Vector3(0.0,0.0,-1.0*CUDAnative.sin(currentTimeStep/10000.0*3.14))\n",
" end\n",
" elseif convert(Float64,N_force.z)==-2.0\n",
" if(currentTimeStep>30000.0||(currentTimeStep>10000&¤tTimeStep<20000))\n",
"# ax=normalizeVector3(N_currPosition)\n",
"# totalForce+=Vector3(-0.1*ax.x,-0.1*ax.y,-0.1*ax.z)\n",
"# totalForce+=Vector3(0.0,0.0,-0.1)\n",
" end\n",
" end\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": 91,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"updateDataAndSaveFEA! (generic function with 1 method)"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function updateDataAndSaveFEA!(setup,fileName)\n",
"# nodes = setup[\"nodes\"]\n",
"# edges = setup[\"edges\"]\n",
" \n",
"# setup[\"animation\"][\"showDisplacement\"]=false\n",
"# voxCount=size(nodes)[1]\n",
"# linkCount=size(edges)[1]\n",
" \n",
"# N_displacement=Array(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",
"\n",
"\n",
"# i=1\n",
"# \tfor edge in edges\n",
"# edge[\"stress\"]=E_stress[i]\n",
"# i=i+1\n",
"\n",
"# end\n",
" \n",
" \n",
"# i=1 \n",
"# \tfor node in nodes\n",
"# node[\"displacement\"][\"x\"]=N_displacement[i].x*100\n",
"# node[\"displacement\"][\"y\"]=N_displacement[i].y*100\n",
"# node[\"displacement\"][\"z\"]=N_displacement[i].z*100\n",
" \n",
"# node[\"angle\"][\"x\"]=N_angle[i].x\n",
"# node[\"angle\"][\"y\"]=N_angle[i].y\n",
"# node[\"angle\"][\"z\"]=N_angle[i].z\n",
"# i=i+1\n",
"\n",
"# end\n",
" \n",
" # pass data as a json string (how it shall be displayed in a file)\n",
" stringdata = JSON.json(setup)\n",
"\n",
" # write the file with the stringdata variable information\n",
" open(fileName, \"w\") do f\n",
" write(f, stringdata)\n",
" end\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"updateDataAndSave! (generic function with 2 methods)"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function updateDataAndSave!(metavoxel,setup,fileName)\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\"]=20.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[\"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",
" node[\"angle\"][\"x\"]=N_angle[i].x\n",
" node[\"angle\"][\"y\"]=N_angle[i].y\n",
" node[\"angle\"][\"z\"]=N_angle[i].z\n",
" i=i+1\n",
"\n",
" end\n",
" \n",
" # pass data as a json string (how it shall be displayed in a file)\n",
" stringdata = JSON.json(setup)\n",
"\n",
" # write the file with the stringdata variable information\n",
" open(fileName, \"w\") do f\n",
" write(f, stringdata)\n",
" end\n",
" \n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"updateDataAndSave! (generic function with 2 methods)"
]
},
"execution_count": 93,
"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\"]=100.0\n",
" setup[\"static\"]=false\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",
" \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": 94,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"runMetavoxelGPU! (generic function with 1 method)"
]
},
"execution_count": 94,
"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",
" # 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",
"# E_currentRestLength[i]=75/sqrt(2)\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",
"\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_posTimeSteps=[]\n",
" # N_angTimeSteps=[]\n",
" N_force=fill(Vector3(),voxCount)\n",
" N_orient=fill(Quaternion(),voxCount)\n",
" N_edgeID=fill(-1,(voxCount,maxNumEdges))\n",
" N_edgeFirst=fill(true,(voxCount,maxNumEdges))\n",
" N_currEdge=fill(1,voxCount)\n",
"\n",
" ############# edges\n",
" E_source=fill(0,linkCount)\n",
" E_target=fill(0,linkCount)\n",
" E_area=fill(0.0f0,linkCount)\n",
" E_density=fill(0.0f0,linkCount)\n",
" E_stiffness=fill(0.0f0,linkCount)\n",
" E_stress=fill(0.0f0,linkCount)\n",
" E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)\n",
" E_currentRestLength=fill(0.0f0,linkCount)\n",
" E_pos2=fill(Vector3(),linkCount)\n",
" E_angle1v=fill(Vector3(),linkCount)\n",
" E_angle2v=fill(Vector3(),linkCount)\n",
" E_angle1=fill(Quaternion(),linkCount)\n",
" E_angle2=fill(Quaternion(),linkCount)\n",
"\n",
" E_intForce1=fill(Vector3(),linkCount)\n",
" E_intMoment1=fill(Vector3(),linkCount) \n",
"\n",
" E_intForce2=fill(Vector3(),linkCount)\n",
" E_intMoment2=fill(Vector3(),linkCount)\n",
" E_damp=fill(false,linkCount)\n",
"\n",
" E_currentTransverseStrainSum=fill(0.0f0,linkCount)# TODO remove ot incorporate\n",
" # E_stressTimeSteps=[]\n",
"\n",
"\n",
" #################################################################\n",
" initialize!(setup)\n",
" #################################################################\n",
"\n",
" ########################## turn to cuda arrays\n",
" ############# nodes\n",
" N_positionGPU= CuArray(N_position) \n",
" N_restrainedGPU= CuArray(N_restrained) \n",
" N_displacementGPU=CuArray(N_displacement) \n",
" N_angleGPU= CuArray(N_angle) \n",
" N_currPositionGPU=CuArray(N_currPosition) \n",
" N_linMomGPU= CuArray(N_linMom) \n",
" N_angMomGPU= CuArray(N_angMom) \n",
" N_intForceGPU= CuArray(N_intForce) \n",
" N_intMomentGPU= CuArray(N_intMoment) \n",
" N_momentGPU= CuArray(N_moment) \n",
" N_forceGPU= CuArray(N_force) \n",
" N_orientGPU= CuArray(N_orient) \n",
" N_edgeIDGPU= CuArray(N_edgeID) \n",
" N_edgeFirstGPU= CuArray(N_edgeFirst) \n",
"\n",
"\n",
" ############# edges\n",
" E_sourceGPU= CuArray(E_source) \n",
" E_targetGPU= CuArray(E_target)\n",
" E_areaGPU= CuArray(E_area) \n",
" E_densityGPU= CuArray(E_density)\n",
" E_stiffnessGPU= CuArray(E_stiffness)\n",
" E_stressGPU= CuArray(E_stress)\n",
" E_axisGPU= CuArray(E_axis) \n",
" E_currentRestLengthGPU= CuArray(E_currentRestLength)\n",
" E_pos2GPU= CuArray(E_pos2)\n",
" E_angle1vGPU= CuArray(E_angle1v)\n",
" E_angle2vGPU= CuArray(E_angle2v)\n",
" E_angle1GPU= CuArray(E_angle1)\n",
" E_angle2GPU= CuArray(E_angle2)\n",
" E_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)\n",
" E_intForce1GPU= CuArray(E_intForce1) \n",
" E_intMoment1GPU= CuArray(E_intMoment1) \n",
" E_intForce2GPU= CuArray(E_intForce2) \n",
" E_intMoment2GPU= CuArray(E_intMoment2)\n",
" E_dampGPU= CuArray(E_damp) \n",
" # E_stressTimeSteps=[]\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",
" \"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",
" )\n",
"\n",
" #########################################\n",
" \n",
"\n",
" dt=0.0251646\n",
" E =2000 # MPa\n",
" s=2.38\n",
" mass=1 \n",
" \n",
" \n",
" \n",
" MaxFreq2=E*s/mass\n",
" dt= 1/(6.283185*sqrt(MaxFreq2))\n",
"# dt=0.01\n",
" println(\"dt: $dt\")\n",
" \n",
" append!(displacements,[Array(metavoxel[\"N_displacementGPU\"])])\n",
" \n",
" t=@timed doTimeStep!(metavoxel,dt,0)\n",
" append!(displacements,[Array(metavoxel[\"N_displacementGPU\"])])\n",
" time=t[2]\n",
" println(\"first timestep took $time seconds\")\n",
" t=@timed simulateParallel!(metavoxel,numTimeSteps-1,dt,returnEvery)\n",
" time=t[2]\n",
" \n",
" if save\n",
" updateDataAndSave!(metavoxel,setup,\"../json/trialJuliaParallelGPU.json\")\n",
" updateDataAndSave!(metavoxel,setup,\"../json/trialJuliaParallelGPUDynamic.json\",displacements)\n",
" end\n",
" println(\"ran latticeSize $latticeSize with $voxCount voxels and $linkCount edges for $numTimeSteps time steps took $time seconds\")\n",
" return\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 95,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"fea (generic function with 1 method)"
]
},
"execution_count": 95,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function fea(setup)\n",
" #######################################################\n",
" function points(element, properties)\n",
" elements = properties[\"elements\"]\n",
" nodes = properties[\"nodes\"]\n",
" degrees_of_freedom = properties[\"degrees_of_freedom\"]\n",
"\n",
" # find the nodes that the lements connects\n",
" fromNode = elements[element][1]\n",
" toNode = elements[element][2]\n",
"\n",
" # the coordinates for each node\n",
" fromPoint = nodes[fromNode]\n",
" toPoint = nodes[toNode]\n",
"\n",
" # find the degrees of freedom for each node\n",
" dofs = degrees_of_freedom[fromNode]\n",
" dofs=vcat(dofs,degrees_of_freedom[toNode])\n",
"\n",
" return fromPoint, toPoint, dofs\n",
" end\n",
"\n",
" function direction_cosine(vec1, vec2)\n",
" return dot(vec1,vec2) / (norm(vec1) * norm(vec2))\n",
" end\n",
"\n",
" function rotation_matrix(element_vector, x_axis, y_axis,z_axis)\n",
" # find the direction cosines\n",
" x_proj = direction_cosine(element_vector, x_axis)\n",
" y_proj = direction_cosine(element_vector, y_axis)\n",
" z_proj = direction_cosine(element_vector, z_axis);\n",
" return [[x_proj y_proj z_proj 0 0 0];[0 0 0 x_proj y_proj z_proj]]\n",
" end\n",
"\n",
" function rotation_matrix(element_vector, x_axis, y_axis,z_axis)\n",
" # find the direction cosines\n",
" L=norm(element_vector)\n",
" l = (element_vector[1])/L\n",
" m = (element_vector[2])/L\n",
" n = (element_vector[3])/L\n",
" D = ( l^2+ m^2+n^2)^0.5\n",
"\n",
" transMatrix=[[l m n 0 0 0 0 0 0 0 0 0];[-m/D l/D 0 0 0 0 0 0 0 0 0 0];[ -l*n/D -m*n/D D 0 0 0 0 0 0 0 0 0];[ 0 0 0 l m n 0 0 0 0 0 0];[ 0 0 0 -m/D l/D 0 0 0 0 0 0 0];[ 0 0 0 -l*n/D -m*n/D D 0 0 0 0 0 0];[ 0 0 0 0 0 0 l m n 0 0 0];[ 0 0 0 0 0 0 -m/D l/D 0 0 0 0];[ 0 0 0 0 0 0 -l*n/D -m*n/D D 0 0 0];[ 0 0 0 0 0 0 0 0 0 l m n];[ 0 0 0 0 0 0 0 0 0 -m/D l/D 0];[ 0 0 0 0 0 0 0 0 0 -l*n/D -m*n/D D]]\n",
"\n",
" return transMatrix\n",
" end\n",
" \n",
" #######################################################\n",
" function get_matrices(setup)\n",
"\n",
" nodes = setup[\"nodes\"]\n",
" edges = setup[\"edges\"]\n",
" ndofs = length(nodes)*6\n",
"\n",
" x_axis = [1 0 0]\n",
" y_axis = [0 1 0]\n",
" z_axis = [0 0 1]\n",
"\n",
" M = zeros((ndofs,ndofs))\n",
" K = zeros((ndofs,ndofs))\n",
" \n",
" \n",
" for edge in edges\n",
" #degrees_of_freedom = properties[\"degrees_of_freedom\"]\n",
"\n",
" element=parse(Int,edge[\"id\"][2:end])\n",
"\n",
" # find the nodes that the lements connects\n",
" fromNode = nodes[edge[\"source\"]+1]\n",
" toNode = nodes[edge[\"target\"]+1]\n",
" \n",
"\n",
" # the coordinates for each node\n",
" fromPoint = [fromNode[\"position\"][\"x\"]*15.0 fromNode[\"position\"][\"y\"]*15.0 fromNode[\"position\"][\"z\"]*15.0]\n",
" toPoint = [toNode[\"position\"][\"x\"]*15.0 toNode[\"position\"][\"y\"]*15.0 toNode[\"position\"][\"z\"]*15.0]\n",
" \n",
" fromPoint = fromPoint .+ [fromNode[\"displacement\"][\"x\"]*15.0 fromNode[\"displacement\"][\"y\"]*15.0 fromNode[\"displacement\"][\"z\"]*15.0]\n",
" toPoint = toPoint .+ [toNode[\"displacement\"][\"x\"]*15.0 toNode[\"displacement\"][\"y\"]*15.0 toNode[\"displacement\"][\"z\"]*15.0]\n",
"\n",
" # find the degrees of freedom for each node\n",
" dofs = convert(Array{Int}, fromNode[\"degrees_of_freedom\"])\n",
" dofs=vcat(dofs,convert(Array{Int}, toNode[\"degrees_of_freedom\"]))\n",
"\n",
" element_vector=toPoint-fromPoint\n",
"\n",
" # find element mass and stifness matrices\n",
" length = norm(element_vector)\n",
" rho = edge[\"density\"]\n",
" area = edge[\"area\"]\n",
" E = edge[\"stiffness\"]# youngs modulus\n",
"\n",
" A = edge[\"area\"]\n",
" G=1.0#todo shear_modulus\n",
" ixx = 1.0#todo section ixx\n",
" iyy = 1.0#todo section.iyy#\n",
" l0=length\n",
" j=1.0;#todo check\n",
" l02 = l0 * l0\n",
" l03 = l0 * l0 * l0\n",
" \n",
" # find element mass and stifness matrices\n",
" length = norm(element_vector)\n",
" rho = edge[\"density\"]\n",
" area = edge[\"area\"]\n",
" E = edge[\"stiffness\"]# youngs modulus\n",
" \n",
" \n",
" \n",
" A = edge[\"area\"]\n",
" G=1.0#todo shear_modulus\n",
" ixx = 1.0#todo section ixx\n",
" iyy = 1.0#todo section.iyy#\n",
" j=1.0;#todo check\n",
" \n",
" \n",
" h = 2.38 # mm\n",
" b = 2.38 # mm\n",
" E =2000 # MPa\n",
" rho = 7.85e-9 / 3 # kg/mm^3\n",
" G = E * 1 / 3 # MPa\n",
" A=h*b\n",
" Q = 1 / 3 - 0.2244 / (min(h / b, b / h) + 0.1607)\n",
" J = Q * min(h * b^3, b * h^3)\n",
" I= b*h^3/12\n",
" ixx=I\n",
" iyy=I\n",
" j=J\n",
" \n",
" l0=length\n",
" l02 = l0 * l0\n",
" l03 = l0 * l0 * l0\n",
"\n",
" # Cm = rho * area * length /6.0\n",
" # Ck= E * area / length \n",
"\n",
" # m = [[2 1];[1 2]]\n",
" # k = [[1 -1];[-1 1]]\n",
"\n",
" k = [[E*A/l0 0 0 0 0 0 -E*A/l0 0 0 0 0 0];[0 12*E*ixx/l03 0 0 0 6*E*ixx/l02 0 -12*E*ixx/l03 0 0 0 6*E*ixx/l02];[0 0 12*E*iyy/l03 0 -6*E*iyy/l02 0 0 0 -12*E*iyy/l03 0 -6*E*iyy/l02 0];[0 0 0 G*j/l0 0 0 0 0 0 -G*j/l0 0 0];[0 0 -6*E*iyy/l02 0 4*E*iyy/l0 0 0 0 6*E*iyy/l02 0 2*E*iyy/l0 0];[0 6*E*ixx/l02 0 0 0 4*E*ixx/l0 0 -6*E*ixx/l02 0 0 0 2*E*ixx/l0];[-E*A/l0 0 0 0 0 0 E*A/l0 0 0 0 0 0];[0 -12*E*ixx/l03 0 0 0 -6*E*ixx/l02 0 12*E*ixx/l03 0 0 0 -6*E*ixx/l02];[0 0 -12*E*iyy/l03 0 6*E*iyy/l02 0 0 0 12*E*iyy/l03 0 6*E*iyy/l02 0];[0 0 0 -G*j/l0 0 0 0 0 0 G*j/l0 0 0];[0 0 -6*E*iyy/l02 0 2*E*iyy/l0 0 0 0 6*E*iyy/l02 0 4*E*iyy/l0 0];[0 6*E*ixx/l02 0 0 0 2*E*ixx/l0 0 -6*E*ixx/l02 0 0 0 4*E*ixx/l0]]\n",
"\n",
" \n",
" ################################\n",
"\n",
" # find rotated mass and stifness matrices\n",
" tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis)\n",
"\n",
" # m_r=transpose(tau)*m*tau\n",
" k_r=transpose(tau)*k*tau\n",
"\n",
" # change from element to global coordinate\n",
" index= dofs.+1\n",
"\n",
" B=zeros((12,ndofs))\n",
" for i in 1:12\n",
" B[i,index[i]]=1.0\n",
" end\n",
"\n",
"\n",
" # M_rG= transpose(B)*m_r*B\n",
" K_rG= transpose(B)*k_r*B\n",
"\n",
" # M += Cm .* M_rG\n",
" # K += Ck .* K_rG\n",
" K += K_rG\n",
"\n",
" end\n",
" \n",
" \n",
" # construct the force vector\n",
" F=zeros(ndofs)\n",
" remove_indices=[];\n",
" for node in nodes\n",
" #insert!(F,i, value);\n",
" #F=vcat(F,value)\n",
" \n",
" \n",
" i=parse(Int,node[\"id\"][2:end])\n",
" f=node[\"force\"]\n",
" \n",
" # println(f)\n",
" F[(i)*6+1]+=f[\"x\"]\n",
" F[(i)*6+2]+=f[\"y\"]\n",
" F[(i)*6+3]+=f[\"z\"]\n",
" F[(i)*6+4]+=0\n",
" F[(i)*6+5]+=0\n",
" F[(i)*6+6]+=0\n",
" Load+=f[\"y\"]\n",
" \n",
" if f[\"z\"]==-1.0\n",
" currPos=Vector3((node[\"position\"][\"x\"]+node[\"displacement\"][\"x\"])*15.0,(node[\"position\"][\"y\"]+node[\"displacement\"][\"y\"])*15.0,(node[\"position\"][\"z\"]+node[\"displacement\"][\"z\"])*15.0)\n",
" currPos=Vector3(0.0,(node[\"position\"][\"y\"]+node[\"displacement\"][\"y\"])*15.0,(node[\"position\"][\"z\"]+node[\"displacement\"][\"z\"])*15.0)\n",
"\n",
" ff=normalizeVector3(currPos)\n",
"# F[(i)*6+1]+=f[\"x\"]*-0.1\n",
"# F[(i)*6+2]+=f[\"y\"]*-0.1\n",
"# F[(i)*6+3]+=f[\"z\"]*-0.1\n",
" F[(i)*6+3]=-0.1\n",
" else\n",
" F[(i)*6+3]=0.0\n",
" end\n",
" \n",
" if (F[(i)*6+2]!=0)\n",
" append!(topNodesIndices,i+1)\n",
" end\n",
" \n",
" dofs = convert(Array{Int}, node[\"degrees_of_freedom\"]).+1\n",
" restrained_dofs=node[\"restrained_degrees_of_freedom\"]\n",
" for (index, value) in enumerate(dofs)\n",
" if restrained_dofs[index]\n",
" append!( remove_indices, value)\n",
" end\n",
" end\n",
" \n",
" end\n",
"\n",
" #println(remove_indices)\n",
" #print(K)\n",
" #print(F)\n",
" \n",
"\n",
" #M = M[setdiff(1:end, remove_indices), :]\n",
" K = K[setdiff(1:end, remove_indices), :]\n",
"\n",
" #M = M[:,setdiff(1:end, remove_indices)]\n",
" K = K[:,setdiff(1:end, remove_indices)]\n",
"\n",
" F = F[setdiff(1:end, remove_indices)]\n",
" \n",
" U=zeros(ndofs)\n",
" \n",
" return M,K,F,U,remove_indices\n",
" end\n",
"\n",
" \n",
" function updateDisplacement(setup, X)\n",
" nodes= setup[\"nodes\"]\n",
" i=0\n",
" for node in nodes\n",
" \n",
"# if !node[\"restrained_degrees_of_freedom\"][2]\n",
" #i=parse(Int,node[\"id\"][2:end])\n",
" node[\"displacement\"][\"x\"]+=X[(i)*6+1]/15\n",
" node[\"displacement\"][\"y\"]+=X[(i)*6+2]/15\n",
" node[\"displacement\"][\"z\"]+=X[(i)*6+3]/15\n",
" node[\"angle\"][\"x\"]+=X[(i)*6+4]\n",
" node[\"angle\"][\"y\"]+=X[(i)*6+5]\n",
" node[\"angle\"][\"z\"]+=X[(i)*6+6]\n",
" append!(displacementFEA,[Vector3(node[\"displacement\"][\"x\"]*15,node[\"displacement\"][\"y\"]*15,node[\"displacement\"][\"z\"]*15)])\n",
" i=i+1\n",
"# else\n",
"# append!(displacementFEA,[Vector3(0,0,0)])\n",
"# end\n",
" end\n",
" end\n",
" \n",
" #######################################################\n",
"\n",
" function get_stresses(setup)\n",
" nodes = setup[\"nodes\"]\n",
" edges = setup[\"edges\"]\n",
" ndofs = length(nodes)*6\n",
"\n",
" x_axis = [1 0 0]\n",
" y_axis = [0 1 0]\n",
" z_axis = [0 0 1]\n",
"\n",
" # find the stresses in each member\n",
" stresses=zeros(length(edges))\n",
" max11=-10e6\n",
" min11=10e6\n",
" for edge in edges\n",
" #degrees_of_freedom = properties[\"degrees_of_freedom\"]\n",
"\n",
" element=parse(Int,edge[\"id\"][2:end])\n",
"\n",
" # find the nodes that the lements connects\n",
" fromNode = nodes[edge[\"source\"]+1]\n",
" toNode = nodes[edge[\"target\"]+1]\n",
"\n",
" # the coordinates for each node\n",
" fromPoint = [fromNode[\"position\"][\"x\"]*15.0 fromNode[\"position\"][\"y\"]*15.0 fromNode[\"position\"][\"z\"]*15.0]\n",
" toPoint = [toNode[\"position\"][\"x\"]*15.0 toNode[\"position\"][\"y\"]*15.0 toNode[\"position\"][\"z\"]*15.0]\n",
" \n",
" fromPoint = fromPoint .+ [fromNode[\"displacement\"][\"x\"]*15.0 fromNode[\"displacement\"][\"y\"]*15.0 fromNode[\"displacement\"][\"z\"]*15.0]\n",
" toPoint = toPoint .+ [toNode[\"displacement\"][\"x\"]*15.0 toNode[\"displacement\"][\"y\"]*15.0 toNode[\"displacement\"][\"z\"]*15.0]\n",
"\n",
"\n",
" # find the degrees of freedom for each node\n",
" dofs = convert(Array{Int}, fromNode[\"degrees_of_freedom\"])\n",
" dofs=vcat(dofs,convert(Array{Int}, toNode[\"degrees_of_freedom\"]))\n",
"\n",
" element_vector=toPoint-fromPoint\n",
"\n",
"\n",
" # find rotated mass and stifness matrices\n",
" tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis)\n",
"\n",
" # i1=parse(Int,fromNode[\"id\"][2:end])\n",
" # i2=parse(Int,toNode[\"id\"][2:end])\n",
"\n",
" # global_displacements=[X[(i1)*6+1] X[(i1)*6+2] X[(i1)*6+3] X[(i1)*6+4] X[(i1)*6+5] X[(i1)*6+6] X[(i2)*6+1] X[(i2)*6+2] X[(i2)*6+3] X[(i2)*6+4] X[(i2)*6+5] X[(i2)*6+6]] # todo change\n",
" global_displacements=[fromNode[\"displacement\"][\"x\"]*15 fromNode[\"displacement\"][\"y\"]*15 fromNode[\"displacement\"][\"z\"]*15 fromNode[\"angle\"][\"x\"] fromNode[\"angle\"][\"y\"] fromNode[\"angle\"][\"z\"] toNode[\"displacement\"][\"x\"]*15 toNode[\"displacement\"][\"y\"]*15 toNode[\"displacement\"][\"z\"]*15 toNode[\"angle\"][\"x\"] toNode[\"angle\"][\"y\"] toNode[\"angle\"][\"z\"]] # todo change\n",
"\n",
" # nodal displacement\n",
"\n",
" q=tau*transpose(global_displacements)\n",
" # println(q)\n",
" # calculate the strain and stresses\n",
" strain =(q[7]-q[1])/norm(element_vector)\n",
" E = edge[\"stiffness\"]# youngs modulus\n",
" E =2000\n",
" stress=E.*strain\n",
" edge[\"stress\"]=stress\n",
" if stress>max11\n",
" max11=stress\n",
" end\n",
" if stress<min11\n",
" min11=stress\n",
" end\n",
" # println(element)\n",
" # println(stress)\n",
" end\n",
"\n",
"\n",
"\n",
" setup[\"viz\"][\"minStress\"]=min11\n",
" setup[\"viz\"][\"maxStress\"]=max11\n",
" setup[\"animation\"][\"exaggeration\"]=20.0\n",
" return stresses\n",
" end\n",
" \n",
" function initialize(setup)\n",
" nodes = setup[\"nodes\"]\n",
" ndofs = length(nodes)*6\n",
" \n",
" i=0\n",
" for node in nodes\n",
" dg=[]\n",
" for ii in 0:5\n",
" append!(dg,i+ii) \n",
" end\n",
" i+=6\n",
" node[\"degrees_of_freedom\"]=dg\n",
" end\n",
" end\n",
"\n",
" #######################################################\n",
" function solveFea(setup)\n",
" // # determine the global matrices\n",
" initialize(setup)\n",
" \n",
" M,K,F,U,ind=get_matrices(setup)\n",
" \n",
" #println(M)\n",
" #println(K)\n",
" #println(F)\n",
"\n",
" #evals=eigvals(K,M)\n",
" #evecs=eigvecs(K,M)\n",
" #frequencies=sqrt.(evals)\n",
" X=inv(K)*F\n",
" U[setdiff(1:end, ind)]=X\n",
"\n",
" updateDisplacement(setup, U)\n",
"\n",
" # determine the stresses in each element\n",
" stresses=get_stresses(setup)\n",
" end\n",
" #######################################################\n",
" displacementFEA=[]\n",
" Load=0\n",
" topNodesIndices=[]\n",
" solveFea(setup)\n",
" return displacementFEA,Load,topNodesIndices\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"getYoungsModulus (generic function with 1 method)"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function getYoungsModulus(latticeSize,voxelSize,disp,Load,topNodesIndices)\n",
" F=-Load\n",
" l0=voxelSize*latticeSize\n",
" A=l0*l0\n",
"\n",
" δl1=-mean( x.y for x in disp[topNodesIndices])\n",
" \n",
" stresses=F/A\n",
" strain=δl1/l0\n",
" println(\"Load=$Load\")\n",
" println(\"stress=$stresses\")\n",
"\n",
" E=stresses/strain \n",
"\n",
" return E\n",
"end\n"
]
},
{
"cell_type": "code",
"execution_count": 97,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"getSetup (generic function with 1 method)"
]
},
"execution_count": 97,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function getSetup(latticeSize)\n",
" setup = Dict()\n",
" name=string(\"../json/setupTestUni$latticeSize\",\".json\")\n",
" name=string(\"../json/canteliver\",\".json\")\n",
" name=string(\"../json/snake\",\".json\")\n",
" \n",
"# open(\"../json/setupValid2.json\", \"r\") do f\n",
"# open(\"../json/setupTest.json\", \"r\") do f\n",
" # open(\"../json/trialJulia.json\", \"r\") do f\n",
"# open(\"../json/setupTestUni4.json\", \"r\") do f\n",
" # open(\"../json/setupChiral.json\", \"r\") do f\n",
"# open(\"../json/setupTestCubeUni10.json\", \"r\") do f\n",
" open(name, \"r\") do f\n",
"# global setup\n",
" dicttxt = String(read(f)) # file information to string\n",
" setup=JSON.parse(dicttxt) # parse and transform data\n",
" end\n",
"\n",
" setup=setup[\"setup\"]\n",
" return setup\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 98,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"44\n"
]
},
{
"data": {
"text/plain": [
"104"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"latticeSize=1\n",
"setup=getSetup(latticeSize)\n",
"displacementFEA,Load,topNodesIndices=fea(setup)\n",
"topNodesIndices\n",
"println(length(setup[\"nodes\"]))\n",
"length(setup[\"edges\"])"
]
},
{
"cell_type": "code",
"execution_count": 99,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"5-element Array{Float64,1}:\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"DDisplacements=[[],[],[],[],[]]\n",
"DDisplacementsFEA=[[],[],[],[],[]]\n",
"Loads=[0.0,0,0,0,0]\n",
"EsFEA=[0.0,0,0,0,0]\n",
"Es=[0.0,0,0,0,0]"
]
},
{
"cell_type": "code",
"execution_count": 100,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Dict{String,Any} with 7 entries:\n",
" \"nodes\" => Any[Dict{String,Any}(\"degrees_of_freedom\"=>Any[0, 1, 2, 3, 4,…\n",
" \"voxelSize\" => 5\n",
" \"hierarchical\" => false\n",
" \"animation\" => Dict{String,Any}(\"speed\"=>3,\"exaggeration\"=>20.0,\"showDisplace…\n",
" \"viz\" => Dict{String,Any}(\"colorMap\"=>0,\"colorMaps\"=>Any[Any[Any[0, An…\n",
" \"edges\" => Any[Dict{String,Any}(\"source\"=>0,\"area\"=>5.6644,\"density\"=>0.…\n",
" \"ndofs\" => 264"
]
},
"execution_count": 100,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"setup"
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"44-element Array{Any,1}:\n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.0, y:0.0, z:0.0 \n",
" x:0.010613496971641401, y:-0.01046280778290931, z:-0.01059129790553449 \n",
" â‹® \n",
" x:-7.762493357730787e-5, y:2.9952846496351806e-15, z:-0.06302666774019734 \n",
" x:-0.010957579507126305, y:0.010158333207868608, z:-0.052980102215182984 \n",
" x:-3.500636791186467e-5, y:6.963266467751947e-5, z:-0.0640476197001775 \n",
" x:-2.3096342910205877e-5, y:2.9823993790924887e-15, z:-0.06309994346837483\n",
" x:0.002222557388210866, y:-0.0033974994957394963, z:-0.06642143213658186 \n",
" x:-0.017667638965804384, y:-0.018868777070104946, z:-0.08187400904630208 \n",
" x:-0.015410075342964763, y:-0.0002086089466902828, z:-0.08403859895720674 \n",
" x:0.002222557388210005, y:0.0033974994957447547, z:-0.06642143213658117 \n",
" x:-0.015582059116631436, y:2.2633004501554508e-15, z:-0.06985538723491842 \n",
" x:-0.01766763896580353, y:0.018868777070110178, z:-0.08187400904630134 \n",
" x:-0.015410075342964766, y:0.00020860894669479622, z:-0.08403859895720513 \n",
" x:-0.015207250272737447, y:2.2496589914265297e-15, z:-0.10102732897327556 "
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"latticeSize=1\n",
"setup=getSetup(latticeSize)\n",
"\n",
"for i in 1:1\n",
" displacementFEA,Load,topNodesIndices=fea(setup)\n",
"end\n",
"# displacementFEA,Load,topNodesIndices=feaDisplacement(setup,latticeSize)\n",
"updateDataAndSaveFEA!(setup,\"../json/trialJuliaParallelGPU.json\")\n",
"displacementFEA"
]
},
{
"cell_type": "code",
"execution_count": 102,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"dt: 0.0023068357674191054\n",
"first timestep took 0.9859363 seconds\n",
"ran latticeSize 1 with 44 voxels and 104 edges for 50000 time steps took 16.158550899 seconds\n",
"FEA displacement= -0.10102732897327556,converged displacement= -0.08487774089695395\n"
]
},
{
"data": {
"image/svg+xml": [
"<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
"<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
"<defs>\n",
" <clipPath id=\"clip9800\">\n",
" <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip9800)\" d=\"\n",
"M0 1600 L2400 1600 L2400 0 L0 0 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<defs>\n",
" <clipPath id=\"clip9801\">\n",
" <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip9800)\" d=\"\n",
"M270.627 1425.62 L2352.76 1425.62 L2352.76 121.675 L270.627 121.675 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<defs>\n",
" <clipPath id=\"clip9802\">\n",
" <rect x=\"270\" y=\"121\" width=\"2083\" height=\"1305\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 325.626,1425.62 325.626,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 718.481,1425.62 718.481,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1111.34,1425.62 1111.34,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1504.19,1425.62 1504.19,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1897.04,1425.62 1897.04,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 2289.9,1425.62 2289.9,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,1248.06 2352.76,1248.06 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,1028.73 2352.76,1028.73 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,809.405 2352.76,809.405 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,590.075 2352.76,590.075 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,370.745 2352.76,370.745 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 270.627,151.415 2352.76,151.415 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,1425.62 2352.76,1425.62 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,1425.62 270.627,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 325.626,1425.62 325.626,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 718.481,1425.62 718.481,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1111.34,1425.62 1111.34,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1504.19,1425.62 1504.19,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1897.04,1425.62 1897.04,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2289.9,1425.62 2289.9,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,1248.06 295.612,1248.06 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,1028.73 295.612,1028.73 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,809.405 295.612,809.405 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,590.075 295.612,590.075 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,370.745 295.612,370.745 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 270.627,151.415 295.612,151.415 \n",
" \"/>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 325.626, 1479.62)\" x=\"325.626\" y=\"1479.62\">0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 718.481, 1479.62)\" x=\"718.481\" y=\"1479.62\">100</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1111.34, 1479.62)\" x=\"1111.34\" y=\"1479.62\">200</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1504.19, 1479.62)\" x=\"1504.19\" y=\"1479.62\">300</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1897.04, 1479.62)\" x=\"1897.04\" y=\"1479.62\">400</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2289.9, 1479.62)\" x=\"2289.9\" y=\"1479.62\">500</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 246.627, 1265.56)\" x=\"246.627\" y=\"1265.56\">-1.00</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 246.627, 1046.23)\" x=\"246.627\" y=\"1046.23\">-0.75</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 246.627, 826.905)\" x=\"246.627\" y=\"826.905\">-0.50</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 246.627, 607.575)\" x=\"246.627\" y=\"607.575\">-0.25</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 246.627, 388.245)\" x=\"246.627\" y=\"388.245\">0.00</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 246.627, 168.915)\" x=\"246.627\" y=\"168.915\">0.25</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:84px; text-anchor:middle;\" transform=\"rotate(0, 1311.69, 73.2)\" x=\"1311.69\" y=\"73.2\">1 Voxel Convergence Study</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(0, 1311.69, 1559.48)\" x=\"1311.69\" y=\"1559.48\">timestep</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:66px; text-anchor:middle;\" transform=\"rotate(-90, 89.2861, 773.647)\" x=\"89.2861\" y=\"773.647\">displacement</text>\n",
"</g>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 329.555,370.745 333.483,370.745 337.412,370.872 341.34,371.472 345.269,372.771 349.197,374.94 353.126,378.087 357.055,382.271 360.983,387.525 364.912,393.861 \n",
" 368.84,401.279 372.769,409.77 376.697,419.32 380.626,429.908 384.554,441.513 388.483,454.111 392.411,467.674 396.34,482.176 400.269,497.589 404.197,513.883 \n",
" 408.126,531.027 412.054,548.99 415.983,567.74 419.911,587.241 423.84,607.459 427.768,628.354 431.697,649.887 435.625,672.016 439.554,694.698 443.483,717.884 \n",
" 447.411,741.527 451.34,765.576 455.268,789.977 459.197,814.674 463.125,839.611 467.054,864.727 470.982,889.961 474.911,915.25 478.839,940.53 482.768,965.735 \n",
" 486.697,990.798 490.625,1015.65 494.554,1040.23 498.482,1064.46 502.411,1088.29 506.339,1111.63 510.268,1134.42 514.196,1156.61 518.125,1178.12 522.053,1198.88 \n",
" 525.982,1218.85 529.911,1237.95 533.839,1256.14 537.768,1273.35 541.696,1289.54 545.625,1304.65 549.553,1318.64 553.482,1331.46 557.41,1343.07 561.339,1353.44 \n",
" 565.267,1362.52 569.196,1370.3 573.125,1376.74 577.053,1381.82 580.982,1385.52 584.91,1387.82 588.839,1388.71 592.767,1388.19 596.696,1386.24 600.624,1382.87 \n",
" 604.553,1378.07 608.481,1371.85 612.41,1364.23 616.339,1355.21 620.267,1344.8 624.196,1333.04 628.124,1319.93 632.053,1305.5 635.981,1289.78 639.91,1272.81 \n",
" 643.838,1254.6 647.767,1235.2 651.695,1214.65 655.624,1192.98 659.553,1170.24 663.481,1146.46 667.41,1121.7 671.338,1095.99 675.267,1069.39 679.195,1041.95 \n",
" 683.124,1013.71 687.052,984.726 690.981,955.049 694.909,924.731 698.838,893.823 702.767,862.379 706.695,830.451 710.624,798.092 714.552,765.355 718.481,732.291 \n",
" 722.409,698.954 726.338,665.394 730.266,631.773 734.195,598.49 738.123,565.823 742.052,533.996 745.981,503.167 749.909,473.444 753.838,444.906 757.766,417.61 \n",
" 761.695,391.599 765.623,366.907 769.552,343.558 773.48,321.574 777.409,300.968 781.337,281.753 785.266,263.936 789.195,247.524 793.123,232.52 797.052,218.924 \n",
" 800.98,206.735 804.909,195.947 808.837,186.554 812.766,178.545 816.694,171.908 820.623,166.624 824.551,162.675 828.48,160.036 832.409,158.681 836.337,158.579 \n",
" 840.266,159.696 844.194,161.995 848.123,165.433 852.051,169.967 855.98,175.549 859.908,182.128 863.837,189.652 867.766,198.063 871.694,207.303 875.623,217.313 \n",
" 879.551,228.029 883.48,239.388 887.408,251.325 891.337,263.774 895.265,276.67 899.194,289.944 903.122,303.531 907.051,317.365 910.98,331.38 914.908,345.512 \n",
" 918.837,359.697 922.765,373.873 926.694,387.981 930.622,401.961 934.551,415.759 938.479,429.321 942.408,442.596 946.336,455.535 950.265,468.094 954.194,480.229 \n",
" 958.122,491.903 962.051,503.078 965.979,513.722 969.908,523.805 973.836,533.303 977.765,542.191 981.693,550.45 985.622,558.066 989.55,565.024 993.479,571.317 \n",
" 997.408,576.937 1001.34,581.883 1005.26,586.153 1009.19,589.752 1013.12,592.686 1017.05,594.963 1020.98,596.596 1024.91,597.597 1028.84,597.984 1032.76,597.776 \n",
" 1036.69,596.992 1040.62,595.656 1044.55,593.793 1048.48,591.429 1052.41,588.59 1056.34,585.306 1060.26,581.607 1064.19,577.525 1068.12,573.089 1072.05,568.334 \n",
" 1075.98,563.291 1079.91,557.994 1083.84,552.475 1087.76,546.769 1091.69,540.909 1095.62,534.926 1099.55,528.853 1103.48,522.723 1107.41,516.565 1111.34,510.411 \n",
" 1115.26,504.29 1119.19,498.229 1123.12,492.257 1127.05,486.398 1130.98,480.678 1134.91,475.12 1138.84,469.746 1142.76,464.575 1146.69,459.628 1150.62,454.922 \n",
" 1154.55,450.471 1158.48,446.29 1162.41,442.391 1166.34,438.786 1170.26,435.483 1174.19,432.49 1178.12,429.812 1182.05,427.454 1185.98,425.419 1189.91,423.708 \n",
" 1193.83,422.32 1197.76,421.253 1201.69,420.505 1205.62,420.07 1209.55,419.944 1213.48,420.118 1217.41,420.585 1221.33,421.335 1225.26,422.358 1229.19,423.642 \n",
" 1233.12,425.177 1237.05,426.948 1240.98,428.942 1244.91,431.144 1248.83,433.541 1252.76,436.117 1256.69,438.856 1260.62,441.743 1264.55,444.76 1268.48,447.893 \n",
" 1272.41,451.125 1276.33,454.439 1280.26,457.82 1284.19,461.251 1288.12,464.717 1292.05,468.201 1295.98,471.69 1299.91,475.168 1303.83,478.62 1307.76,482.034 \n",
" 1311.69,485.396 1315.62,488.692 1319.55,491.913 1323.48,495.045 1327.41,498.078 1331.33,501.003 1335.26,503.81 1339.19,506.491 1343.12,509.039 1347.05,511.446 \n",
" 1350.98,513.708 1354.91,515.817 1358.83,517.771 1362.76,519.566 1366.69,521.198 1370.62,522.666 1374.55,523.969 1378.48,525.107 1382.41,526.079 1386.33,526.886 \n",
" 1390.26,527.531 1394.19,528.015 1398.12,528.342 1402.05,528.515 1405.98,528.538 1409.9,528.416 1413.83,528.154 1417.76,527.757 1421.69,527.232 1425.62,526.585 \n",
" 1429.55,525.823 1433.48,524.952 1437.4,523.98 1441.33,522.914 1445.26,521.763 1449.19,520.534 1453.12,519.234 1457.05,517.872 1460.98,516.455 1464.9,514.993 \n",
" 1468.83,513.491 1472.76,511.96 1476.69,510.405 1480.62,508.835 1484.55,507.258 1488.48,505.68 1492.4,504.108 1496.33,502.549 1500.26,501.009 1504.19,499.496 \n",
" 1508.12,498.014 1512.05,496.568 1515.98,495.165 1519.9,493.809 1523.83,492.505 1527.76,491.255 1531.69,490.065 1535.62,488.937 1539.55,487.875 1543.48,486.88 \n",
" 1547.4,485.955 1551.33,485.102 1555.26,484.322 1559.19,483.615 1563.12,482.983 1567.05,482.425 1570.98,481.942 1574.9,481.533 1578.83,481.196 1582.76,480.932 \n",
" 1586.69,480.738 1590.62,480.613 1594.55,480.554 1598.48,480.56 1602.4,480.627 1606.33,480.754 1610.26,480.937 1614.19,481.173 1618.12,481.459 1622.05,481.791 \n",
" 1625.97,482.167 1629.9,482.581 1633.83,483.031 1637.76,483.512 1641.69,484.022 1645.62,484.555 1649.55,485.109 1653.47,485.679 1657.4,486.261 1661.33,486.852 \n",
" 1665.26,487.448 1669.19,488.045 1673.12,488.64 1677.05,489.229 1680.97,489.809 1684.9,490.377 1688.83,490.929 1692.76,491.462 1696.69,491.975 1700.62,492.463 \n",
" 1704.55,492.926 1708.47,493.36 1712.4,493.763 1716.33,494.133 1720.26,494.469 1724.19,494.77 1728.12,495.033 1732.05,495.258 1735.97,495.443 1739.9,495.589 \n",
" 1743.83,495.694 1747.76,495.758 1751.69,495.78 1755.62,495.762 1759.55,495.702 1763.47,495.602 1767.4,495.461 1771.33,495.281 1775.26,495.062 1779.19,494.805 \n",
" 1783.12,494.511 1787.05,494.181 1790.97,493.817 1794.9,493.42 1798.83,492.991 1802.76,492.533 1806.69,492.046 1810.62,491.533 1814.55,490.994 1818.47,490.433 \n",
" 1822.4,489.851 1826.33,489.25 1830.26,488.632 1834.19,487.998 1838.12,487.351 1842.04,486.692 1845.97,486.024 1849.9,485.349 1853.83,484.668 1857.76,483.983 \n",
" 1861.69,483.296 1865.62,482.608 1869.54,481.922 1873.47,481.239 1877.4,480.561 1881.33,479.888 1885.26,479.223 1889.19,478.566 1893.12,477.92 1897.04,477.285 \n",
" 1900.97,476.662 1904.9,476.052 1908.83,475.456 1912.76,474.875 1916.69,474.31 1920.62,473.761 1924.54,473.228 1928.47,472.713 1932.4,472.215 1936.33,471.734 \n",
" 1940.26,471.272 1944.19,470.827 1948.12,470.4 1952.04,469.99 1955.97,469.599 1959.9,469.224 1963.83,468.866 1967.76,468.525 1971.69,468.2 1975.62,467.89 \n",
" 1979.54,467.595 1983.47,467.315 1987.4,467.048 1991.33,466.794 1995.26,466.553 1999.19,466.322 2003.12,466.102 2007.04,465.892 2010.97,465.691 2014.9,465.497 \n",
" 2018.83,465.311 2022.76,465.131 2026.69,464.956 2030.62,464.786 2034.54,464.619 2038.47,464.455 2042.4,464.292 2046.33,464.131 2050.26,463.97 2054.19,463.808 \n",
" 2058.11,463.645 2062.04,463.48 2065.97,463.312 2069.9,463.14 2073.83,462.965 2077.76,462.785 2081.69,462.6 2085.61,462.409 2089.54,462.212 2093.47,462.009 \n",
" 2097.4,461.799 2101.33,461.581 2105.26,461.357 2109.19,461.124 2113.11,460.884 2117.04,460.636 2120.97,460.379 2124.9,460.115 2128.83,459.843 2132.76,459.562 \n",
" 2136.69,459.274 2140.61,458.977 2144.54,458.673 2148.47,458.362 2152.4,458.043 2156.33,457.718 2160.26,457.385 2164.19,457.047 2168.11,456.702 2172.04,456.352 \n",
" 2175.97,455.997 2179.9,455.637 2183.83,455.272 2187.76,454.904 2191.69,454.532 2195.61,454.157 2199.54,453.78 2203.47,453.401 2207.4,453.02 2211.33,452.639 \n",
" 2215.26,452.257 2219.19,451.875 2223.11,451.493 2227.04,451.113 2230.97,450.734 2234.9,450.357 2238.83,449.982 2242.76,449.61 2246.69,449.242 2250.61,448.877 \n",
" 2254.54,448.516 2258.47,448.16 2262.4,447.809 2266.33,447.463 2270.26,447.122 2274.18,446.787 2278.11,446.459 2282.04,446.136 2285.97,445.82 2289.9,445.512 \n",
" 2293.83,445.21 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9802)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 329.555,459.378 333.483,459.378 337.412,459.378 341.34,459.378 345.269,459.378 349.197,459.378 353.126,459.378 357.055,459.378 360.983,459.378 364.912,459.378 \n",
" 368.84,459.378 372.769,459.378 376.697,459.378 380.626,459.378 384.554,459.378 388.483,459.378 392.411,459.378 396.34,459.378 400.269,459.378 404.197,459.378 \n",
" 408.126,459.378 412.054,459.378 415.983,459.378 419.911,459.378 423.84,459.378 427.768,459.378 431.697,459.378 435.625,459.378 439.554,459.378 443.483,459.378 \n",
" 447.411,459.378 451.34,459.378 455.268,459.378 459.197,459.378 463.125,459.378 467.054,459.378 470.982,459.378 474.911,459.378 478.839,459.378 482.768,459.378 \n",
" 486.697,459.378 490.625,459.378 494.554,459.378 498.482,459.378 502.411,459.378 506.339,459.378 510.268,459.378 514.196,459.378 518.125,459.378 522.053,459.378 \n",
" 525.982,459.378 529.911,459.378 533.839,459.378 537.768,459.378 541.696,459.378 545.625,459.378 549.553,459.378 553.482,459.378 557.41,459.378 561.339,459.378 \n",
" 565.267,459.378 569.196,459.378 573.125,459.378 577.053,459.378 580.982,459.378 584.91,459.378 588.839,459.378 592.767,459.378 596.696,459.378 600.624,459.378 \n",
" 604.553,459.378 608.481,459.378 612.41,459.378 616.339,459.378 620.267,459.378 624.196,459.378 628.124,459.378 632.053,459.378 635.981,459.378 639.91,459.378 \n",
" 643.838,459.378 647.767,459.378 651.695,459.378 655.624,459.378 659.553,459.378 663.481,459.378 667.41,459.378 671.338,459.378 675.267,459.378 679.195,459.378 \n",
" 683.124,459.378 687.052,459.378 690.981,459.378 694.909,459.378 698.838,459.378 702.767,459.378 706.695,459.378 710.624,459.378 714.552,459.378 718.481,459.378 \n",
" 722.409,459.378 726.338,459.378 730.266,459.378 734.195,459.378 738.123,459.378 742.052,459.378 745.981,459.378 749.909,459.378 753.838,459.378 757.766,459.378 \n",
" 761.695,459.378 765.623,459.378 769.552,459.378 773.48,459.378 777.409,459.378 781.337,459.378 785.266,459.378 789.195,459.378 793.123,459.378 797.052,459.378 \n",
" 800.98,459.378 804.909,459.378 808.837,459.378 812.766,459.378 816.694,459.378 820.623,459.378 824.551,459.378 828.48,459.378 832.409,459.378 836.337,459.378 \n",
" 840.266,459.378 844.194,459.378 848.123,459.378 852.051,459.378 855.98,459.378 859.908,459.378 863.837,459.378 867.766,459.378 871.694,459.378 875.623,459.378 \n",
" 879.551,459.378 883.48,459.378 887.408,459.378 891.337,459.378 895.265,459.378 899.194,459.378 903.122,459.378 907.051,459.378 910.98,459.378 914.908,459.378 \n",
" 918.837,459.378 922.765,459.378 926.694,459.378 930.622,459.378 934.551,459.378 938.479,459.378 942.408,459.378 946.336,459.378 950.265,459.378 954.194,459.378 \n",
" 958.122,459.378 962.051,459.378 965.979,459.378 969.908,459.378 973.836,459.378 977.765,459.378 981.693,459.378 985.622,459.378 989.55,459.378 993.479,459.378 \n",
" 997.408,459.378 1001.34,459.378 1005.26,459.378 1009.19,459.378 1013.12,459.378 1017.05,459.378 1020.98,459.378 1024.91,459.378 1028.84,459.378 1032.76,459.378 \n",
" 1036.69,459.378 1040.62,459.378 1044.55,459.378 1048.48,459.378 1052.41,459.378 1056.34,459.378 1060.26,459.378 1064.19,459.378 1068.12,459.378 1072.05,459.378 \n",
" 1075.98,459.378 1079.91,459.378 1083.84,459.378 1087.76,459.378 1091.69,459.378 1095.62,459.378 1099.55,459.378 1103.48,459.378 1107.41,459.378 1111.34,459.378 \n",
" 1115.26,459.378 1119.19,459.378 1123.12,459.378 1127.05,459.378 1130.98,459.378 1134.91,459.378 1138.84,459.378 1142.76,459.378 1146.69,459.378 1150.62,459.378 \n",
" 1154.55,459.378 1158.48,459.378 1162.41,459.378 1166.34,459.378 1170.26,459.378 1174.19,459.378 1178.12,459.378 1182.05,459.378 1185.98,459.378 1189.91,459.378 \n",
" 1193.83,459.378 1197.76,459.378 1201.69,459.378 1205.62,459.378 1209.55,459.378 1213.48,459.378 1217.41,459.378 1221.33,459.378 1225.26,459.378 1229.19,459.378 \n",
" 1233.12,459.378 1237.05,459.378 1240.98,459.378 1244.91,459.378 1248.83,459.378 1252.76,459.378 1256.69,459.378 1260.62,459.378 1264.55,459.378 1268.48,459.378 \n",
" 1272.41,459.378 1276.33,459.378 1280.26,459.378 1284.19,459.378 1288.12,459.378 1292.05,459.378 1295.98,459.378 1299.91,459.378 1303.83,459.378 1307.76,459.378 \n",
" 1311.69,459.378 1315.62,459.378 1319.55,459.378 1323.48,459.378 1327.41,459.378 1331.33,459.378 1335.26,459.378 1339.19,459.378 1343.12,459.378 1347.05,459.378 \n",
" 1350.98,459.378 1354.91,459.378 1358.83,459.378 1362.76,459.378 1366.69,459.378 1370.62,459.378 1374.55,459.378 1378.48,459.378 1382.41,459.378 1386.33,459.378 \n",
" 1390.26,459.378 1394.19,459.378 1398.12,459.378 1402.05,459.378 1405.98,459.378 1409.9,459.378 1413.83,459.378 1417.76,459.378 1421.69,459.378 1425.62,459.378 \n",
" 1429.55,459.378 1433.48,459.378 1437.4,459.378 1441.33,459.378 1445.26,459.378 1449.19,459.378 1453.12,459.378 1457.05,459.378 1460.98,459.378 1464.9,459.378 \n",
" 1468.83,459.378 1472.76,459.378 1476.69,459.378 1480.62,459.378 1484.55,459.378 1488.48,459.378 1492.4,459.378 1496.33,459.378 1500.26,459.378 1504.19,459.378 \n",
" 1508.12,459.378 1512.05,459.378 1515.98,459.378 1519.9,459.378 1523.83,459.378 1527.76,459.378 1531.69,459.378 1535.62,459.378 1539.55,459.378 1543.48,459.378 \n",
" 1547.4,459.378 1551.33,459.378 1555.26,459.378 1559.19,459.378 1563.12,459.378 1567.05,459.378 1570.98,459.378 1574.9,459.378 1578.83,459.378 1582.76,459.378 \n",
" 1586.69,459.378 1590.62,459.378 1594.55,459.378 1598.48,459.378 1602.4,459.378 1606.33,459.378 1610.26,459.378 1614.19,459.378 1618.12,459.378 1622.05,459.378 \n",
" 1625.97,459.378 1629.9,459.378 1633.83,459.378 1637.76,459.378 1641.69,459.378 1645.62,459.378 1649.55,459.378 1653.47,459.378 1657.4,459.378 1661.33,459.378 \n",
" 1665.26,459.378 1669.19,459.378 1673.12,459.378 1677.05,459.378 1680.97,459.378 1684.9,459.378 1688.83,459.378 1692.76,459.378 1696.69,459.378 1700.62,459.378 \n",
" 1704.55,459.378 1708.47,459.378 1712.4,459.378 1716.33,459.378 1720.26,459.378 1724.19,459.378 1728.12,459.378 1732.05,459.378 1735.97,459.378 1739.9,459.378 \n",
" 1743.83,459.378 1747.76,459.378 1751.69,459.378 1755.62,459.378 1759.55,459.378 1763.47,459.378 1767.4,459.378 1771.33,459.378 1775.26,459.378 1779.19,459.378 \n",
" 1783.12,459.378 1787.05,459.378 1790.97,459.378 1794.9,459.378 1798.83,459.378 1802.76,459.378 1806.69,459.378 1810.62,459.378 1814.55,459.378 1818.47,459.378 \n",
" 1822.4,459.378 1826.33,459.378 1830.26,459.378 1834.19,459.378 1838.12,459.378 1842.04,459.378 1845.97,459.378 1849.9,459.378 1853.83,459.378 1857.76,459.378 \n",
" 1861.69,459.378 1865.62,459.378 1869.54,459.378 1873.47,459.378 1877.4,459.378 1881.33,459.378 1885.26,459.378 1889.19,459.378 1893.12,459.378 1897.04,459.378 \n",
" 1900.97,459.378 1904.9,459.378 1908.83,459.378 1912.76,459.378 1916.69,459.378 1920.62,459.378 1924.54,459.378 1928.47,459.378 1932.4,459.378 1936.33,459.378 \n",
" 1940.26,459.378 1944.19,459.378 1948.12,459.378 1952.04,459.378 1955.97,459.378 1959.9,459.378 1963.83,459.378 1967.76,459.378 1971.69,459.378 1975.62,459.378 \n",
" 1979.54,459.378 1983.47,459.378 1987.4,459.378 1991.33,459.378 1995.26,459.378 1999.19,459.378 2003.12,459.378 2007.04,459.378 2010.97,459.378 2014.9,459.378 \n",
" 2018.83,459.378 2022.76,459.378 2026.69,459.378 2030.62,459.378 2034.54,459.378 2038.47,459.378 2042.4,459.378 2046.33,459.378 2050.26,459.378 2054.19,459.378 \n",
" 2058.11,459.378 2062.04,459.378 2065.97,459.378 2069.9,459.378 2073.83,459.378 2077.76,459.378 2081.69,459.378 2085.61,459.378 2089.54,459.378 2093.47,459.378 \n",
" 2097.4,459.378 2101.33,459.378 2105.26,459.378 2109.19,459.378 2113.11,459.378 2117.04,459.378 2120.97,459.378 2124.9,459.378 2128.83,459.378 2132.76,459.378 \n",
" 2136.69,459.378 2140.61,459.378 2144.54,459.378 2148.47,459.378 2152.4,459.378 2156.33,459.378 2160.26,459.378 2164.19,459.378 2168.11,459.378 2172.04,459.378 \n",
" 2175.97,459.378 2179.9,459.378 2183.83,459.378 2187.76,459.378 2191.69,459.378 2195.61,459.378 2199.54,459.378 2203.47,459.378 2207.4,459.378 2211.33,459.378 \n",
" 2215.26,459.378 2219.19,459.378 2223.11,459.378 2227.04,459.378 2230.97,459.378 2234.9,459.378 2238.83,459.378 2242.76,459.378 2246.69,459.378 2250.61,459.378 \n",
" 2254.54,459.378 2258.47,459.378 2262.4,459.378 2266.33,459.378 2270.26,459.378 2274.18,459.378 2278.11,459.378 2282.04,459.378 2285.97,459.378 2289.9,459.378 \n",
" 2293.83,459.378 \n",
" \"/>\n",
"<path clip-path=\"url(#clip9800)\" d=\"\n",
"M1853.56 386.635 L2280.76 386.635 L2280.76 205.195 L1853.56 205.195 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1853.56,386.635 2280.76,386.635 2280.76,205.195 1853.56,205.195 1853.56,386.635 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#009af9; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1877.56,265.675 2021.56,265.675 \n",
" \"/>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2045.56, 283.175)\" x=\"2045.56\" y=\"283.175\">Dynamic</text>\n",
"</g>\n",
"<polyline clip-path=\"url(#clip9800)\" style=\"stroke:#e26f46; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1877.56,326.155 2021.56,326.155 \n",
" \"/>\n",
"<g clip-path=\"url(#clip9800)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2045.56, 343.655)\" x=\"2045.56\" y=\"343.655\">FEA</text>\n",
"</g>\n",
"</svg>\n"
]
},
"execution_count": 102,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"setup=getSetup(latticeSize)\n",
"numTimeSteps=50000 #30000\n",
"displacements=[]\n",
"save=true\n",
"returnEvery=100\n",
"runMetavoxelGPU!(setup,numTimeSteps,latticeSize,displacements,returnEvery,true)\n",
"\n",
"numTimeStepsRecorded=length(displacements)\n",
"d=[]\n",
"dFEA=[]\n",
"j=length(displacements[end])\n",
"step=1\n",
"for i in 1:step:numTimeStepsRecorded\n",
" append!(d,displacements[i][j].z)\n",
" append!(dFEA,displacementFEA[j].z)\n",
"end\n",
"\n",
"\n",
"println(\"FEA displacement= $(displacementFEA[j].z),converged displacement= $(displacements[numTimeStepsRecorded][j].z)\")\n",
"plot(1:step:numTimeStepsRecorded,d,label=\"Dynamic\",xlabel=\"timestep\",ylabel=\"displacement\",title=\"$latticeSize Voxel Convergence Study\")\n",
"plot!(1:step:numTimeStepsRecorded,dFEA,label=\"FEA\")\n",
"# savefig(\"4_voxel_convergence\")\n"
]
},
{
"cell_type": "code",
"execution_count": 222,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n",
"0.03141075907812829\n",
"0.06279051952931337\n",
"0.09410831331851431\n",
"0.12533323356430426\n",
"0.15643446504023087\n",
"0.1873813145857246\n",
"0.21814324139654256\n",
"0.2486898871648548\n",
"0.2789911060392293\n",
"0.3090169943749474\n",
"0.33873792024529137\n",
"0.3681245526846779\n",
"0.3971478906347806\n",
"0.4257792915650727\n",
"0.45399049973954675\n",
"0.4817536741017153\n",
"0.5090414157503713\n",
"0.5358267949789967\n",
"0.5620833778521306\n",
"0.5877852522924731\n",
"0.6129070536529764\n",
"0.6374239897486896\n",
"0.6613118653236518\n",
"0.6845471059286886\n",
"0.7071067811865475\n",
"0.7289686274214116\n",
"0.7501110696304596\n",
"0.7705132427757893\n",
"0.7901550123756903\n",
"0.8090169943749475\n",
"0.8270805742745618\n",
"0.8443279255020151\n",
"0.8607420270039436\n",
"0.8763066800438637\n",
"0.8910065241883678\n",
"0.9048270524660196\n",
"0.9177546256839811\n",
"0.9297764858882513\n",
"0.9408807689542255\n",
"0.9510565162951535\n",
"0.960293685676943\n",
"0.9685831611286311\n",
"0.9759167619387473\n",
"0.9822872507286886\n",
"0.9876883405951378\n",
"0.9921147013144779\n",
"0.99556196460308\n",
"0.9980267284282716\n",
"0.9995065603657316\n",
"1.0\n",
"0.9995065603657316\n",
"0.9980267284282716\n",
"0.99556196460308\n",
"0.9921147013144778\n",
"0.9876883405951377\n",
"0.9822872507286886\n",
"0.9759167619387474\n",
"0.9685831611286312\n",
"0.9602936856769431\n",
"0.9510565162951536\n",
"0.9408807689542255\n",
"0.9297764858882515\n",
"0.9177546256839813\n",
"0.9048270524660195\n",
"0.8910065241883679\n",
"0.8763066800438635\n",
"0.8607420270039436\n",
"0.844327925502015\n",
"0.827080574274562\n",
"0.8090169943749475\n",
"0.7901550123756905\n",
"0.7705132427757893\n",
"0.7501110696304596\n",
"0.7289686274214114\n",
"0.7071067811865476\n",
"0.6845471059286888\n",
"0.6613118653236518\n",
"0.6374239897486899\n",
"0.6129070536529764\n",
"0.5877852522924732\n",
"0.5620833778521305\n",
"0.535826794978997\n",
"0.5090414157503714\n",
"0.4817536741017156\n",
"0.45399049973954686\n",
"0.4257792915650729\n",
"0.3971478906347806\n",
"0.36812455268467814\n",
"0.3387379202452913\n",
"0.3090169943749475\n",
"0.2789911060392291\n",
"0.24868988716485482\n",
"0.21814324139654231\n",
"0.18738131458572502\n",
"0.15643446504023098\n",
"0.12533323356430454\n",
"0.09410831331851435\n",
"0.06279051952931358\n",
"0.031410759078128236\n",
"1.2246467991473532e-16\n"
]
}
],
"source": [
"length(displacements)\n",
"for i in 0:100\n",
" println(sin(i/100*pi))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 156,
"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=\"clip8200\">\n",
" <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip8200)\" 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=\"clip8201\">\n",
" <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<path clip-path=\"url(#clip8200)\" d=\"\n",
"M243.864 1425.62 L2352.76 1425.62 L2352.76 121.675 L243.864 121.675 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<defs>\n",
" <clipPath id=\"clip8202\">\n",
" <rect x=\"243\" y=\"121\" width=\"2110\" height=\"1305\"/>\n",
" </clipPath>\n",
"</defs>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 257.282,1425.62 257.282,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 719.961,1425.62 719.961,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1182.64,1425.62 1182.64,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 1645.32,1425.62 1645.32,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 2108,1425.62 2108,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 243.864,1384.36 2352.76,1384.36 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 243.864,1149.17 2352.76,1149.17 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 243.864,913.988 2352.76,913.988 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 243.864,678.803 2352.76,678.803 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 243.864,443.617 2352.76,443.617 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8202)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
" 243.864,208.432 2352.76,208.432 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,1425.62 2352.76,1425.62 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,1425.62 243.864,121.675 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 257.282,1425.62 257.282,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 719.961,1425.62 719.961,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1182.64,1425.62 1182.64,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1645.32,1425.62 1645.32,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 2108,1425.62 2108,1409.97 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,1384.36 269.171,1384.36 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,1149.17 269.171,1149.17 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,913.988 269.171,913.988 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,678.803 269.171,678.803 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,443.617 269.171,443.617 \n",
" \"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 243.864,208.432 269.171,208.432 \n",
" \"/>\n",
"<g clip-path=\"url(#clip8200)\">\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, 257.282, 1479.62)\" x=\"257.282\" y=\"1479.62\">0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 719.961, 1479.62)\" x=\"719.961\" y=\"1479.62\">10</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 1182.64, 1479.62)\" x=\"1182.64\" y=\"1479.62\">20</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 1645.32, 1479.62)\" x=\"1645.32\" y=\"1479.62\">30</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 2108, 1479.62)\" x=\"2108\" y=\"1479.62\">40</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 219.864, 1401.86)\" x=\"219.864\" y=\"1401.86\">-0.4</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 219.864, 1166.67)\" x=\"219.864\" y=\"1166.67\">-0.3</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 219.864, 931.488)\" x=\"219.864\" y=\"931.488\">-0.2</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 219.864, 696.303)\" x=\"219.864\" y=\"696.303\">-0.1</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 219.864, 461.117)\" x=\"219.864\" y=\"461.117\">0.0</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 219.864, 225.932)\" x=\"219.864\" y=\"225.932\">0.1</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 1298.31, 73.2)\" x=\"1298.31\" y=\"73.2\">Node Displacement</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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, 1298.31, 1559.48)\" x=\"1298.31\" y=\"1559.48\">Node ID</text>\n",
"</g>\n",
"<g clip-path=\"url(#clip8200)\">\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(#clip8202)\" cx=\"303.55\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"349.818\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"396.086\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"442.354\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"488.621\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"534.889\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"581.157\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"627.425\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"673.693\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"719.961\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"766.229\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"812.497\" cy=\"443.617\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"858.765\" cy=\"386.411\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"905.033\" cy=\"529.988\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"951.301\" cy=\"472.725\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"997.569\" cy=\"405.767\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1043.84\" cy=\"348.284\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1090.1\" cy=\"549.237\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1136.37\" cy=\"511.87\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1182.64\" cy=\"635.72\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1228.91\" cy=\"279.302\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1275.18\" cy=\"790.762\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1321.44\" cy=\"597.1\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1367.71\" cy=\"264.563\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1413.98\" cy=\"195.219\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1460.25\" cy=\"775.555\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1506.52\" cy=\"528.717\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1552.78\" cy=\"932.699\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1599.05\" cy=\"197.374\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1645.32\" cy=\"1095.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(#clip8202)\" cx=\"1691.59\" cy=\"698.619\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1737.86\" 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(#clip8202)\" cx=\"1784.12\" cy=\"159.792\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1830.39\" cy=\"1058.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(#clip8202)\" cx=\"1876.66\" cy=\"687.088\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1922.93\" cy=\"1221.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(#clip8202)\" cx=\"1969.19\" cy=\"165.966\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2015.46\" cy=\"1245.27\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2061.73\" cy=\"710.23\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2108\" cy=\"288.088\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2154.27\" cy=\"297.239\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2200.53\" cy=\"1370.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(#clip8202)\" cx=\"2246.8\" cy=\"974.309\" r=\"14\" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2293.07\" 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(#clip8202)\" cx=\"303.55\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"349.818\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"396.086\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"442.354\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"488.621\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"534.889\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"581.157\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"627.425\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"673.693\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"719.961\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"766.229\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"812.497\" cy=\"443.617\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"858.765\" cy=\"468.527\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"905.033\" cy=\"468.527\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"951.301\" cy=\"493.792\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"997.569\" cy=\"468.527\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1043.84\" cy=\"493.079\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1090.1\" cy=\"468.527\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1136.37\" cy=\"493.792\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1182.64\" cy=\"493.079\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1228.91\" cy=\"518.35\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1275.18\" cy=\"518.349\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1321.44\" cy=\"542.529\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1367.71\" cy=\"518.35\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1413.98\" cy=\"543.999\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1460.25\" cy=\"518.349\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1506.52\" cy=\"542.529\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1552.78\" cy=\"543.997\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1599.05\" cy=\"568.135\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1645.32\" cy=\"568.219\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1691.59\" cy=\"594.248\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1737.86\" cy=\"568.135\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1784.12\" cy=\"591.847\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1830.39\" cy=\"568.219\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1876.66\" cy=\"594.248\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1922.93\" cy=\"592.019\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"1969.19\" cy=\"599.831\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2015.46\" cy=\"636.173\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2061.73\" cy=\"641.264\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2108\" cy=\"599.831\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2154.27\" cy=\"607.907\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2200.53\" cy=\"636.173\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2246.8\" cy=\"641.264\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<circle clip-path=\"url(#clip8202)\" cx=\"2293.07\" cy=\"681.219\" r=\"14\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<path clip-path=\"url(#clip8200)\" d=\"\n",
"M1853.56 386.635 L2280.76 386.635 L2280.76 205.195 L1853.56 205.195 Z\n",
" \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
"<polyline clip-path=\"url(#clip8200)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
" 1853.56,386.635 2280.76,386.635 2280.76,205.195 1853.56,205.195 1853.56,386.635 \n",
" \"/>\n",
"<circle clip-path=\"url(#clip8200)\" 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(#clip8200)\">\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",
"<circle clip-path=\"url(#clip8200)\" cx=\"1961.56\" cy=\"326.155\" r=\"21\" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"1\" stroke=\"#000000\" stroke-opacity=\"1\" stroke-width=\"3.2\"/>\n",
"<g clip-path=\"url(#clip8200)\">\n",
"<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 2045.56, 343.655)\" x=\"2045.56\" y=\"343.655\">FEA</text>\n",
"</g>\n",
"</svg>\n"
]
},
"execution_count": 156,
"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].z)\n",
" append!(nFEA,displacementFEA[i].z)\n",
"end\n",
"scatter(1:j,n,label=\"Dyanmic\",xlabel=\"Node ID\",ylabel=\"displacement\",title=\"Node Displacement\")\n",
"scatter!(1:j,nFEA,label=\"FEA\")\n",
"# savefig(\"node_displacement_four_voxel\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"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
}