/*  pdnmesh - a 2D finite element solver
    Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>  
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 2 of the License, or
  (at your option) any later version.
  
  This program is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.
  
  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  $Id: bitree.c,v 1.13 2005/03/10 21:23:44 sarod Exp $
*/


/* something like a 2D search tree with Red Black Tree */
/* bi-tree ?? */
#include <stdio.h>
#include <stdlib.h> /* for malloc(), realloc() */
#include "types.h"

#define NILX &sentinelx /* leaves are sentinel nodes */
node_x sentinelx = { NILX, NILX, NILX, 0, BLACK,0};
#define NILY &sentinely
node_y sentinely = { 0,NILY, NILY, 0, BLACK, 0,0,0,0 };

/* initialized to all 0 */
Mesh M={NILX,0,0,0};
/* (re) initialise the mesh */
void 
BIT_init(Mesh *Ms, int max_points)
{
  if(Ms->maxpoints != 0) { 
   BIT_free(Ms);
  }
  Ms->count = 0;
  Ms->Xtree = NILX;
  Ms->maxpoints= max_points;
  Ms->Narray=(node_y**)calloc(Ms->maxpoints, sizeof(node_y*));
	if (!Ms->Narray) {
		fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
		exit(1);
  }

}


/* procedures neede to keep Red-Black tree balanced */
/* first rotation procedures */
/* note we do not use any global variables */
static void 
rotateLeftx( Mesh *Ms, node_x *x)
{
  /* rotate node x to left */
  node_x *y= x->right;
  /* establish x->right link */
  x->right = y->left;
  if ( y->left != NILX ) y->left->parent = x;

  /* establish y->parent link */
  if ( y != NILX ) y->parent = x->parent;
  if ( x->parent) {
    if ( x == x->parent->left)
      x->parent->left = y;
    else 
      x->parent->right = y;
  } else {
    Ms->Xtree = y;
  }

  /* link x and y */
  y->left = x;
  if ( x != NILX ) x->parent = y;
}

static void 
rotateRightx(  Mesh *Ms, node_x *x)
{
  /* rotate node x to right */
  node_x *y = x->left;

  /* establish x->left link */
  x->left = y->right;
  if ( y->right != NILX ) y->right->parent = x;

  /* establish y->parent link */
  if ( y != NILX ) y->parent = x->parent;
  if ( x->parent ) {
    if ( x == x->parent->right )
      x->parent->right = y;
    else
      x->parent->left = y;
  } else {
    Ms->Xtree = y;
  }

  /* link x and y */
  y->right = x;
  if ( x != NILX ) x->parent = y;
}

static void 
insertFixupx( Mesh *Ms, node_x *x )
{
  /* maintain red-black tree balance after insertion */
  /* check red-black tree properties */
  while ( x != Ms->Xtree && x->parent->colour == RED ) {
    /* we have a violation */
    if ( x->parent == x->parent->parent->left ) {
      node_x *y = x->parent->parent->right;
      if ( y->colour == RED ) {

	/* uncle is red */
	x->parent->colour = BLACK ;
	y->colour = BLACK;
	x->parent->parent->colour = RED;
	x = x->parent->parent;
      } else {

	/* uncle is black */
	if ( x == x->parent->right ) {
	  /* make x a left child */
	  x = x->parent;
	  rotateLeftx(Ms, x);
	}

	/* recolour and rotate */
	x->parent->colour = BLACK;
	x->parent->parent->colour = RED;
	rotateRightx(Ms, x->parent->parent );
      }
    } else {
      /* mirror image of above code */
      node_x *y = x->parent->parent->left;
      if ( y->colour == RED ) {

	/* unlce is red */
	x->parent->colour = BLACK;
	y->colour = BLACK;
	x->parent->parent->colour = RED;
	x = x->parent->parent;
      } else {

	/* uncle is black */
	if ( x == x->parent->left ) {
	  x = x->parent;
	  rotateRightx(Ms, x);
	}
	x->parent->colour = BLACK;
	x->parent->parent->colour = RED;
	rotateLeftx(Ms, x->parent->parent);
      }
    }
  }
  Ms->Xtree->colour = BLACK;
}

/* now a repetition of the above functions for y node types */
static void 
rotateLefty(node_y *x, node_x* xparent)
{
  /* rotate node x to left */
  node_y *y= x->right;
  /* establish x->right link */
  x->right = y->left;
  if ( y->left != NILY ) y->left->parent = x;
  
  /* establish y->parent link */
  if ( y != NILY ) y->parent = x->parent;
  if ( x->parent) {
    if ( x == x->parent->left)
      x->parent->left = y;
    else 
      x->parent->right = y;
  } else {
    xparent->Ytree=y;
  }
  
  /* link x and y */
  y->left = x;
  if ( x != NILY ) x->parent = y;
  
}

