iunction Axisymmetric_Stress_T6_Isoparametric (load_pt) % load_pt = 1 flags that point (ring) load will be input % via file msh_load_pt.tmp on a PER UNIT RADIAN basis %............................................................. % Axisymmetric Stress Analysis (PER UNIT RADIAN) %............................................................. if ( nargin == 0 ) ; % check for optional data load_pt = 0 ; % no point source data end % if from argument count echo_user_remarks % optional use of new file msh_remarks.tmp % Application and element dependent controls n_g = 2 ; % number of DOF per node n_r = 4 ; % number of rows in B_e matrix % Read mesh input data files [n_m, n_s, P, x, y, z] = get_mesh_nodes ; % Extract EBC flags from packed integer flag P [EBC_flag] = get_ebc_flags (n_g, n_m, P) ; % unpack flags EBC_count = sum( sum ( EBC_flag > 0 ) ) ; % # of EBC [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements ; n_d = n_g*n_m ; % system degrees of freedom (DOF) n_i = n_g*n_n ; % number of DOF per element S = zeros (n_d, n_d) ; C = zeros (n_d, 1) ; % initalize sums M = zeros (n_d, n_d) ; ; % initalize sums % Read EBC values, if any if ( EBC_count > 0 ) ; % need EBC data [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data end % if any EBC data expected % Read point (ring) source PER UNIT RADIAN, if any, & insert in C if ( load_pt > 0 ) ; % need point loads data [C] = get_and_add_point_sources (n_g, n_m, C); % add point loads end % if any point source expected C_react = C ; % save before BC for later reaction use % Load the element properties array load msh_properties.tmp ; % one row per element n_prop = size(msh_properties, 1) ; % == 1 if homogeneous fprintf ('\n(Echoing columns of file msh_properties.tmp) \n') if ( n_prop == 1 ) E = msh_properties (1, 1) ; % modulus of elasticity nu = msh_properties (1, 2) ; % Poisson's ratio Omega = msh_properties (1, 3) ; % y angular velocity for x-force Axial_a = msh_properties (1, 4) ; % axial acceleration for y-force Rho = msh_properties (1, 5) ; % mass density % ECHO PROPERTIES fprintf ('\nProperties for all elements \n') fprintf ('Elastic Modulus = %g \n', E) fprintf ('Poisson;s ratio = %g \n', nu) fprintf ('Y angular velocity = %g \n', Omega) fprintf ('Y axial acceleration = %g \n', Axial_a) fprintf ('Mass Density = %g \n', Rho) end % if if ( n_prop > n_e ) error ('ERROR: number of property sets exceeds number of elements') end % if % Tabulate quadrature data for unit triangle n_q = 7 ; % exact for 4-th degree polynomials r_q = [ 0.10128651, 0.47014206, 0.79742699, 0.47014206, ... 0.10128651, 0.05971587, 0.33333333 ] ; s_q = [ 0.10128651, 0.05971587, 0.10128651, 0.47014206, ... 0.79742699, 0.47014206, 0.33333333 ] ; w_q = [ 0.06296959, 0.06619708, 0.06296959, 0.06619708, ... 0.06296959, 0.06619708, 0.11250000 ] ; % GENERATE ELEMENT MATRICES AND ASSYMBLE INTO SYSTEM % Assemble n_d by n_d square matrix terms from n_e by n_e for j = 1:n_e ; % element loop ===>> ===>> ===>> ===>> ===>> ===>> e_nodes = nodes (j, 1:n_n) ; % connectivity % SET ELEMENT PROPERTIES & GEOMETRY if ( n_prop > 1 ) E = msh_properties (j, 1) ; % modulus of elasticity nu = msh_properties (j, 2) ; % Poisson's ratio Omega = msh_properties (j, 3) ; % y angular velocity for x-force Axial_a = msh_properties (j, 4) ; % axial acceleration for y-force Rho = msh_properties (j, 5) ; % mass density % ECHO PROPERTIES fprintf ('\nProperties for element %i \n', j) fprintf ('Elastic Modulus = %g \n', E) fprintf ('Poisson;s ratio = %g \n', nu) fprintf ('Y angular velocity = %g \n', Omega) fprintf ('Y axial acceleration = %g \n', Axial_a) fprintf ('Mass Density = %g \n', Rho) end % if E_e = E / ((1 + nu) * (1 - nu*2)) * ... % 4 x 4 [(1-nu), nu, 0, nu ; ... nu, (1-nu), 0, nu ; ... 0, 0, (1-2*nu)/2, 0 ; ... nu, nu, 0, (1-nu)] ; % stress-strain law % Gather nodal coordinates xy_e (1:n_n, 1) = x(e_nodes) ; % x coord at el nodes xy_e (1:n_n, 2) = y(e_nodes) ; % 6 x 2 % y coord at el nodes % Numerical Intergation of S_e, C_e S_e = zeros (n_i, n_i) ; C_e = zeros (n_i, 1) ; % initalize values B_q = zeros (n_r, n_i) ; % initalize values M_e = zeros (n_i, n_i) ; % initalize values for q = 1:n_q ; % begin loop ---> ---> ---> ---> ---> ---> ---> ---> r = r_q (q) ; s = s_q (q) ; w = w_q (q) ; % recover data % element type interpolation functions, 1 x 6 H_q = [ (1 - 3*r - 3*s + 2*r*r + 4*r*s + 2*s*s), (2*r*r - r), ... (2*s*s - s), 4*(r - r*r - r*s), 4*r*s, 4*(s - r*s - s*s) ] ; DLH_q = [ (-3 + 4 * r + 4 * s), (-1 + 4 * r), 0, ... (4 - 8 * r - 4 * s), (4 * s), (-4 * s) ; ... % dH/dr (-3 + 4 * r + 4 * s), 0, (-1 + 4 * s), ... (-4 * r), (4 * r), (4 - 4 * r - 8 * s) ] ; % dH/ds % Global position and Jacobian xy_q = H_q * xy_e ; % 1 x 2 % global (x, y) Jacobian = DLH_q * xy_e ; % 2 x 2 % d(x,y)/d(r,s) J_det = abs( det (Jacobian)) ; % 1 x 1 % determinant J_Inv = inv (Jacobian) ; % 2 x 2 % inverse % Global derivatives of H DGH_q = J_Inv * DLH_q ; % 2 x 6 % dH/dx & dH/dy E_q = E_e ; % 4 x 4 % properties % Insert gradient terms into the differential operator, 4 x 12 B_q (1, 1:2:11) = DGH_q (1, 1:6) ; % for esplion_x B_q (2, 2:2:12) = DGH_q (2, 1:6) ; % for esplion_y B_q (3, 1:2:11) = DGH_q (2, 1:6) ; % for gamma B_q (3, 2:2:12) = DGH_q (1, 1:6) ; % for gamma R = xy_q (1) ; % radius B_q (4, 1:2:11) = H_q (1:6) / R ; % for hoop esplion % ELEMENT MATRICES UPDATES, 12 x 12 and 12 x 1, dv = 2pi R dA S_e = S_e + (B_q' * E_q * B_q) * R * J_det * w_q(q) ; % square M_e (1:2:11) = M_e (1:2:11) ... + (H_q' * Rho * H_q) * R * J_det * w_q (q) ; % x mass M_e (2:2:12) = M_e (2:2:12) ... + (H_q' * Rho * H_q) * R * J_det * w_q (q) ; % y mass if ( any ([Omega, Axial_a]) ) % then need acceleration loads Body_x = Rho * R * Omega^2 ; % centrifugal load Body_y = -Rho * Axial_a ; % axial load C_e(1:2:11) = C_e(1:2:11) + H_q' * Body_x * R * J_det * w_q (q); % Fx C_e(2:2:12) = C_e(2:2:12) + H_q' * Body_y * R * J_det * w_q (q); % Fy end % if body forces active end % for quadratures <--- <--- <--- <--- <--- <--- <--- <--- <--- % SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS % Insert completed element matrices into system matrices [rows] = get_element_index (n_g, n_n, e_nodes); % eq numbers S (rows, rows) = S (rows, rows) + S_e ; % add to system sq C (rows) = C (rows) + C_e ; % add to sys column end % for elements <<==== <<==== <<==== <<==== <<==== <<==== <<==== % ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY if ( EBC_count > 0 ) ; % reactions occur [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C); end % if essential BC exist (almost always true) % ENFORCE ESSENTIAL BOUNDARY CONDITIONS save_resultant_load_vectors (n_g, C) [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C); % COMPUTE SOLUTION & SAVE T = S \ C ; % Compute temperatures list_save_displacements_results (n_g, n_m, T) % OPTIONAL REACTION RECOVERY & SAVE if ( EBC_count > 0 ) ; % reactions exist ? [EBC_react] = recover_reactions_print_save (n_g, n_d, ... EBC_flag, EBC_row, EBC_col, T); % reaction to EBC end % if EBC exist % POST-PROCESS ELEMENTS (element nodal heat flows) fprintf ('\nPost-prossessing Gradients and Reaction: \n') fid = fopen('el_qp_xyz_fluxes.tmp', 'w') ; % open for writing for j = 1:n_e ; % loop over all elements ====>> ====>> ====>> ====>> e_nodes = nodes (j, 1:n_n) ; % connectivity [rows] = get_element_index (n_g, n_n, e_nodes) ; % eq numbers T_e (1:n_i) = T(rows) ; % gather element displacements % SET ELEMENT PROPERTIES & GEOMETRY if ( n_prop > 1 ) E = msh_properties (j, 1) ; % modulus of elasticity nu = msh_properties (j, 2) ; % Poisson's ratio Omega = msh_properties (j, 3) ; % y angular velocity for x-force Axial_a = msh_properties (j, 4) ; % axial acceleration for y-force Rho = msh_properties (j, 5) ; % mass density end % if E_e = E / ((1 + nu) * (1 - nu*2)) * ... % 4 x 4 [(1-nu), nu, 0, nu ; ... nu, (1-nu), 0, nu ; ... 0, 0, (1-2*nu)/2, 0 ; ... nu, nu, 0, (1-nu)] ; % stress-strain law % Gather nodal coordinates xy_e (1:n_n, 1) = x(e_nodes) ; % x coord at el nodes xy_e (1:n_n, 2) = y(e_nodes) ; % 6 x 2 % y coord at el nodes % Numerical Intergation of S_e, C_e S_e = zeros (n_i, n_i) ; C_e = zeros (n_i, 1) ; % initalize values B_q = zeros (n_r, n_i) ; % initalize values for q = 1:n_q ; % begin loop ---> ---> ---> ---> ---> ---> ---> r = r_q (q) ; s = s_q (q) ; w = w_q (q) ; % recover data % element type interpolation functions H_q = [ (1 - 3*r - 3*s + 2*r*r + 4*r*s + 2*s*s), (2*r*r - r), ... (2*s*s - s), 4*(r - r*r - r*s), 4*r*s, 4*(s - r*s - s*s) ] ; DLH_q = [ (-3 + 4 * r + 4 * s), (-1 + 4 * r), 0, ... (4 - 8 * r - 4 * s), (4 * s), (-4 * s) ; ... % dH/dr (-3 + 4 * r + 4 * s), 0, (-1 + 4 * s), ... (-4 * r), (4 * r), (4 - 4 * r - 8 * s) ] ; % dH/ds % Global position and Jacobian xy_q = H_q * xy_e ; % 1 x 2 % global (x, y) Jacobian = DLH_q * xy_e ; % 2 x 2 % d(x,y)/d(r,s) J_det = abs( det (Jacobian)) ; % 1 x 1 % determinant J_Inv = inv (Jacobian) ; % 2 x 2 % inverse % Global derivatives of H DGH_q = J_Inv * DLH_q ; % 2 x 6 % dH/dx & dH/dy E_q = E_e ; % 3 x 3 % properties % Output point location, and sress components fprintf ('\nElement, Point, Coordinates %i %i %g %g \n', ... j, q, xy_q (1:2)) % Insert gradient terms into the differential operator B_q (1, 1:2:11) = DGH_q (1, 1:6) ; % for esplion_x B_q (2, 2:2:12) = DGH_q (2, 1:6) ; % for esplion_y B_q (3, 1:2:11) = DGH_q (2, 1:6) ; % for gamma B_q (3, 2:2:12) = DGH_q (1, 1:6) ; % for gamma R = xy_q (1) ; % radius B_q (4, 1:2:11) = H_q (1:6) / R ; % for hoop esplion Strains = B_q * T_e' ; % strains at point, 4 x 1 fprintf ('Element, Point, Eps_x, Eps_y, Gamma, Eps_h ... %i %i %g %g %g %g \n', j, q, Strains (1:4)) Stresses = E_q * Strains ; % stresses at point, 4 x 1 fprintf ('Element, Point, Sig_x, Sig_y, Tau, Sig_h ... %i %i %g %g %g %g \n', j, q, Stresses (1:4)) fprintf (fid, '%g %g %g %g %g %g \n', xy_q (1:2), Stresses (1:4)) %save show = 0 ; % skip the showing of element reaction results if ( show ) % list reactions % ELEMENT MATRICES UPDATES, 12 x 12 and 12 x 1 S_e = S_e + (B_q' * E_q * B_q) * t_q * J_det * w_q (q) ; % square if ( any ([Omega, Axial_a]) ) % then need acceleration loads Body_x = Rho * R * Omega^2 ; % centrifugal load Body_y = -Rho * Axial_a ; % axial load C_e(1:2:11) = C_e(1:2:11) + H_q' * Body_x * R * J_det * w_q (q); % Fx C_e(2:2:12) = C_e(2:2:12) + H_q' * Body_y * R * J_det * w_q (q); % Fy end % if body forces active end % if show end % for quadratures <--- <--- <--- <--- <--- <--- <--- <--- <--- if ( show ) % Assign any point source to the first element with that node if ( load_pt ) % then point sources to assign (once) for j = 1:n_i k = rows (j) ; if ( C_react (k) ~= 0 ) C_e = C_e + C_react (k) ; C_react (k) = 0 ; end % if end % for j node end %if load fprintf ('\nNode List for Element %g, \n', j) disp (e_nodes) fprintf ('Given Resultant Sources \n') disp (C_e') % Finally, get the reactions C_m = S_e * T_e' - C_e ; fprintf ('Element Nodal Reactions \n') disp (C_m') end % if show end % for each j element in mesh <<==== <<==== <<==== <<==== <<==== % End finite element calculations. % See /home/mech517/public_html/Matlab_Plots for graphic options % http://www.owlnet.rice.edu/~mech517/help_plot.html for help % end of Axisymmetric_Stress_T6_Isoparametric % +++++++++++++ functions in alphabetical order +++++++++++++++++ function echo_user_remarks file_id = fopen('msh_remarks.tmp', 'r'); % open for read only if ( file_id > 0 ) % then file exists fprintf('\n(Echo of msh_remarks.tmp) \n') while feof(file_id) ==0 % read until an end of file remark_line = fgets(file_id) ; % read a string and append \n fprintf (remark_line) ; % echo the line end % while data exists end % if file supplied % end function echo_user_remarks function [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C) % modify system linear eqs for essential boundary conditions % (by trick to avoid matrix partitions, loses reaction data) n_d = size (C, 1) ; % number of DOF eqs if ( size (EBC_flag, 2) > 1 ) ; % change to vector copy flag_EBC = reshape ( EBC_flag', 1, n_d) ; value_EBC = reshape ( EBC_value', 1, n_d) ; else flag_EBC = EBC_flag ; value_EBC = EBC_value ; end % if for j = 1:n_d % check all DOF for essential BC if ( flag_EBC (j) ) % then EBC here % Carry known columns*EBC to RHS. Zero that column and row. % Insert EBC identity, 1*EBC_dof = EBC_value. EBC = value_EBC (j) ; % recover EBC value C (:) = C (:) - EBC * S (:, j) ; % carry known column to RHS S (:, j) = 0 ; S (j, :) = 0 ; % clear, restore symmetry S (j, j) = 1 ; C (j) = EBC ; % insert identity into row end % if EBC for this DOF end % for over all j-th DOF % end enforce_essential_BC (EBC_flag, EBC_value, S, C) function [C] = get_and_add_point_sources (n_g, n_m, C) load msh_load_pt.tmp ; % node, DOF, value (eq. number) n_u = size(msh_load_pt, 1) ; % number of point sources if ( n_u < 1 ) ; % missing data error ('No load_pt data in msh_load_pt.tmp') end % if user error fprintf ('\nRead %g point sources. \n', n_u) fprintf ('(Echo of file msh_load_pt.tmp) \n') fprintf ('Node, DOF (1=force, 2=couple), Source_value \n') for j = 1:n_u ; % non-zero Neumann pts node = msh_load_pt (j, 1) ; % global node number DOF = msh_load_pt (j, 2) ; % local DOF number value = msh_load_pt (j, 3) ; % point source value fprintf ('%g %g %g \n', node, DOF, value) Eq = n_g * (node - 1) + DOF ; % row in system matrix C (Eq) = C (Eq) + value ; % add to system column matrix end % for each EBC % end get_and_add_point_sources (n_g, n_m, C) function [EBC_flag] = get_ebc_flags (n_g, n_m, P) EBC_flag = zeros(n_m, n_g) ; % initialize for k = 1:n_m ; % loop over all nodes if ( P(k) > 0 ) ; % at least one EBC here [flags] = unpack_pt_flags (n_g, k, P(k)) ; % unpacking EBC_flag (k, 1:n_g) = flags (1:n_g) ; % populate array end % if EBC at node k end % for loop over all nodes % end get_ebc_flags function [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) EBC_value = zeros(n_m, n_g) ; % initialize to zero load msh_ebc.tmp ; % node, DOF, value (eq. number) n_c = size(msh_ebc, 1) ; % number of constraints fprintf ('\nApplied Displacement Boundary Conditions: %g \n', n_c) fprintf ('(Echo of file load msh_ebc.tmp) \n') %b fprintf ('Node, DOF (1=displacement, 2=slope), Value. \n') fprintf ('Node, DOF (1=u, 2=v, 3=r), Value. \n') disp(msh_ebc) ; % echo input for j = 1:n_c ; % loop over ebc inputs node = round (msh_ebc (j, 1)) ; % node in mesh DOF = round (msh_ebc (j, 2)) ; % DOF # at node value = msh_ebc (j, 3) ; % EBC value % Eq = n_g * (node - 1) + DOF ; % row in system matrix EBC_value (node, DOF) = value ; % insert value in array if ( EBC_flag (node, DOF) == 0 ) % check data consistency fprintf ('WARNING: EBC but no flag at node %g & DOF %g. \n', ... node, DOF) % EBC_flag (node, DOF) = 1; % try to recover from data error end % if common user error end % for each EBC EBC_count = sum (sum ( EBC_flag > 0 )) ; % check input data if ( EBC_count ~= n_c ) ; % probable user error fprintf ('\nWARNING: mismatch in bc_flag count & msh_ebc.tmp \n') end % if user error % end get_ebc_values function [rows] = get_element_index (n_g, n_n, e_nodes) % calculate system DOF numbers of element, for gather, scatter rows = zeros (1, n_g*n_n) ; % allow for node = 0 for k = 1:n_n ; % loop over element nodes global_node = round (e_nodes (k)) ; % corresponding sys node for i = 1:n_g ; % loop over DOF at node eq_global = i + n_g * (global_node - 1) ; % sys DOF, if any eq_element = i + n_g * (k - 1) ; % el DOF number if ( eq_global > 0 ) ; % check node=0 trick rows (1, eq_element) = eq_global ; % valid DOF > 0 end % if allow for omitted nodes end % for DOF i % end local DOF loop end % for each element node % end local node loop % end get_element_index function [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements () ; % input file controls (for various data generators) load msh_typ_nodes.tmp ; % el_type, connectivity list (3) n_e = size (msh_typ_nodes,1) ; % number of elements if ( n_e == 0 ) ; % data file missing error ('Error missing file msh_typ_nodes.tmp') end % if error n_n = size (msh_typ_nodes,2) - 1 ; % nodes per element fprintf ('(Echo of file msh_typ_nodes.tmp) \n') fprintf ('Read %g elements with (ignored) type & %g nodes each. \n', ... n_e, n_n) el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1 n_t = max(el_type) ; % number of element types %b fprintf ('Maximum number of element types = %g. \n', n_t) nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n); disp(msh_typ_nodes (:, 1:1+n_n)) % echo data % end get_mesh_elements function [n_m, n_s, P, x, y, z] = get_mesh_nodes () ; % input file controls (for various data generators) % READ MESH AND EBC_FLAG INPUT DATA % specific problem data from MODEL data files (sequential) load msh_bc_xyz.tmp ; % bc_flag, x-, y-, z-coords n_m = size (msh_bc_xyz,1) ; % number of nodal points in mesh if ( n_m == 0 ) ; % data missing ! error ('Error missing file msh_bc_xyz.tmp') end % if error n_s = size (msh_bc_xyz,2) - 1 ; % number of space dimensions fprintf ('Read %g nodes. \n', n_m) fprintf ('(Echo of file msh_bc_xyz.tmp) \n') fprintf ('bc_flag, x-, y-coordinates \n') msh_bc_xyz (:, 1)= round (msh_bc_xyz (:, 1)); P = msh_bc_xyz (1:n_m, 1) ; % integer Packed BC flag x = msh_bc_xyz (1:n_m, 2) ; % extract x column y (1:n_m, 1) = 0.0 ; z (1:n_m, 1) = 0.0 ; % default to zero if (n_s > 1 ) ; % check 2D or 3D y = msh_bc_xyz (1:n_m, 3) ; % extract y column end % if 2D or 3D if ( n_s == 3 ) ; % check 3D z = msh_bc_xyz (1:n_m, 4) ; % extract z column end % if 3D for j = 1:n_m %bfprintf (' %2.2i %g \n', P(j), x(j) ) %bfprintf (' %1.1i %g %g \n', P(j), x(j), y(j) ) fprintf (' %2.2i %g %g \n', P(j), x(j), y(j) ) %bfprintf (' %3.3i %g %g \n', P(j), x(j), y(j) ) %bfprintf (' %2.2i %g %g \n', P(j), x(j), y(j) ) %bfprintf (' %3.3i %g %g %g \n', P(j), x(j), y(j), z(j) ) end % for % end get_mesh_nodes function list_save_beam_displacements (n_g, n_m, T) fprintf('\nCalculated Displacements: \n') fprintf('Node, Y_displacement, Z_rotation at %g nodes \n', n_m) T_matrix = reshape (T, n_g, n_m)' ; % pretty shape % save results (displacements) to MODEL file: node_results.tmp fid = fopen('node_results.tmp', 'w') ; % open for writing for j = 1:n_m ; % node loop, save displ fprintf (fid, '%g %g \n', T_matrix (j, 1:n_g)) ; % to file fprintf ('%g %g %g \n', j, T_matrix (j, 1:n_g)) ; % to screen end % for j DOF % end list_save_beam_displacements (n_g, n_m, T) function list_save_displacements_results (n_g, n_m, T) fprintf('\nCalculated Displacements: \n') fprintf('Node X_disp Y_disp \n') %b fprintf('Node X_disp Y_disp Z_disp \n') %b fprintf('Node U_disp V_disp Rotation \n') T_matrix = reshape (T, n_g, n_m)' ; % pretty shape disp (T_matrix) ; % print displacements % save results (displacements) to MODEL file: node_results.tmp fid = fopen('node_results.tmp', 'w') ; % open for writing for j = 1:n_m ; % save displacements if ( n_g == 1 ) fprintf (fid, '%i %g \n', j, T_matrix (j, 1:n_g)) ; elseif ( n_g == 2 ) fprintf (fid, '%i %g %g \n', j, T_matrix (j, 1:n_g)) ; elseif ( n_g == 3 ) fprintf (fid, '%i %g %g %g \n', j, T_matrix (j, 1:n_g)); elseif ( n_g == 4 ) fprintf (fid, '%i %g %g %g %g \n', j, T_matrix (j, 1:n_g)) ; elseif ( n_g == 5 ) fprintf (fid, '%i %g %g %g %g %g \n', j, T_matrix (j, 1:n_g)) ; elseif ( n_g == 6 ) fprintf (fid, '%i %g %g %g %g %g %g \n', j, T_matrix (j, 1:n_g)) ; else error ('reformat list_save_displacements_results for n_g > 6.') end % if end % for j DOF % end list_save_displacements_results (T) function list_save_temperature_results (T) n_m = size (T, 1) ; % get size fprintf('Temperature at %g nodes \n', n_m) ; % header % save results (temperature) to MODEL file: node_results.tmp fid = fopen('node_results.tmp', 'w') ; % open for writing for j = 1:n_m ; % save temperature fprintf ( fid, '%g \n', T (j)) ; % print fprintf (' %g %g \n', j, T (j)) ; % sequential save end % for j DOF % end list_save_temperature_results (T) function output_T3_heat_flux (n_e, n_g, n_n, n_q, nodes, x, y, T) % POST-PROCESS ELEMENT HEAT FLUX RECOVERY & SAVE fid = fopen('el_qp_xyz_fluxes.tmp', 'w') ; % open for writing fprintf ('\n') ; % blank line fprintf('Elem, X_qp, Y_qp, HeatFlux_x, HeatFlux_y \n');% header for j = 1:n_e ; % loop over elements ====>> e_nodes = nodes (j, 1:n_n) ; % connectivity [a, b, c, center, two_A] = form_T3_geom_constants (x, y, e_nodes); % get DOF numbers for this element, gather solution [rows] = get_element_index (n_g, n_n, e_nodes) ; % eq numbers T_e = T (rows) ; % gather element DOF for q = 1:n_q ; % Loop over element quadrature points ----> % H_i (x,y) = (a_i + b_i*x + c_i*y)/two_A % interpolations B_e (1, 1:3) = b(1:3) / two_A ; % dH/dx B_e (2, 1:3) = c(1:3) / two_A ; % dH/dy % COMPUTE GRADIENT & HEAT FLUX, SAVE LOCATION AND VALUES Gradient = B_e * T_e ; % gradient vector HeatFlux = E_e * Gradient ; % heat flux vector fprintf (fid, '%g %g %g %g \n', center(1:2), HeatFlux(1:2));% save fprintf ('%g %g %g %g %g \n', j, center(1:2), HeatFlux(1:2));% prt end % for loop over n_q element quadrature points <---- end % for each j element in mesh <<==== % end output_T3_heat_flux function [EBC_react] = recover_reactions_print_save (n_g, n_d, ... EBC_flag, EBC_row, EBC_col, T) % get EBC reaction values by using rows of S & C (before EBC) n_d = size (T, 1) ; % number of system DOF % n_c x 1 = n_c x n_d * n_d x 1 + n_c x 1 EBC_react = EBC_row * T - EBC_col ; % matrix reactions (+-) % save reactions (forces) to MODEL file: node_reaction.tmp fprintf ('\nRecovered Reactions at BC: %g \n', ... sum (sum (EBC_flag > 0))) ; % header %b fprintf ('Node, DOF (1=force, 2=couple), Value \n') fprintf ('Node, DOF, Value \n') fid = fopen('node_reaction.tmp', 'w') ; % open for writing if ( size (EBC_flag, 2) > 1 ) ; % change to vector copy flag_EBC = reshape ( EBC_flag', 1, n_d) ; % changed else flag_EBC = EBC_flag ; % original vector end % if Totals = zeros (1, n_g) ; % zero input totals kount = 0 ; % initialize counter for j = 1:n_d ; % extract all EBC reactions if ( flag_EBC(j) ) ; % then EBC here % Output node_number, component_number, value, equation_number kount = kount + 1 ; % copy counter node = ceil(j/n_g) ; % node at DOF j j_g = j - (node - 1)*n_g ; % 1 <= j_g <= n_g React = EBC_react (kount, 1) ; % reaction value fprintf ( fid, '%g %g %g \n', node, j_g, React);% save fprintf ('%g %g %g \n', node, j_g, React); % print Totals (j_g) = Totals (j_g) + React ; % sum all components end % if EBC for this DOF end % for over all j-th DOF fprintf ('Total force and couple = ') ; disp(Totals) ; % echo total % end recover_reactions_print_save (EBC_row, EBC_col, T) function [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C) n_d = size (C, 1) ; % number of system DOF EBC_count = sum (sum (EBC_flag)) ; % count EBC & reactions EBC_row = zeros(EBC_count, n_d) ; % reaction data EBC_col = zeros(EBC_count, 1) ; % reaction data if ( size (EBC_flag, 2) > 1 ) ; % change to vector copy flag_EBC = reshape ( EBC_flag', 1, n_d) ; % changed else flag_EBC = EBC_flag ; % original vector end % if kount = 0 ; % initialize counter for j = 1:n_d % System DOF loop, check for displacement BC if ( flag_EBC (j) ) ; % then EBC here % Save reaction data to be destroyed by EBC solver trick kount = kount + 1 ; % copy counter EBC_row(kount, 1:n_d) = S (j, 1:n_d) ; % copy reaction data EBC_col(kount, 1) = C (j) ; % copy reaction data end % if EBC for this DOF end % for over all j-th DOF % end sys DOF loop % end save_reaction_matrices (S, C, EBC_flag) function save_resultant_load_vectors (n_g, C) % save resultant forces to MODEL file: node_resultants.tmp n_d = size (C, 1) ; % number of system DOF fprintf ('\nResultants of all input sources: \n') fprintf ('Node, DOF (1=force, 2=couple), Value \n') fid = fopen('node_resultant.tmp', 'w'); % open for writing Totals = zeros (1, n_g) ; % zero input totals for j = 1:n_d ; % extract all resultants if ( C (j) ~= 0. ) ; % then source here % Output node_number, component_number, value, equation_number node = ceil(j/n_g) ; % node at DOF j j_g = j - (node - 1)*n_g ; % 1 <= j_g <= n_g value = C (j) ; % resultant value fprintf ( fid, '%g %g %g %g \n', node, j_g, value, j);% save fprintf ('%g %g %g \n', node, j_g, value); % print Totals (j_g) = Totals (j_g) + value ; % sum all inputs end % if non-zero for this DOF end % for over all j-th DOF fprintf ('Totals = ') ; disp(Totals) ; % echo totals % end save_resultant_load_vectors (n_g, n_m, C) function [flags] = unpack_pt_flags (n_g, N, flag) % unpack n_g integer flags from the n_g digit flag at node N % integer flag contains (left to right) f_1 f_2 ... f_n_g full = flag ; % copy integer check = 0 ; % validate input for Left2Right = 1:n_g ; % loop over local DOF at k Right2Left = n_g + 1 - Left2Right ; % reverse direction work = floor (full / 10) ; % work item keep = full - work * 10 ; % work item flags (Right2Left) = keep ; % insert into array full = work ; % work item check = check + keep * 10^(Left2Right - 1) ; % validate end % for each local DOF if ( flag > check ) ; % check for likely error fprintf ('WARNING: bc flag likely reversed at node %g. \n', N) end % if likely user error % end unpack_pt_flags %=================== Running gives =====================================