Home > atmlab > gridcreation > gridselect3D.m

gridselect3D

PURPOSE ^

GRIDSELECT3D Simple selection of 3D grids

SYNOPSIS ^

function [xc,yc,zc,Ac] = gridselect3D( xf, yf, zf, Af, tol, repr, abs_or_rel )

DESCRIPTION ^

 GRIDSELECT3D   Simple selection of 3D grids

    The function needs a set of realisations of the quantity for which the
    grids shall be selected. These realisations are given along dim. 4 of the
    4D variable *Af*, are provided on a sufficient fine grid that 
    representation errors can be neglected and shall span the possible cases 
    that shall be represented. 

    The dimension of Af shall be [ xf, yf, zf, number of realisation ]

    The grid is selected by first including the end points of *xf/yf/zf*. 
    The coarse representation is transformed to the fine one, and a new grid 
    point is included where the maximum deviation is found (where all
    realisations are considered in parallel). The point is included in
    either *xc*, *yc* or *zc*, depending on which that most effectively 
    reduces the deviation. This procedure is repeated until the coarse and 
    representations deviates less then *tol*.

    See also GRIDSELECT1D and GRIDSELECT2D.

 FORMAT   [xc,yc,zc,Ac] = gridselect3D( xf, yf, zf, Af, tol, repr, abs_or_rel)

 OUT   xc          Obtained coarse grid for x dimension.
       yc          Obtained coarse grid for y dimension.
       zc          Obtained coarse grid for z dimension.
       Ac          Af at grid points xc, yc and zc.
 IN    xf          Reference fine grid for x dimension.
       yf          Reference fine grid for y dimension.
       zf          Reference fine grid for y dimension.
       Af          Reference data (a matrix) on fine grid.
       tol         Acceptable tolarance.
       repr        String describing representation used. Possible choices
                   are equal to allowed 'methods' for *interp1*.
 OPT   abs_or_rel  String telling if tolarance is in absolute ('abs') or
                   relative ('rel') units. Default is 'abs'.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

gridselect3D.m

SOURCE CODE ^

