MESH3D Uniform 3D Spectral element mesh on parallelepipedon
Omega=(xa,xb) x (ya,yb) x (za,zb)
V8 _________ V7
/| /|
/ | / |
V5 /________ /V6|
| | | |
|V4|______|__| V3
| / | /
| / |/
|/________/
V1 V2
__________________________
| | | | |
| 3 | 6 | 9 | 12 | Spectral elements
| | | | | ordering in a plane yz, then
__________________________ planes at different x follow.
| | | | |
| 2 | 5 | 8 | 11 |
| | | | |
__________________________
| | | | |
| 1 | 4 | 7 | 10 |
| | | | |
__________________________
[xx,yy,zz,jacx,jacy,jacz,xyz,ww,ifro,nov]=mesh3d(xa,xb,ya,yb,za,zb,...
nex,ney,nez,npdx,npdy,npdz,nov,x,wx,y,wy,z,wz,gammax,gammay,gammaz);
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
za= ordinate of either vertex V1
zb= ordinate of either vertex V5
nex = number of elements along x-direction
ney = number of elements along y-direction
nez = number of elements (equally spaced) along z-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
npdz = number of nodes in each element (the same in every element)
along z-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
z = npdz LGL nodes in [-1,1], previously generated by xwlgl
wz = npdz 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 y-direction. If gammay=[], uniform
decomposition is used.
gammaz = column or row array of length nez-1.
If the deomposition of Omega is not
uniform, gammay is the vector of position of interfaces between
spectral elements along z-direction. If gammay=[], uniform
decomposition is used.
Output: xx = 2-indexes array of size (8,ne)
xx(1:8,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie;...
x_V5_ie;x_V6_ie;x_V7_ie;x_V8_ie]
(ne=nex*ney*nez is the global number of spectral elements)
yy(1:8,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie;...
y_V5_ie;y_V6_ie;y_V7_ie;y_V8_ie]
zz(1:8,ie)=[z_V1_ie;z_V2_ie;z_V3_ie;z_V4_ie;...
z_V5_ie;z_V6_ie;z_V7_ie;z_V8_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
jacz = array (length(jacz)=ne); jacz(ie)= (z_V5_ie-z_V1_ie)/2
xyz = column array with global mesh, length: noe=nov(npdx*npdy*npdz,ne)
ww = column array with global weigths, length: noe=nov(npdx*npdy*npdz,ne)
diag(ww) is the mass matrix associated to SEM discretization
ifro = column array of length noe=nov(npdx*npdy*npdz,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 = 2-index array of local to global map,
size(nov)=[max(npdx*npdy*npdz),ne]
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,zz,jacx,jacy,jacz,xyz,ww,ifro,nov]=mesh3d(xa,xb,ya,yb,za,zb,... 0002 nex,ney,nez,npdx,npdy,npdz,x,wx,y,wy,z,wz,gammax,gammay,gammaz); 0003 % MESH3D Uniform 3D Spectral element mesh on parallelepipedon 0004 % 0005 % Omega=(xa,xb) x (ya,yb) x (za,zb) 0006 % 0007 % V8 _________ V7 0008 % /| /| 0009 % / | / | 0010 % V5 /________ /V6| 0011 % | | | | 0012 % |V4|______|__| V3 0013 % | / | / 0014 % | / |/ 0015 % |/________/ 0016 % V1 V2 0017 % 0018 % 0019 % 0020 % __________________________ 0021 % | | | | | 0022 % | 3 | 6 | 9 | 12 | Spectral elements 0023 % | | | | | ordering in a plane yz, then 0024 % __________________________ planes at different x follow. 0025 % | | | | | 0026 % | 2 | 5 | 8 | 11 | 0027 % | | | | | 0028 % __________________________ 0029 % | | | | | 0030 % | 1 | 4 | 7 | 10 | 0031 % | | | | | 0032 % __________________________ 0033 % 0034 % [xx,yy,zz,jacx,jacy,jacz,xyz,ww,ifro,nov]=mesh3d(xa,xb,ya,yb,za,zb,... 0035 % nex,ney,nez,npdx,npdy,npdz,nov,x,wx,y,wy,z,wz,gammax,gammay,gammaz); 0036 % 0037 % Input: xa= abscissa of either vertex V1 or vertex V4 0038 % xb= abscissa of either vertex V2 or vertex V3 0039 % ya= ordinate of either vertex V1 or vertex V2 0040 % yb= ordinate of either vertex V3 or vertex V4 0041 % za= ordinate of either vertex V1 0042 % zb= ordinate of either vertex V5 0043 % nex = number of elements along x-direction 0044 % ney = number of elements along y-direction 0045 % nez = number of elements (equally spaced) along z-direction 0046 % npdx = number of nodes in each element (the same in every element) 0047 % along x-direction 0048 % npdy = number of nodes in each element (the same in every element) 0049 % along y-direction 0050 % npdz = number of nodes in each element (the same in every element) 0051 % along z-direction 0052 % nov = local -global map, previously generated by cosnov_2d 0053 % x = npdx LGL nodes in [-1,1], previously generated by xwlgl 0054 % wx = npdx LGL weigths in [-1,1], previously generated by xwlgl 0055 % y = npdy LGL nodes in [-1,1], previously generated by xwlgl 0056 % wy = npdy LGL weigths in [-1,1], previously generated by xwlgl 0057 % z = npdz LGL nodes in [-1,1], previously generated by xwlgl 0058 % wz = npdz LGL weigths in [-1,1], previously generated by xwlgl 0059 % gammax = column or row array of length nex-1. 0060 % If the deomposition of Omega is not 0061 % uniform, gammax is the vector of position of interfaces between 0062 % spectral elements along x-direction. If gammax=[], uniform 0063 % decomposition is used. 0064 % gammay = column or row array of length ney-1. 0065 % If the deomposition of Omega is not 0066 % uniform, gammay is the vector of position of interfaces between 0067 % spectral elements along y-direction. If gammay=[], uniform 0068 % decomposition is used. 0069 % gammaz = column or row array of length nez-1. 0070 % If the deomposition of Omega is not 0071 % uniform, gammay is the vector of position of interfaces between 0072 % spectral elements along z-direction. If gammay=[], uniform 0073 % decomposition is used. 0074 % 0075 % Output: xx = 2-indexes array of size (8,ne) 0076 % xx(1:8,ie)=[x_V1_ie;x_V2_ie;x_V3_ie;x_V4_ie;... 0077 % x_V5_ie;x_V6_ie;x_V7_ie;x_V8_ie] 0078 % (ne=nex*ney*nez is the global number of spectral elements) 0079 % yy(1:8,ie)=[y_V1_ie;y_V2_ie;y_V3_ie;y_V4_ie;... 0080 % y_V5_ie;y_V6_ie;y_V7_ie;y_V8_ie] 0081 % zz(1:8,ie)=[z_V1_ie;z_V2_ie;z_V3_ie;z_V4_ie;... 0082 % z_V5_ie;z_V6_ie;z_V7_ie;z_V8_ie] 0083 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0084 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0085 % jacz = array (length(jacz)=ne); jacz(ie)= (z_V5_ie-z_V1_ie)/2 0086 % xyz = column array with global mesh, length: noe=nov(npdx*npdy*npdz,ne) 0087 % ww = column array with global weigths, length: noe=nov(npdx*npdy*npdz,ne) 0088 % diag(ww) is the mass matrix associated to SEM discretization 0089 % ifro = column array of length noe=nov(npdx*npdy*npdz,ne): 0090 % if (x_i,y_i) is internal to Omega then ifro(i)=0, 0091 % if (x_i,y_i) is on \partial\Omega then ifro(i)=1, 0092 % nov = 2-index array of local to global map, 0093 % size(nov)=[max(npdx*npdy*npdz),ne] 0094 % 0095 % 0096 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0097 % "Spectral Methods. Fundamentals in Single Domains" 0098 % Springer Verlag, Berlin Heidelberg New York, 2006. 0099 0100 % Written by Paola Gervasio 0101 % $Date: 2007/04/01$ 0102 0103 mn=npdx*npdy; ldnov=mn*npdz; ne=nex*ney*nez; 0104 xx=zeros(8,ne); yy=zeros(8,ne); zz=zeros(8,ne); 0105 jac=zeros(ne,1); 0106 jacx=zeros(ne,1); 0107 jacy=zeros(ne,1); 0108 jacz=zeros(ne,1); 0109 if sum(size(gammax))==0 0110 Hx=(xb-xa)/nex*ones(nex); 0111 else 0112 Hx=zeros(nex,1); 0113 Hx(1)=gammax(1)-xa; 0114 Hx(2:nex-1)=gammax(2:nex-1)-gammax(1:nex-2); 0115 Hx(nex)=xb-gammax(nex-1); 0116 end 0117 if sum(size(gammay))==0 0118 Hy=(yb-ya)/ney*ones(ney); 0119 else 0120 Hy=zeros(ney,1); 0121 Hy(1)=gammay(1)-ya; 0122 Hy(2:ney-1)=gammay(2:ney-1)-gammay(1:ney-2); 0123 Hy(ney)=yb-gammay(ney-1); 0124 end 0125 if sum(size(gammaz))==0 0126 Hz=(zb-za)/nez*ones(nez); 0127 else 0128 Hz=zeros(nez,1); 0129 Hz(1)=gammaz(1)-za; 0130 Hz(2:nez-1)=gammaz(2:nez-1)-gammaz(1:nez-2); 0131 Hz(nez)=zb-gammaz(nez-1); 0132 end 0133 0134 % construction of xx, yy, zz, jacx, jacy, jacz 0135 0136 0137 % first element 0138 xx(1:4,1)=[xa;xa+Hx(1);xa+Hx(1);xa]; 0139 xx(5:8,1)=xx(1:4,1); 0140 yy(1:2,1)=[ya,ya]; yy(3:4,1)=yy(1:2,1)+Hy(1); 0141 yy(5:6,1)=yy(1:2,1);yy(7:8,1)=yy(3:4,1); 0142 zz(1:4,1)=[za;za;za;za]; 0143 zz(5:8,1)=[za+Hz(1);za+Hz(1);za+Hz(1);za+Hz(1)]; 0144 jacx(1)=Hx(1)*.5;jacy(1)=Hy(1)*.5; jacz(1)=Hz(1)*.5; 0145 0146 % first row of elements along y-direction 0147 for iey=2:ney 0148 ie=iey; 0149 xx(1:8,ie)=xx(1:8,1); 0150 yy(1:2,ie)=yy(3:4,iey-1); yy(3:4,ie)=yy(1:2,ie)+Hy(iey); 0151 yy(5:6,ie)=yy(1:2,ie);yy(7:8,ie)=yy(3:4,ie); 0152 zz(1:4,ie)=zz(1:4,1); zz(5:8,ie)=zz(1:4,ie)+Hz(1); 0153 jacx(ie)=Hx(1)*.5; 0154 jacy(ie)=Hy(iey)*.5; 0155 jacz(ie)=Hz(1)*.5; 0156 end 0157 0158 % other elements on the plane xz 0159 0160 for iez=2:nez 0161 for iey=1:ney 0162 ie=(iez-1)*ney+iey; 0163 0164 xx(1:8,ie)=xx(1:8,1); 0165 yy(1:8,ie)=yy(1:8,iey); 0166 zz(1:4,ie)=zz(5:8,ie-ney); zz(5:8,ie)=zz(1:4,ie)+Hz(iez); 0167 jacx(ie)=Hx(1)*.5; 0168 jacy(ie)=Hy(iey)*.5; 0169 jacz(ie)=Hz(iez)*.5; 0170 end 0171 end 0172 0173 % the whole mesh 0174 0175 ne2=ney*nez; 0176 for iex=2:nex 0177 for iez=1:nez 0178 for iey=1:ney 0179 ie=(iex-1)*ne2+(iez-1)*ney+iey; 0180 xx([1;4;5;8],ie)=xx([2;3;6;7],ie-ne2); 0181 xx([2;3;6;7],ie)=xx([1;4;5;8],ie)+Hx(iex); 0182 yy(1:8,ie)=yy(1:8,iey); 0183 zz(1:8,ie)=zz(1:8,ie-ne2); 0184 jacx(ie)=Hx(iex)*.5; 0185 jacy(ie)=Hy(iey)*.5; 0186 jacz(ie)=Hz(iez)*.5; 0187 end 0188 end 0189 end 0190 0191 % construction of all local meshes 0192 0193 xyz_loc=zeros(ldnov,3,ne); 0194 ifro_loc=zeros(ldnov,ne); 0195 ww_loc=zeros(ldnov,ne); 0196 for ie=1:ne 0197 bpa=(xx(1,ie)+xx(2,ie))*.5; 0198 dpc=(yy(1,ie)+yy(3,ie))*.5; 0199 epf=(zz(1,ie)+zz(5,ie))*.5; 0200 jac=jacx(ie)*jacy(ie)*jacz(ie); 0201 for l=1:npdz 0202 for j=1:npdy 0203 for i=1:npdx 0204 k=(l-1)*mn+(j-1)*npdx+i; 0205 xyz_loc(k,1,ie)=x(i)*jacx(ie)+bpa; 0206 xyz_loc(k,2,ie)=y(j)*jacy(ie)+dpc; 0207 xyz_loc(k,3,ie)=z(l)*jacz(ie)+epf; 0208 ww_loc(k,ie)=wz(l)*wx(i)*wy(j)*jac; 0209 end 0210 end 0211 0212 end 0213 0214 for l=1:npdz 0215 for j=1:npdy 0216 for i=1:npdx 0217 k=(l-1)*mn+(j-1)*npdx+i; 0218 if (l==1 | l==npdz | j==1 | j==npdy | i==1 | i==npdx) 0219 ifro_loc(k,ie)=1; 0220 end 0221 end 0222 end 0223 end 0224 0225 end 0226 0227 0228 % assembling global and unique mesh 0229 xyz=zeros(ldnov*ne,3); 0230 nov=zeros(ldnov,ne); 0231 xyz(1:ldnov,1:3)=xyz_loc(1:ldnov,1:3,1); 0232 nov(1:ldnov,1)=(1:ldnov)'; 0233 ifro(1:ldnov)=ifro_loc(1:ldnov,1); 0234 epsi=1.e-14; 0235 last=ldnov; 0236 for ie=2:ne 0237 for i=1:ldnov 0238 punto_new(1:3)=xyz_loc(i,1:3,ie); 0239 0240 tt=0;j=1; 0241 while (j<=last & tt==0) 0242 if (ifro(j)~=0) 0243 0244 dd=abs(punto_new-xyz(j,1:3)); 0245 if (dd<=epsi) 0246 nov(i,ie)=j; 0247 tt=1; 0248 end 0249 end 0250 j=j+1; 0251 end 0252 if (nov(i,ie)==0) 0253 last=last+1; 0254 xyz(last,1:3)=punto_new(1:3); 0255 ifro(last)=ifro_loc(i,ie); 0256 nov(i,ie)=last; 0257 end 0258 0259 end 0260 end 0261 0262 ifro=ifro'; 0263 noe=last; 0264 xyz=xyz(1:noe,1:3); 0265 ww=zeros(noe,1); 0266 0267 for ie=1:ne 0268 ww(nov(1:ldnov,ie))=ww(nov(1:ldnov,ie))+ww_loc(1:ldnov,ie); 0269 end 0270 clear ww_loc xyz_loc ifro_loc 0271 0272 for i=1:noe 0273 if (ifro(i)~=0) 0274 if((abs(xyz(i,1)-xa) > eps & abs(xyz(i,1)-xb) > eps)... 0275 & (abs(xyz(i,2)-ya) > eps & abs(xyz(i,2)-yb) > eps)... 0276 & (abs(xyz(i,3)-za) > eps & abs(xyz(i,3)-zb) > eps)) 0277 ifro(i)=0; 0278 end 0279 end 0280 end 0281 0282