C PCPRB7.FOR 27 February 1991 C PROGRAM PCPRB7 C C Materially nonlinear problem C C Reference: C C Ivo Babuska, Werner C Rheinboldt, C Reliable Error Estimations and Mesh adaptation for the Finite Element Method, C In: Computation Methods in Nonlinear Mechanics, C Edited by J T Oden, C North Holland Publishing Company, C Amsterdam 1980, Pages 67-108. C C The function: C C The general problem is the two point boundary value problem: C C -D/DX PHI(X,U,DU/DX,T) + PSI(X,U,DU/DX,T) = 0, C C U(X=0)=U0, U(X=1)=U1. C C The particular problem solved here has the form: C C U'' + T*SIN(U+U**2+U**3)=0 C C U(0)=U(1)=0.0 C C U is approximated by dividing the interval into NINT subintervals, C and using the NPOLYS Legendre functions of orders 0 through NPOLYS-1 C on each subinterval. The NINT*NPOLYS coefficients of these basis C functions form the representation of U. C C Moreover, there are NINT*NPOLYS finite element equations set up, C of the form C C INTEGRAL(Subinterval I) -PHI(X,U,U',T)*PL'(J) + PSI(X,U,U',T)*PL(J) = 0 C C There is a complication however, in that Lagrange multipliers are C used to enforce continuity of the solution U at the interface C nodes between neighboring subintervals. These conditions create C new variables, and new equations, and even modify the above standard C finite element equations. C C Options: C C ICHOOZ=1 Piecewise linear, 1 continuity condition C ICHOOZ=2 Piecewise cubic, 1 continuity condition C ICHOOZ=3 Piecewise cubic, 2 continuity conditions C ICHOOZ=4 Piecewise quintic, 1 continuity condition C ICHOOZ=5 Piecewise quintic, 2 continuity conditions C ICHOOZ=6 Piecewise quintic, 3 continuity conditions C C All options use 8 intervals. C INTEGER LIW INTEGER LRW INTEGER MINT INTEGER MVAR C PARAMETER (MINT=8) PARAMETER (MVAR=MINT*6+2+(MINT-1)*3) PARAMETER (LIW=MVAR+29) PARAMETER (LRW=29+(6+MVAR)*MVAR) C EXTERNAL BVAL EXTERNAL DENSLV EXTERNAL FP0011 EXTERNAL FX0011 EXTERNAL PITCON EXTERNAL UVAL C INTRINSIC MOD C REAL BCONE REAL BCZERO REAL DBCODT REAL DBCZDT REAL FPAR(1) REAL GCOEF REAL GPOINT INTEGER I INTEGER ICHOOZ INTEGER IERROR INTEGER IHOLD INTEGER II INTEGER IPAR(1) INTEGER ISKIP INTEGER IWORK(LIW) INTEGER LOUNIT CHARACTER NAME*12 INTEGER NBCZ INTEGER NBCO INTEGER NDRV INTEGER NINT INTEGER NPOLYS INTEGER NVAR INTEGER NVARY INTEGER NVARZ REAL PL REAL PLD REAL RWORK(LRW) REAL TEMP REAL THETA REAL U REAL UPRYM REAL XG REAL XL REAL XR(MVAR) REAL XRR C COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8), 1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO, 2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10) C DATA GCOEF / 1 .1012285363, .2223810345, .3137066459, .3626837834, 2 .3626837834, .3137066459, .2223810345, .1012285363/ DATA GPOINT / 1 .9602898565, .7966664774, .5255324099, .1834346425, 2 -.1834346425, -.5255324099, -.7966664774, -.9602898565/ DATA THETA/1.0, 0.0, 0.0, 0.0, 0.0,0.0, 0.0, 0.0, 0.0, 0.0, 3 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 5 1.0, 3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 7 1.0, 6.0, 15.0, 15.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 9 1.0, 10.0, 45.0, 105.0, 105.0, 0.0, 0.0, 0.0, 0.0, 0.0, 2 1.0, 15.0, 105.0, 420.0, 945.0, 945.0, 0.0, 0.0, 0.0, 0.0, 4 1., 21., 210., 1260., 4725., 10395., 10395., 0., 0., 0., 6 1., 28., 378., 3150., 17325., 62370., 135135., 135135., 0., 0., 8 1., 36., 630., 6930., 51975., * 270270., 945945., 2027025., 2027025., 0., 1 1., 45., 990., 13860., 135135., * 945945., 4729725., 16216200.,34459425.,34459425./ C LOUNIT=6 C C Zero out work arrays C DO 10 I=1,LIW IWORK(I)=0 10 CONTINUE DO 20 I=1,LRW RWORK(I)=0.0 20 CONTINUE C C Set number of variables C ICHOOZ=1 NINT=8 IF(ICHOOZ.EQ.1)NPOLYS=2 IF(ICHOOZ.EQ.2.OR.ICHOOZ.EQ.3)NPOLYS=4 IF(ICHOOZ.EQ.4.OR.ICHOOZ.EQ.5.OR.ICHOOZ.EQ.6)NPOLYS=6 NVARY=NINT*NPOLYS NBCZ=1 IF(ICHOOZ.EQ.1.OR.ICHOOZ.EQ.2.OR.ICHOOZ.EQ.4)NDRV=1 IF(ICHOOZ.EQ.3.OR.ICHOOZ.EQ.5)NDRV=2 IF(ICHOOZ.EQ.6)NDRV=3 NBCO=1 NVARZ=NBCZ+(NINT-1)*NDRV+NBCO NVAR=NVARY+NVARZ+1 C C Starting point C DO 110 I=1,NVAR XR(I)=0.0 110 CONTINUE C C We perturb the zero solution, by setting the 5-th coefficient C to some nonzero value. We will then also require that the program C find an exact solution by correcting this starting point, but C without altering the 5-th coefficient, thus guaranteeing that the C solution is not all zero. C IHOLD=5 XR(IHOLD)=0.001 C C IWORK(1)=0 ; This is a startup C IWORK(2)=IHOLD ; Use index IHOLD for first parameter C IWORK(3)=0 ; Program may choose index C IWORK(4)=0 ; Update jacobian every newton step C IWORK(5)=0 ; Seek no target values. C IWORK(6)=NVAR ; Seek limit points in index NVAR C IWORK(7)=1 ; Moderate amount of output C IWORK(8)=6 ; Output unit to use C IWORK(9)=0 ; Use user's jacobian routine C C IWORK(17)=20 ; (Unusual setting) Increase the allowed number of Newton C iterations to 20. C IWORK(1)=0 IWORK(2)=IHOLD IWORK(3)=0 IWORK(4)=0 IWORK(5)=0 IWORK(6)=NVAR IWORK(7)=1 IWORK(8)=LOUNIT IWORK(9)=0 IWORK(17)=20 C C RWORK(1)=0.0001 ; Absolute error tolerance C RWORK(2)=0.0001 ; Relative error tolerance C RWORK(3)=0.5 ; Minimum stepsize C RWORK(4)=0.5 ; Maximum stepsize C RWORK(5)=0.25 ; Starting stepsize C RWORK(6)=1.0 ; Starting direction C RWORK(7)=0.0 ; Target value C RWORK(1)=0.0001 RWORK(2)=0.0001 RWORK(3)=0.001 RWORK(4)=0.5 RWORK(5)=0.25 RWORK(6)=1.0 RWORK(7)=0.0 C WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'PCPRB7' WRITE(LOUNIT,*)'PITCON sample program' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'The materially nonlinear problem.' WRITE(LOUNIT,*)' ' IF(ICHOOZ.EQ.1)THEN WRITE(LOUNIT,*)'Piecewise linears' ELSEIF(ICHOOZ.EQ.2.OR.ICHOOZ.EQ.3)THEN WRITE(LOUNIT,*)'Piecewise cubics' ELSE WRITE(LOUNIT,*)'Piecewise quintics.' ENDIF IF(ICHOOZ.EQ.6)THEN WRITE(LOUNIT,*)'3 continuity conditions at breakpoints' ELSEIF(ICHOOZ.EQ.5.OR.ICHOOZ.EQ.3)THEN WRITE(LOUNIT,*)'2 continuity conditions at breakpoints' ELSE WRITE(LOUNIT,*)'1 continuity condition at breakpoints.' ENDIF WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Step Type of Point Lambda' WRITE(LOUNIT,*)' ' I=0 NAME='Start point ' WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,XR(NVAR) C C Take another step. C DO 30 I=1,30 C CALL PITCON(FP0011,FPAR,FX0011,IERROR,IPAR,IWORK,LIW, * NVAR,RWORK,LRW,XR,DENSLV) C IF(IWORK(1).EQ.1)THEN NAME='Corrected ' ELSEIF(IWORK(1).EQ.2)THEN NAME='Continuation ' ELSEIF(IWORK(1).EQ.3)THEN NAME='Target point ' ELSEIF(IWORK(1).EQ.4)THEN NAME='Limit point ' ENDIF C WRITE(LOUNIT,'(1X,I3,2X,A12,2X,3G14.6)')I,NAME,XR(NVAR) C C Print out shape of rod every 5 steps C IF(MOD(I,5).EQ.0)THEN WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Current rod coordinates:' WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'X, U(X), U''(X)' WRITE(LOUNIT,*)' ' DO 25 II=1,NINT TEMP=NINT XL=(II-1)/TEMP XRR=(II)/TEMP ISKIP=(II-1)*NPOLYS XG=0.5*(XL+XRR) CALL BVAL(NPOLYS,PL,PLD,XG,XL,XRR) CALL UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,XR,ISKIP) WRITE(LOUNIT,'(1X,3G14.6)')XG,U,UPRYM 25 CONTINUE WRITE(LOUNIT,*)' ' ENDIF C IF(IWORK(1).EQ.3)GO TO 40 C IF(IERROR.NE.0)THEN WRITE(LOUNIT,*)'Iteration halted, IERROR=',IERROR GO TO 40 ENDIF C 30 CONTINUE 40 CONTINUE WRITE(LOUNIT,*)' ' WRITE(LOUNIT,*)'Jacobians ',IWORK(19) WRITE(LOUNIT,*)'Factorizations ',IWORK(20) WRITE(LOUNIT,*)'Solves ',IWORK(21) WRITE(LOUNIT,*)'Functions ',IWORK(22) STOP END SUBROUTINE FX0011(NVAR,FPAR,IPAR,X,FX,IERROR) C C*********************************************************************** C EXTERNAL BD0011 EXTERNAL BVAL EXTERNAL FVAL EXTERNAL UVAL C INTEGER NVAR C REAL BCONE REAL BCZERO REAL DBCODT REAL DBCZDT REAL DTDX REAL DTDXL REAL DTDXR REAL FPAR(*) REAL FX(NVAR) REAL GCOEF REAL GPOINT REAL H2I REAL H2IL REAL H2IR INTEGER I INTEGER IEQN INTEGER IERROR INTEGER INDX INTEGER IPAR(*) INTEGER ISKIP INTEGER IVAR INTEGER J INTEGER K INTEGER KHI INTEGER L INTEGER LHIL INTEGER LHIR INTEGER NBCZ INTEGER NBCO INTEGER NCL INTEGER NCR INTEGER NDRV INTEGER NDSUM INTEGER NEQN INTEGER NINT INTEGER NODES INTEGER NPOLYS INTEGER NPSUM INTEGER NVARY INTEGER NVARZ REAL PHI REAL PHIPT REAL PHIPU REAL PHIPUP REAL PL REAL PLD REAL PSI REAL PSIPT REAL PSIPU REAL PSIPUP REAL S REAL TCON REAL TEMP REAL TERM REAL TERML REAL TERMR REAL THETA REAL U REAL UPRYM REAL X(NVAR) REAL XC REAL XG REAL XL REAL XR C COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8), 1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO, 2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10) C C Zero out FX C TCON=X(NVAR) CALL BD0011 NEQN=NVAR-1 DO 10 J=1,NEQN 10 FX(J)=0.0 C C 1. SET UP THE TERMS (A*Y) INVOLVING THE BIVARIATE FORM C C I FOR EACH INTERVAL C DO 40 I=1,NINT TEMP=NINT XL=(I-1)/TEMP XR=(I)/TEMP ISKIP=(I-1)*NPOLYS C C J FOR EACH GAUSS POINT, EVALUATE THE INTEGRAND C DO 30 J=1,8 XG=XL+(GPOINT(J)+1.0)*(XR-XL)/2.0 CALL BVAL(NPOLYS,PL,PLD,XG,XL,XR) CALL UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,X,ISKIP) CALL FVAL(XG,U,UPRYM,TCON,PHI,PHIPU,PHIPUP,PHIPT, * PSI,PSIPU,PSIPUP,PSIPT) C C L PROJECT ONTO EACH TEST FUNCTION PL(L) AND PLD(L) C IEQN=ISKIP DO 20 L=1,NPOLYS IEQN=IEQN+1 TERM=GCOEF(J)*(XR-XL)*(PSI*PL(L)+PHI*PLD(L))/2.0 FX(IEQN)=FX(IEQN)+TERM 20 CONTINUE 30 CONTINUE 40 CONTINUE C C 2. ADD THE TERMS INVOLVING THE CONTINUITY OF THE TEST FUNCTIONS C WHICH ARE THE TERMS B*Z IN F=A*Y + B*Z C NCL=NVARY NCR=NVARY+NBCZ C C I FOR EACH INTERVAL C DO 110 I=1,NINT TEMP=NINT XL=(I-1)/TEMP XR=(I)/TEMP DTDX=2.0/(XR-XL) C C J FOR THE POLYNOMIALS USED IN APPROXIMATING EACH U, C COUNT CONDITIONS AT LEFT ENDPOINT, LHIL, AND AT RIGHT, LHIR C IF WE ARE IN THE FIRST OR LAST INTERVAL, ONE OF C THESE WILL BE BOUNDARY CONDITIONS C LHIL=NDRV LHIR=NDRV IF (I.EQ.1) LHIL=NBCZ IF (I.EQ.NINT) LHIR=NBCO IF (LHIL.LE.0.AND.LHIR.LE.0) GO TO 100 C C K FOR EACH TEST FUNCTION PL(K) C DO 90 K=1,NPOLYS TERML=0.0 TERMR=0.0 S=(-1.0)**(K+1) IEQN=(I-1)*NPOLYS+K C C L FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0 TO BE ZERO C OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS C AT INTERIOR NODES C OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0 TO BE ZERO C IF (LHIL.LE.0) GO TO 60 H2I=1.0 DO 50 L=1,LHIL S=-S INDX=NCL+L TERML=TERML+S*X(INDX)*H2I*THETA(L,K) H2I=H2I*DTDX 50 CONTINUE 60 IF (LHIR.LE.0) GO TO 80 H2I=1.0 DO 70 L=1,LHIR INDX=NCR+L TERMR=TERMR+H2I*THETA(L,K)*X(INDX) H2I=H2I*DTDX 70 CONTINUE 80 FX(IEQN)=FX(IEQN)+TERML+TERMR 90 CONTINUE IF (LHIL.GT.0) NCL=NCL+LHIL IF (LHIR.GT.0) NCR=NCR+LHIR 100 CONTINUE 110 CONTINUE C C 3. CREATE THE TERMS FOR THE U FUNCTIONS AND THEIR DERIVATIVES C THE MATRIX TERMS ( C*Y ) C ONE EQUATION IS GENERATED FOR COMPONENT AND CONDITION C NPSUM=0 DTDXR=0.0 DTDXL=0.0 C C I FOR EACH NODE C NDSUM=NVARY NODES=NINT+1 DO 210 I=1,NODES TEMP=NINT IF (I.GT.1) XL=(I-2)/TEMP XC=(I-1)/TEMP IF (I.LT.NODES) XR=(I)/TEMP IF (XC.NE.XL) DTDXL=2.0/(XC-XL) IF (XR.NE.XC) DTDXR=2.0/(XR-XC) H2IL=1.0 H2IR=1.0 C C K FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0 C OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS C AT INTERIOR NODES C OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0 C KHI=NDRV IF (I.EQ.1) KHI=NBCZ IF (I.EQ.NODES) KHI=NBCO IF (KHI.LE.0) GO TO 190 DO 180 K=1,KHI S=(-1.0)**(K+1) C C L SET UP THE TERM FROM THE LEFT HAND INTERVAL C TERML=0.0 IF (I.GT.1) GO TO 120 TERML=BCZERO(K) GO TO 140 120 DO 130 L=1,NPOLYS IVAR=NPSUM+L-NPOLYS TERML=TERML+X(IVAR)*H2IL*THETA(K,L) 130 CONTINUE C C L SET UP THE TERM FROM THE RIGHT HAND INTERVAL C 140 IF (I.LT.NODES) GO TO 150 TERMR=-BCONE(K) GO TO 170 150 TERMR=0.0 DO 160 L=1,NPOLYS IVAR=NPSUM+L S=-S TERMR=TERMR+S*X(IVAR)*H2IR*THETA(K,L) 160 CONTINUE 170 IEQN=NDSUM+K FX(IEQN)=TERML+TERMR H2IL=H2IL*DTDXL H2IR=H2IR*DTDXR 180 CONTINUE NDSUM=NDSUM+KHI 190 NPSUM=NPSUM+NPOLYS 200 CONTINUE 210 CONTINUE RETURN END SUBROUTINE FP0011(NVAR,FPAR,IPAR,X,FPRIME,IERROR) C C*********************************************************************** C C JACOBIAN OF MATERIALLY NONLINEAR PROBLEM. C EXTERNAL BD0011 EXTERNAL BVAL EXTERNAL FVAL EXTERNAL UVAL C INTEGER NVAR C REAL BCONE REAL BCZERO REAL DBCODT REAL DBCZDT REAL DTDX REAL FPAR(*) REAL FPRIME(NVAR,NVAR) REAL GCOEF REAL GPOINT REAL H2I INTEGER I INTEGER IEQN INTEGER IERROR INTEGER IPAR(*) INTEGER ISKIP INTEGER IVAR INTEGER J INTEGER K INTEGER KHIL INTEGER KHIR INTEGER L INTEGER LHIL INTEGER LHIR INTEGER N INTEGER NBCZ INTEGER NBCO INTEGER NCL INTEGER NCR INTEGER NDRV INTEGER NINT INTEGER NPOLYS INTEGER NPSUM INTEGER NVARY INTEGER NVARZ REAL PHI REAL PHIPT REAL PHIPU REAL PHIPUP REAL PL REAL PLD REAL PSI REAL PSIPT REAL PSIPU REAL PSIPUP REAL S REAL SUM REAL TCON REAL TEMP REAL TERM REAL TERM1 REAL TERM2 REAL TERM3 REAL TERM4 REAL THETA REAL U REAL UPRYM REAL X(NVAR) REAL XG REAL XL REAL XR C COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8), 1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO, 2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10) C C ZERO OUT THE MATRIX C TCON=X(NVAR) CALL BD0011 DO 20 I=1,NVAR DO 10 J=1,NVAR FPRIME(I,J)=0.0 10 CONTINUE 20 CONTINUE C C 1. SET UP THE TERMS FROM THE BIVARIATE FORM (A*Y ) C C I FOR EACH INTERVAL (XL,XR) C DO 80 I=1,NINT TEMP=NINT XL=(I-1)/TEMP XR=(I)/TEMP ISKIP=(I-1)*NPOLYS C C J FOR EACH GAUSS POINT XG IN THE INTERVAL C NOTE THAT THE LEGENDRE POLYNOMIALS AND DERIVATIVES ARE SET UP C BY THE CALL TO VALUE C DO 70 J=1,8 XG=XL+0.5*(XR-XL)*(GPOINT(J)+1.0) CALL BVAL(NPOLYS,PL,PLD,XG,XL,XR) CALL UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,X,ISKIP) CALL FVAL(XG,U,UPRYM,TCON,PHI,PHIPU,PHIPUP,PHIPT, * PSI,PSIPU,PSIPUP,PSIPT) C C L FOR EACH LEGENDRE POLYNOMIAL COEFFICIENT C IEQN=ISKIP DO 50 L=1,NPOLYS IEQN=IEQN+1 TERM=GCOEF(J)*0.5*(XR-XL)*(PSIPT*PL(L)+PHIPT*PLD(L)) FPRIME(IEQN,NVAR)=FPRIME(IEQN,NVAR)+TERM C C N FOR EACH Y-COEFFICIENT OF A U C DO 30 N=1,NPOLYS IVAR=NPOLYS*(I-1)+N TERM1=PSIPU*PL(N)*PL(L) TERM2=PSIPUP*PLD(N)*PL(L) TERM3=PHIPU*PL(N)*PLD(L) TERM4=PHIPUP*PLD(N)*PLD(L) SUM=TERM1+TERM2+TERM3+TERM4 FPRIME(IEQN,IVAR)=FPRIME(IEQN,IVAR)+ * GCOEF(J)*0.5*(XR-XL)*SUM 30 CONTINUE 50 CONTINUE 70 CONTINUE 80 CONTINUE C C 2. ADD THE TERMS INVOLVING THE CONTINUITY OF THE TEST FUNCTIONS C WHICH ARE THE TERMS B*Z IN F=A*Y+B*Z C C C I FOR EACH INTERVAL, C DO 150 I=1,NINT NCL=NVARY IF (I.GT.1) NCL=NVARY+NBCZ+(I-2)*NDRV NCR=NVARY+NBCZ+(I-1)*NDRV TEMP=NINT XL=(I-1)/TEMP XR=(I)/TEMP DTDX=2.0/(XR-XL) C C J FOR THE POLYNOMIALS USED IN APPROXIMATING EACH U, C COUNT CONDITIONS AT LEFT ENDPOINT, LHIL, AND AT RIGHT, LHIR C IF WE ARE IN THE FIRST OR LAST INTERVAL, ONE OF C THESE WILL BE BOUNDARY CONDITIONS C LHIL=NDRV LHIR=NDRV IF (I.EQ.1) LHIL=NBCZ IF (I.EQ.NINT) LHIR=NBCO IF (LHIL.LE.0.AND.LHIR.LE.0) GO TO 130 C C K FOR EACH TEST FUNCTION PL(K) C DO 120 K=1,NPOLYS S=(-1.0)**(K+1) IEQN=(I-1)*NPOLYS+K C C L FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0 TO BE ZERO C OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS C AT INTERIOR NODES C OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0 TO BE ZERO C C EVALUATE CONTRIBUTION FROM LEFT ENDPOINT C IF (LHIL.LE.0) GO TO 100 H2I=1.0 DO 90 L=1,LHIL S=-S IVAR=NCL+L FPRIME(IEQN,IVAR)=S*H2I*THETA(L,K) H2I=H2I*DTDX 90 CONTINUE C C EVALUATE CONTRIBUTION FROM RIGHT ENDPOINT C 100 IF (LHIR.LE.0) GO TO 120 H2I=1.0 DO 110 L=1,LHIR IVAR=NCR+L FPRIME(IEQN,IVAR)=H2I*THETA(L,K) H2I=H2I*DTDX 110 CONTINUE 120 CONTINUE IF (LHIL.GT.0) NCL=NCL+LHIL IF (LHIR.GT.0) NCR=NCR+LHIR 130 CONTINUE 140 CONTINUE 150 CONTINUE C C 3. CREATE THE TERMS FOR THE U FUNCTIONS AND THEIR DERIVATIVES C THE MATRIX TERMS ( C*Y ) C ONE EQUATION IS GENERATED FOR COMPONENT AND CONDITION C C C I FOR EACH INTERVAL C DO 230 I=1,NINT NCL=NVARY IF (I.GT.1) NCL=NVARY+NBCZ+(I-2)*NDRV NCR=NVARY+NBCZ+(I-1)*NDRV NPSUM=(I-1)*NPOLYS TEMP=NINT XL=(I-1)/TEMP XR=(I)/TEMP DTDX=2.0/(XR-XL) H2I=1.0 C C K FOR DERIVATIVES 0 THRU NBCZ SPECIFIED AT 0.0 C OR DERIVATIVES 0 THRU NDRV SPECIFIED TO BE CONTINUOUS C AT INTERIOR NODES C OR DERIVATIVES 0 THRU NBCO SPECIFIED AT 1.0 C KHIL=NDRV IF (I.EQ.1) KHIL=NBCZ IF (KHIL.LE.0) GO TO 180 DO 170 K=1,KHIL IEQN=NCL+K C C L SET UP THE TERM FROM THE LEFT HAND ENDPOINT C IF (I.EQ.1) FPRIME(IEQN,NVAR)=DBCZDT(K) S=(-1.0)**(K+1) DO 160 L=1,NPOLYS IVAR=NPSUM+L S=-S FPRIME(IEQN,IVAR)=S*H2I*THETA(K,L) 160 CONTINUE H2I=H2I*DTDX 170 CONTINUE NCL=NCL+KHIL 180 CONTINUE H2I=1.0 KHIR=NDRV IF (I.EQ.NINT) KHIR=NBCO IF (KHIR.LE.0) GO TO 210 DO 200 K=1,KHIR IEQN=NCR+K FPRIME(IEQN,NVAR)=0.0 IF (I.EQ.NINT) FPRIME(IEQN,NVAR)=-DBCODT(K) DO 190 L=1,NPOLYS IVAR=NPSUM+L FPRIME(IEQN,IVAR)=H2I*THETA(K,L) 190 CONTINUE H2I=H2I*DTDX 200 CONTINUE NCR=NCR+KHIR 210 NPSUM=NPSUM+NPOLYS 230 CONTINUE RETURN END SUBROUTINE BVAL(NPOLYS,PL,PLD,XG,XL,XR) C C*********************************************************************** C C AUXILLIARY ROUTINE FOR MATERIALLY NONLINEAR PROBLEM. C C EVALUATE LEGENDRE POLYNOMIALS AND DERIVATIVES C REAL A INTEGER I INTEGER NPOLYS REAL PL(8) REAL PLD(8) REAL T REAL XG REAL XL REAL XR C T=-1.0+2.0*(XG-XL)/(XR-XL) PL(1)=1.0 PLD(1)=0.0 PL(2)=T PLD(2)=1.0 A=0.0 DO 10 I=3,8 A=A+1.0 PL(I)=((A+A+1.0)*T*PL(I-1)-A*PL(I-2))/(A+1.0) PLD(I)=((A+A+1.0)*(T*PLD(I-1)+PL(I-1))-A*PLD(I-2))/(A+1.0) 10 CONTINUE C DO 30 I=1,8 PLD(I)=2.0*PLD(I)/(XR-XL) 30 CONTINUE RETURN END SUBROUTINE UVAL(NPOLYS,PL,PLD,XG,U,UPRYM,X,ISKIP) C C SUM ON COEFFICENTS TO EVALUATE U AND UPRYM C INTEGER NPOLYS C INTEGER INDX INTEGER ISKIP INTEGER J REAL PL(NPOLYS) REAL PLD(NPOLYS) REAL U REAL UPRYM REAL X(*) REAL XG C U=0.0 UPRYM=0.0 DO 40 J=1,NPOLYS INDX=ISKIP+J U=U+X(INDX)*PL(J) UPRYM=UPRYM+X(INDX)*PLD(J) 40 CONTINUE RETURN END SUBROUTINE FVAL(XG,U,UPRYM,TCON,PHI,PHIPU,PHIPUP,PHIPT, * PSI,PSIPU,PSIPUP,PSIPT) C C*********************************************************************** C C AUXILLIARY ROUTINE FOR MATERIALLY NONLINEAR PROBLEM. C INTRINSIC COS INTRINSIC SIN C REAL COS REAL PHI REAL PHIPU REAL PHIPUP REAL PHIPT REAL PSI REAL PSIPU REAL PSIPUP REAL PSIPT REAL SIN REAL TCON REAL U REAL UPRYM REAL XG C PHI=-UPRYM PHIPU=0.0 PHIPUP=-1.0 PHIPT=0.0 C PSI=TCON*SIN(U*(1.0+U*(1.0+U))) PSIPU=TCON*(1.0+U*(2.0+3.0*U))*COS(U*(1.0+U*(1.0+U))) PSIPUP=0.0 PSIPT=SIN(U*(1.0+U*(1.0+U))) C RETURN END SUBROUTINE BD0011 C C*********************************************************************** C C AUXILLIARY ROUTINE SETS BOUNDARY VALUES C FOR MATERIALLY NONLINEAR PROBLEM. C REAL BCONE REAL BCZERO REAL DBCODT REAL DBCZDT REAL GCOEF REAL GPOINT INTEGER NBCZ INTEGER NBCO INTEGER NDRV INTEGER NINT INTEGER NPOLYS INTEGER NVARY INTEGER NVARZ REAL PL REAL PLD REAL THETA C COMMON /AUXMEM/ BCZERO(8),DBCZDT(8),BCONE(8),DBCODT(8), 1 NPOLYS,NDRV,NVARY,NVARZ,NBCZ,NBCO, 2 PL(8),PLD(8),GCOEF(8),GPOINT(8),NINT,THETA(10,10) C BCONE(1)=0.0 BCZERO(1)=0.0 DBCODT(1)=0.0 DBCZDT(1)=0.0 RETURN END