function L3_Quintic_BoEF (load_pt) % load_pt = 1 flags that point forces and/or couples will be % input via file msh_load_pt.tmp %............................................................. % quintic beam on elastic foundation, % for point loads & couples, line load % Mech 417, Finite Element Analysis, Rice University % Copyright J.E. Akin 2008. All rights reserved. %............................................................. if ( nargin == 0 ) ; % check for optional data flag load_pt = 0 ; % no point source data end % if from argument count % Application and element dependent controls n_g = 2 ; % number of DOF per node (deflection, slope) n_q = 0 ; % number of quadrature points required n_r = 1 ; % number of rows in B_e matrix % Read mesh nodal 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 if ( EBC_count == 0 ) fprintf ('WARNING: No displacement boundary condition given. \n') else fprintf ('Note: expecting %g displacement BC values. \n', EBC_count) end % if % Read mesh element input data files [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) ; M = zeros (n_d, n_d) ; % stiff & mass C = zeros (n_d, 1) ; % force & moments T = zeros (n_d, 1) ; % displacements % 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 loads & couples data, if any, and 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 fprintf ('(Echoing file msh_properties.tmp) \n') n_prop = size(msh_properties, 1) ; % == 1 if homogeneous F_k (1:n_prop) = msh_properties (1:n_prop, 5) ; % extract foundation k's if ( any ( F_k > 0 )) % BEF pressures exist BEF = 1 ; fprintf ('WARNING: BoEF moment and shear inaccurate for crude meshes \n') else BEF = 0 ; end % if reaction flag if ( n_prop == 1 ) E = msh_properties (1, 1) ; % material modulus I = msh_properties (1, 2) ; % section inertia Line_e (1:2) = msh_properties (1, 3:4) ; % trapezoidal line load k_f = msh_properties (1, 5) ; % foundation stiffness Rho = msh_properties (1, 6) ; % beam mass per unit length % ECHO PROPERTIES fprintf (' \n') fprintf ('Homogeneous Element Properties: \n' ) fprintf ('Elastity modulus (N/m^2) = %g \n', E) fprintf ('Moment of inertia (m^4) = %g \n', I) fprintf ('Line Load (N/m) = [ %g %g ] \n', ... Line_e(1), Line_e(2)) fprintf ('Foundation stiffness (N/m^2) = %g \n', k_f) fprintf ('Mass per unit length (Kg/m) = %g \n', Rho) end % if if ( n_prop > n_e ) error ('ERROR: number of property sets exceeds number of elements') end % if % 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 ; % loop over elements ====>> ====>> ====>> ====>> S_e = zeros (n_i, n_i) ; % clear array NtN = zeros (n_i, n_i) ; % clear array M_e = zeros (n_i, n_i) ; % clear array C_e = zeros (n_i, 1) ; % clear arrays e_nodes = nodes (j, 1:n_n) ; % connectivity % SET ELEMENT PROPERTIES & GEOMETRY L = x(e_nodes(n_n)) - x(e_nodes(1)) ; % beam length Line_e (1:2) = 0 ; if ( n_prop > 1 ) E = msh_properties (j, 1) ; % material modulus I = msh_properties (j, 2) ; % section inertia Line_e (1:2) = msh_properties (j, 3:4) ; % trapezoidal line load k_f = msh_properties (j, 5) ; % foundation stiffness Rho = msh_properties (j, 6) ; % beam mass per unit length % ECHO PROPERTIES fprintf ('\nProperties for element %g \n', j) fprintf ('Elastity modulus (N/m^2) = %g \n', E) fprintf ('Moment of inertia (m^4) = %g \n', I) fprintf ('Line Load (N/m) = [ %g %g ] \n', ... Line_e(1), Line_e(2)) fprintf ('Foundation stiffness (N/m^2) = %g \n', k_f) fprintf ('Mass per unit length (km/m) = %g \n', Rho) end % if % ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES L2 = L^2 ; L3 = L^3 ; % Constant cubic element & foundation stiffnesses & mass matrices S_e = [5092, 1138*L, -3584, 1920*L, -1508, 242*L ; 1138*L, 332*L2, -896*L, 320*L2, -242*L, 38*L2 ; -3584, -896*L, 7168, 0, -3584, 896*L ; 1920*L, 320*L2, 0, 1280*L2, -1920*L, 320*L2 ; -1508, -242*L, -3584, -1920*L, 5092, -1138*L ; 242*L, 38*L2, 896*L, 320*L2, -1138*L, 332*L2 ] ... *E*I/(35 * L^3) ; % stiffness NtN = L*[2092, 114*L, 880, -160*L, 262, -29*L ; 114*L, 8*L2, 88*L, -12*L2, 29*L, -3*L2 ; 880, 88*L, 5632, 0, 880, -88*L ; -160*L, -12*L2, 0, 128*L2, 160*L, -12*L2 ; 262, 29*L, 880, 160*L, 2092, -114*L ; -29*L, -3*L2 -88*L, -12*L2, -114*L, 8*L2 ]/13860; S_e = S_e + k_f * NtN ; % add foundation stiffness before assembly M_e = Rho * NtN ; % mass matrix % Map line load to node forces & moments; C_e = p_To_F * Line_e C_e = zeros (n_i, 1) ; % clear arrays if ( any (Line_e) ) ; % then form forcing vector p_To_F = L * [ 79, 19 ; 5*L, 2*L ; 112, 112 ; -8*L, 8*L ; 19, 79 ; -2*L, -5*L ] / 420 ; % cubic H, linear Line % 6 x 2 * 2 x 1 = 6 x 1 result C_e = p_To_F (1:n_i, 1:2) * Line_e' ; % force moment @ nodes end % if or set up resultant node loads for line load % 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 stiffness M (rows, rows) = M (rows, rows) + M_e ; % add to system mass C (rows) = C (rows) + C_e ; % add to sys force/couples end % for each j element in mesh <<==== <<==== <<==== <<==== <<==== % ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY if ( EBC_count > 0 ) ; % reactions occur [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C); save_resultant_load_vectors (n_g, C) end % if essential BC exist (almost always true) % ENFORCE ESSENTIAL BOUNDARY CONDITIONS if ( EBC_count > 0 ) ; % reactions occur [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C); end % if essential BC exist (almost always true) % COMPUTE SOLUTION & SAVE T = S \ C ; % Compute displacements & rotations list_save_beam_displacements (n_g, n_m, T) ; % save and print % GRAPHS if ( BEF == 1 ) % plot foundation pressure C1_L3_BoEF_pressure (F_k) % graph pause (5) % display for 5 seconds end % if C1_L3_beam_plots(1) % deflection pause (5) % display for 5 seconds C1_L3_beam_plots(2) % slope pause (5) % display for 5 seconds C1_L3_beam_plots(3) % moment diagram pause (5) % display for 5 seconds C1_L3_beam_plots(4) % transverse shear diagram pause (5) % display for 5 seconds % 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 ELEMENT REACTIONS (MEMBER FORCES) %b if ( BEF == 0 ) fprintf ('\nIndividual Element Load and Reaction Summaries: \n') fprintf ('(F_1, M_1, F_2, M_2) \n') for j = 1:n_e ; % loop over 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 L = x(e_nodes(n_n)) - x(e_nodes(1)) ; % beam length Line_e (1:2) = 0 ; if ( n_prop > 1 ) E = msh_properties (j, 1) ; % material modulus I = msh_properties (j, 2) ; % section inertia Line_e (1:2) = msh_properties (j, 3:4) ; % trapezoidal line load k_f = msh_properties (j, 5) ; % foundation stiffness end % if % ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES L2 = L^2 ; L3 = L^3 ; % Constant cubic element & foundation stiffnesses & mass matrices S_e = [5092, 1138*L, -3584, 1920*L, -1508, 242*L ; 1138*L, 332*L2, -896*L, 320*L2, -242*L, 38*L2 ; -3584, -896*L, 7168, 0, -3584, 896*L ; 1920*L, 320*L2, 0, 1280*L2, -1920*L, 320*L2 ; -1508, -242*L, -3584, -1920*L, 5092, -1138*L ; 242*L, 38*L2, 896*L, 320*L2, -1138*L, 332*L2 ] ... *E*I/(35 * L^3) ; % stiffness NtN = L*[2092, 114*L, 880, -160*L, 262, -29*L ; 114*L, 8*L2, 88*L, -12*L2, 29*L, -3*L2 ; 880, 88*L, 5632, 0, 880, -88*L ; -160*L, -12*L2, 0, 128*L2, 160*L, -12*L2 ; 262, 29*L, 880, 160*L, 2092, -114*L ; -29*L, -3*L2 -88*L, -12*L2, -114*L, 8*L2 ]/13860; % Map line load to node forces & moments; C_e = p_To_F * Line_e C_e = zeros (n_i, 1) ; % clear arrays if ( any (Line_e) ) ; % then form forcing vector p_To_F = L * [ 79, 19 ; 5*L, 2*L ; 112, 112 ; -8*L, 8*L ; 19, 79 ; -2*L, -5*L ] / 420 ; % cubic H, linear Line % 6 x 2 * 2 x 1 = 6 x 1 result C_e = p_To_F (1:n_i, 1:2) * Line_e' ; % force moment @ nodes end % if or set up resultant node loads for line load % Assign any point force/couple to the first element with that node if ( j == 1 ) % first element only gets sources from 3 nodes C_e = C_e + C_react (rows) ; else % last 2 nodes C_e (3:6) = C_e (3:6) + C_react (rows(3:6)) ; % quintic C1 beam ONLY !! end % if fprintf ('\nGiven Resultant Loading on Member %g, \n', j) disp (C_e') % Finally, get the reactions C_m = S_e * T_e' - C_e ; if ( BEF ) % add foundation resultants C_f = -k_f * NtN * T_e' ; fprintf ('Resultant Foundation Loading on Member %g, \n', j) disp (C_f') C_m = C_m - C_f ; end % if fprintf ('Net Resultant End Reactions on Member %g, \n', j) disp (C_m') end % for each j element in mesh <<==== <<==== <<==== <<==== <<==== %b end % if no BEF % 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 L3_Quintic_BoEF % +++++++++++++ functions in alphabetical order +++++++++++++++++ function C1_L3_BoEF_pressure (F_k) % Copyright 2008 J.E. Akin. All rights reserved. % ------------------------------------------------------ % Matlab graph of foundation pressurelue at mesh nodes % for a mesh of 3 node quintic Hermite line elements % ------------------------------------------------------ % c_x = x coordinates of nod_per_el line element % msh_typ_nodes = connectivity list , 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 if ( nargin == 0 ) i_p = 1 ; end % if no arguments n_fk = size (F_k, 1) ; % number of different foundations Pts_wide = 2 ; % fat lines format short % 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 %b fprintf ('Read %g mesh coordinates \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 %b fprintf ('Read %g elements connections \n', nt) if ( nod_per_el ~= 3 ) error ('This is not a mesh of 3 node line elements') end % if 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 %b fprintf ('Read %g nodal solution values \n', nr) %b fprintf (' with %g components each \n', max_p) H (3) = 0. ; HC1 (6) = 0. ; DHC1 (6) = 0. ; D2HC1 (6) = 0. ; D3HC1 (6) = 0. ; D4HC1 (6) = 0. ; x (np) = 0. ; % pre-allocate array x y (np) = 0. ; % pre-allocate array y dy (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 t_dy (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 xmax = max (x) ; xmin = min (x) ; y = node_results(:, 1) ; dy = node_results(:, 2) ; clf % clear graphics hold on % hold image for plots xlabel (['X, Node at 45 deg (', int2str(nod_per_el), ... ' per element), Element at 90 deg']) % Are properties needed ? load msh_properties.tmp n_prop = size (msh_properties, 2) ; % Loop over all elements max_all = -1e9 ; min_all = 1e9 ; for it = 1:nt ; % Element foundation stiffness, if any if ( n_fk == 1 ) % then homogeneous k_f = F_k (1) ; else k_f = F_k (it) ; end % if % 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 %b if ( t_x(2) ~= (t_x(1)+t_x(3))/2 ) %b fprintf ('Bad Jaconian element %g \n', it) %b end % if non-constant Jacobian t_y = y (t_nodes) ; % y at those nodes, only t_dy = dy (t_nodes) ; % dy at those nodes, only D (1:2:6) = t_y ; D (2:2:6) = t_dy ; % Loop over local points on the quadratic polynomial element n_poly = ceil ( 95 / nt) ; for k = 1: (n_poly + 1) % points in parametric space % get element parametric interpolation functions R = (k - 1)/n_poly ; % on 0 to 1 X = 2*R - 1 ; % on -1 to 1 % H = ELEMENT SHAPE FUNCTIONS % X = LOCAL COORDINATE OF POINT, -1 TO +1 % LOCAL NODE COORD. ARE -1,0,+1 1-----2-----3 H (1) = 0.5*(X*X - X) ; H (2) = 1. - X*X ; H (3) = 0.5*(X*X + X) ; x_el (k) = H * t_x ; % true x value A = t_x(3) - t_x(1) ; P = X ; P_2 = P * P ; P_3 = P * P_2 ; P_4 = P * P_3 ; P_5 = P * P_4 ; % deflection HC1(1) = (4*P_2 - 5*P_3 - 2*P_4 + 3*P_5) * 0.25 ; HC1(2) = (P_2 - P_3 - P_4 + P_5) * 0.125 * A ; HC1(3) = 1 - 2*P_2 + P_4 ; HC1(4) = (P - 2*P_3 + P_5) * A * 0.5 ; HC1(5) = (4*P_2 + 5*P_3 - 2*P_4 - 3*P_5) * 0.25 ; HC1(6) = (-P_2 - P_3 + P_4 + P_5) * 0.125 * A ; r_el (k) = HC1 * D' ; % true y value end % for k plot points r_el = -r_el* msh_properties(it,5) ; % elem max, min values [V_X, L_X] = max (r_el) ; [V_N, L_N] = min (r_el); if ( V_X > max_all ) max_all = V_X ; end if ( V_N < min_all ) min_all = V_N ; end plot (x_el, r_el, 'b-', 'LineWidth',Pts_wide) end % if zero in connectivity % Plot the element number x_bar = sum (t_x' )/nod_per_el ; y_bar = mean (r_el) ; t_text = sprintf (' (%g)', it); % offset # from pt text (x_bar, y_bar, t_text) % incline end % for over all elements fprintf ('Max foundation pressure value is %g \n', max_all) fprintf ('Min foundation pressure value is %g \n', min_all) null (1:np) = 0.5*(V_N + V_X) ; if ( V_N <= 0 & V_X >= 0 ) null (1:np) = 0. ; end % if % finalize axes ymax = max_all ; ymin = min_all ; diff = abs(ymax-ymin) ; ymax = ymax + abs (diff)/10. ; ymin = ymin - abs (diff)/10. ; axis ([xmin, xmax, ymin, ymax]) % set axes title(['Beam Foundation Pressure from: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) ylabel (['Pressure: max=', num2str(max_all), ... ', min=', num2str(min_all)]) % 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 % -depsc -tiff % for an eps version print ('-dpng', ['BoEF_pressure']) v_text = ['Created BoEF_pressure.png'] ; fprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) hold off % end of C1_L3_BoEF_pressure function C1_L3_beam_plots (i_p) % Copyright 2008 J.E. Akin. All rights reserved. % ------------------------------------------------------ % Matlab graph of i_p-th component value at mesh nodes % for a mesh of 3 node quintic Hermite line elements % i_p = 1 is deflection, = 2 is slope % = 3 is curvature (moment), = 4 is curvature rate (shear) % ------------------------------------------------------ % c_x = x coordinates of nod_per_el line element % msh_typ_nodes = connectivity list , 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 if ( nargin == 0 ) i_p = 1 ; end % if no arguments Pts_wide = 2 ; % draw fat lines format short % 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 %b fprintf ('Read %g mesh coordinates \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 %b fprintf ('Read %g elements connections \n', nt) if ( nod_per_el ~= 3 ) error ('This is not a mesh of 3 node line elements') end % if 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 %b fprintf ('Read %g nodal solution values \n', nr) %b fprintf (' with %g components each \n', max_p) if ( i_p > 4 ) fprintf ('Data requested for component i_p = %g \n', i_p) error ('i_p > available data') end % if error H (3) = 0. ; HC1 (6) = 0. ; DHC1 (6) = 0. ; D2HC1 (6) = 0. ; D3HC1 (6) = 0. ; D4HC1 (6) = 0. ; x (np) = 0. ; % pre-allocate array x y (np) = 0. ; % pre-allocate array y dy (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 t_dy (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 xmax = max (x) ; xmin = min (x) ; y = node_results(:, 1) ; dy = node_results(:, 2) ; clf % clear graphics hold on % hold image for plots xlabel (['X, Node at 45 deg (', int2str(nod_per_el), ... ' per element), Element at 90 deg']) % Are properties needed ? if ( i_p > 2 ) load msh_properties.tmp n_prop = size (msh_properties, 2) ; else n_prop = 0; EI = 1 ; end % if % Loop over all elements max_all = -1e9 ; min_all = 1e9 ; 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 %b if ( t_x(2) ~= (t_x(1)+t_x(3))/2 ) %b fprintf ('Bad Jaconian element %g \n', it) %b end % if non-constant Jacobian t_y = y (t_nodes) ; % y at those nodes, only t_dy = dy (t_nodes) ; % dy at those nodes, only D (1:2:6) = t_y ; D (2:2:6) = t_dy ; if ( i_p == 1 ) plot (t_x, t_y, 'ko') % plot nodal value symbols end % if % Loop over local points on the quadratic polynomial element n_poly = ceil ( 95 / nt) ; for k = 1: (n_poly + 1) % points in parametric space % get element parametric interpolation functions R = (k - 1)/n_poly ; % on 0 to 1 X = 2*R - 1 ; % on -1 to 1 % H = ELEMENT SHAPE FUNCTIONS % X = LOCAL COORDINATE OF POINT, -1 TO +1 % LOCAL NODE COORD. ARE -1,0,+1 1-----2-----3 H (1) = 0.5*(X*X - X) ; H (2) = 1. - X*X ; H (3) = 0.5*(X*X + X) ; x_el (k) = H * t_x ; % true x value A = t_x(3) - t_x(1) ; P = X ; P_2 = P * P ; P_3 = P * P_2 ; P_4 = P * P_3 ; P_5 = P * P_4 ; % select result of interest, r_el if ( i_p == 1 ) % deflection HC1(1) = (4*P_2 - 5*P_3 - 2*P_4 + 3*P_5) * 0.25 ; HC1(2) = (P_2 - P_3 - P_4 + P_5) * 0.125 * A ; HC1(3) = 1 - 2*P_2 + P_4 ; HC1(4) = (P - 2*P_3 + P_5) * A * 0.5 ; HC1(5) = (4*P_2 + 5*P_3 - 2*P_4 - 3*P_5) * 0.25 ; HC1(6) = (-P_2 - P_3 + P_4 + P_5) * 0.125 * A ; r_el (k) = HC1 * D' ; % true y value elseif ( i_p == 2 ) % slope DHC1 (1) = (8*P - 15*P_2 - 8*P_3 + 15*P_4) * 0.5 / A ; DHC1 (2) = (2*P - 3*P_2 - 4*P_3 + 5*P_4) * 0.25 ; DHC1 (3) = (-4*P + 4*P_3) * 2 / A ; DHC1 (4) = (1 - 6*P_2 + 5*P_4) ; DHC1 (5) = (8*P + 15*P_2 - 8*P_3 - 15*P_4) * 0.5 / A ; DHC1 (6) = (-2*P - 3*P_2 + 4*P_3 + 5*P_4) * 0.25 ; r_el (k) = DHC1 * D' ; % true y' value elseif ( i_p == 3 ) % curvature or moment D2HC1 (1) = (8 - 30*P - 24*P_2 + 60*P_3) /A/ A ; D2HC1 (2) = (2 - 6*P - 12*P_2 + 20*P_3) * 0.5/A ; D2HC1 (3) = (-4 + 12*P_2) * 4 / A / A ; D2HC1 (4) = (-12*P + 20*P_3)*2 / A ; D2HC1 (5) = (8 + 30*P - 24*P_2 - 60*P_3) / A /A; D2HC1 (6) = (-2 - 6*P + 12*P_2 + 20*P_3) * 0.5 / A ; r_el (k) = D2HC1 * D' ; % true y'' value elseif ( i_p == 4 ) % curvature rate or shear force D3HC1 (1) = (-30 - 48*P + 180*P_2) * 2 /A/A/ A ; D3HC1 (2) = (-6 - 24*P + 60*P_2) /A /A; D3HC1 (3) = (24*P) * 8 / A / A/A ; D3HC1 (4) = (-12 + 60*P_2)*4 / A/A ; D3HC1 (5) = (30 - 48*P - 180*P_2) * 2 / A /A/A; D3HC1 (6) = (-6 + 24*P + 60*P_2) / A/A ; r_el (k) = D3HC1 * D' ; % true y''' value end % if plot type end % for k plot points if ( n_prop > 1 ) EI = msh_properties(it,1)*msh_properties(it,2) ; r_el = r_el* EI ; elseif ( n_prop == 1 ) EI = msh_properties( 1,1)*msh_properties( 1,2) ; r_el = r_el* EI ; end % if % elem max, min values [V_X, L_X] = max (r_el) ; [V_N, L_N] = min (r_el); if ( V_X > max_all ) max_all = V_X ; end if ( V_N < min_all ) min_all = V_N ; end plot (x_el, r_el, 'b-', 'LineWidth',Pts_wide) end % if zero in connectivity % Plot the element number x_bar = sum (t_x' )/nod_per_el ; y_bar = mean (r_el) ; t_text = sprintf (' (%g)', it); % offset # from pt text (x_bar, y_bar, t_text) % incline end % for over all elements null (1:np) = 0.5*(V_N + V_X) ; if ( V_N <= 0 & V_X >= 0 ) null (1:np) = 0. ; end % if % finalize axes ymax = max_all ; ymin = min_all ; diff = abs(ymax-ymin) ; ymax = ymax + abs (diff)/10. ; ymin = ymin - abs (diff)/10. ; axis ([xmin, xmax, ymin, ymax]) % set axes if ( i_p == 1 ) title(['Qunitic Beam Deflection from: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) ylabel (['Deflection: max=', num2str(max_all), ... ', min=', num2str(min_all)]) fprintf ('Max deflection value is %g \n', max_all) fprintf ('Min deflection value is %g \n', min_all) elseif ( i_p == 2 ) title(['Qunitic Beam Slope from: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) ylabel (['Slope: max=', num2str(max_all), ... ', min=', num2str(min_all)]) fprintf ('Max slope value is %g \n', max_all) fprintf ('Min slope value is %g \n', min_all) elseif ( i_p == 3 ) title(['Qunitic Beam Moment from: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) ylabel (['Moment: max=', num2str(max_all), ... ', min=', num2str(min_all)]) fprintf ('Max moment value is %g \n', max_all) fprintf ('Min moment value is %g \n', min_all) elseif ( i_p == 4 ) title(['Qunitic Beam Shear from: ', ... int2str(nt),' Elements, ', int2str(np),' Nodes']) ylabel (['Transverse Shear: max=', num2str(max_all), ... ', min=', num2str(min_all)]) fprintf ('Max shear value is %g \n', max_all) fprintf ('Min shear value is %g \n', min_all) end % if % 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 % -depsc -tiff % for an eps version if ( i_p == 1 ) print ('-dpng', ['Beam_deflections']) v_text = ['Created Beam_deflections.png'] ; elseif ( i_p == 2 ) print ('-dpng', ['Beam_slope']) v_text = ['Created Beam_slope.png'] ; elseif ( i_p == 3 ) print ('-dpng', ['Beam_moments']) v_text = ['Created Beam_moments.png'] ; elseif ( i_p == 4 ) print ('-dpng', ['Beam_shears']) v_text = ['Created Beam_shears.png'] ; end % if fprintf (1,'%s', v_text) ; fprintf (1, ' \n' ) hold off % end of C1_L3_beam_plots 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. (Echo of msh_load_pt.tmp) \n', n_u) 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 ('Read %g EBC data with Node, DOF, Value. \n', n_c) fprintf ('(Echo of file msh_ebc.tmp) \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 ('WARNING: mismatch in bc_flag count & msh_ebc.tmp') 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 ('Read %g elements with (ignored) type & %g nodes each. \n', ... n_e, n_n) fprintf ('(Echo of file msh_typ_nodes.tmp) \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-coordinate \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 %b disp(msh_bc_xyz (:, 1:1+n_s)) ; % echo data for j = 1:n_m fprintf (' %2.2i %g \n', P(j), x(j) ) %bfprintf (' %2.2i %g %g \n', P(j), x(j), y(j) ) %bfprintf (' %2.2i %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('\nNode, 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 ('\n') ; fprintf('X_disp Y_disp Z_disp at %g nodes \n', n_m) 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, '%g \n', T_matrix (j, 1:n_g)) ; elseif ( n_g == 2 ) fprintf (fid, '%g %g \n', T_matrix (j, 1:n_g)) ; elseif ( n_g == 3 ) fprintf (fid, '%g %g %g \n', T_matrix (j, 1:n_g)) ; elseif ( n_g == 4 ) fprintf (fid, '%g %g %g %g \n', T_matrix (j, 1:n_g)) ; elseif ( n_g == 5 ) fprintf (fid, '%g %g %g %g %g \n', T_matrix (j, 1:n_g)) ; elseif ( n_g == 6 ) fprintf (fid, '%g %g %g %g %g %g \n', 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 [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 ('\nReactions at Displacement BCs \n') fprintf ('Node, DOF, Reaction 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 ('Totals = ') ; disp(Totals) ; % echo totals % 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 ('\n') ; % Skip a line % fprintf ('Node, DOF, Resultant Force (1) or Moment (2) \n') fprintf ('Resultant input sources: \n') fprintf ('Node, DOF, Resultant input sources \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 ======================= % >> L3_Quintic_BoEF(1) % Read 11 nodes with bc_flag & 1 coordinates. % (Echo of file msh_bc_xyz.tmp) % 10 0 % 0 5 % 0 10 % 0 15 % 0 20 % 0 25 % 10 30 % 0 35 % 0 40 % 0 45 % 0 50 % % Note: expecting 2 displacement BC values. % % Read 5 elements with (ignored) type & 3 nodes each. % (Echo of file msh_typ_nodes.tmp) % 1 1 2 3 % 1 3 4 5 % 1 5 6 7 % 1 7 8 9 % 1 9 10 11 % % Read 2 EBC data % (Echo of file msh_ebc.tmp) % Node, DOF, Value. % 1 1 0 % 7 1 0 % % Read 2 point sources. (Echo of msh_load_pt.tmp) % Node, DOF (1=force, 2=couple), Source_value % 5 1 -100 % 9 2 3000 % % (Echoing file msh_properties.tmp) % Properties for element 1 % Elastity modulus (N/m^2) = 1e+09 % Moment of inertia (m^4) = 0.01 % Line Load (N/m) = [ 0 0 ] % Foundation stiffness (N/m^2) = 0 % % Properties for element 2 % Elastity modulus (N/m^2) = 1e+09 % Moment of inertia (m^4) = 0.01 % Line Load (N/m) = [ -10 -10 ] % Foundation stiffness (N/m^2) = 0 % % Properties for element 3 % Elastity modulus (N/m^2) = 1e+09 % Moment of inertia (m^4) = 0.01 % Line Load (N/m) = [ -10 -10 ] % Foundation stiffness (N/m^2) = 0 % % Properties for element 4 % Elastity modulus (N/m^2) = 1e+09 % Moment of inertia (m^4) = 0.01 % Line Load (N/m) = [ 0 0 ] % Foundation stiffness (N/m^2) = 0 % % Properties for element 5 % Elastity modulus (N/m^2) = 1e+09 % Moment of inertia (m^4) = 0.01 % Line Load (N/m) = [ 0 0 ] % Foundation stiffness (N/m^2) = 0 % % Resultant input sources: % Node, DOF, Resultant input sources % 3 1 -23.3333 % 3 2 -16.6667 % 4 1 -53.3333 % 5 1 -146.667 % 6 1 -53.3333 % 7 1 -23.3333 % 7 2 16.6667 % 9 2 3000 % Totals = -300 3000 % % Calculated Displacements % Node, Y_displacement, Z_rotation at 11 nodes % 1 0 -0.00272222 % 2 -0.0131944 -0.00247222 % 3 -0.0238889 -0.00172222 % 4 -0.0296094 -0.000493056 % 5 -0.0281944 0.00111111 % 6 -0.0182899 0.00284028 % 7 0 0.00444444 % 8 0.0259722 0.00594444 % 9 0.0594444 0.00744444 % 10 0.0966667 0.00744444 % 11 0.133889 0.00744444 % Max deflection value is 0.133889 % Min deflection value is -0.0300147 % Created Beam_deflections.png % % Max slope value is 0.00744448 % Min slope value is -0.00272222 % Created Beam_slope.png % % Max moment value is 1749.94 % Min moment value is -0.3 % Created Beam_moments.png % % Max shear value is 50.049 % Min shear value is -24.9615 % Created Beam_shears.png % % Reactions at Displacement BCs % Node, DOF, Reaction Value % 1 1 200 % 7 1 100 % Totals = 300.0000 0 % % Individual Element Load and Reaction Summaries: % (F_1, M_1, F_2, M_2) % Given Resultant Loading on Member 1, % 0 0 0 0 0 0 % Net Resultant End Reactions on Member 1, % 1.0e+03 * % 0.2000 -0.0000 0 0.0000 -0.2000 2.0000 % % Given Resultant Loading on Member 2, % -23.3333 -16.6667 -53.3333 0 -123.3333 16.6667 % Net Resultant End Reactions on Member 2, % 1.0e+03 * % 0.2000 -2.0000 -0.0000 0.0000 -0.0000 3.5000 % % Given Resultant Loading on Member 3, % -23.3333 -16.6667 -53.3333 0 -23.3333 16.6667 % Net Resultant End Reactions on Member 3, % 1.0e+03 * % 0.0000 -3.5000 -0.0000 0 0.1000 3.0000 % % Given Resultant Loading on Member 4, % 0 0 0 0 0 3000 % Net Resultant End Reactions on Member 4, % 1.0e+03 * % -0.0000 -3.0000 0.0000 -0.0000 0.0000 -0.0000 % % Given Resultant Loading on Member 5, % 0 0 0 0 0 0 % Net Resultant End Reactions on Member 5, % 0 0 0 0 0 0 % >> quit