start on OLIM

src/OLIM/OLIM.jl

@@ -0,0 +1,56 @@
+module OLIM
+
+const INFTY = Inf
+const TOL = 1.0e-12
+const NX = 1024
+const NY = 1024
+const K = 22
+const NCURVEMAX = 399
+
+const nx1 = NX - 1
+const ny1 = NY - 1
+const nxy = NX * NY
+
+const KK = K * K
+
+NCURVE = 0
+count = 0
+h = 0.0
+hx = 0.0
+hy = 0.0
+
+ms = fill(0, NX * NY)
+g = fill(INFTY, NX * NY)
+pos = fill(0, NX * NY)
+tree = fill(0, NX * NY)
+
+const neii = [1, NX + 1, NX, NX - 1, -1, -NX - 1, -NX, -NX + 1]
+
+
+Uexact = fill(0.0, NX * NY)
+solinfo = [SolInfo('n', 0, 0, 0.0) for _ in 1:NX*NY]
+NXY = NX * NY
+
+
+NMESH = 0
+x_ipoint = [0.0, 0.0]
+icurve = Vector{MyCurve}(undef, NCURVEMAX)
+
+
+# Variables for the potential
+# // define the computational domain
+XMIN = 0.0
+XMAX = 0.0
+YMIN = 0.0
+YMAX = 0.0
+
+
+
+
+
+include("types.jl")
+include("point.jl")
+include("field.jl")
+
+
+end # module OLIM

src/OLIM/field.jl

@@ -0,0 +1,101 @@
+
+# // choose the vector field b
+# const char chfield='l';
+# // Vector fields with asymptotically stable equilibrium
+# // chfield = 'l': b = [-2, -alin; 2*alin, -1][ x; y]; (in the Matlab notations),
+# // analytic solution  U = 2x^2 + y^2
+# const double alin = 10.0;
+# // chfield = 'q' ->  Maier-Stein model;
+# // chfield = 'r' -> FitzHugh-Nagumo model with a = 1.2
+
+# // Vector fields with stable limit cycle
+# // chfield = 'b' ->  Brusselator;
+# // CHANGE const char *ficurve_name = "circle.txt"; to
+# // const char *ficurve_name = "icurve_brusselator.txt";
+# // chfield = 'c' -> analytic solution U = (1/2)(r^2-1)^2, l = (y,-x)
+
+const alin = 10.0;
+function myfield(x::Vector{Float64}, chfield::Char)::Vector{Float64}
+    v = zeros(2)
+
+    if chfield == 'b'
+        v[1] = 1.0 + x[1]^2 * x[2] - 4.0 * x[1]
+        v[2] = 3 * x[1] - x[1]^2 * x[2]
+    elseif chfield == 'c'
+        aux = 1.0 - x[1]^2 - x[2]^2
+        v[1] = x[2] + x[1] * aux
+        v[2] = -x[1] + x[2] * aux
+    elseif chfield == 'q'
+        v[1] = x[1] - x[1]^3 - 10.0 * x[1] * x[2]^2
+        v[2] = -(1.0 + x[1]^2) * x[2]
+    elseif chfield == 'r'
+        v[1] = x[1] - x[1]^3/3.0 - x[2]
+        v[2] = x[1] + 1.2
+    elseif chfield == 'l'
+        v[1] = -2 * x[1] - alin * x[2]
+        v[2] = 2 * alin * x[1] - x[2]
+    else
+        error("Invalid chfield: $chfield")
+    end
+
+    return v
+end
+
+function exact_solution(x::Vector{Float64}, chfield::Char)::Float64
+    if chfield == 'l'
+        return 2.0 * x[1]^2 + x[2]^2
+    elseif chfield == 'c'
+        return 0.5 * (x[1]^2 + x[2]^2 - 1.0)^2
+    else
+        return 0.0
+    end
+end
+
+function param(chfield::Char)
+    # Define parameters based on the field type
+    x_ipoint = Dict(
+        'q' => [-1.0, 0.0],
+        'r' => [-1.2, -1.2*(1.0-1.2*1.2/3.0)],
+        'l' => [0.0, 0.0]
+    )
+
+    bounds = Dict(
+        'q' => (-2.0, 2.0, -2.0, 2.0),
+        'r' => (-2.5, 2.5, -2.5, 2.5),
+        'l' => (-1.0, 1.0, -1.0, 1.0),
+        'c' => (-2.0, 2.0, -2.0, 2.0),
+        'b' => (0.0, 7.0, 0.0, 7.0)
+    )
+
+    if !haskey(bounds, chfield)
+        error("Invalid chfield: $chfield")
+    end
+
+    XMIN, XMAX, YMIN, YMAX = bounds[chfield]
+
+    println("in param()")
+
+    # Grid parameters
+    hx = (XMAX - XMIN)/(NX-1)
+    hy = (YMAX - YMIN)/(NY-1)
+    h = sqrt(hx^2 + hy^2)
+
+    # Initialize arrays
+    ms = zeros(Int, NX * NY)
+    g = fill(Inf, NX * NY)
+    solinfo = fill('n', NX * NY)
+
+    # Compute exact solution if available
+    Uexact = zeros(NX * NY)
+    if chfield in ['l', 'c']
+        for j in 0:NY-1
+            for i in 0:NX-1
+                ind = i + 1 + NX * j
+                x = [XMIN + i*hx, YMIN + j*hy]
+                Uexact[ind] = exact_solution(x, chfield)
+            end
+        end
+    end
+
+    return hx, hy, h, ms, g, solinfo, Uexact, x_ipoint
+end

