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.
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] = msh2m_gmsh("fiume","scale",1,"clscale",.1); [mesh] = bim2c_mesh_properties(mesh);
Construct an initial guess
For a linear problem only the values at boundary nodes are actually relevant
xu = mesh.p(1,:).'; yu = mesh.p(2,:).'; nelems = columns(mesh.t); nnodes = columns(mesh.p); uin = 3*xu;
Set the coefficients for the problem:
epsilon = .1; alfa = ones(nelems,1); gamma = ones(nnodes,1); eta = epsilon*ones(nnodes,1); beta = xu+yu; delta = ones(nelems,1); zeta = ones(nnodes,1); f = ones(nelems,1); g = ones(nnodes,1);
Construct the discretized operators
AdvDiff = bim2a_advection_diffusion(mesh,alfa,gamma,eta,beta); 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);