1!$Id:$ 2 subroutine vblkn(nr,ns,nt,xl,x,ixl,dr,ds,dt, 3 & ni,ndm,ctype,prt,prth) 4 5! * * F E A P * * A Finite Element Analysis Program 6 7!.... Copyright (c) 1984-2017: Regents of the University of California 8! All rights reserved 9 10!-----[--.----+----.----+----.-----------------------------------------] 11! Purpose: Generate a block of 3-d 8-node brick elements 12 13! Inputs: 14! nr - Number elements in 1-local coordinate dir. 15! ns - Number elements in 2-local coordinate dir. 16! nt - Number elements in 3-local coordinate dir. 17! xl(ndm,*) - Block nodal coordinate array 18! ixl(*) - Block nodal connection list 19! dr - 1-local coordinate increment 20! ds - 2-local coordinate increment 21! dt - 3-local coordinate increment 22! ni - Initial node number for block 23! ndm - Spatial dimension of mesh 24! ctype - Type of block coordinates 25! prt - Output generated data if true 26! prth - Output title/header data if true 27 28! Outputs: 29! x(ndm,*) - Nodal coordinates for block 30!-----[--.----+----.----+----.-----------------------------------------] 31 32 implicit none 33 34 include 'cdata.h' 35 include 'cdat2.h' 36 include 'iofile.h' 37 include 'trdata.h' 38 39 logical prt,prth,phd, pcomp 40 character xh*6, ctype*15 41 integer ni,ndm,nr,ns,nt,i,j,k,l,m,n,mct, ixl(*) 42 real*8 dr,ds,dt, rr, sn2,cn2,sn3,cn3 43 real*8 ss(3),xl(3,*),x(ndm,*),xx(3) 44 45 save 46 47 data xh/' coord'/ 48 49! Check that all corners of brick are defined 50 51 do k = 1,3 52 xx(k) = 0.0d0 53 end do ! k 54 55 do k = 1,8 56 if(ixl(k).ne.k) go to 900 57 end do 58 call bcor3d(ixl,xl) 59 n = ni 60 mct = 0 61 ss(3) = -1.0d0 62 do k = 1,nt 63 ss(2) = -1.0d0 64 do j = 1,ns 65 ss(1) = -1.0d0 66 do i = 1,nr 67 68! Compute coordinates of node 69 70 call xbcor3d(ss,xl, xx) 71 72! Convert coordinates if necessary 73 74 if(pcomp(ctype,'pola',4)) then 75 call pdegree(xx(2), sn2,cn2) 76 rr = xx(1) 77 xx(1) = x0(1) + rr*cn2 78 xx(2) = x0(2) + rr*sn2 79 xx(3) = x0(3) + xx(3) 80 elseif(pcomp(ctype,'sphe',4)) then 81 call pdegree(xx(2), sn2,cn2) 82 call pdegree(xx(3), sn3,cn3) 83 rr = xx(1) 84 xx(1) = x0(1) + rr*cn2*sn3 85 xx(2) = x0(2) + rr*sn2*sn3 86 xx(3) = x0(3) + rr*cn3 87 endif 88 89! Transform to global coordinates 90 91 do m = 1,ndm 92 x(m,n) = xr(m)+tr(m,1)*xx(1)+tr(m,2)*xx(2)+tr(m,3)*xx(3) 93 end do 94 95! Output point 96 97 if(prt) then 98 mct = mct + 1 99 phd = mod(mct,50).eq.1 100 call prtitl(prth.and.phd) 101 if(phd) write(iow,2000) (l,xh,l=1,ndm) 102 write(iow,2001) n,(x(l,n),l=1,ndm) 103 if(ior.lt.0) then 104 if(phd) write(*,2000) (l,xh,l=1,ndm) 105 write(*,2001) n,(x(l,n),l=1,ndm) 106 endif 107 endif 108 n = n + 1 109 ss(1) = ss(1) + dr 110 end do 111 ss(2) = ss(2) + ds 112 end do 113 ss(3) = ss(3) + dt 114 end do 115 116 return 117 118! Error 119 120900 write(iow,3000) k 121 if(ior.lt.0) then 122 write(*,3000) k 123 return 124 endif 125 call plstop(.true.) 126 127! Formats 128 1292000 format(/' N o d a l C o o r d i n a t e s'//6x,'Node',5(i7,a6)) 130 1312001 format(i10,5f13.4) 132 1333000 format(' *ERROR* Block node',i3,' is undefined') 134 135 end 136 137 subroutine bcor3d(ixl,xl) 138 139 implicit none 140 141 integer ixl(27), imid(12),amid(12),bmid(12) 142 real*8 xl(3,27) 143 144 integer i,j 145 146 data imid/9,10,11,12, 13,14,15,16, 18,19,20,21/ 147 data amid/1, 2, 3, 4, 1, 2, 3, 4, 5, 6, 7, 8/ 148 data bmid/5, 6, 7, 8, 2, 3, 4, 1, 6, 7, 8, 5/ 149 150! Mid edge coordinates 151 152 do i = 1,12 153 if(ixl(imid(i)).eq.0) then 154 do j = 1,3 155 xl(j,imid(i)) = 0.5d0*(xl(j,amid(i)) + xl(j,bmid(i))) 156 end do ! j 157 ixl(i) = i 158 endif 159 end do ! i 160 161! Bottom and top 162 163 if(ixl(17).eq.0) then 164 do j = 1,3 165 xl(j,17) = 0.50d0*(xl(j,13) + xl(j,14) + xl(j,15) + xl(j,16)) 166 & - 0.25d0*(xl(j, 1) + xl(j, 2) + xl(j, 3) + xl(j, 4)) 167 end do ! j 168 ixl(17) = 17 169 endif 170 171 if(ixl(22).eq.0) then 172 do j = 1,3 173 xl(j,22) = 0.50d0*(xl(j,18) + xl(j,19) + xl(j,20) + xl(j,21)) 174 & - 0.25d0*(xl(j, 5) + xl(j, 6) + xl(j, 7) + xl(j, 8)) 175 end do ! j 176 ixl(22) = 22 177 endif 178 179! Mid-face 180 181 if(ixl(23).eq.0) then 182 do j = 1,3 183 xl(j,23) = 0.50d0*(xl(j,13) + xl(j, 9) + xl(j,10) + xl(j,18)) 184 & - 0.25d0*(xl(j, 1) + xl(j, 2) + xl(j, 5) + xl(j, 6)) 185 end do ! j 186 ixl(23) = 23 187 endif 188 189 if(ixl(24).eq.0) then 190 do j = 1,3 191 xl(j,24) = 0.50d0*(xl(j,14) + xl(j,10) + xl(j,11) + xl(j,19)) 192 & - 0.25d0*(xl(j, 2) + xl(j, 3) + xl(j, 6) + xl(j, 7)) 193 end do ! j 194 ixl(24) = 24 195 endif 196 197 if(ixl(25).eq.0) then 198 do j = 1,3 199 xl(j,25) = 0.50d0*(xl(j,15) + xl(j,11) + xl(j,12) + xl(j,20)) 200 & - 0.25d0*(xl(j, 3) + xl(j, 4) + xl(j, 7) + xl(j, 8)) 201 end do ! j 202 ixl(25) = 25 203 endif 204 205 if(ixl(26).eq.0) then 206 do j = 1,3 207 xl(j,26) = 0.50d0*(xl(j,16) + xl(j,12) + xl(j, 9) + xl(j,21)) 208 & - 0.25d0*(xl(j, 1) + xl(j, 4) + xl(j, 8) + xl(j, 5)) 209 end do ! j 210 ixl(26) = 26 211 endif 212 213! Center node 214 215 if(ixl(27).eq.0) then 216 do j = 1,3 217 xl(j,27) = 0.25d0*(xl(j,13) + xl(j,14) + xl(j,15) + xl(j,16) 218 & + xl(j,18) + xl(j,19) + xl(j,20) + xl(j,21) 219 & + xl(j,23) + xl(j,24) + xl(j,25) + xl(j,26) 220 & - xl(j, 1) - xl(j, 2) - xl(j, 3) - xl(j, 4) 221 & - xl(j, 5) - xl(j, 6) - xl(j, 7) - xl(j, 8)) 222 223 end do ! j 224 ixl(27) = 27 225 endif 226 227 end 228 229 subroutine xbcor3d(ss,xl, x) 230 231!-----[--.----+----.----+----.-----------------------------------------] 232! Purpose: Compute shape functions and coordinates for each point 233 234! Inputs: 235! ss(3) - Natural coordinates for point 236! xl(3,*) - Nodal coordinates for brick 237 238! Outputs: 239! x(3) - Cartesian coordinates for point 240!-----[--.----+----.----+----.-----------------------------------------] 241 242 implicit none 243 244 integer j,l, ix(27),iy(27),iz(27) 245 real*8 ss(3),xl(3,27),x(3) 246 real*8 lshp(3,3),shp 247 248 data ix/1,3,3,1, 1,3,3,1, 1,3,3,1, 2,3,2,1,2, 2,3,2,1,2, 249 & 2,3,2,1,2/ 250 data iy/1,1,3,3, 1,1,3,3, 1,1,3,3, 1,2,3,2,2, 1,2,3,2,2, 251 & 1,2,3,2,2/ 252 data iz/1,1,1,1, 3,3,3,3, 2,2,2,2, 1,1,1,1,1, 3,3,3,3,3, 253 & 2,2,2,2,2/ 254 255 256 save 257 258 do j = 1,3 259 lshp(1,j) = 0.5d0*ss(j)*(ss(j) - 1.d0) 260 lshp(2,j) = (1.d0 - ss(j)*ss(j)) 261 lshp(3,j) = 0.5d0*ss(j)*(ss(j) + 1.d0) 262 end do ! j 263 264 do j = 1,3 265 x(j) = 0.0d0 266 end do 267 268 do l = 1,27 269 shp = lshp(ix(l),1)*lshp(iy(l),2)*lshp(iz(l),3) 270 do j = 1,3 271 x(j) = x(j) + shp*xl(j,l) 272 end do 273 end do 274 275 end 276