femjr.m

matlab实例

%*****************************************************************************%
% Program FEMjr v2.0 - Finite Element Method framework                        %
% based on Dr. Ken Fyfe's FEMjr code                                          %
%*****************************************************************************%

disp ('*********************************************************************')
disp ('* FEMjr v 2.0 - Finite Element Method junior - An FEM framework     *')
disp ('*********************************************************************')

clear                   % removes all variables from workspace

% Main components of the data structure:
%
% NOTEs: The actual order of the fields in the structures are not necessarily the
%        same as shown below.
%        The sintax used below is not really valid for matlab (although it is very close)
% 
% 
% Geom.title      % title of the problem read from input data file
%     .dim        % problem dimension (1, 2 or 3)
%     .nelemType  % number of element types defined
%     .nelems     % number of elements in the model
%     .nprops     % number of material properties
%     .nnodes     % number of nodes
%     .ndofs      % number of degrees of freedom in the model
%     .ndofs_node % number of dofs per node (makes sense only if all nodes have the same num of dofs!)
%
%     .ETags.spring1d   % tags used to identify element types
%           .truss2d
%
%     .ElemType(i).nnodes        % num. nodes for elements of type i
%                 .ndofs_node    % num. dofs per node for elements of type i
%
%     .ElemData(i).elemID            % user ID for this element (equal i, right now)
%                 .elemType          % type for this element
%                 .propID            % user ID of the property of this element
%                 .connect[1:nnodes] % connectivities for this element
%
%     .Prop(i).propID         % user ID for this property (equal i, right now)
%             .propType       % type for this property. This defines what is stored in properties[]
%             .properties[j]  % property j = K, A, E, etc (problem/element dependent)
%
%     .Nodes(i).nodeID        % user ID for this node (equal i, right now)
%              .coord[1:dim]  % coordinates for node with index i
%
%
% BCs.nnodeConstraint   % number of node constraints
%    .npointLoads       % number of point loads
%    .ndistrLoads       % number of distributed loads
%    .NodeConstraint(i).nodeID    % node to apply this constraint
%                      .dof       % node dof to apply this constraint
%                      .value     % value of the constraint
%
%    .PointLoad(i).nodeID         % node to apply this load
%                 .dof            % node dof to apply this load
%                 .value          
%
%    .DistrLoad(i).elemID           % element to apply this load
%                 .side             % side of the element with load 
%                 .traction[1:dim]  % traction/flux vector to apply
%
%

[Geom.ETags, Geom.ElemType] = define_elements;  % define elements available to FEMjr
Geom.nelemType  = length( Geom.ElemType );      % number of element types defined
%
%Geom

% Could make this a global object.
GaussRule = define_Gauss_rules;

[Geom,BCs, error] = read_input_data(Geom);   % read input data file
if error ~= 0 
    error('FEMjr: Error reading data file. Ending execution');
end
%Geom
%BCs
%explore(error);

print_input_data(Geom, BCs);           % echo data file to screen

% figure(1)
solution_vector = [];
plot_structure(Geom,BCs,solution_vector)     % plot structure

% uncomment if you want to preview model before solving
%button = questdlg('Do you wish to continue?',...
% 'Verify Geometry!','Yes','No','Yes');
%if strcmp(button,'Yes')
%    disp('Creating file')
%elseif strcmp(button,'No')
%    error('FEMjr terminated by user')
%end

kglobal     = zeros(Geom.ndofs);            % initialize global stiffness matrix
load_vector = zeros(Geom.ndofs,1);          % and load vector

for iel = 1:Geom.nelems                % form element stiffness matrices
  switch ( Geom.ElemData(iel).elemType )
    case(Geom.ETags.spring1d)                            % spring1d
      [ke,fe] = spring1d_elem(Geom.ElemData(iel), Geom.Prop);
    case(Geom.ETags.truss2d)                             % truss2d
      [ke,fe] = truss2d_elem(Geom.ElemData(iel), Geom.Prop, Geom.Nodes);
    case(Geom.ETags.pBar2)
      p = 1;
      [ke,fe] = pBar_elem(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule)
    case(Geom.ETags.pBar3)
      p = 2;
      [ke,fe] = pBar_elem(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule)  
    case(Geom.ETags.pBar4)
      p = 3;
      [ke,fe] = pBar_elem(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule)  
    case(Geom.ETags.pBar5)
      p = 4;
      [ke,fe] = pBar_elem(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule)  
  otherwise
       error('FEMjr: Invalid element type for assembly');
  end                                 % then assemble to global stiffness matrix
  [kglobal, load_vector] = assemble(ke, fe, Geom.ElemData(iel), Geom.ElemType, kglobal, load_vector);
end              
fprintf('\n\nStiffness Matrix and Load Vector after assembly')
kglobal
load_vector

% assembly point loads
for load_num = 1:BCs.npointLoads                
    load_vector = loads(BCs.PointLoad(load_num), Geom.ndofs_node, load_vector);
end
fprintf('\nLoad Vector after assembly of point loads')
load_vector

% apply node constraints
for ic = 1:BCs.nnodeConstraint          
  [kglobal,load_vector] = ...
      constraints(BCs.NodeConstraint(ic), Geom.ndofs, Geom.ndofs_node, kglobal,load_vector);
end           

fprintf('\nConstrained Stiffness Matrix and Load Vector\n')
print_stiffness_load(kglobal,load_vector,Geom.ndofs) % print stiffness matrix/load vector

solution_vector = solve_equations(kglobal,load_vector); % solve system of equations

fprintf ('\n<<<<---- Solution Results ---->>>>\n\n')
fprintf ('SOLUTION VECTOR\n dof#      value\n')
for dof_num = 1:Geom.ndofs
  fprintf('%4i %14.5g\n', dof_num, solution_vector(dof_num))
end
energy = 0.5 * solution_vector' * load_vector;
fprintf('\nStrain Energy (assume homogeneous displacement BCs) = %f\n',energy);

fprintf ('ELEMENTAL POSTPROCESSING\n')
energy = 0;
for iel = 1:Geom.nelems      % postprocess elemental results
    switch ( Geom.ElemData(iel).elemType )        % provide some table labels for the results
    case (Geom.ETags.spring1d)        % spring1d
        energy = energy + ...
            spring1d_results(Geom.ElemData(iel), Geom.ElemType, Geom.Prop, solution_vector);
    case(Geom.ETags.truss2d)        % truss2d
        energy = energy + ...
            truss2d_results(Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, solution_vector);
    case(Geom.ETags.pBar2)
        p = 1;
        energy = energy + ...
            pBar_results(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule, solution_vector);
    case(Geom.ETags.pBar3)
        p = 2;
        energy = energy + ...
            pBar_results(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule, solution_vector);
    case(Geom.ETags.pBar4)
        p = 3;
        energy = energy + ...
            pBar_results(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule, solution_vector);
    case(Geom.ETags.pBar5)
        p = 4;
        energy = energy + ...
            pBar_results(p, Geom.ElemData(iel), Geom.ElemType, Geom.Prop, Geom.Nodes, GaussRule, solution_vector);
    otherwise
        error('FEMjr: Invalid element type for postprocessing'); 
    end
end
fprintf('\nStrain Energy = %f\n',energy);

plot_structure(Geom,BCs, solution_vector)     % plot structure and deflections

%