Difference between revisions of "Bim package"

From Octave
Jump to navigation Jump to search
Line 6: Line 6:
 
We want to solve the equation
 
We want to solve the equation
  
<math>  -\mathrm{div}\ ( \varepsilon\ \nabla u(x, y) - \nabla \varphi(x,y)\ u(x, y) ) ) + u(x, y) = 1 </math>
+
<math>  -\mathrm{div}\ ( \varepsilon\ \nabla u(x, y) - \nabla \varphi(x,y)\ u(x, y) ) ) + u(x, y) = 1 \qquad in \Omega</math>
 +
<math>  \varphi(x, y)\ =\ x + y </math>
 +
 
 +
with mixed Dirichlet / Neumann boundary conditions
  
<math> \varphi(x, y)\ =\ x + y </math>
+
<math> u(x, y) = u_d(x, y)\qquad \mbox{ on } \Gamma_D </math>
  
 +
<math> -( \varepsilon\ \nabla u(x, y) - \nabla \varphi(x,y)\ u(x, y) )  \cdot \mathbf{n} \qquad \mbox{ on } \Gamma_N</math>
  
 
<b> Create the mesh and precompute the mesh properties </b>
 
<b> Create the mesh and precompute the mesh properties </b>
Line 34: Line 38:
 
<b> Construct an initial guess</b>
 
<b> Construct an initial guess</b>
  
For a linear problem only the values at boundary nodes are actually relevant<br>
+
We need this even if our problem is linear and stationary
 +
as we are going to use the values at boundary nodes to set
 +
Dirichelet boundary conditions.
 +
 
 +
Get the node coordinates from the mesh structure
  
 
<pre>
 
<pre>
 
xu    = mesh.p(1,:).';
 
xu    = mesh.p(1,:).';
 
yu    = mesh.p(2,:).';
 
yu    = mesh.p(2,:).';
nelems = columns(mesh.t);
+
</pre>
nnodes = columns(mesh.p);
+
 
 +
 
 +
<pre>
 
uin    = 3*xu;
 
uin    = 3*xu;
 
</pre>
 
</pre>
Line 46: Line 56:
 
<b> Set the coefficients for the problem:</b>
 
<b> Set the coefficients for the problem:</b>
  
 +
Get the number of elements and nodes in the mesh
 +
 +
<pre>
 +
nelems = columns(mesh.t);
 +
nnodes = columns(mesh.p);
 +
</pre>
  
 
<pre>
 
<pre>
 
epsilon = .1;
 
epsilon = .1;
alfa    = ones(nelems,1);
+
phi     = xu+yu;
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>
 
</pre>
  
Line 62: Line 71:
  
 
<pre>
 
<pre>
AdvDiff = bim2a_advection_diffusion(mesh,alfa,gamma,eta,beta);
+
AdvDiff = bim2a_advection_diffusion(mesh, epsilon, 1, phi);
 
Mass    = bim2a_reaction(mesh,delta,zeta);
 
Mass    = bim2a_reaction(mesh,delta,zeta);
 
b      = bim2a_rhs(mesh,f,g);
 
b      = bim2a_rhs(mesh,f,g);

Revision as of 14:19, 19 July 2012

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

with mixed Dirichlet / Neumann boundary conditions

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

We need this even if our problem is linear and stationary as we are going to use the values at boundary nodes to set Dirichelet boundary conditions.

Get the node coordinates from the mesh structure

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


uin    = 3*xu;

Set the coefficients for the problem:

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