PROGRAM PATCRK C C This program converts a finite element mesh from the patran neutral file C format to the FRANC new analysis input format. C C Program description: C Unit 20 is used for the input file and unit 21 is used for the output C file. Units 5 and 6 are used for interactive i/o with the user, but C these are assigned to variables and can be easily changed near the top C of the program. C Because this program is intended to translate primarily geometric C data, many of the values input from PATRAN are ignored (ie the difference C between materials and properties, the element configuration, number of C degrees of freedom at each node, etc.). One that is used is the material C pointer. However if its value is zero, it defaults to 1. The element C configuration (type) number is always set to 1 (straight sides). C In FRANC configuration 1 refers to straight sided elements with C the midside nodes at the true midside. FRANC also requires the element C numbering for each element to be counter clockwise ( this produces C a positive area for the element). Both these requirements are handled. C One other feature included is a debug facility. If the program is C compiled with the DEBUG switch on, a line in the input subroutine is C set which allows a debug report to unit IOUT2. This requires that the user C preconnect this unit to a file for output. The program can be easily C modified to include this output as standard. C C Happy Translating! C C Programmer: Mark A. James C Version: 1.00 Date: November 1987 C Version: 1.01 Date: May 1988 C Version: 1.02 C C Changes: Fixed command line input for Apollo. Added prompts for input C if command line input is not present. Put record length in C unformatted OPEN statement. Removed HARRIS system dependencies. C Put in DO ... END DO loops. C C Modified to work for FRANC. C C Note: This program contains some nonstandard FORTRAN statements. C C C INTEGER EMAX, NMAX, MMAX C Emax is element max C Nmax is node max C Mmax is material max PARAMETER (NMAX = 10000, EMAX = 10000) PARAMETER (MMAX = 20) C CHARACTER*80 CMD_LINE ! Command line as input by user. CHARACTER*80 IN_FILE ! Patran neutral format input file name. CHARACTER*80 OUT_FILE ! Cracker output file name. CHARACTER DUM C INTEGER INF, OUTF ! Input and Output file name lengths. INTEGER IIN, IOUT ! Input and Output interactive unit. INTEGER I, J ! Loop counters. C CHARACTER*80 TITLE ! Problem title associated with mesh. LOGICAL EX ! Existance of file for I/O. INTEGER NUMNP ! Number of node points. INTEGER NDMAX ! Maximum number of node points. INTEGER NUMEL ! Number of elements. INTEGER ELMAX ! Maximum number of elements. INTEGER CONN(11,EMAX) ! Connectivity matrix. INTEGER TCONN(10) ! Connectivity matrix. INTEGER NUMMAT ! Number of material property types. INTEGER MAMAX ! Maximum number of materials. INTEGER MAT(MMAX) ! Used to determine # of unique mat's. INTEGER IDEBUG ! Debug switch > 0, debug file 22 is opepatc REAL V1(2), V2(2) C DOUBLE PRECISION X(NMAX),Y(NMAX) ! Global Nodal coordinates. C COMMON IIN, IOUT, IDEBUG C COMMON / NDAT / X, Y COMMON / EDAT / CONN C C C These next few statements are APOLLO insert files. C C%nolist c INCLUDE '/sys/ins/base.ins.ftn' c INCLUDE '/sys/ins/error.ins.ftn' c INCLUDE '/sys/ins/pgm.ins.ftn' C%list INTEGER*4 STATUS ! Error status returned from APOLLO callpa C End Apollo inserts. DATA NDMAX, ELMAX, MAMAX / NMAX, EMAX, MMAX / DATA IDEBUG / 0 / C C If the program is compiled with the debug switch on this line will be compilep C and the READ_PAT_NEUTRAL routine will output information to a debug file. c IDEBUG = 1 C C Zero the connectivity matrix and node location vectors. C DO I=1, 11 DO J=1, EMAX CONN(I,J) = 0 END DO END DO DO I=1,NMAX X(I) = 0.0 Y(I) = 0.0 END DO C C Default terminal I/O units C IIN = 5 IOUT = 6 IOUT2 = 6 C C This is not standard fortran C C Get input file name C INF = 0 c INF = PGM_$GET_ARG( INT2(1), IN_FILE, STATUS, INT2(80) ) ! Apollo systpatcrk C IF (INF.LE.0) THEN 100 CONTINUE WRITE(IOUT,'(x,A$)')'ENTER INPUT FILE NAME: ' READ(IIN,1000) CMD_LINE CALL LEFT(CMD_LINE) INF = INDEX(CMD_LINE,' ') IN_FILE = CMD_LINE(:INF) IF (INF.LE.0) GOTO 100 END IF C C Make sure input file exists C C Hacks because the DECstation doesn't do inquires C C INQUIRE(FILE = IN_FILE(:INF), EXIST = EX ) C IF (.NOT. EX) THEN C WRITE(IOUT,*) ' THE INPUT FILE DOES NOT EXIST!!' C STOP C END IF C C Get output file name C OUTF = 0 c OUTF = PGM_$GET_ARG( INT2(2), OUT_FILE, STATUS, INT2(80) ) ! Apollo systpatcrk C IF (OUTF.LE.0) THEN 200 CONTINUE WRITE(IOUT,'(x,A$)')'ENTER OUTPUT FILE NAME: ' READ(IIN,1000) CMD_LINE CALL LEFT(CMD_LINE) OUTF = INDEX(CMD_LINE,' ') - 1 OUT_FILE = CMD_LINE(:OUTF) IF (OUTF.LE.0) GOTO 200 END IF C C Check if output file exists C if ( index(out_file(:outf),'.inp') .EQ. 0 ) THEN OUT_FILE = OUT_FILE(:OUTF)//'.inp' OUTF = OUTF + 4 end if C Hacks because the DECstation doesn't do inquires C C INQUIRE (FILE=OUT_FILE(:OUTF), EXIST = EX ) C IF (EX) THEN C WRITE(IOUT,1000) ' OUTPUT FILE = '//OUT_FILE(:OUTF)// C _ ' ALREADY EXISTS' C C Check if they want to overwrite it C C300 CONTINUE C WRITE(IOUT,'(A$)') ' OVERWRITE (Y,N)? ' C READ(IIN,1000) DUM C IF (DUM.EQ.'N' .OR. DUM.EQ.'n') STOP C IF ( (DUM.NE.'Y' .AND. DUM.NE.'y') ) GOTO 300 C END IF C C Echo the input and output file names C WRITE(IOUT,1000) ' INPUT FILE : '//IN_FILE(:INF) WRITE(IOUT,1000) ' OUTPUT FILE: '//OUT_FILE(:OUTF) IF (IDEBUG.NE.0) THEN WRITE(IOUT,1000) ' DEBUG FILE : '//OUT_FILE(:OUTF-4)//'.dbg' END IF C C Open the input file and debug file C OPEN (UNIT = 20, FILE = IN_FILE(:INF), FORM = 'FORMATTED', _ ERR = 900, STATUS = 'OLD' ) IF (IDEBUG.NE.0) THEN OPEN (UNIT = 22, FILE = OUT_FILE(:OUTF-4)//'.dbg', _ FORM = 'FORMATTED', ERR = 902, STATUS='NEW') END IF C C Read the neutral file. C CALL READ_PAT_NEUTRAL(CONN,X,Y,NUMEL,NUMNP,NUMMAT,MAT, _ ELMAX,NDMAX,MAMAX,TITLE) C C Open the output file. C OPEN (UNIT = 21, FILE = OUT_FILE(:OUTF), _ FORM = 'FORMATTED', ERR = 901, STATUS='NEW') C C Output the data to FRANC format C WRITE(21,*) TITLE WRITE(21,*) NUMNP, NUMEL, NUMMAT, 2 c materials DO I=1,NUMMAT WRITE(21,'(I5, 4F10.2)') 1, 2900e3, 0.25, 1.0, 1.0 END DO c connectivity: check for clockwise elements and reverse them, but c first put them in franc format. DO K=1,NUMEL IF ( CONN(7,K) .NE. 0 ) THEN TCONN(1) = CONN(9,K) TCONN(2) = CONN(10,K) TCONN(3) = CONN(1,K) TCONN(4) = CONN(5,K) TCONN(5) = CONN(2,K) TCONN(6) = CONN(6,K) TCONN(7) = CONN(3,K) TCONN(8) = CONN(7,K) TCONN(9) = CONN(4,K) TCONN(10) = CONN(8,K) ELSE TCONN(1) = CONN(9,K) TCONN(2) = CONN(10,K) TCONN(3) = CONN(1,K) TCONN(4) = CONN(4,K) TCONN(5) = CONN(2,K) TCONN(6) = CONN(5,K) TCONN(7) = CONN(3,K) TCONN(8) = CONN(6,K) TCONN(9) = 0 TCONN(10) = 0 END IF V1(1) = X(TCONN(7)) - X(TCONN(5)) V1(2) = Y(TCONN(7)) - Y(TCONN(5)) V2(1) = X(TCONN(3)) - X(TCONN(5)) V2(2) = Y(TCONN(3)) - Y(TCONN(5)) XPROD = V1(1)*V2(2) - V2(1)*V1(2) IF ( XPROD .LT. 0.0 ) THEN IF ( TCONN(9) .NE. 0 ) THEN ITMP = TCONN(4) TCONN(4) = TCONN(10) TCONN(10) = ITMP ITMP = TCONN(5) TCONN(5) = TCONN(9) TCONN(9) = ITMP ITMP = TCONN(6) TCONN(6) = TCONN(8) TCONN(8) = ITMP ELSE ITMP = TCONN(4) TCONN(4) = TCONN(8) TCONN(8) = ITMP ITMP = TCONN(5) TCONN(5) = TCONN(7) TCONN(7) = ITMP END IF END IF WRITE(21,'(10I5)') (TCONN(I), I=1, 10) END DO DO I=1,NUMNP WRITE(21,'(I5, 2F10.2)') I, X(i), Y(i) END DO c WRITE(21) TITLE c WRITE(21) NUMEL, NUMNP, NUMMAT c WRITE(21) (X(I), I=1, NUMNP) c WRITE(21) (Y(I), I=1, NUMNP) c WRITE(21) ((CONN(I,J), I=1, 9), J=1, NUMEL) C C Close the input and FRANC files. C CLOSE (20) CLOSE (21) IF (IDEBUG.NE.0) THEN CLOSE (22) END IF C C Tell the user how much stuff was processed C WRITE(IOUT,'(A$)')'PROBLEM TITLE: ' WRITE(IOUT,1000) TITLE WRITE(IOUT,*) WRITE(IOUT,1002) NUMNP WRITE(IOUT,1003) NUMEL WRITE(IOUT,1004) NUMMAT C C This is it! C STOP 900 WRITE(IOUT,1000)' ERROR OPENING INPUT FILE.' STOP 901 WRITE(IOUT,1000)' ERROR OPENING OUTPUT FILE.' STOP 902 WRITE(IOUT,1000)' ERROR OPENING DEBUG FILE.' 1000 FORMAT(A) 1002 FORMAT(' NUMBER OF NODES:',I5) 1003 FORMAT(' NUMBER OF ELEM :',I5) 1004 FORMAT(' NUMBER OF MATER:',I5) END SUBROUTINE LEFT(LINE) CHARACTER*(*) LINE INTEGER I,J IF (LINE.NE.' ') THEN 100 CONTINUE I = INDEX (LINE,' ') IF (I.EQ.1) THEN LINE = LINE (2:) END IF IF ( .NOT.(I.EQ.0 .OR. I.NE.1) ) GOTO 100 END IF RETURN END SUBROUTINE READ_PAT_NEUTRAL(NTABLE,X,Y,NUMEL,NUMNP,NUMMAT,MAT, _ ELMAX, NDMAX, MAMAX, TITLE) INTEGER NUMEL INTEGER NUMNP ! (OUTPUT) CURRENT NUMBER OF NODES INTEGER NUMMAT ! (OUTPUT) TOTAL NUMBER OF MAT TYPES INTEGER ELMAX ! (INPUT) TOTAL ELEMENTS ALLOWED INTEGER NDMAX ! (INPUT) TOTAL NODES ALLOWED INTEGER MAMAX ! (INPUT) TOTAL NUMBER OF UNIQUE MAT'S INTEGER MAT(MAMAX) ! (OUTPUT) UNIQUE MAT TYPES INTEGER NTABLE(11,ELMAX) ! (OUTPUT) TABLE OF ELEMENT NODES (CONN) DOUBLE PRECISION X(NDMAX) ! (OUTPUT) GLOBAL NODE LOCATIONS DOUBLE PRECISION Y(NDMAX) DOUBLE PRECISION X1,Y1,Z1 ! DUMMY USED TO READ COORD IN DOUBLE PRECISION X2, X3, Y2, Y3 ! USED TO MAKE STRAIGHT SIDED ELEMENTS CHARACTER*80 TITLE ! (OUTPUT) MODEL TITLE C** LOCAL VARIABLES *** C** HEADER CARD INTEGER IT ! TYPE OF PACKET THAT WILL FOLLOW INTEGER ID ! IDENTIFICATION NUMBER INTEGER IV ! ADDITIONAL ID INTEGER KC ! CARK COUNT (NUMBER OF DATA CARDS IN THIS PACKET INTEGER N(5) ! SUPPLEMENTAL INTEGER VALUES C** NODE DATA INTEGER ICF ! CONDENSATION FLAG (0=UNREFERENCED) CHARACTER TYPE ! NODE TYPE INTEGER NDF ! NUMBER OF DEGREES OF FREEDOM INTEGER CONFIG ! NODE CONFIGURATION INTEGER CID ! COORDINATE FRAME DO ANALYSIS RESULTS INTEGER M(6) ! 6 PERMANENT SINGLE POINT CONSTRAINT FLAGS INTEGER NODE ! NODE NUMBER DOUBLE PRECISION XNODE,YNODE ! NODAL DATA FOR SIDE NODES C** ELEMENT DATA INTEGER NODES ! TOTAL NUMBER OF NODES DO THIS ELEMENT C CONFIG ! ELEMENT CONFIGURATION INTEGER PID ! PROP ID (+) OR MAT ID (-) INTEGER CEID ! CONGRUENT ELEMENT ID REAL MATOA(3) ! MATERIAL ORIENTATION ANGLES DOUBLE PRECISION DET INTEGER NUM_ND, NUM_EL, NUM_MAT, NUM_PROP, I, J, K, IDEBUG INTEGER MAX_NODE_NUMBER,MAX_ELEMENT_NUMBER INTEGER TWIDTH, IIN, IOUT, IOUT2, II, JJ REAL RADIUS CHARACTER DATE*12, TIME*8, VERSION*9, DUM LOGICAL AT_END COMMON IIN, IOUT, IDEBUG C C Initialize some stuff C AT_END = .FALSE. NUMNP = 0 NUMEL = 0 NUMMAT = 0 MAT(1) = 0 TWIDTH = 11 MAX_NODE_NUMBER = 0 MAX_ELEMENT_NUMBER = 0 C C begin main loop C 5000 CONTINUE C C Read the header card C READ(20,'(I2,8I8)',END=100) IT, ID, IV, KC,(N(I),I=1,5) C C Node data C IF (IT.EQ.1) THEN READ(20,'(3E16.9)',END=200) X1,Y1,Z1 READ(20,'(I1,A1,3I8,2X,6I1)',END=200) ICF,TYPE,NDF,CONFIG, _ CID,(M(I),I=1,6) IF (N(1).NE.0) THEN WRITE(IOUT,*)'ERROR - ' STOP C ELSE IF (CID.NE.0) THEN C WRITE(IOUT,*)'ERROR - BAD COORDINATE FRAME' C STOP C ELSE IF (NDF.NE.1) C WRITE(IOUT,*)'ERROR - TOO MANY DOF AT ELM NODES' C WRITE(IOUT,*)' NODE =',ID C STOP END IF NUMNP = NUMNP + 1 X(ID)=X1 Y(ID)=Y1 IF (ID.GT.MAX_NODE_NUMBER) MAX_NODE_NUMBER = ID C C Element data C ELSE IF (IT.EQ.2) THEN READ (20,'(4I8,3E16.9)',END=200)NODES,CONFIG,PID,CEID, _ (MATOA(I),I=1,3) C C Twidth is 11 in FRANC ==> nodes=6, or nodes=8 C IF ((NODES.NE.6) .AND. (NODES.NE.8)) THEN WRITE(IOUT,*)'ERROR - NUMBER OF NODES IS NOT CORRECT' WRITE(IOUT,*)' ELEM=',ID,' NUMBER OF NODES=',NODES STOP END IF C C Force configuration to straight sided. C CONFIG = 1 C READ(20,'(10I8)',END=200)(NTABLE(I,ID),I=1,NODES) IF (PID.EQ.0) PID = 1 NUMEL = NUMEL + 1 NTABLE(9,ID) = NUMEL NTABLE(10,ID) = IABS(PID) if ( nodes .eq. 6 ) NTABLE(7,ID) = 0 IF (ID.GT.MAX_ELEMENT_NUMBER) MAX_ELEMENT_NUMBER= ID C C Check for a new unique material table value C DO I=1,NUMMAT IF (IABS(PID).EQ.MAT(I)) GO TO 300 END DO IF (NUMMAT.LT.MAMAX) THEN NUMMAT = NUMMAT + 1 MAT(NUMMAT) = IABS(PID) ELSE WRITE(IOUT,*)'ERROR - TOO MANY UNIQUE MATERIALS' WRITE(IOUT,*)' NUMBER OF MAT''S ALLOWED =',MAMAX WRITE(IOUT,*)' MAKE THE PARAMETER MAMAX LARGER' END IF 300 CONTINUE C C Check for associate data C IF (N(1).NE.0) THEN WRITE(IOUT,*)'ERROR - ASSOCIATE DATA SHOULD NOT BE ' _ //' PRESENT FOR ELEMENT READ, ELM=',ID STOP END IF C C Title card C ELSE IF (IT.EQ.25) THEN READ(20,'(A)',END=200) TITLE C C Summary card C ELSE IF (IT.EQ.26) THEN READ(20,'(A12,A8,A9)',END=200) DATE, TIME, VERSION NUM_ND = N(1) NUM_EL = N(2) NUM_MAT = N(3) NUM_PROP = N(4) IF (N(1).GT.NDMAX) THEN WRITE(IOUT,*)'ERROR - NODES TO BE READ IN IS GREATER THAN' _ //' MAX ALLOWED' STOP ELSE IF (N(2).GT.ELMAX) THEN WRITE(IOUT,*)'ERROR - ELM TO BE READ IN IS GREATER THAN' _ //' MAX ALLOWED' STOP END IF C C End of file C ELSE IF(IT.EQ.99) THEN AT_END = .TRUE. C C Data is present that we don't translate C ELSE WRITE(IOUT,*)'Neutral file contains information not used', 1 ' in FRANC' WRITE(IOUT,*)'Skipping the data' WRITE(IOUT,*)' SET TYPE=',IT WRITE(IOUT,*)' LOOK TO PATRAN NATURAL FILE DEF FOR TYPE' DO I=1, KC READ(20,'(A)',END=200) DUM END DO END IF IF ( .NOT. AT_END ) GOTO 5000 C C Check data some C IF (NUMNP.NE.NUM_ND) THEN WRITE(IOUT,*)'ERROR - NUMBER OF NODES READ IN DOES NOT EQUAL', _ ' NUMBER OF NODES IN MODEL' WRITE(IOUT,*)' NDS=',NUMNP,' NUM_NODES=',NUM_ND STOP ELSE IF (NUMEL.NE.NUM_EL) THEN WRITE(IOUT,*)'ERROR - NUMBER OF ELEM READ IN DOES NOT EQUAL', _ ' NUMBER OF ELEM IN MODEL' WRITE(IOUT,*)' ELM=',NUMEL,' NUM_EL=',NUM_EL STOP END IF C C Done reading data. C 100 CONTINUE C C Check if the element and node numbers are compressed. C IF (MAX_NODE_NUMBER.GT.NUMNP.OR.MAX_ELEMENT_NUMBER.GT.NUMEL) THEN WRITE(IOUT2,*)'ERROR Nodes/elements not numbered sequentially.' WRITE(IOUT2,*)' Use COMPRESS command in PATRAN.' STOP END IF C C Make sure all midside nodes cause straight sided elements. C DO I = 1,NUMEL X1 = X(NTABLE(1,I)) X2 = X(NTABLE(2,I)) X3 = X(NTABLE(3,I)) X4 = X(NTABLE(4,I)) Y1 = Y(NTABLE(1,I)) Y2 = Y(NTABLE(2,I)) Y3 = Y(NTABLE(3,I)) Y4 = Y(NTABLE(4,I)) if ( NTABLE(7,I) .EQ. 0 ) then DO JJ = 4,6 IF ((NTABLE(JJ,I).NE.NTABLE(1,I)).AND. 1 (NTABLE(JJ,I).NE.NTABLE(2,I)).AND. 2 (NTABLE(JJ,I).NE.NTABLE(3,I))) THEN IF (JJ.EQ.4) THEN X(NTABLE(4,I)) = (X1+X2)*0.5 Y(NTABLE(4,I)) = (Y1+Y2)*0.5 ELSE IF (JJ.EQ.5) THEN X(NTABLE(5,I)) = (X2+X3)*0.5 Y(NTABLE(5,I)) = (Y2+Y3)*0.5 ELSE X(NTABLE(6,I)) = (X3+X1)*0.5 Y(NTABLE(6,I)) = (Y3+Y1)*0.5 END IF END IF END DO ! DO JJ = 4,6 ELSE DO JJ = 5,8 IF ((NTABLE(JJ,I).NE.NTABLE(1,I)).AND. 1 (NTABLE(JJ,I).NE.NTABLE(2,I)).AND. 2 (NTABLE(JJ,I).NE.NTABLE(3,I)).AND. 3 (NTABLE(JJ,I).NE.NTABLE(4,I))) THEN IF (JJ.EQ.5) THEN X(NTABLE(5,I)) = (X1+X2)*0.5 Y(NTABLE(5,I)) = (Y1+Y2)*0.5 ELSE IF (JJ.EQ.6) THEN X(NTABLE(6,I)) = (X2+X3)*0.5 Y(NTABLE(6,I)) = (Y2+Y3)*0.5 ELSE IF (JJ.EQ.7) THEN X(NTABLE(7,I)) = (X3+X4)*0.5 Y(NTABLE(7,I)) = (Y3+Y4)*0.5 ELSE X(NTABLE(8,I)) = (X4+X1)*0.5 Y(NTABLE(8,I)) = (Y4+Y1)*0.5 END IF END IF END DO ! DO JJ = 4,6 END IF C C Make sure that elem are numbered counter clockwise. C c DET = (X2*Y3-X3*Y2) - (X1*Y3-X3*Y1) + (X1*Y2-X2*Y1) c IF (DET.LT.0.0) THEN c J = NTABLE(2,I) c NTABLE(2,I) = NTABLE(3,I) c NTABLE(3,I) = J c J = NTABLE(4,I) c NTABLE(4,I) = NTABLE(6,I) c NTABLE(6,I) = J c END IF END DO C C Make sure that the side nodes are in the middle of corner nodes. C If the side nodes are the corner nodes (PATRAN has this BUG), then C generate these nodes. C RADIUS = SQRT ( (X(NTABLE(1,1))-X(NTABLE(2,1)))**2 + 1 (Y(NTABLE(1,1))-Y(NTABLE(2,1)))**2 ) * 1.E-6 DO II = 1,NUMEL if ( NTABLE(7,I) .EQ. 0 ) then DO JJ = 4,6 IF ((NTABLE(JJ,II).EQ.NTABLE(1,II)).OR. 1 (NTABLE(JJ,II).EQ.NTABLE(2,II)).OR. 2 (NTABLE(JJ,II).EQ.NTABLE(3,II))) THEN IF (JJ.NE.6) THEN XNODE = (X(NTABLE(JJ-2,II))+X(NTABLE(JJ-3,II)))*0.5 YNODE = (Y(NTABLE(JJ-2,II))+Y(NTABLE(JJ-3,II)))*0.5 ELSE XNODE = (X(NTABLE(1,II))+X(NTABLE(3,II)))*0.5 YNODE = (Y(NTABLE(1,II))+Y(NTABLE(3,II)))*0.5 END IF CALL HIT_NODE_ABS_RAD(X,Y,NUMNP,XNODE, 1 YNODE,NODE,RADIUS) IF ( NODE.EQ.0 ) THEN NUMNP = NUMNP+1 IF ( NUMNP.LT.NDMAX ) THEN NTABLE(JJ,II) = NUMNP X(NUMNP) = XNODE Y(NUMNP) = YNODE ELSE WRITE(IOUT,*)' ERROR - NODE TO BE GENERATED IS ', 1 'GREATER THAN MAX ALLOWED' STOP END IF ELSE NTABLE(JJ,II) = NODE END IF END IF END DO ! of DO JJ = 4,6 ELSE DO JJ = 5,8 IF ((NTABLE(JJ,II).EQ.NTABLE(1,II)).OR. 1 (NTABLE(JJ,II).EQ.NTABLE(2,II)).OR. 2 (NTABLE(JJ,II).EQ.NTABLE(3,II)).OR. 2 (NTABLE(JJ,II).EQ.NTABLE(4,II))) THEN IF (JJ.NE.8) THEN XNODE = (X(NTABLE(JJ-3,II))+X(NTABLE(JJ-4,II)))*0.5 YNODE = (Y(NTABLE(JJ-3,II))+Y(NTABLE(JJ-4,II)))*0.5 ELSE XNODE = (X(NTABLE(1,II))+X(NTABLE(4,II)))*0.5 YNODE = (Y(NTABLE(1,II))+Y(NTABLE(4,II)))*0.5 END IF CALL HIT_NODE_ABS_RAD(X,Y, 1 NUMNP,XNODE,YNODE,NODE,RADIUS) IF ( NODE.EQ.0 ) THEN NUMNP = NUMNP+1 IF ( NUMNP.LT.NDMAX ) THEN NTABLE(JJ,II) = NUMNP X(NUMNP) = XNODE Y(NUMNP) = YNODE ELSE WRITE(IOUT,*)' ERROR - NODE TO BE GENERATED IS ', 1 'GREATER THAN MAX ALLOWED' STOP END IF ELSE NTABLE(JJ,II) = NODE END IF END IF END DO ! of DO JJ = 4,6 END IF END DO ! of DO II = 1,NUMEL C C If compiled with the debug switch on this information will be written C to the default terminal number defined in the calling routine. C IF ( IDEBUG .NE. 0 ) THEN c IOUT2 = 22 ! Iout2 is hardwired in the main program. IOUT2 = 6 WRITE(IOUT2,1000) NDMAX, ELMAX WRITE(IOUT2,*) 'NUMBER OF MATERIALS = ',NUMMAT WRITE(IOUT2,*) TITLE WRITE(IOUT2,*)'NODE LOCATIONS' DO j=1, NUMNP WRITE(IOUT2,'(I8,2X,2(F12.5,2X))')j,X(j),Y(j) END DO WRITE(IOUT2,*)'NODE TABLE' DO k=1,NUMEL WRITE(IOUT2,'(I5,8X,20I8)') k,(NTABLE(J,k),J=1,TWIDTH) END DO END IF RETURN 1000 FORMAT(' MAX NODES ALLOWED=',I5,3X,'MAX ELM ALLOWED=',I5) 200 WRITE(IOUT2,*)'ERROR DURING READ, PROGRAM TERMINATED' STOP END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C SUBROUTINE HIT_NODE_ABS_RAD(X,Y,NUM,XNODE,YNODE,NODE,RADI) patcrk C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C This subroutine determines if the user has hit a node in the mesh. It C C returns the node number if a hit has occured. C C C C Written by: Daniel Swenson - August 1984 C C Modified by: Mark A. James - February 1988. C C Fixed to find closest node to point within the radius. C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C c IMPLICIT NONE C DOUBLE PRECISION X(*),Y(*) ! Coordinates. DOUBLE PRECISION XNODE,YNODE ! Node coordinates. DOUBLE PRECISION RAD ! Radius. REAL RADI C INTEGER NUM ! Number of search nodes. INTEGER NODE ! Hit node. C DOUBLE PRECISION RADIUS,CHKRAD ! Check for proximity. DOUBLE PRECISION MIN ! Look for closest. C INTEGER II ! Counter. C C For this entry point, an absolute radius is give for the search. C RAD = RADI RADIUS = RAD**2 C C Check each node to see if the given point lies within a radius around C a node. Note that it finds the closest node to the point. patcrk  C 5 NODE = 0 MIN = RADIUS DO 10 II = 1,NUM IF(X(II).LE.1.E10) THEN CHKRAD = (X(II)-XNODE)**2+(Y(II)-YNODE)**2 pat IF(CHKRAD.LT.RADIUS) THEN IF ( CHKRAD.LT.MIN ) THEN NODE = II MIN = CHKRAD END IF pa END IF END IF 10 CONTINUE C 20 RETURN END