# Difference between revisions of "Bim package"

Line 3: | Line 3: | ||

This is a short example on how to use <tt>bim</tt> to solve a 2D Diffusion Advection Reaction problem. | This is a short example on how to use <tt>bim</tt> to solve a 2D Diffusion Advection Reaction problem. | ||

The coplete code for this example can is on [[Agora]] at this link | The coplete code for this example can is on [[Agora]] at this link | ||

− | [[http://agora.octave.org/snippet/1bqV| | + | [[http://agora.octave.org/snippet/1bqV|link]]. |

We want to solve the equation | We want to solve the equation |

## Revision as of 05:44, 20 July 2012

## 2D Diffusion Advection Reaction example

This is a short example on how to use `bim` to solve a 2D Diffusion Advection Reaction problem.
The coplete code for this example can is on Agora at this link
[[1]].

We want to solve the equation

with mixed Dirichlet / Neumann boundary conditions

** Create the mesh and precompute the mesh properties **

To define the geometry of the domain we can use gmsh.

the following gmsh input

Point (1) = {0, 0, 0, 0.1}; Point (2) = {1, 1, 0, 0.1}; Point (3) = {1, 0.9, 0, 0.1}; Point (4) = {0, 0.1, 0, 0.1}; Point (5) = {0.3,0.1,-0,0.1}; Point (6) = {0.4,0.4,-0,0.1}; Point (7) = {0.5,0.6,0,0.1}; Point (8) = {0.6,0.9,0,0.1}; Point (9) = {0.8,0.8,0,0.1}; Point (10) = {0.2,0.2,-0,0.1}; Point (11) = {0.3,0.5,0,0.1}; Point (12) = {0.4,0.7,0,0.1}; Point (13) = {0.5,1,0,0.1}; Point (14) = {0.8,0.9,0,0.1}; Line (1) = {3, 2}; Line (2) = {4, 1}; CatmullRom(3) = {1,5,6,7,8,9,3}; CatmullRom(4) = {4,10,11,12,13,14,2}; Line Loop(15) = {3,1,-4,2}; Plane Surface(16) = {15};

will produce the geometry below

we need to load the mesh into Octave and precompute mesh properties check out the tutorial for the msh package for info on the mesh structure

[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)

** Set the coefficients for the problem:**

Get the node coordinates from the mesh structure

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

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, 1, phi); Mass = bim2a_reaction (mesh, 1, 1); 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 we let be composed by sides 1 and 2 and be the rest of the boundary

GammaD = bim2c_unknowns_on_side (mesh, [1 2]); ## DIRICHLET NODES LIST GammaN = bim2c_unknowns_on_side (mesh, [3 4]); ## NEUMANN NODES LIST GammaN = setdiff (GammaN, GammaD); jn = zeros (length (GammaN),1); ## PRESCRIBED NEUMANN FLUXES ud = 3*xu; ## DIRICHLET DATUM Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST

Add = A(GammaD, GammaD); Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!! Adi = A(GammaD, Omega); And = A(GammaN, GammaD); ## shoud be all zeros hopefully!! Ann = A(GammaN, GammaN); Ani = A(GammaN, Omega); Aid = A(Omega, GammaD); Ain = A(Omega, GammaN); Aii = A(Omega, Omega); bd = b(GammaD); bn = b(GammaN); bi = b(Omega);

** Solve for the displacements**

temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; u = ud; u(GammaN) = temp(1:numel (GammaN)); u(Omega) = temp(length(GammaN)+1:end);

** Compute the fluxes through Dirichlet sides**

jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - 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, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi);

** Export data to VTK format**

The resut can be exported to vtk format to visualize with [[2]] or [[3]]

fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1);

you can also plot your data directly in Octave using ` pdesurf `

pdesurf (mesh.p, mesh.t, u)

it will look like this