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

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

parent 822fd2b7
......@@ -10,7 +10,7 @@
https://gitlab.tue.nl/stem/mlfem_nac/tree/v21
* [update] fix code comments throughout (mostly in demos)
* [update] cleaning & check of comments/help in most linear element (support) functions (issue #1)
* [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)
......@@ -19,6 +19,7 @@ https://gitlab.tue.nl/stem/mlfem_nac/tree/v21
* [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
......
......@@ -50,7 +50,7 @@ if ichois == 1
% sanity check
npar = length(par1);
if npar~=nlnodes,
if npar ~= nlnodes
string_err = ['Number of nodes = ' num2str(nlnodes) ...
' for element type ' num2str(ietype) ...
' in ielem ' num2str(ielem) ', should be '];
......
......@@ -86,7 +86,7 @@ end
% perform the element integration
for int = 1:length(intworg);
for int = 1:length(intworg)
if axi == 1 % axi-symmetric
x = n(int, :)*nodcoord(:, 1); % ipts x-coordinate (r)
......@@ -114,4 +114,4 @@ for int = 1:length(intworg);
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
end
......@@ -48,9 +48,9 @@ nedofu = length(nonzeros(vpos(1:2, :)));
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,:));
iu = zeros(1, nedofu);
iu(1:2:nedofu) = nonzeros(vpos(1, :));
iu(2:2:nedofu) = nonzeros(vpos(2, :));
% get nodal displacements
u = sol(pos(ielem, iu));
......@@ -129,4 +129,4 @@ sigma = int2elm_e(ielem, [sig str], coord, top, mat);
elweight = ones(nlnodes, 1);
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
\ No newline at end of file
end
function [Me, Ce, Ke, rhse, dum3] = ...
elefib(ielem, coord, top, mat, pos, sol, soln, dum1, dum2)
function [Me, Ce, Ke, rhse, dummy] = elefib(ielem, coord, top, mat, varargin)
% [Me, Ce, Ke, rhse, dummy] = elefib(ielem, coord, top, mat)
%
% linear elastic, plane stres element, including fibers
% linear elastic (Hookean) plane stress element, including fibers
%
% 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)
% nf = mat(imat, 12); (Number of fibers)
% Ef = mat(imat, 13 + (ifiber-1)*3); (Young's modulus of the i'th fiber
% phi = mat(imat, 14 + (ifiber-1)*3); (angle of the i'th fiber, in radians)
% vf = mat(imat, 15 + (ifiber-1)*3); (Volume fraction of the i'th fiber)
%
% Me : element mass matrix (not used)
% Ce : element damping matrix (not used)
% Ke : element stiffness matrix
% rhse : element right hand side
% coord : global node coordinates array
% top : global 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)
% eltype = mat(imat,11)( Element type specification, see ele_i)
% nf = mat(imat,12); (Number of fibers)
% Ef = mat(imat,13+(ifiber-1)*3); (Young's modulus of the i'th fiber
% phi = mat(imat,14+(ifiber-1)*3); (phi of the i'th fiber, in radians)
% vf = mat(imat,15+(ifiber-1)*3); (Volume fraction of the i'th fiber)
% Output:
% Me element mass matrix
% Ce element damping matrix
% Ke element stiffness matrix
% dummy = [], for consistency of element calls
%
% Remark
% mat.types : 'elefib';
% See also elefib_s, elefib_i, elefib_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
% pick up some material parameters
imat = top(ielem, nlnodes + 1);
E = mat(imat, 1);
nu = mat(imat, 2);
......@@ -45,43 +43,36 @@ axi = mat(imat, 3);
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] = elefib_i( 1, [], etop, mat, 5);
[n, dn, intworg] = elefib_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
% initialize Ke and rhse
Ke = zeros(nedofu, nedofu);
rhse = zeros(nedofu, 1);
% initialize strain-displacement matrix
B = zeros(4, nedofu);
ii = 0:2:nedofu-2;
nfiber = mat(imat, 12);
vf_tot = 0;
Df = zeros(3, 3);
Hf = zeros(3, 3);
for ifiber = 1:nfiber
Ef = mat(imat, 13 + 3*(ifiber-1));
phi = mat(imat, 14 + 3*(ifiber-1));
vf = mat(imat, 15 + 3*(ifiber-1));
vf_tot = vf_tot + vf;
ca = cos(phi); sa = sin(phi);
Df = Df + vf*Ef*[
Hf = Hf + vf*Ef*[
ca^4 ca^2*sa^2 sa*ca^3
ca^2*sa^2 sa^4 sa^3*ca
ca^3*sa ca*sa^3 ca^2*sa^2];
end
if axi == 2 % plane stress
% D matrix currently assumes plane strain
d = (1-vf_tot)*E/(1-nu^2)*[
% H matrix currently assumes plane stress
H = (1-vf_tot)*E/(1-nu^2)*[
1 nu 0
nu 1 0
0 0 (1-nu)/2];
......@@ -90,18 +81,14 @@ else
error('This element only support plane stress, axi = 2');
end
d = d + Df;
% initialize strain-displacement matrix
b = zeros(3 ,nedofu);
ii = 0:2:nedofu-2;
H = H + Hf;
% perform the element integration
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
......@@ -110,14 +97,19 @@ for int = 1:nint
[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)';
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
Ke = Ke + b'*d*b*intwdetj;
Ke = Ke + B'*H*B*intwdetj;
end
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
function [sigma,sig,elweight] = ...
elefib_d(ielem,coord,top,mat,pos,sol,nodaldata,elemdata_in,elemparam);
function [sigma, sig, elweight] = elefib_d(ielem, coord, top, mat, pos, ...
sol, varargin)
% [sigma, sig, elweight] = elefib_d(ielem, coord, top, mat, pos, ...
% sol, varargin)
%
% linear elastic element
% Stresses c du/dx for elefib
%
% Parameters:
%
% input:
%
% ielem : element number
% coord : global node coordinates array
% top : global 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)
% 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);
% mat : material properties
% E = mat(imat, 1); (Young's modulus)
% nu = mat(imat, 2); (Poisson's ratio)
% axi = mat(imat, 3); (2: plane stress)
% nf = mat(imat, 12); (Number of fibers)
% Ef = mat(imat, 13 + (ifiber-1)*3); (Young's modulus of the i'th fiber
% phi = mat(imat, 14 + (ifiber-1)*3); (angle of the i'th fiber, in radians)
% vf = mat(imat, 15 + (ifiber-1)*3); (Volume fraction of the i'th fiber)
% 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 elefib, 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)));
% get the number of nodes on this element
nlnodes = length(nonzeros(top(ielem, 1:end-2)));
% get the coordinates of the nodes of this element
nodcoord=coord(top(ielem,1:nlnodes),:);
mattypes=mat;
mat=mat.mat;
% 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,:);
nodcoord = coord(top(ielem, 1:nlnodes), :);
% 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);
% element topology
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] = elefib_i( 1, [], etop, mat.mat, 5);
[~, dn, intworg] = elefib_s([], [], etop, mat.mat, vshp(1, :));
% number of dofs for the displacements
nedofu = length(nonzeros(vpos(1:2,:)));
% number of integration points
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,:));
% get nodal displacements
u=sol(pos(ielem,iu));
iu = zeros(1, nedofu);
iu(1:2:nedofu) = nonzeros(vpos(1, :));
iu(2:2:nedofu) = nonzeros(vpos(2, :));
% get nodal displacements
u = sol(pos(ielem, iu));
nfiber = mat(imat,12);
vf_tot=0;
Df=zeros(3,3);
for ifiber=1:nfiber
Ef = mat(imat,12+(ifiber-1)*3+1);
angle=mat(imat,12+(ifiber-1)*3+2);
vf = mat(imat,12+(ifiber-1)*3+3);
vf_tot = 0;
Hf = zeros(3, 3);
for ifiber = 1:nfiber
Ef = mat(imat, 12 + (ifiber-1)*3 + 1);
angle = mat(imat, 12 + (ifiber-1)*3 + 2);
vf = mat(imat, 12 + (ifiber-1)*3 + 3);
vf_tot = vf_tot + vf;
ca=cos(angle); sa=sin(angle);
Df = Df + vf*Ef*[
ca = cos(angle); sa = sin(angle);
Hf = Hf + vf*Ef*[
ca^4 ca^2*sa^2 sa*ca^3
ca^2*sa^2 sa^4 sa^3*ca
ca^3*sa ca*sa^3 ca^2*sa^2];
end
if axi==2, % plane strain or axi-symmetric
% D matrix currently assumes plane strain
d= (1-vf_tot)*E/(1-nu^2)*[
1 nu 0
nu 1 0
0 0 (1-nu)/2];
if axi == 2 % plane stress
% H matrix currently assumes plane strain
H = (1-vf_tot)*E/(1-nu^2)*[
1 nu 0
nu 1 0
0 0 (1-nu)/2];
else
error(['This element only support plane stress, axi=1']);
error('This element only support plane stress, axi = 2');
end
d=d+Df;
H = H + Df;
% initialize strain-displacement matrix
b=zeros(3,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)';
strain = b*u;
str(int,:) = [strain ; 0]';
sig(int,:) = [(d*strain) ; 0]';
end
%sigma(ielem,:)=sig(:)';
sigma=int2elm_e(ielem,[sig str],coord,top,mattypes);
%strain=int2elm_e(ielem,str,coord,top,mattypes);
%sigma=[sigma strain];
elweight = ones(nlnodes,1);
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)';
% calculate strain
strain = B*u;
% store strain and sigma
str(int, :) = strain';
sig(int, :) = (H*strain)';
end
% reshape strain and sigma
sigma = int2elm_e(ielem, [sig str], coord, top, mat);
elweight = ones(nlnodes,1);
% part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
end
function [par1,par2]=elefib_i(ielem,coord,top,mat,ichois);
%
% description:
%
%
%
% 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:
%
% 1 - linear triangular element
% 2 - quadratic triangular element
% 3 - bi-linear quadrilateral
% 4 - bi-quadratic quadrilateral element
[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=[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),
% 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
gshp(1,:) = [19 1 3];
end
par1=gshp;
elseif (ichois==5),
if (ietype==1),
% 3-noded triangle with linear pressure
vpos = zeros(2,3);
vpos(1,:) = 1:2:5;
vpos(2,:) = vpos(1,:)+1;
vshp = zeros(2,3);
vshp(1,:) = [23 2 3];
vshp(2,:) = [23 2 3];
elseif (ietype==2),
% 6-noded triangle
vpos = zeros(3,nlnodes);
vpos(1,:) = 1:2:nlnodes*2;
vpos(2,:) = vpos(1,:)+1;
vshp = zeros(2,3);
vshp(1,:) = [26 2 7];
vshp(2,:) = [26 2 7];
elseif (ietype==3),
% 4-noded quad
vpos = zeros(2,nlnodes);
vpos(1,:) = 1:2:nlnodes*2;
vpos(2,:) = vpos(1,:)+1;
vshp = zeros(2,3);
vshp(1,:) = [14 1 2];
vshp(2,:) = [14 1 2];
elseif (ietype==4),
% 9-noded quad
vpos = zeros(2,nlnodes);
vpos(1,:) = 1:2:nlnodes*2;
vpos(2,:) = vpos(1,:)+1;