stencilops.h

数值计算工具库,C语言编写的,可以直接调用.

#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_ARRAYSTENCIL_H
 #error <blitz/array/stencilops.h> must be included via <blitz/array/stencil.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<class T>                              \
  inline _bz_typename T::T_numtype name(T& A)    \
  {

#define BZ_END_STENCIL_OPERATOR   }

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

#define BZ_DECLARE_STENCIL_OPERATOR3(name,A,B,C) \
  template<class 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. * A 
      + 16.*(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. * A 
      + 16.*(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<class T> \
  inline _bz_typename T::T_numtype name(T& A, int dim = firstDim)

#define BZ_DECLARE_MULTIDIFF(name) \
  template<class 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 A.shift(-1,dim) - 2. * A + A.shift(+1,dim);
}

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

BZ_DECLARE_DIFF(central42) {
  return A.shift(-2,dim) + A.shift(2,dim) -4.*(A.shift(-1,dim)+A.shift(+1,dim)) 
    +6.*A.shift(0,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 A.shift(-1,dim)[comp] - 2. * (*A)[comp] + A.shift(+1,dim)[comp];
}

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

BZ_DECLARE_MULTIDIFF(central42) {
  return A.shift(-2,dim)[comp] -4.*A.shift(-1,dim)[comp] 
    +6.*A.shift(0,dim)[comp] -4.*A.shift(1,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 (A.shift(-2,dim) - A.shift(2,dim)) 
     + 8.*(A.shift(1,dim)-A.shift(-1,dim));
}

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

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

BZ_DECLARE_DIFF(central44) {
  return -1.*A.shift(-3,dim)+12.*A.shift(-2,dim)-39.*A.shift(-1,dim)
    +56.*A-39.*A.shift(1,dim)+12.*A.shift(2,dim)-A.shift(3,dim);
}

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

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

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

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

BZ_DECLARE_MULTIDIFF(central44) {
  return -1.*A.shift(-3,dim)[comp]+12.*A.shift(-2,dim)[comp]
    -39.*A.shift(-1,dim)[comp] +56.*(*A)[comp]-39.*A.shift(1,dim)[comp]
    +12.*A.shift(2,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.*A.shift(-1,dim) + A.shift(-2,dim);
}

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

BZ_DECLARE_DIFF(backward41) {
  return A - 4.*A.shift(-1,dim) + 6.*A.shift(-2,dim) -4.*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.*A.shift(-1,dim)[comp] + A.shift(-2,dim)[comp];
}

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

BZ_DECLARE_MULTIDIFF(backward41) {
  return (*A)[comp] - 4.*A.shift(-1,dim)[comp] + 6.*A.shift(-2,dim)[comp]
    -4.*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.*A -4.*A.shift(-1,dim) + A.shift(-2,dim);
}

BZ_DECLARE_DIFF(backward22) {
  return 2.*A -5.*A.shift(-1,dim) + 4.*A.shift(-2,dim) -A.shift(-3,dim);
}

BZ_DECLARE_DIFF(backward32) {
  return 5.*A - 18.*A.shift(-1,dim) + 24.*A.shift(-2,dim) -14.*A.shift(-3,dim) 
    + 3.*A.shift(-4,dim);
}

BZ_DECLARE_DIFF(backward42) {
  return 3.*A -14.*A.shift(-1,dim) + 26.*A.shift(-2,dim) -24.*A.shift(-3,dim) 
    + 11.*A.shift(-4,dim) -2.*A.shift(-5,dim);
}

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

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

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

BZ_DECLARE_MULTIDIFF(backward32) {
  return 5.*(*A)[comp] - 18.*A.shift(-1,dim)[comp] + 24.*A.shift(-2,dim)[comp]
     -14.*A.shift(-3,dim)[comp] + 3.*A.shift(-4,dim)[comp];
}

BZ_DECLARE_MULTIDIFF(backward42) {
  return 3.*(*A)[comp] -14.*A.shift(-1,dim)[comp] + 26.*A.shift(-2,dim)[comp]
    -24.*A.shift(-3,dim)[comp] + 11.*A.shift(-4,dim)[comp] 
    -2.*A.shift(-5,dim)[comp];
}

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

BZ_DECLARE_DIFF(backward12n) { return backward12(A,dim) * recip_2; }
BZ_DECLARE_DIFF(backward22n) { return backward22(A,dim); }