/***********************************************************************/
/* Open Visualization Data Explorer                                    */
/* (C) Copyright IBM Corp. 1989,1999                                   */
/* ALL RIGHTS RESERVED                                                 */
/* This code licensed under the                                        */
/*    "IBM PUBLIC LICENSE - Open Visualization Data Explorer"          */
/***********************************************************************/

#include <dxconfig.h>



#include <dx/dx.h>
#include <math.h>
#include "stream.h"

/*
 ***********************************
 * Reg interface
 ***********************************
 */

typedef struct reg_InstanceVars
{
    struct instanceVars i;
    int 	        indices[3];
    VECTOR_TYPE	        ivector[3];
    double	        weights[9];
    int		        face;
    int		        cp;
    int			ce;
} *Reg_InstanceVars;

typedef struct reg_VectorGrp  *Reg_VectorGrp;
typedef struct reg_VectorPart *Reg_VectorPart;

struct reg_VectorGrp
{
    struct vectorGrp  P;
    Object	      object;
    Matrix  	      mInv;
};

struct reg_VectorPart
{
    struct vectorPart      p;
    int	   	           cnts[3];
    int		           strides[3];
    int		           cstrides[3];
    int		           bumps[8];
    float   	           dels[9];
    float	           org[3];
    Matrix  	           mInv;
    float  	           *vectors;
    int			   constantData;
    int	   	           *partitionNeighbors;
    int		           flag;
    InvalidComponentHandle invElements;
};

       VectorGrp    _dxfReg_InitVectorGrp(Object, char *);
static VectorPart   Reg_InitVectorPart(Field, Reg_VectorGrp, int);
static Error        Reg_DestroyVectorPart(VectorPart);
static int          Reg_FindElement_VectorPart(Reg_InstanceVars, int,
						POINT_TYPE *);
static int          Reg_FindMultiGridContinuation_VectorPart(Reg_InstanceVars,
						int, POINT_TYPE *);
static InstanceVars Reg_NewInstanceVars(VectorGrp);
static Error        Reg_FreeInstanceVars(InstanceVars);
static int          Reg_FindElement(InstanceVars, POINT_TYPE *);
static int          Reg_FindMultiGridContinuation(InstanceVars, POINT_TYPE *);
static int          Reg_Interpolate(InstanceVars, POINT_TYPE *, VECTOR_TYPE *);
static Error        Reg_StepTime(InstanceVars, double, VECTOR_TYPE *, double *);
static Error        Reg_FindBoundary(InstanceVars, 
				POINT_TYPE *, VECTOR_TYPE *, double *);
static int          Reg_Neighbor(InstanceVars, VECTOR_TYPE *);
static Error        Reg_Delete(VectorGrp);
static int          Reg_Weights(InstanceVars, POINT_TYPE *);
static int          Reg_FaceWeights(InstanceVars, POINT_TYPE *);
static Error        Reg_CurlMap(VectorGrp, MultiGrid);

static Matrix	    InvertN(Matrix, int);
static void	    ApplyN(Matrix, int, POINT_TYPE *, POINT_TYPE *);
static void	    ApplyRotationN(Matrix, int, VECTOR_TYPE *, VECTOR_TYPE *);
static Matrix	    GetXYZtoIJKMatrix(Reg_VectorPart, int);


static InstanceVars
Reg_NewInstanceVars(VectorGrp p)
{
    Reg_InstanceVars iI;

    iI = (Reg_InstanceVars)DXAllocate(sizeof(struct reg_InstanceVars));
    memset(iI, -1, sizeof(struct reg_InstanceVars));

    iI->i.currentVectorGrp = p;

    return (InstanceVars)iI;
}

static Error
Reg_FreeInstanceVars(InstanceVars I)
{
    if (I)
	DXFree((Pointer)I);
}

