compute_gdl_velocity2.m

Go to the documentation of this file.

% script computing a velocity field for simple gdl model by solving a
% laplace problem with suitable boundary conditions
%
% three components are put together:
%
%    .p0...______......______......______..p1..
%   |                                         |
%   |                                         |
%    -----------------------------------------
% 
%   |- = neumann 0 conditions
%   ... dirichlet conditions: linearly decreasing pressure 

% Bernard Haasdonk 6.7.2006

%%%%%%%%%%%%%%%%%%%%%%%% compute pressure
params = [];

disp('setting parameters');

params.xrange = [-1250e-6,2250e-6]; % grid x-coordinate range
params.yrange = [0,200e-6];  % grid y-coordinate range
params.xnumintervals = 497;  % number of grid cells in x-direction: 
                             % multiple of 7
params.ynumintervals = 41;   % number of grid cells in y-direction

% set upper boundary to dirichlet, middle of upper and the lower boundary
% to neuman, set left and right to neuman
params.bnd_rect_corner1 = [params.xrange(1) , params.yrange(2) ; ...
                    params.xrange(1)+ 500e-6 , params.yrange(2); ...
                    params.xrange(1)+ 1500e-6 , params.yrange(2); ...
                    params.xrange(1)+ 2500e-6 , params.yrange(2); ...
                    params.xrange(1) , params.yrange(1); ...
                    params.xrange(1) , params.yrange(1); ...
                    params.xrange(2) , params.yrange(1);
                   ];
params.bnd_rect_corner1 = params.bnd_rect_corner1' - eps;
params.bnd_rect_corner2 = [params.xrange(2), params.yrange(2); ...
                    params.xrange(1)+1000e-6, params.yrange(2); ...
                    params.xrange(1)+2000e-6, params.yrange(2); ...
                    params.xrange(1)+3000e-6, params.yrange(2); ...
                    params.xrange(2), params.yrange(1); ...
                    params.xrange(1), params.yrange(2); ...
                    params.xrange(2), params.yrange(2)
                   ];
params.bnd_rect_corner2 = params.bnd_rect_corner2' + eps;

params.bnd_rect_index = [-1, -2, -2, -2, -2, -2, -2];
params.name_dirichlet_values = 'xstripes';
params.c_dir = [4,3,2,1];
% separating borders for changing value
params.dir_borders = [params.xrange(1)+750e-6, ... 
                    params.xrange(1)+1750e-6, ...
                    params.xrange(1)+2750e-6]; 
%params.name_dirichlet_values = 'leftright';
%params.c_dir_left = 1.0; % pressure left
%params.c_dir_right = 0.5; % pressure right
%params.dir_middle = params.xrange(2)/2;

%params.name_dirichlet_values = 'homogeneous';
%params.c_dir = 1.0;
params.name_neuman_values = 'zero';
%params.name_neuman_values = 'homogeneous';
%params.c_neu = 0;

params.show_colorbar = 1;
params.no_lines = 1;
params.verbose = 10;

grid = rectgrid(params);

%tfn = fullfile(getenv('RBMATLABTEMP'),'gdl_pressure2.mat');
tfn = fullfile(rbmatlabhome,'datafunc','data','gdl_pressure2.mat');

if ~exist(tfn,'file')  
  p = fem_laplace(grid,params);

  figure;
  plot_vertex_data(grid,p,params);
  
  %  [px,py] = gradient(reshape(p,params.nx+1,params.ny+1)',...
  %                  params.xrange(2)/params.nx,params.yrange(2)/params.ny);
  %  figure;
  %  quiver(1:(params.nx+1),1:(params.ny+1),-4*px,-4*py);
  
  save(tfn,'p','params');
else
  disp('skipping computation of data, as gdl_pressure2.mat exists')
end;


% plot velocity:

% evaluate velocity in a second grid and plot as piecewise constant
load(tfn,'p','params');
params2 = [];
params2.xnumintervals = 350;
params2.ynumintervals = 50;
params2.xrange = params.xrange;
params2.yrange = params.yrange;
params2.verbose = 0;
params2.show_colorbar = 1;
params2.no_lines = 1;
params2.axis_tight = 1;

grid2 = rectgrid(params2);
[ux, uy] = evaluate_gradient(grid2.CX, grid2.CY, p, params);

factor = 1/1000;
ux = -ux/factor;
uy = -uy/factor;
%plot_element_data_sequence(grid2,[-ux,-uy],params2);
params2.plot = @fv_plot;
plot_sequence([-ux,-uy],grid2,params2);
title('gradient on whole domain');

params2.xrange = [0e-6,1000e-6]; % grid x-coordinate range
params2.yrange = params.yrange;
params2.xnumintervals = 200;
grid2 = rectgrid(params2);
[ux, uy] = evaluate_gradient(grid2.CX, grid2.CY,p,params);

ux = -ux/factor;
uy = -uy/factor;
plot_sequence([-ux,-uy],grid2,params2);
%plot_element_data_sequence(grid2,[-ux,-uy],params2);
title('gradient on part of domain');