/*  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: pdnmesh.c,v 1.57 2005/04/24 05:32:24 sarod Exp $
*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#ifndef WIN32
#include <unistd.h>
#include <gtk/gtk.h>
#include <gdk/gdkkeysyms.h>
#include <gtk/gtkgl.h>
#endif/*WIN32 */
#ifdef WIN32
#include <windows.h>
#endif/*!WIN32*/
#include <GL/gl.h>
#include <GL/glu.h>

#include "types.h"

int ntriangles,mouse_responce_flag,windowid;
DAG_tree dt;
RBT_tree rt;
double zoom_x1,zoom_x2,zoom_y1,zoom_y2;
GLdouble current_mat[16]; /* current zoom transform */
/* for getopt() */
extern char *optarg;
extern int optind, opterr, optopt;


int degree_of_freedom=1; /* degree of freedom for a point */
int requested_degree_of_freedom=0; /* user request, for eigenvalue problems */
int max_iteration_limit; /* limit for iterative solution */
int text_mode_only; /* flag to detect text mode operation */

/* boundary polygons */
extern boundary bound;

/* input cord filename */
char *  input_cord_filename;

/* vector cross product */
#define CROSS_PRODUCT(x1,y1,x2,y2) \
  ((x1)*(y2)-(y1)*(x2))
/* vector dot product */
#define DOT_PRODUCT(x1,y1,x2,y2) \
  ((x1)*(x2)+(y1)*(y2))
/* function to test point p is on a line segment (p1-p2) */
static int
is_point_online_p(point p, point p1,point p2)
{
  /* p = p1 + lambda p2/(1+lambda) */
  double lambda;

  /*if (((p.x == p1.x)&&(p.y==p1.y))||((p.x==p2.x)&&(p.y==p2.y))) { */
  if ((IS_ZERO(p.x - p1.x)&&IS_ZERO(p.y-p1.y))||(IS_ZERO(p.x-p2.x)&&IS_ZERO(p.y-p2.y))) {
    /* points coincide==>on line */
    return(1);
  }
  /* horizontal line */
  if (IS_ZERO(p.x-p1.x)&&IS_ZERO(p.x-p2.x)) {
    return((p.y-p1.y)/(p2.y-p.y) >0.0f);
  }
  /* vertical line */
  if ( IS_ZERO(p.y-p1.y)&&IS_ZERO(p.y-p2.y) ) {
    return((p.x-p1.x)/(p2.x-p.x) >0.0f);
  }
  lambda=(p.x-p1.x)/(p2.x-p.x);
  if ( IS_ZERO(lambda - (p.y-p1.y)/(p2.y-p.y)) ) {
    return(lambda>0.0f);
  }
  /* else  */
  return(0);
}

