Bim package: Difference between revisions

Jump to navigation Jump to search
258 bytes added ,  19 July 2012
no edit summary
No edit summary
No edit summary
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);
43

edits

Navigation menu