pbar_results.m

matlab实例

%
function energy = pBar_results(p, ElemData, ElemType, Prop, Nodes, GaussRule, solution_vector)
% compute and print results for pBar element

% determine elemental displacement values from solution vector
elemType   = ElemData.elemType;
nnodes     = ElemType(elemType).nnodes;
ndofs_node = ElemType(elemType).ndofs_node;

for node = 1:nnodes
    dof_pos = (ElemData.connect(node) - 1) * ndofs_node + 1; % assuming all nodes with same num dofs
    disp_vector(node) = solution_vector(dof_pos);
end

% determine stiffness matrix
[ke, fe] = pBar_elem(p, ElemData, ElemType, Prop, Nodes, GaussRule);

% forces applied to the element
force = ke * disp_vector';                   % form matrix-vector product

% strain energy
energy = 0.5 * disp_vector * force;

% Compute stress at end nodes (nodes 1 and nnodes)
% Adjust the signs: Assume that coord. node 1 < coord. node nnodes
%         + force at node 1     : bar under compres., stress < 0
%         + force at node nnodes: bar under tension., stress > 0
% TBD: Use the same procedure as in truss2d_results to find the sign of
% stress
node = 1; % left node
stress(1) = - force(node)/ Prop( ElemData.propID ).properties(2);
node = nnodes; % right node
stress(2) =   force(node)/ Prop( ElemData.propID ).properties(2);

% Note: The axial force is discontinuous at a node if we have a point force applied there.

fprintf(['  ID         f_i                 Sxx(end nodes)          energy\n'])
fprintf('%4i ', ElemData.elemID)
fprintf('%8.3g ', force),  fprintf('   ')
fprintf('%8.3g', stress), fprintf('   ')
fprintf('%8.3g', energy)
fprintf('\n')

% Another way of computing stress/strain at nodes
% Here, we avoid computing the element stiffness matrix
% We also get the correct sign, regardless of the node numbering used (see
% above)
switch(p)
case (1)
    ksi_nodes = [-1.0, 1.0];
case (2)
    ksi_nodes = [-1.0, 0.0, 1.0];
case (3)
    ksi_nodes = [-1.0, -1/3, 1/3, 1.0];
case (4)
    ksi_nodes = [-1.0, -1/2, 0.0, 1/2, 1.0]; 
otherwise
    error('pBar_results: invalid value for p-order');
end

for node = 1:nnodes
    [psi, dpsi] = pBar_shape_lag( ksi_nodes(node), p);
         
    % Jacobian and global coordinates of this point
    [J, x_pt] = pBar_jacobian( ksi_nodes(node), p, ElemData, Nodes);

    % global derivatives of shape functions, i.e., derivatives w.r.t. x
    dpsi = (1./J) * dpsi;

    % strain at this node =  du/dx = u_i * dN_i/dx
    strain(node) = disp_vector * dpsi';
end    
stress_nd = Prop( ElemData.propID ).properties(1) * strain % Sxx = E * strain