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

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

parent d76ced51
function [Me, Ce, Ke, rhse, dum3] = ...
eler(ielem, coord, top, mat, pos, sol, soln, dum1, dum2)
function [Me, Ce, Ke, rhse, dummy] = eler(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse, dummy] = eler(ielem, coord, top, mat)
%
% linear elastic element, reduced integration
% linear elastic (Hookean) element, reduced integration
%
% 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)
%
% Me : element mass matrix (not used)
% Ce : element damping matrix (not used)
% Ke : element stiffness matrix
% rhse : element right hand side
% coord : coordinates of all nodes
% top : 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)
% fx = mat(imat,4); (Volume force in x-direction)
% fy = mat(imat,5); (Volume force in y-direction)
% eltype = mat(imat,11)(Element type specification, see eler_i)
% Output:
% Me element mass matrix
% Ce element damping matrix
% Ke element stiffness matrix
% dummy = [], for consistency of element calls
%
% Remark
% mat.types : 'eler;
% See also eler_s, eler_i, eler_d, femlin_e
Me = []; Ce = []; dum3 = [];
% generate the empty outputs
Me = []; Ce = []; dummy = [];
% 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)));
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
% get the coordinates of the nodes of this element
nodcoord = coord(top(ielem, 1:nlnodes), :);
% pick up some material parameter
imat = top(ielem, nlnodes+1);
e = mat(imat, 1);
imat = top(ielem, nlnodes + 1);
E = mat(imat, 1);
nu = mat(imat, 2);
axi = mat(imat, 3);
......@@ -43,57 +39,51 @@ axi = mat(imat, 3);
etop = top(ielem, :);
% position of the local degrees of freedom in the global solution array
[vpos, vshp] = eler_i( 1, [], etop, mat, 5);
[n, dn, intworg, intcoord] = eler_s([], [], etop, mat, vshp(1, :));
[vpos, vshp] = eler_i( 1, [], etop, mat, 5);
[n, dn, intworg] = eler_s([], [], etop, mat, vshp(1, :));
% number of dofs for the displacements
nedofu = length(nonzeros(vpos(1:2, :)));
% number of integration points
nint = length(intworg);
% extract the location of the displacement degrees of freedom
iu = zeros(1,nedofu);
iu(1:2:nedofu) = nonzeros(vpos(1, :));
iu(2:2:nedofu) = nonzeros(vpos(2, :));
% initialize Ke and rhs
Ke = zeros(nedofu, nedofu);
rhse = zeros(nedofu, 1);
% initialize strain-displacement matrix
B = zeros(4, nedofu);
ii = 0:2:nedofu-2;
% compute shear modulus and compression modulus from E and Nu
g = e/(2*(1+nu));
kappa = e/(3*(1-2*nu));
G = E/(2*(1+nu));
kappa = E/(3*(1-2*nu));
if axi < 2 % plane strain or axi-symmetric
% D matrix currently assumes plane strain
dd = 1/3*g*[
4 -2 0 -2;
-2 4 0 -2;
0 0 3 0;
-2 -2 0 4];
% sigma = K*tr(epsilon)*I + 2*G*epsilon^d
% epsilon = [xx yy 2xy zz]
% deviatoric part: full integration
Hd = G/3*[
4 -2 0 -2;
-2 4 0 -2;
0 0 3 0;
-2 -2 0 4];
dv = kappa*[
% volumetric part: reduced integration
Hv = kappa*[
1 1 0 1;
1 1 0 1;
0 0 0 0;
1 1 0 1];
1 1 0 1];
end
% initialize strain-displacement matrix
b = zeros(4, nedofu);
ii = 0:2:nedofu-2;
% perform the element integration
% full integration on deviatoric part
for int = 1:nint
for int = 1:length(intworg)
% compute x (=r)-coordinate for axi-symmetric case
x = n(int, :)*nodcoord(:, 1);
if axi == 1 % axi-symmetric
intw = 2*pi*x*intworg;
x = n(int, :)*nodcoord(:, 1); % ipts x-coordinate (r)
intw = 2*pi*x*intworg; % multiply Jacobian with 2*pi*r
else
intw = intworg;
end
......@@ -104,30 +94,28 @@ for int = 1:nint
intwdetj = intw(int)*detj;
% compute the strain displacement matrix
b(1, ii+1) = dndxy(:, 1)';
b(2, ii+2) = dndxy(:, 2)';
b(3, ii+1) = dndxy(:, 2)';
b(3, ii+2) = dndxy(:, 1)';
if axi==1 % axi-symmetric
b(4, ii+1) = n(int, :)/x;
B(1, ii+1) = dndxy(:, 1)';
B(2, ii+2) = dndxy(:, 2)';
B(3, ii+1) = dndxy(:, 2)';
B(3, ii+2) = dndxy(:, 1)';
if axi == 1 % axi-symmetric
B(4, ii+1) = n(int, :)/x;
end
Ke = Ke + b'*dd*b*intwdetj;
% add deviatoric part to stiffness matrix
Ke = Ke + B'*Hd*B*intwdetj;
end
% reduced integration step
[n, dn, intworg, intcoord] = eler_s([], [], etop, mat, vshp(3, :));
% number of integration points
nint = length(intworg);
for int = 1:nint
% compute x (=r)-coordinate for axi-symmetric case
x = n(int, :)*nodcoord(:, 1);
% reduced integration step on volumetric part
[n, dn, intworg] = eler_s([], [], etop, mat, vshp(3, :));
B = zeros(4, nedofu);
for int = 1:length(intworg)
if axi == 1 % axi-symmetric
intw = 2*pi*x*intworg;
x = n(int, :)*nodcoord(:, 1); % ipts x-coordinate (r)
intw = 2*pi*x*intworg; % multiply Jacobian with 2*pi*r
else
intw = intworg;
end
......@@ -138,15 +126,16 @@ for int = 1:nint
intwdetj = intw(int)*detj;
% compute the strain displacement matrix
b(1, ii+1) = dndxy(:, 1)';
b(2, ii+2) = dndxy(:, 2)';
b(3, ii+1) = dndxy(:, 2)';
b(3, ii+2) = dndxy(:, 1)';
B(1, ii+1) = dndxy(:, 1)';
B(2, ii+2) = dndxy(:, 2)';
B(3, ii+1) = dndxy(:, 2)';
B(3, ii+2) = dndxy(:, 1)';
if axi == 1 % axi-symmetric
b(4, ii+1) = n(int, :)/x;
B(4, ii+1) = n(int, :)/x;
end
Ke = Ke + b'*dv*b*intwdetj;
% add volumetric part to stiffness matrix
Ke = Ke + B'*Hv*B*intwdetj;
end
......
function [sigma,sig,elweight] = ...
eler_d(ielem,coord,top,mat,pos,sol,nodaldata,elemdata_in,elemparam);
function [sigma, sig, elweight] = eler_d(ielem, coord, top, mat, pos, ...
sol, varargin)
% [sigma, sig, elweight] = eler_d(ielem, coord, top, mat, pos, ...
% sol, varargin)
%
% linear elastic element
% Stresses c du/dx for ele
%
% Parameters:
% input:
%
% coord : coordinates of all nodes
% top : topology array
% ielem : element number
% coord : global node coordinates array
% top : global topology array
% 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)
% fx = mat(imat,4); (Volume force in x-direction)
% fy = mat(imat,5); (Volume force in y-direction)
% sigma : stresses
% storage:
% sigma(ielem,:)=[sigma_xx sigma_yy sigma_xy sigma_tt .... sigma_xx sigma_yy sigma_xy sigma_tt];
% node_1 node_n
% function [sigma]=ele_d(ielem,coord,top,mat,pos,sol,sigma);
% 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)
% pos : global pos array
% sol : array with calculated displacements (solutions)
%
% output:
% sigma(ielem, :) = [sigma_xx sigma_yy sigma_xy sigma_tt .... sigma_xx sigma_yy sigma_xy sigma_tt];
% node_1 node_n
%
% See also eler, femlin_e
% 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)));
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
% get the coordinates of the nodes of this element
nodcoord=coord(top(ielem,1:nlnodes),:);
nodcoord = coord(top(ielem, 1:nlnodes), :);
mattypes=mat;
mat=mat.mat;
% pick up some material parameters
imat = top(ielem, end-1);
E = mat.mat(imat, 1);
nu = mat.mat(imat, 2);
axi = mat.mat(imat, 3);
% pick up some material parameter
imat=top(ielem,iimat);
e = mat(imat,1);
nu = mat(imat,2);
axi = mat(imat,3);
fx = mat(imat,4);
fy = mat(imat,5);
% element topology
etop = top(ielem,:);
etop = top(ielem, :);
% position of the local degrees of freedom in the global solution array
[vpos,vshp] = ele_i(1,[],etop,mat,5);
[n,dn,intworg,intcoord] = ele_s([],[],etop,mat,vshp(1,:));
[vpos, vshp] = ele_i( 1, [], etop, mat.mat, 5);
[n, dn, intworg] = ele_s([], [], etop, mat.mat, vshp(1, :));
% number of dofs for the displacements
nedofu = length(nonzeros(vpos(1:2,:)));
nedofu = length(nonzeros(vpos(1:2, :)));
% number of integration points
nint = length(intworg);
......@@ -59,86 +53,79 @@ iu(1:2:nedofu) = nonzeros(vpos(1,:));
iu(2:2:nedofu) = nonzeros(vpos(2,:));
% get nodal displacements
u=sol(pos(ielem,iu));
u = sol(pos(ielem, iu));
% compute shear modulus and compression modulus from E and Nu
g =e/(2*(1+nu));
kappa=e/(3*(1-2*nu));
if axi<2, % plane strain or axi-symmetric
% D matrix currently assumes plane strain
d= 1/3*g*[
4 -2 0 -2;
-2 4 0 -2;
0 0 3 0;
-2 -2 0 4];
d=d+kappa*[
1 1 0 1;
1 1 0 1;
0 0 0 0;
1 1 0 1];
elseif axi==2, % plane stress
f0=(-2*g/3+kappa)^2/(4*g/3+kappa);
f1=(4*g/3+kappa)-f0;
f2=(-2*g/3+kappa)-f0;
d=[f1 f2 0 0;
f2 f1 0 0;
0 0 g 0;
0 0 0 0];
G = E/(2*(1+nu));
kappa = E/(3*(1-2*nu));
if axi < 2 % plane strain or axi-symmetric
% sigma = K*tr(epsilon)*I + 2*G*epsilon^d
% epsilon = [xx yy 2xy zz]
H = G/3*[
4 -2 0 -2;
-2 4 0 -2;
0 0 3 0;
-2 -2 0 4];
H = H + kappa*[
1 1 0 1;
1 1 0 1;
0 0 0 0;
1 1 0 1];
elseif axi == 2 % plane stress
% sigma = K*tr(epsilon)*I + 2G*epsilon^d
% epsilon = [xx yy 2xy zz]
% sigma_zz == 0 (see also exercise 18.1.e in the book)
f0 = (-2*G/3 + kappa)^2/(4*G/3 + kappa);
f1 = ( 4*G/3 + kappa) - f0;
f2 = (-2*G/3 + kappa) - f0;
H = [f1 f2 0 0;
f2 f1 0 0;
0 0 G 0;
0 0 0 0];
end
% initialize strain-displacement matrix
b=zeros(4,nedofu);
ii=0:2:nedofu-2;
sig=zeros(nint,4);
str=zeros(nint,4);
B = zeros(4, nedofu);
ii = 0:2:nedofu-2;
% initialise sigma
sig = zeros(nint, 4);
str = zeros(nint, 4);
% perform the element integration
for int=1:nint
% compute x (=r)-coordinate for axi-symmetric case
x=n(int,:)*nodcoord(:,1);
if (axi==1),
intw=2*pi*x*intworg;
else
intw=intworg;
end
% get shape function derivatives with respect to global coordinates
[dndxy,detj]=getderiv(dn,nodcoord,int);
intwdetj=intw(int)*detj;
% compute the strain displacement matrix
b(1,ii+1)=dndxy(:,1)';
b(2,ii+2)=dndxy(:,2)';
b(3,ii+1)=dndxy(:,2)';
b(3,ii+2)=dndxy(:,1)';
if axi==1,
b(4,ii+1)=n(int,:)/x;
end
strain = b*u;
str(int,:) = strain';
sig(int,:) = (d*strain)';
for int = 1:nint
% get shape function derivatives with respect to global coordinates
dndxy = getderiv(dn, nodcoord, int);
% compute the strain displacement matrix
B(1, ii+1) = dndxy(:, 1)';
B(2, ii+2) = dndxy(:, 2)';
B(3, ii+1) = dndxy(:, 2)';
B(3, ii+2) = dndxy(:, 1)';
if axi == 1 % axi-symmetric correction
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
%sigma(ielem,:)=sig(:)';
% reshape strain and sigma
sigma = int2elm_e(ielem, [sig str], coord, top, mat);
sigma=int2elm_e(ielem,[sig str],coord,top,mattypes);
%strain=int2elm_e(ielem,str,coord,top,mattypes);
%sigma=[sigma strain];
elweight = ones(nlnodes, 1);
elweight = ones(nlnodes,1);
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
function [par1,par2]=eler_i(ielem,coord,top,mat,ichois);
function [par1, par2] = eler_i(ielem, ~, top, mat, ichois)
% [par1, par2] = eler_i(ielem, ~, top, mat, ichois)
%
% description:
% Element information for eler, output is determined by ichois
%
% This class of elements solves either the incompressible Stokes problem
% or a (in) compressible elasticity problem
% input:
% ielem : element number
% ~ : not used (coord, nodal coordinates)
% top : mesh topology array
%
% parameters in the mat array
% mat( 1) = g : viscosity or Young's modulus
% mat( 2) = kappa : compression modulus if 0, incompressibility is assumed
% mat( 3) = axi : 0 for plane, and 1 for axi-symmetric problems
% mat(11) = ietype : element type, see below
%
% 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)
%
% 1 - linear triangular element
% 2 - quadratic triangular element
% 3 - bi-linear quadrilateral
% 4 - bi-quadratic quadrilateral element
% See also eler, eler_d, eler_s
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
% get some element info
imat = top(ielem, end-1);
ietype = mat(imat, 11);
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
imat = top(ielem,iimat);
ietype = mat(imat,11);
nlnodes = length(nonzeros(top(ielem,1:maxnlnodes)));
% initialise
par1 = []; par2 = [];
par1=[]; par2=[];
if ichois == 1
% definition of number of degrees of freedom for each node
par1 = 2*ones(1, nlnodes);
elseif ichois == 2
par1 = 8*ones(1, nlnodes);
elseif ichois == 3
% definition of sequence of nodes for plotting purposes
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
% shape function that defines the geometry
if ietype == 1
% linear triangle
gshp = [23 1 2];
elseif ietype == 2
% quadratic triangle
gshp(1, :) = [26 2 7];
elseif ietype == 3
% linear quad
gshp(1, :) = [14 1 3];
elseif ietype == 4
% quadratic quad
gshp(1, :) = [19 1 3];
if (ichois==1),
if (ietype==1),
par1=[2 2 2]; % Linear triangular-element
elseif (ietype==2),
par1=[2 2 2 2 2 2]; % Quadratic element
elseif (ietype==3),
par1=[2 2 2 2 ]; % bilinear element
elseif (ietype==4),
par1=[2 2 2 2 2 2 2 2 2]; % Bi-quadratic 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=8*ones(1,nlnodes);
elseif ietype==2,
par1=8*ones(1,nlnodes);
elseif ietype==3,
par1=8*ones(1,nlnodes);
elseif ietype==4,
par1=8*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),
end
par1 = gshp;
elseif ichois == 5
% definition of location of dofs (vpos) and their shape functions
if ietype == 1
% linear triangle
gshp = [23 1 2];
elseif (ietype==2 ),
% 6-noded triangle
gshp(1,:) = [26 2 7];
elseif (ietype==3),
% 4-noded quad
gshp(1,:) = [14 1 3];
elseif (ietype==4),
% 9-noded quad