Commit 67ca846a authored by Turnhout, M.C. van's avatar Turnhout, M.C. van
Browse files

Merge branch 'testing'

parents 20f1decf 10a73e2f
.gitattributes export-ignore
.gitignore export-ignore
/manual_src export-ignore
mlfem_nac.pdf -diff
......@@ -39,6 +39,7 @@ Thumbs.db
# Linux git ignore
.*
!.gitignore
!.gitattributes
*~
......
# 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 parameters in mat.mat (there is a very good explanation for this, but still...)
# Version 21
https://gitlab.tue.nl/stem/mlfem_nac/tree/v21
* [update] fix code comments throughout (mostly in demos)
* [update] cleaning & check of comments/help in linear element (support) functions (issue #1)
* [update] update names of (element) matrices to reflect terminology in book (NAC I) and lecture notes (NAC II)
* [feature] feedback on which elm1d_c.m or elcd_a.m function is used (fem1d, femlin_cd)
* [feature] calculate flux cdudx in fem1dcd and femlin_cd (closes issue #8)
* [update] update crmesh chapter in manual for more clarity (thanks to Roderick for the feedback!)
* [update] fix some minor issues in manual (typos, typesetting, etc.) and update to v21 (new global matrix names reflect terminology in book)
* [update] update engine chapter in manual with global matrices, equations, and engine output for linear elements
* [update] add references to book exercises and equations to index of manual, and update index lay-out (because I can)
# Version 20-q1
https://gitlab.tue.nl/stem/mlfem_nac/tree/v20-q1
......@@ -11,19 +33,14 @@ https://gitlab.tue.nl/stem/mlfem_nac/tree/v20-q1
* [feature] allow 'options' for sysbuild and 'newtoptions' for sysnewt to be set in the pre-processing phase
* [update] thorough update of the manual
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...)
# 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.
......@@ -33,7 +50,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.
......@@ -41,7 +58,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'.
......
% %%%%%%%%%%%%%% demo of femlin_e %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %
% This script analyses the planar deformation of a thick cylinder at %
% small deformations, using an axi-symmetric model and an %
% incompressible Neo-Hookean material law: %
% small deformations, using an axi-symmetric model and a linear %
% Hookean material law. %
% %
% sigma = -pI + G(B - I) %
% The demo uses femlin_e as engine with the element ele. %
% %
% The demo uses femlin_e as engine with the element elup. %
% %
% See also demo_2ds_axi_cyl for a (compressible) Hookean implementation %
% See also demo_2ds_axi_cylup for an incompressible Neo-Hookean %
% implementation %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, close all %#ok<CLSCR>
......
......@@ -51,10 +51,10 @@ itype = 1; norder = 2; ietype = 4; % quadratic quad
crmesh(curves, subarea, points, norder, itype);
%%%%%%%%%%%%%%%%%% DEFINE THE MATERIAL PROPERTIES %%%%%%%%%%%%%%%%%%%%%%%
E = 1e6; % Youngs modulus
nu = 0.3; % Poisson ratio
mat.mat(1) = E; % E: Youngs modulus
mat.mat(2) = nu; % nu: Poisson ratio
E = 1e6; % Youngs modulus [Pa]
nu = 0.3; % Poisson ratio [-]
mat.mat(1) = 1e6; % E: Youngs modulus
mat.mat(2) = 0.3; % nu: Poisson ratio
mat.mat(3) = 2; % axi: 0-plane strain, 1-axisymmetric, 2-plane stress
mat.mat(11) = ietype;
% define element type(s)
......@@ -73,14 +73,14 @@ bndcon = addbndc(bndcon, coord, top, mat, usercurves, crv, dof, string, ...
iplot);
% apply a downward force of 5e-5 on the upper right corner: points(3, :)
F = -5e-4;
F = -5e-3;
nodfrc = [userpoints(3) 2 F];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create and solve equations %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Solve the problem using the 2D fem program FEMNL %%%%%%%%%%
%%%%%%%%%%%%% Solve the problem using the 2D FEM program FEMLIN_E %%%%%%%
% do the finite element computation
t0 = cputime;
......@@ -96,11 +96,11 @@ disp(['CPU-time consumption in femlin_e : ' num2str(cputime-t0)]);
% plot undeformed and deformed mesh (scale deformations with 1e4 to
% properly show 'em)
figure
pldisp(sol*1e4, coord, top, dest, mat, 1, 1, 2)
pldisp(sol, coord, top, dest, mat, 1e4)
% plot element stresses: s_11
figure
plelmdat(sigma_elm,coord,top,mat,1);
plelmdat(sigma_elm, coord, top, mat, 1);
% calculate displacements for linear beam theory
......@@ -109,9 +109,9 @@ I = B*H^3/12;
u_lin = -((F*L^3)/6 - (F*L^3)/2)/(E*I);
fprintf('Displacement u_lin according to linear theory: %1.3e\n', u_lin);
% get actual displacements (from the upper left corner)
u_a = sol(dest(userpoints(3), 2));
fprintf('Displacement u_a in simulation: %1.3e\n', u_a);
u_s = sol(dest(userpoints(3), 2));
fprintf('Displacement u_sim in simulation: %1.3e\n', u_s);
% ratio?
fprintf('u_a/u_lin: %1.3f\n', u_a/u_lin);
fprintf('u_sim/u_lin: %1.3f\n', u_s/u_lin);
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
......@@ -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
......
function [qme, qde, qke, rhse, dum3] = ...
elcd(ielem, coord, top, mat, pos, sol, soln, dum1, dum2)
function [Me, Ce, Ke, rhse, dummy] = elcd(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse, dummy] = elcd(ielem, coord, top, mat)
%
% convection diffusion element element
%
% du/dt + v*du/dx = d/dx (c du/dx)
%
% Parameters:
% Input:
% coord coordinates of all nodes
% top topology array
% ielem element number
% mat material properties (mat.mat)
% c = mat(imat, 1); (diffusion coefficient)
% vx = mat(imat, 2); (factor for convective velocity in x-direction)
% vy = mat(imat, 3); (factor for convective velocity in y-direction)
% axi = mat(imat, 4); (0: plane strain, 1: axi-symmetric, 2: plane stress)
%
% qme : element mass matrix
% qde : element damping matrix
% qe : element stiffness matrix
% rhse : element right hand side
% coord : coordinates of all nodes
% top : topology array
% ielem : element number
% mat : material properties
% c = mat(imat,1); (diffusion coefficient)
% vx = mat(imat,2); (factor for convective velocity in x-direction)
% vy = mat(imat,3); (factor for convective velocity in y-direction)
% axi = mat(imat,4); (0: plane strain, 1: axi-symmetric, 2: plane stress)
qme = []; qde = []; dum3 = [];
% Output:
% Me element mass matrix
% Ce element damping matrix
% Ke element stiffness matrix
% rhse element right hand side
% dummy = [], for consistency of element calls
%
% See also femlin_cd, elcd_i, elcd_s, elcd_a
% get the number of nodes on this element
[nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ...
sysinfo(size(coord), size(top));
nlnodes = length(nonzeros(top(ielem, 1:maxnlnodes)));
dummy = [];
% get the coordinates of the nodes of this element
nodcoord = coord(top(ielem, 1:nlnodes), :);
nodcoord = coord(top(ielem, 1:end-2), :);
% pick up some material parameter
imat = top(ielem, nlnodes + 1);
imat = top(ielem, end-1);
c = mat(imat, 1);
vx = mat(imat, 2);
vy = mat(imat, 3);
axi = mat(imat, 4);
% element topology
etop = top(ielem, :);
% position of the local degrees of freedom in the global solution array
[vpos, vshp] = elcd_i( 1, [], etop, mat, 5);
[n, dn, intworg, intcoord] = elcd_s([], [], etop, mat, vshp(1, :));
[vpos, vshp] = elcd_i( 1, [], etop, mat, 5);
[n, dn, intworg] = elcd_s([], [], etop, mat, vshp(1, :));
% number of dofs for the temperature
% number of dofs for the element
nedof = length(nonzeros(vpos(1, :)));
% number of integration points
nint = length(intworg);
% extract the location of the displacement degrees of freedom
it = nonzeros(vpos(1, :));
% initialize qe and rhs
qke = zeros(nedof, nedof);
qde = qke;
qme = qke;
% initialize Ke, Ce, Me, and rhse
Ke = zeros(nedof, nedof);
Ce = Ke;
Me = Ke;
rhse = zeros(nedof, 1);
% perform the element integration
for int = 1:nint
% compute x (=r)-coordinate for axi-symmetric case
x = n(int, :)*nodcoord(:, 1);
y = n(int, :)*nodcoord(:, 2);
if axi > 0
intw = 2*pi*x*intworg;
else
intw = intworg;
end
% get the convective velocity
[ax, ay] = elcd_a(x, y, vx, vy);
% get shape function derivatives with respect to global coordinates
[dndxy, detj] = getderiv(dn, nodcoord, int);
intwdetj = intw(int)*detj;
% diffusion part
qke = qke + c*(dndxy(:, 1)*dndxy(:, 1)' + ...
dndxy(:, 2)*dndxy(:, 2)')*intwdetj;
% convection part
qde = qde + (ax*n(int, :)'*dndxy(:, 1)' + ...
ay*n(int, :)'*dndxy(:, 2)')*intwdetj;
% mass matrix
qme = qme + n(int, :)'*n(int, :)*intwdetj;
for int = 1:length(intworg)
% compute x (=r)-coordinate for axi-symmetric case
x = n(int, :)*nodcoord(:, 1);
y = n(int, :)*nodcoord(:, 2);
if axi > 0
intw = 2*pi*x*intworg;
else
intw = intworg;
end
% get the convective velocity
[ax, ay] = elcd_a(x, y, vx, vy);
% get shape function derivatives with respect to global coordinates
[dndxy, detj] = getderiv(dn, nodcoord, int);
intwdetj = intw(int)*detj;
% diffusion part
Ke = Ke + c*(dndxy(:, 1)*dndxy(:, 1)' + ...
dndxy(:, 2)*dndxy(:, 2)')*intwdetj;
% convection part
Ce = Ce + (ax*n(int, :)'*dndxy(:, 1)' + ...
ay*n(int, :)'*dndxy(:, 2)')*intwdetj;
% mass matrix
Me = Me + n(int, :)'*n(int, :)*intwdetj;
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
\ No newline at end of file
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
......@@ -4,7 +4,8 @@ function [ax, ay] = elcd_a(x, y, vx, vy)
% function to create a flow field as a function of x and y
% for the solution of the convection diffusion equation
%
% See also demo_2d_conv_diff_planar, demo_2d_conv_diffusion
% See also elcd, femlin_cd, demo_2d_conv_diff_planar,
% demo_2d_conv_diffusion
ax = vx;
ay = vy;
......
function [par1,par2]=elcd_i(ielem,coord,top,mat,ichois);
%
% description:
%
% This class of elements solves either the incompressible Stokes problem
% or a (in) compressible elasticity problem
%
% ietype:
%
% 1 - linear triangle
% 2 - quadratic triangle
% 3 - bilinear quad
% 4 - biquadratic quad
%
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
imat = top(ielem,iimat);
ietype = mat(imat,11);
nlnodes = length(nonzeros(top(ielem,1:maxnlnodes)));
par1=[]; par2=[];
if (ichois==1),
if ietype==1,
par1=[1 1 1]; % Linear triangular-element
elseif ietype==2,
par1=[1 1 1 1 1 1]; % Quadratic element
elseif ietype==3,
par1=[1 1 1 1]; % Bilinear element
elseif ietype==4,
par1=[1 1 1 1 1 1 1 1 1]; % Biquadratic element
end
npar=length(par1);
if npar~=nlnodes,
string_err = ['Number of nodes = ' num2str(nlnodes) ...
' for element type ' num2str(ietype) ...
' in ielem ' num2str(ielem) ', should be '];
error([string_err num2str(npar)]);
end
elseif (ichois==2),
if (ietype==1),
par1=2*ones(1,nlnodes);
elseif ietype==2,
par1=2*ones(1,nlnodes);
elseif ietype==3,
par1=2*ones(1,nlnodes);
elseif ietype==4,
par1=2*ones(1,nlnodes);
end
elseif (ichois==3),
if (ietype==1),
par1=[1 2; 2 3; 3 1];
elseif (ietype==2),
par1=[1 2; 2 3; 3 4; 4 5; 5 6; 6 1];
elseif (ietype==3),
par1=[1 2; 2 3; 3 4; 4 1];
elseif (ietype==4),
par1=[1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];
end
elseif (ichois==4),
if (ietype==1),
% linear triangle
gshp = [23 2 4];
elseif (ietype==2),
% 6-noded triangle
gshp(1,:) = [26 2 7];
elseif (ietype==3),
% 4-noded quad
gshp(1,:) = [14 1 2];
elseif (ietype==4),
% 9-noded quad
gshp(1,:) = [19 1 3];
end
par1=gshp;
elseif (ichois==5),
if (ietype==1),
% 3-noded triangle with linear pressure
vpos = zeros(1,3);
vpos(1,:) = 1:3;
vshp = zeros(1,3);
vshp(1,:) = [23 2 4];
elseif (ietype==2),
% 6-noded triangle
vpos = zeros(1,nlnodes);
vpos(1,:) = 1:nlnodes;
vshp = zeros(1,3);
vshp(1,:) = [26 2 7];
elseif (ietype==3),
% 4-noded quad
vpos = zeros(1,nlnodes);
vpos(1,:) = 1:nlnodes;
vshp = zeros(1,3);
vshp(1,:) = [14 1 2];
elseif (ietype==4),
% 9-noded quad
vpos = zeros(1,nlnodes);
vpos(1,:) = 1:nlnodes;
vshp = zeros(1,3);
vshp(1,:) = [19 1 3];
end
par1=vpos; par2=vshp;
elseif ichois==6,
if ietype==1,
par1=[1 0; 0 1; 0 0;];
elseif ietype==2,
par1=[1 0; 1/2 1/2; 0 1; 0 1/2; 0 0; 1/2 0];
elseif ietype==3,
par1=[-1 -1; 1 -1; 1 1; -1 1;];
elseif ietype==4,
par1=[-1 -1; 0 -1; 1 -1; 1 0; 1 1; 0 1; -1 1; -1 0; 0 0];
function [par1, par2] = elcd_i(ielem, ~, top, mat, ichois)
% [par1, par2] = elcd_i(ielem, ~, top, mat, ichois)
%
% Element information for elcd, output is determined by ichois
%
% input:
% ielem : element number
% ~ : not used (coord, nodal coordinates)
% top : mesh topology array
% mat : material properties (mat.mat), mat(iimat, 11) = ietype
%
% output:
% ichois = 1 par1: number of degrees of freedom for each node
% par2 = [];
% ichois = 3 par1: sequence of nodes for plotting purposes
% par2 = [];
% ichois = 4 par1: shape function that defines the geometry
% par2 = [];
% ichois = 5 par1: definition of location of dofs (vpos)
% par2: definition of their shape functions (vshp)
%
% See also elcd, elcd_s
% get some element info
imat = top(ielem, end-1);
ietype = mat(imat, 11);
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
% initialise
par1 = []; par2 = [];
if ichois == 1
% definition of number of degrees of freedom for each node
if ietype == 1
% linear triangle
par1 = [1 1 1];
elseif ietype == 2
% quadratic triangle
par1 = [1 1 1 1 1 1];
elseif ietype == 3
% linear quad
par1 = [1 1 1 1];
elseif ietype == 4
% quadratic quad
par1 = [1 1 1 1 1 1 1 1 1];
end
% sanity check
npar = length(par1);
if npar ~= nlnodes
string_err = ['Number of nodes = ' num2str(nlnodes) ...
' for element type ' num2str(ietype) ...
' in ielem ' num2str(ielem) ', should be '];
error([string_err num2str(npar)]);
end
elseif ichois == 2
if ietype == 1
% linear triangle
par1 = 2*ones(1, nlnodes);
elseif ietype == 2
% quadratic triangle
par1 = 2*ones(1, nlnodes);
elseif ietype == 3
% linear quad
par1 = 2*ones(1, nlnodes);
elseif ietype == 4
% quadratic quad
par1 = 2*ones(1, nlnodes);
end
elseif ichois == 3
% definition of sequence of nodes for plotting purposes
if ietype == 1
% linear triangle
par1 = [1 2; 2 3; 3 1];
elseif ietype == 2
% quadratic triangle
par1 = [1 2; 2 3; 3 4; 4 5; 5 6; 6 1];
elseif ietype == 3
% linear quad
par1 = [1 2; 2 3; 3 4; 4 1];
elseif ietype == 4
% quadratic quad
par1 = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 1];
end
elseif ichois == 4
% shape function that defines the geometry
if ietype == 1
% linear triangle
gshp = [23 2 4];
elseif ietype == 2
% quadratic triangle
gshp(1, :) = [26 2 7];
elseif ietype == 3
% linear quad
gshp(1, :) = [14 1 2];
elseif ietype == 4
% quadratic quad
gshp(1, :) = [19 1 3];
end
par1 = gshp;
elseif ichois == 5
% definition of location of dofs (vpos) and their shape functions
if ietype == 1
% linear triangle
vpos = zeros(1, 3);
vpos(1, :) = 1:3;
vshp = zeros(1, 3);
vshp(1, :) = [23 2 4];
elseif ietype == 2
% quadratic triangle
vpos = zeros(1, nlnodes);
vpos(1, :) = 1:nlnodes;
end
elseif ichois==7,