/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>

 
 
#include <dx/dx.h>
#include <ctype.h>
#include <math.h>

/* struct for passing parameters around
 */ 
struct argblk {
    Object toslice;             /* object to process */
    int slicedim;               /* dimension to make slices along */
    int slicecount;             /* number of slices to make */
    int *sliceplace;            /* connection nums at which to make slice */
    int numdim;                 /* dimensionality of input object */
    Object compositeseries;     /* if !NULL, multiple slices through CF */
    int compositecount;         /* number of connections in complete CF */
    float *compositepositions;  /* positions list for complete CF */
    Object outgroup;            /* if !NULL, parent group for slices */
    char *membername;		/* if !NULL, slice name */
    float seriesplace;		/* if series, series position */
    int membernum;		/* if group, enumerated member */
    int goparallel;		/* true/false flag */
    int childtasks;		/* count of parallel child tasks */
    int nprocs;			/* 4 * (number of processors - 1) */
};
 
 
static Object Object_Slice(Object o, struct argblk *a);
static Object Field_Slice(Field o, struct argblk *a);
static Array Array_Slice(Array a, int dim, int pos, int num, int *counts, 
			 int flag);
static Error Slice_Wrapper(Pointer a);
static Error Slice_Wrapper2(Pointer a);
static Error Slice_Wrapper3(Pointer a);

static Error todouble(Type t, Pointer v1, int i1, double *v2);
static Error fromdouble(Type t, double v1, Pointer v2, int i2);

static Error numslices(CompositeField f, struct argblk *a);
static Error fixtype(Series s);
extern void _dxfPermute(int dim, Pointer out, int *ostrides, int *counts, 
		    int size, Pointer in, int *istrides);

/* the source for these is in partreg.c */
extern Field _dxf_has_ref_lists(Field f, int *pos, int *con);
extern Error _dxf_make_map_template(Array a, Array *map);
extern Error _dxf_fix_map_template(Array map, Array newmap);
extern Array _dxf_remap_by_map(Array a, Array map);
 
#define IsEmpty(o)  (DXGetObjectClass(o)==CLASS_FIELD && DXEmptyField((Field)o))

 
#define MAXSHAPE   128        /* XXX - system constant */
#define MAXDIM     8          /* XXX - system constant */


/* 0 = input
 * 1 = dim number
 * 2 = position(s)
 *
 * 0 = output
 */
int
m_Slice(Object *in, Object *out)
{
    char *cp, cdim;
    struct argblk a;

 
    /* initialize things which need it - especially anything which may be 
     *  allocated - it makes the error exit code simpler.  
     */
    a.compositeseries = NULL;		/* object */
    a.compositepositions = NULL;	/* memory */
    a.sliceplace = NULL;		/* memory */
    a.compositecount = 0;		/* integer */
    a.outgroup = NULL;			/* object */
    a.childtasks = 0;			/* integers */
    a.goparallel = 1;
    a.childtasks = 0;
    a.nprocs = (DXProcessors(0) - 1) * 4;

    /* input object */
    if (!in[0])  {
	DXSetError(ERROR_BAD_PARAMETER, "#10000", "input");
	return ERROR;
    }
 
    /* dimension to slice along.  defaults to 0 */
    if (in[1]) {
        /* accept either an integer or the letters x,y or z */
	if (!DXExtractInteger(in[1], &a.slicedim)) {
            if (DXExtractString(in[1], &cp)) {
                cdim = isupper(cp[0]) ? tolower(cp[0]) : cp[0];
                if (cdim < 'x' || cdim > 'z') {
                    DXSetError(ERROR_BAD_PARAMETER, "bad dimension letter");
		    return ERROR;
		}
                a.slicedim = (int)(cdim - 'x');
            } else {
		DXSetError(ERROR_BAD_PARAMETER, "#10620", "dimension number");
		return ERROR;
	    }
        }
        if (a.slicedim < 0) {
            DXSetError(ERROR_BAD_PARAMETER, "#10020", "dimension number");
	    return ERROR;
	}
    } else
	a.slicedim = 0;
    
    
    /* where to make the slice(s).  either a single connection number
     *  or a list of connection numbers along the dimension where the
     *  slices are to be taken.  the default is to pass a count of 0
     *  and a NULL pointer into the worker routine, which will make a
     *  series of slices, one at each connection of the input object.
     */
    if (in[2]) {
	if (!DXQueryParameter(in[2], TYPE_INT, 1, &a.slicecount)) {
	    DXSetError(ERROR_BAD_PARAMETER, "#10050", "positions");
	    return ERROR;
	}
	
	if (a.slicecount <= 0) {
	    DXSetError(ERROR_BAD_PARAMETER, "#10050", "positions");
	    return ERROR;
	}
        
        if (!(a.sliceplace = (int *)DXAllocate(sizeof(int) * a.slicecount)))
            return ERROR;
        
        if (!DXExtractParameter(in[2], TYPE_INT, 1, 
			      a.slicecount, (Pointer)a.sliceplace)) {
            DXSetError(ERROR_BAD_PARAMETER, "#10050", "positions");
            goto error;
        }
        
    } else
	a.slicecount = 0;
    
    
    /* allow empty fields to pass thru unchanged */
    if (IsEmpty(in[0])) {
	out[0] = in[0];
	return OK;
    }

    /* see if there is at least one field in the input somewhere */
    if(!DXExists(in[0], "positions")) {
	DXSetError(ERROR_INVALID_DATA, "#10190", "input");
	return ERROR;
    }
    
    out[0] = Object_Slice(in[0], &a);

    /* if we want to make it so that if the slice location is outside
     *  the object, that it isn't an error, but just returns an empty
     *  object, then add this code here.
     */
#if 1    /* EMPTY_NOT_AN_ERROR */
    if (!out[0] && DXGetError()==ERROR_NONE)
        out[0] = (Object)DXEndField(DXNewField());
#endif

    if (!out[0])
        goto error;

    DXFree((Pointer)a.compositepositions);
    DXFree((Pointer)a.sliceplace);
    return OK;

  error:
    DXDelete((Object)a.compositeseries);
    DXFree((Pointer)a.compositepositions);
    DXFree((Pointer)a.sliceplace);
    DXDelete((Object)a.outgroup);

    return ERROR;
}



/* recursive worker routine.
 *  input is an object and a filled in parameter block.
 */
