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/_irregstream.c,v 1.9 2006/01/03 17:02:21 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 #include "stream.h"
21 #include "vectors.h"
22 #include "_divcurl.h"
23 #include "_partnbrs.h"
24
25 typedef struct neighborsHandle
26 {
27 Array array;
28 int *neighbors;
29 int nPerE, nDim, wrap;
30 int gcounts[10];
31 int nstrides[10];
32 } *NeighborsHandle;
33
34 /*
35 ***********************************
36 * Irreg interface
37 ***********************************
38 */
39
40 typedef struct irreg_InstanceVars
41 {
42 struct instanceVars i;
43 int cp; /* current partition index */
44 int ce; /* current element index */
45 int ct; /* current primitive */
46 int flags; /* flags for zero-area faces */
47 float w[4]; /* current weights */
48 float planes[16]; /* four face planar equations */
49 int face; /* intersection face */
50 float edges[18]; /* primitive edge vectors */
51 ArrayHandle pHandle;
52 ArrayHandle cHandle;
53 NeighborsHandle nHandle;
54 } *Irreg_InstanceVars;
55
56 /*
57 * The following table contains an entry for each face of each of the six
58 * tetras within the element. If the face is external to the element, nelt is
59 * the index of the element neighbor pointer for that element face, -1
60 * otherwise. nprim is the index of the containing tetra within either
61 * the neighbor element (if nelt != -1) or the current tetra (if nelt == -1).
62 */
63 typedef struct {
64 int nelt;
65 int nprim;
66 } NbrTable;
67
68 typedef struct irreg_VectorGrp *Irreg_VectorGrp;
69 typedef struct irreg_VectorPart *Irreg_VectorPart;
70
71 /*
72 * defines for flag field
73 * PARITY_KNOWN the parity of the element list is known
74 * CURRENT_DEGENERATE the current element is degenerate
75 *
76 * the low 16 bits contain flags denoting whether the corresponding
77 * face's planar equation must be reversed.
78 */
79 #define PARITY_KNOWN 0x0100
80 #define CURRENT_DEGENERATE 0x0200
81
82 #define IsParityKnown(i) \
83 (((Irreg_InstanceVars)(i))->flags & PARITY_KNOWN)
84 #define SetParityKnown(i) \
85 (((Irreg_InstanceVars)(i))->flags |= PARITY_KNOWN)
86
87 #define IsCurrentDegenerate(i) \
88 (((Irreg_InstanceVars)(i))->flags & CURRENT_DEGENERATE)
89 #define SetCurrentDegenerate(i) \
90 (((Irreg_InstanceVars)(i))->flags |= CURRENT_DEGENERATE)
91 #define ResetCurrentDegenerate(i) \
92 (((Irreg_InstanceVars)(i))->flags &= ~CURRENT_DEGENERATE)
93
94 #define IsFaceNReversed(i,n) \
95 (((Irreg_InstanceVars)(i))->flags & (1 << n))
96 #define SetFaceNReversed(i,n) \
97 (((Irreg_InstanceVars)(i))->flags |= (1 << n))
98
99 struct irreg_VectorGrp
100 {
101 struct vectorGrp P;
102 int *primTable;
103 NbrTable *nbrTable;
104 int (*InitElement)(Irreg_InstanceVars);
105 int (*FaceWeights)(Irreg_InstanceVars, POINT_TYPE *);
106 int primsPerElement;
107 int nbrsPerElement;
108 int vrtsPerElement;
109 int nbrsPerPrim;
110 int nDimensions;
111 };
112
113 struct irreg_VectorPart
114 {
115 struct vectorPart p;
116 int nElements;
117 int nPositions;
118 int *elements;
119 int *neighbors;
120 int *partitionNeighbors;
121 float *points;
122 int constantData;
123 float *vectors;
124 byte *visited;
125 byte gridConnections;
126 int nd;
127 float origin[3];
128 int pcounts[3];
129 float deltas[9];
130 int gcounts[3];
131 int gstrides[3];
132 int gbumps[8];
133 int nstrides[3];
134 int pstrides[3];
135 int wrapFlags;
136 InvalidComponentHandle invElements;
137 Array pArray, cArray, nArray;
138 };
139
140 #define SetElementVisited(ip, i) \
141 { \
142 if (i < 0 || i > (ip)->nElements) \
143 DXSetError(ERROR_INTERNAL, "invalid visited"); \
144 (ip)->visited[i] = 1; \
145 }
146
147 #define IsElementVisited(ip, i) ((ip)->visited[i] == 1)
148
149 /*
150 ***********************************
151 * Tetrahedra definition
152 ***********************************
153 */
154
155 static int tetrahedra_primTable[] =
156 {
157 0, 1, 2, 3
158 };
159
160 static NbrTable tetrahedra_nbrTable[] =
161 {
162 {0, 0}, {1, 0}, {2, 0}, {3, 0}
163 };
164
165 #define TETRAHEDRA_N_PRIMS 1
166 #define TETRAHEDRA_N_NBRS 4
167 #define TETRAHEDRA_N_VRTS 4
168 #define TETRAHEDRA_N_DIMS 3
169 #define TETRAHEDRA_N_PRIM_NBRS 4
170
171 /*
172 ***********************************
173 * Triangle definition
174 ***********************************
175 */
176
177 static int triangles_primTable[] =
178 {
179 0, 1, 2
180 };
181
182 static NbrTable triangles_nbrTable[] =
183 {
184 {0, 0}, {1, 0}, {2, 0}
185 };
186
187 #define TRIANGLES_N_PRIMS 1
188 #define TRIANGLES_N_NBRS 3
189 #define TRIANGLES_N_VRTS 3
190 #define TRIANGLES_N_DIMS 2
191 #define TRIANGLES_N_PRIM_NBRS 3
192
193 /*
194 ***********************************
195 * Irregular cubes definition
196 ***********************************
197 */
198
199 static int icubes_primTable[] =
200 {
201 0, 3, 5, 1,
202 5, 2, 7, 3,
203 5, 4, 7, 2,
204 5, 3, 0, 2,
205 5, 2, 0, 4,
206 4, 2, 6, 7
207 };
208
209 #define ICUBES_N_PRIMS 6
210 #define ICUBES_N_NBRS 6
211 #define ICUBES_N_VRTS 8
212 #define ICUBES_N_DIMS 3
213 #define ICUBES_N_PRIM_NBRS 4
214
215 static NbrTable icubes_nbrTable[] =
216 {
217 { 5, 4}, { 2, 1}, { 0, 2}, {-1, 3},
218 { 3, 0}, { 5, 5}, {-1, 3}, {-1, 2},
219 {-1, 5}, {-1, 1}, {-1, 4}, { 1, 0},
220 { 0, 5}, {-1, 4}, {-1, 1}, {-1, 0},
221 { 4, 0}, { 2, 5}, {-1, 2}, {-1, 3},
222 { 3, 4}, { 1, 3}, {-1, 2}, { 4, 1}
223 };
224
225 /*
226 ***********************************
227 * Irregular quads definition
228 ***********************************
229 */
230
231 static int iquads_primTable[] =
232 {
233 0, 3, 1,
234 0, 2, 3
235 };
236
237 #define IQUADS_N_PRIMS 2
238 #define IQUADS_N_NBRS 4
239 #define IQUADS_N_VRTS 4
240 #define IQUADS_N_DIMS 2
241 #define IQUADS_N_PRIM_NBRS 3
242
243 static NbrTable iquads_nbrTable[] =
244 {
245 { 3, 1}, { 0, 1}, {-1, 1},
246 { 1, 0}, {-1, 0}, { 2, 0}
247 };
248
249 /*
250 ***********************************
251 * Element definition table
252 ***********************************
253 */
254 typedef struct
255 {
256 char *elementType;
257 int *primTable;
258 NbrTable *nbrTable;
259 int primsPerElement;
260 int nbrsPerElement;
261 int vrtsPerElement;
262 int nDimensions;
263 int nbrsPerPrim;
264 } ElementDef;
265
266 static ElementDef elementTable[] =
267 {
268 {
269 "tetrahedra",
270 tetrahedra_primTable,
271 tetrahedra_nbrTable,
272 TETRAHEDRA_N_PRIMS,
273 TETRAHEDRA_N_NBRS,
274 TETRAHEDRA_N_VRTS,
275 TETRAHEDRA_N_DIMS,
276 TETRAHEDRA_N_PRIM_NBRS
277 },
278 {
279 "triangles",
280 triangles_primTable,
281 triangles_nbrTable,
282 TRIANGLES_N_PRIMS,
283 TRIANGLES_N_NBRS,
284 TRIANGLES_N_VRTS,
285 TRIANGLES_N_DIMS,
286 TRIANGLES_N_PRIM_NBRS
287 },
288 {
289 "cubes",
290 icubes_primTable,
291 icubes_nbrTable,
292 ICUBES_N_PRIMS,
293 ICUBES_N_NBRS,
294 ICUBES_N_VRTS,
295 ICUBES_N_DIMS,
296 ICUBES_N_PRIM_NBRS
297 },
298 {
299 "quads",
300 iquads_primTable,
301 iquads_nbrTable,
302 IQUADS_N_PRIMS,
303 IQUADS_N_NBRS,
304 IQUADS_N_VRTS,
305 IQUADS_N_DIMS,
306 IQUADS_N_PRIM_NBRS
307 },
308 {
309 NULL,
310 NULL,
311 NULL,
312 -1,
313 -1,
314 -1,
315 -1
316 }
317 };
318
319 static VectorPart Irreg_InitVectorPart(Field, Irreg_VectorGrp);
320 static Error Irreg_DestroyVectorPart(VectorPart);
321 static int Irreg_FindElement_VectorPart(Irreg_InstanceVars,
322 int, POINT_TYPE *);
323 static int Irreg_FindMultiGridContinuation_VectorPart(
324 Irreg_InstanceVars, int, POINT_TYPE *);
325 static InstanceVars Irreg_NewInstanceVars(VectorGrp);
326 static Error Irreg_FreeInstanceVars(InstanceVars);
327 static int Irreg_FindElement(InstanceVars, POINT_TYPE *);
328 static int Irreg_FindMultiGridContinuation(InstanceVars, POINT_TYPE *);
329 static int Irreg_Interpolate(InstanceVars, POINT_TYPE *,
330 VECTOR_TYPE *);
331 static Error Irreg_StepTime(InstanceVars, double,
332 VECTOR_TYPE *, double *);
333 static Error Irreg_Walk(InstanceVars, POINT_TYPE *,
334 VECTOR_TYPE *, POINT_TYPE *);
335 static Error Irreg_FindBoundary(InstanceVars, POINT_TYPE *,
336 VECTOR_TYPE *, double *);
337 static int Irreg_Neighbor(InstanceVars, VECTOR_TYPE *);
338 static Error Irreg_Delete(VectorGrp);
339 static Error Irreg_CurlMap(VectorGrp, MultiGrid);
340 static Error Irreg_CurlMap_Task(Pointer);
341 static Error Irreg_CurlMap_Field(Irreg_VectorGrp,
342 Irreg_VectorPart, Group, Field *);
343 static Error Irreg_ResetVectorGrp(VectorGrp);
344 static Error Irreg_ResetVectorPart(VectorPart);
345 static int Irreg_Ghost(InstanceVars I, POINT_TYPE *p);
346
347 static int Tetras_Weights(InstanceVars, POINT_TYPE *);
348 static int Tetras_FaceWeights(InstanceVars, POINT_TYPE *);
349 static int Tetras_InitElement(Irreg_InstanceVars);
350 static int Tris_Weights(InstanceVars, POINT_TYPE *);
351 static int Tris_FaceWeights(InstanceVars, POINT_TYPE *);
352 static int Tris_InitElement(Irreg_InstanceVars);
353
354 static NeighborsHandle
DXCreateNeighborsHandle(Field field)355 DXCreateNeighborsHandle(Field field)
356 {
357 Array cA;
358 NeighborsHandle handle=NULL;
359
360 cA = (Array)DXGetComponentValue(field, "connections");
361 if (! cA)
362 {
363 DXSetError(ERROR_MISSING_DATA, "no connections component");
364 goto error;
365 }
366
367 handle = (NeighborsHandle)DXAllocateZero(sizeof(struct neighborsHandle));
368 if (! handle)
369 goto error;
370
371 handle->wrap = 0;
372 if (DXGetAttribute((Object)cA, "wrap I")) handle->wrap |= 0x1;
373 if (DXGetAttribute((Object)cA, "wrap J")) handle->wrap |= 0x2;
374 if (DXGetAttribute((Object)cA, "wrap K")) handle->wrap |= 0x4;
375
376 if (DXQueryGridConnections(cA, &handle->nDim, handle->gcounts))
377 {
378 int i;
379
380 for (i = 0; i < handle->nDim; i++)
381 handle->gcounts[i] -= 1;
382
383 handle->nstrides[handle->nDim - 1] = 1;
384 for (i = handle->nDim-2; i >= 0; i--)
385 handle->nstrides[i] = handle->nstrides[i+1] * handle->gcounts[i+1];
386 }
387 else
388 {
389 handle->array = DXNeighbors(field);
390 if (! handle->array)
391 goto error;
392
393 DXReference((Object)(handle->array));
394
395 DXGetArrayInfo(handle->array, NULL, NULL, NULL, NULL, &handle->nPerE);
396
397 handle->neighbors = (int *)DXGetArrayData(handle->array);
398 }
399
400 return handle;
401
402 error:
403 if (handle->array)
404 DXDelete((Object)handle->array);
405 DXFree((Pointer)handle);
406 return NULL;
407 }
408
409 static void
DXFreeNeighborsHandle(NeighborsHandle handle)410 DXFreeNeighborsHandle(NeighborsHandle handle)
411 {
412 if (handle)
413 {
414 DXDelete((Object)handle->array);
415 DXFree((Pointer)handle);
416 }
417 }
418
419 #define GET_NEIGHBORS(handle, index, buf) \
420 (handle->neighbors ? \
421 (handle->neighbors + handle->nPerE*(index)) : \
422 _dxfGetNeighbors(handle, index, buf))
423
424 static Error
_dxfGetIndices(int offset,int nDim,int * strides,int * indices)425 _dxfGetIndices(int offset, int nDim, int *strides, int *indices)
426 {
427 int i;
428
429 for (i = 0; i < nDim-2; i++)
430 {
431 indices[i] = offset / strides[i];
432 offset = offset % strides[i];
433 }
434
435 indices[nDim-2] = offset / strides[nDim-2];
436 indices[nDim-1] = offset % strides[nDim-2];
437
438 return OK;
439 }
440
441 static int
_dxfGetOffset(int nDim,int * strides,int * indices)442 _dxfGetOffset(int nDim, int *strides, int *indices)
443 {
444 int i, j;
445
446 for (i = j = 0; i < nDim; i++)
447 j += indices[i]*strides[i];
448
449 return j;
450 }
451
452 int *
_dxfGetNeighbors(NeighborsHandle handle,int index,int * buf)453 _dxfGetNeighbors(NeighborsHandle handle, int index, int *buf)
454 {
455 int i, indices[10];
456
457 _dxfGetIndices(index, handle->nDim, handle->nstrides, indices);
458
459 for (i = 0; i < handle->nDim; i++)
460 {
461 if (indices[i] == 0)
462 if (handle->wrap & (1 << i))
463 buf[i<<1] = index +
464 (handle->gcounts[i]-1)*handle->nstrides[i];
465 else
466 buf[i<<1] = -1;
467 else
468 buf[i<<1] = index - handle->nstrides[i];
469
470 if (indices[i] == (handle->gcounts[i]-1))
471 if (handle->wrap & (1 << i))
472 buf[(i<<1)+1] = index -
473 (handle->gcounts[i]-1)*handle->nstrides[i];
474 else
475 buf[(i<<1)+1] = -1;
476 else
477 buf[(i<<1)+1] = index + handle->nstrides[i];
478 }
479
480 return buf;
481 }
482
483 static InstanceVars
Irreg_NewInstanceVars(VectorGrp p)484 Irreg_NewInstanceVars(VectorGrp p)
485 {
486 Irreg_InstanceVars iI;
487
488 iI = (Irreg_InstanceVars)DXAllocate(sizeof(struct irreg_InstanceVars));
489 if (! iI)
490 return NULL;
491
492 memset(iI, -1, sizeof(struct irreg_InstanceVars));
493 iI->flags = 0;
494 iI->pHandle = NULL;
495 iI->cHandle = NULL;
496 iI->nHandle = NULL;
497
498 iI->i.currentVectorGrp = p;
499 iI->i.currentPartition = NULL;
500 iI->i.isRegular = 0;
501
502 return (InstanceVars)iI;
503 }
504
505 static Error
Irreg_FreeInstanceVars(InstanceVars I)506 Irreg_FreeInstanceVars(InstanceVars I)
507 {
508 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
509
510 if (iI)
511 {
512 if (I->currentPartition)
513 DXDelete(I->currentPartition);
514
515 if (iI->pHandle)
516 {
517 DXFreeArrayHandle(iI->pHandle);
518 iI->pHandle = NULL;
519 }
520 if (iI->cHandle)
521 {
522 DXFreeArrayHandle(iI->cHandle);
523 iI->pHandle = NULL;
524 }
525 if (iI->nHandle)
526 {
527 DXFreeNeighborsHandle(iI->nHandle);
528 iI->nHandle = NULL;
529 }
530
531 DXFree((Pointer)I);
532 }
533 return OK;
534 }
535
536 VectorGrp
_dxfIrreg_InitVectorGrp(Object object,char * elementType)537 _dxfIrreg_InitVectorGrp(Object object, char *elementType)
538 {
539 Irreg_VectorGrp P = NULL;
540 ElementDef *eDef;
541 int i, nP;
542
543 /*
544 * Look in the table to see if we have the info necessary
545 * to operate on this element type
546 */
547 for (eDef = elementTable; eDef->elementType != NULL; eDef ++)
548 if (! strcmp(eDef->elementType, elementType))
549 break;
550
551 if (eDef->elementType == NULL)
552 {
553 DXSetError(ERROR_DATA_INVALID,
554 "unknown element type: %s", elementType);
555 return NULL;
556 }
557
558 if (! _dxfPartitionNeighbors(object))
559 return NULL;
560
561 if (DXGetObjectClass(object) == CLASS_FIELD)
562 nP = 1;
563 else
564 for (nP = 0; DXGetEnumeratedMember((Group)object, nP, NULL); nP++);
565
566 P = (Irreg_VectorGrp)DXAllocate(sizeof(struct irreg_VectorGrp));
567 if (! P)
568 goto error;
569
570 P->P.n = nP;
571
572 P->P.p = (VectorPart *)DXAllocateZero(nP * sizeof(struct irreg_VectorPart));
573 if (! P->P.p)
574 goto error;
575
576 /*
577 * Attach generic irregular element type methods to
578 * return object
579 */
580 P->P.NewInstanceVars = Irreg_NewInstanceVars;
581 P->P.FreeInstanceVars = Irreg_FreeInstanceVars;
582 P->P.FindElement = Irreg_FindElement;
583 P->P.FindMultiGridContinuation = Irreg_FindMultiGridContinuation;
584 P->P.FindElement = Irreg_FindElement;
585 P->P.Interpolate = Irreg_Interpolate;
586 P->P.StepTime = Irreg_StepTime;
587 P->P.FindBoundary = Irreg_FindBoundary;
588 P->P.Neighbor = Irreg_Neighbor;
589 P->P.CurlMap = Irreg_CurlMap;
590 P->P.Delete = Irreg_Delete;
591 P->P.Reset = Irreg_ResetVectorGrp;
592 P->P.Walk = Irreg_Walk;
593 P->P.Ghost = Irreg_Ghost;
594
595 /*
596 * Attach element type-dependent info to return object
597 */
598 P->primTable = eDef->primTable;
599 P->nbrTable = eDef->nbrTable;
600 P->primsPerElement = eDef->primsPerElement;
601 P->nbrsPerElement = eDef->nbrsPerElement;
602 P->vrtsPerElement = eDef->vrtsPerElement;
603 P->nbrsPerPrim = eDef->nbrsPerPrim;
604 P->nDimensions = eDef->nDimensions;
605
606 if ((P->P.nDim = eDef->nDimensions) == 3)
607 {
608 P->InitElement = Tetras_InitElement;
609 P->P.Weights = Tetras_Weights;
610 P->P.FaceWeights = Tetras_FaceWeights;
611 }
612 else
613 {
614 P->InitElement = Tris_InitElement;
615 P->P.Weights = Tris_Weights;
616 P->P.FaceWeights = Tris_FaceWeights;
617 }
618
619 if (DXGetObjectClass(object) == CLASS_FIELD)
620 {
621 P->P.p[0] = Irreg_InitVectorPart((Field)object, P);
622 if (! P->P.p[0])
623 goto error;
624
625 P->P.p[0]->field = (Field)object;
626 }
627 else
628 {
629 Group g = (Group)object;
630 Field f;
631
632 i = 0;
633 while (NULL != (f = (Field)DXGetEnumeratedMember(g, i, NULL)))
634 {
635 P->P.p[i] = Irreg_InitVectorPart(f, P);
636 if (! P->P.p[i] && DXGetError() != ERROR_NONE)
637 goto error;
638
639 if (P->P.p[i])
640 P->P.p[i]->field = f;
641
642 i++;
643 }
644 }
645
646 return (VectorGrp)P;
647
648 error:
649 Irreg_Delete((VectorGrp)P);
650 return NULL;
651 }
652
653 static VectorPart
Irreg_InitVectorPart(Field f,Irreg_VectorGrp iP)654 Irreg_InitVectorPart(Field f, Irreg_VectorGrp iP)
655 {
656 Irreg_VectorPart ip = NULL;
657 Array array;
658 int i;
659
660 if (DXEmptyField(f))
661 return NULL;
662
663 ip = (Irreg_VectorPart)DXAllocate(sizeof(struct irreg_VectorPart));
664 if (! ip)
665 goto error;
666
667 if (! _dxfInitVectorPart((VectorPart)ip, f))
668 goto error;
669
670 ip->visited = NULL;
671
672 ip->cArray = (Array)DXGetComponentValue(f, "connections");
673 if (! ip->cArray)
674 goto error;
675
676 if (DXQueryGridConnections(ip->cArray, &ip->nd, ip->gcounts))
677 {
678 ip->gridConnections = 1;
679
680 ip->gstrides[ip->nd-1] = 1;
681 ip->gcounts[ip->nd-1] -= 1;
682 for (i = ip->nd-2; i >= 0; i--)
683 {
684 ip->gcounts[i] -= 1;
685 ip->gstrides[i] = ip->gstrides[i+1]*ip->gcounts[i+1];
686 }
687 }
688 else
689 ip->gridConnections = 0;
690
691 if (! DXTypeCheck(ip->cArray, TYPE_INT,
692 CATEGORY_REAL, 1, iP->vrtsPerElement))
693 {
694 DXSetError(ERROR_DATA_INVALID, "connections");
695 return NULL;
696 }
697
698 ip->wrapFlags = 0;
699 if (DXGetAttribute((Object)ip->cArray, "wrap I")) ip->wrapFlags |= 0x1;
700 if (DXGetAttribute((Object)ip->cArray, "wrap J")) ip->wrapFlags |= 0x2;
701 if (DXGetAttribute((Object)ip->cArray, "wrap K")) ip->wrapFlags |= 0x4;
702
703 DXGetArrayInfo(ip->cArray, &ip->nElements, NULL, NULL, NULL, NULL);
704
705 if (! DXInvalidateConnections((Object)f))
706 goto error;
707
708 array = (Array)DXGetComponentValue(f, "invalid connections");
709 if (! array)
710 {
711 ip->invElements = NULL;
712 }
713 else
714 {
715 ip->invElements = DXCreateInvalidComponentHandle((Object)f, NULL,
716 "connections");
717
718 if (! ip->invElements)
719 goto error;
720 }
721
722 /*
723 * Also use an invalid component handle to keep track of the
724 * elements the streamline visits. This one will be initialized to
725 * all invalid, so we ignore the incoming.
726 */
727 ip->visited = (byte *)DXAllocateZero(ip->nElements*sizeof(byte));
728 if (! ip->visited)
729 goto error;
730
731 ip->pArray = (Array)DXGetComponentValue(f, "positions");
732 if (! ip->pArray)
733 goto error;
734
735 DXGetArrayInfo(ip->pArray, &ip->nPositions, NULL, NULL, NULL, NULL);
736
737 if (! DXTypeCheck(ip->pArray, TYPE_FLOAT, CATEGORY_REAL, 1, iP->P.nDim))
738 {
739 DXSetError(ERROR_DATA_INVALID, "positions");
740 return NULL;
741 }
742
743 ip->nArray = DXNeighbors(f);
744
745 array = (Array)DXGetComponentValue(f, "partition neighbors");
746 if (array)
747 ip->partitionNeighbors = (int *)DXGetArrayData(array);
748 else
749 ip->partitionNeighbors = NULL;
750
751 array = (Array)DXGetComponentValue(f, "data");
752 if (! array)
753 {
754 DXSetError(ERROR_MISSING_DATA, "no data component found");
755 goto error;
756 }
757
758 if (! DXTypeCheck(array, TYPE_FLOAT, CATEGORY_REAL, 1, iP->P.nDim))
759 {
760 DXSetError(ERROR_DATA_INVALID, "data");
761 return NULL;
762 }
763
764 if (DXQueryConstantArray(array, NULL, NULL))
765 {
766 ip->constantData = 1;
767 ip->vectors = (float *)DXGetConstantArrayData(array);
768 }
769 else
770 {
771 ip->constantData = 0;
772 ip->vectors = (float *)DXGetArrayData(array);
773 }
774
775 return (VectorPart)ip;
776
777 error:
778 if (ip)
779 DXFree((Pointer)ip);
780 return NULL;
781 }
782
783 static Error
Irreg_ResetVectorGrp(VectorGrp P)784 Irreg_ResetVectorGrp(VectorGrp P)
785 {
786 int i;
787
788 for (i = 0; i < P->n; i++)
789 if (! Irreg_ResetVectorPart(P->p[i]))
790 return ERROR;
791
792 return OK;
793 }
794
795 static Error
Irreg_ResetVectorPart(VectorPart p)796 Irreg_ResetVectorPart(VectorPart p)
797 {
798 Irreg_VectorPart ip = (Irreg_VectorPart)p;
799
800 if (ip && ip->visited)
801 memset(ip->visited, 0, ip->nElements);
802
803 return OK;
804 }
805
806 static Error
Irreg_Delete(VectorGrp P)807 Irreg_Delete(VectorGrp P)
808 {
809 Irreg_VectorGrp iP = (Irreg_VectorGrp)P;
810
811 if (iP)
812 {
813 int i;
814
815 for (i = 0; i < iP->P.n; i++)
816 Irreg_DestroyVectorPart(iP->P.p[i]);
817
818 DXFree((Pointer)iP->P.p);
819 DXFree((Pointer)iP);
820 }
821
822 return OK;
823 }
824
825 static Error
Irreg_DestroyVectorPart(VectorPart p)826 Irreg_DestroyVectorPart(VectorPart p)
827 {
828 Irreg_VectorPart ip = (Irreg_VectorPart)p;
829
830 if (ip)
831 {
832 DXFree((Pointer)ip->visited);
833 DXFree((Pointer)ip);
834 }
835
836 return OK;
837 }
838
839 static int
Irreg_FindElement(InstanceVars I,POINT_TYPE * point)840 Irreg_FindElement(InstanceVars I, POINT_TYPE *point)
841 {
842 int i, last, r = 0;
843 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
844 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
845
846 last = iI->cp;
847
848 for (i = 0; i < iP->P.n; i++)
849 if (last != i && iP->P.p[i] != NULL
850 && _dxfIsInBox(iP->P.p[i], point, iP->P.nDim, 0))
851 {
852 if ((r = Irreg_FindElement_VectorPart(iI, i, point)) != 0)
853 break;
854 }
855
856 if (r == -1)
857 return -1;
858 else if (i == iP->P.n)
859 return 0;
860 else
861 {
862 if (I->currentPartition)
863 DXDelete(I->currentPartition);
864 I->currentPartition = DXReference((Object)iP->P.p[i]->field);
865 return 1;
866 }
867 }
868
869 static int
Irreg_FindMultiGridContinuation(InstanceVars I,POINT_TYPE * point)870 Irreg_FindMultiGridContinuation(InstanceVars I, POINT_TYPE *point)
871 {
872 int i, last, r = 0;
873 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
874 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
875
876 last = iI->cp;
877
878 for (i = 0; i < iP->P.n; i++)
879 if (last != i && iP->P.p[i] != NULL
880 && _dxfIsInBox(iP->P.p[i], point, iP->P.nDim, 0))
881 {
882 if ((r = Irreg_FindElement_VectorPart(iI, i, point)) != 0)
883 break;
884 }
885
886 /*
887 * If we don't find an element containing the point, look
888 * for an entry point into the volume - that is, an external
889 * face of an element that the point lies in (within fuzz) and
890 * where the interpolated vector moves into the element.
891 */
892 if (i == iP->P.n)
893 {
894 for (i = 0; i < iP->P.n; i++)
895 if (last != i && iP->P.p[i] != NULL
896 && _dxfIsInBox(iP->P.p[i], point, iP->P.nDim, 1))
897 {
898 r = Irreg_FindMultiGridContinuation_VectorPart(iI, i, point);
899 if (r != 0)
900 break;
901 }
902 }
903
904 if (r == -1)
905 return -1;
906 else if (i == iP->P.n)
907 return 0;
908 else
909 return 1;
910 }
911
912 #define ADD_MINMAX(p, i) \
913 { \
914 if ((p)[i] < min[i]) min[i] = (p)[i]; \
915 if ((p)[i] > max[i]) max[i] = (p)[i]; \
916 }
917
918 static int
Irreg_FindMultiGridContinuation_VectorPart(Irreg_InstanceVars iI,int np,POINT_TYPE * point)919 Irreg_FindMultiGridContinuation_VectorPart(Irreg_InstanceVars iI,
920 int np, POINT_TYPE *point)
921 {
922 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
923 int *elt, i, j, k, l, found;
924 Irreg_VectorPart ip;
925 int ebuf[8], nBuf[8], *nbrs;
926 int nPrimFaces;
927 NbrTable *nptr;
928
929 if (iI->pHandle != NULL)
930 {
931 DXFreeArrayHandle(iI->pHandle);
932 iI->pHandle = NULL;
933 }
934
935 if (iI->cHandle != NULL)
936 {
937 DXFreeArrayHandle(iI->cHandle);
938 iI->cHandle = NULL;
939 }
940
941 if (iI->nHandle)
942 {
943 DXFreeNeighborsHandle(iI->nHandle);
944 iI->nHandle = NULL;
945 }
946
947 iI->cp = np;
948
949 ip = (Irreg_VectorPart)iP->P.p[np];
950 if (! ip)
951 goto error;
952
953 iI->pHandle = DXCreateArrayHandle(ip->pArray);
954 iI->cHandle = DXCreateArrayHandle(ip->cArray);
955 iI->nHandle = DXCreateNeighborsHandle(ip->p.field);
956 if (! iI->pHandle || ! iI->cHandle || ! iI->nHandle)
957 goto error;
958
959 /*
960 * The following is the number of faces of primitive elements
961 * within an element of the current type. There are P.nDim+1
962 * vertices (and therefore faces) if a primitive element of
963 * that dimensionality.
964 */
965 nPrimFaces = iP->P.nDim + 1;
966
967 for (i = 0, found = 0;
968 ! found && i < ip->nElements;
969 i++, elt += iP->vrtsPerElement)
970 {
971 float min[3];
972 float max[3];
973
974 if (ip->invElements && DXIsElementInvalid(ip->invElements, i))
975 continue;
976
977 nbrs = GET_NEIGHBORS(iI->nHandle, i, nBuf);
978 for (j = 0; j < iP->nbrsPerElement; j++)
979 if (nbrs[j] < 0 || (ip->invElements &&
980 DXIsElementInvalid(ip->invElements, nbrs[j])))
981 break;
982
983 if (j == iP->vrtsPerElement)
984 continue;
985
986 elt = (int *)DXGetArrayEntry(iI->cHandle, i, (Pointer)ebuf);
987
988 /*
989 * First check the bounding box with fuzz
990 */
991 min[0] = min[1] = min[2] = DXD_MAX_FLOAT;
992 max[0] = max[1] = max[2] = -DXD_MAX_FLOAT;
993
994
995 for (j = 0; j < iP->vrtsPerElement; j++)
996 {
997 float *p;
998 float pbuf[3];
999
1000 p = (float *)DXGetArrayEntry(iI->pHandle,
1001 elt[j], (Pointer)pbuf);
1002
1003 for (k = 0; k < iP->P.nDim; k++)
1004 ADD_MINMAX(p, k);
1005 }
1006
1007 for (k = 0; k < iP->P.nDim; k++)
1008 {
1009 float f = 0.0001 * (max[k] - min[k]);
1010 if (point[k] < (min[k]-f) || point[k] > (max[k]+f))
1011 break;
1012 }
1013
1014 if (k != iP->P.nDim)
1015 continue;
1016
1017 /*
1018 * Now do it the hard way
1019 */
1020
1021 iI->cp = np;
1022 iI->ce = i;
1023
1024 /*
1025 * for each external face of a primitive element, determine
1026 * whether it lies in the external boundary of the partition.
1027 */
1028 nptr = iP->nbrTable;
1029 for (j = 0; ! found && j < iP->primsPerElement; j++)
1030 {
1031 int first = 1;
1032 iI->ct = j;
1033
1034 /*
1035 * for each face of the primitive, if its external to the
1036 * element and the element has no neighbor in that direction,
1037 * then the primitive face lies in the boundary of the
1038 * partition.
1039 */
1040 for (k = 0; k < nPrimFaces; k++, nptr++)
1041 {
1042 /*
1043 * if the primitive face is internal to the
1044 * element, ignore it.
1045 */
1046 if (nptr->nelt < 0)
1047 continue;
1048
1049 /*
1050 * If there's a valid neighbor in that direction, ignore it.
1051 */
1052 if (nbrs[nptr->nelt] > 0)
1053 if (!ip->invElements ||
1054 !DXIsElementInvalid(ip->invElements, nbrs[nptr->nelt]))
1055 continue;
1056
1057 /*
1058 * Only do the primitive interpolation once.
1059 */
1060 if (first)
1061 {
1062 first = 0;
1063
1064 (*iP->InitElement)(iI);
1065 if (IsCurrentDegenerate(iI))
1066 continue;
1067
1068 (*iP->P.Weights)((InstanceVars)iI, point);
1069 }
1070
1071 if (iI->w[k] < 0.0 && iI->w[k] > -0.01)
1072 {
1073 VECTOR_TYPE v[5];
1074 double dot;
1075 float *plane = iI->planes + 4*k;
1076
1077 iI->w[k] = 0.0;
1078 iI->face = k;
1079
1080 if (! Irreg_Interpolate((InstanceVars)iI, NULL, v))
1081 goto error;
1082
1083 dot = 0.0;
1084 for (l = 0; l < iP->P.nDim; l++)
1085 dot += plane[l] * v[l];
1086
1087 if (dot < 0.0)
1088 found = 1;
1089 }
1090 }
1091 }
1092 }
1093
1094 return found;
1095
1096 error:
1097 return -1;
1098 }
1099
1100 static int
Irreg_FindElement_VectorPart(Irreg_InstanceVars iI,int np,POINT_TYPE * point)1101 Irreg_FindElement_VectorPart(Irreg_InstanceVars iI, int np, POINT_TYPE *point)
1102 {
1103 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1104 int *elt, i, j, k, found;
1105 Irreg_VectorPart ip;
1106 int ebuf[8];
1107
1108 if (iI->pHandle != NULL)
1109 {
1110 DXFreeArrayHandle(iI->pHandle);
1111 iI->pHandle = NULL;
1112 }
1113
1114 if (iI->cHandle != NULL)
1115 {
1116 DXFreeArrayHandle(iI->cHandle);
1117 iI->cHandle = NULL;
1118 }
1119
1120 if (iI->nHandle)
1121 {
1122 DXFreeNeighborsHandle(iI->nHandle);
1123 iI->nHandle = NULL;
1124 }
1125
1126 iI->cp = np;
1127
1128 ip = (Irreg_VectorPart)iP->P.p[np];
1129 if (! ip)
1130 goto error;
1131
1132 iI->pHandle = DXCreateArrayHandle(ip->pArray);
1133 iI->cHandle = DXCreateArrayHandle(ip->cArray);
1134 if (! iI->pHandle || ! iI->cHandle)
1135 goto error;
1136
1137 for (i = 0, found = 0;
1138 ! found && i < ip->nElements;
1139 i++, elt += iP->vrtsPerElement)
1140 {
1141 float min[3];
1142 float max[3];
1143
1144 if (ip->invElements && DXIsElementInvalid(ip->invElements, i))
1145 continue;
1146
1147 elt = (int *)DXGetArrayEntry(iI->cHandle, i, (Pointer)ebuf);
1148
1149 /*
1150 * First check the bounding box
1151 */
1152 min[0] = min[1] = min[2] = DXD_MAX_FLOAT;
1153 max[0] = max[1] = max[2] = -DXD_MAX_FLOAT;
1154
1155
1156 for (j = 0; j < iP->vrtsPerElement; j++)
1157 {
1158 float *p;
1159 float pbuf[3];
1160
1161 p = (float *)DXGetArrayEntry(iI->pHandle,
1162 elt[j], (Pointer)pbuf);
1163
1164 for (k = 0; k < iP->P.nDim; k++)
1165 ADD_MINMAX(p, k);
1166 }
1167
1168 for (k = 0; k < iP->P.nDim; k++)
1169 if (point[k] < min[k] || point[k] > max[k])
1170 break;
1171
1172 if (k != iP->P.nDim)
1173 continue;
1174
1175 /*
1176 * Now do it the hard way
1177 */
1178
1179 iI->cp = np;
1180 iI->ce = i;
1181
1182 for (j = 0; ! found && j < iP->primsPerElement; j++)
1183 {
1184 iI->ct = j;
1185
1186 (*iP->InitElement)(iI);
1187 if (! IsCurrentDegenerate(iI))
1188 if ((*iP->P.Weights)((InstanceVars)iI, point))
1189 found = 1;
1190 }
1191 }
1192
1193 if (found)
1194 {
1195 iI->nHandle = DXCreateNeighborsHandle(ip->p.field);
1196 if (! iI->nHandle)
1197 goto error;
1198 }
1199
1200 return found;
1201
1202 error:
1203 return -1;
1204 }
1205
1206 #define CROSS2PLANE(x, p) \
1207 { \
1208 (x)[3] = -((x)[0]*(p)[0] + (x)[1]*(p)[1] + (x)[2]*(p)[2]); \
1209 }
1210
1211 #define DOT(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]+(b)[2])
1212 #define POINTPLANE(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2]*(b)[2] + a[3])
1213 #define NEGATEPLANE(a) \
1214 { \
1215 (a)[0] = -(a)[0]; \
1216 (a)[1] = -(a)[1]; \
1217 (a)[2] = -(a)[2]; \
1218 (a)[3] = -(a)[3]; \
1219 }
1220
1221 static int
Tetras_Weights(InstanceVars I,POINT_TYPE * pt)1222 Tetras_Weights(InstanceVars I, POINT_TYPE *pt)
1223 {
1224 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
1225 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1226 int e, t;
1227 Irreg_VectorPart ip;
1228 float vPp[3], vPs[3];
1229 float V, nV;
1230 float *p, *s;
1231 float pbuf[3], sbuf[3];
1232 float *xqsr = iI->planes + 0;
1233 float *xrps = iI->planes + 4;
1234 float *xspq = iI->planes + 8;
1235 float *xqpr = iI->planes + 12;
1236 float *w = iI->w;
1237 int i, *elt, *tet;
1238 int ebuf[8];
1239
1240 if (IsCurrentDegenerate(iI))
1241 return 0;
1242
1243 ip = (Irreg_VectorPart)iP->P.p[iI->cp];
1244
1245 /*np = iI->cp;*/
1246 e = iI->ce;
1247 t = iI->ct;
1248
1249 if (ip->invElements && DXIsElementInvalid(ip->invElements, e))
1250 return 0;
1251
1252 elt = (int *)DXGetArrayEntry(iI->cHandle, e, (Pointer)ebuf);
1253
1254 tet = iP->primTable + 4*t;
1255
1256 p = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[0]],(Pointer)pbuf);
1257 /*q = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[1]],(Pointer)qbuf);*/
1258 /*r = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[2]],(Pointer)rbuf);*/
1259 s = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[3]],(Pointer)sbuf);
1260
1261 /*
1262 * DXCross products of first three faces have apexes at p, last has
1263 * apex at s.
1264 */
1265 vPp[0] = pt[0] - p[0]; vPp[1] = pt[1] - p[1]; vPp[2] = pt[2] - p[2];
1266 vPs[0] = pt[0] - s[0]; vPs[1] = pt[1] - s[1]; vPs[2] = pt[2] - s[2];
1267
1268 w[0] = vPs[0]*xqsr[0] + vPs[1]*xqsr[1] + vPs[2]*xqsr[2];
1269 w[1] = vPp[0]*xrps[0] + vPp[1]*xrps[1] + vPp[2]*xrps[2];
1270 w[2] = vPp[0]*xspq[0] + vPp[1]*xspq[1] + vPp[2]*xspq[2];
1271 w[3] = vPp[0]*xqpr[0] + vPp[1]*xqpr[1] + vPp[2]*xqpr[2];
1272
1273 V = w[0] + w[1] + w[2] + w[3];
1274
1275 /*
1276 * If volume is zero, element was degenerate. Otherwise,
1277 * check for a boundary condition: total negative volume must
1278 * be no more than 0.0001 of total volume for the point to be
1279 * considered inside.
1280 */
1281 if (V == 0)
1282 return 0;
1283 else if (V > 0)
1284 {
1285 nV = 0;
1286 for (i = 0; i < 4; i++)
1287 if (w[i] < 0)
1288 nV += w[i];
1289 }
1290 else
1291 {
1292 nV = 0;
1293 for (i = 0; i < 4; i++)
1294 if (w[i] > 0)
1295 nV += w[i];
1296 }
1297
1298 V = 1.0 / V;
1299
1300 for (i = 0; i < 4; i++)
1301 w[i] *= V;
1302
1303 if (nV != 0.0)
1304 return 0;
1305
1306 return 1;
1307 }
1308
1309 static int
Tetras_FaceWeights(InstanceVars I,POINT_TYPE * pt)1310 Tetras_FaceWeights(InstanceVars I, POINT_TYPE *pt)
1311 {
1312 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
1313 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1314 int np, e, t;
1315 Irreg_VectorPart ip;
1316 float w[3];
1317 int i, j, *elt, *tet;
1318 int ebuf[8];
1319 float V, *pts[3], pbuf[3][3];
1320 Vector vecs[3], cross;
1321 float fpt[3];
1322 int npos;
1323
1324 fpt[0] = pt[0];
1325 fpt[1] = pt[1];
1326 fpt[2] = pt[2];
1327
1328 np = iI->cp;
1329 e = iI->ce;
1330 t = iI->ct;
1331
1332 ip = (Irreg_VectorPart)iP->P.p[np];
1333
1334 if (ip->invElements && DXIsElementInvalid(ip->invElements, e))
1335 return 0;
1336
1337 elt = (int *)DXGetArrayEntry(iI->cHandle, e, (Pointer)ebuf);
1338
1339 tet = iP->primTable + 4*t;
1340
1341 for (i = j = 0; i < 4; i++)
1342 {
1343 if (i != iI->face)
1344 {
1345 pts[j] = (float *)DXGetArrayEntry(iI->pHandle,
1346 elt[tet[j]], (Pointer)pbuf[j]);
1347 j++;
1348 }
1349 }
1350
1351 _dxfvector_subtract_3D((Vector *)pts[0], (Vector *)fpt, &vecs[0]);
1352 _dxfvector_subtract_3D((Vector *)pts[1], (Vector *)fpt, &vecs[1]);
1353 _dxfvector_subtract_3D((Vector *)pts[2], (Vector *)fpt, &vecs[2]);
1354
1355 _dxfvector_cross_3D(&vecs[1], &vecs[2], &cross);
1356 w[0] = _dxfvector_length_3D(&cross);
1357
1358 w[1] = _dxfvector_cross_3D(&vecs[2], &vecs[0], &cross);
1359 w[1] = _dxfvector_length_3D(&cross);
1360
1361 w[2] = _dxfvector_cross_3D(&vecs[0], &vecs[1], &cross);
1362 w[2] = _dxfvector_length_3D(&cross);
1363
1364 V = w[0] + w[1] + w[2];
1365
1366 V = 1.0 / V;
1367
1368 for (i = j = npos = 0; i < 4; i++)
1369 {
1370 if (i == iI->face)
1371 iI->w[i] = 0.0;
1372 else
1373 iI->w[i] = w[j++] * V;
1374
1375 if (iI->w[i] >= 0.0)
1376 npos ++;
1377 }
1378
1379 return npos == 0 || npos == 4;
1380 }
1381
1382 static int
Tetras_InitElement(Irreg_InstanceVars iI)1383 Tetras_InitElement(Irreg_InstanceVars iI)
1384 {
1385 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1386 Irreg_VectorPart ip;
1387 float *p, *q, *r, *s, *v;
1388 float pbuf[3], qbuf[3], rbuf[3], sbuf[3];
1389 float *vqp = iI->edges + 0;
1390 float *vrp = iI->edges + 3;
1391 float *vsp = iI->edges + 6;
1392 float *vqs = iI->edges + 9;
1393 float *vrs = iI->edges + 12;
1394 float *vrq = iI->edges + 15;
1395 float *xqsr = iI->planes + 0;
1396 float *xrps = iI->planes + 4;
1397 float *xspq = iI->planes + 8;
1398 float *xqpr = iI->planes + 12;
1399 int i, *elt, ebuf[8];
1400 int np, e, t, *tet;
1401
1402 np = iI->cp;
1403 e = iI->ce;
1404 t = iI->ct;
1405
1406 ip = (Irreg_VectorPart)iP->P.p[np];
1407
1408 elt = (int *)DXGetArrayEntry(iI->cHandle, e, (Pointer)ebuf);
1409
1410 tet = iP->primTable + 4*t;
1411
1412 SetElementVisited(ip, iI->ce);
1413
1414 p = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[0]],(Pointer)pbuf);
1415 q = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[1]],(Pointer)qbuf);
1416 r = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[2]],(Pointer)rbuf);
1417 s = (float*)DXGetArrayEntry(iI->pHandle,elt[tet[3]],(Pointer)sbuf);
1418
1419 /*
1420 * edge vectors
1421 */
1422 vqp[0] = q[0] - p[0]; vqp[1] = q[1] - p[1]; vqp[2] = q[2] - p[2];
1423 vrp[0] = r[0] - p[0]; vrp[1] = r[1] - p[1]; vrp[2] = r[2] - p[2];
1424 vsp[0] = s[0] - p[0]; vsp[1] = s[1] - p[1]; vsp[2] = s[2] - p[2];
1425 vqs[0] = q[0] - s[0]; vqs[1] = q[1] - s[1]; vqs[2] = q[2] - s[2];
1426 vrs[0] = r[0] - s[0]; vrs[1] = r[1] - s[1]; vrs[2] = r[2] - s[2];
1427 vrq[0] = r[0] - q[0]; vrq[1] = r[1] - q[1]; vrq[2] = r[2] - q[2];
1428
1429 /*
1430 * DXCross product terms for the four faces based on vertices p and s
1431 */
1432 xqsr[0] = vrs[1]*vqs[2] - vrs[2]*vqs[1];
1433 xqsr[1] = vrs[2]*vqs[0] - vrs[0]*vqs[2];
1434 xqsr[2] = vrs[0]*vqs[1] - vrs[1]*vqs[0];
1435 CROSS2PLANE(xqsr, s);
1436
1437 xrps[0] = vrp[1]*vsp[2] - vrp[2]*vsp[1];
1438 xrps[1] = vrp[2]*vsp[0] - vrp[0]*vsp[2];
1439 xrps[2] = vrp[0]*vsp[1] - vrp[1]*vsp[0];
1440 CROSS2PLANE(xrps, p);
1441
1442 xspq[0] = vsp[1]*vqp[2] - vsp[2]*vqp[1];
1443 xspq[1] = vsp[2]*vqp[0] - vsp[0]*vqp[2];
1444 xspq[2] = vsp[0]*vqp[1] - vsp[1]*vqp[0];
1445 CROSS2PLANE(xspq, p);
1446
1447 xqpr[0] = vqp[1]*vrp[2] - vqp[2]*vrp[1];
1448 xqpr[1] = vqp[2]*vrp[0] - vqp[0]*vrp[2];
1449 xqpr[2] = vqp[0]*vrp[1] - vqp[1]*vrp[0];
1450 CROSS2PLANE(xqpr, p);
1451
1452 /*
1453 * normalize the edge vectors
1454 */
1455 ResetCurrentDegenerate(iI);
1456 for (i = 0, v = iI->edges; i < 6; i++, v += 3)
1457 {
1458 float d = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
1459
1460 if (d == 0.0)
1461 {
1462 SetCurrentDegenerate(iI);
1463 v[0] = 0.0;
1464 v[1] = 0.0;
1465 v[2] = 0.0;
1466 }
1467 else
1468 {
1469 d = 1.0 / d;
1470 v[0] *= d;
1471 v[1] *= d;
1472 v[2] *= d;
1473 }
1474 }
1475
1476 /*
1477 * Note that the streamline will be found in some non-degenerate element
1478 * before any degenerate elements have to be handled.
1479 */
1480 if (! IsParityKnown(iI) && ! IsCurrentDegenerate(iI))
1481 {
1482 if (POINTPLANE(xqsr, p) > 0.0)
1483 SetFaceNReversed(iI, 0);
1484 if (POINTPLANE(xrps, q) > 0.0)
1485 SetFaceNReversed(iI, 1);
1486 if (POINTPLANE(xspq, r) > 0.0)
1487 SetFaceNReversed(iI, 2);
1488 if (POINTPLANE(xqpr, s) > 0.0)
1489 SetFaceNReversed(iI, 3);
1490
1491 SetParityKnown(iI);
1492 }
1493
1494 if (IsParityKnown(iI))
1495 {
1496 int i;
1497 float *p = iI->planes;
1498
1499 for (i = 0; i < 4; i++, p += 4)
1500 if (IsFaceNReversed(iI, i))
1501 NEGATEPLANE(p);
1502 }
1503
1504 return 1;
1505 }
1506
1507 static int
Tris_Weights(InstanceVars I,POINT_TYPE * pt)1508 Tris_Weights(InstanceVars I, POINT_TYPE *pt)
1509 {
1510 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
1511 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1512 int e, t;
1513 float *p, *q, *r;
1514 float pbuf[3], qbuf[3], rbuf[3];
1515 float x, y, x0, y0, x1, y1, x2, y2;
1516 float A, aInv, nA;
1517 float *w = iI->w;
1518 int i, *elt, *tri, ebuf[8];
1519
1520 if (IsCurrentDegenerate(iI))
1521 return 0;
1522
1523 /*np = iI->cp;*/
1524 e = iI->ce;
1525 t = iI->ct;
1526
1527 /*ip = (Irreg_VectorPart)iP->P.p[np];*/
1528
1529 elt = (int *)DXGetArrayEntry(iI->cHandle, e, (Pointer)ebuf);
1530
1531 tri = iP->primTable + 3*t;
1532
1533 p = (float*)DXGetArrayEntry(iI->pHandle,elt[tri[0]],(Pointer)pbuf);
1534 q = (float*)DXGetArrayEntry(iI->pHandle,elt[tri[1]],(Pointer)qbuf);
1535 r = (float*)DXGetArrayEntry(iI->pHandle,elt[tri[2]],(Pointer)rbuf);
1536
1537 x = pt[0]; y = pt[1];
1538
1539 x0 = p[0]; y0 = p[1];
1540 x1 = q[0]; y1 = q[1];
1541 x2 = r[0]; y2 = r[1];
1542
1543 A = (x0*y1 + x1*y2 + x2*y0 - y0*x1 - y1*x2 - y2*x0) * 0.5;
1544 if (A == 0)
1545 return 0;
1546
1547 aInv = 1.0 / A;
1548
1549 w[0] = 0.5 * (x *y1 + x1*y2 + x2*y - y *x1 - y1*x2 - y2*x ) * aInv;
1550 w[1] = 0.5 * (x0*y + x *y2 + x2*y0 - y0*x - y *x2 - y2*x0) * aInv;
1551 w[2] = 0.5 * (x0*y1 + x1*y + x *y0 - y0*x1 - y1*x - y *x0) * aInv;
1552
1553 nA = 0;
1554 for (i = 0; i < 3; i++)
1555 if (w[i] < 0)
1556 nA -= w[i];
1557
1558 if (nA > 0.0001)
1559 return 0;
1560
1561 return 1;
1562 }
1563
1564 static int
Tris_FaceWeights(InstanceVars I,POINT_TYPE * pt)1565 Tris_FaceWeights(InstanceVars I, POINT_TYPE *pt)
1566 {
1567 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
1568 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1569 int e, t;
1570 float w[2];
1571 int npos, i, j, *elt, *tri, ebuf[8];
1572 float *pts[3], pbuf[3][2];
1573
1574 /*np = iI->cp;*/
1575 e = iI->ce;
1576 t = iI->ct;
1577
1578 /*ip = (Irreg_VectorPart)iP->P.p[np];*/
1579
1580 elt = (int *)DXGetArrayEntry(iI->cHandle, e, (Pointer)ebuf);
1581
1582 tri = iP->primTable + 3*t;
1583
1584 for (i = j = 0; i < 3; i++)
1585 if (i != iI->face)
1586 {
1587 pts[j] = (float *)DXGetArrayEntry(iI->pHandle,
1588 elt[tri[i]], (Pointer)pbuf[j]);
1589 j++;
1590 }
1591
1592 if (pts[0][0] != pts[1][0])
1593 {
1594 w[1] = (pt[0] - pts[0][0]) / (pts[1][0] - pts[0][0]);
1595 w[0] = 1.0 - w[1];
1596 }
1597 else if (pts[0][1] != pts[1][1])
1598 {
1599 w[1] = (pt[1] - pts[0][1]) / (pts[1][1] - pts[0][1]);
1600 w[0] = 1.0 - w[1];
1601 }
1602 else
1603 {
1604 DXSetError(ERROR_INTERNAL, "zero length edge");
1605 return 0;
1606 }
1607
1608 for (i = j = npos = 0; i < 3; i++)
1609 {
1610 if (i == iI->face)
1611 iI->w[i] = 0.0;
1612 else
1613 iI->w[i] = w[j++];
1614
1615 if (iI->w[i] >= 0.0)
1616 npos ++;
1617 }
1618
1619 return npos == 0 || npos == 3;
1620 }
1621
1622 #define POINTLINE(a,b) ((a)[0]*(b)[0] + (a)[1]*(b)[1] + (a)[2])
1623 #define NEGATELINE(a) \
1624 { \
1625 (a)[0] = -(a)[0]; \
1626 (a)[1] = -(a)[1]; \
1627 (a)[2] = -(a)[2]; \
1628 }
1629
1630 static int
Tris_InitElement(Irreg_InstanceVars iI)1631 Tris_InitElement(Irreg_InstanceVars iI)
1632 {
1633 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1634 Irreg_VectorPart ip;
1635 float *p, *q, *r, *v;
1636 float pbuf[3], qbuf[3], rbuf[3];
1637 float *vqp = iI->edges + 0;
1638 float *vrq = iI->edges + 3;
1639 float *vpr = iI->edges + 6;
1640 float *xrq = iI->planes + 0;
1641 float *xpr = iI->planes + 4;
1642 float *xqp = iI->planes + 8;
1643 int i, *elt, ebuf[8];
1644 int np, e, t, *tri;
1645
1646 np = iI->cp;
1647 e = iI->ce;
1648 t = iI->ct;
1649
1650 ip = (Irreg_VectorPart)iP->P.p[np];
1651
1652 elt = (int *)DXGetArrayEntry(iI->cHandle, e, (Pointer)ebuf);
1653
1654 tri = iP->primTable + 3*t;
1655
1656 SetElementVisited(ip, iI->ce);
1657
1658 p = (float*)DXGetArrayEntry(iI->pHandle,elt[tri[0]],(Pointer)pbuf);
1659 q = (float*)DXGetArrayEntry(iI->pHandle,elt[tri[1]],(Pointer)qbuf);
1660 r = (float*)DXGetArrayEntry(iI->pHandle,elt[tri[2]],(Pointer)rbuf);
1661
1662 /*
1663 * edge vectors
1664 */
1665 vqp[0] = q[0] - p[0]; vqp[1] = q[1] - p[1];
1666 vrq[0] = r[0] - q[0]; vrq[1] = r[1] - q[1];
1667 vpr[0] = p[0] - r[0]; vpr[1] = p[1] - r[1];
1668
1669 xqp[0] = vqp[1];
1670 xqp[1] = -vqp[0];
1671 xqp[2] = -(xqp[0]*p[0] + xqp[1]*p[1]);
1672
1673 xrq[0] = vrq[1];
1674 xrq[1] = -vrq[0];
1675 xrq[2] = -(xrq[0]*q[0] + xrq[1]*q[1]);
1676
1677 xpr[0] = vpr[1];
1678 xpr[1] = -vpr[0];
1679 xpr[2] = -(xpr[0]*r[0] + xpr[1]*r[1]);
1680
1681 /*
1682 * normalize the edge vectors
1683 */
1684 ResetCurrentDegenerate(iI);
1685 for (i = 0, v = iI->edges; i < 3; i++, v += 3)
1686 {
1687 float d = v[0]*v[0] + v[1]*v[1];
1688
1689 if (d == 0.0)
1690 {
1691 SetCurrentDegenerate(iI);
1692 v[0] = 0.0;
1693 v[1] = 0.0;
1694 }
1695 else
1696 {
1697 d = 1.0 / d;
1698 v[0] *= d;
1699 v[1] *= d;
1700 }
1701 }
1702
1703 /*
1704 * Note that the streamline will be found in some non-degenerate element
1705 * before any degenerate elements have to be handled.
1706 */
1707 if (! IsParityKnown(iI) && ! IsCurrentDegenerate(iI))
1708 {
1709 if (POINTLINE(xqp, r) > 0.0)
1710 SetFaceNReversed(iI, 0);
1711 if (POINTLINE(xrq, p) > 0.0)
1712 SetFaceNReversed(iI, 1);
1713 if (POINTLINE(xpr, q) > 0.0)
1714 SetFaceNReversed(iI, 2);
1715
1716 SetParityKnown(iI);
1717 }
1718
1719 if (IsParityKnown(iI))
1720 {
1721 int i;
1722 float *p = iI->planes;
1723
1724 for (i = 0; i < 3; i++, p += 4)
1725 if (IsFaceNReversed(iI, i))
1726 NEGATELINE(p);
1727 }
1728
1729 return 1;
1730 }
1731
1732 static Error
Irreg_FindBoundary(InstanceVars I,POINT_TYPE * p,VECTOR_TYPE * v,double * t)1733 Irreg_FindBoundary(InstanceVars I, POINT_TYPE *p, VECTOR_TYPE *v, double *t)
1734 {
1735 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
1736 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1737 double tmin, tface;
1738 int i, j, imin;
1739 float *plane;
1740 double d0, d1, dmax;
1741
1742 plane = iI->planes; imin = -1; tmin = DXD_MAX_FLOAT; dmax = 0.0;
1743 for (i = 0; i < iP->P.nDim+1; i++, plane += 4)
1744 {
1745 d0 = d1 = 0.0;
1746
1747 for (j = 0; j < iP->P.nDim; j++)
1748 {
1749 d0 += (double)plane[j]*v[j];
1750 d1 += (double)plane[j]*p[j];
1751 }
1752
1753 d1 += plane[iP->P.nDim];
1754
1755 if (d0 <= 0.0)
1756 continue;
1757
1758 if (d1 > 0.0)
1759 tface = 0;
1760 else
1761 tface = -d1 / d0;
1762
1763 if (tface >= 0)
1764 if ((tface == tmin && d0 > dmax) || tface < tmin)
1765 {
1766 dmax = d0;
1767 tmin = tface;
1768 imin = i;
1769 }
1770 }
1771
1772
1773 if (imin == -1)
1774 {
1775 DXSetError(ERROR_INTERNAL, "can't find intersection");
1776 return ERROR;
1777 }
1778
1779 iI->face = imin;
1780 *t = tmin;
1781
1782 return OK;
1783 }
1784
1785 static int
Irreg_Interpolate(InstanceVars I,POINT_TYPE * pt,VECTOR_TYPE * vec)1786 Irreg_Interpolate(InstanceVars I, POINT_TYPE *pt, VECTOR_TYPE *vec)
1787 {
1788 Irreg_VectorPart ip;
1789 Irreg_VectorGrp iP;
1790 Irreg_InstanceVars iI;
1791 int *elt, i, j;
1792 float *vdata;
1793 int *prim;
1794 int ebuf[8];
1795
1796 iI = (Irreg_InstanceVars)I;
1797 iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1798
1799 ip = (Irreg_VectorPart)iP->P.p[iI->cp];
1800
1801 if (ip->constantData)
1802 {
1803 for (i = 0; i < iP->P.nDim; i++)
1804 vec[i] = (VECTOR_TYPE)ip->vectors[i];
1805 return OK;
1806 }
1807
1808 if (ip->p.dependency == DEP_ON_POSITIONS)
1809 {
1810 prim = iP->primTable + (iP->P.nDim+1)*iI->ct;
1811
1812 elt = (int*)DXGetArrayEntry(iI->cHandle,iI->ce,(Pointer)ebuf);
1813
1814 vec[0] = 0.0; vec[1] = 0.0; vec[2] = 0.0;
1815
1816 for (i = 0; i < iP->P.nDim+1; i++)
1817 {
1818 vdata = ip->vectors + iP->P.nDim*elt[prim[i]];
1819 for (j = 0; j < iP->P.nDim; j++)
1820 vec[j] += iI->w[i]*vdata[j];
1821 }
1822 }
1823 else
1824 {
1825 vdata = ip->vectors + iP->P.nDim*iI->ce;
1826 for (j = 0; j < iP->P.nDim; j++)
1827 vec[j] = vdata[j];
1828 }
1829
1830 return OK;
1831 }
1832
1833 static Error
Irreg_StepTime(InstanceVars I,double c,VECTOR_TYPE * vec,double * t)1834 Irreg_StepTime(InstanceVars I, double c, VECTOR_TYPE *vec, double *t)
1835 {
1836 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
1837 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
1838 int i, j, nEdges;
1839 double tmin, iprime;
1840 float *edge;
1841
1842 /*
1843 * We are given c, a limit on |t*Vx|/|X| where t is the unknown time step,
1844 * Vx is the projection of V on edge X, and |X| is the length of X.
1845 * For all edges X, find the smallest t such that
1846 *
1847 * c > |t*Vx| / |X| (= t|Vx|/|X|)
1848 *
1849 * For each edge X, we solve for t': the maximum allowable length along
1850 * the edge X:
1851 *
1852 * t' = c|X|/|Vx|
1853 *
1854 * To get Vx, the projection of the vector along the double, we
1855 * scale the vector by the cosine of the angle separating them:
1856 *
1857 * |Vx| = |V| cos theta
1858 *
1859 * where theta is the angle between the edge and the vector. We also
1860 * have that:
1861 *
1862 * cos theta = (V dot X) / |X| |V|
1863 *
1864 * We therefore get:
1865 *
1866 * |Vx| = (V dot X) / |X|
1867 *
1868 * which we plug back into the above:
1869 *
1870 * t' = c|X|/|Vx| = c|X|/((V dot X)/|X|) = c|X|^2 / (V dot X)
1871 *
1872 * But since our edge vectors are normalized,
1873 *
1874 * = c / (V dot X)
1875 */
1876
1877 if (IsCurrentDegenerate(iI))
1878 {
1879 *t = 0.0;
1880 return OK;
1881 }
1882
1883 if (iP->P.nDim == 2)
1884 nEdges = 3;
1885 else
1886 nEdges = 6;
1887
1888 tmin = DXD_MAX_FLOAT;
1889 edge = iI->edges;
1890 for (i = 0; i < nEdges; i++, edge += 3)
1891 {
1892 float dot;
1893
1894 dot = 0;
1895 for (j = 0; j < iP->P.nDim; j++)
1896 dot += vec[j]*edge[j];
1897
1898 if (dot == 0.0)
1899 continue;
1900
1901 if (dot < 0.0)
1902 dot = -dot;
1903
1904 iprime = c / dot;
1905
1906 if (iprime < tmin)
1907 tmin = iprime;
1908 }
1909
1910 *t = tmin;
1911
1912 return OK;
1913 }
1914
1915 static int
_Irreg_Walk(Irreg_InstanceVars iI,Irreg_VectorGrp iP,int * nbrs,int NPerP,int nDim,POINT_TYPE * start,VECTOR_TYPE * vector,POINT_TYPE * target)1916 _Irreg_Walk(Irreg_InstanceVars iI, Irreg_VectorGrp iP,
1917 int *nbrs, /* The element neighbors of elt iI->ce */
1918 int NPerP, /* The number of nbrs of a primitive */
1919 int nDim, /* The dimensionality of the space */
1920 POINT_TYPE *start, VECTOR_TYPE *vector, POINT_TYPE *target)
1921 {
1922 int nbuf[6];
1923 int i, j;
1924 float *plane;
1925 int ne;
1926 NbrTable *nt;
1927 /*
1928 * Check the current element/primitive. If its in, we're done.
1929 */
1930 if ((*iP->P.Weights)((InstanceVars)iI, target))
1931 return WALK_FOUND;
1932
1933 #if 0
1934 plane = iI->planes;
1935 for (i = 0; i < NPerP; i++)
1936 {
1937 float l = 0.0;
1938 for (j = 0; j < NPerP-1; j++)
1939 l += (plane[j]*plane[j]);
1940
1941 if (l == 0)
1942 continue;
1943
1944 l = 1.0 / sqrt(l);
1945
1946 for (j = 0; j < NPerP; j++)
1947 plane[j] *= l;
1948
1949 plane += NPerP;
1950 }
1951 #endif
1952
1953 /*
1954 * Find the exit point for the vector. To do so, intersect
1955 * the vector with the plane of each face, then test the
1956 * intersection point against each other plane.
1957 */
1958 plane = iI->planes;
1959 for (i = 0; i < NPerP; i++)
1960 {
1961 float VDotP = 0.0; /* vector dot plane */
1962 float SDotP = 0.0; /* start dot plane */
1963 float t; /* intersection dist along vector */
1964 POINT_TYPE xit[3]; /* intersection point */
1965
1966 VDotP = 0.0;
1967 for (j = 0; j < nDim; j++)
1968 {
1969 VDotP += (vector[j] * *plane);
1970 SDotP += (start[j] * *plane);
1971 plane ++;
1972 }
1973 SDotP += *plane;
1974 plane++;
1975
1976 /*
1977 * if VDotP <= 0.0 then the vector points away
1978 * from the element face
1979 */
1980 if (VDotP <= 0)
1981 continue;
1982
1983 t = -SDotP / VDotP;
1984 if (t > 1.0)
1985 continue;
1986
1987 for (j = 0; j < nDim; j++)
1988 xit[j] = start[j] + t*vector[j];
1989
1990 /*
1991 * See if the intersection point (xit) actually lies in the
1992 * face.
1993 */
1994 iI->face = i;
1995 if (! (*iP->P.FaceWeights)((InstanceVars)iI, xit))
1996 continue;
1997
1998 if (j == NPerP)
1999 return WALK_ERROR;
2000
2001 nt = iP->nbrTable + NPerP*iI->ct + iI->face;
2002
2003 /*
2004 * If there is no neighbor, then we've found the exit. This
2005 * is the case if the primitive face corresponds to an element
2006 * face (nt->nelt != -1) and there's no corresponding element
2007 * neighbor.
2008 *
2009 * ne is the next element. Its either the same as this one,
2010 * if this is an internal face, or the neighbor if there
2011 * is one. If this face is external to the partition, then
2012 * we'll be returning before its used.
2013 */
2014 ne = iI->ce;
2015 if (nt->nelt != -1 && (ne = nbrs[nt->nelt]) == -1)
2016 {
2017 for (j = 0; j < nDim; j++)
2018 target[j] = xit[j];
2019 return WALK_EXIT;
2020 }
2021
2022 iI->ct = nt->nprim;
2023 if (ne != iI->ce)
2024 {
2025 iI->ce = ne;
2026 nbrs = GET_NEIGHBORS(iI->nHandle, iI->ce, nbuf);
2027 }
2028
2029 (*iP->InitElement)(iI);
2030
2031 if (IsCurrentDegenerate(iI))
2032 {
2033 for (j = 0; j < nDim; j++)
2034 target[j] = xit[j];
2035 return WALK_EXIT;
2036 }
2037
2038 return _Irreg_Walk(iI, iP, nbrs, NPerP, nDim, start, vector, target);
2039 }
2040
2041 /*
2042 * If we fall out, then we didn't find the exit face.
2043 */
2044 return WALK_ERROR;
2045 }
2046
2047 static int
Irreg_Walk(InstanceVars I,POINT_TYPE * start,VECTOR_TYPE * vector,POINT_TYPE * target)2048 Irreg_Walk(InstanceVars I,
2049 POINT_TYPE *start, VECTOR_TYPE *vector, POINT_TYPE *target)
2050 {
2051 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
2052 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
2053 int nbuf[6], *nbrs;
2054
2055 nbrs = GET_NEIGHBORS(iI->nHandle, iI->ce, nbuf);
2056
2057 return _Irreg_Walk(iI, iP, nbrs,
2058 iP->nbrsPerPrim, iP->nDimensions, start, vector, target);
2059 }
2060
2061 static int
Irreg_Neighbor(InstanceVars I,VECTOR_TYPE * v)2062 Irreg_Neighbor(InstanceVars I, VECTOR_TYPE *v)
2063 {
2064 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
2065 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
2066 Irreg_VectorPart ip;
2067 int nbr, nindex, *nbrs, nBuf[6];
2068 float *plane, dot;
2069 int i;
2070
2071 plane = iI->planes + 4*iI->face;
2072
2073 dot = 0;
2074 for (i = 0; i < iP->P.nDim; i++)
2075 dot += *v++ * *plane++;
2076 if (dot < 0)
2077 return 1;
2078
2079 ip = (Irreg_VectorPart)iP->P.p[iI->cp];
2080
2081 nindex = iP->nbrTable[(iP->P.nDim+1)*iI->ct + iI->face].nelt;
2082 iI->ct = iP->nbrTable[(iP->P.nDim+1)*iI->ct + iI->face].nprim;
2083
2084 if (nindex == -1)
2085 {
2086 (*iP->InitElement)(iI);
2087 return 1;
2088 }
2089
2090 nbrs = GET_NEIGHBORS(iI->nHandle, iI->ce, nBuf);
2091 nbr = nbrs[nindex];
2092
2093 if (nbr < 0)
2094 {
2095 nbr = -nbr;
2096
2097 if (ip->partitionNeighbors)
2098 {
2099 Irreg_VectorPart np;
2100
2101 if (ip->gridConnections)
2102 {
2103 /*
2104 * if grid connections, then the partition neighbors
2105 * gives us a pointer to the neighboring partition.
2106 * We figure out the offset from that partitions grid
2107 * strides and the munged indices of the current element.
2108 */
2109 int indices[10], axis, npartition;
2110
2111 npartition = ip->partitionNeighbors[nindex];
2112 if (npartition == -1)
2113 return 0;
2114
2115 /*
2116 * OK. there's a neighboring partition. Determine the
2117 * indices of the current element in the current partition
2118 * and from that determine the indicies in the neighbor.
2119 */
2120 iI->cp = npartition;
2121
2122 np = (Irreg_VectorPart)iP->P.p[npartition];
2123 if (np == NULL)
2124 return 0;
2125
2126 _dxfGetIndices(iI->ce, ip->nd, ip->gstrides, indices);
2127
2128 axis = nindex >> 1;
2129
2130 if (nindex & 0x1)
2131 {
2132 /*
2133 * Then the current element is at the top end of the
2134 * partition's range, so it'll be at the bottom end
2135 * of the neighbor's.
2136 */
2137 indices[axis] = 0;
2138 }
2139 else
2140 {
2141 /*
2142 * Then the current element is at the bottom end of
2143 * the partition's range, so it'll be at the top of the
2144 * neighbor's.
2145 */
2146 indices[axis] = np->gcounts[axis]-1;
2147 }
2148
2149 iI->ce = _dxfGetOffset(ip->nd, np->gstrides, indices);
2150 }
2151 else
2152 {
2153 if (ip->partitionNeighbors[nbr<<1] == -1)
2154 return 0;
2155
2156 iI->cp = ip->partitionNeighbors[nbr<<1];
2157 iI->ce = ip->partitionNeighbors[(nbr<<1)+1];
2158 np = (Irreg_VectorPart)iP->P.p[iI->cp];
2159 }
2160
2161
2162 if (iI->pHandle)
2163 DXFreeArrayHandle(iI->pHandle);
2164 iI->pHandle = DXCreateArrayHandle(np->pArray);
2165
2166 if (iI->cHandle)
2167 DXFreeArrayHandle(iI->cHandle);
2168 iI->cHandle = DXCreateArrayHandle(np->cArray);
2169
2170 if (iI->nHandle)
2171 DXFreeNeighborsHandle(iI->nHandle);
2172 iI->nHandle = DXCreateNeighborsHandle(np->p.field);
2173
2174 if (!iI->pHandle || ! iI->nHandle || ! iI->cHandle)
2175 return -1;
2176
2177 (*iP->InitElement)(iI);
2178 return 1;
2179 }
2180 else
2181 return 0;
2182 }
2183 else
2184 {
2185 if (ip->invElements && DXIsElementInvalid(ip->invElements, nbr))
2186 return 0;
2187
2188 iI->ce = nbr;
2189 (*iP->InitElement)(iI);
2190 return 1;
2191 }
2192 }
2193
2194 typedef struct
2195 {
2196 Irreg_VectorGrp grp;
2197 Irreg_VectorPart part;
2198 CompositeField group;
2199 } CurlMap_Args;
2200
2201 static Error
Irreg_CurlMap(VectorGrp P,MultiGrid mg)2202 Irreg_CurlMap(VectorGrp P, MultiGrid mg)
2203 {
2204 Irreg_VectorGrp iP = (Irreg_VectorGrp)P;
2205 Object curl = NULL;
2206 int i;
2207 CurlMap_Args args;
2208 CompositeField cf = NULL;
2209 Field cfield = NULL;
2210
2211 if (iP->P.n == 1)
2212 {
2213 if (iP->P.p[0]->dependency == DEP_ON_CONNECTIONS)
2214 {
2215 DXSetError(ERROR_DATA_INVALID,
2216 "cannot compute curl when data is dependent on connections");
2217 goto error;
2218 }
2219
2220 cfield = NULL;
2221
2222 if (! Irreg_CurlMap_Field(iP, (Irreg_VectorPart)iP->P.p[0], NULL, &cfield))
2223 goto error;
2224
2225 if (cfield)
2226 {
2227 if (! _dxfDivCurl((Object)cfield, NULL, &curl))
2228 goto error;
2229
2230 if (! curl)
2231 goto error;
2232 }
2233
2234 DXDelete((Object)cfield);
2235 cfield = NULL;
2236 }
2237 else
2238 {
2239 if (! DXCreateTaskGroup())
2240 goto error;
2241
2242 cf = args.group = DXNewCompositeField();
2243 if (! args.group)
2244 goto error;
2245
2246 args.grp = iP;
2247
2248 for (i = 0; i < iP->P.n; i++)
2249 {
2250 args.part = (Irreg_VectorPart)iP->P.p[i];
2251
2252 if (args.part->p.dependency == DEP_ON_CONNECTIONS)
2253 {
2254 DXSetError(ERROR_DATA_INVALID,
2255 "cannot compute curl when data is dependent on connections");
2256 DXAbortTaskGroup();
2257 goto error;
2258 }
2259
2260 if (! DXAddTask(Irreg_CurlMap_Task,
2261 (Pointer)&args, sizeof(args), 1.0))
2262 {
2263 DXAbortTaskGroup();
2264 goto error;
2265 }
2266 }
2267
2268 if (! DXExecuteTaskGroup() || DXGetError() != ERROR_NONE)
2269 goto error;
2270
2271 DXGetMemberCount((Group)cf, &i);
2272
2273 if (i != 0)
2274 {
2275 if (! _dxfDivCurl((Object)cf, NULL, &curl))
2276 goto error;
2277
2278 if (! curl)
2279 goto error;
2280 }
2281
2282 DXDelete((Object)cf);
2283 cf = NULL;
2284 }
2285
2286 if (curl)
2287 if (! DXSetMember((Group)mg, NULL, curl))
2288 goto error;
2289
2290 return OK;
2291
2292 error:
2293
2294 DXDelete((Object)cf);
2295 DXDelete((Object)curl);
2296 DXDelete((Object)cfield);
2297 return ERROR;
2298 }
2299
2300 static Error
Irreg_CurlMap_Task(Pointer ptr)2301 Irreg_CurlMap_Task(Pointer ptr)
2302 {
2303 CurlMap_Args *args = (CurlMap_Args *)ptr;
2304 return Irreg_CurlMap_Field(args->grp, args->part, (Group)args->group, NULL);
2305 }
2306
2307 static Error
Irreg_CurlMap_Field(Irreg_VectorGrp iP,Irreg_VectorPart ip,Group g,Field * f)2308 Irreg_CurlMap_Field(Irreg_VectorGrp iP, Irreg_VectorPart ip, Group g, Field *f)
2309 {
2310 Field newf = NULL;
2311 int *elt;
2312 int i, j;
2313 int ebuf[8];
2314 InvalidComponentHandle activePositions = NULL;
2315 InvalidComponentHandle activeConnections = NULL;
2316 ArrayHandle handle = NULL;
2317
2318 handle = DXCreateArrayHandle(ip->cArray);
2319 if (! handle)
2320 goto error;
2321
2322 newf = (Field)DXCopy((Object)(ip->p.field), COPY_STRUCTURE);
2323 if (! newf)
2324 goto error;
2325
2326 activePositions = DXCreateInvalidComponentHandle((Object)newf,
2327 NULL, "positions");
2328 if (! activePositions)
2329 goto error;
2330
2331 DXSetAllInvalid(activePositions);
2332
2333 activeConnections = DXCreateInvalidComponentHandle((Object)newf,
2334 NULL, "connections");
2335 if (! activeConnections)
2336 goto error;
2337
2338 DXSetAllInvalid(activeConnections);
2339
2340 /*
2341 * Mark the vertices of the alive elements alive
2342 */
2343 for (i = 0; i < ip->nElements; i++)
2344 {
2345 if (IsElementVisited(ip, i))
2346 {
2347 elt = (int *)DXGetArrayEntry(handle, i, (Pointer)ebuf);
2348
2349 for (j = 0; j < iP->vrtsPerElement; j++)
2350 DXSetElementValid(activePositions, elt[j]);
2351
2352 DXSetElementValid(activeConnections, i);
2353 }
2354 }
2355
2356
2357 /*
2358 * For each dead element, if it contains an alive vertex,
2359 * resurrect it UNLESS it was dead to begin with
2360 */
2361 for (i = 0; i < ip->nElements; i++)
2362 {
2363 /*
2364 * If element i is already active, continue
2365 */
2366 if (IsElementVisited(ip, i))
2367 continue;
2368
2369 /*
2370 * If element i is not currently active but not invalid, then
2371 * activate it if it has at least one active position
2372 */
2373 if (!ip->invElements || DXIsElementValid(ip->invElements, i))
2374 {
2375 elt = (int *)DXGetArrayEntry(handle, i, (Pointer)ebuf);
2376
2377 for (j = 0; j < iP->vrtsPerElement; j++)
2378 if (DXIsElementValid(activePositions, elt[j]))
2379 {
2380 DXSetElementValid(activeConnections, i);
2381 break;
2382 }
2383 }
2384 }
2385
2386 /*
2387 * Now go back through, turning on all the vertices of all
2388 * the elements that touch a vertex that an initially alive
2389 * element touches.
2390 */
2391 for (i = 0; i < ip->nElements; i++)
2392 {
2393 if (DXIsElementValid(activeConnections, i))
2394 {
2395 elt = (int *)DXGetArrayEntry(handle, i, (Pointer)ebuf);
2396 for (j = 0; j < iP->vrtsPerElement; j++)
2397 DXSetElementValid(activePositions, elt[j]);
2398 }
2399 }
2400
2401 DXFreeArrayHandle(handle);
2402 handle = NULL;
2403
2404 if (! DXSaveInvalidComponent(newf, activeConnections))
2405 goto error;
2406
2407 if (! DXSaveInvalidComponent(newf, activePositions))
2408 goto error;
2409
2410 if (DXGetComponentValue(newf, "neighbors"))
2411 DXDeleteComponent(newf, "neighbors");
2412
2413 if (! DXEndField(newf))
2414 goto error;
2415
2416 if (! DXCull((Object)newf))
2417 goto error;
2418
2419 if (! DXEmptyField(newf))
2420 {
2421 if (g)
2422 {
2423 if (! DXSetMember((Group)g, NULL, (Object)newf))
2424 goto error;
2425 }
2426 else if (f)
2427 {
2428 *f = newf;
2429 }
2430 else
2431 {
2432 DXSetError(ERROR_INTERNAL, "no place to return the curl field");
2433 DXDelete((Object)newf);
2434 return ERROR;
2435 }
2436 }
2437 else
2438 DXDelete((Object)newf);
2439
2440 if (activePositions)
2441 DXFreeInvalidComponentHandle(activePositions);
2442
2443 if (activeConnections)
2444 DXFreeInvalidComponentHandle(activeConnections);
2445
2446 return OK;
2447
2448 error:
2449 if (activePositions)
2450 DXFreeInvalidComponentHandle(activePositions);
2451 if (activeConnections)
2452 DXFreeInvalidComponentHandle(activeConnections);
2453 DXFreeArrayHandle(handle);
2454 DXDelete((Object)newf);
2455
2456 return ERROR;
2457 }
2458
2459 static int
Irreg_Ghost(InstanceVars I,POINT_TYPE * p)2460 Irreg_Ghost(InstanceVars I, POINT_TYPE *p)
2461 {
2462 Irreg_InstanceVars iI = (Irreg_InstanceVars)I;
2463 Irreg_VectorGrp iP = (Irreg_VectorGrp)(iI->i.currentVectorGrp);
2464 Irreg_VectorPart vp = (Irreg_VectorPart)iP->P.p[iI->cp];
2465 return vp->p.ghosts[iI->ce];
2466 }
2467
2468