/***********************************************************************/ /* 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: /home/gda/dxcvs/dx/src/exec/dxmods/_refineirr.c,v 1.3 1999/05/10 15:45:20 gda Exp $: */ #include #include #include "math.h" #include /*************************** * cubes tables ***************************/ static int cube_edgeTable[] = { 0, 1, 1, 5, 5, 4, 4, 0, 2, 3, 3, 7, 7, 6, 6, 2, 0, 2, 1, 3, 5, 7, 4, 6 }; static int cube_faceTable[] = { 0, 1, 2, 3, 1, 5, 7, 3, 5, 4, 6, 7, 0, 4, 2, 6, 1, 5, 4, 0, 2, 6, 7, 3 }; static int subcubes[] = { 0, 8, 16, 20, 11, 24, 23, 26, 8, 1, 20, 17, 24, 9, 26, 21, 16, 20, 2, 12, 23, 26, 15, 25, 20, 17, 12, 3, 26, 21, 25, 13, 11, 24, 23, 26, 4, 10, 19, 22, 24, 9, 26, 21, 10, 5, 22, 18, 23, 26, 15, 25, 19, 22, 6, 14, 26, 21, 25, 13, 22, 18, 14, 7 }; /*************************** * quads tables ***************************/ static int quad_edgeTable[] = { 0, 1, 1, 3, 3, 2, 2, 0 }; static int subquads[] = { 0, 4, 7, 8, 4, 1, 8, 5, 7, 8, 2, 6, 8, 5, 6, 3 }; /*************************** * lines tables ***************************/ static int sublines[] = { 0, 2, 2, 1 }; /*************************** * tetras tables ***************************/ static int tetra_edgeTable[] = { 0, 1, 1, 3, 1, 2, 2, 3, 0, 2, 0, 3 }; static int subtetras[] = { 0, 4, 8, 9, 4, 1, 6, 5, 7, 6, 8, 2, 9, 5, 7, 3, 5, 4, 7, 6, 4, 8, 7, 6, 4, 7, 8, 9, 4, 5, 7, 9 }; /*************************** * triangles tables ***************************/ static int tri_edgeTable[] = { 0, 1, 1, 2, 2, 0 }; static int subtriangles[] = { 0, 3, 5, 3, 1, 4, 3, 4, 5, 5, 4, 2 }; typedef struct hash { int index; int v[4]; } *Hash; /* * One element structure for all hash elements, but * each will be sized appropriately. */ #define EDGESIZE 3*sizeof(int) #define QUADSIZE 5*sizeof(int) /* * Tables containing the following information are used * to identify the interpolations necessary to create the * additional positional information. The table is a list * of nV * knt vertex indices. For each of knt results, * the data associated with nV vertices are averaged. */ typedef struct { int *table; /* list of vertex indices to be interpolated */ int nV; /* number of vertices to interpolate for each output */ int knt; /* number of interpolations to perform */ } InterpTable; Field _dxfRefineReg(Field, int); static Field RefineICubes(Field, int); static Field RefineIQuads(Field, int); static Field RefineILines(Field, int); static Field RefineITetras(Field, int); static Field RefineITriangles(Field, int); static Array RefineCArray(Array, int); static Array RefineCRefArray(Array, int); static Array RefinePArray(Array, InterpTable *, int); static Array CopyArray(Array); static Error FieldSetup(Field); static int *MakeTable(HashTable, int, int); /* * The following is for keeping track of position derivations so * that invalid positions references are propagated. */ typedef struct listElement { struct listElement *next; int index; } *ListElement; #define VALID(table, index) ((((int **)table)[index]) == NULL) static Error AddHashReferences(HashTable, ListElement *, SegList *, int, int); static Error AddReferences(int *, int, ListElement *, SegList *, int, int); static Error SetupIPTables(Field, int, ListElement **, SegList **); static Error DumpIPTables(Field, int, ListElement **, SegList **); static Error AddQuadToHash(HashTable, int, int, int, int, int *); static int QueryQuad(HashTable, int, int, int, int); static void SetupQuad(int, int, int, int, Hash); static Error AddEdgeToHash(HashTable, int, int, int *); static int QueryEdge(HashTable, int, int); static void SetupEdge(int, int, Hash); static Hash _QueryHash(HashTable, Hash); static int Cmp4(Key, Key); static PseudoKey Hash4(Key); static int Cmp2(Key, Key); static PseudoKey Hash2(Key); Field _dxfRefineIrreg(Field f, int levels) { Object attr; char *str; if (! DXGetComponentValue(f, "connections")) { DXSetError(ERROR_MISSING_DATA, "#10250", "input", "connections"); return NULL; } attr = DXGetComponentAttribute(f, "connections", "element type"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "#10255", "connections", "element type"); return NULL; } if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "#10051", "element type attribute"); return NULL; } str = DXGetString((String)attr); if (! strcmp(str, "cubes")) return RefineICubes(f, levels); else if (! strcmp(str, "quads")) return RefineIQuads(f, levels); else if (! strcmp(str, "lines")) return RefineILines(f, levels); else if (! strcmp(str, "tetrahedra")) return RefineITetras(f, levels); else if (! strcmp(str, "triangles")) return RefineITriangles(f, levels); else { DXSetError(ERROR_INVALID_DATA, "#11380", str); return NULL; } } static Field RefineICubes(Field f, int levels) { int i, j, k, l; HashTable fHash = NULL; HashTable eHash = NULL; InterpTable iTable[4]; Array inP, inC, outC = NULL; int *c, *inElts, *outElts, nElts, nPoints; int eOffset, fOffset, cOffset; int nEdges, nFaces; Object old; char *name; ListElement *invTable = NULL; SegList *listElements = NULL; iTable[0].table = iTable[1].table = NULL; if (! FieldSetup(f)) goto error; for (l = 0; l < levels; l++) { inP = (Array)DXGetComponentValue(f, "positions"); DXGetArrayInfo(inP, &nPoints, NULL, NULL, NULL, NULL); inC = (Array)DXGetComponentValue(f, "connections"); if (! inC) { DXSetError(ERROR_MISSING_DATA, "10250", "input", "connections"); goto error; } inElts = (int *)DXGetArrayData(inC); if (! inElts) goto error; DXGetArrayInfo(inC, &nElts, NULL, NULL, NULL, NULL); outC = (Array)DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 8); if (! outC) goto error; if (! DXAddArrayData(outC, 0, 8*nElts, NULL)) goto error; outElts = (int *)DXGetArrayData(outC); if (! outElts) goto error; /* * Create hash tables */ fHash = DXCreateHash(QUADSIZE, Hash4, Cmp2); if (! fHash) goto error; eHash = DXCreateHash(EDGESIZE, Hash2, Cmp2); if (! eHash) goto error; /* * Put edges and faces into to the hash tables */ c = inElts; nEdges = 0; nFaces = 0; for (i = 0; i < nElts; i++, c += 8) { int *t; t = cube_edgeTable; for (j = 0; j < 12; j++, t += 2) if (! AddEdgeToHash(eHash, c[t[0]], c[t[1]], &nEdges)) goto error; t = cube_faceTable; for (j = 0; j < 6; j++, t += 4) if (! AddQuadToHash(fHash, c[t[0]], c[t[1]], c[t[2]], c[t[3]], &nFaces)) goto error; } /* * for each input cube, create grid containing indices in output * positions component and resulting output cubes. We will be * producing the input vertices first, the edge vertices next, the * face vertices next and the center vertices last. */ eOffset = nPoints; fOffset = eOffset + nEdges; cOffset = fOffset + nFaces; c = inElts; for (i = 0; i < nElts; i++, c += 8) { int grid[27]; int v, *t; int *b = subcubes; v = 0; for (j = 0; j < 8; j++) grid[v++] = c[j]; t = cube_edgeTable; for (j = 0; j < 12; j++, t += 2) { int indx = QueryEdge(eHash, c[t[0]], c[t[1]]); grid[v++] = indx + eOffset; } t = cube_faceTable; for (j = 0; j < 6; j++, t += 4) { int indx = QueryQuad(fHash, c[t[0]], c[t[1]], c[t[2]], c[t[3]]); grid[v++] = indx + fOffset; } grid[v] = cOffset + i; /* * Now the grid is stuffed with the output vertex IDs. * Generate the output cubes. */ for (j = 0; j < 8; j++) for (k = 0; k < 8; k++) *outElts++ = grid[*b++]; } if (! SetupIPTables(f, nPoints, &invTable, &listElements)) goto error; iTable[0].nV = 2; iTable[0].knt = nEdges; iTable[0].table = MakeTable(eHash, nEdges, 2); if (! iTable[0].table) goto error; if (invTable) if (! AddHashReferences(eHash, invTable, listElements, nPoints, 2)) goto error; iTable[1].nV = 4; iTable[1].knt = nFaces; iTable[1].table = MakeTable(fHash, nFaces, 4); if (! iTable[1].table) goto error; if (invTable) if (! AddHashReferences(fHash, invTable, listElements, nPoints+nEdges, 4)) goto error; iTable[2].nV = 8; iTable[2].knt = nElts; iTable[2].table = inElts; if (! iTable[2].table) goto error; if (invTable) if (! AddReferences(inElts, nElts, invTable, listElements, nPoints+nEdges+nFaces, 8)) goto error; iTable[3].table = NULL; DXDestroyHash(eHash); DXDestroyHash(fHash); eHash = fHash = NULL; for (i=0; NULL != (old = DXGetEnumeratedComponentValue(f, i, &name)); i++) { Object attr; Array new = (Array)old; char *str; if (! strcmp(name, "connections")) continue; attr = DXGetComponentAttribute(f, name, "dep"); if (attr) { int boolean; boolean = ! strncmp(name, "invalid", 7); if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "#10051", "dependency attribute"); goto error; } str = DXGetString((String)attr); if (! strcmp(str, "positions")) { new = RefinePArray((Array)old, iTable, boolean); if (! new) goto error; } else if (! strcmp(str, "connections")) { new = RefineCArray((Array)old, 8); if (! new) goto error; } } else { if (! strcmp(name, "invalid connections")) { /* * Must be referential */ new = RefineCRefArray((Array)old, 8); if (! new) goto error; } } /* * It is possible that the refinement process created new * arrays. This is because you cannot extend regular arrays * and they are used for constants. */ if ((Object)old != (Object)new) if (! DXSetComponentValue(f, name, (Object)new)) goto error; } /* * We wait until now to reset the connections component since it * causes deletion of the old one, whose data is used in iTable[2] * to designate the vertices to be interpolated for each cube-center * point. */ if (! DXSetComponentValue(f, "connections", (Object)outC)) goto error; outC = NULL; if (! DumpIPTables(f, nPoints, &invTable, &listElements)) goto error; DXFree((Pointer)iTable[0].table); DXFree((Pointer)iTable[1].table); iTable[0].table = iTable[1].table = NULL; } return DXEndField(f); error: if (listElements) DXDeleteSegList(listElements); if (invTable) DXFree((Pointer)invTable); if (eHash) DXDestroyHash(eHash); if (fHash) DXDestroyHash(fHash); DXFree((Pointer)iTable[0].table); DXFree((Pointer)iTable[1].table); return NULL; } static Field RefineIQuads(Field f, int levels) { int i, j, k, l; HashTable eHash = NULL; InterpTable iTable[3]; Array inP, inC, outC = NULL; int *c, *inElts, *outElts, nElts, nPoints; int eOffset, cOffset; int nEdges; Object old; char *name; ListElement *invTable = NULL; SegList *listElements = NULL; InvalidComponentHandle ich = NULL; iTable[0].table = iTable[1].table = NULL; if (! FieldSetup(f)) goto error; for (l = 0; l < levels; l++) { inP = (Array)DXGetComponentValue(f, "positions"); DXGetArrayInfo(inP, &nPoints, NULL, NULL, NULL, NULL); inC = (Array)DXGetComponentValue(f, "connections"); if (! inC) { DXSetError(ERROR_MISSING_DATA, "10250", "input", "connections"); goto error; } if (DXGetComponentValue(f, "invalid connections")) { ich = DXCreateInvalidComponentHandle((Object)f, NULL, "connections"); if (! ich) goto error; } inElts = (int *)DXGetArrayData(inC); if (! inElts) goto error; DXGetArrayInfo(inC, &nElts, NULL, NULL, NULL, NULL); outC = (Array)DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 4); if (! outC) goto error; if (! DXAddArrayData(outC, 0, 4*nElts, NULL)) goto error; outElts = (int *)DXGetArrayData(outC); if (! outElts) goto error; eHash = DXCreateHash(EDGESIZE, Hash2, Cmp2); if (! eHash) goto error; /* * Put edges into to the hash tables */ c = inElts; nEdges = 0; for (i = 0; i < nElts; i++, c += 4) { int *t; t = quad_edgeTable; for (j = 0; j < 4; j++, t += 2) if (! AddEdgeToHash(eHash, c[t[0]], c[t[1]], &nEdges)) goto error; } /* * for each input quad, create grid containing indices in output * positions component and resulting output quad. We will be * producing the input vertices first, the edge vertices next, the * center vertices last. */ eOffset = nPoints; cOffset = eOffset + nEdges; c = inElts; for (i = 0; i < nElts; i++, c += 4) { int grid[9]; int *t, v; int *b = subquads; v = 0; for (j = 0; j < 4; j++) grid[v++] = c[j]; t = quad_edgeTable; for (j = 0; j < 4; j++, t += 2) { int indx = QueryEdge(eHash, c[t[0]], c[t[1]]); grid[v++] = indx + eOffset; } grid[v] = cOffset + i; /* * Now the grid is stuffed with the output vertex IDs. * Generate the output quads. */ for (j = 0; j < 4; j++) for (k = 0; k < 4; k++) *outElts++ = grid[*b++]; } if (! SetupIPTables(f, nPoints, &invTable, &listElements)) goto error; iTable[0].nV = 2; iTable[0].knt = nEdges; iTable[0].table = MakeTable(eHash, nEdges, 2); if (! iTable[0].table) goto error; if (invTable) if (! AddHashReferences(eHash, invTable, listElements, nPoints, 2)) goto error; iTable[1].nV = 4; iTable[1].knt = nElts; iTable[1].table = inElts; if (! iTable[1].table) goto error; if (invTable) if (! AddReferences(inElts, nElts, invTable, listElements, nPoints+nEdges, 4)) goto error; iTable[2].table = NULL; DXDestroyHash(eHash); eHash = NULL; for (i=0; NULL != (old = DXGetEnumeratedComponentValue(f, i, &name)); i++) { Object attr; Array new = (Array)old; char *str; if (! strcmp(name, "connections")) continue; attr = DXGetComponentAttribute(f, name, "dep"); if (attr) { int boolean; boolean = ! strncmp(name, "invalid", 7); if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10051", "dependency attribute"); goto error; } str = DXGetString((String)attr); if (! strcmp(str, "positions")) { new = RefinePArray((Array)old, iTable, boolean); if (! new) goto error; } else if (! strcmp(str, "connections")) { new = RefineCArray((Array)old, 4); if (! new) goto error; } } else { if (! strcmp(name, "invalid connections")) { /* * Must be referential */ new = RefineCRefArray((Array)old, 4); if (! new) goto error; } } /* * It is possible that the refinement process created new * arrays. This is because you cannot extend regular arrays * and they are used for constants. */ if ((Object)old != (Object)new) if (! DXSetComponentValue(f, name, (Object)new)) goto error; } /* * We wait until now to reset the connections component since it * causes deletion of the old one, whose data is used in iTable[2] * to designate the vertices to be interpolated for each cube-center * point. */ if (! DXSetComponentValue(f, "connections", (Object)outC)) goto error; outC = NULL; if (! DumpIPTables(f, nPoints, &invTable, &listElements)) goto error; DXFree((Pointer)iTable[0].table); iTable[0].table = NULL; } return DXEndField(f); error: if (listElements) DXDeleteSegList(listElements); if (invTable) DXFree((Pointer)invTable); if (eHash) DXDestroyHash(eHash); DXFree((Pointer)iTable[0].table); return NULL; } static Field RefineILines(Field f, int levels) { int i, j, k, l; InterpTable iTable[2]; Array inP, inC, outC = NULL; int *c, *inElts, *outElts, nElts, nPoints; int cOffset; int nEdges; Object old; char *name; ListElement *invTable = NULL; SegList *listElements = NULL; iTable[0].table = iTable[1].table = NULL; if (DXEmptyField(f)) return DXEndField(DXNewField()); if (! FieldSetup(f)) goto error; for (l = 0; l < levels; l++) { inP = (Array)DXGetComponentValue(f, "positions"); DXGetArrayInfo(inP, &nPoints, NULL, NULL, NULL, NULL); inC = (Array)DXGetComponentValue(f, "connections"); if (! inC) { DXSetError(ERROR_MISSING_DATA, "10250", "input", "connections"); goto error; } inElts = (int *)DXGetArrayData(inC); if (! inElts) goto error; DXGetArrayInfo(inC, &nElts, NULL, NULL, NULL, NULL); outC = (Array)DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 2); if (! outC) goto error; if (! DXAddArrayData(outC, 0, 2*nElts, NULL)) goto error; outElts = (int *)DXGetArrayData(outC); if (! outElts) goto error; /* * for each input line, create grid containing indices in output * positions component and resulting output line. We will be * producing the input vertices first and the center vertices last. */ cOffset = nPoints; c = inElts; for (i = 0; i < nElts; i++, c += 2) { int grid[3]; int *t, v; int *b = sublines; v = 0; for (j = 0; j < 2; j++) grid[v++] = c[j]; grid[v] = cOffset + i; /* * Now the grid is stuffed with the output vertex IDs. * Generate the output cubes. */ for (j = 0; j < 2; j++) for (k = 0; k < 2; k++) *outElts++ = grid[*b++]; } if (! SetupIPTables(f, nPoints, &invTable, &listElements)) goto error; iTable[0].nV = 2; iTable[0].knt = nElts; iTable[0].table = inElts; if (! iTable[0].table) goto error; iTable[1].table = NULL; for (i=0; NULL != (old = DXGetEnumeratedComponentValue(f, i, &name)); i++) { Object attr; Array new = (Array)old; char *str; if (! strcmp(name, "connections")) continue; attr = DXGetComponentAttribute(f, name, "dep"); if (attr) { int boolean; boolean = ! strncmp(name, "invalid", 7); if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10051", "dependency attribute"); goto error; } str = DXGetString((String)attr); if (! strcmp(str, "positions")) { new = RefinePArray((Array)old, iTable, boolean); if (! new) goto error; } else if (! strcmp(str, "connections")) { new = RefineCArray((Array)old, 2); if (! new) goto error; } } else { if (! strcmp(name, "invalid connections")) { /* * Must be referential */ new = RefineCRefArray((Array)old, 2); if (! new) goto error; } } /* * It is possible that the refinement process created new * arrays. This is because you cannot extend regular arrays * and they are used for constants. */ if ((Object)old != (Object)new) if (! DXSetComponentValue(f, name, (Object)new)) goto error; } /* * We wait until now to reset the connections component since it * causes deletion of the old one, whose data is used in iTable[2] * to designate the vertices to be interpolated for each cube-center * point. */ if (! DXSetComponentValue(f, "connections", (Object)outC)) goto error; outC = NULL; } return DXEndField(f); error: if (listElements) if (listElements) DXDeleteSegList(listElements); if (invTable) DXFree((Pointer)invTable); DXDeleteSegList(listElements); if (invTable) DXFree((Pointer)invTable); return NULL; } static Field RefineITetras(Field f, int levels) { int i, j, k, l; HashTable eHash = NULL; InterpTable iTable[3]; Array inP, inC, outC = NULL; int *c, *inElts, *outElts, nElts, nPoints; int eOffset, cOffset; int nEdges; Object old; char *name; ListElement *invTable = NULL; SegList *listElements = NULL; iTable[0].table = iTable[1].table = NULL; if (DXEmptyField(f)) return DXEndField(DXNewField()); if (! FieldSetup(f)) goto error; for (l = 0; l < levels; l++) { inP = (Array)DXGetComponentValue(f, "positions"); DXGetArrayInfo(inP, &nPoints, NULL, NULL, NULL, NULL); inC = (Array)DXGetComponentValue(f, "connections"); if (! inC) { DXSetError(ERROR_MISSING_DATA, "10250", "input", "connections"); goto error; } inElts = (int *)DXGetArrayData(inC); if (! inElts) goto error; DXGetArrayInfo(inC, &nElts, NULL, NULL, NULL, NULL); outC = (Array)DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 4); if (! outC) goto error; if (! DXAddArrayData(outC, 0, 8*nElts, NULL)) goto error; outElts = (int *)DXGetArrayData(outC); if (! outElts) goto error; eHash = DXCreateHash(EDGESIZE, Hash2, Cmp2); if (! eHash) goto error; /* * Put edges into to the hash tables */ c = inElts; nEdges = 0; for (i = 0; i < nElts; i++, c += 4) { int *t; t = tetra_edgeTable; for (j = 0; j < 6; j++, t += 2) if (! AddEdgeToHash(eHash, c[t[0]], c[t[1]], &nEdges)) goto error; } /* * for each input element, create grid containing indices in output * positions component and resulting output elements. We will be * producing the input vertices first, the edge vertices next. */ eOffset = nPoints; cOffset = eOffset + nEdges; c = inElts; for (i = 0; i < nElts; i++, c += 4) { int grid[10]; int *t, v; int *b = subtetras; v = 0; for (j = 0; j < 4; j++) grid[v++] = c[j]; t = tetra_edgeTable; for (j = 0; j < 6; j++, t += 2) { int indx = QueryEdge(eHash, c[t[0]], c[t[1]]); grid[v++] = indx + eOffset; } /* * Now the grid is stuffed with the output vertex IDs. * Generate the output tetras. */ for (j = 0; j < 8; j++) for (k = 0; k < 4; k++) *outElts++ = grid[*b++]; } if (! SetupIPTables(f, nPoints, &invTable, &listElements)) goto error; iTable[0].nV = 2; iTable[0].knt = nEdges; iTable[0].table = MakeTable(eHash, nEdges, 2); if (! iTable[0].table) goto error; if (invTable) if (! AddHashReferences(eHash, invTable, listElements, nPoints, 2)) goto error; iTable[1].table = NULL; DXDestroyHash(eHash); eHash = NULL; for (i=0; NULL != (old = DXGetEnumeratedComponentValue(f, i, &name)); i++) { Object attr; Array new = (Array)old; char *str; if (! strcmp(name, "connections")) continue; attr = DXGetComponentAttribute(f, name, "dep"); if (attr) { int boolean; boolean = ! strncmp(name, "invalid", 7); if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10051", "dependency attribute"); goto error; } str = DXGetString((String)attr); if (! strcmp(str, "positions")) { new = RefinePArray((Array)old, iTable, boolean); if (! new) goto error; } else if (! strcmp(str, "connections")) { new = RefineCArray((Array)old, 8); if (! new) goto error; } } else { if (! strcmp(name, "invalid connections")) { /* * Must be referential */ new = RefineCRefArray((Array)old, 8); if (! new) goto error; } } /* * It is possible that the refinement process created new * arrays. This is because you cannot extend regular arrays * and they are used for constants. */ if ((Object)old != (Object)new) if (! DXSetComponentValue(f, name, (Object)new)) goto error; } /* * We wait until now to reset the connections component since it * causes deletion of the old one, whose data is used in iTable[2] * to designate the vertices to be interpolated for each cube-center * point. */ if (! DXSetComponentValue(f, "connections", (Object)outC)) goto error; if (! DumpIPTables(f, nPoints, &invTable, &listElements)) goto error; outC = NULL; DXFree((Pointer)iTable[0].table); iTable[0].table = NULL; } return DXEndField(f); error: if (listElements) DXDeleteSegList(listElements); if (invTable) DXFree((Pointer)invTable); if (eHash) DXDestroyHash(eHash); DXFree((Pointer)iTable[0].table); return NULL; } static Field RefineITriangles(Field f, int levels) { int i, j, k, l; HashTable eHash = NULL; InterpTable iTable[3]; Array inP, inC, outC = NULL; int *c, *inElts, *outElts, nElts, nPoints; int eOffset, cOffset; int nEdges; Object old; char *name; ListElement *invTable = NULL; SegList *listElements = NULL; iTable[0].table = iTable[1].table = NULL; if (DXEmptyField(f)) return DXEndField(DXNewField()); if (! FieldSetup(f)) goto error; for (l = 0; l < levels; l++) { inP = (Array)DXGetComponentValue(f, "positions"); DXGetArrayInfo(inP, &nPoints, NULL, NULL, NULL, NULL); inC = (Array)DXGetComponentValue(f, "connections"); if (! inC) { DXSetError(ERROR_MISSING_DATA, "10250", "input", "connections"); goto error; } inElts = (int *)DXGetArrayData(inC); if (! inElts) goto error; DXGetArrayInfo(inC, &nElts, NULL, NULL, NULL, NULL); outC = (Array)DXNewArray(TYPE_INT, CATEGORY_REAL, 1, 3); if (! outC) goto error; if (! DXAddArrayData(outC, 0, 4*nElts, NULL)) goto error; outElts = (int *)DXGetArrayData(outC); if (! outElts) goto error; eHash = DXCreateHash(EDGESIZE, Hash2, Cmp2); if (! eHash) goto error; /* * Put edges into to the hash tables */ c = inElts; nEdges = 0; for (i = 0; i < nElts; i++, c += 3) { int *t; t = tri_edgeTable; for (j = 0; j < 3; j++, t += 2) if (! AddEdgeToHash(eHash, c[t[0]], c[t[1]], &nEdges)) goto error; } /* * for each input element, create grid containing indices in output * positions component and resulting output elements. We will be * producing the input vertices first, the edge vertices next. */ eOffset = nPoints; cOffset = eOffset + nEdges; c = inElts; for (i = 0; i < nElts; i++, c += 3) { int grid[6]; int *t, v; int *b = subtriangles; v = 0; for (j = 0; j < 3; j++) grid[v++] = c[j]; t = tri_edgeTable; for (j = 0; j < 3; j++, t += 2) { int indx = QueryEdge(eHash, c[t[0]], c[t[1]]); grid[v++] = indx + eOffset; } /* * Now the grid is stuffed with the output vertex IDs. * Generate the output tetras. */ for (j = 0; j < 4; j++) for (k = 0; k < 3; k++) *outElts++ = grid[*b++]; } if (! SetupIPTables(f, nPoints, &invTable, &listElements)) goto error; iTable[0].nV = 2; iTable[0].knt = nEdges; iTable[0].table = MakeTable(eHash, nEdges, 2); if (! iTable[0].table) goto error; if (invTable) if (! AddHashReferences(eHash, invTable, listElements, nPoints, 2)) goto error; iTable[1].table = NULL; DXDestroyHash(eHash); eHash = NULL; for (i=0; NULL != (old = DXGetEnumeratedComponentValue(f, i, &name)); i++) { Object attr; Array new = (Array)old; char *str; if (! strcmp(name, "connections")) continue; attr = DXGetComponentAttribute(f, name, "dep"); if (attr) { int boolean; boolean = ! strncmp(name, "invalid", 7); if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10051", "dependency attribute"); goto error; } str = DXGetString((String)attr); if (! strcmp(str, "positions")) { new = RefinePArray((Array)old, iTable, boolean); if (! new) goto error; } else if (! strcmp(str, "connections")) { new = RefineCArray((Array)old, 4); if (! new) goto error; } } else { if (! strcmp(name, "invalid connections")) { /* * Must be referential */ new = RefineCRefArray((Array)old, 4); if (! new) goto error; } } /* * It is possible that the refinement process created new * arrays. This is because you cannot extend regular arrays * and they are used for constants. */ if ((Object)old != (Object)new) if (! DXSetComponentValue(f, name, (Object)new)) goto error; } /* * We wait until now to reset the connections component since it * causes deletion of the old one, whose data is used in iTable[2] * to designate the vertices to be interpolated for each cube-center * point. */ if (! DXSetComponentValue(f, "connections", (Object)outC)) goto error; if (! DumpIPTables(f, nPoints, &invTable, &listElements)) goto error; outC = NULL; DXFree((Pointer)iTable[0].table); iTable[0].table = NULL; } DXFree((Pointer)iTable[0].table); return DXEndField(f); error: if (listElements) DXDeleteSegList(listElements); if (invTable) DXFree((Pointer)invTable); if (eHash) DXDestroyHash(eHash); DXFree((Pointer)iTable[0].table); return NULL; } static int * MakeTable(HashTable h, int knt, int nVerts) { int *table; Hash hash; table = (int *)DXAllocate(knt*nVerts*sizeof(int)); if (! table) goto error; if (! DXInitGetNextHashElement(h)) goto error; while (NULL != (hash = (Hash)DXGetNextHashElement(h))) { int i, *tPtr = table + hash->index*nVerts; for (i = 0; i < nVerts; i++) *tPtr++ = hash->v[i]; } return table; error: if (table) DXFree((Pointer)table); return NULL; } static Error FieldSetup(Field f) { char *rmv[64]; int nrmv = 0, nkeep = 0, i; Object o; Array a; char *name; for (i = 0; NULL != (o = DXGetEnumeratedComponentValue(f, i, &name)); i++) { /* * if its connections, keep it. */ if (! strcmp(name, "connections")) continue; /* * If it is an invalid array, keep it */ if (! strncmp(name, "invalid", 7)) goto keep; /* * Otherwise, if it refs, get rid of it. */ if (DXGetComponentAttribute(f, name, "ref")) { rmv[nrmv++] = name; continue; } /* * If its derived, get rid of it. */ if (DXGetComponentAttribute(f, name, "der")) { rmv[nrmv++] = name; continue; } keep: /* * Otherwise, replace it with a copy. */ a = CopyArray((Array)o); if (! a) goto error; if (! DXSetComponentValue(f, name, (Object)a)) goto error; } /* * Now delete the components we don't want */ for (i = 0; i < nrmv; i++) DXDeleteComponent(f, rmv[i]); return OK; error: return ERROR; } static Array CopyArray(Array in) { int nD, nIn; Array out; Type type; Category cat; int rank, shape[32], i; if (DXGetObjectClass((Object)in) != CLASS_ARRAY) goto error; DXGetArrayInfo(in, &nIn, &type, &cat, &rank, shape); /* * If the array is constant, it will have to be copied * each time through the loop, anyway, so its not necessary to copy * it here. */ if (DXQueryConstantArray(in, NULL, NULL)) return in; /* * Neither grid positions or constant, go create * an irregular output. */ out = DXNewArrayV(type, cat, rank, shape); if (! out) goto error; if (! DXAddArrayData(out, 0, nIn, DXGetArrayData(in))) goto error; done: return out; error: if (out) DXDelete((Object)out); return NULL; } #define AVERAGE(type, itable) \ { \ type *inPtr = (type *)DXGetArrayData(out); \ type *outPtr; \ int i, j, k, t, l, m; \ \ if (! inPtr) \ goto error; \ \ outPtr = inPtr + nIn*nVars; \ for (t = 0; itable[t].table != NULL; t++) \ { \ int nV = itable[t].nV; \ int knt = itable[t].knt; \ int *tab = itable[t].table; \ float d = 1.0 / itable[t].nV; \ \ for (i = 0; i < knt; i++, tab += nV) \ { \ type *xtab[8]; \ \ for (j = 0; j < nV; j++) \ xtab[j] = inPtr + tab[j]*nVars; \ \ for (k = 0; k < nVars; k++) \ { \ *outPtr = 0.0; \ for (l = 0; l < nV; l++) \ { \ type max = *xtab[l]; \ int maxi = -1; \ for (m = l+1; m < nV; m++) \ if (max < *xtab[m]) \ { \ maxi = m; \ max = *xtab[m]; \ } \ if (maxi != -1) \ { \ type *t; \ t = xtab[maxi]; \ xtab[maxi] = xtab[l]; \ xtab[l] = t; \ } \ *outPtr += max; \ xtab[l] += 1; \ } \ *outPtr++ *= d; \ } \ } \ } \ } #define BOOLEAN(itable) \ { \ ubyte *inPtr = (ubyte *)DXGetArrayData(out); \ ubyte *outPtr; \ int i, j, k, t, l, m; \ \ if (! inPtr) \ goto error; \ \ outPtr = inPtr + nIn; \ for (t = 0; itable[t].table != NULL; t++) \ { \ int nV = itable[t].nV; \ int knt = itable[t].knt; \ int *tab = itable[t].table; \ \ for (i = 0; i < knt; i++, tab += nV, outPtr++) \ { \ *outPtr = 0; \ for (j = 0; j < nV; j++) \ if (inPtr[tab[j]] == 1) { \ *outPtr = 1; \ break; \ } \ } \ } \ } static Array RefinePArray(Array in, InterpTable *table, int boolean) { int nD, nIn; Array out; Type type; Category cat; int rank, shape[32], i; int nNew = 0, nVars, itemSize; if (DXGetObjectClass((Object)in) != CLASS_ARRAY) goto error; DXGetArrayInfo(in, &nIn, &type, &cat, &rank, shape); itemSize = DXGetItemSize(in); for (i = 0; table[i].table != NULL; i ++) nNew += table[i].knt; if (DXQueryConstantArray(in, NULL, NULL)) { out = (Array)DXNewConstantArrayV(nIn+nNew, DXGetConstantArrayData(in), type, cat, rank, shape); } else { out = in; if (! DXAddArrayData(out, nIn, nNew, NULL)) goto error; nVars = itemSize / DXTypeSize(type); if (boolean) { if (type != TYPE_BYTE && type != TYPE_UBYTE) { DXSetError(ERROR_INVALID_DATA, "invalid data flag must be byte or ubyte"); goto error; } if (nVars != 1) { DXSetError(ERROR_INVALID_DATA, "invalid data flag must be scalar or 1-vector"); goto error; } BOOLEAN(table); } else { switch(type) { case TYPE_UBYTE: AVERAGE(ubyte, table); break; case TYPE_BYTE: AVERAGE(byte, table); break; case TYPE_USHORT: AVERAGE(ushort, table); break; case TYPE_SHORT: AVERAGE(short, table); break; case TYPE_UINT: AVERAGE(uint, table); break; case TYPE_INT: AVERAGE(int, table); break; case TYPE_FLOAT: AVERAGE(float, table); break; case TYPE_DOUBLE: AVERAGE(double, table); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "input"); goto error; } } } return out; error: if (out) DXDelete((Object)out); return NULL; } static Array RefineCArray(Array in, int nCopies) { int nD, nIn; Array out; Type type; Category cat; int rank, shape[32], i, j; int itemSize; byte *inPtr, *outPtr; if (DXGetObjectClass((Object)in) != CLASS_ARRAY) goto error; DXGetArrayInfo(in, &nIn, &type, &cat, &rank, shape); if (DXQueryConstantArray(in, NULL, NULL)) { out = (Array)DXNewConstantArrayV(nIn*nCopies, DXGetConstantArrayData(in), type, cat, rank, shape); } else { itemSize = DXGetItemSize(in); out = DXNewArrayV(type, cat, rank, shape); if (! out) goto error; if (! DXAddArrayData(out, 0, nIn*nCopies, NULL)) goto error; inPtr = (byte *)DXGetArrayData(in); outPtr = (byte *)DXGetArrayData(out); if (! inPtr || ! outPtr) goto error; for (i = 0; i < nIn; i++, inPtr += itemSize) for (j = 0; j < nCopies; j++, outPtr += itemSize) memcpy(outPtr, inPtr, itemSize); } return out; error: if (out) DXDelete((Object)out); return NULL; } static Hash _QueryHash(HashTable h, Hash hash) { return (Hash)DXQueryHashElement(h, (Key)hash); } static void SetupEdge(int p, int q, Hash e) { if (p > q) { e->v[0] = p; e->v[1] = q; } else { e->v[0] = q; e->v[1] = p; } } static int QueryEdge(HashTable h, int p, int q) { struct hash hash; Hash hptr; SetupEdge(p, q, &hash); hptr = _QueryHash(h, &hash); if (hptr) return hptr->index; else return -1; } static Error AddEdgeToHash(HashTable h, int p, int q, int *nEdges) { struct hash hash; SetupEdge(p, q, &hash); if (_QueryHash(h, &hash)) return OK; hash.index = *nEdges; *nEdges = *nEdges + 1; return DXInsertHashElement(h, (Element)&hash); } static void SetupQuad(int p, int q, int r, int s, Hash h) { int i; h->v[0] = p; h->v[1] = q; h->v[2] = r; h->v[3] = s; for (i = 0; i < 3; i++) { int max = h->v[i]; int j, k = -1; for (j = i+1; j < 4; j++) if (h->v[j] > max) { k = j; max = h->v[j]; } if (k != -1) { int t = h->v[i]; h->v[i] = h->v[k]; h->v[k] = t; } } } static int QueryQuad(HashTable h, int p, int q, int r, int s) { struct hash hash; Hash hptr; SetupQuad(p, q, r, s, &hash); hptr = _QueryHash(h, &hash); if (hptr) return hptr->index; else return -1; } static Error AddQuadToHash(HashTable h, int p, int q, int r, int s, int *nQuad) { struct hash hash; SetupQuad(p, q, r, s, &hash); if (_QueryHash(h, &hash)) return OK; hash.index = *nQuad; *nQuad = *nQuad + 1; return DXInsertHashElement(h, (Element)&hash); } static PseudoKey Hash4(Key k) { Hash h = (Hash)k; return h->v[0]*17 + h->v[1]*53 + h->v[2]*1907 + h->v[3]*2801; } static int Cmp4(Key k0, Key k1) { Hash h0 = (Hash)k0; Hash h1 = (Hash)k1; return ! (h0->v[0]==h1->v[0] && h0->v[1]==h1->v[1] && h0->v[2]==h1->v[2] && h0->v[3]==h1->v[3]); } static PseudoKey Hash2(Key k) { Hash h = (Hash)k; return h->v[0]*17 + h->v[1]*53; } static int Cmp2(Key k0, Key k1) { Hash h0 = (Hash)k0; Hash h1 = (Hash)k1; return ! (h0->v[0]==h1->v[0] && h0->v[1]==h1->v[1]); } static Error AddHashReferences(HashTable h, ListElement *invTable, SegList *listElements, int base, int nv) { int *table; Hash hash; int n = base; if (! DXInitGetNextHashElement(h)) goto error; while (NULL != (hash = (Hash)DXGetNextHashElement(h))) { int i; for (i = 0; i < nv; i++) { int index = hash->v[i]; if (! VALID(invTable, index)) { ListElement le = DXNewSegListItem(listElements); if (! le) goto error; le->index = n; le->next = invTable[index]; invTable[index] = le; } } n++; } return OK; error: return ERROR; } static Error AddReferences(int *list, int knt, ListElement *invTable, SegList *listElements, int base, int nv) { int i, j; for (i = 0; i < knt; i++) { for (j = 0; j < nv; j++) { int index = *list++; if (! VALID(invTable, index)) { ListElement le = DXNewSegListItem(listElements); if (! le) goto error; le->index = base + i; le->next = invTable[index]; invTable[index] = le; } } } return OK; error: return ERROR; } static Error SetupIPTables(Field f, int nPoints, ListElement **invTable, SegList **listElements) { InvalidComponentHandle iph = NULL; *invTable = NULL; *listElements = NULL; if (DXGetComponentAttribute(f, "invalid positions", "ref")) { int i, index, *iPtr; *invTable = (ListElement *)DXAllocateZero(nPoints * sizeof(Pointer)); if (! *invTable) goto error; iph = DXCreateInvalidComponentHandle((Object)f, NULL, "positions"); if (! iph) goto error; iPtr = (int *)*invTable; DXInitGetNextInvalidElementIndex(iph); while (-1 != (index = DXGetNextInvalidElementIndex(iph))) iPtr[index] = 1; *listElements = DXNewSegList(sizeof(struct listElement), nPoints*2, nPoints); if (! *listElements) goto error; DXFreeInvalidComponentHandle(iph); iph = NULL; DXDeleteComponent(f, "invalid positions"); } return OK; error: if (iph) DXFreeInvalidComponentHandle(iph); if (*invTable) DXFree((Pointer)*invTable); if (*listElements) DXDeleteSegList(*listElements); *invTable = NULL; *listElements = NULL; return ERROR; } static Error DumpIPTables(Field f, int nPoints, ListElement **invTable, SegList **listElements) { InvalidComponentHandle iph = NULL; if (*invTable) { int i; ListElement le; iph = DXCreateInvalidComponentHandle((Object)f, NULL, "positions"); if (! iph) goto error; for (i = 0; i < nPoints; i++) if (NULL != (le = (*invTable)[i])) { if (! DXSetElementInvalid(iph, i)) goto error; do { if (! DXSetElementInvalid(iph, le->index)) goto error; le = le->next; } while (((int)le) != 1); } if (! DXSaveInvalidComponent(f, iph)) goto error; DXFreeInvalidComponentHandle(iph); iph = NULL; DXDeleteSegList(*listElements); *listElements = NULL; DXFree((Pointer)*invTable); *invTable = NULL; } return OK; error: if (iph) DXFreeInvalidComponentHandle(iph); return ERROR; } static Array RefineCRefArray(Array in, int outPerIn) { Type t; Category c; int r, s[32], nIn, i, nRefs, *sPtr, *dPtr; Array out = NULL; if (! DXGetArrayInfo(in, &nIn, &t, &c, &r, s)) goto error; if (t != TYPE_INT && t != TYPE_UINT) { DXSetError(ERROR_INVALID_DATA, "ref components must be int or uint"); goto error; } out = DXNewArrayV(t, c, r, s); if (! out) goto error; if (! DXAddArrayData(out, 0, nIn*outPerIn, NULL)) goto error; nRefs = (DXGetItemSize(in) / sizeof(int)); if (nRefs != 1) { DXSetError(ERROR_INVALID_DATA, "invalid references must be scalar or 1-vector"); goto error; } sPtr = (int *)DXGetArrayData(in); dPtr = (int *)DXGetArrayData(out); for (i = 0; i < nIn; i++, sPtr ++) { int j; for (j = 0; j < outPerIn; j++) *dPtr++ = (*sPtr * outPerIn) + j; } return out; error: DXDelete((Object)out); return NULL; }