43
edits
Bim package: Difference between revisions
→2D Diffusion Advection Reaction example
No edit summary |
|||
| Line 65: | Line 65: | ||
to see the mesh you can use functions from the [[fpl_package|fpl package]] | to see the mesh you can use functions from the [[fpl_package|fpl package]] | ||
{{Code| | {{Code|Visualizing the mesh for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
pdemesh (mesh.p, mesh.e, mesh.t) | pdemesh (mesh.p, mesh.e, mesh.t) | ||
view (2) | view (2) | ||
| Line 80: | Line 80: | ||
xu = mesh.p(1,:).'; | xu = mesh.p(1,:).'; | ||
yu = mesh.p(2,:).'; | yu = mesh.p(2,:).'; | ||
</syntaxhighlight>}} | </syntaxhighlight> | ||
}} | |||
| Line 90: | Line 91: | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
< | {{Code|Set value of coefficients for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
epsilon = .1; | epsilon = .1; | ||
phi = xu + yu; | phi = xu + yu; | ||
</ | </syntaxhighlight>}} | ||
<b> Construct the discretized operators</b> | <b> Construct the discretized operators</b> | ||
| Line 136: | Line 137: | ||
<B> Solve for the | <B> Solve for the tracer density</B> | ||
{{Code| | {{Code|Compute solution of the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; | temp = [Ann Ani ; Ain Aii ] \ [ jn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; | ||
u = ud; | u = ud; | ||
| Line 147: | Line 148: | ||
<b> Compute the fluxes through Dirichlet sides</b><br> | <b> Compute the fluxes through Dirichlet sides</b><br> | ||
{{Code| | {{Code|Boundary fluxes in the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd; | jd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd; | ||
</syntaxhighlight>}} | </syntaxhighlight>}} | ||
| Line 160: | Line 161: | ||
<B> Compute the internal Advection-Diffusion flux</B> | <B> Compute the internal Advection-Diffusion flux</B> | ||
< | {{Code|Total flux for the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
[jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi); | [jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon*ones(nelems, 1), ones(nnodes, 1), ones(nnodes, 1), phi); | ||
</ | </syntaxhighlight>}} | ||
<B> Export data to VTK format</B> | <B> Export data to VTK format</B> | ||
| Line 169: | Line 170: | ||
or [[https://wci.llnl.gov/codes/visit/|visit]] | or [[https://wci.llnl.gov/codes/visit/|visit]] | ||
< | {{Code|Export the solution of the 2D problem to vtk|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1); | fpl_vtk_write_field ("vtkdata", mesh, {u, "Solution"}, {[gx; gy]', "Gradient"}, 1); | ||
</ | </syntaxhighlight>}} | ||
you can also plot your data directly in Octave using <code> pdesurf </code> | you can also plot your data directly in Octave using <code> pdesurf </code> | ||
< | {{Code|Rubbersheet visualization of the solution of the 2D problem|<syntaxhighlight lang="octave" style="font-size:13px"> | ||
pdesurf (mesh.p, mesh.t, u) | pdesurf (mesh.p, mesh.t, u) | ||
</ | </syntaxhighlight>}} | ||
it will look like this | it will look like this | ||