This is a short example on how to use BIM to solve a DAR problem.
The data for this example can be found in the directory BIM/examples/FIUME

Create the mesh and precompute the mesh properties
The geometry of the domain was created using gmsh and is stored in the file fiume.geo

[mesh] = MSH2Mgmsh('fiume',1,.1);
[mesh] = BIM2Cmeshproperties(mesh);

Construct an initial guess
For a linear problem only the values at internal nodes are actually relevant

[xu, yu] = BIM2Cunknowncoord(mesh);
nelem = size(mesh.t,2);
nnodi = size(mesh.p,2);
uin   = 3*xu;

set the coefficients for the problem: -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g

epsilon = .1;
alfa    = ones(nelem,1);
gamma   = ones(nnodi,1);
eta     = epsilon*ones(nnodi,1);
beta    = xu+yu;
delta   = ones(nelem,1);
zeta    = ones(nnodi,1);
f       = ones(nnodi,1);
g       = ones(nelem,1);

Construct the discretized operators

AdvDiff = BIM2Aadvdiff(mesh,alfa,gamma,eta,beta);
Mass    = BIM2Areaction(mesh,delta,zeta);
b       = BIM2Arhs(mesh,f,g);
A       = AdvDiff + Mass;

To Apply Boundary Conditions, partition LHS and RHS
The tags of the sides are assigned by gmsh

Dlist = BIM2Cunknownsonside(mesh, [8 18]); 	   ## DIRICHLET NODES LIST
Nlist = BIM2Cunknownsonside(mesh, [23 24]); 	   ## NEUMANN NODES LIST
Nlist = setdiff(Nlist,Dlist);
Fn    = zeros(length(Nlist),1);           	   ## PRESCRIBED NEUMANN FLUXES
Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## INTERNAL NODES LIST

Add = A(Dlist,Dlist);
Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!!
Adi = A(Dlist,Ilist); 

And = A(Nlist,Dlist); ## shoud be all zeros hopefully!!
Ann = A(Nlist,Nlist);
Ani = A(Nlist,Ilist); 

Aid = A(Ilist,Dlist);
Ain = A(Ilist,Nlist); 
Aii = A(Ilist,Ilist); 

bd = b(Dlist);
bn = b(Nlist); 
bi = b(Ilist); 

ud = uin(Dlist);
un = uin(Nlist); 
ui = uin(Ilist);

Solve for the displacements

temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud];
un   = temp(1:length(un));
ui   = temp(length(un)+1:end);

Compute the fluxes through Dirichlet sides

Fd = Add * ud + Adi * ui + Adn*un - bd;

Compute the gradient of the solution

[gx, gy] = BIM2Cpdegrad(mesh,u);

Compute the internal Advection-Diffusion flux

[jxglob,jyglob] = BIM2Cglobalflux(mesh,u,alfa,gamma,eta,beta);

Visulaization

FPL2pdesurf(mesh,u);
FPL2quiver(mesh,jxglob,jyglob,"sample_density","1000");
FPL2quiver(mesh,gx,gy,"sample_density","1000");