src/OLIM/init.jl

@@ -0,0 +1,160 @@
+function initial_curve()
+    println("in initial_curve()")
+
+    # Read initial curve points from file
+    if !isfile(ficurve_name)
+        error("Please provide a data file for the initial curve")
+    end
+
+    # Initialize curve from file
+    points = Vector{Tuple{Float64,Float64}}()
+    open(ficurve_name, "r") do f
+        for line in eachline(f)
+            x, y = parse.(Float64, split(line))
+            push!(points, (x, y))
+        end
+    end
+
+    n = length(points)
+    if n > NCURVEMAX
+        error("Error: j = $n while NCURVEMAX = $NCURVEMAX")
+    end
+
+    # Create linked list of curve points
+    global NCURVE = n
+    global icurve = Vector{MyCurve}(undef, n)
+
+    for i in 1:n
+        x, y = points[i]
+        next_idx = i == n ? 1 : i + 1
+        prev_idx = i == 1 ? n : i - 1
+        icurve[i] = MyCurve(x, y, nothing, nothing)
+    end
+
+    # Link the nodes
+    for i in 1:n
+        icurve[i] = MyCurve(
+            icurve[i].x,
+            icurve[i].y,
+            i == n ? icurve[1] : icurve[i+1],
+            i == 1 ? icurve[n] : icurve[i-1]
+        )
+    end
+
+    println("start")
+    Nac = 0
+    iac = zeros(Int, NX*100)
+
+    # Process each curve segment
+    for n in 1:NCURVE
+        xc0, yc0 = icurve[n].x, icurve[n].y
+        next_point = icurve[n].next
+        xc1, yc1 = next_point.x, next_point.y
+
+        i0 = min(floor(Int, (xc0 - XMIN)/hx), floor(Int, (xc1 - XMIN)/hx))
+        j0 = min(floor(Int, (yc0 - YMIN)/hy), floor(Int, (yc1 - YMIN)/hy))
+        i1 = max(ceil(Int, (xc0 - XMIN)/hx), ceil(Int, (xc1 - XMIN)/hx))
+        j1 = max(ceil(Int, (yc0 - YMIN)/hy), ceil(Int, (yc1 - YMIN)/hy))
+
+        for i in i0:i1, j in j0:j1
+            ind = i + j*NX
+            if ms[ind] == 0
+                g[ind] = init(getpoint(ind))
+                ms[ind] = 1
+                solinfo[ind] = SolInfo('0', 0, 0, 0.0)
+                addtree(ind)
+                Nac += 1
+                iac[Nac] = ind
+            end
+        end
+    end
+end
+
+function init(x::MyVector)::Float64
+    if chfield == 'l'
+        x = vec_difference(x, x_ipoint)
+        a, b = -2.0, -alin
+        c, d = 2.0*alin, -1.0
+        aux1 = c - b
+        aux2 = a + d
+        aux = aux1*aux1 + aux2*aux2
+        aux1 *= aux2/aux
+        aux2 *= aux2/aux
+        A = -(a*aux2 + c*aux1)
+        B = -(b*aux2 + d*aux1)
+        C = -(d*aux2 - b*aux1)
+        return A*x.x*x.x + 2.0*B*x.x*x.y + C*x.y*x.y
+
+    elseif chfield == 'c'
+        # Find closest point on curve
+        step = max(1, NCURVE÷12)
+        nmin, dmin = 0, INFTY
+
+        # Course search
+        for n in 1:step:NCURVE
+            dtemp = length(icurve[n].x - x.x, icurve[n].y - x.y)
+            if dtemp < dmin
+                dmin = dtemp
+                nmin = n
+            end
+        end
+
+        # Fine search forward
+        imin = nmin
+        step *= 2
+        curr = icurve[nmin]
+        for _ in 1:step
+            curr = curr.next
+            dtemp = length(curr.x - x.x, curr.y - x.y)
+            if dtemp < dmin
+                dmin = dtemp
+                imin = (nmin + n) % NCURVE
+            end
+        end
+
+        # Fine search backward
+        curr = icurve[nmin]
+        for _ in 1:step
+            curr = curr.prev
+            dtemp = length(curr.x - x.x, curr.y - x.y)
+            if dtemp < dmin
+                dmin = dtemp
+                imin = (nmin + NCURVE - n) % NCURVE
+            end
+        end
+
+        # Calculate result
+        x0 = MyVector(icurve[imin].x, icurve[imin].y)
+        x1 = MyVector(icurve[imin].next.x, icurve[imin].next.y)
+        v0 = vec_difference(x0, x)
+        v1 = vec_difference(x1, x)
+
+        if dot_product(v0, vec_difference(v0, v1)) < 0.0
+            x1 = MyVector(icurve[imin].prev.x, icurve[imin].prev.y)
+            v1 = vec_difference(x1, x)
+        end
+
+        v = vec_difference(x1, x0)
+        c = length_vec(v)
+        dmin = abs(v0.x*v1.y - v0.y*v1.x)/c
+        a = sqrt(v0.x*v0.x + v0.y*v0.y - dmin*dmin)
+        aux = a/c
+        v = MyVector(v.x * aux, v.y * aux)
+        x0 = vec_sum(x0, v)
+        x1 = vec_lin_comb(x, x0, 0.5, 0.5)
+        bvec0 = myfield(x0)
+        bvec1 = myfield(x1)
+        bvec = myfield(x)
+        lb0 = length_vec(bvec0)
+        aux = dot_product(bvec0, bvec1)/(lb0*lb0)
+        bvec0.x *= aux
+        bvec0.y *= aux
+        aux = dot_product(bvec0, bvec)/(lb0*lb0)
+        bvec.x *= aux
+        bvec.y *= aux
+        temp = dmin*(4.0*length_vec(vec_difference(bvec1, bvec0)) + length_vec(vec_difference(bvec, bvec0)))/3.0
+        return temp
+    else
+        return 0.0
+    end
+end