VectorGrp
_dxfReg_InitVectorGrp(Object object, char *elementType)
{
    Reg_VectorGrp   P = NULL;
    Reg_VectorPart  p = NULL;
    int                i, nP;

    if (! _dxfPartitionNeighbors(object))
	 return NULL;

    if (DXGetObjectClass(object) == CLASS_FIELD)
	nP = 1;
    else
	for (nP = 0; DXGetEnumeratedMember((Group)object, nP, NULL); nP++);
    
    P = (Reg_VectorGrp)DXAllocate(sizeof(struct reg_VectorGrp));
    if (! P)
	goto error;

    P->object = object;
    
    P->P.n = nP;

    P->P.p = (VectorPart *)DXAllocateZero(nP * sizeof(struct reg_VectorPart));
    if (! P->P.p)
	goto error;
    
    /*
     * Attach generic regular element type methods to 
     * return object
     */
    P->P.NewInstanceVars   		= Reg_NewInstanceVars;
    P->P.FreeInstanceVars  		= Reg_FreeInstanceVars;
    P->P.FindElement       		= Reg_FindElement;
    P->P.FindMultiGridContinuation      = Reg_FindMultiGridContinuation;
    P->P.Interpolate       		= Reg_Interpolate;
    P->P.StepTime          		= Reg_StepTime;
    P->P.FindBoundary      		= Reg_FindBoundary;
    P->P.Neighbor          		= Reg_Neighbor;
    P->P.CurlMap           		= Reg_CurlMap;
    P->P.Weights           		= Reg_Weights;
    P->P.FaceWeights       		= Reg_FaceWeights;
    P->P.Delete            		= Reg_Delete;
    P->P.Reset             		= NULL;

    if (! strcmp(elementType, "cubes"))
	P->P.nDim = 3;
    else if (! strcmp(elementType, "quads"))
	P->P.nDim = 2;
    else
    {
	DXSetError(ERROR_INVALID_DATA, "data");
	goto error;
    }

    if (DXGetObjectClass(object) == CLASS_FIELD)
    {
	P->P.p[0] = Reg_InitVectorPart((Field)object, P, 1);
	if (! P->P.p[0])
	    goto error;

	P->P.p[0]->field = (Field)object;
    }
    else
    {
	Group g = (Group)object;
	Field f;

	i = 0;
	while (NULL != (f = (Field)DXGetEnumeratedMember(g, i, NULL)))
	{
	    P->P.p[i] = Reg_InitVectorPart(f, P, (i == 0));
	    if (! P->P.p[i] && DXGetError() != ERROR_NONE)
	        goto error;
	    
	    P->P.p[i]->field = f;
	    
	    i++;
	}
    }

    return (VectorGrp)P;

error:
    Reg_Delete((VectorGrp)P);
    return NULL;
}