/* function to check point p is the same side as point c of line a-b */
/* return 1 if true, 0 if false, -1 if floating point over/under flow */
static int
same_side(point p, point c, point a, point b) {
 double cp1,cp2;
 cp1=CROSS_PRODUCT(b.x-a.x,b.y-a.y,p.x-a.x,p.y-a.y);
 cp2=CROSS_PRODUCT(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
#ifdef DEBUG
 printf("[cross %lf,%lf",cp1,cp2);
#endif
 /* else */
 /*(cp1*cp2 >= 0) */
 if ((cp1==0.0f)||(cp2==0.0f)
   ||((cp1>0.0f)&&(cp2>0.0f))
   ||((cp1<0.0f)&&(cp2<0.0f))) {
#ifdef DEBUG
  printf("same side]");
#endif
  return(1);
 } else {
#ifdef DEBUG
  printf("NOT same side]");
#endif
  return(0);
 }
}

static int 
inclusion(const void *trig,const void *pt, const void *MM/* Mesh*/){
  triangle *t;
  point *p,a,b,c;
  Mesh *G;
  t=(triangle *)trig;
  p=(point *)pt;
  G=(Mesh*)MM;
  a.x=Mx(t->p0,*G);
  a.y=My(t->p0,*G);
  b.x=Mx(t->p1,*G);
  b.y=My(t->p1,*G);
  c.x=Mx(t->p2,*G);
  c.y=My(t->p2,*G);
  if (is_point_online_p(*p,a,b)
    ||is_point_online_p(*p,b,c)
    ||is_point_online_p(*p,c,a)) {
#ifdef DEBUG
  printf("inclusion: point on boundary\n");
#endif
  return(2);
  }

 /* if point is on line any triangle edge,
    we consider that to be included */
#ifdef DEBUG
 printf("inclusion (%lf,%lf)-(%lf,%lf)-(%lf,%lf) in (%lf,%lf)\n",
     a.x,a.y,b.x,b.y,c.x,c.y,p->x,p->y);
#endif
 /* check for errors */
 if(same_side(*p,c,a,b)
 &&same_side(*p,a,b,c)
 &&same_side(*p,b,c,a)) {
 /* we have floating point over/under flow
   ans switch to alternative method */
#ifdef DEBUG
  printf("inclusion: triangle %d includes the point\n",t->n);
#endif
  return(1);
 }
 /*  else */ 
#ifdef DEBUG
  printf("inclusion: triangle %d NOT include the point\n",t->n);
#endif
  return(0);
}


#ifdef EXTRA_SOURCE
static double
determ( double xx1, double yy1, double xx2, double yy2, double xx3, double yy3 ) {
  return ( xx1*(yy2-yy3)+xx2*(yy3-yy1)+xx3*(yy1-yy2) );
}

static int
ccw(point p0,point p1,point p2)
{
  return(determ(p0.x,p0.y,p1.x,p1.y,p2.x,p2.y)>=0.0f);
}
/* is d inside triangle (a,b,c) ? */
#define INSIDE(a,b,c,d) \
  (ccw(a,b,d)&&ccw(b,c,d)&&ccw(c,a,d))

/* checking inclusion of point in triangle */
static int 
inclusion1(const void *trig,const void *pt, const void *MM)
{
  triangle *t;
  point *p,a,b,c;
  Mesh *G;
  t=(triangle *)trig;
  p=(point *)pt;
  G=(Mesh*)MM;
  a.x=Mx(t->p0,*G);
  a.y=My(t->p0,*G);
  b.x=Mx(t->p1,*G);
  b.y=My(t->p1,*G);
  c.x=Mx(t->p2,*G);
  c.y=My(t->p2,*G);
  if (is_point_online_p(*p,a,b)
    ||is_point_online_p(*p,b,c)
    ||is_point_online_p(*p,c,a)) {
#ifdef DEBUG
  printf("inclusion1: point on boundary\n");
#endif
  return(2);
  }
 /* else */
  return(INSIDE(a,b,c,*p));
}
#endif/* !EXTRA_SOURCE */

/**************************/
/* for triangles */
static int
equ(const void*a,const void *b)
{ 
  triangle *p,*q;
  p=(triangle *)a;
  q=(triangle *)b;

  return(p->n==q->n);
}
static int
lt(const void*a,const void *b)
{ 
  triangle *p,*q;
  p=(triangle *)a;
  q=(triangle *)b;

  return(p->n<q->n);

}

static void
print_detail(const void *a)
{
  triangle *p;
  p=(triangle *)a;
  printf("%d :",p->n);
  printf("%d:(%lf,%lf), %d:(%lf,%lf), %d:(%lf,%lf) ",p->p0,Mx(p->p0,M),My(p->p0,M),
	 p->p1,Mx(p->p1,M),My(p->p1,M),p->p2,Mx(p->p2,M),My(p->p2,M));
  printf("neighbours %d,%d,%d ",p->t0,p->t1,p->t2);

  p=(triangle *)p->dag_p->rec;
  printf("parent in DAG points to %d\n",p->n);
}

static void
delete_triangle(void *a)
{
  triangle *p;
  p=(triangle *)a;
	if (!p) free(p);
}

/* a function to determine whether point p4 is inside */
/* the circumcircle of points p1,p2,p3 */
static int
inside_circle( int p1, int p2, int p3, int p4)
{
  double temp1,temp2,temp3,temp4,a2,b2;
  double A[3][3];
#ifdef DEBUG
  printf("checking (%d,%d,%d) includes %d\n",p1,p2,p3,p4);
#endif
  temp1 = (Mx(p1,M)*Mx(p1,M)-Mx(p2,M)*Mx(p2,M)+My(p1,M)*My(p1,M)-My(p2,M)*My(p2,M));
  temp2 = (Mx(p1,M)*Mx(p1,M)-Mx(p3,M)*Mx(p3,M)+My(p1,M)*My(p1,M)-My(p3,M)*My(p3,M));
  temp3 = (Mx(p1,M)-Mx(p2,M))*(My(p1,M)-My(p3,M))-(Mx(p1,M)-Mx(p3,M))*(My(p1,M)-My(p2,M));
				/* we know this is not zero */
  temp3 = 1/temp3;
  a2 = (temp1*(My(p1,M)-My(p3,M))-temp2*(My(p1,M)-My(p2,M)))*temp3;
  b2 = (temp2*(Mx(p1,M)-Mx(p2,M))-temp1*(Mx(p1,M)-Mx(p3,M)))*temp3;
  temp4 =  Mx(p1,M)*Mx(p1,M)-Mx(p4,M)*Mx(p4,M)+My(p1,M)*My(p1,M)-My(p4,M)*My(p4,M);
  temp3=a2*(Mx(p1,M)-Mx(p4,M))+b2*(My(p1,M)-My(p4,M));
  if ( ABS(temp3 - temp4) > TOL ) {
    /* can make a simple decision */
   if(temp3 < temp4 ) {
#ifdef DEBUG
    printf("Simple case :points are inside circle\n");
#endif
    return(1);
  } else {
#ifdef DEBUG
    printf("Simple case :points are NOT inside circle\n");
#endif
  return(0);
  }
 } else {
  /* this is a close call. evaluate determinant */
  A[0][0]=Mx(p1,M)-Mx(p4,M);A[0][1]=My(p1,M)-My(p4,M);A[0][2]=A[0][0]*A[0][0]+A[0][1]*A[0][1];
  A[1][0]=Mx(p2,M)-Mx(p4,M);A[1][1]=My(p2,M)-My(p4,M);A[1][2]=A[1][0]*A[1][0]+A[1][1]*A[1][1];
  A[2][0]=Mx(p3,M)-Mx(p4,M);A[2][1]=My(p3,M)-My(p4,M);A[2][2]=A[2][0]*A[2][0]+A[2][1]*A[2][1];
  /* evaluate determinant */
  temp3=A[0][0]*(A[1][2]*A[2][2]-A[2][1]*A[1][2])
      -A[0][1]*(A[1][0]*A[2][2]-A[2][0]*A[1][2])
      +A[0][2]*(A[1][0]*A[2][1]-A[2][0]*A[1][1]);
  if(temp3 <=0) {
#ifdef DEBUG
    printf("Hard case :points are inside circle\n");
#endif
    return(1);
  } else {
#ifdef DEBUG
    printf("Hard case :points are NOT inside circle\n");
#endif
  return(0);
  }
 }
 return(0);
}

/* update neighbour link in a triangle */
/* neighbour old will be replaced by neighbour new */
static void
update_neighbours_around(triangle *neighbour,int old,int new)
{
  if ( neighbour ) {
    if ( neighbour->t0 == old ) {
      neighbour->t0 = new;
    } else if ( neighbour->t1 == old ) {
      neighbour->t1 = new;
    } else if ( neighbour->t2 == old ) {
      neighbour->t2 = new;
    }
  }
}

/* tgn - triangle number to refine */
/* ptn - point number of vertex of tgn facing the edge to be refined */
/* assertion - at each step, the remainder of triangles form a valid 
 *  delaunay triangulation */
int 
delaunay_refine(int tgn, int ptn, RBT_tree * rbt,DAG_tree *dag)
{
  triangle tq,*t,*t_n;
  int neighbour,p0,p1,p2,p4;
  int t1,t2,t3,t4; /* neighbours */
  triangle *tp1,*tp2,*tp3,*tp4,*tnew1,*tnew2;

  DAG_node *current_dag_node,*parent_dag_node;
  /* fix everything into this generalized form 
   *                 p2/\
   *                  /| \
   *                 / |  \
   *              t1/  |   \ t4
   *               /   |    \
   *           p0 /    |     \ p4
   *             \  t  | tn  /
   *              \    |    /
   *           t2  \   |   / t3
   *                 \ | /
   *                   p1
   */

  /* get the triangle */
  tq.n=tgn;
  t=RBT_find(&tq, rbt);
  if ( t && ( t->n == tgn )) { /* triangle exists */
    /* first find neighbour facing the edge to be refined */
    if ( ptn == t->p0 ) {
      neighbour=t->t0;
      p0=t->p0;p1=t->p1;p2=t->p2;
      t1=t->t1; t2=t->t2;
    } else if  ( ptn == t->p1 ) {
      neighbour=t->t1;
      p0=t->p1;p1=t->p2;p2=t->p0;
      t1=t->t2; t2=t->t0;
    } else {
      neighbour=t->t2;
      p0=t->p2;p1=t->p0;p2=t->p1;
      t1=t->t0; t2=t->t1;
    }
		/* if the edge to be flipped is a fixed one,
		 * we cannot do anything, we quit */
		/*if ((Mval(p1)==FX) &&(Mval(p2)==FX)){
					printf("cannot refine\n");
					return(1);
		}*/
    /* now find the neighbouring triangle */
    /* it may not exist */
    tq.n=neighbour;
    t_n=RBT_find(&tq, rbt);
    if ( !t_n ) {/* no neighbour */
      return(0);
    }
    /* find the shared edge with neigbour, or rather, the opposing vertex */
    if ( tgn == t_n->t0 ) {
      p4=t_n->p0;
      t3=t_n->t1; t4=t_n->t2;
    } else if  ( tgn == t_n->t1 ) {
      p4=t_n->p1;
      t3=t_n->t2; t4=t_n->t0;
    } else {
      p4=t_n->p2;
      t3=t_n->t0; t4=t_n->t1;
    }
    /* now we have fitted everything to the general framework */
    if ( inside_circle(p0,p2,p1,p4) ) {
      /* we flip the edge (p1,p2) */
      /* this could only be done if p4 can be seen
       * by p0, i.e. (p1,p2) is the largest edge of neighbour 
       * but since we assert a valid initial triangulation,
       * this will hold.
       */

      /* get neighbouring triangles first */
      tq.n=t1;tp1=RBT_find(&tq, rbt);
      tq.n=t2;tp2=RBT_find(&tq, rbt);
      tq.n=t3;tp3=RBT_find(&tq, rbt);
      tq.n=t4;tp4=RBT_find(&tq, rbt);
      /* create 2 new triangles */
      if ((tnew1=(triangle *)malloc(sizeof(triangle)))==0) {  
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
      }
      if ((tnew2=(triangle *)malloc(sizeof(triangle)))==0) {  
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
      }
      /* we transform the setup into 
       *                 p2/\
       *                  /  \
       *                 /    \
       *              t1/      \ t4
       *               /  new1  \
       *           p0 /          \ p4
       *             \- - ------ /
       *              \  new2   /
       *           t2  \       / t3
       *                 \   /
       *                   p1
       */
      tnew1->p0=p0;tnew1->p1=p4;tnew1->p2=p2;
      tnew1->n=ntriangles;
      tnew2->p0=p0;tnew2->p1=p1;tnew2->p2=p4;
      tnew2->n=ntriangles+1;
      tnew1->status=tnew2->status=LIVE;
      t->status=t_n->status=DEAD;
      /* neighbours */
      tnew1->t0=t4;tnew1->t1=t1;tnew1->t2=tnew2->n;
      tnew2->t0=t3;tnew2->t1=tnew1->n;tnew2->t2=t2;
      /* undate the neighbours of the neighbours too */
      update_neighbours_around(tp1,tgn,tnew1->n);
      update_neighbours_around(tp2,tgn,tnew2->n);
      update_neighbours_around(tp3,neighbour,tnew2->n);
      update_neighbours_around(tp4,neighbour,tnew1->n);
      /* insert the 2 new triangles into RB tree */
      RBT_insert((void*)tnew1, rbt);

      parent_dag_node=t->dag_p;
      current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew1); 
      tnew1->dag_p=current_dag_node;
      /* link it to the other parent in dag */
      DAG_link(dag,t_n->dag_p,current_dag_node);

      /* second triangle */
      RBT_insert((void*)tnew2, rbt);

      parent_dag_node=t->dag_p;
      current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew2); 
      tnew2->dag_p=current_dag_node;
      /* link it to the other parent in dag */
      DAG_link(dag,t_n->dag_p,current_dag_node);

      /* increment triangle count */
      ntriangles +=2;
      /* printf("flipped %d neighbour %d  to %d and %d\n",tgn,neighbour,tnew1->n,tnew2->n); */

      /* call this function recursively */
      delaunay_refine(tnew1->n,p0,rbt,dag);
      delaunay_refine(tnew2->n,p0,rbt,dag);
    }


  } else {
    /* error */
    printf("cannot find triangle\n");
    return(-1);
  }
 

  return(0);
}


