# diffusion.tcl --
#     Solve a two-dimensional diffusion problem:
#     - An initial condition consisting of 1.0 in four grid cells
#     - Zero concentration on the boundaries
#     - Constant diffusion
#
#     A small twist: there is an inner boundary that hinders
#     diffusion, though it is modelled too simply. I should
#     suppress the diffusion on the sides on the grid cells
#     rather than suppress the change in the grid cells.
#     Visually the effect is very similar.
#
package require Tk

lappend auto_path .
package require Tensor

# nextTime --
#     Compute the result for the next time
#
# Arguments:
#     None
#
# Side effects:
#     Tensor u updated (du is an auxiliary tensor)
#
proc nextTime {} {

    set ucoeff [expr {$::diffusion * $::deltt / (2.0*$::deltx**2)}]

    set end   [expr {$::number}]
    set endm1 [expr {$::number-1}]
    set endp1 [expr {$::number+1}]

    #
    # Bounds depend on the type of boundary condition - unfortunately
    #
    ::tensor::expr "du(i=1:$end,j=1:$end) = 4.0*u(i=1:$end,j=1:$end) - u(i=0:$endm1,j=1:$end) - u(i=2:$endp1,j=1:$end)
                                                                     - u(i=1:$end,j=0:$endm1) - u(i=1:$end,j=2:$endp1)"

    du *= tensor mask

    ::tensor::expr "u(i=1:$end,j=1:$end) = u(i=1:$end,j=1:$end) - $ucoeff * du(i=1:$end,j=1:$end)"

}

# setBoundaryConditions --
#     Set the boundary conditions - almost trivial in this case
#
# Arguments:
#     time           Time in the computation
#
# Side effects:
#     Tensor u is updated
#
proc setBoundaryConditions {time} {
    set end $::number

    u section [list 0             [list 0 $end]] = scalar 0.0
    u section [list $end          [list 0 $end]] = scalar 0.0
    u section [list [list 0 $end] 0            ] = scalar 0.0
    u section [list [list 0 $end] $end         ] = scalar 0.0
}

# initialisePlot --
#     Initialise the plot
#
# Arguments:
#     None
#
# Side effects:
#     Graph is updated
#
proc initialisePlot {} {

    set width  [.c cget -width]
    set height [.c cget -height]

    set wrect  [expr {$width  / ($::number+2)}]
    set hrect  [expr {$height / ($::number+2)}]

    for { set j 0 } { $j < $::number+1 } { incr j } {
        for { set i 0 } { $i < $::number+1 } { incr i } {
            set xl [expr { $i * $wrect }]
            set xr [expr { $xl + $wrect - 1}]
            set yu [expr { $j * $wrect }]
            set yd [expr { $yu + $wrect - 1}]

            set id [.c create rectangle $xl $yu $xr $yd -fill white]

            cell section [list $i $j] = scalar $id
        }
    }

    .c create text [expr {$width/2}] $height -text "Time: 0.0" -anchor s -tag time
}

# plotData --
#     Plot the result
#
# Arguments:
#     time           Time in the computation
#     values         List of values
#
# Side effects:
#     Graph is updated
#
proc plotData {time} {

    set colours [list darkblue blue  cyan green yellow orange red]
    set values  [list 0.0001   0.003 0.01 0.03  0.1    0.2    0.5]

    .c itemconfigure time -text "Time: [format "%5.2f" $time]"

    for { set j 0 } { $j < $::number+1 } { incr j } {
        for { set i 0 } { $i < $::number+1 } { incr i } {

            set id     [cell section [list $i $j]]
            set value  [u    section [list $i $j]]
            set active [mask section [list $i $j]]

            if { $active == 0.0 } {
                continue
            }
            foreach v $values c $colours {
                if { $value <= $v } {
                    break
                }
            }

            .c itemconfigure $id -fill $c
        }
    }

    after 100 {
        set go 1
    }
    vwait go
}

# main --
#     Create the tensors, store a few
#     computational paramters and then
#     run the time loop
#
set diffusion  0.1
set number    30
set length    10.0
set deltx      [expr {$length/$number}]
set deltt      0.001

set size       [expr {$number + 2}]

tensor::create u    -size [list $size $size] -type double
tensor::create du   -size [list $size $size] -type double

tensor::create cell -size [list $size $size] -type int
tensor::create mask -size [list $size $size] -type double

#
# Initial condition
#
u  = scalar 0.0
du = scalar 0.0

set mid [expr {$number/2}]
u section [list [list $mid [expr {$mid+1}]] [list $mid [expr {$mid+1}]]] = scalar 1.0

#
# Set several cells to inactive
#
mask = scalar 1.0

set inactiveRow [expr {$mid + 4}]
mask section [list $inactiveRow [list [expr {$mid - 4}] [expr {$mid + 4}]]] = scalar 0.0

#
# Loop over time
#
pack [canvas .c -bg white -width 500 -height 500]

tkwait visibility .c

initialisePlot

set t 0
while { $t < 100000} {

    set time [expr {$deltt * $t}]

    if { $t % 200 == 0 } {
        plotData $time
    }

    setBoundaryConditions $time

    nextTime

    incr t
}