function FEM_truss(data,A)
# FEM solver for truss elements
# Author: JV Carstensen, CEE, MIT (JK Guest, Civil Eng, JHU)
# Revised: Aug 22 2017, JVC
# Revised: Amira Abdel-Rahman
# =======================================================================
# ---- READ IN DATA -----------------------------------------------------
# read in external file
E,f,dis,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te= data;
g=zeros(size(dis))
g.=dis
# ---- NUMBER THE EQUATIONS ---------------------------------------------
# line 380
id,neq = number_eq(idb,ndf,nnp);
# ---- FORM THE ELEMENT STIFFNESS MATRICES ------------------------------
# line 324
nee = ndf*nen; # number of element equations
Ke = zeros(nee,nee,nel);
Te = zeros(nen*1,nen*nsd,nel); # *1 is specific to truss
for i = 1:nel
Ke[:,:,i],Te[:,:,i] = Ke_truss(A[i],E[i],ien[:,i],nee,nsd,xn);
end
# display("Element Stiffness Matrix")
# display(Ke)
# ---- PERFORM GLOBAL TO LOCAL MAPPING ----------------------------------
# line 237
LM = zeros(nee,nel);
for i = 1:nel
LM[:,i] = get_local_id(id,ien[:,i],ndf,nee,nen);
end
# ---- IF THERE IS FREE DEGREES OF FREEDOM, THEN SOLVE THE EQUILIBRIUM --
if (neq > 0)
# get global force vector (line 275)
F = globalF(f,g,id,ien,Ke,LM,ndf,nee,nel,nen,neq,nnp);
# solve Kd - F = 0 (line 521)
d,K = solveEQ(F,LM,Ke,nee,nel,neq);
end
# display("F")
# display(F)
# ---- POST-PROCESS THE RESULTS -----------------------------------------
# line 419
dcomp,axial,stress,strain,Fe = post_processing(A,d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp,Te);
# display("dcomp")
# display(dcomp)
# display("axial")
# display(axial)
# display("stress")
# display(stress)
# ---- COMPUTE REACTION FORCES ------------------------------------------
# line 482
Rcomp,idr = reactions(idb,ien,Fe,ndf,nee,nel,nen,nnp);
# display("Rcomp")
# display(Rcomp)
# ---- PLOT THE STRUCTURE -----------------------------------------------
# read in external file
# set the plot factor for the thickness of the truss elements
# plot_fac_bar = 1/A[i];
# A_min = 0;
# plot_truss(A,A_min,f,idr,ien,nel,nnp,nsd,plot_fac_bar,xn);
return K,F,d,stress,dcomp,g
# =======================================================================
end
# =======================================================================
function add_d2dcomp(dcomp,d,id,ndf,nnp)
# function that adds the displacements of the free degrees of freedom to
# the nodal displacements
# -----------------------------------------------------------------------
# dcomp(ndf,nnp) = nodal displacements
# d(neq,1) = displacement at free degrees of freedom
# id(ndf,nnp) = equation numbers of degrees of freedom
# ndf = number of degrees of freedom per node
# nnp = number of nodal points
#------------------------------------------------------------------------
# loop over nodes and degrees of freedom
for n=1:nnp
for i=1:ndf
# if it is a free dof then add the global displacement
if (id[i,n]>0)
dcomp[Int(i),Int(n)] = dcomp[Int(i),Int(n)]+d[Int(id[Int(i),Int(n)])];
end
end
end
return dcomp
end
# =======================================================================
function add_loads_to_force(F,f,id,ndf,nnp)
# function that adds nodal forces to the global force vector
# -----------------------------------------------------------------------
# F(neq,1) = global force vector
# f(ndf,nnp) = prescribed nodal forces
# id(ndf,nnp) = equation numbers of degrees of freedom
# ndf = number of degrees of freedom per node
# nnp = number of nodal points
#------------------------------------------------------------------------
# loop over nodes and degrees of freedom
for n = 1:nnp
for i = 1:ndf
# get the global equation number
M = id[i,n];
# if free degree of freedom, then add nodal load to global force
# vector
if (M > 0)
F[Int(floor(M))] = F[Int(floor(M))] + f[Int(floor(i)),Int(floor(n))];
end
end
end
return F
end
# =======================================================================
function addforce(F,Fe,LM,nee)
# function that adds element forces to the global force vector
# -----------------------------------------------------------------------
# F(neq,1) = global force vector
# Fe(nee,1) = element force vector
# LM(nee,nel) = global to local map for the element
# nee = number of element equations
# =======================================================================
# loop over rows of Fe
for i = 1:nee
# get the global equation number for local equation i
M = LM[i];
# if free dof (eqn number > 0) add to F vector
if (M > 0)
F[Int(M)]=F[Int(M)]+Fe[Int(i)];
end
end
return F
end
# =======================================================================
function get_de_from_dcomp(dcomp,ien,ndf,nen)
# extracts element displacement vector from complete displacement vector
#------------------------------------------------------------------------
# dcomp(ndf,nnp) = nodal displacements
# ien(nen,1) = element connectivity
# ndf = number of degrees of freedom per node
# nen = number of element equations
#
# de(nen,1) = element displacements
#
#------------------------------------------------------------------------
de = zeros((nen-1)*ndf+ndf,1);
# loop over number of element nodes
for i = 1:nen
# loop over number of degrees of freedom per node
for j = 1:ndf
# get the local element number and place displacement in de
leq = (i-1)*ndf+j;
de[leq] = dcomp[j,Int(floor(ien[i]))];
end
end
return de
end
# =======================================================================
function get_local_id(id,ien,ndf,nee,nen)
# functions that performs the global to local mapping of equation numbers
# -----------------------------------------------------------------------
# id(ndf,nnp) = equation numbers of degrees of freedom
# ien(nen,1) = element connectivity
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nen = number of element equations
#
# LM(nee,1) = global to local map for element
# =======================================================================
# initialize global-local mapping matrix
LM = zeros(nee,1);
# initialize local equation number counter
k = 0;
# loop over element nodes
for i = 1:nen
# loop over degrees of freedom at each node
for j = 1:ndf
# update counter and prescribe global equation number
k = k+1;
LM[k] = id[j,Int(floor(ien[i]))];
end
end
return LM
end
# =======================================================================
function globalF(f,g,id,ien,Ke,LM,ndf,nee,nel,nen,neq,nnp)
# function that assembles the global load vector
# -----------------------------------------------------------------------
# id(ndf,nnp) = equation numbers of degrees of freedom
# f(ndf,nnp) = prescribed nodal forces
# g(ndf,nnp) = prescribed nodal displacements
# ien(nen,nel) = element connectivities
# Ke(nee,nee,nel) = element stiffness matrices
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nel = number of elements
# nen = number of element equations
# neq = number of equations
# nnp = number of nodal points
#
# F(neq,1) = global force vector
# =======================================================================
# initialize
F = zeros(neq,1);
# Insert applied loads into F
F = add_loads_to_force(F,f,id,ndf,nnp);
# Compute forces from applied displacements (ds~=0) and insert into F
Fse = zeros(nee,nel);
# loop over elements
for i = 1:nel
# get dse for current element
dse = get_de_from_dcomp(g,ien[:,i],ndf,nen);
# compute element force
Fse[:,i] = -Ke[:,:,i]*dse;
# assemble elem force into global force vector
F = addforce(F,Fse[:,i],LM[:,i],nee);
end
return F
end
# =======================================================================
function Ke_truss(A,E,ien,nee,nsd,xn)
# function that computes the global element stiffness matrix for a truss
# element
# -----------------------------------------------------------------------
# A(1,1) = cross-sectional area of elements
# E(1,1) = Young's modulus of elements
# ien(nen,1) = element connectivity
# nee = number of element equations
# nsd = number of spacial dimensions
# xn(nsd,nnp) = nodal coordinates
#
# Ke(nee,nee,1) = global element stiffness matrix
# Te(nee,nee,1) = global to local transformation matrix for element
# =======================================================================
# form vector along axis of element using nodal coordinates
v = xn[:,Int(floor(ien[2]))]-xn[:,Int(floor(ien[1]))];
# compute the length of the element
Le = norm(v,2);
# rotation of parent domain
# rot=[ cos(theta_x) cos(theta_y) cos(theta_z) ]'
rot = v/Le;
# local element stiffness matrix
ke = E*A/Le*[ 1 -1
-1 1 ];
# Transformation matrix: global to local coordinate system
if (nsd == 2) # 2D case
# truss Te is nen x ndf*nen array
Te = [ rot[1] rot[2] 0 0
0 0 rot[1] rot[2] ];
elseif (nsd == 3) # 3D case
# Truss Te is nen x ndf*nen array
Te = [ rot[1] rot[2] rot[3] 0 0 0
0 0 0 rot[1] rot[2] rot[3] ];
end
# compute the global element stiffness matrix
Ke = zeros(nee,nee);
Ke = Te'*ke*Te;
return Ke,Te
end
# =======================================================================
function number_eq(idb,ndf,nnp)
# function that numbers the unknown degrees of freedom (equations)
# -----------------------------------------------------------------------
# idb(ndf,nnp) = 1 if the degree of freedom is prescribed, 0 otherwise
# ndf = number of degrees of freedom per node
# nnp = number of nodal points
#
# id(ndf,nnp) = equation numbers of degrees of freedom
# neq = number of equations (tot number of degrees of freedom)
# =======================================================================
# initialize id and neq
id = zeros(ndf,nnp);
neq = 0;
# loop over nodes
for n = 1:nnp
# loop over degrees of freedom
for i = 1:ndf
if idb[i,n] == 0
# udate # of equations
neq = neq + 1;
# if no prescribed displacement at dof i of node n
# => give an equation # different from 0
id[i,n] = neq;
end
end
end
return id,neq
end
# =======================================================================
function post_processing(A,d,E,g,id,ien,Ke,ndf,nee,nel,nen,nnp,Te)
# function that performs post processing for truss elements
# -----------------------------------------------------------------------
# A(nel,1) = cross-sectional area of elements
# d(neq,1) = displacements at free degrees of freedom
# E(nel,1) = Young's modulus of elements
# g(ndf,nnp) = prescribed nodal displacements
# id(ndf,nnp) = equation numbers of degrees of freedom
# ien(nen,nel) = element connectivities
# Ke(nee,nee,nel) = element stiffness matrices
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nel = number of elements
# nen = number of element equations
# nnp = number of nodes
# Te(nee,nee,nel) = element transformation matrices
#
# dcomp(ndf,nnp) = nodal displacements
# axial(nel,1) = axial element forces
# stress(nel,1) = element stresses
# strain(nel,1) = element strains
# Fe(nee,nel) = element forces
# =======================================================================
# get the total displacement of the structure in matrix form dcomp(nsd,nnp)
dcomp = add_d2dcomp(g,d,id,ndf,nnp);
# initalize evaluation of global element forces Fe, local element forces
# fe, axial forces, element stresses and strains
Fe = zeros(nee,nel);
de = zeros(nee,nel);
fe = zeros(1*nen,nel); # element local force vector
axial = zeros(nel,1);
stress = zeros(nel,1);
strain = zeros(nel,1); # element axial, stress, strain
# loop over elements
for i=1:nel
# get the element displacaments
de[:,i] = get_de_from_dcomp(dcomp,ien[:,i],ndf,nen);
# compute the element forces
Fe[:,i] = Ke[:,:,i]*de[:,i];
# transform Fe to the local coordinate system
fe[:,i] = Te[:,:,i]*Fe[:,i];
# Compute the axial force, stress, strain
axial[i] = fe[2,i] ; # Use second entry for truss element
stress[i] = axial[i]/A[i];
strain[i] = stress[i]/E[i];
end
# display("Fe")
# display(Fe)
return dcomp,axial,stress,strain,Fe
end
# =======================================================================
function reactions(idb,ien,Fe,ndf,nee,nel,nen,nnp)
# function that computes the reaction forces on the structure
# -----------------------------------------------------------------------
# idb(ndf,nnp) = 1 if the degree of freedom is prescribed, 0 otherwise
# ien(nen,nel) = element connectivities
# Fe(nee,nel) = element forces
# ndf = number of degrees of freedom per node
# nee = number of element equations
# nel = number of elements
# nen = number of element equations
# nnp = number of nodes
#
# Rcomp(ndf,nnp) = nodal reactions
# idbr(ndf,nnp) = 0 if the degree of freedom is prescribed, otherwise
# =======================================================================
# switch BC marker and number the equations for the reaction forces
idbr = 1 .- idb;
idr,neqr = number_eq(idbr,ndf,nnp);
# assemble reactions R from element force vectors Fe
R = zeros(neqr,1);
for i = 1:nel
LMR = get_local_id(idr,ien[:,i],ndf,nee,nen);
R = addforce(R,Fe[:,i],LMR,nee);
end
# organize the reactions in matrix array Rcomp(ndf,nnp)
Rcomp = zeros(ndf,nnp);
Rcomp = add_d2dcomp(Rcomp,R,idr,ndf,nnp);
return Rcomp,idr
end
# =======================================================================
function addstiff(K,Ke,LM,nee)
# function that solves the equilibrium condition
# -----------------------------------------------------------------------
# K(neq,neq) = global stiffness matrix
# Ke(nee,nee,1) = element stiffness matrix
# LM(nee,nel) = global to local map for the element
# nee = number of element equations
# =======================================================================
# loop over rows of Ke
for i =1:nee
# loop over columns of Ke
for j = 1:nee
Mr = LM[i];
Mc = LM[j];
if(Mr > 0 && Mc > 0)
# if equation #'s are non-zero add element contribution to the
# stiffness matrix
K[Int(Mr),Int(Mc)] = K[Int(Mr),Int(Mc)] +Ke[Int(i),Int(j)];
end
end
end
return K
end
# =======================================================================
function solveEQ(F,LM,Ke,nee,nel,neq)
# function that solves the equilibrium condition
# -----------------------------------------------------------------------
# F(neq,1) = global force vector
# LM(nee,nel) = global to local maps
# Ke(nee,nee,nel) = element stiffness matrices
# nee = number of element equations
# nel = number of elements
# neq = number of equations
#
# d(neq,1) = displacements at free degrees of freedom
# =======================================================================
# assemble global stiffness matrix
# K = zeros(neq,neq); # Use 'sparse' for more efficient memory usage
K=spzeros(neq,neq)
for i = 1:nel
K = addstiff(K,Ke[:,:,i],LM[:,i],nee);
end
# display("K")
# display(K)
# solve the equlibrium
d = K\F;
# display("d")
# display(d)
return d,K
end
# =======================================================================
# =======================================================================
function data_truss()
# -----------------------------------------------------------------------
# A(nel,1) = cross-sectional area of elements
# E(nel,1) = Young's modulus of elements
# f(ndf,nnp) = prescribed nodal forces
# g(ndf,nnp) = prescribed nodal displacements
# idb(ndf,nnp) = 1 if the degree of freedom is prescribed, 0 otherwise
# ien(nen,nel) = element connectivities
# ndf = number of degrees of freedom per node
# nel = number of elements
# nen = number of element equations
# nnp = number of nodes
# nsd = number of spacial dimensions
# xn(nsd,nnp) = nodal coordinates
# =======================================================================
# ---- MESH -------------------------------------------------------------
nsd = 2; # number of spacial dimensions
ndf = 2; # number of degrees of freedom per node
nen = 2; # number of element nodes
nel = 3; # number of elements
nnp = 3; # number of nodal points
# ---- NODAL COORDINATES ------------------------------------------------
# xn(i,N) = coordinate i for node N
# N = 1,...,nnp
# i = 1,...,nsd
# -----------------------------------------------------------------------
xn = zeros(nsd,nnp);
xn[1:2,1] = [ 0 0]';
xn[1:2,2] = [ 4/tand(30) 4]';
xn[1:2,3] = [ 4+4/tand(30) 0]';
# ---- NODAL COORDINATES ------------------------------------------------
# ien(a,e) = N
# N: global node number - N = 1,...,nnp
# e: element number - e = 1,...,nel
# a: local node number - a = 1,...,nen
# -----------------------------------------------------------------------
ien = zeros(nen,nel);
ien[1:2,1] = [1 2]' ;
ien[1:2,2] = [2 3]' ;
ien[1:2,3] = [3 1]' ;
# ---- MATERIAL PROPERTIES ----------------------------------------------
# E(e) = E_mat
# e: element number - e = 1,...,nel
# -----------------------------------------------------------------------
E_mat = 200; # Young's modulus # lbf/in^2
E = E_mat*ones(nel,1);
# ---- GEOMETRIC PROPERTIES ---------------------------------------------
# A(e) = A_bar
# e: element number - e = 1,...,nel
# -----------------------------------------------------------------------
A = 1e3*[15; 20; 18]; # mm^2
# ---- BOUNDARY CONDITIONS ----------------------------------------------
# prescribed displacement flags (essential boundary condition)
#
# idb(i,N) = 1 if the degree of freedom i of the node N is prescribed
# = 0 otherwise
#
# 1) initialize idb to 0
# 2) enter the flag for prescribed displacement boundary conditions
# -----------------------------------------------------------------------
idb = zeros(ndf,nnp);
idb[:,1].=1;
idb[2,3]=1;
# ---- BOUNDARY CONDITIONS: PRESCRIBED NODAL DISPLACEMENTS --------------
# g(i,N) = prescribed displacement magnitude
# N = 1,...,nnp
# i = 1,...,nsd
#
# 1) initialize g to 0
# 2) enter the values
# -----------------------------------------------------------------------
g = zeros(ndf,nnp);
# ---- BOUNDARY CONDITIONS: PRESCRIBED NODAL FORCES ---------------------
# f(i,N) = prescribed force magnitude
# N = 1,...,nnp
# i = 1,...,nsd
#
# 1) initialize f to 0
# 2) enter the values
# -----------------------------------------------------------------------
P = 500e3; # lbf
f = zeros(ndf,nnp);
f[1,2] = P*cosd(40);
f[2,2] = P*sind(40);
return A,E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn
end
# =======================================================================
function data_truss2(A)
# -----------------------------------------------------------------------
# A(nel,1) = cross-sectional area of elements
# E(nel,1) = Young's modulus of elements
# f(ndf,nnp) = prescribed nodal forces
# g(ndf,nnp) = prescribed nodal displacements
# idb(ndf,nnp) = 1 if the degree of freedom is prescribed, 0 otherwise
# ien(nen,nel) = element connectivities
# ndf = number of degrees of freedom per node
# nel = number of elements
# nen = number of element equations
# nnp = number of nodes
# nsd = number of spacial dimensions
# xn(nsd,nnp) = nodal coordinates
# =======================================================================
# ---- MESH -------------------------------------------------------------
nsd = 2; # number of spacial dimensions
ndf = 2; # number of degrees of freedom per node
nen = 2; # number of element nodes
nel = 17; # number of elements
nnp = 10; # number of nodal points
# ---- NODAL COORDINATES ------------------------------------------------
# xn(i,N) = coordinate i for node N
# N = 1,...,nnp
# i = 1,...,nsd
# -----------------------------------------------------------------------
L= 30.0/4.0;#*0.3048;
H= 10.0;#*0.3048;
xn = zeros(nsd,nnp);
xn[1:2, 1] = [ 0 0]';# 1
xn[1:2, 2] = [ L 0]';# 2
xn[1:2, 3] = [ 2*L 0]';# 3
xn[1:2, 4] = [ 3*L 0]';# 4
xn[1:2, 5] = [ 4*L 0]';# 5
xn[1:2, 6] = [ 0 -H]';# 6
xn[1:2, 7] = [ L -H]';# 7
xn[1:2, 8] = [ 2*L -H]';# 8
xn[1:2, 9] = [ 3*L -H]';# 9
xn[1:2,10] = [ 4*L -H]';# 10
# ---- NODAL COORDINATES ------------------------------------------------
# ien(a,e) = N
# N: global node number - N = 1,...,nnp
# e: element number - e = 1,...,nel
# a: local node number - a = 1,...,nen
# -----------------------------------------------------------------------
ien = zeros(nen,nel);
ien[1:2,1] = [1 2]' ;#
ien[1:2,2] = [2 3]' ;#
ien[1:2,3] = [3 4]' ;#
ien[1:2,4] = [4 5]' ;#
ien[1:2,5] = [6 7]' ;#
ien[1:2,6] = [7 8]' ;#
ien[1:2,7] = [8 9]' ;#
ien[1:2,8] = [9 10]' ;#
ien[1:2, 9] = [1 6]' ;#
ien[1:2,10] = [2 6]' ;#
ien[1:2,11] = [2 7]' ;#
ien[1:2,12] = [2 8]' ;#
ien[1:2,13] = [3 8]' ;#
ien[1:2,14] = [4 8]' ;#
ien[1:2,15] = [4 9]' ;#
ien[1:2,16] = [4 10]' ;#
ien[1:2,17] = [5 10]' ;#
# ---- MATERIAL PROPERTIES ----------------------------------------------
# E(e) = E_mat
# e: element number - e = 1,...,nel
# -----------------------------------------------------------------------
E_mat = 29000*1000; # Young's modulus # lbf/in^2
E = E_mat*ones(nel,1);
# ---- GEOMETRIC PROPERTIES ---------------------------------------------
# A(e) = A_bar
# e: element number - e = 1,...,nel
# -----------------------------------------------------------------------
#A = ones(nel);
#A[:] .= 400; # mm^2
# ---- BOUNDARY CONDITIONS ----------------------------------------------
# prescribed displacement flags (essential boundary condition)
#
# idb(i,N) = 1 if the degree of freedom i of the node N is prescribed
# = 0 otherwise
#
# 1) initialize idb to 0
# 2) enter the flag for prescribed displacement boundary conditions
# -----------------------------------------------------------------------
idb = zeros(ndf,nnp);
idb[:,6].=1;
idb[2,10]=1;
# ---- BOUNDARY CONDITIONS: PRESCRIBED NODAL DISPLACEMENTS --------------
# g(i,N) = prescribed displacement magnitude
# N = 1,...,nnp
# i = 1,...,nsd
#
# 1) initialize g to 0
# 2) enter the values
# -----------------------------------------------------------------------
g = zeros(ndf,nnp);
# ---- BOUNDARY CONDITIONS: PRESCRIBED NODAL FORCES ---------------------
# f(i,N) = prescribed force magnitude
# N = 1,...,nnp
# i = 1,...,nsd
#
# 1) initialize f to 0
# 2) enter the values
# -----------------------------------------------------------------------
P = 20.0* 1000; # lbf
f = zeros(ndf,nnp);
f[2,7] = -P;
f[2,8] = -P;
f[2,9] = -P;
return [A,E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn]
end
# =======================================================================
function mapp(value, x1, y1, x2, y2)
return (value - x1) * (y2 - x2) / (y1 - x1) + x2;
end
# =======================================================================
function plotTruss(problem,X,scale,threshold=0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
linewidth = X[i]*scale)
# linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
linewidth = 0.0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
p
end
# =======================================================================
function getSetup(fileName)
setup = Dict()
open(fileName, "r") do f
dicttxt = read(f,String) # file information to string
setup=JSON.parse(dicttxt) # parse and transform data
end
return setup
end
# =======================================================================
function plotTrussDeformed(problem,X,scale,threshold=0,exageration=100.0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
xnn=zeros(size(xn))
xnn.=xn .+ exageration.*dcomp
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
p=plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
[xnn[2,Int(ien[1,i])],xnn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,1.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
linewidth = X[i]*scale)
# linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
linewidth = 0.0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
p
end
function plotTrussDeformed3D(problem,X,scale,threshold=0,exageration=100.0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
xnn=zeros(size(xn))
xnn.=xn .+ exageration.*dcomp
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
p=plot!([xnn[1,Int(ien[1,i])],xnn[1,Int(ien[2,i])]],
[xnn[3,Int(ien[1,i])],xnn[3,Int(ien[2,i])]],
[xnn[2,Int(ien[1,i])],xnn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,1.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
linewidth = X[i]*scale)
# linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
linewidth = 0.0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
p
end
function plotTruss3D(problem,X,scale,threshold=0)
nel=length(X)
E,f,g,idb,ien,ndf,nel,nen,nnp,nsd,xn,len,Ke,Te=problem;
K,F,d,stress,dcomp,g=FEM_truss(problem,X);
p=plot(axis=nothing,ticks=nothing, border=nothing, aspect_ratio=:equal)
for i in 1:nel
if X[i]>threshold
if threshold>0
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(0.0,0.0,0.0),
linewidth = 3.0)
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 1, 0),0.0, mapp(stress[i], minimum(stress[:]), maximum(stress[:]), 0, 1)),
linewidth = X[i]*scale)
# linewidth = mapp(X[i], minimum(X[:]), maximum(X[:])+0.0001, scale/10.0, scale))
# p=annotate!((xn[1,Int(ien[1,i])]+xn[1,Int(ien[2,i])])/2.0, (xn[2,Int(ien[1,i])]+xn[2,Int(ien[2,i])])/2.0, "$(floor(stress[i]/1000))*e3", 8)
end
else
p=plot!([xn[1,Int(ien[1,i])],xn[1,Int(ien[2,i])]],
[xn[3,Int(ien[1,i])],xn[3,Int(ien[2,i])]],
[xn[2,Int(ien[1,i])],xn[2,Int(ien[2,i])]],label="",
color=RGB(1.0,1.0,1.0),
linewidth = 0.0)
end
end
# plot!(axis=nothing,ticks=nothing, border=nothing)
p
end