quadgridold.m

MDPSAS工具箱是马里兰大学开发的

function P = quadgrid(varargin)

% quadgrid.m Quadrature grid object constructor.
%            Computes the quadrature points and weights, the differentiation
%            arrays for the first-order differentiation and Laplacian operators,
%            and computes the discrete transorm array corresponding to the
%            appropriate Jacobi polynomial sequence. Call as
%
%      P = quadgrid(geom,ndisc,name,axlim)
%
%          geom  : geometry - 'slab'(rect),'disk'(cyln),'sphe','peri'
%          ndisc : number of discretization points, including endpoints
%          name  : coordinate name, e.g., 'x', 'y','z', 't', 'r', etc.
%          axlim : axis limits in true dimensions
%
% or directly: P = quadgrid(geom,qp,w,Q,name,d,dd,d2,dd2)
%
%      Object fields:
%        P.geom  : Symmetry of coordinate system ('slab','cyln')
%        P.qp    : (ndisc by 1) set of discretization points in
%                  the specified;
%        P.w     : (ndisc by 1) quadrature weight array; note that
%                  the quadrature array contains the factor z.^a;
%        P.Q     : discretized Jacobi polynomials in column format;
%        P.Qinv  : inverse of Q used for interpolation;
%        P.name  : coordinate name;
%        P.d     : (ndisc by ndisc) first order derivative array;
%        P.dd    : (ndisc by ndisc) Laplacian operator array;
%        P.d2    : (ndisc by ndisc) (1/z)(d/dz)z
%        P.dd2   : (ndisc by ndisc) (d/dz)(1/z)(d/dz)z
%
% Use the * operator for quadgrid objects to create higher-dimensional
% quadrature grids.

% Written by Raymond A. Adomaitis as part of MWRtools, version 2.
% Copyright (c) by Raymond A. Adomaitis, 1998-2002.

if nargin == 0
   P = makegrid(' ',[],[],[],[],' ',[],[],[],[]);
   return
   end

if length(varargin) <= 10 & ~isa(varargin{3},'char') % The latter condition fails when varargin{3}=name

% Construct the object with fields specified directly.

   P = makegrid(varargin{1},varargin{2},varargin{3},varargin{4},varargin{5}, ...
                varargin{6},varargin{7},varargin{8},varargin{9});
   return
   end

% Compute quadrature gird points, weights, etc.

geom  = varargin{1};
ndisc = varargin{2};
name  = varargin{3};

if length(varargin) == 4
   axlim = varargin{end};
else
   axlim = [0 1];
   if strcmp(geom,'peri')
      axlim = [0 2*pi];
      end
   end

[z,w,Q,d,dd,d2,dd2] = pd(geom,ndisc,axlim);

P = makegrid(geom,z,w,Q,name,d,dd,d2,dd2);

return

% ------------------------------------------------------ %

function P = makegrid(geom,qp,w,Q,name,d,dd,d2,dd2)

if isa(geom,'cell')
   P.geom = geom;
else
   P.geom = {geom};
   end

if isa(qp,'cell')
   P.qp = qp;
else
   P.qp = {qp};
   end

if isa(w,'cell')
   P.w = w;
else
   P.w = {w};
   end

if isa(Q,'cell')
   P.Q = Q;
   for i = 1:length(Q)
      P.Qinv{i} = inv(Q{i});
      end
else
   P.Q = {Q};
   P.Qinv = {inv(Q)};
   end

if isa(name,'cell')
   P.name = name;
else
   P.name = {name};
   end

if isa(d,'cell')
   P.d = d;
else
   P.d = {d};
   end

if isa(dd,'cell')
   P.dd = dd;
else
   P.dd = {dd};
   end

if isa(d2,'cell')
   P.d2 = d2;
else
   P.d2 = {d2};
   end

if isa(dd2,'cell')
   P.dd2 = dd2;
else
   P.dd2 = {dd2};
   end

P = class(P,'quadgrid');