Efficient calculation of finite element stiffness matrices

Efficient calculation of finite element stiffness matrices I. R. WHITE, K. M O R G A N and R. W. L E W I S Department of Civil Engineering, University of Wales, Singleton Park, Swansea SA2 8PP, UK (Received 14 November 1980)

A method of increasing the efficiency of the construction of stiffness matrices in the finite element method is described. It considers how the efficiency of coding is affected by array referencing and shows how the restricted use of multisubscript arrays can improve the c.p.u, costs. Use is made of the large core storage areas usually available within finite element packages, together with the optimal use of backing store to keep cartesian derivative information which is sequentially loaded into buffers. An estimate of the efficiency of this method is also given which provides a criterion for its use when the global derivatives are known to be required on a given number of occasions during the core life of a program.

from the disadvantage that the numbering of the nodes comprising the mesh is crucial to the efficiency of the solver, and this node numbering becomes difficult in large two- and three-dimensional meshes representing complex structures. This difficulty was overcome in part with the introduction of the 'front' solver ''6 where the machine considered each element of the mesh in turn and performs the forward elimination on a given node as soon as all the relevant details concerning it have been constructed. The solver, however, is not the only portion of a finite element program where savings may be made in machine time. In this paper the authors consider the efficiency of the construction routines often used for calculating the element matrices before they are fed into the solver. CURRENT M E T H O D S O F ANALYSIS Consider the solution of the equilibrium equations:

INTRODUCTION In recent years the increased use of the finite element method has led to a proliferation of suites of software to solve most of the commonly encountered differential equations of engineering and applied mathematics. The consideration of linear problems quickly gave way to the analysis of complex non-linear situations where the quantity of computational effort was greatly increased. There are in the literature some well-tested and documented programs for solving both linear 1'4 and non-linear 2"3 problems: and it is now widely recognised that the programming effort required to implement finite element techniques of analysis is considerably greater than rival methods such as finite differences. These numerical methods are by their nature computer intensive requiring large amounts of machine time and considerable use of on-line peripherals. Some authors have already considered the efficiency of various parts of the programming method. Techniques of solving the systems of linear equations produced by the analyses have been considered in some depth l's'6's, and a wide variety of'solvers' is now available to the analyst when constructing a software system. The original total matrix solvers where the whole matrix is stored in core were quickly observed to be inefficient on core usage and gave way to solvers where the matrix was partitioned into blocks and stored in these blocks on disc or tape. The banded structure of most matrices produced by the finite element method can easily be used to improve the efficiency of the solver where all the non-zero members of the matrix are known to lie within a given bandwidth s. These 'banded' solvers, however, suffer 0141- I 195/81/020077-07$2.00 ©1981 CML Publications

