/***********************************************************************/ /* Open Visualization Data Explorer */ /* (C) Copyright IBM Corp. 1989,1999 */ /* ALL RIGHTS RESERVED */ /* This code licensed under the */ /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */ /***********************************************************************/ /* * $Header: /home/gda/dxcvs/dx/src/exec/dxmods/streamline.c,v 1.3 1999/05/10 15:45:32 gda Exp $: */ #include #include "stream.h" #include "vectors.h" #include "math.h" #include extern VectorGrp _dxfIrreg_InitVectorGrp(Object, char *); extern VectorGrp _dxfReg_InitVectorGrp(Object, char *); typedef enum { STREAM_OK, STREAM_FULL, STREAM_ERROR } StreamFlags; #ifdef STREAM_MAX #undef STREAM_MAX #endif #define STREAM_MAX 25000 #define TOT_MAX 50000 #define ZERO_MAX 100 static Stream NewStream(int); static void DestroyStream(Stream); static Error FlushStream(Stream); static Field StreamToField(Stream); static StreamFlags AddPointToStream(Stream, POINT_TYPE *, VECTOR_TYPE *, double); static VectorField InitVectorField(Object); static Error DestroyVectorField(VectorField); static Error GetElementType(Object, char **); static Array Starts(Object, Object, int); static Error StartsRecurse(Object, Array *); static Error StartsAddPoints(Array *, Array); static Error AlignStartPtsAndTimes(Array, Array); static Error Streamlines(Object, Array, Array, float, float, Group, Object, int); static Error Streamline(VectorField, float *, float, float, float, Group); static Error StreamTask(Pointer); static Error TraceFrame(Pointer); static Error InitFrame(float *, float *, float *); static Error UpdateFrame(float *, float *, float *, float *, float *, float *, float *, float *); static void RotateAroundVector(Vector *, float, float, float *); static int ZeroVector(VECTOR_TYPE *, int nDim); static double VectorDot(VECTOR_TYPE *, VECTOR_TYPE *, int nDim); static int IsRegular(Object); static Error GeometryCheck(Object, int); static InstanceVars FindElement(VectorField, POINT_TYPE *, VectorGrp); static InstanceVars FindMultiGridContinuation(VectorField, POINT_TYPE *, VectorGrp); m_Streamline(Object *in, Object *out) { Array starts = NULL, times = NULL; Object vectors, curls = NULL; int curlFlag = 0; float c; float duration; Group group = NULL; Type type; Category cat; int rank; int nD, n; Class vClass; out[0] = NULL; if (! in[0]) { DXSetError(ERROR_BAD_PARAMETER, "#10000", "data"); goto error; } vectors = in[0]; if ((vClass = DXGetObjectClass(vectors)) == CLASS_GROUP) vClass = DXGetGroupClass((Group)vectors); if (vClass != CLASS_FIELD && vClass != CLASS_COMPOSITEFIELD && vClass != CLASS_MULTIGRID) { DXSetError(ERROR_BAD_PARAMETER, "data must be a field, multigrid, or composite field"); goto error; } if (! DXGetType(vectors, &type, &cat, &rank, &nD)) { DXSetError(ERROR_MISSING_DATA, "#10250", "data", "data"); goto error; } if (type != TYPE_FLOAT || cat != CATEGORY_REAL || (nD != 2 && nD != 3)) { DXSetError(ERROR_BAD_PARAMETER, "#10302", "data"); goto error; } if (vectors) if (! GeometryCheck(vectors, nD)) goto error; /* * Starting points are gathered from arg 1 or centerpoint of vector * field. */ starts = Starts(in[1], in[0], nD); if (starts) DXGetArrayInfo(starts, &n, NULL, NULL, NULL, NULL); if (! starts || n == 0) goto no_streams; if (! DXTypeCheck(starts, TYPE_FLOAT, CATEGORY_REAL, 1, nD)) { DXSetError(ERROR_BAD_PARAMETER, "#11823", "start", "data"); goto error; } /* * Starting times are gathered from arg 2 or are a single 0.0 */ times = Starts(in[2], NULL, 1); /* * Match up the starting points and times */ if (! AlignStartPtsAndTimes(starts, times)) goto error; if (in[3]) { if (! DXExtractFloat(in[3], &duration)) { DXSetError(ERROR_BAD_PARAMETER, "#10080", "head"); goto error; } } else duration = DXD_MAX_FLOAT; if (in[4]) { int nDC; curls = in[4]; if (! DXGetType(curls, &type, &cat, &rank, &nDC)) goto error; if (type != TYPE_FLOAT || cat != CATEGORY_REAL || rank != 1) { DXSetError(ERROR_BAD_PARAMETER, "#10282", "curl"); goto error; } if (nDC != nD) { DXSetError(ERROR_BAD_PARAMETER, "#11400", "curl", "data"); goto error; } curlFlag = 1; } if (in[5]) { if (! DXExtractInteger(in[5], &curlFlag) || (curlFlag < 0 || curlFlag > 1)) { DXSetError(ERROR_BAD_PARAMETER, "#10070", "curlflag"); goto error; } } if (in[6]) { if (! DXExtractFloat(in[6], &c)) { DXSetError(ERROR_BAD_PARAMETER, "#10090", "stepscale"); goto error; } if (c <= 0.0) c = DEFAULT_C; } else c = DEFAULT_C; group = DXNewGroup(); if (! group) goto error; if (! Streamlines(vectors, starts, times, c, duration, group, curls, curlFlag)) goto error; if (! DXGetEnumeratedMember(group, 0, 0)) { DXDelete((Object)group); out[0] = (Object)DXEndField(DXNewField()); } else { DXSetFloatAttribute((Object)group, "fuzz", 4); out[0] = (Object)group; } DXDelete((Object)starts); DXDelete((Object)times); return OK; no_streams: out[0] = NULL; DXDelete((Object)starts); DXDelete((Object)times); return OK; error: DXDelete((Object)starts); DXDelete((Object)times); DXDelete((Object)group); return ERROR; } typedef struct task { VectorField vf; float point[3]; float time; float c, duration; Group g; } Task; static Error Streamlines(Object iVectors, Array starts, Array times, float c, float d, Group g, Object curl, int curlFlag) { Object vectors = NULL; VectorField vf = NULL; int i, j, n; Task task; float *point; float *time; int nDim; Field f; int nP; nP = DXProcessors(0); DXGetArrayInfo(starts, &n, NULL, NULL, NULL, &nDim); if (n < 1) return OK; vectors = DXCopy(iVectors, COPY_STRUCTURE); if (! vectors) return ERROR; vf = InitVectorField(vectors); if (! vf) goto error; point = (float *)DXGetArrayData(starts); time = (float *)DXGetArrayData(times); if (!point || !time) goto error; if (n == 1 || nP == 1) { for (i = 0; i < n; i++) { if (*time <= d) if (! Streamline(vf, point, *time, c, d, g)) goto error; point += nDim; time += 1; } } else { task.vf = vf; task.c = c; task.duration = d; task.g = g; if (! DXCreateTaskGroup()) goto error; for (i = 0; i < n; i++) { for (j = 0; j < nDim; j++) task.point[j] = *point++; task.time = *time++; if (task.time <= task.duration) if (! DXAddTask(StreamTask, (Pointer)&task, sizeof(task), 1.0)) { DXAbortTaskGroup(); goto error; } } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; } /* * Did we generate any streamlines? Count them. */ for (n = 0; DXGetEnumeratedMember(g, n, 0); n++); if (n <= 0) goto no_streams; if (vf->nDim == 3) { if (curlFlag) { Interpolator curlMap; if (curl) curlMap = DXNewInterpolator(curl, INTERP_INIT_PARALLEL, 0.0); else { MultiGrid mg = DXNewMultiGrid(); if (! mg) goto error; for (i = 0; i < vf->nmembers; i++) if (! (*(vf->members[i]->CurlMap))(vf->members[i], mg)) { DXDelete((Object)mg); goto error; } curlMap = DXNewInterpolator((Object)mg, INTERP_INIT_PARALLEL, -1.0); } if (! curlMap && DXGetError() != ERROR_NONE) goto error; if (curlMap) { if (! DXMap((Object)g, (Object)curlMap, "positions", "curl")) { DXDelete((Object)curlMap); goto error; } DXDelete((Object)curlMap); } } if (n == 1 || nP == 1) { for (i = 0; i < n; i++) { f = (Field)DXGetEnumeratedMember(g, i, NULL); if (! TraceFrame((Pointer)f)) goto error; } } else { if (! DXCreateTaskGroup()) goto error; for (i = 0; i < n; i++) { f = (Field)DXGetEnumeratedMember(g, i, NULL); if (! DXAddTask(TraceFrame, (Pointer)f, 0, 1.0)) { DXAbortTaskGroup(); goto error; } } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; } } no_streams: DestroyVectorField(vf); vf = NULL; DXDelete((Object)vectors); return OK; error: DXDelete((Object)vectors); DestroyVectorField(vf); return ERROR; } static Error StreamTask(Pointer p) { Task *t = (Task *)p; return Streamline(t->vf, t->point, t->time, t->c, t->duration, t->g); } static Error Streamline(VectorField vf, /* Vector field to trace streamline in */ float *point, /* streamline origination position */ float time, /* streamline origination time */ float c, /* step criterion */ float duration,/* max. lifetime of stream */ Group g) /* group in which to put the result */ { Stream stream = NULL; VectorGrp vg; POINT_TYPE pb0[3], pb1[3], *ptmp; VECTOR_TYPE vb0[3], vb1[3], *vtmp; POINT_TYPE *p0, *p1; VECTOR_TYPE *v0, *v1; double t, t2, elapsedtime; int nDim; int done, keepPoint; Error (*Interpolate)(); Error (*StepTime)(); Error (*FindBoundary)(); int (*Neighbor)(); int (*Weights)(); int (*FaceWeights)(); Field field = NULL; int i, totKnt, zeroKnt, streamKnt; InstanceVars I = NULL; StreamFlags sFlag; int zero; stream = NULL; elapsedtime = time; /* * Double buffer pointers */ p0 = pb0; p1 = pb1; v0 = vb0; v1 = vb1; /* * Initial conditions */ for (i = 0; i < vf->nDim; i++) p0[i] = point[i]; I = FindElement(vf, p0, NULL); if (! I) goto stream_done; vg = I->currentVectorGrp; Interpolate = vg->Interpolate; StepTime = vg->StepTime; FindBoundary = vg->FindBoundary; Neighbor = vg->Neighbor; Weights = vg->Weights; FaceWeights = vg->FaceWeights; nDim = vg->nDim; if (! (*Interpolate)(I, p0, v0)) goto error; stream = NewStream(nDim); if (! stream) goto error; sFlag = AddPointToStream(stream, p0, v0, elapsedtime); if (sFlag == STREAM_ERROR) goto error; if (ZeroVector(v0, nDim)) goto stream_done; /* * While stream is passing from partition to partition... */ done = 0; totKnt = 1; zeroKnt = 0; streamKnt = 1; while (!done) { /* * Get the point and the appropriate time step */ if (! (*StepTime)(I, c, v0, &t)) goto error; /* * If the time step exceeds the duration, truncate it and mark * the stream complete. */ if (elapsedtime + t > duration) { t = duration - elapsedtime; done = 1; } /* * Get the point a time step away in the direction of the * vector */ for (i = 0; i < nDim; i++) p1[i] = p0[i] + t*v0[i]; if (! (*Weights)(I, p1)) { int found; /* * Segment proceeds from a point inside an element (or on * the boundary of the element) to one outside the element. * Move to the boundary. */ if (! (*FindBoundary)(I, p0, v0, &t2)) goto stream_done; for (i = 0; i < nDim; i++) p1[i] = p0[i] + t2*v0[i]; /* * P1 is the point in the boundary. Interpolate a vector for it * based only on the face. */ if (! (*FaceWeights)(I, p1)) goto stream_done; if (! (*Interpolate)(I, p1, v1)) goto error; found = (*Neighbor)(I, v1); if (found == -1) goto error; else if (! found) if (vf->nmembers != 1) { VectorGrp lastGrp = I->currentVectorGrp; (*vg->FreeInstanceVars)(I); I = FindMultiGridContinuation(vf, p1, lastGrp); if (! I) { done = 1; } else { vg = I->currentVectorGrp; Interpolate = vg->Interpolate; StepTime = vg->StepTime; FindBoundary = vg->FindBoundary; Neighbor = vg->Neighbor; Weights = vg->Weights; FaceWeights = vg->FaceWeights; } } else done = 1; if (done || t2 > 0.05*t) keepPoint = 1; else keepPoint = 0; t = t2; } else { /* * Point is inside the current element. Regardless of whether * p0 is on the boundary, the step is fine and since p1 is not * on the boundary, we clear the boundary flag. */ if (! (*Interpolate)(I, p1, v1)) goto error; keepPoint = 1; } if (ZeroVector(v1, nDim) || VectorDot(v0, v1, nDim) < -0.0) done = 1; elapsedtime += t; if (keepPoint && t != 0.0) { sFlag = AddPointToStream(stream, p1, v1, elapsedtime); if (sFlag == STREAM_ERROR) goto error; else if (sFlag == STREAM_FULL) { DXWarning("streamline point limit exceeded"); done = 1; } streamKnt ++; } zero = 1; for (i = 0; zero && i < nDim; i++) zero = (p0[i] == p1[i]); if (zero) { if (++zeroKnt > ZERO_MAX) { DXWarning("possible infinite loop detected"); done = 1; } } else zeroKnt = 0; if (++totKnt > TOT_MAX) { DXWarning("possible infinite loop detected"); done = 1; } ptmp = p0; p0 = p1; p1 = ptmp; vtmp = v0; v0 = v1; v1 = vtmp; } stream_done: if (stream) { field = StreamToField(stream); if (! field) goto error; DestroyStream(stream); stream = NULL; if (! DXSetMember(g, NULL, (Object)field)) goto error; } if (I) (*vg->FreeInstanceVars)(I); return OK; error: if (I) (*vg->FreeInstanceVars)(I); DestroyStream(stream); DXDelete((Object)field); return ERROR; } static Error TraceFrame(Pointer ptr) { Field field; Array nArray = NULL, bArray = NULL; Array cArray, vArray, pArray; float *points, *vectors, *curls, *normals, *binormals; int i, nVectors, nDim; Object depattr = NULL; float fn[3], fb[3], ft[3]; field = (Field)ptr; if (DXEmptyField(field)) return OK; vArray = (Array)DXGetComponentValue(field, "data"); if (! vArray) goto error; DXGetArrayInfo(vArray, &nVectors, NULL, NULL, NULL, &nDim); if (nDim < 3) goto done; vectors = (float *)DXGetArrayData(vArray); if (! vectors) goto error; pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) goto error; DXGetArrayInfo(pArray, &nVectors, NULL, NULL, NULL, &nDim); if (nDim < 3) goto done; points = (float *)DXGetArrayData(pArray); if (! points) goto error; cArray = (Array)DXGetComponentValue(field, "curl"); if (cArray) { curls = (float *)DXGetArrayData(cArray); if (! curls) goto error; } else curls = NULL; nArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! nArray) goto error; if (! DXAddArrayData(nArray, 0, nVectors, NULL)) goto error; normals = (float *)DXGetArrayData(nArray); if (! normals) goto error; bArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3); if (! bArray) goto error; if (! DXAddArrayData(bArray, 0, nVectors, NULL)) goto error; binormals = (float *)DXGetArrayData(bArray); if (! binormals) goto error; if (! InitFrame(vectors, normals, binormals)) goto error; fn[0] = normals[0]; fn[1] = normals[1]; fn[2] = normals[2]; fb[0] = binormals[0]; fb[1] = binormals[1]; fb[2] = binormals[2]; ft[0] = vectors[0]; ft[1] = vectors[1]; ft[2] = vectors[2]; _dxfvector_normalize_3D((Vector *)fn, (Vector *)fn); _dxfvector_normalize_3D((Vector *)fb, (Vector *)fb); _dxfvector_normalize_3D((Vector *)ft, (Vector *)ft); for (i = 1; i < nVectors; i++) { float avcurl[3]; points += 3; vectors += 3; normals += 3; binormals += 3; if (curls) { avcurl[0] = 0.5 * (curls[0] + curls[3]); avcurl[1] = 0.5 * (curls[1] + curls[4]); avcurl[2] = 0.5 * (curls[2] + curls[5]); curls += 3; if (! UpdateFrame(points, vectors, avcurl, normals, binormals, fn, fb, ft)) goto error; } else { if (! UpdateFrame(points, vectors, NULL, normals, binormals, fn, fb, ft)) goto error; } } depattr = (Object)DXNewString("positions"); if (! depattr) goto error; DXReference(depattr); if (! DXSetComponentValue(field, "normals", (Object)nArray)) goto error; nArray = NULL; if (! DXSetComponentAttribute(field, "normals", "dep", depattr)) goto error; if (! DXSetComponentValue(field, "binormals", (Object)bArray)) goto error; bArray = NULL; if (! DXSetComponentAttribute(field, "binormals", "dep", depattr)) goto error; if (cArray) DXDeleteComponent(field, "curl"); DXDelete(depattr); done: return OK; error: DXDelete((Object)nArray); DXDelete((Object)bArray); return ERROR; } static Array Starts(Object starts, Object vectors, int nD) { Array array = NULL; if (starts) { Object xstarts; xstarts = DXApplyTransform(starts, NULL); if (! xstarts) goto error; if (! StartsRecurse(xstarts, &array)) { if (starts != xstarts) DXDelete((Object)xstarts); goto error; } if (starts != xstarts) DXDelete((Object)xstarts); } else if (vectors) { Point box[8]; float pt[3]; if (! DXBoundingBox(vectors, box)) return NULL; pt[0] = (box[0].x + box[7].x) / 2.0; pt[1] = (box[0].y + box[7].y) / 2.0; pt[2] = (box[0].z + box[7].z) / 2.0; array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nD); if (! array) goto error; if (! DXAddArrayData(array, 0, 1, (Pointer)pt)) goto error; } else { float *ptr; int i; array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nD); if (! array) goto error; if (! DXAddArrayData(array, 0, 1, NULL)) goto error; ptr = (float *)DXGetArrayData(array); if (! ptr) goto error; for (i = 0; i < nD; i++) *ptr++ = 0.0; } return array; error: DXDelete((Object)array); return NULL; } static Error AlignStartPtsAndTimes(Array starts, Array times) { int nStarts, nTimes; DXGetArrayInfo(starts, &nStarts, NULL, NULL, NULL, NULL); DXGetArrayInfo(times, &nTimes, NULL, NULL, NULL, NULL); if (nStarts == nTimes) return OK; else if (nStarts == 1) { int itemSize = DXGetItemSize(starts); int i; unsigned char *src, *dst; if (! DXAddArrayData(starts, 1, nTimes-1, NULL)) goto error; src = (unsigned char *)DXGetArrayData(starts); dst = src + itemSize; for (i = 0; i < nTimes; i++, dst += itemSize) memcpy(dst, src, itemSize); } else if (nTimes == 1) { int itemSize = DXGetItemSize(times); int i; unsigned char *src, *dst; if (! DXAddArrayData(times, 1, nStarts-1, NULL)) goto error; src = (unsigned char *)DXGetArrayData(times); dst = src + itemSize; for (i = 0; i < nStarts; i++, dst += itemSize) memcpy(dst, src, itemSize); } else { DXSetError(ERROR_BAD_PARAMETER, "#11160", "start positions", nStarts, "start times", nTimes); goto error; } return OK; error: return ERROR; } static Error StartsRecurse(Object o, Array *a) { Object c; int i; if (DXGetObjectClass(o) == CLASS_ARRAY) return StartsAddPoints(a, (Array)o); else if (DXGetObjectClass(o) == CLASS_FIELD) { if (! DXEmptyField((Field)o)) return StartsAddPoints(a, (Array)DXGetComponentValue((Field)o, "positions")); else return OK; } else if (DXGetObjectClass(o) == CLASS_GROUP) { i = 0; while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL))) { if (! StartsRecurse(c, a)) goto error; } return OK; } else { DXSetError(ERROR_BAD_PARAMETER, "starts"); return ERROR; } error: DXDelete((Object)*a); *a = NULL; return ERROR; } static Error StartsAddPoints(Array *dst, Array src) { int nSrc, nDst, rank, shape[32]; if (! src || DXGetObjectClass((Object)src) != CLASS_ARRAY) return ERROR; if (*dst == NULL) { DXGetArrayInfo(src, &nSrc, NULL, NULL, &rank, shape); if (nSrc == 0) return OK; *dst = DXNewArrayV(TYPE_FLOAT, CATEGORY_REAL, rank, shape); if (! *dst) return ERROR; if (rank == 0) shape[0] = 1; if (! DXAddArrayData(*dst, 0, nSrc, NULL)) return ERROR; if (! DXExtractParameter((Object)src, TYPE_FLOAT, shape[0], nSrc, DXGetArrayData(*dst))) return ERROR; } else { DXGetArrayInfo(*dst, &nDst, NULL, NULL, &rank, shape); DXGetArrayInfo(src, &nSrc, NULL, NULL, NULL, NULL); if (nSrc == 0) return OK; if (! DXAddArrayData(*dst, nDst, nSrc, NULL)) return ERROR; if (rank == 0) shape[0] = 1; if (! DXExtractParameter((Object)src, TYPE_FLOAT, shape[0], nSrc, (Pointer)(((float *)DXGetArrayData(*dst))+nDst*shape[0]))) return ERROR; } return OK; } static Error DestroyVectorField(VectorField vf) { if (vf) { int i; for (i = 0; i < vf->nmembers; i++) if (vf->members[i]) (*(vf->members[i]->Delete))(vf->members[i]); DXFree((Pointer)vf->members); DXFree((Pointer)vf); } } static VectorField InitVectorField(Object vfo) { char *str; VectorField vf = NULL; Class class; vf = (struct vectorField *)DXAllocateZero(sizeof(struct vectorField)); if (! vf) goto error; class = DXGetObjectClass(vfo); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)vfo); switch(class) { case CLASS_MULTIGRID: { int i, n; char *str; if (! DXGetMemberCount((Group)vfo, &n)) goto error; vf->nmembers = n; vf->members = (VectorGrp *)DXAllocate(n*sizeof(VectorGrp)); if (! vf->members) goto error; for (i = 0; i < n; i++) { Object mo = DXGetEnumeratedMember((Group)vfo, i, NULL); if (!mo) goto error; class = DXGetObjectClass(mo); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)mo); if (class != CLASS_FIELD && class != CLASS_COMPOSITEFIELD) { DXSetError(ERROR_INVALID_DATA, "invalid multigrid member"); goto error; } if (! GetElementType(mo, &str)) { DXSetError(ERROR_INVALID_DATA, "vectors"); goto error; } switch (IsRegular(mo)) { case 1: vf->members[i] = _dxfReg_InitVectorGrp(mo, str); break; case 0: vf->members[i] = _dxfIrreg_InitVectorGrp(mo, str); break; default: { DXSetError(ERROR_INVALID_DATA, "vectors"); goto error; } } if (! vf->members[i]) goto error; } } break; case CLASS_FIELD: case CLASS_COMPOSITEFIELD: { vf->nmembers = 1; vf->members = (VectorGrp *)DXAllocate(sizeof(VectorGrp)); if (! vf->members) goto error; if (! GetElementType(vfo, &str)) { DXSetError(ERROR_INVALID_DATA, "vectors"); goto error; } switch (IsRegular(vfo)) { case 1: vf->members[0] = _dxfReg_InitVectorGrp(vfo, str); break; case 0: vf->members[0] = _dxfIrreg_InitVectorGrp(vfo, str); break; default: { DXSetError(ERROR_INVALID_DATA, "vectors"); goto error; } } if (! vf->members[0]) goto error; } break; default: { DXSetError(ERROR_INVALID_DATA, "invalid object encountered in vector field"); goto error; } } vf->nDim = vf->members[0]->nDim; return vf; error: DestroyVectorField(vf); return NULL; } static Error GetElementType(Object o, char **str) { Object c; *str = NULL; if (DXGetObjectClass(o) == CLASS_FIELD) { if (DXEmptyField((Field)o)) return ERROR; if (! DXGetComponentValue((Field)o, "connections")) { DXSetError(ERROR_INVALID_DATA, "data has no connections"); return ERROR; } c = DXGetComponentAttribute((Field)o, "connections", "element type"); if (! c) { DXSetError(ERROR_MISSING_DATA, "element type attribute"); return ERROR; } if (DXGetObjectClass(c) != CLASS_STRING) { DXSetError(ERROR_INVALID_DATA, "element type attribute"); return ERROR; } *str = DXGetString((String)c); return OK; } else if (DXGetObjectClass(o) == CLASS_GROUP) { int i; Error j = ERROR; i = 0; while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL))) { j = GetElementType(c, str); if (j != ERROR) return OK; } return ERROR; } else { DXSetError(ERROR_BAD_CLASS, "vector field"); return ERROR; } } static Stream NewStream(nDim) { Stream s = NULL; s = (Stream)DXAllocateZero(sizeof(struct stream)); if (! s) goto error; s->points = (float *)DXAllocate(nDim*STREAM_BUF_POINTS*sizeof(float)); s->pArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nDim); s->time = (float *)DXAllocate(STREAM_BUF_POINTS*sizeof(float)); s->tArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0); s->vectors = (float *)DXAllocate(nDim*STREAM_BUF_POINTS*sizeof(float)); s->vArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nDim); if (s->points == NULL || s->time == NULL || s->pArray == NULL || s->tArray == NULL || s->vectors == NULL || s->vArray == NULL) goto error; s->nDim = nDim; s->bufKnt = 0; s->arrayKnt = 0; return s; error: DestroyStream(s); return NULL; } static void DestroyStream(Stream s) { if (s) { DXFree((Pointer)s->points); DXFree((Pointer)s->time); DXFree((Pointer)s->vectors); DXDelete((Object)s->pArray); DXDelete((Object)s->vArray); DXDelete((Object)s->tArray); DXFree((Pointer)s); } } static Error FlushStream(Stream s) { if (! DXAddArrayData(s->pArray, s->arrayKnt, s->bufKnt, (Pointer)s->points)) goto error; if (! DXAddArrayData(s->vArray, s->arrayKnt, s->bufKnt, (Pointer)s->vectors)) goto error; if (! DXAddArrayData(s->tArray, s->arrayKnt, s->bufKnt, (Pointer)s->time)) goto error; s->arrayKnt += s->bufKnt; s->bufKnt = 0; return OK; error: return ERROR; } static StreamFlags AddPointToStream(Stream s, POINT_TYPE *p, VECTOR_TYPE *v, double t) { if (s->bufKnt == STREAM_BUF_POINTS) if (! FlushStream(s)) goto error; if (s->nDim == 2) { float *ptr; ptr = s->points + s->bufKnt*s->nDim; *ptr++ = *p++; *ptr++ = *p++; ptr = s->vectors + s->bufKnt*s->nDim; *ptr++ = *v++; *ptr++ = *v++; } else { float *ptr; ptr = s->points + s->bufKnt*s->nDim; *ptr++ = *p++; *ptr++ = *p++; *ptr++ = *p++; ptr = s->vectors + s->bufKnt*s->nDim; *ptr++ = *v++; *ptr++ = *v++; *ptr++ = *v++; } s->time[s->bufKnt] = t; s->bufKnt ++; if (s->arrayKnt + s->bufKnt > STREAM_MAX) return STREAM_FULL; return STREAM_OK; error: return STREAM_ERROR; } static Field StreamToField(Stream s) { float c[3]; float d[3]; Field f = NULL; Array a = NULL; Object depattr = NULL; int flag = 1; depattr = (Object)DXNewString("positions"); if (! depattr) goto error; DXReference(depattr); c[0] = 0.7; c[1] = 0.7; c[2] = 0.0; d[0] = 0.0; d[1] = 0.0; d[2] = 0.0; if (! FlushStream(s)) goto error; f = DXNewField(); if (! f) goto error; if (! DXSetComponentValue(f, "positions", (Object)s->pArray)) goto error; s->pArray = NULL; if (! DXSetComponentAttribute(f, "positions", "dep", depattr)) goto error; if (! DXSetComponentValue(f, "data", (Object)s->vArray)) goto error; s->vArray = NULL; if (! DXSetComponentAttribute(f, "data", "dep", depattr)) goto error; if (! DXSetComponentValue(f, "time", (Object)s->tArray)) goto error; s->tArray = NULL; if (! DXSetComponentAttribute(f, "time", "dep", depattr)) goto error; a = DXMakeGridConnections(1, s->arrayKnt); if (! a) goto error; if (! DXSetComponentValue(f, "connections", (Object)a)) goto error; a = NULL; a = (Array)DXNewRegularArray(TYPE_FLOAT, 3, s->arrayKnt, (Pointer)c, (Pointer)d); if (! a) goto error; if (! DXSetComponentValue(f, "colors", (Object)a)) goto error; a = NULL; if (! DXSetComponentAttribute(f, "connections", "element type", (Object)DXNewString("lines"))) goto error; if (! DXSetIntegerAttribute((Object)f, "shade", 0)) goto error; if (! DXEndField(f)) goto error; DXDelete(depattr); return f; error: DXDelete(depattr); DXDelete((Object)f); return ERROR; } Error _dxfMinMaxBox(Field f, float *min, float *max) { float box[24], *b, fuzz; int i, j; if (! DXBoundingBox((Object)f, (Point *)box)) return ERROR; min[0] = min[1] = min[2] = DXD_MAX_FLOAT; max[0] = max[1] = max[2] = -DXD_MAX_FLOAT; b = box; for (i = 0; i < 8; i++) for (j = 0; j < 3; j++, b++) { if (*b < min[j]) min[j] = *b; if (*b > max[j]) max[j] = *b; } return OK; } int _dxfIsInBox(VectorPart p, POINT_TYPE *pt, int n, int fuzz) { int i; float f; if (fuzz) { for (i = 0; i < n; i++) { float f = 0.0001 * (p->max[i] - p->min[i]); if (pt[i] < (p->min[i]-f) || pt[i] > (p->max[i]+f)) return 0; } } else { for (i = 0; i < n; i++) if (pt[i] < p->min[i] || pt[i] > p->max[i]) return 0; } return 1; } static int ZeroVector(VECTOR_TYPE *v, int n) { int i; for (i = 0; i < n; i++) if (v[i] != 0.0) return 0; return 1; } static double VectorDot(VECTOR_TYPE *v0, VECTOR_TYPE *v1, int n) { int i; double d = 0; for (i = 0; i < n; i++) d += *v0++ * *v1++; return d; } static Error InitFrame(float *v, float *n, float *b) { n[0] = n[1] = n[2] = 0.0; if (v[0] == 0.0 && v[1] == 0.0 && v[2] == 0.0) { n[0] = 0.0; n[1] = 0.0; n[2] = 0.0; b[0] = 0.0; b[1] = 0.0; b[2] = 0.0; return OK; } if (v[0] == 0.0 && v[1] == 0.0 && v[2] > 0.0) n[0] = 1.0; else if (v[0] == 0.0 && v[1] == 0.0 && v[2] < 0.0) n[0] = -1.0; else if (v[0] == 0.0 && v[1] > 0.0 && v[2] == 0.0) n[2] = 1.0; else if (v[0] == 0.0 && v[1] < 0.0 && v[2] == 0.0) n[2] = -1.0; else if (v[0] > 0.0 && v[1] == 0.0 && v[2] == 0.0) n[1] = 1.0; else if (v[0] < 0.0 && v[1] == 0.0 && v[2] == 0.0) n[1] = -1.0; else { float d; if (v[2] != 0.0) { n[0] = n[1] = 1.0; n[2] = -(v[0] + v[1]) / v[2]; } else if (v[1] != 0.0) { n[0] = n[2] = 1.0; n[1] = -(v[0] + v[2]) / v[1]; } d = 1.0 / sqrt(n[0]*n[0] + n[1]*n[1] + n[2]*n[2]); n[0] *= d; n[1] *= d; n[2] *= d; } _dxfvector_cross_3D((Vector *)n, (Vector *)v, (Vector *)b); _dxfvector_normalize_3D((Vector *)b, (Vector *)b); return OK; } /* * update frame of reference from prior vertex to current vertex based * on incuming and outgoing edge vectors and twist. */ #define VecMat(a, b, c) \ { \ float x = *(a+0); \ float y = *(a+1); \ float z = *(a+2); \ \ *(c+0) = x*(b)[0] + y*(b)[3] + z*(b)[6]; \ *(c+1) = x*(b)[1] + y*(b)[4] + z*(b)[7]; \ *(c+2) = x*(b)[2] + y*(b)[5] + z*(b)[8]; \ } #define MatMat(a, b, c) \ { \ VecMat((a+0), (b), (c+0)); \ VecMat((a+3), (b), (c+3)); \ VecMat((a+6), (b), (c+6)); \ } static Error UpdateFrame(float *p, float *v, float *c, float *n, float *b, float *fn, float *fb, float *ft) { float twist; float mTwist[9], mBend[9], mProduct[9]; float cross[3], len, cA; int bend = 0; float nt[3]; /* * If theres a curl, rotate by the curl around the prior segment vector. */ if (c != NULL) { float step[3]; _dxfvector_subtract_3D((Vector *)p, (Vector *)(p-3), (Vector *)step); /* _dxfvector_scale_3D((Vector *)step, (*t - *(t-1)), (Vector *)step); */ twist = 0.5 * _dxfvector_dot_3D((Vector *)step, (Vector *)c); if (twist != 0.0) { RotateAroundVector((Vector *)ft, cos(twist), sin(twist), mTwist); VecMat(fn, mTwist, fn); VecMat(fb, mTwist, fb); } } /* * Look at the length of the step proceeding from the current * vertex. If its zero length, just leave the frame of reference * alone. */ len = _dxfvector_length_3D((Vector *)v); if (len == 0.0) { n[0] = fn[0]; n[1] = fn[1]; n[2] = fn[2]; b[0] = fb[0]; b[1] = fb[1]; b[2] = fb[2]; return OK; } /* * DXDot the incoming and outgoing tangents to determine * whether a bend occurs. */ _dxfvector_scale_3D((Vector *)v, 1.0/len, (Vector *)nt); cA = _dxfvector_dot_3D((Vector *)nt, (Vector *)ft); /* * If there's a bend angle... */ if (cA < 0.999999) { float sA, angle; bend = 1; /* * DXNormalize the bend axis */ _dxfvector_cross_3D((Vector *)nt, (Vector *)ft, (Vector *)cross); _dxfvector_normalize_3D((Vector *)cross, (Vector *)cross); /* * Rotate the incoming frame of reference by half the angle to * get a vertex FOR */ angle = acos(cA)/2; cA = cos(angle); sA = sin(angle); RotateAroundVector((Vector *)cross, cA, sA, mBend); VecMat(fn, mBend, fn); VecMat(fb, mBend, fb); VecMat(ft, mBend, ft); /* * Now we have the vertex frame of reference */ n[0] = fn[0]; n[1] = fn[1]; n[2] = fn[2]; b[0] = fb[0]; b[1] = fb[1]; b[2] = fb[2]; /* * If there was a bend, perform the second half of the rotation */ if (bend) { VecMat(fn, mBend, fn); VecMat(fb, mBend, fb); } /* * Outgoing tangent is normalized exiting segment vector. Also * normalize the frame normal and binormal. */ _dxfvector_normalize_3D((Vector *)fn, (Vector *)fn); _dxfvector_normalize_3D((Vector *)fb, (Vector *)fb); ft[0] = nt[0]; ft[1] = nt[1]; ft[2] = nt[2]; } else { n[0] = fn[0]; n[1] = fn[1]; n[2] = fn[2]; b[0] = fb[0]; b[1] = fb[1]; b[2] = fb[2]; } return OK; error: return ERROR; } /* * Create a 3x3 matrix defined by the axis v and the sin and cosine * of an angle Theta. */ static void RotateAroundVector(Vector *v, float c, float s, float *M) { float x, y, z; float sx, sy, sz; x = v->x; y = v->y; z = v->z; sx = x*x; sy = y*y; sz = z*z; M[0] = sx*(1.0-c) + c; M[1] = x*y*(1.0-c)-z*s; M[2] = z*x*(1.0-c)+y*s; M[3] = x*y*(1.0-c)+z*s; M[4] = sy*(1.0-c) + c; M[5] = y*z*(1.0-c)-x*s; M[6] = z*x*(1.0-c)-y*s; M[7] = y*z*(1.0-c)+x*s; M[8] = sz*(1.0-c) + c; } static int IsRegular(Object o) { Object c; Array a, b; if (DXGetObjectClass(o) == CLASS_FIELD) { if (DXEmptyField((Field)o)) return -1; a = (Array)DXGetComponentValue((Field)o, "connections"); b = (Array)DXGetComponentValue((Field)o, "positions"); if (! a || ! b) return -1; if (DXQueryGridConnections(a, NULL, NULL) && DXQueryGridPositions(b, NULL, NULL, NULL, NULL)) return 1; else return 0; } else if (DXGetObjectClass(o) == CLASS_GROUP) { int i, j, k; i = 0; k = -1; while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL))) { j = IsRegular(c); if (j == 0) return 0; if (j == 1) k = 1; } return k; } else { DXSetError(ERROR_BAD_CLASS, "vector field"); return -1; } } static Error GeometryCheck(Object vectors, int nD) { Class class = DXGetObjectClass(vectors); if (class == CLASS_GROUP) class = DXGetGroupClass((Group)vectors); switch(class) { case CLASS_FIELD: { Field f; Array a; Object o; char *str; int n; f = (Field)vectors; if (DXEmptyField(f)) return OK; a = (Array)DXGetComponentValue(f, "positions"); if (! DXGetArrayInfo(a, NULL, NULL, NULL, NULL, &n)) goto error; if (n != nD) { DXSetError(ERROR_INVALID_DATA, "vectors are %d-D, vector space is %d-D", nD, n); goto error; } if (! DXGetComponentValue(f, "connections")) { DXSetError(ERROR_MISSING_DATA, "connections"); goto error; } o = DXGetComponentAttribute(f, "connections", "element type"); if (! o) { DXSetError(ERROR_MISSING_DATA, "connections element type"); goto error; } if (DXGetObjectClass(o) != CLASS_STRING) { DXSetError(ERROR_BAD_CLASS, "connections element type attribute"); goto error; } str = DXGetString((String)o); /* switch(n) { case 2: if (strcmp(str, "quads") && strcmp(str, "triangles")) { DXSetError(ERROR_INVALID_DATA, "%d-D vector space requires %d-D elements", n, n); goto error; } break; case 3: if (strcmp(str, "cubes") && strcmp(str, "tetrahedra")) { DXSetError(ERROR_INVALID_DATA, "%d-D vector space requires %d-D elements", n, n); goto error; } break; default: DXSetError(ERROR_INVALID_DATA, "2- or 3-D vectors required"); goto error; } */ break; } case CLASS_MULTIGRID: case CLASS_COMPOSITEFIELD: { Object child; int i; i = 0; while (NULL != (child = DXGetEnumeratedMember((Group)vectors, i++, NULL))) if (! GeometryCheck(child, nD)) return ERROR; break; } default: DXSetError(ERROR_INVALID_DATA, "vectors must be field or composite field"); goto error; } return OK; error: return ERROR; } static InstanceVars FindElement(VectorField vField, POINT_TYPE *point, VectorGrp last) { int i; InstanceVars I; for (i = 0; i < vField->nmembers; i++) { if (vField->members[i] != last) { I = (*(vField->members[i]->NewInstanceVars))(vField->members[i]); if ((*(I->currentVectorGrp->FindElement))(I, point)) { vField->current = vField->members[i]; return I; } (*(I->currentVectorGrp->FreeInstanceVars))(I); } } return 0; } static InstanceVars FindMultiGridContinuation(VectorField vField, POINT_TYPE *point, VectorGrp last) { int i; InstanceVars I; for (i = 0; i < vField->nmembers; i++) { if (vField->members[i] != last) { I = (*(vField->members[i]->NewInstanceVars))(vField->members[i]); if ((*(I->currentVectorGrp->FindMultiGridContinuation))(I, point)) { vField->current = vField->members[i]; return I; } (*(I->currentVectorGrp->FreeInstanceVars))(I); } } return 0; }