function both_error_norms_graph () % Copyright 2000, J.E. Akin. All rights reserved. % ------------------------------------------------------ % Matlab graph of error estimate value at mesh nodes % ------------------------------------------------------ % c_x = x 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 fprintf ('Begin error estimate value graph: \n') % 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 ~= 1) error ('This is not a 1D mesh') end % if not 1D 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 pt_ave_ex_error.tmp nr = size (pt_ave_ex_error, 1); if ( nr == 0 ) error ('Error missing file pt_ave_ex_error.tmp') end % if error max_p = size (pt_ave_ex_error, 2) ; % number of columns fprintf ('Read %g nodal error estimate values \n', nr) load pt_ave_error_est.tmp x (np) = 0. ; % pre-allocate array x y (np) = 0. ; % pre-allocate array y scp_y (np) = 0. ; % pre-allocate array y 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 scpt_y (nod_per_el) = 0 ; % Optional pre-allocation c_x (nod_per_el) = 0 ; % Optional pre-allocation c_y (nod_per_el) = 0 ; % Optional pre-allocation % set constants loop = [1:nod_per_el] ; % default to sequential order % msh_bc_xyz has: pre_p items then: x, y x = msh_bc_xyz (1:np, (pre_p+1)) ; % extract x column y = pt_ave_ex_error ; scp_y = pt_ave_error_est % Cite max, min values [V_X, L_X] = max (scp_y) ; [V_N, L_N] = min (scp_y) ; 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) null (1:np) = V_N ; % Initialize plots xmax = max (x) ; xmin = min (x) ; ymax = max (y) ; ymin = min (y) ; % ymax=50 if ( ymax == ymin ) if ( abs (ymax) > 0 ) ymax = ymax + abs (ymax)/20. ; ymin = ymin - abs (ymin)/20. ; end % if end % if clf % clear graphics ymin=0 ymax=15 axis ([xmin, xmax, ymin, ymax]) % set axes hold on % hold image for plots xlabel ('X, Node number at 45 deg, Element number at 90 deg') title(['Exact (dash), SCP (solid) Nodal Energy Norm Error: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) % ylabel (['SCP Error (max = ', ... % num2str(V_X), ', min = ', num2str(V_N), ')']) ylabel ('Energy Norm Error, % * 100') % Loop over all elements for it = 1:nt ; % Extract element connectivity t_nodes = msh_typ_nodes (it, (pre_e+2):(nod_per_el+pre_e+1)); % Skip point elements, if any if ( all (t_nodes) ) % then valid line % Extract element coordinates & values t_x = x (t_nodes) ; % x at those nodes, only t_y = y (t_nodes) ; % y at those nodes, only scpt_y = scp_y (t_nodes) ; % y at those nodes, only % Plot the element number x_bar = sum (t_x' )/nod_per_el ; t_text = sprintf (' (%g)', it); % offset # from pt text (x_bar, V_N, t_text, 'Rotation', 90) % incline plot (x_bar, V_N, 'k+') %b text (x_bar, 0.0, t_text, 'Rotation', 90) % incline %b plot (x_bar, 0.0, 'k+') % Plot this element (if non-sequential use loop) % c_x = t_x (loop) ; % x for nod_per_el line element % c_y = t_y (loop) ; % values at nodes % plot (c_x, c_y) % plot nod_per_el lines plot (t_x, t_y, 'r--') % plot nod_per_el lines plot (t_x, scpt_y, 'b-') % plot nod_per_el lines end % if zero in connectivity end % for over all elements % plot node points on axis for i = 1:np t_text = sprintf (' %g', i); % offset # from pt text (x(i), null(i), t_text, 'Rotation', 45) % incline end % for all plot (x, null, 'k*') grid % label max min points % plot (x(L_X), y(L_X), 'kx') % plot (x(L_N), y(L_N), 'ko') v_text = sprintf ('---min') ; text (x(L_N), V_N, v_text) ; v_text = sprintf ('---max') ; text (x(L_X), V_X, [v_text]) ; % -depsc -tiff % for an eps version %bprint ('-dpsc', ['both_error_norms_graph']) hold off %bv_text = ['Created both_error_norms_graph.ps'] ; %bfprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) % end of both_error_norms_graph