static VectorPart
Reg_InitVectorPart(Field f, Reg_VectorGrp P, int flag)
{
    Reg_VectorPart    ip = NULL;
    Array 	      array;
    int		      nDim;
    int		      meshOff[3];
    Object	      attr;

    if (DXEmptyField(f))
	return NULL;

    ip = (Reg_VectorPart)DXAllocate(sizeof(struct reg_VectorPart));
    if (! ip)
	goto error;

    attr = DXGetComponentAttribute(f, "data", "dep");
    if (! attr | DXGetObjectClass(attr) != CLASS_STRING)
    {
	DXSetError(ERROR_INVALID_DATA, "data dependency");
	goto error;
    }

    if (!strcmp(DXGetString((String)attr), "positions"))
	ip->p.dependency = DEP_ON_POSITIONS;
    else
	ip->p.dependency = DEP_ON_CONNECTIONS;

    array = (Array)DXGetComponentValue(f, "connections");
    if (! array)
	goto error;

    if (! DXQueryGridConnections(array, NULL, ip->cnts))
    {
	DXSetError(ERROR_INVALID_DATA, "irregular connections");
	goto error;
    }

    if (! DXGetMeshOffsets((MeshArray)array, meshOff))
	goto error;

    array = (Array)DXGetComponentValue(f, "positions");
    if (! array)
	goto error;
    
    if (! DXQueryGridPositions(array, NULL, ip->cnts, ip->org, ip->dels))
    {
	DXSetError(ERROR_INVALID_DATA, "irregular positions");
	goto error;
    }

    if (flag)
	P->mInv = GetXYZtoIJKMatrix(ip, P->P.nDim);

    array = (Array)DXGetComponentValue(f, "partition neighbors");
    if (array)
	ip->partitionNeighbors = (int *)DXGetArrayData(array);
    else
	ip->partitionNeighbors = NULL;
    
    array = (Array)DXGetComponentValue(f, "data");
    if (! array)
    {
	DXSetError(ERROR_MISSING_DATA, "no data component found");
	goto error;
    }

    if (! DXTypeCheck(array, TYPE_FLOAT, CATEGORY_REAL, 1, P->P.nDim))
    {
	DXSetError(ERROR_INVALID_DATA,
	   "dimensionality of positions differs from that of vector field");
	return NULL;
    }

    if (DXQueryConstantArray(array, NULL, NULL))
    {
	ip->vectors = (float *)DXGetConstantArrayData(array);
	ip->constantData = 1;
    }
    else
    {
	ip->vectors = (float *)DXGetArrayData(array);
	ip->constantData = 0;
    }

    if (! _dxfMinMaxBox(f, ip->p.min, ip->p.max))
	goto error;
    
    if (! DXInvalidateConnections((Object)f))
	goto error;
    
    DXGetComponentValue(f, "invalid connections");
    if (DXGetComponentValue(f, "invalid connections"))
    {
	ip->invElements = DXCreateInvalidComponentHandle((Object)f,
						    NULL, "connections");
	if (! ip->invElements)
	    goto error;
    }
    else
	ip->invElements = NULL;
    
    if (P->P.nDim == 3)
    {
	ip->cstrides[0] = (ip->cnts[2] - 1) * (ip->cnts[1] - 1);
	ip->cstrides[1] = (ip->cnts[2] - 1);
	ip->cstrides[2] = 1;
    }
    else
    {
	ip->cstrides[0] = (ip->cnts[1] - 1);
	ip->cstrides[1] = 1;
    }

    if (ip->p.dependency == DEP_ON_CONNECTIONS)
    {
	if (P->P.nDim == 3)
	{
	    ip->strides[0] = (ip->cnts[2] - 1) * (ip->cnts[1] - 1);
	    ip->strides[1] = (ip->cnts[2] - 1);
	    ip->strides[2] = 1;

	    ip->mInv = P->mInv;
	    ip->mInv.b[0] -= meshOff[0];
	    ip->mInv.b[1] -= meshOff[1];
	    ip->mInv.b[2] -= meshOff[2];
	}
	else
	{
	    ip->strides[0] = (ip->cnts[1] - 1);
	    ip->strides[1] = 1;

	    ip->mInv = P->mInv;
	    ip->mInv.b[0] -= meshOff[0];
	    ip->mInv.b[1] -= meshOff[1];
	}
    }
    else
    {
	if (P->P.nDim == 3)
	{
	    ip->strides[0] = ip->cnts[2] * ip->cnts[1];
	    ip->strides[1] = ip->cnts[2];
	    ip->strides[2] = 1;

	    ip->bumps[0] = 0;
	    ip->bumps[1] =                                   ip->strides[0];
	    ip->bumps[2] =                  ip->strides[1];
	    ip->bumps[3] =                  ip->strides[1] + ip->strides[0];
	    ip->bumps[4] = ip->strides[2];
	    ip->bumps[5] = ip->strides[2] +                  ip->strides[0];
	    ip->bumps[6] = ip->strides[2] + ip->strides[1];
	    ip->bumps[7] = ip->strides[2] + ip->strides[1] + ip->strides[0];

	    ip->mInv = P->mInv;
	    ip->mInv.b[0] -= meshOff[0];
	    ip->mInv.b[1] -= meshOff[1];
	    ip->mInv.b[2] -= meshOff[2];
	}
	else
	{
	    ip->strides[0] = ip->cnts[1];
	    ip->strides[1] = 1;

	    ip->bumps[0] = 0;
	    ip->bumps[1] =                  ip->strides[0];
	    ip->bumps[2] = ip->strides[1];
	    ip->bumps[3] = ip->strides[1] + ip->strides[0];

	    ip->mInv = P->mInv;
	    ip->mInv.b[0] -= meshOff[0];
	    ip->mInv.b[1] -= meshOff[1];
	}
    }
	
    ip->flag = 0;

    return (VectorPart)ip;

error:
    if (ip)
	DXFree((Pointer)ip);
    return NULL;
}
    
