-
Amira Abdel-Rahman authoredAmira Abdel-Rahman authored
updateNodes.jl 16.94 KiB
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
# BASED ON https://github.com/jonhiller/Voxelyze
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_fixedDisplacement,N_moment,N_orient,N_edgeID,N_edgeFirst,N_material,N_poissonStrain,
E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,E_axis, E_strain,E_material)
index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
stride = blockDim().x * gridDim().x
## @cuprintln("thread $index, block $stride")
N,M=size(N_edgeID)
for i = index:stride:N
isFloorEnabled=floorEnabled();#todo make these arguments ??
# @inbounds if lengthVector3(Vector3(convert(Float64,N_displacement[i].x),
# convert(Float64,N_displacement[i].y),convert(Float64,N_displacement[i].z))) >1000.0 #todo ?? change this according to ratio
# # @cuprintln("DIVERGEDNODE!!!!!!!!!!")
# return
# end
@inbounds if N_restrained[i]
@inbounds N_linMom[i]=Vector3(0,0,0)
@inbounds N_angMom[i]=Vector3(0,0,0)
@inbounds if(N_material[i].poisson)
@inbounds N_poissonStrain[i]=strain(i,N_material[i].poisson,i, N_edgeID, N_edgeFirst, E_axis, E_strain, E_material, N_restrained[i], N_material[i])
end
return
elseif N_fixedDisplacement[i].x != 0.0 || N_fixedDisplacement[i].y != 0.0 || N_fixedDisplacement[i].z != 0.0
@inbounds N_linMom[i]=Vector3(0,0,0)
@inbounds N_angMom[i]=Vector3(0,0,0)
@inbounds if(N_material[i].poisson)
@inbounds N_poissonStrain[i]=strain(i,N_material[i].poisson,i, N_edgeID, N_edgeFirst, E_axis, E_strain, E_material, N_restrained[i], N_material[i])
end
@inbounds translate=externalDisplacement(currentTimeStep,N_position,N_fixedDisplacement[i])
@inbounds N_currPosition[i]=N_currPosition[i]+translate
@inbounds N_displacement[i]=N_displacement[i]+translate
else
for j in 1:M
@inbounds if (N_edgeID[i,j]!=-1)
#@cuprintln("i $i, j $j, N_edgeID[i,j] $temp")
@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]] )
@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]] )
end
end
FloorStaticFriction=false;
prec=1e8
#get properties
@inbounds E=roundd(N_material[i].E,prec)
@inbounds nomSize=roundd(N_material[i].nomSize,prec)
@inbounds mass=roundd(N_material[i].mass,prec)
@inbounds massInverse=roundd(N_material[i].massInverse,prec)
@inbounds curForce = force(i,N_intForce[i],N_orient[i],N_force[i],N_position[i],currentTimeStep,N_material[i],N_linMom[i])
fricForce = Vector3(curForce.x,curForce.y,curForce.z);
if (isFloorEnabled)
@inbounds curForce,FloorStaticFriction=floorForce!(dt,curForce,N_currPosition[i],
Vector3(N_linMom[i].x,N_linMom[i].y ,N_linMom[i].z),FloorStaticFriction,N_material[i])
end
fricForce = curForce-fricForce;
# if(i==25)
# x=N_intForce[i].x*1e6
# y=N_intForce[i].y*1e6
# z=N_intForce[i].z*1e6
# @cuprintln("N_intForce[i] x $x 1e-6, y $y 1e-6, z $z 1e-6")
# x=curForce.x*1e6
# y=curForce.y*1e6
# z=curForce.z*1e6
# @cuprintln("curForce x $x 1e-6, y $y 1e-6, z $z 1e-6")
# x=N_linMom[i].x*1e6
# y=N_linMom[i].y*1e6
# z=N_linMom[i].z*1e6
# @cuprintln("N_linMom[i] x $x 1e-6, y $y 1e-6, z $z 1e-6")
# x=translate.x*1e0
# y=translate.y*1e0
# z=translate.z*1e0
# @cuprintln("translate x $x, y $y, z $z")
# @cuprintln("")
# end
# x=N_orient[i].x*1e6
# y=N_orient[i].y*1e6
# z=N_orient[i].z*1e6
# @cuprintln("N_orient[i] x $x 1e-6, y $y 1e-6, z $z 1e-6")
# x=N_linMom[i].x*1e6
# y=N_linMom[i].y*1e6
# z=N_linMom[i].z*1e6
# @cuprintln("N_linMom[i] x $x 1e-6, y $y 1e-6, z $z 1e-6")
# x=curForce.x*1e6
# y=curForce.y*1e6
# z=curForce.z*1e6
# @cuprintln("curForce x $x 1e-6, y $y 1e-6, z $z 1e-6")
# @cuprintln("dt $dt")
# @cuprintln("massInverse $massInverse")
#########################################
# @inbounds N_intForce[i]=Vector3(0,0,0) #??
@inbounds N_linMom[i]=N_linMom[i]+curForce*Vector3(dt,dt,dt) #todo make sure right
@inbounds translate=N_linMom[i]*Vector3((dt*massInverse),(dt*massInverse),(dt*massInverse)) # ??massInverse
############################
#we need to check for friction conditions here (after calculating the translation) and stop things accordingly
@inbounds if (isFloorEnabled && floorPenetration( convert(Float64,N_currPosition[i].x),convert(Float64,N_currPosition[i].y),N_material[i].nomSize)>= 0.0 )#we must catch a slowing voxel here since it all boils down to needing access to the dt of this timestep.
work =convert(Float64,fricForce.x*translate.x + fricForce.z*translate.z); #F dot disp
@inbounds hKe = 0.5*massInverse*convert(Float64,(N_linMom[i].x*N_linMom[i].x + N_linMom[i].z*N_linMom[i].z)); #horizontal kinetic energy
if((hKe + work) <= 0.0)
FloorStaticFriction=true; #this checks for a change of direction according to the work-energy principle
end
if(FloorStaticFriction) #if we're in a state of static friction, zero out all horizontal motion
N_linMom[i]=Vector3(0.0 ,convert(Float64,N_linMom[i].y) ,0.0)
translate=Vector3(0.0 ,convert(Float64,translate.y) ,0.0)
end
else
FloorStaticFriction=false
end
@inbounds N_currPosition[i]=N_currPosition[i]+translate
@inbounds N_displacement[i]=N_displacement[i]+translate
# Rotation
@inbounds curMoment = moment(N_intMoment[i],N_orient[i],N_moment[i],N_material[i],N_angMom[i])
@inbounds N_intMoment[i]=Vector3(0,0,0) # do i really need it?
@inbounds N_intForce[i]=Vector3(0,0,0) # do i really need it?
@inbounds N_angMom[i]=N_angMom[i]+curMoment*Vector3(dt,dt,dt)
@inbounds momentInertiaInverse=N_material[i].momentInertiaInverse
@inbounds temp=FromRotationVector(N_angMom[i]*Vector3((dt*momentInertiaInverse),(dt*momentInertiaInverse),(dt*momentInertiaInverse)))
@inbounds N_orient[i]=multiplyQuaternions(temp,N_orient[i])
@inbounds N_angle[i]=ToRotationVector(N_orient[i])
# x=N_orient[i].x*1e0
# y=N_orient[i].y*1e0
# z=N_orient[i].z*1e0
# @cuprintln("N_orient[i] x $x, y $y, z $z ")
# @cuprintln("momentInertiaInverse $momentInertiaInverse")
# if (ext){
# double size = mat->nominalSize();
# if (ext->isFixed(X_TRANSLATE)) {pos.x = ix*size + ext->translation().x; linMom.x=0;}
# if (ext->isFixed(Y_TRANSLATE)) {pos.y = iy*size + ext->translation().y; linMom.y=0;}
# if (ext->isFixed(Z_TRANSLATE)) {pos.z = iz*size + ext->translation().z; linMom.z=0;}
# if (ext->isFixedAnyRotation()){ //if any rotation fixed, all are fixed
# if (ext->isFixedAllRotation()){
# orient = ext->rotationQuat();
# angMom = Vec3D<double>();
# }
# else { //partial fixes: slow!
# Vec3D<double> tmpRotVec = orient.ToRotationVector();
# if (ext->isFixed(X_ROTATE)){ tmpRotVec.x=0; angMom.x=0;}
# if (ext->isFixed(Y_ROTATE)){ tmpRotVec.y=0; angMom.y=0;}
# if (ext->isFixed(Z_ROTATE)){ tmpRotVec.z=0; angMom.z=0;}
# orient.FromRotationVector(tmpRotVec);
# }
# }
# }
# if(i==17)
# x=(N_currPosition[i].x-0.004)*1e6
# y=N_currPosition[i].y*1e6
# z=N_currPosition[i].z*1e6
# @cuprintln("pos 17 x $x 1e-6, y $y 1e-6, z $z 1e-6")
# @cuprintln("")
# end
# if(i==25)
# x=(N_currPosition[i].x)*1e6
# y=N_currPosition[i].y*1e6
# z=N_currPosition[i].z*1e6
# @cuprintln("pos 25 x $x 1e-6, y $y 1e-6, z $z 1e-6") # @cuprintln("")
# end
# if(i==41)
# x=curForce.x*1e6
# y=curForce.y*1e6
# z=curForce.z*1e6
# @cuprintln("curForce x $x 1e-6, y $y 1e-6, z $z 1e-6")
# x=(N_currPosition[i].x-0.004)*1e6
# y=N_currPosition[i].y*1e6
# z=N_currPosition[i].z*1e6
# @cuprintln("pos x $x 1e-6, y $y 1e-6, z $z 1e-6")
# @cuprintln("")
# end
####poisson
@inbounds if(N_material[i].poisson)
@inbounds N_poissonStrain[i]=strain(i,N_material[i].poisson,i, N_edgeID, N_edgeFirst, E_axis, E_strain, E_material, N_restrained[i], N_material[i])
end
end
end
return
end
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_fixedDisplacement,N_moment,N_orient,N_edgeID,N_edgeFirst,N_material,N_poissonStrain,
E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,E_axis, E_strain,E_material)
N=length(N_intForce)
numblocks = ceil(Int, N/256)
CUDA.@sync begin
@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_fixedDisplacement,N_moment,N_orient,N_edgeID,N_edgeFirst,N_material,N_poissonStrain,
E_intForce1,E_intMoment1,E_intForce2,E_intMoment2,E_axis, E_strain,E_material)
end
end
############################
function strain(i,poissonStrain,id,N_edgeID,N_edgeFirst, E_axis, E_strain,E_material,fixed,mat)
###if no connections in the positive and negative directions of a particular axis, strain is zero
###if one connection in positive or negative direction of a particular axis, strain is that strain - ?? and force or constraint?
###if connections in both the positive and negative directions of a particular axis, strain is the average.
# intStrRet= [0.0, 0.0, 0.0] ###intermediate strain return value. axes according to linkAxis enum
# numBondAxis = [0,0,0]; ###number of bonds in this axis (0,1,2). axes according to linkAxis enum
# tension = [false,false,false];
intStrRet1=0.0
intStrRet2=0.0
intStrRet3=0.0
numBondAxis1=0.0
numBondAxis2=0.0
numBondAxis3=0.0
tension1=false
tension2=false
tension3=false
N,M=size(N_edgeID)
for j in 1:M
@inbounds linkID=N_edgeID[id,j]
# @cuprintln("id $id, j $j, linkID $linkID")
if (linkID!=-1)
# if (E_strain[j]!=0.0)
@inbounds axis= normalizeVector3(E_axis[linkID]) @inbounds strain1=axialStrain( N_edgeFirst[id,j],E_strain[linkID],E_material[linkID])
# if (E_strain[j]!=0.0)
# @cuprint("strain1 $strain1")
# end
x= abs(convert(Float64,axis.x))*strain1 #todo change for general axis
y= abs(convert(Float64,axis.y))*strain1 #todo change for general axis
z= abs(convert(Float64,axis.z))*strain1 #todo change for general axis
if (axis.x!=0.0)
intStrRet1+=x
numBondAxis1+=1
# if(i==32)
# @cuprintln("strain1 x $(strain1)")
# end
end
if (axis.y!=0.0)
intStrRet2+=y
numBondAxis2+=1
end
if (axis.z!=0.0)
intStrRet3+=z
numBondAxis3+=1
end
# end
end
end
# @cuprintln("intStrRet 11 x $(intStrRet1*1e6) *1e-6, y $(intStrRet2*1e6) *1e-6, z $(intStrRet3*1e6) *1e-6")
# for i =1:3
if (numBondAxis1>1.0)
intStrRet1=intStrRet1 /numBondAxis1 ###average
end
if (poissonStrain)
tension1 = ((numBondAxis1==2.0) || (fixed&&numBondAxis1==1.0))
#(ext && (numBondAxis[i]==1 && (ext.isFixed((dofComponent)(1<<i)) || ext.force()[i] != 0)))); #if both sides pulling, or just one side and a fixed or forced voxel...
end
if (numBondAxis2>1.0)
intStrRet2=intStrRet2 /numBondAxis2 ###average
end
if (poissonStrain)
tension2 = ((numBondAxis2==2.0) || (fixed&&numBondAxis2==1.0))
#(ext && (numBondAxis[i]==1 && (ext.isFixed((dofComponent)(1<<i)) || ext.force()[i] != 0)))); #if both sides pulling, or just one side and a fixed or forced voxel...
end
if (numBondAxis3>1.0)
intStrRet3=intStrRet3 /numBondAxis3 ###average
end
if (poissonStrain)
tension3 = ((numBondAxis3==2.0) || (fixed&&numBondAxis3==1.0))
#(ext && (numBondAxis[i]==1 && (ext.isFixed((dofComponent)(1<<i)) || ext.force()[i] != 0)))); #if both sides pulling, or just one side and a fixed or forced voxel...
end
# end
# if(i==32)
# @cuprintln("intStrRet 1 x $(intStrRet1) , y $(intStrRet2) , z $(intStrRet3)")
# @cuprintln("numBondAxis x $(numBondAxis1) , y $(numBondAxis2) , z $(numBondAxis3) ")
# end
if (poissonStrain)
if (!(tension1 && tension2 && tension3 )) ###if at least one isn't in tension
add = 0.0;
# for i =1:3
if (tension1)
add=add+intStrRet1
end
if (tension2)
add=add+intStrRet2
end
if (tension3)
add=add+intStrRet3
end # end
# @cuprintln("add x $(add*1e6)")
# @cuprintln("mat.nu x $(mat.nu)")
# value = pow((1.0 + add),(-mat.nu))-1.0 #((1.0 + add)^(-mat.nu)) - 1.0;
value = CUDA.pow( 1.0 + add, -mat.nu) - 1.0;
# if(i==32)
# @cuprintln("value x $(value)")
# @cuprintln("mat.nu x $(mat.nu)")
# @cuprintln("add x $(add)")
# end
# value =((1.0 + add)^(-roundd(mat.nu,10e5))) - 1.0;
# @cuprintln("value x $(value*1e6)")
# value=0.0
# for i =1:3
if (!tension1)
intStrRet1=value;
end
if (!tension2)
intStrRet2=value;
end
if (!tension3)
intStrRet3=value;
end
# end
end
end
# if(i==32)
# @cuprintln("intStrRet 2 x $(intStrRet1*1e6) *1e-6, y $(intStrRet2*1e6) *1e-6, z $(intStrRet3*1e6) *1e-6")
# end
intStrRet= Vector3(intStrRet1,intStrRet2,intStrRet3); ###intermediate strain return value. axes according to linkAxis enum
# @cuprintln("intStrRet x $(intStrRet1*1e6) *1e-6, y $(intStrRet2*1e6) *1e-6, z $(intStrRet3*1e6) *1e-6")
# @cuprintln("numBondAxis x $(numBondAxis1) , y $(numBondAxis2) , z $(numBondAxis3) ")
# if tension1
# tension1=1.0
# else
# tension1=0.0
# end
# if tension2
# tension2=1.0
# else
# tension2=0.0
# end
# if tension3
# tension3=1.0
# else
# tension3=0.0
# end
# @cuprintln("tension x $(tension1) , y $(tension2) , z $(tension3) ")
return intStrRet
end
function axialStrain( positiveEnd,strain,mat)
#strainRatio = pVPos->material()->E/pVNeg->material()->E;
strainRatio=convert(Float64,mat.strainRatio)
# positiveEnd=true
return positiveEnd ? 2.0 *strain*strainRatio/(1.0+strainRatio) : 2.0*strain/(1.0+strainRatio)
end