MESH2D Uniform 2D SEM mesh on rectangular domain Omega=(xa,xb) x (ya,yb)
side 3
V4 __________ V3
| |
| |
side 4 | Omega | side 2
| |
|__________|
V1 V2
side 1
[xx,yy,jacx,jacy,xy,ww,ifro,xn]=mesh2d(xa,xg,ya,yb,nex,ney,npdx,npdy,...
nov,x,wx,y,wy)
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
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
ney = number of elements along y-direction
npdx = number of nodes in each element (the same in every element)
along x-direction
npdy = number of nodes in each element (the same in every element)
along y-direction
nov = local -global map, previously generated by cosnov_2d
x = npdx LGL nodes in [-1,1], previously generated by xwlgl
wx = npdx LGL weigths in [-1,1], previously generated by xwlgl
y = npdy LGL nodes in [-1,1], previously generated by xwlgl
wy = npdy LGL weigths in [-1,1], previously generated by xwlgl
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 x-direction. If gammay=[], uniform
decomposition is used.
Output: 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)
yy = 2-indexes array of size (4,ne):
yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie]
jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2
jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2
xy = column array with global mesh, length: noe=nov(npdx*npdy,ne)
ww = column array with global weigths, length: noe=nov(npdx*npdy,ne)
diag(ww) is the mass matrix associated to SEM discretization
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 and Dirichlet boundary
condition is imposed there, then ifro(i)=1,
if (x_i,y_i) is on \partial\Omega and Neumann boundary
condition is imposed there, then ifro(i)=31,
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[xx,yy,jacx,jacy,xy,ww,ifro]=mesh2d(xa,xb,ya,yb,cb,nex,ney,npdx,npdy,... 0002 nov,x,wx,y,wy,gammax,gammay) 0003 % MESH2D Uniform 2D SEM mesh on rectangular domain Omega=(xa,xb) x (ya,yb) 0004 % 0005 % side 3 0006 % V4 __________ V3 0007 % | | 0008 % | | 0009 % side 4 | Omega | side 2 0010 % | | 0011 % |__________| 0012 % V1 V2 0013 % side 1 0014 % 0015 % [xx,yy,jacx,jacy,xy,ww,ifro,xn]=mesh2d(xa,xg,ya,yb,nex,ney,npdx,npdy,... 0016 % nov,x,wx,y,wy) 0017 % 0018 % Input: xa= abscissa of either vertex V1 or vertex V4 0019 % xb= abscissa of either vertex V2 or vertex V3 0020 % ya= ordinate of either vertex V1 or vertex V2 0021 % yb= ordinate of either vertex V3 or vertex V4 0022 % cb= a string of four elements, each for any side of \partial\Omega. 0023 % cb(i) refers to side number "i" of \partial\Omega 0024 % 'd' character stands for Dirichlet boundary conditions on that side 0025 % 'n' character stands for Neumann boundary conditions on that side 0026 % nex = number of elements along x-direction 0027 % ney = number of elements along y-direction 0028 % npdx = number of nodes in each element (the same in every element) 0029 % along x-direction 0030 % npdy = number of nodes in each element (the same in every element) 0031 % along y-direction 0032 % nov = local -global map, previously generated by cosnov_2d 0033 % x = npdx LGL nodes in [-1,1], previously generated by xwlgl 0034 % wx = npdx LGL weigths in [-1,1], previously generated by xwlgl 0035 % y = npdy LGL nodes in [-1,1], previously generated by xwlgl 0036 % wy = npdy LGL weigths in [-1,1], previously generated by xwlgl 0037 % gammax = column or row array of length nex-1. 0038 % If the deomposition of Omega is not 0039 % uniform, gammax is the vector of position of interfaces between 0040 % spectral elements along x-direction. If gammax=[], uniform 0041 % decomposition is used. 0042 % gammay = column or row array of length ney-1. 0043 % If the deomposition of Omega is not 0044 % uniform, gammay is the vector of position of interfaces between 0045 % spectral elements along x-direction. If gammay=[], uniform 0046 % decomposition is used. 0047 % 0048 % Output: xx = 2-indexes array of size (4,ne) 0049 % xx(1:4,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie] 0050 % (ne=nex*ney is the global number of spectral elements) 0051 % yy = 2-indexes array of size (4,ne): 0052 % yy(1:4,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie] 0053 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0054 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0055 % xy = column array with global mesh, length: noe=nov(npdx*npdy,ne) 0056 % ww = column array with global weigths, length: noe=nov(npdx*npdy,ne) 0057 % diag(ww) is the mass matrix associated to SEM discretization 0058 % ifro = column array of length noe=nov(npdx*npdy,ne): 0059 % if (x_i,y_i) is internal to Omega then ifro(i)=0, 0060 % if (x_i,y_i) is on \partial\Omega and Dirichlet boundary 0061 % condition is imposed there, then ifro(i)=1, 0062 % if (x_i,y_i) is on \partial\Omega and Neumann boundary 0063 % condition is imposed there, then ifro(i)=31, 0064 % 0065 % 0066 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0067 % "Spectral Methods. Fundamentals in Single Domains" 0068 % Springer Verlag, Berlin Heidelberg New York, 2006. 0069 0070 % Written by Paola Gervasio 0071 % $Date: 2007/04/01$ 0072 0073 ldnov=npdx*npdy; ne=nex*ney; 0074 noe=nov(ldnov,ne); nx=npdx-1; 0075 xx=zeros(4,ne); yy=zeros(4,ne); 0076 xy=zeros(noe,2); ww=zeros(noe,1); 0077 jac=zeros(ne,1); 0078 jacx=zeros(ne,1); 0079 jacy=zeros(ne,1); 0080 if sum(size(gammax))==0 0081 Hx=(xb-xa)/nex*ones(nex); 0082 else 0083 Hx=zeros(nex,1); 0084 Hx(1)=gammax(1)-xa; 0085 Hx(2:nex-1)=gammax(2:nex-1)-gammax(1:nex-2); 0086 Hx(nex)=xb-gammax(nex-1); 0087 end 0088 if sum(size(gammay))==0 0089 Hy=(yb-ya)/ney*ones(ney); 0090 else 0091 Hy=zeros(ney,1); 0092 Hy(1)=gammay(1)-ya; 0093 Hy(2:ney-1)=gammay(2:ney-1)-gammay(1:ney-2); 0094 Hy(ney)=yb-gammay(ney-1); 0095 end 0096 0097 for iex=1:nex 0098 ie1=(iex-1)*ney+1; 0099 xx(1,ie1)=xa+sum(Hx(1:iex-1)); 0100 xx(2,ie1)=xx(1,ie1)+Hx(iex); 0101 xx(3,ie1)=xx(2,ie1); xx(4,ie1)=xx(1,ie1); 0102 jacx(ie1)=Hx(iex)*.5; 0103 for iey=2:ney 0104 ie=ie1-1+iey; xx(:,ie)=xx(:,ie1); 0105 jacx(ie)=Hx(iex)*.5; 0106 end 0107 end 0108 0109 0110 for iey=1:ney 0111 yy(1,iey)=ya+sum(Hy(1:iey-1)); 0112 yy(2,iey)=yy(1,iey); 0113 yy(3,iey)=yy(1,iey)+Hy(iey); yy(4,iey)=yy(3,iey); 0114 jacy(iey)=Hy(iey)*.5; 0115 for iex=2:nex 0116 ie=(iex-1)*ney+iey; yy(:,ie)=yy(:,iey); 0117 jacy(ie)=Hy(iey)*.5; 0118 end 0119 end 0120 0121 0122 for ie=1:ne 0123 bpa=(xx(1,ie)+xx(2,ie))*.5; bpc=(yy(1,ie)+yy(3,ie))*.5; 0124 0125 for j=1:npdy 0126 for i=1:npdx 0127 k=(j-1)*npdx+i; ii=nov(k,ie); 0128 xy(ii,1)=x(i)*jacx(ie)+bpa; xy(ii,2)=y(j)*jacy(ie)+bpc; 0129 ww(ii)=ww(ii)+wx(i)*jacx(ie)*wy(j)*jacy(ie); 0130 end 0131 end 0132 0133 end 0134 0135 0136 ifro=zeros(noe,1); 0137 0138 % Loop on elements 0139 for ie=1:ne 0140 0141 % side 1 0142 for i=1:npdx 0143 k=i; ii=nov(k,ie); 0144 if(abs(xy(ii,2)-ya)<=1.e-12 ) 0145 if cb(1)=='d' 0146 ifro(ii)=1; 0147 elseif cb(1)=='n' 0148 ifro(ii)=31; 0149 end 0150 end 0151 end 0152 0153 % side 2 0154 for j=1:npdy 0155 k=j*npdx; ii=nov(k,ie); 0156 if(abs(xy(ii,1)-xb)<=1.e-12 ) 0157 if cb(2)=='d' 0158 ifro(ii)=1; 0159 elseif cb(2)=='n' 0160 ifro(ii)=31; 0161 end 0162 end 0163 end 0164 0165 % side 3 0166 nm=npdx*(npdy-1); 0167 for i=1:npdx 0168 k=i+nm; ii=nov(k,ie); 0169 if(abs(xy(ii,2)-yb)<=1.e-12 ) 0170 if cb(3)=='d' 0171 ifro(ii)=1; 0172 elseif cb(3)=='n' 0173 ifro(ii)=31; 0174 end 0175 end 0176 end 0177 0178 % side 4 0179 for j=1:npdy 0180 k=(j-1)*npdx+1; ii=nov(k,ie); 0181 if(abs(xy(ii,1)-xa)<=1.e-12 ) 0182 if cb(4)=='d' 0183 ifro(ii)=1; 0184 elseif cb(4)=='n' 0185 ifro(ii)=31; 0186 end 0187 end 0188 end 0189 % check vertices: if vertex V belongs to two sides with different (N/D) 0190 %boundary conditions, dirichlet b.c. prevails. 0191 0192 % V1 0193 ip=nov(1,ie); 0194 if ifro(ip)>0 & cb(1)~=cb(4) 0195 ifro(ip)=1; 0196 end 0197 % V2 0198 ip=nov(npdx,ie); 0199 if ifro(ip)>0 & cb(1)~=cb(2) 0200 ifro(ip)=1; 0201 end 0202 % V3 0203 ip=nov(ldnov,ie); 0204 if ifro(ip)>0& cb(2)~=cb(3) 0205 ifro(ip)=1; 0206 end 0207 % V4 0208 ip=nov(ldnov-nx,ie); 0209 if ifro(ip)>0& cb(4)~=cb(3) 0210 ifro(ip)=1; 0211 end 0212 0213 0214 % end loop for ie=1:ne 0215 end