static Matrix
GetXYZtoIJKMatrix(Reg_VectorPart p, int n)
{
    Matrix         m;

    if (n == 3)
    {
	memcpy((float *)m.A, p->dels, 9*sizeof(float));
	memcpy((float *)m.b, p->org, 3*sizeof(float));

	return InvertN(m, 3);
    }
    else
    {
	m.A[0][0] = p->dels[0];
	m.A[0][1] = p->dels[1];
	m.A[0][2] = 0.0;
	m.A[1][0] = p->dels[2];
	m.A[1][1] = p->dels[3];
	m.A[1][2] = 0.0;
	m.A[2][0] = 0.0;
	m.A[2][1] = 0.0;	
	m.A[2][2] = 0.0;

	m.b[0] = p->org[0];
	m.b[1] = p->org[1];
	m.b[2] = 0.0;

	return InvertN(m, 2);

    }
}
    
static Error
Reg_Delete(VectorGrp P)
{
    Reg_VectorGrp  iP = (Reg_VectorGrp)P;
    Reg_VectorPart ip;

    if (iP)
    {
	int i;

	for (i = 0; i < iP->P.n; i++)
	    if (iP->P.p[i])
		Reg_DestroyVectorPart(iP->P.p[i]);

	DXFree((Pointer)iP->P.p);
	DXFree((Pointer)iP);
    }

    return OK;
}

static Error 
Reg_DestroyVectorPart(VectorPart p)
{
    Reg_VectorPart ip = (Reg_VectorPart)p;

    if (ip)
    {
	if (ip->invElements)
	    DXFreeInvalidComponentHandle(ip->invElements);

	DXFree((Pointer)ip);
    }

    return OK;
}

static int
Reg_FindElement(InstanceVars I, POINT_TYPE *point)
{
    int i;
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)I->currentVectorGrp;

    for (i = 0; i < iP->P.n; i++)
	if (_dxfIsInBox(iP->P.p[i], point, iP->P.nDim, 0))
	{
	    if (Reg_FindElement_VectorPart(iI, i, point))
		break;
	}
    
    if (i == iP->P.n)
	return 0;
    else
	return 1;
}

static int
Reg_FindMultiGridContinuation(InstanceVars I, POINT_TYPE *point)
{
    int i;
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)I->currentVectorGrp;

    for (i = 0; i < iP->P.n; i++)
	if (_dxfIsInBox(iP->P.p[i], point, iP->P.nDim, 1))
	{
	    if (Reg_FindMultiGridContinuation_VectorPart(iI, i, point))
		break;
	}
    
    if (i == iP->P.n)
	return 0;
    else
	return 1;
}

