Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RFLU_ModOctree.F90
Go to the documentation of this file.
1 ! *********************************************************************
2 ! * Rocstar Simulation Suite *
3 ! * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4 ! * *
5 ! * Illinois Rocstar LLC *
6 ! * Champaign, IL *
7 ! * www.illinoisrocstar.com *
8 ! * sales@illinoisrocstar.com *
9 ! * *
10 ! * License: See LICENSE file in top level of distribution package or *
11 ! * http://opensource.org/licenses/NCSA *
12 ! *********************************************************************
13 ! *********************************************************************
14 ! * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15 ! * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16 ! * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17 ! * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18 ! * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19 ! * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20 ! * Arising FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21 ! * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22 ! *********************************************************************
23 !*******************************************************************************
24 !
25 ! Purpose: Suite of routines to carry out octree operations.
26 !
27 ! Description: None.
28 !
29 ! Notes:
30 ! 1. The routines contained in this module are from Tim Baker. They are
31 ! modified so that they no longer use COMMON statements and follow the
32 ! design philosophy of Rocflu. Some of the routines have been renamed:
33 !
34 ! The routine RFLU_BuildOctree was called OCT_SETUP by Tim Baker, and
35 ! RFLU_QueryOctree was called OCT_SEARCH.
36 !
37 ! To create and use a octree, one has to take the following steps:
38 ! 1. RFLU_CreateOctree
39 ! 2. RFLU_BuildOctree
40 ! 3. RFLU_QueryOctree
41 ! 4. RFLU_DestroyOctree
42 !
43 !*******************************************************************************
44 !
45 ! $Id: RFLU_ModOctree.F90,v 1.9 2008/12/06 08:44:23 mtcampbe Exp $
46 !
47 ! Copyright: (c) 2002-2004 by the University of Illinois
48 !
49 !*******************************************************************************
50 
51 MODULE rflu_modoctree
52 
53  USE modglobal, ONLY: t_global
54  USE moddatatypes
55  USE modparameters
56  USE moderror
57 
58  IMPLICIT NONE
59 
60  PRIVATE
61  PUBLIC :: rflu_createoctree,rflu_buildoctree,rflu_queryoctree, &
62  rflu_destroyoctree
63 
64  SAVE
65 
66 ! ******************************************************************************
67 ! Declarations and definitions
68 ! ******************************************************************************
69 
70  INTEGER, PARAMETER :: MTEST = 5000, &
71  MNUMP = 100
72 
73  INTEGER :: MOCTR,NNODE
74 
75  INTEGER, DIMENSION(:), ALLOCATABLE :: ICHK,NLINK
76  INTEGER, DIMENSION(:,:), ALLOCATABLE :: NOCTR
77 
78  REAL(RFREAL) :: DISMIN
79  REAL(RFREAL), DIMENSION(2) :: XFAR,YFAR,ZFAR
80  REAL(RFREAL), DIMENSION(:), ALLOCATABLE :: X,Y,Z
81 
82 ! ******************************************************************************
83 ! Routines
84 ! ******************************************************************************
85 
86  CONTAINS
87 
88 ! ==============================================================================
89 ! Create octree
90 ! ==============================================================================
91 
92  SUBROUTINE rflu_createoctree(global,nPoints)
93 
94  IMPLICIT NONE
95 
96  INTEGER, INTENT(IN) :: npoints
97  TYPE(t_global), POINTER :: global
98 
99  INTEGER :: errorflag
100 
101  CALL registerfunction(global,'RFLU_CreateOctree',&
102  'RFLU_ModOctree.F90')
103 
104  nnode = npoints
105  moctr = 2*nnode
106 
107  ALLOCATE(ichk(nnode),stat=errorflag)
108  global%error = errorflag
109  IF ( global%error /= err_none ) THEN
110  CALL errorstop(global,err_allocate,__line__,'ICHK')
111  END IF ! global%error
112 
113  ALLOCATE(nlink(nnode),stat=errorflag)
114  global%error = errorflag
115  IF ( global%error /= err_none ) THEN
116  CALL errorstop(global,err_allocate,__line__,'NLINK')
117  END IF ! global%error
118 
119  ALLOCATE(noctr(2,moctr),stat=errorflag)
120  global%error = errorflag
121  IF ( global%error /= err_none ) THEN
122  CALL errorstop(global,err_allocate,__line__,'NOCTR')
123  END IF ! global%error
124 
125  ALLOCATE(x(nnode),stat=errorflag)
126  global%error = errorflag
127  IF ( global%error /= err_none ) THEN
128  CALL errorstop(global,err_allocate,__line__,'X')
129  END IF ! global%error
130 
131  ALLOCATE(y(nnode),stat=errorflag)
132  global%error = errorflag
133  IF ( global%error /= err_none ) THEN
134  CALL errorstop(global,err_allocate,__line__,'Y')
135  END IF ! global%error
136 
137  ALLOCATE(z(nnode),stat=errorflag)
138  global%error = errorflag
139  IF ( global%error /= err_none ) THEN
140  CALL errorstop(global,err_allocate,__line__,'Z')
141  END IF ! global%error
142 
143  CALL deregisterfunction(global)
144 
145  END SUBROUTINE rflu_createoctree
146 
147 
148 ! ==============================================================================
149 ! Build octree
150 ! ==============================================================================
151 
152  SUBROUTINE rflu_buildoctree(XI,YI,ZI,XLOW,XUPP,YLOW,YUPP,ZLOW,ZUPP)
153 !
154 ! ******************************************************************
155 ! * *
156 ! * GIVEN A MESH DEFINED AT THE SET OF POINTS *
157 ! * (X(N),Y(N),Z(N), N=1,NNODE), STORE THEM IN AN OCTREE DATA *
158 ! * STRUCTURE FOR USE IN FUTURE SEARCH QUERIES. *
159 ! * *
160 ! ******************************************************************
161 ! ******************************************************************
162 ! * *
163 ! * COPYRIGHT (C) TIM BAKER 2002 *
164 ! * *
165 ! ******************************************************************
166 !
167 ! ******************************************************************
168 !
169 
170 ! --- arguments
171 
172  REAL(RFREAL), INTENT(IN) :: xlow,xupp,ylow,yupp,zlow,zupp
173  REAL(RFREAL), INTENT(IN) :: xi(nnode),yi(nnode),zi(nnode)
174 
175 ! --- locals
176 
177  INTEGER :: n
178 !
179 ! ******************************************************************
180 !
181 ! --- copy coordinates into internal coordinate array
182 !
183  DO n=1,nnode
184  x(n) = xi(n)
185  y(n) = yi(n)
186  z(n) = zi(n)
187  END DO ! N
188 !
189 ! DETERMINE MINIMUM AND MAXIMUM EXTENT OF ORIGINAL DOMAIN
190 !
191  ichk(1:nnode) = 0
192 
193  xfar(1) = xlow
194  xfar(2) = xupp
195  yfar(1) = ylow
196  yfar(2) = yupp
197  zfar(1) = zlow
198  zfar(2) = zupp
199 !
200 ! INSERT ORIGINAL MESH POINTS INTO OCTREE DATA STRUCTURE
201 !
202  nlink(1) = 0
203  noctr(1,1) = 1
204  noctr(2,1) = 0
205 
206  CALL octree
207 
208  END SUBROUTINE rflu_buildoctree
209 
210 
211 
212 ! =============================================================================
213 ! Construct octree
214 ! =============================================================================
215 
216  SUBROUTINE octree
217 !
218 ! ******************************************************************
219 ! * *
220 ! * SET UP OCTREE DATA STRUCTURE FOR THE POINTS (X(N),Y(N),Z(N)) *
221 ! * WHERE N RANGES FROM 1 TO NNODE *
222 ! * *
223 ! ******************************************************************
224 ! ******************************************************************
225 ! * *
226 ! * COPYRIGHT (C) TIM BAKER 2002 *
227 ! * *
228 ! ******************************************************************
229 !
230 ! ******************************************************************
231 !
232 ! --- locals
233 
234  INTEGER :: i,ioctr,iroot,j,k,l,ltest,n,next,ntest,nxsgn,nysgn,nzsgn
235  INTEGER :: jj(8),nstore(8)
236  REAL(RFREAL) :: tol, &
237  xhalf,xhigh,xlow,xshift,xsize,xsgn, &
238  yhalf,yhigh,ylow,yshift,ysize,ysgn, &
239  zhalf,zhigh,zlow,zshift,zsize,zsgn
240  REAL(RFREAL) :: xroot(2),yroot(2),zroot(2)
241 !
242 ! ******************************************************************
243 !
244  tol = 1.000001_rfreal
245  ioctr = 1
246 !
247 ! LOOP OVER LIST OF POINTS
248 !
249  DO 150 n=2,nnode
250  i = 1
251  xlow = xfar(1)
252  xhigh = xfar(2)
253  ylow = yfar(1)
254  yhigh = yfar(2)
255  zlow = zfar(1)
256  zhigh = zfar(2)
257 
258  IF (noctr(1,i).LT.0) go to 80
259  l = 0
260  next = noctr(1,i)
261  6 next = nlink(next)
262  l = l +1
263  IF (next.NE.0) go to 6
264  10 iroot = i
265  xroot(1) = xlow
266  xroot(2) = xhigh
267  yroot(1) = ylow
268  yroot(2) = yhigh
269  zroot(1) = zlow
270  zroot(2) = zhigh
271  IF (l.EQ.8) go to 30
272  l = l +1
273  nlink(n) = 0
274  IF (l.GT.1) go to 15
275  noctr(1,i) = n
276  go to 150
277  15 next = noctr(1,i)
278  20 ntest = nlink(next)
279  IF (ntest.EQ.0) go to 25
280  next = ntest
281  go to 20
282  25 nlink(next) = n
283  go to 150
284 !
285 ! OCTANT NOW CONTAINS NINE POINTS. FORM EIGHT NEW OCTANTS
286 ! AND AMEND TREE STRUCTURE.
287 !
288  30 DO 35 k=1,8
289  ioctr = ioctr +1
290  IF (ioctr.GT.moctr) go to 210
291  noctr(2,ioctr) = i
292  35 noctr(1,ioctr) = 0
293 !
294 ! ASSIGN POINTS TO NEW OCTANTS
295 !
296  next = noctr(1,i)
297  DO 40 l=1,8
298  nstore(l) = next
299  next = nlink(next)
300  40 jj(l) = 0
301  DO 45 k=1,8
302  j = nstore(k)
303  xshift = x(j) -.5_rfreal*(xlow +xhigh)
304  xsize = max(1.e-9_rfreal,abs(xshift))
305  nxsgn = (int(tol*xshift/xsize) +1)/2*2
306  yshift = y(j) -.5_rfreal*(ylow +yhigh)
307  ysize = max(1.e-9_rfreal,abs(yshift))
308  nysgn = (int(tol*yshift/ysize) +1)/2*2
309  zshift = z(j) -.5_rfreal*(zlow +zhigh)
310  zsize = max(1.e-9_rfreal,abs(zshift))
311  nzsgn = (int(tol*zshift/zsize) +1)/2*2
312  l = 1 +nxsgn/2 +nysgn +2*nzsgn
313  jj(l) = jj(l) +1
314  nlink(j) = 0
315  IF (jj(l).GT.1) go to 41
316  noctr(1,ioctr-8+l) = j
317  go to 45
318  41 next = noctr(1,ioctr-8+l)
319  42 ntest = nlink(next)
320  IF (ntest.EQ.0) go to 43
321  next = ntest
322  go to 42
323  43 nlink(next) = j
324  45 CONTINUE
325 !
326 ! CHECK WHETHER ALL EIGHT POINTS LIE IN ONE OCTANT
327 !
328  noctr(1,i) = 7 -ioctr
329  ltest = 0
330  DO 60 l=1,8
331  IF (jj(l).EQ.8) ltest = l
332  60 CONTINUE
333  IF (ltest.EQ.0) go to 70
334  i = ioctr -8 +ltest
335  xhalf = .5_rfreal*(xlow +xhigh)
336  xsgn = mod(ltest,2)
337  xlow = xsgn*xlow +(1._rfreal -xsgn)*xhalf
338  xhigh = xsgn*xhalf +(1._rfreal -xsgn)*xhigh
339  yhalf = .5_rfreal*(ylow +yhigh)
340  ysgn = isign(1,2*mod(ltest-1,4)-3)
341  ylow = .5_rfreal*((1._rfreal-ysgn)*ylow +(1._rfreal+ysgn)*yhalf)
342  yhigh = .5_rfreal*((1._rfreal-ysgn)*yhalf+(1._rfreal+ysgn)*yhigh)
343  zhalf = .5_rfreal*(zlow +zhigh)
344  zsgn = isign(1,2*ltest-9)
345  zlow = .5_rfreal*((1._rfreal-zsgn)*zlow +(1._rfreal+zsgn)*zhalf)
346  zhigh = .5_rfreal*((1._rfreal-zsgn)*zhalf+(1._rfreal+zsgn)*zhigh)
347  go to 30
348 !
349 ! LOCATE SUB-OCTANT IN WHICH POINT LIES
350 !
351  70 i = iroot
352  xlow = xroot(1)
353  xhigh = xroot(2)
354  ylow = yroot(1)
355  yhigh = yroot(2)
356  zlow = zroot(1)
357  zhigh = zroot(2)
358  80 xhalf = .5_rfreal*(xlow +xhigh)
359  yhalf = .5_rfreal*(ylow +yhigh)
360  zhalf = .5_rfreal*(zlow +zhigh)
361  xshift = x(n) -xhalf
362  xsize = max(1.e-9_rfreal,abs(xshift))
363  nxsgn = (int(tol*xshift/xsize) +1)/2*2
364  yshift = y(n) -yhalf
365  ysize = max(1.e-9_rfreal,abs(yshift))
366  nysgn = (int(tol*yshift/ysize) +1)/2*2
367  zshift = z(n) -zhalf
368  zsize = max(1.e-9_rfreal,abs(zshift))
369  nzsgn = (int(tol*zshift/zsize) +1)/2*2
370  l = 1 +nxsgn/2 +nysgn +2*nzsgn
371  i = -noctr(1,i) +l -1
372  xsgn = mod(l,2)
373  xlow = xsgn*xlow +(1._rfreal -xsgn)*xhalf
374  xhigh = xsgn*xhalf +(1._rfreal -xsgn)*xhigh
375  ysgn = isign(1,2*mod(l-1,4)-3)
376  ylow = .5_rfreal*((1._rfreal-ysgn)*ylow +(1._rfreal+ysgn)*yhalf)
377  yhigh = .5_rfreal*((1._rfreal-ysgn)*yhalf+(1._rfreal+ysgn)*yhigh)
378  zsgn = isign(1,2*l-9)
379  zlow = .5_rfreal*((1._rfreal-zsgn)*zlow +(1._rfreal+zsgn)*zhalf)
380  zhigh = .5_rfreal*((1._rfreal-zsgn)*zhalf+(1._rfreal+zsgn)*zhigh)
381  IF (noctr(1,i).LT.0) go to 80
382  l = 0
383  IF (noctr(1,i).EQ.0) go to 10
384  next = noctr(1,i)
385  85 next = nlink(next)
386  l = l +1
387  IF (next.NE.0) go to 85
388  go to 10
389  150 CONTINUE
390  RETURN
391  210 WRITE (6,610)
392 
393  stop
394 
395  610 FORMAT(5x,'DIMENSION OF NOCTR ARRAY EXCEEDED IN ROUTINE OCTREE.'/ &
396  5x,'INCREASE SIZE OF MOCTR')
397 
398  END SUBROUTINE octree
399 
400 
401 ! ==============================================================================
402 ! Query octree
403 ! ==============================================================================
404 
405  SUBROUTINE rflu_queryoctree(XPT,YPT,ZPT,NUMP,NEIGHP)
406 !
407 ! ******************************************************************
408 ! * *
409 ! * FIND THE NUMP POINTS IN THE OCTREE THAT ARE CLOSEST TO THE *
410 ! * POINT (XPT,YPT,ZPT). THE ADDRESSES OF THE CLOSEST POINTS ARE *
411 ! * GIVEN BY (NEIGHP(M), M = 1, NUMP). *
412 ! * *
413 ! ******************************************************************
414 ! ******************************************************************
415 ! * *
416 ! * COPYRIGHT (C) TIM BAKER 2002 *
417 ! * *
418 ! ******************************************************************
419 !
420 ! ******************************************************************
421 !
422 
423 ! --- arguments
424 
425  INTEGER, INTENT(IN) :: nump
426  INTEGER, INTENT(INOUT) :: neighp(nump)
427  REAL(RFREAL), INTENT(IN) :: xpt,ypt,zpt
428 
429 ! --- locals
430 
431  INTEGER :: i,ic,iflag,ii,itry,istart,j,jj,k,kc,l,lflag,m, &
432  nxsgn,nysgn,nzsgn
433  INTEGER :: nsrch(mtest),ksrch(mtest)
434  REAL(RFREAL) :: dist,dmin,tol,tolpt, &
435  xhalf,xhigh,xhigh0,xl,xlow,xlow0,xshift,xsgn,xsize,xu, &
436  yhalf,yhigh,yhigh0,yl,ylow,ylow0,yshift,ysgn,ysize,yu, &
437  zhalf,zhigh,zhigh0,zl,zlow,zlow0,zshift,zsgn,zsize,zu
438  REAL(RFREAL) :: xoctr(2,mtest),yoctr(2,mtest),zoctr(2,mtest), &
439  xhold(2,mtest),yhold(2,mtest),zhold(2,mtest), &
440  xkeep(2),ykeep(2),zkeep(2)
441 !
442 ! ******************************************************************
443 !
444  tolpt = 1.e-9_rfreal
445  tol = 1.000001_rfreal
446  IF (xpt.LT.xfar(1)-tolpt.OR.xpt.GT.xfar(2)+tolpt) go to 200
447  IF (ypt.LT.yfar(1)-tolpt.OR.ypt.GT.yfar(2)+tolpt) go to 200
448  IF (zpt.LT.zfar(1)-tolpt.OR.zpt.GT.zfar(2)+tolpt) go to 200
449 !
450 ! FIRST FIND OCTANT WHICH CONTAINS POINT (XPT,YPT,ZPT)
451 !
452  m = 0
453  i = 1
454  xlow = xfar(1)
455  xhigh = xfar(2)
456  ylow = yfar(1)
457  yhigh = yfar(2)
458  zlow = zfar(1)
459  zhigh = zfar(2)
460  xlow0 = xlow
461  xhigh0 = xhigh
462  ylow0 = ylow
463  yhigh0 = yhigh
464  zlow0 = zlow
465  zhigh0 = zhigh
466  istart = 1
467  IF (noctr(1,i).GT.0) go to 50
468 !
469 ! LOCATE SUB-OCTANT IN WHICH POINT LIES
470 !
471  20 xhalf = .5_rfreal*(xlow +xhigh)
472  yhalf = .5_rfreal*(ylow +yhigh)
473  zhalf = .5_rfreal*(zlow +zhigh)
474  xshift = xpt -xhalf
475  xsize = max(1.e-9_rfreal,abs(xshift))
476  nxsgn = (int(tol*xshift/xsize) +1)/2*2
477  yshift = ypt -yhalf
478  ysize = max(1.e-9_rfreal,abs(yshift))
479  nysgn = (int(tol*yshift/ysize) +1)/2*2
480  zshift = zpt -zhalf
481  zsize = max(1.e-9_rfreal,abs(zshift))
482  nzsgn = (int(tol*zshift/zsize) +1)/2*2
483  l = 1 +nxsgn/2 +nysgn +2*nzsgn
484  i = -noctr(1,i) +l -1
485  xsgn = mod(l,2)
486  xlow = xsgn*xlow +(1. -xsgn)*xhalf
487  xhigh = xsgn*xhalf +(1. -xsgn)*xhigh
488  ysgn = isign(1,2*mod(l-1,4)-3)
489  ylow = .5_rfreal*((1._rfreal-ysgn)*ylow +(1._rfreal+ysgn)*yhalf)
490  yhigh = .5_rfreal*((1._rfreal-ysgn)*yhalf+(1._rfreal+ysgn)*yhigh)
491  zsgn = isign(1,2*l-9)
492  zlow = .5_rfreal*((1._rfreal-zsgn)*zlow +(1._rfreal+zsgn)*zhalf)
493  zhigh = .5_rfreal*((1._rfreal-zsgn)*zhalf+(1._rfreal+zsgn)*zhigh)
494  IF (noctr(1,i).LT.0) go to 20
495  xlow0 = xlow
496  xhigh0 = xhigh
497  ylow0 = ylow
498  yhigh0 = yhigh
499  zlow0 = zlow
500  zhigh0 = zhigh
501  istart = i
502 !
503 ! SEARCH FOR POINT NEIGHP(M) THAT IS NEAREST TO (XPT,YPT,ZPT)
504 !
505  50 m = m +1
506  neighp(m) = 1
507  dismin = (xfar(2) -xfar(1))**2 +(yfar(2) -yfar(1))**2 &
508  +(zfar(2) -zfar(1))**2
509  xlow = xlow0
510  xhigh = xhigh0
511  ylow = ylow0
512  yhigh = yhigh0
513  zlow = zlow0
514  zhigh = zhigh0
515  i = istart
516  IF (noctr(1,i).EQ.0) go to 65
517  j = noctr(1,i)
518  60 dist = (xpt -x(j))**2 +(ypt -y(j))**2 &
519  +(zpt -z(j))**2
520  IF (dist.LT.dismin.AND.ichk(j).EQ.0) THEN
521  neighp(m) = j
522  dismin = dist
523  ENDIF
524  j = nlink(j)
525  IF (j.NE.0) go to 60
526  65 dmin = sqrt(dismin)
527  xl = xpt -dmin
528  xu = xpt +dmin
529  yl = ypt -dmin
530  yu = ypt +dmin
531  zl = zpt -dmin
532  zu = zpt +dmin
533  70 IF (xl.GT.xlow.AND.xu.LT.xhigh.AND. &
534  yl.GT.ylow.AND.yu.LT.yhigh.AND. &
535  zl.GT.zlow.AND.zu.LT.zhigh) THEN
536  IF (m.EQ.nump) THEN
537  DO 72 m=1,nump
538  ichk(neighp(m)) = 0
539  72 CONTINUE
540  RETURN
541  ELSE
542  ichk(neighp(m)) = 1
543  go to 50
544  ENDIF
545  ENDIF
546  IF (i.EQ.1) THEN
547  IF (m.EQ.nump) THEN
548  DO 74 m=1,nump
549  ichk(neighp(m)) = 0
550  74 CONTINUE
551  RETURN
552  ELSE
553  ichk(neighp(m)) = 1
554  go to 50
555  ENDIF
556  ENDIF
557  iflag = i
558  i = noctr(2,i)
559  lflag = iflag +noctr(1,i) +1
560  xsgn = mod(lflag,2)
561  xoctr(1,1) = (2._rfreal -xsgn)*xlow -(1._rfreal -xsgn)*xhigh
562  xoctr(2,1) = (1._rfreal +xsgn)*xhigh -xsgn*xlow
563  ysgn = isign(1,2*mod(lflag-1,4)-3)
564  yoctr(1,1) = .5_rfreal*((3._rfreal+ysgn)*ylow -(1._rfreal+ysgn)*yhigh)
565  yoctr(2,1) = .5_rfreal*((3._rfreal-ysgn)*yhigh -(1._rfreal-ysgn)*ylow)
566  zsgn = isign(1,2*lflag-9)
567  zoctr(1,1) = .5_rfreal*((3._rfreal+zsgn)*zlow -(1._rfreal+zsgn)*zhigh)
568  zoctr(2,1) = .5_rfreal*((3._rfreal-zsgn)*zhigh -(1._rfreal-zsgn)*zlow)
569  xkeep(1) = xoctr(1,1)
570  xkeep(2) = xoctr(2,1)
571  ykeep(1) = yoctr(1,1)
572  ykeep(2) = yoctr(2,1)
573  zkeep(1) = zoctr(1,1)
574  zkeep(2) = zoctr(2,1)
575 !
576 ! EXAMINE PARENT OF OCTANT AND ALL ITS OFFSPRING
577 !
578  kc = 1
579  ksrch(1) = i
580  75 ic = 0
581  DO 90 j=1,kc
582  ii = ksrch(j)
583  DO 90 k=1,8
584  itry = -noctr(1,ii) +k -1
585  IF (itry.EQ.iflag) go to 90
586  xhalf = .5_rfreal*(xoctr(1,j) +xoctr(2,j))
587  xsgn = mod(k,2)
588  xlow = xsgn*xoctr(1,j) +(1._rfreal -xsgn)*xhalf
589  xhigh = xsgn*xhalf +(1._rfreal -xsgn)*xoctr(2,j)
590  yhalf = .5_rfreal*(yoctr(1,j) +yoctr(2,j))
591  ysgn = isign(1,2*mod(k-1,4)-3)
592  ylow = .5_rfreal*((1._rfreal-ysgn)*yoctr(1,j) +(1._rfreal+ysgn)*yhalf)
593  yhigh = .5_rfreal*((1._rfreal-ysgn)*yhalf +(1._rfreal+ysgn)*yoctr(2,j))
594  zhalf = .5_rfreal*(zoctr(1,j) +zoctr(2,j))
595  zsgn = isign(1,2*k-9)
596  zlow = .5_rfreal*((1._rfreal-zsgn)*zoctr(1,j)+(1._rfreal+zsgn)*zhalf)
597  zhigh = .5_rfreal*((1._rfreal-zsgn)*zhalf +(1._rfreal+zsgn)*zoctr(2,j))
598  IF (xl.GT.xhigh.OR.xu.LT.xlow.OR. &
599  yl.GT.yhigh.OR.yu.LT.ylow.OR. &
600  zl.GT.zhigh.OR.zu.LT.zlow) go to 90
601  IF (noctr(1,itry).GE.0) go to 80
602  ic = ic +1
603  IF (ic.GT.mtest) go to 210
604  xhold(1,ic) = xlow
605  xhold(2,ic) = xhigh
606  yhold(1,ic) = ylow
607  yhold(2,ic) = yhigh
608  zhold(1,ic) = zlow
609  zhold(2,ic) = zhigh
610  nsrch(ic) = itry
611  go to 90
612  80 IF (noctr(1,itry).EQ.0) go to 90
613  jj = noctr(1,itry)
614  85 dist = (xpt -x(jj))**2 +(ypt -y(jj))**2 &
615  +(zpt -z(jj))**2
616  IF (dist.LT.dismin.AND.ichk(jj).EQ.0) THEN
617  neighp(m) = jj
618  dismin = dist
619  ENDIF
620  jj = nlink(jj)
621  IF (jj.NE.0) go to 85
622  dmin = sqrt(dismin)
623  xl = xpt -dmin
624  xu = xpt +dmin
625  yl = ypt -dmin
626  yu = ypt +dmin
627  zl = zpt -dmin
628  zu = zpt +dmin
629  90 CONTINUE
630  IF (ic.EQ.0) THEN
631  xlow = xkeep(1)
632  xhigh = xkeep(2)
633  ylow = ykeep(1)
634  yhigh = ykeep(2)
635  zlow = zkeep(1)
636  zhigh = zkeep(2)
637  go to 70
638  ENDIF
639  kc = ic
640  DO 95 j=1,kc
641  ksrch(j) = nsrch(j)
642  xoctr(1,j) = xhold(1,j)
643  xoctr(2,j) = xhold(2,j)
644  yoctr(1,j) = yhold(1,j)
645  yoctr(2,j) = yhold(2,j)
646  zoctr(1,j) = zhold(1,j)
647  zoctr(2,j) = zhold(2,j)
648  95 CONTINUE
649  go to 75
650  200 WRITE (6,600) xpt,ypt,zpt
651  stop
652  210 WRITE (6,610)
653  stop
654 
655  600 FORMAT(5x,'POINT XPT = ',f12.4,' YPT = ',f12.4,' ZPT = ',f12.4/ &
656  5x,'LIES OUTSIDE CONVEX HULL'/ &
657  5x,'PROGRAM STOPPED IN ROUTINE OCT_SEARCH')
658  610 FORMAT(5x,'DIMENSION OF NSRCH AND KSRCH ARRAYS EXCEEDED'/ &
659  5x,'IN ROUTINE OCT_SEARCH. INCREASE SIZE OF MTEST')
660 
661  END SUBROUTINE rflu_queryoctree
662 
663 
664 ! ==============================================================================
665 ! Destroy octree
666 ! ==============================================================================
667 
668  SUBROUTINE rflu_destroyoctree(global)
669 
670  IMPLICIT NONE
671 
672  TYPE(t_global), POINTER :: global
673 
674  INTEGER :: errorflag
675 
676  CALL registerfunction(global,'RFLU_DestroyOctree',&
677  'RFLU_ModOctree.F90')
678 
679  DEALLOCATE(ichk,stat=errorflag)
680  global%error = errorflag
681  IF ( global%error /= err_none ) THEN
682  CALL errorstop(global,err_deallocate,__line__,'ICHK')
683  END IF ! global%error
684 
685  DEALLOCATE(nlink,stat=errorflag)
686  global%error = errorflag
687  IF ( global%error /= err_none ) THEN
688  CALL errorstop(global,err_deallocate,__line__,'NLINK')
689  END IF ! global%error
690 
691  DEALLOCATE(noctr,stat=errorflag)
692  global%error = errorflag
693  IF ( global%error /= err_none ) THEN
694  CALL errorstop(global,err_deallocate,__line__,'NOCTR')
695  END IF ! global%error
696 
697  DEALLOCATE(x,stat=errorflag)
698  global%error = errorflag
699  IF ( global%error /= err_none ) THEN
700  CALL errorstop(global,err_deallocate,__line__,'X')
701  END IF ! global%error
702 
703  DEALLOCATE(y,stat=errorflag)
704  global%error = errorflag
705  IF ( global%error /= err_none ) THEN
706  CALL errorstop(global,err_deallocate,__line__,'Y')
707  END IF ! global%error
708 
709  DEALLOCATE(z,stat=errorflag)
710  global%error = errorflag
711  IF ( global%error /= err_none ) THEN
712  CALL errorstop(global,err_deallocate,__line__,'Z')
713  END IF ! global%error
714 
715  CALL deregisterfunction(global)
716 
717  END SUBROUTINE rflu_destroyoctree
718 
719 
720 
721 ! ******************************************************************************
722 ! End
723 ! ******************************************************************************
724 
725 END MODULE rflu_modoctree
726 
727 
728 ! ******************************************************************************
729 !
730 ! RCS Revision history:
731 !
732 ! $Log: RFLU_ModOctree.F90,v $
733 ! Revision 1.9 2008/12/06 08:44:23 mtcampbe
734 ! Updated license.
735 !
736 ! Revision 1.8 2008/11/19 22:17:34 mtcampbe
737 ! Added Illinois Open Source License/Copyright
738 !
739 ! Revision 1.7 2004/02/16 23:49:32 haselbac
740 ! Changed setting of MOCTR (for PLAG testing) and allocation of NOCTR
741 !
742 ! Revision 1.6 2002/11/27 22:34:06 haselbac
743 ! Changed estimate of MOCTR to resolve occasional error
744 !
745 ! Revision 1.5 2002/10/08 15:49:21 haselbac
746 ! {IO}STAT=global%error replaced by {IO}STAT=errorFlag - SGI problem
747 !
748 ! Revision 1.4 2002/10/05 19:15:22 haselbac
749 ! Some cleaning up and fixed bug
750 !
751 ! Revision 1.3 2002/09/09 15:11:46 haselbac
752 ! global now under regions, bug fixes by me and by Tim Baker, pass bounding box
753 !
754 ! Revision 1.2 2002/07/25 14:53:03 haselbac
755 ! Added improved estimate for MOCTR based on Tim Bakers email
756 !
757 ! Revision 1.1 2002/06/27 15:48:16 haselbac
758 ! Initial revision
759 !
760 ! *****************************************************************************
761 
762 
763 
764 
765 
766 
767 
768 
kc
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE kc
Definition: RFLO_FindSourceCell.h:27
Aff_transformation_2::m
FT m(int i, int j) const
Definition: Aff_transformation_2.h:228
rflu_modoctree::rflu_queryoctree
subroutine, public rflu_queryoctree(XPT, YPT, ZPT, NUMP, NEIGHP)
Definition: RFLU_ModOctree.F90:405
k
j indices k indices k
Definition: Indexing.h:6
max
Vector_n max(const Array_n_const &v1, const Array_n_const &v2)
Definition: Vector_n.h:354
moderror::registerfunction
subroutine registerfunction(global, funName, fileName)
Definition: ModError.F90:449
ic
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ic
Definition: RFLO_FindSourceCell.h:27
sqrt
double sqrt(double d)
Definition: double.h:73
rflu_modoctree::rflu_buildoctree
subroutine, public rflu_buildoctree(XI, YI, ZI, XLOW, XUPP, YLOW, YUPP, ZLOW, ZUPP)
Definition: RFLU_ModOctree.F90:152
rflu_modoctree::rflu_destroyoctree
subroutine, public rflu_destroyoctree(global)
Definition: RFLU_ModOctree.F90:668
i
blockLoc i
Definition: read.cpp:79
n
const NT & n
Definition: rational_rotation.h:72
rflu_modoctree
Definition: RFLU_ModOctree.F90:51
j
j indices j
Definition: Indexing.h:6
modglobal::t_global
Definition: ModGlobal.F90:50
rflu_modoctree::octree
subroutine octree
Definition: RFLU_ModOctree.F90:216
moderror::errorstop
subroutine errorstop(global, errorCode, errorLine, addMessage)
Definition: ModError.F90:483
Element_node_enumerator::next
void next()
Go to the next element within the connectivity tables of a pane.
Definition: Element_accessors.C:79
dist
long double dist(long double *coord1, long double *coord2, int size)
Definition: utilities/rc_cmp/util.cpp:305
moderror::deregisterfunction
subroutine deregisterfunction(global)
Definition: ModError.F90:469
rflu_modoctree::rflu_createoctree
subroutine, public rflu_createoctree(global, nPoints)
Definition: RFLU_ModOctree.F90:92