static Object Object_Slice(Object o, struct argblk *a)
{
    Object newo, newo2;
    Object subo, subo2;
    char *name;
    Matrix m;
    float position;
    int fixed, z;
    int i, newi, numparts = 0;
    int goparallel = a->goparallel;

    if (a->childtasks >= a->nprocs)
	goparallel = 0;
 
    switch(DXGetObjectClass(o)) {
 
      case CLASS_GROUP:
	switch (DXGetGroupClass((Group)o)) {
	  case CLASS_COMPOSITEFIELD:   /* partitioned data */

            /* three different cases:
             *  number of slices == 0 (create slices at all locations)
             *  number of slices == 1 (create a single slice)
             *  number of slices > 1 (multiple slices at user selected locations)
             */


	    if (!numslices((CompositeField)o, a))
		return NULL;


            /* compute the number of connections in the composite field,
             *  and their actual location.  input is an object divided
             *  into fields based on the locations of the positions of the
             *  points.  output is a series object, where each field is a
             *  partitioned slice.  the series position is the original 
             *  location of the origin of the slice.  eventually this will
             *  run in parallel, so create a series group with initially empty
             *  composite fields to hold the slices here.
             */

            if (a->slicecount == 1)
                goto normalgroup;


            a->compositeseries = (Object)DXNewSeries();
            if (!a->compositeseries)
                return NULL;

            /* for each new slice */
	    for (i=0; i < a->slicecount; i++) {
                newo = DXCopy(o, COPY_ATTRIBUTES);
                if (!newo)
                    return NULL;

                if (!DXSetSeriesMember((Series)a->compositeseries, i,
                                     (double)a->compositepositions[i], newo)) {
                    DXDelete(newo);
                    return NULL;
                }
            }
	    
	    if (goparallel) {
		if (!DXCreateTaskGroup())
		    return NULL;

		/* add the total number of subtasks which will be added to the
		 *  task queue here, so when the child tasks decide whether
		 *  they can go parallel, they know how many previous tasks
		 *  have been started.
		 */
		DXGetMemberCount((Group)o, &numparts);
		a->childtasks += numparts;
	    }


            /* for each original partition */
            for (i=0; subo = DXGetEnumeratedMember((Group)o, i, &name); i++) {
                if (DXEmptyField((Field)subo))
                    continue;
		
		if (goparallel) {
		    a->toslice = subo;
		    if (!DXAddTask(Slice_Wrapper, (Pointer)a, 
				 sizeof(struct argblk), 1.0)) {
			DXAbortTaskGroup();
			return NULL;
		    }

		} else if (!Field_Slice((Field)subo, a))
		    return NULL;
		    
	    }

	    if (goparallel) {
		if (!DXExecuteTaskGroup())
		    return NULL;

		a->childtasks -= numparts;
	    }
		
	    /* and now bubble the data type up from the first non-empty
             *  composite field and set the series group type.
	     */
	    if (!fixtype((Series)a->compositeseries))
		return NULL;

            o = a->compositeseries;
	    break;
 
	  default: 	    /* normal groups */

          normalgroup:
            newo2 = DXCopy(o, COPY_ATTRIBUTES);
            if (!newo2)
                return NULL;

            /* for each member */
	    if (goparallel) {
		a->outgroup = newo2;
		if (!DXCreateTaskGroup())
		    return NULL;

		/* add the total number of subtasks which will be added to the
		 *  task queue here, so when the child tasks decide whether
		 *  they can go parallel, they know how many previous tasks
		 *  have been started.
		 */
		DXGetMemberCount((Group)o, &numparts);
		a->childtasks += numparts;
	    }

            newi = 0;
	    for (i=0; subo = DXGetEnumeratedMember((Group)o, i, &name); i++) {
		if (goparallel) {
		    a->toslice = subo;
		    a->membername = name;
		    
		    if (!DXAddTask(Slice_Wrapper2, (Pointer)a, 
				 sizeof(struct argblk), 1.0)) {
			DXAbortTaskGroup();
			return NULL;
		    }

		} else {
		    if (!(newo = Object_Slice(subo, a))) {
			if (DXGetError()==ERROR_NONE)
			    continue;
			
			DXDelete(newo2);
			return NULL;
		    }
		    
		    /* add new object to copy of input, with the same name if
		     *  input is named - else add as an unnamed member.
		     */
		    name ? DXSetMember((Group)newo2, name, newo) :
			   DXSetEnumeratedMember((Group)newo2, newi, newo);

                newi++;
		}
	    }

	    if (goparallel) {
		if (!DXExecuteTaskGroup())
		    return NULL;
		a->childtasks -= numparts;
	    }
            o = newo2;
	    break;
 
	  case CLASS_SERIES:
            newo2 = DXCopy(o, COPY_HEADER);
            if (!newo2)
                return NULL;

            /* for each member */
	    if (goparallel) {
		a->outgroup = newo2;
		if (!DXCreateTaskGroup())
		    return NULL;

		/* add the total number of subtasks which will be added to the
		 *  task queue here, so when the child tasks decide whether
		 *  they can go parallel, they know how many previous tasks
		 *  have been started.
		 */
		DXGetMemberCount((Group)o, &numparts);
		a->childtasks += numparts;
	    }

            /* for each member */
	    for (i=0; subo = DXGetSeriesMember((Series)o, i, &position); i++) {
		if (goparallel) {
		    a->toslice = subo;
		    a->seriesplace = position;
		    a->membernum = i;
		    
		    if (!DXAddTask(Slice_Wrapper3, (Pointer)a, 
				 sizeof(struct argblk), 1.0)) {
			DXAbortTaskGroup();
			return NULL;
		    }

		} else {
		    if (!(newo = Object_Slice(subo, a))) {
			if (DXGetError()==ERROR_NONE)
			    /* for a series, have to create a placeholder */
			    newo = (Object)DXEndField(DXNewField());
			else {
			    /* error code was set.  delete copy, return error.
			     */
			    DXDelete(newo2);
			    return NULL;
			}
			
		    }

		    if(!DXSetSeriesMember((Series)newo2, i, position, newo)) {
			DXDelete(newo2);
			return NULL;
		    }
		}
		
	    }
	    if (goparallel) {
		if (!DXExecuteTaskGroup())
		    return NULL;
		a->childtasks -= numparts;
	    }
	    o = newo2;
	    break;
	}
	
        /* end of group objects */
	break;
	
 
      case CLASS_FIELD:
	if (DXEmptyField((Field)o))
	    break;
	
	o = Field_Slice((Field)o, a);
	break;
 
 
      case CLASS_XFORM:
	if (!DXGetXformInfo((Xform)o, &subo, &m))
	    return NULL;
 
	if (!(newo=Object_Slice(subo, a)))
	    return NULL;
 
	o = (Object)DXNewXform(newo, m);
	break;
 
 
      case CLASS_SCREEN:
	if (!DXGetScreenInfo((Screen)o, &subo, &fixed, &z))
	    return NULL;
 
	if (!(newo=Object_Slice(subo, a)))
	    return NULL;
 
	o = (Object)DXNewScreen(newo, fixed, z);
	break;
 
 
      case CLASS_CLIPPED:
	if (!DXGetClippedInfo((Clipped)o, &subo, &subo2))
	    return NULL;
 
	if (!(newo=Object_Slice(subo, a)))
	    return NULL;
 
	o = (Object)DXNewClipped(newo, subo2);
	break;
 
	
      default:
	/* any remaining objects can't contain fields, so we can safely
	 *  pass them through unchanged - just make an extra copy.
	 */
        o = DXCopy(o, COPY_STRUCTURE);
	break;
 
    }
 
    return o;
}
 
 
/*
 */
static Error Slice_Wrapper(Pointer a)
{
    Object o;

    o = Field_Slice((Field)((struct argblk *)a)->toslice, (struct argblk *)a);
    return (o ? OK : ERROR);
}

/*
 */
static Error Slice_Wrapper2(Pointer a)
{
    struct argblk *ap = (struct argblk *)a;
    struct argblk a2;
    Object o;

    memcpy((char *)&a2, (char *)a, sizeof(struct argblk));

    o = Object_Slice(a2.toslice, &a2);

    if (!o)
	return ((DXGetError() == ERROR_NONE) ? OK : ERROR);

    if (IsEmpty(o))
	return OK;
    
    DXSetMember((Group)ap->outgroup, ap->membername, o);
    
    return OK;
}