static int
Reg_FindElement_VectorPart(Reg_InstanceVars iI, int np, POINT_TYPE *point)
{
    Reg_VectorGrp  iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart ip = (Reg_VectorPart)iP->P.p[np];
    POINT_TYPE	   cpoint[3];
    int		   nd, i, base;
    int		   *strd = ip->cstrides;

    if (! ip)
	return -1;
    
    ApplyN(ip->mInv, (nd = iP->P.nDim), point, cpoint);

    for (i = 0; i < nd; i++)
	if (cpoint[i] < 0.0 || cpoint[i] >= (ip->cnts[i]-1))
	    break;
    
    if (i != nd)
	return 0;
	
    for (i = 0, base = 0; i < nd; i++)
    {
	iI->indices[i] = (int)cpoint[i];
	base += iI->indices[i]*(*strd++);
    }
	
    if (ip->invElements && DXIsElementInvalid(ip->invElements, base))
	    return 0;
    
    iI->cp = np;
    iI->ce = base;

    ((Reg_VectorPart)(iP->P.p[np]))->flag = 1;
    Reg_Weights((InstanceVars)iI, point);
    return 1;
}

static int
Reg_FindMultiGridContinuation_VectorPart(Reg_InstanceVars iI,
					int np, POINT_TYPE *point)
{
    Reg_VectorGrp  iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart ip = (Reg_VectorPart)iP->P.p[np];
    POINT_TYPE	   cpoint[3];
    int		   nd, i, base;
    int		   *strd = ip->cstrides;

    if (! ip)
	return -1;
    
    ApplyN(ip->mInv, (nd = iP->P.nDim), point, cpoint);

    /*
     * If we passed the bounding box test, we clamp to
     * the grid's boundaries
     */
    for (i = 0; i < nd; i++)
	if (cpoint[i] < 0.0)
	    cpoint[i] = 0.0;
	else if (cpoint[i] >= (ip->cnts[i]-1))
	    cpoint[i] = (ip->cnts[i]-1.00001);
	
    if (i != nd)
	return 0;
	
    for (i = 0, base = 0; i < nd; i++)
    {
	iI->indices[i] = (int)cpoint[i];
	base += iI->indices[i]*(*strd++);
    }
	
    if (ip->invElements && DXIsElementInvalid(ip->invElements, base))
	    return 0;
    
    iI->cp = np;
    iI->ce = base;

    ((Reg_VectorPart)(iP->P.p[np]))->flag = 1;
    Reg_Weights((InstanceVars)iI, point);
    return 1;
}

static int
Reg_Weights(InstanceVars I, POINT_TYPE *pt)
{
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart   ip = (Reg_VectorPart)iP->P.p[iI->cp];
    POINT_TYPE	     dels[3], ipoint[3];
    int		     nd, i;
    int		     *strd = ip->cstrides;
    int		     base  = 0;
    int		     indices[3];

    if (! ip)
	return -1;
    
    ApplyN(ip->mInv, (nd = iP->P.nDim), pt, ipoint);

    for (i = 0; i < nd; i++)
    {
	if (ipoint[i] < -0.0001 || ipoint[i] > (ip->cnts[i]-1) + 0.0001)
	    break;

	if (ipoint[i] < 0.0)
	    ipoint[i] = 0.0;
	
	if (ipoint[i] >= (ip->cnts[i]-1))
	{
	    ipoint[i] = (ip->cnts[i]-1);
	    indices[i] = (ip->cnts[i]-2);
	    dels[i] = 1.0;
	}
	else
	{
	    indices[i] = (int)(ipoint[i]);
	    dels[i] = ipoint[i] - indices[i];
	}
    }

    if (i < nd)
	return 0;
    
    for (i = 0; i < nd; i++)
	base += ((int)ipoint[i]) * *strd++;
    
    if (base != iI->ce)
	return 0;

    if (nd == 3)
    {
	double *w = iI->weights;
	double dx = dels[2];
	double dy = dels[1];
	double dz = dels[0];
	double A  = dx * dy;
	double B  = dx * dz;
	double C  = dy * dz;
	double D  = dx * dy * dz;

	*w++ = 1.0 - dz - dy - dx + A + B + C - D;
	*w++ = dz - B - C + D;
	*w++ = dy - C - A + D;
	*w++ = C - D;
	*w++ = dx - B - A + D;
	*w++ = B - D;
	*w++ = A - D;
	*w++ = D;
    }
    else if (nd == 2)
    {
	double dx   = dels[0];
	double dy   = dels[1];
	double dxdy = dx * dy;
	double *w   = iI->weights;

	*w++ = 1 - dx - dy + dxdy;
	*w++ = dx - dxdy;
	*w++ = dy - dxdy;
	*w++ = dxdy;
    }
    else
    {
	double dx   = dels[0];
	double *w   = iI->weights;

	*w++ = 1 - dx;
	*w++ = dx;
    }
    
    return 1;
}

