PROGRAM SETMBC2
C=======================================================================
C DATA GENERATING PROGRAM FOR FEM4Q.FOR
C DOMAIN : TLX BY TLY
C BOUNDARY CONDITIONS
C ON FACES OF -X : H1
C ON FACES OF -Y: Q=0.
C ON FACE OF +X: H2 AND SEEPAGE FACE
C ON FACE OF +Y: MOVING WITH Q=0.
C EIJI FUKUMORI NOVEMBER 28, 1995
C FEBRUARY 06, 2013
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER (ND=4,MXE=4000, MXN=5000,MXB=1000 )
PARAMETER ( TLX = 10., TLY = 5. )
PARAMETER ( INTEPT=2, MXMVG=100,MXGEACH=100 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDT(MXB),
* BVT(MXB), ITYPET(MXB),NMOVE(MXMVG), MOVEND(MXMVG,MXGEACH)
DIMENSION XBOTTOM(MXMVG), XTOP(MXMVG)
CHARACTER PROJECT*12,EXFILE*3
LOGICAL YES
C=======================================================================
DATA PROJECT / 'MOVINGBC.DAT' /
C=======================================================================
C NEY: NEMBER OF VERTICAL ELEMENTS (NUMBER OF NODES: NEY+1)
C NEX: NEMBER OF HORIZONTAL ELEMENTS (NUMBER OF NODES: NEX+1)
C=======================================================================
MAXIT=700
ERALLOW=0.0005D0
C
NEX = 30
NEY = 30
C
DX = TLX / NEX
DY = TLY / NEY
NDX=NEX+1
NDY=NEY+1
EXX = 1.D0
EXY = 0.D0
EYY = 1.D0
C
C---------------NDMAXH2 = # OF NODES IN Y DOWN-STREAM HEAD IS ASSIGNED
H1 = TLY
NDMAXH2 = 0.25D0*NEY+1
C=======================================================================
C ELEMENT CREATION
NE = 0
DO I = 1 , NEX
DO J = 1 , NEY
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = NDY*(I-1) + J
NODEX(NE,2) = NODEX(NE,1) + NDY
NODEX(NE,3) = NODEX(NE,2) + 1
NODEX(NE,4) = NODEX(NE,3) - NDY
END DO
END DO
C=======================================================================
C NODAL COORDINATE CREATION
C------XCOORDINATE VALUES ALONG BOTTOM AND TOP LINES
TLXTOP = TLX*0.7
DXTOP = TLXTOP/NEX
DO I = 1 , NDX
XBOTTOM(I) = DX*(I-1)
XTOP(I) = TLX/2-TLXTOP/2+DXTOP*(I-1)
END DO
C
NNODE = 0
DO I = 1 , NDX
DO J = 1 , NDY
NNODE = NNODE + 1
IF ( NNODE .GT. MXN ) STOP 'NNODE > MXN'
NODE = NDY*(I-1) + J
XCOORD(NODE) = F1(XBOTTOM(I),XTOP(I),NDY,J)
YCOORD(NODE) = DY*(J-1)
END DO
END DO
C=======================================================================
NBT = 0
C--------- FACE OF -X
DO J = 1 , NDY
NBT = NBT + 1
IBNDT(NBT) = J
ITYPET(NBT) = 1
BVT(NBT) = H1
END DO
C--------- FACE OF -Y
C NO FLOW; Q=0
C--------- FACE OF +X
C== HYDRAULIC HEAD = CONSTANT NODE
H2 = YCOORD(NDY*NEX+NDMAXH2)
DO J = 1 , NDMAXH2
NBT = NBT + 1
IBNDT (NBT) = NDY*NEX + J
ITYPET(NBT) = 1
BVT(NBT ) = H2
END DO
C== SEEPAGE FACE <---------------- ITYPET(NBT) = 3 -----------------
DO J = NDMAXH2+1, NEY
NBT = NBT + 1
IBNDT (NBT) = NDY*NEX + J
ITYPET(NBT) = 3
BVT(NBT ) = YCOORD (IBNDT(NBT))
END DO
C--------- FACE OF +Y
C NO FLOW; Q=0 AND FREE SURFACE
C
C---------- MOVING NODES (GROUPING)
C*******************************************************************
NGROUP = 0
DO I = 2 , NDX
II = I-1
NGROUP = NGROUP + 1
NMOVEI = 0
ISTART = 1
IF ( I .EQ. NDX ) ISTART = NDMAXH2
IF ( I .EQ. NDX-1 ) ISTART = NDMAXH2-2
IF ( I .EQ. NDX-2 ) ISTART = NDMAXH2-3
IF ( ISTART .LE. 0 ) ISTART = 1
DO J = ISTART, NDY
NMOVEI = NMOVEI + 1
JJ = J - ISTART+1
MOVEND(II,JJ) = NDY*(I-1)+J
END DO
NMOVE(II) = NMOVEI
END DO
C=======================================================================
WRITE (*,*) ' NUMBER OF ELEMENTS (NE) = ',NE
WRITE (*,*) ' NUMBER OF NODAL POINTS (NNODE) = ',NNODE
C=======================================================================
C DATA FILE INQUIRY
EXFILE = 'NEW'
INQUIRE ( FILE = PROJECT, EXIST = YES )
IF ( YES ) EXFILE='OLD'
C=======================================================================
C MAKING DATA FILES
C---------- 'PROJECT'.JNK
OPEN ( 1, FILE=PROJECT, STATUS = EXFILE )
C---------- MATERIAL PARAMETER
WRITE (1,*) EXX, EXY, EYY
C---------- CONVERGENCE PARAMETER
WRITE (1,*) MAXIT,ERALLOW
C---------- ELEMENT INFORMATION
WRITE (1,*) NE
DO IEL = 1 , NE
WRITE (1,*) IEL,(NODEX(IEL,J),J=1,ND)
END DO
C---------- NODAL COORDINATE INFOMATION
WRITE (1,*) NNODE
DO NODE = 1 , NNODE
WRITE (1,*) NODE,XCOORD(NODE),YCOORD(NODE)
END DO
C---------- BOUNDARY CONDITIONS
WRITE (1,*) NBT
DO I = 1 , NBT
WRITE (1,*) IBNDT(I), ITYPET(I), BVT(I)
END DO
C---------- MOVING NODE INFORMATION
WRITE (1,*) NGROUP
DO I = 1 , NGROUP
NNMOVE = NMOVE(I)
WRITE (1,*) NNMOVE, (MOVEND(I,J), J=1,NNMOVE)
END DO
CLOSE (1)
C=======================================================================
CALL BANDWID ( ND, MXE, NE, NODEX , NBW )
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
WRITE (1,*) ' PARAMETER ( ND=',ND,' )'
WRITE (1,*) ' PARAMETER ( INTEPT=',INTEPT,')'
WRITE (1,*) ' PARAMETER ( MXE=',NE,')'
WRITE (1,*) ' PARAMETER ( MXN=',NNODE,')'
WRITE (1,*) ' PARAMETER ( MXW=',NBW,')'
WRITE (1,*) ' PARAMETER ( MXB=',NBT,')'
WRITE (1,*) ' PARAMETER ( MXMVG=',NGROUP,')'
NMOVEMX = 0
DO I = 1 , NGROUP
NMOVEMX = MAX0 ( NMOVEMX , NMOVE(I) )
END DO
WRITE (1,*) ' PARAMETER ( MXMVEACH=',NMOVEMX,')'
CLOSE (1)
C=======================================================================
STOP "NORMAL TERMINATION"
END
C
C
SUBROUTINE BANDWID ( ND, MXE, NE, NODEX , NBW )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND)
NBW = 0
DO I = 1 , NE
DO J = 1 , ND-1
DO K = J+1 , ND
NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
END DO
END DO
END DO
NBW = NBW + 1
WRITE(*,*) 'HALF BANDWIDTH + 1 =', NBW
RETURN
END
C
C
FUNCTION F1(XBOTTOM,XTOP,NDY,J)
IMPLICIT REAL*8 ( A-H , O-Z )
T = DFLOAT(J-1)/DFLOAT(NDY-1)
F1 = (1.-T)*XBOTTOM + T*XTOP
RETURN
END