static void 
rotateRighty( node_y *x, node_x* xparent )
{
  /* rotate node x to right */
  node_y *y = x->left;
  
  /* establish x->left link */
  x->left = y->right;
  if ( y->right != NILY ) y->right->parent = x;
  
  /* establish y->parent link */
  if ( y != NILY ) y->parent = x->parent;
  if ( x->parent ) {
    if ( x == x->parent->right )
      x->parent->right = y;
    else
      x->parent->left = y;
  } else {
    xparent->Ytree=y;
  }
  
  /* link x and y */
  y->right = x;
  if ( x != NILY ) x->parent = y;
  
}

static void 
insertFixupy( node_y *x, node_x* xparent )
{
  /* maintain red-black tree balance after insertion */
  /* check red-black tree properties */
  while ( x != xparent->Ytree && x->parent->colour == RED ) {
    /* we have a violation */
    if ( x->parent == x->parent->parent->left ) {
      node_y *y = x->parent->parent->right;
      if ( y->colour == RED ) {
	
	/* uncle is red */
	x->parent->colour = BLACK ;
	y->colour = BLACK;
	x->parent->parent->colour = RED;
	x = x->parent->parent;
      } else {
	
	/* uncle is black */
	if ( x == x->parent->right ) {
	  /* make x a left child */
	  x = x->parent;
	  rotateLefty( x, xparent );
	}

	/* recolour and rotate */
	x->parent->colour = BLACK;
	x->parent->parent->colour = RED;
	rotateRighty( x->parent->parent, xparent );
      }
    } else {
	 /* mirror image of above code */
      node_y *y = x->parent->parent->left;
      if ( y->colour == RED ) {

	/* unlce is red */
	x->parent->colour = BLACK;
	y->colour = BLACK;
	x->parent->parent->colour = RED;
	x = x->parent->parent;
      } else {
	
	/* uncle is black */
	if ( x == x->parent->left ) {
	     x = x->parent;
	     rotateRighty( x, xparent );
	}
	x->parent->colour = BLACK;
	x->parent->parent->colour = RED;
	rotateLefty( x->parent->parent, xparent );
      }
    }
  }
  xparent->Ytree->colour = BLACK;
}

