path: root/src/GetLocalLine.c

#include <stdio.h>
#include <stdlib.h>
#include <math.h>  

#include "cctk.h"
#include "CactusPUGH/PUGH/src/include/pugh.h"

static char *rcsid = "$Header$";

CCTK_FILEVERSION(CactusPUGH_PUGHSlab_GetLocalLine_c)

#define MAX_DIM 3
#define MAX_FACE 6
#define BAD -42
#define EPS 1e-10

#define ABS(a) ((a)<0 ? -(a) : (a))
#define SIGN(a) ((a)<0.0 ? (-1.0) : (1.0))
#define LIN_INDEX3D(extras,P) ((P[0])+ (extras)->nsize[0]*((P[1]) + (extras)->nsize[1]*(P[2])))

#define MIN(a,b) ((a)<(b) ? (a) : (b))
#define MAX(a,b) ((a)>(b) ? (a) : (b))

/*$#define LINE_DEBUG$*/

int IsPointOnProc_INT (int dim, int *gownership, const int *pnt);
int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt);

int IsPointOnProc_INT(int dim, int *gownership, const int *pnt)
{
  int retval=1; /*assume YES */
  int idim;
  for (idim=0;idim<dim;idim++)
  {
    /* out of range -> set retval to zero;*/
    if ((gownership[2*idim  ]>pnt[idim]) || 
	(gownership[2*idim+1]<pnt[idim]))
      retval&=0;  /* incremental AND on the return value */
  }
  return(retval);
}

int IsPointOnProc_REAL(int dim, int *gownership, CCTK_REAL *pnt)
{
  int retval=1; /*assume YES */
  int idim;
  for (idim=0;idim<dim;idim++)
  {
    /* out of range -> set retval to zero;*/
    if ((gownership[2*idim  ]>pnt[idim]) || 
	(gownership[2*idim+1]<pnt[idim]))
      retval&=0;  /* incremental AND on the return value */
  }
  return(retval);
}

