PROGRAM SET4QUARTERMODEL
C=======================================================================
C POISSON'S EQUATION FOR TORSION PROBLEM
C D2(PHI)/DXDX + D2(PHI)/DYDY +2*GM*THETA = 0
C 1/4 MODEL
C
C ^Y
C I
C I Φ=0
C Y=+TLY I-------<--------
C I I
C I I
C V IΦ=0
C DΦ/DN =0 I I
C I I
C I I X=+TLX
C 0*------->-------------------->X
C 0 DΦ/DN =0
C
C
C *=START AND END POINT FOR BOUNDARY INTEGRAL
C GM=SHEAR MODULUS
C THETA = RATE OF TWIST
C (PHI) = STRESS FUNCTION
C DATA GENERATING PROGRAM FOR TORSION4Q.FOR
C AREA INTEGRAL: FOUR-NODED ISOPARAMETRIC ELEMENT
C BEM BOUNDARY INTEGRAL: TWO-NODED LINEAR ELEMENT
C EIJI FUKUMORI DECEMBER 29, 1993, 04/NOV./2024
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=4,MXE=50000,MXN=52000 )
PARAMETER ( NDBEM=2,MXEBEM=50000,MXNBEM=50000 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN)
DIMENSION NODEXBEM(MXEBEM,NDBEM),XCOORDBEM(MXNBEM),
* YCOORDBEM(MXNBEM), IELTYPE(MXEBEM), BV(MXEBEM,NDBEM),
* XI(MXN), YI(MXN)
C=======================================================================
C NEX = NUMBER OF ELEMENTS IN X-COORDINATE
C NEY = NUMBER OF ELEMENTS IN Y-COORDINATE
C BASIC PARAMETERS
PI = 4.D0* DATAN( 1.D0)
GM = 26538461.5384615D0
THETA = 0.0000670055862315403D0
CCCCCCCCCCCCCCCCCC GM = 1.D0
CCCCCCCCCCCCCCCCCC THETA = 1.D0
TLX = 1.D0
TLY = 1.D0
NEX = 10
NEY = 10
C=======================================================================
C DOMAIN ELEMENT CREATION
DX = TLX / NEX
DY = TLY / NEY
NDX=NEX+1
NDY=NEY+1
NE = 0
DO I = 1 , NEY
DO J = 1 , NEX
NE = NE + 1
IF ( NE .GT. MXE ) STOP 'NE > MXE'
NODEX(NE,1) = NDX*(I-1) + J
NODEX(NE,2) = NODEX(NE,1) + 1
NODEX(NE,3) = NODEX(NE,2) + NDX
NODEX(NE,4) = NODEX(NE,1) + NDX
END DO
END DO
C=======================================================================
C NORDAL COORDINATES CREATION FOR DOMAIN INTEGRAL
NNODE = 0
DO I = 1 , NDY
DO J = 1 , NDX
NNODE = NNODE + 1
XCOORD(NNODE) = DX*(J-1)
YCOORD(NNODE) = DY*(I-1)
END DO
END DO
C=======================================================================
C LINEAR ELEMENT IS USED FOR BOUNDARY INTEGRATION
C NODAL NUMBERS CREATION FOR THE BOUNDARY ELEMENTS
NEBEM = 0
C--------- FACE OF -Y
DO I = 1 , NEX
NEBEM = NEBEM + 1
NODEXBEM(NEBEM,1) = I
NODEXBEM(NEBEM,2) = I+1
IELTYPE(NEBEM) = 2
BV(NEBEM,1) = 0.D0
BV(NEBEM,2) = 0.D0
END DO
C--------- FACE OF +X
DO I = 1 , NEY
NEBEM = NEBEM + 1
NODEXBEM(NEBEM,1) = NEX + I
NODEXBEM(NEBEM,2) = NEX + I+1
IELTYPE(NEBEM) = 1
BV(NEBEM,1) = 0.D0
BV(NEBEM,2) = 0.D0
END DO
C--------- FACE OF +Y
DO I = 1 , NEX
NEBEM = NEBEM + 1
NODEXBEM(NEBEM,1) = NEX + NEY + I
NODEXBEM(NEBEM,2) = NEX + NEY + I+1
IELTYPE(NEBEM) = 1
BV(NEBEM,1) = 0.D0
BV(NEBEM,2) = 0.D0
END DO
C--------- FACE OF -X
DO I = 1 , NEY
NEBEM = NEBEM + 1
NODEXBEM(NEBEM,1) = 2*NEX + NEY + I
IF ( I .EQ. NEY ) THEN
NODEXBEM(NEBEM,2) = 1
ELSE
NODEXBEM(NEBEM,2) = 2*NEX + NEY + I+1
END IF
IELTYPE(NEBEM) = 2
BV(NEBEM,1) = 0.D0
BV(NEBEM,2) = 0.D0
END DO
C=======================================================================
C NODAL COORDINATES CREATION FOR BOUNDARY INTEGRATION
NNODEBEM = 0
C--------- FACE OF -Y
DO I = 1 , NDX
NNODEBEM = NNODEBEM + 1
XCOORDBEM(NNODEBEM) = DX*(I-1)
YCOORDBEM(NNODEBEM) = 0.D0
END DO
C--------- FACE OF +X
DO I = 2 , NDY
NNODEBEM = NNODEBEM + 1
XCOORDBEM(NNODEBEM) = DX*NEX
YCOORDBEM(NNODEBEM) = DY*(I-1)
END DO
C--------- FACE OF +Y
DO I = 2 , NDX
NNODEBEM = NNODEBEM + 1
XCOORDBEM(NNODEBEM) = DX*NEX - DX*(I-1)
YCOORDBEM(NNODEBEM) = DY*NEY
END DO
C--------- FACE OF -X
DO I = 2 , NDY
NNODEBEM = NNODEBEM + 1
XCOORDBEM(NNODEBEM) = 0.D0
YCOORDBEM(NNODEBEM) = DY*NEY - DY*(I-1)
END DO
C=======================================================================
C------- INTERNAL POINTS: STRESS FUNCTION VALUE TO BE EVALUATED --------
NIPT = NNODE
DO I = 1 , NIPT
XI(I) = XCOORD(I)
YI(I) = YCOORD(I)
END DO
C=======================================================================
C MAKING DATA FILES FOR DOMAIN INTEGRAL
OPEN ( 1, FILE='DOMAIN.DAT', STATUS = 'UNKNOWN' )
WRITE(1,*) NE
C------ ELEMENTS CONECTIVITY INFORMATION
DO I = 1 , NE
WRITE (1,*) I,(NODEX(I,J),J=1,ND)
END DO
C------ NODAL COORDINATES
WRITE(1,*) NNODE
DO I = 1 , NNODE
WRITE(1,*) I, XCOORD(I), YCOORD(I)
END DO
CLOSE (1)
C=======================================================================
WRITE(*,*) 'NE(FEM)=',NE
WRITE(*,*) 'NNODE(FEM)=',NNODE
C=======================================================================
C MAKING DATA FILES FOR BOUNDARY INTEGRAL
OPEN ( 1, FILE='BEM1.DAT', STATUS = 'UNKNOWN' )
WRITE(1,*) GM, THETA
WRITE(1,*) NEBEM
DO I = 1 , NEBEM
WRITE (1,*)I,(NODEXBEM(I,J),J=1,NDBEM),IELTYPE(I),
* (BV(I,J),J=1,NDBEM)
END DO
DO I = 1 , NEBEM
WRITE(1,*) I, XCOORDBEM(I), YCOORDBEM(I)
END DO
WRITE(1,*) NIPT
DO I = 1 , NIPT
WRITE(1,*) I, XI(I), YI(I)
END DO
CLOSE (1)
C=======================================================================
WRITE(*,*) 'NE(BEM)=',NEBEM
WRITE(*,*) 'NNODE(BEM)=',NNODEBEM
C=======================================================================
C DOMAIN Element Mapping
C=======================================================================
OPEN ( 1, FILE='ELEMENT4.OUT', STATUS = 'UNKNOWN' )
DO IEL = 1 , NE
DO J = 1 , ND
NODE = NODEX(IEL,J)
WRITE(1,*) XCOORD(NODE),YCOORD(NODE)
END DO
NODE = NODEX(IEL,1)
WRITE(1,*) XCOORD(NODE),YCOORD(NODE)
WRITE(1,*)
END DO
CLOSE (1)
C=======================================================================
C Boundary element mapping
OPEN ( 1, FILE='ELEMENT2.OUT', STATUS = 'UNKNOWN' )
DO IEL = 1 , NEBEM
DO J = 1 , NDBEM
NODE = NODEXBEM(IEL,J)
WRITE(1,*) XCOORDBEM(NODE),YCOORDBEM(NODE)
END DO
WRITE(1,*)
END DO
CLOSE (1)
C=======================================================================
C------ CREATION OF PARAMETER FILE TO BE USED IN INCLUDE STATEMENT
OPEN ( 1, FILE='PARAM.DAT', STATUS='UNKNOWN' )
C------------ BOUNDARY ANALYSIS DIMENSIONS
INTEPT = 2
WRITE (1,*) ' PARAMETER ( ND=',NDBEM,')'
WRITE (1,*) ' PARAMETER ( INTEPT=',INTEPT,')'
WRITE (1,*) ' PARAMETER ( MXE=',NEBEM,')'
WRITE (1,*) ' PARAMETER ( MXN=',NEBEM,')'
WRITE (1,*) ' PARAMETER ( MXI=',NNODE,')'
C------------ DOMAIN INTERAL ANALYSIS DIMENSIONS
C------------ DOMAIN INTERAL ANALYSIS DIMENSIONS
INTEPTFEM = 2
INTEPTM = INTEPTFEM + 1
WRITE (1,*) ' PARAMETER ( NDFEM=',ND,')'
WRITE (1,*) ' PARAMETER ( INTEPTFEM=',INTEPTFEM,')'
WRITE (1,*) ' PARAMETER ( MXEFEM=',NE,')'
WRITE (1,*) ' PARAMETER ( MXNFEM=',NNODE,')'
WRITE (1,*) ' PARAMETER ( INTEPTM=',INTEPTM,')'
CLOSE (1)
C=======================================================================
STOP "NORMAL TERMINATION"
END