Changes

Jump to navigation Jump to search
939 bytes added ,  08:40, 20 July 2012
no edit summary
to see the mesh you can use functions from the [[fpl_package|fpl package]]
{{Code|Plot mesh in the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
pdemesh (mesh.p, mesh.e, mesh.t)
view (2)
</presyntaxhighlight> }}
[[File:fiume_msh.png]]
Get the node coordinates from the mesh structure
{{Code|Get mesh coordinates in the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
xu = mesh.p(1,:).';
yu = mesh.p(2,:).';
</presyntaxhighlight>}}
Get the number of elements and nodes in the mesh
{{Code|Get number of elements in the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
nelems = columns (mesh.t);
nnodes = columns (mesh.p);
</presyntaxhighlight>}}
<pre>
<b> Construct the discretized operators</b>
{{Code|Discretized operators for the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi);
Mass = bim2a_reaction (mesh, 1, 1);
b = bim2a_rhs (mesh,f,g);
A = AdvDiff + Mass;
</presyntaxhighlight>}}
<b> To Apply Boundary Conditions, partition LHS and RHS</b>
and <math> \Gamma_D </math> be the rest of the boundary
{{Code|Boundary conditions for the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
GammaD = bim2c_unknowns_on_side (mesh, [1 2]); ## DIRICHLET NODES LIST
GammaN = bim2c_unknowns_on_side (mesh, [3 4]); ## NEUMANN NODES LIST
ud = 3*xu; ## DIRICHLET DATUM
Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST
</pre>
 
<pre>
Add = A(GammaD, GammaD);
Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!!
bn = b(GammaN);
bi = b(Omega);
</presyntaxhighlight>}} 
<B> Solve for the displacements</B>
{{Code|Displacement in the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)];
u = ud;
u(GammaN) = temp(1:numel (GammaN));
u(Omega) = temp(length(GammaN)+1:end);
</presyntaxhighlight>}}
<b> Compute the fluxes through Dirichlet sides</b><br>
{{Code|Fluxes at sides in the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd;
</presyntaxhighlight>}}
<B> Compute the gradient of the solution </B>
{{Code|Gradient of solution in the 2D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
[gx, gy] = bim2c_pde_gradient (mesh, u);
</presyntaxhighlight>}}
<B> Compute the internal Advection-Diffusion flux</B>
<math> u = \exp (- \left(\frac{x-.2}{.2}\right)^2 - \left(\frac{y-.2}{.2}\right)^2 - \left(\frac{z-.2}{.2}\right)^2)</math>
{{Code|Define the 3D problem|<presyntaxhighlight lang="octave" style="font-size:13px">
pkg load bim
fpl_vtk_write_field (name, msh, {U(ii,:)', 'u'}, {}, 1);
endfor
</presyntaxhighlight>}}
[http://youtu.be/2E6Z_AcV8CQ This is a video] showing the .3 isosurface of the solution.
657

edits

Navigation menu