gridspec_api/regrid/nccf_compute_regrid_weights.c

#include <math.h>
#include <netcdf.h>

#include <nccf_coord.h>
#include <nccf_grid.h>
#include <nccf_utility_functions.h>
#include "nccf_regrid.h"

#define MULTIPLE_TRIALS 1
#define SNAKE 0
#define EXHAUSTIVE_SEARCH 0

struct CFLISTITEM *CFLIST_REGRID;

int nccf_is_forbidden(int ndims, const double dIndices[], 
                      struct CFLISTITEM *box_lohi) {
  int res = 0;
  int in_this_box;
  int i, id;
  int *lohi;
  nccf_li_begin(&box_lohi);
  while (nccf_li_next(&box_lohi)) {
    id = nccf_li_get_id(&box_lohi);
    lohi = nccf_li_find(&box_lohi, id);
    in_this_box = 1;
    for (i = 0; i < ndims; ++i) {
      in_this_box *= (dIndices[i] >= lohi[i] && dIndices[i] <= lohi[i + ndims]);
    }
    res += in_this_box;
  }
  return res;
}

int nccf_find_indices_exhaustive(int ndims, const int oriDims[], 
                                 const double ** coordOriData, 
                                 const int is_periodic[], 
                                 const double xyz[], 
                                 int nitermax, double tolpos, 
                                 double dIndices[]) {
  int status = 1;
  int k, i;
  int ntot = 1;
  int ijk[ndims];
  for (i = 0; i < ndims; ++i) {
    ntot *= oriDims[i];
  }
  for (k = 0; k < ntot; ++k) {
    nccf_get_multi_index(ndims, oriDims, k, ijk);
    for (i = 0; i < ndims; ++i) {
      dIndices[i] = (double) ijk[i];
    }
    status = nccf_find_indices_double(ndims, oriDims,
                                      coordOriData,
                                      is_periodic,
                                      xyz, nitermax, tolpos,
                                      NULL, dIndices); 
    if (!status) {
      break;
    }
  }
  return status;
}

int nccf_find_indices_from_corners(int ndims, const int oriDims[], 
                                   const double ** coordOriData, 
                                   const int is_periodic[], 
                                   const double xyz[], 
                                   int nitermax, double tolpos, 
                                   double dIndices[]) {
  int status = 0;
  int numSides = 1;
  int j, i;
  int cornerVector[ndims];
  for (i = 0; i < ndims; ++i) {
    numSides *= 3;
  }
  for (j = 0; j < numSides; ++j) {
    nccf_index_to_corner_vector(j, ndims, cornerVector);
    for (i = 0; i < ndims; ++i) {
      dIndices[i] = dIndices[i] = (cornerVector[i] + 1) * (oriDims[i] - 1) / 2;
    }
    status = nccf_find_indices_double(ndims, oriDims,
                                      (const double **) coordOriData,
                                      is_periodic,
                                      xyz, nitermax, tolpos,
                                      NULL, dIndices);
    if (!status) {
      /* success */
      break;
    }
  }
  return status;
}