/*
 */
static Error Slice_Wrapper3(Pointer a)
{
    struct argblk *ap = (struct argblk *)a;
    struct argblk a2;
    Object o;

    memcpy((char *)&a2, (char *)a, sizeof(struct argblk));

    o = Object_Slice(a2.toslice, &a2);

    if (!o && (DXGetError() != ERROR_NONE))
	return ERROR;

    if (!o || (IsEmpty(o)))
	o = (Object)DXEndField(DXNewField());
    
    DXSetSeriesMember((Series)ap->outgroup, ap->membernum, ap->seriesplace, o);
    
    return OK;
}

/*
 * 
 */
static Object Field_Slice(Field o, struct argblk *a)
{
    Series s = NULL;
    Object newo, subo;
    Array ac = NULL;
    Array ap = NULL;
    Array ao;
    Field f = NULL;
    char *name, *dattr, *rattr;
    int i, j, k, n;
    int justone = 0;
    int start;
    int rank = 0;
    int stride[MAXDIM];
    int counts[MAXDIM], ccounts[MAXDIM], shape[MAXDIM];
    int gridoffset[MAXDIM], ngridoffset[MAXDIM];
    float origins[MAXDIM], deltas[MAXDIM*MAXDIM];
    float worldposition, scratch[MAXDIM];
    float *fp;
    ArrayHandle ah;
    int refpos, refcon;
    Array posmap = NULL;
    Array conmap = NULL;  /* temp for now - should use algorithm */
    
 
    /* have to do positions & connections first, so we know the
     *  original shape.  then other components can be sliced.
     *  do a quick check on connections to be sure they're regular
     *  before we start.
     */
    ao = (Array)DXGetComponentValue(o, "connections");
    if (!ao) {  
	DXSetError(ERROR_MISSING_DATA, "#10240", "connections");
	return NULL;
    }
    
    if (!DXQueryGridConnections(ao, &a->numdim, ccounts)) {
	DXSetError(ERROR_BAD_TYPE, "#10610", "input");
	return NULL;
    }
    
    if (a->slicedim >= a->numdim) {
	DXSetError(ERROR_BAD_PARAMETER, "#11360", a->slicedim, "input", 
		 a->numdim-1);
	return NULL;
    }
    
    /* all fields have offsets.  if the input isn't partitioned, it still
     *  returns 0 - but if you extract a field from a composite field, it
     *  may still have offsets.  set them to zero if you didn't get a CF.
     */
    if (a->compositepositions) {
	if (!DXGetMeshOffsets((MeshArray)ao, gridoffset))
	    return NULL;
    } else
	memset((char *)gridoffset, '\0', sizeof(int)*MAXDIM);

    /* if the input isn't partitioned, set the limits from this single field.
     */
    if (a->compositecount == 0)
	a->compositecount = ccounts[a->slicedim];

    
    /* if they gave us a position list, check it for out of range values 
     */
    for (i=0; i<a->slicecount; i++) {
	if (a->sliceplace[i] < 0 || a->sliceplace[i] >= a->compositecount) {
	    if (a->slicecount > 1) {
		DXSetError(ERROR_BAD_PARAMETER, "#10041",
			 "position item", i+1, 0, a->compositecount - 1);
		return NULL;
	    } else {
		DXSetError(ERROR_BAD_PARAMETER, "#10040", 
			 "position", 0, a->compositecount - 1);
		return NULL;
	    }
	}
	
    }
    
    /* if they didn't give us a list, initialize a positions list.
     *  making a slice at every connection.
     */
    if (a->slicecount == 0) {
	
	a->slicecount = a->compositecount;
	
	if (!(a->sliceplace = (int *)DXAllocate(a->slicecount * sizeof(int))))
	    return NULL;
	
	for (i=0; i<a->slicecount; i++)
	    a->sliceplace[i] = i;
    }
 
 
    /* if the input is partitioned, then an output field has already been
     *  initialized.  if there is more than one slice, then the output
     *  will be a new series, with each slice in a separate field.  otherwise
     *  the input is a simple field and the output is a simple field - make
     *  a copy of the input a little further down and work on it.
     */
    if (!a->compositeseries) {
	if (a->slicecount != 1) {
	    if (!(s = DXNewSeries()))
		return NULL;
	} else
	    justone++;
    }
 
    /* loop here for each position */
    for (i=0; i<a->slicecount; i++) {
	
	/* check here to see if the slice is inside this partition.
	 *  if slices, each partition excludes the starting face, except
         *  for the very first partition, to prevent creating 2 slices at
         *  the same location because of internal partition boundaries, but
         *  includes the ending face. 
	 */

	/* these are start and end position number, relative to this
	 *  partition.
	 */	
	start = a->sliceplace[i] - gridoffset[a->slicedim];

	/* completely outside?
         */
	if (start < 0 ||  (start >= ccounts[a->slicedim]))
	    continue;
	
	/* don't include partition minimum.
	 */
	if (a->sliceplace[i] != 0 && start <= 0)
	    continue;


	/* if there is only 1 slice, use a copy of the existing field.  
	 *  else, create a new Field to be added to the Series.
	 */
	if (justone)
	    f = (Field)DXCopy((Object)o, COPY_STRUCTURE);
	else {
	    f = DXNewField();
	    DXCopyAttributes((Object)f, (Object)o);
	}
	if (!f)
	    goto error;
	
	/*
	 * do the connections 
	 */
	for (j=0; j<a->numdim; j++) {
	    counts[j] = ccounts[j];
	    ngridoffset[j] = gridoffset[j];
	}

	/* object goes down 1 dimension in connection shape.
         *  treat lines going to a point as a special case - points don't
	 *  have a connections component.
	 */
	if (a->numdim > 1) {
	    for (j=a->slicedim+1; j<a->numdim; j++) {
		counts[j-1] = ccounts[j];
		ngridoffset[j-1] = gridoffset[j];
	    }
	    
	    if (!(ac = DXMakeGridConnectionsV(a->numdim-1, counts))
		|| !DXSetStringAttribute((Object)ac, "ref", "positions")
		|| !DXSetMeshOffsets((MeshArray)ac, ngridoffset)) {
		DXDelete((Object)ac);
		goto error;
	    }
	}
	
	/* if the new connections type would be points, either remove the
	 *  existing connections component if just one field, or for a
	 *  series, don't add a connections component.
	 */
	if (justone)
	    if (!DXDeleteComponent(f, "connections"))
		goto error;
	
	if (a->numdim > 1) {
	    if (!DXSetComponentValue(f, "connections", (Object)ac))
		goto error;
	    
	    if (justone)
		if (!DXChangedComponentValues(f, "connections"))
		    goto error;
	}
	
    
	/* now check the positions component */
	ap = (Array)DXGetComponentValue(o, "positions");
	if (!ap) {
	    DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
	    return NULL;
	}
    
    
	/* new positions.  get the actual position value for the starting
	 *  connection number.  then slice the positions array.
	 */
 
	/* if regular, compute position.  if irregular, index into data
	 *  array to find it.
	 */
	if (DXQueryGridPositions(ap, &n, NULL, origins, deltas)) {
	    if (n != a->numdim) {
		DXSetError(ERROR_INVALID_DATA, "#11161", 
			   "connection dimensions", a->numdim,
			   "position dimensions", n);
		return NULL;
	    }
	    worldposition = origins[a->slicedim];
	    for (j=0; j<n; j++)
	        worldposition += start * deltas[n*j + a->slicedim];
	} else {
            DXGetArrayInfo(ap, NULL, NULL, NULL, &rank, shape);
	    if (rank != 1 || shape[0] != a->numdim) {
		DXSetError(ERROR_INVALID_DATA, "#11161", 
			   "connection dimensions", a->numdim,
			   "position dimensions", shape[0]);
		return NULL;
	    }
	    k = 1;
	    for (j=a->numdim-1; j>=0; j--) {
		stride[j] = k;
		k *= ccounts[j];
	    }

	    ah = DXCreateArrayHandle(ap);
	    if (!ah)
		goto error;
            fp = (float *)DXGetArrayEntry(ah, start*stride[a->slicedim], 
					  (Pointer)scratch);
	    worldposition = fp[a->slicedim];

#if 0  /* very broken, and in dx 3.1.0 */
	    worldposition = *(float *)DXGetArrayEntry(ah, stride[a->slicedim],
						      (Pointer)scratch);
#endif
	    DXFreeArrayHandle(ah);
	}
	
	/* check for 1d point */
	if (a->numdim > 1) {
	    ap = Array_Slice(ap, a->slicedim, start, a->numdim, ccounts, 1);
	    if (!ap) {
		DXDelete((Object)ac);
		goto error;
	    }
	} else {
	    ap = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 1);
	    if (!DXAddArrayData(ap, 0, 1, (Pointer)&worldposition)) {
		DXDelete((Object)ac);
		goto error;
	    }
	}
 
	
	if (!DXSetComponentValue(f, "positions", (Object)ap)
	||  (justone && !DXChangedComponentValues(f, "positions"))) {
	    DXDelete((Object)ac);
	    DXDelete((Object)ap);
	    goto error;
	}



	/* check to see if we will need to make a map for any ref components.
	 *  build a map for positions.
	 */
	
	_dxf_has_ref_lists(o, &refpos, &refcon);

	/* can't do invalid connections - can't slice connection dep data.
	 */
	if (refcon) {
	    DXSetError(ERROR_NOT_IMPLEMENTED, 
		       "cannot slice components which `ref' connections");
	    goto error;
	}

	/* make a refpos map.
	 * the map is made by making a regular array of 0-N items and then
	 * permuting them the same way we are going to permute the positions.
	 */
	if (refpos) {
	    Array tmap = NULL;
	    
	    ap = (Array)DXGetComponentValue(o, "positions");
	    if (!ap)
		goto error;
	    
	    if (!_dxf_make_map_template(ap, &posmap))
		goto error;
	    
	    tmap = Array_Slice(posmap, a->slicedim, start,
			       a->numdim, ccounts, 0);
	    
	    if (!tmap)
		goto error;
	    
	    if (!_dxf_fix_map_template(posmap, tmap)) {
		DXDelete((Object) tmap);
		goto error;
	    }
	    
	    DXDelete((Object) tmap);
	    /* posmap now captures where the positions have moved to */
	}
 
	
	/* for each of the remaining components:
	 *   skip connections & positions - we've already sliced them
	 *   if it doesn't dep on positions or connections;
         *    for a single slice (therefore a single output field), we've
	 *    already made a copy of all the components, so just continue.
	 *    for a series, each slice is a separate field, so the component
	 *    must be copied to the new field.
	 *   if it does dep on positions or connections:
	 *    slice the new component and add it to the field.  for single
	 *    slices, this replaces the old component.  for series of slices,
	 *    this adds the component for the first time.
	 */
	for (j=0; subo=DXGetEnumeratedComponentValue(o, j, &name); j++) {
	    
	    if (!strcmp("positions", name) || !strcmp("connections", name))
		continue;
	    
	    /* skip derived components */
	    if (DXGetStringAttribute(subo, "der", NULL) != NULL) {
		if (justone)
		    DXDeleteComponent(f, name);
		
		continue;
	    }
	    
	    /* if dep on positions, or ref positions, slice.
	     */
	    if (!DXGetStringAttribute(subo, "dep", &dattr))
		dattr = NULL;
	    if (!DXGetStringAttribute(subo, "ref", &rattr))
		rattr = NULL;

	    if (dattr && !strcmp("positions", dattr)) {
		
		newo = (Object)Array_Slice((Array)subo, a->slicedim, start,
					   a->numdim, ccounts, 0);
		
	    }
	    else if (dattr && !strcmp("connections", dattr)) {
		
                DXSetError(ERROR_NOT_IMPLEMENTED, 
			   "cannot slice connection-dependent %s", name);
		goto error;
		
	    }
	    else if (rattr && !strcmp("positions", rattr)) {
		
		newo = (Object) _dxf_remap_by_map((Array)subo, posmap);

	    } else if (rattr && !strcmp("connections", rattr)) {

		DXSetError(ERROR_NOT_IMPLEMENTED, 
			   "cannot slice components which `ref' connections");
		goto error;

	    } else {

		if (!justone)
		    f = DXSetComponentValue(f, name, subo);
		
		continue;
		
	    }
	    
	    if (!newo
	    || !DXSetComponentValue(f, name, newo)
	    || !DXChangedComponentValues(f, name)
	    || !DXCopyAttributes(subo, newo)) {
		DXDelete(newo);
		goto error;
	    }

	}
	
	
	if (!DXEndField(f))
	    goto error;

	/* if the output is a series of composite fields, get the matching
         *  group member and add this new field to it.
	 * else if the input is a simple field and there is more than one
         *  slice, add this field as part of a series.
         * else the input is a simple field and the output is only a single
	 *  slice - there is nothing extra to do here.
	 */
	if (a->compositeseries) {
	    subo = DXGetSeriesMember((Series)a->compositeseries, i, 
				                               &worldposition);
	    if (!subo)
		goto error;

	    if (!DXSetMember((Group)subo, NULL, (Object)f))
		goto error;

	} else if (!justone) {
	    if (!DXSetSeriesMember(s, i, worldposition, (Object)f))
		goto error;
 
	}
    }

    
    if (a->compositeseries)
	return (Object)a->compositeseries;
    else if (!justone)
	return (Object)s;
    else
	return (Object)f;
 
