N_pde.h

Go to the documentation of this file.

/*****************************************************************************
 *
 * MODULE:       Grass PDE Numerical Library
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
 *                 soerengebbert <at> gmx <dot> de
 *
 * PURPOSE:      This file contains definitions of variables and data types
 *
 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS
 *               for details.
 *
 *****************************************************************************/
 
#include <grass/gis.h>
#include <grass/raster3d.h>
#include <grass/glocale.h>
#include <grass/gmath.h>
 
#ifndef _N_PDE_H_
#define _N_PDE_H_
 
#define N_NORMAL_LES        0
#define N_SPARSE_LES        1
/*!
 * Boundary conditions for cells
 */
#define N_CELL_INACTIVE     0
#define N_CELL_ACTIVE       1
#define N_CELL_DIRICHLET    2
#define N_CELL_TRANSMISSION 3
/*!
 * \brief the maximum number of available cell states (eg: boundary condition,
 * inactiven active)
 * */
#define N_MAX_CELL_STATE    20
 
#define N_5_POINT_STAR      0
#define N_7_POINT_STAR      1
#define N_9_POINT_STAR      2
#define N_27_POINT_STAR     3
 
#define N_MAXIMUM_NORM      0
#define N_EUKLID_NORM       1
 
#define N_ARRAY_SUM         0 /* summ two arrays */
#define N_ARRAY_DIF         1 /* calc the difference between two arrays */
#define N_ARRAY_MUL         2 /* multiply two arrays */
#define N_ARRAY_DIV                                          \
    3 /* array division, if div with 0 the NULL value is set \
       */
 
#define N_UPWIND_FULL   0 /*full upwinding stabilization */
#define N_UPWIND_EXP    1 /*exponential upwinding stabilization */
#define N_UPWIND_WEIGHT 2 /*weighted upwinding stabilization */
 
/* *************************************************************** */
/* *************** LINEARE EQUATION SYSTEM PART ****************** */
/* *************************************************************** */
 
/*!
 * \brief The linear equation system (les) structure
 *
 * This structure manages the Ax = b system.
 * It manages regular quadratic matrices or
 * sparse matrices. The vector b and x are normal one dimensional
 * memory structures of type double. Also the number of rows
 * and the matrix type are stored in this structure.
 * */
typedef struct {
    double *x;             /*the value vector */
    double *b;             /*the right side of Ax = b */
    double **A;            /*the normal quadratic matrix */
    G_math_spvector **Asp; /*the sparse matrix */
    int rows;              /*number of rows */
    int cols;              /*number of cols */
    int quad;              /*is the matrix quadratic (1-quadratic, 0 not) */
    int type;              /*the type of the les, normal == 0, sparse == 1 */
} N_les;
 
extern N_les *N_alloc_les_param(int cols, int rows, int type, int param);
extern N_les *N_alloc_les(int rows, int type);
extern N_les *N_alloc_les_A(int rows, int type);
extern N_les *N_alloc_les_Ax(int rows, int type);
extern N_les *N_alloc_les_Ax_b(int rows, int type);
extern N_les *N_alloc_nquad_les(int cols, int rows, int type);
extern N_les *N_alloc_nquad_les_A(int cols, int rows, int type);
extern N_les *N_alloc_nquad_les_Ax(int cols, int rows, int type);
extern N_les *N_alloc_nquad_les_Ax_b(int cols, int rows, int type);
extern void N_print_les(N_les *les);
extern void N_free_les(N_les *les);
 
/* *************************************************************** */
/* *************** GEOMETRY INFORMATION ************************** */
/* *************************************************************** */
 
/*!
 * \brief Geometric information about the structured grid
 * */
typedef struct {
    int planimetric; /*If the projection is not planimetric (0), the array
                        calculation is different for each row */
    double *area;    /* the vector of area values for non-planimetric projection
                        for each row */
    int dim;         /* 2 or 3 */
 
    double dx;
    double dy;
    double dz;
 
    double Az;
 
    int depths;
    int rows;
    int cols;
 
} N_geom_data;
 
extern N_geom_data *N_alloc_geom_data(void);
extern void N_free_geom_data(N_geom_data *geodata);
extern N_geom_data *N_init_geom_data_3d(RASTER3D_Region *region3d,
                                        N_geom_data *geodata);
extern N_geom_data *N_init_geom_data_2d(struct Cell_head *region,
                                        N_geom_data *geodata);
extern double N_get_geom_data_area_of_cell(N_geom_data *geom, int row);
 
/* *************************************************************** */
/* *************** READING RASTER AND VOLUME DATA **************** */
/* *************************************************************** */
 
typedef struct {
    int type; /* which raster type CELL_TYPE, FCELL_TYPE, DCELL_TYPE */
    int rows, cols;
    int rows_intern, cols_intern;
    int offset;         /*number of cols/rows offset at each boundary */
    CELL *cell_array;   /*The data is stored in an one dimensional array
                           internally */
    FCELL *fcell_array; /*The data is stored in an one dimensional array
                           internally */
    DCELL *dcell_array; /*The data is stored in an one dimensional array
                           internally */
} N_array_2d;
 
extern N_array_2d *N_alloc_array_2d(int cols, int rows, int offset, int type);
extern void N_free_array_2d(N_array_2d *data_array);
extern int N_get_array_2d_type(N_array_2d *array2d);
extern void N_get_array_2d_value(N_array_2d *array2d, int col, int row,
                                 void *value);
extern CELL N_get_array_2d_c_value(N_array_2d *array2d, int col, int row);
extern FCELL N_get_array_2d_f_value(N_array_2d *array2d, int col, int row);
extern DCELL N_get_array_2d_d_value(N_array_2d *array2d, int col, int row);
extern void N_put_array_2d_value(N_array_2d *array2d, int col, int row,
                                 char *value);
extern void N_put_array_2d_c_value(N_array_2d *array2d, int col, int row,
                                   CELL value);
extern void N_put_array_2d_f_value(N_array_2d *array2d, int col, int row,
                                   FCELL value);
extern void N_put_array_2d_d_value(N_array_2d *array2d, int col, int row,
                                   DCELL value);
extern int N_is_array_2d_value_null(N_array_2d *array2d, int col, int row);
extern void N_put_array_2d_value_null(N_array_2d *array2d, int col, int row);
extern void N_print_array_2d(N_array_2d *data);
extern void N_print_array_2d_info(N_array_2d *data);
extern void N_copy_array_2d(N_array_2d *source, N_array_2d *target);
extern double N_norm_array_2d(N_array_2d *array1, N_array_2d *array2, int type);
extern N_array_2d *N_math_array_2d(N_array_2d *array1, N_array_2d *array2,
                                   N_array_2d *result, int type);
extern int N_convert_array_2d_null_to_zero(N_array_2d *a);
extern N_array_2d *N_read_rast_to_array_2d(char *name, N_array_2d *array);
extern void N_write_array_2d_to_rast(N_array_2d *array, char *name);
extern void N_calc_array_2d_stats(N_array_2d *a, double *min, double *max,
                                  double *sum, int *nonzero, int withoffset);
 
typedef struct {
    int type; /* which raster type FCELL_TYPE, DCELL_TYPE */
    int rows, cols, depths;
    int rows_intern, cols_intern, depths_intern;
    int offset;          /*number of cols/rows/depths offset at each boundary */
    float *fcell_array;  /*The data is stored in an one dimensional array
                            internally */
    double *dcell_array; /*The data is stored in an one dimensional array
                            internally */
} N_array_3d;
 
extern N_array_3d *N_alloc_array_3d(int cols, int rows, int depths, int offset,
                                    int type);
extern void N_free_array_3d(N_array_3d *data_array);
extern int N_get_array_3d_type(N_array_3d *array3d);
extern void N_get_array_3d_value(N_array_3d *array3d, int col, int row,
                                 int depth, void *value);
extern float N_get_array_3d_f_value(N_array_3d *array3d, int col, int row,
                                    int depth);
extern double N_get_array_3d_d_value(N_array_3d *array3d, int col, int row,
                                     int depth);
extern void N_put_array_3d_value(N_array_3d *array3d, int col, int row,
                                 int depth, char *value);
extern void N_put_array_3d_f_value(N_array_3d *array3d, int col, int row,
                                   int depth, float value);
extern void N_put_array_3d_d_value(N_array_3d *array3d, int col, int row,
                                   int depth, double value);
extern int N_is_array_3d_value_null(N_array_3d *array3d, int col, int row,
                                    int depth);
extern void N_put_array_3d_value_null(N_array_3d *array3d, int col, int row,
                                      int depth);
extern void N_print_array_3d(N_array_3d *data);
extern void N_print_array_3d_info(N_array_3d *data);
extern void N_copy_array_3d(N_array_3d *source, N_array_3d *target);
extern double N_norm_array_3d(N_array_3d *array1, N_array_3d *array2, int type);
extern N_array_3d *N_math_array_3d(N_array_3d *array1, N_array_3d *array2,
                                   N_array_3d *result, int type);
extern int N_convert_array_3d_null_to_zero(N_array_3d *a);
extern N_array_3d *N_read_rast3d_to_array_3d(char *name, N_array_3d *array,
                                             int mask);
extern void N_write_array_3d_to_rast3d(N_array_3d *array, char *name, int mask);
extern void N_calc_array_3d_stats(N_array_3d *a, double *min, double *max,
                                  double *sum, int *nonzero, int withoffset);
 
/* *************************************************************** */
/* *************** MATRIX ASSEMBLING METHODS ********************* */
/* *************************************************************** */
/*!
 * \brief Matrix entries for a mass balance 5/7/9 star system
 *
 * Matrix entries for the mass balance of a 5 star system
 *
 * The entries are center, east, west, north, south and the
 * right side vector b of Ax = b. This system is typically used in 2d.
 
 \verbatim
 N
 |
 W-- C --E
 |
 S
 \endverbatim
 
 * Matrix entries for the mass balance of a 7 star system
 *
 * The entries are center, east, west, north, south, top, bottom and the
 * right side vector b of Ax = b. This system is typically used in 3d.
 
 \verbatim
 T N
 |/
 W-- C --E
 /|
 S B
 \endverbatim
 
 * Matrix entries for the mass balance of a 9 star system
 *
 * The entries are center, east, west, north, south, north-east, south-east,
 * north-wast, south-west and the
 * right side vector b of Ax = b. This system is typically used in 2d.
 
 \verbatim
 NW  N  NE
 \ | /
 W-- C --E
 / | \
 SW  S  SE
 \endverbatim
 
 * Matrix entries for the mass balance of a 27 star system
 *
 * The entries are center, east, west, north, south, north-east, south-east,
 * north-wast, south-west, same for top and bottom and the
 * right side vector b of Ax = b. This system is typically used in 2d.
 
 \verbatim
 top:
 NW_T N_Z NE_T
 \ | /
 W_T-- T --E_T
 / | \
 SW_T S_T SE_T
 
 center:
 NW  N  NE
 \ | /
 W-- C --E
 / | \
 SW  S  SE
 
 bottom:
 NW_B N_B NE_B
 \ | /
 W_B-- B --E_B
 / | \
 SW_B S_B SE_B
 \endverbatim
 
 */
typedef struct {
    int type;
    int count;
    double C, W, E, N, S, NE, NW, SE, SW, V;
    /*top part */
    double T, W_T, E_T, N_T, S_T, NE_T, NW_T, SE_T, SW_T;
    /*bottom part */
    double B, W_B, E_B, N_B, S_B, NE_B, NW_B, SE_B, SW_B;
} N_data_star;
 