/* function to test point is on a line segment */
static double
is_point_online(int p, int p1,int p2)
{
  /* p = p1 + lambda p2/(1+lambda) */
  double x,x1,x2,y,y1,y2,lambda;
 
  x=Mx(p,M);
  x1=Mx(p1,M);
  x2=Mx(p2,M);
  y=My(p,M);
  y1=My(p1,M);
  y2=My(p2,M);

 /*printf("(%lf,%lf)(%lf,%lf)(%lf,%lf)\n",x,y,x1,y1,x2,y2); */
  if ((IS_ZERO(x - x1)&&IS_ZERO(y-y1))||(IS_ZERO(x-x2)&&IS_ZERO(y-y2))) {
    /* points coincide */
    return(-1);
  }
  if (IS_ZERO(x-x1)&&IS_ZERO(x-x2)) {
    return((y-y1)/(y2-y));
  }
  if ( IS_ZERO(y-y1)&&IS_ZERO(y-y2) ) {
    return((x-x1)/(x2-x));
  }
  lambda=(x-x1)/(x2-x);
  if ( IS_ZERO(lambda - (y-y1)/(y2-y)) ) {
    return(lambda);
  } 
  /*printf("l1=%lf,l2=%lf\n",(x-x1)/(x2-x),(y-y1)/(y2-y)); */
  /* else  */
  return(-1);
}

