structopt.jl

    if run_own
        #update_area(elements, A0)
        λ = [1.0]
        for k = 1:n_iters
            update_limits!(L, U, α, β, μ, s, k, xmax, xmin, x_k, x_k1, x_k2)
            L = zeros(length(elements))
            ∇g0 = zeros(n_eles)
            g0 = objective(x, ∇g0)
            r0, q0 = compute_mma(elements, x, u_full, g0, ∇g0, L)
            φ = (λ) -> compute_φ(λ, Float64[], x, elements, r0, Vmax, q0, L, α, β)
            φgra = (λ, gra) -> compute_φ(λ, gra, x, elements, r0, Vmax, q0, L, α, β)
            d = DifferentiableFunction(φ, φgra, φgra)
            results = fminbox(d, λ, [0.0], [Inf], xtol=1e-9, ftol=0.0)
            λ = results.minimum
        end
        figure()
        plot_deformed(elements, u_full, x, 0.1)
        print_stuff(x)

    #for (i, ele) in enumerate(elements)
    #    println(stress(ele, x[i], u_full))
    #end

    end



plot_undeformed(elements)


    return
end




function steepest_ascent(g, λ)
    λs = Float64[]
    fs = Float64[]
    rel_tol = 1e-9
    ∇f = [0.0]

    max_iters = 50
    c = 0.5
    τ = 0.5
    γ = 0.5
    done = false
    i = 0
    for i = 1:max_iters
        f = g(λ, ∇f)
        # Line search
        n_ls = 1
        # Armijo - Goldstein
        while true

            γ *= 1.5
            debug("γ = $γ")
            n_ls += 1
            λ_new = λ - γ * ∇f
            f_new = g(λ_new, Float64[])
            #return
            if -f_new > -f + c * γ * norm(∇f)^2
                #λ += γ * ∇f
                λ = λ_new
             #   push!(λs, λ)
             #   push!(fs, f_new)
                break
            else
                γ *= 0.5
            end

            #debug("Updating γ from $(2*γ) to $γ")
            if (norm(λ_new - λ) / norm(λ) ) < rel_tol
                #debug("n_iters = $i, λ = $λ_new, f = $f_new")
                done = true
                break
            end
        end
        if done
            break
        end
    end
   # plot(λs, fs, "*-")
    return λ

end


function update_limits!(L, U, α, β, μ, s, k, xmax, xmin, x_k2, x_k1, x_k)
    for j = 1:length(x_k)
        if k == 1 || k == 2
            L[j] = x_k[j] - (xmax[j] - xmin[j])
            U[j] = x_k[j] + (xmax[j] - xmin[j])
        else
            if sign(x_k[j] - x_k1[j]) != sign(x_k1[j] - x_k2[j])
                L[j] = x_k[j] - s * (x_k1[j] - L[j])
                U[j] = x_k[j] + s * (U[j] - x_k1[j])
            else
                L[j] = x_k[j] - (x_k1[j] - L[j]) / s
                U[j] = x_k[j] + (U[j] - x_k1[j]) / s
            end
        end
        #α[j] = max(xmin[j], L[j] + μ * (x_k[j] - L[j]))
        #β[j] = min(xmax[j], U[j] - μ * (U[j] - x_k[j]))
        α[j] = max(0.9 * L[j] + 0.1 * x_k1[j], xmin[j])
        β[j] = min(0.9 * U[j] + 0.1 * x_k1[j], xmax[j])
    end
end

function compute_areas!(λ, x, bars, L, q0, α, β)
    for (j, element) in enumerate(bars)
        lj = length(element)
        x_opt = L[j] + sqrt(q0[j] / (λ*lj))
      #  println("j: $j")
      #  println("xopt: $x_opt")
        if x_opt < α[j]
           x[j] = α[j]
        elseif x_opt > β[j]
            x[j] = β[j]
        else
            x[j] = x_opt
        end
    end
end

function compute_φ(λ, δφ_δλ, x, bars, r0, Vmax, q0, L, α, β)
    if λ[1] < 0.0 # Bro you don't wanna be here...
        if length(δφ_δλ) > 0
            δφ_δλ[1] = 10e8
        end
        return 10e9
    end
    φ = r0 - λ[1] * Vmax
    if length(δφ_δλ) > 0
        δφ_δλ[1] = - Vmax
    end
    compute_areas!(λ[1], x, bars, L, q0, α, β)
  #  println("x_nub: $x")
  #  println("q0_nub: $q0")
    for (j, element) in enumerate(bars)
        lj = length(element)
        φ += q0[j] / (x[j] - L[j]) + λ[1]*lj*x[j]
        if length(δφ_δλ) > 0
            δφ_δλ[1] += lj * x[j]
        end
    end
   scale!(δφ_δλ, -1.0)
    return -φ
end


function compute_mma(elements::Vector{Bar}, x, u, g0, ∇g0, L)
    q0 = zeros(length(elements))
    r0 = g0
    for (j, element) in enumerate(elements)
        q0[j] = (x[j] - L[j])^2 * (-∇g0[j])
        r0 -= (x[j] - L[j]) * (-∇g0[j])
    end
    return r0, q0
end