aij.j = .[i

in an N-dimensional region fl, where a~j is the stress tensor at a point,f~ is the body force in direction i and the subscript following the comma denotes the direction of differentiation. Now if u~ is the displacement in the ith direction, N m the shape function associated with node m and uT' is the displacement in direction i at node m then the strain tensor is given by:

1

~i j = ~( U i.j 21- li j,i)

(1)

and the displacement field is given by: ui = Nmum Substituting into equation (1) the expression: --

m

m

~ij -- ~ijkLlk

(2)

is obtained where:

and/~j is the Kronecker delta. The stresses and strains at a point in the interior of f~ are linked by the material property tensor ~' to give % = g'ij, lek~

Ado. Eng. Software, 1981, Vol. 3, No. 2

(3)

77

The finite element analysis of the behaviour of f~ under various loading conditions rests on the substitution of equations (2) and (3) into the equilibrium equations to obtain the following:

buffer to store the global derivatives used in the calculation of , ~ .

,~/" i~"u"i= L'f

When a Fortran array is referenced from an object code the position of the data within the core is calculated by the machine from a knowledge of the subscript values and the sizes of the array as given in its declarator (usually a D I M E N S I O N or C O M M O N ) . As a result the time taken to perform a given numerical task may not be a function solely of the operations count of that task. As an example consider the simple task of multiplying together two square matrices B and C and storing the result in A. If the matrices are all n × n in size this task requires n 2 multiplications and n 2 additions. To perform this task as it would be done on paper would result in the code given below:

(4)

where L m represents the loads in direction i applied to node m. ,~g'7~"is a tensor, usually referred to as the stiffness matrix, and is given by:

-.-f U --

~uklN~Nn.I

d~

(5)

fl

The tensor ~auu.is, for linear elastic material behaviour, comprised only of constants. However, this paper is concerned with the evaluation of fft"~"~"where ~Uk~ is non-constant. In these situations the stiffness matrix j,(-7j, is evaluated a large number of times due to the fact that g'Uk~ varies with strain, strain rate, stress etc. 9"1°. A variety of iterative solution methods have been developed for the solving of equation (4) from a given initial estimate of u~. In many of these methods the repeated calculation of J ( ' ~ " takes u'p a considerable portion of the run-time with the numerical integration of equation (5). This means that the rapid calculation of the integrand in equation (5) is of crucial importance to the efficiency of the program. Even in a linear problem the evaluation of°//"m" •..,.. ij may take up a significant proportion of the program's c.p.u. usage. This is mainly due to the fact that the integral in equation (5) is evaluated numerically by mapping each element e (e = 1..... %) onto a prototype element of the same type and of a standard size positioned at or near the origin. Thus the quantity actually computed in equation (5) is:

,yUU ={ Inn

.... l I ,~uuN.kN.,,,,{°MIj,,,

M U L T I S U B S C R I P T ARRAYS

SUBROUTINETRYI (NSIZEp MSIZEr AMATX~BMATx,C H A T X ) DIMENSIONAMATX(MSIZE,MSIZE), BMATX(MSIZE,MSIZE), CHATX(MSIZE,MSIZE)

DO 30 IROWS : I,NSIZE DO 20 ICOLS = ItN$IZE AHATX(IROWS~ICOLS) = 0.0 DO 10 KTERM = 1,NSIZE

AHATX(IROWStICOLS) = AHATX(IROWS,ICOLS) + BHATX(IROWS,KTERH) • CHATX(KTERM,ICOLS) 10 "CONTINU E 20 CONTINUE 30 CONTINUE RETURN END

In this method each element of the matrix A is initialized at zero and the products from B and C are added in term by term. This code, in addition to the necessary n 2 multiplications and n 2 additions, contains 4n 3 + n 2 references to 2-subscript arrays. The c.p.u, time taken to execute this code is shown as method 1 in Fig. 2 as a function of n. The quantity of array referencing may be reduced by declaring a vector W O R K A as a work area in which one row of the matrix B is stored.

where o~ is the integration weight applied to sampling point g (g = 1..... %) in element e and IJI is the determinant of the Jacobian of the transform from element e to its prototype evaluated at the sampling point. In a program which solves only linear problems the calculation offft"Tj" needs to be carried out once only but when the material property tensor varies the amount of computation involves is significantly increased. This repeated calculation of ~ involves re-evaluating {N.~}g e and {oglJl}g, repeatedly. The derivative buffering method described below aims at reducing this computational load to a minimum. The reduction is achieved by: (a) multiple use of existing core storage areas within the program, (b) the avoidance of extensive use of multisubscript arrays, (c) optimal use of backing store.

Shape function Evaluation

Numerical integration orion of global dm'ivofives of J~, on disk "

Adv. Eng. Software, 1981, Vol. 3, No. 2

Cek:ulofion of [D][B]

STFFNESS HATRIX J~

M U L T I P L E USE O F S T O R A G E AREAS

78

~ Y~

SETTING UP

Forwerd Elimination

The first of these factors stems from the observation that most finite element programs contain a relatively large array within the'solve routine which is used as a work area. When control is not in the solver this work area remains untouched and contains information superfluous to the current stage of the program. The derivative buffering method described in this study uses this area as a

/Cofculation of / Jocobions and / their inverses

mab'lces

CALCULATION

FRONT SOLVER ~ L

_

~

CONSIRUCTIONOF LOAO VECTOR

DATA INPUT

INITIAL PROCESSING

Figure 1.

-

/

/

°:I

--o.--

Method 1: 93W Method 2: 107W Method 3: 102W

METH09 2

The percentage increase in speed of the last two methods over the first is shown in Fig. 3. From this it can be seen that the improvement is roughly constant as n increases and that the avoidance of extensive use of multisubscript arrays has resulted in a doubling of the speed of execution of this simple task. DERIVATIVE B U F F E R I N G M E T H O D 5 a.' -

'Jo-

i

O. OORDER OF HATRICES I

20.

-~'0

,'o.

5'0

Figure 2.

10

2O 3O

SUBROUTINE TRY2 (NSIZE, MSIZE, N4ATX, BHATX, CHATX, WORKA) DIMENSION AHATX(HSIZE,MSIZE), BHATX(MSIZE,MSIZE), CMATX(MSIZE, MSIZE), WORKA(MSIZE) "DO 30 IROWS = I,NSIZE DO 10 KTERM : I,NSIZE WORKA(KTERM) : BMATX(IROWS,KTERM) DO 30 ICOLS : 1, NSIZE SAVED = 0.0 DO 20 KTERM = 1, NSIZE Q4A SAVED : SAVED + WORKA(KTERM)* TX(KTERM,ICOLS) AMATX(IROWS,ICOLS) = SAVED RETURN END

The addition of the term is then performed in a temporary location SAVED to be placed in A when complete. This second method, although faster than the first, suffers from the disadvantage that it has a larger core requirement due to the presence of the work area. However, the improvement in speed can be seen from Fig 2. The third method shown below makes use of the oi'dering in which the matrices B and C are presented to the subroutine. A, B and C are declared as vectors of unit length and the machine calculates explicitly the address of each piece of data in the array as it is required. The declaring of A, B and C as vectors is made possible in Fortran due to the fact that the data are passed to the subroutine by location and not by value ? and this optimization obviates the need to reference any of the arrays as 2-subscript arrays.

The derivative buffering method stems from the observation that, using conventional techniques, the terms of [ N ~ g e and Ucoldl]ge in equation (6) are recalculated o n each occasion that they are required. The proportion of the total run-time of a finite element program that the calculation of ~ ( o c c u p i e s can be seen from Fig. 1. The program used here was that published in ref. 1 (with some minor modifications) and it can be seen that the calculation of ,~"occupies nearly half of the total run-time and uses more than twice the c.p.u, time of the front solver. The problem solved to produce Fig. 1 consisted of twelve quadrilateral elements in an array producing a frontwidth of 16. It should be recognised that in larger problems involving hundreds of elements the solver will take a larger proportion of the total run time due to the much larger frontwidth. The method described here takes the values of I-N~g e and I~olJl-lgeand stares them optimally on backing store. These values are then accessed only when they are required and set up in a way which avoids extensive use of multisubscript arrays but leaves them in an ordering that is conventional for many finite element packages. The data handling is achieved by constructing a series of 3-level data trees as shown in Fig. 4, the third level of which comprises the buffer to be written onto backing store. Each buffer contains data relating to a small patch of elements, the size of the patch being determined by the core space available for a work area. Every element in a patch is indexed incrementally and this patch index number indicates the position in the block indices of the control parameters for the element under consideration. For the ith element in a patch the control parameters consist of d(i), the number of coordinate dimensions in the

SUBROUTINE TRY3 (NSIZE, MSIZE r AMATX, BHATX, CHATX) DIMENSION AMATX(I ), BHATX( I ), CMATX( I ) IASTA = 1 ICSTA = 1 DO 3 0 JCOLS = 1,NSIZE IAPOS : IASTA DO 20 IROWS = I,NSIZE IBPOS = IROWS

J _<

---0---0--

ME1HO0 2 .EIHO0 3

ICFIN = ICSTA + NSIZE - 1 SAVE]) = 0.0

10 20 30

DO 10 ICPOS = ICSTA,ICFIN SAVED = SAVED + BHATX(IBPOS)*CMATX(ICPOS) IBPOS = IBPOS + MSIZE AMATX(IAPOS) = SAVED IAPOS = IAP(~ + I IASTA : IASTA + MSIZE ICSTA = ICSTA + HSIZE RETURN END

go" o - . . . ~ o

o"

~

..o.____.o-~o

o.

o

z

:

:

o

o

:

:

o

o"

Although long in terms of lines of code the decrease in c.p.u, usage over the previous two methods can again be seen from Fig. 2. The increase in object code size is small. On an ICL 1904S these sizes are as follows:

zg"

,'0

ORQER 0 ¢ MATRICES

~'0

9'0

,'0

~0

Figure 3.

Adv. Eng. Software, 1981, Vol. 3, No. 2 79

I, ~w~rnent index h;n I~tc,h

-:-1 element mmt~-

i i

Figure 4. element, h(i), the number of nodes in the element and s(i), the number of integration sampling points in the element. Clearly these parameters are bounded above by maxima determined by the dimensions of the program. Ifm d and m,~ are the upper bounds of d(k) and if(k) for elements k in the patch then the position in the buffer of the data for element i can be evaluated from these parameters. The position P~ of the first member of the [~olJI],~ list is given by:

parts of the method. The picture is further complicated by the use of backing store albeit in as economical a fashion as possible. As a consequence, the question of whether this method is computationally cheaper in terms of c.p.u, time use is machine dependent and can only be accurately assessed after the user has carried out some preliminary experimentation on the behaviour of a few simple machine operations. In order to test its efficiency the authors have derived a predictive method of assessing whether or not derivative buffering is an improvement in terms of c.p.u, usage over the procedure of evaluating the global derivatives on each occasion they are required within the program. The fact that the software for this method is capable of analysing mixed meshes in arbitrary derivative buffer lengths means that the c.p.u, time taken to evaluate N'.7 may vary from element to element. Let te. k be the c.p.u. time taken to evaluate N ~.,,in element k, tw.k be the time to write the kth buffer onto backing store (of whatever medium), tr,. be the time to execute a rewind, t,. k be the time to read block k from backing store and t~.k be the time to set up the derivatives from the buffer into the relevant arrays for use within the program. Then if there are n elements requiring n b buffers we may define the quantity y by: tlb

i=1

i-I

P1 = 1 + mdmr, ~ 's(k) k=l

and the position P2 of the first member of the cartesian derivative list is given by:

P2=s(i)+Pt

k=l

It can be shown (see Appendix A) that the method is always less efficient if y ~<0 regardless of the number of times the global derivatives are required by the program. However, when y > 0, as was the case in all of the author's tests, the method is more efficient if m, the number of times the global derivatives are required, exceeds a critical value mc~, (not necessarily an integer) given by:

where

'f(k) = k=l

f 0

, ~ k~

n

I=0

n

/>0 In the 12-element example considered above the value of Incrit using an ICL 1904S mainframe and disc backing store was 3.62. So the derivative buffering method is more efficient than explicit calculation if the global derivatives are required on four or more occasions during the run of the entire program. This is certainly the case in all but the simplest non-linear analyses.

The values of P~ and P2 may be calculated incrementally within the program as the data are transferred to the buffer. When an element is encountered that would overload the buffer (the broken line in Fig. 4) the entire data tree is written on to backing store and the construction of another is commenced. By processing the data in this way transfers to and from backing store are reduced to a minimum for a given core work area. This work area should be made as large as possible within the bounds of the program and is normally overlaid on the 'global stiffness array' or another large work area used within the solver. At the end of the construction phase the disc file created by the rountines contains the three block indices and one buffer on each record. The buffer is comprised of lists of the weighted Jacobians, wlJlg e, and the global derivatives stored sequentially with the members of the block indices determining the position of the data for a given element as described above.

A method has been described which enables finite element stiffn6ss matrices to be evaluated with greater efficiency. This is achieved by buffering values of N~ as part of a data tree and storing a sequence of these trees optimally on backing store. The effects of multiple dimension array referencing on efficiency have been examined in order to optimize both the in-core data structure and c.p.u, usage. Practical efficiency bounds have also been derived so that a prospective user may assess whether or not this method improves on conventional techniques in a particular case.

PRACTICAL EFFICIENCY BOUNDS

APPENDIX A

When using techniques such as the one under consideration here the efficiency cannot be assessed purely by adding the operations counts for the various constituent

In comparing the c.p.u, times taken by the two methods it is assumed that the number of occasions the global derivatives are required by the program is m. Thus for

80

Adv. Eng. Software, 1981, Vol. 3, No. 2

CONCLUSIONS

increased efficiency we must have: & m ~ t~.;> (time to set up N ".~, sort into buffers and write i=1

onto backing store)+ m (time to read and dismantle the buffers) which leads to nb

m r > ~.t~.i+ E tw.k>O i=1

k=l

where y is as defined in the text. Consequently if 7 < 0 we obtain m 0 and so efficiency was improved provided m exceeded /~Zcrit. APPENDIX B The coding for the derivative buffering technique is given below. The data are passed to the processing segments through argument lists to enable these subroutines to be set up in program libraries and thus avoid repeated code handling and compilation costs. A description of the contents of the arrays and the other subroutines accessed is also given to enable the user to incorporate this technique into existing finite element packages.

Cartesian derivative evaluation The evaluation of the cartesian derivatives throughout an element of any type in up to three dimensions is performed by the subroutine CARTD2. It evaluates the terms of tensors N ~. number within ,1 where r e = n o d e element, i = direction of integration. This is carried out by solving the equations:

ment when sampled at Gauss point IGAUS is stored in DERIV (IDIME, INPEL, IGAUS). This array is vectored within CARTD2 and should be assembled before the subroutine is entered. ELCOD Real array with dimensions (MDIME, MNPEL). This array contains the coordinates of the nodes within the element. The coordinate in direction I D I M E of node I N P E L should be pre-stored in E L C O D (IDIME, INPEL) before the subroutine is called. TOLER The value of the determinant tolerance ~. If the determinant of J lies in the range (-e,e) it is treated as singular and the program is halted. NDIME The number of coordinate dimensions this element contains. MDIME The maximum number ofcoordinate dimensions allowed within the problem. NGAUS The number of Gauss points used in this element. MGAUS The maximum number of Gauss points allowed within the problem. NNPEL The number of nodes in this element. MNPEL The maximum number of nodes allowed in any element within the problem. MCART The maximum number of cartesian derivatives to be calculated. This should be set to the size of the CARTD array and is thus given by M D I M E * M N P E L , M G A U S . The function of the CARTD2 subroutine is to set up the arrays CARTD and DETJA for use either in CARTD3 or elsewhere in the program. All the other data must be set up prior to calling the routine. C l C A ] ~ I ' D 2 | l | | | t l l l l l l I l l J J|lttl~ll l l l l l l l l l l l | l l l l i l

DETJA

DERIV

Real array with dimensions (MDIME, M N P E L , MGAUS). The global derivative in the ith coordinate direction of the shape function associated with n o d e j in an element sampled at sampling point (or 'Gauss point') k is assembled in CARTD (ij, k). This array is vectored within the subroutines. Real array with dimensions (MGAUS). The determinant of the Jacobian at each Gauss point is assembled in this array. Real array with dimensions (MDIME, M N P E L , MGAUS). The local derivatives within the prototype element are stored in this array. The derivative in the coordinate direction I D I M E of the shape function associated with node I N P E L within the ele-

Jl|i~

J ! |MiUllltlll

F~)IMEI NGAUS, HGAUSI NNPEL, HNPEL, MCART) C C C C C C C

C C

SUBROUTINE TO CALCULATE THE CARTESIAN DERIVATIVES AT ALL GAUSS POINTS WITHIN AN ELE~4ENT IAN WHITE

IKrRODUCED AT IANLIB VERSION cR6q (JULY 1980) UPGRADED AT VERSION GO75 (AUGUST 1980) UPGRADED AT VERSION G082 (AUGUST 1980) UPGRADED AT VERSION GO84 (AUGUST 1980) LOGICAL SINGA REAL JACOB(313) DIHENSION CARTD(HCANT), DERIV(HCANT), ELCOD(PDIHE, HNPEL), DETJA(HGAUS) EOUIVALENCE (JACOB( 1 , 1 ) , XJA11 ) , (JACOB(2,1 ) , XJ~21 ) j (JACOB( 3 , 1 ) , XJA31 ) (JACOB(1,2),XJA12), (JACOB(2,2),XJA22), (JACOB(3,2),XJA32), (JACOB(1,3),XJA13), (JACOB(2,3),XJA23), (JACOB(3,3),XJA33) COHHON/IANLIB/LCHAN ICART = 1

Parameter list CARTD

! li,llli

SUBROUTINE CARTD2 (CARTD, DETJA, DERIV, ELCOD, TOLER I NDIMEI

C

at each Guass point..N~ are the local derivatives of N ~ in the prototype element. The determinant of the Jacobian J is also evaluated and tested to ensure that it is bounded away from zero within the tolerance of the machine.

li,

C

C C

SCAN THE GAUSS POINTS AND EVALUATE THE JACOBIAN

C DO 90 IGAUS = I, NGAUS DO 20 IDIHE = I,NDIHE

10 20 C

C C

DO 20 JDIHE = 1,NDIHE XJACO = O. DO 10 INPEL = 1,NNPEL IDERI = IDIHE + FQ)IHE='(INPEL-I+HNPELZ'(IGAUS-1)) XJACO = XJACO + DERIV(IDERI)*ELCOD(JDIHE,INPEL) JACOB(IDIHE,JDIHE) = XJACO

