stencilops.h

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

BZ_DECLARE_DIFF(backward32n) { return backward32(A,dim) * recip_2; }
BZ_DECLARE_DIFF(backward42n) { return backward42(A,dim); }

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

BZ_DECLARE_MULTIDIFF(backward12n) { return backward12(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(backward22n) { return backward22(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(backward32n) { return backward32(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(backward42n) { return backward42(A,comp,dim); }

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

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

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

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

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

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

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

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

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

BZ_DECLARE_MULTIDIFF(forward41) {
  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];
}

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

BZ_DECLARE_DIFF(forward11n) { return forward11(A,dim); }
BZ_DECLARE_DIFF(forward21n) { return forward21(A,dim); }
BZ_DECLARE_DIFF(forward31n) { return forward31(A,dim); }
BZ_DECLARE_DIFF(forward41n) { return forward41(A,dim); }

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

BZ_DECLARE_MULTIDIFF(forward11n) { return forward11(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward21n) { return forward21(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward31n) { return forward31(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward41n) { return forward41(A,comp,dim); }

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

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

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

BZ_DECLARE_DIFF(forward32) {
  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(forward42) {
  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);
}

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

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

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

BZ_DECLARE_MULTIDIFF(forward32) {
  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(forward42) {
  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] 
    + 11.*A.shift(5,dim)[comp];
}


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

BZ_DECLARE_DIFF(forward12n) { return forward12(A,dim) * recip_2; }
BZ_DECLARE_DIFF(forward22n) { return forward22(A,dim); }
BZ_DECLARE_DIFF(forward32n) { return forward32(A,dim) * recip_2; }
BZ_DECLARE_DIFF(forward42n) { return forward42(A,dim); }

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

BZ_DECLARE_MULTIDIFF(forward12n) { return forward12(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(forward22n) { return forward22(A,comp,dim); }
BZ_DECLARE_MULTIDIFF(forward32n) { return forward32(A,comp,dim) * recip_2; }
BZ_DECLARE_MULTIDIFF(forward42n) { return forward42(A,comp,dim); }

/****************************************************************************
 * Gradient operators
 ****************************************************************************/

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central12(A,firstDim),
    central12(A,secondDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central14(A,firstDim),
    central14(A,secondDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central12(A,firstDim),
    central12(A,secondDim),
    central12(A,thirdDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central14(A,firstDim),
    central14(A,secondDim),
    central14(A,thirdDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2Dn(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central12n(A,firstDim),
    central12n(A,secondDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> grad2D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central14n(A,firstDim),
    central14n(A,secondDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3Dn(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central12n(A,firstDim),
    central12n(A,secondDim),
    central12n(A,thirdDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> grad3D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central14n(A,firstDim),
    central14n(A,secondDim),
    central14n(A,thirdDim));
}

/****************************************************************************
 * Grad-squared operators
 ****************************************************************************/

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central22(A,firstDim),
    central22(A,secondDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central24(A,firstDim),
    central24(A,secondDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3D(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central22(A,firstDim),
    central22(A,secondDim),
    central22(A,thirdDim));
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3D4(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central24(A,firstDim),
    central24(A,secondDim),
    central24(A,thirdDim));
}

/****************************************************************************
 * Grad-squared operators (normalized)
 ****************************************************************************/

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2Dn(T& A) {
  return gradSqr2D(A);
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,2> gradSqr2D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,2>(
    central24(A,firstDim) * recip_12,
    central24(A,secondDim) * recip_12);
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3Dn(T& A) {
  return gradSqr3D(A);
}

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> gradSqr3D4n(T& A) {
  return TinyVector<_bz_typename T::T_numtype,3>(
    central24(A,firstDim) * recip_12,
    central24(A,secondDim) * recip_12,
    central24(A,thirdDim) * recip_12);
}

/****************************************************************************
 * Gradient operators on a vector field
 ****************************************************************************/

template<class T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
    T::T_numtype>::T_element, 3, 3>
Jacobian3D(T& A)
{
    const int x=0, y=1, z=2;
    const int u=0, v=1, w=2;

    TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
        T::T_numtype>::T_element, 3, 3> grad;

    grad(u,x) = central12(A,u,x);
    grad(u,y) = central12(A,u,y);
    grad(u,z) = central12(A,u,z);
    grad(v,x) = central12(A,v,x);
    grad(v,y) = central12(A,v,y);
    grad(v,z) = central12(A,v,z);
    grad(w,x) = central12(A,w,x);
    grad(w,y) = central12(A,w,y);
    grad(w,z) = central12(A,w,z);

    return grad;
}

template<class T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
    T::T_numtype>::T_element, 3, 3>
Jacobian3Dn(T& A)
{
    const int x=0, y=1, z=2;
    const int u=0, v=1, w=2;

    TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
        T::T_numtype>::T_element, 3, 3> grad;
    
    grad(u,x) = central12n(A,u,x);
    grad(u,y) = central12n(A,u,y);
    grad(u,z) = central12n(A,u,z);
    grad(v,x) = central12n(A,v,x);
    grad(v,y) = central12n(A,v,y);
    grad(v,z) = central12n(A,v,z);
    grad(w,x) = central12n(A,w,x);
    grad(w,y) = central12n(A,w,y);
    grad(w,z) = central12n(A,w,z);

    return grad;
}

template<class T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename
    T::T_numtype>::T_element, 3, 3>
Jacobian3D4(T& A)
{
    const int x=0, y=1, z=2;
    const int u=0, v=1, w=2;

    TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
        T::T_numtype>::T_element, 3, 3> grad;
    
    grad(u,x) = central14(A,u,x);
    grad(u,y) = central14(A,u,y);
    grad(u,z) = central14(A,u,z);
    grad(v,x) = central14(A,v,x);
    grad(v,y) = central14(A,v,y);
    grad(v,z) = central14(A,v,z);
    grad(w,x) = central14(A,w,x);
    grad(w,y) = central14(A,w,y);
    grad(w,z) = central14(A,w,z);

    return grad;
}

template<class T>
inline TinyMatrix<_bz_typename multicomponent_traits<_bz_typename
    T::T_numtype>::T_element, 3, 3>
Jacobian3D4n(T& A)
{
    const int x=0, y=1, z=2;
    const int u=0, v=1, w=2;

    TinyMatrix<_bz_typename multicomponent_traits<_bz_typename 
        T::T_numtype>::T_element, 3, 3> grad;
    
    grad(u,x) = central14n(A,u,x);
    grad(u,y) = central14n(A,u,y);
    grad(u,z) = central14n(A,u,z);
    grad(v,x) = central14n(A,v,x);
    grad(v,y) = central14n(A,v,y);
    grad(v,z) = central14n(A,v,z);
    grad(w,x) = central14n(A,w,x);
    grad(w,y) = central14n(A,w,y);
    grad(w,z) = central14n(A,w,z);

    return grad;
}

/****************************************************************************
 * Curl operators
 ****************************************************************************/

// O(h^2) curl, using central difference

template<class T>
inline TinyVector<_bz_typename T::T_numtype,3> 
curl(T& vx, T& vy, T& vz) {
  const int x = firstDim, y = secondDim, z = thirdDim;

  return TinyVector<_bz_typename T::T_numtype,3>(
    central12(vz,y)-central12(vy,z),
    central12(vx,z)-central12(vz,x),
    central12(vy,x)-central12(vx,y));
}

// Normalized O(h^2) curl, using central difference