Category : Science and Education
Archive   : SURF1311.ZIP
Filename : FEM2SURF.BAS

10 ' ***************************************
20 ' ************ READFEM ******************
30 ' ***************************************
40 ' Read Finite Element output data and convert to format for SURFMODL.
50 ' Copyright 1986 by Kenneth Van Camp
60 '
80 DEFINT I-N
100 MXNODE=1000: MXSURF=500: MXELEM=1000: NN=8: NC=4
120 DIM XYZ(MXNODE,3), IELEM(MXELEM,NN), NFLAG(MXNODE), N(MXSURF,NC), ISURF(6,4)
140 'INITIALIZATION
150 FOR I=1 TO 6
160 FOR J=1 TO 4
170 READ ISURF(I,J)
180 NEXT J
190 NEXT I
200 ' DATA FOR SURFACES ON EACH ELEMENT
210 DATA 1,2,3,4, 8,7,6,5, 2,6,7,3, 3,7,8,4, 4,8,5,1, 1,5,6,2
220 '
230 INPUT "Enter file name for FEM data"; IFILE$
240 'Convert first char to upper case
250 OPEN "I", #1, IFILE$
260 '
270 ' Read the finite element data
280 '
290 INPUT #1, NNODE, NELEM
300 FOR I=1 TO NNODE
310 INPUT #1, XYZ(I,1), XYZ(I,2), XYZ(I,3)
320 NEXT I
330 FOR I=1 TO NELEM
340 INPUT #1, IELEM(I,1), IELEM(I,2), IELEM(I,3), IELEM(I,4), IELEM(I,5), IELEM(I,6), IELEM(I,7), IELEM(I,8)
350 NEXT I
360 CLOSE #1
370 '
380 PRINT "Data for";NNODE;"nodes and";NELEM;"elements read."
390 INPUT "Do you wish to make a slice thru the solid [Y or N]"; A$
400 IF A$="N" OR A$="n" THEN 540
410 INPUT "Enter 1, 2 or 3 to cut X, Y or Z plane respectively"; IPLANE
420 IF IPLANE<0 OR IPLANE>3 THEN 390
430 INPUT "Cutoff value"; CUTOFF
440 INPUT "Enter 1 to keep nodes > cutoff or -1 to keep nodes < cutoff"; IORIENT
450 IF IORIENT<>1 AND IORIENT<>-1 THEN 440
460 FOR I=1 TO NNODE
470 IF IORIENT < 0 THEN 500
480 IF XYZ(I,IPLANE) >= CUTOFF THEN COUNT(I)=0 ELSE COUNT(I)=8
490 GOTO 520
500 ' Else
510 IF XYZ(I,IPLANE) <= CUTOFF THEN COUNT(I)=0 ELSE COUNT(I)=8
520 ' Endif
530 NEXT I
540 'Endif
550 PRINT "Finding surface nodes ..."
560 '
570 'Count the #elements that share each node. If a node is shared by 8
580 ' or more elements, it is assumed to be an interior node. If it is
590 ' shared by fewer than 8 elements, it is assumed to be exterior.
600 ' An exterior surface is one that has all four of its corner nodes
610 ' on the exterior.
620 '
630 FOR I=1 TO NELEM
640 FOR J=1 TO 8
650 ' Warning: Degenerate elements will cause a node to be incorrectly
660 ' considered an interior node!
670 COUNT(IELEM(I,J)) = COUNT(IELEM(I,J)) + 1
680 NEXT J
690 NEXT I
700 '
710 PRINT "Finding exterior surfaces ..."
720 NSURF=0
730 FOR I=1 TO NELEM
740 FOR J=1 TO 6
750 FOR K=1 TO 4
760 IF COUNT(IELEM(I,ISURF(J,K))) > 7 THEN 830
770 NEXT K
780 ' This is an exterior surface
790 NSURF = NSURF + 1
800 FOR K=1 TO 4
810 N(NSURF,K) = ISURF(J,K)
820 NEXT K
830 NEXT J
840 NEXT I
850 '
860 'Count the surface nodes
870 '
880 NSNODE=0
890 FOR I=1 TO NNODE
900 IF COUNT(I) < 8 THEN NSNODE = NSNODE + 1
910 NEXT I
920 INPUT "Enter file name for SURFMODL data"; OFILE$
930 INPUT "Enter title for plot";TITL$
940 OPEN "O", #1, OFILE$
950 PRINT#1, TITL$
955 PRINT #1,4
960 PRINT #1,1;NSNODE;NSURF;4;1
970 PRINT #1, 1;0;0;3;.1
980 FOR I=1 TO NNODE
990 IF COUNT(I) < 8 THEN PRINT #1, XYZ(I,1); XYZ(I,2); XYZ(I,3)
1000 NEXT I
1020 FOR I=1 TO NSURF
1030 PRINT #1, 4; 1; N(I,1); N(I,2); N(I,3); N(I,4)
1040 NEXT I
1060 CLOSE #1
1065 PRINT "FINISHED CREATING ";OFILE$
1070 END
 PRINT #1, 4; 1; N(I,1); N(I,2); N(I,3); N(I,4)
1040 NEXT I
1060 CLOSE #1