31
32
33
34
35
36
37
38#include "implicit_f.inc"
39
40
41
42 INTEGER, INTENT(IN) :: NMN, NSN, NRT, NUMNOD
43 INTEGER, INTENT(IN) :: IRECT(4,*), NSEG(*), NSV(*), MSR(*)
44 INTEGER, INTENT(OUT) :: LCSEG(*), ILOC(*)
46
47
48
49 INTEGER II, I, JJ, J, K, L, N, SEG1
51 . cms, dms,xmax_m,ymax_m,zmax_m,xmax_s,ymax_s,zmax_s,
52 . xmin_m,ymin_m,zmin_m,xmin_s,ymin_s,zmin_s,
53 . xmin,ymin,zmin,xmax,
ymax,zmax,xs,ys,zs,dist,mindist,
54 . x1,y1,z1,aaa
55 INTEGER IIX(NSN),IIY(NSN),IIZ(NSN)
56 INTEGER LAST_NOD(NSN),NEXT_NOD(NSN)
57 INTEGER NBX,NBY,NBZ,NE,FOUND
58 INTEGER FIRST,LAST,IX1,IY1,IZ1,NN,IX,IY,IZ,IX2,IY2,IZ2
59 INTEGER, DIMENSION(:,:,:),ALLOCATABLE :: VOXEL
60 INTEGER, DIMENSION(:), ALLOCATABLE :: KNOD2SEG
61 INTEGER, DIMENSION(:), ALLOCATABLE :: NOD2SEG
62
63 IF(nrt==0) RETURN
64
65 next_nod = 0
66 iix = 0
67 iiy = 0
68 iiz = 0
69
70 xmax_m=-ep30
71 ymax_m=-ep30
72 zmax_m=-ep30
73 xmin_m= ep30
74 ymin_m= ep30
75 zmin_m= ep30
76
77 DO i=1,nmn
78 j=msr(i)
79 xmax_m=
max(xmax_m,x(1,j))
80 ymax_m=
max(ymax_m,x(2,j))
81 zmax_m=
max(zmax_m,x(3,j))
82 xmin_m=
min(xmin_m,x(1,j))
83 ymin_m=
min(ymin_m,x(2,j))
84 zmin_m=
min(zmin_m,x(3,j))
85 ENDDO
86
87 xmin_s= ep30
88 xmax_s=-ep30
89 ymin_s= ep30
90 ymax_s=-ep30
91 zmin_s= ep30
92 zmax_s=-ep30
93
94 DO i=1,nsn
95 j=nsv(i)
96 xmin_s=
min(xmin_s,x(1,j))
97 ymin_s=
min(ymin_s,x(2,j))
98 zmin_s=
min(zmin_s,x(3,j))
99 xmax_s=
max(xmax_s,x(1,j))
100 ymax_s=
max(ymax_s,x(2,j))
101 zmax_s=
max(zmax_s,x(3,j))
102 ENDDO
103 xmin=
min(xmin_m,xmin_s)
104 ymin=
min(ymin_m,ymin_s)
105 zmin=
min(zmin_m,zmin_s)
106 xmax=
max(xmax_m,xmax_s)
108 zmax=
max(zmax_m,zmax_s)
109
110
111 aaa = sqrt(nsn /
112 . ((xmax-xmin)*(
ymax-ymin)
113 . +(
ymax-ymin)*(zmax-zmin)
114 . +(zmax-zmin)*(xmax-xmin)))
115
116 aaa = three_over_4*aaa
117
118 nbx = nint(aaa*(xmax-xmin))
119 nby = nint(aaa*(
ymax-ymin))
120 nbz = nint(aaa*(zmax-zmin))
124
125
126
127 ALLOCATE(voxel(nbx,nby,nbz))
128 voxel=0
129 DO i=1,nsn
130 iix(i)=0
131 iiy(i)=0
132 iiz(i)=0
133 j=nsv(i)
134
135 iix(i)=int(nbx*(x(1,j)-xmin)/(xmax-xmin))
136 iiy(i)=int(nby*(x(2,j)-ymin)/(
ymax-ymin))
137 iiz(i)=int(nbz*(x(3,j)-zmin)/(zmax-zmin))
138
139 iix(i)=
max(1,
min(nbx,iix(i)))
140 iiy(i)=
max(1,
min(nby,iiy(i)))
141 iiz(i)=
max(1,
min(nbz,iiz(i)))
142
143 first = voxel(iix(i),iiy(i),iiz(i))
144 IF(first == 0)THEN
145
146 voxel(iix(i),iiy(i),iiz(i)) = i
147 next_nod(i) = 0
148 last_nod(i) = 0
149 ELSEIF(last_nod(first) == 0)THEN
150
151
152 next_nod(first) = i
153 last_nod(first) = i
154 next_nod(i) = 0
155 ELSE
156
157
158 last = last_nod(first)
159 next_nod(last) = i
160 last_nod(first) = i
161 next_nod(i) = 0
162 ENDIF
163 ENDDO
164
165
166
167 DO ne=1,nmn
168 mindist = ep30
169 x1 = x(1,msr(ne))
170 y1 = x(2,msr(ne))
171 z1 = x(3,msr(ne))
172 ix1=int(nbx*(x1-xmin)/(xmax-xmin))
173 iy1=int(nby*(y1-ymin)/(
ymax-ymin))
174 iz1=int(nbz*(z1-zmin)/(zmax-zmin))
175
179
180 ix2=ix1
181 iy2=iy1
182 iz2=iz1
183
184 found = 0
185 DO WHILE(found == 0)
186 DO iz = iz1,iz2
187 DO iy = iy1,iy2
188 DO ix = ix1,ix2
189
190 jj = voxel(ix,iy,iz)
191
192 DO WHILE(jj /= 0)
193 nn=nsv(jj)
194 xs = x(1,nn)
195 ys = x(2,nn)
196 zs = x(3,nn)
197 dist = (x1-xs)**2+(y1-ys)**2+(z1-zs)**2
198 IF( dist < mindist ) THEN
199 iloc(ne)=jj
200 mindist = dist
201 ENDIF
202 found = 1
203 jj = next_nod(jj)
204 ENDDO
205 ENDDO
206 ENDDO
207 ENDDO
208 ix1 = ix1-1
209 iy1 = iy1-1
210 iz1 = iz1-1
214 ix2 = ix2+1
215 iy2 = iy2+1
216 iz2 = iz2+1
220 ENDDO
221
222
223 x1 = x(1,msr(ne))
224 y1 = x(2,msr(ne))
225 z1 = x(3,msr(ne))
226 ix1=int(nbx*(x1-mindist-xmin)/(xmax-xmin))
227 iy1=int(nby*(y1-mindist-ymin)/(
ymax-ymin))
228 iz1=int(nbz*(z1-mindist-zmin)/(zmax-zmin))
229
233
234 ix2=int(nbx*(x1+mindist-xmin)/(xmax-xmin))
235 iy2=int(nby*(y1+mindist-ymin)/(
ymax-ymin))
236 iz2=int(nbz*(z1+mindist-zmin)/(zmax-zmin))
237
241
242 DO iz = iz1,iz2
243 DO iy = iy1,iy2
244 DO ix = ix1,ix2
245
246 jj = voxel(ix,iy,iz)
247
248 DO WHILE(jj /= 0)
249 nn=nsv(jj)
250 xs = x(1,nn)
251 ys = x(2,nn)
252 zs = x(3,nn)
253 dist = (x1-xs)**2+(y1-ys)**2+(z1-zs)**2
254 IF( dist < mindist ) THEN
255 iloc(ne)=jj
256 mindist = dist
257 ENDIF
258 jj = next_nod(jj)
259 ENDDO
260 ENDDO
261 ENDDO
262 ENDDO
263 ENDDO
264 DEALLOCATE(voxel)
265
266
267
268 ALLOCATE(knod2seg(numnod+1))
269 ALLOCATE(nod2seg(4*nrt))
270 nod2seg(1:4*nrt) = 0
271 knod2seg(1:numnod+1) = 0
272 DO i=1,nrt
273 DO k=1,4
274 n = irect(k,i)
275 knod2seg(n) = knod2seg(n) + 1
276 END DO
277 END DO
278
279 DO i=1,numnod
280 knod2seg(i+1) = knod2seg(i+1) + knod2seg(i)
281 END DO
282
283 DO n=numnod,1,-1
284 knod2seg(n+1)=knod2seg(n)
285 END DO
286 knod2seg(1)=0
287
288 DO i=1,nrt
289 DO k=1,4
290 n = irect(k,i)
291 knod2seg(n) = knod2seg(n) + 1
292 nod2seg(knod2seg(n)) = i
293 END DO
294 END DO
295
296 DO n=numnod,1,-1
297 knod2seg(n+1)=knod2seg(n)
298 END DO
299 knod2seg(1)=0
300
301
302
303 ii=0
304 DO i=1,nsn
305 k = nsv(i)
306 DO j=knod2seg(k)+1,knod2seg(k+1)
307 seg1 = nod2seg(j)
308 DO l=1,4
309 IF (irect(l,seg1) /= k)THEN
310 ii=ii+1
311 lcseg(ii)=seg1
312 EXIT
313 ENDIF
314 ENDDO
315 ENDDO
316 ENDDO
317
318 DEALLOCATE(knod2seg)
319 DEALLOCATE(nod2seg)
320 RETURN
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)