Skip to content

GitLab

Commit 03febfb1 authored by Andreas Zilian's avatar Andreas Zilian
Browse files

Initial commit.

parents
Hide whitespace changes
Inline Side-by-side
+mafe/Beam1D.m 0 → 100644
View file @ 03febfb1
classdef Beam1D < mafe.Li2Element1D
% beam finite element in axial direction only (1D)
properties
%
end
methods
%% constructor
function obj = Beam1D(varargin)
obj = obj@mafe.Li2Element1D(varargin{:});
end
% ======================================================================
%% activation of dofs at the nodes of the element
function activateDofTypes(self)
ineed = [ mafe.DofType.Disp2 mafe.DofType.Rota3 ];
for i = 1:length(self.node)
self.node(i).dof = [self.node(i).dof ineed];
end
end
%% provide element index vector of dof
function idx = getDofSysIndices(self)
disp2 = find( [self.node.dof] == mafe.DofType.Disp2 );
rota3 = find( [self.node.dof] == mafe.DofType.Rota3 );
ndidx = [ self.node.idx ];
idx = [ ndidx(disp2(1)), ndidx(rota3(1)), ndidx(disp2(2)), ndidx(rota3(2)) ];
end
% ======================================================================
%% calculate transformation matrix local->global
function [L, T] = transform(self)
L = 2.*det( self.J(0.0) );
T = eye(4); % not needed in fact
end
% ======================================================================
%% calculation of element stiffness matrix
function ele_mat = stiffness(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
EI = self.sect.E * self.sect.Iy;
% element matrix
h1 = 2*EI/L;
h2 = 6*EI/L^2;
h3 = 12*EI/L^3;
ele_mat = [ +h3, +h2, -h3, +h2 ; ...
+h2, 2*h1, -h2, +h1 ; ...
-h3, -h2, +h3, -h2 ; ...
+h2, +h1, -h2, 2*h1 ];
end
%% calculation of element mass matrix
function ele_mat = mass(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
rhoA = self.sect.rho * self.sect.A;
% local element matrix
h1 = rhoA/420 * L^3;
h2 = rhoA/420 * L^2;
h3 = rhoA/420 * L;
ele_mat = [ 156*h3, 22*h2, 54*h3, -13*h2 ; ...
22*h2, 4*h1, 13*h2, -3*h1 ; ...
54*h3, 13*h2, 156*h3, -22*h2 ; ...
-13*h2, -3*h1, -22*h2, 4*h1 ];
end
%% calculation of element damping matrix
function ele_mat = damping(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
d_w = self.sect.d_w;
% local element matrix
h1 = d_w/420 * L^3;
h2 = d_w/420 * L^2;
h3 = d_w/420 * L;
ele_mat = [ 156*h3, 22*h2, 54*h3, -13*h2 ; ...
22*h2, 4*h1, 13*h2, -3*h1 ; ...
54*h3, 13*h2, 156*h3, -22*h2 ; ...
-13*h2, -3*h1, -22*h2, 4*h1 ];
% provide damping matrix based on Rayleigh assumption
% get parameters
d_alpha = self.sect.d_alpha;
d_beta = self.sect.d_beta;
% D = alpha * M + beta * K
ele_mat = ele_mat + d_alpha * self.mass() + d_beta * self.stiffness();
end
%% calculation of element force vector
function ele_vec = force(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
pX([1, 3]) = self.load.p(mafe.DofType.Disp2,:);
pX([2, 4]) = self.load.p(mafe.DofType.Rota3,:);
% local element vector
ele_vec = L/6./210 * [ 468, -630/L, 162, -630/L; ...
L*66, 126, L*39, -126; ...
162, 630/L, 468, 630/L; ...
-L*39, -126, -L*66, 126 ] * pX';
end
%------------------------------------------------------------------
%% graphical output: plot element
function plot(self, config, scale)
% my nodes
n1 = self.node(1);
n2 = self.node(2);
% coordinates and displacements
x1 = [ n1.ref(1), n2.ref(1) ];
uu = [ n1.lhs( find(n1.dof == mafe.DofType.Disp2) ), ...
n1.lhs( find(n1.dof == mafe.DofType.Rota3) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Disp2) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Rota3) ) ];
%
function H = interpolate_hermite(z)
[L, T] = self.transform();
H = [ 0.25*(1.-z).^2 .*( 2.+z), ...
0.25*(1.-z).^2 .*( 1.+z) * L/2., ...
0.25*(1.+z).^2 .*( 2.-z), ...
0.25*(1.+z).^2 .*(-1.+z) * L/2. ];
end
NN = 20;
xx = linspace(x1(1),x1(2),NN);
zz = linspace(-1,+1,NN);
ww = interpolate_hermite(zz') * uu';
%
% defaults
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
lwidth = 0.5;
mwidth = 10.;
ecolor = [0.4 0.4 0.4];
marker = '.';
%
switch config
case 'deformed'
XX1 = xx;
XX2 = ww*scale;
lcolor = [43.9, 50.6, 0.0]/100;
mcolor = [43.9, 50.6, 0.0]/100;
falpha = 0.6;
lwidth = 2.0;
mwidth = 20.;
otherwise
plot@mafe.Li2Element1D(self, config, scale);
return;
end
mafe.patchline(XX1, XX2, ...
'EdgeColor', lcolor, ...
'LineWidth', lwidth, ...
'EdgeAlpha', falpha);
plot(XX1([1 end]), XX2([1 end]), ...
marker, ...
'color', mcolor, ...
'markersize', mwidth);
end
%% print element data
function str = print(self)
str = sprintf('%+4.3e ', [ +self.internal(1), ... % left Q
-self.internal(2), ... % left M
-self.internal(3), ... % right Q
+self.internal(4) ]); % right M
end
end
end
+mafe/Beam2D.m 0 → 100644
View file @ 03febfb1
classdef Beam2D < mafe.Li2Element2D
% beam finite element in the plane (2D)
properties
%
end
methods
%% constructor
function obj = Beam2D(varargin)
obj = obj@mafe.Li2Element2D(varargin{:});
end
%% activation of dofs at the nodes of the element
function activateDofTypes(self)
ineed = [ mafe.DofType.Disp1 mafe.DofType.Disp2 mafe.DofType.Rota3 ];
for i = 1:length(self.node)
self.node(i).dof = [self.node(i).dof ineed];
end
end
%% provide element index vector of dof
function idx = getDofSysIndices(self)
disp1 = find( [self.node.dof] == mafe.DofType.Disp1 );
disp2 = find( [self.node.dof] == mafe.DofType.Disp2 );
rota3 = find( [self.node.dof] == mafe.DofType.Rota3 );
ndidx = [ self.node.idx ];
idx = [ ndidx(disp1(1)), ndidx(disp2(1)), ndidx(rota3(1)), ndidx(disp1(2)), ndidx(disp2(2)), ndidx(rota3(2)) ];
end
%% calculate transformation matrix local->global
function [L, T] = transform(self)
dx1 = self.node(2).ref(1) - self.node(1).ref(1);
dx2 = self.node(2).ref(2) - self.node(1).ref(2);
L = sqrt( dx1^2 + dx2^2 );
T = 1/L*[ dx1, dx2, 0, 0, 0, 0; ...
-dx2, dx1, 0, 0, 0, 0; ...
0, 0, L, 0, 0, 0; ...
0, 0, 0, dx1, dx2, 0; ...
0, 0, 0,-dx2, dx1, 0; ...
0, 0, 0, 0, 0, L ];
end
%% calculation of element stiffness matrix
function ele_mat = stiffness(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
EI = self.sect.E * self.sect.Iy;
% local element matrix
h1 = 2*EI/L;
h2 = 6*EI/L^2;
h3 = 12*EI/L^3;
ele_mat_local = [ 0, 0, 0, 0, 0, 0 ;
0, +h3, +h2, 0, -h3, +h2 ;
0, +h2, 2*h1, 0, -h2, +h1 ;
0, 0, 0, 0, 0, 0 ;
0, -h3, -h2, 0, +h3, -h2 ;
0, +h2, +h1, 0, -h2, 2*h1 ];
% global element matrix
ele_mat = T' * ele_mat_local * T;
end
%% calculation of element mass matrix
function ele_mat = mass(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
rhoA = self.sect.rho * self.sect.A;
% local element matrix
h1 = rhoA/420 * L^3;
h2 = rhoA/420 * L^2;
h3 = rhoA/420 * L;
ele_mat_local = [ 0, 0, 0, 0, 0, 0 ; ...
0, 156*h3, 22*h2, 0, 54*h3, -13*h2 ; ...
0, 22*h2, 4*h1, 0, 13*h2, -3*h1 ; ...
0, 0, 0, 0, 0, 0 ; ...
0, 54*h3, 13*h2, 0, 156*h3, -22*h2 ; ...
0, -13*h2, -3*h1, 0, -22*h2, 4*h1 ];
% global element matrix
ele_mat = T' * ele_mat_local * T;
end
%% calculation of element damping matrix
function ele_mat = damping(self)
% provide damping matrix based on Rayleigh assumption
% get parameters
d_alpha = self.sect.d_alpha;
d_beta = self.sect.d_beta;
% D = alpha * M + beta * K
ele_mat = d_alpha * self.mass() + d_beta * self.stiffness();
end
%% calculation of element force vector
function ele_vec = force(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
px1 = self.load.p(mafe.DofType.Disp1,:);
px2 = self.load.p(mafe.DofType.Disp2,:);
pxr = self.load.p(mafe.DofType.Rota3,:);
% local element vector
pX([1,4]) = T(1,1)*px1 + T(1,2)*px2;
pX([2,5]) = T(2,1)*px1 + T(2,2)*px2;
pX([3,6]) = pxr;
% local element vector
ele_vec_local = L/6./210 * [ 420, 0, 0, 210, 0, 0; ...
0, 468, -630/L, 0, 162, -630/L; ...
0, L*66, 126, 0, L*39, -126; ...
210, 0, 0, 420, 0, 0; ...
0, 162, 630/L, 0, 468, 630/L; ...
0, -L*39, -126, 0, -L*66, +126 ] * pX';
% global element vector
ele_vec = T' * ele_vec_local;
end
%------------------------------------------------------------------
%% graphical output: plot element
function plot(self, config, scale)
% my nodes
n1 = self.node(1);
n2 = self.node(2);
% coordinates and displacements
[L, T] = self.transform();
x1 = [ n1.ref(1), n2.ref(1) ];
x2 = [ n1.ref(2), n2.ref(2) ];
uu = [ n1.lhs( find(n1.dof == mafe.DofType.Disp1) ), ...
n1.lhs( find(n1.dof == mafe.DofType.Disp2) ), ...
n1.lhs( find(n1.dof == mafe.DofType.Rota3) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Disp1) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Disp2) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Rota3) )];
%
function N = interpolate_hermite(z)
[L, T] = self.transform();
N = [ (0.25*(1.-z).^2) .*(2.+z), ...
(0.25*(1.-z).^2) .*( 1.+z) * L/2., ...
(0.25*(1.+z).^2) .*(2.-z), ...
(0.25*(1.+z).^2) .*(-1.+z) * L/2. ];
end
NN = 20;
xx1 = linspace(x1(1),x1(2),NN)';
xx2 = linspace(x2(1),x2(2),NN)';
zz = linspace(-1,+1,NN);
ww = interpolate_hermite(zz') * uu([2,3,5,6])';
ww1 = T(2,1) * ww;
ww2 = T(2,2) * ww;
%
% defaults
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
lwidth = 0.5;
mwidth = 10.;
ecolor = [0.4 0.4 0.4];
marker = '.';
%
switch config
case 'deformed'
XX1 = xx1+ww1*scale;
XX2 = xx2+ww2*scale;
lcolor = [43.9, 50.6, 0.0]/100;
mcolor = [43.9, 50.6, 0.0]/100;
falpha = 0.6;
lwidth = 2.0;
mwidth = 20.;
otherwise
plot@mafe.Li2Element2D(self, config, scale);
return;
end
mafe.patchline(XX1, XX2, ...
'EdgeColor', lcolor, ...
'LineWidth', lwidth, ...
'EdgeAlpha', falpha);
plot(XX1([1 end]), XX2([1 end]), ...
marker, ...
'color', mcolor, ...
'markersize', mwidth);
end
%% print element data
function str = print(self)
str = sprintf('%+4.3e ', [ +self.internal(2), ... % left Q
-self.internal(3), ... % left M
-self.internal(5), ... % right Q
+self.internal(6) ]); % right M
end
end
end
+mafe/Constraint.m 0 → 100644
View file @ 03febfb1
classdef Constraint < matlab.mixin.Heterogeneous & handle
% Constraint on degree of freedom at node(s)
properties
cdof; % the dof type to constrain
type; % the constraint type
tfun = mafe.TimeFunction.Static; % default time function of constraint value
end
methods (Abstract)
getDofSysIndices(obj)
value(obj)
end
end
+mafe/ConstraintEdge.m 0 → 100644
View file @ 03febfb1
classdef ConstraintEdge < mafe.Constraint
% Constraint on degree of freedom along a polygonal edge composed of many nodes
properties
node ; % the nodes in correct order (polygon)
cval ; % the given contraint value(s), can be scalar or vectorial
% (assumed to act in the local edge coordinate system, not in the global one)
end
methods
%% constructor
function obj = ConstraintEdge(node, cdof, type, cval, tfun)
if nargin >= 4
obj.node = node;
obj.cdof = cdof;
obj.type = type;
obj.cval = cval;
end
if nargin >= 5
obj.tfun = tfun;
end
% consistency checks and convenience
if length(node) < 2
disp('ConstraintEdge: Error with given number of nodes < 2 !');
end
if length(cval) == 1 && length(type) == 1 % single scalar given
obj.cval = cval * ones(1,length(node));
end
if length(cval) == 2 && length(obj.cdof) == 2 % single vector given
obj.cval = [ cval(1) * ones(1,length(obj.node)), cval(2) * ones(1,length(obj.node)) ];
end
if length(obj.cval) ~= length(obj.node)*length(obj.cdof)
disp('ConstraintEdge: Error with given number of values != num nodes * num doftypes!')
end
end
% ----------------------------------------------------------------------
%% provide constraint index vector of dof
function idx = getDofSysIndices(self)
idx = [];
% loop over dof types
for cdof = self.cdof
% get nodal index positions
dof = find( [self.node.dof] == cdof );
% check
if length(dof) ~= length(self.node)
disp('ConstraintEdge:getDofSysIndices() --> Node missing requested dof!');
end
% get global dof indices
ndi = [ self.node.idx ];
idx = [ idx ndi(dof) ];
end
% we have: idx = [ idx doftype1 for all nodes, idx dotype2 for all nodes ]
end
%% provide constraint values
function vec = value(self)
% IMPORTANT NOTICE:
% All ConstraintEdge values are assumed to act in the **local**
% edge coordinate system, not in the global one
% This applies to both, Dirichlet and Neumann constraint types!
%
vec = zeros( 1, length(self.node)*length(self.cdof) );
% loop over all edges = 2-node pairs
for ii = 1:length(self.node)-1
% get node pair
n1 = self.node(ii );
n2 = self.node(ii+1);
% transformation matrix (note order of dof !)
dx1 = n2.ref(1) - n1.ref(1);
dx2 = n2.ref(2) - n1.ref(2);
L = sqrt( dx1^2 + dx2^2 );
T = 1/L*[ dx1, 0, dx2, 0; ...
0, dx1, 0, dx2; ...
-dx2, 0, dx1, 0; ...
0,-dx2, 0, dx1 ];
% convenvience: handle case of consecutive identical nodes
if ( L < eps )
continue;
end
% get parameters (in local coordinate system)
pT = self.cval( [ii, ii+1] );
pN = self.cval( [ii, ii+1]+length(self.node) );
% local element vector
vec_local = L/6.*[ +2, +1, 0, 0; ...
+1, +2, 0, 0; ...
0, 0, +2, +1; ...
0, 0, +1, +2 ] * [pT, pN]';
% global vector
vec_global = T' * vec_local;
%
idx = [ii, ii+1, ii+length(self.node), ii+1+length(self.node)];
vec( idx ) = vec( idx ) + vec_global';
end
end
end
end
+mafe/ConstraintNode.m 0 → 100644
View file @ 03febfb1
classdef ConstraintNode < mafe.Constraint
% Constraint on degree of freedom at single or multiple node(s)
properties
node ; % the involved node(s)
cval ; % the given contraint value(s), can be scalar or vectorial
end
methods
%% constructor
function obj = ConstraintNode(node, cdof, type, cval, tfun)
if nargin >= 4
obj.node = node;
obj.cdof = cdof;
obj.type = type;
obj.cval = cval;
end
if nargin >= 5
obj.tfun = tfun;
end
% consistency checks and convenience
if length(cval) == 1 && length(type) == 1 % single scalar given
obj.cval = cval * ones(1,length(node));
end
if length(cval) == 2 && length(obj.cdof) == 2 % single vector given
obj.cval = [ cval(1) * ones(1,length(obj.node)), cval(2) * ones(1,length(obj.node)) ]
end
if length(obj.cval) ~= length(obj.node)*length(obj.cdof)
disp('ConstraintNode: Error with given number of values != num nodes * num doftypes!')
end
end
% ----------------------------------------------------------------------
%% provide constraint index vector of dof
function idx = getDofSysIndices(self)
idx = [];
% loop over dof types
for cdof = self.cdof
% get nodal index positions
dof = find( [self.node.dof] == cdof );
% check
if length(dof) ~= length(self.node)
disp('ConstraintEdge:getDofSysIndices() --> Node missing requested dof!');
end
% get global dof indices
ndi = [ self.node.idx ];
idx = [ idx ndi(dof) ];
end
% we have: idx = [ idx doftype1 for all nodes, idx doftype2 for all nodes ]
end
%% provide constraint values
function vec = value(self)
vec = [ self.cval ];
end
end
end
+mafe/ConstraintType.m 0 → 100644
View file @ 03febfb1
classdef ConstraintType < uint8
% Enumeration of available constraint types
enumeration
Dirichlet (1), % constraint of Dirichlet type, e.g. displacement
Neumann (2), % constraint of Neumann type, e.g. force
end
end
+mafe/DofType.m 0 → 100644
View file @ 03febfb1
classdef DofType < uint8
% Enumeration of available dof types
enumeration
Disp1 ( 1), % displacement in y1
Disp2 ( 2), % displacement in y2
Disp3 ( 3), % displacement in y3
Rota1 ( 4), % rotation about y1
Rota2 ( 5), % rotation about y2
Rota3 ( 6), % rotation about y3
%
Velo1 ( 7), % velocity in y1
Velo2 ( 8), % velocity in y2
Velo3 ( 9), % velocity in y3
Omeg1 (10), % rotation velocity about y1
Omeg2 (11), % rotation velocity about y2
Omeg3 (12), % rotation velocity about y3
%
Press (13), % pressure
%
BoundaryStress1 (14), % boundary stress in y1 = Lagrange multiplier
BoundaryStress2 (15), % boundary stress in y2
BoundaryStress3 (16) % boundary stress in y3
end
end
+mafe/EleLoad.m 0 → 100644
View file @ 03febfb1
classdef EleLoad < handle
% An element load descriptor
properties
% nodal values of external distributed load
p = zeros(6,2); % stored in a matrix, rows = doftypes, cols = nodal vals
% time function
tfun = mafe.TimeFunction.Static;
end
methods
function obj = EleLoad(tdof, p, tfun)
if nargin >= 2
if length(p) == 1
obj.p(tdof,:) = [p, p];
else
obj.p(tdof,:) = p(1:2);
end
end
if nargin >= 3
obj.tfun = tfun;
end
end
end
end
+mafe/Element.m 0 → 100644
View file @ 03febfb1
classdef Element < matlab.mixin.Heterogeneous & handle
% A generic finite element
properties
internal = [];
end
methods
%------------------------------------------------------------------
%% post-process (compute internal forces)
function postprocess(self, u, du, ddu)
if ~exist('ddu', 'var')
ddu = zeros(length(u),1);
end
if ~exist('du', 'var')
du = zeros(length(u),1);
end
% transformation matrix
[L, T] = self.transform();
% perform post-processing for local element forces
f = self.stiffness() * u + self.damping() * du + self.mass() * ddu;
% consider the effect of the element-local external forces
self.internal = T * ( f - self.force() ); % TODO: only static elememt load considered here, time fucntions?!
end
%------------------------------------------------------------------
%% placeholder -> graphical output: plot element
function plot(self, config)
% do nothing here
end
%% placeholder -> print element data
function str = print(self)
str = '';
end
end
end
+mafe/FeAnalysis.m 0 → 100644
View file @ 03febfb1
classdef FeAnalysis < handle
% An base class for any analysis of problems discretised with FE
properties
fep; % the description of the finite element problem
end
methods
%
end
end
+mafe/FeAnalysisDynamicFD.m 0 → 100644
View file @ 03febfb1
classdef FeAnalysisDynamicFD < mafe.FeAnalysis
% dynamic analysis of finite element problem in frequency domain
properties
X = []; % matrix of modes columns-wise
L = []; % vector of eigenvalues
A = []; % vector of coefficients (initial conditions)
Z = []; % column-matrix of particular solutions
P = []; % column-matrix of forcing vectors
M = []; % mass matrix
D = []; % damping matrix
K = []; % stiffness matrix
tfuns = [ mafe.TimeFunction.Static ]; % vector of time functions used
end
methods
%% constructor
function obj = FeAnalysisDynamicFD(fep, tfuns)
if nargin >= 1
obj.fep = fep;
end
if nargin >= 2
obj.tfuns = [ mafe.TimeFunction.Static, tfuns ];
end
end
%% calculate dynamic response
function [u, lambda] = analyse(self)
% make sure the problem is initialised
self.fep.init();
% get system mass and stiffness matrix and force vector
self.M = self.fep.assembleSystemMassMatrix();
self.D = self.fep.assembleSystemDampingMatrix();
self.K = self.fep.assembleSystemStiffnessMatrix();
% identify indices of prescribed Dirchlet constraints
[~,dbc,act] = self.fep.applyConstraintsDirichlet(self.K, zeros(self.fep.ndofs,1));
% solve eigenproblem -----------------------------------------------
% (a) compute eigenvalues lambda (including complex conjugate)
sfac = 1./max( [ normest(self.M,1e-4), normest(self.D,1e-4), normest(self.K,1e-4) ] ); % scale just for polyeig
self.L = polyeig(self.K(act,act)*sfac, self.D(act,act)*sfac, self.M(act,act)*sfac);
% (b) compute eigenvectors
self.X = zeros(self.fep.ndofs,length(self.L));
for ii = 1:length(self.L)
lambda = self.L(ii);
A = full( lambda^2 * self.M(act,act) + lambda * self.D(act,act) + self.K(act,act) );
% compute eigenvector by solution of linear system A(lambda_k)*x_k=e
% b = -A(:,1);
% A(:,1) = 0;
% A(1,:) = 0;
% A(1,1) = 1; b(1) = 1;
% self.X(act,ii) = A\b;
% compute eigenvector by singular value decomposition
% TODO: check if this is generalisable
[U,S,V] = svd(A);
self.X(act,ii) = V(:,end);
% scale computed eigenvector
self.X(:,ii) = 1./sqrt(self.X(:,ii)'*self.M*self.X(:,ii)) * self.X(:,ii); % normalize mode wrt M
end
% % (xx) compute eigenvalues of the first order system
% idx1 = [1:self.fep.ndofs]; idx2 = idx1 + self.fep.ndofs;
% size = self.fep.ndofs*2;
% A = zeros(size,size);
% B = zeros(size,size);
% A(idx1,idx1) = -eye(self.fep.ndofs,self.fep.ndofs);
% A(idx2,idx2) = +M;
% B(idx1,idx2) = +eye(self.fep.ndofs,self.fep.ndofs);
% B(idx2,idx1) = +K;
% B(idx2,idx2) = +D;
% aact = [act act+self.fep.ndofs];
% self.L = eig(B(aact,aact),-A(aact,aact));
% self.X = zeros(self.fep.ndofs,length(self.L));
% for ii = 1:length(self.L)
% lambda = self.L(ii);
% C = lambda * A(aact,aact) + B(aact,aact);
% b = -C(:,1);
% C(:,1) = 0;
% C(1,:) = 0;
% C(1,1) = 1; b(1) = 1;
% tempx = C\b;
% self.X(act,ii) = tempx(1:end/2);
% %norm( ( lambda * A(aact,aact) + B(aact,aact) ) * tempx, 2 )
% %norm( ( lambda^2 * M(act,act) + lambda * D(act,act) + K(act,act) ) * self.X(act,ii), 2 )
% self.X(:,ii) = 1./sqrt(self.X(:,ii)'*M*self.X(:,ii)) * self.X(:,ii); % normalize mode wrt M
% end
% (c) sort eigenvalues/eigenvectors in ascending order
[Lsorted, permutation] = sort(self.L);
self.X = self.X(:,permutation);
self.L = Lsorted;
self.A = ones(length(self.L),1); % put factor=1 as default
% solve particular problem(s)
self.Z = zeros(self.fep.ndofs,length(self.tfuns)); % particular solution vectors
% store forcing vector(s)
self.P = zeros(self.fep.ndofs,length(self.tfuns)); % particular forcing vectors
%
ii = 1;
for tf = self.tfuns
% compute load vector associated with time function tf
p = self.fep.assembleSystemForceVector(tf);
% apply Neumann constraints associated with time function tf
f = self.fep.applyConstraintsNeumann(p, tf);
% apply Dirichlet constraints associated with time function tf
fa = self.fep.applyConstraintsDirichlet(self.M, zeros(self.fep.ndofs,1), tf);
fv = self.fep.applyConstraintsDirichlet(self.D, zeros(self.fep.ndofs,1), tf);
fu = self.fep.applyConstraintsDirichlet(self.K, zeros(self.fep.ndofs,1), tf);
% multiply load vector with evaluated time function
fa = fa * ( -1*tf.Omega^2 );
fv = fv * ( 1i*tf.Omega );
% compose load vector
ff = f + fa + fv + fu;
% calculate system response associated with time function tf
self.Z(act,ii) = ( -tf.Omega^2 * self.M(act,act) + 1i * tf.Omega * self.D(act,act) + self.K(act,act) ) \ ff(act);
% put kinematic constraint values to solution
self.Z(dbc,ii) = fu(dbc);
% put local force vector
self.P(:,ii) = f(:);
%
ii = ii + 1;
end
end
%% apply initial conditions
function applyInitialConditions(self, u0, v0, t0)
u0 = u0;
%
if ~exist('v0', 'var')
v0 = zeros(self.fep.ndofs,1);
end
if ~exist('t0', 'var')
t0 = 0.0;
end
%
C = zeros(2*self.fep.ndofs, length(self.L));
%
for ii = 1:length(self.L)
C(1+0*self.fep.ndofs:1*self.fep.ndofs,ii) = self.X(:,ii) * exp(self.L(ii)*t0);
C(1+1*self.fep.ndofs:2*self.fep.ndofs,ii) = self.X(:,ii) * self.L(ii) * exp(self.L(ii)*t0);
end
self.A = C\( [u0; v0] - [self.solInhomogeneous(t0,0); self.solInhomogeneous(t0,1)] );
%
% plot the amplitude spectrum of complex conjugate pairs
% clf; hold on;
% plot( 0.5*real(self.A(1:2:end)), '*b' );
% plot( 0.5*imag(self.A(1:2:end)), 'xr' );
% pause;
end
%% compute total solution to homogeneous equation of motion at given time t
function xh = solHomogeneous(self, time, difforder)
if ~exist('time', 'var')
time = 0.0;
end
if ~exist('difforder', 'var')
difforder = 0;
end
xh = zeros(self.fep.ndofs,1);
switch(difforder)
case 0
for ii = 1:length(self.L)
% component ii of the homogeneous solution x_h
% is typically a complex vector
xh = xh + self.A(ii) * exp(self.L(ii)*time) * self.X(:,ii);
end
case 1
for ii = 1:length(self.L)
% component ii of the homogeneous solution dot x_h
% is typically a complex vector
xh = xh + self.A(ii) * self.L(ii) *exp(self.L(ii)*time) * self.X(:,ii);
end
case 2
for ii = 1:length(self.L)
% component ii of the homogeneous solution ddot x_h
% is typically a complex vector
xh = xh + self.A(ii) * self.L(ii)^2*exp(self.L(ii)*time) * self.X(:,ii);
end
otherwise
disp('Warning: Wrong choice for time diff order!')
end
% 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!')
end
% return the real part only
xh = real(xh);
end
%% compute solution to particular equation of motion at given time t
function xp = solInhomogeneous(self, time, difforder)
if ~exist('time', 'var')
time = 0.0;
end
if ~exist('difforder', 'var')
difforder = 0;
end
xp = zeros(self.fep.ndofs,1);
switch(difforder)
case 0
for ii = 1:length(self.tfuns)
tf = self.tfuns(ii);
xp = xp + self.Z(:,ii) * tf.eval_f(time);
end
case 1
for ii = 1:length(self.tfuns)
tf = self.tfuns(ii);
xp = xp + self.Z(:,ii) * tf.eval_dotf(time);
end
case 2
for ii = 1:length(self.tfuns)
tf = self.tfuns(ii);
xp = xp + self.Z(:,ii) * tf.eval_ddotf(time);
end
otherwise
disp('Warning: Wrong choice for time diff order!')
end
% 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);
% --- POST-PROCESSING ---
% get local force vector at specific time
p = zeros(self.fep.ndofs,1);
for ii = 1:length(self.tfuns)
tf = self.tfuns(ii);
p = p + self.P(:,ii) * tf.eval_f(time);
end
p = real(p); % transform complex notation to real
% compute reaction forces (including dynamic effects)
r = self.M * a + self.D * v + self.K * u - p;
% feedback of solution values (u and r) to nodal dofs
self.fep.updateNodalDofs(u, r);
% postprocess elements
self.fep.postprocessElements(u,v,a);
end
%% put selected mode to dof solution
function putModeToDof(self, mode, component)
if nargin < 1
mode = 1;
end
if nargin < 2
component = 'real';
end
lambda = self.L(mode)
a = self.A(mode)
switch component
case 'imag'
u = imag( a * self.X(:,mode) );
case 'real'
u = real( a * self.X(:,mode) );
otherwise
u = a * self.X(:,mode);
end
u = 1./norm(u,2) * u
% compute reaction forces
K = self.fep.assembleSystemStiffnessMatrix();
f = self.fep.assembleSystemForceVector();
r = K * u - f;
% feedback of solution values (u and r) to nodal dofs
self.fep.updateNodalDofs(u, r);
% postprocess elements
self.fep.postprocessElements(u);
end
end
end
+mafe/FeAnalysisDynamicTD.m 0 → 100644
View file @ 03febfb1
classdef FeAnalysisDynamicTD < mafe.FeAnalysis
% dynamic analysis of finite element problem in time domain
properties
L = []; % lower part of decomposition of the effective system matrix
U = []; % upper part of decomposition of the effective system matrix
M = [];
D = [];
K = [];
Mf = [];
Df = [];
Kf = [];
u = []; % vector of displacement at t_{n+1}
v = []; % vector of velocity at t_{n+1}
a = []; % vector of acceleration at t_{n+1}
tn = 0.; % time instant at t_{n}
dt = 1.; % time step size
scheme = mafe.TimeIntegrationType.CrankNicholson; % default scheme
params; % parameter set of the time integration scheme
tfuns = [ mafe.TimeFunction.Static ]; % vector of time functions used
end
methods
%% constructor
function obj = FeAnalysisDynamicTD(fep, tfuns)
if nargin >= 1
obj.fep = fep;
end
if nargin >= 2
obj.tfuns = [ mafe.TimeFunction.Static, tfuns ];
end
end
%% initialise the calculation of the dynamic response
function [] = initialise(self, dt, scheme, params)
if nargin >= 2
self.dt = dt;
end
if nargin >= 3
self.scheme = scheme;
end
if nargin >= 4
self.params = params;
end
% ------------------------------------------------------------------
% scheme identification and settings
switch (self.scheme)
case mafe.TimeIntegrationType.CrankNicholson
%
beta = 0.25;
gamma = 0.50;
%
self.params.theta_ = [ 1.0, 1.0, 1.0 ];
self.params.gamma_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
self.params.beta_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
%
case mafe.TimeIntegrationType.CentralDifference
%
beta = 0.00;
gamma = 0.50;
%
self.params.theta_ = [ 1.0, 1.0, 1.0 ];
self.params.gamma_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
self.params.beta_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
%
case mafe.TimeIntegrationType.LinearAcceleration
%
beta = 1./6;
gamma = 0.50;
%
self.params.theta_ = [ 1.0, 1.0, 1.0 ];
self.params.gamma_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
self.params.beta_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
%
case mafe.TimeIntegrationType.FoxGoodwin
%
beta = 1/12;
gamma = 0.50;
%
self.params.theta_ = [ 1.0, 1.0, 1.0 ];
self.params.gamma_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
self.params.beta_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
%
case mafe.TimeIntegrationType.GeneralisedNewmark
%
beta = self.params.beta;
gamma = self.params.gamma;
%
self.params.theta_ = [ 1.0, 1.0, 1.0 ];
self.params.gamma_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
self.params.beta_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
%
case mafe.TimeIntegrationType.GeneralisedAlpha
%
beta = self.params.beta;
gamma = self.params.gamma;
alpha_f = self.params.alpha_f;
alpha_m = self.params.alpha_m;
%
self.params.theta_ = [ 1.0-alpha_f, 1.0-alpha_f, 1.0-alpha_f ];
self.params.gamma_ = [ 6.0*beta, 2.0*gamma, (1.0-alpha_m)/(1.0-alpha_f) ];
self.params.beta_ = [ 6.0*beta, 2.0*gamma, 1.0 ];
%
case mafe.TimeIntegrationType.Bathe
%
% TODO --> check if this scheme can be fit into the above!?
%
otherwise
disp('FeAnalysisDynamicTD::initialise() : Unknown scheme!');
return;
end
% ------------------------------------------------------------------
% make sure the problem is initialised
self.fep.init();
% get system mass and stiffness matrix and force vector
M = self.fep.assembleSystemMassMatrix();
D = self.fep.assembleSystemDampingMatrix();
K = self.fep.assembleSystemStiffnessMatrix();
f = self.fep.assembleSystemForceVector();
% store system matrices
self.Mf = M;
self.Df = D;
self.Kf = K;
% apply constraints
[~, ~, ~, M] = self.fep.applyConstraintsDirichlet(M, zeros(self.fep.ndofs,1));
[~, ~, ~, D] = self.fep.applyConstraintsDirichlet(D, zeros(self.fep.ndofs,1));
[~, ~, ~, K] = self.fep.applyConstraintsDirichlet(K, zeros(self.fep.ndofs,1));
% ------------------------------------------------------------------
A = ( self.params.gamma_(3)*self.params.theta_(1) ) * M;
A = A + ( self.params.gamma_(2)*self.params.theta_(2)*(self.dt)^1/2. ) * D;
A = A + ( self.params.gamma_(1)*self.params.theta_(3)*(self.dt)^2/6. ) * K;
%
[self.L, self.U] = lu(A);
% ------------------------------------------------------------------
% set initial conditions to default zero
self.u = zeros(self.fep.ndofs,1);
self.v = zeros(self.fep.ndofs,1);
self.a = zeros(self.fep.ndofs,1);
% ------------------------------------------------------------------
% store system matrices
self.M = M;
self.D = D;
self.K = K;
end
%% perform one time integration step
function solveTimeStep(self)
% compute estimates
da_tilde = self.a * 0.0;
u_tilde = self.u;
u_tilde = u_tilde + ( self.params.theta_(1)*(self.dt) ) * self.v;
u_tilde = u_tilde + ( self.params.theta_(2)*(self.dt)^2/2. ) * self.a;
u_tilde = u_tilde + ( self.params.theta_(3)*(self.dt)^3/6.*self.params.gamma_(1) ) * da_tilde / self.dt;
v_tilde = self.v;
v_tilde = v_tilde + ( self.params.theta_(1)*(self.dt) ) * self.a;
v_tilde = v_tilde + ( self.params.theta_(2)*(self.dt)^2/2.*self.params.gamma_(2) ) * da_tilde / self.dt;
a_tilde = self.a;
a_tilde = a_tilde + ( self.params.theta_(1)*(self.dt)^2/1.*self.params.gamma_(3) ) * da_tilde / self.dt;
% compute forcing vector at current time instant (using time fcts)
fp = zeros(self.fep.ndofs,1);
fa = zeros(self.fep.ndofs,1);
fv = zeros(self.fep.ndofs,1);
fu = zeros(self.fep.ndofs,1);
% current and next time instant
time = self.tn;
tnext = time + self.dt;
% identify indices of prescribed Dirchlet constraints
[~,dbc] = self.fep.applyConstraintsDirichlet(self.Kf, fu);
% load contributions from stationary and non-stationary loads and constraints
for tf = self.tfuns
% compute load vector associated with time function tf
fp_ = self.fep.assembleSystemForceVector(tf);
% apply Neumann constraints associated with time function tf
fp_ = self.fep.applyConstraintsNeumann(fp_, tf);
% apply Dirichlet constraints associated with time function tf
fa_ = self.fep.applyConstraintsDirichlet(self.Mf, zeros(self.fep.ndofs,1), tf);
fv_ = self.fep.applyConstraintsDirichlet(self.Df, zeros(self.fep.ndofs,1), tf);
fu_ = self.fep.applyConstraintsDirichlet(self.Kf, zeros(self.fep.ndofs,1), tf);
% multiply load vector with evaluated time function
fa = fa + fa_ * ( (1.0-self.params.theta_(1))*( tf.eval_cos_ddotf(time) + tf.eval_sin_ddotf(time) ) + self.params.theta_(1)*( tf.eval_cos_ddotf(tnext) + tf.eval_sin_ddotf(tnext) ) );
fv = fv + fv_ * ( (1.0-self.params.theta_(1))*( tf.eval_cos_dotf(time) + tf.eval_sin_dotf(time) ) + self.params.theta_(1)*( tf.eval_cos_dotf(tnext) + tf.eval_sin_dotf(tnext) ) );
fu = fu + fu_ * ( (1.0-self.params.theta_(1))*( tf.eval_cos_f(time) + tf.eval_sin_f(time) ) + self.params.theta_(1)*( tf.eval_cos_f(tnext) + tf.eval_sin_f(tnext) ) );
fp = fp + fp_ * ( (1.0-self.params.theta_(1))*( tf.eval_cos_f(time) + tf.eval_sin_f(time) ) + self.params.theta_(1)*( tf.eval_cos_f(tnext) + tf.eval_sin_f(tnext) ) );
end
% compose load vector
p = fp + fa + fv + fu;
% compute resudial force vector
r = p - ( self.M * a_tilde + self.D * v_tilde + self.K * u_tilde );
% forward/backward solve re-using factorisation from initialisation
y = self.L\r;
x = self.U\y;
% update state vectors
da = da_tilde + x;
self.u = self.u + ( (self.dt) ) * self.v;
self.u = self.u + ( (self.dt)^2/2. ) * self.a;
self.u = self.u + ( self.params.beta_(1)*(self.dt)^3/6. ) * da / self.dt;
self.v = self.v + ( (self.dt) ) * self.a;
self.v = self.v + ( self.params.beta_(2)*(self.dt)^2/2. ) * da / self.dt;
self.a = self.a + ( self.params.beta_(3)*(self.dt)^1/1. ) * da / self.dt;
% include given dirichlet values
self.u(dbc) = fu(dbc);
self.v(dbc) = fv(dbc);
self.a(dbc) = fa(dbc);
% update time
self.tn = self.tn + self.dt;
% --- POST-PROCESSING ---
% compute reaction forces (including dynamic effects)
r = self.Mf * self.a + self.Df * self.v + self.Kf * self.u - fp;
% feedback of solution values (u and r) to nodal dofs
self.fep.updateNodalDofs(self.u, r);
% postprocess elements
self.fep.postprocessElements(self.u, self.v, self.a);
end
%% apply initial conditions
function applyInitialConditions(self, u0, v0, t0)
self.u = u0;
%
if ~exist('v0', 'var')
self.v = zeros(self.fep.ndofs,1);
end
if ~exist('t0', 'var')
self.tn = 0.0;
end
end
end
end
+mafe/FeAnalysisStatic.m 0 → 100644
View file @ 03febfb1
classdef FeAnalysisStatic < mafe.FeAnalysis
% static analysis of finite element problem
properties
%
end
methods
%% constructor
function obj = FeAnalysisStatic(fep)
if nargin > 0
obj.fep = fep;
end
end
%% calculate static response
function analyse(self)
% make sure the problem is initialised
self.fep.init();
% get system stiffness matrix and force vector
K = self.fep.assembleSystemStiffnessMatrix();
p = self.fep.assembleSystemForceVector();
% apply constraints
[f ] = self.fep.applyConstraintsNeumann(p);
[f, dbc, act] = self.fep.applyConstraintsDirichlet(K, f);
% allocate solution vector
u = zeros(self.fep.ndofs,1);
% solve system of equations (active indices only)
u(act) = K(act,act)\f(act);
u(dbc) = f(dbc);
% compute reaction forces
r = K * u - p;
% feedback of solution values (u and r) to nodal dofs
self.fep.updateNodalDofs(u, r);
% postprocess elements
self.fep.postprocessElements(u);
end
end
end
+mafe/FeAnalysisSteadyNonlinear.m 0 → 100644
View file @ 03febfb1
classdef FeAnalysisSteadyNonlinear < mafe.FeAnalysis
% steady state nonlinear analysis of finite element problem
properties
%
end
methods
%% constructor
function obj = FeAnalysisSteadyNonlinear(fep)
if nargin > 0
obj.fep = fep;
end
end
%% calculate static response
function analyse(self)
% make sure the problem is initialised
self.fep.init();
% allocate solution vector
x = zeros(self.fep.ndofs,1);
dx = zeros(self.fep.ndofs,1);
% get system tangent matrix and force vector
self.fep.updateNodalDofs(x, 0*x);
L = self.fep.assembleSystemMatrix_Linear_dx();
p = -L*x;
% apply constraints
[f ] = self.fep.applyConstraintsNeumann(p);
[f, dbc, act] = self.fep.applyConstraintsDirichlet(L, f);
ii = 0;
for lambda = linspace(0.0, 1.0, 10)
x(dbc) = f(dbc)*lambda;
self.fep.updateNodalDofs(x, 0*x);
for iter = 1:12
% get system tangent matrix and force vector
tic
A = L + self.fep.assembleSystemMatrix_Nonlinear_dx();
toc
tic
p = -L*x - self.fep.assembleSystemForceVector();
toc
% apply constraints
[f ] = self.fep.applyConstraintsNeumann(p);
[f, dbc, act] = self.fep.applyConstraintsDirichlet(A, f);
% solve system of equations (active indices only)
dx(act) = A(act,act)\f(act);
if iter == 1
dx(dbc) = f(dbc)*lambda;
else
dx(dbc) = 0.0;
end
% update solution
x = x + dx;
% feedback of solution values (x and r) to nodal dofs
self.fep.updateNodalDofs(x, A*dx-p);
% cla; clf;
% fig = figure(1); subplot('Position',[0.05 0.05 0.90 0.90]); axis equal; hold on;
% self.fep.plotSystem('pressure');
% self.fep.plotSystem('velocity', 1.e-4);
% saveas(gcf,sprintf('nscyl_%04d',ii),'png');
ii = ii + 1;
% check
n = norm(dx);
%n = dx'*p;
fprintf( sprintf('lambda = %5.5f | iter = %2d | norm = %4.6e \n', lambda, iter, n) );
if (n < 1.e-9)
break;
end
end
end
% postprocess elements
self.fep.postprocessElements(x);
end
end
end
+mafe/FeProblem.m 0 → 100644
View file @ 03febfb1
classdef FeProblem < handle
% A generic FE problem description
properties
nodes = mafe.Node.empty(0,0); % the vector of nodes of the system
elems = mafe.Element.empty(0,0); % the vector of elements of the system
const = mafe.Constraint.empty(0,0); % the vector of constraints of the system
ndofs = -1; % the total number of degree of freedom in the system
end
methods
%% constructor
function obj = FeProblem(nodes, elems, const)
if nargin >= 3
obj.nodes = nodes;
obj.elems = elems;
obj.const = const;
end
% call initialisation of finite element problem
obj.init();
end
%% initialise finite element problem
function init(self)
% loop over all elements and setup dof types
for elem = self.elems
elem.activateDofTypes();
end
% loop over all nodes and setup system indices of dofs
self.ndofs = 0;
for node = self.nodes
node.dof = unique(node.dof,'sorted'); %sort(unique(node.dof));
ndof = length(node.dof);
node.idx = self.ndofs + [1:ndof];
self.ndofs = self.ndofs + ndof;
end
%sprintf('feproblem ndofs = %d', self.ndofs)
end
%% assemble system K matrix
function A = assembleSystemStiffnessMatrix(self)
tic
% 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) );
toc
return
tic
% create system stiffness matrix
id1 = zeros(0,1);
id2 = zeros(0,1);
mat = zeros(0,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();
% add element contribution to system matrix
id1 = [ id1, repmat( e_idx, 1, length(e_idx) ) ]; % TODO: clean up here
id2 = [ id2; reshape( repmat( e_idx, length(e_idx), 1 ), [], 1 ) ];
mat = [ mat; reshape( e_mat', [], 1 ) ];
end
A = sparse( id1, id2, mat );
toc
tic
% create system mass matrix
A = sparse(self.ndofs, self.ndofs);
% loop over all elements and collect element mass matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.stiffness();
% add element contribution to system matrix
A(e_idx,e_idx) = A(e_idx,e_idx) + e_mat;
end
toc
end
%% assemble system M matrix
function A = assembleSystemMassMatrix(self)
% create system mass matrix
A = sparse(self.ndofs, self.ndofs);
% loop over all elements and collect element mass matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.mass();
% add element contribution to system matrix
A(e_idx,e_idx) = A(e_idx,e_idx) + e_mat;
end
end
%% assemble system D matrix
function A = assembleSystemDampingMatrix(self)
% create system damping matrix
A = sparse(self.ndofs, self.ndofs);
% loop over all elements and collect element damping matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.damping();
% add element contribution to system matrix
A(e_idx,e_idx) = A(e_idx,e_idx) + e_mat;
end
end
%% assemble system f vector
function f = assembleSystemForceVector(self, tfun)
% create system force vector
f = zeros(self.ndofs, 1);
% loop over all elements and collect element mass matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element vector
if ~exist('tfun', 'var')
e_vec = elem.force(); % tfun not given, collect all
elseif isprop(elem,'load')
if elem.load.tfun == tfun
e_vec = elem.force(); % collect matching tfun
else
e_vec = elem.force() * 0.0; % collect zero
end
else
e_vec = zeros(length(e_idx),1);
end
% add element contribution to system vector
f(e_idx) = f(e_idx) + e_vec;
end
end
%% apply constraints of type Neumann to the rhs
function [b, on, off] = applyConstraintsNeumann(self, b, tfun)
% apply constraints - dynamic/static
idx = [];
val = [];
for con = self.const
if con.type == mafe.ConstraintType.Neumann
idx = [ idx, con.getDofSysIndices() ];
if ~exist('tfun', 'var')
% tfun not given -> time-constant dynamic/static constraint
val = [ val, con.value() ];
elseif con.tfun == tfun
% collect matching tfun value
val = [ val, con.value() ];
else
% collect zero (vector value)
val = [ val, 0.0*con.value() ];
end
end
end
% add values and consider repeated indices:
if length(idx) > 0 % check for zero-length idx because of accumarray
b = b + accumarray(idx',val', [length(b) 1]);
end
% create index vector of dof system indices without Neumann constraints
off = setdiff([1:self.ndofs], idx);
% create index vector of dof system indices with Neumann constraints
on = idx;
end
%% apply constraints of type Dirichlet to coefficent matrix and rhs
function [b, on, off, A] = applyConstraintsDirichlet(self, A, b, tfun)
% apply constraints - kinematic
idx = [];
val = [];
for con = self.const
if con.type == mafe.ConstraintType.Dirichlet
idx = [ idx, con.getDofSysIndices() ];
if ~exist('tfun', 'var')
% tfun not given -> time-constant kinematic constraint
val = [ val, con.value() ];
elseif con.tfun == tfun
% collect matching tfun value
val = [ val, con.value() ];
else
% collect zero, treat this as homogeneous constraint
val = [ val, 0.0*con.value() ];
end
end
end
% check for repeated indices in kinematic constraints
if length(idx) ~= length(unique(idx))
%disp('WARNING: FeProblem::applyConstraints() : Repeated indices in Dirichlet constraint !');
% give preference to first occurence
[~,occurence,~] = unique(idx, 'first');
val = val(occurence);
idx = idx(occurence);
end
% impose inhomogeneous Dirchlet conditions
b = b - A(:,idx) * val'; % assuming distinct indices here
A(idx,:) = 0.0; A(:,idx) = 0.0; A(idx,idx) = eye(length(idx));
b(idx) = val;
% create index vector of dof system indices without Dirichlet constraints
off = setdiff([1:self.ndofs], idx);
% create index vector of dof system indices with Dirichlet constraints
on = idx;
end
%% update solution values to nodal dofs
function updateNodalDofs(self, lhs, rhs)
for node = self.nodes
node.lhs = lhs(node.idx);
node.rhs = rhs(node.idx);
end
end
%% postprocess elements
function postprocessElements(self, u, v, a)
if ~exist('v', 'var')
v = zeros(self.ndofs,1);
end
if ~exist('a', 'var')
a = zeros(self.ndofs,1);
end
for elem = self.elems
idx = elem.getDofSysIndices();
elem.postprocess( u(idx), v(idx), a(idx) );
end
end
%% plot the system
function plotSystem(self, config, scale)
if ~exist('scale', 'var')
scale = 1.0;
end
if ~exist('config', 'var')
config = 'reference';
end
% loop over elements and draw them
for elem = self.elems
elem.plot(config, scale);
end
end
%% print the system
function printSystem(self)
% loop over nodes and print them
disp('nodal data:')
for node = self.nodes
fprintf( [ sprintf('node %3d: ', i), node.print(), '\n' ] );
end
% loop over elements and print them
disp('element data:')
for elem = self.elems
fprintf( [ sprintf('elem %3d: ', i), elem.print(), '\n' ] );
end
end
end
end
+mafe/FeProblemLinearised.m 0 → 100644
View file @ 03febfb1
classdef FeProblemLinearised < mafe.FeProblem
% A generic FE problem description
properties
%
end
methods
%% constructor
function obj = FeProblemLinearised(varargin)
obj = obj@mafe.FeProblem(varargin{:});
end
%% assemble system matrix related to dx
function A = assembleSystemMatrix_Linear_dx(self)
% % create system tangent matrix
% id1 = zeros(0,1);
% id2 = zeros(0,1);
% mat = zeros(0,1);
% % loop over all elements and collect element matrices
% for elem = self.elems
% % get element index vector
% e_idx = elem.getDofSysIndices();
% % get element matrix
% e_mat = elem.elemat_dx( x(e_idx) );
% % add element contribution to system matrix
% id1 = [ id1, repmat( e_idx, 1, length(e_idx) ) ]; % TODO: clean up here
% id2 = [ id2; reshape( repmat( e_idx, length(e_idx), 1 ), [], 1 ) ];
% mat = [ mat; reshape( e_mat', [], 1 ) ];
% end
% A = sparse( id1, id2, mat );
% create system damping matrix
A = sparse(self.ndofs, self.ndofs);
% loop over all elements and collect element damping matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.elemat_linear_dx();
% add element contribution to system matrix
A(e_idx,e_idx) = A(e_idx,e_idx) + e_mat;
end
end
%% assemble system matrix related to dx
function A = assembleSystemMatrix_Nonlinear_dx(self)
% % create system tangent matrix
% id1 = zeros(0,1);
% id2 = zeros(0,1);
% mat = zeros(0,1);
% % loop over all elements and collect element matrices
% for elem = self.elems
% % get element index vector
% e_idx = elem.getDofSysIndices();
% % get element matrix
% e_mat = elem.elemat_dx( x(e_idx) );
% % add element contribution to system matrix
% id1 = [ id1, repmat( e_idx, 1, length(e_idx) ) ]; % TODO: clean up here
% id2 = [ id2; reshape( repmat( e_idx, length(e_idx), 1 ), [], 1 ) ];
% mat = [ mat; reshape( e_mat', [], 1 ) ];
% end
% A = sparse( id1, id2, mat );
% create system damping matrix
A = sparse(self.ndofs, self.ndofs);
% loop over all elements and collect element damping matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.elemat_nonlinear_dx();
% add element contribution to system matrix
A(e_idx,e_idx) = A(e_idx,e_idx) + e_mat;
end
end
%% assemble system matrix related to dot_dx
function A = assembleSystemMatrix_dot_dx(self, x_lastiter)
% create system tangent matrix
id1 = zeros(0,1);
id2 = zeros(0,1);
mat = zeros(0,1);
% loop over all elements and collect element matrices
for elem = self.elems
% get element index vector
e_idx = elem.getDofSysIndices();
% get element matrix
e_mat = elem.elemat_dot_dx( x_lastiter(e_idx) );
% add element contribution to system matrix
id1 = [ id1, repmat( e_idx, 1, length(e_idx) ) ]; % TODO: clean up here
id2 = [ id2; reshape( repmat( e_idx, length(e_idx), 1 ), [], 1 ) ];
mat = [ mat; reshape( e_mat', [], 1 ) ];
end
A = sparse( id1, id2, mat );
end
%% postprocess elements
function postprocessElements(self, x, dot_x, ddot_x)
if ~exist('x', 'var')
x = zeros(self.ndofs,1);
end
if ~exist('dot_x', 'var')
dot_x = zeros(self.ndofs,1);
end
if ~exist('ddot_x', 'var')
ddot_x = zeros(self.ndofs,1);
end
for elem = self.elems
idx = elem.getDofSysIndices();
elem.postprocess( x(idx), dot_x(idx), ddot_x(idx) );
end
end
end
end
+mafe/GaussCoordWeights.m 0 → 100644
View file @ 03febfb1
%% Gauss integration parameters
function [GC, GW] = GaussCoordWeights(num)
% Gauss integration parameters, local coordinate system z = [-1,+1]
switch(num)
case 1
GC = [ 0. ]; % 1D Gauss coordinates: up to P(1)
GW = [ 2. ]; % 1D Gauss weights: sum=2
case 2
GC = [-sqrt(1./3.), +sqrt(1./3.)]; % up to P(3)
GW = [ 1., 1.]; %
case 3
GC = [-sqrt(3./5.), 0., +sqrt(3./5.)]; % up to P(5)
GW = [ 5./9., 8./9., 5./9.]; %
otherwise
disp('Quadrature.GaussCoordWeights(num) : Order is not implemented!')
end
end
+mafe/Li2Element1D.m 0 → 100644
View file @ 03febfb1
classdef Li2Element1D < mafe.Element
% basis for linear 2-node elements in 1D
properties
node = mafe.Node1D.empty(0,2); % the two nodes
sect = mafe.Section1D(); % the section descriptor
load = mafe.EleLoad(); % the element load descriptor
end
methods
%% initialisation
function obj = Li2Element1D(nodes, section, eleload)
if exist('nodes','var')
obj.node = nodes;
end
if exist('section','var')
obj.sect = section;
end
if exist('eleload','var')
obj.load = eleload;
end
end
% ======================================================================
%% matrix of nodal coordinates
function X = X(self)
X = reshape( [ self.node.ref ], 1, 2 )'; % (2x1) = [ x1_i ]
end
%% linear ansatz function in local coordinates z
function N = N(~, z)
% (1x2) = [N_1, N_2]
N = 0.5 * [ 1-z(1), 1+z(1) ];
end
%% 1st derivative of ansatz wrt local coordinate z
function Nz = Nz(~, z)
% (1x2) = [N_1,z1, N_2,z1 ]
Nz = 0.5 * [ -1, +1 ];
end
%% Hermite ansatz function in local coordinates z
function H = H(~, z, L)
% (1x4) = [H_1, H_2, H_3, H_4]
N = 0.25 * [ (1-z(1))^2 * (2+z(1)), (1-z(1))^2 * (1+z(1)) * L/2, (1+z(1))^2 * (2-z(1)), (1+z(1))^2 * (-1+z(1)) * L/2 ];
end
%% 2nd derivative of Hermite ansatz wrt local coordinate z
function Hzz = Hzz(~, z, L)
% (1x4) = [H_1,z1, H_2,z1, H_3,z1, H_4,z1 ]
Nz = 0.125 * [ 12*z(1), L*(6*z(1)-2), -12*z(1), L*(6*z(1)+2) ] ;
end
%% Jacobi matrix J = J_ij = d x_j / d z_i
function J = J(self, z)
J = self.Nz(z) * self.X();
end
% ======================================================================
%-----------------------------------------------------------------------
%% graphical output: plot element
function plot(self, config, scale)
% my nodes
X = self.X();
% coordinates
x1 = X(:,1);
% plot defaults
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
lwidth = 0.5;
mwidth = 10.;
ecolor = [0.4 0.4 0.4];
marker = '.';
%
switch config
otherwise
XX1 = x1;
XX2 = x1*0;
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
falpha = 1.0;
lwidth = 0.5;
mwidth = 12.;
end
mafe.patchline(XX1, XX2, ...
'EdgeColor', lcolor, ...
'LineWidth', lwidth, ...
'EdgeAlpha', falpha);
plot(XX1([1 end]), XX2([1 end]), ...
marker, ...
'color', mcolor, ...
'markersize', mwidth);
end
%% print element data
function str = print(self)
%
str = [ 'No local data ! (Not implemented?)' ];
end
end
end
+mafe/Li2Element2D.m 0 → 100644
View file @ 03febfb1
classdef Li2Element2D < mafe.Element
% basis for linear 2-node elements in 2D
properties
node = mafe.Node2D.empty(0,2); % the two nodes
sect = mafe.Section1D(); % the section descriptor
load = mafe.EleLoad(); % the element load descriptor
end
methods
%% initialisation
function obj = Li2Element2D(nodes, section, eleload)
if exist('nodes','var')
obj.node = nodes;
end
if exist('section','var')
obj.sect = section;
end
if exist('eleload','var')
obj.load = eleload;
end
end
% ======================================================================
%% matrix of nodal coordinates
function X = X(self)
X = reshape( [ self.node.ref ], 2, 2 ); % (2x2) = [ x1_i ]
end
%% linear ansatz function in local coordinates z
function N = N(~, z)
% (1x2) = [N_1, N_2]
N = 0.5 * [ 1-z(1), 1+z(1) ];
end
%% 1st derivative of ansatz wrt local coordinate z
function Nz = Nz(~, z)
% (1x2) = [N_1,z1, N_2,z1 ]
Nz = 0.5 * [ -1, +1 ];
end
%% Hermite ansatz function in local coordinates z
function H = H(~, z, L)
% (1x4) = [H_1, H_2, H_3, H_4]
N = 0.25 * [ (1-z(1))^2 * (2+z(1)), (1-z(1))^2 * (1+z(1)) * L/2, (1+z(1))^2 * (2-z(1)), (1+z(1))^2 * (-1+z(1)) * L/2 ];
end
%% 2nd derivative of Hermite ansatz wrt local coordinate z
function Hzz = Hzz(~, z, L)
% (1x4) = [H_1,z1, H_2,z1, H_3,z1, H_4,z1 ]
Nz = 0.125 * [ 12*z(1), L*(6*z(1)-2), -12*z(1), L*(6*z(1)+2) ] ;
end
%% Jacobi matrix J = J_ij = d x_j / d z_i
function J = J(self, z)
J = self.Nz(z) * self.X();
end
% ======================================================================
%-----------------------------------------------------------------------
%% graphical output: plot element
function plot(self, config, scale)
% my nodes
X = self.X();
% coordinates
x1 = X(1,:);
x2 = X(2,:);
% plot defaults
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
lwidth = 0.5;
mwidth = 10.;
ecolor = [0.4 0.4 0.4];
marker = '.';
%
switch config
otherwise
XX1 = x1;
XX2 = x2;
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
falpha = 1.0;
lwidth = 0.5;
mwidth = 12.;
end
mafe.patchline(XX1, XX2, ...
'EdgeColor', lcolor, ...
'LineWidth', lwidth, ...
'EdgeAlpha', falpha);
plot(XX1([1 end]), XX2([1 end]), ...
marker, ...
'color', mcolor, ...
'markersize', mwidth);
end
%% print element data
function str = print(self)
%
str = [ 'No local data ! (Not implemented?)' ];
end
end
end
+mafe/Member2D.m 0 → 100644
View file @ 03febfb1
classdef Member2D < mafe.Li2Element2D
% member (= beam + truss) finite element in the plane (2D)
properties
%
end
methods
%% constructor
function obj = Member2D(varargin)
obj = obj@mafe.Li2Element2D(varargin{:});
end
%% activation of dofs at the nodes of the element
function activateDofTypes(self)
ineed = [ mafe.DofType.Disp1 mafe.DofType.Disp2 mafe.DofType.Rota3 ];
for i = 1:length(self.node)
self.node(i).dof = [self.node(i).dof ineed];
end
end
%% provide element index vector of dof
function idx = getDofSysIndices(self)
disp1 = find( [self.node.dof] == mafe.DofType.Disp1 );
disp2 = find( [self.node.dof] == mafe.DofType.Disp2 );
rota3 = find( [self.node.dof] == mafe.DofType.Rota3 );
ndidx = [ self.node.idx ];
idx = [ ndidx(disp1(1)), ndidx(disp2(1)), ndidx(rota3(1)), ndidx(disp1(2)), ndidx(disp2(2)), ndidx(rota3(2)) ];
end
%% calculate transformation matrix local->global
function [L, T] = transform(self)
dx1 = self.node(2).ref(1) - self.node(1).ref(1);
dx2 = self.node(2).ref(2) - self.node(1).ref(2);
L = sqrt( dx1^2 + dx2^2 );
T = 1/L*[ dx1, dx2, 0, 0, 0, 0; ...
-dx2, dx1, 0, 0, 0, 0; ...
0, 0, L, 0, 0, 0; ...
0, 0, 0, dx1, dx2, 0; ...
0, 0, 0,-dx2, dx1, 0; ...
0, 0, 0, 0, 0, L ];
end
%% calculation of element stiffness matrix
function ele_mat = stiffness(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
EA = self.sect.E * self.sect.A ;
EI = self.sect.E * self.sect.Iy;
% local element matrix
h1 = 2*EI/L;
h2 = 6*EI/L^2;
h3 = 12*EI/L^3;
h4 = EA/L;
ele_mat_local = [ +h4, 0, 0, -h4, 0, 0 ;
0, +h3, +h2, 0, -h3, +h2 ;
0, +h2, 2*h1, 0, -h2, +h1 ;
-h4, 0, 0, +h4, 0, 0 ;
0, -h3, -h2, 0, +h3, -h2 ;
0, +h2, +h1, 0, -h2, 2*h1 ];
% global element matrix
ele_mat = T' * ele_mat_local * T;
end
%% calculation of element mass matrix
function ele_mat = mass(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
rhoA = self.sect.rho * self.sect.A;
% local element matrix
h1 = rhoA/420 * L^3;
h2 = rhoA/420 * L^2;
h3 = rhoA/420 * L;
h4 = rhoA/6 * L;
ele_mat_local = [ 2*h4, 0, 0, 1*h4, 0, 0 ; ...
0, 156*h3, 22*h2, 0, 54*h3, -13*h2 ; ...
0, 22*h2, 4*h1, 0, 13*h2, -3*h1 ; ...
1*h4, 0, 0, 2*h4, 0, 0 ; ...
0, 54*h3, 13*h2, 0, 156*h3, -22*h2 ; ...
0, -13*h2, -3*h1, 0, -22*h2, 4*h1 ];
% global element matrix
ele_mat = T' * ele_mat_local * T;
end
%% calculation of element damping matrix
function ele_mat = damping(self)
% provide damping matrix based on Rayleigh assumption
% get parameters
d_alpha = self.sect.d_alpha;
d_beta = self.sect.d_beta;
% D = alpha * M + beta * K
ele_mat = d_alpha * self.mass() + d_beta * self.stiffness();
end
%% calculation of element force vector
function ele_vec = force(self)
% transformation matrix
[L, T] = self.transform();
% get parameters
px1 = self.load.p(mafe.DofType.Disp1,:);
px2 = self.load.p(mafe.DofType.Disp2,:);
pxr = self.load.p(mafe.DofType.Rota3,:);
% local element vector
pX([1,4]) = T(1,1)*px1 + T(1,2)*px2;
pX([2,5]) = T(2,1)*px1 + T(2,2)*px2;
pX([3,6]) = pxr;
% local element vector
ele_vec_local = L/6./210 * [ 420, 0, 0, 210, 0, 0; ...
0, 468, -630/L, 0, 162, -630/L; ...
0, L*66, 126, 0, L*39, -126; ...
210, 0, 0, 420, 0, 0; ...
0, 162, 630/L, 0, 468, 630/L; ...
0, -L*39, -126, 0, -L*66, 126 ] * pX';
% global element vector
ele_vec = T' * ele_vec_local;
end
%------------------------------------------------------------------
%% graphical output: plot element
function plot(self, config, scale)
% my nodes
n1 = self.node(1);
n2 = self.node(2);
% coordinates and displacements
[L, T] = self.transform();
x1 = [ n1.ref(1), n2.ref(1) ];
x2 = [ n1.ref(2), n2.ref(2) ];
aa = [ n1.lhs( find(n1.dof == mafe.DofType.Disp1) ), ...
n1.lhs( find(n1.dof == mafe.DofType.Disp2) ), ...
n1.lhs( find(n1.dof == mafe.DofType.Rota3) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Disp1) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Disp2) ), ...
n2.lhs( find(n2.dof == mafe.DofType.Rota3) )];
aa = (T*aa')';
%
function N = interpolate_hermite(z)
[L, T] = self.transform();
N = [ (0.25*(1.-z).^2) .*(2.+z), ...
(0.25*(1.-z).^2) .*( 1.+z) * L/2., ...
(0.25*(1.+z).^2) .*(2.-z), ...
(0.25*(1.+z).^2) .*(-1.+z) * L/2. ];
end
function N = interpolate_linear(z)
N = [ 0.5*(1.-z), 0.5*(1.+z) ];
end
NN = 20;
zz = linspace(-1,+1,NN);
xx1 = interpolate_linear(zz') * x1';
xx2 = interpolate_linear(zz') * x2';
uu = interpolate_linear(zz') * aa([1,4])';
ww = interpolate_hermite(zz') * aa([2,3,5,6])';
ww1 = T(1,1) * uu + T(2,1) * ww;
ww2 = T(1,2) * uu + T(2,2) * ww;
%
% defaults
lcolor = [0.6 0.6 0.6];
mcolor = [0.2 0.2 0.2];
lwidth = 0.5;
mwidth = 10.;
ecolor = [0.4 0.4 0.4];
marker = '.';
%
switch config
case 'deformed'
XX1 = xx1+ww1*scale;
XX2 = xx2+ww2*scale;
lcolor = [43.9, 50.6, 0.0]/100;
mcolor = [43.9, 50.6, 0.0]/100;
falpha = 0.6;
lwidth = 1.0;
mwidth = 20.;
case 'intensity'
XX1 = xx1;
XX2 = xx2;
normalforce = 0.5 * ( -self.internal(1) + self.internal(4) );
if normalforce < 0.0
lcolor = 'b';%[15.3, 31.0, 54.5]/100;
elseif normalforce > 0.0
lcolor = 'r';%[82.7, 28.4, 16.9]/100;
else
lcolor = 'g';%[58.8, 65.5, 7.8]/100;
end
mcolor = [0.0 0.0 0.0 0.0];
falpha = 0.6;
lwidth = 2.0;
mwidth = 0.1;
otherwise
plot@mafe.Li2Element2D(self, config, scale);
return;
end
mafe.patchline(XX1, XX2, ...
'EdgeColor', lcolor, ...
'LineWidth', lwidth, ...
'EdgeAlpha', falpha);
plot(XX1([1 end]), XX2([1 end]), ...
marker, ...
'color', mcolor, ...
'markersize', mwidth);
end
%% print element data
function str = print(self)
str = sprintf('%+4.3e ', [ -self.internal(1), ... % left N
self.internal(2), ... % left Q
-self.internal(3), ... % left M
self.internal(4), ... % right N
-self.internal(5), ... % right Q
self.internal(6) ]); % right M
end
end
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment