STIFF_3D_SP Computes 3D local stiffness SEM matrix associated to (nabla(phi_j), nabla(phi_i))_N
(nabla(phi_j), nabla(phi_i))_N
[A]=stiff_3d_sp(wx,dx,jacx,wy,dy,jacy,wz,dz,jacz);
produces the matrix
A_{ij}=(nabla(phi_j), nabla(phi_i))_N
of size
(mn,mn) where mn=npdx*npdy is the local number of d.o.f.
Input :
wx = npdx LGL weigths in [-1,1],
(produced by calling [x,wx]=xwlgl(npdx))
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
wz = npdz LGL weigths in [-1,1],
(produced by calling [z,wz]=xwlgl(npdz))
dz =first derivative LGL matrix (by calling dzderlgl(z,npdz))
jacz= array (length(jacy)=ne); jacz(ie)= (y_V5_ie-y_V1_ie)/2
ww = column array with local weigths, length: npdx*npdy*npdz
Output: A = matrix (npdx*npdy*npdz,npdx*npdy*npdz)
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 [A]=stiff_3d_sp(wx,dx,jacx,wy,dy,jacy,wz,dz,jacz); 0002 % STIFF_3D_SP Computes 3D local stiffness SEM matrix associated to (nabla(phi_j), nabla(phi_i))_N 0003 % 0004 % (nabla(phi_j), nabla(phi_i))_N 0005 % 0006 % 0007 % [A]=stiff_3d_sp(wx,dx,jacx,wy,dy,jacy,wz,dz,jacz); 0008 % produces the matrix 0009 % A_{ij}=(nabla(phi_j), nabla(phi_i))_N 0010 % of size 0011 % (mn,mn) where mn=npdx*npdy is the local number of d.o.f. 0012 % 0013 % Input : 0014 % wx = npdx LGL weigths in [-1,1], 0015 % (produced by calling [x,wx]=xwlgl(npdx)) 0016 % dx =first derivative LGL matrix (by calling dx=derlgl(x,npdx)) 0017 % jacx = array (length(jacx)=ne); jacx(ie)= (x_V2_ie-x_V1_ie)/2 0018 % wy = npdy LGL weigths in [-1,1], 0019 % (produced by calling [y,wy]=xwlgl(npdy)) 0020 % dy =first derivative LGL matrix (by calling dy=derlgl(y,npdy)) 0021 % jacy = array (length(jacy)=ne); jacy(ie)= (y_V3_ie-y_V1_ie)/2 0022 % wz = npdz LGL weigths in [-1,1], 0023 % (produced by calling [z,wz]=xwlgl(npdz)) 0024 % dz =first derivative LGL matrix (by calling dzderlgl(z,npdz)) 0025 % jacz= array (length(jacy)=ne); jacz(ie)= (y_V5_ie-y_V1_ie)/2 0026 % ww = column array with local weigths, length: npdx*npdy*npdz 0027 % 0028 % Output: A = matrix (npdx*npdy*npdz,npdx*npdy*npdz) 0029 % 0030 % 0031 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0032 % "Spectral Methods. Fundamentals in Single Domains" 0033 % Springer Verlag, Berlin Heidelberg New York, 2006. 0034 0035 % Written by Paola Gervasio 0036 % $Date: 2007/04/01$ 0037 0038 0039 npdx=length(wx); 0040 npdy=length(wy); 0041 npdz=length(wz); 0042 jac=jacx*jacy*jacz; 0043 jacyx=jacy/jacx; 0044 jacxy=jacx/jacy; 0045 mn=npdx*npdy; 0046 ldnov=mn*npdz; 0047 A=sparse(ldnov,ldnov); 0048 0049 % * dx' * diag(wx) * dx 0050 sx = dx' * diag(wx) * dx; 0051 0052 % * dy' * diag(wy) * dy 0053 sy = dy' * diag(wy) * dy; 0054 0055 % * dz' * diag(wz) * dz 0056 sz = dz' * diag(wz) * dz; 0057 0058 for li=1:npdz; 0059 il=(li-1)*mn; 0060 0061 for ki=1:npdy; 0062 ik=(ki-1)*npdx; 0063 coefx=wy(ki)*wz(li)*jacz*jacy/jacx; 0064 0065 for i=1:npdx; 0066 ii=il+ik+i; 0067 coefy=wx(i)*wz(li)*jacz*jacx/jacy; 0068 coefz=wx(i)*wy(ki)*jacy*jacx/jacz; 0069 0070 i1=il+ik+1; 0071 i2=il+ik+npdx; 0072 A(ii,i1:i2)=A(ii,i1:i2)+sx(i,1:npdx)*coefx; 0073 0074 i1=il+i;i2=il+mn; 0075 A(ii,i1:npdx:i2)=A(ii,i1:npdx:i2)+sy(ki,1:npdy)*coefy; 0076 0077 i1=ik+i; 0078 A(ii,i1:mn:ldnov)=A(ii,i1:mn:ldnov)+sz(li,1:npdz)*coefz; 0079 end 0080 end 0081 end 0082 0083 0084 % A=sparse(ldnov,mn);A=0; 0085 % B=sparse(mn,mn);B=0; 0086 % for ki=1:npdy 0087 % inde=((ki-1)*npdx+1:ki*npdx); 0088 % A(inde,inde)=dx'*(diag(ww(inde))*dx)*jacyx; 0089 % end 0090 % 0091 % for i=1:npdx 0092 % inde=(i:npdx:mn); 0093 % B(inde,inde)=dy'*(diag(ww(inde))*dy)*jacxy; 0094 % end 0095 % A=A+B; 0096 %