int
insert_point_to_mesh(point p)
{
  DAG_node *current_dag_node,*parent_dag_node;
  triangle *tg,*tg1,*tg2,*tg3,*tempt,tq;
  triangle *tg4,*t0,*t1,*t2,*t3,*neigh;
  int status,p0,p1,p2,p3,point_number,neighbour_number;
  double lambda1,lambda2,lambda3;


  point_number=BIT_insert(&M,p);
#ifdef DEBUG
	printf("insert_point_to_mesh: new point %d\n",point_number);
	printf("insert_point_to_mesh: Insert point count=%d, maximum=%d\n",M.count,M.maxpoints);
#endif
current_dag_node=DAG_find(&dt, (void *)&p,(void*)&M);
    if ( current_dag_node ) {
      /* point is included in this triangle */
      tg=(triangle *)current_dag_node->rec;
#ifdef DEBUG
       printf("insert_p_..:point inside %d->(%d,%d,%d)\n",tg->n,tg->p0,tg->p1,tg->p2); 
#endif

      /* insert point into  bitree */
      lambda1=is_point_online(point_number,tg->p0,tg->p1);
      lambda2=is_point_online(point_number,tg->p1,tg->p2);
      lambda3=is_point_online(point_number,tg->p2,tg->p0);
#ifdef DEBUG
      printf("insert_p_..:lambdas %lf,%lf,%lf\n",lambda1,lambda2,lambda3); 
#endif
      if ( (lambda1<0)&&(lambda2<0)&&(lambda3<0)) {

	/* split the current triangle into 3 */
	/*    p0 /\                  p0 /|\
	 *      /  \                   / | \
	 *     /    \      to         /1 |  \ t1
	 *    /      \           t2  /  / \3 \
	 *   /        \             /  / pn\  \
	 *  /          \           / /       \ \
	 * / p1         \ p2      //    2      \\
	 * --------------         ---------------
	 *                       p1     t0       p2
	 */
	if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
	tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles;
	tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1;
	tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2;
	/* arrange neighbours */
	tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2;
	tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0;
	tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1;
	/* new triangles are alive */
	tg1->status=tg2->status=tg3->status=LIVE;
	/* old triangle is dead */
	tg->status=DEAD;
                
	/* update links from neighbours too */
	tq.n=tg->t1;tempt=RBT_find(&tq, &rt);
	update_neighbours_around(tempt,tg->n,tg3->n);
	tq.n=tg->t2;tempt=RBT_find(&tq, &rt);
	update_neighbours_around(tempt,tg->n,tg1->n);
	tq.n=tg->t0;tempt=RBT_find(&tq, &rt);
	update_neighbours_around(tempt,tg->n,tg2->n);


	parent_dag_node=current_dag_node;
	current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); 
	tg1->dag_p=current_dag_node;
	current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); 
	tg2->dag_p=current_dag_node;
	current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3); 
	tg3->dag_p=current_dag_node;

	status = RBT_insert((void*)tg1, &rt);
	status = RBT_insert((void*)tg2, &rt);
	status = RBT_insert((void*)tg3, &rt);
	ntriangles+=3;

	/* now Delaunay refine the neighbouring triangles */
	delaunay_refine(tg1->n,point_number,&rt,&dt);
	delaunay_refine(tg2->n,point_number,&rt,&dt);
	delaunay_refine(tg3->n,point_number,&rt,&dt);
      } else {
	/* one lambda is positive */
	/* split to 4 */
	/* general form */
				/*                      p0               p0
				 *                      /\               /\
				 *                     /  \  t3         / |\
				 *               t0   /    \           /  | \
				 *                   /  t   \         /  1| 2\
				 *              p1  /____p___\ p3  p1/____|___\ p3
				 *                  \        /       \    | 4 /
				 *                   \  nei  /        \  3|  /
				 *                    \    /           \  | /
				 *               t1    \  /  t2         \ |/
				 *                      \/               \/
				 *                      p2               p2
				 */
	/*printf("split 4\n"); */
	neigh=0;
	if ( lambda1 > 0 ) {
	  p0=tg->p2;
	  p1=tg->p0;
	  p3=tg->p1;
	  tq.n=tg->t2;
neighbour_number=tq.n;
	  neigh=RBT_find(&tq, &rt);

	  tq.n=tg->t1;
	  t0=RBT_find(&tq, &rt);
	  tq.n=tg->t0;
	  t3=RBT_find(&tq, &rt);
             
	} else if ( lambda2 > 0 ) {
	  p0=tg->p0;
	  p1=tg->p1;
	  p3=tg->p2;
	  tq.n=tg->t0;
neighbour_number=tq.n;
	  neigh=RBT_find(&tq, &rt);

	  tq.n=tg->t2;
	  t0=RBT_find(&tq, &rt);
	  tq.n=tg->t1;
	  t3=RBT_find(&tq, &rt);
             
	} else { /* lambda3 > 0 */
	  p0=tg->p1;
	  p1=tg->p2;
	  p3=tg->p0;
	  tq.n=tg->t1;
neighbour_number=tq.n;
	  neigh=RBT_find(&tq, &rt);

	  tq.n=tg->t0;
	  t0=RBT_find(&tq, &rt);
	  tq.n=tg->t2;
	  t3=RBT_find(&tq, &rt);

	}
				/* if neighbour exists */
	if ( neigh ) {
	  /*printf("neighbour exists %d\n",neigh->n); */
	  if ( neigh->t0 == tg->n ) {
	    p2=neigh->p0;
	    tq.n=neigh->t1;
	    t1=RBT_find(&tq,&rt);
	    tq.n=neigh->t2;
	    t2=RBT_find(&tq,&rt);
	  } else if (neigh->t1 == tg->n ) {
	    p2=neigh->p1;
	    tq.n=neigh->t2;
	    t1=RBT_find(&tq,&rt);
	    tq.n=neigh->t0;
	    t2=RBT_find(&tq,&rt);
	  } else { 
	    p2=neigh->p2;
	    tq.n=neigh->t0;
	    t1=RBT_find(&tq,&rt);
	    tq.n=neigh->t1;
	    t2=RBT_find(&tq,&rt);
	  }

	/*	printf("generalized %d %d %d %d\n",p0,p1,p2,p3); */
	/* split to 4 */
if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
	  tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
	  tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
	  tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2;
	  tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3;

	  /* alive, dead triangles */
	  /* new triangles are alive */
	  tg1->status=tg2->status=tg3->status=tg4->status=LIVE;
	  /* old triangle is dead */
	  tg->status=neigh->status=DEAD;

	  /* arrange neighbours */
	  tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0 ? t0->n:-1);
	  tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
	  tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n;
	  tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n;
	  /* update old neighbours */
	  update_neighbours_around(t0,tg->n,tg1->n);
	  update_neighbours_around(t3,tg->n,tg2->n);
	  update_neighbours_around(t1,neigh->n,tg3->n);
	  update_neighbours_around(t2,neigh->n,tg4->n);

	  /* insert into RBT */
	  RBT_insert((void*)tg1, &rt);
	  RBT_insert((void*)tg2, &rt);
	  RBT_insert((void*)tg3, &rt);
	  RBT_insert((void*)tg4, &rt);
	  ntriangles+=4;

	  /* update DAG */
	  parent_dag_node=tg->dag_p;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); 
	  tg1->dag_p=current_dag_node;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); 
	  tg2->dag_p=current_dag_node;
	  parent_dag_node=neigh->dag_p;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3); 
	  tg3->dag_p=current_dag_node;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg4); 
	  tg4->dag_p=current_dag_node;

  delaunay_refine(tg1->n,point_number,&rt,&dt);
	delaunay_refine(tg2->n,point_number,&rt,&dt);
	delaunay_refine(tg3->n,point_number,&rt,&dt);
	delaunay_refine(tg4->n,point_number,&rt,&dt);

	 /* printf("split to %d %d %d %d\n",tg1->n,tg2->n,tg3->n,tg4->n);
      DAG_traverse(&dt); */
	} else { /* no neighbour */
if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}

	  tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
	  tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;

	  /* alive, dead triangles */
	  /* new triangles are alive */
	  tg1->status=tg2->status=LIVE;
	  /* old triangle is dead */
	  tg->status=DEAD;

	  /* arrange neighbours */
	  tg1->t0=neighbour_number;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
	  tg2->t0=neighbour_number;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
	  /* update old neighbours */
	  update_neighbours_around(t0,tg->n,tg1->n);
	  update_neighbours_around(t3,tg->n,tg2->n);

	  /* insert into RBT */
	  RBT_insert((void*)tg1, &rt);
	  RBT_insert((void*)tg2, &rt);
	  ntriangles+=2;

	  /* update DAG */
	  parent_dag_node=tg->dag_p;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); 
	  tg1->dag_p=current_dag_node;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); 
	  tg2->dag_p=current_dag_node;
  delaunay_refine(tg1->n,point_number,&rt,&dt);
	delaunay_refine(tg2->n,point_number,&rt,&dt);

	}
      }

    }
	return(point_number);
}