int GetLocalLine(cGH *GH, pGExtras *extras,
		 int vardim, const int *S_l, const int *u_l, int ds, int len_l, 
		 int *locstart, int *locend, int *nlocpoints)
{

  const int n_s[MAX_DIM][MAX_FACE] = 
  {{ -1, 1, 0, 0, 0, 0},
   { 0, 0,-1, 1, 0, 0},
   { 0, 0, 0, 0,-1, 1}};


  /* variables for proc. surfaces */
  CCTK_REAL lmd_s[MAX_FACE], lmd[2], t_lmd;
  int P_s[MAX_DIM][MAX_FACE];
  int gownership[MAX_FACE];

  int notcoll[MAX_FACE];
  CCTK_REAL t_pnt1[MAX_DIM];

  int prochasline=0;

  int tmp1,tmp2;
  CCTK_REAL rtmp1;

  int idim, iface, linsec;
  int myproc;



#ifdef LINE_DEBUG
  printf("\n+++++++++ GET LOCAL LINE \n");
  printf("S_l [XYZ] : %d %d %d \n",
	 S_l[0], S_l[1], S_l[2]);
  printf("u_l [XYZ] : %d %d %d    ds %d \n",
	 u_l[0], u_l[1], u_l[2], ds);

  printf("n_s: face 0/1: [%d %d %d], [%d %d %d] \n",
	 n_s[0][0],n_s[1][0],n_s[2][0],
	 n_s[0][1],n_s[1][1],n_s[2][1]);
  printf("n_s: face 2/3: [%d %d %d], [%d %d %d] \n",
	 n_s[0][2],n_s[1][2],n_s[2][2],
	 n_s[0][3],n_s[1][3],n_s[2][3]);
  printf("n_s: face 4/5: [%d %d %d], [%d %d %d] \n",
	 n_s[0][4],n_s[1][4],n_s[2][4],
	 n_s[0][5],n_s[1][5],n_s[2][5]);
#endif

  myproc = CCTK_MyProc (GH);

  for (idim=0;idim<vardim;idim++)
  {
    /* GLOBAL OWNERSHIP: lower/upperbnd are INCLUSIVE */
    gownership[2*idim]  = extras->lb[myproc][idim];
    if (GH->cctk_bbox[2*idim]==0) 
      gownership[2*idim]  += extras->nghostzones[idim];

    gownership[2*idim+1]= extras->ub[myproc][idim];
    if (GH->cctk_bbox[2*idim+1]==0) 
      gownership[2*idim+1]-= (extras->nghostzones[idim]);    
  }

  
  for (idim=0;idim<vardim;idim++)
  {
    P_s[idim][0] = gownership[2*idim];    /* Aufpunkt, minus faces */
    P_s[idim][2] = gownership[2*idim];    /* Aufpunkt, minus faces */
    P_s[idim][4] = gownership[2*idim];    /* Aufpunkt, minus faces */
    P_s[idim][1] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
    P_s[idim][3] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
    P_s[idim][5] = gownership[2*idim+1];  /* Aufpunkt, plus faces */
  }
  
#ifdef LINE_DEBUG
  printf("P_s (face 0): %d %d %d \n", 
	 P_s[0][0],P_s[1][0],P_s[2][0]);
  printf("P_s (face 1): %d %d %d \n\n",
	 P_s[0][1],P_s[1][1],P_s[2][1]);
#endif

  /* Lamda for each surface */
  for (iface=0;iface<2*vardim;iface++)
  {
    tmp1=0.0; tmp2=0.0;

    for (idim=0;idim<vardim;idim++) 
    {
      tmp1+= (S_l[idim]-P_s[idim][iface])*n_s[idim][iface];
      tmp2+= (ds*u_l[idim]*n_s[idim][iface]);
    }
    if (tmp2==0) 
    {
      lmd_s[iface]=BAD;
      notcoll[iface]=0;
    }
    else
    {
      notcoll[iface]=1;
      rtmp1         = (CCTK_REAL)(-1.0)*(tmp1)/(tmp2);
      if (rtmp1<0.0) lmd_s[iface]=BAD;
      else           lmd_s[iface]=rtmp1;
#ifdef LINE_DEBUG
      printf("lmd_s (face%d) %f \n",iface,rtmp1);
#endif
    }
  }

   /* Find 2 valid intersection with the proc surfaces */
  /* tricky: duplicate points corners and the startpoint */
  /* We have to restrict the line to positive lambda, otherwise
     it is getting extremly ugly in finding out the local
     points */

  lmd[0] = BAD;
  lmd[1] = BAD;
  linsec = BAD;

#ifdef LINE_DEBUG
  printf("testing for local Startpt: %d %d %d",
	 S_l[0],S_l[1],S_l[2]);
#endif

  if (IsPointOnProc_INT(vardim,gownership,S_l)==1)
  {  
#ifdef LINE_DEBUG
    printf(" ... YES ... valid first\n");
#endif
    lmd[0]      = 0;
    linsec      = LIN_INDEX3D(extras,S_l);
  }
  else
  {
#ifdef LINE_DEBUG
    printf(" ... NO \n");
#endif
  }
  

  for (iface=0;iface<MAX_FACE; iface++)
  {
#ifdef LINE_DEBUG
    printf("checking face %d: coll: %d  lmd_s %d ", 
	   iface, notcoll[iface], lmd_s[iface]);
#endif
    if (notcoll[iface]==1 && lmd_s[iface]!=BAD) 
    {
      for (idim=0;idim<vardim;idim++)
	t_pnt1[idim]=S_l[idim]+lmd_s[iface]*ds*u_l[idim];

#ifdef LINE_DEBUG
      printf("  testing for local intersec. (face %d lamda %+4.2f):  %4.2f %4.2f %4.2f\n",
	     iface, lmd_s[iface], t_pnt1[0],t_pnt1[1],t_pnt1[2]);
#endif
      
      /* check if point is local */
      if (IsPointOnProc_REAL(vardim,gownership,t_pnt1)==1)
      {
#ifdef LINE_DEBUG
	printf("    ... YES ... ");
#endif

	/* CASE1: first valid face: this point is good */
	if (linsec<0)
	{
#ifdef LINE_DEBUG
	  printf(" valid first \n");
#endif
	  lmd[0]      = lmd_s[iface];
	  linsec      = LIN_INDEX3D(extras,t_pnt1);
	}            
	/* CASE2: we already have intersec. 
	   check for duplicates: eg. a diagonal through the n-space, with 
	   one exception: make sure that we allow the endpoint to be located 
	   on the startpoint if lmd == 0 (which is the case if we start at the last
	   point of a proc) */
	else if (linsec==LIN_INDEX3D(extras,t_pnt1) && lmd_s[iface]>EPS)
	{
#ifdef LINE_DEBUG
	  printf(" but duplicate, skip \n");
#endif
	}
	else
	{
#ifdef LINE_DEBUG
	  printf(" valid second \n");
#endif
	  lmd[1]      =lmd_s[iface];
	}
      }
      else
      {
#ifdef LINE_DEBUG
	printf("     ... NO \n");
#endif
      }
    }
    else
    {
#ifdef LINE_DEBUG
    printf(" ... skipping face: %d %4.2f \n",iface,lmd_s[iface]);
#endif
    }
  }
 
  /* flag if a proc has intersections, order the intersections 
     if we have any  */
  if (lmd[0]!=BAD && lmd[1]!=BAD)
  {   
    prochasline=1;

    /* switch: lmd[0] has to be the closer one to S_l */
    if (lmd[0]>lmd[1])
    {
      t_lmd  = lmd[0];
      lmd[0] = lmd[1];
      lmd[1] = t_lmd;
    }
    lmd[0]=ceil (lmd[0]);
    lmd[1]=floor(lmd[1]);
#ifdef LINE_DEBUG
    printf("LOCAL INTERSECTIONS: (ordered in lmd) %4.2f %4.2f \n",lmd[0],lmd[1]);
#endif
  }
  else if (lmd[0]==BAD && lmd[1]==BAD)
  {   
    prochasline = 0;
#ifdef LINE_DEBUG
    printf("NO INTERSECTION ON THIS PROC\n");
#endif
  }
  else if ((lmd[0]==BAD && lmd[1]!=BAD) ||
	   (lmd[1]==BAD && lmd[0]!=BAD))
  {
    printf("ERROR: only one local intersection: you are screwed. lmd 0/1  %4.2f %4.2f\n",
	   lmd[0],lmd[1]);
  }
  
  /* patch if lambdas are greater specified line length */
  if (len_l>=0) 
  {
    if (len_l>lmd[0] && len_l<lmd[1])
      lmd[1]=floor(len_l);
    else if (len_l<lmd[0])
    {
      lmd[0]=BAD;
      lmd[1]=BAD;
      prochasline=0;
    }
  }

  /* Calculate the intersections points in global index space*/
  if (prochasline) 
  {
    nlocpoints[0]= ABS(lmd[0]-lmd[1])+1;
    for (idim=0;idim<vardim;idim++)
    {
      locstart[idim] = (S_l[idim]+lmd[0]*ds*u_l[idim]);
      locend[idim]   = (S_l[idim]+lmd[1]*ds*u_l[idim]);
    }
  } 
  else
  {
    nlocpoints[0]  = 0;
    for (idim=0;idim<vardim;idim++)
    {
      locstart[idim] = BAD;
      locend[idim]   = BAD;  
    }
  }

#ifdef LINE_DEBUG
  printf("\n>>>>> GLOBAL START/END: (%d %d %d) -> (%d %d %d)  = %d\n",
	 locstart[0],locstart[1],locstart[2],
	 locend[0],locend[1],locend[2],nlocpoints[0]);
#endif

  return(0);
  
}