Fix anastomosis (#70)

* add solved flag

* 2.6.1

* ifelse -> max/min

* cleanup

Project.toml

@@ -1,7 +1,7 @@
 name = "openBF"
 uuid = "e815b1a4-10eb-11ea-25f1-272ff651e618"
 authors = ["alessandro <alessandro_melis@rocketmail.com>"]
-version = "2.6.0"
+version = "2.6.1"
 
 [deps]
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"

src/anastomosis.jl

@@ -26,7 +26,7 @@ getUan(v1::Vessel, v2::Vessel, v3::Vessel) = SVector{6, Float64}(
 function solveAnastomosis(v1::Vessel, v2::Vessel, v3::Vessel)
     k = (sqrt(1.5*v1.gamma[end]), sqrt(1.5*v2.gamma[end]), sqrt(1.5*v3.gamma[1]))
     U = getUan(v1, v2, v3)
-    W = (U[1] + 4k[1]* U[4], U[2] + 4k[2] * U[5], U[3] + 4k[3] * U[6])
+    W = (U[1] + 4k[1]* U[4], U[2] + 4k[2] * U[5], U[3] - 4k[3] * U[6])
     J = getJan(v1, v2, v3, U, k)
     F = getFan(v1, v2, v3, U, k, W)
 
@@ -56,7 +56,7 @@ function getFan(v1::Vessel, v2::Vessel, v3::Vessel, U, k, W)
     v2.beta[end] * (U[5]^2 / sqrt(v2.A0[end]) - 1) - (v3.beta[1] * ((U[6]^2) / sqrt(v3.A0[1]) - 1)))
 end
 
-function getJan(v1::Vessel, v2::Vessel, v3::Vessel, U::SArray, k::SArray)
+function getJan(v1::Vessel, v2::Vessel, v3::Vessel, U, k)
     J::Array{Float64, 2} = zeros(Float64, 6, 6)
 
     J[1, 1] = 1.0
@@ -89,5 +89,5 @@ function NRan(U, W, J, F, k, v1, v2, v3)
         F = getFan(v1, v2, v3, U, k, W)
         J = getJan(v1, v2, v3, U, k)
     end
-    J
+    U
 end

src/bifurcations.jl

@@ -24,7 +24,6 @@ function getJbif(v1::Vessel, v2::Vessel, v3::Vessel, U, k)
     J[1, 1] = 1.0
     J[2, 2] = 1.0
     J[3, 3] = 1.0
-    J[4, 4] = 1.0
 
     J[1, 4] = 4k[1]
     J[2, 5] = -4k[2]

src/solver.jl

@@ -22,52 +22,82 @@ function calculateΔt(n::Network)
         maxspeed = 0.0
         for i=eachindex(v.u)
             @inbounds speed = abs(v.u[i] + wave_speed(v.A[i], v.gamma[i+1]))
-            maxspeed = ifelse(speed>maxspeed, speed, maxspeed)
+            maxspeed = max(maxspeed, speed)
         end
         Δt = v.dx / maxspeed
-        minΔt = ifelse(Δt < minΔt, Δt, minΔt)
+        minΔt = min(minΔt, Δt)
     end
     minΔt*n.Ccfl
 end
 
+function skipthis(n::Network, s, t)
+    if n.vessels[(s, t)].solved
+        return true
+    end
+
+    if s == 1
+        return false
+    end
+
+    indeg::Int64 = Graphs.indegree(n.graph, s)
+
+    if indeg == 1
+        parent_src = first(Graphs.inneighbors(n.graph, s))
+        return ~n.vessels[parent_src, s].solved
+    end
+
+    if indeg == 2
+        parent_src = Graphs.inneighbors(n.graph, s)
+        return ~(n.vessels[first(parent_src), s].solved && n.vessels[last(parent_src), s].solved)
+    end
+
+    # TODO: this should error
+
+    return false
+end
+
 function solve!(n::Network, dt::Float64, current_time::Float64)::Nothing
-    for edge in n.edges
-        s = Graphs.src(edge)
-        t = Graphs.dst(edge)
-
-        # inlet
-        s == 1 && inbc!(n.vessels[(s, t)], current_time, dt, n.heart)
-
-        # TODO: multiple inlets
-
-        # vessel
-        muscl!(n.vessels[(s, t)], dt, n.blood)
-
-        # downstream
-        outdeg::Int64 = Graphs.outdegree(n.graph, t)
-        if outdeg == 0 # outlet
-            outbc!(n.vessels[(s, t)], dt, n.blood.rho)
-        elseif outdeg == 1
-            indeg::Int64 = Graphs.indegree(n.graph, t)
-            d::Int64 = first(Graphs.outneighbors(n.graph, t))
-            if indeg == 1 # conjunction
-                join_vessels!(n.vessels[(s, t)], n.vessels[(t, d)], n.blood.rho)
+    for (k, v) in n.vessels
+        v.solved = false
+    end
+
+    while ~all(v.solved for v in values(n.vessels))
+        for edge in n.edges
+            s = Graphs.src(edge)
+            t = Graphs.dst(edge)
 
-            # TODO: test ANASTOMOSIS!!!
-            elseif indeg == 2 # anastomosis
-                ps::Vector{Int64} = Graphs.inneighbors(n.graph, t)
-                if t == max(ps[1], ps[2])
+            skipthis(n, s, t) && continue
+
+            # inlet
+            s == 1 && inbc!(n.vessels[(s, t)], current_time, dt, n.heart)
+
+            # TODO: multiple inlets
+
+            # vessel
+            muscl!(n.vessels[(s, t)], dt, n.blood)
+
+            # downstream
+            outdeg::Int64 = Graphs.outdegree(n.graph, t)
+            if outdeg == 0 # outlet
+                outbc!(n.vessels[(s, t)], dt, n.blood.rho)
+            elseif outdeg == 1
+                indeg::Int64 = Graphs.indegree(n.graph, t)
+                d::Int64 = first(Graphs.outneighbors(n.graph, t))
+                if indeg == 1 # conjunction
+                    join_vessels!(n.vessels[(s, t)], n.vessels[(t, d)], n.blood.rho)
+                elseif indeg == 2 # anastomosis
+                    ps::Vector{Int64} = Graphs.inneighbors(n.graph, t)
                     solveAnastomosis(
                         n.vessels[(ps[1], t)],
                         n.vessels[(ps[2], t)],
                         n.vessels[(t, d)],
                     )
                 end
+            elseif outdeg == 2 # bifurcation
+                ds::Vector{Int64} = Graphs.outneighbors(n.graph, t)
+                join_vessels!(n.vessels[(s, t)], n.vessels[t, ds[1]], n.vessels[t, ds[2]])
             end
-
-        elseif outdeg == 2 # bifurcation
-            ds::Vector{Int64} = Graphs.outneighbors(n.graph, t)
-            join_vessels!(n.vessels[(s, t)], n.vessels[t, ds[1]], n.vessels[t, ds[2]])
+            n.vessels[(s, t)].solved = true
         end
     end
 end

src/vessel.jl

@@ -28,6 +28,7 @@ function Blood(config::Dict)
 end
 
 mutable struct Vessel
+    solved::Bool
     label::SubString{String}
     tosave::Bool
 
@@ -264,7 +265,10 @@ function Vessel(config::Dict{Any,Any}, b::Blood, jump::Int64, tokeep::Vector{Str
         waveforms[q] = zeros(Float64, jump, 6)
     end
 
+    solved = false
+
     Vessel(
+        solved,
         vessel_name,
         tosave,
         sn,