/*
* $Id: hydram.c,v 1.1.1.1 2005/09/18 22:05:49 dhmunro Exp $
* workers to convert hydra bnd specs into HX_boundary arrays
*/
/* Copyright (c) 2005, The Regents of the University of California.
* All rights reserved.
* This file is part of yorick (http://yorick.sourceforge.net).
* Read the accompanying LICENSE file for details.
*/
#include "hex.h"
#include "hydram.h"
/* ------------------------------------------------------------------------ */
static int get_2dims(int pn, long blo[], long sdims[]);
static long set_scratch(long n, long *ndx, long bnd[], long blo[],
long *scratch, long *sdims);
static int orient_compute(long bnds[], long bndr[], long *p2dp, long *p2dq,
long pq0, long *p2p, long *p2q);
long hydra_blks(long nblks, long (*blo)[4])
{
long offset= 0;
long lscratch= 0;
long b, ni, nj, nk, nn;
for (b=0 ; b<nblks ; b++) {
blo[b][0]= offset;
ni= blo[b][1];
nj= blo[b][2];
nk= blo[b][3];
offset+= (blo[b][3]*= (blo[b][2]*= ni));
if (ni<nj) nn= nk<ni? ni*nj : nj*nk;
else nn= nk<nj? ni*nj : nk*ni;
if (nn>lscratch) lscratch= nn;
}
return lscratch;
}
long hydra_mrk(int ibnd, long (*tbound)[3], long *blo, long *bnd,
long n, long *ndx)
{
int pn= (int)bnd[0];
long jp= bnd[1];
long kp= bnd[2];
long jq= blo[1]; /* always jp-3?? */
long kq= blo[2]; /* always kp-3*jp?? */
long i, px0, px1, nijk[3];
long ijk= -1;
int p0, p1;
if (pn<0) pn= -pn;
pn--;
p0= !pn;
p1= pn^p0^3;
nijk[0]= jq;
if (pn!=1) nijk[1]= kq/jq;
if (pn!=2) nijk[2]= blo[3]/kq;
px0= nijk[p0];
px1= nijk[p1];
for (i=0 ; i<n ; i++) {
/* use hydra in-memory strides from bnd[] to extract separate
* i, j, and k indices from ndx[], then adjust these by subtracting
* the 2 ghost planes hydra uses
* -- could make -2 a parameter to this function... */
nijk[0]= ndx[i]%jp - 2;
nijk[1]= (ndx[i]%kp)/jp - 2;
nijk[2]= ndx[i]/kp - 2;
if (nijk[pn]<0) return -13;
if (ijk!=nijk[pn]) {
/* all the ndx[] must lie in a single plane of constant pn */
if (i) return -1;
ijk= nijk[pn];
}
if (nijk[p0]<0 || nijk[p1]<0) {
/* hydra may undershoot BC by one node */
if (nijk[p0]<-1 || nijk[p1]<-1) return -11;
continue;
}
if (nijk[p0]>=px0 || nijk[p1]>=px1) {
/* hydra may overshoot BC by one node */
if (nijk[p0]>px0 || nijk[p1]>px1) return -12;
continue;
}
tbound[nijk[0]+nijk[1]*jq+nijk[2]*kq][pn]= ibnd;
}
return ijk; /* index of plane containing all of ndx[] */
}
long hydra_adj(long (*bound)[3], long (*tbound)[3],
long *blo, long n, long *tcheck)
{
long start= -1;
long strides[4];
long i, j, k, is, js, ix, jx, face;
int pn, p0, p1;
strides[0]= 1;
strides[1]= blo[1];
strides[2]= blo[2];
strides[3]= blo[3];
for (bound+=blo[0] ; n-- > 0 ; tcheck+=2) {
pn= tcheck[0];
if (!pn) continue;
face= (pn>0);
if (!face) pn= -pn;
pn--;
p0= !pn;
p1= pn^p0^3;
if (start<0) {
face|= (pn<<1);
if (!(face&1)) face+= 6*strides[pn];
}
k= tcheck[1]*strides[pn];
is= strides[p0];
ix= strides[p0+1];
js= strides[p1];
jx= k+strides[p1+1];
for (j=k+js ; j<jx ; j+=js) {
for (i=0 ; i<ix-is ;) {
for ( ; i<ix-is ; i+=is)
if (tbound[i+j-js][pn] && tbound[i+j][pn]) break;
for (i+=is ; i<ix ; i+=is) {
if (!tbound[i+j-js][pn] || !tbound[i+j][pn]) break;
bound[i+j][pn]= tbound[i+j][pn];
if (start<0 && tbound[i+j][pn]==-1)
start= 6*(i+j+blo[0]) + face;
}
}
}
}
return start;
}
long hydra_bnd(long ibnd, long (*bound)[3], long *scratch,
long *blos, long *blor, long *bnds, long *bndr,
long n, long *ndxs,
long *ndxr, HX_blkbnd *mbnds, long rblock)
{
long i, j, ii, jj, sdims[2];
/* determine the face strides and orientations */
int pns= get_2dims((int)bnds[0], blos, sdims);
long js= sdims[0]; /* scratch 2nd dimension stride */
long ks= sdims[1]; /* scratch size */
long ib= pns? 1 : blos[1];
long jb= (pns==2)? blos[1] : blos[2];
long start;
/* mark scratch space with nodes in ndxs */
long splane= set_scratch(n, ndxs, bnds, blos, scratch, sdims);
if (splane<0) return splane;
/* scan to find first point */
for (j=0 ; j<ks-js-1 ; j++)
if (scratch[j] && scratch[j+1] && scratch[j+js] && scratch[j+js+1]) break;
/* first quad face is j, j+1, j+js, j+js+1 */
if (j>=ks-js) return -2;
/* offset bound to the proper plane of the proper block
* note that in the pns direction, the index is a nodal index, while
* in the other two directions, the mark will be at a cell index */
start= splane*(pns?blos[pns]:1) + blos[0] + ib + jb;
bound+= start;
i= j%js;
j-= i;
ii= i*ib;
jj= (j/js)*jb;
if (!ndxr) {
/* THIS DOES NOT WORK AS INTENDED
* the hydra_mrk and hydra_adj routines must be
* used instead to mark BCs */
#ifdef NEVER_DEFINE_THIS
for (; j<ks-js ; i=ii=0,j+=js,jj+=jb) { /* loop on rows */
do { /* loop on segments in this row */
for (; i<js-1 ; i++,ii+=ib) /* scan to first face in segment */
if (scratch[j+i] && scratch[j+i+js]) break;
for (i++ ; i<js ; i++,ii+=ib) { /* scan to end of segment */
if (!scratch[j+i] || !scratch[j+i+js]) break;
bound[jj+ii][pns]= ibnd; /* mark the face */
}
} while (i<js);
}
/* need to return face index (6*cell+face) of first face in this bc */
if (bnds[0]>0) {
pns= (pns<<1) | 1; /* 0-5 face index, no cell adjustment */
} else {
start+= (pns?blos[pns]:1); /* cell index */
pns<<= 1; /* 0-5 face index */
}
return 6*start + pns;
#endif
return -4;
} else {
/* bnds, bndr, and the first quad determine
* the relative orientation of the s and r blocks
* - note that scratch markers are indices+1 into ndxs, ndxr */
long pq0= ndxr[scratch[i+j]-1];
long dp= ndxr[scratch[i+j+1]-1] - pq0;
long dq= ndxr[scratch[i+j+js]-1] - pq0;
long p= blor[1];
long q= blor[2];
long jsmi = js - i;
int orient= orient_compute(bnds, bndr, &dp, &dq, pq0, &p, &q);
if (orient<0) return -3;
q+= blor[0]+p; /* q is global index of first cell in r block */
j += i; /* like q, combine initial i into j ... */
jj += ii; /* ... and ii into jj */
for (i=ii=p=0 ; j<ks-jsmi ; i=ii=p=0,j+=js,jj+=jb,q+=dq) {
do {
for (; i<jsmi-1 ; i++,ii+=ib,p+=dp)
if (scratch[j+i] && scratch[j+i+js]) break;
for (i++ ; i<jsmi ; i++,ii+=ib,p+=dp) {
if (!scratch[j+i] || !scratch[j+i+js]) break;
bound[jj+ii][pns]= ibnd+1;
mbnds[ibnd].block= rblock;
mbnds[ibnd].cell= p+q;
mbnds[ibnd].orient= orient;
ibnd++;
}
} while (i<jsmi);
}
return ibnd;
}
}
static int get_2dims(int pn, long blo[], long sdims[])
{
long jq= blo[1];
long kq= blo[2];
if (pn<0) pn= -pn;
pn--;
if (pn==0) {
sdims[0]= kq/jq;
sdims[1]= blo[3]/jq;
} else if (pn==1) {
sdims[0]= jq;
sdims[1]= jq*(blo[3]/kq);
} else {
sdims[0]= jq;
sdims[1]= kq;
}
return pn;
}
static long set_scratch(long n, long *ndx, long bnd[], long blo[],
long *scratch, long *sdims)
{
int pn= (int)bnd[0];
long jp= bnd[1];
long kp= bnd[2];
long jq= blo[1]; /* always jp-3?? */
long kq= blo[2]; /* always kp-3*jp?? */
long j0= sdims[0];
long ns= sdims[1];
long i, px0, px1, nijk[3];
long ijk= -1;
int p0, p1;
if (pn<0) pn= -pn;
pn--;
p0= !pn;
p1= pn^p0^3;
nijk[0]= jq;
if (pn!=1) nijk[1]= kq/jq;
if (pn!=2) nijk[2]= blo[3]/kq;
px0= nijk[p0];
px1= nijk[p1];
for (i=0 ; i<ns ; i++) scratch[i]= 0;
for (i=0 ; i<n ; i++) {
/* use hydra in-memory strides from bnd[] to extract separate
* i, j, and k indices from ndx[], then adjust these by subtracting
* the 2 ghost planes hydra uses
* -- could make -2 a parameter to this function... */
nijk[0]= ndx[i]%jp - 2;
nijk[1]= (ndx[i]%kp)/jp - 2;
nijk[2]= ndx[i]/kp - 2;
if (nijk[pn]<0) return -13;
if (ijk!=nijk[pn]) {
/* all the ndx[] must lie in a single plane of constant pn */
if (i) return -1;
ijk= nijk[pn];
}
if (nijk[p0]<0 || nijk[p1]<0) {
/* hydra may undershoot BC by one node */
if (nijk[p0]<-1 || nijk[p1]<-1) return -11;
continue;
}
if (nijk[p0]>=px0 || nijk[p1]>=px1) {
/* hydra may overshoot BC by one node */
if (nijk[p0]>px0 || nijk[p1]>px1) return -12;
continue;
}
scratch[nijk[p0]+nijk[p1]*j0]= i+1;
}
return ijk; /* index of plane containing all of ndx[] */
}
static int orient_compute(long bnds[], long bndr[], long *p2dp, long *p2dq,
long pq0, long *p2p, long *p2q)
{
long i, j, k;
int f0, f2;
long dp= *p2dp;
long dq= *p2dq;
int sns= bnds[0]<0;
int pns= (int)(sns? -bnds[0] : bnds[0]) - 1;
/* int pps= !pns; */
/* int pqs= pns^pps^3; */
int snr= bndr[0]<0;
int pnr= (int)(snr? -bndr[0] : bndr[0]) - 1;
int ppr, pqr;
/* direction pnr in r block corresponds to direction pns in s block,
* or opposite to pns if sns==snr
* dp is r block stride corresponding to + step in direction pps
* dq is r block stride corresponding to + step in direction pqs
*/
int sp= dp<0;
int sq= dq<0;
int sr= sns==snr;
if (sp) dp= -dp;
if (sq) dq= -dq;
/* check consistency of r block orientation */
if (pnr) {
if (dp==1) {
ppr= 0;
pqr= 3^pnr;
if (dq!=bndr[pqr]) return -1;
} else {
ppr= 3^pnr;
pqr= 0;
if (dq!=1 || dp!=bndr[ppr]) return -1;
}
} else {
if (dp==bndr[1]) {
ppr= 1;
pqr= 2;
if (dq!=bndr[pqr]) return -1;
} else {
ppr= 2;
pqr= 1;
if (dp!=bndr[ppr] || dq!=bndr[pqr]) return -1;
}
}
if ((pns==1) ^ (((pnr==2)?0:pnr+1)!=ppr) ^ sp^sq^sr) return -1;
/* find which faces in r correspond to faces 0 and 2 in s */
if (pns) {
/* face 0s is +- dp direction (pps==0 if pns!=0) */
f0= (ppr<<1) | sp;
if (pns==1) {
/* face 2s is +- bndr[0] direction */
f2= (pnr<<1) | sr;
} else {
/* face 2s is +- dq direction */
f2= (pqr<<1) | sq;
}
} else {
/* face 0s is +- bndr[0] direction, face 2s is +-dp direction */
f0= (pnr<<1) | sr;
f2= (ppr<<1) | sp;
}
/* compute separate (p,q) for first cell
* corners are pq0, pq0+dp, pq0+dq, pq0+dp+dq
* if bndr[0]>0, the boundary is a min face,
* so the nodal index pq0 must be shifted to become a cell index
* similar shifts must be applied in the other two directions
* after all this, pq0 is the cell index in a mesh with ghosts,
* so the next step is to split it into its ijk and remove ghosts
* note that *p2p=blor[1]=stride for j, *p2q=blor[2]=stride for k
* on input
*/
if (snr) pq0+= pnr? bndr[pnr] : 1;
if (!sp) pq0+= dp;
if (!sq) pq0+= dq;
i= pq0%bndr[1] - 2;
j= (pq0%bndr[2])/bndr[1] - 2;
k= pq0/bndr[2] - 2;
if (i<0 || j<0 || k<0) return -1;
j*= (*p2p);
k*= (*p2q);
if (ppr==0) {
*p2dp= sp? -1 : 1;
*p2dq= (pqr==1)? (*p2p) : (*p2q);
if (sq) *p2dq= -(*p2dq);
*p2p= i;
*p2q= j+k;
} else if (ppr==1) {
*p2dp= sp? -(*p2p) : (*p2p);
*p2dq= (pqr==0)? 1 : (*p2q);
if (sq) *p2dq= -(*p2dq);
*p2p= j;
*p2q= k+i;
} else {
*p2dp= sp? -(*p2q) : (*p2q);
*p2dq= (pqr==0)? 1 : (*p2p);
if (sq) *p2dq= -(*p2dq);
*p2p= k;
*p2q= i+j;
}
/* return orientation */
f2^= ((f0&4)? f0-4 : f0+2);
if (f2&4) f2^= 6;
return (f0<<2) | f2;
}
/* ------------------------------------------------------------------------ */
#ifdef DEGHOST_ADJUSTS_NDX
static long deghost(long n, long *ndx, long bnd[], long blo[], long *sdims)
{
int pn= (int)bnd[0];
long jp= bnd[1];
long kp= bnd[2];
long jq= blo[1]; /* always jp-3?? */
long kq= blo[2]; /* always kp-3*jp?? */
long i, nijk[3];
long ijk= -1;
int p0, p1;
if (pn<0) pn= -pn;
pn--;
p0= !pn;
p1= pn^p0^3;
for (i=0 ; i<n ; i++) {
/* use hydra in-memory strides from bnd[] to extract separate
* i, j, and k indices from ndx[], then adjust these by subtracting
* the 2 ghost planes hydra uses
* -- could make -2 a parameter to this function... */
nijk[0]= ndx[i]%jp - 2;
nijk[1]= (ndx[i]%kp)/jp - 2;
nijk[2]= ndx[i]/kp - 2;
if (nijk[0]<0 || nijk[1]<0 || nijk[2]<0) return -1;
/* all the ndx[] must lie in a single plane of constant abs(pn) */
if (ijk!=nijk[pn]) {
if (i) return -1;
ijk= nijk[pn];
}
/* reconstruct ndx[] adjusted to index into arrays without
* any ghost planes, as in hex */
ndx[i]= nijk[0] + jq*nijk[1] + kq*nijk[2];
}
return ijk; /* index of plane containing all of ndx[] */
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1