Attachment 104674 Details for Bug 73707 – test case

test case

gb8414.f (text/plain), 12.07 KB, created by Richard Henderson on 2004-10-03 04:07:11 UTC

(hide)

Description:

Filename:

MIME Type:

Creator: Richard Henderson

Created: 2004-10-03 04:07:11 UTC

Size: 12.07 KB

patch

obsolete

>C LINEARIZE SIDECONDITION
>      SUBROUTINE GB8414(SCD,PI,DT,CT,NT,MT,C0,
>     *                  P,T1,T2,K,M,N)
>      INTEGER PI,NT,MT,DT(*),M,N
>      REAL SCD(*),CT(MT),C0,P(M,N),T1(M,N),T2(M,N),K(N)
>C
>C A linear sidecondition is converted into a form
>C that contains only degrees of freedom and their coefficients
>C and a constant term.
>C
>C SCD,PI: sidecondition  (INPUT)
>C         SCD: definition of sideconditions
>C         PI : pointer to scd
>C DT,CT,NT,MT: sidecondition expressed by the dofs.  (OUTPUT)
>C              DT: dof. numbers
>C              CT: coefficients of the dofs.
>C              NT: number of dofs. involved
>C              MT: max. of NT  (INPUT)
>C P,T1,T2,K,M,N: points, directions and parameters (INPUT)
>C
>C HSe
>C ****************************************************************
>C
>      INTEGER NI,G2,J,IP,IT,IC,TYP,GT(4),L,AX1(3),AX2(3),
>     *        I1,I2,IS,IQ,IP1,IER
>      REAL Q,QT(10),AUX,RT(4),LL,UL,DK,D1,D2,D3,AUX1,AUX2,
>     *     DK1,DK2,DK3,DU(2),DV(2),A,DA,QC,Q1,Q2,DA1
>      LOGICAL LGI
>C
>      DATA LL /0.333333/
>      DATA UL /0.666666/
>      DATA AX1 /2,1,1/
>      DATA AX2 /3,3,2/
>      DATA QC /1.4/
>C
>      NT=0
>      C0=0
>      TYP=SCD(PI)
>      IF (TYP.EQ.1.OR.TYP.EQ.31.OR.TYP.EQ.61) THEN
>C        linear sidecondition
>C        --------------------
>         NI=SCD(PI+1)
>         C0=SCD(PI+NI+NI+3)
>         DO 110 J=1,NI
>            G2=SCD(PI+1+J)
>            IF (G2.GT.0) THEN
>C              g2: dof. number
>               NT=NT+1
>               IF (NT.GT.MT) GOTO 9000
>               DT(NT)=G2
>               CT(NT)=SCD(PI+2+NI+J)
>            ELSE
>C              pointer to tables p,t1,t2 etc.
>               G2=-G2
>               IP=G2/100
>               IT=(G2-100*IP)/10
>               IC=MOD(G2,10)
>               IF (IT.EQ.1) THEN
>                  Q=P(IC,IP)
>               ELSEIF (IT.EQ.2) THEN
>                  Q=T1(IC,IP)
>               ELSEIF (IT.EQ.3) THEN
>                  Q=T2(IC,IP)
>               ELSE
>                  Q=-9999
>               ENDIF
>               IF (Q.LT.-10000) THEN
>C                 degree of freedom
>                  G2=NINT(-Q-10000)
>                  NT=NT+1
>                  IF (NT.GT.MT) GOTO 9000
>                  DT(NT)=G2
>                  CT(NT)=SCD(PI+2+NI+J)
>               ELSE
>                  IQ=NINT(Q)
>                  IF (IQ.EQ.-999.OR.IQ.EQ.-9999) THEN
>C                    unknown variable -> sidecondition cannot be used
>                     NT=0
>                     GOTO 9000
>                  ENDIF
>                  C0=C0+Q*SCD(PI+2+NI+J)
>               ENDIF
>            ENDIF
>110      CONTINUE
>C
>      ELSEIF (TYP.EQ.2) THEN
>C        quadratic sidecondition
>C        -----------------------
>         NI=SCD(PI+1)
>         C0=SCD(PI+3*NI+3)
>         DO 210 J=1,NI
>CC            GT(1)=SCD(PI+1+J)
>CC            GT(2)=SCD(PI+1+J+NI)
>            GT(1)=SCD(PI+2*J)
>            GT(2)=SCD(PI+2*J+1)
>            DO 220 L=1,2
>               G2=GT(L)
>               IF (G2.GT.0) THEN
>C                 degree of freedom
>                  QT(L)=0
>               ELSE
>C                 pointer to tables p,t1,t2 etc.
>                  G2=-G2
>                  IP=G2/100
>                  IT=(G2-100*IP)/10
>                  IC=MOD(G2,10)
>                  IF (IT.EQ.1) THEN
>                     Q=P(IC,IP)
>                  ELSEIF (IT.EQ.2) THEN
>                     Q=T1(IC,IP)
>                  ELSEIF (IT.EQ.3) THEN
>                     Q=T2(IC,IP)
>                  ELSE
>C                    sidecondition cannot be used
>                     NT=0
>                     GOTO 9000
>                  ENDIF
>                  IF (Q.LT.-10000) THEN
>C                    degree of freedom
>                     G2=NINT(-Q-10000)
>                     GT(L)=G2
>                     QT(L)=0
>                  ELSE
>                     GT(L)=0
>                     QT(L)=Q
>                  ENDIF
>               ENDIF
>220         CONTINUE
>            AUX=SCD(PI+2+2*NI+J)
>            IF (GT(1).GT.0.AND.GT(2).GT.0) THEN
>               NT=NT+1
>               L=2*NT
>               DT(L-1)=GT(1)
>               DT(L)=GT(2)
>               CT(NT)=AUX
>            ELSEIF (GT(1).GT.0) THEN
>               NT=NT+1
>               L=2*NT
>               DT(L-1)=GT(1)
>               DT(L)=0
>               CT(NT)=AUX*QT(2)
>            ELSEIF (GT(2).GT.0) THEN
>               NT=NT+1
>               L=2*NT
>               DT(L-1)=GT(2)
>               DT(L)=0
>               CT(NT)=AUX*QT(1)
>            ELSE
>               C0=C0+AUX*QT(1)*QT(2)
>            ENDIF
>210      CONTINUE
>C
>      ELSEIF (TYP.EQ.3.OR.TYP.EQ.33.OR.TYP.EQ.63) THEN
>C        Elimination of extranous inflection points.
>C        -------------------------------------------
>C        There are no inflection points if 1/3 <zpar <2/3.
>C        Format: '3,  segment, coordinate, zpar, weight'
>C                     IP       IC          AUX
>         IP=SCD(PI+1)
>         IF (IP.GE.N) GOTO 9000
>         IC=SCD(PI+2)
>         AUX=SCD(PI+3)
>         IF (NINT(AUX).EQ.-999) THEN
>C           evaluate zval (an approximation)
>            D2=0
>            DK2=K(IP+1)-K(IP)
>            IF (DK2.GT.1E-3) D2=(P(IC,IP+1)-P(IC,IP))/DK2
>            IF (IP.GT.1) THEN
>               D1=0
>               DK1=K(IP)-K(IP-1)
>               IF (DK1.GT.1E-3) D1=(P(IC,IP)-P(IC,IP-1))/DK1
>            ELSE
>               D1=D2
>            ENDIF
>            IF (IP.LT.N-1) THEN
>               D3=0
>               DK3=K(IP+2)-K(IP+1)
>               IF (DK3.GT.1E-3) D3=(P(IC,IP+2)-P(IC,IP+1))/DK3
>            ELSE
>               D3=D2
>            ENDIF
>            QT(1)=P(IC,IP)
>            QT(2)=0.5*(D1+D2)
>            QT(3)=P(IC,IP+1)
>            QT(4)=0.5*(D2+D3)
>            DK=DK2
>            AUX1=DK*QT(2)+QT(1)-QT(3)
>            AUX2=DK*QT(4)+QT(1)-QT(3)
>            AUX=AUX2-AUX1
>            IF (ABS(AUX).GT.1E-5*ABS(AUX2)) THEN
>               AUX=AUX2/AUX
>               IF (ABS(AUX-UL).LT.ABS(AUX-LL)) THEN
>                  AUX=UL
>               ELSE
>                  AUX=LL
>               ENDIF
>               SCD(PI+3)=AUX
>            ELSE
>               AUX=-998
>            ENDIF
>         ENDIF
>CC       IF (AUX.GT.LL.AND.AUX.LT.UL) GOTO 9000
>         IF (NINT(AUX).EQ.-998) GOTO 9000
>         QT(1)=P(IC,IP)
>         QT(2)=T2(IC,IP)
>         QT(3)=P(IC,IP+1)
>         QT(4)=T1(IC,IP+1)
>         DK=K(IP+1)-K(IP)
>         RT(1)=1
>         RT(2)=AUX*DK
>         RT(3)=-1
>         RT(4)=(1-AUX)*DK
>         C0=0
>         DO 310 J=1,4
>            IF (QT(J).LT.-10000) THEN
>               NT=NT+1
>               DT(NT)=NINT(-QT(J)-10000)
>               CT(NT)=RT(J)
>            ELSE
>               C0=C0+QT(J)*RT(J)
>            ENDIF
>310      CONTINUE
>C
>      ELSEIF (TYP.EQ.4.OR.TYP.EQ.34.OR.TYP.EQ.64) THEN
>C        convexity at nodepoint
>C        ---------------------
>C        format: 'typ, point, side, plane, weight'
>C        typ=4: equality requirement
>C        typ=64: apply convexity (option <)
>C        typ=34: break convexity (option >)
>         IP=SCD(PI+1)
>         IF (IP.LT.2.OR.IP.GT.N-1) GOTO 9000
>         IS=SCD(PI+2)
>         IC=SCD(PI+3)
>         I1=AX1(IC)
>         I2=AX2(IC)
>         QT(1)=P(I1,IP-1)
>         QT(2)=P(I2,IP-1)
>         QT(3)=P(I1,IP)
>         QT(4)=P(I2,IP)
>         QT(5)=P(I1,IP+1)
>         QT(6)=P(I2,IP+1)
>         QT(7)=T1(I1,IP)
>         QT(8)=T1(I2,IP)
>         QT(9)=T2(I1,IP)
>         QT(10)=T2(I2,IP)
>         DO 410 J=1,6
>            IF (QT(J).LT.-10000) THEN
>C              degree of freedom
>               GOTO 9000
>            ENDIF
>410      CONTINUE
>         DU(1)=QT(3)-QT(1)
>         DV(1)=QT(4)-QT(2)
>         DU(2)=QT(5)-QT(3)
>         DV(2)=QT(6)-QT(4)
>         AUX=DU(1)*DV(2)-DV(1)*DU(2)
>         D1=SQRT(DU(1)*DU(1)+DV(1)*DV(1))
>         D2=SQRT(DU(2)*DU(2)+DV(2)*DV(2))
>         IF (IS.EQ.1) THEN
>C           before the point
>            I1=7
>            I2=8
>         ELSE
>C           after the point
>            I1=9
>            I2=10
>         ENDIF
>         IF (ABS(AUX).LT.D1*D2*1E-3.and..false.) THEN
>C           zero angle
>         ELSE
>C           nonzero angle
>            IF (QT(I1).LT.-10000) THEN
>               NT=NT+1
>               DT(NT)=NINT(-QT(I1)-10000)
>               CT(NT)=-DV(IS)
>            ELSE
>               C0=C0-DV(IS)*QT(I1)
>            ENDIF
>            IF (QT(I2).LT.-10000) THEN
>               NT=NT+1
>               DT(NT)=NINT(-QT(I2)-10000)
>               CT(NT)=DU(IS)
>            ELSE
>               C0=C0+DU(IS)*QT(I2)
>            ENDIF
>            LGI=AUX.LT.0.AND.IS.EQ.1.OR.AUX.GT.0.AND.IS.EQ.2
>            IF (LGI.AND.TYP.EQ.34 .OR. .NOT.LGI.AND.TYP.EQ.64) THEN
>               DO 420 J=1,NT
>                  CT(J)=-CT(J)
>420            CONTINUE
>               C0=-C0
>            ENDIF
>         ENDIF
>C
>      ELSEIF (TYP.EQ.5.OR.TYP.EQ.35.OR.TYP.EQ.65) THEN
>C        inequality constraint iabs(case)=1 of gb812
>C        -------------------------------------------
>C        '5, point, side, plane, weight'
>         IP=SCD(PI+1)
>         IS=SCD(PI+2)
>         IC=SCD(PI+3)
>         I1=AX1(IC)
>         I2=AX2(IC)
>         IF (IS.EQ.1) THEN
>            IP1=IP-1
>            CALL MA1811(P(1,IP1),P(1,IP),A,IC,IER)
>            IF (IER.GT.0) GOTO 9000
>         ELSE
>            IP1=IP+1
>            CALL MA1811(P(1,IP),P(1,IP1),A,IC,IER)
>            IF (IER.GT.0) GOTO 9000
>         ENDIF
>         CALL MA1810(P(1,IP-1),P(1,IP),P(1,IP+1),DA,IC,IER)
>         IF (IER.GT.0) GOTO 9000
>         CALL MA1810(P(1,IP1-1),P(1,IP1),P(1,IP1+1),DA1,IC,IER)
>         IF (IER.GT.0) GOTO 9000
>C        a: angle of a segment
>C        da: angle between connecting chords at ip
>C        da1: angle between connecting chords at ip1
>         A=A+QC*ABS(DA1)
>         D1=COS(A)
>         D2=SIN(A)
>         IF (IS.EQ.1) THEN
>            Q1=T1(I1,IP)
>            Q2=T1(I2,IP)
>         ELSE
>            Q1=T2(I1,IP)
>            Q2=T2(I2,IP)
>         ENDIF
>         IF (Q1.LT.-10000) THEN
>            NT=NT+1
>            DT(NT)=NINT(-Q1-10000)
>            CT(NT)=-D2
>         ELSE
>            C0=C0-D2*Q1
>         ENDIF
>         IF (Q2.LT.-10000) THEN
>            NT=NT+1
>            DT(NT)=NINT(-Q2-10000)
>            CT(NT)=D1
>         ELSE
>            C0=C0+D1*Q2
>         ENDIF
>         LGI=DA.GT.0.AND.IS.EQ.1 .OR. DA.LT.0.AND.IS.EQ.2
>         IF (LGI.AND.TYP.EQ.35 .OR. .NOT.LGI.AND.TYP.EQ.65) THEN
>            DO 510 J=1,NT
>               CT(J)=-CT(J)
>510         CONTINUE
>            C0=-C0
>         ENDIF
>C
>      ELSEIF (TYP.EQ.6.OR.TYP.EQ.36.OR.TYP.EQ.66) THEN
>C        inequality constraint iabs(case)=2 of gb812
>C        -------------------------------------------
>C        '6, point, side, plane, weight'
>         IP=SCD(PI+1)
>         IS=SCD(PI+2)
>         IC=SCD(PI+3)
>         I1=AX1(IC)
>         I2=AX2(IC)
>         IF (IS.EQ.1) THEN
>            IP1=IP-1
>            CALL MA1811(P(1,IP1-1),P(1,IP1),A,IC,IER)
>            IF (IER.GT.0) GOTO 9000
>         ELSE
>            IP1=IP+1
>            CALL MA1811(P(1,IP1),P(1,IP1+1),A,IC,IER)
>            IF (IER.GT.0) GOTO 9000
>         ENDIF
>         CALL MA1810(P(1,IP-1),P(1,IP),P(1,IP+1),DA,IC,IER)
>         IF (IER.GT.0) GOTO 9000
>         CALL MA1810(P(1,IP1-1),P(1,IP1),P(1,IP1+1),DA1,IC,IER)
>         IF (IER.GT.0) GOTO 9000
>C        a: angle of a segment
>C        da: angle between connecting chords at ip
>C        da1: angle between connecting chords at ip1
>         A=A+QC*DA1
>         D1=COS(A)
>         D2=SIN(A)
>         IF (IS.EQ.1) THEN
>            Q1=T2(I1,IP1)
>            Q2=T2(I2,IP1)
>         ELSE
>            Q1=T1(I1,IP1)
>            Q2=T1(I2,IP1)
>         ENDIF
>         IF (Q1.LT.-10000) THEN
>            NT=NT+1
>            DT(NT)=NINT(-Q1-10000)
>            CT(NT)=-D2
>         ELSE
>            C0=C0-D2*Q1
>         ENDIF
>         IF (Q2.LT.-10000) THEN
>            NT=NT+1
>            DT(NT)=NINT(-Q2-10000)
>            CT(NT)=D1
>         ELSE
>            C0=C0+D1*Q2
>         ENDIF
>         LGI=DA.GT.0.AND.IS.EQ.1 .OR. DA.LT.0.AND.IS.EQ.2
>         IF (LGI.AND.TYP.EQ.66 .OR. .NOT.LGI.AND.TYP.EQ.36) THEN
>            DO 610 J=1,NT
>               CT(J)=-CT(J)
>610         CONTINUE
>            C0=-C0
>         ENDIF
>C
>      ELSEIF (TYP.EQ.7.OR.TYP.EQ.37.OR.TYP.EQ.67) THEN
>C        curvature sidecondition
>C        -----------------------
>C        '7, segment, parameter, radius, weight'
>C        ...
>C
>      ENDIF
>C
>9000  CONTINUE
>      END

Actions: View

Attachments on bug 73707: 104674