/*  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: input.c,v 1.15 2005/02/15 06:14:20 sarod Exp $
*/


/* IO rutines */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "types.h"

poly_edge *edge_array;
int nedges;
boundary bound; /* boundaries */
MY_DOUBLE g_xoff,g_yoff; /* x,y offsets */
MY_DOUBLE g_xscale,g_yscale; /* x,y scales */
MY_DOUBLE g_xrange,g_yrange; /* x,y ranges, used for contour plot */

/* global no for current boundary */
MY_INT  current_boundary_number;

/* skips comment lines */
static int
skip_lines(FILE *fin)
{

		int c;
		do {
		if ( ( c = getc(fin) ) == EOF )
		return(-1);
	/* handle empty lines */
	if ( c == '\n' )
	continue; /* next line */
	if ( (c != '#') ) {
	ungetc(c,fin);
	return(0);
	} else { /* skip this line */
	do {
	if ( ( c = getc(fin) ) == EOF )
	return(-1);
	} while (  c != '\n') ;
	}
	} while( 1 );
}

/* skips rest of line */
static int
skip_restof_line(FILE *fin)
{
	int c;
	do {
	if ( ( c = getc(fin) ) == EOF )
	return(-1);
	} while (  c != '\n') ;
	return(1);
}

int
read_input_file(point **input_points, const char *infile)
{
 /* read points from infile and store coords in x,y, arrays */
 FILE *infd;
 int i,j,count;
	MY_DOUBLE xmax,xmin,ymax,ymin;
	int npoint, input_degree_of_freedom=1;
 poly_edge e; 
int buff_len,k,c;
char *buff;

 infd =fopen( infile, "r" );
 if ( infd == NULL ) {
   fprintf(stderr, "%s: %d: could not open file %s\n", __FILE__,__LINE__,infile);
   exit(1);
 }
	skip_lines(infd);
	/* numer of points */
 if ( fscanf(infd, "%d" , &npoint ) == EOF ) {
    fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    exit(1);
 }
	skip_restof_line(infd);
	skip_lines(infd);
#ifdef DEBUG
  printf("input: reading %d points\n",npoint);
#endif

 /* allocate memory for points */
 if ((*input_points=(point*)malloc((size_t)(npoint+3)*sizeof(point)))==0) {
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	xmax=-1000;xmin=1000;
	ymax=-1000;ymin=1000;

 for (i=3; i<npoint+3; i++ ) {
		/* allocate memory for z array */
	if (input_degree_of_freedom) {
  if (((*input_points)[i].z=(MY_DOUBLE*)malloc((size_t)(input_degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }
	}
	/* read point data according to degree of freedom */
	 if ( input_degree_of_freedom==0 ) { /* wave equation, do not need this */
    if ( fscanf(infd,"%d "MDF" "MDF"  %d",&count, &((*input_points)[i].x),&((*input_points)[i].y),&j) == EOF ) {
      fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    }
   } else if ( input_degree_of_freedom==1 ) { /* potentials */
    if ( fscanf(infd,"%d "MDF" "MDF" "MDF" %d",&count, &((*input_points)[i].x),&((*input_points)[i].y),&((*input_points)[i].z[0]),&j) == EOF ) {
      fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    }
  } else { /* input_degree_of_freedom==2, elasticity problems etc. */
    if ( fscanf(infd,"%d "MDF" "MDF" "MDF" "MDF" %d",&count, &((*input_points)[i].x),&((*input_points)[i].y),&((*input_points)[i].z[0]),&((*input_points)[i].z[1]), &j) == EOF ) {
      fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    }
  }

#ifdef DEBUG
  printf("input: read point %d as %d ("MDF", "MDF"): input degree_of_freedom=%d ",i,count,(*input_points)[i].x, (*input_points)[i].y,input_degree_of_freedom);
	if ( input_degree_of_freedom==0 ) {
  printf("no z values\n");
  } else if (input_degree_of_freedom==1) {
   printf("z[0]="MDF", j=%d\n",(*input_points)[i].z[0],j);
	} else {
   printf("z[0]="MDF", z[1]="MDF",j=%d\n",(*input_points)[i].z[0], (*input_points)[i].z[1], j);
	}
#endif
 skip_restof_line(infd);
 skip_lines(infd);

 if ( xmin > (*input_points)[i].x ) {
								xmin=(*input_points)[i].x;
	}
	if ( xmax < (*input_points)[i].x ) {
								xmax=(*input_points)[i].x;
	}
 if ( ymin > (*input_points)[i].y ) {
								ymin=(*input_points)[i].y;
	}
	if ( ymax < (*input_points)[i].y ) {
								ymax=(*input_points)[i].y;
	}
/* point type */
	if (j==1) {
					/* fixed point */
					(*input_points)[i].val=FX;
	} else {
					/* variable point */
					(*input_points)[i].val=VAR;
	}

 }
 g_xoff=-(xmax+xmin)*0.5;
 g_yoff=-(ymax+ymin)*0.5;
 g_xrange=(xmax-xmin)*0.5;
 g_yrange=(ymax-ymin)*0.5;

	 /* now find the maximum value of coordinates */
  ymax = ABS(ymax);
		xmax = ABS(xmax);
		xmin = ABS(xmin);
		ymin = ABS(ymin);
		if ( xmax < ymax )
			xmax = ymax;
		if ( xmax < ymin )
			xmax = ymin;
		if ( xmax < xmin  )
			xmax = xmin;
  g_xscale=xmax;
  g_yscale=xmax;
		/* first 3 points (-3M,-3M),(3M,0),(0,3M) */
		 for (i=3;i<npoint+3;i++) {
          (*input_points)[i].x+=(g_xoff);
          (*input_points)[i].y+=(g_yoff);

          (*input_points)[i].x/=(xmax);
          (*input_points)[i].y/=(xmax);
		 }
			(*input_points)[0].x=-3;(*input_points)[0].y=-3;
			(*input_points)[1].x=3;(*input_points)[1].y=0;
			(*input_points)[2].x=0;(*input_points)[2].y=3;


	/* now the edges */
	/* first number of edges */
   if ( fscanf(infd, "%d" , &nedges) == EOF ) {
    fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    exit(1);
 }
 skip_restof_line(infd);
	skip_lines(infd);


#ifdef DEBUG
  printf("input: reading %d edges\n",nedges);
#endif

 /* allocate memory for edge array */
 if ((edge_array =(poly_edge *)malloc((size_t)nedges*sizeof(poly_edge)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

 for (i=0; i< nedges; i++) {
 /* FIXME - edges need not be in number order */
  if ( fscanf(infd, "%d %d %d %d" , &count, &(edge_array[i].p1), &(edge_array[i].p2), &j) == EOF ) {
	      fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
				exit(1);
		}
 skip_restof_line(infd);
	skip_lines(infd);


		/* account for the 3 outer points, so increment point numbers */
   edge_array[i].p1+=3;edge_array[i].p2+=3;
			/* dirichlet or neumann ? */
			if (j==1) {
				edge_array[i].type=DR;
			} else {
				edge_array[i].type=NU;
			}
#ifdef DEBUG
  printf("input: edge %d (%d,%d) ?%d\n",i,edge_array[i].p1,edge_array[i].p2,j);
#endif

 }

	/* read the boundaries */
 /* first number of boundaries */
   if ( fscanf(infd, "%d" , &bound.npoly) == EOF ) {
    fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    exit(1);
 }
 skip_restof_line(infd);
	skip_lines(infd);


	/* allocate memory for boundaries */
 if ((bound.poly_array=(polygon*)malloc((size_t)bound.npoly*sizeof(polygon)))==0){
   fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
   exit(1);
 }

for (count=0; count< bound.npoly; count++) {
 skip_lines(infd);
	/* now read the outer boundary */
	/* first, number of edges, hollow? 1/mu rho */
	if ( fscanf(infd, "%d %d "MDF" " , &i, 
							&bound.poly_array[count].hollow,
							&bound.poly_array[count].mu) == EOF ) {
    fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
    exit(1);
 }
/* since we are reading 1/mu get the reciprocal */
 if (bound.poly_array[count].mu >TOL ){
  bound.poly_array[count].mu=1/bound.poly_array[count].mu;
 } else {
    fprintf(stderr, "%s: %d: Warning: invalid value for 1/mu or epsilon in file %s\n", __FILE__,__LINE__,infile);
 }
 /* set epsilon to 1.0 - free space */
 bound.poly_array[count].epsilon=1.0;
	/* now read the expression or constant for rho */
	/* worst case we assume it is an expression -> read it to a buffer */
	/* then check the buffer for any alphabetical chars */

	  /* now read expression only if boundary not hollow ? FIXME*/
	    buff_len = 30; k = 0;
			buff = (char*)malloc(sizeof(char)*(size_t)(buff_len+2));
			while ( ((c = fgetc(infd)) != EOF ) && k < buff_len) {
						buff[k++] = c;
						if ( c == '\n')
						break;
			}
			if ( c != '\n') {
					buff_len += 30;
					buff = (char*)realloc((void*)buff,sizeof(char)*(size_t)(buff_len+2));
					buff[k++] = c;
					while ( ((c = fgetc(infd)) != EOF ) && k < buff_len ) {
							buff[k++] = c;
							if ( c == '\n')
							break;
					}
			}
			buff[k++] = ';';
			buff[k++] = '\0';

			
			/* now see if all chars in the buffer are +,-,. or digits */
			c=0; /* flag */
			for (k=0;k<(int)strlen(buff)-1;k++) {
							if ( isalpha(buff[k]) || buff[k]=='(' 
								|| buff[k]==')' || buff[k]=='*'
								|| buff[k]=='/' || buff[k]=='^') {
									c=1;/* is an expression */
									break;
								}
			 }
			if (c ) {
      /* make global boundary this -- used in parser */
			current_boundary_number = count;
			/* now parse equation */
#ifdef DEBUG
			printf("input: parsing rho expression for "MIF" using %s\n",current_boundary_number,buff);
#endif
        bound.poly_array[count].rho=0;
				equation_parse(buff);
			} else {
#ifdef DEBUG
			printf("input: parsing rho value for "MIF" using %s\n",current_boundary_number,buff);
#endif
        bound.poly_array[count].rho=(MY_DOUBLE)strtod(buff,0);
				bound.poly_array[count].rhoexp=0;
			}
			free(buff);
			 
 skip_lines(infd);


#ifdef DEBUG
   printf("input: boundary %d has %d edges, hollow=%d, 1/mu="MDF", rho="MDF"\n",
																			count,i,bound.poly_array[count].hollow,
																			bound.poly_array[count].mu,bound.poly_array[count].rho);
#endif
  /* initialize polygon */
	init_poly(&bound.poly_array[count]);
/* read edge numbers per line */
	for (; i; i--){
    if ( fscanf(infd,"%d",&j) == EOF ) {
      fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
			exit(1);
    }
	e.p1 =edge_array[j].p1; e.p2=edge_array[j].p2; 
	e.type=edge_array[j].type;
	insert_edge_to_poly(&bound.poly_array[count],&e);
#ifdef DEBUG
   printf("input: %d of %d is edge %d :(%d  %d)\n",i,bound.poly_array[count].nedges,j,e.p1,e.p2);
#endif

	}
	/* arrange polygon edges */
#ifdef DEBUG
	printf("input: rearranging polygon %d\n",count);
#endif
 arrange_polgon_edges(&bound.poly_array[count]);
}

   fclose(infd);
   return(npoint+3);
}


/* function to get the expression for rho from input and build parse tree */
/* here is how its called: 1: equation_parse(buffer)->call yyparse() with buffer
 * as input->parses buffer and builds tree->calls insert_expression to 
 * attach tree to current boundary globally */
int
insert_expression(exp_nodetype *p, MY_INT boundary_no)
{
				      bound.poly_array[boundary_no].rhoexp = p;
							return(0);
}