Bim package: Difference between revisions

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

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

view (2)

[[File:fiume_msh.png]]

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

<pre>

<b> Construct the discretized operators</b>

AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi);

Mass    = bim2a_reaction (mesh, 1, 1);

b      = bim2a_rhs (mesh,f,g);

A      = AdvDiff + Mass;

<b> To Apply Boundary Conditions, partition LHS and RHS</b>

and <math> \Gamma_D </math> 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

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!!

bn = b(GammaN);  

bi = b(Omega);  

<B> Solve for the displacements</B>

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

<b> Compute the fluxes through Dirichlet sides</b><br>

jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd;

<B> Compute the gradient of the solution </B>

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

<B> Compute the internal Advection-Diffusion flux</B>

<math> u = \exp (- \left(\frac{x-.2}{.2}\right)^2 - \left(\frac{y-.2}{.2}\right)^2 - \left(\frac{z-.2}{.2}\right)^2)</math>

pkg load bim

   fpl_vtk_write_field (name, msh, {U(ii,:)', 'u'}, {}, 1);

endfor

[http://youtu.be/2E6Z_AcV8CQ This is a video] showing the .3 isosurface of the solution.