C PROGRAM NUMBER C (VERSION OF 11 NOVEMBER 1999) C C BY PETER BIRD C DEPARTMENT OF EARTH AND SPACE SCIENCES C UNIVERSITY OF CALIFORNIA C LOS ANGELES, CA 90095-1567 C C C (C) COPYRIGHT 1993, 1995, 1997, 1999 BY PETER BIRD AND THE C REGENTS OF THE UNIVERSITY OF CALIFORNIA. C C PROGRAM TO READ A FINITE ELEMENT GRID OF 6-NODE ELEMENTS, C RENUMBER IT FOR GREATER EFFICIENCY, AND WRITE THE RESULT. C ACCEPTS EITHER -FAULTS- OR -PLATES- FORMATS. C C *** NOTE: SET VARIABLE "REDEEM" TO TRUE OR FALSE. C (IN INTERACTIVE VERSION OF PROGRAM, THIS IS DONE BY SETTING C INTEGER VARIABLE IMODE TO 1 OR 2.) C IF TRUE, IT MEANS THAT FAKE NODES (BOUNDARY NODES WITH NO C MECHANICAL DEGREES OF FREEDOM) WILL BE CONVERTED TO REAL C NODES DURING THE RENUMBERING. THIS OPTION IS USEFUL TO C PREPARE A RENUMBERING FOR USE BY THE GRAPHICS PROGRAMS. C IF "REDEEM" IS FALSE, FAKE NODES ARE NOT RENUMBERED, AND C WILL APPEAR AT THE END OF THE LIST AS BEFORE. C "REDEEM" ALSO CONTROLS THE TYPE OF OUTPUT FILE; IF TRUE, C THEN ONLY A LIST OF OLD NODE NUMBERS, NEW NODE NUMBERS, AND C A CROSS-REFERENCE (OLD NUMBERS LISTED IN ORDER OF NEW) IS C OUTPUT ON UNIT "UNITO". IF "REDEEM" IS FALSE, HOWEVER, A C COMPLETE, RESORTED FINITE ELEMENT GRID FILE IS WRITTEN. C THIS REFLECTS THE PRACTICAL REALITY THAT ONE USUALLY C CHOOSES "REDEEM=.TRUE." WHEN PREPARING TO DO GRAPHICS, BUT C CHOOSES "REDEEM=.FALSE." WHEN OPTIMIZING A HAND-EDITED GRID C PRIOR TO FINITE ELEMENT CALCULATIONS. C PARAMETER (MAXNOD=12000,MAXEL=4000,MAXFEL=2000,MAXGSI=4000000) C LOGICAL*1 CHECKE,CHECKF,CHECKN,DONEXT,EDGEFS,EDGETS, + INSIDE,LGRAPH LOGICAL BRIEF,REDEEM INTEGER*2 NODTYP CHARACTER*80 TITLE1 C DIMENSION ATNODE(MAXNOD), + CHECKE(MAXEL),CHECKF(MAXFEL),CHECKN(MAXNOD), + DONEXT(MAXNOD),DQDTDA(MAXNOD),ELEV(MAXNOD), + EDGEFS(2,MAXFEL),EDGETS(3,MAXEL), + FAZ(2,MAXFEL),FDIP(3,MAXFEL), + IALIAS(MAXNOD),IDIAG(MAXNOD),INSIDE(MAXNOD), + IUSER(MAXNOD),LGRAPH(MAXGSI),LIST(MAXNOD), + LISTTK(MAXNOD),MAXLST(MAXNOD), + NODEF(6,MAXFEL),NODES(6,MAXEL),NODTYP(MAXNOD), + OFFSET(MAXFEL),TLNODE(MAXNOD), + XNODE(MAXNOD),YNODE(MAXNOD),ZMNODE(MAXNOD) C DATA IUNITI /1/ DATA IUNITO /2/ DATA IUNITT /6/ C MXNODE=MAXNOD MXEL=MAXEL MXFEL=MAXFEL MXGSIZ=MAXGSI C WRITE (IUNITT,1) 1 FORMAT (//' Program NUMBER' + /' (version for grids of 6-node continuum triangles,' + /' and, optionally, 6-node fault elements, such as' + /' those produced by interactive program DRAWGRID)' + //' Its function is to renumber the nodes for lower' + /' bandwidth, allowing more economical solution of' + /' any kind of linear system. Always apply NUMBER' + /' after editing with DRAWGRID, and before plotting' + /' with FAULTS2AI or PLATES2AI, and before solving' + /' neotectonic problems with FAULTS or PLATES.' + /' Peter Bird, UCLA, 4/97') WRITE (IUNITT,2) 2 FORMAT (/' NUMBER WORKS IN TWO MODES:' + /' 1: gives a second renumbered .FEG file as output.' + /' 2: writes only a table of new node numbers (aliases)' + /' (needed for some types of graphical plots).' + /' Note that in mode 1 the fake boundary nodes are kept' + /' distinct, but in mode 2 they are mixed in.' + /' ENTER YOUR CHOICE OF MODE: '\) 3 READ (*,*) IMODE IF (IMODE.EQ.1) THEN REDEEM=.FALSE. ELSE IF (IMODE.EQ.2) THEN REDEEM=.TRUE. ELSE WRITE (IUNITT,4) 4 FORMAT (' ERROR: ENTER ONLY 1 OR 2: '\) GO TO 3 END IF C CALL GETNET (INPUT,IUNITI,IUNITT, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,N1000,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE, + WORK,CHECKE,CHECKF,CHECKN) C BRIEF=.TRUE. C CALL SQUARE (INPUT,BRIEF,FAZ,FDIP,IUNITT, + MXEL,MXFEL,MXNODE, + NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000, + MODIFY,XNODE,YNODE, + OUTPUT,LENGTH,LIST,NODTYP, + WORK,CHECKN,EDGEFS,EDGETS,LISTTK) C IF (REDEEM) THEN CALL NUMBRE (INPUT,IUNITT,LENGTH,LIST, + MXEL,MXFEL,MXGSIZ,MXNODE, + NFL,NODEF,NODES,NODTYP,NUMNOD, + NUMEL,NUMNOD,XNODE,YNODE, + OUTPUT,IALIAS,IUSER,MAXDIF,MINDIF, + WORK,DONEXT,IDIAG,INSIDE,LGRAPH, + MAXLST) ELSE CALL NUMBRE (INPUT,IUNITT,LENGTH,LIST, + MXEL,MXFEL,MXGSIZ,MXNODE, + NFL,NODEF,NODES,NODTYP,NREALN, + NUMEL,NUMNOD,XNODE,YNODE, + OUTPUT,IALIAS,IUSER,MAXDIF,MINDIF, + WORK,DONEXT,IDIAG,INSIDE,LGRAPH, + MAXLST) ENDIF C C REPORT RESULTS TO USER ON DEVICE "IUNITT": C WRITE (IUNITT,15) MAXDIF 15 FORMAT (/ /' BANDWIDTH BEFORE RENUMBERING = ',I5) WRITE (IUNITT,18) MINDIF 18 FORMAT (/ /' BANDWIDTH AFTER RENUMBERING = ',I5) C C WRITE (IUNITT,20) 20 FORMAT (/' READY TO WRITE OUTPUT FILE:') IF (IMODE.EQ.1) WRITE (IUNITT,21) 21 FORMAT (' NOTE: A filename extension of .FEG is suggested!') IF (REDEEM) THEN DO 23 I=1,NUMNOD WRITE (IUNITO,22) I, IALIAS(I), IUSER(I) 22 FORMAT (3I10) 23 CONTINUE ELSE IF (MINDIF.LT.MAXDIF) THEN C C REORDER THE VALUES IN THE SIMPLE FLOATING-POINT VECTOR ARRAYS: C CALL MIXUPX (INPUT,IALIAS,MXNODE,NUMNOD, + MODIFY,XNODE, + WORK,ATNODE) CALL MIXUPX (INPUT,IALIAS,MXNODE,NUMNOD, + MODIFY,YNODE, + WORK,ATNODE) CALL MIXUPX (INPUT,IALIAS,MXNODE,NUMNOD, + MODIFY,ELEV, + WORK,ATNODE) CALL MIXUPX (INPUT,IALIAS,MXNODE,NUMNOD, + MODIFY,DQDTDA, + WORK,ATNODE) CALL MIXUPX (INPUT,IALIAS,MXNODE,NUMNOD, + MODIFY,ZMNODE, + WORK,ATNODE) CALL MIXUPX (INPUT,IALIAS,MXNODE,NUMNOD, + MODIFY,TLNODE, + WORK,ATNODE) C C MODIFY THE VALUES IN THE 6-ENTRY INTEGER ARRAYS: C CALL MIXUPI (INPUT,IALIAS,MXEL,MXNODE,NUMEL, + MODIFY,NODES) CALL MIXUPI (INPUT,IALIAS,MXFEL,MXNODE,NFL, + MODIFY,NODEF) C C OUTPUT MODIFIED FILE ON DEVICE "IUNITO": C CALL PUTNET (INPUT,IUNITO, + BRIEF,DQDTDA,ELEV,FAZ,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE) ELSE WRITE (IUNITT,99) 99 FORMAT (/ /' SINCE BANDWIDTH COULD NOT BE IMPROVED,', + ' GRID WAS NOT MODIFIED, AND NO OUTPUT FILE', + ' WAS WRITTEN.') ENDIF ENDIF C STOP ' ' END PROGRAM NUMBER C C C SUBROUTINE NUMBRE (INPUT,IUNITT,LENGTH,LIST, + MXEL,MXFEL,MXGSIZ,MXNODE, + NFL,NODEF,NODES,NODTYP,NREALN, + NUMEL,NUMNOD,XNODE,YNODE, + OUTPUT,IALIAS,IUSER,MAXDIF,MINDIF, + WORK,DONEXT,IDIAG,INSIDE,LGRAPH, + MAXLST) C C FIND A GENERAL RENUMBERING OF THE NODES OF THE FINITE-ELEMENT C GRID WHICH IS (NEARLY) OPTIMAL IN TERMS OF MINIMIZING THE C BANDWIDTH OF THE CONNECTIVITY MATRIX (A SURROGATE C FOR THE STIFFNESS MATRIX, BUT SYMMETRICAL, AND WITH ONLY ONE C ROW AND ONE COLUMN PER NODE). C C IN THE SECTION WHICH FILLS THE CONNECTIVITY MATRIX C "LGRAPH", IT IS ASSUMED THAT THE GRID IS MADE UP OF C A SET OF "NUMEL" 6-NODE TRIANGULAR ELEMENTS (WHOSE OLD C NODE NUMBERS ARE STORED AS "NODES(K=1..6,I=1..NUMEL)"), AND C "NFL" 6-NODE FAULT ELEMENTS (WHOSE OLD NODE NUMBERS ARE C STORED AS "NODEF(K=1..6,I=1..NFL)" ). C NFL MAY BE ZERO. C THIS SECTION CAN BE EASILY MODIFIED FOR OTHER TYPES OF ELEMENTS. C IT DOES NOT MATTER WHETHER OR NOT THE 6 NUMBERS DEFINING C EACH ELEMENT OCCUR IN ANY PARTICULAR ORDER, OR NOT. C HOWEVER, ARRAY "NODTYP" MUST CONTAIN 1 IF THE NODE IS A CORNER C OR END NODE, BUT 2 IF IT IS A MIDPOINT NODE. C DEGENERATE ELEMENTS WITH TWO OR MORE IDENTICAL NODE NUMBERS C ARE ALSO PERMITTED. C C THE TOTAL NUMBER OF NODES IN THE GRID IS INPUT PARAMETER C "NUMNOD". THE FIRST "NREALN" OF THESE WILL BE RENUMBERED. C USUALLY, NUMNOD WILL EQUAL NREALN. OBVIOUSLY, NREALN CANNOT C BE GREATER THAN NUMNOD. THE DISTINCTION IS MADE IN CASE C YOU ARE USING A FINITE ELEMENT PROGRAM LIKE "FAULTS" IN C WHICH FAKE NODES (THOSE WITH NO DEGREES OF FREEDOM) ARE C LISTED LAST (WITH NUMBERS FROM NREALN+1 TO NUMNOD), AND C YOU WANT TO KEEP THIS GROUP AT THE END, AND NOT RENUMBER C THEM. C C THE ALGORITHM USED IS AN ORIGINAL VARIANT OF THE C CUTHILL-MCKEE ALGORITHM C DESCRIBED BY H. R. SCHWARZ IN "FINITE ELEMENT METHODS", C "COMPUTATIONAL MATHEMATICS AND APPLICATIONS" SERIES, C ACADEMIC PRESS/HARCOURT-BRACE-JOVANOVITCH, 1988, PP.167-185. C C INSTEAD OF USING DEGREE TO DECIDE WHICH NODES WILL BE NUMBERED C FIRST, THIS VARIANT NUMBERS NODES WITHIN ONE LEVEL ON THE BASIS C OF THEIR AZIMUTH FROM THE INITIAL NODE. C C BECAUSE OF THIS CHOICE OF ALGORITHM, THIS CODE CAN NOT BE EXPECTED C TO WORK WELL FOR GRIDS WHICH ARE NOT SIMPLY CONNECTED (I.E., C THOSE WITH HOLES IN THEM), AND IT CANNOT BE EXPECTED TO WORK C WELL FOR GRIDS WHICH HAVE EXTREMELY CONCAVE BOUNDARIES. C (FORTUNATELY, BOTH OF THESE CASES DO NOT ARISE IN EARTH SCIENCES, C WHERE THE SIDE BOUNDARIES OF THE GRID ARE ALWAYS ARTIFICIAL, AND C ARE ALMOST ALWAYS CHOSEN AS A CONVEX POLYHEDRON.) C C THE PRODUCTS OF THIS SUBPROGRAM ARE: C *"IALIAS" = AN INTEGER VECTOR OF NEW NODE NUMBERS, ARRANGED C IN THE ORDER OF THE OLD NODE NUMBERS. C *"IUSER" = AN INTEGER VECTOR OF OLD NODE NUMBERS, ARRANGED C IN THE ORDER OF THE NEW NODE NUMBERS. C ("IALIAS" AND "IUSER" ARE MUTUAL INVERSES, AS SHOWN BY C THE FOLLOWING IDENTITIES: C IUSER(IALIAS(I))=I C IALIAS(IUSER(J))=J ) C * MAXDIF = HALF BANDWIDTH BEFORE RENUMBERING. C * MINDIF = HALF BANDWIDTH AFTER RENUMBERING. C (THE DEFINITION OF HALF-BANDWIDTH USED HERE C DOES NOT INCLUDE THE DIAGONAL, AND CAN BE DEFINED AS THE C MAXIMUM OF THE ABSOLUTE VALUE OF THE DIFFERENCE OF THE C (NEW) NUMBERS IN ALL POSSIBLE PAIRS OF CONNECTED NODES.) C C NOTE: THE FOLLOWING TYPE DECLARATION IS VERY IMPORTANT FOR C ECONOMIZING STORAGE, AS LGRAPH WILL HAVE A SIZE OF C APPROXIMATELY NREALN*(NREALN-2)/2 LOGICAL ENTRIES. LOGICAL*1 DONEXT,INSIDE,LGRAPH C INTEGER*2 NODTYP C LOGICAL FOUND,LESSOK,MOREOK C DIMENSION DONEXT(MXNODE),IALIAS(MXNODE),IDIAG(MXNODE), + INSIDE(MXNODE), + IUSER(MXNODE),LGRAPH(MXGSIZ),LIST(MXNODE), + MAXLST(MXNODE), + NODEF(6,MXFEL),NODES(6,MXEL), + NODTYP(MXNODE), + XNODE(MXNODE),YNODE(MXNODE) C C STATEMENT FUNCTION, GIVING STORAGE LOCATION (VECTOR INDEX) OF THE C ELEMENT OF THE GRAPH MATRIX AT ROW "IROW" AND COLUMN "JCOL": INDEXN(IROW,JCOL)=IDIAG(MIN(IROW,JCOL))+ABS(JCOL-IROW) C IF (NREALN.GT.NUMNOD) THEN WRITE (IUNITT,5) NREALN,NUMNOD 5 FORMAT (/ /' SUBPROGRAM NUMBER WAS ASKED TO RENUMBER', + ' THE FIRST ',I6,' NODES,'/ + ' OUT OF A TOTAL OF ',I6,'.'/ + ' THIS IS IMPOSSIBLE.'/ + ' INPUT PARAMETER NREALN MUST NOT EXCEED', + ' INPUT PARAMETER NUMNOD.') STOP ENDIF C C COMPUTE THE LOCATIONS OF THE DIAGONAL ELEMENTS OF C THE CONNECTIVITY MATRIX "LGRAPH". C (NOTE THAT THE GRAPH IS SYMMETRIC, SO WE ONLY STORE THE C UPPER TRIANGLE, AND WE ACCESS THE MATRIX C INDIRECTLY USING STATEMENT FUNCTION "INDEXN".) C NUSED=0 DO 10 I=1,NREALN IDIAG(I)=NUSED NINROW=NREALN-I NUSED=NUSED+NINROW 10 CONTINUE IF (NUSED.GT.MXGSIZ) THEN WRITE (IUNITT,11) MXGSIZ, NUSED 11 FORMAT (/ /' INCREASE PARAMETER MAXGSI ', + '(CURRENTLY ',I12,') TO AT LEAST ',I12/ + ' AND RECOMPILE THIS PROGRAM.') STOP ENDIF C C CREATE THE GRAPH MATRIX C DO 20 I=1,MXGSIZ LGRAPH(I)=.FALSE. 20 CONTINUE C DO 50 I=1,NUMEL DO 40 J=1,5 DO 30 K=J+1,6 IROW=NODES(J,I) JCOL=NODES(K,I) IF ((IROW.LT.1).OR.(IROW.GT.NUMNOD)) THEN WRITE (IUNITT,21) IROW,I 21 FORMAT (/ /' ERR0R IN INPUT DATA:'/ + ' NODE NUMBER ',I5,' IN DEFINITION OF', + ' TRIANGULAR ELEMENT ',I5,' IS ILLEGAL.') STOP ENDIF IF ((JCOL.LT.1).OR.(JCOL.GT.NUMNOD)) THEN WRITE (IUNITT,21) JCOL,I STOP ENDIF IF ((IROW.LE.NREALN).AND.(JCOL.LE.NREALN) + .AND.(IROW.NE.JCOL)) + LGRAPH(INDEXN(IROW,JCOL))=.TRUE. 30 CONTINUE 40 CONTINUE 50 CONTINUE C DO 150 I=1,NFL DO 140 J=1,5 DO 130 K=J+1,6 IROW=NODEF(J,I) JCOL=NODEF(K,I) IF ((IROW.LT.1).OR.(IROW.GT.NUMNOD)) THEN WRITE (IUNITT,51) IROW,I 51 FORMAT (/ /' ERR0R IN INPUT DATA:'/ + ' NODE NUMBER ',I5,' IN DEFINITION OF', + ' FAULT ELEMENT ',I5,' IS ILLEGAL.') STOP ENDIF IF ((JCOL.LT.1).OR.(JCOL.GT.NUMNOD)) THEN WRITE (IUNITT,51) JCOL,I STOP ENDIF IF ((IROW.LE.NREALN).AND.(JCOL.LE.NREALN) + .AND.(IROW.NE.JCOL)) + LGRAPH(INDEXN(IROW,JCOL))=.TRUE. 130 CONTINUE 140 CONTINUE 150 CONTINUE C C DETERMINE THE ORIGINAL BANDWIDTH C MAXDIF=0 DO 152 I=1,NREALN-1 DO 151 J=I+1,NREALN IF (LGRAPH(INDEXN(I,J))) MAXDIF=MAX(MAXDIF,ABS(J-I)) 151 CONTINUE 152 CONTINUE C M1=1 M9=LENGTH C C=================================================================== C 200 DO 10000 M=M1,M9 C C CREATE A RENUMBERING SCHEME, USING INITIAL NODE "LIST(M)" C C (NECESSARY INITIALIZATION; RESET NUMBERING SCHEME TO NULL) DO 201 I=1,NREALN IALIAS(I)=0 IUSER(I)=0 201 CONTINUE C C (BEGIN REAL WORK) C X1=XNODE(LIST(M)) Y1=YNODE(LIST(M)) IF (M.LT.LENGTH) THEN MMORE=M+1 ELSE MMORE=1 ENDIF IF (M.GT.1) THEN MLESS=M-1 ELSE MLESS=LENGTH ENDIF XMORE=XNODE(LIST(MMORE)) YMORE=YNODE(LIST(MMORE)) XLESS=XNODE(LIST(MLESS)) YLESS=YNODE(LIST(MLESS)) DXMORE=XMORE-X1 DYMORE=YMORE-Y1 DXLESS=XLESS-X1 DYLESS=YLESS-Y1 MOREOK=(ABS(DXMORE)+ABS(DYMORE)).GT.0. LESSOK=(ABS(DXLESS)+ABS(DYLESS)).GT.0. IF (MOREOK.AND.LESSOK) THEN ARGMOR=ATAN2F(DYMORE,DXMORE) ARGLES=ATAN2F(DYLESS,DXLESS) ELSE IF (MOREOK) THEN ARGMOR=ATAN2F(DYMORE,DXMORE) ARGLES=ARGMOR-3.14159 IF (ARGLES.LT.-3.14159) ARGLES=ARGLES+6.28318 ELSE IF (LESSOK) THEN ARGLES=ATAN2F(DYLESS,DXLESS) ARGMOR=ARGLES+3.14159 IF (ARGMOR.GT.3.14159) ARGMOR=ARGMOR-6.28318 ELSE WRITE (IUNITT,202) LIST(M) 202 FORMAT (/ /' PROBLEM IN SUBPROGRAM NUMBER:'/ + ' BOUNDARY AZIMUTH IS UNDEFINED AT NODE ',I5/ + ' BECAUSE BOTH NEIGHBORING BOUNDARY NODES ARE LOCATED'/ + ' AT EXACTLY THE SAME POINT.') STOP ENDIF ARGCUT=0.5*(ARGMOR+ARGLES) XM=COS(ARGMOR) YM=SIN(ARGMOR) XC=COS(ARGCUT) YC=SIN(ARGCUT) CROSS=XC*YM-YC*XM IF (CROSS.LT.0.) THEN ARGCUT=ARGCUT+3.14159 IF (ARGCUT.GT.3.14159) ARGCUT=ARGCUT-6.28318 ENDIF CUTDEG=ARGCUT*57.298 C C NUMBER THE FIRST NODE C IF (LIST(M).LE.NREALN) THEN ISTART=LIST(M) ELSE ISTART=1 FOUND=.FALSE. DO 205 I=1,NUMEL DO 204 K=1,6 IF (NODES(K,I).EQ.LIST(M)) THEN DO 203 J=1,6 IF (NODES(J,I).LE.NREALN) THEN FOUND=.TRUE. ISTART=NODES(J,I) ENDIF 203 CONTINUE ENDIF 204 CONTINUE 205 CONTINUE IF (.NOT.FOUND) THEN DO 215 I=1,NFL DO 214 K=1,6 IF (NODEF(K,I).EQ.LIST(M)) THEN DO 213 J=1,6 IF (NODEF(J,I).LE.NREALN) THEN FOUND=.TRUE. ISTART=NODEF(J,I) ENDIF 213 CONTINUE ENDIF 214 CONTINUE 215 CONTINUE ENDIF ENDIF IALIAS(ISTART)=1 IUSER(1)=ISTART NDONE=1 LEVEL=0 LAST1=1 LAST9=1 C C----------------- BEGIN INDEFINATE LOOP ON LEVEL -------------------- C C MARK THE NODES TO BE NUMBERED IN THE NEXT LEVEL (UNNUMBERED NODES C CONNECTED TO NODES THAT HAVE ALREADY BEEN NUMBERED) C 220 LEVEL=LEVEL+1 DO 250 I=1,NREALN DONEXT(I)=.FALSE. 250 CONTINUE DO 300 I=LAST1,LAST9 NODE=IUSER(I) IF (NODE.GT.1) THEN DO 260 J=1,NODE-1 DONEXT(J)=DONEXT(J).OR. + ( (IALIAS(J).EQ.0) .AND. LGRAPH(INDEXN(J,NODE)) ) 260 CONTINUE ENDIF IF (NODE.LT.NREALN) THEN DO 270 J=NODE+1,NREALN DONEXT(J)=DONEXT(J).OR. + ( (IALIAS(J).EQ.0) .AND. LGRAPH(INDEXN(J,NODE)) ) 270 CONTINUE ENDIF 300 CONTINUE C C HOW MANY ARE TO BE NUMBERED AT THIS LEVEL? C NTODO=0 DO 310 I=1,NREALN IF (DONEXT(I)) NTODO=NTODO+1 310 CONTINUE IF (NTODO.EQ.0) THEN IF (NDONE.EQ.NREALN) THEN GO TO 8000 ELSE WRITE (IUNITT,311) ISTART,NDONE,IUSER(NDONE),NREALN 311 FORMAT (/ /' ERR0R IN GRID TOPOLOGY:'/ + ' AFTER BEGINNING NUMBERING AT NODE ',I5/ + ' AND NUMBERING ',I5,' NODES UP THROUGH ', + 'NODE ',I5,','/' NO CONNECTED NEIGHBORING', + ' NODES COULD BE FOUND,'/' SO IT IS NOT', + ' POSSIBLE TO NUMBER ',I5,' NODES.') STOP ENDIF ENDIF C C FIRST, DIVIDE NODES TO BE NUMBERED AT THIS LEVEL INTO "INSIDE" C AND OUTSIDE NODES. (INSIDE NODES WILL GET LOWER NUMBERS.) C DO 313 I=1,NREALN INSIDE(I)=DONEXT(I).AND.(NODTYP(I).EQ.2) IF (INSIDE(I)) THEN DO 312 J=1,NREALN IF (J.NE.I) THEN IF (LGRAPH(INDEXN(I,J))) THEN INSIDE(I)=INSIDE(I).AND. + ((IALIAS(J).GT.0).OR.DONEXT(J)) IF (.NOT.INSIDE(I)) GO TO 313 ENDIF ENDIF 312 CONTINUE ENDIF 313 CONTINUE NINSID=0 DO 314 I=1,NREALN IF (INSIDE(I)) NINSID=NINSID+1 314 CONTINUE C C LOOP ON INSIDE NODES, TAKING THEM IN ORDER OF THEIR ARGUMENT C ANGLES FROM THE INITIAL NODE. C DO 400 N=1,NINSID C C FIND LOWEST-ARGUMENT NODE LEFT TO BE DONE C ARGMIN=999. DO 350 I=1,NREALN IF (INSIDE(I)) THEN ARG=ATAN2F(YNODE(I)-Y1,XNODE(I)-X1) IF (ARG.LT.ARGCUT) ARG=ARG+6.28318 IF (ARG.LT.ARGMIN) THEN ARGMIN=ARG NODE=I ENDIF ENDIF 350 CONTINUE C C NUMBER IT ! ! ! C NDONE=NDONE+1 IUSER(NDONE)=NODE IALIAS(NODE)=NDONE INSIDE(NODE)=.FALSE. DONEXT(NODE)=.FALSE. 400 CONTINUE C C LOOP ON OUTSIDE NODES, TAKING THEM IN ORDER OF THEIR ARGUMENT C ANGLES FROM THE INITIAL NODE. C NOUT=0 DO 401 I=1,NREALN IF (DONEXT(I)) NOUT=NOUT+1 401 CONTINUE DO 500 N=1,NOUT C C FIND LOWEST-ARGUMENT NODE LEFT TO BE DONE C ARGMIN=999. DO 450 I=1,NREALN IF (DONEXT(I)) THEN ARG=ATAN2F(YNODE(I)-Y1,XNODE(I)-X1) IF (ARG.LT.ARGCUT) ARG=ARG+6.28318 IF (ARG.LT.ARGMIN) THEN ARGMIN=ARG NODE=I ENDIF ENDIF 450 CONTINUE C C NUMBER IT ! ! ! C NDONE=NDONE+1 IUSER(NDONE)=NODE IALIAS(NODE)=NDONE DONEXT(NODE)=.FALSE. 500 CONTINUE C C SUMMARIZE THE LAST LEVEL DONE C LAST1=LAST9+1 LAST9=LAST1+NTODO-1 C C LOOP TO NEXT LEVEL C IF (NDONE.LT.NREALN) GO TO 220 C C------------------- END INDEFINATE LOOP ON LEVEL ----------------- C C DETERMINE THE NEW HALF-BANDWIDTH C 8000 MINDIF=0 DO 9020 I=1,NREALN-1 DO 9010 J=I+1,NREALN IF (LGRAPH(INDEXN(I,J))) THEN IALDIF=ABS(IALIAS(I)-IALIAS(J)) MINDIF=MAX(MINDIF,IALDIF) IF ((M1.EQ.M9).AND.(IALDIF.EQ.MINMAX)) THEN WRITE (IUNITT,9005) I,J,IALIAS(I),IALIAS(J), + MINMAX 9005 FORMAT (' OLD NODES ',I5,' AND ',I5,' HAVE', + ' NEW NUMBERS ',I5,' AND ',I5/ + ' WHOSE DIFFERENCE OF ',I5,' DEFINES', + ' THE NEW BANDWIDTH.') ENDIF ENDIF 9010 CONTINUE 9020 CONTINUE MAXLST(M)=MINDIF IF (ISTART.EQ.LIST(M)) THEN WRITE (IUNITT,9031) M, LIST(M),MINDIF ELSE WRITE (IUNITT,9032) M, LIST(M),ISTART,MINDIF ENDIF 9031 FORMAT(' ATTEMPT NUMBER ',I5,' BEGINS WITH BOUNDARY NODE',I5,'.' + /' WHEN THIS NODE IS #1, THE BANDWIDTH IS ',I5,'.') 9032 FORMAT(' ATTEMPT NUMBER ',I5,' BEGINS WITH BOUNDARY NODE',I5,',' + /' AND NEAREST REAL NODE TO THAT IS NUMBER ',I5,'.'/ + ' WHEN THIS NODE IS #1, THE BANDWIDTH IS ',I5,'.') C C 10000 CONTINUE C================================================================= C C REPEAT THE NUMBERING FROM THE BEST STARTING POINT C IF (M9.GT.M1) THEN MINMAX=NREALN DO 10001 M=1,LENGTH IF (MAXLST(M).LT.MINMAX) THEN MINMAX=MAXLST(M) MSTART=M ENDIF 10001 CONTINUE M1=MSTART M9=MSTART GO TO 200 ENDIF C IF (NUMNOD.GT.NREALN) THEN DO 90000 I=NREALN+1,NUMNOD IALIAS(I)=I IUSER(I)=I 90000 CONTINUE ENDIF C RETURN END SUBROUTINE NUMBRE C C C SUBROUTINE GETNET (INPUT,IUNIT7,IUNIT8, + MXEL,MXFEL,MXNODE, + OUTPUT,BRIEF,DQDTDA,ELEV,FAZ,FDIP, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,N1000, + OFFSET,TITLE1,TLNODE,XNODE,YNODE, + ZMNODE, + WORK,CHECKE,CHECKF,CHECKN) C C READ FINITE ELEMENT GRID FROM UNIT "IUNIT7". C ECHO THE IMPORTANT VALUES TO A PRINT DATASET ON UNIT "IUNIT8". C C READ NUMBER OF NODES, AND HOW MANY ARE "REAL" VERSUS "FAKE". C INPUT NODAL LOCATIONS (X,Y), ELEVATIONS, HEAT-FLOW, AND ISOSTATIC C GRAVITY ANOMALIES (ONE RECORD PER NODE). C CHARACTER*80 TITLE1 LOGICAL ALLOK,BRIEF C C NOTE: FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL*1 CHECKE,CHECKF,CHECKN C DIMENSION CHECKE(MXEL),CHECKF(MXFEL),CHECKN(MXNODE), + DQDTDA(MXNODE),ELEV(MXNODE), + FAZ(2,MXFEL),FDIP(3,MXFEL), + NODEF(6,MXFEL),NODES(6,MXEL), + OFFSET(MXFEL),TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(3),IFN(6),VECTOR(7) C WRITE (IUNIT8,1) IUNIT7 1 FORMAT (/' ATTEMPTING TO READ FINITE ELEMENT GRID FROM UNIT',I3/) READ (IUNIT7,2,IOSTAT=IOS) TITLE1 2 FORMAT (A80) WRITE (IUNIT8,3) TITLE1 3 FORMAT(/' TITLE OF FINITE ELEMENT GRID ='/' ',A80) C C READ NODES, AND SUBDIVIDE INTO 2 CLASSES C (OPTION "BRIEF" SUPPRESSES MOST OUTPUT) C READ (IUNIT7,*) NUMNOD,NREALN,NFAKEN,N1000,BRIEF C IF (NUMNOD.NE.(NREALN+NFAKEN)) THEN WRITE (IUNIT8,5) 5 FORMAT (/' INCONSISTENT DATA:'/ + ' NUMBER OF NODES SHOULD EQUAL TOTAL OF REAL', + ' NODES AND FAKE NODES.') STOP ENDIF IF (NUMNOD.GT.MXNODE) THEN WRITE (IUNIT8,10) NUMNOD 10 FORMAT(/' INCREASE PARAMETER MAXNOD TO BE AT LEAST' + /' THE NUMBER OF NODES (',I6,') AND RECOMPILE.') STOP ENDIF N2=2*NUMNOD IF (NREALN.GT.N1000) THEN WRITE (IUNIT8,20) NREALN 20 FORMAT (/' INCREASE THE DATA VALUE N1000 TO BE GREATER' + /' OR EQUAL TO NREALN (',I6,') AND RECOMPILE.') STOP ENDIF C NBASE=N1000+1 NTOP=N1000+NFAKEN IF (BRIEF) THEN WRITE (IUNIT8,35) 35 FORMAT(/' (SINCE OPTION BRIEF=.TRUE., GRID WILL NOT BE ', + 'ECHOED HERE. BE CAREFUL!!!)') ELSE WRITE (IUNIT8,40) NUMNOD,NREALN,NREALN 40 FORMAT (/' THERE ARE',I5,' NODES IN THE GRID:'/ ' ',I5, + ' OF THESE ARE NUMBERED 1-',I4,' AND THESE REAL NODES', + ' HAVE TWO VELOCITY VARIABLES UNLESS CONSTRAINED.') IF (NFAKEN.GT.0) WRITE (IUNIT8,42) NFAKEN,NBASE,NTOP 42 FORMAT(' ',I5, + ' OF THESE ARE NUMBERED ',I6,'-',I6,' AND THESE ARE', + ' ARTIFICIAL;'/' THEIR VELOCITIES MUST BE', + ' COMPLETELY SPECIFIED.') WRITE (IUNIT8,49) 49 FORMAT(/ + ' (NOTE: X AND Y COORDINATES MAY BE ZERO FOR NODES WHICH' + ,' WILL BE AT MIDPOINTS OF ELEMENT SIDES AND/OR FAULTS.'/ + ' THE PROGRAM WILL INTERPOLATE VALUES FOR THESE.)') WRITE (IUNIT8,50) 50 FORMAT (/' ', + ' CRUSTAL MANTLE'/ + ' NODE X Y ELEVATION', + ' HEAT-FLOW THICKNESS THICKNESS'/) ENDIF DO 90 K=1,NUMNOD CHECKN(K)=.FALSE. 90 CONTINUE DO 100 K=1,NUMNOD CALL READN (INPUT,IUNIT7,IUNIT8,7, + OUTPUT,VECTOR) INDEX=VECTOR(1)+0.5 XI=VECTOR(2) YI=VECTOR(3) ELEVI=VECTOR(4) QI=VECTOR(5) ZMI=VECTOR(6) TLI=VECTOR(7) C (NODES NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) IF (INDEX.LE.NREALN) THEN I=INDEX ELSE I=INDEX-N1000+NREALN ENDIF CHECKN(I)=.TRUE. ELEV(I)=ELEVI DQDTDA(I)=QI IF (QI.LT.0.) THEN WRITE (IUNIT8,91) 91 FORMAT (' NEGATIVE HEAT-FLOW IS NON-PHYSICAL.') STOP ENDIF XNODE(I)=XI YNODE(I)=YI IF (ZMI.LT.0.) THEN WRITE (IUNIT8,92) 92 FORMAT(' NEGATIVE CRUSTAL THICKNESS IS NON-PHYSICAL.') STOP ENDIF ZMNODE(I)=ZMI IF (TLI.LT.0.) THEN WRITE (IUNIT8,93) 93 FORMAT(' NEGATIVE MANTLE LITHOSPHERE THICKNESS', + ' IS NON-PHYSICAL.') STOP ENDIF TLNODE(I)=TLI IF (.NOT.BRIEF) THEN WRITE (IUNIT8,95) INDEX,XI,YI,ELEVI,QI,ZMI,TLI 95 FORMAT (' ',I5,1P,2E11.3,4E10.2) ENDIF 100 CONTINUE ALLOK=.TRUE. DO 101 I=1,NUMNOD ALLOK=ALLOK.AND.CHECKN(I) 101 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,102) 102 FORMAT(' THE FOLLOWING NODES WERE NEVER READ:') DO 104 I=1,NUMNOD IF (.NOT.CHECKN(I)) WRITE(IUNIT8,103)I 103 FORMAT (' ',36X,I6) 104 CONTINUE STOP ENDIF C C READ TRIANGULAR ELEMENTS C READ (IUNIT7,*) NUMEL IF (NUMEL.GT.MXEL) THEN WRITE (IUNIT8,108) NUMEL 108 FORMAT(/' INCREASE PARAMETER MAXEL TO BE AT LEAST EQUAL' + /' TO THE NUMBER OF ELEMENTS (',I6,') AND RECOMPILE.') STOP ENDIF DO 109 K=1,NUMEL CHECKE(K)=.FALSE. 109 CONTINUE IF (.NOT.BRIEF) WRITE (IUNIT8,110) NUMEL 110 FORMAT(/' THERE ARE ',I6,' TRIANGULAR CONTINUUM ELEMENTS.'/ + ' (NODE NUMBERS FOR EACH ARE GIVEN CORNERS-FIRST, COUNTER', + 'CLOCKWISE; THEN'/' MIDPOINTS, COUNTERCLOCKWISE, BEGINNING' + ,' WITH THE MIDPOINT BETWEEN CORNER #1 AND CORNER #2)'/ / + ' ELEMENT C1 C2 C3 M1 M2', + ' M3') DO 200 K=1,NUMEL C (ELEMENTS NEED NOT BE INPUT IN ORDER, BUT MUST ALL BE PRESENT.) READ (IUNIT7,*) I,(IFN(J),J=1,6) CHECKE(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,120) I,(IFN(J),J=1,6) 120 FORMAT (' ',I6,':',6I6) DO 130 J=1,6 N=IFN(J) IF ((N.LE.0).OR.(N.GT.NTOP).OR. + ((N.GT.NREALN).AND.(N.LE.N1000))) THEN WRITE (IUNIT8,125) N 125 FORMAT (' NODE NUMBER ',I6,' IS ILLEGAL.') STOP ENDIF IF (N.GT.NREALN) N=N-N1000+NREALN NODES(J,I)=N 130 CONTINUE 200 CONTINUE ALLOK=.TRUE. DO 201 I=1,NUMEL ALLOK=ALLOK.AND.CHECKE(I) 201 CONTINUE IF (.NOT.ALLOK) THEN WRITE (IUNIT8,202) 202 FORMAT (' THE FOLLOWING ELEMENTS WERE NEVER READ:') DO 204 I=1,NUMEL IF (.NOT.CHECKE(I)) WRITE(IUNIT8,203)I 203 FORMAT (' ',39X,I6) 204 CONTINUE STOP ENDIF C C READ FAULT ELEMENTS C READ(IUNIT7,*) NFL IF (NFL.GT.MXFEL) THEN WRITE (IUNIT8,220)NFL 220 FORMAT (/' INCREASE PARAMETER MAXFEL TO BE AT LEAST EQUAL' + /' TO THE NUMBER OF FAULTS (',I6,') AND RECOMPILE.') STOP ENDIF DO 222 I=1,NFL CHECKF(I)=.FALSE. 222 CONTINUE IF (.NOT.BRIEF) WRITE(IUNIT8,230) NFL 230 FORMAT(/ /' THERE ARE ',I6,' CURVILINEAR FAULT ELEMENTS.') IF ((.NOT.BRIEF).AND.(NFL.GT.0)) WRITE(IUNIT8,231) 231 FORMAT (/' (THE 6 NODE NUMBERS DEFINING EACH ELEMENT MUST BE', + ' IN A COUNTERCLOCKWISE ORDER:'/ + ' N1, N2, AND N3 ARE IN LEFT-TO-RIGHT SEQUENCE ON THE', + ' NEAR SIDE,'/ + ' THEN N4 IS OPPOSITE N3, N5 IS OPPOSITE N2, AND ', + 'N6 IS OPPOSITE N1.)'/' (FAULT DIPS ARE GIVEN AT N1, N2, ', + 'AND N3, IN DEGREES FROM HORIZONTAL;'/ + ' POSITIVE DIPS ARE TOWARD N1, N2, AND N3, RESPECTIVELY, '/ + ' WHILE NEGATIVE DIPS ARE TOWARD N6, N5, AND N4.)'/ + ' (THE ARGUMENT OF THE FAULT TRACE IS GIVEN AT N1 AND N3,'/ + ' IN DEGREES COUNTERCLOCKWISE FROM THE X AXIS.)'/ / + ' ELEMENT N1 N2 N3 N4 N5 N6', + ' DIP1 DIP2 DIP3 ARG1 ARG3 OFFSET'/) 240 FORMAT (' ',I5,6I5,3F5.0,2F7.1,E10.2) DO 300 K=1,NFL READ (IUNIT7,*) I,(IFN(J),J=1,6),(DIPS(L),L=1,3), + AZ1,AZ3,OFF CHECKF(I)=.TRUE. IF (.NOT.BRIEF) WRITE (IUNIT8,240) I,(IFN(J),J=1,6), + (DIPS(L),L=1,3),AZ1,AZ3,OFF DO 250 J=1,6 N=IFN(J) IF ((N.LE.0).OR.(N.GT.NTOP).OR. + ((N.GT.NREALN).AND.(N.LE.N1000))) THEN WRITE (IUNIT8,125) N STOP ENDIF IF (N.GT.NREALN) N=N-N1000+NREALN NODEF(J,I)=N 250 CONTINUE DO 260 L=1,3 IF (ABS(DIPS(L)).GT.90.) THEN WRITE(IUNIT8,252) DIPS(L) 252 FORMAT(' ILLEGAL DIP OF ',F10.4,'; SHOULD BE IN', + ' RANGE OF -90. TO +90. DEGREES.'/ + ' (NOTE: ALL DIPS ARE IN DEGREES FROM THE', + ' HORIZONAL;'/ + ' A + PREFIX (OR NONE) INDICATES A DIP', + ' TOWARD THE N1-N2-N3 SIDE;'/ + ' A - PREFIX INDICATES A DIP TOWARD', + ' THE N6-N5-N4 SIDE.)') STOP ENDIF IF (DIPS(L).LT.0.) DIPS(L)=180.+DIPS(L) FDIP(L,I)=DIPS(L)*0.017453293 260 CONTINUE IF ((ABS(AZ1).GT.361.).OR.(ABS(AZ3).GT.361.)) THEN WRITE(IUNIT8,272) AZ1,AZ3 272 FORMAT(' ILLEGAL ARGUMENT OF ',F10.4,' OR ',F10.4, + '; SHOULD BE IN RANGE -360. TO +360. DEGREES.') STOP ENDIF FAZ(1,I)=AZ1*0.017453293 FAZ(2,I)=AZ3*0.017453293 OFFSET(I)=OFF 300 CONTINUE ALLOK=.TRUE. DO 301 I=1,NFL ALLOK=ALLOK.AND.CHECKF(I) 301 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,302) 302 FORMAT(' THE FOLLOWING FAULTS WERE NEVER READ:') DO 304 I=1,NFL IF (.NOT.CHECKF(I)) WRITE(IUNIT8,303)I 303 FORMAT(' ',36X,I6) 304 CONTINUE STOP ENDIF IF (.NOT. BRIEF) WRITE (IUNIT8,999) 999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END SUBROUTINE GETNET C C C SUBROUTINE PUTNET (INPUT,IUNITO, + BRIEF,DQDTDA,ELEV,FAZ,FDIP, + MXEL,MXFEL,MXNODE,N1000, + NFAKEN,NFL,NODEF,NODES, + NREALN,NUMEL,NUMNOD,OFFSET, + TITLE1,TLNODE,XNODE,YNODE,ZMNODE) C C WRITES FINITE ELEMENT GRID TO UNIT "IUNITO". C CHARACTER*80 TITLE1 LOGICAL BRIEF C DIMENSION DQDTDA(MXNODE),ELEV(MXNODE), + FAZ(2,MXFEL),FDIP(3,MXFEL),NP(6), + NODEF(6,MXFEL),NODES(6,MXEL), + OFFSET(MXFEL),TLNODE(MXNODE), + XNODE(MXNODE),YNODE(MXNODE),ZMNODE(MXNODE) DIMENSION DIPS(3),ARGS(2) C WRITE (IUNITO,1) TITLE1 1 FORMAT (A80) C WRITE (IUNITO,2) NUMNOD, NREALN,NFAKEN,N1000,BRIEF 2 FORMAT(4I8,L8,' (NUMNOD, NREALN, NFAKEN, N1000, BRIEF)') C DO 100 I=1,NUMNOD IF (I.LE.NREALN) THEN IPRINT=I ELSE IPRINT=N1000+(I-NREALN) ENDIF X=XNODE(I) Y=YNODE(I) IF ((X.NE.0.).AND.(Y.NE.0.)) THEN WRITE (IUNITO,91) I,X,Y,ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 91 FORMAT (I5,1P,2E13.5,4E10.2) ELSE IF (X.NE.0.) THEN WRITE (IUNITO,92) I,X, ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 92 FORMAT (I5,1P, E13.5,' 0. ',4E10.2) ELSE IF (Y.NE.0.) THEN WRITE (IUNITO,93) I, Y,ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 93 FORMAT (I5,1P,' 0. ',E13.5,4E10.2) ELSE WRITE (IUNITO,94) I, ELEV(I),DQDTDA(I),ZMNODE(I), + TLNODE(I) 94 FORMAT (I5,1P,' 0. ',' 0. ', + 4E10.2) ENDIF 100 CONTINUE C WRITE (IUNITO,110) NUMEL 110 FORMAT (I10,' (NUMEL = NUMBER OF TRIANGULAR CONTINUUM ELEMENTS)') C DO 200 I=1,NUMEL DO 150 K=1,6 IF (NODES(K,I).LE.NREALN) THEN NP(K)=NODES(K,I) ELSE NP(K)=N1000+(NODES(K,I)-NREALN) ENDIF 150 CONTINUE WRITE (IUNITO,160) I,(NP(K),K=1,6) 160 FORMAT (I5,6I5) 200 CONTINUE C WRITE (IUNITO,210) NFL 210 FORMAT (I10,' (NFL = NUMBER OF CURVILINEAR FAULT ELEMENTS)') C DO 300 I=1,NFL DO 220 K=1,6 IF (NODEF(K,I).LE.NREALN) THEN NP(K)=NODEF(K,I) ELSE NP(K)=N1000+(NODEF(K,I)-NREALN) ENDIF 220 CONTINUE ARGS(1)=FAZ(1,I)*57.2958 ARGS(2)=FAZ(2,I)*57.2958 DO 230 K=1,3 DIPS(K)=FDIP(K,I) DIPS(K)=DIPS(K)*57.2958 IF (DIPS(K).GT.90.01) DIPS(K)=DIPS(K)-180. 230 CONTINUE WRITE (IUNITO,250) I,(NP(K),K=1,6), + (DIPS(K),K=1,3),(ARGS(K),K=1,2),OFFSET(I) 250 FORMAT (I5,6I6,3F5.0,2F7.1,1P,E10.2) 300 CONTINUE C RETURN END SUBROUTINE PUTNET C C C SUBROUTINE MIXUPX (INPUT,INEW,MXNODE,NUMNOD, + MODIFY,DATA, + WORK,ATNODE) C C SORT FLOATING-POINT SCALAR VALUES IN ARRAY "DATA" SO THAT THE C ENTRY CURRENTLY IN POSITION I ENDS UP IN POSITION "INEW(I)". C DIMENSION ATNODE(MXNODE),INEW(MXNODE), + DATA(MXNODE) C DO 10 I=1,NUMNOD ATNODE(I)=DATA(I) 10 CONTINUE C DO 20 I=1,NUMNOD DATA(INEW(I))=ATNODE(I) 20 CONTINUE C RETURN END SUBROUTINE MIXUPX C C C SUBROUTINE MIXUPI (INPUT,INEW,MXEL,MXNODE,NUMEL, + MODIFY,NODES) C C CHANGE 6-ENTRY INTEGER SETS IN ARRAY "NODES" SO THAT THE C INTEGER WHICH IS CURRENTLY "I" IS CHANGED TO "INEW(I)". C DIMENSION INEW(MXNODE), + NODES(6,MXEL) C DO 20 K=1,6 DO 10 J=1,NUMEL NODES(K,J)=INEW(NODES(K,J)) 10 CONTINUE 20 CONTINUE C RETURN END SUBROUTINE MIXUPI C C C SUBROUTINE NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) C C DETERMINE WHETHER THERE ARE MORE ELEMENTS ADJACENT TO SIDE J OF C TRIANGULAR CONTINUUM ELEMENT #I. C J = 1 MEANS THE SIDE OPPOSITE NODE # NODES(1,I). C J = 2 MEANS THE SIDE OPPOSITE NODE # NODES(2,I). C J = 3 MEANS THE SIDE OPPOSITE NODE # NODES(3,I). C IF A FAULT ELEMENT IS ADJACENT, ITS NUMBER IS KFAULT; OTHERWISE, C KFAULT IS SET TO ZERO. C IF ANOTHER TRIANGULAR CONTINUUM ELEMENT IS ADJACENT (EVEN ACROSS C FAULT ELEMENT KFAULT!) THEN ITS NUMBER IS KELE; OTHERWISE, KELE = 0. C LOGICAL FOUNDE,FOUNDF DIMENSION NODEF(6,MXFEL),NODES(6,MXEL) C C THREE NODE NUMBERS ALONG THE SIDE OF INTEREST, COUNTERCLOCKWISE: N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J, 3)+4,I) N3=NODES(MOD(J+1,3)+1,I) C CHECK FOR ADJACENT FAULT ELEMENT FIRST: FOUNDF=.FALSE. DO 10 K=1,NFL M1=NODEF(1,K) M2=NODEF(2,K) M3=NODEF(3,K) M4=NODEF(4,K) M5=NODEF(5,K) M6=NODEF(6,K) IF (((M1.EQ.N3).AND.(M2.EQ.N2).AND.(M3.EQ.N1)).OR. + ((M4.EQ.N3).AND.(M5.EQ.N2).AND.(M6.EQ.N1))) THEN FOUNDF=.TRUE. KFAULT=K GO TO 11 ENDIF 10 CONTINUE 11 IF (.NOT.FOUNDF) KFAULT=0 C IF THERE WAS A FAULT, REPLACE 3 NODE NUMBERS THAT WE SEARCH FOR: IF (FOUNDF) THEN IF (M2.EQ.N2) THEN N1=M4 N2=M5 N3=M6 ELSE N1=M1 N2=M2 N3=M3 ENDIF ENDIF C SEARCH FOR ADJACENT TRIANGULAR CONTINUUM ELEMENT: FOUNDE=.FALSE. DO 20 K=1,NUMEL IF (K.NE.I) THEN DO 15 L=1,3 M1=NODES(MOD(L, 3)+1,K) M2=NODES(MOD(L, 3)+4,K) M3=NODES(MOD(L+1,3)+1,K) IF ((M3.EQ.N1).AND.(M2.EQ.N2).AND.(M1.EQ.N3)) THEN FOUNDE=.TRUE. KELE=K GO TO 21 ENDIF 15 CONTINUE ENDIF 20 CONTINUE 21 IF (.NOT.FOUNDE) KELE=0 RETURN END SUBROUTINE NEXT C C C REAL FUNCTION ATAN2F (Y,X) C C CORRECTS FOR PROBLEM OF TWO ZERO ARGUMENTS C IF ((Y.NE.0.).OR.(X.NE.0.)) THEN ATAN2F=ATAN2(Y,X) ELSE ATAN2F=0. ENDIF RETURN END FUNCTION ATAN2F C C C SUBROUTINE SQUARE (INPUT,BRIEF,FAZ,FDIP,IUNIT8, + MXEL,MXFEL,MXNODE, + NFL,NODEF,NODES,NREALN, + NUMEL,NUMNOD,N1000, + MODIFY,XNODE,YNODE, + OUTPUT,NCOND,NODCON,NODTYP, + WORK,CHECKN,EDGEFS,EDGETS,LIST) C C CHECK, CORRECT, AND COMPLETE THE GEOMETRY OF THE GRID C LOGICAL AGREED,ALLOK,BRIEF,MATCH C C NOTE: THE FOLLOWING TYPE COULD BE "LOGICAL*1" IN IBM VS-FORTRAN: LOGICAL*1 CHECKN,EDGEFS,EDGETS C C NOTE: THE FOLLOWING COULD BE MADE "INTEGER*2" IN VS-FORTRAN: INTEGER*2 NODTYP C DOUBLE PRECISION FPHI,FPOINT,FGAUSS COMMON /SFAULT/ FPOINT COMMON /FPHIS/ FPHI COMMON /FGLIST/ FGAUSS DIMENSION FPHI(6,7),FPOINT(7),FGAUSS(7) DIMENSION CHECKN(MXNODE), + EDGEFS(2,MXFEL),EDGETS(3,MXEL),FDIP(3,MXFEL), + FAZ(2,MXFEL), + LIST(MXNODE),NODCON(MXNODE), + NODEF(6,MXFEL),NODES(6,MXEL),NODTYP(MXNODE), + XNODE(MXNODE),YNODE(MXNODE) C C (1) CHECK THAT ALL REAL NODES ARE CONNECTED TO AT LEAST ONE C CONTINUUM (TRIANGULAR) ELEMENT; C DO 110 I=1,NREALN CHECKN(I)=.FALSE. 110 CONTINUE DO 130 I=1,NUMEL DO 120 J=1,6 CHECKN(NODES(J,I))=.TRUE. 120 CONTINUE 130 CONTINUE DO 132 I=1,NFL DO 131 J=1,6 CHECKN(NODEF(J,I))=.TRUE. 131 CONTINUE 132 CONTINUE ALLOK=.TRUE. DO 140 I=1,NREALN ALLOK=ALLOK.AND.CHECKN(I) 140 CONTINUE IF (.NOT.ALLOK) THEN WRITE(IUNIT8,150) 150 FORMAT(' BAD GRID TOPOLOGY: FOLLOWING REAL NODES DO NOT'/ 1 ' BELONG TO ANY CONTINUUM ELEMENT OR FAULT:') DO 160 I=1,NREALN IF (.NOT.CHECKN(I)) WRITE (IUNIT8,155) I 155 FORMAT (' ',43X,I6) 160 CONTINUE STOP ENDIF C C (2) CHECK THAT EVERY NODE IS EITHER A CORNER OR A MIDPOINT NODE, C BUT NOT BOTH. C DO 210 I=1,NUMNOD NODTYP(I)=0 210 CONTINUE NTOFIX=0 ALLOK=.TRUE. DO 250 I=1,NUMEL DO 240 J=1,6 IF (J.LE.3) THEN ITYPE=1 ELSE ITYPE=2 ENDIF N=NODES(J,I) IF (NODTYP(N).EQ.0) THEN NODTYP(N)=ITYPE IF (ITYPE.EQ.2) THEN IF ((XNODE(N).EQ.0.).AND.(YNODE(N).EQ.0.)) + NTOFIX=NTOFIX+1 ENDIF ELSE IF (NODTYP(N).NE.ITYPE) THEN ALLOK=.FALSE. WRITE (IUNIT8,220) N 220 FORMAT(' BAD GRID TOPOLOGY: NODE ',I6, + ' CANNOT BE AN ELEMENT CORNER AND AN', + ' ELEMENT SIDE-MIDPOINT AT THE SAME', + ' TIME.') ENDIF ENDIF 240 CONTINUE 250 CONTINUE DO 290 I=1,NFL DO 280 J=1,6 IF ((J.EQ.2).OR.(J.EQ.5)) THEN ITYPE=2 ELSE ITYPE=1 ENDIF N=NODEF(J,I) IF (NODTYP(N).EQ.0) THEN NODTYP(N)=ITYPE ELSE IF (NODTYP(N).NE.ITYPE) THEN ALLOK=.FALSE. WRITE (IUNIT8,220) N ENDIF ENDIF 280 CONTINUE 290 CONTINUE IF (.NOT.ALLOK) STOP C C (3) CHECK THAT EACH FAULT SIDE WITH REAL NODES ALONG IT SHARES C THOSE SAME 3 NODES WITH A TRIANGULAR CONTINUUM ELEMENT. C ALLOK=.TRUE. DO 390 I=1,NFL DO 380 J=2,5,3 N=NODEF(J,I) IF (N.LE.NREALN) THEN DO 320 K=1,NUMEL DO 310 L=4,6 IF (NODES(L,K).EQ.N) THEN LP=L-2 IF (LP.EQ.4) LP=1 LM=L-3 MATCH=((NODEF(J-1,I).EQ.NODES(LP,K)) + .OR.(NODEF(J-1,I).GT.NREALN)) + .AND.((NODEF(J+1,I).EQ.NODES(LM,K)) + .OR.(NODEF(J+1,I).GT.NREALN)) IF (.NOT.MATCH) THEN ALLOK=.FALSE. WRITE(IUNIT8,305) I,K 305 FORMAT(' BAD GRID TOPOLOGY:', + ' FAULT ',I6,' IS NOT PROPERL' + ,'Y CONNECTED TO ELEMENT ',I6) ELSE GO TO 380 ENDIF ENDIF 310 CONTINUE 320 CONTINUE ENDIF 380 CONTINUE 390 CONTINUE IF (.NOT.ALLOK) STOP C C (4) AVERAGE TOGETHER THE COORDINATES OF ALL NODES AT ONE "POINT" C DO 410 I=1,NUMNOD CHECKN(I)=.FALSE. C (MEANS "NOT YET INVOLVED IN AVERAGING') 410 CONTINUE DO 490 I=1,NFL DO 480 J1=1,3,2 NJ1=NODEF(J1,I) C (FAULT ENDS ARE THE ONLY PLACES THAT CAN HAVE PROBLEMS) IF (.NOT.CHECKN(NJ1)) THEN LIST(1)=NJ1 CHECKN(NJ1)=.TRUE. C BEGIN LIST OF NEIGHBORS WITH PAIRED NODE J2=7-J1 NJ2=NODEF(J2,I) LIST(2)=NJ2 CHECKN(NJ2)=.TRUE. NINSUM=2 C FIND SHORTEST FAULT CONNECTED TO EITHER ONE SHORT=SQRT( + (XNODE(NODEF(1,I))-XNODE(NODEF(3,I)))**2+ + (YNODE(NODEF(1,I))-YNODE(NODEF(3,I)))**2) DO 470 K=1,NFL NL1=NODEF(1,K) NL3=NODEF(3,K) NL4=NODEF(4,K) NL6=NODEF(6,K) IF ((NJ1.EQ.NL1).OR.(NJ2.EQ.NL1).OR. + (NJ1.EQ.NL3).OR.(NJ2.EQ.NL3).OR. + (NJ1.EQ.NL4).OR.(NJ2.EQ.NL4).OR. + (NJ1.EQ.NL6).OR.(NJ2.EQ.NL6)) THEN TEST=SQRT( + (XNODE(NL1)-XNODE(NL3))**2+ + (YNODE(NL1)-YNODE(NL3))**2) SHORT=MIN(SHORT,TEST) ENDIF 470 CONTINUE C COLLECT ALL CORNER NODES WITHIN 10% OF THIS TOLER=SHORT/10. T2=TOLER**2 DO 471 K=1,NUMNOD IF (.NOT.CHECKN(K)) THEN IF (NODTYP(K).EQ.1) THEN R2=(XNODE(NJ1)-XNODE(K))**2+ + (YNODE(NJ1)-YNODE(K))**2 IF (R2.LT.T2) THEN NINSUM=NINSUM+1 LIST(NINSUM)=K CHECKN(K)=.TRUE. ENDIF ENDIF ENDIF 471 CONTINUE C (QUICK EXIT IF ALL NODES IN SAME PLACE) AGREED=.TRUE. DO 472 K=2,NINSUM AGREED=AGREED.AND. + (XNODE(K).EQ.XNODE(1)).AND. + (YNODE(K).EQ.YNODE(1)) 472 CONTINUE IF (AGREED) GO TO 480 XSUM=0. YSUM=0. DO 473 K=1,NINSUM XSUM=XSUM+XNODE(LIST(K)) YSUM=YSUM+YNODE(LIST(K)) 473 CONTINUE XMEAN=XSUM/NINSUM YMEAN=YSUM/NINSUM RMAX=0. DO 474 K=1,NINSUM R=SQRT((XNODE(LIST(K))-XMEAN)**2+ + (YNODE(LIST(K))-YMEAN)**2) RMAX=MAX(RMAX,R) 474 CONTINUE DO 475 K=1,NINSUM XNODE(LIST(K))=XMEAN YNODE(LIST(K))=YMEAN 475 CONTINUE IF (.NOT.BRIEF) THEN IF (RMAX.GT.0.) THEN WRITE(IUNIT8,476) NINSUM, + (LIST(N),N=1,NINSUM) 476 FORMAT(/ + ' AVERAGING TOGETHER THE POSITIONS OF', + ' THESE ',I6,' NODES:',(/' ',12I6)) WRITE (IUNIT8,477) RMAX 477 FORMAT (' MAXIMUM CORRECTION TO ', + 'ANY POSITION IS',1P,E10.2/ + ' YOU ARE RESPONSIBLE FOR ', + ' DECIDING WHETHER THIS IS A', + ' SERIOUS ERR0R!') ENDIF ENDIF ENDIF 480 CONTINUE 490 CONTINUE C C (6) COMPUTE COORDINATES OF MIDPOINT NODES THAT WERE NOT INPUT. C C FIRST, FAULTS: IF ((.NOT.BRIEF).AND.(NFL.GT.0)) WRITE (IUNIT8,540) 540 FORMAT(/ /' FOLLOWING FAULT MID-POINT POSITIONS WERE COMPUTED:'/ + /' FAULT NODE2 NODE5 X Y'/) DO 550 I=1,NFL I1=NODEF(1,I) I2=NODEF(2,I) I3=NODEF(3,I) I5=NODEF(5,I) DX= XNODE(I3)- XNODE(I1) DY= YNODE(I3)- YNODE(I1) AZ=ATAN2(DY,DX) PHI1=FAZ(1,I)-AZ PHI1=MOD(PHI1+1.570796,3.14159265)-1.570796 PHI2=AZ-FAZ(2,I) PHI2=MOD(PHI2+1.570796,3.14159265)-1.570796 IF ((ABS(PHI1).GT.0.001).OR.(ABS(PHI2).GT.0.001)) THEN T1=TAN(PHI1) T2=TAN(PHI2) IF (ABS(T2-T1).GE.ABS(T1+T2)) THEN FACTOR=0.99*ABS(T1+T2)/ABS(T2-T1) IF (ABS(T1).GT.ABS(T2)) THEN T2=T1+FACTOR*(T2-T1) ELSE T1=T2+FACTOR*(T1-T2) ENDIF ENDIF PARRAL=(T2-T1)/(4.*(T1+T2)) PERPEN= T1*T2 /(2.*(T1+T2)) XNODE(I2)=XNODE(I1)+DX/2.+PARRAL*DX-PERPEN*DY YNODE(I2)=YNODE(I1)+DY/2.+PERPEN*DX+PARRAL*DY ELSE XNODE(I2)=(XNODE(I1)+XNODE(I3))/2. YNODE(I2)=(YNODE(I1)+YNODE(I3))/2. ENDIF XNODE(I5)= XNODE(I2) YNODE(I5)= YNODE(I2) NTOFIX=NTOFIX-1 IF (.NOT.BRIEF) WRITE (IUNIT8,549) I,I2,I5,XNODE(I2), 1 YNODE(I2) 549 FORMAT(' ',I6,2I10,1P,2E12.4) 550 CONTINUE C C NEXT, OTHER ELEMENT SIDES, IF NEEDED: IF ((.NOT.BRIEF).AND.(NTOFIX.GT.0)) WRITE (IUNIT8,551) 551 FORMAT(/ /' FOLLOWING MID-POINTS OF CONTINUUM ELEMENT SIDES', + ' THAT WERE 0.0 IN THE' + / ' INPUT DATASET ARE NOW COMPUTED, AS FOLLOWS:' + / / ' ELEMENT NODE X Y'/) DO 590 I=1,NUMEL DO 580 J=4,6 N=NODES(J,I) IF ((XNODE(N).EQ.0.).AND.(YNODE(N).EQ.0.)) THEN JP=J-2 IF (J.EQ.6) JP=1 JM=J-3 XNODE(N)=0.5* + (XNODE(NODES(JP,I))+XNODE(NODES(JM,I))) YNODE(N)=0.5* + (YNODE(NODES(JP,I))+YNODE(NODES(JM,I))) IF (.NOT.BRIEF) + WRITE (IUNIT8,579) I,N,XNODE(N),YNODE(N) 579 FORMAT(' ',I6,I10,1P,2E12.4) ENDIF 580 CONTINUE 590 CONTINUE C C (9) MAKE A LIST OF NODES THAT ARE ON THE BOUNDARY AND REQUIRE C BOUNDARY CONDITIONS (NODCON); THESE ARE IN COUNTERCLOCKWISE C ORDER. ALSO MAKE A LISTS OF ELEMENT SIDES WHICH CONTAIN THESE C NODES: EDGETS AND EDGEFS. C WRITE (IUNIT8,800) 800 FORMAT (/' ----------------------------------------' + / /' COMPILING AN ORDERED LIST OF BOUNDARY NODES...'/) NCOND=0 DO 801 I=1,NUMNOD CHECKN(I)=.FALSE. 801 CONTINUE DO 802 I=1,NFL EDGEFS(1,I)=.FALSE. EDGEFS(2,I)=.FALSE. 802 CONTINUE DO 810 I=1,NUMEL DO 809 J=1,3 CALL NEXT (INPUT,I,J,MXEL,MXFEL,NFL,NODEF,NODES,NUMEL, + OUTPUT,KFAULT,KELE) IF (KELE.GT.0) THEN C (ORDINARY INTERIOR SIDE) EDGETS(J,I)=.FALSE. ELSE IF (KFAULT.EQ.0) THEN C (EXTERIOR SIDE) EDGETS(J,I)=.TRUE. N1=NODES(MOD(J, 3)+1,I) N2=NODES(MOD(J, 3)+4,I) N3=NODES(MOD(J+1,3)+1,I) IF (.NOT.CHECKN(N1)) THEN NCOND=NCOND+1 CHECKN(N1)=.TRUE. ENDIF IF (.NOT.CHECKN(N2)) THEN NCOND=NCOND+1 CHECKN(N2)=.TRUE. ENDIF IF (.NOT.CHECKN(N3)) THEN NCOND=NCOND+1 CHECKN(N3)=.TRUE. ENDIF ELSE C (TRIANGULAR ELEMENT HAS AN EXTERIOR FAULT ELEMENT C ADJACENT TO IT) EDGETS(J,I)=.FALSE. N2=NODES(MOD(J, 3)+4,I) IF (NODEF(2,KFAULT).EQ.N2) THEN EDGEFS(2,KFAULT)=.TRUE. DO 806 K=4,6 N=NODEF(K,KFAULT) IF (.NOT.CHECKN(N)) THEN NCOND=NCOND+1 CHECKN(N)=.TRUE. ENDIF 806 CONTINUE ELSE EDGEFS(1,KFAULT)=.TRUE. DO 808 K=1,3 N=NODEF(K,KFAULT) IF (.NOT.CHECKN(N)) THEN NCOND=NCOND+1 CHECKN(N)=.TRUE. ENDIF 808 CONTINUE ENDIF ENDIF 809 CONTINUE 810 CONTINUE IF (NUMNOD.GT.NREALN) THEN DO 824 I=NREALN+1,NUMNOD IF (.NOT.CHECKN(I)) THEN IO=N1000+I-NREALN WRITE(IUNIT8,822) IO 822 FORMAT(' BAD GRID TOPOLOGY; FAKE NODES ARE NOT', + ' PERMITTED IN THE INTERIOR.'/' CHECK NODE ',I6) STOP ENDIF 824 CONTINUE ENDIF C BEGIN CIRCUIT WITH LOWEST-NUMBERED BOUNDARY NODE DO 830 I=1,NUMNOD IF (CHECKN(I)) GO TO 831 830 CONTINUE 831 NODCON(1)=I NDONE=1 NLEFT=NCOND C BEGINNING OF INDEFINATE LOOP WHICH TRACES AROUND THE PERIMETER. C EACH TIME, IT PROGRESSES BY ONE OF 3 STEPS: C -2 NODES AT A TIME ALONG A TRIANGLE SIDE, OR C -2 NODES AT A TIME ALONG A FAULT ELEMENT SIDE, OR C -BY FINDING ANOTHER (CORNER) NODE WHICH SHARES THE SAME LOCATION. C FIRST, BE SURE THAT WE ARE NOT STARTING ON A MIDPOINT: IF (NODTYP(I).EQ.2) THEN DO 833 K=1,NUMEL DO 832 L=1,3 IF (EDGETS(L,K)) THEN N2=NODES(MOD(L, 3)+4,K) IF (N2.EQ.I) THEN J=NODES(MOD(L+1,3)+1,K) GO TO 839 ENDIF ENDIF 832 CONTINUE 833 CONTINUE DO 835 K=1,NFL IF (EDGEFS(1,K)) THEN IF (NODEF(2,K).EQ.I) THEN J=NODEF(3,K) GO TO 839 ENDIF ELSE IF (EDGEFS(2,K)) THEN IF (NODEF(5,K).EQ.I) THEN J=NODEF(6,K) GO TO 839 ENDIF ENDIF 835 CONTINUE 839 NDONE=2 NODCON(2)=J NLEFT=NCOND-1 ENDIF C BEGINNING OF MAIN INDEFINATE LOOP: 840 NODE=NODCON(NDONE) X=XNODE(NODE) Y=YNODE(NODE) C LOOK FOR AN ADJACENT TRIANGULAR ELEMENT USING THIS NODE. DO 844 I=1,NUMEL DO 842 J=1,3 IF (EDGETS(J,I)) THEN N1=NODES(MOD(J,3)+1,I) IF (N1.EQ.NODE) GO TO 846 ENDIF 842 CONTINUE 844 CONTINUE GO TO 850 846 N2=NODES(MOD(J,3)+4,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. N3=NODES(MOD(J+1,3)+1,I) NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N3 CHECKN(N3)=.FALSE. NLEFT=NLEFT-2 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR AN ADJACENT FAULT ELEMENT USING THIS NODE. 850 DO 854 I=1,NFL IF (EDGEFS(1,I)) THEN IF (NODEF(1,I).EQ.NODE) THEN N2=NODEF(2,I) N3=NODEF(3,I) GO TO 856 ENDIF ELSE IF (EDGEFS(2,I)) THEN IF (NODEF(4,I).EQ.NODE) THEN N2=NODEF(5,I) N3=NODEF(6,I) GO TO 856 ENDIF ENDIF 854 CONTINUE GO TO 860 856 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N2 CHECKN(N2)=.FALSE. NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=N3 CHECKN(N3)=.FALSE. NLEFT=NLEFT-2 IF (NLEFT.GT.0) THEN GO TO 840 ELSE GO TO 870 ENDIF C ELSE, LOOK FOR ANOTHER EXTERIOR CORNER NODE AT SAME LOCATION. 860 DO 865 I=1,NUMNOD IF ((I.NE.NODE).AND.CHECKN(I)) THEN IF ((NODTYP(I).EQ.1).AND. + ((XNODE(I).EQ.X).AND.(YNODE(I).EQ.Y)))GO TO 867 ENDIF 865 CONTINUE WRITE(IUNIT8,866) NDONE, NODE, X, Y 866 FORMAT(/' AFTER CONNECTING ',I6,' NODES AROUND THE', + ' PERIMETER, '/ + ' PROCESS WAS STOPPED BY BAD GRID TOPOLOGY;'/ + ' COULD NOT FIND ANY WAY TO CONTINUE FROM NODE ',I6/ + ' AT (X=',1P,E10.3,',Y=',E10.3,')'/ + ' EITHER THROUGH SHARED BOUNDARY ELEMENTS, OR'/ + ' THROUGH OTHER BOUNDARY NODES SHARING THE SAME ', + 'POSITION.') STOP 867 NDONE=NDONE+1 IF (NDONE.LE.NCOND) NODCON(NDONE)=I CHECKN(I)=.FALSE. NLEFT=NLEFT-1 IF (NLEFT.GT.0) GO TO 840 C END OF INDEFINATE LOOP WHICH TRACES AROUND PERIMETER. 870 IF (.NOT.BRIEF) THEN WRITE(IUNIT8,880) 880 FORMAT(/ /' HERE FOLLOWS A LIST, IN CONSECUTIVE ORDER,'/ + ' OF THE NODES WHICH DEFINE THE PERIMETER'/ + ' OF THE MODEL; THESE NODES REQUIRE BOUNDARY', + ' CONDITIONS:'/' BC# NODE') DO 890 I=1,NCOND N=NODCON(I) IF (N.GT.NREALN) N=N1000+N-NREALN WRITE(IUNIT8,882) I, N 882 FORMAT(' ',2I6) 890 CONTINUE N=NODCON(1) IF (N.GT.NREALN) N=N1000+N-NREALN WRITE (IUNIT8,892) N 892 FORMAT(' (NOTE: NODE ',I6,' COMPLETES THE LOOP, BUT WILL', + ' NOT BE LISTED TWICE.)') ENDIF C IF (.NOT. BRIEF) WRITE (IUNIT8,9999) 9999 FORMAT (' --------------------------------------------------', + '-----------------------------') RETURN END SUBROUTINE SQUARE C C C SUBROUTINE READN (INPUT,IUNITP,IUNITT,N, + OUTPUT,VECTOR) C C A UTILITY ROUTINE DESIGNED TO PERMIT -FAULTS- INPUT FILES C TO ALSO BE USED BY -PLATES-, WHICH EXPECTS MORE NUMBERS C IN SOME RECORDS. C THIS ROUTINE ATTEMPTS TO READ 'N' FLOATING-POINT VALUES C (USING * FORMAT) FROM THE NEXT RECORD ON DEVICE 'IUNITP'. C IF ANYTHING GOES WRONG, THE MISSING VALUES ARE SET TO ZERO. C CHARACTER*1 C CHARACTER*80 LINE INTEGER IOS LOGICAL ANYIN,DOTTED,EXPON,SIGNED DIMENSION VECTOR(N) C LINE=' '// + ' ' READ (IUNITP,1,IOSTAT=IOS) LINE 1 FORMAT (A80) C NUMBER=0 ANYIN=.FALSE. EXPON=.FALSE. SIGNED=.FALSE. DOTTED=.FALSE. DO 10 I=1,80 C=LINE(I:I) IF ((C.EQ.' ').OR.(C.EQ.',').OR.(C.EQ.'/')) THEN SIGNED=.FALSE. EXPON=.FALSE. DOTTED=.FALSE. IF (ANYIN) THEN NUMBER=NUMBER+1 ANYIN=.FALSE. ENDIF ELSE IF ((C.EQ.'+').OR.(C.EQ.'-')) THEN IF (SIGNED) THEN GO TO 50 ELSE SIGNED=.TRUE. ENDIF ELSE IF ((C.EQ.'E').OR.(C.EQ.'D').OR. + (C.EQ.'e').OR.(C.EQ.'d')) THEN IF (EXPON) THEN GO TO 50 ELSE EXPON=.TRUE. SIGNED=.FALSE. DOTTED=.TRUE. ENDIF ELSE IF (C.EQ.'.') THEN IF (DOTTED) THEN GO TO 50 ELSE DOTTED=.TRUE. ENDIF ELSE IF ((C.EQ.'0').OR.(C.EQ.'1').OR.(C.EQ.'2').OR. + (C.EQ.'3').OR.(C.EQ.'4').OR.(C.EQ.'5').OR. + (C.EQ.'6').OR.(C.EQ.'7').OR.(C.EQ.'8').OR. + (C.EQ.'9')) THEN SIGNED=.TRUE. ANYIN=.TRUE. ELSE GO TO 50 ENDIF 10 CONTINUE IF (ANYIN) NUMBER=NUMBER+1 C 50 IF (NUMBER.EQ.0) THEN WRITE (IUNITT,91) N,LINE 91 FORMAT (/' ERR0R: A LINE OF PARAMETER INPUT WHICH', + ' WAS SUPPOSED TO CONTAIN 1-',I2,' NUMBERS'/ + ' COULD NOT BE INTERPRETED. LINE FOLLOWS:'/ + ' ',A80) STOP ELSE NUMBER=MIN(NUMBER,N) BACKSPACE IUNITP READ (IUNITP,*) (VECTOR(I),I=1,NUMBER) IF (NUMBER.LT.N) THEN DO 99 I=NUMBER+1,N VECTOR(I)=0. 99 CONTINUE ENDIF ENDIF RETURN END SUBROUTINE READN