cosmos_cantilever_twoside.m

%% Script parameters: change as needed

% Copyright 2006, The MathWorks, Inc.

%
% Directory in which to find the data files
%
datadir='/local/user/vchudnov/Private/TMW/AEU/FlexPaper/data';
datadir='/home/vchudnov/Playground/SimMechanics/FlexibleBodies/Paper/MLNN2006Q1/AEU/FlexPaper/data';
if ispc
    datadir='\Playground\SimMechanics\FlexibleBodies\Paper\MLNN2006Q1\AEU\FlexPaper\data';
end

%
% Name prefix for the data files
%
namePrefix='cantilever_part';

%
% In this example, we use uniform damping (the same for all the modes).
% This is the coefficient that multiplies all the omegas in the damping
% matrix and is equal to 2*xi.
%
uniformDampingRatio = .15;

%
% Whether we should attempt to trim the mesh when obtaining the modes, so
% that the trimmed mesh only contains the input and output nodes. When the
% FEA mesh is very fine, this can make the MATLAB data structures smaller
% by throwing out data we don't use.
%
doTrimMesh = true;

%
% A scale factor for the force (e.g., to convert between systems of units).
% In this example, both the FEA model and the SimMechanics model use MKS
% units, so this is unity.
%
forcescale=1;

%
% Should we model the input as being distributed along the entire right
% edge of the beam or just at the first point we find on the right edge?
%
distributeinput=true;

%
% Which modes to include on the analysis
%
whichModes = [ 1 ];

%% Obtaining the FEA data

%
% Get the original mesh
% 
[ md , theMesh ] = cosmos2m( datadir , namePrefix , true );
origMesh = theMesh;
figure(1);
subplot(3,1,1);plot3( theMesh(:,2), theMesh(:,3), theMesh(:,4) , 'b.');
title('Original Mesh'); xlabel('x');ylabel('y');zlabel('z');
fprintf( 1 , '==>Original mesh has %i points\n' , length( theMesh) );

% Select the nodes of particular interest
rightEdgeXCoord = max( theMesh( : , 2 ) );
leftEdgeXCoord = min( theMesh( : , 2 ) );


%
% Trim the mesh if requested
%
rightMeshNodes = ( theMesh( : , 2 ) == rightEdgeXCoord ) ;
leftMeshNodes  = ( theMesh( : , 2 ) == leftEdgeXCoord ) ;
if doTrimMesh
    trimMeshNodes = (  rightMeshNodes | leftMeshNodes );
else
    trimMeshNodes = ones( [ 1 length( theMesh ) ] ); 
end

trimMeshIndeces = find( trimMeshNodes );
trimMesh = theMesh( trimMeshIndeces , : );

trimMeshNum = length( trimMeshIndeces );

 
%
% Process the FEA data at the selected nodes only
%
[ md , theMesh ] = cosmos2m( datadir , namePrefix , [] , trimMeshIndeces , whichModes );
fprintf( 1 , '==>New mesh (x==%f,%f) has %i points.\n' , rightEdgeXCoord , leftEdgeXCoord , trimMeshNum );
subplot(3,1,2);plot3( theMesh(:,2), theMesh(:,3), theMesh(:,4) , 'b.');
title('Trimmed Mesh');xlabel('x');ylabel('y');zlabel('z');


%% Process the FEA data for the state-space model

%
% Make sure to identify the nodes at the right edge (necessary if we did
% not trim the mesh earlier )
%

rightMeshNodes = ( theMesh( : , 2 ) == rightEdgeXCoord ) ;
leftMeshNodes  = ( theMesh( : , 2 ) == leftEdgeXCoord ) ;
interestMeshNodes = (  rightMeshNodes | leftMeshNodes );

interestMeshIndeces = find( interestMeshNodes );
interestEdgeMeshNum = sum( interestMeshNodes );
interestEdgeMesh    = theMesh( interestMeshIndeces , : );

rightMeshIndeces = find( rightMeshNodes );
rightEdgeMeshNum = sum( rightMeshNodes );
rightEdgeMesh    = theMesh( rightMeshIndeces , : );

leftMeshIndeces = find( leftMeshNodes );
leftEdgeMeshNum = sum( leftMeshNodes );
leftEdgeMesh    = theMesh( leftMeshIndeces , : );


subplot(3,1,3);plot3( interestEdgeMesh(:,2), interestEdgeMesh(:,3), interestEdgeMesh(:,4) , 'b.');
title('Interest Edge Mesh');xlabel('x');ylabel('y');zlabel('z');

%
% Set up the input nodes and their scaling magnitudes
%
systemInputs = zeros( [ 2 trimMeshNum ] );
if distributeinput
    rightEdgeDistributedForce = forcescale / rightEdgeMeshNum;
    leftEdgeDistributedForce = forcescale / leftEdgeMeshNum;
    systemInputs( 1 , rightMeshIndeces  ) = rightEdgeDistributedForce;
    systemInputs( 2 , leftMeshIndeces   ) = leftEdgeDistributedForce;
else
    rightEdgeDistributedForce = forcescale;
    leftEdgeDistributedForce = forcescale;
    systemInputs( 1 , rightMeshIndeces( 1 )  ) = rightEdgeDistributedForce;
    systemInputs( 2 , leftMeshIndeces( 1 )   ) = leftEdgeDistributedForce;
end



%
% Set up the output nodes
% 
whichOutput = [ rightMeshIndeces( 1 ) leftMeshIndeces( 1 ) ];
outputCoordinate = interestEdgeMesh( whichOutput , : );
systemOutputs = zeros( [ 2 trimMeshNum ] );
systemOutputs( 1 , whichOutput(1) ) = 1;
systemOutputs( 2 , whichOutput(2) ) = 1;

    
%
% Encapsulate the state-space data
%
ssData.modeList = 1:md.num_modes;
ssData.nodeList = 1:trimMeshNum;
ssData.modalMatrix = md.nodal_displacements_Y( ssData.nodeList , ssData.modeList );
ssData.eigenvalues = md.frequency_angular';
ssData.modalDamping = uniformDampingRatio * ones( size( ssData.modeList ) );
ssData.inputNodes = systemInputs';
ssData.outputNodes = systemOutputs;




%% Create the state-space system
[ ss_cantilever_twoside_pva , ss_cantilever_twoside_p ] = fea_to_statespace( ssData.modalMatrix , ...
                             ssData.eigenvalues , ssData.modalDamping , ...
                             ssData.inputNodes , ssData.outputNodes );


%
% Save the state-space system
%
save ss_cantilever_twoside_pva ss_cantilever_twoside_pva
disp( '==> Saved ss_cantilever_pva' )