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

Education version 2016

parent e0bcecce
Version 16, Education version 2016
==================================
https://gitlab.tue.nl/stem/mlfem/tree/v16
The original scripts of version 01 weer adapted for educational use over the years.
Most of those adaptations were made by Cees Oomens and a team of student assistants. Sadly, history has not been so kind as to preserve the identities of those students for prosperity.
This is the state of the code in the academic year 2016/2017.
Version 01, Original scripts by Frank Baaijens
==============================================
......
......@@ -76,7 +76,7 @@ for int=1:nint
r0=n(int,:)*nodcoord0(:,1);
if (axi==1),
intw=r*intworg;
intw=2*pi*r*intworg;
else
intw=intworg;
end
......
function [qem,qed,qe,rhse,Fhist]=elcneo(ielem,coord,top,mat,pos,sol,dum1,Fhistn,dum2);
%
% neohookean material behaviour
%
function [qem,qed,qe,rhse,Fhist]=elcneo(ielem,coord,top,mat,pos,sol,dum1,Fhistn,dum2);
%
% neohookean material behaviour
%
% mat(1) = G: modulus
% mat(2) = kappa : compression modulus
% mat(3) = 0: plane strain, 1: axi-symmetric
%
% set mass matrix and damping matrix to zero
qem=[]; qed=[]; tau=[];
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top));
nlnodes=length(nonzeros(top(ielem,1:maxnlnodes)));
% get material parameters
imat = top(ielem,iimat);
G = mat(imat,1);
kappa= mat(imat,2);
% mat(2) = kappa : compression modulus
% mat(3) = 0: plane strain, 1: axi-symmetric
%
% set mass matrix and damping matrix to zero
qem=[]; qed=[]; tau=[];
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top));
nlnodes=length(nonzeros(top(ielem,1:maxnlnodes)));
% get material parameters
imat = top(ielem,iimat);
G = mat(imat,1);
kappa= mat(imat,2);
axi = mat(imat,3);
etop = top(ielem,:);
[vpos,vshp] = elcneo_i(1,[],etop,mat,5);
[n,dn,intworg,intcoord] = elcneo_s([],[],etop,mat,vshp(1,:));
nedof = length(nonzeros(vpos(1:2,:)));
nint = length(intworg);
iu = zeros(1,nedof);
iu(1:2:nedof) = nonzeros(vpos(1,:));
iu(2:2:nedof) = nonzeros(vpos(2,:));
% * get incremental displacements and current pressure per element
u=sol(pos(ielem,iu));
% * get element coordinates at reference configuration
nodcoord0=zeros(nlnodes,2);
nodcoord0(:,1)=coord(top(ielem,1:nlnodes),1);
nodcoord0(:,2)=coord(top(ielem,1:nlnodes),2);
% get current nodal coordinates
nodcoord=zeros(nlnodes,2);
ii=0:2:(nlnodes-1)*2;
nodcoord(:,1)=nodcoord0(:,1)+u(ii+1);
nodcoord(:,2)=nodcoord0(:,2)+u(ii+2);
% * set element matrix and vector to zero
ndof= length(n(1,:))*2;
qe = zeros(ndof,ndof);
rhse = zeros(ndof,1);
% * initialize the deformation tensor f
F=zeros(3,3);
% define the unit matrix
%
unit=eye(3);
b=zeros(4,ndof);
etop = top(ielem,:);
[vpos,vshp] = elcneo_i(1,[],etop,mat,5);
[n,dn,intworg,intcoord] = elcneo_s([],[],etop,mat,vshp(1,:));
nedof = length(nonzeros(vpos(1:2,:)));
nint = length(intworg);
iu = zeros(1,nedof);
iu(1:2:nedof) = nonzeros(vpos(1,:));
iu(2:2:nedof) = nonzeros(vpos(2,:));
% * get incremental displacements and current pressure per element
u=sol(pos(ielem,iu),1);
% * get element coordinates at reference configuration
nodcoord0=zeros(nlnodes,2);
nodcoord0(:,1)=coord(top(ielem,1:nlnodes),1);
nodcoord0(:,2)=coord(top(ielem,1:nlnodes),2);
% get current nodal coordinates
nodcoord=zeros(nlnodes,2);
ii=0:2:(nlnodes-1)*2;
nodcoord(:,1)=nodcoord0(:,1)+u(ii+1);
nodcoord(:,2)=nodcoord0(:,2)+u(ii+2);
% * set element matrix and vector to zero
ndof= length(n(1,:))*2;
qe = zeros(ndof,ndof);
rhse = zeros(ndof,1);
% * initialize the deformation tensor f
F=zeros(3,3);
% define the unit matrix
%
unit=eye(3);
b=zeros(4,ndof);
nint=length(intworg);
Fhist=zeros(1,9*nint);
for int=1:nint
r =n(int,:)*nodcoord(:,1);
r0=n(int,:)*nodcoord0(:,1);
if (axi==1),
intw=r*intworg;
else
intw=intworg;
end
% * get shape function derivatives with respect to global coordinates
[dndxy,detj] = getderiv(dn,nodcoord,int);
[dndxy0] = getderiv(dn,nodcoord0,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(4,ii+2)=dndxy(:,1)';
%b(5,ii+1)=axi*n(int,:)/r;
%
% get deformation tensor
F(1,1)=dndxy0(:,1)'*nodcoord(:,1);
F(1,2)=dndxy0(:,2)'*nodcoord(:,1);
F(2,1)=dndxy0(:,1)'*nodcoord(:,2);
F(2,2)=dndxy0(:,2)'*nodcoord(:,2);
F(3,3)=axi*(r/r0);
for int=1:nint
r =n(int,:)*nodcoord(:,1);
r0=n(int,:)*nodcoord0(:,1);
if (axi==1),
intw=2*pi*r*intworg;
else
intw=intworg;
end
% * get shape function derivatives with respect to global coordinates
[dndxy,detj] = getderiv(dn,nodcoord,int);
[dndxy0] = getderiv(dn,nodcoord0,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(4,ii+2)=dndxy(:,1)';
%b(5,ii+1)=axi*n(int,:)/r;
%
% get deformation tensor
F(1,1)=dndxy0(:,1)'*nodcoord(:,1);
F(1,2)=dndxy0(:,2)'*nodcoord(:,1);
F(2,1)=dndxy0(:,1)'*nodcoord(:,2);
F(2,2)=dndxy0(:,2)'*nodcoord(:,2);
F(3,3)=axi*(r/r0);
if axi==0, F(3,3)=1; end
ihist=(int-1)*9+(1:9);
Fn=reshape(Fhistn(ihist),3,3);
F=F*Fn;
Fhist(ihist) = F(:)';
% ihist=(int-1)*9+(1:9)
% Fn=reshape(Fhistn(ihist),3,3)
%
% F=F*Fn
% Fhist(ihist) = F(:)';
% volume ratio
J=det(F);
......@@ -115,7 +115,7 @@ for int=1:nint
B=F*F';
% Cauchy stress tensor
sigma = kappa *(J-1)*eye(3) + (G/J)*(B-J^(2/3)*eye(3));
sigma = kappa *(J-1)*eye(3) + (G/J)*(B-J^(2/3)*eye(3));
M = (kappa*J + (G/3)*J^(-1/3))*eye(3) - (G/J)*B;
......@@ -144,7 +144,7 @@ for int=1:nint
B(1,2) B(1,2) B(2,2) B(1,1)
];
% Element stiffness matrix
qe = qe + b'*(DF + Ds + DJ)*b*intw(int)*detj;
......@@ -153,11 +153,9 @@ for int=1:nint
rhse = rhse - b'*s*intw(int)*detj;
end
end
s=s;
function [qem,qed,qe,rhse,Fhist]=elcneo(ielem,coord,top,mat,pos,sol,dum1,Fhistn,dum2);
%
% neohookean material behaviour
%
% mat(1) = G: modulus
% mat(2) = kappa : compression modulus
% mat(3) = 0: plane strain, 1: axi-symmetric
%
% set mass matrix and damping matrix to zero
qem=[]; qed=[]; tau=[];
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top));
nlnodes=length(nonzeros(top(ielem,1:maxnlnodes)));
% get material parameters
imat = top(ielem,iimat);
G = mat(imat,1);
kappa= mat(imat,2);
axi = mat(imat,3);
etop = top(ielem,:);
[vpos,vshp] = elcneo_i(1,[],etop,mat,5);
[n,dn,intworg,intcoord] = elcneo_s([],[],etop,mat,vshp(1,:));
nedof = length(nonzeros(vpos(1:2,:)));
nint = length(intworg);
iu = zeros(1,nedof);
iu(1:2:nedof) = nonzeros(vpos(1,:));
iu(2:2:nedof) = nonzeros(vpos(2,:));
% * get incremental displacements and current pressure per element
u=sol(pos(ielem,iu),1);
% * get element coordinates at reference configuration
nodcoord0=zeros(nlnodes,2);
nodcoord0(:,1)=coord(top(ielem,1:nlnodes),1);
nodcoord0(:,2)=coord(top(ielem,1:nlnodes),2);
% get current nodal coordinates
nodcoord=zeros(nlnodes,2);
ii=0:2:(nlnodes-1)*2;
nodcoord(:,1)=nodcoord0(:,1)+u(ii+1);
nodcoord(:,2)=nodcoord0(:,2)+u(ii+2);
% * set element matrix and vector to zero
ndof= length(n(1,:))*2;
qe = zeros(ndof,ndof);
rhse = zeros(ndof,1);
% * initialize the deformation tensor f
F=zeros(3,3);
% define the unit matrix
%
unit=eye(3);
b=zeros(4,ndof);
nint=length(intworg);
Fhist=zeros(1,9*nint);
for int=1:nint
r =n(int,:)*nodcoord(:,1);
r0=n(int,:)*nodcoord0(:,1);
if (axi==1),
intw=2*pi*r*intworg;
else
intw=intworg;
end
% * get shape function derivatives with respect to global coordinates
[dndxy,detj] = getderiv(dn,nodcoord,int);
[dndxy0,detj0] = getderiv(dn,nodcoord0,int);
% * compute the strain displacement matrix , current coordinates
b(1,ii+1)=dndxy(:,1)';
b(2,ii+2)=dndxy(:,2)';
b(3,ii+1)=dndxy(:,2)';
b(4,ii+2)=dndxy(:,1)';
%b(5,ii+1)=axi*n(int,:)/r;
% * compute the strain displacement matrix , undeformed coordinates
b0(1,ii+1)=dndxy0(:,1)';
b0(2,ii+2)=dndxy0(:,2)';
b0(3,ii+1)=dndxy0(:,2)';
b0(4,ii+2)=dndxy0(:,1)';
%b(5,ii+1)=axi*n(int,:)/r;
%
% get deformation tensor
F(1,1)=dndxy0(:,1)'*nodcoord(:,1);
F(1,2)=dndxy0(:,2)'*nodcoord(:,1);
F(2,1)=dndxy0(:,1)'*nodcoord(:,2);
F(2,2)=dndxy0(:,2)'*nodcoord(:,2);
F(3,3)=axi*(r/r0);
if axi==0, F(3,3)=1; end
% ihist=(int-1)*9+(1:9);
% Fn=reshape(Fhistn(ihist),3,3);
%
% F=F*Fn;
% Fhist(ihist) = F(:)';
% volume ratio
J=det(F);
% Finger tensor
B=F*F';
% Cauchy stress tensor
sigma = kappa *(J-1)*eye(3) + (G/J)*(B-J^(2/3)*eye(3));
MF=inv(F)*sigma*J;
DF = - [
MF(1,1) 0 MF(1,2) 0
0 MF(2,2) 0 MF(2,1)
MF(2,1) 0 MF(2,2) 0
0 MF(1,2) 0 MF(1,1)];
DJ = [
MF(1,1) MF(1,1) 0 0
MF(2,2) MF(2,2) 0 0
MF(1,2) MF(1,2) 0 0
MF(2,1) MF(2,1) 0 0];
M=inv(F)*J*((kappa*J + (G/3)*J^(-1/3))*eye(3) - (G/J)*B);
M2=inv(F)*2*B;
Ds = [
M(1,1) M(1,1) 0 0
M(2,2) M(2,2) 0 0
M(2,1) M(2,1) 0 0
M(1,2) M(1,2) 0 0] ...
+ ...
G*[
M2(1,1) 0 M2(1,2) 0
0 M2(2,2) 0 M2(2,1)
M2(2,1) 0 M2(2,2) 0
0 M2(1,2) 0 M2(1,1)
];
% Element stiffness matrix
qe = qe + b0'*(DF + Ds + DJ)*b0*intw(int)*detj0;
% right hand side array
%s=[sigma(1,1) sigma(2,2) sigma(1,2) sigma(1,2)]';
t=inv(F)*sigma;
s=[t(1,1) t(2,2) t(1,2) t(1,2)]';
%rhse = rhse - b'*s*intw(int)*detj;
rhse = rhse - b0'*s*intw(int)*J*detj0;
end
function [par1,par2]=elcneo_i(ielem,coord,top,mat,ichois);
%
% description:
%
% This class of elements solves either the incompressible Stokes problem
% or a (in) compressible elasticity problem
%
% paramters 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 - 4-noded quad with constant pressure
% 2 - 9-noded quad with discontinuous pressure
% 3 - 9-noded quad with continuous bi-linear pressure
% 4 - Mini-element extended linear vel, linear continuous pressure
% 5 - Crouzeix-Raviart element, extended quadratic vel, line, disc. pressure
% 6 - Taylor-Hood element, quadratic u, linear continuous pressure
[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=[3 3 3]; % linear triangle
elseif (ietype==2),
par1=[2 2 2 2 2 2 2 2 2]; % bi-quadratic element
elseif (ietype==4),
par1=[2 2 2 2 2 2]; %quadratic triangle
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),
error(['ichois = 2 not yet supported in elneo_i']);
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 7; 7 8; 8 1];
elseif (ietype==4 ),
par1=[1 2; 2 3; 3 4; 4 5; 5 6; 6 1];
end
elseif (ichois==4),
if (ietype==1),
% 3-noded triangle
gshp = [24 1 2];
elseif (ietype==2),
% 9-noded quad
gshp(1,:) = [19 1 3];
elseif (ietype==4 ),
% 6-noded triangle
gshp(1,:) = [26 2 7];
end
par1=gshp;
elseif (ichois==5),
if (ietype==1),
% 3-noded triangle
vpos = zeros(2,4);
vpos(1,:) = 1:2:5;
vpos(2,:) = vpos(1,:)+1;
vshp = zeros(2,3);
vshp(1,:) = [23 2 4];
vshp(2,:) = [23 2 4];
elseif (ietype==2),
% 9-noded quad
vpos = zeros(2,nlnodes);
vpos(1,:) = 1:2:nlnodes*2;
vpos(2,:) = vpos(1,:)+1;
vshp = zeros(2,3);
vshp(1,:) = [19 1 3];
vshp(2,:) = [19 1 3];
elseif (ietype==4),
% 6-noded triangle with continuous linear pressure
vpos = zeros(2,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];
end
par1=vpos; par2=vshp;
elseif ichois==6,
if ietype==1,
par1=[1 0; 0 1; 0 0; 1/3 1/3];
elseif ietype==2,
par1=[-1 -1; 0 -1; 1 -1; 1 0; 1 1; 0 1; -1 1; -1 0; 0 0];
elseif ietype==4,
par1=[1 0; 1/2 1/2; 0 1; 0 1/2; 0 0; 1/2 0];
end
elseif ichois==7,
if ietype==1 | ietype==2 | ietype==4 ,
nvar=9; nmode=1;
if ietype==1, nint=4; end
if ietype==2, nint=9; end
if ietype==4, nint=7; end
end
par1=[nint nvar nmode];
elseif ichois==12, % Definition of the edges
if ietype==1,
par1=[1 2; 2 3; 3 4; 4 1];