/* Last edited on 2011-09-24 21:48:13 by stolfilocal */

bz_patch_ddim_t bz_patch_get_ddim(bz_patch_t *P)
  { assert(P->C.ds.na = P->d + P->r);
    return P->d;
  }
  
bz_patch_rinds_t bz_patch_get_rinds(bz_patch_t *P)
  { assert(P->C.ds.na = P->d + P->r);
    return P->r;
  }
  
void bz_patch_compute_steps
  ( bz_patch_t *P,  /* A Bézier patch. */
    int step[],     /* (OUT) control point index increment for each axis. */
    int *sz         /* (OUT) total number of control points. */
  )
  {
    int tsz = 1, i;
    bz_degree_t *g = &(P->g[0]);
    for (i = P->d-1; i >= 0; i--) { step[i] = tsz; tsz *= (g[i] + 1); }
    (*sz) = tsz;
  }

bz_patch_t bz_patch_facet_new(bz_patch_t *P, box_axis_index_t j)
  { 
    
    bz_degree_t fg[bz_patch_MAX_DDIM];
       
    ix_size_t fsz[];
    int bi, fi;
    for (bi = 0, fi = 0; bi < d; bi++)
      { if (bi != j) { fsz[fi] = g[bi]; fi++; } }
    return bz_patch_new(P->d-1, fsz, );
  }

void bz_patch_set_face
  ( bz_patch_t *P, 
    box_signed_dir_t dir[], 
    bz_patch_t *t
  )
  { int tm = t->d;
    int bm = P->d;
    const bz_degree_t *bg = &(P->g[0]);
    const bz_degree_t *tg = &(t->g[0]);
    affirm(t->n == P->n, "range dimension mismatch");
    { int it, ib;
      for (ib = 0,it = 0; ib < P->d; ib++)
        { if (dir[ib] == 00) 
            { affirm(tg[it] == bg[ib], "degree mismatch");
              it++; 
            }
        }
      affirm(tm == it, "domain dimension mismatch");
    }
    { int n = P->n;
      static int bstep[bz_patch_MAX_DDIM];
      static int tstep[bz_patch_MAX_DDIM];
      int bsz, tsz, bix;
      
      bz_patch_compute_steps(P, bstep, &bsz);
      bz_patch_compute_steps(t, tstep, &tsz);
      /* Modify control points: */
      for (bix = 0; bix < bsz; bix++)
        { /* Compute {tix} such that {t->c[tix]} is  the control point of
            {t} closest to {P->c[ix]}; compute the index {btix} such that
            {P->c[btix]} is the corresponding control point in {P};
            and compute the influence factor {dfac} of that control 
            point on {P->c[bix]}: */
          int bb = bix, tix = 0, btix = 0, it, ib;
          double dfac = 1.0;
          for (ib = bm-1, it = tm-1; ib >= 0; ib--) 
            { int bnumi = bg[ib] + 1, bei = bb % bnumi;
              double xi = ((double)bei)/((double)bg[ib]);
              bb /= bnumi;
              if (dir[ib] == -1)
                { dfac *= (1.0 - xi); /* btix += 0*bstep[ib]; */ }
              else if (dir[ib] == +1)
                { dfac *= xi; btix += bg[ib]*bstep[ib]; }
              else
                { btix += bei*bstep[ib]; tix += bei*tstep[it]; it--; }
            }
          if (dfac > 0.0)
            { /* Modify {P->c[bix]} by {t->c[tix]-P->c[btix]} times {dfac}: */
              double *bc = &(P->c[n*bix]);
              double *tc = &(t->c[n*tix]);
              int j;
              if (dfac == 1.0)
                { affirm(bix == btix, "indexing bug");
                  for (j = 0; j < n; j++,tc++,bc++) { (*bc) = (*tc); }
                }
              else
                { double *btc = &(P->c[n*btix]);
                  affirm(bix != btix, "indexing bug");
                  for (j = 0; j < n; j++,tc++,bc++,btc++)
                    { double delta = dfac*((*tc) - (*btc));
                      (*bc) += delta;
                    } 
                }
              
            }
        }
    }
  }
  
void bz_patch_get_face ( bz_patch_t *P, box_signed_dir_t dir[], bz_patch_t *T )
  { int tm = T->d;
    int bm = P->d;
    const bz_degree_t *bg = &(P->g[0]);
    const bz_degree_t *tg = &(t->g[0]);
    affirm(t->n == P->n, "range dimension mismatch");
    { int it, ib;
      for (ib = 0,it = 0; ib < P->d; ib++)
        { if (dir[ib] == 00) 
            { affirm(it < tm, "dimension too low for location");
              affirm(tg[it] == bg[ib], "degree mismatch");
              it++; 
            }
        }
      affirm(tm == it, "dimension inconsistent with location");
    }
    { int n = P->n;
      static int bstep[bz_patch_MAX_DDIM];
      static int tstep[bz_patch_MAX_DDIM];
      int bsz, tsz, tix;
      
      bz_patch_compute_steps(P, bstep, &bsz);
      bz_patch_compute_steps(t, tstep, &tsz);
      /* Copy control points: */
      for (tix = 0; tix < tsz; tix++)
        { /* Compute index {bix} such that {P->c[bix]} corresponds to {t->c[tix]}: */
          int bix = 0, tt = tix, it, ib;
          for (ib = bm-1, it = tm-1; ib >= 0; ib--) 
            { if (dir[ib] == -1) 
                { /* bix += 0*bstep[ib]; */ }
              else if (dir[ib] == +1) 
                { bix += bg[ib]*bstep[ib]; }
              else
                { int tnumi = tg[it]+1, tei = tt % tnumi;
                  bix += tei*bstep[ib];
                  tt /= tnumi; it--;
                }
            }
          /* Copy the coefficient vector: */
          { double *tc = &(t->c[n*tix]), *bc = &(P->c[n*bix]); 
            int j;
            for (j = 0; j < n; j++,tc++,bc++) { (*tc) = (*bc); }
          }
        }
    }
  }