diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index f74740f31..155b34bc6 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,4 +1,4 @@
-name: ci
+name: linux-macos-windows
 on:
   pull_request:
     branches:
@@ -17,7 +17,7 @@ jobs:
         version:
           - 'lts' # Replace this with the minimum Julia version that your package supports. E.g. if your package requires Julia 1.5 or higher, change this to '1.5'.
           - '1'   # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia.
-          - 'nightly'
+#          - 'pre'
         os:
           - ubuntu-latest
           - windows-latest
@@ -31,7 +31,7 @@ jobs:
             arch: aarch64
           - os: windows-latest
             arch: aarch64
-          - os: macOS-13 
+          - os: macOS-13
             arch: aarch64
           - os: macOS-14
             arch: x64
@@ -79,4 +79,3 @@ jobs:
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
-          
\ No newline at end of file
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 5561fd28b..ecfa462d2 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,7 +4,10 @@
   - use `precs` based linear solver API, see https://github.com/SciML/LinearSolve.jl/pull/514
   - stop re-exporting ForwardDiff.value
   - try to remove type piracies
-  
+## v2.1.0 October 23, 2024
+  - Features:
+    - Add VoronoiFVMExtendableFEMBaseExt  extension allowing to use velocity fields caculated  with [ExtendableFEM](https://github.com/chmerdon/ExtendableFEM.jl) in convection-diffusion problems
+
 ## v2.0.2 October 12, 2024
   - Bugfixes
     - replace internally used `values` and `vec` methods by new `dofs`
@@ -37,7 +40,7 @@
 ### Modified
   - Sparse solver default for Float64 in all space dimensions is now  UMFPACKFactorization. For all other
     number types it is SparspakFactorization. Before it was KLU for 1D, Sparspak for 2D an UMFPACK for 3D.
-    
+
 ## v1.25.0, Sept 9, 2024
 - [`VoronoiFVM.Physics`](@ref) callbacks (`flux`, `storage`, etc.) without `data` argument are now deprecated
 
diff --git a/README.md b/README.md
index a56e9d8c0..b998198b7 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,7 @@
 VoronoiFVM.jl
 ===============
 
-[![Build status](https://github.com/j-fu/VoronoiFVM.jl/workflows/linux-macos-windows/badge.svg)](https://github.com/j-fu/VoronoiFVM.jl/actions)
+[![Build status](https://github.com/j-fu/VoronoiFVM.jl/actions/workflows/ci.yml/badge.svg?branch=master)](https://github.com/j-fu/VoronoiFVM.jl/actions/workflows/ci.yml?query=branch%3Amaster)
 [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://j-fu.github.io/VoronoiFVM.jl/stable)
 [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://j-fu.github.io/VoronoiFVM.jl/dev)
 [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3529808.svg)](https://doi.org/10.5281/zenodo.3529808)
@@ -33,7 +33,7 @@ VoronoiFVM.jl and most of these packages are  part of the meta package [PDELib.j
 - [ExtendableFEM.jl](https://github.com/chmerdon/ExtendableFEM.jl): finite element library implementing gradient robust FEM
   from the same package base by Ch. Merdon
 - [SkeelBerzins.jl](https://github.com/gregoirepourtier/SkeelBerzins.jl): a Julian variation on Matlab's `pdepe` API
-- [Trixi.jl](https://github.com/trixi-framework/Trixi.jl):  numerical simulation framework for hyperbolic conservation laws 
+- [Trixi.jl](https://github.com/trixi-framework/Trixi.jl):  numerical simulation framework for hyperbolic conservation laws
 - [GridAP.jl](https://github.com/gridap/Gridap.jl) Grid-based approximation of partial differential equations in Julia
 - [Ferrite.jl](https://github.com/Ferrite-FEM/Ferrite.jl) Finite element toolbox for Julia
 - [FinEtools.jl](https://github.com/PetrKryslUCSD/FinEtools.jl)  Finite element tools for Julia
diff --git a/test/xgrid/Example120_ThreeRegions1D.jl b/test/xgrid/Example120_ThreeRegions1D.jl
deleted file mode 100644
index 1e2708756..000000000
--- a/test/xgrid/Example120_ThreeRegions1D.jl
+++ /dev/null
@@ -1,120 +0,0 @@
-# # 120: Differing species sets in regions, 1D
-# ([source code](SOURCE_URL))
-
-module Example120_ThreeRegions1D
-
-using Printf
-using VoronoiFVM
-
-
-
-
-function main(;n=30,Plotter=nothing,plot_grid=false, verbose=false,unknown_storage=:sparse)
-    h=3.0/(n-1)
-    X=collect(0:h:3.0)
-    grid=VoronoiFVM.Grid(X)
-    cellmask!(grid,[0.0],[1.0],1)
-    cellmask!(grid,[1.0],[2.1],2)
-    cellmask!(grid,[1.9],[3.0],3)
-
-    subgrid1=subgrid(grid,[1])   
-    subgrid2=subgrid(grid,[1,2,3])
-    subgrid3=subgrid(grid,[3])
-    
-    if isplots(Plotter)&&plot_grid
-        p=Plotter.plot()
-        VoronoiFVM.plot(Plotter,grid,p=p)
-        Plotter.gui(p)
-        return
-    end
-    
-    eps=[1,1,1]
-    k=[1,1,1]
-
-    physics=VoronoiFVM.Physics(
-    num_species=3,
-    reaction=function(f,u,node)
-        if node.region==1
-            f[1]=k[1]*u[1]
-            f[2]=-k[1]*u[1]
-        elseif node.region==3
-            f[2]=k[3]*u[2]
-            f[3]=-k[3]*u[2]
-        else
-            f[1]=0
-        end
-    end,
-    
-    flux=function(f,u,edge)   
-        if edge.region==1
-            f[1]=eps[1]*(u[1,1]-u[1,2])
-            f[2]=eps[2]*(u[2,1]-u[2,2])
-        elseif edge.region==2
-            f[2]=eps[2]*(u[2,1]-u[2,2])
-        elseif edge.region==3
-            f[2]=eps[2]*(u[2,1]-u[2,2])
-            f[3]=eps[3]*(u[3,1]-u[3,2])
-        end
-    end,
-    
-    source=function(f,node)
-        if node.region==1
-           f[1]=1.0e-4*(3.0-node[1])
-        end
-    end,
-    
-    storage=function(f,u,node)
-        f.=u
-    end
-    )
-
-    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
-
-    enable_species!(sys,1,[1])
-    enable_species!(sys,2,[1,2,3])
-    enable_species!(sys,3,[3])
-
-    boundary_dirichlet!(sys,3,2,0.0)
-    
-    inival=unknowns(sys)
-    U=unknowns(sys)
-    inival.=0
-    
-    control=VoronoiFVM.NewtonControl()
-    control.verbose=verbose
-    tstep=0.01
-    time=0.0
-    istep=0
-    tend=10
-    testval=0
-    while time<tend
-        time=time+tstep
-        solve!(U,inival,sys,control=control,tstep=tstep)
-        inival.=U
-        if verbose
-            @printf("time=%g\n",time)
-        end
-        tstep*=1.0
-        istep=istep+1
-        testval=U[2,15]
-        if isplots(Plotter)
-            Plots=Plotter
-            p=Plots.plot()
-            U1=view(U[1,:],subgrid1)
-            U2=view(U[2,:],subgrid2)
-            U3=view(U[3,:],subgrid3)
-            plot(subgrid1, U1,label="spec1", color=(0.5,0,0),p=p,show=false, Plotter=Plots)
-            plot(subgrid2, U2,label="spec2", color=(0.0,0.5,0),p=p,show=false, Plotter=Plots)
-            plot(subgrid3, U3,label="spec3", color=(0.0,0.0,0.5),p=p,show=false, Plotter=Plots)
-            Plots.gui(p)
-        end
-    end
-    return testval
-end
-
-function test()
-    main(unknown_storage=:sparse) ≈ 0.00039500514567080265 &&
-        main(unknown_storage=:dense) ≈ 0.00039500514567080265
-end
-
-end
diff --git a/test/xgrid/Example121_PoissonPointCharge1D.jl b/test/xgrid/Example121_PoissonPointCharge1D.jl
deleted file mode 100644
index 1b1dfe442..000000000
--- a/test/xgrid/Example121_PoissonPointCharge1D.jl
+++ /dev/null
@@ -1,117 +0,0 @@
-#=
- # 121: 1D Poisson with point charge
- ([source code](SOURCE_URL))
-
-Solve a Poisson equation
-```math
-- \Delta u = 0
-```
-
-in $\Omega=(-1,1)$
-with a point charge $Q$ at $x=0$. 
-=#
-
-module Example121_PoissonPointCharge1D
-
-using Printf
-
-using VoronoiFVM
-
-function main(;nref=0,Plotter=nothing, verbose=false, unknown_storage=:sparse, brea=false)
-    
-    ## Create grid in (-1,1) refined around 0
-    hmax=0.2/2.0^nref
-    hmin=0.05/2.0^nref
-    X1=geomspace(-1.0,0.0, hmax,hmin)
-    X2=geomspace(0.0,1.0, hmin,hmax)
-    X=glue(X1,X2)
-    grid=VoronoiFVM.Grid(X)
-
-    ## Edit default region numbers:
-    ##   additional boundary region 3 at 0.0
-    bfacemask!(grid, [0.0],[0.0],3)
-    ## Material 1 left of 0
-    cellmask!(grid, [-1.0],[0.0],1)
-    ## Material 2 right of 0
-    cellmask!(grid, [0.0],[1.0],2)
-    
-
-    Q=0.0
-
-    function flux!(f,u,edge)
-        f[1]=u[1]-u[2]
-    end
-    function storage!(f,u,node)
-        f[1]=u[1]
-    end
-
-    ## Define boundary reaction defining charge
-    ## Note that the term  is written on  the left hand side, therefore the - sign
-    function breaction!(f,u,node)
-        if node.region==3
-            f[1]=-Q
-        end
-    end
-
-   ## Create physics
-    physics=VoronoiFVM.Physics(
-        flux=flux!,
-        storage=storage!,
-        breaction=breaction!
-    )
-
-    ## Create system
-    sys=VoronoiFVM.System(grid,physics,unknown_storage=:dense)
-
-
-    ##  put potential into both regions
-    enable_species!(sys,1,[1,2])
-
-    ## Set boundary conditions
-
-    boundary_dirichlet!(sys,1,1,1.0)
-    boundary_dirichlet!(sys,1,2,0.0)
-
-    ## Create a solution array
-    inival=unknowns(sys)
-    U=unknowns(sys)
-    inival.=0
-
-    ## Create solver control info
-    control=VoronoiFVM.NewtonControl()
-    control.verbose=verbose
-    if ispyplot(Plotter)
-        Plotter.clf()
-    end
-
-    ## Solve and plot for several values of charge
-    for q in [0.1,0.2,0.4,0.8,1.6]
-        
-        if brea
-            ## Charge in reaction term
-            Q=q
-        else
-            ## Charge as boundary condition
-            sys.boundary_values[1,3]=q
-        end
-        solve!(U,inival,sys, control=control)
-        
-        ## Plot data
-        if ispyplot(Plotter)
-            Plotter.grid()
-            coord=grid[Coordinates]
-            Plotter.plot(coord[1,:],U[1,:],label=@sprintf("Q=%.2f",q))
-            Plotter.xlabel("x")
-            Plotter.ylabel("\$\\phi\$")
-            Plotter.legend(loc="upper right")
-            Plotter.pause(1.0e-10)
-        end
-    end
-    return sum(U)
-end
-
-function test()
-    main()≈20.254591679579015
-end
-end 
-
diff --git a/test/xgrid/Example202_Laplace2D_Unstructured.jl b/test/xgrid/Example202_Laplace2D_Unstructured.jl
deleted file mode 100644
index 01b816d90..000000000
--- a/test/xgrid/Example202_Laplace2D_Unstructured.jl
+++ /dev/null
@@ -1,82 +0,0 @@
-#=
-
-# 202: 2D Laplace equation on unstructured mesh 
-([source code](SOURCE_URL))
-
-=#
-module Example202_Laplace2D_Unstructured
-
-using VoronoiFVM
-using LinearAlgebra
-
-
-# Flux function which describes the flux
-# between neighboring control volumes $\omega_k$ and $\omega_l$
-function flux!(f,u,edge)
-    f[1]=u[1,1]-u[1,2]
-end
-
-function source!(f,edge)
-    f[1]=1
-end
-
-
-function main(;Plotter=nothing, plot_grid=false,nref=0)
-
-    nspecies=1
-    ispec=1
-    nrad=10*2^nref
-    grid=VoronoiFVM.Grid(points=reduce(hcat,[ [cos(2*pi*i/nrad), sin(2*pi*i/nrad)] for i=1:nrad]),
-                     bfaces=reduce(hcat,vcat([ [i,i+1] for i=1:nrad-1],[[nrad,1]])),                  
-                     bfaceregions=ones(nrad),
-                     regionpoints=[0.0 0.0;],
-                     regionnumbers=[1],
-                     regionvolumes=[0.1*2.0^(-2*nref)])
-
-    if plot_grid
-        p=plot(Plotter=Plotter,grid)
-        return
-    end
-
-    
-    # Create a physics structure
-    physics=VoronoiFVM.Physics(num_species=nspecies,flux=flux!,source=source!)
-
-    # Create a finite volume system with dense storage of unknowns
-    sys=VoronoiFVM.System(grid,physics,unknown_storage=:dense)
-
-    # Enable species 1 in region 1
-    enable_species!(sys,ispec,[1])
-
-    # Set boundary conditions
-    # Dirichlet boundary conditions are marked by setting a corresponding value of the
-    # boundary factor
-    for i=1:num_bfaceregions(grid)
-        boundary_dirichlet!(sys,ispec,i,0.0)
-    end
-    
-    # Create & initialize array for solution and initial value
-    inival=unknowns(sys)
-    solution=unknowns(sys)
-    inival.=0
-    solution.=0
-
-    # Solve stationary problem
-    solve!(solution,inival,sys)
-
-    if Plotter!=nothing
-        plot(grid, solution[1,:], Plotter=Plotter)
-    end
-    
-    # Return test value
-    return norm(solution,Inf)
-end
-
-#
-# Called by unit test
-#
-function test()
-    main() ≈ 0.23422190956113936
-end
-
-end
diff --git a/test/xgrid/Example220_NonlinearPoisson2D_BoundarySpecies.jl b/test/xgrid/Example220_NonlinearPoisson2D_BoundarySpecies.jl
deleted file mode 100644
index 15ac0dfac..000000000
--- a/test/xgrid/Example220_NonlinearPoisson2D_BoundarySpecies.jl
+++ /dev/null
@@ -1,113 +0,0 @@
-# # 220: 2D Nonlinear Poisson with boundary reaction and boundary species
-# ([source code](SOURCE_URL))
-
-module Example220_NonlinearPoisson2D_BoundarySpecies
-
-using Printf
-using VoronoiFVM
-
-
-function main(;n=10,Plotter=nothing,verbose=false,unknown_storage=:sparse)
-    
-    
-    h=1.0/convert(Float64,n)
-    X=collect(0.0:h:1.0)
-    Y=collect(0.0:h:1.0)
-    
-
-
-    grid=VoronoiFVM.Grid(X,Y)
-    
-    
-    k=1.0
-    eps=1.0
-    physics=VoronoiFVM.Physics(
-    num_species=3,
-    breaction=function(f,u,node)
-        if  node.region==2
-            f[1]=k*(u[1]-u[3])
-            f[3]=k*(u[3]-u[1])+ k*(u[3]-u[2])
-            f[2]=k*(u[2]-u[3])
-        end
-    end,
-    
-    bstorage=function(f,u,node)
-        if  node.region==2
-            f[3]=u[3]
-        end
-    end,
-    
-    
-    flux=function(f,u,edge)
-        f[1]=eps*(u[1,1]-u[1,2])
-        f[2]=eps*(u[2,1]-u[2,2])
-    end,
-    
-    source=function(f,node)
-        x1=node[1]-0.5
-        x2=node[2]-0.5
-        f[1]=exp(-20.0*(x1^2+x2^2))
-    end,
-    
-    storage=function(f,u,node)
-        f[1]=u[1]
-        f[2]=u[2]
-    end
-    )
-    
-    sys=VoronoiFVM.System(grid,physics,unknown_storage=unknown_storage)
-    
-    enable_species!(sys,1,[1])
-    enable_species!(sys,2,[1])
-    enable_boundary_species!(sys,3,[2])
-
-    
-    function tran32!(a,b)
-        a[1]=b[2]
-    end
-    
-    bgrid2=subgrid(grid,[2],boundary=true,transform=tran32!)
-   
-    inival=unknowns(sys)
-    inival.=0.0
-    U=unknowns(sys)
-
-    eps=1.0e-2
-    
-    control=VoronoiFVM.NewtonControl()
-    control.verbose=verbose
-    control.tol_linear=1.0e-5
-    control.tol_relative=1.0e-5
-    control.max_lureuse=0
-    tstep=0.01
-    time=0.0
-    istep=0
-    u5=0
-    while time<1
-        time=time+tstep
-        solve!(U,inival,sys,control=control,tstep=tstep)
-        inival.=U
-        if verbose
-            @printf("time=%g\n",time)
-        end
-        tstep*=1.0
-        istep=istep+1
-        U_bound=view(U[3,:],bgrid2)
-        u5=U_bound[5]
-        if isplots(Plotter)
-            p1=Plotter.contourf(X,Y,reshape(U[1,:],length(X),length(Y)),levels=collect(0:0.1:0.6),clim=(0,0.6),colorbar=:right,color=:viridis,title=@sprintf("max1=%g max2=%g maxb=%g\n",maximum(U[1,:]),maximum(U[2,:]),maximum(U_bound[3,:])))
-            p2=Plotter.contourf(X,Y,reshape(U[2,:],length(X),length(Y)),levels=collect(0:0.0001:0.002),clim=(0,0.002), colorbar=:right,color=:viridis)
-            p3=Plotter.plot(grid=true,ylims=(0,0.0025))
-            plot(bgrid2,U_bound,p=p3,show=false,Plotter=Plotter)
-            p=Plotter.plot(p1,p2,p3,layout=(3,1) )
-            Plotter.gui(p)
-        end
-    end
-    return u5
-end
-
-function test()
-    main(unknown_storage=:sparse) ≈ 0.0020781361856598
-    main(unknown_storage=:dense) ≈ 0.0020781361856598
-end
-end