/*!
 * \brief callback structure for 3d matrix assembling
 * */
typedef struct {
    N_data_star *(*callback)(void *, N_geom_data *, int, int, int);
} N_les_callback_3d;
 
/*!
 * \brief callback structure for 2d matrix assembling
 * */
typedef struct {
    N_data_star *(*callback)(void *, N_geom_data *, int, int);
} N_les_callback_2d;
 
extern void N_set_les_callback_3d_func(
    N_les_callback_3d *data,
    N_data_star *(*callback_func_3d)(void *, N_geom_data *, int, int, int));
extern void N_set_les_callback_2d_func(
    N_les_callback_2d *data,
    N_data_star *(*callback_func_2d)(void *, N_geom_data *, int, int));
extern N_les_callback_3d *N_alloc_les_callback_3d(void);
extern N_les_callback_2d *N_alloc_les_callback_2d(void);
extern N_data_star *N_alloc_5star(void);
extern N_data_star *N_alloc_7star(void);
extern N_data_star *N_alloc_9star(void);
extern N_data_star *N_alloc_27star(void);
extern N_data_star *N_create_5star(double C, double W, double E, double N,
                                   double S, double V);
extern N_data_star *N_create_7star(double C, double W, double E, double N,
                                   double S, double T, double B, double V);
extern N_data_star *N_create_9star(double C, double W, double E, double N,
                                   double S, double NW, double SW, double NE,
                                   double SE, double V);
extern N_data_star *
N_create_27star(double C, double W, double E, double N, double S, double NW,
                double SW, double NE, double SE, double T, double W_T,
                double E_T, double N_T, double S_T, double NW_T, double SW_T,
                double NE_T, double SE_T, double B, double W_B, double E_B,
                double N_B, double S_B, double NW_B, double SW_B, double NE_B,
                double SE_B, double V);
extern N_data_star *N_callback_template_3d(void *data, N_geom_data *geom,
                                           int col, int row, int depth);
extern N_data_star *N_callback_template_2d(void *data, N_geom_data *geom,
                                           int col, int row);
extern N_les *N_assemble_les_3d(int les_type, N_geom_data *geom,
                                N_array_3d *status, N_array_3d *start_val,
                                void *data, N_les_callback_3d *callback);
extern N_les *N_assemble_les_3d_active(int les_type, N_geom_data *geom,
                                       N_array_3d *status,
                                       N_array_3d *start_val, void *data,
                                       N_les_callback_3d *callback);
extern N_les *N_assemble_les_3d_dirichlet(int les_type, N_geom_data *geom,
                                          N_array_3d *status,
                                          N_array_3d *start_val, void *data,
                                          N_les_callback_3d *callback);
extern N_les *N_assemble_les_3d_param(int les_type, N_geom_data *geom,
                                      N_array_3d *status, N_array_3d *start_val,
                                      void *data, N_les_callback_3d *callback,
                                      int cell_type);
extern N_les *N_assemble_les_2d(int les_type, N_geom_data *geom,
                                N_array_2d *status, N_array_2d *start_val,
                                void *data, N_les_callback_2d *callback);
extern N_les *N_assemble_les_2d_active(int les_type, N_geom_data *geom,
                                       N_array_2d *status,
                                       N_array_2d *start_val, void *data,
                                       N_les_callback_2d *callback);
extern N_les *N_assemble_les_2d_dirichlet(int les_type, N_geom_data *geom,
                                          N_array_2d *status,
                                          N_array_2d *start_val, void *data,
                                          N_les_callback_2d *callback);
extern N_les *N_assemble_les_2d_param(int les_type, N_geom_data *geom,
                                      N_array_2d *status, N_array_2d *start_val,
                                      void *data, N_les_callback_2d *callback,
                                      int cell_Type);
extern int N_les_pivot_create(N_les *les);
int N_les_integrate_dirichlet_2d(N_les *les, N_geom_data *geom,
                                 N_array_2d *status, N_array_2d *start_val);
