Update _initialize! to reflect new boundary condition API

[Da-Be-Ru]

As of the last API update, all boundary condition specifications can be pooled into one function and passed to the System constructor using the bcondition keyword argument.

However, the old boundary condition handling using boundary_factor and boundary_value is of course still present not only in eval_and_assemble! among

VoronoiFVM.jl/src/vfvm_solver.jl

Line 561 in 70739a1

F[ispec,K]+=boundary_factor*(U[ispec,K]-boundary_value)

but also in _initialize! e.g.

VoronoiFVM.jl/src/vfvm_system.jl

Line 931 in 70739a1

if system.boundary_factors[ispec,ibreg]≈ Dirichlet

which is invoked in every solve! call through eval_and_assemble!.

It might be sensible to commit fully to the new boundary condition handling (which pleases me very much btw, thank you for that!) because fields like boundary_factor which are assigned once within the System construcutor can lead to problems for some narrow applications where grid-dependent boundary conditions are needed.

---

Example:

I want to evaluate a convection-diffusion system on a simple simplexgrid, but set Neumann boundary conditions on one boundary face of one cell of the grid. If I want to cycle through different boundary faces in some loop, I would set up the system once where I assign a new bregion, for a 2D rectangle say ibreg=5, and specify my Neumann condition among others in some bcond function:

function bconditions(y,u,bnode)
       
        boundary_neumann!(y,u,bnode,species=1,region=5,value=1)

        # some other bconds...
        
end 

If I then make a solve! call for this system somewhere else, VoronoiFVM protests since boundary_factor will be out of bounds, having been assigned fixed when the System constructor was called.

I currently avoid this behaviour by simply setting one entry in grid[BFaceRegions] to 5 before calling the constructor and resetting it afterward.
I am sure this is not how it is intended to be used. I have not tested this further, but if the _initialize! functions were to simply invoke the breaction, being an alias for bcondition, and boundary_factor, boundary_value were made obsolete, that would at least avoid the bounds error here and make the code more consistent.

---

Apologies for the wall of text 🙃