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