NORMAH1_2D Computes H1-norm in 2D
[err_h1]=normah1_2d(fdq,nq,errtype,x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,...
xy,ww,nov,un,u,uex,uex_x,uex_y);
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
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 (4,ne)
xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_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 (4,ne):
yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie]
jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
xy = column array with global mesh, length: noe=nov(npdx*npdy,ne)
ww = column array with global weigths, length: noe=nov(npdx*npdy,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)[uex(x,y)], with .*, .^, ./)
uex_x = exact first x-derivative (uex_x=@(x,y)[uex_x(x,y)],
with .*, .^, ./)
uex_y = exact first y-derivative (uex_y=@(x,y)[uex_y(x,y)],
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_2d(fdq,nq,errtype,x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0002 xy,ww,nov,un,u,uex,uex_x,uex_y); 0003 % NORMAH1_2D Computes H1-norm in 2D 0004 % 0005 % [err_h1]=normah1_2d(fdq,nq,errtype,x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0006 % xy,ww,nov,un,u,uex,uex_x,uex_y); 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 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 (4,ne) 0021 % xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_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 (4,ne): 0028 % yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie] 0029 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0030 % xy = column array with global mesh, length: noe=nov(npdx*npdy,ne) 0031 % ww = column array with global weigths, length: noe=nov(npdx*npdy,ne) 0032 % diag(ww) is the mass matrix associated to SEM discretization 0033 % nov = local to global map. 2-indexes array, size(nov)=[nov,ne] 0034 % un = column array with the numerical solution 0035 % u = column array with the evaluation of exact solution 0036 % uex = exact solution (uex=@(x,y)[uex(x,y)], with .*, .^, ./) 0037 % uex_x = exact first x-derivative (uex_x=@(x,y)[uex_x(x,y)], 0038 % with .*, .^, ./) 0039 % uex_y = exact first y-derivative (uex_y=@(x,y)[uex_y(x,y)], 0040 % with .*, .^, ./) 0041 % 0042 % Output: err_h1 = error in H1-norm 0043 % 0044 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0045 % "Spectral Methods. Fundamentals in Single Domains" 0046 % Springer Verlag, Berlin Heidelberg New York, 2006. 0047 0048 % Written by Paola Gervasio 0049 % $Date: 2007/04/01$ 0050 0051 0052 0053 npdx=length(x); npdy=length(y); 0054 [ldnov,ne]=size(nov); noe=nov(ldnov,ne); 0055 num=0; den=0; 0056 err_h1=0; 0057 0058 if fdq==0 0059 0060 % Legendre Gauss nodes and weigths 0061 0062 [xg,wxg] = xwlg(nq,-1,1); 0063 [yg,wyg] = xwlg(nq,-1,1); 0064 [wxg1,wyg1]=meshgrid(wxg,wyg); wxyg=wxg1.*wyg1;wxyg=wxyg'; wxyg=wxyg(:); 0065 clear wxg1 wyg1; 0066 0067 % Evaluation of Lagrange basis polynomials at quadrature nodes xg 0068 % 0069 [phix]= intlag_lgl (x,xg); 0070 [phiy]= intlag_lgl (y,yg); 0071 0072 % Loop on spectral elements 0073 0074 for ie=1:ne 0075 un_loc=un(nov(1:ldnov,ie)); 0076 [norma_loc1,norma_loc2]=normah1_ex_loc(errtype,nq,uex,uex_x,uex_y,un_loc,... 0077 wxyg, x,dx,jacx(ie),xx(1:4,ie),xg,phix,... 0078 y,dy,jacy(ie),yy(1:4,ie),yg,phiy); 0079 num=num+norma_loc1; 0080 den=den+norma_loc2; 0081 end 0082 0083 elseif fdq==1 0084 0085 for ie=1:ne 0086 u_loc=u(nov(1:ldnov,ie)); 0087 un_loc=un(nov(1:ldnov,ie)); 0088 xy_loc=[xy(nov(1:ldnov,ie),1),xy(nov(1:ldnov,ie),2)]; 0089 [norma_loc1,norma_loc2]=normah1_loc(errtype,uex_x,uex_y,u_loc,un_loc,... 0090 xy_loc, wx,dx,jacx(ie),wy,dy,jacy(ie)); 0091 num=num+norma_loc1; 0092 den=den+norma_loc2; 0093 end 0094 0095 end 0096 if errtype==0 0097 err_h1=sqrt(num); 0098 elseif errtype==1 0099 if abs(den)>1.d-14; err_h1=sqrt(num/den); end 0100 end 0101 0102 return 0103 0104 function [norma_loc1,norma_loc2]=normah1_ex_loc(errtype,nq,uex,uex_x,uex_y,... 0105 un,wxyg,x,dx,jacx,xx,xg,phix,... 0106 y,dy,jacy,yy,yg,phiy); 0107 0108 % High degree Legendre Gaussian formulas to compute H1-norm error 0109 npdx=length(dx); 0110 npdy=length(dy); 0111 0112 0113 % mapping quadrature nodes on element ie 0114 0115 xxg=xg*jacx+(xx(2)+xx(1))*.5; 0116 yyg=yg*jacy+(yy(3)+yy(1))*.5; 0117 [xg1,yg1]=meshgrid(xxg,yyg); xg1=xg1'; yg1=yg1'; xyg=[xg1(:),yg1(:)]; 0118 clear xg1; clear yg1; 0119 0120 0121 % evaluation of exact solution at quadrature nodes. 0122 0123 [U]=uex(xyg(:,1),xyg(:,2)); 0124 [UX]=uex_x(xyg(:,1),xyg(:,2)); 0125 [UY]=uex_y(xyg(:,1),xyg(:,2)); 0126 0127 % Compute partial numerical derivative 0128 uloc=reshape(un,npdx,npdy); 0129 ux=dx*uloc/jacx; 0130 uy=(dy*uloc')'/jacy; 0131 0132 0133 % evaluate numerical solution and its derivative at quadrature nodes. 0134 0135 uloc_i=phix*(phiy*uloc')'; 0136 ux_i=phix*(phiy*ux')'; 0137 uy_i=phix*(phiy*uy')'; 0138 0139 uloc_i=uloc_i(:);ux_i=ux_i(:); uy_i=uy_i(:); 0140 0141 0142 % compute the sum 0143 0144 norma_loc1=sum(((U-uloc_i).^2+(UX-ux_i).^2+(UY-uy_i).^2).*wxyg)*jacx*jacy; 0145 0146 if errtype==0 0147 norma_loc2=0; 0148 else 0149 norma_loc2=sum((U.^2+UX.^2+UY.^2).*wxyg)*jacx*jacy; 0150 end 0151 0152 return 0153 0154 function [norma_loc1,norma_loc2]=normah1_loc(errtype,uex_x,uex_y,U,un_loc,... 0155 xy_loc, wx,dx,jacx,wy,dy,jacy); 0156 0157 % LGL quadrature formulas on npdx nodes to compute H1-norm error 0158 npdx=length(dx); 0159 npdy=length(dy); 0160 [wx1,wy1]=meshgrid(wx,wy); ww=wx1.*wy1;ww=ww(:); clear wx1 wy1; 0161 0162 0163 0164 % evaluation of exact solution at quadrature nodes. 0165 0166 UX=uex_x(xy_loc(:,1),xy_loc(:,2)); 0167 UY=uex_y(xy_loc(:,1),xy_loc(:,2)); 0168 0169 % Compute numerical derivative 0170 0171 un_loc=reshape(un_loc,npdx,npdy); 0172 ux=dx*un_loc/jacx; ux=ux(:); 0173 uy=(dy*un_loc')'/jacy; uy=uy(:); 0174 un_loc=un_loc(:); 0175 0176 0177 % compute the sum 0178 0179 norma_loc1=sum(((U-un_loc).^2+(UX-ux).^2+(UY-uy).^2).*ww)*jacx*jacy; 0180 if errtype==0 0181 norma_loc2=0; 0182 else 0183 norma_loc2=sum((U.^2+UX.^2+UY.^2).*ww)*jacx*jacy; 0184 end 0185 0186 return 0187