LAP_2D Numerical solution of the 2D b.v.p. -Delta u + gam u = f
-Delta u + gam u = f x in Omega
+ Dirichlet bc
by SEM Numerical Integration with LGL quadrature formulas.
[xy,un,D,param]=lap_2d(xa,xb,ya,yb,gam,...
uex,uex_x,uex_y,ff,g,h,nex,nx,ney,ny,gammax,gammay,param);
Omega=(xa,xb) x (ya,yb)
side 3
V4 __________ V3
| |
| |
side 4 | Omega | side 2
| |
|__________|
V1 V2
side 1
__________________________
| | | | |
| 3 | 6 | 9 | 12 | Omega and spectral elements
| | | | | ordering
__________________________
| | | | |
| 2 | 5 | 8 | 11 |
| | | | |
__________________________
| | | | |
| 1 | 4 | 7 | 10 |
| | | | |
__________________________
Input: xa= abscissa of either vertex V1 or vertex V4
xb= abscissa of either vertex V2 or vertex V3
ya= ordinate of either vertex V1 or vertex V2
yb= ordinate of either vertex V3 or vertex V4
gam = coefficient of zeroth order term (constant>=0)
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 .*, .^, ./)
ff = r.h.s. solution (ff=@(x,y)[ff(x,y)], with .*, .^, ./)
g = @(x,y)[....] function handle to the expression of Dirichlet
boundary data
h =@(x,y)[....] function handle to the expression of Neumann
boundary data. It is a vector of 4 functions, each for any
side.
cb= a string of four elements, each for any side of \partial\Omega.
cb(i) refers to side number "i" of \partial\Omega
'd' character stands for Dirichlet boundary conditions on that side
'n' character stands for Neumann boundary conditions on that side
nex = number of elements along x-direction
nx = polynomial degree in each element (the same in each element)
along x-direction
ney = number of elements along y-direction
ny = polynomial degree in each element (the same in each element)
along y-direction
gammax = column or row array of length nex-1.
If the deomposition of Omega is not
uniform, gammax is the vector of position of interfaces between
spectral elements along x-direction. If gammax=[], uniform
decomposition is used.
gammay = column or row array of length ney-1.
If the deomposition of Omega is not
uniform, gammay is the vector of position of interfaces between
spectral elements along y-direction. If gammay=[], uniform
decomposition is used.
param = array of parameters
param(1) = 1 : SEM-NI approach, see Sect. 5.3.3, CHQZ3
= 2 : Patching approach: See Sect. 5.13, CHQZ3
param(2) = 0 : no reordering of the matrix, see pag. 194, CHQZ2
1 : Cuthill-McKee ordering (by symrcm function of
Matlab)
2 : Symmetric minimum degree ordering (by symamd function of
Matlab)
ONLY used if param(3)==1
param(3) = 1 : Solve the differential problem with a direct method
(Cholesky) and compute errors on the exact solution
2 : Compute extrema eigenvalues of stiffness matrix
3 : Compute interface Schur complement, solve the diff.
problem by Schur complement
4: Compute Schur complement and its extrema eigenvalues
5: Compute and plot all the eigenvalues of A (stiffness
matrix)
6: Compute extrema eigenvalues of M^{-1}A, i.e.
the generalized eigenvalues of Av=lambda* M v,
where M is the SEM-NI mass matrix
param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm)
on the exact solution
2: no errors are computed
param(5) = 0: LG quadrature formulas with high precision degree are
used to compute norms (exact norms)
1: LGL quadrature formulas with npdx,npdy nodes are
used to compute norms (discrete norms)
(used only if param(4) == 1)
param(6) = number of nodes for high degree quadrature formula,
(used only if param(5) == 0 & param(4) == 1)
param(7) = 0: absolute errors are computed
1: relative errors are computed
(used only if param(4) == 1)
param(8) = 0: not plot the solution
1: plot the solution (mesh)
2: plot the solution (surf)
3: plot the solution (contour)
param(9) = number of nodes in each element and along each direction
for plotting the solution
(used only if param(8) > 0)
Output: xy = 2-indexes array wiht coordinates of 2D LGL mesh
un = numerical solution
D = (if param(3)==5) eigenvalues of A
param(21) = (if param(3)==1) cputime for reordering stiffness
matrix befor solving
param(22) = (if param(3)==1) cputime for solving linear system
by Choleski factorization
param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of
stiffness matrix A
param(25,26) = (if param(3)==4) extrema eigenvalues (max,min) of
interface Schur complement
param(27,28) = (if param(3)==6) extrema generalized eigenvalues (max,min) of
A v = lambda M v, where A is the stiffness matrix and M is the
mass matrix
param(29,30,31) = (if param(4)==1 & (param(3)==1 | param(3)==3))
errors on the exact solution:
L^inf-norm, L2-norm, 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 [xy,un,D,param]=lap_2d(xa,xb,ya,yb,gam,... 0002 uex,uex_x,uex_y,ff,g,h,cb,nex,nx,ney,ny,gammax,gammay,param); 0003 % LAP_2D Numerical solution of the 2D b.v.p. -Delta u + gam u = f 0004 % 0005 % -Delta u + gam u = f x in Omega 0006 % 0007 % + Dirichlet bc 0008 % 0009 % by SEM Numerical Integration with LGL quadrature formulas. 0010 % 0011 % [xy,un,D,param]=lap_2d(xa,xb,ya,yb,gam,... 0012 % uex,uex_x,uex_y,ff,g,h,nex,nx,ney,ny,gammax,gammay,param); 0013 % 0014 % Omega=(xa,xb) x (ya,yb) 0015 % 0016 % side 3 0017 % V4 __________ V3 0018 % | | 0019 % | | 0020 % side 4 | Omega | side 2 0021 % | | 0022 % |__________| 0023 % V1 V2 0024 % side 1 0025 % 0026 % 0027 % __________________________ 0028 % | | | | | 0029 % | 3 | 6 | 9 | 12 | Omega and spectral elements 0030 % | | | | | ordering 0031 % __________________________ 0032 % | | | | | 0033 % | 2 | 5 | 8 | 11 | 0034 % | | | | | 0035 % __________________________ 0036 % | | | | | 0037 % | 1 | 4 | 7 | 10 | 0038 % | | | | | 0039 % __________________________ 0040 % 0041 % 0042 % Input: xa= abscissa of either vertex V1 or vertex V4 0043 % xb= abscissa of either vertex V2 or vertex V3 0044 % ya= ordinate of either vertex V1 or vertex V2 0045 % yb= ordinate of either vertex V3 or vertex V4 0046 % gam = coefficient of zeroth order term (constant>=0) 0047 % uex = exact solution (uex=@(x,y)[uex(x,y)], with .*, .^, ./) 0048 % uex_x = exact first x-derivative (uex_x=@(x,y)[uex_x(x,y)], with .*, .^, ./) 0049 % uex_y = exact first y-derivative (uex_y=@(x,y)[uex_y(x,y)], with .*, .^, ./) 0050 % ff = r.h.s. solution (ff=@(x,y)[ff(x,y)], with .*, .^, ./) 0051 % g = @(x,y)[....] function handle to the expression of Dirichlet 0052 % boundary data 0053 % h =@(x,y)[....] function handle to the expression of Neumann 0054 % boundary data. It is a vector of 4 functions, each for any 0055 % side. 0056 % cb= a string of four elements, each for any side of \partial\Omega. 0057 % cb(i) refers to side number "i" of \partial\Omega 0058 % 'd' character stands for Dirichlet boundary conditions on that side 0059 % 'n' character stands for Neumann boundary conditions on that side 0060 % nex = number of elements along x-direction 0061 % nx = polynomial degree in each element (the same in each element) 0062 % along x-direction 0063 % ney = number of elements along y-direction 0064 % ny = polynomial degree in each element (the same in each element) 0065 % along y-direction 0066 % gammax = column or row array of length nex-1. 0067 % If the deomposition of Omega is not 0068 % uniform, gammax is the vector of position of interfaces between 0069 % spectral elements along x-direction. If gammax=[], uniform 0070 % decomposition is used. 0071 % gammay = column or row array of length ney-1. 0072 % If the deomposition of Omega is not 0073 % uniform, gammay is the vector of position of interfaces between 0074 % spectral elements along y-direction. If gammay=[], uniform 0075 % decomposition is used. 0076 % param = array of parameters 0077 % param(1) = 1 : SEM-NI approach, see Sect. 5.3.3, CHQZ3 0078 % = 2 : Patching approach: See Sect. 5.13, CHQZ3 0079 % param(2) = 0 : no reordering of the matrix, see pag. 194, CHQZ2 0080 % 1 : Cuthill-McKee ordering (by symrcm function of 0081 % Matlab) 0082 % 2 : Symmetric minimum degree ordering (by symamd function of 0083 % Matlab) 0084 % ONLY used if param(3)==1 0085 % param(3) = 1 : Solve the differential problem with a direct method 0086 % (Cholesky) and compute errors on the exact solution 0087 % 2 : Compute extrema eigenvalues of stiffness matrix 0088 % 3 : Compute interface Schur complement, solve the diff. 0089 % problem by Schur complement 0090 % 4: Compute Schur complement and its extrema eigenvalues 0091 % 5: Compute and plot all the eigenvalues of A (stiffness 0092 % matrix) 0093 % 6: Compute extrema eigenvalues of M^{-1}A, i.e. 0094 % the generalized eigenvalues of Av=lambda* M v, 0095 % where M is the SEM-NI mass matrix 0096 % param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm) 0097 % on the exact solution 0098 % 2: no errors are computed 0099 % param(5) = 0: LG quadrature formulas with high precision degree are 0100 % used to compute norms (exact norms) 0101 % 1: LGL quadrature formulas with npdx,npdy nodes are 0102 % used to compute norms (discrete norms) 0103 % (used only if param(4) == 1) 0104 % param(6) = number of nodes for high degree quadrature formula, 0105 % (used only if param(5) == 0 & param(4) == 1) 0106 % param(7) = 0: absolute errors are computed 0107 % 1: relative errors are computed 0108 % (used only if param(4) == 1) 0109 % param(8) = 0: not plot the solution 0110 % 1: plot the solution (mesh) 0111 % 2: plot the solution (surf) 0112 % 3: plot the solution (contour) 0113 % param(9) = number of nodes in each element and along each direction 0114 % for plotting the solution 0115 % (used only if param(8) > 0) 0116 % 0117 % Output: xy = 2-indexes array wiht coordinates of 2D LGL mesh 0118 % un = numerical solution 0119 % D = (if param(3)==5) eigenvalues of A 0120 % param(21) = (if param(3)==1) cputime for reordering stiffness 0121 % matrix befor solving 0122 % param(22) = (if param(3)==1) cputime for solving linear system 0123 % by Choleski factorization 0124 % param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of 0125 % stiffness matrix A 0126 % param(25,26) = (if param(3)==4) extrema eigenvalues (max,min) of 0127 % interface Schur complement 0128 % param(27,28) = (if param(3)==6) extrema generalized eigenvalues (max,min) of 0129 % A v = lambda M v, where A is the stiffness matrix and M is the 0130 % mass matrix 0131 % param(29,30,31) = (if param(4)==1 & (param(3)==1 | param(3)==3)) 0132 % errors on the exact solution: 0133 % L^inf-norm, L2-norm, H1-norm 0134 % 0135 % 0136 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0137 % "Spectral Methods. Fundamentals in Single Domains" 0138 % Springer Verlag, Berlin Heidelberg New York, 2006. 0139 0140 % Written by Paola Gervasio 0141 % $Date: 2007/04/01$ 0142 0143 0144 0145 xy=[]; un=[]; D=[]; 0146 0147 npdx=nx+1; npdy=ny+1; ldnov=npdx*npdy; mn=ldnov; ne=nex*ney; 0148 0149 [x,wx]=xwlgl(npdx); [dx]=derlgl(x,npdx); 0150 [y,wy]=xwlgl(npdy); [dy]=derlgl(y,npdy); 0151 % 0152 % nov construction 0153 % 0154 [nov]=cosnov_2d(npdx,nex,npdy,ney); 0155 noe=nov(ldnov,ne); 0156 0157 % Mesh generation 0158 0159 [xx,yy,jacx,jacy,xy,ww,ifro]=mesh_2d(xa,xb,ya,yb,cb,nex,ney,npdx,npdy,... 0160 nov,x,wx,y,wy,gammax,gammay); 0161 0162 % Matrix assembling 0163 0164 A=stiff_2d_se(npdx,nex,npdy,ney,nov,wx,dx,jacx,wy,dy,jacy); 0165 M=spdiags(ww,0,noe,noe); 0166 0167 if gam ~=0 0168 A=A+gam*M; 0169 end 0170 0171 % Rigth Hand Side 0172 f=ff(xy(:,1),xy(:,2)).*ww; 0173 0174 % Neumann boundary condition 0175 [f]=neumannbc(f,h,jacx,jacy,wx,wy,xy,ifro,nov); 0176 0177 % 0178 % Patching approach 0179 % 0180 if (param(1)==2) 0181 [A,f]=patch_se(A,f,ifro,nov,dx,jacx,dy,jacy); 0182 end 0183 0184 % Lists of internal, boundary, interface nodes. All lists are referred to 0185 % global ordering on Omega 0186 % 0187 % lbor: list of Dirichlet boundary nodes 0188 % lint: list of internal nodes 0189 % lintint: list of internal nodes which are internal to spectral elements 0190 % lgamma: list of internal nodes which are on the interface between spectral elements 0191 0192 0193 0194 [ldir,lint,lintint,lgamma,ifro]=liste(ifro,nov); 0195 0196 0197 % Setting Dirichlet boundary conditions on both matrix and r.h.s 0198 ub=zeros(noe,1); 0199 for i=1:noe 0200 if(ifro(i)==1) 0201 ub(i)=uex(xy(i,1),xy(i,2)); 0202 end 0203 end 0204 0205 if length(ldir)>0 0206 f(lint)=f(lint)-A(lint,ldir)*ub(ldir); 0207 A=A(lint,lint); M=M(lint,lint); 0208 f=f(lint); 0209 noei=length(lint); 0210 else 0211 noei=noe; 0212 end 0213 0214 % selection on param(3) 0215 0216 0217 if (param(3)==1) 0218 0219 if (param(1)==1) 0220 % SEM-NI discretization A is symmetric 0221 % Solve the linear system by Cholesky factorization 0222 0223 if param(2)==0 0224 % no reordering 0225 t_ord=0; 0226 t=cputime; L=chol(A); un=L\(L'\f);t_solv=cputime-t; 0227 0228 elseif param(2)==1 0229 % Cuthill-McKee ordering 0230 t=cputime; pp=symrcm(A);t_ord=cputime-t; 0231 t=cputime; L=chol(A(pp,pp)); un(pp)=L\(L'\f(pp));t_solv=cputime-t; 0232 0233 elseif param(2)==2 0234 % Symmetric Minimum degree ordering 0235 t=cputime; pp=symamd(A);t_ord=cputime-t; 0236 t=cputime; L=chol(A(pp,pp)); un(pp)=L\(L'\f(pp));t_solv=cputime-t; 0237 end 0238 elseif (param(1)==2) 0239 % Patching discretization A is not symmetric 0240 % Solve the linear system by LU factorization 0241 0242 if param(2)==0 0243 % no reordering 0244 t_ord=0; 0245 t=cputime; un=A\f;t_solv=cputime-t; 0246 0247 elseif param(2)==1 0248 % Cuthill-McKee ordering 0249 t=cputime; pp=symrcm(A);t_ord=cputime-t; 0250 t=cputime; [L,U,PP]=lu(A(pp,pp)); un(pp)=U\(L\(PP*f(pp)));t_solv=cputime-t; 0251 0252 elseif param(2)==2 0253 % Symmetric Minimum degree ordering 0254 t=cputime; pp=symamd(A);t_ord=cputime-t; 0255 t=cputime; [L,U,PP]=lu(A(pp,pp)); un(pp)=U\(L\(PP*f(pp)));t_solv=cputime-t; 0256 end 0257 end 0258 0259 param(21)=t_ord; 0260 param(22)=t_solv; 0261 0262 % 0263 ub(lint)=un; 0264 un=ub; 0265 0266 0267 0268 elseif(param(3)==2) 0269 0270 % Compute extrema eigenvalues of the stiffness matrix 0271 keigs=1; 0272 if(param(1)==1) 0273 % Symmetric problem 0274 opts.issym=1;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0275 lambda_max=eigs(A,keigs,'lm',opts); 0276 lambda_min=eigs(A,1,'sm',opts); 0277 else 0278 opts.issym=0;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0279 lambda_max=eigs(A,keigs,'lm',opts); 0280 lambda_min=eigs(A,1,'sm',opts); 0281 end 0282 param(23)=lambda_max; 0283 param(24)=lambda_min; 0284 0285 0286 elseif(param(3)==3) 0287 0288 % Compute interface Schur complement, solve the differential problem by 0289 % the Schur complement 0290 0291 Sigma=A(lgamma,lgamma)-A(lgamma,lintint)*inv(A(lintint,lintint))*A(lintint,lgamma); 0292 0293 % Solve Schur linear system and internal linear systems by \ of matlab 0294 0295 fsigma=f(lgamma)-A(lgamma,lintint)*inv(A(lintint,lintint))*f(lintint); 0296 xsigma=Sigma\fsigma; 0297 xint=A(lintint,lintint)\(f(lintint)-A(lintint,lgamma)*xsigma); 0298 u2=zeros(noei,1); 0299 u2(lgamma)=xsigma; 0300 u2(lintint)=xint; 0301 ub(lint)=u2; 0302 un=ub; 0303 0304 0305 elseif(param(3)==4) 0306 0307 % Compute interface Schur complement and its extrema eigenvalues 0308 0309 Sigma=A(lgamma,lgamma)-A(lgamma,lintint)*inv(A(lintint,lintint))*A(lintint,lgamma); 0310 0311 clear A; clear f; 0312 opts.issym=1;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0313 lambda_max=eigs(Sigma,1,'lm',opts); 0314 lambda_min=eigs(Sigma,1,'sm',opts); 0315 param(25)=lambda_max; 0316 param(26)=lambda_min; 0317 0318 0319 elseif(param(3)==5) 0320 % compute and plot eigenvalues of A 0321 D=eig(full(A)); 0322 fig=figure(..., 0323 'Name','Spectrum of SEM-NI stiffness matrix',... 0324 'Visible','on') 0325 plot(real(D),imag(D),'o') 0326 set(gca,'Fontname','Times','Fontsize',18) 0327 xlabel('Re'); ylabel('Im') 0328 if param(1)==1 0329 set(gca,'Fontname','Times','Fontsize',18) 0330 title(['N=',num2str(nx),' H=',num2str(2/nex), ' SEM-NI']) 0331 else 0332 set(gca,'Fontname','Times','Fontsize',18) 0333 title(['N=',num2str(nx),' H=',num2str(2/nex), ' Patching']) 0334 end 0335 0336 0337 elseif(param(3)==6) 0338 keigs=1; 0339 % compute extrema eigenvalues of M^{-1}A, i.e. 0340 % the generalized eigenvalues of Av=lambda* M v 0341 if(param(1)==1) 0342 opts.issym=1;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0343 lambda_max=eigs(A,M,keigs,'lm',opts); 0344 lambda_min=eigs(A,M,keigs,'sm',opts); 0345 param(27)=lambda_max; 0346 param(28)=lambda_min; 0347 0348 elseif(param(2)==2) 0349 disp('No eigenvalues are computed') 0350 end 0351 0352 end 0353 0354 if (param(4)==1 & (param(3)==1 | param(3)==3)) 0355 % compute errors on the exact solution 0356 [err_inf,err_h1,err_l2]=errors_2d(x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0357 xy,ww,nov,un,uex,uex_x,uex_y,param); 0358 param(29:31)=[err_inf;err_h1;err_l2]; 0359 end 0360 0361 if (param(8)>0 & (param(3)==1 | param(3)==3)) 0362 if param(8)==1 0363 command='mesh'; 0364 elseif param(8)==2 0365 command='surf'; 0366 elseif param(8)==3; 0367 command='contour' 0368 end 0369 % plot the solution 0370 fig=figure(..., 0371 'Name','SEM-NI solution of -Delta u + gam u=f, Dirichlet bc',... 0372 'Visible','on'); 0373 [ha]=plot_sem_2d(fig,command,nex,ney,x,xx,jacx,y,yy,jacy,xy,ww,nov,un,param(9)); 0374 0375 end 0376 0377 return