error:
    DXDelete((Object)s);
    DXDelete((Object)f);
    return NULL;
}
 
/* slice an array.  args are:
 *  array to slice
 *  connection dimension number along which to slice
 *  connection offset along that dimension at which to start slicing
 *  number of original dimensions in the connections
 *  number of counts along the original dimensions in the connections
 *  if flag set, slice contents of array also (e.g. for positions)
 */
static Array 
Array_Slice(Array a, int dim, int pos, int num, int *ccounts, int flag)
{
    Array na = NULL, na2 = NULL, terms[MAXDIM];
    int nitems, itemsize;
    Type t;
    Category c;
    double o_x[MAXDIM], d_x[MAXDIM*MAXDIM];
    Pointer origins = o_x, deltas = d_x;
    float *ofp, *ifp;
    int rank, shape[MAXSHAPE], counts[MAXDIM], ncounts[MAXDIM];
    int ostride[MAXDIM], nstride[MAXDIM];
    int dims[MAXDIM], tempdims[MAXDIM];
    double tvalue1, tvalue2;
    Pointer op, np;
    int i, j, k, l;
    
    
    /* check here for fully regular array and shortcut the following code */
    
    if (DXQueryGridPositions(a, &i, counts, 
			     (float *)origins, (float *)deltas)) {
	
	int temp = 0;
	int dimtemp = dim;
	
	/* dim of this array less than connections? like quads in 1 space? 
	 * expand it to get the right correspondance.  same if the slice
	 * dim is already greater than the dim of this array.
	 */
	if (i < num || dim >= i)
	    goto expanded;

	/* if you are going to slice the contents of each array item in
	 * addition to slicing the total number of items, you can't use
	 * this code.
	 */
	if (flag == 1 && i != num)
	    goto expanded;

        /* if dim of array doesn't match dim of connections, expand them
         *  to be sure you get the right ones.
         */
        if (i > num) {
	    for (j=0; j<i; j++) {
		if (counts[j] <= 1) {
		    temp++;
		    if (dimtemp >= j) 
			dimtemp++;
		} else {
		    if (counts[j] != ccounts[j-temp])
			goto expanded;
		}
	    }
	    
	    /* make sure the dimensions match before you try to slice
	     * something which might not really be the same shape.
	     */
	    if (i != num+temp)
		goto expanded;
	    
	    dim = dimtemp;
	} else {
	    /* even if the dims match, check the counts to be sure the
	     * shapes are the same.
	     */
	    for (j=0; j<i; j++)
		if (counts[j] != ccounts[j])
		    goto expanded;
	}
	
	   
	/* it should be safe here to assume this array matches the connections
	 * and we can leave it regular.
	 *
	 * remove the sliced dimension 
	 */
	for (j=dim; j<i-1; j++) {
	    counts[j] = counts[j+1];
	    ((float *)origins)[j] = ((float *)origins)[j+1];
	}
	
	if (flag == 1) {
	    int realdelta;
	    int whichdim;

	    /*
	     * if only one axis varies in this dim, cut it out.
	     * (this is new code and may not match exactly what
	     *  happens if the data is irregularized first, but
	     *  it is a requested function, so i'll try it out...)
	     */
	    realdelta = 0;
	    whichdim = -1;
	    for (j=0; j<i; j++) {
		if (((float *)deltas)[j*i + dim] != 0.0) {
		    if (realdelta > 0)
			goto expanded;
		    whichdim = j;
		    realdelta++;
		}
	    }
	    if (whichdim < 0)
		whichdim = dim;
	    for (j=0; j<i; j++) {
		if (j == dim)
		    continue;
		if(((float *)deltas)[whichdim*i + j] != 0.0)
		    goto expanded;
	    }
	    
	    ifp = ofp = (float *)deltas;
	    for(j=0; j<i; j++) {
		if(j == whichdim)
		    ofp += i;
		
		else for(k=0; k<i; k++)
		    if(k == dim) 
			ofp++;
		    else
			*ifp++ = *ofp++;
	    }
	}
	
	na = DXMakeGridPositionsV(i-1, counts, (float *)origins, (float *)deltas);
	DXCopyAttributes((Object)na, (Object)a);
	return na;
    }
    
    
    /* not a regular grid */
    
    switch(DXGetArrayClass(a)) {
	
      default:
      case CLASS_ARRAY:
	/* normal expanded array */
expanded:
	if (!DXGetArrayInfo(a, &nitems, &t, &c, &rank, shape))
	    return NULL;

	itemsize = DXGetItemSize(a);
	op = DXGetArrayData(a);

	na = DXNewArrayV(t, c, rank, shape);
	DXCopyAttributes((Object)na, (Object)a);
	
	np = DXGetArrayData(DXAddArrayData(na, 0, nitems/ccounts[dim], NULL));
	if(!np)
	    goto error;

	if (num <= 1) {
	    /* don't bother with Permute */
	    op = (Pointer)((char *)op + pos * itemsize);
	    memcpy(np, op, itemsize);
	    return na;
	}

	/* try for one general routine, and be prepared to optimize
	 * on slicedim == numdim, and itemsize == sizeof(float)
	 */
	
	k = 0;
	for(i=0; i<dim; i++)
	    dims[k++] = i;
	for(i=dim+1; i<num; i++)
	    dims[k++] = i;
	
	
	/* old strides - w/ sliced dim */
	k = 1;
	for (i=num-1; i>=0; i--) {
	    tempdims[i] = k;
	    k *= ccounts[i];
	}
	
	/* rearrange dims */
	for (i=num-2; i>=0; i--) {
	    ostride[i] = tempdims[dims[i]];	
	    ncounts[i] = ccounts[dims[i]];
	}
	
	/* new strides - w/o sliced dim */
	k = 1;
	for (i=num-1; i>dim; i--) {
	    nstride[i-1] = k;
	    k *= ncounts[i-1];
	}
	for (i=dim-1; i>=0; i--) {
	    nstride[i] = k;
	    k *= ncounts[i];
	}
	
	op = (Pointer)((char *)DXGetArrayData(a) + pos*tempdims[dim]*itemsize);
	np = (Pointer)DXGetArrayData(na);
	if (!op || !np)
	    goto error;

#if 0
	DXMarkTimeLocal("begin permute");
#endif
	_dxfPermute(num-1, np, nstride, ncounts, itemsize, op, ostride);
#if 0
	DXMarkTimeLocal("end permute");
#endif
	
	if(!flag)
	    return na;

	/* take care of shape here, because i don't think the other
	 *  cases need to worry about it.  this is inefficient, since
	 *  we are traversing the data array twice - once to cut out
	 *  the items in the sliced dimension, and once to remove the
	 *  sliced dimension from each data item.
	 */
	DXGetArrayInfo(na, &nitems, &t, &c, &rank, shape);
	if(rank != 1) {
	    DXSetError(ERROR_INVALID_DATA, 
		"array contents shape must be 1");
	    goto error;
	}
        if (shape[0] != num) {
	    DXSetError(ERROR_INVALID_DATA,
		 "array contents dimensions must match connection dimensions");
	    goto error;
	}
        if (shape[0] <= 1) {
	    DXSetError(ERROR_INVALID_DATA,
		 "can't slice scalar array contents");
	    goto error;
	}

	--shape[0];
	
	na2 = DXNewArrayV(t, c, rank, shape);
	DXCopyAttributes((Object)na2, (Object)na);
	if(!na2)
	    goto error;
	
	if (t != TYPE_FLOAT) {
	    DXSetError(ERROR_NOT_IMPLEMENTED, "only float contents");
	    goto error;
	}

	ifp = (float *)DXGetArrayData(na);
	ofp = (float *)DXGetArrayData(DXAddArrayData(na2, 0, nitems, NULL));
	if(!ofp)
	    goto error;
	
	
	/* delete a column from the array (copy, not in place) */
	for(k=0; k < nitems; k++) {
            
	    for(l=0; l<dim; l++)
		ofp[l] = ifp[l];
	    
	    for(l=dim; l<shape[0]; l++)
		ofp[l] = ifp[l+1];
	    
	    ofp += shape[0];
	    ifp += shape[0] + 1;
	}
	
	DXDelete((Object)na);
	na = na2;

	break;

	
      case CLASS_REGULARARRAY:
	/* the output slice can't still be regular unless
         *  you are slicing along the slowest varying dimension...
	 */
	if (dim != 0)
	    goto expanded;

	/* need shape[0] from this */
	DXGetArrayInfo(a, NULL, &t, NULL, NULL, shape);
	DXGetRegularArrayInfo((RegularArray)a, &nitems, origins, deltas);

        if (shape[0] != num && flag == 1) {
	    DXSetError(ERROR_INVALID_DATA,
		 "array contents dimensions must match connection dimensions");
	    return NULL;
	}

	/* compute the new origin and correct the counts 
	 */
	k = 1;
	for (i=num-1; i>0; --i)
	    k *= ccounts[i];
	k *= pos;
	
	for (j=0; j<shape[0]; j++) {
	    todouble(t, deltas, j, &tvalue1);
	    todouble(t, origins, j, &tvalue2);
	    
	    tvalue2 += tvalue1 * k;
	    
	    fromdouble(t, tvalue2, origins, j);
	}
	nitems /= ccounts[dim];

	if (flag == 0) {
	    na = (Array)DXNewRegularArray(t, shape[0], nitems, origins, deltas);
	    DXCopyAttributes((Object)na, (Object)a);
	} else {
	    if (shape[0] == 1) {
		DXSetError(ERROR_INVALID_DATA, "can't slice scalar array");
		return NULL;
	    }
	    
#define Plus(type, var) ((type *)var)[j] = ((type *)var)[j+1];

	    switch (t) {
	      case TYPE_UBYTE:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(ubyte, deltas); Plus(ubyte, origins);
		}
		break;
		
	      case TYPE_BYTE:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(byte, deltas); Plus(byte, origins);
		}
		break;
		
	      case TYPE_USHORT:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(ushort, deltas); Plus(ushort, origins);
		}
		break;

	      case TYPE_SHORT:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(short, deltas); Plus(short, origins);
		}
		break;

	      case TYPE_UINT:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(uint, deltas); Plus(uint, origins);
		}
		break;

	      case TYPE_INT:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(int, deltas); Plus(int, origins);
		}
		break;

	      case TYPE_FLOAT:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(float, deltas); Plus(float, origins);
		}
		break;

	      case TYPE_DOUBLE:
		for (j=dim; j<shape[0]-1; j++) {
		    Plus(double, deltas); Plus(double, origins);
		}
		break;

	      default:
		DXSetError(ERROR_NOT_IMPLEMENTED, "unrecognized data type");
		return NULL;
	    }

	    na = (Array)DXNewRegularArray(t, shape[0]-1, nitems, origins, 
					  deltas);
	    DXCopyAttributes((Object)na, (Object)a);
	}
	break;
 
      case CLASS_CONSTANTARRAY:
	DXGetArrayInfo(a, &nitems, &t, &c, &rank, shape);

        if ((rank != 1 || shape[0] != num) && flag == 1) {
	    DXSetError(ERROR_INVALID_DATA,
		 "array contents dimensions must match connection dimensions");
	    return NULL;
	}

	nitems /= ccounts[dim];

	if (flag == 0) {
	    na = (Array)DXNewConstantArrayV(nitems, DXGetConstantArrayData(a),
                                           t, c, rank, shape);
	    DXCopyAttributes((Object)na, (Object)a);
            break;
	} else {
	    if (rank != 1) {
		DXSetError(ERROR_INVALID_DATA, 
                    "array contents dimensions must match connection dimensions");
		return NULL;
	    }


	    memcpy(origins, DXGetConstantArrayData(a), DXGetItemSize(a));

	    /* inherit Plus macro from above */
	    switch (t) {
	      case TYPE_UBYTE:
		for (j=dim; j<shape[0]-1; j++) Plus(ubyte, origins); break;
		
	      case TYPE_BYTE:
		for (j=dim; j<shape[0]-1; j++) Plus(byte, origins); break;
		
	      case TYPE_USHORT:
		for (j=dim; j<shape[0]-1; j++) Plus(ushort, origins); break;

	      case TYPE_SHORT:
		for (j=dim; j<shape[0]-1; j++) Plus(short, origins); break;

	      case TYPE_UINT:
		for (j=dim; j<shape[0]-1; j++) Plus(uint, origins); break;

	      case TYPE_INT:
		for (j=dim; j<shape[0]-1; j++) Plus(int, origins); break;

	      case TYPE_FLOAT:
		for (j=dim; j<shape[0]-1; j++) Plus(float, origins); break;

	      case TYPE_DOUBLE:
		for (j=dim; j<shape[0]-1; j++) Plus(double, origins); break;

	      default:
		DXSetError(ERROR_NOT_IMPLEMENTED, "unrecognized data type");
		return NULL;
	    }

	    shape[0]--;

	    na = (Array)DXNewConstantArrayV(nitems, origins, t, c, 
                                            rank, shape);
	    DXCopyAttributes((Object)na, (Object)a);
	}
	break;
 
      case CLASS_PRODUCTARRAY:
	if (flag == 1)
	    goto expanded;

	/* before slicing this like the connections, make sure it is
         *  shaped exactly like the connections.  otherwise, go to the
	 *  irregular case to be sure we get the right items.
	 */
	DXGetProductArrayInfo((ProductArray)a, &j, terms);
	if (j != num || dim > j) 
	    goto expanded;
	for (i=0; i<j; i++) {
	    if (!DXGetArrayInfo(terms[i], &nitems, NULL, NULL, NULL, NULL))
		return NULL;
	    
	    if (nitems != ccounts[i] + 1)
		goto expanded;
	}

	switch(j) {
	  case 0:
	    return NULL;
	  case 1:
	    na = Array_Slice(terms[0], dim, pos, num, ccounts, flag);
	    break;
	  case 2:
	    na = terms[dim];
	    break;
	  default:
	    for (i=dim; i<num-1; i++)
		terms[i] = terms[i+1];

	    na = (Array)DXNewProductArrayV(j-1, terms);
	    DXCopyAttributes((Object)na, (Object)a);
	    break;
	}
	break;
	
	
      case CLASS_MESHARRAY:
	if (flag == 1)
	    goto expanded;

	/* see the comment above in the product array case */
	DXGetMeshArrayInfo((MeshArray)a, &j, terms);
	if (j != num || dim > j) 
	    goto expanded;
	for (i=0; i<j; i++) {
	    if (!DXGetArrayInfo(terms[i], &nitems, NULL, NULL, NULL, NULL))
		return NULL;
	    
	    if (nitems != ccounts[i] + 1)
		goto expanded;
	}

	switch(j) {
	  case 0:
	    return NULL;
	  case 1:
	    na = Array_Slice(terms[0], dim, pos, num, ccounts, flag);
	    break;
	  case 2:
	    na = terms[dim];
	    break;
	  default:
	    for (i=dim; i<num-1; i++)
		terms[i] = terms[i+1];
	    
	    na = (Array)DXNewMeshArrayV(j-1, terms);
	    DXCopyAttributes((Object)na, (Object)a);
	    break;
	}
	break;
	
      case CLASS_PATHARRAY:
	na = (Array)DXNewPathArray(1);   /* single point */
	DXCopyAttributes((Object)na, (Object)a);
	break;
	
	
    }
    
    return na;
    
  error:
    DXDelete((Object)na);
    DXDelete((Object)na2);
    return NULL;
}
 

