Editing Bim package
Jump to navigation
Jump to search
The edit can be undone. Please check the comparison below to verify that this is what you want to do, and then publish the changes below to finish undoing the edit.
Latest revision | Your text | ||
Line 1: | Line 1: | ||
This is a short example on how to use <tt>bim</tt> to solve a DAR problem.<br> | |||
The data for this example can be found in the doc directory inside the | |||
bim installation directory. | |||
This is a short example on how to use <tt>bim</tt> to solve a | |||
<b> Create the mesh and precompute the mesh properties </b> | <b> Create the mesh and precompute the mesh properties </b> | ||
The geometry of the domain was created using gmsh and is stored in the file <tt>fiume.geo</tt> | |||
created with [http://gmsh.geuz.org gmsh] | |||
[[File:fiume.png]] | |||
<pre> | <pre> | ||
[mesh] = msh2m_gmsh("fiume","scale",1,"clscale",.1); | |||
[mesh] = bim2c_mesh_properties(mesh); | |||
</pre> | </pre> | ||
<b> Construct an initial guess</b> | |||
For a linear problem only the values at boundary nodes are actually relevant<br> | |||
<pre> | |||
xu = mesh.p(1,:).'; | xu = mesh.p(1,:).'; | ||
yu = mesh.p(2,:).'; | yu = mesh.p(2,:).'; | ||
</ | nelems = columns(mesh.t); | ||
nnodes = columns(mesh.p); | |||
uin = 3*xu; | |||
</pre> | |||
<b> Set the coefficients for the problem:</b> | |||
<math> -div ( \alpha \gamma ( \eta \nabla u - \beta u ) )+ \delta \zeta u = f g </math> | |||
<pre> | |||
epsilon = .1; | 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); | |||
</pre> | |||
<b> Construct the discretized operators</b> | <b> Construct the discretized operators</b> | ||
<pre> | |||
AdvDiff = bim2a_advection_diffusion (mesh, | AdvDiff = bim2a_advection_diffusion(mesh,alfa,gamma,eta,beta); | ||
Mass = bim2a_reaction (mesh, | Mass = bim2a_reaction(mesh,delta,zeta); | ||
b = bim2a_rhs (mesh,f,g); | b = bim2a_rhs(mesh,f,g); | ||
A = AdvDiff + Mass; | A = AdvDiff + Mass; | ||
</ | </pre> | ||
<b> To Apply Boundary Conditions, partition LHS and RHS</b> | <b> To Apply Boundary Conditions, partition LHS and RHS</b> | ||
The tags of the sides are assigned by gmsh | The tags of the sides are assigned by gmsh | ||
<pre> | |||
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 | |||
</pre> | |||
Add = A( | <pre> | ||
Adn = A( | Add = A(Dlist,Dlist); | ||
Adi = A( | Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!! | ||
Adi = A(Dlist,Ilist); | |||
And = A( | And = A(Nlist,Dlist); ## shoud be all zeros hopefully!! | ||
Ann = A( | Ann = A(Nlist,Nlist); | ||
Ani = A( | Ani = A(Nlist,Ilist); | ||
Aid = A( | Aid = A(Ilist,Dlist); | ||
Ain = A( | Ain = A(Ilist,Nlist); | ||
Aii = A( | Aii = A(Ilist,Ilist); | ||
bd = b( | bd = b(Dlist); | ||
bn = b( | bn = b(Nlist); | ||
bi = b( | bi = b(Ilist); | ||
ud = uin(Dlist); | |||
un = uin(Nlist); | |||
ui = uin(Ilist); | |||
</pre> | |||
<B> Solve for the | <B> Solve for the displacements</B> | ||
<pre> | |||
temp = [Ann Ani ; Ain Aii ] \ [ | 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; | |||
</pre> | |||
<b> Compute the fluxes through Dirichlet sides</b><br> | <b> Compute the fluxes through Dirichlet sides</b><br> | ||
<pre> | |||
Fd = Add * ud + Adi * ui + Adn*un - bd; | |||
</pre> | |||
<B> Compute the gradient of the solution </B><BR> | |||
<pre> | |||
[gx, gy] = bim2c_pde_gradient(mesh,u); | |||
</pre> | |||
<B> Compute the | <B> Compute the internal Advection-Diffusion flux</B><BR> | ||
<pre> | |||
[jxglob,jyglob] = bim2c_global_flux(mesh,u,alfa,gamma,eta,beta); | |||
</pre> | |||
<B> Save data for later visualization</B><BR> | |||
<pre> | |||
fpl_dx_write_field("dxdata",mesh,[gx; gy]',"Gradient",1,2,1); | |||
fpl_vtk_write_field ("vtkdata", mesh, {}, {[gx; gy]', "Gradient"}, 1); | |||
<B> | </pre> | ||
</ | |||
</ | |||
[[Category: | [[Category:OctaveForge]] | ||
[[Category:Packages]] |