Programming Finite Element Computations - Programming the Finite Element Method - page 97

the velocity (vector) quantity and pressure (scalar) quantity and this is what is done in
programs in Chapters 9 and 12.
The velocity shape functions are designated as fun and the pressure shape functions as
funf . Similarly, the velocity derivatives are deriv and the pressure derivatives derivf .
The arrays nd1 , nd2 , ndf1 , ndf2 , nfd1 ,and nfd2 hold the results of cross products
between the velocity and pressure shape functions and their derivatives as shown below.
Thus, the element integrals which have to be evaluated numerically from equation (2.115)
are of the form

dtkd=MATMUL(MATMUL(TRANSPOSE(deriv),kay),deriv)

(3.103)

CALL cross_product(fun,deriv(1,:),nd1)

CALL cross_product(fun,deriv(2,:),nd2)

followed by

nip

c11 (=c33) = det i *weights(i)*dtkd i

i = 1

nip

+ ubar i (3.104)

det i *weights(i)*nd1 i

i = 1

nip

+ vbar i

det i *weights(i)*nd2 i

i = 1

In these equations, deriv(1,:) signifies the first row of deriv andsooninthe
usual Fortran 95 style. The diagonal terms in kay represent the reciprocal of the Reynolds
number. Note the identity of the first term of c11 with (3.63) for uncoupled flow.
The remaining submatrices are formed as follows:

CALL cross_product(fun,derivf(1,:),ndf1)

CALL cross_product(fun,derivf(2,:),ndf2)

(3.105)

CALL cross_product(funf,deriv(1,:),nfd1)

CALL cross_product(funf,deriv(2,:),nfd2)

followed by

nip

1
rho

c12 = det i *weights(i)*ndf1 i

i = 1

nip

1
rho

c32 =

det i *weights(i)*ndf2 i

i = 1