/* the following will insert the given point into the mesh if is already is */
/* not there, else it will do nothing but in both cases, the point number will be returned */
unsigned int 
BIT_insert(Mesh *Ms, point p )
{/* Mesh M is the global variable */
  node_x *currentx, *parentx, *xx;
  node_y *currenty, *parenty, *yy;
  
  currentx = Ms->Xtree;
  parentx=0;
     
  /* find the x coordinate first */
  while ( currentx != NILX ) {
    if ( p.x == currentx->x ) {
#ifdef DEBUG
      /* x cordinates match, find y coordinate */
     printf("x match\n"); 
#endif
      currenty = currentx->Ytree;
      parenty=0;
      while ( currenty != NILY ) {
	if ( p.y == currenty->y ) {
#ifdef DEBUG
	  printf("y match as well\n"); 
#endif
	  return( currenty->n );
	}
	parenty = currenty;
	currenty = ( p.y < currenty->y ) ? currenty->left : currenty->right;
      }
      /* if we reach here, we know x corrdinate exist but y does not */
      /* so add new y node to Ytree */
      /* set up new Y node */
      if ( ( yy =( node_y *) malloc( sizeof(node_y) )) == 0 ) {
	     fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	     exit(1);
      }
      yy->y = p.y;
      yy->z = p.z;
      yy->val = p.val;
      yy->parent = parenty;
      yy->left = NILY;
      yy->right = NILY;
      yy->colour = RED;
      yy->n = Ms->count; /* increment node count */
      Ms->Narray[Ms->count] = yy; /* point it to this */
      yy->xp = &(currentx->x); /* x coordinate */
      /* now insert it in y tree */
      if ( parenty ) {
	if ( p.y < parenty->y )
	  parenty->left = yy;
	else
	  parenty->right = yy;
      } else {
	currentx->Ytree = yy;
	/* certainly we will not reach this step */
	/* because there will be no node with just */
	/* an x value */
      }
#ifdef DEBUG
      printf("new y added\n"); 
#endif
      /* modify global variable */
      insertFixupy( yy, currentx );
 
      /* if root has changed in the function need to do the same here */
      /* if there was no change, just as fine because no change is done */
      if ( Ms->count >= Ms->maxpoints-1 ) {
#ifdef DEBUG	
	printf("point number exceeded, resizing\n");
#endif	
	/* now realloc the array to accomodate bigger size */
	Ms->maxpoints +=RESIZE; /* thousand more is enough */
        Ms->Narray=(node_y**)realloc((void*)Ms->Narray,Ms->maxpoints*sizeof(node_y*));
	if ( Ms->Narray == 0 ) {
	  fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	  exit(2);
	}
#ifdef DEBUG
	printf("new size %d\n",Ms->maxpoints);
#endif      
      }
      /* now increment count */
      Ms->count++;
      return(yy->n);
    }
    parentx = currentx;
    currentx=( p.x < currentx->x )? currentx->left: currentx->right;
  }
  /* now we know that x coordinate is not in the tree, and */
  /* obviously the y coordinate */
  /* set up new x node */
  if ( (xx =(node_x*) malloc( sizeof(node_x))) == 0 ){
		fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
    exit(1);
  }
  xx->x = p.x;
  xx->parent= parentx;
  xx->left = NILX;
  xx->right= NILX;
  xx->colour = RED;
  /* now insert y node to this x node */
  if ( ( yy =(node_y*) malloc( sizeof(node_y) )) == 0 ) {
		fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
    exit(1);
  }
  yy->y = p.y;
  yy->z = p.z;
  yy->val = p.val;
  yy->parent = 0; /* spacial case--should be null-- */
  yy->left = NILY;
  yy->right = NILY;
  yy->colour = RED; /* can make this BLACK and omit calling insertFixupy() */
  yy->n = Ms->count; /* increment node count */
  Ms->Narray[Ms->count] = yy; /* point it to this */
  /* we have a limit to the number of points */
  /* check this before incrementing */
  /* do it at the end */
  yy->xp = &(xx->x); /* x coordinate */
  xx->Ytree = yy;
#ifdef DEBUG
  printf("new x and y added\n"); 
#endif
  /* now insert xx in the X tree */
  if ( parentx ) {
    if ( p.x < parentx->x )
      parentx->left = xx;
    else 
      parentx->right = xx;
  } else {
    Ms->Xtree = xx; /* the first point */
  }
  insertFixupx(Ms, xx ); /* we pass globally */
  insertFixupy( xx->Ytree, xx); 
   /* it took me about two weeks to figure out to include the above */
   /* insertFixupy() for the tree to properly work :)               */
   /* this teaches us a lesson: the bugs in a program are NOT where */
   /* you  expect them to be */
   /* note we need not call this function either, we can modify line 364 */
  if ( Ms->count >= Ms->maxpoints-1 ) {  
#ifdef DEBUG
    printf("point number exceeded, resizing\n");
#endif    
    /* now realloc the array to accomodate bigger size */
    Ms->maxpoints +=RESIZE; /* thousand more is enough */
    Ms->Narray=(node_y**)realloc((void*)Ms->Narray,Ms->maxpoints*sizeof(node_y*));
    if (!Ms->Narray) {
			fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
      exit(2);
    }
#ifdef DEBUG
    printf("new size %d\n",Ms->maxpoints);
#endif 
  }
  Ms->count++;
 
  return( Ms->count-1 );

}

/* initialise the mesh */
static void
destroy_Xtree( node_x* xx )
{
if ( xx->right != NILX )
  destroy_Xtree( xx->right );
if ( xx->left != NILX )
  destroy_Xtree( xx->left );
free( xx);
}


/* free the mesh */
void 
BIT_free(Mesh *Ms)
{ 
  /* free up memory */
 unsigned int i;
 for ( i =0; i < Ms->count; i++) {
      if (Ms->Narray[i] != 0) {
#ifdef DEBUG
 printf("bit_free %d "MDF" "MDF" "MDF"\n",i,*((Ms->Narray[i])->xp),
		 (Ms->Narray[i])->y,*(Ms->Narray[i])->z);
#endif
   if ((Ms->Narray[i])->z != 0) free((Ms->Narray[i])->z); /* z array */
    free( Ms->Narray[i] ); /* free up the y nodes */
	  }
 }
 free( Ms->Narray );
 Ms->count=0;
 Ms->maxpoints=0;

 if ((Ms->Xtree != NILX) && ( Ms->Xtree != 0) ) destroy_Xtree( Ms->Xtree );
}