NORMAH1_1D Computes H1-norm in 1D [err_h1]=normah1_1d(fdq, nq, errtype, u, uex, uexx,... x, wx, dx, xx, jacx, xy,nov) 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 on the nodes in which the discrete function is knonwn. nq = nodes (in each elemetn) 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|| u = numerical solution, column array of length noe uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./) uexx = first derivative of exact solution (uexx=@(x)[uex(xx)], with .*, .^, ./) x = column array with LGL nodes in [-1,1] wx= column array with LGL weigths in [-1,1] dx= derivative matrix produced with derlgl xx = 2-indexes array of size (2,ne): xx(1:2,ie)=[xa_ie;xb_ie] jacx = column array of length = ne, containing jacobians of the maps F_ie:[-1,1]---->[xa_ie,xb_ie] xy = mesh, column array of length noe nov = local to global map. 2-indexes array, size(nov)=[nov,ne] 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_1d(fdq, nq, errtype, u, uex, uexx,... 0002 x, wx, dx, xx, jacx, xy,nov) 0003 % NORMAH1_1D Computes H1-norm in 1D 0004 % 0005 % [err_h1]=normah1_1d(fdq, nq, errtype, u, uex, uexx,... 0006 % x, wx, dx, xx, jacx, xy,nov) 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 on the 0011 % nodes in which the discrete function is knonwn. 0012 % nq = nodes (in each elemetn) for GL quadrature formulas. Not used if 0013 % fdq=1 0014 % errtype = 0 for absolute error ||u-u_ex|| 0015 % 1 for relative error ||u-u_ex||/||u_ex|| 0016 % u = numerical solution, column array of length noe 0017 % uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./) 0018 % uexx = first derivative of exact solution 0019 % (uexx=@(x)[uex(xx)], with .*, .^, ./) 0020 % x = column array with LGL nodes in [-1,1] 0021 % wx= column array with LGL weigths in [-1,1] 0022 % dx= derivative matrix produced with derlgl 0023 % xx = 2-indexes array of size (2,ne): xx(1:2,ie)=[xa_ie;xb_ie] 0024 % jacx = column array of length = ne, containing 0025 % jacobians of the maps F_ie:[-1,1]---->[xa_ie,xb_ie] 0026 % xy = mesh, column array of length noe 0027 % nov = local to global map. 2-indexes array, size(nov)=[nov,ne] 0028 % 0029 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0030 % "Spectral Methods. Fundamentals in Single Domains" 0031 % Springer Verlag, Berlin Heidelberg New York, 2006. 0032 0033 % Written by Paola Gervasio 0034 % $Date: 2007/04/01$ 0035 0036 0037 npdx=length(x); [ldnov,ne]=size(nov); noe=nov(npdx,ne); 0038 num=0; den=0; 0039 err_h1=0; 0040 0041 if fdq==0 0042 0043 % Legendre Gauss nodes and weigths 0044 0045 [xg,wxg] = xwlg(nq,-1,1); 0046 0047 % Evaluation of Lagrange basis polynomials at quadrature nodes xg 0048 % 0049 [phix]= intlag_lgl (x,xg); 0050 0051 % Loop on spectral elements 0052 0053 for ie=1:ne 0054 u_loc=u(nov(1:npdx,ie)); 0055 [norma_loc1,norma_loc2]=normah1_ex_loc(errtype,npdx,uex,uexx,... 0056 u_loc,x,dx,jacx(ie),xx(1:2,ie),xg,wxg,nq,phix); 0057 num=num+norma_loc1; 0058 den=den+norma_loc2; 0059 end 0060 0061 elseif fdq==1 0062 0063 for ie=1:ne 0064 u_loc=u(nov(1:npdx,ie)); 0065 xy_loc=xy(nov(1:npdx,ie)); 0066 [norma_loc1,norma_loc2]=normah1_loc(errtype,npdx,uex,uexx,... 0067 u_loc,xy_loc, wx,dx,jacx(ie)); 0068 num=num+norma_loc1; 0069 den=den+norma_loc2; 0070 end 0071 0072 end 0073 if errtype==0 0074 err_h1=sqrt(num); 0075 elseif errtype==1 0076 if abs(den)>1.d-14; err_h1=sqrt(num/den); end 0077 end 0078 0079 return 0080 0081 function [norma_loc1,norma_loc2]=normah1_ex_loc(errtype,npdx,uex,uexx,... 0082 u, x,dx,jacx,xx,xg,wxg,nq,phix); 0083 0084 % High degree Legendre Gaussian formulas to compute H1-norm error 0085 0086 % mapping quadrature nodes on element ie 0087 0088 xxg=xg*jacx+(xx(2)+xx(1))*.5; 0089 0090 % evaluation of exact solution at quadrature nodes. 0091 0092 U=uex(xxg); 0093 UX=uexx(xxg); 0094 0095 % Compute numerical derivative 0096 ux=dx*u/jacx; 0097 0098 % evaluate numerical solution and its derivative at quadrature nodes. 0099 0100 u_i=phix*u; 0101 ux_i=phix*ux; 0102 % figure(1) 0103 % mesh(xxg,yyg,reshape(U,npx,npy)); 0104 % figure(2) 0105 % mesh(xxg,yyg,uloc_i) 0106 0107 % compute the sum 0108 0109 norma_loc1=sum(((U-u_i).^2+(UX-ux_i).^2).*wxg)*jacx; 0110 0111 if errtype==0 0112 norma_loc2=0; 0113 else 0114 norma_loc2=sum((U.^2+UX.^2).*wxg)*jacx; 0115 end 0116 0117 return 0118 0119 function [norma_loc1,norma_loc2]=normah1_loc(errtype,npdx,uex,uexx, u_loc,... 0120 xy_loc, wx,dx, jacx); 0121 0122 % LGL quadrature formulas on npdx nodes to compute H1-norm error 0123 0124 0125 % evaluation of exact solution at quadrature nodes. 0126 0127 U=uex(xy_loc); 0128 UX=uexx(xy_loc); 0129 0130 % Compute numerical derivative 0131 ux=dx*u_loc/jacx; 0132 0133 % compute the sum 0134 0135 norma_loc1=sum(((U-u_loc).^2+(UX-ux).^2).*wx)*jacx; 0136 if errtype==0 0137 norma_loc2=0; 0138 else 0139 norma_loc2=sum((U.^2+UX.^2).*wx)*jacx; 0140 end 0141 0142 return 0143