PROGRAM BEM8LINQTORSIONM
C======================================================================
C BOUNDARY ELEMENT METHOD
C APPLIED TO
C POISSON'S EQUATION FOR SOLVING TORSION PROBLEM
C UNKNOWN VARIABLE: STRESS FUCTION (PHI)
C EQUATION: D2(PHI)/DXDX + D2(PHI)/DYDY + 2*G*THETA = 0
C G = SHEAR MODULUS
C THETA = RATE OF TWISTE
C ELEMENT TYPE: LINEAR ELEMENT
C== ALONG WITH DIRECT FORMATION OF [A]{X}={RHS} ==
C PROGRAMMED BY EIJI FUKUMORI // SUNY AT BUFFALO // 1984 SPRING
C revised 07/NOV/2024
C======================================================================
INCLUDE 'PARAM.DAT'
IMPLICIT REAL*8 ( A-H , O-Z )
CCCCCCCCCCC PARAMETER (MXE=300, MXN=MXE, MXI=2000,INTEPT=2, ND=2)
CCCC PARAMETER (MXEFEM=10000, MXNFEM=10000, INTEPTFEM=2, NDFEM=4)
CCCC PARAMETER ( INTEPTM=2 )
C======================== DIMENSION FOR BEM ===========================
DIMENSION XE(ND),YE(ND),GE(ND),FE(ND),
* SAI(INTEPT),W(INTEPT), NODEX(MXE,ND),IELTYPE(MXE),BV(MXE,ND),
* A(MXN,MXN), C(MXN), X(MXN),Y(MXN),QN(MXN),H(MXN),RHS(MXN),
* NDTYPE(MXN), XI(MXI),YI(MXI),HI(MXI),CI(MXI),ICHECK(MXI)
C======================== DIMENSION FOR FEM ===========================
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM), SAIFEM(INTEPTFEM),WFEM(INTEPTFEM)
DIMENSION BPPFEM(2,NDFEM,INTEPTFEM,INTEPTFEM),
* SFFEM(NDFEM,INTEPTFEM,INTEPTFEM), SS(NDFEM,NDFEM),
* XDUMMY(NDFEM), YDUMMY(NDFEM)
C======================== DIMENSION FOR MOMENT ========================
DIMENSION SAIM(INTEPTM),WM(INTEPTM)
DIMENSION BPPM(2,NDFEM,INTEPTM,INTEPTM),
* SFM(NDFEM,INTEPTM,INTEPTM)
C========================= INITIAL SET-UPS ============================
C
PI = 4.D0 * DATAN( 1.D0)
C1 = - 1.D0/ ( 2.D0 * PI )
CALL GRULE ( INTEPT, SAI, W )
CALL GRULE ( INTEPTFEM, SAIFEM, WFEM )
CALL DERIV ( NDFEM, INTEPTFEM, XDUMMY, YDUMMY, SAIFEM, BPPFEM )
CALL SHAPEF( NDFEM, INTEPTFEM, XDUMMY, SAIFEM, SFFEM )
C=============== INETEGARTION FOR MOMENT CALCULATION ==================
CALL GRULE ( INTEPTM, SAIM, WM )
CALL DERIV ( NDFEM, INTEPTM, XDUMMY, YDUMMY, SAIM, BPPM )
CALL SHAPEF( NDFEM, INTEPTM, XDUMMY, SAIM, SFM )
C========================= READING IN DATA ============================
C BOUNDARY ELEMENT DATA
CALL INPUT (GM,THETA,ND,MXE,MXI,NE,NODEX,IELTYPE,BV,X,Y,
* NIPT,XI,YI)
C ...TYPE(I)=1 ---> STRESS FUNCTION(PHI) PRESCRIBED
C ...TYPE(I)=2 ---> D(PHI)/DN PRESCRIBED
NNODE = NE
C======================================================================
C FEM DOMAIN DATA
CALL INPUTFEM ( NDFEM,MXEFEM,MXNFEM,NEFEM,NNODEFEM,
* NODEXFEM,XCOORDFEM,YCOORDFEM )
C======================================================================
C================= ANALYSIS ON NODAL CONDITION ========================
C
CALL NDANALYS (ND,MXE,MXN,NE,IELTYPE,NODEX,NDTYPE )
C====================== DIRICHLET B.C. TO H(I) ========================
C
CALL BVARRANG ( ND, MXE,MXN, NE,NODEX,IELTYPE, BV, H )
C=============== FORMATION OF MATRIX A AND VECTOR RHS ==============
C
CALL MTXFORM (ND,MXE,MXN,C1,INTEPT,SAI,W,NE,NODEX,X,Y,
* NDTYPE, XE,YE, GE,FE, C ,IELTYPE,BV,RHS,H, A,
* INTEPTFEM, MXEFEM,MXNFEM,NDFEM, NEFEM,XCOORDFEM, YCOORDFEM,
* GM, THETA, NODEXFEM,BPPFEM, SFFEM,WFEM )
C===================== READY TO SOLVE A . X = C ===================
C
CALL SYSTEM ( MXE , NNODE, A , RHS )
C====================== SORTING SOLUTION ==============================
C
CALL SORTSOLN (ND,MXE,MXN,NE,IELTYPE,NODEX,RHS,NDTYPE,
* H,QN,BV )
C======================= INTERNAL POINTS ============================
C
CALL CHECK ( MXN, MXI, NIPT, NE, XI, YI, X, Y, ICHECK )
C
CALL DOMAIN ( INTEPT,ND,MXE,MXN,MXI,C1,NIPT,NE,
* GE,FE, SAI,W,XI,YI,NODEX,X,Y,H,BV,HI, CI,XE,YE,
* INTEPTFEM, MXEFEM,MXNFEM,NDFEM, NEFEM,XCOORDFEM, YCOORDFEM,
* GM, THETA, NODEXFEM,BPPFEM, SFFEM,WFEM,ICHECK , C )
C
C==================== GRAPHIC FILES MAKING ============================
CALL GRAPHIC ( NDFEM,MXEFEM,MXI,NIPT,XI,YI,HI,GM,THETA,NEFEM,
* NODEXFEM )
C
C================= PRINTING RESULTS =================================
CALL ECHOSOLN( MXE,MXN,MXI,NE,ND,NIPT,NODEX,X,Y,IELTYPE,
* H, BV, XI, YI, HI,CI, C )
CALL MOMENT ( MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM,NEFEM,
* SFFEM,XCOORDFEM,YCOORDFEM,NODEXFEM, HI,TORQUE )
C
C
CALL MOMENT2 ( INTEPT,ND,MXE,MXN,C1,NE,
* GE,FE, SAI,W,NODEX,X,Y,H,BV,XE,YE,
* GM, THETA,MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM,
* NEFEM,SFFEM, XCOORDFEM,YCOORDFEM,NODEXFEM,
* INTEPTM, SAIM, WM, BPPM, SFM, TORQUE2 )
C
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE NDANALYS (ND,MXE,MXN,NE,IELTYPE,NODEX,NDTYPE )
DIMENSION IELTYPE(MXE), NDTYPE(MXN), NODEX(MXE,ND)
C------- THIS EVALUATES NODAL BOUNDARY CONDITION -------
C NDTYPE(NODE)=1 IMPLIES H IS KNOWN AT THE NODE.
C NDTYPE(NODE)=2 IMPLIES QN IS KNOWN AT THE NODE.
NNODE = NE
DO I = 1 , NNODE
NDTYPE(I) = 0
END DO
C------- CORNER(NODAL) BOUNDARY CONDITION EVALUATION -------
DO IEL = 1 , NE
DO J = 1 , ND
NDTYPE(NODEX(IEL,J)) = NDTYPE(NODEX(IEL,J)) + IELTYPE(IEL)
END DO
END DO
C-------
DO I = 1 , NNODE
IF ( NDTYPE(I) .EQ. 2 ) J = 1
IF ( NDTYPE(I) .EQ. 3 ) J = 3
IF ( NDTYPE(I) .EQ. 4 ) J = 2
NDTYPE(I) = J
END DO
RETURN
END
C
C
SUBROUTINE SORTSOLN (ND,MXE,MXN,NE,IELTYPE,NODEX,RHS,NDTYPE,
* H,QN,BV )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION IELTYPE(MXE),RHS(MXN),QN(MXN),H(MXN),NDTYPE(MXN),
* BV(MXE,ND),NODEX(MXE,ND)
NNODE = NE
DO I = 1 , NNODE
IF ( NDTYPE(I) .EQ. 2 ) THEN
H (I) = RHS(I)
ELSE
QN(I) = RHS(I)
END IF
END DO
C------- REDISTRIBUTION OF QN(I) INTO BV(I,J)
DO IEL = 1 , NE
IF ( IELTYPE(IEL) .EQ. 1 ) THEN
DO J = 1 , ND
NODE = NODEX(IEL,J)
BV(IEL,J) = QN(NODE)
END DO
END IF
END DO
RETURN
END
C
C
SUBROUTINE BVARRANG ( ND, MXE,MXN, NE,NODEX,IELTYPE, BV, H )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION IELTYPE(MXE),BV(MXE,ND),H(MXN),NODEX(MXE,ND)
C------ CONVERT STRESS FUNCTION VALUE ALONG ELEMENT TO THE NODEL
DO I = 1 , NE
DO J = 1 , ND
IF ( IELTYPE(I) .EQ. 1 ) THEN
H(NODEX(I,J)) = BV(I,J)
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE MTXFORM (ND,MXE,MXN,C1,INTEPT,SAI,W,NE,NODEX,X,Y,
* NDTYPE, XE,YE, GE,FE, C ,IELTYPE,BV,RHS,H, A,
* INTEPTFEM, MXEFEM,MXNFEM,NDFEM, NEFEM,XCOORDFEM, YCOORDFEM,
* GM, THETA, NODEXFEM,BPPFEM, SFFEM,WFEM )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION C(MXN),X(MXN),Y(MXN),GE(ND),FE(ND),
* XE(ND),YE(ND),SAI(INTEPT),W(INTEPT) , NODEX(MXE,ND),
* IELTYPE(MXE),BV(MXE,ND), H(MXN),RHS(MXN),A(MXN,MXN),
* NDTYPE(MXN)
C============================= FEM DIMENSION ==========================
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM), WFEM(INTEPTFEM),
* BPPFEM(2,NDFEM,INTEPTFEM,INTEPTFEM),
* SFFEM(NDFEM,INTEPTFEM,INTEPTFEM)
C======================================================================
C=============== FORMATION OF MATRIX G,F AND VECTOR C
C---------- CLEAR MATRIX A(I,J) AND RHS(I) AND C(I)
NNODE = NE
DO I = 1 , NNODE
C(I) = 0.D0
RHS(I) = 0.D0
DO J = 1 , NNODE
A(I,J) = 0.D0
END DO
END DO
C ....... (XP,YP) = COORDINATES OF OBSERVATION POINT.
DO ISOURCE = 1 , NNODE
XP = X(ISOURCE)
YP = Y(ISOURCE)
DO IEL = 1 , NE
DO I = 1 , ND
XE(I) = X(NODEX(IEL,I))
YE(I) = Y(NODEX(IEL,I))
END DO
C------------------ INTEGRAL ON AN ELEMENT
NDDIFF1 = ISOURCE-NODEX(IEL,1)
NDDIFF2 = ISOURCE-NODEX(IEL,2)
IF ( NDDIFF1*NDDIFF2 .EQ. 0 ) THEN
CALL FINE ( ND,NDDIFF2, C1, XE, YE, GE, FE )
ELSE
CALL INTE ( INTEPT,ND,XP,YP,C1,XE,YE,SAI,W, GE, FE )
END IF
C------------------ MATRIX FORMATION
C.............NON-FREE TERMS
DO L = 1 , ND
ICURREN = NODEX(IEL,L)
IF (IELTYPE(IEL) .EQ. 1 ) THEN
RHS(ISOURCE) = RHS(ISOURCE) + FE(L) * H(ICURREN)
A(ISOURCE,ICURREN) = A(ISOURCE,ICURREN) + GE(L)
END IF
IF (IELTYPE(IEL) .EQ. 2 ) THEN
RHS(ISOURCE) = RHS(ISOURCE) - GE(L) * BV(IEL,L)
IF ( NDTYPE(ICURREN) .EQ. 2 ) THEN
A(ISOURCE,ICURREN) = A(ISOURCE,ICURREN) - FE(L)
ELSE
RHS(ISOURCE) = RHS(ISOURCE) + FE(L) * H(ICURREN)
END IF
END IF
C(ISOURCE) = C(ISOURCE) + FE(L)
END DO
END DO
C.............FREE TERM
IF (NDTYPE(ISOURCE) .EQ. 2 ) THEN
A(ISOURCE,ISOURCE) = A(ISOURCE,ISOURCE) + C(ISOURCE)
ELSE
RHS(ISOURCE) = RHS(ISOURCE) - C(ISOURCE) * H(ISOURCE)
END IF
C.............. BELOW DO ISOURCE = 1 , NNODE...........................
CALL GDM ( MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM, NEFEM,
* SFFEM,XCOORDFEM,YCOORDFEM,NODEXFEM,GM,THETA,C1, XP,YP, AJ )
RHS(ISOURCE) = RHS(ISOURCE) + AJ
END DO
RETURN
END
C
C=======================================================================
C=======================================================================
SUBROUTINE INPUT (GM,THETA,ND,MXE,MXI,NE,NODEX,IELTYPE,BV,X,Y,
* NIPT,XI,YI)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXE),Y(MXE),IELTYPE(MXE),BV(MXE,ND),NODEX(MXE,ND),
* XI(MXI),YI(MXI)
OPEN ( 1, FILE='BEM1.DAT', STATUS='UNKNOWN' )
READ (1,*) GM, THETA
C WRITE (*,*) GM, THETA
READ (1,*) NE
C WRITE (1,*) NE
DO IEL = 1 , NE
READ (1,*) I,(NODEX(I,J),J=1,ND),IELTYPE(I), (BV(I,J),J=1,ND)
END DO
NNODE = NE
DO I = 1 , NNODE
READ (1,*) NODE, X(NODE), Y(NODE)
END DO
READ (1,*) NIPT
IF ( NIPT .GE. 1 ) THEN
DO J = 1 , NIPT
READ (1,*) I, XI(I), YI(I)
END DO
END IF
CLOSE (1)
RETURN
END
C=======================================================================
C=======================================================================
SUBROUTINE INPUTFEM ( NDFEM,MXEFEM,MXNFEM,NEFEM,NNODEFEM,
* NODEXFEM,XCOORDFEM,YCOORDFEM )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM)
IR = 1
OPEN ( IR, FILE='DOMAIN.DAT', STATUS='UNKNOWN' )
READ (IR,*) NEFEM
DO I = 1 , NEFEM
READ (IR,*) IEL,( NODEXFEM(IEL,J),J=1,NDFEM )
END DO
READ (IR,*) NNODEFEM
DO I = 1 , NNODEFEM
READ (IR,*) NODE,XCOORDFEM(NODE),YCOORDFEM(NODE)
END DO
CLOSE (IR)
RETURN
END
C
C
SUBROUTINE FINE ( ND,NDDIFF2, C1, XE, YE, GE, FE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XE(ND),YE(ND), GE(ND),FE(ND)
FE(1) = 0.D0
FE(2) = 0.D0
DX = XE(2) - XE(1)
DY = YE(2) - YE(1)
DS = DSQRT ( DX*DX + DY*DY )
GE(1) = C1 * (DS/2.D0) *( DLOG(DS) - 1.5D0 )
GE(2) = C1 * (DS/2.D0) *( DLOG(DS) - 0.5D0 )
IF ( NDDIFF2 .EQ. 0 ) THEN
TEMP = GE(1)
GE(1) = GE(2)
GE(2) = TEMP
END IF
RETURN
END
C
C
SUBROUTINE INTE (INTEPT, ND, XP, YP, C1, XE, YE, SAI,W, GE,FE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XE(ND),YE(ND), SAI(INTEPT),W(INTEPT), GE(ND),FE(ND)
DO I = 1 , ND
GE(I) = 0.D0
FE(I) = 0.D0
END DO
DX = XE(2) - XE(1)
DY = YE(2) - YE(1)
DS = DSQRT ( DX*DX + DY*DY )
DETJ = DS/2.D0
XM = ( XE(2) + XE(1) ) /2.D0
YM = ( YE(2) + YE(1) ) /2.D0
C--------GAUSS INTEGRATION
DO IGAUSS = 1 , INTEPT
XGAUSS = DX/2.D0*SAI(IGAUSS) + XM
YGAUSS = DY/2.D0*SAI(IGAUSS) + YM
RX = XGAUSS - XP
RY = YGAUSS - YP
R = DSQRT ( RX*RX + RY*RY )
SF1 = 0.5D0 * ( 1.D0 - SAI(IGAUSS) )
SF2 = 0.5D0 * ( 1.D0 + SAI(IGAUSS) )
C-------- INTEGRATION OF G(R)
TEMP = DLOG(R) * W(IGAUSS)
GE(1) = GE(1) + TEMP * SF1
GE(2) = GE(2) + TEMP * SF2
C-------- INTEGRATION OF F(R)
TEMP = (RX*DY-RY*DX) / (R*R) * W(IGAUSS)
FE(1) = FE(1) + TEMP * SF1
FE(2) = FE(2) + TEMP * SF2
END DO
GE(1) = C1 * DETJ * GE(1)
GE(2) = C1 * DETJ * GE(2)
FE(1) = -C1 * DETJ /DS * FE(1)
FE(2) = -C1 * DETJ /DS * FE(2)
RETURN
END
C
C
SUBROUTINE DOMAIN ( INTEPT,ND,MXE,MXN,MXI,C1,NIPT,NE,
* GE,FE, SAI,W,XI,YI,NODEX,X,Y,H,BV,HI, CI,XE,YE,
* INTEPTFEM, MXEFEM,MXNFEM,NDFEM, NEFEM,XCOORDFEM, YCOORDFEM,
* GM, THETA, NODEXFEM,BPPFEM, SFFEM,WFEM,ICHECK, C )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
* XI(MXI),YI(MXI), X(MXN),Y(MXN),H(MXN), C(MXN), BV(MXE,ND),
* XE(ND),YE(ND),HI(MXI),CI(MXI),ICHECK(MXI),
* GE(ND),FE(ND)
C============================= FEM DIMENSION ==========================
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM), WFEM(INTEPTFEM),
* BPPFEM(2,NDFEM,INTEPTFEM,INTEPTFEM),
* SFFEM(NDFEM,INTEPTFEM,INTEPTFEM)
C======================================================================
IF ( NIPT .EQ. 0 ) RETURN
C----------------------------------------------------------------------
DO INSIDE = 1 , NIPT
IF ( ICHECK(INSIDE) .EQ. 0 ) THEN
XP = XI(INSIDE)
YP = YI(INSIDE)
SUM = 0.D0
CP = 0.D0
C======================================================================
DO IEL = 1 , NE
DO I = 1 , ND
XE(I) = X( NODEX(IEL,I) )
YE(I) = Y( NODEX(IEL,I) )
END DO
CALL INTE ( INTEPT, ND, XP, YP, C1, XE, YE, SAI,W,GE,FE )
DO J = 1 , ND
CP = CP + FE(J)
SUM = SUM + FE(J)*H(NODEX(IEL,J)) - GE(J)*BV(IEL,J)
END DO
END DO
C======================================================================
CALL GDM ( MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM, NEFEM,
* SFFEM,XCOORDFEM,YCOORDFEM,NODEXFEM,GM,THETA,C1, XP,YP, AJ )
CI(INSIDE) = CP
HI(INSIDE) = SUM + AJ
ELSE
HI(INSIDE) = H( ICHECK(INSIDE) )
CI(INSIDE) = C( ICHECK(INSIDE) )
END IF
END DO
C----------------------------------------------------------------------
RETURN
END
C
C
SUBROUTINE ECHOSOLN( MXE,MXN,MXI,NE,ND,NIPT,NODEX,X,Y,IELTYPE,
* H, BV, XI, YI, HI,CI, C )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN), Y(MXN), XI(MXI), YI(MXI), HI(MXI),CI(MXI),
* NODEX(MXE,ND), IELTYPE(MXE), H(MXN), BV(MXE,ND), C(MXN)
CHARACTER*9 BC
OPEN ( 1, FILE='BEM.SOL', STATUS='UNKNOWN' )
C---------INPUT COORDINATES AND B.C.
WRITE (1,*)"ELEMENT# NODE(1) NODE(2) B.C. B-VALUE(1) B-VALUE(2)"
DO IEL = 1 , NE
IF ( IELTYPE(IEL) .EQ. 1 ) THEN
BV1 = H(NODEX(IEL,1))
BV2 = H(NODEX(IEL,2))
BC ='DIRICHLET'
END IF
IF ( IELTYPE(IEL) .EQ. 2 ) THEN
BV1 = BV(IEL,1)
BV2 = BV(IEL,2)
BC ='NUEMANN'
END IF
WRITE (1,*) IEL, NODEX(IEL,1), NODEX(IEL,2), BC, BV1, BV2
END DO
C
WRITE(1,*)
WRITE(1,*) "NODAL# XCOORD YCOORD STRESS-FUNCTION FREE-TERM"
NNODE = NE
DO I = 1 , NNODE
WRITE (1,*) I, X(I) , Y(I) , H(I) , C(I)
END DO
I = 1
WRITE (1,*) I, X(I) , Y(I) , H(I) , C(I)
C
C------------- D(PHI)/DN'S ARE STORED IN BV(IEL,1) AND BV(IEL,2).------
WRITE(1,*)
WRITE(1,*) "ELEMENT# X1 Y1 X2 Y2 D(PHI)/DN(1) D(PHI)/DN(2)"
DO IEL = 1 , NE
WRITE(1,*) IEL, X(NODEX(IEL,1)), Y(NODEX(IEL,1)),
* X(NODEX(IEL,2)), Y(NODEX(IEL,2)),
* BV(IEL,1), BV(IEL,2)
END DO
CLOSE (1)
C----------- MAKING FILE FOR SOLUTION ON INTERNAL POINTS --------
IF ( NIPT .NE. 0 ) THEN
OPEN ( 3, FILE='INTERNAL.SOL', STATUS='UNKNOWN' )
WRITE (3,*) "INTERNALPOINT# XCRD YCRD STRESS-FUNCTION FREE-TERM"
DO I = 1 , NIPT
WRITE(3,*) I, XI(I), YI(I), HI(I), CI(I)
END DO
CLOSE (3)
END IF
C----------- CREATION OF FILE FOR POST PROCESS --------
OPEN ( 2, FILE='POSTPROC.DAT', STATUS='UNKNOWN' )
DO I = 1 , NNODE
WRITE (2,*) H(I)
END DO
DO I = 1 , NIPT
WRITE(2,*) HI(I)
END DO
DO IEL = 1 , NE
WRITE(2,*) IEL, BV(IEL,1), BV(IEL,2)
END DO
CLOSE (2)
RETURN
END
C
C
SUBROUTINE SYSTEM ( MXN , N, A , C )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A (MXN,MXN) , C (MXN)
N1 = N - 1
DO K = 1, N1
L = K + 1
DO I = L , N
P = A (I,K) / A (K,K)
IF ( P .NE. 0. ) THEN
DO J = L , N
A (I,J) = A (I,J) - P * A ( K , J )
END DO
C ( I ) = C ( I) - P * C ( K )
END IF
END DO
END DO
C---- BACK SUBSTITUTION
C (N) = C (N) / A (N,N)
DO K = 1 , N1
I = N - K
L = I + 1
P = C ( I )
DO J = L , N
P = P - C (J) * A (I,J)
END DO
C ( I ) = P / A (I,I)
END DO
RETURN
END
C
C
SUBROUTINE GRULE ( N , SAI , W )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(N) , W(N), TABLESAI(8,4), TABLEW(8,4)
IF ( N .LT. 1 ) STOP'N<1'
IF ( N .GT. 8 ) STOP'N>8'
C-------- 1-POINT GAUSS-LEGENDRE
TABLESAI(1,1) = 0.D0
TABLEW(1,1) = 2.D0
C-------- 2-POINT GAUSS-LEGENDRE
TABLESAI(2,1) = DSQRT(3.D0)/3.D0
TABLEW(2,1) = 1.D0
C-------- 3-POINT GAUSS-LEGENDRE
TABLESAI(3,1) = DSQRT(15.D0)/5.D0
TABLESAI(3,2) = 0.D0
TABLEW(3,1) = 5.D0/ 9.D0
TABLEW(3,2) = 8.D0/ 9.D0
C-------- 4-POINT GAUSS-LEGENDRE
TABLESAI(4,1) = 0.33998104358485626480D0
TABLESAI(4,2) = 0.86113631159405257522D0
TABLEW(4,1) = 0.6521451548625461426D0
TABLEW(4,2) = 0.34785484513745385737D0
C-------- 5-POINT GAUSS-LEGENDRE
TABLESAI(5,1) = 0.90617984593866399279D0
TABLESAI(5,2) = 0.53846931010568309103D0
TABLESAI(5,3) = 0.D0
TABLEW(5,1) = 0.23692688505618908751D0
TABLEW(5,2) = 0.47862867049936646804D0
TABLEW(5,3) = 5.12D0 / 9.D0
C-------- 6-POINT GAUSS-LEGENDRE
TABLESAI(6,1) = 0.23861918608319690863D0
TABLESAI(6,2) = 0.66120938646626451366D0
TABLESAI(6,3) = 0.93246951420315202781D0
TABLEW(6,1) = 0.46791393457269104738D0
TABLEW(6,2) = 0.36076157304813860756D0
TABLEW(6,3) = 0.17132449237917034504D0
C-------- 7-POINT GAUSS-LEGENDRE
TABLESAI(7,1) = 0.94910791234275852452D0
TABLESAI(7,2) = 0.74153118559939443986D0
TABLESAI(7,3) = 0.40584515137739716690D0
TABLESAI(7,4) = 0.D0
TABLEW(7,1) = 0.12948496616886969327D0
TABLEW(7,2) = 0.27970539148927666790D0
TABLEW(7,3) = 0.38183005050511894495D0
TABLEW(7,4) = 0.41795918367346938775D0
C-------- 8-POINT GAUSS-LEGENDRE
TABLESAI(8,1) = 0.96028985649753623168D0
TABLESAI(8,2) = 0.79666647741362673959D0
TABLESAI(8,3) = 0.52553240991632898581D0
TABLESAI(8,4) = 0.18343464249564980493D0
TABLEW(8,1) = 0.10122853629037625915D0
TABLEW(8,2) = 0.22238103445337447054D0
TABLEW(8,3) = 0.31370664587788728733D0
TABLEW(8,4) = 0.36268378337836198296D0
C----- SAI BETWEEN 0 AND +1
NFRONT = N / 2
NREAR = N - NFRONT
DO I = 1, NREAR
SAI(I) = TABLESAI(N,I)
W(I) = TABLEW(N,I)
END DO
C----- FOR A CASE OF N=1
IF ( NFRONT .EQ. 0 ) THEN
RETURN
END IF
C----- SAI BETWEEN -1 AND 0
DO I = 1 , NFRONT
J = N - I + 1
SAI(J) = - TABLESAI(N,I)
W(J) = TABLEW(N,I)
END DO
RETURN
END
C
C
SUBROUTINE GDM (MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM,NEFEM,
* SFFEM,XCOORDFEM,YCOORDFEM,NODEXFEM,GM,THETA,C1, XP,YP, AJ )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM)
DIMENSION BPPFEM(2,NDFEM,INTEPTFEM,INTEPTFEM),WFEM(INTEPTFEM),
* SFFEM(NDFEM,INTEPTFEM,INTEPTFEM)
C---------------------------------------
CONST = 2.D0*GM*THETA
AJ = 0.D0
C---------------------------------------
DO IEL = 1 ,NEFEM
C
DO K = 1 , INTEPTFEM
DO L = 1 , INTEPTFEM
WEIGHT = WFEM(K) * WFEM(L)
YAC11 = 0.D0
YAC12 = 0.D0
YAC21 = 0.D0
YAC22 = 0.D0
DO I = 1 , NDFEM
YAC11 = YAC11 + BPPFEM(1,I,K,L) * XCOORDFEM(NODEXFEM(IEL,I))
YAC12 = YAC12 + BPPFEM(1,I,K,L) * YCOORDFEM(NODEXFEM(IEL,I))
YAC21 = YAC21 + BPPFEM(2,I,K,L) * XCOORDFEM(NODEXFEM(IEL,I))
YAC22 = YAC22 + BPPFEM(2,I,K,L) * YCOORDFEM(NODEXFEM(IEL,I))
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
BETA = WEIGHT * DETJ
C
X = 0.D0
Y = 0.D0
DO I = 1 , NDFEM
X = X + SFFEM(I,K,L)*XCOORDFEM(NODEXFEM(IEL,I))
Y = Y + SFFEM(I,K,L)*YCOORDFEM(NODEXFEM(IEL,I))
END DO
DX = X - XP
DY = Y - YP
R = DSQRT ( DX*DX + DY*DY )
IF ( R .EQ. 0.D0 ) R=1.D-12
AJ = AJ + C1*DLOG(R)*BETA*CONST
END DO
END DO
C
END DO
RETURN
END
C
C
SUBROUTINE DERIV ( ND, INTEPT, F0, F1, SAI, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT),F0(ND),F1(ND)
DX = 0.5D0
DO K = 1 , INTEPT
E1 = SAI (K)
DO L = 1 , INTEPT
E2 = SAI (L)
COMPUTATION OF BP(J) = DN(J) / DETA1
CALL ISOPARA ( ND, E1+DX , E2 , F1 )
CALL ISOPARA ( ND, E1-DX , E2 , F0 )
DO I = 1 , ND
BPP(1,I,K,L) = F1(I) - F0(I)
END DO
COMPUTATION OF BP(J)= DN(J)/DETA2
CALL ISOPARA ( ND, E1 , E2+DX , F1 )
CALL ISOPARA ( ND, E1 , E2-DX , F0 )
DO I = 1 , ND
BPP(2,I,K,L) = F1(I) - F0(I)
END DO
END DO
END DO
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND, E1 , E2 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) = 0.25D0*(1.D0 - E1)*(1.D0 - E2)
F(2) = 0.25D0*(1.D0 + E1)*(1.D0 - E2)
F(3) = 0.25D0*(1.D0 + E1)*(1.D0 + E2)
F(4) = 0.25D0*(1.D0 - E1)*(1.D0 + E2)
RETURN
END
C
C
SUBROUTINE SHAPEF ( ND, INTEPT, F, SAI, SF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND) , SAI(INTEPT) , SF(ND,INTEPT,INTEPT)
DO K = 1 , INTEPT
E1 = SAI (K)
DO L = 1 , INTEPT
E2 = SAI( L )
CALL ISOPARA ( ND, E1 , E2 , F )
DO I = 1 , ND
SF(I,K,L) = F(I)
END DO
END DO
END DO
RETURN
END
C
C
SUBROUTINE CHECK( MXN, MXI, NIPT, NE, XI, YI, X, Y, ICHECK)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION ICHECK(MXI), XI(MXI),YI(MXI), X(MXN),Y(MXN)
C======================================================================
IF ( NIPT .EQ. 0 ) RETURN
C----------------------------------------------------------------------
C- ICHECK(I)=NODE# IMPLIES THE INTERNAL POINT(I) ON THE BOUNDARY.
C- ICHECK(I)=0 IMPLIES THE INTERNAL POINT(I) NOT ON THE BOUNDARY.
DO I = 1 , NIPT
ICHECK(I) = 0
END DO
C
CALL RADIUSMX ( MXN,NE,X,Y,RMAX )
C
DO INTERNAL = 1 , NIPT
XP = XI(INTERNAL)
YP = YI(INTERNAL)
DO NODE = 1 , NE
DX = XP - X(NODE)
DY = YP - Y(NODE)
R = DSQRT ( DX*DX + DY*DY )
IF ( R/RMAX .LE. 1.D-10 ) THEN
ICHECK(INTERNAL) = NODE
EXIT
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE GRAPHIC ( NDFEM,MXEFEM,MXI,NIPT,XI,YI,HI,GM,THETA,
* NEFEM,NODEXFEM )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XI(MXI),YI(MXI), HI(MXI), NODEXFEM(MXEFEM,NDFEM)
C======================================================================
IF ( NIPT .EQ. 0 ) RETURN
C----------------------------------------------------------------------
OPEN ( 1, FILE = 'TRS4DATA.DAT', STATUS = 'UNKNOWN' )
WRITE (1,*) GM, THETA
WRITE (1,*) NEFEM
DO IEL = 1 , NEFEM
WRITE (1,*) IEL,( NODEXFEM(IEL,J), J=1,NDFEM )
END DO
WRITE (1,*) NIPT
DO I = 1 , NIPT
WRITE (1,*) I, XI(I), YI(I)
END DO
CLOSE (1)
C=======================================================================
C========> FILENAME SOLUTION4.BIN
OPEN (1,FILE="SOLUTION4.BIN",STATUS='UNKNOWN',FORM='UNFORMATTED')
WRITE (1) ( HI(I) , I = 1 , NIPT )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE MOMENT ( MXE,MXN,INTEPT,ND,BPP,W,NE,SF,
* XCOORD,YCOORD,NODEX, RHS,TORQUE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RHS(MXN),
* SS(ND,ND),BPP(2,ND,INTEPT,INTEPT),W(INTEPT),SF(ND,INTEPT,INTEPT)
TORQUE = 0.D0
C
DO IEL = 1 ,NE
C
DO K = 1 , INTEPT
DO L = 1 , INTEPT
WEIGHT = W(K) * W(L)
YAC11 = 0.D0
YAC12 = 0.D0
YAC21 = 0.D0
YAC22 = 0.D0
DO I = 1 , ND
YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODEX(IEL,I))
YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODEX(IEL,I))
YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODEX(IEL,I))
YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODEX(IEL,I))
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
BETA = WEIGHT * DETJ
DO I = 1 , ND
TORQUE = TORQUE + SF(I,K,L) * BETA*RHS(NODEX(IEL,I))
END DO
END DO
END DO
C
END DO
C
TORQUE = 2.D0*TORQUE
OPEN (1,FILE="TORQUEBEM.OUT",STATUS='UNKNOWN' )
WRITE (1,*) 'MOMENT =', TORQUE
CLOSE (1)
RETURN
END
C
C
SUBROUTINE MOMENT2 ( INTEPT,ND,MXE,MXN,C1,NE,
* GE,FE, SAI,W,NODEX,X,Y,H,BV,XE,YE,
* GM, THETA, MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM,
* NEFEM,SFFEM, XCOORDFEM,YCOORDFEM,NODEXFEM,
* INTEPTM, SAIM, WM, BPPM, SFM, TORQUE )
IMPLICIT REAL*8 ( A-H , O-Z )
C----------------------------- BEM DIMENSION --------------------------
DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
* X(MXN),Y(MXN),H(MXN),BV(MXE,ND),XE(ND),YE(ND),GE(ND),FE(ND)
C------------------- FEM DIMENSION --------------------------------
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM),
* BPPFEM(2,NDFEM,INTEPTFEM,INTEPTFEM),WFEM(INTEPTFEM),
* SFFEM(NDFEM,INTEPTFEM,INTEPTFEM)
C============== DIMENSION FOR MOMENT CALCULATION ======================
DIMENSION SAIM(INTEPTM),WM(INTEPTM)
DIMENSION BPPM(2,NDFEM,INTEPTM,INTEPTM),
* SFM(NDFEM,INTEPTM,INTEPTM)
C
TORQUE = 0.D0
C
DO IEL = 1 ,NEFEM
C
DO K = 1 , INTEPTM
DO L = 1 , INTEPTM
WEIGHT = WM(K) * WM(L)
YAC11 = 0.D0
YAC12 = 0.D0
YAC21 = 0.D0
YAC22 = 0.D0
DO I = 1 , NDFEM
YAC11 = YAC11 + BPPM(1,I,K,L) * XCOORDFEM(NODEXFEM(IEL,I))
YAC12 = YAC12 + BPPM(1,I,K,L) * YCOORDFEM(NODEXFEM(IEL,I))
YAC21 = YAC21 + BPPM(2,I,K,L) * XCOORDFEM(NODEXFEM(IEL,I))
YAC22 = YAC22 + BPPM(2,I,K,L) * YCOORDFEM(NODEXFEM(IEL,I))
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
BETA = WEIGHT * DETJ
C------------- COORDINATE OF XP WITH YP IS SOURCE PONIT.
XP = 0.D0
YP = 0.D0
DO I = 1 , NDFEM
XP = XP + SFM(I,K,L)*XCOORDFEM(NODEXFEM(IEL,I))
YP = YP + SFM(I,K,L)*YCOORDFEM(NODEXFEM(IEL,I))
END DO
C
CALL INTPOINT ( INTEPT,ND,MXE,MXN,C1,NE,
* GE,FE, SAI,W,NODEX,X,Y,H,BV,XE,YE,
* INTEPTFEM, MXEFEM,MXNFEM,NDFEM, NEFEM,XCOORDFEM, YCOORDFEM,
* GM, THETA, NODEXFEM,BPPFEM, SFFEM,WFEM ,XP,YP, SUM )
C
TORQUE = TORQUE + BETA*SUM
END DO
END DO
C
END DO
C
TORQUE = 2.D0*TORQUE
OPEN (1,FILE="TORQUEBEM2.OUT",STATUS='UNKNOWN' )
WRITE (1,*) 'MOMENT2 =', TORQUE
CLOSE (1)
RETURN
END
C
C
SUBROUTINE INTPOINT ( INTEPT,ND,MXE,MXN,C1,NE,
* GE,FE, SAI,W,NODEX,X,Y,H,BV,XE,YE,
* INTEPTFEM, MXEFEM,MXNFEM,NDFEM, NEFEM,XCOORDFEM, YCOORDFEM,
* GM, THETA, NODEXFEM,BPPFEM, SFFEM,WFEM ,XP,YP, SUM )
IMPLICIT REAL*8 ( A-H , O-Z )
C----------------------------- BEM DIMENSION --------------------------
DIMENSION NODEX(MXE,ND), SAI(INTEPT), W(INTEPT),
* X(MXN),Y(MXN),H(MXN),BV(MXE,ND),XE(ND),YE(ND),GE(ND),FE(ND)
C============================= FEM DIMENSION ==========================
DIMENSION NODEXFEM(MXEFEM,NDFEM),XCOORDFEM(MXNFEM),
* YCOORDFEM(MXNFEM), WFEM(INTEPTFEM),
* BPPFEM(2,NDFEM,INTEPTFEM,INTEPTFEM),
* SFFEM(NDFEM,INTEPTFEM,INTEPTFEM)
C======================================================================
SUM = 0.D0
C======================================================================
DO IEL = 1 , NE
DO I = 1 , ND
XE(I) = X( NODEX(IEL,I) )
YE(I) = Y( NODEX(IEL,I) )
END DO
CALL INTE ( INTEPT, ND, XP, YP, C1, XE, YE, SAI,W,GE,FE )
DO J = 1 , ND
SUM = SUM + FE(J)*H(NODEX(IEL,J)) - GE(J)*BV(IEL,J)
END DO
END DO
C======================================================================
CALL GDM ( MXEFEM,MXNFEM,INTEPTFEM,NDFEM,BPPFEM,WFEM, NEFEM,
* SFFEM,XCOORDFEM,YCOORDFEM,NODEXFEM,GM,THETA,C1, XP,YP, AJ )
SUM = SUM + AJ
C----------------------------------------------------------------------
RETURN
END
C
C
SUBROUTINE RADIUSMX ( MXN,NE,X,Y,RMAX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION X(MXN),Y(MXN)
C======================================================================
RMAX = 0.D0
C======================================================================
DO I = 1 , NE-1
XS = X(I)
YS = Y(I)
DO J = I+1, NE
DX = X(J) - XS
DY = Y(J) - YS
R = DSQRT ( DX*DX + DY*DY )
IF ( R .GT. RMAX ) RMAX = R
END DO
END DO
RETURN
END