PFV get pressure 2
From CFD-Wiki
Matlab function GetPres4424W.m to retrieve consistent pressure from velocity.
function [P,Px,Py] = GetPres44243W(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBCPr,nu) %GETPRES44243G - Compute continuous simple cubic pressure from (quartic-based) % velocity field on general quadrilateral grid (bilinear geometric mapping). % % Usage % [P,PX,PY]=Getpres44243W(NumNod,NodNdx,Elcon,nn2nft,Xgrid,Ygrid,Q,EBDp,nu); % Inputs % NumNod - number of nodes % NodNdx - nodal index into Xgrid and Ygrid % Elcon - element connectivity, nodes in element % nn2nft - global number and type of (non-pressure) DOF at each node % Xgrid - array of nodal x-coordinates % Ygrid - array of nodal y-coordinates % Q - array of DOFs for cubic velocity elements % EBCp - essential pressure boundary condition structure % EBCp.nodn - node number of fixed pressure node % EBCp.val - value of pressure % nu - diffusion coefficient % Outputs % P - pressure % Px - x-derivative of pressure % Py - y-derivative of pressure % Uses % ilu_gmres_with_EBC - to solve the system with essential/Dirichlet BCs % GQ3, GQ4, GQ5 - quadrature rules. % Jonas Holdeman, January 2007, revised July 2011 % Constants and fixed data nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order nnd = 4; % Number of nodes in elements nd = 6; ND=1:nd; % Number DOFs in velocity fns (bicubic-derived) np = 3; % Number DOFs in pressure fns (simple cubic) % Parameters for GMRES solver GMRES.Tolerance=1.e-9; GMRES.MaxIterates=8; GMRES.MaxRestarts=6; DropTol=1e-7; % Drop tolerence for ilu preconditioner % Define 3-point quadrature data once, on first call (if needed). % Gaussian weights and absissas to integrate 5th degree polynomials exactly. global GQ3; if (isempty(GQ3)) % Define 3-point quadrature data once, on first call. Aq=[-.774596669241483, .000000000000000,.774596669241483]; %Abs Hq=[ .555555555555556, .888888888888889,.555555555555556]; %Wts GQ3.size=9; GQ3.xa=[Aq;Aq;Aq]; GQ3.ya=GQ3.xa'; wt=[Hq;Hq;Hq]; GQ3.wt=wt.*wt'; end % Define 4-point quadrature data once, on first call (if needed). % Gaussian weights and absissas to integrate 7th degree polynomials exactly. global GQ4; if (isempty(GQ4)) % Define 4-point quadrature data once, on first call. Aq=[-.861136311594053,-.339981043584856,.339981043584856, .861136311594053]; %Abs Hq=[ .347854845137454, .652145154862546,.652145154862546, .347854845137454]; %Wts GQ4.size=16; GQ4.xa=[Aq;Aq;Aq;Aq]; GQ4.ya=GQ4.xa'; wt=[Hq;Hq;Hq;Hq]; GQ4.wt=wt.*wt'; % 4x4 end % Define 5-point quadrature data once, on first call (if needed). % Gaussian weights and absissas to integrate 9th degree polynomials exactly. global GQ5; if (isempty(GQ5)) % Has 5-point quadrature data been defined? If not, define arguments & weights. Aq=[-.906179845938664,-.538469310105683, .0, .538469310105683, .906179845938664]; Hq=[ .236926885056189, .478628670499366, .568888888888889, .478628670499366, .236926885056189]; GQ5.size=25; GQ5.xa=[Aq;Aq;Aq;Aq;Aq]; GQ5.ya=GQ5.xa'; wt=[Hq;Hq;Hq;Hq;Hq]; GQ5.wt=wt.*wt'; % 5x5 end % Set IsAffine=1 for affine elements, IsAffine=0 for general quad IsAffine=0; % -------------- end fixed data ------------------------ NumEl=size(Elcon,2); % Number of elements [NumNy,NumNx]=size(Xgrid); NumNod=NumNy*NumNx; % Number of nodes MxVDof=nd*NumNod; % Max number velocity DOFs MxPDof=np*NumNod; % Max number pressure DOFs % We can use the same nodal connectivities (Elcon) for pressure elements, % but need new pointers from nodes to pressure DOFs nn2nftp=zeros(2,NumNod); % node number -> pressure nf & nt nf=-np+1; nt=1; for n=1:NumNod nf=nf+np; % all nodes have 3 dofs nn2nftp(:,n)=[nf;nt]; % dof number & type (all nodes type 1) end; % Allocate space for pressure matrix, velocity data pMat = spalloc(MxPDof,MxPDof,36*MxPDof); % allocate P mass matrix Vdata = zeros(MxPDof,1); % allocate velocity data Qp=zeros(MxPDof,1); % Nodal pressure DOFs % Begin essential boundary conditions, allocate space MaxPBC = 1; EBCp.Mxdof=MxPDof; % Essential boundary condition for pressure EBCp.dof = nn2nftp(1,EBCPr.nodn(1)); % Degree-of-freedom index EBCp.val = EBCPr.val; % Pressure Dof value % partion out essential (Dirichlet) dofs p_vec = [1:EBCp.Mxdof]'; % List of all dofs EBCp.p_vec_undo = zeros(1,EBCp.Mxdof); % form a list of non-diri dofs EBCp.ndro = p_vec(~ismember(p_vec, EBCp.dof)); % list of non-diri dofs % calculate p_vec_undo to restore Q to the original dof ordering EBCp.p_vec_undo([EBCp.ndro;EBCp.dof]) = [1:EBCp.Mxdof]; %p_vec'; Qp(EBCp.dof(1))=EBCp.val(1); Vdof = zeros(nd,nnd); % Nodal velocity DOFs Xe = zeros(2,nnd); % BEGIN GLOBAL MATRIX ASSEMBLY for ne=1:NumEl for k=1:nnd ki=NodNdx(:,Elcon(k,ne)); Xe(:,k)=[Xgrid(ki(2),ki(1));Ygrid(ki(2),ki(1))]; end % Get stream function and velocities for n=1:nnd Vdof(ND,n)=Q((nn2nft(1,Elcon(n,ne))-1)+ND); % Loop over local nodes of element end [pMmat,Rndx,Cndx] = pMassMat(Xe,Elcon(:,ne),nn2nftp); % Pressure "mass" matrix pMat=pMat+sparse(Rndx,Cndx,pMmat,MxPDof,MxPDof); % Global mass assembly [CDat,RRndx] = PCDat(Xe,Elcon(:,ne),nn2nftp,Vdof); % Convective data term Vdata([RRndx]) = Vdata([RRndx])-CDat(:); [DDat,RRndx] = PDDatL(Xe,Elcon(:,ne),nn2nftp,Vdof); % Diffusive data term Vdata([RRndx]) = Vdata([RRndx]) + nu*DDat(:); % +nu*DDat; end; % Loop ne % END GLOBAL MATRIX ASSEMBLY Vdatap=Vdata(EBCp.ndro)-pMat(EBCp.ndro,EBCp.dof)*EBCp.val; pMatr=pMat(EBCp.ndro,EBCp.ndro); Qs=Qp(EBCp.ndro); % Non-Dirichlet (active) dofs Pr=ilu_gmres_with_EBC(pMatr,Vdatap,[],GMRES,Qs,DropTol); Qp=[Pr;EBCp.val]; % Augment active dofs with esential (Dirichlet) dofs Qp=Qp(EBCp.p_vec_undo); % Restore natural order Qp=reshape(Qp,np,NumNod); P= reshape(Qp(1,:),NumNy,NumNx); Px=reshape(Qp(2,:),NumNy,NumNx); Py=reshape(Qp(3,:),NumNy,NumNx); return; % >>>>>>>>>>>>> End pressure recovery <<<<<<<<<<<<< % -------------------- function pMassMat ---------------------------- function [MM,Rndx,Cndx]=pMassMat(Xe,Elcon,nn2nftp) % Simple cubic gradient element "mass" matrix % % -------------- Constants and fixed data ----------------- global GQ4; xa=GQ4.xa; ya=GQ4.ya; wt=GQ4.wt; Nq=GQ4.size; nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order nnd=4; np=3; np4=nnd*np; NP=1:np; % global ZG3412pm; if (isempty(ZG3412pm)|size(ZG3412pm,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZG3412pm=cell(nnd,Nq); for k=1:Nq for m=1:nnd ZG3412pm{m,k}=Gr(nn(m,:),xa(k),ya(k)); end end end % if(isempty(*)) % --------------------- end fixed data ----------------- TGi=cell(nnd); for m=1:nnd % Loop over corner nodes J=(Xe*GBL(nn(:,:),nn(m,1),nn(m,2)))'; % GBL is gradient of bilinear function TGi{m} = blkdiag(1,J); end % Loop m MM=zeros(np4,np4); G=zeros(2,np4); % Preallocate arrays for k=1:Nq % Initialize functions and derivatives at the quadrature point (xa,ya). J=(Xe*GBL(nn(:,:),xa(k),ya(k)))'; % Jacobian at (xa,ya) Det=J(1,1)*J(2,2)-J(1,2)*J(2,1); % Determinant of Jt & J Ji=[J(2,2),-J(1,2); -J(2,1),J(1,1)]/Det; mm = 0; for m=1:nnd G(:,mm+NP)=Ji*ZG3412pm{m,k}*TGi{m}; mm = mm+np; end % loop m MM = MM + G'*G*(wt(k)*Det); end % end loop k (quadrature pts) gf=zeros(np4,1); % array of global freedoms m=0; for n=1:4 % Loop over element nodes gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; % Get global freedoms m=m+np; end Rndx=repmat(gf,1,np4); % Row indices Cndx=Rndx'; % Column indices MM = reshape(MM,1,np4*np4); Rndx=reshape(Rndx,1,np4*np4); Cndx=reshape(Cndx,1,np4*np4); return; % --------------------- funnction PCDat ----------------------- function [PC,Rndx]=PCDat(Xe,Elcon,nn2nftp,Vdof) % Evaluate convective force data % % ----------- Constants and fixed data --------------- global GQ5; xa=GQ5.xa; ya=GQ5.ya; wt=GQ5.wt; Nq=GQ5.size; nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order nnd=4; % number of nodes np = 3; np4=4*np; NP=1:np; nd = 6; nd4=4*nd; ND=1:nd; % global ZS4424pc; global ZSX4424pc; global ZSY4424pc; global ZG4424pc; if (isempty(ZS4424pc)|size(ZS4424pc,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZS4424pc=cell(nnd,Nq); ZSX4424pc=cell(nnd,Nq); ZSY4424pc=cell(nnd,Nq); ZG4424pc=cell(nnd,Nq); for k=1:Nq for m=1:nnd ZS4424pc{m,k} =Sr(nn(m,:),xa(k),ya(k)); ZSX4424pc{m,k}=Sxr(nn(m,:),xa(k),ya(k)); ZSY4424pc{m,k}=Syr(nn(m,:),xa(k),ya(k)); ZG4424pc{m,k} =Gr(nn(m,:),xa(k),ya(k)); end % loop m over nodes end % loop k over Nq end % if(isempty(*)) % ----------------- end fixed data ------------------ Ti=cell(nnd); TGi=cell(nnd); for m=1:nnd % Loop over corner nodes Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2)); JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ... % alt Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ... Jt(2,1)^2, 2*Jt(2,1)*Jt(2,2), Jt(2,2)^2]; T6=blkdiag(1,JtiD,T2); Bxy=Xe*BLxy(nn(:,:),nn(m,1),nn(m,2)); % Second cross derivatives T6(5,2:3)=Bxy([2,1])'; Ti{m}=T6; TGi{m} = blkdiag(1,Jt'); end % Loop m over corner nodes PC=zeros(np4,1); S=zeros(2,nd4); Sx=zeros(2,nd4); Sy=zeros(2,nd4); G=zeros(2,np4); for k=1:Nq % Loop over quadrature points Jt=Xe*GBL(nn(:,:),xa(k),ya(k)); % transpose of Jacobian at (xa,ya) Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J Jtd=Jt/Det; JiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; Ji=JiD/Det; for m=1:nnd % Loop over element nodes mm=nd*(m-1); S(:,mm+ND) =Jtd*ZS4424pc{m,k}*Ti{m}; Sx(:,mm+ND)=Jtd*(Ji(1,1)*ZSX4424pc{m,k}+Ji(2,1)*ZSY4424pc{m,k})*Ti{m}; Sy(:,mm+ND)=Jtd*(Ji(1,2)*ZSX4424pc{m,k}+Ji(2,2)*ZSY4424pc{m,k})*Ti{m}; mm=np*(m-1); G(:,mm+NP)=Ji*ZG4424pc{m,k}*TGi{m}; end % end loop over element nodes % Compute the fluid velocities at the quadrature point. U = S*Vdof(:); Ux = Sx*Vdof(:); Uy = Sy*Vdof(:); UgU = U(1)*Ux+U(2)*Uy; % U dot grad U PC = PC + G'*UgU*(wt(k)*Det); end % end loop over Nq quadrature points gf=zeros(1,np4); % array of global freedoms m=0; for n=1:4 % Loop over element nodes gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; % Get global freedoms m=m+np; end Rndx=gf; PC = reshape(PC,1,np4); return; % ----------------- function PDDatL ------------------------- function [PD,Rndx]=PDDatL(Xe,Elcon,nn2nftp,Vdof) % Evaluate diffusive force data (Laplacian form) % --------------- Constants and fixed data ------------------- global GQ3; xa=GQ3.xa; ya=GQ3.ya; wt=GQ3.wt; Nq=GQ3.size; nnd=4; % number of nodes nn=[-1 -1; 1 -1; 1 1; -1 1]; % defines local nodal order np = 3; npdf=nnd*np; NP=1:np; nd = 6; nd4=nnd*nd; ND=1:nd; global ZSXX4424pd; global ZSXY4424pd; global ZSYY4424pd; global ZG4424pd; if (isempty(ZSXX4424pd)|size(ZSXX4424pd,2)~=Nq) % Evaluate and save/cache the set of shape functions at quadrature pts. ZSXX4424pd=cell(nnd,Nq); ZSXY4424pd=cell(nnd,Nq); ZSYY4424pd=cell(nnd,Nq); ZG4424pd=cell(nnd,Nq); for k=1:Nq for m=1:nnd ZSXX4424pd{m,k}=Sxxr(nn(m,:),xa(k),ya(k)); ZSXY4424pd{m,k}=Sxyr(nn(m,:),xa(k),ya(k)); ZSYY4424pd{m,k}=Syyr(nn(m,:),xa(k),ya(k)); ZG4424pd{m,k}=Gr(nn(m,:),xa(k),ya(k)); end % loop m over nodes end % loop k over Nq end % if(isempty(*)) % ------------ end fixed data ------------------- % Ti=cell(nnd); TGi=cell(nnd); for m=1:nnd % Loop over corner nodes Jt=Xe*GBL(nn(:,:),nn(m,1),nn(m,2)); JtiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; % det(J)*inv(J) T2=[Jt(1,1)^2, 2*Jt(1,1)*Jt(1,2), Jt(1,2)^2; ... % alt Jt(1,1)*Jt(2,1), Jt(1,1)*Jt(2,2)+Jt(1,2)*Jt(2,1), Jt(1,2)*Jt(2,2); ... Jt(2,1)^2, 2*Jt(2,1)*Jt(2,2), Jt(2,2)^2]; T6=blkdiag(1,JtiD,T2); Bxy=Xe*BLxy(nn(:,:),nn(m,1),nn(m,2)); % Second cross derivatives T6(5,2:3)=Bxy([2,1])'; Ti{m}=T6; TGi{m} = blkdiag(1,Jt'); end % Loop m over corner nodes PD=zeros(npdf,1); Sxx=zeros(2,nd4); Syy=zeros(2,nd4); G=zeros(2,npdf); for k=1:Nq % Loop over quadrature points Jt=Xe*GBL(nn(:,:),xa(k),ya(k)); % transpose of Jacobian at (xa,ya) Det=Jt(1,1)*Jt(2,2)-Jt(1,2)*Jt(2,1); % Determinant of Jt & J Jtd=Jt/Det; JiD=[Jt(2,2),-Jt(1,2); -Jt(2,1),Jt(1,1)]; Ji=JiD/Det; for m=1:nnd % Loop over element nodes mm=nd*(m-1); % This transform is approximate !! % Sxx(:,mm+ND)=Ji(1,1)^2*Jtd*ZSXX4424pd{m,k}*Ti{m}; % Syy(:,mm+ND)=Ji(2,2)^2*Jtd*ZSYY4424pd{m,k}*Ti{m}; Sxx(:,mm+ND)=Jtd*(Ji(1,1)^2*ZSXX4424pd{m,k}+2*Ji(1,1)*Ji(1,2)*ZSXY4424pd{m,k}+Ji(1,2)^2*ZSYY4424pd{m,k})*Ti{m}; Syy(:,mm+ND)=Jtd*(Ji(2,1)^2*ZSXX4424pd{m,k}+2*Ji(2,1)*Ji(2,2)*ZSXY4424pd{m,k}+Ji(2,2)^2*ZSYY4424pd{m,k})*Ti{m}; mm=np*(m-1); G(:,mm+NP) =Ji*ZG4424pd{m,k}*TGi{m}; end % end loop over element nodes LapU = (Sxx+Syy)*Vdof(:); PD = PD+G'*LapU*(wt(k)*Det); end % end loop over quadrature points gf=zeros(1,npdf); % array of global freedoms m=0; K=1:np; for n=1:nnd % Loop over element nodes nfn=nn2nftp(1,Elcon(n))-1; % Get global freedom gf(m+NP)=(nn2nftp(1,Elcon(n))-1)+NP; m=m+np; end Rndx=gf; PD = reshape(PD,1,npdf); return; % ------------------------------------------------------------------------------ function gv=Gr(ni,q,r) %GR Gradient of cubic Hermite basis functions for pressuret. qi=ni(1); q0=q*ni(1); ri=ni(2); r0=r*ni(2); % Cubic Hermite gradient gv=[1/8*qi*(1+r0)*(r0*(1-r0)+3*(1-q^2)), -1/8*(1+r0)*(1+q0)*(1-3*q0), ... -1/8*qi/ri*(1-r^2)*(1+r0); ... 1/8*ri*(1+q0)*(q0*(1-q0)+3*(1-r^2)), -1/8/qi*ri*(1-q^2)*(1+q0), ... -1/8*(1+q0)*(1+r0)*(1-3*r0)]; return; function gx=Gxr(ni,q,r) %GRX - x-derivative of cubic Hermite gradient basis functions for pressure. qi=ni(1); q0=q*ni(1); ri=ni(2); r0=r*ni(2); % x-derivative of irrotational vector gx=[-3/4*qi^2*q0*(1+r0), 1/4*qi*(1+r0)*(1+3*q0), 0; ... 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)]; return; function gy=Gyr(ni,q,r) %GRY - y-derivative of cubic Hermite gradient basis functions for pressuret. qi=ni(1); q0=q*ni(1); ri=ni(2); r0=r*ni(2); % y-derivative of irrotational vector 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); ... -3/4*ri^2*r0*(1+q0), 0, 1/4*ri*(1+q0)*(1+3*r0)]; return; % ------------------------------------------------------------------------------ function S=Sr(ni,q,r) %SR Array of solenoidal basis functions (quartic). qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors S=[3/64*ri*(1+q0)^2*(1-r^2)*((1+q0)*(8-9*q0+3*q^2)+(2-q0)*(1-5*r^2)), ... -1/64*(1+q0)^2*(1+r0)^2*(2-q0)*(7-26*r0+15*r^2), ... 3/64*ri/qi*(1+q0)^3*(1-r^2)*(1-q0)*(5-3*q0), ... 3/64*ri/qi^2*(1+q0)^3*(1-r^2)*(1-q0)^2, ... 1/16/qi*(1+q0)^2*(1+r0)*(1-q0)*(1-3*r0), ... 1/64/ri*(1+q0)^2*(1+r0)^2*(1-r0)*(1-5*r0)*(2-q0); ... -3/64*qi*(1+r0)^2*(1-q^2)*((1+r0)*(8-9*r0+3*r^2)+(2-r0)*(1-5*q^2)), ... 3/64*qi/ri*(1+r0)^3*(1-q^2)*(1-r0)*(5-3*r0), ... -1/64*(1+r0)^2*(1+q0)^2*(2-r0)*(7-26*q0+15*q^2), ... -1/64/qi*(1+r0)^2*(1+q0)^2*(1-q0)*(1-5*q0)*(2-r0), ... -1/16/ri*(1+r0)^2*(1+q0)*(1-r0)*(1-3*q0), ... -3/64*qi/ri^2*(1+r0)^3*(1-q^2)*(1-r0)^2]; return; function S=Sxr(ni,q,r) %SXR Array of x-derivatives of solenoidal basis functions. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors S=[9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), -3/64*qi*(1-q^2)*(1+r0)^2*(1-3*r0)*(7-5*r0), ... 3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), 3/64*ri/qi*(1+q0)^2*(1-r^2)*(1-q0)*(1-5*q0), ... 1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), 3/64*qi/ri*(1-q^2)*(1+r0)^2*(1-r0)*(1-5*r0); ... 3/32*qi^2*q0*(1+r0)^2*(r0*(10*q^2+3*r^2-7)+20-20*q^2-6*r^2), -3/32*qi^2/ri*q0*(1+r0)^3*(1-r0)*(5-3*r0),... 3/16*qi*(1-q^2)*(1+r0)^2*(2-r0)*(1+5*q0), 1/16*(1+q0)*(1+r0)^2*(2-r0)*(1+2*q0-5*q^2), ... 1/8*qi/ri*(1+r0)^2*(1-r0)*(1+3*q0), 3/32*qi^2/ri^2*q0*(1+r0)^3*(1-r0)^2]; return; function s=Syr(ni,q,r) %SYR Array of y-derivatives of solenoidal basis functions. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors s=[-3/32*ri^2*r0*(1+q0)^2*(q0*(10*r^2+3*q^2-7)+20-20*r^2-6*q^2), 3/16*ri*(1+q0)^2*(1-r^2)*(2-q0)*(1+5*r0),... -3/32*ri^2*r0/qi*(1+q0)^2*(1-q^2)*(5-3*q0), -3/32*ri^2/qi^2*r0*(1+q0)*(1-q^2)^2, ... -1/8*ri/qi*(1+q0)*(1-q^2)*(1+3*r0), -1/16*(1+q0)^2*(1+r0)*(2-q0)*(1+2*r0-5*r^2); ... -9/64*qi*ri*(1-q^2)*(1-r^2)*(6-5*q^2-5*r^2), 3/64*qi*(1+r0)^2*(1-q^2)*(1-3*r0)*(7-5*r0), ... -3/64*ri*(1+q0)^2*(1-r^2)*(1-3*q0)*(7-5*q0), -3/64*ri/qi*(1+q0)*(1-q^2)*(1-r^2)*(1-5*q0), ... -1/16*(1+q0)*(1+r0)*(1-3*q0)*(1-3*r0), -3/64*qi/ri*(1-q^2)*(1-r^2)*(1+r0)*(1-5*r0)]; return; function s=Sxxr(ni,q,r) %SXXR Array of second x-derivatives of solenoidal basis functions. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors s=[-9/32*qi^2*ri*q0*(1-r^2)*(11-10*q^2-5*r^2), 3/32*qi^2*q0*(1+r0)^2*(1-3*r0)*(7-5*r0) , ... -9/16*qi*ri*(1-r^2)*(1-q^2)*(1+5*q0), -3/16*ri*(1+q0)*(1-r^2)*(1+2*q0-5*q^2), ... -1/8*qi*(1+r0)*(1+3*q0)*(1-3*r0), -3/32*qi^2/ri*q0*(1+r0)^2*(1-r0)*(1-5*r0); ... 3/32*qi^3*(1+r0)^2*(6+(r0-2)*(30*q^2+3*r^2-7)), -3/32*qi^3/ri*(1+r0)^3*(1-r0)*(5-3*r0), ... 3/16*qi^2*(1+r0)^2*(2-r0)*(5-2*q0-15*q^2), 3/16*qi*(1+r0)^2*(2-r0)*(1-2*q0-5*q^2), ... 3/8*qi^2/ri*(1+r0)^2*(1-r0), 3/32*qi^3/ri^2*(1+r0)^3*(1-r0)^2]; return; function s=Syyr(ni,q,r) %SYYR Array of second y-derivatives of solenoidal basis functions. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors s=[-3/32*ri^3*(1+q0)^2*(6+(q0-2)*(30*r^2+3*q^2-7)), 3/16*ri^2*(1+q0)^2*(2-q0)*(5-2*r0-15*r^2), ... -3/32*ri^3/qi*(1+q0)^2*(1-q^2)*(5-3*q0), -3/32*ri^3/qi^2*(1+q0)*(1-q^2)^2, ... -3/8*ri^2/qi*(1+q0)*(1-q^2), -3/16*ri*(1+q0)^2*(2-q0)*(1-2*r0-5*r^2); ... 9/32*qi*ri^2*r0*(1-q^2)*(11-5*q^2-10*r^2), -9/16*qi*ri*(1-q^2)*(1-r^2)*(1+5*r0), ... 3/32*ri^2*r0*(1+q0)^2*(1-3*q0)*(7-5*q0), 3/32*ri^2/qi*(1+q0)*(1-q^2)*r0*(1-5*q0), ... 1/8*ri*(1+q0)*(1-3*q0)*(1+3*r0), 3/16*qi*(1-q^2)*(1+r0)*(1+2*r0-5*r^2)]; return; function s=Sxyr(ni,q,r) %SXYR Array of second (cross) xy-derivatives of solenoidal basis functions. qi=ni(1); q0=q*ni(1); q1=1+q0; ri=ni(2); r0=r*ni(2); r1=1+r0; % array of solenoidal vectors s=[9/32*qi*ri^2*r0*(1-q^2)*(5*q^2-11+10*r^2), 9/16*qi*ri*(1-q^2)*(1-r^2)*(5*r0+1), ... -3/32*ri^2*(q0+1)^2*r0*(-1+3*q0)*(-7+5*q0), 3/32*ri^2/qi*(1-q^2)*(1+q0)*r0*(-1+5*q0), ... 1/8*ri*(q0+1)*(-1+3*q0)*(1+3*r0), 3/16*qi*(1-q^2)*(1+r0)*(5*r^2-2*r0-1); ... -9/32*qi^2*ri*q0*(1-r^2)*(5*r^2-11+10*q^2), -3/32*qi^2*q0*(3*r0-1)*(5*r0-7)*(1+r0)^2, ... 9/16*qi*ri*(1-q^2)*(1-r^2)*(1+5*q0), -3/16*ri*(q0+1)*(1-r^2)*(5*q^2-2*q0-1), ... -1/8*qi*(1+r0)*(3*r0-1)*(1+3*q0), -3/32*qi^2/ri*q0*(5*r0-1)*(1-r^2)*(1+r0)]; return; function G=GBL(ni,q,r) % Transposed gradient (derivatives) of scalar bilinear mapping function. % The parameter ni can be a vector of coordinate pairs. G=[.25*ni(:,1).*(1+ni(:,2).*r), .25*ni(:,2).*(1+ni(:,1).*q)]; return; function D=BLxy(ni,q,r) % Transposed second cross-derivative of scalar bilinear mapping function. % The parameter ni can be a vector of coordinate pairs. D=[.25*ni(:,1).*ni(:,2)]; return;