EVALUATE THE DETERHINANT OF d AND TEST FOR NEAR SINGULARITY IF (NDIHE. NE. I ) GOTO 30 DETJB = XJA11 GOTO 5O

30

IF (NDIHE.NE.2) GOTO qO DETJB = XJA11'=XJA22 - XJA211XJA12 GOTO 50

qo

IF (NDIHE.NE.3) CALL ERROR ('CARTD2*'p 11 LCHAN1 .TRUE.) DETJB = XJA11*(XJA22*XJA33 - XJA23*XJA32) - XJA12*(XJA21*XJA33 - XJA23*XJA31) + XJA13*(XJA21'~XJA32 - XJA22*XJA31)

Adv. Eng. Software, 1981, Vol. 3, No. 2 81

50 52 53 55 C C C

IF (ADS(DRTJB).GT.TOLER) GOTO55 IF (LCHAN.LE.O) GOTO53 DO 52 IDIME = I,NDIME WRITE (LCHAN,IO00) (JACOB(IDIME,JDIHE), JDIHE = I,NDIME) WRITE (LCHAN,1000) DETJB, TOLER CAll ERROR('CARTD2i ' , 2, LCHAN, .TRUE.) DETJA(IGAUS) -- DETJB

SOLVE FOR CARTESIAN DERIVATIVES

60

DETJI = I./DETJB DO 80 INPEL = I,NNPEL IF (NDIME. NE. I) GOTO 60 CARTD(ICART) = DERIV(ICART)IDETJI GOTO 80 JCART : ICART + I DERII = DERIV(ICART) DERI2 = DERIV(JCART) IF (NDIME.NE.2) GOTO 70 CARTD(ICART) = ( D E R I I m X J A 2 2 - DERI2aXJA12)WDETJI CARTD(JCART) = (DERI21XJA11 -DERIIUXJA21)IDETJI

70

GOTO 80 KCART= JCART + I DERI3 = DERIV(KCART)

80 90

SDETI : XJA22~DERI3 - XJA32*DERI2 SDET2 = XJA21mDERI3 - XJA31~DERI2 SDET3 = XJA23mDERI3 - XJA33tDERI2 DIX33 : DERII"XJA33 DIX32 = DERIImXJA32 DIX31 : DERIImXJA31 CARTD(ICART) = ( X J A 1 2 J S D E T 3 - XJAI3~SDET1 • + DIX33~XJA22 - DIX321XJA23)UDETJI CARTD(JCART) = ( X J A 1 3 ~ S D E T 2 - XJA11~SDET3 • + DIX31mXJA23 - DIX33tXJA21)IDETJI CARTD(KCART) = (XJA11~SDET1 - XJAI2~SDET2 + DIX32~XJA21 - DIX31~XJA22)~DETJI ICART = ICART + ~O)IM£ ICART = ICART + (HNPEL - NNPEL)IPDIME

RETURN 10OO FORMAT (3G15.7) END

Buffer assembly and storage The assembly and storage of the buffers is performed in subroutine CARTD3. This process is carried out by setting up the shape function values and local derivatives in the arrays SHAPE and DERIV and then evaluating the global derivatives with a call to CARTD2. The block indices are then assembled and an assessment is made on the current state of the buffer. If there is sufficient space to store CARTD then no writing to disc will be performed and another element will be added to the current patch. If, however, the buffer is already full it will be written on to backing store and a new buffer commenced with the derivatives of the current element.

Parameterlist COORD NCONN

NGCOD

JCHAN KCHAN MDIME MELEB NELEM MELEM MGAUS

82

Real array with dimensions (MDIME, MNODE) containing the global coordinates of each node. Integer array with dimensions (MNPEL, MELEM) which contains the nodal connectivity. The number of node I N P E L in element IELEM should be pre-stored in N C O N N ( I N P E L , TELEM). Integer array of size (MELEM) containing the Gauss code for each element. This Gauss code determines the element type and Gaussian integration rule to be used in a given element (see note on subroutine GAUSS2). See note on GAUSS2. The channel number to which the disc or tape is attached where the buffers are stored. See parameter list of CARTD2. The maximum number of elements allowed per buffer. See parameter list of CARTD2. See parameter list of CARTD2. See parameter list of CARTD2.

Adr. Eng. So[tware, 1981, Vol. 3, No. 2

MNODE MNPEL MSIZE

The maximum number of nodes in the mesh. See parameter list of CARTD2. The size of the work area where each buffer is stored in core. See parameter list of CARTD2. See parameter list of CARTD2. See parameter fist of CARTD2.

CARTD DERIV DETJA ELCOD See parameter list of CARTD2. GAUSW Real array of size (MGAUS) which is used within CARTD3 to contain the integration weighting factors for each sampling point. NDIMB Integer array of size (MELEB) that is used within the current buffer. NGAUB Integer array of size (MELEB) that is used within CARTD3 to contain the number of Gauss points in each element within the current buffer. NNPEB Integer array of size (MELEB) that is used within CARDT3 to contain the number of nodal points in each element within the current buffer. SHAPE Real array of dimensions (MNPEL, MGAUS) used internally within CARTD3 to contain the values of shape functions within an element. WORKA Real array of size (MSIZE) used as the work area where each buffer is stored until written onto backing store. MCART See parameter list of CARTD2. Of these parameters the arrays CARTD, DERIV, DETJA, ELCOD, GAUSW, NDIMB, NGUAB, NNPEB, SHAPE, W O R K A are used within the subroutine only. Thus they may be overwritten with other material in the calling segment and may be treated as work areas. All the remaining data must be set up in the machine prior to calling CARTD3. CRCA~311tlltllttlllllllllllt

lllti

llll

ltllllllllll

Jlllll

lltlllil

liliilJ

C

SUBROUTINE CARTD3 (COORD, NCONN, NGCOD, JCHAN, KCHAN, bg)IME, MELEB, NELEM, MELEM, MGAUS, MNDDE, MNPEL, MSIZE, CARTD, DERIV, DETJA, ELCOD, GAUSW, NDIMB, NGAUB, NNPEB, SHAPE, WORKA, MCART) C C C C C C C C C

SUBROUTINE TO CALCULATE THE CARTESIAN DERIVATIVES OF ALL SHAPEFUNCTIONS AND BUFFERTHEMTOGETHERPRIOR TO WRITING THEMON TO CHANNELKCHAN INTRODUCED AT VERSIONGO60 UPGRADED AT VERSIONG075 (AUGUST 1980) IAN WHITE - APRIL 1980 DIMENSION GAUSW(MGAUS), SHAPE(MNPEL,MGAUS), DERIV(NDIME, MNPEL,blGAUS), NCONN(MNPEL,MELE/4), COORD(PDIME,MNDDE) , DETJA(HGAUS), ELCOD(FDIME,MNPEL), WORKA(MSIZE), NGCOD(HELEM), NGAUB(MELEB), NDIHB(MELEB), NNPEB(HELEB), CARTD(MCART) REAL JACOB(3,3) LOGICAL GCALL CO~ON/IANLIB/LCHAN DATA GCALL/.TRUE./

C C C C

INITIALISE COUNTERS, TAPE, SET UP GAUSS CODE

SCAN ELEMENTS AND

REWIND KCHAN WRITE (KCHAN) I~)IME, MNPEL TOLER = I.E-08 KPOSB = I KELEB = I DO 50 IELEM = 1,NELEI4 IF (GCALL) CALL GAUSS2 (NGCOD(IELEM), JCHAN, ~IME, HNPEL, HGAUS, NDIME, NNPEL, NGAUS, GAUSW, SHAPE, DERIV) GCALL = NGCOD(HOD(IELEM,NELEM)+1). NE.NGCOD(IELEM) C C C

