1 /***********************************************************************/
2 /* Open Visualization Data Explorer */
3 /* (C) Copyright IBM Corp. 1989,1999 */
4 /* ALL RIGHTS RESERVED */
5 /* This code licensed under the */
6 /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */
7 /***********************************************************************/
8 /*
9 * $Header: /src/master/dx/src/exec/dxmods/streakline.c,v 1.9 2006/06/10 16:33:58 davidt Exp $:
10 */
11
12 #include <dxconfig.h>
13
14 #if defined(HAVE_STRING_H)
15 #include <string.h>
16 #endif
17
18 #include <dx/dx.h>
19 #include <math.h>
20
21 #include "stream.h"
22 #include "vectors.h"
23
24 #define NAME_ARG 0
25 #define DATA_ARG 1
26 #define POINTS_ARG 2
27 #define TIMES_ARG 3
28 #define HEAD_ARG 4
29 #define FRAME_ARG 5
30 #define CURL_ARG 6
31 #define FLAG_ARG 7
32 #define SCALE_ARG 8
33 #define HOLD_ARG 9
34
35 typedef enum
36 {
37 STREAK_OK,
38 STREAK_FULL,
39 STREAK_ERROR
40 } StreakFlags;
41
42 #ifdef STREAM_MAX
43 #undef STREAM_MAX
44 #endif
45
46 #define STREAK_MAX 25000
47 #define TOT_MAX 50000
48 #define ZERO_MAX 100
49
50 /*
51 * Streak status flags: streaks awaiting their starting time,
52 * Streaks currently alive (eg. as of the ending time of the last
53 * interval) and streaks that have either exited the vector field.
54 */
55 #define STREAK_NEW 1
56 #define STREAK_ALIVE 2
57 #define STREAK_DEAD 3
58 #define STREAK_PENDING 4
59
60 typedef struct cacheStreak *CacheStreak;
61 typedef struct cacheObject *CacheObject;
62
63 struct cacheStreak
64 {
65 Series streak; /* If the streak is alive, this is it. */
66 double time; /* If not, this is the time to start it. */
67 POINT_TYPE point[3]; /* And this is the place to start it. */
68 VECTOR_TYPE norm[3]; /* And this is the propagated normal */
69 VECTOR_TYPE binorm[3]; /* And this is the propagated binormal */
70 VECTOR_TYPE tan[3]; /* And this is the propagated tangent */
71 int status; /* status: waiting, alive or dead */
72 int start; /* offset of newly added vertices */
73 };
74
75 struct cacheObject
76 {
77 int frame;
78 int nDim;
79 float t0, t1;
80 Object v0, v1;
81 Interpolator c0, c1;
82 VectorField vf0, vf1;
83 int nStreaks;
84 CacheStreak streaks;
85 float *frameTimes;
86 };
87
88 typedef struct streakVars
89 {
90 VectorField vf0;
91 VectorField vf1;
92 InstanceVars I0;
93 InstanceVars I1;
94 float t0, t1, t;
95 int inside0, inside1;
96 int nDim;
97 } StreakVars;
98
99 static Array Starts(Object, Object, int);
100 static Array Times(Object, Object);
101 static Error StartsRecurse(Object, Array *);
102 static Error StartsAddPoints(Array *, Array);
103 static Error AlignStartPtsAndTimes(Array, Array);
104 static Error Streaklines(CacheObject, float, int);
105 static Error Streakline(CacheObject, int, float);
106 static Error StreakTask(Pointer);
107 static Error GetFrameData(Object, Object, CacheObject);
108 static Error TraceFrame(Pointer);
109 static Error InitFrame(Vector *, Vector *, Vector *);
110 static Error UpdateFrame(float *, float *, float *,
111 float *, float *, float *,
112 Vector *, Vector *, Vector *);
113 static void RotateAroundVector(Vector *, float, float, float *);
114 static int ZeroVector(VECTOR_TYPE *, int nDim);
115 static int IsRegular(Object);
116 static Error GetElementType(Object, char **);
117 static Error GetCacheObject(Object *, CacheObject *, Object *);
118 static Object SetCacheObject(Object *, CacheObject);
119 static void GetCacheKey(Object *, Object *, int *);
120 static Error FreeCacheObject(Pointer);
121 static Error NewCacheObject(Array, Array, CacheObject *);
122 static Stream NewStreakBuf(int);
123 static void DestroyStreakBuf(Stream);
124 static Field StreakToField(Stream);
125 static Error SetupStreakVars(StreakVars *, CacheObject);
126 static Error FreeStreakVars(StreakVars *);
127 static int Streak_FindElement(StreakVars *, POINT_TYPE *);
128 static int Streak_FindMultiGridContinuation(StreakVars *, POINT_TYPE *);
129 static Error Streak_Interpolate(StreakVars *, POINT_TYPE *,
130 double, VECTOR_TYPE *);
131 static Error Streak_StepTime(StreakVars *, double, VECTOR_TYPE *, double *);
132 static int Streak_FaceWeights(StreakVars *, POINT_TYPE *);
133 static int Streak_InsideWeights(StreakVars *, POINT_TYPE *);
134 static Error Streak_FindBoundary(StreakVars *, POINT_TYPE *,
135 POINT_TYPE *, double *);
136 static int Streak_Neighbor(StreakVars *, VECTOR_TYPE *);
137 static StreakFlags AddPointToStreak(Stream, POINT_TYPE *, VECTOR_TYPE *, double);
138 static Error GetTail(Object, char *, Pointer);
139 static int IsSeries(Object o);
140 static Error ResetVectorField(VectorField);
141 static Error GeometryCheck(Object, int);
142
143 static VectorField InitVectorField(Object);
144 static Error DestroyVectorField(VectorField);
145
146 Error
m_Streakline(Object * in,Object * out)147 m_Streakline(Object *in, Object *out)
148 {
149 Object vectors, curls = NULL;
150 int curlFlag = 0;
151 float c;
152 int head;
153 int frame;
154 Type type;
155 Category cat;
156 int rank;
157 int nD, n;
158 CacheObject cstreaks;
159 int thisframe;
160 int i, f;
161 Group group = NULL;
162 CacheStreak streak;
163 Array starts = NULL;
164 Array times = NULL;
165 Class vClass;
166 Object cachedObject = NULL;
167 float termTime, seriesEnd=0;
168
169 out[0] = NULL;
170
171 if (! in[DATA_ARG])
172 {
173 DXSetError(ERROR_MISSING_DATA, "#10000", "data");
174 goto error;
175 }
176
177 vectors = in[DATA_ARG];
178
179 if ((vClass = DXGetObjectClass(vectors)) == CLASS_GROUP)
180 vClass = DXGetGroupClass((Group)vectors);
181
182 if (vClass != CLASS_FIELD &&
183 vClass != CLASS_COMPOSITEFIELD &&
184 vClass != CLASS_MULTIGRID &&
185 vClass != CLASS_SERIES)
186 {
187 DXSetError(ERROR_DATA_INVALID,
188 "data must be field, multigrid, composite field or series");
189 goto error;
190 }
191
192 if (! DXGetType(vectors, &type, &cat, &rank, &nD))
193 {
194 DXSetError(ERROR_MISSING_DATA, "#10240", "data");
195 goto error;
196 }
197
198 if (type != TYPE_FLOAT || cat != CATEGORY_REAL || (nD != 2 && nD != 3))
199 {
200 DXSetError(ERROR_BAD_PARAMETER, "#10302", "data");
201 goto error;
202 }
203
204 if (! GeometryCheck(vectors, nD))
205 goto error;
206
207 if (IsSeries(vectors))
208 {
209 int i;
210 float f;
211
212 for (i = 0; DXGetSeriesMember((Series)vectors, i, &f); i++)
213 seriesEnd = f;
214
215 head = i - 1;
216 }
217
218 /*
219 * If head is given, then OK.
220 */
221 if (in[HEAD_ARG])
222 {
223 if (! DXExtractInteger(in[HEAD_ARG], &head))
224 {
225 DXSetError(ERROR_BAD_PARAMETER, "head");
226 goto error;
227 }
228 }
229
230 /*
231 * Regarding frame: if its not given, then we run to head time.
232 */
233 if (in[FRAME_ARG])
234 {
235 if (! DXExtractInteger(in[FRAME_ARG], &frame))
236 {
237 DXSetError(ERROR_BAD_PARAMETER, "frame");
238 goto error;
239 }
240 }
241 else
242 frame = head;
243
244 /*
245 * If the input is a series, then the termination time should
246 * be the series position of the frame'th member of the
247 * series. If the input is NOT a series, then its found
248 * in the series position attribute.
249 */
250 if (IsSeries(vectors))
251 {
252 if (! DXGetSeriesMember((Series)vectors, frame, &termTime))
253 termTime = seriesEnd;
254 }
255 else
256 {
257 Object attr = DXGetAttribute((Object)vectors, "series position");
258 if (! attr)
259 {
260 DXSetError(ERROR_BAD_PARAMETER,
261 "vector series member must have series position attribute");
262 goto error;
263 }
264
265 if (! DXExtractFloat(attr, &termTime))
266 {
267 DXSetError(ERROR_DATA_INVALID, "series position attribute");
268 goto error;
269 }
270 }
271
272 if (in[CURL_ARG])
273 {
274 int nDC;
275
276 curls = in[CURL_ARG];
277 if (! DXGetType(curls, &type, &cat, &rank, &nDC))
278 goto error;
279
280 if (IsSeries(vectors) && !IsSeries(curls))
281 {
282 DXSetError(ERROR_BAD_PARAMETER,
283 "since data is series curl must be series");
284 goto error;
285 }
286
287 if (type != TYPE_FLOAT || cat != CATEGORY_REAL || rank != 1)
288 {
289 DXSetError(ERROR_BAD_PARAMETER, "curls must be float/real vectors");
290 goto error;
291 }
292
293 if (nDC != nD)
294 {
295 DXSetError(ERROR_BAD_PARAMETER,
296 "curl vectors must match data in dimensionality");
297 goto error;
298 }
299
300 curlFlag = 1;
301 }
302
303 if (in[FLAG_ARG])
304 {
305 if (!DXExtractInteger(in[FLAG_ARG], &curlFlag) ||
306 (curlFlag < 0 || curlFlag > 1))
307 {
308 DXSetError(ERROR_BAD_PARAMETER, "curl flag must be 0 or 1");
309 goto error;
310 }
311 }
312
313 if (in[SCALE_ARG])
314 {
315 if (! DXExtractFloat(in[SCALE_ARG], &c))
316 {
317 DXSetError(ERROR_BAD_PARAMETER, "c");
318 goto error;
319 }
320
321 if (c <= 0.0)
322 c = DEFAULT_C;
323 }
324 else
325 c = DEFAULT_C;
326
327 /*
328 * Look in the cache to see if there is an object representing earlier
329 * calls to streakline.
330 */
331 if (! GetCacheObject(in, &cstreaks, &cachedObject))
332 goto error;
333
334 if (cstreaks)
335 {
336 /*
337 * if the current termination time is earlier than the
338 * termination time of the streaks already computed, determine
339 * the frame number of the termination.
340 */
341 if (cstreaks->frame > 0)
342 {
343 if (termTime <= cstreaks->frameTimes[cstreaks->frame])
344 {
345 for (frame = 0; frame <= cstreaks->frame; frame++)
346 if (termTime == cstreaks->frameTimes[frame])
347 break;
348
349 if (frame == cstreaks->frame+1)
350 {
351 DXSetError(ERROR_BAD_PARAMETER,
352 "out-of-sequence series member");
353 goto error;
354 }
355 }
356 }
357 }
358 else
359 {
360 /*
361 * If there isn't a streak object, then we create a new one.
362 *
363 * Starting points are gathered from argument or centerpoint of
364 * first member of vector series
365 */
366 starts = Starts(in[POINTS_ARG], in[DATA_ARG], nD);
367 if (starts)
368 DXGetArrayInfo(starts, &n, NULL, NULL, NULL, NULL);
369 if (! starts || n == 0)
370 goto no_streaks;
371
372 if (! DXTypeCheck(starts, TYPE_FLOAT, CATEGORY_REAL, 1, nD))
373 {
374 DXSetError(ERROR_BAD_PARAMETER, "starts");
375 goto error;
376 }
377
378 /*
379 * Starting points are gathered from argument, or are extracted
380 * from the series position of the first member of the vector series.
381 */
382 times = Times(in[TIMES_ARG], in[DATA_ARG]);
383 if (! times)
384 goto error;
385
386 /*
387 * Match up the starting points and times
388 */
389 if (! AlignStartPtsAndTimes(starts, times))
390 goto error;
391
392 if (! NewCacheObject(starts, times, &cstreaks))
393 goto error;
394
395 DXDelete((Object)starts); starts = NULL;
396 DXDelete((Object)times); times = NULL;
397 }
398
399 if (frame <= cstreaks->frame)
400 goto past_head;
401
402 thisframe = cstreaks->frame + 1;
403
404 /*
405 * If a frame argument was given, then make sure its in sequence.
406 * If frame argument was NOT given, then if the input is a series, frame
407 * is the same as head, which is either given or the number of series
408 * members. If the input is NOT a series, then frame is the current
409 * frame only.
410 */
411 if (in[FRAME_ARG])
412 {
413 if (!IsSeries(vectors) && frame != thisframe)
414 {
415 DXSetError(ERROR_BAD_PARAMETER,
416 "frame parameter given is out of sequence");
417 goto error;
418 }
419
420 if (frame > head)
421 frame = head;
422 }
423 else if (IsSeries(vectors))
424 frame = head;
425 else
426 frame = thisframe;
427
428 /*
429 * The following loop updates the frame and generates streaks.
430 */
431 for (f = thisframe; f <= frame; f++)
432 {
433 if (! GetFrameData(vectors, curls, cstreaks))
434 goto error;
435
436 if (! Streaklines(cstreaks, c, curlFlag))
437 goto error;
438 }
439
440 if (!cachedObject)
441 if (NULL == (cachedObject = SetCacheObject(in, cstreaks)))
442 goto error;
443
444 past_head:
445
446 streak = cstreaks->streaks;
447 for (i = 0; i < cstreaks->nStreaks; i++, streak++)
448 {
449 Series s = NULL;
450 Field f;
451 int j, k;
452 float t;
453
454 k = 0;
455 for (j = 0; j <= frame; j++)
456 {
457 f = (Field)DXGetSeriesMember(streak->streak, j, &t);
458 if (!f)
459 break;
460
461 if (! DXEmptyField(f))
462 {
463 if (s == NULL)
464 {
465 s = DXNewSeries();
466 if (! s)
467 goto error;
468 }
469
470 if (! DXSetSeriesMember(s, k++, t, (Object)f))
471 {
472 DXDelete((Object)s);
473 goto error;
474 }
475 }
476 }
477
478 if (s)
479 {
480 if (group == NULL)
481 {
482 group = DXNewGroup();
483 if (! group)
484 {
485 DXDelete((Object)s);
486 goto error;
487 }
488 }
489
490 if (! DXSetMember(group, NULL, (Object)s))
491 {
492 DXDelete((Object)s);
493 goto error;
494 }
495 }
496 }
497
498 if (! group)
499 out[0] = (Object)DXEndField(DXNewField());
500 else
501 {
502 DXSetFloatAttribute((Object)group, "fuzz", 4);
503 out[0] = (Object)group;
504 }
505
506 DXDelete(cachedObject);
507
508 return OK;
509
510 no_streaks:
511 DXDelete(cachedObject);
512 DXDelete((Object)group);
513 DXDelete((Object)starts);
514 DXDelete((Object)times);
515 out[0] = NULL;
516 return OK;
517
518 error:
519 DXDelete(cachedObject);
520 DXDelete((Object)starts);
521 DXDelete((Object)times);
522 DXDelete((Object)group);
523 out[0] = NULL;
524 return ERROR;
525 }
526
527 typedef struct streakargs
528 {
529 CacheObject cstreak;
530 int i;
531 float c;
532 } StreakArgs;
533
534 typedef struct frameargs
535 {
536 int curlFlag;
537 Interpolator i0, i1;
538 float t0, t1;
539 CacheStreak streak;
540 } FrameArgs;
541
542 static Error
Streaklines(CacheObject cstreak,float c,int curlFlag)543 Streaklines(CacheObject cstreak, float c, int curlFlag)
544 {
545 int i, n;
546 StreakArgs stask;
547 FrameArgs ftask;
548 int nP;
549 Interpolator curlMap0 = NULL;
550 Interpolator curlMap1 = NULL;
551 MultiGrid cObject0 = NULL;
552 MultiGrid cObject1 = NULL;
553
554
555 /*
556 * If this is the 0-th frame, put an empty field into the streaks
557 * Otherwise, put any necessary streaks.
558 */
559 if (cstreak->frame == 0)
560 {
561 CacheStreak streak = cstreak->streaks;
562 Field eField = DXEndField(DXNewField());
563
564 for (i = 0; i < cstreak->nStreaks; i++, streak++)
565 {
566 if (! DXSetSeriesMember(streak->streak, 0, cstreak->t1,
567 (Object)eField))
568 goto error;
569 }
570 }
571 else
572 {
573 if (cstreak->t1 <- cstreak->t0)
574 cstreak->t1 = cstreak->t0;
575
576 nP = DXProcessors(0);
577
578 n = cstreak->nStreaks;
579 if (n < 1)
580 return OK;
581
582 if (n == 1 || nP == 1)
583 {
584 for (i = 0; i < n; i++)
585 {
586 if (! Streakline(cstreak, i, c))
587 goto error;
588 }
589 }
590 else
591 {
592 stask.cstreak = cstreak;
593 stask.c = c;
594
595 if (! DXCreateTaskGroup())
596 goto error;
597
598 for (i = 0; i < n; i++)
599 {
600 stask.i = i;
601
602 if (! DXAddTask(StreakTask, (Pointer)&stask, sizeof(stask), 1.0))
603 {
604 DXAbortTaskGroup();
605 goto error;
606 }
607 }
608
609 if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
610 goto error;
611 }
612
613 if (cstreak->nDim == 3)
614 {
615 if (curlFlag)
616 {
617 if (cstreak->c0)
618 curlMap0 = cstreak->c0;
619 else
620 {
621 int i;
622
623 cObject0 = DXNewMultiGrid();
624 if (! cObject0)
625 goto error;
626
627 for (i = 0; i < cstreak->vf0->nmembers; i++)
628 {
629 VectorGrp vg = cstreak->vf0->members[i];
630
631 if (! (*(vg->CurlMap))(vg, cObject0))
632 goto error;
633 }
634
635 curlMap0 = DXNewInterpolator((Object)cObject0,
636 INTERP_INIT_IMMEDIATE, -1.0);
637
638 if (! curlMap0)
639 goto error;
640 }
641
642 DXReference((Object)curlMap0);
643
644 if (cstreak->c1)
645 curlMap1 = cstreak->c1;
646 else
647 {
648 int i;
649
650 cObject1 = DXNewMultiGrid();
651 if (! cObject1)
652 goto error;
653
654 for (i = 0; i < cstreak->vf1->nmembers; i++)
655 {
656 VectorGrp vg = cstreak->vf1->members[i];
657
658 if (! (*(vg->CurlMap))(vg, cObject1))
659 goto error;
660 }
661
662 curlMap1 = DXNewInterpolator((Object)cObject1,
663 INTERP_INIT_IMMEDIATE, 0.0);
664
665 if (! curlMap1)
666 goto error;
667 }
668
669 DXReference((Object)curlMap1);
670
671 if (! curlMap0 || ! curlMap1)
672 {
673 DXSetError(ERROR_INTERNAL, "missing curl map");
674 goto error;
675 }
676
677 ftask.curlFlag = 1;
678 ftask.i0 = curlMap0;
679 ftask.i1 = curlMap1;
680 }
681 else
682 {
683 ftask.curlFlag = 0;
684 ftask.i0 = NULL;
685 ftask.i1 = NULL;
686 }
687
688 ftask.t0 = cstreak->t0;
689 ftask.t1 = cstreak->t1;
690
691 if (n == 1 || nP == 1)
692 {
693 for (i = 0; i < n; i++)
694 {
695 ftask.streak = cstreak->streaks + i;
696 if (! TraceFrame((Pointer)&ftask))
697 goto error;
698 }
699 }
700 else
701 {
702 if (! DXCreateTaskGroup())
703 goto error;
704
705 for (i = 0; i < n; i++)
706 {
707 ftask.streak = cstreak->streaks + i;
708 if (! DXAddTask(TraceFrame, (Pointer)&ftask,
709 sizeof(ftask), 1.0))
710 {
711 DXAbortTaskGroup();
712 goto error;
713 }
714 }
715
716 if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
717 goto error;
718 }
719
720 DXDelete((Object)curlMap0);
721 DXDelete((Object)curlMap1);
722 curlMap0 = curlMap1 = NULL;
723 }
724 }
725
726 return OK;
727
728 error:
729 if (cObject0)
730 DXDelete((Object)cObject0);
731 if (cObject1)
732 DXDelete((Object)cObject1);
733 if (curlMap0)
734 DXDelete((Object)curlMap0);
735 if (curlMap1)
736 DXDelete((Object)curlMap1);
737 return ERROR;
738 }
739
740 static Error
StreakTask(Pointer p)741 StreakTask(Pointer p)
742 {
743 StreakArgs *t = (StreakArgs *)p;
744
745 return Streakline(t->cstreak, t->i, t->c);
746 }
747
748 static Error
Streakline(CacheObject cstreak,int which,float c)749 Streakline(CacheObject cstreak, /* accumulated streak info */
750 int which, /* which to operate on */
751 float c) /* step criterion */
752 {
753 double elapsedtime = 0;
754 float head = cstreak->t1;
755 CacheStreak streak = cstreak->streaks + which;
756 int nDim = cstreak->nDim;
757 Stream streakBuf = NULL;
758 POINT_TYPE pb0[3], pb1[3];
759 VECTOR_TYPE vb0[3], vb1[3];
760 POINT_TYPE *p0, *p1, *ptmp;
761 VECTOR_TYPE *v0, *v1, *vtmp;
762 double t, t2;
763 int done, keepPoint;
764 Field field = NULL;
765 int i, totKnt;
766 StreakFlags sFlag;
767 StreakVars svars;
768 int zeroVector;
769
770 memset(&svars, 0, sizeof(svars));
771
772 /*
773 * Double buffer pointers
774 */
775 p0 = pb0; p1 = pb1;
776 v0 = vb0; v1 = vb1;
777
778 /*
779 * If streak has already died, well, OK.
780 */
781 if (streak->status == STREAK_DEAD)
782 return OK;
783
784 /*
785 * If streak is awaiting birth and its start time is after the current
786 * interval, well, OK. But got to put an empty field marker in.
787 */
788 if (streak->status == STREAK_PENDING && streak->time > head)
789 goto streak_done;
790
791 /*
792 * If streak is alive, then we need to get its current termination point.
793 * If the streak is newborn, then the point is already in the streak
794 * header. Otherwise, we need to get it from the fragment already
795 * computed.
796 */
797 if (streak->status == STREAK_ALIVE)
798 {
799 float fbuf[3];
800 GetTail((Object)streak->streak,"positions",(Pointer)&fbuf[0]);
801 for (i = 0; i < nDim; i++)
802 streak->point[i] = fbuf[i];
803 }
804 else
805 streak->status = STREAK_ALIVE;
806
807 for (i = 0; i < nDim; i++)
808 p0[i] = streak->point[i];
809
810 if (streak->time > cstreak->t0)
811 elapsedtime = streak->time;
812 else
813 elapsedtime = cstreak->t0;
814
815 if (! SetupStreakVars(&svars, cstreak))
816 goto error;
817
818 nDim = cstreak->nDim;
819
820 /*
821 * If we can't find the start, we can't extend the streak.
822 */
823 if (! Streak_FindElement(&svars, p0))
824 {
825 streak->status = STREAK_DEAD;
826 goto streak_done;
827 }
828
829 if (! Streak_Interpolate(&svars, p0, elapsedtime, v0))
830 goto error;
831
832 streakBuf = NewStreakBuf(nDim);
833 if (! streakBuf)
834 goto error;
835
836 sFlag = AddPointToStreak(streakBuf, p0, v0, elapsedtime);
837 if (sFlag == STREAK_ERROR)
838 goto error;
839
840 /*
841 * While streak is passing from partition to partition...
842 */
843 done = 0; totKnt = 1;
844 while (!done)
845 {
846 zeroVector = ZeroVector(v0, nDim);
847
848 if (zeroVector)
849 {
850 /*
851 * Interpolate this point at t1. If its zero there too,
852 * then it'll be zero everywhere in between (due to linear
853 * interpolation in time.
854 */
855
856 if (! Streak_Interpolate(&svars, p0, svars.t1, v1))
857 goto error;
858
859 if (ZeroVector(v1, nDim))
860 {
861 sFlag = AddPointToStreak(streakBuf, p0, v0, svars.t1);
862 if (sFlag == STREAK_ERROR)
863 goto error;
864
865 streak->status = STREAK_ALIVE;
866 goto streak_done;
867 }
868 else
869 {
870 /*
871 * If it isn't zero at t1, then take a small step
872 * in time to get off zero.
873 */
874 while (zeroVector)
875 {
876 elapsedtime += 0.025 * (svars.t1 - elapsedtime);
877
878 if (! Streak_Interpolate(&svars, p0, elapsedtime, v0))
879 goto error;
880
881 zeroVector = ZeroVector(v0, nDim);
882 }
883 }
884 }
885
886 if (! Streak_StepTime(&svars, c, v0, &t))
887 goto error;
888
889 /*
890 * If the time step exceeds the head, truncate it and mark
891 * the streak complete.
892 */
893 if (elapsedtime + t > svars.t1)
894 {
895 t = svars.t1 - elapsedtime;
896 streak->status = STREAK_ALIVE;
897 done = 1;
898 }
899
900 /*
901 * Get the point a time step away in the direction of the
902 * vector
903 */
904 if (! zeroVector)
905 for (i = 0; i < nDim; i++)
906 {
907 p1[i] = p0[i] + t*v0[i];
908 if (p1[i] == p0[i])
909 v0[i] = 0.0;
910 }
911
912 /*
913 * If the new point is not in the current element, find the point
914 * at which it exits the element and iterponlate the point there.
915 * Then move on to the neighbor, if one exists.
916 *
917 * Otherwise, just interpolate the endpoint.
918 */
919 if (! Streak_InsideWeights(&svars, p1))
920 {
921 int found;
922
923 /*
924 * If the new sample isn't in the previous element, find the
925 * exit point of the streak and interpolate in the current
926 * element there.
927 */
928 if (! Streak_FindBoundary(&svars, p0, v0, &t2))
929 goto error;
930
931 elapsedtime += t2;
932
933 for (i = 0; i < nDim; i++)
934 {
935 double d = ((double)p0[i]) + ((double)t2)*((double)v0[i]);
936 p1[i] = d;
937 }
938
939 if (! Streak_FaceWeights(&svars, p1))
940 goto error;
941
942 if (! Streak_Interpolate(&svars, p1, elapsedtime, v1))
943 goto error;
944
945 found = Streak_Neighbor(&svars, v1);
946 if (found == -1)
947 goto error;
948 if (! found)
949 found = Streak_FindMultiGridContinuation(&svars, p0);
950 if (! found)
951 {
952 streak->status = STREAK_DEAD;
953 done = 1;
954 }
955 else
956 done = 0;
957
958 if (done || t2 > 0)
959 keepPoint = 1;
960 else
961 keepPoint = 0;
962
963 t = t2;
964 }
965 else
966 {
967 keepPoint = 1;
968 elapsedtime += t;
969
970 if (! Streak_Interpolate(&svars, p1, elapsedtime, v1))
971 goto error;
972 }
973
974 if (keepPoint && t > 0.0)
975 {
976 sFlag = AddPointToStreak(streakBuf, p1, v1, elapsedtime);
977 if (sFlag == STREAK_ERROR)
978 goto error;
979 else if (sFlag == STREAK_FULL)
980 done = 1;
981 }
982
983 if (++totKnt > TOT_MAX)
984 {
985 DXWarning("possible infinite loop detected");
986 done = 1;
987 }
988
989 ptmp = p0;
990 p0 = p1;
991 p1 = ptmp;
992
993 vtmp = v0;
994 v0 = v1;
995 v1 = vtmp;
996 }
997
998 streak_done:
999
1000 if (streakBuf)
1001 {
1002 field = StreakToField(streakBuf);
1003 if (! field)
1004 goto error;
1005
1006 DestroyStreakBuf(streakBuf);
1007 streakBuf = NULL;
1008 }
1009 else
1010 field = DXEndField(DXNewField());
1011
1012 if (! DXSetSeriesMember(streak->streak, cstreak->frame,
1013 svars.t0, (Object)field))
1014 goto error;
1015
1016 FreeStreakVars(&svars);
1017
1018 return OK;
1019
1020 error:
1021 FreeStreakVars(&svars);
1022 DestroyStreakBuf(streakBuf);
1023 DXDelete((Object)field);
1024
1025 return ERROR;
1026 }
1027
1028 static Error
TraceFrame(Pointer ptr)1029 TraceFrame(Pointer ptr)
1030 {
1031 Interpolator i0 = ((FrameArgs *)ptr)->i0;
1032 Interpolator i1 = ((FrameArgs *)ptr)->i1;
1033 float t0 = ((FrameArgs *)ptr)->t0;
1034 float t1 = ((FrameArgs *)ptr)->t1;
1035 CacheStreak cs = ((FrameArgs *)ptr)->streak;
1036 int flag = ((FrameArgs *)ptr)->curlFlag;
1037 Series s = cs->streak;
1038 Array cArray0 = NULL;
1039 Array cArray1 = NULL;
1040 Field field = NULL;
1041 Object next;
1042 Field prev;
1043 Array nArray = NULL, bArray = NULL;
1044 Array vArray, tArray, pArray;
1045 float *points, *vecarr, *curls0=NULL, *curls1=NULL,
1046 *normals, *binormals, *time;
1047 int i, nVectors, nDim, seg;
1048 Object dattr = NULL;
1049 Vector fn, fb, ft, vectors;
1050
1051 field = NULL; prev = NULL;
1052 for (seg = 0; NULL != (next = DXGetSeriesMember(s, seg, NULL)); seg++)
1053 {
1054 if (field && ! DXEmptyField(field))
1055 prev = field;
1056 field = (Field)next;
1057 }
1058
1059 if (DXEmptyField(field))
1060 return OK;
1061
1062 pArray = (Array)DXGetComponentValue(field, "positions");
1063 if (! pArray)
1064 goto error;
1065
1066 points = (float *)DXGetArrayData(pArray);
1067
1068 if (flag)
1069 {
1070 i0 = (Interpolator)DXCopy((Object)i0, COPY_STRUCTURE);
1071 if (! i0)
1072 goto error;
1073
1074 cArray0 = DXMapArray(pArray, i0, NULL);
1075 if (! cArray0)
1076 goto error;
1077
1078 DXDelete((Object)i0);
1079
1080 i1 = (Interpolator)DXCopy((Object)i1, COPY_STRUCTURE);
1081 if (! i1)
1082 goto error;
1083
1084 cArray1 = DXMapArray(pArray, i1, NULL);
1085 if (! cArray1)
1086 goto error;
1087
1088 DXDelete((Object)i1);
1089
1090 curls0 = (float *)DXGetArrayData(cArray0);
1091 curls1 = (float *)DXGetArrayData(cArray1);
1092 if (! curls0 || ! curls1)
1093 goto error;
1094 }
1095
1096 vArray = (Array)DXGetComponentValue(field, "data");
1097 if (! vArray)
1098 goto error;
1099
1100 DXGetArrayInfo(vArray, &nVectors, NULL, NULL, NULL, &nDim);
1101 if (nDim < 3)
1102 goto done;
1103
1104 vecarr = (float *)DXGetArrayData(vArray);
1105 if (! vecarr)
1106 goto error;
1107 vectors.x = vecarr[0];
1108 vectors.y = vecarr[1];
1109 vectors.z = vecarr[2];
1110
1111 tArray = (Array)DXGetComponentValue(field, "time");
1112 if (! tArray)
1113 goto error;
1114
1115 time = (float *)DXGetArrayData(tArray);
1116 if (! time)
1117 goto error;
1118
1119 nArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
1120 if (! nArray)
1121 goto error;
1122
1123 if (! DXAddArrayData(nArray, 0, nVectors, NULL))
1124 goto error;
1125
1126 normals = (float *)DXGetArrayData(nArray);
1127 if (! normals)
1128 goto error;
1129
1130 bArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
1131 if (! bArray)
1132 goto error;
1133
1134 if (! DXAddArrayData(bArray, 0, nVectors, NULL))
1135 goto error;
1136
1137 binormals = (float *)DXGetArrayData(bArray);
1138 if (! binormals)
1139 goto error;
1140
1141 /*
1142 * If this is not the first segment of the stream, we need to get
1143 * the propagated frame of reference and initial vertex normal and
1144 * binormal from the cached stuff. Otherwise, we initialize the
1145 * process.
1146 */
1147 if (prev)
1148 {
1149 fn.x = cs->norm[0]; fn.y = cs->norm[1]; fn.z = cs->norm[2];
1150 fb.x = cs->binorm[0]; fb.y = cs->binorm[1]; fb.z = cs->binorm[2];
1151 ft.x = cs->tan[0]; ft.y = cs->tan[1]; ft.z = cs->tan[2];
1152
1153 GetTail((Object)prev, "normals", (Pointer)normals);
1154 GetTail((Object)prev, "binormals", (Pointer)binormals);
1155 }
1156 else
1157 {
1158 if (! InitFrame(&vectors, &fn, &fb))
1159 goto error;
1160
1161 vecarr[0] = vectors.x; vecarr[1] = vectors.y; vecarr[2] = vectors.z;
1162
1163 ft.x = vectors.x; ft.y = vectors.y; ft.z = vectors.z;
1164
1165 _dxfvector_normalize_3D(&fn, &fn);
1166 _dxfvector_normalize_3D(&fb, &fb);
1167 _dxfvector_normalize_3D(&ft, &ft);
1168
1169 normals[0] = fn.x; normals[1] = fn.y; normals[2] = fn.z;
1170 binormals[0] = fb.x; binormals[1] = fb.y; binormals[2] = fb.z;
1171 }
1172
1173 for (i = 1; i < nVectors; i++)
1174 {
1175 points += 3; time += 1; vecarr += 3; normals += 3; binormals += 3;
1176
1177 if (flag)
1178 {
1179 float avcurl[3], pcurl[3], ncurl[3], d;
1180
1181 d = (*time - t0) / (t1 - t0);
1182
1183 pcurl[0] = curls0[0] + d*(curls1[0] - curls0[0]);
1184 pcurl[1] = curls0[1] + d*(curls1[1] - curls0[1]);
1185 pcurl[2] = curls0[2] + d*(curls1[2] - curls0[2]);
1186
1187 ncurl[0] = curls0[3] + d*(curls1[3] - curls0[3]);
1188 ncurl[1] = curls0[4] + d*(curls1[4] - curls0[4]);
1189 ncurl[2] = curls0[5] + d*(curls1[5] - curls0[5]);
1190
1191 avcurl[0] = 0.5 * (pcurl[0] + ncurl[0]);
1192 avcurl[1] = 0.5 * (pcurl[1] + ncurl[1]);
1193 avcurl[2] = 0.5 * (pcurl[2] + ncurl[2]);
1194
1195 curls0 += 3;
1196 curls1 += 3;
1197
1198 if (! UpdateFrame(points, vecarr, avcurl, time, normals,
1199 binormals, &fn, &fb, &ft))
1200 goto error;
1201 }
1202 else
1203 {
1204 if (! UpdateFrame(points, vecarr, NULL, time, normals,
1205 binormals, &fn, &fb, &ft))
1206 goto error;
1207 }
1208 }
1209
1210 cs->norm[0] = fn.x; cs->norm[1] = fn.y; cs->norm[2] = fn.z;
1211 cs->binorm[0] = fb.x; cs->binorm[1] = fb.y; cs->binorm[2] = fb.z;
1212 cs->tan[0] = ft.x; cs->tan[1] = ft.y; cs->tan[2] = ft.z;
1213
1214 dattr = (Object)DXNewString("positions");
1215 if (! dattr)
1216 goto error;
1217 DXReference(dattr);
1218
1219 if (! DXSetComponentValue(field, "normals", (Object)nArray))
1220 goto error;
1221 nArray = NULL;
1222
1223 if (! DXSetComponentAttribute(field, "normals", "dep", dattr))
1224 goto error;
1225
1226 if (! DXSetComponentValue(field, "binormals", (Object)bArray))
1227 goto error;
1228 bArray = NULL;
1229
1230 if (! DXSetComponentAttribute(field, "binormals", "dep", dattr))
1231 goto error;
1232
1233 DXDelete(dattr);
1234 DXDelete((Object)cArray0);
1235 DXDelete((Object)cArray1);
1236
1237 done:
1238 return OK;
1239
1240 error:
1241 DXDelete((Object)nArray);
1242 DXDelete((Object)bArray);
1243 DXDelete((Object)cArray0);
1244 DXDelete((Object)cArray1);
1245
1246 return ERROR;
1247 }
1248
1249 static Array
Starts(Object starts,Object vectors,int nD)1250 Starts(Object starts, Object vectors, int nD)
1251 {
1252 Array array = NULL;
1253
1254 if (starts)
1255 {
1256 Object xstarts = DXApplyTransform(starts, NULL);
1257 if (! xstarts)
1258 goto error;
1259
1260 if (! StartsRecurse(xstarts, &array))
1261 {
1262 if (starts != xstarts)
1263 DXDelete((Object)xstarts);
1264 goto error;
1265 }
1266
1267 if (starts != xstarts)
1268 DXDelete((Object)xstarts);
1269 }
1270 else if (vectors)
1271 {
1272 Point box[8];
1273 float pt[3];
1274
1275 if (! DXBoundingBox(vectors, box))
1276 return NULL;
1277
1278 pt[0] = (box[0].x + box[7].x) / 2.0;
1279 pt[1] = (box[0].y + box[7].y) / 2.0;
1280 pt[2] = (box[0].z + box[7].z) / 2.0;
1281
1282 array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nD);
1283 if (! array)
1284 goto error;
1285
1286 if (! DXAddArrayData(array, 0, 1, (Pointer)pt))
1287 goto error;
1288
1289 }
1290 else
1291 {
1292 float *ptr;
1293 int i;
1294
1295 array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nD);
1296 if (! array)
1297 goto error;
1298
1299 if (! DXAddArrayData(array, 0, 1, NULL))
1300 goto error;
1301
1302 ptr = (float *)DXGetArrayData(array);
1303 if (! ptr)
1304 goto error;
1305
1306 for (i = 0; i < nD; i++)
1307 *ptr++ = 0.0;
1308 }
1309
1310 return array;
1311
1312 error:
1313 DXDelete((Object)array);
1314 return NULL;
1315 }
1316
1317 static Array
Times(Object times,Object vectors)1318 Times(Object times, Object vectors)
1319 {
1320 Array array = NULL;
1321
1322 if (times)
1323 {
1324 if (! StartsRecurse(times, &array))
1325 goto error;
1326 }
1327 else if (vectors)
1328 {
1329 float time;
1330
1331 if (IsSeries(vectors))
1332 DXGetSeriesMember((Series)vectors, 0, &time);
1333 else
1334 {
1335 Object attr = DXGetAttribute(vectors, "series position");
1336
1337 if (attr == NULL)
1338 DXSetError(ERROR_MISSING_DATA,
1339 "series position attribute on data");
1340
1341 if (! DXExtractFloat(attr, &time))
1342 DXSetError(ERROR_DATA_INVALID,
1343 "series position attribute on data");
1344 }
1345
1346 array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0);
1347 if (! array)
1348 goto error;
1349
1350 if (! DXAddArrayData(array, 0, 1, (Pointer)&time))
1351 goto error;
1352
1353 }
1354 else
1355 {
1356 float *ptr;
1357
1358 array = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1);
1359 if (! array)
1360 goto error;
1361
1362 if (! DXAddArrayData(array, 0, 1, NULL))
1363 goto error;
1364
1365 ptr = (float *)DXGetArrayData(array);
1366 if (! ptr)
1367 goto error;
1368
1369 *ptr = 0.0;
1370 }
1371
1372 return array;
1373
1374 error:
1375 DXDelete((Object)array);
1376 return NULL;
1377 }
1378
1379 static Error
AlignStartPtsAndTimes(Array starts,Array times)1380 AlignStartPtsAndTimes(Array starts, Array times)
1381 {
1382 int nStarts, nTimes;
1383
1384 DXGetArrayInfo(starts, &nStarts, NULL, NULL, NULL, NULL);
1385 DXGetArrayInfo(times, &nTimes, NULL, NULL, NULL, NULL);
1386
1387 if (nStarts == nTimes)
1388 return OK;
1389 else if (nStarts == 1)
1390 {
1391 int itemSize = DXGetItemSize(starts);
1392 int i;
1393 unsigned char *src, *dst;
1394
1395 if (! DXAddArrayData(starts, 1, nTimes-1, NULL))
1396 goto error;
1397
1398 src = (unsigned char *)DXGetArrayData(starts);
1399 dst = src + itemSize;
1400 for (i = 0; i < nTimes; i++, dst += itemSize)
1401 memcpy(dst, src, itemSize);
1402 }
1403 else if (nTimes == 1)
1404 {
1405 int itemSize = DXGetItemSize(times);
1406 int i;
1407 unsigned char *src, *dst;
1408
1409 if (! DXAddArrayData(times, 1, nStarts-1, NULL))
1410 goto error;
1411
1412 src = (unsigned char *)DXGetArrayData(times);
1413 dst = src + itemSize;
1414 for (i = 0; i < nStarts; i++, dst += itemSize)
1415 memcpy(dst, src, itemSize);
1416 }
1417 else
1418 {
1419 DXSetError(ERROR_BAD_PARAMETER, "#11160",
1420 "start positions", nStarts, "start times", nTimes);
1421 goto error;
1422 }
1423
1424 return OK;
1425
1426 error:
1427 return ERROR;
1428 }
1429
1430 static Error
StartsRecurse(Object o,Array * a)1431 StartsRecurse(Object o, Array *a)
1432 {
1433 Object c;
1434 int i;
1435
1436 if (DXGetObjectClass(o) == CLASS_ARRAY)
1437 return StartsAddPoints(a, (Array)o);
1438 else if (DXGetObjectClass(o) == CLASS_FIELD)
1439 {
1440 if (! DXEmptyField((Field)o))
1441 return StartsAddPoints(a,
1442 (Array)DXGetComponentValue((Field)o, "positions"));
1443 else
1444 return OK;
1445 }
1446 else if (DXGetObjectClass(o) == CLASS_GROUP)
1447 {
1448 i = 0;
1449 while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL)))
1450 {
1451 if (! StartsRecurse(c, a))
1452 goto error;
1453 }
1454 return OK;
1455 }
1456 else
1457 {
1458 DXSetError(ERROR_BAD_PARAMETER, "starts");
1459 return ERROR;
1460 }
1461
1462 error:
1463 DXDelete((Object)*a);
1464 *a = NULL;
1465 return ERROR;
1466 }
1467
1468 static Error
StartsAddPoints(Array * dst,Array src)1469 StartsAddPoints(Array *dst, Array src)
1470 {
1471 int nSrc, nDst, rank, shape[32];
1472
1473 if (! src || DXGetObjectClass((Object)src) != CLASS_ARRAY)
1474 return ERROR;
1475
1476 if (*dst == NULL)
1477 {
1478 DXGetArrayInfo(src, &nSrc, NULL, NULL, &rank, shape);
1479 if (nSrc == 0)
1480 return OK;
1481
1482 *dst = DXNewArrayV(TYPE_FLOAT, CATEGORY_REAL, rank, shape);
1483 if (! *dst)
1484 return ERROR;
1485
1486 if (rank == 0)
1487 shape[0] = 1;
1488
1489 if (! DXAddArrayData(*dst, 0, nSrc, NULL))
1490 return ERROR;
1491
1492 if (! DXExtractParameter((Object)src, TYPE_FLOAT, shape[0], nSrc,
1493 DXGetArrayData(*dst)))
1494 return ERROR;
1495 }
1496 else
1497 {
1498 DXGetArrayInfo(*dst, &nDst, NULL, NULL, &rank, shape);
1499 DXGetArrayInfo(src, &nSrc, NULL, NULL, NULL, NULL);
1500 if (nSrc == 0)
1501 return OK;
1502
1503 if (! DXAddArrayData(*dst, nDst, nSrc, NULL))
1504 return ERROR;
1505
1506 if (rank == 0)
1507 shape[0] = 1;
1508
1509 if (! DXExtractParameter((Object)src, TYPE_FLOAT, shape[0], nSrc,
1510 (Pointer)(((float *)DXGetArrayData(*dst))+nDst*shape[0])))
1511 return ERROR;
1512 }
1513
1514 return OK;
1515 }
1516
1517 static Error
DestroyVectorField(VectorField vf)1518 DestroyVectorField(VectorField vf)
1519 {
1520 if (vf)
1521 {
1522 int i;
1523 for (i = 0; i < vf->nmembers; i++)
1524 if (vf->members[i])
1525 (*(vf->members[i]->Delete))(vf->members[i]);
1526
1527 DXFree((Pointer)vf->members);
1528 DXFree((Pointer)vf);
1529 }
1530 return OK;
1531 }
1532
1533 static VectorField
InitVectorField(Object vfo)1534 InitVectorField(Object vfo)
1535 {
1536 char *str;
1537 VectorField vf = NULL;
1538 Class class;
1539
1540 vf = (struct vectorField *)DXAllocateZero(sizeof(struct vectorField));
1541 if (! vf)
1542 goto error;
1543
1544 class = DXGetObjectClass(vfo);
1545 if (class == CLASS_GROUP)
1546 class = DXGetGroupClass((Group)vfo);
1547
1548 switch(class)
1549 {
1550 case CLASS_MULTIGRID:
1551 {
1552 int i, n;
1553 char *str;
1554
1555 if (! DXGetMemberCount((Group)vfo, &n))
1556 goto error;
1557
1558 vf->nmembers = n;
1559
1560 vf->members = (VectorGrp *)DXAllocate(n*sizeof(VectorGrp));
1561 if (! vf->members)
1562 goto error;
1563
1564 for (i = 0; i < n; i++)
1565 {
1566 Object mo = DXGetEnumeratedMember((Group)vfo, i, NULL);
1567 if (!mo)
1568 goto error;
1569
1570 class = DXGetObjectClass(mo);
1571 if (class == CLASS_GROUP)
1572 class = DXGetGroupClass((Group)mo);
1573
1574 if (class != CLASS_FIELD && class != CLASS_COMPOSITEFIELD)
1575 {
1576 DXSetError(ERROR_DATA_INVALID, "invalid multigrid member");
1577 goto error;
1578 }
1579
1580 if (! GetElementType(mo, &str))
1581 {
1582 DXSetError(ERROR_DATA_INVALID, "vectors");
1583 goto error;
1584 }
1585
1586 switch (IsRegular(mo))
1587 {
1588 case 1:
1589 vf->members[i] = _dxfReg_InitVectorGrp(mo, str);
1590 break;
1591 case 0:
1592 vf->members[i] = _dxfIrreg_InitVectorGrp(mo, str);
1593 break;
1594
1595 default:
1596 DXSetError(ERROR_DATA_INVALID, "vectors");
1597 goto error;
1598 }
1599
1600 if (! vf->members[i])
1601 goto error;
1602 }
1603 }
1604 break;
1605
1606 case CLASS_FIELD:
1607 case CLASS_COMPOSITEFIELD:
1608 {
1609 vf->nmembers = 1;
1610
1611 vf->members = (VectorGrp *)DXAllocate(sizeof(VectorGrp));
1612 if (! vf->members)
1613 goto error;
1614
1615 if (! GetElementType(vfo, &str))
1616 {
1617 DXSetError(ERROR_DATA_INVALID, "vectors");
1618 goto error;
1619 }
1620
1621 switch (IsRegular(vfo))
1622 {
1623 case 1:
1624 vf->members[0] = _dxfReg_InitVectorGrp(vfo, str);
1625 break;
1626 case 0:
1627 vf->members[0] = _dxfIrreg_InitVectorGrp(vfo, str);
1628 break;
1629
1630 default:
1631 DXSetError(ERROR_DATA_INVALID, "vectors");
1632 goto error;
1633 }
1634
1635 if (! vf->members[0])
1636 goto error;
1637 }
1638 break;
1639
1640 default:
1641 {
1642 DXSetError(ERROR_DATA_INVALID,
1643 "invalid object encountered in vector field");
1644 goto error;
1645 }
1646 }
1647
1648 vf->nDim = vf->members[0]->nDim;
1649
1650 return vf;
1651
1652 error:
1653 DestroyVectorField(vf);
1654 return NULL;
1655 }
1656
1657 static Error
GetElementType(Object o,char ** str)1658 GetElementType(Object o, char **str)
1659 {
1660 Object c;
1661
1662 *str = NULL;
1663
1664 if (DXGetObjectClass(o) == CLASS_FIELD)
1665 {
1666 if (DXEmptyField((Field)o))
1667 return OK;
1668
1669 if (! DXGetComponentValue((Field)o, "connections"))
1670 {
1671 DXSetError(ERROR_DATA_INVALID, "data has no connections");
1672 return ERROR;
1673 }
1674
1675 c = DXGetComponentAttribute((Field)o, "connections", "element type");
1676 if (! c)
1677 {
1678 DXSetError(ERROR_MISSING_DATA, "element type attribute");
1679 return ERROR;
1680 }
1681
1682 if (DXGetObjectClass(c) != CLASS_STRING)
1683 {
1684 DXSetError(ERROR_DATA_INVALID, "element type attribute");
1685 return ERROR;
1686 }
1687
1688 *str = DXGetString((String)c);
1689 return OK;
1690 }
1691 else if (DXGetObjectClass(o) == CLASS_GROUP)
1692 {
1693 int i;
1694 Error j = ERROR;
1695
1696 i = 0;
1697 while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL)))
1698 {
1699 j = GetElementType(c, str);
1700 if (j != ERROR)
1701 return OK;
1702 }
1703
1704 return ERROR;
1705 }
1706 else
1707 {
1708 DXSetError(ERROR_BAD_CLASS, "vector field");
1709 return ERROR;
1710 }
1711 }
1712
1713 static Stream
NewStreakBuf(int nDim)1714 NewStreakBuf(int nDim)
1715 {
1716 Stream s = NULL;
1717
1718 s = (Stream)DXAllocateZero(sizeof(struct stream));
1719 if (! s)
1720 goto error;
1721
1722 s->pArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nDim);
1723 s->tArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 0);
1724 s->vArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, nDim);
1725
1726 if (s->pArray == NULL || s->vArray == NULL || s->tArray == NULL)
1727 goto error;
1728
1729 s->arrayKnt = 0;
1730
1731 s->points = (float *)DXAllocate(nDim*STREAM_BUF_POINTS*sizeof(float));
1732 s->time = (float *)DXAllocate(STREAM_BUF_POINTS*sizeof(float));
1733 s->vectors = (float *)DXAllocate(nDim*STREAM_BUF_POINTS*sizeof(float));
1734 if (s->points == NULL || s->time == NULL || s->vectors == NULL)
1735 goto error;
1736
1737 s->nDim = nDim;
1738 s->bufKnt = 0;
1739
1740 return s;
1741
1742 error:
1743 DestroyStreakBuf(s);
1744
1745 return NULL;
1746 }
1747
1748 static void
DestroyStreakBuf(Stream s)1749 DestroyStreakBuf(Stream s)
1750 {
1751 if (s)
1752 {
1753 DXDelete((Object)s->pArray);
1754 DXDelete((Object)s->vArray);
1755 DXDelete((Object)s->tArray);
1756 DXFree((Pointer)s->points);
1757 DXFree((Pointer)s->time);
1758 DXFree((Pointer)s->vectors);
1759 DXFree((Pointer)s);
1760 }
1761 }
1762
1763 static Error
FlushStreak(Stream s)1764 FlushStreak(Stream s)
1765 {
1766 if (! DXAddArrayData(s->pArray, s->arrayKnt, s->bufKnt, (Pointer)s->points))
1767 goto error;
1768
1769 if (! DXAddArrayData(s->vArray, s->arrayKnt, s->bufKnt, (Pointer)s->vectors))
1770 goto error;
1771
1772 if (! DXAddArrayData(s->tArray, s->arrayKnt, s->bufKnt, (Pointer)s->time))
1773 goto error;
1774
1775 s->arrayKnt += s->bufKnt;
1776 s->bufKnt = 0;
1777
1778 return OK;
1779
1780 error:
1781 return ERROR;
1782 }
1783
1784 static StreakFlags
AddPointToStreak(Stream s,POINT_TYPE * p,VECTOR_TYPE * v,double t)1785 AddPointToStreak(Stream s, POINT_TYPE *p, VECTOR_TYPE *v, double t)
1786 {
1787 if (s->bufKnt == STREAM_BUF_POINTS)
1788 if (! FlushStreak(s))
1789 goto error;
1790
1791 if (s->nDim == 2)
1792 {
1793 float *ptr;
1794
1795 ptr = s->points + s->bufKnt*s->nDim;
1796 *ptr++ = *p++;
1797 *ptr++ = *p++;
1798
1799 ptr = s->vectors + s->bufKnt*s->nDim;
1800 *ptr++ = *v++;
1801 *ptr++ = *v++;
1802 }
1803 else
1804 {
1805 float *ptr;
1806
1807 ptr = s->points + s->bufKnt*s->nDim;
1808 *ptr++ = *p++;
1809 *ptr++ = *p++;
1810 *ptr++ = *p++;
1811
1812 ptr = s->vectors + s->bufKnt*s->nDim;
1813 *ptr++ = *v++;
1814 *ptr++ = *v++;
1815 *ptr++ = *v++;
1816 }
1817
1818 s->time[s->bufKnt] = t;
1819
1820 s->bufKnt ++;
1821
1822 if (s->bufKnt + s->arrayKnt >= STREAK_MAX)
1823 return STREAK_FULL;
1824
1825 return STREAK_OK;
1826
1827 error:
1828 return STREAK_ERROR;
1829 }
1830
1831 static Field
StreakToField(Stream buf)1832 StreakToField(Stream buf)
1833 {
1834 float c[3];
1835 float d[3];
1836 Array a = NULL;
1837 Object dattr = NULL;
1838 Field field;
1839
1840 if (! FlushStreak(buf))
1841 goto error;
1842
1843 c[0] = 0.7; c[1] = 0.7; c[2] = 0.0;
1844 d[0] = 0.0; d[1] = 0.0; d[2] = 0.0;
1845
1846 dattr = (Object)DXNewString("positions");
1847 if (! dattr)
1848 goto error;
1849
1850 DXReference(dattr);
1851
1852 field = DXNewField();
1853 if (! field)
1854 goto error;
1855
1856 if (! DXSetComponentValue(field, "positions", (Object)buf->pArray))
1857 goto error;
1858 buf->pArray = NULL;
1859
1860 if (! DXSetComponentAttribute(field, "positions", "dep", dattr))
1861 goto error;
1862
1863 if (! DXSetComponentValue(field, "data", (Object)buf->vArray))
1864 goto error;
1865 buf->vArray = NULL;
1866
1867 if (! DXSetComponentAttribute(field, "data", "dep", dattr))
1868 goto error;
1869
1870 if (! DXSetComponentValue(field, "time", (Object)buf->tArray))
1871 goto error;
1872 buf->tArray = NULL;
1873
1874 if (! DXSetComponentAttribute(field, "time", "dep", dattr))
1875 goto error;
1876
1877 DXDelete(dattr);
1878
1879 a = DXMakeGridConnections(1, buf->arrayKnt);
1880 if (! a)
1881 goto error;
1882
1883 if (! DXSetComponentValue(field, "connections", (Object)a))
1884 goto error;
1885 a = NULL;
1886
1887 if (! DXSetComponentAttribute(field, "connections", "element type",
1888 (Object)DXNewString("lines")))
1889 goto error;
1890
1891 a = (Array)DXNewRegularArray(TYPE_FLOAT, 3,
1892 buf->arrayKnt, (Pointer)c, (Pointer)d);
1893 if (! a)
1894 goto error;
1895
1896 if (! DXSetComponentValue(field, "colors", (Object)a))
1897 goto error;
1898 a = NULL;
1899
1900 if (! DXSetIntegerAttribute((Object)field, "shade", 0))
1901 goto error;
1902
1903 if (! DXEndField(field))
1904 goto error;
1905
1906 return field;
1907
1908 error:
1909 DXDelete(dattr);
1910 return ERROR;
1911 }
1912
1913 static int
ZeroVector(VECTOR_TYPE * v,int n)1914 ZeroVector(VECTOR_TYPE *v, int n)
1915 {
1916 int i;
1917
1918 for (i = 0; i < n; i++)
1919 if (v[i] != 0.0) return 0;
1920
1921 return 1;
1922 }
1923
1924 static Error
InitFrame(Vector * v,Vector * n,Vector * b)1925 InitFrame(Vector *v, Vector *n, Vector *b)
1926 {
1927 n->x = n->y = n->y = 0.0;
1928
1929 if (v->x == 0.0 && v->y == 0.0 && v->z == 0.0)
1930 {
1931 n->x = 0.0; n->y = 0.0; n->z = 0.0;
1932 b->x = 0.0; b->y = 0.0; b->z = 0.0;
1933 return OK;
1934 }
1935
1936 if (v->x == 0.0 && v->y == 0.0 && v->z > 0.0)
1937 n->x = 1.0;
1938 else if (v->x == 0.0 && v->y == 0.0 && v->z < 0.0)
1939 n->x = -1.0;
1940 else if (v->x == 0.0 && v->y > 0.0 && v->z == 0.0)
1941 n->z = 1.0;
1942 else if (v->x == 0.0 && v->y < 0.0 && v->z == 0.0)
1943 n->z = -1.0;
1944 else if (v->x > 0.0 && v->y == 0.0 && v->z == 0.0)
1945 n->y = 1.0;
1946 else if (v->x < 0.0 && v->y == 0.0 && v->z == 0.0)
1947 n->y = -1.0;
1948 else
1949 {
1950 float d;
1951
1952 if (v->z != 0.0)
1953 {
1954 n->x = n->y = 1.0;
1955 n->z = -(v->x + v->y) / v->z;
1956 }
1957 else if (v->y != 0.0)
1958 {
1959 n->x = n->z = 1.0;
1960 n->y = -(v->x + v->z) / v->y;
1961 }
1962
1963 d = 1.0 / sqrt(n->x*n->x + n->y*n->y + n->z*n->z);
1964
1965 n->x *= d;
1966 n->y *= d;
1967 n->z *= d;
1968 }
1969
1970 _dxfvector_cross_3D(v, n, b);
1971 _dxfvector_normalize_3D(b, b);
1972
1973 return OK;
1974 }
1975
1976 /*
1977 * update frame of reference from prior vertex to current vertex based
1978 * on incuming and outgoing edge vectors and twist.
1979 */
1980
1981 #define VecMat(a, b, c) \
1982 { \
1983 float x = a->x; \
1984 float y = a->y; \
1985 float z = a->z; \
1986 \
1987 c->x = x*(b)[0] + y*(b)[3] + z*(b)[6]; \
1988 c->y = x*(b)[1] + y*(b)[4] + z*(b)[7]; \
1989 c->z = x*(b)[2] + y*(b)[5] + z*(b)[8]; \
1990 }
1991
1992 #define MatMat(a, b, c) \
1993 { \
1994 VecMat((a+0), (b), (c+0)); \
1995 VecMat((a+3), (b), (c+3)); \
1996 VecMat((a+6), (b), (c+6)); \
1997 }
1998
1999 static Error
UpdateFrame(float * p,float * v,float * c,float * t,float * n,float * b,Vector * fn,Vector * fb,Vector * ft)2000 UpdateFrame(float *p, float *v, float *c, float *t, float *n, float *b,
2001 Vector *fn, Vector *fb, Vector *ft)
2002 {
2003 float twist;
2004 float mTwist[9], mBend[9];
2005 float len, cA;
2006 int bend = 0;
2007 Vector v1, v2, vv, cross, nt;
2008
2009
2010 /*
2011 * If theres a curl, rotate by the curl around the prior segment vector.
2012 */
2013 if (c != NULL)
2014 {
2015 Vector step;
2016 Vector vc;
2017
2018 vc.x = c[0]; vc.y = c[1]; vc.z = c[2];
2019
2020 v1.x = p[0]; v1.y = p[1]; v1.z = p[2];
2021 v2.x = (p-3)[0]; v2.y = (p-3)[1]; v2.z = (p-3)[2];
2022 _dxfvector_subtract_3D(&v1, &v2, &step);
2023
2024 twist = 0.5 * _dxfvector_dot_3D(&step, &vc);
2025 if (twist != 0.0)
2026 {
2027 float fna[3], fba[3];
2028 fna[0] = fn->x; fna[1] = fn->y; fna[2] = fn->z;
2029 fba[0] = fb->x; fba[1] = fb->y; fba[2] = fb->z;
2030
2031 RotateAroundVector(ft, cos(twist), sin(twist), mTwist);
2032
2033 VecMat(fn, mTwist, fn);
2034 VecMat(fb, mTwist, fb);
2035 }
2036 }
2037
2038 /*
2039 * Look at the length of the step proceeding from the current
2040 * vertex. If its zero length, just leave the frame of reference
2041 * alone.
2042 */
2043 vv.x = v[0]; vv.y = v[1]; vv.z = v[2];
2044 len = _dxfvector_length_3D(&vv);
2045 if (len == 0.0)
2046 {
2047 n[0] = fn->x; n[1] = fn->y; n[2] = fn->z;
2048 b[0] = fb->x; b[1] = fb->y; b[2] = fb->z;
2049 return OK;
2050 }
2051
2052 /*
2053 * DXDot the incoming and outgoing tangents to determine
2054 * whether a bend occurs.
2055 */
2056 _dxfvector_scale_3D(&vv, 1.0/len, &nt);
2057 cA = _dxfvector_dot_3D(&nt, ft);
2058
2059 /*
2060 * If there's a bend angle...
2061 */
2062 if (cA < 0.999999)
2063 {
2064 float sA, angle;
2065
2066 bend = 1;
2067
2068 /*
2069 * DXNormalize the bend axis
2070 */
2071 _dxfvector_cross_3D(&nt, ft, &cross);
2072 _dxfvector_normalize_3D(&cross, &cross);
2073
2074 /*
2075 * Rotate the incoming frame of reference by half the angle to
2076 * get a vertex FOR
2077 */
2078 angle = acos(cA)/2;
2079 cA = cos(angle);
2080 sA = sin(angle);
2081
2082 RotateAroundVector(&cross, cA, sA, mBend);
2083
2084 VecMat(fn, mBend, fn);
2085 VecMat(fb, mBend, fb);
2086 VecMat(ft, mBend, ft);
2087
2088 /*
2089 * Now we have the vertex frame of reference
2090 */
2091 n[0] = fn->x; n[1] = fn->y; n[2] = fn->z;
2092 b[0] = fb->x; b[1] = fb->y; b[2] = fb->z;
2093
2094 /*
2095 * If there was a bend, perform the second half of the rotation
2096 */
2097 if (bend)
2098 {
2099 VecMat(fn, mBend, fn);
2100 VecMat(fb, mBend, fb);
2101 }
2102
2103 /*
2104 * Outgoing tangent is normalized exiting segment vector. Also
2105 * normalize the frame normal and binormal.
2106 */
2107 _dxfvector_normalize_3D(fn, fn);
2108 _dxfvector_normalize_3D(fb, fb);
2109 ft->x = nt.x; ft->y = nt.y; ft->z = nt.z;
2110 }
2111 else
2112 {
2113 n[0] = fn->x; n[1] = fn->y; n[2] = fn->z;
2114 b[0] = fb->x; b[1] = fb->y; b[2] = fb->z;
2115 }
2116
2117 return OK;
2118 }
2119
2120 /*
2121 * Create a 3x3 matrix defined by the axis v and the sin and cosine
2122 * of an angle Theta.
2123 */
2124 static void
RotateAroundVector(Vector * v,float c,float s,float * M)2125 RotateAroundVector(Vector *v, float c, float s, float *M)
2126 {
2127 float x, y, z;
2128 float sx, sy, sz;
2129
2130 x = v->x; y = v->y; z = v->z;
2131 sx = x*x; sy = y*y; sz = z*z;
2132
2133 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;
2134 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;
2135 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;
2136 }
2137
2138 static int
IsRegular(Object o)2139 IsRegular(Object o)
2140 {
2141 Object c;
2142 Array a, b;
2143
2144 if (DXGetObjectClass(o) == CLASS_FIELD)
2145 {
2146 if (DXEmptyField((Field)o))
2147 return -1;
2148
2149 a = (Array)DXGetComponentValue((Field)o, "connections");
2150 b = (Array)DXGetComponentValue((Field)o, "positions");
2151
2152 if (! a || ! b)
2153 return -1;
2154
2155 if (DXQueryGridConnections(a, NULL, NULL) &&
2156 DXQueryGridPositions(b, NULL, NULL, NULL, NULL))
2157 return 1;
2158 else
2159 return 0;
2160 }
2161 else if (DXGetObjectClass(o) == CLASS_GROUP)
2162 {
2163 int i, j, k;
2164
2165 i = 0; k = -1;
2166 while(NULL != (c = DXGetEnumeratedMember((Group)o, i++, NULL)))
2167 {
2168 j = IsRegular(c);
2169 if (j == 0)
2170 return 0;
2171 if (j == 1)
2172 k = 1;
2173 }
2174
2175 return k;
2176 }
2177 else
2178 {
2179 DXSetError(ERROR_BAD_CLASS, "vector field");
2180 return -1;
2181 }
2182 }
2183
2184 static Error
NewCacheObject(Array pArray,Array tArray,CacheObject * co)2185 NewCacheObject(Array pArray, Array tArray, CacheObject *co)
2186 {
2187 CacheStreak cs = NULL;
2188 float *points, *times;
2189 int nDim, nStreaks, i, j;
2190
2191 *co = NULL;
2192
2193 points = (float *)DXGetArrayData(pArray);
2194 if (! points)
2195 goto error;
2196
2197 times = (float *)DXGetArrayData(tArray);
2198 if (! times)
2199 goto error;
2200
2201 if (! DXGetArrayInfo(pArray, &nStreaks, NULL, NULL, NULL, &nDim))
2202 goto error;
2203
2204 *co = (CacheObject)DXAllocateZero(sizeof(struct cacheObject));
2205 if (! *co)
2206 goto error;
2207
2208 cs = (CacheStreak)DXAllocateZero(nStreaks*sizeof(struct cacheStreak));
2209 if (!cs)
2210 goto error;
2211
2212 (*co)->streaks = cs;
2213 (*co)->nDim = nDim;
2214
2215 for (i = 0; i < nStreaks; i++, cs++)
2216 {
2217 for (j = 0; j < nDim; j++)
2218 cs->point[j] = *points++;
2219 cs->time = *times++;
2220 cs->status = STREAK_PENDING;
2221 cs->start = 0;
2222 cs->streak = DXNewSeries();
2223 if (cs->streak == NULL)
2224 goto error;
2225 }
2226
2227 /*
2228 * Note: no frame data is currently contained in the following
2229 * structure
2230 */
2231 (*co)->frame = -1;
2232 (*co)->nStreaks = nStreaks;
2233 (*co)->frameTimes = NULL;
2234
2235
2236 return OK;
2237
2238 error:
2239 FreeCacheObject(*co);
2240 return ERROR;
2241 }
2242
2243 static Error
FreeCacheObject(Pointer ptr)2244 FreeCacheObject(Pointer ptr)
2245 {
2246 int i;
2247 CacheObject co = ptr;
2248
2249 if (co)
2250 {
2251 if (co->streaks)
2252 {
2253 for (i = 0; i < co->nStreaks; i++)
2254 if (co->streaks[i].streak)
2255 DXDelete((Object)co->streaks[i].streak);
2256 DXFree((Pointer)co->streaks);
2257 }
2258
2259 if (co->frameTimes)
2260 DXFree(co->frameTimes);
2261
2262 if (co->v0)
2263 DXDelete(co->v0);
2264 if (co->c0)
2265 DXDelete((Object)co->c0);
2266 if (co->vf0)
2267 DestroyVectorField(co->vf0);
2268
2269 if (co->v1)
2270 DXDelete(co->v1);
2271 if (co->c1)
2272 DXDelete((Object)co->c1);
2273 if (co->vf1)
2274 DestroyVectorField(co->vf1);
2275
2276 DXFree((Pointer)co);
2277 }
2278
2279 return OK;
2280 }
2281
2282 static void
GetCacheKey(Object * in,Object * keys,int * nk)2283 GetCacheKey(Object *in, Object *keys, int *nk)
2284 {
2285 int n = 0;
2286
2287 if (in[DATA_ARG] && IsSeries(in[DATA_ARG]))
2288 keys[n++] = in[DATA_ARG];
2289
2290 if (in[NAME_ARG]) keys[n++] = in[NAME_ARG];
2291 if (in[POINTS_ARG]) keys[n++] = in[POINTS_ARG];
2292 if (in[TIMES_ARG]) keys[n++] = in[TIMES_ARG];
2293 if (in[HEAD_ARG]) keys[n++] = in[HEAD_ARG];
2294 if (in[FLAG_ARG]) keys[n++] = in[FLAG_ARG];
2295 if (in[SCALE_ARG]) keys[n++] = in[SCALE_ARG];
2296
2297 *nk = n;
2298 }
2299
2300 static Error
GetCacheObject(Object * in,CacheObject * co,Object * obj)2301 GetCacheObject(Object *in, CacheObject *co, Object *obj)
2302 {
2303 int nk;
2304 Object keys[8];
2305 Private pco;
2306 int hold = 0;
2307
2308
2309 if (in[HOLD_ARG])
2310 if (!DXExtractInteger(in[HOLD_ARG], &hold) || hold < 0 || hold > 1)
2311 {
2312 DXSetError(ERROR_BAD_PARAMETER, "hold must be integer flag");
2313 return ERROR;
2314 }
2315
2316 if (hold)
2317 {
2318 int frame;
2319
2320 if (! in[FRAME_ARG])
2321 {
2322 DXSetError(ERROR_BAD_PARAMETER, "if hold is set, frame is required");
2323 return ERROR;
2324 }
2325
2326 if (!DXExtractInteger(in[FRAME_ARG], &frame))
2327 {
2328 DXSetError(ERROR_BAD_PARAMETER, "frame must be integer");
2329 return ERROR;
2330 }
2331
2332 if (frame != 0)
2333 {
2334 pco = (Private)DXGetCacheEntry("Streakline", 0, 1, in[NAME_ARG]);
2335 if (! pco)
2336 {
2337 DXSetError(ERROR_BAD_PARAMETER,
2338 "hold set and frame not 0, but no cached state is present");
2339 return ERROR;
2340 }
2341 }
2342 else
2343 {
2344 DXSetCacheEntry(NULL, CACHE_PERMANENT,
2345 "Streakline", 0, 1, in[NAME_ARG]);
2346 pco = NULL;
2347 }
2348 }
2349 else
2350 {
2351 GetCacheKey(in, keys, &nk);
2352 pco = (Private)DXGetCacheEntryV("Streakline", 0, nk, keys);
2353 }
2354
2355 if (! pco)
2356 {
2357 *co = NULL;
2358 *obj = NULL;
2359 }
2360 else
2361 {
2362
2363 *co = (CacheObject)DXGetPrivateData(pco);
2364 *obj = DXReference((Object)pco);
2365 }
2366
2367 return OK;
2368 }
2369
2370 static Object
SetCacheObject(Object * in,CacheObject co)2371 SetCacheObject(Object *in, CacheObject co)
2372 {
2373 Private pco = NULL;
2374 Object keys[8];
2375 int nk;
2376
2377 pco = DXNewPrivate((Pointer)co, FreeCacheObject);
2378 if (! pco)
2379 return ERROR;
2380
2381 DXReference((Object)pco);
2382
2383 if (in[HOLD_ARG])
2384 {
2385 if (! DXSetCacheEntry((Object)pco, CACHE_PERMANENT,
2386 "Streakline", 0, 1, in[NAME_ARG]))
2387 return ERROR;
2388 }
2389 else
2390 {
2391 GetCacheKey(in, keys, &nk);
2392
2393 if (! DXSetCacheEntryV((Object)pco, CACHE_PERMANENT,
2394 "Streakline", 0, nk, keys))
2395 goto error;
2396 }
2397
2398 return (Object)pco;
2399
2400 error:
2401 DXDelete((Object)pco);
2402 return NULL;
2403 }
2404
2405 static Error
GetTail(Object o,char * name,Pointer l)2406 GetTail(Object o, char *name, Pointer l)
2407 {
2408 Field f;
2409 Object next;
2410 int n, i;
2411 byte *ptr;
2412 int size;
2413 Array a;
2414
2415 if (DXGetObjectClass(o) == CLASS_FIELD)
2416 {
2417 f = (Field)o;
2418 }
2419 else /* object is a series */
2420 {
2421 f = NULL;
2422 for (i = 0; NULL != (next = DXGetSeriesMember((Series)o, i, NULL)); i++)
2423 f = (Field)next;
2424 }
2425
2426 if (! f || DXEmptyField(f))
2427 {
2428 DXSetError(ERROR_INTERNAL, "GetTail called with bad field");
2429 return ERROR;
2430 }
2431
2432 a = (Array)DXGetComponentValue(f, name);
2433 if (! a)
2434 {
2435 DXSetError(ERROR_INTERNAL, "GetTail: no component %s in field", name);
2436 return ERROR;
2437 }
2438
2439 size = DXGetItemSize(a);
2440 DXGetArrayInfo(a, &n, NULL, NULL, NULL, NULL);
2441
2442 ptr = ((byte *)DXGetArrayData(a)) + (n-1)*size;
2443 memcpy(l, ptr, size);
2444
2445 return OK;
2446 }
2447
2448 static Error
SetupStreakVars(StreakVars * svars,CacheObject cstreak)2449 SetupStreakVars(StreakVars *svars, CacheObject cstreak)
2450 {
2451 svars->I0 = svars->I1 = NULL;
2452
2453 svars->vf0 = cstreak->vf0;
2454 svars->vf1 = cstreak->vf1;
2455
2456 svars->t0 = cstreak->t0;
2457 svars->t1 = cstreak->t1;
2458
2459 svars->inside0 = svars->inside1 = 0;
2460
2461 svars->nDim = cstreak->nDim;
2462
2463 return OK;
2464 }
2465
2466 static Error
FreeStreakVars(StreakVars * svars)2467 FreeStreakVars(StreakVars *svars)
2468 {
2469 if (svars->I0)
2470 (*(svars->I0->currentVectorGrp->FreeInstanceVars))(svars->I0);
2471 if (svars->I1)
2472 (*(svars->I1->currentVectorGrp->FreeInstanceVars))(svars->I1);
2473 svars->I0 = svars->I1 = NULL;
2474 return OK;
2475 }
2476
2477 static Error
GetFrameData(Object vectors,Object curls,CacheObject co)2478 GetFrameData(Object vectors, Object curls, CacheObject co)
2479 {
2480 Object curl = NULL;
2481
2482 co->frame ++;
2483
2484 if (co->v0)
2485 DXDelete((Object)co->v0);
2486 co->v0 = co->v1;
2487
2488 if (co->c0)
2489 DXDelete((Object)co->c0);
2490 co->c0 = co->c1;
2491
2492 if (co->vf0)
2493 {
2494 DestroyVectorField(co->vf0);
2495 co->vf0 = NULL;
2496 }
2497
2498 if (co->vf1)
2499 {
2500 co->vf0 = co->vf1;
2501 if (! ResetVectorField(co->vf0))
2502 goto error;
2503 }
2504
2505 co->t0 = co->t1;
2506
2507 if (IsSeries(vectors))
2508 {
2509 co->v1 = DXGetSeriesMember((Series)vectors, co->frame, &co->t1);
2510 if (! co->v1)
2511 {
2512 DXSetError(ERROR_MISSING_DATA, "data: empty series group");
2513 goto error;
2514 }
2515
2516 co->v1 = DXReference(DXCopy(co->v1, COPY_STRUCTURE));
2517 if (! co->v1)
2518 goto error;
2519
2520 if (curls)
2521 {
2522 float t;
2523
2524 curl = DXGetSeriesMember((Series)curls, co->frame, &t);
2525 if (! curl)
2526 {
2527 DXSetError(ERROR_MISSING_DATA, "curl: empty series group");
2528 goto error;
2529 }
2530
2531 if (t != co->t1)
2532 {
2533 DXSetError(ERROR_DATA_INVALID,
2534 "curl and data series positions do not match");
2535 goto error;
2536 }
2537 }
2538 }
2539 else
2540 {
2541 Object attr;
2542
2543 attr = DXGetAttribute(vectors, "series position");
2544 if (! attr)
2545 {
2546 DXSetError(ERROR_MISSING_DATA, "series position");
2547 goto error;
2548 }
2549
2550 if (! DXExtractFloat(attr, &co->t1))
2551 {
2552 DXSetError(ERROR_DATA_INVALID, "series position");
2553 goto error;
2554 }
2555
2556 co->v1 = DXReference(DXCopy(vectors, COPY_STRUCTURE));
2557 if (! co->v1)
2558 goto error;
2559
2560 curl = curls;
2561
2562 }
2563
2564 if (curl)
2565 {
2566 co->c1 = DXNewInterpolator(curl, INTERP_INIT_PARALLEL, 0.0);
2567 if (! co->c1)
2568 goto error;
2569
2570 DXReference((Object)co->c1);
2571 }
2572
2573 co->vf1 = InitVectorField(co->v1);
2574 if (! co->vf1)
2575 goto error;
2576
2577 if (co->frameTimes)
2578 co->frameTimes =
2579 (float *)DXReAllocate(co->frameTimes,
2580 (co->frame+1)*sizeof(float));
2581 else
2582 co->frameTimes =
2583 (float *)DXAllocate((co->frame+1)*sizeof(float));
2584
2585 if (! co->frameTimes)
2586 goto error;
2587
2588 co->frameTimes[co->frame] = co->t1;
2589
2590 return OK;
2591
2592 error:
2593 return ERROR;
2594 }
2595
2596 static InstanceVars
Streak_FindElement_1(VectorField vf,POINT_TYPE * point,VectorGrp last)2597 Streak_FindElement_1(VectorField vf, POINT_TYPE *point, VectorGrp last)
2598 {
2599 InstanceVars I;
2600 int i;
2601
2602 for (i = 0; i < vf->nmembers; i++)
2603 if (vf->members[i] != last)
2604 {
2605 I = (*(vf->members[i]->NewInstanceVars))(vf->members[i]);
2606
2607 if ((*(I->currentVectorGrp->FindElement))(I, point))
2608 {
2609 vf->current = vf->members[i];
2610 return I;
2611 }
2612
2613 (*(I->currentVectorGrp->FreeInstanceVars))(I);
2614 }
2615
2616 return 0;
2617 }
2618
2619 static InstanceVars
Streak_FindMultiGridContinuation_1(VectorField vf,POINT_TYPE * point,VectorGrp last)2620 Streak_FindMultiGridContinuation_1(VectorField vf, POINT_TYPE *point, VectorGrp last)
2621 {
2622 InstanceVars I;
2623 int i;
2624
2625 for (i = 0; i < vf->nmembers; i++)
2626 if (vf->members[i] != last)
2627 {
2628 I = (*(vf->members[i]->NewInstanceVars))(vf->members[i]);
2629
2630 if ((*(I->currentVectorGrp->FindMultiGridContinuation))(I, point))
2631 {
2632 vf->current = vf->members[i];
2633 return I;
2634 }
2635
2636 (*(I->currentVectorGrp->FreeInstanceVars))(I);
2637 }
2638
2639 return 0;
2640 }
2641
2642 static int
Streak_FindElement(StreakVars * svars,POINT_TYPE * p)2643 Streak_FindElement(StreakVars *svars, POINT_TYPE *p)
2644 {
2645 if (! svars->inside0)
2646 {
2647 VectorGrp last;
2648
2649 if (svars->I0)
2650 {
2651 last = svars->I0->currentVectorGrp;
2652 (*(last->FreeInstanceVars))(svars->I0);
2653 }
2654 else
2655 last = NULL;
2656
2657 svars->I0 = Streak_FindElement_1(svars->vf0, p, last);
2658 if (! svars->I0)
2659 return 0;
2660 }
2661
2662 if (! svars->inside1)
2663 {
2664 VectorGrp last;
2665
2666 if (svars->I1)
2667 {
2668 last = svars->I1->currentVectorGrp;
2669 (*(last->FreeInstanceVars))(svars->I1);
2670 }
2671 else
2672 last = NULL;
2673
2674 svars->I1 = Streak_FindElement_1(svars->vf1, p, last);
2675 if (! svars->I1)
2676 return 0;
2677 }
2678 return 1;
2679 }
2680
2681
2682 static int
Streak_FindMultiGridContinuation(StreakVars * svars,POINT_TYPE * p)2683 Streak_FindMultiGridContinuation(StreakVars *svars, POINT_TYPE *p)
2684 {
2685 if (! svars->inside0)
2686 {
2687 VectorGrp last;
2688
2689 if (svars->I0)
2690 {
2691 last = svars->I0->currentVectorGrp;
2692 (*(last->FreeInstanceVars))(svars->I0);
2693 }
2694 else
2695 last = NULL;
2696
2697 svars->I0 = Streak_FindMultiGridContinuation_1(svars->vf0, p, last);
2698 if (! svars->I0)
2699 return 0;
2700 }
2701
2702 if (! svars->inside1)
2703 {
2704 VectorGrp last;
2705
2706 if (svars->I1)
2707 {
2708 last = svars->I1->currentVectorGrp;
2709 (*(last->FreeInstanceVars))(svars->I1);
2710 }
2711 else
2712 last = NULL;
2713
2714 svars->I1 = Streak_FindMultiGridContinuation_1(svars->vf1, p, last);
2715 if (! svars->I1)
2716 return 0;
2717 }
2718 return 1;
2719 }
2720
2721
2722 static Error
Streak_Interpolate(StreakVars * svars,POINT_TYPE * p,double t,VECTOR_TYPE * v)2723 Streak_Interpolate(StreakVars *svars, POINT_TYPE *p, double t, VECTOR_TYPE *v)
2724 {
2725 double d;
2726 int i;
2727 VECTOR_TYPE v0[3], v1[3];
2728
2729 if (! (*(svars->I0->currentVectorGrp->Interpolate))(svars->I0, p, v0))
2730 return ERROR;
2731
2732 if (! (*(svars->I1->currentVectorGrp->Interpolate))(svars->I1, p, v1))
2733 return ERROR;
2734
2735 d = (t - svars->t0) / (svars->t1 - svars->t0);
2736
2737 for (i = 0; i < svars->nDim; i++)
2738 v[i] = v0[i] + d*(v1[i] - v0[i]);
2739
2740 return OK;
2741 }
2742
2743 static Error
Streak_StepTime(StreakVars * svars,double c,VECTOR_TYPE * v,double * t)2744 Streak_StepTime(StreakVars *svars, double c, VECTOR_TYPE *v, double *t)
2745 {
2746 double t0, t1;
2747
2748 if (! (*(svars->I0->currentVectorGrp->StepTime))(svars->I0, c, v, &t0))
2749 return ERROR;
2750
2751 if (! (*(svars->I1->currentVectorGrp->StepTime))(svars->I1, c, v, &t1))
2752 return ERROR;
2753
2754 if (t0 < t1)
2755 *t = t0;
2756 else
2757 *t = t1;
2758
2759 return OK;
2760 }
2761
2762 static int
Streak_InsideWeights(StreakVars * svars,POINT_TYPE * p)2763 Streak_InsideWeights(StreakVars *svars, POINT_TYPE *p)
2764 {
2765
2766 svars->inside0 = (*(svars->I0->currentVectorGrp->Weights))(svars->I0, p);
2767 svars->inside1 = (*(svars->I1->currentVectorGrp->Weights))(svars->I1, p);
2768
2769 if (svars->inside0 && svars->inside1)
2770 return 1;
2771 else
2772 return 0;
2773 }
2774
2775 #if 0
2776 static int
2777 Streak_Weights(StreakVars *svars, POINT_TYPE *p)
2778 {
2779
2780 (*(svars->I0->currentVectorGrp->Weights))(svars->I0, p);
2781 (*(svars->I1->currentVectorGrp->Weights))(svars->I1, p);
2782
2783 return 1;
2784 }
2785 #endif
2786
2787 static int
Streak_FaceWeights(StreakVars * svars,POINT_TYPE * p)2788 Streak_FaceWeights(StreakVars *svars, POINT_TYPE *p)
2789 {
2790
2791 if (! svars->inside0)
2792 (*(svars->I0->currentVectorGrp->FaceWeights))(svars->I0, p);
2793
2794 if (! svars->inside1)
2795 (*(svars->I0->currentVectorGrp->FaceWeights))(svars->I1, p);
2796
2797 return 1;
2798 }
2799
2800 static Error
Streak_FindBoundary(StreakVars * svars,POINT_TYPE * p,VECTOR_TYPE * v,double * t)2801 Streak_FindBoundary(StreakVars *svars, POINT_TYPE *p, VECTOR_TYPE *v, double *t)
2802 {
2803 double t0, t1;
2804
2805 if (svars->inside0)
2806 t0 = DXD_MAX_FLOAT;
2807 else
2808 if (! (*(svars->I0->currentVectorGrp->FindBoundary))
2809 (svars->I0, p, v, &t0))
2810 return ERROR;
2811
2812 if (svars->inside1)
2813 t1 = DXD_MAX_FLOAT;
2814 else
2815 if (! (*(svars->I1->currentVectorGrp->FindBoundary))
2816 (svars->I1, p, v, &t1))
2817 return ERROR;
2818
2819 if (t1 < t0)
2820 *t = t1;
2821 else
2822 *t = t0;
2823
2824 if (*t == DXD_MAX_FLOAT)
2825 {
2826 DXSetError(ERROR_INTERNAL, "error in Streak_FindBoundary");
2827 return ERROR;
2828 }
2829
2830 return OK;
2831 }
2832
2833 static int
Streak_Neighbor(StreakVars * svars,VECTOR_TYPE * v)2834 Streak_Neighbor(StreakVars *svars, VECTOR_TYPE *v)
2835 {
2836 int found;
2837
2838 if (! svars->inside0)
2839 {
2840 found = (*(svars->I0->currentVectorGrp->Neighbor))(svars->I0, v);
2841 if (found != 1)
2842 return found;
2843 }
2844
2845 if (! svars->inside1)
2846 {
2847 found = (*(svars->I1->currentVectorGrp->Neighbor))(svars->I1, v);
2848 if (found != 1)
2849 return found;
2850 }
2851
2852 return 1;
2853 }
2854
2855 static int
IsSeries(Object o)2856 IsSeries(Object o)
2857 {
2858 if (DXGetObjectClass(o) != CLASS_GROUP)
2859 return 0;
2860 if (DXGetGroupClass((Group)o) != CLASS_SERIES)
2861 return 0;
2862 return 1;
2863 }
2864
2865 static Error
ResetVectorField(VectorField vf)2866 ResetVectorField(VectorField vf)
2867 {
2868 int i;
2869
2870 for (i = 0; i < vf->nmembers; i++)
2871 {
2872 VectorGrp vg = vf->members[i];
2873
2874 if (vg->Reset)
2875 if (! (*(vg->Reset))(vg))
2876 return ERROR;
2877 }
2878
2879 return OK;
2880 }
2881
2882 static Error
GeometryCheck(Object vectors,int nD)2883 GeometryCheck(Object vectors, int nD)
2884 {
2885 Class class = DXGetObjectClass(vectors);
2886 if (class == CLASS_GROUP)
2887 class = DXGetGroupClass((Group)vectors);
2888
2889 switch(class)
2890 {
2891 case CLASS_FIELD:
2892 {
2893 Field f;
2894 Array a;
2895 Object o;
2896 char *str;
2897 int n;
2898
2899 f = (Field)vectors;
2900 if (DXEmptyField(f))
2901 return OK;
2902
2903 a = (Array)DXGetComponentValue(f, "positions");
2904 if (! DXGetArrayInfo(a, NULL, NULL, NULL, NULL, &n))
2905 goto error;
2906
2907 if (n != nD)
2908 {
2909 DXSetError(ERROR_DATA_INVALID,
2910 "vectors are %d-D, vector space is %d-D", nD, n);
2911 goto error;
2912 }
2913
2914 if (! DXGetComponentValue(f, "connections"))
2915 {
2916 DXSetError(ERROR_MISSING_DATA, "connections");
2917 goto error;
2918 }
2919
2920 o = DXGetComponentAttribute(f, "connections", "element type");
2921 if (! o)
2922 {
2923 DXSetError(ERROR_MISSING_DATA, "connections element type");
2924 goto error;
2925 }
2926
2927 if (DXGetObjectClass(o) != CLASS_STRING)
2928 {
2929 DXSetError(ERROR_BAD_CLASS, "connections element type attribute");
2930 goto error;
2931 }
2932
2933 str = DXGetString((String)o);
2934
2935 switch(n)
2936 {
2937 case 2:
2938 if (strcmp(str, "quads") && strcmp(str, "triangles"))
2939 {
2940 DXSetError(ERROR_DATA_INVALID,
2941 "%d-D vector space requires %d-D elements", n, n);
2942 goto error;
2943 }
2944 break;
2945
2946 case 3:
2947 if (strcmp(str, "cubes") && strcmp(str, "tetrahedra"))
2948 {
2949 DXSetError(ERROR_DATA_INVALID,
2950 "%d-D vector space requires %d-D elements", n, n);
2951 goto error;
2952 }
2953 break;
2954
2955 default:
2956 DXSetError(ERROR_DATA_INVALID, "2- or 3-D vectors required");
2957 goto error;
2958 }
2959 break;
2960 }
2961
2962 case CLASS_SERIES:
2963 case CLASS_MULTIGRID:
2964 case CLASS_COMPOSITEFIELD:
2965 {
2966 Object child;
2967 int i;
2968
2969 i = 0;
2970 while (NULL != (child = DXGetEnumeratedMember((Group)vectors, i++, NULL)))
2971 if (! GeometryCheck(child, nD))
2972 return ERROR;
2973
2974 break;
2975 }
2976
2977 default:
2978 DXSetError(ERROR_DATA_INVALID,
2979 "vectors must be field or composite field");
2980 goto error;
2981 }
2982
2983 return OK;
2984
2985 error:
2986 return ERROR;
2987 }
2988