/***************************/
/* remove (not literally) triangles outside a polygon */
static void
remove_exterior_and_update_triangles(boundary *b)
{
	triangle *rec;
	point p;
 int i;
#ifdef DEBUG
 int j;
 triangle brec;
#endif
	/* remove outer boundary triangles */
	 DAG_traverse_list_reset(&dt);
  rec=DAG_traverse_prune_list(&dt);
	while (rec) {
     if (rec &&rec->status==LIVE) {
       p.x=(Mx(rec->p0,M)+Mx(rec->p1,M)+Mx(rec->p2,M))/3; 
       p.y=(My(rec->p0,M)+My(rec->p1,M)+My(rec->p2,M))/3; 
		   if (!inside_poly_point(&(b->poly_array[0]),p)) {
			     rec->status=HIDDEN;
#ifdef DEBUG
		printf("remove_ext: updated triangle %d(%d,%d,%d) status %d\n",rec->n,rec->p0,rec->p1,rec->p2,rec->status);
#endif


			 }
		 }

   rec=DAG_traverse_prune_list(&dt);
 }

		/* update other boundary triangles */
  for (i=0;i<b->npoly;i++) {
   DAG_traverse_list_reset(&dt);
  rec=DAG_traverse_prune_list(&dt);
	while (rec) {
     if (rec &&rec->status==LIVE) {
       p.x=(Mx(rec->p0,M)+Mx(rec->p1,M)+Mx(rec->p2,M))/3; 
       p.y=(My(rec->p0,M)+My(rec->p1,M)+My(rec->p2,M))/3; 
		   if (inside_poly_point(&(b->poly_array[i]),p)) {
			     rec->boundary=i;
			 /* if hollow boundary, make triangle hidden*/
			    if ( (b->poly_array[i]).hollow==0 ) {
								rec->status=HIDDEN;
					}
			 }
		 }
   rec=DAG_traverse_prune_list(&dt);
 }
		}

#ifdef DEBUG
	printf("===============\n");
		for (j=0;j<ntriangles;j++) {
			brec.n=j;
    rec= RBT_find(&brec, &rt);

     if (rec->status!=DEAD) {
		printf("remove_ext: check triangle %d(%d,%d,%d) status %d\n",rec->n,rec->p0,rec->p1,rec->p2,rec->status);
					}
		} 
#endif
}