/* if input is a CF and there will be more than 1 slice, you need to
 *  compute the positions for the entire output series now, so you can
 *  create the output series before you start creating any individual
 *  slices.  this is because you don't want to iterate through each incoming
 *  partition once per input slice - you want to compute all slices which
 *  cross this partition at the same time.  this is a bit awkward, since
 *  the partial slices each go into a different output composite field, so
 *  at the end, each outgoing CF will contain a complete single slice.
 *  all the positions are computed here, in order, and then the slices
 *  are added to the already existing CF's as they are computed.
 */
 
static Error 
numslices(CompositeField cf, struct argblk *ap)
{
    int i, j, k, l, m;
    Field f;
    Array arrc, arrp;
    int shape[MAXSHAPE]; 
    int stride[MAXDIM], counts[MAXDIM], gridoffsets[MAXDIM];
    float origins[MAXDIM], deltas[MAXDIM*MAXDIM], *pos;
    int startconn;
    int *sorted = NULL;
    struct gridlist {
	int seenoffset;
	int shape;
	int partno;
    } *gridlist = NULL;
    int beenallocated = 0;
    int seen = 0;
 
 
    /* clear these, so in case we return with an error, they are set */
    ap->compositecount = 0;
    ap->compositepositions = NULL;
 
 
    /* do some initial checks on the first partition member.
     */
    f = (Field)DXGetEnumeratedMember((Group)cf, 0, NULL);
    if (!f || DXGetObjectClass((Object)f) != CLASS_FIELD)
	DXErrorReturn(ERROR_INVALID_DATA, "compositefield member must be field");

    
    arrc = (Array)DXGetComponentValue(f, "connections");
    
    if (!DXQueryGridConnections(arrc, &ap->numdim, shape)) {
	DXSetError(ERROR_BAD_PARAMETER, "#10610", "input");
	return NULL;
    }

    if (ap->slicedim >= ap->numdim) {
	DXSetError(ERROR_BAD_PARAMETER, "#11360", ap->slicedim, "input", 
		 ap->numdim-1);
	return NULL;
    }
    
    /* before we compute the slice positions, we have to collect all the 
     *  partitions along the slice axis and compute the total number of slices.
     */
    
    for (i=0; f = (Field)DXGetEnumeratedMember((Group)cf, i, NULL); i++) {
	
	arrc = (Array)DXGetComponentValue(f, "connections");
	
	if (!DXQueryGridConnections(arrc, &ap->numdim, shape)) {
	    DXSetError(ERROR_BAD_PARAMETER, "#10610", "input");
	    return NULL;
	}
	
	if (!DXGetMeshOffsets((MeshArray)arrc, gridoffsets))
	    return NULL;
	
	/* only keep the partitions which are along the origin in the
	 *  other dimensions, so that later when the positions are
	 *  computed, if the data is deformed, the position will be
	 *  the value on the ap->slicedim axis.
	 */
	for (j=0; j<ap->numdim; j++) {
	    if (j == ap->slicedim)
		continue;
	    
	    if (gridoffsets[j] != 0)
		break;
	}
	if (j != ap->numdim)
	    continue;
	
	for (j=0; j<seen; j++)
	    if (gridoffsets[ap->slicedim] == gridlist[j].seenoffset)
		break;
	
	if (j != seen)
	    continue;
	
	if (seen >= beenallocated) {
	    beenallocated = (int)(16+seen * 1.5);
	    
	    gridlist = (struct gridlist *)DXReAllocate((Pointer)gridlist, 
				     beenallocated * sizeof(struct gridlist));
	    if (!gridlist)
		return NULL;
	    
	}
	
	gridlist[seen].seenoffset = gridoffsets[ap->slicedim];
	gridlist[seen].shape = shape[ap->slicedim];
	gridlist[seen].partno = i;
	seen++;
	/* allow for duplicated faces */
	ap->compositecount += shape[ap->slicedim] - 1;
		
    }
    /* take into account the last face which isn't duplicated.
     *  this is now a positions count.  num of elements equals positions-1.
     */
    ap->compositecount++;
    
 
    /* if the user gave an explicit list, go through the components and
     *  find the starting positions which correspond to those slices.
     */
    if (ap->slicecount >= 1) {
 
	ap->compositepositions = (float *)
	                         DXAllocateZero(ap->slicecount * sizeof(float));
	if (!ap->compositepositions)
	    goto error;
 
	for (i=0; f = (Field)DXGetEnumeratedMember((Group)cf, i, NULL); i++) {
 
	    arrc = (Array)DXGetComponentValue(f, "connections");
 
	    if (!DXQueryGridConnections(arrc, &ap->numdim, shape)) {
		DXSetError(ERROR_BAD_PARAMETER, "#10610", "input");
		goto error;
	    }
    
	    if (!DXGetMeshOffsets((MeshArray)arrc, gridoffsets))
		goto error;
	    
	    /* only keep the partitions which are along the origin in the
             *  other dimensions, so that later when the positions are
             *  computed, if the data is deformed, the position will be the
	     *  value on the ap->slicedim axis.
             */
	    for (j=0; j<ap->numdim; j++) {
		if (j == ap->slicedim)
		    continue;
		
		if (gridoffsets[j] != 0)
		    break;
	    }
	    if (j != ap->numdim)
		continue;
	    
	    /* if any slices are inside this partition, set their actual
             *  positions.
	     */
	    arrp = NULL;
	    for (j=0; j<ap->slicecount; j++) {

		startconn = ap->sliceplace[j];

		if (startconn < 0 || startconn >= ap->compositecount) {
		    DXSetError(ERROR_INVALID_DATA, "#10040", 
			           "slice position", 0, ap->compositecount-1);
		    goto error;
		}
		
		/* if this slice is outside this partition, skip it.
		 */
		if ((startconn < gridoffsets[ap->slicedim])
		|| (startconn > gridoffsets[ap->slicedim]+shape[ap->slicedim]))
		    continue;


		if (!arrp)
		    arrp = (Array)DXGetComponentValue(f, "positions");
		
		/* if regular, compute position.  if irregular, index into data
		 *  array to find it.
		 */
		if (DXQueryGridPositions(arrp, &k, NULL, origins, deltas))
		    ap->compositepositions[j] = origins[ap->slicedim] + 
			startconn * deltas[k*ap->slicedim + ap->slicedim];
		else {
		    k = 1;
		    for (l=ap->slicedim; l>=0; l--) {
			stride[l] = k;
			k *= counts[l];
		    }
		    pos = (float *)DXGetArrayData(arrp);
		    if (!pos)
			goto error;
		    ap->compositepositions[j] = 
			              pos[startconn * stride[ap->slicedim]];
		}
	    }
	}
 
	goto done;
    }
 
    
    /* if you get here, the rangelist was not set - the default is to 
     *  make a slice at each position.
     */

    if (ap->compositecount <= 0) {
	DXSetError(ERROR_BAD_PARAMETER, "#10030", "number of slices");
	goto error;
    }

    ap->slicecount = ap->compositecount;
    
    if (!(ap->sliceplace = (int *)DXAllocate(ap->slicecount * sizeof(int))))
	return NULL;
    
    for (i=0; i<ap->slicecount; i++)
	ap->sliceplace[i] = i;

    ap->compositepositions = (float *)
	                     DXAllocateZero(ap->slicecount * sizeof(float));
    if (!ap->compositepositions)
	goto error;
 

	
    /* sort the partitions so the slices/slices will be created in numeric
     *  order (if the user gave an explicit positions list, that will
     *  define the order of the output slices and this isn't necessary)
     */
    sorted = (int *)DXAllocate(seen * sizeof(int));
    if (!sorted)
	goto error;
    
    for (i=0; i < seen; i++)
	sorted[i] = i;
    
    for (i=0; i < seen-1; i++) {
	for (j=i+1; j < seen; j++)
	    if (gridlist[sorted[j]].seenoffset < 
		gridlist[sorted[i]].seenoffset) {
		k = sorted[i];
		sorted[i] = sorted[j];
		sorted[j] = k;
	    }
    }
    
    /* compute the actual positions.
     */
    k = -1;
    f = NULL;

    for (i=0; i < ap->slicecount; i++) {

	if (!f || ap->sliceplace[i] > 
	               gridoffsets[ap->slicedim] + gridlist[sorted[k]].shape) {

	    k++;    
	    f = (Field)DXGetEnumeratedMember((Group)cf, 
					   gridlist[sorted[k]].partno, NULL);
 
	    arrc = (Array)DXGetComponentValue(f, "connections");
	    if (!DXGetMeshOffsets((MeshArray)arrc, gridoffsets))
		goto error;
	    
	    arrp = (Array)DXGetComponentValue(f, "positions");
	}
	
	startconn = ap->sliceplace[i] - gridoffsets[ap->slicedim];

	/* if regular, compute position.  if irregular, index into data
	 *  array to find it.
	 */
	if (DXQueryGridPositions(arrp, &l, NULL, origins, deltas))
	    ap->compositepositions[i] = origins[ap->slicedim] + 
		startconn * deltas[l*ap->slicedim + ap->slicedim];
	else {
	    l = 1;
	    for (m=ap->slicedim; m>=0; m--) {
		stride[m] = l;
		l *= counts[m];
	    }
	    pos = (float *)DXGetArrayData(arrp);
	    if (!pos)
		goto error;
	    ap->compositepositions[i] = pos[startconn * stride[ap->slicedim]];
	}
	
    }


 done: 
    DXFree((Pointer)gridlist);
    DXFree((Pointer)sorted);
    return OK;
 
  error:
    DXFree((Pointer)gridlist);
    DXFree((Pointer)sorted);
    return ERROR;
}


