1 /*****************************************************************************
2 *
3 * Elmer, A Finite Element Software for Multiphysical Problems
4 *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program (in file fem/GPL-2); if not, write to the
19 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
20 * Boston, MA 02110-1301, USA.
21 *
22 *****************************************************************************/
23
24 /*******************************************************************************
25 *
26 * File reading + MATC model setting command
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02101 Espoo, Finland
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 27 Sep 1995
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 ******************************************************************************/
46
47
48 #include "../elmerpost.h"
49 #include <tcl.h>
50
51 extern double XMin,XMax,YMin,YMax,ZMin,ZMax;
52 extern double xmin,xmax,ymin,ymax,zmin,zmax;
53 extern int CurrentTimeStep,KeepScale;
54 static int E[50],code,NV,NE,NF,NT;
55
56 #define BUFFER_SIZE 8192
57
epReadFile(ClientData cl,Tcl_Interp * interp,int argc,char ** argv)58 static int epReadFile( ClientData cl,Tcl_Interp *interp,int argc,char **argv )
59 {
60 int i,j,k,n,t = 0,total = 0,NamesGiven;
61
62 FILE *fp;
63
64 static vertex_t vertex;
65 element_type_t *EL;
66
67 static char *str,name[512],*ptr;
68
69 double s,*NodeArray,*Velo,*Vabs,*Temp,*Pres,fdummy;
70
71 double *Vector[1000], *Scalar[1000], *Times,*tvar;
72
73 VARIABLE *Var1;
74
75 group_t *grp;
76
77 struct {
78 int type;
79 char name[256];
80 } variable[64];
81
82 int groupid,gid,StartTime=1,EndTime=1,IncTime=1,ToRead,last,gotit;
83 static char groupname[128];
84
85 group_t *group;
86 extern Tcl_Interp *TCLInterp;
87
88 if ( argc < 2 ) return TCL_ERROR;
89
90 fp = fopen( argv[1], "r" );
91 if ( !fp )
92 {
93 sprintf( interp->result, "ReadModel: can't open file [%s]\n",argv[1] );
94 return TCL_ERROR;
95 }
96
97 if ( argc > 2 ) StartTime=atoi(argv[2]);
98 if ( argc > 3 ) EndTime=atoi(argv[3]);
99 if ( argc > 4 ) IncTime=atoi(argv[4]);
100
101 str = (char *)malloc( BUFFER_SIZE*sizeof(char) );
102 if ( !str ) {
103 fprintf( stderr, "ERROR: ElmerPost: memory allocation error.\n" );
104 exit(0);
105 }
106
107 NV = NE = NT = 0;
108 fgets( str, BUFFER_SIZE-1, fp );
109 sscanf( str, "%d %d %d %d", &NV,&NE,&NF,&NT );
110
111 if ( NV <= 0 || NE <=0 )
112 {
113 Tcl_SetResult( interp, "Bad element model file.\n",TCL_STATIC );
114 return TCL_ERROR;
115 }
116
117 ptr = str;
118 for( i=0; i<4; i++ )
119 {
120 while( *ptr && isspace(*ptr) ) ptr++;
121 while( *ptr && !isspace(*ptr) ) ptr++;
122 }
123 while( *ptr && isspace(*ptr) ) ptr++;
124
125 i = 0;
126 while ( *ptr )
127 {
128 if ( sscanf( ptr, "vector:%s", name ) == 1 )
129 {
130 variable[i].type = 1;
131 #if 0
132 variable[i].name = (char *)malloc( strlen(name)+1 );
133 #endif
134 strcpy( variable[i].name, name );
135
136 i++;
137 } else if ( sscanf( ptr, "scalar:%s", name ) == 1 )
138 {
139 variable[i].type = 2;
140 #if 0
141 variable[i].name = (char *)malloc( strlen(name)+1 );
142 #endif
143 strcpy( variable[i].name, name );
144
145 i++;
146 }
147
148 while( *ptr && !isspace(*ptr) ) ptr++;
149 while( *ptr && isspace(*ptr) ) ptr++;
150
151 while( *ptr && !isspace(*ptr) ) ptr++;
152 while( *ptr && isspace(*ptr) ) ptr++;
153 }
154 NamesGiven = i;
155
156 if ( !CurrentObject->ElementModel )
157 CurrentObject->ElementModel = (element_model_t *)calloc( 1,sizeof(element_model_t) );
158
159 if ( CurrentObject == &VisualObject )
160 CurrentObject->ElementModel->NodeArray = NodeArray =
161 MATR( var_new( "nodes", TYPE_DOUBLE, 3, NV ) );
162 else
163 {
164 sprintf( str, "nodes_%s", CurrentObject->Name );
165
166 CurrentObject->ElementModel->NodeArray = NodeArray =
167 MATR( var_new( str, TYPE_DOUBLE, 3, NV ) );
168 }
169
170 Tcl_LinkVar( TCLInterp, "NumberOfTimesteps", (char *)&CurrentObject->ElementModel->NofTimesteps, TCL_LINK_INT );
171
172 xmin = ymin = zmin = DBL_MAX;
173 xmax = ymax = zmax = -DBL_MAX;
174 for( i=0; i<NV; i++ )
175 {
176 fgets(str,BUFFER_SIZE-1,fp);
177 if ( *str == '#' ) { i--; continue; }
178
179 sscanf( str, "%lf %lf %lf", &NodeArray[i],&NodeArray[NV+i],&NodeArray[2*NV+i] );
180
181 xmin = MIN( xmin, NodeArray[i] );
182 ymin = MIN( ymin, NodeArray[NV+i] );
183 zmin = MIN( zmin, NodeArray[2*NV+i] );
184
185 xmax = MAX( xmax, NodeArray[i] );
186 ymax = MAX( ymax, NodeArray[NV+i] );
187 zmax = MAX( zmax, NodeArray[2*NV+i] );
188 }
189
190 if (CurrentObject->ElementModel->Elements )
191 {
192 for( i=0; i<CurrentObject->ElementModel->NofElements; i++ )
193 if (CurrentObject->ElementModel->Elements[i].Topology )
194 {
195 free( CurrentObject->ElementModel->Elements[i].Topology );
196 }
197 free(CurrentObject->ElementModel->Elements );
198 }
199 CurrentObject->ElementModel->Elements = Elements = (element_t *)calloc( NE,sizeof(element_t) );
200
201 geo_free_groups(CurrentObject->ElementModel->Groups );
202 CurrentObject->ElementModel->Groups = NULL;
203
204 for( i=0; i<NE; i++ )
205 {
206 fgets(str,BUFFER_SIZE-1,fp);
207
208 if ( *str == '#' )
209 {
210 if ( strncmp( str, "#group ", 7 ) == 0 )
211 {
212 sscanf( str, "#group %s", groupname );
213 groupid = geo_group_id( CurrentObject->ElementModel,groupname,1 );
214 } else if ( strncmp( str, "#endgroup ",10 ) == 0 )
215 {
216 sscanf( str, "#endgroup %s", groupname );
217 groupid = geo_group_id( CurrentObject->ElementModel,groupname,0 );
218 }
219 i--;
220 continue;
221 }
222
223 for( gid=0; gid<MAX_GROUP_IDS; gid++ )CurrentObject->ElementModel->Elements[i].GroupIds[gid] = -1;
224
225 n = sscanf( str, "%s %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d",
226 groupname, &code, &E[0],&E[1],&E[2],&E[3],&E[4],&E[5],&E[6],&E[7],&E[8],&E[9],
227 &E[10],&E[11],&E[12],&E[13],&E[14],&E[15],&E[16],&E[17],&E[18],&E[19],
228 &E[20],&E[21],&E[22],&E[23],&E[24],&E[25],&E[26] );
229
230 n = n - 2;
231
232 while( n < (code - 100 * (code / 100)) )
233 {
234 fgets( str,BUFFER_SIZE-1,fp );
235 n = n + sscanf( str, "%d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d %d",
236 &E[n],&E[n+1],&E[n+2],&E[n+3],&E[n+4],&E[n+5],&E[n+6],&E[n+7],&E[n+8],&E[n+9],
237 &E[n+10],&E[n+11],&E[n+12],&E[n+13],&E[n+14],&E[n+15],&E[n+16],&E[n+17],&E[n+18],
238 &E[n+19],&E[n+20],&E[n+21],&E[n+22],&E[n+23],&E[n+24],&E[n+25],&E[n+26] );
239 }
240
241
242 groupid = geo_group_id( CurrentObject->ElementModel,groupname,1 );
243
244 for( EL=ElementDefs.ElementTypes; EL!=NULL; EL=EL->Next )
245 if ( code == EL->ElementCode )
246 {
247 Elements[total].ElementType = EL;
248 Elements[total].DisplayFlag = TRUE;
249
250 for( gid=0; gid<MAX_GROUP_IDS; gid++ ) if ( Elements[total].GroupIds[gid] < 0 ) break;
251
252 groupid = 0;
253 for( grp=CurrentObject->ElementModel->Groups; grp != NULL; grp=grp->Next, groupid++ )
254 if ( grp->Open ) if ( gid < MAX_GROUP_IDS ) Elements[total].GroupIds[gid++] = groupid;
255
256 Elements[total].Topology = (int *)malloc( EL->NumberOfNodes*sizeof(int) );
257
258 if ( !Elements[total].Topology )
259 {
260 Tcl_SetResult( interp,"FATAL: can't alloc element connection tables.\n",TCL_STATIC );
261 return TCL_ERROR;
262 }
263
264 for( j=0; j<EL->NumberOfNodes; j++ ) Elements[total].Topology[j] = E[j];
265
266 total++;
267
268 break;
269 }
270
271 groupid = geo_group_id( CurrentObject->ElementModel,groupname,0 );
272
273 if ( EL == NULL )
274 {
275 if ( code != 101 )
276 fprintf( stderr,"Unknown element type: [%d]. Skipping Element.\n", code );
277 }
278 }
279
280 CurrentObject->ElementModel->NofElements = NE = total;
281 CurrentObject->ElementModel->NofNodes = NV;
282 CurrentObject->ElementModel->NofTimesteps = (EndTime-StartTime+IncTime) / IncTime;
283 ToRead = CurrentObject->ElementModel->NofTimesteps;
284
285 last = 0;
286 t = 0;
287 if ( NF>0 )
288 {
289 Times = malloc( 3*ToRead*sizeof(double) );
290 for( i=0; i<ToRead; i++ )
291 {
292 Times[0*ToRead+i] = i+1;
293 Times[1*ToRead+i] = i+1;
294 Times[2*ToRead+i] = i+1;
295 }
296
297 if ( !NamesGiven ) {
298 if ( NF>=3 )
299 {
300 Vabs = MATR( var_new( "vabs", TYPE_DOUBLE, 1, ToRead*NV ) );
301 Velo = MATR( var_new( "velo", TYPE_DOUBLE, 3, ToRead*NV ) );
302 }
303 Pres = MATR( var_new( "pres", TYPE_DOUBLE, 1, ToRead*NV ) );
304 Temp = MATR( var_new( "temp", TYPE_DOUBLE, 1, ToRead*NV ) );
305
306 if ( NF>=3 )
307 {
308 t = 0;
309 for( i=0; i<EndTime; i++ )
310 {
311 if ( feof(fp) ) goto exit_loop;
312
313 if ( i<StartTime-1 || ((i-StartTime+1)%IncTime) ) {
314 for( k=0; k<NV; k++ ) {
315 fgets( str,BUFFER_SIZE-1,fp );
316 if ( feof(fp) ) goto exit_loop;
317 if ( *str == '#' ) k--;
318 }
319 } else {
320 for( k=0; k<NV; k++ )
321 {
322 fgets( str,BUFFER_SIZE-1,fp );
323 if ( feof( fp ) ) goto exit_loop;
324 if ( *str == '#' ) { k--; continue; }
325 sscanf( str, "%lf %lf %lf %lf %lf",&Velo[t*NV+k], &Velo[NV*(t+ToRead)+k],
326 &Velo[NV*(t+2*ToRead)+k], &Pres[t*NV+k],&Temp[t*NV+k] );
327 }
328 t++;
329 }
330 }
331
332 for( i=0; i<NV*ToRead; i++ )
333 {
334 Vabs[i] = sqrt( Velo[i]*Velo[i]+Velo[NV*ToRead+i]*Velo[NV*ToRead+i]+
335 Velo[2*NV*ToRead+i]*Velo[2*NV*ToRead+i] );
336 }
337 } else {
338 t = 0; last = 0;
339 for( i=0; i<EndTime; i++ )
340 {
341 if ( feof(fp) ) goto exit_loop;
342 if ( i < StartTime-1 || ((i-StartTime+1)%IncTime) )
343 {
344 for( k=0; k<NV; k++ ) {
345 fgets( str,BUFFER_SIZE-1,fp );
346 if ( feof( fp ) ) goto exit_loop;
347 if ( *str == '#' ) k--;
348 }
349 } else {
350 for( k=0; k<NV; k++ )
351 {
352 fgets( str,BUFFER_SIZE-1,fp );
353 if ( feof( fp ) ) goto exit_loop;
354 if ( *str == '#' ) { k--; continue; }
355 sscanf( str, "%lf %lf", &Pres[t*NV+k],&Temp[t*NV+k] );
356 }
357 t++;
358 }
359 }
360 }
361 } else {
362 for( i=0; i<NamesGiven; i++ ) {
363 if ( variable[i].type == 1 ) {
364 Vector[i] = MATR( var_new( variable[i].name, TYPE_DOUBLE, 3, ToRead*NV ) );
365 strcpy( str, variable[i].name ); strcat( str, "_abs" );
366 Scalar[i] = MATR( var_new( str, TYPE_DOUBLE, 1, ToRead*NV ) );
367 } else {
368 Scalar[i] = MATR( var_new( variable[i].name, TYPE_DOUBLE, 1, ToRead*NV ) );
369 }
370 }
371
372 t = 0;
373 for( i=0; i<EndTime; i++ )
374 {
375 #if 0
376 if ( feof(fp) ) goto exit_loop;
377 if ( i<StartTime-1 || ((i-StartTime+1)%IncTime) )
378 {
379 for( j=0; j<NV*NF; j++ )
380 {
381 if ( feof(fp) ) goto exit_loop;
382 while ( 1 != fscanf( fp, "%lf", &fdummy ) )
383 {
384 fgets( str, BUFFER_SIZE-1, fp );
385 if ( feof( fp ) ) goto exit_loop;
386 if ( strncmp( str, "#time ", 6 ) == 0 )
387 {
388 double t1,t2,t3;
389
390 sscanf( str, "#time %lf %lf %lf", &t1,&t2,&t3 );
391 Times[0*ToRead+t] = t1;
392 Times[1*ToRead+t] = t2;
393 Times[2*ToRead+t] = t3;
394 }
395 }
396 }
397 } else
398 #endif
399 {
400 for( k=0; k<NV; k++ )
401 {
402 if ( feof( fp ) ) goto exit_loop;
403 for( j=0; j<NamesGiven; j++ )
404 {
405 if ( feof(fp) ) goto exit_loop;
406
407 if ( variable[j].type == 1 )
408 {
409 if ( feof(fp) ) goto exit_loop;
410 while ( 1 != fscanf( fp, "%lf",&Vector[j][t*NV+k] ) ) {
411 fgets( str,BUFFER_SIZE-1,fp );
412 if ( feof( fp ) ) goto exit_loop;
413 if ( strncmp( str, "#time ", 6 ) == 0 )
414 {
415 double t1,t2,t3;
416
417 sscanf( str, "#time %lf %lf %lf", &t1,&t2,&t3 );
418 Times[0*ToRead+t] = t1;
419 Times[1*ToRead+t] = t2;
420 Times[2*ToRead+t] = t3;
421 }
422 }
423
424 if ( feof(fp) ) goto exit_loop;
425 while ( 1 != fscanf( fp, "%lf",&Vector[j][(t+ToRead)*NV+k] ) ) {
426 fgets( str,BUFFER_SIZE-1,fp );
427 if ( feof(fp) ) goto exit_loop;
428 if ( strncmp( str, "#time ", 6 ) == 0 )
429 {
430 double t1,t2,t3;
431
432 sscanf( str, "#time %lf %lf %lf", &t1,&t2,&t3 );
433 Times[0*ToRead+t] = t1;
434 Times[1*ToRead+t] = t2;
435 Times[2*ToRead+t] = t3;
436 }
437 }
438
439 if ( feof(fp) ) goto exit_loop;
440 while ( 1 != fscanf( fp, "%lf",&Vector[j][(t+2*ToRead)*NV+k] ) ) {
441 fgets( str,BUFFER_SIZE-1,fp );
442 if ( feof(fp) ) goto exit_loop;
443 if ( strncmp( str, "#time ", 6 ) == 0 )
444 {
445 double t1,t2,t3;
446
447 sscanf( str, "#time %lf %lf %lf", &t1,&t2,&t3 );
448 Times[0*ToRead+t] = t1;
449 Times[1*ToRead+t] = t2;
450 Times[2*ToRead+t] = t3;
451 }
452 }
453
454 Scalar[j][t*NV+k] = sqrt(
455 Vector[j][t*NV+k]*Vector[j][t*NV+k]+
456 Vector[j][NV*(t+ToRead)+k]*Vector[j][NV*(t+ToRead)+k] +
457 Vector[j][NV*(t+2*ToRead)+k]*Vector[j][NV*(t+2*ToRead)+k] );
458
459 } else {
460 if ( feof(fp) ) goto exit_loop;
461 while ( 1 != fscanf( fp, "%lf",&Scalar[j][t*NV+k] ) )
462 {
463 fgets( str,BUFFER_SIZE-1,fp );
464 if ( feof(fp) ) goto exit_loop;
465 if ( strncmp( str, "#time ", 6 ) == 0 )
466 {
467 double t1,t2,t3;
468
469 sscanf( str, "#time %lf %lf %lf", &t1,&t2,&t3 );
470 Times[0*ToRead+t] = t1;
471 Times[1*ToRead+t] = t2;
472 Times[2*ToRead+t] = t3;
473 }
474 }
475 }
476 }
477 }
478 last = 0;
479 if ( i>=StartTime-1 && !((i-StartTime+1)%IncTime) ) { last=1; t++; }
480 if ( t >= ToRead ) goto exit_loop;
481 }
482 }
483 }
484 }
485 exit_loop:
486
487 if ( t == 0 && ToRead > 0 && i > 0 ) t++;
488
489 fclose( fp );
490
491 if ( t > 0 ) {
492 CurrentObject->ElementModel->NofTimesteps = t;
493 if ( CurrentObject == &VisualObject )
494 tvar = MATR( var_new( "times", TYPE_DOUBLE, 3, t ) );
495 else
496 {
497 sprintf( str, "times_%s", CurrentObject->Name );
498 tvar = MATR( var_new( str, TYPE_DOUBLE, 3, t ) );
499 }
500
501 for( i=0; i<t; i++ ) {
502 tvar[0*t+i] = Times[0*ToRead+i];
503 tvar[1*t+i] = Times[1*ToRead+i];
504 tvar[2*t+i] = Times[2*ToRead+i];
505 }
506
507 free( Times );
508 }
509
510 if ( t != ToRead ) {
511 fprintf( stderr,"WARNING: ElmerPost: Not enough data for all timesteps"
512 " requested. REQUEST: %d, GOT: %d\n", ToRead, t );
513
514 for( i=0; i<NamesGiven; i++ )
515 {
516 if ( t > 0 ) {
517 if ( variable[i].type == 1 )
518 {
519 Var1 = lst_find( VARIABLES, variable[i].name );
520 #if 0
521 NCOL(Var1) = t*NV;
522 #endif
523
524 strcpy( str, variable[i].name ); strcat( str, "_abs" );
525 Var1 = lst_find( VARIABLES, str );
526 #if 0
527 NCOL(Var1) = t*NV;
528 #endif
529 } else {
530 Var1 = lst_find( VARIABLES, variable[i].name );
531 #if 0
532 NCOL(Var1) = t*NV;
533 #endif
534 }
535 } else {
536 if ( variable[i].type == 1 )
537 {
538 var_delete( variable[i].name );
539 strcpy( str, variable[i].name ); strcat( str, "_abs" );
540 var_delete( str );
541 } else {
542 var_delete( variable[i].name );
543 }
544 }
545 }
546 }
547
548 groupid = 0;
549 for( group =CurrentObject->ElementModel->Groups; group!=NULL; group=group->Next,groupid++ )
550 {
551 sprintf( name, "Groups(%d)", groupid );
552 Tcl_SetVar( TCLInterp, name, group->Name, TCL_GLOBAL_ONLY );
553 }
554
555 sprintf( name, "NumberOfGroups" );
556 sprintf( str, "%d", groupid );
557 Tcl_SetVar( TCLInterp, name, str, TCL_GLOBAL_ONLY );
558
559 if ( CurrentObject->Geometry )
560 {
561 geo_free_edge_tables( CurrentObject->Geometry );
562 geo_free_vertex_face_tables( CurrentObject->Geometry );
563 } else {
564 CurrentObject->Geometry = (geometry_t *)calloc( 1,sizeof(geometry_t) );
565 }
566
567 CurrentObject->Geometry->VertexCount = 0;
568
569 if ( !KeepScale ) {
570 s = MAX( MAX( xmax-xmin, ymax-ymin ), zmax-zmin );
571 } else {
572 s = CurrentObject->Geometry->Scale;
573 xmin = CurrentObject->Geometry->MinMax[0].x[0];
574 ymin = CurrentObject->Geometry->MinMax[0].x[1];
575 zmin = CurrentObject->Geometry->MinMax[0].x[2];
576
577 xmax = CurrentObject->Geometry->MinMax[1].x[0];
578 ymax = CurrentObject->Geometry->MinMax[1].x[1];
579 zmax = CurrentObject->Geometry->MinMax[1].x[2];
580 }
581 XMin = YMin = ZMin = DBL_MAX;
582 XMax = YMax = ZMax = -DBL_MAX;
583
584 for( i=0; i<NV; i++ )
585 {
586 vertex.x[0] = ( 2.0 * ( NodeArray[i] - xmin) - (xmax - xmin)) / s;
587 vertex.x[1] = ( 2.0 * ( NodeArray[NV+i] - ymin) - (ymax - ymin)) / s;
588 vertex.x[2] = ( 2.0 * ( NodeArray[2*NV+i] - zmin) - (zmax - zmin)) / s;
589
590 XMin = MIN( XMin,vertex.x[0] );
591 YMin = MIN( YMin,vertex.x[1] );
592 ZMin = MIN( ZMin,vertex.x[2] );
593
594 XMax = MAX( XMax,vertex.x[0] );
595 YMax = MAX( YMax,vertex.x[1] );
596 ZMax = MAX( ZMax,vertex.x[2] );
597
598 vertex.ElementModelNode = TRUE;
599 geo_add_vertex( CurrentObject->Geometry, &vertex );
600 }
601
602 CurrentObject->Geometry->Scale = s;
603 CurrentObject->Geometry->MinMax[0].x[0] = xmin;
604 CurrentObject->Geometry->MinMax[0].x[1] = ymin;
605 CurrentObject->Geometry->MinMax[0].x[2] = zmin;
606
607 CurrentObject->Geometry->MinMax[1].x[0] = xmax;
608 CurrentObject->Geometry->MinMax[1].x[1] = ymax;
609 CurrentObject->Geometry->MinMax[1].x[2] = zmax;
610
611 CurrentObject->Geometry->Edges = (edge_t *)calloc( CurrentObject->Geometry->VertexCount, sizeof(edge_t) );
612
613 if ( !CurrentObject->Geometry->Edges )
614 {
615 fprintf( stderr, "Can't alloc edge array.\n" );
616 exit( 0 );
617 }
618
619 CurrentObject->Geometry->TriangleCount = 0;
620
621 for( i=0; i<NE; i++ )
622 {
623 EL = CurrentObject->ElementModel->Elements[i].ElementType;
624 (*EL->Triangulate)( CurrentObject->Geometry,(element_t *)&Elements[i],(element_t *)&Elements[i] );
625 }
626
627 geo_vertex_normals( CurrentObject->Geometry,50.0 );
628
629 CurrentTimeStep = 0;
630
631 UpdateObject( 0,NULL,0,NULL );
632 DrawItSomeTimeWhenIdle();
633
634 free( str );
635
636 return TCL_OK;
637 }
638
639 #if 0
640 static VARIABLE *SetModel( VARIABLE *ptr )
641 {
642 static vertex_t vertex;
643 element_type_t *EL;
644
645 double s,*NodeArray,*Velo,*Vabs,*Temp,*Pres;
646
647 int i,j,k,total;
648
649 double *Topology = MATR(NEXT(ptr));
650 double *Type = MATR(NEXT(NEXT(ptr)));
651
652 NV = NCOL(ptr);
653
654 NE = NROW(NEXT(ptr));
655
656 if ( NEXT(NEXT(NEXT(ptr))) )
657 NT = M(NEXT(NEXT(NEXT(ptr))),0,0);
658 else
659 NT = 1;
660
661 CurrentObject->ElementModel->NodeArray = NodeArray = MATR( ptr );
662
663 xmin = ymin = zmin = DBL_MAX;
664 xmax = ymax = zmax = -DBL_MAX;
665 for( i=0; i<NV; i++ )
666 {
667 xmin = MIN( xmin, NodeArray[i] );
668 ymin = MIN( ymin, NodeArray[NV+i] );
669 zmin = MIN( zmin, NodeArray[2*NV+i] );
670
671 xmax = MAX( xmax, NodeArray[i] );
672 ymax = MAX( ymax, NodeArray[NV+i] );
673 zmax = MAX( zmax, NodeArray[2*NV+i] );
674 }
675
676 if (CurrentObject->ElementModel.Elements )
677 {
678 for( i=0; i<ElementModel.NofElements; i++ )
679 if (CurrentObject->ElementModel.Elements[i].Topology )
680 {
681 free( CurrentObject->ElementModel.Elements[i].Topology );
682 }
683 free(CurrentObject->ElementModel.Elements );
684 }
685 CurrentObject->ElementModel.Elements = Elements = (element_t *)calloc( NE,sizeof(element_t) );
686
687 total = 0;
688 for( i=0; i<NE; i++ )
689 {
690 for( EL=ElementDefs.ElementTypes; EL != NULL; EL=EL->Next )
691 if ( Type[i] == EL->ElementCode )
692 {
693 Elements[total].ElementType = EL;
694
695 Elements[total].Topology = (int *)malloc( EL->NumberOfNodes*sizeof(int) );
696
697 if ( !Elements[total].Topology )
698 {
699 error( "FATAL: can't alloc element connection tables.\n" );
700 }
701
702 for( j=0; j<EL->NumberOfNodes; j++ ) Elements[total].Topology[j] = Topology[i*NCOL(NEXT(ptr))+j];
703 total++;
704
705 break;
706 }
707
708 if ( EL == NULL )
709 {
710 fprintf( stderr, "Unknown element type: [%d]. Skipping Element.\n", Type[i] );
711 }
712 }
713
714 CurrentObject->ElementModel.NofElements = NE = total;
715 CurrentObject->ElementModel.NofNodes = NV;
716 CurrentObject->ElementModel.NofTimesteps = NT;
717
718 geo_free_edge_tables( &Geometry );
719 geo_free_vertex_face_tables( &Geometry );
720
721 Geometry.VertexCount = 0;
722
723 s = MAX( MAX( xmax-xmin, ymax-ymin ), zmax-zmin );
724 XMin = YMin = ZMin = DBL_MAX;
725 XMax = YMax = ZMax = -DBL_MAX;
726
727 for( i=0; i<NV; i++ )
728 {
729 vertex.x[0] = ( 2.0 * ( NodeArray[i] - xmin) - (xmax - xmin)) / s;
730 vertex.x[1] = ( 2.0 * ( NodeArray[NV+i] - ymin) - (ymax - ymin)) / s;
731 vertex.x[2] = ( 2.0 * ( NodeArray[2*NV+i] - zmin) - (zmax - zmin)) / s;
732
733 XMin = MIN( XMin,vertex.x[0] );
734 YMin = MIN( YMin,vertex.x[1] );
735 ZMin = MIN( ZMin,vertex.x[2] );
736
737 XMax = MAX( XMax,vertex.x[0] );
738 YMax = MAX( YMax,vertex.x[1] );
739 ZMax = MAX( ZMax,vertex.x[2] );
740
741 vertex.ElementModelNode = TRUE;
742 geo_add_vertex( &Geometry, &vertex );
743 }
744
745 Geometry.MinMax[0].x[0] = xmin;
746 Geometry.MinMax[0].x[1] = ymin;
747 Geometry.MinMax[0].x[2] = zmin;
748
749 Geometry.MinMax[1].x[0] = xmax;
750 Geometry.MinMax[1].x[1] = ymax;
751 Geometry.MinMax[1].x[2] = zmax;
752
753 Geometry.Edges = (edge_t *)calloc( Geometry.VertexCount, sizeof(edge_t) );
754
755 if ( !Geometry.Edges )
756 {
757 error( "SetModel: FATAL: Can't alloc edge array.\n" );
758 }
759
760 Geometry.TriangleCount = 0;
761
762 for( i=0; i<NE; i++ )
763 {
764 EL = Elements[i].ElementType;
765 (*EL->Triangulate)( &Geometry,(element_t *)&Elements[i],(element_t *)&Elements[i] );
766 }
767
768 geo_vertex_normals( &Geometry,50.0 );
769
770 CurrentTimeStep = 0;
771
772 fprintf( stderr, "TRIANGLES: %d %d, TS: %d\n", Geometry.TriangleCount,Geometry.MaxTriangleCount,NT );
773
774 UpdateObject( 0,NULL,0,NULL );
775 DrawItSomeTimeWhenIdle();
776
777 return (VARIABLE *)NULL;
778 }
779 #endif
780
Readfile_Init(Tcl_Interp * interp)781 int Readfile_Init( Tcl_Interp *interp )
782 {
783 Tcl_CreateCommand( interp,"cReadFile",epReadFile,(ClientData)NULL,(Tcl_CmdDeleteProc *)NULL);
784 #if 0
785 com_init
786 (
787 "model", FALSE, FALSE,SetModel,3,4,
788 "Usage: model(node-array,topology,element-type,[timesteps])\nSet current element model."
789 );
790 #endif
791
792 return TCL_OK;
793 }
794