src/OLIM/point.jl

@@ -0,0 +1,31 @@
+
+
+function getpoint(ind)
+    if ind < NXY
+        x = hx * (mod(ind, NX)) + XMIN
+        y = hy * (div(ind, NX)) + YMIN
+        return (x=x, y=y)  # Returns a NamedTuple as a vector-like structure
+    end
+end
+
+
+function ipoint()
+    # Array for surface indices offset
+    isur = [0, 1, NX+1, NX]
+
+    # Calculate initial indices
+    i = floor(Int, (x_ipoint.x - XMIN)/hx)
+    j = floor(Int, (x_ipoint.y - YMIN)/hy)
+    ind0 = i + j*NX
+
+    # Process the four surrounding points
+    for m in 1:4
+        ind = ind0 + isur[m]
+        x = getpoint(ind)
+        gtemp = init(x)
+        g[ind] = gtemp
+        ms[ind] = 1
+        addtree(ind)
+        solinfo[ind].type = 0
+    end
+end

src/OLIM/types.jl

@@ -0,0 +1,19 @@
+
+struct MyCurve
+    x::Float64
+    y::Float64
+    next::Union{MyCurve, Nothing}
+    prev::Union{MyCurve, Nothing}
+end
+
+struct MySol
+    g::Float64
+    c::Char
+end
+
+struct SolInfo
+    type::Char # solution type: 1 = 1ptupdate, 2 = 2ptupdate, 0 = initialization, 'n' = never reached
+    ind0::Int
+    ind1::Int
+    s::Float64
+end