int nccf_compute_regrid_weights(int regrid_id,
                                int nitermax, 
                                double tolpos, 
                                const int is_periodic[]) {

  int totErr = NC_NOERR;
  int status;
  int dims[] = {0, 0}, indom_dim[] = {0};
  const char *dim_names[] = {CF_DIMNAME_NTARGETS, CF_DIMNAME_NNODES};
  const char *indom_name[] = {CF_DIMNAME_NNODES};

  struct nccf_regrid_type *self;
  self = nccf_li_find(&CFLIST_REGRID, regrid_id);

  /* Fix */
  nitermax = (nitermax <= 0? 1: nitermax);
  tolpos = (tolpos < 1.e-8? 1.e-8: tolpos);
  
  int coordTargetIds[self->ndims];
  status = nccf_get_struct_grid_coord_ids(self->tgt_grid_id, coordTargetIds);
  totErr += abs(status);

  /* assume ndims >= 1 */
  int targetDims[self->ndims];
  status = nccf_get_coord_dims(coordTargetIds[0], targetDims);
  totErr += abs(status);

  
  int i, j, k;
  int numNodes = 1;
  self->ntargets = 1;
  for (i = 0; i < self->ndims; ++i) {
    self->ntargets *= targetDims[i];
    numNodes *= 2;
  }

  /* Get the original grid coordinates */
  int coordOriIds[self->ndims];
  status = nccf_get_struct_grid_coord_ids(self->ori_grid_id, coordOriIds);
  totErr += abs(status);
  
  int oriDims[self->ndims];
  status = nccf_get_coord_dims(coordOriIds[0], oriDims);
  totErr += abs(status);  

  double **coordOriData;
  double **coordTargetData;
  coordOriData = (double **) malloc(self->ndims * sizeof(double *));
  coordTargetData = (double **) malloc(self->ndims * sizeof(double *));

  for (i = 0; i < self->ndims; ++i) {
    double *dataPtr;
    status = nccf_get_coord_data_pointer(coordOriIds[i], &dataPtr);
    totErr += abs(status);
    coordOriData[i] = dataPtr;
    status = nccf_get_coord_data_pointer(coordTargetIds[i], &dataPtr);
    totErr += abs(status);
    coordTargetData[i] = dataPtr;
  }  

  /* Locate the cells and compute the weights */
  double dIndices[self->ndims];
  for (i = 0; i < self->ndims; ++i) {
    /* initial guess */
    dIndices[i] = 0.5; 
  }

  self->weights = (double **) malloc(self->ntargets * sizeof(double *));
  self->indices = (int **) malloc(self->ntargets * sizeof(int *));
  self->inside_domain = (int *) malloc(self->ntargets * sizeof(int));

  int displ[self->ndims];
  int indx[self->ndims];

  int prodDims[self->ndims];
  prodDims[self->ndims - 1] = 1; 
  for (i = self->ndims - 2; i >= 0; --i) {
    prodDims[i] = prodDims[i + 1] * 2;
  }

  double xyz[self->ndims];
  int countFirstHit = 0;
  int countMultipleTrials = 0;
#if SNAKE == 1
  int ijk[self->ndims];
  for (i = 0; i < self->ndims; ++i) {
    ijk[i] = 0;
  }
  
  k = 0;
  int inside_domain = 1;
  while (inside_domain) {
    k = nccf_get_flat_index(self->ndims, targetDims, ijk);
#else
  for (k = 0; k < self->ntargets; ++k) {
#endif // SNAKE
    for (i = 0; i < self->ndims; ++i) {
      xyz[i] = coordTargetData[i][k];
    }
    status = nccf_find_indices_double(self->ndims, oriDims, 
                                      (const double **) coordOriData, is_periodic, 
                                      xyz, nitermax, tolpos, 
                                      NULL, dIndices);
    if (!status && nccf_is_forbidden(self->ndims, dIndices, self->box_lohi)) {
      status = -1;
    }
    if (!status) {
      countFirstHit++;
    }
#if MULTIPLE_TRIALS == 1
    /*
       A non-zero value indicates failure to find the set of (float) indices.
       This may be due to the target position being outside of the domain or
       the position was found to be in the forbidden region.
    */
    if (status) {
      status = nccf_find_indices_from_corners(self->ndims, oriDims,
                                              (const double **) coordOriData,
                                              is_periodic,
                                              xyz, nitermax, tolpos,
                                              dIndices);

#if EXHAUSTIVE_SEARCH == 1
      if (status) {
        status = nccf_find_indices_exhaustive(self->ndims, oriDims,
                                              (const double **) coordOriData,
                                              is_periodic,
                                              xyz, nitermax, tolpos,
                                              dIndices);
      }
#endif
      if (!status && nccf_is_forbidden(self->ndims, dIndices, self->box_lohi)) {
          status = -1;
      }
      if (!status) {
        /* success */
        countMultipleTrials++;
      }
    }
#endif // MULTIPLE_TRIALS

    self->inside_domain[k] = 0;
    if (status == 0) {
      self->inside_domain[k] = 1;
      self->n_nonmasked++;
    }

    // compute weights
    self->weights[k] = (double *) malloc(numNodes * sizeof(double));
    status = nccf_get_linear_interp_weights(self->ndims, oriDims, 
                                            dIndices, self->weights[k]);
    totErr += abs(status);

    // compute indices
    self->indices[k] = (int *) malloc(numNodes * sizeof(int));
    for (j = 0; j < numNodes; ++j) {

      // fill in displ,  is a vector of zeros and ones
      for (i = 0; i < self->ndims; ++i) {
        displ[i] = j / prodDims[i] % 2;
      }

      // compute the index set of the node
      for (i = 0; i < self->ndims; ++i) {
        indx[i] = (int) floor(dIndices[i]) + displ[i];
        // Fix index to prevent memory errors. Note this should not
        // influence the interpolation because the weights should be 
        // zero for nodal values outside the domain. 
        indx[i] = (indx[i] > oriDims[i] - 1? oriDims[i] - 1: indx[i]);
      }

      self->indices[k][j] = nccf_get_flat_index(self->ndims, oriDims, indx);
    }
#if SNAKE == 1
    inside_domain = nccf_find_next_adjacent(self->ndims, targetDims, ijk);
#endif
  }

/* Populate the data structures initialized in nccf_def_regrid */
  dims[0] = self->ntargets;
  dims[1] = numNodes;
  self->nnodes  =  numNodes;
  indom_dim[0] = self->ntargets;
  int ddims = 2;  /* Number of dimensions for writing. NOT ndims */

  nccf_varSetDims( &self->inside_domain_stt, 1, indom_dim, indom_name );
  nccf_varSetDataPtr( &self->inside_domain_stt, NC_INT, self->inside_domain );
  nccf_varSetDims( &self->weights_stt, ddims, dims, dim_names );
  nccf_varSetDataPtr( &self->weights_stt, NC_DOUBLE, self->weights );
  nccf_varSetDims( &self->indices_stt, ddims, dims, dim_names );
  nccf_varSetDataPtr( &self->indices_stt, NC_INT, self->indices );

  /* Clean up */
  free(coordOriData);
  free(coordTargetData);
  
  return totErr;
}