C	FORTRAN PROGRAM TO GENERATE A FINITE ELEMENT MESH
C	FOR 2-DIMENSIONAL PROBLEMS
C
C
	DIMENSION IDBLK(400),NSD(20),NWD(20)
        DIMENSION NNAR(1000)
	DIMENSION XB(450,3),SR(450,3),WR(450,3),MERG(30,4)
	DIMENSION X(450,3),XP(8,3),NOC(1000,4),MAT(1000),SH(8)
        COMMON/AA/NGCN(450)
	NDIM=2
	WRITE(*,1001)
1001	FORMAT(1X,'CHOOSE ELEMENT TYPE <1=TRIANGLE,
     +  2=QUADRILATERAL>:',$)
	READ(*,*)NTMP
	IF(NTMP.EQ.1) THEN
	 NEN=3
	 ELSE 
	  NEN=4
	ENDIF
	WRITE(*,1002)
1002	FORMAT(1X,' NUMBER OF S-SPANS=?',$)
	READ(*,*)NS
	WRITE(*,1003)
1003	FORMAT(1X,' NUMBER OF W-SPANS=?',$)
	READ(*,*)NW
	NSW=NS*NW
	NGN=(NS+1)*(NW+1)
	NM=1
C	BLOCK IDENTIFIERS FOR BLOCKS
C
	DO 10 I=1,NSW
	 IDBLK(I)=1
10	CONTINUE
	WRITE(*,1004)
1004	FORMAT(1X,' NO. OF VOID BLOCKS OR BLOCKS WITH
     +  MATERIAL # NOT 1 ?'$) 
	READ(*,*)NON1BLK
        IF(NON1BLK.NE.0) THEN
	WRITE(*,*)'GIVE THE BLOCK AND MATERIAL # FOR EACH 
     + VOID BLOCK OR '
	WRITE(*,*) ' BLOCK WITH MATERIAL # NOT=1' 
	DO 20 I=1,NON1BLK
15	 READ(*,*)J,IDBLK(J)
	 IF(NM.LT.IDBLK(J)) NM=IDBLK(J)
20	CONTINUE
	ENDIF

C	SPAN DIVISIONS
C
	NNS=1
	NNW=1
	WRITE(*,*)' INPUT THE NO. OF DIVISIONS IN EACH S-SPAN' 
	DO 30 KS=1,NS
	 WRITE(*,1005)KS
1005	 FORMAT(1X,'S-SPAN # =',I3,'   NO. OF DIVISIONS=?',$)
25	 READ(*,*) NSD(KS)
	 NNS=NNS+NSD(KS)
30	CONTINUE

	WRITE(*,*)' INPUT THE NO. OF DIVISIONS IN EACH W-SPAN'
	DO 40 KW=1,NW
	 WRITE(*,1006) KW
1006	 FORMAT(1X,'W-SPAN # =',I3,'   NO.OF DIVISIONS =?',$)
35	 READ(*,*)NWD(KW)
	 NNW=NNW+NWD(KW)
40	CONTINUE
C
C	BLOCK CORNER AND SIDE DATA
C
	NSR=NS*(NW+1)
	NWR=NW*(NS+1)
	WRITE(*,*)' INPUT THE NO. OF CORNER NODES'
	READ(*,*)NOCN
C.......Z-COORDINATES ARE SET TO ZERO , BUT MAY BE READ IN HERE !
	WRITE(*,*)' GIVE CORNER NODE NO. AND X , Y COORDINATES
     + FOR EACH NODE'
	DO 50 I=1,NOCN
	 READ(*,*) N,XB(I,1),XB(I,2)
         XB(I,3)=0.0
50	CONTINUE
C
C	EVALUATE MID-POINTS OF S-SIDES
C
	DO 60 I=1,NW+1
	 DO 60 J=1,NS
	  IJ=(I-1)*NS+J
	  SR(IJ,1)=.5*(XB(IJ+I-1,1)+XB(IJ+I,1))
	  SR(IJ,2)=.5*(XB(IJ+I-1,2)+XB(IJ+I,2))
	  SR(IJ,3)=.5*(XB(IJ+I-1,3)+XB(IJ+I,3))
60	CONTINUE
C
C	EVALUATE MID-POINTS OF W-SIDES
C
	DO 70 I=1,NW
	 DO 70 J=1,NS+1
	  IJ=(I-1)*(NS+1)+J
	  WR(IJ,1)=.5*(XB(IJ,1)+XB(IJ+NS+1,1))
	  WR(IJ,2)=.5*(XB(IJ,2)+XB(IJ+NS+1,2))
	  WR(IJ,3)=.5*(XB(IJ,3)+XB(IJ+NS+1,3))
70	CONTINUE
C
C	READ CO-ORDINATES OF CURVED SIDES OR SIDES WITH
C	DISPLACED MID-NODES
C
	WRITE(*,1015)
1015	FORMAT(1X,' INPUT THE NUMBER OF S-SIDES WHICH ARE
     + CURVED OR THOSE WHICH ARE STRAIGHT') 
	WRITE(*,1025)
1025	FORMAT(1X,'BUT HAVE DISPLACED MID-POINTS :',$)
	READ(*,*)NSDP
	IF(NSDP.EQ.0)GOTO 81
	WRITE(*,*)'GIVE THE S-SIDE# AND THE X,Y
     + COORDINATES FOR EACH SIDE'
	DO 80 I=1,NSDP
75	 READ(*,*)NSS,SR(NSS,1),SR(NSS,2)
         SR(NSS,3) = 0.0
80	CONTINUE
81	WRITE(*,1007)
1007	FORMAT(1X,' INPUT THE NUMBER OF W-SIDES WHICH ARE
     + CURVED OR WITH DISPLACED MID-POINTS :',$)
	READ(*,*)NWDP
	IF(NWDP.EQ.0)GOTO 91
	WRITE(*,*)' GIVE THE W-SIDE# & THE X,Y COORDINATES FOR
     + EACH W-SIDE'
	DO 90 I=1,NWDP
85	 READ(*,*)NWW,WR(NWW,1),WR(NWW,2)
         WR(NWW,3) = 0.0
90	CONTINUE
C
C	MERGING SIDES
C
 91	WRITE(*,1008)
1008	FORMAT(1X,'NUMBER OF PAIRS OF SIDES TO BE MERGED=?',$)
	READ(*,*)NSJ
	IF(NSJ.EQ.0)GOTO 107
	DO 100 I=1,NSJ
	 WRITE(*,1009)
1009	 FORMAT(1X,'END NODES OF SIDE 1 < SIDE 1 IS THE ONE
     + WITH LOWER CORNER NODE NUMBERS>=?',$)
95	 READ(*,*)L1,L2
	 CALL TEST(NNS,L1,L2,II1)
	 WRITE(*,1010)
1010	 FORMAT(1X,' INPUT END NODES OF MERGING SIDE',$)
	 READ(*,*)L3,L4
	 CALL TEST(NNS,L3,L4,II2)
	 IF(II1.NE.II2) THEN
	  WRITE(*,*)'DIVISIONS ON MERGING SIDES DONT MATCH'
	  WRITE(*,*)'ERROR IN INPUT DATA;PLEASE CHECK'
	  STOP
	 ENDIF
	 MERG(I,1)=L1
	 MERG(I,2)=L2
	 MERG(I,3)=L3
	 MERG(I,4)=L4
100	CONTINUE
C
C	GLOBAL NODE LOCATIONS OF CORNER NODES
C
 107	NTMPI=1
	DO 110 I=1,NW+1
	 IF(I.EQ.1) THEN
	  IINC=0
	 ELSE
	  IINC=NNS*NWD(I-1)
	 ENDIF
	 NTMPI=NTMPI+IINC
	 NTMPJ=0
	 DO 110 J=1,NS+1
	  IJ=(NS+1)*(I-1)+J
	  IF(J.EQ.1) THEN
	     JINC=0
	   ELSE
	     JINC=NSD(J-1)
	  ENDIF
	  NTMPJ=NTMPJ+JINC
	  NGCN(IJ)=NTMPI+NTMPJ
110	CONTINUE
C
C	NODE POINT ARRAY
C
	NNT=NNS*NNW
	DO 120 I=1,NNT
	 NNAR(I)=-1
120	CONTINUE
C
C	ZERO NON EXISTING NODE LOCATIONS
C
	DO 130 KW=1,NW
	 DO 130 KS=1,NS
	  KSW=NS*(KW-1)+KS
	  IF(IDBLK(KSW).GT.0) GOTO 130
C	
C	OPERATIONS WITHIN AN EMPTY BLOCK
C
          K1=(KW-1)*(NS+1)+KS
	  N1=NGCN(K1)
	  NS1=2
	  IF(KS.EQ.1) NS1=1
	  NW1=2
	  IF(KW.EQ.1) NW1=1
	  NS2=NSD(KS)+1
 	  IF(KS.EQ.NS) GOTO 123
 	  IF(IDBLK(KSW+1).GT.0) NS2=NSD(KS)
123	  NW2=NWD(KW)+1
	  IF(KW.EQ.NW) GOTO 126
	  IF(IDBLK(KSW+NS).GT.0) NW2=NWD(KW)
126	  DO 128 I=NW1,NW2
	    IN1=N1+(I-1)*NNS
	    DO 128 J=NS1,NS2
	     IJ=IN1+J-1
	     NNAR(IJ)=0
128       CONTINUE
	  IF((NS2.EQ.NSD(KS)).OR.(NW2.EQ.NWD(KW))) GOTO 130
	  IF((KS.EQ.NS).OR.(KW.EQ.NW)) GOTO 130
	  IF(IDBLK(KSW+NS+1).GT.0) NNAR(IJ)=-1
130	CONTINUE
C
C	NODE IDENTIFICATION FOR SIDE MERGING
C
	IF(NSJ.EQ.0) GOTO 141
	DO 140 I=1,NSJ
	 I1=MERG(I,1)
	 I2=MERG(I,2)
	 CALL TEST(NNS,I1,I2,IDIV)
	 IA1=NGCN(I1)
	 IA2=NGCN(I2)
	 IASTP=(IA2-IA1)/IDIV
	 I1=MERG(I,3)
	 I2=MERG(I,4)
	 CALL TEST(NNS,I1,I2,IDIV)
	 IB1=NGCN(I1)
	 IB2=NGCN(I2)
	 IBSTP=(IB2-IB1)/IDIV
	 IAA=IA1-IASTP
	  DO 140 IBB=IB1,IB2,IBSTP
	   IAA=IAA+IASTP
	   IF(IBB.EQ.IAA)THEN 
	       NNAR(IAA)=-1
	      ELSE 
	       NNAR(IBB)=IAA
	   ENDIF
140	CONTINUE
C
C	FINAL NODE NUMBERS IN THE ARRAY
C
 141	NODE=0
	DO 150 I=1,NNT
	 IF(NNAR(I).EQ.0) GOTO 150
	 IF(NNAR(I).GT.0) GOTO 145
	 NODE=NODE+1
	 NNAR(I)=NODE
	 GOTO 150
145 	 II=NNAR(I)
	 NNAR(I)=NNAR(II)
150	CONTINUE
C
C	NODAL COORDINATES
C
	NN=NODE
	NELM=0
	DO 160 KW=1,NW
	 DO 160 KS=1,NS
	  KSW=NS*(KW-1)+KS
	  IF(IDBLK(KSW).EQ.0) GOTO 159
C
C	  EXTRACTION OF BLOCK DATA
C
	  NODW=NGCN(KSW+KW-1)-NNS-1
	  DO 160 JW=1,(NWD(KW)+1)
	   ETA=-1.+(2.*FLOAT(JW-1)/FLOAT(NWD(KW)))
	   NODW=NODW+NNS
	   NODS=NODW
	   DO 160 JS=1,(NSD(KS)+1)
	    XI=-1.+2.*(FLOAT(JS-1)/FLOAT(NSD(KS)))
	    NODS=NODS+1
	    NODE=NNAR(NODS)
	    CALL COORD(NS,KSW,KW,XB,SR,WR,XP)
	    CALL SHAPE(SH,XI,ETA)
	    DO 155 J=1,3
	     C1=0.
	     DO 153 I=1,8
	      C1=C1+SH(I)*XP(I,J)
153	     CONTINUE
	     X(NODE,J)=C1
155	    CONTINUE
C
C	    CONNECTIVITY
C
	    IF((JS.EQ.(NSD(KS)+1)).OR.(JW.EQ.(NWD(KW)+1)))
     +       GOTO 159
	    N1=NODE
	    N2=NNAR(NODS+1)
	    N3=NNAR(NODS+NNS+1)
	    N4=NNAR(NODS+NNS)
	    NELM=NELM+1
	    IF(NEN.EQ.3) GOTO 2600
C
C	   QUADRILATERAL ELEMENTS
C
	   NOC(NELM,1)=N1
	   NOC(NELM,2)=N2
	   NOC(NELM,3)=N3
	   NOC(NELM,4)=N4
	   MAT(NELM)=IDBLK(KSW)
	   GOTO 159

C
C	   TRIANGULAR ELEMENTS
C
2600	   NOC(NELM,1)=N1
	   NOC(NELM,2)=N2
	   NOC(NELM,3)=N3
	   MAT(NELM)=IDBLK(KSW)
	   NELM=NELM+1
	   NOC(NELM,1)=N3
	   NOC(NELM,2)=N4
	   NOC(NELM,3)=N1
	   MAT(NELM)=IDBLK(KSW)
 159     CONTINUE
160	CONTINUE
	IF(NEN.EQ.4) GOTO 170
	NE=NELM
C
C	READJUSMENT FOR TRIANGULAR CONNECTIVITY
C

	NE2=NE/2
	DO 165 I=1,NE2
	I1=2*I-1
	N1=NOC(I1,1)
	N2=NOC(I1,2)
	N3=NOC(I1,3)
	N4=NOC(2*I,2)
	X13=X(N1,1)-X(N3,1)
	Y13=X(N1,2)-X(N3,2)
	Z13=X(N1,3)-X(N3,3)
	X24=X(N2,1)-X(N4,1)
	Y24=X(N2,2)-X(N4,2)
	Z24=X(N2,3)-X(N4,3)
	XFACT=(X13**2+Y13**2+Z13**2)/
     +        (X24**2+Y24**2+Z24**2)
	XFACT=SQRT(XFACT)
	IF(XFACT.LE.0.9) GOTO 165
	NOC(I1,3)=N4
	NOC(2*I,3)=N2
165	CONTINUE	
170	CALL OUTPUT(NN,NELM,X,NOC,MAT,NEN,NDIM,NM)
	CLOSE(5)
	END


C	SUBROUTINE TO CALCULATE THE NUMBER OF DIVISIONS
C	FOR SIDES I1 AND I2
C
	SUBROUTINE TEST(NNS,I1,I2,IDIV)
	COMMON/AA/ NGCN(450)
	IMIN=I1
	IMAX=I2
	IF(I1.EQ.I2) WRITE(*,*)'*****ERROR;I1=I2'
	IF(IMIN.GT.I2) THEN
	 IMIN=I2
	 IMAX=I1
	ENDIF
	IF((IMAX-IMIN).EQ.1) THEN
	  IDIV=NGCN(IMAX)-NGCN(IMIN)
	  RETURN
	ENDIF
	IDIV=(NGCN(IMAX)-NGCN(IMIN))/NNS
	RETURN
	END


C
C	SUBROUTINE TO EXTRACT THE COORDINATES OF EACH OF
C	THE EIGHT NODES OF THE BLOCK FROM THE BULK DATA
C
	SUBROUTINE COORD(NS,KSW,KW,XB,SR,WR,XP)
	DIMENSION XB(450,3),SR(450,3),WR(450,3),XP(8,3)
	N1=KSW+KW-1
	XP(1,1)=XB(N1,1)
	XP(1,2)=XB(N1,2)
	XP(1,3)=XB(N1,3)
	XP(3,1)=XB(N1+1,1)
	XP(3,2)=XB(N1+1,2)
	XP(3,3)=XB(N1+1,3)
	XP(5,1)=XB(N1+NS+2,1)
	XP(5,2)=XB(N1+NS+2,2)
	XP(5,3)=XB(N1+NS+2,3)
	XP(7,1)=XB(N1+NS+1,1)
	XP(7,2)=XB(N1+NS+1,2)
	XP(7,3)=XB(N1+NS+1,3)
	XP(2,1)=SR(KSW,1)
	XP(2,2)=SR(KSW,2)
	XP(2,3)=SR(KSW,3)
	XP(6,1)=SR(KSW+NS,1)
	XP(6,2)=SR(KSW+NS,2)
	XP(6,3)=SR(KSW+NS,3)
	XP(8,1)=WR(N1,1)
	XP(8,2)=WR(N1,2)
	XP(8,3)=WR(N1,3)
	XP(4,1)=WR(N1+1,1)
	XP(4,2)=WR(N1+1,2)
	XP(4,3)=WR(N1+1,3)
	RETURN
	END


C
C	SUBROUTINE TO EVALUATE THE SHAPE FUNCTIONS
C
	SUBROUTINE SHAPE(SH,XI,ETA)
	DIMENSION SH(8)
	SH(1)=-(1.-XI)*(1.-ETA)*(1.+XI+ETA)/4.
	SH(2)=(1.-XI*XI)*(1.-ETA)/2.
	SH(3)=-(1.+XI)*(1.-ETA)*(1.-XI+ETA)/4.
	SH(4)=(1.-ETA*ETA)*(1.+XI)/2.
	SH(5)=-(1.+XI)*(1.+ETA)*(1.-XI-ETA)/4.
	SH(6)=(1.-XI*XI)*(1.+ETA)/2.
	SH(7)=-(1.-XI)*(1.+ETA)*(1.+XI-ETA)/4.
	SH(8)=(1.-ETA*ETA)*(1.-XI)/2.
	RETURN
	END

C
C	SUBROUTINE TO SAVE GENERATED MESH DATA
C
	SUBROUTINE OUTPUT(NN,NELM,X,NOC,MAT,NEN,NDIM,NM)
	DIMENSION X(450,3),NOC(1000,4),MAT(1000)
        CHARACTER*12 OUTFIL
        LOUT = 11
        write(6,107)
 107    FORMAT(1X,'# CONSTRAINED DOF"s, # OF LOAD COMPONENTS = ',$)
        READ(5,*) ND,NL
        WRITE(6,109)
 109    FORMAT(1X,'FILE NAME FOR OUTPUT : ',$)
        READ(5,108)OUTFIL
 108    FORMAT(A)
	OPEN(11,FILE=OUTFIL,STATUS='UNKNOWN')
	WRITE(LOUT,172) NN,NELM,ND,NL,NM,NDIM,NEN
172	FORMAT(16I5)
	DO 200 I=1,NN
	 WRITE(LOUT,199)(X(I,J),J=1,NDIM)
199	 FORMAT(6(1X,E12.5))
200	CONTINUE
	DO 210 I=1,NELM
	 WRITE(LOUT,172)(NOC(I,J),J=1,NEN)
210	CONTINUE
         WRITE(LOUT,172) (MAT(I), I=1,NELM)
        CLOSE(UNIT= 11)
	RETURN
	END