SET UP COORDSOF NODESIN ELEMENT 'IELE~' DO 10 INPEL = I,NNPEL NODEN = IABS(NCONN(INPEL,IELEM))

I0

DO 10 IDIME : I,NDIME ELCOD(IDIME,INPEL) : COORD(IDIME,NODEN)

CICA~rDLjJII

SET UP CARTESIAN DERIVATIVES

NCART = NDIHEmHNPELmNGAUS NEXPO = NGAUS + NCART NEEDS = KPOSB + NEXPO IF (NEEDS.LE.MSIZE) GOTO 20 IF (IELEM. EO. I) CALL ERROR ( 'CARTD3 j' ~ I t LCHANp .TRUE. ) NELEB = KELEB - 1 NPOSB = KPOSB - 1 WRITE (KCHAN) HSIZEp NPOSB, NELEBj (NGAUB(IELEB)p IELEB = 1,NELEB), (NDIHB(IELEB), IELEB = 1,NELEB)t (NNPEB(IELEB), IELEB = 1, NELEB), • (WORKA(IPOSB), IPOSB = 1,NPOSB) KPOSB = 1 KELEB = 1 C C C 2O 30 qo

5O C C C

C C C C C C C C

Accessing and dismantling the buffers T h e s u b r o u t i n e C A R T D 4 is called within the stiffness m a t r i x a s s e m b l y r o u t i n e in o r d e r to set u p the cartesian derivatives. C h e c k s are m a d e o n the m a x i m u m n u m b e r of c o o r d i n a t e d i m e n s i o n s a n d n o d e s per element when the file is initialized. This is necessary in o r d e r to ensure t h a t the d i m e n s i o n s of the C A R T D a r r a y are identical on b o t h the calls to the C A R T D 3 a n d C A R T D 4 subroutines. C A R T D 4 s h o u l d be accessed sequentially, once for each element a n d a test is included to ensure t h a t this is done. T h e r o u t i n e then reads in a new buffer w h e n e v e r a new p a t c h of elements is to be processed a n d sets up the cartesian derivatives from the w o r k area.

