Some code to test and play around

The following code works as it is. It takes as parameters a function $u$, a field $\kappa$ and the dimension dim. It then calculates $-\nabla(\kappa\nabla u)=:f$ and creates discrete versions of $\kappa$ and $f$ to generate a numerical solution $u$. Finally it compares the numerical and the exact solution in $L^2$ and plots the numerical result in case dim==2.

The boundary conditions below are set periodic in dimension 1 , Dirichlet at $x_2=1$ and Neumann at $x_2=0$ but you can change this according to your whishes.

The nodes distribution is as standard iid but you can provide a density. The code is implemented to generate 1000 nodes but you can change this as well.

using HighVoronoi
using SparseArrays
using IterativeSolvers
using NearestNeighbors
using LinearAlgebra
using StaticArrays


##########################################################################################

## Derivatives

##########################################################################################

function ∂_k(f,x,k;dim=length(x),vec2=MVector{dim}(zeros(Float64,dim)))
    h = 0.00245
    vec2 .= x
    vec2[k] += h
    f1 = f(vec2)
    vec2[k] += h
    f2 = f(vec2)
    vec2 .= x
    vec2[k] -= h
    f3 = f(vec2)
    vec2[k] -= h
    f4 = f(vec2)
    return ( 8*(f1-f3) + (f4-f2) ) / ( 12*h ) # five-point stencil
end

function ∇(f::Function,dim)
    vec2=MVector{dim}(zeros(Float64,dim))
    return x->map(k->∂_k(f,x,k,dim=dim,vec2=vec2),1:dim)
end

function ∇_buffered(f::Function,dim,vec=MVector{dim}(zeros(Float64,dim)),base=HighVoronoi.empty_local_Base(dim))
    vec = MVector{dim}(zeros(Float64,dim))
    vec2 = MVector{dim}(zeros(Float64,dim))
    return x->map!(k->∂_k(f,x,k,dim=dim,vec2=vec2),vec,1:dim)
end

function ∇cdot(f::Function,dim)
    function sum_partials(f,x,dim,vec2)
        f_sum = 0.0
        for k in 1:dim
            f_sum += ∂_k(y->f(y)[k],x,k,dim=dim,vec2=vec2)
        end
        return f_sum
    end
    vec = MVector{dim}(zeros(Float64,dim))
    return x->sum_partials(f,x,dim,vec)
end

function neumann_bc(flux,domain,x)
    function plane_of_x(domain,x)
        k = 0
        dist = 10.0
        ldp = length(domain.planes)
        for i in 1:ldp
            d = dot(domain.planes[i].base-x,domain.planes[i].normal)
            if d<dist
                k=i
            end
        end
        return domain.planes[k].normal
    end
    normal = plane_of_x(domain,x)
    return dot(normal,flux(x))
end


##########################################################################################

## Solving -∇⋅(κ∇u) = RHS    on  'domain'
## keep track of signs!!!

##########################################################################################

# returns all necessary data to perform a numerical calculation to solve 
# -∇⋅(κ∇u) = RHS   using HighVoronoi tools
# calculates RHS := -∇⋅(κ∇u) and if needed Neumann condition
# provides u_exact:=u and κ
function make_set(u::Function,κ::Function,dim;periodic=[], 
                    dirichlet_boundary=collect(1:(2*dim)), 
                    neumann_boundary=nothing, density=nothing, 
                    number_of_nodes=1000)
    my_domain = cuboid(dim,periodic=periodic)
    ∇u = ∇(u,dim)
    rhs = ∇cdot(x->-1.0*κ(x)*∇u(x),dim)
    flux(x) = -κ(x)*∇u(x)
    neumann_flux(x) = -κ(x)*∇_buff_u(x) # faster but with internal buffer
    if density!=nothing
        return (u_exact = u,  κ=κ, RHS = rhs, 
                domain=my_domain, dim=dim,
                dirichlet_boundary=dirichlet_boundary, 
                neumann_boundary=neumann_boundary, 
                number_of_nodes = number_of_nodes,
                neumann = x->neumann_bc(flux,my_domain,x))
    else
        return (u_exact = u,  κ=κ, RHS = rhs, 
                domain=my_domain, dim=dim,
                density=density,
                dirichlet_boundary=dirichlet_boundary, 
                neumann_boundary=neumann_boundary, 
                number_of_nodes = number_of_nodes,
                neumann = x->neumann_bc(flux,my_domain,x))
    end
