# Треугольный элемент. Функции формы для элементов высокого порядка. Вычисление производных функций формы. Программа 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. Бреббиа и Конор  “Метод конечных элементов для решения задач жесткости”