Bim package

From Octave
Revision as of 16:31, 19 July 2012 by Carlo (talk | contribs)
Jump to navigation Jump to search
The printable version is no longer supported and may have rendering errors. Please update your browser bookmarks and please use the default browser print function instead.

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


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

Fiume.png

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

Fiume msh.png

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