/***********************************************************************/ /* 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 /* * Header: * * Log: */ #include #include /* * We will create a hash table to associate neighbor partitions. * This table associates each partition with each of its corners in * IJK space. Given this table, a partition process can have access * to its neighbors by using its corners as keys into this hash table. * * The table will be keyed by the IJK point of a partition corner * and a tag indicating which corner this is. This tag contains one bit * for each of the dimensions of the data that indicates whether the * corresponding delta has been added onto the origin to get the corner * point. Thus the origin corner of a 3-D partition is tagged 000, the * corner adjacent along the 0-th dimension is tagged 001, the corner * adjacent in the 1-st dimension is tagged 010, the opposite corner in * the (0,1) plane is 011 and so on. */ struct partHash { int base[3]; int tag; Field partition; }; /* * Structure holding the overlap lengths along an axis */ struct lohi { int lo, hi; }; typedef struct lohi *LoHi; /* * Looping structure for multi-dimensional looping */ struct loop { int length; /* iteration limit */ int inc; /* iteration count */ int srcStart; /* index of start in source along this dimension */ int srcSkip; /* source bump size for this dimension */ byte *srcPtr; /* source current pointer */ int dstStart; /* index of start in destination along this dimension */ int dstSkip; /* destination bump size for this dimension */ byte *dstPtr; /* destination current pointer */ }; static PseudoKey PartHash(Key); static int PartCmp(Key, Key); static Error AddPartition(HashTable, Field); static Error AddPartitionEntry(HashTable, int *, int, Field); static Error GetIJKBox(Field, int *, int *); static Field FindNeighbor(HashTable, int *, int); static Error DestroyPartitionHash(HashTable); static Array GrowRegularGrid(Array, LoHi); static Array GrowRegularArray(Array, int *, int, int); static Array InitGrowIrregArray(Array, int *, int, int); static Error GrowPartition(Pointer); static Field GrowPartition1(HashTable, Field, int, char **, Pointer); static Field GrowPartition2(HashTable, Field, int, char **, Pointer); static Field GrowPartition3(HashTable, Field, int, char **, Pointer); static Error AddOverlapData(Field, int *, int *, Field, int *, LoHi, LoHi, LoHi, int, char **); static Error FillEmptyOverlap(Field, int *, LoHi, int, char **, Pointer); static Error ShrinkPartition(Pointer); static Array ShrinkArray(Array, int *, int *, int *, int, int); static Array ShrinkRegularArray(RegularArray, int, int, int); static Array ShrinkIrregularArray(Array, int *, int *, int *, int, int); static Array ShrinkProductArray(ProductArray, int *, int *, int *, int, int); extern Array _dxfReRef(Array, int); static int DoIt(char *, char **); static Error GetDepRef(Array, int *, int *); extern Error _dxf_RemoveDupReferences(Field); #define DEP_ON_POSITIONS 1 #define DEP_ON_CONNECTIONS 2 #define NOT_DEP 3 #define REFS_POSITIONS 1 #define REFS_CONNECTIONS 2 #define NOT_REFS 3 Field DXQueryOriginalMeshExtents(Field, int *, int *); #define GetEArray(f,i,n) \ (Array)DXGetEnumeratedComponentValue((f), (i), (n)) #define GetAttr(f,n,a) \ DXGetComponentAttribute((f), (n), (a)) ? \ DXGetString((String)DXGetComponentAttribute((f), (n), (a))) : \ NULL; #define GetEAttr(f,i,n,a) \ DXGetEnumeratedComponentAttribute((f), (i), (n), (a))) ? \ DXGetString((String)DXGetEnumeratedComponentAttribute((f),(i),(n),(a))) : \ NULL; #define SetOrigName(buf, name) sprintf((buf), "original %s", (name)) struct regGrowTask { Field partition; HashTable hashTable; int nRings; Array compArray; Pointer fill; }; Object _dxfRegGrow(Object object, int nRings, Pointer fill, Array compArray) { int i; Group group; Field child; HashTable hashTable; struct regGrowTask task; Object hashCopy; hashCopy = NULL; hashTable = NULL; if (DXGetObjectClass(object) != CLASS_GROUP) { DXSetError(ERROR_INVALID_DATA, "_dxfRegGrow called with object of wrong class"); return NULL; } group = (Group)object; hashCopy = DXCopy(object, COPY_STRUCTURE); if (! hashCopy) goto error; hashTable = DXCreateHash(sizeof(struct partHash), PartHash, PartCmp); if (! hashTable) return NULL; i = 0; while (NULL != (child = (Field)DXGetEnumeratedMember((Group)hashCopy, i++, NULL))) if (! AddPartition(hashTable, child)) goto error; group = (Group)object; if (! DXCreateTaskGroup()) goto error; task.hashTable = hashTable; task.nRings = nRings; task.compArray = compArray; task.fill = fill; i = 0; while (NULL != (child = (Field)DXGetEnumeratedMember(group, i++, NULL))) { task.partition = (Field)child; if (! DXAddTask(GrowPartition, (Pointer)&task, sizeof(task), 1.0)) { DXAbortTaskGroup(); goto error; } } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; DXDelete(hashCopy); DestroyPartitionHash(hashTable); return (Object)object; error: DXDelete(hashCopy); DestroyPartitionHash(hashTable); return NULL; } struct regShrinkTask { Field partition; }; Object _dxfRegShrink(Object object) { int i; Group group; Field child; struct regShrinkTask task; switch(DXGetObjectClass(object)) { case CLASS_FIELD: if (ShrinkPartition((Pointer)&object) == OK) return object; else return NULL; case CLASS_GROUP: if (DXGetGroupClass((Group)object) == CLASS_COMPOSITEFIELD) break; default: { DXSetError(ERROR_INVALID_DATA, "_dxfRegGrow called with object of wrong class"); return NULL; } } group = (Group)object; if (! DXCreateTaskGroup()) goto error; i = 0; while (NULL != (child = (Field)DXGetEnumeratedMember(group, i++, NULL))) { task.partition = (Field)child; if (! DXAddTask(ShrinkPartition, (Pointer)&task, sizeof(task), 1.0)) { DXAbortTaskGroup(); goto error; } } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; return (Object)object; error: return NULL; } static Error DestroyPartitionHash(HashTable ht) { struct partHash *ph; if (ht == NULL) return OK; DXInitGetNextHashElement(ht); while (NULL != (ph = (struct partHash *)DXGetNextHashElement(ht))) DXDelete((Object)(ph->partition)); DXDestroyHash(ht); return OK; } static Error GrowPartition(Pointer ptr) { HashTable ht; Field partition; int nRings; Array array; Object attr; char *str; char *components[256]; int i, j, doit; char *name, origName[256]; Pointer fill; Array compArray; ht = ((struct regGrowTask *)ptr)->hashTable; partition = ((struct regGrowTask *)ptr)->partition; nRings = ((struct regGrowTask *)ptr)->nRings; compArray = ((struct regGrowTask *)ptr)->compArray; fill = ((struct regGrowTask *)ptr)->fill; if (compArray) { char *s; int i, nComponents, len; DXGetArrayInfo(compArray, &nComponents, NULL, NULL, NULL, &len); s = (char *)DXGetArrayData(compArray); for (i = 0; i < nComponents; i++) { components[i] = s; s += len; } components[i] = NULL; } else components[0] = NULL; if (DXEmptyField(partition)) return OK; attr = DXGetComponentAttribute(partition, "connections", "element type"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "no element type specified"); goto error; } str = DXGetString((String)attr); if (! str) { DXSetError(ERROR_INVALID_DATA, "bad element type specification"); goto error; } if (! strcmp(str, "lines")) { if (! GrowPartition1(ht, partition, nRings, components, fill)) goto error; } else if (! strcmp(str, "quads")) { if (! GrowPartition2(ht, partition, nRings, components, fill)) goto error; } else if (! strcmp(str, "cubes")) { if (! GrowPartition3(ht, partition, nRings, components, fill)) goto error; } else { DXSetError(ERROR_NOT_IMPLEMENTED, "Cannot grow regular connections of element type %s", str); goto error; } if (! _dxf_RemoveDupReferences(partition)) goto error; if (! DXEndField(partition)) goto error; return OK; error: return ERROR; } static int cornerIndices1[3] = {0, 0, 4}; static int tags1[3] = {4, 0, 0}; static Field GrowPartition1(HashTable ht, Field field, int nRings, char **components, Pointer fill) { int nCount, oCount; Array nArray, oArray; int i, j; int box[8][3]; char *name; Field neighbors[3]; struct lohi growth, overlap, extra; int cornerIndex, tag; int *corner; int meshOrigin; int index[1]; if (DXEmptyField(field)) return field; if (! GetIJKBox(field, (int *)box, NULL)) return NULL; oArray = (Array)DXGetComponentValue(field, "connections"); if (! DXQueryGridConnections(oArray, NULL, &oCount)) { DXSetError(ERROR_INTERNAL, "regular connections expected"); return NULL; } /* * Get pointers to neighbors */ overlap.lo = overlap.hi = 0; for (i = 0; i < 3; i++) { cornerIndex = cornerIndices1[i]; corner = (int *)(box + cornerIndex); tag = tags1[i]; neighbors[i] = FindNeighbor(ht, corner, tag); if (neighbors[i]) { Array cA; int count; cA = (Array)DXGetComponentValue(neighbors[i], "connections"); if (! cA) return ERROR; if (! DXQueryGridConnections(cA, NULL, &count)) return ERROR; if (i == 0) { if (count-1 < nRings) overlap.lo = count-1; else overlap.lo = nRings; } else if (i == 2) { if (count-1 < nRings) overlap.hi = count-1; else overlap.hi = nRings; } } } /* * If fill is not NULL, then we generate complete growth area. * Otherwise, only generate overlap area where it overlaps a neighbor. */ if (fill) { extra.lo = nRings - overlap.lo; extra.hi = nRings - overlap.hi; growth.lo = nRings; growth.hi = nRings; } else { extra.lo = 0; extra.hi = 0; growth.lo = overlap.lo; growth.hi = overlap.hi; } nCount = oCount + (growth.lo + growth.hi); meshOrigin = box[0][0] - growth.lo; index[0] = 1; if (! AddOverlapData(field, &nCount, &meshOrigin, neighbors[1], index, &growth, &overlap, &extra, 1, components)) goto GrowPartition1_error; neighbors[1] = NULL; /* * Now for each neighboring partition, add in the irregular data from that * partition */ for (i = 0; i < 3; i++) { if (neighbors[i]) if (! AddOverlapData(field, &nCount, &meshOrigin, neighbors[i], &i, &growth, &overlap, &extra, 1, components)) goto GrowPartition1_error; } /* * If fill was requested, now's the time. */ if (fill != NULL && fill != GROW_NOFILL) { if (extra.lo) { i = 0; FillEmptyOverlap(field, &i, &extra, 1, components, fill); } if (extra.hi) { i = 2; FillEmptyOverlap(field, &i, &extra, 1, components, fill); } } DXDeleteComponent(field, "data statistics"); return field; GrowPartition1_error: return NULL; } static int cornerIndices2[3][3] = { {0, 0, 2}, {0, 0, 2}, {4, 4, 6} }; static int tags2[3][3] = { {6, 4, 4}, {2, 0, 0}, {2, 0, 0} }; static Field GrowPartition2(HashTable ht, Field field, int nRings, char **components, Pointer fill) { int nCounts[2], oCounts[2]; int indices[2]; Array nArray, oArray; int i, j; int box[8][3]; char *name; Field neighbors[3][3]; struct lohi growth[2], overlap[2], extra[2]; int cornerIndex, tag; int *corner; int meshOrigins[2]; if (DXEmptyField(field)) return field; if (! GetIJKBox(field, (int *)box, NULL)) return NULL; oArray = (Array)DXGetComponentValue(field, "connections"); if (! DXQueryGridConnections(oArray, NULL, oCounts)) { DXSetError(ERROR_INTERNAL, "regular connections expected"); return NULL; } /* * Get pointers to neighbors */ overlap[0].lo = overlap[0].hi = 0; overlap[1].lo = overlap[1].hi = 0; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) { cornerIndex = cornerIndices2[(i)][(j)]; corner = (int *)(box + cornerIndex); tag = tags2[(i)][(j)]; neighbors[i][j] = FindNeighbor(ht, corner, tag); if (neighbors[i][j]) { Array cA; int counts[2]; cA = (Array)DXGetComponentValue(neighbors[i][j], "connections"); if (! cA) return ERROR; if (! DXQueryGridConnections(cA, NULL, counts)) return ERROR; if (i == 0) { if (counts[0]-1 < nRings) overlap[0].lo = counts[0]-1; else overlap[0].lo = nRings; } else if (i == 2) { if (counts[0]-1 < nRings) overlap[0].hi = counts[0]-1; else overlap[0].hi = nRings; } if (j == 0) { if (counts[1]-1 < nRings) overlap[1].lo = counts[1]-1; else overlap[1].lo = nRings; } else if (j == 2) { if (counts[1]-1 < nRings) overlap[1].hi = counts[1]-1; else overlap[1].hi = nRings; } } } /* * If fill is not NULL, then we generate complete growth area. * Otherwise, only generate overlap area where it overlaps a neighbor. */ if (fill) { for (i = 0; i < 2; i++) { extra[i].lo = nRings - overlap[i].lo; extra[i].hi = nRings - overlap[i].hi; growth[i].lo = nRings; growth[i].hi = nRings; } } else { for (i = 0; i < 2; i++) { extra[i].lo = 0; extra[i].hi = 0; growth[i].lo = overlap[i].lo; growth[i].hi = overlap[i].hi; } } for (i = 0; i < 2; i++) { nCounts[i] = oCounts[i] + (growth[i].lo + growth[i].hi); meshOrigins[i] = box[i][i] - growth[i].lo; } indices[0] = indices[1] = 1; if (! AddOverlapData(field, nCounts, meshOrigins, neighbors[1][1], indices, growth, overlap, extra, 2, components)) goto GrowPartition2_error; neighbors[1][1] = NULL; /* * Now for each neighboring partition, add in the irregular data from that * partition */ for (i = 0; i < 3; i++) { indices[0] = i; for (j = 0; j < 3; j++) { indices[1] = j; if ((i != 1 || j != 1) && neighbors[i][j]) if (! AddOverlapData(field, nCounts, meshOrigins, neighbors[i][j], indices, growth, overlap, extra, 2, components)) goto GrowPartition2_error; } } /* * If fill was requested, now's the time. */ if (fill != NULL && fill != GROW_NOFILL) { if (extra[0].lo) { indices[0] = 0; indices[1] = 1; FillEmptyOverlap(field, indices, extra, 2, components, fill); } if (extra[0].hi) { indices[0] = 2; indices[1] = 1; FillEmptyOverlap(field, indices, extra, 2, components, fill); } if (extra[1].lo) { indices[0] = 1; indices[1] = 0; FillEmptyOverlap(field, indices, extra, 2, components, fill); } if (extra[1].hi) { indices[0] = 1; indices[1] = 2; FillEmptyOverlap(field, indices, extra, 2, components, fill); } } DXDeleteComponent(field, "data statistics"); return field; GrowPartition2_error: return NULL; } static int cornerIndices3[3][3][3] = { { {0, 0, 1}, {0, 0, 1}, {2, 2, 3} }, { {0, 0, 1}, {0, 0, 1}, {2, 2, 3} }, { {4, 4, 5}, {4, 4, 5}, {6, 6, 7} } }; static int tags3[3][3][3] = { { {7, 6, 6}, {5, 4, 4}, {5, 4, 4} }, { {3, 2, 2}, {1, 0, 0}, {1, 0, 0} }, { {3, 2, 2}, {1, 0, 0}, {1, 0, 0} } }; static Field GrowPartition3(HashTable ht, Field field, int nRings, char **components, Pointer fill) { int nCounts[3], oCounts[3]; int indices[3]; Array nArray, oArray; int i, j, k; int box[8][3]; char *name; Field neighbors[3][3][3]; struct lohi growth[3], overlap[3], extra[3]; int cornerIndex, tag; int *corner; int meshOrigins[3]; if (DXEmptyField(field)) return field; if (! GetIJKBox(field, (int *)box, NULL)) return NULL; oArray = (Array)DXGetComponentValue(field, "connections"); if (! DXQueryGridConnections(oArray, NULL, oCounts)) { DXSetError(ERROR_INTERNAL, "regular connections expected"); return NULL; } /* * Get pointers to neighbors */ overlap[0].lo = overlap[0].hi = 0; overlap[1].lo = overlap[1].hi = 0; overlap[2].lo = overlap[2].hi = 0; for (i = 0; i < 3; i++) for (j = 0; j < 3; j++) for (k = 0; k < 3; k++) { cornerIndex = cornerIndices3[i][j][k]; corner = (int *)(box + cornerIndex); tag = tags3[i][j][k]; neighbors[i][j][k] = FindNeighbor(ht, corner, tag); if (neighbors[i][j][k]) { Array cA; int counts[3]; cA = (Array)DXGetComponentValue(neighbors[i][j][k], "connections"); if (! cA) return ERROR; if (! DXQueryGridConnections(cA, NULL, counts)) return ERROR; if (i == 0) { if (counts[0]-1 < nRings) overlap[0].lo = counts[0]-1; else overlap[0].lo = nRings; } else if (i == 2) { if (counts[0]-1 < nRings) overlap[0].hi = counts[0]-1; else overlap[0].hi = nRings; } if (j == 0) { if (counts[1]-1 < nRings) overlap[1].lo = counts[1]-1; else overlap[1].lo = nRings; } else if (j == 2) { if (counts[1]-1 < nRings) overlap[1].hi = counts[1]-1; else overlap[1].hi = nRings; } if (k == 0) { if (counts[2]-1 < nRings) overlap[2].lo = counts[2]-1; else overlap[2].lo = nRings; } else if (k == 2) { if (counts[2]-1 < nRings) overlap[2].hi = counts[2]-1; else overlap[2].hi = nRings; } } } /* * If fill is not NULL, then we generate complete growth area. * Otherwise, only generate overlap area where it overlaps a neighbor. */ if (fill) { for (i = 0; i < 3; i++) { extra[i].lo = nRings - overlap[i].lo; extra[i].hi = nRings - overlap[i].hi; growth[i].lo = nRings; growth[i].hi = nRings; } } else { for (i = 0; i < 3; i++) { extra[i].lo = 0; extra[i].hi = 0; growth[i].lo = overlap[i].lo; growth[i].hi = overlap[i].hi; } } for (i = 0; i < 3; i++) { nCounts[i] = oCounts[i] + (growth[i].lo + growth[i].hi); meshOrigins[i] = box[i][i] - growth[i].lo; indices[i] = 1; } if (! AddOverlapData(field, nCounts, meshOrigins, neighbors[1][1][1], indices, growth, overlap, extra, 3, components)) goto GrowPartition3_error; neighbors[1][1][1] = NULL; /* * Now for each neighboring partition, add in the irregular data from that * partition */ for (i = 0; i < 3; i++) { indices[0] = i; for (j = 0; j < 3; j++) { indices[1] = j; for (k = 0; k < 3; k++) { indices[2] = k; if (neighbors[i][j][k]) if (! AddOverlapData(field, nCounts, meshOrigins, neighbors[i][j][k], indices, growth, overlap, extra, 3, components)) goto GrowPartition3_error; } } } /* * If fill was requested, now's the time. */ if (fill != NULL && fill != GROW_NOFILL) { if (extra[0].lo) { indices[0] = 0; indices[1] = 1; indices[2] = 1; FillEmptyOverlap(field, indices, extra, 3, components, fill); } if (extra[0].hi) { indices[0] = 2; indices[1] = 1; indices[2] = 1; FillEmptyOverlap(field, indices, extra, 3, components, fill); } if (extra[1].lo) { indices[0] = 1; indices[1] = 0; indices[2] = 1; FillEmptyOverlap(field, indices, extra, 3, components, fill); } if (extra[1].hi) { indices[0] = 1; indices[1] = 2; indices[2] = 1; FillEmptyOverlap(field, indices, extra, 3, components, fill); } if (extra[2].lo) { indices[0] = 1; indices[1] = 1; indices[2] = 0; FillEmptyOverlap(field, indices, extra, 3, components, fill); } if (extra[2].hi) { indices[0] = 1; indices[1] = 1; indices[2] = 2; FillEmptyOverlap(field, indices, extra, 3, components, fill); } } DXDeleteComponent(field, "data statistics"); return field; GrowPartition3_error: return NULL; } #define INDICES(ref, indices, strides, ndim) \ { \ int _i, _j = (ref); \ for (_i = 0; _i < (ndim)-1; _i++) { \ (indices)[_i] = _j / (strides)[_i]; \ _j = _j % (strides)[_i]; \ } \ (indices)[_i] = _j; \ } #define REFERENCE(indices, ref, strides, ndim) \ { \ int _i, _j = 0; \ for (_i = 0; _i < (ndim)-1; _i++) \ _j += (indices)[_i]*(strides)[_i]; \ (ref) = _j + (indices)[(ndim)-1]; \ } #define OFFSET(indices, offsets, ndim) \ { \ int _i; \ for (_i = 0; _i < (ndim); _i++) \ (indices)[_i] += (offsets)[_i]; \ } #define IS_IN(in, indices, min, max, ndim) \ { \ int _i; \ for ((in) = 1, _i = 0; (in) && _i < (ndim); _i++) \ (in) = ((indices)[_i] >= (min)[_i]) && \ ((indices)[_i] <= (max)[_i]); \ } static Error ShrinkPartition(Pointer ptr) { Field partition; Array array, oArray, nArray = NULL, rArray = NULL; int rOffsets[8], oOffsets[8], pSizes[8], gSizes[8]; int oLength, length; int nPositions, nConnections; char *name; int i, nDim; char origName[128]; char *attr; int scstrides[3], spstrides[3], *sstrides; int dcstrides[3], dpstrides[3], *dstrides; int pmaxi[3], cmaxi[3], *maxi; int dep, ref; char *toShrink[256], *toReplace[256], nS, nR; partition = ((struct regShrinkTask *)ptr)->partition; if (DXEmptyField(partition)) return OK; array = (Array)DXGetComponentValue(partition, "connections"); if (! array) { DXSetError(ERROR_MISSING_DATA, "connections missing"); goto error; } if (! DXQueryGridConnections(array, &nDim, gSizes)) { DXSetError(ERROR_MISSING_DATA, "connections should be regular"); goto error; } if (! DXQueryOriginalMeshExtents(partition, oOffsets, pSizes)) return OK; if (! DXQueryOriginalSizes(partition, &nPositions, &nConnections)) return NULL; scstrides[nDim-1] = 1; spstrides[nDim-1] = 1; dcstrides[nDim-1] = 1; dpstrides[nDim-1] = 1; for (i = nDim-2; i >= 0; i --) { scstrides[i] = scstrides[i+1]*(gSizes[i+1]-1); spstrides[i] = spstrides[i+1]*gSizes[i+1]; dcstrides[i] = dcstrides[i+1]*(pSizes[i+1]-1); dpstrides[i] = dpstrides[i+1]*pSizes[i+1]; } for (i = 0; i < nDim; i++) { cmaxi[i] = pSizes[i] - 1; pmaxi[i] = pSizes[i]; rOffsets[i] = -oOffsets[i]; } /* * Go through the components of the input field. For each component * determine if we are to shrink it, replace it with the pre-grown * original, or ignore it altogether. */ nS = nR = 0; for (i = 0; NULL != (array = GetEArray(partition, i, &name)); i++) if (strncmp(name, "original", 8)) { SetOrigName(origName, name); if (DXGetComponentValue(partition, origName)) { toReplace[nR++] = name; } else if (DXGetAttribute((Object)array, "ref") || DXGetAttribute((Object)array, "dep")) { toShrink[nS++] = name; } } /* * Replace the ones that have "originals" - IF the "original" has * non-zero entries. */ for (i = 0; i < nR; i++) { int nItems; SetOrigName(origName, toReplace[i]); array = (Array)DXGetComponentValue(partition, origName); if (! array) goto error; DXGetArrayInfo(array, &nItems, NULL, NULL, NULL, NULL); if (nItems == 0) array = NULL; DXSetComponentAttribute(partition, toReplace[i], "ref", NULL); DXSetComponentAttribute(partition, toReplace[i], "dep", NULL); if (! DXSetComponentValue(partition, toReplace[i], (Object)array)) goto error; DXDeleteComponent(partition, origName); } /* * Now shrink the one that require it. */ for (i = 0; i < nS; i++) { name = toShrink[i]; oArray = (Array)DXGetComponentValue(partition, name); if (! oArray) { DXSetError(ERROR_INTERNAL, "component not found"); goto error; } nArray = (Array)DXReference((Object)oArray); attr = GetAttr(partition, name, "dep"); if (! attr) { dep = NOT_DEP; } else { if (! strcmp(attr, "positions")) { oLength = nPositions; dep = DEP_ON_POSITIONS; } else if (! strcmp(attr, "connections")) { oLength = nConnections; dep = DEP_ON_CONNECTIONS; } else { DXSetError(ERROR_INVALID_DATA, "illegal dependency: %s", DXGetString((String)attr)); goto error; } } attr = GetAttr(partition, name, "ref"); if (! attr) { ref = NOT_REFS; } else { if (! strcmp(attr, "positions")) { sstrides = spstrides; dstrides = dpstrides; maxi = pmaxi; ref = REFS_POSITIONS; } else if (! strcmp(attr, "connections")) { sstrides = scstrides; dstrides = dcstrides; maxi = cmaxi; ref = REFS_CONNECTIONS; } else ref = NOT_REFS; } DXGetArrayInfo(oArray, &length, NULL, NULL, NULL, NULL); /* * If the lengths mismatch, truncate the array. Note that this * creates a new array. */ if (dep != NOT_DEP) if (length != oLength) { Array tmp; tmp = ShrinkArray(nArray, oOffsets, pSizes, gSizes, nDim, dep); if (! tmp) goto error; DXDelete((Object)nArray); nArray = (Array)DXReference((Object)tmp); } /* * If the component references something, we need to * re-reference it to reflect to portion of the referenced * array that has been removed. */ if (ref != NOT_REFS) { int nItemsIn; Type t; Category c; int r, s[32]; int *sPtr, *dPtr; int nRefsPerElt; DXGetArrayInfo(nArray, &nItemsIn, &t, &c, &r, s); /* * Two possibilities: either the input is dep on * soneting, in which case we just renumber it, or it * doesn't, in which case we remove elements that contain * no references into the original portion of the partition. */ if (dep == DEP_ON_POSITIONS) { length = nPositions; } else if (dep == DEP_ON_CONNECTIONS) { length = nConnections; } else { /* * Count the ones that will remain. */ int i, j, in; nRefsPerElt = DXGetItemSize(nArray) / sizeof(int); sPtr = (int *)DXGetArrayData(nArray); length = 0; for (i = 0; i < nItemsIn; i++, sPtr += nRefsPerElt) { for (in = 1, j = 0; in && j < nRefsPerElt; j++) { int k, indices[3]; INDICES(*sPtr, indices, sstrides, nDim); OFFSET(indices, rOffsets, nDim); for (k = 0; in && k < nDim; k++) if (indices[k] < 0 || indices[k] >= maxi[k]) in = 0; if (in) length++; } } } /* * If the length is zero, we delete the nArray reference and * null it out. Then the SetComponentValue will zip it * out of the field. Otherwise, we shrink it. */ if (length == 0) { DXDelete((Object)nArray); nArray = NULL; } else { /* * If the new array is not the original array and the * counts match, we can avoid creating a new array. * Otherwise... */ if (nArray == oArray || length != nItemsIn) { rArray = DXNewArrayV(t, c, r, s); if (! rArray) goto error; if (! DXAddArrayData(rArray, 0, length, NULL)) goto error; } else { rArray = (Array)DXReference((Object)nArray); } sPtr = (int *)DXGetArrayData(nArray); dPtr = (int *)DXGetArrayData(rArray); for (i = 0; i < nItemsIn; i++) { int j, k, first, keep; int refs[32]; for (keep = 0, j = 0; j < nRefsPerElt; j++, sPtr++) { int indices[3], first, in; INDICES(*sPtr, indices, sstrides, nDim); OFFSET(indices, rOffsets, nDim); for (in = 1, k = 0; in && k < nDim; k++) if (indices[k] < 0 || indices[k] >= maxi[k]) in = 0; if (in) { REFERENCE(indices, refs[j], dstrides, nDim); keep = 1; } else refs[j] = -1; } if (keep) { memcpy((char *)dPtr, (char *)refs, nRefsPerElt*sizeof(int)); dPtr += nRefsPerElt; } } DXDelete((Object)nArray); nArray = (Array)DXReference((Object)rArray); } } if (! DXSetComponentValue(partition, name, (Object)nArray)) goto error; DXDelete((Object)nArray); nArray = NULL; } if (! DXEndField(partition)) goto error; return OK; error: DXDelete((Object)nArray); return ERROR; } static Array ShrinkArray(Array array, int *oOffs, int *oCnts, int *gCnts, int nDim, int dep) { Array a; /* * The possibilities are: regular, product, and irregular. */ if (DXGetArrayClass(array) == CLASS_REGULARARRAY) { a = ShrinkRegularArray((RegularArray)array, *oOffs, *oCnts, dep); } else if (DXGetArrayClass(array) == CLASS_PRODUCTARRAY) { a = ShrinkProductArray((ProductArray)array, oOffs, oCnts, gCnts, nDim, dep); } else { a = ShrinkIrregularArray(array, oOffs, oCnts, gCnts, nDim, dep); } return a; } static Array ShrinkProductArray(ProductArray array, int *oOffs, int *oCnts, int *gCnts, int nDims, int dep) { int nTerms, i; Array inTerms[32], outTerms[32]; DXGetProductArrayInfo(array, &nTerms, inTerms); /* * If we have one term per dimension, we can retain product array-hood. */ if (nTerms != nDims) { return ShrinkIrregularArray((Array)array, oOffs, oCnts, gCnts, nDims, dep); } else { for (i = 0; i < nTerms; i++) { outTerms[i] = ShrinkArray(inTerms[i], oOffs+i, oCnts+i, gCnts+i, 1, dep); if (NULL == outTerms[i]) { for (i--; i >= 0; i--) DXDelete((Object)outTerms[i]); return NULL; } } } return (Array)DXNewProductArrayV(nTerms, outTerms); } static Array ShrinkRegularArray(RegularArray array, int off, int len, int dep) { Type t; double aOrg[32], aDelta[32]; int i, shape; DXGetArrayInfo((Array)array, NULL, &t, NULL, NULL, &shape); DXGetRegularArrayInfo(array, NULL, (Pointer)aOrg, (Pointer)aDelta); if (dep == DEP_ON_CONNECTIONS) len -= 1; /* * We need to bump origin by offset times the delta. */ if (t == TYPE_FLOAT) { float *fOrg, *fDelta; int i; fOrg = (float *)aOrg; fDelta = (float *)aDelta; for (i = 0; i < shape; i++) *fOrg++ += off * *fDelta++; } else if (t == TYPE_DOUBLE) { double *dOrg, *dDelta; dOrg = (double *)aOrg; dDelta = (double *)aDelta; for (i = 0; i < shape; i++) *dOrg++ += off * *dDelta++; } return (Array)DXNewRegularArray(t, shape, len, (Pointer)aOrg, (Pointer)aDelta); } static Array ShrinkIrregularArray(Array sArray, int *offs, int *oCnts, int *gCnts, int nDim, int dep) { int i, j, k, itemSize, nItems; struct loop loop[3]; Array dArray; byte *srcData, *dstData; Type type; Category cat; int rank, shape[32]; int oK[8], gK[8]; if (dep == DEP_ON_POSITIONS) for (i = 0; i < nDim; i++) { oK[i] = oCnts[i]; gK[i] = gCnts[i]; } else for (i = 0; i < nDim; i++) { oK[i] = oCnts[i] - 1; gK[i] = gCnts[i] - 1; } dArray = NULL; nItems = 1; for (i = 0; i < nDim; i++) nItems *= oK[i]; /* * Set up generic looping. Note that we invert the looping order * here. */ for (i = 0; i < nDim; i++) { loop[i].srcStart = offs[(nDim-1)-i]; loop[i].dstStart = 0; loop[i].length = oK[(nDim-1)-i]; if (i == 0) { loop[i].srcSkip = 1; loop[i].dstSkip = 1; } else { loop[i].srcSkip = gK[nDim-i] * loop[i-1].srcSkip; loop[i].dstSkip = oK[nDim-i] * loop[i-1].dstSkip; } } DXGetArrayInfo(sArray, NULL, &type, &cat, &rank, shape); dArray = DXNewArrayV(type, cat, rank, shape); if (! dArray) return NULL; if (! DXAddArrayData(dArray, 0, nItems, NULL)) { DXDelete((Object)dArray); return NULL; } itemSize = DXGetItemSize(sArray); srcData = (byte *)DXGetArrayData(sArray); dstData = (byte *)DXGetArrayData(dArray); /* * Initialize looping. Note: dst start is always [0][0]...[0] * For the source, we need to skip the overlap data */ for (j = 0; j < nDim; j++) srcData += (loop[j].srcStart * loop[j].srcSkip * itemSize); for (j = 0; j < nDim; j++) { loop[j].inc = 0; loop[j].srcPtr = srcData; loop[j].dstPtr = dstData; } /* * Now loop */ for ( ;; ) { memcpy(loop[0].dstPtr, loop[0].srcPtr, loop[0].length*itemSize); for (j = 1; j < nDim; j++) { loop[j].srcPtr += (loop[j].srcSkip * itemSize); loop[j].dstPtr += (loop[j].dstSkip * itemSize); loop[j].inc ++; if (loop[j].inc < loop[j].length) break; } if (j == nDim) break; for (k = (j - 1); k >= 0; k--) { loop[k].inc = 0; loop[k].srcPtr = loop[j].srcPtr; loop[k].dstPtr = loop[j].dstPtr; } } return dArray; } /* * Note that it is critical to call AddOverlapData with the central * partition first so that the new positions and connections components * get set up correctly. */ static Error AddOverlapData(Field dstField, int *dstCounts, int *meshOffsets, Field srcField, int *indices, LoHi growth, LoHi overlap, LoHi extra, int nDim, char **components) { int i, j, k, itemSize; struct loop pdep_loop[3]; struct loop cdep_loop[3]; struct loop *loop; Array cArray, srcArray, dstArray, origArray = NULL; int srcCounts[3]; char *name; int doit; char **cmp; Class class; Object attr; int spmin[3], spmax[3]; int scmin[3], scmax[3]; int dpmin[3], dpmax[3]; int dcmin[3], dcmax[3]; int *dmin, *dmax, *smin, *smax; int scstride[3], spstride[3], *sstride; int dcstride[3], dpstride[3], *dstride; int offset[3]; int ddep, dref, nRefs, doPositions; int sdep, sref; char origName[64]; if (DXEmptyField(srcField)) return OK; /* * First deal with connections. Get the original and stick it into * "original connections". Then create the new one, install the * approproate attributes and stick it in "connections" */ srcArray = (Array)DXGetComponentValue(srcField, "connections"); if (! srcArray) return OK; if (! DXQueryGridConnections(srcArray, NULL, srcCounts)) { DXSetError(ERROR_MISSING_DATA, "irregular connections"); return ERROR; } if (! DXGetComponentValue(dstField, "original connections")) { if (! DXSetComponentValue(dstField, "original connections", (Object)srcArray)) goto error; dstArray = DXMakeGridConnectionsV(nDim, dstCounts); if (! dstArray) return ERROR; if (! DXSetMeshOffsets((MeshArray)dstArray, meshOffsets)) { DXDelete((Object)dstArray); goto error; } if (! DXCopyAttributes((Object)dstArray, (Object)srcArray)) { DXDelete((Object)dstArray); goto error; } if (! DXSetComponentValue(dstField, "connections", (Object)dstArray)) { DXDelete((Object)dstArray); goto error; } } /* * Given a reference in the source component, we need to * * a) determine whether it lies in the overlap region, and * b) what the index in the grown field will be. * * In the following loop we determine the index ranges * of the included portions of the source field, with * spmin[i] and spmax[i] representing the min and max * position indices along the i-th dimension and, similarly, * scmin[i] and scmax[i] representing the min and max * connection indices along the i-th dimension. Later * this info will enable us to determine whether an element * referred to in the source will be copied into the destination. * * We also get the offset of the source grid into the destination * grid. Given this we can translate a source reference's * indices into a destination's indices, from where we can get * the destinations reference. */ for (i = 0; i < nDim; i++) { switch(indices[i]) { case 0: spmin[i] = (srcCounts[i] - 1) - overlap[i].lo; spmax[i] = (srcCounts[i] - 1); offset[i] = -((srcCounts[i] - 1) - overlap[i].lo); break; case 1: spmin[i] = 0; spmax[i] = (srcCounts[i] - 1); offset[i] = overlap[i].lo; break; case 2: spmin[i] = 0; spmax[i] = overlap[i].hi; offset[i] = (dstCounts[i] - 1) - overlap[i].hi; break; } scmin[i] = spmin[i]; scmax[i] = spmax[i] - 1; dpmin[i] = 0; dpmax[i] = dstCounts[i] - 1; dcmin[i] = 0; dcmax[i] = dpmax[i] - 1; } /* * Now determine the positions and connections strides in both the * source and destination */ scstride[nDim-1] = 1; spstride[nDim-1] = 1; dcstride[nDim-1] = 1; dpstride[nDim-1] = 1; for (i = nDim-2; i >= 0; i--) { scstride[i] = scstride[i+1]*(srcCounts[i+1]-1); spstride[i] = spstride[i+1]*srcCounts[i+1]; dcstride[i] = dcstride[i+1]*(dstCounts[i+1]-1); dpstride[i] = dpstride[i+1]*dstCounts[i+1]; } /* * Set up generic looping. Note that we invert the looping order * here. */ for (i = 0; i < nDim; i++) { switch(indices[(nDim-1)-i]) { case 0: pdep_loop[i].srcStart = (srcCounts[(nDim-1)-i] - overlap[(nDim-1)-i].lo) - 1; pdep_loop[i].dstStart = extra[(nDim-1)-i].lo; pdep_loop[i].length = overlap[(nDim-1)-i].lo; cdep_loop[i].srcStart = ((srcCounts[(nDim-1)-i]-1) - overlap[(nDim-1)-i].lo); cdep_loop[i].dstStart = extra[(nDim-1)-i].lo; cdep_loop[i].length = overlap[(nDim-1)-i].lo; break; case 1: pdep_loop[i].srcStart = 0; pdep_loop[i].dstStart = (extra[(nDim-1)-i].lo + overlap[(nDim-1)-i].lo); pdep_loop[i].length = dstCounts[(nDim-1)-i] - ((overlap[(nDim-1)-i].lo + overlap[(nDim-1)-i].hi) + (extra[(nDim-1)-i].lo + extra[(nDim-1)-i].hi)); cdep_loop[i].srcStart = 0; cdep_loop[i].dstStart = (extra[(nDim-1)-i].lo + overlap[(nDim-1)-i].lo); cdep_loop[i].length = (dstCounts[(nDim-1)-i]-1) - ((overlap[(nDim-1)-i].lo + overlap[(nDim-1)-i].hi) + (extra[(nDim-1)-i].lo + extra[(nDim-1)-i].hi)); break; case 2: pdep_loop[i].srcStart = 1; pdep_loop[i].dstStart = dstCounts[(nDim-1)-i] - (overlap[(nDim-1)-i].hi + extra[(nDim-1)-i].hi); pdep_loop[i].length = overlap[(nDim-1)-i].hi; cdep_loop[i].srcStart = 0; cdep_loop[i].dstStart = (dstCounts[(nDim-1)-i]-1) - (overlap[(nDim-1)-i].hi + extra[(nDim-1)-i].hi); cdep_loop[i].length = overlap[(nDim-1)-i].hi; break; } if (i == 0) { pdep_loop[i].srcSkip = 1; pdep_loop[i].dstSkip = 1; cdep_loop[i].srcSkip = 1; cdep_loop[i].dstSkip = 1; } else { pdep_loop[i].srcSkip = srcCounts[nDim-i] * pdep_loop[i-1].srcSkip; pdep_loop[i].dstSkip = dstCounts[nDim-i] * pdep_loop[i-1].dstSkip; cdep_loop[i].srcSkip = (srcCounts[nDim-i]-1) * cdep_loop[i-1].srcSkip; cdep_loop[i].dstSkip = (dstCounts[nDim-i]-1) * cdep_loop[i-1].dstSkip; } } /* * Now deal with positions */ srcArray = (Array)DXGetComponentValue(srcField, "positions"); dstArray = (Array)DXGetComponentValue(dstField, "positions"); if (dstArray) { /* * If its regular, then we handle it separately. Otherwise, * handle it just like any other irregular component. */ if (DXQueryGridPositions(dstArray, NULL, NULL, NULL, NULL)) { doPositions = 0; /* * If there isn't an "original" version yet, this is the * first time through and we go ahead and grow it. */ if (! DXGetComponentValue(dstField, "original positions")) { if (! DXSetComponentValue(dstField, "original positions", (Object)dstArray)) goto error; dstArray = GrowRegularGrid(srcArray, growth); if (! dstArray) goto error; if (! DXCopyAttributes((Object)dstArray, (Object)srcArray)) { DXDelete((Object)dstArray); goto error; } if (! DXSetComponentValue(dstField, "positions", (Object)dstArray)) { DXDelete((Object)dstArray); goto error; } } } else doPositions = 1; } i = 0; while (NULL != (srcArray = GetEArray(srcField, i++, &name))) { if (! DoIt(name, components)) continue; if (! strcmp(name, "connections")) continue; if (! strcmp(name, "positions") && ! doPositions) continue; if (! GetDepRef(srcArray, &sdep, &sref)) goto error; if (sdep == NOT_DEP && sref == NOT_REFS) continue; /* * Get the destination array */ dstArray = (Array)DXGetComponentValue(dstField, name); if (! dstArray) { /* * OK. There's no equivalent in the destination field to * a component in the source field. Assume that a referential * component will suffice. Better be an invalid component. */ if (strncmp(name, "invalid", 7)) { DXSetError(ERROR_INVALID_DATA, "mismatching components: %s", name); goto error; } dstArray = DXNewArray(TYPE_INT, CATEGORY_REAL, 0); if (sdep == DEP_ON_POSITIONS || sref == REFS_POSITIONS) { if (! DXSetAttribute((Object)dstArray, "ref", (Object)DXNewString("positions"))) { DXDelete((Object)dstArray); goto error; } dref = REFS_POSITIONS; } else { if (! DXSetAttribute((Object)dstArray, "ref", (Object)DXNewString("connections"))) { DXDelete((Object)dstArray); goto error; } dref = REFS_CONNECTIONS; } ddep = NOT_DEP; if (! DXSetComponentValue(dstField, name, (Object)dstArray)) { DXDelete((Object)dstArray); goto error; } } else { if (! GetDepRef(dstArray, &ddep, &dref)) goto error; } if (ddep == DEP_ON_POSITIONS || sdep == DEP_ON_POSITIONS) loop = pdep_loop; else if (ddep == DEP_ON_CONNECTIONS || sdep == DEP_ON_CONNECTIONS) loop = cdep_loop; else loop = NULL; /* * Make sure there's already a saved copy */ SetOrigName(origName, name); if (! DXGetComponentValue(dstField, origName)) { Type type; Category cat; int knt, rank, shape[32]; /* * If not, this is the first time we've encountered * this component. Place dstArray in "original..." and * create a copy in "..." */ origArray = dstArray; DXGetArrayInfo(origArray, &knt, &type, &cat, &rank, shape); dstArray = InitGrowIrregArray(origArray, dstCounts, nDim, sdep); if (! dstArray) goto error; if (! DXCopyAttributes((Object)dstArray, (Object)origArray)) { DXDelete((Object)dstArray); dstArray = NULL; goto error; } if (! DXSetComponentValue(dstField, origName, (Object)origArray)) goto error; if (! DXSetComponentValue(dstField, name, (Object)dstArray)) goto error; origArray = NULL; } if (sref == REFS_POSITIONS || dref == REFS_POSITIONS) { sstride = spstride; dstride = dpstride; smin = spmin; smax = spmax; dmin = dpmin; dmax = dpmax; } else if (sref == REFS_CONNECTIONS || dref == REFS_CONNECTIONS) { sstride = scstride; dstride = dcstride; smin = scmin; smax = scmax; dmin = dcmin; dmax = dcmax; } itemSize = DXGetItemSize(srcArray); nRefs = itemSize / sizeof(int); /* just in case */ if (sdep != NOT_DEP && ddep != NOT_DEP) /* case 0 */ { byte *srcData = (byte *)DXGetArrayData(srcArray); byte *dstData = (byte *)DXGetArrayData(dstArray); /* * Both are dependent. Destination array must already exist. * * Initialize looping */ for (j = 0; j < nDim; j++) { srcData += (loop[j].srcStart * loop[j].srcSkip) * itemSize; dstData += (loop[j].dstStart * loop[j].dstSkip) * itemSize; } for (j = 0; j < nDim; j++) { loop[j].inc = 0; loop[j].srcPtr = srcData; loop[j].dstPtr = dstData; } /* * Now loop */ for ( ;; ) { if (sref == NOT_REFS) { memcpy(loop[0].dstPtr, loop[0].srcPtr, loop[0].length*itemSize); } else { int i, in; int *sPtr = (int *)loop[0].srcPtr; int *dPtr = (int *)loop[0].dstPtr; for (i = 0; i < loop[0].length*nRefs; i++, sPtr++, dPtr++) { int indices[3], dRef, in; if (*sPtr == -1) *dPtr = -1; else { INDICES(*sPtr, indices, sstride, nDim); IS_IN(in, indices, smin, smax, nDim); if (in) { OFFSET(indices, offset, nDim); REFERENCE(indices, *dPtr, dstride, nDim); } else *dPtr = -1; } } } for (j = 1; j < nDim; j++) { loop[j].srcPtr += (loop[j].srcSkip * itemSize); loop[j].dstPtr += (loop[j].dstSkip * itemSize); loop[j].inc ++; if (loop[j].inc < loop[j].length) break; } if (j == nDim) break; for (k = (j - 1); k >= 0; k--) { loop[k].inc = 0; loop[k].srcPtr = loop[j].srcPtr; loop[k].dstPtr = loop[j].dstPtr; } } } else if (sdep == NOT_DEP && ddep != NOT_DEP) /* case 1 */ { /* * source component contains references while the destination * component is dependent. Must be an invalid component. * Dependent destination component must already exist. */ int *srcData = (int *)DXGetArrayData(srcArray); byte *dstData = (byte *)DXGetArrayData(dstArray); int nRefs, dIndex; if ((ddep == DEP_ON_CONNECTIONS && sref != REFS_CONNECTIONS) || (ddep == DEP_ON_POSITIONS && sref != REFS_POSITIONS )) { DXSetError(ERROR_INVALID_DATA, "depref mismatch"); goto error; } DXGetArrayInfo(srcArray, &nRefs, NULL, NULL, NULL, NULL); dstData = (byte *)DXGetArrayData(dstArray); for (j = 0; j < nRefs; j++) { int indices[10]; int dIndex, in; INDICES(srcData[j], indices, sstride, nDim); IS_IN(in, indices, smin, smax, nDim); if (in) { OFFSET(indices, offset, nDim); REFERENCE(indices, dIndex, dstride, nDim); dstData[dIndex] = 1; } } } else if (sdep != NOT_DEP && ddep == NOT_DEP) /* case 2 */ { /* * destination component contains references while the source * component is dependent. Must be an invalid component. */ byte *srcData; int *dstData; int nOldRefs, nNewRefs, sIndex; if ((sdep == DEP_ON_CONNECTIONS && dref != REFS_CONNECTIONS) || (sdep == DEP_ON_POSITIONS && dref != REFS_POSITIONS )) { DXSetError(ERROR_INVALID_DATA, "depref mismatch"); goto error; } /* * Initialize looping to count references */ srcData = (byte *)DXGetArrayData(srcArray); for (j = 0; j < nDim; j++) srcData += loop[j].srcStart * loop[j].srcSkip; for (j = 0; j < nDim; j++) { loop[j].inc = 0; loop[j].srcPtr = srcData; } /* * Now loop, counting references */ nNewRefs = 0; for ( ;; ) { byte *ptr = loop[0].srcPtr; for (j = 0; j < loop[0].length; j++) if (*ptr++) nNewRefs++; for (j = 1; j < nDim; j++) { loop[j].srcPtr += loop[j].srcSkip; loop[j].inc ++; if (loop[j].inc < loop[j].length) break; } if (j == nDim) break; for (k = (j - 1); k >= 0; k--) { loop[k].inc = 0; loop[k].srcPtr = loop[j].srcPtr; } } DXGetArrayInfo(dstArray, &nOldRefs, NULL, NULL, NULL, NULL); if (! DXAddArrayData(dstArray, nOldRefs, nNewRefs, NULL)) goto error; dstData = ((int *)DXGetArrayData(dstArray)) + nOldRefs; for (j = 0; j < nDim; j++) { loop[j].inc = 0; loop[j].srcPtr = srcData; } for ( ;; ) { byte *ptr = loop[0].srcPtr; for (j = 0; j < loop[0].length; j++) if (*ptr++) { int indices[32], k; for (k = 0; k < nDim-1; k++) indices[k] = loop[(nDim-1)-k].dstStart + loop[(nDim-1)-k].inc; indices[nDim-1] = loop[0].dstStart + j; REFERENCE(indices, k, dstride, nDim); *dstData++ = k; } for (j = 1; j < nDim; j++) { loop[j].srcPtr += loop[j].srcSkip; loop[j].inc ++; if (loop[j].inc < loop[j].length) break; } if (j == nDim) break; for (k = (j - 1); k >= 0; k--) { loop[k].inc = 0; loop[k].srcPtr = loop[j].srcPtr; } } } else if (sdep == NOT_DEP && ddep == NOT_DEP) { /* * Not a dependent array, but contains refs */ Type type; Category cat; int rank, shape[32]; int *dPtr, *sPtr = (int *)DXGetArrayData(srcArray); int nDstItems = 0, newDstItems = 0, nSrcItems, in; DXGetArrayInfo(srcArray, &nSrcItems, &type, &cat, &rank, shape); /* * Count refs that will be copied into the overlap region */ for (j = 0, newDstItems = 0; j < nSrcItems; j++, sPtr += nRefs) { for (k = 0, in = 0; !in && k < nRefs; k++) { int indices[3]; if (sPtr[k] == -1) continue; INDICES(sPtr[k], indices, sstride, nDim); IS_IN(in, indices, smin, smax, nDim); } if (in) newDstItems ++; } if (newDstItems) { if (! dstArray) { dstArray = DXNewArrayV(type, cat, rank, shape); if (! dstArray) goto error; if (dref == REFS_POSITIONS) { if (! DXSetAttribute((Object)dstArray, "ref", (Object)DXNewString("positions"))) goto error; } else { if (! DXSetAttribute((Object)dstArray, "ref", (Object)DXNewString("connections"))) goto error; } if (! DXSetComponentValue(dstField, name, (Object)dstArray)) { DXDelete((Object)dstArray); goto error; } nDstItems = 0; } else DXGetArrayInfo(dstArray, &nDstItems, NULL,NULL,NULL,NULL); if (! DXAddArrayData(dstArray, nDstItems, newDstItems, NULL)) return NULL; sPtr = (int *)DXGetArrayData(srcArray); dPtr = ((int *)DXGetArrayData(dstArray))+nDstItems; /* * Count refs that will be copied into the overlap region */ for (j = 0, newDstItems = 0; j < nSrcItems; j++, sPtr += nRefs) { for (k = 0, in = 0; !in && k < nRefs; k++) { int indices[3]; if (sPtr[k] == -1) continue; INDICES(sPtr[k], indices, sstride, nDim); IS_IN(in, indices, smin, smax, nDim); } if (in) { for (k = 0, in = 0; !in && k < nRefs; k++) { int indices[3]; if (sPtr[k] == -1) dPtr[k] = -1; INDICES(sPtr[k], indices, sstride, nDim); IS_IN(in, indices, smin, smax, nDim); if (in) { OFFSET(indices, offset, nDim); REFERENCE(indices, dPtr[k], dstride, nDim); } } dPtr += nRefs; } } } } } return OK; error: DXDelete((Object)origArray); return ERROR; } static Array InitGrowIrregArray(Array oArray, int *nCounts, int nDim, int dep) { Type type; Category cat; int rank, shape[32]; int i, nItems; Array nArray; if (! DXGetArrayInfo(oArray, NULL, &type, &cat, &rank, shape)) return NULL; nArray = DXNewArrayV(type, cat, rank, shape); if (! nArray) return NULL; if (dep != NOT_DEP) { if (dep == DEP_ON_POSITIONS) { nItems = 1; for (i = 0; i < nDim; i++) nItems *= nCounts[i]; } else { nItems = 1; for (i = 0; i < nDim; i++) nItems *= (nCounts[i]-1); } if (! DXAddArrayData(nArray, 0, nItems, NULL)) { DXDelete((Object)nArray); return NULL; } memset(DXGetArrayData(nArray), 0, nItems); } return nArray; } static Array GrowRegularGrid(Array oArray, LoHi growth) { int oCounts[3], nCounts[3]; float oOrigin[3], nOrigin[3]; float deltas[9], *d; Array nArray; int nDim, oDim, i, j; DXQueryGridPositions(oArray, &nDim, oCounts, oOrigin, deltas); /* * Origin is original, minus growth[i].lo*deltas; */ for (i = 0; i < nDim; i++) nOrigin[i] = oOrigin[i]; d = deltas; oDim = 0; for (i = 0; i < nDim; i++) { if (oCounts[i] > 1) { /* * DXAdd in counts for growth at either end of this axis, as * required by the growth array */ nCounts[i] = oCounts[i] + (growth[oDim].lo+growth[oDim].hi); /* * If we grew along the lower end of this axis, we need to * back the origin off appropriately. */ if (growth[oDim].lo) for (j = 0; j < nDim; j++) nOrigin[j] -= growth[oDim].lo * *d++; else d += nDim; oDim ++; } else { nCounts[i] = 1; d += nDim; } } nArray = DXMakeGridPositionsV(nDim, nCounts, nOrigin, deltas); return nArray; } static Array GrowRegularArray(Array oArray, int *counts, int nDim, int dep) { int len, i, n; int count; float o[32]; float d[32]; Type type; if (! DXGetRegularArrayInfo((RegularArray)oArray, NULL, (Pointer)o, (Pointer)d)) return NULL; DXGetArrayInfo(oArray, NULL, &type, NULL, NULL, &len); if (type != TYPE_FLOAT) return NULL; for (i = 0; i < len; i++) if (d[i] != 0.0) return NULL; if (dep == DEP_ON_POSITIONS) { n = 1; for (i = 0; i < nDim; i++) n *= counts[i]; } else { n = 1; for (i = 0; i < nDim; i++) n *= (counts[i]-1); } return (Array)DXNewRegularArray(type, len, n, (Pointer)o, (Pointer)d); } static Error AddPartition(HashTable ht, Field partition) { int i, nDim, inc; int box[8][3]; inc = 0; if (DXEmptyField(partition)) return OK; if (!GetIJKBox(partition, (int *)box, &nDim)) return ERROR; if (nDim == 1) inc = 4; else if (nDim == 2) inc = 2; else if (nDim == 3) inc = 1; for (i = 0; i < 8; i += inc) { if (! AddPartitionEntry(ht, (int *)(box + i), i, partition)) return ERROR; } return OK; } static Error GetIJKBox(Field partition, int *box, int *nDim) { int i, j, n; int counts[3], offsets[3]; Array array; array = (Array)DXGetComponentValue(partition, "connections"); if (! array) { DXSetError(ERROR_MISSING_DATA, "no connections component"); return ERROR; } if (!DXQueryGridConnections(array, &n, counts) || n > 3 || !DXGetMeshOffsets((MeshArray)array, offsets)) { DXSetError(ERROR_INVALID_DATA,"regular connections of dim <= 3 expected"); return ERROR; } for (i = n; i < 3; i++) counts[i] = offsets[i] = 0; /* * Generate (i,j,k) coordinates of corners of partition based on * the origin of the original partitioned field. Use the three bits * of the corner index i to select which counts we add in. Theres a * reversal happening here: if i == 0x1, we add in the Z count while * if i == 0x4, we dd in X count */ for (i = 0; i < 8; i++, box += 3) { box[0] = offsets[0]; box[1] = offsets[1]; box[2] = offsets[2]; for (j = 0; j < 3; j++) if (i & (1<partition; } static PseudoKey PartHash(Key key) { struct partHash *p; register long l; p = (struct partHash *)key; l = p->base[0]*17 + p->base[1]*107 + p->base[2]*557 + p->tag*1039; return (PseudoKey) l; } static int PartCmp(Key k0, Key k1) { struct partHash *p0, *p1; p0 = (struct partHash *)k0; p1 = (struct partHash *)k1; if (p0->base[0] != p1->base[0]) return 1; if (p0->base[1] != p1->base[1]) return 1; if (p0->base[2] != p1->base[2]) return 1; if (p0->tag != p1->tag) return 1; return 0; } Field DXQueryOriginalMeshExtents(Field f, int *offsets, int *sizes) { int oOffsets[32], grownOffsets[32]; Array orig, grown; int nDim, i; orig = (Array)DXGetComponentValue(f, "original connections"); if (! orig) return NULL; if (sizes || offsets) { if (! DXQueryGridConnections(orig, &nDim, sizes)) { DXSetError(ERROR_INVALID_DATA, "regular connections required"); return NULL; } } if (offsets) { grown = (Array)DXGetComponentValue(f, "connections"); if (! grown) { DXSetError(ERROR_INVALID_DATA, "grown connections component not found"); return NULL; } if (!DXGetMeshOffsets((MeshArray)orig, oOffsets)) { DXSetError(ERROR_INVALID_DATA, "error retrieving original mesh offset"); return NULL; } if (!DXGetMeshOffsets((MeshArray)grown, grownOffsets)) { DXSetError(ERROR_INVALID_DATA, "error retrieving grown mesh offset"); return NULL; } for (i = 0; i < nDim; i++) offsets[i] = oOffsets[i] - grownOffsets[i]; } return f; } static Error FillEmptyOverlap(Field field, int *indices, LoHi growth, int nDim, char **comp, Pointer fill) { int i, j, k, itemSize; struct loop cdep_loop[3]; struct loop pdep_loop[3]; struct loop *loop; Array cArray, array; int counts[3], cdep_sizes[3], pdep_sizes[3]; byte *srcData, *dstData; char *name; int doit; char **cmp; int cstFill; Object attr; char *dep; cstFill = (fill > (Pointer)0x3); cArray = (Array)DXGetComponentValue(field, "connections"); if (! cArray) return NULL; if (! DXQueryGridConnections(cArray, NULL, counts)) { DXSetError(ERROR_MISSING_DATA, "FillEmptyOverlap: irregular connections"); return ERROR; } cdep_sizes[nDim-1] = pdep_sizes[nDim-1] = 1; for (i = nDim-2; i >= 0; i--) { cdep_sizes[i] = cdep_sizes[i+1] * (counts[i+1] - 1); pdep_sizes[i] = pdep_sizes[i+1] * counts[i+1]; } /* * Set up generic looping. */ for (i = 0, j = 0; i < nDim; i++) { switch(indices[i]) { case 0: pdep_loop[nDim-1].srcStart = growth[i].lo * pdep_sizes[i]; pdep_loop[nDim-1].dstStart = 0; pdep_loop[nDim-1].length = growth[i].lo; pdep_loop[nDim-1].srcSkip = pdep_sizes[i]; cdep_loop[nDim-1].srcStart = growth[i].lo * cdep_sizes[i]; cdep_loop[nDim-1].dstStart = 0; cdep_loop[nDim-1].length = growth[i].lo; cdep_loop[nDim-1].srcSkip = cdep_sizes[i]; break; case 1: pdep_loop[j].srcStart = 0; pdep_loop[j].dstStart = 0; pdep_loop[j].length = counts[i]; pdep_loop[j].srcSkip = pdep_sizes[i]; cdep_loop[j].srcStart = 0; cdep_loop[j].dstStart = 0; cdep_loop[j].length = (counts[i] - 1); cdep_loop[j].srcSkip = cdep_sizes[i]; j++; break; case 2: pdep_loop[nDim-1].srcStart = ((counts[i] - growth[i].hi) - 1) * pdep_sizes[i]; pdep_loop[nDim-1].dstStart = (counts[i] - growth[i].hi) * pdep_sizes[i]; pdep_loop[nDim-1].length = growth[i].hi; pdep_loop[nDim-1].srcSkip = pdep_sizes[i]; cdep_loop[nDim-1].srcStart = (((counts[i] - 1) - growth[i].hi) - 1) * cdep_sizes[i]; cdep_loop[nDim-1].dstStart = ((counts[i] - 1) - growth[i].hi) * cdep_sizes[i]; cdep_loop[nDim-1].length = growth[i].hi; cdep_loop[nDim-1].srcSkip = cdep_sizes[i]; break; } } i = 0; while (NULL != (array = GetEArray(field, i++, &name))) { doit = 0; if (! strcmp(name, "positions")) doit = 1; if (comp) for (cmp = comp; !doit && *cmp; cmp++) if (! strcmp(*cmp, name)) { doit = 1; break; } if (! doit) continue; attr = DXGetComponentAttribute(field, name, "dep"); if (! attr) continue; if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_MISSING_DATA, "invalid dependency attribute"); return ERROR; } dep = DXGetString((String)attr); if (! strcmp(dep, "connections")) loop = cdep_loop; else if (! strcmp(dep, "positions")) loop = pdep_loop; else { DXSetError(ERROR_MISSING_DATA, "invalid dependency attribute"); return ERROR; } if (DXQueryGridPositions(array, NULL, NULL, NULL, NULL)) continue; /* * If its a regular array and we are replicating the edge value, * then we are already done. If its regular, then if we are inserting * a fill value, we may need to irregularize the array. */ if (DXGetArrayClass(array) == CLASS_REGULARARRAY) { Type type; Category cat; int rank, shape[32]; int nItems; int size, ok; if (! cstFill) continue; size = DXGetItemSize(array); DXGetArrayInfo(array, &nItems, &type, &cat, &rank, shape); #define CHECK(type, ok) \ { \ type *o = NULL, *d = NULL; \ int i; \ \ o = (type *)DXAllocate(size); \ d = (type *)DXAllocate(size); \ if (! o|| ! d) \ { \ DXFree((Pointer)o); \ DXFree((Pointer)d); \ goto error; \ } \ \ DXGetRegularArrayInfo((RegularArray)array, NULL, (Pointer)o, (Pointer)d); \ \ for (i = 0, ok = 1; i < shape[0] && ok; i++) \ { \ if (d[i] != (type)0) \ ok = 0; \ if (o[i] != ((type *)fill)[i]) \ ok = 0; \ } \ \ DXFree((Pointer)o); \ DXFree((Pointer)d); \ } switch(type) { case TYPE_DOUBLE: CHECK(double, ok); break; case TYPE_FLOAT: CHECK(float, ok); break; case TYPE_INT: CHECK(int, ok); break; case TYPE_SHORT: CHECK(short, ok); break; case TYPE_UBYTE: CHECK(ubyte, ok); break; } if (ok) continue; #define IRREGULARIZE(dx_type, type) \ { \ type *o = NULL, *d = NULL, *ptr; \ Array nArray = NULL; \ int i, j; \ \ nArray = DXNewArrayV(dx_type, cat, rank, shape); \ if (! nArray) \ goto error; \ \ if (! DXAddArrayData(nArray, 0, nItems, NULL)) \ { \ DXDelete((Object)nArray); \ goto error; \ } \ \ ptr = (type *)DXGetArrayData(nArray) + loop[nDim-1].srcStart*itemSize; \ if (! ptr) \ { \ DXDelete((Object)nArray); \ goto error; \ } \ \ o = (type *)DXAllocate(size); \ d = (type *)DXAllocate(size); \ if (! o || ! d) \ { \ DXFree((Pointer)o); \ DXFree((Pointer)d); \ goto error; \ } \ \ DXGetRegularArrayInfo((RegularArray)array, NULL, (Pointer)o, (Pointer)d); \ \ for (i = 0; i < nItems - loop[nDim-1].srcStart; i++) \ for (j = 0; j < shape[0]; j++) \ *ptr++ = o[j] + i*d[j]; \ \ \ DXSetComponentValue(field, name, (Object)nArray); \ \ DXFree((Pointer)o); \ DXFree((Pointer)d); \ } switch(type) { case TYPE_DOUBLE: IRREGULARIZE(TYPE_DOUBLE, double); break; case TYPE_FLOAT: IRREGULARIZE(TYPE_FLOAT, float); break; case TYPE_INT: IRREGULARIZE(TYPE_INT, int); break; case TYPE_SHORT: IRREGULARIZE(TYPE_SHORT, short); break; case TYPE_UBYTE: IRREGULARIZE(TYPE_UBYTE, ubyte); break; } } itemSize = DXGetItemSize(array); dstData = srcData = (byte *)DXGetArrayData(array); dstData += loop[nDim-1].dstStart * itemSize; if (cstFill) srcData = (byte *)fill; else srcData += loop[nDim-1].srcStart * itemSize; for (j = 0; j < nDim; j++) { loop[j].inc = 0; loop[j].srcPtr = srcData; loop[j].dstPtr = dstData; } /* * Now loop */ for ( ;; ) { if (! cstFill) srcData = loop[0].srcPtr; dstData = loop[0].dstPtr; for (j = 0; j < loop[0].length; j++) { memcpy(dstData, srcData, itemSize); if (! cstFill) srcData += loop[0].srcSkip * itemSize; dstData += loop[0].srcSkip * itemSize; } for (j = 1; j < nDim; j++) { if (j < (nDim-1)) loop[j].srcPtr += (loop[j].srcSkip * itemSize); loop[j].dstPtr += (loop[j].srcSkip * itemSize); loop[j].inc ++; if (loop[j].inc < loop[j].length) break; } if (j == nDim) break; for (k = (j - 1); k >= 0; k--) { loop[k].inc = 0; loop[k].srcPtr = loop[j].srcPtr; loop[k].dstPtr = loop[j].dstPtr; } } } return OK; error: return ERROR; } static int DoIt(char *name, char **components) { if (! strcmp(name, "positions") || ! strcmp(name, "invalid positions") || ! strcmp(name, "connections") || ! strcmp(name, "invalid connections")) return 1; if (components) while (*components) { if (! strcmp(*components, name)) return 1; components ++; } return 0; } static Error GetDepRef(Array array, int *dep, int *ref) { Object attr; attr = DXGetAttribute((Object)array, "ref"); if (attr) { if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10200", "ref attribute"); return ERROR; } if (!strcmp(DXGetString((String)attr), "positions")) *ref = REFS_POSITIONS; else if (! strcmp(DXGetString((String)attr), "connections")) *ref = REFS_CONNECTIONS; else *ref = NOT_REFS; } else *ref = NOT_REFS; attr = DXGetAttribute((Object)array, "dep"); if (attr) { if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10200", "dep attribute"); return ERROR; } if (!strcmp(DXGetString((String)attr), "positions")) *dep = DEP_ON_POSITIONS; else if (! strcmp(DXGetString((String)attr), "connections")) *dep = DEP_ON_CONNECTIONS; else *dep = NOT_DEP; } else *dep = NOT_DEP; return OK; } typedef struct { int index; int limit; int base; int stride; } Loop; static Error DupTask(Pointer p) { Field field = (Field)p; int mo[32], mc[32], stride[32], i, j, k, d, l, nd, done; InvalidComponentHandle ich = NULL; Loop loop[32]; Array a; a = (Array)DXGetComponentValue(field, "invalid positions"); if (a) { if (! DXSetComponentValue(field, "saved invalid positions", (Object)a)) goto error; } a = (Array)DXGetComponentValue(field, "connections"); if (! a) return OK; if (! DXGetMeshOffsets((MeshArray)a, mo)) goto error; if (! DXQueryGridConnections(a, &nd, mc)) goto error; stride[nd-1] = 1; for (i = nd-2; i >= 0; i--) stride[i] = stride[i+1]*mc[i+1]; ich = DXCreateInvalidComponentHandle((Object)field, NULL, "positions"); if (! ich) goto error; if (nd == 1) { if (mc[0] != 0) if (! DXSetElementInvalid(ich, 0)) goto error; } else { for (i = 0; i < nd; i++) { if (mo[i] != 0) { for (j = k = 0; j < nd; j++) if (i != j) { loop[k].index = 0; loop[k].limit = mc[j]; loop[k].base = 0; loop[k].stride = stride[j]; k++; } done = 0; while (! done) { int b = loop[0].base; for (j = 0; j < loop[0].limit; j++, b += loop[0].stride) if (! DXSetElementInvalid(ich, b)) goto error; for (j = 1; j < (nd-1); j++) if (++loop[j].index < loop[j].limit) { loop[j].base += loop[j].stride; for (k = j-1; k >= 0; k--) { loop[k].base = loop[j].base; loop[k].index = 0; } break; } done = (j == (nd-1)); } } } } if (! DXSaveInvalidComponent(field, ich)) goto error; DXFreeInvalidComponentHandle(ich); return OK; error: if (ich) DXFreeInvalidComponentHandle(ich); return ERROR; } Error _dxfRegInvalidateDupBoundary(CompositeField cf) { Field f; int i; if (! DXCreateTaskGroup()) goto error; i = 0; while (NULL != (f = (Field)DXGetEnumeratedMember((Group)cf, i++, NULL))) if (! DXAddTask(DupTask, f, 0, 1.0)) goto error; if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; return OK; error: return ERROR; }