Bim package: Difference between revisions

Jump to navigation Jump to search
181 bytes added ,  19 July 2012
no edit summary
No edit summary
No edit summary
Line 55: Line 55:


<pre>
<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>
</pre>


Line 68: Line 68:
[[File:fiume_msh.png]]
[[File:fiume_msh.png]]


<b> Construct an initial guess</b>


<b> Set the coefficients for the problem:</b>
<b> Set the coefficients for the problem:</b>
Line 83: Line 82:


<pre>
<pre>
nelems = columns(mesh.t);
nelems = columns (mesh.t);
nnodes = columns(mesh.p);
nnodes = columns (mesh.p);
</pre>
</pre>


Line 95: Line 94:


<pre>
<pre>
AdvDiff = bim2a_advection_diffusion(mesh, epsilon, 1, phi);
AdvDiff = bim2a_advection_diffusion (mesh, epsilon, 1, 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);
A      = AdvDiff + Mass;
A      = AdvDiff + Mass;
</pre>
</pre>
Line 103: Line 102:
<b> To Apply Boundary Conditions, partition LHS and RHS</b>
<b> To Apply Boundary Conditions, partition LHS and RHS</b>


The tags of the sides are assigned by gmsh
The tags of the sides are assigned by gmsh we let <math> \Gamma_D </math> be composed by sides 1 and 2
and <math> \Gamma_D </math> be the rest of the boundary


<pre>
<pre>
GammaD = bim2c_unknowns_on_side(mesh, [1 2]);   ## DIRICHLET NODES LIST
GammaD = bim2c_unknowns_on_side (mesh, [1 2]);           ## DIRICHLET NODES LIST
GammaN = bim2c_unknowns_on_side(mesh, [3 4]);   ## NEUMANN NODES LIST
GammaN = bim2c_unknowns_on_side (mesh, [3 4]);           ## NEUMANN NODES LIST
Corners = setdiff(GammaD,GammaN);
GammaN = setdiff (GammaN, GammaD);
jn    = zeros(length(GammaN),1);             ## PRESCRIBED NEUMANN FLUXES
 
ud    = 3*xu;                                     ## DIRICHLET DATUM
jn    = zeros (length (GammaN),1);                     ## PRESCRIBED NEUMANN FLUXES
Ilist = setdiff(1:length(uin),union(Dlist,Nlist)); ## INTERNAL NODES LIST
ud    = 3*xu;                                             ## DIRICHLET DATUM
Omega = setdiff (1:length (uin), union (GammaD, GammaN)); ## INTERIOR NODES LIST
</pre>
</pre>




<pre>
<pre>
Add = A(Dlist,Dlist);
Add = A(GammaD, GammaD);
Adn = A(Dlist,Nlist); ## shoud be all zeros hopefully!!
Adn = A(GammaD, GammaN); ## shoud be all zeros hopefully!!
Adi = A(Dlist,Ilist);  
Adi = A(GammaD, Omega);  


And = A(Nlist,Dlist); ## shoud be all zeros hopefully!!
And = A(GammaN, GammaD); ## shoud be all zeros hopefully!!
Ann = A(Nlist,Nlist);
Ann = A(GammaN, GammaN);
Ani = A(Nlist,Ilist);  
Ani = A(GammaN, Omega);  


Aid = A(Ilist,Dlist);
Aid = A(Ilist, GammaD);
Ain = A(Ilist,Nlist);  
Ain = A(Ilist, GammaN);  
Aii = A(Ilist,Ilist);  
Aii = A(Ilist, Omega);  


bd = b(Dlist);
bd = b(GammaD);
bn = b(Nlist);  
bn = b(GammaN);  
bi = b(Ilist);
bi = b(Omega);  
 
ud = uin(Dlist);
un = uin(Nlist);
ui = uin(Ilist);  
</pre>
</pre>


Line 140: Line 137:


<pre>
<pre>
temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud ; bi-Aid*ud];
temp = [Ann Ani ; Ain Aii ] \ [ Fn+bn-And*ud(GammaD) ; bi-Aid*ud(GammaD)];
un   = temp(1:length(un));
u(Nlist)   = temp(1:numel (GammaN));
ui   = temp(length(un)+1:end);
u(Omega)   = temp(length(un)+1:end);
u(Dlist) = ud;
u(Ilist) = ui;
u(Nlist) = un;
</pre>
</pre>


<b> Compute the fluxes through Dirichlet sides</b><br>
<b> Compute the fluxes through Dirichlet sides</b><br>
<pre>
<pre>
Fd = Add * ud + Adi * ui + Adn*un - bd;
Fd = [Add Adi Adn] * u([GammaD; Omega; GammaN]) - bd;
</pre>
</pre>




<B> Compute the gradient of the solution </B><BR>
<B> Compute the gradient of the solution </B>
 
<pre>
<pre>
[gx, gy] = bim2c_pde_gradient(mesh,u);
[gx, gy] = bim2c_pde_gradient (mesh, u);
</pre>
</pre>


<B> Compute the internal Advection-Diffusion flux</B><BR>
<B> Compute the internal Advection-Diffusion flux</B>
 
<pre>
<pre>
[jxglob,jyglob] = bim2c_global_flux(mesh,u,alfa,gamma,eta,beta);
[jxglob, jyglob] = bim2c_global_flux (mesh, u, epsilon, 1, phi);
</pre>
</pre>


<B> Save data for later visualization</B><BR>
<B> Save data for later visualization</B><BR>
The resut can be exported to vtk format to visualize with [[http://www.paraview.org|paraview]]
or [[https://wci.llnl.gov/codes/visit/|visit]]
<pre>
<pre>
fpl_dx_write_field("dxdata",mesh,[gx; gy]',"Gradient",1,2,1);
fpl_vtk_write_field ("vtkdata", mesh, {}, {[gx; gy]', "Gradient"}, 1);
fpl_vtk_write_field ("vtkdata", mesh, {}, {[gx; gy]', "Gradient"}, 1);
</pre>
</pre>
43

edits

Navigation menu