NORMAH1_3D Computes H1-norm in 3D
[err_h1]=normah1_3d(fdq,nq,errtype,x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,...
z,wz,dz,zz,jacz,xyz,ww,nov,un,u,uex,uex_x,uex_y,uex_z);
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 = 2-indexes array of size (8,ne)
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 .*, .^, ./)
uex_x = exact first x-derivative (uex_x=@(x,y,z)[uex_x(x,y,z)],
with .*, .^, ./)
uex_y = exact first y-derivative (uex_y=@(x,y,z)[uex_y(x,y,z)],
with .*, .^, ./)
uex_z = exact first z-derivative (uex_z=@(x,y,z)[uex_y(x,y,z)],
with .*, .^, ./)
Output: err_h1 = error in H1-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_h1]=normah1_3d(fdq,nq,errtype,x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0002 z,wz,dz,zz,jacz,xyz,ww,nov,un,u,uex,uex_x,uex_y,uex_z); 0003 % NORMAH1_3D Computes H1-norm in 3D 0004 % 0005 % [err_h1]=normah1_3d(fdq,nq,errtype,x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0006 % z,wz,dz,zz,jacz,xyz,ww,nov,un,u,uex,uex_x,uex_y,uex_z); 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 = 2-indexes array of size (8,ne) 0021 % xx(1:8,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie;... 0022 % x_V5_ie;x_V6_ie;x_V7_ie;x_V8_ie] 0023 % (ne=nex*ney is the global number of spectral elements) 0024 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0025 % y = npdy LGL nodes in [-1,1], previously generated by xwlgl 0026 % wy = npdy LGL weigths in [-1,1], previously generated by xwlgl 0027 % dy= derivative matrix produced with derlgl 0028 % yy = 2-indexes array of size (8,ne): 0029 % yy(1:8,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie;... 0030 % y_V5_ie;y_V6_ie;y_V7_ie;y_V8_ie] 0031 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0032 % z = npdy LGL nodes in [-1,1], previously generated by xwlgl 0033 % wz = npdy LGL weigths in [-1,1], previously generated by xwlgl 0034 % dz= derivative matrix produced with derlgl 0035 % zz = 2-indexes array of size (8,ne): 0036 % zz(1:8,ie)=[z_V1_ie;z_V2_ie;z_V3_ie;z_V4_ie;... 0037 % z_V5_ie;z_V6_ie;z_V7_ie;z_V8_ie] 0038 % jacz = array (length(jacz)=ne); jacz(ie)= (z_V5_ie-z_V1_ie)/2 0039 % xyz = column array with global mesh, length: noe=nov(npdx*npdy*npdz,ne) 0040 % ww = column array with global weigths, length: noe=nov(npdx*npdy*npdz,ne) 0041 % diag(ww) is the mass matrix associated to SEM discretization 0042 % nov = local to global map. 2-indexes array, size(nov)=[nov,ne] 0043 % un = column array with the numerical solution 0044 % u = column array with the evaluation of exact solution 0045 % uex = exact solution (uex=@(x,y,z)[uex(x,y,z)], with .*, .^, ./) 0046 % uex_x = exact first x-derivative (uex_x=@(x,y,z)[uex_x(x,y,z)], 0047 % with .*, .^, ./) 0048 % uex_y = exact first y-derivative (uex_y=@(x,y,z)[uex_y(x,y,z)], 0049 % with .*, .^, ./) 0050 % uex_z = exact first z-derivative (uex_z=@(x,y,z)[uex_y(x,y,z)], 0051 % with .*, .^, ./) 0052 % 0053 % Output: err_h1 = error in H1-norm 0054 % 0055 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0056 % "Spectral Methods. Fundamentals in Single Domains" 0057 % Springer Verlag, Berlin Heidelberg New York, 2006. 0058 0059 % Written by Paola Gervasio 0060 % $Date: 2007/04/01$ 0061 0062 0063 0064 npdx=length(x); npdy=length(y); npdz=length(wz); 0065 mn=npdx*npdy; 0066 [ldnov,ne]=size(nov); noe=nov(ldnov,ne); 0067 num=0; den=0; 0068 err_h1=0; 0069 0070 if fdq==0 0071 0072 % Legendre Gauss nodes and weigths 0073 0074 [xg,wxg] = xwlg(nq,-1,1); 0075 yg=xg;zg=xg; 0076 wyg=wxg;wzg=wxg; 0077 wxyg=zeros(ldnov,1); 0078 nq2=nq*nq;nq3=nq*nq2; 0079 for l=1:nq 0080 for j=1:nq 0081 for i=1:nq 0082 wxyg((l-1)*nq2+(j-1)*nq+i)=wxg(i)*wyg(j)*wzg(l); 0083 end 0084 end 0085 end 0086 0087 0088 % Evaluation of Lagrange basis polynomials at quadrature nodes xg 0089 % 0090 [phix]= intlag_lgl (x,xg); 0091 [phiy]= intlag_lgl (y,yg); 0092 [phiz]= intlag_lgl (z,zg); 0093 0094 % Loop on spectral elements 0095 0096 for ie=1:ne 0097 un_loc=un(nov(1:ldnov,ie)); 0098 [norma_loc1,norma_loc2]=normah1_ex_loc(errtype,nq,uex,uex_x,uex_y,uex_z,... 0099 un_loc,wxyg,... 0100 x,dx,jacx(ie),xx(1:8,ie),xg,phix,... 0101 y,dy,jacy(ie),yy(1:8,ie),yg,phiy,... 0102 z,dz,jacz(ie),zz(1:8,ie),zg,phiz); 0103 num=num+norma_loc1; 0104 den=den+norma_loc2; 0105 end 0106 0107 elseif fdq==1 0108 0109 for ie=1:ne 0110 nn=nov(1:ldnov,ie); 0111 u_loc=u(nn); 0112 un_loc=un(nn); 0113 xyz_loc=[xyz(nn,1),xyz(nn,2),xyz(nn,3)]; 0114 [norma_loc1,norma_loc2]=normah1_loc(errtype,uex_x,uex_y,uex_z,u_loc,un_loc,... 0115 xyz_loc, wx,dx,jacx(ie),wy,dy,jacy(ie),wz,dz,jacz(ie)); 0116 num=num+norma_loc1; 0117 den=den+norma_loc2; 0118 end 0119 0120 end 0121 if errtype==0 0122 err_h1=sqrt(num); 0123 elseif errtype==1 0124 if abs(den)>1.d-14; err_h1=sqrt(num/den); end 0125 end 0126 0127 return 0128 0129 function [norma_loc1,norma_loc2]=normah1_ex_loc(errtype,nq,uex,uex_x,uex_y,... 0130 uex_z,un,wxyg,x,dx,jacx,xx,xg,phix,... 0131 y,dy,jacy,yy,yg,phiy,z,dz,jacz,zz,zg,phiz); 0132 0133 % High degree Legendre Gaussian formulas to compute H1-norm error 0134 npdx=length(dx); 0135 npdy=length(dy); 0136 npdz=length(dz); 0137 mn=npdx*npdy; 0138 0139 % mapping quadrature nodes on element ie 0140 0141 xxg=xg*jacx+(xx(2)+xx(1))*.5; 0142 yyg=yg*jacy+(yy(3)+yy(1))*.5; 0143 zzg=zg*jacz+(zz(5)+zz(1))*.5; 0144 ldnov=mn*npdz; 0145 0146 nq2=nq*nq;nq3=nq2*nq; 0147 for l=1:nq 0148 for j=1:nq 0149 for i=1:nq 0150 k=(l-1)*nq2+(j-1)*nq+i; 0151 xyg(k,1:3)=[xxg(i),yyg(j),zzg(l)]; 0152 end 0153 end 0154 end 0155 jac=jacx*jacy*jacz; 0156 0157 % evaluation of exact solution at quadrature nodes. 0158 0159 [U]=zeros(nq3,1)+uex(xyg(:,1),xyg(:,2),xyg(:,3)); 0160 [UX]=zeros(nq3,1)+uex_x(xyg(:,1),xyg(:,2),xyg(:,3)); 0161 [UY]=zeros(nq3,1)+uex_y(xyg(:,1),xyg(:,2),xyg(:,3)); 0162 [UZ]=zeros(nq3,1)+uex_z(xyg(:,1),xyg(:,2),xyg(:,3)); 0163 0164 % Compute partial numerical derivative 0165 0166 uloc=reshape(un,npdx,npdy,npdz); 0167 for l=1:npdz 0168 ux(:,:,l)=dx*uloc(:,:,l)/jacx; 0169 uy(:,:,l)=(dy*(uloc(:,:,l))')'/jacy; 0170 uloc_z(l,:,:)=uloc(:,:,l); 0171 end 0172 0173 for j=1:npdy 0174 uz_z(:,:,j)=dz*uloc_z(:,:,j)/jacz; 0175 end 0176 0177 for l=1:npdz 0178 ux_z(l,:,:)=ux(:,:,l); 0179 uy_z(l,:,:)=uy(:,:,l); 0180 end 0181 0182 % evaluate numerical solution and its derivative at quadrature nodes. 0183 0184 for j=1:npdy 0185 tt_z(:,:,j)=phiz*uloc_z(:,:,j); 0186 end 0187 for l=1:nq 0188 tt(:,:,l)=tt_z(l,:,:); 0189 end 0190 for l=1:nq 0191 uloc_i(:,:,l)=phix*(phiy*(tt(:,:,l))')'; 0192 end 0193 0194 for j=1:npdy 0195 tt_z(:,:,j)=phiz*ux_z(:,:,j); 0196 end 0197 for l=1:nq 0198 tt(:,:,l)=tt_z(l,:,:); 0199 end 0200 for l=1:nq 0201 ux_i(:,:,l)=phix*(phiy*(tt(:,:,l))')'; 0202 end 0203 0204 for j=1:npdy 0205 tt_z(:,:,j)=phiz*uy_z(:,:,j); 0206 end 0207 for l=1:nq 0208 tt(:,:,l)=tt_z(l,:,:); 0209 end 0210 for l=1:nq 0211 uy_i(:,:,l)=phix*(phiy*(tt(:,:,l))')'; 0212 end 0213 0214 for j=1:npdy 0215 tt_z(:,:,j)=phiz*uz_z(:,:,j); 0216 end 0217 for l=1:nq 0218 tt(:,:,l)=tt_z(l,:,:); 0219 end 0220 for l=1:nq 0221 uz_i(:,:,l)=phix*(phiy*(tt(:,:,l))')'; 0222 end 0223 0224 0225 0226 0227 uloc_i=uloc_i(:);ux_i=ux_i(:); uy_i=uy_i(:); uz_i=uz_i(:); 0228 0229 0230 % compute the sum 0231 0232 norma_loc1=sum(((U-uloc_i).^2+(UX-ux_i).^2+(UY-uy_i).^2+(UZ-uz_i).^2).*wxyg)*jac; 0233 0234 if errtype==0 0235 norma_loc2=0; 0236 else 0237 norma_loc2=sum((U.^2+UX.^2+UY.^2+UZ.^2).*wxyg)*jac; 0238 end 0239 0240 return 0241 0242 function [norma_loc1,norma_loc2]=normah1_loc(errtype,uex_x,uex_y,uex_z,U,... 0243 un_loc,xyz_loc, wx,dx,jacx,wy,dy,jacy,wz,dz,jacz); 0244 0245 % LGL quadrature formulas on npdx nodes to compute H1-norm error 0246 npdx=length(dx); 0247 npdy=length(dy); 0248 npdz=length(dz); 0249 mn=npdx*npdy; 0250 ldnov=mn*npdz; 0251 ww=zeros(ldnov,1); 0252 for l=1:npdz 0253 for j=1:npdy 0254 for i=1:npdx 0255 ww((l-1)*mn+(j-1)*npdx+i)=wx(i)*wy(j)*wz(l); 0256 end 0257 end 0258 end 0259 jac=jacx*jacy*jacz; 0260 0261 0262 % evaluation of exact solution at quadrature nodes. 0263 0264 UX=uex_x(xyz_loc(:,1),xyz_loc(:,2),xyz_loc(:,3)); 0265 UY=uex_y(xyz_loc(:,1),xyz_loc(:,2),xyz_loc(:,3)); 0266 UZ=uex_z(xyz_loc(:,1),xyz_loc(:,2),xyz_loc(:,3)); 0267 0268 % Compute numerical derivative 0269 0270 un_loc=reshape(un_loc,npdx,npdy,npdz); 0271 for l=1:npdz 0272 ux(:,:,l)=dx*un_loc(:,:,l)/jacx; 0273 uy(:,:,l)=(dy*(un_loc(:,:,l))')'/jacy; 0274 un_loc_z(l,:,:)=un_loc(:,:,l); 0275 end 0276 for j=1:npdy 0277 uz_z(:,:,j)=dz*un_loc_z(:,:,j)/jacz; 0278 end 0279 0280 for l=1:npdz 0281 uz(:,:,l)=uz_z(l,:,:); 0282 end 0283 ux=ux(:); 0284 uy=uy(:); 0285 uz=uz(:); 0286 un_loc=un_loc(:); 0287 0288 0289 % compute the sum 0290 0291 norma_loc1=sum(((U-un_loc).^2+(UX-ux).^2+(UY-uy).^2+(UZ-uz).^2).*ww)*jac; 0292 if errtype==0 0293 norma_loc2=0; 0294 else 0295 norma_loc2=sum((U.^2+UX.^2+UY.^2+UZ.^2).*ww)*jac; 0296 end 0297 0298 return 0299