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

original code by Frank Baaijens

parent b788f72f
......@@ -3,3 +3,4 @@ Version 01, Original scripts by Frank Baaijens
https://gitlab.tue.nl/stem/mlfem/tree/v01
These are the original scripts that Frank Baaijens wrote for his personal use at Philips.
......@@ -2,7 +2,9 @@ mlfem
=====
https://gitlab.tue.nl/STEM/mlfem
Software belonging to the text book
Software belonging to the text book
Biomechanics: Concepts and Computation, by Cees Oomens, Marcel Brekelmans and Frank Baaijens (Cambridge University Press, 2009).
using & downloading
......
99.01.22 16:24 B C:\mlfiles\elemfnct\elembld.m --> sg1 /users/wfw2/roelb/mlfiles/elemfnct elembld.m
99.01.22 16:24 B C:\mlfiles\elemfnct\eleminfo.m --> sg1 /users/wfw2/roelb/mlfiles/elemfnct eleminfo.m
99.01.22 16:24 B C:\mlfiles\elemfnct\elemshp.m --> sg1 /users/wfw2/roelb/mlfiles/elemfnct elemshp.m
99.01.29 09:43 B C:\mlfiles\elemfnct\elembld.m --> sg26 /users/sg26/baaijens/mlfiles/elemfnct elembld.m
99.01.29 09:43 B C:\mlfiles\elemfnct\eleminfo.m --> sg26 /users/sg26/baaijens/mlfiles/elemfnct eleminfo.m
99.01.29 09:43 B C:\mlfiles\elemfnct\elemshp.m --> sg26 /users/sg26/baaijens/mlfiles/elemfnct elemshp.m
99.01.29 09:43 B C:\mlfiles\elemfnct\WS_FTP.LOG --> sg26 /users/sg26/baaijens/mlfiles/elemfnct WS_FTP.LOG
99.06.24 13:46 B C:\mlfiles\elemfnct\WS_FTP.LOG --> hc96 /home/rens/mlfiles/elemfnct WS_FTP.LOG
99.06.24 13:46 B C:\mlfiles\elemfnct\elembld.m --> hc96 /home/rens/mlfiles/elemfnct elembld.m
99.06.24 13:46 B C:\mlfiles\elemfnct\elemdrv.m --> hc96 /home/rens/mlfiles/elemfnct elemdrv.m
99.06.24 13:46 B C:\mlfiles\elemfnct\eleminfo.m --> hc96 /home/rens/mlfiles/elemfnct eleminfo.m
99.06.24 13:46 B C:\mlfiles\elemfnct\elemshp.m --> hc96 /home/rens/mlfiles/elemfnct elemshp.m
function [elmass,eldamp,elstif,elrhs,elemdataout]= ...
elembld(ielem,coord,top,mat,pos,sol,nodaldata,elemdatain,param);
%
if (nargin<5),
error('Illegal number of input arguments')
end
if (nargin<6),
sol=[];
end
if (nargin<7),
nodaldata=[];
end
if (nargin<8),
elemdatain=[];
end
if (nargin<9),
param=[];
end
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]= sysinfo(size(coord),size(top));
eltype = top(ielem,iitype);
ntypes=length(mat.types(:,1));
if eltype>ntypes,
error('This type not allowed');
end
type=deblank(mat.types(eltype,:));
[elmass,eldamp,elstif,elrhs,elemdataout]= ...
feval(type,ielem,coord,top,mat.mat,pos,sol,nodaldata,elemdatain,param);
\ No newline at end of file
function [elmdat1,elmdat2,elmdat3]= ...
elemdrv(ielem,coord,top,mat,pos,sol,nodaldata,elemdatain,param);
%
if (nargin<5),
error('Illegal number of input arguments')
end
if (nargin<6),
sol=[];
end
if (nargin<7),
nodaldata=[];
end
if (nargin<8),
elemdatain=[];
end
if (nargin<9),
param=[];
end
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]= sysinfo(size(coord),size(top));
eltype = top(ielem,iitype);
ntypes=length(mat.types(:,1));
if eltype>ntypes,
error('This type not allowed');
end
type=deblank(mat.types(eltype,:));
type=[type '_d'];
[elmdat1,elmdat2,elmdat3]= ...
feval(type,ielem,coord,top,mat,pos,sol,nodaldata,elemdatain,param);
\ No newline at end of file
function [elinfo1,elinfo2]= eleminfo(ielem,coord,top,mat,ichois);
% [elinfo1,elinfo2]=eleminfo(ielem,coord,top,mat,ichois);
%
%
% ichois : 1 - ldof
% 2 - ddof
% 3 - elgraph
% 4 - gshape
% 5 - vshape
% 6 - dshape
% 7:10 - reserved
% >11 - for user
%
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]= sysinfo(size(coord),size(top));
eltype = top(ielem,iitype);
ntypes=length(mat.types(:,1));
if eltype>ntypes,
error('This type not allowed');
end
type=deblank(mat.types(eltype,:));
type=[type '_i'];
[elinfo1,elinfo2]=feval(type,ielem,coord,top,mat.mat,ichois);
function [n,dn,weight,intcoord,detj]=elemshp(ielem,coord,top,mat,ichois,intcoord,idet);
%
%
if nargin<6 ,
intcoord=[];
idet=[];
elseif nargin<7
idet=[];
end
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]= sysinfo(size(coord),size(top));
eltype = top(ielem,iitype);
ntypes=length(mat.types(:,1));
if eltype>ntypes,
error('This type not allowed');
end
type=deblank(mat.types(eltype,:));
type=[type '_s'];
[n,dn,weight,intcoord,detj]= feval(type,ielem,coord,top,mat.mat,ichois,intcoord,idet);
99.01.18 09:33 B C:\mlfiles\elmlib\defdi.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib defdi.m
99.01.18 09:33 B C:\mlfiles\elmlib\elneo.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo.m
99.01.18 09:33 B C:\mlfiles\elmlib\elneo.m_ --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo.m_
99.01.18 09:33 B C:\mlfiles\elmlib\elneo_i.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo_i.m
99.01.18 09:33 B C:\mlfiles\elmlib\elneo_s.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo_s.m
99.01.18 09:33 B C:\mlfiles\elmlib\elucm.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elucm.m
99.01.18 09:33 B C:\mlfiles\elmlib\elucm_i.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elucm_i.m
99.01.18 09:33 B C:\mlfiles\elmlib\elucm_s.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elucm_s.m
99.01.18 09:33 B C:\mlfiles\elmlib\elup.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elup.m
99.01.18 09:33 B C:\mlfiles\elmlib\elup_i.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elup_i.m
99.01.18 09:33 B C:\mlfiles\elmlib\elup_s.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elup_s.m
99.01.18 09:33 B C:\mlfiles\elmlib\int2elm.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib int2elm.m
99.01.22 16:24 B C:\mlfiles\elmlib\defdi.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib defdi.m
99.01.22 16:24 B C:\mlfiles\elmlib\elneo.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elneo.m
99.01.22 16:24 B C:\mlfiles\elmlib\elneo.m_ --> sg1 /users/wfw2/roelb/mlfiles/elmlib elneo.m_
99.01.22 16:24 B C:\mlfiles\elmlib\elneo_i.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elneo_i.m
99.01.22 16:24 B C:\mlfiles\elmlib\elneo_s.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elneo_s.m
99.01.22 16:24 B C:\mlfiles\elmlib\elucm.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elucm.m
99.01.22 16:24 B C:\mlfiles\elmlib\elucm_i.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elucm_i.m
99.01.22 16:24 B C:\mlfiles\elmlib\elucm_s.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elucm_s.m
99.01.22 16:24 B C:\mlfiles\elmlib\elup.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elup.m
99.01.22 16:24 B C:\mlfiles\elmlib\elup_i.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elup_i.m
99.01.22 16:24 B C:\mlfiles\elmlib\elup_s.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib elup_s.m
99.01.22 16:24 B C:\mlfiles\elmlib\int2elm.m --> sg1 /users/wfw2/roelb/mlfiles/elmlib int2elm.m
99.01.22 16:24 B C:\mlfiles\elmlib\WS_FTP.LOG --> sg1 /users/wfw2/roelb/mlfiles/elmlib WS_FTP.LOG
99.01.29 09:43 B C:\mlfiles\elmlib\defdi.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib defdi.m
99.01.29 09:43 B C:\mlfiles\elmlib\elneo.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo.m
99.01.29 09:43 B C:\mlfiles\elmlib\elneo.m_ --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo.m_
99.01.29 09:43 B C:\mlfiles\elmlib\elneo_i.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo_i.m
99.01.29 09:43 B C:\mlfiles\elmlib\elneo_s.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elneo_s.m
99.01.29 09:43 B C:\mlfiles\elmlib\elucm.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elucm.m
99.01.29 09:43 B C:\mlfiles\elmlib\elucm_i.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elucm_i.m
99.01.29 09:43 B C:\mlfiles\elmlib\elucm_s.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elucm_s.m
99.01.29 09:43 B C:\mlfiles\elmlib\elup.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elup.m
99.01.29 09:43 B C:\mlfiles\elmlib\elup_i.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elup_i.m
99.01.29 09:43 B C:\mlfiles\elmlib\elup_s.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib elup_s.m
99.01.29 09:43 B C:\mlfiles\elmlib\int2elm.m --> sg26 /users/sg26/baaijens/mlfiles/elmlib int2elm.m
99.01.29 09:43 B C:\mlfiles\elmlib\WS_FTP.LOG --> sg26 /users/sg26/baaijens/mlfiles/elmlib WS_FTP.LOG
99.04.29 08:42 B C:\mlfiles\elmlib\elneo.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo.m
99.04.29 08:49 B C:\mlfiles\elmlib\elneo.m <-- sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo.m
99.04.29 13:45 B C:\mlfiles\elmlib\elneo.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo.m
99.05.11 16:12 B C:\mlfiles\elmlib\WS_FTP.LOG --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib WS_FTP.LOG
99.05.11 16:12 B C:\mlfiles\elmlib\defdi.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib defdi.m
99.05.11 16:12 B C:\mlfiles\elmlib\elcd.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elcd.m
99.05.11 16:12 B C:\mlfiles\elmlib\elcd_i.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elcd_i.m
99.05.11 16:12 B C:\mlfiles\elmlib\elcd_s.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elcd_s.m
99.05.11 16:12 B C:\mlfiles\elmlib\ele.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib ele.m
99.05.11 16:12 B C:\mlfiles\elmlib\ele_d.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib ele_d.m
99.05.11 16:12 B C:\mlfiles\elmlib\ele_i.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib ele_i.m
99.05.11 16:12 B C:\mlfiles\elmlib\ele_s.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib ele_s.m
99.05.11 16:12 B C:\mlfiles\elmlib\elneo.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo.m
99.05.11 16:12 B C:\mlfiles\elmlib\elneo_i.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo_i.m
99.05.11 16:12 B C:\mlfiles\elmlib\elneo_org.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo_org.m
99.05.11 16:12 B C:\mlfiles\elmlib\elneo_s.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo_s.m
99.05.11 16:12 B C:\mlfiles\elmlib\elucm.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elucm.m
99.05.11 16:12 B C:\mlfiles\elmlib\elucm_i.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elucm_i.m
99.05.11 16:13 B C:\mlfiles\elmlib\elucm_s.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elucm_s.m
99.05.11 16:13 B C:\mlfiles\elmlib\elup.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elup.m
99.05.11 16:13 B C:\mlfiles\elmlib\elup_i.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elup_i.m
99.05.11 16:13 B C:\mlfiles\elmlib\elup_s.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elup_s.m
99.05.11 16:13 B C:\mlfiles\elmlib\int2elm.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib int2elm.m
99.05.11 16:13 B C:\mlfiles\elmlib\int2elm_e.m --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib int2elm_e.m
99.05.11 16:13 B C:\mlfiles\elmlib\elneo.m_ --> sg26.wfw.wtb.tue.nl /users/sg26/baaijens/mlfiles/elmlib elneo.m_
99.06.24 13:46 B C:\mlfiles\elmlib\WS_FTP.LOG --> hc96 /home/rens/mlfiles/elmlib WS_FTP.LOG
99.06.24 13:46 B C:\mlfiles\elmlib\defdi.m --> hc96 /home/rens/mlfiles/elmlib defdi.m
99.06.24 13:46 B C:\mlfiles\elmlib\elcd.m --> hc96 /home/rens/mlfiles/elmlib elcd.m
99.06.24 13:46 B C:\mlfiles\elmlib\elcd_i.m --> hc96 /home/rens/mlfiles/elmlib elcd_i.m
99.06.24 13:46 B C:\mlfiles\elmlib\elcd_s.m --> hc96 /home/rens/mlfiles/elmlib elcd_s.m
99.06.24 13:46 B C:\mlfiles\elmlib\ele.m --> hc96 /home/rens/mlfiles/elmlib ele.m
99.06.24 13:46 B C:\mlfiles\elmlib\ele_d.m --> hc96 /home/rens/mlfiles/elmlib ele_d.m
99.06.24 13:46 B C:\mlfiles\elmlib\ele_i.m --> hc96 /home/rens/mlfiles/elmlib ele_i.m
99.06.24 13:46 B C:\mlfiles\elmlib\ele_s.m --> hc96 /home/rens/mlfiles/elmlib ele_s.m
99.06.24 13:46 B C:\mlfiles\elmlib\elneo.m --> hc96 /home/rens/mlfiles/elmlib elneo.m
99.06.24 13:46 B C:\mlfiles\elmlib\elneo_i.m --> hc96 /home/rens/mlfiles/elmlib elneo_i.m
99.06.24 13:46 B C:\mlfiles\elmlib\elneo_org.m --> hc96 /home/rens/mlfiles/elmlib elneo_org.m
99.06.24 13:46 B C:\mlfiles\elmlib\elneo_s.m --> hc96 /home/rens/mlfiles/elmlib elneo_s.m
99.06.24 13:46 B C:\mlfiles\elmlib\elucm.m --> hc96 /home/rens/mlfiles/elmlib elucm.m
99.06.24 13:46 B C:\mlfiles\elmlib\elucm_i.m --> hc96 /home/rens/mlfiles/elmlib elucm_i.m
99.06.24 13:46 B C:\mlfiles\elmlib\elucm_s.m --> hc96 /home/rens/mlfiles/elmlib elucm_s.m
99.06.24 13:46 B C:\mlfiles\elmlib\elup.m --> hc96 /home/rens/mlfiles/elmlib elup.m
99.06.24 13:46 B C:\mlfiles\elmlib\elup_i.m --> hc96 /home/rens/mlfiles/elmlib elup_i.m
99.06.24 13:46 B C:\mlfiles\elmlib\elup_s.m --> hc96 /home/rens/mlfiles/elmlib elup_s.m
99.06.24 13:46 B C:\mlfiles\elmlib\int2elm.m --> hc96 /home/rens/mlfiles/elmlib int2elm.m
99.06.24 13:46 B C:\mlfiles\elmlib\int2elm_e.m --> hc96 /home/rens/mlfiles/elmlib int2elm_e.m
99.06.24 13:46 B C:\mlfiles\elmlib\elneo.m_ --> hc96 /home/rens/mlfiles/elmlib elneo.m_
function [dp,d1,d2,d3,d4,d5]=defdi(press,s,b,axi)
[m,n]=size(s);
if m==3, s=[s(1,1) s(2,2) s(1,2) s(3,3)]; end
[m,n]=size(b);
if m==3, b=[b(1,1) b(2,2) b(1,2) b(3,3)]; end
d1=-[s(1) 0 s(3)/2 s(3)/2 0;
0 s(2) s(3)/2 -s(3)/2 0;
s(3)/2 s(3)/2 (s(1)+s(2))/4 (-s(1)+s(2))/4 0;
-s(3)/2 s(3)/2 (s(1)-s(2))/4 -(s(1)+s(2))/4 0;
0 0 0 0 s(4);
];
d2=[s(1) 0 s(3)/2 s(3)/2 0;
0 s(2) s(3)/2 -s(3)/2 0;
s(3)/2 s(3)/2 (s(1)+s(2))/4 (-s(1)+s(2))/4 0;
-s(3)/2 s(3)/2 (s(1)-s(2))/4 (s(1)+s(2))/4 0;
0 0 0 0 s(4)];
d3=[s(1) 0 s(3)/2 s(3)/2 0;
0 s(2) s(3)/2 -s(3)/2 0;
s(3)/2 s(3)/2 (s(1)+s(2))/4 (-s(1)+s(2))/4 0;
s(3)/2 -s(3)/2 (-s(1)+s(2))/4 (s(1)+s(2))/4 0;
0 0 0 0 s(4)
];
d4=[b(1) 0 b(3)/2 -b(3)/2 0;
0 b(2) b(3)/2 b(3)/2 0;
b(3)/2 b(3)/2 (b(1)+b(2))/4 (b(1)-b(2))/4 0;
-b(3)/2 b(3)/2 (b(1)-b(2))/4 (b(1)+b(2))/4 0;
0 0 0 0 b(4)
];
d5=[b(1) 0 b(3)/2 -b(3)/2 0;
0 b(2) b(3)/2 b(3)/2 0;
b(3)/2 b(3)/2 (b(1)+b(2))/4 (b(1)-b(2))/4 0;
b(3)/2 -b(3)/2 (-b(1)+b(2))/4 -(b(1)+b(2))/4 0;
0 0 0 0 b(4)] ;
dp= press*[1 0 0 0 0;
0 1 0 0 0;
0 0 .5 0 0;
0 0 0 -.5 0;
0 0 0 0 axi];
function [qme,qde,qke,rhse,dum3]=elcd(ielem,coord,top,mat,pos,sol,soln,dum1,dum2);
%
% convection diffusion element element
%
% Parameters:
%
% qme : element mass matrix
% qde : element damping matrix
% qe : 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)
% fx = mat(imat,2); (factor for convective velocity in x-direction)
% fy = mat(imat,2); (factor for convective velocity in y-direction)
% axi = mat(imat,4); (0: plane strain, 1: axi-symmetric, 2: plane stress)
qme=[]; qde=[]; dum3=[];
% get the number of nodes on this elemen
[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 elemen
nodcoord=coord(top(ielem,1:nlnodes),:);
% pick up some material parameter
imat=top(ielem,nlnodes+1);
c = mat(imat,1);
fx = mat(imat,2);
fy = 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,:));
% number of dofs for the temperature
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 qe and rhs
qke=zeros(nedof,nedof);
qde=qke;
qme=qke;
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==1),
intw=2*pi*x*intworg;
else
intw=intworg;
end
% get the convective velocity
[ax,ay]=elcd_a(x,y,fx,fy);
%* get shape function derivatives with respect to global coordinates
[dndxy,detj]=getderiv(dn,nodcoord,int);
intwdetj=intw(int)*detj;
% diffusion part
qke = qke + c*(dndxy(:,1)*dndxy(:,1)' + dndxy(:,2)*dndxy(:,2)')*intwdetj;
% convection part
qde = qde + (ax*n(int,:)'*dndxy(:,1)' + ay*n(int,:)'*dndxy(:,2)')*intwdetj;
% mass matrix
qme = qme + n(int,:)'*n(int,:)*intwdetj;
end
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];
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
function [n,dn,intw,intcoord,detj]=elcd_s(ielem,coord,top,mat,shp,intcoord,idet);
%
% input:
% ielem : element number
% coord : nodal coordinates
%
if (nargin<6), intcoord=[]; end
if (nargin<7), idet=[]; end
if isempty(ielem), ielem=1; end
[nnodes,nsd,nelem,maxnlnodes,iimat,iitype]=sysinfo(size(coord),size(top));
imat = top(ielem,iimat);
shfunction = shp(1);
intrule = shp(2);
nint = shp(3);
nsd = 2;
if length(shp)==4,
nsd=shp(4);
end
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);
end
[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
end
function [qme,qde,qke,rhse,dum3]=elcd(ielem,coord,top,mat,pos,sol,soln,dum1,dum2);
%