/***********************************************************************/ /* 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/_divcurl.c,v 1.3 1999/05/10 15:45:17 gda Exp $ */ #include #include #include #include static Object _dxfDivCurlObject(Object, int); static Error _dxfDivCurlField(Field, int); static Error _dxfDivCurlTask(Pointer ptr); static Error _dxfDivCurlRegular(Field, int *, float *, int); static Error _dxfDivCurlCubesRegular(Field, int *, float *); static Error _dxfDivCurlIrregular(Field); static Error _dxfRecursiveSetGroupTypeV(Object, Type, Category, int, int *); static void _dxfTetrasDivCurl(ArrayHandle, Pointer, Type, float *, float *, byte *, int *); static void _dxfCubesDivCurl(ArrayHandle, Pointer, Type, float *, float *, byte *, int *); typedef struct { Field field; int flags; } Task; #define DO_DIV 0x01 #define DO_CURL 0x02 #define TRUE 1 #define FALSE 0 Object _dxfDivCurl(Object object, Object *div, Object *curl) { int flags = 0; Object copy = NULL; if (div) flags |= DO_DIV; if (curl) flags |= DO_CURL; copy = DXCopy(object, COPY_STRUCTURE); if (! copy) goto error; if (! _dxfDivCurlObject(copy, flags)) goto error; if (div) { *div = DXCopy(copy, COPY_STRUCTURE); if (*div == NULL) goto error; if (! DXRename(*div, "div", "data")) goto error; if (curl) if (! DXRemove(*div, "curl")) goto error; } if (curl) { *curl = DXCopy(copy, COPY_STRUCTURE); if (*curl == NULL) goto error; if (! DXRename(*curl, "curl", "data")) goto error; if (div) if (! DXRemove(*curl, "div")) goto error; } if (copy != object) DXDelete((Object)copy); return object; error: if (div && *div) { DXDelete(*div); *div = NULL; } if (curl && *curl) { DXDelete(*curl); *curl = NULL; } if (copy != object) DXDelete(copy); return NULL; } static Object _dxfDivCurlObject(Object object, int flags) { int i; Field field; Object child; Class class; Type type; Category cat; int rank, shape[32]; Task task; class = DXGetObjectClass(object); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)object); switch (class) { case CLASS_FIELD: if (! _dxfDivCurlField((Field)object, flags)) return NULL; else return object; case CLASS_COMPOSITEFIELD: if (! DXGrow(object, 1, NULL, "data", NULL)) return NULL; DXCreateTaskGroup(); task.flags = flags; i = 0; while(NULL != (field = DXGetPart(object, i++))) { task.field = field; if (! DXAddTask(_dxfDivCurlTask, (Pointer)&task, sizeof(Task), 1.0)) { DXAbortTaskGroup(); return NULL; } } if (!DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) return NULL; if (! DXShrink(object)) return NULL; return object; case CLASS_MULTIGRID: case CLASS_SERIES: case CLASS_GROUP: i = 0; while (NULL != (child=DXGetEnumeratedMember((Group)object,i++,NULL))) if (! _dxfDivCurlObject(child, flags)) return NULL; return object; case CLASS_XFORM: if (! DXGetXformInfo((Xform)object, &child, 0)) return NULL; if (! _dxfDivCurlObject(child, flags)) return NULL; else return object; case CLASS_CLIPPED: if (! DXGetClippedInfo((Clipped)object, &child, 0)) return NULL; if (! _dxfDivCurlObject(child, flags)) return NULL; else return object; default: DXSetError(ERROR_INVALID_DATA, "#11381"); return NULL; } } static Error _dxfDivCurlTask(Pointer ptr) { Task *task = (Task *)ptr; return _dxfDivCurlField(task->field, task->flags); } static Error _dxfDivCurlField(Field field, int flags) { Type type; Category cat; int dim, rank, shape[32]; Array positions; Array connections; Array data; float del[9]; int counts[3]; int grid; int nPoints, pDim; Object attr; if (DXEmptyField(field)) return OK; data = (Array)DXGetComponentValue(field, "data"); if (! data) { DXSetError(ERROR_MISSING_DATA, "#10240", "data"); return ERROR; } DXGetArrayInfo(data, NULL, &type, &cat, &rank, shape); if (cat != CATEGORY_REAL) { DXSetError (ERROR_INVALID_DATA, "#11150", "data"); return ERROR; } if (rank != 1 || shape[0] != 3) { DXSetError (ERROR_INVALID_DATA, "#10230", "data", 3); return ERROR; } positions = (Array)DXGetComponentValue(field, "positions"); connections = (Array)DXGetComponentValue(field, "connections"); DXGetArrayInfo(positions, &nPoints, NULL, NULL, NULL, &pDim); if (rank != 1 && shape[0] == pDim) { DXSetError (ERROR_INVALID_DATA, "#10230", "data", pDim); return ERROR; } attr = DXGetComponentAttribute(field, "data", "dep"); if (! attr) { DXSetError(ERROR_INVALID_DATA, "#10255", "data", "dep"); return ERROR; } else if (strcmp(DXGetString((String)attr), "positions")) { DXSetError(ERROR_INVALID_DATA, "#11251", "data"); return ERROR; } DXReference(attr); /* * Set up div and curl arrays */ if (flags & DO_DIV) { Array array; array = (Array)DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0); if (! array) return ERROR; if (! DXAddArrayData(array, 0, nPoints, NULL)) { DXDelete((Object)array); return ERROR; } if (! DXSetComponentValue(field, "div", (Object)array)) { DXDelete((Object)array); return ERROR; } if (! DXSetComponentAttribute(field, "div", "dep", attr)) return ERROR; } if (flags & DO_CURL) { Array array; array = (Array)DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! array) return ERROR; if (! DXAddArrayData(array, 0, nPoints, NULL)) { DXDelete((Object)array); return ERROR; } if (! DXSetComponentValue(field, "curl", (Object)array)) { DXDelete((Object)array); return ERROR; } if (! DXSetComponentAttribute(field, "curl", "dep", attr)) return ERROR; } DXDelete(attr); /* * We need to ensure that there are no invalid positions referenced by * a valid connection element. */ if (! DXInvalidateUnreferencedPositions((Object)field)) return ERROR; /* * If the grid is regular and axis-aligned, we can use a faster method. * Unless there are any invalid positions. */ if (DXQueryGridConnections(connections, NULL, NULL) && DXQueryGridPositions(positions, &dim, counts, NULL, del) && ! DXGetComponentValue(field, "invalid positions")) { switch(dim) { case 1: grid = 1; break; case 2: grid = del[0] != 0.0 && del[1] == 0.0 && del[2] == 0.0 && del[3] != 0.0; break; case 3: grid = del[0] != 0.0 && del[1] == 0.0 && del[2] == 0.0 && del[3] == 0.0 && del[4] != 0.0 && del[5] == 0.0 && del[6] == 0.0 && del[7] == 0.0 && del[8] != 0.0; break; default: grid = 0; } } else grid = 0; if (grid) { if (! _dxfDivCurlRegular(field, counts, del, dim)) return ERROR; } else { if (! _dxfDivCurlIrregular(field)) return ERROR; } if (DXGetComponentValue(field, "original data")) DXDeleteComponent(field, "original data"); return OK; } static Error _dxfRecursiveSetGroupTypeV (Object parent, Type type, Category cat, int rank, int *shape) { Object child; int i; if (DXGetObjectClass(parent) == CLASS_GROUP) { for (i = 0; NULL != (child = DXGetEnumeratedMember((Group) parent, i, NULL)); i++) { if (! _dxfRecursiveSetGroupTypeV (child, type, cat, rank, shape)) return ERROR; } if (! DXSetGroupTypeV ((Group) parent, type, cat, rank, shape)) return ERROR; } return (OK); } static Error _dxfDivCurlRegular(Field field, int *counts, float *deltas, int nDim) { Object attr; /* * Make sure data is dep on positions */ attr = DXGetComponentAttribute(field, "data", "dep"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "#10255", "data", "dependency"); return ERROR; } if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10200", "dependency attribute"); return ERROR; } if (strcmp(DXGetString((String)attr), "positions")) { DXSetError(ERROR_INVALID_DATA, "#11251", "data"); return ERROR; } switch(nDim) { case 3: return _dxfDivCurlCubesRegular(field, counts, deltas); default: DXSetError(ERROR_INVALID_DATA, "#10230", "positions", 3); return ERROR; } } #define REG_CURL(dx, dy, dz, curl) \ { \ curl[0] = (dy[2] - dz[1]); \ curl[1] = (dz[0] - dx[2]); \ curl[2] = (dx[1] - dy[0]); \ } #define REG_DIV(dx, dy, dz, div) \ { \ *(div) = (dx[0] + dy[1] + dz[2]); \ } static Error _dxfDivCurlCubesRegular(Field field, int *counts, float *deltas) { float dxInv = 1.0 / deltas[0]; float dyInv = 1.0 / deltas[4]; float dzInv = 1.0 / deltas[8]; int xKnt = counts[0]; int yKnt = counts[1]; int zKnt = counts[2]; float *div, *curl; Array inArray, dArray, cArray; Type type; int i, j, k; int lx, ly, lz, nx, ny, nz; float dx, dy, dz; int nItems; inArray = (Array)DXGetComponentValue(field, "data"); DXGetArrayInfo(inArray, &nItems, &type, NULL, NULL, NULL); dArray = (Array)DXGetComponentValue(field, "div"); cArray = (Array)DXGetComponentValue(field, "curl"); if (dArray) { div = (float *)DXGetArrayData(dArray); if (!div || DXGetError()) goto error; } else div = NULL; if (cArray) { curl = (float *)DXGetArrayData(cArray); if (!curl || DXGetError()) goto error; } else curl = NULL; #define REG_CUBES_DIVCURL(type) \ { \ type *data, *last, *next; \ int xSkip, ySkip, zSkip; \ float dx[3], dy[3], dz[3]; \ float xStep, yStep, zStep; \ \ data = (type *)DXGetArrayData(inArray); \ if (! data || DXGetError()) \ goto error; \ \ zSkip = 3; \ ySkip = zSkip * counts[2]; \ xSkip = ySkip * counts[1]; \ \ for (i = 0; i < xKnt; i++) \ { \ xStep = dxInv; \ if (i == 0) \ { \ lx = 0; \ nx = xSkip; \ xStep = dxInv; \ } \ else if (i == (xKnt - 1)) \ { \ lx = xSkip; \ nx = 0; \ xStep = dxInv; \ } \ else \ { \ lx = xSkip; \ nx = xSkip; \ xStep = 0.5 * dxInv; \ } \ \ for (j = 0; j < yKnt; j++) \ { \ yStep = dyInv; \ if (j == 0) \ { \ ly = 0; \ ny = ySkip; \ yStep = dyInv; \ } \ else if (j == (yKnt - 1)) \ { \ ly = ySkip; \ ny = 0; \ yStep = dyInv; \ } \ else \ { \ ly = ySkip; \ ny = ySkip; \ yStep = 0.5 * dyInv; \ } \ \ for (k = 0; k < zKnt; k++) \ { \ zStep = dzInv; \ if (k == 0) \ { \ lz = 0; \ nz = zSkip; \ zStep = dzInv; \ } \ else if (k == (zKnt - 1)) \ { \ lz = zSkip; \ nz = 0; \ zStep = dzInv; \ } \ else \ { \ lz = zSkip; \ nz = zSkip; \ zStep = 0.5 * dzInv; \ } \ \ last = data - lx; \ next = data + nx; \ \ dx[0] = ((float)*next++ - (float)*last++) * xStep; \ dx[1] = ((float)*next++ - (float)*last++) * xStep; \ dx[2] = ((float)*next - (float)*last ) * xStep; \ \ last = data - ly; \ next = data + ny; \ \ dy[0] = ((float)*next++ - (float)*last++) * yStep; \ dy[1] = ((float)*next++ - (float)*last++) * yStep; \ dy[2] = ((float)*next - (float)*last ) * yStep; \ \ last = data - lz; \ next = data + nz; \ \ dz[0] = ((float)*next++ - (float)*last++) * zStep; \ dz[1] = ((float)*next++ - (float)*last++) * zStep; \ dz[2] = ((float)*next - (float)*last ) * zStep; \ \ if (div) \ { \ REG_DIV(dx, dy, dz, div); \ div ++; \ } \ \ if (curl) \ { \ REG_CURL(dx, dy, dz, curl); \ curl += 3; \ } \ \ data += 3; \ } \ } \ } \ } switch(type) { case TYPE_DOUBLE: REG_CUBES_DIVCURL(double); break; case TYPE_FLOAT: REG_CUBES_DIVCURL(float); break; case TYPE_INT: REG_CUBES_DIVCURL(int); break; case TYPE_UINT: REG_CUBES_DIVCURL(uint); break; case TYPE_SHORT: REG_CUBES_DIVCURL(short); break; case TYPE_USHORT: REG_CUBES_DIVCURL(ushort); break; case TYPE_BYTE: REG_CUBES_DIVCURL(byte); break; case TYPE_UBYTE: REG_CUBES_DIVCURL(ubyte); break; default: DXSetError(ERROR_INVALID_DATA, "#10320", "data"); goto error; } return OK; error: return ERROR; } static Error _dxfDivCurlIrregular(Field field) { Array pA = NULL, cA = NULL, dA = NULL, curlA = NULL, divA = NULL; void (*elementMethod)() = NULL; Object attr; int pDim, nVperE, nPoints, nElements, eDim = 0; Type dataType; char *eltType; byte *counts = NULL; float *divs = NULL; float *curls = NULL; Pointer data = NULL; int *e, i, j; float *dPtr, *cPtr; byte *c; int curlsInLocal, divsInLocal; ArrayHandle pH = NULL, cH = NULL; InvalidComponentHandle icH = NULL; if (DXEmptyField(field)) goto done; /* * Make sure invalid connections are up to data */ if (! DXInvalidateConnections((Object)field)) goto error; /* * Make sure data is dep on positions */ attr = DXGetComponentAttribute(field, "data", "dep"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "#10255", "data", "dependency"); goto error; } if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10200", "dependency attribute"); goto error; } if (strcmp(DXGetString((String)attr), "positions")) { DXSetError(ERROR_INVALID_DATA, "#11251", "data"); goto error; } /* * Get the element type */ attr = DXGetComponentAttribute(field, "connections", "element type"); if (! attr) { DXSetError(ERROR_INVALID_DATA, "#10255", "connections", "element type"); goto error; } if (DXGetObjectClass(attr) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "#10200", "element type attribute"); goto error; } eltType = DXGetString((String)attr); /* * Get various arrays, pointers to their contents and other * related information */ pA = (Array)DXGetComponentValue(field, "positions"); if (! pA) { DXSetError(ERROR_MISSING_DATA, "#10240", "data", "positions"); goto error; } pH = DXCreateArrayHandle(pA); if (! pH) goto error; DXGetArrayInfo(pA, &nPoints, NULL, NULL, NULL, &pDim); if (DXGetError() != ERROR_NONE) goto error; cA = (Array)DXGetComponentValue(field, "connections"); if (! cA) { DXSetError(ERROR_MISSING_DATA, "#10240", "data", "connections"); goto error; } cH = DXCreateArrayHandle(cA); if (! cH) goto error; DXGetArrayInfo(cA, &nElements, NULL, NULL, NULL, &nVperE); if (DXGetError() != ERROR_NONE) goto error; dA = (Array)DXGetComponentValue(field, "data"); data = DXGetArrayDataLocal(dA); DXGetArrayInfo(dA, NULL, &dataType, NULL, NULL, NULL); /* * Get element-type dependent info: the dimensionality of the elements * and a pointer to a routine which will determine element div and curl */ if (! strcmp(eltType, "tetrahedra")) { eDim = 3; elementMethod = _dxfTetrasDivCurl; } else if (! strcmp(eltType, "cubes")) { eDim = 3; elementMethod = _dxfCubesDivCurl; } else { DXSetError(ERROR_NOT_IMPLEMENTED, "#11380", eltType); goto error; } /* * Dimensionality of elements and the space in which they are embedded * had better agree. */ if (eDim != pDim) { DXSetError(ERROR_INVALID_DATA, "#11360", pDim, eltType, eDim); goto error; } /* * Allocation a counts array, recalling the number of element * samples that have been accumulated for each position. */ counts = (byte *)DXAllocateLocalZero(nPoints); if (! counts) { DXResetError(); counts = (byte *)DXAllocateZero(nPoints); } if (! counts) goto error; /* * Ideally, use a local memory buffer for the output buffers * If insufficient memory is available, use the output array directly * to avoid an unnecessary copy. */ curlA = (Array)DXGetComponentValue(field, "curl"); if (curlA) { curls = (float *)DXAllocateLocalZero(nPoints * pDim * sizeof(float)); if (! curls) { DXResetError(); curls = (float *)DXGetArrayData(curlA); if (! curls || DXGetError()) goto error; curlsInLocal = FALSE; } else curlsInLocal = TRUE; } else curls = NULL; divA = (Array)DXGetComponentValue(field, "div"); if (divA) { divs = (float *)DXAllocateLocalZero(nPoints * sizeof(float)); if (! divs) { DXResetError(); divs = (float *)DXGetArrayData(divA); if (! divs || DXGetError()) goto error; divsInLocal = FALSE; } else divsInLocal = TRUE; } else divs = NULL; /* * Get the invalid component handle for connections */ icH = DXCreateInvalidComponentHandle((Object)field, NULL, "connections"); if (! icH) goto error; /* * For each element, compute the local redults and add the results * to the slots of the buffers that correspond to the vertices * of the element. Bump a counter for the affected vertices. */ for (i = 0; i < nElements; i++) { if (DXIsElementValid(icH, i)) { int ebuf[16]; e = DXGetArrayEntry(cH, i, (Pointer)ebuf); (*elementMethod)(pH, data, dataType, divs, curls, counts, e); } } DXFreeInvalidComponentHandle(icH); icH = NULL; /* * For each output datum, divide through by the number of elements * that contributed to it. */ dPtr = divs; cPtr = curls; c = counts; for (i = 0; i < nPoints; i++, c++) { if (*c) { float d = 1.0 / counts[i]; if (dPtr) *dPtr++ *= d; if (cPtr) for (j = 0; j < pDim; j++) *cPtr++ *= d; } else { if (dPtr) dPtr ++; if (cPtr) cPtr ++; } } DXFree((Pointer)counts); counts = NULL; /* * If we were able to create local-memory output buffers, we * now have to copy out the local buffers. * Otherwise, the arrays already exist and the data's already there. */ if (divA && divsInLocal) { if (! DXAddArrayData(divA, 0, nPoints, (Pointer)divs)) goto error; DXFree((Pointer)divs); divs = NULL; } if (curlA && curlsInLocal) { if (! DXAddArrayData(curlA, 0, nPoints, (Pointer)curls)) goto error; DXFree((Pointer)curls); curls = NULL; } DXFreeArrayHandle(pH); DXFreeArrayHandle(cH); if (dA && data) DXFreeArrayDataLocal(dA, (Pointer)data); done: return OK; error: DXFree((Pointer)counts); if (icH) DXFreeInvalidComponentHandle(icH); DXFreeArrayHandle(pH); DXFreeArrayHandle(cH); if (dA && data) DXFreeArrayDataLocal(dA, (Pointer)data); if (divA && divsInLocal) DXFreeArrayDataLocal(divA, (Pointer)divs); if (curlA && curlsInLocal) DXFreeArrayDataLocal(curlA, (Pointer)curls); return ERROR; } #define TETRA_WEIGHTS(x, y, z, w0, w1, w2, w3) \ { \ float x0t, y0t, z0t; \ float xt1, yt1, zt1; \ float xt2, yt2, zt2; \ float xt3, yt3, zt3; \ float yzt12, yzt13, yzt23, yz0t2, yz0t3, yz01t, yz02t; \ float wc, wx, wy, wz; \ \ x0t = (x) - x0; y0t = (y) - y0; z0t = (z) - z0; \ xt1 = x1 - (x); yt1 = y1 - (y); zt1 = z1 - (z); \ xt2 = x2 - (x); yt2 = y2 - (y); zt2 = z2 - (z); \ xt3 = x3 - (x); yt3 = y3 - (y); zt3 = z3 - (z); \ \ yzt12 = yt1 * zt2 - yt2 * zt1; \ yzt13 = yt1 * zt3 - yt3 * zt1; \ yzt23 = yt2 * zt3 - yt3 * zt2; \ yz0t2 = y0t * z02 - y02 * z0t; \ yz0t3 = y0t * z03 - y03 * z0t; \ yz01t = y01 * z0t - y0t * z01; \ yz02t = y02 * z0t - y0t * z02; \ \ w0 = (xt1 * yzt23 - xt2 * yzt13 + xt3 * yzt12) * invVolume; \ w1 = (x0t * yz023 - x02 * yz0t3 + x03 * yz0t2) * invVolume; \ w2 = (x01 * yz0t3 - x0t * yz013 + x03 * yz01t) * invVolume; \ w3 = (x01 * yz02t - x02 * yz01t + x0t * yz012) * invVolume; \ } #define TETRA_VECTOR(type, p, q, r, s, wp, wq, wr, ws, vector) \ { \ type *dp, *dq, *dr, *ds; \ \ /* \ * Get the data values at the vertices \ */ \ dp = ((type *)data) + 3*p; \ dq = ((type *)data) + 3*q; \ dr = ((type *)data) + 3*r; \ ds = ((type *)data) + 3*s; \ \ /* \ * DXInterpolate the data values for the \ * center point and the points a step away \ * along the three axes \ */ \ vector[0] = wp*(*dp++)+wq*(*dq++)+wr*(*dr++)+ws*(*ds++); \ vector[1] = wp*(*dp++)+wq*(*dq++)+wr*(*dr++)+ws*(*ds++); \ vector[2] = wp*(*dp) +wq*(*dq) +wr*(*dr) +ws*(*ds) ; \ } #define DELTA_VECTORS(vc, vx, vy, vz, dx, dy, dz) \ { \ dx[0]=vx[0]-vc[0]; dx[1]=vx[1]-vc[1]; dx[2]=vx[2]-vc[2]; \ dy[0]=vy[0]-vc[0]; dy[1]=vy[1]-vc[1]; dy[2]=vy[2]-vc[2]; \ dz[0]=vz[0]-vc[0]; dz[1]=vz[1]-vc[1]; dz[2]=vz[2]-vc[2]; \ } #define IRREG_CURL(dx, dy, dz, curl) \ { \ curl[0] = (dy[2] - dz[1]) * invStep; \ curl[1] = (dz[0] - dx[2]) * invStep; \ curl[2] = (dx[1] - dy[0]) * invStep; \ } #define IRREG_DIV(dx, dy, dz, div) \ { \ div = (dx[0] + dy[1] + dz[2]) * invStep; \ } static void _dxfTetrasDivCurl(ArrayHandle pH, Pointer data, Type type, float *divs, float *curls, byte *counts, int *tetra) { float *pt; int p, q, r, s; float x0, x1, x2, x3; float y0, y1, y2, y3; float z0, z1, z2, z3; float x01, x02, x03; float y01, y02, y03; float z01, z02, z03; float wc0, wc1, wc2, wc3; /* weights for 4 vertices for centroid */ float wx0, wx1, wx2, wx3; /* weights for 4 vertices for x step */ float wy0, wy1, wy2, wy3; /* weights for 4 vertices for y step */ float wz0, wz1, wz2, wz3; /* weights for 4 vertices for z step */ float volume, invVolume; float step, invStep; float yz012, yz013, yz023; float cx, cy, cz; /* centroid */ float curl[3], div; float vc[3], vx[3], vy[3], vz[3]; float dc[3], dx[3], dy[3], dz[3]; float *gPtr; float pbuf[3]; p = *tetra++; q = *tetra++; r = *tetra++; s = *tetra; pt = DXGetArrayEntry(pH, p, (Pointer)pbuf); x0 = *pt++; y0 = *pt++; z0 = *pt; pt = DXGetArrayEntry(pH, q, (Pointer)pbuf); x1 = *pt++; y1 = *pt++; z1 = *pt; pt = DXGetArrayEntry(pH, r, (Pointer)pbuf); x2 = *pt++; y2 = *pt++; z2 = *pt; pt = DXGetArrayEntry(pH, s, (Pointer)pbuf); x3 = *pt++; y3 = *pt++; z3 = *pt; x01 = x1 - x0; y01 = y1 - y0; z01 = z1 - z0; x02 = x2 - x0; y02 = y2 - y0; z02 = z2 - z0; x03 = x3 - x0; y03 = y3 - y0; z03 = z3 - z0; /* * Need to get a handle on the size of the tetrahedra. Use manhattan * length to estimate edge length; average three edges incident on p0 * to get an average edge length; take 0.01 of that. */ step = ((x01 + y01 + z01) + (x02 + y02 + z02) + (x03 + y03 + z03)) * 0.0333; if (step == 0.0) return; /* * Need to divide data delta by step length */ invStep = 1.0 / step; /* * A few intermediates we will be re-using */ yz012 = y01 * z02 - y02 * z01; yz013 = y01 * z03 - y03 * z01; yz023 = y02 * z03 - y03 * z02; /* * volume of tetrahedra */ volume = x01 * yz023 - x02 * yz013 + x03 * yz012; if (volume == 0.0) return; /* * Need to normalize weights later */ invVolume = 1.0 / volume; /* * Get a point inside the tetrahedra */ cx = (x0 + x1 + x2 + x3) * 0.25; cy = (y0 + y1 + y2 + y3) * 0.25; cz = (z0 + z1 + z2 + z3) * 0.25; /* * evalate weights for center point and for * points one step up along each axis */ TETRA_WEIGHTS(cx, cy, cz, wc0, wc1, wc2, wc3); TETRA_WEIGHTS(cx+step, cy, cz, wx0, wx1, wx2, wx3); TETRA_WEIGHTS(cx, cy+step, cz, wy0, wy1, wy2, wy3); TETRA_WEIGHTS(cx, cy, cz+step, wz0, wz1, wz2, wz3); switch(type) { case TYPE_DOUBLE: TETRA_VECTOR(double, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(double, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(double, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(double, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_FLOAT: TETRA_VECTOR(float, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(float, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(float, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(float, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_INT: TETRA_VECTOR(int, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(int, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(int, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(int, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_UINT: TETRA_VECTOR(uint, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(uint, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(uint, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(uint, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_SHORT: TETRA_VECTOR(short, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(short, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(short, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(short, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_USHORT: TETRA_VECTOR(ushort, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(ushort, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(ushort, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(ushort, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_BYTE: TETRA_VECTOR(byte, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(byte, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(byte, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(byte, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; case TYPE_UBYTE: TETRA_VECTOR(ubyte, p, q, r, s, wc0, wc1, wc2, wc3, vc); TETRA_VECTOR(ubyte, p, q, r, s, wx0, wx1, wx2, wx3, vx); TETRA_VECTOR(ubyte, p, q, r, s, wy0, wy1, wy2, wy3, vy); TETRA_VECTOR(ubyte, p, q, r, s, wz0, wz1, wz2, wz3, vz); break; default: return; } DELTA_VECTORS(vc, vx, vy, vz, dx, dy, dz); if (curls) IRREG_CURL(dx, dy, dz, curl); if (divs) IRREG_DIV(dx, dy, dz, div); if (divs) { divs[p] += div; divs[q] += div; divs[r] += div; divs[s] += div; } if (curls) { float *cPtr; cPtr = curls + 3*p; *cPtr++ += curl[0]; *cPtr++ += curl[1]; *cPtr += curl[2]; cPtr = curls + 3*q; *cPtr++ += curl[0]; *cPtr++ += curl[1]; *cPtr += curl[2]; cPtr = curls + 3*r; *cPtr++ += curl[0]; *cPtr++ += curl[1]; *cPtr += curl[2]; cPtr = curls + 3*s; *cPtr++ += curl[0]; *cPtr++ += curl[1]; *cPtr += curl[2]; } counts[p] ++; counts[q] ++; counts[r] ++; counts[s] ++; return; } static void _dxfCubesDivCurl(ArrayHandle pH, Pointer data, Type type, float *div, float *curl, byte *counts, int *cube) { int tetra[4]; tetra[0] = cube[0]; tetra[1] = cube[6]; tetra[2] = cube[5]; tetra[3] = cube[3]; _dxfTetrasDivCurl(pH, data, type, div, curl, counts, tetra); tetra[0] = cube[1]; tetra[1] = cube[2]; tetra[2] = cube[4]; tetra[3] = cube[7]; _dxfTetrasDivCurl(pH, data, type, div, curl, counts, tetra); }