/*************************/
/* generate mesh */
void
generate_mesh(const char *filename) {

  DAG_node *current_dag_node,*parent_dag_node;
  triangle *tg,*tg1,*tg2,*tg3,*tempt,tq;
  triangle *tg4,*t0,*t1,*t2,*t3,*neigh;
  int status,p0,p1,p2,p3,total,i,point_number;
  point p;
		point *point_array;
  double lambda1,lambda2,lambda3;
  polygon temp_polygon;
  poly_edge e;

  /* read input */
  total=read_input_file(&point_array,filename);
  
  if ((tg=(triangle *)malloc(sizeof(triangle)))==0){
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
  }
  p.x=point_array[0].x;p.y=point_array[0].y; 
  if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 p.z[0]=0;
		p.val=FX; /* we dont care they are FX or not */
#ifdef DEBUG 
    printf("inserting point %d "MDF","MDF"\n",0,p.x,p.y);
#endif
  tg->p0=BIT_insert(&M, p );
  p.x=point_array[1].x;p.y=point_array[1].y; 
		if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
  p.z[0]=0;
	p.val=FX;
#ifdef DEBUG 
    printf("inserting point %d "MDF","MDF"\n",1,p.x,p.y);
#endif
  tg->p1=BIT_insert(&M, p );
  p.x=point_array[2].x;p.y=point_array[2].y; 
		if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
  p.z[0]=0;
		p.val=FX;
#ifdef DEBUG 
    printf("inserting point %d "MDF","MDF"\n",2,p.x,p.y);
#endif
  tg->p2=BIT_insert(&M, p );

  tg->n=0;
  /* triangle is alive */
  tg->status=LIVE;
  /* make neighbours -1 */
  tg->t0=tg->t1=tg->t2=-1;
  p0=tg->p0;p1=tg->p1;p2=tg->p2;
  ntriangles=1;

  status = RBT_insert((void*)tg, &rt);
  parent_dag_node=dt.root;
  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg); 
  /* set pointer to dag */
  tg->dag_p=current_dag_node;
  /* we have inserted the outer triangle */
  /* now process the remaining points */
  for (i=3;i<total;i++) {
    p.x=point_array[i].x;p.y=point_array[i].y;p.val=point_array[i].val;
    p.z=point_array[i].z;
    point_number=BIT_insert(&M,p);
#ifdef DEBUG 
    printf("do_intial_triangulation: inserting point %d "MDF","MDF" z="MDF" as %d\n",i,p.x,p.y,*p.z,point_number);
#endif
				if ( point_number != i) {
									fprintf(stderr,"***: insertion of new point %d returned an old point %d.\n your input file is screwed up. you may have duplicate points. aborting..\n",i,point_number);
									abort();
				}
    current_dag_node=DAG_find(&dt, (void *)&p,(void*)&M);
    if ( current_dag_node ) {
      /* point is included in this triangle */
      tg=(triangle *)current_dag_node->rec;
#ifdef DEBUG 
      printf("do_intial_triangulation: point inside %d->(%d,%d,%d)\n",tg->n,tg->p0,tg->p1,tg->p2);
#endif
      /* insert point into  bitree */
      lambda1=is_point_online(point_number,tg->p0,tg->p1);
      lambda2=is_point_online(point_number,tg->p1,tg->p2);
      lambda3=is_point_online(point_number,tg->p2,tg->p0);
#ifdef DEBUG
      printf("lambdas %lf,%lf,%lf\n",lambda1,lambda2,lambda3);
#endif
      if ( (lambda1<0)&&(lambda2<0)&&(lambda3<0)) {

	/* split the current triangle into 3 */
	/*    p0 /\                  p0 /|\
	 *      /  \                   / | \
	 *     /    \      to         /1 |  \ t1
	 *    /      \           t2  /  / \3 \
	 *   /        \             /  / pn\  \
	 *  /          \           / /       \ \
	 * / p1         \ p2      //    2      \\
	 * --------------         ---------------
	 *                       p1     t0       p2
	 */
	if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
	tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles;
	tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1;
	tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2;
	/* arrange neighbours */
	tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2;
	tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0;
	tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1;
	/* new triangles are alive */
	tg1->status=tg2->status=tg3->status=LIVE;
	/* old triangle is dead */
	tg->status=DEAD;
                
	/* update links from neighbours too */
	tq.n=tg->t1;tempt=RBT_find(&tq, &rt);
	update_neighbours_around(tempt,tg->n,tg3->n);
	tq.n=tg->t2;tempt=RBT_find(&tq, &rt);
	update_neighbours_around(tempt,tg->n,tg1->n);
	tq.n=tg->t0;tempt=RBT_find(&tq, &rt);
	update_neighbours_around(tempt,tg->n,tg2->n);


	parent_dag_node=current_dag_node;
	current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); 
	tg1->dag_p=current_dag_node;
	current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); 
	tg2->dag_p=current_dag_node;
	current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3); 
	tg3->dag_p=current_dag_node;

	status = RBT_insert((void*)tg1, &rt);
	status = RBT_insert((void*)tg2, &rt);
	status = RBT_insert((void*)tg3, &rt);
	ntriangles+=3;

	/* now Delaunay refine the neighbouring triangles */
	delaunay_refine(tg1->n,point_number,&rt,&dt);
	delaunay_refine(tg2->n,point_number,&rt,&dt);
	delaunay_refine(tg3->n,point_number,&rt,&dt);
      } else {
	/* one lambda is positive */
	/* split to 4 */
	/* general form */
				/*                      p0               p0
				 *                      /\               /\
				 *                     /  \  t3         / |\
				 *               t0   /    \           /  | \
				 *                   /  t   \         /  1| 2\
				 *              p1  /____p___\ p3  p1/____|___\ p3
				 *                  \        /       \    | 4 /
				 *                   \  nei  /        \  3|  /
				 *                    \    /           \  | /
				 *               t1    \  /  t2         \ |/
				 *                      \/               \/
				 *                      p2               p2
				 */
	/*printf("split 4\n"); */
	neigh=0;
	if ( lambda1 > 0 ) {
	  p0=tg->p2;
	  p1=tg->p0;
	  p3=tg->p1;
	  tq.n=tg->t2;
	  neigh=RBT_find(&tq, &rt);

	  tq.n=tg->t1;
	  t0=RBT_find(&tq, &rt);
	  tq.n=tg->t0;
	  t3=RBT_find(&tq, &rt);
             
	} else if ( lambda2 > 0 ) {
	  p0=tg->p0;
	  p1=tg->p1;
	  p3=tg->p2;
	  tq.n=tg->t0;
	  neigh=RBT_find(&tq, &rt);

	  tq.n=tg->t2;
	  t0=RBT_find(&tq, &rt);
	  tq.n=tg->t1;
	  t3=RBT_find(&tq, &rt);
             
	} else { /* lambda3 > 0 */
	  p0=tg->p1;
	  p1=tg->p2;
	  p3=tg->p0;
	  tq.n=tg->t1;
	  neigh=RBT_find(&tq, &rt);

	  tq.n=tg->t0;
	  t0=RBT_find(&tq, &rt);
	  tq.n=tg->t2;
	  t3=RBT_find(&tq, &rt);

	}
				/* if neighbour exists */
	if ( neigh ) {
	  /*printf("neighbour exists %d\n",neigh->n); */
	  if ( neigh->t0 == tg->n ) {
	    p2=neigh->p0;
	    tq.n=neigh->t1;
	    t1=RBT_find(&tq,&rt);
	    tq.n=neigh->t2;
	    t2=RBT_find(&tq,&rt);
	  } else if (neigh->t1 == tg->n ) {
	    p2=neigh->p1;
	    tq.n=neigh->t2;
	    t1=RBT_find(&tq,&rt);
	    tq.n=neigh->t0;
	    t2=RBT_find(&tq,&rt);
	  } else { 
	    p2=neigh->p2;
	    tq.n=neigh->t0;
	    t1=RBT_find(&tq,&rt);
	    tq.n=neigh->t1;
	    t2=RBT_find(&tq,&rt);
	  }

	/*	printf("generalized %d %d %d %d\n",p0,p1,p2,p3); */
	/* split to 4 */
if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
	  tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
	  tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
	  tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2;
	  tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3;

	  /* alive, dead triangles */
	  /* new triangles are alive */
	  tg1->status=tg2->status=tg3->status=tg4->status=LIVE;
	  /* old triangle is dead */
	  tg->status=neigh->status=DEAD;

	  /* arrange neighbours */
	  tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0 ? t0->n:-1);
	  tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
	  tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n;
	  tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n;
	  /* update old neighbours */
	  update_neighbours_around(t0,tg->n,tg1->n);
	  update_neighbours_around(t3,tg->n,tg2->n);
	  update_neighbours_around(t1,neigh->n,tg3->n);
	  update_neighbours_around(t2,neigh->n,tg4->n);

	  /* insert into RBT */
	  RBT_insert((void*)tg1, &rt);
	  RBT_insert((void*)tg2, &rt);
	  RBT_insert((void*)tg3, &rt);
	  RBT_insert((void*)tg4, &rt);
	  ntriangles+=4;

	  /* update DAG */
	  parent_dag_node=tg->dag_p;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); 
	  tg1->dag_p=current_dag_node;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); 
	  tg2->dag_p=current_dag_node;
	  parent_dag_node=neigh->dag_p;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3); 
	  tg3->dag_p=current_dag_node;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg4); 
	  tg4->dag_p=current_dag_node;

  delaunay_refine(tg1->n,point_number,&rt,&dt);
	delaunay_refine(tg2->n,point_number,&rt,&dt);
	delaunay_refine(tg3->n,point_number,&rt,&dt);
	delaunay_refine(tg4->n,point_number,&rt,&dt);

	  /*printf("split to %d %d %d %d\n",tg1->n,tg2->n,tg3->n,tg4->n);
      DAG_traverse(&dt); */
	} else { /* no neighbour */
if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}
if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
	fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	exit(1);
	}

	  tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
	  tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;

	  /* alive, dead triangles */
	  /* new triangles are alive */
	  tg1->status=tg2->status=LIVE;
	  /* old triangle is dead */
	  tg->status=DEAD;

	  /* arrange neighbours */
	  tg1->t0=tg->t2;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
	  tg2->t0=tg->t2;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
	  /* update old neighbours */
	  update_neighbours_around(t0,tg->n,tg1->n);
	  update_neighbours_around(t3,tg->n,tg2->n);

	  /* insert into RBT */
	  RBT_insert((void*)tg1, &rt);
	  RBT_insert((void*)tg2, &rt);
	  ntriangles+=2;

	  /* update DAG */
	  parent_dag_node=tg->dag_p;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); 
	  tg1->dag_p=current_dag_node;
	  current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); 
	  tg2->dag_p=current_dag_node;
  delaunay_refine(tg1->n,point_number,&rt,&dt);
	delaunay_refine(tg2->n,point_number,&rt,&dt);

	}
      }

    }
  }

	  /* free point array */
		free(point_array);
		/* note: we dont free point_array[i].z values 
		 * because they are stored in the tree */