0001 % GRIDSELECT3D   Simple selection of 3D grids
0002 %
0003 %    The function needs a set of realisations of the quantity for which the
0004 %    grids shall be selected. These realisations are given along dim. 4 of the
0005 %    4D variable *Af*, are provided on a sufficient fine grid that
0006 %    representation errors can be neglected and shall span the possible cases
0007 %    that shall be represented.
0008 %
0009 %    The dimension of Af shall be [ xf, yf, zf, number of realisation ]
0010 %
0011 %    The grid is selected by first including the end points of *xf/yf/zf*.
0012 %    The coarse representation is transformed to the fine one, and a new grid
0013 %    point is included where the maximum deviation is found (where all
0014 %    realisations are considered in parallel). The point is included in
0015 %    either *xc*, *yc* or *zc*, depending on which that most effectively
0016 %    reduces the deviation. This procedure is repeated until the coarse and
0017 %    representations deviates less then *tol*.
0018 %
0019 %    See also GRIDSELECT1D and GRIDSELECT2D.
0020 %
0021 % FORMAT   [xc,yc,zc,Ac] = gridselect3D( xf, yf, zf, Af, tol, repr, abs_or_rel)
0022 %
0023 % OUT   xc          Obtained coarse grid for x dimension.
0024 %       yc          Obtained coarse grid for y dimension.
0025 %       zc          Obtained coarse grid for z dimension.
0026 %       Ac          Af at grid points xc, yc and zc.
0027 % IN    xf          Reference fine grid for x dimension.
0028 %       yf          Reference fine grid for y dimension.
0029 %       zf          Reference fine grid for y dimension.
0030 %       Af          Reference data (a matrix) on fine grid.
0031 %       tol         Acceptable tolarance.
0032 %       repr        String describing representation used. Possible choices
0033 %                   are equal to allowed 'methods' for *interp1*.
0034 % OPT   abs_or_rel  String telling if tolarance is in absolute ('abs') or
0035 %                   relative ('rel') units. Default is 'abs'.
0036 
0037 % 2004-09-23   Created by Patrick Eriksson.
0038 
0039 
0040 function [xc,yc,zc,Ac] = gridselect3D( xf, yf, zf, Af, tol, repr, abs_or_rel )
0041 
0042 if nargin < 7
0043   abs_or_rel = 'abs';
0044 end
0045 
0046 
0047 nx = length( xf );
0048 ny = length( yf );
0049 nz = length( zf );
0050 
0051 
0052 %=== Check input
0053 %
0054 if ~isvector( xf )
0055   error( 'Input argument *xf* must be a vector.' );
0056 end
0057 %
0058 if ~isvector( yf )
0059   error( 'Input argument *yf* must be a vector.' );
0060 end
0061 %
0062 if ~isvector( yf )
0063   error( 'Input argument *zf* must be a vector.' );
0064 end
0065 %
0066 if ~( ndims( Af ) == 3  |  ndims( Af ) == 4 )  |  ~isnumeric( Af )
0067   error( 'Input argument *Af* must be a 4D matrix.' );
0068 end
0069 %
0070 if size( Af, 1 ) ~= nx  &&  size( Af, 2 ) ~= ny  &&  size( Af, 3 ) ~= nz
0071   error( 'Sizes of *xf*/*yf* and *Af* do not match.' );
0072 end
0073 %
0074 if strcmp( abs_or_rel, 'rel' )  &  any( Af == 0 )
0075    error( 'With *abs_or_rel* = ''rel'', *Af can not contain zeros.' );
0076 end
0077 
0078 xind = [1 nx];
0079 yind = [1 ny];
0080 zind = [1 nz];
0081 
0082 xf = vec2row( xf );
0083 yf = vec2row( yf );
0084 zf = vec2row( zf );
0085 
0086 Ac   = zeros( size( Af ) );
0087 
0088 Ac = do_interp( xf(xind), yf(yind), zf(zind), Af(xind,yind,zind,:), ...
0089                                                             xf, yf, zf, repr );
0090 [e,xip,yip,zip] = get_max_dev( Af, Ac, abs_or_rel );
0091 
0092 while e > tol
0093 
0094   %= Test to put in xip
0095   if ~any( xip == xind )
0096     xindtmp = sort( [ xind xip ] );
0097     Ac = do_interp( xf(xindtmp), yf(yind), zf(zind), ...
0098                                    Af(xindtmp,yind,zind,:), xf, yf, zf, repr );
0099     [e1,xip1,yip1,zip1] = get_max_dev( Af, Ac, abs_or_rel );
0100   else
0101     e1 = Inf;
0102   end
0103 
0104   %= Test to put in yip
0105   if ~any( yip == yind )
0106     yindtmp = sort( [ yind yip ] );
0107     Ac = do_interp( xf(xind), yf(yindtmp), zf(zind), ...
0108                                    Af(xind,yindtmp,zind,:), xf, yf, zf, repr );
0109     [e2,xip2,yip2,zip2] = get_max_dev( Af, Ac, abs_or_rel );
0110   else
0111     e2 = Inf;
0112   end
0113 
0114   %= Test to put in zip
0115   if ~any( zip == zind )
0116     zindtmp = sort( [ zind zip ] );
0117     Ac = do_interp( xf(xind), yf(yind), zf(zindtmp), ...
0118                                    Af(xind,yind,zindtmp,:), xf, yf, zf, repr );
0119     [e3,xip3,yip3,zip3] = get_max_dev( Af, Ac, abs_or_rel );
0120   else
0121     e3 = Inf;
0122   end
0123   
0124   if e1 <= e2  &  e1 <= e3
0125     xind = xindtmp;
0126     xip  = xip1;
0127     yip  = yip1;
0128     zip  = zip1;
0129     e    = e1;
0130   elseif e2 <= e3
0131     yind = yindtmp;
0132     xip  = xip2;
0133     yip  = yip2;
0134     zip  = zip2;
0135     e    = e2;
0136   else
0137     zind = zindtmp;
0138     xip  = xip3;
0139     yip  = yip3;
0140     zip  = zip3;
0141     e    = e3;
0142   end
0143 end
0144 
0145 xc = xf(xind);
0146 yc = yf(yind);
0147 zc = zf(zind);
0148 
0149 if nargout > 3
0150   Ac = Af(xind,yind,zind,:);
0151 end
0152 
0153 
0154 return
0155 
0156 
0157 function Ac = do_interp( xc, yc, zc, A, xf, yf, zf, repr )
0158   %
0159   for i = 1 : size(A,4)
0160     Ac(:,:,:,i) = interp3( yc, xc, zc, A(:,:,:,i), yf, xf, zf', repr );
0161   end
0162   %
0163 return
0164 
0165 
0166 function [e,xip,yip,zip] = get_max_dev( Af, Ac, abs_or_rel );
0167   %
0168   if strcmp( abs_or_rel, 'abs' )
0169     [e,xip] = max( max( max( max( abs( Ac - Af ), [], 4 ), [], 3 ), [], 2 ) );
0170     [e,yip] = max( max( max( max( abs( Ac - Af ), [], 4 ), [], 3 ), [], 1 ) );
0171     [e,zip] = max( max( max( max( abs( Ac - Af ), [], 4 ), [], 2 ), [], 1 ) );
0172   elseif strcmp( abs_or_rel, 'rel' )
0173     [e,xip] = max( max( max( max( abs( (Ac - Af)./Af ), [], 4 ), [], 3 ), [], 2 ) );
0174     [e,yip] = max( max( max( max( abs( (Ac - Af)./Af ), [], 4 ), [], 3 ), [], 1 ) );
0175     [e,zip] = max( max( max( max( abs( (Ac - Af)./Af ), [], 4 ), [], 2 ), [], 1 ) );
0176   else
0177     error( ...
0178       sprintf('Unknown error type (%s). Can only be ''abs'' or ''rel''.', ...
0179                                                                 abs_or_rel ) );
0180   end
0181   %
0182 return

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005