Commit 7f8003c2 authored by Turnhout, M.C. van's avatar Turnhout, M.C. van
Browse files

rename global matrices for consistency with book

parent 09b58ee5
......@@ -3,12 +3,9 @@
https://gitlab.tue.nl/stem/mlfem_nac/tree/v20-q2
* [update] fix code comments throughout (mostly in demos)
* [update] update names of (element) matrices to reflect terminology in book (NAC I) and lecture notes (NAC II)
Known issues (see also https://gitlab.tue.nl/stem/mlfem_nac/-/issues ):
* [**severe**] the implementation of the **updated/total Lagrange** method is **broken** (there are no error messages, but there are no sensible results either).
* [mild] time-dependent pulsatile flow is hard-coded in the engine femnlt
* [minor] femlt puts some, but all, istat parametters in mat.mat (there is a very good explanation for this, but still...)
Known issues: see *Version 20-q1*
# Version 20-q1
......@@ -27,16 +24,16 @@ https://gitlab.tue.nl/stem/mlfem_nac/tree/v20-q1
Known issues (see also https://gitlab.tue.nl/stem/mlfem_nac/-/issues ):
* [**severe**] the implementation of the **updated/total Lagrange** method is **broken** (there are no error messages, but there are no sensible results either).
* [mild] time-dependent pulsatile flow is hard-coded in the engine femnlt
* [minor] femlt puts some, but all, istat parametters in mat.mat (there is a very good explanation for this, but still...)
* [minor] femlt puts some, but all, istat parameters in mat.mat (there is a very good explanation for this, but still...)
# Version 19, The legacy of Cees Oomens
https://gitlab.tue.nl/stem/mlfem_nac/tree/v19
In the last year(s) before his pension, Cees Ooomens took it upon himself to go debug, stream line and improve the entire code once again for our educational purposes.
In the last year(s) before his pension, Cees Oomens took it upon himself to go debug, stream line and improve the entire code once again for our educational purposes.
This is the state of the code at Cees' pension date (June 25th, 2020).
This is the state of the code at Cees' pension date (June 25th, 2020).
Includes the source for the documentation.
......@@ -46,7 +43,7 @@ https://gitlab.tue.nl/stem/mlfem_nac/tree/v16
The original scripts of version 01 were adapted for educational use over the years.
Most of those adaptations were made by Cees Oomens and a team of student assistants. Sadly, history has not been so kind as to preserve the identities of those students for prosperity.
Most of those adaptations were made by Cees Oomens and a team of student assistants. Sadly, history has not been so kind as to preserve the identities of those students for prosperity.
This is the state of the code in the academic year 2016/2017.
......@@ -54,7 +51,7 @@ This is the state of the code in the academic year 2016/2017.
https://gitlab.tue.nl/stem/mlfem_nac/tree/v01
These are the original scripts that Frank Baaijens wrote for his personal use at Philips.
These are the original scripts that Frank Baaijens wrote for his personal use at Philips.
He took the scripts with him when he joined Eindhoven University of Technology.
This appears to have been around 1994, making this 'version 00.94' which rounds off to 'version 01'.
......
......@@ -126,7 +126,7 @@ femlin_e
% extract the total (reaction) shear force on the top plate (crv = 3)
% first get the reaction forces
f = q*sol; % q is the (global) stiffness matrix
f = K*sol; % K is the (global) stiffness matrix
% then select the nodes of the top plate
nodes = nonzeros(usercurves(3, :));
% get the locations of the first degree of freedom (x) in sol
......
% pick up some global variables
% nnodes: number of nodes
% nsd : space dimension (1 for 1d problems, 2 for 2d problems, etc)
% nelem : number of elements
% maxnlnodes: maximum number of nodes in an element
[nnodes, nsd, nelem, maxnlnodes] = sysinfo(size(coord), size(top));
if ~exist('nodfrc', 'var'), nodfrc = []; end
% compute the pos and dest array
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
% global K and f
disp('* start of matrix-building process ...')
% ndof: the total number of degrees of freedom
ndof = max(dest(:));
% initialize q and rhs;
% q: the global stiffness matrix (spalloc allocated space for a sparse matrix)
% K: the global stiffness matrix (spalloc allocated space for a sparse matrix)
% rhs: the global right-hand-side matrix
ldof = length(pos(1, :));
q = spalloc(ndof, ndof, sum(ldof)^2);
K = spalloc(ndof, ndof, sum(ldof)^2);
rhs = zeros(ndof,1);
% assemble q and rhs
......@@ -26,17 +26,18 @@ for ielem = 1:nelem
% assemble the global stiffness matrix and the right-hand-side array
ii = pos(ielem, :);
q(ii, ii) = q(ii, ii) + Ke;
K(ii, ii) = K(ii, ii) + Ke;
rhs(ii) = rhs(ii) + rhse;
end
% solve the system of equations
if ~exist('nodfrc', 'var'), nodfrc = []; end
% solvestat takes care of boundary conditions and partitioning of system of equations
sol = solvestat(q, rhs, bndcon, nodfrc, dest);
% sum & partition and solve
disp('* solving the equations ...')
% solvestat takes care of boundary conditions and partitioning of system of equations
sol = solvestat(K, rhs, bndcon, nodfrc, dest);
% compute the rhs including the reaction forces
rhs = q*sol;
rhs = K*sol;
%* compute the stresses
sigma = zeros(nelem, maxnlnodes);
......
......@@ -4,29 +4,6 @@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pick up some global variables
[nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ...
sysinfo(size(coord), size(top));
% define position array
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
% ndof: the total number of degrees of freedom
ndof = max(dest(:));
%d efine number of unknows for each node
ldof = ones(1, maxnlnodes);
% initialize q and rhs;
% qk: the global stiffness matrix
% qm: the global mass matrix
% qd: the global damping matrix
% rhs: the global right-hand-side matrix
qk = spalloc(ndof, ndof, sum(ldof)^2);
qm = qk;
qd = qk;
rhs = zeros(ndof, 1);
if ~exist('sol', 'var')
% initialise sol
ndof = max(dest(:));
......@@ -38,6 +15,29 @@ if ~exist('sol', 'var')
end
if ~exist('nodfrc', 'var'), nodfrc = []; end
% pos and dest
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
% global M, C, K and f
disp('* start of matrix-building process ...')
% ndof: the total number of degrees of freedom
ndof = max(dest(:));
% define number of unknows for each node
ldof = ones(1, size(top, 2)-2);
% initialize q and rhs;
% K: the global stiffness matrix
% M: the global mass matrix
% C: the global damping matrix
% rhs: the global right-hand-side matrix
K = spalloc(ndof, ndof, sum(ldof)^2);
M = K;
C = K;
rhs = zeros(ndof, 1);
% assemble q and rhs
for ielem = 1:nelem
% compute the element stiffness matrix and right-hand-side column
......@@ -45,26 +45,29 @@ for ielem = 1:nelem
% assemble the global stiffness matrix and the right-hand-side colomn
ii = pos(ielem, :);
qk(ii, ii) = qk(ii, ii) + Ke;
qm(ii, ii) = qm(ii, ii) + Me;
qd(ii, ii) = qd(ii, ii) + Ce;
K(ii, ii) = K(ii, ii) + Ke;
M(ii, ii) = M(ii, ii) + Me;
C(ii, ii) = C(ii, ii) + Ce;
rhs(ii) = rhs(ii) + rhse;
end
% sum & partition and solve
disp('* solving the equations ...')
if istat == 1
% solve the system of equations, stationairy problem
q = qk + qd;
q = K + C;
sol = solvestat(q, rhs, bndcon, [], dest);
elseif istat == 2
% unsteady solution
q = qm/dt + theta*qd + theta*qk;
q = M/dt + theta*C + theta*K;
soln = sol(:, 1);
for itime = 1:ntime
fprintf('** Start of time step %i\n', itime);
rhsn = rhs + (qm/dt - (1 - theta)*(qd + qk))*soln;
rhsn = rhs + (M/dt - (1 - theta)*(C + K))*soln;
sol(:, itime + 1) = solvestat(q, rhsn, bndcon, [], dest);
soln = sol(:, itime+1);
end
......
% obtain some global information
% nnodes : total number of nodes
% nsd : spatial dimension (1, 2, or 3D)
% nelem : number of elements
% maxnlnodes : maximum number of nodes of each element
% iimat : column of the top array that contains the material
% property identifiers
% iitype : column of the top array that contains the element
% type identifiers
[nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ...
sysinfo(size(coord), size(top));
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
options = [1 1 1 1 1];
% options: used in sysbuild to build systems of equations
% options(1): build mass matrix
% (2): build damping matrix
% (3): build stiffness matrix
% (4): build rhs
% (5): pass only one element topology etc to element function
if ~exist('options', 'var')
options = [1 1 1 1 1];
% options: used in sysbuild to build systems of equations
% options(1): build mass matrix
% (2): build damping matrix
% (3): build stiffness matrix
% (4): build rhs
% (5): pass only one element topology etc to element function
end
nodaldata = [];
elemdata_in = [];
......@@ -30,7 +16,7 @@ bndcon_org = bndcon;
if ~exist('sol', 'var')
% initialise sol
ndof = max(dest(:));
if istat==1 % steady state solution
if istat == 1 % steady state solution
sol = zeros(ndof, 1);
elseif istat == 2 % unsteady solution
sol = zeros(ndof, ntime + 1);
......@@ -38,24 +24,31 @@ if ~exist('sol', 'var')
end
if ~exist('nodfrc', 'var'), nodfrc = []; end
% pos and dest
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
% global M, C, K and f
disp('* start of matrix-building process ...')
[qm, qd, qk, rhs] = sysbuild(coord, top, mat, pos, dest, options, ...
[M, C, K, rhs] = sysbuild(coord, top, mat, pos, dest, options, ...
[], [], [], []);
% sum & partition and solve
disp('* solving the equations ...')
if istat == 1
% solve the system of equations, steady state problem
q = qk + qd;
q = K + C;
sol = solvestat(q, rhs, bndcon, [], dest);
elseif istat == 2
% unsteady solution
q = qm/dt + theta*qd + theta*qk;
q = M/dt + theta*C + theta*K;
soln = sol(:, 1);
for itime = 1:ntime
fprintf('** Start of time step %i\n', itime);
rhsn = rhs + (qm/dt - (1 - theta)*(qd + qk))*soln;
rhsn = rhs + (M/dt - (1 - theta)*(C + K))*soln;
sol(:, itime + 1) = solvestat(q, rhsn, bndcon, [], dest);
soln = sol(:, itime + 1);
end
......
......@@ -5,11 +5,6 @@
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ...
sysinfo(size(coord), size(top));
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
if ~exist('options', 'var')
options = [1 1 1 1 0];
......@@ -21,6 +16,10 @@ if ~exist('options', 'var')
% (5): pass only one element topology etc to element function
end
% pos and dest
disp('* create equation numbers ...');
[pos, dest] = equatnr(coord, top, mat);
if ~exist('sol', 'var')
% initialise sol
ndof = max(dest(:));
......@@ -28,12 +27,15 @@ if ~exist('sol', 'var')
end
if ~exist('nodfrc', 'var'), nodfrc = []; end
% global K and f
disp('* start of matrix-building process ...')
[qm, qc, q, rhs] = sysbuild(coord, top, mat, pos, dest, options, sol, ...
[~, ~, K, rhs] = sysbuild(coord, top, mat, pos, dest, options, sol, ...
[], [], []);
% solution
sol = solvestat(q, rhs, bndcon, nodfrc, dest);
% partition and solve
disp('* solving the equations ...')
sol = solvestat(K, rhs, bndcon, nodfrc, dest);
% compute derived quantities
[sigma_elm, sigma_int, elweight] = sysdrv(coord, top, mat, pos, dest, ...
sol, [], [], []);
......
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