PROGRAM SET4WARPING
C=======================================================================
C LAPLACE EQUATION FOR TORSION PROBLEM
C WARPING
C D2(PSI)/DXDX + D2(PSI)/DYDY = 0
C FULL MODEL
C
C ^Y
C I I
C I Dψ/DN=-X IY=+TLY/2
C Y=+TLY I-------<-----------------------
C I I I
C I I I
C V I I
C Dψ/DN =Y I I I Dψ/DN =Y
C I I I
C I 0I I X=+TLX/2
C X=-TLX/2 0I--------------0------------------->X
C I I I
C I I I
C V I
C Dψ/DN =Y I I I Dψ/DN =Y
C I I I
C I I I
C *------->------------------------
C Dψ/DN=-X Y=-TLY/2
C
C *=START AND END POINT FOR BOUNDARY INTEGRAL
C GM=SHEAR MODULUS
C THETA = RATE OF TWIST
C (PSI) = WARPING FUNCTION
C DATA GENERATING PROGRAM FOR BEM8LINQ.FOR
C BEM BOUNDARY INTEGRAL: TWO-NODED LINEAR ELEMENT
C EIJI FUKUMORI DECEMBER 29, 1993, 01/DEC/2024
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=2, MXE=5000, MXN=5000 )
C=======================================================================
DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN),
* IELTYPE(MXE), BV(MXE,ND), XI(MXN), YI(MXN)
C=======================================================================
F(X)=-(0.383141762452107)*X*(X-1.9)*(X+1.9)
CCCCCCCCCCCCCCC F(X)=-(1./3.)*X*(X-2.)*(X+2.)
CCCCCCCCCCCCCCC F(X)=-(0.445428571428571)*X*(X-1.8)*(X+1.8)
C=======================================================================
C NEX = NUMBER OF ELEMENTS IN X-COORDINATE
C NEY = NUMBER OF ELEMENTS IN Y-COORDINATE
C=======================================================================
C BASIC PARAMETERS
PI = 4.D0* DATAN( 1.D0 )
GM = 26538461.5384615D0
THETA = 0.0000670055862315403D0
C=======================================================================
TLX = 2.D0
TLY = 2.D0
NEX = 20
NEY = 20
C=======================================================================
DX = TLX / NEX
DY = TLY / NEY
NDX=NEX+1
NDY=NEY+1
EPS = 0.001D0
HX = TLX/2.D0
HY = TLY/2.D0
C=======================================================================
C LINEAR ELEMENT FOR BOUNDARY INTEGRATION
C NODAL COORDINATES CREATION FOR BOUNDARY INTEGRATION
NNODE = 0
C--------- FACE OF -Y
DO I = 1 , NDX
NNODE = NNODE + 1
YCOORD(NNODE) = -TLY/2.D0
XCOORD(NNODE) = DX*(I-1) - TLX/2.D0
XCOORD(NNODE) = HX*F(XCOORD(NNODE)/HX)
IF ( I .EQ. 2 ) XCOORD(NNODE) = EPS*DX - TLX/2.D0
IF ( I .EQ. NDX-1 ) XCOORD(NNODE) = -EPS*DX + TLX/2.D0
END DO
C--------- FACE OF +X
DO I = 2 , NDY
NNODE = NNODE + 1
XCOORD(NNODE) = TLX/2.D0
YCOORD(NNODE) = DY*(I-1)-TLY/2.D0
YCOORD(NNODE) = HY*F(YCOORD(NNODE)/HY)
IF ( I .EQ. 2 ) YCOORD(NNODE) = EPS*DY - TLY/2.D0
IF ( I .EQ. NDY-1 ) YCOORD(NNODE) = -EPS*DY + TLY/2.D0
END DO
C--------- FACE OF +Y
DO I = 2 , NDX
NNODE = NNODE + 1
YCOORD(NNODE) = TLY/2.D0
XCOORD(NNODE) = TLX/2.D0 - DX*(I-1)
XCOORD(NNODE) = HX*F(XCOORD(NNODE)/HX)
IF ( I .EQ. 2 ) XCOORD(NNODE) = TLX/2.D0 - EPS*DX
IF ( I .EQ. NDX-1 ) XCOORD(NNODE) = EPS*DX - TLX/2.D0
END DO
C--------- FACE OF -X
DO I = 2 , NDY
NNODE = NNODE + 1
XCOORD(NNODE) = -TLX/2.D0
YCOORD(NNODE) = TLY/2.D0 - DY*(I-1)
YCOORD(NNODE) = HY*F(YCOORD(NNODE)/HY)
IF ( I .EQ. 2 ) YCOORD(NNODE) = -EPS*DY + TLY/2.D0
IF ( I .EQ. NDY-1 ) YCOORD(NNODE) = EPS*DY - TLY/2.D0
END DO
C=======================================================================
C ELEMENT NODES AND BOUNDARY CONDITIONS
NE = 0
C--------- FACE OF -Y BOTH END ELEMENTS: DIRICHILET, OTHER: NEUMANN
DO I = 1 , NEX
NE = NE + 1
NODEX(NE,1) = I
NODEX(NE,2) = I+1
IF ( (I .EQ. 1) .OR. (I .EQ. NEX) ) THEN
IELTYPE(NE) = 1
BV(NE,1) = 0.D0
BV(NE,2) = 0.D0
ELSE
IELTYPE(NE) = 2
BV(NE,1) = -XCOORD(NODEX(NE,1))
BV(NE,2) = -XCOORD(NODEX(NE,2))
END IF
END DO
C--------- FACE OF +X
DO I = 1 , NEY
NE = NE + 1
NODEX(NE,1) = NEX + I
NODEX(NE,2) = NEX + I+1
IF ( (I .EQ. 1) .OR. (I .EQ. NEY) ) THEN
IELTYPE(NE) = 1
BV(NE,1) = 0.D0
BV(NE,2) = 0.D0
ELSE
IELTYPE(NE) = 2
BV(NE,1) = -YCOORD(NODEX(NE,1))
BV(NE,2) = -YCOORD(NODEX(NE,2))
END IF
END DO
C--------- FACE OF +Y
DO I = 1 , NEX
NE = NE + 1
NODEX(NE,1) = NEX + NEY + I
NODEX(NE,2) = NEX + NEY + I+1
IF ( (I .EQ. 1) .OR. (I .EQ. NEX) ) THEN
IELTYPE(NE) = 1
BV(NE,1) = 0.D0
BV(NE,2) = 0.D0
ELSE
IELTYPE(NE) = 2
BV(NE,1) = XCOORD(NODEX(NE,1))
BV(NE,2) = XCOORD(NODEX(NE,2))
END IF
END DO
C--------- FACE OF -X
DO I = 1 , NEY
NE = NE + 1
NODEX(NE,1) = 2*NEX + NEY + I
IF ( I .EQ. NEY ) THEN
NODEX(NE,2) = 1
ELSE
NODEX(NE,2) = 2*NEX + NEY + I+1
END IF
IF ( (I .EQ. 1) .OR. (I .EQ. NEY) ) THEN
IELTYPE(NE) = 1
BV(NE,1) = 0.D0
BV(NE,2) = 0.D0
ELSE
IELTYPE(NE) = 2
BV(NE,1) = YCOORD(NODEX(NE,1))
BV(NE,2) = YCOORD(NODEX(NE,2))
END IF
END DO
C=======================================================================
C------- INTERNAL POINTS: STRESS FUNCTION VALUE TO BE EVALUATED --------
NIPT = 0
C=======================================================================
C=======================================================================
C=======================================================================
C MAKING DATA FILES FOR BOUNDARY INTEGRAL
OPEN ( 1, FILE='BEM1.DAT', STATUS = 'UNKNOWN' )
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
WRITE(1,*) GM, THETA
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
WRITE(1,*) NE
DO I = 1 , NE
WRITE (1,*)I,(NODEX(I,J),J=1,ND),IELTYPE(I),
* (BV(I,J),J=1,ND)
END DO
DO I = 1 , NE
WRITE(1,*) I, XCOORD(I), YCOORD(I)
END DO
WRITE(1,*) NIPT
IF ( NIPT .GT. 0 ) THEN
DO I = 1 , NIPT
WRITE(1,*) I, XI(I), YI(I)
END DO
END IF
CLOSE (1)
C=======================================================================
WRITE(*,*) 'NE(BEM)=',NE
WRITE(*,*) 'NNODE(BEM)=',NNODE
C=======================================================================
C=======================================================================
C Boundary element mapping
OPEN ( 1, FILE='ELEMENT2.OUT', STATUS = 'UNKNOWN' )
DO IEL = 1 , NE
DO J = 1 , ND
NODE = NODEX(IEL,J)
WRITE(1,*) XCOORD(NODE),YCOORD(NODE)
END DO
WRITE(1,*)
END DO
CLOSE (1)
C=======================================================================
CCCCC PARAMETER (MXE=300, MXN=MXE, MXI=2000,INTEPT=2, ND=2)
CCCCC PARAMETER (MXEFEM=10000, MXNFEM=10000, INTEPTFEM=2, NDFEM=4)
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=',ND,')'
WRITE (1,*) ' PARAMETER ( INTEPT=',INTEPT,')'
WRITE (1,*) ' PARAMETER ( MXE=',NE,')'
WRITE (1,*) ' PARAMETER ( MXN=',NE,')'
WRITE (1,*) ' PARAMETER ( MXI=',NNODE,')'
CLOSE (1)
C=======================================================================
STOP "NORMAL TERMINATION"
END