Mathematical Physics - Volume II - Numerical Methods

3.6 Program for solving of elliptical problems

91

COMMON /cmatrx/ gk(maxn,maxband) , gf(maxn) COMMON /cint1 / etaq(3) , ww(3) COMMON /cint2 / xiq(2,9) , xit(2,7) , wq(9) , wt(7)

c c

OPEN (UNIT=5,FILE=’fem.in’) OPEN (UNIT=6,FILE=’fem.out’) CALL setint

100 CALL prep CALL proc CALL post GOTO 100 END

SUBROUTINE setint COMMON /cint1 / etaq(3) , ww(3) COMMON /cint2 / xiq(2,9) , xit(2,7) , wq(9) , wt(7)

c nine point quadrature rule for quadrilaterals sq35 = sqrt(3./5.)

wq(1) = 25./81. wq(2) = 40./81. wq(3) = 25./81. wq(4) = 40./81. wq(5) = 64./81. wq(6) = 40./81. wq(7) = 25./81. wq(8) = 40./81. wq(9) = 25./81. xiq(1,1) = -sq35 xiq(2,1) = -sq35 xiq(1,2) = 0. xiq(2,2) = -sq35 xiq(1,3) = sq35 xiq(2,3) = -sq35 xiq(1,4) = -sq35 xiq(2,4) = 0. xiq(1,5) = 0. xiq(2,5) = 0. xiq(1,6) = sq35 xiq(2,6) = 0. xiq(1,7) = -sq35 xiq(2,7) = sq35 xiq(1,8) = 0.

xiq(2,8) = sq35 xiq(1,9) = sq35 xiq(2,9) = sq35 c seven point rule for triangulars sq15 = sqrt(15.0) a = (155.0+sq15)/2400.0