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

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

parent 6209ecc9
......@@ -3,7 +3,7 @@ function [sigma, sig, elweight] = ele_d(ielem, coord, top, mat, pos, ...
% [sigma, sig, elweight] = ele_d(ielem, coord, top, mat, pos, ...
% sol, varargin)
%
% Stresses in a linear elastic Hookean element
% Stresses c du/dx for ele
%
% input:
%
......@@ -19,7 +19,7 @@ function [sigma, sig, elweight] = ele_d(ielem, coord, top, mat, pos, ...
% sigma(ielem, :) = [sigma_xx sigma_yy sigma_xy sigma_tt .... sigma_xx sigma_yy sigma_xy sigma_tt];
% node_1 node_n
%
% See also ele
% See also ele, femlin_e
% get the number of nodes on this element
......
function [par1, par2] = ele_i(ielem, coord, top, mat, ichois)
function [par1, par2] = ele_i(ielem, ~, top, mat, ichois)
% [par1, par2] = ele_i(ielem, ~, top, mat, ichois)
%
% description:
% Element information for ele, output is determined by ichois
%
% 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
% input:
% ielem : element number
% ~ : not used (coord, nodal coordinates)
% top : mesh topology array
%
% 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)
%
% 1 - linear triangular element
% 2 - quadratic triangular element
% 3 - bi-linear quadrilateral
% 4 - bi-quadratic quadrilateral element
% See also ele, ele_d, ele_s
imat = top(ielem, end-1);
ietype = mat(imat, 11);
......
function [n, dn, intw, intcoord, detj] = ele_s(ielem, coord, top, ~, ...
shp, intcoord, idet)
% [n, dn, intw, intcoord, detj] = ele_s(ielem, coord, top, ~, shp, ...
% intcoord, idet)
%
% input:
% ielem : element number
......@@ -18,26 +20,31 @@ function [n, dn, intw, intcoord, detj] = ele_s(ielem, coord, top, ~, ...
% n : rows with shape function values
% dn : columns with shape function derivative values
% intw : weighing factors for the integration points
% intcoord : (local) coordiantes for the integration points
% detj : Jacobian
% intcoord : (local) coordinates for the integration points
% detj : value of Jacobian per integration point
%
% See also ele, shape, gauss, gausstri, gaussl, getderiv
% 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);
% overwrite nsd
nsd = shp(4);
else
nsd = 2;
end
% initialise detj & intw
detj = []; intw = [];
%
if (intrule == 1 && isempty(intcoord))
[intw, intcoord] = gauss(nint, nsd);
elseif (intrule == 2 && isempty(intcoord))
......
......@@ -18,10 +18,7 @@ function [Me, Ce, Ke, rhse] = elm1dcd(ielem, coord, top, mat)
% rhse : element right hand side
% 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), :);
......@@ -36,8 +33,8 @@ fsupg = mat(imat, 4); % if fsupg ~= 0, use supg, otherwise galerkin
etop = top(ielem, :);
% position of the local degrees of freedom in the global solution array
[vpos, vshp] = elm1dcd_i( 1, [], etop, mat, 5);
[n, dndxi, intw, intcoord] = elm1dcd_s([], [], etop, mat, vshp(1, :));
[vpos, vshp] = elm1dcd_i( 1, [], etop, mat, 5);
[n, dndxi, intw] = elm1dcd_s([], [], etop, mat, vshp(1, :));
nint = length(intw);
ndofu = length(nonzeros(vpos));
......@@ -54,6 +51,7 @@ if fsupg > 0
alpha = coth(beta)-1/beta;
end
% loop over integration points
for int = 1:nint
% get shape function derivatives with respect to global coordinates
......@@ -87,3 +85,4 @@ for int = 1:nint
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
function [sigma]=elm1dcd_d(ielem,coord,top,mat,pos,sol);
function cdudx = elm1dcd_d(ielem, coord, top, mat, pos, sol)
% cdudx = elm1dcd_d(ielem, coord, top, mat, pos, sol, varargin)
%
% linear elastic element
% Flux c du/dx for elm1dcd
%
% Parameters:
% input:
%
% ielem : element number
% Ke : element stiffness matrix
% rhse : element right hand side
% coord : coordinates of all nodes
% top : topology array
% mat : material properties
% sol : solution array
% pos : pos array, see mpos
% f = mat(imat, 1); (source term)
% c = mat(imat, 2); (diffusion constant)
%
% output:
% cdudx(ielem, :, itime) = row with flux c du/dx for each node of the
% element
%
% See also elm1dcd, elm1dcd_i, elm1dcd_s, fem1dcd
% 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)));
% 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 parameters
imat=top(ielem,nlnodes+1);
c = mat(imat,1);
f = mat(imat,3);
imat = top(ielem, end-1);
c = mat(imat, 1);
% element topology
etop=top(ielem,:);
etop = top(ielem,:);
% position of the local degrees of freedom in the global solution array
[vpos,vshp] = elm1dcd_i(1,[],etop,mat,5);
[n,dndxi,intw,intcoord] = elm1dcd_s([],[],etop,mat,vshp(1,:));
nint=length(intw);
ndofu=length(nonzeros(vpos));
[~, vshp] = elm1dcd_i( 1,[], etop, mat, 5);
[n, dndxi, intw] = elm1dcd_s([],[], etop, mat, vshp(1, :));
% initialize Ke and rhs
Ke=zeros(ndofu,ndofu);
rhse=zeros(ndofu,1);
% nodal solutions
unod = sol(pos(ielem, :));
% get the nodal solution
unod=sol(pos(ielem,:));
% stresses at the integration points
sigmaint=zeros(2,1);
for int = 1:nint
% allocate for flux at the integration points
cdudxint = zeros(2, 1);
for int = 1:length(intw)
% get shape function derivatives with respect to global coordinates
[dndxy, detj] = getderiv(dndxi, nodcoord, int);
dndxy = getderiv(dndxi, nodcoord, int);
% compute the strain displacement matrix
b = dndxy(:, 1)';
% stress (c du/dx) at this integration point
sigmaint(int) = c*b*unod;
cdudxint(int) = c*b*unod;
end
% sigmaint can be extrapolated to the nodes
sigma = (n\sigmaint)';
% cdudxint can be extrapolated to the nodes
cdudx = (n\cdudxint)';
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
function [par1,par2]=elm1dcd_i(ielem,coord,top,mat,ichois);
function [par1, par2] = elm1dcd_i(ielem, ~, top, mat, ichois)
% [par1, par2] = elm1dcd_i(ielem, ~, top, mat, ichois)
%
% description:
% Element information for elm1dcd, 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
%
% 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)
%
% 1 - 2-noded linear element
% 2 - 3-noded quadratic element
% See also elm1dcd, elm1dcd_d, elm1dcd_s
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
% get some element info
imat = top(ielem, end-1);
ietype = mat(imat ,5);
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
imat = top(ielem,iimat);
ietype = mat(imat,5);
nlnodes = length(nonzeros(top(ielem,1:maxnlnodes)));
% initialise
par1 = []; par2 = [];
par1=[]; par2=[];
% definition of number of degrees of freedom for each node
if (ichois==1),
if (ietype==1),
par1=[1 1]; % Linear element
elseif (ietype==2),
par1=[1 1 1]; % quadratic element
end
npar=length(par1);
if npar~=nlnodes,
string_err = ['Number of nodes = ' num2str(nlnodes) ...
if ichois == 1
% definition of number of degrees of freedom for each node
if ietype == 1
% linear element
par1 = [1 1];
elseif ietype == 2
% quadratic element
par1 = [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),
par1=3*ones(1,nlnodes);
elseif ietype==2,
par1=4*ones(1,nlnodes);
elseif ietype==4,
par1=4*ones(1,nlnodes);
end
elseif (ichois==3),
% definition of sequence of nodes for plotting purposes
if (ietype==1),
par1=[1 2];
elseif (ietype==2),
par1=[1 2 3];
end
elseif (ichois==4),
% shape function that defines the geometry
if (ietype==1),
% linear
gshp = [2 1 2];
elseif (ietype==2),
% quadratic
gshp(1,:) = [3 1 3];
end
par1=gshp;
elseif (ichois==5),
% defintion of location of dofs (vpos) and their shape functions
if (ietype==1),
% linear element
vpos = zeros(1,2);
vpos(1,:) = 1:2;
vshp = zeros(1,3);
vshp(1,:) = [2 1 2];
elseif (ietype==2),
% quadratic element
vpos = zeros(1,nlnodes);
vpos(1,:) = 1:3;
vshp = zeros(1,3);
vshp(1,:) = [3 1 3]
end
par1=vpos; par2=vshp;
error([string_err num2str(npar)]);
end
elseif ichois == 2
if ietype == 1
% linear element
par1 = 3*ones(1, nlnodes);
elseif ietype== 2
% quadratic element
par1 = 4*ones(1, nlnodes);
elseif ietype == 4
par1 = 4*ones(1, nlnodes);
end
elseif ichois == 3
% definition of sequence of nodes for plotting purposes
if ietype == 1
% linear element
par1 = [1 2];
elseif ietype == 2
% quadratic element
par1 = [1 2 3];
end
elseif ichois == 4
% shape function that defines the geometry
if ietype == 1
% linear element
gshp = [2 1 2];
elseif ietype == 2
% quadratic element
gshp(1, :) = [3 1 3];
end
par1 = gshp;
elseif ichois == 5
% defintion of location of dofs (vpos) and their shape functions
if ietype == 1
% linear element
vpos = zeros(1, 2);
vpos(1, :) = 1:2;
vshp = zeros(1, 3);
vshp(1, :) = [2 1 2];
elseif ietype == 2
% quadratic element
vpos = zeros(1, nlnodes);
vpos(1, :) = 1:3;
vshp = zeros(1, 3);
vshp(1, :) = [3 1 3];
end
par1=vpos; par2=vshp;
else
error([' Ichois = ' num2str(ichois) ' not supported ']);
error([' Ichois = ' num2str(ichois) ' not supported ']);
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
function [n,dn,intw,intcoord,detj]=elm1dcd_s(ielem,coord,top,mat,shp,intcoord,idet);
function [n, dn, intw, intcoord, detj] = elm1dcd_s(ielem, coord, top, ...
~, shp, intcoord, idet)
% [n, dn, intw, intcoord, detj] = elm1dcd_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 elm1dcd, shape, gauss, 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 = 1;
if length(shp)==4,
nsd=shp(4);
if length(shp) == 4
% overwrite nsd
nsd = shp(4);
else
nsd = 1;
end
% initialise detj & intw
detj = [];
intw = [];
detj=[]; intw=[];
if (intrule==1 & isempty(intcoord)),
[intw,intcoord]=gauss(nint,nsd);
if intrule == 1 && isempty(intcoord)
% get integration points
[intw, intcoord] = gauss(nint, nsd);
end
% get shape functions and derivatives
[n, dn] = shape(intcoord, shfunction);
[n,dn]=shape(intcoord,shfunction);
if ~isempty(idet),
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
nlnodes = length(nonzeros(top(ielem,1:maxnlnodes)));
nodcoord=coord(top(ielem,1:nlnodes),:);
nint=length(intcoord(:,1));
detj=zeros(1,nint);
for int=1:nint
[dndxy,detj(int)]=getderiv(dn,nodcoord,int);
end
% get Jacobian
if ~isempty(idet)
% nodal coordinates
nodcoord = coord(top(ielem, 1:end-2), :);
% number of integration points
nint = length(intcoord(:, 1));
% allocate
detj = zeros(1, nint);
for int = 1:nint
[~, detj(int)] = getderiv(dn, nodcoord, int);
end
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
......@@ -57,9 +57,9 @@ if istat == 1
sol = solvestat(q, rhs, bndcon, [], dest);
% calculate flux (c du/dx)
sigma = zeros(size(top, 1), size(top, 2)-2);
cdudx = zeros(size(top, 1), size(top, 2)-2);
for ielem = 1:size(top, 1)
sigma(ielem, :) = elm1dcd_d(ielem, coord, top, mat.mat, ...
cdudx(ielem, :) = elm1dcd_d(ielem, coord, top, mat.mat, ...
pos, sol);
end
......@@ -70,9 +70,9 @@ elseif istat == 2
soln = sol(:, 1);
% calculate flux (c du/dx)
sigma = zeros(size(top, 1), size(top, 2)-2, ntime + 1);
cdudx = zeros(size(top, 1), size(top, 2)-2, ntime + 1);
for ielem = 1:size(top, 1)
sigma(ielem, :, 1) = elm1dcd_d(ielem, coord, top, ...
cdudx(ielem, :, 1) = elm1dcd_d(ielem, coord, top, ...
mat.mat, pos, soln);
end
......@@ -84,7 +84,7 @@ elseif istat == 2
% calculate flux (c du/dx)
for ielem = 1:size(top, 1)
sigma(ielem, :, itime + 1) = elm1dcd_d(ielem, coord, top, ...
cdudx(ielem, :, itime + 1) = elm1dcd_d(ielem, coord, top, ...
mat.mat, pos, soln);
end
end
......
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