Editing Bim package
Jump to navigation
Jump to search
The edit can be undone. Please check the comparison below to verify that this is what you want to do, and then publish the changes below to finish undoing the edit.
Latest revision | Your text | ||
Line 1: | Line 1: | ||
This is a short example on how to use <tt>bim</tt> to solve a DAR problem.<br> | |||
The data for this example can be found in the doc directory inside the | |||
bim installation directory. | |||
This is a short example on how to use <tt>bim</tt> to solve a | |||
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 \qquad | <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> | <math> \varphi(x, y)\ =\ x + y </math> | ||
Line 58: | Line 54: | ||
on the mesh structure | on the mesh structure | ||
<pre> | |||
[mesh] = msh2m_gmsh ("fiume","scale",1,"clscale",.1); | [mesh] = msh2m_gmsh ("fiume","scale",1,"clscale",.1); | ||
[mesh] = bim2c_mesh_properties (mesh); | [mesh] = bim2c_mesh_properties (mesh); | ||
</ | </pre> | ||
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]] | ||
<pre> | |||
pdemesh (mesh.p, mesh.e, mesh.t) | pdemesh (mesh.p, mesh.e, mesh.t) | ||
view (2) | view (2) | ||
</ | </pre> | ||
[[File:fiume_msh.png]] | [[File:fiume_msh.png]] | ||
Line 77: | Line 73: | ||
Get the node coordinates from the mesh structure | Get the node coordinates from the mesh structure | ||
<pre> | |||
xu = mesh.p(1,:).'; | xu = mesh.p(1,:).'; | ||
yu = mesh.p(2,:).'; | yu = mesh.p(2,:).'; | ||
</ | </pre> | ||
Get the number of elements and nodes in the mesh | Get the number of elements and nodes in the mesh | ||
<pre> | |||
nelems = columns (mesh.t); | nelems = columns (mesh.t); | ||
nnodes = columns (mesh.p); | nnodes = columns (mesh.p); | ||
</ | </pre> | ||
<pre> | |||
epsilon = .1; | epsilon = .1; | ||
phi = xu + yu; | phi = xu + yu; | ||
</ | </pre> | ||
<b> Construct the discretized operators</b> | <b> Construct the discretized operators</b> | ||
<pre> | |||
AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi); | AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 1, phi); | ||
Mass = bim2a_reaction (mesh, 1, 1); | Mass = bim2a_reaction (mesh, 1, 1); | ||
b = bim2a_rhs (mesh,f,g); | b = bim2a_rhs (mesh,f,g); | ||
A = AdvDiff + Mass; | A = AdvDiff + Mass; | ||
</ | </pre> | ||
<b> To Apply Boundary Conditions, partition LHS and RHS</b> | <b> To Apply Boundary Conditions, partition LHS and RHS</b> | ||
Line 110: | Line 105: | ||
and <math> \Gamma_D </math> be the rest of the boundary | and <math> \Gamma_D </math> be the rest of the boundary | ||
<pre> | |||
GammaD = bim2c_unknowns_on_side (mesh, [1 2]); | GammaD = bim2c_unknowns_on_side (mesh, [1 2]); ## DIRICHLET NODES LIST | ||
GammaN = bim2c_unknowns_on_side (mesh, [3 4]); | GammaN = bim2c_unknowns_on_side (mesh, [3 4]); ## NEUMANN NODES LIST | ||
GammaN = setdiff (GammaN, GammaD); | GammaN = setdiff (GammaN, GammaD); | ||
jn = zeros (length (GammaN),1); | jn = zeros (length (GammaN),1); ## PRESCRIBED NEUMANN FLUXES | ||
ud = 3*xu; | ud = 3*xu; ## DIRICHLET DATUM | ||
Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST | Omega = setdiff (1:nnodes, union (GammaD, GammaN)); ## INTERIOR NODES LIST | ||
</pre> | |||
<pre> | |||
Add = A(GammaD, GammaD); | Add = A(GammaD, GammaD); | ||
Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!! | Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!! | ||
Line 127: | Line 125: | ||
Ani = A(GammaN, Omega); | Ani = A(GammaN, Omega); | ||
Aid = A( | Aid = A(Ilist, GammaD); | ||
Ain = A( | Ain = A(Ilist, GammaN); | ||
Aii = A( | Aii = A(Ilist, Omega); | ||
bd = b(GammaD); | bd = b(GammaD); | ||
bn = b(GammaN); | bn = b(GammaN); | ||
bi = b(Omega); | bi = b(Omega); | ||
</ | </pre> | ||
<B> Solve for the displacements</B> | |||
< | <pre> | ||
temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)]; | |||
u(Nlist) = temp(1:numel (GammaN)); | |||
temp = [Ann Ani ; Ain Aii ] \ [ | u(Omega) = temp(length(un)+1:end); | ||
</pre> | |||
u( | |||
u(Omega) = temp(length( | |||
</ | |||
<b> Compute the fluxes through Dirichlet sides</b><br> | <b> Compute the fluxes through Dirichlet sides</b><br> | ||
<pre> | |||
Fd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd; | |||
</ | </pre> | ||
<B> Compute the gradient of the solution </B> | <B> Compute the gradient of the solution </B> | ||
<pre> | |||
[gx, gy] = bim2c_pde_gradient (mesh, u); | [gx, gy] = bim2c_pde_gradient (mesh, u); | ||
</ | </pre> | ||
<B> Compute the internal Advection-Diffusion flux</B> | <B> Compute the internal Advection-Diffusion flux</B> | ||
<pre> | |||
[jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon | [jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon, 1, phi); | ||
</ | </pre> | ||
<B> Export data to VTK format</B> | <B> Export data to VTK format</B> | ||
Line 169: | Line 165: | ||
The resut can be exported to vtk format to visualize with [[http://www.paraview.org|paraview]] | The resut can be exported to vtk format to visualize with [[http://www.paraview.org|paraview]] | ||
or [[https://wci.llnl.gov/codes/visit/|visit]] | or [[https://wci.llnl.gov/codes/visit/|visit]] | ||
<pre> | |||
fpl_vtk_write_field ("vtkdata", mesh, {}, {[gx; gy]', "Gradient"}, 1); | |||
</pre> | |||
[[Category:OctaveForge]] | |||
[[Category:Packages]] | |||
[[ | |||
[[Category: |