How do i write this finite element code for a 2D truss beam in matlab? You need
ID: 1841937 • Letter: H
Question
How do i write this finite element code for a 2D truss beam in matlab?
You need to write a Finite Element Code for 2D truss structures using 1D Two-Noded Bar Elements. • Main.m file containing a function. You need to write the input to your code and test it.
% Call to function to import nodes, elements, Young’s moduli and cross section areas
[Nodes,Elements, E, A] = ReadFiles(‘elements.txt’,’nodes.txt’,’youngs.txt’,’areas.txt’);
% Plot undeformed configuration
PlotMesh(Nodes,Elements);
% Define here the boundary conditions
% Your code here: arrays F, Ux, Uy must contain your BCs in this way:
% F = [NodeID Fx Fy]
% Ux = [NodeID value: 0 if fixed, any number otherwise]
& Uy = [NodeID value] see Verification Example 1.
% Call to StiffnessMatrix function
K = StiffnessMatrix(Nodes,Elements,E,A);
% Call to Solve function
[d,R] = Solve(K,F,Ux,Uy);
% Post-processing – Plot Deformed configuration and compute stresses in each bar
Stress = PlotDeformed(Nodes,Elements,E,A,d);
Explanation / Answer
Matlab program:
function Modular_2D_Truss (load_pt)
%.............................................................
% Classic planar truss for point loads (& line load soon)
% XY COORDINATES CLOSED FORM INTEGRALS
%.............................................................
% pre_e = # of dummy items before el_type & connectivity
% pre_p = # of dummy items before BC_flag % coordinates
pre_e = 0 ; pre_p = 0 ; % default, consistent with plots
if ( nargin == 0 ) ; % check for optional data
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 (axial_d, transverse_d)
n_q = 0 ; % number of quadrature points required
n_r = 1 ; % number of rows in B_e matrix
% Read mesh input data files
[n_m, n_s, P, x, y, z] = get_mesh_nodes (pre_p) ;
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements (pre_e) ;
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
% 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
% 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 or moments, 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
% ============== ASSUMING HOMOGENOUS PROPERTIES =================
% GENERATE ELEMENT MATRICES AND ASSYMBLE INTO SYSTEM
% Assemble n_d by n_d square matrix terms from n_e by n_e
for i = 1:n_e % loop over elements ====>> ====>> ====>> ====>>
S_e = zeros (n_i, n_i) ; M_e = zeros (n_i, n_i) ; % sys arrays
C_p = zeros (n_i, 1) ; C_e = zeros (n_i, 1) ; % sys arrays
s_L = zeros (n_i, n_i) ; m_L = zeros (n_i, n_i) ; % loc arrays
t_L = zeros (n_i, n_i) ; c_L = zeros (n_i, 1) ; % loc arrays
e_nodes = nodes (i, 1:n_n) ; % connectivity
% SET ELEMENT PROPERTIES & GEOMETRY
Option = 1 ; % select analysis case
[A, E, Line_e, Rho] = set_2D_truss_properties (n_n, Option) ;
%--> find member length and direction cosines
dx = x(e_nodes(2)) - x(e_nodes(1)) ; % x length
dy = y(e_nodes(2)) - y(e_nodes(1)) ; % y length
L_e = sqrt (dx * dx + dy * dy) ; % total length
cx = dx / L_e ; cy = dy / L_e ; % direction cosines
% ELEMENT CONDUCTION AND INTERNAL SOURCE MATRICES
% for q = 1:n_q ; % Loop over quadrature points ----> ---->
% Linear axial bar and cubic bending. DOF = u, v, u, v
% Form arrays in local axes, transform. 1 2 3 4
% stiffness
s_L = [ 1, 0, -1, 0 ;
0, 0, 0, 0 ;
-1, 0, 1, 0 ;
0, 0, 0, 0 ] * E * A / L_e ;
% Map line load to node forces ONLY
if ( any (Line_e) ) % then form forcing vector
fprintf ('WARNING: line loads not yet active. ')
end % if for set up line load nodal resultants
% Optional local mass matrix
if ( Rho > 0 )
m_L = [140, 0, 70, 0 ;
0, 0, 0, 0 ;
70, 0, 140, 0 ;
0, 0, 0, 0 ] * Rho * A * L_e ;
end % if mass requested
% end % for loop over n_q element quadrature points <---- <----
% Define local to system DOF transformation matrix
t_L = [ cx cy 0 0 ; % (inverse t_L = transpose t_L)
-cy cx 0 0 ;
0 0 cx cy ; % joint 2
0 0 -cy cx ] ;
% Transform from local to system
S_e = t_L' * s_L * t_L ; M_e = t_L' * m_L * t_L ;
C_e = t_L' * c_L ;
% 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 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);
end % if essential BC exist (almost always true)
% ECHO PROPERTIES (add Rho)
fprintf ('Application properties are: ')
fprintf ('Elastity modulus = %g ', E)
fprintf ('Cross-sectional area = %g ', A)
fprintf ('Line Load = [ %g %g ] ', Line_e(1), Line_e(2))
% ENFORCE ESSENTIAL BOUNDARY CONDITIONS
save_resultant_load_vectors (n_g, C)
if ( all ( C == 0 ) ) % then null solution
fprintf ('WARNING: No loads or support movement applied. ')
end % if null solution
[S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C);
% COMPUTE SOLUTION & SAVE
T = S C ; % Compute displacements & rotations
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 ELEMENT REACTIONS (MEMBER FORCES)
% output_2D_truss_el_reactions (n_e, n_g, n_n, n_q, nodes, x, y, T)
% +++++++++++++ functions in alphabetical order +++++++++++++++++
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 p = 1:n_d % check all DOF for essential BC
if ( flag_EBC (p) ) % 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 (p) ; % recover EBC value
C (:) = C (:) - EBC * S (:, p) ; % carry known column to RHS
S (:, p) = 0 ; S (p, :) = 0 ; % clear, restore symmetry
S (p, p) = 1 ; C (p) = 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)
end
function [a, b, c, center, two_A] = form_T3_geom_constants (x, y, e_nodes)
% Planar 3 node triangle geometry: H_i (x,y) = (a_i + b_i*x + c_i*y)/two_a
% define nodal coordinates, ccw: i, j, k
x_e = x(e_nodes) ; y_e = y(e_nodes) ; % coord at el nodes
x_i = x_e(1) ; x_j = x_e(2) ; x_k = x_e(3) ; % change notation
y_i = y_e(1) ; y_j = y_e(2) ; y_k = y_e(3) ; % change notation
% define centroid coordinates (quadrature point)
center (1) = (x_i + x_j + x_k)/3 ;
center (2) = (y_i + y_j + y_k)/3 ;
% geometric parameters: H_i (x,y) = (a_i + b_i*x + c_i*y)/two_a
a_i = x_j * y_k - x_k * y_j ; b_i = y_j - y_k ; c_i = x_k - x_j ;
a_j = x_k * y_i - x_i * y_k ; b_j = y_k - y_i ; c_j = x_i - x_k ;
a_k = x_i * y_j - x_j * y_i ; b_k = y_i - y_j ; c_k = x_j - x_i ;
a (1:3) = [a_i, a_j, a_k] ;
b (1:3) = [b_i, b_j, b_k] ;
c (1:3) = [c_i, c_j, c_k] ;
% calculate twice element area
two_A = a_i + a_j + a_k ; % = b_j*c_k - b_k*c_j also
% end form_T3_geom_constants (x, y, e_nodes)
end
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 ('Read %g point sources. ', n_u)
fprintf ('Node DOF Source_value ')
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 ', 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
fprintf (' ')
% end get_and_add_point_sources (n_g, n_m, C)
end
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
end
function [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag)
EBC_count = sum( sum ( EBC_flag > 0 ) ) ; % # of EBC
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 essential boundary condition sets. ', n_c)
fprintf ('Read %g essential boundary condition sets. ', n_c)
if ( n_c ~= EBC_count ) % then probable user error
fprintf ('WARNING, expected %g EBC sets. ', EBC_count)
end % if error expected
fprintf ('Node DOF Value_EBC ')
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
fprintf ('%g %g %g ', node, DOF, value)
if ( EBC_flag (node, DOF) == 0 ) % check data consistency
fprintf ('WARNING: EBC but no flag at node %g & DOF %g. ', ...
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
fprintf (' ')
% end get_ebc_values
end
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
end
function [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements (pre_e) ;
% MODEL input file controls (for various data generators)
if (nargin == 0) ; % default to no proceeding items in data
pre_e = 0 ; % Dummy items before el_type & connectivity
end % if
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) - pre_e - 1 ; % nodes per element
el_type = round (msh_typ_nodes(:, pre_e+1)); % el type number >= 1
n_t = max(el_type) ; % number of element types
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, (pre_e+2:pre_e+1+n_n));
fprintf ('Read %g elements. ', n_e)
fprintf ('Elem, Type, Connectivity_List ')
for j = 1:n_e
if ( n_n == 1 )
fprintf ('%g %g %g ', j, el_type(j), nodes(j,:));
elseif ( n_n == 2 )
fprintf ('%g %g %g %g ', j, el_type(j), nodes(j,:));
elseif ( n_n == 3 )
fprintf ('%g %g %g %g %g ', j, el_type(j), nodes(j,:));
elseif ( n_n == 4 )
fprintf ('%g %g %g %g %g %g ', j, el_type(j), nodes(j,:));
elseif ( n_n == 5 )
fprintf ('%g %g %g %g %g %g %g ', j, el_type(j), nodes(j,:));
elseif ( n_n == 6 )
fprintf ('%g %g %g %g %g %g %g %g ', j, el_type(j), nodes(j,:));
else
fprintf ('%g %g ', j, el_type(j)); dips( nodes(j,:));
end % if
end % for each element
fprintf ('Maximum element type number = %g ', n_t)
% end get_mesh_elements
end
function [n_m, n_s, P, x, y, z] = get_mesh_nodes (pre_p) ;
% MODEL input file controls (for various data generators)
if (nargin == 0) % set usual default
pre_p = 0 ; % Dummy items before BC_flag % coordinates
end % if
% 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) - pre_p - 1 ; % number of space dimensions
msh_bc_xyz (:, (pre_p+1))= round (msh_bc_xyz (:, (pre_p+1)));
P = msh_bc_xyz (1:n_m, (pre_p+1)) ; % integer Packed BC flag
x = msh_bc_xyz (1:n_m, (pre_p+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, (pre_p+3)) ; % extract y column
end % if 2D or 3D
if ( n_s == 3 ) ; % check 3D
z = msh_bc_xyz (1:n_m, (pre_p+4)) ; % extract z column
end % if 3D
%b if ( pre_p ~= 1 ) % not given node number, sequential data
fprintf ('Read %g nodes. ', n_m)
fprintf ('Node, BC_Flag, Coordinates ')
for j = 1:n_m ; % list nodes
fprintf ('%g %g %g %g %g ', j, P(j), x(j), y(j), z(j)) ;
end % for j DOF
fprintf (' ')
% end get_mesh_nodes
end
function list_save_beam_displacements (n_g, n_m, T)
fprintf (' ') ;
fprintf('Node Y_displacement Z_rotation at %g nodes ', 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 ', T_matrix (j, 1:n_g)) ; % to file
fprintf ('%g %g %g ', j, T_matrix (j, 1:n_g)) ; % to screen
end % for j DOF
% end list_save_beam_displacements (n_g, n_m, T)
end
function list_save_displacements_results (n_g, n_m, T)
fprintf('Computed nodal displacements at %g nodes ', n_m)
fprintf('Node DOF_1 DOF_2 DOF_3 DOF_4 DOF_5 DOF_6 ')
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 x = 1:n_m ; % save displacements
if ( n_g == 1 )
fprintf (fid, '%g ', T_matrix (x, 1:n_g)) ;
fprintf ('%g %g ', j, T_matrix (x, 1:n_g)) ;
elseif ( n_g == 2 )
fprintf (fid, '%g %g ', T_matrix (x, 1:n_g)) ;
fprintf ('%g %g %g ', j, T_matrix (x, 1:n_g)) ;
elseif ( n_g == 3 )
fprintf (fid, '%g %g %g ', T_matrix (x, 1:n_g)) ;
fprintf ('%g %g %g %g ', j, T_matrix (x, 1:n_g)) ;
elseif ( n_g == 4 )
fprintf (fid, '%g %g %g %g ', T_matrix (x, 1:n_g)) ;
fprintf ('%g %g %g %g %g ', j, T_matrix (x, 1:n_g)) ;
elseif ( n_g == 5 )
fprintf (fid, '%g %g %g %g %g ', T_matrix (x, 1:n_g)) ;
fprintf ('%g %g %g %g %g %g ', j, T_matrix (x, 1:n_g)) ;
elseif ( n_g == 6 )
fprintf (fid, '%g %g %g %g %g %g ', T_matrix (x, 1:n_g)) ;
fprintf ('%g %g %g %g %g %g %g ', j, T_matrix (x, 1:n_g)) ;
else
error('reformat list_save_displacements_results for n_g > 6.')
end % if
end % for j DOF
fprintf (' ') ;
% end list_save_displacements_results (T)
end
function list_save_temperature_results (T)
n_m = size (T, 1) ; % get size
fprintf('Temperature at %g nodes ', 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 ', T (j)) ; % print
fprintf (' %g %g ', j, T (j)) ; % sequential save
end % for j DOF
% end list_save_temperature_results (T)
end
function output_PlaneStress_stresses(n_e, n_g, n_n, n_q, nodes, x,y,T)
% POST-PROCESS ELEMENT STRESS RECOVERY & SAVE
fid = fopen('el_qp_xyz_fluxes.tmp', 'w') ; % open for writing
fprintf (' ') ; % blank line
fprintf('Elem, QP, X_qp, Y_qp ') ;% header
fprintf('Elem, QP, Stress_qp: xx yy xy ');% 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);
[t_e, Body_e, E_e] = set_constant_plane_stress_prop; % properties
% 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:2:5) = b (1:3)/two_A ; B_e (2, 2:2:6) = c (1:3)/two_A ;
B_e (3, 1:2:5) = c (1:3)/two_A ; B_e (3, 2:2:6) = b (1:3)/two_A ;
% COMPUTE GRADIENT & HEAT FLUX, SAVE LOCATION AND VALUES
Strains = B_e * T_e ; % mechanical strain
Stress = E_e * Strains ; % mechanical stress
fprintf (fid,'%g %g %g %g %g ', center(1), center(2), ...
Stress(1), Stress(2), Stress(3));% save
fprintf ('%g %g %g %g ', j, q, center(1:2));% prt
fprintf ('%g %g %g %g %g ', j, q, Stress(1:3));% prt
fprintf (' ') ;% prt
end % for loop over n_q element quadrature points <----
end % for each j element in mesh
% end output_PlaneStress_stresses (n_e, n_g, n_n, n_q, nodes, x, y, T)
end
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 (' ') ; % blank line
fprintf('Elem, X_qp, Y_qp, HeatFlux_x, HeatFlux_y ');% header
end
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);
[t_e, Body_e, E_e] = set_constant_plane_stress_prop; % properties
% 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 ', center(1:2), HeatFlux(1:2));% save
fprintf ('%g %g %g %g %g ', 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 (n_e, n_g, n_n, n_q, nodes, x, y, T)
end
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 ('Recovered %g Reactions ', ...
sum (sum (EBC_flag > 0))) ; % header
fprintf ('Node, DOF, Value of reaction ')
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 ', node, j_g, React);% save
fprintf ('%g %g %g ', 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 ('%g DOF Totals = ', n_g) ; disp(Totals) ; % echo totals
fprintf (' ') ; % Skip a line
% end recover_reactions_print_save (EBC_row, EBC_col, T)
end
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)
end
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 (' ') ; % Skip a line
% fprintf ('Node, DOF, Resultant Force (1) or Moment (2) ')
fprintf ('Node, DOF, Resultant input sources ')
fid = fopen('node_resultant.tmp', 'w'); % open for writing
Totals = zeros (1, n_g) ; % zero input totals
for j = 1:n_d ; % extract resultants
if ( C (j) ~= 0. ) ; % then source here
% Output node_number, component_number, value
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 ', node,j_g,value,j);% save
fprintf ('%g %g %g ', 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 ('%g DOF Totals = ', n_g) ; disp(Totals) ; % echo totals
% end save_resultant_load_vectors (n_g, n_m, C)
end
function [I, E, Rho, Line_e]=set_constant_beam_prop (n_n, Option);
if ( nargin == 1 )
Option = 1 ;
elseif ( nargin == 0 )
n_n = 2 ; Option = 1 ;
end % if problem Option number
Line_e = zeros (n_n, 1) ; % default line load at nodes
switch Option
case 1 % Propped cantilever with uniform load, L, L/4
% Wall reactions: V=37*Line*L/64 M=7*Line*L^2/64
% Roller reaction: R= 43*Line*L/64
% Total vertical load: 5*Line*L/4
% *-----(1)-----*-----(2)-----*--(3)--* EI
% Fixed@1 L/2 2 Roller@3 L/4 4
I = 1.0 ; E = 1.0 ; Rho = 0.0 ; Line_e = [1.0; 1.0] ;
case 2 % cantilever with uniform load, L
I = 1.0 ; E = 1.0 ; Rho = 0.0 ; Line_e = [1.0; 1.0] ;
otherwise
I = 1.0 ; E = 1.0 ; Rho = 0.0; % default shape & material
end % switch
% end set_constant_beam_prop;
end
function [t_e, Body_e, E_e] = set_constant_plane_stress_prop ;
t_e = 1 ; Body_e (1:2) = 0. ; % defaults
% case 1
t_e = 5e-3 ; % thickness
Body_e (1:2) = [5e5, 0.] ; % components
E = 15e9 ; % Elastic modulus
nu = 0.25 ; % Poisson's ratio
% plane stress
E_v = E/(1 - nu^2) ; % constant
E_e (1, 1) = E_v ; E_e (1, 2) = E_v * nu ; % non-zero term
E_e (2, 1) = E_v * nu ; E_e (2, 2) = E_v ; % non-zero term
E_e (3, 3) = E_v * (1 - nu) / 2 ; % non-zero term
%end set_constant_plane_stress_prop
end
function [t_e, Q_e, E_e] = set_constant_2D_conduction_prop
% Manually set constant element properties (Fig 11.9 text)
q_e = 0. ; t_e = 1. ; % defaults
% case 1
Kx = 8. ; Ky = 8. ; Kxy = 0. ; % thermal conductivity
% case 2
kx = 1. ; Ky = 1. ; Kxy = 0. ;
% insert
E_e = zeros (2, 2) ; % constitutive matrix
E_e (1, 1) = Kx ; E_e (1, 2) = Kxy ; % non-zero term
E_e (2, 1) = Kxy ; E_e (2, 2) = Ky ; % non-zero term
% end set_constant_2D_conduction_prop
end
function [A, E, I, Line_e, Rho] = set_2D_frame_properties (n_n, Option)
if ( nargin == 1 )
Option = 1 ;
elseif ( nargin == 0 )
n_n = 2 ; Option = 1 ;
end % if problem Option number
Line_e = zeros (n_n, 1) ; % default line load at nodes
switch Option
case 1
% Weaver plane frame example X-------- F_y=-32 K,
% E=10000ksi, A=10 in sq 2 (1) 1 M_z=-1050 in-K,
% I=1000 in^4, L_1=100 in (2) at node 1.
% L_2x=100 in, L_2y=-75 in [no line load]
% Node 1 disp: -0.02026 in, -0.09936 in, X
% and -0.001798 radians 3
% Reactions, node 2: 20.261 K, 1.1378 K, 236.65 in K
% node 3:-20.261 K, 30.862 K, -639.52 in-K
A=10; I=1e3; E=1e4; Rho=0; Line_e = [0.0; 0.0] ;
case 2 % cantilever with uniform load, L
A = 1 ; I = 1 ; E = 1 ; Rho = 0 ; Line_e = [1.0; 1.0] ;
otherwise
A = 1 ; I = 1 ; E = 1 ; Rho = 0; % default shape & material
end % switch
% end set_2D_frame_properties (n_n, Option)
end
function [A, E, Line_e, Rho] = set_2D_truss_properties (n_n, Option)
if ( nargin == 1 )
Option = 1 ;
elseif ( nargin == 0 )
n_n = 2 ; Option = 1 ;
end % if problem Option number
Line_e = zeros (n_n, 1) ; % default line load at nodes
switch Option
case 1
% 2 4 5 Meek's Example 7.2 truss
% *--(10)-*--(11)-* E = 30,000 ksi, A = 1 in^2
% |(4) /|(7) /| Two 10 inch bays
% | / | / | Max vertical deflection @ 3
% (3) X (6) X (9) is -4.4013E-03 inches
% | / | / |
% |/(5) |/(8) | Reactions are 5K each at 1, 6
% 1 #--(1)--*--(2)--o 6
% Pin 3| Roller
% v P=10K
A=1; E=30e3 ; Rho=0; Line_e = [0.0; 0.0] ;
otherwise
A = 1 ; E = 1 ; Rho = 0; % default shape & material
end % switch
% end set_2D_truss_properties (n_n, Option)
end
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
%b size(flag)
%b size(full)
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)
end % if likely user error
% end unpack_pt_flags
end
Related Questions
drjack9650@gmail.com
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.