# Bim package

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 doc directory inside the
bim installation directory.

We want to solve the equation

with mixed Dirichlet / Neumann boundary conditions

** 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`
created with gmsh

[mesh] = msh2m_gmsh("fiume","scale",1,"clscale",.1); [mesh] = bim2c_mesh_properties(mesh);

to see the mesh you can use functions from the [fpl] package

pdemesh (mesh.p, mesh.e, mesh.t) view (2)

** Construct an initial guess**

We need this even if our problem is linear and stationary as we are going to use the values at boundary nodes to set Dirichelet boundary conditions.

Get the node coordinates from the mesh structure

xu = mesh.p(1,:).'; yu = mesh.p(2,:).';

uin = 3*xu;

** Set the coefficients for the problem:**

Get the number of elements and nodes in the mesh

nelems = columns(mesh.t); nnodes = columns(mesh.p);

epsilon = .1; phi = xu+yu;

** Construct the discretized operators**

AdvDiff = bim2a_advection_diffusion(mesh, epsilon, 1, phi); Mass = bim2a_reaction(mesh,delta,zeta); b = bim2a_rhs(mesh,f,g); A = AdvDiff + Mass;

** To Apply Boundary Conditions, partition LHS and RHS**

The tags of the sides are assigned by gmsh

Dlist = bim2c_unknowns_on_side(mesh, [8 18]); ## DIRICHLET NODES LIST Nlist = bim2c_unknowns_on_side(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); u(Dlist) = ud; u(Ilist) = ui; u(Nlist) = un;

** Compute the fluxes through Dirichlet sides**

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

** Compute the gradient of the solution **

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

** Compute the internal Advection-Diffusion flux**

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

** Save data for later visualization**

fpl_dx_write_field("dxdata",mesh,[gx; gy]',"Gradient",1,2,1); fpl_vtk_write_field ("vtkdata", mesh, {}, {[gx; gy]', "Gradient"}, 1);