Mathematical Physics - Volume II - Numerical Methods

Chapter 3. Comparison of finite element method and finite difference method

98

DIMENSION psi(9) , dpsi(9,2)

c initilaze element arrays DO 100 i = 1 , n ef(i) = 0. DO 50 j = 1 , n ek(i,j) = 0. 50 CONTINUE 100 CONTINUE CALL getmat(xk,xb,xf,mat) c begin integration point loop DO 300 l = 1 , nl

CALL shape2(xi(1,l),xi(2,l),n,psi,dpsi)

c

calculate dxds DO 150 i = 1 , 2

DO 120 j = 1 , 2 dxds(i,j) = 0.

DO 110 k = 1 , n dxds(i,j) = dxds(i,j) + dpsi(k,j)*x(i,k)

110 120

CONTINUE

CONTINUE

150 CONTINUE c

calculate dsdx detj = dxds(1,1)*dxds(2,2) - dxds(1,2)*dxds(2,1)

IF ( detj.LE.0. ) GOTO 500 dsdx(1,1) = dxds(2,2)/detj dsdx(2,2) = dxds(1,1)/detj dsdx(1,2) = -dxds(1,2)/detj dsdx(2,1) = -dxds(2,1)/detj

c

calculate d(psi)/dx DO 200 i = 1 , n

dpsix(i) = dpsi(i,1)*dsdx(1,1) + dpsi(i,2)*dsdx(2,1) dpsiy(i) = dpsi(i,1)*dsdx(1,2) + dpsi(i,2)*dsdx(2,2)

200 CONTINUE c

accumulate integration point value of integrals fac = detj*w(l) DO 250 i = 1 , n ef(i) = ef(i) + xf*psi(i)*fac DO 220 j = i , n ek(i,j) = ek(i,j) &

+ fac*(xk*(dpsix(i)*dpsix(j)+dpsiy(i)*dpsiy(j)

&

)+xb*psi(i)*psi(j))

220

CONTINUE

250 CONTINUE 300 CONTINUE c calculate lower symmetric part of ek DO 400 i = 1 , n DO 350 j = 1 , i ek(i,j) = ek(j,i) 350 CONTINUE