stencilops.h

著名的数学计算类库

// -*- C++ -*-
/***************************************************************************
 * blitz/array/stencilops.h  Stencil operators
 *
 * Copyright (C) 1997-2001 Todd Veldhuizen <[e-mail]>
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * Suggestions:          [e-mail]
 * Bugs:                 [e-mail]
 *
 * For more information, please see the Blitz++ Home Page:
 *    http://oonumerics.org/blitz/
 *
 ****************************************************************************/
#ifndef BZ_ARRAYSTENCILOPS_H
#define BZ_ARRAYSTENCILOPS_H

// NEEDS_WORK: need to factor many of the stencils in terms of the
// integer constants, e.g. 16*(A(-1,0)+A(0,-1)+A(0,1)+A(1,0))

#ifndef BZ_ARRAYSTENCILS_H
 #error <blitz/array/stencilops.h> must be included via <blitz/array/stencils.h>
#endif

#ifndef BZ_GEOMETRY_H
 #include <blitz/array/geometry.h>
#endif

#ifndef BZ_TINYMAT_H
 #include <blitz/tinymat.h>
#endif

BZ_NAMESPACE(blitz)

#define BZ_DECLARE_STENCIL_OPERATOR1(name,A)     \
  template<typename T>                              \
  inline _bz_typename T::T_numtype name(T& A)    \
  {

#define BZ_END_STENCIL_OPERATOR   }

#define BZ_DECLARE_STENCIL_OPERATOR2(name,A,B)       \
  template<typename T>                                  \
  inline _bz_typename T::T_numtype name(T& A, T& B)  \
  {

#define BZ_DECLARE_STENCIL_OPERATOR3(name,A,B,C) \
  template<typename T>                              \
  inline _bz_typename T::T_numtype name(T& A, T& B, T& C)    \
  {

// These constants are accurate to 45 decimal places = 149 bits of mantissa
const double recip_2 = .500000000000000000000000000000000000000000000;
const double recip_4 = .250000000000000000000000000000000000000000000;
const double recip_6 = .166666666666666666666666666666666666666666667;
const double recip_8 = .125000000000000000000000000000000000000000000;
const double recip_12 = .0833333333333333333333333333333333333333333333;
const double recip_144 = .00694444444444444444444444444444444444444444444;

/****************************************************************************
 * Laplacian Operators
 ****************************************************************************/

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D, A)
  return -4.0 * (*A)
    + A.shift(-1,0) + A.shift(1,0)
    + A.shift(-1,1) + A.shift(1,1);
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian3D, A)
  return -6.0 * (*A) 
    + A.shift(-1,0) + A.shift(1,0) 
    + A.shift(-1,1) + A.shift(1,1)
    + A.shift(-1,2) + A.shift(1,2);
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D4, A)
  return -60.0 * (*A) 
    + 16.0 * (A.shift(-1,0) + A.shift(1,0) + A.shift(-1,1) + A.shift(1,1))
    -        (A.shift(-2,0) + A.shift(2,0) + A.shift(-2,1) + A.shift(2,1));
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D4n, A)
  return Laplacian2D4(A) * recip_12;
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian3D4, A)
  return -90.0 * (*A) 
    + 16.0 * (A.shift(-1,0) + A.shift(1,0) + A.shift(-1,1) + A.shift(1,1) +
              A.shift(-1,2) + A.shift(1,2))
    -        (A.shift(-2,0) + A.shift(2,0) + A.shift(-2,1) + A.shift(2,1) +
              A.shift(-2,2) + A.shift(2,2));
BZ_END_STENCIL_OPERATOR

BZ_DECLARE_STENCIL_OPERATOR1(Laplacian3D4n, A)
  return Laplacian3D4(A) * recip_12;
BZ_END_STENCIL_OPERATOR

/****************************************************************************
 * Derivatives
 ****************************************************************************/

#define BZ_DECLARE_DIFF(name)  \
  template<typename T> \
  inline _bz_typename T::T_numtype name(T& A, int dim = firstDim)

#define BZ_DECLARE_MULTIDIFF(name) \
  template<typename T> \
  inline _bz_typename multicomponent_traits<_bz_typename     \
     T::T_numtype>::T_element name(T& A, int comp, int dim)

/****************************************************************************
 * Central differences with accuracy O(h^2)
 ****************************************************************************/

BZ_DECLARE_DIFF(central12) {
  return A.shift(1,dim) - A.shift(-1,dim);
}

BZ_DECLARE_DIFF(central22) {
  return -2.0 * (*A) + A.shift(1,dim) + A.shift(-1,dim);
}

