/***********************************************************************/ /* 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/streakline.c,v 1.3 1999/05/10 15:45:31 gda Exp $: */ #include #include "stream.h" #include "vectors.h" #include "math.h" #define NAME_ARG 0 #define DATA_ARG 1 #define POINTS_ARG 2 #define TIMES_ARG 3 #define HEAD_ARG 4 #define FRAME_ARG 5 #define CURL_ARG 6 #define FLAG_ARG 7 #define SCALE_ARG 8 #define HOLD_ARG 9 extern VectorGrp _dxfIrreg_InitVectorGrp(Object, char *); extern VectorGrp _dxfReg_InitVectorGrp(Object, char *); typedef enum { STREAK_OK, STREAK_FULL, STREAK_ERROR } StreakFlags; #ifdef STREAM_MAX #undef STREAM_MAX #endif #define STREAK_MAX 25000 #define TOT_MAX 50000 #define ZERO_MAX 100 /* * Streak status flags: streaks awaiting their starting time, * Streaks currently alive (eg. as of the ending time of the last * interval) and streaks that have either exited the vector field. */ #define STREAK_NEW 1 #define STREAK_ALIVE 2 #define STREAK_DEAD 3 #define STREAK_PENDING 4 typedef struct cacheStreak *CacheStreak; typedef struct cacheObject *CacheObject; struct cacheStreak { Series streak; /* If the streak is alive, this is it. */ double time; /* If not, this is the time to start it. */ POINT_TYPE point[3]; /* And this is the place to start it. */ VECTOR_TYPE norm[3]; /* And this is the propagated normal */ VECTOR_TYPE binorm[3]; /* And this is the propagated binormal */ VECTOR_TYPE tan[3]; /* And this is the propagated tangent */ int status; /* status: waiting, alive or dead */ int start; /* offset of newly added vertices */ }; struct cacheObject { int frame; int nDim; float t0, t1; Object v0, v1; Interpolator c0, c1; VectorField vf0, vf1; int nStreaks; CacheStreak streaks; float *frameTimes; }; typedef struct streakVars { VectorField vf0; VectorField vf1; InstanceVars I0; InstanceVars I1; float t0, t1, t; int inside0, inside1; int nDim; } StreakVars; static Array Starts(Object, Object, int); static Array Times(Object, Object); static Error StartsRecurse(Object, Array *); static Error StartsAddPoints(Array *, Array); static Error AlignStartPtsAndTimes(Array, Array); static Error Streaklines(CacheObject, float, int); static Error Streakline(CacheObject, int, float); static Error StreakTask(Pointer); static Error GetFrameData(Object, Object, CacheObject); static Error TraceFrame(Pointer); static Error InitFrame(float *, float *, float *); static Error UpdateFrame(float *, float *, float *, float *, float *, float *, float *, float *, float *); static void RotateAroundVector(Vector *, float, float, float *); static int ZeroVector(VECTOR_TYPE *, int nDim); static int IsRegular(Object); static Error GetElementType(Object, char **); static Error GetCacheObject(Object *, CacheObject *, Object *); static Object SetCacheObject(Object *, CacheObject); static void GetCacheKey(Object *, Object *, int *); static Error FreeCacheObject(Pointer); static Error NewCacheObject(Array, Array, CacheObject *); static Stream NewStreakBuf(int); static void DestroyStreakBuf(Stream); static Field StreakToField(Stream); static Error SetupStreakVars(StreakVars *, CacheObject); static Error FreeStreakVars(StreakVars *); static int Streak_FindElement(StreakVars *, POINT_TYPE *); static int Streak_FindMultiGridContinuation(StreakVars *, POINT_TYPE *); static Error Streak_Interpolate(StreakVars *, POINT_TYPE *, double, VECTOR_TYPE *); static Error Streak_StepTime(StreakVars *, double, VECTOR_TYPE *, double *); static int Streak_Weights(StreakVars *, POINT_TYPE *); static int Streak_FaceWeights(StreakVars *, POINT_TYPE *); static int Streak_InsideWeights(StreakVars *, POINT_TYPE *); static Error Streak_FindBoundary(StreakVars *, POINT_TYPE *, POINT_TYPE *, double *); static int Streak_Neighbor(StreakVars *, VECTOR_TYPE *); static StreakFlags AddPointToStreak(Stream, POINT_TYPE *, VECTOR_TYPE *, double); static Error GetTail(Object, char *, Pointer); static int IsSeries(Object o); static Error ResetVectorField(VectorField); static Error GeometryCheck(Object, int); static VectorField InitVectorField(Object); static Error DestroyVectorField(VectorField); Error m_Streakline(Object *in, Object *out) { Object vectors, curls = NULL; int curlFlag = 0; float c; int head; int frame; Type type; Category cat; int rank; int nD, n; CacheObject cstreaks; int thisframe; int i, f; Group group = NULL; CacheStreak streak; Array starts = NULL; Array times = NULL; Class vClass; Object cachedObject = NULL; float termTime, seriesEnd; out[0] = NULL; if (! in[DATA_ARG]) { DXSetError(ERROR_MISSING_DATA, "#10000", "data"); goto error; } vectors = in[DATA_ARG]; if ((vClass = DXGetObjectClass(vectors)) == CLASS_GROUP) vClass = DXGetGroupClass((Group)vectors); if (vClass != CLASS_FIELD && vClass != CLASS_COMPOSITEFIELD && vClass != CLASS_MULTIGRID && vClass != CLASS_SERIES) { DXSetError(ERROR_INVALID_DATA, "data must be field, multigrid, composite field or series"); goto error; } if (! DXGetType(vectors, &type, &cat, &rank, &nD)) { DXSetError(ERROR_MISSING_DATA, "#10240", "data"); goto error; } if (type != TYPE_FLOAT || cat != CATEGORY_REAL || (nD != 2 && nD != 3)) { DXSetError(ERROR_BAD_PARAMETER, "#10302", "data"); goto error; } if (! GeometryCheck(vectors, nD)) goto error; if (IsSeries(vectors)) { int i; float f; for (i = 0; DXGetSeriesMember((Series)vectors, i, &f); i++) seriesEnd = f; head = i - 1; } /* * If head is given, then OK. */ if (in[HEAD_ARG]) { if (! DXExtractInteger(in[HEAD_ARG], &head)) { DXSetError(ERROR_BAD_PARAMETER, "head"); goto error; } } /* * Regarding frame: if its not given, then we run to head time. */ if (in[FRAME_ARG]) { if (! DXExtractInteger(in[FRAME_ARG], &frame)) { DXSetError(ERROR_BAD_PARAMETER, "frame"); goto error; } } else frame = head; /* * If the input is a series, then the termination time should * be the series position of the frame'th member of the * series. If the input is NOT a series, then its found * in the series position attribute. */ if (IsSeries(vectors)) { if (! DXGetSeriesMember((Series)vectors, frame, &termTime)) termTime = seriesEnd; } else { Object attr = DXGetAttribute((Object)vectors, "series position"); if (! attr) { DXSetError(ERROR_BAD_PARAMETER, "vector series member must have series position attribute"); goto error; } if (! DXExtractFloat(attr, &termTime)) { DXSetError(ERROR_INVALID_DATA, "series position attribute"); goto error; } } if (in[CURL_ARG]) { int nDC; curls = in[CURL_ARG]; if (! DXGetType(curls, &type, &cat, &rank, &nDC)) goto error; if (IsSeries(vectors) && !IsSeries(curls)) { DXSetError(ERROR_BAD_PARAMETER, "since data is series curl must be series"); goto error; } if (type != TYPE_FLOAT || cat != CATEGORY_REAL || rank != 1) { DXSetError(ERROR_BAD_PARAMETER, "curls must be float/real vectors"); goto error; } if (nDC != nD) { DXSetError(ERROR_BAD_PARAMETER, "curl vectors must match data in dimensionality"); goto error; } curlFlag = 1; } if (in[FLAG_ARG]) { if (!DXExtractInteger(in[FLAG_ARG], &curlFlag) || (curlFlag < 0 || curlFlag > 1)) { DXSetError(ERROR_BAD_PARAMETER, "curl flag must be 0 or 1"); goto error; } } if (in[SCALE_ARG]) { if (! DXExtractFloat(in[SCALE_ARG], &c)) { DXSetError(ERROR_BAD_PARAMETER, "c"); goto error; } if (c <= 0.0) c = DEFAULT_C; } else c = DEFAULT_C; /* * Look in the cache to see if there is an object representing earlier * calls to streakline. */ if (! GetCacheObject(in, &cstreaks, &cachedObject)) goto error; if (cstreaks) { /* * if the current termination time is earlier than the * termination time of the streaks already computed, determine * the frame number of the termination. */ if (cstreaks->frame > 0) { if (termTime <= cstreaks->frameTimes[cstreaks->frame]) { for (frame = 0; frame <= cstreaks->frame; frame++) if (termTime == cstreaks->frameTimes[frame]) break; if (frame == cstreaks->frame+1) { DXSetError(ERROR_BAD_PARAMETER, "out-of-sequence series member"); goto error; } } } } else { /* * If there isn't a streak object, then we create a new one. * * Starting points are gathered from argument or centerpoint of * first member of vector series */ starts = Starts(in[POINTS_ARG], in[DATA_ARG], nD); if (starts) DXGetArrayInfo(starts, &n, NULL, NULL, NULL, NULL); if (! starts || n == 0) goto no_streaks; if (! DXTypeCheck(starts, TYPE_FLOAT, CATEGORY_REAL, 1, nD)) { DXSetError(ERROR_BAD_PARAMETER, "starts"); goto error; } /* * Starting points are gathered from argument, or are extracted * from the series position of the first member of the vector series. */ times = Times(in[TIMES_ARG], in[DATA_ARG]); if (! times) goto error; /* * Match up the starting points and times */ if (! AlignStartPtsAndTimes(starts, times)) goto error; if (! NewCacheObject(starts, times, &cstreaks)) goto error; DXDelete((Object)starts); starts = NULL; DXDelete((Object)times); times = NULL; } if (frame <= cstreaks->frame) goto past_head; thisframe = cstreaks->frame + 1; /* * If a frame argument was given, then make sure its in sequence. * If frame argument was NOT given, then if the input is a series, frame * is the same as head, which is either given or the number of series * members. If the input is NOT a series, then frame is the current * frame only. */ if (in[FRAME_ARG]) { if (!IsSeries(vectors) && frame != thisframe) { DXSetError(ERROR_BAD_PARAMETER, "frame parameter given is out of sequence"); goto error; } if (frame > head) frame = head; } else if (IsSeries(vectors)) frame = head; else frame = thisframe; /* * The following loop updates the frame and generates streaks. */ for (f = thisframe; f <= frame; f++) { if (! GetFrameData(vectors, curls, cstreaks)) goto error; if (! Streaklines(cstreaks, c, curlFlag)) goto error; } if (!cachedObject) if (NULL == (cachedObject = SetCacheObject(in, cstreaks))) goto error; past_head: streak = cstreaks->streaks; for (i = 0; i < cstreaks->nStreaks; i++, streak++) { Series s = NULL, o = NULL; Field f; int j, k; float t; streak->streak; k = 0; for (j = 0; j <= frame; j++) { f = (Field)DXGetSeriesMember(streak->streak, j, &t); if (!f) break; if (! DXEmptyField(f)) { if (s == NULL) { s = DXNewSeries(); if (! s) goto error; } if (! DXSetSeriesMember(s, k++, t, (Object)f)) { DXDelete((Object)s); goto error; } } } if (s) { if (group == NULL) { group = DXNewGroup(); if (! group) { DXDelete((Object)s); goto error; } } if (! DXSetMember(group, NULL, (Object)s)) { DXDelete((Object)s); goto error; } } } if (! group) out[0] = (Object)DXEndField(DXNewField()); else { DXSetFloatAttribute((Object)group, "fuzz", 4); out[0] = (Object)group; } DXDelete(cachedObject); return OK; no_streaks: DXDelete(cachedObject); DXDelete((Object)group); DXDelete((Object)starts); DXDelete((Object)times); out[0] = NULL; return OK; error: DXDelete(cachedObject); DXDelete((Object)starts); DXDelete((Object)times); DXDelete((Object)group); out[0] = NULL; return ERROR; } typedef struct streakargs { CacheObject cstreak; int i; float c; } StreakArgs; typedef struct frameargs { int curlFlag; Interpolator i0, i1; float t0, t1; CacheStreak streak; } FrameArgs; static Error Streaklines(CacheObject cstreak, float c, int curlFlag) { Object member = NULL; int i, j, n; StreakArgs stask; FrameArgs ftask; int nDim; Field f; int nP; Interpolator curlMap0 = NULL; Interpolator curlMap1 = NULL; MultiGrid cObject0 = NULL; MultiGrid cObject1 = NULL; /* * If this is the 0-th frame, put an empty field into the streaks * Otherwise, put any necessary streaks. */ if (cstreak->frame == 0) { CacheStreak streak = cstreak->streaks; Field eField = DXEndField(DXNewField()); for (i = 0; i < cstreak->nStreaks; i++, streak++) { if (! DXSetSeriesMember(streak->streak, 0, cstreak->t1, (Object)eField)) goto error; } } else { if (cstreak->t1 <- cstreak->t0) cstreak->t1 = cstreak->t0; nP = DXProcessors(0); n = cstreak->nStreaks; if (n < 1) return OK; if (n == 1 || nP == 1) { for (i = 0; i < n; i++) { if (! Streakline(cstreak, i, c)) goto error; } } else { stask.cstreak = cstreak; stask.c = c; if (! DXCreateTaskGroup()) goto error; for (i = 0; i < n; i++) { stask.i = i; if (! DXAddTask(StreakTask, (Pointer)&stask, sizeof(stask), 1.0)) { DXAbortTaskGroup(); goto error; } } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; } if (cstreak->nDim == 3) { if (curlFlag) { if (cstreak->c0) curlMap0 = cstreak->c0; else { int i; cObject0 = DXNewMultiGrid(); if (! cObject0) goto error; for (i = 0; i < cstreak->vf0->nmembers; i++) { VectorGrp vg = cstreak->vf0->members[i]; if (! (*(vg->CurlMap))(vg, cObject0)) goto error; } curlMap0 = DXNewInterpolator((Object)cObject0, INTERP_INIT_IMMEDIATE, -1.0); if (! curlMap0) goto error; } DXReference((Object)curlMap0); if (cstreak->c1) curlMap1 = cstreak->c1; else { int i; cObject1 = DXNewMultiGrid(); if (! cObject1) goto error; for (i = 0; i < cstreak->vf1->nmembers; i++) { VectorGrp vg = cstreak->vf1->members[i]; if (! (*(vg->CurlMap))(vg, cObject1)) goto error; } curlMap1 = DXNewInterpolator((Object)cObject1, INTERP_INIT_IMMEDIATE, 0.0); if (! curlMap1) goto error; } DXReference((Object)curlMap1); if (! curlMap0 || ! curlMap1) { DXSetError(ERROR_INTERNAL, "missing curl map"); goto error; } ftask.curlFlag = 1; ftask.i0 = curlMap0; ftask.i1 = curlMap1; } else { ftask.curlFlag = 0; ftask.i0 = NULL; ftask.i1 = NULL; } ftask.t0 = cstreak->t0; ftask.t1 = cstreak->t1; if (n == 1 || nP == 1) { for (i = 0; i < n; i++) { ftask.streak = cstreak->streaks + i; if (! TraceFrame((Pointer)&ftask)) goto error; } } else { if (! DXCreateTaskGroup()) goto error; for (i = 0; i < n; i++) { ftask.streak = cstreak->streaks + i; if (! DXAddTask(TraceFrame, (Pointer)&ftask, sizeof(ftask), 1.0)) { DXAbortTaskGroup(); goto error; } } if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE) goto error; } DXDelete((Object)curlMap0); DXDelete((Object)curlMap1); curlMap0 = curlMap1 = NULL; } } no_streaks: return OK; error: if (cObject0) DXDelete((Object)cObject0); if (cObject1) DXDelete((Object)cObject1); if (curlMap0) DXDelete((Object)curlMap0); if (curlMap1) DXDelete((Object)curlMap1); return ERROR; } static Error StreakTask(Pointer p) { StreakArgs *t = (StreakArgs *)p; return Streakline(t->cstreak, t->i, t->c); } static Error Streakline(CacheObject cstreak, /* accumulated streak info */ int which, /* which to operate on */ float c) /* step criterion */ { double elapsedtime = 0; float head = cstreak->t1; CacheStreak streak = cstreak->streaks + which; int nDim = cstreak->nDim; Stream streakBuf = NULL; POINT_TYPE pb0[3], pb1[3]; VECTOR_TYPE vb0[3], vb1[3]; POINT_TYPE *p0, *p1, *ptmp; VECTOR_TYPE *v0, *v1, *vtmp; double t, t2; int done, keepPoint; Field field = NULL; int i, totKnt, zero, zeroKnt; StreakFlags sFlag; StreakVars svars; int zeroVector; memset(&svars, 0, sizeof(svars)); /* * Double buffer pointers */ p0 = pb0; p1 = pb1; v0 = vb0; v1 = vb1; /* * If streak has already died, well, OK. */ if (streak->status == STREAK_DEAD) return OK; /* * If streak is awaiting birth and its start time is after the current * interval, well, OK. But got to put an empty field marker in. */ if (streak->status == STREAK_PENDING && streak->time > head) goto streak_done; /* * If streak is alive, then we need to get its current termination point. * If the streak is newborn, then the point is already in the streak * header. Otherwise, we need to get it from the fragment already * computed. */ if (streak->status == STREAK_ALIVE) { float fbuf[3]; GetTail((Object)streak->streak,"positions",(Pointer)&fbuf[0]); for (i = 0; i < nDim; i++) streak->point[i] = fbuf[i]; } else streak->status = STREAK_ALIVE; for (i = 0; i < nDim; i++) p0[i] = streak->point[i]; if (streak->time > cstreak->t0) elapsedtime = streak->time; else elapsedtime = cstreak->t0; if (! SetupStreakVars(&svars, cstreak)) goto error; nDim = cstreak->nDim; /* * If we can't find the start, we can't extend the streak. */ if (! Streak_FindElement(&svars, p0)) { streak->status = STREAK_DEAD; goto streak_done; } if (! Streak_Interpolate(&svars, p0, elapsedtime, v0)) goto error; streakBuf = NewStreakBuf(nDim); if (! streakBuf) goto error; sFlag = AddPointToStreak(streakBuf, p0, v0, elapsedtime); if (sFlag == STREAK_ERROR) goto error; /* * While streak is passing from partition to partition... */ done = 0; totKnt = 1; zeroKnt = 0; while (!done) { zeroVector = ZeroVector(v0, nDim); if (zeroVector) { /* * Interpolate this point at t1. If its zero there too, * then it'll be zero everywhere in between (due to linear * interpolation in time. */ if (! Streak_Interpolate(&svars, p0, svars.t1, v1)) goto error; if (ZeroVector(v1, nDim)) { sFlag = AddPointToStreak(streakBuf, p0, v0, svars.t1); if (sFlag == STREAK_ERROR) goto error; streak->status = STREAK_ALIVE; goto streak_done; } else { /* * If it isn't zero at t1, then take a small step * in time to get off zero. */ while (zeroVector) { elapsedtime += 0.025 * (svars.t1 - elapsedtime); if (! Streak_Interpolate(&svars, p0, elapsedtime, v0)) goto error; zeroVector = ZeroVector(v0, nDim); } } } if (! Streak_StepTime(&svars, c, v0, &t)) goto error; /* * If the time step exceeds the head, truncate it and mark * the streak complete. */ if (elapsedtime + t > svars.t1) { t = svars.t1 - elapsedtime; streak->status = STREAK_ALIVE; done = 1; } /* * Get the point a time step away in the direction of the * vector */ if (! zeroVector) for (i = 0; i < nDim; i++) { p1[i] = p0[i] + t*v0[i]; if (p1[i] == p0[i]) v0[i] = 0.0; } /* * If the new point is not in the current element, find the point * at which it exits the element and iterponlate the point there. * Then move on to the neighbor, if one exists. * * Otherwise, just interpolate the endpoint. */ if (! Streak_InsideWeights(&svars, p1)) { int found; /* * If the new sample isn't in the previous element, find the * exit point of the streak and interpolate in the current * element there. */ if (! Streak_FindBoundary(&svars, p0, v0, &t2)) goto error; elapsedtime += t2; for (i = 0; i < nDim; i++) { double d = ((double)p0[i]) + ((double)t2)*((double)v0[i]); p1[i] = d; } if (! Streak_FaceWeights(&svars, p1)) goto error; if (! Streak_Interpolate(&svars, p1, elapsedtime, v1)) goto error; found = Streak_Neighbor(&svars, v1); if (found == -1) goto error; if (! found) found = Streak_FindMultiGridContinuation(&svars, p0); if (! found) { streak->status = STREAK_DEAD; done = 1; } else done = 0; if (done || t2 > 0) keepPoint = 1; else keepPoint = 0; t = t2; } else { keepPoint = 1; elapsedtime += t; if (! Streak_Interpolate(&svars, p1, elapsedtime, v1)) goto error; } if (keepPoint && t > 0.0) { sFlag = AddPointToStreak(streakBuf, p1, v1, elapsedtime); if (sFlag == STREAK_ERROR) goto error; else if (sFlag == STREAK_FULL) done = 1; } if (++totKnt > TOT_MAX) { DXWarning("possible infinite loop detected"); done = 1; } ptmp = p0; p0 = p1; p1 = ptmp; vtmp = v0; v0 = v1; v1 = vtmp; } streak_done: if (streakBuf) { field = StreakToField(streakBuf); if (! field) goto error; DestroyStreakBuf(streakBuf); streakBuf = NULL; } else field = DXEndField(DXNewField()); if (! DXSetSeriesMember(streak->streak, cstreak->frame, svars.t0, (Object)field)) goto error; FreeStreakVars(&svars); return OK; error: FreeStreakVars(&svars); DestroyStreakBuf(streakBuf); DXDelete((Object)field); return ERROR; } static Error TraceFrame(Pointer ptr) { Interpolator i0 = ((FrameArgs *)ptr)->i0; Interpolator i1 = ((FrameArgs *)ptr)->i1; float t0 = ((FrameArgs *)ptr)->t0; float t1 = ((FrameArgs *)ptr)->t1; CacheStreak cs = ((FrameArgs *)ptr)->streak; int flag = ((FrameArgs *)ptr)->curlFlag; Series s = cs->streak; Array cArray0 = NULL; Array cArray1 = NULL; Field field = NULL; Object next; Field prev; Array nArray = NULL, bArray = NULL; Array vArray, tArray, pArray; float *points, *vectors, *curls0, *curls1, *normals, *binormals, *time; int i, nVectors, nDim, seg; Object dattr = NULL; float fn[3], fb[3], ft[3]; field = NULL; prev = NULL; for (seg = 0; NULL != (next = DXGetSeriesMember(s, seg, NULL)); seg++) { if (field && ! DXEmptyField(field)) prev = field; field = (Field)next; } if (DXEmptyField(field)) return OK; pArray = (Array)DXGetComponentValue(field, "positions"); if (! pArray) goto error; points = (float *)DXGetArrayData(pArray); if (flag) { i0 = (Interpolator)DXCopy((Object)i0, COPY_STRUCTURE); if (! i0) goto error; cArray0 = DXMapArray(pArray, i0, NULL); if (! cArray0) goto error; DXDelete((Object)i0); i1 = (Interpolator)DXCopy((Object)i1, COPY_STRUCTURE); if (! i1) goto error; cArray1 = DXMapArray(pArray, i1, NULL); if (! cArray1) goto error; DXDelete((Object)i1); curls0 = (float *)DXGetArrayData(cArray0); curls1 = (float *)DXGetArrayData(cArray1); if (! curls0 || ! curls1) goto error; } 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; tArray = (Array)DXGetComponentValue(field, "time"); if (! tArray) goto error; time = (float *)DXGetArrayData(tArray); if (! time) goto error; 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 this is not the first segment of the stream, we need to get * the propagated frame of reference and initial vertex normal and * binormal from the cached stuff. Otherwise, we initialize the * process. */ if (prev) { fn[0] = cs->norm[0]; fn[1] = cs->norm[1]; fn[2] = cs->norm[2]; fb[0] = cs->binorm[0]; fb[1] = cs->binorm[1]; fb[2] = cs->binorm[2]; ft[0] = cs->tan[0]; ft[1] = cs->tan[1]; ft[2] = cs->tan[2]; GetTail((Object)prev, "normals", (Pointer)normals); GetTail((Object)prev, "binormals", (Pointer)binormals); } else { if (! InitFrame(vectors, fn, fb)) goto error; 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); normals[0] = fn[0]; normals[1] = fn[1]; normals[2] = fn[2]; binormals[0] = fb[0]; binormals[1] = fb[1]; binormals[2] = fb[2]; } for (i = 1; i < nVectors; i++) { points += 3; time += 1; vectors += 3; normals += 3; binormals += 3; if (flag) { float avcurl[3], pcurl[3], ncurl[3], d; d = (*time - t0) / (t1 - t0); pcurl[0] = curls0[0] + d*(curls1[0] - curls0[0]); pcurl[1] = curls0[1] + d*(curls1[1] - curls0[1]); pcurl[2] = curls0[2] + d*(curls1[2] - curls0[2]); ncurl[0] = curls0[3] + d*(curls1[3] - curls0[3]); ncurl[1] = curls0[4] + d*(curls1[4] - curls0[4]); ncurl[2] = curls0[5] + d*(curls1[5] - curls0[5]); avcurl[0] = 0.5 * (pcurl[0] + ncurl[0]); avcurl[1] = 0.5 * (pcurl[1] + ncurl[1]); avcurl[2] = 0.5 * (pcurl[2] + ncurl[2]); curls0 += 3; curls1 += 3; if (! UpdateFrame(points, vectors, avcurl, time, normals, binormals, fn, fb, ft)) goto error; } else { if (! UpdateFrame(points, vectors, NULL, time, normals, binormals, fn, fb, ft)) goto error; } } cs->norm[0] = fn[0]; cs->norm[1] = fn[1]; cs->norm[2] = fn[2]; cs->binorm[0] = fb[0]; cs->binorm[1] = fb[1]; cs->binorm[2] = fb[2]; cs->tan[0] = ft[0]; cs->tan[1] = ft[1]; cs->tan[2] = ft[2]; dattr = (Object)DXNewString("positions"); if (! dattr) goto error; DXReference(dattr); if (! DXSetComponentValue(field, "normals", (Object)nArray)) goto error; nArray = NULL; if (! DXSetComponentAttribute(field, "normals", "dep", dattr)) goto error; if (! DXSetComponentValue(field, "binormals", (Object)bArray)) goto error; bArray = NULL; if (! DXSetComponentAttribute(field, "binormals", "dep", dattr)) goto error; DXDelete(dattr); DXDelete((Object)cArray0); DXDelete((Object)cArray1); done: return OK; error: DXDelete((Object)nArray); DXDelete((Object)bArray); DXDelete((Object)cArray0); DXDelete((Object)cArray1); return ERROR; } static Array Starts(Object starts, Object vectors, int nD) { Array array = NULL; if (starts) { Object 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 Array Times(Object times, Object vectors) { Array array = NULL; if (times) { if (! StartsRecurse(times, &array)) goto error; } else if (vectors) { float time; if (IsSeries(vectors)) DXGetSeriesMember((Series)vectors, 0, &time); else { Object attr = DXGetAttribute(vectors, "series position"); if (attr == NULL) DXSetError(ERROR_MISSING_DATA, "series position attribute on data"); if (! DXExtractFloat(attr, &time)) DXSetError(ERROR_INVALID_DATA, "series position attribute on data"); } array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0); if (! array) goto error; if (! DXAddArrayData(array, 0, 1, (Pointer)&time)) goto error; } else { float *ptr; int i; array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1); if (! array) goto error; if (! DXAddArrayData(array, 0, 1, NULL)) goto error; ptr = (float *)DXGetArrayData(array); if (! ptr) goto error; *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 OK; 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 NewStreakBuf(int nDim) { Stream s = NULL; s = (Stream)DXAllocateZero(sizeof(struct stream)); if (! s) goto error; s->pArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nDim); s->tArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0); s->vArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nDim); if (s->pArray == NULL || s->vArray == NULL || s->tArray == NULL) goto error; s->arrayKnt = 0; s->points = (float *)DXAllocate(nDim*STREAM_BUF_POINTS*sizeof(float)); s->time = (float *)DXAllocate(STREAM_BUF_POINTS*sizeof(float)); s->vectors = (float *)DXAllocate(nDim*STREAM_BUF_POINTS*sizeof(float)); if (s->points == NULL || s->time == NULL || s->vectors == NULL) goto error; s->nDim = nDim; s->bufKnt = 0; return s; error: DestroyStreakBuf(s); return NULL; } static void DestroyStreakBuf(Stream s) { if (s) { DXDelete((Object)s->pArray); DXDelete((Object)s->vArray); DXDelete((Object)s->tArray); DXFree((Pointer)s->points); DXFree((Pointer)s->time); DXFree((Pointer)s->vectors); DXFree((Pointer)s); } } static Error FlushStreak(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 StreakFlags AddPointToStreak(Stream s, POINT_TYPE *p, VECTOR_TYPE *v, double t) { if (s->bufKnt == STREAM_BUF_POINTS) if (! FlushStreak(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->bufKnt + s->arrayKnt >= STREAK_MAX) return STREAK_FULL; return STREAK_OK; error: return STREAK_ERROR; } static Field StreakToField(Stream buf) { float c[3]; float d[3]; Array a = NULL; Object dattr = NULL; Field field; int flag = 1; if (! FlushStreak(buf)) goto error; c[0] = 0.7; c[1] = 0.7; c[2] = 0.0; d[0] = 0.0; d[1] = 0.0; d[2] = 0.0; dattr = (Object)DXNewString("positions"); if (! dattr) goto error; DXReference(dattr); field = DXNewField(); if (! field) goto error; if (! DXSetComponentValue(field, "positions", (Object)buf->pArray)) goto error; buf->pArray = NULL; if (! DXSetComponentAttribute(field, "positions", "dep", dattr)) goto error; if (! DXSetComponentValue(field, "data", (Object)buf->vArray)) goto error; buf->vArray = NULL; if (! DXSetComponentAttribute(field, "data", "dep", dattr)) goto error; if (! DXSetComponentValue(field, "time", (Object)buf->tArray)) goto error; buf->tArray = NULL; if (! DXSetComponentAttribute(field, "time", "dep", dattr)) goto error; DXDelete(dattr); a = DXMakeGridConnections(1, buf->arrayKnt); if (! a) goto error; if (! DXSetComponentValue(field, "connections", (Object)a)) goto error; a = NULL; if (! DXSetComponentAttribute(field, "connections", "element type", (Object)DXNewString("lines"))) goto error; a = (Array)DXNewRegularArray(TYPE_FLOAT, 3, buf->arrayKnt, (Pointer)c, (Pointer)d); if (! a) goto error; if (! DXSetComponentValue(field, "colors", (Object)a)) goto error; a = NULL; if (! DXSetIntegerAttribute((Object)field, "shade", 0)) goto error; if (! DXEndField(field)) goto error; return field; error: DXDelete(dattr); return ERROR; } 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 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 *)v, (Vector *)n, (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 *t, 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], vt[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); 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 NewCacheObject(Array pArray, Array tArray, CacheObject *co) { CacheStreak cs = NULL; float *points, *times; int nDim, nStreaks, i, j; *co = NULL; points = (float *)DXGetArrayData(pArray); if (! points) goto error; times = (float *)DXGetArrayData(tArray); if (! times) goto error; if (! DXGetArrayInfo(pArray, &nStreaks, NULL, NULL, NULL, &nDim)) goto error; *co = (CacheObject)DXAllocateZero(sizeof(struct cacheObject)); if (! *co) goto error; cs = (CacheStreak)DXAllocateZero(nStreaks*sizeof(struct cacheStreak)); if (!cs) goto error; (*co)->streaks = cs; (*co)->nDim = nDim; for (i = 0; i < nStreaks; i++, cs++) { for (j = 0; j < nDim; j++) cs->point[j] = *points++; cs->time = *times++; cs->status = STREAK_PENDING; cs->start = 0; cs->streak = DXNewSeries(); if (cs->streak == NULL) goto error; } /* * Note: no frame data is currently contained in the following * structure */ (*co)->frame = -1; (*co)->nStreaks = nStreaks; (*co)->frameTimes = NULL; return OK; error: FreeCacheObject(*co); return ERROR; } static Error FreeCacheObject(Pointer ptr) { int i; CacheObject co = ptr; if (co) { if (co->streaks) { for (i = 0; i < co->nStreaks; i++) if (co->streaks[i].streak) DXDelete((Object)co->streaks[i].streak); DXFree((Pointer)co->streaks); } if (co->frameTimes) DXFree(co->frameTimes); if (co->v0) DXDelete(co->v0); if (co->c0) DXDelete((Object)co->c0); if (co->vf0) DestroyVectorField(co->vf0); if (co->v1) DXDelete(co->v1); if (co->c1) DXDelete((Object)co->c1); if (co->vf1) DestroyVectorField(co->vf1); DXFree((Pointer)co); } return OK; } static void GetCacheKey(Object *in, Object *keys, int *nk) { int n = 0; if (in[DATA_ARG] && IsSeries(in[DATA_ARG])) keys[n++] = in[DATA_ARG]; if (in[NAME_ARG]) keys[n++] = in[NAME_ARG]; if (in[POINTS_ARG]) keys[n++] = in[POINTS_ARG]; if (in[TIMES_ARG]) keys[n++] = in[TIMES_ARG]; if (in[HEAD_ARG]) keys[n++] = in[HEAD_ARG]; if (in[FLAG_ARG]) keys[n++] = in[FLAG_ARG]; if (in[SCALE_ARG]) keys[n++] = in[SCALE_ARG]; *nk = n; } static Error GetCacheObject(Object *in, CacheObject *co, Object *obj) { int nk; Object keys[8]; Private pco; int hold = 0; if (in[HOLD_ARG]) if (!DXExtractInteger(in[HOLD_ARG], &hold) || hold < 0 || hold > 1) { DXSetError(ERROR_BAD_PARAMETER, "hold must be integer flag"); return ERROR; } if (hold) { int frame; if (! in[FRAME_ARG]) { DXSetError(ERROR_BAD_PARAMETER, "if hold is set, frame is required"); return ERROR; } if (!DXExtractInteger(in[FRAME_ARG], &frame)) { DXSetError(ERROR_BAD_PARAMETER, "frame must be integer"); return ERROR; } if (frame != 0) { pco = (Private)DXGetCacheEntry("Streakline", 0, 1, in[NAME_ARG]); if (! pco) { DXSetError(ERROR_BAD_PARAMETER, "hold set and frame not 0, but no cached state is present"); return ERROR; } } else { DXSetCacheEntry(NULL, CACHE_PERMANENT, "Streakline", 0, 1, in[NAME_ARG]); pco = NULL; } } else { GetCacheKey(in, keys, &nk); pco = (Private)DXGetCacheEntryV("Streakline", 0, nk, keys); } if (! pco) { *co = NULL; *obj = NULL; } else { *co = (CacheObject)DXGetPrivateData(pco); *obj = DXReference((Object)pco); } return OK; } static Object SetCacheObject(Object *in, CacheObject co) { Private pco = NULL; Object keys[8]; int nk; pco = DXNewPrivate((Pointer)co, FreeCacheObject); if (! pco) return ERROR; DXReference((Object)pco); if (in[HOLD_ARG]) { if (! DXSetCacheEntry((Object)pco, CACHE_PERMANENT, "Streakline", 0, 1, in[NAME_ARG])) return ERROR; } else { GetCacheKey(in, keys, &nk); if (! DXSetCacheEntryV((Object)pco, CACHE_PERMANENT, "Streakline", 0, nk, keys)) goto error; } return (Object)pco; error: DXDelete((Object)pco); return NULL; } static Error GetTail(Object o, char *name, Pointer l) { Field f; Object next; int n, nd, i; byte *ptr; int size; Array a; if (DXGetObjectClass(o) == CLASS_FIELD) { f = (Field)o; } else /* object is a series */ { f = NULL; for (i = 0; NULL != (next = DXGetSeriesMember((Series)o, i, NULL)); i++) f = (Field)next; } if (! f || DXEmptyField(f)) { DXSetError(ERROR_INTERNAL, "GetTail called with bad field"); return ERROR; } a = (Array)DXGetComponentValue(f, name); if (! a) { DXSetError(ERROR_INTERNAL, "GetTail: no component %s in field", name); return ERROR; } size = DXGetItemSize(a); DXGetArrayInfo(a, &n, NULL, NULL, NULL, NULL); ptr = ((byte *)DXGetArrayData(a)) + (n-1)*size; memcpy(l, ptr, size); } static Error SetupStreakVars(StreakVars *svars, CacheObject cstreak) { svars->I0 = svars->I1 = NULL; svars->vf0 = cstreak->vf0; svars->vf1 = cstreak->vf1; svars->t0 = cstreak->t0; svars->t1 = cstreak->t1; svars->inside0 = svars->inside1 = 0; svars->nDim = cstreak->nDim; return OK; error: FreeStreakVars(svars); return ERROR; } static Error FreeStreakVars(StreakVars *svars) { if (svars->I0) (*(svars->I0->currentVectorGrp->FreeInstanceVars))(svars->I0); if (svars->I1) (*(svars->I1->currentVectorGrp->FreeInstanceVars))(svars->I1); svars->I0 = svars->I1 = NULL; return OK; } static Error GetFrameData(Object vectors, Object curls, CacheObject co) { Object curl = NULL; co->frame ++; if (co->v0) DXDelete((Object)co->v0); co->v0 = co->v1; if (co->c0) DXDelete((Object)co->c0); co->c0 = co->c1; if (co->vf0) { DestroyVectorField(co->vf0); co->vf0 = NULL; } if (co->vf1) { co->vf0 = co->vf1; if (! ResetVectorField(co->vf0)) goto error; } co->t0 = co->t1; if (IsSeries(vectors)) { co->v1 = DXGetSeriesMember((Series)vectors, co->frame, &co->t1); if (! co->v1) { DXSetError(ERROR_MISSING_DATA, "data: empty series group"); goto error; } co->v1 = DXCopy(co->v1, COPY_STRUCTURE); if (! co->v1) goto error; if (curls) { Object c; float t; curl = DXGetSeriesMember((Series)curls, co->frame, &t); if (! curl) { DXSetError(ERROR_MISSING_DATA, "curl: empty series group"); goto error; } if (t != co->t1) { DXSetError(ERROR_INVALID_DATA, "curl and data series positions do not match"); goto error; } } } else { Object attr; attr = DXGetAttribute(vectors, "series position"); if (! attr) { DXSetError(ERROR_MISSING_DATA, "series position"); goto error; } if (! DXExtractFloat(attr, &co->t1)) { DXSetError(ERROR_INVALID_DATA, "series position"); goto error; } co->v1 = DXCopy(vectors, COPY_STRUCTURE); if (! co->v1) goto error; curl = curls; } if (curl) { co->c1 = DXNewInterpolator(curl, INTERP_INIT_PARALLEL, 0.0); if (! co->c1) goto error; DXReference((Object)co->c1); } co->vf1 = InitVectorField(co->v1); if (! co->vf1) goto error; if (co->frameTimes) co->frameTimes = (float *)DXReAllocate(co->frameTimes, (co->frame+1)*sizeof(float)); else co->frameTimes = (float *)DXAllocate((co->frame+1)*sizeof(float)); if (! co->frameTimes) goto error; co->frameTimes[co->frame] = co->t1; return OK; error: return ERROR; } static InstanceVars Streak_FindElement_1(VectorField vf, POINT_TYPE *point, VectorGrp last) { InstanceVars I; int i; for (i = 0; i < vf->nmembers; i++) if (vf->members[i] != last) { I = (*(vf->members[i]->NewInstanceVars))(vf->members[i]); if ((*(I->currentVectorGrp->FindElement))(I, point)) { vf->current = vf->members[i]; return I; } (*(I->currentVectorGrp->FreeInstanceVars))(I); } return 0; } static InstanceVars Streak_FindMultiGridContinuation_1(VectorField vf, POINT_TYPE *point, VectorGrp last) { InstanceVars I; int i; for (i = 0; i < vf->nmembers; i++) if (vf->members[i] != last) { I = (*(vf->members[i]->NewInstanceVars))(vf->members[i]); if ((*(I->currentVectorGrp->FindMultiGridContinuation))(I, point)) { vf->current = vf->members[i]; return I; } (*(I->currentVectorGrp->FreeInstanceVars))(I); } return 0; } static int Streak_FindElement(StreakVars *svars, POINT_TYPE *p) { if (! svars->inside0) { VectorGrp last; if (svars->I0) { last = svars->I0->currentVectorGrp; (*(last->FreeInstanceVars))(svars->I0); } else last = NULL; svars->I0 = Streak_FindElement_1(svars->vf0, p, last); if (! svars->I0) return 0; } if (! svars->inside1) { VectorGrp last; if (svars->I1) { last = svars->I1->currentVectorGrp; (*(last->FreeInstanceVars))(svars->I1); } else last = NULL; svars->I1 = Streak_FindElement_1(svars->vf1, p, last); if (! svars->I1) return 0; } return 1; } static int Streak_FindMultiGridContinuation(StreakVars *svars, POINT_TYPE *p) { if (! svars->inside0) { VectorGrp last; if (svars->I0) { last = svars->I0->currentVectorGrp; (*(last->FreeInstanceVars))(svars->I0); } else last = NULL; svars->I0 = Streak_FindMultiGridContinuation_1(svars->vf0, p, last); if (! svars->I0) return 0; } if (! svars->inside1) { VectorGrp last; if (svars->I1) { last = svars->I1->currentVectorGrp; (*(last->FreeInstanceVars))(svars->I1); } else last = NULL; svars->I1 = Streak_FindMultiGridContinuation_1(svars->vf1, p, last); if (! svars->I1) return 0; } return 1; } static Error Streak_Interpolate(StreakVars *svars, POINT_TYPE *p, double t, VECTOR_TYPE *v) { double d; int i; VECTOR_TYPE v0[3], v1[3]; if (! (*(svars->I0->currentVectorGrp->Interpolate))(svars->I0, p, v0)) return ERROR; if (! (*(svars->I1->currentVectorGrp->Interpolate))(svars->I1, p, v1)) return ERROR; d = (t - svars->t0) / (svars->t1 - svars->t0); for (i = 0; i < svars->nDim; i++) v[i] = v0[i] + d*(v1[i] - v0[i]); return OK; } static Error Streak_StepTime(StreakVars *svars, double c, VECTOR_TYPE *v, double *t) { double t0, t1; if (! (*(svars->I0->currentVectorGrp->StepTime))(svars->I0, c, v, &t0)) return ERROR; if (! (*(svars->I1->currentVectorGrp->StepTime))(svars->I1, c, v, &t1)) return ERROR; if (t0 < t1) *t = t0; else *t = t1; return OK; } static int Streak_InsideWeights(StreakVars *svars, POINT_TYPE *p) { svars->inside0 = (*(svars->I0->currentVectorGrp->Weights))(svars->I0, p); svars->inside1 = (*(svars->I1->currentVectorGrp->Weights))(svars->I1, p); if (svars->inside0 && svars->inside1) return 1; else return 0; } static int Streak_Weights(StreakVars *svars, POINT_TYPE *p) { (*(svars->I0->currentVectorGrp->Weights))(svars->I0, p); (*(svars->I1->currentVectorGrp->Weights))(svars->I1, p); return 1; } static int Streak_FaceWeights(StreakVars *svars, POINT_TYPE *p) { if (! svars->inside0) (*(svars->I0->currentVectorGrp->FaceWeights))(svars->I0, p); if (! svars->inside1) (*(svars->I0->currentVectorGrp->FaceWeights))(svars->I1, p); return 1; } static Error Streak_FindBoundary(StreakVars *svars, POINT_TYPE *p, VECTOR_TYPE *v, double *t) { double t0, t1; if (svars->inside0) t0 = DXD_MAX_FLOAT; else if (! (*(svars->I0->currentVectorGrp->FindBoundary)) (svars->I0, p, v, &t0)) return ERROR; if (svars->inside1) t1 = DXD_MAX_FLOAT; else if (! (*(svars->I1->currentVectorGrp->FindBoundary)) (svars->I1, p, v, &t1)) return ERROR; if (t1 < t0) *t = t1; else *t = t0; if (*t == DXD_MAX_FLOAT) { DXSetError(ERROR_INTERNAL, "error in Streak_FindBoundary"); return ERROR; } return OK; } static int Streak_Neighbor(StreakVars *svars, VECTOR_TYPE *v) { int found; if (! svars->inside0) { found = (*(svars->I0->currentVectorGrp->Neighbor))(svars->I0, v); if (found != 1) return found; } if (! svars->inside1) { found = (*(svars->I1->currentVectorGrp->Neighbor))(svars->I1, v); if (found != 1) return found; } return 1; } static int IsSeries(Object o) { if (DXGetObjectClass(o) != CLASS_GROUP) return 0; if (DXGetGroupClass((Group)o) != CLASS_SERIES) return 0; return 1; } static Error ResetVectorField(VectorField vf) { int i; for (i = 0; i < vf->nmembers; i++) { VectorGrp vg = vf->members[i]; if (vg->Reset) if (! (*(vg->Reset))(vg)) return ERROR; } return OK; } 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_SERIES: 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; }