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

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

parent f90df6fa
function [Me, Ce, Ke, rhse, dum3] = ...
elcd(ielem, coord, top, mat, pos, sol, soln, dum1, dum2)
function [Me, Ce, Ke, rhse, dum3] = elcd(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse] = 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)
%
% Me : element mass matrix
% Ce : element damping matrix
% Ke : 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)
Me = []; Ce = []; dum3 = [];
% Output:
% Me element mass matrix
% Ce element damping matrix
% Ke element stiffness matrix
% rhse element right hand side
%
% 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)));
dum3 = [];
% 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 Ke and rhs
% 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
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;
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: defintion of location of dofs (vpos)
% par2: defintion 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
% defintion 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,
if ietype==1 | ietype==2 | ietype==3 | ietype==4 | ietype==5,
nvar=4; nmode=1;
if ietype==1, nint=4; end
if ietype==2, nint=9; end
if ietype==3, nint=9; end
if ietype==4, nint=7; end
if ietype==5, nint=7; end
end
par1=[nint nvar nmode];
else
error([' Ichois = ' num2str(ichois) ' not supported ']);
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
vshp = zeros(1, 3);
vshp(1, :) = [26 2 7];
elseif ietype == 3
% linear quad
vpos = zeros(1, nlnodes);
vpos(1, :) = 1:nlnodes;
vshp = zeros(1,3);
vshp(1, :) = [14 1 2];
elseif ietype == 4
% quadratic 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
% 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];
end
elseif ichois == 7
nvar=4; nmode=1;
if ietype == 1, nint = 4; end
if ietype == 2, nint = 9; end
if ietype == 3, nint = 9; end
if ietype == 4, nint = 7; end
if ietype == 5, nint = 7; end
par1 = [nint nvar nmode];
else
error([' Ichois = ' num2str(ichois) ' not supported ']);
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
function [n,dn,intw,intcoord,detj]=elcd_s(ielem,coord,top,mat,shp,intcoord,idet);
function [n, dn, intw, intcoord, detj] = elcd_s(ielem, coord, top, ...
~, shp, intcoord, idet)
% [n, dn, intw, intcoord, detj] = elcd_s(ielem, coord, top, ~, shp, ...
% intcoord, idet)
%
% input:
% input:
% ielem : element number
% coord : nodal coordinates
%
if (nargin<6), intcoord=[]; end
if (nargin<7), idet=[]; end
if isempty(ielem), ielem=1; end
% top : mesh topology array
% ~ : mat structure (not used)
% shp : array with shape function options
% shp(1) inttype: integration type parameter for shape
% shp(2) integration rule (generates ipts for missing intcoord)
% shp(3) nint: number of integration points
% shp(4) (optional) overwrites nsd = 2
% intcoord : (optional) local coordinates of integration points
% idet : (optional) Jacobian
%
% output:
% n : rows with shape function values
% dn : columns with shape function derivative values
% intw : weighing factors for the integration points
% intcoord : (local) coordinates for the integration points
% detj : value of Jacobian per integration point
%
% See also elcd, shape, gauss, gausstri, gaussl, getderiv
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
imat = top(ielem,iimat);
% parse arguments
if nargin < 6, intcoord = []; end
if nargin < 7, idet = []; end
if isempty(ielem), ielem = 1; end
% parse shp
shfunction = shp(1);
intrule = shp(2);
nint = shp(3);
nsd = 2;
if length(shp)==4,
nsd=shp(4);
if length(shp) == 4
% overwrite nsd
nsd = shp(4);
else
nsd = 2;
end
% initialise detj & intw
detj = [];
intw = [];
detj=[]; intw=[];
if (intrule==1 & isempty(intcoord)),
[intw,intcoord]=gauss(nint,nsd);
elseif (intrule==2 & isempty(intcoord)),
[intw,intcoord]=gausstri(nint); intw=intw/2;
elseif (intrule==3 & isempty(intcoord));
[intw,intcoord]=gaussl(nint,nsd);
if isempty(intcoord)
% get integration points
if intrule == 1
[intw, intcoord] = gauss(nint, nsd);
elseif intrule == 2
[intw, intcoord] = gausstri(nint);
intw = intw/2;
elseif intrule == 3
[intw, intcoord] = gaussl(nint, nsd);
end
end
[n, dn] = shape(intcoord, shfunction);
[n,dn]=shape(intcoord,shfunction);