begin
    import Pkg as _Pkg
    haskey(ENV, "PLUTO_PROJECT") && _Pkg.activate(ENV["PLUTO_PROJECT"])
    using Revise
    using Test
    using SimplexGridFactory, Triangulate, ExtendableGrids, VoronoiFVM
    using PlutoUI, GridVisualize
    using CairoMakie
    CairoMakie.activate!(; type = "png", visible = false)
    GridVisualize.default_plotter!(CairoMakie)
end;

Grid

function swiss_cheese_2d()
    function circlehole!(builder, center, radius; n = 20)
        points = [point!(builder, center[1] + radius * sin(t), center[2] + radius * cos(t))
                  for t in range(0, 2π; length = n)]
        for i = 1:(n - 1)
            facet!(builder, points[i], points[i + 1])
        end
        facet!(builder, points[end], points[1])
        holepoint!(builder, center)
    end

    builder = SimplexGridBuilder(; Generator = Triangulate)
    cellregion!(builder, 1)
    maxvolume!(builder, 0.1)
    regionpoint!(builder, 0.1, 0.1)

    p1 = point!(builder, 0, 0)
    p2 = point!(builder, 10, 0)
    p3 = point!(builder, 10, 10)
    p4 = point!(builder, 0, 10)

    facetregion!(builder, 1)
    facet!(builder, p1, p2)
    facet!(builder, p2, p3)
    facet!(builder, p3, p4)
    facet!(builder, p4, p1)

    holes = [1.0 2.0
             8.0 9.0
             2.0 8.0
             8.0 4.0
             9.0 1.0
             3.0 4.0
             4.0 6.0
             7.0 9.0
             4.0 7.0
             7.0 5.0
             2.0 1.0
             4.0 1.0
             4.0 8.0
             3.0 6.0
             4.0 9.0
             6.0 9.0
             3.0 5.0
             1.0 4.0]'

    radii = [
        0.15,
        0.15,
        0.1,
        0.35,
        0.2,
        0.3,
        0.1,
        0.4,
        0.1,
        0.4,
        0.2,
        0.2,
        0.2,
        0.35,
        0.15,
        0.25,
        0.15,
        0.25,
    ]

    for i = 1:length(radii)
        facetregion!(builder, i + 1)
        circlehole!(builder, holes[:, i], radii[i])
    end

    simplexgrid(builder)
end

swiss_cheese_2d (generic function with 1 method)

grid = swiss_cheese_2d()

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       2
   nnodes =    1434
   ncells =    2443
  nbfaces =     459

System + solution

mutable struct Params
    val11::Float64
end

params = Params(5)

Params(5.0)

flux(y, u, edge, data) = y[1] = u[1, 1] - u[1, 2];

function bc(y, u, bnode, data)
    boundary_dirichlet!(y, u, bnode; region = 2, value = 10.0)
    boundary_dirichlet!(y, u, bnode; region = 3, value = 0.0)
    boundary_dirichlet!(y, u, bnode; region = 11, value = data.val11)
end

bc (generic function with 1 method)

system = VoronoiFVM.System(grid; flux = flux, species = 1, bcondition = bc, data = params)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=2, nnodes=1434, ncells=2443,
  nbfaces=459),  
  physics = Physics(data=Main.var"workspace#3".Params, flux=flux,
  storage=default_storage, breaction=bc, ),  
  num_species = 1)

begin
    params.val11 = val11
    sol = solve(system)
end;

@test params.val11 != 5.0 || isapprox(sum(sol), 7842.2173682050525; rtol = 1.0e-12)

Test Passed

Flux reconstruction

nf = nodeflux(system, sol)

2×1×1434 Array{Float64, 3}:
[:, :, 1] =
  0.05468367090633706
 -0.05468367090633706

[:, :, 2] =
 0.01852257911924369
 0.014818063295393813

[:, :, 3] =
 0.0
 0.0

;;; … 

[:, :, 1432] =
 0.020523328431052094
 0.036034112502081016

[:, :, 1433] =
 0.028906869669800477
 0.012529162794698063

[:, :, 1434] =
 0.03545781788403497
 0.0342143068249023

@test params.val11 != 5.0 || isapprox(sum(nf), 978.000534849034; rtol = 1.0e-14)

Test Passed

vis = GridVisualizer(; dim = 2, resolution = (400, 400))

begin
    scalarplot!(vis, grid, sol[1, :]; levels = 9, colormap = :summer, clear = true)
    vectorplot!(vis, grid, nf[:, 1, :]; clear = false, vscale = 1.5)
    reveal(vis)
end

src(x) = exp(-x^2 / 0.01)

src (generic function with 1 method)

source1d(y, node, data) = y[1] = src(node[1])

source1d (generic function with 1 method)

flux1d(y, u, edge, data) = y[1] = u[1, 1]^2 - u[1, 2]^2

flux1d (generic function with 1 method)

function bc1d(y, u, bnode, data)
    boundary_dirichlet!(y, u, bnode; region = 1, value = 0.01)
    boundary_dirichlet!(y, u, bnode; region = 2, value = 0.01)
end

bc1d (generic function with 1 method)

grid1d = simplexgrid(-1:0.01:1)

ExtendableGrids.ExtendableGrid{Float64, Int32}
      dim =       1
   nnodes =     201
   ncells =     200
  nbfaces =       2

sys1d = VoronoiFVM.System(grid1d;
                          flux = flux1d,
                          bcondition = bc1d,
                          source = source1d,
                          species = [1],)

VoronoiFVM.System{Float64, Float64, Int32, Int64, Matrix{Int32}}(
  grid = ExtendableGrids.ExtendableGrid{Float64, Int32}(dim=1, nnodes=201, ncells=200,
  nbfaces=2),  
  physics = Physics(flux=flux1d, storage=default_storage, source=source1d,
  breaction=bc1d, ),  
  num_species = 1)

sol1d = solve(sys1d; inival = 0.1)

1×201 VoronoiFVM.DenseSolutionArray{Float64, 2}:
 0.01  0.0314043  0.0432719  0.0525231  …  0.0525231  0.0432719  0.0314043  0.01

nf1d = nodeflux(sys1d, sol1d)

1×1×201 Array{Float64, 3}:
[:, :, 1] =
 -0.08862269254527583

[:, :, 2] =
 -0.08862269254527581

[:, :, 3] =
 -0.08862269254527581

;;; … 

[:, :, 199] =
 0.08862269254527581

[:, :, 200] =
 0.08862269254527581

[:, :, 201] =
 0.08862269254527583

let
    vis1d = GridVisualizer(; dim = 1, resolution = (500, 250), legend = :lt)
    scalarplot!(vis1d, grid1d, map(src, grid1d); label = "rhs", color = :blue)
    scalarplot!(vis1d, grid1d, sol1d[1, :]; label = "solution", color = :red, clear = false)
    vectorplot!(vis1d, grid1d, nf1d[:, 1, :]; label = "flux", clear = false, color = :green)
    reveal(vis1d)
end

@test nf1d[1, 1, 101] ≈ 0.0

Test Passed

@test nf1d[1, 1, 1] ≈ -nf1d[1, 1, end]

Test Passed