43
edits
Bim package: Difference between revisions
Jump to navigation
Jump to search
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> | <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> | ||
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,:).'; | ||
</pre> | |||
<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; | ||
phi = xu+yu; | |||
</pre> | </pre> | ||
| Line 62: | Line 71: | ||
<pre> | <pre> | ||
AdvDiff = bim2a_advection_diffusion(mesh, | 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); | ||