int N_les_integrate_dirichlet_3d(N_les *les, N_geom_data *geom,
                                 N_array_3d *status, N_array_3d *start_val);
 
/* *************************************************************** */
/* *************** GPDE STANDARD OPTIONS ************************* */
/* *************************************************************** */
 
/*! \brief Standard options of the gpde library
 * */
typedef enum {
    N_OPT_SOLVER_SYMM,     /*! solver for symmetric, positive definite linear
                              equation systems */
    N_OPT_SOLVER_UNSYMM,   /*! solver for unsymmetric linear equation systems */
    N_OPT_MAX_ITERATIONS,  /*! Maximum number of iteration used to solver the
                              linear equation system */
    N_OPT_ITERATION_ERROR, /*! Error break criteria for the iterative solver
                              (jacobi, sor, cg or bicgstab) */
    N_OPT_SOR_VALUE, /*! The relaxation parameter used by the jacobi and sor
                        solver for speedup or stabilizing */
    N_OPT_CALC_TIME  /*! The calculation time in seconds */
} N_STD_OPT;
 
extern struct Option *N_define_standard_option(int opt);
 
/* *************************************************************** */
/* *************** GPDE MATHEMATICAL TOOLS *********************** */
/* *************************************************************** */
 
extern double N_calc_arith_mean(double a, double b);
extern double N_calc_arith_mean_n(double *a, int size);
extern double N_calc_geom_mean(double a, double b);
extern double N_calc_geom_mean_n(double *a, int size);
extern double N_calc_harmonic_mean(double a, double b);
extern double N_calc_harmonic_mean_n(double *a, int size);
extern double N_calc_quad_mean(double a, double b);
extern double N_calc_quad_mean_n(double *a, int size);
 
/* *************************************************************** */
/* *************** UPWIND STABILIZATION ALGORITHMS *************** */
/* *************************************************************** */
 
extern double N_full_upwinding(double sprod, double distance, double D);
extern double N_exp_upwinding(double sprod, double distance, double D);
 
/* *************************************************************** */
/* *************** METHODS FOR GRADIENT CALCULATION ************** */
/* *************************************************************** */
/*!
   \verbatim
 
   ______________
   |    |    |    |
   |    |    |    |
   |----|-NC-|----|
   |    |    |    |
   |   WC    EC   |
   |    |    |    |
   |----|-SC-|----|
   |    |    |    |
   |____|____|____|
 
 
   |  /
   TC NC
   |/
   --WC-----EC--
   /|
   SC BC
   /  |
 
   \endverbatim
 
 */
 
/*! \brief Gradient between the cells in X and Y direction */
typedef struct {
 
    double NC, SC, WC, EC;
 
} N_gradient_2d;
 
/*! \brief Gradient between the cells in X, Y and Z direction */
typedef struct {
 
    double NC, SC, WC, EC, TC, BC;
 
} N_gradient_3d;
 
/*!
   \verbatim
 
   Gradient in X direction between the cell neighbours
   ____ ____ ____
   |    |    |    |
   |   NWN  NEN   |
   |____|____|____|
   |    |    |    |
   |   WN    EN   |
   |____|____|____|
   |    |    |    |
   |   SWS  SES   |
   |____|____|____|
 
   Gradient in Y direction between the cell neighbours
   ______________
   |    |    |    |
   |    |    |    |
   |NWW-|-NC-|-NEE|
   |    |    |    |
   |    |    |    |
   |SWW-|-SC-|-SEE|
   |    |    |    |
   |____|____|____|
 
   Gradient in Z direction between the cell neighbours
   /______________/
   /|    |    |    |
   | NWZ| NZ | NEZ|
   |____|____|____|
   /|    |    |    |
   | WZ | CZ | EZ |
   |____|____|____|
   /|    |    |    |
   | SWZ| SZ | SEZ|
   |____|____|____|
   /____/____/____/
 
 
   \endverbatim
 */
 
