component RelaxationTest
export Executable

run(args:String...):() = do
    (* The purpose of the main program is to test the constitutive driver
       in a relaxation test. The output will be in MATLAB syntax for easy
       plotting with octave/gnuplot. *)

    println("% Relaxation of viscoelastic material" )
    (* Constants *)
    E       :RR64 = 210.0 TIMES 10^9
    n       :RR64 = 3.0
    t_star  :RR64 = 1000.0
    t_end   :RR64 = 1000000.0
    m       :ZZ32 = 100

    (* Variables *)
    var sigma :Array[\RR64,ZZ32\] := array[\RR64\](m+1)
    var t     :Array[\RR64,ZZ32\] := array[\RR64\](m+1)
    var E_tang:RR64

    (* First load is applied very quickly, so the response is fully elastic *)
    sigma[0] := 0.1 E
    t[0]     := 0.0
    (* Now calculate the stress response by keeping the strain rate at zero *)
    println "data = ["
    for i <- seq(0#m) do
        (sigma[i+1],E_tang) := ViscoElastic1D(sigma[i], 0.0, t_end/m, E, n, t_star)
        t[i+1] := t[i]+t_end/m
        println t[i] ", " sigma[i] ";"
    end
    println t[m] ", " sigma[m] "];"
    println "plot(data(:,1),data(:,2))"
end

(* Strain driven constitutive driver *)
ViscoElastic1D(sigma_old, DELTA_epsilon, DELTA_t, E, n, t_star):(Number,Number) 
    requires{ DELTA_t > 0, t_star > 0, n >= 1, E > 0 } = do

    var sigma :Number
    var E_tang:Number

    if (n = 1) then
        sigma  := (E DELTA_epsilon + sigma_old/DELTA_t) / (1/DELTA_t + 1/t_star)
        E_tang := E
    else
        iter_max  = 20
        tolerance = 10^(-3)
        temp1     = DELTA_t/(t_star E^(n-1))
        temp2     = E DELTA_epsilon
        (* There is no point in making this look parallell *)
        sigma         := sigma_old
        var h :Number := 2.0 tolerance
        var i :ZZ     := 0
        while (i < iter_max AND abs(h) > tolerance sigma ) do
            i     += 1
            temp3  = temp1 (abs(sigma))^(n-1)
            f      = sigma - sigma_old + temp3 sigma - temp2;
            f'     = 1 + n temp3
            h     := f/f'
            sigma -= h
        end
        E_tang := E/(1 + n DELTA_t/t_star (abs(sigma)/E)^(n-1))
    end
    (sigma,E_tang)
end

(* I believe this should already exist, as |x| *)
abs(x:Number):(Number) = if (x >= 0) then x else -x end

end