end


using Plots
# plotting the results if dimension is 2
function plot_2d_surface(nodes, values)
    # The following two lines are necessary in order for the plot to look nicely
    func = StepFunction(nodes,values)
    new_nodes = vcat([VoronoiNode([k/10,j*1.0]) for k in 0:10, j in 0:1], [VoronoiNode([j*1.0,k/10]) for k in 1:9, j in 0:1])
    append!(nodes,new_nodes)
    append!(values,[func(n) for n in new_nodes])
    
    x = [node[1] for node in nodes]
    y = [node[2] for node in nodes]
    
    p = Plots.surface(x, y, values, legend=false)
    xlabel!("X")
    ylabel!("Y")
    zlabel!("Values")
    title!("2D Surface Graph")
    
    display(p)
end


# Flux function passed as a parameter to HighVoronoi
function SQRA_flux(;para_i,para_j,mass_ij,normal,kwargs...) 
    # kwargs... collects all additional parameters which are not used in the current function.
    weight = norm(normal)^(-1) * mass_ij * sqrt(para_i[:κ]*para_j[:κ])
    return weight, weight
end

# RHS passed to HighVoronoi
myRHS(;para_i,mass_i,kwargs...) = mass_i * para_i[:f] 

# performs numerical calculations to solve -∇⋅(κ∇u) = RHS
function simulation(set)
    # adjust RHS for periodic domain 
    new_RHS = HighVoronoi.PeriodicFunction(set.RHS,set.domain)
    set = (neumann = x->0.0, set..., RHS=new_RHS)
    # get nodes
    nodes = nothing
    if isdefined(set,:density)
        nodes = VoronoiNodes( set.number_of_nodes;density=set.density, 
                    domain=cuboid(set.dim,periodic=[]), silence=false)
    else
        nodes = VoronoiNodes(rand(set.dim,set.number_of_nodes))
    end
    # Voronoi Geometry and integration. We could also set up VG_κ directly...
    VG_basis = VoronoiGeometry(nodes,set.domain,integrator=HighVoronoi.VI_GEOMETRY)
    VG_κ = VoronoiGeometry(VG_basis, integrator=HighVoronoi.VI_POLYGON, integrand=x->[set.κ(x),set.RHS(x)])
    vd = VoronoiData(VG_κ) # needed for volumes in the final calculations
    # set up fluxes and RHS
    vfvp = VoronoiFVProblem(VG_κ,  
                              integralfunctions = (κ = set.κ, f = set.RHS, ), 
                              fluxes = ( j1 = SQRA_flux, ),
                              rhs_functions = (F = myRHS,) )
    # define functions that can be applied as boundary conditions
    harmonic = FVevaluate_boundary(x->0.0) # turn a function into the format HighVoronoi needs
    compatibility = FVevaluate_boundary(x->set.u_exact(x))
    neumann = FVevaluate_boundary(x->set.neumann(x))
    # construct linear system from fluxes, RHS and boundary conditions
    r,c,v,f = linearVoronoiFVProblem(vfvp, flux = :j1, rhs = :F, 
            Dirichlet = set.dirichlet_boundary!=nothing ? (set.dirichlet_boundary,compatibility) : nothing, 
            Neumann = set.neumann_boundary!=nothing ? (set.neumann_boundary,neumann) : nothing)
    A = sparse(r,c,v) # a sparse matrix with rows `r`, coloumns `c` and values `v`
    # solve linear system using IterativeSolvers-Package
    solution_u = cg(A,f) # conjugate gradients
    # print out approximate L²-error between exact and numerical solutions 
    println("Approximate L²-error: ",sqrt(sum(map(k->abs2(solution_u[k]-set.u_exact(nodes[k]))*vd.volume[k],1:length(nodes)))))
    return nodes, solution_u  # return nodes and values for plotting...
end


###################################################################################

## Putting together the above pieces

###################################################################################

# create parameters
my_set = make_set(x->sin(x[1]*2*π)^2 * sin(x[2]*2*π)^2, x->1.0, 2, 
                periodic=[1], dirichlet_boundary=3, neumann_boundary=4)
                # reminder: periodic=[1] identifies boundary 1 with boundary 2
# perform simulation
nodes, values = simulation(my_set)
# plot results
my_set.dim==2 && plot_2d_surface(nodes, values)