/*! \brief Gradient between the cell neighbours in X direction */
typedef struct {
 
    double NWN, NEN, WC, EC, SWS, SES;
 
} N_gradient_neighbours_x;
 
/*! \brief Gradient between the cell neighbours in Y direction */
typedef struct {
 
    double NWW, NEE, NC, SC, SWW, SEE;
 
} N_gradient_neighbours_y;
 
/*! \brief Gradient between the cell neighbours in Z direction */
typedef struct {
 
    double NWZ, NZ, NEZ, WZ, CZ, EZ, SWZ, SZ, SEZ;
 
} N_gradient_neighbours_z;
 
/*! \brief Gradient between the cell neighbours in X and Y direction */
typedef struct {
 
    N_gradient_neighbours_x *x;
    N_gradient_neighbours_y *y;
 
} N_gradient_neighbours_2d;
 
/*! \brief Gradient between the cell neighbours in X, Y and Z direction */
typedef struct {
 
    N_gradient_neighbours_x *xt; /*top values */
    N_gradient_neighbours_x *xc; /*center values */
    N_gradient_neighbours_x *xb; /*bottom values */
 
    N_gradient_neighbours_y *yt; /*top values */
    N_gradient_neighbours_y *yc; /*center values */
    N_gradient_neighbours_y *yb; /*bottom values */
 
    N_gradient_neighbours_z *zt; /*top-center values */
    N_gradient_neighbours_z *zb; /*bottom-center values */
 
} N_gradient_neighbours_3d;
 
/*! Two dimensional gradient field */
typedef struct {
 
    N_array_2d *x_array;
    N_array_2d *y_array;
    int cols, rows;
    double min, max, mean, sum;
    int nonull;
 
} N_gradient_field_2d;
 
/*! Three dimensional gradient field */
typedef struct {
 
    N_array_3d *x_array;
    N_array_3d *y_array;
    N_array_3d *z_array;
    int cols, rows, depths;
    double min, max, mean, sum;
    int nonull;
 
} N_gradient_field_3d;
 
extern N_gradient_2d *N_alloc_gradient_2d(void);
extern void N_free_gradient_2d(N_gradient_2d *grad);
extern N_gradient_2d *N_create_gradient_2d(double NC, double SC, double WC,
                                           double EC);
extern int N_copy_gradient_2d(N_gradient_2d *source, N_gradient_2d *target);
extern N_gradient_2d *N_get_gradient_2d(N_gradient_field_2d *field,
                                        N_gradient_2d *gradient, int col,
                                        int row);
extern N_gradient_3d *N_alloc_gradient_3d(void);
extern void N_free_gradient_3d(N_gradient_3d *grad);
extern N_gradient_3d *N_create_gradient_3d(double NC, double SC, double WC,
                                           double EC, double TC, double BC);
extern int N_copy_gradient_3d(N_gradient_3d *source, N_gradient_3d *target);
extern N_gradient_3d *N_get_gradient_3d(N_gradient_field_3d *field,
                                        N_gradient_3d *gradient, int col,
                                        int row, int depth);
extern N_gradient_neighbours_x *N_alloc_gradient_neighbours_x(void);
extern void N_free_gradient_neighbours_x(N_gradient_neighbours_x *grad);
extern N_gradient_neighbours_x *
N_create_gradient_neighbours_x(double NWN, double NEN, double WC, double EC,
                               double SWS, double SES);
extern int N_copy_gradient_neighbours_x(N_gradient_neighbours_x *source,
                                        N_gradient_neighbours_x *target);
extern N_gradient_neighbours_y *N_alloc_gradient_neighbours_y(void);
extern void N_free_gradient_neighbours_y(N_gradient_neighbours_y *grad);
extern N_gradient_neighbours_y *
N_create_gradient_neighbours_y(double NWW, double NEE, double NC, double SC,
                               double SWW, double SEE);
