|
Octave-Forge - Extra packages for GNU Octave |
| Home · Packages · Developers · Documentation · FAQ · Bugs · Mailing Lists · Links · Code |
|
|
Octave-Forge - Extra packages for GNU Octave |
| Home · Packages · Developers · Documentation · FAQ · Bugs · Mailing Lists · Links · Code |
Build the Scharfetter-Gummel stabilized OSC stiffness matrix for a diffusion-advection problem.
For details on the Orthogonal Subdomain Collocation (OSC) method see: M.Putti and C.Cordes, SIAM J.SCI.COMPUT. Vol.19(4), pp.1154-1168, 1998.
The equation taken into account is:
- div (alpha ( grad (u) - grad (v) u)) = f
where v is a piecewise linear continuous scalar functions and alpha is a piecewise constant scalar function.
See also: bim3a_rhs, bim3a_osc_laplacian, bim3a_reaction, bim3a_laplacian, bim3c_mesh_properties.
The following code
gmsh_input = [["Point(1) = {0, 0, 0, .1}; \n"], ...
["Point(2) = {1, 0, 0, .1}; \n"], ...
["Point(3) = {0, -.3, 0, .1}; \n"], ...
["Point(4) = {0, +.3, 0, .1}; \n"], ...
["Point(5) = {1, -.3, 0, .1}; \n"], ...
["Point(6) = {1, 0.3, 0, .1}; \n"], ...
["Point(7) = {0, 0, -.3, .1}; \n"], ...
["Point(8) = {0, 0, +.3, .1}; \n"], ...
["Point(9) = {1, 0, -.3, .1}; \n"], ...
["Point(10) = {1, 0, 0.3, .1}; \n"], ...
["Circle(1) = {4, 1, 7}; \n"], ...
["Circle(2) = {7, 1, 3}; \n"], ...
["Circle(3) = {3, 1, 8}; \n"], ...
["Circle(4) = {8, 1, 4}; \n"], ...
["Circle(5) = {6, 2, 9}; \n"], ...
["Circle(6) = {9, 2, 5}; \n"], ...
["Circle(7) = {5, 2, 10}; \n"], ...
["Circle(8) = {10, 2, 6}; \n"], ...
["Line(9) = {4, 6}; \n"], ...
["Line(10) = {3, 5}; \n"], ...
["Line(11) = {8, 10}; \n"], ...
["Line(12) = {7, 9}; \n"], ...
["Line Loop(13) = {4, 1, 2, 3}; \n"], ...
["Plane Surface(14) = {13}; \n"], ...
["Line Loop(15) = {5, 6, 7, 8}; \n"], ...
["Plane Surface(16) = {15}; \n"], ...
["Line Loop(17) = {9, -8, -11, 4}; \n"], ...
["Ruled Surface(18) = {17}; \n"], ...
["Line Loop(19) = {12, -5, -9, 1}; \n"], ...
["Ruled Surface(20) = {19}; \n"], ...
["Line Loop(21) = {12, 6, -10, -2}; \n"], ...
["Ruled Surface(22) = {21}; \n"], ...
["Line Loop(23) = {11, -7, -10, 3}; \n"], ...
["Ruled Surface(24) = {23}; \n"], ...
["Surface Loop(25) = {18, 20, 22, 16, 24, 14}; \n"], ...
["Volume(26) = {25}; \n"]];
fname = tmpnam ();
[fid, msg] = fopen (strcat (fname, ".geo"), "w");
if (fid < 0); error (msg); endif
fputs (fid, gmsh_input);
fclose (fid);
msh = bim3c_mesh_properties (msh3m_gmsh (fname, "clscale", ".25"));
x = msh.p (1, :).';
u = x;
bnd = bim3c_unknowns_on_faces (msh, [14, 16]);
int = setdiff (1:columns (msh.p), bnd);
Mosc = bim3a_osc_advection_diffusion (msh, 1, msh.p(1,:)'*0);
Mgal = bim3a_advection_diffusion (msh, 1, msh.p(1,:)'*0);
u(int) = Mosc(int, int) \ ( - Mosc(int, bnd) * u(bnd));
uosc = u;
u(int) = Mgal(int, int) \ ( - Mgal(int, bnd) * u(bnd));
ugal = u;
fname_out = tmpnam ();
printf ("saving results to %s \n", strcat (fname_out, ".vtu"));
fpl_vtk_raw_write_field (fname_out, msh, {uosc, "u_osc"; ugal, "u_galerkin"}, {});
unlink (fname);
Produces the following output
Generating mesh... Processing gmsh data... Creating PDE-tool like mesh... Check for hanging nodes... Deleting temporary files... saving results to /var/tmp/oct-LgwrG1.vtu
The following code
gmsh_input = [["Point(1) = {0, 0, 0, .1}; \n"], ...
["Point(2) = {1, 0, 0, .1}; \n"], ...
["Point(3) = {0, -.3, 0, .1}; \n"], ...
["Point(4) = {0, +.3, 0, .1}; \n"], ...
["Point(5) = {1, -.3, 0, .1}; \n"], ...
["Point(6) = {1, 0.3, 0, .1}; \n"], ...
["Point(7) = {0, 0, -.3, .1}; \n"], ...
["Point(8) = {0, 0, +.3, .1}; \n"], ...
["Point(9) = {1, 0, -.3, .1}; \n"], ...
["Point(10) = {1, 0, 0.3, .1}; \n"], ...
["Circle(1) = {4, 1, 7}; \n"], ...
["Circle(2) = {7, 1, 3}; \n"], ...
["Circle(3) = {3, 1, 8}; \n"], ...
["Circle(4) = {8, 1, 4}; \n"], ...
["Circle(5) = {6, 2, 9}; \n"], ...
["Circle(6) = {9, 2, 5}; \n"], ...
["Circle(7) = {5, 2, 10}; \n"], ...
["Circle(8) = {10, 2, 6}; \n"], ...
["Line(9) = {4, 6}; \n"], ...
["Line(10) = {3, 5}; \n"], ...
["Line(11) = {8, 10}; \n"], ...
["Line(12) = {7, 9}; \n"], ...
["Line Loop(13) = {4, 1, 2, 3}; \n"], ...
["Plane Surface(14) = {13}; \n"], ...
["Line Loop(15) = {5, 6, 7, 8}; \n"], ...
["Plane Surface(16) = {15}; \n"], ...
["Line Loop(17) = {9, -8, -11, 4}; \n"], ...
["Ruled Surface(18) = {17}; \n"], ...
["Line Loop(19) = {12, -5, -9, 1}; \n"], ...
["Ruled Surface(20) = {19}; \n"], ...
["Line Loop(21) = {12, 6, -10, -2}; \n"], ...
["Ruled Surface(22) = {21}; \n"], ...
["Line Loop(23) = {11, -7, -10, 3}; \n"], ...
["Ruled Surface(24) = {23}; \n"], ...
["Surface Loop(25) = {18, 20, 22, 16, 24, 14}; \n"], ...
["Volume(26) = {25}; \n"]];
fname = tmpnam ();
[fid, msg] = fopen (strcat (fname, ".geo"), "w");
if (fid < 0); error (msg); endif
fputs (fid, gmsh_input);
fclose (fid);
msh = bim3c_mesh_properties (msh3m_gmsh (fname, "clscale", ".25"));
x = msh.p (1, :).';
u = x;
bnd = bim3c_unknowns_on_faces (msh, [14, 16]);
int = setdiff (1:columns (msh.p), bnd);
Mosc = bim3a_osc_advection_diffusion (msh, 1, msh.p(1,:)'*0);
Mgal = bim3a_advection_diffusion (msh, 1, msh.p(1,:)'*0);
f = bim3a_rhs (msh, 10, 1);
u(int) = Mosc(int, int) \ (f(int) - Mosc(int, bnd) * u(bnd));
uosc = u;
u(int) = Mgal(int, int) \ (f(int) - Mgal(int, bnd) * u(bnd));
ugal = u;
fname_out = tmpnam ();
printf ("saving results to %s \n", strcat (fname_out, ".vtu"));
fpl_vtk_raw_write_field (fname_out, msh, {uosc, "u_osc"; ugal, "u_galerkin"}, {});
unlink (fname);
Produces the following output
Generating mesh... Processing gmsh data... Creating PDE-tool like mesh... Check for hanging nodes... Deleting temporary files... saving results to /var/tmp/oct-m7e2KT.vtu
Package: bim