C C C

JJJJJ J

IF WE'RE STARTING WITH ELEMENT I SEToUP THE TAPE IF (IELEM.GT.I) GOTO 5 REWIND KCHAN READ (KCHAN) KDIHE, KNPEL IF (g2)IHE.NE.~IME) CALL ERROR ('CARTD4*'~ 3, LCHAN, .TRUE.) IF (KNPEL.NE.MNPEL) CALL ERROR ('CARTDq* ' t q, LCHANj .TRUE.)

C C C

IF A NEW BUFFER IS REOUIRED READ ONE IN 5

10 C C C 20 C C C

30

IF (KELEB.LE.NELEB.AND.IELEM.NE.I) GOTO 10 READ (KCHAN) NSIZE, NPOSB, NELEBj (NGAUB(IELEB), IELEB = I t NELEB) f (NDIHB(IELEB), IELEB=I,NELEB) (NNPEB(IELE3)j IELEB = I,NELEB)p • (WORKA(IPOSB), IPOSB= I.NPOSB) IF (NPOSB.GT.HSIZE) CALL ERROR ('CARTDq it, It LCHANm .TRUE.) IF (NELEB.GT.MELEB) CALL ERROR ('CARTDq m', 2, LCHANf .TRUE.) KPOSB = 1 KELEB = 1 NGAUS= NGAUB(KELEB) NDIME = NDIMB(KELED) NNPEL = NNPEB(KELEB) SET UP DVOLU

DO 20 IGAUS = 1,NGAUS DVOLU(IGAUS) = WORKA(KP(~B) KPOSB= KPOSB + 1 SET UP CARTD • NCART = Pg)IHE*HNPEL*NGAUS IF (NCART.GT.MCART) CALL ERROR (~CARTDq IT, 5, LCAHN, .TRUE.) DO 30 ICART = I,NCART CANTD(ICART) = WORKA(KPO~B) KPOSB = KPOSB + I KELEB = KELEB + 1 RETURN END

REFERENCES 1 2 3 4

Parameter list CARTD A real a r r a y as d e s c r i b e d in c o n n e c t i o n with C A R T D 2 which, on return, will c o n t a i n the cartesian derivatives in direction i of s h a p e function j at i n t e g r a t i o n s a m p l i n g p o i n t k in C A R T D (ij,k). DVOLU A real a r r a y of size ( M G A U S ) in which the values of [o~glJl-lg will be stores g = 1. . . . . NGAUS. IELEM T h e n u m b e r of the element c u r r e n t l y being processed. All o t h e r p a r a m e t e r s have a l r e a d y been d e s c r i b e d in c o n n e c t i o n with s u b r o u t i n e s C A R T D 2 a n d C A R T D 3 .

J Jill

IAN WHITE IKI'RODUCED AT IANLIB VERSION G064 (JULY 1980) UPGRADED AT VERSION GO?5 (AUOUST 1980) DIMENSION ROAUB(HELEB)~ NDIHB(HELEB), NNPEB(HELEB), WOREA(HSIZE) ~ DVOLU(HGAUS), • CANTD(MCANT) COWHON/IANLIB/LCHAN DATA KELEB/O/~ NELEB/O/

WRITE REMAINDER OF BUFFER NELEB = KELE8 - I NPOSB = KPOSB - 1 WRITE (KCHAN) NSIZE, NPOSB, NELEB~ (NGAUB(IELEM), IELEM = 1,NELEB), (NDIHB(IELEB), IELER = 1, NELEB), (NNPEB(IELEB), IELEB = 1, NELEB), (WORKA(IPOSB), ~POSB = 1,~POSB) "RETURN END

JlJJJJJJl

SUBROUTINE TO READ THE CARTESIAN DERIVATIVES FROM CHANNEL 'KCHAN' AND DISMANTLE THE BUFFER APPROPRIATELY

PUT CURRENT ELEMENT INTO THE BUFFER DO 30 IGAUS = 1, NGAUS WORKA(KPOSB) = DETJA(IGAUS)*GAUSW(IGAUS) KPOSB = KPG~B + 1 DO 40 ICART = 1,NCANT WORKA(KPOSB) = CARTD(ICART) KPOSB = KPOSB + t NGAUB(KELEB) = NGAUS NDIMB(KELEB) = NDIHE NNPEB(KELEB) = NNPEL KELEB = KELEB + 1

Jl JJ J J JJJJ J| JJJJ JJJJlJJJJJJllJJi

SUBROUTINE CARTD~ (CANTD~ DVOLU~ NNPEL, NGAUS, NDIME, IELEM, KCHAN, MNPELj HGAUS, FU)IMEt MELEB MCANT~ HSIZE, NGAUBp NDIMB, NNPEB, WORKAI

CALL CARTD2 (CARTD, DETJA, DERIV~ ELCOD~ TOLER, NDIHE, HDIME, NGAUSr HGAUS~ NNPELp HNPELt HCART) EVALUATE HOW MUCH MORE SPACE THIS ELEMENT WILL OCCUPY II THE BUFFER AND WRITE OUT THE BUFFER IF THERE ISN'T ENOUGH ROOM

JJJlJllJJ

C

5 6 7 8 9 10

Hinton, E. and Owen, D. R. J. Finite Element Programming, Academic Press, New York and London, 1977 Hinton, E. and Owen, D. R. J. Finite Elements in Plasticity, Pineridge Press, Swansea, 1980 Zienkiewicz, O. C. The Finite Element Method, McGraw-Hill, New York, 1977 Brebbia, C. A. and Connor, J. J. Fundamentals of Finite Element Techniques, Butterworths, London, 1973 Connor, J. J. and Brebbia, C. A. Finite Element Techniquesfor Fluid Flow, Butterworths, London, 1976 Irons, B. M. A frontal solution program, Int. J. Num. Math. Eng. 1970, 2, 5 ICL Extended Fortran, ICL Software Distribution, Reading, Technical Pub. 4269 Bathe, K.-J. and Wilson, E. L. Numerical Methods in Finite Element Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1976 Gupta, A. K. and Mohraz, B. A method for computing numerically integrated stiffness matrices, Int. J. Num. Meth. Eng. 1972,5, 83 Damrath, R. and Pahl, P. J. Data structures for a system of finite element programs, Proc. US-German Syrup. Formulations and Algorithms in Fhdte Element Analysis, MIT Press, 1977, Chapter 6

Adv. Eng. Softlvare, 1981, Vol. 3, No. 2

83