function exact_3d_graph (i_p) % Copyright 2000, J.E. Akin. All rights reserved. % ------------------------------------------------------ % Matlab line plot of i_p-th component value, in z-direction, % on top of line mesh (e.g. truss) at z=0. Mesh is black, % component is red. If i_p = 0, show RMS value. % ------------------------------------------------------ % c_x = x coordinates of nod_per_el line element % c_y = y coordinates of nod_per_el line element % msh_typ_nodes = connectivity list for elements, nt x nod_per_el % loop = corners for nod_per_el line element % 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 fprintf ('Begin component value 3D graph: \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 exact_node_solution.tmp nrr = size (exact_node_solution, 1); max_pp = size (exact_node_solution, 2) ; % number of columns if ( nrr == 0 ) error ('Error missing file exact_node_solution.tmp') end % if error fprintf ('Read %g exact nodal value sets \n', nrr) fprintf (' with %g components each \n', max_pp) if ( i_p > max_pp ) 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 t_z (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 c_z (nod_per_el + 1) = 0 ; % Optional pre-allocation % set constants [loop] = get_El_Loop (nod_per_el) ; %b fprintf ('Begin element plots:') % 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 = exact_node_solution(:, i_p) ; else % i_p = 0, get root mean sq for k = 1:np z (k) = sqrt ( sum (exact_node_solution (k, 1:max_pp).^2)) ; end % for k end % if get RMS value % Cite max, min values [V_X, L_X] = max (z(1:np)) ; [V_N, L_N] = min (z(1:np)) ; fprintf ('Max component value is %g at node %g \n', V_X, L_X) fprintf ('Min component value is %g at node %g \n', V_N, L_N) z_t (1:nod_per_el+1) = V_N ; % Fake z element coordinate null(1:np) = V_N ; % Fake z element coordinate % node to value connect x2 (1:2) = V_N; y2 (1:2) = V_N; z2 (1:2) = V_N; % Initialize plots xmax = max (x) ; xmin = min (x) ; ymax = max (y) ; ymin = min (y) ; zmax = V_X ; zmin = V_N ; if ( zmax == zmin ) zmax = zmin + 1 ; if ( ymax == ymin ) ymax = ymin + 1 ; end % if end % if planar or axial clf % clear graphics axis ([xmin, xmax, ymin, ymax, zmin, zmax]) % set axes % Get view azimuth and elevation [az, el] = view ; hold on % hold image for plots xlabel ('X, nodes at z=0, try view([0,1,0])') % add label ylabel ('Y, nodes at z=0, try view([1,0,0])') % add label if ( i_p >= 1 ) title(['Exact Solution Component\_', int2str(i_p), ... ' at ', int2str(np),' Nodes (', int2str(nt),' Elements)']) %title(['Exact Solution Component\_', int2str(i_p), ... %' Evaluated at the Nodes']) zlabel (['Component ', int2str(i_p), ' (max = ', ... num2str(V_X), ', min = ', num2str(V_N), ')']) else % i_p = 0, get root mean sq title(['Exact Solution RMS Value', ... ' at ', int2str(np),' Nodes (', int2str(nt),' Elements)']) %title(['Exact Solution RMS Value Evaluated at the Nodes']) zlabel (['RMS Value (max = ', ... num2str(V_X), ', min = ', num2str(V_N), ')']) 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)); % 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 % Plot this element c_x = t_x (loop) ; % x for nod_per_el line element c_y = t_y (loop) ; % y for nod_per_el line element c_z = t_z (loop) ; plot3 (c_x, c_y, z_t, 'g-') % plot mesh lines plot3 (c_x, c_y, c_z, 'k-') % plot result lines end % for over all elements plot3(c_x, c_y, c_z, 'k-'), grid % add grid to last one % Mark node points plot3 (x, y, null,'k*') plot3 (x, y, z, 'r*') % 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]) ; % Label node numbers at value location for i = 1:np x2 (1:2) = x(i) ; y2 (1:2) = y(i) ; z2 (1) = V_N ; z2(2) = z(i) ; plot3 (x2, y2, z2, 'r--') % connector % p_text = sprintf (' %g', i); % offset # from pt % text (x(i), y(i), z(i), p_text) end % for nodes % -depsc -tiff % for an eps version %bprint ('-dpsc', ['exact_', int2str(i_p), '_3d_graph']) hold off %bv_text = ['Created exact_', int2str(i_p), '_3d_graph.ps'] ; %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) % end of exact_3d_graph