EIG_SCHUR_2D Eigenvalues computation Schur complement matrix
Eigenvalues computation for the matrix associated to
the 2D boundary value problem
-Delta u + gam u = f x in Omega
+ Dirichlet bc
by Schur Complement Matrix (Sect. 6.4.3, CHQZ3)
and by SEM Numerical Integration with LGL quadrature formulas.
[param]=eig_schur_2d(xa,xb,ya,yb,gam,...
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)
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
It works only if cb='dddd';
nex = number of elements (equally spaced) along x-direction
nx = polynomial degree in each element (the same in each element)
along x-direction
ney = number of elements (equally spaced) 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 : No preconditioner (P=I)
2 : P= Neumann-Neumann (Alg. 6.4.3, CHQZ3)
3 : P= balancing Neumann-Neumann (Alg. 6.4.4, CHQZ3)
param(2) : maximum number of iterations in eigs
param(3) : tolerance for stopping test in eigs
Output:
param(11,12) extrema eigenvalues (max,min) of
the preconditioned Schur complement
param(13) iterative condition number (|lam_max|/|lam_min|) of
the preconditioned Schur complement
References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Fundamentals in Single Domains"
Springer Verlag, Berlin Heidelberg New York, 2006.
CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Evolution to Complex Geometries
and Applications to Fluid DynamicsSpectral Methods"
Springer Verlag, Berlin Heidelberg New York, 2007.0001 function [param]=eig_schur_2d(xa,xb,ya,yb,gam,... 0002 cb,nex,nx,ney,ny,gammax,gammay,param); 0003 % EIG_SCHUR_2D Eigenvalues computation Schur complement matrix 0004 % 0005 % Eigenvalues computation for the matrix associated to 0006 % the 2D boundary value problem 0007 % 0008 % -Delta u + gam u = f x in Omega 0009 % 0010 % + Dirichlet bc 0011 % 0012 % by Schur Complement Matrix (Sect. 6.4.3, CHQZ3) 0013 % and by SEM Numerical Integration with LGL quadrature formulas. 0014 % 0015 % [param]=eig_schur_2d(xa,xb,ya,yb,gam,... 0016 % nex,nx,ney,ny,gammax,gammay,param); 0017 % 0018 % Omega=(xa,xb) x (ya,yb) 0019 % 0020 % side 3 0021 % V4 __________ V3 0022 % | | 0023 % | | 0024 % side 4 | Omega | side 2 0025 % | | 0026 % |__________| 0027 % V1 V2 0028 % side 1 0029 % 0030 % 0031 % __________________________ 0032 % | | | | | 0033 % | 3 | 6 | 9 | 12 | Omega and spectral elements 0034 % | | | | | ordering 0035 % __________________________ 0036 % | | | | | 0037 % | 2 | 5 | 8 | 11 | 0038 % | | | | | 0039 % __________________________ 0040 % | | | | | 0041 % | 1 | 4 | 7 | 10 | 0042 % | | | | | 0043 % __________________________ 0044 % 0045 % 0046 % Input: xa= abscissa of either vertex V1 or vertex V4 0047 % xb= abscissa of either vertex V2 or vertex V3 0048 % ya= ordinate of either vertex V1 or vertex V2 0049 % yb= ordinate of either vertex V3 or vertex V4 0050 % gam = coefficient of zeroth order term (constant>=0) 0051 % cb= a string of four elements, each for any side of \partial\Omega. 0052 % cb(i) refers to side number "i" of \partial\Omega 0053 % 'd' character stands for Dirichlet boundary conditions on that side 0054 % 'n' character stands for Neumann boundary conditions on that side 0055 % It works only if cb='dddd'; 0056 % nex = number of elements (equally spaced) along x-direction 0057 % nx = polynomial degree in each element (the same in each element) 0058 % along x-direction 0059 % ney = number of elements (equally spaced) along y-direction 0060 % ny = polynomial degree in each element (the same in each element) 0061 % along y-direction% gammax = column or row array of length nex-1. 0062 % If the deomposition of Omega is not 0063 % uniform, gammax is the vector of position of interfaces between 0064 % spectral elements along x-direction. If gammax=[], uniform 0065 % decomposition is used. 0066 % gammay = column or row array of length ney-1. 0067 % If the deomposition of Omega is not 0068 % uniform, gammay is the vector of position of interfaces between 0069 % spectral elements along y-direction. If gammay=[], uniform 0070 % decomposition is used. 0071 % param = array of parameters 0072 % param(1) = 1 : No preconditioner (P=I) 0073 % 2 : P= Neumann-Neumann (Alg. 6.4.3, CHQZ3) 0074 % 3 : P= balancing Neumann-Neumann (Alg. 6.4.4, CHQZ3) 0075 % param(2) : maximum number of iterations in eigs 0076 % param(3) : tolerance for stopping test in eigs 0077 % 0078 % Output: 0079 % param(11,12) extrema eigenvalues (max,min) of 0080 % the preconditioned Schur complement 0081 % param(13) iterative condition number (|lam_max|/|lam_min|) of 0082 % the preconditioned Schur complement 0083 % 0084 % 0085 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0086 % "Spectral Methods. Fundamentals in Single Domains" 0087 % Springer Verlag, Berlin Heidelberg New York, 2006. 0088 % CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0089 % "Spectral Methods. Evolution to Complex Geometries 0090 % and Applications to Fluid DynamicsSpectral Methods" 0091 % Springer Verlag, Berlin Heidelberg New York, 2007. 0092 0093 % Written by Paola Gervasio 0094 % $Date: 2007/04/01$ 0095 0096 xy=[]; un=[]; 0097 0098 npdx=nx+1; npdy=ny+1; ldnov=npdx*npdy; mn=ldnov; ne=nex*ney; 0099 0100 [x,wx]=xwlgl(npdx); [dx]=derlgl(x,npdx); 0101 [y,wy]=xwlgl(npdy); [dy]=derlgl(y,npdy); 0102 0103 % Generation of nov: its implements the 0104 % restriction map R_m from \overline\Omega to overline\Omega_m 0105 % (CHQZ, pag. 394) 0106 0107 [nov]=cosnov_2d(npdx,nex,npdy,ney); 0108 noe=nov(ldnov,ne); 0109 0110 % Mesh generation 0111 0112 [xx,yy,jacx,jacy,xy,ww,ifro]=mesh_2d(xa,xb,ya,yb,cb,nex,ney,npdx,npdy,... 0113 nov,x,wx,y,wy,gammax,gammay); 0114 0115 0116 % Generation of lists of internal, boundary, interface nodes. All lists are referred to 0117 % global ordering on Omega 0118 % 0119 % ldir: list of Dirichlet boundary nodes 0120 % lint: list of internal nodes 0121 % lintint: list of internal nodes which are internal to spectral elements 0122 % lgamma: list of internal nodes which are on the interfaces between spectral elements 0123 0124 [ldir,lint,lintint,lgamma,ifro]=liste(ifro,nov); 0125 [novi,nvli]=cosnovi(nov,ifro,lint); 0126 0127 0128 noei=length(lint); 0129 ifroi=ifro(lint); 0130 xyi=[xy(lint,1),xy(lint,2)]; 0131 wwi=ww(lint); 0132 0133 % Generation of novg: its implements the 0134 % restriction map R_{\Gamma_m} from \Gamma to \Gamma_m (= 0135 % \partial\Omega_m\cap\Gamma) (CHQZ, pag. 394) 0136 0137 ngamma=length(lgamma); 0138 [novg]=cosnovg(xyi,noei,ifroi,lgamma,ldnov,novi,nvli); 0139 0140 % Construction of LGG, to set Rgamma 0141 0142 LGG=cell(ne,1); 0143 for ie=1:ne 0144 ifro_l=ifro(nov(1:mn,ie)); 0145 [lbor,lint]=liste1(ifro_l); 0146 ifroi=ifro_l(lint); 0147 [lbor,lint,lintint,lg]=liste1(ifroi); 0148 LGG(ie)={lg}; 0149 end 0150 0151 % Rgamma = matrix of size (ne, ngamma) (defined in CHQZ3, pag. 399) 0152 % (R_\Gamma)_ij: = 1/n_j if x_j \in \Gamma_i 0153 % = 0 otherwise 0154 0155 [Rgamma]=cosrgam(novg,LGG,ne,ngamma); 0156 0157 % Computes the diagonal weighting matrix D relative to the interface 0158 % unknowns: 0159 % D_ii=1/n_i where n_i is the number of subdomains x_i belongs to. 0160 % D_m=D(novg(LGG_m,m)) is the restriction of D to Gamma_m 0161 % (used in (6.4.40), pag. 397, CHQZ3 0162 0163 0164 if (param(1)==1) 0165 D=[]; 0166 else 0167 D=partition(Rgamma); 0168 end 0169 0170 % Generation of both Schur complement matrix and preconditioner 0171 0172 [Sigma,PNN]=schur_matrix(ifro,nov,wx,dx,jacx,... 0173 wy,dy,jacy,nvli,gam,novg,lint,lgamma,D,Rgamma,param); 0174 0175 0176 if(param(1)~=1) % Neumann Neumann o balancing Neumann-Neumann 0177 Sigma=PNN*Sigma; 0178 end 0179 0180 if param(2)==0 0181 param(2)=1000; 0182 end 0183 if param(3)==0 0184 param(3)=1.d-12; 0185 end 0186 opts.issym=1;opts.disp=0; 0187 opts.maxit=param(2); opts.tol=param(3); 0188 lambda_max=eigs(Sigma,1,'lm',opts); 0189 lambda_min=eigs(Sigma,1,'sm',opts); 0190 param(11)=lambda_max; 0191 param(12)=lambda_min; 0192 param(13)=lambda_max/lambda_min; 0193 0194 0195 return