/* if the result of slice is a series of composite fields, the data type
 *  information is set on the composite fields when the first member
 *  with data is added, but this information isn't propagated up to the
 *  series.   this routine looks for the first non-empty field in the
 *  series and sets the group type based on it.
 */
static Error
fixtype(Series s)
{
    int i, j;
    Object subo, subo2;
    Type t;
    Category c;
    int rank;
    int shape[MAXDIM];

    for (i=0; subo = DXGetEnumeratedMember((Group)s, i, NULL); i++) {
	if (DXGetObjectClass(subo) != CLASS_GROUP ||
	    DXGetGroupClass((Group)subo) != CLASS_COMPOSITEFIELD) {
	    DXErrorReturn(ERROR_INTERNAL, "missing composite field");
	}

	for (j=0; subo2 = DXGetEnumeratedMember((Group)subo, j, NULL); j++) {
	    if (DXGetObjectClass(subo2) != CLASS_FIELD || 
		DXEmptyField((Field)subo2))
		continue;
	    
	    /* no data? */
	    if (!DXGetType(subo2, &t, &c, &rank, shape))
		continue;

	    if (!DXSetGroupTypeV((Group)s, t, c, rank, shape))
		return ERROR;

	    return OK;
	}
    }

    /* if there are no fields in the series, i guess it's ok to return
     *  without setting the type - because there isn't one.
     */
    return OK;
}






