% NORMAL2_3D Computes L2-norm in 3D
[err_l2]=normal2_2d(fdq,nq,errtype,x,wx,xx,jacx,y,wy,yy,jacy,...
z,wz,zz,jacz,xyz,ww,nov,un,u,uex);
Input : fdq = 0 uses Legendre Gauss quadrature formulas with nq nodes
in each element (exactness degree = 2*nq+1)
= 1 uses Legendre Gauss Lobatto quadrature formulas with npdx/y/z
quadrature nodes in each element.
Quadrature nodes are the nodes of the mesh.
nq = nodes (in each element and along each direction)
for GL quadrature formulas. Not used if fdq == 1
errtype = 0 for absolute error ||u-u_ex||
1 for relative error ||u-u_ex||/||u_ex||
x = column array with npdx LGL nodes in [-1,1]
wx= column array with npdx LGL weigths in [-1,1]
dx= derivative matrix produced with derlgl
xx(1:8,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie;...
x_V5_ie;x_V6_ie;x_V7_ie;x_V8_ie]
(ne=nex*ney is the global number of spectral elements)
jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
y = npdy LGL nodes in [-1,1], previously generated by xwlgl
wy = npdy LGL weigths in [-1,1], previously generated by xwlgl
dy= derivative matrix produced with derlgl
yy = 2-indexes array of size (8,ne):
yy(1:8,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie;...
y_V5_ie;y_V6_ie;y_V7_ie;y_V8_ie]
jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
z = npdy LGL nodes in [-1,1], previously generated by xwlgl
wz = npdy LGL weigths in [-1,1], previously generated by xwlgl
dz= derivative matrix produced with derlgl
zz = 2-indexes array of size (8,ne):
zz(1:8,ie)=[z_V1_ie;z_V2_ie;z_V3_ie;z_V4_ie;...
z_V5_ie;z_V6_ie;z_V7_ie;z_V8_ie]
jacz = array (length(jacz)=ne); jacz(ie)= (z_V5_ie-z_V1_ie)/2
xyz = column array with global mesh, length: noe=nov(npdx*npdy*npdz,ne)
ww = column array with global weigths, length: noe=nov(npdx*npdy*npdz,ne)
diag(ww) is the mass matrix associated to SEM discretization
nov = local to global map. 2-indexes array, size(nov)=[nov,ne]
un = column array with the numerical solution
u = column array with the evaluation of exact solution
uex = exact solution (uex=@(x,y,z)[uex(x,y,z)], with .*, .^, ./)
Output: err_l2 = error in L2-norm
Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Fundamentals in Single Domains"
Springer Verlag, Berlin Heidelberg New York, 2006.0001 function [err_l2]=normal2_3d(fdq,nq,errtype,x,wx,xx,jacx,y,wy,yy,jacy,... 0002 z,wz,dz,zz,jacz,xyz,ww,nov,un,u,uex); 0003 %% NORMAL2_3D Computes L2-norm in 3D 0004 % 0005 % [err_l2]=normal2_2d(fdq,nq,errtype,x,wx,xx,jacx,y,wy,yy,jacy,... 0006 % z,wz,zz,jacz,xyz,ww,nov,un,u,uex); 0007 % 0008 % Input : fdq = 0 uses Legendre Gauss quadrature formulas with nq nodes 0009 % in each element (exactness degree = 2*nq+1) 0010 % = 1 uses Legendre Gauss Lobatto quadrature formulas with npdx/y/z 0011 % quadrature nodes in each element. 0012 % Quadrature nodes are the nodes of the mesh. 0013 % nq = nodes (in each element and along each direction) 0014 % for GL quadrature formulas. Not used if fdq == 1 0015 % errtype = 0 for absolute error ||u-u_ex|| 0016 % 1 for relative error ||u-u_ex||/||u_ex|| 0017 % x = column array with npdx LGL nodes in [-1,1] 0018 % wx= column array with npdx LGL weigths in [-1,1] 0019 % dx= derivative matrix produced with derlgl 0020 % xx(1:8,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie;... 0021 % x_V5_ie;x_V6_ie;x_V7_ie;x_V8_ie] 0022 % (ne=nex*ney is the global number of spectral elements) 0023 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0024 % y = npdy LGL nodes in [-1,1], previously generated by xwlgl 0025 % wy = npdy LGL weigths in [-1,1], previously generated by xwlgl 0026 % dy= derivative matrix produced with derlgl 0027 % yy = 2-indexes array of size (8,ne): 0028 % yy(1:8,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie;... 0029 % y_V5_ie;y_V6_ie;y_V7_ie;y_V8_ie] 0030 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0031 % z = npdy LGL nodes in [-1,1], previously generated by xwlgl 0032 % wz = npdy LGL weigths in [-1,1], previously generated by xwlgl 0033 % dz= derivative matrix produced with derlgl 0034 % zz = 2-indexes array of size (8,ne): 0035 % zz(1:8,ie)=[z_V1_ie;z_V2_ie;z_V3_ie;z_V4_ie;... 0036 % z_V5_ie;z_V6_ie;z_V7_ie;z_V8_ie] 0037 % jacz = array (length(jacz)=ne); jacz(ie)= (z_V5_ie-z_V1_ie)/2 0038 % xyz = column array with global mesh, length: noe=nov(npdx*npdy*npdz,ne) 0039 % ww = column array with global weigths, length: noe=nov(npdx*npdy*npdz,ne) 0040 % diag(ww) is the mass matrix associated to SEM discretization 0041 % nov = local to global map. 2-indexes array, size(nov)=[nov,ne] 0042 % un = column array with the numerical solution 0043 % u = column array with the evaluation of exact solution 0044 % uex = exact solution (uex=@(x,y,z)[uex(x,y,z)], with .*, .^, ./) 0045 % 0046 % Output: err_l2 = error in L2-norm 0047 % 0048 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0049 % "Spectral Methods. Fundamentals in Single Domains" 0050 % Springer Verlag, Berlin Heidelberg New York, 2006. 0051 0052 % Written by Paola Gervasio 0053 % $Date: 2007/04/01$ 0054 0055 0056 0057 npdx=length(x); npdy=length(y);npdz=length(wz); 0058 mn=npdx*npdy; 0059 0060 [ldnov,ne]=size(nov); noe=nov(ldnov,ne); 0061 num=0; den=0; 0062 err_l2=0; 0063 0064 if fdq==0 0065 0066 % Legendre Gauss nodes and weigths 0067 0068 [xg,wxg] = xwlg(nq,-1,1); 0069 yg=xg;zg=xg; 0070 wyg=wxg;wzg=wxg; 0071 wxyg=zeros(ldnov,1); 0072 nq2=nq*nq; 0073 nq3=nq2*nq; 0074 for l=1:nq 0075 for j=1:nq 0076 for i=1:nq 0077 wxyg((l-1)*nq2+(j-1)*nq+i)=wxg(i)*wyg(j)*wzg(l); 0078 end 0079 end 0080 end 0081 0082 % Evaluation of Lagrange basis polynomials at quadrature nodes xg 0083 % 0084 [phix]= intlag_lgl (x,xg); 0085 [phiy]= intlag_lgl (y,yg); 0086 [phiz]= intlag_lgl (z,zg); 0087 0088 % Loop on spectral elements 0089 0090 for ie=1:ne 0091 un_loc=un(nov(1:ldnov,ie)); 0092 [norma_loc1,norma_loc2]=normal2_ex_loc(errtype,nq,uex,un_loc,wxyg,... 0093 jacx(ie),xx(1:8,ie),xg,phix,... 0094 jacy(ie),yy(1:8,ie),yg,phiy,... 0095 jacz(ie),zz(1:8,ie),zg,phiz); 0096 num=num+norma_loc1; 0097 den=den+norma_loc2; 0098 end 0099 0100 elseif fdq==1 0101 0102 num=sum((u-un).^2.*ww); 0103 den=sum(u.^2.*ww); 0104 0105 end 0106 0107 0108 if errtype==0 0109 err_l2=sqrt(num); 0110 elseif errtype==1 0111 if abs(den)>1.d-14; err_l2=sqrt(num/den); end 0112 end 0113 0114 return 0115 0116 function [norma_loc1,norma_loc2]=normal2_ex_loc(errtype,nq,uex,... 0117 un,wxyg,jacx,xx,xg,phix,... 0118 jacy,yy,yg,phiy,jacz,zz,zg,phiz); 0119 0120 % High degree Legendre Gaussian formulas to compute H1-norm error 0121 0122 [nq,npdx]=size(phix); 0123 [nq,npdy]=size(phiy); 0124 [nq,npdz]=size(phiz); 0125 0126 % mapping quadrature nodes on element ie 0127 0128 xxg=xg*jacx+(xx(2)+xx(1))*.5; 0129 yyg=yg*jacy+(yy(3)+yy(1))*.5; 0130 zzg=zg*jacz+(zz(5)+zz(1))*.5; 0131 0132 nq2=nq*nq;nq3=nq2*nq; 0133 for l=1:nq 0134 for j=1:nq 0135 for i=1:nq 0136 k=(l-1)*nq2+(j-1)*nq+i; 0137 xyg(k,1:3)=[xxg(i),yyg(j),zzg(l)]; 0138 end 0139 end 0140 end 0141 jac=jacx*jacy*jacz; 0142 0143 0144 % evaluation of exact solution at quadrature nodes. 0145 0146 [U]=uex(xyg(:,1),xyg(:,2),xyg(:,3)); 0147 0148 0149 % evaluate numerical solution at quadrature nodes. 0150 0151 un=reshape(un,npdx,npdy,npdz); 0152 for l=1:npdz 0153 un_z(l,:,:)=un(:,:,l); 0154 end 0155 0156 for j=1:npdy 0157 tt_z(:,:,j)=phiz*un_z(:,:,j); 0158 end 0159 for l=1:nq 0160 tt(:,:,l)=tt_z(l,:,:); 0161 end 0162 for l=1:nq 0163 un_i(:,:,l)=phix*(phiy*(tt(:,:,l))')'; 0164 end 0165 0166 un_i=un_i(:); 0167 0168 0169 % compute the sum 0170 0171 norma_loc1=sum((U-un_i).^2.*wxyg)*jac; 0172 0173 if errtype==0 0174 norma_loc2=0; 0175 else 0176 norma_loc2=sum(U.^2.*wxyg)*jac; 0177 end 0178 0179 return 0180 0181