Commit 10ca8ebb by Turnhout, M.C. van

### rename element matrices for consistency with book

parent 48a25687
 function [qme, qde, qke, rhse, dum3] = ... function [Me, Ce, Ke, rhse, dum3] = ... elcd(ielem, coord, top, mat, pos, sol, soln, dum1, dum2) % % convection diffusion element element ... ... @@ -7,9 +7,9 @@ function [qme, qde, qke, rhse, dum3] = ... % % Parameters: % % qme : element mass matrix % qde : element damping matrix % qe : element stiffness matrix % 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 ... ... @@ -20,7 +20,7 @@ function [qme, qde, qke, rhse, dum3] = ... % vy = mat(imat,3); (factor for convective velocity in y-direction) % axi = mat(imat,4); (0: plane strain, 1: axi-symmetric, 2: plane stress) qme = []; qde = []; dum3 = []; Me = []; Ce = []; dum3 = []; % get the number of nodes on this element [nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ... ... ... @@ -53,10 +53,10 @@ 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; % initialize Ke and rhs Ke = zeros(nedof, nedof); Ce = Ke; Me = Ke; rhse = zeros(nedof, 1); % perform the element integration ... ... @@ -80,15 +80,15 @@ for int = 1:nint intwdetj = intw(int)*detj; % diffusion part qke = qke + c*(dndxy(:, 1)*dndxy(:, 1)' + ... Ke = Ke + c*(dndxy(:, 1)*dndxy(:, 1)' + ... dndxy(:, 2)*dndxy(:, 2)')*intwdetj; % convection part qde = qde + (ax*n(int, :)'*dndxy(:, 1)' + ... Ce = Ce + (ax*n(int, :)'*dndxy(:, 1)' + ... ay*n(int, :)'*dndxy(:, 2)')*intwdetj; % mass matrix qme = qme + n(int, :)'*n(int, :)*intwdetj; Me = Me + n(int, :)'*n(int, :)*intwdetj; end ... ...
 function [qme,qde,qke,rhse,dum3]=elcdu(ielem,coord,top,mat,pos,sol,solu,posu,matu); function [Me,Ce,Ke,rhse,dum3]=elcdu(ielem,coord,top,mat,pos,sol,solu,posu,matu); % % linear elastic element % % Parameters: % % qme : element mass matrix (not used) % qde : element damping matrix (not used) % qe : element stiffness matrix % Me : element mass matrix (not used) % Ce : element damping matrix (not used) % Ke : element stiffness matrix % rhse : element right hand side % coord : coordinates of all nodes % top : topology array ... ... @@ -18,7 +18,7 @@ function [qme,qde,qke,rhse,dum3]=elcdu(ielem,coord,top,mat,pos,sol,solu,posu,mat % fx = mat(imat,4); (Volume force in x-direction) % fy = mat(imat,5); (Volume force in y-direction) qme=[]; qde=[]; dum3=[]; Me=[]; Ce=[]; dum3=[]; % get the number of nodes on this element [nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top)); ... ... @@ -65,10 +65,10 @@ nedof = length(nonzeros(vpos(1,:))); % 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; % initialize Ke and rhs Ke=zeros(nedof,nedof); Ce=Ke; Me=Ke; rhse=zeros(nedof,1); ... ... @@ -100,13 +100,13 @@ for int=1:nint intwdetj=intw(int)*detj; % diffusion part qke = qke + c*(dndxy(:,1)*dndxy(:,1)' + dndxy(:,2)*dndxy(:,2)')*intwdetj; Ke = Ke + 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; Ce = Ce + (ax*n(int,:)'*dndxy(:,1)' + ay*n(int,:)'*dndxy(:,2)')*intwdetj; % mass matrix qme = qme + n(int,:)'*n(int,:)*intwdetj; Me = Me + n(int,:)'*n(int,:)*intwdetj; end ... ...
 function [qem,qed,qe,rhse,eldat_out]=elchypo(ielem,coord,top,mat,pos,sol,dum1,eldat_in,dum2); function [Me,Ce,Ke,rhse,eldat_out]=elchypo(ielem,coord,top,mat,pos,sol,dum1,eldat_in,dum2); % % neohookean material behaviour ... ... @@ -9,7 +9,7 @@ function [qem,qed,qe,rhse,eldat_out]=elchypo(ielem,coord,top,mat,pos,sol,dum1,el % % set mass matrix and damping matrix to zero qem=[]; qed=[]; tau=[]; Me=[]; Ce=[]; tau=[]; [nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top)); ... ... @@ -53,7 +53,7 @@ nodcoord(:,2)=nodcoord0(:,2)+u(ii+2); % set element matrix and vector to zero ndof= length(n(1,:))*2; qe = zeros(ndof,ndof); Ke = zeros(ndof,ndof); rhse = zeros(ndof,1); % initialize the deformation tensor f ... ... @@ -153,7 +153,7 @@ for int=1:nint ]; % Element stiffness matrix qe = qe + b'*(DF + Ds + DJ)*b*intw(int)*detj; Ke = Ke + b'*(DF + Ds + DJ)*b*intw(int)*detj; % right hand side array s=[sigma(1,1) sigma(2,2) sigma(1,2) sigma(1,2)]'; ... ...
 function [qem, qed, qe, rhse, Fhist] = ... function [Me, Ce, Ke, rhse, Fhist] = ... elcneo(ielem, coord, top, mat, pos, sol, dum1, Fhistn, dum2) % % neohookean material behaviour ... ... @@ -9,7 +9,7 @@ function [qem, qed, qe, rhse, Fhist] = ... % % set mass matrix and damping matrix to zero qem = []; qed = []; tau = []; Me = []; Ce = []; tau = []; [nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ... ... ... @@ -52,7 +52,7 @@ nodcoord(:, 2) = nodcoord0(:, 2) + u(ii+2); % set element matrix and vector to zero ndof = length(n(1, :))*2; qe = zeros(ndof, ndof); Ke = zeros(ndof, ndof); rhse = zeros(ndof, 1); % initialize the deformation tensor f ... ... @@ -138,7 +138,7 @@ for int = 1:nint % Element stiffness matrix qe = qe + b'*(DF + Ds + DJ)*b*intw(int)*detj; Ke = Ke + b'*(DF + Ds + DJ)*b*intw(int)*detj; % right hand side array s = [sigma(1, 1) sigma(2, 2) sigma(1, 2) sigma(1, 2)]'; ... ...
 function [qem,qed,qe,rhse,tau]=elcneofib(ielem,coord,top,mat,pos,sol,dum1,tau0,dum2); function [Me,Ce,Ke,rhse,tau]=elcneofib(ielem,coord,top,mat,pos,sol,dum1,tau0,dum2); % % neohookean material behaviour ... ... @@ -9,7 +9,7 @@ function [qem,qed,qe,rhse,tau]=elcneofib(ielem,coord,top,mat,pos,sol,dum1,tau0,d % % set mass matrix and damping matrix to zero qem=[]; qed=[]; tau=[]; Me=[]; Ce=[]; tau=[]; [nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top)); ... ... @@ -53,7 +53,7 @@ nodcoord(:,2)=nodcoord0(:,2)+u(ii+2); % set element matrix and vector to zero ndof= length(n(1,:))*2; qe = zeros(ndof,ndof); Ke = zeros(ndof,ndof); rhse = zeros(ndof,1); % initialize the deformation tensor f ... ... @@ -140,7 +140,7 @@ for int=1:nint [Dfib,sigmafib]=elcneofib_D(F,mat(imat,:)); % Element stiffness matrix qe = qe + b'*(DF + Ds + DJ + Dfib)*b*intw(int)*detj; Ke = Ke + b'*(DF + Ds + DJ + Dfib)*b*intw(int)*detj; % right hand side array s=[sigma(1,1) sigma(2,2) sigma(1,2) sigma(1,2)]'; ... ...
 function [qem, qed, qe, rhse, Fhist] = ... function [Me, Ce, Ke, rhse, Fhist] = ... elcneotot(ielem, coord, top, mat, pos, sol, dum1, Fhistn, dum2) % % neohookean material behaviour ... ... @@ -9,7 +9,7 @@ function [qem, qed, qe, rhse, Fhist] = ... % % set mass matrix and damping matrix to zero qem = []; qed = []; tau = []; Me = []; Ce = []; tau = []; [nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ... ... ... @@ -52,7 +52,7 @@ nodcoord(:,2)=nodcoord0(:,2)+u(ii+2); % set element matrix and vector to zero ndof= length(n(1,:))*2; qe = zeros(ndof,ndof); Ke = zeros(ndof,ndof); rhse = zeros(ndof,1); % initialize the deformation tensor f ... ... @@ -158,7 +158,7 @@ for int=1:nint % Element stiffness matrix qe = qe + b0'*(DF + Ds + DJ)*b0*intw(int)*detj0; Ke = Ke + 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)]'; ... ...
 function [qme, qde, qe, rhse, dummy] = ele(ielem, coord, top, mat, ... varargin) % [qme, qde, qe, rhse, dummy] = ele(ielem, coord, top, mat, ... % varargin) function [Me, Ce, Ke, rhse, void] = ele(ielem, coord, top, mat, varargin) % [Me, Ce, Ke, rhse, dummy] = ele(ielem, coord, top, mat, varargin) % % linear elastic (Hookean) element % % Parameters: % % qme : element mass matrix (not used, empty) % qde : element damping matrix (not used, empty) % qe : element stiffness matrix % Me : element mass matrix (not used, empty) % Ce : element damping matrix (not used, empty) % Ke : element stiffness matrix % rhse : element right hand side (not calculated, zeros) % coord : global coordinate array for all nodes % top : topology array % top : mesh topology array % ielem : element number % mat : material properties % E = mat(imat, 1) (Young's modulus) ... ... @@ -20,13 +18,12 @@ function [qme, qde, qe, rhse, dummy] = ele(ielem, coord, top, mat, ... % axi = mat(imat, 3) (0: plane strain, 1: axi-symmetric, 2: plane stress) % eltype = mat(imat, 11) (Element type specification, see ele_i) % % Remark % mat.types : 'ele; % void = [] and varargin are for consistenty of element calls % % See also ele_s, ele_i, ele_d, femlin_e % generate the empty outputs qme = []; qde = []; dummy = []; Me = []; Ce = []; void = []; % get the number of nodes on this element nlnodes = length(nonzeros(top(ielem, 1:end-2))); ... ... @@ -49,8 +46,8 @@ etop = top(ielem, :); % number of dofs for the displacements nedofu = length(nonzeros(vpos(1:2, :))); % initialize qe and rhse qe = zeros(nedofu, nedofu); % initialize Ke and rhse Ke = zeros(nedofu, nedofu); rhse = zeros(nedofu, 1); % initialize strain-displacement matrix B = zeros(4, nedofu); ... ... @@ -114,7 +111,7 @@ for int = 1:length(intworg); B(4, ii+1) = n(int, :)/x; end qe = qe + B'*H*B*intwdetj; Ke = Ke + B'*H*B*intwdetj; end ... ...
 function [qme, qde, qe, rhse, dum3] = ... function [Me, Ce, Ke, rhse, dum3] = ... elefib(ielem, coord, top, mat, pos, sol, soln, dum1, dum2) % % linear elastic, plane stres element, including fibers % % Parameters: % % qme : element mass matrix (not used) % qde : element damping matrix (not used) % qe : element stiffness matrix % Me : element mass matrix (not used) % Ce : element damping matrix (not used) % Ke : element stiffness matrix % rhse : element right hand side % coord : coordinates of all nodes % top : topology array ... ... @@ -25,7 +25,7 @@ function [qme, qde, qe, rhse, dum3] = ... % Remark % mat.types : 'elefib'; qme = []; qde = []; dum3 = []; Me = []; Ce = []; dum3 = []; % get the number of nodes on this element [nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ... ... ... @@ -59,8 +59,8 @@ iu = zeros(1, nedofu); iu(1:2:nedofu) = nonzeros(vpos(1, :)); iu(2:2:nedofu) = nonzeros(vpos(2, :)); % initialize qe and rhs qe = zeros(nedofu, nedofu); % initialize Ke and rhs Ke = zeros(nedofu, nedofu); rhse = zeros(nedofu, 1); ... ... @@ -116,7 +116,7 @@ for int = 1:nint b(3,ii+1) = dndxy(:, 2)'; b(3,ii+2) = dndxy(:, 1)'; qe = qe + b'*d*b*intwdetj; Ke = Ke + b'*d*b*intwdetj; end ... ...
 function [qme, qde, qe, rhse, dum3] = ... function [Me, Ce, Ke, rhse, dum3] = ... eler(ielem, coord, top, mat, pos, sol, soln, dum1, dum2) % % linear elastic element, reduced integration % % Parameters: % % qme : element mass matrix (not used) % qde : element damping matrix (not used) % qe : element stiffness matrix % Me : element mass matrix (not used) % Ce : element damping matrix (not used) % Ke : element stiffness matrix % rhse : element right hand side % coord : coordinates of all nodes % top : topology array ... ... @@ -23,7 +23,7 @@ function [qme, qde, qe, rhse, dum3] = ... % Remark % mat.types : 'eler; qme = []; qde = []; dum3 = []; Me = []; Ce = []; dum3 = []; % get the number of nodes on this element [nnodes, nsd, nelem, maxnlnodes, iimat, iitype] = ... ... ... @@ -57,8 +57,8 @@ iu = zeros(1,nedofu); iu(1:2:nedofu) = nonzeros(vpos(1, :)); iu(2:2:nedofu) = nonzeros(vpos(2, :)); % initialize qe and rhs qe = zeros(nedofu, nedofu); % initialize Ke and rhs Ke = zeros(nedofu, nedofu); rhse = zeros(nedofu, 1); % compute shear modulus and compression modulus from E and Nu ... ... @@ -112,7 +112,7 @@ for int = 1:nint b(4, ii+1) = n(int, :)/x; end qe = qe + b'*dd*b*intwdetj; Ke = Ke + b'*dd*b*intwdetj; end ... ... @@ -146,7 +146,7 @@ for int = 1:nint b(4, ii+1) = n(int, :)/x; end qe = qe + b'*dv*b*intwdetj; Ke = Ke + b'*dv*b*intwdetj; end ... ...
 function [qem,qed,qe,rhse,Fhist]=elineo(ielem,coord,top,mat,pos,sol,dum1,Fhistn,dum2); function [Me,Ce,Ke,rhse,Fhist]=elineo(ielem,coord,top,mat,pos,sol,dum1,Fhistn,dum2); % % incompressible neohookean material behaviour ... ... @@ -8,7 +8,7 @@ function [qem,qed,qe,rhse,Fhist]=elineo(ielem,coord,top,mat,pos,sol,dum1,Fhistn, % % set mass matrix and damping matrix to zero qem=[]; qed=[]; Me=[]; Ce=[]; [nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top)); ... ... @@ -59,7 +59,7 @@ nodcoord(:,2)=nodcoord0(:,2)+u(ii+2); ndofu= length(n(1,:))*2; ndofp= length(np(1,:)); qe = zeros(ndofu+ndofp,ndofu+ndofp); Ke = zeros(ndofu+ndofp,ndofu+ndofp); rhse = zeros(ndofu+ndofp,1); % initialize the deformation tensor f ... ... @@ -145,12 +145,12 @@ for int=1:nint % element matrix qe(iu,iu)=qe(iu,iu)+b'*(DF+Dtau)*b*intw(int)*detj; Ke(iu,iu)=Ke(iu,iu)+b'*(DF+Dtau)*b*intw(int)*detj; % part of the stiffness matrix that belongs to the % incompressibility constraint qe(ip,iu(ii+1))=qe(ip,iu(ii+1))-np(int,:)'*(dndxy(:,1)'+axi*n(int,:)/r)*intw(int)*detj; qe(ip,iu(ii+2))=qe(ip,iu(ii+2))-np(int,:)'*dndxy(:,2)'*intw(int)*detj; Ke(ip,iu(ii+1))=Ke(ip,iu(ii+1))-np(int,:)'*(dndxy(:,1)'+axi*n(int,:)/r)*intw(int)*detj; Ke(ip,iu(ii+2))=Ke(ip,iu(ii+2))-np(int,:)'*dndxy(:,2)'*intw(int)*detj; % rhs of the incompressibility constraint rhse(ip)=rhse(ip)+ np(int,:)'*(J-1)/J*intw(int)*detj; ... ... @@ -167,6 +167,6 @@ for int=1:nint end % take care of incomp. constraint qe(iu,ip)=qe(ip,iu)'; Ke(iu,ip)=Ke(ip,iu)'; % part of mlfem_nac: https://gitlab.tue.nl/STEM/mlfem_nac
 function [qe, rhse] = elm1d(ielem, coord, top, mat) function [Ke, rhse] = elm1d(ielem, coord, top, mat) % % One dimensional diffusion problem % ... ... @@ -13,7 +13,7 @@ function [qe, rhse] = elm1d(ielem, coord, top, mat) % mat : material properties % % Output: % qe : element stiffness matrix % Ke : element stiffness matrix % rhse : element right hand side % get the number of nodes on this element ... ... @@ -39,8 +39,8 @@ etop = top(ielem, :); nint = length(intw); ndofu = length(nonzeros(vpos)); % initialize qe and rhse qe = zeros(ndofu, ndofu); % initialize Ke and rhse Ke = zeros(ndofu, ndofu); rhse = zeros(ndofu, 1); for int = 1:nint ... ... @@ -56,7 +56,7 @@ for int = 1:nint c_int = elm1d_c(c, x); % contribution to the element coefficient matrix qe = qe + b'*c_int*b*intw(int)*detj; Ke = Ke + b'*c_int*b*intw(int)*detj; % apply source term rhse = rhse + n(int, :)'*f*intw(int)*detj; ... ...
 ... ... @@ -38,8 +38,8 @@ etop=top(ielem,:); nint=length(intw); ndofu=length(nonzeros(vpos)); % initialize qe and rhs qe=zeros(ndofu,ndofu); % initialize Ke and rhs Ke=zeros(ndofu,ndofu); rhse=zeros(ndofu,1); % get the nodal solution ... ...
 function [qme, qde, qe, rhse] = elm1dcd(ielem, coord, top, mat) function [Me, Ce, Ke, rhse] = elm1dcd(ielem, coord, top, mat) % % linear convection diffusion element % ... ... @@ -12,9 +12,9 @@ function [qme, qde, qe, rhse] = elm1dcd(ielem, coord, top, mat) % mat : material properties % % Output % qme : element mass matrix % qde : element damping matrix % qe : element stifness matrix % Me : element mass matrix % Ce : element damping matrix % Ke : element stifness matrix % rhse : element right hand side % get the number of nodes on this element ... ... @@ -42,10 +42,10 @@ etop = top(ielem, :); nint = length(intw); ndofu = length(nonzeros(vpos)); % initialize qe, qme and qde and rhse qe = zeros(ndofu, ndofu); qde = qe; qme = qe; % initialize Ke, Me and Ce and rhse Ke = zeros(ndofu, ndofu); Ce = Ke; Me = Ke; rhse = zeros(ndofu, 1); if fsupg > 0 ... ... @@ -63,22 +63,22 @@ for int = 1:nint b = dndxy(:, 1)'; % diffusion part qe = qe + b'*c*b*intw(int)*detj; Ke = Ke + b'*c*b*intw(int)*detj; % convection part if fsupg > 0 qde = qde + (n(int, :)' + ... Ce = Ce + (n(int, :)' + ... h*alpha*a/(2*abs(a))*dndxy(:, 1))*a*dndxy(:, 1)'*intw(int)*detj; else qde = qde + n(int, :)'*a*dndxy(:, 1)'*intw(int)*detj; Ce = Ce + n(int, :)'*a*dndxy(:, 1)'*intw(int)*detj; end % mass matrix if fsupg > 0 qme = qme + (n(int, :)' + ... Me = Me + (n(int, :)' + ... h*alpha*a/(2*abs(a))*dndxy(:, 1))*n(int, :)*intw(int)*detj; else qme = qme + n(int, :)'*n(int, :)*intw(int)*detj; Me = Me + n(int, :)'*n(int, :)*intw(int)*detj; end % apply volume forces ... ...
 ... ... @@ -5,7 +5,7 @@ function [sigma]=elm1dcd_d(ielem,coord,top,mat,pos,sol); % Parameters: % % ielem : element number % qe : element stiffness matrix % Ke : element stiffness matrix % rhse : element right hand side % coord : coordinates of all nodes % top : topology array ... ... @@ -36,8 +36,8 @@ etop=top(ielem,:); nint=length(intw); ndofu=length(nonzeros(vpos)); % initialize qe and rhs qe=zeros(ndofu,ndofu); % initialize Ke and rhs Ke=zeros(ndofu,ndofu); rhse=zeros(ndofu,1); % get the nodal solution ... ...
 function [qme,qde,qe,rhse,dum3]=elme(ielem,coord,top,mat,pos,sol,soln,dum1,dum2); function [Me,Ce,Ke,rhse,dum3]=elme(ielem,coord,top,mat,pos,sol,soln,dum1,dum2); % % linear elastic element % % Parameters: % % qme : element mass matrix % qde : element damping matrix (not used) % qe : element stiffness matrix % Me : element mass matrix % Ce : element damping matrix (not used) % Ke : element stiffness matrix % rhse : element right hand side % coord : coordinates of all nodes % top : topology array ... ... @@ -22,7 +22,7 @@ function [qme,qde,qe,rhse,dum3]=elme(ielem,coord,top,mat,pos,sol,soln,dum1,dum2) % Remark % mat.types : 'elme'; qme=[]; qde=[]; dum3=[]; Me=[]; Ce=[]; dum3=[]; % get the number of nodes on this element [nnodes,nsd,nelem,maxnlnodes,iimat,iitype] = sysinfo(size(coord),size(top)); ... ... @@ -58,9 +58,9 @@ iu = zeros(1,nedofu); iu(1:2:nedofu) = nonzeros(vpos(1,:)); iu(2:2:nedofu) = nonzeros(vpos(2,:)); % initialize qe and rhs qe=zeros(nedofu,nedofu);