1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
|
*DECK XELLSR
SUBROUTINE XELLSR( NDIM, NCP, NSUR, MAXREM, REMESH,
> INDEL, MINDIM, MAXDIM, ICOORD, ICUR, INCR,
> TRKORI, TRKDIR, TRKCUT, NSCUT, NCROS,
> TOTLEN)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Find the beginning and ending surfaces crossed by a track.
*
*Copyright:
* Copyright (C) 1990 Ecole Polytechnique de Montreal
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version
*
*Author(s): R. Roy
*
*Parameters: input
* NDIM number of dimension (2 or 3).
* NCP number of cylindres of a type + 3.
* NSUR number of surfaces.
* MAXREM max number of real mesh values in REMESH.
* REMESH real mesh values (rect/cyl).
* INDEL numbering of surfaces and zones.
* MINDIM min index values for all axes (rect/cyl).
* MAXDIM max index values for all axes (rect/cyl).
* ICOORD principal axes direction (X/Y/Z) for meshes.
* ICUR current zonal location for a track segment.
* INCR increment direction for next track segment.
* TRKORI origin of a track.
* TRKDIR direction of a track in all axes.
*
*Parameters: output
* TRKCUT points where track cut the domain.
* NSCUT surface where the track begins/ends.
* NCROS number of surface crossing.
* TOTLEN total length of the track.
*
*-----------------------------------------------------------------------
*
IMPLICIT NONE
*
INTEGER NDIM, NCP, NSUR, MAXREM, NCROS
REAL TRKCUT(3,2), REMESH(MAXREM), TRKDIR(3), TRKORI(3), TOTLEN,
> TKBEG1, TKBEG2, R2BEG
INTEGER MINDIM(NCP), MAXDIM(NCP), ICUR(NCP), INCR(NCP),
> ICOORD(NCP), IFACUT(2), ISFCUT(2),
> IORD(4), NSCUT(2), INDEL(4,*)
INTEGER IOUT
PARAMETER ( IOUT=6 )
REAL ANORM2, CENTRE, A, B, XYZP2, CONBEG, CONEND, CON,
> XYZP1
INTEGER I, J, NUM, NCRBEG, NCREND, NP1, NUMP1, NP2, NUMP2,
> N, NUMP0, K, IBEGIN, KELSUR, KWW, IDM
*
ANORM2(A,B)= A*A + B*B
CENTRE(I,J)= REMESH( MAXDIM(I-1) + J )
NUM(J)= J + 1 - NSUR
NUMP2= 0
IFACUT(:2)=0
ISFCUT(:2)=0
*
IF( NDIM.EQ.2 )THEN
NP2= 3
NUMP2 = ICOORD(NP2)
XYZP2= 0.0
ENDIF
*
* IF THERE ARE NO CYLINDER AT ALL
NSCUT(1)= 0
NSCUT(2)= 0
NCRBEG= 0
NCREND= 0
CONBEG=+1.0E+36
CONEND=-1.0E+36
*
* FING BEGINNING AND ENDING POINTS OF THE TRACK
DO 75 N = 1, NDIM
NUMP0 = ICOORD(N )
IF( INCR(NUMP0).EQ.0 ) GO TO 75
NP1 = MOD(N ,NDIM)+1
NUMP1 = ICOORD(NP1)
IF( NDIM.EQ.3 )THEN
NP2 = MOD(N+1 ,NDIM)+1
NUMP2 = ICOORD(NP2)
ENDIF
DO 70 IDM = MINDIM(N), MAXDIM(N), MAXDIM(N)-MINDIM(N)
CON = (REMESH(IDM) - TRKORI(NUMP0)) / TRKDIR(NUMP0)
XYZP1 = TRKORI(NUMP1) + CON * TRKDIR(NUMP1)
IF( XYZP1.LT.REMESH(MINDIM(NP1)).OR.
> XYZP1.GT.REMESH(MAXDIM(NP1))) GO TO 70
IF( NDIM.EQ.3 )THEN
XYZP2 = TRKORI(NUMP2) + CON * TRKDIR(NUMP2)
IF( XYZP2.LT.REMESH(MINDIM(NP2)).OR.
> XYZP2.GT.REMESH(MAXDIM(NP2))) GO TO 70
ENDIF
IF( CON.LT.CONBEG )THEN
NCRBEG=1
NCREND=MAX(1,NCREND)
IFACUT(1)= NUMP0
ISFCUT(1)= IDM
IF( IDM.EQ.MINDIM(N) ) ISFCUT(1)= ISFCUT(1)-1
CONBEG=CON
TRKCUT(NUMP0,1)= REMESH(IDM)
TRKCUT(NUMP1,1)= XYZP1
TRKCUT(NUMP2,1)= XYZP2
ENDIF
IF( CON.GT.CONEND )THEN
NCREND=2
NCRBEG=MIN(2,NCRBEG)
IFACUT(2)= NUMP0
ISFCUT(2)= IDM
IF( IDM.EQ.MINDIM(N) ) ISFCUT(2)= ISFCUT(2)-1
CONEND=CON
TRKCUT(NUMP0,2)= REMESH(IDM)
TRKCUT(NUMP1,2)= XYZP1
TRKCUT(NUMP2,2)= XYZP2
ENDIF
70 CONTINUE
75 CONTINUE
NCROS = NCREND + NCRBEG
TOTLEN= CONEND - CONBEG
IF( NCROS.EQ.0 ) GO TO 1000
NCROS = NCREND + 1 - NCRBEG
*
* FIND BEGINNING AND ENDING SURFACE NUMBERS
DO 900 K= NCRBEG, NCREND
DO 90 I = 1, NDIM
N = ICOORD(I)
ICUR(I)= MINDIM(I)
DO 80 J = MINDIM(I), MAXDIM(I)-1
IF(TRKCUT(N,K).GE.REMESH(J)) ICUR(I)= J
80 CONTINUE
90 CONTINUE
ICUR(IFACUT(K))= ISFCUT(K)
IBEGIN= MAXDIM(3) + 3
DO 110 I = 4, NCP
N = ICOORD(I)
NP1 = MOD(N ,3) + 1
NP2 = MOD(N+1,3) + 1
TKBEG1= CENTRE(I,1) - TRKCUT(NP1,K)
TKBEG2= CENTRE(I,2) - TRKCUT(NP2,K)
R2BEG = ANORM2(TKBEG1,TKBEG2)
ICUR(I) = IBEGIN - 1
DO 100 J = IBEGIN, MAXDIM(I)
IF( R2BEG .GE. REMESH(J) )ICUR(I)= J
100 CONTINUE
IBEGIN= MAXDIM(I) + 3
110 CONTINUE
*
* FIND IORD(4) FOR LOCATION IN THE INDEX VECTOR
DO 115 I= 1,NCP
IORD(MIN(4,I))= ICUR(I)
IF( I.GT.3.AND.ICUR(I).LT.MAXDIM(I)) GOTO 116
115 CONTINUE
IORD(4)= 0
116 CONTINUE
*
* FIND NSCUT=BEGINNING/ENDING SURFACE #S
KELSUR= NSUR
INDEL(1,NUM(0))= IORD(1)
INDEL(2,NUM(0))= IORD(2)
INDEL(3,NUM(0))= IORD(3)
INDEL(4,NUM(0))= IORD(4)
880 CONTINUE
IF( IORD(1).EQ.INDEL(1,NUM(KELSUR)).AND.
> IORD(2).EQ.INDEL(2,NUM(KELSUR)).AND.
> IORD(3).EQ.INDEL(3,NUM(KELSUR)).AND.
> IORD(4).EQ.INDEL(4,NUM(KELSUR)) ) GO TO 890
KELSUR= KELSUR + 1
GO TO 880
890 NSCUT(K)= KELSUR
IF( KELSUR.EQ.0 )THEN
WRITE(IOUT,*) ' BAD SURFACE IDENTIFICATION'
WRITE(IOUT,*) ' NSCUT=', NSCUT(K)
WRITE(IOUT,*) 'TRKCUT=', (TRKCUT(KWW,K),KWW=1,3)
WRITE(IOUT,*) ' IORD=', IORD
WRITE(IOUT,*) ' ICUR=', (ICUR(KWW),KWW=1,NCP)
CALL XABORT('XELLSR: BAD SURFACE IDENTIFICATION')
ENDIF
900 CONTINUE
*
1000 RETURN
END
|
