rbmatlab  1.16.09
 All Classes Namespaces Files Functions Variables Modules Pages
stokes_plot_sim_data.m
Go to the documentation of this file.
1 function p = stokes_plot_sim_data(model, model_data, sim_data, plot_params)
2 %function p = stokes_plot_sim_data(model, model_data, sim_data, plot_params)
3 % Visualizes solution of stokes problem
4 %
5 % plot_params.modes is a cell array of strings. Possible mode-strings:
6 % - 'pressure'
7 % - 'velocity'
8 % - 'velocity_xcomp'
9 % - 'velocity_ycomp'
10 % - 'velocity_abs'
11 % - 'velocity_vec'
12 % - 'velocity_vec_plus_pressure'
13 %
14 % if plot_params.axes_handle is given, it is used and returned, else a cell
15 % array holding the figure handles is returned.
16 
17 % IM 24.10.2012
18 
19 if nargin < 4
20  plot_params = [];
21 end
22 
23 if ~isfield(plot_params, 'axis_equal')
24  plot_params.axis_equal = 1;
25 end
26 
27 if ~isfield(plot_params, 'title_on')
28  plot_params.title_on = 1;
29 end
30 
31 if ~isfield(plot_params, 'modes')
32  plot_params.modes = {'pressure', 'velocity'};
33 end
34 
35 if model.has_geometry_transformation
36  plot_params.global_transformation = @(glob) model.geometry_transformation(glob, model);
37 end
38 
39 xrange = model.xrange;
40 yrange = model.yrange;
41 
42 if model.has_geometry_transformation
43  corners = [xrange(1) yrange(1); xrange(1) yrange(2); xrange(2) yrange(1); xrange(2) yrange(2)];
44  corners = plot_params.global_transformation(corners);
45 
46  xrange(1) = min(corners([1 2], 1));
47  xrange(2) = max(corners([3 4], 1));
48  yrange(1) = min(corners([1 2], 2));
49  yrange(2) = max(corners([3 4], 2));
50 end
51 
52 if any(strncmp('velocity_vec', plot_params.modes, 12) + ...
53  strcmp('velocity', plot_params.modes))
54 
55  if ~isfield(plot_params, 'vector_distance')
56  plot_params.vector_distance = ...
57  max([xrange(2)-xrange(1), yrange(2)-yrange(1)]) / 30;
58  end
59 
60  if ~isfield(plot_params, 'vector_color')
61  plot_params.vector_color = 0;
62  end
63 end
64 
65 axes_handle_given = isfield(plot_params, 'axes_handle');
66 
67 p = {};
68 v = [xrange(1) xrange(2) yrange(1) yrange(2)];
69 
70 if axes_handle_given
71  axes(plot_params.axes_handle);
72  p = plot_params.axes_handle;
73  v = axis;
74  v(1) = min(v(1), xrange(1));
75  v(2) = max(v(2), xrange(2));
76  v(3) = min(v(3), yrange(1));
77  v(4) = max(v(4), yrange(2));
78 end
79 
80 
81 vel_ind = get_index(model_data.df_info, 'velocity');
82 pre_ind = get_index(model_data.df_info, 'pressure');
83 
84 vh = get_component(sim_data.uh, vel_ind);
85 ph = get_component(sim_data.uh, pre_ind);
86 
87 for i = 1:length(plot_params.modes)
88 
89  if isfield(plot_params, 'rbgui') && plot_params.rbgui
90  subplot(length(plot_params.modes)+1,1,i)
91  else
92  if ~axes_handle_given
93 
94  cf = figure;
95  p = [p(:); cf];
96  else
97 
98  hold on
99  end
100  end
101  switch (plot_params.modes{i})
102 
103  case ('pressure')
104 
105  plot(ph, plot_params);
106  if plot_params.title_on
107  title('pressure solution');
108  end
109 
110  case ('velocity_abs')
111 
112  % df_info's of components are assumed to be identical
113  vh_x = get_component(vh, 1);
114  vh_y = get_component(vh, 2);
115 
116  vh_abs = Fem.DiscFunc([], vh_x.df_info);
117  vh_abs.dofs = sqrt(vh_x.dofs.^2 + vh_y.dofs.^2);
118 
119  plot(vh_abs, plot_params);
120  if plot_params.title_on
121  title('velocity solution (absolute values)');
122  end
123 
124  case ('velocity_xcomp')
125 
126  vh_x = get_component(vh, 1);
127 
128  plot(vh_x, plot_params);
129  if plot_params.title_on
130  title('velocity solution (first component)');
131  end
132 
133  case ('velocity_ycomp')
134 
135  vh_y = get_component(vh, 2);
136 
137  plot(vh_y, plot_params);
138  if plot_params.title_on
139  title('velocity solution (second component)');
140  end
141 
142  case {'velocity_vec', 'velocity', 'velocity_vec_plus_pressure'}
143 
144  vh_vals = evaluate(vh, 1:vh.grid.nelements, [1/3, 1/3]);
145 
146  dist = plot_params.vector_distance;
147 
148  X = xrange(1)+dist/2 : dist : xrange(2);
149  Y = yrange(1)+dist/2 : dist : yrange(2);
150 
151  [X, Y] = meshgrid(X,Y);
152 
153  CXY = [model_data.grid.CX, model_data.grid.CY];
154  ESX = model_data.grid.ESX;
155  ESY = model_data.grid.ESY;
156 
157  if isfield(plot_params, 'global_transformation')
158  CXY = plot_params.global_transformation(CXY);
159  ESXY = plot_params.global_transformation([ESX(:), ESY(:)]);
160 
161  ESX = reshape(ESXY(:, 1), size(model_data.grid.ESX));
162  ESY = reshape(ESXY(:, 2), size(model_data.grid.ESY));
163  end
164 
165  xind = round((CXY(:, 1) - xrange(1)) / dist) + 1;
166  yind = round((CXY(:, 2) - yrange(1)) / dist) + 1;
167 
168  xind = xind - (xind > size(X, 2));
169  yind = yind - (yind > size(X, 1));
170 
171  ind = sub2ind(size(X), yind, xind);
172  X = X(:);
173  Y = Y(:);
174 
175  delete_ind = [];
176  for j = 1:length(ind)
177  [min_val, el_ind] = min((X(ind(j)) - CXY(:, 1)).^2 ...
178  + (Y(ind(j)) - CXY(:, 2)).^2);
179  min_val = min_val(1);
180  el_ind = el_ind(1);
181  edge_ind = find(model_data.grid.NBI(el_ind, :) < 0);
182  for k = 1:length(edge_ind)
183  X_virt = 2*ESX(el_ind, edge_ind(k)) ...
184  - CXY(el_ind, 1);
185  Y_virt = 2*ESY(el_ind, edge_ind(k)) ...
186  - CXY(el_ind, 2);
187  if (X(ind(j)) - X_virt)^2 + (Y(ind(j)) - Y_virt)^2 < min_val - 10*eps
188  delete_ind = [delete_ind; j];
189  break
190  end
191  end
192  end
193 
194  U = zeros(size(X));
195  V = zeros(size(X));
196 
197  U(ind) = vh_vals(:,1);
198  V(ind) = vh_vals(:,2);
199  U(ind(delete_ind)) = 0;
200  V(ind(delete_ind)) = 0;
201 
202  plot_title = 'velocity solution (vectors)';
203 
204  if strcmp(plot_params.modes{i}, 'velocity')
205 
206  % df_info's of components are assumed to be identical
207  vh_x = get_component(vh, 1);
208  vh_y = get_component(vh, 2);
209 
210  vh_abs = Fem.DiscFunc([], vh_x.df_info);
211  vh_abs.dofs = sqrt(vh_x.dofs.^2 + vh_y.dofs.^2);
212 
213  plot(vh_abs, plot_params);
214  plot_title = 'velocity solution';
215  v = axis;
216 
217  elseif strcmp(plot_params.modes{i}, 'velocity_vec_plus_pressure')
218 
219  plot(ph, plot_params);
220  plot_title = 'pressure solution and velocity vectors';
221  v = axis;
222  end
223 
224  if isfield(plot_params, 'vector_scale')
225  U = plot_params.vector_scale*U;
226  V = plot_params.vector_scale*V;
227  quiver(X ,Y, U, V, 'color', [0 0 plot_params.vector_color], ...
228  'AutoScale', 'off');
229  else
230  quiver(X, Y, U, V, 1, 'color', [0 0 plot_params.vector_color]);
231  end
232 
233  if plot_params.title_on
234  title(plot_title);
235  end
236  if plot_params.axis_equal
237  axis equal
238  axis(v);
239  end
240 
241  otherwise
242 
243  disp('plot mode unknown');
244  end
245 end
246 
247 end
248 
function p = stokes_plot_sim_data(model, model_data, sim_data, plot_params)
Visualizes solution of stokes problem.
dofs
DOF vector.
Definition: DiscFunc.m:70
represents a continous piecewise polynomial function of arbitrary dimension. Can be used for all fini...
Definition: DiscFunc.m:18