3 % Visualizes solution of stokes problem
5 % plot_params.modes is a cell array of strings. Possible mode-strings:
12 % -
'velocity_vec_plus_pressure'
14 %
if plot_params.axes_handle is given, it is used and returned,
else a cell
15 % array holding the figure handles is returned.
23 if ~isfield(plot_params,
'axis_equal')
24 plot_params.axis_equal = 1;
27 if ~isfield(plot_params, 'title_on')
28 plot_params.title_on = 1;
31 if ~isfield(plot_params, 'modes')
32 plot_params.modes = {
'pressure',
'velocity'};
35 if model.has_geometry_transformation
36 plot_params.global_transformation = @(glob) model.geometry_transformation(glob, model);
39 xrange = model.xrange;
40 yrange = model.yrange;
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);
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));
52 if any(strncmp('velocity_vec', plot_params.modes, 12) + ...
53 strcmp('velocity', plot_params.modes))
55 if ~isfield(plot_params, 'vector_distance')
56 plot_params.vector_distance = ...
57 max([xrange(2)-xrange(1), yrange(2)-yrange(1)]) / 30;
60 if ~isfield(plot_params, 'vector_color')
61 plot_params.vector_color = 0;
65 axes_handle_given = isfield(plot_params, 'axes_handle');
68 v = [xrange(1) xrange(2) yrange(1) yrange(2)];
71 axes(plot_params.axes_handle);
72 p = plot_params.axes_handle;
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));
81 vel_ind = get_index(model_data.df_info, 'velocity');
82 pre_ind = get_index(model_data.df_info, 'pressure');
84 vh = get_component(sim_data.uh, vel_ind);
85 ph = get_component(sim_data.uh, pre_ind);
87 for i = 1:length(plot_params.modes)
89 if isfield(plot_params, 'rbgui') && plot_params.rbgui
90 subplot(length(plot_params.modes)+1,1,i)
101 switch (plot_params.modes{i})
105 plot(ph, plot_params);
106 if plot_params.title_on
107 title(
'pressure solution');
110 case (
'velocity_abs')
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);
116 vh_abs = Fem.DiscFunc([], vh_x.df_info);
117 vh_abs.dofs = sqrt(vh_x.dofs.^2 + vh_y.dofs.^2);
119 plot(vh_abs, plot_params);
120 if plot_params.title_on
121 title('velocity solution (absolute values)
');
124 case ('velocity_xcomp
')
126 vh_x = get_component(vh, 1);
128 plot(vh_x, plot_params);
129 if plot_params.title_on
130 title('velocity solution (first component)
');
133 case ('velocity_ycomp
')
135 vh_y = get_component(vh, 2);
137 plot(vh_y, plot_params);
138 if plot_params.title_on
139 title('velocity solution (second component)
');
142 case {'velocity_vec
', 'velocity
', 'velocity_vec_plus_pressure
'}
144 vh_vals = evaluate(vh, 1:vh.grid.nelements, [1/3, 1/3]);
146 dist = plot_params.vector_distance;
148 X = xrange(1)+dist/2 : dist : xrange(2);
149 Y = yrange(1)+dist/2 : dist : yrange(2);
151 [X, Y] = meshgrid(X,Y);
153 CXY = [model_data.grid.CX, model_data.grid.CY];
154 ESX = model_data.grid.ESX;
155 ESY = model_data.grid.ESY;
157 if isfield(plot_params, 'global_transformation
')
158 CXY = plot_params.global_transformation(CXY);
159 ESXY = plot_params.global_transformation([ESX(:), ESY(:)]);
161 ESX = reshape(ESXY(:, 1), size(model_data.grid.ESX));
162 ESY = reshape(ESXY(:, 2), size(model_data.grid.ESY));
165 xind = round((CXY(:, 1) - xrange(1)) / dist) + 1;
166 yind = round((CXY(:, 2) - yrange(1)) / dist) + 1;
168 xind = xind - (xind > size(X, 2));
169 yind = yind - (yind > size(X, 1));
171 ind = sub2ind(size(X), yind, xind);
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);
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)) ...
185 Y_virt = 2*ESY(el_ind, edge_ind(k)) ...
187 if (X(ind(j)) - X_virt)^2 + (Y(ind(j)) - Y_virt)^2 < min_val - 10*eps
188 delete_ind = [delete_ind; j];
197 U(ind) = vh_vals(:,1);
198 V(ind) = vh_vals(:,2);
199 U(ind(delete_ind)) = 0;
200 V(ind(delete_ind)) = 0;
202 plot_title = 'velocity solution (vectors)
';
204 if strcmp(plot_params.modes{i}, 'velocity
')
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);
211 vh_abs.
dofs = sqrt(vh_x.dofs.^2 + vh_y.dofs.^2);
213 plot(vh_abs, plot_params);
214 plot_title =
'velocity solution';
217 elseif strcmp(plot_params.modes{i},
'velocity_vec_plus_pressure')
219 plot(ph, plot_params);
220 plot_title = 'pressure solution and velocity vectors';
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], ...
230 quiver(X, Y, U, V, 1, 'color', [0 0 plot_params.vector_color]);
233 if plot_params.title_on
236 if plot_params.axis_equal
243 disp('plot mode unknown');
function p = stokes_plot_sim_data(model, model_data, sim_data, plot_params)
Visualizes solution of stokes problem.
represents a continous piecewise polynomial function of arbitrary dimension. Can be used for all fini...