/usr/lib/x86_64-linux-gnu/scilab-getfem++/demos/demo_stokes_3D_tank_alt.sce is in scilab-getfem++ 4.2.1~beta1~svn4635~dfsg-3+b1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 | // old example, which uses the deprecated gf_solve function
// see demo_stokes_3D_tank.m for a "modern" version
lines(0);
stacksize('max');
path = get_absolute_file_path('demo_stokes_3D_tank_alt.sce');
if getos()=='Windows' then
// Under Windows, all the trace messages are available in the dos console
// Under Linuxs, all the trace messages are redirected to the Scilab console
consolebox('on');
end
gf_util('trace level',3);
gf_util('warning level',3);
gf_workspace('clear all');
printf('demo stokes_3D_tank_alt started\n');
disp('3D stokes demonstration on a quadratic mesh -- 512MB of memory needed for the solve!!');
compute = input(' 1:compute the solution\n 0:load a previously computed solution\n ? ');
global verbosity; verbosity=1;
pde = init_pde();
pde('type') = 'stokes'; // YC:
pde('viscos') = 1.0;
pde('F') = list(list(0,0,0));
pde = add_empty_bound(pde);
pde('bound')($)('type') = 'Dirichlet';
pde('bound')($)('R') = list(list('9-(y.^2+(z-6.0).^2)',0,0));
pde = add_empty_bound(pde);
pde('bound')($)('type') = 'Dirichlet';
pde('bound')($)('R') = list(list('9-(y.^2+(z-6.0).^2)',0,0));
pde = add_empty_bound(pde);
//pde('bound')($)('type') = 'Mixed';
pde('bound')($)('type') = 'Dirichlet';
pde('bound')($)('R') = list(0,0,0);
pde('bound')($)('H') = list(list(0,0,0),list(0,0,0),list(0,0,1));
pde('bound')($)('G') = list(list(0,0,0));
pde = add_empty_bound(pde);
pde('bound')($)('type') = 'Dirichlet';
pde('bound')($)('R') = list(list(0,0,0));
m = gf_mesh('import','GiD', path + '/data/tank_quadratic_2500.GiD.msh');
mfulag = gf_mesh_fem(m,3);
pde('mf_u') = [];
pde('mf_p') = [];
pde('mf_d') = [];
pde('mim') = [];
pde('mf_u') = gf_mesh_fem(m,3);
pde('mf_p') = gf_mesh_fem(m,1);
pde('mf_d') = gf_mesh_fem(m,1);
pde('mim') = gf_mesh_im(m, gf_integ('IM_TETRAHEDRON(5)'));
// this is a good example of the usefulness of the cubic bubble
// -> if not used, the pression has strange values
//gf_mesh_fem_set(pde('mf_u'),'fem',gf_fem('FEM_PK_WITH_CUBIC_BUBBLE(3,2)')
gf_mesh_fem_set(pde('mf_u'),'fem',gf_fem('FEM_PK(3,2)'));
gf_mesh_fem_set(pde('mf_d'),'fem',gf_fem('FEM_PK(3,2)'));
//gf_mesh_fem_set(pde('mf_p'),'fem',gf_fem('FEM_PK_DISCONTINUOUS(3,1)'))
gf_mesh_fem_set(pde('mf_p'),'fem',gf_fem('FEM_PK(3,1)'));
// we use a P3 mesh fem for interpolation of the U field, since
// because of its cubic bubble function, the pde.mf_u is not lagrangian
gf_mesh_fem_set(mfulag,'fem',gf_fem('FEM_PK(3,1)'));
all_faces = gf_mesh_get(m, 'outer faces', gf_mesh_get(m, 'cvid'));
P = gf_mesh_get(m,'pts');
INpid = find(abs(P(1,:)+25) < 1e-4);
OUTpid = find(abs(P(1,:)-25) < 1e-4);
TOPpid = find(abs(P(3,:)-20) < 1e-4);
INfaces = gf_mesh_get(m, 'faces from pid', INpid);
OUTfaces = gf_mesh_get(m, 'faces from pid', OUTpid);
TOPfaces = gf_mesh_get(m, 'faces from pid', TOPpid);
gf_mesh_set(m, 'boundary', 1, INfaces);
gf_mesh_set(m, 'boundary', 2, OUTfaces);
gf_mesh_set(m, 'boundary', 3, TOPfaces);
gf_mesh_set(m, 'boundary', 4, _setdiff(all_faces',union(union(INfaces',OUTfaces','r'),TOPfaces','r'),'rows')');
disp(sprintf('nbdof: mf_u=%d, mf_p=%d',gf_mesh_fem_get(pde('mf_u'),'nbdof'),gf_mesh_fem_get(pde('mf_p'),'nbdof')));
if (compute) then
// unfortunately, the basic stokes solver is very slow...
// on this problem, the fastest way is to reduce to a (full) linear system on the pression...
// drawback: scilab will be killed if you don't have 512MB of memory
pde('solver') = 'brute_stokes';
tic;
[U,P] = gf_solve(pde);
disp(sprintf('solve done in %.2f sec', toc));
save(path + '/demo_stokes_3D_tank_UP.mat',U,P);
disp('[the solution has been saved in ''demo_stokes_3D_tank_UP.mat'']');
else
load(path + '/demo_stokes_3D_tank_UP.mat');
end
mfu = pde('mf_u');
mfp = pde('mf_p');
mfd = pde('mf_d');
disp('Got a solution, now you can call demo_stokes_3D_tank_draw to generate graphics');
printf('demo stokes_3D_tank_alt terminated\n');
|