1 /*  pdnmesh - a 2D finite element solver
2     Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>
3   This program is free software; you can redistribute it and/or modify
4   it under the terms of the GNU General Public License as published by
5   the Free Software Foundation; either version 2 of the License, or
6   (at your option) any later version.
7 
8   This program is distributed in the hope that it will be useful,
9   but WITHOUT ANY WARRANTY; without even the implied warranty of
10   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11   GNU General Public License for more details.
12 
13   You should have received a copy of the GNU General Public License
14   along with this program; if not, write to the Free Software
15   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
16   $Id: input.c,v 1.20 2005/04/26 05:16:47 sarod Exp $
17 */
18 
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22 
23 /* IO rutines */
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <ctype.h>
28 #include "types.h"
29 
30 poly_edge *edge_array;
31 int nedges;
32 boundary bound; /* boundaries */
33 MY_DOUBLE g_xoff,g_yoff; /* x,y offsets */
34 MY_DOUBLE g_xscale,g_yscale; /* x,y scales */
35 MY_DOUBLE g_xrange,g_yrange; /* x,y ranges, used for contour plot */
36 
37 
38 /* skips comment lines */
39 static int
skip_lines(FILE * fin)40 skip_lines(FILE *fin)
41 {
42 
43 		int c;
44 		do {
45 		if ( ( c = getc(fin) ) == EOF )
46 		return(-1);
47 	/* handle empty lines */
48 	if ( c == '\n' )
49 	continue; /* next line */
50 	if ( (c != '#') ) {
51 	ungetc(c,fin);
52 	return(0);
53 	} else { /* skip this line */
54 	do {
55 	if ( ( c = getc(fin) ) == EOF )
56 	return(-1);
57 	} while (  c != '\n') ;
58 	}
59 	} while( 1 );
60 }
61 
62 /* skips rest of line */
63 static int
skip_restof_line(FILE * fin)64 skip_restof_line(FILE *fin)
65 {
66 	int c;
67 	do {
68 	if ( ( c = getc(fin) ) == EOF )
69 	return(-1);
70 	} while (  c != '\n') ;
71 	return(1);
72 }
73 
74 /* reads the next string (isalphanumeric() contiguous set of characters)
75   separated by spaces, tabs or a newline. If the last character read is newline
76   1 is returned, else 0 returned. */
77 /* buffer is automatically adjusted is length is not enough */
78 static int
read_next_string(char * buff,int * buff_len,FILE * infd)79 read_next_string(char *buff, int *buff_len, FILE *infd) {
80    int k,c,flag;
81    k = 0;
82    /* intialize buffer */
83    buff[0]='\0';
84    /* skip leading white space */
85    do {
86    c=fgetc(infd);
87    if(c=='\n' || c==EOF) return 1;
88    } while(c != EOF && isblank(c));
89    if(c=='\n' || c==EOF) return 1;
90    /* now we have read a non whitespace character */
91    buff[k++]=c;
92   if (k==*buff_len) {
93     /* now we have run out of buffer */
94     *buff_len += 30;
95     if ((buff = (char*)realloc((void*)buff,sizeof(char)*(size_t)(*buff_len)))==NULL) {
96   fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
97   exit(1);
98   }
99   }
100    flag=0;
101    while ( ((c = fgetc(infd)) != EOF ) && k < *buff_len) {
102       if ( c == '\n') {  flag=1; break; }
103       if ( isblank(c) ) {  break; }/* not end of line */
104       buff[k++] = c;
105       if (k==*buff_len) {
106        /* now we have run out of buffer */
107        *buff_len += 30;
108        if((buff = (char*)realloc((void*)buff,sizeof(char)*(size_t)(*buff_len)))==NULL) {
109     fprintf(stderr,"%s: %d: no free memory\n",__FILE__,__LINE__);
110     exit(1);
111        }
112       }
113   }
114   /* now c == blank , \n or EOF */
115   if (k==*buff_len-2) {
116     /* now we have run out of buffer */
117     *buff_len += 3;
118     if((buff = (char*)realloc((void*)buff,sizeof(char)*(size_t)(*buff_len)))==NULL) {
119     fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
120     exit(1);
121    }
122   }
123 
124   /* add trailing ';' to buffer and '\0' */
125   buff[k++]=';';
126   buff[k++]='\0';
127   return flag;
128 }
129 
130 
131 /* read input coordinate file */
132 int
read_input_file(point ** input_points,const char * infile)133 read_input_file(point **input_points, const char *infile)
134 {
135  /* read points from infile and store coords in x,y, arrays */
136  FILE *infd;
137  int i,j,count;
138 	MY_DOUBLE xmax,xmin,ymax,ymin;
139 	int npoint, input_degree_of_freedom=1;
140  poly_edge e;
141 int buff_len,c;
142 char *buff;
143  exp_nodetype *pp=0; /* pointer to get root of parse tree from lex/yacc */
144 
145  infd =fopen( infile, "r" );
146  if ( infd == NULL ) {
147    fprintf(stderr, "%s: %d: could not open file %s\n", __FILE__,__LINE__,infile);
148    exit(1);
149  }
150 	skip_lines(infd);
151 	/* numer of points */
152  if ( fscanf(infd, "%d" , &npoint ) == EOF ) {
153     fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
154     exit(1);
155  }
156 	skip_restof_line(infd);
157 	skip_lines(infd);
158 #ifdef DEBUG
159   printf("input: reading %d points\n",npoint);
160 #endif
161 
162  /* allocate memory for points */
163  if ((*input_points=(point*)malloc((size_t)(npoint+3)*sizeof(point)))==0) {
164    fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
165    exit(1);
166  }
167 
168 	xmax=-1000;xmin=1000;
169 	ymax=-1000;ymin=1000;
170 
171  for (i=3; i<npoint+3; i++ ) {
172 		/* allocate memory for z array */
173 	if (input_degree_of_freedom) {
174   if (((*input_points)[i].z=(MY_DOUBLE*)malloc((size_t)(input_degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
175    fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
176    exit(1);
177  }
178 	}
179 	/* read point data according to degree of freedom */
180 	 if ( input_degree_of_freedom==0 ) { /* wave equation, do not need this */
181     if ( fscanf(infd,"%d "MDF" "MDF"  %d",&count, &((*input_points)[i].x),&((*input_points)[i].y),&j) == EOF ) {
182       fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
183     }
184    } else if ( input_degree_of_freedom==1 ) { /* potentials */
185     if ( fscanf(infd,"%d "MDF" "MDF" "MDF" %d",&count, &((*input_points)[i].x),&((*input_points)[i].y),&((*input_points)[i].z[0]),&j) == EOF ) {
186       fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
187     }
188   } else { /* input_degree_of_freedom==2, elasticity problems etc. */
189     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 ) {
190       fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
191     }
192   }
193 
194 #ifdef DEBUG
195   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);
196 	if ( input_degree_of_freedom==0 ) {
197   printf("no z values\n");
198   } else if (input_degree_of_freedom==1) {
199    printf("z[0]="MDF", j=%d\n",(*input_points)[i].z[0],j);
200 	} else {
201    printf("z[0]="MDF", z[1]="MDF",j=%d\n",(*input_points)[i].z[0], (*input_points)[i].z[1], j);
202 	}
203 #endif
204  skip_restof_line(infd);
205  skip_lines(infd);
206 
207  if ( xmin > (*input_points)[i].x ) {
208 								xmin=(*input_points)[i].x;
209 	}
210 	if ( xmax < (*input_points)[i].x ) {
211 								xmax=(*input_points)[i].x;
212 	}
213  if ( ymin > (*input_points)[i].y ) {
214 								ymin=(*input_points)[i].y;
215 	}
216 	if ( ymax < (*input_points)[i].y ) {
217 								ymax=(*input_points)[i].y;
218 	}
219 /* point type */
220 	if (j==1) {
221 					/* fixed point */
222 					(*input_points)[i].val=FX;
223 	} else {
224 					/* variable point */
225 					(*input_points)[i].val=VAR;
226 	}
227  }
228  g_xoff=-(xmax+xmin)*0.5;
229  g_yoff=-(ymax+ymin)*0.5;
230  g_xrange=(xmax-xmin)*0.5;
231  g_yrange=(ymax-ymin)*0.5;
232 
233 	 /* now find the maximum value of coordinates */
234   ymax = ABS(ymax);
235 		xmax = ABS(xmax);
236 		xmin = ABS(xmin);
237 		ymin = ABS(ymin);
238 		if ( xmax < ymax )
239 			xmax = ymax;
240 		if ( xmax < ymin )
241 			xmax = ymin;
242 		if ( xmax < xmin  )
243 			xmax = xmin;
244   g_xscale=xmax;
245   g_yscale=xmax;
246 		/* first 3 points (-3M,-3M),(3M,0),(0,3M) */
247 		 for (i=3;i<npoint+3;i++) {
248           (*input_points)[i].x+=(g_xoff);
249           (*input_points)[i].y+=(g_yoff);
250 
251           (*input_points)[i].x/=(xmax);
252           (*input_points)[i].y/=(xmax);
253 		 }
254 			(*input_points)[0].x=-3;(*input_points)[0].y=-3;
255 			(*input_points)[1].x=3;(*input_points)[1].y=0;
256 			(*input_points)[2].x=0;(*input_points)[2].y=3;
257 
258 
259 	/* now the edges */
260 	/* first number of edges */
261    if ( fscanf(infd, "%d" , &nedges) == EOF ) {
262     fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
263     exit(1);
264  }
265  skip_restof_line(infd);
266 	skip_lines(infd);
267 
268 
269 #ifdef DEBUG
270   printf("input: reading %d edges\n",nedges);
271 #endif
272 
273  /* allocate memory for edge array */
274  if ((edge_array =(poly_edge *)malloc((size_t)nedges*sizeof(poly_edge)))==0){
275    fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
276    exit(1);
277  }
278 
279  for (i=0; i< nedges; i++) {
280  /* FIXME - edges need not be in number order */
281   if ( fscanf(infd, "%d %d %d %d" , &count, &(edge_array[i].p1), &(edge_array[i].p2), &j) == EOF ) {
282 	      fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
283 				exit(1);
284 		}
285  skip_restof_line(infd);
286 	skip_lines(infd);
287 
288 
289 		/* account for the 3 outer points, so increment point numbers */
290    edge_array[i].p1+=3;edge_array[i].p2+=3;
291 			/* dirichlet or neumann ? */
292 			if (j==1) {
293 				edge_array[i].type=DR;
294 			} else {
295 				edge_array[i].type=NU;
296 			}
297 #ifdef DEBUG
298   printf("input: edge %d (%d,%d) ?%d\n",i,edge_array[i].p1,edge_array[i].p2,j);
299 #endif
300 
301  }
302 
303 	/* read the boundaries */
304  /* first number of boundaries */
305    if ( fscanf(infd, "%d" , &bound.npoly) == EOF ) {
306     fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
307     exit(1);
308  }
309  skip_restof_line(infd);
310 	skip_lines(infd);
311 
312 
313 	/* allocate memory for boundaries */
314  if ((bound.poly_array=(polygon*)malloc((size_t)bound.npoly*sizeof(polygon)))==0){
315    fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
316    exit(1);
317  }
318 
319 /* setup buffer to read expressions */
320  buff_len = 32;
321  if((buff = (char*)malloc(sizeof(char)*(size_t)(buff_len)))==NULL) {
322         fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
323         exit(1);
324  }
325 
326 for (count=0; count< bound.npoly; count++) {
327   skip_lines(infd);
328   /* initialize polygon */
329 	init_poly(&bound.poly_array[count]);
330 
331 	/* now read the outer boundary */
332  /* boundary format :
333    number of edges, hollow?, mu, epsilon, rho
334     OR
335    number of edge, hollow?, 1/mu or epsilon, rho - old format */
336 	/* first, number of edges, hollow?*/
337 	if ( fscanf(infd, "%d %d", &i,
338 							&bound.poly_array[count].hollow) == EOF ) {
339     fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
340     exit(1);
341  }
342 
343  /* next read mu */
344  memset(buff,0,buff_len);
345  bound.poly_array[count].muexp=0;
346  c=read_next_string(buff,&buff_len,infd);
347 #ifdef DEBUG
348  printf("%d, mu=%s\n",count, buff);
349 #endif
350  /* if c==1, '\n' was detected and we cannot read any more */
351  if(c != 1) {
352   /* parse buffer and extract mu */
353   /* mu is an expression */
354    pp=0;
355    bound.poly_array[count].mu=1.0;
356 	 equation_parse(buff,&pp);
357 	 bound.poly_array[count].muexp= pp;
358    if((bound.poly_array[count].mustr
359        =(char*)malloc(sizeof(char)*(size_t)(strlen(buff)+1)))==NULL) {
360        fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
361        exit(1);
362    }
363    strcpy(bound.poly_array[count].mustr, buff);
364  /* if we are reading 1/mu get the reciprocal (mu cannot be below 1) */
365 
366  /* next read epsilon or rho. we do not know this at this point,
367   but we assume it is epsilon */
368  memset(buff,0,buff_len);
369  bound.poly_array[count].epsexp=0;
370  c=read_next_string(buff,&buff_len,infd);
371 #ifdef DEBUG
372  printf("%d, eps=%s\n",count, buff);
373 #endif
374   /* parse buffer and extract mu */
375   /* mu is an expression */
376    pp=0;
377    bound.poly_array[count].epsilon=1.0;
378 	 equation_parse(buff,&pp);
379 	 bound.poly_array[count].epsexp= pp;
380    if((bound.poly_array[count].epsstr=(char*)malloc(sizeof(char)*(size_t)(strlen(buff)+1)))==NULL) {
381      fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
382      exit(1);
383    }
384    strcpy(bound.poly_array[count].epsstr, buff);
385    /* if c==1, '\n' was detected and we cannot read any more */
386    if(c != 1) {
387 
388   /* now read the expression for rho */
389     memset(buff,0,buff_len);
390     bound.poly_array[count].rhoexp=0;
391     c=read_next_string(buff,&buff_len,infd);
392 #ifdef DEBUG
393  printf("%d, rho=%s\n",count,buff);
394 #endif
395         pp=0;
396         bound.poly_array[count].rho=0;
397 				equation_parse(buff,&pp);
398 				bound.poly_array[count].rhoexp= pp;
399        if((bound.poly_array[count].rhostr=(char*)malloc(sizeof(char)*(size_t)(strlen(buff)+1)))==NULL) {
400          fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
401         exit(1);
402        }
403        strcpy(bound.poly_array[count].rhostr, buff);
404  } else {
405    /* c==1 means end of line. so we have read rho, NOT epsilon */
406    /* copy back everything to rho */
407    bound.poly_array[count].rho=bound.poly_array[count].epsilon;
408    bound.poly_array[count].rhoexp=bound.poly_array[count].epsexp;
409    bound.poly_array[count].rhostr=bound.poly_array[count].epsstr;
410 
411    bound.poly_array[count].epsilon=1.0;
412    bound.poly_array[count].epsexp=0;
413    /*bound.poly_array[count].epsstr="1.0;"; */
414    if((bound.poly_array[count].epsstr=(char*)malloc(sizeof(char)*(size_t)(5)))==NULL) {
415      fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
416      exit(1);
417    }
418    strcpy(bound.poly_array[count].epsstr, "1.0;");
419 
420 
421  }
422  } else {
423     fprintf(stderr,"%s: %d: invalid numeric value for mu %s\n",__FILE__,__LINE__,buff);
424    bound.poly_array[count].mu=1.0;
425    bound.poly_array[count].epsilon=1.0;
426    bound.poly_array[count].rho=0.0;
427 
428    bound.poly_array[count].muexp=0;
429    bound.poly_array[count].epsexp=0;
430    bound.poly_array[count].rhoexp=0;
431 
432    /*bound.poly_array[count].mustr="1.0;"; */
433    if((bound.poly_array[count].mustr=(char*)malloc(sizeof(char)*(size_t)(5)))==NULL) {
434      fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
435      exit(1);
436    }
437    strcpy(bound.poly_array[count].mustr, "1.0;");
438    /*bound.poly_array[count].epsstr="1.0;"; */
439    if((bound.poly_array[count].epsstr=(char*)malloc(sizeof(char)*(size_t)(5)))==NULL) {
440      fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
441      exit(1);
442    }
443    strcpy(bound.poly_array[count].epsstr, "1.0;");
444    /*bound.poly_array[count].rhostr="1.0;"; */
445    if((bound.poly_array[count].rhostr=(char*)malloc(sizeof(char)*(size_t)(5)))==NULL) {
446      fprintf(stderr,"%s: %d: No free memory\n",__FILE__,__LINE__);
447      exit(1);
448    }
449    strcpy(bound.poly_array[count].rhostr, "1.0;");
450  }
451 #ifdef DEBUG
452  printf("boundary %d, mu=%s, epsilon=%s, rho=%s\n",count, bound.poly_array[count].mustr,
453   bound.poly_array[count].epsstr,bound.poly_array[count].rhostr);
454 #endif
455 
456  skip_lines(infd);
457 
458 #ifdef DEBUG
459    printf("input: boundary %d has %d edges, hollow=%d, 1/mu="MDF", rho="MDF"\n",
460 																			count,i,bound.poly_array[count].hollow,
461 																			bound.poly_array[count].mu,bound.poly_array[count].rho);
462 #endif
463 /* read edge numbers per line */
464 	for (; i; i--){
465     if ( fscanf(infd,"%d",&j) == EOF ) {
466       fprintf(stderr, "%s: %d: could not read file %s\n", __FILE__,__LINE__,infile);
467 			exit(1);
468     }
469 	e.p1 =edge_array[j].p1; e.p2=edge_array[j].p2;
470 	e.type=edge_array[j].type;
471 	insert_edge_to_poly(&bound.poly_array[count],&e);
472 #ifdef DEBUG
473    printf("input: %d of %d is edge %d :(%d  %d)\n",i,bound.poly_array[count].nedges,j,e.p1,e.p2);
474 #endif
475 
476 	}
477 	/* arrange polygon edges */
478 #ifdef DEBUG
479 	printf("input: rearranging polygon %d\n",count);
480 #endif
481  arrange_polgon_edges(&bound.poly_array[count]);
482 }
483 
484 	 free(buff);
485    fclose(infd);
486    return(npoint+3);
487 }
488 
489 
490 /* function to get the expression for rho from input and build parse tree */
491 /* here is how its called: 1: equation_parse(buffer)->call yyparse() with buffer
492  * as input->parses buffer and builds tree->calls insert_expression to
493  * attach tree to current boundary globally */
494 /*int
495 insert_expression(exp_nodetype *p, MY_INT boundary_no)
496 {
497 				      bound.poly_array[boundary_no].rhoexp = p;
498 							return(0);
499 }*/
500