BZ_DECLARE_DIFF(central32) {
  return -2.0 * (A.shift(1,dim) - A.shift(-1,dim))
    +           (A.shift(2,dim) - A.shift(-2,dim));
}

BZ_DECLARE_DIFF(central42) {
  return 6.0 * (*A)
    - 4.0 * (A.shift(1,dim) + A.shift(-1,dim))
    +       (A.shift(2,dim) + A.shift(-2,dim));
}

/****************************************************************************
 * Central differences with accuracy O(h^2)  (multicomponent versions)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(central12) {
  return A.shift(1,dim)[comp] - A.shift(-1,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(central22) {
  return -2.0 * (*A)[comp]
    + A.shift(1,dim)[comp] + A.shift(-1,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(central32) {
  return -2.0 * (A.shift(1,dim)[comp] - A.shift(-1,dim)[comp])
    +           (A.shift(2,dim)[comp] - A.shift(-2,dim)[comp]);
}

BZ_DECLARE_MULTIDIFF(central42) {
  return 6.0 * (*A)[comp]
    -4.0 * (A.shift(1,dim)[comp] + A.shift(-1,dim)[comp])
    +      (A.shift(2,dim)[comp] + A.shift(-2,dim)[comp]);
}

/****************************************************************************
 * Central differences with accuracy O(h^2)  (normalized versions)
 ****************************************************************************/

BZ_DECLARE_DIFF(central12n) {
  return central12(A,dim) * recip_2;
}

BZ_DECLARE_DIFF(central22n) {
  return central22(A,dim);
}

BZ_DECLARE_DIFF(central32n) {
  return central32(A,dim) * recip_2;
}

BZ_DECLARE_DIFF(central42n) {
  return central42(A,dim);
}

/****************************************************************************
 * Central differences with accuracy O(h^2)  (normalized multicomponent)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(central12n) {
  return central12(A,comp,dim) * recip_2;
}

BZ_DECLARE_MULTIDIFF(central22n) {
  return central22(A,comp,dim);
}

BZ_DECLARE_MULTIDIFF(central32n) {
  return central32(A,comp,dim) * recip_2;
}

BZ_DECLARE_MULTIDIFF(central42n) {
  return central42(A,comp,dim);
}

/****************************************************************************
 * Central differences with accuracy O(h^4)  
 ****************************************************************************/

BZ_DECLARE_DIFF(central14) {
  return 8.0 * (A.shift(1,dim) - A.shift(-1,dim))
    -          (A.shift(2,dim) - A.shift(-2,dim));
}

BZ_DECLARE_DIFF(central24) {
  return -30.0 * (*A)
    + 16.0 * (A.shift(1,dim) + A.shift(-1,dim))
    -        (A.shift(2,dim) + A.shift(-2,dim));
}

BZ_DECLARE_DIFF(central34) {
  return -13.0 * (A.shift(1,dim) - A.shift(-1,dim))
    +      8.0 * (A.shift(2,dim) - A.shift(-2,dim))
    -            (A.shift(3,dim) - A.shift(-3,dim));
}

BZ_DECLARE_DIFF(central44) {
  return 56.0 * (*A)
    - 39.0 * (A.shift(1,dim) + A.shift(-1,dim))
    + 12.0 * (A.shift(2,dim) + A.shift(-2,dim))
    -        (A.shift(3,dim) + A.shift(-3,dim));
}

/****************************************************************************
 * Central differences with accuracy O(h^4)  (multicomponent versions)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(central14) {
  return 8.0 * (A.shift(1,dim)[comp] - A.shift(-1,dim)[comp])
    -          (A.shift(2,dim)[comp] - A.shift(-2,dim)[comp]);
}

BZ_DECLARE_MULTIDIFF(central24) {
  return - 30.0 * (*A)[comp]
    + 16.0 * (A.shift(1,dim)[comp] + A.shift(-1,dim)[comp])
    -        (A.shift(2,dim)[comp] + A.shift(-2,dim)[comp]);
}

BZ_DECLARE_MULTIDIFF(central34) {
  return -13.0 * (A.shift(1,dim)[comp] - A.shift(-1,dim)[comp])
    +      8.0 * (A.shift(2,dim)[comp] - A.shift(-2,dim)[comp])
    -            (A.shift(3,dim)[comp] - A.shift(-3,dim)[comp]);
}

BZ_DECLARE_MULTIDIFF(central44) {
  return 56.0 * (*A)[comp]
    - 39.0 * (A.shift(1,dim)[comp] + A.shift(-1,dim)[comp])
    + 12.0 * (A.shift(2,dim)[comp] + A.shift(-2,dim)[comp])
    -        (A.shift(3,dim)[comp] + A.shift(-3,dim)[comp]);
}

/****************************************************************************
 * Central differences with accuracy O(h^4)  (normalized)
 ****************************************************************************/

