PRECOASC Solves the linear system P_as z= r where P_as is the additive Schwarz prec.
P_as is the additive Schwarz preconditioner with overlapping elements and
coarse mesh (Sect. 6.3.3, CHQZ3)
[z]=precoasc(r,param,noei,lint,p_unity,...
xy, ww, nov, x,wx,y,wy,xx,jacx,yy,jacy,...
Aq1, wwq1, linte, nove, nvle,...
Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc);
Input : r = column array, r.h.s. of the linear system to solve
param = array of parameters
noei = number of nodes internal to Omega
lint = list of internal nodes (set in liste.n)
p_untity = unity partition for extended decomposition, set in partition_e.m
xy = 2-indexes array wiht coordinates of 2D LGL mesh
ww = mass matrix (in column array)
nov = local -global map, previously generated by cosnov_2d
x = npdx LGL nodes in [-1,1],
(produced by calling [x,wx]=xwlgl(npdx))
wx= npdx LGL weigths in [-1,1],
(produced by calling [x,wx]=xwlgl(npdx))
y = npdx LGL nodes in [-1,1],
(produced by calling [y,wy]=xwlgl(npdy))
wy= npdx LGL weigths in [-1,1],
(produced by calling [y,wy]=xwlgl(npdy))
xx = 2-indexes array of size (4,ne)
xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie]
(ne=nex*ney is the global number of spectral elements)
jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
yy = 2-indexes array of size (4,ne):
yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie]
jacy = array (length(jacy)=ne); jacy(ie)= (x_V3_ie-x_V1_ie)/2
Aq1 = cell array with Q1 stiffness matrices (internal/internal nodes)
on extended elements
wwq1 = cell array with Q1 mass matrices (internal nodes)
on extended elements
linte = cell array with list of internal nodes of extended
elements
nove = 2-index array of "extended local" to global map, permuted.
nvle = 2-index array:
nvle (:,1)=number of nodes of the extendend elements
nvle (:,2)=number of Q1 elements of the extendend elements
along x-direction
nvle (:,3)=number of Q1 elements of the extendend elements
along y-direction
Ac = Cholesky factor of the Q1 stiffness matrix
(internal/internal nodes) on the coarse mesh
Acb = Q1 stiffness matrix (internal/boundary nodes)
on the coarse mesh
wwc = Q1 mass matrix (internal nodes)
on the coarse mesh
r0t = matrix R_H^T
lista_coarse = list of nodes of Omega wich are nodes of the coarse mesh
noec = number of nodes of the coarse mesh
novc(4,ne) = it maps each macro-Q1 element on the coarse mesh Q1
lintc = lista of internal nodes of the coarse mesh
ldirc = lista of boudanry nodes of the coarse mesh
Output: z= column array with the solution of P_as z= r0001 function [z]=precoasc(r,param,noei,lint,p_unity,... 0002 xy, ww, nov, x,wx,y,wy,xx,jacx,yy,jacy,... 0003 Aq1, wwq1, linte, nove, nvle,... 0004 Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0005 % PRECOASC Solves the linear system P_as z= r where P_as is the additive Schwarz prec. 0006 % 0007 % P_as is the additive Schwarz preconditioner with overlapping elements and 0008 % coarse mesh (Sect. 6.3.3, CHQZ3) 0009 % 0010 % [z]=precoasc(r,param,noei,lint,p_unity,... 0011 % xy, ww, nov, x,wx,y,wy,xx,jacx,yy,jacy,... 0012 % Aq1, wwq1, linte, nove, nvle,... 0013 % Ac, Acb, wwc, r0t,lista_coarse,noec,novc,lintc,ldirc); 0014 % 0015 % Input : r = column array, r.h.s. of the linear system to solve 0016 % param = array of parameters 0017 % noei = number of nodes internal to Omega 0018 % lint = list of internal nodes (set in liste.n) 0019 % p_untity = unity partition for extended decomposition, set in partition_e.m 0020 % xy = 2-indexes array wiht coordinates of 2D LGL mesh 0021 % ww = mass matrix (in column array) 0022 % nov = local -global map, previously generated by cosnov_2d 0023 % x = npdx LGL nodes in [-1,1], 0024 % (produced by calling [x,wx]=xwlgl(npdx)) 0025 % wx= npdx LGL weigths in [-1,1], 0026 % (produced by calling [x,wx]=xwlgl(npdx)) 0027 % y = npdx LGL nodes in [-1,1], 0028 % (produced by calling [y,wy]=xwlgl(npdy)) 0029 % wy= npdx LGL weigths in [-1,1], 0030 % (produced by calling [y,wy]=xwlgl(npdy)) 0031 % xx = 2-indexes array of size (4,ne) 0032 % xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie] 0033 % (ne=nex*ney is the global number of spectral elements) 0034 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0035 % yy = 2-indexes array of size (4,ne): 0036 % yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie] 0037 % jacy = array (length(jacy)=ne); jacy(ie)= (x_V3_ie-x_V1_ie)/2 0038 % Aq1 = cell array with Q1 stiffness matrices (internal/internal nodes) 0039 % on extended elements 0040 % wwq1 = cell array with Q1 mass matrices (internal nodes) 0041 % on extended elements 0042 % linte = cell array with list of internal nodes of extended 0043 % elements 0044 % nove = 2-index array of "extended local" to global map, permuted. 0045 % nvle = 2-index array: 0046 % nvle (:,1)=number of nodes of the extendend elements 0047 % nvle (:,2)=number of Q1 elements of the extendend elements 0048 % along x-direction 0049 % nvle (:,3)=number of Q1 elements of the extendend elements 0050 % along y-direction 0051 % Ac = Cholesky factor of the Q1 stiffness matrix 0052 % (internal/internal nodes) on the coarse mesh 0053 % Acb = Q1 stiffness matrix (internal/boundary nodes) 0054 % on the coarse mesh 0055 % wwc = Q1 mass matrix (internal nodes) 0056 % on the coarse mesh 0057 % r0t = matrix R_H^T 0058 % lista_coarse = list of nodes of Omega wich are nodes of the coarse mesh 0059 % noec = number of nodes of the coarse mesh 0060 % novc(4,ne) = it maps each macro-Q1 element on the coarse mesh Q1 0061 % lintc = lista of internal nodes of the coarse mesh 0062 % ldirc = lista of boudanry nodes of the coarse mesh 0063 % 0064 % Output: z= column array with the solution of P_as z= r 0065 0066 % 0067 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0068 % "Spectral Methods. Fundamentals in Single Domains" 0069 % Springer Verlag, Berlin Heidelberg New York, 2006. 0070 % CHQZ3 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0071 % "Spectral Methods. Evolution to Complex Geometries 0072 % and Applications to Fluid DynamicsSpectral Methods" 0073 % Springer Verlag, Berlin Heidelberg New York, 2007. 0074 0075 % Written by Paola Gervasio 0076 % $Date: 2007/04/01$ 0077 0078 nx=param(13);npdx=nx+1; 0079 ny=param(14); npdy=ny+1; 0080 mn=npdx*npdy; 0081 nex=param(15); ney=param(16);ne=nex*ney; 0082 noe=length(xy); 0083 0084 rb=zeros(noe,1); 0085 rb(lint)=r; 0086 z=zeros(noei,1); 0087 0088 % local solver 0089 for ie=1:ne 0090 listaie=linte{ie}; mne=nvle(ie,1); 0091 Aloc=Aq1{ie}; 0092 0093 r0=rb(nove(1:mne,ie)); 0094 r0=r0(listaie); 0095 z0=(Aloc\(Aloc'\r0)); 0096 0097 % solution expansion on the global mesh 0098 zc=zeros(mne,1); 0099 zc(listaie)=z0; 0100 zb=zeros(noe,1); 0101 zb(nove(1:mne,ie))=zc; 0102 z=z+zb(lint); 0103 end 0104 0105 clear Aloc; 0106 0107 % solver on the coarse mesh 0108 rb=zeros(noe,1); 0109 rb(lint)=r; 0110 rb=rb.*p_unity; 0111 r0=r0t'*rb; 0112 r0i=r0(lintc)-(Acb*r0(ldirc)); 0113 z0i=(Ac\(Ac'\(r0i))); 0114 z0=r0;z0(lintc)=z0i; 0115 zc=r0t*z0; 0116 zc=zc.*p_unity; 0117 z=z+zc(lint); 0118 return