1!$Id:$ 2 subroutine int1dg(l,sw) 3 4! * * F E A P * * A Finite Element Analysis Program 5 6!.... Copyright (c) 1984-2017: Regents of the University of California 7! All rights reserved 8 9!-----[--.----+----.----+----.-----------------------------------------] 10! Purpose: Gauss quadrature for 1-d element 11 12! Inputs: 13! l - Number of points 14 15! Outputs: 16! sw(1,*) - Gauss points 17! sw(2,*) - Gauss weights 18!-----[--.----+----.----+----.-----------------------------------------] 19 implicit none 20 21 integer l 22 real*8 sw(2,*), t 23 24 save 25 26 if(l.eq.1) then 27 28 sw(1,1) = 0.0d0 29 sw(2,1) = 2.0d0 30 31 elseif(l.eq.2) then 32 33 sw(1,1) = -1.d0/sqrt(3.d0) 34 sw(1,2) = -sw(1,1) 35 sw(2,1) = 1.0d0 36 sw(2,2) = 1.0d0 37 38 elseif(l.eq.3) then 39 40 sw(1,1) = -sqrt(0.6d0) 41 sw(1,2) = 0.0d0 42 sw(1,3) = -sw(1,1) 43 sw(2,1) = 5.d0/9.d0 44 sw(2,2) = 8.d0/9.d0 45 sw(2,3) = sw(2,1) 46 47 elseif(l.eq.4) then 48 49 t = sqrt(4.8d0) 50 sw(1,1) = -sqrt((3.d0+t)/7.d0) 51 sw(1,2) = -sqrt((3.d0-t)/7.d0) 52 sw(1,3) = -sw(1,2) 53 sw(1,4) = -sw(1,1) 54 t = 1.d0/3.d0/t 55 sw(2,1) = 0.5d0 - t 56 sw(2,2) = 0.5d0 + t 57 sw(2,3) = sw(2,2) 58 sw(2,4) = sw(2,1) 59 60 elseif(l.eq.5) then 61 62 t = sqrt(1120.0d0) 63 64 sw(1,1) = (70.d0+t)/126.d0 65 sw(1,2) = (70.d0-t)/126.d0 66 67 t = 1.d0/(15.d0 * (sw(1,2) - sw(1,1))) 68 69 sw(2,1) = (5.0d0*sw(1,2) - 3.0d0)*t/sw(1,1) 70 sw(2,2) = (3.0d0 - 5.0d0*sw(1,1))*t/sw(1,2) 71 sw(2,3) = 2.0d0*(1.d0 - sw(2,1) - sw(2,2)) 72 sw(2,4) = sw(2,2) 73 sw(2,5) = sw(2,1) 74 75 sw(1,1) = -sqrt(sw(1,1)) 76 sw(1,2) = -sqrt(sw(1,2)) 77 sw(1,3) = 0.0d0 78 sw(1,4) = -sw(1,2) 79 sw(1,5) = -sw(1,1) 80 81! Compute points and weights 82 83 else 84 85 call gausspw(l,sw) 86 87 endif 88 89 end 90 91 subroutine gausspw (nn,sw) 92 93!-----[--.----+-!--.----+----.----+------------------------------------] 94! Input: 95! nn = Number of Gauss Points to compute 96 97! Outputs: 98! sw(1,*) = Gauss Point Coordinates 99! sw(2,*) = Gauss Point Weights 100!-----[--.----+-!--.----+----.----+------------------------------------] 101 implicit none 102 103 integer nn, n 104 real*8 sw(2,*) 105 real*8 fn, beta, cc, xt, dpn,pn1, flgama 106 107 fn = dble(nn) 108 beta = exp(2.d0*flgama(1.d0) - flgama(2.d0)) 109 cc = 2.d0*beta 110 do n = 2,nn 111 cc = cc*4.d0*dble(n-1)**4 112 & / (dble(2*n-1)*dble(2*n-3)*dble(2*n-2)**2) 113 end do ! n 114 115 do n = 1,nn 116 117! Largest zero 118 119 if(n.eq.1) then 120 xt = 1.d0 - 2.78d0/(4.d0 + fn*fn) 121 122! Second zero 123 124 elseif(n.eq.2) then 125 xt = xt - (4.1d0 + 0.246d0*(fn - 8.d0)/fn)*(1.d0 - xt) 126 127! Third zero 128 129 elseif(n.eq.3) then 130 xt = xt - (1.67d0 + 0.3674d0*(fn - 8.d0)/fn)*(sw(1,1) - xt) 131 132! Second last zero 133 134 elseif(n.eq.nn-1) then 135 xt = xt + (xt - sw(1,n-2))/0.766d0/(1.d0 + 0.639d0*(fn-4.d0) 136 & /(1.d0 + 0.710d0*(fn-4.d0))) 137 138! Last zero 139 140 elseif(n.eq.nn) then 141 xt = xt + (xt - sw(1,n-2))/1.67d0/(1.d0 + 0.22d0*(fn-8.d0)/fn) 142 143! Intermediate roots 144 145 else 146 xt = 3.d0*sw(1,n-1) - 3.d0*sw(1,n-2) + sw(1,n-3) 147 endif 148 149! Find root using xt-value 150 151 call root (xt,nn,dpn,pn1) 152 sw(1,n) = xt 153 sw(2,n) = cc/(dpn*pn1) 154 155 end do ! n 156 157! Reverse order of points 158 159 do n = 1, nn 160 sw(1,n) = - sw(1,n) 161 end do ! n 162 163 end 164 165 subroutine root (x,nn,dpn,pn1) 166 167!-----[--.----+-!--.----+----.----+------------------------------------] 168! Improve approximate root x; in addition we also obtain 169! dpn = derivative of p(n) at x 170! pn1 = value of p(n-1) at x 171!-----[--.----+-!--.----+----.----+------------------------------------] 172 implicit none 173 174 logical notconv 175 integer nn, iter 176 real*8 x,dpn,pn1, d,p,dp 177 178 real*8 eps 179 data eps / 1.d-39 / 180 181 iter = 0 182 183 notconv = .true. 184 do while(notconv .and. iter.lt.50) 185 iter = iter + 1 186 call recur (p,dp,pn1,x,nn) 187 d = p/dp 188 x = x - d 189 if(abs(d).le.eps) then 190 notconv = .false. 191 endif 192 end do ! while 193 dpn = dp 194 195 end 196 197 subroutine recur (pn,dpn,pn1,x,nn) 198 199 implicit none 200 201 integer nn, n 202 real*8 pn,dpn,pn1,x, c 203 real*8 p1, p,dp, dp1, q,dq 204 205 p1 = 1.d0 206 p = x 207 dp1 = 0.d0 208 dp = 1.d0 209 do n = 2,nn 210 c = 4.d0*dble(n-1)**4 211 & / (dble(2*n-1)*dble(2*n-3)*dble(2*n-2)**2) 212 q = x*p - c*p1 213 dq = x*dp - c*dp1 + p 214 p1 = p 215 p = q 216 dp1 = dp 217 dp = dq 218 end do ! n 219 pn = p 220 dpn = dp 221 pn1 = p1 222 223 end 224 225 function flgama (w) 226 227 implicit none 228 229 integer m, i 230 real*8 flgama, w, pi,x, p, fk, y, z,zz 231 232 pi = acos(-1.d0) 233 x = w 234 fk = -1.d0 235 236! w less eq 0.5 237 238 if (x .lt. 0.5d0) then 239 m = 1 240 p = pi/sin(x*pi) 241 x = 1.d0 - x 242 else 243 m = 0 244 p = 0.0d0 245 endif 246 247 do while(x + fk - 6.d0 .le.0.0d0) 248 fk = fk + 1.d0 249 end do ! while 250 251 z = x + fk 252 zz = z*z 253 254 y = (z - 0.5d0)*log(z) - z + 0.9189385332047d0 + (((((-4146.d0/zz 255 & + 1820.d0)/zz - 1287.d0)/zz + 1716.d0)/zz - 6006d0)/zz 256 & + 180180.d0)/z/2162160.d0 257 258 if(fk.gt.0.0d0) then 259 do i = 1,int(fk) 260 fk = fk - 1.d0 261 y = y - log(x + fk) 262 end do ! i 263 endif 264 265 if(m.ne.0) then 266 if(p.le.0.0d0) then 267 write (*,2000) w 268 y = 0.d0 269 else 270 y = log(p) - y 271 endif 272 endif 273 flgama = y 274 2752000 format (2x,'gamma(',e11.4,') is negative') 276 277 end 278