BZ_DECLARE_DIFF(central14n) {
  return central14(A,dim) * recip_12;
}

BZ_DECLARE_DIFF(central24n) {
  return central24(A,dim) * recip_12;
}

BZ_DECLARE_DIFF(central34n) {
  return central34(A,dim) * recip_8;
}

BZ_DECLARE_DIFF(central44n) {
  return central44(A,dim) * recip_6;
}

/****************************************************************************
 * Central differences with accuracy O(h^4)  (normalized, multicomponent)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(central14n) {
  return central14(A,comp,dim) * recip_12;
}

BZ_DECLARE_MULTIDIFF(central24n) {
  return central24(A,comp,dim) * recip_12;
}

BZ_DECLARE_MULTIDIFF(central34n) {
  return central34(A,comp,dim) * recip_8;
}

BZ_DECLARE_MULTIDIFF(central44n) {
  return central44(A,comp,dim) * recip_6;
}

/****************************************************************************
 * Backward differences with accuracy O(h)
 ****************************************************************************/

BZ_DECLARE_DIFF(backward11) {
  return (*A) - A.shift(-1,dim);
}

BZ_DECLARE_DIFF(backward21) {
  return (*A) - 2.0 * A.shift(-1,dim) + A.shift(-2,dim);
}

BZ_DECLARE_DIFF(backward31) {
  return (*A) - 3.0 * A.shift(-1,dim) + 3.0 * A.shift(-2,dim)
    - A.shift(-3,dim);
}

BZ_DECLARE_DIFF(backward41) {
  return (*A) - 4.0 * A.shift(-1,dim) + 6.0 * A.shift(-2,dim)
    - 4.0 * A.shift(-3,dim) + A.shift(-4,dim);
}

/****************************************************************************
 * Backward differences with accuracy O(h) (multicomponent versions)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(backward11) {
  return (*A)[comp] - A.shift(-1,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward21) {
  return (*A)[comp] - 2.0 * A.shift(-1,dim)[comp] + A.shift(-2,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward31) {
  return (*A)[comp] - 3.0 * A.shift(-1,dim)[comp] + 3.0 * A.shift(-2,dim)[comp]
    - A.shift(-3,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward41) {
  return (*A)[comp] - 4.0 * A.shift(-1,dim)[comp] + 6.0 * A.shift(-2,dim)[comp]
    - 4.0 * A.shift(-3,dim)[comp] + A.shift(-4,dim)[comp];
}

/****************************************************************************
 * Backward differences with accuracy O(h)  (normalized)
 ****************************************************************************/

BZ_DECLARE_DIFF(backward11n) { return backward11(A,dim); }
BZ_DECLARE_DIFF(backward21n) { return backward21(A,dim); }
BZ_DECLARE_DIFF(backward31n) { return backward31(A,dim); }
BZ_DECLARE_DIFF(backward41n) { return backward41(A,dim); }

/****************************************************************************
 * Backward differences with accuracy O(h)  (normalized, multicomponent)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(backward11n) { return backward11(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(backward21n) { return backward21(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(backward31n) { return backward31(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(backward41n) { return backward41(A,comp,dim); }

/****************************************************************************
 * Backward differences with accuracy O(h^2)
 ****************************************************************************/

BZ_DECLARE_DIFF(backward12) {
  return 3.0 * (*A) - 4.0 * A.shift(-1,dim) + A.shift(-2,dim);
}

BZ_DECLARE_DIFF(backward22) {
  return 2.0 * (*A) - 5.0 * A.shift(-1,dim) + 4.0 * A.shift(-2,dim)
    - A.shift(-3,dim);
}

BZ_DECLARE_DIFF(backward32) {
  return 5.0 * (*A) - 18.0 * A.shift(-1,dim) + 24.0 * A.shift(-2,dim)
    - 14.0 * A.shift(-3,dim) + 3.0 * A.shift(-4,dim);
}

BZ_DECLARE_DIFF(backward42) {
  return 3.0 * (*A) - 14.0 * A.shift(-1,dim) + 26.0 * A.shift(-2,dim)
    - 24.0 * A.shift(-3,dim) + 11.0 * A.shift(-4,dim) - 2.0 * A.shift(-5,dim);
}

/****************************************************************************
 * Backward differences with accuracy O(h^2) (multicomponent versions)
 ****************************************************************************/

BZ_DECLARE_MULTIDIFF(backward12) {
  return 3.0 * (*A)[comp] - 4.0 * A.shift(-1,dim)[comp]
    + A.shift(-2,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward22) {