PROGRAM TORSION9Q
C=======================================================================
C 2-DIM FEM PROGRAM FOR SOLVING TORSION EQUATION
C USING UPPER HALF BANDED MATRIX
C EQUATION: DFAI/DX + DFAI/DY + 2*THETA*G = 0.
C ELEMENT : 9-NODED ISO-PARAMETERIC
C SEQUENCE: (-1,-1),(0,-1),(+1,-1),(+1,0),(+1,+1),(0,+1),
C (-1,+1),(-1,0),(0,0)
C ORIGINAL:1984 EIJI FUKUMORI BUFFALO NY & REVISED: 1994 ACHI
C=======================================================================
IMPLICIT REAL*8 ( A-H , O-Z )
PARAMETER ( ND=9,MXE=20000,MXN=20000,MXB=3000,INTEPT=3,MXW=500 )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), IBND(MXB),
* ITYPE(MXB), BVALUE(MXB),RK(MXN,MXW), RHS(MXN),SK(ND,ND),
* X(ND), Y(ND),BX(2,ND),SAI(INTEPT), W(INTEPT),SS(ND,ND),
* BPP(2,ND,INTEPT,INTEPT), SF(ND,INTEPT,INTEPT),
* DUMMY0(ND), DUMMY1(ND), XS(ND), YS(ND), BPPX(2,ND,ND)
DIMENSION SOURCE(MXN),TAUZX(MXN),TAUZY(MXN),NDREPEAT(MXN)
C=======================================================================
WRITE (*,*)' TORSIN9Q.FOR SOLVER'
C=======================================================================
CALL GRULE ( INTEPT, SAI, W )
CALL DERIV ( ND, INTEPT, X, Y, SAI, BPP )
CALL SHAPEF( ND, INTEPT, X, SAI, SF )
C=======================================================================
WRITE(*,*)' INPUT FILE NAME =TRS9DATA.DAT'
CALL INPUT ( ND,MXE,MXN,MXB,NE,NNODE,NB,GM,THETA, NODEX,
* XCOORD,YCOORD,IBND,ITYPE,BVALUE )
C=======================================================================
CALL BANDWID ( MXE, ND, NE,NODEX, NBW )
IF ( NBW .GT. MXW ) STOP 'Error ... NBW > MXW'
C=======================================================================
CALL DAREA ( MXE,MXN,INTEPT,ND,BPP,W,NE,XCOORD,YCOORD,
* NODEX, AREA )
C=======================================================================
CALL GSM ( INTEPT,ND,NBW ,MXE,MXN,MXW,NE,NNODE,BPP,W,BX,SS,
* NODEX,XCOORD,YCOORD,RK )
C=======================================================================
CALL SRCTRM ( MXE,MXN,INTEPT,ND,NNODE,BPP,W,NE,GM,THETA,
* SF,XCOORD,YCOORD,NODEX,SS, SOURCE )
C=======================================================================
CALL FORM (MXN,MXB,MXW,NNODE,NB,NBW,RK,RHS,ITYPE,BVALUE,IBND,
* SOURCE)
C=======================================================================
CALL SYSTEM ( MXN, MXW, NNODE, NBW, RK, RHS )
C=======================================================================
CALL MOMENT ( MXE,MXN,INTEPT,ND,BPP,W,NE,SF,
* XCOORD,YCOORD,NODEX, RHS,TORQUE )
C=======================================================================
CALL STRESS ( ND,MXE,MXN,NE,NNODE,BPPX,NODEX,XCOORD,YCOORD,
* XS, YS, DUMMY0, DUMMY1, RHS, NDREPEAT,TAUZX,TAUZY )
C=======================================================================
CALL PRINT ( MXN, NNODE , RHS, TAUZX,TAUZY )
CALL PRINTG ( ND,MXE,MXN, NNODE,NE, XCOORD, YCOORD,
* GM,THETA, NODEX,MXB, NB, IBND,ITYPE,BVALUE,
* RHS,TORQUE,TAUZX,TAUZY, AREA )
C=======================================================================
STOP 'NORMAL TERMINATION'
END
C
C
SUBROUTINE ISOPARA ( ND, E1 , E2 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) = 0.25D0*E1*(1.D0- E1)*E2*(1.D0- E2)
F(2) = -0.50D0*(1.D0- E1*E1)*E2*(1.D0- E2)
F(3) = -0.25D0*E1*(1.D0+ E1)*E2*(1.D0- E2)
F(4) = 0.50D0*E1*(1.D0+ E1)*(1.D0- E2*E2)
F(5) = 0.25D0*E1*(1.D0+ E1)*E2*(1.D0+ E2)
F(6) = 0.50D0*(1.D0- E1*E1)*E2*(1.D0+ E2)
F(7) = -0.25D0*E1*(1.D0- E1)*E2*(1.D0+ E2)
F(8) = -0.50D0*E1*(1.D0- E1)*(1.D0- E2*E2)
F(9) = (1.D0- E1*E1)*(1.D0- E2*E2)
RETURN
END
C
C
SUBROUTINE GSM ( INTEPT,ND,NBW ,MXE,MXN,MXW,NE,NNODE,BPP,W,BX,SS,
* NODEX,XCOORD,YCOORD,RK )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RK(MXN,MXW),
* SS(ND,ND), BX(2,ND),W(INTEPT), BPP(2,ND,INTEPT,INTEPT)
C====== RESET GLOBAL MATRIX RK ======
DO I = 1, NNODE
DO J = 1, NBW
RK(I,J) = 0.D0
END DO
END DO
C===== INTEGRATION =======
DO 500 IEL = 1 ,NE
DO I = 1 , ND
DO J = I , ND
SS(I,J) = 0.D0
END DO
END DO
DO 400 K = 1 , INTEPT
DO 300 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
Y11 = YAC22 / DETJ
Y12 = - YAC12 / DETJ
Y21 = - YAC21 / DETJ
Y22 = YAC11 / DETJ
DO J = 1 , ND
BX(1,J) = Y11 * BPP(1,J,K,L) + Y12 * BPP(2,J,K,L)
BX(2,J) = Y21 * BPP(1,J,K,L) + Y22 * BPP(2,J,K,L)
END DO
BETA = WEIGHT * DETJ
DO I = 1 , ND
DO J = I , ND
SS(I,J) = SS(I,J) + ( BX(1,I)*BX(1,J)+BX(2,I)*BX(2,J) )*BETA
END DO
END DO
300 CONTINUE
400 CONTINUE
DO I = 1 , ND
DO J = I , ND
SS(J,I) = SS(I,J)
END DO
END DO
C----ASSEMBLY
DO K = 1 , ND
I = NODEX(IEL,K)
DO L = 1 , ND
J = NODEX(IEL,L) - I + 1
IF ( J .GT. 0 ) RK(I,J) = RK(I,J) - SS(K,L)
END DO
END DO
500 CONTINUE
RETURN
END
C
C
SUBROUTINE INPUT ( ND,MXE,MXN,MXB,NE,NNODE,NB,GM,THETA,NODEX,
* XCOORD,YCOORD,IBND,ITYPE,BVALUE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBND(MXB),
* ITYPE(MXB), BVALUE(MXB)
IR = 1
OPEN ( IR, FILE='TRS9DATA.DAT', STATUS='UNKNOWN' )
READ (IR,*) GM, THETA
READ (IR,*) NE
DO I = 1 , NE
READ (IR,*) IEL,( NODEX(IEL,J),J=1,ND )
END DO
READ (IR,*) NNODE
READ (IR,*) ( NODE,XCOORD(NODE),YCOORD(NODE),J=1,NNODE )
READ (IR,*) NB
READ (IR,*) (IBND(I), ITYPE(I), BVALUE(I),I=1,NB)
CLOSE (IR)
RETURN
END
C
C
SUBROUTINE PRINT ( MXN, NNODE , RHS, TAUZX,TAUZY )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION RHS(MXN), TAUZX(MXN),TAUZY(MXN)
IW = 1
OPEN (IW,FILE='SOLUTION9.BIN',STATUS='UNKNOWN',FORM='UNFORMATTED')
WRITE (IW) (RHS(I),I=1,NNODE)
WRITE (IW) (TAUZX(I),I=1,NNODE)
WRITE (IW) (TAUZY(I),I=1,NNODE)
CLOSE (IW)
RETURN
END
C
C
SUBROUTINE PRINTG ( ND,MXE,MXN, NNODE,NE, XCOORD, YCOORD,
* GM,THETA, NODEX,MXB, NB, IBND,ITYPE,BVALUE,
* RHS,TORQUE,TAUZX,TAUZY, AREA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION RHS(MXN), XCOORD(MXN), YCOORD(MXN),NODEX(MXE,ND),
* IBND(MXB),ITYPE(MXB), BVALUE(MXB), TAUZX(MXN),TAUZY(MXN)
IW = 1
OPEN ( IW, FILE='SOLUTION9.FEM', STATUS='UNKNOWN' )
WRITE (IW,*) 'GM,THETA=', GM,THETA
WRITE (IW,*)
WRITE (IW,*) 'NOMBER-OF-ELEMENTS:', NE
DO I = 1, NE
WRITE (IW,*) I, (NODEX(I,J),J=1,ND)
END DO
WRITE (IW,*)
WRITE (IW,*) 'NOMBER-OF-NODES:',NNODE
DO I = 1, NNODE
WRITE (IW,*) I, XCOORD(I), YCOORD(I)
END DO
WRITE (IW,*)
WRITE (IW,*) 'RESULTS:AREA-OF-DOMAIN =', AREA
WRITE (IW,*)
WRITE (IW,*) 'NOMBER-OF-BOUNDARY-NODES:',NB
DO I = 1, NB
WRITE (IW,*) IBND(I),ITYPE(I), BVALUE(I)
END DO
WRITE (IW,*)
WRITE (IW,*) 'RESULTS:STRESS-FUNCTION'
WRITE (IW,*) 'NODE X Y STRESS-FUNCTION'
DO I = 1, NNODE
WRITE (IW,*) I, XCOORD(I), YCOORD(I), RHS(I)
END DO
WRITE (IW,*)
WRITE (IW,*) 'RESULTS:TORQUE=', TORQUE
WRITE (IW,*)
WRITE (IW,*) 'RESULTS: NODE X Y TAUZX TAUZY'
DO I = 1, NNODE
WRITE (IW,*) I, XCOORD(I), YCOORD(I),TAUZX(I), TAUZY(I)
END DO
CLOSE (IW)
RETURN
END
C
C
SUBROUTINE BANDWID ( MXE, ND, NE, NODEX , NBW )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND)
NBW = 0
DO 30 I = 1 , NE
DO 20 J = 1 , ND-1
DO 10 K = J+1 , ND
NBW = MAX0(NBW,IABS(NODEX(I,J)-NODEX(I,K)))
10 CONTINUE
20 CONTINUE
30 CONTINUE
NBW = NBW + 1
write(6,*) ' bandwidth =', nbw
RETURN
END
C
C
SUBROUTINE SYSTEM ( MXN, MXW, NUMNP, MBAND, A, B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXW) , B(MXN)
C---------- ELIMINATION ------------------
DO N = 1 , NUMNP
DO L = 2 , MBAND
C = A(N,L) / A(N,1)
I = N + L - 1
IF ( I .LE. NUMNP ) THEN
J = 0
DO K = L , MBAND
J = J + 1
A(I,J) = A(I,J) - C * A(N,K)
END DO
A(N,L) = C
B(I) = B(I) - A(N,L) * B(N)
ENDIF
END DO
B(N) = B(N) / A(N,1)
END DO
C---------- BACKSUBSTITUTION -------------
N = NUMNP
DO WHILE ( N .GT. 0 )
DO K = 2 , MBAND
L = N + K - 1
IF ( L .LE. NUMNP ) B(N) = B(N)-A(N,K)*B(L)
END DO
N = N - 1
END DO
RETURN
END
C
C
SUBROUTINE FORM (MXN,MXB,MXW,NNODE,NB,NBW,A,RHS,IBC,BV,IBND,
* SOURCE)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION RHS(MXN),A(MXN,MXW),IBC(MXB),BV(MXB),IBND(MXB),
* SOURCE(MXN)
C BOUNDARY CONDITIONS ARE COUPLED HERE.
C ITYPE(I) = 1 ---> DIRICHLET, ITYPE(I) = 2 ---> NEUMANN
C----- REFORMATION OF VECTOR RHS DUE TO FIRST KIND B.C.
C------ CLEAR RIGHT HAND SIDE.
DO I = 1 , NNODE
RHS (I) = -SOURCE(I)
END DO
DO 50 K = 1 , NB
IF ( IBC(K) .EQ. 2 ) GO TO 50
I = IBND(K)
DO 20 J = 2 , NBW
I = I - 1
IF ( I.LE. 0 ) GO TO 30
RHS(I) = RHS(I) - BV(K) * A(I,J)
20 CONTINUE
30 I = IBND(K)
DO 40 J = 2 , NBW
I = I + 1
IF ( I .GT. NNODE ) GO TO 50
RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
40 CONTINUE
50 CONTINUE
C-----REFORMATION OF MATRIX A.
DO 70 K = 1 , NB
I = IBND (K)
IF ( IBC(K) .EQ. 2 ) GO TO 65
RHS(I) = BV(K)
A(I,1) = 1.
DO 60 J = 2 , NBW
L = I - J + 1
A(I,J) = 0.
IF ( L .LE. 0 ) GO TO 60
A( L ,J) = 0.
60 CONTINUE
GO TO 70
65 RHS (I) = RHS(I) - BV(K)
70 CONTINUE
RETURN
END
C
C
SUBROUTINE GRULE ( N , SAI , W )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(N) , W(N)
IF ( N .LT. 2 ) STOP'N<2'
IF ( N .GT. 6 ) STOP'N>6'
IF ( N .EQ. 2 ) THEN
SAI(1) = DSQRT(3.D0)/3.D0
W(1) = 1.D0
SAI(2) = - SAI(1)
W(2) = W(1)
RETURN
END IF
IF ( N .EQ. 3 ) THEN
SAI(1) = DSQRT(15.D0)/5.D0
SAI(2) = 0.D0
W(1) = 5.D0/ 9.D0
W(2) = 8.D0/ 9.D0
SAI(3) = - SAI(1)
W(3) = W(1)
RETURN
END IF
IF ( N .EQ. 4 ) THEN
SAI(1) = 0.33998104358485D0
SAI(2) = 0.86113631159405D0
W(1) = 0.65214515486254D0
W(2) = 0.34785484513745D0
SAI(3) = -SAI(2)
SAI(4) = -SAI(1)
W(3) = W(2)
W(4) = W(1)
RETURN
END IF
IF ( N .EQ. 5 ) THEN
SAI(1) = 0.90617984593866D0
SAI(2) = 0.53846931010568D0
SAI(3) = 0.D0
W(1) = 0.23692688505619D0
W(2) = 0.47862867049937D0
W(3) = 5.12D0 / 9.D0
SAI(4) = -SAI(2)
SAI(5) = -SAI(1)
W(4) = W(2)
W(5) = W(1)
RETURN
END IF
IF ( N .EQ. 6 ) THEN
SAI(1) = 0.23861918608320D0
SAI(2) = 0.66120938646626D0
SAI(3) = 0.93246951420315D0
W(1) = 0.46791393457269D0
W(2) = 0.36076157304814D0
W(3) = 0.17132449237917D0
SAI(4) = -SAI(1)
SAI(5) = -SAI(2)
SAI(6) = -SAI(3)
W(4) = W(1)
W(5) = W(2)
W(6) = W(3)
END IF
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)
DO 40 K = 1 , INTEPT
E1 = SAI (K)
DO 30 L = 1 , INTEPT
E2 = SAI (L)
COMPUTATION OF BP(J) = DN(J) / DETA1
CALL ISOPARA ( ND, E1+0.5D0 , E2 , F1 )
CALL ISOPARA ( ND, E1-0.5D0 , 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+0.5D0 , F1 )
CALL ISOPARA ( ND, E1 , E2-0.5D0 , F0 )
DO I = 1 , ND
BPP(2,I,K,L) = F1(I) - F0(I)
END DO
30 CONTINUE
40 CONTINUE
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 40 K = 1 , INTEPT
E1 = SAI (K)
DO 30 L = 1 , INTEPT
E2 = SAI( L )
CALL ISOPARA ( ND, E1 , E2 , F )
DO 20 I = 1 , ND
SF(I,K,L) = F(I)
20 CONTINUE
30 CONTINUE
40 CONTINUE
RETURN
END
C
C
SUBROUTINE SRCTRM ( MXE,MXN,INTEPT,ND,NNODE,BPP,W,NE,GM,THETA,
* SF,XCOORD,YCOORD,NODEX,SS, SOURCE )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* SOURCE(MXN),
* SS(ND,ND),BPP(2,ND,INTEPT,INTEPT),W(INTEPT),SF(ND,INTEPT,INTEPT)
C
TRSCONST = 2.*GM*THETA
DO I = 1, NNODE
SOURCE(I) = 0.D0
END DO
DO 500 IEL = 1 ,NE
DO I = 1 , ND
DO J = 1 , ND
SS(I,J) = 0.D0
END DO
END DO
DO 400 K = 1 , INTEPT
DO 300 L = 1 , INTEPT
WEIGHT = W(K) * W(L) * TRSCONST
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
DO J = 1 , ND
SS(I,J) = SS(I,J) + SF(I,K,L) * SF(J,K,L) * BETA
END DO
END DO
300 CONTINUE
400 CONTINUE
C
DO K = 1 , ND
SUM = 0.D0
DO L = 1 , ND
SUM = SUM + SS(K,L)
END DO
SOURCE(NODEX(IEL,K)) = SOURCE(NODEX(IEL,K)) + SUM
END DO
500 CONTINUE
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
DO 500 IEL = 1 ,NE
DO 400 K = 1 , INTEPT
DO 300 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
300 CONTINUE
400 CONTINUE
500 CONTINUE
TORQUE = 2.D0*TORQUE
RETURN
END
C
C
SUBROUTINE ELNDPT ( ND, XS, YS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XS(ND), YS(ND)
IF ( ND .LT. 4 ) STOP'ND<4'
IF ( ND .GT. 9 ) STOP'ND>9'
IF ( ND .EQ. 4) THEN
XS(1) = -1.
XS(2) = 1.
XS(3) = 1.
XS(4) = -1.
YS(1) = -1.
YS(2) = -1.
YS(3) = 1.
YS(4) = 1.
END IF
IF ( ND .EQ. 9) THEN
XS(1) = -1.
XS(2) = 0.
XS(3) = 1.
XS(4) = 1.
XS(5) = 1.
XS(6) = 0.
XS(7) = -1.
XS(8) = -1.
XS(9) = 0.
YS(1) = -1.
YS(2) = -1.
YS(3) = -1.
YS(4) = 0.
YS(5) = 1.
YS(6) = 1.
YS(7) = 1.
YS(8) = 0.
YS(9) = 0.
END IF
RETURN
END
C
C
SUBROUTINE STRESS ( ND ,MXE,MXN,NE,NNODE, BPPX,NODEX,
* XCOORD,YCOORD,
* XS, YS, DUMMY0, DUMMY1, RHS, NDREPEAT,TAUZX,TAUZY )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),RHS(MXN),
* BPPX(2,ND,ND),XS(ND),YS(ND),DUMMY0(ND),DUMMY1(ND),
* TAUZX(MXN),TAUZY(MXN), NDREPEAT(MXN)
C====================================
CALL ELNDPT ( ND, XS, YS )
CALL DERIVX ( ND, DUMMY0, DUMMY1, XS, YS, BPPX )
C===== INTEGRATION =======
DO I = 1, NNODE
NDREPEAT(I) = 0
TAUZX(I) = 0.D0
TAUZY(I) = 0.D0
END DO
DO IEL = 1 ,NE
DO K = 1 , ND
YAC11 = 0.D0
YAC12 = 0.D0
YAC21 = 0.D0
YAC22 = 0.D0
DO I = 1 , ND
YAC11 = YAC11 + BPPX(1,I,K) * XCOORD(NODEX(IEL,I))
YAC12 = YAC12 + BPPX(1,I,K) * YCOORD(NODEX(IEL,I))
YAC21 = YAC21 + BPPX(2,I,K) * XCOORD(NODEX(IEL,I))
YAC22 = YAC22 + BPPX(2,I,K) * YCOORD(NODEX(IEL,I))
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
Y11 = YAC22 / DETJ
Y12 = - YAC12 / DETJ
Y21 = - YAC21 / DETJ
Y22 = YAC11 / DETJ
DTDX = 0.D0
DTDY = 0.D0
DO J = 1 , ND
DNDX = Y11*BPPX(1,J,K)+Y12*BPPX(2,J,K)
DNDY = Y21*BPPX(1,J,K)+Y22*BPPX(2,J,K)
DTDX = DTDX + DNDX*RHS(NODEX(IEL,J))
DTDY = DTDY + DNDY*RHS(NODEX(IEL,J))
END DO
NODE = NODEX(IEL,K)
TAUZX(NODE) = TAUZX(NODE) + DTDY
TAUZY(NODE) = TAUZY(NODE) - DTDX
NDREPEAT(NODE) = NDREPEAT(NODE) + 1
END DO
END DO
DO I = 1 , NNODE
TAUZX(I) = TAUZX(I)/NDREPEAT(I)
TAUZY(I) = TAUZY(I)/NDREPEAT(I)
END DO
RETURN
END
C
C
SUBROUTINE DERIVX ( ND, DUMMY0, DUMMY1, XS, YS, BPPX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XS(ND),YS(ND),BPPX(2,ND,ND),DUMMY0(ND),DUMMY1(ND)
DO 40 K = 1 , ND
E1 = XS (K)
E2 = YS (K)
COMPUTATION OF BP(J) = DN(J) / DETA1
CALL ISOPARA ( ND, E1+0.5D0 , E2 , DUMMY1 )
CALL ISOPARA ( ND, E1-0.5D0 , E2 , DUMMY0 )
DO I = 1 , ND
BPPX(1,I,K) = DUMMY1(I) - DUMMY0(I)
END DO
COMPUTATION OF BP(J)= DN(J) / DETA2
CALL ISOPARA ( ND, E1 , E2+0.5D0 , DUMMY1 )
CALL ISOPARA ( ND, E1 , E2-0.5D0 , DUMMY0 )
DO I = 1 , ND
BPPX(2,I,K) = DUMMY1(I) - DUMMY0(I)
END DO
40 CONTINUE
RETURN
END
C
C
SUBROUTINE DAREA ( MXE,MXN,INTEPT,ND,BPP,W,NE,XCOORD,YCOORD,
* NODEX, AREA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* BPP(2,ND,INTEPT,INTEPT),W(INTEPT)
AREA = 0.D0
DO IEL = 1 ,NE
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
AREA = AREA + WEIGHT * DETJ
END DO
END DO
END DO
RETURN
END