LAP_3D Numerical solution of the 3D b.v.p. -Delta u + gam u = f
-Delta u + gam u = f x in Omega
+ Dirichlet bc
by SEM Numerical Integration with LGL quadrature formulas.
[xy,un,D,param]=lap_3d(xa,xb,ya,yb,za,zb,gam,...
uex,uex_x,uex_y,uex_z,ff,nex,nx,ney,ny,nez,nz,gammax,gammay,gammaz,...
param);
Input: xa= abscissa of either vertex V1
xb= abscissa of either vertex V2
ya= ordinate of either vertex V1
yb= ordinate of either vertex V3
za= ordinate of either vertex V1
zb= ordinate of either vertex V5
gam = coefficient of zeroth order term (constant>=0)
uex = exact solution (uex=@(x,y,z)[uex(x,y,z)], with .*, .^, ./)
uex_x = exact first x-derivative (uex_x=@(x,y,z)[uex_x(x,y,z)],
with .*, .^, ./)
uex_y = exact first y-derivative (uex_y=@(x,y,z)[uex_y(x,y,z)],
with .*, .^, ./)
uex_z = exact first z-derivative (uex_z=@(x,y,z)[uex_z(x,y,z)],
with .*, .^, ./)
ff = r.h.s. solution (ff=@(x,y,z)[ff(x,y,z)], with .*, .^, ./)
nex = number of elements (equally spaced) along x-direction
nx = polynomial degree in each element (the same in each element)
along x-direction
ney = number of elements (equally spaced) along y-direction
ny = polynomial degree in each element (the same in each element)
along y-direction
nez = number of elements (equally spaced) along z-direction
nz = polynomial degree in each element (the same in each element)
along z-direction
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.
param = array of parameters
param(1) = 1 : SEM-NI approach, see Sect. 5.3.3, CHQZ3
param(2) = 0 : no reordering of the matrix, see pag. 194, CHQZ2
1 : Cuthill-McKee ordering (by symrcm function of
Matlab)
2 : Symmetric minimum degree ordering (by symamd function of
Matlab)
ONLY used if param(3)==1
param(3) = 1 : Solve the differential problem with a direct method
(Cholesky) and compute errors on the exact solution
2 : Compute extrema eigenvalues of stiffness matrix
3 : Compute interface Schur complement, solve the diff.
problem by Schur complement
4: Compute Schur complement and its extrema eigenvalues
5: Compute and plot all the eigenvalues of A (stiffness
matrix)
6: Compute extrema eigenvalues of M^{-1}A, i.e.
the generalized eigenvalues of Av=lambda* M v,
where M is the SEM-NI mass matrix
param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm)
on the exact solution
2: no errors are computed
param(5) = 0: LG quadrature formulas with high precision degree are
used to compute norms (exact norms)
1: LGL quadrature formulas with npdx,npdy nodes are
used to compute norms (discrete norms)
(used only if param(4) == 1)
param(6) = number of nodes for high degree quadrature formula,
(used only if param(5) == 0 & param(4) == 1)
param(7) = 0: absolute errors are computed
1: relative errors are computed
(used only if param(4) == 1)
param(8) = 0: not plot the solution
1: plot the solution (mesh)
2: plot the solution (surf)
3: plot the solution (contour)
param(9) = number of nodes in each element and along each direction
for plotting the solution
(used only if param(8) > 0)
Output: xyz = 3-indexes array wiht coordinates of 3D LGL mesh
un = numerical solution
D = (if param(3)==5) eigenvalues of A
param(21) = (if param(3)==1) cputime for reordering stiffness
matrix befor solving
param(22) = (if param(3)==1) cputime for solving linear system
by Choleski factorization
param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of
stiffness matrix A
param(25,26) = (if param(3)==4) extrema eigenvalues (max,min) of
interface Schur complement
param(27,28) = (if param(3)==6) extrema generalized eigenvalues
(max,min) of
A v = lambda M v, where A is the stiffness matrix and M is the
mass matrix
param(29,30,31) = (if param(4)==1 & (param(3)==1 | param(3)==3))
errors on the exact solution:
L^inf-norm, L2-norm, H1-norm
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 |
| | | | |
__________________________0001 function [xyz,un,D,param]=lap_3d(xa,xb,ya,yb,za,zb,gam,... 0002 uex,uex_x,uex_y,uex_z,ff,nex,nx,ney,ny,nez,nz,gammax,gammay,gammaz,... 0003 param); 0004 % LAP_3D Numerical solution of the 3D b.v.p. -Delta u + gam u = f 0005 % 0006 % -Delta u + gam u = f x in Omega 0007 % 0008 % + Dirichlet bc 0009 % 0010 % by SEM Numerical Integration with LGL quadrature formulas. 0011 % 0012 % [xy,un,D,param]=lap_3d(xa,xb,ya,yb,za,zb,gam,... 0013 % uex,uex_x,uex_y,uex_z,ff,nex,nx,ney,ny,nez,nz,gammax,gammay,gammaz,... 0014 % param); 0015 % 0016 % Input: xa= abscissa of either vertex V1 0017 % xb= abscissa of either vertex V2 0018 % ya= ordinate of either vertex V1 0019 % yb= ordinate of either vertex V3 0020 % za= ordinate of either vertex V1 0021 % zb= ordinate of either vertex V5 0022 % gam = coefficient of zeroth order term (constant>=0) 0023 % uex = exact solution (uex=@(x,y,z)[uex(x,y,z)], with .*, .^, ./) 0024 % uex_x = exact first x-derivative (uex_x=@(x,y,z)[uex_x(x,y,z)], 0025 % with .*, .^, ./) 0026 % uex_y = exact first y-derivative (uex_y=@(x,y,z)[uex_y(x,y,z)], 0027 % with .*, .^, ./) 0028 % uex_z = exact first z-derivative (uex_z=@(x,y,z)[uex_z(x,y,z)], 0029 % with .*, .^, ./) 0030 % ff = r.h.s. solution (ff=@(x,y,z)[ff(x,y,z)], with .*, .^, ./) 0031 % nex = number of elements (equally spaced) along x-direction 0032 % nx = polynomial degree in each element (the same in each element) 0033 % along x-direction 0034 % ney = number of elements (equally spaced) along y-direction 0035 % ny = polynomial degree in each element (the same in each element) 0036 % along y-direction 0037 % nez = number of elements (equally spaced) along z-direction 0038 % nz = polynomial degree in each element (the same in each element) 0039 % along z-direction 0040 % gammax = column or row array of length nex-1. 0041 % If the deomposition of Omega is not 0042 % uniform, gammax is the vector of position of interfaces between 0043 % spectral elements along x-direction. If gammax=[], uniform 0044 % decomposition is used. 0045 % gammay = column or row array of length ney-1. 0046 % If the deomposition of Omega is not 0047 % uniform, gammay is the vector of position of interfaces between 0048 % spectral elements along y-direction. If gammay=[], uniform 0049 % decomposition is used. 0050 % gammaz = column or row array of length nez-1. 0051 % If the deomposition of Omega is not 0052 % uniform, gammay is the vector of position of interfaces between 0053 % spectral elements along z-direction. If gammay=[], uniform 0054 % decomposition is used. 0055 % param = array of parameters 0056 % param(1) = 1 : SEM-NI approach, see Sect. 5.3.3, CHQZ3 0057 % param(2) = 0 : no reordering of the matrix, see pag. 194, CHQZ2 0058 % 1 : Cuthill-McKee ordering (by symrcm function of 0059 % Matlab) 0060 % 2 : Symmetric minimum degree ordering (by symamd function of 0061 % Matlab) 0062 % ONLY used if param(3)==1 0063 % param(3) = 1 : Solve the differential problem with a direct method 0064 % (Cholesky) and compute errors on the exact solution 0065 % 2 : Compute extrema eigenvalues of stiffness matrix 0066 % 3 : Compute interface Schur complement, solve the diff. 0067 % problem by Schur complement 0068 % 4: Compute Schur complement and its extrema eigenvalues 0069 % 5: Compute and plot all the eigenvalues of A (stiffness 0070 % matrix) 0071 % 6: Compute extrema eigenvalues of M^{-1}A, i.e. 0072 % the generalized eigenvalues of Av=lambda* M v, 0073 % where M is the SEM-NI mass matrix 0074 % param(4) = 1: compute errors (L^inf-norm, L2-norm, H1-norm) 0075 % on the exact solution 0076 % 2: no errors are computed 0077 % param(5) = 0: LG quadrature formulas with high precision degree are 0078 % used to compute norms (exact norms) 0079 % 1: LGL quadrature formulas with npdx,npdy nodes are 0080 % used to compute norms (discrete norms) 0081 % (used only if param(4) == 1) 0082 % param(6) = number of nodes for high degree quadrature formula, 0083 % (used only if param(5) == 0 & param(4) == 1) 0084 % param(7) = 0: absolute errors are computed 0085 % 1: relative errors are computed 0086 % (used only if param(4) == 1) 0087 % param(8) = 0: not plot the solution 0088 % 1: plot the solution (mesh) 0089 % 2: plot the solution (surf) 0090 % 3: plot the solution (contour) 0091 % param(9) = number of nodes in each element and along each direction 0092 % for plotting the solution 0093 % (used only if param(8) > 0) 0094 % 0095 % Output: xyz = 3-indexes array wiht coordinates of 3D LGL mesh 0096 % un = numerical solution 0097 % D = (if param(3)==5) eigenvalues of A 0098 % param(21) = (if param(3)==1) cputime for reordering stiffness 0099 % matrix befor solving 0100 % param(22) = (if param(3)==1) cputime for solving linear system 0101 % by Choleski factorization 0102 % param(23,24) = (if param(3)==2) extrema eigenvalues (max,min) of 0103 % stiffness matrix A 0104 % param(25,26) = (if param(3)==4) extrema eigenvalues (max,min) of 0105 % interface Schur complement 0106 % param(27,28) = (if param(3)==6) extrema generalized eigenvalues 0107 % (max,min) of 0108 % A v = lambda M v, where A is the stiffness matrix and M is the 0109 % mass matrix 0110 % param(29,30,31) = (if param(4)==1 & (param(3)==1 | param(3)==3)) 0111 % errors on the exact solution: 0112 % L^inf-norm, L2-norm, H1-norm 0113 % 0114 % 0115 % Omega=(xa,xb) x (ya,yb) x (za,zb) 0116 % 0117 % V8 _________ V7 0118 % /| /| 0119 % / | / | 0120 % V5 /________ /V6| 0121 % | | | | 0122 % |V4|______|__| V3 0123 % | / | / 0124 % | / |/ 0125 % |/________/ 0126 % V1 V2 0127 % 0128 % 0129 % 0130 % __________________________ 0131 % | | | | | 0132 % | 3 | 6 | 9 | 12 | Spectral elements 0133 % | | | | | ordering in a plane yz then 0134 % __________________________ planes at different x follow. 0135 % | | | | | 0136 % | 2 | 5 | 8 | 11 | 0137 % | | | | | 0138 % __________________________ 0139 % | | | | | 0140 % | 1 | 4 | 7 | 10 | 0141 % | | | | | 0142 % __________________________ 0143 % 0144 0145 xyz=[]; un=[]; D=[]; 0146 0147 npdx=nx+1; npdy=ny+1; npdz=nz+1; 0148 mn=npdx*npdy; ldnov=mn*npdz; ne=nex*ney*nez; 0149 0150 0151 [x,wx]=xwlgl(npdx); [dx]=derlgl(x,npdx); 0152 [y,wy]=xwlgl(npdy); [dy]=derlgl(y,npdy); 0153 [z,wz]=xwlgl(npdz); [dz]=derlgl(z,npdz); 0154 0155 % Mesh generation (with nov construction) 0156 0157 [xx,yy,zz,jacx,jacy,jacz,xyz,ww,ifro,nov]=mesh_3d(xa,xb,ya,yb,za,zb,nex,ney,nez,... 0158 npdx,npdy,npdz,x,wx,y,wy,z,wz,gammax,gammay,gammaz); 0159 0160 % Matrix assembling 0161 noe=length(ww); 0162 A=stiff_3d_se(npdx,nex,npdy,ney,npdz,nez,nov,wx,dx,jacx,wy,dy,jacy,wz,dz,jacz); 0163 M=spdiags(ww,0,noe,noe); 0164 0165 if gam ~=0 0166 A=A+gam*M; 0167 end 0168 0169 % Rigth Hand Side 0170 f=ff(xyz(:,1),xyz(:,2),xyz(:,3)).*ww; 0171 0172 0173 0174 % Lists of internal, boundary, interface nodes. All lists are referred to 0175 % global ordering on Omega 0176 % 0177 % lbor: list of boundary nodes 0178 % lint: list of internal nodes 0179 % lintint: list of internal nodes which are internal to spectral elements 0180 % lgamma: list of internal nodes which are on the interface between spectral elements 0181 0182 0183 0184 [ldir,lint,lintint,lgamma,ifro]=liste(ifro,nov); 0185 0186 0187 % Setting Dirichlet boundary conditions on both matrix and r.h.s 0188 ub=ifro.*uex(xyz(:,1),xyz(:,2),xyz(:,3)); 0189 0190 f(lint)=f(lint)-A(lint,ldir)*ub(ldir); 0191 A=A(lint,lint); M=M(lint,lint); 0192 f=f(lint); 0193 noei=length(lint); 0194 0195 % selection on param(3) 0196 0197 0198 if (param(3)==1) 0199 0200 % Solve the linear system by a direct method 0201 0202 if param(2)==0 0203 % no reordering 0204 t_ord=0; 0205 t=cputime; un=A\f;t_solv=cputime-t; 0206 0207 elseif param(2)==1 0208 % Cuthill-McKee ordering 0209 t=cputime; pp=symrcm(A);t_ord=cputime-t; 0210 t=cputime; un(pp)=A(pp,pp)\f(pp);t_solv=cputime-t; 0211 0212 elseif param(2)==2 0213 % Symmetric Minimum degree ordering 0214 t=cputime; pp=symamd(A);t_ord=cputime-t; 0215 t=cputime; un(pp)=A(pp,pp)\f(pp);t_solv=cputime-t; 0216 end 0217 0218 0219 param(21)=t_ord; 0220 param(22)=t_solv; 0221 0222 % 0223 ub(lint)=un; 0224 un=ub; 0225 0226 0227 0228 elseif(param(3)==2) 0229 0230 % Compute extrema eigenvalues of the stiffness matrix 0231 keigs=1; 0232 if(param(1)==1) 0233 % Symmetric problem 0234 opts.issym=1;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0235 lambda_max=eigs(A,keigs,'lm',opts); 0236 lambda_min=eigs(A,1,'sm',opts); 0237 else 0238 opts.issym=0;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0239 lambda_max=eigs(A,keigs,'lm',opts); 0240 lambda_min=eigs(A,1,'sm',opts); 0241 end 0242 param(23)=lambda_max; 0243 param(24)=lambda_min; 0244 0245 0246 elseif(param(3)==3) 0247 0248 % Compute interface Schur complement, solve the differential problem by 0249 % the Schur complement 0250 0251 Sigma=A(lgamma,lgamma)-A(lgamma,lintint)*inv(A(lintint,lintint))*A(lintint,lgamma); 0252 0253 % Solve Schur linear system and internal linear systems by \ of matlab 0254 0255 fsigma=f(lgamma)-A(lgamma,lintint)*inv(A(lintint,lintint))*f(lintint); 0256 xsigma=Sigma\fsigma; 0257 xint=A(lintint,lintint)\(f(lintint)-A(lintint,lgamma)*xsigma); 0258 u2=zeros(noei,1); 0259 u2(lgamma)=xsigma; 0260 u2(lintint)=xint; 0261 ub(lint)=u2; 0262 un=ub; 0263 0264 0265 elseif(param(3)==4) 0266 0267 % Compute interface Schur complement and its extrema eigenvalues 0268 0269 Sigma=A(lgamma,lgamma)-A(lgamma,lintint)*inv(A(lintint,lintint))*A(lintint,lgamma); 0270 0271 clear A; clear f; 0272 opts.issym=1;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0273 lambda_max=eigs(Sigma,1,'lm',opts); 0274 lambda_min=eigs(Sigma,1,'sm',opts); 0275 param(25)=lambda_max; 0276 param(26)=lambda_min; 0277 0278 0279 elseif(param(3)==5) 0280 % compute and plot eigenvalues of A 0281 D=eig(full(A)); 0282 fig=figure(..., 0283 'Name','Spectrum of SEM-NI stiffness matrix',... 0284 'Visible','on') 0285 plot(real(D),imag(D),'o') 0286 set(gca,'Fontname','Times','Fontsize',18) 0287 xlabel('Re'); ylabel('Im') 0288 if param(1)==1 0289 set(gca,'Fontname','Times','Fontsize',18) 0290 title(['N=',num2str(nx),' H=',num2str(2/nex), ' SEM-NI']) 0291 else 0292 set(gca,'Fontname','Times','Fontsize',18) 0293 title(['N=',num2str(nx),' H=',num2str(2/nex), ' Patching']) 0294 end 0295 0296 0297 elseif(param(3)==6) 0298 keigs=1; 0299 % compute extrema eigenvalues of M^{-1}A, i.e. 0300 % the generalized eigenvalues of Av=lambda* M v 0301 if(param(1)==1) 0302 opts.issym=1;opts.disp=0;opts.maxit=1000; opts.tol=1.d-12; 0303 lambda_max=eigs(A,M,keigs,'lm',opts); 0304 lambda_min=eigs(A,M,keigs,'sm',opts); 0305 param(27)=lambda_max; 0306 param(28)=lambda_min; 0307 0308 elseif(param(2)==2) 0309 disp('No eigenvalues are computed') 0310 end 0311 0312 end 0313 0314 if (param(4)==1 & (param(3)==1 | param(3)==3)) 0315 % compute errors on the exact solution 0316 [err_inf,err_h1,err_l2]=errors_3d(x,wx,dx,xx,jacx,y,wy,dy,yy,jacy,... 0317 z,wz,dz,zz,jacz,... 0318 xyz,ww,nov,un,uex,uex_x,uex_y,uex_z,param); 0319 param(29:31)=[err_inf;err_h1;err_l2]; 0320 end 0321 0322 if (param(8)>0 & (param(3)==1 | param(3)==3)) 0323 if param(8)==1 0324 command='mesh'; 0325 elseif param(8)==2 0326 command='surf'; 0327 elseif param(8)==3; 0328 command='contour' 0329 end 0330 % plot the solution 0331 fig=figure(..., 0332 'Name','SEM-NI solution of -Delta u + gam u=f, Dirichlet bc',... 0333 'Visible','on'); 0334 %[hf,ha]=plot_sem(fig,command,nex,ney,x,xx,jacx,y,yy,jacy,xy,ww,nov,un,param); 0335 0336 end