transfer: allow grids that do not start at coordinate zero

2 files changed, 42 insertions, 9 deletions

diff --git a/transfer.c b/transfer.c
index 14b9367..40acd0a 100644
--- a/transfer.c
+++ b/transfer.c
@@ -122,13 +122,13 @@ static int transfer_factors_lagrange(GridTransferContext *ctx, unsigned int dim)
         return -ENOMEM;
 
     for (ptrdiff_t idx_dst = 0; idx_dst < ctx->dst.size[dim]; idx_dst++) {
-        const ptrdiff_t idx_src = (ptrdiff_t)(scale * idx_dst) - (priv->stencil / 2 - 1);
-        const double coord_dst  = idx_dst * step_dst;
+        const ptrdiff_t idx_src = (ptrdiff_t)(scale * (idx_dst + ctx->dst.start[dim])) - ctx->src.start[dim] - (priv->stencil / 2 - 1);
+        const double coord_dst  = (idx_dst + ctx->dst.start[dim]) * step_dst;
 
         double *fact = priv->fact[dim] + priv->stencil * idx_dst;
 
         for (int i = 0; i < priv->stencil; i++)
-            coord_src[i] = (idx_src + i) * step_src;
+            coord_src[i] = (idx_src + ctx->src.start[dim] + i) * step_src;
 
         for (int i = 0; i < priv->stencil; i++) {
             fact[i] = 1.0;
@@ -275,13 +275,16 @@ typedef struct GridTransferFW {
 
 static int fw_transfer_task(void *arg, unsigned int job_idx, unsigned int thread_idx)
 {
-    GridTransferFW *fw = arg;
-    NDArray       *dst = fw->dst;
-    const NDArray *src = fw->src;
+    GridTransferContext *ctx = arg;
+    GridTransferFW       *fw = ctx->priv;
+    NDArray             *dst = fw->dst;
+    const NDArray       *src = fw->src;
 
     const ptrdiff_t ss0 = src->stride[0];
     const ptrdiff_t ss1 = src->stride[1];
-    const double    *sd = src->data + (job_idx * 2 - fw->order) * ss0 - fw->order * ss1;
+    const double    *sd = src->data
+                          + ss0 * ((job_idx + ctx->dst.start[0]) * 2 - ctx->src.start[0] - fw->order)
+                          + ss1 * (           ctx->dst.start[1]  * 2 - ctx->src.start[1] - fw->order);
 
     double          *dd = dst->data + job_idx * dst->stride[0];
     const ptrdiff_t ds1 = dst->stride[1];
@@ -309,7 +312,7 @@ static int transfer_fw_transfer(GridTransferContext *ctx,
     fw->src = src;
     fw->dst = dst;
 
-    tp_execute(ctx->tp, ctx->dst.size[0], fw_transfer_task, fw);
+    tp_execute(ctx->tp, ctx->dst.size[0], fw_transfer_task, ctx);
 
     fw->src = NULL;
     fw->dst = NULL;
@@ -321,6 +324,12 @@ static int transfer_fw_init(GridTransferContext *ctx)
 {
     GridTransferFW *fw = ctx->priv;
 
+    for (int i = 0; i < 2; i++) {
+        /* the destination must have twice the stepsize of the source */
+        if (fabs(ctx->dst.step[i] / ctx->src.step[i] - 2.0) > 1e-14)
+            return -EINVAL;
+    }
+
     switch (ctx->op) {
     case GRID_TRANSFER_FW_1:
         fw->order   = 1;
@@ -403,7 +412,22 @@ int mg2di_gt_init(GridTransferContext *ctx)
     if (!ctx->tp)
         return -EINVAL;
 
-    /* TODO check that the destination grid is not outside the source one */
+    /* check that the destination grid is inside the source one */
+    for (int i = 0; i < 2; i++) {
+        double src_start, src_end, dst_start, dst_end;
+
+        if (ctx->src.step[i] <= 0.0 || ctx->dst.step[i] <= 0.0 ||
+            ctx->src.size[i] < 2 || ctx->dst.size[i] < 1)
+            return -EINVAL;
+
+        src_start = ctx->src.start[i] * ctx->src.step[i];
+        src_end   = (ctx->src.start[i] + ctx->src.size[i] - 1) * ctx->src.step[i];
+        dst_start = ctx->dst.start[i] * ctx->dst.step[i];
+        dst_end   = (ctx->dst.start[i] + ctx->dst.size[i] - 1) * ctx->dst.step[i];
+
+        if (src_start >= src_end || dst_start < src_start || dst_end > src_end)
+            return -EINVAL;
+    }
 
     if (t->init) {
         ret = t->init(ctx);
diff --git a/transfer.h b/transfer.h
index bb9a035..d000dcd 100644
--- a/transfer.h
+++ b/transfer.h
@@ -34,8 +34,17 @@ enum GridTransferOperator {
     GRID_TRANSFER_FW_3,
 };
 
+/**
+ * A 2-dimensional grid with regular spacing.
+ */
 typedef struct RegularGrid2D {
     /**
+     * Indices of the grid origin, relative to coordinate origin.
+     * (o0, o1) = (start[0] * step[0], start[1] * step[1])
+     */
+    ptrdiff_t           start[2];
+
+    /**
      * Number of points in the domain.
      */
     size_t              size[2];