static int
Reg_FaceWeights(InstanceVars I, POINT_TYPE *pt)
{
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart   ip = (Reg_VectorPart)iP->P.p[iI->cp];
    POINT_TYPE	     dels[3], ipoint[3];
    int		     nd, i;

    if (! ip)
	return -1;
    
    ApplyN(ip->mInv, (nd = iP->P.nDim), pt, ipoint);

    for (i = 0; i < nd; i++)
    {
	if (iI->face == (i<<1)+0)
	{
	    dels[i]        = 0.0;
	}
	else if (iI->face == (i<<1)+1)
	{
	    dels[i]        = 1.0;
	}
	else
	{
	    if (ipoint[i] < 0.0)
	    {
		dels[i] = 0.0;
	    } 
	    else if (ipoint[i] >= (ip->cnts[i]-1))
	    {
		dels[i] = 1.0;
	    }
	    else
	    {
		dels[i] = ipoint[i] - iI->indices[i];
	    }
	}
    }

    if (i < nd)
	return 0;

    if (nd == 3)
    {
	double *w = iI->weights;
	double dx = dels[2];
	double dy = dels[1];
	double dz = dels[0];
	double A  = dx * dy;
	double B  = dx * dz;
	double C  = dy * dz;
	double D  = dx * dy * dz;

	*w++ = 1.0 - dz - dy - dx + A + B + C - D;
	*w++ = dz - B - C + D;
	*w++ = dy - C - A + D;
	*w++ = C - D;
	*w++ = dx - B - A + D;
	*w++ = B - D;
	*w++ = A - D;
	*w++ = D;
    }
    else if (nd == 2)
    {
	double dx   = dels[0];
	double dy   = dels[1];
	double dxdy = dx * dy;
	double *w   = iI->weights;

	*w++ = 1 - dx - dy + dxdy;
	*w++ = dx - dxdy;
	*w++ = dy - dxdy;
	*w++ = dxdy;
    }
    else
    {
	double dx   = dels[0];
	double *w   = iI->weights;

	*w++ = 1 - dx;
	*w++ = dx;
    }
    
    return 1;
}

static Error
Reg_FindBoundary(InstanceVars I, POINT_TYPE *p, POINT_TYPE *v, double *t)
{
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart   ip = (Reg_VectorPart)iP->P.p[iI->cp];
    int   i, nd;
    POINT_TYPE tmin, ipoint[3];
    int   face = -1, f;

    ApplyN(ip->mInv, (nd = iP->P.nDim), p, ipoint);

    tmin = DXD_MAX_FLOAT;
    for (i = 0; i < nd; i++)
    {
	double iv = iI->ivector[i];
	double l, t0;

	l  = ipoint[i] - iI->indices[i];

	if (l < 0.0) l = 0.0;
	else if (l > 1.0) l = 1.0;

	if (iv < 0)
	{
	    iv = -iv;
	    f  = (i << 1);
	}
	else if (iv > 0)
	{
	    l = 1.0 - l;
	    f  = (i << 1) + 1;
	}
	else
	    continue;

	t0 = l / iv;
	if (t0 < tmin)
	{
	    tmin = t0;
	    face = f;
	}
    }

    if (tmin == DXD_MAX_FLOAT || tmin < 0.0 || face == -1)
    {
        tmin = 0.0;
	iI->face = 0;
	return ERROR;
    }

    iI->face = face;
    *t = tmin;

    return OK;
}
	