extern int N_copy_gradient_neighbours_y(N_gradient_neighbours_y *source,
                                        N_gradient_neighbours_y *target);
extern N_gradient_neighbours_z *N_alloc_gradient_neighbours_z(void);
extern void N_free_gradient_neighbours_z(N_gradient_neighbours_z *grad);
extern N_gradient_neighbours_z *
N_create_gradient_neighbours_z(double NWZ, double NZ, double NEZ, double WZ,
                               double CZ, double EZ, double SWZ, double SZ,
                               double SEZ);
extern int N_copy_gradient_neighbours_z(N_gradient_neighbours_z *source,
                                        N_gradient_neighbours_z *target);
extern N_gradient_neighbours_2d *N_alloc_gradient_neighbours_2d(void);
extern void N_free_gradient_neighbours_2d(N_gradient_neighbours_2d *grad);
extern N_gradient_neighbours_2d *
N_create_gradient_neighbours_2d(N_gradient_neighbours_x *x,
                                N_gradient_neighbours_y *y);
extern int N_copy_gradient_neighbours_2d(N_gradient_neighbours_2d *source,
                                         N_gradient_neighbours_2d *target);
extern N_gradient_neighbours_2d *
N_get_gradient_neighbours_2d(N_gradient_field_2d *field,
                             N_gradient_neighbours_2d *gradient, int col,
                             int row);
extern N_gradient_neighbours_3d *N_alloc_gradient_neighbours_3d(void);
extern void N_free_gradient_neighbours_3d(N_gradient_neighbours_3d *grad);
extern N_gradient_neighbours_3d *N_create_gradient_neighbours_3d(
    N_gradient_neighbours_x *xt, N_gradient_neighbours_x *xc,
    N_gradient_neighbours_x *xb, N_gradient_neighbours_y *yt,
    N_gradient_neighbours_y *yc, N_gradient_neighbours_y *yb,
    N_gradient_neighbours_z *zt, N_gradient_neighbours_z *zb);
extern int N_copy_gradient_neighbours_3d(N_gradient_neighbours_3d *source,
                                         N_gradient_neighbours_3d *target);
extern void N_print_gradient_field_2d_info(N_gradient_field_2d *field);
extern void N_calc_gradient_field_2d_stats(N_gradient_field_2d *field);
extern N_gradient_field_2d *N_alloc_gradient_field_2d(int cols, int rows);
extern void N_free_gradient_field_2d(N_gradient_field_2d *field);
extern int N_copy_gradient_field_2d(N_gradient_field_2d *source,
                                    N_gradient_field_2d *target);
extern N_gradient_field_2d *
N_compute_gradient_field_2d(N_array_2d *pot, N_array_2d *weight_x,
                            N_array_2d *weight_y, N_geom_data *geom,
                            N_gradient_field_2d *gradfield);
extern void N_compute_gradient_field_components_2d(N_gradient_field_2d *field,
                                                   N_array_2d *x_comp,
                                                   N_array_2d *y_comp);
extern void N_print_gradient_field_3d_info(N_gradient_field_3d *field);
extern void N_calc_gradient_field_3d_stats(N_gradient_field_3d *field);
extern N_gradient_field_3d *N_alloc_gradient_field_3d(int cols, int rows,
                                                      int depths);
extern void N_free_gradient_field_3d(N_gradient_field_3d *field);
extern int N_copy_gradient_field_3d(N_gradient_field_3d *source,
                                    N_gradient_field_3d *target);
extern N_gradient_field_3d *
N_compute_gradient_field_3d(N_array_3d *pot, N_array_3d *weight_x,
                            N_array_3d *weight_y, N_array_3d *weight_z,
                            N_geom_data *geom, N_gradient_field_3d *gradfield);
extern void N_compute_gradient_field_components_3d(N_gradient_field_3d *field,
                                                   N_array_3d *x_comp,
                                                   N_array_3d *y_comp,
                                                   N_array_3d *z_comp);
 
#endif