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
188
189
190
191
192
193
|
*DECK NXTIHA
FUNCTION NXTIHA(POSHEX ,PINPOS,VOLINT)
*
*----------
*
*Purpose:
* Compute the volume of intersection between
* a 2--D hexagon and an annular pin.
*
*Copyright:
* Copyright (C) 2010 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): G. Marleau.
*
*Parameters: input
* POSHEX spatial description of the hexagon with:
* POSHEX(0) the dimension of one of its sides;
* POSHEX(1) the $X$ position of hexagon center;
* POSHEX(2) the $Y$ position of hexagon center.
* PINPOS spatial description of the annular pin region with:
* PINPOS(0) the radius of the annular pin;
* PINPOS(1) the $X$ position of the annular pin center;
* PINPOS(2) the $Y$ position of the annular pin center.
*
*Parameters: output
* NXTIHA type of intersection between haxagon and annular pin, where:
* = 0 means that there is no intersection
* between the two regions;
* = 1 means that the hexagon
* is all located inside the annular pin;
* = 2 means that the annular pin
* is all located inside the hexagon;
* =-1 means that the intersection between
* the hexagon and the annular pin is partial.
* VOLINT 2-D volume of intersection (area) between hexagon and
* annular pin.
*
*Reference:
* G. Marleau,
* New Geometries Processing in DRAGON: The NXT: Module,
* Report IGE-260, Polytechnique Montreal,
* Montreal, 2005.
*
*----
*
IMPLICIT NONE
*----
* Subroutine arguments
*----
INTEGER NXTIHA
DOUBLE PRECISION POSHEX(0:2),PINPOS(0:2)
DOUBLE PRECISION VOLINT
*----
* Local parameters
*----
INTEGER IOUT
CHARACTER NAMSBR*6
PARAMETER (IOUT=6,NAMSBR='NXTIHA')
INTEGER IPRLOC
PARAMETER (IPRLOC=10)
DOUBLE PRECISION DCUTOF
PARAMETER (DCUTOF=1.0D-8)
DOUBLE PRECISION DZERO,DONE,DHALF,DSQ3O2
PARAMETER (DZERO=0.0D0,DONE=1.0D0,
> DHALF=0.5D0,DSQ3O2=0.86602540378444D0)
*----
* Functions
*----
INTEGER NXTITA
DOUBLE PRECISION XDRCST,PI
*----
* Local variables
*----
INTEGER IDIR,ICORN,IFACE,NFPINS,NCIN,ITRI,INTTRI
DOUBLE PRECISION POSTRI(2,3),RADC,DISTF,DISTC,RADP2,
> VTPINT
*----
* Data
*----
DOUBLE PRECISION CORNRH(2,6),DIRFAC(2,6)
SAVE CORNRH,DIRFAC
DATA CORNRH / 0.86602540378444D0,-0.5D0,
> 0.86602540378444D0, 0.5D0,
> 0.0D0 , 1.0D0,
> -0.86602540378444D0, 0.5D0,
> -0.86602540378444D0,-0.5D0,
> 0.0D0 ,-1.0D0/
DATA DIRFAC /-1.0D0 , 0.0D0,
> -0.5D0 ,-0.86602540378444D0,
> 0.5D0 ,-0.86602540378444D0,
> 1.0D0 , 0.0D0,
> 0.5D0 , 0.86602540378444D0,
> -0.5D0 , 0.86602540378444D0/
*----
* Print header if required
*----
IF(IPRLOC .GE. 200) THEN
WRITE(IOUT,6000) NAMSBR
WRITE(IOUT,6010) (POSHEX(IFACE),IFACE=0,2)
WRITE(IOUT,6011) (PINPOS(IFACE),IFACE=0,2)
ENDIF
*----
* Initialize PI, NXTIHA and VOLINT
*----
PI=XDRCST('Pi',' ')
NXTIHA=0
VOLINT=DZERO
*----
* Evaluate distance from FACES to pin center
*----
RADP2=PINPOS(0)**2
NFPINS=0
NCIN=0
ICORN=1
POSTRI(1,ICORN)=POSHEX(1)
POSTRI(2,ICORN)=POSHEX(2)
DO IFACE=1,6
DISTF=DZERO
RADC=DZERO
* write(6,*) 'CORNRH',(CORNRH(IDIR,IFACE),IDIR=1,2)
* write(6,*) 'DIRFAC',(DIRFAC(IDIR,IFACE),IDIR=1,2)
DO IDIR=1,2
DISTC=PINPOS(IDIR)-(POSHEX(IDIR)+POSHEX(0)*CORNRH(IDIR,IFACE))
DISTF=DISTF+DISTC*DIRFAC(IDIR,IFACE)
RADC=RADC+DISTC**2
ENDDO
* write(6,*) DISTC,DISTF,RADC,PINPOS(0)
IF(DISTF .LT. -PINPOS(0)) THEN
*----
* Pin outside hexagon
* Return
*----
NXTIHA=0
VOLINT=DZERO
RETURN
ELSE IF(DISTF .GE. PINPOS(0)) THEN
NFPINS=NFPINS+1
ENDIF
IF(RADC .LT. RADP2) THEN
NCIN=NCIN+1
ENDIF
ENDDO
IF(NFPINS .EQ. 6) THEN
NXTIHA=2
VOLINT=PI*RADP2
ELSE IF(NCIN .EQ. 6) THEN
NXTIHA=1
VOLINT=3.0D0*DSQ3O2*POSHEX(0)**2
ELSE
NXTIHA=-1
*----
* First five triangles
*----
DO ITRI=1,5
DO ICORN=2,3
DO IDIR=1,2
POSTRI(IDIR,ICORN)=POSTRI(IDIR,ICORN)
> +POSHEX(0)*CORNRH(IDIR,IFACE+ICORN-2)
ENDDO
ENDDO
INTTRI=NXTITA(POSTRI,PINPOS,VTPINT)
VOLINT=VOLINT+VTPINT
ENDDO
*----
* Last five triangles
*----
DO IDIR=1,2
POSTRI(IDIR,ICORN)=POSTRI(IDIR,ICORN)
> +POSHEX(0)*CORNRH(IDIR,6)
POSTRI(IDIR,ICORN)=POSTRI(IDIR,ICORN)
> +POSHEX(0)*CORNRH(IDIR,1)
ENDDO
ITRI=NXTITA(POSTRI,PINPOS,VTPINT)
VOLINT=VOLINT+VTPINT
ENDIF
IF(IPRLOC .GE. 200) THEN
WRITE(IOUT,6012) NAMSBR,NXTIHA,VOLINT
WRITE(IOUT,6001) NAMSBR
ENDIF
RETURN
*----
* Output formats
*----
6000 FORMAT('(* Output from --',A6,'-- follows ')
6001 FORMAT(' Output from --',A6,'-- completed *)')
6010 FORMAT('POSHEX={',2(F20.10,','),F20.10,'};')
6011 FORMAT('PINPOS={',2(F20.10,','),F20.10,'};')
6012 FORMAT(A6,'={',I5,',',F20.10,'};')
END
|
