schur_matrix

 SCHUR_MATRIX Computes  Schur complement matrix Sigma and its preconditioner

 Computes  Schur complement matrix Sigma and, if param(1)~=1, 
 the preconditioner for Sigma

    [Sigma,PNN]=schur_matrix(ifro,nov,wx,dx,jacx,...
    wy,dy,jacy,nvli,gam,novg,lint,lgamma,D,Rgamma,param);

 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
               subdomain
        novg = restriction map from (\cup_i \Omega_i) to \Gamma
                 set in cosnovg.m
        lint= list of internal nodes
        lgamma= list of internal nodes which are on the interface between
                spectral elements
        D, column array of size ngamma (global number of interface nodes)
        Rgamma = matrix contructed in cosrgam.m
        param = array of parameters
        
 Output: Sigma = Schur complement matrix
         PNN = empty array if param(1)==1,
               Neumann-Neumann preconditioner for Sigma, if param(1)==2,
               bal Neumann-Neumann preconditioner for Sigma, if param(1)==3
       

 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 [Sigma,PNN]=schur_matrix(ifro,nov,wx,dx,jacx,...
0002     wy,dy,jacy,nvli,gam,novg,lint,lgamma,D,Rgamma,param);
0003 % SCHUR_MATRIX Computes  Schur complement matrix Sigma and its preconditioner
0004 %
0005 % Computes  Schur complement matrix Sigma and, if param(1)~=1,
0006 % the preconditioner for Sigma
0007 %
0008 %    [Sigma,PNN]=schur_matrix(ifro,nov,wx,dx,jacx,...
0009 %    wy,dy,jacy,nvli,gam,novg,lint,lgamma,D,Rgamma,param);
0010 %
0011 % Input:
0012 %         ifro = column array of length noe=nov(npdx*npdy,ne):
0013 %            if (x_i,y_i) is internal to Omega then ifro(i)=0,
0014 %            if (x_i,y_i) is on \partial\Omega then ifro(i)=1,
0015 %        nov = local -global map, previously generated by cosnov_2d
0016 %        wx = npdx LGL weigths in [-1,1], previously generated by xwlgl
0017 %        dx =first derivative LGL matrix (by calling dx=derlgl(x,npdx))
0018 %        jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
0019 %        wy = npdy LGL weigths in [-1,1],
0020 %            (produced by calling [y,wy]=xwlgl(npdy))
0021 %        dy =first derivative LGL matrix (by calling dy=derlgl(y,npdy))
0022 %        jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
0023 %        nvli = column array with number of internal nodes of each
0024 %               subdomain
0025 %        novg = restriction map from (\cup_i \Omega_i) to \Gamma
0026 %                 set in cosnovg.m
0027 %        lint= list of internal nodes
0028 %        lgamma= list of internal nodes which are on the interface between
0029 %                spectral elements
0030 %        D, column array of size ngamma (global number of interface nodes)
0031 %        Rgamma = matrix contructed in cosrgam.m
0032 %        param = array of parameters
0033 %
0034 % Output: Sigma = Schur complement matrix
0035 %         PNN = empty array if param(1)==1,
0036 %               Neumann-Neumann preconditioner for Sigma, if param(1)==2,
0037 %               bal Neumann-Neumann preconditioner for Sigma, if param(1)==3
0038 %
0039 %
0040 % References: CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0041 %                    "Spectral Methods. Evolution to Complex Geometries
0042 %                     and Applications to Fluid DynamicsSpectral Methods"
0043 %                    Springer Verlag, Berlin Heidelberg New York, 2007.
0044 
0045 %   Written by Paola Gervasio
0046 %   $Date: 2007/04/01$
0047 
0048 [ldnov,ne]=size(nov);
0049 noe=nov(ldnov,ne); 
0050 npdx=length(wx); npdy=length(wy); mn=npdx*npdy;
0051 n=length(lgamma);
0052 Sigma=sparse(n,n);PSH=[];
0053 if(param(1)==1)
0054     PNN=[];
0055 else
0056     PNN=sparse(n,n);
0057 end
0058 
0059 [wx1,wy1]=meshgrid(wx,wy); w_ie=wx1.*wy1; w_ie=w_ie'; w_ie=w_ie(:); clear wx1 wy1;
0060 if(param(1)==1) % P=I
0061     
0062 for ie=1:ne
0063 Al=sparse(mn,mn);
0064 [Al]=stiff_2d_sp(wx,dx,jacx(ie),wy,dy,jacy(ie));
0065 jac=jacx(ie)*jacy(ie);
0066 
0067 Al=Al+jac*gam*spdiags(w_ie,0,mn,mn);
0068 ifro_l=ifro(nov(1:mn,ie));
0069 [lbor,lint]=liste1(ifro_l);
0070 
0071 Ali=Al(lint,lint);
0072 
0073 ifroi=ifro_l(lint);
0074 [lbor,lint,lintint,lg]=liste1(ifroi);
0075 
0076 
0077 Amm=Ali(lintint,lintint); 
0078 Agg=Ali(lg,lg);
0079 Agm=Ali(lg,lintint);
0080 Sigma_m=Agg-Agm*inv(Amm)*Agm';
0081 Sigma(novg(lg,ie),novg(lg,ie))=Sigma(novg(lg,ie),novg(lg,ie))+Sigma_m;
0082 end
0083 
0084 elseif(param(1)==2)  %P= Neumann -Neumann preconditioner
0085 
0086 for ie=1:ne
0087 [Al]=stiff_2d_sp(wx,dx,jacx(ie),wy,dy,jacy(ie));
0088 jac=jacx(ie)*jacy(ie);
0089 
0090 Al=Al+jac*gam*spdiags(w_ie,0,mn,mn);
0091 
0092 ifro_l=ifro(nov(1:mn,ie));
0093 [lbor,lint]=liste1(ifro_l);
0094 Ali=Al(lint,lint);
0095 ifroi=ifro_l(lint);
0096 [lbor,lint,lintint,lg]=liste1(ifroi);
0097 
0098 Amm=Ali(lintint,lintint); 
0099 Agg=Ali(lg,lg);
0100 Agm=Ali(lg,lintint);
0101 
0102 if nvli(ie)==mn
0103 % Assembling preconditioner on the element ie.
0104 % Mass/H^2 is added if the boundary of element ie  and boundary of
0105 % Omega have null intersection
0106 % ww is not multiplied by jac, since Mass/H^2= diag(ww*jac)/jac
0107  Amp=Ali+spdiags(w_ie,0,mn,mn);
0108 else
0109 Amp=Ali;
0110 end
0111 
0112 Sigma_m=Agg-Agm*inv(Amm)*Agm';
0113 Sigma(novg(lg,ie),novg(lg,ie))=Sigma(novg(lg,ie),novg(lg,ie))+Sigma_m;
0114 clear Sigma_m;
0115 
0116 Sigmap_m=inv(Amp(lg,lg)-Amp(lg,lintint)*inv(Amp(lintint,lintint))*Amp(lintint,lg));
0117 D_m=D(novg(lg,ie));
0118 PNN(novg(lg,ie),novg(lg,ie))=...
0119     PNN(novg(lg,ie),novg(lg,ie))+diag(D_m)*Sigmap_m*diag(D_m);
0120 end
0121 
0122 elseif(param(1)==3)  %P= balancing Neumann -Neumann preconditioner
0123 
0124 for ie=1:ne
0125 [Al]=stiff_2d_sp(wx,dx,jacx(ie),wy,dy,jacy(ie));
0126 jac=jacx(ie)*jacy(ie);
0127 
0128 Al=Al+jac*gam*spdiags(w_ie,0,mn,mn);
0129 
0130 ifro_l=ifro(nov(1:mn,ie));
0131 [lbor,lint]=liste1(ifro_l);
0132 Ali=Al(lint,lint);
0133 ifroi=ifro_l(lint);
0134 [lbor,lint,lintint,lg]=liste1(ifroi);
0135 
0136 Amm=Ali(lintint,lintint); 
0137 Agg=Ali(lg,lg);
0138 Agm=Ali(lg,lintint);
0139 
0140 if nvli(ie)==mn
0141 % Assembling preconditioner on the element ie.
0142 % Mass/H^2 is added if the boundary of element ie  and boundary of
0143 % Omega have null intersection
0144 % ww is not multiplied by jac, since Mass/H^2= diag(ww*jac)/jac
0145  Amp=Ali+spdiags(w_ie,0,mn,mn);
0146 else
0147 Amp=Ali;
0148 end
0149 nl=novg(lg,ie);
0150 Sigma_m=Agg-Agm*inv(Amm)*Agm';
0151 Sigma(nl,nl)=Sigma(nl,nl)+Sigma_m;
0152 clear Sigma_m;
0153 
0154 Sigmap_m=inv(Amp(lg,lg)-Amp(lg,lintint)*inv(Amp(lintint,lintint))*Amp(lintint,lg));
0155 D_m=D(nl);
0156 PNN(nl,nl)=PNN(nl,nl)+diag(D_m)*Sigmap_m*diag(D_m);
0157 end
0158 
0159 % A_H=Rgamma*Sigma*Rgamma'
0160 % Sigma_H^+
0161 
0162 PSH=Rgamma'*pinv(full(Rgamma*Sigma*Rgamma'))*Rgamma;
0163 PNN=(speye(n)-PSH*Sigma)*PNN*(speye(n)-Sigma*PSH)+PSH;
0164 
0165 end  % end loop on ie
0166 return
0167