# // Amira Abdel-Rahman
# // (c) Massachusetts Institute of Technology 2020
using LinearAlgebra
import JSON
function rotation_matrix(element_vector, x_axis, y_axis,z_axis)
# find the direction cosines
L=norm(element_vector)
l = (element_vector[1])/L
m = (element_vector[2])/L
n = (element_vector[3])/L
D = ( l^2+ m^2+n^2)^0.5
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]]
return transMatrix
end
function get_stresses(setup)
nodes = setup["nodes"]
edges = setup["edges"]
ndofs = length(nodes)*6
x_axis = [1 0 0]
y_axis = [0 1 0]
z_axis = [0 0 1]
# find the stresses in each member
stresses=zeros(length(edges))
max11=-10e6
min11=10e6
for edge in edges
#degrees_of_freedom = properties["degrees_of_freedom"]
element=parse(Int,edge["id"][2:end])
# find the nodes that the lements connects
fromNode = nodes[edge["source"]+1]
toNode = nodes[edge["target"]+1]
# the coordinates for each node
fromPoint = [fromNode["position"]["x"] fromNode["position"]["y"] fromNode["position"]["z"]]
toPoint = [toNode["position"]["x"] toNode["position"]["y"] toNode["position"]["z"]]
# find the degrees of freedom for each node
dofs = convert(Array{Int}, fromNode["degrees_of_freedom"])
dofs=vcat(dofs,convert(Array{Int}, toNode["degrees_of_freedom"]))
element_vector=toPoint-fromPoint
# find rotated mass and stifness matrices
tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis)
# i1=parse(Int,fromNode["id"][2:end])
# i2=parse(Int,toNode["id"][2:end])
# 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
global_displacements=[fromNode["displacement"]["x"] fromNode["displacement"]["y"] fromNode["displacement"]["z"] fromNode["angle"]["x"] fromNode["angle"]["y"] fromNode["angle"]["z"] toNode["displacement"]["x"] toNode["displacement"]["y"] toNode["displacement"]["z"] toNode["angle"]["x"] toNode["angle"]["y"] toNode["angle"]["z"]] # todo change
# nodal displacement
q=tau*transpose(global_displacements)
# println(q)
# calculate the strain and stresses
strain =(q[7]-q[1])/norm(element_vector)
E = edge["stiffness"]# youngs modulus
stress=E.*strain
edge["stress"]=stress
if stress>max11
max11=stress
end
if stress<min11
min11=stress
end
# println(element)
# println(stress)
end
setup["viz"]["minStress"]=min11
setup["viz"]["maxStress"]=max11
return stresses
end
function get_matrices(setup)
nodes = setup["nodes"]
edges = setup["edges"]
ndofs = length(nodes)*6
x_axis = [1 0 0]
y_axis = [0 1 0]
z_axis = [0 0 1]
M = zeros((ndofs,ndofs))
K = zeros((ndofs,ndofs))
for edge in edges
#degrees_of_freedom = properties["degrees_of_freedom"]
element=parse(Int,edge["id"][2:end])
# find the nodes that the lements connects
fromNode = nodes[edge["source"]+1]
toNode = nodes[edge["target"]+1]
# the coordinates for each node
fromPoint = [fromNode["position"]["x"] fromNode["position"]["y"] fromNode["position"]["z"]]
toPoint = [toNode["position"]["x"] toNode["position"]["y"] toNode["position"]["z"]]
# find the degrees of freedom for each node
dofs = convert(Array{Int}, fromNode["degrees_of_freedom"])
dofs=vcat(dofs,convert(Array{Int}, toNode["degrees_of_freedom"]))
element_vector=toPoint-fromPoint
# find element mass and stifness matrices
length = norm(element_vector)
rho = edge["density"]
area = edge["area"]
E = edge["stiffness"]# youngs modulus
A = edge["area"]
G=1.0#todo shear_modulus
ixx = 1.0#todo section ixx
iyy = 1.0#todo section.iyy#
l0=length
j=1.0;#todo check
l02 = l0 * l0
l03 = l0 * l0 * l0
# Cm = rho * area * length /6.0
# Ck= E * area / length
# m = [[2 1];[1 2]]
# k = [[1 -1];[-1 1]]
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]]
# find rotated mass and stifness matrices
tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis)
# m_r=transpose(tau)*m*tau
k_r=transpose(tau)*k*tau
# change from element to global coordinate
index= dofs.+1
B=zeros((12,ndofs))
for i in 1:12
B[i,index[i]]=1.0
end
# M_rG= transpose(B)*m_r*B
K_rG= transpose(B)*k_r*B
# M += Cm .* M_rG
# K += Ck .* K_rG
K += K_rG
end
# construct the force vector
F=zeros(ndofs)
remove_indices=[];
for node in nodes
#insert!(F,i, value);
#F=vcat(F,value)
i=parse(Int,node["id"][2:end])
f=node["force"]
# println(f)
F[(i)*6+1]=f["x"]
F[(i)*6+2]=f["y"]
F[(i)*6+3]=f["z"]
F[(i)*6+4]=0
F[(i)*6+5]=0
F[(i)*6+6]=0
dofs = convert(Array{Int}, node["degrees_of_freedom"]).+1
restrained_dofs=node["restrained_degrees_of_freedom"]
for (index, value) in enumerate(dofs)
if restrained_dofs[index]
append!( remove_indices, value)
end
end
end
# println(remove_indices)
# print(K)
# print(F)
# M = M[setdiff(1:end, remove_indices), :]
K = K[setdiff(1:end, remove_indices), :]
# M = M[:,setdiff(1:end, remove_indices)]
K = K[:,setdiff(1:end, remove_indices)]
F = F[setdiff(1:end, remove_indices)]
return M,K,F
end
function updateDisplacement(setup, X)
nodes= setup["nodes"]
i=0
for node in nodes
if !node["restrained_degrees_of_freedom"][1]
#i=parse(Int,node["id"][2:end])
node["displacement"]["x"]=X[(i)*6+1]
node["displacement"]["y"]=X[(i)*6+2]
node["displacement"]["z"]=X[(i)*6+3]
node["angle"]["x"]=X[(i)*6+4]
node["angle"]["y"]=X[(i)*6+5]
node["angle"]["z"]=X[(i)*6+6]
i=i+1
end
end
end
function solveFea(setup)
// # determine the global matrices
M,K,F=get_matrices(setup)
#println(M)
#println(K)
#println(F)
#evals=eigvals(K,M)
#evecs=eigvecs(K,M)
#frequencies=sqrt.(evals)
X=inv(K)*F
#updateDisplacement(displacements);
updateDisplacement(setup, X)
# determine the stresses in each element
stresses=get_stresses(setup)
end
###################
### Read data #####
###################
setup = Dict()
open("./json/setupChiral3.json", "r") do f
global setup
dicttxt = String(read(f)) # file information to string
setup=JSON.parse(dicttxt) # parse and transform data
end
setup=setup["setup"]
###################
### Solve FEA data #####
###################
solveFea(setup)
###################
### Write data #####
###################
# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(setup)
# write the file with the stringdata variable information
open("./json/trialJulia.json", "w") do f
write(f, stringdata)
end