#ifdef DEBUG 
  DAG_traverse(&dt);
#endif
	/* copy boundary to a temp one */
  /*copy_polygon(&bound,&temp_polygon); */
  /* create a temporary polygon with all edges */
	init_poly(&temp_polygon);
  for (i=0; i<nedges; i++){
	 e.p1 =edge_array[i].p1; e.p2=edge_array[i].p2; 
		e.type=edge_array[i].type;
#ifdef DEBUG
		printf("do_initial_triangulation: building plygon edge (%d,%d)\n",e.p1,e.p2);
#endif
	 insert_edge_to_poly(&temp_polygon,&e);
	}

  remove_polygon_mesh_intersections(&temp_polygon);

  destroy_polygon(&temp_polygon);

	/* remove exterior triangles */
  remove_exterior_and_update_triangles(&bound);

		/* subdivide triangles */
  subdivide_all_triangles();
#ifdef DEBUG
		for (i=0; i<M.count;i++) {
		  printf("triangulate: %d (%lf,%lf) = %lf\n",i,Mx(i,M),My(i,M),Mz(i,M));
		}
#endif

}

/* function to generate mesh and solve */
int 
solve_problem_from_scratch(const char *filename)
{

	/* set the degree of freedom  of points */
	degree_of_freedom=1;
/* (re) generate mesh */
  generate_mesh(filename);

		/* solve mesh */
  if ( solve_equation!=POISSON && solve_equation!=POISSON_SPARSE ) {
  /* number of eigenvectors to be found */
	if (requested_degree_of_freedom) {
	 degree_of_freedom=requested_degree_of_freedom; /* global change */
  } else {
   degree_of_freedom=3;
  }
 if((eig_array= (MY_DOUBLE*) realloc((void*)eig_array,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
 if ( solve_equation==HELMHOLTZ) {
   re_number_and_solve_helmholtz();
	} else if ((solve_equation == HELMHOLTZ_INHOMO) 
     || (solve_equation == HELMHOLTZ_FREQ)
     || (solve_equation == HELMHOLTZ_BETA)) {
   buildup_hash_table_and_equation(M,dt);
  }
  } else if (solve_equation==POISSON_SPARSE) {
   re_number_and_solve_poisson_sparse();
  }  else { /* POISSON */
   re_number_and_solve_poisson();
  }

  return(0);
}


/* free all allocated variables */
void
free_everything(void)
{
	int i;
	/* we dont free anything if M.count==0 */
	/* because we assume we have not allocated any memory */
	if (M.count) {
	/* Mesh */
	BIT_free(&M);
	/* RBT */
	RBT_free(&rt);
	/* DAG */
  DAG_free(&dt);

	/* free edge array */
	free(edge_array);
	edge_array=0;
  /* free boundaries */
	for (i=0;i<bound.npoly;i++) {
			destroy_polygon(&(bound.poly_array[i]));
	}
	free(bound.poly_array);
	bound.poly_array=0;
	bound.npoly=0;
 }
  /* initialize bitree for 100 points intially */
  BIT_init(&M,100);
  /* initialize DAG */
  DAG_init(&dt,&inclusion,&print_detail);
  /* initialize RBT with triangles */
  RBT_init(&rt,&equ,&lt,&print_detail,&delete_triangle);

}

/* function to print help */
static void 
print_help(void)
{
  printf("usage: "PACKAGE" [options] \n\n");
  printf("options are:\n");
  printf("-i (FILE) : input file name is 'FILE'\n");
  printf("-a (val)  : set min area of each triangle to 'val'\n");
  printf("-A (val)  : set max area of each triangle to 'val'\n");
  printf("-h, -?    : help, print info\n");
/*  printf("-v    : be verbose\n");                                         
  printf("-b (val) : set max badness of each triangle to 'val'\n");            
  printf("-s (val) : set splitting priority to 'val'\n");                 
  printf("-l (val) : set max nodes to 'val'\n");
  printf("-r (val) : limit iterations to 'val'\n");
  printf("-conj    : use preconditioned conjugate gradient method to solve\n");
  printf("-do      : output data file of mesh as \"data.in\"\n");     */
  printf("-c (val)  : set contour levels to 'val'\n");                     
	printf("-d (contour|grad|fill): display contours,gradient,colourfill by default\n");
	printf("-e (poisson|helmholtz|helmin|helmk|helmbeta): type of equation to solve\n");
  printf("-s        : use sparse solver (works only for Poisson Equation)\n");
  printf("-x        : do not go into display mode, dump output to file 'FILE.eps'\n");
  printf("-V        : print version\n");
}

/* function to print version */
static void
print_version(void)
{
  printf(PACKAGE"  "VERSION", Copyright (C) 2001-2005 Sarod Yatawatta\n");
  printf(PACKAGE" comes with ABSOLUTELY NO WARRANTY\n"); 
  printf("This is free software, and you are welcome to redistribute\n");
  printf("it under certain conditions given by GNU GENERAL PUBLIC LICENSE\n");
  printf("$Id: pdnmesh.c,v 1.57 2005/04/24 05:32:24 sarod Exp $\n"); 
	printf("This is a beta release. Please report bugs to\n");
	printf(" <"PACKAGE_BUGREPORT">. Thank you.\n");
}


#ifndef WIN32
void 
initialize_global_variables(int argc, char *argv[])
{

	int c;
	opterr=0;

contour_levels=40; /* no of contour levels */
current_plotting_contour=0; /* number in z array < degree_of_freedom */
g_badness_limit=2; /* max triangle badness */
g_area_floor=0.001; /* min triangle area */
g_area_ceil=0.005; /* max triangle area */
max_iteration_limit=20;

/* default waveguide parameters */
global_wave_k0=1.0;
global_wave_beta=10.0;

g_gradient_limit=1000;

solve_equation = POISSON; /* POISSON, HELMHOLTZ */

/* default plotting */
 plot_mesh=0;
 plot_cont=0;
 plot_grad=0;
 plot_legend=0;
 text_mode_only=0;

 mesh_already_present=0;

 update_status_bar("Starting intialization...");
/* getopt loop */
 while ((c=getopt(argc,argv,"i:e:d:a:A:c:xsvhV"))!=-1) {
		switch (c) {
						case 'i':
        if (optarg) {
  input_cord_filename=(char*)realloc((void*)input_cord_filename,sizeof(char)*(size_t)(strlen((char*)optarg)+1));
	if ( input_cord_filename == 0 ) {
	  fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
	  exit(1);
	}
  strcpy(input_cord_filename,(char*)optarg);
 } else {
 input_cord_filename=0;
 }
			break;
			      case 'e':
               if (!strcmp(optarg,"helmholtz")) {
							   solve_equation=HELMHOLTZ;
	              } else if( !strcmp(optarg,"poisson")) {
							   solve_equation=POISSON;
	              } else if( !strcmp(optarg,"helmin")) {
							   solve_equation=HELMHOLTZ_INHOMO;
	              } else if( !strcmp(optarg,"helmk")) {
							   solve_equation=HELMHOLTZ_FREQ;
	              } else if( !strcmp(optarg,"helmbeta")) {
							   solve_equation=HELMHOLTZ_BETA;
	              } else {
								fprintf(stderr,"Unknown equation type `%s`.\n",optarg);
							  }
							break;
			      case 'd':
               if (!strcmp(optarg,"contour")) {
							   plot_cont=1;
	              } else if( !strcmp(optarg,"grad")) {
							   plot_grad=1;
	              } else if( !strcmp(optarg,"fill")) {
							   plot_fill=1;
	              } else {
								fprintf(stderr,"Unknown display option `%s`.\n",optarg);
							  }
							break;
            case 'a':
               g_area_floor=strtod(optarg,0);
               if ( g_area_floor == 0.0f ) {
								fprintf(stderr,"Invalid value for area `%s`.\n",optarg);
                g_area_floor=0.001;
               }
             break;
            case 'A':
               g_area_ceil=strtod(optarg,0);
               if ( g_area_ceil== 0.0f ) {
								fprintf(stderr,"Invalid value for area `%s`.\n",optarg);
                g_area_ceil=0.005;
               }
             break;
            case 'c':
               contour_levels=strtol(optarg,0,0);
               if ( !contour_levels ) {
								fprintf(stderr,"Invalid value for contours `%s`.\n",optarg);
                contour_levels=40;
               }
             break;
           	case 's':
					  if(solve_equation==POISSON) solve_equation=POISSON_SPARSE;
            break;
           	case 'x':
					   text_mode_only=1;
            break;
						case 'h':
			      print_help();
						exit(0);
            case 'V':
			      print_version();
						exit(0);

						case '?':
             if (isprint(optopt))
								fprintf(stderr,"Unknown option `-%c`.\n",optopt);
						 else
							  fprintf(stderr,
										"Unknown option character `\\x%x`.\n",optopt);
						exit(1);
					 default:
						break;
		}
	}

 if (!plot_cont && !plot_grad && !plot_fill) {
   plot_cont=1;
 }
  /* initialize bitree for 100 points intially */
  BIT_init(&M,100);
  /* initialize DAG */
  DAG_init(&dt,&inclusion,&print_detail);
  /* initialize RBT with triangles */
  RBT_init(&rt,&equ,&lt,&print_detail,&delete_triangle);



 /* set eig_array to null */
 if((eig_array= (MY_DOUBLE*) calloc((size_t)degree_of_freedom,sizeof(MY_DOUBLE)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

 update_status_bar("Done intialization.");
}
#endif /*WIN32 */
/******************************************/

#ifndef WIN32
int 
main(int argc, char **argv)
{
/* gtkglext */
 GtkWidget *window;
 GdkGLConfig *glconfig;
 char *output_filename;

/* use getopt to parse commandline */
initialize_global_variables(argc,argv);

if ( input_cord_filename ) {
    solve_problem_from_scratch(input_cord_filename);
}
/*free_everything(); */

if(text_mode_only && input_cord_filename) {
  /* allocate memory for output filename */
 if((output_filename= (char*) malloc((strlen(input_cord_filename)+5)*sizeof(char)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
  strcpy(output_filename,input_cord_filename);
  strcat(output_filename,".eps");
  /* save output file */
  if (plot_cont) {
   print_contour_all(output_filename,1,M);  
  }
  if (plot_grad) {
   print_gradient(output_filename,M);  
  }

  free(output_filename);
  return(0);
 }
 /* intialiaze gtk */
 gtk_init(&argc,&argv);
 /* intialize GtkGlExt */
 gtk_gl_init(&argc,&argv);
 /* configure OpenGL framebuffer */
 glconfig=configure_gl();

 /* create and show the application window */
 window=create_window(glconfig);
 gtk_widget_show(window);
 gtk_main();
  return(0);
}
#endif/*WIN32 */