topologyOptimization.jl

                        append!(izg,k)
                        append!(iyg,j)
                        append!(izg,k)
                    end
                    if(I[j,i,k,3]>threshold)
                        append!(ixb,nelx-i+nelx-1)
                        append!(ixb,i)
                        append!(iyb,j)
                        append!(izb,k)
                        append!(iyb,j)
                        append!(izb,k)
                        
                    end
                end
            end
        end
    end
    p=scatter(ixr,izr,iyr,color="red",label="",markersize =4, aspect_ratio=:equal,markerstrokealpha = 0.2,markeralpha = 0.6,markershape = :square,camera = (30, 60),xlim=(0,nelx),zlim=(0,nely),ylim=(-10,10))#,markershape = :square
    p=scatter!(ixg,izg,iyg,color="green",label="",markersize =4, aspect_ratio=:equal,markerstrokealpha = 0.2,markeralpha = 0.6,markershape = :square,camera = (30, 60),xlim=(0,nelx),zlim=(0,nely),ylim=(-10,10))#,markershape = :square
    p=scatter!(ixb,izb,iyb,color="blue",label="",markersize =4, aspect_ratio=:equal,markerstrokealpha = 0.2,markeralpha = 0.6,markershape = :square,camera = (30, 60),xlim=(0,nelx),zlim=(0,nely),ylim=(-10,10))#,markershape = :square
    return p

end

function plotBoundaryConditions3D(nelx,nely,nelz,U,F,freedofs,z=1)
    display("Supports")
    UU=U[:,1]
    UU[freedofs].=1.0
    UUU=reshape(UU,3,nely+1,nelx+1,nelz+1)[1,:,:,z]
    display(heatmap(UUU,aspect_ratio=:equal))
    UUU=reshape(UU,3,nely+1,nelx+1,nelz+1)[2,:,:,z]
    display(heatmap(UUU,aspect_ratio=:equal))
    UUU=reshape(UU,3,nely+1,nelx+1,nelz+1)[3,:,:,z]
    display(heatmap(UUU,aspect_ratio=:equal))


    UU=U[:,2]
    UU[freedofs].=1.0
    UUU=reshape(UU,3,nely+1,nelx+1,nelz+1)[1,:,:,z]
    display(heatmap(UUU,aspect_ratio=:equal))
    UUU=reshape(UU,3,nely+1,nelx+1,nelz+1)[2,:,:,z]
    display(heatmap(UUU,aspect_ratio=:equal))
    UUU=reshape(UU,3,nely+1,nelx+1,nelz+1)[3,:,:,z]
    display(heatmap(UUU,aspect_ratio=:equal))

    display("Forces")
    display("din")
    FF=F[:,1]
    FFF=Array(reshape(FF,3,nely+1,nelx+1,nelz+1)[1,:,:,z])
    display(heatmap(FFF,aspect_ratio=:equal))
    FFF=Array(reshape(FF,3,nely+1,nelx+1,nelz+1)[2,:,:,z])
    display(heatmap(FFF,aspect_ratio=:equal))
    FFF=Array(reshape(FF,3,nely+1,nelx+1,nelz+1)[3,:,:,z])
    display(heatmap(FFF,aspect_ratio=:equal))

    display("dout")
    FF=F[:,2]
    FFF=Array(reshape(FF,3,nely+1,nelx+1,nelz+1)[1,:,:,z])
    display(heatmap(FFF,aspect_ratio=:equal))
    FFF=Array(reshape(FF,3,nely+1,nelx+1,nelz+1)[2,:,:,z])
    display(heatmap(FFF,aspect_ratio=:equal))
    FFF=Array(reshape(FF,3,nely+1,nelx+1,nelz+1)[3,:,:,z])
    display(heatmap(FFF,aspect_ratio=:equal))
end
#####################################Utils####################################################################
function getIndex(i,j,nelx,nely)
    return (i-1)*(nely+1)+(j-1)+1
end

# get index for 3d compliant
function getIndex3d(i,j,k,nelx,nely,nelz)
    return Int((nelx*nely)*(k-1) + nely*(i-1) + j);
end

# get index for 3d compliant
function getIndex3d1(i,j,k,dof,nelx,nely,nelz,ndof)
    return Int((ndof*nelx*nely)*(k-1) + ndof*nely*(i-1) + ndof*(j-1)+dof);
end

function getIndex3d(i,j,k,dof,nelx,nely,nelz,ndof)
    return Int.((ndof.*nelx.*nely).*(k.-1.0) .+ ndof.*nely.*(i.-1.0) .+ ndof.*(j.-1.0).+dof);
end