using LinearAlgebra; using SparseArrays; using Plots; using ArrayFire; #Structs struct Material E::Float64 #Young's Modulus v::Float64 #Poisson's Ratio Area::Float64 THIC::Float64 #Thickness end struct Element Type::Int64 #Plane Stress/Strain Node::Array{Float64,2} #Nodes in Element Material::Int64 #Material Number NINT::Int64 #Gauss Numerical Integration Order LM::Array{Int64,1} #DOF Indices NN::Array{Int64,1} #Node Numbers end struct ElementGroup NUME::Int64 #Number of Elements Type::Int64 #Quads/OFE Ele::Array{Element,1} #Array of Elements end #Setting up FEM function OFESolve(IIN; verbose = false, tradmode = false, useGPU = false) global Verbose = verbose; global line = 1; f = open(IIN); global text = readlines(f); close(f); global p = plot(0,0); global usingOFE = false; global TradMode = tradmode; global UseGPU = useGPU; global NumMonomials = 6; #3 = Linear, 4 = Bilinear, 6 = Quadratic if (~verbose) print("Control Information\t\t\t"); end @time RCI(); if (~verbose) print("Nodal Inputs\t\t\t\t"); end @time ReadInput(); if (~verbose) print("Resolve Number of Independent Equations\t"); end @time IDDOF(); if (~verbose) print("Read Load Information\t\t\t"); end @time Loads(); if (~verbose) print("Read Element Information\t\t"); end @time ReadElement(); if (~verbose) print("Assemble Global Stiffness Matrix\t"); end @time K = Assem(); if (~verbose) print("Assemble Load Vector\t\t\t"); end @time UpdateLoad(); if (~verbose) print("Solve DOF\t\t\t\t"); end @time Solve(K); if (~verbose) print("Update Displacements if OFE is used\t") end @time UpdateNodes(); print("Total Time\t\t\t\t"); end #Read Line by Line function RL() global line += 1; return parse(Float64,text[line-1]); end #Gauss-Quadrature Functions function Legendre(n,x) #Algorithm taken from https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Legendre_quadrature if (n == 0) return 1; elseif (n == 1) return x; else return ((2*n-1)*x/n)*Legendre(n-1,x) - ((n-1)/n)*Legendre(n-2,x); end end function GaussPoints(n) points = zeros(1,n); #Algorithm taken from https://math.stackexchange.com/questions/12160/roots-of-legendre-polynomial if (n == 1) points = [0.0]; elseif (n == 2) points = [-0.5773502691896 0.5773502691896]; elseif (n == 3) points = [-0.7745966692415 0.0 0.7745966692415]; elseif (n == 4) points = [0.8611363115941 -0.3399810435849 0.3399810435849 0.8611363115941]; else for i = 1:n points[n+1-i] = (1 - (1/(8*(n^2))) + (1/(8*(n^3))))*cos(pi*(4*i-1)/(4*n + 2)); end end return points; end function GaussWeights(n) weights = zeros(1,n); x = GaussPoints(n); #Algorithm taken from http://mathworld.wolfram.com/Legendre-GaussQuadrature.html if (n == 1) weights = [2.0]; elseif (n == 2) weights = [1.0 1.0]; elseif (n == 3) weights = [5/9 8/9 5/9]; elseif (n == 4) weights = [0.3478548451375 0.6521451548625 0.6521451548625 0.3478548451375]; else for i = 1:n weights[i] = 2*(1 - (x[i]^2))/(((n+1)^2)*(Legendre(n+1,x[i])^2)); end end return weights; end function λ1(x,y) ans = (-x - (y-1)); return ans; end function λ2(x,y) ans = x; return ans; end function λ3(x,y) return 1 - λ1(x,y) - λ2(x,y); end function Barycentric(x,y) return [λ1(x,y), λ2(x,y), λ3(x,y)]; end function TriangleGaussPoints() #Data taken from https://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html GP = [0.33333333333333331 0.33333333333333331 2.06349616025259287E-002 0.48968251919873701 0.48968251919873701 2.06349616025259287E-002 0.48968251919873701 0.48968251919873701 0.12582081701412900 0.43708959149293553 0.43708959149293553 0.12582081701412900 0.43708959149293553 0.43708959149293553 0.62359292876193562 0.18820353561903219 0.18820353561903219 0.62359292876193562 0.18820353561903219 0.18820353561903219 0.91054097321109406 4.47295133944529688E-002 4.47295133944529688E-002 0.91054097321109406 4.47295133944529688E-002 4.47295133944529688E-002 0.74119859878449801 3.68384120547362581E-002 0.74119859878449801 0.22196298916076573 3.68384120547362581E-002 0.74119859878449801 3.68384120547362581E-002 0.22196298916076573 0.22196298916076573 0.74119859878449801 0.22196298916076573 3.68384120547362581E-002]; ans = zeros(length(GP[:,1]),3); for i = 1:length(GP[:,1]) ans[i,:] = Barycentric(GP[i,1],GP[i,2]); end return ans; end function TriangleGaussWeights() return [9.71357962827961025E-002 3.13347002271398278E-002 3.13347002271398278E-002 3.13347002271398278E-002 7.78275410047754301E-002 7.78275410047754301E-002 7.78275410047754301E-002 7.96477389272090969E-002 7.96477389272090969E-002 7.96477389272090969E-002 2.55776756586981006E-002 2.55776756586981006E-002 2.55776756586981006E-002 4.32835393772893970E-002 4.32835393772893970E-002 4.32835393772893970E-002 4.32835393772893970E-002 4.32835393772893970E-002 4.32835393772893970E-002]; end #Manager for Element types function Elemnt(Type, IND, K) if (Type == 1) Quads(IND, K); elseif (Type == 2) TriOFE2D(IND, K); else #Error end end function plotmesh() #global Mesh = Array{Float64,1}(undef,0); global plt = plot(0,0, legend = :none, line = :blue); for i = 1:NUMEG Elemnt(EG[i].Type, 4, ~); end plt end #Jacobian Solver for 2D Quadrilateral Elements function STDM(R, S, Ele) RP = 1+R; SP = 1+S; RM = 1-R; SM = 1-S; P = 0.25*[SP -SP -SM SM; RP RM -RM -RP]; QNodes = Ele.Node[:,1:2]; XJ = P*QNodes; M = inv(XJ)*P; B = zeros(3,8); for i = 1:8 j = Int(1.5 + 0.5*((-1)^i)); k = Int(ceil(i/2)); B[j,i] = M[j,k]; B[3,i] = M[(j%2)+1,k]; end return [B,det(XJ)]; end #FEM Solver for 2D Quadrilateral Elements function Quads(IND, K) if (IND == 1) if (Verbose) println("Quads"); end NUMMAT = Int(RL()); if (NUMMAT == 0) NUMMAT = 1; end NUME = Int(RL()); E = zeros(1,NUMMAT); THIC = zeros(1,NUMMAT); v = zeros(1,NUMMAT); resize!(Mat, length(Mat) + NUMMAT); for i = 1:NUMMAT N = Int(RL()); if (Verbose) println("\tMaterial ", N, "/", NUMMAT); end E[N] = RL(); if (Verbose) println("\t\tYoung's Modulus:\t", E[N]); end THIC[N] = RL(); if (Verbose) println("\t\tThickness:\t\t", THIC[N]); end v[N] = RL(); if (Verbose) println("\t\tPoisson's Ratio:\t", v[N]); end Mat[length(Mat)-NUMMAT+i] = Material(E[N],v[N],0,THIC[N]); end Ele = Array{Element,1}(undef, 0); resize!(Ele, length(Ele) + NUME); NN = vec(zeros(1,4)); for N = 1:NUME LM = Array{Int64,1}(undef, 0); qNodes = zeros(4,3); resize!(LM, 12); M = Int(RL()); if (Verbose) println("\tElement ", M, "/", NUME); end if (Verbose) print("\t\tNodes:\t"); end for i = 1:4 qN = Int(RL()); if (Verbose) print(qN, " "); end for j = 1:3 qNodes[i,j] = Node[qN,j]; LM[(i-1)*3 + j] = ID[qN,j]; end NN[i] = qN; FENodes[qN] = 1; end if (Verbose) println("\n\t\tLM:\t", LM); end MTYP = Int(RL()); if (Verbose) println("\t\tMaterial Type:\t", MTYP); end STYP = Int(RL()); if (Verbose) println("\t\tStress Type:\t", STYP); end NINT = Int(RL()); if (Verbose) println("\t\tGauss Order:\t", NINT); end Ele[N] = Element(STYP,qNodes,MTYP,NINT,LM,NN); end resize!(EG, length(EG) + 1); EG[length(EG)] = ElementGroup(NUME,1,Ele); elseif (IND == 2) findEG = 1; while (EG[findEG].Type != 1) findEG += 1; end EleGroup = EG[findEG]; if (usingOFE) global MeshStep = dist(EleGroup.Ele[1].Node[1,1:2], EleGroup.Ele[1].Node[2,1:2]); end for N = 1:EleGroup.NUME Ele = EleGroup.Ele[N]; MTYPE = Ele.Material; F = Mat[MTYPE].E/(1+Mat[MTYPE].v); G = F*Mat[MTYPE].v/(1-2*Mat[MTYPE].v); H = F+G; D = [H G 0; G H 0; 0 0 F/2]; if (Ele.Type == 1) XBAR = 1; elseif (Ele.Type == 2) D -= [(G^2)/H (G^2)/H 0; (G^2)/H (G^2)/H 0; 0 0 0]; XBAR = Mat[MTYPE].THIC; end #If not running OFE, use this if (!usingOFE) S = zeros(8,8); Pts = GaussPoints(Ele.NINT); Wts = GaussWeights(Ele.NINT); for LX = 1:Ele.NINT RI = Pts[LX]; for LY = 1:Ele.NINT SI = Pts[LY]; Res = STDM(RI,SI,Ele); B = Res[1]; DET = Res[2]; WT = Wts[LX]*Wts[LY]*XBAR*DET; S += B'*D*B*WT; end end LM = zeros(1,8); for i = 1:4 for j = 1:2 LM[(i-1)*2+j] = Ele.LM[(i-1)*3+j]; end end for i = 1:8 I = Int(LM[i]); for j = 1:8 J = Int(LM[j]); if (I != 0)&&(J != 0) K[I,J] += S[i,j]; end end end #Else run this code instead else S = zeros(8,8); Pts = GaussPoints(Ele.NINT); Wts = GaussWeights(Ele.NINT); for LX = 1:Ele.NINT RI = Pts[LX]; for LY = 1:Ele.NINT SI = Pts[LY]; Res = STDM(RI,SI,Ele); B = Res[1]; DET = Res[2]; WT = Wts[LX]*Wts[LY]*XBAR*DET; S += B'*D*B*WT; end end LM = zeros(1,8); for i = 1:4 for j = 1:2 LM[(i-1)*2+j] = OFEID[Ele.NN[i],(j-1)*NumMonomials+1]; end end for i = 1:8 I = Int(LM[i]); for j = 1:8 J = Int(LM[j]); if (I != 0)&&(J != 0) K[I,J] += S[i,j]; #Threads.atomic_add!(K[I,J], S[i,j]); end end end end end elseif (IND == 3) #Code for Stress Recovery elseif (IND == 4) findEG = 1; while (EG[findEG].Type != 1) findEG += 1; end EleGroup = EG[findEG]; for N = 1:EleGroup.NUME Ele = EleGroup.Ele[N]; xplt = [Ele.Node[:,1]; Ele.Node[1,1]]; yplt = [Ele.Node[:,2]; Ele.Node[1,2]]; plot!(plt, xplt, yplt, legend = :none, line = :blue); end elseif (IND == 5) findEG = 1; while (EG[findEG].Type != 1) findEG += 1; end EleGroup = EG[findEG]; for N = 1:EleGroup.NUME Ele = EleGroup.Ele[N]; xplt = zeros(5,1); yplt = zeros(5,1); for i = 1:5 xplt[i] = displacedNodes[Ele.NN[(i-1)%4+1],1]; yplt[i] = displacedNodes[Ele.NN[(i-1)%4+1],2]; end plot!(solnplt, xplt, yplt, legend = :none, line = :blue); end else #Error end end #2D TriOFE Functions function W(s::Float64) if (s >= 1) return 0; else return 1 - 6*s^2 + 8*s^3 - 3*s^4; end end function dist(I,J) return sqrt(sum((I-J).^2)); end function Φ(I,J,i,Nd1,Nd2,Nd3,Rad) S = zeros(3,6); ans = 0; alpha = 1; EP = zeros(3,2); EP[1,1] = (Nd1[1]+Nd2[1])/2; EP[1,2] = (Nd1[2]+Nd2[2])/2; EP[2,1] = (Nd2[1]+Nd3[1])/2; EP[2,2] = (Nd2[2]+Nd3[2])/2; EP[3,1] = (Nd3[1]+Nd1[1])/2; EP[3,2] = (Nd3[2]+Nd1[2])/2; if (i == J) ans = 1; elseif ((i <= 3)&&(J <= 3)) ans = 0; else S[1,4] = dist(EP[1,1:2],Nd1[1:2])/(alpha*Rad[1]); S[1,5] = dist(EP[2,1:2],Nd1[1:2])/(alpha*Rad[1]); S[1,6] = dist(EP[3,1:2],Nd1[1:2])/(alpha*Rad[1]); S[2,4] = dist(EP[1,1:2],Nd2[1:2])/(alpha*Rad[2]); S[2,5] = dist(EP[2,1:2],Nd2[1:2])/(alpha*Rad[2]); S[2,6] = dist(EP[3,1:2],Nd2[1:2])/(alpha*Rad[2]); S[3,4] = dist(EP[1,1:2],Nd3[1:2])/(alpha*Rad[3]); S[3,5] = dist(EP[2,1:2],Nd3[1:2])/(alpha*Rad[3]); S[3,6] = dist(EP[3,1:2],Nd3[1:2])/(alpha*Rad[3]); if (I == 1) if (i == 4) ans = W(S[J,i])/(W(S[1,i]) + W(S[2,i])); if (J == 3) ans = 0; end elseif (i == 5) ans = W(S[J,i])/(W(S[1,i]) + W(S[2,i]) + W(S[3,i])); elseif (i == 6) ans = W(S[J,i])/(W(S[1,i]) + W(S[3,i])); if (J == 2) ans = 0; end end elseif (I == 2) if (i == 4) ans = W(S[J,i])/(W(S[1,i]) + W(S[2,i])); if (J == 3) ans = 0; end elseif (i == 5) ans = W(S[J,i])/(W(S[2,i]) + W(S[3,i])); if (J == 1) ans = 0; end elseif (i == 6) ans = W(S[J,i])/(W(S[1,i]) + W(S[2,i]) + W(S[3,i])); end elseif (I == 3) if (i == 4) ans = W(S[J,i])/(W(S[1,i]) + W(S[2,i]) + W(S[3,i])); elseif (i == 5) ans = W(S[J,i])/(W(S[2,i]) + W(S[3,i])); if (J == 1) ans = 0; end elseif (i == 6) ans = W(S[J,i])/(W(S[1,i]) + W(S[3,i])); if (J == 2) ans = 0; end end end end return ans; end function A(Nodes) #Area Nds = [Nodes; Nodes[1,:]']; ans = 0; for i = 1:3 ans += 0.5*(Nds[i,1]*Nds[i+1,2] - Nds[i+1,1]*Nds[i,2]); end return ans; end function h(i, I, J, K, Nd, NINT) ans = 0; j = i%3 + 1; k = (i+1)%3 + 1; Nodes = [I'; J'; K']; Area = A(Nodes); #r = (j+1)%3 + 1; #s = (j+2)%3 + 1; s = j; r = k; ai = (Nodes[s,1]*Nodes[r,2] - Nodes[r,1]*Nodes[s,2])/(2*Area); bi = (Nodes[s,2] - Nodes[r,2])/(2*Area); ci = (Nodes[r,1] - Nodes[s,1])/(2*Area); Li = ai + bi*Nd[1] + ci*Nd[2]; Lj = 0; if (NINT == 1) ans = Li; elseif (NINT == 2) ans = bi; #This is for dh/dx (1st order shape function h) elseif (NINT == 3) ans = ci; #This is for dh/dy (1st order shape function h) elseif (NINT == 4) if (i <= 3) ans = Li*(2*Li - 1); else s = j%3 + 1; r = k%3 + 1; aj = (Nodes[s,1]*Nodes[r,2] - Nodes[r,1]*Nodes[s,2])/(2*Area); bj = (Nodes[s,2] - Nodes[r,2])/(2*Area); cj = (Nodes[r,1] - Nodes[s,1])/(2*Area); Lj = aj + bj*Nd[1] + cj*Nd[2]; ans = 4*Li*Lj; end elseif (NINT == 5) #This is for dh/dx (2nd order shape function h) if (i <= 3) ans = 4*bi*Li - bi; else s = j%3 + 1; r = k%3 + 1; aj = (Nodes[s,1]*Nodes[r,2] - Nodes[r,1]*Nodes[s,2])/(2*Area); bj = (Nodes[s,2] - Nodes[r,2])/(2*Area); cj = (Nodes[r,1] - Nodes[s,1])/(2*Area); Lj = aj + bj*Nd[1] + cj*Nd[2]; ans = 4*(bi*Lj + bj*Li); end elseif (NINT == 6) #This is for dh/dy (2nd order shape function h) if (i <= 3) ans = 4*ci*Li - ci; else s = j%3 + 1; r = k%3 + 1; aj = (Nodes[s,1]*Nodes[r,2] - Nodes[r,1]*Nodes[s,2])/(2*Area); bj = (Nodes[s,2] - Nodes[r,2])/(2*Area); cj = (Nodes[r,1] - Nodes[s,1])/(2*Area); Lj = aj + bj*Nd[1] + cj*Nd[2]; ans = 4*(ci*Lj + cj*Li); end end return ans; end function OFERadius(Nd, Connections) Radius = 0; for i = 1:length(Connections) if (Connections[i] != 0) if (Radius == 0) Radius = dist(Node[Nd,:], Node[i,:]); else if (dist(Node[Nd,:], Node[i,:]) > Radius) Radius = dist(Node[Nd,:], Node[i,:]); end end end end if (FENodes[Nd] > 1) Radius = max(Radius,sqrt(2)*MeshStep); end return Radius; end function OFE2DSTDM(Nodes, GP, Nd) B = zeros(3,6*NumMonomials); Rad = vec(zeros(1,3)); for i = 1:3 Rad[i] = OFERadius(Nd[i], OFEM[Nd[i],:]); end Pm = zeros(3,6); Pm[:,1] = ones(3,1); Pm[:,2] = (GP[1]*ones(3,1) - Nodes[:,1]); Pm[:,3] = (GP[2]*ones(3,1) - Nodes[:,2]); Pm[:,4] = ((GP[1]*ones(3,1) - Nodes[:,1]).*(GP[2]*ones(3,1) - Nodes[:,2])); Pm[:,5] = (GP[1]*ones(3,1) - Nodes[:,1]).*(GP[1]*ones(3,1) - Nodes[:,1]); Pm[:,6] = (GP[2]*ones(3,1) - Nodes[:,2]).*(GP[2]*ones(3,1) - Nodes[:,2]); dPmdx = zeros(3,6); dPmdy = zeros(3,6); for i = 1:3 dPmdx[i,:] = [0, 1, 0, GP[2]-Nodes[i,2], 2*(GP[1]-Nodes[i,1]), 0]; dPmdy[i,:] = [0, 0, 1, GP[1]-Nodes[i,1], 0, 2*(GP[2]-Nodes[i,2])]; end ρsum = 0; for I = 1:3 tm1 = 0; tm2x = 0; tm3x = 0; tm2y = 0; tm3y = 0; for k = 1:3 ϕIK = 0; dϕIKdx = 0; dϕIKdy = 0; for i = 1:6 hi = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,4); dhidx = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,5); dhidy = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,6); ΦIK = Φ(k,I,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕIK += hi*ΦIK; dϕIKdx += dhidx*ΦIK; dϕIKdy += dhidy*ΦIK; end hK = h(k,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhKdx = h(k,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhKdy = h(k,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); tm1 += hK*ϕIK; tm2x += dhKdx*ϕIK; tm3x += hK*dϕIKdx; tm2y += dhKdy*ϕIK; tm3y += hK*dϕIKdy; end dρImdx = zeros(4,1); dρImdy = zeros(4,1); ρI = tm1; dρIdx = tm2x + tm3x; dρIdy = tm2y + tm3y; ρsum += ρI; #At high alpha's, these should be equal #= println(ρI, "\t", h(I,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1)); println(dρIdx, "\t", h(I,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2)); println(dρIdy, "\t", h(I,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3)); =# #= println(abs(ρI - h(I,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1)) > 1e-8); println(abs(dρIdx - h(I,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2)) > 1e-8); println(abs(dρIdy - h(I,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3)) > 1e-8); =# dρImdx = dPmdx[I,:]*ρI + Pm[I,:]*dρIdx; dρImdy = dPmdy[I,:]*ρI + Pm[I,:]*dρIdy; for m = 1:NumMonomials B[1,(I-1)*2*NumMonomials+2*m-1] = dρImdx[m]; B[2,(I-1)*2*NumMonomials+2*m] = dρImdy[m]; B[3,(I-1)*2*NumMonomials+2*m-1] = dρImdy[m]; B[3,(I-1)*2*NumMonomials+2*m] = dρImdx[m]; end end return B; end function Coupling2DSTDM(Nodes, GP, Nd, Case) B = zeros(3,6*NumMonomials); Rad = vec(zeros(1,3)); for i = 1:3 Rad[i] = OFERadius(Nd[i], OFEM[Nd[i],:]); end #println(Rad, Nd) Pm = zeros(3,6); Pm[:,1] = ones(3,1); Pm[:,2] = (GP[1]*ones(3,1) - Nodes[:,1]); Pm[:,3] = (GP[2]*ones(3,1) - Nodes[:,2]); Pm[:,4] = ((GP[1]*ones(3,1) - Nodes[:,1]).*(GP[2]*ones(3,1) - Nodes[:,2])); Pm[:,5] = (GP[1]*ones(3,1) - Nodes[:,1]).*(GP[1]*ones(3,1) - Nodes[:,1]); Pm[:,6] = (GP[2]*ones(3,1) - Nodes[:,2]).*(GP[2]*ones(3,1) - Nodes[:,2]); dPmdx = zeros(3,6); dPmdy = zeros(3,6); for i = 1:3 dPmdx[i,:] = [0, 1, 0, GP[2]-Nodes[i,2], 2*(GP[1]-Nodes[i,1]), 0]; dPmdy[i,:] = [0, 0, 1, GP[1]-Nodes[i,1], 0, 2*(GP[2]-Nodes[i,2])]; end if (Case == 1) #I is Traditional, L and M are MFS ϕIL = 0; dϕILdx = 0; dϕILdy = 0; ϕIM = 0; dϕIMdx = 0; dϕIMdy = 0; ϕLL = 0; dϕLLdx = 0; dϕLLdy = 0; ϕLM = 0; dϕLMdx = 0; dϕLMdy = 0; ϕML = 0; dϕMLdx = 0; dϕMLdy = 0; ϕMM = 0; dϕMMdx = 0; dϕMMdy = 0; for i = 1:6 hi = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,4); dhidx = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,5); dhidy = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,6); ΦIL = Φ(2,1,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕIL += hi*ΦIL; dϕILdx += dhidx*ΦIL; dϕILdy += dhidy*ΦIL; ΦLL = Φ(2,2,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕLL += hi*ΦLL; dϕLLdx += dhidx*ΦLL; dϕLLdy += dhidy*ΦLL; ΦML = Φ(2,3,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕML += hi*ΦML; dϕMLdx += dhidx*ΦML; dϕMLdy += dhidy*ΦML; ΦIM = Φ(3,1,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕIM += hi*ΦIM; dϕIMdx += dhidx*ΦIM; dϕIMdy += dhidy*ΦIM; ΦLM = Φ(3,2,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕLM += hi*ΦLM; dϕLMdx += dhidx*ΦLM; dϕLMdy += dhidy*ΦLM; ΦMM = Φ(3,3,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕMM += hi*ΦMM; dϕMMdx += dhidx*ΦMM; dϕMMdy += dhidy*ΦMM; end hI = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhIdx = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhIdy = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); hL = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhLdx = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhLdy = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); hM = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhMdx = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhMdy = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); CM = hI + hL*ϕIL + hM*ϕIM; dCMdx = dhIdx + dhLdx*ϕIL + hL*dϕILdx + dhMdx*ϕIM + hM*dϕIMdx; dCMdy = dhIdy + dhLdy*ϕIL + hL*dϕILdy + dhMdy*ϕIM + hM*dϕIMdy; dρIdx = dhIdx*CM + hI*dCMdx; dρIdy = dhIdy*CM + hI*dCMdy; dρLdx = Pm[2,:]*(dhLdx*ϕLL + hL*dϕLLdx + dhMdx*ϕLM + hM*dϕLMdx) + dPmdx[2,:]*(hL*ϕLL + hM*ϕLM); dρLdy = Pm[2,:]*(dhLdy*ϕLL + hL*dϕLLdy + dhMdy*ϕLM + hM*dϕLMdy) + dPmdy[2,:]*(hL*ϕLL + hM*ϕLM); dρLdx[1] += dhLdx*CM + hL*dCMdx; dρLdy[1] += dhLdy*CM + hL*dCMdy; dρMdx = Pm[3,:]*(dhLdx*ϕML + hL*dϕMLdx + dhMdx*ϕMM + hM*dϕMMdx) + dPmdx[3,:]*(hM*ϕMM + hL*ϕML); dρMdy = Pm[3,:]*(dhLdy*ϕML + hL*dϕMLdy + dhMdy*ϕMM + hM*dϕMMdy) + dPmdy[3,:]*(hM*ϕMM + hL*ϕML); dρMdx[1] += dhMdx*CM + hM*dCMdx; dρMdy[1] += dhMdy*CM + hM*dCMdy; B[1,1] = dρIdx; B[2,2] = dρIdy; B[3,1] = dρIdy; B[3,2] = dρIdx; for i = 1:NumMonomials B[1,2*NumMonomials+2*i-1] = dρLdx[i]; B[2,2*NumMonomials+2*i] = dρLdy[i]; B[3,2*NumMonomials+2*i-1] = dρLdy[i]; B[3,2*NumMonomials+2*i] = dρLdx[i]; B[1,4*NumMonomials+2*i-1] = dρMdx[i]; B[2,4*NumMonomials+2*i] = dρMdy[i]; B[3,4*NumMonomials+2*i-1] = dρMdy[i]; B[3,4*NumMonomials+2*i] = dρMdx[i]; end elseif (Case == 2) #I is MFS, L and M are Traditional ϕLI = 0; ϕMI = 0; dϕLIdx = 0; dϕLIdy = 0; dϕMIdx = 0; dϕMIdy = 0; ϕII = 0; dϕIIdx = 0; dϕIIdy = 0; for i = 1:6 hi = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,4); dhidx = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,5); dhidy = h(i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,6); ΦII = Φ(1,1,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕII += hi*ΦII; dϕIIdx += dhidx*ΦII; dϕIIdy += dhidy*ΦII; ΦLI = Φ(1,2,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕLI += hi*ΦLI; dϕLIdx += dhidx*ΦLI; dϕLIdy += dhidy*ΦLI; ΦMI = Φ(1,3,i,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],Rad); ϕMI += hi*ΦMI; dϕMIdx += dhidx*ΦMI; dϕMIdy += dhidy*ΦMI; end hI = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhIdx = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhIdy = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); hL = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhLdx = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhLdy = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); hM = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,1); dhMdx = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhMdy = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); CM = hL + hM + hI*(ϕLI + ϕMI); dCMdx = dhLdx + dhMdx + dhIdx*(ϕLI + ϕMI) + hI*(dϕLIdx + dϕMIdx); dCMdy = dhLdy + dhMdy + dhIdy*(ϕLI + ϕMI) + hI*(dϕLIdy + dϕMIdy); dρIdx = Pm[1,:]*(dhIdx*ϕII + hI*dϕIIdx) + dPmdx[1,:]*hI*ϕII; dρIdy = Pm[1,:]*(dhIdy*ϕII + hI*dϕIIdy) + dPmdy[1,:]*hI*ϕII; dρIdx[1] += dhIdx*CM + hI*dCMdx; dρIdy[1] += dhIdy*CM + hI*dCMdy; dρLdx = dhLdx*CM + hL*dCMdx; dρLdy = dhLdy*CM + hL*dCMdy; dρMdx = dhMdx*CM + hM*dCMdx; dρMdy = dhMdy*CM + hM*dCMdy; for i = 1:NumMonomials B[1,2*i-1] = dρIdx[i]; B[2,2*i] = dρIdy[i]; B[3,2*i-1] = dρIdy[i]; B[3,2*i] = dρIdx[i]; end B[1,2*NumMonomials+1] = dρLdx; B[2,2*NumMonomials+2] = dρLdy; B[3,2*NumMonomials+1] = dρLdy; B[3,2*NumMonomials+2] = dρLdx; B[1,4*NumMonomials+1] = dρMdx; B[2,4*NumMonomials+2] = dρMdy; B[3,4*NumMonomials+1] = dρMdy; B[3,4*NumMonomials+2] = dρMdx; elseif (Case == 3) dhIdx = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhIdy = h(1,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); dhLdx = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhLdy = h(2,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); dhMdx = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,2); dhMdy = h(3,Nodes[1,1:2],Nodes[2,1:2],Nodes[3,1:2],GP,3); B[1,1] = dhIdx; B[2,2] = dhIdy; B[3,1] = dhIdy; B[3,2] = dhIdx; B[1,2*NumMonomials+1] = dhLdx; B[2,2*NumMonomials+2] = dhLdy; B[3,2*NumMonomials+1] = dhLdy; B[3,2*NumMonomials+2] = dhLdx; B[1,4*NumMonomials+1] = dhMdx; B[2,4*NumMonomials+2] = dhMdy; B[3,4*NumMonomials+1] = dhMdy; B[3,4*NumMonomials+2] = dhMdx; end return B; end #FEM Solver for TriOFE2D function TriOFE2D(IND, K) if (IND == 1) if (Verbose) println("TriOFE"); end NUMMAT = Int(RL()); if (NUMMAT == 0) NUMMAT = 1; end NUME = Int(RL()); E = zeros(1,NUMMAT); THIC = zeros(1,NUMMAT); v = zeros(1,NUMMAT); resize!(Mat, length(Mat) + NUMMAT); for i = 1:NUMMAT N = Int(RL()); if (Verbose) println("\tMaterial ", N, "/", NUMMAT); end E[N] = RL(); if (Verbose) println("\t\tYoung's Modulus:\t", E[N]); end THIC[N] = RL(); if (Verbose) println("\t\tThickness:\t\t", THIC[N]); end v[N] = RL(); if (Verbose) println("\t\tPoisson's Ratio:\t", v[N]); end Mat[length(Mat)-NUMMAT+i] = Material(E[N],v[N],0,THIC[N]); end Ele = Array{Element,1}(undef, 0); OFECircleNodes = Array{Int64,1}(undef, 0); OFECircles = Array{typeof(OFECircleNodes),1}(undef, 0); #Connected Nodes CI = vec(zeros(1,6*NUME)); CJ = vec(zeros(1,6*NUME)); CV = vec(zeros(1,6*NUME)); #Elements of a Node EI = Array{Int64,1}(undef, 0); EJ = Array{Int64,1}(undef, 0); EV = Array{Int64,1}(undef, 0); resize!(Ele, length(Ele) + NUME); for N = 1:NUME qNodes = zeros(3,3); LM = Array{Int64,1}(undef, 0); resize!(LM, 9); M = Int(RL()); if (Verbose) println("\tElement ", M, "/", NUME); end if (Verbose) print("\t\tNodes:\t"); end NN = vec(zeros(1,3)); for i = 1:3 qN = Int(RL()); if (Verbose) print(qN, " "); end for j = 1:3 qNodes[i,j] = Node[qN,j]; LM[(i-1)*3 + j] = ID[qN,j]; end CI[(N-1)*6 + 2*i-1] = qN; CI[(N-1)*6 + 2*i] = qN; CV[(N-1)*6 + 2*i-1] = 1; CV[(N-1)*6 + 2*i] = 1; resize!(EI, length(EI) + 1); resize!(EJ, length(EJ) + 1); resize!(EV, length(EV) + 1); EI[length(EI)] = qN; EJ[length(EI)] = N; EV[length(EI)] = 1; NN[i] = qN; end n1 = 2; n2 = 3; for i = 1:3 CJ[(N-1)*6 + 2*i-1] = CI[(N-1)*6 + 2*n1]; CJ[(N-1)*6 + 2*i] = CI[(N-1)*6 + 2*n2]; n1 = n1%3 + 1; n2 = n2%3 + 1; end MTYP = Int(RL()); if (Verbose) println("\n\t\tMaterial Type:\t", MTYP); end STYP = Int(RL()); if (Verbose) println("\t\tStress Type:\t", STYP); end NINT = Int(RL()); if (Verbose) println("\t\tGauss Order:\t", NINT); end Ele[N] = Element(STYP,qNodes,MTYP,NINT,LM,NN); end resize!(EG, length(EG) + 1); EG[length(EG)] = ElementGroup(NUME,2,Ele); global OFEM = sparse([CI; NUMNP],[CJ; NUMNP],[CV; 0]); #Which OFE nodes are connected to each other global OFEID = [ID zeros(length(ID[:,1]),2*NumMonomials-3)]; #DOF of all nodes expanded into OFE basis global OFENodes = Array{Int64,1}(undef, 0); #Collection of Nodes that are part of OFE global OFEEle = sparse(EI,EJ,EV); # OFEIDCount = 0; for i = 1:length(OFEID[:,1]) temp = OFEID[i,2]; OFEID[i,2] = 0; OFEID[i,NumMonomials+1] = temp; if (sum(OFEM[i,:]) > 0) if (FENodes[i] > 0) #OFEM[i,:] = zeros(1,NUMNP); #OFEM[:,i] = zeros(1,NUMNP); if (OFEID[i,1] > 0) OFEIDCount += 1; OFEID[i,1] = OFEIDCount; end if (OFEID[i,NumMonomials+1] > 0) OFEIDCount += 1; OFEID[i,NumMonomials+1] = OFEIDCount; end else if (OFEID[i,1] > 0) for j = 1:NumMonomials OFEIDCount += 1; OFEID[i,j] = OFEIDCount; end else for j = 2:NumMonomials OFEIDCount += 1; OFEID[i,j] = OFEIDCount; end end if (OFEID[i,NumMonomials+1] > 0) for j = NumMonomials+1:2*NumMonomials OFEIDCount += 1; OFEID[i,j] = OFEIDCount; end else for j = NumMonomials+2:2*NumMonomials OFEIDCount += 1; OFEID[i,j] = OFEIDCount; end end resize!(OFENodes, length(OFENodes) + 1); OFENodes[length(OFENodes)] = i; end else if (OFEID[i,1] > 0) OFEIDCount += 1; OFEID[i,1] = OFEIDCount; end if (OFEID[i,NumMonomials+1] > 0) OFEIDCount += 1; OFEID[i,NumMonomials+1] = OFEIDCount; end end end if (OFEIDCount > NEQ) global usingOFE; usingOFE = true; end global NEQ = OFEIDCount; elseif (IND == 2) global newR = zeros(NEQ,1); findEG = 1; while (EG[findEG].Type != 2) findEG += 1; end EleGroup = EG[findEG]; LMArray = zeros(EleGroup.NUME,6*NumMonomials); SArray = zeros(EleGroup.NUME,6*NumMonomials,6*NumMonomials); for N = 1:EleGroup.NUME Ele = EleGroup.Ele[N]; MTYPE = EleGroup.Ele[N].Material; E = Mat[MTYPE].E; v = Mat[MTYPE].v; NINT = Ele.NINT; if (Ele.Type == 1) XBAR = 1; else XBAR = Mat[MTYPE].THIC; end F = E/(1+v); G = F*v/(1-2*v); H = F+G; D = [H G 0; G H 0; 0 0 F/2]; if (Ele.Type == 2) D -= [(G^2)/H (G^2)/H 0; (G^2)/H (G^2)/H 0; 0 0 0]; end Nd = Ele.NN[1]; Nd2 = Ele.NN[2]; Nd3 = Ele.NN[3]; S = zeros(6*NumMonomials,6*NumMonomials); Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); #Coordinates Area = A(Nds); if (Area < 0) Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd3,:]) vec(Node[Nd2,:])]'); temp = Nd2; Nd2 = Nd3; Nd3 = temp; Area = -Area; end countTraditional = 0; bin = zeros(3,1); if (FENodes[Nd] > 0) countTraditional += 1; bin[1] = 1; end if (FENodes[Nd2] > 0) countTraditional += 1; bin[2] = 1; end if (FENodes[Nd3] > 0) countTraditional += 1; bin[3] = 1; end #println(bin, " ", countTraditional) if (countTraditional == 1) if (bin[1] == 1) elseif (bin[2] == 1) temp = Nd; Nd = Nd2; Nd2 = Nd3; Nd3 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); else temp = Nd; Nd = Nd3; Nd3 = Nd2; Nd2 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); end elseif (countTraditional == 2) if (bin[1] == 0) elseif (bin[2] == 0) temp = Nd; Nd = Nd2; Nd2 = Nd3; Nd3 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); else temp = Nd; Nd = Nd3; Nd3 = Nd2; Nd2 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); end end #Gauss Points #= GP = zeros(12,length(Node[Nd,:])); TriWt = [0.050844906370207*ones(3,1); 0.116786275726379*ones(3,1); 0.082851075618374*ones(6,1)]; AreaCoord = zeros(12,3); AreaCoord[1,:] = [0.873821971016996 0.063089014491502 0.063089014491502]; AreaCoord[2,:] = [0.063089014491502 0.873821971016996 0.063089014491502]; AreaCoord[3,:] = [0.063089014491502 0.063089014491502 0.873821971016996]; AreaCoord[4,:] = [0.501426509658179 0.249286745170911 0.249286745170911]; AreaCoord[5,:] = [0.249286745170911 0.501426509658179 0.249286745170911]; AreaCoord[6,:] = [0.249286745170911 0.249286745170911 0.501426509658179]; AreaCoord[7,:] = [0.636502499121399 0.310352451033785 0.053145049844816]; AreaCoord[8,:] = [0.636502499121399 0.053145049844816 0.310352451033785]; AreaCoord[9,:] = [0.310352451033785 0.636502499121399 0.053145049844816]; AreaCoord[10,:] = [0.310352451033785 0.053145049844816 0.636502499121399]; AreaCoord[11,:] = [0.053145049844816 0.636502499121399 0.310352451033785]; AreaCoord[12,:] = [0.053145049844816 0.310352451033785 0.636502499121399]; =# #= GP = zeros(13,length(Node[Nd,:])); TriWt = [-0.149570044467670; 0.175615257433204*ones(3,1); 0.053347235608839*ones(3,1); 0.077113760890257*ones(6,1)]; AreaCoord = zeros(length(GP[:,1]),3); AreaCoord[1,:] = [1/3 1/3 1/3]; AreaCoord[2,:] = [0.479308067841923 0.260345966079038 0.260345966079038]; AreaCoord[3,:] = [0.260345966079038 0.479308067841923 0.260345966079038]; AreaCoord[4,:] = [0.260345966079038 0.260345966079038 0.479308067841923]; AreaCoord[5,:] = [0.869739794195568 0.065130102902216 0.065130102902216]; AreaCoord[6,:] = [0.065130102902216 0.869739794195568 0.065130102902216]; AreaCoord[7,:] = [0.065130102902216 0.065130102902216 0.869739794195568]; AreaCoord[8,:] = [0.638444188569809 0.312865496004875 0.048690315425316]; AreaCoord[9,:] = [0.638444188569809 0.048690315425316 0.312865496004875]; AreaCoord[10,:] = [0.312865496004875 0.638444188569809 0.048690315425316]; AreaCoord[11,:] = [0.312865496004875 0.048690315425316 0.638444188569809]; AreaCoord[12,:] = [0.048690315425316 0.638444188569809 0.312865496004875]; AreaCoord[13,:] = [0.048690315425316 0.312865496004875 0.638444188569809]; =# AreaCoord = TriangleGaussPoints(); TriWt = TriangleGaussWeights(); GP = zeros(length(AreaCoord[:,1]),length(Node[Nd,:])); for i = 1:length(GP[:,1]) GP[i,1] = sum(Nds[:,1].*AreaCoord[i,:]); GP[i,2] = sum(Nds[:,2].*AreaCoord[i,:]); if (countTraditional == 0) B = OFE2DSTDM(Nds, GP[i,1:2], [Nd Nd2 Nd3]); else B = Coupling2DSTDM(Nds, GP[i,1:2], [Nd Nd2 Nd3], countTraditional); end S += TriWt[i]*B'*D*B*Area*XBAR; end LM = []; for i = 1:NumMonomials LM = [LM; OFEID[Nd,i]; OFEID[Nd,NumMonomials+i]]; end for i = 1:NumMonomials LM = [LM; OFEID[Nd2,i]; OFEID[Nd2,NumMonomials+i]]; end for i = 1:NumMonomials LM = [LM; OFEID[Nd3,i]; OFEID[Nd3,NumMonomials+i]]; end LMArray[N,:] = Int.(LM); SArray[N,:,:] = S; end for N = 1:EleGroup.NUME for i = 1:6*NumMonomials I = Int(LMArray[N,i]); for j = 1:6*NumMonomials J = Int(LMArray[N,j]); if ((I != 0)&&(J != 0)) K[I,J] += SArray[N,i,j]; end end end end elseif (IND == 3) #Stress Recovery (Work in Progress) findEG = 1; while (EG[findEG].Type != 2) findEG += 1; end EleGroup = EG[findEG]; LMArray = zeros(EleGroup.NUME,6*NumMonomials); SArray = zeros(EleGroup.NUME,6*NumMonomials,6*NumMonomials); N = 0; for CompareStress = 1:sum(OFEEle[EvalNode,:]) N += 1; while (OFEEle[EvalNode,N] == 0) N += 1; end Ele = EleGroup.Ele[N]; MTYPE = EleGroup.Ele[N].Material; E = Mat[MTYPE].E; v = Mat[MTYPE].v; NINT = Ele.NINT; if (Ele.Type == 1) XBAR = 1; else XBAR = Mat[MTYPE].THIC; end F = E/(1+v); G = F*v/(1-2*v); H = F+G; D = [H G 0; G H 0; 0 0 F/2]; if (Ele.Type == 2) D -= [(G^2)/H (G^2)/H 0; (G^2)/H (G^2)/H 0; 0 0 0]; end Nd = Ele.NN[1]; Nd2 = Ele.NN[2]; Nd3 = Ele.NN[3]; σ = zeros(3,1); Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); #Coordinates Area = A(Nds); if (Area < 0) Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd3,:]) vec(Node[Nd2,:])]'); temp = Nd2; Nd2 = Nd3; Nd3 = temp; Area = -Area; end countTraditional = 0; bin = zeros(3,1); if (FENodes[Nd] > 0) countTraditional += 1; bin[1] = 1; end if (FENodes[Nd2] > 0) countTraditional += 1; bin[2] = 1; end if (FENodes[Nd3] > 0) countTraditional += 1; bin[3] = 1; end if (countTraditional == 1) if (bin[1] == 1) elseif (bin[2] == 1) temp = Nd; Nd = Nd2; Nd2 = Nd3; Nd3 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); else temp = Nd; Nd = Nd3; Nd3 = Nd2; Nd2 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); end elseif (countTraditional == 2) if (bin[1] == 0) elseif (bin[2] == 0) temp = Nd; Nd = Nd2; Nd2 = Nd3; Nd3 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); else temp = Nd; Nd = Nd3; Nd3 = Nd2; Nd2 = temp; Nds = Matrix([vec(Node[Nd,:]) vec(Node[Nd2,:]) vec(Node[Nd3,:])]'); end end EvalCoord = Node[EvalNode,:]; if (countTraditional == 0) B = OFE2DSTDM(Nds, EvalCoord[1:2], [Nd Nd2 Nd3]); else B = Coupling2DSTDM(Nds, EvalCoord[1:2], [Nd Nd2 Nd3], countTraditional); end displacement = zeros(6*NumMonomials,1); for i = 1:NumMonomials if (OFEID[Nd,i] == 0) displacement[2*i-1] = 0; else displacement[2*i-1] = α[Int(OFEID[Nd,i])]; end end for i = 1:NumMonomials if (OFEID[Nd,NumMonomials+i] == 0) displacement[2*i] = 0; else displacement[2*i] = α[Int(OFEID[Nd,NumMonomials+i])]; end end for i = 1:NumMonomials if (OFEID[Nd2,i] == 0) displacement[2*NumMonomials+2*i-1] = 0; else displacement[2*NumMonomials+2*i-1] = α[Int(OFEID[Nd2,i])]; end end for i = 1:NumMonomials if (OFEID[Nd2,NumMonomials+i] == 0) displacement[2*NumMonomials+2*i] = 0; else displacement[2*NumMonomials+2*i] = α[Int(OFEID[Nd2,NumMonomials+i])]; end end for i = 1:NumMonomials if (OFEID[Nd3,i] == 0) displacement[4*NumMonomials+2*i-1] = 0; else displacement[4*NumMonomials+2*i-1] = α[Int(OFEID[Nd3,i])]; end end for i = 1:NumMonomials if (OFEID[Nd3,NumMonomials+i] == 0) displacement[4*NumMonomials+2*i] = 0; else displacement[4*NumMonomials+2*i] = α[Int(OFEID[Nd3,NumMonomials+i])]; end end σ = D*B*displacement; println(σ); end elseif (IND == 4) findEG = 1; while (EG[findEG].Type != 2) findEG += 1; end EleGroup = EG[findEG]; for N = 1:EleGroup.NUME Ele = EleGroup.Ele[N]; xplt = [Ele.Node[:,1]; Ele.Node[1,1]]; yplt = [Ele.Node[:,2]; Ele.Node[1,2]]; plot!(plt, xplt, yplt, legend = :none, line = :blue); end elseif (IND == 5) findEG = 1; while (EG[findEG].Type != 2) findEG += 1; end EleGroup = EG[findEG]; for N = 1:EleGroup.NUME Ele = EleGroup.Ele[N]; xplt = zeros(4,1); yplt = zeros(4,1); for i = 1:4 xplt[i] = displacedNodes[Ele.NN[(i-1)%3+1],1]; yplt[i] = displacedNodes[Ele.NN[(i-1)%3+1],2]; end plot!(solnplt, xplt, yplt, legend = :none, line = :blue); end else #Error end end #Read Control Information function RCI() global line = 1; global NUMNP = Int(RL()); if (Verbose) println("Number of Nodal Points:\t\t", NUMNP); end global NUMEG = Int(RL()); if (Verbose) println("Number of Element Groups:\t", NUMEG); end global NLCASE = Int(RL()); if (Verbose) println("Number of Load Cases:\t\t", NLCASE); end global MODEX = Int(RL()); if (Verbose) println("Mode of Execution:\t\t", MODEX); end end #Read Input function ReadInput() global line = 5; global ID = zeros(Int64,NUMNP,3); global Node = zeros(NUMNP,3); if (TradMode) global FENodes = ones(NUMNP,1); else global FENodes = zeros(NUMNP,1); #FENodes[2:4] = ones(3,1); end N = 0; while (N != NUMNP) N = Int(RL()); for i = 1:3 ID[N,i] = Int(RL()); end for i = 1:3 Node[N,i] = RL(); end KN = RL(); if (Verbose) println("Node ", N); end if (Verbose) println("BC:\t\t", ID[N,1], "\t", ID[N,2], "\t", ID[N,3]);end if (Verbose) println("Coordinates:\t", Node[N,1], "\t", Node[N,2], "\t", Node[N,3]);end end end #Identifying Degrees of Freedom function IDDOF() global NEQ = 0; for N = 1:NUMNP for i = 1:3 if (ID[N,i] == 0) NEQ += 1; ID[N,i] = NEQ; else ID[N,i] = 0; end end end end #Calculate and Store Load Vectors function Loads() global line = 5 + 8*NUMNP; global R = zeros(NEQ,1); for L = 1:NLCASE LL = Int(RL()); NLOAD = Int(RL()); if (Verbose) println("Load Case ", L, "/", NLCASE); end for l = 1:NLOAD if (Verbose) println("\tLoad ", l ,"/", NLOAD); end NOD = Int(RL()); if (Verbose) println("\t\tNode:\t\t\t", NOD); end IDIRN = Int(RL()); if (Verbose) println("\t\tForce Direction:\t", (IDIRN == 1) ? "x" : ((IDIRN == 2) ? "y" : "z")); end FLOAD = RL(); if (Verbose) println("\t\tForce:\t\t\t", FLOAD, "N"); end R[ID[NOD,IDIRN]] += FLOAD; end end end #Read Element Data function ReadElement() global EG = Array{ElementGroup,1}(undef,0); global Mat = Array{Material,1}(undef,0); #logg = [] for i = 1:NUMEG Type = Int(RL()); if (Verbose) print("Element ", i, "/", NUMEG, ":\t"); end x = Elemnt(Type, 1, ~); #push!(logg,x) end #logg end #Assemble Stiffness Matrix function Assem() logg = [] K = zeros(NEQ,NEQ); for i = 1:NUMEG Elemnt(EG[i].Type, 2, K); #push!(logg, Elemnt(EG[i].Type, 2, K)); end return sparse(K); #return sparse(K), logg; end #Update Load Vector function UpdateLoad() if (usingOFE) for i = 1:length(R) if (R[i] != 0) j = 0; l = 1; findNode = false; while ((!findNode)&&(j <= NUMNP)) j += 1; for k = 1:2 if (ID[j,k] == i) findNode = true; l = k; break; end end end if (j > NUMNP) println("Unknown Error Encountered") break; end newR[Int(OFEID[j,(l-1)*NumMonomials+1])] = R[i]; end end end #newR end function Solve(K) global α = 0; global u = zeros(NUMNP,2); global Stiffness = K; if (usingOFE) if (UseGPU) α = Array(AFArray(Matrix(K))\AFArray(newR)); #α = Array(AFArray(K)\AFArray(newR)); #ArrayFire does not work with Sparse Matrices yet else α = K\newR; end for i = 1:NUMNP for j = 1:2 if (OFEID[i,(j-1)*NumMonomials + 1] > 0) u[i,j] = α[Int(OFEID[i,(j-1)*NumMonomials + 1])]; end end end else if (UseGPU) α = Array(AFArray(Matrix(K))\AFArray(R)); else α = K\R; end for i = 1:NUMNP for j = 1:2 if (ID[i,j] > 0) u[i,j] = α[Int(ID[i,j])]; end end end end end function UpdateNodes() global displacedNodes = Node[:,1:length(u[1,:])]; for i = 1:NUMNP displacedNodes[i,:] += u[i,:]; end end function plotsolution() global solnplt = plot(0,0, legend = :none, line = :blue); for i = 1:NUMEG Elemnt(EG[i].Type, 5, ~) end solnplt end #Input file has to be consistent with the STAP input procedure #tradmode = true sets the solver to use only traditional Finite Elements @time OFESolve("PlateWithHole.txt", verbose = false, tradmode = false, useGPU = false) #Original Geometry plotmesh() #Deformed Geometry plotsolution() #Evaluate Stresses at an OFE Node (Work in Progress) #Returns [σ_xx, σ_yy, σ_xy] once for each triangular element the node is part of EvalNode = 6464 TriOFE2D(3,Stiffness) #Global Stiffness Matrix Matrix(Stiffness) #Force vectors for the various DOFs #= for i = 1:NEQ if (i < 100) print("DOF ", i , ":\t\t", newR[i], "N\n"); else print("DOF ", i , ":\t", newR[i], "N\n"); end end =# #DOF indexing for each node #= for i = 1:NUMNP if (i < 10) print("Node ", i, ":\t\t"); else print("Node ", i, ":\t"); end for j = 1:2*NumMonomials print(Int(OFEID[i,j]), "\t") end print("\n"); end =# #Nodal Displacements [["Node";range(1,stop=length(u[:,1]),length(u[:,1]))] ["u_x" "u_y";u]]