function ave_est_exact_norms () % 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 % 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 x (np) = 0. ; % pre-allocate array x y (np) = 0. ; % pre-allocate array y yE(np) = 0. ; % pre-allocate array y % 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 exact error norms \n', nr) y = pt_ave_ex_error ; % Cite max, min values [V_X, L_X] = max (y) ; [V_N, L_N] = min (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 ; load pt_ave_error_est.tmp nr = size (pt_ave_error_est, 1); if ( nr == 0 ) error ('Error missing file pt_ave_error_est.tmp') end % if error max_p = size (pt_ave_error_est, 2) ; % number of columns fprintf ('Read %g nodal estimate error norms \n', nr) yE = pt_ave_error_est ; % Cite max, min values [VE_X, LE_X] = max (yE) ; [VE_N, LE_N] = min (yE) ; fprintf ('Max value is %g at node %g \n', VE_X, LE_X) fprintf ('Min value is %g at node %g \n', VE_N, LE_N) 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) = 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 % Initialize plots xmax = max (x) ; xmin = min (x) ; ymax = max (y) ; ymin = min (y) ; Emax = max (yE) ; Emin = min (yE) ; ymax = max (Emax, ymax) ; ymin = min (Emin, ymin) ; % 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 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(['FEA Estimated & Exact Nodal Energy Norm Error, % * 100: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) ylabel (['Error Norm (o - Exact max = ', ... num2str(V_X), ', min = ', num2str(V_N), ')']) % 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 % 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+') % Plot this element (if non-sequential use loop) plot (t_x, t_y) % plot nod_per_el lines t_y = yE (t_nodes) ; plot (t_x, t_y, 'k--') plot (t_x, t_y, 'ko') 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 %b print ('-dpsc', ['ave_est_exact_norms']) hold off %b v_text = ['Created ave_est_exact_norms.ps'] ; %b fprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) % end of ave_est_exact_norms