#define TYPE_FROM(type, in1, out1, i1) \
{ type *p1;                   \
  p1 = (type *)out1;          \
  p1[i1] = (type) in1;      \
}

#define TYPE_TO(type, in1, i1, out1) \
{ type *p1;		      \
  p1 = (type *)in1;           \
  *out1 = (double) p1[i1];   \
}

#define TYPE_EQ(type, in1, i1, rc) \
{ type *p1;                    \
  p1 = (type *)in1;            \
  rc = (double) p1[i1] == 0.0; \
}


/* add convert a double to a value of any type
 */
static Error fromdouble(Type t, double v1, Pointer v2, int i2)
{
    switch (t) {
      case TYPE_UBYTE:  TYPE_FROM(ubyte,  v1, v2, i2);  return OK;
      case TYPE_BYTE:   TYPE_FROM(byte,   v1, v2, i2);  return OK;
      case TYPE_USHORT: TYPE_FROM(ushort, v1, v2, i2);  return OK;
      case TYPE_SHORT:  TYPE_FROM(short,  v1, v2, i2);  return OK;
      case TYPE_UINT:   TYPE_FROM(uint,   v1, v2, i2);  return OK;
      case TYPE_INT:    TYPE_FROM(int,    v1, v2, i2);  return OK;
      case TYPE_FLOAT:  TYPE_FROM(float,  v1, v2, i2);  return OK;
      case TYPE_DOUBLE: TYPE_FROM(double, v1, v2, i2);  return OK;
    }

    DXSetError(ERROR_NOT_IMPLEMENTED, "unrecognized data type");
    return ERROR;
}

