matrices_1d

 MATRICES_1D   Assembles SEM stiffness and mass matrices for 1D b.v.p.

     M_{ij}=(phi_j,phi_i)_N
     A_{ij}=(grad phi_j, grad phi_i)_N
    
 SEM  Numerical Integration with LGL quadrature formulas.

  [A,M]=matrices_1d(xa,xb,cb,ne,nx);

 Input: xa, xb = extrema of computational domain Omega=(xa,xb)
      cb = string containing boundary conditions, for example
           cb='dn' imposes Dirichlet in xa and Neumann in xb 
      ne = number of elements (equally spaced)
      nx = polynomial degree in each element (the same in each element)
             to be set only if p=4, otherwise, nx=p;

 Output: A = stiffness matrix
         M = Mass matrix

 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,M]=matrices_1d(xa,xb,cb,ne,nx);
0002 % MATRICES_1D   Assembles SEM stiffness and mass matrices for 1D b.v.p.
0003 %
0004 %     M_{ij}=(phi_j,phi_i)_N
0005 %     A_{ij}=(grad phi_j, grad phi_i)_N
0006 %
0007 % SEM  Numerical Integration with LGL quadrature formulas.
0008 %
0009 %  [A,M]=matrices_1d(xa,xb,cb,ne,nx);
0010 %
0011 % Input: xa, xb = extrema of computational domain Omega=(xa,xb)
0012 %      cb = string containing boundary conditions, for example
0013 %           cb='dn' imposes Dirichlet in xa and Neumann in xb
0014 %      ne = number of elements (equally spaced)
0015 %      nx = polynomial degree in each element (the same in each element)
0016 %             to be set only if p=4, otherwise, nx=p;
0017 %
0018 % Output: A = stiffness matrix
0019 %         M = Mass matrix
0020 %
0021 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0022 %                    "Spectral Methods. Fundamentals in Single Domains"
0023 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0024 
0025 %   Written by Paola Gervasio
0026 %   $Date: 2007/04/01$
0027 
0028 
0029 %
0030 npdx=nx+1; 
0031 
0032 [x,wx]=xwlgl(npdx); dx=derlgl(x,npdx);
0033 
0034 % nov
0035 ldnov=npdx; nov=zeros(ldnov,ne);
0036 [nov]=cosnov1d(npdx,ne,nov);
0037 noe=nov(npdx,ne);
0038 
0039 % Uniform decomposition of  Omega in ne elements
0040 
0041 [xx,jacx,xy,ww]=mesh1d(xa,xb,ne,npdx,nov,x,wx);
0042 
0043 % ww contains the diagonal of the mass matrix
0044 
0045 % Stiffness Matrix assembling
0046 A=stiff_1d_se(npdx,ne,nov,wx,dx,jacx); 
0047 M=spdiags(ww,0,noe,noe);
0048 
0049 % Setting boundary conditions on the matrix
0050 
0051 if cb(1)=='d' & cb(2)=='d'; 
0052 lint=(2:noe-1); 
0053 elseif cb(1)=='n' & cb(2)=='d';
0054 lint=(1:noe-1);
0055 elseif cb(1)=='d' & cb(2)=='n';
0056 lint=(2:noe);
0057 elseif cb(1)=='n' & cb(2)=='n';
0058 lint=(1:noe);
0059 end
0060 
0061 A=A(lint,lint); M=M(lint,lint);
0062 
0063 return