static int
Reg_Interpolate(InstanceVars I, POINT_TYPE *pt, VECTOR_TYPE *vec)
{
    Reg_InstanceVars iI    = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP    = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart   ip    = (Reg_VectorPart)iP->P.p[iI->cp];
    int		     nDim  = iP->P.nDim;
    int		     nVrts = 1 << nDim;
    int		     *ind  = iI->indices;
    int		     *strd = ip->strides;
    int		     *bump = ip->bumps;
    double	     *w    = iI->weights;
    int		     base;
    float	     *vector;
    int		     i, j;

    if (ip->constantData)
    {
	for (i = 0; i < nDim; i++)
	    vec[i] = ip->vectors[i];
	return OK;
    }

    if (ip->p.dependency == DEP_ON_POSITIONS)
    {
	base = 0;
	for (i = 0; i < nDim; i++)
	    base += *ind++ * *strd++;
	
	vec[0] = vec[1] = vec[2] = 0.0;
	for (i = 0; i < nVrts; i++, bump++, w++)
	{
	    vector = ip->vectors + nDim*(base + *bump);
	    for (j = 0; j < nDim; j++)
		vec[j] += *w * *vector++;
	}
    }
    else
    {
	base = 0;
	for (i = 0; i < nDim; i++)
	    base += *ind++ * *strd++;
    
	vector = ip->vectors + nDim*base;
	for (j = 0; j < nDim; j++)
	    vec[j] = *vector++;
    }

    return OK;
}

static Error
Reg_StepTime(InstanceVars I, double c, VECTOR_TYPE *vec, double *t)
{
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart   ip = (Reg_VectorPart)iP->P.p[iI->cp];
    int              i, nd;
    double 	     l;

    /*
     * determine the step in time that results in a step along
     * the given step vector of length c.
     */
    
    ApplyRotationN(ip->mInv, (nd = iP->P.nDim), vec, iI->ivector);

    l = 0.0;
    for (i = 0; i < nd; i++)
	l += (iI->ivector[i] * iI->ivector[i]);
    
    if (l == 0.0)
	*t = 0.0;
    else
	*t = c / sqrt(l);

    return OK;
}
	
static int
Reg_Neighbor(InstanceVars I, VECTOR_TYPE *v)
{
    Reg_InstanceVars iI = (Reg_InstanceVars)I;
    Reg_VectorGrp    iP = (Reg_VectorGrp)iI->i.currentVectorGrp;
    Reg_VectorPart   ip = (Reg_VectorPart)iP->P.p[iI->cp];
    int 	     face, axis, base;
    int 	     *ind;
    int 	     *strd;
    int 	     i, nd;

    axis = iI->face >> 1;
    face = iI->face & 0x1;

    /*
     * Look to see whether this is an internal boundary arising
     * from an invalidated connection or an external boundary.
     */
    if ((face == 0 && iI->indices[axis] == 0) ||
        (face == 1 && iI->indices[axis] == (ip->cnts[axis]-2)))
    {
	int nbr;

	if (! ip->partitionNeighbors)
	    return 0;

	nbr = ip->partitionNeighbors[iI->face];

	if (nbr < 0)
	{
	    return 0;
	}

	iI->cp = nbr;

	ip = (Reg_VectorPart)iP->P.p[iI->cp];

	/*
	 * Find the element of the neighbor we are moving into
	 */
	ind  = iI->indices;
	strd = ip->cstrides;
	nd   = iP->P.nDim;

	if (face == 0)
	    ind[axis] = ip->cnts[axis] - 2;
	else
	    ind[axis] = 0;

	for (base = 0, i = 0; i < nd; i++)
	    base += *ind++ * *strd++;
	    
	/*
	 * Check to see if we will be entering an invalidated element
	 */
	if (ip->invElements)
	    if (DXIsElementInvalid(ip->invElements, base))
		return 0;

	iI->ce = base;
	ip->flag = 1;
	return 1;
    }
    else
    {
	ind  = iI->indices;
	strd = ip->cstrides;
	nd   = iP->P.nDim;

	for (base = 0, i = 0; i < nd; i++)
	    base += *ind++ * *strd++;
	
	if (face == 0)
	{
	    base -= ip->cstrides[axis];
	    iI->indices[axis] -= 1;
	}
	else
	{
	    base += ip->cstrides[axis];
	    iI->indices[axis] += 1;
	}

	if (ip->invElements != NULL)
	    if (DXIsElementInvalid(ip->invElements, base))
		return 0;
	
	iI->ce = base;
	return 1;
    }
}

