https://cfd-online.com/W/index.php?title=PFV_get_pressure&feed=atom&action=history
PFV get pressure - Revision history
2024-03-28T23:15:51Z
Revision history for this page on the wiki
MediaWiki 1.16.5
https://cfd-online.com/W/index.php?title=PFV_get_pressure&diff=17536&oldid=prev
Jonas Holdeman: Update to MatLab version R2012b and generalize to handle rectangular and triangular elements. This is a major revision.
2013-03-15T16:10:53Z
<p>Update to MatLab version R2012b and generalize to handle rectangular and triangular elements. This is a major revision.</p>
<a href="https://cfd-online.com/W/index.php?title=PFV_get_pressure&diff=17536&oldid=13021">Show changes</a>
Jonas Holdeman
https://cfd-online.com/W/index.php?title=PFV_get_pressure&diff=13021&oldid=prev
Jonas Holdeman: Created page with "Matlab function '''GetPresW.m''' to retrieve consistent pressure from velocity. <pre> function [P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu) %GETPRESW..."
2011-06-29T04:01:18Z
<p>Created page with "Matlab function '''GetPresW.m''' to retrieve consistent pressure from velocity. <pre> function [P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu) %GETPRESW..."</p>
<p><b>New page</b></p><div>Matlab function '''GetPresW.m''' to retrieve consistent pressure from velocity. <br />
<br />
<pre><br />
function [P,Px,Py] = GetPresW(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu)<br />
%GETPRESW - Compute continuous simple cubic pressure and derivatives from (simple-cubic) <br />
% velocity field on general quadrilateral grid (bilinear geometric mapping). <br />
%<br />
% Inputs<br />
% NumNod - number of nodes <br />
% NodNdx - nodal index into Xgrid and Ygrid <br />
% Elcon - element connectivity, nodes in element <br />
% nn2nft - global number and type of (non-pressure) DOF at each node <br />
% Xgrid - array of nodal x-coordinates <br />
% Ygrid - array of nodal y-coordinates <br />
% Q - array of DOFs for cubic velocity elements <br />
% EBCp - essential pressure boundary condition structure<br />
% EBCp.nodn - node number of fixed pressure node <br />
% EBCp.val - value of pressure <br />
% nu - diffusion coefficient <br />
% Outputs <br />
% P - pressure <br />
% Px - x-derivative of pressure <br />
% Py - y-derivative of pressure <br />
% Uses <br />
% ilu_gmres_with_EBC - to solve the system with essential/Dirichlet BCs <br />
% GQ3, GQ4, GQ5 - quadrature rules.<br />
<br />
% Jonas Holdeman, January 2007, revised June 2011 <br />
<br />
% Constants and fixed data<br />
nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order<br />
nnd = 4; % Number of nodes in elements<br />
nd = 3; ND=1:nd; % Number DOFs in velocity fns (bicubic-derived)<br />
np = 3; % Number DOFs in pressure fns (simple cubic)<br />
% Parameters for GMRES solver <br />
GMRES.Tolerance=1.e-9;<br />
GMRES.MaxIterates=8; <br />
GMRES.MaxRestarts=6;<br />
DropTol = 1e-7; % Drop tolerence for ilu preconditioner <br />
<br />
% Define 3-point quadrature data once, on first call (if needed). <br />
% Gaussian weights and absissas to integrate 5th degree polynomials exactly. <br />
global GQ3;<br />
if (isempty(GQ3)) % Define 3-point quadrature data once, on first call. <br />
Aq=[-.774596669241483, .000000000000000,.774596669241483]; %Abs<br />
Hq=[ .555555555555556, .888888888888889,.555555555555556]; %Wts<br />
GQ3.size=9; GQ3.xa=[Aq;Aq;Aq]; GQ3.ya=GQ3.xa'; <br />
wt=[Hq;Hq;Hq]; GQ3.wt=wt.*wt';<br />
end<br />
% Define 4-point quadrature data once, on first call (if needed). <br />
% Gaussian weights and absissas to integrate 7th degree polynomials exactly. <br />
global GQ4;<br />
if (isempty(GQ4)) % Define 4-point quadrature data once, on first call. <br />
Aq=[-.861136311594053,-.339981043584856,.339981043584856, .861136311594053]; %Abs<br />
Hq=[ .347854845137454, .652145154862546,.652145154862546, .347854845137454]; %Wts<br />
GQ4.size=16; GQ4.xa=[Aq;Aq;Aq;Aq]; GQ4.ya=GQ4.xa'; <br />
wt=[Hq;Hq;Hq;Hq]; GQ4.wt=wt.*wt'; % 4x4 <br />
end<br />
% Define 5-point quadrature data once, on first call (if needed). <br />
% Gaussian weights and absissas to integrate 9th degree polynomials exactly. <br />
global GQ5;<br />
if (isempty(GQ5)) % Has 5-point quadrature data been defined? If not, define arguments & weights. <br />
Aq=[-.906179845938664,-.538469310105683, .0, .538469310105683, .906179845938664];<br />
Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189];<br />
GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa'; <br />
wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt'; % 5x5 <br />
end<br />
% -------------- end fixed data ------------------------<br />
<br />
NumEl=size(Elcon,2); % Number of elements <br />
[NumNy,NumNx]=size(Xgrid);<br />
NumNod=NumNy*NumNx; % Number of nodes <br />
MxVDof=nd*NumNod; % Max number velocity DOFs <br />
MxPDof=np*NumNod; % Max number pressure DOFs <br />
<br />
% We can use the same nodal connectivities (Elcon) for pressure elements, <br />
% but need new pointers from nodes to pressure DOFs <br />
nn2nftp=zeros(2,NumNod); % node number -> pressure nf & nt<br />
nf=-np+1; nt=1;<br />
for n=1:NumNod<br />
nf=nf+np; % all nodes have 3 dofs <br />
nn2nftp(:,n)=[nf;nt]; % dof number & type (all nodes type 1)<br />
end;<br />
<br />
% Allocate space for pressure matrix, velocity data <br />
pMat = spalloc(MxPDof,MxPDof,30*MxPDof); % allocate P mass matrix<br />
Vdata = zeros(MxPDof,1); % allocate for velocity data (RHS)<br />
Qp=zeros(MxPDof,1); % Nodal pressure DOFs <br />
<br />
% Begin essential boundary conditions, allocate space <br />
MaxPBC = 1;<br />
EBCp.Mxdof=MxPDof;<br />
% Essential boundary condition for pressure <br />
EBCp.dof = nn2nftp(1,EBCPr.nodn(1)); % Degree-of-freedom index<br />
EBCp.val = EBCPr.val; % Pressure Dof value <br />
<br />
% partion out essential (Dirichlet) dofs<br />
p_vec = [1:EBCp.Mxdof]'; % List of all dofs<br />
EBCp.p_vec_undo = zeros(1,EBCp.Mxdof);<br />
% form a list of non-diri dofs<br />
EBCp.ndro = p_vec(~ismember(p_vec, EBCp.dof)); % list of non-diri dofs<br />
% calculate p_vec_undo to restore Q to the original dof ordering<br />
EBCp.p_vec_undo([EBCp.ndro;EBCp.dof]) = [1:EBCp.Mxdof]; %p_vec';<br />
<br />
Qp(EBCp.dof(1))=EBCp.val(1);<br />
<br />
Vdof = zeros(nd,nnd); % Nodal velocity DOFs <br />
Xe = zeros(2,nnd);<br />
<br />
% BEGIN GLOBAL MATRIX ASSEMBLY<br />
for ne=1:NumEl <br />
for k=1:4<br />
ki=NodNdx(:,Elcon(k,ne));<br />
Xe(:,k)=[Xgrid(ki(2),ki(1));Ygrid(ki(2),ki(1))]; <br />
end<br />
% Get stream function and velocities<br />
for n=1:nnd <br />
Vdof(ND,n)=Q((nn2nft(1,Elcon(n,ne))-1)+ND); % Loop over local nodes of element<br />
end<br />
[pMmat,Rndx,Cndx] = pMassMat(Xe,Elcon(:,ne),nn2nftp); % Pressure "mass" matrix <br />
pMat=pMat+sparse(Rndx,Cndx,pMmat,MxPDof,MxPDof); % Global mass assembly <br />
<br />
[CDat,RRndx] = PCDat(Xe,Elcon(:,ne),nn2nftp,Vdof); % Convective data term<br />
Vdata([RRndx]) = Vdata([RRndx])-CDat(:);<br />
<br />
[DDat,RRndx] = PDDatL(Xe,Elcon(:,ne),nn2nftp,Vdof); % Diffusive data term <br />
Vdata([RRndx]) = Vdata([RRndx]) + nu*DDat(:); % +nu*DDat;<br />
end; % Loop ne <br />
% END GLOBAL MATRIX ASSEMBLY<br />
<br />
Vdatap=Vdata(EBCp.ndro)-pMat(EBCp.ndro,EBCp.dof)*EBCp.val;<br />
pMatr=pMat(EBCp.ndro,EBCp.ndro);<br />
Qs=Qp(EBCp.ndro); % Non-Dirichlet (active) dofs <br />
<br />
Pr=ilu_gmres_with_EBC(pMatr,Vdatap,[],GMRES,Qs,DropTol);<br />
<br />
Qp=[Pr;EBCp.val]; % Augment active dofs with esential (Dirichlet) dofs<br />
Qp=Qp(EBCp.p_vec_undo); % Restore natural order<br />
Qp=reshape(Qp,np,NumNod);<br />
P= reshape(Qp(1,:),NumNy,NumNx); <br />
Px=reshape(Qp(2,:),NumNy,NumNx); <br />
Py=reshape(Qp(3,:),NumNy,NumNx); <br />
return;<br />
% >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<<<br />
<br />
% -------------------- function pMassMat ----------------------------<br />
<br />
function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp)<br />
% Simple cubic gradient element "mass" matrix <br />
% -------------- Constants and fixed data -----------------<br />
global GQ4;<br />
xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size;<br />
nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order<br />
nnd=4; <br />
np=3; np4=nnd*np; NP=1:np; <br />
%<br />
global ZG3412pm; <br />
if (isempty(ZG3412pm)|size(ZG3412pm,2)~=Nq)<br />
% Evaluate and save/cache the set of shape functions at quadrature pts. <br />
ZG3412pm=cell(nnd,Nq); <br />
for k=1:Nq<br />
for m=1:nnd<br />
ZG3412pm{m,k}=Gr(nn(m,:),xa(k),ya(k));<br />
end<br />
end<br />
end % if(isempty(*))<br />
% --------------------- end fixed data -----------------<br />
<br />
TGi=cell(nnd);<br />
for m=1:nnd % Loop over corner nodes <br />
J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))'; % GBL is gradient of bilinear function<br />
TGi{m} = blkdiag(1,J);<br />
end % Loop m <br />
<br />
MM=zeros(np4,np4); G=zeros(2,np4); % Preallocate arrays<br />
for k=1:Nq <br />
% Initialize functions and derivatives at the quadrature point (xa,ya).<br />
J=(Xe*GBL(nn(:,:),xa(k),ya(k)))'; % transpose of Jacobian J at (xa,ya)<br />
Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of J <br />
Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det; % inverse of J<br />
<br />
mm = 0;<br />
for m=1:nnd <br />
G(:,mm+NP) = Ji*ZG3412pm{m,k}*TGi{m};<br />
mm = mm+np;<br />
end % loop m<br />
MM = MM + G'*G*(wt(k)*Det);<br />
end % end loop k (quadrature pts)<br />
<br />
gf=zeros(np4,1); % array of global freedoms <br />
m=0; <br />
for n=1:4 % Loop over element nodes <br />
gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; % Get global freedoms<br />
m=m+np;<br />
end<br />
<br />
Rndx=repmat(gf,1,np4); % Row indices<br />
Cndx=Rndx'; % Column indices<br />
<br />
MM = reshape(MM,1,np4*np4);<br />
Rndx=reshape(Rndx,1,np4*np4);<br />
Cndx=reshape(Cndx,1,np4*np4); <br />
return;<br />
<br />
% --------------------- funnction PCDat -----------------------<br />
<br />
function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof)<br />
% Evaluate convective force data <br />
% ----------- Constants and fixed data ---------------<br />
global GQ5;<br />
xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size;<br />
nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order<br />
nnd=4; % number of nodes <br />
np = 3; np4=4*np; NP=1:np;<br />
nd = 3; nd4=4*nd; ND=1:nd;<br />
%<br />
global ZS3412pc; global ZSX3412pc; global ZSY3412pc; global ZG3412pc; <br />
if (isempty(ZS3412pc)|size(ZS3412pc,2)~=Nq)<br />
% Evaluate and save/cache the set of shape functions at quadrature pts. <br />
ZS3412pc=cell(nnd,Nq); ZSX3412pc=cell(nnd,Nq); <br />
ZSY3412pc=cell(nnd,Nq); ZG3412pc=cell(nnd,Nq); <br />
for k=1:Nq<br />
for m=1:nnd<br />
ZS3412pc{m,k} =Sr(nn(m,:),xa(k),ya(k));<br />
ZSX3412pc{m,k}=Sxr(nn(m,:),xa(k),ya(k));<br />
ZSY3412pc{m,k}=Syr(nn(m,:),xa(k),ya(k));<br />
ZG3412pc{m,k}=Gr(nn(m,:),xa(k),ya(k));<br />
end % loop m over nodes <br />
end % loop k over Nq<br />
end % if(isempty(*))<br />
% ----------------- end fixed data ------------------<br />
<br />
Ti=cell(nnd); TGi=cell(nnd);<br />
for m=1:nnd % Loop over corner nodes <br />
J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))'; % Jacobian at (xa,ya)<br />
Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of J & Jt <br />
JiD=[J(2,2),-J(1,2); -J(2,1),J(1,1)]; % inv(J)*det(J)<br />
Ti{m} = blkdiag(1,JiD');<br />
TGi{m} = blkdiag(1,J);<br />
end % Loop m over corner nodes<br />
<br />
PC=zeros(np4,1);<br />
S=zeros(2,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4); G=zeros(2,np4);<br />
<br />
for k=1:Nq % Loop over quadrature points <br />
J=(Xe*GBL(nn(:,:),xa(k),ya(k)))'; % Jacobian at (xa,ya)<br />
Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of J & Jt <br />
Jtbd=(J/Det)'; % transpose(J/D)<br />
JiD=[J(2,2),-J(1,2); -J(2,1),J(1,1)]; % inv(J)*det(J)<br />
Ji=JiD/Det; % inverse(J)<br />
for m=1:4 % Loop over element nodes <br />
mm=nd*(m-1);<br />
S(:,mm+ND) =Jtbd*ZS3412pc{m,k}*Ti{m};<br />
Sx(:,mm+ND)=Jtbd*(Ji(1,1)*ZSX3412pc{m,k}+Ji(1,2)*ZSY3412pc{m,k})*Ti{m};<br />
Sy(:,mm+ND)=Jtbd*(Ji(2,1)*ZSX3412pc{m,k}+Ji(2,2)*ZSY3412pc{m,k})*Ti{m};<br />
mm=np*(m-1);<br />
G(:,mm+NP)=Ji*ZG3412pc{m,k}*TGi{m};<br />
end % end loop over element nodes<br />
<br />
% Compute the fluid velocities at the quadrature point.<br />
U = S*Vdof(:);<br />
Ux = Sx*Vdof(:);<br />
Uy = Sy*Vdof(:);<br />
UgU = U(1)*Ux+U(2)*Uy; % U dot grad U <br />
PC = PC + G'*UgU*(wt(k)*Det);<br />
end % end loop over Nq quadrature points<br />
<br />
gf=zeros(1,np4); % array of global freedoms <br />
m=0; <br />
for n=1:4 % Loop over element nodes <br />
gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; % Get global freedoms<br />
m=m+np;<br />
end<br />
Rndx=gf;<br />
PC = reshape(PC,1,np4);<br />
return;<br />
<br />
% ----------------- function PDDatL -------------------------<br />
<br />
function [PD,Rndx]=PDDatL(Xe,Elcon,nn2nftp,Vdof)<br />
% Evaluate diffusive force data (Laplacian form) <br />
% --------------- Constants and fixed data -------------------<br />
global GQ3;<br />
xa=GQ3.xa; ya=GQ3.ya; wt=GQ3.wt; Nq=GQ3.size;<br />
nnd=4; % number of nodes <br />
nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order<br />
np = 3; npdf=nnd*np; NP=1:np;<br />
nd = 3; nd4=nnd*nd; ND=1:nd;<br />
global ZSXX3412pd; global ZSXY3412pd; global ZSYY3412pd; global ZG3412pd; <br />
if (isempty(ZSXX3412pd)|size(ZSXX3412pd,2)~=Nq)<br />
% Evaluate and save/cache the set of shape functions at quadrature pts. <br />
ZSXX3412pd=cell(nnd,Nq); ZSXY3412pd=cell(nnd,Nq); <br />
ZSYY3412pd=cell(nnd,Nq); ZG3412pd=cell(nnd,Nq);<br />
global ZSYY3412pd;<br />
for k=1:Nq<br />
for m=1:nnd<br />
ZSXX3412pd{m,k}=Sxxr(nn(m,:),xa(k),ya(k));<br />
ZSXY3412pd{m,k}=Sxyr(nn(m,:),xa(k),ya(k));<br />
ZSYY3412pd{m,k}=Syyr(nn(m,:),xa(k),ya(k));<br />
ZG3412pd{m,k}=Gr(nn(m,:),xa(k),ya(k));<br />
end % loop m over nodes <br />
end % loop k over Nq<br />
end % if(isempty(*))<br />
% ------------ end fixed data -------------------<br />
%<br />
Ti=cell(nnd); TGi=cell(nnd);<br />
for m=1:nnd % Loop over corner nodes <br />
J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))'; % Jacobian at (xa,ya)<br />
Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of J & Jt <br />
JiD=[J(2,2),-J(1,2); -J(2,1),J(1,1)]; % inv(J)*det(J)<br />
Ti{m} = blkdiag(1,JiD');<br />
TGi{m} = blkdiag(1,J);<br />
end % Loop m over corner nodes<br />
<br />
PD=zeros(npdf,1);<br />
Sxx=zeros(2,nd4); Syy=zeros(2,nd4); G=zeros(2,npdf);<br />
for k=1:Nq % Loop over quadrature points <br />
Jt=(Xe*GBL(nn(:,:),xa(k),ya(k))); % transpose of Jacobian at (xa,ya)<br />
Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J <br />
Jtd=Jt/Det;<br />
JiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)];<br />
Ji=JiD/Det;<br />
for m=1:nnd % Loop over element nodes <br />
mm=nd*(m-1); % This transform is approximate !!<br />
Sxx(:,mm+ND)=Jtd*(Ji(1,1)^2*ZSXX3412pd{m,k}+2*Ji(1,1)*Ji(1,2)*ZSXY3412pd{m,k}+Ji(1,2)^2*ZSYY3412pd{m,k})*Ti{m};<br />
Syy(:,mm+ND)=Jtd*(Ji(2,1)^2*ZSXX3412pd{m,k}+2*Ji(2,1)*Ji(2,2)*ZSXY3412pd{m,k}+Ji(2,2)^2*ZSYY3412pd{m,k})*Ti{m};<br />
mm=np*(m-1);<br />
G(:,mm+NP) =Ji*ZG3412pd{m,k}*TGi{m};<br />
end % end loop over element nodes<br />
<br />
LapU = (Sxx+Syy)*Vdof(:);<br />
PD = PD+G'*LapU*(wt(k)*Det);<br />
end % end loop over quadrature points<br />
<br />
gf=zeros(1,npdf); % array of global freedoms <br />
m=0; K=1:np;<br />
for n=1:nnd % Loop over element nodes <br />
nfn=nn2nftp(1,Elcon(n))-1; % Get global freedom <br />
gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP;<br />
m=m+np;<br />
end<br />
Rndx=gf;<br />
PD = reshape(PD,1,npdf);<br />
return;<br />
<br />
% ------------------------------------------------------------------------------<br />
function gv=Gr(ni,q,r)<br />
%GR Cubic Hermite gradient basis functions for pressure gradient.<br />
qi=ni(1); q0=q*ni(1); <br />
ri=ni(2); r0=r*ni(2); <br />
% Cubic Hermite gradient <br />
gv=[1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), -1/8*(1+r0)*(1+q0)*(1-3*q0), ...<br />
-1/8*qi/ri*(1-r^2)*(1+r0); ...<br />
1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), -1/8/qi*ri*(1-q^2)*(1+q0), ...<br />
-1/8*(1+q0)*(1+r0)*(1-3*r0)];<br />
return;<br />
<br />
function gx=Gxr(ni,q,r)<br />
%GRX - Cubic Hermite gradient basis functions for pressure gradient.<br />
qi=ni(1); q0=q*ni(1); <br />
ri=ni(2); r0=r*ni(2); <br />
% x-derivative of irrotational vector <br />
gx=[-3/4*qi^2*q0*(1+r0), 1/4*qi*(1+r0)*(1+3*q0), 0; ...<br />
1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*ri*(1+q0)*(1-3*q0), -1/8*qi*(1+r0)*(1-3*r0)];<br />
return;<br />
<br />
function gy=Gyr(ni,q,r)<br />
%GRY - Cubic Hermite gradient basis functions for pressure gradient.<br />
qi=ni(1); q0=q*ni(1); <br />
ri=ni(2); r0=r*ni(2); <br />
% y-derivative of irrotational vector <br />
gy=[1/8*qi*ri*(4-3*q0^2-3*r0^2), -1/8*ri*(1+q0)*(1-3*q0), -1/8*qi*(1+r0)*(1-3*r0); ...<br />
-3/4*ri^2*r0*(1+q0), 0, 1/4*ri*(1+q0)*(1+3*r0)];<br />
return;<br />
<br />
% ------------------------------------------------------------------------------<br />
function S=Sr(ni,q,r)<br />
%SR Array of solenoidal basis functions.<br />
qi=ni(1); q0=q*ni(1); q1=1+q0; <br />
ri=ni(2); r0=r*ni(2); r1=1+r0;<br />
% array of solenoidal vectors <br />
S=[ .125*ri*q1*(q0*(1-q0)+3*(1-r^2)), .125*q1*r1*(3*r0-1), .125*ri/qi*q1^2*(1-q0); ...<br />
-.125*qi*r1*(r0*(1-r0)+3*(1-q^2)), .125*qi/ri*r1^2*(1-r0), .125*q1*r1*(3*q0-1)];<br />
return;<br />
<br />
function S=Sxr(ni,q,r)<br />
%SXR Array of x-derivatives of solenoidal basis functions.<br />
qi=ni(1); q0=q*ni(1); q1=1+q0; <br />
ri=ni(2); r0=r*ni(2); r1=1+r0;<br />
% array of solenoidal vectors <br />
S=[.125*qi*ri*(4-3*q^2-3*r^2), .125*qi*r1*(3*r0-1), -.125*ri*q1*(3*q0-1); ...<br />
.75*qi^2*r1*q0, 0, .25*qi*r1*(3*q0+1)];<br />
return;<br />
<br />
function s=Syr(ni,q,r)<br />
%SYR Array of y-derivatives of solenoidal basis functions.<br />
qi=ni(1); q0=q*ni(1); q1=1+q0; <br />
ri=ni(2); r0=r*ni(2); r1=1+r0;<br />
% array of solenoidal vectors <br />
s=[-.75*ri^2*q1*r0, .25*ri*q1*(3*r0+1), 0 ; ...<br />
-.125*qi*ri*(4-3*q^2-3*r^2), .125*qi*r1*(1-3*r0), .125*ri*q1*(3*q0-1)];<br />
return;<br />
<br />
function s=Sxxr(ni,q,r)<br />
%SXXR Array of second x-derivatives of solenoidal basis functions.<br />
qi=ni(1); q0=q*ni(1); q1=1+q0; <br />
ri=ni(2); r0=r*ni(2); r1=1+r0;<br />
% array of solenoidal vectors <br />
s=[-3/4*qi^2*ri*q0, 0, -1/4*ri*qi*(1+3*q0); 3/4*qi^3*r1, 0, 3/4*qi^2*r1 ];<br />
return;<br />
<br />
function s=Syyr(ni,q,r)<br />
%SYYR Array of second y-derivatives of solenoidal basis functions.<br />
qi=ni(1); q0=q*ni(1); q1=1+q0; <br />
ri=ni(2); r0=r*ni(2); r1=1+r0;<br />
% array of solenoidal vectors <br />
s=[-3/4*ri^3*q1, 3/4*ri^2*q1, 0; 3/4*qi*ri^2*r0, -1/4*qi*ri*(1+3*r0), 0 ];<br />
return;<br />
<br />
function s=Sxyr(ni,q,r)<br />
%SXYR Array of second (cross) xy-derivatives of solenoidal basis functions.<br />
qi=ni(1); q0=q*ni(1); q1=1+q0; <br />
ri=ni(2); r0=r*ni(2); r1=1+r0;<br />
% array of solenoidal vectors <br />
s=[-3/4*qi*ri^2*r0, 1/4*qi*ri*(1+3*r0), 0; 3/4*qi^2*ri*q0, 0, 1/4*qi*ri*(1+3*q0)];<br />
return;<br />
<br />
function G=GBL(ni,q,r)<br />
% Transposed gradient (derivatives) of scalar bilinear mapping function. <br />
% The parameter ni can be a vector of coordinate pairs. <br />
G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)];<br />
return;<br />
</pre></div>
Jonas Holdeman