Commit 18b557c3 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Resolved merge conflicts.

parents 40d1781a 7354cd3e
......@@ -186,8 +186,8 @@ classdef FeAnalysisDynamicFD < mafe.FeAnalysis
% due to the appearance of xh_ii as complex conjugates,
% the homogeneous solution should consist of real values only
norm_imag = norm( imag(xh), 2 );
if norm_imag > 1.e-8
disp('Warning: Solution has imaginary component!')
if norm_imag > 1.e-6
disp(sprintf('Warning: Solution has imaginary component! (norm = %2.3e)', norm_imag));
end
% return the real part only
xh = real(xh);
......@@ -223,11 +223,26 @@ classdef FeAnalysisDynamicFD < mafe.FeAnalysis
% return the real part only
xp = real(xp);
end
%% put xh+xp to dof solution
function putSolToDof(self, time)
u = self.solHomogeneous(time,0) + self.solInhomogeneous(time,0);
v = self.solHomogeneous(time,1) + self.solInhomogeneous(time,1);
a = self.solHomogeneous(time,2) + self.solInhomogeneous(time,2);
%% put a specific solution component to dof
function putSolToDof(self, time, component)
if ~exist('component', 'var')
component = 'full';
end
% extract solution
switch(component)
case 'full'
u = self.solHomogeneous(time,0) + self.solInhomogeneous(time,0);
v = self.solHomogeneous(time,1) + self.solInhomogeneous(time,1);
a = self.solHomogeneous(time,2) + self.solInhomogeneous(time,2);
case 'homogeneous'
u = self.solHomogeneous(time,0);
v = self.solHomogeneous(time,1);
a = self.solHomogeneous(time,2);
case 'particular'
u = self.solInhomogeneous(time,0);
v = self.solInhomogeneous(time,1);
a = self.solInhomogeneous(time,2);
end
% --- POST-PROCESSING ---
% get local force vector at specific time
p = zeros(self.fep.ndofs,1);
......@@ -251,8 +266,8 @@ classdef FeAnalysisDynamicFD < mafe.FeAnalysis
if nargin < 2
component = 'real';
end
lambda = self.L(mode)
a = self.A(mode)
lambda = self.L(mode);
a = self.A(mode);
switch component
case 'imag'
u = imag( a * self.X(:,mode) );
......@@ -261,7 +276,7 @@ classdef FeAnalysisDynamicFD < mafe.FeAnalysis
otherwise
u = a * self.X(:,mode);
end
u = 1./norm(u,2) * u
u = 1./norm(u,2) * u;
% compute reaction forces
K = self.fep.assembleSystemStiffnessMatrix();
f = self.fep.assembleSystemForceVector();
......
......@@ -34,37 +34,20 @@ classdef FeProblem < handle
node.idx = self.ndofs + [1:ndof];
self.ndofs = self.ndofs + ndof;
end
%sprintf('feproblem ndofs = %d', self.ndofs)
% initialise state to zero (global dof values)
z = zeros(self.ndofs,1);
% initialise nodal vectors lhs and rhs
self.updateNodalDofs(z,z)
% initialise element internal quantities (forces, stresses)
self.postprocessElements(z)
end
%% assemble system K matrix
function A = assembleSystemStiffnessMatrix(self)
% create system stiffness matrix
est_entries = round( length(self.elems) * 12*12 * 1.10 );
num_entries = 0;
id1 = zeros(est_entries,1);
id2 = zeros(est_entries,1);
mat = zeros(est_entries,1);
% loop over all elements and collect element stiffness matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.stiffness();
% determine number of entries
e_size = length(e_idx)^2;
% add element contribution to system matrix
id1(1+num_entries:num_entries+e_size) = repmat( e_idx, 1, length(e_idx) );
id2(1+num_entries:num_entries+e_size) = reshape( repmat( e_idx, length(e_idx), 1 ), [], 1 );
mat(1+num_entries:num_entries+e_size) = reshape( e_mat', [], 1 );
% add up to total number of entries
num_entries = num_entries + e_size;
end
A = sparse( id1(1:num_entries), id2(1:num_entries), mat(1:num_entries) );
%
%A = self.assembleSystemStiffnessMatrixElastic() + self.assembleSystemStiffnessMatrixGeometric();
function A = assembleSystemStiffnessMatrix(self)
A = self.assembleSystemStiffnessMatrixElastic() + self.assembleSystemStiffnessMatrixGeometric();
end
%% assemble system K matrix: elastic
function A = assembleSystemStiffnessMatrixElastic(self)
% create system stiffness matrix
est_entries = round( length(self.elems) * 12*12 * 1.10 );
num_entries = 0;
......
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 1 (task 2.1-2.2):
%
% definition of cross section and material properties
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
% show the objects created above and the variables they hold
s1
s2
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 2 (task 2.3):
%
% definition of nodes and beam/point elements
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 3 (task 2.4):
%
% definition of constraints (kinematic and static/dynamic)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 1000. ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 4 (task 3.1):
%
% create finite element problem (and plot the undeformed configuration)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 1000. ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
% visualise the system
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference');
% print the current state of the system
fep.printSystem();
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 5 (task 3.2-3.6):
%
% create and perform a static finite element analysis (and plot result)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 1000. ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: static
ana = mafe.FeAnalysisStatic(fep);
%% Calculate response
ana.analyse();
%% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference');
fep.plotSystem('deformed');
fep.printSystem();
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 6 (task 4.1-4.5):
%
% create and perform a dynamic finite element analysis (and extract results)
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
s1.E = 2.0e11; % [N/m^2]
s1.rho = 8.0e+3; % [kg/m^3]
s2 = mafe.Section0D(); % for the point mass (m, but no k)
s2.m = 50.0; % [kg]
s2.k = 0.00; % [N]
s3 = mafe.Section0D(); % for the point spring (k, but no m)
s3.m = 0.00; % [kg]
s3.k = 1000; % [N]
%% Definition of nodes in the system
% Node 1
n1 = mafe.Node2D( [0.0, 0.0] );
% Node 2
n2 = mafe.Node2D( [0.0, 2.5] );
% Node 3
n3 = mafe.Node2D( [2.0, 2.5] );
% Node 4
n4 = mafe.Node2D( [4.0, 2.5] );
% Vector of nodes
nodes = [ n1, n2, n3, n4 ];
%% Definition of elements in the system
% Element 1
e1 = mafe.Member2D( [n1, n2], s1 );
% Element 2
e2 = mafe.Member2D( [n2, n3], s1 );
% Element 3
e3 = mafe.Member2D( [n3, n4], s1 );
% Element 4
e4 = mafe.Point( n3, s2, mafe.DofType.Disp1 );
% Element 5
e5 = mafe.Point( n3, s2, mafe.DofType.Disp2 );
% Element 6
e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 1000. ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% Select an analysis type: dynamic
ana = mafe.FeAnalysisDynamicFD( fep );
%% Calculate response
ana.analyse();
%% Visualise system response (individual modes of vibration = eigenvectors)
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference'); % plot the undeformed system
ana.putModeToDof(1, 'real'); % select the first mode (complex-conjugate eigenstates)
fep.plotSystem('deformed'); % plot the mode
ana.putModeToDof(3, 'real'); % select the third mode (complex-conjugate eigenstates)
fep.plotSystem('deformed'); % plot the mode
ana.putModeToDof(5, 'real'); % select the fifth mode (complex-conjugate eigenstates)
fep.plotSystem('deformed'); % plot the mode
%% Extract results for eigenvalues and eigenvectors
eigenvalues = ana.L([1 3 5])
eigenvectors = ana.X(:,[1 3 5])
% get the first three (distinct) eigenfrequencies
angular_frequency_in_rad_per_second = abs(imag(eigenvalues)) % omega
natural_frequency_in_hertz = abs(imag(eigenvalues)) / (2*pi) % f
%% Tutorial 1: Structural Dynamics
% A frame system with point mass and spring - static and dynamic analysis
%
% step 7 (task 5.1-5.3):
%
% provide initial conditions and plot system response at certain time instances
%
addpath([ '..' filesep '..' filesep '..' filesep ]);
%% Definition of sections
% -----------------------------------------------------------DEFINE FE ITEMS---
%% Definition of sections (containing geometrical data and material data)
s1 = mafe.Section1D(); % for the members
s1.A = 5.0e-4; % [m^2] = 0.01m * 0.05m
s1.Iy = 1.0e-7; % [m^4] = 0.01m * (0.05m)^3 / 12
......@@ -40,64 +46,33 @@ e6 = mafe.Point( n4, s3, mafe.DofType.Disp1 );
% vector of elements
elems = [ e1, e2, e3, e4, e5, e6 ];
%% Definition of constraints in the system
% Constraint 1 and 2 (fixed vertically and horizontally at base)
% Constraint 1 and 2 (fixed vertically and horizontally at lower left base point)
c1 = mafe.ConstraintNode( n1, mafe.DofType.Disp1, mafe.ConstraintType.Dirichlet, 0.0 );
c2 = mafe.ConstraintNode( n1, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 3 (fixed vertically right end)
% Constraint 3 (fixed vertically upper right end)
c3 = mafe.ConstraintNode( n4, mafe.DofType.Disp2, mafe.ConstraintType.Dirichlet, 0.0 );
% Constraint 4 (single load)
% Constraint 4 (horizontal point force at the corner)
c4 = mafe.ConstraintNode( n2, mafe.DofType.Disp1, mafe.ConstraintType.Neumann, 1000. ); % 1000 [N]
% vector of constraints
const = [ c1, c2, c3, c4 ];
% --------------------------------------------DEFINE FE PROBLEM AND ANALYSIS---
%% Setup of the finite element problem
fep = mafe.FeProblem( nodes, elems, const );
%% -----------------------------------------------------------------------------
%% Select an analysis type: static
ana = mafe.FeAnalysisStatic(fep);
%% Calculate response
ana.analyse();
%% Visualise system response
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
fep.plotSystem('reference');
fep.plotSystem('deformed');
fep.printSystem();
%% -----------------------------------------------------------------------------
disp('press key to proceed'); pause;
%% -----------------------------------------------------------------------------
%% Select a n analysis type: dynamic
%% Select an analysis type: dynamic
ana = mafe.FeAnalysisDynamicFD( fep );
%% Calculate response
ana.analyse();
%% Visualise system response (individual modes of vibration = eigenvectors)
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
ylim([-1 3])
ana.putModeToDof(1, 'real');
fep.plotSystem('reference');
fep.plotSystem('deformed');
%% -----------------------------------------------------------------------------
disp('press key to proceed'); pause;
%% -----------------------------------------------------------------------------
%% Visualise system response (homogeneous solution for given initial conditions)
%% Visualise system response (overall solution for given initial conditions)
cla; clf;
fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
% apply initial conditions
modeshape = real(ana.X(:,1:2:fep.ndofs)); % get some selected (or all) modes
u0 = 1.00 * modeshape(:,1) + 0.00 * modeshape(:,2);
u0 = 0.00 * modeshape(:,1) + 1.00 * modeshape(:,2);
v0 = 0.00 * modeshape(:,1) + 0.00 * modeshape(:,2);
ana.applyInitialConditions( u0, v0 );
% evaluate and plot for given time instants
% first undamped eigenfrequency is oemga = ~5.9 [rad/s] -> T = 2pi/omega = 1.1 [s]
for tt = 0:0.02:1.1*2
% first undamped eigenfrequency is omega = ~5.9 [rad/s] -> T = 2pi/omega = 1.1 [s]
for tt = 0:0.02:1.1*5
ana.putSolToDof(tt);
clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
xlim([-1 5])
......@@ -106,3 +81,18 @@ for tt = 0:0.02:1.1*2
fep.plotSystem('deformed');
drawnow;
end
% % The same as above, but with a video being generated (*.avi format)
% video = VideoWriter( mfilename('fullpath') );
% open(video);
% for tt = 0:0.02:1.1*5
% ana.putSolToDof(tt);
% clf; subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
% xlim([-1 5])
% ylim([-1 3])
% fep.plotSystem('reference');
% fep.plotSystem('deformed');
% writeVideo(video, getframe);
% drawnow;