static Matrix
InvertN(Matrix mIn, int nDim)
{
    if (nDim == 3)
	return DXInvert(mIn);
    else
    {
	Matrix mOut;
	float d;

	d = 1.0 / (mIn.A[0][0]*mIn.A[1][1] - mIn.A[0][1]*mIn.A[1][0]);

	mOut.A[0][0] =  mIn.A[1][1] * d;
	mOut.A[0][1] = -mIn.A[1][0] * d;
	mOut.A[1][0] = -mIn.A[0][1] * d;
	mOut.A[1][1] =  mIn.A[0][0] * d;

	mOut.b[0] = -(mIn.b[0]*mOut.A[0][0] + mIn.b[1]*mOut.A[0][1]);
	mOut.b[1] = -(mIn.b[0]*mOut.A[1][0] + mIn.b[1]*mOut.A[1][1]);

	return mOut;
    }
}

static void
ApplyN(Matrix m, int n, POINT_TYPE *in, POINT_TYPE *out)
{
    int i, j;

    for (i = 0; i < n; i++)
    {
	out[i] = m.b[i];
	for (j = 0; j < n; j++)
	    out[i] += in[j] * m.A[j][i];
	
    }

    return;
}

static void
ApplyRotationN(Matrix m, int n, VECTOR_TYPE *in, VECTOR_TYPE *out)
{
    int i, j;

    for (i = 0; i < n; i++)
    {
	out[i] = 0.0;
	for (j = 0; j < n; j++)
	    out[i] += in[j] * m.A[j][i];
    }

    return;
}

extern Object _dxfDivCurl(Object, Object *, Object *);

static Error
Reg_CurlMap(VectorGrp P, MultiGrid mg)
{
    Reg_VectorGrp  iP = (Reg_VectorGrp)P;
    int 	   i;
    Object 	   curl;

    if (P->n == 1)
    {
	if (! _dxfDivCurl((Object)P->p[0]->field,
					NULL, &curl))
	    goto error;

	if (! curl)
	    goto error;

	if (! DXSetMember((Group)mg, NULL, curl))
	    goto error;
    }
    else
    {
	CompositeField cf = DXNewCompositeField();
	if (! cf)
	    goto error;

	for (i = 0; i < P->n; i++)
	{
	    Object attr;

	    attr = DXGetComponentAttribute(P->p[i]->field, "data", "dep");
	    if (! attr || DXGetObjectClass(attr) != CLASS_STRING)
	    {
		DXSetError(ERROR_INVALID_DATA, "data dependency attribute");
		DXDelete((Object)cf);
		goto error;
	    }

	    if (!strcmp(DXGetString((String)attr), "connections"))
	    {
		DXSetError(ERROR_INVALID_DATA, 
		    "cannot compute curl when data is dep on connections");
		DXDelete((Object)cf);
		goto error;
	    }

	    if (! DXSetMember((Group)cf, NULL, (Object)P->p[i]->field))
	    {
		DXDelete((Object)cf);
		goto error;
	    }
	}

	if (! _dxfDivCurl((Object)cf, NULL, &curl))
	    goto error;

	if (! curl)
	    goto error;

	if (! DXSetMember((Group)mg, NULL, curl))
	    goto error;
	
	DXDelete((Object)cf);

	if (! DXSetMember((Group)mg, NULL, curl))
	    goto error;
    }

    return OK;

error:
    return ERROR;
}