/***********************************************************************/ /* 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" */ /***********************************************************************/ /* * $Header: /home/gda/dxcvs/dx/src/exec/dxmods/slab.c,v 1.3 1999/05/10 15:45:31 gda Exp $ */ #include #include #include #include #define argblk slabargblk /* struct for passing parameters around */ struct argblk { Object toslab; /* object to process */ int slabdim; /* dimension to make slabs along */ int slabcount; /* number of slabs to make */ int slabthickness; /* thickness (in connection numbers) */ int *slabplace; /* connection numbers at which to make slabs */ int numdim; /* dimensionality of input object */ Object compositeobj; /* if !NULL, multiple slabs 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_Slab(Object o, struct argblk *a); static Object Field_Slab(Field o, struct argblk *a); static Array Array_SlabP(Array a, int dim, int pos, int thick, int num, int *counts); static Array Array_SlabC(Array a, int dim, int pos, int thick, int num, int *counts); static Error Slab_Wrapper(Pointer a); static Error Slab_Wrapper2(Pointer a); static Error Slab_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 numslabs(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 MAXSHAPE 128 /* XXX - system constant */ #define MAXDIM 8 /* XXX - system constant */ #define IsEmpty(o) (DXGetObjectClass(o)==CLASS_FIELD && DXEmptyField((Field)o)) /* 0 = input * 1 = dim number * 2 = thickness * 3 = position(s) * * 0 = output */ int m_Slab(Object *in, Object *out) { char *cp, cdim; struct argblk a; /* initialize anything which will be allocated eventually. * it makes the error exit code simpler. */ a.compositeobj = NULL; /* object */ a.compositepositions = NULL; /* memory */ a.slabplace = 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 slab along. defaults to 0 */ if (in[1]) { /* accept either an integer or the letters x,y or z */ if (!DXExtractInteger(in[1], &a.slabdim)) { 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.slabdim = (int)(cdim - 'x'); } else { DXSetError(ERROR_BAD_PARAMETER, "#10620", "dimension number"); return ERROR; } } if (a.slabdim < 0) { DXSetError(ERROR_BAD_PARAMETER, "#10020", "dimension number"); return ERROR; } } else a.slabdim = 0; /* where to make the slab(s). either a single connection number * or a list of connection numbers along the dimension where the * slabs 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 slabs, one at each connection of the input object. */ if (in[2]) { if (!DXQueryParameter(in[2], TYPE_INT, 1, &a.slabcount)) { DXSetError(ERROR_BAD_PARAMETER, "#10050", "positions"); return ERROR; } if (a.slabcount <= 0) { DXSetError(ERROR_BAD_PARAMETER, "#10050", "positions"); return ERROR; } if (!(a.slabplace = (int *)DXAllocate(sizeof(int) * a.slabcount))) return ERROR; if (!DXExtractParameter(in[2], TYPE_INT, 1, a.slabcount, (Pointer)a.slabplace)) { DXSetError(ERROR_BAD_PARAMETER, "#10050", "positions"); goto error; } } else a.slabcount = 0; /* the thickness of each slab. defaults to 0 if user gives a position, * else defaults to 1. */ if (in[3]) { if (!DXExtractInteger(in[3], &a.slabthickness)) { DXSetError(ERROR_BAD_PARAMETER, "#10030", "thickness"); goto error; } if (a.slabthickness < 0) { DXSetError(ERROR_BAD_PARAMETER, "#10030", "thickness"); goto error; } } else a.slabthickness = in[2] ? 0 : 1; /* allow empty fields to pass thru unchanged */ if (IsEmpty(in[0])) { out[0] = in[0]; return OK; } /* make sure there is a valid field in the input somewhere */ if(!DXExists(in[0], "positions")) { DXSetError(ERROR_INVALID_DATA, "#10190", "input"); return ERROR; } /* ok, try to slab */ out[0] = Object_Slab(in[0], &a); /* if we want to make it so that if the slab 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.slabplace); return OK; error: DXDelete((Object)a.compositeobj); DXFree((Pointer)a.compositepositions); DXFree((Pointer)a.slabplace); DXDelete((Object)a.outgroup); return ERROR; } /* recursive worker routine. * input is an object and a filled in parameter block. */ static Object Object_Slab(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 slabs == 0 (slab at all locations) * number of slabs == 1 (create a single slab) * number of slabs > 1 (multiple slabs at user selected locations) */ if (!numslabs((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 slab. the series position is the original * location of the origin of the slab. eventually this will * run in parallel, so create a series group with initially empty * composite fields to hold the slabs here. */ if (a->slabcount == 1) goto normalgroup; a->compositeobj = (Object)DXNewSeries(); if (!a->compositeobj) return NULL; /* for each new slab */ for (i=0; i < a->slabcount; i++) { newo = DXCopy(o, COPY_ATTRIBUTES); if (!newo) return NULL; if (!DXSetSeriesMember((Series)a->compositeobj, 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->toslab = subo; if (!DXAddTask(Slab_Wrapper, (Pointer)a, sizeof(struct argblk), 1.0)) { DXAbortTaskGroup(); return NULL; } } else if (!Field_Slab((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->compositeobj)) return NULL; o = a->compositeobj; 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->toslab = subo; a->membername = name; if (!DXAddTask(Slab_Wrapper2, (Pointer)a, sizeof(struct argblk), 1.0)) { DXAbortTaskGroup(); return NULL; } } else { if (!(newo = Object_Slab(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; 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->toslab = subo; a->seriesplace = position; a->membernum = i; if (!DXAddTask(Slab_Wrapper3, (Pointer)a, sizeof(struct argblk), 1.0)) { DXAbortTaskGroup(); return NULL; } } else { if (!(newo = Object_Slab(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_Slab((Field)o, a); break; case CLASS_XFORM: if (!DXGetXformInfo((Xform)o, &subo, &m)) return NULL; if (!(newo=Object_Slab(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_Slab(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_Slab(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 Slab_Wrapper(Pointer a) { Object o; o = Field_Slab((Field)((struct argblk *)a)->toslab, (struct argblk *)a); return (o ? OK : ERROR); } /* */ static Error Slab_Wrapper2(Pointer a) { struct argblk *ap = (struct argblk *)a; struct argblk a2; Object o; memcpy((char *)&a2, (char *)a, sizeof(struct argblk)); o = Object_Slab(a2.toslab, &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 Slab_Wrapper3(Pointer a) { struct argblk *ap = (struct argblk *)a; struct argblk a2; Object o; memcpy((char *)&a2, (char *)a, sizeof(struct argblk)); o = Object_Slab(a2.toslab, &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_Slab(Field o, struct argblk *a) { Series s = NULL; Object newo, subo; Array ac = NULL; Array ap = NULL; Array ao; Field f = NULL; char *name, *attr; int i, j, k, n; int justone = 0; int start, slab; int stride[MAXDIM]; int counts[MAXDIM], ccounts[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 slabbed. * 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->slabdim >= a->numdim) { DXSetError(ERROR_BAD_PARAMETER, "#11360", a->slabdim, "input", a->numdim-1); return NULL; } /* all fields have offsets. if the input isn't partitioned, it still * returns 0. but if this field was extracted from a composite field, * then the mesh offsets don't make sense - set them to zero. */ 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->slabdim]; /* check the thickness against the actual number of connections */ if ((a->slabthickness+1) < 0) { DXSetError(ERROR_BAD_PARAMETER, "#10040", "thickness", 0, a->compositecount - 1); return NULL; } /* if they gave us a position list, check it for out of range values */ for (i=0; islabcount; i++) { if (a->slabplace[i] < 0 || a->slabplace[i] >= a->compositecount) { if (a->slabcount > 1) { DXSetError(ERROR_BAD_PARAMETER, "#10040", "position item %d", 0, a->compositecount - 1, i+1); return NULL; } else { DXSetError(ERROR_BAD_PARAMETER, "#10040", "position", 0, a->compositecount - 1); return NULL; } } /* check start+thick to see if it extends out of volume */ if ((a->slabplace[i]+a->slabthickness+1) > a->compositecount) { slab = a->compositecount - a->slabplace[i] - 1; if (a->slabcount > 1) DXWarning("slab %d extends out of data, thickness set to %d", i+1, slab); else DXWarning("slab extends out of data, thickness set to %d", slab); } } /* if they didn't give us a list, initialize a positions list, * making a position every N==thickness connections. */ if (a->slabcount == 0) { if (a->slabthickness > 0) a->slabcount = (a->compositecount-1 + a->slabthickness-1) / a->slabthickness; else a->slabcount = a->compositecount; if (!(a->slabplace = (int *)DXAllocate(a->slabcount * sizeof(int)))) return NULL; if (a->slabthickness > 0) for (i=0, slab=0; islabcount; i++, slab += a->slabthickness) a->slabplace[i] = slab; else for (i=0; islabcount; i++) a->slabplace[i] = i; } /* if the input is partitioned, then an output field has already been * initialized. if there is more than one slab/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->compositeobj) { if (a->slabcount != 1) { if (!(s = DXNewSeries())) return NULL; } else justone++; } /* loop here for each position */ for (i=0; islabcount; i++) { /* check here to see if the slab/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. * if slabs, each partition includes all endpoints, but slabs * cannot begin at the exact end of a partition. */ /* these are start and end position number, relative to this * partition. */ start = a->slabplace[i] - gridoffset[a->slabdim]; slab = start + a->slabthickness; /* completely outside? */ if ((start < 0 && slab < 0) || (start >= ccounts[a->slabdim] && slab >= ccounts[a->slabdim])) continue; /* if slices, don't include partition minimum. if slabs, don't * start a slab at the last position in a partition, or end a * slab at the first position in a partition. */ if (a->slabthickness == 0 && a->slabplace[i] != 0 && start <= 0 && slab <= 0) continue; if (a->slabthickness > 0 && start+1 >= ccounts[a->slabdim]) continue; if (a->slabthickness > 0 && slab <= 0) continue; /* if there is only 1 slab, 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; jnumdim; j++) { counts[j] = ccounts[j]; ngridoffset[j] = gridoffset[j]; } /* reset slab so it is now the thickness instead of the ending * position number. */ slab = a->slabthickness + 1; /* if the start of this slab is before this partition, adjust the * thickness first (start is negative, so this makes slab smaller), * and then set the start to the beginning of this partition. */ if (start < 0) { slab += start; start = 0; } /* if a slab, and thickness extends past this partition, decrease * the thickness. */ if (a->slabthickness > 0 && start+slab > ccounts[a->slabdim]) slab = ccounts[a->slabdim]-start; counts[a->slabdim] = slab; if (ngridoffset[a->slabdim] > 0) { if (ngridoffset[a->slabdim] <= a->slabplace[i]) ngridoffset[a->slabdim] = 0; else ngridoffset[a->slabdim] -= a->slabplace[i]; } /* if thickness is 0, then object goes down 1 dimension in * connection shape, but stays 3D in positions. treat lines going to * a point as a special case. points have no connections component. */ if (slab == 1) { if (a->numdim > 1) { for (j=a->slabdim+1; jnumdim; j++) { counts[j-1] = ccounts[j]; ngridoffset[j-1] = gridoffset[j]; } /* makegridconnections sets elemtype now */ if (!(ac = DXMakeGridConnectionsV(a->numdim-1, counts)) || !DXSetStringAttribute((Object)ac, "ref", "positions") || !DXSetMeshOffsets((MeshArray)ac, ngridoffset)) { DXDelete((Object)ac); goto error; } } } else { if(!(ac = DXMakeGridConnectionsV(a->numdim, 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 (slab != 1 || a->numdim != 1) { if (!DXSetComponentValue(f, "connections", (Object)ac)) goto error; /* do this so it doesn't get deleted twice on error */ ac = NULL; 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 slab the positions array. */ /* if regular, compute position. if irregular, index into data * array to find it. if the dimensionality of the positions doesn't * match the dimensionality of the connections (e.g. quads in * three space), then try to account for the missing dimension. */ if (DXQueryGridPositions(ap, &n, NULL, origins, deltas)) { k = a->slabdim; if (n > a->numdim) for (j=0; j= j) k++; worldposition = origins[k]; for (j=0; jnumdim-1; j>=0; j--) { stride[j] = k; k *= ccounts[j]; } /* borrow the counts array temporarily, to look at shape[0] */ if (!DXGetArrayInfo(ap, NULL, NULL, NULL, NULL, counts)) return NULL; k = a->slabdim; if (counts[0] > a->numdim) for (j=0; j= j) k++; #if 1 /* new */ ah = DXCreateArrayHandle(ap); if (!ah) goto error; fp = (float *)DXGetArrayEntry(ah, start*stride[k], (Pointer)scratch); worldposition = fp[k]; #if 0 /* very wrong, and in the dx 3.1.0 release. sigh */ worldposition = *(float *)DXGetArrayEntry(ah, stride[k], (Pointer)scratch); #endif DXFreeArrayHandle(ah); #else /* bad. this expands partially regular data in the input object */ worldposition = ((float *)DXGetArrayData(ap))[stride[k]]; #endif } ap = Array_SlabP(ap, a->slabdim, start, slab, a->numdim, ccounts); if (!ap) { 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. * since we know that the connections are regular (you can't get into * this code otherwise), we don't have to build a connections map, * but the positions may not be regular, so we have to build a map * for them. (the following code does make a map for the connections * anyway, as a quick hack to get the function into the system). */ _dxf_has_ref_lists(o, &refpos, &refcon); /* make a refpos map (and a refcon map for now until i get a * map routine which keeps the map compact) if needed. * 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_SlabP(posmap, a->slabdim, start, slab, a->numdim, ccounts); 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 */ } if (refcon) { Array tmap = NULL; ac = (Array)DXGetComponentValue(o, "connections"); if (!ac) goto error; if (!_dxf_make_map_template(ac, &conmap)) goto error; tmap = Array_SlabC(conmap, a->slabdim, start, slab, a->numdim, ccounts); if (!tmap) goto error; if (!_dxf_fix_map_template(conmap, tmap)) { DXDelete((Object) tmap); goto error; } DXDelete((Object) tmap); /* conmap now captures where the connections have moved to */ } /* for each of the remaining components: * skip connections & positions - we've already slabbed them * either skip or delete derived components * if it doesn't dep on positions or connections; * for a single slab (therefore a single output field), we've * already made a copy of all the components, so just continue * for a series, each slab is a separate field, so the component * must be copied to the new field. * if it does dep on positions or connections: * slab the new component and add it to the field. for single * slabs, this replaces the old component. for series of slabs, * 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; } /* try dep first */ if (DXGetStringAttribute(subo, "dep", &attr)) { if (!strcmp("positions", attr)) { newo = (Object)Array_SlabP((Array)subo, a->slabdim, start, slab, a->numdim, ccounts); } else if (!strcmp("connections", attr)) { if (slab <= 1) { DXSetError(ERROR_NOT_IMPLEMENTED, "cannot make 0 thickness slab in connection dependent data"); goto error; } newo = (Object)Array_SlabC((Array)subo, a->slabdim, start, slab, a->numdim, ccounts); } else { /* dep on something we don't understand - copy thru */ if (!justone) DXSetComponentValue(f, name, subo); continue; } } else if (DXGetStringAttribute(subo, "ref", &attr)) { if (!strcmp(attr, "positions")) { newo = (Object) _dxf_remap_by_map((Array)subo, posmap); } else if (!strcmp(attr, "connections")) { newo = (Object) _dxf_remap_by_map((Array)subo, conmap); } else { /* ref on something we don't understand - copy thru */ if (!justone) DXSetComponentValue(f, name, subo); continue; } } else { /* not dep or ref anything - pass a copy thru */ if (!justone) DXSetComponentValue(f, name, subo); continue; } if (!newo) goto error; if (!DXSetComponentValue(f, name, newo) || !DXChangedComponentValues(f, name) || !DXCopyAttributes(newo, subo)) { 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 * slab, 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->compositeobj) { subo = DXGetSeriesMember((Series)a->compositeobj, 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; } } DXDelete((Object)posmap); DXDelete((Object)conmap); if (a->compositeobj) return (Object)a->compositeobj; else if (!justone) return (Object)s; else return (Object)f; error: DXDelete((Object)posmap); DXDelete((Object)conmap); DXDelete((Object)s); DXDelete((Object)f); return NULL; } /* slab an array dep on positions. args are: * array to slab * connection dimension number along which to slab * connection offset along that dimension at which to start slabbing * number of connections to include in the slab * number of original dimensions in the connections * number of counts along the original dimensions in the connections */ static Array Array_SlabP(Array a, int dim, int pos, int thick, int num, int *ccounts) { Array na = 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; int rank, shape[MAXSHAPE], counts[MAXDIM], ncounts[MAXDIM]; int ostride[MAXDIM], nstride[MAXDIM]; double tvalue1, tvalue2; Pointer op, np; int i, j, k; /* 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 slab * dim is already greater than this array dim. */ if (i < num || dim > i) goto expanded; /* dim of this array larger than connections (like quads in 3 space). * try to see if the number of flat dimensions plus the number of * connection dimensions equals this number of dimensions. if so, * skip over the flat dimensions. make sure the remaining connections * match the connection dimensions in length - otherwise expand them. */ if (i > num) { for (j=0; j= j) dimtemp++; } else { if (counts[j] != ccounts[j-temp]) goto expanded; } } /* make sure the dimensions match before you try to slab * 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=0; i--) { ostride[i] = j; j *= ccounts[i]; } for (i=0; i=0; i--) { nstride[i] = j; j *= ncounts[i]; } nitems = j; na = DXNewArrayV(t, c, rank, shape); if (!na) return NULL; np = DXGetArrayData(DXAddArrayData(na, 0, nitems, NULL)); if (!np) { DXDelete((Object)na); return NULL; } /* call permute to move the data */ _dxfPermute(num, np, nstride, ncounts, itemsize, op, ostride); break; case CLASS_REGULARARRAY: /* the output slab can't still be regular unless * you are slabbing along the slowest varying dimension... */ if (dim != 0) goto expanded; /* need type and shape[0] from this */ if (!DXGetArrayInfo(a, NULL, &t, NULL, &rank, shape) || !DXGetRegularArrayInfo((RegularArray)a, &nitems, origins, deltas)) 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 j) goto expanded; for (i=0; i j) goto expanded; for (i=0; i=0; i--) { ostride[i] = j; j *= ccounts[i] - 1; } for (i=0; i=0; i--) { nstride[i] = j; j *= ncounts[i]; } nitems = j; na = DXNewArrayV(t, c, rank, shape); if (!na) return NULL; np = DXGetArrayData(DXAddArrayData(na, 0, nitems, NULL)); if (!np) { DXDelete((Object)na); return NULL; } /* call permute to move the data */ _dxfPermute(num, np, nstride, ncounts, itemsize, op, ostride); break; case CLASS_REGULARARRAY: /* the output slab can't still be regular unless * you are slabbing along the slowest varying dimension... */ if (dim != 0) goto expanded; /* need type and shape[0] from this */ if (!DXGetArrayInfo(a, NULL, &t, NULL, &rank, shape) || !DXGetRegularArrayInfo((RegularArray)a, &nitems, origins, deltas)) return NULL; /* compute the new origin and correct the counts */ k = 1; for (i=num-1; i>0; --i) k *= (ccounts[i]-1); k *= pos; for (j=0; j j) goto expanded; for (i=0; i j) goto expanded; for (i=0; icompositecount = 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) { DXSetError(ERROR_INVALID_DATA, "compositefield member must be field"); goto error; } arrc = (Array)DXGetComponentValue(f, "connections"); if (!DXQueryGridConnections(arrc, &ap->numdim, ccounts)) { DXSetError(ERROR_BAD_PARAMETER, "#10610", "input"); goto error; } /* 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, ccounts)) { 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->slabdim axis. */ for (j=0; jnumdim; j++) { if (j == ap->slabdim) continue; if (gridoffsets[j] != 0) break; } if (j != ap->numdim) continue; for (j=0; jslabdim] == 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) goto error; } gridlist[seen].seenoffset = gridoffsets[ap->slabdim]; gridlist[seen].ccounts = ccounts[ap->slabdim]; gridlist[seen].partno = i; seen++; /* allow for duplicated faces */ ap->compositecount += ccounts[ap->slabdim] - 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->slabcount >= 1) { ap->compositepositions = (float *) DXAllocateZero(ap->slabcount * 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, ccounts)) { 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->slabdim axis. */ for (j=0; jnumdim; j++) { if (j == ap->slabdim) 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; jslabcount; j++) { startconn = ap->slabplace[j]; endconn = ap->slabplace[j] + ap->slabthickness; if (startconn < 0 || startconn >= ap->compositecount) { DXSetError(ERROR_INVALID_DATA, "#10040", "slab position", 0, ap->compositecount-1); goto error; } /* if this slab is entirely outside this partition, skip it. */ if (((startconn < gridoffsets[ap->slabdim]) && (endconn < gridoffsets[ap->slabdim])) || ((startconn > gridoffsets[ap->slabdim]+ccounts[ap->slabdim]) && (endconn > gridoffsets[ap->slabdim]+ccounts[ap->slabdim]))) continue; if (startconn < gridoffsets[ap->slabdim]) startconn = 0; 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->slabdim] + startconn * deltas[k*ap->slabdim + ap->slabdim]; else { k = 1; for (l=ap->slabdim; l>=0; l--) { stride[l] = k; k *= ccounts[l]; } pos = (float *)DXGetArrayData(arrp); if (!pos) goto error; ap->compositepositions[j] = pos[startconn * stride[ap->slabdim]]; } } } 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 slabs"); goto error; } if (ap->slabthickness > 0) ap->slabcount = (ap->compositecount-1 + ap->slabthickness-1) / ap->slabthickness; else ap->slabcount = ap->compositecount; if (!(ap->slabplace = (int *)DXAllocate(ap->slabcount * sizeof(int)))) goto error; if (ap->slabthickness > 0) for (i=0, m=0; islabcount; i++, m += ap->slabthickness) ap->slabplace[i] = m; else for (i=0; islabcount; i++) ap->slabplace[i] = i; ap->compositepositions = (float *) DXAllocateZero(ap->slabcount * sizeof(float)); if (!ap->compositepositions) goto error; /* sort the partitions so the slabs/slices will be created in numeric * order (if the user gave an explicit positions list, that will * define the order of the output slabs 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->slabcount; i++) { if (!f || ap->slabplace[i] > gridoffsets[ap->slabdim] + gridlist[sorted[k]].ccounts) { 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->slabplace[i] - gridoffsets[ap->slabdim]; /* 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->slabdim] + startconn * deltas[l*ap->slabdim + ap->slabdim]; else { l = 1; for (m=ap->slabdim; m>=0; m--) { stride[m] = l; l *= ccounts[m]; } pos = (float *)DXGetArrayData(arrp); if (!pos) goto error; ap->compositepositions[i] = pos[startconn * stride[ap->slabdim]]; } } done: DXFree((Pointer)gridlist); DXFree((Pointer)sorted); return OK; error: DXFree((Pointer)gridlist); DXFree((Pointer)sorted); return ERROR; } /* if the result of slab 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 * group. this routine looks for the first non-empty field in the * group and sets the group type based on it. */ static Error fixtype(Series s) { int i; 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) return OK; for (i=0; subo2 = DXGetEnumeratedMember((Group)subo, i, NULL); i++) { 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, return without setting anything. */ 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]; \ } /* 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; }