function ave_node_flux_shrink (i_p) % Copyright 2000, J.E. Akin. All rights reserved. % ------------------------------------------------------ % Matlab carpet plot of i_p-th component value, % of the SCP average fluxes at mesh node locations % If i_p = 0, show RMS value % ------------------------------------------------------ % c_x = x coordinates of nod_per_el line polygon % c_y = y coordinates of nod_per_el line polygon % msh_typ_nodes = connectivity list for elements, nt x nod_per_el % loop = corners for nod_per_el line polygon % nod_per_el = Nodes per element % np = Number of Points % nt = Number of elements % pre_e = Element items before connectivity list pre_e = 0 ; % pre_p = Nodal items before coordinates pre_p = 1; % msh_bc_xyz = Nodal coordinates (with preceeding data) % S_1 = Element shrink ratio, 1 = none S_1 = 0.85 ; % S_2 = Gap size, = 1 - S_1 S_2 = 1 - S_1 ; % t_x = x coordinates of nod_per_el corners % t_y = y coordinates of nod_per_el corners % fprintf ('Begin flux component value shrink plots: \n') if ( nargin == 0 ) i_p = 0 ; end % if no arguments % Read coordinate file and connectivity file % integer bc code, real xy pairs for np points (pre_p = 1) load msh_bc_xyz.tmp ; % Set control data: number of points np = size (msh_bc_xyz,1) ; % number of nodal points fprintf ('Read %g mesh coordinate pairs \n', np) ns = size (msh_bc_xyz,2) - pre_p ; % space dimension if ( ns < 2 ) error ('This is not a 2D mesh') end % if not 2D data % Set control data: number elements load msh_typ_nodes.tmp ; % nod_per_el nodes per element nt = size (msh_typ_nodes,1) ; % number of elements in mesh nod_per_el = size (msh_typ_nodes,2) - pre_e - 1 ; % nodes per elem fprintf ('Read %g elements connections \n', nt) load scp_node_ave_fluxes.tmp nr = size (scp_node_ave_fluxes, 1); if ( nr == 0 ) error ('Error missing file scp_node_ave_fluxes.tmp') end % if error max_p = size (scp_node_ave_fluxes, 2) ; % number of columns fprintf ('Read %g nodal values of fluxes \n', nr) fprintf (' with %g components each \n', max_p) if ( i_p > max_p ) error ('i_p > available data') end % if error x (np) = 0. ; % pre-allocate array x y (np) = 0. ; % pre-allocate array y z (np) = 0. ; % pre-allocate array z t_nodes (nod_per_el) = 0 ; % Optional pre-allocation t_x (nod_per_el) = 0 ; % Optional pre-allocation t_y (nod_per_el) = 0 ; % Optional pre-allocation c_x (nod_per_el + 1) = 0 ; % Optional pre-allocation c_y (nod_per_el + 1) = 0 ; % Optional pre-allocation loop (nod_per_el + 1) = 0 ; % Optional pre-allocation % set constants [loop] = get_El_Loop (nod_per_el) ; % msh_bc_xyz has: pre_p items then: x, y x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column y = msh_bc_xyz (1:np, (pre_p+2)) ; % extract y column if ( i_p >= 1 ) z = scp_node_ave_fluxes(:, i_p) ; else % i_p = 0, get root mean sq for k = 1:np z (k) = sqrt ( sum (scp_node_ave_fluxes (k, 1:max_p).^2)) ; end % for k end % if get RMS value % Cite max, min values [V_X, L_X] = max (z) ; [V_N, L_N] = min (z) ; fprintf ('Max value is %g at node %g \n', V_X, L_X) fprintf ('Min value is %g at node %g \n', V_N, L_N) % Initialize plots xmax = max (x) ; xmin = min (x) ; ymax = max (y) ; ymin = min (y) ; zmax = max (z) ; zmin = min (z) ; clf % clear graphics axis ([xmin, xmax, ymin, ymax, zmin, zmax]) % set axes hold on % hold image for plots xlabel ('X') % add label ylabel ('Y') % add label if ( i_p >= 1 ) zlabel (['Flux ', int2str(i_p), ' (max = ', ... num2str(V_X), ', min = ', num2str(V_N), ')']) else % i_p = 0, get root mean sq zlabel (['RMS Value (max = ', ... num2str(V_X), ', min = ', num2str(V_N), ')']) end % if get RMS value if ( i_p >= 1 ) title(['Average FEA SCP Flux Component\_', int2str(i_p), ... ' at the ', int2str(np), 'Nodes (', int2str(nt),' Elements)']) % title(['Average FEA SCP Flux Component\_', int2str(i_p),' at the Nodes']) %title(['Average FEA SCP Flux Component\_', int2str(i_p),' at the Nodes']) else title(['Average FEA SCP Flux RMS Value at the ', ... int2str(np), ' Nodes (', int2str(nt),' Elements)']) %title(['Average FEA SCP Flux RMS Value at the Nodes']) end % if % Loop over all elements for it = 1:nt ; % Extract corner connectivity t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e + 1)); % Extract corner coordinates t_x = x (t_nodes) ; % x at those nodes, only t_y = y (t_nodes) ; % y at those nodes, only t_z = z (t_nodes) ; % z at those nodes, only % Get the centroid x_bar = sum (t_x' )/nod_per_el ; y_bar = sum (t_y' )/nod_per_el ; z_bar = sum (t_z' )/nod_per_el ; % Plot this polygon c_x = t_x (loop) ; % x for nod_per_el line polygon c_y = t_y (loop) ; % y for nod_per_el line polygon c_z = t_z (loop) ; c_x = S_1 * c_x + S_2 * x_bar ; % shrink x c_y = S_1 * c_y + S_2 * y_bar ; % shrink y c_z = S_1 * c_z + S_2 * z_bar ; % shrink z plot3 (c_x, c_y, c_z) % plot nod_per_el lines end % for over all elements plot3(c_x, c_y, c_z), grid % add grid to last one % label max min points % plot3 (x(L_X), y(L_X), V_X, 'kx') % plot3 (x(L_N), y(L_N), V_N, 'ko') v_text = sprintf ('----min') ; text (x(L_N), y(L_N), V_N, v_text) ; v_text = sprintf ('----max') ; text (x(L_X), y(L_X), V_X, [v_text]) ; ; % end % if show labels % -depsc -tiff % for an eps version %bprint ('-dpsc', ['ave_node_flux_', int2str(i_p), '_shrink_plot']) hold off %bv_text = ['Created ave_node_flux_', int2str(i_p), '_shrink_plot.ps']; %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) % end of ave_node_flux_shrink