function color_result_brick (i_p) % Copyright 2000, J.E. Akin. All rights reserved. % ------------------------------------------------------ % Matlab carpet plot of i_p-th component value, % at mesh node locations % If i_p = 0, show RMS value % ------------------------------------------------------ %? global az el ? these are not changing with new view % 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) % t_x = x coordinates of nod_per_el corners % t_y = y coordinates of nod_per_el corners % Drawing connection order for 6 faces % f1 1 2 3 4 & 1 % f2 1 5 6 2 & 1 % f3 2 6 7 3 & 2 % f4 3 7 8 4 & 3 % f5 1 4 8 5 & 1 % f6 5 8 7 6 & 5 fprintf ('Begin component value carpet 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 ~= 3 ) error ('This is not a 3D mesh') end % if not 3D 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) if ( nod_per_el ~= 8 ) error ('This is not a hexahedra mesh') end % if H8 face_nodes = 4; faces=6; load node_results.tmp nr = size (node_results, 1); if ( nr == 0 ) error ('Error missing file node_results.tmp') end % if error max_p = size (node_results, 2) ; % number of columns fprintf ('Read %g nodal solution values \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 R (np) = 0. ; % pre-allocate array R 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 t_z (nod_per_el) = 0 ; % Optional pre-allocation t_R (nod_per_el) = 0 ; % Optional pre-allocation c_x (face_nodes + 1) = 0 ; % Optional pre-allocation c_y (face_nodes + 1) = 0 ; % Optional pre-allocation c_z (face_nodes + 1) = 0 ; % Optional pre-allocation c_R (face_nodes + 1) = 0 ; % Optional pre-allocation loop (face_nodes + 1) = 0 ; % Optional pre-allocation % set constants % 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 z = msh_bc_xyz (1:np, (pre_p+3)) ; % extract z column if ( i_p >= 1 ) R = node_results(:, i_p) ; else % i_p = 0, get root mean sq for k = 1:np R (k) = sqrt ( sum (node_results (k, 1:max_p).^2)) ; end % for k end % if get RMS value % Cite max, min values [V_X, L_X] = max (R) ; [V_N, L_N] = min (R) ; 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) caxis([V_N V_X]) ; % set full color ranges % 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 zlabel ('Z') if ( i_p >= 1 ) title(['FEA Hexahedra Component\_', int2str(i_p), ... ' (', num2str(V_N), ' to ' num2str(V_X),'): ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) else % i_p = 0, get root mean sq title(['FEA Hexahedra RMS Value: ', ... ' (', num2str(V_N), ' to ' num2str(V_X),'): ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) end % if get RMS value % 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)); if ( all (t_nodes > 0) ) % skips face and edge elements % 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 t_R = R (t_nodes) ; % R at those nodes, only % Loop over each face for jf = 1:faces % Extract face nodes [loop] = get_H8_Face_Loop (jf) ; % continue if normal vector is visible % 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_R = t_R (loop) ; color= c_R ; fill3 (c_x, c_y, c_z, color) % plot nod_per_el lines end % for faces end % if all end % for over all elements fill3(c_x, c_y, c_z, color), grid % add grid to last one colorbar % label max min points 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', ['color_result_', int2str(i_p), '_brick']) hold off %bv_text = ['Created color_result_', int2str(i_p), '_brick.ps'] ; %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) % end of color_result_brick