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

cleaning & check comments/help in element (support) functions up to ele (issue #1)

parent f2aed5f8
......@@ -10,7 +10,7 @@
https://gitlab.tue.nl/stem/mlfem_nac/tree/v20-q2
* [update] fix code comments throughout (mostly in demos)
* [update] cleaning & check of comments/help in element (support) functions (closes issue #1)
* [update] cleaning & check of comments/help in most 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)
......
function [Me, Ce, Ke, rhse, dum3] = elcd(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse] = elcd(ielem, coord, top, mat)
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
%
......@@ -20,10 +20,11 @@ function [Me, Ce, Ke, rhse, dum3] = elcd(ielem, coord, top, mat, varargin)
% 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
dum3 = [];
dummy = [];
% get the coordinates of the nodes of this element
nodcoord = coord(top(ielem, 1:end-2), :);
......
function [par1, par2] = elcd_i(ielem, ~,top, mat, ichois)
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
......@@ -168,7 +168,7 @@ elseif ichois == 6
elseif ichois == 7
nvar=4; nmode=1;
nvar = 4; nmode = 1;
if ietype == 1, nint = 4; end
if ietype == 2, nint = 9; end
......
function [Me, Ce, Ke, rhse, void] = ele(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse, dummy] = ele(ielem, coord, top, mat, varargin)
function [Me, Ce, Ke, rhse, dummy] = ele(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse, dummy] = ele(ielem, coord, top, mat)
%
% linear elastic (Hookean) element
%
% Parameters:
% Input:
% coord coordinates of all nodes
% top topology array
% ielem element number
% mat material properties (mat.mat)
% E = mat(imat, 1); (Youngs modulus)
% nu = mat(imat, 2); (Poisson ratio)
% axi = mat(imat, 3); (0: plane strain, 1: axi-symmetric, 2: plane stress)
%
% Me : element mass matrix (not used, empty)
% Ce : element damping matrix (not used, empty)
% Ke : element stiffness matrix
% rhse : element right hand side (not calculated, zeros)
% coord : global coordinate array for all nodes
% top : mesh topology array
% ielem : element number
% mat : material properties
% E = mat(imat, 1) (Young's modulus)
% nu = mat(imat, 2) (Poisson's ratio)
% axi = mat(imat, 3) (0: plane strain, 1: axi-symmetric, 2: plane stress)
% eltype = mat(imat, 11) (Element type specification, see ele_i)
%
% void = [] and varargin are for consistenty of element calls
% Output:
% Me element mass matrix
% Ce element damping matrix
% Ke element stiffness matrix
% dummy = [], for consistency of element calls
%
% See also ele_s, ele_i, ele_d, femlin_e
% generate the empty outputs
Me = []; Ce = []; void = [];
Me = []; Ce = []; dummy = [];
% get the number of nodes on this element
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
......
......@@ -97,17 +97,8 @@ str = zeros(nint, 4);
% perform the element integration
for int = 1:nint
if axi == 1 % axi-symmetric
x = n(int, :)*nodcoord(:, 1); % ipts x-coordinate (r)
intw = 2*pi*x*intworg; % multiply Jacobian with 2*pi*r
else
intw = intworg;
end
% get shape function derivatives with respect to global coordinates
[dndxy, detj] = getderiv(dn, nodcoord, int);
intwdetj = intw(int)*detj;
dndxy = getderiv(dn, nodcoord, int);
% compute the strain displacement matrix
B(1, ii+1) = dndxy(:, 1)';
......@@ -118,15 +109,15 @@ for int = 1:nint
B(4, ii+1) = n(int, :)/x;
end
% calculate strain
strain = B*u;
% store strain and sigma
str(int, :) = strain';
sig(int, :) = (H*strain)';
end
% reshape
% reshape strain and sigma
sigma = int2elm_e(ielem, [sig str], coord, top, mat);
%strain=int2elm_e(ielem,str,coord,top,mattypes);
......
......@@ -20,14 +20,17 @@ function [par1, par2] = ele_i(ielem, ~, top, mat, ichois)
%
% See also ele, ele_d, ele_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 triangular-element
par1 = [2 2 2];
......@@ -41,6 +44,7 @@ if ichois == 1
par1 = [2 2 2 2];
end
% sanity check
npar = length(par1);
if npar ~= nlnodes
string_err = ['Number of nodes = ' num2str(nlnodes) ...
......@@ -66,6 +70,7 @@ elseif ichois == 2
elseif ichois == 3
% definition of sequence of nodes for plotting purposes
if ietype == 1
par1 = [1 2; 2 3; 3 1];
......@@ -80,26 +85,26 @@ elseif ichois == 3
end
elseif ichois == 4
% shape function that defines the geometry
if ietype == 1
% linear triangle
gshp = [23 1 2];
elseif ietype == 2
% 6-noded triangle
% quadratic triangle
gshp(1, :) = [26 2 7];
elseif ietype == 3
% 4-noded quad
% linear quad
gshp(1, :) = [14 1 3];
elseif ietype == 4
% 9-noded quad
% quadratic quad
gshp(1, :) = [19 1 3];
elseif ietype == 5
% 4-noded triangle
% extended (4-noded) triangle
gshp(1, :) = [24 1 3];
end
......@@ -108,9 +113,10 @@ elseif ichois == 4
elseif ichois == 5
% defintion of location of dofs (vpos) and their shape functions
if ietype == 1
% 3-noded triangle with linear pressure
% linear triangle
vpos = zeros(2, 3);
vpos(1, :) = 1:2:5;
vpos(2, :) = vpos(1, :) + 1;
......@@ -120,7 +126,7 @@ elseif ichois == 5
vshp(2, :) = [23 2 3];
elseif ietype == 2
% 6-noded triangle
% quadratic triangle
vpos = zeros(3, nlnodes);
vpos(1, :) = 1:2:nlnodes*2;
vpos(2, :) = vpos(1, :) + 1;
......@@ -130,7 +136,7 @@ elseif ichois == 5
vshp(2, :) = [26 2 7];
elseif ietype == 3
% 4-noded quad
% linear quad
vpos = zeros(2, nlnodes);
vpos(1, :) = 1:2:nlnodes*2;
vpos(2, :) = vpos(1, :) + 1;
......@@ -140,7 +146,7 @@ elseif ichois == 5
vshp(2, :) = [14 1 2];
elseif ietype == 4
% 9-noded quad
% quadratic quad
vpos = zeros(2, nlnodes);
vpos(1, :) = 1:2:nlnodes*2;
vpos(2, :) = vpos(1, :) + 1;
......@@ -150,7 +156,7 @@ elseif ichois == 5
vshp(2, :) = [19 1 3];
elseif ietype == 5
% 4-noded triangle with linear pressure
% extended (4-noded) triangle
vpos = zeros(2, 4);
vpos(1, :) = 1:2:7;
vpos(2, :) = vpos(1, :) + 1;
......@@ -165,40 +171,37 @@ elseif ichois == 5
elseif ichois == 6
if ietype == 1
% linear triangle
par1 = [1 0; 0 1; 0 0];
elseif ietype == 2
% quadratic triangle
par1 = [1 0; 1/2 1/2; 0 1; 0 1/2; 0 0; 1/2 0];
elseif ietype == 3
% linear quad
par1 = [-1 -1; 1 -1; 1 1; -1 1;];
elseif ietype == 4
% quadratic quad
par1 = [-1 -1; 0 -1; 1 -1; 1 0; 1 1; 0 1; -1 1; -1 0; 0 0];
elseif ietype == 5
% extended (4-noded) triangle
par1 = [1 0; 0 1; 0 0; 1/3 1/3];
end
elseif ichois == 7
if ietype < 6
nvar = 8; nmode = 1;
if ietype == 1, nint = 3; end
if ietype == 2, nint = 7; end
if ietype == 3, nint = 4; end
if ietype == 4, nint = 9; end
if ietype == 5, nint = 3; end
end
nvar = 8; nmode = 1;
if ietype == 1, nint = 3; end
if ietype == 2, nint = 7; end
if ietype == 3, nint = 4; end
if ietype == 4, nint = 9; end
if ietype == 5, nint = 3; end
par1 = [nint nvar nmode];
else
......
......@@ -43,8 +43,7 @@ end
% initialise detj & intw
detj = []; intw = [];
%
% get integration points
if (intrule == 1 && isempty(intcoord))
[intw, intcoord] = gauss(nint, nsd);
elseif (intrule == 2 && isempty(intcoord))
......@@ -53,14 +52,16 @@ elseif (intrule == 2 && isempty(intcoord))
elseif (intrule == 3 && isempty(intcoord))
[intw, intcoord] = gaussl(nint, nsd);
end
[n, dn] = shape(intcoord, shfunction);
if ~isempty(idet),
% [nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
% get Jacobian
if ~isempty(idet)
% nodal coordinates
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
nodcoord = coord(top(ielem, 1:nlnodes), :);
% number of integration points
nint = length(intcoord(:, 1));
% allocate
detj = zeros(1, nint);
for int = 1:nint
[~, detj(int)] = getderiv(dn, nodcoord, int);
......
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