/* add convert a value of any type to double.
 */
static Error todouble(Type t, Pointer v1, int i1, double *v2)
{
    switch (t) {
      case TYPE_UBYTE:  TYPE_TO(ubyte,  v1, i1, v2);  return OK;
      case TYPE_BYTE:   TYPE_TO(byte,   v1, i1, v2);  return OK;
      case TYPE_USHORT: TYPE_TO(ushort, v1, i1, v2);  return OK;
      case TYPE_SHORT:  TYPE_TO(short,  v1, i1, v2);  return OK;
      case TYPE_UINT:   TYPE_TO(uint,   v1, i1, v2);  return OK;
      case TYPE_INT:    TYPE_TO(int,    v1, i1, v2);  return OK;
      case TYPE_FLOAT:  TYPE_TO(float,  v1, i1, v2);  return OK;
      case TYPE_DOUBLE: TYPE_TO(double, v1, i1, v2);  return OK;
    }

    DXSetError(ERROR_NOT_IMPLEMENTED, "unrecognized data type");
    return ERROR;
}

/* check if value1 when cast to double equals 0.0.  1 is true, 0 is false
 */
static int my_iszero (Type t, Pointer v1, int i1)
{
    int rc; 

    switch (t) {
      case TYPE_UBYTE:  TYPE_EQ(ubyte,  v1, i1, rc);  return rc;
      case TYPE_BYTE:   TYPE_EQ(byte,   v1, i1, rc);  return rc;
      case TYPE_USHORT: TYPE_EQ(ushort, v1, i1, rc);  return rc;
      case TYPE_SHORT:  TYPE_EQ(short,  v1, i1, rc);  return rc;
      case TYPE_UINT:   TYPE_EQ(uint,   v1, i1, rc);  return rc;
      case TYPE_INT:    TYPE_EQ(int,    v1, i1, rc);  return rc;
      case TYPE_FLOAT:  TYPE_EQ(float,  v1, i1, rc);  return rc;
      case TYPE_DOUBLE: TYPE_EQ(double, v1, i1, rc);  return rc;
    }

    return 0;
}