Bim package: Difference between revisions

678 bytes added ,  19 July 2012
no edit summary
No edit summary
No edit summary
Line 13: Line 13:
<math> u(x, y) = u_d(x, y)\qquad \mbox{ on } \Gamma_D </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>
<math> -( \varepsilon\ \nabla u(x, y) - \nabla \varphi(x,y)\ u(x, y) )  \cdot \mathbf{n} = j_N(x, y)\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>


The geometry of the domain was created using gmsh and is stored in the file <tt>fiume.geo</tt>
To define the geometry of the domain we can use [http://gmsh.geuz.org gmsh].
created with [http://gmsh.geuz.org gmsh]
 
the following gmsh input
 
<pre>
Point (1)  = {0, 0, 0, 0.1};
Point (2)  = {1, 1, 0, 0.1};
Point (3)  = {1, 0.9, 0, 0.1};
Point (4)  = {0, 0.1, 0, 0.1};
Point (5) = {0.3,0.1,-0,0.1};
Point (6) = {0.4,0.4,-0,0.1};
Point (7) = {0.5,0.6,0,0.1};
Point (8) = {0.6,0.9,0,0.1};
Point (9) = {0.8,0.8,0,0.1};
Point (10) = {0.2,0.2,-0,0.1};
Point (11) = {0.3,0.5,0,0.1};
Point (12) = {0.4,0.7,0,0.1};
Point (13) = {0.5,1,0,0.1};
Point (14) = {0.8,0.9,0,0.1};
 
Line (1)  = {3, 2};
Line (2) = {4, 1};
 
CatmullRom(3) = {1,5,6,7,8,9,3};
CatmullRom(4) = {4,10,11,12,13,14,2};
Line Loop(15) = {3,1,-4,2};
Plane Surface(16) = {15};
</pre>
 
will produce the geometry below


[[File:fiume.png]]
[[File:fiume.png]]
we need to load the mesh into Octave and precompute mesh properties
check out the tutorial for the [[msh_package|msh package]] for info
on the mesh structure


<pre>
<pre>
Line 27: Line 59:
</pre>
</pre>


to see the mesh you can use functions from the [fpl] package
to see the mesh you can use functions from the [[fpl_package|fpl package]]


<pre>
<pre>
Line 38: Line 70:
<b> Construct an initial guess</b>
<b> Construct an initial guess</b>


We need this even if our problem is linear and stationary
<b> Set the coefficients for the problem:</b>
as we are going to use the values at boundary nodes to set
Dirichelet boundary conditions.


Get the node coordinates from the mesh structure
Get the node coordinates from the mesh structure
Line 49: Line 79:
</pre>
</pre>


<pre>
uin    = 3*xu;
</pre>
<b> Set the coefficients for the problem:</b>


Get the number of elements and nodes in the mesh
Get the number of elements and nodes in the mesh
Line 65: Line 89:
<pre>
<pre>
epsilon = .1;
epsilon = .1;
phi    = xu+yu;
phi    = xu + yu;
</pre>
</pre>


Line 82: Line 106:


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

edits