Треугольный элемент. Функции формы для элементов высокого порядка. Вычисление производных функций формы. Программа MCHB, страница 3

     SUBROUTINE MCHB(R,A,M,N,MUD,IOP,EPS,IER)

C

C

      DIMENSION R(1),A(1)

      DOUBLE PRECISION TOL,SUM,PIV

C

C        TEST ON WRONG INPUT PARAMETERS

      IF(IABS(IOP)-3)1,1,43

    1 IF(MUD)43,2,2

    2 MC=MUD+1

      IF(M-MC)43,3,3

    3 MR=M-MUD

      IER=0

C

C        MC IS THE MAXIMUM NUMBER OF ELEMENTS IN THE ROWS OF ARRAY A

C        MR IS THE INDEX OF THE LAST ROW IN ARRAY A WITH MC ELEMENTS

C

C     ******************************************************************

C

C        START FACTORIZATION OF MATRIX A

      IF(IOP)24,4,4

    4 IEND=0

      LLDST=MUD

      DO 23 K=1,M

      IST=IEND+1

      IEND=IST+MUD

      J=K-MR

      IF(J)6,6,5

    5 IEND=IEND-J

    6 IF(J-1)8,8,7

    7 LLDST=LLDST-1

    8 LMAX=MUD

      J=MC-K

      IF(J)10,10,9

    9 LMAX=LMAX-J

   10 ID=0

      TOL=A(IST)*EPS

C

C        START FACTORIZATION-LOOP OVER K-TH ROW

      DO 23 I=IST,IEND

      SUM=0.D0

      IF(LMAX)14,14,11

C

C        PREPARE INNER LOOP

   11 LL=IST

      LLD=LLDST

C

C        START INNER LOOP

      DO 13 L=1,LMAX

      LL=LL-LLD

      LLL=LL+ID

      SUM=SUM+A(LL)*A(LLL)

      IF(LLD-MUD)12,13,13

   12 LLD=LLD+1

   13 CONTINUE

C        END OF INNER LOOP

C

C        TRANSFORM ELEMENT A(I)

   14 SUM=DBLE(A(I))-SUM

      IF(I-IST)15,15,20

C

C        A(I) IS DIAGONAL ELEMENT. ERROR TEST.

   15 IF(SUM)43,43,16

C

C        TEST ON LOSS OF SIGNIFICANT DIGITS AND WARNING

   16 IF(SUM-TOL)17,17,19

   17 IF(IER)18,18,19

   18 IER=K-1

C

C        COMPUTATION OF PIVOT ELEMENT

   19 PIV=DSQRT(SUM)

      A(I)=PIV

      PIV=1.D0/PIV

      GO TO 21

C

C        A(I) IS NOT DIAGONAL ELEMENT

   20 A(I)=SUM*PIV

C

C        UPDATE ID AND LMAX

   21 ID=ID+1

      IF(ID-J)23,23,22

   22 LMAX=LMAX-1

   23 CONTINUE

C

C        END OF FACTORIZATION-LOOP OVER K-TH ROW

C        END OF FACTORIZATION OF MATRIX A

C

C     ******************************************************************

C

C        PREPARE MATRIX DIVISIONS

      IF(IOP)24,44,24

   24 ID=N*M

      IEND=IABS(IOP)-2

      IF(IEND)25,35,25

C

C     ******************************************************************

C

C        START DIVISION BY TRANSPOSE OF MATRIX TU (TU IS STORED IN

C        LOCATIONS OF A)

   25 IST=1

      LMAX=0

      J=-MR

      LLDST=MUD

      DO 34 K=1,M

      PIV=A(IST)

      IF(PIV)26,43,26

   26 PIV=1.D0/PIV

C

C        START BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R

      DO 30 I=K,ID,M

      SUM=0.D0

      IF(LMAX)30,30,27

C

C        PREPARE INNER LOOP

   27 LL=IST

      LLL=I

      LLD=LLDST

C

C        START INNER LOOP

      DO 29 L=1,LMAX

      LL=LL-LLD

      LLL=LLL-1

      SUM=SUM+A(LL)*R(LLL)

      IF(LLD-MUD)28,29,29

   28 LLD=LLD+1

   29 CONTINUE

C        END OF INNER LOOP

C

C        TRANSFORM ELEMENT R(I)

   30 R(I)=PIV*(DBLE(R(I))-SUM)

C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R

C

C        UPDATE PARAMETERS LMAX, IST AND LLDST

      IF(MC-K)32,32,31

   31 LMAX=K

   32 IST=IST+MC

      J=J+1

      IF(J)34,34,33

   33 IST=IST-J

      LLDST=LLDST-1

   34 CONTINUE

C

C        END OF DIVISION BY TRANSPOSE OF MATRIX TU

C

C     ******************************************************************

C

C        START DIVISION BY MATRIX TU (TU IS STORED ON LOCATIONS OF A)

      IF(IEND)35,35,44

   35 IST=M+(MUD*(M+M-MC))/2+1

      LMAX=0

      K=M

   36 IEND=IST-1

      IST=IEND-LMAX

      PIV=A(IST)

      IF(PIV)37,43,37

   37 PIV=1.D0/PIV

      L=IST+1

C

C        START BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R

      DO 40 I=K,ID,M

      SUM=0.D0

      IF(LMAX)40,40,38

   38 LLL=I

C

C        START INNER LOOP

      DO 39 LL=L,IEND

      LLL=LLL+1

   39 SUM=SUM+A(LL)*R(LLL)

C        END OF INNER LOOP

C

C        TRANSFORM ELEMENT R(I)

   40 R(I)=PIV*(DBLE(R(I))-SUM)

C        END OF BACKSUBSTITUTION-LOOP FOR K-TH ROW OF MATRIX R

C

C

C        UPDATE PARAMETERS LMAX AND K

      IF(K-MR)42,42,41

   41 LMAX=LMAX+1

   42 K=K-1

      IF(K)44,44,36

C

C        END OF DIVISION BY MATRIX TU

C

C     ******************************************************************

C

C        ERROR EXIT IN CASE OF WRONG INPUT PARAMETERS OR PIVOT ELEMENT

C        LESS THAN OR EQUAL TO ZERO

   43 IER=-1

   44 RETURN

      END

Литература

1. Сегерлинд , Ларри Дж  “Применение метода конечных элементов” М: “Мир”, 1997г.

2. Бреббиа и Конор  “Метод конечных элементов для решения задач жесткости”