feat: added option to modify singular boundary elements

julia/TUG/src/Boundary.jl

@@ -289,3 +289,37 @@ function setBoundarySideConstant!(
 
     boundary.boundaries[Int(side)] = [BoundaryElement{T}(value) for _ = 1:n]
 end
+
+"""
+    setBoundarySideElement(boundary::Boundary{T}, side::SIDE, index::Int, be::BoundaryElement{T})::BoundaryElement{T} where {T}
+
+Sets the boundary element at the specified side and index.
+
+# Arguments
+- `boundary::Boundary{T}`: The boundary.
+- `side::SIDE`: The side of the boundary.
+- `index::Int`: The index of the boundary element.
+- `be::BoundaryElement{T}`: The boundary element.
+
+# Returns
+The boundary element at the specified side and index.
+"""
+function setBoundarySideElement(
+    boundary::Boundary{T},
+    side::SIDE,
+    index::Int,
+    be::BoundaryElement{T},
+)::BoundaryElement{T} where {T}
+    if boundary.dim == 1 && (side == BOTTOM || side == TOP)
+        throw(
+            ArgumentError(
+                "For the one-dimensional case, only the left and right borders exist.",
+            ),
+        )
+    end
+    if index < 1 || index > (side in (LEFT, RIGHT) ? boundary.rows : boundary.cols)
+        throw(ArgumentError("Index out of bounds!"))
+    end
+
+    boundary.boundaries[Int(side)][index] = be
+end

julia/TUG/src/Core/BTCS.jl

@@ -37,7 +37,6 @@ function createCoeffMatrix(
     rowIndex::Int,
     sx::T,
 )::Tridiagonal{T} where {T}
-    # Determine left side boundary coefficients based on boundary condition
     centerCoeffTop, rightCoeffTop = calcBoundaryCoeff(
         alpha[rowIndex, 1],
         alpha[rowIndex, 2],
@@ -45,7 +44,6 @@ function createCoeffMatrix(
         getType(bcLeft[rowIndex]),
     )
 
-    # Determine right side boundary coefficients based on boundary condition
     centerCoeffBottom, leftCoeffBottom = calcBoundaryCoeff(
         alpha[rowIndex, numCols],
         alpha[rowIndex, numCols-1],
@@ -95,7 +93,6 @@ function calcExplicitConcentrationsBoundaryClosed(
     return (sy .* alpha .+ (1 .- sy .* alpha)) .* conc_center
 end
 
-# creates a solution vector for next time step from the current state of concentrations
 function writeSolutionVector!(
     sv::Vector{T},
     concentrations::Matrix{T},
@@ -112,7 +109,6 @@ function writeSolutionVector!(
     numRows = size(concentrations, 1)
     length = size(sv, 1)
 
-    # First row
     if rowIndex == 1
         @inbounds for i = 1:length
             if getType(bcTop[i]) == CONSTANT
@@ -136,14 +132,13 @@ function writeSolutionVector!(
         end
     end
 
-    # Inner rows
     if rowIndex > 1 && rowIndex < numRows
         @inbounds for i = 1:length
             alpha_here_below =
                 calcAlphaIntercell(alphaY[rowIndex, i], alphaY[rowIndex+1, i])
             alpha_here_above =
                 alphaY[rowIndex+1, i] == alphaY[rowIndex-1, i] ? alpha_here_below :
-                calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i]) # calcAlphaIntercell is symmetric, so we can use it for both directions
+                calcAlphaIntercell(alphaY[rowIndex-1, i], alphaY[rowIndex, i])
             sv[i] =
                 sy * alpha_here_below * concentrations[rowIndex+1, i] +
                 (1 - sy * (alpha_here_below + alpha_here_above)) *
@@ -152,7 +147,6 @@ function writeSolutionVector!(
         end
     end
 
-    # Last row
     if rowIndex == numRows
         @inbounds for i = 1:length
             if getType(bcBottom[i]) == CONSTANT
@@ -176,12 +170,10 @@ function writeSolutionVector!(
         end
     end
 
-    # First column - additional fixed concentration change from perpendicular dimension in constant BC case
     if getType(bcLeft[rowIndex]) == CONSTANT
         sv[1] += 2 * sx * alphaX[rowIndex, 1] * getValue(bcLeft[rowIndex])
     end
 
-    # Last column - additional fixed concentration change from perpendicular dimension in constant BC case
     if getType(bcRight[rowIndex]) == CONSTANT
         sv[end] += 2 * sx * alphaX[rowIndex, end] * getValue(bcRight[rowIndex])
     end

julia/TUG/src/Core/FTCS.jl

@@ -14,7 +14,6 @@ function calcHorizontalChange(
     concentration_prev::T,
     concentration_current::T,
 ) where {T}
-
     intercellAlpha_next = calcAlphaIntercell(alphaX_next, alphaX_current)
     intercellAlpha_prev = calcAlphaIntercell(alphaX_prev, alphaX_current)
 
@@ -134,7 +133,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
 
     row = 1
 
-    # inner cells
     for col = 2:colMax-1
         concentrations_t1[row, col] +=
             sx * calcHorizontalChange(
@@ -147,7 +145,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             )
     end
 
-    # left boundary
     leftBoundaryType = getBoundaryElementType(bc, LEFT, row)
     leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row)
     concentrations_t1[row, 1] +=
@@ -160,7 +157,6 @@ function FTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             leftBoundaryValue,
         )
 
-    # right boundary
     rightBoundaryType = getBoundaryElementType(bc, RIGHT, row)
     rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row)
     concentrations_t1[row, colMax] +=
@@ -187,12 +183,6 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
     concentrations = getConcentrations(grid)
     concentrations_t1 = copy(concentrations)
 
-    # Currently boundary sides can only be set for the whole side, not for each cell individually
-    leftBoundaryType = getBoundaryElementType(bc, LEFT, 1)
-    rightBoundaryType = getBoundaryElementType(bc, RIGHT, 1)
-    topBoundaryType = getBoundaryElementType(bc, TOP, 1)
-    bottomBoundaryType = getBoundaryElementType(bc, BOTTOM, 1)
-
     Threads.@threads for row = 2:rowMax-1
         for col = 2:colMax-1
             concentrations_t1[row, col] +=
@@ -214,17 +204,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
                 )
         end
 
-        # Boundary conditions for each row
-        # Left boundary without corners
-        leftBoundaryValue = getBoundaryElementValue(bc, LEFT, row)
         concentrations_t1[row, 1] +=
             sx * calcHorizontalChangeLeftBoundary(
-                leftBoundaryType,
+                getBoundaryElementType(bc, LEFT, row),
                 alphaX[row, 2],
                 alphaX[row, 1],
                 concentrations[row, 2],
                 concentrations[row, 1],
-                leftBoundaryValue,
+                getBoundaryElementValue(bc, LEFT, row),
             ) +
             sy * calcVerticalChange(
                 alphaY[row+1, 1],
@@ -235,16 +222,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
                 concentrations[row, 1],
             )
 
-        # Right boundary without corners
-        rightBoundaryValue = getBoundaryElementValue(bc, RIGHT, row)
         concentrations_t1[row, colMax] +=
             sx * calcHorizontalChangeRightBoundary(
-                rightBoundaryType,
+                getBoundaryElementType(bc, RIGHT, row),
                 alphaX[row, colMax-1],
                 alphaX[row, colMax],
                 concentrations[row, colMax-1],
                 concentrations[row, colMax],
-                rightBoundaryValue,
+                getBoundaryElementValue(bc, RIGHT, row),
             ) +
             sy * calcVerticalChange(
                 alphaY[row+1, colMax],
@@ -256,18 +241,15 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             )
     end
 
-    # Handle top/bottom boundaries
     Threads.@threads for col = 2:colMax-1
-        # Top boundary
-        topBoundaryValue = getBoundaryElementValue(bc, TOP, col)
         concentrations_t1[1, col] +=
             sy * calcVerticalChangeTopBoundary(
-                topBoundaryType,
+                getBoundaryElementType(bc, TOP, col),
                 alphaY[2, col],
                 alphaY[1, col],
                 concentrations[2, col],
                 concentrations[1, col],
-                topBoundaryValue,
+                getBoundaryElementValue(bc, TOP, col),
             ) +
             sx * calcHorizontalChange(
                 alphaX[1, col+1],
@@ -278,16 +260,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
                 concentrations[1, col],
             )
 
-        # Bottom boundary
-        bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, col)
         concentrations_t1[rowMax, col] +=
             sy * calcVerticalChangeBottomBoundary(
-                bottomBoundaryType,
+                getBoundaryElementType(bc, BOTTOM, col),
                 alphaY[rowMax, col],
                 alphaY[rowMax-1, col],
                 concentrations[rowMax, col],
                 concentrations[rowMax-1, col],
-                bottomBoundaryValue,
+                getBoundaryElementValue(bc, BOTTOM, col),
             ) +
             sx * calcHorizontalChange(
                 alphaX[rowMax, col+1],
@@ -299,9 +279,6 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             )
     end
 
-    # Handle corners
-    # Top left corner
-    topBoundaryValue = getBoundaryElementValue(bc, TOP, 1)
     concentrations_t1[1, 1] +=
         sx * calcHorizontalChange(
             alphaX[1, 2],
@@ -312,16 +289,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             concentrations[1, 1],
         ) +
         sy * calcVerticalChangeTopBoundary(
-            topBoundaryType,
+            getBoundaryElementType(bc, TOP, 1),
             alphaY[2, 1],
             alphaY[1, 1],
             concentrations[2, 1],
             concentrations[1, 1],
-            topBoundaryValue,
+            getBoundaryElementValue(bc, TOP, 1),
         )
 
-    # Top right corner
-    topBoundaryValue = getBoundaryElementValue(bc, TOP, colMax)
     concentrations_t1[1, colMax] +=
         sx * calcHorizontalChange(
             alphaX[1, colMax-1],
@@ -332,16 +307,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             concentrations[1, colMax],
         ) +
         sy * calcVerticalChangeTopBoundary(
-            topBoundaryType,
+            getBoundaryElementType(bc, TOP, colMax),
             alphaY[2, colMax],
             alphaY[1, colMax],
             concentrations[2, colMax],
             concentrations[1, colMax],
-            topBoundaryValue,
+            getBoundaryElementValue(bc, TOP, colMax),
         )
 
-    # Bottom left corner
-    bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, 1)
     concentrations_t1[rowMax, 1] +=
         sx * calcHorizontalChange(
             alphaX[rowMax, 2],
@@ -352,16 +325,14 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             concentrations[rowMax, 1],
         ) +
         sy * calcVerticalChangeBottomBoundary(
-            bottomBoundaryType,
+            getBoundaryElementType(bc, BOTTOM, 1),
             alphaY[rowMax, 1],
             alphaY[rowMax-1, 1],
             concentrations[rowMax, 1],
             concentrations[rowMax-1, 1],
-            bottomBoundaryValue,
+            getBoundaryElementValue(bc, BOTTOM, 1),
         )
 
-    # Bottom right corner
-    bottomBoundaryValue = getBoundaryElementValue(bc, BOTTOM, colMax)
     concentrations_t1[rowMax, colMax] +=
         sx * calcHorizontalChange(
             alphaX[rowMax, colMax-1],
@@ -372,12 +343,12 @@ function FTCS_2D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
             concentrations[rowMax, colMax],
         ) +
         sy * calcVerticalChangeBottomBoundary(
-            bottomBoundaryType,
+            getBoundaryElementType(bc, BOTTOM, colMax),
             alphaY[rowMax, colMax],
             alphaY[rowMax-1, colMax],
             concentrations[rowMax, colMax],
             concentrations[rowMax-1, colMax],
-            bottomBoundaryValue,
+            getBoundaryElementValue(bc, BOTTOM, colMax),
         )
 
     setConcentrations!(grid, concentrations_t1)

julia/TUG/src/TUG.jl

@@ -24,11 +24,12 @@ export Boundary,
     BoundaryElement, TYPE, SIDE, BC_TYPE_CLOSED, BC_TYPE_CONSTANT, LEFT, RIGHT, TOP, BOTTOM
 export getType,
     getValue,
-    setBoundarySideClosed!,
-    setBoundarySideConstant!,
     getBoundarySide,
     getBoundaryElementType,
-    getBoundaryElementValue
+    getBoundaryElementValue,
+    setBoundarySideClosed!,
+    setBoundarySideConstant!,
+    setBoundarySideElement
 
 include("AbstractSimulation.jl")