SCHUR_LOC Computes local matrices and lists for implementing the Schur method
Computes local matrices and lists for implementing Schur complement algorithm
Input : ifro = column array of length noe=nov(npdx*npdy,ne):
if (x_i,y_i) is internal to Omega then ifro(i)=0,
if (x_i,y_i) is on \partial\Omega then ifro(i)=1,
nov = local -global map, previously generated by cosnov_2d
wx = npdx LGL weigths in [-1,1], previously generated by xwlgl
dx =first derivative LGL matrix (by calling dx=derlgl(x,npdx))
jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
wy = npdy LGL weigths in [-1,1],
(produced by calling [y,wy]=xwlgl(npdy))
dy =first derivative LGL matrix (by calling dy=derlgl(y,npdy))
jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
nvli = column array with number of internal nodes of each
gam = coefficient of zeroth order term (constant>=0)
f = r.h.s. column array
ub = column array with the solution on the boundary
param = array of parameters. only param(1) is used.
Output: AGG = structure of local matrices A_{\Gamma_m,\Gamma_m}
Amm = structure of local matrices (A_{m,m})^{-1}
AGm = structure of local matrices A_{\Gamma_m,m}
Lmm = structure containing maps from internal{\Omega_m to \Omega_m
set in schur_local.m
LGG = structure with maps from \Gamma_i to \Gamma
set in schur_local.m
novg = restriction map from (\cup_i \Omega_i) to \Gamma
set in cosnovg.m
Rgamma = matrix of size (ne, ngamma) (defined in CHQZ3, pag. 399)
(R_\Gamma)_ij: = 1/n_j if x_j \in \Gamma_i
= 0 otherwise
References: 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 [AGG,Amm,AGm,Lmm,LGG,Am,f]=schur_loc(ifro,nov,... 0002 wx,dx,jacx,wy,dy,jacy,nvli,gam,f,ub,param); 0003 % SCHUR_LOC Computes local matrices and lists for implementing the Schur method 0004 % 0005 % Computes local matrices and lists for implementing Schur complement algorithm 0006 % 0007 % Input : ifro = column array of length noe=nov(npdx*npdy,ne): 0008 % if (x_i,y_i) is internal to Omega then ifro(i)=0, 0009 % if (x_i,y_i) is on \partial\Omega then ifro(i)=1, 0010 % nov = local -global map, previously generated by cosnov_2d 0011 % wx = npdx LGL weigths in [-1,1], previously generated by xwlgl 0012 % dx =first derivative LGL matrix (by calling dx=derlgl(x,npdx)) 0013 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0014 % wy = npdy LGL weigths in [-1,1], 0015 % (produced by calling [y,wy]=xwlgl(npdy)) 0016 % dy =first derivative LGL matrix (by calling dy=derlgl(y,npdy)) 0017 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0018 % nvli = column array with number of internal nodes of each 0019 % gam = coefficient of zeroth order term (constant>=0) 0020 % f = r.h.s. column array 0021 % ub = column array with the solution on the boundary 0022 % param = array of parameters. only param(1) is used. 0023 % 0024 % Output: AGG = structure of local matrices A_{\Gamma_m,\Gamma_m} 0025 % Amm = structure of local matrices (A_{m,m})^{-1} 0026 % AGm = structure of local matrices A_{\Gamma_m,m} 0027 % Lmm = structure containing maps from internal{\Omega_m to \Omega_m 0028 % set in schur_local.m 0029 % LGG = structure with maps from \Gamma_i to \Gamma 0030 % set in schur_local.m 0031 % novg = restriction map from (\cup_i \Omega_i) to \Gamma 0032 % set in cosnovg.m 0033 % Rgamma = matrix of size (ne, ngamma) (defined in CHQZ3, pag. 399) 0034 % (R_\Gamma)_ij: = 1/n_j if x_j \in \Gamma_i 0035 % = 0 otherwise 0036 % 0037 % References: CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0038 % "Spectral Methods. Evolution to Complex Geometries 0039 % and Applications to Fluid DynamicsSpectral Methods" 0040 % Springer Verlag, Berlin Heidelberg New York, 2007. 0041 0042 % Written by Paola Gervasio 0043 % $Date: 2007/04/01$ 0044 0045 [ldnov,ne]=size(nov); 0046 noe=nov(ldnov,ne); 0047 npdx=length(wx); npdy=length(wy); mn=npdx*npdy; 0048 0049 0050 Amm=cell(ne,1); AGG=cell(ne,1); AGm=cell(ne,1); 0051 Lmm=cell(ne,1);LGG=cell(ne,1); 0052 Am=cell(ne,1); 0053 [wx1,wy1]=meshgrid(wx,wy); w_ie=wx1.*wy1; w_ie=w_ie'; w_ie=w_ie(:); clear wx1 wy1; 0054 0055 for ie=1:ne 0056 Al=sparse(mn,mn); 0057 [Al]=stiff_2d_sp(wx,dx,jacx(ie),wy,dy,jacy(ie)); 0058 jac=jacx(ie)*jacy(ie); 0059 0060 Al=Al+jac*gam*spdiags(w_ie,0,mn,mn); 0061 0062 % Ali is the block of Al relative to internal nodes 0063 0064 ifro_l=ifro(nov(1:mn,ie)); 0065 ub_loc=ub(nov(1:mn,ie)); 0066 [lbor,lint]=liste1(ifro_l); 0067 if(isempty(lbor)==0) 0068 f(nov(lint,ie))=f(nov(lint,ie))-Al(lint,lbor)*ub_loc(lbor); 0069 end 0070 0071 Ali=Al(lint,lint); 0072 ifroi=ifro_l(lint); 0073 [lbor,lint,lintint,lg]=liste1(ifroi); 0074 Lmm(ie)={lintint}; LGG(ie)={lg}; 0075 ngammam=length(lg); 0076 AGG(ie)={Ali(lg,lg)}; 0077 AGm(ie)={Ali(lg,lintint)}; 0078 Amm(ie)={inv(Ali(lintint,lintint))}; 0079 0080 % Assemble preconditioner on the element ie. 0081 % Mass/H^2 is added if the boundary of element ie and boundary of 0082 % Omega have null intersection 0083 0084 if(param(1)~=1) 0085 if nvli(ie)==mn 0086 0087 % ww is not multiplied by jac, since Mass/H^2= diag(ww*jac)/jac 0088 0089 Ali=Ali+spdiags(w_ie,0,mn,mn); 0090 end 0091 Am(ie)={sparse(chol(Ali))}; 0092 end 0093 0094 % end loop on ie 0095 end 0096