You cannot select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1461 lines
56 KiB
Fortran

!IPK LAST UPDATE SEP 23 2015 ADD TESTING FOR CHNAGED ELEMENTS/NODES
! last update Sept 20 1999
! Last change: IPK 13 Jan 98 10:05 am
!ipk last update Nov 18 1997
!ipk last update Oct 24 1996
SUBROUTINE REFB
!
! Routines to control refinement of elements
!
USE BLK1MOD
INCLUDE 'BFILES.I90'
! INCLUDE 'BLK1.COM'
!
CHARACTER*1 ANS,ANSW(10)
DATA ANSW/'f','l','s','t','v','n',' ','m',' ','q'/
!
! Draw box around selections
!
100 CONTINUE
NHTP=8
NMESS=0
NBRR=0
CALL HEDR
!
! Get answer
!
!ipk jan98
210 continue
call wrtbox(idelv)
call xyloc(XPT,YPT,ANS,IBOX)
IF(IRMAIN .EQ. 1) RETURN
!ipk jan98 add option for deleting elevation on move
IF(IBOX .EQ. 7 .or. ANS .eq. 'e') THEN
IDELV=MOD(IDELV+1,2)
GO TO 210
ENDIF
IF(ANS .EQ. 'c') THEN
if(ibox .eq. 0) go to 210
ANS=ANSW(IBOX)
ENDIF
!
! Element generation
!
IF (ANS .EQ. 'f') THEN
!
! Refine elements by four
!
IECHG=0
!IPK MAY03
ICHG=0
CALL REFIN(0)
IRDONE=0
IF(IRMAIN .EQ. 1) RETURN
!
ELSEIF (ANS .EQ. 'l') THEN
!
! Refine elements by two long
!
IECHG=0
!IPK MAY03
ICHG=0
CALL REFIN(1)
IRDONE=0
IF(IRMAIN .EQ. 1) RETURN
!
ELSEIF (ANS .EQ. 's') THEN
!
! Refine elements by two short
!
IECHG=0
!IPK MAY03
ICHG=0
CALL REFIN(2)
IRDONE=0
IF(IRMAIN .EQ. 1) RETURN
!
!
ELSEIF (ANS .EQ. 't') THEN
!
! Refine elements by splitting quads
!
IECHG=0
!IPK MAY03
ICHG=0
CALL REFIN(3)
IRDONE=0
IF(IRMAIN .EQ. 1) RETURN
!
!
ELSEIF (ANS .EQ. 'v') THEN
!
! Reverse element diagonals for quads
!
IECHG=0
!IPK MAY03
ICHG=0
CALL REFIN(4)
IRDONE=0
IF(IRMAIN .EQ. 1) RETURN
!
ELSEIF (ANS .EQ. 'n') THEN
!
! Clean up element refinement
!
IECHG=0
!IPK MAY03
ICHG=0
CALL CLENUP(0)
IRDONE=0
IF(IRMAIN .EQ. 1) RETURN
!
ELSEIF (ANS .EQ. 'm') THEN
IF(IRMAIN .EQ. 1) RETURN
!
! simplify layout
!
IECHG=0
!IPK MAY03
ICHG=0
CALL SMFY !
IRDONE=0
ELSEIF (ANS .EQ. 'q') THEN
CALL WRTOUT(0) !
RETURN
!
! Look again
!
ENDIF
GO TO 100
END
!
SUBROUTINE REFIN(ITYPT)
!
! Routine to refine elements
!
USE BLK1MOD
! INCLUDE 'BLK1.COM'
INCLUDE 'TXFRM.COM'
!IPK MAY02 COMMON /TXFRM/ XS, YS, TXSCAL
!
DIMENSION NTRAN(9),IELGB(8)
CHARACTER*1 IFLAG
DIST(N1,N2)=SQRT((CORD(N1,1)-CORD(N2,1))**2 &
& +(CORD(N1,2)-CORD(N2,2))**2)
!
ITYP=ITYPT
IF(NEFL .GT. 0) GO TO 150
!ipk may94 change so that refine does not change display
! DO 2 I=1,9
! IPSW(I)=0
! 2 CONTINUE
! IPSW(4)=1
! CALL PLOTOT
!ipk may94 end changes
3 CONTINUE
NHTP=0
NMESS=12
NBRR=3
CALL HEDR
!
! Write out
!
NEFL=0
4 CONTINUE
IBOX=1
CALL PROX(XC,YC,NE,XX,YY,IELEM,IFLAG,IESKP,IBOX)
IF(IRMAIN .EQ. 1) RETURN
!
IF (IFLAG .EQ. 'c') THEN
NEFL=NEFL+1
NEFLAG(NEFL)=IELEM
CALL FILLEM(IELEM)
ELSEIF (IFLAG .EQ. 'U') THEN
NEFLAG(NEFL)=0
NEFL=NEFL-1
CALL PLOTOT(1)
CALL HEDR
DO IELEM=1,NEFL
CALL FILLEM(NEFLAG(IELEM))
ENDDO
!
! ELSEIF(IFLAG .EQ. 'r') THEN
! CALL PLOTS(0)
! CALL PLOTOT
! GO TO 4
ELSEIF(IFLAG .EQ. 'q' .OR. IFLAG .EQ. 'e') THEN
GO TO 152
!
ELSE
!IPK JAN98 WRITE(*,*) CHAR(7),CHAR(7)
ENDIF
!
GOTO 4
!
!
150 CONTINUE
! IPSWO=IPSW
! IPSW=4
! CALL PLOTS(0)
!ipk oct96 DO 151 I=1,9
!ipk oct96 IPSW(I)=0
!ipk oct96 151 CONTINUE
!ipk oct96 IPSW(4)=1
!ipk nov97 add (1)
CALL PLOTOT(1)
! IPSW=IPSWO
!
! Define NEF and process elements
!
152 CONTINUE
DO N=1,NE
DO M=1,8
NOPSV(N,M)=NOP(N,M)
ENDDO
IMATSV(N)=IMAT(N)
ENDDO
NPUNDO=0
NEUNDO=0
NESAV=NE
NEFSAV=NENTRY
IF(NENTRY .GT. 0) THEN
DO N=1,NENTRY
DO M=1,3
NEFSV(N,M)=NEF(N,M)
ENDDO
ENDDO
ENDIF
ITYPSV=ITYP
DO 250 NN=1,NEFL
ITYP=ITYPSV
N=NEFLAG(NN)
IF(IMAT(N) .GT. 900 .AND. IMAT(N) .LT. 904) GO TO 250
! IF(IMAT(N) .EQ. 999) ITYP=1
NCN=NCORN(N)
!
! Split a one-dimensional element in two
!
IF(NCN .EQ. 3) THEN
N1=NOP(N,1)
N2=NOP(N,2)
N3=NOP(N,3)
IF(NOP(N,2) .EQ. 0) THEN
CALL GETNOD(N2)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=N2
ELSEIF(INEW(N2) .EQ. 1) THEN
GO TO 153
ENDIF
CORD(N2,1)=(CORD(N1,1)+CORD(N3,1))/2.
CORD(N2,2)=(CORD(N1,2)+CORD(N3,2))/2.
IF(LOCK(N1) .EQ. 1 .AND. LOCK(N3) .EQ. 1) LOCK(N2)=1
XUSR(N2) = CORD(N2,1)*TXSCAL - XS
YUSR(N2) = CORD(N2,2)*TXSCAL - YS
INEW(N2) = 1
INSKP(N2) =0
153 CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
NOP(NEM,3)=N3
NOP(N,2)=0
NOP(N,3)=N2
NOP(NEM,1)=N2
NOP(NEM,2)=0.
NOP(NEM,3)=N3
IMAT(NEM)=IMAT(N)
IESKP(NEM)=0
NCORN(NEM)=3
!ipk jan98
IF(IDELV .EQ. 1) then
WD(N2)=-9999.
WIDTH(N2)=0.
SS1(N2)=0.
SS2(N2)=0.
WIDS(N2)=0.
ELSE
WD(N2)=(WD(N1)+WD(N3))/2.
WIDTH(N2)=(WIDTH(N1)+WIDTH(N3))/2.
SS1(N2)=(SS1(N1)+SS1(N3))/2.
SS2(N2)=(SS2(N1)+SS2(N3))/2.
WIDS(N2)=(WIDS(N1)+WIDS(N3))/2.
IF(ICRIN .EQ. 23) CALL COMPWGT
ENDIF
GO TO 250
ENDIF
!
! Setup for each type of refinement
!
!ipk jan08
IF(ITYP .EQ. 0) THEN
!
! Full refinement all nodes are eligible
!
! IF(imat(n) .eq. 999) then
! IELGB(2)=2
! IELGB(4)=0
! IELGB(6)=2
! IELGB(8)=0
! ELSE
DO M=2,NCN
IELGB(M)=1
ENDDO
! ENDIF
ELSEIF(ITYP .EQ. 1 .OR. ITYP .EQ. 2) THEN
!
! Setup for long or short side refinement
!
IF(ITYP .EQ. 1) THEN
DISTLL=0.
DISTL=0.
ELSE
DISTLL=-VOID
DISTL=-VOID
ENDIF
!
! Sort out longest or shortest sides
!
DO 165 M=2,NCN,2
IELGB(M)=0
N1=NOP(N,M-1)
N2=MOD(M,NCN)+1
N2=NOP(N,N2)
DSEP=DIST(N1,N2)
IF(ITYP .EQ. 1) THEN
IF(DISTLL .LT. DSEP) THEN
! Separation greater DISTLL
IF(DISTLL .GT. 0.) THEN
! DISTLL already exists then move it down the line
DISTL=DISTLL
NDS=NDSS
ENDIF
! Save separation
DISTLL=DSEP
NDSS=M
GO TO 165
ELSEIF(DISTL .LT. DSEP) THEN
! 2nd longest
DISTL=DSEP
NDS=M
ENDIF
ELSE
IF(DSEP .LT. DISTLL) THEN
! Separation less than DISTLL
IF(DISTLL .LT. -VDX) THEN
! DISTLL already exists then move it up the line
DISTL=DISTLL
NDS=NDSS
ENDIF
DISTLL=DSEP
NDSS=M
GO TO 165
ELSEIF(DSEP .LT. DISTL) THEN
! 2nd shortest
DISTL=DSEP
NDS=M
ENDIF
ENDIF
165 CONTINUE
IELGB(NDSS)=2
IELGB(NDS)=2
ELSEIF(ITYP .EQ. 3) THEN
!ipk jan98 IF(NCN .EQ. 8) CALL SPLIT(N)
IF(NCN .GT. 5) CALL SPLIT(N)
GO TO 250
ELSEIF(ITYP .EQ. 4) THEN
NPL=NEFLAG(NN+1)
CALL REVERS(N,NPL)
GO TO 255
ENDIF
!
! Loop through element sides
!
DO 200 M=2,NCN,2
IF(IELGB(M) .EQ. 0) GO TO 200
N1=NOP(N,M-1)
N3=MOD(M+1,NCN)
N3=NOP(N,N3)
!
! Search table for N1
!
IF(NENTRY .EQ. 0) GO TO 182
DO 180 J=1,NENTRY
IF(N1 .EQ. NEF(J,3) .AND. N3 .EQ. NEF(J,1)) THEN
!
! We have found match so use this info
!
NOP(N,M)=NEF(J,2)
!
! For regular ops remove value in NEF(J,1) so that it seems blank and s
! otherwise set value negative
IF(IELGB(M) .EQ. 1) THEN
NEF(J,1)=0
ELSE
NEF(J,1)=-NEF(J,1)
ENDIF
GO TO 200
ENDIF
180 CONTINUE
182 CONTINUE
!
! Define a node, enter it, initialize it, and make entry in NEF
!
IF(IMAT(N) .EQ. 999 .AND. (M .EQ. 4 .OR. M .EQ. 8)) GO TO 200
IF(NOP(N,M) .EQ. 0) THEN
CALL GETNOD(N2)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=N2
NOP(N,M)=N2
CORD(N2,1)=(CORD(N1,1)+CORD(N3,1))/2.
CORD(N2,2)=(CORD(N1,2)+CORD(N3,2))/2.
IF(LOCK(N1) .EQ. 1 .AND. LOCK(N3) .EQ. 1) LOCK(N2)=1
XUSR(N2) = CORD(N2,1)*TXSCAL - XS
YUSR(N2) = CORD(N2,2)*TXSCAL - YS
INEW(N2) = 1
INSKP(N2) =0
ELSE
N2=NOP(N,M)
IF(INEW(N2) .NE. 1) THEN
CORD(N2,1)=(CORD(N1,1)+CORD(N3,1))/2.
CORD(N2,2)=(CORD(N1,2)+CORD(N3,2))/2.
XUSR(N2) = CORD(N2,1)*TXSCAL - XS
YUSR(N2) = CORD(N2,2)*TXSCAL - YS
INEW(N2) = 1
INSKP(N2) =0
ENDIF
ENDIF
!ipk jan98
IF(IDELV .EQ. 1) then
WD(N2)=-9999.
ELSE
WD(N2)=(WD(N1)+WD(N3))/2.
ENDIF
IF(M .EQ. 2 .AND. IMAT(N) .EQ. 999) THEN
WIDTH(N2)=(WIDTH(N1)+WIDTH(N3))/2.
SS1(N2)=(SS1(N1)+SS1(N3))/2.
SS2(N2)=(SS2(N1)+SS2(N3))/2.
WIDS(N2)=(WIDS(N1)+WIDS(N3))/2.
ELSE
WIDTH(N2)=0.
SS1(N2)=0.
SS2(N2)=0.
WIDS(N2)=0.
ENDIF
NENTRY=NENTRY+1
NEF(NENTRY,1)=N1
NEF(NENTRY,2)=N2
NEF(NENTRY,3)=N3
200 CONTINUE
IF(ITYP .GT. 0) GO TO 250
!
! Copy NOP into temporary NTRAN for processing then delete element
!
DO 220 K=1,8
NTRAN(K)=NOP(N,K)
NOP(N,K)=0
220 CONTINUE
NRMAT=IMAT(N)
IMAT(N)=0
IESKP(N)=-1
NTYP=1
NELAST= MIN(NELAST,N)
IF(NCN .EQ. 8) THEN
IF(NRMAT .EQ. 999) THEN
IF(NTRAN(2) .EQ. 0) THEN
CALL GETNOD(N2)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=N2
N1=NTRAN(1)
N3=NTRAN(3)
CORD(N2,1)=(CORD(N1,1)+CORD(N3,1))/2.
CORD(N2,2)=(CORD(N1,2)+CORD(N3,2))/2.
INEW(N2) = 1
IF(LOCK(N1) .EQ. 1 .AND. LOCK(N3) .EQ. 1) LOCK(N2)=1
NTRAN(2)=N2
WD(N2)=(WD(N1)+WD(N3))/2.
WIDTH(N2)=(WIDTH(N1)+WIDTH(N3))/2.
SS1(N2)=(SS1(N1)+SS1(N3))/2.
SS2(N2)=(SS2(N1)+SS2(N3))/2.
WIDS(N2)=(WIDS(N1)+WIDS(N3))/2.
ENDIF
IF(NTRAN(6) .EQ. 0) THEN
CALL GETNOD(N6)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=N6
N5=NTRAN(5)
N7=NTRAN(7)
CORD(N6,1)=(CORD(N5,1)+CORD(N7,1))/2.
CORD(N6,2)=(CORD(N5,2)+CORD(N7,2))/2.
INEW(N6) = 1
IF(LOCK(N5) .EQ. 1 .AND. LOCK(N7) .EQ. 1) LOCK(N6)=1
NTRAN(6)=N6
WD(N6)=(WD(N5)+WD(N7))/2.
WIDTH(N6)=(WIDTH(N5)+WIDTH(N7))/2.
SS1(N6)=(SS1(N5)+SS1(N7))/2.
SS2(N6)=(SS2(N5)+SS2(N7))/2.
WIDS(N6)=(WIDS(N5)+WIDS(N7))/2.
ENDIF
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
NOP(NEM,1)=NTRAN(1)
NOP(NEM,3)=NTRAN(2)
NOP(NEM,5)=NTRAN(6)
NOP(NEM,7)=NTRAN(7)
IMAT(NEM)=999
IESKP(NEM)=0
NCORN(NEM)=8
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
NOP(NEM,1)=NTRAN(2)
NOP(NEM,3)=NTRAN(3)
NOP(NEM,5)=NTRAN(5)
NOP(NEM,7)=NTRAN(6)
IMAT(NEM)=999
IESKP(NEM)=0
NCORN(NEM)=8
ELSE
CALL GETNOD(N2)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=N2
CORD(N2,1)=(CORD(NTRAN(1),1)+CORD(NTRAN(3),1) &
& +CORD(NTRAN(5),1)+CORD(NTRAN(7),1))/4.
CORD(N2,2)=(CORD(NTRAN(1),2)+CORD(NTRAN(3),2) &
& +CORD(NTRAN(5),2)+CORD(NTRAN(7),2))/4.
INEW(N2) = 1
IF(LOCK(NTRAN(1)) .EQ. 1 .AND. LOCK(NTRAN(3)) .EQ. 1 .AND. &
& LOCK(NTRAN(5)) .EQ. 1 .AND. LOCK(NTRAN(7)) .EQ. 1) LOCK(N2)=1
!ipk jan98
IF(IDELV .EQ. 1) then
WD(N2)=-9999.
ELSE
WD(N2) =(WD(NTRAN(1))+WD(NTRAN(3)) &
& +WD(NTRAN(5))+WD(NTRAN(7)))/4.
ENDIF
WIDTH(N2)=0.
SS1(N2)=0.
SS2(N2)=0.
WIDS(N2)=0.
XUSR(N2) = CORD(N2,1)*TXSCAL - XS
YUSR(N2) = CORD(N2,2)*TXSCAL - YS
NTRAN(9)=N2
INSKP(N2)=0
CALL RGEN(NTRAN,NTYP,NRMAT)
ENDIF
ELSE
CALL TGEN(NTRAN,NTYP,NRMAT)
ENDIF
IF(MOD(NN,20) .EQ. 0) THEN
!
! Compress NEF for later use
!
NCT=0
DO 245 N=1,NENTRY
IF(NEF(N,1) .NE. 0) THEN
NCT=NCT+1
NEF(NCT,1)=NEF(N,1)
NEF(NCT,2)=NEF(N,2)
NEF(NCT,3)=NEF(N,3)
ENDIF
245 CONTINUE
NENTRY=NCT
ENDIF
250 END DO
255 CONTINUE
IF(ITYP .GT. 2) THEN
!ipk nov97 add (1)
call plotot(1)
NEFL=0
RETURN
ENDIF
!
! Process the ITYP = 1 or 2 situation
!
IF(ITYP .GT. 0) THEN
CALL CLENUP(ITYP)
ENDIF
!
! Search for left over entries NEF
!
DO 600 I=1,NENTRY
DO 500 N=1,NE
IF(IMAT(N) .EQ. 0) GO TO 500
NCN=NCORN(N)
!ipk sep99 add test for line element
if(ncn .eq. 3) then
if(nef(i,2) .eq. nop(n,2)) go to 600
go to 500
endif
!
! Loop on sides
!
DO 400 K=2,NCN,2
IF(NOP(N,K-1) .EQ. NEF(I,3)) THEN
KP=MOD(K+1,NCN)
IF(NOP(N,KP) .EQ. ABS(NEF(I,1))) THEN
!
! We have a match, quit search for this entry
!
GO TO 600
ENDIF
ENDIF
400 CONTINUE
500 CONTINUE
!
! No match this must be a boundary eliminate NEF value
!
NEF(I,1)=0
NEF(I,3)=0
600 END DO
!
! Now compress remaining NEF for later use
!
NCT=0
DO 700 N=1,NENTRY
IF(NEF(N,1) .GT. 0) THEN
NCT=NCT+1
NEF(NCT,1)=NEF(N,1)
NEF(NCT,2)=NEF(N,2)
NEF(NCT,3)=NEF(N,3)
ENDIF
700 END DO
NENTRY=NCT
NEFL=0
RETURN
END
!
SUBROUTINE CLENUP(ITYP)
!
! Clean up transitions on the boundary of the refined area
!
!
USE BLK1MOD
! INCLUDE 'BLK1.COM'
INCLUDE 'TXFRM.COM'
!IPK MAY02 COMMON /TXFRM/ XS, YS, TXSCAL
!
DIMENSION NTEMP(9),NTRAN(9),NSWT(8)
!
! First loop through elements looking for transitions
!
IF(ITYP .EQ. 0) THEN
NEO=NE
ELSE
NEO=NEFL
ENDIF
! DO KN=1,NEO
! WRITE(234,*) KN,NEFLAG(KN),NEF(KN,1),NEF(KN,2),NEF(KN,3)
! ENDDO
DO 500 KN=1,NEO
IF(ITYP .EQ. 0) THEN
N=KN
IF(IMAT(N) .EQ. 0) GO TO 500
ELSE
N=NEFLAG(KN)
ENDIF
NCN=NCORN(N)
!ipk sep99 add test for line element
if(ncn .eq. 3) then
do i=1,nentry
if(nop(n,2) .eq. nef(i,2)) then
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
nop(nem,1)=nef(I,2)
nop(nem,3)=nef(I,3)
imat(nem)=imat(n)
ncorn(nem)=3
IESKP(NEM)=0
IERC=0
CALL PLTELM(NEM,IERC)
nop(n,2)=0
nop(n,3)=nef(I,2)
go to 500
endif
enddo
go to 500
endif
!
! Loop on sides
!
IFND=0
NSWT(8)=0
DO 400 K=2,NCN,2
!
! Search for left over entry in NEF
!
DO 350 I=1,NENTRY
IF(NOP(N,K-1) .EQ. NEF(I,3)) THEN
KP=MOD(K+1,NCN)
IF(NOP(N,KP) .EQ. ABS(NEF(I,1))) THEN
!
! We have a match, start building TEMP
!
NTEMP(K-1)=NEF(I,3)
NTEMP(K)=NEF(I,2)
NSWT(K)=1
IFND=1
GO TO 400
ENDIF
ENDIF
IF(ITYP .GT. 0) THEN
IF(NOP(N,K-1) .EQ. ABS(NEF(I,1))) THEN
KP=MOD(K+1,NCN)
IF(NOP(N,KP) .EQ. NEF(I,3)) THEN
!
! We have a match, start building TEMP
!
NTEMP(K-1)=ABS(NEF(I,1))
NTEMP(K)=NEF(I,2)
NSWT(K)=1
IFND=1
GO TO 400
ENDIF
ENDIF
ENDIF
350 CONTINUE
!
! No match copy old values
!
NTEMP(K-1)=NOP(N,K-1)
NTEMP(K)=NOP(N,K)
NSWT(K)=0
400 CONTINUE
IF(IFND .EQ. 0) GO TO 500
!
! Now test for match
!
NTOT=NSWT(2)+NSWT(4)+NSWT(6)+NSWT(8)
IF(NTOT .EQ. 0) GO TO 500
!
! Delete element
!
DO 420 K=1,8
NOP(N,K)=0
420 CONTINUE
NRMAT=IMAT(N)
IMAT(N)=0
NELAST=MIN(NELAST,N)
!
! Work with triangles first
!
IF(NCN .EQ. 6) THEN
!
! Determine transition type and prepare to rotate connections
!
IF(NTOT .EQ. 1) THEN
NTYP=3
IF(NSWT(2) .EQ. 1) THEN
ISHIFT=0
ELSEIF(NSWT(4) .EQ. 1) THEN
ISHIFT=2
ELSEIF(NSWT(6) .EQ. 1) THEN
ISHIFT=4
ENDIF
ELSEIF(NTOT .EQ. 2) THEN
NTYP=2
IF(NSWT(2) .EQ. 0) THEN
ISHIFT=2
ELSEIF(NSWT(4) .EQ. 0) THEN
ISHIFT=4
ELSEIF(NSWT(6) .EQ. 0) THEN
ISHIFT=0
ENDIF
ELSE
NTYP=1
ISHIFT=0
ENDIF
!
! Now rotate so that first mid node is refined
!
DO 430 K=1,NCN
KS=MOD(K+ISHIFT,NCN)
IF(KS .EQ. 0) KS=NCN
NTRAN(K)=NTEMP(KS)
430 CONTINUE
!
! Now generate transition refined elements
!
CALL TGEN(NTRAN,NTYP,NRMAT)
!
! Now work on quadrilateral elements
!
ELSE
!
! Determine transition type and prepare to rotate connections
!
IF(NTOT .EQ. 1) THEN
NTYP=2
IF(NSWT(2) .EQ. 1) THEN
ISHIFT=0
ELSEIF(NSWT(4) .EQ. 1) THEN
ISHIFT=2
ELSEIF(NSWT(6) .EQ. 1) THEN
ISHIFT=4
ELSEIF(NSWT(8) .EQ. 1) THEN
ISHIFT=6
ENDIF
ELSEIF(NTOT .EQ. 2) THEN
IF(NSWT(2) .EQ. 1) THEN
IF(NSWT(4) .EQ. 1) THEN
NTYP=3
ISHIFT=0
ELSEIF(NSWT(6) .EQ. 1) THEN
NTYP=4
ISHIFT=0
ELSE
NTYP=3
ISHIFT=6
ENDIF
ELSEIF(NSWT(4) .EQ. 1) THEN
IF(NSWT(6) .EQ. 1) THEN
NTYP=3
ISHIFT=2
ELSEIF(NSWT(8) .EQ. 1) THEN
NTYP=4
ISHIFT=2
ENDIF
ELSE
NTYP=3
ISHIFT=4
ENDIF
ELSEIF(NTOT .EQ. 3) THEN
NTYP=5
IF(NSWT(2) .EQ. 0) THEN
ISHIFT=2
ELSEIF(NSWT(4) .EQ. 0) THEN
ISHIFT=4
ELSEIF(NSWT(6) .EQ. 0) THEN
ISHIFT=6
ELSEIF(NSWT(8) .EQ. 0) THEN
ISHIFT=0
ENDIF
ELSE
NTYP=1
ISHIFT=0
ENDIF
!
! Now rotate so that first mid node is refined
!
DO 450 K=1,NCN
KS=MOD(K+ISHIFT,NCN)
IF(KS .EQ. 0) KS=NCN
NTRAN(K)=NTEMP(KS)
450 CONTINUE
!
IF(NTYP .EQ. 1 .OR. NTYP .EQ. 5) THEN
!
! If appropriate define a new node at the centroid
!
CALL GETNOD(N2)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=N2
CORD(N2,1)=(CORD(NTEMP(1),1)+CORD(NTEMP(3),1) &
& +CORD(NTEMP(5),1)+CORD(NTEMP(7),1))/4.
CORD(N2,2)=(CORD(NTEMP(1),2)+CORD(NTEMP(3),2) &
& +CORD(NTEMP(5),2)+CORD(NTEMP(7),2))/4.
IF(LOCK(NTEMP(1)) .EQ. 1 .AND. LOCK(NTEMP(3)) .EQ. 1 .AND. &
& LOCK(NTEMP(5)) .EQ. 1 .AND. LOCK(NTEMP(7)) .EQ. 1) LOCK(N2)=1
INEW(N2) = 1
!ipk jan98
IF(IDELV .EQ. 1) then
WD(N2)=-9999.
ELSE
WD(N2)= (WD(NTEMP(1))+WD(NTEMP(3)) &
& +WD(NTEMP(5))+WD(NTEMP(7)))/4.
ENDIF
WIDTH(N2)=0.
SS1(N2)=0.
SS2(N2)=0.
WIDS(N2)=0.
XUSR(N2) = CORD(N2,1)*TXSCAL - XS
YUSR(N2) = CORD(N2,2)*TXSCAL - YS
NTRAN(9)=N2
INSKP(N2)=0
!
! Now generate transition refined elements
!
ENDIF
CALL RGEN(NTRAN,NTYP,NRMAT)
ENDIF
500 END DO
IF(ITYP .EQ. 0) THEN
NENTRY=0
ELSE
DO 600 I=1,NENTRY
IF(NEF(I,1) .LT. 0) NEF(I,1)=0
600 CONTINUE
ENDIF
RETURN
END
!
SUBROUTINE RGEN(NTRAN,NTYP,NRMAT)
!
! Routine to refine quadrilateral elements
!
USE BLK1MOD
! INCLUDE 'BLK1.COM'
!
! IRGEN contains pointers to the various connections
!
INTEGER*2 IRGEN
DIMENSION NTRAN(9),IRGEN(8,5,5)
!
! DATA IRGEN /1,0,2,0,9,0,8,0,3,0,4,0,9,0,2,0,5,0,6,0,9,0,4,0, &
! & 7,0,8,0,9,0,6,0,8*0, &
DATA IRGEN /1,0,2,0,9,0,8,0,2,0,3,0,4,0,9,0,9,0,4,0,5,0,6,0, &
& 8,0,9,0,6,0,7,0,8*0, &
& 1,0,2,0,7,8,0,0,3,4,5,0,2,0,0,0,5,6,7,0,2,0,0,0,16*0, &
& 1,0,2,0,7,8,0,0,3,0,4,0,2,0,0,0,5,6,7,0,4,0,0,0, &
& 7,0,2,0,4,0,0,0,8*0, &
& 1,0,2,0,6,0,7,8,2,0,3,4,5,0,6,0,24*0, &
& 1,0,2,0,9,0,0,0,3,0,4,0,9,0,2,0,5,0,6,0,9,0,4,0, &
& 7,0,9,0,6,0,0,0,7,8,1,0,9,0,0,0/
!
DO 300 N=1,5
IF(IRGEN(1,N,NTYP) .EQ. 0) GO TO 310
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
DO 250 K=1,7,2
INN=IRGEN(K,N,NTYP)
INP=IRGEN(K+1,N,NTYP)
IF(INN .GT. 0) INN=NTRAN(INN)
IF(INP .GT. 0) INP=NTRAN(INP)
NOP(NEM,K)=INN
NOP(NEM,K+1)=INP
250 CONTINUE
IF(NOP(NEM,7) .EQ. 0) THEN
NCORN(NEM)=6
ELSE
NCORN(NEM)=8
ENDIF
IMAT(NEM)=NRMAT
IESKP(NEM)=0
!IPK JAN98
IERC=0
CALL PLTELM(NEM,IERC)
300 END DO
310 CONTINUE
RETURN
END
!
SUBROUTINE TGEN(NTRAN,NTYP,NRMAT)
!
! Routine to refine triangular elements
!
USE BLK1MOD
! INCLUDE 'BLK1.COM'
!
! ITGEN contains pointers to the various connections
!
INTEGER*2 ITGEN
DIMENSION NTRAN(9),ITGEN(8,4,3)
!
DATA ITGEN /1,0,2,0,6,0,0,0,3,0,4,0,2,0,0,0, &
& 5,0,6,0,4,0,0,0,2,0,4,0,6,0,0,0, &
& 1,0,2,0,4,0,5,6,2,0,3,0,4,0,0,0,16*0, &
& 1,0,2,0,5,6,0,0,3,4,5,0,2,0,0,0,16*0/
!
DO 300 N=1,4
IF(ITGEN(1,N,NTYP) .EQ. 0) GO TO 310
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
DO 250 K=1,7,2
INN=ITGEN(K,N,NTYP)
INP=ITGEN(K+1,N,NTYP)
IF(INN .GT. 0) INN=NTRAN(INN)
IF(INP .GT. 0) INP=NTRAN(INP)
NOP(NEM,K)=INN
NOP(NEM,K+1)=INP
250 CONTINUE
IF(NOP(NEM,7) .EQ. 0) THEN
NCORN(NEM)=6
ELSE
NCORN(NEM)=8
ENDIF
IMAT(NEM)=NRMAT
IESKP(NEM)=0
IERC=0
CALL PLTELM(NEM,IERC)
300 END DO
310 CONTINUE
RETURN
END
SUBROUTINE SPLIT(N)
!
! Routine to split quadrilateral elements in two
!
USE BLK1MOD
! INCLUDE 'BLK1.COM'
INCLUDE 'TXFRM.COM'
!IPK MAY02 COMMON /TXFRM/ XS, YS, TXSCAL
!
DIST(N1,N2)=SQRT((CORD(N1,1)-CORD(N2,1))**2 &
& +(CORD(N1,2)-CORD(N2,2))**2)
if(nop(n,7) .eq. 0) go to 100
!
! Loop around element looking for longest diagonal
!
L1=NOP(N,1)
L5=NOP(N,5)
D15=DIST(L1,L5)
L3=NOP(N,3)
L7=NOP(N,7)
D37=DIST(L3,L7)
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
IF(D15 .LT. D37) THEN
NOP(NEM,1)=L1
NOP(NEM,2)=0
NOP(NEM,3)=L5
NOP(NEM,4)=NOP(N,6)
NOP(NEM,5)=L7
NOP(NEM,6)=NOP(N,8)
IMAT(NEM)=IMAT(N)
IESKP(NEM)=0
NCORN(NEM)=6
NOP(N,6)=0
NOP(N,7)=0
NOP(N,8)=0
NCORN(N)=6
ELSE
NOP(NEM,1)=L3
NOP(NEM,2)=NOP(N,4)
NOP(NEM,3)=L5
NOP(NEM,4)=NOP(N,6)
NOP(NEM,5)=L7
NOP(NEM,6)=0
IMAT(NEM)=IMAT(N)
IESKP(NEM)=0
NCORN(NEM)=6
NOP(N,4)=0
NOP(N,5)=L7
NOP(N,6)=NOP(N,8)
NOP(N,7)=0
NOP(N,8)=0
NCORN(N)=6
ENDIF
! call plotot
RETURN
100 continue
!
! triangle split
!
l1=nop(n,1)
l3=nop(n,3)
l5=nop(n,5)
d13=dist(l1,l3)
d35=dist(l3,l5)
d51=dist(l5,l1)
CALL GETELM(NEM)
NEUNDO=NEUNDO+1
IELDEL(NEUNDO)=NEM
IMAT(NEM)=IMAT(N)
IESKP(NEM)=0
NCORN(NEM)=6
write(90,*) l1,l3,l5,d13,d35,d51,nentry
if(d13 .gt. d35) then
if(d13 .gt. d51) then
!
! Search table for L1
!
IF(NENTRY .NE. 0) THEN
DO J=1,NENTRY
IF(L1 .EQ. NEF(J,3) .AND. L3 .EQ. NEF(J,1)) THEN
!
! We have found match so use this info
!
NOP(N,2)=NEF(J,2)
NEWND=NEF(J,2)
!
! For regular ops remove value in NEF(J,1) so that it seems blank and s
! otherwise set value negative
! IF(IELGB(2) .EQ. 1) THEN
! NEF(J,1)=0
! ELSE
NEF(J,1)=-NEF(J,1)
! ENDIF
GO TO 200
ENDIF
ENDDO
ENDIF
!
! Define a node, enter it, initialize it, and make entry in NEF
!
IF(NOP(N,2) .EQ. 0) THEN
CALL GETNOD(NEWND)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=NEWND
NOP(N,2)=NEWND
CORD(NEWND,1)=(CORD(L1,1)+CORD(L3,1))/2.
CORD(NEWND,2)=(CORD(L1,2)+CORD(L3,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
IF(LOCK(L1) .EQ. 1 .AND. LOCK(L3) .EQ. 1 ) LOCK(NEWND)=1
INSKP(NEWND) =0
ELSE
NEWND=NOP(N,2)
IF(INEW(NEWND) .NE. 1) THEN
CORD(NEWND,1)=(CORD(L1,1)+CORD(L3,1))/2.
CORD(NEWND,2)=(CORD(L1,2)+CORD(L3,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
INSKP(NEWND) =0
ENDIF
ENDIF
!ipk jan98
IF(IDELV .EQ. 1) then
WD(NEWND)=-9999.
ELSE
WD(NEWND)=(WD(L1)+WD(L3))/2.
ENDIF
WIDTH(NEWND)=0.
SS1(NEWND)=0.
SS2(NEWND)=0.
WIDS(NEWND)=0.
NENTRY=NENTRY+1
NEF(NENTRY,1)=L1
NEF(NENTRY,2)=NEWND
NEF(NENTRY,3)=L3
200 CONTINUE
nop(nem,1)=nop(n,1)
nop(nem,3)=newnd
nop(nem,5)=nop(n,5)
nop(nem,6)=nop(n,6)
nop(n,1)=newnd
nop(n,2)=0
nop(n,6)=0
else
!
! Search table for L5
!
IF(NENTRY .NE. 0) THEN
DO J=1,NENTRY
IF(L5 .EQ. NEF(J,3) .AND. L1 .EQ. NEF(J,1)) THEN
!
! We have found match so use this info
!
NOP(N,2)=NEF(J,2)
NEWND=NEF(J,2)
!
! For regular ops remove value in NEF(J,1) so that it seems blank and s
! otherwise set value negative
! IF(IELGB(2) .EQ. 1) THEN
! NEF(J,1)=0
! ELSE
NEF(J,1)=-NEF(J,1)
! ENDIF
GO TO 300
ENDIF
ENDDO
ENDIF
!
! Define a node, enter it, initialize it, and make entry in NEF
!
IF(NOP(N,6) .EQ. 0) THEN
CALL GETNOD(NEWND)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=NEWND
NOP(N,6)=NEWND
CORD(NEWND,1)=(CORD(L5,1)+CORD(L1,1))/2.
CORD(NEWND,2)=(CORD(L5,2)+CORD(L1,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
IF(LOCK(L1) .EQ. 1 .AND. LOCK(L5) .EQ. 1) LOCK(NEWND)=1
INSKP(NEWND) =0
ELSE
NEWND=NOP(N,6)
IF(INEW(NEWND) .NE. 1) THEN
CORD(NEWND,1)=(CORD(L5,1)+CORD(L1,1))/2.
CORD(NEWND,2)=(CORD(L5,2)+CORD(L1,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
INSKP(NEWND) =0
ENDIF
ENDIF
!ipk jan98
IF(IDELV .EQ. 1) then
WD(NEWND)=-9999.
ELSE
WD(NEWND)=(WD(L5)+WD(L1))/2.
ENDIF
WIDTH(NEWND)=0.
SS1(NEWND)=0.
SS2(NEWND)=0.
WIDS(NEWND)=0.
NENTRY=NENTRY+1
NEF(NENTRY,1)=L5
NEF(NENTRY,2)=NEWND
NEF(NENTRY,3)=L1
300 CONTINUE
nop(nem,1)=nop(n,1)
nop(nem,2)=nop(n,2)
nop(nem,3)=nop(n,3)
nop(nem,5)=newnd
nop(n,1)=newnd
nop(n,2)=0
nop(n,6)=0
endif
elseif(d35 .gt. d51) then
!
! Search table for L3
!
IF(NENTRY .NE. 0) THEN
DO J=1,NENTRY
IF(L3 .EQ. NEF(J,3) .AND. L5 .EQ. NEF(J,1)) THEN
!
! We have found match so use this info
!
NOP(N,4)=NEF(J,2)
NEWND=NEF(J,2)
!
! For regular ops remove value in NEF(J,1) so that it seems blank and s
! otherwise set value negative
! IF(IELGB(2) .EQ. 1) THEN
! NEF(J,1)=0
! ELSE
NEF(J,1)=-NEF(J,1)
! ENDIF
GO TO 400
ENDIF
ENDDO
ENDIF
!
! Define a node, enter it, initialize it, and make entry in NEF
!
IF(NOP(N,4) .EQ. 0) THEN
CALL GETNOD(NEWND)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=NEWND
NOP(N,4)=NEWND
CORD(NEWND,1)=(CORD(L3,1)+CORD(L5,1))/2.
CORD(NEWND,2)=(CORD(L3,2)+CORD(L5,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
IF(LOCK(L3) .EQ. 1 .AND. LOCK(L5) .EQ. 1) LOCK(NEWND)=1
INSKP(NEWND) =0
ELSE
NEWND=NOP(N,4)
IF(INEW(NEWND) .NE. 1) THEN
CORD(NEWND,1)=(CORD(L3,1)+CORD(L5,1))/2.
CORD(NEWND,2)=(CORD(L3,2)+CORD(L5,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
INSKP(NEWND) =0
ENDIF
ENDIF
!ipk jan98
IF(IDELV .EQ. 1) then
WD(NEWND)=-9999.
ELSE
WD(NEWND)=(WD(L3)+WD(L5))/2.
ENDIF
WIDTH(NEWND)=0.
SS1(NEWND)=0.
SS2(NEWND)=0.
WIDS(NEWND)=0.
NENTRY=NENTRY+1
NEF(NENTRY,1)=L3
NEF(NENTRY,2)=NEWND
NEF(NENTRY,3)=L5
400 CONTINUE
nop(nem,1)=nop(n,1)
nop(nem,2)=nop(n,2)
nop(nem,3)=nop(n,3)
nop(nem,5)=newnd
nop(n,3)=newnd
nop(n,2)=0
nop(n,4)=0
else
!
! Search table for L5
!
IF(NENTRY .NE. 0) THEN
DO J=1,NENTRY
IF(L5 .EQ. NEF(J,3) .AND. L1 .EQ. NEF(J,1)) THEN
!
! We have found match so use this info
!
NOP(N,2)=NEF(J,2)
NEWND=NEF(J,2)
!
! For regular ops remove value in NEF(J,1) so that it seems blank and s
! otherwise set value negative
! IF(IELGB(2) .EQ. 1) THEN
! NEF(J,1)=0
! ELSE
NEF(J,1)=-NEF(J,1)
! ENDIF
GO TO 500
ENDIF
ENDDO
ENDIF
!
! Define a node, enter it, initialize it, and make entry in NEF
!
IF(NOP(N,6) .EQ. 0) THEN
CALL GETNOD(NEWND)
NPUNDO=NPUNDO+1
NODDEL(NPUNDO)=NEWND
NOP(N,6)=NEWND
CORD(NEWND,1)=(CORD(L5,1)+CORD(L1,1))/2.
CORD(NEWND,2)=(CORD(L5,2)+CORD(L1,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
IF(LOCK(L1) .EQ. 1 .AND. LOCK(L5) .EQ. 1) LOCK(NEWND)=1
INSKP(NEWND) =0
ELSE
NEWND=NOP(N,6)
IF(INEW(NEWND) .NE. 1) THEN
CORD(NEWND,1)=(CORD(L5,1)+CORD(L1,1))/2.
CORD(NEWND,2)=(CORD(L5,2)+CORD(L1,2))/2.
XUSR(NEWND) = CORD(NEWND,1)*TXSCAL - XS
YUSR(NEWND) = CORD(NEWND,2)*TXSCAL - YS
INEW(NEWND) = 1
INSKP(NEWND) =0
ENDIF
ENDIF
!ipk jan98
IF(IDELV .EQ. 1) then
WD(NEWND)=-9999.
ELSE
WD(NEWND)=(WD(L5)+WD(L1))/2.
ENDIF
WIDTH(NEWND)=0.
SS1(NEWND)=0.
SS2(NEWND)=0.
WIDS(NEWND)=0.
NENTRY=NENTRY+1
NEF(NENTRY,1)=L5
NEF(NENTRY,2)=NEWND
NEF(NENTRY,3)=L1
500 CONTINUE
nop(nem,1)=nop(n,1)
nop(nem,2)=nop(n,2)
nop(nem,3)=nop(n,3)
nop(nem,5)=newnd
nop(n,1)=newnd
nop(n,2)=0
nop(n,6)=0
endif
return
END
SUBROUTINE REVERS(N1,N2)
!
! Routine to reverse diagonal of two quadrilateral elements
!
USE BLK1MOD
USE BLK2MOD
! INCLUDE 'BLK1.COM'
INCLUDE 'TXFRM.COM'
!IPK MAY02 COMMON /TXFRM/ XS, YS, TXSCAL
!
! Search for common nodes on the elements
!
DO 300 M=1,NCORN(N1),2
J=NOP(N1,M)
DO 250 MM=1,NCORN(N2),2
JJ=NOP(N2,MM)
IF(JJ .EQ. J) THEN
!
! We have a match find the other nodes around element
!
MID1=M+1
JMID1=NOP(N1,MID1)
write(90,*) n1,mid1,jmid1
MID2=M+3
IF(M .EQ. 5) MID2=2
JMID2=NOP(N1,MID2)
MID3=M+5
IF(MID3 .GT. 6) MID3=MID3-6
JMID3=NOP(N1,MID3)
!
! Now find the other node
!
M1=M+2
IF(M1 .GT. 6) M1=1
J1=NOP(N1,M1)
MM1=MM-2
IF(MM1 .LT. 1) MM1=5
JJ1=NOP(N2,MM1)
IF(J1 .EQ. JJ1) THEN
!
! We have the other match find nodes around the element
!
MID4=MM+1
JMID4=NOP(N2,MID4)
MID5=MM+3
IF(MM .EQ. 5) MID5=2
JMID5=NOP(N2,MID5)
M2=9-M-M1
MM2=9-MM-MM1
J2=NOP(N1,M2)
JJ2=NOP(N2,MM2)
NOP(N1,1)=J2
NOP(N1,2)=JMID3
NOP(N1,3)=J
NOP(N1,4)=JMID4
NOP(N1,5)=JJ2
NOP(N1,6)=JMID1
NOP(N2,1)=JJ2
NOP(N2,2)=JMID5
NOP(N2,3)=J1
NOP(N2,4)=JMID2
NOP(N2,5)=J2
NOP(N2,6)=JMID1
write(90,*) (nop(n1,i),i=1,6)
write(90,*) (nop(n2,i),i=1,6)
if(jmid1 .gt. 0) then
CORD(JMID1,1) = (CORD(J2,1)+CORD(JJ2,1))/2.0
CORD(JMID1,2) = (CORD(J2,2)+CORD(JJ2,2))/2.0
XUSR(JMID1) = CORD(JMID1,1)*TXSCAL - XS
YUSR(JMID1) = CORD(JMID1,2)*TXSCAL - YS
IF(NECON(JMID2,1) .EQ. N1) NECON(JMID2,1)=N2
IF(NECON(JMID2,2) .EQ. N1) NECON(JMID2,2)=N2
IF(NECON(JMID4,1) .EQ. N2) NECON(JMID4,1)=N1
IF(NECON(JMID4,2) .EQ. N2) NECON(JMID4,2)=N1
endif
GO TO 350
ENDIF
ENDIF
250 CONTINUE
300 END DO
350 CONTINUE
! CALL PLOTOT
RETURN
END