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